綜合EGM2008模型和SRTM／DTM2006.0剩余地形模型的GPS高程轉換方法

張興福，劉 成

1.廣東工業大學測繪工程系，廣東廣州510006；2.鐵道第三勘察設計院集團有限公司測繪分院，天津300251

1 引 言

隨著衛星重力計劃的逐步實施和完善，利用衛星重力資料并結合地面重力數據獲得的全球重力場模型的精度和分辨率在逐步提高，特別是EGM2008模型的公布［1］，使得利用重力場模型進行GPS高程轉換的精度和可靠性獲得了大幅度提升。EGM2008模型的階次完全至2159（另外球諧系數的階擴展至2190，次為2159），相當于模型的空間分辨率為5′（約9km）。EGM2008模型采用了GRACE衛星跟蹤數據（ITG－GRACE03S位系數信息以及相應的協方差信息）、衛星測高數據和地面重力數據等，該模型無論在精度還是在分辨率方面均取得了巨大進步［2－3］。文獻［4］利用我國大陸GPS／水準實測高程異常對EGM2008地球重力場模型進行了外部精度測試，結果表明：EGM2008模型高程異常精度在我國大陸的整體精度為20cm，且與全球精度相當。

地球重力場模型展開到一定的階次，對應于一定的空間分辨率，EGM2008模型的階數高達2190，也只有約9km的空間分辨率，很難表示更高頻段的高程異常。為了提高GPS高程轉換的精度，須借助于高分辨率的地形模型。文獻［5］推導出一套適合厘米級似大地水準面精化的Molodensky解實用簡便的算法公式，并結合實例測試驗證了厘米級地形影響的技術特征，文獻［6—8］對RTM技術用于補償EGM2008信號截斷誤差進行了研究，獲得了良好的效果。筆者基于剩余地形模型理論［9－11］，利用SRTM以及DTM2006.0全球地形模型構建RTM數據，并計算剩余地形模型高程異常。任意GPS點的總高程異常可表示為地球重力場模型高程異常、剩余地形模型高程異常以及殘余高程異常3部分總和的形式，通過計算可獲得所有GPS點總高程異常中的前兩項，然后選擇少量GPS／水準點的實測高程異常，扣除EGM2008模型以及SRTM與DTM2006.0模型求得的剩余模型高程異常，對殘余高程異常進行擬合，從而提高GPS高程轉換的精度。利用2個不同測區的高精度GPS／水準數據進行實例分析，結果表明，該方法能顯著提高GPS高程轉換精度。

2 原理與方法

2.1 重力場模型高程異常計算

根據Bruns公式，地球表面上任意點P的模型高程異常可由下式獲得［4］

2.2 剩余地形模型（RTM）高程異常計算

RTM數據表示真實地形表面和參考地形表面的差值，見圖1。一般情況下，真實地形表面用分辨率較高的DTM數據表示，而參考地形表面用分辨率較低的DTM數據表示，該數據是對表示真實地形表面的DTM經過一定平滑處理后獲得的。在地形起伏較大測區，僅僅依靠一定階次的地球重力場模型以及少量的GPS／水準數據很難對短波長的重力場信號進行模型化，而這種方法又是非常經濟實用的GPS高程轉換方法，理論上RTM技術有能力表示由局部地形變化引起的重力場短波長部分。

圖1 剩余地形模型（RTM）Fig.1 Residual terrain model

RTM高程異常計算的方法很多，本文探討利用棱柱積分法將RTM數據轉換到剩余高程異常，每個棱柱（格網）對應的引力位為［13－14］

式中，G為引力常數；r為坐標原點到點（x，y，z）的距離；x1、x2、y1、y2、z1、z2為棱柱體的邊界角點坐標，z2－z1表示剩余高程；ρ0為標準的地形質量密度，取ρ0＝2670kg／m3；式（2）是平面近似，計算中需考慮地球曲率的影響［11］；其他變量的含義及具體計算過程請參考相關文獻［11，13－14］，將棱柱引力位轉換為對應的高程異常，公式為

