唐新功，侯帅，李丹丹，2

1.油气资源与勘探技术教育部重点实验室(长江大学)，湖北 武汉 430100

2.武汉大学中国南极测绘研究中心，湖北 武汉 430079

作为最成熟的地球物理位场勘探方法之一，重磁勘探在油气勘探、固体矿产勘探、地下水、工程物探、军事、大地测量、深部构造以及地震预报中都发挥着不可或缺的重要作用[1]。在重磁勘探中，实测重磁数据包含了地下不同规模、不同深度各种地质异常体组合产生的叠加异常，为了对目标地质异常体进行勘探，必须在叠加重磁异常中提取出目标异常。因此，重磁资料处理中异常分离是不可或缺的一步，其中使用重磁数据进行异常边界识别也是位场资料处理的重要过程，能否分辨出地下断裂构造或密度体边界是衡量重磁勘探方法应用效果的重要指标。随着勘探精度的不断提高，直接从布格重力异常或磁异常中得到地质边界信息的难度越来越大，因此边界增强或突出技术显得尤为重要。边界增强方法的基本原理是利用导数的零值点或极值点位置来判定地质异常体的边界，在地质异常体边界附近，重磁异常的变化率或梯度较大。边界识别方法众多，主要可分为数理统计和数值计算两大类。

数理统计类重磁边界识别方法主要包括小子域滤波和归一化标准差两种。杨高印[2]最早提出了小子域滤波法，该方法是在滑动平均法的基础上进行改进的，其基本原理是将滑动平均的窗口划分成位于中心的不同侧面的8个小子域，然后将均方差最小的子域的平均值作为窗口中心的值。段晓旭[3]在小子域滤波的基础上提出了滑动子域滤波法，改进了拐点处的数据处理效果，但是存在的折线走样现象并没有得到很好的解决。另一种数理统计类方法是归一化标准差方法，该方法通过计算一个滑动窗口内的垂向一阶导数的标准差与X、Y、Z方向的标准差之和的比值作为窗口中心的值。数理统计类方法的优点是可以根据实际情况的不同，选择不同的窗口大小来消除噪音干扰，提高边界识别的分辨率。

相对于数理统计方法，数值计算类方法更多地被用于重磁资料的处理中。该类方法主要包括垂向导数法、总水平方向导数法、解析信号法、倾斜角法和Theta图法等[4-9]。数值计算类方法的边界识别原理主要是通过对重磁异常进行求取垂向导数、水平导数等处理之后再根据结果的极大值位置或者零值位置来识别出地质异常体的边界。此外，数值计算类方法还衍生出了一些处理方法，主要包括解析信号倾斜角法、广义导数算子法、增强型均衡滤波器法、平面全张量梯度角度法和改进二阶倾斜角法等[10-14]。

上述边界识别方法的适用条件不同，处理效果也不同，了解其适用条件及如何选择哪种方法或哪些方法更加合理至关重要。目前，针对不同处理方法适用条件及其优缺点的详细探讨和对比的文献较少。因此，笔者使用自行编制的重磁处理软件对上述边界识别方法的优缺点进行了详细的研究和对比，对于帮助提高重磁勘探的应用效果具有重要的理论和实用价值。

1 边界识别方法

1)一阶导数类边界识别方法。主要包括垂向一阶导数法、总水平方向导数法、解析信号法、倾斜角法、Theta图法、归一化标准差法等几种处理方法，也是重磁勘探中最基本和最常用的方法。其中，垂向一阶导数法、总水平方向梯度法以及解析信号法是对原函数直接求取一阶导数，而倾斜角法和Theta图法则是在前3种方法的基础上经过求取比值而获得的[4-10]。

2)二阶导数类边界识别方法。主要包括了解析信号倾斜角法、广义导数算子法、增强型均衡滤波器法、平面全张量梯度角度法、改进的二阶倾斜角法等[11-15]。解析信号倾斜法与倾斜角法利用零值位置识别边界不同的是，该方法利用极大值位置来识别地质异常体的边界位置。

