地球物理学进展  2016, Vol. 31 Issue (3): 1164-1172   PDF    
引用本文
王林飞, 郭灿灿, 薛典军, 熊盛青. 2016. 磁梯度张量解析信号分析法及其在场源位置识别中的应用[J]. 地球物理学进展, 31(3): 1164-1172, doi: 10.6038/pg20160333  
WANG Lin-fei, GUO Can-can, XUE Dian-jun, XIONG Sheng-qing. 2016. Analytic signals of magnetic gradient tensor and their application to estimate source location. Progress in Geophysics, 31(3): 1164-1172, doi: 10.6038/pg20160333   
磁梯度张量解析信号分析法及其在场源位置识别中的应用
王林飞1, 郭灿灿2, 薛典军1, 熊盛青1    
1. 中国国土资源航空物探遥感中心, 北京 100083;
2. 中国地质大学(北京)地球物理与信息技术学院, 北京 100083
收稿日期: 2015-08-21 修回日期: 2015-12-02.
基金项目: 国家高技术研究发展计划“863计划”(2013AA063905)、航空物探方法技术与数据处理解释系统研究(1212011087010)和海洋地质保障工程综合数据库建设(航遥中心)(GZH201200513)联合资助.
作者简介: 王林飞,女,1977年生,浙江宁波人,博士,高级工程师,现就职于中国国土资源航空物探遥感中心,主要从事航空物探方法技术研究、软件开发及信息系统建设工作.(E-mail:wlfei4103@126.com)
摘要: 相对于传统磁异常数据,磁梯度张量数据可以提供有关异常体更多的信息,分辨率更高,具有广阔的应用前景.本文给出了基于磁梯度张量解析信号进行边界增强的新方法,该方法是通过对异常垂向一阶导数与z方向解析信号的比值做归一化得到的,其零值线位置可以很准确的反映出地质体边界位置,而且可以有效降低噪声的干扰.文中还证明了异常垂向一阶导数与x、y、z三个方向解析信号的比值同样满足欧拉齐次方程,且在计算过程中不需考虑构造指数N的影响,避免了因构造指数不当而引起的反演误差.通过对单一地质体及组合地质体模型的实验证明:与常规欧拉反褶积法相比,该反演方法能够更好的得到地质体界面及深度信息,所得的解更集中.将其应用到保定实测数据中,获得了更精确的场源信息.
关键词: 磁梯度张量     欧拉反褶积     解析信号     边界识别    
Analytic signals of magnetic gradient tensor and their application to estimate source location
WANG Lin-fei1, GUO Can-can2, XUE Dian-jun1, XIONG Sheng-qing1    
1. China Aero Geophysical Survey & Remote Sensing Center for Land and Resources, Beijing 100083, China;
2. School of Geophysics and Information Technology, China University of Geosciences, Beijing 100083, China
Abstract: Magnetic gradient tensor data can provide more information about the anomalies with higher resolution in contrast to the traditional gravity anomaly data, having vast application prospect. In this paper, we present a new edge identification method based on the analytic signal of magnetic gradient tensor, which can be derived by normalization of the ratio of first derivative and analytic signal in z directions. This method can reflect the edge of causative bodies exactly and reduce the noise interference effectively. I also demonstrate in this paper that the ratio of first derivative and analytic signal in x, y, z directions are homogeneous and satisfy Euler's homogeneity equation, thus a linear equation can be acquired. This method allows the estimation of location parameters without any prior information about the nature (structural index) of the source, avoiding inversion errors that conventional Euler deconvolution method have caused by improper choice of structural index N. The presented method is tested on single and synthetic magnetic anomalies, and the inversion results show that the method can provide the location parameters more accurately compared with conventional Euler deconvolution method. And we obtained more accurate information sources by applying it to the real magnetic data in Baoding area.
Key words: magnetic gradient tensor     Euler deconvolution     analytic signal     edge identification    
0 引 言

目前，测量地磁场矢量三个分量的变化率，是各国地球物理学家研究和开发的热点.张量元素受地磁场倾向、偏角影响小，由他们计算得到的张量不变量不需要额外的处理就可以很好的描述磁场源，且磁梯度张量反演能够更精确的描述场源体的磁化方向和几何形态(张昌达，2006).相对于磁总场梯度，磁梯度张量有两个优点：

(1) 梯度张量算出的不变量等值线易于解释.

(2) 可利用偶极子-追踪算法精确确定磁偶极子的深度和水平位置.Wynn等(1975) 最早提出利用磁力梯度张量数据的偶极子-追踪算法，可精确确定磁偶极子的深度和水平位置.