式中，γQ表示正常重力。GPS點P所對應的RTM高程異常可由計算區域內所有單個棱柱對應的高程異常總合表示，公式為

為了提高棱柱積分的計算速度，在實際應用中，離計算點稍遠的區域采用分辨率稍低DTM格網數據，離計算點較近的區域采用分辨率較高的DTM格網數據，這種實用的剩余地形計算的過程稱為粗糙／詳細DTM格網系統。在計算過程中，一般會給兩個積分計算半徑R1和R2，當積分半徑小于R1時用詳細格網數據，當積分半徑在R1和R2時用粗糙格網數據。另外在計算點P周圍3×3格網的地形數據用雙三次樣條函數進行加密處理，以獲得更光滑的地形數據，這種加密也可避免計算點P處于格網邊界時對計算結果產生影響［11］，粗糙／詳細DTM格網數據使用情況見圖2。

圖2 粗糙／詳細DTM格網使用情況圖Fig.2 Use of coarse and detailed DTM grid

2.3 剩余地形模型（RTM）數據構建方法

利用RTM方法計算高程異常的關鍵是構建RTM數據。目前可用于RTM數據構建的全球DTM數據有3″的SRTM、1″的ASTER GDEM以及局部高精度的DTM數據等。在我國要獲得局部高精度的DTM數據難度較大，而SRTM和GDEM數據是公開的，SRTM（V4.1）的高程精度要優于GDEM。本文RTM的構建方法為：①采用免費的3″SRTM（V4.1）數據來表示地形表面，詳細及粗糙DTM數據均由該數據獲得；②一般情況下，參考面DTM可通過對3″SRTM（V4.1）數據進行平滑、降低分辨率而獲得，本文為了與EGM2008模型更匹配，參考面DTM數據采用EGM2008研發課題組完成的DTM2006.0模型；③利用①和②數據構建RTM數據，即zRTM＝z2－z1＝HSRTM－HDTM2006.0。

DTM2006.0模型是為滿足新一代超高階地球重力場模型（如：EGM2008）研發所涉及的與地形有關的重力量計算而完成的全球高分辨率DTM模型。該模型階次為2190，容納了大約240萬對正規化的高程系數，地面上任意一點P的模型高程可用式（5）求得［15］

2.4 綜合EGM2008模型和SRTM／DTM2006.0剩余地形模型的高程異常擬合

基于EGM2008模型，利用式（1）可獲得任意GPS點的重力場模型高程異常ζEGM2008，由式（4）可獲得剩余地形模型高程異常ζRTM，則任意點P總的高程異常為

式中，ζres表示殘余高程異常。為提高GPS高程轉換精度，在實際數據處理過程中，一般可采用如下方法：①基于EGM2008模型計算所有GPS點重力場模型高程異常ζEGM2008；②根據RTM數據計算所有GPS點的RTM高程異常ζRTM；③采用優化選擇法選擇少量GPS／水準點，計算這些點的實測高程異常，然后扣除模型高程異常ζEGM2008以及RTM高程異常ζRTM，獲得殘余高程異常ζres，利用一定的擬合方法建立殘余高程異常擬合模型，從而推求待求GPS點殘余高程異常，通過式（6）即可獲得待求點的總高程異常。其中ζRTM和ζres可表示為EGM2008模型的截斷誤差。選擇少量GPS／水準點的目的有兩個：①建立殘余高程異常擬合模型；②實現高程基準轉換，即將EGM2008模型所對應的高程基準轉換到水準點對應的高程基準。

對于面狀測區，待求點的殘余高程異常可通過平面或曲面等擬合函數獲得。而對于公路、鐵路等線路測區，由于GPS點一般沿線路延伸方向呈近似直伸形式分布，若采用平面或曲面擬合函數，則需要在線路兩側聯測足夠多的GPS／水準點，否則容易造成觀測方程病態。一般可選擇某一GPS點作為坐標原點，以線路延伸方向作為坐標軸的橫軸，過坐標原點且與橫軸垂直的方向為縱軸建立線路獨立坐標系，然后忽略垂直線路方向的高程異常變化，并利用直線、二次曲線或三次曲線擬合函數。在實際計算中選擇不同的GPS／水準點，可能會得到不同精度的高程轉換結果，因此有必要對GPS／水準點進行優化選擇，可采用逐步剔除法［16］。

