重磁异常是地下物质密度和磁性不均匀分布的综合反映^[1]，在研究地质目标体的横向不均匀性，特别是地质目标体边缘位置时，重磁位场有其独特的优势^[2]。利用重磁位场研究推断地质体边界位置之前，需要对重磁位场数据进行滤波处理，以便从综合的异常数据中提取能够突出目标体的异常。但传统的滤波方法会模糊异常之间的界限，降低异常分辨率^[1]，为此，杨高印^[1]提出一种非线性滤波方法——小子域滤波法，该方法在进行低通滤波的同时，能以较高的分辨率保留异常中的区域特征；严良俊等^[3]称其为梯级带滤波增强技术；马涛等^[4]对其改进并称之为网格数据保持梯度滤波法；张凤旭等^[5-7]利用小子域滤波法分别检测3个方向偏移处理结果中的梯级带，称为三方向小子域滤波法；肖锋等^[8]、蒋甫玉等^[9]、马国庆等^[10]、段晓旭^[11]、蔡钟等^[12]均对小子域滤波法的子域划分方式进行了改进，使其滤波输出更合理、稳定。夏玲燕等^[13]、许海红等^[14]提出小子域滤波与Tilt梯度或总水平导数联合处理的方法，取得了良好的边界识别效果。

根据前人资料，小子域滤波法自提出以来，一直有学者对该方法进行改进及应用研究，但对于改进后应用效果的系统整理和对比研究较少，这给使用者带来诸多不便。本文以5×5数据点窗口为例，采用小子域滤波法已有的10种类型分别对2组理论模型的正演重力异常数据进行处理，对比分析10种类型滤波处理的效果，讨论滤波窗口、迭代次数等参数对小子域滤波处理结果的影响，根据试验分析优选出滤波效果最佳的子域划分方式及滤波参数组合，并用于实际资料中，获得了较好的地质效果。

1.   小子域滤波法分类

小子域滤波法是基于滑动平均法原理进行改进的^[1]。它是将一个给定窗口分解成位于中心点不同侧面的子域^[1]，通过对各个“子域”进行计算，采用均方差^[1]或水平导数值^[10]判别比较，并根据判别结果将均方差或水平导数值最小者对应的子域的异常平均值记为给定窗口中心点的滤波输出值，依次类推完成的滤波过程。根据“窗口”划分方式及“子域”的形状，将目前使用的小子域滤波方法分为以下10种类型：“米”字形子域滤波法^{[1, 3]}、Tomita-Tsuji形子域滤波法^[15]、Nagao-Matsuyama形子域滤波法^[15]、对称形子域滤波法^[4]，“田”字形子域滤波法^[8]、“十”字形子域滤波法^[10]、“×”字形子域滤波法^[10]、平行四边形子域滤波法^[9]、滑动形子域滤波法^[11]、五六边形子域滤波法^[12]。

2.   模型试验及效果分析

为更好地显示10种小子域滤波方法的处理效果，设计A和B两组模型，以5×5数据点窗口为例分别采用上述10种方法对2组模型的理论重力异常进行处理，并将处理结果成图对比。模型A、B参数见表 1，模型分布情况及其理论重力异常见图 1。

表  1  理论模型参数

Table  1.  Parameters of the theoretical model

模型编号	地质体编号	x，y角点坐标/m	顶面埋深/m	底面埋深/m	密度差/(g· cm-3)
	A1	（3000，9000）; （5000，9000）; （5000，11000）; （3000，11000）	100	2100	0.7
A	A2	（9000，9000）; （11000，9000）; （11000，11000）; （9000，11000）	500	2500	0.7
	A3	（15000，9000）; （17000，9000）; （17000，11000）; （15000，11000）	1100	3100	0.7
	B1	（7125，7485）; （8050，7485）; （8050，8415）; （7125，8415）	100	200	0.7
	B2	（3855，3640）; （1960，2640）; （5160，855）; （6250，1320）	400	1000	0.7
