位场垂向梯度最佳自比值的边界检测技术

王彦国，张凤旭，刘 财，王祝文，张 瑾

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

2东华理工大学核工程与地球物理学院，南昌 330013

1 引 言

应用位场数据识别地质体边界是地质解释的一项重要工作.位场异常中包含有场源边界的信息，但边界信息的提取需要进行相关的数据处理.为此，地球物理工作者提出了许多基于位场梯度的边界检测方法.这些方法按边界检测的计算模式，可分为导数分析和数理统计两类［1］.

导数分析法主要有垂向导数分析法［2］、水平总梯 度［3］、tilt angle［4］及 其 水 平 总 梯 度［5］、Theta map［6］等方法.垂向导数分析法是Evjen提出的，主要是利用垂向一、二阶导数的零值线确定地质体边界［2］；Cordell提出了水平梯度法，利用水平梯度模极大值检测场源边界［3］；这两种方法有效地提高了对原异常信息的识别能力，但难以同时反映出深度不同的地质体边界［7］.Miller等提出了tilt angle法，该法识别地质体边界的原理与垂向导数类似［4］；Verduzco在tilt angle法基础上提出了tilt angle水平总梯度法，用计算结果的极大值检测地质体边界，但在计算过程中采用了异常的二阶导数，更易受高频干扰的影响［5］.Wijns等采用水平总梯度与一阶解析信号比值（Theta map）来分析磁源边界，对低纬度磁性体有较好的边界检测能力，但检测的边界与深部地质体的实际边界位置有偏差［6］.

数理统计法主要有Euler反褶积［8－10］、小子域滤波［11］、导数归一化标准差［12］、均方差比归一化垂向梯度［13］等方法.Euler反褶积法是将位场异常和其梯度联系在一起的Euler齐次方程，该方法可用于地质体边界的识别，但方程式灵敏度较高，导数计算结果稍有偏差，便会在计算结果中产生较大的偏离，从而影响反演的精度［14］.小子域滤波法是杨高印针对异常分离提出的，可用于边界的识别，但只对原异常中以梯级带形式显示的边界位置有效［11］；张凤旭等［15－16］在小子域滤波基础上，提出了三方向小子域滤波法，获得了较为详细的断裂构造信息，但计算结果与实际地质体边界也有一定的偏差.Cooper等［12］借助于异常导数的标准差，提出了可以检测不同深度地质体边界的归一化标准差法；王彦国等［13］在2012年提出了均方差比归一化垂向梯度的边界检测方法，对不同深度的地质体边界也有较好的探测效果，但这两种方法在计算的过程中同样采用了导数计算，因此也会受噪声等高频干扰.

显然，上述边界检测方法大都是以位场梯度计算为基础的，方法运算简单，物理意义也明确，但在实际资料处理中都存在着受高频干扰的影响［14］.在处理这个问题上，如果直接采用低通滤波进行降噪，则会损失部分有用信号.

在利用位场异常进行场源边界识别中，一阶导数分析是目前常用的一类方法.不过笔者认为，导数的阶次越高，对场源的边界识别能力和定位精度就越强，然而高阶导数易受噪声等干扰的影响，而且导数的阶次越高，高频干扰的影响就越为严重.如果能采用某种数学处理手段，既能尽量减少导数异常中有效边界信息的损失，又具有抗干扰能力的话，势必可以提高场源的边界识别能力.为此，本文提出了位场垂向梯度最佳自比值的边界检测方法.

2 方法原理分析

2.1 自比值的定义

对于位场异常数据f（i，j），设其垂向m阶导数为fz（m），那么可以将导数异常fz（m）窗口内数据的均值与其均方根的比值定义为fz（m）的一次自比值κ，即

2.2 自比值数值分析

构造如下二次函数Fij（x）：

2.3 位场垂向梯度自比值的物理意义

设垂向导数fz（m）在所选定的计算窗口内数据平均值为：

2.4 最佳自比次数和最佳窗口长度的确定

依据文献［17－19］，可定义垂向m阶导数的第n次自比值与第n＋1次自比值的归一化互相关系数R为：