2 软件编制

根据重磁正演理论方法，使用C#语言编写了三维可视化重磁正演程序。该程序能够计算球体和直立长方体的重磁正演异常，也可以进行重磁数据的初步处理。在程序启动界面可以选择计算重力异常或磁异常，然后设置测网大小和点距。

为了验证本软件的正确性，设计并计算了一个球体和一个长方体模型的磁异常响应。球体的半径是200 m，球心坐标X是300 m，Y是0 m，Z是500 m，磁化强度1 A/m，磁化倾角45°，磁化偏角40°。长方体的长为500 m，宽为300 m，高为200 m，中心坐标X是-300 m，Y是200 m，Z是500 m，磁化强度是1 A/m，磁化倾角45°，磁化偏角20°。测量区域网格大小X方向坐标是-1 000～1 000 m，Y方向的坐标是-1 000～1 000 m，网格点数X方向是201个，Y方向是201个，计算结果如图1(a)所示。为了进一步验证本软件的正确性，使用RGIS商业软件的重磁模型库模块计算了相同模型的磁异常正演，并对比了两者的计算结果，误差如图1(b)所示。可以看出，两个软件的计算结果非常接近，误差极小，表明了本软件的正确性。

图1 作者编制软件的磁异常正演计算结果及其与RGIS商业软件计算结果的误差

3 不同地质异常体边界识别方法的计算与对比

为了检验上述不同地质异常体边界识别方法在不同地质条件下对边界增强效果的优劣及不同特征，设置了3种理论重磁模型(地下单一地质异常体模型，地下两个规模相同但埋深不同的地质异常体模型，添加随机噪音的单一地质异常体模型)，分别使用不同地质异常体边界识别方法对得到的重磁响应进行处理和对比。

3.1 地下单一地质异常体模型

在均匀半空间中设置一个长、宽、高分别为800 m、800 m和400 m的三维地质异常体模型，顶面埋深100 m、剩余密度为300 kg/m3；磁化倾角45°，磁化偏角0°，测量区域网格大小X方向坐标为-1 000～1 000 m，Y方向坐标为-1 000～1 000 m，网格点数X方向为201个，Y方向为201个(见图2(a))。通过计算其重磁响应以及使用不同地质异常体边界识别方法处理后的效果对比不同方法的边界识别能力的优劣。

图2 均匀半空间中地下单一地质异常体模型产生的重力异常及使用不同地质异常体识别方法处理后的结果图

图2是均匀半空间中地下单一地质异常体模型产生的重力异常及使用不同地质异常体识别方法处理后的结果图。由图2可知，垂向一阶导数法利用其零值位置范围可以确定出模型的边界位置，只是比模型的范围略大一点(见图2(b))；解析信号法的极大值位置集中在模型的中部，基本能够反映模型的位置，但是边界识别不清晰(见图2(c))；总水平方向导数法的极大值位置与模型边界位置基本重合，基本可以识别出地质异常体边界，但在角落处误差较大(见图2(d))；倾斜角法和垂向一阶导数法一样都利用零值位置识别边界，其效果与垂向一阶导数法类似，但是识别出的边界有偏差，边界识别不清晰(见图2(e))；Theta图法利用极大值位置识别边界，其极大值的位置和垂向一阶导数法零值位置重合，但是极大值的峰值较宽，对边界识别的分辨率不足(见图2(f))；归一化标准差法与总水平方向导数法结果相似，其极大值与模型边界基本重合，对边界识别的分辨率较好(见图3(g))；解析信号倾斜角法的识别结果与解析信号法基本一致(见图3(h))；广义导数算子法计算时使用的仰角是90°此时方位角未发挥作用，利用其零值位置识别的边界效果和垂向一阶导数法的结果相同(见图2(i))；增强型倾斜角法的极大值位置与模型边界重合，且极大值边界尖锐，有利于识别边界(见图2(j))；平面全张量倾斜角法的极大值位置较模型边界更大一些，识别结果没有总水平方向导数法好(见图2(k))；改进二阶倾斜角法的计算结果非常锐利清晰，能够最佳地突出地质异常体的边界位置(见图2(l))。

