一种改进的重力梯度全张量数据的边界识别方法

2021-09-06杜威吴贺宇陈祥忠刘俊成张阳阳惠梦琳

地球物理学报 2021年9期

杜威，吴贺宇，陈祥忠，刘俊成，张阳阳，惠梦琳

1 中国科学院地球化学研究所矿床地球化学重点研究实验室，贵阳 550081 2 吉林大学地球探测科学与技术学院，长春 130026 3 北京桔灯地球物理勘探股份有限公司，北京 102200 4 安徽省勘查技术院，合肥 230031

0 引言

随着重力梯度测量技术的发展，越来越多的重力梯度全张量数据被应用到数据处理中.重力梯度全张量数据是重力位的二阶偏导数，由9个张量数据组成，较重力数据包含有更高频率的信息.对重力梯度全张量数据进行处理，可以得到更全面，更准确的地质信息.

边界识别在位场数据处理中有着重要的作用(Cooper and Cowan，2006；马国庆等，2011；Archibald et al.，1999；鲁宝亮等，2020)，其方法种类有很多，常用的有：垂向导数(Evjen，1936)、总水平导数法(THDR)(Verduzco et al.，2004)、解析信号(ASM)(Hsu et al.，1996)、倾斜角法(Miller and Singh，1994)、倾斜角水平导数法(王明等，2013)、Theta图法(Wijns et al.，2005)以及其他方法(段杰翔和徐亚，2020；王丁丁等，2021)，这些方法都是在传统的重力数据基础上提出的.随着重磁勘探方法的多样化和高精度勘探仪器的发展，梯度仪可以测量所有的重磁张量分量.许多经典的边缘检测方法被应用于梯度张量数据(郑强等，2019).Zhang等(2000)提出张量-欧拉反褶积方法，提高了地质体水平位置的分辨率(Mikhailov et al.，2007).Beiki(2010a，b)利用重力张量数据的方向解析信号振幅来确定地质体位置.朱自强等(2015)通过对重力梯度张量解析信号应用欧拉反褶积，提高了反演结果的可靠性.Yuan和Yu(2015)定义了二阶方向解析信号，提出了几种基于水平方向解析信号和二阶水平方向解析信号的新边界识别方法.Oruç和Keskinsezer(2008)使用倾斜角法，利用重力梯度张量数据对简单模型的边缘进行检测.另外，重力梯度张量的特征值也被用来描绘地质体的边缘.Beiki和Pedersen(2010)利用重力梯度张量数据的特征向量来估计场源的位置和走向.Sertcelik和Kafadar(2012)提出了一种利用结构张量特征值来识别地质体边缘的方法.Yuan等(2014a)提出了三种平衡结构张量最大特征值的方法，并重新定义了结构张量.Zhou等(2017)提出了基于三维结构张量方向总水平导数的边缘检测新方法.Oruç等(2013)利用重力梯度张量数据曲率的最大特征值圈定了地质构造的边缘.一些学者基于Oruç 等提出的方法进行了相关改进(Zhou et al.，2013；Wang et al.，2015；Zuo and Hu，2015).Yuan等(2016)提出了方向Theta方法，并利用水平方向Theta组合(ThetaX与ThetaY相加)的最大值定义了一种新的边缘检测方法(ED)来描绘重力梯度张量数据的边界.当地质体同时存在正异常和负异常时，该方法可以避免产生虚假边缘，但该方法受地质体走向的影响.本文针对其受地质体走向影响的问题，提出了一种改进的边界识别方法(IED)，替换原方法中的加法运算，通过选择合适的阈值来削弱虚假边界的识别.模型试验表明，该方法不再受地质体走向的影响，在实际数据处理中取得了良好的效果.

1 ED方法及存在问题

重力梯度全张量(GGT)数据是重力位U的二阶偏导数，其矩阵形式为

(1)

由位场理论可知，地质体外的引力位满足拉普拉斯方程(Beiki et al., 2010)：

gxx+gyy+gzz=0，

(2)

T为对称矩阵，它只包括5个独立张量，gxy、gxz、gyz、gxx和gzz.方向解析信号的振幅为

(3)

Yuan和Geng(2014b)、Yuan和Yu(2015)定义了方向总水平导数.基于Theta方法的定义，提出了GGT数据的方向Theta方法(Yuan et al.，2016)，可以表示为

(4)

并在此基础上，将ThetaX与ThetaY得到的结果进行叠加，提出了ED方法，该方法用极大值识别边界：

ED=ThetaX+ThetaY.

