重力数据LTHG 均衡边界识别方法对比分析及其应用

彭莉红，张伟盟，程莎莎，孙栋华，骆燕，陈伟

1 核工业航测遥感中心，河北 石家庄 050002

2 中核集团铀资源地球物理勘查技术（重点实验室），河北 石家庄 050002

3 河北省航空探测与遥感技术重点实验室，河北 石家庄 050002

边界识别方法在重力数据解释中起着至关重要的作用，国内外学者对它进行了广泛研究，提出多种边界识别方法［1-16］。众多边界识别方法以导数分析类为主，而导数分析类中最基本的方法是总水平梯度（THG）［1-3］、垂向导数（VDR）［4-6］和解析信号（AS）［7-9］，但上述三种方法均存在异常幅度随地质体埋深的增加而迅速衰减的情况，因此对于埋深较深、异常幅度较小的地质体，边界识别能力稍显不足。近年来，多种通过求取重力数据导数比值来均衡不同幅度异常的边界识别方法相继被提了出来。Miller 和Singh［17］提出了均衡不同幅度异常的边界识别方法——倾斜角法（TDR），该方法对于场源体的深度不敏感，对埋深不同的多场源异常体的边界识别效果较好；Verduzco 等［18］针对TDR 的处理结果峰值位于场源体中心，识别结果不直接的问题，提出了倾斜角的总水平导数法（THDR），该方法采用极大值识别场源边界，结果相较于TDR 简单、直接，且具有更高的分辨率，然而也存在随场源深度的增加、异常幅度降低的情况［19］；Wijns 等［20］提出另一种均衡不同振幅的识别方法，利用解析信号与总水平导数的比值来进行地质体的边界识别，即Theta图法（TM），该方法相对于解析信号法的分辨能力增强，缺点是所识别出的边界较模糊。Ferreira 等人［19］提出了总水平梯度的倾斜角法（TTHG），该方法能较好地识别出不同埋深的场源体边界信息；Cooper［21］提出了另一种改进的倾斜角法，称为解析信号的倾斜角法（TAS），虽然TAS 可以更好地识别出埋深较深的场源体边界信息，但存在与AS 法同样的不足，即随埋深的增加，其极大值位置会向实际边界的内部收敛。王万银等［22］提出了归一化总水平导数垂向导数边界识别方法（NVDR_THDR），该方法采用零阈值技术消除地质体边界，结果直观，易于识别，是均衡不同埋深场源边界识别效果的有效方法。马国庆等提出增强型均衡滤波器，该滤波器利用不同阶导数之间的组合识别地质体边界。该方法能有效识别深、浅部地质体边界，但高阶导数的引入使噪声不可避免地被放大，从而影响整体边界识别的稳定性［23-24］。Luan 等提出基于Logistic 函数及总水平导数的重力数据边界增强方法（LTHG 法），认为该方法处理结果简单、直接且对于场源边界具有更高的分辨能力。LTHG 方法是目前最新提出的可用于重力数据的边界识别方法，但还未有其与主流边界识别方法效果对比的相关研究成果，也缺少实际应用方面的公开论述。

为此，本文基于理论模型及实测数据，对比分析了不同情况下LTHG 法与常用的THG、AS、TDR、THDR、TTHG、TAS、TM 和NVDR_THDR8 种方法对不同埋深场源体的均衡能力、分辨能力及抗噪能力，进而验证LTHG 法的适用性及应用效果。

1 LTHG 方法简介

LTHG 法是Luan 等（2019）新提出的一种基于Logistic 函数及总水平导数的位场数据均衡边界识别方法，其核心原理是基于Logistic 函数，通过对总水平导数的一阶垂向导数和水平导数的比值进行Logistic 函数运算来达到突出显示不同深度的场源体边界的效果。

Logistic 函数是机器学习领域中常用的一种激活函数，定义域（-∞，+∞），值域为（0，1），形式如下：

Logistic 函数（函数曲线形式见图1）在定义域内单调连续，呈现出先缓慢增长，然后加速增长，最后逐渐稳定的趋势，与常用于识别重力数据边界的arctan 函数相似。

