大地电磁测深利用“丰富”的天然场源可以探测地下几百米到几十公里深的电性结构，在国内外都得到大力的发展和广泛的应用.反演是地球物理的核心内容，目前，有关大地电磁测深的反演方法从最早的一维地电结构假设已过渡到二维反演甚至三维反演阶段.从反演方法本身可将其分为线性反演和非线性反演两类.就实用性而言，以线性反演方法为主的算法应用较广，而非线性反演算法由于普遍存在计算效率低的问题，目前尚处于研究阶段.因此，本文针对线性反演算法就其理论和反演效果进行比较.由于地下电性结构的复杂性和数据的不准确性(数据本身带有误差和数据的有限性)以及算法本身的限制，任何一种单纯的反演算法都不可能真正再现地下电性结构，对于前两者，反演本身是无能为力的，只能通过算法本身的约束来实现最大可能的近似.国内外学者经过多年的研究和实践，也由此提出了许多优秀的算法.

在一维反演方面：Bostic作为一种半定量化反演方法由于能较好地反映地电断面的基本特征，被广泛应用于初步解释和作为其它反演方法的初始模型；上世纪80年代中国地质大学王家映教授等(1985)借助于常规地震勘探处理中的某些基本理论提出大地电磁拟地震解释法；90年代徐世浙院士和刘斌(1995)提出一维连续介质反演的曲线对比法.其它一些一维反演算法如：马夸特法、广义逆矩阵法等也都得到了很好应用.但由于地下电性结构的复杂性，很少存在真正的一维情况，目前一维反演方法常用于一些辅助性解释.在二维反演方面：Constable等(1987)基于最小可能构造下寻找拟合数据的光滑模型开发了OCCAM反演；Smith和Booker(1991)为提高反演效率，基于电性主要以垂向变化为主开发了RRI反演；Rodi和Mackie(2001)抛弃了线性反演的框架，提出了非线性共轭梯度法(其本质任为线性反演)；Smith等人(1999)基于如何更清晰地突出电性体界面提出陡边界反演(Sharp boundary inversion)，Myeong-Jong Yi等(2003)提出Smoothness constrained least-squares反演.国内学者戴世坤(1997)提出了MT二维和三维连续介质快速反演；严良俊和胡文宝(2004)提出MT资料二次函数逼近非线性反演方法；胡祖志等(2006)提出MT非线性共轭梯度拟三维反演； 鉴于最小构造和最大平滑准则思想的反演方法所得地电断面图像比较模糊，同济大学刘小军等(2007)根据Protniaguine提出的聚焦反演思想实现了MT聚焦反演法(Focus).本文就普遍应用的线性反演算法从算法理论本身到反演效果进行对比分析，并给出实际工作中的大地电磁剖面反演解释实例. 1 反演理论

1.1 反演目标函数

反演问题就是求解目标函数的最优化问题，由于地球物理反演固有的非唯一性与不稳定性，直接拟合数据所求得的模型往往有无穷多个，吉洪诺夫提出正则化目标函数，大地电磁正则化反演目标函数构制一般分为两部分，即：数据拟合项与模型约束项:

令：F为正演算子，d为观测数据向量，σ _j为第j个数据的标准差，m ₀为参考模型(包含先验信息)，W _m为模型约束矩阵，W_d为利用实测数据标准差进行规一化的矩阵.

由于实测数据不可避免地存在误差，因此实际反演中并不是拟合差越小越好，而应保持在由数据误差水平决定的一定的范围内，此时：

对于模型结构项，由于解得非唯一性，为了求得最优解，通常以模型范数的大小进行衡量，一般认为“最小模型”保留了真实地电模型的最基本特征，去除了冗余的复杂构造，即：

不同的模型结构项代表不同的反演算法，其相应的反演效果也存在一定的差别.

由于MT反演是非线性问题，因此在线性反演时，将非线性问题在初始模型 m _k附近线性化，并保留一次项，忽略高次项，即：

