利用加强水平方向总水平导数对位场全张量数据进行边界识别

袁园， 黄大年， 余青露

1 吉林大学地球探测科学与技术学院， 长春 130026 2 国家海洋局第二海洋研究所， 杭州 310012 3 国家海洋局海底科学重点实验室， 杭州 310012 4 中国石化石油物探技术研究院， 南京 211103



利用加强水平方向总水平导数对位场全张量数据进行边界识别

袁园1,2,3， 黄大年1*， 余青露4

1 吉林大学地球探测科学与技术学院， 长春 130026 2 国家海洋局第二海洋研究所， 杭州 310012 3 国家海洋局海底科学重点实验室， 杭州 310012 4 中国石化石油物探技术研究院， 南京 211103

位场全张量梯度数据以其信息量大、含有更高频的信号成分，能更好地描述小的异常特征等优点在地球物理领域中得到广泛应用.边界检测是位场解释中不可缺少的任务，需要新的边界探测器来处理位场梯度张量数据.为了充分利用位场梯度张量数据的多信息成分，本文定义了方向总水平导数和加强方向总水平导数，并利用其定义新的边界检测器.为了能同时显示不同振幅大小异常的边界，本文对其进行了归一化处理.通过模型试验，证明了归一化方法能更加清晰准确地显示浅部和深部的地质体边界信息.最后将该边界检测方法用于加拿大圣乔治湾实际测得全张量重力梯度数据和中国朱日和地区的磁异常数据中，并得到了较好的边界检测结果.

边界检测; 位场梯度张量; 加强方向总水平导数; 归一化

1 引言

边界检测是位场数据解释中的重要组成部分.目前，有许多方法用于边界检测和边界加强.水平导数和垂向导数经常用于加强这些潜在特征.垂向导数通常被用于描述密度体或磁性体的边界(Evjen, 1936).后来人们发现总水平导数和解析信号振幅的最大值直接描述密度和磁化强度边界的检测(Cordell and Grauch, 1985; Cordell, 1979; Nabighian, 1984).Hsu等(1996)提出加强解析信号方法进行到高阶导数以增加该方法的求解能力.Cooper和Cowan (2004)使用不同阶的垂向导数进行增强地球物理异常的细节.

Miller和Singh(1994)使用垂向导数与总水平导数的比来进行边界检测.它是第一个均衡滤波器，能均衡不同振幅强度的异常信息.Rajagopalan和Milligan (1995)使用自动增益控制技术均衡不同大小的振幅.然而，这种方法与窗口大小的选择相关.Verduzco等(2004)用倾斜角的总水平导数(THDR)来进行边界检测.Wijns等(2005)提出了Theta图滤波器进行边界探测，它是用解析信号振幅来对总水平导数进行归一化处理.这些均衡边界探测器虽然能同时显示不同振幅大小的异常的边界信息，但识别出来的边界分辨率较低.马国庆等(2012)利用高阶导数来定义均衡滤波器，提高了对地质体边界的分辨率.但是，当测量的数据中同时包含正负异常体的情况下，这些均衡探测器可能会引入错误的边界信息.

近十几年，位场梯度张量数据被广泛使用.位场梯度张量数据可以通过测量或者数值计算得到，它提供了9个张量成分，包含了比传统位场数据更高频的信号成分，因此，对于梯度张量数据的解释能得到更高分辨率和更加详细的地质构造.这要求新的方法来解释梯度张量数据，尤其是对小地质构造的边界检测.Cooper 和Cowan (2006)以及Oruç和Keskinsezer(2008)提出了方向倾斜角来探测边界，但是只有垂直方向的倾斜角，即传统的倾斜角方法，能得到较好的检测结果.Beiki(2010)和Mikhailov等(2007)定义了方向解析信号来处理重力梯度张量数据.Beiki(2010)利用方向解析信号进行边界检测，比传统的解析信号方法能得到更好的结果，但是该方法也不能同时描述大振幅和小振幅异常体.近年来，一些学者利用重力梯度张量数据的特征值来进行边界识别(Oruç et al., 2013; Sertcelik and Kafadar, 2012; Zhou et al., 2013).

Marson和Klingele(1993)指出总水平导数比解析信号有更高的分辨率.因此，本文定义了方向总水平导数和加强方向总水平导数，并根据其定义不同的边界探测器.为了同时显示不同振幅大小异常的边界，给出了归一化方法.

2 方向总水平导数

位场梯度张量数据为位场U在空间上的二阶偏导数，张量矩阵的形式：

(1)

