赵建宇， 明彦伯， 孙成城， 张志东， 刘晓甲

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

重磁(梯度)张量数据边界识别方法研究

赵建宇， 明彦伯， 孙成城， 张志东， 刘晓甲

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

边界识别是重力数据解释的常规任务之一，但地质体的边界不能很好地与重力原始异常数据对应，而与重磁异常水平导数极大值、垂直导数零值相对应，因此大多利用该性质完成重磁边界识别任务。研究发现，现有边界识别方法会出现虚假边界，针对现有边界识别方法进行总结，并针对多余边界问题进行改进，提出新的边界识别，通过模型试验和实际数据证明，该方法可有效地去除多余边界，且降低了噪声的干扰，为重磁数据边界识别方法提供新的思路。

边界识别； 重磁(梯度)张量数据； 重磁勘探

0 引言

重磁勘探是一种可以有效地圈定异常，划分构造的物理勘探方法，由于它的经济、快速、范围广等特点而得到广泛应用。但是在实际生产工作中能够发现，重磁原始异常数据和地质体的边界不能很好地对应[1]，在异常边界识别过程中可能会出现虚假的异常，给研究和生产造成不必要的麻烦。于是学者们开始对重磁异常数据的边界识别进行研究，垂直导数[2-3]是最早被用来进行边界识别的方法，后来又发现水平导数的极值位于密度和磁化率变化的位置[4-5]，从这点出发，同样证明了水平导数对于边界识别的作用。但是由于导数随地质体埋深增加，其数值会迅速衰减[6]，因此，垂直导数和水平导数不能很好地识别深部异常。为了解决这一问题，学者将目光转向对均衡滤波器的研究与开发，Miller等[7]提出了第一个均衡边界识别滤波器—tilt angle；Verduzco等[8]提出利用tilt angle的总水平导数(THDR)进行边界的识别；Wijns等[9]利用总水平导数与解析信号的比值(theta map)来进行边界识别操作。这些方法有的可以均衡强弱异常之间的幅值，有的可以进行边界识别。

重磁张量的研究是近年来兴起的一种高精度勘探手段，目前国内、外只有针对原始异常边界识别方法的研究，并未开展张量数据的研究，以及张量数据加权方式对边界识别结果精度影响的讨论，而利用重磁张量边界识别技术可准确地获得地质水平位置信息；利用延拓或者相关方法来进行张量数据的边界识别，获得地质体的空间位置信息。通过研究重磁张量数据的处理方法，进一步开发适用于全张量重磁数据的处理手段，并解决现有边界识别方法所获得的地质体边界范围较大，以及不能获得深度信息的问题，提出新思路。在这一研究背景下，笔者分析了各种边界识别方法，比较优缺点，并加以改进，整理出优化方案，最后得到新方法的成果图像。

1 重磁张量数据的特征

张量表示为G的二阶导数，可采用矩阵的形式表示为：

(1)

根据拉普拉斯方程可知gxy=gyx、gxz=gzx、gyz=gzy，张量测量仅有五个是独立的。即：

(2)

根据棱柱体重力张量异常(图1)分析张量属性：gxx能有效地增强地质体的南北向边界；gyy能增强异常的东西向边界；gxy能描述地质体的角点信息；gzx的最大值能识别出地质体的南北边界；gzy能识别地质体的东西边界，不同的张量异常能反映出地质体的不同边界信息。

2 边界识别方法

常用重磁数据边界识别方法如下：

1)Cordell[10]提出总水平导数(THD)。

图1 棱柱体重力张量异常Fig.1 Gravity tensor abnormity of prismoid

(3)

2)Miller and Singh[7]倾斜角法。

(4)

3)Verduzco等[8]提出利用tilt angle的总水平导数(THDR)进行边界的识别。

(5)

4)Wijns等[9]利用总水平导数与解析信号的比值(theta map)来进行边界识别。

(6)

5)Cooper等[11]提出正则倾斜角法来进行异常体边界的识别。

(7)

采用这五种边界识别方法分别进行一正一负密度模型正演后再进行边界识别处理(图2)，通过比较可以看出这五种方法各有利弊。

