数据空间磁异常模量三维反演

李泽林, 姚长利, 郑元满, 孟小红, 张聿文

地下信息探测技术与仪器教育部重点实验室, 地质过程与矿产资源国家重点实验室，中国地质大学(北京)， 北京 100083



数据空间磁异常模量三维反演

李泽林, 姚长利*, 郑元满, 孟小红, 张聿文

地下信息探测技术与仪器教育部重点实验室, 地质过程与矿产资源国家重点实验室，中国地质大学(北京)， 北京 100083

强剩磁的存在通常导致了总磁化强度方向未知，进而影响了磁异常的反演和解释.磁异常模量是一种受磁化方向影响小的转换量，可以在强剩磁条件下通过反演三维磁化强度大小分布来推测场源分布状态.我们提出了一种数据空间磁异常模量反演算法来减少剩磁的影响.与标准的模型空间L2范数正则化反演方法相比，我们的方法有两个优点：一是无需搜索正则化参数(需要反复求解非线性反演问题)，因而可以减少计算时间；二是反演结果更加聚焦，深度分辨率更高，我们对此进行了原因分析.通过模型和实测数据测试证明了该算法的有效性和更好的反演效果.

数据空间； 磁异常模量； 剩磁； 三维反演

1 引言

三维磁化率反演(Li and Oldenburg，1996，2003；Pilkington，1997，2009；Portniaguine and Zhdanov，2002；姚长利等，2003，2007)逐步成为磁异常定量解释中的一种重要方法，已经成功地应用到了矿产和油气资源勘探以及区域地质构造解释中.传统的三维磁化率反演方法通常假设场源的磁化方向已知且与地磁场方向一致.然而，由于实际地质情况的复杂性，当强剩磁存在时，场源的磁化方向往往与地磁场方向大不相同，此时采用传统的反演方法会得到错误的结果.

近年来，针对上述问题，国内外学者已进行了相当深入的研究.在已有的研究成果中，减少反演中剩磁影响的方法主要有三种：第一种方法是首先估计磁异常的磁化方向(Lourenco and Morrison，1973；Phillips，2005；Dannemiller and Li，2006；Gerovska et al.，2009)，然后将估计的磁化方向用于三维反演，条件理想时，这种方法的反演精度较高，但不足之处是其仅适用于孤立场源情况.第二种方法是直接反演磁化强度矢量.王妙月等(2004)提出了二维层状模型的磁化强度矢量反演方法；Lelièvre和Oldenburg(2009)以及Ellis等(2012)提出了直接反演三维磁化强度矢量的方法；刘双等(2013)提出了基于磁异常模量的井中二维磁化强度矢量反演方法.这些方法可以同时反演磁化强度的大小和方向，但通常需要添加特定的约束来减少反演参数增多带来的更加明显的非唯一性(Lelièvre and Oldenburg，2009；Li S L and Li Y G，2014)，目前的应用还很少，研究还需进一步深化.第三种方法是首先将磁异常数据转化为受磁化方向影响小的转换量，然后再对这个转换量进行反演.Shearer(2005)和Li等(2010)通过反演磁异常模量来得到磁化强度大小的三维分布；Pilkington和Beiki(2013)通过反演归一化磁源强度来减少剩磁的影响.这种基于转换量的反演方法可以用于多场源复杂磁异常的反演.但是，进一步的分析也发现，由于转换量消除了原有磁异常中包含的相位信息，因此有时难以恢复地质体的准确产状，只能反演出大概位置(Li et al.，2010)，这方面的研究工作还需要不断深化.

本文主要的工作属于第三种方法中磁异常模量反演的范畴.因为比较而言，这种方法需要最少的先验信息，且可以用于多场源复杂磁异常的反演(Li S L and Li Y G，2014).磁异常模量反演已经成功应用到盆地火山岩勘探(Li et al.，2012b)、金属矿勘探(Santos et al.，2012；Ribeiro et al.，2013)以及海洋基底构造研究(Li et al.，2012a)，具有广泛的应用.由于磁异常模量与地下场源磁化率之间是非线性关系，Shearer(2005)和Li等(2010)提出的反演方法需要通过反复求解非线性反演问题来搜索正则化因子，这个过程非常耗时；此外，由于该算法属于L2范数正则化反演方法，因此反演结果较为模糊.针对上述两个问题，我们将数据空间反演算法(Pilkington，2009)扩展到了磁异常模量反演中，在数据空间中，无需搜索正则化因子，可大幅提高计算效率；此外由于模型范数定义为加权磁化率平方根向量的平方和，反演结果中的高磁化率部分和相邻单元较大的磁化率变化将会被保留(Lelièvre and Oldenburg，2006)，反演结果更加聚焦，可在一定程度上改善深部分辨率.

