柳建新 ，汤文武 ，童孝忠

(1. 中南大学 地球科学与信息物理学院，湖南 长沙，410083；2. 有色资源与地质灾害探查湖南省重点实验室，湖南 长沙，410083)

2.5 维数值模拟指的是源为三维的、正演模型为二维的地球物理数值模拟问题[1−3]。众所周知，2.5维数值模拟问题有着重要的实用价值，这是因为将三维问题简化为2.5维模型后只需要对剖面而不是整个体积空间进行离散处理，这样就大大减小了矩阵的尺寸[4]。

2.5维频率域可控源电磁法(frequency domain controlled-source electromagnetic, FCSEM)的有限元正演模拟，国内外都进行了研究。Unsworth等[5]将一次场与二次场分离，利用有限元对二维地电模型下水平电偶源的电磁场响应进行了模拟计算；Mitsuhata[6]借鉴地震模拟中的伪delta函数模拟电偶源，只需一次性计算总场，并利用等参有限元实现了带地形的2.5维频率域电磁正演模拟；底青云等[7−8]对2.5维频率域电磁进行了有限元模拟并给出了一组波数，对层状模型的模拟计算效果较好；薛东川等[9−10]对2.5维频率域电磁的吸收边界条件进行了讨论，发现吸收边界条件能够有效压制边界处的反射。在此，本文作者在前人研究基础上，重点研究基于双二次插值的FCSEM有限元正演模拟，并对波数选取进行研究，给出一组兼顾计算速度及模拟精度的波数。

1 有限元方程的推导

图1所示为2.5维FCSEM正演模拟示意图：其y方向为走向方向，水平电偶源Idy沿走向方向布置。电磁法勘探是基于电磁场理论的一类地球物理勘探方法，2.5维频率域可控源电磁(FCSEM)也应满足麦克斯韦方程。假设时谐因子为eiωt，则麦克斯韦方程表示如下[3]：

式中：E为电场强度；H为磁场强度；Js为外加电偶源；阻抗率=iμω；导纳率=σ+εω。对式(1)和(2)中的所有电场分量及磁场分量都进行傅里叶变换得：

可以得到如下的微分控制方程[6]：

图1 2.5-D FCSEM模拟示意图Fig. 1 Simulation schematic diagram of 2.5-D FCSEM

首先利用有限元方法求解式(4)和(5)可得到傅氏域中的电场分量y及磁场分量y，再利用式(6)～(9)求得其余4个分量，最后通过傅氏逆变换可求得真实电磁场。这里直接给出由伽辽金有限元方法推导得到的有限元方程如下[6]：

2 有限元求解过程

先对模拟区域进行剖分。由于同等计算量下双二次插值一般较双线性插值精度高，因此，本文采用双二次插值以提高计算精度；接着对各个单元进行分析、计算，得到各个小型刚度矩阵；最后，对所有单元总体集成得到总体刚度矩阵[11−15]。

2.1 网格剖分及插值

首先在研究区域内进行矩形单元剖分。为了提高计算精度及减小利用式(6)～(9)求解时带来的误差，本文采用精度较高的双二次插值进行模拟。图2所示为坐标变换关系及局部节点编号，其中图 2(a)所示为母单元，图2(b)所示为子单元；a和b分别表示子单元的宽度及高度；(x0，z0)为矩形单元中心点的坐标。文中采用的插值函数表示如下：

图2 矩形单元剖分示意图Fig. 2 Schematic diagram of rectangle meshing element

子单元与母单元之间的坐标变换关系式表示如下：

2.2 单元分析

首先对式(10)进行单元分析。

记U1=(u1，u2，…，u8)T，V1=(v1，v3，…，v8)T，则

式(10)中左边第1项：

式(10)中左边第2项：

式(10)中左边第3项：

式(10)中左边第4项：

式(10)中左边第5项：

再对式(11)进行单元分析。

式(11)中左边第1项：

