徐义贤，蒋礼，岳瑞永，杨波，刘营，张秉政，孔亚丽

徐义贤，蒋礼，岳瑞永，杨波，刘营，张秉政，孔亚丽_大人地球物理科学

徐义贤^1a，b，蒋礼^1，2，岳瑞永^1a，3，杨波^1a，刘营^1a，张秉政^1a，孔亚丽^1a

1.中国地质大学，a.地球物理与空间信息学院；b.地质过程与矿产资源国家重点实验室，武汉，430074；

2.华北水利水电学院，郑州，450011；3.大连测控技术研究所，大连，116013

摘　要：中国地质大学地球物理学科在大地电磁法的理论研究与应用方面有40多年的历史，做出了许多重要贡献。本研究是对近10年来笔者所在团队在大地电磁观测与方法研究方面所取得的主要进展进行了综述，对每项研究的背景进行了简要介绍，对大地电磁法的发展前景进行了概述。

关键词：大地电磁法，远参考技术，多站叠加技术，电各向异性，大地电磁正反演

前　言

大地电磁法（MT）在中国地质大学已有40多年的研究历史，最早可以追溯到陈乐寿老师20世纪60年代末期在匈牙利留学期间的工作。从20世纪70年代开始，陈乐寿、吴广耀、王家映、魏文博、张胜业、吴农、戴联其、方胜、董浩斌、魏胜等先后加入到MT的教学和研究工作中。至2002年，中国地质大学在该领域贡献了多本广为采用的教材和广为引用的学术专著。在中国地质大学该领域众多的学术成果中，笔者认为至少有3个方面的贡献值得提及。王家映老师在MT的定量反演方面做出了重要贡献，尤其是在拟地震解释和非线性反演方法上的诸多努力^［1，2］，为推动MT在中国的推广应用发挥了重要作用。陈乐寿老师在MT资料的处理和解释方面有很多建树，尤其是他领导了第一期中美合作的INDEPTH项目，该项目首次获得了青藏高原喜马拉雅造山带的深部电性结构^［3］，为推动青藏高原的地球物理研究工作做出了贡献。魏文博老师领导了第二期至第四期的INDEPTH项目，第一次根据电性预测了青藏高原中下地壳广泛存在流体^［4］。同时，他领导了国家专项《深部探测技术与实验研究》中“大陆电磁参数标准网实验研究”项目（SinoProbe-1），为推动MT方法在中国大陆地球物理格架构建工作中的应用做出了贡献。

笔者于1993—1996年间追随王家映老师攻读博士学位，学习了MT方法的基本原理和反演方法，完成了《大地电磁多尺度反演》的博士学位论文。期间随当时的MT组张胜业老师等参加了柴达木盆地西缘的MT大剖面工作，熟悉了MT的野外工作方法技术，学习了MT的数据处理方法。2002年以来，虽然主要精力已不在MT方法的研究方面，但仍情被所系而无法放弃，其间借助于承担的1项国家自然科学基金项目、国家专项SinoProbe-01-03课题及中石化和国土资源部先后资助的“塔里木盆地深部流体探测”专题，指导2名博士和5名硕士研究生开展了MT野外观测方法技术、电各向异性介质对MT的影响及三维MT正反演方法等研究。在庆祝建校60周年之际，笔者与在读或已工作的从事过MT研究工作的研究生们一起，将2002年以来的研究成果或摘编、或改写、或新著述，结成此文，献给我们深爱的中国地质大学，向培育了我们的恩师们汇报。

1　MT野外观测方法技术

除了人文噪声，磁场信号中的噪声来源较为简单，磁棒的方位误差和地震动是两种很难避免的磁场噪声来源^［5］。地震动引起的噪声主要集中在MT信号的“死带”^［6］，因而使得任何测点估计的MT“死带”视电阻率和相位曲线均出现跳跃。大量文献资料表明，当远参考站设置在20公里至几百公里范围内时，远参考技术可以极大地改善MT“死带”的信号质量。笔者认为一种可能的解释是常时微动的主要周期为1～3s，因此在通常的远参考距离上，地震动噪声是不相关的。磁棒的方位误差引起的噪声影响程度可以和近地表电性不均匀体所造成的畸变效应相当^［6］，特别是短磁棒的采用可能会加大磁棒方位误差，因此必须在野外磁棒布设时加倍小心。裸露的长金属管线可以产生很强的非平面波磁场噪声^［7］，但可以在选点时尽量加以避免。一般地，测点需要远离裸露的长金属管线至少500m以上。其他类型的人文磁场噪声在相关的规范中均已提及并作出了严格的技术处理规定。

电场信号中的噪声来源相当复杂，除了人文噪声外，电极对的方位误差、地震动、极差的稳定性及温度效应等是需要考虑的几个因素。相对于磁棒，由于测量电极对（偶极子）的间距较大，因此方位误差较小。地震动对磁场测量信号的影响是直接的，而对电场信号的影响则可能主要来自于孔隙介质中产生的“动电效应”，这与电极对的高差及地下水的流动方向也产生了联系，因而在理论和实践上都是一个难以准确预测的误差来源。极差的稳定性极大地影响长周期MT观测数据的质量。控制实验表明：AgCl电极的灵敏度最高，而PbCl₂电极的长期稳定性较好^［8］。除电极选择之外，接地条件包括温度效应是野外最需要认真处理的问题。在我国西部干旱地区进行长周期MT测量时，保持接地条件的长期一致性是很困难的，需要采取多重措施，如挖掘很深的电极坑、底部放置湿润的泥土、电极包裹盐水浸湿的棉絮、采取措施防止电极周边温度变化太快、以淋滤方式长期加入盐水至电极坑等。

