0 引言

运用位场数据识别地质体边界是对地球物理资料进行地质解释的一项重要内容。近年来，在利用位场资料推算断裂构造的水平位置及深度方面已出现了许多方法。其中基于位场梯度的边界识别方法算法^[1]简单、易于实现，是一类常用的方法(如水平总梯度^[2-3]、垂向导数^[4-5]、Theta map^[6]、Tilt angle^[7]、归一化标准差^[8]等)。但是这类方法对区域边界信息不敏感，且易受随机干扰的影响。Euler反褶积法^[9-11]也可以用于边界识别，该方法是通过计算位场一阶导数和点位坐标实现的，导数计算结果或计算点位稍有偏差，便会在反演中产生较大的偏离，从而极大地影响边界反演精度。利用重力归一化总梯度及其相位同样可以确定场源深度和边界水平位置^[12], 但此法受谐波数、向下延拓的影响，场源深度和水平位置的估计受到了一定限制，因此精度不高。

小子域滤波法是杨高印^[13]提出的一项位场数据处理技术，是对滑动平均法的一种改进，该方法能够增强区域异常的界限特征。夏玲燕等^[14]将小子域滤波和Tilt梯度相结合，许海红等^[15]将小子域滤波和水平总梯度相结合，均获得了较良好的边界识别效果。马涛等^[16]针对小子域滤波存在偏心问题对梯度保持滤波法进行改进，使滤波结果更加稳定可靠；肖锋等^[17]采用“田”字型子域划分模式，使计算更便捷；马国庆等^[18]针对小子域滤波存在异常曲线扭曲的缺陷，提出了优化小子域滤波，更好地展示了区域异常间的界限。段晓旭^[19]采用加密点距的方式解决小子域存在异常曲线扭曲的问题，并提出点位移动法增强异常梯度带。蔡钟等^[20]针对小子域滤波存在边界异常失真的现象，提出了五/六边形子域滤波法，在一定程度上改善了小子域滤波的效果。许海红等^[21]对上述多种小子域滤波方法进行了对比分析，给出了不同方法的优缺点及适用条件。但是上述改进方法仅能对具有明显梯级带特征的断裂构造进行识别，而无法检测具有弱梯级带或非梯级带特征的场源边界。为此，张凤旭等^[22-23]提出了三方向小子域滤波进行断裂构造识别，提高了边界识别能力；Jiang等^[24]对小子域滤波的剖分模式进行了改进，并实现了基于梯度张量的小子域滤波，同样获得了丰富的边界信息；Tai等^[25]在三方向小子域滤波基础上，提出了加权小子域滤波法，使断裂构造的划分更快捷、更精确。但这三种改进方法均采用了差分或导数计算，因此计算更易受到干扰的影响。

针对上述方法的优缺点，本文提出了新的子域剖分模式，使其可以更准确地检测出不同方向的异常界限; 同时为了消除干扰对小子域滤波的影响和提高异常分辨率，给出了一种新的滤波技术——基于迭代差分的稳定增强滤波。

1 改进小子域滤波 1.1 传统小子域滤波

传统小子域滤波^[13]是将滑动窗口按中心点的不同侧面划分出8个子域，并将8个子域均方差最小子域的平均值作为滤波输出，图 1为传统小子域滤波子域剖分模式的示意图。该方法可以较好地保留异常间的界限，但还存在着以下不足：①虽然对位场异常进行了低通滤波处理，然而不同子域对噪声的敏感程度不同，从而导致滤波结果中出现明显的“虚假”异常^[26]；②8个子域的重心均偏离窗口中心点，易造成异常形态不规则扭曲，从而产生虚假信息^{[15, 18]}；③不能较好地检测出梯级带不明显或非梯级带特征的异常界限，且检测结果存在较大偏差^[22]。

1.2 改进型小子域滤波

