GOCE＋EGM08＋RTM 組合模型計算高程異常

徐新強

1 西安測繪信息技術(shù)總站，西安市西影路36號，710054

GOCE任務(wù)的基本目標(biāo)是在以前所未有的精度恢復(fù)地球中高階重力場、實現(xiàn)空間分辨率100 km 的情況下，大地水準(zhǔn)面精度達(dá)到1cm。截止2014－07，歐空局已發(fā)布5代GOCE 重力場模型，第5代模型采用了4a的GOCE 衛(wèi)星數(shù)據(jù)，包括時域法模型GO＿CONS＿GCF＿2＿TIM＿R5（簡稱Tim5）和直接法模型GO＿CONS＿GCF＿2＿DIR＿R5（簡稱Dir5）兩個產(chǎn)品。Tim5只采用了GOCE衛(wèi)星數(shù)據(jù)，即重力梯度數(shù)據(jù)和GPS高低衛(wèi)衛(wèi)跟蹤軌道數(shù)據(jù)；而Dir5除采用GOCE 衛(wèi)星數(shù)據(jù)外，還采用了25a的LAGEOS 1／2 激光測衛(wèi)數(shù)據(jù)、10a的GRACE低低跟蹤衛(wèi)星數(shù)據(jù)和K 波段距離變率數(shù)據(jù)。

在缺乏重力數(shù)據(jù)的地區(qū)，通常用超高階重力模型（如EGM08）來計算重力高程異常，然后擬合到已知GPS／水準(zhǔn)點上。EGM08 模型構(gòu)建過程中，對重力數(shù)據(jù)專利區(qū)或無數(shù)據(jù)區(qū)，采用如下方法進(jìn)行填充［1］：用GRACE 衛(wèi)星跟蹤數(shù)據(jù)建立GGM02S模型的前60 階，加上EGM96 模型的61～360階球諧系數(shù)，再加上剩余地形重力模型361～2 160階系數(shù)，計算得到5′×5′的重力異常。而EGM96模型的重力異常又是根據(jù)已有30′重力異常和5′GEOSAT 海洋大地水準(zhǔn)面，采用Forsberg協(xié)方差模型的最小二乘配置法內(nèi)插得到。盡管中國大陸屬于數(shù)據(jù)專利區(qū)，但其EGM08模型高程異常整體精度與全球精度相當(dāng)，為20cm［2］。

從EGM08模型的構(gòu)建過程來看，在缺乏重力數(shù)據(jù)的地區(qū)，至少可以從兩方面提高模型高程異常的計算精度：1）由高程異常的能量譜可知，高程異常主要由重力場模型的中低頻信號決定，因此可以用GOCE梯度信息來提高EGM08模型中低頻段的精度；2）EGM08模型完全階次至2 159，相應(yīng)的空間分辨率為5′×5′，舍去分辨率高于5′×5′的甚高頻信號，因此可以借助于高分辨率的地形數(shù)據(jù)來彌補這一損失。關(guān)于第一方面，文獻(xiàn)［3］用實測GPS／水準(zhǔn)數(shù)據(jù)證明，GOCE 模型相對于EGM08在中頻更有優(yōu)勢，尤其是像西藏這樣的地形復(fù)雜地區(qū)。關(guān)于第二方面，文獻(xiàn)［4－5］用SRTM 數(shù)字地形模型減去DTM2006.0模型計算的地形高，構(gòu)造剩余地形模型RTM，并用RTM高程異常彌補EGM08模型的甚高頻信號。文獻(xiàn)［4］在阿爾卑斯山高山區(qū)的實驗表明，EGM08＋RTM 高程異常精度比單獨使用EGM08時提高了2cm，相對精度提高了49%，同樣，文獻(xiàn)［5］中組合模型高程異常精度提高了1.3cm，相對精度提高了42%，說明EGM08＋RTM 組合模型能大幅提高山區(qū)垂線偏差精度。文獻(xiàn)［6］在阿爾卑斯山的垂線偏差實驗表明，重力場高頻分量對垂線偏差影響很大，RTM 對垂線偏差精度的貢獻(xiàn)非常明顯，子午分量從3.5″提高到0.8″，卯酉分量從3.2″提高到0.7″。因此，本文嘗試采用GOCE＋EMG08＋RTM 的組合模型來計算重力高程異常。

