剩余地形模型高程異常計算的積分法及精度分析

張永毅　張興福　周波陽　魏德宏　丘廣新

1　廣東工業大學土木與交通工程學院，廣州市外環西路100號，510006 2　廣州市城市規劃勘測設計研究院，廣州市建設大馬路10號，510006

?

剩余地形模型高程異常計算的積分法及精度分析

張永毅1張興福1周波陽1魏德宏1丘廣新2

1廣東工業大學土木與交通工程學院，廣州市外環西路100號，510006 2廣州市城市規劃勘測設計研究院，廣州市建設大馬路10號，510006

以構建剩余地形模型高程異常的數字地形模型的分辨率及其參考面的選擇為研究對象，系統分析兩者對剩余地形模型高程異常計算效率及精度的影響。實驗結果表明：1)將DTM2006.0模型作為參考面時，在海岸帶區域產生較大的誤差，而在陸地與RET2012和RET2014模型的計算結果相差不大；2)在構建我國東部地區剩余地形模型高程異常時，為保證計算效率及精度，計算時內外圈的積分半徑分別取50 km和200 km，SRTM數據的分辨率分別采用7.5″和15″，參考面模型使用RET2012。

地球重力場模型；剩余地形模型；數字地形模型；參考面；高程異常

由于重力場模型階次有限，存在模型截斷誤差，故通常采用剩余地形模型(RTM)對其進行補償[1-3]。RTM數據表示的是真實地形表面與參考地形表面的差值[4]。實際計算中，不同分辨率的DTM數據以及參考面DTM數據對所構建的RTM數據會產生不同程度的影響，因此兩者的選擇至關重要。本文從影響RTM計算精度的兩個因素出發，著重討論利用不同分辨率的地形數據與不同參考面地形數據分別構建RTM數據，并計算RTM高程異常，最后利用實測高精度GPS/水準數據對計算結果進行檢核。

1　原理與方法

1.1RTM構建

采用SRTM數據[5]以及對應的參考面DTM數據進行計算。RTM計算公式為：

(1)

式中，HRTM(i)為第i個格網點的剩余高程，HSRTM(i)為第i個格網點的SRTM高程，Href(i)為第i個格網點的參考面高程。

目前有兩種方法可以計算參考地形表面。方法一是利用一定階次的球諧地形模型數據計算，如DTM2006.0、RET2012及RET2014模型，其中DTM2006.0和RET2012模型為2 160階次的球諧地形模型，RET2014模型為10 800階次的球諧地球地形模型。計算公式為[6]:

(2)

方法二是對高分辨率的SRTM地形數據進行平滑(一般可取n×n個格網平均值)或結合低通濾波器獲得較為光滑的參考面。

若選擇與某一地球重力場模型分辨率一致的參考面，則意味著RTM數據可以表示比該地球重力場模型分辨率更高的高頻重力信息，使得該地球重力場模型全波段信號得以恢復。

1.2RTM高程異常計算

根據Forsberg[4]的研究，可以利用棱柱積分法或快速傅立葉變換將RTM轉換成RTM高程異常。由于快速傅立葉變換是近似計算，雖然計算速度較快，但精度相對較差[3]，因此，本文僅探討利用棱柱積分法將RTM數據轉換為RTM高程異常。每個棱柱(格網)對應的引力位為[4]：

V=Gρ0|||xyln(z+r)+yzln(x+r)+

(3)

式中，G為引力常數，r為坐標原點到點(x,y,z)的距離，(x1,x2,y1,y2,z1,z2)為棱柱體的邊界角點坐標，z2-z1表示剩余高程，ρ0=2 670 kg/m3為標準的地形質量密度。式(3)是平面近似，計算中需考慮地球曲率的影響[4]。其他變量的含義及具體計算過程可參考文獻[4,7]。將棱柱引力位轉換為對應的高程異常：

(4)

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

(5)

由于陸地巖石和海水密度不一樣，因此在陸海交界區域利用式(3)～(5)計算RTM高程異常時存在諸多不便，可通過壓縮因子將海水深度壓縮成等效的巖石高[8]：

(6)

式中，H′為壓縮后的海水深度，H為原海水深度，ρw=1 030 kg/m3為海水密度，ρ0=2 670 kg/m3為地形質量密度。利用式(6)可以將海岸帶區域計算RTM高程異常所需要的海洋水深數據轉換為等效的巖石高，方便利用式(3)～(5)進行相關計算。