(5)

本文设计了模型试验(图1).图1a、b分别是模型的平面图和三维视图.测区范围为100 km×100 km，点距为0.5 km，模型参数见表1.

图1 模型位置图 (a) 平面图； (b) 立体图.Fig.1 Positions of models (a) Planar view; (b) 3D view.

表1 图1模型的参数Table 1 Parameters of the models in Fig.1

在模型试验中，当重力梯度张量用于边界识别时，可能会受到地质体走向的影响，因此模型体设计成与坐标轴成一定角度.此外，模型的设计考虑了不同深度地质体及正负异常相间对边界识别的影响.图2和图3分别为模型的重力异常图和理论重力张量图.

图2 模型重力异常图Fig.2 Gravity anomalies of models

图3 模型理论张量 (a)gxx; (b) gyy; (c) gzz; (d) gxz; (e) gyz; (f) gxy.Fig.3 Synthetic gravity gradient tensors of models

本文将总水平导数(THDR)方法的边界识别效果(图4)与ED方法的边界识别效果(图5)进行对比，其中THDR方法为重磁数据处理中常用的边界识别方法，该方法受深度影响.图4中，模型2较模型1、3深度大，所以模型2的边界识别效果较模型1、3差，边界模糊.

图4 THDR方法边界识别结果Fig.4 Edge detection resultby THDR

图5是ED方法边界识别效果，Yuan等(2016)等已经验证该方法不受正负异常梯度带的影响.与THDR方法识别结果进行比较，其中模型2的边界识别效果与模型1的相同，并且比THDR方法识别出的边界清晰，说明该方法不受深度的影响.图5中模型3与模型1的深度相同，但是识别出的边界清晰程度不同，考虑到只有模型3边界垂直于x轴或y轴，所以推测该方法受地质体走向的影响.

图5 ED方法边界识别结果Fig.5 Edge detection result by ED

从ThetaX和ThetaY的等值线图(图6)中可以看出，模型体1，2的边界在ThetaX和ThetaY中分别对应极大值，且大小相近.而模型体3，在ThetaX的图中只能识别出垂直于x轴的边界，在ThetaY的图中只能识别出垂直于y轴的边界，对ThetaX和ThetaY进行叠加后，模型体1,2边界的数值绝对值会达到模型体3边界的2倍左右，因此模型体3的边界较为模糊.因此，ED方法识别边界的结果中，当待识别区域内存在垂直或平行于坐标轴的地质体和与坐标轴成一定角度的地质体时，与坐标轴垂直或平行的地质体的边界较模糊，不利于边界信息的提取.

图6 边界识别结果 (a) ThetaX； (b) ThetaY.Fig.6 Edge detection results

2 方法改进

为了减小地质体走向对ED方法的影响，本文提出IED方法，表达式为

IED(i,j)=max(ThetaX(i,j),ThetaY(i,j)),(6)

式中，i、j分别为点线号，将ThetaX和ThetaY逐点对比，取较大的值赋给IED.这样可以避免与坐标轴成一定角度的边界叠加两次.但这样运算会使ThetaX和ThetaY中的较弱的虚假异常凸显出来.为了减弱这种效应，本文根据ThetaX和ThetaY各自的最小值TXmin和TYmin，定义了阈值T：

T=α·max(TXmin,TYmin).

(7)

由于ThetaX和ThetaY都为负值，T也为负值，α取值的范围为0～1.设定判定条件：当测点的ThetaX和ThetaY的值都小于T，则IED取二者较小的值；其他情况，IED取二者较大的值，即：

IED(i,j)=

(8)

将其应用到模型试验中，令α分别为0.9、0.6、0.3，并与EI方法及常用的斜导数法、Theta法进行对比(图7).

由图7，对比几种常用方法及α分别为0.9、0.6、0.3的IED方法边界识别结果，可以看出，斜导数法和Theta方法都会在正负异常之间产生虚假边界，从而影响边界识别效果.IED方法对三个模型体的边界识别效果都较明显，一定程度上改善了ED方法中模型体3相比于模型体1、2边界模糊的现象，而且随着α值的减小，虚假异常也逐渐减弱，但是当α值减小到一定程度，边界的有用信息也开始有一定的损失，通过实验对比分析得到当α值取0.6～0.8时效果较好.根据图1中虚线位置提取斜导数方法，Theta方法，ED方法及α值为0.6的IED方法边界识别结果剖面图如图8所示.

