有限范围的重力层间改正算法

马健^1,2, 魏子卿^2,3, 武丽丽⁴, 杨正辉⁵     

1. 信息工程大学地理空间信息学院, 河南 郑州 450052; 
2. 地理信息工程国家重点实验室, 陕西 西安 710054; 
3. 西安测绘研究所, 陕西 西安 710054; 
4. 信息工程大学导航与空天目标工程学院, 河南 郑州 450052; 
5. 长安大学地质工程与测绘学院, 陕西 西安 710054

收稿日期：2016-04-19； 修回日期：2016-12-29

基金项目：国家自然科学基金（41674025；41674082）；地理信息工程国家重点实验室开放研究基金（SKLGIE2016-M-1-5）

第一作者简介： 马健(1988-), 女, 博士生, 研究方向为物理大地测量。

E-mail:

摘要：层间改正是重力归算的一项重要内容，传统的平面层层间改正、球面层层间改正与地形改正的范围不一致，因此均存在远区虚拟地形引入的近似误差，且计算点高度越高，此误差越大。本文提出使用有限范围的层间改正进行重力归算的方法，使其区域范围与地形改正的范围一致。然后给出了有限范围层间改正的简便计算方法，该算法与通过地形改正严密积分法演化来的算法具有较好的一致性。内插试验说明当计算点地形高于1000 m时，内插应使用基于有限范围层间改正的重力归算方法。

关键词： 传统层间改正     有限范围层间改正     简便公式     地形改正严密积分演化公式    

The Bouguer Correction Algorithm for Gravity with Limited Range

MA Jian^1,2, WEI Ziqing^2,3, WU Lili⁴, YANG Zhenghui⁵     

Abstract: The Bouguer correction is an important item in gravity reduction, while the traditional Bouguer correction, whether the plane Bouguer correction or the spherical Bouguer correction, exists approximation error because of far-zone virtual terrain. The error grows as the calculation point gets higher. Therefore gravity reduction using the Bouguer correction with limited range, which was in accordance with the scope of the topographic correction, was researched in this paper. After that, a simplified formula to calculate the Bouguer correction with limited range was proposed. The algorithm, which is innovative and has the value of mathematical theory to some extent, shows consistency with the equation evolved from the strict integral algorithm for topographic correction. The interpolation experiment shows that gravity reduction based on the Bouguer correction with limited range is prior to unlimited range when the calculation point is taller than 1000 m.

Key words: the traditional Bouguer correction     the Bouguer correction with limited range     simplified formula     evolution equation from the strict integral algorithm for topographic correction    

通过重力归算量化地形质量对地表重力的引力作用，是离散重力数据格网化以及局部大地水准面精化的重要基础工作^[-]。重力归算具体包括空间改正、层间改正^[-]、地形改正(有的文献称为局部地形改正)^[-]3个部分。本文重点对层间改正展开研究。传统层间改正的布格层分平面层和球面层^[]两种类型。文献[-]给出了球面布格层层间改正与平面布格层层间改正相差约2倍的结论。布格层是高度为计算点高程的等高度地形层，无论是无限延伸的平面层还是覆盖全球的球壳的形式，由于其区域范围与地形改正的范围不一致，因此两种类型的层间改正中均含有远区虚拟地形引入的误差。

基于此，本文提出使用布格层的区域范围与地形改正的范围一致的层间改正的方法，推导了有限范围层间改正的简便计算公式，该公式与由地形改正严密积分算法演化的公式具有一致性，并且在数学上具有一定的创新意义和理论价值。最后通过试验得出球面布格层不适宜内插，计算点地形超过1000 m时，应使用基于有限范围层间改正的重力归算进行内插的结论。

1 传统层间改正算法--平面层、球面层

传统层间改正所使用的布格层有平面和球面两种类型。

平面布格层是计算点所在的平面与该点在大地水准面上对应点所在平面之间的无限延伸的地形层。当地面计算点位于大地水准面上方时，平面布格层位于计算点之下，若除去这部分地形质量的引力作用，重力值将比实测重力值小，因此此时层间改正符号为负。在以计算点为原点的柱空间坐标系下，积分半径为a、高度为Δh(计算点位于大地水准面之上时，Δh符号取正)的平面布格地形层对计算点的地形改正公式的近似结果为^[]

 (1)

式中，G为地球引力常数；ρ为地壳密度。实际计算时通常假设平面布格层是一个无限延伸的地形层，即a取无限大，因此平面无限布格层层间改正的实用公式为

 (2)

球面布格层是计算点所在的球面与该点在大地水准面上对应点所在的球面之间的球壳状的地形层。球面布格层的层间改正计算公式为^[]

 (3)

式中，R为地球平均半径。比较式(2)、式(3)可知，球面布格层改正在量级上约为平面布格层的2倍。

重力归算中的地形影响包含层间改正与地形改正两部分，由于实际计算时，地形改正仅采用有限区域范围，这使层间改正与地形改正计算范围不一致，即平面布格层和球面布格层层间改正均包含未进行地形改正的远区虚拟地形产生的误差。由于布格层的高度为计算点的高程，因此该误差对不同高度的计算点也不相同：计算点高度越大，其层间改正中远区虚拟地形产生的误差也越大。

