EGM2008和加權(quán)整體最小二乘在GPS高程擬合中的應(yīng)用

范昆飛，易桂軒，孔建，劉立

(1.南寧市國土測繪地理信息中心，廣西 南寧 530021； 2.武漢市測繪研究院，湖北 武漢 430000；3.武漢大學(xué)，湖北 武漢 430079； 4.浙江省第一測繪院，浙江 杭州 310012)

EGM2008和加權(quán)整體最小二乘在GPS高程擬合中的應(yīng)用

范昆飛1*，易桂軒2，孔建3，劉立4

(1.南寧市國土測繪地理信息中心，廣西 南寧 530021； 2.武漢市測繪研究院，湖北 武漢 430000；3.武漢大學(xué)，湖北 武漢 430079； 4.浙江省第一測繪院，浙江 杭州 310012)

在利用地球重力場模型EGM2008結(jié)合GPS/水準(zhǔn)數(shù)據(jù)計算得到剩余高程異常的基礎(chǔ)上，利用加權(quán)整體最小二乘(WTLS)方法對剩余高程異常進(jìn)行三次曲面擬合，獲得更加精確的GPS高程擬合模型。通過對CORS系統(tǒng)GPS/水準(zhǔn)數(shù)據(jù)的擬合和外部精度檢核表明，綜合EGM2008和WTLS方法能夠顯著提高GPS高程異常擬合精度。

EGM2008；WTLS；高程異常；三次多項式

1 引 言

GPS觀測技術(shù)具有高效性、全球性、高精度等優(yōu)點，在道路交通、工程建設(shè)等諸多領(lǐng)域獲得了廣泛應(yīng)用，尤其在測繪領(lǐng)域引起了革命性的變化。傳統(tǒng)水準(zhǔn)測量獲得高程的方法費時費力，且容易受到觀測環(huán)境的影響。GPS測量輕便、高效，觀測點之間相互獨立，沒有誤差積累，相對于傳統(tǒng)水準(zhǔn)測量有很大優(yōu)勢。但是GPS觀測獲得大地高為相對于地球橢球的幾何高程，需要將其轉(zhuǎn)換為正常高[1，2]。最常用的高程擬合方法為將高程異常描述為曲面函數(shù)，利用GPS/水準(zhǔn)點的高程異常數(shù)據(jù)和平面坐標(biāo)，采用最小二乘方法進(jìn)行曲面函數(shù)系數(shù)求解[3～5]。但是傳統(tǒng)最小二乘方法只對觀測值進(jìn)行改正，當(dāng)觀測向量和系數(shù)矩陣由于函數(shù)模型、儀器誤差等影響含有誤差時，最后求解結(jié)果是有偏的。加權(quán)整體最小二乘(Weighted Total Least Squares)以觀測值和系數(shù)矩陣殘差平方和最小為準(zhǔn)則，同時對觀測向量和系數(shù)矩陣進(jìn)行改進(jìn)，能夠改正函數(shù)模型存在誤差的問題。同時通過迭代的方法確定觀測向量和系數(shù)矩陣的權(quán)，有效解決了GPS觀測數(shù)據(jù)和水準(zhǔn)數(shù)據(jù)精度等級不同的問題[6～9]。

國際上先后實施了CHAMP、GRACE、 GOCE衛(wèi)星重力計劃。高精度、高時空分辨率的衛(wèi)星重力數(shù)據(jù)極大地完善了地球重力場數(shù)據(jù)。美國國家地理空間情報局(NGA)綜合利用衛(wèi)星重力數(shù)據(jù)、地面重力數(shù)據(jù)、數(shù)字地面模型等，研制新一代超高階地球重力場模型EGM2008。利用高精度地球重力場模型能夠精確計算重力場短波分量，彌補GPS水準(zhǔn)擬合在反映地形起伏方面的不足。綜合利用加權(quán)整體最小二乘和EGM2008，對實際GPS/水準(zhǔn)數(shù)據(jù)進(jìn)行區(qū)域似大地水準(zhǔn)面擬合，擬合結(jié)果精度表明該方法能夠有效提高GPS高程精度。

