球坐标系下基于全球DEM 数据的重力远二区高精度地改算法

张 伟 廖国忠 张秋冬 李 华

1 中国地质调查局成都地质调查中心，成都市一环路北三段2号，610082

2 成都理工大学地球物理学院，成都市二仙桥东三1号，610059

随着重力传感器制造工艺的不断发展，μGal级的重力测量工作已成为可能。文献［1－2］对影响重力地改精度的主要原因进行探讨，认为误差精度取决于数字地形高程与真实地形地貌的逼近程度。数字地形高程的网格距越小，高程值越精确，其地形改正值的计算精度就越高［3－4］。而对于球坐标系下重力异常的计算，文献［5，6］从不同角度推导了相同的解析式。本文在球坐标系下重力远区地形改正理论模型的基础上，提出使用公开的全球高分辨率DEM 数据计算远二区地改。实例表明，该方法能够大大提高远二区地改的计算精度，达到μGal量级。

1 球坐标系下的远二区地改理论

1.1 球坐标下的远区地改计算公式

球坐标系下的球体理论模型如图1所示，地球球心O为原点，N、S点分别为地球北、南极，P、N、S在同一个过球心O的平面上，点M在空间曲面上，以OP的延长线方向为Z轴，OX轴垂直于PNS平面，OY轴垂直于XOZ平面，X、Y、Z轴满足右手坐标系。α1、α2、α3可看成以P为北极的“经度”，β1、β2为其“纬度”，空间曲面上的点M可表示为Μ（α，β，h）。在球极坐标系下，随地形起伏的质量体Ω可描述为地壳扇形块h）dv，其底面是以地球平均半径R所作的球面，顶面是地形的自由起伏面。当积分区域为Ω＝｛（α，β，h）｜α1≤α≤α2，β1≤β≤β2，0≤h≤H｝时，质量体Ω相对于测点P的地形改正值为：

图1 球坐标系下的重力远二区地改模型Fig.1 Far zone II’s terrain correction model in spherical coordinates

式中，R为地球平均半径，取大地水准面；z、H为P、M点的海拔高程；G为引力常数；ρ为地壳平均密度。

式（1）中的积分函数对变量α而言是常量积分，当α∈［0，2π］时，ΔgΩ变为环带（β1，β2）相对于测点P的 重力地改值Δg［β1，β2］。将式（1）中的变量h设为常量，取h＝，为积分区域的平均高程，此时式（1）可简化为［6］：

再令

于是，

1.2 地壳扇形“质量体”的高程

计算远二区的地形改正值时，采用的高程资料是全国5′×5′高程节点网，以每个地壳扇形“质量体”落入5′×5′高程数据个数相加求和取平均值作为该扇形体的平均高程。但由于这类数据网格距较大，式（2）中与真实地貌间必然存在较大误差。而全球ASTER GDEM 数据的分辨率能达到3″×3″，从而大大提高远二区地改的精度。ASTER GDEM 和SRTM3 DEM 数据基于WGS－84坐标系，该坐标系下的重力远区地改模型如图2所示。当β∈［β1，β2］，α∈［0，2π］时，M点的“轨迹”是球面上以OP为轴的环带，可将式（5）离散化成式（6），即将环带等分为M个网格，M值越大，环带上的网格划分越精细，对应式（2）中的越逼近真实地形。当M→∞时，式（6）便变为理论解析式（2）。

图2 WGS－84坐标下的重力远区地改模型Fig.2 Far zone II’s terrain correction model in WGS－84coordinates

M为某环上的网格数，而P点远区地改值可表示为二重离散积分公式。式（7）可由计算机编程实现：

N为20～166.7km间的环带数。

代入式（7）前的一个关键问题是将WGS－84坐标系下的经纬度坐标（θ，Φ）转换为球坐标下的坐标（α，β）。根据球面三角形公式［7］，两坐标系之间的转换关系可由式（8）和式（9）表示。如图2所示，α的取值以PN方向为起始0°，令方位角按相对于球心O的顺时针方向增加。当M点位于平面POY前方时，即θ0－θ＞0，α的取值范围为180°～360°；当M点位于平面POY后方时，即θ0－θ≤0，α的取值范围为0～180°。

2 计算方法与实现

目前，基本覆盖全球的DEM 数据主要有：ASTER GDEM V2数据、SRTM3 V4.1数据、GTOPO30数据以及我国的资源三号测绘卫星立体影像数据等。其中，ASTER GDEM V2数据的分辨率最高（1″×1″，约30m×30m），SRTM3次之（3″×3″，约90m×90m）。两种数据历经几次版本升级后，质量日趋精良，绝对高度误差小于16m，相对高度误差小于10m，满足重力远区地改的高程精度要求。

式（7）中，M、N越大，其环带和环带上的网格剖分越精细，计算结果越逼近解析值，但计算量也越大。区域重力调查规范推荐的分环和网格参数如图3和表1所示。

表1 区域重力调查规范中远二区计算的环带和方位数（R＝6 371.025km）Tab.1 The subdivision of zone and azimuth in regional gravity survey specification（R＝6 371.025km）

图3 区域重力调查规范中环带及其网格划分模型Fig.3 Zone subdivision model and its mesh in regional gravity survey specification

M、N和网格剖分方案确定后，即可根据α、β计算出某一网格4个节点在WGS－84坐标系下的经纬度值。根据球面三角形边余弦定理：