图3 均匀半空间中单一地质异常体模型产生的磁异常及使用不同地质异常体识别方法处理后的结果图

通过图2不同地质异常体识别方法的处理效果可以看出，对目标体边界位置识别最准确和最清晰的是改进二阶倾斜角法和增强型倾斜角法；其次是总水平方向导数法与归一化标准差法，二者基本能够识别出模型的边界位置；平面全张量倾斜角法相比总水平方向导数法的边界位置略微向外偏移，对边界识别效果不佳；垂向一阶导数法、倾斜角法、Theta图法和广义导数算子法的边界识别结果相似，虽然基本能够识别模型的边界位置，但是得到的边界位置都与模型边界有偏差，识别效果一般；解析信号法和解析信号倾斜角法对边界的识别效果最差。

图3是均匀半空间中单一地质异常体模型产生的磁异常及使用不同地质异常体识别方法处理后的结果图。为了考察倾斜磁化的影响，这里均不做化极处理。由图3可知，垂向一阶导数的零值位置与模型边界位置相差较大，不利于边界识别(见图3(b))；解析信号法的极大值与模型边界位置基本重合，能够较为准确地指示边界位置，但是不够完整，只能识别与磁化方向垂直的边界(见图3(c))；总水平导数法的极大值存在双峰现象，也能在一定程度上识别出模型边界位置(见图3(d))；倾斜角法的零值位置和垂向一阶导数法的零值位置重合，同样与模型边界位置相差较大(见图3(e))；Theta图法的极大值位置与垂向一阶导数法的零值位置基本相同，识别效果较好，但在边界识别中具有一定误差(见图3(f))；归一化标准差法的极大值位置识别出的边界相对完整，与磁化方向平行的边界位置能够准确地识别出来，但与磁化方向垂直的边界识别会出现双峰现象(见图3(g))；解析信号倾斜角法的极大值位置与模型边界位置基本重合，能够准确且较完整地识别出边界(见图3(h))；广义导数算子法的识别结果与垂向一阶导数基本相同，无法准确指示出模型边界位置(见图3(i))；增强型倾斜角法的边界识别结果比较准确，但是存在较多的假边界(见图3(j))；平面全张量倾斜角法也能准确的指示出部分边界的位置，但是存在较多的假边界(见图3(k))；改进二阶倾斜角法的极大值位置清晰，边界识别结果也较准确，但是同样也存在多余的假边界现象(见图3(l))。

综上所述，使用磁异常数据进行地下地质异常体边界识别时，解析信号倾斜角法的边界识别效果最好；其次是解析信号法和Theta图法，受到倾斜磁化的影响也较小；总水平导数法的识别结果较好，能够显示边界位置，但是由于倾斜磁化的影响，对边界的位置识别得不太准确，且会出现双峰现象；归一化标准差法、增强型倾斜角法和改进二阶倾斜角法都可以识别出边界，但是在正常场区域会存在较多的假边界；其他方法受到斜磁化的影响较大，边界识别的误差也较大，不太适合直接用于边界的识别。

3.2 地下两个规模相同但埋深不同的地质异常体模型

为了考察上述不同方法对相邻多个规模相同但埋深不同地质异常体模型的分辨能力，在均匀半空间中设置了两个规模相同(长宽高规模均分别为500 m、500 m和400 m)、剩余密度相同(300 kg/m3)、但埋深不同(左右两个地质异常体的顶面埋深分别为100 m和300 m)的地质异常体模型。磁化强度1 A/m，磁化倾角45°，磁化偏角0°。左边长方体中心位置为坐标X=-400 m，Y=0 m，Z=300 m；右边长方体中心位置为坐标X=400 m，Y=0 m，Z=500 m。