B	B3	（1425，7292）; （1425，5278）; （3554，5278）; （2986，6920）	600	1000	0.7
	B4	（0，8750）; （3100，8150）; （5800，4150）; （10000，2733）; （10000，3200）; （6150，4500）; （3450，8500）; （0，9170）	200	400	0.7

下载: 导出CSV 

| 显示表格

图  1  模型AB及其理论重力异常（参考文献[9]修改）

a—模型A及其分布；b—模型A理论重力异常；c—模型B及其分布；d—模型B理论重力异常

Figure  1.  The model of AB and its theoretical gravity anomalies

下载: 全尺寸图片 幻灯片

2.1   模型A重力异常及其处理结果

如图 1-a所示，模型A由3个地质体A1、A2、A3组成，三者均为立方体且平面间距相同，但埋深不同（表 1）。图 1-b为模型A的重力异常平面图，成图采用的网格间距为80m × 80m，等值线间距为0.8mgal。由图 1-b可见，模型A的重力异常呈南北向对称分布，其中A1地质体的异常幅值最大，A2次之，A3最小，即由西向东异常幅值逐渐减小，且地质体的边界位置在图 1-b中不易直接识别，为此，采用前述10种方法对模型A的重力异常进行处理，分别得到10种处理结果（图 2）。

图  2  模型A重力异常的小子域滤波处理结果

Figure  2.  Results of model A by use of small subdomain filtering

下载: 全尺寸图片 幻灯片

由图 2可见，与图 1-b相比，“米”字形子域滤波结果将地质体A1、A2、A3的边界异常形态不同程度得到加强（图 2-a），但该方法对地质体拐角边界位置的反映失真，各地质体的拐角向地质体中心方向发生变形，且随着地质体埋深不断增大，这种变形程度逐渐加大。Tomita-Tsuji形子域滤波处理后A1、A2、A3地质体的边界形态均得到加强（图 2-b），且该方法对拐角位置的刻画较准确。Nagao-Matsuyama形子域滤波处理后也将地质体A1、A2、A3的边界特征增强（图 2-c），但该方法处理结果与“米”字形子域滤波处理结果类似，存在对地质体拐角边界位置刻画失真的问题，只是该方法处理后各地质体拐角向外部发生变形。对称形子域滤波处理结果（图 2-d）与图 1-b相比，地质体边界未见明显增强或变化。“田”字形子域滤波结果将地质体边界增强（图 2-e），且对地质体拐角边界位置的刻画准确，克服了“米”字形和Nagao-Matsuyama形子域滤波法的缺点，对地质体边界形态增强的整体效果较明确。“+”字形（图 2-f）和“×”字形（图 2-g）处理结果虽然也使地质体边界形态加强，且拐角边界没有发生明显的失真变形，但其边界的锐化程度不佳，如图 2-f中A2地质体边界形态不及图 2-e中A2地质体处理结果，且由图 2-f和图 2-g可见，“×”字形子域滤波的边界锐化效果不及“十”字形子域滤波的结果。平行四边形子域滤波法对地质体边界形态的增强作用较明显，但是也存在对地质体拐角边界位置刻画失真变形的问题（图 2-h）。滑动形子域滤波法（图 2-i）与“田”字形子域滤波法（图 2-e）处理结果相近，在增强地质体边界形态的同时对拐角边界的刻画也较准确。五六边形子域滤波法（图 2-j）与Nagao-Matsuyama形子域滤波（图 2-c）处理结果相似，虽然增强了地质体边界特征，但是也存在对地质体拐角边界位置刻画失真的问题。

2.2   模型B重力异常及其处理结果

为进一步探讨各种方法在不同边界延伸方向地质体中的处理效果，设计了模型B，如图 1-c所示，模型B为B1、B2、B3和B4四个地质体叠加成的一个组合模型。该模型中，4个地质体的形状大小和埋深各不相同，其中B1和B4地质体埋深较小，B2和B3地质体埋深较大（表 1），模型B引起的重力异常如图 1-d所示，成图采用的网格间距为50m× 50m，等值线间距为0.2mgal。由图 1-d可以看出，B1、B4地质体引起的重力异常幅值较小，B2、B3地质体引起的重力异常幅值较大，且B3与B4地质体邻近，它们引起的重力异常叠加，各地质体边界位置在图 1-d中不易确定。

