1 引言

最小二乘配置最初是指组合各种不同类型的数据资料去研究地球形状和重力场的一种数学方法(Borre，2006).Moritz(1972)对最小二乘配置进行了进一步的深入研究.地壳运动不仅仅存在整体的趋势形变，同时存在局部复杂的点位变化信息.传统的三角形内插法(Eringen，1980)、有限元法(蒋志浩等，2009)、多面函数法(黄立人等，1993；刘经南等，2001)等方法均不能很好的表达局部点位变化信息.为顾及形变场的局部点位变化信息，藤井洋一郎首次用最小二乘配置方法处理地壳形变分析问题.El-Fiky和Kato(1998)利用最小二乘配置方法获得研究区域的形变信号.Wu等(2011)通过不断加大观测数据噪声探讨了最小二乘配置、多面函数以及球谐函数的抗差性能，通过试验得出最小二乘配置的抗差性较为良好.陈光保等(2010)针对形变场存在各向异性的问题，首先利用传统最小二乘配置方法给出地壳形变监测点趋势项估值和梯度，以梯度调整监测点之间的距离，再通过迭代的方法获得更接近实际情况的协方差函数.杨元喜等(2011)针对最小二乘配置易受到先验协方差函数误差影响的问题，利用方差分量估计构造自适应因子调整随机信号的先验方差阵，自适应地调整观测向量与信号向量的先验权比，使得随机信号和观测误差对参数估计的影响得到平衡.刘念(2001)针对最小二乘配置在重力异常内插中易受到观测值中的粗差点的影响，提出协方差函数抗差拟合方法以达到抗差的目的.柴洪洲等(2009)考虑到地球物理模型不受粗差的影响，利用研究区域的地球物理模型进行观测值中心化，以文献(江在森等，2003)的方法进行协方差函数参数k的拟合，通过这两种方法达到抗差的目的.赵丽华等(2011)针对研究区域可能存在的明显的局部变形或观测异常的测站，提出高崩溃污染率抗差估计法进行处理.邹镇宇等(2015)提出一种改进的拟准检定初选指标方法去剔除观测数据中的粗差点，然后再进行协方差函数的拟合以达到抗差的目的.目前，最小二乘配置已被广泛应用于重力异常(Jarmołowski，2015)、高程异常(Yang et al., 2009)、坐标转换(You and Hwang, 2006)以及地壳形变分析(柴洪洲等，2009；江在森和刘经南，2010；杨元喜等，2011；赵丽华等，2011)中.但最小二乘配置结果的合理性，关键在于经验协方差函数的拟合.当观测数据存在异常点或是粗差点时，常规的拟合方法给出的协方差函数不能精确的表征其统计特性，且观测数据中的粗差点或异常点导致以经验方式给出的权重所确定的方差并不准确，把不具抗差性的协方差函数计算的信号协方差及不准确的观测值方差阵代入到最小二乘配置计算公式中，获得的最小二乘配置参数解并不准确，最终给出点位拟合推估结果的合理性值得商榷.

针对上述问题，本文在文献(刘念，2001；柴洪洲等，2009；赵丽华等，2011；谢曦霖等，2017)的基础上，提出两种抗差最小二乘配置迭代方法，并利用这两种方法与传统的抗差最小二乘配置方法、常规最小二乘配置方法分别对芦山地震以及意大利L′Aqulia地震数据进行处理分析，验证本文方法的有效性.

2 地壳形变分析模型的构建 2.1 传统最小二乘配置方法

地壳形变分析包含趋势项以及随机项两部分，根据最小二乘配置原理，建立模型为(Yang，1992；崔希璋等，2009)：

(1)

式中，L为观测向量；GY为趋势项，其中Y为非随机参数，G为对应的系数矩阵；B_ZZ为随机项，其中Z为随机信号，包含了已测点信号X和未测点信号X′，B_Z=[B, 0]^T，B为q阶单位阵，0为l－q阶零矩阵；n为观测噪声.

根据广义最小二乘准则，得到非随机参数和随机信号计算表达式为(崔希璋等，2009)：

(2)

式中，D_X为已测点信号协方差，D_n为观测误差自协方差，D_X′X为未测点信号与已测点信号的协方差阵.其中，D_X和D_X′X由协方差函数给出.

