1 引言

大地电磁法(MT)已广泛应用于地球物理探测各个领域.我国的山区面积比较大，在起伏地形条件下开展MT 工作时，地形的影响产生电磁场畸变，这给资料处理和解释工作带来很大困难，因此有必要研究带地形的MT 数值模拟算法.

目前，在二维起伏地形和地电模型条件下，国内外学者采用过多种有限元剖分方式的二维MT 有限元数值模拟.在有限元电性参数(如电导率、传播系数等)分块均匀的前提下，Wannamaker^[1]、Mauriello等^[2]采用矩形三角网格剖分、双线性插值有限元进行数值模拟，陈小斌^[3-4]、谢飞^[5]用自适应地形三角网格剖分、双线性插值有限元数值模拟，周熙襄^[6]、王绪本^[7]、刘云等^[8]采用自适应地形四边形网格剖分、双二次插值有限元数值模拟，对斜坡型地形产生的地形影响做模拟和分析;认为电性参数分块均匀的模拟方式，在单元块内的电性参数分布是均匀的，单元块与单元块之间的电性参数分布则不一定是连续的.但是在野外实际中，岩石、矿物的电性参数往往是变化的，因此采用电性参数分块连续变化的数值模拟方法，更符合野外实际.徐世浙、李予国等^[9-10]采用矩形网格剖分，进行电性参数分块连续变化、双线性插值、双二次插值等有限元数值模拟;对于矩形网格剖分，如果在有限元矩形单元块内或单元节点处填入空气介质，从而形成台阶型地形模拟地形影响^[11-12]，是有其局限性的^{[4， 13]}.

阮百尧^[14]在起伏地形二维直流电阻率法正演数值模拟中，采用矩形内三角网格剖分、电导率分块连续变化、双线性插值有限元数值模拟.本文在此基础上，导出了起伏地形三角单元网格剖分电性参数分块连续变化的二维MT 有限元数值模拟方法.

2 二维MT 变分问题

如图 1所示，假设三维地质模型的电、磁等参数沿构造方向(x方向)上无变化，这时可将三维模型问题转化为二维模型问题进行处理.此时，与二维MT 边值问题相应的变分问题^[13]是

(1)

式(1)中，u为E_x或H_x，Ω 为二维研究区域.

对于TE 模式(E_x型):上边界AB处离地面足够远.则

对于TM 模式(H_x型):上边界AB在地面上.则

其中，σ 和μ 分别为介质电导率和磁导率，ω 为圆频率，ε 为介电常数.下边界CD处离目标区域足够远，这时σ 为CD处均匀介质电导率;左、右边界AC、BD离目标区域足够远.

3 有限单元法

用有限单元法求解式(1)的二维MT 变分问题，具体步骤如下:

3.1 网格剖分

对于TE 模式，在高频时，空气介质中位移电流的影响不可忽略^[3]，因此研究区域为空气和地下;对于TM 模式，场值H_x在空气中近似为常数(在y和z方向的偏导数近似为零)，即不考虑空气介质中位移电流的影响，此时AB边界设在地面上，研究区域为地下.

如图 1所示，研究区域Ω用三角单元进行网格剖分.如图 1a所示，当模拟水平地形时，考虑到地下介质体形状的任意性，首先对研究区域进行矩形网格剖分，再在每一个矩形网格内剖分出4 个三角单元网格.这样，一方面，避免了三角网格过于尖锐的情况;另一方面，可利用三角形的斜边模拟地形线和任意介质体倾斜的界面;如图 1b所示，当模拟起伏地形时，在图 1a网格剖分的基础上，地空边界采用沿实际地形线进行网格剖分.

3.2 线性插值

如图 2所示，假设三角单元e内电磁场u和电性参数τ，λ，k呈线性变化，则在每一个单元中:

如果三角单元e的23边位于CD边界上时，