提高大地电磁（MT）观测资料的信噪比是获得高质量视电阻率和相位曲线的关键。传统上，除了选择电磁干扰小的观测点和采用抗噪能力强的数据处理方法^［9，10］外，延长观测时间和同步设置远参考站^{［11～13］}是两项用于提高MT信噪比的主要技术。这两项技术无疑都以更多的人力和物力投入为代价，而且并不总是奏效，如强电磁干扰环境下延长观测时间并不能提高电磁场信号的信噪比，当参考站与观测点电磁场噪声相关时远参考技术也无能为力。

1.1　远参考技术

自从20世纪70年代末期Goubau等1978年提出远参考技术^［11］以来，国内外众多学者进行了理论和实验研究。Gamble^{［12，13］}等（1979a、1979b）通过设置4.8km的远参考站，在0.01～100s周期段获得了明显改善的MT视电阻率和相位曲线。Jones等^［14］报道了远参考站分别为30km和135km的Robust处理结果，在接近6 000s时，远参考站距离为135km的结果优于30km。Tarakura^［15］等和黄哲^［16］进行了类似的实验对比，验证了远参考站的合适距离与需要改善的MT视电阻率和相位的周期有关，高频段需要的距离相对小，低频段需要的距离相对大。陈清礼等^［17］的实验研究给出了有意思的结果：在中频段和低频段，近参考（2km）和远参考（1 000km）的处理结果基本一致，仅在高频段有差异。杨生等^［18］通过实际资料的处理对比，建议在保证噪音不相关的条件下尽量靠近实测点，并采用均匀半空间模型的分析结果认为参考距大于14倍勘探深度为宜。Shalivahan^［19］等在实验中将参考距分别设置为80m、115m、215km，3个参考距获得的结果均对视电阻率和相位曲线的连续性有明显改善，而其中215km参考距对MT“死带”的改善结果最佳。从前人的研究成果中可以肯定的是：参考距的选择与频带有关，与参考站的电磁环境有关（也就是与记录的信噪比有关），最远距离需满足MT信号的均匀场假设。

由于远参考技术建立在MT有效输入信号一致性（高相干）基础上，与参考站下方的地电结构无关，而输入信号的一致性与天然电磁场的场源空间结构有关。因此，参考站的选择除了电磁环境噪声方面的考虑之外，最重要的应该是测点与远参考点的输入信号（入射到地球表面的电磁场水平分量）具有几乎相同的频率含量。这样，选择与测点处地球磁场的水平分量在观测频带内高相干且电磁干扰小的地点设置远参考站就成为唯一的原则。但是，设置一个远参考站一般难以满足超长周期MT观测的需要，因为单独的参考站所观测的电磁场甚至仅仅磁场也难以在超宽频带上都与测点的有效输入信号一致。

在海南岛中部将等距离（40km）参考站分别设置在测点正南和正西方向，实验结果如图1所示。该实验表明南北向远参考优于东西向远参考，40km左右的远参考对改善大地电磁“死带”的质量有很大帮助，但对于1 000s以上的频带来说，带远参考的处理结果更差，这至少说明远参考的距离需要考虑改善的信号频带。图2为湖北荆门实验测点的结果，由于北向远参考点7天资料中因仪器故障只有5天可用，东向远参考7天资料全部可用。实际结果表明，不论是北向远参考还是东向远参考，联合估计的MT阻抗质量均有提高，而且测点与北向远参考点联合估计的质量要高于测点与东向远参考点联合估计的质量。

上述两个观测实例说明：①远参考技术对于提高MT阻抗的估计质量十分有效，但成本成倍增加；②在我国中部和东南部地区，经线方向远参考效果优于纬度向远参考；③远参考点的距离需要考虑信号频带。一般地，参考站应设置在需要改善的信号频带对应的趋肤深度相当的距离范围。

1.2　多站叠加技术

在提高MT阻抗估计质量的野外观测方法和技术中，远参考技术和延长观测时间的策略，均存在成本较高的缺陷。为此，我们提出并经过理论和实验证明了类似于反射地震中检波器组合的MT多站叠加技术的有效性^{［20，21］}。

1.2.1　多站叠加原理

假设M套观测仪器对大地电磁信号进行同步记录，x_i（k）为第i套仪器的实测信号，s（k）为有效信号（为讨论方便假设其功率为1），a_i为该仪器的有效增益，n_i（k）为该仪器所接收的零均值白噪声（其方差为σ²_i），且设定各路噪声存在相关性，n_i（k）和n_j（k）的互相关系数为r_ij。当记录时长为K时，多站观测系统所记录的时间序列可记为：

图1　不同方向等距离远参考实验结果（测点位于海南岛中部）

（a）测点与其同一纬度西向40km同步观测的远参考站联合估计的阻抗；（b）测点独立数据估计的MT阻抗；（c）测点与其同一经度南向40km同步观测的远参考站联合估计的阻抗。除远参考外，处理中的其他参数均相同

图2　不同方向等距离远参考实验（湖北荆门测点）

（a）测点与其同一纬度东向220km远参考站联合估计的阻抗；（b）测点独立数据估计的MT阻抗；c）测点与其同一经度北向220km远参考站联合估计的阻抗。除远参考外，处理中的其他参数均相同

该式可缩记为：

设大地电磁多站叠加的加权系数为W＝［w₁，…，w_M］^T，各台站观测数据的单站最大信噪比为SNR＝a/σ，那么通过多站叠加方法所能提高的最大增益为：

