改进的均衡滤波器在位场数据边界识别中的应用

马国庆，黄大年，于 平，李丽丽

吉林大学地球探测科学与技术学院，长春 130021

1 引 言

位场主要是指重力和磁力场，重磁勘探具有快速、经济、范围大的优点，是进行构造划分和异常圈定快速而有效的方法，但是原始异常与地质体边界之间不存在良好的对应关系，所以直接根据重磁异常不易识别出地质体的边界.重磁异常水平导数的极大值和垂直导数的零值与地质体边界相对应，该性质常被用来进行地质体边界的识别.垂直导数（1936）是最早被用来进行地质体边界识别的方法［1-2］；后来人们发现水平导数的极值位于密度与磁化率发生突变的位置，该方法被证明是一种非常有效的边界识别工具［3-4］；解析信号的极大值与磁性体的边界也存在良好的对应关系［5］，Hsu等利用高阶解析信号使磁性体的边界更加清晰［6］，但上述边界识别方法均存在不能同时显示浅部和深部地质体边界的缺点［7］，这是由于导数随地质体埋深的衰减速度较快，因此更易凸显出浅部地质体的效应.为了解决这一问题，人们开始致力于均衡滤波器的研究.Miller和Singh在1994年提出了第一个均衡边界识别滤波器—tilt angle［8］，该方法能很好地均衡强弱异常之间的幅值，但不是一个很好的边界识别工具；Rajagopalan和 Milligan在1995年提出利用自动增益控制来进行磁源体边界的识别［9］，该方法的识别结果依赖于窗口的尺寸，灵活性较差；Verduzco等在2004年提出利用tilt angle的总水平导数（THDR）进行地质体边界的识别［10］，取得不错的效果；Wijns等在2005年利用总水平导数与解析信号的比值（theta map）来进行边界识别［11］，取得了很好的效果；Cooper和Cowan在2006总结了多种不同形式的边界识别工具［12］，并提出了其它形式的倾斜角进行异常体边界的识别；Cooper和Cowan 2008年利用水平与垂直导数的均方差进行异常体边界的识别［13］，该方法的结果与窗口的尺寸有关；后来人们利用异常的希尔伯特变换来进行地质体边界的识别［7，14］，该方法可有效地降低噪声的干扰，但所识别出的地质体边界比较模糊；马国庆等在2011年提出利用总水平导数与垂直导数的相关系数进行地质体边界的识别［15］，取得了不错的应用效果，但是该方法对小范围异常体的识别效果不是很好.Ma和Li在2012年采用归一化总水平导数法进行地质体边界的识别［16］，有效地提高了输出结果的稳定性，但识别出的边界存在一定的发散.

本文提出增强型均衡滤波器来进行地质体边界的识别，该滤波器是由不同阶导数构建的，其能更加清晰和准确地识别出地质体的边界.

2 增强型均衡滤波器

现有的边界识别滤波器均是同阶导数之间的组合，增强型滤波器是利用不同阶导数之间的非线性组合，其输出结果更加准确和清晰.本文给出了三种不同形式的增强型均衡滤波器.

2.1 增强型tilt angle

Miller和Singh在1994年提出采用tilt angle进行地质体边界的识别，其表达式为：

其中，f为原始重力或磁异常.tilt angle能很好地均衡不同深度异常之间的幅度，但是该方法并不能很好地进行地质体边界的识别［15］，因此Verduzco等在2004年提出利用tilt angle的总水平导数（THDR）进行地质体的边界识别工作.

本文提出了另一种形式的tilt angle，称之为增强型tilt angle（ETA），其表达式为：

其中，n表示垂直导数的阶数，也以此表示滤波器的阶数，

系数cn是为了使公式（3）满足数学运算法则，并可均衡不同阶导数之间的强度.增强型tilt angle为垂直导数与其水平导数的组合，其分子与分母是不同阶的，这是其与常规tilt angle的主要区别，增强型tilt angle的最大值与地质体的边界相对应.

