径向基函数插值中形状参数的选取方法

Selection of Shape Parameters in Radial Basis Function Interpolation

摘要: 径向基函数的形状参数对插值精度的影响很大。如何选取形状参数使得插值误差最小的问题，受到国内外学者的广泛关注。结合径向基函数插值的误差理论，本文围绕形状参数的选取方法进行归纳总结，并通过数值实验，对现有的方法对比研究。为了提高变参数径向基函数的插值精度，提出用拉格朗日法和径向基函数法相结合的方式加以改进。

Abstract: The shape parameters of radial basis function have a great influence on interpolation accuracy. How to select the shape parameters to minimize the interpolation error has been widely concerned by scholars at home and abroad. Combined with the error theory of radial basis function interpolation, the selection methods of shape parameters are summarized in this paper. Through numerical experiments, the existing methods are compared. In order to improve the interpolation accuracy of variable parameter radial basis function, the combination of Lagrange method and radial basis function method is proposed.

1. 引言

Frank在文献 [1] 中比较了27种数据处理方法，最后得出结论：径向基函数插值法在所有数值方法当中综合性能最好。到目前为止，径向基函数的理论研究已经基本完善。众多研究和实际应用表明，径向基函数插值中形状参数对插值精度有较大的影响，因此形状参数的选取是至关重要的，国内外众多学者从不同角度对该问题展开了讨论。本文拟对径向基函数插值中形状参数的选取方法进行归纳总结，结合径向基函数插值误差理论，并通过数值实验对现有的方法对比研究。为了提高变参数径向基函数的插值精度，提出用两种插值方法相结合的方式加以改进。

2. 径向基函数插值

E.M. Stein和G. Weiss [2] 对径向基函数是这样定义的：

径向基函数是一个取值仅与离原点距离有关的实值函数 $\varphi$，即 $\Phi \left(x\right)=\varphi \left(\Vert x\Vert \right)$。如果满足： $\Vert x_1\Vert =\Vert x_2\Vert$，那么 $\varphi \left(x_1\right)=\varphi \left(x_2\right)$，其中， $\Vert ·\Vert$ 是标准欧式范数，常用的径向基函数(吴宗敏教授在文献 [3])：

Kriging方法的Gauss分布函数： $\varphi \left(\Vert x_k-x_j\Vert \right)={\text{e}}^{-c^2{\Vert x_k-x_j\Vert }^2}$ ；

Hardy的MQ函数： $\varphi \left(\Vert x_k-x_j\Vert \right)={\left(c^2+{\Vert x_k-x_j\Vert }^2\right)}^{\beta }$ ；

Hardy的逆MQ函数： $\varphi \left(\Vert x_k-x_j\Vert \right)={\left(c^2+{\Vert x_k-x_j\Vert }^2\right)}^{-\beta }$ (其中β是正实数)。

其中，c可以确定基函数的形状，称之为形状参数。

径向基函数插值的定义为：

对于给定的n个样本点 ${\left\{x_j,f_j\right\}}_{j=1}^n\in R^n\otimes R$。选取径向基函数 $\varphi :R_{+}\to R$ 构造径向基函数空间 ${\left\{\varphi \left(\Vert x-x_j\Vert \right)\right\}}_{j=1}^n$，并寻找形如 $S\left(x\right)={\displaystyle \underset{j=1}{\overset{n}{\sum }}{\lambda }_j}\varphi \left(\Vert x-x_j\Vert \right)$ ( ${\lambda }_j$ 为插值系数)的插值函数S(x)，使其满足条件 $S\left(x_j\right)=f_j,j=1,2,\cdots ,n$。

3. 径向基函数插值误差估计

当选定合适的基函数之后，进一步，需要考察径向基函数的拟合效果，即考察验证样本点的误差估计。取m个验证样本点 ${\left\{x_i,f_i\right\}}_{i=1}^m$，S(x_i)为径向基函数在x_i处的预测值。常用的误差估计方法有：

均方根误差 [4] (Root Mean Square Error)： $RMSE=\sqrt{\frac{{\displaystyle \underset{i=1}{\overset{m}{\sum }}{\left[f_i-S\left(x_i\right)\right]}^2}}{m}}$。

最大误差 [5] (Maximum Error)： $ME=\max \left|f_i-S\left(x_i\right)\right|,i=1,2,\cdots ,m$。