式(11)中左边第2项：

式(11)中左边第3项：

式(11)中左边第4项：

式(11)中左边第5项：

观察式(15)～(24)，借助符号计算工具maple可求得：

2.3 总体合成

将式(15)～(19)中各项累加，从而可将式(10)化为

其中：J代表电流源项。将式(20)～(24)中各项累加，从而可将式(11)化为

再对式(30)及(31)化简，令

则式(30)和(31)可记为：

将式(36)等价为

将式(37)中左边的刚度矩阵记作

经分析可知：M矩阵为复对称矩阵。求解线性方程组(38)即可得到傅氏域的 2个电磁场分量，再利用式(6)～(9)求得傅氏域中的其余4个电磁场分量；最后，利用傅里叶逆变换求得各个电磁场分量。

3 波数的选取

2.5 维数值模拟是先在波数域里计算电磁场，再通过傅里叶逆变换求得空间域的电磁场。由于在数值模拟时不可能取无穷多个波数来还原真实的电磁场，因此，需要选取若干个合适的波数以尽可能、快速地求得真实的电磁场。如何选取波数及确定波数的选取范围便是一个关键性问题。

本文通过实际计算及分析，选取一组兼顾计算速度与精度的波数。在收发距为4 km，频率f为2 Hz和1 024 Hz时，分别进行正演模拟计算并绘制了及随波数的变化曲线图，如图3所示。其中，图3(a)所示为随波数的变化曲线图，图3(b)所示为随波数的变化曲线图。从图3(a)可知：当波数较小时，基本保持不变；随着波数继续增大时，开始下降；当波数约为0.01时到达波谷，而后又稍有回升。由图 3(b)可知：当波数较小时，很小；随着波数增大，开始先上升，到达波峰后以较陡倾角下降；同样地，当波数约为0.01时到达低谷，而后又稍有回升，但变化曲线没有规律。据文献[6]，波谷右侧的曲线是不可靠的。事实上，曲线到达波谷后继续往下递减，且由于此时取值已经很小，对傅里叶逆变换中的积分结果的影响可以忽略，因此，只需要在波谷左边部分选取足够的能够精确计算积分的一组波数即可。本文通过不断进行正演计算及对比分析，选取由20个波数构成的一组波数，如表1所示。

图3 傅氏域中电场和磁场随波数变化示意图Fig. 3 Schematic diagram of electric field and magnetic field in Fourier domain of changing with wave number

表1 波数值选取Table 1 Choice of wavenumbers

事实验证，这一组波数能够有效且快速、精确计算2.5维FCSEM的电磁场响应。本文后续的工作都是基于这一组波数进行的。

4 数值算例

4.1 算法验证

为了验证本文算法的正确性，对1个电阻率为100Ω·m的均匀半空间模型进行正演模拟计算及与解析解的对比分析。图4所示为解析解与数值解对比结果图(其中，x轴代表的是收发距)。由图4可知：电磁场分量数值解与解析解相吻合，且频率低时精度更高。计算结果表明：对于电场分量，当f=2 Hz时，最小相对误差为 0.106%，最大相对误差为 0.639%，平均相对误差为0.460%；当f=1 024 Hz时，最小相对误差为0.451%，最大相对误差为 1.763%，平均相对误差为0.920%。对于磁场分量，当f=2 Hz时，最小相对误差为 0.016%，最大相对误差为 0.513%，平均相对误差为0.220%；当f=1 024 Hz时，最小相对误差为0.545%，最大相对误差为 1.921%，平均相对误差为 1.028%。可见，完全满足正演计算的要求。

4.2 模型计算

4.2.1 低阻体模型

