磁总场梯度数据三维反演方法及应用

孙石达,张青杉,陈 超,王浩然,陈海弟,戴继舒

(1. 中国冶金地质总局地球物理勘查院，河北保定 071000；2. 中国地质大学(武汉)地球物理与空间信息学院，地球内部多尺度成像湖北省重点实验室，湖北武汉 430074)



磁总场梯度数据三维反演方法及应用

孙石达1,2,张青杉1,2,陈 超2,王浩然2,陈海弟1,戴继舒1

(1. 中国冶金地质总局地球物理勘查院，河北保定 071000；2. 中国地质大学(武汉)地球物理与空间信息学院，地球内部多尺度成像湖北省重点实验室，湖北武汉 430074)

磁总场三方位梯度数据相对于总场数据包含更丰富的异常信息，将梯度数据应用于三维磁化率反演中，可以更准确地描述异常体。本文采用最小模型结构反演方法进行三维磁化率反演成像，并采用对数障碍法对磁化率的反演取值范围进行约束。通过模型试验，对磁总场异常数据及其三轴梯度数据进行单独反演、联合反演，结果表明三轴梯度联合反演结果可以更好地刻画异常体形态，更有效地分辨邻近异常体，反演物性更合乎实际。将此反演技术应用于大冶铁矿高精度航磁数据反演解释，取得了较好的应用效果。

磁总场梯度 磁化率 最小模型结构反演 对数障碍法

Sun Shi-da, Zhang Qing-shan, Chen Chao, Wang Hao-ran, Chen Hai-di, Dai Ji-shu. 3D inversion of gradients of total field magnetic anomalies and its application[J]. Geology and Exploration, 2015, 51(6): 1016-1024.

1 引言

磁法测量的对象除地磁场的总场强度外，还包括总场空间变化率(梯度)、总场三分量及其空间变化率(梯度张量)。其中，总场强度的梯度测量始于20世纪80年代，并已取得了较好的应用效果(Hood,1981；Hoodetal.,1989)。由于地磁场三分量测量受测量设备精度及姿态影响严重，目前尚未得到较为理想的测量结果，Schmidtetal.(1997,1998)证实由观测磁总场强度进而推算其三轴分量更加切实可行。磁梯度张量数据由于其包含了更为丰富的信息而被认为是磁法勘探的下一个突破点(Christensenetal.,2000)，Schmidtetal.(2006)介绍了梯度张量测量相比分量测量和总场梯度测量的优势，并讨论了梯度张量数据的应用前景。目前，国内航空磁测仍以观测磁总场强度为主，梯度测量方面垂直梯度测量开展相对较早(姚长利等，1996)，并且近几年也实现了航磁三方向梯度测量工作(李晓禄等，2009；郭华等，2013)。中国冶金地质总局物勘院于20世纪80年代开始使用核旋磁力仪开展直升机吊舱式垂向磁总场梯度测量工作，效果较为理想；2007年后，该院使用光泵磁力仪开展三方位磁总场梯度测量工作，工作模式为直升机吊舱式及固定翼软补偿式，取得了较好成果。

针对国内磁法勘探现状，本文将磁总场梯度数据作为主要研究对象。磁总场梯度测量较总场测量而言，其原始资料具备如下主要优点(Hoodetal.,1989；秦葆瑚，1985；张昌达，2006；Schmidtetal.,2006；张青杉等，2010)：

(1) 数据质量较高，各探头受外界磁干扰比较一致，且梯度观测仪器的精度要求高于总场观测，一定程度上提高了磁异常图的质量；

(2) 数据类别及数量有较大幅度提升，而且梯度数据包含垂直测线方向上一定距离的信息，更利于欠采样观测数据的解释，提高了磁测效能，同时不同梯度数据间可相互印证，更多的数据处理方法更利于推断解释工作；

(3) 梯度测量可以消除区域磁场和随时间变化磁场(地磁日变及磁暴)的影响，因此可以不设置日变站及进行日变校正。

就数据解译而言，总场梯度数据尤其水平梯度对于场源体的边界位置具有更高的分辨能力(张青杉，2009)，许多处理方法如欧拉反褶积法(Thompson,1982；王明等，2012)、总梯度模法(Roestetal.,1992；黄临平等，1998；张季生，2000)均需使用梯度信息推算磁性体边界位置或估算顶面边缘埋深等，而且总梯度模法还能够在一定程度上消除剩磁影响；石磊等(2012)对磁总场及其垂向梯度数据采用相关成像技术进行处理，结果表明梯度数据效果更好。

