球坐标系下磁化强度矢量反演方法及其应用

杨 力， 周可法， 王金林

(1.中国科学院 新疆生态与地理研究所，新疆矿产资源研究中心，乌鲁木齐 830011；2.中国科学院大学，北京 100049)

0 引言

剩磁改变了总磁化强度的大小和方向，导致磁异常的幅值和形态产生畸变，从而影响磁测资料的反演解释。磁化强度矢量反演方法因其能同时反演得到岩石的磁化强度的大小和方向，消除岩石剩磁的影响，成为研究的热点。目前，研究磁异常强度矢量反演的方法主要有两种：①直接反演磁化强度矢量(如王妙月等[1-2]通过反演磁化强度矢量的三个分量最终获得磁化强度大小和方向，李泽林等[3-4]通过磁异常模量约束反演磁化强度矢量，Lelilèvre等[5]利用有限差分法、刘双[6]利用有限元方法重建地下高磁化率磁性体的三维磁化率分布并取得较好效果)；②利用地面化极异常垂直梯度和总梯度的互相关关系来估计磁化强度的倾角和偏角[7]。磁场强度矢量反演方法适用于笛卡尔直角坐标系和球坐标系[5]，大多数的反演方法都是基于笛卡尔直角坐标系展开的，但是磁化强度矢量反演将会产生三倍于常规磁化率反演方法的变量。另外，磁化强度矢量反演方法基于正则化方程计算磁化强度和磁化方向的三维分布，反演结果较为平滑，地质体界限模糊，将反演结果直接用于地质解释会导致很多信息丢失。在球坐标系中，地质信息的稀疏化约束和正则化方程具有较好的耦合关系，反演结果能够更加真实的反映地质体分布信息。

近年来，我国加紧对老矿山深部及外围的矿产资源进行勘查，促进了高精度航磁测量仪器和观测技术的发展。航空高精度磁测具有探测精度高、受地形影响小、快速高效的特点[8]，因此在球坐标系下对高精度磁测数据的反演解释是必要的。笔者首先推导由笛卡尔直角坐标系下磁化强度矢量反演的公式。然后将其转换到球坐标系下进行计算并建立试验模型，最后将该方法应用于对东天山卡拉塔格矿区，分析处理其1∶50 000高精度航磁测量数据[9]，对反演结果进行实测数据验证，且将反演结果与笛卡尔直角系下的反演结果进行对比，并进一步分析该矿床磁化强度异常的展布特征，综合地质和岩石磁性资料，探讨了红海矿床的找矿潜力。

1 磁化强度矢量反演方法

1.1 直角坐标系下的反演方法和原理

磁感应强度大小与观测距离之间的关系定义为式(1)。

(1)

将积分方程离散化后得到，b=TM，其中：T是二阶张量，M=[(MxMyMz)]T，b=[(bxbybz)]T。每一个离散单元的磁化强度为：

M=κH+MNRM

(2)

其中:κ是磁化率；H是磁场强度；MNRM是天然剩磁。磁感应强度观测值与m的关系可通过算子F定义为：

b=F(m)

(3)

其中:m是三维空间的磁化强度矢量模型。磁化强度矢量反演问题转化为求解m的问题，反问题可以写为：

m=F(m)-1dobs

(4)

其中:dobs是观测到的磁感应场强度大小，F(m)-1是算子F(m)的逆。根据Tikhonov[10]正则化方程，磁化强度矢量反演的反问题可以转化为求解正则化方程的形式：

φ(m)=φd+βφm

(5)

其中:Wd是标准差的倒数矩阵；W是空间距离加权矩阵，β是权重系数。对φ(m)求关于m的偏导数为“0”，可得：

βWTW(m-mref)=0

(6)

其中:雅克比矩阵J是预测观测数据对m的偏导数。易得式中数据拟合差项依赖于雅克比矩阵，通过空间距离矩阵W加入先验信息对模型范数项约束可以得到更加聚焦的解，模型范数项可定义为：

(7)

其中:Gx、Gy、Gz是三个方向上梯度的有限差分算子，[Ws,Wi,i=x,y,z]和[Rs,Ri,i=x,y,z]是稀疏化对角矩阵。敏感度矩阵Wr的选取能够对反演结果产生重要的影响。建立迭代更新的敏感度权重方程[11]，可以写为：

(8)