式中：Θ＝E ［NN^H］，表示求解叠加噪声的协方差矩阵^［22］。

假设M个观测台站同步记录数据，且每个台站的信噪比完全相等（即a_i/σ_i＝a/σ）；而且对于所有台站而言，任意两组台站所测得随机噪声的互相关系数也完全相等（即r_ij＝r）。依据式（3），此时M站叠加的最大信噪比增益为^［23］：

从式（4）中可以看出，多站叠加增益不仅与进行叠加的观测台站数目有关，而且与各台站间的噪声相关性有关。M个观测台站进行叠加所能获得的最大信噪比增益与成正比。

由于式（3）不要求各站噪声的互相关系数相等，因此更具一般性。而式（4）对于认识多站叠加对信噪比的改善作用是有帮助的，它与地震检波器组合中导出的公式^［24］相同。

1.2.2　多站叠加技术对近地表电性不均匀体引起的静位移的压制作用

多站叠加能够克服EMAP仅能压制沿测线方向电场畸变的弱点，对于由地表电性不均匀体引起的任意方向的电场异常有较好的压制作用。如图3所示，假设O点为标准测点，在以O点为圆心，r为半径的圆周上均匀布设多个观测台站用以进行多站叠加处理。为了论述方便，此处仅选择A、B、C3个台站的观测数据进行叠加，这3个观测台站呈等边三角形分布，且A站位于y轴上。假设半径为R的电性不均匀体的中心位于D点，此时D点距离标准测点O和辅助测点A、B、C的距离分别为：r_O、r_A、r_B和r_C。假设背景场E₀的方向为x轴的正方向，并设定x轴的正半轴与r_O的逆时针夹角为φ₀。

如果没有异常体存在，在地表任意位置的观测站都将测量到：E_x＝E₀，E_y＝0这种真实的背景场值。但由于电性异常体的存在，致使背景场发生畸变，此时地表的不同测站将由于方位和距离的不同，测量到不同的场值^［25］。定义不均匀体和地下半空间的电导率之比为K，

即：。

（1）假设P点位于不均匀体的内部，则该点处在x和y两个方向的电场值分别为：

分别去掉两个方向上的真实背景场值，可以得到这两个方向上的异常场值：

（2）假设P点位于不均匀体的外部，则该点处在x和y两个方向的电场值分别为：

同理，可以得到这两个方向上的异常场值：

式（7）和式（8）中：x_P、y_P是以D点为坐标原点的P点横、纵坐标值，r_P为P点距坐标原点的距离。

假设电性不均匀体的半径R为20m，辅助台站距离标准测站的距离r为100m，且电性不均匀体与地下半空间的电导率之比K等于10。分析的过程由不均匀体的中心点与标准测点重合开始，随着r_O的增加不均匀体逐渐远离标准测站而接近辅助测站。每当r_O确定之后，将不均匀体逆时针旋转360°，以便对其进行全方位的分析。为了作图美观，设定角度变量φ为r_O和y轴正半轴之间的夹角（图3），即：；为了作图方便，式（5）～式（8）全部除以背景场值E₀，即进行归一化处理。

图3　多站观测时均匀半空间含半球形电性不均匀体模型

（1）当电性不均匀体的中心位于O点，即：r_O＝0。

因为r》R，所以在x方向上，三站叠加的异常场值远小于标准点的异常场值。

在y方向上，ΔE_y^O＝0。又因ΔE_y^A＝0、ΔE_y^B＝-ΔE_y^C，所以三站叠加的异常场ΔE_y＝0。此时在y方向上，标准点与多站叠加的结果都未受到异常场的干扰。

（2）当O点依旧处于电性不均匀体的内部，即：0＜r_O＜20m。

将上述模型的参数值代入式（5）中，可以求出标准点的异常场值为ΔE_x^O＝0.75，而由式（7）可以求出多站叠加的异常场值，其最大值maxΔEx＜＜0.01，可参见图4。可见此时，在x方向上三站叠加的异常场值也远小于标准点的异常场值。

此时ΔE_y^O＝0，但ΔE_y≠0，尽管其异常场值非常微弱（图4）。所以在y方向上标准点的结果优于多站叠加的结果。

（3）当O点处于电性不均匀体的外部，即：r_O＞20m。

标准点与多站叠加的异常场值对比结果如图4所示，根据式（7）可以绘制出x方向上标准点的异常场值图［图4（a）］和多站叠加的异常场值图［图4（b）］；根据式（8）可以绘制出y方向上标准点的异常场值图［图4（c）］和多站叠加的异常场值图［图4（d）］。在图中横轴表示角度φ（单位为弧度），纵轴表示相对异常场值（已除以E₀），不同的线型代表着异常体距标准测点的不同距离。

对比图4（a）与图4（b）及图4（c）与图4（d）可以看出，就整体趋势而言，随着r_O的逐渐增大，不均匀体渐渐远离中心点，O点处的异常场值越来越小；同时不均匀体距辅助台站越来越近，多站叠加的异常场值逐渐增大。但总体而言，图4（b）和图4（d）的异常场值明显小于图4（a）和图4（c）的异常场值。可见，只要r_O≤50m时，三站叠加的异常场值整体上要小于中心点的异常场值。

图4　中心点与多站叠加观测到的由于电性不均匀体产生的附加异常场随方位的变化

1.2.3　多站叠加技术对横向传播的平面电磁波的压制效应