这里，u_i为三角单元节点处的电磁场，τi、λ_i、k_i为三角单元节点处的电性参数，(i)、1、2和3为三角单元的三角节点号.N_i=1/2Δ(ai_y+b_iz+c_i)是关于y和z的线性插值形函数(i=1，2，3)，其中Δ =1/2(a₁b₂ -a₂b₁)是三角单元的面积，其中

其中，(y₁，z₁)、(y₂，z₂)和(y₃，z₃)是三角单元节点号的坐标.

3.3 单元分析

在整个研究区域Ω 中，将泛函(1)式离散化，表示为所有三角单元e的线性组合，即

(2)

3.3.1 单元分析1

u对y求偏导数，有

其中，，.

显然，.同理可得，.

面积分，Δ 为三角单元的面积^[13].因此，对于(2)式中第一个积分项，有

式中，.

3.3.2 单元分析2

面积分，其中a，b，c，为非负整数^[13].

因此，对于(2)式中第二个积分项，有

式中，.

其中，c₁₁ = (6 2 2)λ，c₁₂ = (2 2 1)λ，c₁₃ =(2 1 2)λ，c₂₁= (2 2 1)λ，c₂₂= (2 6 2)λ，c₂₃= (1 2 2)λ，c₃₁= (2 1 2)λ，c₃₂= (1 2 2)λ，c₃₃= (2 2 6)λ，λ= (λ₁ λ₂ λ₃)^T.

3.3.3 单元分析3

线积分，其中a，b为非负整数，l为单元底边的长度^[13].

对于(2)式中第三个积分项只对CD边界进行.当单元的23边落在CD边界上时，线积分

式中，，l为23边长度.其中，d₂₂ = (12 3 3 2)τk，d₃₂ = (3 2 2 3)τk，d₃₃ = (2 3 3 12)τk，d₂₃ =d₃₂，τk= (τ₂k₂τ₂k₃ τ₃k₂ τ₃k₃)^T.

3.4 系数矩阵总体合成及求变分

将研究区域Ω 内每个单元系数矩阵ke按照总体节点号进行扩展，得到k_e，相加得到总体系数矩阵K，即

其中，K矩阵的数学表现形式为(i，j=1，…，NOD)，NE为总体单元数，NOD为总体节点号，e为单元号，i，j为总体节点号变量.显然，K是一个对称、正定的大型稀疏矩阵，可按照变带宽方式存储^{[13， 15]}.令F(u)的变分为零，得线性方程组Ku=0，代入第一类边界条件(u|_AB=1)，得Ku=P，解此线性方程组^{[13， 16]}，便可得各节点u值，即E_x或H_x值.

如图 3所示，采用矩形内三角网格剖分方式，比较矩形网格剖分来说，增加了一个中间节点5，大大增加了计算量.但是节点5 只与该矩形内4 个角点有关系，而与其它不相邻的矩形网格没有任何直接联系，所以在单元分析中，运用高斯消元法事先可消去这个节点.节点5只是一个虚设的节点，并不包含在总体节点数中.到解方程结束后，亦可根据矩形网格内4个角点的已知场值和节点系数之间的关系，直接计算出节点5的场值.这样，一方面不增加节点总数，节省了计算量，另一方面每个三角单元仍具有各自的不同物性.具体做法为:节点5 的消除，将k_i，j(i，j= 1，2，3，4，5)替换为k′_i，j，k′_i，j= k_i，j-;节点5 的恢复，u₅ =.

4 辅助场、视电阻率和阻抗相位的计算 4.1 辅助场的计算

辅助场是通过计算主场值沿某方向的方向导数得到的^[17-19].当各节点主场值求出后，由3.2 节可知，单元内场值u可表示为，则u对y求偏导数，得

其中a₁、a₂、a₃，是与三角单元的三顶点坐标有关的常数，所以在单元内，u的偏导数是常数.

同理可得