相关研究表明，在地壳形变分析中，高斯函数作为协方差函数较为合适(江在森和刘经南，2010)，因此本文选取高斯函数作为经验协方差函数，协方差函数模型为

(3)

式中，d表示两点间的距离，C₀为距离为0的方差，k为待定系数.

假定协方差只与任何两点间的距离有关，跟方向与点位无关.用以下方法计算C₀和k，首先利用最小二乘法计算出非随机参数，公式为

(4)

根据式(1) 去除趋势项，得到信号X为

(5)

对信号进行统计，得到C_L(0) 和C_d为

(6)

(7)

式中，n_d为任何两个距离为d的点对的总个数.分别把C_d和对应的距离d代入式(3)，利用最小二乘方法得到协方差函数参数C₀和k.

2.2 基于观测值中位数初值的抗差最小二乘配置迭代解法 2.2.1 初始观测值权阵

在观测数据不含异常点或粗差时，利用常规的最小二乘配置获得研究区域的运动模型具有高度可靠性.顾及到观测数据存在异常点或粗差点时，常规最小二乘配置解并非是一个最优解.同时传统的抗差方法中，给出的参数初值一般为最小二乘解.由于观测数据存在异常点或粗差点，最小二乘解易受到粗差点或异常点的影响，使得解产生较大的偏离，从而导致后续去除趋势项后，给出的信号不准确，等价权也随之失效(赵丽华等，2011).相关研究表明，中位数的抗差性能较好，其崩溃污染率为50%(刘念，2001).因此，采用中位数进行观测值中心化具备良好的抗差效果，利用中位数进行中心化可表达为

(8)

根据给出的观测值中心化，计算初始方差因子为(Yang, 1999; 赵丽华等，2011)：

(9)

后续选择合适的等价权方案进行观测值的抗差计算.值得注意的是，本文方法在此处暂不需要进行迭代计算.在更新观测值权阵后，初始化协方差函数，计算出最小二乘配置的初始解.

2.2.2 协方差函数的抗差拟合

获得准确的经验协方差函数是求解最小二乘配置问题的关键.由上述2.2.1节得到最小二乘配置的初始解后，把初始解代入式(5) 计算出信号X，并根据信号X得到C_d，计算方法为：

设有s+1个等间距的d_l=l×μ，式中l=0, 1, 2, …s，μ为分段的间隔距离.由于数据点分布不规则，根据实际情况选取一个常量Δd，利用满足|d_ij－d_l| < Δd的数据点对的信号X计算经验协方差的函数值C_dl，为了进一步减弱观测数据中的粗差点或异常点对后续计算的经验协方差参数解的影响，通过应用等价权抗差估计得到观测数据的一系列权因子，C_dl的计算公式为

(10)

式中，s_l为满足|d_ij－d_l| < Δd点对的总对数.

对式(3) 的高斯函数两边进行求对数处理，则有：

(11)

计算式(11) 的最小二乘解，得到经验协方差函数的抗差估计.基于观测值中位数初值的抗差最小二乘配置迭代方法具体计算步骤如下：

(1) 初始化协方差函数，给出其参数初始值C₀⁽⁰⁾=1，k⁽⁰⁾=0，利用观测值中位数进行观测值中心化，并应用合适的等价权方案更新观测值方差阵，给出最小二乘配置初始解.

(2) 由得到的最小二乘配置初始解，去除趋势项，应用等价权方案更新观测值方差阵，进行信号协方差统计，得到新的协方差函数参数估值C₀⁽ⁱ⁾和k⁽ⁱ⁾.

(3) 由更新的协方差函数参数估值以及观测值方差阵，计算信号协方差D_X.把D_X及更新后的观测值方差阵D_n代入式(2)，计算出最小二乘配置解，根据步骤(2) 的流程更新协方差函数参数估值以及观测值方差阵.

(4) 重复步骤(3)，直至满足迭代条件‖C₀⁽ⁱ⁺¹⁾－C₀⁽ⁱ⁾‖ < ε&&‖k⁽ⁱ⁺¹⁾－k⁽ⁱ⁾‖ < δ(文中设为ε=δ=10^－10)时，迭代结束，否则回到步骤(3)，继续迭代计算，最终获得稳健的协方差函数参数估值及最小二乘配置解.

