非線性模型的抗差最小二乘法及其在圓心擬合中的應(yīng)用

李春雷,王亞男,葛仁磊,陳珣

(海洋石油工程股份有限公司,山東 青島 266555)

圓形及圓柱是最常見的自然曲面之一,在各類工程中得到了非常廣泛的應(yīng)用,使用測量數(shù)據(jù)計算結(jié)構(gòu)物的圓心坐標也是工程及工業(yè)測量中最基本的任務(wù)之一。在實際測量工作中,粗差往往不可避免,但在圓心計算中,即使錯誤數(shù)據(jù)只占所有測量數(shù)據(jù)很小的一部分,也會導致計算結(jié)果產(chǎn)生較大的偏離。

解決此問題最常用的方法是對測量數(shù)據(jù)進行人工檢查,剔除其中的粗差點,然后再進行后續(xù)計算。但此過程由于需要人為干預(yù),導致效率較低,而且由于不同的數(shù)據(jù)處理人員的主觀意識影響,也會使不同的數(shù)據(jù)處理人員得到的數(shù)據(jù)處理結(jié)果各不相同,難以統(tǒng)一。

抗差最小二乘法就是為了解決粗差數(shù)據(jù)對處理結(jié)果影響過大的問題而提出的?？共钭钚《朔ㄗ罨镜脑硎窃谧钚《朔ǖ幕A(chǔ)上,使用測量數(shù)據(jù)自動產(chǎn)生等價權(quán),并根據(jù)不同的原則調(diào)整后驗權(quán),使粗差點對結(jié)果的影響減弱甚至消除精差點的影響。

1 選權(quán)迭代法

最小二乘原則在測量平差中得到了非常廣泛的應(yīng)用。但是,當觀測值存在粗差時,最小二乘原則就會出現(xiàn)明顯的不適應(yīng)性,即單個測值的偏差即可能導致結(jié)果面目全非[1]。基于此,抗差最小二乘法應(yīng)運而生?？共钭钚《朔ㄊ窃谧钚《朔ǖ幕A(chǔ)上,通過改變觀測值的權(quán)重,以減少甚至消除觀測值中粗差的影響。

進行測量數(shù)據(jù)平差時最常用的方法是間接平差法,間接平差的誤差方程為[2]:

V=AX-L

根據(jù)條件平差理論,不僅可以對未知參數(shù)X進行估計,還可以對未知參數(shù)X的后驗權(quán)進行估計,X及其后驗權(quán)陣為PXX的估計值為[2]:

X=(ATPA)-1ATPL

PXX=ATPA

假設(shè)觀測值中有一個粗差i,根據(jù)中誤差σi與權(quán)pi的關(guān)系[2]:

在一次計算中,單位權(quán)中誤差σ0固定,也就是說觀測值的方差越大,則其權(quán)就越小。也就是說使用后驗權(quán)陣PXX當作權(quán)陣再進行一次平差,粗差i對結(jié)果的影響就會變小,未知參數(shù)X的估計就會更接近真值。如果我們反復迭代,那么粗差的影響會越來越小,未知參數(shù)X的估計也會越來越接近真值。

上述做法就叫做選權(quán)迭代法[3-4],其選權(quán)原則為最小二乘法原則,即:

上式中:pi為PXX的第i個對角素;νi為向量V中第i個元素。

顯然,對于最小二乘原則:

wi=1

2 常用選權(quán)原則

最小二乘原則是抗差最小二乘法中最常用的一種原則,第1節(jié)中已經(jīng)作了詳細介紹,這里就不作贅述。下面就繼續(xù)為大家介紹兩種常用原則。

2.1 殘差絕對和最小原則

殘差絕對和最小原則的估計準則為:

相應(yīng)的權(quán)因子為[3]:

在實際計算中,為了避免由于νi≈0時造成計算錯誤,一般使用下式來計算:

c一般取比較小的值,在實際使用中常取c=0.001。

殘差絕對和最小原則的抗差性非常好,其崩潰污染率高達50%,是所有抗差估計類中,是崩潰污染率最大的一種。但其精度并不高,所以常用于初始值的估計[3]。

2.2 IGG Ⅲ選權(quán)原則[4]

IGG為中國科學院測量與地理研究所的英文縮寫,IGG Ⅲ方案是基于測量誤差的有界性提出來的,它對測量抗差估計比較有效。其等價權(quán)因子取為[4]:

式中,σ0為單位權(quán)中誤差。由圖1wi與νi/σ0的關(guān)系可以很明確地看出k0和k1的幾何意義。

圖1 wi與νi/σ0的關(guān)系圖

在抗差最小二乘法中最常用的取值為k0=1.5,k1=2.5。但其取值并不是固定不變的,在精密工程測量與工業(yè)測量中,可以根據(jù)實際工程需要對k0和k1進行設(shè)置。而且當k0=k1=3時,此等價權(quán)設(shè)計方案就變成常用的萊因達準則。

3 三步抗差方案[5]

學者利用上述幾種選權(quán)方案的優(yōu)點,提出了“三步抗差方案”,之前的討論多在線性模型的基礎(chǔ)上進行的[1,5]。但在實際數(shù)據(jù)處理工作中,非線性模型則更為常見。非線性模型的引入,會加大大增加三步抗差的復雜度。在本文根據(jù)工程測量實際,對三步抗差方案進行適當?shù)暮喕?并介紹非線性模型的三步抗差方案。