图4是在均匀半空间中放置了两个埋深与规模相同的地质异常体产生的重力异常及使用不同地质异常体识别方法处理的结果图，可以看出，不同的处理方法对于多个地质异常体边界的刻画能力与单个地质异常体情况很相似，对边界位置刻画最准确、分辨率最高的仍是改进二阶倾斜角法和增强型倾斜角法；总水平方向导数法、归一化标准差法的极值均与模型边界基本重合，识别效果较好，但总水平方向导数法受到相邻模型的影响较大，导致相邻两个模型之间的边界突出不够；平面全张量倾斜角法受相邻模型的影响也较大，甚至在二者之间产生了虚假异常；垂向一阶导数法对边界识别能力较好，只是边界略微偏大一点；倾斜角法、Theta图法和广义导数算子法的边界识别结果与模型边界均有偏差，分辨精度不足；解析信号法和解析信号倾斜角法的极值大致能与两个模型位置相对应，并指示出模型的中心位置，但是清晰展现边界的能力稍显不足。

图4 均匀半空间中两个埋深不同地质异常体模型产生的重力异常及使用不同地质异常体识别方法处理后的结果图

图5是两个埋深不同的长方体模型产生磁异常图和不同地质异常体识别方法的处理结果图。同样为了考察倾斜磁化的影响，所有图件也均不做化极处理。由图5可知，垂向一阶导数法的处理结果由于受到倾斜磁化和不同埋深的地质异常体的影响，其零值位置对比复杂，埋深较大的地质异常体反映不明显，不能准确指示地质异常体边界位置(见图5(b))；解析信号法能够准确地识别浅部地质异常体的部分边界位置，但是却无法显示深部地质异常体的边界位置(见图5(c))；总水平导数法同样受到不同埋深的影响，无法有效地识别埋深较大地质异常体的边界(见图5(d))；倾斜角法受埋深的影响较小，能够平衡不同埋深的异常，但是其零值位置与垂向一阶导数的零值位置一样，存在较大误差，对识别边界位置不利(见图5(e))；Theta图法也能够平衡不同埋深的异常，其边界识别结果同样受到倾斜磁化的影响，存在较大误差(见图5(f))；归一化标准差法受地质异常体埋深不同的影响较小，能够平衡不同埋深的异常值，浅部地质异常体边界识别结果较准确，深部的识别结果存在一定误差(见图5(g))；解析信号倾斜角法不易受到埋深的影响，边界识别结果较准确，出现的假边界也较少(见图5(h))；广义导数算子法能够平衡不同埋深的异常，但是其零值位置与地质异常体边界位置相差较大(见图5(i))；增强型倾斜角法的识别结果对比准确，深部地质异常体的边界识别存在一定误差，而且存在假边界(见图5(j))；平面全张量倾斜角法能够平衡不同埋深的异常，埋深较浅的地质异常体边界识别结果较好，深部识别误差较大(见图5(k))；改进二阶倾斜角法的结果与增强型倾斜角类似，但是对边界识别的结果更加清晰(见图5(l))。

图5 均匀半空间中两个埋深不同地质异常体模型产生的磁异常及使用不同地质异常体识别方法处理后的结果图

综上所述，解析信号倾斜角法、归一化标准差法和倾斜角法等使用了比值原理进行归一化的方法受埋深影响较小，能够较好地同时识别出不同埋深的地质异常体边界；没有进行归一化的方法由于受埋深影响，对深部地质异常体边界的识别效果较差。

3.3 添加随机噪音的单一地质异常体模型

实测的重磁异常数据往往包含着各种噪音，在进行识别之前要经过滤波处理以消除噪音干扰。但是滤波等处理通常仅能够消除部分噪音，剩余的噪音往往会影响边界识别的效果。为了测试上述地质异常体识别方法在存在噪音情况下对边界的分辨能力，设置了与3.1完全相同的地质异常体模型，在正演的重磁异常值中添加了5%的随机噪音，以此研究上述方法在未处理干净所有噪音情况下的边界识别效果。

