有限单元法与有限差分法在MT一维正演模拟中的对比分析

王 涛，柳建新，童孝忠＊，曹创华，谭神湘

（1.中国海洋大学 海洋地球科学学院，青岛 266100；2.中南大学 地球科学与信息物理学院，湖南 长沙 410083）

0 前言

大地电磁测深法（Magnetotelluric Sounding，简称MT）是一种重要的地球物理勘探方法，由苏联学者 Tikhonov［1］和Cagniard［2］于上世纪五十年代初期提出。与其它地球物理勘探方法一样，正演问题是大地电磁测深法的理论基础，同时也是我们认识各种地电条件下大地电磁场响应特征的良好途径，其研究一直受到广泛关注。通过对不同地质模型的正演研究，我们总结出在不同地质条件下大地电磁场的分布规律；同时通过正演，也能准确了解在不同的地形起伏情况下，大地电磁场的分布特点。

大地电磁测深正演模拟的数值方法主要分为三种：①有限差分法；②积分方程法；③有限单元法。作者主要针对有限单元法和有限差分法进行讨论。有限单元法（Finite Element Method或简称FEM）是将要分析的连续场分割为很多较小的区域。Coggon［3］首先将有限单元法应用在电磁法正演模拟中；七十年代末，朱伯芳［4］将有限单元法引入国内。有限差分法（Finite Difference Method或简称FDM）是一种经典的数值模拟方法。地球物理工作者从八十年代开始研究有限差分法正演计算问题，为有限差分法的发展做了相当大的贡献（周熙襄等［5］；罗延钟等［6－7］。根据长期的研究，有限单元法和有限差分法在实际应用中均有各自的优势和不足，作者在本文中主要对有限单元法和有限差分法进行对比分析。

1 数值模拟

1.1 大地电磁测深正演的基本理论

麦克斯韦方程组是电磁场必须遵从的微分方程组，利用傅里叶变换可将任意随时间变化的电磁场分解为一系列谐变场的组合，通常我们以e－iωt表示谐变场的时间因子（即以负谐时表示），电场强度和磁场强度可表示为式（1）与式（2）。

大地电磁测深所讨论的电磁场频率是极低的，故在大地介质中可忽略位移电流对场分布的影响，即在大地电磁测深正演中研究的是似稳电磁场问题。于是，导电介质低频谐变场的麦克斯韦方程组为式（3）～式（6）。

式中 ▽·E＝0是因为导电介质内部体电荷密度实际上为零，公式中时间因子都隐含在电场E和磁场H 中，方程组（3）～方程组（6）是大地电磁测深理论研究的出发点。

1.2 一维有限单元正演算法

电场满足的微分方程为公式（7）。

根据边值问题所对应的变分问题为

采用有限单元法进行计算，首先应将区域剖分为若干单元。在单元内电导率（或电阻率）必须是连续的，也就是说，电导率的间断点不能在单元内。对于高频，电磁波衰减非常迅速，如果采用线性插值，必须取很小的单元，从而增加了计算工作量。

（1）单元积分为式（9）。

（2）单元积分为式（10）。

然后将各单元的扩展矩阵相加，则得K＝∑K1e－∑K2e＋K3。再根据δF（Ex）＝0，可得

将电场所满足的上边界条件和下边界条件代入上面的方程组，则有

求解上述线性方程组式（12）即可得到节点处的电场值，从而也可以进一步计算模型响应的视电阻率和相位。

1.3 一维有限差分正演算法

根据电场所满足的微分方程，将一维地电模型离散化，如图1所示。

图1 地电模型离散化Fig.1 Geoelectric model discretization

然后对电场进行离散处理，采用有限差分近似计算一阶、二阶偏导数。写成矩阵形式如式（13）。

根据下边界条件，电场衰减为“0”，当上边界条件取Hy＝1，则（E1＋1／2－E1－1／2）＝iωμ，得到公式（14）。

求解方程组式（14）可得节点处电场值，从而可以进一步计算模型响应的视电阻率和相位。

2 模型算例

2.1 均匀半空间模型

设计一个均匀半空间模型，ρ＝10Ω·m，分别用有限单元法和有限差分法进行正演模拟。图2和图3分别显示有限单元法和有限差分法模拟的视电阻率的相对误差。

由图2和图3可知，两幅曲线图的变化趋势基本一致，在低频段，这两种方法的计算结果与实际值基本一致；而在高频段，则出现随着频率升高视电阻率的相对误差增大的现象。从整体上看，这两种方法的计算结果都基本符合正演的要求，但是有限差分法正演结果的相对误差比有限单元法正演结果的相对误差要小，说明有限差分法在均匀半空间模型的正演效果比有限单元法略好。通过以上模拟的计算结果对比，证明了两种方法正演算法的正确性，以及所编写程序的正确性。

