航空重力測量測線網平差中的粗差處理*

周波陽 羅志才 林 旭 周 浩

(1)武漢大學測繪學院，武漢 430079

2)武漢大學地球空間環(huán)境與大地測量教育部重點實驗室，武漢430079)

航空重力測量測線網平差中的粗差處理*

周波陽1)羅志才1，2)林 旭1)周 浩1)

(1)武漢大學測繪學院，武漢 430079

2)武漢大學地球空間環(huán)境與大地測量教育部重點實驗室，武漢430079)

航空重力測量測線網平差的主要任務是進行系統(tǒng)誤差的補償。為了削弱粗差觀測值的影響，提高補償效果，引入閾值法和抗差估計，并采用模擬數據驗證了這兩種方法的有效性和可靠性。通過比較發(fā)現(xiàn)，抗差估計在航空重力測線網平差中具有更好的實用價值。

航空重力測量;測線網平差;閾值法;抗差估計;權函數

1 引言

在航空重力測量中，為了盡可能地構成重復觀測，正確評價測量精度，測線往往布設成交叉網狀，在交叉點處兩條測線上的重力值理論上應該是相等的，但由于測量平臺在高空作業(yè)，處于不斷高速的運動狀態(tài)，觀測環(huán)境較為復雜，各類觀測值不可避免地受到“污染”，從而導致觀測值在交叉點處存在不符值。采用低通數字濾波等技術消除高頻噪聲后，觀測值中仍存在各種誤差源的干擾和影響。這些干擾和影響按其表現(xiàn)形式可分為偶然誤差、系統(tǒng)誤差以及粗差。偶然誤差的原因十分復雜，空氣擾動、機器零件的摩擦、氣壓、氣溫的變化等，以及它們的綜合影響均能產生;系統(tǒng)誤差主要由厄特弗斯效應、姿態(tài)誤差、衛(wèi)星星歷誤差等引起;GPS、加速度計等的工作狀態(tài)反生異常時觀測值中容易出現(xiàn)粗差。

航空重力測量測線網平差中的粗差處理通常和系統(tǒng)誤差的補償同步進行，忽略粗差直接采用最小二乘來估計系統(tǒng)誤差模型參數時，所建立的誤差模型并不能真實反映系統(tǒng)誤差的變化規(guī)律。平差后，大粗差的影響也會被分配到許多正常觀測值上，從而扭曲平差結果，參數估值的效率和可靠性就會受到嚴重的損害［1］。有效地消除或削弱觀測數據中粗差對于參數估值的不良影響，是提高航空重力測量網平差精度的關鍵環(huán)節(jié)。一般來說，對粗差的處理方式主要有3類:第一類即計算機輔助圖形可視化法，通過人機交互的方式來實現(xiàn)粗差的探測和剔除。當前數據處理都趨于自動化，人為干預排除粗差的方式將被淘汰［2］;第二類是在平差前采用某些方法對粗差定位并剔除，Dixon準則［3］、Grubbs準則［3］、χ2檢驗［4］、τ標準檢驗［4］、F檢驗［5］、t檢驗［6］、小波分析方法［3，7］以及閾值法［3］可以對粗差定位。其中Dixon準則無需估算樣本均值、樣本標準差和殘余誤差，在將所有觀測值按照從小到大的順序排列后，Dixon準則認為粗差觀測值僅存在于重排序之后的數據兩端。Dixon準則僅適合于小樣本(3＜n＜30)，采用Dixon準則進行小樣本數據的粗差探測能夠得到較嚴密的結果［8］。Grubbs檢驗、F檢驗等統(tǒng)計檢驗法都派生于最小二乘估計，對于偏離正態(tài)分布的數據是不可靠的，它們幾乎不可能把剔除粗差與估計量的計算分開，利用不可靠的估計量來進行統(tǒng)計檢驗就很難保證這種統(tǒng)計檢驗方法的正確性［2］。小波分析方法實質上是對觀測值信號進行濾波，利用小波變換的高低頻分離的特點，可在不丟失原始信號重要信息成分的前提下，對信號進行放大并對其進行特征提取，它也是一種較為有效的粗差探測方法，但算法較為復雜。閾值法計算方便簡單，適用于觀測數據的快速粗差探測;第三類則是基于多維M估計建立參數模型的抗差解式，采用等價權思想將其轉化為最小二乘解形式［9］。采用迭代的抗差M估計，既保留了最小二乘法的優(yōu)越性，又保證被估參數既能抵制模型誤差又能抗拒粗差的影響［1，10］。本文將抗差M估計和閾值法應用于航空重力測量測線網平差中的粗差處理，并通過模擬計算對這兩種方法進行評價。

