重力梯度之地形改正方法与研究

冯进凯，王庆宾，黄佳喜，张 超

(信息工程大学，河南 郑州 450000)

引用著录：冯进凯，王庆宾，黄佳喜,等.重力梯度之地形改正方法与研究[J].测绘工程，2017，26(12)：50-54.

DOI:10.19349/j.cnki.issn1006-7949.2017.12.010

重力梯度之地形改正方法与研究

冯进凯，王庆宾，黄佳喜，张 超

(信息工程大学，河南 郑州 450000)

简要介绍地形起伏对地面点重力梯度的影响，推导棱柱法计算重力梯度改正的严格积分公式并和积分核累加法对比验证了其正确性；将积分区域划分为中央区和远区，中央区采用严格的棱柱法进行计算，远区采用线性积分的方式逼近真实值，简化了积分过程；根据重力梯度应用的要求针对不同类型地形设定了中央区划分准则。实验表明：根据文中提出的方法和分区准则，相比于传统的积分法和棱柱法，计算速度得到提高，且精度要求相比于积分核累加法优于0.01 E。

重力梯度；地形改正；分区准则

重力梯度是重力沿某一方向的变化率，它的空间分布和变化与地下构造、分布密度有着很强的联系[1]，广泛用于矿产探测、地下军事目标识别和重力梯度匹配导航等领域[2]。特别在航空重力梯度测量中，需要高精度的地面控制点进行控制和联合平差。高精度高分辨率的地形数据为充分利用地形信息进行重力梯度改正提供了可能。在完成大范围的梯度地形改正过程中，过去通常采用直接积分法[3]，即将地形格网进行分割，然后对积分核累加求和，这种方法耗时严重，难以在短时间内完成。例如：使用30″×30″的地形数据，积分范围取2°×2°，在CD-4360小型计算机上完成百万图幅格网点梯度地形改正所需时间超过400 h[1]。近些年来，国内外学者将快速Fourier变换(FFT)应用其中。Tziavos利用FFT算法进行了梯度地形改正的计算[3]；Kass也推进了这一领域的研究[4]，使得计算速度成百倍的提高，但是精度略有损失，特别是在计算机性能高速发展的今天，其优势不太明显；边少锋等将圆柱体模型和棱柱模型应用于梯度改正，确定了梯度垂直分量Tzz的值[5]； 武凛等也进行了相应的研究[6]；LiZhi Zhu等将重力梯度正演模型计算方法进行了归类总结和比较[7]。

本文通过分析不同DEM数据特点的基础上，将积分区域划分为中央区和远区，并分别利用不同的方法进行计算。中央区域采用棱柱积分的方法；远区采用线性卷积分的方法来逼近真实值。通过分析实验结果并针对不同地形数据提出了中央区划分准则，并在此基础上进行梯度地形改正计算。实验结果表明：此种改进方法可以使得计算速度大幅提高，且精度与直接积分法相当。

1 地面点重力梯度局部地形改正原理

1.1 传统重力梯度地形改正的算法

梯度地形改正就是地面点(计算点)周围地形起伏对该点梯度的影响[6]，以计算点为原点，建立局部笛卡尔坐标系，x轴正方向指向北，y轴正方向指向东，为表示方便x1=x，x2=y，x3=z。图1所示在计算点P(ip,jp)处，重力梯度的地形改正可以表示为

z.

(1)

在计算中，通常将其化为累加形式，积分核累加形式如下：

(2)

其中:N，M分别为积分区域中积分单元的行数、列数，Q为高度方向划分的单元格数。

图1 局部地形重力梯度改正

1.2 改进1：棱柱法求解梯度分量

根据式(1)可得梯度三方向的分量如下[6]：

(3)

现将式(3)中的函数直接积分求得原函数，以格网的边界为上下限，每个地形格网只计算一次，从而避免大量的累加运算。

地形网格数据步长近似为Δx=RΔB，Δy=RcosBmΔL，Bm为积分区域的平均纬度。ΔL和ΔB分别为球面网格的分辨率。选择积分区域为N×M，N为行数，M为列数,为避免积分中的强奇异性问题，在此将计算点选在格网中心。