本文首先回顾磁异常模量和数据空间反演方法的基本原理，然后再推导基于数据空间的磁异常模量反演方法的迭代公式，最后通过模型试验和实测数据验证本方法的有效性和计算效率.

2 数据空间磁异常模量反演方法

2.1 磁异常模量

由于实际勘探所用的磁力仪测量的是磁性体产生的磁异常的某个投影分量，其结果受磁化方向影响很大.前人的研究发现磁异常模量具有弱敏感于磁化方向和与场源平面位置对应关系更好的特点，可以减少磁测数据解释中的剩磁影响(Stavrev and Gerovska，2000；Shearer，2005；Li et al.，2010；Liu et al.，2013；Li S L and Li Y G，2014).在三维情况下，磁异常模量的定义如下：

(1)

其中Ta是磁性体产生的磁异常矢量Ta的模量，Hax、Hay、Za分别为磁异常矢量Ta在x、y、z三个方向的投影分量.

由于现在磁力仪(如质子磁力仪)实测数据通常为总场异常ΔT(ΔT相当于是磁异常矢量Ta在地磁场方向上的投影分量)，因此需要利用ΔT换算磁异常模量.对于平面网格数据，可通过频率域转换的方法(Pedersen，1978)将总场异常转化为异常三分量，然后投影合成得到模量数据.对于起伏地形观测数据，常规分量转换不易实现，可以采用空间域的等效源方法(Dampney，1969)，进一步得到模量数据，然后再进行基于模量数据的反演.由于磁异常模量与磁异常的量级一致，因此在转换的过程中不会放大噪声，保持了长波长信息.

2.2 数据空间的优势

在磁异常模量反演中，假设模型和数据个数分别为M和N，由于M通常是远大于N的，特别是三维反演情况.所以，该反演属于纯欠定问题.在模型空间需要求解一个M×M的方程组，而在数据空间中只需要求解一个N×N的方程组(Tarantola，1987)，因此，数据空间反演算法将在一定程度上提高计算效率，减少计算时间.此外，由于模型空间中包含庞大的零空间(至少M-N维)，在模型空间中求解需要施加正则化来改善系数矩阵的病态程度，而正则化因子的搜索是一个非常耗时的过程(通常是需要反复求解非线性反演的问题)；在数据空间中则不存在零空间，因此大大提高了系数矩阵的稳定性，同时因为无需搜索正则化参数，所以，计算量得到进一步减少.为此，这里我们立足于数据空间进行磁异常模量反演(关于模型空间和数据空间的差异，在附录A-1中有介绍).

另外，我们经过深入分析认为，模型空间中改善系数矩阵的病态程度的正则化过程，也有其副作用，就是增加了反演结果图像的模糊性.而这恰恰是基于数据空间不需要的，因而，基于数据空间的反演具有更好的模型聚焦度，也即反演结果图像会更加清晰.

2.3 模型定义与磁化率正约束

在实际矿产资源勘探中，绝大多数岩石的磁化率均为正值，在反演中加入磁化率正约束可以有效地恢复地质体的产状信息，得到更有地质意义的反演结果(Li and Oldenburg，1996).首先定义有效磁化率κ=μ0M/T0，其中，μ0=4π×10-7H·m-1为真空中的磁导率，M为有效磁化强度，单位为A·m-1，T0为地磁场强度，单位为T.在数据空间反演算法中，利用平方根算子m=κ1/2对模型进行投影变换实现正约束(Lelièvre and Oldenburg，2006)是一种简单的技巧，其中m为模型向量，κ为有效磁化率向量.即模型向量可以在(-∞,∞)变化，但有效磁化率向量只能在(0,∞)之间变化.这样做的优势有两个：一是反演所需的雅可比矩阵与常规方式比较，具有一种更简单的形式；二是反演结果中的高磁化率部分和相邻单元较大的磁化率变化将会被保留，反演结果更加聚焦(Lelièvre and Oldenburg，2006)，这个优点也是我们需要重视的.

