剩余地形模型高程异常计算的积分法及精度分析

张永毅　张兴福　周波阳　魏德宏　丘广新

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。