2.2 增强型theta map

Wijns等在2005年提出theta map来进行地质体边界的识别，其具体公式为：

增强型theta map（ETM）的具体表达式为：

增强型theta map的最大值也与地质体边界相对应.

2.3 增强型cosine function

第三种是利用不同阶导数之间比值的余弦函数来进行异常体边界的识别，称之为增强型cosine function（ECF），其具体表达式为：

其中，R表示取实部运算，这是因为反余弦函数的运算结果中会出现虚数.增强型cosine function的最大值与地质体边界相对应.为了显示增强型滤波器的可行性，给出一系列重力异常剖面（图1）.

图1a为模型的原始异常，2个模型埋深分别为12m和15m.图1b表示原始异常的THDR计算结果，该方法不能清晰地给出深部地质体的边界，且会产生多余的干扰异常；图1c表示异常的theta map计算结果，该方法能同时显示浅部与深部异常的边界，但是所识别出的边界较发散；图1d—1f分别表示ECF1、ETA1和ETM1的计算结果；图1g—1i分别表示ECF2、ETA2和ETM2的计算结果；一阶和二阶增强型均衡滤波器均能很清晰地识别出地质体的边界，边界非常集中.图1j—1l分别表示ECF3、ETA3和ETM3的计算结果，三阶增强型均衡滤波器的识别结果产生了多余的峰值，不利于识别出地质体的边界，无法对异常进行解释.因此仅可以利用一阶和二阶增强型均衡滤波器来进行地质体边界的识别.

从增强型均衡滤波器的表达式中可以看出，二阶增强型滤波器需要计算二阶垂直导数，二阶垂直导数的计算会明显增大噪声的干扰.为了减小噪声的干扰，采用Laplace方程［17］来完成二阶垂直导数的计算：

图1 不同方法边界识别结果的剖面图（a）原始重力异常；（b）异常的THDR结果；（c）异常的theta map结果；（d）异常的ETA1 结果；（e）异常的ETM1 结果；（f）异常的ECF1 结果；（g）异常的ETA2 结果；（h）异常的ETM2 结果；（i）异常的ECF2 结果；（j）异常的ETA3 结果；（k）异常的ETM3 结果；（l）异常的ECF3 结果.Fig.1 Profiles showing different edge identifications from original gravity data（a）Original gravity anomaly；（b）THDR of the data in（a）；（c）theta map of the data in（a）；（d）ETA1of the data in（a）；（e）ETM1of thedata in（a）；（f）ECF1of the data in（a）；（g）ETA2of the data in（a）；（h）ETM2of the data in（a）；（i）ECF2of the data in（a）；（j）ETA3of the data in（a）；（k）ETM3of the data in（a）；（l）ECF3of the data in（a）.

从（7）式中可以看出，二阶垂直导数可以采用两个二阶水平导数之和来完成，水平导数采用空间域方法来计算得到，不会增大噪声的干扰.公式中所涉及的其它低阶垂直导数的计算采用傅里叶变换来完成.

3 理论模型试验

为了试验增强型均衡滤波器的可行性，建立如下地质模型：地下布设了2个中心埋深分别为15m和20m、与围岩密度差为ρ＝2g／cm3、半边长为10m的棱柱体.图2a为模型理论重力异常，白虚线标识地质体的真实水平位置.分别利用上述的方法对重力异常进行处理（图2）.

图2b表示异常THDR的结果，THDR法不能很清晰地给出地质体的边界信息；图2c表示异常theta map的结果，theta map能够清晰地识别出地质体的边界，但是边界比较发散，导致两个地质体的边界已经相连.高阶导数对地质体的边界具有更高的分辨率，因此可以通过增加导数的阶次来使地质体边界更加清晰，考虑采用一阶垂直导数构建theta map滤波器，其具体表达式为：

