穩(wěn)健回歸分析在GPS高程擬合中的應(yīng)用

劉 振，吳長悅

(華北理工大學(xué) 礦業(yè)工程學(xué)院，河北 唐山 063210)

GPS是隨著現(xiàn)代衛(wèi)星技術(shù)在測繪領(lǐng)域的普及而產(chǎn)生的一種空間定位技術(shù)，具有很高的測量精度，不受時間、空間地域的限制，成本較低等優(yōu)點，被廣泛地應(yīng)用在地質(zhì)工程、交通運輸、測量等領(lǐng)域[1]。通常情況下，利用全球定位系統(tǒng)(GPS)得到的海拔是WGS-84下的大地高，在實際項目施工時，GPS測得的大地高是不能直接利用的，需要通過高程異常值求得海拔[2]，高程的異常值是大地高與正高之間作差而得到的，對于大部分工程項目，由于施工區(qū)域面積較小，高程異常的變化比較緩慢，利用普通的多項式擬合就能夠滿足施工精度的需求[3]。高程擬合的實際機理就是把高程的異常值與地面坐標(biāo)近似地看作一個多項式內(nèi)部的關(guān)系，運用已知公共點的大地高和正常高構(gòu)成誤差方程，根據(jù)最小二乘法的原理(LS)求解多項式的系數(shù)。最小二乘法規(guī)定所求殘差的平方和最小，它重點是針對觀測過程中偶然產(chǎn)生的誤差起作用[4]。在實際施工測量中不可避免地存在較大誤差,由于傳統(tǒng)的最小二乘法不具備抵抗較大誤差的能力，所以當(dāng)樣本存在較大誤差時，一定會對整個回歸方程的系數(shù)估計構(gòu)成嚴(yán)重影響，因抗差估計具有抵抗較大誤差的特性，故本文把抗差估計引入GPS高程擬合，以達到提高精度的目的。

1 高程擬合的原理

我國在施工測量中，通常是利用正高作為高程值，地球表面上的點沿鉛垂方向到與地球近似的橢球表面的長度叫作正常高Hr，在本實驗中，正常高是利用二等水準(zhǔn)得到的[5]；地表上的點沿法線的方向延伸到參考橢球面的長度稱為大地高H84，是利用GPS點連結(jié)成網(wǎng)，經(jīng)過平差處理得到的以WGS-84為橢球基準(zhǔn)面的高，這就有必要找出GPS 點的大地高H84與正常高Hr之間的關(guān)系，并將H84轉(zhuǎn)換為Hr。大地高與正常高的關(guān)系如圖1所示，其中，高程異常值用似大地水準(zhǔn)面到橢球面之間的高差ξ來表示[6]。由圖1可知，要求得某點的GPS正常高Hr，只需該點的高程異常值ξ和它的大地高H84。

圖1 大地高與正常高的關(guān)系

Hr=H84-ξ

(1)

把水準(zhǔn)測量和GPS測量相結(jié)合，是當(dāng)前GPS施工測量中頻繁使用的手段之一，當(dāng)前世界上用GPS和水準(zhǔn)相結(jié)合求得高程的主流方法為曲面擬合法[7]。

2 曲面擬合法原理

曲面擬合法作為一種多項式擬合的方法，把高程的異常值ξ近似看作因變量y,測區(qū)內(nèi)的各個坐標(biāo)看作自變量x，然后依據(jù)測量范圍內(nèi)已知的平面坐標(biāo)(x,y)和ξ值組成方程組，通過解方程就可以求得函數(shù)的擬合系數(shù)，從而擬合出所測區(qū)域內(nèi)的似大地水準(zhǔn)面，再內(nèi)插未知點平面坐標(biāo)解得未知點的ξ，求出未知點的正常高Hr[8]。

設(shè)測站點的高程異常值ξi與坐標(biāo)(xi,yi)之間存在以下函數(shù)關(guān)系:

ξi=f(xi,yi)+εi

(2)

式中,(xi,yi)為ξi的近似值；εi為擬合殘差。

