EGM2008和加權整體最小二乘在GPS高程擬合中的應用

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

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

EGM2008和加權整體最小二乘在GPS高程擬合中的應用

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

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

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

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

1 引 言

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

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

2 GPS高程擬合與加權整體最小二乘

2.1 三次曲面函數

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

(1)

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

(2)

設在該區域內有n個GPS水準點，則曲面函數的矩陣形式為：

l=Aβ

(3)

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

矩陣的具體形式為：

(4)

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

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

l-ε=(A-σ)β

(5)

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

QA=Q0?Qx

(6)

WTLS的迭代步驟為：

3 EGM2008地球重力場模型

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

(7)

利用已知GPS/水準點數據可以計算得到水準高程異常ζ，EGM2008可以計算得到GPS/水準點的模型高程異常ζG，則剩余高程異常

ζres=ζ-ζG

(8)

通過三次曲面函數采用加權整體最小二乘方法對剩余高程異常ζres進行擬合，得到區域剩余高程異常模型，與EGM2008計算得到的模型高程異常ζG結合，最終獲得該區域高程異常模型。

4 實例分析

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

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

圖1 GPS/水準點分布

含有粗差的GPS/水準數據 表1

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

曲面擬合參數 表2

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

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

圖2 LS方法GPS高程異常擬合結果

圖3 WTLS方法GPS高程異常擬合結果

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

精度統計 表3

利用LS、WTLS、EGM2008+WTLS方法得到區域高程異常模型后，將9個檢核點的平面坐標帶入模型，獲得這些點的擬合高程異常，與已知的水準高程異常做差得到模型的擬合誤差。表3給出了模型擬合的精度統計結果。從表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 結 論

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

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

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

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

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

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

[6] 趙輝，張書畢，張秋昭. 基于加權總體最小二乘法的GPS高程擬合[J]. 大地測量與地球動力學，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] 丁海勇，孫景領. GPS高程轉換的總體最小二乘方法研究[J]. 大地測量與地球動力學，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地形數據質量的評估[J]. 武漢大學學報·信息科學版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高程轉換方法[J]. 測繪學報，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維護與應用。

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