0 引言

重力场分离是重力数据反演、地质解释与推断的前提^[1-3]，是业内公认的重点研究课题和技术难题。早在1952年，Dobrin等^[4]就提出重力解释中两个最重要的问题：一个是区域场与局部场的分离；另一个是基于分离得到的异常获取目标地质体的空间位置和密度信息。程方道等^[5]指出，与磁异常相比，重力异常普遍存在强大的区域背景场，往往会造成局部异常发生畸变，甚至难于辨识，要对重力资料进行有效的地质解释，首先必须分离重力区域场与局部场；管志宁等^[6]强调，正确划分不同规模的区域场与局部场是面积性测量资料解释能否获得良好地质效果的一个重要因素；曾华霖^[7]认为，在密度界面反演中，只要平均深度取值合理，界面之间的密度差异比较接近实际情况，无论采用何种反演方法，都能得到达到一定的精度及相似的结果，但前提是用于反演的异常必须单纯由目标体引起；徐世浙等^[8]认为重力场三维反演是个难度极大的课题，提出异于常规解决方案的开拓性思路，主要包括三个方面：对平面上的重力场进行不同深度的分离；用大深度位场向下延拓的迭代法将各个切割层在地面的重力场下延至相应深度；将顶面的重力场反演为视密度。现今，位场大深度的向下延拓技术和层源异常的反演技术已经得到了很好发展，位场分离方法是目前影响反演方法效果的关键技术，继续加强对异常分离方法的研究具有重要意义。

重力场分离方法大致可分为空间域方法和频率域(波数域)方法。频率域方法主要包括维纳滤波法^[9]、匹配滤波法^[10-12]、小波多尺度分析法^[13-17]、优选滤波法^[18-19]和优化滤波法^[20]等，这些方法各有特点和优势，但同时又有不同的局限性。维纳滤波法是维纳滤波器在信号与干扰不相关条件下的一个特例，适用于白噪声的去除，不适合于场的分离，目前相关研究与应用不多；匹配滤波法对于埋深和尺寸相差较大的场源分离效果较好，但是若异常的功率谱不具明显的线性特征时，如何更加有效地确定深、浅源场的频段分布，是该方法在实际资料处理过程中的一个难点；小波多尺度分析法是一种新兴的位场分离方法，在区域性深部地壳结构研究方面取得了一些成果，但是存在小波母函数选择和阶次选择较难等问题；优选滤波法实现了指定频段或尺度的局部异常的分离，但实际应用中延拓高度往往很难正确选择；优化滤波法异常分离不受延拓高度的困扰，但其前提是目标层与其他深度层场源信息互不相关。实践表明，空间域方法有一定的技术优势。管志宁等^[6]基于多年的重磁数据处理经验认为，由于受测区范围有限、数据离散化，以及不同埋深、不同位置场源体的重磁异常频谱特征部分重叠等诸多因素的影响，在波数域变换过程中会产生明显的误差和畸变，只能用作初步的定性估计；Spector^[21]也认为，波数域方法只能对深源(区域场)和浅源(局部场)异常作定性的划分，不能用于定量研究，无法取代空间域中异常曲线拟合类方法。因而，本文以空间域方法为研究对象。

空间域方法主要包括趋势面拟合法^[22-24]、经验模态分解法^[25]、最小曲率法^[26-28]、曲率滤波法^[29-30]、有限单元法^[31-37]、卡尔曼滤波法^[38-39]、矩阵底秩分解法^[40-41]和插值切割法^[42-53]等。目前，空间域方法中较为前沿且应用较广的重力场分离方法是插值切割法，该方法基于多次切割法发展而来。程方道等^[5]首次提出多次切割法，其基本原理是以四点圆周平均值为切割算子，通过连续切割得到切割区域场，并将其从总异常中减去，便得到切割局部场。为了压制多次切割法中的切割发散效应，文百红等^[42]基于新型算子首次提出了插值切割法；之后，秦葆瑚^[43]、段本春等^[44]、徐世浙等^[45]、邢怡等^[46]、刘东甲^[47]、葛粲等^[48-49]围绕优化切割算子、探寻切割半径与异常埋深的关系展开了系统的研究工作，形成了以插值切割法为核心技术的方法系列，汪炳柱等^[50]、杨金玉等^[51]、万胜荣等^[52]和赵文举等^[53]将该方法广泛应用于油气勘探、矿产勘查、基础地学研究等诸多领域。

