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，文中所提方案有利于更大范圍的重力數據處理工作，可為區域或國家的重力基準建立，提升測繪基準的現勢性提供參考。