针对传统小子域滤波存在的问题，本文对子域剖分模式进行改进，改进型小子域滤波与传统方法的区别在于剖分子域的个数不同。以5×5窗口为例，传统剖分是按中心点的不同侧面划分出8个子域(图 1)，而本文方法的剖分是将5×5窗口剖分出21个子域，每个子域包含9个数据点(图 2)。从图 2中可以看出，改进后的子域划分方式可以检测任意方向上的异常变化，能更好地体现出不同走向上的构造，获得更为准确、丰富的异常界限。该方法的计算流程与传统方法一致，即将窗口内21个子域均方差最小的子域数据平均值作为窗口中心点的滤波结果，以紧缩梯级带作为场源边界的识别依据。

对于N×N窗口来说，子域剖分个数为$\frac{(N+1)(N+9)}{4}$，每个子域的数据点为$\frac{(N+1)^{2}}{4}$。也就是说，随着窗口增大，子域剖分个数也随之增加，而传统小子域滤波及其他改进方法的子域剖分个数大多是固定的，与窗口大小无关。因此本文提出的子域剖分模式则更能有效地检测异常的细节变化。

2 异常增强滤波技术

传统滤噪技术一般是低通滤波。低通滤波在消除噪声时，对低频有效信号同样也有一定的消减，从而导致异常分辨率降低。因此，在原异常中表现为非梯度带的异常界限，小子域滤波仍不能有效识别。鉴于此，本文提出了一种稳定增强滤波技术，该方法在滤除噪声的同时，可以提高对异常的分辨能力。

2.1 基本原理

设位场异常u的波谱为U，低通滤波算子为φ，0≤φ≤1，滤波后的波谱可以表示为

$ U_{\varphi}=U \varphi $

(1)

虽然上述低通滤波可以滤除随机干扰，但对低频有效信息也存在一定的消减作用。因此，采用迭代补偿模式完成对低频成分的弥补，迭代过程如下：

(1) 将U与U_φ的差值利用滤波算子φ进行修正，修正结果为

$ U_{\varphi}^{(1)}=\left(U-U_{\varphi}\right) \varphi+U_{\varphi}=U \varphi[1+(1-\varphi)] $

(2)

(2) 将U与U_φ⁽¹⁾的差值再次利用φ进行修正，修正结果为

$ \begin{aligned} U_{\varphi}^{(2)} &=\left(U-U_{\varphi}^{(1)}\right) \varphi+U_{\varphi}^{(1)} \\ &=U \varphi\left[1+(1-\varphi)+(1-\varphi)^{2}\right] \end{aligned} $

(3)

(3) 按照如上迭代补偿模式，第n次修正结果为

$ U_{\varphi}^{(n)}=\left(U-U_{\varphi}^{(n-1)}\right) \varphi+U_{\varphi}^{(n-1)} $

(4)

根据式(2)~式(4)，得到迭代通式

$ \begin{aligned} U_{\varphi}^{(n)}=& U_{\varphi}\left[1+(1-\varphi)+(1-\varphi)^{2}+\cdots+\right.\\ &\left.(1-\varphi)^{n}\right]=U\left[1-(1-\varphi)^{n+1}\right] \end{aligned} $

(5)

上述迭代模式事实上与侯重初^[27]提出的补偿圆滑滤波在原理上是一致的。由式(5)可以看出，采用迭代法得到的波谱U_φ⁽ⁿ⁾等价于在波谱U中加入了一个与迭代次数n相关的低通滤波1-(1-φ)ⁿ⁺¹，因此迭代法仍是起到压制干扰的作用，而无法提高异常的分辨率。

为了提高异常分辨率，对迭代结果进行如下处理：(1)设初始迭代次数为n的滤波结果u_φ⁽ⁿ⁾=F^-1[U_φ⁽ⁿ⁾]，迭代次数间隔为Δn的滤波结果u_φ^(n+Δn)=F^-1[U_φ^(n+Δn)]，则一阶迭代差分的结果为