近年来，随着工程技术及超导量子干涉磁力仪(SQUID)的发展，磁梯度张量数据得到了越来越多的应用(曾邵发和王者江，2003).因此，基于磁梯度张量分析的数据处理与解释方法也得到了不断的发展与改进(吴招才和刘天佑，2008; 张朝阳等，2009; Beiki et al.，2011; 马国庆等，2012).

如何准确的确定场源体的位置及规模是位场数据处理与解释中的一项重要内容，目前常用的自动反演方法有温纳反褶积法(Werner, 1953)，小波变换法(Moreau et al.，1999)，欧拉反褶积法(Thompson, 1982)以及解析信号法(Nabighian, 1972).解析信号法的最大优势在于不需要提供异常体的先验信息，自提出以后获得了广泛应用(Keating and Sailhac, 2004;Ma et al.，2012)并得到了诸多学者的进一步研究与改进(Thurston and Smith, 1997;Bastani and Pedersen, 2001).Nabighian(1984) 又将解析信号理论扩展到三维位场中，提出位场水平两个方向上的梯度与垂向梯度间转换的希尔伯特变换关系，但没有进行实际应用研究；骆遥等(2011) 利用希尔伯特变换提出位场解释的直接解析信号法，并将其应用于航空位场数据的处理.Roest等(1992) 提出利用三维解析信号振幅—总梯度模的极大值识别异常体边界位置.Hsu等(1996) 将Nabighian(1974) 提出的二维数据的增强型解析信号推广到三维，给出了三维数据的n阶增强型解析信号的定义，并用来实现边界检测.Bournas和Baker(2001) 提出利用解析信号振幅总水平导数方法 .Beiki(2010) 提出了利用两个水平方向解析信号的垂向一阶导数的平方和(ED)来进行边界增强，并根据三个方向解析信号的欧拉方程的矩阵形式反演场源位置.以上提出的边界增强方法，均是在原解析信号的基础上求取一阶或多阶导数得到，这势必会对噪声产生放大作用.

本文利用位场垂向一阶导数(即Uzz)与三个方向的解析信号做比值，然后对z方向的比值做归一化得到一种新的边界增强方法，这里我们称之为NR法.同时，本文还证明了垂向一阶导数与磁张量数据的方向解析信号的比值同样满足欧拉齐次方程，该方程消除掉构造指数N的影响，相对于常规欧拉反演算法具有更高的精度.可以很好的用来估算场源的位置及埋深.

2 方法原理

三个相互垂直的轴上的三个磁场分量的空间变化率，共9个值，是磁梯度张量g的要素，公式为

式中：U为磁性体的磁位.从式(1) 中可以看出，磁力梯度张量是对称的，即Uij=Uji.并且，由于在无源区域U遵守Laplace方程，因此G的对角元素之和为零，即Uxx+Uyy+Uzz=0.从这两个性质中可以知道，梯度张量只有5个元素是独立的.

在三维数据中，可以针对矩阵G的每一行定义其解析信号，称为x, y, z方向解析信号(Beiki，2010)，并可写成矩阵形式为

其方向解析信号的幅值可以表示为
我们用rα表示位场垂向一阶导数与方向解析信号的比值(Ratio of first derivative and analytic signal, RDAS)，即
其中，α代表x, y, z三个方向，本文提出对z方向的比值rz做归一化(Normalized ratio of first derivative and analytic signal)处理来边界识别，具体形式为
常规欧拉反褶积是利用位场及其三个方向的导数估算场源的位置，其具体公式(Thompson，1982)为
其中x0，y0，z0是场源坐标，B是背景场，N是构造指数，不同地质体对应不同的N值.

Huang等(1995) 证明，如果函数的构造指数为N, 则相应的解析信号的齐次为N+1，满足方程

因此对于方向解析信号有
我们对方程(8) 求z方向偏导数有
然后对公式(10) 两边乘以fz，对公式(11) 两边乘以Aα，并用前者减去后者有

我们知道，重力异常的垂向一阶导数与重力位的垂向二阶导数相同，所以公式(6) 中的Uzz可以用fz来代替，然后我们对公式(6) 求x, y, z三个方向的偏导数，公式为

把以上三式带入到公式(12) 中有
写成矩阵的形式为
采用该方法可对磁梯度张量数据反演时，不需要考虑磁化方向的干扰且不需要事先知道场源体性质，可自动完成场源位置的计算.