1 基本原理

1.1 模型高程異常

重力場模型高程異常的計算詳見文獻(xiàn)［7］，軟件可采用美國地理空間情報局NGA 提供的hsynth＿WGS84.f。不同組織發(fā)布的重力場模型對應(yīng)的參考橢球長半軸、地球引力常量和GRS80橢球都不一致，因此，首先需要將重力場位系數(shù)統(tǒng)一到GRS80橢球上來：

1.2 RTM 高程異常

DTM2006.0是隨EGM08一起發(fā)布的5′分辨率的球諧系數(shù)模型，在EGM08構(gòu)建中使用該模型。DTM2006.0是將3″×3″的SRTM 數(shù)字地形模型平滑成5′×5′的格網(wǎng)，再調(diào)和生成球諧系數(shù)，因此DTM2006.0相當(dāng)于低通濾波器，濾掉了地形球諧模型中高于2 160階的甚高頻信號。由SRTM－DTM2006.0生成的RTM 引力位可表示比EGM08更高的甚高頻信息。采用下式計算模型地形高：

每個六棱柱RTM 在平面近似下對計算點的引力位為［8］：

式中，G為引力常數(shù)，ρ0為標(biāo)準(zhǔn)地形質(zhì)量密度，r為積分流動點（x，y，z）到計算點的距離，x1、x2、y1、y2、z1、z2為六棱柱邊界角點坐標(biāo)，z2－z1表示剩余高程。對計算點周圍所有的六棱柱引力位求和，可得到整個剩余地形的引力位，然后應(yīng)用Bruns公式得到由RTM 引起的計算點高程異常：

式中，γQ為計算所對應(yīng)的近似地形面上的正常重力。RTM 高程異常的計算可采用Forsberg提供的TC.f軟件［9］。

2 數(shù)值實驗

2.1 重力場模型的選擇

按文獻(xiàn)［7］給出的重力場模型功率譜分析，繪制4個GOCE模型和EGM08模型的大地水準(zhǔn)面累積誤差圖（圖1）。從圖1 可看出，隨著GOCE重力梯度數(shù)據(jù)的不斷積累，解算的第5代模型精度明顯比第4代高，且在中波段的精度GOCE 重力場模型都高于EGM08模型。在200 階處，即相應(yīng)空間分辨率為100km 時，Dir5累積誤差為0.8cm，而Tim5累積誤差為2.2cm。

圖1 不同模型的累積大地水準(zhǔn)面誤差Fig.1 Cumulative geoid errors for different geopotential models

Dir系列模型精度高于同代Tim 的原因是：1）Dir系列模型除了使用GOCE 數(shù)據(jù)外，還使用了LAGEOS激光測衛(wèi)數(shù)據(jù)和GRACE衛(wèi)星數(shù)據(jù)，也就是說Dir系列模型充分利用了3種衛(wèi)星重力測量技術(shù)在不同頻段恢復(fù)重力場的優(yōu)勢。衛(wèi)星激光測距主要用于確定地球動力學(xué)扁率J2項，衛(wèi)星跟蹤衛(wèi)星技術(shù)主要用于恢復(fù)70階以下的重力場，衛(wèi)星重力梯度技術(shù)在恢復(fù)70階以上的地球重力場中更有優(yōu)勢。2）Dir模型構(gòu)造過程中不依賴任何先驗地球重力場模型，理論框架非常嚴(yán)密，因此反演地球重力場精度較高［10］。

目前GOCE 模型精度在200 階之后還不太高，累積大地水準(zhǔn)面誤差不太穩(wěn)定，而EGM08模型在200階之后的累積大地水準(zhǔn)面誤差逐漸趨于穩(wěn)定?？紤]到Dir5 模型和EGM08 模型在不同頻段的優(yōu)勢，本文采用200階Dir5模型加上201～2 160階的EGM08模型來計算重力高程異常。在北緯27°～29°、東經(jīng)86°～88°的珠峰地區(qū)，計算200階次的Dir5 和同階次EGM08 高程異常差（圖2）。從圖2中看出，兩個模型在該地區(qū)差異非常明顯，最大達(dá)到1.4 m，說明Dir5在確定高程異常的長波項中具有很大潛力。