3 計算與結果分析

為對本文提出的方法進行精度分析，現選擇華南地區兩個地形起伏相對較大測區GPS／水準數據作為測試數據：算例A（帶狀GPS／水準數據）和算例B（面狀GPS／水準數據）。

3.1 計算數據準備

為滿足實際計算需要，本文共收集或計算了如下幾組數據。

（1）3″的SRTM數據：算例A（25°N－30°N，107°E－115°E）和算例B（23°N－27°N，111°E－116°E）3″的SRTM數據，該數據直接從SRTM官方網站獲得，用于表示測區詳細地形。

（2）15″的SRTM數據：15″SRTM數據是通過對3″SRTM數據進行平滑處理獲得，用于表示測區的粗糙地形。

（3）參考面DTM數據：根據DTM2006.0模型利用式（5）計算了這兩個測區相應范圍的分辨率為30″的參考面DTM數據，模型階次取2160。

（4）高精度的GPS／水準點：算例A共收集到高精度GPS／水準點109個，算例B共收集到高精度GPS／水準點41個，兩個算例的GPS點坐標基準均為CGCS2000大地坐標系，GPS點大地高精度均優于2cm，水準測量等級不低于三等。

兩個測區地形情況（15″的SRTM數據）和GPS／水準點概略點位見圖3，其中十字絲表示GPS／水準點。測區RTM高程數據見圖4（SRTM高程與DTM2006.0高程之差，分辨率為30″）。

圖3 測區地形（SRTM）及GPS／水準點Fig.3 SRTM terrain and GPS／leveling points

圖4 測區RTM高程數據（SRTM－DTM2006.0）Fig.4 RTM elevations（SRTM－DTM2006.0）

3.2 算例A

算例A為某一鐵路客運專線GPS／水準數據，線路跨度近500km，共有GPS／水準點109個，整個線路西部地勢較東部高，東部GPS／水準點間距略比西部大，整個線路平均點間距約4.5km。為對殘余高程異常進行有效擬合，提高GPS高程轉換精度，本文將整個線路GPS／水準數據分為4段：A1、A2、A3和A4，分段的基本過程和原則為：①鑒于EGM2008地球重力場模型具有很高的精度，本文首先利用2160階次的EGM2008地球重力場模型計算所有GPS／水準點的模型高程異常；② 沿線路走向繪制線路模型高程異常變化圖；③ 在保證每段高程異常變化趨勢基本一致或變化不大前提下，盡量延長線路長度，且每個測段宜近似直線走向。A1、A2、A3和A4測段模型高程異常變化圖見圖5，每個圖橫軸表示該測段點序號，整個測區高程異常變化約9m。

圖5 算例A測區高程異常變化圖Fig.5 The height anomaly variability for case A

圖6是利用3″的SRTM數據（詳細DTM）、15″的SRTM數據（粗糙DTM）、30″的DTM2006.0數據（參考面DTM）以及GPS數據獲得的各點的RTM高程異常。在計算過程中，采用詳細／粗糙DTM數據組合的目的是為了提高計算效率。隨著積分半徑R2的增大，計算點P的RTM高程異常值會逐步趨于穩定。在積分半徑R2一定的前提下，理論上隨著積分半徑R1的增大，計算點P的RTM高程異常值的精度會越高，但計算效率越低。同時RTM高程異常值的變化范圍越小，若給定一個可容許的變化值（如1cm），則可從中確定積分半徑R1及R2。在實際計算中本文取R1＝40km，R2＝200km。圖6結果顯示，在該區域RTM高程異常在±0.05m范圍內，RTM高程異常最大值為0.042m，最小值為－0.044m，且東部GPS／水準點RTM高程異常變化較小，這與東部地形起伏相對較小相吻合（見圖3）。