如图5所示，假设E₀为平面电磁波，箭头所示为平面波的传播方向；A、B为两个辅助测站，同步记录大地电磁数据，A、B间连线距离的长度为l；假设电磁波的传播方向与l的夹角为α，而波前面与l的夹角为β。

图5　两站叠加对于平面波的压制作用

根据图5，可以方便地计算出电磁波路程差为：

Δ＝l·cosα＝l·sinβ

将路程差直接转换为相位差：

式中：σ为电导率；μ为磁导率。

由相位差的知识可知：当Δφ为π的偶数倍时，叠加之后对平面波的压制效果最小；而当Δφ为π的奇数倍时，叠加之后对平面波的压制效果最大。

在野外实际的大地电磁多站叠加数据采集过程中，如果各辅助台站的位置、平面波的传播方向和地下介质的电性参数全部确定下来后，相位差就仅与频率的选择有关。依然以两站叠加为例，假设台站间距l为100m，波前面与l的夹角β为π/6（图5），磁导率μ等于μ₀为4π×10^-7 H/m，而电导率σ为0.1S/m。将以上参数代入式（9）可计算出相位差为：

从式（10）中可以看出，相位差与频率的平方根成正比。假定长周期大地电磁观测的频带范围从320Hz到10 000s，进行两站叠加的其他参数不变，利用式（10）可以计算出相位差介于0.18π至10^-4π之间。如果对于高导介质，电导率σ为1S/m，那么对于同样频段而言，相位差将介于0.57π至3.16×10^-4π之间。显而易见，多站叠加对于平面波的压制在低频段并无太大影响，但对于高频信号尤其在低阻介质的情况下却不可忽略。

根据上述理论和实验研究，我们认为MT多站叠加技术在增加成本较小的情况下可以有效提高MT阻抗的估计质量。带来的另外一个好处是：多站叠加类似于设置了一个二维空间滤波器，可以压制测点附近电性不均匀体产生的静位移，但对沿地下传播的平面电磁波的衰减却很小。

2　一维电各向异性介质对MT和CSAMT的影响

电阻率各向异性的可能起因包括岩石中某些矿物的优势分布、岩石颗粒的定向伸长、微裂隙或微结构的定向排列、断层、薄互层等。观测到地球介质的电阻率各向异性最早可追溯到20世纪30年代^［26］（Schlumberger，1934）。Hill^［27］（1972）在实验室测量变质岩样品时发现存在强烈的电各向异性，并且随着频率的降低而增强。Herwanger^［28］等（2004）对采自英格兰西南部Cornwall的细层状沉积岩样品进行了电阻率测量，发现电阻率最大到最小的变化范围为1 240～330Ω·m。在野外尺度的观测数据中，大量的研究者发现了岩石圈存在电各向异性^{［29～31］}（如Eisel、Haak，1999；Bahr，2000；Wannamaker，2005），含水层也存在电各向异性^{［32，33］}（如Darboux-Afouda、Louis，1989；Ritzi、Andolsek，1992）。Weidelt^［34］强调了在三维电磁正反演中研究电各向异性的重要性。Wannamaker^［31］全面回顾了描述电各向异性的理论，产生电各向异性的地质因素，研究电各向异性的意义。我们检查了CSAMT近场区传统定义在一维电各向异性介质中的不适用性，研究了一维大地电磁反演中用电阻率各向同性模型拟合电各向异性模型的误差。

2.1　描述电各向异性介质的一般模型及常用模型的一维MT阻抗特征

对于各向异性问题，可以把电导率看作一个张量。一般来说，它是一个3×3的调和且正定的矩阵：

将电导率张量旋转到3个主电导率方向上，参考图6：

其中R是单位旋转矩阵，

图6 一般各向异性介质坐标旋转示意图

（据Pek，2002^［35］，有修改）

对于实际地球物理问题来说，可以将电各向异性介质分为：横各向同性介质（TI介质）、方位各向异性介质、TTI介质以及一般各向异性介质。

2.1.1　含横各向同性介质的MT阻抗

对于一般沉积盆地，若每一层水平方向电阻率相等且与其在垂直方向上的电阻率不等，则称为电横向各向同性介质。对于一维TI介质，正演的ρ_xy和ρ_yx视电阻率曲线总是重合的，而ρ_xx和ρ_yy的值为零，说明大地电磁测深对一维TI介质的电各向异性是无法分辨的。根据大地电磁法原理，由于场源是垂直入射的平面电磁波，一维TI介质的对称轴与场源的入射方向一致，水平平面内不存在电性主轴，因此无论如何改变纵向与横向电导率的差异，都不能改变ρ_xy和ρ_yx视电阻率曲线总是重合的事实。

2.1.2　含方位各向异性介质的MT阻抗

近20余年来，人们利用横波窗内观测的地震波S波分裂资料，先后在全球许多地区发现在约20km深度以上的上地壳内比较普遍地存在地震波传播的各向异性现象，并将其物理机理归结为上部地壳岩石的膨胀各向异性（EDA）。上地壳内，特别是对于一些构造比较活动的地区，最大主应力和最小主应力都近似位于水平面内，形成了沿最大主应力方向优势排列的近于直立的裂隙，这些裂缝一般显示出方位各向异性性质^［36］。对于表1给出的模型，其第二层为方位各向异性介质，正演所得的视电阻率张量如图7所示。

表1　方位各向异性介质模型

图7　夹层为电阻率方位各向异性模型（表1）MT正演的视电阻率曲线

（a）ρ_xy和ρ_yx；（b）ρ_xx和ρ_yy