2　計算與分析

為對影響RTM構建的數字地形模型和參考面模型進行分析，收集以下幾組數據(范圍均為18°～32°N,107°～121°E)：1)由SRTM官方提供的7.5″和15″ SRTM 數據；2)根據式(6)并利用2 160 階的DTM2006.0、RET2012和RET2014等地形模型分別獲得的15″分辨率的參考面模型。

2.1RTM構建中參考面的影響

為了分析參考面對RTM構建的影響，利用15″的SRTM數據分別與15″的DTM2006.0、RET2012和RET2014模型數據(參考面DTM)構建RTM數據，并通過數值積分計算求得RTM高程異常值。計算區域為20°～30°N、109°～119°E范圍內的陸地，計算點的分辨率為1′，結果見圖1。

圖1　RTM高程異常值Fig.1　RTM height anomaly

從圖1可看出，利用RET2012和RET2014模型數據作為參考面計算所得的RTM高程異常值相當，而以DTM2006.0模型數據為參考面計算所得的RTM高程異常值在海岸帶區域與其他兩個模型對應的結果相差較大。為進一步分析參考面對RTM構建的影響，分別選擇一定范圍的陸地和海岸帶區域的RTM高程異常數據進行分析。在陸地上選取范圍為25°～29°N、111°～116°E的區域，同時在海岸帶附近選取范圍為20°～23°N、109°～116°E的區域分別進行分析，并以利用RET2012模型數據為參考面所求得的RTM高程異常值為參考，分別與利用DTM2006.0和RET2014模型數據為參考面所求得的RTM高程異常值進行對比分析，結果見表1和表2。

表1　RTM高程異常比較結果(陸地)

表2　RTM高程異常比較結果(海岸帶)

從表1看出，在陸地區域以DTM2006.0、RET2012和RET2014模型數據為參考面計算得到的RTM高程異常數據相差不大。因此，在我國東部陸地區域構建RTM數據時可任意選擇3個參考面中的一個。從表2看出，在海岸帶區域分別以DTM2006.0和RET2012模型數據為參考面計算所得到的RTM高程異常數據相差較大，最大值為38 cm，最小值為-1.55 cm；而分別以RET2014和RET2012模型數據為參考面所構建的RTM高程異常數據相差較小，最大值為0.07 cm，最小值為-0.40 cm。結合圖1和表2發現，以DTM2006.0模型數據為參考面所構建的RTM高程異常值在海岸帶與實際地形相差較大，而以RET2012和RET2014模型數據所構建的RTM高程異常值與實際地形比較吻合，這與3個參考面模型構建時所采用的數據源不同有關(在海岸帶區域構建RTM數據可選擇RET2012和RET2014兩個參考面中的一個)。

2.2RTM數據構建的效率及精度分析

利用棱柱積分法將RTM數據轉換到剩余高程異常時，為提高棱柱積分法的計算速度，距離計算點較遠的區域可采用分辨率稍低的DTM格網數據，距離計算點較近的區域采用分辨率較高的DTM格網數據[9]。實際計算過程中，一般會提供兩個積分計算半徑r1和r2，當積分半徑小于r1時用分辨率較高的詳細格網數據，當積分半徑在r1和r2之間時用分辨率較低的粗糙格網數據。隨著積分半徑r2的增大，計算點的RTM高程異常值會逐步趨于穩定。在積分半徑r2一定的前提下，

隨著積分半徑r1的增大，計算點的RTM高程異常值的精度會提高，但計算效率會降低。因此，有必要在計算過程中選擇一個合適的積分半徑組合方案，使得計算精度和計算效率同時得到保證。

構建分辨率為5′的RTM數據進行分析。對于構建分辨率不高于7.5″的RTM數據，可采用7.5″的SRTM數據作為詳細DTM格網數據，15″的SRTM數據作為粗糙DTM格網數據。采用7.5″與15″的SRTM數據和15″的RET2012模型數據(參考面DTM)，并結合多種不同的積分半徑組合方案(表3)進行RTM高程異常值的計算，計算區域為20°～30°N、109°～119°E的陸地，結果見圖2。

表3　積分半徑組合方案

理論上，方案5計算的RTM高程異常數據精度最高但計算效率最低。為分析其他積分半徑組合方案對RTM高程異常構建的效率及精度的影響，以方案5為參考標準，其他方案所計算的結果均與其作差比較，比較結果見圖3。其中方案A、B、C、D分別為方案1、2、3、4與方案5的RTM高程異常值的較差結果，見表4。