圖6 RTM高程異常值（算例A）Fig.6 RTM height anomaly（case A）

表1為采用2160階次的EGM2008模型以及RTM數據計算得到的高程異常與GPS／水準點實測高程異常的比較結果。從表1中可以看到，2160階的EGM2008模型在該區域具有較高的精度，最差精度為6.0cm，最高精度為3.1cm。EGM2008高程異常與GPS／水準高程異常差值的平均值不為零，則暗示兩者存在系統偏差。該系統偏差主要是由EGM2008模型所定義的高程基準和我國國家高程基準不一致造成的。綜合EGM2008模型高程異常及RTM剩余高程異常后，GPS高程轉換的精度都有不同程度的提高，4個測段高程異常精度提高量分別為5%、10%、42%和2%，說明RTM技術能提高GPS高程轉換精度。

為了分析RTM高程異常及殘余高程異常擬合技術對GPS高程轉換的影響，在測段A1、A2、A3和A4中各選擇3～5個GPS／水準點。選點的原則有兩個：①要求每一測段中GPS／水準點平均點間距不小于25km；②采用逐步剔除法選擇最優的GPS／水準點。由于每個測段是近似直伸形式，故擬合模型采用二次曲線和三次曲線，即在每一測段都建立線路獨立坐標，并忽略垂直于線路方向的高程異常變化。每一測段所選擇的GPS／水準點個數及高程轉換精度統計結果見表2，高差轉換精度統計結果見表3。

表1 EGM2008／RTM高程異常與GPS／水準高程異常比較結果（算例A）Tab.1 The statistical results of GPS height anomaly differences between EGM2008／RTM and GPS／leveling（case A） cm

從表1和表2可以看出：

（1）與用EGM2008模型直接轉換相比，通過對僅扣除EGM2008模型高程異常的殘余高程異常進行擬合，對于二次曲線擬合，A1和A3測段精度略微降低，而A2和A4測段精度提高較大，約提高50%；對于三次曲線擬合，A3測段精度略微降低，而A1、A2和A4測段精度提高較大，除A1測段外，三次曲線擬合與二次曲線擬合精度基本相當。

（2）通過對扣除EGM2008／RTM模型高程異常的殘余高程異常進行擬合，精度能明顯提高，與未顧及RTM高程異常情況相比，若選擇擬合方法為二次曲線，則A1、A2、A3和A44個測段精度分別提高了19%、20%、71%和10%；若選擇擬合方法為三次曲線，則A1、A2、A3和A44個測段精度分別提高了33%、65%、76%和21%；顧及RTM高程異常后，無論采用二次曲線還是三次曲線擬合方法，精度都會得到不同程度的提高，最大可提高76%。測段A4精度提高量相對較小，這和測段A4及周圍地形起伏相對較小有關（見圖3），圖6結果也顯示該測段RTM高程異常變化較小。

（3）與利用EGM2008／RTM高程異常直接轉換結果相比，對殘余高程異常進行擬合，轉換精度也獲得了不同程度的提高。整體上，顧及RTM高程異常后，三次曲線擬合比二次曲線擬合精度略高，整體精度可以達到2.0cm，在A2和A3測段，精度甚至可達到1.0cm。

表2 基于EGM2008模型及RTM的GPS高程轉換精度統計結果（算例A）Tab.2 The statistical results of GPS height transformation accuracy based on EGM2008model and RTM（case A）cm

為了進一步分析本文方法獲得的高差精度，將表2中顧及RTM高程異常后的三次曲線擬合高程結果轉換為高差（沿線路坐標依次求高差），并和三等、四等水準測量限差要求進行比較，結果見表3。

表3結果顯示：除測段A1外，其余3個測段平均90%的高差滿足三等水準測量精度要求。測段A1、A2和A3中所有高差均滿足四等水準測量限差要求，A4測段中95%高差滿足四等水準測量限差要求，經分析發現有1個GPS點高程轉換精度較差，若扣除該點，則剩余高差均滿足四等水準測量限差要求。

