由Hardy^[15-16]提出的基本方案可知，任意函数f可以写成N个连续可微的基础函数的组合：

式中，g_j为核函数；α_j为待定系数；N为核函数个数。核函数一般可选为：

式中，(x_j, y_j)为节点坐标；δ为光滑因子；k决定了核函数的形状。陶本藻等^[19]研究表明，在推估中k值可取为0.5，即采用正双曲面型，δ取0。在推估中，将研究区域视为一数学表面，根据多面函数的原理可利用多个核函数的组合来逼近该数学表面。利用已知离散点数据计算得到所有待定系数α_j的估值，即可利用式(1)推估得到所有未测点的估值。

标准的最小二乘配置函数模型为：

式中，B为设计矩阵； l为观测向量；Δ为观测噪声向量；G_y为倾向部分； x为随机信号向量。最小二乘配置的原理是通过已测点数据和已测点信号与未知点信号之间的协方差来推估未知点的值，其核心问题在于已测点信号与未知点信号之间协方差的确定。式(3)的最佳线性无偏解为^[6]：

式中，l_x为信号的均值，通常取l_x = 0；D_x为已测点信号之间的协方差矩阵；D_x'x为推估点信号与已测点信号之间的协方差阵；$\mathit{\boldsymbol{\hat x}}{\rm{'}}$为待推估的信号值。

在最小二乘配置的框架下，本文将多面函数法融入到最小二乘配置的数学模型中，即用多面函数法对区域的倾向项进行估计，然后利用拟合推估方法对空间随机变化项进行补偿，综合两种方法的优点以达到更精确的推估效果。改进最小二乘配置的函数模型可以表示为：

式中，l为n×1观测向量；Q为n×m核函数阵，m为节点数目；a为m×1待定系数向量；Δ为n×1误差向量。与式(4)同理, 得到待定系数a以及信号x的最佳线性无偏解：

式中，l_x为信号的均值，通常取l_x =0；D_x为已测点信号之间的协方差矩阵；D_x'x为推估点信号与已测点信号之间的协方差阵；$\mathit{\boldsymbol{\hat x}}{\rm{'}}$为待推估的信号值。

由式(4)可知，多面函数待定参数的估计需要已知的信号协方差阵，而信号的协方差阵则需要首先从观测值中将趋势项去除才能统计拟合得到。本文使用陈家军等^[21]提出的多次迭代的方法解决这一矛盾，获取稳健的经验协方差函数和多面函数待定系数值。对于地表垂直形变，可以近似认为是各项同性的，用高斯模型来拟合经验协方差值^[1]。高斯模型的表达式如下：

式中，d为两点间的距离；C(0)与γ为待定系数, C(0)表示信号的方差, γ为信号的空间相关尺度。改进方法具体的计算步骤如下：

1) 初始化数据协方差阵为单位阵(D_x= E)，即C(0)_i=1，γ_i=0，i=0；

2) 利用协方差函数模型参数C(0)_i和γ_i得到协方差阵D_x。由于B、Q已知，通过式(6)得到待定系数的估值a_i；

3) 根据式(5)去除观测数据l中的趋势项Qa_i，得到信号x_i；

4) 根据文献[1]中的方法对信号x_i进行实验协方差统计，得到不同距离段的协方差值(d, C(d))，继而利用最小二乘拟合得到协方差函数模型参数C(0)_i+1和γ_i+1，更新协方差阵D_x；

5) 将更新协方差阵D_x代入式(6)得到新的待定系数值a_i+1；

6) 判断条件：|C(0)_i+1-C(0)_i| < 10^-6且|γ_i+1-γ_i| < 10^-6，不满足则i=i+1，回到第2步。否则继续；

