0 引 言

可控源音频电磁方法(CSAMT)的野外装置分为两个子系统，即发射子系统和接收子系统，这两个子系统相距几公里.接收电偶极和发射源的布设方向一致，接收由发射源发出的、耦合了地电信息的水平电场信号，同时接收的还有地面水平磁场，磁场的方向与电偶极的方向垂直.用观测到的磁场和电场计算出卡尼亚视电阻率，卡尼亚视电阻率便是最终用于处理解释的观测资料之一.

从CSAMT的工作装置看出，CSAMT本身是三维的方法.探测的目标体如金属矿脉、含水溶洞等都具有三维特性，这些因素决定了对CSAMT方法进行三维研究的必要性.

三维电磁正演模拟早在上世纪70年代便已有学者研究(Raiche，1974).模拟方法有积分方程法(Wannamaker et al., 1984; Wannamaker et al., 1991; Xiong，1992; 殷长春，1994；Zhdanov et al., 2000; Garbon et al., 2002)，有限差分法(Adhidjaja et al., 1989; Mackie et al., 1993; Wang et al., 1996; 谭捍东，2000；沈金松，2003)和有限元方法(FEM).由于有限元方法相对于有限差分方法和积分方程法能更有效地处理三维复杂问题，所以许多学者致力于有限元在电磁法中的应用研究并取得了较大进展(Coggon，1971; Pridmore et al., 1981; 徐世浙，1994；金建铭，1998；阎述等，2000；黄临平等，2002；王若等，2007；汤井田等，2007；张继锋等，2009；徐志锋等，2010).

作者曾对CSAMT三维三分量(3D3C)有限元正演做了初步研究(王若等，2007)，由于其计算量大，花费时间长，不利于应用于反演中，所以本次研究以提高计算精度，缩短计算时间为目的.

边界条件在有限元研究中必不可少.一阶吸收边界条件需要计算区域范围很大(王若等，2006)，剖分网格较多，虽可通过加大网格间距来减小网格数量，但无疑会降低计算精度.本文试图应用边界场值不为零的第一类边界条件，借此来减小研究区域，从而在网格剖分数不变的情况下，提高计算精度.

由于计算区的边界场值未知，需先用其它的方法计算得到.文献(徐志锋等，2010)做了这种尝试，作者用标量势与磁矢量位的3D有限元方法计算了3D CSAMT电场，取得了较大进展.然而该文作者在计算边界上的标量势和磁矢量位时，采用了均匀半空间的解析公式来确定，然后将计算得到的边界上的值作为第一类边界条件.这种方法只有当计算区的电性结构缓变时，解才较可靠，从普通意义上说尚有一定的局限性.本文采用一种更为普通适用的方法求取边界场值，即用3D积分方程法中层状介质的格林函数计算边界场值，应用这种方法的好处是可得到层状介质的边界条件，若获得成功，可以模拟层状介质背景中的异常体.

除此之外，本文还将三维三分量退化为三维单分量(3D1C)来研究.由于有限元结点上自由度的减少，计算速度大大提高.一般来说，三维CSAMT方法需要研究三维三分量，然而，当电性介质的变化不是很大时，电场分量之间的耦合可以忽略，从而理论场的正演可退化为三维单分量问题.

目前国内常用的是标量CSAMT方法，所以三维单分量的研究有一定的实际意义. 1 三维单分量电磁场有限元正演 1.1 三维单分量电磁场边值问题

当电性结构缓变时，电场的散度为零.从麦克斯韦方程可推导出

在边界上满足为

其中，E_x为电场的x分量，σ是介质电导率，ω是圆频率，μ是介质的磁导率，J_cx是外加源，J_cx=Idsδ(x)δ(y)δ(z)，其中I为外加电流强度，ds为偶极源长度，δ()为δ函数.Γ为人工截断边界，E_xb是边界上背景电场.

1.2 第一类边界条件

方程(2)为第一类边界条件.除此之外，有限元分析中还可应用边界值为零的第一类边界条件、场的导数为零的第二类边界条件和第三类边界条件(即一阶吸收边界条件).但应用这些边界条件，需将边界设置得离一次源较远，而边界设置得越远，研究区域就会越大.对于有限元方法来说，研究区域扩大就意味着剖分的单元数量、内存消耗量和计算量都会大幅增加.如果想减少计算量和内存消耗量，就要使研究区域尽量小.应用边界条件不为零的第一类边界条件，既使边界位置不会离源太远，又能保证边界上的值准确可靠.

