IGG III选权迭代法的重力网抗差估计

李 琦，武 威

（1.西安测绘总站，陕西 西安 710054）

重力加密测量为重力基准的建立、油气资源的勘探、地质构造及地下水资源的开发利用提供了重要的基础数据和手段[1-3]。我国2000国家重力基本网已经使用了近20 a，随着经济社会的发展部分点位遭到破坏[4]，以及地球内部及表面的质量调整，进行加密重力测量弥补2000国家重力基本网在部分地区的不足就尤为必要。相对重力仪的非线性零点漂移及格值系数变化、重力场时变信号改正不尽完善，另外野外作业人员操作失误等因素造成观测数据中可能包含有较大误差、甚至是粗差[5-6]。同时，加密重力测量也是一项费时、费力的精细工作，获取的数据资料珍贵[7]，故对加密重力测量数据进行高精度处理就显得非常关键。

在2000国家重力基本网平差中，采用“弱基准”的方式，对绝对观测量及相对观测量将赋以适当的权进行平差，以减弱两者精度不匹配的影响[8]。在沿海省市的重力网数据处理中，用抗差等价权来调整相对重力的权、重力仪参数的取舍[9]。

抗差估计的核心是按照权值将观测值进行正常观测值、可利用观测值和粗差观测值，权值大小不变的为正常观测值。对残差较大的观测值进行降权处理，为可利用观测值，而权值降为零的观测值为粗差观测值，从而避免粗差对结果的影响[10]。最小截断二乘法（least trimmed squares，LTS）的抗差方案中，处理含有粗差的大样本量数据处理效率不高[11]。最小二乘估计不具有抗干扰性[12-14]，抗差最小二乘估计通过等价权把抗差估计与最小二乘法结合起来。等价权的选取有L1法、German-McClure法、Danish法、Huber法IGG法和IGG III等多种方法[15-16]。为此，以某区域加密重力网为例，在平差中采用IGG III的权因子函数[16]，进行权值迭代平差，获取可靠的参数估值。

1 重力网的抗差估计方案

依据绝对重力观测值，其绝对观测误差方程为：

式中，gi为i点的平差重力值；为i点的观测重力值；pi为根据绝对重力观测值与相对重力观测值的实际精度确定的权。

经过固体潮改正、海潮改正、气压改正、仪器高改正、零漂改正后，对一台仪器在第i点和j点之间的段差观测值的误差方程为：

式中，gi、gj分别为测站i、j点平差后的重力值；gRZi、gRZj分别为测站i、j点经过四项改正的相对联测最后观测值；Ri、Rj分别为仪器在测站i、j点的观测读数；Ck为重力仪的K次格值改正因子；M为格值因子的最高次幂，一般M=1，只有个别仪器M=2；pij为相对观测量的权；Xn、Yn为周期误差的参数；Tn为周期误差的周期。

采用相对测量观测量与绝对测量观测量，依据误差方程式（1）、（2），形成误差方程：

式中，V为残差向量；A为系数矩阵；X为未知向量；L为观测向量；P为权矩阵。最小二乘解为：

其未知数的协因数阵：

单位权方差：

为了便于选权方法的阐述，以第i个观测值为例

式中，ai为矩阵A第i行元素组成的向量；Li为第i个观测值；vi为第i个观测值的改正数；pi为对应观测值的权。于是，当观测值之间相互独立时，定义准则函数为：

式中，ρ(vi)为vi的凸函数；令ρ′(vi)=ψ(vi)，可得

IGG III[16]的等价权定义如下：

式中，ko=1.0～1.5；k1=3.0～4.5。

可以得到如下未知参数估计及其协方差矩阵为：

2 重力网平差验后单位权方差的检验

接受H0假设，检验通过。反之，则拒绝H0，接受H1。未通过检验，认为观测值中含有粗差或验前权选取不合适。

在重力网的最小二乘平差完成后，进行χ2分布检验，若未通过检验，则需进行抗差估计。在计算完成后，再一次进行χ2检验。若不满足式（15），一方面需要对观测值进行检查，数据中可能存在粗差，另一方面需要对验前权矩阵进行检查，即重新选择合适的验前权阵。

3 数据与平差流程

3.1 实验数据

以某地区重力加密网为例，测区地理特征以平原、高原、山区为主，地区交通条件较好，气候的季风性特征显著。在数据采集过程中，采用了14台CG5/6仪器，由147条相对重力测量边段组成，获得了1 156条相对观测数据，平差的数据均通过闭合差检核。加密网的目标是实现对该地区的地面重力基准进行更新，网型结构如图1所示。

图1 重力加密网分布图

3.2 观测数据平差流程

略去对观测数据进行气压改正、仪器高改正、潮汐改正及零漂改正等处理过程的讨论。依据最小二乘法进行处理时，文献中多有叙述[18-19]。在本次实验项目中，采用迭代抗差估计计算，主要流程如下：

1）由初始解求得X^(0)，V(0)，并计算相应的权P(0)；

2）如果定义第k次迭代的残差为V(k)，等价权P(k)则通过式（12）计算得到；若pˉk=0，则认为该观测值含有粗差，不参与平差计算。

3）求解第k+1次迭代参数估值X^(k+1)和残差:

4 数据处理与分析

为了比较最小二乘法、Huber法、IGG III对含有粗差数据的处理效果，通过检验稳健估计结果对比分析，确定哪种方法抗差效果更好。在加密网中均匀挑选8条边，分别在每一边的终点重力观测值中加入0.1×10-5ms-2作为粗差进行平差处理。粗差数据在图1中加粗黑线位置处。

从表1中看出，最小二乘计算结果容易受到数据中粗差的影响。在选权迭代时，Huber方案的阈值为中误差的1.35倍作为计算值，对加入100微伽粗差数据处理时，点值平均中误差较最小二乘法有明显提高，权值更新时迭代了4次[20]，仅对其中6条粗差数据降权。IGGIII方案比Huber方案的处理精度也有一定的提高，权值迭代了5次，并将8条粗差数据的权值降为零。

表1 不同情况下各方案的平差结果对/（1×10-8 ms-2）

对未包含100微伽的数据处理时，检查IGGIII方案解的可靠性，对相对观测量的残差进行统计，从图2中可以看出，残差统计量符合正态分布的要求[10]。

图2 残差修正的预测曲线

图2 相对观测量数据平差后残差统计图

5 结语

根据多次试算结果分析，选用IGG III权因子函数，选择ko=1.5，k1=4.5，使平差计算的降权过程有效、平稳。抗差估计模型对粗差较为敏感，可准确的判断出粗差的位置和大小。

2000国家重力网的259个点位的重力点值平均中误差为7.4×10-8ms-2，而本次加密网139个点结果为6.5×10-8ms-2，文中所提方案有利于更大范围的重力数据处理工作，可为区域或国家的重力基准建立，提升测绘基准的现势性提供参考。