磁张量梯度是磁场矢量在不同方向的导数，能更准确地描述场源体的几何形态(Jekeli, 1993; Bell et al., 1997; Schmidt and Clark, 2000)，提高了对磁源体的水平分辨率.随着超导量子干涉器件(Superconducting Quantum Interference Devices，SQUID)的发展，21世纪初，德国耶拿光子技术研究所(Institute of Photonic Technology，IPHT)率先实现了航磁全张量测量，为矿产资源勘探提供了一种新的技术手段(舒晴和周坚鑫，2006；Luo et al., 2015).我国航磁全张量测量技术研究起步于“十二五”，近年来，中国科学院上海微系统所和吉林大学先后研制出具有自主知识产权的航磁张量梯度测量系统，并在部分地区开展了试验飞行和应用示范.

获取场源体水平分布和深度信息是航磁张量梯度解释的主要目的，现有磁张量解释技术大多是从磁总场异常解释技术基础上发展而来(Gamey et al., 2004; Doll et al., 2006).Mikhailov等(2007)提出张量欧拉反褶积法进行全张量重磁数据的解释工作，但该方法需事先知道场源体的类型(构造指数).Schmidt和Clark (2006)讨论了磁张量各分量的性质及特征，以及张量不变量与地质体形状间的对应关系.Qruç(2010a)利用磁张量不变量的最大值计算规则地质体的深度，但针对实际数据难以开展；Qruç(2010b)利用磁张量数据的解析信号进行简单地质体的反演，但需给定构造指数；Beiki等(2011)利用特征向量法实现了航磁张量数据的解释，但该方法在应用中限定条件较多.马国庆等(2012)提出位场全张量数据解释的张量局部波数法，试验了张量数据不同组合方式的应用效果.

磁张量数据的解析信号具有降低磁化方向干扰的优点(Roest et al., 1992)，能准确圈定地质体的范围，但解析信号随着深度的增加而幅值降低，因此对较深地质体的应用效果不明显.本文提出的基于方向解析信号比值函数的磁张量数据均衡边界识别方法，可有效均衡不同深度地质体响应，从而同时突出不同深度地质体的边界，提高深部地质体的识别能力，对深部地球资源的开发与利用具有重要的意义.针对磁源体埋深，笔者提出了磁张量数据深度成像技术，利用实测张量数据与假定模型异常的相关系数确定地质体的深度，该方法不需进行复杂的反演运算，是大数据量解释的有效方法.通过理论模型试验和实际数据验证了所提出方法的准确性和有效性，并应用于河北迁安磁张量数据的解释，取得了较好的应用效果.

磁张量梯度是磁场矢量在x, y, z方向的变化率，可以写成如下矩阵形式：

(1)

针对矩阵M定义每一行磁梯度数据的模为方向解析信号，称作x, y, z方向解析信号，其表达式为：

(2)

(3)

(4)

为了解磁张量数据方向解析信号的平面特征，并为后续方法试验提供基础数据，选用棱柱体模型开展了正演模拟.图 1为水平位置分别为(30 m，30 m)和(70 m, 70 m)的两个埋深分别为15 m和20 m棱柱体所产生的磁张量异常及其方向解析信号，地磁倾角为60°，地磁偏角为0°，地质体磁化方向与地磁方向一致，地质体水平位置在图 1a中用黑框线标出.

图 1

Fig. 1

图 1 棱柱体磁张量异常及其方向解析信号 Fig. 1 Magnetic tensor anomalies and directional analytic signal of a prism

图 1结果表明磁张量数据能更好地描述地质体的分布特征，但仍会受到倾斜磁化影响而发生变形.方向解析信号A_x，A_y，A_z受倾斜磁化的影响较小，异常形态规则，且A_x, A_y的极值能准确地识别浅部地质体的位置信息，但是对深部地质体的分辨能力较差，边界较模糊.为了更加清晰地获得较深地质体的边界信息，提出张量数据的均衡边界识别法，其为方向解析信号的比值函数，表达式为：

(5)

根据图 1中所示方向解析信号的平面特征可知该方法的极值对应于地质体边界.

磁张量数据的深度成像技术是利用实测数据与假定地质体异常的相关性来获取地质体的深度信息.为了计算方便，假定地质体的位置按照如下网格分布(图 2).

图 2

Fig. 2

图 2 地质体分布方式 Fig. 2 Grid distribution of a geologic body

任意地质体均可以通过多个点源组合而成，因此j点处地质体以点源方式存在，计算点源所产生的张量异常与实测数据的相关系数.当模型假定位置与真实地质体的位置相一致时相关系数取得最大值，从而判定场源体的分布.张量数据存在9个分量，其中只有5个分量是独立的，因此磁张量数据的深度成像法的表达式为：

(6)

其中，M_{r* *}代表实测磁张量异常，M_{s* *}代表假定模型计算得到的磁张量异常；cov[M_{r* *}, M_{s* *}]= , , N为窗口大小，磁张量数据的相关系数是无量纲的.张量数据深度成像法采用加和的方式来获得最终相关系数的平均值，通过多分量综合利用可有效降低解释结果的不确定性，且加和的方式可有效地避免单个分量误差大而降低反演精度.此外，该方法无需进行复杂的反演运算，有效地提高了多参数张量数据的反演效率，是大数据量反演的一种有效手段.点源磁源体的张量异常表达式为：

