IGG III選權(quán)迭代法的重力網(wǎng)抗差估計

李 琦，武 威

（1.西安測繪總站，陜西 西安 710054）

重力加密測量為重力基準(zhǔn)的建立、油氣資源的勘探、地質(zhì)構(gòu)造及地下水資源的開發(fā)利用提供了重要的基礎(chǔ)數(shù)據(jù)和手段[1-3]。我國2000國家重力基本網(wǎng)已經(jīng)使用了近20 a，隨著經(jīng)濟(jì)社會的發(fā)展部分點位遭到破壞[4]，以及地球內(nèi)部及表面的質(zhì)量調(diào)整，進(jìn)行加密重力測量彌補(bǔ)2000國家重力基本網(wǎng)在部分地區(qū)的不足就尤為必要。相對重力儀的非線性零點漂移及格值系數(shù)變化、重力場時變信號改正不盡完善，另外野外作業(yè)人員操作失誤等因素造成觀測數(shù)據(jù)中可能包含有較大誤差、甚至是粗差[5-6]。同時，加密重力測量也是一項費(fèi)時、費(fèi)力的精細(xì)工作，獲取的數(shù)據(jù)資料珍貴[7]，故對加密重力測量數(shù)據(jù)進(jìn)行高精度處理就顯得非常關(guān)鍵。

在2000國家重力基本網(wǎng)平差中，采用“弱基準(zhǔn)”的方式，對絕對觀測量及相對觀測量將賦以適當(dāng)?shù)臋?quán)進(jìn)行平差，以減弱兩者精度不匹配的影響[8]。在沿海省市的重力網(wǎng)數(shù)據(jù)處理中，用抗差等價權(quán)來調(diào)整相對重力的權(quán)、重力儀參數(shù)的取舍[9]。

抗差估計的核心是按照權(quán)值將觀測值進(jìn)行正常觀測值、可利用觀測值和粗差觀測值，權(quán)值大小不變的為正常觀測值。對殘差較大的觀測值進(jìn)行降權(quán)處理，為可利用觀測值，而權(quán)值降為零的觀測值為粗差觀測值，從而避免粗差對結(jié)果的影響[10]。最小截斷二乘法（least trimmed squares，LTS）的抗差方案中，處理含有粗差的大樣本量數(shù)據(jù)處理效率不高[11]。最小二乘估計不具有抗干擾性[12-14]，抗差最小二乘估計通過等價權(quán)把抗差估計與最小二乘法結(jié)合起來。等價權(quán)的選取有L1法、German-McClure法、Danish法、Huber法IGG法和IGG III等多種方法[15-16]。為此，以某區(qū)域加密重力網(wǎng)為例，在平差中采用IGG III的權(quán)因子函數(shù)[16]，進(jìn)行權(quán)值迭代平差，獲取可靠的參數(shù)估值。

1 重力網(wǎng)的抗差估計方案

依據(jù)絕對重力觀測值，其絕對觀測誤差方程為：

式中，gi為i點的平差重力值；為i點的觀測重力值；pi為根據(jù)絕對重力觀測值與相對重力觀測值的實際精度確定的權(quán)。

經(jīng)過固體潮改正、海潮改正、氣壓改正、儀器高改正、零漂改正后，對一臺儀器在第i點和j點之間的段差觀測值的誤差方程為：

式中，gi、gj分別為測站i、j點平差后的重力值；gRZi、gRZj分別為測站i、j點經(jīng)過四項改正的相對聯(lián)測最后觀測值；Ri、Rj分別為儀器在測站i、j點的觀測讀數(shù)；Ck為重力儀的K次格值改正因子；M為格值因子的最高次冪，一般M=1，只有個別儀器M=2；pij為相對觀測量的權(quán)；Xn、Yn為周期誤差的參數(shù)；Tn為周期誤差的周期。

采用相對測量觀測量與絕對測量觀測量，依據(jù)誤差方程式（1）、（2），形成誤差方程：

式中，V為殘差向量；A為系數(shù)矩陣；X為未知向量；L為觀測向量；P為權(quán)矩陣。最小二乘解為：

其未知數(shù)的協(xié)因數(shù)陣：

單位權(quán)方差：

為了便于選權(quán)方法的闡述，以第i個觀測值為例

式中，ai為矩陣A第i行元素組成的向量；Li為第i個觀測值；vi為第i個觀測值的改正數(shù)；pi為對應(yīng)觀測值的權(quán)。于是，當(dāng)觀測值之間相互獨立時，定義準(zhǔn)則函數(shù)為：

式中，ρ(vi)為vi的凸函數(shù)；令ρ′(vi)=ψ(vi)，可得

IGG III[16]的等價權(quán)定義如下：

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