2 GPS高程擬合與加權(quán)整體最小二乘

2.1 三次曲面函數(shù)

在對區(qū)域GPS高程進(jìn)行擬合時，通常采用二次多項式函數(shù)對GPS點的高程異常ζ與平面坐標(biāo)(x，y)進(jìn)行描述，具體函數(shù)關(guān)系為

(1)

而文獻(xiàn)[10]研究表明，三次多項式曲面函數(shù)在GPS高程擬合中精度更高。在現(xiàn)有計算機(jī)計算能力條件下，二次多項式和三次多項式計算耗費時間都十分微小，因此本文采用三次多項式對GPS高程進(jìn)行擬合，具體函數(shù)表達(dá)式為：

(2)

設(shè)在該區(qū)域內(nèi)有n個GPS水準(zhǔn)點，則曲面函數(shù)的矩陣形式為：

l=Aβ

(3)

其中l(wèi)為高程異常ζ觀測向量，A為GPS/水準(zhǔn)點的平面坐標(biāo)(xi，yi)組成的設(shè)計矩陣，β為三次多項式待求系數(shù)。

矩陣的具體形式為：

(4)

2.2 加權(quán)整體最小二乘(WTLS)

由于觀測方程中，觀測向量l和設(shè)計矩陣A中的誤差不可避免，因此引入EIV模型[11]：

l-ε=(A-σ)β

(5)

其中l(wèi)為含有隨機(jī)誤差ε的高程異常觀測向量；A為含有隨機(jī)誤差σ的設(shè)計矩陣；β為三次多項式待求系數(shù)。加權(quán)整體最小二乘計算中，設(shè)計矩陣A的權(quán)陣為：

QA=Q0?Qx

(6)

WTLS的迭代步驟為：

3 EGM2008地球重力場模型

美國空間情報局2009年發(fā)布了超高階地球重力場模型EGM2008，該模型階次達(dá)到了 2 159階，空間分辨率約為5′×5′。由地球重力場模型EGM2008可以計算地球上任意一點的模型高程異常：

(7)

利用已知GPS/水準(zhǔn)點數(shù)據(jù)可以計算得到水準(zhǔn)高程異常ζ，EGM2008可以計算得到GPS/水準(zhǔn)點的模型高程異常ζG，則剩余高程異常

ζres=ζ-ζG

(8)

通過三次曲面函數(shù)采用加權(quán)整體最小二乘方法對剩余高程異常ζres進(jìn)行擬合，得到區(qū)域剩余高程異常模型，與EGM2008計算得到的模型高程異常ζG結(jié)合，最終獲得該區(qū)域高程異常模型。

4 實例分析

本文選取某市CORS站點數(shù)據(jù)對綜合加權(quán)整體最小二乘和EGM2008方法的有效性和精度進(jìn)行驗證。通過后期增加水準(zhǔn)聯(lián)測工作，該區(qū)域內(nèi)聯(lián)測水準(zhǔn)的GPS控制點共25個，經(jīng)過平面坐標(biāo)標(biāo)準(zhǔn)化處理后，分布情況如圖1所示。為檢核上述方法的精度，選取其中9個點作為外部檢核點(黑色三角)，利用剩余16個點(空心圓圈)進(jìn)行GPS高程擬合。從圖1可以看出，該區(qū)域的GPS/水準(zhǔn)點分布較為均勻，檢核點主要分布在區(qū)域內(nèi)部，不會因為檢核點位置影響所得高程異常擬合模型檢核精度。

為了檢驗加權(quán)整體最小二乘對設(shè)計矩陣A含有誤差時的有效性，在擬合點中，D11和D24平面坐標(biāo)加入了粗差，具體坐標(biāo)和高程異常值如表1所示。

圖1 GPS/水準(zhǔn)點分布

含有粗差的GPS/水準(zhǔn)數(shù)據(jù) 表1