其中M、N分别为测线数和每条测线的测点数.

事实上，在选定的窗口下，计算的自比值互相关系数与自比次数的关系曲线中存在一个明显极大值（后文将进一步阐述），现将这个极大值点对应的自比次数定义为最佳次数，将该次数对应的自比值称为最佳自比值.

由于在选定的初始窗口下，最佳自比值未必是确定地质体边界位置的最佳结果，因此需要选择最佳窗口的长度.最佳窗口长度的选取方案是：首先设定一个置信度R0（其取值一般介于0.96和0.99之间），如果在初始窗口下，自比值互相关系数的极大值max（R）大于给定的置信度R0，则认为最佳自比值去噪能力强，结果是可信的，那么初始窗口即为最佳窗口；若max（R）不大于R0，则认为初始窗口下的自比值结果可信度较低，需调整计算窗口的长度，重新计算自比值互相关系数与自比次数的关系，直到某一窗口下满足max（R）＞R0时，终止计算，这时所对应的窗口大小便是最佳计算窗口长度.

2.5 最佳自比值计算步骤

由以上分析知，本文提出的垂向梯度最佳自比值法在给定置信度R0之后，可以自适应地寻找到最佳窗口下的最佳自比次数，然后通过最佳自比值进行场源边界的识别.为了方便起见，在此给出算法的步骤：

（1）采用快速傅里叶变换计算位场数据的垂向梯度fz（m）；

（2）给定初始窗口长度（一般选取D＝5），利用公式（1）和（2）计算fz（m）的n次自比值

（4）给定置信度R0，将初始窗口下的互相关系数极大值max（R）与R0进行比较.若max（R）＞R0，则输出最佳自比值；若max（R）≤R0，则调整窗口长度，重新计算（2）、（3）步骤，直到满足max（R）＞R0为止.

值得说明的是：为了进一步提高地质体边界的定位精度，可采用杨高印提出的小子域滤波［11］对最佳自比值进行滤波输出.

3 理论模型分析

为了验证本文方法的边界识别能力，设计了一个包含4个不同参数的长方体组合模型（模型体参数见表1，位置见图1a），其产生的理论重力异常见图1b.同时为了验证方法消除噪声的能力，在叠加模型体产生的重力异常中加入了随机噪声（图1c）.

从模型重力异常（图1b）中可以看出，地质体A由于埋深和规模较大，其产生的异常在边界位置表现出的梯级带范围较宽，因此利用该异常所显示的边界特征是无法精确确定地质体边界位置的；地质体B、C、D规模较小，受地质体A的影响，其产生的异常在叠加异常中表现为等值线同形扭曲，显然只依靠这种异常特征难以确定出这三个地质体的边界；而且在加入噪声后，异常图（图1c）在三个地质体位置所显示的等值线同形扭曲特征变的模糊不清，这进一步增加了边界识别的难度.

表1 模型参数Table 1 The model parameters

在模型试验中，笔者先用tilt angle、水平总梯度、Theta map、导数归一化标准差、小子域滤波和均方差比归一化垂向梯度6种已有的边界识别方法对含噪声的重力异常分别进行了计算.各种方法计算的结果见图2，从中可以看出，由于受噪声和异常叠加的影响，6种方法的计算结果均无法清晰地显示模型体A的边界，更无法有效地识别其它3个较试验表明，该极值点就是确定最佳窗口大小和自比次数的重要依据.

在计算最佳窗口和最佳自比次数（图4）时，首先选择初始窗口D＝5，由于该窗口条件下垂向一阶和二阶导数的自比值互相关系数极大值均大于可信度R0（文中选取R0＝0.98），因此可确定垂向一阶、二阶导数的最佳窗口与初始窗口相同（D＝5），相对应的最佳自比次数分别为n＝3和n＝4.对于垂向三阶导数，由于在初始窗口下的互相关系数极大值max（Rz（3））＝0.956，小于R0，因此将计算窗口长度增加至D＝7，此窗口下的max（Rz（3））＝0.994，大于R0，即垂向三阶导数自比值的最佳窗口为D＝7，其对应的最佳自比次数为n＝4.