$ \begin{aligned} \mathrm{d} u_{\varphi}^{1 , n} &=u_{\varphi}^{(n+\Delta n)}-u_{\varphi}^{(n)}=\mathrm{F}^{-1}\left[U_{\varphi}^{(n+\Delta n)}-U_{\varphi}^{(n)}\right] \\ &=\mathrm{F}^{-1}\left\{U(1-\varphi)^{n+1}\left[1-(1-\varphi)^{\Delta n}\right]\right\} \end{aligned} $

(6)

式中F^-1表示反傅里叶变换。

(2) 令迭代次数间隔为2Δn的滤波结果为u_φ^(n+2Δn)=F^-1[U_φ^(n+2Δn)]，则二阶迭代差分的结果为

$ \begin{aligned} \mathrm{d} u_{\varphi}^{2, n} &=2 u_{\varphi}^{(n+\Delta n)}-\left[u_{\varphi}^{(n)}+u_{\varphi}^{(n+2 \Delta n)}\right] \\ &=\mathrm{F}^{-1}\left\{U(1-\varphi)^{n+1}\left[1-(1-\varphi)^{\Delta n}\right]^{2}\right\} \end{aligned} $

(7)

(3) 仿照上述流程，第m阶迭代差分结果为

$ \mathrm{d} u_{\varphi}^{m , n}=\mathrm{F}^{-1}\left\{U(1-\varphi)^{n+1}\left[1-(1-\varphi)^{\Delta n}\right]^{m}\right\} $

(8)

由式(8)可知，迭代差分m次的滤波因子Φ=(1-φ)ⁿ⁺¹[1-(1-φ)^Δn]^m事实上是一个低通滤波Φ_L=[1-(1-φ)^Δn]^m与高通滤波Φ_H=(1-φ)ⁿ⁺¹的组合。

2.2 滤波特性

令滤波算子φ=exp(-wh)，m=1，Δn=10，其中w表示波数，h=10Δx，Δx=0.1km。图 3a给出了上文提及的滤波算子的滤波曲线。可以看出，低通滤波φ(向上延拓算子)虽然对高频成分具有较好的压制作用，但对中、低频的有用信号一定程度上也有消减作用；相对于低通滤波φ，增强滤波器中的低通滤波Φ_L可以更好地保留低频成分；稳定增强滤波Φ在中、低频则趋近于高通滤波算子Φ_H，而在高频处则趋于低通滤波Φ_L，即滤波算子Φ为一带通滤波器，对中频成分起到了放大作用，而对高频成分起到了压制作用。也就是说，稳定增强滤波后的位场异常在理论上不仅具有较强的稳定性，而且还可以提高原异常的分辨能力，因此，将这种迭代差分滤波称为稳定增强滤波。

图 3b~图 3d是不同参数的增强滤波响应曲线。可以看出，随着初始迭代次数n的增加(图 3b)，滤波算子的主峰向高频移动，但幅值却随之降低；随着迭代次数间隔Δn的增加(图 3c)，增强滤波算子频带变宽，主峰同时向高频移动；随着迭代差分阶数m的增加(图 3d)，滤波算子频带变窄且主峰向低频移动。因此选择较少的初始迭代次数n、较小的迭代次数间隔Δn和较大的迭代差分次数m，会使滤波结果更稳定。

理论上，稳定增强滤波器中的低通滤波Φ_L、迭代初始次数n、迭代间隔Δn及迭代差分阶数m的选择具有随意性，但为确保计算结果的稳定性及异常分辨能力，建议φ选用向上延拓不少于5倍点距的延拓算子，选取n=1、Δn < 50及m < 5。后文中的模型及实例部分低通滤波算子均是选取向上延拓10倍的延拓算子，m=n=1，Δn=10。

3 理论模型分析

