一種基于丹麥法的改進(jìn)型雙步M 估計(jì)

林國(guó)標(biāo) 劉立龍 蔡成輝 黎峻宇 黃良珂

1 廣西礦冶與環(huán)境科學(xué)實(shí)驗(yàn)中心，桂林市建干路12號(hào)，541004

2 桂林理工大學(xué)測(cè)繪地理信息學(xué)院，桂林市建干路12號(hào)，541004

3 廣西空間信息與測(cè)繪重點(diǎn)實(shí)驗(yàn)室，桂林市建干路12號(hào)，541004

抗差估計(jì)的抗差性及有效性主要取決于初值的準(zhǔn)確性與權(quán)函數(shù)的合理性［1］。目前，一般選取最小二乘估計(jì)的結(jié)果作為抗差估計(jì)的初值［2］，但最小二乘估計(jì)對(duì)粗差具有均衡性和不敏感性，致使粗差觀測(cè)值的殘差并非最大，導(dǎo)致選權(quán)迭代時(shí)出現(xiàn)錯(cuò)誤的判斷。同時(shí)，合理的權(quán)函數(shù)應(yīng)該分3段：正常段、可疑段、淘汰段，其臨界值和均方差因子應(yīng)根據(jù)實(shí)際情況選取，具有可變性［1，3］。在傳統(tǒng)的抗差估計(jì)方法中，Andrews法、Tukey 法不含正常段，Huber 法、丹麥法沒有淘汰段，IGG法、Hampel法雖然具有正常段、可疑段及淘汰段，但其臨界值往往取某一固定的先驗(yàn)值。這些方法僅采用一步抗差估計(jì)［4］，具有一定的局限性。

針對(duì)以上問題，本文提出一種基于丹麥法的改進(jìn)型雙步M 估計(jì)方法。采用抗差性較強(qiáng)的L1－范估計(jì)的平差值作為抗差估計(jì)的初值，在丹麥法等價(jià)權(quán)函數(shù)的基礎(chǔ)上引入淘汰域，將改進(jìn)的丹麥法權(quán)函數(shù)作為抗差權(quán)，選取可變動(dòng)的臨界值，通過雙步M 估計(jì)，提高抗差估計(jì)的精度。

1 改進(jìn)型雙步M 估計(jì)

1.1 初值選取

由于LS估計(jì)不具抗差性，LS估計(jì)的殘差不能有效反映粗差的位置。而L1－范估計(jì)的平差值具有較好的抗差性，能較為準(zhǔn)確地體現(xiàn)真實(shí)殘差［5］，本文將其作為雙步M 估計(jì)的初值。

首先，列誤差方程：

解得殘差：

然后，以L1－范權(quán)函數(shù)重新定權(quán)，再進(jìn)行迭代計(jì)算，收斂時(shí)得殘差初值：

式中，Pi為觀測(cè)權(quán)，為等價(jià)權(quán)，V0為初始?xì)埐?，l為常數(shù)項(xiàng)，c為一個(gè)很小的正常數(shù)。

1.2 等價(jià)權(quán)函數(shù)選擇

丹麥法的權(quán)函數(shù)如下［3］：

丹麥法權(quán)函數(shù)缺乏淘汰段，導(dǎo)致抗差效率降低。加入淘汰段，對(duì)丹麥法進(jìn)行改進(jìn)：

1.3 第一步——綜合抗差階段

該階段有兩個(gè)關(guān)鍵問題需要解決：一是確定臨界值k0、k1，二是選取等價(jià)權(quán)。k0、k1一般分別取某一固定值［3］，k0＝1～2.5，k1＝3～6。臨界值在k0、k1基礎(chǔ)上再乘以一個(gè)可變因子（r為多余度），可解決臨界值取值不當(dāng)?shù)膯栴}，實(shí)現(xiàn)系數(shù)陣空間與觀測(cè)空間的同時(shí)抗差［6］。r的平均值為（n－m）／n（n、m分別為觀測(cè)值個(gè)數(shù)與參數(shù)個(gè)數(shù)），故可變因子為。此時(shí)臨界值會(huì)隨具體問題而變化，增強(qiáng)了抗差估計(jì)的靈活性。本文取k0＝2d，k1＝4d。

