基于拉格朗日中值定理的航空重力异常向下延拓方法研究

胡 起 徐新禹 赵永奇 丁子凡

1 武汉大学测绘学院，武汉市珞喻路129号，430079

重力异常向下延拓是物理大地测量的经典问题。长期以来，国内外开展了很多研究来解决向下延拓中的不适定问题，其中以求解球外Dirichlet问题的逆Poisson积分法应用最为广泛。迭代法是解算逆泊松积分的常用方法。Martine与Huang[1-2]利用迭代方式将加拿大落基山脉的航空重力值进行向下延拓；Kingdon、Ma与Zeng等[3-5]对迭代法进行了更深入的研究。正则化方法也是解决病态问题的常用方法。王兴涛等[6]采用Tikhonov正则化方法对5′×5′的航空重力异常进行向下延拓。蒋涛等[7]则将岭估计与广义岭估计应用于航空重力向下延拓。赵启龙等[8]使用半参数估计与最小二乘配置方法对台湾地区的航空重力进行向下延拓。在频率域进行重力延拓也是常用的方法。将重力异常从空间域利用FFT算法转换到频率域，在频率域进行向下延拓，最后将延拓后的重力异常从频率域转换回空间域[9]。宁津生等[10]使用小波分析研究了多尺度边缘约束的重力场向下延拓。黄谟涛等[11]提出基于向上延拓的向下解析延拓方法，利用向上延拓与垂直梯度进行航空重力向下延拓。Zhang也提出利用积分中值定理进行航空重力向下延拓的数值方法[12]。不同于逆Poisson积分法，还有一些方法是利用物理信息进行重力场延拓的研究，如直接代表法、虚拟点质量法、联合重力场模型与地形进行重力场向下延拓等。

本文结合拉格朗日中值定理以及重力异常垂直梯度在垂向上的分布特点，改进使用重力异常垂直梯度进行重力异常延拓的方法，回避了向下延拓的不适定性。针对经典积分求解重力异常垂直梯度的奇异性，本文使用差分的方式计算航空重力高度面上的垂直梯度，并在台湾地区进行模拟与实测航空重力异常向下延拓实验。

1 基本原理与方法

1.1 利用重力异常垂直梯度进行向下延拓

重力异常向下延拓中，最直接的方法是利用垂直梯度进行向下延拓。根据拉格朗日中值定理，有：

(1)

将式(1)变形得到：

(2)

重力异常垂直梯度(Δgr)在垂向上的变化较为缓慢，且呈现出与高度较强的线性关系。为了说明重力异常垂直梯度在垂向的分布特点，本文在台湾地区随机选取100个点，采用2 160阶次的EGM2008重力场模型模拟每个点从海平面到海平面以上10 000 m每隔100 m的重力异常，点位分布见图1。图2为各点位上每隔100 m模拟重力异常垂直梯度与高度的散点图。

图1 Δgr垂向分布实验的点位分布Fig.1 Distribution of points validating the radial pattern of Δgr

图2 各点位上Δgr与高度的散点图Fig.2 Scatter point plots of vertical gradients and height at each test point

由图2可以看出，多数点位上的模拟重力异常垂直梯度与高度呈现出强线性关系。为进一步说明这种线性关系，求得各点位上每隔100 m模拟重力异常垂直梯度与高度的线性相关系数(图3)。在87个点位上，线性相关系数的绝对值大于0.85，重力异常垂直梯度在垂向上与高度呈现出极强的线性关系。

图3 模拟点位上Δgr与高度的线性相关系数Fig.3 Linear correlation coefficients between Δgr and height at each test point

在垂向上，Δgr的变化整体较为缓慢。表1为各个点位上高度7 500 m与5 000 m、7 500 m与2 500 m以及5 000 m与2 500 m的Δgr差值统计。可以看出，三者Δgr差异的绝对值平均值分别为3.3 E、8.1 E、4.8 E。按5 000 m延拓高差计算，由梯度近似造成的相应的误差绝对值平均值分别为1.7 mGal、4.0 mGal以及2.4 mGal，对当前的航空重力观测精度以及向下延拓的精度而言，这个误差在可接受范围之内。

表1 不同高度面Δgr模型值的差异Tab.1 Differences of Δgr on different heights

为了进一步说明前文中重力异常垂直梯度呈强线性分布的结论，采用线性外推计算2 500 m高度处的Δgr：

图4 Δgr线性外推示意图Fig.4 Demonstration of Δgr linear extrapolation

(3)

1.2 重力异常向上延拓与垂直梯度

球面上某一点重力异常的垂直梯度可以用该高度面上的重力异常表示：

(4)

利用重力异常进行重力异常向上延拓的积分公式为：

