大地电磁测深法(Magnetotelluric，简称MT)是一种重要的地球物理勘探方法，具有不受高阻屏蔽和对低阻层反映灵敏的优点，探测深度可以到达下地壳和上地幔.自从20世纪50年代初由A.N.Tikhonov和L.Cagnird分别提出来后，其在国内外发展迅速.

有限单元法是20世纪50年代发展起来的一种求解偏微分方程有力的数值方法.其最早可追溯到Courant(1943)在一个子域上采用逐段连续函数来确定接近未知函数，使用了一组三角形单元和最小势能原理去分析扭转问题.将研究区域剖分成许多子单元，用古典变分原理和多项式插值的方式进行数值分析.有限元法的整体刚度矩阵是一个维度很大的矩阵，使得计算量很大，所以开始并没有得到广泛的应用.直到1960年以后，随着电子数值计算机的广泛应用和发展，有限元法的发展速度才显著加快(王勖成和邵敏，1997).

1971年Coggon提出用有限元法进行大地电磁场的数值模拟(Coggon，1971).他从电磁场最小能量原理出发推导出了赫姆霍兹方程的有限元格式，实现了二维大地电磁有限元正演计算.此后，有许多的学者和专家在大地电磁有限元正演方面进行了研究.Rodi将有限元法引入到MT正演计算中，采用矩形单元进行二维大地电测正演问题(Rodi，1976).Rijo用通用性网格剖分的有限元法，进行大地电磁正演模拟，给出了相应的程序(Rijo，1977).1981年，Pridmore发表了三维大地电磁有限元正演的研究成果(Pridmore et al.，1981).Wanamaker采用三角网格剖分进行了带地形的有限元正演(Wannamaker et al.，1986)，另外还给出了能直接解出平行走向的二次场有限元算法(Wannamaker et al.，1987).Franke等用自适应网格进行了各向异性介质大地电磁二维有限元正演(Franke et al.，2004).自从有限元法在大地电磁中的应用引入到国内后(徐世浙，1982)，越来越多的学者也在这方面进行了研究(徐世浙和刘斌，1995；史明娟等，1997；陈小斌等，2000；冯德山和王珣，2013；毕武和刘云，2014)；带地形的二维大地电磁正演问题的研究(王绪本等，1999；刘树才等，2005；刘云和王绪本，2010)；基于计算二次场响应的方法，进行二维大地电磁数值模拟(刘小军等，2007).随后许多学者开始逐渐研究如何提高有限元正演过程中大型线性方程组的计算精度和效率.如陈小斌研究出的有限元直接迭代算法，避开了传统有限元算法中的刚度矩阵的形成(陈小斌，1999；陈小斌和胡文宝，2002；陈小斌和赵国泽，2004)；使得正演过程速度快、稳定性好、正演结果精度更高的不完全LU分解预处理的BICGSTAB算法(柳建新等，2009；徐凌华等，2009).

这些有关文献中，很少见到二维大地电磁有限元正演的推导过程，以及节点编码的详细描述.本文基于前人的研究，介绍了二大地电磁有限元正演的原理、方法过程，详细反复推导，特别是形象描述了单元刚度矩阵合成整体刚度矩阵的过程，使得方法原理和过程实现变得更加清晰.

1 二维大地电磁的基本原理

二维模型是指介质的电学性质除深度外还沿一个水平坐标轴方向有变化，而沿另一个水平坐标方向保持不变(走向).本文中取x轴为二维介质的走向方向，z轴为深度方向，y轴为电学性质变化的方向(倾向).

在大地电磁测深方法中，考虑到应用的观测频率范围一般为10^-4～10³ Hz，构成地壳浅部的介质的电阻率一般取1～1000 Ω·m，所以可以忽略位移电流的影响(陈乐寿和王光锷，1990).其麦克斯韦方程组(取时域中的谐变因子为e^-iωt)为

将上式展开得到二维地电模型的偏微分方程

可知在忽略位移电流的情况下，空气中的H_x为常量.即使考虑位移电流，由于对大地电磁来说，有$\frac{\sigma }{\omega \varepsilon }\gg 1$，H_x也近似为常量.所以TM模式的研究区域为地下.对于TE模式，在高频时，不能忽略空气中的位移电流，因此研究区域为空气和地下.为方便起见本文取TM模式进行有限元法的正演数值模拟，令