選取改進(jìn)后的丹麥法等價(jià)權(quán)作為抗差權(quán)，單位權(quán)中誤差σ0采用第一次或第二次抗差估計(jì)所求得的觀測(cè)殘差來計(jì)算，能夠獲得比較準(zhǔn)確的均方差因子［7］：

1.4 第二步——優(yōu)化階段

將第一步所得的估計(jì)值作為該階段的初值，利用截尾LS法重新平差計(jì)算。取權(quán)函數(shù)［8］：

臨界值取k＝3。根據(jù)此時(shí)的v與σ2，將大于3倍均方差的觀測(cè)值刪去，然后利用剩余的觀測(cè)值進(jìn)行加權(quán)LS 平差，得到最終的參數(shù)解。精度評(píng)定采用：

式中，n、m、t分別為總觀測(cè)數(shù)、必要觀測(cè)數(shù)和刪去的粗差觀測(cè)數(shù)；為最終權(quán)。

2 算例分析

本文結(jié)合一個(gè)測(cè)邊網(wǎng)算例和一個(gè)水準(zhǔn)網(wǎng)算例［9］，對(duì)比本文提出的改進(jìn)型雙步M 估計(jì)與丹麥法、Huber法、L1法、IGGⅠ法及兩步抗差估計(jì)5種方法的平差結(jié)果。將無粗差時(shí)的LS法平差值當(dāng)作真值，傳統(tǒng)抗差估計(jì)方法和兩步抗差估計(jì)的臨界值取先驗(yàn)值，均方差因子通過初次平差的殘差由中位數(shù)法計(jì)算，即0.674 5。

算例1 如圖1，測(cè)邊網(wǎng)中A、B、C、D為已知點(diǎn)，P1、P2、P3、P4為未知點(diǎn)，觀測(cè)13條邊長(zhǎng)，測(cè)距精度σS＝3mm＋1×10－6S。起算數(shù)據(jù)及觀測(cè)邊長(zhǎng)見表1。在第3、12條邊上分別附上－30dm、50 dm 的粗差，各種抗差估計(jì)方法的計(jì)算結(jié)果見表2。

圖1 測(cè)邊網(wǎng)圖Fig.1 Trilateration network

算例中，測(cè)邊網(wǎng)的污染率為15.4%，屬于較嚴(yán)重污染。由表2可知：

1）在參數(shù)估計(jì)方面，幾種傳統(tǒng)抗差估計(jì)方法的估計(jì)結(jié)果與真值之差較大：丹麥法、Huber法、L1法和IGGⅠ法最大的差值分別為11.2cm、30.4cm、20.3cm、43.1cm；兩步抗差估計(jì)的結(jié)果與真值之差較小，最大為－3.2cm；本文方法與真值之差較小，最大差值為－2.5cm，整體比兩步抗差估計(jì)小一些。

2）在最終權(quán)方面，P3、P12為粗差權(quán)，其余為正常觀測(cè)值的最終權(quán)，傳統(tǒng)抗差估計(jì)方法的P3、P12未能同時(shí)降到最小或0，如丹麥法、Huber法、L1法和IGGⅠ法的P3分別為0.03、0.11、0.03、0.13，L1法正常觀測(cè)值的最終權(quán)受粗差影響波動(dòng)很大；兩步抗差估計(jì)的P3、P12均降為0，而P6、P8、P13分別降為0.42、0.40、0.73；本文方法的P3、P12均降為0，其余權(quán)值均為觀測(cè)權(quán)。

若在上述粗差基礎(chǔ)上再加上一個(gè)粗差，即在第3、10、12條邊上分別附上－30dm、－20dm、50dm 的3個(gè)粗差，再一次進(jìn)行平差計(jì)算，則各種抗差估計(jì)方法的計(jì)算結(jié)果見表3。此時(shí)測(cè)邊網(wǎng)的污染率為23.1%。由表3可知：

表1 測(cè)邊網(wǎng)觀測(cè)數(shù)據(jù)及起算數(shù)據(jù)Tab.1 Observed value and known data

表2 各種抗差估計(jì)方法的計(jì)算結(jié)果Tab.2 Results by various robust estimation methods

表3 各種抗差估計(jì)方法的計(jì)算結(jié)果Tab.3 Results by various robust estimation methods

