李 政，李永樹，楚 彬

（西南交通大學 地球科學與環(huán)境工程學院，四川 成都610031）

穩(wěn)健估計法能夠保證所估的參數(shù)不受或少受模型誤差（首先指的是粗差）的影響，這種方法主要用來發(fā)現(xiàn)粗差和對粗差進行定位［7］。然而一般穩(wěn)健估計權函數(shù)多通過經(jīng)驗法選取，而且權被表示成改正數(shù)的函數(shù)。由于改正數(shù)僅是真誤差的可見部分，所以經(jīng)驗權函數(shù)未顧及平差的幾何條件。事實上，粗差可以視為來自期望為零、方差很大的正態(tài)母體的子樣［8］。根據(jù)整體最小二乘的驗后方差估計，求出觀測值的驗后方差，通過方差檢驗可找出方差異常大（即含粗差）的觀測值；根據(jù)經(jīng)典權與觀測值方差成反比的性質(zhì)賦予它一個相應小的權進行下一步迭代平差，逐步實現(xiàn)粗差定位。

1 驗后方差估計定權

大地測量和攝影測量平差往往包含多組觀測值，每組包含的觀測值具有相同的精度。若假定觀測值互不相關，則可按Forst ner［9］提出的方法求解各組觀測值的驗后方差：

其中

為了發(fā)現(xiàn)每組觀測值內(nèi)的粗差，對第i組內(nèi)任一觀測值li，j求其方差估值＾和相應的多余觀測分量ri，j。根據(jù)式（1）、式（2）可得

由此，建立下列統(tǒng)計量來檢驗該方差是否異常，即相應的觀測值是否包含粗差［10］：

H0假設：E（）＝E），

注:Rpei,Rpej—高低壓繞組單元與油箱之間的絕緣電阻;Cii,Cij—高低壓繞組單元的對地電容;hi,lj—高低壓側的節(jié)點;Ki,Kj—高低壓繞組各單元間的串聯(lián)電容;Rpi,Rpj—高低壓繞組各單元之間的絕緣電阻;Rsi,Rsj—高低壓繞組各單元的歐姆電阻;Chlij—高低壓繞組之間的電容;L—自感;M—互感。uhi-1,ulj-1—hi-1和lj-1節(jié)點的節(jié)點電壓;ihi,ilj—流過高低壓繞組單元阻抗支路的支路電流。

假設觀測值li，j不含粗差，即H0假設成立，則統(tǒng)計量Ti，j近似為自由度為1和ri的中心F分布。若Ti，j＞Fa，1，ri，則表明該觀測值方差與該組觀測值方差有顯著差異，它很可能包含粗差。于是按下列權函數(shù)計算迭代平差中觀測值的權：

由于第一次迭代中，含粗差觀測值的權是不準確的，因此所有的殘差及＾δ0均會受影響，即此時的估計是有偏差的。由于驗后方差在迭代過程中不斷變化，將粗差觀測值的權逐步減小，直至接近于零，最終不影響平差結果，從而使估計從有偏走向無偏。

2 穩(wěn)健整體最小二乘解

在選定合適的權函數(shù)之后，根據(jù)Krar up［10］提出的最小二乘穩(wěn)健估計方法對觀測向量的協(xié)因數(shù)陣Qy進行Cholesky分解：Qy＝Gy，Gy為上三角矩陣。由于在EIV模型中，系數(shù)矩陣中的元素可能為常數(shù)或可確定的元素，其協(xié)方差陣QA的某些行或列的元素可能全為0，因此QA并非為正定矩陣，所以不能采用Cholesky分解。因此本文引入Schur分解定理［11］：假設A為n×n維實對稱矩陣，那么則有n×n維正交矩陣S和對角矩陣U U中對角線上的元素為A的特征值），使得STAS＝U，STS＝In。

將QA進行Schur分解可得QA＝SUST。結合文獻［15］中加權整體最小二乘解（WTLS）和文獻［10］的最小二乘穩(wěn)健估計方法，可得穩(wěn)健整體最小二乘解（RTLS）迭代過程：

1）構造系數(shù)矩陣的協(xié)方差矩陣QA，

2）Qy＝Gy，QA＝SUST，