目前，航空磁测数据仍然以总磁场强度数据(TMI)为主，为了能够利用大量磁测数据，本文选取总场强度数据进行计算。总场强度可以定义为：

dTMI=F(κ)=PTH0κ=Fκ

dTMI∈Rn,F∈Rn×n,κ∈Rn

(9)

其中:H0是地磁场强度；P是将磁场转换为总场强度的投影矩阵。在迪卡尔直角坐标系中，磁感应强度大小可以写为：

b=T[(HpHsHt]κ

(10)

其中:Hp、Hs、Ht是笛卡尔右手直角坐标系中的磁场强度分块矩阵；p轴指向磁感应场方向；κ是等效磁化率；等效磁化率κe={κpκsκt} 可以定义为：

κp=MP/H0

κs=MP/H0

κt=MP/H0

(11)

式(9)中的总场强度可以写成：

dTMI=Fpκp+Fsκs+Ftκt

(12)

通过一个简单的立方体模型来说明直角坐标系下磁化强度矢量反演方法对结果的影响。立方体的体积为15 m×15 m×15 m，磁性体模型的的磁化率为0.05，磁倾角为-45°，磁偏角为215°。背景磁场的磁场强度设定为5 000 nT，并且磁倾角为90°，磁偏角为0°。构造的观测数据为总磁场强度数据，并且数据均带有1 nT的高斯白噪声用于模拟实际情况。具体模型的反演结果见图1。

图1 直角坐标系下的模型水平切片和模型垂直切片Fig.1 Plan section through the recovered magnetization vector model using the cartesian formulation(a)模型水平切片；(b)模型垂直切片

图2 球坐标系下的模型水平切片和模型垂直切片Fig.2 Plan section through the recovered magnetization vector model using the spherical formulation(a)模型水平切片；(b)模型垂直切片

1.2 球坐标系下的反演方法和原理

在球坐标系下，磁场强度大小{α}和磁化方向{θ,φ}能够较好地解耦。与直角坐标系下的反演相比具有两大优势，①野外磁测数据和实验室磁测数据都会给出球坐标系下的参数值，将其转换到直角坐标系会增加运算的繁琐程度；②球坐标系下能分别对磁化强度大小和磁化方向实施稀疏化约束。

对模型范数项求解关于{α,θ,φ}的偏微分，其形式可以写成：

(13)

J=FS

(14)

其中:矩阵S可以写成：

(15)

在球坐标系下，敏感度矩阵是非线性的[5]。在使用高斯-牛顿迭代方法进行计算时，模型范数项、磁化强度大小和磁化方向收敛速度很快。所以，敏感度矩阵也需要适应于这种变化。式(8)中的空间加权矩阵可以定义为：

(16)

由于磁倾角的取值范围为[-π/2,π/2]，磁偏角的取值范围为[-π,π]。磁化强度α的变化范围很大，通过选取最大值约束进行迭代更新。权重ωθ，ωφ可以定义为：

采用与直角坐标系下反演时同样的模型进行反演，图2(b)是球坐标系下磁场强度矢量反演结果。对比直角坐标系下磁化强度矢量反演得到模型的水平切片和垂直切片(图1)和球坐标系下的反演结果(图2)，可以看出与直角坐标系下的反演结果相比，在水平位置以及深度数值方面，球坐标系下的反演结果都能得到更加聚焦的解。

2 航磁异常在勘查区的应用

2.1 矿区地质特征

红海-黄土坡铜锌多金属矿床分布于卡拉塔格隆起区西南边缘的斜坡地带(图3)。红海铜锌多金属矿与黄土坡多金属矿分属不同公司的勘探区，但在地质上属于同一矿床。矿区发育的地层为晚奥陶-早志留统荒草坡群大柳沟组，为一套活动陆源性质的海相火山岩-火山碎屑岩建造，自下自上可分为四个岩性段。矿区主要出露第二及第四岩性段，第二岩性段岩性主要为英安岩、英安质凝灰岩、沉凝灰岩、砂砾岩，第四岩性段岩性主要为流纹岩、流纹斑岩、含砾流纹质沉凝灰岩、角砾凝灰熔岩。矿区中发育一套形成于早古生代活动大陆边缘弧构造环境的英云闪长岩-花岗闪长岩侵入组合，主要出露于矿区中南部，另有少量的北西-北北西向的石英闪长岩脉、二长花岗岩脉和流纹斑岩脉产出。地表发育北西向-北西西向-近东西向的断裂系统，它们都是继承早期基底构造发展而成的长期多次火山活动断裂，铜锌多金属矿就位于早期构造控制的火山构造盆地中[14]。

图3 卡拉塔格地区地质简图 [12] Fig.3 Geological map of the Kalatage region

图4 磁测数据ΔT等值线图Fig.4 ΔT anomaly contour map of aeromagnetic in Kalatage district

2.2 航磁ΔT异常特征

通过磁异常ΔT等值线图(图4)和地质图(图3)对比分析。其中，主要的航磁异常有M1、M2、M3、M4，其分布见图4。区内酸性岩(如流纹岩、花岗岩等)多具中低磁性，可引起低磁异常；中性岩(闪长岩)具有较高磁性，可引起较高的磁异常特征。研究区内沉积碎屑岩磁性较低，火山岩类随岩石基性程度有增强趋势，且磁性变化范围大(表1)。

航磁ΔT磁异常值变化范围大，磁异常主要沿北西西-近东西向和北西向呈带状分布，与区域构造方向一致。M1、M2、M4磁异常带比M3高，对应于研究区中的中性-基性岩石。M3低磁异常可能是被破碎带、碎屑岩覆盖的区域，由于ΔT值不能反映磁化率的三维分布结构，因此限制了对磁异常体的进一步分析。

表1 红海矿区岩矿石磁化率[13] Tab.1 Susceptibility parameter of rocks and minerals in Honghai deposit

图5 磁化率模型Fig.5 Magnetic susceptibility model through the recovered MVI model(a)直角坐标系下反演结果;(b)球坐标系下反演结果

3 反演结果

磁化强度矢量反演方法能够得到准确的磁化率三维分布。在反演中，将地下半空间剖分为30 m×30 m×10 m的单元，为了减少边界效应，模型的边界范围比研究区大，模型核部包含183×178×80个长方体单元。针对该区域对222 m地下深部进行切片分析。在第一节中通过模型试验验证了利用球坐标系下磁化强度矢量反演方法的优势。对红海-黄土坡矿区的实测数据进行反演，显示模型核部的反演结果。

直角坐标系下，反演得到的结果总是发散的，不能真实地反映地质体的展布特征(图5(a))。而在球坐标系下，能够得到聚焦的解(图5(b))。通过球坐标系下三维磁化率模型的切片图(图5(b))，可以划分出F1，F2，F3，F4四条断层，与图5a相比具有更清晰的边界特征。对比ΔT等值线图可以发现F2与F3之间的高磁异常区可能是磁性矿物的聚集区。红海-黄土坡铜多金属矿床属于隐伏矿床，通过磁化率反演得出深部F2、F3断裂控制了火山裂陷沉积中心的发育，并成为红海-黄土坡铜锌多金属矿矿质的热液通道和主要富集沉淀地。

图6 磁化方向的计算结果Fig.6 Magnetic inclination and declination through the recovered MVI model(a)磁偏角分布图; (b)磁倾角分布图

另外，通过磁化强度矢量反演得到的磁化方向分布，能够指示岩体和沉积物的形成时代，为矿体的圈定提供依据。从磁化方向计算结果中可以看到(图6(a),图6(b))，M3磁异常带具有强的偶极磁化异常，表明在F2与F3断裂中间沉积了大量的碎屑物质，而这些碎屑物质为VMS矿床的形成提供了良好的就位空间。一般火山碎屑沉积越厚，铜锌矿化越发育[15]。因此研究区深部有较大的找矿-成矿潜力，具有大型矿床规模远景。对比磁化方向和磁异常ΔT等值线图发现，M1和M2虽然在地表表现出较强的磁异常，但在深部表现出低的磁化率，表明M1和M2不具有成矿潜力。M1、M2、M3磁异常带，具有相同的磁化方向所以很有可能都是同一时期火山活动的产物，而M4磁异常区，也具有较强的偶极异常，同时在观测值中呈现出高的磁异常，也具有一定的找矿潜力。

4 结论

1)球坐标系下磁化强度矢量反演不仅比笛卡尔直角坐标系下反演的速度更快，而且得到的三维磁化率模型比笛卡尔直角坐标系下反演得到的结果更加聚焦。

2)矿区发育的M1、M2、M3磁异常带与区内断层发育方向一致，表明断裂控制了该区域磁异常分布。F2和F3之间的磁异常M3是可能的成矿远景区。