1）在參數(shù)估計(jì)方面，傳統(tǒng)抗差估計(jì)方法和兩步抗差估計(jì)的結(jié)果與真值之差比較大，如丹麥法、Huber法、L1法、IGGⅠ法和兩步抗差估計(jì)的最大差值分別為109.0、175.6、115.1、121.8、121.8 cm；本文方法的結(jié)果與真值之差較小，最大僅為－2.7cm。

2）在最終權(quán)方面，P3、P10、P12為粗差權(quán)，其余為正常觀測(cè)值的最終權(quán)。傳統(tǒng)抗差估計(jì)方法和兩步抗差估計(jì)的P3、P10仍為觀測(cè)權(quán)，分別為1.86、14.2，L1法的粗差權(quán)雖然最后有所降低，但是正常觀測(cè)權(quán)之比波動(dòng)很大；本文方法的P3、P10、P12均降為0，其余權(quán)仍為觀測(cè)權(quán)。

綜上所述，當(dāng)測(cè)邊網(wǎng)含有2個(gè)粗差時(shí)，傳統(tǒng)抗差估計(jì)方法的抗差效果較差，兩步抗差估計(jì)和本文方法的抗差效果較好，而本文方法的整體抗差效果優(yōu)于兩步抗差估計(jì)；當(dāng)測(cè)邊網(wǎng)含有3個(gè)粗差時(shí)，傳統(tǒng)抗差估計(jì)方法和兩步抗差估計(jì)基本喪失了抗差估計(jì)的能力，而本文方法的抗差效果依然比較顯著。這是由于傳統(tǒng)抗差估計(jì)方法和兩步抗差估計(jì)往往采用對(duì)粗差具有均衡性的LS法的平差結(jié)果作為抗差初值，不能準(zhǔn)確反映粗差的位置。而前者臨界值是固定的，不能有效適應(yīng)實(shí)際問題，因此它們抵御多維粗差的效果并不顯著。

算例2 如圖2，水準(zhǔn)網(wǎng)中A、B為已知點(diǎn)，HA＝8.016m，HB＝9.016m，點(diǎn)P1、P2、P3為待定點(diǎn)。進(jìn)行7條水準(zhǔn)路線觀測(cè)，見表4。在第1、4路線上分別附加－20mm、40mm 的粗差，各種方法的計(jì)算結(jié)果見表5。

圖2 水準(zhǔn)網(wǎng)圖Fig.2 Leveling network

表4 水準(zhǔn)網(wǎng)觀測(cè)高差值及路線長(zhǎng)Tab.4 Observed value and length

表5 各種抗差估計(jì)方法的計(jì)算結(jié)果Tab.5 Results by various robust estimation methods

算例中，多余觀測(cè)數(shù)為4，污染率為28.6%，屬于比較嚴(yán)重的污染，因此本算例只進(jìn)行包含2個(gè)粗差的比較。由表5可知：

1）在參數(shù)估計(jì)方面，傳統(tǒng)抗差估計(jì)方法和兩步抗差估計(jì)的估計(jì)結(jié)果與真值之差較大，丹麥法、Huber法、L1法、IGGⅠ法和兩步抗差估計(jì)的最大差值分別為－8.7 mm、－8.0 mm、－6.7 mm、－8.3mm、－8.8mm；本文方法的結(jié)果與真值之差較小，最大僅為－0.3mm。

2）在最終權(quán)方面，P1、P4為粗差權(quán)，其余為正常觀測(cè)值的最終權(quán)，傳統(tǒng)抗差估計(jì)和兩步抗差估計(jì)的P1、P4未能同時(shí)降為0，如丹麥法、Huber法、IGGⅠ法及兩步抗差估計(jì)的P1觀測(cè)權(quán)均為2.73，L1法的P1為0.21，L1法的P3、P7有所降低；本文方法的P1、P4均降為0，其余最終權(quán)維持不變。

和算例1同理，因?yàn)閭鹘y(tǒng)抗差估計(jì)和兩步抗差估計(jì)常常采用不能準(zhǔn)確反映粗差位置的LS法平差結(jié)果作為抗差初值，前者臨界值是固定的，所以往往不能抵抗多維粗差的干擾。而本方法以L1－范估計(jì)結(jié)果作為初值，第一步的抗差權(quán)采用改進(jìn)的丹麥法等價(jià)權(quán)，臨界值可變；第二步采用截尾LS法，取得了較好的估計(jì)結(jié)果。