以上提及的插值切割法及其改进算法很少涉及重力异常产生的物理属性，多针对异常曲线本身进行插值切割计算。如果忽视异常产生的地质属性，往往会造成分析问题不够全面，甚至产生一定的计算误差。故本文在分析插值切割法基本原理基础上，证明了插值切割法的切割算子仅考虑了局部异常与区域异常弯曲方向相同的叠合情况，不能完全体现实际复杂的地下地质情况。因此，本文结合重力局部异常与区域异常的四种叠合关系，提出了一种基于异常约束条件的插值切割法，并将其命名为约束插值切割法。通过典型理论模型对比试算和实际测量数据应用，验证了约束插值切割法的准确性和优越性，可为大区域的基础地学研究提供技术支撑，在具备局部高精度重力资料的情况下，也可为油气远景区的预测提供一定技术参考。

1 插值切割法

插值切割法是一种简单、有效的空间域非线性滤波方法，以当前计算点场值与周围多点平均值的插值运算为切割算子，通过连续切割得到切割区域场，并将其从位场异常中减去，便得到切割局部场^[42]。根据插值切割算子的不同，形成了基于四点圆周插值切割算子^[42]、八点窗口插值切割算子^[46]和动态改进型插值切割算子^[49]的插值切割方法。

1.1 计算思路

插值切割法虽衍生出不同的切割算子，但其计算原理是一致的。下面以经典的四点圆周插值切割算子为例介绍插值切割法的计算思路。

设重力异常Z(x, y)由区域场R(x, y)和局部场L(x, y)组成，即

$ Z(x, y)=R(x, y)+L(x, y) $

(1)

令M(x, y)表示与点(x, y)距离为r的4个点的重力异常平均值，即

$ \begin{gathered} M(x, y)=\frac{1}{4}[Z(x+r, y)+Z(x-r, y)+ \\ Z(x, y+r)+Z(x, y-r)] \end{gathered} $

(2)

式中r也称为切割半径。区域场R(x, y)是Z(x, y)与M(x, y)的加权平均，即

$ R(x, y)=a M(x, y)+b Z(x, y) $

(3)

式中a和b是加权系数，且a+b=1。令

$ \overline{\Delta H_x}=Z(x, y)-\frac{1}{2}[Z(x+r, y)+Z(x-r, y)] $

(4)

$ \overline{\Delta H_y}=Z(x, y)-\frac{1}{2}[Z(x, y+r)+Z(x, y-r)] $

(5)

$ \Delta H_x=Z(x+r, y)+Z(x-r, y) $

(6)

$ \Delta H_y=Z(x, y+r)+Z(x, y-r) $

(7)

$ \begin{array}{*{20}{c}} {c=\frac{\left(\Delta H_x\right)^2}{\left(\overline{\Delta H_x}\right)^2+\left(\Delta H_x\right)^2}}\\ {\left(0 \leqslant c \leqslant 1\right., 且当 \Delta H_x=\overline{\Delta H_x}=0 时, \left.c=1\right)} \end{array} $

(8)

$ \begin{array}{*{20}{c}} {d=\frac{\left(\Delta H_y\right)^2}{\left(\overline{\Delta H_y}\right)^2+\left(\Delta H_y\right)^2}}\\ {\left(0 \leqslant d \leqslant 1\right., 且当 \Delta H_y=\overline{\Delta H_y}=0 时, \left.d=1\right)} \end{array} $