其中，Py（i）和PA（i）由式（5）確定，（·）－0．5為矩陣·的算術平方根逆根。在迭代過程中，對統(tǒng)計量Ti，j進行假設檢驗，然后根據(jù)式（5）更新下一次迭代的權函數(shù)。通過迭代，含粗差觀測值的權函數(shù)元素值會逐步趨近于0，不含粗差觀測值的權函數(shù)元素變化不大。因此，此方法不僅可以對粗差進行定位，而且所估參數(shù)受粗差影響較小，具有穩(wěn)健性。

3 實驗與討論

在坐標轉換實驗中，由于觀測方程等號兩邊都為觀測值，不可避免會存在隨機誤差，因此采用整體最小二乘方法求解待估參數(shù)更為合理［12］。然而，當觀測值混入粗差時，普通整體最小二乘方法解得的各待估參數(shù)受粗差影響較大，精度較低。為了檢驗本文穩(wěn)健整體最小二乘方法（RTLS）的可行性，分別采用一般最小二乘（LS）、文獻［15］中提出的加權整體最小二乘（WTLS）以及本文的穩(wěn)健整體最小二乘（RTLS）三種方法求解待估參數(shù)，并評價其精度。實驗步驟流程見圖1。

3．1 實驗數(shù)據(jù)

坐標轉換的基本原理是在兩套坐標系下通過足夠數(shù)量同名點求取坐標轉換參數(shù)，然后根據(jù)坐標轉換參數(shù)將某一坐標系下坐標轉換到另一坐標系下。因此，坐標轉換的關鍵就是求解坐標轉換參數(shù)。假設兩套坐標系有以下轉換關系：

假設有n組觀測值時，觀測方程如下：

圖1 實驗步驟流程圖

實驗以文獻［16］數(shù)據(jù)為例，假設有13組觀測值，見表1，通過參數(shù)a1＝0．9，b1＝－0．8，c1＝1，a2＝0．6，b2＝0．7，c2＝5從某一坐標系轉換至另一坐標系。

表1 模擬觀測數(shù)據(jù)

3．2 觀測值僅含有偶然誤差

本文對x0和y0添加ε1∈［－0．3～0．3］的隨機誤差得到組成系數(shù)矩陣的觀測值，對xt和yt添加ε2∈［－0．2～0．2］的隨機誤差構成與系數(shù)矩陣不同精度的觀測向量的值，然后通過編寫程序，分別采用一般最小二乘（LS）、加權整體最小二乘（WTLS）、穩(wěn)健整體最小二乘（RTLS）求得坐標轉換參數(shù)，求解結果見表2。

表2 待估參數(shù)統(tǒng)計表（觀測值僅含偶然誤差）

由表2可知，當觀測向量和系數(shù)矩陣同時含有偶然誤差時，采用加權整體最小二乘法（WTLS）和穩(wěn)健整體最小二乘法（RTLS）所求待估參數(shù)的精度要高于一般最小二乘法（LS）。由于觀測值不含粗差，因此在假設檢驗過程中，H0假設成立，迭代過程中其權函數(shù)并未發(fā)生變化，所以采用加權整體最小二乘法（WTLS）和穩(wěn)健整體最小二乘法（RTLS）獲得的結果相同。

3．3 觀測值含有偶然誤差與粗差

為了驗證穩(wěn)健整體最小二乘法（RTLS）是否有正確定位粗差并且降低粗差對其他觀測數(shù)據(jù)影響的能力，對模擬觀測數(shù)據(jù)加入隨機誤差的同時在第3組數(shù)據(jù)x0與y0中混入3δ0粗差，第7組數(shù)據(jù)xt與yt中混入10δ0粗差，然后分別采用一般最小二乘法（LS）、加權整體最小二乘法（WTLS）和穩(wěn)健整體最小二乘法（RTLS）求解坐標轉換參數(shù)，見表3。

表3 待估參數(shù)統(tǒng)計表（觀測值含偶然誤差和粗差）

由表3可知，在觀測數(shù)據(jù)混入粗差時，無論采用一般最小二乘法（LS）還是加權整體最小二乘法（WTLS）求得的待估參數(shù)都偏離真值較遠。并且所含粗差越大，對其他觀測數(shù)據(jù)影響越大。采用穩(wěn)健整體最小二乘法（RTLS）求解待估參數(shù)，雖然不能完全剔除粗差的影響，但是無論所含粗差大小，其他觀測值受其影響較小 因此本文所述方法具有穩(wěn)健性。為了探究穩(wěn)健整體最小二乘法（RTLS）如何抵抗粗差影響，將迭代過程中含粗差的觀測值權值變化情況繪制成曲線，如圖2所示。

