使用径向基函数进行曲面插值

你有没有遇到过这样一种情况：在笛卡尔平面内，有一组非均匀分布的点，用于对面高度进行采样，例如地图等高线上的点或表示某些材料属性数据的数据点？你可能还想过，要在这些点之间重建或插值出一个连续的光滑曲面。对于这种想法，你可以使用 COMSOL Multiphysics 的核心功能，即径向基函数来构建这样的曲面。让我们来看看应该如何操作。

径向基函数简介

径向基函数(RBF) 通常是指根据到某点的距离（或半径）定义的任意函数：

(1)

式中，w是径向基函数的权重；\mathbf{c}=(x_c,y_c)是该点的坐标或中心；r是从 XY 平面内的任意其他点到此中心的距离。

径向基函数本身可以是许多不同类型函数的一种。多谐波样条系列通常用于插值，特别是薄板样条函数，这种基函数的形式为：

(2)

仅仅这样一个径向基函数并不是那么吸引人，但我们可以在有限数量的具有不同权重的不同中心上进行求和，并选择性地添加具有权重的线性多项式项a_0, a_1, a_2，从而得到函数：

(3)

如果有足够的中心点，那么一组基函数的总和可以用来表示非常复杂的单值函数。当使用薄板样条基函数时，还有一个额外的优点，就是这个函数在任何地方都是平滑的，并且可以无限微分。

现在让我们来看看如何使用这些径向基函数来插值得到光滑曲面。如果我们获得一组有限的中心点位置\mathbf{c}_1, \mathbf{c}_2, … \mathbf{c}_N，以及它们对应的已知高度z_1, z_2, … z_N，那么可以写一个线性方程组：

(4)

其中，系统矩阵项\phi_{i,j}= \phi ( || \mathbf{c}_i- \mathbf{c}_j||)是各中心点之间计算的径向基函数。

当使用薄板样条基函数时，几乎所有的非对角线项都是非零的，因此这个系统矩阵非常密集。线性系统可以求解所有权重，我们可以计算 xy 平面中任何其他点的加权径向基函数之和，从而得到平滑插值函数。现在，让我们看看如何使用 COMSOL Multiphysics 的核心功能来计算这些权重并将插值函数可视化。

使用 COMSOL Multiphysics 中的径向基函数进行曲面插值

我们从建立一个包含三维组件的模型开始，该模型使用无量纲单位系统，可以在组件 1的设置中选择单位系统。如果我们的数据代表材料属性而不是几何形状，则无量纲单位系统使用起来更简单。

模型中的几何图形由两个特征组成。首先，点特征用于定义点集（可以从文本文件复制坐标列表）。累积选择用于定义所有这些点的命名选择，如下面的屏幕截图所示。此外，还有一个长方体特征，其尺寸略大于数据点的范围，并被定位为包住所有数据点。

显示要插值的数据点的定义和累积选择定义的屏幕截图。

数据点和边界块。

定义几何图形后，我们在刚刚创建的点上定义一个积分组件耦合算子。由于积分是在一组点上完成的，因此这个运算符等效于获取在一组点上计算的表达式总和。接下来，我们定义三个变量，如下面的屏幕截图所示。

首先，变量r = eps+sqrt((dest(x)-x)^2+(dest(y)-y)^2)将用于计算所有中心之间的距离。请注意运算符dest()的用法，这会强制在目标点而不是源点上计算运算符中的表达式。添加一个非常小的非零项（eps是机器精度），以使这个表达式永远不会为零。

接下来，变量phi = r^2*log(r)是方程(2)中的薄板样条基函数。请注意，当半径为零时，此函数应该收敛为零，但由于对数函数的原因，我们确实需要使半径不能完全为零，以便可以在零处也能计算基函数。还值得一提的是，这个函数可以更改为任何其他所需的基函数。

最后，定义RBF = intop1(w*phi)+a0+a1*x+a2*y是方程(3)，即插值曲面本身，其权重尚未计算。请记住，积分组件耦合算子用于获取算子中的表达式在这些点上的总和。

显示变量定义的屏幕截图。

现在几何已经设置好了，所有变量都已定义，我们已经准备好求解径向基函数和多项式项的权重了。这是通过点常微分和微分代数方程接口在我们要插值的点上进行定义而完成的，如下面的屏幕截图所示。我们可以将所有单位设置为无量纲，因为点位置也是无量纲的。这些设置定义了一组未知数w，每个点将具有不同的值。

点常微分方程和微分代数方程接口的设置。

在这个物理场接口中，只需要修改两个特征。首先，需要调整分布式常微分方程的设置，如下图所示。源项定义为z-RBF。由于求解稳态研究时方程中的所有其他项均为零，因此该项表示在所有选定点处RBF=z。使用这个功能，我们可以定义等式(4)的第 1 行到第 N 行。

每个点的源项设置。

其次，我们需要定义等式(4)的最后三行。这是通过全局方程功能完成的，如下面的屏幕截图所示，这三个方程求解权重a0，a1和a2。同样，积分耦合运算符用于获取所有选定点的表达式之和。有了这两个特征，问题的定义就完整了，差不多可以进行求解了。

显示用于多项式权重的全局方程的屏幕截图。

求解这个模型要求我们在所有点上都有一个网格，因此我们在边界框上应用了自由四面体网格，然后使用稳态求解器求解。完成求解后，我们就可以绘制变量RBF的插值函数，如下所示。

穿过所有数据点的平滑且可微分的插值曲面。

将功能打包到一个仿真 App 中

如果你想直接使用这个功能，而不用在自己的模型中设置所有这些功能，欢迎你从 COMSOL案例库中下载我们的演示仿真 App，其中可以从一个用逗号分隔的文件中获取xyz数据点并计算插值曲面。使用此演示 App 最多可以插值 5000 个数据点。

除了对曲面进行计算之外，这个 App 还可以导出描述曲面的完整分析函数，以及它的 COMSOL 格式的 CAD 文件，所有这些都是 COMSOL Multiphysics 的核心功能。CAD 数据是NURBS 表面，因此仅近似表示函数，但对于相当光滑的曲面，精度非常高。该仿真 App 用户界面的屏幕截图如下所示。

计算了插值函数并导出函数和 CAD 曲面文件的仿真 App 屏幕截图。

自己动手尝试

点击下方按钮，下载仿真 App，尝试使用仿真 App 计算一组点之间的插值：

更多资源

如果你有兴趣了解有关模型开发器 和 COMSOL Server™ 的更多信息(它们可用于开发和运行此仿真 App)，请查看以下资源。

• 阅读App 开发器简介手册中的完整说明摘要
• 观看视频，了解如何使用 App 开发器和 COMSOL Server™
• 阅读博客文章，了解如何在各种应用中使用仿真 App

你想使用 COMSOL Multiphysics 做什么？你对径向基函数的功能还有其他疑问吗？欢迎联系我们寻求帮助。