图1 Logistic 函数变化曲线Fig. 1 The Logistic function curve

基于Logistic 函数，LTHG 法的计算公式如下：

式（2）中： THG—重力数据的总水平导数；a—常数。LTHG 法求解后得到一个比值，是无量纲，其平面峰值表征了地质体的边界。

图2 为二维情况下，采用LTHG 法a分别取值为0.5、1.0、5.0、10、20、50 和100 时直立长方体密度模型的边界识别结果，可以看出a 取不同参数时，均可以较好的识别出模型的边界位置，且边界位置处异常幅值基本相同，但不同取值计算的结果还是存在一定的差异，当a＞10 时，计算结果当中，边界位置处异常峰型会变的更加尖锐，但当a取值大于10 时，处理结果中会因存在畸变点而出现空值的情况，当a＜1 时，模型边界位置处异常峰型平缓，1≤a＜10 时LTHG 法应用效果较好。

图2 不同参数LTHG 法边界识别结果Fig. 2 The LTHG filter of gravity anomaly due to the 2D block with different a.

2 边界识别方法对比

为说明LTHG 法识别边界的有效性，选取八种常用的重力数据的边界识别方法进行对比，即总水平梯度（THG）、解析信号（AS）、倾斜角法（TDR）、倾斜角的总水平导数法（THDR）、Theta 图法（TM）、总水平梯度倾斜角法（TTHG）、解析信号倾斜角法（TAS）和归一化总水平导数垂向导数（NVDR_THDR）法，各方法定义如下：

THG 法［25］：

AS 法［25］：

TDR 法［25］：

THDR 法［25］：

TM 法［25］：

TAS 法［25］：

TTHG 法［25］：

NVDR_THDR 法［22］：

计算总水平导数（THG）的n阶垂向导数（VDRn），n—垂向导数阶数，n=1，2，3…。阶数越大，横向分辨力越高，通常对于重磁异常，阶数取2 时较为合适；

使用阈值大于0 来计算总水平导数峰值（PTHDR）：

计算总水平导数峰值及总水平导数的比值：

计算总水平导数垂向导数的最大值（VDR - THDRmax），并使用最大值进行总水平导数垂向导数归一化，最终得到归一化总水平导数垂向导数（NVDR-THDR）：

建立了由5 个不同规模、不同埋深的长方体组成的组合模型，模型三维立体图见图3，参数见表1，正演计算了3 种情况下组合模型的理论重力异常。情况1 为所有长方体的剩余密度均为正值，长方体物体1 和2 剩余密度为0.5×10-3kg·m-3，3、4、5 体的剩余密度为0.4×10-3kg·m-3，计算的重力异常如图4A 所示。情况2 为长方体1、4 剩余密度为负值，体1 剩余密度为-0.5×10-3kg·m-3，体4 剩余密度为-0.4×10-3kg·m-3，体2、3、5 剩余密度与第一种情况保存一致，计算的理论重力异常如图5A 所示。情况3 为在情况2 数据中加入了2 %的高斯噪声，结果如图6A 所示。选用包括THG、AS、TDR、THDR、TM、TAS、TTHG、NVDR_THDR 及LTHG在内的九种方法对以上三种情况下的理论重力异常数据进行处理，对比分析各方法的不同振幅重力异常的均衡能力、边界分辨能力及抗干扰能力。

表1 组合模型参数Table 1 Parameters of the synthetic model

图3 模型三维立体图Fig.3 A 3D view of the synthetic model

图4 模型理论重力异常（重力正异常）及边界识别结果Fig. 4 Synthetic gravity anomaly of the first model and the edges detected by different methods

图5 模型理论重力异常（重力正、负异常）及边界识别结果Fig. 5 Synthetic gravity anomaly of the second model and the edges detected by different methods

图6 加入噪声的模型理论重力异常（重力正、负异常）及边界识别结果Fig. 6 Synthetic gravity anomaly of the third model and the edges detected by different methods