(9)

$ e=c+d \quad(0 \leqslant e \leqslant 2) $

(10)

取$a=\frac{1}{2} e$，b=1－a，并代入式(3)得到

$ R(x, y)=\frac{1}{2} e Z(x, y)+\left(1-\frac{1}{2} e\right) M(x, y) $

(11)

用上述方法得到的区域场称为第1次切割的区域场，用R₁(x, y)表示；对R₁(x, y)重复使用以上方法，得到第2次切割的区域场R₂(x, y)；依次迭代，当

$ \mathop {\lim }\limits_{n \to \infty } \left| {{R_n}(x,y) - {R_{n - 1}}(x,y)} \right| \to 0 $

(12)

有

$ R(x, y) = \mathop {\lim }\limits_{n \to \infty } {R_n}(x, y) $

(13)

上式得到的区域场即为切割半径的区域场，将其与重力异常Z(x, y)相减，得到局部场

$ L(x, y)=Z(x, y)-R(x, y) $

(14)

一般地，切割半径r越大，得到的局部场所反映的密度体的深度和规模就越大。根据前人大量模拟实验和实测资料^[42-53]分析，以切割半径r得到的局部场可近似代表深度r以上地层的重力效应，但尚未得出切割半径r与密度体埋深的确切计算公式。由于切割算子中的加权系数引入了半二阶差分量，它们与负二阶水平导数即曲率成正比，因此插值切割对不同非线性异常有不同的作用，可提高对不同特征异常场的分离效果。

1.2 误差传递分析

插值切割法采用当前计算点值与四点圆周平均值之间的插值关系作为切割算子G，将切割算子G作用于重力场Z(x, y)，即可得到区域场

$ \begin{aligned} R(x, y)=& G[Z(x, y)]=\frac{1}{2} e Z(x, y)+\\ &\left(1-\frac{1}{2} e\right) M(x, y) \end{aligned} $

(15)

设原始观测数据各点的均方误差δZ(x, y)=ε₀，根据误差传递公式^[54]得到区域场的误差传递关系为

$ \begin{aligned} \varepsilon^2=&\left(\frac{1}{2} e\right)^2[\delta Z(x, y)]^2+\left(1-\frac{1}{2} e\right)^2 \times \\ & {[\delta M(x, y)]^2 } \\ =&\left(\frac{5}{16} e^2-4 e+\frac{1}{4}\right) \varepsilon_0^2 \end{aligned} $

(16)

令误差传递系数$k=\frac{5}{16} e^2-4 e+\frac{1}{4}$，因0≤e≤2，根据一元二次方程最值公式，当e=2时，得到最大值k_max=1；当e=2/5时，得到最小值k_min=1/5。

那么，第1次切割后区域场的误差传递系数k₁有

$ \frac{\sqrt{5}}{5} \leqslant k_1 \leqslant 1 $

(17)

第i次切割后区域场的误差传递系数k_i有

$ \left(\frac{\sqrt{5}}{5}\right)^i \leqslant k_i \leqslant 1 $

(18)

由上式可见，多次切割后区域场的传递误差不会放大。

2 改进的插值切割法

分析插值切割法基本理论公式，该方法主要对异常曲线本身进行插值切割计算，而对重力异常产生的地质属性考虑不足。忽略产生异常的地质属性得到的计算结果难免产生一定的误差。下面首先在理论层面分析插值切割法的不足之处，再从实际地质属性出发，结合重力场局部异常与区域异常的四种叠合关系，对插值切割法加以改进，提出一种基于异常约束条件的插值切割法，即约束插值切割法。

2.1 插值切割法不足的理论分析