作者利用在Utah大学访问学习的机会，采用CEMI组研发的3D积分方程法，算出了层状介质(可将空气和大地均匀半空间介质各看作一层，也可将地下介质分为多层)边界上的电磁场三分量场值.用计算的结果作为第一类边界条件应用到三维有限元正演分析中.

1.3 有限元分析

将研究区域用六面体划分成一系列的小单元，(ξ，ζ，η)为小单元内的局域坐标，(x，y，z)为总体坐标；设单元边长分别为a，b，c(如图 1所示)，各点在两种坐标系下的对应关系为

图 1

图 1 六面体单元的编号及坐标 Fig. 1 Number and coordination of hexahedron element

采用局域坐标后，函数的积分变为

设小单元内的电场可用角点的电场线性表示，如果定义线性插值时所用的形函数N_i(i =1，2，3，…，8)为

则单元中的电场u可用形函数N_i表示为

其中，u_i是电场在单元的八个角点的待定场值.

由方程(1)，定义余量为

以N_i 为权函数，使余量的加权积分为零：

根据(6)式，可将单元内的E_x用形函数与结点处的电场值表示为

E _x= N T E，(9)

其中： N T=(N₁ N₂ N₃ N₄ N₅ N₆ N₇ N₈)，

E =(E₁ E₂ E₃ E₄ E₅ E₆ E₇ E₈)T.

将(9)式代入(7)式，有

将(10)式代入(8)式，经整理可写出有限元方程

KE=P，(11)

其中，

(14)式中，当单元边界为内边界时，在单元边界S上的面积分与相邻不同单元的贡献相互抵消，当单元边界为外边界的一部分时，其贡献不能取为零，而是通过外边界上已知场值计算得到.

对于(11)式右侧的源项，有

集成后的有限元方程为

通过解方程组(16)可得到电场E_x，磁场可由麦克斯韦尔方程计算，磁场和卡尼亚视电阻率的公式为

其中，ρ_s为视电阻率，H_y为y方向的磁场，(17)式的可行性将单独讨论.

1.4 程序可靠性检验

为了检验程序的可靠性，用电阻率是100 Ωm的均匀半空间模型来做试验.

将坐标原点置于发射源中心，x，y轴位于地平面上，x轴东西向放置，东为正向，z轴垂直地面向下，x，y，z轴的方向遵循右手螺旋法则.发射偶极源的长度为500 m，电流为20 A，沿x轴方向布设.将研究区域剖分成20×24×25的网格，相关参数见表 1.表中计算区域指网格剖分所覆盖的区域，所给各值为边界位置.勘探区域是指展示计算结果的区域，相当于野外的勘探区.从表中可以看出，边界离源的距离较近，最远不过6600 m，边界离勘探区域则更近.

表 1(Table 1)

表 1 计算参数表 Table 1 Calculation parameters 网格距 计算区域(m) 勘探区域(m) x方向 100 m，边界附近较大 -3600~3600 -1000~1000 y方向 源附近100 m，其余300 m 0~6600 4000~6000 z方向 空气层200 m，地下不均匀 -1000～3700 地面z=0

表 1 计算参数表 Table 1 ExamplesCalculation parameters

图 2为在地面上测点(100，4800，0)处用有限元计算得到的电磁场、视电阻率与解析解的对照曲线.在这三幅图中，横轴均为观测频率，纵轴分别为电场幅值、磁场幅值和视电阻率值.实线为解析解曲线，离散点为相应频率用有限元计算得到的结果.从图中可以看出电磁场的有限元结果和解析结果比较接近，但是并没有完全吻合，视电阻率和解析结果的吻合程度相对较好.这是因为磁场是由电场计算得到的，二者有着相近的数值模拟误差，而视电阻率是用电场和磁场的比值得到的，可以部分消除掉数值模拟误差的影响.电场、磁场和视电阻率与解析解的相对百分比误差分别为9.47%，10.76%和5.15％.以上信息说明正演程序基本正确.在文献(王若等，2007)中，由三维三分量有限元计算的视电阻率误差为13.23%，可见，通过应用场值不为零的第一类边界条件，不但缩小了计算区域，而且提高了计算精度.