联立式（9）和（10），可求得网格4个节点的经纬度坐标。当网格较小时，可将4个节点视为同平面，通过平面插值方法，利用待插值点邻近四周的DEM 数据得到节点的高程值。将4 个节点的高程值取平均，便可得到该网格的平均高程。代入式（7），得到测点在球坐标系下的远区地形改正值。由于式（9）存在除数为0时的奇异值，因此在程序编写时需要进行特殊处理，当cosΦ＝0时，Φ＝90°，M点处于北极或者南极，θ＝0°；当cosΦ0＝0时，Φ0＝90°，P点处于北极或者南极，θ＝α。

3 算例分析

3.1 理论模型

选用的理论模型为高程＝1 000m 的有限球壳，且计算点高程z＝1 000m。将、z、β1（β1＝20km／R）、β2（β2＝166.7km／R）代入式（2），可求得有限球壳模型的理论地改值为－3.749 57mGal。

在有限球壳模型下，将不同的环数和方位数代入式（5），计算值与理论值之间的差值如表2所示。可知，在不考虑高程误差时，式（5）的计算精度主要取决于环数，与方位数无关，环数越大，其计算精度越高，当N＞4 000时，式（5）的计算精度达到μGal级。ASTER GDEM V2 数据网格30m×30 m，其环数最大取值可到5 500，满足μGal级远二区地改要求。

表2 有限球壳理论模型下不同环数及方位的试算结果对比Tab.2 Comparative results within different subdivisions in limited spherical shell model

3.2 误差精度

为进一步分析高程误差对本文计算方法的影响，在上述有限球壳模型（＝1 000m，z＝1 000 m）基础上，在参与计算的每个高程网格节点上加入高斯随机噪声进行试算，高斯噪声的均值为0，标准差为20，其噪声水平与ASTER GDEM V2 DEM 的数据质量相当，计算结果如表3所示。从表3可知，在该噪声水平下，加入噪声和没有加入噪声时的计算结果差值在μGal量级以下，可以忽略。因此，ASTER GDEM V2和SRTM3数据的质量满足于μGal级远二区地改的工作要求。如表3第10行所示，当N＝5 500、M＝35 000时，有噪声与无噪声时的计算结果差值小于μGal数量级。

表3 有限球壳理论模型在高斯噪声条件下的试算结果对比Tab.3 Comparative results within different subdivisions under Gauss noise environment

［1］李东汉.试论区域重力测量远区地改的精度［J］.物探与化探，1984，8（2）：99－104（Li Donghan.A Discussion on the Accuracy of Far Sector Topographic Correction in Regional Gravity Survey［J］.Geophysical ＆ Geochemical Exploration，1984，8（2）：99－104）

［2］张俊，张宝松，邸兵叶，等.高程数据网格间距对重力中区地形改正精度的影响［J］.物探与化探，2014，38（1）：157－161（Zhang Jun，Zhang Baosong，Di Bingye，et al.The Effect of the Grid Spacing of Elevation on the Accurcy of Median Region Terrain Correction of Gravity［J］.Geophysical ＆Geochemical Exploration，2014，38（1）：157－161）

［3］赵海涛，张兵，左正立，等.中国及周边区域ASTER GDEM 与SRTM DEM 高程对比分析及互补修复［J］.测绘科学，2012，37（1）：8－11（Zhao Haitao，Zhang Bing，Zuo Zhengli，et al.Elevation Comparison and Complementary Repair of ASTER GDEM and SRTM DEM in China and Surrounding Areas［J］.Science of Surveying and Mapping，2012，37（1）：8－11）

［4］杜小平，郭华东，范湘涛，等.基于ICESat／GLAS数据的中国典型区域SRTM 与ASTER GDEM 高程精度评价［J］.地球科学—中国地质大学学报，2013，38（4）：887－897（Du Xiaoping，Guo Huadong，Fan Xiangtao，et al.Verical Accuracy Assesment of SRTM and ASTER GDEM over Typical Regions of China Using ICESat／GLAS［J］.Earth Science－Journal of China University of Geosciences，2013，38（4）：887－897）

［5］杜劲松，陈超，梁青，等.球冠体积分的重力异常正演方法及其与Tesseriod单元体泰勒级数展开方法的比较［J］.测绘学报，2012，41（3）：339－346（Du Jinsong，Chen Chao，Liang Qing，et al.Gravity Anomaly Calculation Based on Volume Integral in Spherical Cap and Comparision with the Tesserio－Taylor Series Expansion Approach［J］.Acta Geodaetica et Cartographica Sinca，2012，41（3）：339－346）

［6］安玉林，张明华，黄金明，等.纯球坐标系内各项重力校正值计算方案和计算过程［J］.物探与化探，2010，34（6）：697－705（An Yulin，Zhang Minghua，Huang Jinming，et al.The Computation Scheme and Computation Process for Gravity Correction Values Within the Pure Spherical Coordinate System［J］.Geophysical ＆Geochemical Exploration，2010，34（6）：697－705）

［7］周雪娟，褚善东.球面三角形公式及其应用［J］.浙江国际海运职业技术学院学报，2008，4（1）：59－62（Zhou Xuejuan，Chu Shandong.Formula of Spherical Triangle and Its Application［J］.Journal of Zhejiang International Maritime College，2008，4（1）：59－62）