图 3为模型B重力异常的小子域滤波处理结果，与图 1-d相比，“米”字形子域滤波处理结果（图 3-a）使各地质体边界特征得到增强，其中B4地质体边界增强效果最佳，B1、B3地质体均在拐角位置发生变形，B2地质体仅有局部边界变形。TomitaTsuji形子域滤波法处理后各地质体的边界特征均得到增强（图 3-b），且地质体拐角位置形态未发生失真。Nagao-Matsuyama形子域滤波处理结果（图 3-c）使B1、B4地质体边界增强效果明显，B2、B3地质体边界有一定程度的增强，但对边界位置的刻画不够精准，且地质体拐角边界向外变形。对称形子域滤波处理后（图 3-d），重力异常形态未见明显变化。“田”字形子域滤波处理后，4个地质体边界特征都得到增强（图 3-e），且该方法处理结果对地质体拐角的刻画较准确。“十”字形子域滤波处理后（图 3-f）各地质体边界特征得到一定程度的增强，但对地质体边界的锐化程度和平面位置反映的准确程度不够精确。由“×”字形子域滤波法处理结果（图 3-g）可见，地质体边界特征未得到明显增强，且地质体异常局部存在变形。平行四边形子域滤波处理后各地质体边界均不同程度得到增强（图 3-h），但B2、B3地质体拐角局部变形失真较明显。滑动形子域滤波处理结果（图 3-i）与“田”字形子域滤波法处理结果（图 3-e）基本一致，各地质体边界特征均得到增强，且地质体拐角位置未发生变形失真。五六边形子域滤波法（图 3-j）与Nagao-Matsuyama形子域滤波（图 3-c）处理结果相似，虽然将B1、B4地质体边界特征增强，但是也存在对B2、B3地质体拐角边界位置刻画失真的问题。

图  3  模型B重力异常的小子域滤波处理结果

Figure  3.  Results of model B by use of small subdomain filtering

下载: 全尺寸图片 幻灯片

由模型A、B重力异常小子域滤波5×5数据点窗口处理结果对比可见，“米”字形子域滤波法与平行四边形子域滤波法处理结果相似，其处理结果使地质体边界异常得到增强，但地质体拐角位置发生了变形失真。Nagao-Matsuyama形子域滤波法、“十”字形子域滤波法、“×”字形子域滤波法、五六边形子域滤波法处理结果类似，这4种方法处理后地质体边界增强作用不强或会出现异常边界局部模糊变形的现象。Tomita-Tsuji形子域滤波法、田字形子域滤波法和滑动形子域滤波法的处理效果相当，这3种方法对地质体边界增强和拐角位置的刻画作用均较好。

3.   滤波参数对比试验

3.1   对称形子域滤波法参数试验

根据5×5数据点窗口滤波处理结果可知，对称形子域滤波法对模型A、B理论重力异常的处理结果（图 2-d、图 3-d）与模型的原始异常形态异常特征等（图 1-b、d）无明显变化，为了进一步查明对称形子域滤波法的滤波效果，采用变化滤波窗口及迭代次数等参数对模型B的重力异常进行处理比较（图 4）。由图 4及图 3-d可见，随着滤波窗口的增大，模型B重力异常的增强效果（图 4-a、b）未见有明显改善；随着迭代次数的增加，其处理结果（图 4-c、d）也没有显著的变化，且B1、B2地质体的异常等值线有向中心变形的现象。

图  4  不同滤波窗口、迭代次数对称形子域滤波法对模型B重力异常的处理结果

a—9×9数据点滤波窗口、迭代1次；b—13×13数据点滤波窗口、迭代1次；c—9×9数据点滤波窗口、迭代2次；d—13×13数据点滤波窗口、迭代2次

Figure  4.  Results of model B by use of symmetrical shape small subdomain filtering

下载: 全尺寸图片 幻灯片