本文针对磁化率参数研究适用最小模型结构反演方法，并通过对数障碍法对磁化率范围施加约束。通过多组模型试验，对(单或多)梯度数据及总场数据的反演结果进行对比，可以看出两者均有准确可靠的表现。在磁性体形态描述的精细性、邻近磁性体的分离能力以及反演磁化率与其真值的偏离度等方面，梯度(尤其是多梯度)数据反演结果更为优越。将上述反演技术应用于大冶地区高精度航磁数据的反演解释，效果较为理想。

2 方法原理

2.1 磁总场及其梯度正演

磁总场及其梯度正演公式采用骆遥等(2007)提出的长方体磁场及其梯度无解析奇点表达式。假设地下空间剖分为M个规则网格排列的长方体，观测面上存在N个数据观测点，则第j个长方体单元在第i个观测点处的磁总场或其梯度异常均可简化为：

(1)

其中，mj为第j个单元的磁化率，Gii为根据正演公式计算得到的核系数，其只与正常场的磁化强度及方向、长方体单元的剩磁以及与观测点的相对位置有关。

由于位场的叠加性，整个模型空间对第i个观测点的异常可通过所有单元在该点产生异常值的线性叠加，因此可将所有数据点的正演计算写为矩阵形式：

T=Gm

(2)

其中G为N×M维矩阵，可称为核矩阵或敏感度矩阵，其中第i行第j列元素即为式(1)中核系数，T和m分别为正演结果和模型磁化率构成的N维向量和M维向量。

2.2 最小模型结构反演

磁化率反演的目的是根据磁测数据求取比较符合实际地质情况的地下空间磁化率分布。为使反演结果尽量精细，一般情况下地下空间剖分单元个数远多于观测数据点数，反演过程需要求解欠定方程组，同时由于观测数据中存在噪声和误差，仅仅通过精确拟合观测数据所获得的反演结果未必理想，因此需要加入某些确定的模型约束来控制解模型空间的结构，改善反演的病态程度，降低多解性，这种反演方法统称为最小模型结构反演，目前已被广泛应用于地球物理数据反演中(Lietal.,1998；Farquharson,2008；Sunetal.,2014)，而加入模型结构约束的过程也被称为正则化(Tikhonovaetal.,1977)。

最终的反演目标函数通常包括两部分：数据拟合差函数Φd和模型目标函数Φm，前者体现的是反演结果正演数据与观测数据的拟合程度，后者则表示反演结果的模型复杂程度和光滑程度，两者相对重要性由正则化参数β来确定。因此目标函数形式可写为：

(3)

目前常用的确定最优正则化参数的方法主要有L曲线法和广义交叉验证法(Asteretal.,2005)，并且均被较为广泛的应用到了重磁单数据反演中(Li and Oldenburg,2003；Farquharson and Oldenburg,2004)，但是由于后者目前较难应用于多数据联合反演的情况，因此本文中选用L曲线法来确定最优正则化参数，具体方法可参见Asteretal.(2005)，这里不再赘述。

本文以磁化率为反演参数，可将Φd和Φm表述为常规的计算矢量长度形式，式(3)演化为：

(4)

其中dobs为观测数据，G为依据所剖分地下空间及数据观测点位正演计算得到的核矩阵，m为磁化率模型空间，Gm即为所求取数据的正演表达，Wd为依据观测数据标准差确定的对角矩阵；Ws和Wi(i=x,y,z)分别是测算模型空间m模型复杂度的单位矩阵和三方向粗糙度的差分矩阵(Lietal.,1996)；mref为参考模型。Φg表示测度矢量长度的函数，不同的测度函数具有不同的性质(Farquharson,2008；杜劲松，2014；Sunetal.,2014)，可导致反演结果呈现不同的分布特征，本文采用常规的L2范数测度形式，式(4)推演如下：

(5)

这是应用较多的一种形式，得到的物性结果较为平滑，并且由于磁正演是线性问题，因此采用二阶范数求解时反演过程也是线性的，最终求解公式为：

(6)