2 系統(tǒng)誤差補償的兩步法

航空重力測量的測線網平差實際上是系統(tǒng)誤差的補償問題，進行系統(tǒng)誤差補償的一般方法為自檢校測線網平差法，但自檢校平差會導致法方程奇異，為了解決這個問題，將自檢校測線網平差過程進行簡化，將誤差理論模型簡化為實用模型，即為誤差補償的兩步法［11-13］。

2.1 交叉點平差

第一階段使用條件平差法對測線交叉點進行平差。在主測線、副測線的交叉點處可建立如下的誤差方程:

式中，Δgik與Δgjk分別為主測線i、副測線j在交叉點k處的重力異常，Δgik-Δgjk表示交叉點不符值，Vik、Vjk為改正數。對于具有多個交叉點的測線網，可以寫出交叉點條件方程的矩陣形式為

其中，V由每一條測線上交叉點處的改正數組成，系數陣B是由1、-1和大量的0構成，W為交叉點不符值。依據最小二乘平差理論可解得

一般來說，測線上各個測點為獨立等精度觀測，權陣P為單位陣。

2.2 測線系統(tǒng)誤差估計

以觀測時間t為自變量，以式得到的交叉點不符值的改正數V為虛擬觀測值，遍歷所有測線對每個交叉點建立如下形式的誤差方程

vi為V中的第i個元素，f(t)為系統(tǒng)誤差模型，Δi為噪聲，m為該測線上的交叉點個數。f(t)有如下兩種形式:

1)一般多項式模型

2)三角多項式形式

則誤差方程的矩陣形式為

3 粗差處理方法

3.1 閾值法

Huber［2］在1977年采用Tukey污染分布模擬了一些含有不同數量粗差的正態(tài)分布數據，其污染正態(tài)分布模式為

其中Φ是標準正態(tài)分布，ε為污染率，是一小量，μ和σ分別為母體均值和標準差。式(9)表明，與均值之差的絕對值大于3σ的觀測值被認為是粗差。而實際觀測中，被測量的真實值大多是未知，μ和σ不能精確知道，通常以樣本均值和標準差S代替。閾值法的實用形式為

在實際操作過程中，把交叉點平差后的不符值的改正數V作為觀測值向量X，當觀測值xi(xi實際上為交叉點不符值的改正數vi)滿足式(10)時，被認作是粗差。

3.2 抗差M估計

經典最小二乘估計是利用一組來自于母體為正態(tài)分布的觀測值，求定母體參數的一種參數估計，它的最優(yōu)統(tǒng)計特性(不能脫離觀測值母體)是正態(tài)分布，即觀測值中僅含有偶然誤差這一前提［14］。因此，最小二乘估計對偏離正態(tài)分布的數據是不可靠的。與最小二乘估計相比，抗差估計應具備兩大特點:一是它能消除和削弱粗差對估值的影響;二是它基本上具備經典估計的優(yōu)良特性。

按M估計原理，取極值函數為

其中，ai為系數陣A的第i行。對X求導，并令其為0，同時記φ(Vi)=?ρ/?Vi，則有

令φ(Vi)/Vi=Wi(權因子)=piiWi為等價權元素，則式(12)可寫為

由此，可得到抗差M估值

抗差估計采用迭代解法，每一次迭代相當于進行一次最小二乘估計，但應以等價權代替先驗權。式(7)的第k+1步迭代解為

抗差估計的結果取決于所選取的等價權函數，權函數不一樣，抗差估計的模型就不同，抗差的效果也不一樣［15］。常見的等價權模型有Huber權函數、Hampel權函數、Tukey權函數和 Andrews權函數等［9］。本文僅以Huber等價權模型為例進行了抗差估計，其具體形式為:

式中c一般取2.0。

4 算例及結果分析

本文模擬了某航空重力測區(qū)，它由10條橫向測線和10條縱向測線組成，飛機起伏飛行，飛行時離地面高度約為1.8 km。應用兩步處理法對觀測數據進行處理，第一階段交叉點平差后主、副測線驗后單位權中誤差見表1～2。為了更好地檢驗閾值法和抗差M估計的實際效果，在這些交叉點不符值的改正數中隨機挑選10個加入絕對值分別為(10～18)×10-5ms-2和(6～10)×10-5ms-2的粗差，形成兩組觀測值data1和data2，其交叉點的分布如圖1。在兩步處理法第二階段中分別采用閾值法和抗差估計來消除粗差的影響，系統(tǒng)誤差模型為2階和3階一般多項式，表3～9給出了相應的最小二乘估計(LSE)和抗差估計(RE)結果。

圖1 交叉點的分布Fig.1 Distribution of crosspoints