将(2)(3)(4)式代入(1)式，并根据函数的极值条件，求目标函数 φ(m，μ)关于模型向量m或m ^T的偏导数，并令其等于零，可得不同反演算法第k次迭代模型的求解公式.假设实测数据误差相互独立，则：W_d=diag{1/σ₁，1/σ₂，…，1/σ_N}.不同的反演算法其模型约束项 W _m或模型结构项φ_mod的表现形式不同.方便起见，令W_d=W，并定义R为模型光滑矩阵，C为模型聚焦矩阵，S为尖锐边界矩阵，其具体的定义形式及模型迭代求解公式见下表 1.

表 1(Table 1)

表 1 MT二维线性反演算法模型约束及其模型求解公式 Table 1 The structural penalty of MT 2D linear inversion algorithms

表 1 MT二维线性反演算法模型约束及其模型求解公式 Table 1 The structural penalty of MT 2D linear inversion algorithms

快速松弛反演(RRI)与表 1所列二维情况下反演方法不同的是：为了避免计算雅克比矩阵导致计算效率降低的问题，它通过解与一维相近的反演问题，来计算在每个测点位置下方的电阻率扰动，把二维问题转化为一系列一维反演.其反演迭代过程为(谭捍东等，2003)：给定初始模型，正演计算得到模型响应及需要的场值；对每个测点计算Frechet偏导数，并进行单点反演得到每个测点下的地电模型；内插形成新的模型；再正演得到新模型的响应和场值，直到满足终止条件.

OCCAM反演是利用模型参数梯度的最小范数求在一定拟合误差下的纵向(一维情况下)或纵向和横向(二维)的最光滑解，属于最平缓模型约束反演，但不利于突出电性体边界段.Focus反演和Sharp boundary反演的模型目标函数利用最小梯度支撑泛函，相比OCCAM反演能获得较清晰的电性界面图像，具有突出电性体边界的优点，但Sharp boundary反演由于把模型表示为有限数量的均匀电导率层(这些层的深度在横向上是变化的，反演的参数是各层的电导率以及所有测点下各层的深度)，在块体建模方面存在缺陷(张罗磊，2009).Smoothness constrained least-squares反演与马夸特反演的迭代求解公式相似，即：模型结构项的构制都基于模型修改量 Δ=m _k+1－ m k来进行的，而非模型本身；所不同的是马夸特法利用一个单位矩阵来改善系数矩阵的条件，使得病态或奇异线性方程组的解趋于稳定，而Smoothness constrained least-squares将单位矩阵替换为一个模型光滑矩阵的内积，这样不仅改善了系数矩阵的条件数，也使得求解模型光滑.

为了便于比较不同算法的反演效果，给出一个六层的一维地电断面模型(吴小平，1998)具体模型参数见表 2，并分别基于Focus、Smoothness constrained least-squares和OCCAM 进行反演.从图 1可以看出：无论Focus、Smoothness constrained least-squares还是OCCAM反演，它们对低阻层的反映相似，与真实低阻层接近；但在对高阻层的反映方面，OCCAM由于求模型的光滑解，使得对高阻反映不足，只是在趋势上与真实模型一致，对于复杂的二维/三维模型，有时甚至会使得高阻体完全淹没；Smoothness constrained least-squares反演则不同，由于它即保持了马夸特反演的优势(即对模型的修改量进行约束限制，而非模型本身)，同时还尽量使模型修改量光滑，一般而言由于初始模型为半空间或无剧烈的地电变化，因此求解结果与真实地电断面模型更为接近，但通过多个理论模型的反演结果，相比OCCAM反演它有时会出现冗余的构造；Focus反演由于利用最小梯度支撑泛函约束模型，反演结果对高阻层(第四层)出现了大于真实模型一个数量级的电阻率值.

表 2(Table 2)

表 2 一维地电断面模型参数 Table 2 The 1D model parameters of six lays 1 2 3 4 5 6 Resistivity(Ωm) 100 1000 10 1000 10 1000 Thickness(m) 100 100 50 400 150