根据插值切割法的基本理论可知，切割区域场计算公式(式(15))中，0≤e≤2，不难证明R(x, y)介于Z(x, y)与M(x, y)之间，这就意味着经插值切割法计算得到的区域场值不能超出总场异常值与周围点平均值的取值区间。这种情形仅代表局部异常与区域异常弯曲方向相同的叠合情况，并不能完全代表复杂多变的实际地下地质情况，下面进一步分析局部异常与区域异常不同叠合关系的情形。

2.2 局部异常与区域异常叠合关系分类

在勘探地球物理领域，一般认为区域异常是由埋藏较深、水平延伸范围较大的地质体或构造所产生的，在异常图中表现出形态宽缓、幅值较大、范围较广的特征；局部异常则是由埋藏较浅、水平延伸范围较小的地质体所产生的，在异常图中表现出形态陡窄、幅值较小、异常范围小的特征。在遵循这一基本原则的前提下，本文针对四种实际地质情形，总结了相应的局部异常与区域异常的四种叠合关系分类情形，详见表 1和图 1。不难发现，式(11)仅适用于分类情形1和2，即区域异常介于总场异常值与周围点平均值的取值区间。然而，对于区域异常与局部异常曲线弯曲方向相反的情形，即分类情形3和4，如仍使用式(11)计算切割区域场，计算结果虽可无限向周围点平均值逼近，但依然忽略了图 1c和图 1d中黑色竖线所示部分的异常，这部分异常代表计算得到切割区域场与理论区域场的差异。

2.3 算法改进方案

为弥补传统插值切割法的不足，本文提出了一种基于异常约束条件的改进方案，其基本思路是：首先，给定切割半径r，分析区域异常和局部异常曲线弯曲方向；其次，对区域异常与局部异常曲线弯曲方向相同的情形(即分类情形1和2)，采用式(11)计算切割区域场，对弯曲方向相反的情形(即分类情形3和4)，在式(11)基础上，增加切割区域异常修正量ΔZ(x, y)(图 2)，即

$ \begin{aligned} R(x, y)=& \frac{1}{2} e Z(x, y)+\\ &\left(1-\frac{1}{2} e\right) M(x, y)+\Delta Z(x, y) \end{aligned} $

(19)

首先考虑x方向异常剖面，假设局部异常附近的区域异常是曲率半径为N、曲率中心为O的曲率圆的部分弧段(图 2中红色曲线与蓝色曲线重合部分)，则x方向切割区域的异常修正量为

$ \Delta Z(x)=N-\overline{O C}=N-\sqrt{N^2-\overline{A C^2}} $

(20)

式中OC、AC分别为线段OC和AC的长度，其中AC可根据勾股定理求得

$ \overline{A C}=\frac{1}{2} \sqrt{\overline{A B}_x{ }^2+{\overline{A B_z}}^2} $

(21)

其中AB_x、AB_z分别为线段AB在x轴和z轴的投影长度。曲率半径N可根据曲率圆的曲率K求得

$ N=\frac{1}{K}=\frac{\left(1+Z_{D, x}^{\prime 2}\right)^{\frac{3}{2}}}{Z_{D, x}^{\prime \prime}} $

(22)

式中Z′_D，x、Z″_D，x分别为Z(x, y)在D点对x的一阶导数和二阶导数。用一阶差分替代一阶导数，并忽略一阶余量，得到

$ Z_{D, x}^{\prime} \approx \frac{Z(B)-Z(A)}{\overline{A B}} $

(23)

式中：Z(A)、Z(B)分别为Z(x, y)在A点和B点的值；AB表示线段AB的长度。用二阶差分替代二阶导数，并忽略二阶余量，则

$ Z_{D, x}^{\prime \prime} \approx \frac{Z(A)+Z(E)-2 Z(B)}{\overline{A B}_x} $

(24)

同理，可得到y方向切割区域的异常修正量ΔZ(y)，并取ΔZ(x, y)=[ΔZ(x)+ΔZ(y)]/2。

约束插值切割法详细计算步骤见附录A。以上算法基于MATLAB(R2008a)编程平台实现程序和界面编制。