对式(3)积分后去离散求和形式如下：

(4)

其中：(i,j)为计算点编号，(ip,jp)为积分点编号，至此将复杂的三重积分化为了简单形式，即在一定分辨率的地形数据下，每一个地形格网数据只需计算一次，避免了式(2)中复杂的累加运算，将原来的N×M×Q次累加运算降低为N×M运算。

1.3 改进2：划分近远区求解梯度分量

在上述计算中，所有的距离量都是相对于积分点而言的，在距离计算点较远区域，h(i,j)-h(ip,jp)相对于距离量来说为小值。为加快计算速度，将

积分区σ划分为中央区σ0和远区σ1=σ-σ0， 最后结果为中央区与远区地形改正之和：

(5)

(6)

2 数值计算分析

为研究不同地形数据对中央区划分准则的影响以及划分准则对计算速度的影响，本文搜集了3块区域(1°×1°)，3块区域的高程统计特性见表1。

本次采集数据为DEM数据，分辨率为3″×3″，选取了3块地形特征明显的区域。平原地区地形较为平缓，高程较低；山地地区地形落差大，高程普遍较高，地形起伏明显。

表1 3块区域高程统计信息

2.1 方法正确性检验

积分核累加法是严格的数值计算方法，其正确性已经得到广泛检验。考虑到积分核累加法对时间的要求，本实验采取固定积分区域。本文提出算法直接使用中央区积分法。在积分核累加法中，将地形网格划分为底面积1 m×1 m，高度为1 m的正方体，计算点区域为3′×3′，使用两种方法计算得到差值见表2。

表2 两种方法计算结果对比

根据式(3)严格的数学推导和上述实验两种方法的差值比较，本文方法的正确性得到检验。

2.2 中央区划分准则

通过数值初步估算，距离计算点20 km以外积分点对计算点的梯度影响几乎为0，因此固定积分半径为20 km，密度为2.67 g/cm3，中央区域采用棱柱法积分，远区采用简化积分，通过改变中央积分区域大小，得到在3种地形数据下，重力梯度3个方向的值随中央区半径的变化。

图2 区域1X 方向分量随中央区的变化

图3 区域1Y 方向分量随中央区的变化

图4 区域1Z 方向分量随中央区的变化

图5 区域2X 方向分量随中央区的变化

图6 区域2Y 方向分量随中央区的变化

图7 区域2Z 方向分量随中央区的变化

图8 区域3X 方向分量随中央区的变化

图9 区域3Y 方向分量随中央区的变化

图10 区域3Z 方向分量随中央区的变化

图2—图4、图5—图7、图8—图10，分别对应区域1(丘陵)、区域2(平原)、区域3(山地)的重力梯度三方向的改正值。实验数据表明：当中央区积分半径过小时，不能将高差看做距离的小值，所以计算结果存在一定的误差。为突出最后收敛结果特征，将中央区半径从1 km开始计算。

从图中分析可得：随着中央区的增大，重力梯度3个方向分量值呈收敛状态，趋近真实值；当积分半径较小时，可以清楚看出，地形对梯度值的影响较为剧烈，而且不同地形的收敛节点不同，因此设定不同地形的不同中央区准则是合理的；每一组三方向分量的计算值满足拉普拉斯方程，即Txx+Tyy+Tzz=0，从侧面证明了计算结果的正确性；三区域梯度值对比可以得出平原地区梯度值的收敛速度最快，丘陵地区次之，山区收敛速度最慢；要满足梯度控制点和矿产探测的精度要求(±0.01 E)，现提出平原地区中央区积分半径不小于5 km，丘陵地区中央区积分半径不小于7 km，山地区域中央区积分半径不小于9 km的中央区划分准则。这样的划分准则，既满足了点位的精度要求，又在一定程度上加快了计算速度；通过对比以上数据，不难发现，与地形高低起伏最为相关的是重力的垂直方向的绝对值，垂直梯度反映了重力在径向的变化速率，因此有正有负。

2.3 基于以上算法和准则的大范围计算点实验