为了验证本文方法的有效性和优越性，建立了由4个埋深不同、尺度不同、密度不同的重力异常体(表 1)组成的复杂模型进行试验。同时，在模型的正演重力异常中添加了随机干扰，以检验方法的稳定性。图 4是该模型产生的重力异常和分别添加1%、3%噪声的重力异常及利用本文方法对图 4c进行稳定增强滤波处理的结果。可以看出，理论重力异常(图 4a)可以通过梯度带大致识别异常体A、B、C的边界位置，但梯度带较宽，边界位置不明确；模型体D的异常受叠加场影响，圈闭极值与模型体中心不对应，在异常图主要以等值线同向扭曲为主。在重力异常含1%随机噪声时(图 4b)，异常体D所对应的异常更加模糊。当重力异常含3%噪声时(图 4c)，难以通过异常梯度带识别所有模型边界。稳定增强滤波后的重力异常(图 4d)不仅等值线相对于图 4c更加圆滑，且异常分辨率也得到明显提高，尤其对规模小、埋深大的异常体D。

图 5是无噪数据在不同窗口下的传统小子域滤波与本文改进方法滤波处理结果。可以看出，当计算窗口较小时(如3×3)，传统小子域(图 5a上)与改进小子域滤波(图 5a下)仅对埋深较浅的异常体C的边界异常进行了有效紧缩；随着窗口尺度的增加，无论是传统小子域滤波还是改进小子域滤波，地质体边界异常被紧缩得更加明显，如9×9窗口时，除异常体D缺失部分边界，其他三个异常体的边界异常均得到了有效增强。由图 5可以明显发现，传统小子域滤波结果在模型体角点处存在明显的等值线扭曲现象，而改进型小子域滤波的处理结果并不存在此现象，能更好地识别异常体的边界。

图 6是不同噪声水平的重力异常在不同窗口下的传统小子域滤波与改进小子域滤波处理结果对比。可以看出，由于随机噪声对不同方向子域的影响程度不同，导致了两种方法的滤波结果中均存在着多条不连续分布的异常紧缩带，且无法确定哪些紧缩带是有效的；而且，噪声水平越高，小子域滤波处理结果中的虚假信息就越多，处理结果的可靠性越差。

图 7是含3%噪声重力异常经过稳定增强滤波后的传统小子域滤波和改进型小子域滤波处理结果对比。相对于图 6来说，稳定增强滤波后的小子域滤波更加稳定，不但不包含明显的虚假紧缩带，更能有效地检测出地质体边界，尤其异常体D的边界也有明显的展示。

4 实际重力资料应用

为了验证稳定增强异常改进型小子域滤波对实际资料的处理效果，选取吉林省南部鸭绿江盆地实测重力数据(点距和线距均为1km)进行试验。鸭绿江盆地位于中朝板块东北缘，二级大地构造单元隶属于辽东台隆区，盆地主体为太子河—浑江坳陷，盆地内除志留系、泥盆系地层缺失外，其他时代的岩石均有出露(图 8)。研究区存在一系列逆冲推覆构造，使得地层关系十分复杂、连续性极差及地层倒置现象非常普遍^[28]，这给地质体边界的确定甚至地质解释都带来了较大难度。研究区内各个时代的岩石出露较全，奥陶系与石炭系间存在明显的密度差(表 2)，这为方法试验的可靠性提供了较好的物性基础。