分別利用傳統(tǒng)LS方法、WTLS、EGM2008+WTLS三種方法對16個GPS/水準(zhǔn)點進(jìn)行高程擬合，得到該區(qū)域的高程異常模型，然后將9個檢核點平面坐標(biāo)帶入各自模型進(jìn)行模型精度檢驗。表2給出了三種擬合方法所得的三次多項式擬合參數(shù)。

曲面擬合參數(shù) 表2

圖2和圖3分別給出了LS和WTLS方法GPS高程擬合結(jié)果的三維視圖，圖中標(biāo)出了粗差點D11和D24的位置。從圖中可以看出，D11由于加入了人為粗差，位置偏離到了邊緣位置，而該位置周圍點的高程異常值都非常小。從圖2與圖3對比中可以看出，在利用LS方法進(jìn)行GPS高程擬合時，由于粗差點D11的存在，該區(qū)域的高程異常被明顯向上拉升而偏離了周圍GPS/水準(zhǔn)點的高程異常值，導(dǎo)致了擬合結(jié)果與實際情況有較大偏差。粗差點D24同樣導(dǎo)致了擬合結(jié)果的偏差。而WTLS方法綜合考慮了GPS/水準(zhǔn)點的高程異常和平面坐標(biāo)的誤差，對GPS/水準(zhǔn)點進(jìn)行迭代定權(quán)，消除了粗差對擬合結(jié)果的影響，建立了更加合理的擬合模型。從圖3中可以看出，粗差點D11和D24并沒有影響周圍高程異常的擬合結(jié)果，并且其他部分的高程異常也更加貼合實際高程異常值。

圖4給出了EGM2008+WTLS的GPS高程擬合結(jié)果的三維視圖，從圖中可以看出，與WTLS擬合結(jié)果相比，EGM2008+WTLS也沒有受到粗差點的影響，并且由于EGM2008為檢核點提供了更加精確的模型高程異常，使得擬合模型更加符合區(qū)域地形情況。其中區(qū)域中部起伏更加明顯，整體趨勢與檢核點更加貼合。

圖2 LS方法GPS高程異常擬合結(jié)果

圖3 WTLS方法GPS高程異常擬合結(jié)果

圖4 WTLS+EGM2008方法GPS高程異常擬合結(jié)果

精度統(tǒng)計 表3

利用LS、WTLS、EGM2008+WTLS方法得到區(qū)域高程異常模型后，將9個檢核點的平面坐標(biāo)帶入模型，獲得這些點的擬合高程異常，與已知的水準(zhǔn)高程異常做差得到模型的擬合誤差。表3給出了模型擬合的精度統(tǒng)計結(jié)果。從表3中可以看出，LS方法的擬合誤差較大，最大值為 0.060 1 m，最小值為 0.027 3 m，檢核點中誤差為 0.048 3 m；WTLS擬合精度相對于LS方法有較大提高，EGM2008+WTLS擬合精度最高，其中檢核點誤差最大值為 0.047 1 m，最小值為 0.005 4 m，中誤差為 0.025 3 m。

5 結(jié) 論

傳統(tǒng)最小二乘方法(LS)只考慮了觀測向量中的誤差，當(dāng)設(shè)計矩陣A中含有誤差時，所得結(jié)果通常是與實際情況有較大偏差的。整體最小二乘雖然考慮了設(shè)計矩陣中的誤差，但是將其與觀測向量誤差等精度處理，忽略了兩種數(shù)據(jù)的不等精度的現(xiàn)實。本文首先采用新一代的超高階地球重力場模型，計算區(qū)域模型高程異常，獲得剩余高程異常后，利用加權(quán)整體最小二乘方法，通過引入EIV模型，并迭代計算各個GPS/水準(zhǔn)點數(shù)據(jù)的權(quán)，獲得更加合理的擬合模型。實際高精度CORS系統(tǒng)GPS/水準(zhǔn)數(shù)據(jù)擬合結(jié)果的外部檢核精度表明，EGM2008+WTLS方法GPS高程擬合模型精度優(yōu)于LS和WTLS方法。

[1] 魏子卿，王剛. 用地球位模型和GPS/水準(zhǔn)數(shù)據(jù)確定我國大陸似大地水準(zhǔn)面[J]. 測繪學(xué)報，2003，32(1)：1～5.