由上述试验可知，对称形子域滤波法在地质体边界异常增强方面效果不明显。根据马涛等^[4]的研究，在对各种观测的地球物理场数据进行滤波处理的过程中，常需要既保持原始场的正常变化梯度，又要能够消除随机噪声干扰。为此，对模型B重力异常（图 1-d）加入异常数值5%大小的随机噪声，然后采用5×5数据点窗口对称形子域滤波法对其进行滤波处理（图 5）。由图 5可见，对模型B重力异常（图 1-d）添加噪声后（图 5-a），模型B中B1、B2、B3、B4地质体异常边界特征被淹没，几乎不可辨认，由对称形子域滤波进行1次滤波处理后（图 5-b），大部分随机噪声被滤除掉，在1次滤波处理的基础上再进行1次滤波处理，得到对称形子域2次滤波处理结果（图 5-c），相比1次滤波处理结果，图 5-b中的随机噪声不仅得到进一步的滤除，而且数据变化趋势可以较清晰地分辨并与原始重力异常场（图 1-d）的变化趋势基本保持一致，各地质体引起的异常等值线形态没有发生变形或失真的现象。为深入讨论滤波窗口和迭代次数的影响情况，采用不同的处理参数（9×9、13×13数据点窗口，迭代1次、迭代2次等）对模型B含噪重力异常进行处理（图 6），由图 6及图 5-b、图 5-c可见，随着滤波窗口和迭代次数的增加，其滤波去噪效果未见明显提升，且滤波窗口和迭代次数太大时处理结果会出现变形失真的现象，如图 6-b、图 6-d中B2地质体异常向中心位置发生变形“偏心”现象。通过对各种滤波参数处理结果的对比，结合对称形子域2次（图 5-c）、3次（图 5-d）滤波处理结果可得，采用5×5数据点滤波窗口，进行2次迭代处理的对称形子域滤波处理结果即可达到处理目的，即能够在滤除随机噪声干扰的同时较好地保持原始场的正常变化梯度特征。

图  5  模型B含噪重力异常及对称形子域滤波法处理结果

a—对模型B重力异常加入异常数值5%的随机噪声后的异常图；b—对称形子域滤波法对图 5-a进行1次处理的结果；c—对称形子域滤波法对图 5-a进行2次处理的结果；d—对称形子域滤波法对图 5-a进行3次处理的结果

Figure  5.  Gravity anomaly of model B with noise and results by use of symmetrical shape small subdomain filtering

下载: 全尺寸图片 幻灯片

图  6  不同滤波窗口、迭代次数对称形子域滤波法对模型B含噪重力异常的处理结果

a—9×9数据点滤波窗口、迭代1次；b—13×13数据点滤波窗口、迭代1次；c—9×9数据点滤波窗口、迭代2次；d—13×13数据点滤波窗口、迭代2次

Figure  6.  Results of noised model B by use of symmetrical shape small subdomain filtering with different filter windows and different iterations

下载: 全尺寸图片 幻灯片

3.2   其他类型子域滤波参数试验

由5×5数据点窗口滤波处理结果（图 2、图 3）可知，“米”字形子域滤波法、Nagao-Matsuyama形子域滤波法、“十”字形子域滤波法、“×”字形子域滤波法、平行四边形子域滤波法、五六边形子域滤波法6种方法处理结果使地质体边界异常得到不同程度的增强，但也使地质体拐角位置发生了变形失真，或出现局部异常边界模糊变形或发散的现象；而Tomita-Tsuji形子域滤波法、“田”字形子域滤波法和滑动形子域滤波法3种方法对地质体边界增强和拐角位置的刻画作用均较好，因此选择这3种类型的子域滤波方法深入探究滤波参数对处理效果的影响。首先，变化滤波窗口（以9×9和13×13为例）对模型B重力异常（图 1-d）进行处理，由处理结果（图 3、图 7）可见，随着滤波窗口的增大，这3种方法对地质体异常边界的压缩程度逐渐增强，但是当窗口增加到一定尺寸时，增强效果不再有明显改善，且大窗口的处理结果会使异常等值线失真变形（图 7中B2、B3地质体异常边界），窗口越大变形情况越明显，相比而言，TomitaTsuji形子域滤波结果变形最大（图 7-b）、“田”字形子域滤波结果变形中等（图 7-d）、滑动形子域滤波结果变形最小（图 7-f）。由图 7-a、图 7-c、图 7-e与图 3-b、图 3-e和图 3-i对比可见，9×9数据点窗口比5×5数据点窗口滤波处理结果对地质体边界异常有一定的增强，但是增强效果不明显，且大窗口滤波处理结果存在使异常等值线产生失真变形的情况，因此选择5×5数据点滤波处理窗口即可。

