重力异常向下延拓是物理大地测量的经典问题。长期以来，国内外开展了很多研究来解决向下延拓中的不适定问题，其中以求解球外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 利用重力异常垂直梯度进行向下延拓

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

$\frac{\Delta g_{P^{\prime}}-\Delta g_{P}}{H}=\left.\frac{\partial \Delta g}{\partial r}\right|_{P^{\prime\prime}} $

(1)

式中，Δg_P′为P′点的重力异常，Δg_P为P点的重力异常，H为P′点与P点间的高差，$\left.\frac{\partial \Delta g}{\partial r}\right|_{P^{\prime \prime}}$为P″点(P″点为P′点与P点间的某一点)处的重力异常垂直梯度。

将式(1)变形得到：

$\Delta g_{P}=\Delta g_{P^{\prime}}-\left.\frac{\partial \Delta g}{\partial r}\right|_{P^{\prime \prime}} H $

(2)

式中，Δg_P′与H为已知量，但是无法准确得到$\left.\frac{\partial \Delta g}{\partial r}\right|_{P^{\prime \prime}}$，原因在于：1)无法确定P″的位置；2)没有理论支持计算已知重力异常所在高度面以下的重力异常垂直梯度。利用式(2)进行向下延拓，需要对$\left.\frac{\partial \Delta g}{\partial r}\right|_{P^{\prime \prime}}$进行近似求解。

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

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

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

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

$\Delta g_{r_{-} 2500}^{e x}=2 \times \Delta g_{r_{-} 5000}-\Delta g_{r_{-} 7500} $

式中，Δg_{r_2500}^ex为2 500 m高度处外推值。Δg_{r_2500}^ex与Δg_{r_2500}的差异统计见表 1的最后一行，此时Δg_{r_2500}^ex-Δg_{r_2500}的绝对值平均值为1.8 E。同样按5 000 m的延拓高差计算，此时由梯度近似造成的误差仅为0.9 mGal。

利用图 4对式(2)中$\left.\frac{\partial \Delta g}{\partial r}\right|_{P^{\prime \prime}}$的近似进行说明。在图 4中，H为延拓高度，P^2H为将航空重力高度面观测值向上延拓H的点位，则P_M^2H为P^2H与P′中的某一点。根据线性外推有：

$\left.\frac{\partial \Delta g}{\partial r}\right|_{P^{\prime \prime}} \approx\left(2 \times\left.\frac{\partial \Delta g}{\partial r}\right|_{P^{\prime}}-\left.\frac{\partial \Delta g}{\partial r}\right|_{P_{M}^{2 H}}\right) $

(3)

式中，$\left.\frac{\partial \Delta g}{\partial r}\right|_{P^{\prime \prime}}$、$\left.\frac{\partial \Delta g}{\partial r}\right|_{P^{\prime}}$、$\left.\frac{\partial \Delta g}{\partial r}\right|_{P^{2H}_M}$分别为P″、P′与P_M^2H点处的重力异常垂直梯度。

$\left.\frac{\partial \Delta g}{\partial r}\right|_{P^{\prime}}$可利用航空重力高度面上的重力异常积分得到，但是在离散形式下通常需要克服奇异性的问题^[13]。本文同样使用向上延拓后差分的方式对航空重力高度面上垂直梯度进行近似。

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

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

$\left.\frac{\partial \Delta g}{\partial r}\right|_{P}=-\frac{1}{R} \sum\limits_{n=0}^{\infty}(n+2) \Delta g_{n}=-\frac{1}{R} \sum\limits_{n=0}^{\infty} n \Delta g_{n}-\\ \frac{2}{R} \Delta g=\frac{R^{2}}{2 \pi} \iint\limits_{\sigma} \frac{\Delta g-\Delta g_{P}}{l_{0}^{3}} \mathrm{~d} \sigma-\frac{2}{R} \Delta g_{P} $

(4)

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

$\Delta g_{P^{\prime}}=\frac{R^{2}}{4 \pi r} \iint\limits_{\sigma} \frac{r^{2}-R^{2}}{l^{3}} \Delta g \mathrm{~d} \sigma $

(5)

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

$\begin{array}{*{20}{c}} \left.\frac{\partial \Delta g}{\partial r}\right|_{P}=\lim\limits_{H \rightarrow 0} \frac{\Delta g_{P^{\prime}}-\Delta g_{P}}{H}=\lim\limits_{H \rightarrow 0} \frac{\frac{R^{2}}{4 \pi r} \iint \frac{r^{2}-R^{2}}{l^{3}} \Delta g \mathrm{~d} \sigma-\frac{R^{2}}{r^{2}} \Delta g_{P}+\frac{R^{2}}{r^{2}} \Delta g_{P}-\Delta g_{P}}{H}=\\ \lim\limits_{H \rightarrow 0} \frac{\frac{R^{2}}{4 \pi r} \iint \frac{r^{2}-R^{2}}{l^{3}}\left(\Delta g-\Delta g_{P}\right) \mathrm{d} \sigma+\left(\frac{R^{2}}{r^{2}}-1\right) \Delta g_{P}}{H}=\\ \lim\limits_{H \rightarrow 0} \frac{\frac{R^{2}}{4 \pi(R+H)} \iint \frac{(R+H)^{2}-R^{2}}{l^{3}}\left(\Delta g-\Delta g_{P}\right) \mathrm{d} \sigma+\left(\frac{R^{2}}{(R+H)^{2}}-1\right) \Delta g_{P}}{H}=\\ \lim\limits_{H \rightarrow 0} \frac{\frac{R^{2}}{4 \pi(R+H)} \iint \frac{2 R H+H^{2}}{l^{3}}\left(\Delta g-\Delta g_{P}\right) \mathrm{d} \sigma-\left(\frac{2 R H+H^{2}}{(R+H)^{2}}\right) \Delta g_{P}}{H}=\\ \lim\limits_{H \rightarrow 0} \frac{R^{2}}{4 \pi(R+H)} \iint \frac{2 R+H}{l^{3}}\left(\Delta g-\Delta g_{P}\right) \mathrm{d} \sigma-\left(\frac{2 R+H}{(R+H)^{2}}\right) \Delta g_{P}=\frac{R^{2}}{2 \pi} \iint \frac{\left(\Delta g-\Delta g_{P}\right)}{l^{3}} \mathrm{~d} \sigma-\frac{2}{R} \Delta g_{P} \end{array} $

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