2.4 总目标函数

磁异常模量反演和磁异常分量反演的目标函数形式上是一样的，带参考模型约束的总目标函数有如下形式：

S(m)= (d-F(m))TD-1(d-F(m))

+(m-m0)TW-1(m-m0),

(2)

2.5 目标函数的优化

对(2)式的求解通常是在模型空间进行的.这里我们参考Tarantola(1987)和Pilkington(2009)的处理方式，建立起磁异常模量基于数据空间反演的目标函数求解公式.经过推导在数据空间中极小化(2)式可以得到模型m的不动点迭代方程(推导过程见附录A-1)：

m=m0+WJT(D+JWJT)-1(d-F(m)

+J(m-m0)),

(3)

其中J代表F(m)的雅可比矩阵(维数为N×M，具体形式见附录A-2).由于(3)式直接对模型进行迭代更新，其求解的稳定性和收敛性会很差，所以，我们采取将每次迭代得到的模型m作为下次迭代的参考模型，得到第k次迭代时数据空间的近似迭代公式：

Δmk=mk+1-mk

(4)

(4)式对模型的修正量进行迭代更新，稳定性和收敛性较好.其中α为迭代步长，实际计算中，初始选择1，计算拟合差，若不减小，则逐次减少为原来的1/3，直至拟合差减小为止.

将(4)式改写为紧凑形式：

(5)

fk=Akbk,

(6)

2.6 实际应用中需要注意的问题

在实际应用中，还需要注意以下几个方面的问题：

首先，在实测数据中，可能存在异常不完整或分布于数据体边缘的情况，需要在模型剖分时对水平方向进行适当的扩边，数值模拟表明，这样做不但可以减少边界效应，还可以加速算法的收敛.

其次，在反演迭代的过程中，建议将共轭梯度迭代的残差设置为期望的噪声水平(Pilkington，2009)，这种不完全共轭梯度迭代有两个优点：一是通过降低迭代的次数来减少反演时间，二是可以起到类似正则化的作用，提高反演算法的收敛性和稳定性.

另外，因为雅可比矩阵计算的要求，初始模型m0不能给0，通常可以给一个较小的值，例如0.001.

3 模型实验

我们采用两个简单模型和一个复杂模型来测试反演算法的有效性.简单模型包括立方体模型和倾斜板状体模型，复杂模型为组合倾斜板状体模型.这三种模型已被多位学者用于重磁异常物性反演的测试(Li and Oldenburg，1996，1998；Portniaguine and Zhdanov，2002；姚长利等，2007).

观测数据均为21×21的网格数据，网格间距均为50 m.我们首先利用最小曲率法对磁异常进行扩边，然后利用频率域方法将其转化为模量，最后分别采用模型空间反演方法(Li et al., 2010)和数据空间反演方法对模量进行反演.在反演中，模型剖分为立方体单元，边长与观测数据的网格间距一致，为了减少边界效应，模型剖分为31×31×10的立方体，同时为了便于对比，只显示中间的21×21×10的数据.在模型空间反演方法中，为了减少模型的光滑效应，我们将αx、αy、αz均设置为0，使模型保持最高的分辨率.在模型空间和数据空间反演方法中，初始模型均设置为0.001.在主频为2.90 GHz的计算机上，利用两种反演算法对上述三个模型进行了反演测试.

3.1 简单模型

模型一是一个边长为200 m的立方体，立方体的中心位于地下250 m，其有效磁化率为0.05SI，地磁场强度为50000 nT，地磁场倾角和偏角分别为90°和0°.磁化倾角和偏角分别为60°和30°.立方体模型产生的磁异常如图1a所示，其中包含了标准差为2%数据绝对值加上0.5 nT的高斯噪声.由磁异常换算得到的模量如图1b所示，可以看出，模量与立方体在水平面上的投影位置对应关系良好.

图2a、2b和2c、2d分别为采用模型空间反演算法和数据空间反演算法的反演结果切片图.可以看出，模型空间的反演结果较为光滑，场源边界模糊一些，且磁化率值整体偏低.而数据空间的反演结果聚焦程度高，场源边界清晰且磁化率值更接近模型的真实磁化率.数据空间反演共需185次共轭梯度(以下简称CG)迭代，耗时2.2 s.模型空间反演共需802次CG迭代，耗时14.3 s.