图3 含噪声重力异常垂向导数（a）垂向一阶导数；（b）垂向二阶导数；（c）垂向三阶导数.Fig.3 Vertical derivatives of gravity anomaly including random noise（a）First－order vertical derivative；（b）Second－order vertical derivative；（c）Third－order vertical derivative.

图4 自比值互相关系数与自比次数的关系曲线Fig.4 The relationship of cross－correlation coefficient of two successive auto－ratios and the number of auto－ratio

以上对比分析充分证明，本文提出的边界检测方法不但具有较高的精度，而且计算过程更为稳定.

4 实际重力资料1）应用

为了验证垂向梯度最佳自比值对实际资料的处理效果，选取了吉林省南部鸭绿江盆地2）实测重力数据进行试验.

鸭绿江盆地位于中朝板块东北缘，二级大地构造单元隶属于辽东台隆区，主体位于太子河—浑江坳陷内（图6）.在盆地内，除志留纪、泥盆纪和早石提供了较好的对比依据.

从研究区布格重力异常（图8a）中可以看出，反映地质体边界位置或断裂构造的梯级带主要分布在二道江—新立屯—板石镇一线、四道江—六道江—白山—孙家堡子一线、孤砬子—红土崖—三道湖—石人一线、公益村—大路村一线以及蚂蚁河—闹枝沟屯圈闭区等位置，这几组异常梯级带可能是大型地质体边界（或规模较大的断裂构造）的反映.但这些异常梯度带较平缓，边界位置难以直接从重力异常图中精确定位.另外，较小型地质体边界受区域异常影响较大，在异常图中主要表现为异常等值线突然变宽或变窄以及同形扭曲等非梯级带特征，这类地质体边界位置确定难度更大.

图8b为结合区域地质和前人地质－地球物理工作成果勾划的构造分区图3）.图中给出的大型断裂位置、燕山期花岗岩出露区、中新生代凹陷分布范围以及浑江煤田工作区等信息，可为文中方法的有效性提供佐证.

图8c为研究区垂向三阶导数自比值互相关系数与自比次数的关系曲线，仿照图4做法，可检测出最佳窗口长度D＝5和最佳自比次数n＝4，进而获得了垂向三阶导数最佳自比值的计算结果（图8d）.从中可以看出，最佳自比值不仅清晰地反映出龙岗隆起和浑江坳陷之间以及浑江坳陷和老岭推覆区之间的大型断裂，还较好地显示出了不同岩性间的接触带（图7）.另外，自比值圈定的负异常较好地反映出了研究区中－新生代地层和燕山期花岗岩等相对低密度岩石的分布，也客观地反映出了浑江煤田的工作范围.这些结果与已知地质资料以及前人的工作成果符合较好，有力地佐证了方法的有效性.

图7 鸭绿江盆地地质图Fig.7 Geology map of Yalujiang basin

5 结 论

有效地利用位场相应阶次的导数异常，能够提高地质构造解释的分辨率，但高阶导数换算对干扰的放大作用一直困扰着地球物理工作者.本文在前人工作的基础上，另辟蹊径，提出了位场垂向梯度最佳自比值的边界检测方法.在文中，笔者对方法进行了数值分析，并阐述了该算法检测地质体边界的物理机制.模型试验和实际算例表明，垂向梯度最佳自比值算法不但能够有效地消除导数中的干扰成分，而且能够精细地识别出地质体的边界.

本文实现了高阶导数在位场边界识别中的应用.由于自比值可以有效地压制干扰，因此可以提供更丰富、精确的边界信息，为高阶导数在位场数据处理中的有效应用提供了新的思路.

图8 （a）鸭绿江盆地布格重力异常（单位10－5 m·s－2）；（b）研究区构造分区图；（c）垂向三阶导数自比值互相关系数与自比值次数关系曲线；（d）垂向三阶导数最佳自比值Fig.8 （a）Bouguer gravity anomaly of Yalujiang basin（unit：10－5 m·s－2）；（b）The tectonic division of research area；（c）The relationship of cross－correlation coefficient of two successive auto－ratios of third－order vertical derivative and the number of auto－ratio；（d）Optimal auto－ratio of third－order vertical derivative