$\Delta g_{P^{\prime}}=\frac{R}{r} \Delta g_{P}+\frac{R}{4 \pi} \iint\limits_{\sigma} \frac{r^{2}-R^{2}}{l^{3}}\left(\frac{R}{r} \Delta g-\Delta g_{P}\right) \mathrm{d} \sigma $

(6)

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

$\left.\frac{\partial \Delta g}{\partial r}\right|_{P_{M}^{2 H}}=\frac{\left(\Delta g_{P^{2 H}}-\Delta g_{P^{\prime}}\right)}{H} $

(7)

式中，Δg_P^2H可由重力异常向上延拓稳定得到。

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

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

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

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

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

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

考虑到台湾地区地形复杂，使用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。

在模拟实验中，仅利用向上延拓H后差分得到的梯度进行向下延拓的残差标准差为6.3 mGal。值得注意的是，该模拟实验的$\left.\frac{\partial \Delta g}{\partial r}\right|_{P^{\prime}}$由EGM2008模型模拟得到，直接使用该模拟梯度进行向下延拓，残差标准差为3.2 mGal；而在进行梯度外推后向下延拓，残差标准差提升到2.3 mGal。证明了本文利用梯度进行向下延拓的理论有效性，同时也提供了§1.1中重力异常梯度在垂向上强线性分布的间接证据。

上述结果中的$\left.\frac{\partial \Delta g}{\partial r}\right|_{P^{\prime}}$由模型得到，在实际应用中往往并不具备该观测值，而需要根据重力异常观测值进行计算。使用§1.2中的步骤再次对模拟数据进行延拓，同时对实际观测值进行处理。在对航空重力观测值进行向上延拓的过程中，采用200阶的GO_CONS_GCF_2_TIM_R5重力场模型进行移去-恢复以提升向上延拓的精度^[18]。延拓结果统计见表 3。

表 3中，第1行的结果为使用模型值进行延拓的统计结果，第2行为实测航空重力向下延拓的结果。与表 2的模拟数据延拓结果相比较可以发现，仅使用$\left.\frac{\partial \Delta g}{\partial r}\right|_{P^{\prime}}$延拓的残差标准差由3.2 mGal增大到6.7 mGal，原因在于：积分公式在离散形式下在较低延拓高度时存在较高的奇异性，$\left.\frac{\partial \Delta g}{\partial r}\right|_{P^{\prime}}$的计算精度不高，而使用差分近似同样存在精度损失。$\left.\frac{\partial \Delta g}{\partial r}\right|_{P^{\prime \prime}}$由$\left.\frac{\partial \Delta g}{\partial r}\right|_{P_{M}^{2 H}}$与$\left.\frac{\partial \Delta g}{\partial r}\right|_{P^{\prime}}$外推得到，所以$\left.\frac{\partial \Delta g}{\partial r}\right|_{P^{\prime}}$的误差同时会影响到使用外推梯度进行延拓的结果。在此种情形下，梯度计算误差影响超过外推近似，得到精度提升，所以仅使用$\left.\frac{\partial \Delta g}{\partial r}\right|_{P^{2H}_M}$梯度进行延拓的精度反而最高。在实测航空重力的实验中有同样的精度变化，仅使用$\left.\frac{\partial \Delta g}{\partial r}\right|_{P^{2H}_M}$的残差标准差为10.1 mGal, 仅使用$\left.\frac{\partial \Delta g}{\partial r}\right|_{P^{\prime}}$的残差标准差为11.8 mGal，而使用外推后的梯度进行延拓的残差标准差为15.4 mGal，其中$\left.\frac{\partial \Delta g}{\partial r}\right|_{P^{\prime}}$的计算误差影响更大。

特别指出的是，使用外推梯度进行延拓具有更高的理论严密性，理应获得更高的延拓精度，这一点在表 2的结果中得到印证。而在实际计算中，计算误差的影响超过了理论严密的作用，仅使用$\left.\frac{\partial \Delta g}{\partial r}\right|_{P^{2H}_M}$得到了最佳的精度，这是在实际应用中的妥协，也就导致了表 3的结果。如果能得到航空重力高度面上的高精度垂直梯度值，比如在航空重力观测时同时进行重力梯度观测，或改进梯度计算方法得到更高精度的$\left.\frac{\partial \Delta g}{\partial r}\right|_{P^{\prime}}$，可以预期，使用外推梯度进行向下延拓结果更好。

3 结语

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

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

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

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

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

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