由表1可見，用閾值法進行粗差的剔除后，數據的標準偏差變小，測線網精度有所提高。這說明閾值法能有效對粗差進行識別和定位，是一種可靠的粗差探測手段，但它的缺點也很明顯，其對絕對值較小的粗差不敏感，且在剔除粗差觀測值之后，會導致交叉點個數變少，這將給后續(xù)研究帶來不利影響。

從表4～7可以看出，當部分觀測值含有粗差時，主、副測線的最小二乘估計驗后單位權中誤差均大于抗差估計驗后單位權中誤差，這說明當觀測值含有粗差時，從整體精度上，采用抗差估計進行系統(tǒng)誤差補償的成果明顯優(yōu)于采用最小二乘估計進行系統(tǒng)誤差補償的成果。但我們也不能單純的以交叉點不符值的均方差大小作為衡量抗差效果唯一標準。當信噪比一定時，計算結果很大程度上取決于誤差模型的合理選擇，從本算例來看，二階多項式的抗差效果更好。反之，當系統(tǒng)誤差模型一定時，信噪比越低，觀測值含粗差的絕對值較大時，抗差效果越明顯。

表1 不含粗差時主測線系統(tǒng)誤差模型驗后單位權中誤差(單位：10-5ms-2)Tab.1 Posteriori unit weight mean errors of system error model on main line of data without outliers(unit：10-5ms-2)

表2 不含粗差時副測線系統(tǒng)誤差模型驗后單位權中誤差(單位：10-5ms-2)Tab.2 Posteriori unit weight mean errors of system error model on vice line of data without outliers(unit：10-5ms-2)

表3 閾值法處理結果(單位：10-5ms-2)Tab.3 Results of Thresholding method(unit：10-5ms-2)

表4 data1主測線系統(tǒng)誤差模型參數驗后單位權中誤差(單位：10-5ms-2)Tab.4 Posteriori unit weight mean errors of system error model on main line of data1(unit：10-5ms-2)

表5 data1副測線系統(tǒng)誤差模型參數驗后單位權中誤差(單位：10-5ms-2)Tab.5 Posteriori unit weight mean errors of system error model on vice line of data1(unit：10-5ms-2)

表6 data2主測線系統(tǒng)誤差模型參數驗后單位權中誤差(單位：10-5ms-2)Tab.6 Posteriori unit weight mean errors of system error model on main line of data2(unit：10-5ms-2)

表7 data2副測線系統(tǒng)誤差模型參數驗后單位權中誤差(單位：10-5ms-2)Tab.7 Posteriori unit weight mean errors of system error model on deputy line of data2(unit：10-5ms-2)

表8 data1交叉點最小二乘殘差與抗差估計殘差(單位：10-5ms-2)Tab.8 Residuals of least squares estimation and robust estimation of data1(unit：10-5ms-2)

表9 data2交叉點最小二乘殘差與抗差估計殘差(單位：10-5ms-2)Tab.9 Residuals of least squares estimation and robust estimation of data2(unit：10-5ms-2)

表8～9列出了不含粗差時最小二乘殘差、加粗差后最小二乘殘差以及加粗差后抗差M估計的殘差(系統(tǒng)誤差模型為二階多項式)。不難發(fā)現(xiàn)，加粗差后抗差M估計的殘差比最小二乘估計的殘差更接近于粗差值本身，在異常觀測值上抗差M估計具有明顯的抵御粗差干擾的能力，這是因為在迭代計算時，抗差M估計利用等價權模型對異常觀測值實施了降權操作。通過對殘差進行分析對比，抗差估計不僅能對粗差的位置進行識別，而且還能對粗差的大小進行估計。上述結果還表明，當粗差絕對值較小時，抗差估計仍然能夠對其進行識別和定位。從這方面來講，抗差估計比閾值法更加有效和可靠。

5 結語

對于各測線上的所有采樣點的系統(tǒng)誤差補償來說，閾值法和抗差估計均具有一定的抗差效果。閾值法是直接剔除異常觀測值，而抗差M估計能在拒絕和接受一個觀測值之間起一個平滑的作用。相比較而言，抗差估計可靠性更高，更加有效。值得注意的是，任何一種抗差估計方法的實際抗差能力都是有一定條件的，當粗差個數過多、絕對值較小時，閾值法和抗差估計對粗差的處理結果往往不夠理想。后續(xù)研究將考慮更加復雜的系統(tǒng)誤差模型(如三角多項式)和更為豐富的權函數模型以及對系統(tǒng)誤差模型的參數進行顯著性檢驗。

1 楊元喜.抗差估計理論及其應用［M］.北京:八一出版社，1993.(Yang Yuanxi.Theory of robust estimation and its application［M］.Beijing:Bayi Press，1993)