相对误差 [6] (Relative Error)： $RE=\frac{\left|f_i-S\left(x_i\right)\right|}{f_i},i=1,2,\cdots ,m$。

交叉验证误差 [7] (Cross Validation Error)：在样本点 $X=\left\{x_{\text{1}},x_{\text{2}},\cdots ,x_m\right\}$ 中，去掉X中的x_i，得到 $X_{-i}=\left\{x_{\text{1}},x_{\text{2}},\cdots ,x_i{}_{-\text{1}},x_i{}_{+\text{1}},\cdots ,x_m\right\}$，用构造径向基模型 ${\tilde{S}}_{-i}\left(x_i\right)$，并用x_i作为验证点，在x_i处的误差可以表示为： $CV{E}_i=\left|S\left(x_i\right)-{\tilde{S}}_{-i}\left(x_i\right)\right|,i=1,2,\cdots ,m$。

相关系数 [8] (Correlation Coefficient)： $r\left(F,S\right)=\frac{Cov\left(F,S\right)}{\sqrt{Var\left(F\right)Var\left(S\right)}}$。其中， $F={\left[f_{\text{1}},f_{\text{2}},\cdots ,f_m\right]}^{\text{T}}$， $S={\left[S\left(x_{\text{1}}\right),S\left(x_{\text{2}}\right),\cdots ,S\left(x_m\right)\right]}^{\text{T}}$， $Cov\left(F,S\right)$ 为F与S的协方差，Var(F)为F的方差，Var(S)为S的方差。比较好的模型预测值和真实值有较高的相关系数，最大值为1。

拟合优度 [7] (Goodness of Fit)： $R^2=1-\frac{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left[f_i-S\left(x_i\right)\right]}^2}}{S{S}_T}$。其中，SS_T称为平方和，即 $S{S}_T={\displaystyle \underset{i=1}{\overset{n}{\sum }}{f}_i^2}-\frac{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{f}_i^2}}{n}$。R²在0到1范围内取值，R²越接近1，模型的精度就越高。

4. 径向基函数插值中形状参数的选取方法

对于径向基函数，形状参数c是一个自由参数。在实际应用过程中c的取值对计算结果有很大的影响，如何选取形状参数使得插值误差最小，一直是研究人员关注的课题。目前，形状参数c的选取有两种观点：一种观点认为参数c是常数，与样本点无关；另一种观点则认为c在每一个样本点处是可变的。以下将分别对这两类径向基函数归纳总结，并通过数值实验比较这些方法的优缺点。

4.1. 常参数径向基函数中形状参数的选取方法

Hardy [9] 早在1971年，就提出了MQ函数中形状参数的计算公式，对其研究的地形问题有很好的拟合能力。常参数径向基函数中，形状参数的选取方法可分为以下两种：

1) 优化误差确定法

显然误差是与形状参数c有关的，记为E(c)，优化误差法的基本思想：把形状参数的选取问题转化为优化误差的问题，即通过最小化误差，为给定的径向基函数选择最优形状参数。

Rippa [4] 定义了代价函数cost(c)来表示插值函数和实际未知函数之间的均方根误差，通过最小化代价函数，得到较好的形状参数。并从插值矩阵的条件数及计算精度、样本点的数量和分布、径向基函数类型三个方面证明了该选取方法的有效性。魏月兴，许林，陈小前 [7] 提出在控制Runge现象 [10] 的同时，最小化交叉验证误差获得形状参数的算法。该算法先确定c的初始值(c_init = mean(d_j)，这里，d_j是x_j和其他样本点之间的最小距离)和步长(l = m/n)。魏月兴，陈小前，许林 [11] 利用MQ函数插值Rosenbroke函数，以插值矩阵条件数小于10¹⁵作为约束条件，最小化RMSE(c)，得到c的最优值。Roque和Ferreira [12] 在MQ形状参数的选择中，应用留一交叉验证(LOOCV) [13] 优化技术模拟给定边值问题的误差函数，并证明该方法适用于一维和二维情况。牛瑞萍 [5] 在利用逆MQ函数求解混合介质中Cauchy热传导反问题时，提出了一种自动选择算法，该算法通过最小化近似温度和给定温度值的残差选择合适的形状参数，并利用Tikhonov正则化 [14] 方法保证问题获得精确稳定的解。