上式适用于针对磁总场或磁梯度单一观测参数的三维反演，对于三方位磁总场梯度数据的联合反演问题，反演目标函数可表述为将针对三个参量的反演目标函数加权叠加(由于三个方位磁梯度正演核函数随深度衰减特性较为一致，本文采用相同的深度加权函数，以便于目标函数的实现)，即将式(6)中数据拟合差部分扩展为三梯度数据拟合差即可：

式中μl为各梯度数据反演权重，其余参数同式(6)。

(7)

2.3 物性约束

为避免出现有悖于地质理论的物性反演结果，提高反演结果的可靠性，合理描述地质体形态，对反演模型的物性范围进行合理约束是必要的。本文采用对数障碍法(Logarithmic barrier method)来实现物性约束。该方法最早用于解决不等值约束条件下的线性和二次规划问题(e.g.Gill et al.,1991；Rooset al.,1997)，是一种近似全局精确的范围约束方法(白富生等，2000)，也被广泛应用在多种重磁位场数据反演中并取得了较好的物性范围约束效果(e.g.LiandOldenburg,2003；Farquharson,2008；Martinezet al.,2013)。假设反演磁化率值约束在(m-,m+)的范围内，按照对数障碍法把该不等值约束转换为式(9)所示的障碍函数并将其作为反演目标函数的一部分，因此式(3)扩展为：

(8)

其中λ为障碍参数，Xbf为障碍函数：

(9)

(10)

使用对数障碍法约束反演物性时，给定的初始模型值m(0)需全部处于物性约束范围(m-,m+)内，而初始障碍参数λ(0)可通过式(10)求得。由于加入障碍函数使得式(5)的线性反演问题非线性化，因此本文采用加权迭代最小二乘方法(Iteratively Reweighted Least Squares，IRLS)求解此非线性问题，求解思路为用第(n-1)次迭代的观测数据拟合差计算第n次迭代的物性模型修改量。求解物性模型修改量Δm的公式推导如下(以单参量反演为例，对于多参量联合反演，只需比照式(7)将数据拟合差部分扩展为多参量即可)。

设初始物性模型为m(0)；设第(n-1)次迭代的物性模型为m(n-1)，其正演结果为d(n-1)(等同于Gm(n-1))，则第n次迭代可视为用物性模型变化量Δm拟合dobs与d(n-1)的差值；第n次迭代得到的物性模型即为m(n)=m(n-1)+Δm；因此第n次迭代反演目标函数由(8)式演化为：

(11)

在仅考虑单一参量反演的情况下，求取第n次Δm的最终表达式由式(6)变为：

(12)

式中，X1和X2均为对角矩阵，分别由障碍函数对物性模型向量的一阶和二阶求导所得，λ(n)为第n次障碍参数，通过上式即可求得第n次Δm。

为使物性模型保持在约束范围内，第n次物性结果可表示为：

(13)

其中γ为处于(0,1)区间内近于1的常数，η为最大允许步长，用于将m(n)约束在(m-,m+)内，其值可由(14)式求得。

(14)

为使反演过程稳定收敛，应使障碍函数对反演目标函数的影响逐步降低，故将下一次迭代的障碍参数设定为λ(n+1)=[1-min(η,γ)]λ(n)，可见障碍参数随着迭代次数越来越小，而障碍函数对整个目标函数的作用也相应越来越弱，当达到规定的终止条件时(如反演目标函数变化量小于整体1%、障碍函数小于反演目标函数1%)即可终止迭代，求得最终物性模型结果。

3 模型试验

为验证上述反演方法的可靠性，测试磁总场及其梯度数据的反演效果，以下采用倾斜体模型及两组组合模型进行试验。通过对磁总场及其梯度数据进行单参量反演、多参量联合反演，对反演结果在异常体形态描述准确性、邻近异常体区分能力以及与真实物性值接近程度等方面的表现特征进行分析。

图1 倾斜体模型三维显示图及其正演磁异常图Fig.1 Perspective view of tilting model and its magnetic anomalies of forward modelinga-三维模型；b-磁总场异常ΔT；c-磁总场北向梯度Tx；d-磁总场 东向梯度Ty；e-磁总场垂向梯度Tza-3D model;b-total-field magnetic anomalies ΔT;c-Tx,d-Ty;e-Tzare the total magnetic gradient in x,y and z directions, respectively

3.1 倾斜体模型