需要指出的是，由于三维磁性体磁异常模量依然受到磁化方向的影响，所以两种方法的反演结果都略微偏离了真实模型的位置.但是，这是模量反演中在磁化方向未知的情况下得到的结果.同样条件下，如果是磁异常的反演，则磁化方向的选择误差会造成很大的结果误差，甚至是完全错误的结果，这样的情况，我们在下面的组合模型例子中得到了更明显的体现.

图1 (a)立方体模型产生的磁异常；(b)由磁异常转换得到的模量Fig.1 (a) The total field anomaly produced by the cube model； (b) The amplitude data computed from the total field anomaly in Fig.1a

图2 图1b所示的磁异常模量的反演结果，黑框代表真实模型的位置(仅显示有效磁化率为0.0015SI以上的值) (a)和(b)利用模型空间反演算法的Z=-225 m和N=500 m的反演结果切片图； (c)和(d)利用数据空间反演算法的Z=-225 m和N=500 m的反演结果切片图.Fig.2 Inversion results of the amplitude data in Fig.1b. The black lines indicate the position of the true model. Only values above 0.0015SI are displayed(a，b) Model-space inversion results, slice at 225 m depth and 500 m north； (c，d) Data-space inversion results, slice at 225 m depth and 500 m north.

模型二是一个45°倾角的倾斜板状体，其有效磁化率为0.05SI，地磁场倾角和偏角分别为65°和-25°，磁化倾角和偏角分别为45°和75°.倾斜板状体模型产生的磁异常如图3a所示，其中包含了标准差为2%数据绝对值加上2 nT的高斯噪声.由磁异常换算得到的模量如图3b所示，可以看出，模量与板状体在水平面上的投影位置对应关系良好.

图4a、4b和4c、4d分别为采用模型空间反演算法和数据空间反演算法的反演结果切片图.可以看出，模型空间的反演结果场源边界模糊，产状不明显，但磁化率极大值接近真实场源磁化率.而数据空间的反演结果场源边界清晰，产状明显，不足之处是磁化率极大值大于模型的真实磁化率，这也是聚焦性体现的一个附带现象.

本例的数据空间反演共需89次CG迭代，耗时1.6 s.对应的模型空间反演共需832次CG迭代，耗时15.2 s.

3.2 复杂模型

模型三为组合倾斜板状体模型，由两个长宽、磁化率、延伸长度、磁化方向不同，倾向相反但走向相同的倾斜板状体组成，主要用于测试干扰场源的影响.地磁场倾角和偏角分别为65°和-25°.西侧板状体的有效磁化率为0.04SI，磁化倾角和偏角分别为-25°和-25°.东侧板状体的有效磁化率为0.05SI，磁化倾角和偏角分别为45°和75°.该组合模型产生的磁异常如图5a所示，其中包含了标准差为2%数据绝对值加上2 nT的高斯噪声.由磁异常换算得到的模量如图5b所示.

图3 (a) 倾斜板状体模型产生的磁异常；(b) 由磁异常转换得到的模量Fig.3 (a) The total field anomaly produced by the dipping slap model； (b) The amplitude data computed from the total field anomaly in Fig.3a

图4 图3b所示的磁异常模量的反演结果，黑框代表真实模型的位置(仅显示有效磁化率为0.0066SI以上的值)(a)和(b)利用模型空间反演算法的Z=-225 m和N=500 m的反演结果切片图；(c)和(d)利用数据空间反演算法的Z=-225 m和N=500 m的反演结果切片图.Fig.4 Inversion results of the amplitude data in Fig.3b. The black lines indicate the position of the true model. Only values above 0.0066SI are displayed(a，b) Model-space inversion results, slice at 225 m depth and 500 m north； (c，d) Data-space inversion results, slice at 225 m depth and 500 m north.

对于叠加异常，为了体现模量反演和常规磁异常反演的差异，我们首先假设磁化方向与地磁场方向一致，利用标准磁异常三维反演算法对图5a所示磁总场异常进行了反演.图6a、6b和6c、6d分别为采用模型空间反演算法(Li and Oldenburg，1996，2003)和数据空间反演算法(Pilkington，2009)的反演结果切片图.可以看出，标准磁异常三维反演算法的反演结果几乎是完全错误的.