［1］ 王万银，邱之云，杨永等.位场边缘识别方法研究进展.地球物理学进展，2010，25（1）：196－210.Wang W Y，Qiu Z Y，Yang Y，et al.Some advances in the edge recognition of the potential field.ProgressinGeophys.（in Chinese），2010，25（1）：196－210.

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

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

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

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

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

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

［8］ Thompson D T.EULDPH：A new technique for making computer－assisted depth estimates from magnetic data.Geophysics，1982，47（1）：31－37.

［9］ Reid A B，Allsop J M，Granser H，et al.Magnetic interpretation in three dimensions using Euler deconvolution.Geophysics，1990，55（1）：80－91.

［10］ Hansen R O，Laura S.Multiple－source Euler deconvolution.Geophysics，2002，67（2）：525－535.

［11］ 杨高印.位场数据处理的一项新技术——小子域滤波法.石油地球物理勘探，1995，30（2）：240－244.

Yang G Y.A new technique for potential－field data processing：small subdomain filtering.OilGeophysical Prospecting（in Chinese），1995，30（2）：240－244.

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

［13］ 王彦国，王祝文，张凤旭等.基于均方差比归一化垂向梯度法的位场边界检测.中国石油大学学报（自然科学版），2012，36（2）：86－90.

Wang Y G，Wang Z W，Zhang F X，et al.Edge detection of potential field based on normalized vertical gradient of mean square error ratio.JournalofChinaUniversityof Petroleum（EditionofNaturalScience），2012，36（2）：86－90.

［14］ 张恒磊，刘天佑，杨宇山.各向异性标准化方差计算重磁源边界.地球物理学报，2011，54（7）：1921－1927.

Zhang H L，Liu T Y，Yang Y S.Calculation of gravity and magnetic source boundary based on anisotropy normalized variance.ChineseJ.Geophys.（in Chinese），2011，54（7）：1921－1927.

［15］ 张凤旭，张凤琴，刘财等.断裂构造精细解释技术——三方向小子域滤波.地球物理学报，2007，50（5）：1543－1550.

Zhang F X，Zhang F Q，Liu C，et al.A technique for elaborate explanation of faulted structures：three－directional small subdomain filtering.ChineseJ.Geophys.（in Chinese），2007，50（5）：1543－1550.

［16］ 张凤旭，张兴洲，张凤琴等.中国东北地区重力场研究——利用改进的三方向小子域滤波划分主构造线及大地构造单元.地球物理学报，2010，53（6）：1475－1485.

Zhang F X，Zhang X Z，Zhang F Q，et al.Study of gravity field in Northeastern China area：Classification of main structure lines and tectonic units using the improved threedirectional small subdomain filtering.ChineseJ.Geophys.（in Chinese），2010，53（6）：1475－1485.

［17］ 高光珠，李忠武，余理富等.归一化互相关系数在图像序列目标检测中的应用.计算机工程与科学，2005，27（3）：38－40.

Gao G Z，Li Z W，Yu L F，et al.Application of the normalized cross correlation coefficient in image sequence object detection.ComputerEngineeringandScience（in Chinese），2005，27（3）：38－40.

［18］ 曾华霖，许德树.最佳向上延拓高度的估计.地学前缘，2002，9（2）：499－504.

Zeng H L，Xu D S.Estimation of optimum upward continuation height.EarthScienceFrontiers（in Chinese），2002，9（2）：499－504.

［19］ 孟小红，刘国锋，陈召曦等.基于剩余异常相关成像的重磁物性反演方法.地球物理学报，2012，55（1）：304－309.

Meng X H，Liu G F，Chen Z X，et al.3－D gravity and magnetic inversion for physical properties based on residual anomaly correlation.ChineseJ.Geophys.（in Chinese），2012，55（1）：304－309.