經(jīng)過研究大量有關(guān)GPS高程擬合的相關(guān)資料，最終得到一致的結(jié)論：在所有的擬合方法中二次曲面擬合的結(jié)果精度要高于其余的擬合方法，所以本文直接利用二次曲面擬合。

二次曲面擬合方程可表述為：

ξ=a0a1x+a2y+a3x2+a4xy+a5y2

(3)

相應(yīng)的誤差方程為：

(4)

如果有n個水準(zhǔn)測量點，用矩陣可表示為：

b=Aβ

(5)

式中,

b=[ξ1,ξ2,...,ξn]T

β=[a0,a1,a2,a3,a4,a5]T

由上式可知，系數(shù)矩陣A除第一列是常數(shù)項之外，其它項都是由觀測值構(gòu)成，所以，誤差一定會存在于觀測向量b和系數(shù)矩陣A之中，使用傳統(tǒng)的最小二乘法已經(jīng)不能滿足需求，引入穩(wěn)健回歸分析[9]：

b-ε=(A-σ)β

(6)

式中，b為含隨機誤差ε的觀測向量；A為含隨機誤差δ的n×m維系數(shù)矩陣；β為m維待估參數(shù)向量。

3 實驗分析

3.1 數(shù)據(jù)介紹

實驗數(shù)據(jù)選取華北理工大學(xué)18級研究生空間大地測量學(xué)第四和第二測量小組的靜態(tài)GPS數(shù)據(jù)和二等水準(zhǔn)數(shù)據(jù)。將11個已知點分為兩組，其中8個點作為擬合點，如表1所示，進行二次曲面的多項式擬合，剩下3個作為檢核點，如表2所示，用來對擬合結(jié)果精度的檢核，從而判斷所建立回歸方程的準(zhǔn)確性和有效性。

表1 擬合點數(shù)據(jù)/m

表2 檢核點數(shù)據(jù)/m

3.2 算法流程

高程擬合的具體流程如圖2所示。

圖2 高程擬合流程圖

3.3 在MATLAB中進行驗證

3.3.1 高斯-牛頓迭代法在MATLAB中的實現(xiàn)

在MATLAB中，提供了nlinfit函數(shù)用于求解高斯-牛頓法的非線性最小二乘數(shù)據(jù)擬合[7]。函數(shù)的調(diào)用格式如下：

[beta,r,J,COVB,mse]=nlinfit(X,y,fun,beta0):同時返回的beta為擬合系數(shù)，r為殘差，J為jacobi矩陣，COVB為評估的協(xié)方差矩陣，mse為誤差的方差。輸入?yún)?shù)beta0為初始預(yù)測值[11]。

[...]=(X,y,fun,beta0,options):指定控制參數(shù)后返回值。參數(shù)options包括MaxIter、TolFun、TolX、Display、DerivStep等。

當(dāng)X為矩陣時，則X的每一列為自變量的取值，y是一個相應(yīng)的列向量。如果fun中使用了@，則表示函數(shù)的句柄。

核心代碼如下所示：

function F=myfunn1(beta,x)

z1=p2-q2

x=x1-mean(x1);y=y1-mean(y1);

xi=x2-mean(x2);yi=y2-mean(y2);

myfunn1=inline('beta(1)+beta(2)*x(:,1)+beta(3)*x(:,2)+beta(4)*x(:,1).*x(:,1)+beta(5)*x(:,1).*x(:,2)+beta(6)*x(:,2).*x(:,2)','beta','x');

z=p1-q1;

X=[x y];

[beta,r,J,COVB,mse]=nlinfit(X,z,myfunn1,[1 1 1 1 1 1])

zi=beta(1)+beta(2)*xi+beta(3)*yi+beta(4)*xi.*xi+beta(5)*xi.*yi+beta(6)*yi.*yi

zj=beta(1)+beta(2)*x+beta(3)*y+beta(4)*x.*x+beta(5)*x.*y+beta(6)*y.*y

z2=[z1-zi]

z3=z-zj

K=z3'*z3

V=z2'*z2

N=sqrt(K/7)

M=sqrt(V/2)

3.3.2 抗差估計在MATLAB中的實現(xiàn)

可以利用MATLAB中已有的robustfit函數(shù)來做抗差估計[12]，一般使用方式為：