Beiki(2010)和Mikhailov等(2007)定义了方向解析信号，其表达式为

(2)

(3)

(4)

其中，下标x,y和z表示方向.

这里，我们在x，y和z方向上定义方向总水平导数，其表达式分别为

(5)

(6)

(7)

其中，下标x,y和z表示方向.THDx的最大值描绘地质体的N-S向边界；THDy的最大值描绘地质体的E-W向边界；垂直方向上的总水平导数THDz与传统使用的总水平导数相同，振幅的最大值描绘地质体的边界.

图1比较分析了方向解析信号振幅和方向总水平导数的边界识别能力.图1a为一垂直棱柱体产生的重力异常.图1b—1d为x,y和z方向上的方向解析信号振幅；图1e—1g为x,y和z方向上的方向总水平导数.通过比较可知，方向总水平导数相对于方向解析信号振幅对于边界识别有着更高的分辨率.因此，根据水平方向的总水平导数来定义边界探测器有着更高的分辨率.利用THDx和THDy定义地质体边界探测器：

(8)

EDT的最大值定义地质体的边界.

3 加强方向总水平导数

Hsu等(1996)应用高阶导数定义加强解析信号振幅来提高边界的分辨率.为了提高边界检测的分辨率，我们使用高阶导数来定义方向总水平导数，叫做加强方向总水平导数.

x，y和z方向上的加强方向总水平导数分别为

(9)

(10)

(11)

其中，下标x,y和z表示方向.上标n表示阶数，大于等于1.当n=1时，与方向总水平导数相同.

同样，加强方向总水平导数ETHDx的最大值描绘地质体的N-S向边界；ETHDy的最大值描绘地质体的E-W向边界.因此，我们可以根据加强方向总水平导数来定义地质体边界探测器：

(12)

EEDTn的最大值定义地质体的边界.通常情况下，我们只取n=1,2.因为位场数据高阶导数的计算将放大噪声的影响.为了识别出来的边界更加准确，我们应避免高阶导数的计算.图1h—1j为x,y和z方向上的二阶加强方向总水平导数.可以看出，相对于方向总水平导数有着更高的分辨率.这里，我们用EEDT来表示由二阶加强方向总水平导数定义的边界探测器EEDT2.在EEDT中，需要计算二阶加强方向总水平的导数，然而，垂向导数的直接计算会增加噪声的影响.我们用如下表达式来计算二阶ETHDx和ETHDy导数：

(13)

(14)

(15)

4 归一化方法

边界探测器EDT和EEDT都为非均衡滤波器，不能很好地同时显示不同振幅大小的异常的边界.为了避免这一缺点，对其采用了归一化方法.

本文采用位场数据的一阶垂向导数来对EDT进行归一化处理，表达式为

(16)

这里引入常数p，一般大小在0.001到0.1.常数p的引入是为了避免产生额外错误的边界.大的p值将降低均衡效果，小的p值将加强弱振幅异常的边界.

边界探测器EEDT的归一化表达式为NEEDT=

(17)

同样，参数p与(13)式相同.这里，用如下表达式来计算Uzzz：

(18)

为了验证本文方法的可行性，我们选用两种常用的方法进行对比.它们是倾角总水平导数(THDR)和Theta图(Verduzcoetal., 2004;Wijnsetal., 2005).

5 模型试验

建立如下重力异常模型一：模型中包含4个垂直棱柱体，各个棱柱体的上顶埋深分别为1km，1km，2km，2km.图2显示了模型体的平面图和3维图.各个棱柱体的厚度都为1km，剩余密度为200kg·m-3.图3为合成的重力异常模型.重力异常模型的梯度张量数据可以通过数值计算和实际测量得到.这里，合成模型的重力梯度张量值通过理论公式计算得到(Forsberg, 1984)(图4a—4c).图4d—4i分别为EDT，EEDT，NEDT，NEEDT，THDR和Theta图的边界结果.可以看出，EDT和EEDT为非均衡滤波器，对于深部异常体的边界检测能力较弱，显示结果不清晰.NEDT和NEEDT相对于两种常用的边界探测器THDR和Theta图法更加清晰和准确.

为了进一步验证方法的有效性，建立更加复杂的重力异常模型二，其中包含正负异常，更加符合实际地质情况.模型体的大小和位置与图2相同.唯一不同的是棱柱体1和3的剩余密度为200kg·m-3；棱柱体2和4的剩余密度为-200kg·m-3.合成的重力异常见图5.图6为计算的梯度张量和相应的边界检测结果.可以看出，当同时包含正负异常时，传统的边界探测器THDR和Theta图虽然能很好地识别地质体的边界，但是引入了一些额外的错误边界.然而，本文提出的均衡探测器能很好地识别异常体的边界，没有引入错误的边界信息，能更好地完成边界识别工作.其中，NEDT中的p值为0.1；NEEDT中的p值为0.08.