2.3 基于中位参数法的抗差最小二乘配置迭代解法

从式(1) 中可以看出，有q组观测方程，p个非随机参数，在保证可以进行精度评定的基础上，能够获得C_q^p+1组最小二乘配置解.其中，第i个参数解组成C_q^p+1维向量，公式为

(12)

根据式(12) 可得到第i个参数解的中位数median()，则非随机参数的中位数为

(13)

计算各组解与非随机参数中位数的差值，即：

(14)

取对应差值向量的最小范数对应的非随机参数解作为最后的参数解，即：

(15)

相关研究表明，中位参数法的崩溃污染率可达50%(杨玲等，2011).基于中位参数法的抗差最小二乘配置迭代解法步骤如下：

(1) 初始协方差函数，给出其参数初始值C₀⁽⁰⁾=1，k⁽⁰⁾=0，计算出初始信号协方差D_X⁽⁰⁾.

(2) 利用中位参数法给出最小二乘配置初始解，去除趋势项，应用合适的等价权方案更新观测值方差阵，进行信号协方差统计，得到新的协方差函数参数估值C₀⁽ⁱ⁾和k⁽ⁱ⁾.

(3) 根据步骤(2) 更新的协方差函数模型及观测值方差阵D_n⁽ⁱ⁾，计算出信号协方差D_X⁽ⁱ⁾，由中位参数法给出最小二乘配置解，根据步骤(2) 流程更新协方差函数参数估值以及观测值方差阵，直至满足迭代条件：‖C₀⁽ⁱ⁺¹⁾－C₀⁽ⁱ⁾‖ < ε&&‖k⁽ⁱ⁺¹⁾－k⁽ⁱ⁾‖ < δ(文中设为ε=δ=10^－10)时，迭代结束，否则回到步骤(3)，继续迭代计算，最终获得稳健的协方差函数参数估值及最小二乘配置解.

3 地震实例分析 3.1 芦山地震

2013年4月20日，四川芦山发生7.0级地震，这次地震造成了人们的巨大损失.本文所用的33个GPS垂直形变以及相应的中误差数据来源于文献(Jiang et al., 2014)，如图 1所示.其中随机选取20个GPS点作为已测点，剩余13个GPS点作为外部检核点.设计4种计算方案对数据进行处理.

方案1：常规最小二乘配置方法.

方案2：传统抗差最小二乘配置方法，即文献(王海栋等，2010)方法.

方案3：基于中位参数法的抗差最小二乘配置迭代方法.

方案4：基于观测值中位数初值的抗差最小二乘配置迭代方法.

值得注意的是，经过试验表明，方案2、3和4的等价权函数取huber等价权方案时较为合适，根据下式确定观测值权阵及方差阵，公式为

(16)

c一般取值2.四种方案已测点和外部检核点的残差统计表分别如表 1和2所示，协方差统计如图 2所示，已测点与外部检核点残差分布图如图 3所示.

为了进一步验证本文提出的方法的抗差效果，遵循粗差点在观测数据中一般不超过10%的原则，在已测点中的其中两个点分别混入粗差，并不断增加该两个点的粗差大小，根据粗差大小的变化，四种方案的已测点和外部检核点的RMS变化如图 4所示.

3.2 意大利地震

2009年4月6日，意大利L′Aqulia发生了6.3级地震，这次地震给当地人们造成了极大的伤害，包括300多人死亡，1000多人受伤，2万多人流离失所，无家可归(温扬茂等，2012).鉴于最小二乘配置是用于空间域离散点观测信号的拟合推估方法，基于被观测信号理论上是空间连续分布的模型假设的前提下.由于同震位移场本身在发震断层处存在显著的不连续分布，利用最小二乘配置处理发震断层边缘的点必然会产生较大的模型残差，因此用最小二乘配置整体处理同震位移场数据并不合适.本文通过选取某一部分连续分布区域的Envisat卫星Track129干涉得到的650个InSAR点作为实验数据，其形变值中误差均为6.8 mm.随机选取600个点作为已测点，其余50个点作为外部检核点，如图 5所示.

由于我们需要的是垂直方向的形变值，而实际获取的是InSAR视线向的形变值，根据垂直形变与视线向形变之间的转换关系，将650个视线向形变值转换成垂直方向的形变值，具体转换公式为