圖2　RTM高程異常值Fig.2　RTM height anomaly

圖3　RTM高程異常比較結果Fig.3　The comparison results of RTM height anomaly

由表4可知，B、C、D(r2為200 km，r1分別為50、75、100 km)3種方案的RTM高程異常值相差較小，但計算效率相差較大，其計算效率提高量分別為46.7%、34.8%、19.6%。說明r1大于50 km后計算點的RTM高程異常值的精度逐步趨于穩定，但計算效率有著不同程度的降低，所以選擇方案B(r1=50 km，r2=200 km)既能保證計算精度，也能提高計算效率。若對計算精度的要求較低，可以選擇方案A(r1=25 km，r2=200 km)，這將使計算效率大幅提高。

表4　4種方案RTM高程異常比較結果

2.3RTM高程異常計算結果檢核分析

為檢核所構建的RTM高程異常數據，收集某地區的44個高精度GPS/水準點，分布范圍為24.6°～24.9°N、113.4°～113.7°E，其中高程異常最大值為-8.291 m，最小值為-9.180 m，平均值為-8.741 m。為了驗證§2.1的分析結果，分別以DTM2006.0、RET2012和RET2014模型數據為參考面計算得到RTM高程異常數據，根據§2.2的分析結果，進行計算時選擇r1=50 km，r2=200 km的積分半徑組合。 檢核步驟如下：

1)構建RTM高程異常數據。利用7.5″的SRTM數據(詳細DTM)和15″的SRTM數據(粗糙DTM)分別與15″的DTM2006.0、RET2012和RET2014模型數據(參考面DTM)構建RTM數據，并通過積分計算求得RTM高程異常值(下文分別以RTM_a、RTM_b、RTM_c表示)，計算范圍為20°～30°N、109°～119°E，積分半徑為r1=50 km，r2=200 km，計算點的分辨率為1′。

2)利用重力場模型計算GPS/水準點的高程異常。分別利用EIGEN-6C4(2 160階)和EGM2008(2 160階)模型計算GPS/水準點的高程異常，并與實測高程異常進行比較分析。

3)綜合重力場模型及RTM數據得到的高程異常與GPS/水準點實測高程異常進行比較分析。

檢核結果見表5和表6。表5為采用2 160階次的EIGEN-6C4模型以及RTM數據計算得到的高程異常與GPS/水準點實測高程異常的比較結果。從表5可以看出，2 160階的EIGEN-6C4模型高程異常在該區域的精度為2.3 cm。結合EIGEN-6C4模型高程異常及RTM高程異常后，GPS高程轉換的精度有著不同程度的提高，其中結合RTM_a(以DTM2006.0為參考面)、RTM_b(以RET2012為參考面)和RTM_c(以RET2014為參考面)后的高程異常精度提高量分別為26%、30%、30%。

表5　EIGEN-6C4/RTM高程異常與GPS/水準高程異常比較結果

表6　EGM2008/RTM高程異常與GPS/水準高程異常比較結果

表6為采用2 160階次的EGM2008模型以及RTM數據計算得到的高程異常與GPS/水準點實測高程異常的比較結果。從表6看出，2 160階的EGM2008模型高程異常在該區域的精度為2.3 cm。結合EGM2008模型高程異常及RTM高程異常后，GPS高程轉換的精度有著不同程度的提高，其中結合RTM_a、RTM_b、RTM_c后的高程異常精度提高量分別為43%、39%、39%。

綜合表5和表6看出：1)綜合重力場模型高程異常及RTM高程異常能提高GPS高程轉換的精度；2)在研究區域內，EGM2008模型高程異常結合RTM高程異常對GPS高程轉換的精度較EIGEN-6C4模型高程異常結合RTM高程異常稍高；3)在研究區域內，RTM_a、RTM_b和RTM_c高程異常的精度基本相當。

3　結　語

1)在研究區域內，以DTM2006.0、RET2012和RET2014模型數據為參考面計算得到的RTM高程異常值在陸地區域相差不大。在海岸帶，以DTM2006.0為參考面的計算結果較差，以RET2012和RET2014為參考面的計算結果與實際比較符合，這與3種模型構建的數據源有較大關系。