B=robustfit(X,y)，輸入?yún)?shù)X是已知點坐標(biāo)所組成的設(shè)定好的矩陣，它是一個n*p的矩陣，輸入?yún)?shù)y為已知點的高程異常值所組成的矩陣，是一個n*1的列向量[13]，通常采用如下調(diào)用格式：

[b,statas]=robustfit(X,y,wfun,tune):用參數(shù)wfun指定加權(quán)函數(shù)，用參數(shù)tune指定調(diào)節(jié)常數(shù)[14]。

核心代碼如下：

x=x1-mean(x1);y=y1-mean(y1);% 數(shù)據(jù)中心化

xi=x2-mean(x2);yi=y2-mean(y2);% 數(shù)據(jù)中心化

z=p1-q1; % 已知點高程異常矩陣

X=[x y x.*x x.*y y.*y];

[b stats]=robustfit(X,z)

B=[1 1 1 1 1 1 1 1]'

C=horzcat(B,X)

z2=C*b

z3=z-z2%擬合殘差

D=[1 1 1]'

E=[xi yi xi.*xi xi.*yi yi.*yi]

X1=horzcat(D,E)

z4=X1*b

z5=z1-z4%檢核殘差

subplot(2,1,2);plot(z5,'g')

ylabel('殘差(m)')

title('圖(b)檢核點殘差圖')

grid on;

subplot(2,1,1);plot(z3,'g')

xlabel('點')

ylabel('殘差(m)')

title('圖(a)擬合點殘差圖')

grid on;

N=sqrt((z3'*z3)/7)

M=sqrt((z5'*z5)/2)

3.4 實驗結(jié)果及精度評定

通過MATLAB程序的運行可以得到殘差圖，如圖3、圖4所示。

圖3 高斯-牛頓迭代法得到的殘差圖

圖4 穩(wěn)健回歸分析得到的殘差圖

實驗結(jié)果如表3所示。A0～A5是擬合得到的二次曲面的系數(shù)。為了更好的進行精度評定，表中列出兩種方法的擬合點和檢核點的最大最小殘差來進行對比，并引入內(nèi)外符合精度作為評價擬合效果的指標(biāo)。

表3 實驗結(jié)果

3.5 實驗結(jié)果分析

從殘差圖來看，高斯-牛頓迭代法與穩(wěn)健回歸分析法的擬合點殘差都有波動，但波動范圍較小，最大殘差均不超過1.5 cm,檢核點殘差均不超過5 cm，但穩(wěn)健回歸分析的最大殘差相對較小，不超過3 cm，高斯-牛頓迭代法的最大殘差已經(jīng)超過4 cm。由表3可知，用高斯-牛頓迭代法進行二次曲面擬合與用穩(wěn)健回歸分析擬合法對比如下：擬合系數(shù)差異較大，前者A0是-354.220 465，后者A0是0；擬合點的殘差總體相差不大；檢核點的殘差，用穩(wěn)健回歸分析的最大殘差和最小殘差都相對較??；內(nèi)符合精度相差不大；外符合精度穩(wěn)健回歸數(shù)值是0.021 1，相對較小，說明精度較高。穩(wěn)健回歸相比高斯-牛頓迭代法在高程擬合方面誤差相對較小，具有一定的優(yōu)勢，也在一定程度上說明穩(wěn)健回歸分析用于高程擬合的現(xiàn)實可行性。

4 結(jié)論與不足

本文采用高斯-牛頓迭代最小二乘法與穩(wěn)健回歸分析法進行高程擬合，通過精度檢驗和對比分析，不僅說明在地勢比較平坦的地區(qū)用GPS水準(zhǔn)進行高程擬合基本可以滿足精度要求，還證明了在觀測數(shù)據(jù)存在粗差的情況下，用穩(wěn)健回歸能夠在一定程度上提高擬合精度，具有現(xiàn)實可行性。但是由于條件的限制，實驗數(shù)據(jù)量較少，高程擬合的測區(qū)范圍也不大，實驗結(jié)果存在一定的局限性，可以通過大區(qū)域和更多的數(shù)據(jù)來進行更深入的研究。