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

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

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

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

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

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

EGM2008；WTLS；高程異常；三次多項(xiàng)式

1 引 言

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

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

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

2.1 三次曲面函數(shù)

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

(1)

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

(2)

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

l=Aβ

(3)

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

矩陣的具體形式為：

(4)

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

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

l-ε=(A-σ)β

(5)

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

QA=Q0?Qx

(6)

WTLS的迭代步驟為：

3 EGM2008地球重力場(chǎng)模型

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

(7)

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

ζres=ζ-ζG

(8)

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

4 實(shí)例分析

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

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

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

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

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

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

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

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

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

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

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

精度統(tǒng)計(jì) 表3

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

5 結(jié) 論

傳統(tǒng)最小二乘方法(LS)只考慮了觀測(cè)向量中的誤差，當(dāng)設(shè)計(jì)矩陣A中含有誤差時(shí)，所得結(jié)果通常是與實(shí)際情況有較大偏差的。整體最小二乘雖然考慮了設(shè)計(jì)矩陣中的誤差，但是將其與觀測(cè)向量誤差等精度處理，忽略了兩種數(shù)據(jù)的不等精度的現(xiàn)實(shí)。本文首先采用新一代的超高階地球重力場(chǎng)模型，計(jì)算區(qū)域模型高程異常，獲得剩余高程異常后，利用加權(quán)整體最小二乘方法，通過引入EIV模型，并迭代計(jì)算各個(gè)GPS/水準(zhǔn)點(diǎn)數(shù)據(jù)的權(quán)，獲得更加合理的擬合模型。實(shí)際高精度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]. 測(cè)繪學(xué)報(bào)，2003，32(1)：1～5.

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

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

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

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

[6] 趙輝，張書畢，張秋昭. 基于加權(quán)總體最小二乘法的GPS高程擬合[J]. 大地測(cè)量與地球動(dòng)力學(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]. 大地測(cè)量與地球動(dòng)力學(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] 陳俊勇，對(duì)SRTM3和GTOPO30地形數(shù)據(jù)質(zhì)量的評(píng)估[J]. 武漢大學(xué)學(xué)報(bào)·信息科學(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]. 測(cè)繪學(xué)報(bào)，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—)，男，工程師，主要研究方向?yàn)镃ORS維護(hù)與應(yīng)用。

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