图2 不同方法的边界识别结果（a）原始重力异常；（b）异常的THDR结果；（c）异常的theta map结果；（d）异常的一阶垂直导数的theta map结果；（e）异常的ETA1 结果；（f）异常的ETM1 结果；（g）异常的ECF1 结果；（h）异常的ETA2 结果；（i）异常的ETM2 结果；（j）异常的ECF2 结果.Fig.2 Identifications of edges by different methods（a）Original gravity anomaly；（b）THDR of the data in（a）；（c）theta map of the data in（a）；（d）theta map of the first vertical derivative of the data in（a）；（e）ETA1of the data in（a）；（f）ETM1of the data in（a）；（g）ECF1of the data in（a）；（h）ETA2of the data in（a）；（i）ETM2of the data in（a）；（j）ECF2of the data in（a）.

图2d表示一阶垂直导数的theta map（Theta2），该方法使边界更加收敛，但是在高值外侧存在明显的低值，为异常的解释工作带来了困难.一阶垂直导数与其水平导数组成的边界识别滤波器（图2e，2f，2g）能很好地完成边界识别工作，且未产生干扰异常，识别出的边界很清晰，采用本文方法构建的边界识别工具避免了采用高阶导数直接构建边界识别滤波器所带来的干扰.二阶垂直导数与其水平导数组成的边界识别工具（图2h，2i，2j）能更加准确地描述地质体的水平位置，且也不存在干扰异常.从该试验中可以看出，增强型均衡滤波器能准确地识别出地质体的水平位置，识别出的边界非常清晰，能很好地完成异常的解释工作.

为了进一步试验方法的有效性，建立如下较为复杂的模型：重力异常由四个长方体组成，埋深分别为20m（1），30m（2），30m（3），5m（4）.图3a为理论重力异常，图中白虚线表示地质体的真实水平位置，分别利用上述边界识别方法对该异常进行处理（图3）.

图3b表示异常THDR的结果，该方法能大致地给出地质体的边界信息，但是埋深较深的地质体边界不是很清晰；图3c表示异常theta map的结果，该方法能给出较大的地质体边界，而不能给出较小地质体边界，且边界比较发散.一阶增强型均衡滤波器（图3d，3e，3f）能很清晰地给出较大的地质体的边界，而对于较小的地质体的边界比较模糊.二阶增强型均衡滤波器（图3g，3h，3i）能更加准确地给出地质体的边界信息，且较小的地质体的边界反映也十分清晰.从图3中可以看出，本文方法能更好地完成边界识别工作.从图3d，3e，3f与3g，3h，3i的对应中可以看出，随着所使用的导数阶次的增加能识别出更多的细节信息.

为了试验方法的稳定性，在图3a所示的重力异常中加入信噪比为50的高斯白噪声，图4a为加入噪声后的重力异常，利用边界识别方法对该异常进行处理（图4）.

图3 不同方法的边界识别结果（a）原始重力异常；（b）异常的THDR结果；（c）异常的theta map结果；（d）异常的ETA1 结果；（e）异常的ETM1 结果；（f）异常的ECF1 结果；（g）异常的ETA2 结果；（h）异常的ETM2结果；（i）异常的ECF2 结果.Fig.3 Identifications of edges by different methods（a）Original gravity anomaly；（b）THDR of the data in（a）；（c）theta map of the data in（a）；（d）ETA1of the data in（a）；（e）ETM1of the data in（a）；（f）ECF1of the data in（a）；（g）ETA2of the data in（a）；（h）ETM2of the data in（a）；（i）ECF2of the data in（a）.