在大地电磁法的数值模拟中，常以地面三角单元边上任意点或三角单元的节点为实际测量点，此时要考虑到相邻单元间场值偏导数的不连续性问题^[13].对于彼此相邻的单元、相同单元边或节点上场值的偏导数，可以在所属不同单元中予以分别计算，再取平均值.本文以地面三角单元的节点为实际测量点.不难发现，对于TE 模式，需要计算同一测点的地表和地表上一空气层所有单元的场值偏导数，再取平均值;而对于TM 模式，则只需求出同一测点的地表所有单元的场值偏导数，再取平均值即可.

4.2 视电阻率定义和阻抗相位的计算

对于TE 模式:u=E_x，辅助场.在二维起伏地形条件下，E_x、H_y在野外实际中均可测得，因此，视电阻率的定义仍采用普遍形式的阻抗定义方式，即

其中，偏导数按4.1 节方式求出.所以，ρ_E= ρ ，;对于TM 模式:u=H_x，在二维起伏地形条件下，H_x可测得，E_y不可测得，只能测得沿地形线切线方向上的E_l.这时，辅助场E_l可以分解成E_y和E_z两个分量(x方向为二维构造方向)^[4]，即.当为水平地形情况时，E_z=0.因此，视电阻率形式为

其中，偏导数按4.1 节方式求出.所以，ρ_H = ρ ，

5 模型计算

为验证本文正演算法的有效性，首先与水平层状一维连续介质理论模型的解析法、以及前人已有地形模型做计算分析，之后用该算法对倾斜界面异常体模型做计算分析.

在以下所有模型的计算中，有限元网格的剖分方法如下:如图 1 所示，在横向网格中，目标区域采用等间隔网格，间隔距等于点距.左、右稀疏网格数各为18个，网格大小等于点距乘以系数，其中18个系数依次分别为1，1，1，1，1，2，2，2，2，4，4，4，8，8，16，32，64，128;在纵向网格中，空气稀疏网格坐标(TM 模式无空气网格，单位为m)14 个，依次分别为100000，50000，10000，5000，1000，500，200，100，55，35，25，15，10，5.地下网格的剖分，从1~10层的网格大小(单位为m)依次分别为5，5，5，5，10，10，10，15，15，20，第1层网格大小的默认值为5 m.当然，也可以任意设置第1 层网格大小，这样，从1~10层的网格大小则相应地等比例放大或缩小.从11层及以下的网格剖分为等比网格，该层的网格大小等于上一层网格大小乘以等比系数，等比系数一般为1.1或1.2.这样，横向网格数则取决于测点数的多少，纵向网格数则取决于探测深度的选择;当模拟起伏地形时，上、下地表附近还要加上地形网格，地形网格数等于所有测点的高程数(对于重复的高程，只计1个)，地形网格大小等于相邻高程之差.

5.1 一维连续介质模型

如图 4所示，三层一维连续介质地电断面模型，第一层电阻率为100Ωm、层厚度为1000m;第二层为连续介质层，介质电阻率随深度线性减小，其变化关系如图中所示，层厚度为1000 m;第三层电阻率为1Ωm;测点区域为0~4km, 点距为100m, 共41个测点，频率范围为10³~10^-3 Hz, 以10为底的对数间隔采样，共51个频点.

在一维层状介质的解析法中，将中间的连续介质层剖分为10层进行模拟.在二维有限元数值模拟中，取测线中心处第21 号点测点为研究对象.则解析法与二维数值解法的结果对比如图 5所示.

在TE、TM 视电阻率曲线图和阻抗相位曲线图上，解析法和数值解曲线形态、值的大小基本一致，这表明本文方法能够有效地模拟连续介质模型.正演结果的数据统计表明，视电阻率、阻抗相位的解析法结果和数值模拟结果均方根误差均小于1%.本文的数值模拟结果与解析法结果之间的误差分析主要为:一方面，与有限元数值模拟的网格剖分有关系，如高频网格和低频网格的剖分方式是不一样的;另一方面，解析法是基于电性参数分层均匀的计算方式，而本文的数值模拟方法是基于电性参数分块连续变化的模拟.

5.2 地形模型