(5)

式中，P′点为与式(4)中P点同一向径且高出P点H的点。式(5)与式(4)在数学上是等价的，证明如下。

在离散形式下通常使用式(6)进行重力异常向上延拓，以避免中心区域导致积分结果发散[14]:

(6)

在图4中，由中值定理有：

(7)

式中，ΔgP2H可由重力异常向上延拓稳定得到。

本文利用重力异常垂直梯度进行延拓的过程分为3步：

1)根据向下延拓高度H，将航空重力异常向上延拓H，利用式(7)得到航空重力高度面上方某点的垂直梯度；

2)近似计算航空重力高度面的垂直梯度。同样采用向上延拓后差分的方式，首先将重力异常向上延拓1/5H(1/5H为本文差分近似的经验值)，然后差分得到垂直梯度近似值；

3)利用式(3)外推航空重力高度面下方某点的垂直梯度，代入式(2)进行向下延拓。

2 台湾地区航空重力异常向下延拓

本文在台湾地区进行航空重力向下延拓实验。为说明本文方法的有效性与实用性，首先采用EGM2008模型模拟观测值进行模拟实验[15]，然后对实测航空重力数据进行处理。台湾地区有丰富的航空重力观测与地面重力测量数据，本文采用5 156 m高度的航空重力数据进行向下延拓，以及3 360个地面重力观测点进行检核[8]。图5为台湾地区实测航空重力格网值，图6为该地区地面重力值分布及相应重力异常。

图5 台湾地区航空重力格网值Fig.5 Airborne gravity anomaly grid values in Taiwan area

图6 地面重力观测点位分布及观测值大小Fig.6 Ground gravity points and their gravity anomaly

考虑到台湾地区地形复杂，使用GRAVSOFT软件对航空重力观测值及地面重力值进行以15′ 分辨率为平滑半径的残差地形改正[16-17]。考虑到边缘效应，将航空重力在22°～25°N、120°～122°E区域内进行格网化，然后向下延拓，格网分辨率为2′。航空重力异常最大值为253.72 mGal，最小值为-159 mGal。在该区域内有地面重力观测值3 360个。地面重力异常最低点高程为0，最高点高程为3 942.61 m；重力异常最小值为-52.54 mGal，最大值为271.95 mGal，标准差为78.50 mGal。在向下延拓实验中，将航空重力异常分别向下延拓至各地面重力点，再将延拓值与各点重力观测值进行比较。

为了说明本文方法的有效性，首先使用同样位置分布的EGM2008模型的模拟值进行延拓，向下延拓精度见表2。

表2 EGM2008模拟数据向下延拓精度统计Tab.2 Statistics of downward continuation results of EGM2008 simulated gravity anomaly

表3 航空重力向下延拓精度统计Tab.3 Statistics of downward continuation results of airborne gravity observations

3 结 语

本文联合拉格朗日中值定理以及重力异常垂直梯度径向分布的特点，改进了利用重力异常垂直梯度进行航空重力异常向下延拓的方法，利用向上延拓后进行差分的方法获取重力异常垂直梯度，然后外推进行向下延拓。主要结论如下：

1)与常用的泰勒级数一阶展开方式相比，使用拉格朗日中值定理进行垂直梯度外推然后进行向下延拓具有更高的理论严密性。

2)重力异常垂直梯度在台湾地区在海平面到海平面以上10 000 m的空间范围内，在垂向分布上与高度具有强线性关系。在本文进行模拟实验随机选取的100个点位中，87个点位的重力异常垂直梯度与高度的线性相关系数绝对值超过了0.85。

3)在使用EGM2008模拟实验中，直接模拟航空重力高度面上的垂直梯度，利用外推得到航空重力高度面以下的垂直梯度，然后进行向下延拓，精度达到2.3 mGal，较仅利用航空重力高度面垂直梯度或向上延拓差分梯度进行向下延拓精度有较大提升。

4)在差分计算航空重力高度面垂直梯度中，航空重力高度面上垂直梯度计算的不准确会导致利用外推梯度延拓的结果精度反而下降，这是理论严密与数值精度之间的矛盾。将航空重力向上延拓相应点位处的向下延拓高差，将向上延拓前后重力异常与高度进行差分得到差分梯度值，仅利用该差分梯度进行向下延拓，精度达到10.1 mGal，该结果亦优于传统延拓结果[8,19]。这也说明了使用垂直梯度进行向下延拓的优越性。

可以预期，在能得到航空重力高度面上较高精度垂直梯度的情况下，使用本文方法将进一步提升向下延拓的精度。下一步研究将集中在优化利用重力异常计算重力异常垂直梯度的方法上。