2)在研究區域內構建分辨率不高于7.5″RTM高程異常數據時，可采用7.5″和15″SRTM數據來構建RTM數據，選擇r1=50 km，r2=200 km的積分半徑組合可同時保證計算精度及效率。

3)利用實測高精度GPS/水準點高程異常數據對RTM高程異常數據進行檢核，結果表明，綜合不同重力場模型高程異常及不同RTM高程異常對GPS高程轉換的精度有著不同程度的提高，其中EGM2008模型高程異常結合RTM高程異常對GPS高程轉換的精度較EIGEN-6C4模型高程異常結合RTM高程異常稍高。

[1]Hirt C, Featherstone 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

[2]Hirt C. RTM Gravity Forward-Modeling Using Topography/Bathymetry Data to Improve High-Degree Global Geopotential Models in the Coastal Zone[J]. Marine Geodesy, 2013, 36(2):183-202

[3]Hirt C, Kuhn M, Claessens S, et al. Study of the Earth’s Short-Scale Gravity Field Using the ERTM2160 Gravity Model[J]. Computers & Geosciences, 2014, 73: 71-80

[4]Forsberg R. Terrain Effects in Geoid Computations[R].International School for the Determination and Use of the Geoid,Milano,1994

[5]Hirt C, Filmer M S, Featherstone W E. Comparison and Validation of the Recent Freely Available ASTER-GDEM Ver1, SRTM Ver4.1 and GEODATA DEM-9S Ver3 Digital Elevation Models over Australia[J]. Australian Journal of Earth Sciences, 2010, 57(3): 337-347

[6]Pavlis N K, Factor J K, Holmes S A. Terrain-Related Gravimetric Quantities Computed for the Next EGM[C]. The 1st International Symposium of the International Gravity Field Service (IGFS), Istanbul, 2007

[7]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

[8]Hirt C, Kuhn M, Featherstone W E, et al. Topographic/Isostatic Evaluation of New-Generation GOCE Gravity Field Models[J]. Journal of Geophysical Research: Solid Earth, 2012, 117(B5)

[9]Farr T G, Rosen P A, Caro E, et al. The Shuttle Radar Topography Mission[J]. Reviews of Geophysics, 2000, 45(2): 37-55

Foundation support:National Natural Science Foundation of China, No. 41104002, 41504013；Basic Research Fund of Geomatics of Key Laboratory of Geospace Environment and Geodesy, Ministry of Education, No. 10-01-07, 14-01-03.

About the first author:ZHANG Yongyi, postgraduate, majors in surveying data processing and soft developing, E-mail: yyzhang199107@163.com.

The Integral Method and Accuracy Analysis of Residual Terrain Model Height Anomaly

ZHANGYongyi1ZHANGXingfu1ZHOUBoyang1WEIDehong1QIUGuangxin2

1School of Civil and Transportation Engineering, Guangdong University of Technology, 100 West-Waihuan Road, Guangzhou 510006,China 2Guangzhou Urban Planning & Design Survey Research Institute, 10 Jianshedama Road, Guangzhou 510006, China

This paper studies the resolution of the digital terrain model, the selection of reference surface, and their influence on the calculation of speed and precision by RTM. The results show that:1) When using DTM2006.0 as a reference surface, RTM and RET2012 have similar results for mainland areas, while the precision of RTM is not as fine as that of RET2012 and RET2014 for coastal areas; 2) Considering computational efficiency and accuracy, we can use the 7.5″ SRTM for the inner zone with a radius of 50 km, the 15″ SRTM for the outer zone with a radius of 200 km, and choose the RET2012 model data as a reference surface to construct the residual terrain model at the east area of our country.

earth gravity field model; residual terrain model; digital terrain model; reference surface; height anomaly

ZHANG Xingfu, PhD, associate professor, majors in satellite gravity and GNSS data processing, E-mail: xfzhang77@163.com.

2015-09-11

張興福，博士，副教授，主要研究方向為衛星重力、GNSS數據處理，E-mail: xfzhang77@163.com。

10.14075/j.jgg.2016.09.005

1671-5942(2016)09-0770-05

P223

A

項目來源：國家自然科學基金(41104002；41504013)；地球空間環境與大地測量教育部重點實驗室測繪基礎研究基金(10-01-07,14-01-03)。第一作者簡介：張永毅，碩士生，主要研究方向為測量數據處理及軟件開發，E-mail：yyzhang199107@163.com。