2) 数值实验确定法

到目前为止，在用径向基函数解决具体插值问题或方程问题时，如何确定形状参数使得插值精度尽可能高，还没有普适性结论。因此，很多学者通过改变形状参数c的值，进行重复性插值数值实验，得出误差值，在这些误差值中选择最小误差，其对应的形状参数即为最优。这种确定形状参数最优值的方法比较机械，且缺乏普适性，但数值实验得出的部分结论，却对后面的研究探索有指导性意义。

Fasshauer [15] 用MQ函数对基于网格的非线性偏微分方程进行数值求解时，得出形状参数的选取公式 $c=2/\sqrt{n}$。聂鑫 [6] 用MQ插值求解边值问题时，指出插值精度和稳定性与形状参数相关，一维时，条件数对数与形状参数呈线性关系。形状参数较小时插值误差较大，而过大时又会发生剧烈振荡；随着节点数的增加，形状参数选取范围逐渐减小；只有参数在特定范围时才能取得较好的计算精度。齐静 [16] 用MQ函数、Gauss函数对一元函数 $f\left(x\right)=\sin \left(2x\right)$ 和二元函数 $f\left(x\right)=\sin \left(3x\right)\sin \left(2y\right)$ 作插值的数值实验，均得到：c的取值越小，对应的误差也越小，因此在实际应用中可适当减小c的值。陈风雷 [17] 用MQ函数对函数 $y=\sin x$，分别模拟了p阶导数插值与p重积分插值的数值实验(p取1, 2, 3, 4)，给出了形状参数在积分插值方法中最适宜的取值范围是(0, 1/n)，在导数插值方法中最适宜的取值范围是(1/n, 3)，显然对于导数形状参数的取值要比积分形式稍微大；多重积分插值相对高阶导数插值更稳定、精度高，对于形状参数的选择更灵活。

4.2. 变参数径向基函数中形状参数的选取方法

Kansa [18] 最早发现，在使用MQ插值时，如果使形状参数随空间变化，则插值精度有提高的可能，对于变参数径向基函数 [19]，形状参数的选取方法可分为以下两种：

1) 公式确定法

根据形状参数就是该点对样本空间影响能力的意义，将该样本点与其他样本点之间的最小距离作为该样本点的形状参数，即

$c_j=\min r_{jk},j=1,2,\cdots ,n,j\ne k$ (1)

其中，r_jk表示第j个样本点与第k个样本点的距离。由于样本点间的距离并不相等，每个样本点有不同的影响范围，为了很好地对原始函数进行拟合预测，当采样点个数 $n\to \infty$ 的情况下，Kitayama [20] 提出了如公式(2)所示的形状参数确定方法。

$c_j=\frac{d_{j,\max }}{\sqrt{\alpha }\sqrt[\alpha ]{n-1}}$ (2)

式中， $d_j{}_{,\max }$ 为第j个样本点到其他样本点间的最大距离，α为原函数所含变量的个数。

2) 样本点局部密度确定法

武泽平 [8] [21] 选择Gauss函数作为序列近似优化算法的插值模型，指出形状参数c_j的本质物理意义是样本点对整个样本空间的影响范围大小，据此提出一种基于样本点局部密度确定形状参数的方法。该方法的基本思路为：计算每个样本点的局部密度ρ_j和影响体积分数v_j (所有样本点的影响体积分数之和为1)。通过确定样本点的影响体积总和V_t，计算每个样本点的影响体积V_j (V_j应该与其密度成反比)，变量V_j的n次方根即为形状参数c_j。

4.3. 数值实验

对函数 $f=\frac{x\sin \left(2\pi x\right)}{1000{\text{e}}^x}$ 在区间 $x\in \left[0,\text{1}\right]$ 上随机地取8个样本点： $\left(0.1,5.3185\times {10}^{-5}\right)$ ； $\left(0.2,1.5573\times {10}^{-4}\right)$ ； $\left(0.49,1.8849\times {10}^{-5}\right)$ ； $\left(0.56,-1.1775\times {10}^{-4}\right)$ ； $\left(0.7,-3.3060\times {10}^{-4}\right)$ ； $\left(0.78,-3.5122\times {10}^{-4}\right)$ ； $\left(0.8,-3.4187\times {10}^{-4}\right)$ ； $\left(0.95,-1.1353\times {10}^{-4}\right)$。

