1 引言

大地电磁法是研究地壳及上地幔电性结构的主要手段，近年来，我国学者取得了显著的进展^[1-4]. Wanamaker^[5]在2005年详细阐述了电性各向异性在地电模型和地球动力学解释中的重要作用，同时说明了电性各向异性现象在地球内部是普遍存在的.二维大地电磁电性各向异性正演模拟是通过解偏微分方程或积分方程的数值计算方法实现的，被广泛使用的两种方法是：有限元法和有限差分法^[6].早在1975年，Reddy ^[7]等采用有限元法对偏微分方程做了数值计算，较早的开展了二维电性各向异性问题的理论研究工作.虽然在其后的二十多年间（直到2000年之后才有学者将其理论改进后用于大地电磁实测资料处理）未能广泛运用在大地电磁实测资料处理中，但其所做研究在理论意义上的影响却尤为深远.徐世浙等^[8]用有限元法计算了二维电性各向异性地电断面的大地电磁场，是我国较早见于公开发表的理论研究文章.李予国^[9]公布了用于计算二维电性各向异性结构感应电磁场的有限元公式.该公式采用有限元法求解表示平行于走向的场分量E_x和H_x的一组耦合方程组得到了有限元的线性方程，使用预条件共轭梯度法求得这个线性方程的结果；通过样条插值对E_x和H_x进行数值微分求得垂直于走向的场分量E_y和H_y.随后，与Pek等^[10]数值计算的结果做了对比，证明了理论的正确性.李予国等^[11]采用自适应网格剖分法对之前的有限元法正演计算的二维电性各向异性结构的大地电磁响应做了一次成功改进，并取得了很好的效果. Pek等^[10]在Reedy等^[7]研究基础上，将前人（Haak^[12]，Brewitt-Taylor & Weaver^[13]，Červ & Praus^[14]等）在二维电性各向同性介质中使用的有限差分法引入求解二维电性各向异性问题中，用有限差分法对偏微分方程做了数值模拟计算，是正演模拟计算二维电性各向异性介质中大地电磁响应问题影响较为深远的一种求解方法.

2 大地电磁电性各向异性二维正演理论

本文采用有限差分法对大地电磁二维电性各向异性正演问题进行研究，参照了Pek等^[10]的研究方法.

假设在笛卡尔坐标系中，有一个二维电性结构，它的构造走向平行于x轴，y轴垂直于x轴，z轴垂直于xy平面，且正向向下.模型由一个有限区域组成，在其中的电性各向异性区域是二维有限结构.假设地表为水平，对应z=0.在地表上空（z＜0）为绝缘空气层.来自z→-∞的平面波垂直向下传播，ω=2π/T，T为周期，单位为秒（s）.

谐变电磁场的麦克斯韦方程组为：

(1)

(2)

上式为描述电磁场分布的微分方程.e^-iωt为时谐因子，σ为张量电导率.x轴平行于构造走向，由2-D性质∂/∂x=0，公式（1）-（2）分解为：

(3)

(4)

(5)

(6)

(7)

(8)

在公式（3）-（8）中，如果E_x和H_x已知，那么E_y，E_z，H_y和H_z均可求得；将公式简化为只含E_x和H_x的一组二阶偏微分方程.

电场沿走向时（即TE模式）：

(9)

磁场沿走向时（即TM模式）：

(10)

公式（9）-（10）中：

A=（σ_xyσ_yz-σ_xzσ_yy）/D，B=（σ_xyσ_yz-σ_xzσ_zz）/D，D=σ_yyσ_zz-σ_yz².

用有限差分法求解公式（9）-（10）：首先将二维结构投影到一个数值网格（是一个有限网格区域），在电性各向异性区域，再分为一个矩形单元的电性均匀体系.网格通常是不规则的，首先应该符合模型的几何形式，同时满足数值运算的一般规则，即设计的数值网格应该满足数值模拟的需要.

使用Červ和Praus^[14]的方法，将积分插值法用于导出各个网格节点的有限差分方程.由图 1所示，黑色圆实点代表节点（j，k），阴影区域代表该节点四周的矩形积分单元G_jk，通过该积分单元对方程（9）和（10）求积分.矩形积分单元G_jk表示为：

求解这些积分，需用各自剖分网格四周的四个节点做近似计算.

公式（6）对应的积分形式为：

(11)

公式（11）中Γ（G_jk）是积分单元G_jk的定向边界，dg是沿着此边界的积分路径的元素，dG=dGe_x，e_x是在x方向的单位矢量，dG是G_jk的区域元素，表示为：

(12)