圖2 Dir5和EGM08模型在珠峰地區(qū)的高程異常差Fig.2 Differences of height anomalies between Dir5and EGM08geopotential models in Mount Everest area

2.2 RTM 高程異常

地形效應(yīng)主要與地形起伏有關(guān)，與地形絕對海拔高度無關(guān)。在珠峰附近選取一個1°×1°區(qū)域的SRTM 模型，該地區(qū)地形落差較大，最低高程193m，最高高程8 771m，平均高程2 755m，標(biāo)準(zhǔn)差1 607 m。按式（3）計算的RTM 高程見圖3，RTM 高程最低－1 430m，最高2 706m，平均11m，標(biāo)準(zhǔn)差434m。

計算RTM 引力位非常耗時，在實際計算中離計算點較近的區(qū)域通常采用高分辨率的剩余地形數(shù)據(jù)，離計算點較遠(yuǎn)的區(qū)域則采用粗網(wǎng)格數(shù)據(jù)。按文獻(xiàn)［4－5］的建議，本文內(nèi)區(qū)精細(xì)網(wǎng)格積分半徑選40km，外區(qū)粗網(wǎng)格積分半徑為160km。實際上TC.f軟件就是這么處理的，它需要準(zhǔn)備3 個地形文件作為輸入，本文準(zhǔn)備的文件是：1）3″的SRTM 精細(xì)網(wǎng)格文件；2）15″的粗網(wǎng)格文件，由3″的SRTM 模型直接平滑獲得；3）3″的參考面地形模型，由DTM2006.0模型按式（2）計算得到。

在實驗區(qū)域按7.5″的間隔逐點計算480×480個網(wǎng)格點的RTM 高程異常（圖4），最小值－6.5cm，最大值9.7cm，平均值1.2cm，標(biāo)準(zhǔn)差2.3cm，且RTM 高程異常輪廓與圖3基本一致。從圖4看出，高精度大地水準(zhǔn)面計算中，RTM 對高程異常的影響不容忽視。

圖3 剩余地形高（SRTM－DTM2006.0）Fig.3 Residual terrain elevations（SRTM－DTM2006.0）

圖4 剩余地形高程異常Fig.4 Height anomalies about residual terrain model

2.3 GPS／水準(zhǔn)數(shù)據(jù)比較

在藏南地區(qū)收集228 個C 級GPS／水準(zhǔn)數(shù)據(jù)，數(shù)據(jù)分布呈東西走向，東西寬800km，南北寬200km，平均海拔約4 500 m，GPS／水準(zhǔn)主要沿道路兩側(cè)布設(shè)，地形起伏不大。用GPS／水準(zhǔn)點上的幾何高程異常，減去模型計算的高程異常，統(tǒng)計結(jié)果見表1。

表1 GPS／水準(zhǔn)與組合模型高程異常之差的統(tǒng)計表Tab.1 Statistics of the height anomalies differences between GPS／leveling and different combined geopotential models

從表1 可看出，EGM08 模型在藏南精度較差，使用GOCE 重力場模型替換EGM08中低頻段的組合模型Dir5＋EGM08的精度提高了0.26 m，相對精度提升43%，而Dir5＋EGM08＋RTM組合模型精度提高了0.27 m，相對精度提升45%，即是說在該地區(qū)RTM 只貢獻(xiàn)了1cm 的精度。這是因為高程異常主要受重力場中低頻支配，且該地區(qū)雖然整體海拔高，但地形起伏不大，RTM 影響不明顯。

3 結(jié) 語

1）圖1表明，截至2014－07，綜合SLR、GRACE和GOCE數(shù)據(jù)的Dir5模型代表了GOCE 重力場發(fā)展的最高水平。

2）圖2和表1表明，在缺少重力數(shù)據(jù)的地區(qū)，尤其是高海拔地區(qū)，GOCE 大地水準(zhǔn)面展現(xiàn)了明顯的優(yōu)勢。