该模型因其倾斜特性多被应用于磁物性反演试验中(Lietal.,1996；Boulangeretal.,2001)。将地下空间剖分为21×21×10=4410个网格单元，网格单元边长均为50 m，倾斜体模型三维立体显示如图1a，其倾斜角度为45°，倾向东，倾斜体磁化率为0.1 SI，背景空间磁化率为0 SI(本文中所有模型试验均假定磁性体无剩磁，自退磁效应忽略不计)；地磁场倾角45°，偏角10°，强度为57000 nT；正演计算得到磁总场异常ΔT及其梯度异常如图1所示。为计算反演结果正演后与原始模型正演数据的偏差，在此未对异常数据加入任何噪声。

对上述ΔT及其梯度异常数据进行磁化率反演，反演过程中仅加入物性范围约束，约束范围为0～0.1 SI，正则化参数通过L曲线方法确定，得到最终磁化率反演结果(图2)。可以看出，ΔT及其三轴梯度单独反演结果大体相当，Tx、Ty反演结果分别在x及y方向表现略好；而三轴梯度联合反演结果则表现较好，既继承了单梯度数据反演结果的优点，能够对异常体边界进行较准确的刻画，同时对深部描述也更好。从反演模型及原始模型分别正演后标准偏差统计结果(表1)看，ΔT反演结果不能很好地拟合梯度数据，单梯度反演结果也不能很好地拟合ΔT及其它梯度数据，而三轴梯度联合反演结果对所有数据均拟合较好，具备真实刻画实际异常体的能力。

图2 倾斜体模型正演数据反演结果切片对比图Fig.2 Comparison of slices for forward and inversion results recovered from the anomalies in Fig.1a-理论模型；b-ΔT反演结果；c-Tx反演结果；d-Ty反演结果；e-Tz 反演结果；f-三方向梯度联合反演结果Slices through (a)the model and results recovered from (b)ΔT,(c) Tx,(d) Ty,(e) Tz and (f) joint three gradients

3.2 组合模型1

为进一步测试上述反演方法的效果，采用如图3a的组合模型，使水平x、y方向均存在高-高磁与高-低磁组合异常体，以便考察三轴梯度反演结果对邻近异常体的区分分辨能力。将地下空间剖分为31×31×10= 9610个网格单元，网格单元边长均为50 m；1、2、3号异常体磁化率均为0.1 SI，4号异常体磁化率为0.05 SI，背景空间磁化率0 SI；假定地磁场倾角45°，偏角10°，强度为57000 nT；正演计算得到磁总场异常ΔT及其梯度异常(均加入异常幅值2%的高斯噪声)如图3b～e所示。

表1 倾斜体下反演模型与原始模型正演结果标准差Table 1 Standard deviations between inversion and forward modeling for a tilting body

图3 组合模型1三维显示图及其正演磁异常图Fig. 3 Same as Fig. 1 but on combined model 1a-三维模型；b-磁总场异常ΔT；c-磁总场北向梯度Tx；d- 磁总场东向梯度Ty；e-磁总场垂向梯度Tza-3D model;b-total-field magnetic anomaly ΔT;c-Tx;d-Ty;e-Tz are the total magnetic gradient in x,y and z directions,re- spectively

反演过程同3.1节。为便于查看横向分辨力，本文给出深度200 m处切片图(图4)。可以看出，ΔT反演结果(图4b)的异常体边缘刻画能力弱于梯度数据反演结果，当异常体靠近时，ΔT反演结果中两者之间的边界或空隙渐渐模糊，容易形成两者连为一体的假象；Tx、Ty单独反演结果(图4c～d)在各自对应方向上能够较为清晰地区分异常体；Tz单独反演(图4e)与三梯度联合反演(图4f)结果在切片深度具有较为一致的异常体分辨能力，能够清晰地分辨所有异常体。

图4 组合模型1正演数据反演结果切片对比图Fig. 4 Same as Fig. 2 but on combined model 1a-理论模型；b-ΔT反演结果；c-Tx反演结果；d-Ty反演结果；e-Tz 反演结果；f-三方向梯度联合反演结果Slices through (a)the model and results recovered from (b)ΔT,(c) Tx,(d) Ty,(e) Tz and (f) joint three gradients

3.3 组合模型2