公式（12）中是指函数穿过积分单元边界KL的值，相应的，，意义可类推. 是积分单元点K、L之间的距离（从x正向看，沿逆时针方向为正），S（G_jk）是积分单元的面积，j_x是电流密度的x分量，j_x=（σE）_x=σ_xxE_x+σ_xyE_y+σ_xzE_z，是j_x穿过积分单元G_jk的值.

基于公式（12）的近似计算，得到了一个在网格节点（j，k）对于TE模式的近似差分方程，这个近似差分方程代表了公式（9）的有限差分形式：

(13)

公式（13）中p可以是j-1，j和j+1，q可以是k-1，k和k+1.C表示各个剖分节点上是否有E_x和H_x分量，其上标的第一个E表示为TE模式，第二个E和H分别为电场分量和磁场分量；且当（p，q）∈ {（j-1，k-1），（j-1，k+1），（j+1，k-1），（j+1，k+1）}时，有：

对应的TM模式积分形式，由公式（3）在积分单元G_jk上的积分表示为：

(14)

同样可以获得在网格节点（j，k）对于H极化下的线性差分方程：

(15)

与TE模式的线性差分公式（13）相比，在公式（15）中，当

（p，q）∈ {（j-1，k-1），（j-1，k+1），（j+1，k-1），（j+1，k+1）}时，只有：

在这些节点，相应的磁场系数C_jk^HH（p，q）一般不为零，它代表了电性各向异性的影响.系数矩阵C_jk^αβ（p，q）的元素并不是每一个都有数值，也就是说在任意情况下，都不用计算全部9个节点的系数元素（18个值）.对于TE模式，有5个节点（10个值）；对于TM模式，有全部9个节点，但只有14个非零系数元素.在全电性各向异性情况下，系数参数的个数会达到最大值.但随着网格单元包含的电性各向异性的减少，非零的系数参数也会逐渐减少.

用有限差分法对公式（9）和公式（10）做数值近似后，须将线性代数公式（13）和公式（15）安排在一个系统，以便做进一步的处理.因为包含了E_x和H_x两个变量，且这两个变量不在相同的网格节点，TE模式的方程需要在全部网格节点做近似计算，而TM模式的方程只需在导电底层内部做近似计算.

由图 2可知，系数矩阵是一个对称、带宽限定矩阵，可以使用高斯消元对其进行特定的修改，从而达到求解有限差分线性代数方程组的目的.为了存储矩阵的上半频带，需要对高斯消元公式做一次改进来实现，复数（N_STOR）需被放置在内存中，N_STOR为：

(16)

公式（16）中N是水平网格步长的数量，M_E、M_A分别是导电地层和空气中垂向网格步长的数量.

3 二维电性各向异性模型计算与分析

大地电磁法被普遍应用于地壳与上地幔大尺度构造研究中，在一般情况下，地形的起伏相对于测量剖面的长度较小，所以用水平层状二维电性各向异性结构的正演模拟，可满足多数地质条件下电磁传播特性的研究目的.

在两层电性各向同性介质中放置一个二维电性各向异性介质，对地电模型做正演计算，分析研究介质中大地电磁场的形态特征和分布情况.第一层称为介质1，为各向同性，在第二层各向同性介质（称为介质3）中放置各向异性介质2.模型参数为：介质1：ρ₁=300 Ωm；介质2：ρ₂¹=10 Ωm，ρ₂²=100 Ωm，ρ ₂³=100 Ωm，α=30°，其中α表示测量轴与电性主轴的夹角；介质3：ρ₃=1000 Ωm.介质1深度范围是0 km到-4 km；介质2的大小为20×4 km²，在水平方向上的位置为24 km至44 km之间、深度方向上的位置为-6 km至-10 km之间.如图 3所示：

正演计算时将模型剖分为45×35的剖分单元（垂向35个剖分单元中包含10个空气层的剖分单元），测点数为45个（即在地表的网格节点处），采用21个周期点（0.1，0.3，0.5，0.7，1，3，5，7，10，30，50，70，100，300，500，700，1000，1300，1500，1700，2000，单位为s）.计算结果如图 4所示（图中视电阻率拟断面图单位为Ωm，相位拟断面图单位为°）：

图 4展示了正演计算得到的大地电磁响应结果.正演计算的结果能够很好的反映出电性各向异性介质的存在范围和区域：TE模式的拟断面图清晰地反映了电性各向异性介质在纵向上的顶部埋深和延展范围；而TM模式的拟断面图反映电性各向异性介质在横向上的延续范围和长度比较准确.图 4也很好的反映了地电模型的层状性质，在频率log（f）=-1以上部分，视电阻率值表现出了覆盖层的相对低阻特性；在频率log（f）=-1以下部分，视电阻率值表现出了背景场值的高阻特性.

