1 引言

磁力勘探是所有地球物理方法中发展最早、应用最广泛的一种方法，在区域地质调查、矿产资源和油气勘探、工程地质和考古等方面都可不同程度地发挥作用.实测磁异常的反演解释则是磁力勘探解决各种实际地质问题的重要环节.常用的磁异常反演方法有物性反演法和欧拉反褶积法.

磁异常物性反演方法是基于反演理论、在最小二乘意义下使目标函数达到极小的线性或非线性反演.在引入足够约束条件的情况下，能够给出接近于实际地质情况的物性分布和几何形态.几十年来，国内外众多学者^[1~7]对磁异常反演方法和算法作了较多的研究与应用.由于这种反演方法在数学上呈病态，具有固有的多解性，需引入先验约束条件，计算维数巨大、计算量庞大，这些问题都制约着反演效果和实际应用.近年来，Portniaguine^[8]、Li^[9]、姚长利^{[10, 11]}、Pilkington^[12]等提出了一些改进方法和措施，在一定程度上改善了反演效果和效率.

磁异常欧拉反褶积方法^[13]以欧拉齐次方程为基础，运用磁异常、空间导数以及各种地质体具有特定的“构造指数”来确定异常场源的位置.它无需场源物性的先验知识，计算简便快速，适合于大面积磁测资料的分析解释.但是实际场源构造指数的未知性、反褶积解的发散性等实际问题制约着欧拉反褶积法的实用性.2007年Salem^[14]在常规欧拉反褶积法的基础上提出了磁倾角欧拉反褶积方法，它不需要已知场源的构造指数就能估计出场源的位置分布，进一步可估计出场源的构造指数.这将促进欧拉反褶积方法的发展和应用.

概率成像方法是Patella^[15]首次提出用于自然电位异常的解释，Mauriello^[16]将这种方法推广到了重力领域.Iuliano^[17]在任意标量和矢量地球物理场概率成像的基础上，又定义了联合概率成像的概念，并利用重力和磁力、重力和自然电位、磁力和自然电位对维苏威火山进行了联合三维概率成像研究，清晰地勾画出了火山岩浆通道系统.概率成像法无需先验信息约束，计算简便稳定，但它只对简单孤立异常成像效果较好，对复杂叠加异常成像效果差，分辨率低.2009年笔者^[18]在此基础上，提出了重力和重力梯度数据三维相关成像方法及基于平均值滤波分离异常的三维相关成像算法.由于重力梯度异常具有比重力异常本身高的分辨率，因此重力梯度数据三维相关成像的分辨率要比重力异常三维相关成像的分辨率明显提高.复杂模型数据试验验证了笔者的方法和算法能改善成像效果和提高成像分辨率.

本文将三维相关成像方法进一步推广到磁力勘探领域，提出磁异常ΔT的三维相关成像方法和基于熵滤波分离异常的三维相关成像算法，利用组合模型磁异常数据和实际磁测资料进行试验和验证.

2 方法原理 2.1 磁异常ΔT三维相关成像

在地形起伏的测区，笛卡儿坐标的（x，y）平面在基准面上，z取垂直向下为正，假设测区地下任意坐标为（x_q，y_q，z_q）、体积为v_q的第q个均匀磁化球体的磁化强度为J_q，磁矩M_q=J_qv_q，地磁场倾角为I₀、偏角为A'₀，磁化强度倾角为I、偏角为A'.由于均匀磁化球体外部磁场与位于球心的偶极子磁场等效，因此，这里把均匀磁化球体视为偶极子.则第q个偶极子在测区上任意点（x，y，z）处的总强度磁异常ΔT_q（x，y，z）^[19]可表示为

(1)

其中，μ₀为真空磁导率，

定义测区上实测ΔT异常与第q个偶极子测区上的ΔT_q异常的皮氏（Pearson）互相关系数为

(2)

其中，ΔT（x_i，y_i，z_i）为（x_i，y_i，z_i）观测点的实测ΔT异常，ΔT_q（x_i，y_i，z_i）为地下第q个偶极子在该观测点的ΔT异常，N_s为测区观测点总数.

假设J_q>0，则将公式（1）代入公式（2）可得

(3)

其中，B_q（x_i，y_i，z_i）为第q个偶极子的基函数，

(4)

根据柯西（Cauchy）不等式可知，

(5)

因此，公式（3）的相关系数C_q的数值范围为-1≤C_q≤+1.