3 理论模型试算 3.1 理论模型设计

本文理论模型设计参考了程方道等^[5]、Kaftan等^[31]、Zhu等^[40]及刘东甲^[47]的设计思路，即深部大尺度密度体产生的异常视为区域异常，浅部小尺度密度体产生的异常视为局部异常，不同之处是在同一个理论模型中综合考虑了局部异常与区域异常的四种叠合关系。鉴于2.2部分中分类情形1、3分别与2、4呈镜像对称关系，可合并处理，故本文设计了图 3所示的理论模型。三个不同深度的浅部球体(编号①、②、③)产生的是局部重力异常，其中浅部球体①、②产生局部重力高，浅部球体③产生局部重力低；三个不同深度的深部球体(编号④、⑤、⑥)产生的是区域重力异常，其中球体④、⑥引起的是区域重力高异常，球体⑤引起的是重力低异常。

3.2 理论模型正演

理论模型的重力异常为各球体剩余质量的引力位沿重力方向(定义为z方向)的导数之和^[55]，即

$ \Delta g(x, y, z) = \sum\limits_{i = 1}^n {\frac{{\partial {V_i}}}{{\partial z}}} = \sum\limits_{i = 1}^n {\frac{{f{m_i}\left( {{z_{i0}} - z} \right)}}{{{r_i}}}} $

(25)

式中：V_i为球体i产生的引力位；n表示异常体数量，此例n=6；f=6.672×10^-11m³·kg^-1·s^-2为万有引力常数；m_i=43πR_i³σ_i为球体i的剩余质量，这里R_i和σ_i分别为球体i的半径和剩余密度；r_i=$\sqrt{\left(x_{i 0}-x\right)^2+\left(y_{i 0}-y\right)^2+\left(z_{i 0}-z\right)^2}$为测点(x，y，z)与球体i的球心的距离，这里(x_i0, y_i0, z_i0)表示球体i的球心坐标。

基于式(25)，通过正演计算可得模型3的正演重力异常(图 4、图 5)。由图 4可见，局部异常基本淹没在区域异常之中。由图 5可见，1、2、3号异常分别对应2.2部分中分类情形1、4、3。

3.3 理论模型重力场分离

分别采用插值切割法、约束插值切割法、小波多尺度分解法^[13]对理论模型正演结果进行重力场分离，分离得到的区域场和局部异常见图 6。基于切割半径与密度体埋深近似相等的对应关系，计算过程中以r=400m为初始值，以100m为增减量，采用多个切割半径r，经对比筛选，确定取r=500m。根据小波多尺度分解的尺度场源深度转换公式^[14]，因网格间距为100m，故采用1~4阶小波细节代表浅部球体产生的局部重力异常，用4阶小波逼近代表深部球体产生的区域重力异常。对比约束插值切割法和小波多尺度分解法得到的局部异常，可见前者更接近于理论局部异常(图 4c)。由于区域异常幅值较大，基本淹没了局部异常，故通过三种方法分离出的区域异常(图 6上)与理论区域异常(图 4b)基本一致。为了更直观地分析本文方法的优越性，绘制了沿y=0剖面的三种方法分离得到的局部异常曲线(图 7)。由图可见，由于1号异常为区域重力高叠加局部重力高(分类情形1)，三种方法的结果基本一致，都能较准确地体现分类情形1所产生的局部重力异常；关于2号和3号异常，本文方法得到的分离结果明显改善了插值切割法的计算结果，且与小波多尺度分解的结果吻合。三种方法均在局部异常两翼附近产生了小幅度假异常，这是因为方法本身存在一定的计算误差。以上分析证明了本文提出的约束插值切割法基本理论的正确性和分离效果的正确性。

为了定量对比分析插值切割法、约束插值切割法和小波多尺度分解法的重力场分离效果，用下式计算局部异常Δg_Ⅰ与理论局部异常Δg_Ⅱ的相关系数