图  7  不同滤波窗口3种类型小子域滤波法对模型B重力异常的处理结果

a—Tomita-Tsuji形子域9×9数据点滤波窗口、迭代1次；b—Tomita-Tsuji形子域13×13数据点滤波窗口、迭代1次；c— “田”字形子域9×9数据点滤波窗口、迭代1次；d— “田”字形子域13×13数据点滤波窗口、迭代1次；e—滑动形子域9×9数据点滤波窗口、迭代1次；

Figure  7.  Results of model B by use of 3 different shapes of small subdomain filtering with different filter windows

下载: 全尺寸图片 幻灯片

迭代次数也是影响滤波处理效果的一个因素，当滤波效果达不到预期目的时可以通过增加迭代次数进一步增强效果^[11]。根据上述不同滤波窗口的试验结果，在此选择5×5数据点窗口，变化迭代次数对模型B重力异常进行处理，由处理结果（图 3、图 8）可见，随着迭代次数的增加，各地质体边界异常等值线越来越密集，迭代次数越多效果越明显，但同时会产生一定程度的失真^[11]，且迭代次数越多失真情况越严重（如图 8中B2、B3地质体边界异常特征尤为明显）。

图  8  不同迭代次数3种类型子域滤波法对模型B重力异常的处理结果

a—Tomita-Tsuji形子域5×5数据点滤波窗口、迭代2次；b—Tomita-Tsuji形子域5×5数据点滤波窗口、迭代3次；c— “田”字形子域5×5数据点滤波窗口、迭代2次；d— “田”字形子域5×5数据点滤波窗口、迭代3次；e—滑动形子域5×5数据点滤波窗口、迭代2次；f—滑动形子域5×5数据点滤波窗口、迭代3次

Figure  8.  Results of model B by use of 3 different shapes of small subdomain filtering with different iterations

下载: 全尺寸图片 幻灯片

依据重力异常研究或推断地质构造和地质体边界，即对异常的解释一般是从“读图”或异常识别开始的^[16]，因而对于处理结果中出现的异常等值线失真或变形会影响解释精度，另一方面当这种失真变形较严重时甚至会误导异常的解释结果。由图 8、图 3对比可见，迭代次数越多，对地质体边界异常的增强效果越明显，但同时会产生异常等值线失真变形的情况；再由图 3与图 1-d对比可知，迭代1次的处理结果能够将地质体的边界异常得到增强，可以较清晰地反映出异常特征且没有出现明显的失真或变形，因此迭代次数选择1次即可。

4.   应用实例

模型A、B重力异常的处理对比结果表明，小子域滤波法在位场数据处理中有明确的有效性：一方面，经过对称形子域滤波法的2次迭代处理后，原始异常中的随机噪声干扰基本可以被滤除掉，而且原始场的正常梯度变化特征保持不变；另一方面，因滑动形子域滤波的子域划分方式更合理，其处理结果对地质体边界位置的刻画程度更准确。下面以SH区块重力异常为例说明小子域滤波法在资料处理中的实用性。