(17)

式中，d_ui为第i点视线向形变值分解到垂直方向的形变值，d_losi为第i点雷达视线向形变值，θ_i为第i点雷达波入射角.

然后，根据式(17) 和误差传播定律，得到各点垂直形变的中误差.设计3种计算方案对意大利L′Aqulia地震的数据进行处理.

方案1：常规最小二乘配置方法.

方案2：传统抗差最小二乘配置方法.

方案3：基于观测值中位数初值的抗差最小二乘配置迭代方法.

值得注意的是，经过试验表明，方案2和3的等价权函数采用IGGI方案时较为合适，根据下式确定观测值方差阵，公式为

(18)

在这里，k₀取值1.5.

利用三种方案进行计算，已测点和外部检核点残差统计表分别如表 3和4所示，协方差统计结果如图 6所示，残差统计图如图 7所示.

同样的，为了进一步验证本文提出的方法的抗差效果，遵循粗差点在观测数据中一般不超过10%的原则，在已测点中的其中60个点分别混入粗差，并不断增加这60个点的粗差大小，根据粗差大小的变化，三种方案的已测点和外部检核点RMS变化如图 8所示.

3.3 结果分析

从表 1中可以看出，抗差方法均比常规方法的拟合推估效果要好.在与传统的抗差方法的比较中，本文提出的基于观测值中位数初值的抗差最小二乘配置方法拟合效果最好.在芦山地震实例中，在标准差上，基于观测值中位数初值的抗差最小二乘配置方法与传统抗差方法相当、比基于中位参数法的抗差最小二乘配置法提高了7.6%，比常规最小二乘配置方法提高了76.2%.另一方面，四种方案在推估效果上，从表 2可以看出基于观测值中位数初值的抗差最小二乘配置方法推估效果较优，与传统抗差方法推估效果相当.在标准差上，基于观测值中位数初值的抗差最小二乘配置方法比基于中位参数法的抗差最小二乘配置法提高了35.7%、比常规最小二乘配置提高了31.1%.综上所述，基于观测值中位数初值的抗差最小二乘配置方法在拟合和推估效果上均较为良好，这是因为常规最小二乘配置方法受到观测数据的粗差点或异常点的影响较大，导致计算的结果与实际情况存在较大的偏差；基于观测值中位数初值的抗差最小二乘配置方法应用合适的等价权抗差估计，并通过迭代的形式自适应性地调整观测值方差阵，以抗差的方法计算协方差函数参数C₀和k，高效的利用了所有的观测数据进行计算，而基于中位参数法的抗差最小二乘配置方法虽然具备一定的抗差性，但实际上中位参数法最终给出的最小二乘配置解仅仅是由p+1个观测数据给出的解，缺少某些特征点参与最后的计算，且该方法计算量较大，导致了计算效率较低，完成整个计算过程耗时2339 s，传统抗差方法耗时0.09 s，常规最小二乘配置方法仅耗时0.03 s，观测值中位数初值的最小二乘配置方法仅耗时0.07 s.

考虑到基于中位参数法的抗差最小二乘配置方法在计算效率以及抗差效果相比于基于观测值中位数初值的抗差最小二乘配置较差的情况，在意大利L′Aqulia地震实例中，仅设计常规最小二乘配置、传统抗差最小二乘配置方法和基于观测值中位数初值的抗差最小二乘配置方法这三种方案进行计算.在拟合精度上，从表 3中可以看出，基于观测值中位数初值的抗差最小二乘配置方法的标准差与传统抗差最小二乘配置方法相当，比常规最小二乘配置精度提高了19.9%.在推估精度上，从表 4可以看出，基于观测值中位数初值的抗差最小二乘配置方法的标准差与传统抗差最小二乘配置方法相当，比常规最小二乘配置精度提高了24.6%.在计算效率上，基于观测值中位数初值的抗差最小二乘配置方法耗时10.1 s、传统抗差最小二乘配置方法耗时9.8 s、常规最小二乘配置方法耗时9.5 s.进一步验证基于观测值中位数初值的抗差最小二乘配置的良好抗差效果.