$ {\mathop{\rm Re}\nolimits} \left( {\Delta {g_{\rm{I}}}, \Delta {g_{{\rm{II}}}}} \right) = \frac{{{\mathop{\rm Cor}\nolimits} \left( {\Delta {g_{\rm{I}}}, \Delta {g_{{\rm{II}}}}} \right)}}{{\sqrt {{\mathop{\rm Var}\nolimits} \left( {\Delta {g_{\rm{I}}}} \right){\mathop{\rm Var}\nolimits} \left( {\Delta {g_{{\rm{II}}}}} \right)} }} $

(26)

式中：Cor(·)表示协方差；Var(·)表示均方差。经计算，插值切割法、约束插值切割法和小波多尺度分解法得到的局部异常与理论局部异常的相关系数分别为0.42、0.75、0.70，可见本文提出的约束插值切割法得到的局部异常与理论局部异常具有最强的相关性。

另外，为了评价约束插值切割法对插值切割法的改善效果，引入评价航磁动态补偿效果的“改善率”(Improvement Ratio，IR)^[56]评价改善效果，其计算公式为

$ \mathrm{IR}=\frac{\operatorname{Std}\left(\Delta g_{\mathrm{cz}}-\Delta g_{\mathrm{II}}\right)}{\operatorname{Std}\left(\Delta g_{\mathrm{ys}}-\Delta g_{\mathrm{II}}\right)} $

(27)

式中：Δg_cz表示插值切割法得到的局部异常；Δg_ys表示约束插值切割法得到的局部异常；Std(·)表示求标准差。

因为本文提出的新方法仅对2号、3号异常进行了优化，故仅对2号、3号局部异常的改善率进行评估。经计算，约束插值切割法对插值切割法的改善率IR为1.48，可见相比于插值切割法，新方法计算的局部异常效果改善了48%，优化效果明显。

4 应用实例 4.1 中国大陆及周边区域重力数据分离

本文以中国大陆及周边区域的重力数据为试验对象，范围为东经63°~136°，北纬15°~58°。重力数据来自美国国家地理空间情报局地球重力场研发小组发布的超高阶地球重力场模型EGM2008(Earth Gravitational Model 2008)，该模型的阶次完全至2159次，空间分辨率约为5′，可用于小比例尺重力编图和构造研究^[57]。本实例使用的布格重力异常的网格间距为4km。

分别利用插值切割法和本文提出的约束插值切割法对布格重力异常进行分离。为了展示新方法的正确性和改进效果，以葛粲等^[48]的研究结果为参照。根据徐世浙等^[8]、葛粲等^[48]提出的“切割半径应近似等于异常体埋深”的层切割思想，本例的切割半径r以网格间距4km为基准，并以其整数倍的形式逐步递增，每次切割得到的区域场作为下次切割的输入，并结合前人经验总结的“切割半径r得到的局部场可以近似代表深度r以上的地层重力效应”，当切割半径为19倍网格间距(r=76km)时，计算得到的局部异常可视为72~76km参考深度层的重力异常值，并将其视为上地幔顶部附近介质产生的重力异常。插值切割法和约束插值切割法计算得到的上地幔顶部附近重力异常分别见图 8和图 9。

对比图 8与葛粲等^[48]的计算结果，可见重力异常分布特征基本一致，验证了本文编程实现的插值切割法的正确性。对比插值切割法(图 8)与约束插值切割法(图 9)得到的上地幔顶部附近重力异常，可见整体上重力异常形态相似，但约束插值切割法得到的正值重力异常幅度和范围相对更大，说明在负值背景下叠加正异常的情况下，约束插值切割法对异常分布进行了合理调整，这与约束插值切割法的基本原理和理论模型结论相一致。

