0 引言

重力归一化总梯度法由前苏联学者Березкин^[1]于1967年提出，因其先对重力异常进行向下解析延拓，再计算归一化总梯度而得名。它是一种利用高精度重力异常确定场源、断裂位置及密度界面的方法。肖一鸣^[2]首次将该方法引入中国，给出了利用泰勒级数展开计算重力归一化总梯度的方法，并探讨了影响该方法应用效果的因素，如谐波数、随机噪声、测线长度等。随后该方法在中国得到了广泛应用。肖一鸣等^[3]将该方法成功应用于油气勘探；王家林等^[4]利用该方法分析断层和密度分界面；吴燕冈等^[5]提出将重力归一化总梯度与相位叠置处理，用于确定深大断裂的空间位置；张凤旭等^[6-7]利用Hilbert变换改进重力归一化总梯度法，并提出在位场转换的向下延拓和求导过程中分别引入圆滑滤波因子，提高了该方法的分辨率、可靠性和稳定性，增大了延拓深度；肖鹏飞等^[8]利用向下延拓的迭代算法替代波数域的直接向下延拓算子，提高了重力归一化总梯度法的稳定性；张凤琴等^[9]采用DCT变换计算重力归一化总梯度，增加了下延深度；王彦国等^[10]和郭灿灿等^[11]提出了基于泰勒级数迭代法的快速稳定向下延拓的重力归一化总梯度法；苏超等^[12]对归一化函数的分母计算几何平均，提高了重力归一化总梯度法的抗噪能力；王选平等^[13]将正则化因子引入重力场的下延计算，提出了利用正则化方法计算重力归一化总梯度的方法；王彦国等^[14]提出基于幂次平均的离散归一化总梯度法，可以有效识别叠加场源的位置信息。由此可见，对重力归一化总梯度法的改进主要集中在向下延拓和求导这两步计算过程，其稳定性和精度直接影响重力归一化总梯度法的稳定性和精度。

张冲等^[15-17]针对向下延拓问题，提出了基于求解微分方程的Adams-Bashforth公式法和基于Simpson求积公式的Milne法。前者下延深度能达到5倍点距，后者甚至可超过10倍点距。Milne法利用观测面上的重力垂向一阶导数、向上延拓不同高度的重力场和重力垂向一阶导数求解下延深度的重力场。向下延拓Milne法具有下延深度大、相对误差小的特点，能获得稳定准确的结果。但该方法需要计算垂向导数，且求导的精度直接影响向下延拓的精度。

早在2001年，Fedi等^[18]基于重力位的二阶导数在场源外满足拉普拉斯方程的原理，提出了对随机噪声有压制作用的积分垂向二阶导数法(Integrated Second Vertical Derivative, ISVD)。之后，崔瑞华等^[19]将该方法引入中国，首先对重力异常积分，计算得到重力位，然后计算重力位的水平二阶导数，最后根据拉普拉斯方程计算重力异常的垂向一阶导数。相对于常规的波数域求垂向导数的方法，ISVD法由于进行了积分运算，对随机噪声有抑制作用，因而更稳定。

本文将这两种稳定算法，即向下延拓Milne法和ISVD法引入重力归一化总梯度的计算，通过理论模型试算和实际资料处理，验证了方法的有效性。

1 方法原理 1.1 重力归一化总梯度法的基本原理

以重力剖面观测数据为例，归一化总梯度的表达式^[2]为

$ {G^{\rm{H}}}(x, z) = \frac{{\sqrt {V_{zx}^2(x, z) + V_{zz}^2(x, z)} }}{{\frac{1}{M}\sum\limits_{i = 1}^M {\sqrt {V_{zx}^2\left( {{x_i}, z} \right) + V_{zz}^2\left( {{x_i}, z} \right)} } }} $

(1)

式中：G^H(x, z)为计算点(x, z)处的重力归一化总梯度值；M为重力剖面上的测点总数；V_zx(x, z)和V_zz(x, z)为计算点的重力位二阶偏导数，同时也是计算点(x, z)处重力异常沿x方向和z方向的一阶偏导数。

1.2 向下延拓Milne法的基本原理

重力场向下延拓四阶Milne公式^[15]为

$ \begin{array}{c}{V_{z}\left(x, z_{h}\right)=V_{z}\left(x, z_{-3 h}\right)+\frac{4 h}{3}\left[2 V_{z z}\left(x, z_{-2 h}\right)-\right.} \\ {\left.V_{z z}\left(x, z_{-h}\right)+2 V_{z z}\left(x, z_{0}\right)\right]}\end{array} $

(2)

式中：h为向下延拓深度；V_z(x, z_h)表示由观测面z₀向下延拓h深度的重力场；V_z(x, z_-3h)表示由观测面向上延拓3h高度的重力场；V_zz(x, z_-h)和V_zz(x, z_-2h)分别表示由观测面向上延拓h和2h高度的重力垂向一阶导数；V_zz(x, z₀)表示观测面上的重力垂向一阶导数。由式(2)可知，计算下延重力场时需要利用向上延拓所得到的重力场及其垂向一阶导数。