3 結(jié) 語(yǔ)

基于丹麥法的改進(jìn)型雙步M 估計(jì)是一種以L1－范估計(jì)的平差值作為初值、以改進(jìn)后的丹麥等價(jià)權(quán)作為綜合抗差階段權(quán)因子的抗差估計(jì)方法，比較有效地解決了抗差估計(jì)初值與臨界值選取不當(dāng)?shù)膯栴}，提高了抗差估計(jì)的精度。算例表明，無論是測(cè)邊網(wǎng)平差，還是水準(zhǔn)網(wǎng)平差，該方法都能有效抵御多維粗差的干擾，獲得優(yōu)于傳統(tǒng)抗差估計(jì)和兩步抗差估計(jì)的平差結(jié)果，具有較強(qiáng)的抗差性。

［1］楊元喜，吳富梅.臨界值可變的抗差估計(jì)等價(jià)權(quán)函數(shù)［J］.測(cè)繪科學(xué)技術(shù)學(xué)報(bào)，2006（5）：320－324（Yang Yuanxi，Wu Fumei.Modified Equivalent Weight Function with Variable Criterion for Robust Estimation［J］.Journal of Geomatics Science and Technology，2006，（5）：320－324）

［2］邱衛(wèi)寧.具有穩(wěn)健初值的選權(quán)迭代法［J］.武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2003（4）：452－454（Qiu Weining.Method for Selecting Weight Iteration with Robust Initial Value［J］.Geomatics and Information Science of Wuhan University，2003，（4）：452－454）

［3］李浩軍，唐詩(shī)華，黃杰.經(jīng)典選權(quán)迭代法研究與兩步抗差估計(jì)的提出［J］.海洋測(cè)繪，2007（1）：17－20（Li Haojun，Tang Shihua，Huang Jie.The Discussion about Some Kinds of Selecting Weight Iteration Method［J］.Hydrographic Surveying and Charting，2007（1）：17－20）

［4］靳奉祥.抗差估計(jì)理論與方法研究［J］.山東科技大學(xué)學(xué)報(bào)：自然科學(xué)版，2003（4）：1－6（Jin Fengxiang.Study on Robust Estimation Theory and Method［J］.Journal of Shandong University of Science and Technology：Natural Science，2003（4）：1－6）

［5］王振杰，歐吉坤，曲國(guó)慶.抗差估計(jì)的初值選擇［J］.地殼形變與地震，2001（3）：32－35（Wang Zhenjie，Ou Jikun，Qu Guoqing.Choice of Initial Value for Robust Estimation［J］.Crustal Deformation and Earthquake，2001，（3）：32－35）

［6］歐吉坤.一種三步抗差方案的設(shè)計(jì)［J］.測(cè)繪學(xué)報(bào)，1996（3）：173－179（Ou Jikun.Design of a New Scheme of Robust Estimation by Three Steps［J］.Acta Geodaetica et Cartographica Sinica，1996（3）：173－179）

［7］高曉，戴吾蛟.抗差Helmert方差分量估計(jì)在GPS／BDS組合定位中的應(yīng)用［J］.大地測(cè)量與地球動(dòng)力學(xué)，2014，34（1）：173－176（Gao Xiao，Dai Wujiao.Application of Robusthelmert Variance Component Estimation to Position in Combination of GPS and BDS［J］.Journal of Geodesy and Geodynamics，2014，34（1）：173－176）

［8］王福麗，成英燕，韋鋮，等.利用抗差多項(xiàng)式擬合法探測(cè)修復(fù)GNSS周跳［J］.大地測(cè)量與地球動(dòng)力學(xué)，2013，33（3）：129－132（Wang Fuli，Cheng Yingyan，Wei Cheng，et al.GNSS Cycle Slip Detection and Correction Using Robust Polynomial Fitting［J］.Journal of Geodesy and Geodynamics，2013，33（3）：129－132）

［9］武漢大學(xué)測(cè)繪學(xué)院測(cè)量平差學(xué)科組.誤差理論與測(cè)量平差基礎(chǔ)［M］.武漢：武漢大學(xué)出版社，2009（Research Group of Surveying Adjustment in School of Geodesy and Geomatics of Wuhan University.Error Theory and Fundation of Surveying Adjustment［M］.Wuhan：Wuhan University Press，2009）