然后我们利用磁异常模量反演算法对图5b所示磁异常模量进行了反演，图7a、7b和7c、7d分别为采用模型空间反演算法和数据空间反演算法的反演结果切片图.可以看出，磁异常模量反演算法的反演结果较好.比较而言，模型空间的反演结果难以区分两个板状体，而数据空间的反演结果场源边界则更加清晰.当然，两种反演方法均不能反演出此模型的深部构造特征，原因主要是深部构造产生的磁异常较小，且被西侧板状体产生的异常所掩盖.

图5 (a)组合倾斜板状体模型产生的磁异常；(b) 由磁异常转换得到的模量Fig.5 (a) The total field anomaly produced by the model composed of two dipping slaps；(b) The amplitude data computed from the total field anomaly in Fig.5a

图6 图5a所示的磁总场异常的反演结果，黑框代表真实模型的位置(仅显示磁化率为0.0046SI以上的值)(a)和(b)利用模型空间反演算法的Z=-125 m和N=500 m的反演结果切片图；(c)和(d)利用数据空间反演算法的Z=-125 m和N=500 m的反演结果切片图.Fig.6 Inversion results of the total field anomaly data in Fig.5a. The black lines indicate the position of the true model. Only values above 0.0046SI are displayed(a，b) Model-space inversion results, slice at 125 m depth and 500 m north；(c，d) Data-space inversion results, slice at 125 m depth and 500 m north.

反演计算对比中，数据空间反演共需189次CG迭代，耗时2.5 s.模型空间反演共需1435次CG迭代，耗时20.1 s.

从上述模型的反演结果来看，数据空间反演算法会得到一个更加聚焦的反演结果，同时计算量和计算时间更少.

4 实际资料试验

图8a为澳大利亚某地区的实测航磁数据，数据大小为53×71，网格距为50 m.该地区地磁场强度为58240 nT，地磁倾角和偏角分别为-64.7°和2.1°.

图7 图5b所示的磁异常模量的反演结果，黑框代表真实模型的位置(仅显示磁化率为0.0042SI以上的值)(a)和(b)利用模型空间反演算法的Z=-125 m和N=500 m的反演结果切片图；(c)和(d)利用数据空间反演算法的Z=-125 m和N=500 m的反演结果切片图.Fig.7 Inversion results of the amplitude data in Fig.5b. The black lines indicate the position of the true model. Only values above 0.0042SI are displayed(a，b) Model-space inversion results, slice at 125 m depth and 500 m north;(c，d) Data-space inversion results, slice at 125 m depth and 500 m north.

图8 (a) 实测航磁数据；(b) 由磁异常转换得到的模量Fig.8 (a) The measured aeromagnetic data；(b) The amplitude data computed from the total field anomaly in Fig.8a

可以看出，图中存在三个磁化方向明显不同的叠加异常，采用单一磁化方向对该数据进行反演显然是不合适的.由磁异常换算得到的模量如图8b所示.

在模量反演中，为了减少边界效应，地下剖分为63×81×27的立方体，即模型的边界范围比采集的数据边界范围稍大.同时为了便于对比，仅显示中间53×71×27的反演结果.图9a、9b和9c、9d分别为采用模型空间反演算法和数据空间反演算法的反演结果切片图.从图9a和9c可以看出，两种方法均显示了三个场源的平面位置，数据空间反演的结果要更加聚焦一些.对于北侧异常而言，模型空间显示为一个场源，而数据空间则显示为两个.由于没有钻孔数据，还无法进行直接验证，但该反演结果更细致地展示了场源分布的细节.图9b和9d均显示了北侧场源的产状，模型空间反演的结果模糊一些，而数据空间反演的结果则更加聚焦，产状也更加明显.

该实测数据的数据空间反演共需187次CG迭代，耗时600.7 s.对应的模型空间反演共需2402次CG迭代，耗时4509.3 s.

图9 实测数据的反演结果，虚线代表切片的位置(仅显示磁化率为0.005SI以上的值)(a)和(b)利用模型空间反演算法的Z=-425 m和N=6493825 m的反演结果切片图； (c)和(d)利用数据空间反演算法的Z=-425 m和N=6493825 m的反演结果切片图.Fig.9 Inversion results of the amplitude data in Fig.8b. The dash lines indicate the position of the slices. Only values above 0.005SI are displayed(a，b) Model-space inversion results, slice at 425 m depth and 6493825 m north； (c，d) Data-space inversion results, slice at 425 m depth and 6493825 m north.