图2 模型体平面图(a)和三维图(b)Fig.2 Plan view (a) and 3D view (b) of synthetic model

图3 合成重力异常图,棱柱体的剩余密度为200 kg·m-3Fig.3 Synthetic gravity anomaly model. The contrasted densities of prisms are 200 kg·m-3

为了验证方法的稳定性，对图6中模型二的各个重力梯度值增加5%的高斯白噪声，图7a—7c为加入噪声后的重力梯度值.图7d—7i为噪声数据的边界检测结果.其中，NEDT中的p值为0.1，NEEDT中的p值为0.0.001.可以看出，EEDT相对于EDT受噪声的影响较大，主要是由于EEDT利用了位场数据的三阶导数，放大了噪声的影响；由于噪声的干扰，THDR已经无法给出地质体的边界，而NEEDT和Theta图法受到噪声的影响较大.相对而言，NEDT抗噪能力强，更加稳定.

图5 合成的重力异常，棱柱体1和3的剩余密度为 200 kg·m-3，棱柱体2和4的剩余密度为-200 kg·m-3Fig.5 Synthetic gravity anomaly model. The contrasted densities of prisms 1 and 3 are 200 kg·m-3, the contrasted densities of prisms 2 and 4 are -200 kg·m-3

图6 同图4，但为模型二的重力梯度和不同方法的边界识别结果Fig.6 Same as Fig.4, but for gravity gradients and the edges detected by different methods of model two

图7 同图4，但为加入噪声后的模型二的重力梯度和不同方法的边界识别结果Fig.7 Same as Fig.4, but for gravity gradients and the edges detected by different methods of model two add 5% Gaussian noise

图8 同图4，但为加拿大圣乔治湾地区全张量重力梯度值边界检测结果Fig.8 Same as Fig. 4, but for edges of full tensor gravity gradients anomalies in St. Georges Bay, Canada

图9 朱日和地区磁异常边界检测结果 (a)化极后的磁异常数据；(b)异常的EDT结果；(c)异常的EEDT结果；(d)异常的NEDT结果； (e)异常的NEEDT结果；(f)异常的THDR结果；(g)异常的Theta图结果.Fig.9 Edge results of magnetic anomalies in the Zhurihe area (a) Reduction to pole of magnetic data; (b) EDT of the data in (a); (c) EEDT of the data in (a); (d) NEDT of the data in (a); (e) NEEDT of the data in (a); (f) THDR of the data in (a); (g) Theta map of the data in (a).

利用本文提出的方法对磁异常数据进行处理前，需要对磁异常数据进行化极处理，因为磁数据及其导数均受磁化方向的干扰.因此采用化极后的数据进行解释所获得的结果将更加准确(Li, 2006).

6 应用到实际数据

为了验证方法在实际中的应用效果，分别对加拿大圣乔治湾的全张量重力梯度张量数据和中国朱日和地区的磁异常数据进行边界检测.

图8a—8c为Air-FTG全张量重力梯度仪在圣乔治湾测得的重力梯度张量数据，飞行线距为500 m，联络线距为5000 m，飞行的高度为100 m；图8d—8i为识别的边界结果.可以看出，非均衡边界探测器EDT和EEDT的识别能力较差，识别的边界结果不够清晰，且不能同时显示不同振幅大小的异常的边界.常用的均衡探测器THDR的显示能力较差，不能显示地质体边界.Theta图法可以很好地显示地质体的边界，但是识别出来的边界并不十分清晰.而NEDT和NEEDT识别出来的边界相对于Theta图的结果更加清晰、准确，尤其是NEEDT识别出来的边界更清晰准确.其中，NEDT中的p值为0.005，NEEDT中的p值为0.001.

图9a为朱日和地区化极后的磁异常数据，勘探区域面积为73 km×117 km，采样间隔为20 m，飞行高度为1500 m，飞行线距为500 m.勘探区域为新元古代超层序，主要由大陆沉积物组成，除一些富含铁的沙堰堤坝外基本不含任何铁磁性物质.磁异常数据主要由SE-NW向的堤坝产生的近线性异常和铁异常构成.图中的黑色点线为已知的构造.相应的梯度值通过数值计算得到(Minkus and Hinojosa, 2001).图9b—9g为相应的边界检测结果.可以看出，NEDT和NEEDT的边界识别效果明显好于其他方法的，能更加清晰地显示边界信息.其中，NEDT和NEEDT中的p值为0.005.

