基于EGM2008地球重力場模型-Griddata函數組合模型的似大地水準面精化探討

付來迎，周世健，許亞男，王奉偉，周 青

(1.東華理工大學測繪工程學院，330013，南昌；2.江西省科學院，330096，南昌)

基于EGM2008地球重力場模型-Griddata函數組合模型的似大地水準面精化探討

付來迎1，周世健2，許亞男1，王奉偉1，周 青1

(1.東華理工大學測繪工程學院，330013，南昌；2.江西省科學院，330096，南昌)

介紹了EGM2008地球重力場模型以及MATLAB提供的Griddata函數，結合工程實例，探討基于EGM2008地球重力場模型-Griddata函數組合模型對似大地水準面進行精化，在Griddata函數提供的4種不同插值方法下進行分析比較，并根據4種不同插值的特點進行改進，最后得出基于EGM2008地球重力場模型-Griddata函數法中利用四格點樣條插值法替代3次多項式插值法中凸包點法的擬合效果最好，具有很好的實際推廣意義。

EGM2008；Griddata函數；移去-恢復法；似大地水準面精化

0 引言

由GPS測量定位得到的基線向量，經平差后可得到高精度的大地高程。若網中有一點或多點具有精確的WGS-84大地坐標系的大地高程，則在GPS網平差后，可得各個GPS點的WGS-84大地坐標系的大地高程。GPS相對定位高程方面的相對精度一般可達(2×3)×10-6，在絕對精度方面，對于10 km以下的基線邊長，可達幾個厘米，如果在觀測和計算時采用一些消除誤差的措施，其精度將優于±1 cm。但在實際應用中，地面點一般采用正常高程系統[1]。因此，為了盡可能減少外業工作量、廣泛應用新技術，充分發揮GPS測量方便、快捷、精度高、成本低等優點，就必須知道WGS-84參考橢球面與似大地水準面之間的差距，即高程異常。知道了每個GPS點的高程異常ζ，就能通過式(1)，計算出相應的正常高，實現用GPS測量代替水準測量。

H大地高=H正常高+ζ

(1)

確定高程異常方法大致可分為幾何曲面擬合法和重力法兩大類。幾何曲面擬合法就是根據區域內若干個公共點上的高程異常值，構造某一種曲面逼近似大地水準面，由于所構造的曲面不同，計算方法也不同，其主要方法有：平面擬合法、曲面擬合法、多面函數擬合法、樣條函數法等[2]。重力法是指利用局部或者全球重力數據求解高程異常，將GPS大地高轉換為正常高。重力法充分利用了地球物理場，但由于重力數據缺乏和重力場模型精度對于工程應用而言較低，所以工程建設一般都不采用重力法直接進行GPS高程轉換[3]。在區域似大地水準面精化中，最嚴密、有效的方法是利用由地球重力場模型、地面重力數據和DEM數據獲得的該區域剩余重力異常，采用移去-恢復法確定重力似大地水準面，再用GPS水準數據對重力似大地水準面進行擬合，求得與國家或地方高程系統定義一致的似大地水準面[4]。由于EGM2008模型的精度明顯好于其他模型[5]，而且它的相對精度比絕對精度高[6]，所以將應用EGM2008模型計算的高程異常差和GPS水準數據確定局部似大地水準面的方法，并通過算例對其應用的可行性進行分析。

1 EGM2008全球重力場模型簡介

EGM2008是美國國家地理空間情報局(National Geospatia1-Intelligence Agency)，簡稱NGA，經過多年的研究和總結，在以往構建地球重力場模型的經驗和理論基礎上，采用最先進的建模技術與算法，以PGM2007B(PGM2007A的變種模型)為參考模型，利用GRACE衛星采集的重力數據和全球5′×5′的重力異常數據，TOPEX衛星測高數據以及現勢性好、分辨率高的地形數據，結合精度高，覆蓋面廣的地面重力數據完成的最新一代全球重力模型[7]。該模型的階次完全至2 159(另外球諧系數的階擴展至2 190次)，相當于模型的空間分辨率約為5′(約9 km)。地球重力場球諧函數模型EGM2008由一系列完全正常化的球諧函數系數組成，根據球諧函數連續疊加可計算擾動位T，利用Bruns公式可求得地面上任意點的高程異常[8]。

(2)

2 EGM2008-Griddata組合模型原理

由于EGM2008模型定義的全球高程基準與我國采用的國家高程基準不同，兩者之間存在著一個系統差，因而，首先將高程異常可如式(3)可分成兩部分

ζ=ζGM+△ζ

(3)

式中：ζGM為EGM2008地球重力場模型求得的高程異常；△ζ為實測高程異常與EGM2008地球重力場高程異常的殘差。從而，當利用EGM2008模型計算高程異常進而推算正常高時，式(1)相應改為：

H大地高=H正常高+ζGM+△ζ

(4)

利用“移去-恢復”法，首先，將基于已有的m個GPS水準聯測點，可通過式(1)的變形式：ζi=H大地高i-H正常高i,(i=1，2，3，…，k)，計算出這k個點的高程異常值ζi；再利用地球重力場模型EGM2008計算這些點的高程異常的近似值ζGMi，然后根據式(3)計算出這k個點的高程異常殘差△ζi。第2，將求出的k個高程異常殘差值△ζi作為已知數據，用MATLAB所提供的Griddata函數插值方法進行擬合計算，最后再內插出未知點的高程異常殘差△ζ。最后，利用地球重力場模型EGM2008求出未聯測水準的點上的高程異常近似值ζGM，再與擬合出的該點的高程異常殘差值△ζ相加，得出最終的高程異常值ζ，再通過利用GPS測得的大地高數據由式(1)計算出所有待求點的正常高。

3 實例分析