图5所示为 1个低阻体模型的示意图，其中，x方向标出的刻度代表垂直收发距。图6和图7所示为正演模拟结果。其中，图6(a)和图6(b)所示分别为TE模式下的视电阻率及相位等值线图，这2个图中都反映出了异常，且异常的中心位置与理论模型一致，异常都有往右下方拉长的趋势，这正是场源阴影效应的影响；视电阻率等值线图异常最低为30 Ω·m，接近异常体的真实值。图7(a)和图7(b)所示分别为TM模式下的视电阻率及相位等值线图，这2个图中也都反映出了异常，且相比 TE模式异常垂向拉长明显。这是由于TM模式下测定的是电场的法向分量，界面积累电荷效应使得受静态效应明显，造成垂向拉长，这提高了方法的灵敏度，有利于发现小规模低阻体，故TM模式发现异常的能力比TE模式的强。

图4 均匀半空间模型解析解与数值解对比结果Fig. 4 Comparison between analytic solution and numerical solution of homogenous half space

图5 低阻体模型示意图Fig. 5 Schematic diagram of low resistivity body

图6 TE模式下低阻体模拟结果Fig. 6 Schematic simulation result of low resistivity body of TE mode

4.2.2 高阻体模型

图7 TM模式下低阻体模拟结果Fig. 7 Simulation results of low resistivity body of TM mode

图8 高阻体模型示意图Fig. 8 Schematic diagram of high resistivity body

图8 所示为1个高阻体模型的示意图，图9和图10所示为正演模拟结果。其中，图9(a)和(b)所示分别为TE模式下的视电阻率及相位等值线图，这2个图中都反映出了异常，视电阻率等值线图中异常最高为130 Ω·m，而真实异常为500 Ω·m，说明频率域可控源电磁对高阻体反映不灵敏，这是电磁法类勘探方法普遍存在的一个问题。图 10(a)和 10(b)所示分别为 TM模式下的视电阻率及相位等值线图，这2个图中都反映出异常，视电阻率等值线图中异常最高为200 Ω·m，说明TM模式较TE模式对高阻体反映较灵敏。

图9 TE模式下高阻体模拟结果Fig. 9 Simulation results of high resistivity body of TE mode

图10 TM模式下高阻体模拟结果Fig. 10 Simulation results of high resistivity body of TM mode

5 结论

(1) 通过均匀半空间模型下低频及高频时的有限元数值解与解析解对比，结果显示有限元数值解与解析解基本吻合，从而从一定程度上验证了基于双二次插值的2.5维FCSEM正演模拟算法的正确性。

(2) 本文2.5维FCSEM有限元数值模拟结果能够正确反映二维模型的浅部异常及构造，进一步验证了该算法的可靠性，可以有效地模拟2.5维地电模型。

(3) TM模式下低阻及高阻异常体响应均有往下拉长的趋势，对低阻及高阻异常体反应较TE模式明显，且横向分辨能力更强。

[1]底青云, 王若. 可控源音频大地电磁数据正反演及方法应用[M]. 北京: 科学出版社, 2008: 74−84.DI Qingyun, WANG Ruo. Forward modeling and inversion of CSAMT data and its application in earth[M]. Beijing: Science Press, 2008: 74−84.

[2]汤井田, 何继善. 可控源音频大地电磁法及其应用[M]. 长沙:中南大学出版社, 2005: 9−33.TANG Jingtian, HE Jishan. Audio-frequency of controlled source electromagnetic method and its application in earth[M].Changsha: Central South University Press, 2005: 9−33.

[3]米萨克 N 纳比吉安. 勘查地球物理电磁法. 1卷: 理论)[M].赵经祥, 王艳君, 译. 北京: 地质出版社, 1992: 330−333.Nabighian M. Electromagnetic methods in applied geophysics:Ⅰ. Theory[M]. ZHAO Jingxiang, WANG Yanjun, trans. Beijing:Geological Publishing House, 1992: 330−333.

[4]熊彬. 关于瞬变电磁法2.5维正演中的几个问题[J]. 物探化探计算技术, 2005, 28(2): 124−128.XIONG Bin. Several problems related to 2.5 dimensional forward modeling of transient electromagnetic method[J].Computing Techniques for Geophysical and Geochemical Exploration, 2005, 28(2): 124−128.