表3 基于EGM2008模型及RTM的GPS高差轉換精度統計結果（算例A）Tab.3 The statistical results of GPS height difference transformation accuracy based on EGM2008 model and RTM（case A） cm

3.3 算例B

算例B是某一城市GPS／水準數據，控制面積600km2多，平均點間距約4km，整個GPS控制區域地形起伏相對不大，但周圍地形起伏很大，共有GPS／水準點41個。

圖7是利用SRTM數據、DTM2006.0數據以及GPS數據獲得的各點的RTM高程異常，計算方法同3.2節。圖7結果顯示，在該區域RTM高程異常在±0.03m范圍內，RTM高程異常最大值為0.03m，最小值為－0.021m。

表4為采用2160階次的EGM2008模型以及RTM數據計算得到的高程異常與GPS／水準點實測高程異常的比較結果。從表4中可以看到，2160階的EGM2008模型在該區域精度可達2.4cm。綜合EGM2008模型高程異常及RTM剩余高程異常后，GPS高程轉換的精度獲得了提高，其提高量為38%，說明RTM技術能提高GPS高程轉換精度。

同樣對于該測區，為了分析RTM高程異常及殘余高程異常擬合技術對GPS高程轉換的影響。根據測區情況，選擇平面和曲面函數建立殘余高程異常擬合模型，采用GPS／水準點優化選擇法確定已知點，個數分別為3個和6個。平面擬合和曲面擬合中約束點間距約25km和10km，計算結果見表5。

表4和表5結果顯示：與用EGM2008模型直接轉換相比，僅通過對扣除EGM2008模型高程異常的殘余高程異常進行擬合，則平面和曲面擬合的精度分別為2.5cm和2.3cm，而EGM2008模型精度為2.4cm，精度基本相當。顧及RTM高程異常后，GPS高程轉換精度顯著提高，平面和曲面殘余高程異常擬合方法的精度分別為1.2cm和0.9cm，精度分別提高了52%和61%。與利用EGM2008／RTM模型高程異常直接轉換結果相比，該方法精度可分別提高20%和40%。

表4 EGM2008／RTM高程異常與GPS／水準高程異常比較結果（算例B）Tab.4 The statistical results of GPS height anomaly differences between EGM2008／RTM and GPS／leveling（case B）cm

表5 基于EGM2008模型及RTM的GPS高程轉換精度統計結果（算例B）Tab.5 The statistical results of GPS height transformation accuracy based on EGM2008model and RTM（case B）cm

4 結 論

（1）EGM2008地球重力場模型所定義的高程基準和我國國家高程基準存在系統性偏差。在局部區域，2160階次的EGM2008地球重力場模型具有很高的精度。在算例A中，4個測段的精度從3.1～6.0cm不等，而在算例B中精度可達到2.4cm。

（2）利用SRTM／DTM2006.0構建RTM數據，并將其轉換為RTM高程異常，能在一定程度上對2160階次的EGM2008模型的截斷誤差進行補償。在算例A中，4個測段的精度提高量從2%～42%不等，在算例B中，精度提高約38%。

（3）若不考慮RTM高程異常，僅通過對扣除EGM2008模型高程異常后的殘余高程異常進行擬合，則能在一定程度上提高高程異常的計算精度。而通過對扣除EGM2008／RTM模型高程異常后的殘余高程異常進行擬合，高程異常擬合精度明顯提高。與未顧及RTM高程異常情況相比，算例A中高程精度平均提高量超過40%，高程平均精度優于2.0cm（三次曲線擬合），而算例B中高程精度平均提高量超過50%，高程平均精度約1.0cm。在算例A中平均90%的高差滿足三等水準測量精度要求，近乎100%的高差滿足四等水準測量精度要求。與利用EGM2008／ RTM模型高程異常直接轉換結果相比，由該方法獲得的高程異常的精度也獲得了較大程度的提高，精度提高的主要原因是RTM高程異常表示高程異常中的高頻部分，扣除該部分后，有利于殘余高程異常的建模。