由图7（a）可见，ρ_xy和ρ_yx在曲线首支和尾支重合，与设置的模型第一层和第三层均为100Ω·m的各向同性层吻合。在曲线中间部分ρ_xy和ρ_yx发生了分离，这是模型第二层是方位各向异性层的响应。从图7（b）可以看出，ρ_xx和ρ_yy两条视电阻率曲线是重合的，而且它们的值非常小；当α²_S＝0时，即观测方向和电性主轴方向重合时，ρ_xx和ρ_yy的值为零，因此两者值的大小受α²_S的影响。

2.1.3　含TTI介质的MT阻抗

在一些沉积盆地，沉积地层由于受构造作用的改造，或者泻湖沉积、近海岸沉积等的作用，造成TI介质的对称轴不再垂直向下，而是倾斜的，称为电TTI介质。这种情况下，沿倾斜面方向的电阻率相同，而其与垂直于倾斜面方向的电阻率不同。对于表2给出的电性TTI模型，其正演曲线如图8所示。

表2　电性TTI介质模型

图8　表层为电阻率TTI介质模型（表2）MT正演视电阻率曲线

（a）ρ_xy和ρ_yx；（b）ρ_xx和ρ_yy

从图8（a）可以看出，ρ_xy和ρ_yx在曲线尾支基本重合，与设置的模型第三层为各向同性的情况是相对应的。曲线首支和中段的ρ_xy、ρ_yx发生了分离，这是模型第一层为倾斜各向异性层的响应。ρ_xy和ρ_yx首支的渐近线值并非等于100Ω·m，而是有所偏差，这主要是因为观测方向和电性主轴方向存在角度α¹、α¹所造成的。从图8（b）可以看出，ρ和ρ两条视电阻率曲线_SDxxyy是重合的，而且它们的值比较小，两者的值也同时受到α¹、α¹的影响。_SD

2.1.4　含一般各向异性介质的MT阻抗

在一些复杂地质构造背景下，可能出现任意的电阻率各向异性介质，比如TTI介质又发生了旋转运动，此时3个角度都会有值。对于表3给出的一般各向异性介质模型，其MT正演曲线如图9所示。

表3　一般电阻率各向异性介质模型

图9　夹层为电阻率一般各向异性介质模型（表3）MT正演视电阻率曲线

（a）ρ_xy和ρ_yx；（b）ρ_xx和ρ_yy

从图9（a）可以看出，ρ_xy和ρ_yx在曲线首支和尾支重合，与设置的模型第一层和第三层是各向同性的情况相符合。ρ_xy和ρ_yx在曲线中间部分发生了分离，这是因为模型第二层是电阻率一般各向异性的响应，决定于3个方向电阻率的值与α²、α²、α²的组合。从图9（b）可以看出，_SDLρ和ρ两条视电阻率曲线是重合的，而且它们的值非常小，值的大小也同时受α²、α²、α²的_xxyySDL影响。

通过分析上述3个不同电阻率各向异性模型的MT正演曲线结果，可以得到如下结论：①MT无法鉴别是否存在对称轴垂直于地面的TI介质；②除TI介质外，电阻率各向异性层的存在使得MT视电阻率张量中两非对角线元素在相应的频段产生明显偏离，而两对角线元素保持重合，这是识别是否存在电阻率各向异性层的有用标志；③虽然电阻率各向异性层的存在会改变MT视电阻率张量中两对角线元素，但一般情况下两对角线元素的值均很小。

2.2　一维方位电各向异性对CSAMT近场距离的影响

近场效应是CSAMT法必须要考虑的问题。在近场区域，场值以1/r²规律衰减的磁场变得饱和，即不再随频率和地下介质的电阻率的变化而变化；电场虽仍然是电阻率的函数但不随频率而变化，场值按照1/r³规律衰减，这里r指收发距。在此区域，勘探深度不再与频率有关而仅与收发距有关，以发射点为中心，这一区域的半径大约为趋肤深度的3～5倍^［37］（Zonge、Hughes，1991）。由于前人的结论是从电各向同性均匀半空间模型出发的，笔者从电方位各向异性均匀半空间模型出发研究了CSAMT法的近场区域大小^［38］。

选择坐标系：x轴指向构造走向并与最大电导率方向平行，y轴垂直于构造走向，z轴垂直向下。对于电方位各向异性介质，电导率张量有形式：

各向异性系数和平均电导率则分别定义为

σ_n。此时，CSAMT在远场区域的ρ_xy和ρ_yx与大地电磁TE模式的视电阻率ρ_TE和TM模式的视电阻率ρ_TM应该一致，当它们之间的相对差异超过5%时，定义为近场区域。由于存在两种极化模式，因而有两种判断标准，分别称为TE（η₁）准则和TM（η₂）准则：

图10　不同电方位各向异性系数情况下根据TE准则计算的近场区域分布

假定σ_t＝0.01s/m，λ分别为1、1.5、2。分别计算CSAMT的视电阻率^［39］和大地电磁的视电阻率^［40］并比较它们的相对误差，可以定量地画出近场区域的大小。根据TE准则给出的近场区域分布见图10。如果测线沿x轴方向，定义趋肤深度δ_x≈503.29，则当λ为1.5时的趋肤深度与x轴方向的近场距离FFD_x有定量关系^［38］：FFD_x＝5.665δ_x。这说明当存在电方位各向异性，且测量方向与最大电导率方向平行时，近场距离超出了传统的估计范围^［37］（Zonge、Hughes，1991）。