图 2和图 6分别为芦山地震和意大利L′Aqulia地震的协方差统计图.值得注意的是，芦山地震的点距单位为km，而意大利L′Aqulia地震的点距单位为m.这样选取的原因是在意大利L′Aqulia地震实例中，点距以m为单位时观测值中位数初值的抗差最小二乘配置方法的计算效率更为理想，且点距单位的不同对协方差函数的参数k存在着这样的一个关系：k_km=1000 k_m，式中k_km为点距以km为单位的协方差函数参数k的估值，k_m为点距以m为单位的协方差函数参数估k值；而对于另一个协方差函数参数C₀，点距单位对其没有直接影响.从图 2和图 6可以看出，芦山地震的相关尺度为375 km，而意大利L′Aqulia地震的相关尺度为31.5 km；在连续分布场中，空间位置靠近的点的观测值理论上较为接近，观测点之间观测信号随着点距的增大而衰减，若不具备这样的特性，则表明数据中的信噪比较低，难以识别出有效的变形信号(江在森和刘经南，2010).从图 2c可以看出，其信噪比较低，基于中位参数法的抗差最小二乘配置方法难以识别出有效的变形信号，而其他方案的观测点之间观测信号随着点距的增大而衰减较为明显，能够识别出有效的形变信息.而从图 6可以看出，意大利L′Aqulia地震相关尺度仅为31.5 km，表明依然能从观测噪声中提取有效的形变信息.不同方案获得的协方差函数参数估计存在差异，抗差方法之间也略有差别.这是因为在协方差函数拟合过程中，不同的抗差方法均可以抑制粗差在计算过程中的影响，同时能够在一定程度上抑制异常运动点的影响，而基于观测值中位数初值的抗差最小二乘配置方法通过迭代计算给出更为合理的协方差函数参数估值，给出的协方差函数更为客观的反映点位之间的相关性.

从图 3和图 7可以看出，基于观测值中位数初值的抗差最小二乘配置方法与最小二乘配置方法的已测点和外部检核点残差相比常规最小二乘配置方法较小，反映了抗差方法拟合推估的形变值更接近实际情况，在计算过程中能够很好的抑制粗差或异常点的影响，标准差明显的降低.

为了进一步检验抗差方法的抗差性能，在芦山地震和意大利L′Aqulia地震均人为的进一步引入了粗差，两个地震实例的标准差变化图分别如图 4和图 8所示.通过粗差大小变化及其对应的标准差变化情况可以看出，常规的最小二乘配置方法抗差能力较差，传统抗差方法虽然具备一定的抗差能力，但在粗差增大到某一程度后，其抗差能力失效，而基于观测值中位数初值的抗差最小二乘配置方法抗差性能最好，当粗差增大到某一程度时，依然保持良好的抗差性能.

4 结论

在地壳形变分析中，观测数据不可避免的存在粗差点或异常点的情况，常规最小二乘配置方法不具备抗粗差的能力，在处理含有粗差点或异常点的观测数据时，给出的解并非是一个最优解，不能正确的反映研究区域的实际情况；本文提出两种抗差最小二乘配置方法，分别利用观测值中位数和结合中位参数法，均在给出协方差函数参数初始值下，通过迭代计算的方法，自适应地调整观测值方差阵，以抗差方法计算协方差函数参数，最终获得稳健的协方差函数参数估计及最小二乘配置解.通过地震实例表明，相对于传统的最小二乘配置方法和抗差最小二乘配置方法，本文提出的基于观测值中位数初值抗差最小二乘配置方法具备良好的抗差性能，能够获得较好的拟合推估效果.

抗差方法在一定程度上提高了拟合推估的精度，但在提高解算精度的同时，降低了计算效率，特别是基于中位参数法的抗差最小二乘配置方法计算效率明显较差，而基于观测值中位数初值的抗差最小二乘配置方法很好的在精度与计算效率之间得到了平衡，与传统的抗差最小二乘配置的计算效率相当，计算效率在可接受的范围之中.

本文提出的抗差方法同样适用于重力异常、高程异常和坐标转换等研究领域.

致谢

感谢审稿专家对本文提出的宝贵意见，感谢武汉大学温扬茂博士提供意大利L′Aqulia地震数据.本文芦山地震数据来源于文献(Jiang et al., 2014)，在此一并致谢.