2 黃幼才.數據探測與抗差估計［M］.北京:測繪出版社，1990.(Huang Youcai.Data detection and robust estimation［M］.Beijing:Mapping Press，1993)

3 Kern M，et al.Outlier detection algorithms and their performance in GOCE gravity field Processing［J］.Journal of Geodesy，2005，78(9):509-519.

4 Hwang C，Wang C G and Lee L H.Adjustment of relative gravity measurements using weighted and datum-free constraints［J］.Computeramp;Geoscience，2002，28:1 005-1 015.

5 李德仁.誤差處理和可靠性理論［M］.北京:測繪出版社，1988.(Li Deren.Error processing and theory of reliablity［M］.Beijing:Mapping Press，1988)

6 丁士俊，陶本藻.半參數模型的統(tǒng)計診斷量與粗差檢驗［J］.大地測量與地球動力學，2005，(3):24-28.(Ding Shijun and Tao Benzao.Measures of statistical diagnosis and gross test［J］.Journal of Geodesy and Geodynamics，2005，(3):24-28)

7 Wang Y.Jump and sharp cusp detection by wavelets［J］.Biometrika，1995，82(2):385-397.

8 楊茂興.小樣本容量測量數據中粗差的剔除［J］.計量與測試技術，2005，(1):27-28.(Yang Maoxing.The coarse error eliminating in small sample space［J］.Metrology and Measurement Technique，2005(1):27-28)

9 楊元喜.參數平差模型的抗差最小二乘解［J］.測繪通報，1994，(6):33-35.(Yang Yuanxi.Robust estimation solutions of parameter adjustment model［J］.Bulletin of Surveying and Mapping，1994，(6):33-35)

10 周江文.經典誤差理論與抗差估計［J］.測繪學報，1989，18(2):115-120.(Zhou Jiangwen.Classical theory of error and robust estimation［J］.Acta Geodaetica et Carto graphica Sinica，1989，18(2):115-120)

11 黃謨濤，等.海洋重力測量系統(tǒng)誤差補償的兩步處理法［J］.武漢大學學報(信息科學版)，2002，27(3):251 -255.(Huang Motao，et al.Two-step processing for compensating the systemic errors in marine gravity measurements［J］.Geomatics and Information Science of Wuhan University，2002，27(3):251-255)

12 李海.航空重力測量測線網平差的理論與方法［D］.解放軍信息工程大學，2002.(Li Hai.Theory and methods of network adjustment of airborne gravimetry［D］.The PLA Information Engineering University，2002)

13 蔡劭琨.航空重力測量網絡平差方法研究［D］.國防科學技術大學，2009.(Cai Shaokun.Research on the network adjustment of airborne gravimetry［D］.National University of Defense Technology，2009)

14 余學詳，呂偉才.抗差估計在粗差探測和平差計算中的運用［J］.測繪工程，1998，8(3):40-44.(Yu Xuexiang and Lv Weicai.The application of robust estimation to the gross error detection and adjustment caculation［J］.Engineering of Surveying and Mapping，1998，8(3):40-44)

15 余學祥，呂偉才.基于標準化的相關觀測抗差估計模型［J］.武漢測繪科技大學學報，1999，24(1):72-78.(Yu Xuexiang and Lü Weicai.Robust estimation model for correlated observations based on standardized residuals［J］.Journal of Wuhan Technical University of Surveying and Mapping，1999，24(1):72-78)

PROCESSING METHODS FOR OUTLIERS IN NETWORK ADJUSTMENT OF AIRBORNE GRAVIMETRY

Zhou Boyang1)，Luo Zhicai1，2)，Lin Xu1)and Zhou Hao1)

(1)School of Geodesy and Geomatics，Wuhan University，Wuhan 430079
2)Key Laboratory of Geospace Environment and Geodesy，Wuhan University，Wuhan 430079)

The main task of network adjustment of airborne gravimetry is to compensate systemic errors.In order to weaken the influence of outliers and improve the effect of compensation，this paper brings in thresholding and robust estimation，and confirms these two methods’validity and reliability in simulated data processing.By means of comparison，a conclusion is reached that the robust estimation has more practical value in the network adjustment of airborne gravimetry.

airborne gravimetry;network adjustment;thresholding;robust estimation;weight function

1671-5942(2012)02-0110-05

2011-12-25

國家自然科學基金(41174062);中國博士后科學基金(20110491189);中央高校基本科研業(yè)務費專項資金(111110)

周波陽，1983年生，博士研究生，目前從事航空重力測量數據處理研究.E-mail:byzhou@whu.edu.cn

P207

A