TM极化模式外边界条件(区域示意见图 1)：

图 1

Fig. 1

图 1 节点编号和单元编号示意图 1 ，2，…代表节点号 ①，②，…代表单元号 Fig. 1 Node number and unit number

(1) 上边界AB直接取在地面上，并以该处的u为1单位，则有

(2) 下边界CD以下为均质的岩石，局部不均匀体所产生的异常场在CD上为零，电磁波在CD以下的传播方程为

(3) 取左右边界AD、BC离局部不均匀体足够远，电磁场在AD、BC上左右对称，其上的边界条件是

上述大地电磁的边界问题归纳为

根据边值问题与变分问题的关系，(12)式与下列问题等价：

二维有限元的网格剖分常用三角形单元或者矩形单元，本文采用矩形单元对整个区域进行剖分，整体节点编号和单元编号的顺序示意见图 1.

对于每个矩形单元内部，取8个点(4个顶点和4条边的中点)，进行双二次插值.单元内局部节点编号和坐标如图 2所示.

图 2

Fig. 2

图 2 母单元与子单元(a)母单元；(b)子单元. Fig. 2 Mother unit and sub unit

在图 2a所示的正方形单元(母单元)上构造如下的形函数：

形函数(14)式满足：

子单元上的u可表示为

由(15)可知

母单元和子单元间的坐标变换关系为

有限元法研究区域的场值问题的重要一步就将整个区域的积分分解个各单元的积分之和，因此在对区域进行单元剖分后，将(13)的F(u)的积分进行分解，公式为

对(22)式分别进行单元积分，公式为

同理可得：

单元积分为

其中

式(22)最后一项线积分只对边界单元进行，当单元的 $\overline{263}$ 边落在无穷远边界CD上时，线积分为

以上过程得到的K_1e、K_2e和K_3e为各单元的刚度矩阵，之后需要将各个单元的刚度矩阵进行合成，得到整体刚度矩阵.各个单元中的元素在整体刚度矩阵中的位置，是由网格中的整体编码和单元中局部节点编码的关系决定的.

以图 1中的①单元为例，其局部编码如图 2a所示.即单元①局部节点编码1、2、3、4、5、6、7和8在整体结构中的编码分别为1、3、14、12、2、9、13和8.①单元的单元刚度矩阵中的元素在整体刚度矩阵中的位置如下图所示.

图 3

Fig. 3

图 3 局部编码和总体编码的关系 Fig. 3 The relationship between local coding and global coding

图中k₁₁^①表示①单元的单元刚度矩阵中的元素.将K_1e、K_2e和K_3e扩展成全部节点组成的矩阵，再将3个扩展矩阵相加得

将上边界值代入，接线性方程组即得节点的u，它代表个节点的H_x.

MATLAB软件是一种功能强大，运算效率高的数字工具软件.它允许用户以数学形式的语言编写程序，使用上比C、Fortran和Basic等高级语言都要方便得多.

本文用MATLAB软件，编程实现上述二维有限元正演过程，给出了两组模型的正演响应.模型一和其响应如图 4所示，求解区域为4 km，异常体和围岩的电阻率分别是 10 Ω·m和100 Ω·m.模型二和其响应如图 5所示.

图 4

Fig. 4

图 4 模型一有限元正演响应 (a)试验模型一;(b)TM模式视电阻率响应;(c)TM模式相位响应. Fig. 4 The finite elementforward response of the model1

图 5

Fig. 5

图 5 模型二有限元正演响应 (a)试验模型二;(b)TM模式视电阻率响应;(c)TM模式相位响应. Fig. 5 The finite elementforward response of the model2

有限元法是一种强有力的数值模拟方法，用有限元法对大地电磁二维问题进行正演，其原理清晰易懂，运算结果稳定可靠.本文在详细推导了大地电磁二维问题的正演过程后，给出了正确的形函数、单元刚度矩阵并将局部编码和整体编码的关系用图表的方式展现出来，完整地表现了整个正演的原理、方法和步骤.

致 谢 感谢审稿专家提出的修改意见和编辑部的大力支持！