为检测反演技术的综合能力，本文模拟某地实际地质情况建立模型(图5a)，地下空间剖分为31×31×10 = 9610个网格单元，网格单元边长均为50 m；1号异常体磁化率0.05 SI， 2号异常体磁化率0.1 SI， 3号异常体磁化率0.15 SI；地磁场倾角45°，偏角10°，强度为57000 nT；令异常观测面为0.1 m，正演计算得到ΔT及其梯度异常(均加入异常幅值2%的高斯噪声)如图5b～e所示。

图5 组合模型2三维显示图及其正演磁异常图Fig. 5 Same as Fig.3 but on combined model 2a-三维模型；b-磁总场异常ΔT；c-磁总场北向梯度Tx；d-磁 总场东向梯度Ty；e-磁总场垂向梯度Tza-3D model;b-total-field magnetic anomaly ΔT;c-Tx;d-Ty;e-Tz are the total magnetic gradient in x,y and z direction respec tively

在此仅分析比对ΔT和三梯度联合反演结果(图6b～c)，可以看出，后者能量更为集中，也对异常体边界描述更准确，反演磁化率更接近实际。

图6 组合模型2正演数据反演结果切片对比图Fig. 6 Same as Fig.4 but on combined model 2a-理论模型；b-ΔT反演结果；c-三方向梯度联合反演结果Slices through (a)the model and results recovered from (b)ΔT and (c) joint three gradients

通过上述理论模型试验，三梯度联合反演结果与ΔT反演结果相比较，能够更为准确地刻画异常体形态，横向上可更好地区分邻近异常体，反演磁化率更接近实际。

4 实例应用

大冶铁(铜)矿床位处铁山岩体南缘中段与三叠系下统大冶组大理岩接触带，顶底板围岩分别为中细粒含石英闪长斑岩及黑云母透辉石闪长岩。据岩矿石物性资料，矿区铁矿石实测磁化率平均值约0.8 SI，黑云母透辉石闪长岩和闪长斑岩磁化率稍低，大理岩磁化率最低(朱永刚等，2006；高宝龙，2010；陶德益等，2011)。矿区1∶1万航磁ΔT数据经匹配滤波剔除区域场后用于反演(图7a)，图7b～d为推算的三方位梯度异常图。

图7 大冶矿区航磁总场异常及换算三方位梯度异常图Fig.7 Total-field magnetic anomalies in the Daye Mine and its derived gradients a-磁总场ΔT局部异常图；b-磁总场北向梯度Tx；c-磁总场东向梯度Ty；d-磁总场垂向梯度Tza-total-field magnetic anomaly ΔT；b-Tx;c-Ty;d-Tz are the total mag-netic gradient in x,y and z direction respectively derived from ΔT

对上述数据进行ΔT反演及三梯度联合反演(不考虑剩磁影响)，将地下模型空间剖分为56×36×20个单元，每个单元边长均为100 m；矿区地磁场倾角45.6，偏角-3；磁化率约束范围0～0.8 SI；反演结果如图8。ΔT反演及三梯度联合反演结果中的高值异常自西北向东南分别与铁门坎、龙洞、象鼻山、狮子山四个矿区对应良好，并且对四个异常体都能较为清晰的刻画。相比较而言，后者结果聚焦更好，且反演磁化率值更接近实际，可以更好的描述异常体浅部以及深部的形态，能够更为清晰的标示出可能存在矿体的高磁性异常体的空间分布范围。结果中(图8b)磁化率高值区南西侧存在北西西-南东东走向的连续低值带，由于该区矿产主要赋存于燕山期含石英闪长岩与三叠系下统大冶群灰岩的接触带上，灰岩在接触带附近已变质为大理岩及白云质大理岩(高宝龙等，2010)，因此推测该低值带为无磁或低磁性的大理岩或灰岩，同时三维立体图(图8f)中显示高磁性体的走向与当地主要构造方向即北西西向相同，其倾向以及与接触带的接触面整体向南西方向倾斜，与薛清泼等(2006)结论相同。

图8 ΔT反演结果及三梯度联合反演结果水平切片图、断面图及三维显示图Fig. 8 Horizontal and vertical sections and perspective view of the inversion results of ΔT and joint inversion of three gradientsΔT反演结果：a-550 m深度切片图；c-AA’断面；e-三维图；三梯度联合反演结果：b-550 m深度切片图；d-AA'断面；f- 三维图ΔT inversion result: a-horizontal sections at 550m depth; c-vertical section AA’; e-3D result; three-dimensional joint inversion results: b-horizontal sections at 550m depth; d-ertical section AA'; e-vertical section AA'; f-3D result