可以得到如下未知參數(shù)估計及其協(xié)方差矩陣為：

2 重力網(wǎng)平差驗后單位權(quán)方差的檢驗

接受H0假設(shè)，檢驗通過。反之，則拒絕H0，接受H1。未通過檢驗，認(rèn)為觀測值中含有粗差或驗前權(quán)選取不合適。

在重力網(wǎng)的最小二乘平差完成后，進(jìn)行χ2分布檢驗，若未通過檢驗，則需進(jìn)行抗差估計。在計算完成后，再一次進(jìn)行χ2檢驗。若不滿足式（15），一方面需要對觀測值進(jìn)行檢查，數(shù)據(jù)中可能存在粗差，另一方面需要對驗前權(quán)矩陣進(jìn)行檢查，即重新選擇合適的驗前權(quán)陣。

3 數(shù)據(jù)與平差流程

3.1 實驗數(shù)據(jù)

以某地區(qū)重力加密網(wǎng)為例，測區(qū)地理特征以平原、高原、山區(qū)為主，地區(qū)交通條件較好，氣候的季風(fēng)性特征顯著。在數(shù)據(jù)采集過程中，采用了14臺CG5/6儀器，由147條相對重力測量邊段組成，獲得了1 156條相對觀測數(shù)據(jù)，平差的數(shù)據(jù)均通過閉合差檢核。加密網(wǎng)的目標(biāo)是實現(xiàn)對該地區(qū)的地面重力基準(zhǔn)進(jìn)行更新，網(wǎng)型結(jié)構(gòu)如圖1所示。

圖1 重力加密網(wǎng)分布圖

3.2 觀測數(shù)據(jù)平差流程

略去對觀測數(shù)據(jù)進(jìn)行氣壓改正、儀器高改正、潮汐改正及零漂改正等處理過程的討論。依據(jù)最小二乘法進(jìn)行處理時，文獻(xiàn)中多有敘述[18-19]。在本次實驗項目中，采用迭代抗差估計計算，主要流程如下：

1）由初始解求得X^(0)，V(0)，并計算相應(yīng)的權(quán)P(0)；

2）如果定義第k次迭代的殘差為V(k)，等價權(quán)P(k)則通過式（12）計算得到；若pˉk=0，則認(rèn)為該觀測值含有粗差，不參與平差計算。

3）求解第k+1次迭代參數(shù)估值X^(k+1)和殘差:

4 數(shù)據(jù)處理與分析

為了比較最小二乘法、Huber法、IGG III對含有粗差數(shù)據(jù)的處理效果，通過檢驗穩(wěn)健估計結(jié)果對比分析，確定哪種方法抗差效果更好。在加密網(wǎng)中均勻挑選8條邊，分別在每一邊的終點重力觀測值中加入0.1×10-5ms-2作為粗差進(jìn)行平差處理。粗差數(shù)據(jù)在圖1中加粗黑線位置處。

從表1中看出，最小二乘計算結(jié)果容易受到數(shù)據(jù)中粗差的影響。在選權(quán)迭代時，Huber方案的閾值為中誤差的1.35倍作為計算值，對加入100微伽粗差數(shù)據(jù)處理時，點值平均中誤差較最小二乘法有明顯提高，權(quán)值更新時迭代了4次[20]，僅對其中6條粗差數(shù)據(jù)降權(quán)。IGGIII方案比Huber方案的處理精度也有一定的提高，權(quán)值迭代了5次，并將8條粗差數(shù)據(jù)的權(quán)值降為零。

表1 不同情況下各方案的平差結(jié)果對/（1×10-8 ms-2）

對未包含100微伽的數(shù)據(jù)處理時，檢查IGGIII方案解的可靠性，對相對觀測量的殘差進(jìn)行統(tǒng)計，從圖2中可以看出，殘差統(tǒng)計量符合正態(tài)分布的要求[10]。

圖2 殘差修正的預(yù)測曲線

圖2 相對觀測量數(shù)據(jù)平差后殘差統(tǒng)計圖

5 結(jié)語

根據(jù)多次試算結(jié)果分析，選用IGG III權(quán)因子函數(shù)，選擇ko=1.5，k1=4.5，使平差計算的降權(quán)過程有效、平穩(wěn)?？共罟烙嬆Ｐ蛯Υ植钶^為敏感，可準(zhǔn)確的判斷出粗差的位置和大小。

2000國家重力網(wǎng)的259個點位的重力點值平均中誤差為7.4×10-8ms-2，而本次加密網(wǎng)139個點結(jié)果為6.5×10-8ms-2，文中所提方案有利于更大范圍的重力數(shù)據(jù)處理工作，可為區(qū)域或國家的重力基準(zhǔn)建立，提升測繪基準(zhǔn)的現(xiàn)勢性提供參考。