图2 不同边界识别方法比较Fig.2 Compare with different boundary identification methods(a)原始异常数据;(b)总水平导数;(c)倾斜角法;(d)THDR;(e)theta map;(f)TDX

由图2可以发现：图2(a)原始异常数据对边界识别效果不好，边界非常模糊；图2(b)总水平导数法较原始数据有所加强，但是边界识别效果仍然有些模糊；图2(c)倾斜角法、图2(d)THDR 、图2(e)theta map、图2(f)正则化倾斜角法这四种方法会产生多余边界。

3 改进的边界识别方法

由上述结论发现：Cooper[11]提出正则倾斜角法来进行异常体边界的识别(TDX)对于异常体的边界呈现效果更好，但是却在一正一负密度体模型成像时出现虚假异常边界。所以，可以针对TXD方法进行改进，消除虚假异常边界。

目前，对边界识别方法的改进方面，马国庆等[3]利用垂直导数的水平导数构成均衡边界识别滤波器进行地质体的边界识别；王万银等[2]提出归一化总水平导数垂向导数进行地质体的边界识别，并分析了位场总水平导数极值和垂直导数零值的空间位置变化，总结得到总水平导数极大值随着地质体埋深的加大与地质体边界偏差越大，利用垂直导数的总水平导数进行地质体的边界识别可减小偏差。在这些成果的启发下，需将TDX的公式进行适当的调整 (由于TDX只在一正一负密度模型中的边界识别结果有虚假异常，而在同为正或负的密度模型中没有，因此只讨论前一种情况。)。

首先，分析TXD方法产生多余边界的原因：

(7)

其中f为重磁异常值。如式(7)所示，该方法采用极大值法来进行地质体边界的识别，即当垂直导数为零时取得极大值对应地质体边界。而在一正一负密度模型中的垂直导数Vzz为“0”的点有5处(图3)，其中四处在地质体边界处，另外一处因为正密度地质体Vzz的负值和负密度地质体的Vzz正值相叠加产生的“0”，即为多余边界产生的原因(其中Vzz为∂f/∂z，即垂直导数)。

图3 重力垂直导数剖面图Fig.3 Vertical derivative of gravity data profile

对TDX进行改进，提出一种基于TDX改进的新方法CTDX(Change TDX)：

(8)

由式(8)可以看出：CTDX是在TDX的基础上在其分母加一个Δ值。使原有方法产生多余边界处的总水平导数与式(8)分母的比值趋近于零值；而在真实的地质体边界处总水平导数与式(8)分母的比值不为零值。这样，可以使TDX中产生的多余边界和地质体边界分离开。

对TDX与CTDX方法进行比较(图4)。

图4 TDX和CTDX比较Fig.4 Comparison between TDX and CTDX(a)TDX ；(b)CTDX

从图4可以看出，CTDX没有TDX里的虚假异常边界，并且其边界识别能力也没有得到影响。说明CTDX的方法可行(其中小数Δ取值0.01)，但是对于小数Δ的取值还有待商榷。

对于小变量或定量Δ，不同的取值会对边界的识别产生不同的效果。因此对相同模型不同Δ取值进行了实验：

Δ分别取0、0.000 1、0.001、0.01、0.1、1，以及总水平导数(图5)。

从图5可以得出结论,Δ太小依然会产生多余边界，Δ增大会使结果逐渐趋近与总水平导数，而总水平导数对较深地质体效果不好。Δ的取值将受不同地质条件产生变化，所以Δ值的确定是目前研究的难点和重点。

以上Δ取值为具体数值，具有一定的经验性，为了将Δ的取值问题应用到实际的数据处理过程中，进行了不同Δ计算方式的模型试算：

图5 不同Δ取值结果比较图Fig.5 Comparison with different Δ(a)Δ=0；(b)Δ=0.0001；(c)Δ=0.001；(d)Δ=0.01；(e)Δ=0.1；(f)Δ=1；(g)总水平导数

图6 不同Δ计算方式的模型试算图Fig.6 Comparison with different models of compute Δ(a)采用TDX对模型进行计算;(b)～(h)使用取垂直导数Vzz的最大值和最小值之差的1/2、1/5、1/10、1/20、1/30、1/40、1/50;(i)模型的垂直导数