为消除层间改正中远区虚拟地形产生的误差，本文提出使用有限范围的层间改正进行重力归算，令布格层的范围与归算时地形改正的经纬度范围一致。虽然式(1)将布格层限定为圆柱也是有限区域范围，但由于实际地形数据都是经纬度划分的网格形式，因而使用圆柱计算层间改正不能完全与地形改正的范围吻合，存在近似误差。本文将层间改正的区域设定为以经纬度网格划定的区域，较圆柱形式的布格层层间改正更为严密。

2 有限范围的层间改正算法推导

本文对有限范围的层间改正算法的推导是在地形改正算法的基础上进行的。不顾及地形改正符号时，空间直角坐标系下的地面点地形改正δΔg的计算公式为

 (4)

式中，(x_P, y_P, z_P)、(x, y, z)分别为计算点和流动点的三维直角坐标；为计算点与流动点之间的空间距离。对式(4)的z轴积分可得式(5)

 (5)

对于地形改正而言，不同流动点的z坐标是不同的，而层间改正可视作流动点高程均等于计算点高程的一种特殊的地形改正，因此当式(5)用于计算层间改正时，则需将流动点与计算点之差z_P-z固定为常数Δh。考虑到层间改正的符号，有限范围的层间改正δΔg_limit的算法公式为

 (6)

式中，Ω′表示有限布格层的区域范围，不是任意指定的，需与重力归算中地形改正的范围保持一致，即式(5)中E与式(6)中的Ω′范围相同，均可表达为式(7)

 (7)

式(6)的计算可直接由地形改正中的严格积分公式^[, ]转化而来

 (8)

为进一步简化有限范围的层间改正的计算，本文推导了一种更为简便的计算公式。推导的基本原理是有限范围的层间改正等于无限平面层层间改正减去远区平面虚拟地形(有限范围以外的地形)产生的层间改正，推导过程中多次运用了极限的数学方法。

首先将无限平面布格层分为有限区域和远区两部分，顾及式(2)，得到式(9)

 (9)

式中，Ω-Ω′表示远区范围；Ω′表示有限区域范围，根据式(7)，Ω-Ω′的坐标范围为式(10)

 (10)

将式(9)右端第2项远区层间改正标记为δΔg_Ω-Ω′，标明y轴积分域为式(11)的形式

 (11)

根据积分关系式(12)

 (12)

得到式(11)的积分核对y轴的不定积分

 (13)

当取负极限时，式(13)简化为式(14)

 (14)

当取正极限时，式(13)简化为式(15)

 (15)

将式(14)、式(15)代入式(11)，得

 (16)

利用积分关系式(17)

 (17)

及关系式(18)

 (18)

对式(16)进一步简化

然后利用下列几项极限等式

得到δΔg_Ω-Ω′的计算结果为式(19)

 (19)

综合式(9)、式(19)，得到有限范围的层间改正δΔg_limit的计算公式为

 (20)

对比式(8)、式(20)，得下列等式

 (21)

将式(21)代入式(8)，得到有限范围的层间改正的简单算式(22)

 (22)

一旦给定计算区域的范围，利用式(22)即可快速、严格地计算有限范围的层间改正。

3 试验与分析3.1 算法分析

为了验证式(8)与式(22)的等效性，令计算点的纬度为30°，高度为500 m，并且积分范围(近似为平面)是经纬度等度数的矩形区域。将两种算法计算的有限范围层间改正及式(2)计算的无限平面层层间改正随积分范围的变化表示于，其中横坐标为计算点距离积分区域边缘的度数。

从可以看出，利用地形改正严格积分算法直接演化而来的式(8)与本文推导的简便公式(22)计算的有限范围的层间改正结果的一致性较好，两者的差异随着积分范围的越大而越小，当积分区域范围为20′时，两算法计算的层间改正相差0.11×10^-5 m/s²(即0.11 mGal)。同时，积分范围越大，有限范围的层间改正与传统的无限范围平面层层间改正的差异越小，当积分区域范围为20′时，两算法计算的层间改正与平面层层间改正分别相差0.37×10^-5 m/s²、0.26×10^-5 m/s²。

上文的分析中曾指出传统无限布格层的层间改正中远区虚拟地形产生的误差随计算点高度的增加而增大，为了说明这一误差的大小，下面令计算点的纬度为30°，高度分别取100 m、1000 m、2000 m，将式(22)计算的有限范围的层间改正与式(2)的传统无限范围平面层层间改正的差值表示于。

从可以分析出，有限范围与无限范围层间改正的差异随着计算点高度的增加而增大。当积分区域范围为20′时，对于100 m、1000 m、2000 m高度的计算点，该差异分别为0.01×10^-5 m/s²、1.03×10^-5 m/s²、4.14×10^-5 m/s²。当计算点高程超过1000 m时，有限范围层间改正与传统平面层层间改正相差超过1.0×10^-5 m/s²，此时将有限范围层间改正代替传统的无限范围层间改正是非常有必要的。