2.2 层状介质模型

2.2.1 二层层状介质

设计一个二层层状介质模型，ρ1＝10Ω·m，ρ2＝100Ω·m，h1＝1 000m。

（1）解析法正演结果。正演计算的曲线如图4所示。

（2）有限单元法正演结果。有限单元法的网格剖分如表1所示，有限单元法计算结果与解析法计算结果比较如图5所示。

公式（13）如下：

公式（14）如下：

（3）有限差分法正演结果。有限差分法的网格剖分如表2所示，有限差分法计算结果与解析法计算结果比较如图6所示。

2.2.2 三层层状介质

（1）解析法正演结果。正演计算的曲线如图7所示。

（2）有限单元法正演结果。有限单元法的网格剖分如表3所示，有限单元法计算结果与解析法计算结果比较如图8所示。

（3）有限差分法正演结果。有限差分法的网格剖分如表4所示，有限差分法计算结果与解析法计算结果比较如图9所示。

图4 二层层状介质大地电磁响应曲线Fig.4 Two－tier media magnetotelluric response curve

表1 有限单元法网格剖分参数Tab.1 Finite element method meshing parameters

表2 有限差分法网格剖分参数Tab.2 Finite difference method meshing parameters

图7 三层层状介质大地电磁响应曲线Fig.7 Three－tier media magnetotelluric response curve

表3 有限单元法网格剖分参数Tab.3 Finite element method meshing parameters

表4 有限差分法网格剖分参数Tab.4 Finite difference method meshing parameters

（1）由图5和图6可得，在二层层状介质模型中，当频率较高时，有限单元法与有限差分法计算的视电阻率曲线与理论值曲线吻合非常好；而当频率较低时，两种方法计算的视电阻率曲线逐渐地脱离了理论的视电阻率曲线，在尾端出现了上升的趋势。

（2）由图8和图9可得，在三层层状介质模型中，有限单元法与有限差分法计算的视电阻率曲线与理论值曲线基本上一致，但二者在曲线的极大值处附近与理论值有所偏离。①有限单元法计算的视电阻率曲线与理论值曲线对比的趋势为：开始偏离后渐渐吻合，再过极大值点偏离，后渐渐吻合；②有限差分法计算的视电阻率曲线与理论值曲线对比的趋势为：开始吻合，过极小值点后偏离，再过极大值后吻合。

3 结论

作者在本文以大地电磁测深的基本理论为基础，从麦克斯韦方程组出发，以有限单元法和有限差分法为计算工具，研究了大地电磁测深的一维正演问题。通过对一些特定的模型进行正演模拟，验证了有限单元法和有限差分法在大地电磁测深正演计算中的有效性。从以上的二层层状介质模型和三层层状模型对比分析中，我们可以看出，有限单元法和有限差分法的计算结果与解析法的计算结果（即理论值）基本一致，两种方法的正演结果均真实地反映了模型的地电参数，说明了有限单元正演算法和有限差分正演算法的正确性，同时也说明了编写程序的正确性。

［1］TIKNONOV A N.On determining electrical characteristics of the deep layers of the earth’s crust［J］.Deki Akud Nuck，1950（73）：295－297.

［2］CAGNIARD L.Basic theory of the magnetotelluric methods of geophysical prospecting［J］.Geophysics，1953（18）：605－635.

［3］COGGON J H.Electromagnetic and electrical modeling by the finite element method［J］.Geophysics，1971，36（2）：132－151.

［4］朱伯芳.有限单元法原理与应用［M］.北京：水利出版社，1979.

［5］周熙襄，钟本善，江东玉.点源二维电阻率法有限差分正演计算［J］.物探化探电子计算技术，1983，5（3）：34－38.

［6］罗延钟，万乐.二维地形不平条件下外电场的有限差分模拟［J］.物探化探计算技术，1984，6（4）：15－26.

［7］罗延钟，万乐.用数值模拟方法构组保角变换坐标网［J］.物探化探计算技术，1986，8（1）：23－31.

［8］陈小斌.有限元直接迭代算法［J］.物探化探计算技术，1999，21（2）：165－171.

［9］徐世浙.地球物理中的有限单元法［M］.北京：科学出版社，1994.

［10］陈乐寿.有限元法在大地电磁测深正演计算中的应用与改进［J］.石油物探，1981，20（3）：84－103.

［11］陈乐寿，孙必俊.有限元法在大地电磁测深中正演计算的改进［J］.石油地球物理勘探，1982，17（3）：69－72.