为了定量对比约束插值切割法和插值切割法得到的重力异常，计算了两者的差值(图 10)，即前文所述的修正量。由图可见，在准噶尔盆地、喜马拉雅地块等局部地区，本文方法的修正达50%以上。

根据图 9可得到初步的地质认识，即在喜马拉雅地块边界局部地区(图中绿色虚线椭圆区域)，深部高密度支撑物质强烈北移侵入，而柴达木盆地和准噶尔盆地正高值范围的扩大(图中白色虚线所示)，则意味着稳定盆地下方可能存在高密度物质。

杨文采等^[58]利用小波多尺度分解法对中国陆域及周边布格重力异常进行了分离，获取了包括上地幔顶部异常(D7+D8)在内的4个深度界面所产生的重力异常。对比约束插值切割法与小波多尺度分解法对上地幔顶部异常的分离结果，可见在青藏高原等区域的异常分布特征与杨文采等^[58]的研究成果相近，但约束插值切割法得到的异常细节更加丰富，主要是因为小波多尺度分解法仅将布格重力异常分离为4层，而约束插值切割法可分离为20层，对异常的分辨率更高^[48]。

由于篇幅限制，结合地震、地质等资料的综合地质解释将另文详细阐述。据本例分析可见，约束插值切割法在实际应用中实现了对插值切割法计算结果的修正，可得到更准确的重力异常、更明确的重力异常梯度带。

4.2 中国南海万安盆地及周边区域重力数据分离

为了进一步验证约束插值切割法的实际应用效果，选取中国南海万安盆地及周边区域布格重力异常数据为研究对象，范围为东经104°~111°，北纬5°~11°，数据来源同4.1部分重力数据。陈玲等^[59]认为南海万安盆地新生界沉积层厚度最大为12.5km，冯旭亮等^[60]的研究结果显示该区域最大厚度大于10km，故本例约束插值切割法中选取切割半径为3倍网格间距(12km)，经计算得到万安盆地及周边区域新生界沉积层的重力异常(图 11a)。图 11b为选取的一条地震剖面所对应的新生界沉积层重力异常，图 11c为该地震剖面构造解释剖面^[60]。

为了说明本文方法计算结果的可靠性，与冯旭亮等^[60]采用最小曲率位场分离方法得到的计算结果(文献[60]中的图 7、图 8)进行了对比，可以发现，图 11a中重力异常的整体分布特征与文献[60]中的图 7较吻合，图 11b中的重力异常曲线走势与文献[60]中的图 8a大致相同，但本文的计算结果更能反映沉积基底的起伏变化细节，特别是更加直观地显示了F1断层右侧沉积基底的局部凸起(图 11c中的虚线圈所示)。因此，基于本文方法分离得到的局部重力异常，可更清晰地解释新生界沉积层底界面的深度变化，为该区域新生界深度反演研究提供了较可靠的数据基础，同时也验证了约束插值切割法的实用性。

5 结论

本文通过基本理论推导，并结合实际重力异常的叠合分类情形，实现了对插值切割法的优化和改进，提出了基于异常约束条件的插值切割法(约束插值切割法)，并开展了理论模型试验和卫星重力数据处理、分析。得出以下几点认识。

(1) 通过对插值切割法切割算子的公式推导，结合四种情形的异常叠合关系，证明了插值切割法的切割算子仅考虑了局部异常与区域异常弯曲方向相同的叠合情况。基于此，给出了插值切割算子的补偿方法，提出了基于异常约束的插值切割法，制定了新方法的计算流程，并基于Matlab平台编程实现。

(2) 理论模型试算结果显示，插值切割法、约束插值切割法、小波多尺度分解法中约束插值切割法得到的局部异常与理论局部异常具有最强的相关性(相关系数为0.75)，与插值切割法相比，对局部异常效果的改善达48%，优化效果明显。