图 9-a为SH区块布格重力异常，根据地质资料可知，该区重力高异常和重力低异常分别是凸起和凹陷构造的反映，为查明该区断裂体系，落实地质体构造边界位置，采用小子域滤波法对布格重力异常进行处理；选择5×5数据点窗口，首先，采用对称形子域滤波法对布格重力异常（图 9-a）进行2次迭代处理，根据处理结果（图 9-b）可见，原始布格重力异常中的随机噪声干扰基本被消除，异常等值线更加光滑，且图 9-b较好地保持了原始重力异常场的变化特征。然后使用滑动形子域滤波法对图 9-b进行处理（1次迭代），与原始布格重力异常相比，滑动形子域滤波法使区内重力高异常与重力低异常之间的界线得到增强、重力异常梯级带得到压缩（图 9-c），地质构造的边界位置更明显。最后，求取滑动形子域滤波法处理结果（图 9-c）的总水平导数，并以伪彩色图像显示其结果（图 9-d），由图 9-d可见，小子域滤波法与总水平导数联合处理后，将地质体界线异常等值线转换为重力异常极大值连线，与重力异常平面等值线图相比，伪彩色图像显示的方法更形象、更直观地表达了研究区的断裂格局及构造界线位置。

图  9  SH区块重力异常及小子域滤波法处理结果

a—SH区块布格重力异常；b—对称形子域滤波法对图 9-a进行2次迭代处理的结果；c—滑动形子域滤波法对图 9-b处理的结果；d—对图 9-c求取总水平导数的结果

Figure  9.  Gravity anomaly of block SH and the results by use of small subdomain filtering

下载: 全尺寸图片 幻灯片

根据小子域滤波及与总水平导数的联合处理，刻画了SH区块平面断裂位置，推测了SH区块平面断裂体系的基本构造格架（图 10-a）。由图 10-a可见，该区断裂较发育，其中6条主要断裂是控制凸起和凹陷构造的边界断裂，次级断裂若干且其控制了凸起/凹陷构造内部的局部构造。为了验证SH区块断裂解释结果的有效性，收集了该区完成的2条地震剖面和1条电法剖面进行对比分析，剖面位置见图 10-a。图 10-b是A-A′地震剖面及其确定的断裂，该剖面在图 10-a中与F₃、F₄断裂相交，通过对比地震剖面中断裂的产状，以及这2条断裂在重力异常中的位置，可以确认A-A′地震剖面解释的2条断裂，即为图 10-a中划分的F3、F4断裂，它们是控制该区中央凹陷的边界断裂。图 10-c是B-B′地震剖面及确定的断裂，在图 10-a中与F₃、F₄断裂相交，该剖面与A-A′剖面类似，通过对比可见，该剖面确定的2条断裂即为重力资料推断的F₃、F₄断裂。图 10-d是C′-C电阻率断面图，由断面图可见该剖面上存在明显的电性差异，根据电阻率等值线密集的垂向高低阻过渡带推测了该剖面的主要断裂，与图 10-a对比，可以确定该2条断裂即为F₁、F₂断裂，重力和地震、电法资料解释结果一致。

图  10  SH区块断裂体系及利用地震、电法剖面确定的断裂

a—断裂平面分布；b—A-A′地震剖面；c—B-B′地震剖面；d—C′-C电法视电阻率剖面

Figure  10.  Fault structures of block SH and the evidence by seismic and electrical sections

下载: 全尺寸图片 幻灯片

5.   结论

（1）理论模型和实际资料处理结果表明，“米”字形子域滤波法与平行四边形子域滤波法使地质体边界得到增强，但地质体拐角位置会发生变形失真。Nagao-Matsuyama形子域滤波法、“十”字形子域滤波法、“×”字形子域滤波法和五六边形子域滤波法处理后地质体边界增强作用不明显，且会出现异常边界局部模糊变形的现象。Tomita-Tsuji形子域滤波法、“田”字形子域滤波法和滑动形子域滤波法对地质体边界增强和拐角位置的刻画作用均较好，其中，滑动形子域滤波法增强效果最好。

（2）根据滤波法参数对比试验可见，对称形子域滤波法的处理效果主要体现在，能够滤除随机噪声的同时较好地保持了原始异常场的梯度变化特征，滤波窗口选择5×5数据点，迭代次数选择2次较好；滑动形子域滤波主要表现在对地质体边界异常的增强、锐化等方面，这种情况下滤波窗口选择5×5数据点，迭代次数选择1次即可。