付來(lái)迎，周世健，許亞男，王奉偉，周 青

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

基于EGM2008地球重力場(chǎng)模型-Griddata函數(shù)組合模型的似大地水準(zhǔn)面精化探討

付來(lái)迎1，周世健2，許亞男1，王奉偉1，周 青1

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

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

EGM2008；Griddata函數(shù)；移去-恢復(fù)法；似大地水準(zhǔn)面精化

0 引言

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

H大地高=H正常高+ζ

(1)

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

1 EGM2008全球重力場(chǎng)模型簡(jiǎn)介

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

(2)

2 EGM2008-Griddata組合模型原理

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

ζ=ζGM+△ζ

(3)

式中：ζGM為EGM2008地球重力場(chǎng)模型求得的高程異常；△ζ為實(shí)測(cè)高程異常與EGM2008地球重力場(chǎng)高程異常的殘差。從而，當(dāng)利用EGM2008模型計(jì)算高程異常進(jìn)而推算正常高時(shí)，式(1)相應(yīng)改為：

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

(4)

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

3 實(shí)例分析

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

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

圖1 點(diǎn)位分布情況

圖2 4種方法4個(gè)學(xué)習(xí)樣本擬合殘差對(duì)照?qǐng)D

圖3 4種方法6個(gè)學(xué)習(xí)樣本擬合殘差對(duì)照?qǐng)D

圖4 4種方法8個(gè)學(xué)習(xí)樣本擬合殘差對(duì)照?qǐng)D

通過(guò)以上3個(gè)圖對(duì)比可以看出：較于最近點(diǎn)插值法替代，四格點(diǎn)樣條插值法替代線性插值法和3次多項(xiàng)式插值法中出現(xiàn)的“凸包”點(diǎn)的效果更好。同時(shí)，從以上3幅圖可以看出“凸包”點(diǎn)替代后的3次多項(xiàng)式插值法插值殘差值最小，替代后的線性插值法效果次之，四格點(diǎn)樣條插值法效果稍差，最近點(diǎn)插值法殘差波動(dòng)性較大且殘差值較大；并且4種方法都會(huì)隨著學(xué)習(xí)樣本數(shù)量的增加，插值的殘差減小且波動(dòng)性會(huì)降低。

通過(guò)4種方法不同數(shù)量學(xué)習(xí)樣本得到的擬合成果還可以分別通過(guò)內(nèi)符合精度和外符合精度進(jìn)行評(píng)定，其中內(nèi)符合精度可按式(5)計(jì)算求得，外符合精度可按式(6)求得

(5)

(6)

其中,vi為各點(diǎn)的擬合殘差，m為學(xué)習(xí)集樣本個(gè)數(shù)，n為總體樣本數(shù)量。其結(jié)果如表1所示。

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

表1 不同插值方法不同數(shù)量學(xué)習(xí)樣本精度統(tǒng)計(jì)表

4 總結(jié)

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

[1] 張福榮,田倩.GPS測(cè)量技術(shù)與應(yīng)用[M].成都:西南交通大學(xué)出版社,2013:177-178.

[2]張勤,李家權(quán),等.GPS測(cè)量原理及應(yīng)用[M].北京:科學(xué)出版社,2005:214-219.

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

[4]馮林剛,郅軍義,寶因?yàn)趿?應(yīng)用EGM2008模型和GPS/水準(zhǔn)數(shù)據(jù)確定局部似大地水準(zhǔn)面[J].測(cè)繪通報(bào),2011(1):18-20.

[5]章傳銀,郭春喜,陳俊勇,等.EGM2008地球重力場(chǎng)模型在中國(guó)大陸適用性分析[J].測(cè)繪學(xué)報(bào),2009,38(4):283-289.

[6]張興福,劉成,劉紅新.利用GPS冰準(zhǔn)數(shù)據(jù)檢核EGM2008重力場(chǎng)模型的精度[J].測(cè)繪通報(bào),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高程轉(zhuǎn)換的新方法研究[J].工程勘察,2009,37(3):60-62.

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

[10]吳守亮.基于Matlab的三維數(shù)字地形模擬及空間分析[J].測(cè)繪工程,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

付來(lái)迎(1990-)，男，山東萊蕪人，碩士研究生，主要研究方向?yàn)楝F(xiàn)代測(cè)量數(shù)據(jù)處理。

國(guó)家自然科學(xué)基金資助項(xiàng)目(41374007)；2014年江西省研究生創(chuàng)新專項(xiàng)資金項(xiàng)目(YC2014-S317)。

10.13990/j.issn1001-3679.2014.06.012

P223

A

1001-3679(2014)06-0794-04