同線性模型一樣,三步抗差方案分為初值選取階段、粗差剔除階段及結(jié)果優(yōu)化階段。

3.1 初值選取階段

此階段的目標是選擇抗差性能優(yōu)良的初始值。根據(jù)前文的探討,此階段使用殘差絕對和最小原則進行初值的選取。選擇過程如圖2所示。

圖2 初值選取過程

3.2 粗差剔除階段

此階段的目標是使粗差充分暴露。一般通過選權(quán)迭代法將粗差的權(quán)重逐漸降低?？梢赃x用IGG Ⅲ選權(quán)原則進行此階段的計算。計算過程如圖3所示。

圖3 粗差剔除過程

3.3 最優(yōu)結(jié)果階段

此階段的目的是重新進行權(quán)矩陣評估,進行精度估計,得到最優(yōu)結(jié)果。過程如圖4所示。

圖4 最優(yōu)結(jié)果階段過程

在實際測量數(shù)據(jù)中,一般只將數(shù)據(jù)分成準確數(shù)據(jù)和粗差數(shù)據(jù)兩部分。準確數(shù)據(jù)直接使用先驗權(quán),粗差數(shù)據(jù)直接剔除。在這一步初廣泛接受的是萊因達準則,本文也使用萊因達準則。

4 圓心擬合

在工程測量數(shù)據(jù)處理過程中,可以將空間圓心擬合轉(zhuǎn)化為平面圓心擬合[7]。不失一般性,本文所述的圓心擬合都指平面圓心擬合。

以(x0,y0)為圓心,以R為半徑的圓可以表示為[8]:

(x-x0)2+(y-y0)2=R2

對于觀測點(xi,yi)列出其誤差方程為:

由圓的幾何意義可知,νi相當于點i與圓弧之間的距離。

根據(jù)文獻[2]將其線性化:

為了便于對比,引用文獻[9]圓柱體底部各點坐標數(shù)據(jù),經(jīng)三步平差方案擬合后,圓心為(49.986 2,60.012 8),半徑為0.430 4,單位權(quán)中誤差σ0=0.001 89。每點殘差及中誤差如表1所示。

表1 圓柱體底部各點坐標 單位:m

與文獻[9]比較可以發(fā)現(xiàn),本文方法與其穩(wěn)健估計的結(jié)果基本一致。都探測出了兩個粗差點(4號點、11號點)。

為了進一步測試三步平差方案的穩(wěn)定性及健壯性,使用程序自動生成圓心為(0,0),半徑為0.5的測試數(shù)據(jù),并用三步平差方案進行擬合,以便對抗差方案進行分析。測量數(shù)據(jù)(xi,yi)生成方法如下:

αi=i×18+εi

ri=0.5+Δi

上式中:α為圓心到圓上的點與x軸的夾角,i∈[0,1,2,…,17];ε為從-4～+5的隨機整數(shù);Δ為均值為0,標準差為0.001的,符合正態(tài)分布的隨機數(shù)。

為了模擬測量數(shù)據(jù)中存在的粗差,在上述的生成的測量數(shù)據(jù)中增加α為0～2π的隨機數(shù),半徑r為0～4.996或0.504～1的隨機數(shù)的粗差。

用上述方法共生成60 000個測試數(shù)據(jù)(每一組用數(shù)據(jù)分別有1～30個粗差,共200組數(shù)據(jù)),對所有數(shù)據(jù)進行三步平差,結(jié)果統(tǒng)計如表2所示(只截取部分反映趨勢的結(jié)果)。

表2 結(jié)果統(tǒng)計

表2中第一階段污染表示其中誤差大于等于0.01,第一階段崩潰表示其中誤差大于等于0.02。第二階段污染表示其中誤差大于等于0.005,第二階段崩潰表示其中誤差大于等于0.01。結(jié)果污染表示其中誤差大于等于0.003,結(jié)果段崩潰表示其中誤差大于等于0.006。由上述結(jié)果可知,此方法的崩潰污染率接近50%(當污染率為50%時,崩潰率只有0.5%),但當污染率超過50%時,崩潰率迅速增加。隨著崩潰率的增加,會逐漸出現(xiàn)結(jié)果不收斂的情況。這是由于非線性模型的線性化是基于非線性函數(shù)的泰勒級數(shù)展開并忽略其高階無窮小。但由于第一階段的崩潰,無法提供給后續(xù)數(shù)據(jù)較好的初值,使舍去的高階無窮小量足以影響計算結(jié)果,導致算法不收斂。

5 結(jié)論

抗差最小二乘法能夠很好地減少測量數(shù)據(jù)中粗差的影響。如果使用三步抗差方案,在測量數(shù)據(jù)良好的情況下,崩潰污染率也高達50%,可以很好地應(yīng)對工程測量數(shù)據(jù)處理中出現(xiàn)的各種問題。而且,經(jīng)過第三階段的結(jié)果優(yōu)化,測量數(shù)據(jù)的精度也可以得到保證。