用Gauss函数对其进行插值，分别利用上述方法，选取形状参数，并画出误差图像，比较这些方法的拟合效果( $\alpha =1,n=8,m=101$ )。

4.3.1. 常参数径向基函数中形状参数的选取方法

1) 优化误差确定法

解决如下优化问题：

$\{\begin{array}{l}\min E_{\max }\left(c\right)=\min \left(\max {\left|S\left(x,c\right)-f\left(x\right)\right|}_x\right)\\ Cond\left(A\right)\le {10}^8,c\ge 0\end{array}$ (3)

其中，A为径向基函数插值矩阵，利用Matlab软件编程，得到利用Gauss函数对原函数f(x)作插值的最优形状参数为 $c_{opt}=\text{4}.\text{581}$，且当c取该值时， $\min \left(E_{\max }\right)=2.2405\times {10}^{-5}$。图1、图2显示最优形状参数 $c_{opt}=\text{4}.\text{581}$ 时的插值效果图及误差图。

2) 数值实验确定法

对形状参数c分别取不同的值，用Matlab软件做插值实验，在矩阵非奇异的条件下，计算不同形状参数所对应的插值最大误差，部分结果如表1所示。最大误差值的最小值min(E_max)对应最优形状参数。

Table 1. Shape parameters and corresponding error values

表1. 各形状参数及对应误差值

从表1可以看出，利用Gauss函数对原函数f(x)作插值的最优形状参数为 $c_{opt}=1.73$，且当c取该值时， $\min \left(E_{\max }\right)=1.2261\times {10}^{-7}$。图3、图4显示最优形状参数 $c_{opt}=1.73$ 时的插值效果图及误差图。

图1、图3中插值曲线与原函数曲线基本吻合，说明选取的形状参数都有很好的插值拟合效果，但图4中误差达到10⁻⁷，而图2的误差达到10⁻⁵。说明方法(2)选取的形状参数使得插值误差更小。另外，图4的误差整体波动较图2要小得多，说明方法(1)选取的形状参数使得误差变化幅度大。经计算，方法(1)的插值矩阵条件数 $Cond\left(A\right)=3.0044\times {10}^4$，方法(2)的插值矩阵条件数 $Cond\left(A\right)=5.3486\times {10}^9$，说明方法(1)选取的形状参数使得插值稳定性 [19] 高。图5显示了验证样本点的最大误差随形状参数的变化趋势，显然，整体上，最大误差E_max随形状参数c的增大先减小后增大。但当 $c\in \left[0.\text{7},\text{4}.\text{6}\right]$ 时最大误差有两次较大的波动，最优的形状参数正好对应两次波动的最低点。

Figure 1. Original function f(x) and its interpolation function S(x)

图1. 原函数f(x)及其插值函数S(x)

Figure 2. Interpolation error graph E(x)

图2. 插值误差图E(x)

Figure 3. Original function f(x) and its interpolation function S(x)

图3. 原函数f(x)及其插值函数S(x)

Figure 4. Interpolation error graph E(x)

图4. 插值误差图E(x)

Figure 5. Maximum error E_max with shape parameter c

图5. 最大误差E_max随形状参数c的变化图像

4.3.2. 变参数径向基函数中形状参数的选取方法

1) 公式确定法

根据公式(1)，计算出在各个样本点处的形状参数如表2所示。

利用表2中形状参数的值，确定径向基函数S(x)，编写Matlab程序，对原函数f(x)进行插值，得到插值曲线及误差图像如图6、图7所示。

Table 2. Shape parameters of each sample point

表2. 各样本点的形状参数

2) 样本点局部密度确定法

根据武泽平给出的方法步骤，计算出在各个样本点处的形状参数如表3所示。

利用表3中形状参数的值，确定径向基函数S(x)，编写Matlab程序，对原函数f(x)进行插值，得到插值曲线及误差图像如图8、图9所示。

Table 3. Related data to shape parameters of each sample point (take Vt = 2)

表3. 各样本点形状参数相关数据(取V_t = 2)

Figure 6. Original function f(x) and its interpolation function S(x)

图6. 原函数f(x)及其插值函数S(x)

Figure 7. Interpolation error graph E(x)