图4 加入噪声后不同方法的边界识别结果（a）原始重力异常；（b）异常的THDR结果；（c）异常的theta map结果；（d）异常的ETA1 结果；（e）异常的ETM1 结果；（f）异常的ECF1 结果；（g）由方程（7）计算的异常的ETA2结果；（h）由方程（7）计算的异常的ETM2结果；（i）由方程（7）计算的异常的ECF2结果；（j）由Fourier变换计算的异常的ETA2结果；（k）由Fourier变换计算的异常的ETM2结果；（l）由Fourier变换计算的异常的ECF2结果.Fig.4 Identifications of edges by different methods after adding noise（a）Original gravity anomaly；（b）THDR of the data in（a）；（c）theta map of the data in（a）；（d）ETA1of the data in（a）；（e）ETM1of the data in（a）；（f）ECF1of the data in（a）；（g）ETA2of the data in（a）computed by the Eq.（7）；（h）ETM2of the data in（a）computed by the Eq.（7）；（i）ECF2of the data in（a）computed by the Eq.（7）；（j）ETA2of the data in（a）computed by Fourier transformation；（k）ETM2of the data in（a）computed by Fourier transformation；（l）ECF2of the data in（a）computed by Fourier transformation.

图4b表示异常THDR的结果，由于噪声的干扰，该方法已经无法给出地质体的边界；图4c表示异常theta map的结果，该方法依旧可以给出地质体的边界.一阶增强型均衡滤波器（图4d，4e，4f）能很清晰地给出地质体的边界，受噪声影响较小.图4g—4i为采用方程（7）计算得到的异常的ETA2、ETM2和ECF2结果，该方法计算出的二阶增强型均衡滤波器仍能很清晰地给出地质体的边界.图4j—4l为采用常规Fourier变换计算得到的ETA2、ETM2和ECF2结果，但是采用该方法计算出的二阶增强型均衡滤波器根本无法给出地质体的边界，边界已经被损坏.因此采用方程（7）进行高阶垂直导数的计算能有效地提高输出结果的稳定性.从对比中可以看出，本文提出的增强型均衡滤波器的稳定性不比现有的仅由一阶导数组成的边界识别滤波器差，即使在存在噪声的情况下依旧能得到稳定的结果.

在应用本文方法进行磁异常处理前，要对磁异常进行化极处理，因为磁数据及其导数均受磁化方向的干扰［18］，因此采用化极异常进行解释所获得的结果将更加准确.

4 实际应用效果

为了验证方法在实际中的应用效果，对四川地区重力异常和朱日和地区磁异常进行处理.图5a为四川地区重力异常，重力数据采自国家测绘总局编绘的1∶100万布格重力异常图，并用黑虚线在图中标识出已知断裂的水平位置.利用上述方法对重力数据进行处理，来得到该地区的线性构造和地层之间的界线（图5）.

图5 同图3，但为四川地区重力异常边界识别结果Fig.5 Same as Fig.3，but for edge identification of gravity anomalies in Sichuan

图6 同图3，但为朱日和地区磁力异常边界识别结果Fig.6 Same as Fig.3，but for edge identification of magnetic anomalies in the Zhurihe area

THDR法（图5b）与theta map法（图5c）划分断裂的能力较差，仅能给出四川地区部分已知断裂的位置且不是十分清晰.一阶增强型均衡滤波器（图5d，5e，5f）能清晰地给出断裂的分布特征，根据该结果可以很容易地划分出断裂的位置及走势，且与已知断裂对应较好.二阶增强型均衡滤波器的结果（图5g，5h，5i）与已知断裂也存在较好的对应，此外还发现了更多的小的断裂及小范围的异常，能得到更多的细节信息.

下面试验本文方法在磁异常数据处理中的作用，对中国西北部朱日和地区的磁异常进行处理（图6）.图6a为该地区化极后磁异常，并根据前期的研究工作划定出已知成矿带的大致范围.

从图6中可以看出，增强型均衡滤波器的边界识别效果要明显好于其它边界识别技术，能更加清晰地显示出成矿带的界线，且与前期划定的成矿带范围拟合较好，为进一步开发提供了有力的保证.二阶增强型均衡滤波器能发现更多小范围的异常，很好地描述了异常的细节特征.

5 结 论