5 结论

本文将基于对数障碍法进行物性范围约束的最小模型结构反演方法应用于磁总场及其梯度数据反演中，通过模型数据验证了方法的可靠性，并将该反演方法应用于大冶矿区航磁数据反演，反演结果与实际地质情况较为一致。总体而言，ΔT数据反演结果对异常体边界刻画稍显模糊，异常体横向区分能力偏弱，反演磁化率与实际偏差稍大；水平单梯度数据反演结果偏重于刻画异常体相应边界，具备较好的横向分辨能力；垂向梯度数据反演结果在深部表现较好；三轴梯度数据联合反演结果保持了单梯度反演结果的优点，能够更准确地描述异常体的整体形态。值得指出的是，梯度数据反演结果会出现反演物性值远大于实际值情况，需要更准确的物性约束以确保反演结果真实可靠。

本文中模型试验及实际资料处理均未考虑剩磁，由于总梯度模对磁性体总磁化方向不敏感，后续研究中将补充梯度数据与总梯度模数据的联合反演，使本文反演技术能够适用于强剩磁情况。

Aster R C,Borchers B,Thurber C H.2011.Parameter Estimation and Inverse Problems [M].The second edition.Elsevier: 89-109

Bai Fu-sheng,Zhang Liang-sheng.2000.Approximate global exact bnarrier method for logarithmic barrier function [J].OR Transactions,4(3): 13-18

Boulanger O,Chouteau M.2001.Constraints in 3D gravity inversion [J].Geophysical Prospecting,49: 265-280

Christensen A,Rajagopalan S.2000.The magnetic vector and gradient tensor in mineral and oil exploration [J].Preview,84: 77

Du Jin-song.2014.Study on processing,forward modeling and inversion algorithms of satellite magnetic anomaly data in spherical coordinate system [D].Wuhan: China University of Geosciences(Wuhan): 84-85(in Chinese)

Farquharson C G.2008.Constructing piecewise-constant models in multidimensional minimum-structure inversions [J].Geophysics,73(1): K1-K9

Farquharson C G,Oldenburg D W.2004.A comparison of automatic techniques for estimating the regularization parameter in non-Linear inverse problems [J].Geophysical Journal International,156(3):411-425

Gao Bao-long,Tao De-yi,Zhan Ying-lin,Fan Zhi-xiong,Zhou Kui.2010.Application of aero-surface and borehole magnetic exploration to the prospecting of exhausted mines in the Daye iron mine [J].Geology and Exploration,46(3): 0483-0490(in Chinese with English abstract)

Gill P E,Murray W,Poncele D B.1991.Solving reduced KKT systems in barrier methods for linear and quadratic programming[Z].Technical Report SOL 91-7,Stanford University,Stanford,CA: 1-23

Guan Zhi-ning,Hou Jun-sheng,Yao Chang-li.1996.Application of aeromagnetic gradient data in geological mapping and metallogenetic prognosis of gold deposits [J].Geoscience Journal of Graduate School,China University of Geosciences,10(02): 239-249(in Chinese with English abstract)

Guo Hua,Wang Ping,Xie Ru-kuan.2013.A study of geological interpretation with the tri-axial aeromagnetic gradients [J].Progress in Geophys,28(05): 2688-2692(in Chinese with English abstract)

Huang Lin-ping,Guan Zhi-ning.1998.The determination of magnetic causative boundaries using total gradient modules of magnetic anomalies [J].Journal of East China Geological Institute,21(2): 44-51(in Chinese with English abstract)

Hood PJ. 1980. Aeromagnetic gradiometry: A superior geological mapping tool for mineral exploration programs. In: Weinstock H, and Overton WC. (eds.), SQUID applications to geophysics. Society of Exploration Geophysicists: 72-77

Hood P J,Teskey D J.1989.Aeromagnetic gradiometer program of the Geological Survey of Canada [J].Geophysics,54(8): 1012-1022

Li Xiao-lu,Chang Shu-shuai.2009.Aeromagnetic gradient survey and elementary application in sandstone type uranium deposits prospecting [J].Uranium Geology,25(06): 355-360(in Chinese with English abstract)