图 2

图 2 有限元结果与解析结果的对比 实线为解析结果，离散点为有限元计算结果 (a)电场;(b)磁场;(c)视电阻率. Fig. 2 Comparison of FEM results and analytic solution Solid line is analytic solution，disperse point is FEM results (a)Electric field;(b)Magnetic field;(c)Apparent resistivity.

本文展示了两个数值模型的有限元模拟结果.第一个是单个低阻模型，第二个是平行x轴方向排列的两个低阻模型.

这两个模型的坐标系统、发射源的位置和有限元网格的剖分方式与均匀半空间模型中所用的相同.每个模型计算共发射13个频率的场，频率范围为2¹~2¹³ Hz.

1.5.1 单个低阻模型

(1)模型设计

在100 Ωm 的均匀半空间中赋存着一个10 Ωm的低阻异常体，低阻异常体在x，y，z方向的尺寸分别为200 m，300 m，200 m.异常体的上界面位于地下200 m处.其中心在地面的投影点位于y轴上，且与源中心的距离为4950 m.模型示意图如图 3所示.

图 3

图 3 单个低阻体模型示意图 均匀半空间：100 Ωm，低阻体：10 Ωm，空气层：1011 Ωm AB=1000 m，I=10 A，方框为正演结果显示区域，相当于测区，图中给出了测区两个角点的坐标，下同. Fig. 3 Sigle low resistivity anomaly model half space: 100 Ωm，low resistivity anomaly: 10 Ωm，air: 1011 Ωm. AB=1000 m，I=10 A. the rectangle is survey area，two corner coordirations are show in the figure，same in the following figure.

(2)正演结果

图 4为异常体所在区域的电磁场和视电阻率在观测面上(地面上)的有限元分析结果，每一幅图的横坐标为x方向的测点，纵坐标为y方向的测点.图 4包括三行，每一行代表一个频率的计算结果，共给出三个频率的计算结果，从上到下依次对应频率8192 Hz、128 Hz和4 Hz.每一行图件从左到右的物理量分别为电场、磁场和视电阻率.

图 4

图 4 地面上单个低阻模型的有限元正演结果 (a)(b)(c)分别为f=8192 Hz时的电场、磁场和视电阻率正演结果; (d)(e)(f)分别为f=128 Hz时的电场、磁场和视电阻率正演结果 (g)(h)(i)分别为f=4 Hz时的电场、磁场和视电阻率正演结果. Fig. 4 Forward modeling result of single low resistivity anomaly (a)(b)(c)Are electric field，magnetic field and apparent resistivity，respectively，at 8192 Hz; (d)(e)(f)Are electric field，magnetic field and apparent resistivity，respectively，at 128 Hz; (g)(h)(i)Are electric field，magnetic field and apparent resistivity，respectively，at 4 Hz.

当频率是8192 Hz时，在电磁场和视电阻率平面等值线图上看不到异常信息.这是因为高频时电磁波在浅层传播，主要体现的是浅层均匀背景场的信息.从视电阻率值也能看出，电阻率值在96~101 Ωm之间，非常接近背景值100 Ωm.当频率为128 Hz时，电场和磁场都出现了明显异常，异常在视电阻率等值线平面图上显示得更加清楚，异常中心的电阻率值在76~81 Ωm之间，比背景值低，但与真值10 Ωm尚存在一定差异.至4 Hz时，电场平面图上没有明显的异常反映，但在磁场图上还有较弱的异常存在，电阻率平面图上虽有异常显示，但与128 Hz相比，异常有弱化现象.这是因为电磁场具有体积效应，频率越低，地面观测到的电磁场所携带的信息是较大体积内地电介质的综合反映，所以4 Hz时异常信息减弱，但异常的信息还未被周围介质的信息完全平滑. 1.5.2 双低阻模型_沿与x轴平行方向排列

(1)模型设计

在100 Ωm的均匀半空间中有两个10 Ωm的低阻异常体，低阻异常体的尺寸大小与埋深都和单个异常体模型中的相同.两个异常体中心连线的中点在地面的投影位于y轴上，与y轴交于4950 m处.两个异常体中心在x轴上的投影分别为-300 m和300 m.模型示意图如图 5所示.