表 2 一维地电断面模型参数 Table 2 The 1D model parameters of six lays

图 1

Fig. 1

图 1 一维地电断面不同反演方法结果 (a)原始视电阻率相位曲线图；(b)反演结果图；其中：Focus Rms=0.059；Smoothness constrained least-squares Rms=0.048；OCCAM Rms=0.061.Fig. 1 The inversion results of 1D model with six lays (a)Apparent resistivity and phase curves;(b)Resistivity curves of inversion results; Focus Rms=0.059；Smoothness constrained least-squares Rms=0.048；OCCAM Rms=0.061.

我们曾针对多个理论模型的数据进行了反演效果的对比(包括Sasaki模型、COMMMEM2D-3、COMMEMI2D-4模型等)，限于篇幅，这里给出一个针对金属矿产勘查较典型的一个二维地电模型：其地表自西向东分别为出露高阻岩体(中基性或中酸性)；背景为含粉砂质蚀变岩、大理岩化碎屑岩等中低阻体；东部为低阻残坡积物；且在剖面中部有一西倾含水破碎带(断裂)，剖面东向中深部为高阻侵入岩体；具体模型图见图 2所示.正演模拟计算时：频率取1000 Hz到1 Hz 共20个频点，测点点距为0.25 km，共20个测点，并将正演数据加入20%的随机误差(考虑到实际工作中高频数据质量总体大于低频这一现实，在数据加入噪声是也保持这一特点).为了便于比较反演计算效果，初始模型均为50 Ωm的半空间，且模型网格剖分一致(地表网格厚度为5 m，并以1:1.1的厚度递增).

图 2

Fig. 2

图 2 加入20%的随机相对误差反演结果图 (a)RRI反演；(b)Smoothness constrained least-squares反演；(c)OCCAM反演；(d)Focus反演.Fig. 2 The inversion results with adding 20% r and om noise to data (a)RRI inversion result;(b)Smoothness constrained least-squares inversion result;(c)OCCAM inversion result;(d)Focus inversion result; The initial model is show in a.

由图 2可以看出，四种反演方法(RRI、Smoothness constrained least-squares、OCCAM、Focus)都能较好还原理论模型的轮廓，其共同特征都为对低阻体的反应灵敏度大于高阻体(这也是MT方法本身的特点)，在电性异常的边界部位，低阻体有所扩大，对高阻体的反应能力显的不足.所不同的是Smoothness constrained least-squares反演结果的电性体边界部位更为清晰，且对高阻体的反映更突出一些，但同时也出现多余的构造；而OCCAM和Focus则表现得较好；RRI方法对倾斜低阻体反映为近似直立的低阻体，它对剖面东向中深部高阻反映则与OCCAM、Foucs类似都表现不足.

不同反演算法由于目标函数构制中模型约束机制不同，其反演结果存在一定的差异.总体而言，OCCAM反演结果平滑，能体现出地下电性体的基本特征，但有可能会产生模糊的地电图像，对于精细勘探不利于解释；Focus和Sharp boundary反演能够较清晰地突地下电性体的边界段，但Sharp boundary存在块体建模方面的缺陷，而纯粹的Focus反演电性变化剧烈也不利于定量解释；Smoothness constrained least-squares反演由于对模型修改量进行光滑限制而非模型本身，能够较好地反映地电模型，但同时也会产生冗余的构造；RRI反演效率高，但不易控制横向光滑因子，有时会将倾斜的地质体反映为直立的地质体.因此，在实测数据反演时，须考虑以上各种反演方法的优缺点，根据具体目标任务有选择性地制定反演方案或结合几种反演结果综合解译. 3 实际反演效果