从盆地内布格重力异常图(图 9a)可以看出，反映地质体边界位置或断裂构造的梯级带异常梯度较平缓，因而地质体边界不易直接根据重力异常进行精确解释。较小型地质体的边界受区域异常的影响较大，在异常图中主要表现为异常等值线突然变宽或变窄以及同形扭曲等非梯级带特征，地质体的边界位置确定难度较大。从图 9b可以看出，经稳定增强处理的重力异常图的分辨率获得了明显提高，不仅异常梯级带连续性增强，且异常呈现更明显的NE走向，与地质图(图 8)中的岩石分布呈现出了很强的一致性。图 10是布格重力异常的传统小子域和改进型小子域滤波结果，可以看出，对于原始布格重力异常，传统小子域滤波(图 10上)对大型异常梯度带进行了较好的紧缩，然而紧缩带的连续性较差，且存在明显的无规律性弯曲；而改进型小子域滤波(图 10下)同样也反映了大型断裂的位置，且边界异常信息更丰富，紧缩带的连续性也更强。对于稳定增强异常，无论是经传统小子域滤波(图 11上)还是改进型小子域滤波(图 11下)处理的重力数据对边界的识别能力显然都获得了大幅度的提升，边界异常信息更丰富，而改进型小子域滤波结果显示出更丰富的构造边界信息，且紧缩带的走向和连续性也更好。需要指出的是，5×5窗口的小子域滤波(图 11a)数据中的一些紧缩带在9×9窗口的小子域滤波中力图(图 11b)上没有显示。如在大安镇的东南侧，5×5窗口的小子域滤波数据中存在两条NE走向的异常紧缩带，但9×9窗口的小子域滤波数据中并未被有效识别出来，且检测出的紧缩带连续性也较差。这说明在实际资料处理中，需要结合不同窗口的小子域滤波数据进行地质解释。另外，新方法获得的异常紧缩带与不同时代的岩性界限(图 8)具有良好的对应性；增强异常的改进小子域滤波数据的负异常分布也较好地反映了研究区中—新生代地层和燕山期花岗岩等低密度岩石的分布。

图 12a是结合区域地质及前人工作成果绘制的构造分区图。图中给出了大型断裂的位置、低密度岩石的分布区(燕山期花岗岩及石炭—白垩纪地层)以及浑江煤田工作区等信息，以此验证本文方法的可靠性。图 12b是基于图 11下解释的3条大型断裂和17条中小型断裂以及在浑江坳陷内圈定的负异常区域，其中F1为二道江—江源断裂，是龙岗隆起与浑江坳陷的界限断裂；F2为石人—大青沟断裂，是浑江坳陷与老岭隆起的分界断裂；F3为鸭绿江断裂，在研究区东部则是老岭隆起与四道沟坳陷的界限。这三条大型断裂在重力异常中均以大型梯级带为异常特征，其余17条中、小型断裂在重力异常图上表现为弱异常梯级带或等值线突然变宽或变窄以及同形扭曲等异常特征。需要指出的是，在浑江坳陷内，地表出露的低密度岩石分布区仅占30%左右的面积，而本文方法识别出的负异常分布范围可达50%以上，也就是说浑江坳陷内有大面积的低密度体隐伏于高密度岩石之下，其中勘探程度较高的浑江煤田工作区中的一些飞来峰下的煤层(主要分布在石炭、二叠、侏罗纪地层中)已被钻探或采掘所证实^[28-29]，这也再一次证实了本文方法的实用性，同时该方法也为研究区后续煤田勘查有利区(尤其勘探程度较低的六道沟—三道湖—石人镇一带)和地质—地球物理综合解释的有效手段之一。

5 结论

针对传统小子域滤波存在异常曲线扭曲的缺陷，本文提出了一种包含更多子域的新子域剖分模式，可以更准确地检测不同走向构造的界限。针对小子域滤波存在易受干扰影响和难以识别非梯级带特征的断裂构造的缺陷，提出了基于迭代差分的稳定增强滤波技术。该技术在滤除噪声干扰的同时，还可以提高异常的分辨率。模型试验和实例应用表明，改进型小子域滤波比传统算法效果更佳，异常的界限可以得到更加准确、精细的刻画；稳定增强滤波后的重力异常不仅稳定性强，且异常分辨率明显提高；增强异常改进型小子域滤波则能够更加准确、精细地刻画密度异常体的界限，同时还可以通过重力正、负异常的分布进一步提高方法的解释能力。