MATLAB是MathWorks公司推出的一套高性能的數值計算和可視化軟件，集數值分析、矩陣運算與圖形顯示于一體，可方便地應用于數學計算、數據分析和工程繪圖，所采用的數值計算方法均采用公認的先進、可靠的算法，其程序均由世界一流專家編制并經高度優化[9]。MATLAB中的Griddata函數可以將位于同一空間坐標系下的散點插值為規則格網，提供了4種插值方法：線性插值(linear)、3次多項式插值(cubic)、最近點插值(nearest)、四格點樣條插值(v4)，可以方便地實現結合鄰近離散點分布特征的光滑曲面擬合。其中線性插值和最近點插值結果構成的曲面不光滑不連續，3次多項式插值和四格點樣條插值結果構成的曲面較光滑[10]。

為對Griddata函數提供的4種插值方法進行合理的選擇，本文以某地區16個GPS水準數據為例，該測區南北跨度約20 km，東西跨度約15 km，地形較為復雜，測區起伏較大，最大高差約為170 m之多，其點位分布情況如圖1所示。

圖1 點位分布情況

圖2 4種方法4個學習樣本擬合殘差對照圖

圖3 4種方法6個學習樣本擬合殘差對照圖

圖4 4種方法8個學習樣本擬合殘差對照圖

通過以上3個圖對比可以看出：較于最近點插值法替代，四格點樣條插值法替代線性插值法和3次多項式插值法中出現的“凸包”點的效果更好。同時，從以上3幅圖可以看出“凸包”點替代后的3次多項式插值法插值殘差值最小，替代后的線性插值法效果次之，四格點樣條插值法效果稍差，最近點插值法殘差波動性較大且殘差值較大；并且4種方法都會隨著學習樣本數量的增加，插值的殘差減小且波動性會降低。

通過4種方法不同數量學習樣本得到的擬合成果還可以分別通過內符合精度和外符合精度進行評定，其中內符合精度可按式(5)計算求得，外符合精度可按式(6)求得

(5)

(6)

其中,vi為各點的擬合殘差，m為學習集樣本個數，n為總體樣本數量。其結果如表1所示。

由表1可以看出，隨著樣本數量的增多，4種不同的插值方法插值效果都有所改善，最大偏差值都在減小，外符合精度不斷提高；同時，在線性插值法與3次多項式插值法插值中利用四格點樣條插值法進行凸包點替代的外符合精度高于最近點插值法進行替代的外符合精度。

表1 不同插值方法不同數量學習樣本精度統計表

4 總結

通過工程實例采用基于EGM2008重力場模型的Griddata函數的4種不同方法進行插值比較，利用四格點樣條插值法進行替代的線性插值法與3次多項式插值法有著較好的插值效果，基于插值效果與擬和曲面光滑性的綜合考慮，利用四格點樣條插值法進行凸包點替代的3次多項式插值法在采用基于EGM2008重力場模型的Griddata函數模型的似大地水準面擬合中有更好的適用性，并在實際生產中具有較好的實際推廣意義。

[1] 張福榮,田倩.GPS測量技術與應用[M].成都:西南交通大學出版社,2013:177-178.

[2]張勤,李家權,等.GPS測量原理及應用[M].北京:科學出版社,2005:214-219.

[3]谷延超,范東明.顧及EGM2008和殘差地形模型的GPS高程轉換方法研究[J].測繪工程,2013,22(2):26-29.

[4]馮林剛,郅軍義,寶因烏力吉.應用EGM2008模型和GPS/水準數據確定局部似大地水準面[J].測繪通報,2011(1):18-20.

[5]章傳銀,郭春喜,陳俊勇,等.EGM2008地球重力場模型在中國大陸適用性分析[J].測繪學報,2009,38(4):283-289.

[6]張興福,劉成,劉紅新.利用GPS冰準數據檢核EGM2008重力場模型的精度[J].測繪通報,2009(2):7-9.

[7]Pavlis N K,Holmes S A,Kenyon S C,etal.An Earth Gravitational Model to Degree 2160:EGM2008[C].Vienna:Presented at the 2008 General Assembly of the European Geosciences Union,2008.

[8]游為,范東明,付淑娟,等.GPS高程轉換的新方法研究[J].工程勘察,2009,37(3):60-62.

[9]求是科技.Matlab7.0從入門到精通[M].北京:人民郵電出版社,2006.

[10]吳守亮.基于Matlab的三維數字地形模擬及空間分析[J].測繪工程,2011,20(3):54-57.

DiscussionofRefiningtheGeoidBasedonEGM2008-GriddataFunctionCombinationModel

FU Laiying1,ZHOU Shijian2,XU Yanan1,WANG Fengwei1,ZHOU Qing1

(1.East China Institute of Technology.Faculty of Geomatics,330013,Nanchang,PRC;2.Jiangxi Academy of Science,330096,Nanchang,PRC)

The paper introduced the EGM2008 gravity field model and the Griddata function that was provided by MATLAB.At the same time,the paper combining with the engineering example presented a new method which was the EGM2008 gravity field model-Griddata function model to refine the geoid.Meanwhile,Analysis and comparison was made with four kinds of interpolation methods provided by Griddata function.At the last,we got a conclusion that the new method had a good fitting effect on the geoid refinement by analyzing the result of the date,and it had an important significance in the actual production.

EGM2008;Griddata function;remove-restore method;refine the geoid

2014-10-12;

2014-11-21

付來迎(1990-)，男，山東萊蕪人，碩士研究生，主要研究方向為現代測量數據處理。

國家自然科學基金資助項目(41374007)；2014年江西省研究生創新專項資金項目(YC2014-S317)。

10.13990/j.issn1001-3679.2014.06.012

P223

A

1001-3679(2014)06-0794-04