进一步假定测量总是沿x轴方向（也是最大电导率方向），根据TE准则和TM准则，可以获得不同的近场区域分布模式。笔者分别计算了λ从1～2变化范围内，FFD_x在两种准则下的变化曲线（图11）。显见对于TE模式，当方位各向异性系数大于1.5，最大的近场距离将逐渐超出传统的估计；而对于TM模式，几乎只要存在电方位各向异性，最大的近场距离就会超出传统的估计。因此，在进行CSAMT野外工作时，为了数据的有效性起见，最好预先用电测深方法进行不同方位的电阻率测量，以确定最大和最小的电导率方向及各向异性系数，从而确定最小的收发距。

图11　测线沿最大电导率方向时近场距离随方位各向异性系数的变化情况

（a）TE准则；（b）TM准则给出的结果

3　三维大地电磁模拟与反演

近年来，大地电磁数据的反演与解释技术已逐渐向三维发展。交错网格的有限差分由于其简洁、高效的优点，从早期到现今，一直都是三维大地电磁正演模拟的主流方法^{［41～43］}。在反演方面，共轭梯度法是最早用于解决三维大地电磁数据反演的方法^［41］、非线性共轭梯度法^［44］和数据域的Occam^{［45，46］}反演算法受到了广泛的关注，已逐步走向了实用化。

3.1　三维大地电磁正演模拟原理

在三维介质情况下，电场E满足：

将式（1）按电场的3个分量进行分离，则可得

上式中，k＝√iωμσ。这是3个互相耦合的方程组。采用图12所示的交错网格采样方案将式（16）～（18）离散，则可以形成一个大型的对称复线性方程组：K·E＝S

（19）

式（19）中：矩阵K为一复对称稀疏矩阵，每一行的非零元最多不超过13个；E为网格结点上待求的电场各分量；S为离散的边界条件。

应用Krylov子空间迭代方法（如QMR、BICG等）来求解方程式（19），再通过插值和辅助场计算得到测站处的电场和磁场值，从而可以计算阻抗、视电阻率和相位等参数。

3.2　三维大地电磁数据反演原理

定义目标函数为：

Φ（m）＝Φ_d＋λ^-1Φ_m

图12　Yee网格中的电磁场各分量采样方案

式中：d^obs为观测数据；d^pre为正演预测数据；m为待求的模型；m_ref为参考模型；β为正则化因子；C_d为数据协方差矩阵；C_m为模型协方差矩阵。将第k＋1次迭代的模型m_k＋1代入式（6）并将正演算子进行泰勒展开，可以得到模型域的Occam反演的迭代式为：

式中：J_k为第k次迭代所得模型m_k处的数据灵敏度矩阵。对式（21）进行一些数学运算可以得到其等价变换为：

在式（21）中，求逆部分的矩阵的阶数为模型参数的个数m，而在式（22）中，求逆部分的矩阵阶数等于数据个数n，故将式（22）给出的迭代式称为数据域Occam反演方法。在三维大地电磁数据的反演中，通常情况下都有n远小于m，因而使用数据域的Occam反演可以大大减小反演迭代过程中大型矩阵求逆的计算量，使得三维的大地电磁数据反演在PC上得以实现。

图13　低阻块体模型三维透视图

3.3　合成数据的反演

为检验反演算法的性能，设计一个电性模型为一个低阻块体赋存于50Ω·m的半空间中。低阻体的电阻率为1Ω·m，规模为16 000m×16 000m×5 000m（图13～图15）。对该低阻块体模型进行正演计算，在低阻体上方所在区域设置6×6的二维测网，每个测点观测5个频点，周期分别是0.1s、1s、10s、100s、1 000s。在正演数据中加入5%的噪声以模型实际数据，再将加入噪声的正演数据作为观测数据输入反演程序进行反演计算，反演结果如图16～图18所示。

图14　低阻块体模型东西方向切片图

图15　低阻块体模型顶面切片图

图16　低阻块体反演恢复模型三维透视图（电阻率小于10Ω·m的区域）

图17　低阻块体反演恢复模型东西方向切片图

图18　低阻块体反演恢复模型顶面切片图（深度100m处）

反演的初始模型为50Ω·m的半空间，经过5次迭代后，数据拟合差达到目标值，共耗时约5.5小时。

对比反演结果与真实模型可以看出：反演结果基本恢复了真实模型的主要构造信息，低阻区域的位置、形态和规模均与真实模型吻合良好。不过，虽然反演结果中的低阻区域比较光滑，但还是可以明显看出在低阻块体的4个角点附近存在着多余的构造，这与模型目标函数的选取相关。虽然本例所选取的数据点较少，但仍然较好地恢复了真实模型，这反映了Occam反演的确非常稳定可靠。同时，也因为问题的模型很小，所以，反演的计算耗时比较合理。但如果问题的规模稍微加大，则计算时间会显著增加，因而，对于计算量巨大的问题依然是三维大地电磁数据反演所面临的挑战。

4　展　望

从国内外检索的文献资料来看，MT法本身呈现如下发展趋势^［47］：①空间高密度采集^{［10，48］}；②全张量信息的深入分析和利用^［31］；③三维多参数成像和多数据源联合成像。从仪器设备和数据处理方面来看，该方法已发展到相当成熟阶段，但三维带地形的反演仍是一个难题。近年来三维反演方法越来越受到重视，这对减小包括近地表电性不均匀体引起的电磁场畸变在内的体积效应作用巨大。大地电磁数据三维反演中最本质的困难来自于计算资源的限制和反演结果的可靠性评价。在计算资源足够的情况下，反演方法本身比低维反演不存在更多的技术难点，但由于解空间的急剧增大，造成反演结果的非唯一性更为严重，对解的评价就显得非常重要。改善三维反演唯一性的途径只能是充分利用各种已知信息构建尽可能接近于真实模型的反演初始模型或者对解模型施加约束。三维反演解的评价问题迄今无人涉足，仍是一项难度很大的挑战。