5 结论

强剩磁的存在改变了磁化强度的方向，进而影响了磁异常的反演和解释.磁异常模量是一种受磁化方向影响小的转换量，且由磁异常转换模量的过程中能达到不放大噪声，保持了长波长信息.我们将数据空间反演算法扩展到了磁异常模量反演中，模型试验表明，与模型空间反演算法相比，其具有计算量小，成像结果聚焦的优势.但该算法同时也存在反演结果的磁化率极大值可能偏大的缺点，这是在实际应用中需要注意的地方.

致谢 本文研究工作得到了加拿大地质调查局Mark Pilkington的帮助，使用的实测数据来源于澳大利亚地质调查局网站公布的资料，两位审稿专家提出了宝贵的修改建议，在此一并感谢.

附录A 数据空间磁异常模量反演算法

1 数据空间反演算法迭代公式的推导

磁性反演计算，就是对目标函数(2)式求解极小值的过程.在目标函数的极小值处，(2)式的梯度一定为0，故可得到：

=0,

(A1)

其中，J为F(m)的雅可比矩阵.我们对(A1)式进一步转化，可以得到适合计算的具体公式.由(A1)式整理得W-1(m-m0)=JTD-1(d-F(m))，两边加上JTD-1J(m-m0)则得到：

(JTD-1J+W-1)(m-m0)=

JTD-1(d-F(m)+J(m-m0)),

(A2)

再由(A2)式两边乘上(JTD-1J+W-1)-1，得到：

m=m0+(JTD-1J+W-1)-1JTD-1(d-F(m)

+J(m-m0)).

(A3)

在(A3)式中，JTD-1J+W-1为M×M阶矩阵，即存在于模型空间，该矩阵的计算如求逆等，计算量通常会是非常大的.为此，设法对其进一步转换.

由于(JTD-1J+W-1)和JTD-1满足下列恒等式：

JTD-1(D+JWJT) =JT+JTD-1JWJT

=(JTD-1J+W-1)WJT,

(A4)

且(D+JWJT)和(JTD-1J+W-1)均可逆，由(A4)式则可以得出：

(JTD-1J+W-1)-1JTD-1=WJT(D+JWJT)-1,

(A5)

进一步将(A5)式代入(A3)式中进行等式替换即得到了方程在数据空间中的不动点形式：m=m0+WJT(D+JWJT)-1(d-F(m)+J(m-m0)),

(A6)

其中，D+JWJT为N×N阶矩阵，即存在于数据空间.与(A3)式比较，该矩阵求逆等计算的计算量将大大减少.

2 雅可比矩阵J的推导

由于磁异常模量与地下模型的磁化率(或磁化强度)之间是非线性关系，因此需要在迭代的过程中计算其雅可比矩阵.首先将磁异常三分量的正演过程表达为：

Hax=Gxκ,

Hay=Gyκ,

Za=Gzκ,

(A7)

其中，Hax、Hay、Za分别为磁异常向量的x、y、z分量.Gx、Gy、Gz是磁异常三分量正演的灵敏度矩阵.κ=m2为有效磁化率向量.

在迭代的过程中，Ta的第i个分量Tai对第j个模型mj的偏导数为：

×2mj,

(A8)

由(A8)式可以得到雅可比矩阵J：

×diag(2m),

(A9)

其中，diag(Hax/Ta)表示主对角由Hax/Ta向量构成的对角矩阵.diag(Hay/Ta)，diag(Za/Ta)以及diag(2m)的定义类似.

这样在反演迭代的过程中，仅需存储Gx，Gy，Gz三个矩阵即可快速完成相关的矩阵向量运算.计算这三个矩阵需要正确的磁化方向，但真实的磁化方向是未知的，由于磁异常模量受磁化方向的影响小，因此在计算这三个矩阵时，可以采用地磁场方向作为假设的磁化方向，这是模量反演中需要的(Li et al.，2010).显然，这个磁化方向的选择已经不再像常规磁异常分量反演那么重要了.

3 数据空间反演步骤

为便于了解在数据空间进行反演的计算细节，这里将具体计算步骤列出如下：

(6)mi+1=mi+Δmi；

Dampney C N G. 1969. The equivalent source technique.Geophysics, 34(1): 39-53.Dannemiller N, Li Y G. 2006. A new method for determination of magnetization direction.Geophysics, 71(6): L69-L73.