对图 3的地电模型做一个微小的调整，将介质2的参数改为：ρ₂¹=10 Ωm，ρ₂²=100 Ωm，ρ₂³=100 Ωm，α=0°.即：令测量轴和电性主轴相互重合，但介质电阻率值保持不变.通过正演计算，研究两种地电模型中大地电磁场分布的差异，计算结果如图 5所示（图中视电阻率拟断面图单位为Ωm，相位拟断面图单位为°）：

图 5为修改后模型正演计算得到的大地电磁响应结果.图 5与图 4在形态上大体相似，但细节处又有明显不同，主要表现在：TE模式下图 5（a）的低阻区域较图 4（a）有明显扩大，图 5（c）的高值区域有明显扩大，而低值区域的延展范围有了明显缩小；TM模式下图 5（b）的低阻区域较图 4（b）在水平方向上略有收窄，图 5（d）中高值区域的延展范围有了明显缩小；TE模式下图 5（a）的视电阻率最小值较图 4（a）有了明显降低，而TM模式下图 5（b）的视电阻率最小值较图 4（b）有了明显升高.这些细微的不同表明了观测角对于大地电磁响应结果的影响和作用：观测角的存在会令两个电性主轴的电阻率值有平均的效果，即具有低值电性主轴表现出来的视电阻率值不会太低，具有高值电性主轴表现出来的视电阻率值不会太高；但当电性主轴为0时，具有低值电性主轴表现出来的视电阻率值会降低，具有高值电性主轴表现出来的视电阻率值会升高.

对图 3的地电模型再做一个微小的调整，将介质2的参数改为：ρ₂¹=10 Ωm，ρ₂²=1000 Ωm，ρ₂³=1000 Ωm，α=0°.即：不仅令测量轴和电性主轴相互重合，而且使ρ₂²、ρ₂³的电阻率值等于背景场值.计算结果如图 6所示（图中视电阻率拟断面图单位为Ωm，相位拟断面图单位为°）.

由图 6的正演计算结果可知：在TM模式下，无法识别该模型的二维结构，只表现出了层状介质的电性各向同性性质.原因是当观测角为0时，测量轴与电性主轴方向是重合的，在电阻率值和背景值是相同的情况下，该方向具有层状同性电性特征.在TE模式下，与图 5一样，可以明显的看到一个低阻异常体的存在，并且，低阻异常的范围较初始模型计算得到的结果（图 4（a））有了非常明显的扩大，特别是在深度方向上异常的范围比水平方向上的异常范围要延展很多.

通过对二维电性各向异性介质的正演计算，并对正演计算结果进行研究分析，能够做出以下结论：在电性各向异性主轴同测量轴存在夹角的情况下，TE、TM模式都可以有效的分辨出异常区域，并对二维结构在顶部埋深、深度和宽度的延展上都有非常好的反映；观测角的存在会令两个电性主轴的电阻率值有平均的效果，在电性各向异性主轴同测量轴不存在夹角的情况下，若测量轴其中的一个方向上的电阻率值同背景场值相同，必然会造成对应的极化模式无法分辨出二维电性各向异性结构.

4 大地电磁实测资料解释及电性结构分析

选取的大地电磁测区位于新疆北部，该测线（图 7所示）位于准噶尔盆地的古尔班通古特沙漠边缘，地势平坦，整条测线高程起伏较小，海拔高度在300~400 m.测点均在沙漠边缘，干扰较小.中心测点（C）位置为N45°00′，E86°00′，西一测点（W1）距中心点约15 km，西二测点（W2）距中心点约40 km，东一测点（E1）距中心点约20 km，东二测点（E2）距中心点约40 km.

测线各个测点的大地电磁视电阻率和阻抗相位的曲线形态都非常平滑，仅在个别频点上存在飞点，只需做简单的圆滑处理.对此测线的大地电磁数据选取了频率为11.2~0.0005 Hz之间的数据进行分析研究.该测线的视电阻率拟断面图和阻抗相位拟断面图如图 8所示（图中视电阻率拟断面图单位为Ωm，相位拟断面图单位为°）.

因为测点较少，必须通过插值成图，造成了阻抗相位拟断面图的成图效果较差，但视电阻率拟断面图的成图效果还是能够满足研究要求的.在以往的研究中，基本上是通过对整条测线的大地电磁数据做TE、TM或TE+TM反演，从而研究测区的地下电性结构，反演方法也多种多样（Occam^[15-16]、RRI^[17]、NLCG^[18-19]等），这些研究的理论基础都是建立在电性各向同性理论上.本文尝试通过建立一个二维电性各向异性地电模型，基于前文所阐述的电性各向异性理论，对其做数值模拟计算，将计算结果同图 8做对比，从而分析该测区的电性结构.建立的二维电性各向异性地电模型如图 9所示.

