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

張永毅　張興?！≈懿?yáng)　魏德宏　丘廣新

1　廣東工業(yè)大學(xué)土木與交通工程學(xué)院，廣州市外環(huán)西路100號(hào)，510006 2　廣州市城市規(guī)劃勘測(cè)設(shè)計(jì)研究院，廣州市建設(shè)大馬路10號(hào)，510006

?

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

張永毅1張興福1周波陽(yáng)1魏德宏1丘廣新2

1廣東工業(yè)大學(xué)土木與交通工程學(xué)院，廣州市外環(huán)西路100號(hào)，510006 2廣州市城市規(guī)劃勘測(cè)設(shè)計(jì)研究院，廣州市建設(shè)大馬路10號(hào)，510006

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

地球重力場(chǎng)模型；剩余地形模型；數(shù)字地形模型；參考面；高程異常

由于重力場(chǎng)模型階次有限，存在模型截?cái)嗾`差，故通常采用剩余地形模型(RTM)對(duì)其進(jìn)行補(bǔ)償[1-3]。RTM數(shù)據(jù)表示的是真實(shí)地形表面與參考地形表面的差值[4]。實(shí)際計(jì)算中，不同分辨率的DTM數(shù)據(jù)以及參考面DTM數(shù)據(jù)對(duì)所構(gòu)建的RTM數(shù)據(jù)會(huì)產(chǎn)生不同程度的影響，因此兩者的選擇至關(guān)重要。本文從影響RTM計(jì)算精度的兩個(gè)因素出發(fā)，著重討論利用不同分辨率的地形數(shù)據(jù)與不同參考面地形數(shù)據(jù)分別構(gòu)建RTM數(shù)據(jù)，并計(jì)算RTM高程異常，最后利用實(shí)測(cè)高精度GPS/水準(zhǔn)數(shù)據(jù)對(duì)計(jì)算結(jié)果進(jìn)行檢核。

1　原理與方法

1.1RTM構(gòu)建

采用SRTM數(shù)據(jù)[5]以及對(duì)應(yīng)的參考面DTM數(shù)據(jù)進(jìn)行計(jì)算。RTM計(jì)算公式為：

(1)

式中，HRTM(i)為第i個(gè)格網(wǎng)點(diǎn)的剩余高程，HSRTM(i)為第i個(gè)格網(wǎng)點(diǎn)的SRTM高程，Href(i)為第i個(gè)格網(wǎng)點(diǎn)的參考面高程。

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

(2)

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

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

1.2RTM高程異常計(jì)算

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

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

(3)

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

(4)

式中，γQ表示正常重力。點(diǎn)P所對(duì)應(yīng)的RTM高程異?？捎捎?jì)算區(qū)域內(nèi)所有單個(gè)棱柱對(duì)應(yīng)的高程異常總和表示：

(5)

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

(6)

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

2　計(jì)算與分析

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

2.1RTM構(gòu)建中參考面的影響

為了分析參考面對(duì)RTM構(gòu)建的影響，利用15″的SRTM數(shù)據(jù)分別與15″的DTM2006.0、RET2012和RET2014模型數(shù)據(jù)(參考面DTM)構(gòu)建RTM數(shù)據(jù)，并通過(guò)數(shù)值積分計(jì)算求得RTM高程異常值。計(jì)算區(qū)域?yàn)?0°～30°N、109°～119°E范圍內(nèi)的陸地，計(jì)算點(diǎn)的分辨率為1′，結(jié)果見(jiàn)圖1。

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

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

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

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

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

2.2RTM數(shù)據(jù)構(gòu)建的效率及精度分析

利用棱柱積分法將RTM數(shù)據(jù)轉(zhuǎn)換到剩余高程異常時(shí)，為提高棱柱積分法的計(jì)算速度，距離計(jì)算點(diǎn)較遠(yuǎn)的區(qū)域可采用分辨率稍低的DTM格網(wǎng)數(shù)據(jù)，距離計(jì)算點(diǎn)較近的區(qū)域采用分辨率較高的DTM格網(wǎng)數(shù)據(jù)[9]。實(shí)際計(jì)算過(guò)程中，一般會(huì)提供兩個(gè)積分計(jì)算半徑r1和r2，當(dāng)積分半徑小于r1時(shí)用分辨率較高的詳細(xì)格網(wǎng)數(shù)據(jù)，當(dāng)積分半徑在r1和r2之間時(shí)用分辨率較低的粗糙格網(wǎng)數(shù)據(jù)。隨著積分半徑r2的增大，計(jì)算點(diǎn)的RTM高程異常值會(huì)逐步趨于穩(wěn)定。在積分半徑r2一定的前提下，

隨著積分半徑r1的增大，計(jì)算點(diǎn)的RTM高程異常值的精度會(huì)提高，但計(jì)算效率會(huì)降低。因此，有必要在計(jì)算過(guò)程中選擇一個(gè)合適的積分半徑組合方案，使得計(jì)算精度和計(jì)算效率同時(shí)得到保證。

構(gòu)建分辨率為5′的RTM數(shù)據(jù)進(jìn)行分析。對(duì)于構(gòu)建分辨率不高于7.5″的RTM數(shù)據(jù)，可采用7.5″的SRTM數(shù)據(jù)作為詳細(xì)DTM格網(wǎng)數(shù)據(jù)，15″的SRTM數(shù)據(jù)作為粗糙DTM格網(wǎng)數(shù)據(jù)。采用7.5″與15″的SRTM數(shù)據(jù)和15″的RET2012模型數(shù)據(jù)(參考面DTM)，并結(jié)合多種不同的積分半徑組合方案(表3)進(jìn)行RTM高程異常值的計(jì)算，計(jì)算區(qū)域?yàn)?0°～30°N、109°～119°E的陸地，結(jié)果見(jiàn)圖2。

表3　積分半徑組合方案

理論上，方案5計(jì)算的RTM高程異常數(shù)據(jù)精度最高但計(jì)算效率最低。為分析其他積分半徑組合方案對(duì)RTM高程異常構(gòu)建的效率及精度的影響，以方案5為參考標(biāo)準(zhǔn)，其他方案所計(jì)算的結(jié)果均與其作差比較，比較結(jié)果見(jiàn)圖3。其中方案A、B、C、D分別為方案1、2、3、4與方案5的RTM高程異常值的較差結(jié)果，見(jiàn)表4。

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

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

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

表4　4種方案RTM高程異常比較結(jié)果

2.3RTM高程異常計(jì)算結(jié)果檢核分析

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

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

2)利用重力場(chǎng)模型計(jì)算GPS/水準(zhǔn)點(diǎn)的高程異常。分別利用EIGEN-6C4(2 160階)和EGM2008(2 160階)模型計(jì)算GPS/水準(zhǔn)點(diǎn)的高程異常，并與實(shí)測(cè)高程異常進(jìn)行比較分析。

3)綜合重力場(chǎng)模型及RTM數(shù)據(jù)得到的高程異常與GPS/水準(zhǔn)點(diǎn)實(shí)測(cè)高程異常進(jìn)行比較分析。

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

表5　EIGEN-6C4/RTM高程異常與GPS/水準(zhǔn)高程異常比較結(jié)果

表6　EGM2008/RTM高程異常與GPS/水準(zhǔn)高程異常比較結(jié)果

表6為采用2 160階次的EGM2008模型以及RTM數(shù)據(jù)計(jì)算得到的高程異常與GPS/水準(zhǔn)點(diǎn)實(shí)測(cè)高程異常的比較結(jié)果。從表6看出，2 160階的EGM2008模型高程異常在該區(qū)域的精度為2.3 cm。結(jié)合EGM2008模型高程異常及RTM高程異常后，GPS高程轉(zhuǎn)換的精度有著不同程度的提高，其中結(jié)合RTM_a、RTM_b、RTM_c后的高程異常精度提高量分別為43%、39%、39%。

綜合表5和表6看出：1)綜合重力場(chǎng)模型高程異常及RTM高程異常能提高GPS高程轉(zhuǎn)換的精度；2)在研究區(qū)域內(nèi)，EGM2008模型高程異常結(jié)合RTM高程異常對(duì)GPS高程轉(zhuǎn)換的精度較EIGEN-6C4模型高程異常結(jié)合RTM高程異常稍高；3)在研究區(qū)域內(nèi)，RTM_a、RTM_b和RTM_c高程異常的精度基本相當(dāng)。

3　結(jié)　語(yǔ)

1)在研究區(qū)域內(nèi)，以DTM2006.0、RET2012和RET2014模型數(shù)據(jù)為參考面計(jì)算得到的RTM高程異常值在陸地區(qū)域相差不大。在海岸帶，以DTM2006.0為參考面的計(jì)算結(jié)果較差，以RET2012和RET2014為參考面的計(jì)算結(jié)果與實(shí)際比較符合，這與3種模型構(gòu)建的數(shù)據(jù)源有較大關(guān)系。

2)在研究區(qū)域內(nèi)構(gòu)建分辨率不高于7.5″RTM高程異常數(shù)據(jù)時(shí)，可采用7.5″和15″SRTM數(shù)據(jù)來(lái)構(gòu)建RTM數(shù)據(jù)，選擇r1=50 km，r2=200 km的積分半徑組合可同時(shí)保證計(jì)算精度及效率。

3)利用實(shí)測(cè)高精度GPS/水準(zhǔn)點(diǎn)高程異常數(shù)據(jù)對(duì)RTM高程異常數(shù)據(jù)進(jìn)行檢核，結(jié)果表明，綜合不同重力場(chǎng)模型高程異常及不同RTM高程異常對(duì)GPS高程轉(zhuǎn)換的精度有著不同程度的提高，其中EGM2008模型高程異常結(jié)合RTM高程異常對(duì)GPS高程轉(zhuǎn)換的精度較EIGEN-6C4模型高程異常結(jié)合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

張興福，博士，副教授，主要研究方向?yàn)樾l(wèi)星重力、GNSS數(shù)據(jù)處理，E-mail: xfzhang77@163.com。

10.14075/j.jgg.2016.09.005

1671-5942(2016)09-0770-05

P223

A

項(xiàng)目來(lái)源：國(guó)家自然科學(xué)基金(41104002；41504013)；地球空間環(huán)境與大地測(cè)量教育部重點(diǎn)實(shí)驗(yàn)室測(cè)繪基礎(chǔ)研究基金(10-01-07,14-01-03)。第一作者簡(jiǎn)介：張永毅，碩士生，主要研究方向?yàn)闇y(cè)量數(shù)據(jù)處理及軟件開(kāi)發(fā)，E-mail：yyzhang199107@163.com。