这里选择在波数域中进行向上延拓。向上延拓是一种稳定算法，对高波数成分(包括噪声)有压制作用。将观测面上的重力异常向上延拓nh高度的表达式为

$ \widetilde{V}_{z}\left(u, z_{-n h}\right)=\widetilde{V}_{z}\left(u, z_{0}\right) \mathrm{e}^{-2 \pi n h \sqrt{u^{2}}} $

(3)

式中：$\widetilde{V}_{z}\left(u, z_{-n h}\right)$和$\widetilde{V}_{z}\left(u, z_{0}\right)$分别表示由观测面向上延拓高度为nh的重力场的波谱和观测面上的重力场的波谱；$\mathrm{e}^{-2 \pi n h \sqrt{u^{2}}}$为波数域向上延拓nh高度的转换因子；u为x方向的波数。

式(2)右端涉及向上延拓和垂向导数两个计算过程，其中向上延拓具有压制随机噪声的作用，而垂向求导有提高分辨率的作用，因此该算法比较稳定，也具有较高的精度。根据张冲等^[15]的理论模型试验结果，下延超过5倍以上点距深度时，Milne法的计算误差显著低于积分迭代法和传统波数域向下延拓方法，即使下延10倍点距时，下延结果也不会发散。

式(2)中的V_zz(x, z_-2h)和V_zz(x, z_-h)是重力场经过向上延拓后再计算其垂向一阶导数，因而对随机噪声不敏感；V_zz(x, z₀)是观测面的重力垂向一阶导数，一般情况下，它是由观测面上的重力异常直接求垂向导数所得，因此受随机噪声影响最大。故本文选择对噪声有一定压制作用的垂向导数算法，即ISVD法。

1.3 ISVD法的基本原理

二度体引起的重力位二阶导数，在场源外满足拉普拉斯方程

$ V_{x x}+V_{x}=0 $

(4)

式中：V_xx表示重力位的x方向二阶导数；V_zz表示重力位的垂向二阶导数。通常实测数据是重力异常，为了获得重力位，可先在波数域对重力异常V_z积分得到重力位V

$ \widetilde{V}=\frac{1}{2 \pi \sqrt{u^{2}}} \widetilde{V}_{z} $

(5)

式中：$\widetilde{V}$表示重力位的波谱；$\widetilde{V}_{z}$表示重力异常的波谱；$\frac{1}{2 \pi \sqrt{u^{2}}}$为波数域中的积分算子。然后在空间域计算V_xx。为了对重力数据中的随机噪声进行压制，用三点滑动窗口内的二次曲线拟合重力异常，再求取其水平导数^[20-21]。空间域水平一阶导数的公式为

$ V_{x(0)}=\frac{1}{2 \Delta x}\left(V_{(1)}-V_{(-1)}\right) $

(6)

式中：V_x(0)表示三点滑动窗口内中心点重力位的水平一阶导数；V₍₁₎和V_(-1)分别表示滑动窗口右端点和左端点的重力位；Δx表示点距。将第一次使用式(6)计算得到的V_x再次代入式(6)，可求得重力位V在x方向的二阶导数V_xx。最后将V_xx代入式(4)即得到重力异常垂向一阶导数V_zz。由于该方法先进行积分运算，在一定程度上平滑了重力异常中的随机噪声。值得注意的是，该方法并不能完全抑制数据中的非随机噪声，例如浅层地质体引起的干扰。当重力异常中含有“毛刺”等明显干扰时，建议先去除高波数的干扰成分，然后再做后续处理。

1.4 基于向下延拓Milne法的重力归一化总梯度法的实现步骤

① 以向下延拓深度h为例，利用式(3)在波数域对观测剖面上重力异常的波谱$\widetilde{V}_{z}\left(u, z_{0}\right)$)进行不同高度(h、2h、3h)的向上延拓，得到剖面以上不同高度重力异常的波谱$\widetilde{V}_{z}\left(x, z_{-h}\right), \widetilde{V}_{z}\left(x, z_{-2 h}\right), \widetilde{V}_{z}(x, \left. {{z_{ - 3h}}} \right)$。

② 根据式(5)，在波数域对$\widetilde{V}_{z}\left(u, z_{0}\right), \widetilde{V}_{z}(u$, ${z_{- h}}), {\tilde V_z}(u, {z_{ - 2h}})$进行积分，获得重力位的波谱$\widetilde{V}\left(u, z_{0}\right), \widetilde{V}\left(u, z_{-h}\right), \widetilde{V}\left(u, z_{-2 h}\right)$。两次应用式(6)，在空间域分别计算其重力位沿x方向的二阶导数。再将结果代入式(4)，依次得到重力位的垂向二阶导数V_zz(x, z₀)、V_zz(x, z_-h)、V_zz(x, z_-2h)。

③ 将第①步和第②步计算得到的结果代入式(2)，求出下延深度为h的重力异常V_z(x, z_h)。

