2017, Vol. 60
磁电法是一种测量磁场的传导类电法勘探方法，磁场由直流或者伪直流的人工导线源所激发 (Chen et al., 2002).不同于传统电阻率法中测量地面上两点之间的电位差，磁电法中改为观测各点的磁场，通过磁场的空间变化从而研究地下电导率变化 (邓靖武，2005).其中异常磁场由流经地下目标体的电流所产生，是其几何形状和相对电导率的函数 (Bouchedda et al., 2015).磁电法可以分为磁电阻率法 (MMR)(Edwards，1974；Edwards and Nabighian, 1991；Chen et al., 2002；Chen and Oldenburg, 2004)、磁激发极化法 (MIP)(Seigel，1974；Howland-Rose et al., 1980a, b；Chen and Oldenburg, 2003；李建华等，2010) 和观测总磁场强度的亚音频磁测法 (SAM)(Cattach et al., 1993；Fathianpour and Cattach, 1995；Fathianpour et al., 2005a, 2005b).MMR除了地面观测之外，在井中勘探 (Acoata and Worthington, 1983；Bishop et al., 1997；Godber and Bishop, 2007) 和海洋勘探 (Nobes et al., 1986；Yang，2005；翁爱华等，2009；Seama et al., 2013) 中都有较多应用，而MIP和SAM方法则主要被应用在地面勘探.地面磁电阻率法不仅被应用在传统的地质结构调查和找矿领域 (Edwards and Howell, 1976；Meju，2002)，近年来在工程勘察和污染物监测等领域也取得了良好的应用效果 (Zhu and Yang, 2008；Kofoed et al., 2011；Maxwell et al., 2013).
磁电阻率法的反演技术发展比较缓慢.有学者用典型曲线法对一些简单的二维目标体进行定性分析 (Howland-Rose et al., 1980b)，通过重力和磁电阻率之间的相似性用试错法进行数据解释 (Szarka，1987).Asten (1988)也利用这种关系，用标准二维重力模拟的方法处理井中磁电阻率数据.Chen等 (2002)采用正则化的反演思想率先实现了对于磁电阻率数据的三维反演，对于目标函数的求解采用共轭梯度最小二乘法 (CGLS).LaBrecque等 (2003)采用Occam方法联合反演由同一低频源激发的ERT和MMR数据.对于地面磁电阻率数据反演而言，最大的问题在于异常体的深度分辨能力较差，异常形态浮于地表，而浅部异常体分辨能力较高 (Chen et al., 2002).同样，在低频可控源资料反演中，只采用准静态磁场Hy数据反演浅部效果好，深部效果较差 (翁爱华等，2015).Li和Oldenburg (1996)认为地面磁场数据对异常体深度分辨能力的不足是因为核函数随深度的快速衰减造成的，他们引入了权函数w(z) 来抵消核函数的几何扩散从而提高深部分辨能力.随后这种思路又被用在了重力场数据反演中 (Li and Oldenburg, 1998)，Chen等 (2002)稍作变形引入到地面磁电阻率数据反演中，但难以从核函数中分离出明显的衰减因子，并且权函数的参数设定依据经验估计与数值模拟实验，不同的参数会影响反演结果 (Li and Oldenburg, 1996，1998；Chen et al., 2002).
针对以上这些问题，本文提出基于频率域Helmholtz方程的地面磁电阻率数据三维非线性共轭梯度反演技术，在反演过程中结合印模法以提高地面磁电阻率法的深部分辨能力.印模法是叶涛等 (2013)在二维大地电磁资料反演时为压制静态效应和初始模型对于反演的影响提出的.他们将初始模型的某个深度 (叶涛等称其为印模深度) 以下设置为均匀半空间 (辅模型)，印模深度以上将已有的反演结果 (元模型) 和均匀半空间采用加权求和来确定电阻率值，并在此基础上进行多次的迭代重构反演，直到达到理想的效果.孔志召等 (2015)在AMT数据的二维反演中也采用印模法，他们的方法为根据一维TM模式的平均电阻率构建元模型和辅模型，根据理论探测深度确定印模深度.印模法反演结果既保留了已有反演结果的主要轮廓，又在深部通过平均作用消除了“挂面条”现象.考虑到地面磁电阻率法反演深部效果差、浅部效果好，与大地电磁反演纵向分辨能力具有对称性的特点，本文借鉴印模法和迭代重构的反演思想以解决地面磁电阻率数据反演中对异常体深度分辨能力不足的问题.在应用过程中，提出了新的初始模型重构方式，并明确了迭代重构的终止条件，即将上次重构反演最后一次迭代模型变化量非0的区域的纵向分界面作为印模深度.印模深度以上替换成均匀半空间，以下为上一次的反演结果，作为下一次重构的初始模型.当迭代重构结束时模型变化量为0或者两次迭代重构结束时的模型变化量的差值小于阈值时迭代重构结束.简单来讲，本文印模法的使用目的是通过初始模型对反演结果加以影响，使得反演效果符合实际情况，这是与大地电磁中应用印模法最本质的区别.
Chen和Oldenburg (2004)提出，在充分低频的情况下，可用一维频率域长导线源来解决MMR的正演问题.基于此，本文讨论的MMR磁场数据利用低频交流导线源磁场近似.因此，本文的研究思路如下：发射源用频率域低频载流导线代替直流导线，三维正演基于考虑导电效应的Helmholtz方程，数值上采用三维有限差分算法，从而建立起磁场与地面电导率的直接关系；反演采用三维非线性共轭梯度 (NLCG) 反演技术，反演过程中引入印模法以提高深度分辨能力；最后，将上述技术应用到合成数据与实际数据的反演计算中.
2 正演模拟 2.1 磁场计算方法的选择在测量电场的直流电法中，我们可以通过地面测得的电位差得到地下介质的视电阻率.但对于磁电阻率法，测得的磁场数据只与介质间的电导率之比有关，而与绝对电导率无关.磁电阻率法中磁场的计算无论是从毕奥-萨伐尔定律 (Edwards et al，1978；Yang and Tseng, 1992) 或者静电场的麦克斯韦方程组入手 (Chen et al., 2002；Chaladgarn and Yooyuanyong, 2013)，都将其看作是静磁场问题.由附录A可见，假设任意两种不同电导率的模型满足σ1(x, y, z)=χσ2(x, y, z)，其中为比例系数.由式 (A3) 可见ϕ2=ϕ1/χ，也就是σ1∇ϕ1=σ2∇ϕ2.根据式 (A4) 可见，不同的电导率模型对于磁场强度并没有影响.但如果从频率域Helmholtz方程出发 (见式 (A5))，我们可以得到磁场与电导率之间的关系，从而已知磁场就可以求得电导率值，而不是电导率之比 (见式 (A6)).而且在实际工作中，磁电阻率法也是采用低频作为发射信号 (Chen et al., 2002；李建华等，2010)，考虑频率信息更符合实际情况.所以，我们由频率域Helmholtz方程出发对磁场进行计算.
2.2 控制方程与分量选择三维磁场计算通过异常场满足的Helmholtz方程求解，公式为
|
(1) |
其中En为层状背景介质产生的基本 (一次) 场，Ea为三维异常体产生的异常 (二次) 场，kn2是背景介质复波数，k2是含异常后介质复波数.导线源的计算采用虚界面方法由偶极子迭加得到 (Weng et al., 2016).然后采用交错网格有限差分技术离散方程 (1)(Yee，1966；Alumbaugh et al., 1996).这样，总场E可以表示为
|
(2) |
首先计算电场，磁场由电场计算得到.而剖分单元面中心处磁场B则由该单元四周边界中心总电场差分得到，公式为
|
(3) |
其中Db是差分算子，而地面任意观测点处的磁场Bobs由B经过空间差值得到 (Egbert and Kelbert, 2012)，在实际测量中由磁场三分量数据组成 (以源方向为x方向)(Kofoed et al., 2011)，公式为
|
(4a) |
|
(4b) |
|
(4c) |
其中Bn是背景导电介质中电流形成的背景磁场，Ba是由局部异常扰动背景电流密度引起的异常磁场，Bw为导线源中的电流在自由空气中激发的磁场.
从异常场的场值大小来讲，垂直分量Baz与y方向分量Bay相差不多，但都大于x方向分量Bax，而z分量由于受地形影响噪声水平要大于y分量 (Chen et al., 2002)，所以本文采用磁场的y分量进行正反演计算.采用y分量还可以克服Bw的影响.在传统的地面磁电法测量中多采用U型布源 (Chen et al., 2002；Fathianpour et al., 2005b；李建华等，2010)，主要目的是为了减小线场对于测区的干扰.但如 (4) 式可见，线场只对观测数据的垂直分量有影响，观测到的磁场水平分量可以不加处理，垂直分量减去线场的解析解 (Asten, 1988) 即可得到异常场.因此，可以采用长直导线源作为磁电阻率法发射源，既省去实际施工中布设U型源的不便之处，同时也方便计算.
2.3 精度验证为了验证本文算法模拟精度，我们设计了如图 1所示的三维仿球体模型.将其磁场的计算结果与球体解析结果 (傅良魁，1984) 对比.等效球体半径为30 m，中心点距离地面60 m，球体电阻率为10 Ωm，背景电阻率为100 Ωm.场源为x方向的60 m长电偶极子，位于异常体正上方，发射频率采用0.3 Hz.对计算区域采用长方体单元进行离散，网格数为40×40×30.图 2所示为球体正上方x与y方向主剖面的精度验证结果.由图中可见，模拟结果与解析解的结果基本一致.主要误差出现在异常体正上方，产生的原因在于三维仿球体模型无法精确模拟球体形状.因此，由低频的频率域发射源模拟磁电阻率法发射源是可行的.
|
图 1 正演结果验证模型 (a) 球体解析解模型；(b) 三维仿球体模型. Fig. 1 Models to verify the precision of forward modeling results (a) Analytic solution sphere model; (b) 3D sphere-simulation model. |
|
图 2 正演验证结果 (a) x轴主剖面；(b) y轴主剖面. Fig. 2 Verification of forward modeling results (a) x axis principal section; (b) y axis principal section. |
Newman和Alumbaugh (2000)率先采用非线性共轭反演技术对三维大地电磁数据进行反演.其后陆续被应用到地面可控源 (林昌洪等，2012；翁爱华等, 2012, 2015) 和航空电磁 (刘云鹤和殷长春，2013) 的三维反演中，都取得了良好的效果.本文将该反演技术用于反演地面磁电阻率数据.具体计算流程如下 (Newman and Alumbaugh, 2000):
(1) 令i=1，mi=mref并计算ri=－∇φ(mi)，mref是初始模型；
(2) 如果‖ri‖满足要求则停止计算，否则令ui=Mi-1ri，M为预处理算子，u是指定的共轭搜索方向；
(3) 搜索步长αi使得φ(mi+αui) 最小化.当步长αi太小时减小正则化因子λ，当λ达到阈值的时候则反演终止；
(4) 令mi+1=mi+αiui，ri+1=－∇φ(mi+1)；
(5) 如果‖ri+1‖满足要求则停止计算，否则跳到 (6)；
(6) 令
(7) 令
(8) 令i=i+1，返回到步骤 (3).
3.2 地面磁场数据的深度分辨能力设计图 3所示的单一异常体模型检验地面磁电阻率数据非线性共轭梯度反演效果.针对磁场y分量进行反演.发射源采用1200 m长直导线，用中梯方式观测，导线沿x方向敷设，测区范围为源中心1/3，400 m×400 m，共21条测线，441个测点，线距和点距都为20 m.发射频率为0.3 Hz，电流1 A.计算区域被离散成长方体网格，网格数为40×40×30.在理论数据中加入5%的高斯噪声来进行反演研究.
|
图 3 单一三维异常体模型 (a) 过异常中心剖面；(b) 俯视图. Fig. 3 Model with a single 3D conductive cube buried in a half space. (a) Cross section through cube center; (b) Overlook of the model. |
图 4a为反演结果.从结果可以看出，基于Helmholtz方程，用地面准静磁场数据反演，相对于基于静磁场正演进行反演的方法，可以获得电阻率值的分布；但是反演的结果主要集中在地面，难以获得异常体的埋深，和前人的研究结果一样 (Li and Oldenburg, 1996，2000；Chen et al., 2002；Zhdanov，2002).
由毕奥-萨伐尔定律计算磁场B公式为
|
(5) |
|
图 4 印模法初始模型重构示意图 (xoz剖面) (a) 元模型；(b) 辅模型；(c) 重构后的初始模型. Fig. 4 Starting model reconstruction of impressing method (xoz profile) (a) Meta model; (b) Auxiliary model; (c) Starting model after reconstruction. |
其中μ是磁导率，J是电流密度.从式 (5) 可见，该方程属于第一类Fredholm积分方程 (IFKs)，核函数为∇(1/r)，随着深度的增加，核函数迅速地衰减 (Aster et al., 2005)，导致静磁场的深度分辨能力不佳.而且，根据磁场高斯定律可知，通过任意闭合曲面的磁通量恒等于0.这样，对于地面静磁场数据而言，如果只是知道边界曲面上场的分布，我们是无法确定曲面内具体的异常源分布情况.除此之外，磁场观测数据的有限性和精度同样影响深度分辨能力 (Zhdanov，2002).
3.3 适用于地面磁电阻率数据三维反演的印模法及迭代重构过程为了解决图 3所示异常体的深度分辨能力的问题，在反演中引入印模法的思想，在深部 (印模深度以下) 我们保留上次反演结果的深部信息，将其看做先验信息保留在下一次反演的初始模型中；在浅部 (印模深度以上)，用均匀半空间替代效果不好的反演结果.这样已有的反演结果在下，均匀半空间在上，组合出新的反演初始模型，得到反演结果后再按此方法组合，作为下一次反演的初始模型，如此进行多次反演, 重构的过程如图 4所示.
下面介绍该方法的数学原理.根据Egbert和Kelbert (2012)的阐述，电磁反演目标函数ϕ(m, d) 可以表示为
|
(6) |
其中m是电导率反演结果向量，mref是电导率初始模型向量，d是数据向量，f(m) 是正演算子.Cd和Cm分别是数据方差矩阵和模型方差矩阵，λ是正则化因子.对 (6) 式进行拟线性化反演，使其最小化的解为
|
(7) |
其中G是正演数据核矩阵.(7) 式为初始模型mref与反演结果m之间的关系.
由式 (7) 可见，对于给定的数据向量d，当使用相同的正反演方法时，反演结果m与初始模型mref相关.为了方便表述，我们定义两个关于数据和模型的平衡矩阵D和M为
|
(8) |
|
(9) |
式 (7) 可以表述为
|
(10) |
其中mref1为第一次反演的初始模型，在缺乏先验信息的情况下可以设置为均匀半空间，m1是第一次反演后得到结果.我们将第一次反演的结果和初始模型根据印模深度 (下节讨论) 进行组合，组合的过程可以用函数Ψ1表示.这样我们就得到了用于第二次反演的初始模型mref2为
|
(11) |
因为第一次反演结果m1的加入，此时的初始模型mref1与最开始的初始模型mref2相比相当于添加了先验地质信息 (当然如果确实有先验信息的话就可以在Ψ1中添加).我们使用初始模型mref1同样进行如式 (10) 的反演，可以得到第二次的反演结果m2，我们将第二次反演的结果和初始模型再根据印模深度进行组合，组合的过程认为是函数Ψ2.这样的过程称之为迭代重构.综上所述，印模法的具体流程为
|
(12) |
|
(13) |
其中mj是第j次迭代重构反演后的结果，mrefj是第j次迭代重构所用的初始模型，Ψj是用于求得第j+1次反演初始模型mrefj+1的函数.具体到地面磁电阻率数据的反演，因为主要目的是为了提高深度分辨能力，函数Ψj是将印模深度h以上设置为均匀半空间，印模深度h以下为已有的反演结果，具体为
|
(14) |
其中a、b和c分别为三维交错网格有限差分法形成的电导率立体数据块x、y和z三个方向上的数据量.因此，在印模法的反演过程中，每次迭代都涉及到三个模型.如图 4所示，(a) 为已有的反演结果，称之为元模型，(b) 为均匀半空间模型，称之为辅模型，(c) 为根据式 (14) 组合形成的用于下一次反演的初始模型.印模深度的选择和迭代重构的终止条件通过模型变化量来确定，在下节具体讨论.现行的三维反演的主要局限性在于计算时间成本过大 (汤井田等，2007)，但是相对于可控源中的多频反演，磁电阻率法的发射源只是直流或者较低的单一频率人工源，反演受时间的限制并不大，这是可以进行迭代重构反演的一个重要前提.
3.4 印模深度和迭代重构终止条件的确定印模法反演的一个关键问题是印模深度的确定.叶涛等 (2013)利用过渡区的概念来确定印模深度.他们认为过渡区是受数据约束程度由强变弱再到无的区域，在反演结果中表现为电阻率逐渐向初始模型值靠近的区域.当元模型的深部电阻率值大于真实模型的基底电阻率值时选择过渡区的极小值点作为印模深度，反之则选择极大值点.与二维大地电磁中印模法构建初始模型正好相反，本文中三维地面磁电阻率数据反演中重构的初始模型浅层为均匀半空间，深部为已有的反演结果.所以印模深度的选择方法必然也不相同.本文通过模型变化量来选择印模深度.
下面以非线性共轭梯度反演为例说明其原理.根据2.1节步骤 (4) 中所示：
|
(15) |
模型变化量Δm为每次更新步长和共轭梯度方向之后两次反演结果之间的差值.可见，当一次反演结束时，模型变化量Δm较大的区域对应的是模型修正空间较大的区域.此时虽然反演结束，但是还有进一步反演的空间，这个区域也是我们下一步通过重构模型将要改造的区域.Δm为0的区域我们认为是拟合效果已经不错的区域.在地面磁电阻率数据反演中，主要是浅层的模型变化量较大，深部较小.因此所以我们选择印模深度时的原则是保留Δm为0的区域，将Δm非0的区域替换为均匀半空间，在下一次迭代重构时用地面数据对其重新拟合.
随着迭代重构次数的增加，虽然反演的效果越来越好，但是其提升的程度越来越小 (叶涛等，2013).当Δm全部为0时可以认为模型没有再进一步修正的空间了，迭代重构可以终止.当Δm不为0时，我们通过前后两次迭代重构终止时Δm的差值来确定迭代重构是否终止.如果差值小于给定的阈值ε，则迭代重构同样终止，如式 (16) 所示：
|
(16) |
其中Δmj是第j次重构初始模型后反演终止时的模型变化量；Δmj+1是第j+1次重构初始模型后反演终止时的模型变化量.
综上两节所述，使用印模法及迭代重构反演流程如图 5所示.
|
图 5 反演流程图 Fig. 5 Flowchart of inversion |
图 6所示为图 3模型的磁电阻率地面数据y分量的初始反演及每次迭代重构后的反演结果，图 7为相对应的每次重构结束时的模型变化量.基于上文的分析，我们选择模型变化量非0区域的底界面处为印模深度 (见图 7).由2.1节的步骤3和步骤4可见，每次反演的终止条件取决于目标函数的梯度‖ri‖和正则化因子λ是否达到阈值.为了基于同一的标准比较每次迭代重构后的反演效果，我们设定每次迭代重构时‖ri‖和λ的阈值不变.这样，对于相同的地面数据而言，反演结果只与初始模型有关.从图 6a—h中可见，每一次迭代重构都使得异常体的深度位置有所下降.由图 7可见，对于单一异常体而言，随着迭代重构的进行，每次反演结束时的模型变化量越来越小，当第7次迭代重构结束时模型变化量为0(见图 7h)，此时迭代重构终止，异常体的埋深与模型情况良好符合 (见图 6h).并且，随着迭代重构的进行，低阻异常体的视电阻率也逐渐下降.初始反演结束时异常体中心的视电阻率为150 Ωm (见图 6a)，第7次迭代重构后，异常体的中心电阻率下降到55 Ωm (见图 6h).图 6i为第7次迭代重构后xoz剖面图，图 4a为对应的初始反演xoz剖面图.与图 6h对比可见异常体在埋深上同样反映良好，差别在于x方向上的横向分辨能力不如y方向上的横向分辨能力.可见，低频磁场y分量在y方向上的横向分辨能力要明显好于x方向上的横向分辨能力，不过用印模法同样可以准确确定异常体的埋深位置.综上所述，地面磁电阻率数据经过印模法构造初始模型并进行迭代重构反演后，异常体埋深的位置更加准确，视电阻率的分辨能力也有明显提高.
|
图 6 单一异常体初始反演及每次迭代重构后的反演结果 (a) 初始反演yoz剖面；(b) 迭代重构Ⅰ，印模深度120 m，yoz剖面；(c) 迭代重构Ⅱ，印模深度120 m，yoz剖面；(d) 迭代重构Ⅲ，印模深度120 m，yoz剖面；(e) 迭代重构Ⅳ，印模深度120 m，yoz剖面；(f) 迭代重构Ⅴ，印模深度100 m，yoz剖面；(g) 迭代重构Ⅵ，印模深度100 m，yoz剖面；(h) 迭代重构Ⅶ，印模深度80 m，yoz剖面；(i) 迭代重构Ⅷ，印模深度80 m，xoz剖面. Fig. 6 Anomalous body′s inversion results of initial inversion and each iteration (a) Initial inversion, yoz; (b) 1st iteration, h=120 m, yoz; (c) 2nd iteration, h=120 m, yoz; (d) 3rd iteration, h=120 m, yoz; (e) 4th iteration, h=120 m, yoz; (f) 5th iteration, h=100 m, yoz; (g) 6th iteration, h=100 m, yoz; (h) 7th iteration, h=80 m, yoz; (i) 7th inversion, h=80 m, xoz. |
|
图 7 初始反演及每次迭代重构结束时的模型变化量 (yoz剖面，x=0 m) (a) 初始反演；(b) 迭代重构Ⅰ，印模深度120 m；(c) 迭代重构Ⅱ，印模深度120 m；(d) 迭代重构Ⅲ，印模深度120 m；(e) 迭代重构Ⅳ，印模深度120 m；(f) 迭代重构Ⅴ，印模深度100 m；(g) 迭代重构Ⅵ，印模深度100 m；(h) 迭代重构Ⅶ，印模深度80 m. Fig. 7 Model variations of initial inversion and each iteration finishes (yoz profile, x=0 m) (a) 1st inversion; (b) 1st iterative, h=120 m; (c) 2nd iteration, h=120 m; (d) 3rd iteration, h=120 m; (e) 4th iteration, h=120 m; (f) 5th iteration, h=100 m; (g) 6th iteration, h=100 m; (h) 7th inversion, h=80 m. |
为了进一步验证本方法的有效性，本文设计了如图 8a和图 8b所示的含浅层不均匀导电体的复杂异常模型, 其中L1和L2剖面分别为过浅层导电体A和B中心的yoz方向剖面.观测方式、发射频率和电流、网格剖分和噪声水平等设置与图 3单一异常体模型相一致.设计这样的模型出于两方面的考虑：首先，印模法中浅层的初始模型采用均匀半空间可以有效的压制假异常，提升深部分辨能力，但要确定其是否压制了浅层的电性异常；其次，在可控源与大地电磁勘探中，静态效应的存在是影响勘探效果的一个重要因素 (雷达，2010；李鑫等，2016)，磁电阻率法测磁场而不是测电场，受静态效应的影响可能会较小.
|
图 8 含浅层不均匀导电体的异常体模型及反演结果 (印模深度=120 m) (a) 模型xoz剖面，y=0 m；(b) 模型xoy平面，z=0 m；(c) 初始反演L1剖面；(d) 初始反演L2剖面；(e) 初始反演xoz剖面，y=0 m；(f) 迭代重构Ⅷ L1剖面；(g) 迭代重构Ⅶ L2剖面；(h) 迭代重构Ⅶ yoz剖面，y=0 m. Fig. 8 Anomalous prism model with shallow nonuniform conductor and inversion results (h=120 m) (a) xoz profile, y=0 m; (b) xoy profile, z=0 m; (c) 1st inversion, L1 profile; (d) 1st inversion, L2 profile; (e) 1st inversion, xoz, y=0 m; (f) 7th iteration, L1 profile; (g) 7th iteration, L2 profile; (h) 7th iteration, xoz, y=0 m. |
图 8c—e为初次反演结果，图 8f—h最后一次迭代重构的反演结果，因篇幅所限，本文略去中间迭代重构的反演结果.由图 8c—e可见，初始反演后异常仍然浮于地表，不过从异常形态上可以区分出浅层电性不均匀导电体A和B，而且靠近异常体的导电体B的异常更为明显，且有所下移.因为此时的反演效果为浅层导电体和深部异常体相叠加形成的，二者距离越近，这种叠加效果越明显.这也就解释了为什么埋深相对较深的导电体B的异常形态反而比距离测点更近的导电体A的异常形态明显.图 8f—h即为第7次迭代重构后的反演结果.对比初始反演结果，在浅层不均匀导电体的影响下依然能准确确定异常体的埋深，而且对于浅层导电体A和B，也能确定二者的位置.与单一异常体一样，y方向低频磁场数据同样是y方向上的横向分辨率要好于x方向上的横向分辨率.对比图 6图 8可见，浅层不均匀导电体对于深部异常体的影响主要是在不均匀导电体的走向上，而且不均匀导电体距离异常体越近则二者之间的相互影响越明显.这种相互影响使得反演结果中导电体B的埋深较实际要深 (见图 7g)，深部异常体在横向上向距离较近的导电体方向偏移 (见图 8h).图 9a和图 9b分别为初始反演和最后一次迭代重构结束时的模型变化量 (过导电体B中心的yoz剖面图)，同样略去中间迭代重构的过程.可见，因为浅层电性不均匀体的存在，模型变化量在深部上变化不大，所以在迭代重构中印模深度始终保持不变，选择模型变化量非0区域的底界面120 m作为印模深度.图 9c为第6次迭代重构与第7次迭代重构结束时模型变化量的差值，可见几乎为0，根据式 (16)，此时迭代重构结束.总体而言，本文所采用的方法在提高异常体的深度分辨能力的同时不会减弱对于浅层异常体的分辨能力，深部异常体受浅层不均匀导电体的影响也相对较小.
|
图 9 图 7模型初始反演及迭代重构结束时的模型变化量 (印模深度=120 m，L2剖面) (a) 初始反演；(b) 迭代重构Ⅶ；(c)︱迭代重构Ⅵ模型变化量-迭代重构Ⅶ模型变化量︱. Fig. 9 Model variations of first inversion and last iteration (h=120 m, L2 profile) (a) 1st inversion; (b) 7th iteration; (c) Model variation I of 6th iteration-model variation of 7th iteration. |
为了验证本文算法对于实测资料的有效性，对澳大利亚Mons Cupri矿区的地面MMR实测数据进行三维反演.此矿区分布有典型的侵染状Cu-Zn-Pb矿化带，埋深在20 m到200 m范围内 (Chen et al., 2002；Huston，2006).采用0.3 Hz磁场y分量进行反演.图 10为实测数据MMR响应的平面图.这里需要注意的是，一般所指的MMR异常为异常场Ba对测区中心点处的背景场Bn0归一化得到的百分比 (Chen et al., 2002；李建华等，2010)，公式为
|
(17) |
|
图 10 Mons Cupri矿区MMR响应 (百分比) 平面图 Fig. 10 MMR response (per cent) over the Mons Cupri deposit |
反演网格为30×30×25个 (不包括扩编网格)，在反演过程中采用自适应方式选择正则化因子 (翁爱华等，2012)，初始正则化因子为0.8.
图 11为实测资料的反演结果.其中图a—c为采用400 Ωm均匀半空间为初始模型的反演结果 (即初始反演)，图d—f为经过三次迭代重构后的反演结果，三次迭代重构的印模深度分别为60 m、60 m和30 m.图 11中A所示为低阻异常目标矿化带.从图中可见，利用印模法进行3次迭代重构后，与初始反演相比，高低阻异常更加明晰，与背景电阻率差异更大.深部分辨能力明显增强，不再像初始反演那样异常只是集中于地表.低阻异常体视电阻率为40 Ωm至90 Ωm.矿体空间位置、视电阻率与实际地质情况基本相符 (Chen et al., 2002；Huston，2006).
|
图 11 实测资料反演结果 (a) 初始反演，z=-80 m平面；(b) 初始反演，x=1150 m剖面；(c) 初始反演，y=950 m剖面；(d) 迭代重构Ⅲ，z=-80 m平面；(e) 迭代重构Ⅲ，x=1150 m剖面；(f) 迭代重构Ⅲ，y=950 m剖面，平面图中虚线为剖面所在位置，xoz剖面图中标出了目标矿化带A的位置. Fig. 11 Inversion results of measured data (a) Initial inversion, z=-80 m; (b) Initial inversion, x=1150 m; (c) Initial inversion, y=950 m; (d) 3rd iteration, z=-80 m; (e) 3rd iteration, x=1150 m; (f) 3rd iteration, y=950 m. Dashed lines indicate the location of profiles. A is location of mineralized zone. |
为了在地面磁电阻率法中反演出模型的绝对电阻率值，本文从考虑导电性的低频麦克斯韦方程出发，将非线性共轭梯度算法应用于地面MMR数据反演中，并通过引入印模法解决地面MMR数据分辨率不足的问题.通过合成数据与实际数据反演，得到了以下结论：
(1) 从频率域Helmholtz方程出发反演磁电阻法问题，可以得到绝对电导率值，而不仅仅是相对电导率值.
(2) NLCG与印模法相结合，通过迭代重构可用来对地面MMR数据进行三维反演，提高异常体的深度分辨能力.初始模型的重构方案为印模深度以上设置成均匀半空间，印模深度以下为上一次的反演结果.
(3) 根据每次重构反演结束时最后一次迭代的模型变化量在数学上提出确定印模深度的方法.原则是将反演结果中模型变化量非0的区域替换成均匀半空间，与原结果相组合，作为下一次重构的初始模型.当迭代重构结束时模型变化量为0或者两次迭代重构结束时的模型变化量的差值小于阈值时迭代重构结束.
(4) 采用本文的方法对地面MMR数据进行三维反演，在提升异常体的深度分辨能力的同时兼顾浅部的信息.异常体的电阻率随着迭代重构的次数增加而更接近反演模型的电阻率值.
本文中采用的是均匀半空间作为辅模型，如存在先验信息，也可以将其作为辅模型进行迭代重构.文中每次重构中采用的是非线性共轭梯度方法，但是本文的方法也可以与其他反演方法相结合.本文选用磁场单分量对异常体模型进行了反演研究，多源多分量的反演研究是未来工作的方向.因为磁电阻率法在我国相对应用较少，实测资料的采集与反演也是将来工作的重点.基于静磁场与重力场在数学形式上的相似性 (Chen et al., 2002)，在重力数据反演中也可以借鉴本文方法提高异常体的深度分辨能力.
附录A 磁场与电导率之间的关系推导我们分别从静电场麦克斯韦方程组和频率域麦克斯韦方程组两方面出发推导磁场与电导率之间的关系.
静电场 (频率为0) 的麦克斯韦方程组为
|
(A1a) |
|
(A1b) |
|
(A1c) |
其中E是电场强度，H是磁场强度，Jcs是电流密度，σ是电导率, μ是磁导率.根据式 (A1a)，存在一个电场的标量势ϕ，使得：
|
(A2) |
对式 (A1b) 取散度，就可以得到静电场的泊松方程为
|
(A3) |
将式 (A2) 带入到式 (A1b) 可以就可以得到磁场与电导率之间的关系为
|
(A4) |
如果从频率域麦克斯韦方程组出发，引入阻抗率
|
(A5a) |
|
(A5b) |
其中
|
(A6) |
其中F和A分别是引入的电场和磁场的矢量势函数，电导率σ包含在
感谢中国地震局地质研究所的陈小斌研究员对于本文提出的建议，感谢审稿专家的修改意见.
| Acosta J E, Worthington M H. 1983. A borehole magnetometric resistivity experiment. Geophysical Prospecting, 31(5): 800-809. DOI:10.1111/j.1365-2478.1983.tb01087.x | |
| Alumbaugh D L, Newman G A, Prevost L, et al. 1996. Three-dimensional wideband electromagnetic modeling on massively parallel computers. Radio Science, 31(1): 1-23. DOI:10.1029/95RS02815 | |
| Asten M W. 1988. The down-hole magnetometric resistivity (DHMMR) method. Exploration Geophysics, 19(2): 12-16. DOI:10.1071/EG988012 | |
| Aster R C, Borchers B, Thurber C H. Parameter Estimation and Inverse Problems. Amsterdam: Elsevier Academic Press, 2005. | |
| Bishop J, Carroll N, Asten M, et al. 1997. Finding sphalerite at Broken Hill with drillhole magnetometric resistivity. Exploration Geophysics, 28(2): 6-10. DOI:10.1071/EG997006 | |
| Bouchedda A, Giroux B, Allard M. 2015. Down-hole Magnetometric Resistivity inversion for zinc and lead lenses localization at Tobermalug, County Limerick, Irland.//2015 SEG New Orleans Annual Meeting. New Orleans, Louisiana: SEG, 2007—2011, doi: 10.1190/segam2015-5849363.1. | |
| Cattach M K, Stanley J M, Lee S J, et al. 1993. Sub-Audio Magnetics (SAM)—a high resolution technique for simultaneously mapping electrical and magnetic properties. Exploration Geophysics, 24(4): 387-400. DOI:10.1071/EG993387 | |
| Chaladgarn T, Yooyuanyong S. 2013. Mathematical model of magnetometric resistivity sounding for a conductive host with a bulge overburden. Applied Mathematical Sciences, 7(5-8): 335-348. | |
| Chen J, Oldenburg D W. 2003. 3-D inversion of magnetic induced polarization data. Exploration Geophysics, 37(3): 245-253. DOI:10.1071/EG06245 | |
| Chen J, Oldenburg D W. 2004. Magnetic and electric fields of direct currents in a layered earth (short note). Exploration Geophysics, 35(2): 157-163. DOI:10.1071/EG04157 | |
| Chen J P, Haber E, Oldenburg D W. 2002. Three-dimensional numerical modelling and inversion of magnetometric resistivity data. Geophysical Journal International, 149(3): 679-697. DOI:10.1046/j.1365-246X.2002.01688.x | |
| Deng J W. 2005. Forward Theory Study of Magnetometric Resistivity and Induced Polarization Methods[Ph. D. thesis](in Chinese). Beijing: China University of Geosciences. | |
| Edwards R N. 1974. The magnetometric resistivity method and its application to the mapping of a fault. Can. J. Earth Sci., 11(8): 1136-1156. DOI:10.1139/e74-108 | |
| Edwards R N, Howell E C. 1976. A field test of the magnetometric resistivity (MMR) method. Geophysics, 41(6): 1170-1183. DOI:10.1190/1.2035911 | |
| Edwards R N, Lee H, Nabighian M N. 1978. On the theory of magnetometric resistivity (MMR) methods. Geophysics, 43(6): 1176-1203. DOI:10.1190/1.1440887 | |
| Edwards R N, Nabighian M N. 1991. The magnetometric resistivity method.//Electromagnetic Methods in Applied Geophysics: Investigations in Geophysics. Society of Exploration Geophysicists, 47-104. | |
| Egbert G D, Kelbert A. 2012. Computational recipes for electromagnetic inverse problems. Geophysical Journal International, 189(1): 251-267. DOI:10.1111/j.1365-246X.2011.05347.x | |
| Fathianpour N, Cattach M K. 1995. Analytical solutions for the total field magnetometric resistivity (TFMMR) technique. Exploration Geophysics, 26(3): 158-166. DOI:10.1071/EG995158 | |
| Fathianpour N, Heinson G, White A. 2005a. The total field magnetometric resistivity (TFMMR) method: Part Ⅰ: theory and 2. 5D forward modelling. Exploration Geophysics, 36(2): 181-188. DOI:10.1071/EG05181 | |
| Fathianpour N, Heinson G, White A. 2005b. The total field magnetometric resistivity (TFMMR) method: Part Ⅱ: 2D resistivity inversion of data from the Flying Doctor Deposit, Broken Hill, Australia. Exploration Geophysics, 36(2): 189-197. DOI:10.1071/EG05189 | |
| Fu L K. The Principle of Magnetometric Resistivity Method. Beijing: Geological Publishing House, 1984. | |
| Godber K E, Bishop J R. 2007. DHMMR: Coming of age.//Milkereit B. Exploration in the new millennium: Proceeding of the 5th Decennial International Conference on Mineral Exploration. Keystone: AGU-SEG, 1119-1123. | |
| Howland-Rose A W, Linford G, Pitcher D H, et al. 1980a. Some recent magnetic induced-polarization developments—Part Ⅰ: Theory. Geophysics, 45(1): 37-54. DOI:10.1190/1.1441038 | |
| Howland-Rose A W, Linford G, Pitcher D H, et al. 1980b. Some recent magnetic induced-polarization developments-Part Ⅱ: Survey results. Geophysics, 45(1): 55-74. DOI:10.1190/1.1441040 | |
| Huston D L. 2006. Mineralization and regional alteration at the Mons Cupri stratiform Cu-Zn-Pb deposit, Pilbara Craton, Western Australia. Mineralium Deposita, 41(1): 17-32. DOI:10.1007/s00126-005-0036-4 | |
| Kofoed V O, Jessop M L, Wallace M J, et al. 2011. Unique applications of MMR to track preferential groundwater flow paths in dams, mines, environmental sites, and leach fields. Leading Edge, 30(2): 192-204. DOI:10.1190/1.3555330 | |
| Kong Z Z, Shan K S, Wu Y, et al. 2015. The improvement and applications of the impression method in AMT data processing. Geophysical and Geochemical Exploration, 39(2): 416-420. | |
| LaBrecque D, Sharpe R, Casale D, et al. 2003. Combined electrical and magnetic resistivity tomography: synthetic model study and inverse modeling. Journal of Environmental & Engineering Geophysics, 8(4): 251-262. | |
| Lei D. 2010. Studies and applications of 2-D CSAMT modeling and inversion with a dipole source and topography. Chinese J. Geophys., 53(4): 982-993. DOI:10.3969/j.issn.0001-5733.2010.04.023 | |
| Li J H, Lin P R, Guo P. 2010. A trial study of magnetic induced polarization. Seismology and Geology, 32(3): 492-499. DOI:10.3969/j.issn.0253-4967.2010.03.016 | |
| Li X, Bai D H, Yan Y L, et al. 2016. Three-dimensional inversion of magnetotelluric resistivity model with galvanic distortion. Chinese J. Geophys., 59(6): 2302-2315. DOI:10.6038/cjg20160632 | |
| Li Y G, Oldenburg D W. 1996. 3-D inversion of magnetic data. Geophysics, 61(2): 394-408. DOI:10.1190/1.1443968 | |
| Li Y G, Oldenburg D W. 1998. 3-D inversion of gravity data. Geophysics, 63(1): 109-119. DOI:10.1190/1.1444302 | |
| Li Y G, Oldenburg D W. 2000. Joint inversion of surface and three-component borehole magnetic data. Geophysics, 65(2): 540-552. DOI:10.1190/1.1444749 | |
| Lin C H, Tan H D, Shu Q, et al. 2012. Three-dimensional conjugate gradient inversion of CSAMT data. Chinese J. Geophys., 55(11): 3829-3838. DOI:10.6038/j.issn.0001-5733.2012.11.030 | |
| Liu Y H, Yin C C. 2013. 3D inversion for frequency-domain HEM data. Chinese J. Geophys., 56(12): 4278-4287. DOI:10.6038/cjg20131230 | |
| Maxwell M, Eso R, Oldenburg D. 2013. Using 2D and 3D electrical resistivity and magnetometric resistivity techniques for investigating dam and dike soil conditions for leak detection — field examples and forward modelling.//Smith B D ed. Symposium on the Application of Geophysics to Engineering and Environmental Problems 2013. Denver: Environment and Engineering Geophysical Society, 711, doi: 10.4133/sageep2013-264.1. | |
| Meju M A. 2002. Geoelectromagnetic exploration for natural resources: models, case studies and challenges. Surveys in Geophysics, 23(2-3): 133-206. | |
| Nabighian M N. 1987. Electromagnetic Methods in Applied Geophysics. Vol. 1. Tulsa: SEG. | |
| Newman G A, Alumbaugh D L. 2000. Three-dimensional magnetotelluric inversion using non-linear conjugate gradients. Geophysical Journal International, 140(2): 410-424. DOI:10.1046/j.1365-246x.2000.00007.x | |
| Nobes D C, Law L K, Edwards R N. 1986. The determination of resistivity and porosity of the sediment and fractured basalt layers near the Juan de Fuca Ridge. Geophysical Journal International, 86(2): 289-317. DOI:10.1111/j.1365-246X.1986.tb03830.x | |
| Seama N, Tada N, Goto T N, et al. 2013. A continuously towed vertical bipole source for marine magnetometric resistivity surveying. Earth, Planets & Space, 65(8): 883-891. DOI:10.5047/eps.2013.03.007 | |
| Seigel H O. 1974. The magnetic induced polarization (MIP) method. Geophysics, 39(3): 321-339. DOI:10.1190/1.1440431 | |
| Szarka L. 1987. Geophysical mapping by stationary electric and magnetic field components: a combination of potential gradient mapping and magnetometric resistivity (MMR) methods. Geophysical Prospecting, 35(4): 424-444. DOI:10.1111/gpr.1987.35.issue-4 | |
| Tang J T, Ren Z Y, Hua X R. 2007. The forward modeling and inversion in geophysical electromagnetic field. Progress in Geophysics, 22(4): 1181-1194. | |
| Weng A H, Li D J, Li Y B, et al. 2015. Selection of parameter types in Controlled Source Electromagnetic Method. Chinese J. Geophys., 58(2): 697-708. DOI:10.6038/cjg20150230 | |
| Weng A H, Liu Y H, Gao L J, et al. 2009. Magnetoelectric resistivity response characteristics for natural gas hydrate. OGP, 44(Supplement1): 158-161. | |
| Weng A H, Liu Y H, Jia D Y, at al. 2012. Three-dimensional controlled source electromagnetic inversion using non-linear conjugate gradients. Chinese J. Geophys., 55(10): 3506-3515. DOI:10.6038/j.issn.0001-5733.2012.10.034 | |
| Weng A H, Liu Y H, Yin C C, et al. 2016. Singularity-free Green's function for EM sources embedded in a stratified medium. Applied Geophysics, 13(1): 25-36. DOI:10.1007/s11770-016-0549-x | |
| Yang C H, Tseng H W. 1992. Topographic responses in magnetometric resistivity modeling. Geophysics, 57(11): 1409-1418. DOI:10.1190/1.1443208 | |
| Yang J. 2005. Geo-electrical responses associated with hydrothermal fluid circulation in oceanic crust: feasibility of magnetometric and electrical resistivity methods in mapping off-axis convection cells. Exploration Geophysics, 36(3): 281-286. DOI:10.1071/EG05281 | |
| Ye T, Chen X B, Yan L J. 2013. Refined techniques for data processing and two-dimensional inversion in magnetotelluric (Ⅲ): using the Impressing Method to construct starting model of 2D magnetotelluric inversion. Chinese J. Geophys., 56(10): 3596-3606. DOI:10.6038/cjg20131034 | |
| Yee K S. 1966. Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media. IEEE Trans. Ant. Prop., 14(3): 302-307. DOI:10.1109/TAP.1966.1138693 | |
| Zhdanov M S. Geophysical inverse theory and regularization problems. Amsterdam: Elsevier Academic Press, 2002. | |
| Zhu K G, Yang J W. 2008. Time-dependent magnetometric resistivity anomalies of groundwater contamination: Synthetic results from computational hydro-geophysical modeling. Applied Geophysics, 5(4): 322-330. DOI:10.1007/s11770-008-0041-3 | |
| 邓靖武. 2005. 磁电法正演理论研究[博士论文]. 北京: 中国地质大学. | |
| 傅良魁. 磁电勘探法原理. 北京: 地质出版社, 1984. | |
| 孔志召, 山科社, 吴勇, 等. 2015. 印模法在AMT数据处理中的改进与应用. 物探与化探, 39(2): 416–420. DOI:10.11720/wtyht.2015.2.34 | |
| 雷达. 2010. 起伏地形下CSAMT二维正反演研究与应用. 地球物理学报, 53(4): 982–993. DOI:10.3969/j.issn.0001-5733.2010.04.023 | |
| 李建华, 林品荣, 郭鹏. 2010. 磁激电方法技术试验研究. 地震地质, 32(3): 492–499. DOI:10.3969/j.issn.0253-4967.2010.03.016 | |
| 李鑫, 白登海, 闫永利, 等. 2016. 考虑电流型畸变的大地电磁三维反演. 地球物理学报, 59(6): 2302–2315. DOI:10.6038/cjg20160632 | |
| 林昌洪, 谭捍东, 舒晴, 等. 2012. 可控源音频大地电磁三维共轭梯度反演研究. 地球物理学报, 55(11): 3829–3838. DOI:10.6038/j.issn.0001-5733.2012.11.030 | |
| 刘云鹤, 殷长春. 2013. 三维频率域航空电磁反演研究. 地球物理学报, 56(12): 4278–4287. DOI:10.6038/cjg20131230 | |
| 汤井田, 任政勇, 化希瑞. 2007. 地球物理学中的电磁场正演与反演. 地球物理学进展, 22(4): 1181–1194. | |
| 翁爱华, 李大俊, 李亚彬, 等. 2015. 数据类型对三维地面可控源电磁勘探效果的影响. 地球物理学报, 58(2): 697–708. DOI:10.6038/cjg20150230 | |
| 翁爱华, 刘云鹤, 高丽娟, 等. 2009. 天然气水合物的磁电阻率响应特征. 石油地球物理勘探, 44(增刊1): 158–161. | |
| 翁爱华, 刘云鹤, 贾定宇, 等. 2012. 地面可控源频率测深三维非线性共轭梯度反演. 地球物理学报, 55(10): 3506–3515. DOI:10.6038/j.issn.0001-5733.2012.10.034 | |
| 叶涛, 陈小斌, 严良俊. 2013. 大地电磁资料精细处理和二维反演解释技术研究 (三)——构建二维反演初始模型的印模法. 地球物理学报, 56(10): 3596–3606. DOI:10.6038/cjg20131034 | |