由前文所述，整条测线的高程起伏较小，故图 9中的地电模型不包含地形起伏.图 9中该地电模型包含6个区域（分别由圆圈中的数字所指示），其中2、3、6、8区域为电性各向异性区域，4、5为电性各向同性区域.模型参数如下所示：介质2：ρ₂¹=36 Ωm，ρ₂²=30 Ωm，ρ₂³=36 Ωm，α=20°；介质3：ρ₃¹=30 Ωm，ρ₃²=20 Ωm，ρ₃³=30 Ωm，α=30°；介质6：ρ₆¹=2 Ωm，ρ₆²=100000 Ωm，ρ₆³=2 Ωm，α=0°；介质8：ρ₈¹=10 Ωm，ρ₈²=7 Ωm，ρ₈³=10 Ωm，α=0°；介质4：ρ₄=2 Ωm；介质5：ρ₅=15 Ωm.各个介质的区域范围如图 9所示，这里不再一一叙述.正演计算时将模型剖分为81×55的剖分单元（垂向55个剖分单元中包含10个空气层的剖分单元），测点数为81个（即在地表的网格节点处），采用28个周期点（0.1，0.3，0.5，0.7，1，3，5，7，10，20，30，40，50，60，70，80，90，100，200，400，500，600，700，900，1000，1300，1500，2000，单位为s）.计算结果如图 10所示（图中视电阻率拟断面图单位为Ωm，相位拟断面图单位为°）.

为方便对比，图 10使用了同图 8一致的色标，可以直观看到，不论是TE模式还是TM模式，视电阻率拟断面图都有很好的拟合；阻抗相位拟断面图在整体上差强人意，主要的差别出现在log（f）≥0和50~70 km之间的区域.对于整个低阻背景场（2 Ωm）来说，存在高阻异常区域的3个主要部分分别是：log（f）≥0和0~50 km之间的区域；log（f）≥0和75~80 km之间的区域；-1≤log（f）＜0和20~50 km之间的区域.无论是高阻异常区域的视电阻率值，还是异常区域的空间位置以及在水平和深度方向的延伸，图 10中的视电阻率拟断面图中都有很好的反映.由此可以判断，所建立的二维电性各向异性地电模型能够反映出测区的电性结构.由图 9可知，该区域的电性结构相当复杂，存在多个电性各向异性区域，各个电性区域间的相互影响也非常大，特别是在-1≤log（f）＜0和20~50 km之间出现的高阻异常区域，此区域在TE模式中没有任何异常出现，而在TM模式中变为一个非常明显的高阻异常区（相对于背景场值），这一现象本文在前文中有过详细阐述，属于电性各向异性结构的特例.在电磁法中，低阻背景场值会掩盖很多有用的电性信息，对高阻区域的影响尤为明显.在建立模型时，需要适当的扩大各个高阻区域的范围.在数值模拟计算时，对于此类背景场值和各个区域的电阻率值等级相差很大的情况，需要细分剖分网格，并加密频率，才能够将计算结果在成图后，很好的展现出各个不同电性区域的边界和范围.

对于此测区所做的二维电性各向异性结构的研究，不仅验证了二维电性各向异性基础理论及算法的正确性，也说明了在大地电磁资料中普遍存在电性各向异性现象，同时为今后的实际大地电磁资料处理提供了一种新的思路和方法.

5 结论

本文基于Pek等（1997）^[10]的研究成果，系统而详实的阐述了二维电性各向异性正演理论；并将电性各向异性理论成功的应用于大地电磁实测资料的处理解释中，为今后大地电磁资料的处理解释提供了新的思路和方法，主要取得了以下成果：

1）从麦克斯韦方程组出发，采用矩形剖分方法，使用有限差分法为正演研究工具，对二维大地电磁电性各向异性正演问题进行了系统研究.通过对二维电性各向异性地电模型的正演计算，研究了二维电性各向异性结构对观测大地电磁场的影响，为今后大地电磁实测资料的处理解释提供理论依据.

2）以本文的研究成果为基础，将电性各向异性理论引入对实测大地电磁资料的处理解释中，通过对新疆古尔班通古特沙漠边缘大地电磁资料做二维正演拟合解释，说明了电性各向异性现象是普遍存在的，也验证了理论的正确性及算法的实用性，能够为今后分析和解释大地电磁资料中的电性各向异性现象提供理论依据和技术指导，开拓了对大地电磁实测资料处理解释的思路和方法.

致谢

感谢捷克科学院地球物理研究所Josef Pek教授对本文的帮助.