④ 再利用第②步方法，计算V_z(x, z_h)的垂向一阶导数V_zz(x, z_h)，并应用式(6)直接在空间域计算V_z(x, z_h)沿x方向的一阶导数V_zx(x, z_h)。

⑤ 将V_zz(x, z_h)和V_zx(x, z_h)代入式(1)计算G^H(x, z_h)。

⑥ 修改下延深度，重复以上步骤，获得不同深度的重力归一化总梯度值。

⑦ 最后绘制G^H(x, z)等值线剖面，极大值点即是异常体的中心位置。

2 理论模型实验

设计一个无限长水平圆柱体模型，剩余密度为-0.5×10³kg/m³，中心位置为x=100m、z=25m，观测剖面长度为200m，测量点距为1m，共201个观测点。计算该模型的重力异常，分别利用本文方法和基于泰勒级数展开的重力归一化总梯度法计算其中心位置，结果如图 1所示。

由图 1可知，采用两种不同的重力归一化总梯度法，其极大值点(图 1b和图 1c中的红色圆点)均能准确地指示出水平圆柱体模型的中心位置，说明这两种方法处理理论模型数据的效果较好。

实测重力数据中通常包含区域场和随机误差。为检验该方法的抗干扰能力，对模型数据加入5%的随机噪声和一阶趋势背景场，处理结果如图 2所示。

由图 2可知，基于向下延拓Milne法的重力归一化总梯度法所得到的模型中心位置为x=100m、z=28m(图 2b)，与理论坐标位置相比，计算结果偏深3m，水平位置无偏移，说明随机噪声和背景场仍然对该方法有一定的影响。同样，基于泰勒级数展开的重力归一化总梯度法也受影响，计算埋深偏大1m，水平位置向左偏移2m(图 2c)。对比图 2b与图 2c可见，后者出现很多等值线圈闭。在实际资料处理中，这些虚假异常圈闭容易引起对真实场源位置的误判。引起这些圈闭的原因是基于泰勒级数展开的重力归一化总梯度法在计算过程中放大了高波数信号。相比而言，基于向下延拓Milne法的重力归一化总梯度法在一定程度上压制了高波数噪声，计算结果稳定，虚假异常圈闭相对较少。

选取图 2a剖面两端的数据(0~20m及180~200m)进行一阶多项式拟合，用拟合值作为区域异常；再用重力异常减去区域异常得到局部异常，并对其进行圆滑滤波处理，得到图 3a中虚线所示的异常曲线；最后分别用基于向下延拓Milne法的重力归一化总梯度法和基于泰勒级数展开的重力归一化总梯度法得到图 3b和图 3c所示结果。由图 3b可见，计算得到的模型中心位置为x=98m、z=25m，中心埋深与理论模型一致，水平方向有2m的左向偏移；由图 3c可见，计算得到的模型中心位置为x=98m、z=26m，中心埋深比实位置偏大1m，水平方向向左偏移2m。这两种方法均出现水平位置的偏差，分析应该与区域背景场分离不彻底有关，因为分离后的局部重力异常的极小值也有微小的左向偏移。

由模型试验可知，利用本文方法计算重力归一化总梯度具有中心埋深计算准确、虚假异常圈闭相对较少等优点；但若应用于含噪、含背景场干扰的重力数据处理时，需要进行背景场的剥离和去噪处理，这样才能得出比较准确的场源中心位置。

3 实测重力剖面数据应用

将本文方法应用于辽宁省葫芦岛市杨家杖子镇经济技术开发区黑鱼沟A矿洞的实测重力剖面数据。A矿洞位于黑鱼沟西山寒武系下方，是铅锌矿开采矿洞。矿洞断面为马蹄形，水平延伸近50m，可近似为一个无限长的水平圆柱体。矿洞中心埋深为5m，高和宽均约2m。在地面采集重力数据，重力观测剖面垂直于矿洞走向。

观测布格重力异常如图 4a中的实线所示。首先对该数据进行异常分离，这里采用非线性滤波分离区域(背景)场(图 4a中虚线)；去掉背景场后的局部重力异常如图 4b中实线所示，可见局部重力异常含有明显的高波数干扰成分，故对其进行去噪处理。选用多次圆滑方法压制随机噪声，去噪后的结果如图 4b中的虚线所示，对其分别利用本文方法和基于泰勒级数展开得到如图 4c和图 4d所示的重力归一化总梯度剖面。对比图 4c与图 4d可见，根据本文方法计算结果解释的矿洞深度与实际情况基本一致，而根据图 4d解释的矿洞埋深较真实位置偏大1.5m。对比其他方法计算的该矿洞中心埋深(5.0~5.6m)^[22-25]可以证实，本文方法处理结果的准确度较高。

4 结论与建议

理论模型和实际数据计算结果表明，基于向下延拓Milne法的重力归一化总梯度法具有算法稳定、准确性较高的优点。相对于其他重力归一化总梯度法，不需要设置最佳谐波数、正则化参数及迭代次数等参数，减少了数据处理中参数选取引入的误差。但原始重力数据中的随机噪声和背景场对本方法有一定的影响，故在实际应用中需首先进行区域场的分离和去噪处理。