7 结论

本文针对位场梯度张量数据包含多个信息成分的优势，提出方向总水平导数和加强方向总水平导数.根据方向总水平导数和加强方向总水平导数建立了两个新的边界探测器，并对其进行归一化处理，以便能同时检测不同振幅大小异常的边界.当异常数据同时包含正负异常时，本文提出的新探测器能很好地避免产生错误的边界信息，而传统的THDR和Theta图将产生额外的错误边界信息.通过模型试验和实际数据，证明了基于方向总水平导数和加强方向总水平导数定义的边界探测器能很好地识别地质体的边界，且识别出来的边界更加清晰、准确，避免了错误边界信息的引入.

致谢 感谢加拿大自然资源部(Natural Resources Canada)提供圣乔治湾(St. Georges Bay)的Air-FTG全张量重力梯度测量数据，并允许发表数据处理结果.感谢为评审本文所付出努力的专家.

Beiki M. 2010. Analytic signals of gravity gradient tensor and their application to estimate source location.Geophysics, 75(6): I59-I74.

Cooper G R J, Cowan D R. 2004. Filtering using variable order vertical derivatives.Computer&Geosciences, 30(5): 455-459.

Cooper G R J, Cowan D R. 2006. Enhancing potential field data using filters based on the local phase.Computer&Geosciences, 32(10): 1585-1591.

Cordell L. 1979. Gravimetric expression of Graben faulting in Santa Fe country and the Espanola Basin, New Mexico. New Mexico Geol. Soc. Guidebook, 30thField Conf., 1979: 59-64.

Cordell L, Grauch V J S. 1985. Mapping basement magnetization zones from aeromagnetic data in the San Juan basin, New Mexico. ∥Hinzc W J ed. The Utility of Regional Gravity and Magnetic Anomaly. Society of Exploration Geophysics, 181-197.

Evjen H M. 1936. The place of the vertical gradient in gravitational interpretations.Geophysics, 1(1): 127-136.

Forsberg R. 1984. A study of terrain reductions, density anomalies and geophysical inversion methods in gravity field modelling. Dept. of Geodetic Science and Surveying, Report 355, OhioState University.

Hsu S H, Sibuet J C, Shyu C T. 1996. High-resolution detection of geologic boundaries from potential-field anomalies: An enhanced analytic signal technique.Geophysics, 61(2): 373-386. Li X. 2006. Understanding 3D analytic signal amplitude.Geophysics, 71(2): L13-L16. Ma G Q, Huang D N, Yu P, et al. 2012. Application of improved balancing filters to edge identification of potential field data.ChineseJ.Geophys. (in Chinese), 55(12): 4288-4295, doi: 10.6038/j.issn.0001-5733.2012.12.040.

Marson I, Klingele E E. 1993. Advantages of using the vertical gradient of gravity for 3-D interpretation.Geophysics, 58(11): 1588-1595. Mikhailov V, Pajot G, Diament M, et al. 2007. Tensor deconvolution: A method to locate equivalent sources from full tensor gravity data.Geophysics, 72(5): 161-169.

Miller H G, Singh V. 1994. Potential field tilt—a new concept for location of potential field sources.JournalofAppliedGeophysics, 32(2-3): 213-217. Minkus K L, Hinojosa J H. 2001. The complete gravity gradient tensor derived from the vertical component of gravity: a Fourier transform technique.JournalofAppliedGeophysics, 46(3): 159-174.

Nabighian M N. 1984. Toward a three-dimensional automatic interpretation of potential field data via generalized Hilbert transforms: Fundamental relations.Geophysics, 49(6): 780-786.

Oruç B, Keskinsezer A. 2008. Structural setting of the northeastern

Biga Peninsula (Turkey) from tilt derivatives of gravity gradient tensors and magnitude of horizontal gravity components.PureandAppliedGeophysics, 165(9-10): 1913-1927. Oruç B, Sertçelik I, Kafadar Ö, et al. 2013. Structural interpretation of the Erzurum Basin, eastern Turkey, using curvature gravity gradient tensor and gravity inversion of basement relief.JournalofAppliedGeophysics, 88: 105-113. Rajagopalan S, Milligan P. 1994. Image enhancement of aeromagnetic data using automatic gain control.ExplorationGeophysics, 25(4): 173-178.

Sertcelik I, Kafadar O. 2012. Application of edge detection to potential field data using eigenvalue analysis of structure tensor.JournalofAppliedGeophysics, 84: 86-94.