本文提出一种新的边界识别滤波器，称为增强型均衡滤波器，其利用不同阶导数之间的组合来进行地质体边界的识别.通过理论模型试验证明增强型均衡滤波器相对于现有的边界识别方法具有更高的分辨率，能更加清晰和准确地识别出地质体的边界，且随着增强型均衡滤波器阶次的增加会得到更多的细节信息.在计算过程中为了降低二阶垂直导数计算所带来的噪声干扰，引入了一种计算二阶垂直导数的稳定算法，使滤波器的输出结果更加稳定.最后将其应用于四川盆地重力异常及朱日和地区磁异常的解释中，获得了区域地质构造的水平位置，并发现了更多的小范围构造.增强型均衡滤波器是一种非常有效的边界识别工具，具有良好的应用前景.

（References）

［1］Evjen H M.The place of the vertical gradient in gravitational interpretations.Geophysics，1936，1（1）：127－136.

［2］Hood P J，Teskey D J.Aeromagnetic gradiometer program of the Geological Survey of Canada.Geophysics，1989，54（8）：1012－1022.

［3］Cordell L.Gravimetric expression of Graben faulting in Santa Fe country and the Espanola Basin，New Mexico.New Mexico Geol.Soc.Guidebook，30thField Conf.，1979：59－64.

［4］Cordell L，Grauch V J S.Mapping basement magnetization zones from aeromagnetic data in the San Juan basin，New Mexico.／／Hinzc W J ed.The Utility of Regional Gravity and Magnetic Anomaly.Society of Exploration Geophysicists，1985：181－197.

［5］Roest W R，Verhoef J，Pilkington M.Magnetic interpretation using the 3－D analytic signal.Geophysics，1992，57（1）：116－125.

［6］Hsu S，Sibuet J C，Shyu C.High－resolution detection of geologic boundaries from potential field anomalies：an enhanced analytic signal technique.Geophysics，1996，61（2）：373－386.

［7］Cooper G R J.Balancing images of potential－field data.Geophysics，2009，74（3）：L17－L20.

［8］Miller H G，Singh V.Potential field tilt—a new concept for location of potential field sources.Journal of Applied Geophysics，1994，32（2－3）：213－217.

［9］Rajagopalan S，Milligan P.Image enhancement of aeromagnetic data using automatic gain control.Exploration Geophysics，1995，25（4）：173－178.

［10］Verduzco B，Fairhead J D，Green C M，et al.New insights into magnetic derivatives for structural mapping.The Leading Edge，2004，23（2）：116－119.

［11］Wijns C，Perez C，Kowalczyk P.Theta map：edge detection in magnetic data.Geophysics，2005，70（4）：L39－L43.

［12］Cooper G R J，Cowan D R.Enhancing potential field data using filters based on the local phase.Computers ＆Geosciences，2006，32（10）：1585－1591.

［13］Cooper G R J，Cowan D R.Edge enhancement of potentialfield data using normalized statistics.Geophysics，2008，73（3）：H1－H4.

［14］骆遥，王明，罗锋等.重磁场二维希尔伯特变换——直接解析信号解释方法.地球物理学报，2011，54（7）：1912－1920.Luo Y，Wang M，Luo F，et al.Direct analytic signal interpretation of potential field data using 2－D Hilbert transform.Chinese J.Geophys.（in Chinese），2011，54（7）：1912－1920.

［15］马国庆，杜晓娟，李丽丽.利用水平与垂直导数的相关系数进行位场数据的边界识别.吉林大学学报（地球科学版），2011，41（S1）：345－348.Ma G Q，Du X J，Li L L.Edge detection of potential field data using correlation coefficients of horizontal and vertical derivatives.Journal of Jilin University （Earth Science Edition）（in Chinese），2011，41（S1）：345－348.

［16］Ma G，Li L.Edge detection in potential fields with the normalized total horizontal derivative.Computers ＆Geosciences，2012，41：83－87.

［17］Fedi M，Florio G.Detection of potential fields source boundaries by enhanced horizontal derivative method.Geophysical Prospecting，2001，49（1）：40－58.

［18］Li X.Understanding 3Danalytic signal amplitude.Geophysics，2006，71（2）：L13－L16.