对大地电磁成像结果的解释方面，令地球科学家们较为关注的是中下地壳高导层和上地幔高导层的解释。对于中下地壳存在区域性分布的高导层，其成因机制包括^［49］：①高传导性的矿物如石墨、金属硫化物等；②相互联通的空隙中充满含盐流体；③部分熔融；④岩石力学强度软弱带。一般认为上地幔高导层对应于地幔矿物（玄武岩、硫化物等）的部分熔融。利用地幔矿物的电阻率随温度升高而降低的特性，可以利用大地电磁成像结果粗略估计地幔的温度分布。另外，由于引起电各向异性和速度各向异性的机制不同，地震各向异性主要决定于具不同力学性质基质矿物的定向排列，而电各向异性主要决定于岩石空隙内传导物质分布的优势方位，因而综合利用地震成像和电磁成像的结果，可以约束地壳和上地幔的组构和热力学状态的解释。综合地球物理多分支学科的成果，重视高温高压岩石物理实验证据，从而对深部电性结构作出合理的解释是地球电磁学家们长期努力的方向。

地球物理观测、数据分析及成像方法技术的进步是其在地球科学研究和应用领域取得成功的3个关键要素。地球物理学的发展历史给出的基本逻辑是：基于地球物理场理论的探测方法和观测技术的提出—仪器装备的产生和更新—数据分析和成像方法推陈出新。通俗地说，新的观测技术的提出和推广使用，将催生新一代仪器装备的诞生，相应地将带动新的数据分析和成像方法的进步，MT法也不例外。因此，有理由相信，带更多通道的分布式可同步观测和数据传输的MT仪器的推出，与流动地震台站观测给地震学带来的影响一样，会对MT法产生巨大的推动作用，如本文提出的多站叠加技术可以演变为多道叠加技术或者聚束处理技术。

MT法自诞生至今已有60多年，与中国地质大学的建校历史相当。前辈们为MT法的发展做出了重要贡献，吾辈也需继续努力，继往开来。更期待更多的青年才俊投身这一激动人心的、“绿色”的（对人类和环境无影响）、应用日益广泛的地球物理方法研究中来。

参考文献

［1］王家映.大地电磁拟地震解释法［M］.北京：石油工业出版社，1995.

［2］王家映.地球物理反演理论（第二版）［M］.高等教育出版社，2002.

［3］Chen L S，Booker J R，Jones A G，et al.Electrically conductive crust in southern Tibet from INDEPTH magnetotelluric surveying［J］.Science，1996，274：1 694-1 696.

［4］Wei W B，Unsworth M，Jones A G，et al.Detection of widespread fluids in the Tibetan crust by magnetotelluric studies［J］.Science，2001，292：716-719.

［5］Nichols E A，Morrison H F，Clarke J.Signals and noise in measurements of low-frequency geomagnetic fields［J］.J.Geophys.Res.，1988，93：13 743-13 754.

［6］Pedersen L B.Some aspects of magnetotelluric field procedures［J］.Surv.Geophys.，1988，9：245-257.

［7］Brasse H，Junge A.The influence of geomagnetic-variations on pipelines and an application for large-scale magnetotelluric depth sounding［J］.Journal of Geophysics-ZEITSCHRIFI FUR GEOPHYSIK，1984，55（1）：31-36.

［8］Petiau G，Dupis A.Noise temperature coefficient，and long time stability of electrodes for telluric observations［J］.Geoph Prosp，28：792-804.

［9］王家映，徐义贤，张胜业，等.国外大地电磁响应函数估计方法［J］.地学前缘，1998，5（1-2）：217-222.

［10］Egbert G D.Processing and interpretation electromagnetic induction array data［J］.Surv.Geophys.，2002，23：207-249.

［11］Goubau W M，Gamble T D，Clarke J.Magnetotelluric data analysis：removal of bias［J］.Geophysics，1978，43：1 157-1 166.

［12］Gamble T C，Goubau W M，Clarke J.Magnetotellurics with a remote magnetic reference［J］.Geophysics，1979a，44：53-68.

［13］Gamble T C，Goubau W M，Clarke J.Error analysis for remote reference magnetotellurics［J］.Geophysics，1979b，44：959-968.

［14］Jones A G，Chave A D，Egbert G，et al.A comparison of techniques for magnetotelluric impedance estimation［J］.J.Geophys.Res.，1989，94：14 201-14 213.

［15］Tarakura S，Takeda M，Matsu K.Efferts of regional noise of magretotellurics and tleir removal by far remote reference method（in Japanese with English abstract）.［J］.Butsui-Tansa，1994，47：24-35.

［16］黄哲.大地电磁测深远参考处理原理及效果分析［J］.地震地质，2001，23（2）：212-216.

［17］陈清礼，胡文宝，苏朱刘，等.长距离远参考大地电磁测深试验研究［J］.石油地球物理勘探，2002，37（2）：145-148.

［18］杨生，鲍光淑，张全胜.远参考大地电磁测深法应用研究［J］.物探与化探，2002，26（1）：27-31.

［19］Shalivahan B，Bhattacharya B.How remote can the far remote reference site for magnetotelluric measurements be？［J］.J.Geophys.Res.，2002，107，doi：10.1029/2000JB000119.