（4）在算例A和算例B中，計算了GPS／水準點平均點間距不低于25km的情況，為利用少量GPS／水準點進行較大范圍的局部似大地水準面精化提供一種參考。本文計算所利用的數據源均可以免費獲得（EGM2008模型、SRTM數據及DTM2006.0模型）。利用少量的GPS／水準點，即可獲得較為滿意的高程轉換結果。

［1］ PAVLIS N K，HOLMES S A，KENYON S C，et al.An Earth Gravitational Model to Degree 2160：EGM2008［R］.Vienna：EGU General Assembly 2008，2008.

［2］ GRUBER T.Evaluation of the EGM2008Gravity Field by Means of GPS Leveling and Sea Surface Topography Solutions［J］.Newton’s Bulletin，2009（4）：3－17.

［3］ JEKELI C，YANH H J，KWON J H.Evaluation of EGM08—Globally and Locally in South Korea［J］.Newton’s Bulletin，2009（4）：38－49.

［4］ ZHANG Chuanyin，GUO Chunxi，CHEN Junyong，et al.EGM2008and Its Application Analysis in Chinese Mainland［J］.Acta Geodaetica et Cartographica Sinica，2009，38（4）：283－289.（章傳銀，郭春喜，陳俊勇，等.EGM2008地球重力場模型在中國大陸適用性分析［J］.測繪學報，2009，38（4）：283－289.）

［5］ ZHANG Chuanyin，CHAO Dingbo，DING Jian，et al.Arithmetic and Characters Analysis of Terrain Effects for cm Order Precision Height Anomaly［J］.Acta Geodaetica et Cartographica Sinica，2006，35（4）：308－314.（章傳銀，晁定波，丁劍，等.厘米級高程異常地形影響的算法及特征分析［J］.測繪學報，2006，35（4）：308－314.）

［6］ HIRT C，FEATHERSTONE W E，MARTI U.Combining EGM2008and SRTM／DTM2006.0Residual Terrain Model Data to Improve Quasigeoid Computations in Mountainous Areas Devoid of Gravity Data［J］.Journal of Geodesy，2010，84（9）：557－567.

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

［8］ HIRT C，MARTI U，BüRKI B.Assessment of EGM2008 in Europe Using Accurate Astrogeodetic Vertical Deflections and Omission Error Estimates from SRTM／DTM2006.0Residual Terrain Model Data［J］.Journal of Geophysical Research，2010，115（B10404）：1－13.

［9］ FORSBERG R，TSCHERNING C C.The Use of Height Data in Gravity Field Approximation by Collocation［J］.Journal of Geophysical Research，1981，86（B9）：7843－7854.

［10］ FORSBERG R.A study of Terrain Reductions，Density Anomalies and Geophysical Inversion Methods in Gravity Field Modelling［R］.Report 355.Columbus：Ohio State University，1984.

［11］ FORSBERG R.Terrain Effects in Geoid Computations［R］.Copenhagen：National Survey and Cadastre，1994.

［12］ HOLMES S A，FEATHERSTONE W E.A Unified Approach to the Clenshaw Summation and the Recursive Computation of Very High Degree and Order Normalised Associated Legendre Functions［J］.Journal of Geodesy，2002，76（5）：279－299.

［13］ 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.

［14］ NAGY D，PAPP G，BENEDEK J.Corrections to“The Gravitational Potential and Its Derivatives for the Prism”［J］.Journal of Geodesy，2002，76（8）：475.

［15］ PAVLIS N K，FACTOR J K，HOLMES S A.Terrainrelated Gravimetric Quantities Computed for the Next EGM［C］∥Proceedings of the 1st International Symposium of the International Gravity Field Service：18.Istanbul：［s.n.］，2007：318－323.

［16］ SHEN Yunzhong，GAO Dakai，ZHANG Xingfu.Optimal Choice of GPS Leveling Points［J］.Geotechnical Investigation ＆Surveying，2004（6）：49－51.（沈云中，高達凱，張興福.GPS水準點優化選擇法［J］.工程勘察，2004（6）：49－51.）