2019, Vol. 62
磁张量梯度为磁矢量在不同方向的变化率，对于磁源体形态特征具有更高的分辨率和精度，能更清晰区分叠加异常和细节特征(Jekeli, 1993; Bell et al., 1997; Schmidt and Clark, 2000)，在矿产资源勘探中的应用越来越广泛.20世纪90年代国外就开始采用航磁张量测量技术开展金属矿勘探方面的工作(Blakely, 1995; Smith and Bracken, 2004; Gamey et al., 2004; Doll et al., 2006; Schmidt and Clark, 2006)，能更加准确地圈定地质体的分布，取得了非常显著的实际应用效果，此外，磁张量梯度数据还能有效提高水平分辨率以及去除深层背景异常的干扰.2000年开始我国也陆续开展了航磁梯度测量技术方面的研究，并采用磁总场测量设备搭建航磁梯度测量系统，在部分地区也开展了航磁梯度数据的采集实验，但航磁张量梯度测量系统还处于研发阶段.
Schmidt和Clark(2000, 2006)系统地分析了磁张量梯度数据的分布特征，并分析了磁张量矩阵的特征值、特征向量与磁源体分布的对应关系.起初，磁张量数据的解释方法主要是将原始磁异常的解释方法进行推广(Gamey et al., 2004; Doll, 2006)，利用磁张量梯度数据的多参量特征相互约束来提高结果的精度.Zhang等(2000)和Mikhailov等(2007)将欧拉反褶积法进行扩展应用于重磁张量梯度数据的解释.Qruç(2010a, b)利用重磁张量矩阵的特征向量实现地质体深度的计算；Beiki等(2011)利用磁张量矩阵的特征向量在给定构造指数及倾斜磁化参数的条件下实现深度的反演.Zhdanov等(2012)实现了磁张量数据的磁化率分布的反演；马国庆等(2012)提出位场全张量数据解释的局部波数法，综合利用数据特征来完成地质体的反演.马国庆等(2014)利用磁张量数据的解析信号实现了磁张量数据的解释；Yin等(2016)利用频率域正则化方法通过磁总场来转换计算磁梯度张量数据，获得了更加准确的结果.Zhou等(2015)利用交叉梯度反演实现了剩磁情况下的磁张量数据的反演.Luo等(2015)和张光等(2016)讨论了磁张量梯度测量技术的补偿方法.Yin等(2016)和Ren等(2017)提出了复杂模型的重磁张量梯度异常的快速正演算法.磁张量梯度测量在地下管线或水下导弹的定位中也得到良好的应用(刘继昊等，2017).磁张量梯度数据会受到倾斜磁化的干扰，而使异常形态不规则.为了降低磁化方向不固定和剩磁的干扰，磁异常解释大多采用其解析信号来完成(Hsu et al., 1996, 1998; Ma and Du, 2012; Cooper and Gordon, 2014, 2015; Ma et al., 2017)，能有效地降低倾斜磁化的干扰(Roest et al., 1992).
本文提出基于磁张量梯度数据的解析信号比值的均衡边界识别和空间位置反演技术.均衡边界识别技术为张量数据不同方向解析信号的比值函数，其能有效地均衡不同深度地质体的响应，在计算过程中针对高阶垂直导数的计算采用Laplace方程来进行.空间位置反演是通过解析信号比值与地质体位置参数的线性方程来获得地质体的位置和深度，并利用解析信号比值函数的分布特征作为筛选条件来获得更加准确的结果.通过理论模型试验和实际数据验证方法的准确性和实用性.
1 方法基本原理磁张量梯度数据是磁场M不同方向分量的变化率，可写为如下的矩阵形式：
|
(1) |
为了体现磁张量数据的分布特征，设计相互叠加的磁异常模型.图 1为1个埋深为15 m(左上)和3个埋深为10 m棱柱体所产生的磁张量梯度异常，磁倾角为70°，黑框线表示地质体的水平位置.
|
图 1 棱柱体所产生的磁张量梯度异常 Fig. 1 Magnetic tensor gradient anomalies generated by rectangular |
图 1中可以看出磁张量梯度数据能很好地凸显地质体的边界、角点的特征，尤其是水平方向的分量能有效突出地质体的分布特征.由于倾斜磁化的干扰，水平方向张量梯度的极值与地质体位置并不是良好对应.
解析信号为磁梯度数据的模，具有降低倾斜磁化干扰的优点.磁张量矩阵M每一行数据的平方和分别称为x, y, z方向解析信号(Beiki, 2010)，可表示为
|
(2) |
|
(3) |
|
(4) |
x和y方向解析信号的极值对应于地质体的边界.Beiki(2010)利用水平方向解析信号模(ED)和水平方向解析信号垂直导数的模(ED2)进行磁张量梯度数据的边界识别.
|
(5) |
|
(6) |
根据两种方法的极值对应于地质体的边界.利用棱柱体模型所产生的磁张量异常来试验边界识别方法的应用效果(图 2).图 2a为原始磁异常，可以看出斜磁化的干扰使规则形体所产生的磁异常发生变形.图 2b为水平方向解析信号模的计算结果，可以看出该方法的识别结果受倾斜磁化影响小，但是水平分辨率较低，且对于较深地质体的识别效果较差.图 2c为水平方向解析信号垂直导数模的计算结果，结果水平分辨率较高，边界更加收敛，但较深地质体的边界不清晰.为了提高边界识别方法对深部地质体的分辨能力，提出利用方向解析信号的比值构建均衡边界识别滤波器(Balanced Analytic Signal, BAS)：
|
图 2 不同方法的边界识别结果 (a)棱柱体的原始磁异常；(b)水平方向解析信号模的边界识别结果；(c)水平方向解析信号垂直导数模的边界识别结果；(d)不同方向解析信号比值的边界识别结果；(e)解析信号导数比值的边界识别结果；(f)阶跃解析信号导数比值的边界识别结果. Fig. 2 Edge detection results of different methods (a) Original magnetic anomaly of rectangular; (b) Edge detection results of the amplitude of horizontal direction analytic signals; (c) Edge detection results of the amplitude of the vertical derivatives of horizontal direction analytic signal; (d) Edge detection results of the ratio of analytic signal in different directions; (e) Edge detection results of the ratio of derivatives of analytic signal; (f) Edge detection results of the ratio of step derivatives of analytic signal. |
|
(7) |
均衡边界识别滤波器通过比值的方式可有效地均衡不同深度地质体的响应，从而有效地获得较深地质体的边界.图 2d为式(7)计算结果，可以看出该方法能有效地均衡不同深度地质体的响应，有效地提高了对深部地质体的分辨率，但是边界相对发散，对于直接来划定地质体的边界存在一定的困难.为了提高边界识别方法的水平分辨率，提出基于解析信号导数的均衡边界识别滤波器(Gradient Balanced Analytic signal, GBA)：
|
(8) |
为了获得更加准确、收敛的地质体边界，依据不同阶导数特征点的分布特征(Ma et al., 2016)，建立阶跃解析信号导数比值的均衡边界识别滤波器(Step Balanced Analytic signal, SBA)：
|
(9) |
由于位场数据的方向解析信号、总水平导数等函数并不是调和函数，对于式(8)和(9)所涉及方向解析信号垂直导数的计算，不能直接采用位场现有公式计算，应先将其展开为位场的函数，如
|
(10) |
分别计算公式(10)中每项的结果，然后求和来获得最终的导数计算结果，对于位场高阶垂直导数的计算采用Laplace方程来完成(马国庆等，2014)，该方法能将垂直导数的计算转变为水平导数的计算，水平导数采用空间域方法来进行，从而有效降低噪声的干扰.图 2e和2f分别为式(8)和式(9)的边界识别结果，可以看出两种方法的极值能清晰地获得地质体的边界，而阶跃解析信号导数比值所识别边界更加收敛，且与实际地质体边界吻合更好.
采用含噪声异常来试验阶跃解析信号导数比值边界识别方法的稳定性，图 3a为在图 2a模型异常中加入信噪比为40的高斯白噪声.图 3b为水平方向解析信号模的识别结果，该结果受噪声干扰较小，但是边界结果的分辨率较低.图 3c为水平方向解析信号导数模的计算结果，可以看出该方法水平分辨率较高，但较深地质体的边界不清晰，且由于噪声的干扰，边界信息较模糊.
|
图 3 含噪声磁异常的不同方法边界识别结果 (a)原始磁异常加入信噪比为40高斯白噪声后磁异常；(b)水平方向解析信号模的边界识别结果；(c)水平方向解析信号垂直导数模的边界识别结果；(d)不同方向解析信号比值的边界识别结果；(e)解析信号导数比值的边界识别结果；(f)阶跃解析信号导数比值的边界识别结果. Fig. 3 Edge detection results of different methods for noise-contaminated magnetic anomaly (a) Magnetic anomaly contaminated by Gaussian white noise with the signal-to-noise ratio is 40; (b) Edge detection results of the amplitude of horizontal direction analytic signals; (c) Edge detection results of the amplitude of the vertical derivatives of horizontal direction analytic signal; (d) Edge detection results of the ratio of analytic signal in different directions; (e) Edge detection results of the ratio of derivatives of analytic signal; (f) Edge detection results of the ratio of step derivatives of analytic signal. |
图 3d为不同方向解析信号比值的边界识别结果，该方法提高了对深部地质体的分辨率，且受噪声的影响较小.图 3e和图 3f分别为解析信号导数比值和阶跃导数比值的边界识别结果，可以看出由于引入拉普拉斯方程来计算二阶垂直导数，并未明显地增大噪声的效应，且方法所识别结果分辨率更高.
磁张量数据的不同方向解析信号满足如下的欧拉线性方程(马国庆等，2014)：
|
(11) |
|
(12) |
|
(13) |
不同方向解析信号所体现的地质体特征不同，且会造成部分地质体位置的响应与背景异常相同，因此单一分量反演并不能得到准确的结果(马国庆等，2014)，利用方程(11)—(13)组合来获得构造指数N与位置参数(x0，y0，z0).上述线性方程组求解会造成结果的发散效应，需通过设定相应的筛选准则来获得准确的结果.首先是根据数值的大小来判断是否为异常区，主要是利用参与反演函数的平均值作为门限值，当计算点上的数值小于0.5倍平均值时该点不参与反演.解析信号的幅值随着埋藏深度的增加急剧衰减，采用该准则进行反演时无法识别较深地质体，大大降低了对深部目标体的分辨能力.
利用解析信号均衡比值函数建立位置参数的线性反演方程，从而提高对深部地质体的分辨率，式(11)×Ax加上式(12)×Ay可得
|
(14) |
式(13)
|
(15) |
解析信号比值函数BAS在x，y，z方向导数分别为：
|
(16) |
|
(17) |
|
(18) |
将式(16)—(18)代入式(15)可得
|
(19) |
利用解析信号比值函数所构建的线性方程可估算出场源体的位置坐标x0, y0和z0.解析信号比值本身具有均衡不同深度地质体响应的优势，因此在利用其作为计算筛选条件时不会去除深部地质体，具体准则如下：
(1) 计算点解析信号比值小于0.5倍平均值时，该点不参与反演；
(2) 当求解结果与窗口中心的距离大于窗口半径时该解无效，因为当窗口在场源附近时解的准确性较高；
(3) 利用一个较小窗口滤除一些孤立的解，以求解结果为中心，计算窗口为半径，该范围内无其他结果时，该结果无效；
(4) 经过以上的筛选后，解已经比较集中.以结果为中心，窗口为半径的范围内解的个数小于3时，该结果无效，因为当解的个数小于3时不被认定为有效异常.
2 模型试验采用上述棱柱体模型所产生的磁张量梯度异常来试验解析信号比值线性方程的应用效果.图 4a为地质体所产生的原始磁异常，由于倾斜磁化的存在使异常形态发生了一定的变形.首先利用解析信号比值的线性方程来计算地质体的分布，图 4b为反演得到的地质体水平位置的分布，可以看出反演得到的水平与地质体实际位置相吻合，且形态也相对比较规则，降低了倾斜磁化的影响，反演结果中产生多余的干扰异常，因此本文方法所设定的筛选准则可有效地去除异常叠加所产生的干扰.
|
图 4 磁异常解析信号比值函数反演结果 (a)棱柱体原始磁异常；(b)式(19)反演得到的地质体位置信息；(c)反演结果的三维显示；(d)反演得到深度的统计. Fig. 4 Inversion results of analytic signal ratio function of magnetic anomaly (a) Original magnetic anomaly of rectangular; (b) Location of geological body inverted by the Eq.(19);(c) 3D view of inversion results; (d) Statistic of inverted depth. |
图 4c为反演结果的三维显示，由于地质体之间的相互干扰，反演深度存在小幅度的误差.图 4d为反演深度的统计，反演深度的平均值分别为9.7 m和14.6 m，与地质体实际埋深相吻合，因此本文方法能准确地计算得到地质体的位置信息.
在实际数据解释中噪声是不可避免的干扰因素.图 5a为在图 3a所示异常中加入信噪比为30的高斯噪声后的含噪磁异常.图 5b为噪声异常的阶跃解析信号导数比值的边界识别结果，可以看出该方法所识别边界与地质体实际水平位置相一致，且并未由于导数的计算而严重增大噪声的干扰.
|
图 5 含噪声磁异常的解析信号比值函数反演结果 (a)含噪磁异常；(b)含噪磁异常的阶跃解析信号导数比值的边界识别结果；(c)式(19)反演结果的水平分布；(d)反演结果的三维显示. Fig. 5 Inversion results of analytic signal ratio function of noise-contaminated magnetic anomaly (a) Noise-contaminated magnetic anomaly; (b) Edge detection results of the ratio of step derivatives of analytic signal of noise-contaminated magnetic anomaly; (c) Horizontal distribution of inversion results computed by the Eq.(19); (d) 3D view of inversion results. |
图 5c为利用式(19)反演结果的水平位置分布，可以看出噪声的加入降低了反演结果的精度，产生了额外干扰异常，但是反演得到的位置与地质体实际位置依旧吻合良好.图 5d为反演结果的三维显示，可以看出反演得到的深度范围与真实值相接近，并未由于噪声的加入而出现较大的偏差，因此方法具有较高的稳定性.
3 实际数据应用2013年7月我国首次引入国外航磁张量梯度测量系统，并在河北迁安铁矿区开展了磁张量梯度数据的采集，平均飞行高度为100 m，采样点距为20 m，总场数据测量精度为4.6 nT，经调平后实测磁张量梯度数据如图 6所示.
|
图 6 测区实测磁张量异常 Fig. 6 Measured magnetic tensor anomaly of studied area |
从图 6中可以看出测区北部异常均呈现北东向条带状分布，且连续性较好，因此具有较好的资源潜力，异常在局部存在圈闭，表明存在矿体的局部上涌.测区南部异常大多表现为圈闭形态，呈现北西向分布，因此大多为孤立的块状矿体.采用张量梯度数据的解析信号比值的均衡边界识别和空间位置反演技术获得地质体的分布(图 7).
|
图 7 本文方法反演得到的位置和深度 (a)磁张量梯度数据的阶跃解析信号导数比值的边界识别结果；(b)式(19)反演得到的地质体位置和深度. Fig. 7 Inverted location and depth by the method in this paper (a) Edge detection results of the ratio of step derivatives of analytic signal of magnetic tensor gradient data; (b) Horizontal location and depth computed by the Eq.(19). |
根据图 7a的边界识别结果可以圈定出5个较集中的铁矿分布区域(白虚线)，矿体以条带状分布为主，呈现为北东和北西走向.图 7b为地质体位置和场源的深度，反演得到的位置范围与边界识别方法相一致，深度结果表明5个铁矿区的深度集中在200~600m范围.矿体位置和深度具体信息如表 1所示.
|
表 1 矿体反演结果统计表 Table 1 Statistic inversion results of mineral source |
在图 7所示5号矿体的钻孔见矿深度为414 m，矿体底深为603 m，与实际反演得到的矿体深度范围相一致，因此本文方法结果具有良好的可信度.
4 结论基于磁张量梯度数据的分布特征提出解析信号比值的均衡边界识别及空间位置反演技术.解析信号比值的均衡识别技术为不同方向解析信号的不同阶导数的比值，其能有效地均衡不同深度地质体的响应，从而可同时显示不同深度地质体的水平范围，所识别的边界更加清晰和准确，且降低了倾斜磁化的干扰.解析信号比值的空间位置反演技术是通过解析信号比值与位置参数的线性方程来获得地质体的位置和深度，且根据数据特征设定了相应的筛选准则来提高计算效率和精度.将本文方法应用于河北迁安矿区实测磁张量梯度数据的解释，圈定出5个较大型矿藏，埋深集中在200~600 m，呈现为条带状分布，为矿产资源勘查提供了良好的基础资料.
Beiki M. 2010. Analytic signals of gravity gradient tensor and their application to estimate source location. Geophysics, 75(6): I59-I74. DOI:10.1190/1.3493639 |
Beiki M, Pedersen L B, Nazi H. 2011. Interpretation of aeromagnetic data using eigenvector analysis of pseudogravity gradient tensor. Geophysics, 76(3): L1-L10. DOI:10.1190/1.3555343 |
Bell R E, Anderson R N, Pratson L F. 1997. Gravity gradiometry resurfaces. The Leading Edge, 16(1): 55-60. DOI:10.1190/1.1437431 |
Blakely R J. 1995. Potential Theory in Gravity and Magnetic Applications. Cambridge: Cambridge University Press.
|
Cooper, Gordon R J. 2014. The automatic determination of the location and depth of contacts and dykes from aeromagnetic data. Pure and Applied Geophysics, 171(9): 2417-2423. DOI:10.1007/s00024-014-0789-8 |
Cooper, Gordon R J. 2015. Using the analytic signal amplitude to determine the location and depth of thin dikes from magnetic data. Geophysics, 80(1): J1-J6. DOI:10.1190/geo2014-0061.1 |
Doll W E, Gamey T J, Beard L P, et al. 2006. Airborne vertical magnetic gradient for near-surface applications. The Leading Edge, 25(1): 50-53. DOI:10.1190/1.2164755 |
Gamey T J, Starr T, Doll W E, et al. 2004. Initial design and testing of a full-tensor airborne SQUID magnetometer for detection of unexploded ordnance. //2004 SEG Annual Meeting. Denver, Colorado: SEG, 798-801.
|
Hsu S K, Sibuet J C, Shyu C T. 1996. High-resolution detection of geologic boundaries from potential anomalies: An enhanced analytic signal technique. Geophysics, 61(2): 373-386. DOI:10.1190/1.1443966 |
Hsu S K, Coppens D, Shyu C T. 1998. Depth to magnetic source using the generalized analytic signal. Geophysics, 63(6): 1947-1957. DOI:10.1190/1.1444488 |
Jekeli C. 1993. A review of gravity gradiometer survey system data analyses. Geophysics, 58(4): 508-514. DOI:10.1190/1.1443433 |
Liu J H, Li X H, Zeng X N. 2017. Online magnetic target location method based on the magnetic gradient tensor of two points. Chinese J. Geophys (in Chinese), 60(10): 3995-4004. DOI:10.6038/cjg20171026 |
Luo Y, Wu M, Wang P, et al. 2015. Full magnetic gradient tensor from triaxial aeromagnetic gradient measurements: Calculation and application. Applied Geophysics, 12(3): 283-291. DOI:10.1007/s11770-015-0508-y |
Ma G Q, Du X J. 2012. An improved analytic signal technique for the depth and structural index from 2D magnetic anomaly data. Pure and Applied Geophysics, 169(12): 2193-2200. DOI:10.1007/s00024-012-0484-6 |
Ma G Q, Du X J, Li L L. 2012. Comparison of the tensor local wavenumber method with the conventional local wavenumber method for interpretation of total tensor data of potential fields. Chinese J. Geophys (in Chinese), 55(7): 2450-2461. DOI:10.6038/j.issn.0001-5733.2012.07.029 |
Ma G Q, Huang D N, Li L L, et al. 2014. A normalized local wavenumber method for interpretation of gravity and magnetic anomalies. Chinese J. Geophys (in Chinese), 57(4): 1300-1309. DOI:10.6038/cjg20140427 |
Ma G Q, Huang D N, Liu C. 2016. Step-edge detection filters for the interpretation of potential field data. Pure and Applied Geophysics, 173(3): 795-803. DOI:10.1007/s00024-015-1053-6 |
Ma G Q, Liu C, Xu J S, et al. 2017. Correlation imaging method based on local wavenumber for interpreting magnetic data. Journal of Applied Geophysics, 138: 17-22. DOI:10.1016/j.jappgeo.2017.01.003 |
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): I61-I69. DOI:10.1190/1.2749317 |
Qruç B. 2010a. Depth estimation of simple causative sources from gravity gradient tensor invariants and vertical component. Pure and Applied Geophysics, 167(10): 1259-1272. DOI:10.1007/s00024-009-0021-4 |
Qruç B. 2010b. Location and depth estimation of point-dipole and line of dipoles using analytic signals of the magnetic gradient tensor and magnitude of vector components. Journal of Applied Geophysics, 70(1): 27-37. |
Ren Z Y, Chen C J, Tang J T, et al. 2017. Closed-form formula of magnetic gradient tensor for a homogeneous polyhedral magnetic target: A tetrahedral grid example. Geophysics, 82(6): WB21-WB28. DOI:10.1190/geo2016-0470.1 |
Roest W R, Verhoef J, Pilkington M. 1992. Magnetic interpretation using the 3-D analytic signal. Geophysics, 57(1): 116-125. DOI:10.1190/1.1443174 |
Schmidt P W, Clark D A. 2000. Advantages of measuring the magnetic gradient tensor. Preview, 85: 26-30. |
Schmidt P W, Clark D A. 2006. The magnetic gradient tensor: Its properties and uses in source characterization. The Leading Edge, 25(1): 75-78. DOI:10.1190/1.2164759 |
Smith D V, Bracken R E. 2004. Field experiments with the tensor magnetic gradiometer system at Yuma Proving Ground, Arizona. // Proceedings of the 17th EEGS Symposium on the Application of Geophysics to Engineering and Environmental Problems. EEGS, 1675-1691.
|
Yin G, Zhang Y T, Mi S L, et al. 2016. Calculation of the magnetic gradient tensor from total magnetic anomaly field based on regularized method in frequency domain. Journal of Applied Geophysics, 134: 44-54. DOI:10.1016/j.jappgeo.2016.08.010 |
Zhang C Y, Mushayandevu M F, Reid A B, et al. 2000. Euler deconvolution of gravity tensor gradient data. Geophysics, 65(2): 512-520. DOI:10.1190/1.1444745 |
Zhang G, Zhang Y T, Yin G, et al. 2016. Magnetic tensor compensation method for the carrier of the magnetic tensor detection system. Chinese J. Geophys (in Chinese), 59(1): 311-317. DOI:10.6038/cjg20160126 |
Zhdanov M S, Cai H Z, Wilson G A. 2012. 3D inversion of SQUID magnetic tensor data. J. Geol. Geosci, 1: 104. |
Zhou J J, Meng X H, Guo L H, et al. 2015. Three-dimensional cross-gradient joint inversion of gravity and normalized magnetic source strength data in the presence of remanent magnetization. Journal of Applied Geophysics, 119: 51-60. DOI:10.1016/j.jappgeo.2015.05.001 |
刘继昊, 李夕海, 曾小牛. 2017. 基于两点磁梯度张量的磁性目标在线定位方法. 地球物理学报, 60(10): 3995-4004. DOI:10.6038/cjg20171026 |
马国庆, 杜晓娟, 李丽丽. 2012. 解释位场全张量数据的张量局部波数法及其与常规局部波数法的比较. 地球物理学报, 55(7): 2450-2461. DOI:10.6038/j.issn.0001-5733.2012.07.029 |
马国庆, 黄大年, 李丽丽, 等. 2014. 重磁异常解释的归一化局部波数法. 地球物理学报, 57(4): 1300-1309. DOI:10.6038/cjg20140427 |
张光, 张英堂, 尹刚, 等. 2016. 一种磁张量探测系统载体的磁张量补偿方法. 地球物理学报, 59(1): 311-317. DOI:10.6038/cjg20160126 |