馬 琳,鄧文彬,張廣泰
(1.新疆大學 建筑工程學院,新疆 烏魯木齊 830046)
抗差估計在變形監測數據處理中的應用
馬 琳1,鄧文彬1,張廣泰1
(1.新疆大學 建筑工程學院,新疆 烏魯木齊 830046)

對抗差估計及其迭代初值進行了研究分析,并且針對不同估計方法所得到的迭代初值,利用Matlab對其進行模擬數據實驗,對其“抗差性”進行比較,結合模擬結果選取“抗差性”較好的迭代初值進行下一步的抗差估計。最后利用某變形監測數據,同時采用最小二乘法和抗差估計兩種方法進行比較,驗證了抗差估計對粗差的“抗干擾性”。
抗差估計;M估計;變形監測;迭代初值;ρ函數
在變形監測工作中,不論采用何種精密測量儀器或測量方法,都有可能出現粗差[1,2]。粗差的出現,會使數據處理分析結果出現偏差,進而給后期的監測工作帶來一定的困擾。在變形監測數據處理中,由于所測變形值很小,與觀測中的誤差極易混淆,這就要求在對數據進行處理分析時,盡可能避免或減小粗差對整體數據解算的干擾,以免影響對變形的正確分析[3]。
抗差估計與傳統的統計方法相比,可以在數據中有粗差存在的情況下,減少粗差影響,從而得到具有“抗干擾性”的解算結果。本文以某工程項目的變形監測數據為基礎,對其進行抗差估計分析,驗證其對粗差的“抗干擾性”。
抗差估計是指在有粗差存在的情況下,選擇合適的估計方法,盡可能降低粗差對未知量估值的影響,得出最優估值[1]。抗差估計基本上可分為3大類:M估計、L估計和R估計。本文中所用到的是M估計。
M估計是經典的極大似然估計的推廣,稱為廣義極大似然型估計。



在對變形監測數據進行處理分析時,常采用線性回歸法,若將抗差M估計與線性回歸法結合在一起,可以有效抵抗粗差的干擾。


令φ(Vi)/Vi=Wi(權因子)為等價權元素,則有,將誤差方程帶入式(5),則有:

對式(6)的求解采用選權迭代法,選擇合適的迭代初值解算參數第一次估值,由其解算出誤差,進而確定新的等價權,得到下一次的參數估值,以此類推,直至前后兩次解的差值符合一定的限差要求,即得到最終的參數估值[1,2,4]。
抗差M估計的關鍵是選擇合適的迭代初值和ρ函數。
2.1 迭代初值
2.1.1 最小二乘估計
目前一般都是采用最小二乘估計值來作為迭代初值,最小二乘估計的原理是最小化殘差平方總和,即但由于最小二乘估計對粗差具有均衡作用,且對粗差不敏感,往往會降低抗差估計的抗差性[5,6]。
2.1.2 一次范數最小估計的線性規劃算法
線性規劃是研究線性約束條件下線性目標函數的極值問題的數學理論和方法,其數學模型為:


把X、V均視為待求參數,由于線性規劃要求所有參數均為非負,而X、V可正可負,故設:

X+與X-、V+與V-互不獨立,且不能同時存在非零解,則得到數學模型為:

即,

求解式(11),則可以得到X的估值[4]。
2.1.3 最小中位數平方
最小中位數平方的原理是最小化殘差平方的中位數,即Median()= min。最小中位數平方估值可以利用重復抽樣算法來求解,其基本步驟為:從n個誤差方程中隨機抽取一個容量為p+1的子樣本,其中p為自變量的維數,解算參數X,并計算p+1組的殘差及其平方的中位數,選擇殘差平方中位數最小值所對應的參數估值為最終的抗差估值[7-9]。
2.2ρ函數
抗差M估計的對粗差的“抗干擾性”主要取決于迭代計算時所選取的ρ函數。一般情況下,ρ函數應盡量滿足以下條件:
1)對稱函數,即ρ(-v)= ρ(v)。
2)ρ(0)=0;ρ(v)在(-∞,0)區間上非增;ρ(v)在(0,-∞)區間上非降。
3)ρ(v)在(-∞,∞)區間上處處連續。
IGG法屬于有淘汰區的M估計,權因子之間變化較平緩,同時這種估計方案充分考慮了測量數據的實際情況,是一種適合處理測量數據的抗差方案[1,10]。本文采用IGG法,其ρ函數及權函數為:

式中,u為標準化殘差(ui= vi/σ);ρ(u)為ρ函數;w (u)為權函數;b、c為調和系數,選取時可以參考有關文獻推薦值[11];d為常數;k為很小的數,避免u為0時出現計算問題[12]。
3.1 數據實驗與分析
針對上文提到的3種迭代初值的計算方法,采用線性回歸模型,利用Matlab進行模擬比較,選取較為合適的迭代初值。
采用模型yi= a1+a2xi+ε,令a1=5.0,a2=-2.9,模擬自變量x服從正態分布N (2,22),ε服從正態分布N(0,0.22),針對模型計算系數a1、a2,對于粗差,模擬分布yi~N(-3,0.52),j為粗差個數,實驗數據為100個,粗差為30個,其模擬結果見表1。

表1 3種不同迭代初值計算結果比較
從表1可以看出,3種迭代初值計算方法中,最小二乘法對粗差的抗干擾能力比較弱,其余兩種方法均可以較好地抵抗粗差。本文采用最小中位數平方估值作為迭代初值,進行下一步抗差估計。
3.2 變形觀測數據的處理及分析
本文選取某工程的變形觀測數據,對其進行抗差線性回歸分析,原始數據如表2。

表2 原始數據
利用Matlab畫出散點圖,依據散點圖選用指數模型y=a×eb/x對其進行抗差線性回歸分析。由于指數模型本身不是線性模型,故先將指數模型線性化:y=a×eb/x?(lny)=(lna)+b×(1/x),然后利用后面線性化后的模型進行數據處理與分析。
若在第4、8、12、16天的下沉量上各加入2、5、7、9 mm的粗差,采用常用的最小二乘法對其進行線性回歸,同時采用以最小中位數平方估值作為迭代初值,以IGG法作為權函數的抗差估計與其進行對比,結果如圖1所示。

圖1 最小二乘法與抗差估計繪制的擬合結果圖
從圖1可以直接看出,在有粗差存在的情況下,采用的最小二乘法所繪制出來的擬合圖已經完全偏離了真實數據,偏向有粗差的數據;而抗差估計計算繪制出來的擬合圖則沒有受到粗差數據的影響,正確預測出了變形數據的趨勢。從表3中的各項指標數據中也可以看出,最小二乘法的解算結果沒有抗差估計的解算結果好。

表3 數據處理結果
常用的最小二乘法在有粗差存在的情況下,所得的估值有一定的偏差,而抗差估計則可以在有粗差存在的情況下,仍然對其進行無偏估計。所以將抗差估計應用在變形監測的數據處理中,不但可以對其進行分析預測,還可以使其在有粗差的情況下,避免粗差對分析預測結果產生干擾。但是在變形監測數據處理中,影響因素往往有多個,而本文中所用到的變形監測數據只有一個影響因素,若將抗差估計應用在有多個影響因素的變形監測數據處理中,則抗差估計的結果可能會更好,這還有待進一步的驗證研究。
[1] 周江文,黃幼才,楊元喜.抗差最小二乘法[M].武漢∶華中理工大學出版社,1997
[2] 楊元喜.抗差估計理論及其應用[M].北京∶八一出版社,1993
[3] 康世英,張宏偉.變形觀測數據處理粗差的定位與剔除[J].桂林工學院學報,2003,23(3)∶310-313
[4] 劉大杰,陶本藻.實用測量數據處理方法[M].北京∶測繪出版社,2000
[5] Andersen R.現代穩健回歸方法[M].上海∶格致出版社,2012
[6] 邱衛寧.具有穩健初值的選權迭代法[J].武漢大學學報∶信息科學版,2003,28(4)∶452-454
[7] 王海娜.線性回歸模型的若干穩健估計方法及應用實例[D].濟南∶山東大學,2013
[8] 楊玲,沈云中,樓立志.基于中位參數初值的等價權抗差估計方法[J].測繪學報,2011,40(1)∶28-31
[9] 王彤,王琳娜,何大衛.基于LMS回歸的一步M估計與加權最小二乘估計[J].現代預防醫學,1996,26(3)∶281-283
[10] 周江文.經典誤差理論與抗差估計[J].測繪學報, 1989,18(2)∶115-120
[11] 李浩軍,唐詩華,黃杰.抗差估計中幾種選權迭代法常數選取的探討[J].測繪科學,2006,31(6)∶70-72
[12] 賈超.穩健回歸分析方法在變形監測中的應用[D].太原∶太原理工大學,2011
P258
B
1672-4623(2016)01-0089-03
10.3969/j.issn.1672-4623.2016.01.026
馬琳,碩士,主要研究方向為大地測量、工程測量。
2014-12-31。
項目來源:國家自然科學基金資助項目(51368056);武漢大學精密工程與工業測量國家測繪地理信息局重點實驗室開放基金資助項目(PF2012-20);新疆大學校院聯合基金資助項目(XY110135)。