图4B、C 为THG、AS 的处理结果，可以看出当存在多个埋深不同的场源体时，THG 和AS不能均衡不同埋深的场源体边界信息，随着深度增加，场源体边界变得越来越不清晰，同时虽然薄板1 和2 埋深较浅，但AS 仍无法很好的显示出它们的边界信息。图4D～I 分别为TDR、THDR、TM、TAS、TTHG 和NVDR_THDR的处理结果，虽然TDR 均衡了不同埋深的场源体的边界信息，但处理结果的峰值位于场源体中心，处理结果不直接。THDR 法克服了TDR法的不足，利用极大值来识别场源体的边界，但对于埋深较深的场源体，边界分辨能力仍稍显不足。TM 法识别的边界位置不够准确，且存在假的边界信息。TAS 法同样利用极大值来识别场源体边界，然而对于较深的场源，其峰值会从实际的边界位置向内偏移，使识别到的边界位于真实边界内；对于较薄场源，TAS 的峰值直接位于场源体中心。TTHG 法较好地均衡了不同埋深的场源体边界信息，识别出的边界信息接近真实的场源体边界，但随着场源埋深的增加，识别的边界信息也存在一定程度的扩散。NVDR_THDR 法可直接用峰值识别出边界，效果比上述各方法更好。图4J 为 LTHG 法的处理结果，可以看出LTHG 很好地克服了导数计算随深度增加而迅速衰减的问题，均衡了不同埋深场源体的边界信息，结果与NVDR_THDR 相似，但LTHG 方法对边界的分辨能力要更强。对比分析结果表明：LTHG 不仅能更清晰、准确地识别出场源边界，而且能给出较高的边界分辨率。

为对比分析上述各方法同时识别正负异常边界的能力，选用包括THG、AS、TDR、THDR、TM、TAS、TTHG、NVDR_THDR 及LTHG在内的9 种方法对第2 种情况下的理论重力异常数据进行了处理。图5B、C 为THG、AS 的处理结果，同样THG 法和AS 法不能同时均衡振幅不同的异常。图5D～I 分别为TDR、THDR、TM、TAS、TTHG 和NVDR_THDR 的处理结果，由图5D 可知，当同时存在正负异常时，TDR 法无法准确识别场源的边界，且结果中出现了假边界信息。THDR 法也存在假边界信息。TM法除不能很好地识别出埋深较深或薄的场源体边界外假边界信息更加显著。TAS 法能避免假边界信息的产生，但识别到的边界位于真实边界内。TTHG 能避免假边界信息的产生，但识别的边界信息随着场源埋深的增加有一定程度的扩散。NVDR_THDR 法识别出的边界信息轮廓清晰、位置准确，无虚假边界信息。图5J 为 LTHG 法的处理结果，可以看出LTHG 法识别的边界位置接近真实场源边界，无假边界信息，相较于NVDR_THDR，LTHG 方法对边界的分辨能力要更强。

为对比分析各方法的抗噪能力，在第二种情况的数据中加入幅度为异常振幅2 %的高斯噪声。图6B～J 分别为THG、AS、TDR、THDR、TM、TAS、TTHG、NVDR_THDR 及LTHG 的处理结果。THG、AS 法对浅源体3 和4 的边界识别效果相对较好，但对深源体5 的边界识别效果却相对模糊，同时对于薄板1 和2 ，AS 法也不能很好地发挥作用。TDR 和TM 依然无法对场源体边界进行准确的识别，受噪声的影响THDR、TAS 法识别的场源边界变的模糊，分辨率降低，TTHG、NVDR_THDR 和LTHG 三种方法依然可以较准确的识别出场源体边界位置，且NVDR_THDR 和LTHG 受深度变化的影响小，未出现明显的边界向内或外侧偏移的情况，相对而言LTHG 法识别的位置更清晰，分辨率更高。

3 应用实例

埃塞俄比亚欧加登盆地是非洲东海岸油气勘探前景最好的盆地之一，是发育在古生代结晶基底之上的一个古内陆裂谷—被动大陆边缘盆地，根据基底构造特征，可进一步划分为6 个次一级构造单元：北部隆起、北部斜坡、西部斜坡、中部坳陷、南部斜坡和南部隆起（图7），目前在中部坳陷带发现了C、H 和 D 三个气田，展示了盆地良好的勘探潜力［26-30］。但盆地的勘探还处于早期评价阶段，对断裂的整体空间延伸及展布规律还认识不清。对重力数据及其转换参数进行分析，可以为研究盆地断裂构造发育特征及其与油气成藏的关系提供重要依据［31-33］。