[5]Unsworth M J, Travis B J, Chave A D. Electromagnetic induction by a finite electric dipole source over a 2-D earth[J].Geophysics, 1993, 58(2): 198−214.

[6]Mitsuhata Y. 2-D electromagnetic modeling by finite-element method with a dipole source and topography[J]. Geophysics,2000, 65(2): 465−475.

[7]底青云, Unsworth M, 王妙月. 复杂介质有限元法2.5维可控源音频大地电磁法数值模拟[J]. 地球物理学报, 2004, 47(4):723−730.DI Qingyun, Unsworth M, WANG Miaoyue. 2.5 D CSAMT modeling with the finite element method over 2D complex earth media[J]. Journal of Geophysics, 2004, 47(4): 723−730.

[8]底青云, Unsworth M, 王妙月. 有限元法2.5维CSAMT数值模拟[J]. 地球物理学进展, 2004, 19(2): 317−324.DI Qingyun, Unsworth M, WANG Miaoyue. 2.5 D CSAMT modeling with finite element method[J]. Progress in Geophysics,2004, 19(2): 317−324.

[9]薛东川, 戴世坤. 频率域2.5维电磁测深有限元模拟中的吸收边界条件[J]. 中国石油大学学报(自然科学版), 2008, 32(6):57−61.XUE Dongchuan, DAI Shikun. Absorbing boundary condition for simulating 2.5-D electromagnetic sounding in frequency domain by finite element method[J]. Journal of China University o f Petroleum (Natural Science Edition), 2008, 32(6): 57−61.

[10]王华军, 罗延钟. 中心回线瞬变电磁法2.5维有限单元算法[J].地球物理学报, 2003, 46(6): 855−862.WANG Huajun, LUO Yanzhong. Algorithm of a 2.5 dimensional finite element method for transient electromagnetic with a central loop[J]. Journal of Geophysics, 2003, 46(6):855−862.

[11]TONG Xiaozhong, LIU Jianxin, XIE Wei, et al.Three-dimensional forward modeling for magnetotelluric sounding by finite element method[J]. Journal of Central South University of Technology, 2009, 16(1): 136−142.

[12]柳建新, 蒋鹏飞, 童孝忠, 等. 不完全 LU 分解预处理的BICGSTAB算法在大地电磁二维正演模拟中的应用[J]. 中南大学学报(自然科学版), 2009, 40(2): 484−491.LIU Jianxin, JIANG Pengfei, TONG Xiaozhong, et al.Application of BICGSTAB algorithm with incomplete LU decomposition preconditioning to two-dimensional magnetotelluric forward modeling[J]. Journal of Central South University (Science and Technology), 2009, 40(2): 484−491.

[13]TANG Jingtian, WANG Feiyan, REN Zheng-yong, et al. 3-D direct current resistivity forward modeling by adaptive multigrid finite element method[J]. Journal of Central South University of Technology, 2010, 17(3): 587−592.

[14]刘长生, 汤井田, 任政勇, 等. 基于非结构化网格的三维大地电磁自适应矢量有限元模拟[J]. 中南大学学报(自然科学版),2010, 41(5): 1855−1859.LIU Changsheng, TANG Jingtian, REN Zhengyong, et al.Three-dimension magnetotellurics modeling by adaptive edge finite-element using unstructured meshes[J]. Journal of Central South University (Science and Technology), 2010, 41(5):1855−1859.

[15]张继锋, 汤井田, 王烨, 等. 基于预处理共轭梯度的大地电磁快速正演[J]. 中南大学学报(自然科学版), 2010, 41(5):1877−1882.ZHANG Jifeng, TANG Jingtian, WANG Ye, et al.Magnetotellurics fast forward based on preconditioning conjugate gradient[J]. Journal of Central South University(Science and Technology), 2010, 41(5): 1855−1859.