图7. 插值误差图E(x)

Figure 8. Original function f(x) and its interpolation function S(x)

图8. 原函数f(x)及其插值函数S(x)

Figure 9. Interpolation error graph E(x)

图9. 插值误差图E(x)

图6、图8中，当 $x>0.09$ 时，插值曲线与原函数曲线基本吻合，说明两种方法选取的形状参数在除起始端点的很小区间 $x\in \left[0,0.09\right]$ 外，都有比较好的插值拟合效果，方法(1)的 $E_{\max }=8.7201\times {10}^{-5}$ (图7)，而方法(2)的 $E_{\max }=2.3055\times {10}^{-4}$ (图9)，即优先选择方法(1)确定变参数径向基函数的形状参数，可使得插值误差更小。另外，除插值区间的两侧端点 $x\in \left[0,0.2\right]$ 及 $\left[0.9,1\right]$ 外，图7、图9的误差波动都很小。说明两种方法选取的形状参数使得误差变化幅度都很小。经计算，方法(1)的插值矩阵条件数 $Cond\left(A\right)=1.8036\times {10}^{17}$ ；方法(2)的插值矩阵条件数 $Cond\left(A\right)=1.7653\times {10}^{17}$ ；说明两种方法的插值稳定性相当。

4.3.3. 方法的比较与改进

对上述数值实验的结果进行比较发现，对于要求稳定性的热传导等问题，应该优先选用常参数径向基函数的方法(1)。而对于要求精确性的飞机外形设计等问题，应该优先选用常参数径向基函数的方法(2)。变参数径向基函数的两种方法都能保证插值误差波动较小，但插值精度都相对较低。从图6、图8可以看出，导致这种结果的原因是，在区间端点附近，插值误差明显很大。据此，我们需要提高端点附近的插值精度。

拉格朗日法是常用的数据插值方法，其公式结构整齐紧凑，在不增加插值节点的情况下，有比较好的拟合效果。因此，我们用下述分段函数(4)作为插值函数，结合Matlab软件编程拟合f(x)。得到插值曲线如图10、图11所示。

$\{\begin{array}{ll}{L}_1\left(x\right)={\displaystyle \underset{j=1}{\overset{r_1}{\sum }}{P}_j\left(x\right)f_j}\hfill & x\in \left[0,0.2\right]\hfill \\ S\left(x\right)={\displaystyle \underset{j=1}{\overset{n}{\sum }}{\lambda }_j}\varphi \left(\Vert x-x_j\Vert \right)\hfill & x\in \left[0.2,0.9\right]\hfill \\ L_2\left(x\right)={\displaystyle \underset{j=1}{\overset{r_2}{\sum }}{P}_j\left(x\right)f_j}\hfill & x\in \left[0.9,1\right]\hfill \end{array}$ (4)

其中， $P_j\left(x\right)={\displaystyle \underset{k=1,k\ne j}{\overset{r}{\prod }}\frac{x-x_k}{x_j-x_k}}$ 为拉格朗日基本多项式(r为插值节点的个数， $r_{\text{1}}=\text{5}$， $r_{\text{2}}=\text{3}$ )。

Figure 10. The improved method (1) interpolation effect map

图10. 改进后的方法(1)插值效果图

Figure 11. The improved method (2) interpolation effect map

图11. 改进后的方法(2)插值效果图

图10、图11中，插值曲线与原函数曲线基本吻合，说明改进后的两种方法在整个区间 $x\in \left[0,\text{1}\right]$ 上都有很好的插值拟合效果。通过计算，改进后的方法(1)的拟合优度R²从0.9299变为0.9915；改进后的方法(2)的拟合优度R²从0.8165变为0.9803，这比常参数径向基函数中方法(1)的拟合优度0.9686都高，说明此改进方法可以有效提高插值精度。在实际工程问题中，我们同样可以使用上述方式对散乱数据进行处理，得到更好的插值拟合效果。

5. 结论

本文结合径向基函数插值误差理论，对现有的形状参数的选取方法进行归纳总结，并通过数值实验，比较各方法的优缺点。针对变参数径向基函数区间端点附近插值精度不高的问题，利用拉格朗日插值法进行改进，Matlab软件的运行结果表明，改进后的插值方法能够得到更高的精度。

基金项目

国家自然科学基金项目资助(11601151)。