从公式（3）可见，测区上实测ΔT异常与第q个偶极子在测区上的ΔT_q异常的归一化互相关实质上是实测ΔT异常与第q个偶极子的基函数的归一化互相关.相关系数C_q表征实测ΔT异常与偶极子基函数的相关性程度，其值的意义是实测ΔT异常是由第q个偶极子产生的可能性.C_q值为正，说明该偶极子磁性（磁化率或磁化强度）高；反之，C_q值为负，说明该偶极子磁性（磁化率或磁化强度）低. C_q值越接近于1，表明该偶极子磁性高的可能性越大；反之，C_q值越接近于-1，表明该偶极子磁性低的可能性越大.

将地下待成像的空间剖分成均匀网格，然后根据公式（3）从地下浅层到深层计算所有网格节点的相关系数，即实现对实测ΔT异常的三维相关成像，成像结果表征地下异常地质体的空间赋存状态和等效磁性分布.

2.2 基于熵滤波分离异常的三维相关成像算法

笔者在重力和重力梯度数据三维相关成像^[18]文中提出了利用基于平均值滤波分离异常的三维相关成像算法来提高成像分辨率.然而采用平均值滤波法分离异常，其剩余异常往往伴生着虚假异常，这样三维相关成像算法不仅对剩余异常中的正常异常计算成像出正常场源，同时也对虚假异常计算成像出虚假场源，影响了成像可信度.因此，异常分离方法的分离效果将一定程度上影响着基于异常分离的三维相关成像算法的成像效果.

当前位场异常分离方法众多，但大部分方法都很难从叠加异常中精细准确地分离出某个或某些地质体引起的异常，异常分离仍然是当今重磁资料处理解释中没有很好解决的一个难题.这里笔者研究Nikitin的地球物理异常分离的熵滤波法^[20]，进而提出基于熵滤波分离异常的三维相关成像算法来提高成像分辨率和改善成像效果.

熵滤波法是Nikitin^[20]在熵的理论基础上提出的一种空间域滑动窗口滤波方法.设任意待滤波点p的异常ΔT_p，熵滤波法所用的滑动窗口的中心为该点，窗口大小自定义，落入滑动窗口内的磁异常为ΔT_i（i=1，2，…，N），则待滤波点p的磁异常熵滤波分离的实现步骤如下.

①计算磁异常ΔT_i（i=1，2，…，N）的平均值ΔT_mean、最小值ΔT_min和最大值ΔT_max；

②计算异常ΔT_i（i=1，2，…，N）与ΔT_mean之差的绝对值的最大值，r_max=max{|ΔT_i-ΔT_mean |}；

③按照下式计算各点（i=1，2，…，N）的熵值，

(6)

④按照下式计算各点（i=1，2，…，N）的熵权系数，

(7)

⑤根据熵权系数h_i，按照下式对待滤波点p进行熵滤波，

(8)

⑥待滤波点p熵滤波后的剩余异常为

(9)

可见，当滑动窗口内某点磁异常的熵值越小，表明它的变异程度越大，在熵滤波中所起的作用越大，其熵权重值越大；反之，某点磁异常的熵值越大，表明它的变异程度越小，在熵滤波中所起的作用越小，其熵权重值越小.根据计算得到的滑动窗口内各点的熵权重值大小对待滤波点异常进行熵滤波，分离出剩余异常.

基于熵滤波分离异常的三维相关成像算法是将磁异常熵滤波分离和三维相关成像动态结合起来.在成像地下某一深度层时，首先根据深度大小选择熵滤波的滑动窗口大小，一般地，随着成像深度的增加滑动窗口大小逐步变大；然后对原始磁异常进行滑动窗口熵滤波，分离出剩余异常；最后对剩余异常进行该深度层的三维相关成像.按照此算法，从地表到深部逐层进行基于熵滤波分离异常的三维相关成像计算，成像的结果表征地下异常地质体的空间赋存状态和等效磁性分布.下面的组合模型数据和实际资料试验将验证这种算法的有效性.

3 组合模型数据计算

组合模型由处于三个不同深度层的、尺寸大小和磁化强度各异的11个垂直磁化直立长方体组成，各长方体几何参数和磁化强度见表 1所示.正演计算该模型在深度为0km处、测网为101×101、纵横向网格间距均为0.2km的平面网格上的磁异常ΔT.图 1a显示了各直立长方体在水平面的投影分布，其中形体最大的长方体A1、A2和A3处于最深层，长方体B1、B2、B3、B4和B5处于中间层，形体最小的长方体C1、C2和C3处于最浅层.图 1b显示了理论磁异常ΔT图，其中A1、A2、A3、B3、B4、B5和C3各长方体引起的异常呈孤立或相对孤立分布，而B1、B2、C1和C2各长方体引起的异常均叠加在A1和A2长方体引起的异常之上.