图7 欧加登盆地构造单元划分及研究区位置（据参考文献［30］）Fig. 7 Location and simplified geological maps of the study area

核工业航测遥感中心和中国石油集团东方地球物理有限责任公司于2018 年在欧加登盆地中部共同开展了1：10 万航空重力测量，获得一套高精度航空重力数据。本文为验证LTHG 法实际应用效果，选取欧加登盆地中部坳陷带作为研究区（图7），对研究区的航空布格重力数据进行LTHG 处理，并与THG、AS、TDR、THDR、TM、TAS、TTHG 和NVDR_THDR在内的几种方法进行对比分析，为降低数据噪声的影响，对布格重力异常数据进行了向上延拓，延拓高度为2 km，经延拓处理后识别的场源体边界会更加光滑，但场源体的形态不会发生改变；LTHG 法a取值为5（a为5 时，识别结果足够尖锐且无空值现象）。从各方法的识别结果可以看出，研究区断裂走向整体呈NE 和NEE 向，与地质资料相吻合，但各方法的识别结果及识别精度上存在一定差异。由图8A、B可见，THG 法和AS 法对幅值较大的异常识别效果较好，对幅值较低的识别效果稍显不足；图8C 显示，TDR 虽然均衡了不同埋深的场源体的边界信息，但无法产生锐利的边界信息；图8D、E 和F 显示THDR、TM 和TAS 识别效果不甚理想，假边界信息较多。图8G、H 和I 显示TTHG、NVDR_THDR 和LTHG 三种方法的识别效果较好，处理得到的重力异常线性特征更加明显，同时又显示更多的地质细节，但相对而言LTHG 识别的结果更加直观，分辨率更高。

图8 研究区航空布格重力异常边界识别结果Fig.8 The edge recognition results of airborne gravity data in the study area

结合研究区的区域地质资料和其他物探资料，以LTHG 方法为主，TTHG、NVDR_THDR等方法为辅助，推断了研究区内的断裂共10 条（表2），以近NE 向和近NW 向为主，其中，NE 向断裂5 条，分别为 F1、F2、F3、F4、F5，NW 向断裂5 条，为 F6、F7、F8、F9、F10。 将本次的边界识别结果与前人文献［30］进行对比分析，如图9 所示。本次推断的NE 向断裂与前人研究成果基本一致，同时新推断了一组NW 向断裂，从布格重力数据及边界识别结果来看，NW 向断裂基本错切了NE 向断裂；前人研究成果中，中部凹陷、中央凸起及西部凹陷北界断裂与本次推断的F2 断裂一致，南界断裂与本次推断F3 断裂一致。

表2 研究区断裂推断成果表Table 2 The results of inference fault in the study area

图9 研究区推断断裂与前人解释成果对比图Fig.9 Comparison of interpretation faults and previous interpretation results in the study area

4 结 论

1）理论模型的合成重力数据处理结果显示，LTHG 法能很好的识别出不同埋深的场源体边界，提供更丰富的场源体信息。与其他常用方法进行比，LTHG 法的识别的结果准确、清晰，分辨率更高，且具有一定的抗干扰能力；

2）理论模型试验表明，采用LTHG 法进行边界识别时，边界位置不受a取值的影响，但当a＞10 时，识别结果会出现数据缺失的情况，a＜1 时，模型边界位置处异常峰型平缓，因此1≤a＜10 时，LTHG 法应用效果较好；

3）对实测航空重力数据处理得到的重力异常总体上与研究区的断裂有较好的相关性，同时显示了更多的地质细节，最终以LTHG 法处理的结果为主，并综合其他方法对研究区的断裂进行了推断解释，推断NE 向断裂与前人研究成果基本一致，新推断一组NW 向断裂，且NW 向断裂错切了NE 向断裂。