将图6(b)、图6(c)、图6(d)、图6(e)和图6(i)进行比较发现，若Δ取值太大，对于深处的地质体的边界刻画能力较弱，但是比较图6(f)、图6(g)、图6(h)可以发现，Δ取值太小，对多余边界的削弱程度不够，仍然会有多余边界出现。综合上述对试算模型的结果分析，可以将Δ计算方式定义为式(9)。

Δ=|(Vzzmax-Vzzmin)|/30

(9)

此时的CTDX既可以对较深地质体边界进行较好地识别，也没有产生多余边界，而对两个正密度体模型也有着较好的效果(图7)。

图7 CTDX对不同模型的效果Fig.7 CTDX for different models(a)一正一负密度体模型； (b)两个正密度体模型

4 噪声压制后的效果

通过正演理论模型试验本文方法的抗噪性，在图5正演模型中加入了随机噪声，得到了两个正密度体加噪声模型的处理结果(图8)，一正一负密度体加噪声模型的处理结果(图9)。

图8 两个正密度体加噪声Fig.8 Two positive density body models with noise(a)TDX ；(b)CTDX

图9 一正一负密度体加噪声Fig.9 One positive density body models and one negative density body models with noise(a)TDX ；(b)CTDX

通过对比处理结果可以得出结论：新方法与传统方法相比，能较好地压制噪声和提高信噪比的效果。在随机噪声的干扰下，仍能较好地反映出异常体边界。

5 实际数据的处理

根据研究区四川地区的重磁异常值(图10、图11)，可以发现，虽然该方法可以识别大体上的一些构造单元，但是对于有一些构造单元会产生多余的处理，加大实际解释时的难度。将TDX方法与改进后的CTDX方法(图12)行对比不难发现，TDX方法在处理磁异常数据时效果不是很好，对于异常的识别显得非常杂乱，如果想要详细划分解释必然要经过更多地处理。

分别经过TDX、CTDX边界识别方法处理后，可以得到如图12～图15所示的结果。

图10 研究区原始重力数据图Fig.10 The original gravity data of the area

图11 研究区磁数据图Fig.11 The original magnetic data of the area

图12 TDX方法处理研究区重力数据Fig.12 TDX of the gravity data in research area

图13 TDX方法处理研究区重力数据Fig.13 TDX of the gravity data in research area

图14 TDX方法处理研究区磁数据Fig.14 TDX of the magnetic data in research area

图15 CTDX方法处理研究区磁数据Fig.15 CTDX of the magnetic data in research area

图13相较图12少了一些多余的边界(红色框中的阴影)，但对于已经可以识别的构造单元能更好地识别其边界，并且有效地增强了深度场源体的边界，使构造划分工作变得简单。

比较图15与图14，两者的效果非常明显，可以看出CTDX相较于TDX的处理结果可以更清楚地划分出构造单元。从图15与图13结果结合来看，研究区内应存在一定数量的北东走向地质构造体，依据研究区的地质资料，推测是一系列的断裂。且CTDX方法对于重力数据的改进效果有限，对磁数据的改进效果较为可观。

以上从重、磁两各方面印证了改进的正则化倾斜角法，对于边界识别问题的效果明显，说明改进的正则化倾斜角法可以在一定程度上优化边界识别问题，去除多余边界。

6 结论

总结了前人对异常体边界识别的方法，分析了各个方法的优缺点，并在正则倾斜角法基础上进行改进，从虚假异常边界产生的原因入手，引入Δ=|(Vzzmax-Vzzmin)|/30，经过数据处理后，能够很好地去除掉虚假异常边界，并针对噪声干扰进行模型实验，发现该方法(CTDX)对噪声的压制效果也很可观。通过对实际数据进行处理，可以发现改进后的方法可以去除多余边界，增强深度场源体的边界效果。

[1] 马国庆, 黄大年, 于平, 等.改进的均衡滤波器在位场数据边界识别中的应用[J].地球物理学报, 2012,55(12)：4288-4295.

MA G Q, HUANG D N, YU P, et al. Application of improved balancing filter to edge identification of potential field data[J].Chinese J Geophys, 2012,55 (12): 4288-4295.(In Chinese)

[2] 王万银，潘玉，邱之云.位场数据归一化总水平导数垂向导数边缘识别方法(英文)[J];Applied Geophysics,2009(03)：226 - 233.