这里先对熵滤波法分离异常的效果进行检验.对组合模型磁异常ΔT分别进行平均值滤波法分离异常和熵滤波法分离异常，滤波窗口大小纵横向均为5个点距，分离出的剩余异常分别见图 2a和2b所示.图 2c显示了这两种方法在Y=8km处的剩余异常剖面对比，其中蓝色为平均值滤波法结果，红色为熵滤波法结果，黑色为理论磁异常曲线.可见，平均值滤波法的剩余异常中B层和C层各长方体的异常周围均伴生了明显的负虚假异常；而熵滤波法的剩余异常中的负虚假异常明显弱多了，尤其在B层各长方体的异常周围.因此，熵滤波法分离异常效果优于平均值滤波法，采用基于熵滤波分离异常的三维相关成像算法将减少或削弱虚假异常源的发生，改善成像效果.

下面对组合模型磁异常ΔT进行三种三维相关成像试验，第一种是常规的直接对ΔT进行三维相关成像，第二种是对ΔT进行基于平均值滤波分离异常的三维相关成像，第三种是对ΔT进行基于熵滤波分离异常的三维相关成像，这三种成像的深度范围均为0~3km，深度步长均为0.12km.图 3a和3b分别显示了磁异常ΔT在Y=8km处和在Y=14km处的剖面，图 3（c，d）分别显示了第一种常规三维相关成像结果在Y=8km处和Y=14km处的深度切片，图 3（e，f）分别显示了第二种基于平均值滤波分离异常的三维相关成像结果在Y=8km处和Y=14km处的深度切片.图 3（g，h）分别显示了第三种基于熵滤波分离异常的三维相关成像结果在Y=8km处和Y=14km处的深度切片.图中，相关值为正表明磁性高，正相关值越大表明磁性高的可能性越大，相关值为负表明磁性低，负相关值的绝对值越大表明磁性低的可能性越大.

从图 3（c，d）可见，常规三维相关成像对长方体A1、A2、A3、B3、B4和B5引起的异常成像效果较好，能成像出这些长方体的空间和等效磁性分布，其中B3、B4和B5的等效磁性的强度比模型实际偏弱些，但对长方体B1、B2、C1、C2和C3引起的异常成像效果差，或未能分辨出.这说明了常规三维相关成像对简单孤立异常成像效果较好，而对复杂叠加异常成像效果差，分辨率低.从图 3（e~h）可见，基于平均值滤波分离异常的三维相关成像和基于熵滤波分离异常的三维相关成像都能较好成像出模型A层、B层和C层各长方体的空间和等效磁性分布，横向和纵向分辨率都较常规三维相关成像有了明显提高，其中B1和B2的横向分辨率较好、纵向分辨率略低些.进一步对比图 3（e~h）可见，由于平均值滤波法分离异常伴生了明显的虚假异常，造成基于平均值滤波分离异常的三维相关成像在各长方体周围成像出了明显的磁性低的虚假异常源，而基于熵滤波分离异常的三维相关成像则削弱了虚假异常源的发生.因此，采用基于熵滤波分离异常的三维相关成像算法能有效提高成像分辨率、改善成像效果.

4 实际资料计算

图 4a是某地面磁测并化磁极处理后的化极磁异常ΔT等值线图，测区范围为6000m×5000m，平面数据网格为61×51，纵横向网格间距均为100m.对此化极磁异常数据进行基于熵滤波分离异常的三维相关成像试验，成像深度范围为0~2000m，深度步长为60m，共34个深度层.这里选取测区的AB剖面，其位置见图 4a的白色直线，绘制它的成像结果. 图 4b显示了AB剖面的化极磁异常ΔT曲线，图 4c显示了基于熵滤波分离异常的三维相关成像结果在此剖面的深度切片，其中，相关值为正表明磁性高，正相关值越大表明磁性高的可能性越大，相关值为负表明磁性低，负相关值的绝对值越大表明磁性低的可能性越大.从图可见，基于熵滤波分离异常的三维相关成像能成像出地下地质体的空间和等效磁性分布，具有良好的横向和纵向分辨率，其中，0~500m浅层磁性场源形体小，磁性高低相间分布，500m以下深层存在两处明显的磁性高的场源，形体较大，水平位置分别位于400~1200m和3900~4700m，深度均延深到1000 m以下.基于熵滤波分离异常的三维相关成像结果将为测区的地质解释提供参考.

5 结论

本文将重力和重力梯度数据三维相关成像方法推广到磁力勘探领域，推导并建立了磁异常ΔT的三维相关成像方法，同时提出了用基于熵滤波分离异常的三维相关成像算法来提高成像分辨率.本文方法无需常规物性反演中给出初始模型和正反演的迭代拟合过程，成像过程中不需要先验信息约束，方法简单，计算稳定.组合模型磁异常数据和实际化极磁异常资料试验分析表明，本文方法能有效成像出地下磁性体的空间赋存状态和等效磁性分布，具有良好的分辨率.