(7)

(8)

(9)

(10)

(11)

(12)

其中，r_x=x-x₀，r_y=y-y₀，z₀为中心埋深，p₀为地质体磁化强度，为测点距场源中心的距离，l=cosIcosD，m=cosIsinD，n=sinI.I, D分别为磁化倾角和偏角.

现有张量边界识别方法包括张量不变量法、分量组合法和总水平解析信号法.张量不变量法的表达式为：

(13)

磁张量不同分量组合法的表达式为：

(14)

总水平解析信号法的表达式为：

(15)

为了检验磁张量梯度数据均衡边界识别方法的应用效果，利用图 1所示张量异常数据开展了模型试验，不同方法的边界识别结果如图 3所示.

图 3

Fig. 3

图 3 棱柱体张量异常边界识别结果 (a)张量不变量法；(b)分量组合法；(c)总水平解析信号法；(d)方向解析信号的均衡边界识别法. Fig. 3 Edge detection results of prism tensor anomalies (a) Tensor invariant method; (b) Components combination method; (c) Total horizontal derivative of analytic signal method; (d) Balanced edge detection method of directional analytic signal.

图 3a和图 3b分别为磁张量数据张量不变量和分量组合边界识别方法计算结果，可以看出张量不变量和分量组合边界识别法所识别出的边界均受到倾斜磁化的影响存在一定程度的变形，且较深地质体分辨能力较低.图 3c为张量数据总水平方向解析信号的计算结果，可以看出该方法能很好地完成对浅部地质体边界的识别，且几乎不受磁化方向的干扰，但是该方法对于较深地质体边界的分辨率仍然较低.图 3d为张量数据的均衡边界识别结果，其可清晰地显示不同深度地质体的边界，且减少了倾斜磁化的干扰，与地质体实际边界吻合.

为了试验均衡边界识别方法的稳定性，在磁张量异常中加入高斯白噪声，数据信噪比为30 dB，各方法的边界识别结果如图 4所示.

图 4

Fig. 4

图 4 含噪棱柱体张量异常边界识别结果 (a)张量不变量法；(b)分量组合法；(c)总水平解析信号法；(d)方向解析信号的均衡边界识别法. Fig. 4 Edge detection results of noise-bearing prism tensor anomalies (a) Tensor invariant method; (b) Components combination method; (c) Total horizontal derivative of analytic signal method; (d) Balanced edge detection method of directional analytic signal.

从图 4a和图 4b结果中可以看出，张量不变量和分量组合边界识别法由于噪声的干扰降低了分辨率，已经难以识别出较深地质体.图 4c为张量数据总水平方向解析信号的计算结果，该方法在噪声干扰下依旧能清晰识别浅部地质体的边界，但不能很好地识别较深地质体的边界.图 4d为方向解析信号的均衡边界识别结果，该方法在存在噪声情况下仍然可以清晰地识别出不同深度地质体的边界，稳定性更强.

采用磁张量数据深度成像法来计算地质体的深度信息(图 5).图 5a为计算得到的三维相关系数结果，可以看出存在两个相关系数高值区，给定了地质体的位置和分布.为了清晰地获得地质体的深度，截取如图 5a所示黑线位置的剖面数据，箭头代表剖面线方向.图 5b和图 5c分别为剖面位置的磁异常和解析信号()结果. 图 5d为截取剖面所在位置的深度成像结果，极值采用“+”标识，通过极值点可以得到地质体埋深分别为15 m和20 m，与模型实际埋深一致，且所识别位置与解析信号的极值点具有良好的对应关系.

图 5

Fig. 5

图 5 棱柱体模型数据相关成像结果 (a)张量数据的相关成像结果；(b)剖面磁异常；(c)剖面解析信号；(d)剖面位置相关成像结果. Fig. 5 Results of correlation imaging on prism model (a) Correlation imaging result of tensor data; (b) Magnetic anomaly of profile; (c) Analytic signal of profile; (d) Correlation imaging result of profile.

为了进一步检验相关成像技术的实用性，开展了叠加模型试验.模型采用板状体，水平位置分别位于-30，-20，0，30 m处，埋深5 m，磁倾角为60°，磁偏角为0°，地质体顺层磁化，磁张量异常及深度成像结果如图 6所示.图 6a和图 6c分别为无噪声垂直方向磁张量M_zz和其方向解析信号A_z，地质体的位置在图 6a中用“+”标识，可以看出垂直方向解析信号能有效地降低倾斜磁化的干扰，且极值点对应于地质体的水平位置.

图 6

Fig. 6