Ellis R G, de Wet B, Macleod I N. 2012. Inversion of magnetic data from remanent and induced sources. ∥ 22nd International Geophysical Conference and Exhibition. Brisbane, Australia.

Gerovska D, Aráuzo-Bravo M J, Stavrev P. 2009. Estimating the magnetization direction of sources from southeast Bulgaria through correlation between reduced-to-the-pole and total magnitude anomalies.GeophysicalProspecting, 57(4): 491-505.

Lelièvre P G, Oldenburg D W. 2006. Magnetic forward modelling and inversion for high susceptibility.GeophysicalJournalInternational, 166(1): 76-90.

Lelièvre P G, Oldenburg D W. 2009. A 3D total magnetization inversion applicable when significant, complicated remanence is present.Geophysics, 74(3): L21-L30.

Li S L, Li Y G, Meng X H. 2012a. The 3D magnetic structure beneath the continental margin of the northeastern South China Sea.AppliedGeophysics, 9(3): 237-246.

Li S L, Li Y G. 2014. Inversion of magnetic anomaly on rugged observation surface in the presence of strong remanent magnetization.Geophysics, 79(2): J11-J19.

Li Y G, Oldenburg D W. 1996. 3-D inversion of magnetic data.Geophysics, 61(2): 394-408.

Li Y G, Oldenburg D W. 1998. 3-D inversion of gravity data.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.GeophysicalJournalInternational, 152(2): 251-265.

Li Y G, Shearer S E, Haney M M, et al. 2010. Comprehensive approaches to 3D inversion of magnetic data affected by remanent magnetization.Geophysics, 75(1): L1-L11.

Li Y G, He Z X, Liu Y X. 2012b. Application of magnetic amplitude inversion in exploration for volcanic units in a basin environment.Geophysics, 77(5): B219-B225.

Liu S, Feng J, Gao W L, et al. 2013. 2D inversion for borehole magnetic data in the presence of significant remanence and demagnetization.ChineseJ.Geophys. (in Chinese), 56(12): 4297-4309, doi: 10.6038/cjg20131232.

Liu S, Hu X Y, Liu T Y, et al. 2013. Magnetization vector imaging for borehole magnetic data based on magnitude magnetic anomaly.Geophysics, 78(6): D429-D444.

Lourenco J S, Morrison H F. 1973. Vector magnetic anomalies derived from measurements of a single component of the field.Geophysics, 38(2): 359-368.

Nocedal J, Wright S. 2006. Numerical Optimization. New York: Springer.

Pedersen L B. 1978. Wavenumber domain expressions for potential fields from arbitrary 2-, 21/2-, and 3-dimensional bodies.

Geophysics, 43(3): 626-630.

Phillips J D. 2005. Can we estimate total magnetization directions from aeromagnetic data using Helbig′s integrals?Earth,PlanetsandSpace, 57(8): 681-689. Pilkington M. 1997. 3-D magnetic imaging using conjugate gradients.Geophysics, 62(4): 1132-1142. Pilkington M. 2009. 3D magnetic data-space inversion with sparseness constraints.Geophysics, 74(1): L7-L15. Pilkington M, Beiki M. 2013. Mitigating remanent magnetization effects in magnetic data using the normalized source strength.Geophysics, 78(3): J25-J32.

Portniaguine O, Zhdanov M S. 2002. 3D magnetic inversion with data compression and image focusing.Geophysics, 67(5): 1532-1541.

Santos M L, Li Y G, Moraes R. 2012. Application of 3D magnetic amplitude inversion to Fe oxide-Cu-Au deposits at low magnetic latitudes: A case study from Carajás Mineral Province, Brazil.∥ 2012 SEG Annual Meeting. Society of Exploration Geophysicists. 1-5.

Shearer S E. 2005. Three-dimensional inversion of magnetic data in the presence of remanent magnetization [Master′s thesis]. Colorado: Colorado School of Mines.

Stavrev P, Gerovska D. 2000. Magnetic field transforms with low sensitivity to the direction of source magnetization and high centricity.GeophysicalProspecting, 48(2): 317-340.

Tarantola A. 1987. Inverse Problem Theory. Methods for Data Fitting and Model Parameter Estimation. New York: Elsevier.