如图 6 下部所示为圆柱凹陷地形模型，Wannamaker等^[1]和徐世浙等^[20]在文献中对这个模型做过模拟.测点区域为0~0.4km, 点距为10m, 共41个测点，圆半径为50 m, 频率为0.01 Hz, 计算TM 模式下视电阻率.

本文数值模拟结果如图 6上部曲线图所示，对照文献^{[1， 20]}中模拟结果，在0.01 Hz频点处，TM模式下，三种模拟方法与地形测点对应的视电阻率起伏形态和值的大小基本一致.

图 7所示为山脊地形模型，在文献^{[1， 20]}中对这个模型做过模拟.测点区域为0~4km, 点距为50m, 共81个测点，起伏落差为100m, 山底宽度为2.4km, 频率为10Hz.

本文数值模拟结果如图 8 所示.对照文献^{[1， 20]}中的模拟结果，在10 Hz频点处，TE、TM 模式下，三种模拟方法与地形测点对应的视电阻率、阻抗相位起伏形态和值的大小基本一致.

通过对以上两例均匀大地起伏地形模型的模拟，得到了与前人研究相符的结论.分块均匀法认为地面或地下单元内节点之间的电性参数是均匀分布的，而本文的分块连续法则认为地面或地下单元内节点之间的电性参数是连续变化的(空气介质除外).显然，分块均匀法是分块连续法的特殊形表现形式，分块连续法更具有普遍性.用三角网格模拟起伏地形、分块连续变化的剖分方式且更加符合野外实际地质情况.

5.3 倾斜界面异常体模型

如图 9所示为“W"字形倾斜界面异常体模型，背景电阻率为100Ωm, 异常体电阻率为20Ωm, 其上顶面距地面100m, 下底面距地面500m, 异常体上顶面外宽4800m, 下底面外宽2400m.测点区域为0~8km, 点距为100m, 共81个测点，频率为103~10-1Hz, 以10为底的对数间隔采样，共41个频点.

有限元数值模拟网格剖分为:最大探测深度选择为10km, 等比网格系数选择为1.1;TE 模式，网格数为116×66(横向网格数为116，左、右稀疏网格数各为18;纵向网格数为66，其中包括空气稀疏网格数为14);TM 模式，网格数为116×52(横向网格数同TE 模式;纵向网格数为52，没有空气网格).大型稀疏矩阵的存储方式均按文献^{[13， 15]}中变带宽方式存储，用不带平方根的Cholesky分解法求解线性方程组^{[13， 16]}.

本文数值模拟方法的结果剖面图如图 10所示，计算每一个频点的时间为0.78s(TE 和TM 模式同时计算).在TE、TM 两种极化模式下，视电阻率、阻抗相位剖面上都出现明显低阻异常，“W"形异常形态清晰可见，表明采用三角网格模拟任意倾斜界面异常体的方法是有效的.

6 结论

本文提出三角单元网格剖分、电性参数分块连续变化二维MT 有限元数值模拟方法，比较矩形单元剖分、电性参数分块均匀的模拟方法，这样既不增加节点，节省了计算时间，又能很好模拟野外实际地形和任意地电体模型.通过对一维连续介质模型模拟计算，与解析法计算结果对比，验证了本文算法的可靠性;对前人地形模型模拟计算，结果与前人结果相符;对任意倾斜界面异常体模型进行模拟，结果表明本方法是有效的.同时指出，本方法仍然是考虑各向同性介质情况，对各向异性介质分块连续变化的数值模拟方法将是本文的后续研究工作.

在本文程序的编制、有限元网格剖分方法以及辅助场的计算中，陈小斌研究员和作者进行多次坦诚交流和讨论，并提出许多宝贵意见，谨在此表示衷心感谢.

致谢

在本文程序的编制、有限元网格剖分方法以及辅助场的计算中，陈小斌研究员和作者进行多次坦诚交流和讨论，并提出许多宝贵意见，谨在此表示衷心感谢.