图 6 叠加模型的深度成像结果 (a)无噪声垂直方向磁张量异常；(b)含噪的垂直方向磁张量异常；(c)无噪声磁张量异常的z方向解析信号；(d)含噪磁张量异常的z方向解析信号；(e)无噪声磁张量数据的深度相关结果；(f)含噪磁张量数据的深度相关结果. Fig. 6 Depth imaging results of superimosed model (a) Noise-free magnetic tensor anomaly Mzz in the z-direction; (b) Noise-corrupted magnetic tensor anomaly Mzz in the z-direction; (c) Analytic signal of noise-free magnetic tensor anomalies in the z-direction; (d)Analytic signal of noise-corrupted magnetic tensor anomalies in the z-direction; (e) Depth correlation imaging result of noise-free magnetic tensor data; (f) Depth correlation imaging result of noise-corrupted magnetic tensor data.

图 6e为无噪声张量数据的深度成像结果，极值点所标识的埋深分别为4.5 m, 4.6 m, 5 m和5 m.对于水平距离较小的地质体，由于异常之间的相互严重干扰，反演结果的精度较低，因此在实际数据解释时可根据实际情况进行场分离操作.异常中加入信号比为30 dB的高斯白噪声，图 6b和6d分别为含噪声(信号比30 dB)垂直方向磁张量异常和z方向解析信号.图 6f为含噪磁张量异常的深度成像结果，深度分别为4.6 m, 4.6 m, 5 m和5 m，可以看出噪声的加入并未降低方法的反演精度，多张量数据的联合反演可有效地降低噪声的干扰.

2014年7月中国完成首套航磁张量梯度测量系统，并应用于河北迁安地区铁矿潜力评价.测区范围39.8°N—40.2°N，118.4°E—118.7°E，当地磁倾角为58.1°，磁偏角为-7°，平均飞行高度100 m，采样点距20 m，磁场数据和磁梯度数据测量精度分别为1.6 nT和2.0 nT/m，实测磁张量数据如图 7所示.

图 7

Fig. 7

图 7 测区实测磁张量异常 Fig. 7 Measured magnetic gradient tensor anomalies

从图 7中可以看出测区北部矿体呈现北东条带状分布，测区南部大多为孤立的块状体，呈现北西走向.采用张量数据均衡边界识别技术计算得到的地质体的水平分布如图 8所示，白色虚线圈定出5个较为集中的铁矿分布区域，铁矿以条带状分布为主，呈现北东和北西走向.

图 8

Fig. 8

图 8 磁张量数据的边界识别结果 Fig. 8 Edge detection result of magnetic tensor data

图 9所示为5个铁矿集中区域欧拉反褶积和深度成像技术计算得到的矿体深度.图 9a为磁张量数据欧拉反褶积结果，底图为原始磁异常，可以看出矿体的深度范围在200~600 m之间.图 9b为Line1剖面的深度成像结果，图 9c为Line2剖面的深度成像结果.从图 8和图 9中可以看出1号和2号矿藏呈现为北东向展布，3号和4号主要表现为北西向延伸，5号矿藏为北东向展布，但形态接近于球体.具体位置和信息如表 1所示，深度相关成像结果与张量欧拉反褶积法一致.

图 9

Fig. 9

图 9 实测张量数据欧拉及相关成像结果 （a）张量欧拉反褶积结果；（b）Line1剖面位置相关成像结果；（c）Line2剖面位置相关成像结果. Fig. 9 Euler deconvolution and correlation imaging results of measured data (a) Tensor deconvolution results; (b) Correlation imaging result of the Line1 profile; (c) Correlation imaging result of the Line2 profile.

表 1 Table 1

表 1 反演结果统计表 Table 1 Statistics of inversion results 序号 分布特征 水平范围(km) 中心埋深(m)深度成像 深度范围(m)张量欧拉结果 走向 1 x(500~505) y(4440~4445) 426 m 435±91 m 北西 2 x(502~510) y(4435~4443) 207 m 241±174 m 北西 3 x(503~508) y(4426~4433) 435 m 435±57 m 北东 4 x(503~506) y(4425~4427) 267 m 260±143 m 北东 5 x(503~504) y(4421~4423.5) 623 m 609±132 m 北东

表 1 反演结果统计表 Table 1 Statistics of inversion results

深度相关成像结果与张量欧拉反褶积法相一致，因此本文方法是准确和实用的.

本文提出的磁张量数据均衡边界识别和深度成像技术综合利用磁张量数据实现了地质体水平位置和深度的反演.磁张量数据均衡边界识别利用方向解析信号的比值来均衡不同深度地质体的响应，从而同时突出了浅部和深部地质体的边界，模型试验表明，均衡边界识别技术可清晰和准确识别不同深度地质体边界，且能有效地压制倾斜磁化的干扰，提高了对较深地质体的分辨能力.张量数据深度成像技术综合运用多分量数据的相关系数来获得地质体的深度信息，避免了复杂的反演运算，提高了大数据量张量数据解释的效率，模型反演结果表明，该方法能准确地获得地质体的深度，且具有较强的抗噪性.应用上述张量数据解释方法完成河北地区实测磁张量数据解释，取得了较好的应用效果，验证了方法有效性和实用性，为地区矿产资源勘查提供了有力的支撑.