Verduzco B, Fairhead J D, Green C M, et al. 2004. New insights into magnetic derivatives for structural mapping.TheLeadingEdge, 23(2): 116-119.

Wijns C, Perez C, Kowalczyk P. 2005. Theta map: edge detection in magnetic data.Geophysics, 70(4): L39-L43.

Zhou W N, Du X J, Li J Y. 2013. The limitation of curvature gravity gradient tensor for edge detection and a method for overcoming it.JournalofAppliedGeophysics, 98: 237-242.

附中文参考文献

马国庆, 黄大年, 于平等. 2012. 改进的均衡滤波器在位场数据边界识别中的应用. 地球物理学报, 55(12): 4288-4295, doi: 10.6038/j.issn.0001-5733.2012.12.040.

(本文编辑 胡素芳)

Using enhanced directional total horizontal derivatives to detect the edges of potential-field full tensor data

YUAN Yuan1,2,3, HUANG Da-Nian1*, YU Qing-Lu4

1CollegeofGeoExplorationScienceandTechnology,JilinUniversity,Changchun130026,China2SecondInstituteofOceanography,StateOceanicAdministration,Hangzhou310012,China3KeyLaboratoryofSubmarineGeoscience,StateOceanicAdministration,Hangzhou310012,China4SinopecGeophysicalResearchInstitute,Nanjing211103,China

With the ongoing development of the full potential field gradient techniques, more and more gradient tensor data have been widely used in geophysical exploration for its large amount of information and containing higher frequency signals than potential field data. The full potential field gradient technique can simultaneously measure six gradient components. Each component has its own geophysical meaning. The high frequency gradient tensor data can be used to delineate small scale anomalies.Edge detection is required for interpretation of potential field data, and has been used in exploration technology for discovery of mineral resources. The main geological edges are fault lines and the borders of geological or rock bodies of different density, magnetic nature, etc. On account of the higher-frequency signals, their interpretation enables high-resolution and detailed investigation of geological structures. Development of new methods is required to enable interpretation of these data, especially for edge detection for small geological structures.In order to make use of multiple component information, we define directional total horizontal derivatives and enhanced directional total horizontal derivatives and use them to define new edge detectors. However, the new defined edge detectors can balance the detected edge signals of different amplitude anomalies. In order to balance the edge signal amplitude, we present a normalization method, where we divide the new edge detectors by the vertical potential field gradient data. Besides, we introduce a constant parameter in the denominator of the normalization method, which can effectively avoid bringing some additional false edges when real geological bodies contain positive and negative anomalies simultaneously.These methods have been tested with synthetic data to verify that the new methods can delineate the edges of different amplitude anomalies clearly. The results show that new defined directional total horizontal derivatives and enhanced directional total horizontal derivatives have higher resolution than directional analytic signal. The normalized edge detectors can display the edges of large and small amplitude anomalies simultaneously, and avoid introducing additional false edges. To further test the stability, we demonstrate the new edge detectors with the model data corrupted with 5% Gaussian noise. Finally, we apply these methods to real full gravity gradient tensor data in St. Georges Bay, Canada and magnetic anomalies in Zhurihe area, China, and get good edge results.

Edge detection; Potential field gradient tensor; Enhanced directional total horizontal derivative; Normalization

10.6038/cjg20150730.

国家高技术研究发展计划(863计划)课题(2014AA06A613)，国家重点基础研究发展计划项目(973项目)西南印度洋洋中脊热液成矿过程与硫化物矿区预测的第五课题“硫化物矿区特征和找矿标志”(2012CB417305)，大洋十二五项目西南印度洋脊合同区多金属硫化物资源评价联合资助.

袁园，男，1988年生，助理研究员，主要从事位场数据处理及解释方面的研究.E-mail：yuanyuan12@mails.jlu.edu.cn

*通讯作者 黄大年，男，1958年生，教授，博士生导师，主要从事移动平台探测数据处理与解释及一体化软件平台开发. E-mail:dnhuang@jlu.edu.cn

10.6038/cjg20150730

P631

2015-02-22，2015-06-16收修定稿

袁园， 黄大年， 余青露. 2015. 利用加强水平方向总水平导数对位场全张量数据进行边界识别.地球物理学报，58(7):2556-2565,

Yuan Y, Huang D N, Yu Q L. 2015. Using enhanced directional total horizontal derivatives to detect the edges of potential-field full tensor data.ChineseJ.Geophys. (in Chinese),58(7):2556-2565,doi:10.6038/cjg20150730.