3）圖4表明，在地形落差較大的地區(qū)，RTM能有效地彌補重力場的甚高頻信號。圖中采用3″×3″分辨率的RTM 網(wǎng)格數(shù)據(jù)，某種程度上相當(dāng)于將重力場擴展到了216 000 階，因此RTM 導(dǎo)致的高程異常不容忽視。

4）表1表明，GOCE＋EGM08＋RTM 組合模型從兩端（中低頻和甚高頻）改造EGM08重力場模型，提高了無重力數(shù)據(jù)地區(qū)高程異常計算的精度。

［1］楊金玉，張訓(xùn)華，張菲菲，等.EGM2008地球重力模型數(shù)據(jù)在中國大陸地區(qū)的精度分析［J］.地球物理學(xué)進(jìn)展，2012，27（4）：1 298－1 306（Yang Jinyu，Zhang Xunhua，Zhang Feifei.On the Accuracy of EGM2008Earth Gravitational Model in Chinese Mainland［J］.Progress in Geophysics，2012，27（4）：1 298－1 306）

［2］章傳銀，郭春喜，陳俊勇，等.EGM2008 地球重力場模型在中國大陸適用性分析［J］.測繪學(xué)報，2009，38（4）：283－289（Zhang Chuanyin，Guo Chunxi，Chen Junyong.EGM 2008and Its Application Analysis in Chinese Mainland［J］.Acta Geodaetica et Cartographica Sinica，2009，38（4）：283－289）

［3］高春春，陸洋，張子占，等.GOCE 地球重力場模型在青藏地區(qū)的評價及分析［J］.測繪學(xué)報，2014，43（1）：21－29（Gao Chunchun，Lu Yang，Zhang Zizhan.Evaluation and Analysis of GOCE Earth Gravity Field Model in Qinghai－Tibet Region［J］.Acta Geodaetica et Cartographica Sinica，2014，43（1）：21－29）

［4］Hirt C，F(xiàn)eatherstone W E，Marti U.Combining EGM2008 and SRTM／DTM2006.0 Residual Terrain Model Data to Improve Quasigeoid Computations in Mountainous Areas Devoid of Gravity Data［J］.Journal of Geodesy，2010，84（9）：557－567

［5］張興福，劉成.綜合EGM2008模型和SRTM／DTM2006.0剩余地形模型的GPS高程轉(zhuǎn)換方法［J］.測繪學(xué)報，2012，41（1）：25－32（Zhang Xingfu，Liu Cheng.The Approach of GPS Height Transformation Based on EGM2008 and SRTM／DTM2006.0 Residual Terrain Model［J］.Acta Geodaetica et Cartographica Sinica，2012，41（1）：25－32）

［6］Hirt C.Prediction of Vertical Deflections from High－Degree Spherical Harmonic Synthesis and Residual Terrain Model Data［J］.Journal of Geodesy，2010，84（3）：179－190

［7］王正濤，黨亞民，晁定波.超高階地球重力位模型確定的理論與方法［M］.北京：測繪出版社，2011（Wang Zhengtao，Dang Yamin，Chao Dingbo.Theory and Methodology of Ultra－High－Degree Geopotential Model Determination［M］.Beijing：Surveying and Mapping Press，2011）

［8］Nagy D，Papp G，Benedek J.The Gravitational Potential and Its Derivatives for the Prism［J］.Journal of Geodesy，2000，74（7－8）：552－560

［9］Forsberg R.A Study of Terrain Reductions，Density Anomalies and Geophysical Inversion Methods in Gravity Field Modeling［R］.Report 355，Department of Geodetic Science and Surveying，Columbus，1984

［10］鄭偉，許厚澤，鐘敏，等.衛(wèi)星重力梯度反演研究進(jìn)展［J］.大地測量與地球動力學(xué)，2014，34（4）：1－8（Zheng Wei，Xu Houze，Zhong Min.Research Progress in Satellite Gravity Gradiemetry Recovery［J］.Journal of Geodesy and Geodynamics，2014，34（4）：1－8）