WANG WY, PAN Y , QIU ZY . A new edge recognition technology based on the normalized vertical derivative of the total horizontal derivative for potential field data[J]. Applied Geophysics, 2009(3):226 - 233. (In Chinese)

[3] 马国庆, 杜晓娟, 李丽丽. 利用水平与垂直导数的相关系数进行位场数据的边界识别[J]. 吉林大学学报(地球科学版), 2011, 41(S1) 345-348.

MA G Q，DU X J，LI L L. Edge detection of potential filed data using correlation coefficients of horizontal and vertical derivatives[J].Journal of Jilin University( Earth Science Edition), 2011, 41(S1): 345-348.(In Chinese)

[4] 刘艳,李桐林,黄旭钊,等.航磁异常轴向特征提取算法研究及其应用[J].地球物理学进展,2002(1)： 141-144.

LIU Y， LI TL， HUANG X Z, et al. An automatic method of detection of aeromagnetic axes and its application[J].Progress in Geophysics, 2002 (01) :141-144.(In Chinese)

[5] 柏冠军,吴汉宁,赵希刚,等.重力资料识别鄂尔多斯盆地线性构造方法研究[J].地球物理学进展,2007(05):1386-1392.

BAI G J,WU H N,ZHAO X G,et al.Research on recognition of linear structures using gravity data in ordos basin[J].Progress in Geophysics,2007(05): 1386-1392. (In Chinese)

[6] 郭华,吴燕冈,高铁.重力斜导数方法在时间域中的理论模型与研究[J].吉林大学学报(地球科学版),2006(S2): 9-14.

GUO H, WU Y G, GAO T.The research of theories model in time area with the method of title derivative in gravity[J].Journal of Jilin University( Earth Science Edition),2006 (S2): 9-14.(In Chinese)

[7] MILLER H G, SINGH V. Potential field tilt—a new concept for location of potential field sources[J].Journal of Applied Geophysics, 1994, 32: 213-217.

[8] VERDUZCO B., FAIRHEAD, J. D. ,Green, C. M. New insights into magnetic derivatives for structural mapping[J].The Leading Edge, 2004, 23(2):116-119.

[9] EVJEN H M. The place of the vertical gradient in gravitational interpretations[J].Geophysics,1936(1):127-136.

[10] CORDELL L. Gravimetric expression of graben faulting in Santa Fe Country and the Espanola Basin, New Mexico[J]. New Mexico Geol. Soc. Guidebook, 30th Field Conf,1979:59 - 64.

[11] COOPER G R J, Cowan D R. Enhancing potential field data using filters based on the local phase[J]. Computers &Geosciences, 2006, 32: 1585-1591.

Investigationaboutedgerecognitionusing(gradient)tensordataofmagneticandgravity

ZHAO Jianyu , MING Yanbo , SUN Chengcheng, ZHANG Zhidong , LIU Xiaojia

(College of Geoexploration Science And Technology Jilin University , Changchun 130021, China)

The boundary identification is one of the regular tasks of gravity and magnetic data interpretation, but gravity and magnetic raw anomaly data does not well correspond to the boundary of geological body. The maximum value of the horizontal derivative of gravity and magnetic anomalies are corresponding to the zero value of vertical derivative, which is used mostly to complete the assignment of gravity and magnetic boundary identification. We find that false boundary is caused by using current boundary identification method. This article summarizes the current boundary identification method and corrects the problem of excessive boundaries to present new boundary identification method. This new method can effectively remove excessive boundaries and reduces the disturbance of noise, which provides a new thinking way for gravity and magnetic data boundary identification method.

boundary identification; (gradient) tensor data of magnetic and gravity; gravity and magnetic exploration

2016-10-21 改回日期： 2017-03-03

国家重点研发项目(2016YFC0600505)

赵建宇(1993-)，男，硕士，主要研究方向为固体地球物理学，E-mai：lzjy849891627@163.com。

明彦伯(1994-)，男，硕士，主要研究方向为固体地球物理学，E-mai：ybming1994@qq.com。

1001-1749(2017)06-0748-07

P 631.2

A

10.3969/j.issn.1001-1749.2017.06.06