3.2 基于重力归算的内插试验

层间改正是重力归算的一项重要内容，而重力归算的一个重要应用方向是离散重力数据内插，因此下面通过分析采用不同层间改正方法的重力归算的内插精度，说明有限范围层间改正的实用意义。试验区范围为107.5°-114.5°E、27.5°-32.5°N，通过地面重力测量获得的地表重力离散点48 390个，重力实测数据的特点是分布不均匀，地势平坦的区域重力点较为密集，精度相对较高，地势起伏较大的区域重力点较为稀疏，精度相对较低。重力离散点高程最大2870 m、最小-34 m，试验区地形分辨率为5″×5″，分布如所示。

首先分别计算球面布格层、平面布格层、有限范围布格层(1.5°、1°，此为实践中通常采用的地形改正的区域范围)的层间改正，将结果统计于。

表 1 不同布格层层间改正Tab. 1 Statistics of Bouguer correction of different Bouguer layers

从可以看出，球面层层间改正的量级约为平面层的2倍，有限范围层间改正与平面层层间改正的数值相差较小，且计算区域的范围越大，差值越小。进一步将有限范围层间改正与传统无限范围平面层层间改正的差值示于。

表 2 有限范围层间改正与平面层差值Tab. 2 The correction difference between the limited Bouguer layer and the plane layer

中1.5°区域范围时二者差异最大值为1.5×10^-5 m/s²，均方根为0.05×10^-5 m/s²。为了说明两者差值大小的具体分布情况，将1°区域范围下的差值分布显示于。

体现了有限范围层间改正与传统平面层层间改正的差异与数据点的高程有关，高程越大该差异也越大。这与前文的分析一致。综合、可以得出，有限范围层间改正与平面层层间改正的差值与重力点的高程、计算范围有关，高度越大、范围越小，此项差异越大。

下面开始内插试验，以通过内插精度的提高说明有限范围层间改正的有效性。首先将试验区48 390个离散重力点筛选出3229个点作为检验点，其余点作为内插点，以检验点内插结果与原数据的差值作为评价内插精度的指标。内插算法为Shepard反距离内插法，重力归算时分别采用不同种层间改正方法，特别强调的是有限范围层间改正的区域范围与地形改正的区域范围一致。为插值统计结果。

表 3 不同种层间改正的插值结果Tab. 3 Interpolation results based on different Bouguer corrections

由可以看出，球面布格层的插值效果很差，不宜内插。表中平面层与有限范围布格层插值均方根相差较小，这是因为3229个检测点多数位于平坦的海拔低的区域，高程低的数据处的平面层层间改正与有限范围层间改正本身的差异较小。当计算范围为1°、1.5°时，有限范围层间改正的内插精度较平面层均方差均提升了3×10^-4×10^-5 m/s²。对于平面层而言，当地形改正范围从1°扩充到1.5°时，均方根提升了6×10^-4×10^-5 m/s²，参照此量级，说明此时有限范围层间改正对内插精度的提高作用不应该被忽略。为了说明有限范围层间改正相对于传统平面层间改正内插的改善效果，将1°区域范围下两者差异显示于。

说明了有限范围层间改正较传统无限范围平面层层间改正对内插精度的提高主要在海拔较高的计算点上，对于地形高度为2000 m的计算点而言，1°区域范围下有限范围层间改正内插精度将提高约1.5×10^-5 m/s²。结合上文的分析，在计算点地形高度超过1000 m时，应使用基于有限范围层间改正的重力归算进行内插。

4 结论

本文在传统平面层、球面层层间改正的基础上，提出了使用有限范围层间改正进行重力归算的方法，其范围与重力归算中地形改正的范围一致。然后通过推导给出了有限范围层间改正的简便计算方法，该式具有一定的数学理论价值。通过试验证明该算法与通过地形改正严密积分法演化来的算法具有一致性。最后通过基于重力归算的内插试验，说明球面层层间改正不适宜内插；当计算点地形高大于1000 m时，应使用基于有限范围层间改正的重力归算进行内插。

【引文格式】马健，魏子卿，武丽丽，等。 有限范围的重力层间改正算法[J]. 测绘学报，2017，46(1)：26-33. DOI: 10.11947/j.AGCS.2017.20160173 

更多精彩内容：

知识|地图先问钓鱼岛：钓鱼岛地图、地理信息与《钓鱼岛图志》

博士论文摘要|巩翼龙：地表覆盖与建筑物目标变化检测技术研究

书讯 | 杨元喜院士：自适应动态导航定位（第二版）

行业动态 | 2017中国航天三大任务世界瞩目：“北斗”全球布局

博士论文摘要|魏钜杰：极化SAR影像海面船只目标自动检测方法研究

博士论文摘要|陈秋杰：基于改进短弧积分法的GRACE重力反演理论、方法及应用

院士论坛|杨元喜：我国海洋大地测量基准与海洋导航技术研究进展与展望

论文推荐|姚宜斌：2011年电离层和太阳活动指数的准21.5天振荡分析

论文推荐|曾安敏：顾及框架点坐标误差的三维基准转换严密模型