某地锌多金属矿床主要赋存于下侏罗统日当组第四岩性段含碳钙质板岩夹褐黄色钙质砂岩中，矿体的产出严格受控于构造破碎带，呈层状、似层状或透镜状产出，倾向265°~295°，倾角40°~70°.区内地表覆盖严重，为直接根据地表矿化和蚀变找矿带来困难，且石墨化(碳化)围岩造成大面积的区内低阻高极化，也给激电直接找矿带来严重的地质干扰.但根据地表出露断裂带以及钻孔揭示的断裂，其含碳量相对围岩明显减少，具有石英脉填充、不含水等特征.物性测试结果显示：相对围岩，矿化(矿)为中高阻、中高极化(3%~5%)特征，因此物探工作的主要兴趣为寻找相对的中高阻体，且其极化率约为4%左右.众所周知，音频大地电磁测深设备轻便，探测深度可达1~2公里甚至更深，但对低阻体的反映灵敏度大于高阻体，因此在反演方法上如果利用传统的模型光滑约束反演，极有可能使得本就对高阻反映差的特征变得雪上加霜.通过前面的论述和试验，我们发现Smoothness constrained least-squares反演在对高阻体的反映方面相比直接对模型进行光滑约束的反演有所改善，反演结果稳定，因此在反演方案上以传统的OCCAM反演确定大致地电断面，利用Smoothness constrained least-squares反演来突出高阻体，两种反演的结果同时又互相印证.

图 3给出40线中段AMT测深以及激电测深反演结果，由AMT资料的RRI反演结果和Smoothness constrained least-squares反演结果可以看出，两种反演结果的地电断面形态大致一致，浅地表结果与都与地质资料对应较好，但在一些细节方面存在差异，如AMT11号点到AMT18号点间的电性结构，前者反映为一个低阻体由西向东逐渐过渡到中高阻体，而后者则反映为两个紧密相邻的低阻体和高阻体，也就是说后者更容易突出高阻体的边界.两种结果的AMT 22号测点下方都为大块低阻体，据以往面积性激电测量结果以及激电测深结果均显示该段为高极化特征(极化率大于8%)，因此推断为碳质板岩所致，ZK4001钻孔结果也验证了这一推断；AMT 18-19号测点下方等值线变化密集，为低阻到中高阻的梯度带，激电测深反映为中低极化率，据遥感解译区内存在的南北向断裂经过18号测点，部分地质露头也验证了断裂的存在，综合旁侧AMT测量结果以及化探异常特征，排除高阻低极化体，推测AMT 18号测点下方赋存西倾含矿断裂，根据设计位置在ZK4002孔深约200米处见6米厚铅锌矿体.

图 3

Fig. 3

图 3 某地AMT实测资料反演解释图 (a)测点地质底图；(b)等值线为RRI反演；(c)等值线为Smoothness constraine least-squares反演；其中颜色填充底图为激电测量反演结果.Fig. 3 The example of inversion results with observation data (a)geological map of observation posts； b、c：The fill map is polarization sounding inversion results，and the contours is MT 2D inversion results;(b)RRI inversion;(c)Smoothness constrained least-squares inversion.

基于正则化目标函数的大地电磁反演同时考虑到数据拟合和模型约束机制，反演结果稳定，其反演效果依据不同的模型约束机制对模型的光滑程度或电性边界段的反映不同，各有优缺点.总而言之，在同一误差水平条件下，反演效果依赖于其算法的数学表达式.在实际工作中，反演并非一蹴而就，应依据目标任务结合具体反演算法的优缺点制定反演方案或反演组合，不断反复求证.

随着计算机硬件水平的不断提高，近年来有关各向异性、非线性和三维反演成为反演发展的趋势.据有关文献报道在上世纪90年代，由于受当时的计算机硬件条件限制，OCCAM反演一条剖面需要几十个小时(MacⅡ)，但在计算机硬件技术高速发展的现在，一条剖面的反演往往仅需一个小时左右，反演效率提高了几十倍.因此，可以预见基于非线性反演方法、三维反演以及各向异性反演将会随着计算机计算效率的提高成为反演方向发展的主流.