图6是均匀半空间中单个地质异常体产生的重力响应添加5%随机噪音后再滤波的结果及不同地质异常体边界识别方法处理后的结果图，可以看出，垂向一阶导数法、解析信号法和总水平方向导数法受噪音影响最小，非异常区域产生的假异常幅值与异常区域幅值差异明显，能够分辨出地质异常体的边界；倾斜角法、广义导数算子法和平面全张量倾斜角法受噪音影响相对较小，可以分辨出原始地质异常体的位置和形状；Theta图法容易受到噪音的影响，产生了较多虚假异常，虽可大致看出地质异常体边界，但边界不够清晰；归一化标准差法通过选取窗口计算标准差，虽在一定程度上放大了噪音的干扰，产生了较多虚假异常，但对边界识别的效果尚可；解析信号倾斜角法、增强型倾斜角法和改进二阶倾斜角法在计算中均使用了高阶导数，可以提高边界识别的分辨率，但同时也变得对噪音更加敏感，产生了较多的虚假边界，对识别出边界的形状反而不利。

图6 均匀半空间中单个地质异常体产生的重力响应添加5%随机噪音后再滤波的结果及不同地质异常体边界识别方法处理后的结果图

图7是均匀半空间中单个磁异常模型正演数据添加5%噪音并滤波后的结果及不同地质异常体边界识别方法处理后的结果图。同样地均未做化极处理。由图7可知，垂向一阶导数法、解析信号法和总水平导数法受噪音影响较小，处理结果与未添加噪音的结果相差不大，极值位置的幅值与其他区域的幅值相差明显，能够较清晰地识别边界的位置；倾斜角法、Theta图法、归一化标准差法、解析信号倾斜角法、广义导数算子法和平面全张量倾斜角法均对噪音较敏感，边界位置出现了变形，但也可以基本识别出边界；增强型倾斜角法和改进二阶倾斜角法受噪音影响更严重，非边界位置产生了许多虚假干扰，无法对边界进行准确识别。

图7 均匀半空间中单个磁异常模型正演数据添加5%噪音并滤波后的结果及不同地质异常体边界识别方法处理后的结果图

综上所述，在存在噪音的情况下，总水平方向导数法和平面全张量倾斜角法受到的影响较小，能够较好地识别出地质异常体的边界；垂向一阶导数法、倾斜角法和广义导数算子法对噪音不太敏感，抗噪能力较强，虽然其零值位置比边界略大，但也基本能够标识出地质异常体边界；其他方法在存在明显噪音情况下则难以准确识别出地质异常体的边界。

此外，根据上述研究还可以看出，在重磁异常边界识别中，总水平导数法虽然定位边界的位置对比准确，但却不利于识别埋深较大的地质异常体的边界；归一化标准差法和增强型倾斜角法都能够平衡不同埋深的异常，识别出更多的边界，其中增强型倾斜角的识别结果更加清晰。在磁异常边界的识别中，总水平导数法和解析信号法不利于识别埋深较大的地质异常体的边界；归一化标准差法能够识别出更多的边界位置，解析信号倾斜角法的效果最好，不仅能够定位更多的边界位置，而且结果也更加清晰。

4 结束语

本文通过理论模型计算，对重磁勘探中多种地质异常体边界识别方法进行了详细的分析和对比。通过本文的对比研究发现，垂向一阶导数法、倾斜角法和广义导数算子法的极小值位置基本相同，与Theta图法的极大值位置也基本相同。因此上述基于相似数学原理的处理方法对边界的识别效果也相似，在实际应用中只需要选取一种方法即可。其中广义导数算子法相较于其他方法，既能够平衡不同深度的异常，而且受噪音的影响也更小，相对来说处理效果更好；解析信号倾斜角法、增强型倾斜角法和改进二阶倾斜角法使用了高阶导数的方法，虽然分辨率更高，但同时对噪音也更加敏感，使用之前应该首先对原始数据进行充分的滤波等去噪处理；总水平方向导数法的边界识别结果较准确，受噪音影响也较小，但其也存在无法平衡不同埋深异常的缺点。通过本文研究，总体上认为，总水平方向导数法、增强型倾斜角法和归一化标准差法相较于其他地质异常体边界识别方法效果更好。