Wang M Y, Di Q Y, Xu K, et al. 2004. Magnetization vector inversion equations and 2D forward and inversed model study.ChineseJ.Geophys. (in Chinese), 47(3): 528-534, doi: 10.3321/j.issn:0001-5733.2004.03.025.

Yao C L, Hao T Y, Guan Z N, et al. 2003. High-speed computation and efficient storage in 3-D gravity and magnetic inversion based on genetic algorithms.ChineseJ.Geophys. (in Chinese), 46(2): 252-258.

Yao C L, Zheng Y M, Zhang Y W. 2007. 3-D gravity and magnetic inversion for physical properties using stochastic subspaces.ChineseJ.Geophys. (in Chinese), 50(5): 1576-1583.

附中文参考文献

刘双, 冯杰, 高文利等. 2013. 强剩磁强退磁条件下的二维井中磁测反演. 地球物理学报, 56(12): 4297-4309, doi: 10.6038/cjg20131232.

王妙月, 底青云, 许琨等. 2004. 磁化强度矢量反演方程及二维模型正反演研究. 地球物理学报, 47(3): 528-534, doi: 10.3321/j.issn:0001-5733.2004.03.025.

姚长利, 郝天珧, 管志宁等. 2003. 重磁遗传算法三维反演中高速计算及有效存储方法技术. 地球物理学报, 46(2): 252-258.

姚长利, 郑元满, 张聿文. 2007. 重磁异常三维物性反演随机子域法方法技术. 地球物理学报, 50(5): 1576-1583.

(本文编辑 何燕)

3D data-space inversion of magnetic amplitude data

LI Ze-Lin, YAO Chang-Li*, ZHENG Yuan-Man, MENG Xiao-Hong, ZHANG Yu-Wen

KeyLaboratoryofGeo-detection(ChinaUniversityofGeosciences,Beijing),MinistryofEducation,StateKeyLaboratoryofGeologicalProcessesandMineralResources,Beijing100083,China

3D magnetic inversion for susceptibility distribution is a powerful tool in quantitative interpretation of magnetic data and has been successfully applied to exploration of mineral and oil gas resources and to interpretation of regional geologic structure. The traditional inversion algorithms require knowledge of magnetization direction, which means that one should assume there is no remanent magnetization and self-demagnetization. Consequently, the direction of magnetization is assumed to be the same as the current geomagnetic field direction. However, strong remanent magnetization always exists and the total magnetization direction can be significantly different from that of the geomagnetic field direction due to complicated geological conditions, and in this case the traditional inversion algorithms become ineffective and the inversion result may be wrong.Magnetic amplitude data are less sensitive to the total magnetization direction and can be used to invert for 3D magnetization strength distribution to delineate the positions of causative bodies in the presence of strong remanent magnetization. We present a data-space magnetic amplitude inversion algorithm to reduce the effects of remanent magnetization. We also give a detail formula derivation of the proposed algorithm. In the data space, the matrix dimensions are equal toN×N(the number of data) rather thanM×M(the number of model parameters), whereNis far less thanM. So the computational efficiency is improved. The computational time is further reduced because this method does not need to search for a regularization parameter by using an incomplete conjugate gradient method. Moreover, a square root operator is used to impose a positivity constraint on the effective susceptibility and likewise to focus the inversion results.Three synthetic data examples and a real data example are used to verify the proposed data-space algorithm. Tests on these examples show that the proposed method can focus the inversion results and likewise increase solution speed by up to an order of magnitude compared with the standard model-space, L2-norm regularized inversion.

Data-space; Magnetic amplitude data; Remanence; 3D inversion

10.6038/cjg20151030.

Li Z L, Yao C L, Zheng Y M, et al. 2015. 3D data-space inversion of magnetic amplitude data.ChineseJ.Geophys. (in Chinese),58(10):3804-3814,doi:10.6038/cjg20151030.

国家高技术研究发展计划(863计划)项目(2014AA06A613)和国家自然科学基金项目(41304104)资助．

李泽林，主要从事重磁处理反演研究. E-mail：zelin.lee@foxmail.com

*通讯作者 姚长利，从事重磁勘探理论与方法技术研究.E-mail：clyao@cugb.edu.cn

10.6038/cjg20151030

P631

2014-10-24，2015-03-05收修定稿

≤≥≪ ≫李泽林, 姚长利, 郑元满等. 2015. 数据空间磁异常模量三维反演.地球物理学报，58(10):3804-3814,