图7 几种常用方法及不同α值的IED方法边界识别结果 (a) 斜导数法； (b) Theta角法； (c) ED法； (d) IED of α=0.9； (e) IED of α=0.6； (f) IED of α=0.3.Fig.7 Edge detection results of several commonly used methods and IED with different values of α (a) Tilt derivative method； (b) Theta angle method； (c) ED method； (d) IED with α=0.9； (e) IED with α=0.6； (f) IED with α=0.3.

图8a是斜导数方法边界识别结果的剖面图，该方法用零值线(图中虚线)识别边界，图中零值点与地质体边界对应较好，但是在正负异常之间有多余的零值点d，该点显示为虚假边界.图8b是Theta角法边界识别结果的剖面图，与图8a类似，在正负异常之间也有虚假边界存在(f点).图8c是ED方法边界识别结果的剖面图，从图中可以看出该方法能够较好的识别出地质体边界，且没有虚假边界产生，但是左侧异常体边界处的值是右侧异常体边界处值的2倍以上，如果地下地质构造情况较复杂，有些地质体的边界可能无法识别出来.图8d是IED方法边界识别结果的剖面图，可以看出该方法能够较好的识别出地质体边界，没有虚假异常产生，且异常体边界处的值几乎相等.为进一步讨论该方法的实用性，设计了加干扰的复杂模型，模型参数如表中所示.

图8 几种方法边界识别结果剖面图 (a) 斜导数法； (b) Theta角法； (c) ED法； (d) IED α=0.6.Fig.8 Edge detection resultsby methods shown by profile (a) Tilt derivative method； (b) Theta angle method； (c) ED method； (d) IED with α=0.6.

图9为复杂模型的示意图，表2为复杂模型参数.在模型试验中加入异常最大幅值1%的随机噪声.用IED方法进行边界识别，结果如图10所示.

图9 复杂模型示意图Fig.9 Schematic diagram of complex model

从图10的处理结果可以看出，对于加入噪声的复杂模型，IED方法也能较好地识别出边界.边界识别结果在模型体边界位置受噪声影响最小，则该方法在地质体边界处对干扰的压制作用最强，这一特点在边界识别中是有利的.

图10 IED方法边界识别结果(α=0.6)Fig.10 Edge detection result of IED(α=0.6)

3 实际数据处理

本文实际数据是加拿大自然资源部(Natural Resources Canada)提供的圣乔治湾的(St.George′s Bay)全张量重力梯度测量数据(Miller and Singh，1994)，包含88条测线，线距为400 m，总飞行长度为6602 km，测区覆盖面积为3056 km2，飞行高度为80 m，联络线距为5000 m.该地区测得的重力全张量数据，如图11所示图，其中黑色实线是已知的地质边界.

图11 圣乔治湾地区重力全张量 (a) gxx; (b) gyy; (c) gzz; (d) gxz; (e) gyz; (f) gxy.Fig.11 Fullgravity tensors in St. George′s Bay, Canada

测区位于圣劳伦斯在纽芬兰岛西南海岸的海湾，测区包含了圣乔治湾的大面积石炭系凹陷，石炭系和上泥盆统的沉积岩在北、东、南三个方向出露，覆盖在前寒武纪基底上.测区中分布着不同时代的地层，可以通过重力梯度全张量数据对不同密度的地质边界进行边界识别.图12为利用ED方法(图12a、b)和α值为0.7时的IED方法(图12c、d)进行边界识别的图像.

对比图12a、c，可以看出两种方法都能识别出该地区大的构造边界.ED法和IED法在对断裂F1、F2、F4、F6、F7的识别效果是相近的，因为这5条断裂与水平方向存在一定的夹角，而对于断裂F3与F5，对比图中红色虚线圈闭框(B和C)标示出来的区域(即F5断裂)，会发现ED方法一些边界由于受走向的影响在图中只有较弱的显示，边界不连续，点状分布，而IED方法能清楚地分辨出边界的位置；对比红色虚线框(A)标示出来的区域(即断裂F3位置)，可以看出一些边界通过ED方法甚至无法识别出来，但通过IED方法却能得到清楚连续的边界，可以对F3断裂进行更好的解释.通过对实际数据的处理发现，与ED方法对比，IED方法在实际数据处理中边界识别更加清晰，细节处更加明显，处理结果更加可靠.

图12 圣乔治湾重力梯度全张量数据处边界识别效果 (a)和(b) ED方法； (c)和(d) IED方法.Fig.12 Edge detection results of full tensor gravity gradient data in the St. George′s Bay (a)and(b) ED method； (c)and(d) IED method.