［20］Jiang L，Xu Y X.Comparisons of the magnetotelluric experimental results by multisite remote reference and multisite superposition［M］.Beijing：International Workshop on Gravity，Electrical ＆Magnetic Methods and Their Applications，by China Geophysical Society and Society of Exploration Geophysicists，2011.

［21］Xu Y X，Jiang L.Theoretical and measurement investigations of magnetotelluric multistation super-position method［D］.The International Symposium on Deep Exploration into the Lithosphere，Scientific Proceedings，Chinese Academy of Geological Sciences，Beijing，2011，P30.

［22］Horn R A，Johnson C R.Matrix analysis［M］.Cambridge University Press，1985.

［23］张贤达.矩阵分析与应用［M］.北京：清华大学出版社，2004.

［24］陆基孟.地震勘探原理（上册）［M］.东营：石油大学出版社，1993.

［25］Groom R W，Bailey R C.Analytic investigations of the effects of near-surface three dimensional galvanic scatterers on MT tensor decompositions［J］.Geophysics，1991，56（4）：496-518.

［26］Schlumberger C，Schlumberger M，Leonardon E G.Some observations concerning electrical measure ments in anisotropic media and their interpretation［J］.Tronsactions of the American institute of Mining Engineers，1934，100：159-182.

［27］Hill D G.Alaboratory investigation of electrical anisotropy in precambrian rocks［J］.Geophysics，1972，37：1 022-1 038.

［28］Herwanger J V，Pain C C，Binley A.Amisotropic resitivity tomogrophy［J］.Geophys Jint，2004，158：409-425.

［29］Eisel M and Haak V.Macro anisotropy of the electrical conductivity of the crust：A magnetotelluric study of the German Continental Deep Drilling Site（KTB）［J］.Geophys.J.Int.，1999，136：109-122.

［30］Bahr K，Bantin M，Schneider E，Storz W.Electrical anisotropy from electromagnetic array data：Implications for the conduction mechanism and for distortion at long periods［J］.Phys.Earth Plane.Interi.，2000，119，237-257.

［31］Wannamaker P E.Anisotropy versus heterogeneity in continental solid earth electromagnetic studies：Fundamental response characteristics and implications for physicochemical state［J］.Surveys in Geophysics，2005，26：733-765.

［32］Darboux-Afouda R and Louis P.Contributino of electrical anisotrophy to research in froctwed aquifers in the crystalline region of Benin［J］.Geophys.Prosp.，1989，37：91-106.

［33］Ritzi R，Andolsek R.Relation between anisotropic transmissivity and azimuthal resistivity surreys in shallow，fractured，carbonate flow systems［J］.Ground Water，1992，30：774-780.

［34］Weidelt P.3-D Conductivity models：implications of electrical anisotropy：in Three-dimensional electro magnetics［J］.Oristaglio M and Spies B.SEG，1999：119-137.

［35］Pek J，Santos F A M.Magnetotelluric impedances and parametric sensitivities for 1-D anisotropic layered media［J］.Computers ＆Geosciences，2002，28：939-950.

［36］林长佑，杨长福，武玉霞，等.各向异性介质的大地电磁反演及其有关地震前兆场的讨论［J］.地震学报.1996，18（3）：326-332.

［37］Zonge K L and Hughes L J.Controlled source andio-frquency magretotelluric in Electromagnetic methods in applied geophysics［J］.Nabighian，M.N（ed.）SEG.Tulsa，1991，2（B）：713-809.

［38］Xu Y X，Yue R.Preliminary understanding of near-field effect of CSAMT in electric azimuthally anisotropic half-space［J］.Journal of Environmental and Engineering Geophysics，2006，11（2）：67-72.

［39］Li X and Pedersen L B.The electronagnetic response of an azimuthally anisotropic half-space［J］.Geo-physics 1991，56：1 462-1 473.

［40］Yin Changchun.Inherent nan-uniqueness in magnetetolluric inversion for 10anisotropic model［J］.Geophysics，2003，68：138-146.

［41］Mackie R，Madden T.Three-dimensional magnetotelluric inversion using conjugate gradients［J］.Geophys.J.Int.，1993，115：215-229.

［42］Mackie R，Madden T，Wanamaker P E.Three-dimensional magnetotelluric modeling using difference equations-theory and comparisons to integral equation solutions［J］.Geophysics，1993，58（2）：215-226.

［43］Mackie R，Smith T，Madden T.Three-dimensional electromagnetic modeling using finite difference equation：the magnetotelluric example［J］.Radio Science，1994，29（4）：923-935.

［44］Newman G，Alumbaugh D.Three-dimensional magnetotelluric inversion using nonlinear conjugate gradients［J］.Geophys.J.Int.，2000，140：410-424.

［45］Siripunvaraporn W，Egbert G.An efficient data-subspace inversion method for 2-D magnetotelluric data［J］.Geophysics，2000，65（3）：791-803.

［46］Siripunvaraporn W，Egbert G，Lenbury Y，et al.Three-dimensional magnetotelluric inversion：data-space method［J］.Physics of the Earth and Planetary Interiors，2005，150：3-14.

［47］徐义贤.西北与华南阵列式区域大地电磁场标准网控制格架示范性实验研究［D］.《深部探测技术与实验研究》国家专项课题设计书，2008，55.

［48］Patro P K，Egbert G D.Regional conductivity structure of Cascadia：preliminary results from 3D inversion of USArray transportable array magnetotelluric data［J］.Geophys.Res.Lett.，2008，35，doi：10.1029/2008GL035326.

［49］徐义贤.中下地壳高导层成因研究综述［J］.地质科技情报，1995，14（3）：15-22.