为保证积分法能在有限时间内计算完成，所以本次试验区地形数据分辨率为1′×1′，选取3块不同地形区域(2°×3°)，现通过积分核累加法和本文提出算法进行计算，二者总积分区域半径选为20 km，积分核累加法中将累加单元设为底面积边长为300 m，高为1 m的正方体，密度为2.67 g/cm3。本文提出算法中，中央区化半径分别为三块区域设定的划分准则(平原5 km，丘陵7 km，山地9 km)，计算所用时间对比如表3所示。

表3 不同计算方法所用时间对比

积分核累加法虽然形式简单且具有严格的数学意义，但是小单元累加中，长方体划分越小，累加次数越多，计算速度越慢，山区尤为明显。而本文介绍的棱柱法，每个网格数据中只需计算一次，采用本文介绍的方法，使得计算时间大大缩短，而且基于上述提出的分区准则，计算时间还会进一步缩减。

3 结 论

本文主要针对地形起伏对地表点的重力梯度的影响，采用不同方法进行试验，得到如下结论：

1)在以计算点为原点的局部坐标系中推导了梯度地形改正积分公式，大大简化了计算流程，且试验验证了该方法的正确性，两种方法差值均值绝对值不超过±0.01 E，且计算速度呈数十倍的提高。

2)推导了远区的积分公式，即在远区用线性卷积分的方法去逼近近似值，实验表明，采用该种分区方法计算时间会进一步缩短。

3)提出基于不同地形条件下的中央区分区准则，对以后采用更高分辨率地形数据细化分区准则有一定的指导意义。

[1] YULE D E. Microgravity investigation of foundation conditions[J]. Geophysics, 1998, 63(1).

[2] 张赤军.用地形数据确定重力异常垂直梯度[J].科学通报,1999,44(6)：656-661.

[3] TZIAVO I N. The effect of the terrain on airborne gravity and gradiometry [J].Journal of Geophysical Research,1988,93 (B8):173-186.

[4] KASS M A, LI Y. Practical aspects of terrain correction in airborne gravity gradiometry surveys[J]. Exploration Geophysics，2008, 39: 198-203.

[5] 边少锋,张赤军.地面扰动重力垂直梯度的确定[J].地球物理学进展,2005,20(4):969-973.

[6] 武凛,马杰,周瑶,等.重力匹配导航的全张量重力梯度基准图模拟[J].系统仿真学报,2009,21(22):7037-7041.

[7] CHRISTOPHER J, ZHU L Z. Comparison of methods to model the gravitational gradients from topographic data bases[J]. Geophys.J.Int,2006,166:999-1014.

[8] DAVID H．Evaluation of a full tensor gravity gradiometer for kimberlite exploration[A]. Airborne Gravity 2014一 Abstracts from the ASEG-PESA Airborne Gravity 2014 Workshop: Geoscience Australia Record[C]. 2014,18:73-79.

[9] 边少锋,张赤军.地形起伏对重力垂直梯度影响的计算[J]. 物探化探计算技术,1999, 21(2): 133-140.

[10] PAJOT G, DE V O,DIAMENT M, et al. LequentrecLalaneette4, and V. Mikhailov. Noise reduction through joint processing of gravity and gravity gradient data[J]. Geophysics, 2008, 73(3):123-134.

[责任编辑：刘文霞]

Methodsofterraincorrectionongravitygradiometry

FENG Jinkai, WANG Qingbin, HUANG Jiaxi, ZHANG Chao

(Information Engineering University,Zhengzhou 45000,China)

This paper presents the influence of gravity gradient caused by topographic relief on ground points briefly. Strict integral formula of gradient correction is deduced in Performance Prism, whose correctness is verified compared with Integral Nuclear Accumulation Method. The integration zone is divided in two parts, one of which is the center zone, and the other is the outside zone. The integration process is simplified using linear integral approximation to the real values in outside zone and time is shorten using the rules of partition for different types of terrain. Experiments show that computing time is shorter and accuracy is better than the 0.01 E compared with Integral Nuclear Accumulation Method based on the method and the partition criterion proposed in this paper.

gravity gradient；terrain correction；partition

P631

A

1006-7949(2017)12-0050-05

2016-10-22

冯进凯(1992-)，男，硕士研究生.