3 模型实验

为了检验本文方法的应用效果，我们首先设计了一个长方体模型，模型中心位置(40 km, 40 km)，顶面埋深1 km, 底面埋深6 km, 磁化强度6 A/m.利用本文边界识别方法(NR法)对该模型产生的磁异常(图 1a)进行处理，并将其与常规边界识别方法，即HGA法(Cordell and Grauch, 1985)以及ED法(Beiki, 2010)进行比较，二者公式为

式中，Ax, z和Ay, z分别代表方向解析信号Ax和Ay的垂向一阶导数.

图 1 长方体解析信号特征及边界识别结果
(a)模型异常;(b)垂向一阶导数与x方向解析信号比值rx号;(c)垂向一阶导数与y方向解析信号比值ry号；(d)垂向一阶导数与z方向解析信号比值rz 号；(e)HGA法识别结果；(f)ED法识别结果；(g)NR法识别结果. Fig. 1 Analytic signal characteristic and edge detection result of single rectangular prism model

图 1b-d显示了垂向一阶导数与方向解析信号比值的特征，图 1e给出了HGA法的识别效果，可以看出，HGA法的极值幅值较宽，不利于边界的确定，ED法(图 1f)和本文方法(图 1g)都可以更精确的显示出边界位置.

为了验证NR法的稳定性，在图 1a异常中分别加入0.1%(图 2a)和0.5%(图 2e)的高斯随机噪声，并对三种方法的检测结果进行对比分析.可以看出，当加入0.1%的的噪声之后，三者的检测结果并没有太大的变化，但是当加入0.5%的噪声之后，本文方法显现出明显的优势，即受噪声干扰较小，仍然可以清晰的反映出边界的位置，随着噪声幅值的增加，并没有受到明显的干扰.

图 2 加入噪声后三种方法对比
(a)加入0.1%的高斯噪声后的异常；(b)加入0.1%后HGA法识别结果；(c)加入0.1%后ED法识别结果；(d)加入0.1%后NR法识别结果；(e)加入0.5%的高斯噪声后的异常；(f)加入0.5%后HGA法识别结果；(g)加入0.5%后ED法识别结果；(h)加入0.5%后NR法识别结果. Fig. 2 Comparison of three methods using anomaly with noise

分别采用常规欧拉反褶积法和本文方法对上述单一模型异常进行反演，在用常规欧拉反褶积法反演时构造指数N值选取0，结果如图 3所示.从图中可以看出：常规欧拉反褶积法得到的解比较分散，难以准确的确定出异常体的边界位置，而且反演的深度与理论值偏差较大；而本文方法对异常体边界位置反映清晰，且反演得到的深度解较为均匀的分布在顶底界面之间，反演的深度均值3.3 km与理论深度均值3.5 km更为接近.

图 3 长方体异常反演结果
(a)常规欧拉反褶积法反演结果；(b)常规欧拉反褶积法的解的分布情况；(c)张量方向解析信号法反演结果；(d)张量方向解析信号法的解的分布情况. Fig. 3 Inversion results of single rectangular prism

下面通过长方体和直立薄板的组合模型验证来本文方法对不同地质体模型的适用性.长方体模型中心坐标(30，40)km, x方向宽度10 km, y方向长度20 km, 顶面埋深2 km, 底面埋深6 km, 直立薄板顶面埋深3 km, 底面埋深10 km, x方向延伸宽度0.4 km, 两个模型的磁化强度均为6 A/m.组合模型磁异常如图 4a.

图 4 长方体与岩脉组合模型异常及反演结果
(a)组合模型异常;(b)组合模型3D视图;(c)常规欧拉反褶积法反演结果;(d)图 4c的3D视图;(e)张量方向解析信号法反演结果;(f)图 4e的3D视图. Fig. 4 Anomaly and inversion results of a rectangular prism and a dike

分别采用常规欧拉法和本文方法对图 4a中所示异常进行反演.由于直立薄板的构造指数N同样为0，因此对组合模型异常做常规欧拉反演时，N取0值，反演结果如图 4c-f.从反演得到的平面位置(图 4c、图 4e)来看，本文方法比常规欧拉法得到的长方体边界位置更为准确，且分散解较少，而且常规欧拉法不能清晰的反映出薄板的位置.从深度图(图 4d、图 4f)来看，本文方法得到解更集中，更为接近模型中心深度，可以清楚的反演出薄板的位置，而常规欧拉法得到的薄板的深度值较为分散，基本得不到薄板的位置信息.