Li Y G,Oldenburg D W.1996.3-D inversion of magnetic data[J].Geophysics,61(2): 394-408

Li Y G,Oldenburg D W.1998.3-D inversion of gravity data[J].Geophysics,63(1): 109-119

Li Y G,Oldenburg D W.2003.Fast inversion of large-scale magnetic data using wavelet transforms and a logarithmic barrier method [J].Geophysical Journal International,152(2): 251-265

Luo Yao,Yao Chang-li.2007.Theoretical study on cuboid magnetic field and its gradient expression without analytic singular point [J].OGP,42(6): 714-719(in Chinese with English abstract)

Martinez C,Li Y,Krahenbuhl R.2013.3D inversion of airborne gravity gradiometry data in mineral exploration: A case study in the Quadrilátero Ferrífero,Brazil [J].Geophysics,78(1): B1-B11

Roest W R,Verhoef J,Pilkington M.1992.Magnetic interpretation using the 3-D analytic signal [J].Geophysics,57(1):116-125

Roos C,Terlaky T,Vial J.1997.Theory and algorithms for linear optimization: An interior point approach [M].Wiley Chichester:109-144

Schmidt P W,Clark D A.1997.Directions of magnetization and vector anomalies derived from total field surveys [J].Preview,70: 30-32

Schmidt P W,Clark D A.1998.The calculation of magnetic components and moments from TMI: A case study from the Tuckers igneous complex,Queensland [J].Exploration Geophysics,29: 609-614

Schmidt P W,Clark D A.2006.The magnetic gradient tensor: Its properties and uses in source characterization [J].The Leading Edge,25(1): 75

Shi Lei,Guo Liang-hui,Meng Xiao-hong.2012.3D correlation imaging of the vertical gradient of magnetic total field anomaly [J].Progress in Geophys,27(4): 1609-1614(in Chinese with English abstract)

Sun J J,Li Y G.2014.Adaptive Lp inversion for simultaneous recovery of both blocky and smooth features in a geophysical model [J].Geophysical Journal International,197: 882-899

Tao De-yi,Yang Hai-yan,Xiao Ming-yao,Gao Bao-long,Shu Xiu-feng.2011.The application effect of borehole magnetic exploration in the prospecting of the Daye iron mine [J].Resources Environment & Engineering,25(4): 358-363(in Chinese with English abstract)

Thompson D T.1982.EULDPH: A new technique for making computer-assisted depth estimated from magnetic data [J].Geophysics,47(1): 31-37

Tikhonov A V,Arsenin V Y.1977.Solution of Ill-Posed Problems[M].John Wiley and Sons:93-127

Wang Ming,Luo Yao,Luo Feng,Tian Song.2012.The application and development of Euler deconvolution in gravity and magnetic field [J].Geophysical & Geochemical Exploration,36(5): 834-841(in Chinese with English abstract)

Zhang Chang-da.2006.Airborne tensor magnetic gradiometry—The latest progress of airborne magnetometric technology [J].Chinese Journal of Engineering Geophysics,3(5): 354-361(in Chinese with English abstract)

Zhang Ji-sheng.2000.The application of 3D analytical signal technique to the treatment of aeromagnetic anomalies in south China [J].Geophysical & Geochemical Exploration,24(03): 190-196(in Chinese with English abstract)

Zhang Qing-shan,Ma Feng-lin,Xu Li-yun.2010.Study on the 3D gradient aeromagnetic survey [J].Geology and Exploration,46(6): 1087-1091(in Chinese with English abstract)

Zhu Yong-gang,Yu Chang-chun.2007.Analysis of magnetic properties of Daye iron deposit in Hubei province [J].Contributions to Geology and Mineral Resources Research,21(sup): 155-159(in Chinese with English abstract)

[附中文参考文献]

白富生,张连生.2000.近似全局精确障碍方法：对数障碍函数法[J].运筹学学报,4(3): 13-18

杜劲松.2014.基于球坐标系的卫星磁异常数据处理与正反演方法研究[D].武汉：中国地质大学(武汉)：84-85

高宝龙,陶德益,詹应林,范志雄,周 逵,2010.大冶铁矿接替资源勘查项目中“空、地、井”磁法测量的应用[J].地质与勘探,46(03):483-490