[2] 高原，張恒璟，趙春江. 多項式曲面模型在GPS高程擬合中的應(yīng)用[J]. 測繪科學(xué)，2011，36(3)：179～181.

[3] 王解先. 工業(yè)測量中一種二次曲面的擬合方法[J]. 武漢大學(xué)學(xué)報·信息科學(xué)版，2007，32(1)：47～50.

[4] 吳良才，胡振琪. GPS高程轉(zhuǎn)換方法和正常高計算[J]. 測繪學(xué)院學(xué)報，2004，21(4)：256～258.

[5] 高偉，徐紹銓. GPS高程分區(qū)擬合轉(zhuǎn)換正常高的研究[J]. 武漢大學(xué)學(xué)報·信息科學(xué)版，2004，29(10)：908～911.

[6] 趙輝，張書畢，張秋昭. 基于加權(quán)總體最小二乘法的GPS高程擬合[J]. 大地測量與地球動力學(xué)，2011(5)：88～91.

[7] VAN HUFFEL S，VANDEWALLE.The Total Least Squares and Least Squares Techniq ues in the Presence of Errors on all Data[J]. SIAM Journal on Numerical Anal-ysis，1991，25(5)：765～769.

[8] SCHAFFRIN B，ANDREAS W.On Weighted Total Least squares Adjustment for Linear Regression[J].Journal of Geodesy，2008，82(7)：415～421.

[9] Felus Y A，Schaffrin B.Performinng Similarity Transformations Using the Error-in-variables Model[C]. ASPRS 2005 Annual Conference Baltimore，Maryland，2005.

[10] 丁海勇，孫景領(lǐng). GPS高程轉(zhuǎn)換的總體最小二乘方法研究[J]. 大地測量與地球動力學(xué)，2013，33(3)：52～55.

[11] GOLUB G H，VAN LOAN CH.An Analysis of the Total Least Squares Problem [J]. SIAMJournalonNu-merical Analysis，1980，17(6)：883～893.

[12] 陳俊勇，對SRTM3和GTOPO30地形數(shù)據(jù)質(zhì)量的評估[J]. 武漢大學(xué)學(xué)報·信息科學(xué)版2005，30(11)：941～942.

[13] ?gren J. Evaluation of EGM2008 and PGM2007A over Sweden[J]. Newton’s Bull，2009，4：99～109.

[14] 張興福，劉成. 綜合EGM2008模型和SRTM/DTM2006.0剩余地形模型的GPS高程轉(zhuǎn)換方法[J]. 測繪學(xué)報，2012，41(1)：25～32.

Application of Integrated EGM2008 and Weighted Total Least Squares in GPS Height Fitting

Fan Kunfei1，Yi Guixuan2，Kong Jian3，Liu Li4

(1.Nanning Land Surveying and Mapping Geographic Information Center，Nanning 530021，China；2.Wuhan Geomatics Institute，Wuhan 430000，China；3.Wuhan University，Wuhan 430079，China；4.The first surveying and Mapping Institute of Zhejiang Province，Hangzhou 310012，China)

Based on the calculation of the residual height anomaly by using the earth gravity model EGM2008 and GPS/ leveling data，the weighted total least squares(WTLS)method is employed to carry out the three curve fitting of the residual height anomaly and obtain more accurate GPS height fitting model. The fitting of the CORS network GPS/ level data and the external precision test show that，comprehensive EGM2008 and WTLS methods can significantly improve the accuracy of GPS height anomaly fitting.

EGM2008；WTLS；height anomaly；three degree polynomial

1672-8262(2016)06-88-05

P228

A

2016—09—02

范昆飛(1986—)，男，工程師，主要研究方向為CORS維護(hù)與應(yīng)用。

南寧市第三批特聘專家科研項目；南寧市人才小高地資助項目。 獲獎項目：本論文的研究工作獲得中國地理信息產(chǎn)業(yè)協(xié)會地理信息科技進(jìn)步獎二等獎。