4 应用实例

2011年，中国国土资源航空物探遥感中心在河北保定地区的试验生产中，针对同一测区分别进行了测量高度为200 m(离地高度)和400 m(离地高度)的航磁三轴梯度测量，2个高度层的测量的比例尺均为1:25000(测线间距 250 m).该测区磁场信息丰富，不同高度梯度测量所反映的地质信息高度一致，本文即选取该次试验中低高度测量数据进行试验.该梯度资料覆盖面积约720 km2，测线长35 km, 飞行方向90°或270°，测量的平均离地高度为197 m.我们选取其中一块异常较复杂的区域进行反演计算(图 5)，该区域数据的网格化点距和线距离均为150 m, 可以看出梯度资料分辨的地质细节明显较ΔT丰富.由于梯度测量系统中探头的垂向间距仅为2 m, 故垂向梯度Tz(图 5d)相对噪声水平偏高.在利用两种方法反演之前，均做上延0.2倍点距进行去噪处理.利用磁梯度数据与磁张量数据之间的关系(Blakely, 1995)获得该测区的张量数据(图 6).

图 5 保定测区实测航磁异常及梯度异常
(a)实测磁异常;(b)x方向梯度异常;(c)y方向梯度异常;(d)z方向梯度异常. Fig. 5 Aeromagnetic and gradient anomaly of Baoding area

图 6 保定测区梯度张量异常
(a)磁梯度张量Mxx分量;(b)磁梯度张量Mxy分量;(c)磁梯度张量Mxz分量;(d)磁梯度张量Myy分量; (e)磁梯度张量Myz分量;(f)磁梯度张量Mzz分量. Fig. 6 Magnetic gradient tensor of Baoding area

利用本文方法对保定测区实测异常进行反演，并与常规欧拉反演结果进行对比(图 7).根据异常的形态特征，在进行常规欧拉反褶积反演时构造指数选取0，从反演结果可以看出，常规欧拉法(图 7a)分散解较多且发散，不利于构造划分；而张量解析信号法(图 7b)分散解较少，可以更清晰的反映出该盆地内的构造格局.并且利用本文反演方法计算，不用考虑构造指数的选取不当对反演结果造成的影响.当对测区构造信息掌握不足时，本文方法无疑更为便捷有效.

图 7 反演结果
(a)常规欧拉反褶积法反演结果;(b)张量方向解析信号法反演结果. Fig. 7 Inversion results
5 结 论

本文根据磁梯度张量数据的方向解析信号以及垂向一阶导数的欧拉齐次方程，提出了基于张量数据的方向解析信号欧拉反演方法，该反演在反演公式中消去了构造指数N, 避免了构造指数N的选取不当带来的反演误差.模型试验中，验证了本文方法在位场反演中的有效性，其反演精度明显高于常规欧拉反褶积法.将该方法应用于保定测区实测磁梯度数据，可以得到更准确的断裂信息.可为利用位场数据精确定位地质体边界的平面位置和深度以及划分地质构造单元提供一种新的研究思路.

致 谢 感谢审稿专家提出的修改意见和编辑部的大力支持！