(3) 对中国陆域及周边区域的布格重力异常进行了分离，获取了72~76km参考深度层产生的重力异常，约束插值切割法获取了该层位更真实可靠的重力异常值信息，在准噶尔盆地、喜马拉雅地块等局部地区的修正达50%以上，可为上地幔顶部构造特征研究提供参考。

(4) 利用约束插值切割法分离得到了万安盆地及周边区域新生界沉积层重力异常，与地震剖面显示的基底起伏特征基本一致，说明该方法可用于基底深度变化的反演研究。

总之，通过理论模型数据和实际重力数据分析，约束插值切割法计算结果均明显优化了插值切割法的计算结果，验证了前者理论的正确性及实用性，可为实际重力数据的地质解释提供更准确、可靠的重力异常信息。本文提出的新方法也可用于磁场数据的异常分离。

附录A 约束插值切割法计算步骤

(1) 根据目标层位深度信息，对重力异常Z(x, y)设置切割半径r；

(2) 利用式(2)计算与点(x, y)相距r的某4个点重力异常的平均值M(x, y)；

(3) 利用式(10)计算加权参数e;

(4) 比较区域异常与局部异常曲线弯曲方向，具体可细分为2步：

① 计算下列参数

$ \begin{aligned} &I_{x 1}=Z(x+r, y+r)-\frac{Z(x+2 r, y+r)+Z(x, y+r)}{2} \\ &I_{y 1}=Z(x+r, y+r)-\frac{Z(x+r, y+2 r)+Z(x, y+r)}{2} \\ &I_{x 2}=Z(x+2 r, y+r)-\frac{Z(x+4 r, y+r)+Z(x, y+r)}{2} \\ &I_{y 2}=Z(x+r, y+2 r)-\frac{Z(x+r, y+4 r)+Z(x+r, y)}{2} \end{aligned} $

② 当比较参数不能同时满足下列条件时，表明区域异常和局部异常弯曲方向相同(即分类情形1和2)

$ \begin{aligned} &I_{x 1} I_{y 1}>0 \\ &I_{y 1} I_{y 2}<0 \\ &I_{x 1} I_{x 2}<0 \\ &I_{y 2} I_{x 2}>0: \end{aligned} $

当比较参数同时满足下列条件时，表明区域异常和局部异常弯曲方向相反(即分类情形3和4)

$ \begin{aligned} &I_{x 1} I_{y 1}>0 \\ &I_{y 1} I_{y 2}<0 \\ &I_{x 1} I_{x 2}<0 \\ &I_{y 2} I_{x 2}>0: \end{aligned} $

(5) 对于弯曲方向相同的情形，区域场等于Z(x, y)与M(x, y)的加权平均，采用式(11)计算切割区域场；

(6) 对于弯曲方向相反的情形，区域场为Z(x, y)与M(x, y)的加权平均的基础上增加切割区域异常修正量ΔZ(x, y)，此种情况可细分为以下4步：

① 首先考虑x方向异常剖面，假设局部异常附近的区域异常为曲率半径为N、曲率中心为O的曲率圆的部分弧段，利用式(20)计算x方向切割区域异常修正量ΔZ(x)。

② 同理可得y方向切割区域异常修正量ΔZ(y)；

③ 取ΔZ(x, y)=[ΔZ(x)+ΔZ(y)]/2；

④ 计算修正量ΔZ(x, y)后，采用式(19)计算切割区域场；

(7) 用上述方法得到第1次切割的区域场R₁(x, y)，对R₁(x, y)重复使用步骤(1)~(6)，得到第2次切割的区域场R₂(x, y)，依次迭代，直到满足式(12)，迭代终止，得到切割区域场；

(8) 利用式(14)计算最终的局部场R(x，y)。

在构思本文的初期，与插值切割法的原创者——中国石油勘探开发研究院文百红老师进行了有益探讨，特别是在插值切割法的细节上得到了文百红老师的悉心指导，在此深表谢意；感谢中国自然资源航空物探遥感中心骆遥教授为文章的修改提出了建设性的意见！