圖2 含粗差觀測值權變化曲線圖

圖2 中含有粗差的觀測值權函數(shù)數(shù)值變化較快，隨著迭代過程，其值越來越小。平差過程就是誤差分配的過程，穩(wěn)健整體最小二乘法（RTLS）根據(jù)觀測值的驗后方差，通過假設檢驗的方法檢測方差異常大的觀測值，并在下次迭代平差中賦予一個相應小的權值，逐步實現(xiàn)粗差定位。由于含粗差觀測值的權值較小，其對待估參數(shù)求解“貢獻”也較小，因此待估參數(shù)能盡可能少的受到來自粗差的污染。

4 結 論

由上述實驗可知，當觀測向量與系數(shù)矩陣同時存在偶然誤差時，采用加權整體最小二乘法（WTLS）以及本文的穩(wěn)健整體最小二乘法（RTLS）求解待估參數(shù)效果好于一般最小二乘法（LS）；當觀測向量和系數(shù)矩陣同時混入粗差，采用一般最小二乘法（LS）以及加權整體最小二乘法（WTLS）求解待估參數(shù)時，改正數(shù)受粗差影響較大，所求參數(shù)嚴重偏離真值；而采用驗后方差估計進行定權，能保證權函數(shù)盡可能少受粗差污染，在得到合適的權函數(shù)之后采用穩(wěn)健整體最小二乘法（RTLS）求解待估參數(shù)，其參數(shù)估值受粗差影響較小，具有穩(wěn)健性。

［1］ 陶本藻，邱衛(wèi)寧．誤差理論與測量平差［M］．武漢：武漢大學出版社，2012．

［2］ 劉經(jīng)南 曾文憲 徐培亮．整體最小二乘估計的研究進展［J］．武漢大學學報：信息科技版，2013，38（5）：505－512．

［3］ 梁景輝．基于IGGⅢ穩(wěn)健估計法在大壩沉降監(jiān)測中的應用［J］．測繪與空間地理信息，2014，37（7）：196－197．

［4］ GOLUB G H，VAN L CH．An Analysis of the Total Least Squares Problem［J］．SIA M Jour nal on Nu merical Analysis，1980，17（6）：883－893．

［5］ SCHAFFRIN B，ANDREAS W．On Weighted Total Least squares Adjust ment f or Linear Regression［J］．Journal of Geodesy，2008，82（7）：415－421．

［6］ NEITZEL F．Generalization of Total Least－Squares on Example of Unweighted and Weighted 2D Si milarity Transfor mation［J］．Jour nal of Geodesy，2010，84（12）：751－762．

［7］ CHOI Y J，KI M J Y，SUMG K M．A Robust Algorith m of Total Least Squares Met hod［J］．IEICETRANS，F(xiàn)undamentals，1997，7（E80－A）：1336－1339．

［8］ 李德仁，袁修孝．誤差處理與可靠性理論［M］．武漢：武漢大學出版社，2002．

［9］ FORSTNER W．Uncertain Geo metric Relationships and Their Use f or Object Location in Digital Images［M］．The Tutorial Notes f or The Second Inter national Wor kshop on Robust Co mputer Vision，Bonn，1992．

［10］KRARUP T，JUHL J，KUBIK K．Gotterdammerung over least squares adjust ment［J］．Pr oceedings of 14th Congress of the International Society of Photogrammetry，1980，B3：369－378．

［11］卜長江，羅躍生．矩陣論［M］．哈爾濱：哈爾濱工程大學出版社，2007．

［12］楚彬，范東明．基于比例整體最小二乘的GPS高程擬合［J］．測繪工程，2014，23（4）：37－39．

［13］李成仁，岳東杰，袁豹，等．基于最小二乘配置的九參數(shù)模型在三維坐標轉換中的應用［J］．測繪與空間地理信息，2014，37（7）：193－195．

［14］馮劍橋，黃張裕，徐秀杰，等．總體最小二乘法在坐標轉換中的應用［J］．測繪與空間地理信息，2014，37（7）：205－206．

［15］MAHBOUB V．On weighted total least－squares f or geodetic transfor mations［J］．Jour nal of Geodesy，2012，86（5）：359－367．

［16］MAHBOUB V，A MIRI－SI MKOOEI A R，SHARIFI M A．Iteratively reweighted total least squares：a robust esti mation in errors－in－variables models［J］．Survey Review，2013，45（329）：92－99．