参考文献
[1] Beiki M. 2010. Analytic signals of gravity gradient tensor and their application to estimate source location[J]. Geophysics, 75(6): 159-174.
[2] Beiki M, Pedersen L B. 2010. Eigenvector analysis of gravity gradient tensor to locate geologic bodies[J]. Geophysics, 75(6): 137-149.
[3] Beiki M, Pedersen L B, Nazi H. 2011. Interpretation of aeromagnetic data using eigenvector analysis of pseudo gravity gradient tensor[J]. Geophysics, 76(3): L1-L10.
[4] Blakely R J. 1995. Potential Theory in Gravity and Magnetic Applications[M]. Cambridge: Cambridge University Press.
[5] Bournas N, Baker H A. 2001. Interpretation of magnetic anomalies using the horizontal gradient analytic signal[J]. Annali DiGeofisica, 44(3): 505-526.
[6] Cordell L, Grauch V J S. 1985. Mapping basement magnetization zones from aeromagnetic data in the San Juan basin, New Mexico[C].//Hinze W J ed. The Utility of Regional Gravity and Magnetic Anomaly Maps. SEG, 181-197.
[7] Hsu S K, Sibuet J C, Shyu C T. 1996. High-resolution detection of geologic boundaries from potential-field anomalies: An enhanced analytic signal technique[J]. Geophysics, 61(2): 373-386.
[8] Huang D, Gubbins D, Clark R A, et al. 1995. Combined study of Euler's homogeneity equation for gravity and magnetic field[C].//57th EAGE Conference and Exhibition. Glasgow, UK: EAGE, 144.
[9] Keating P, Sailhac P. 2004. Use of the analytic signal to identify magnetic anomalies due to kimberlite pipes[J]. Geophysics, 69(1): 180-190.
[10] Luo Y, Wang M, Luo F, et al. 2011. Direct analytic signal interpretation of potential field data using 2-D Hilbert transform[J]. Chinese Journal of Geophysics (in Chinese), 54(7): 1912-1920, doi: 10.3969/j.issn.0001-5733.2011.07.025.
[11] Ma G Q, Du X J, Li L L, et al. 2012. Interpretation of magnetic anomalies by horizontal and vertical derivatives of the analytic signal[J]. Applied Geophysics, 9(4): 468-474.
[12] Ma G Q, Li L L, Du X J. 2012. Boundary detection and interpretation by magnetic gradient tensor data[J]. Oil Geophysical Prospecting (in Chinese), 47(5): 815-821.
[13] Moreau F, Gibert D, Holschneider M, et al. 1999. Identification of sources of potential fields with the continuous wavelet transform: Basic theory[J]. Journal of Geophysical Research, 104(B3): 5003-5013.
[14] Nabighian M N. 1972. The analytic signal of two-dimensional magnetic bodies with polygonal cross-section: Its properties and use for automated anomaly interpretation[J]. Geophysics, 37(3): 507-517.
[15] Thompson D T. 1982. EULDPH: a new technique for making computer-assisted depth estimates from magnetic data[J]. Geophysics, 47(1): 31-37.
[16] Thurston J B, Smith R S. 1997. Automatic conversion of magnetic data to depth, dip, and susceptibility contrast using the SPI (TM) method[J]. Geophysics, 62(3): 807-813.
[17] Werner S. 1953. Interpretation of magnetic anomalies of sheet-like bodies[M]. Sveriges Geologiska Underskning.
[18] Wu Z C, Liu T Y. 2008. Magnetic gradient tensor: its properties and uses in geophysics[J]. Geological Science and Technology Information (in Chinese), 27(3): 107-110.
[19] Wynn W N, Frahm C P, Carroll P J, et al. 1975. Advanced superconducting gradiometer/magnetometer arrays and a novel signal processing technique[J]. IEEE Transactions on Magnetics, 11(2): 701-707.
[20] 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[J]. Chinese J. Geophys. (in Chinese), 58(7): 2556-2565, doi: 10.6038/cjg20150730.
[21] Zhang C D. 2006. Airborne tensor magnetic gradiometry——the latest progress of airborne magnetometric technology [J]. Chinese Journal of Engineering Geophysics (in Chinese), 3(5): 354-361.
[22] Zhang Z Y, Xiao C H, Yan H. 2009. Localization of a magnetic object based on magnetic gradient tensor at a single point [J]. Journal of Detection & Control (in Chinese), 31(4): 44-48.
[23] Zeng Z F, Wang Z J. 2003. SQUID and its application in geophysics[J]. Progress in Geophysics (in Chinese), 18(4): 608-613.
[24] 骆遥, 王明, 罗锋,等. 2011. 重磁场二维希尔伯特变换—直接解析信号解释方法[J]. 地球物理学报, 54(7): 1912-1920, doi: 10.3969/j.issn.0001-5733.2011.07.025.
[25] 马国庆, 李丽丽, 杜晓娟. 2012. 磁张量数据的边界识别和解释方法[J]. 石油地球物理勘探, 47(5): 815-821.
[26] 吴招才, 刘天佑. 2008. 磁力梯度张量测量及应用[J]. 地质科技情报, 27(3): 107-110.
[27] 袁圆, 黄大年, 余青露. 2015. 利用加强水平方向总水平导数对位场全张量数据进行边界识别[J]. 地球物理学报, 58(7): 2556-2565, doi: 10.6038/cjg20150730.
[28] 张昌达. 2006. 航空磁力梯度张量测量—航空磁测技术的最新进展[J]. 工程地球物理学报, 3(5): 354-361.
[29] 张朝阳, 肖昌汉, 阎辉. 2009. 磁性目标的单点磁梯度张量定位方法[J]. 探测与控制学报, 31(4): 44-48.
[30] 曾邵发, 王者江. 2003. SQUID及在地球物理中的应用[J]. 地球物理学进展, 18(4): 608-613.