图 5

图 5 双低阻体模型(平行于x轴)示意图 Fig. 5 Double low resistivity anomaly model-parallel to the x-axis

(2)正演结果

图 6给出了地面上异常体所在区域的电磁场和视电阻率的观测结果，各图的意义同图 4.

从图中可以看出两个异常体在128 Hz与4 Hz时均有显示，除此之外，还可看出和图 4相似的信息，即：频率是8192 Hz时，电磁场和视电阻率图看不到低阻异常体产生的异常信息；当频率是128 Hz时，电场和磁场都出现了异常，而且很明显是两个异常，视电阻率平面图上显示得更加清楚；至4 Hz时，电场平面图上没有明显的异常反映，但在磁场图上还有较弱的异常存在，电阻率平面图上也有弱异常显示.其中的机理不再阐述.

图 6

图 6 地面上双低阻模型(平行于x轴)的有限元正演结果 (a)(b)(c)分别为f=8192 Hz时的电场、磁场和视电阻率正演结果; (d)(e)(f)分别为f=128 Hz时的电场、磁场和视电阻率正演结果; (g)(h)(i)分别为f=4 Hz时的电场、磁场和视电阻率正演结果. Fig. 6 Forward modeling result of double low resistivity anomaly parallel to x-axis (a)(b)(c)Are electric field，magnetic field and apparent resistivity，respectively，at 8192 Hz; (d)(e)(f)Are electric field，magnetic field and apparent resistivity，respectively，at 128 Hz; (g)(h)(i)Are electric field，magnetic field and apparent resistivity，respectively，at 4 Hz.

在式(17)中，磁场的计算只应用了电场的一个分量，而在完整公式中，磁场由两个分量得到为(王若，2006)

其中，E_x、E_z分别是电场的 x分量和z分量.

通过计算发现，在地面上z 分量电场幅值非常小，这是因为电场在远区是垂直入射到地中的，电场只存在和入射方向垂直的水平分量.三维积分方程的计算结果也证明，在地面上垂直电场E_z的量级为10^-11，相对于E_x的量级10^-5和E_y的量级10^-6而言，量值小到可以忽略.经过计算电场的导数 E_x z和E_z x，发现对于本文中所用的模型来说，E_x z 的量级在8192 Hz为10^-7，在128 Hz为10^-8，在4 Hz为10^-9，而 E_z x 的量级从高频到低频均为10^-12，相对于 E_x z 也是可以忽略的.鉴于此，当用(19)式计算磁场时，可以简化为只用x方向的电场近似求出y方向的磁场，即公式(17)是可行的.说明针对本文所用的模型可以把三维三分量电磁场的有限元分析退化为三维单分量电磁场来研究.

但值得引起注意的是：当把低阻模型顶部置于地面时，E_z x 的量级从高频到低频均为10^-8，相对于 E_x z 来说便不可忽略. 3 结 论

通过本文的研究，可以得到以下认识：

(1)应用场值不为零的第一类边界条件使三维三分量的有限元计算结果有了很大改善.将这种边界条件应用于 三维单分量有限元分析中，使得在离源5 km左右得到的视电阻率与解析解的误差达到5%左右.说明所用边界条件有效，应用场值不为零的第一类边界条件后不但缩小了计算区域，而且提高了计算精度.

(2)用三维单分量来做有限元分析，相对于三维三分量电场正演减少了每个结点上自由度，从而可以将计算效率大大提高，继而可将3D正演应用于反演中.

(3)在勘探区，只由电场单个分量计算磁场的公式只适用于异常体位于地下的情形，当异常体出露于地表时，磁场的计算应该仍用两个电场分量计算得到.

(4)通过数值模拟结果可以看出，有限元分析结果明显反映出了电磁波的传播规律及在地中传播时固有的体积效应，表明CSAMT单分量有限元分析的计算结果正确.

虽然本文得到了一些可喜结果，但计算误差有待进一步降低，特别是在接近于一次源的区域误差还较大，算法仍需进一步改进.

致 谢 感谢在犹他大学访问时，CEMI学科组给予的帮助.