7) 根据待定系数a_i+1，利用式(6)得到未测点的信号值$\mathit{\boldsymbol{\hat x}}{\rm{'}}$，最后综合未测点的趋势项与信号值得到最终估值：$\mathit{\boldsymbol{\hat l}}' = \mathit{\boldsymbol{\hat x}}'\mathit{ + }\mathit{\boldsymbol{Q}}{\mathit{\boldsymbol{a}}_{i + 1}}$。

为了验证改进方法的有效性，本文首先通过一个模拟的地震滑动数据来比较改进方法与常规方法的推估效果。考虑一个位于均匀的，各项同性的弹性半空间中的倾向断层，其长度为75 km，宽为60 km，走向角为70°，倾角为15°，模拟的滑动分布见图 1。

图  1  模拟的断层滑动分布图

Figure 1.  Fault Geometry and Slip Distribution Used in the Synthetic Test

利用Okada公式^[22]计算模拟断层滑动引起的覆盖研究区域的374个均匀格网点(格网点X与Y方向间隔均为4 km)的地表垂直形变，分布模式见图 2。此外，分别按1 km和2 km进行采样得到的结果与4 km采样的结果类似，鉴于篇幅所限，文中主要展示了点位最为稀疏的采样(即4 km)。

图  2  模拟垂直形变值与数据点位

Figure 2.  Vertical Surface Displacement and

从生成的374个已知形变值的点中随机选择200个点作为计算点，另外174个点作为外部检核点。对计算点加入均值为0 mm，标准差为3 mm的高斯误差。多面函数计算中，选取了全部计算点数作为节点数，且均匀覆盖研究区域。为了进行比较，使用趋势项为线性多项式的最小二乘配置与多面函数法进行了计算。实验中对计算点的随机选取共模拟了500次，几种不同方法的外部检核残差统计见表 1。通过对计算结果以及拟合残差均方根(RMS)统计结果(表 1)的比较分析可以发现：多面函数法得到的结果外部检核点的RMS为最大，以线性多项式为趋势的常规最小二乘配置的RMS次之，改进方法的RMS是最小的，表明改进方法得到的结果最为精确。通过比较这几种方法RMS的标准差可知，改进方法受计算点选取的影响最小，得到的结果更加稳定。

为了评估不同噪声水平对推估结果的影响，本文还对模拟数据加入均值为0 mm，标准差分别为5 mm和10 mm的高斯噪声进行试验。当加入标准差为5 mm的高斯噪声后，多面函数法得到RMS均值为13.9 mm，常规最小二乘配置为8.7 mm，改进方法为7.5 mm。当加入标准差为10 mm的高斯噪声后，多面函数为22.8 mm，常规最小二乘配置为10.5 mm，改进方法为10.1 mm。随着噪声方差值的增加，多面函数得到的结果RMS增加最迅速，而改进方法与常规最小二乘配置得到的RMS趋于一致。这是因为随着外部噪声水平的增加，当噪声的数值与信号数值相当时(图 3)，无法有效地将信号从噪声中提取出来，因此改进方法无法发挥其优势。

图  3  噪声均值为0 mm, 标准差为3 mm时利用改进算法得到的信号项

Figure 3.  Signal Distribution Estimated by Refined Least Squares Collocation with 0 mm Mean and 3 mm Standard Deviation Gaussian Error and Trend Distribution Estimated

2009年4月6日，意大利中部Abruzzo地区的L’Aquila发生Mw 6.3级地震, 致使300人死亡。为了对改进方法的效果进行检验，本文使用了意大利L’Aquila地震InSAR和GPS的同震形变数据来进行计算分析。Envisat卫星Track 129的干涉经降采样得到的1 282个点沿视线向的形变值来自文献[23]，56个GPS水平形变值与26个GPS垂直形变值来自文献[24]，点位分布见图 4。InSAR点视线向形变中误差为6.8 mm，GPS水平形变中误差为3.9 mm，GPS垂直形变中误差为7.8 mm。由于InSAR数据得到的是沿视线向的形变值，并不能得到垂直方向的形变。通过InSAR成像几何可以得到视线向形变值与东、北和垂直向形变值的转换关系：

图  4  L’Aquila地震震中区的InSAR点与垂直GPS点分布图

Figure 4.  Distribution of InSAR and Vertical GPS Observations

式中，[d_e d_n d_u]为地面点在东、北和垂直向的形变量；α_h为雷达飞行方向的方位角；θ为雷达的侧视角。等式右边第一项表示视线向形变值方向单位矢量在东、北和垂直向的分量。利用GPS东、北向离散形变值插值得到所有InSAR点东、北方向形变值，再利用式(8)即可将雷达视线向形变d_los转换到垂直向形变d_u。

为了检核推估效果，本文将1 282个InSAR数据点作为已测点，将26个GPS点垂直向形变值作为外部检核。计算中，由于InSAR点数量较多，且空间分布较为接近，导致信号的协方差阵的复共线性严重，从而导致配置解的不稳定。为了保证得到稳定的结果，需要对其计算模型进行正则化修正。由于推估点与计算点是同类型数据，即设计矩阵为单位阵B = E。根据文献[25]的方法可将式(6)改写为：

式中，K = D_x + D_Δ。令：

假设G为稳定的，则只需要得到稳定的q，即可确保可以获得稳定的待估量${\mathit{\boldsymbol{\hat a}}}$。将式(10)改写为：

对式(12)引入Tikhonov正则化法，其正则化准则为^[26]：

式中，α为正则化参数；‖q为q‖范数。通过式(13)可得q的正则化解：

通过式(14)可知协方差阵K的逆的稳定形式可以用(K^TK+αI)^-1 K ^T表示。考虑G为不稳定的，令：

则式(9)可改写为：

同理，要得到稳定的待估量${\mathit{\boldsymbol{\hat a}}}$，则必须得到稳定的H。通过相似的推导可以得到${\mathit{\boldsymbol{\hat a}}}$的正则化解：

式中，β为正则化参数。

利用上述方法进行正则化修正得到信号的正则化解为：

上文中正则化参数的估计可以采用Hansen^[27]提出的L曲线法求得，其核心是定位L曲线上曲率最大的那个点的正则化参数值。改进方法的推估流程与§2中的流程一致，只是用式(17)~式(19)替代式(6)。

为了对计算结果进行比较分析，文中还同时使用了正则化修正的多面函数和以线性多项式为趋势项的最小二乘配置来计算。不同方法的外部检核残差统计结果(表 2)表明, 改进方法相对于另外两种常规方法，各项指标都是最优的。实测数据计算得到了与模拟数据计算相同的结论，这表明将多面函数引入后能够更真实地反映出该区域系统性的形变，将系统部分去除得到的随机场也更加平稳，因此能够通过推估的方法得到准确的信号值。但是由于InSAR数据处理中，视线向形变值残留了轨道误差、DEM误差以及大气误差等的影响，而且GPS的垂直向精度也不是很高，这会导致两者之间残差RMS较大。

在利用最小二乘配置对非平稳空间场进行推估时，趋势去除法是最为广泛的策略，但是在实际应用中由于趋势项数学模型选择的不准确导致无法完全去除区域的趋势。本文在分析多面函数法的优良数学特性后，提出了一种改进的最小二乘配置方法来解决上述问题。该方法利用多面函数法来准确拟合区域的倾向项，利用迭代的方法获得稳定的待定系数值与协方差函数系数值，最后综合趋势项与信号项得到了待估点的估值。

本文首先通过对一组模拟数据进行计算验证了改进方法的有效性。改进方法得到的结果要优于多面函数法以及以线性多项式为趋势的最小二乘配置法。此外，还利用该方法对2009年意大利L’Aquila地震的同震形变数据进行计算。针对InSAR数据量多，且分布密集的特点，引入了正则化方法对改进方法解的表达式进行修正以解决病态对结果带来的影响。文中使用了多面函数与常规最小二乘配置法作为对比，通过外部检核GPS点残差的统计结果表明，改进方法同样能够得到最佳的推估值。

本文所使用的协方差函数模型只考虑了各项同性的情况，如果区域内信号有明显的各向异性，那么在进行推估时需要引入各向异性的协方差函数模型。由于改进方法是基于多面函数的，那么在利用多面函数进行推估时遇到的问题，改进方法也同样会遇到，即节点数量以及位置的选择、平滑因子的选择和核函数形状的选择。至于用多面函数法进行趋势拟合后对信号的影响，有待于更进一步研究。