管志宁,侯俊胜,姚长利.1996.航磁梯度资料在金矿地质填图和成矿预测中的应用[J].现代地质,10(2): 239-249

郭 华,王 平,谢汝宽.2013.航磁全轴梯度数据地质解释优势研究[J].地球物理学进展,28(05): 2688-2692

黄临平,管志宁.1998.利用磁异常总梯度模确定磁源边界位置[J].华东地质学院学报,21(2): 44-51

李晓禄,常树帅.2009.航磁梯度测量及其在砂岩型铀矿勘查中的应用初探[J].铀矿地质,25(06): 355-360

骆 遥,姚长利.2007.长方体磁场及其梯度无解析奇点表达式理论研究[J].石油地球物理勘探,42(6): 714-719

秦葆瑚.1985.国外磁场梯度测量的新进展.国外地质勘探技术[J].(5): 18-24

石 磊,郭良辉,孟小红.2012.磁总场异常垂直梯度三维相关成像[J].地球物理学进展,27(4): 1609-1614

陶德益,杨海燕,肖明尧,高宝龙,舒秀峰.2011.井中磁测在大冶铁矿深部勘查中的应用效果[J].资源环境与工程,25(4): 358-363

王 明,骆 遥,罗 锋,田 嵩.2012.欧拉反褶积在重磁位场中应用与发展[J].物探与化探,36(5): 834-841

薛清泼,徐九华,陈伟,石教波,谢玉玲,2006.大冶铁矿尖林山—狮子山矿段接触带形态及控矿规律分析[J].有色金属(矿山部分),58(05): 14-16

张昌达.2006.航空磁力梯度张量测量——航空磁测技术的最新进展[J].工程地球物理学报,3(5): 354-361

张季生.2000.用三维解析信号技术处理华南航磁异常[J].物探与化探,24(3): 190-196

张青杉.2009.发展三维磁梯度勘查技术，提升深部铁矿勘查能力[C].中国地质学会2009年学术年会论文摘要汇编: 594-596

张青杉.2010.三维磁梯度测量的高精度方案及其应用前景[C].中国冶金地质总局科技大会论文集：412-423

张青杉,麻丰林,许丽云.2010.航空三维磁梯度测量方案研究[J].地质与勘探,46(5): 1087-1091

朱永刚,于长春.2007.湖北省大冶铁矿区内矿石磁性特征分析[J].资源环境与工程,21(sup): 155-159

3D Inversion of Gradients of Total Field Magnetic Anomalies and its Application

SUN Shi-da2,ZHANG Qing-shan1,2,CHEN Chao2,WANG Hao-ran2,CHEN Hai-di1,DAI Ji-shu1

(1.GeophysicalExplorationAcademyofChinaMetallurgicalGeologyBureau,Baoding,Hebei071000;2.HubeiSubsurfaceMulti-scaleImagingKeyLaboratory,InstituteofGeophysicsandGeomatics,ChinaUniversityofGeosciences,Wuhan,Hubei430074)

The 3-direction gradients of total-field magnetic anomalies (ΔT) contain more information compared with ΔT.Applying such data to 3D inversion of susceptibility allows us to describe anomaly bodies more accurately.This work employs the minimum-structure inversion method to conduct such inversion and uses the logarithmic barrier method to constrain the susceptibility limits.We compare the recovered models from inverting data ΔT,their gradients independently and jointly,separately.The results indicate that the inversion gradients,especially by the joint inversion,are capable of distinguishing nearby sources and can describe the sources more completely and exactly,and the inverted physical properties are closer to the true ones.And we apply this method to the aeromagnetic data acquired in Daye for inversion and interpretation,and obtain good results.

magnetic total-field gradients,susceptibility,minimum structure inversion,logarithmic barrier method

2015-03-28；

2015-07-09；[责任编辑]郝情情。

国家“973”计划课题(编号：2012CB416805)和科研院所基本科研资金资助项目(编号： WHS201211)联合资助。

孙石达(1990年-)，男，在读博士生，主要从事重磁数据处理、反演与解释方面研究。E-mail: shidasun.cug@gmail.com。

张青杉(1968年-)，男，中国冶金总局物勘院副总工程师，主要从事地球物理勘探与研究工作。E-mail: qingshan-zhang@163.com。

P631

A

0495-5331(2015)06-1016-9