大地电磁测深法(MT)是利用大地中频率范围很宽广泛分布的天然变化交变电磁场研究地球电性结构的一种地球物理勘探方法.由于它不用人工供电，成本低，工作方便，所以在实际中被广泛应用，因此对于大地电磁的正演模拟也就会非常重要.

本论文的研究内容就是应用矢量有限元方法进行三维大地电磁正演模拟研究，这种有限单元法使用矢量插值函数，将要求的场值建立在单元的棱边而不是节点，这种基于单元棱边的矢量基函数，在保证边界上磁场连续性的同时也保证了正切电场的连续性，因此采用矢量有限元的方法对于强加边界条件非常容易，另外由于矢量基函数的散度为零，因此也会避免节点有限元中伪解的出现.矢量有限元方法在大地电磁正演方面取得了一定的成果，Xueming Shi等人研究了结合散度校正的三维大地电磁矢量有限元正演模拟(Shi，2004)；随后，韩国学者Myung等人研究了利用矢量有限元方法进行带地形三维大地电磁的正演(Nam et al.，2007)；刘长生、任正勇、汤井田等人应用了基于非结构化四面体网格的H型自适应矢量有限单元法对大地电磁法进行了三维正演(刘长生等，2008)；中南大学的王烨博士毕业论文研究了基于矢量有限元的高频大地电磁法三维正演(王烨，2008).

本文在前人工作基础上，以矢量有限元方法作为数值模拟工具，采用结构化六面体网格，在目标区网格均匀剖分，然后向两边逐渐扩大网格间距.直接从电场的双旋度方程出发，利用伽辽金方法推导了关于电场的变分形式；在传统的三维大地电磁边界条件基础上，给求解区域四个垂直侧面加上第一类边界条件；把不完全Cholesky分解和双复共轭梯度法相结合求解方程，能够保证计算精度的同时，高效、稳定的收敛；最后对典型地电模型进行了正演计算，并与解析解进行了对比，精度满足实际要求.

设初始大地电磁场是平面波场，初始电场E的偏振方向沿x轴，如图 1所示，选取足够大的六面体区域Ω，使地下三维异常体产生在区域Γ边界上的异常电磁场为零，则电磁场的边界条件：

图 1

Fig. 1

图 1 三维电性结构和坐标示意图 Fig. 1 Three-dimensional electrical schematic structure and coordinate

(1)在区域顶界面(ABCD)上，E_x=1，E_y=0，E_z>=0 ；

(2)在四个垂直的边界面上，电磁场传播方向垂直向下，与边界面的法线方向垂直，E×H⊥Γ；

(3)在区域底界面(EFGH)上，电磁场按指数规律向下传播，E_x=Ce^－χz，其中C是常数，χ=，其中σ是EFGH面以下的均质电导率.

以上是传统的三维大地电磁边界条件，但是在后面模型验证的过程中发现利用矢量有限元方法进行大地电磁正演时，如果计算区域四个垂直侧面不加边界条件计算精度低，正演速度慢，低频时甚至不收敛，针对此情况，在上面边界条件的基础上给定新的边界条件，前后面：E_z=0，左右面：E_y=0，E_z=0，底界面：E _y=0.

以e^－iwt表示谐变场的时间因子.导电介质中的麦克斯韦方程组变为

w是角频率，u(4π×10^－7H/m)是磁导率，ε(8.85×10^-12 F/m)是介电常数，σ是介质电导率.对第一式两边求旋度，将第二式代入(1)，得到关于电场的双旋度方程为

三维情况下，E=E_xe_x+E_ye_y+E_ze_z，由广义变分原理可以得到(3)式的泛函为

(4)式便为基于电场矢量波动方程的大地电磁变分公式，第三项是剖分单元边界落在计算区域底界面时的边界条件公式.

本文将研究区域用六面体单元进行剖分，向右为X方向，向前为Y方向，向下为Z方向，如图 2，为了实际模拟的需要将中心区域等距离剖分并且需要剖分的尺寸较小，远处可适当加大剖分单元的尺寸.每个单元和每个单元的各条棱边均有自己的编号，x、y、z方向的边长分别为A，B，C，其中心坐标为(x_c^e，y_c^e，z_c^e)，经过推导，每个六面体单元中的内部电场分量可用六面体的12条棱边的场值分量通过插值求出：

图 2

Fig. 2

图 2 单元中的棱边号次序 Fig. 2 The edge order number in the cell edge

将(6)式结合(5)式代入(4)式中，可以求得任意一个子单元内部的矩阵分块方程组：

从矢量有限元数值模拟中计算出的是六面体棱边上的电场值，我们还必须计算出地表所求处的视电阻率，视电阻率的表达式为

本项研究基于矩形六面体剖分，整个模型区域(63250 m×63250 m×168850 m)沿三个方向剖分块数为N_x×N_y×N_z=45×45×42，其中空气层为10层.

首先给出了均匀半空间模型，大地的电阻率为100 Ω·m，表 1为数值模拟结果：

表 1(Table 1)

表 1 矢量有限元解与解析解对比 Table 1 Comparison of vector finite element solution with analytical solution 频点 电阻率(ρx)/Ω·m 相位(φ)/(°) 矢量有限元解 解析解 误差(%) 矢量有限元解 解析解 误差(%) 10000 97.42975 100 2.570254 46.46353 45 3.252292 5000 98.29741 100 1.702594 45.56092 45 1.246496 1000 99.40596 100 0.59404 45.04157 45 0.092373 500 99.70895 100 0.291051 45.06803 45 0.151184 100 99.95241 100 0.047592 45.02291 45 0.050913 50 100.0072 100 0.00723 45.00362 45 0.008048 10 99.83015 100 0.169849 44.95041 45 0.110205 5 99.75044 100 0.249558 44.91745 45 0.183441 1 99.38591 100 0.614092 44.93628 45 0.141609 0.5 99.30442 100 0.695584 44.92528 45 0.166046 0.1 99.31773 100 0.682268 45.01517 45 0.033717 0.05 99.23363 100 0.76637 44.93809 45 0.137581 0.01 96.84162 100 3.15838 44.41066 45 1.309641 0.005 92.59314 100 7.406861 47.59291 45 5.762017 0.001 216.7257 100 116.7257 47.27401 45 5.053361

表 1 矢量有限元解与解析解对比 Table 1 Comparison of vector finite element solution with analytical solution

这里分别给出了两个三层模型(H型和K型)，在不同频点计算了电阻率值和相位值，与解析解进行对比.

模型1(K型)：水平均匀层状模型，分三层，第一层电阻率为100 Ω·m，深度为300 m；第二层电阻率为1000 Ω·m，深度为600 m；基底电阻率为10 Ω·m.

模型2(H型)：水平均匀层状模型，分三层，第一层电阻率为1000 Ω·m，深度为300 m；第二层电阻率为10 Ω·m，深度为600 m；基底电阻率为100 Ω·m.

图 3

Fig. 3

图 3 K型模型视电阻率曲线 Fig. 3 Curve of resistivity for K layer earth

图 4

Fig. 4

图 4 K型模型相位曲线 Fig. 4 Curve of phase for K layer earth

图 5

Fig. 5

图 5 H型模型视电阻率曲线 Fig. 5 Curve of resistivity for H layer earth

图 6

Fig. 6

图 6 H型模型相位曲线 Fig. 6 Curve of phase for H layer earth

均匀半空间以及层状模型的数值模拟结果与解析解很接近，可以说明矢量有限元方法及程序的正确性.

模型1：在地下均匀半空间(100 Ω·m)中存在一低阻异常体，电阻率为10 Ω·m，异常体埋深200 m，长200 m，宽200 m，高50 m，位于-75 m至125 m处.

模型2：在地下均匀半空间(100 Ω·m)中存在一高阻异常体，电阻率为500 Ω·m，异常体埋深200 m，长200 m，宽200 m，高50 m，位于-75 m至125 m处.

矢量有限元数值模拟的视电阻率和相位误差从高频到低频呈现“大—小—大”的规律，电阻率最小误差为0.00723%(50 Hz)，相位最小误差为0.008048%(50 Hz)，0.05~5000 Hz视电阻率和相位的误差均在3%以下，当频率小于0.01 Hz时，视电阻率和相位的误差随着频率的减小越来越大，视电阻率在0.001 Hz的误差甚至达到了116.7257%，而相位的最大误差为5.762017%.经过分析认为低频段误差大的原因是Z向剖分不够造成的，低频段电磁波探测深度大，分辨率降低，Z向剖分过小时大地电磁场在地下的传播并不满足底界面边界条件，这样插值计算出来的场值误差就会比较大.为了验证我的猜想，建立如下四个Z向不同剖分尺寸大小的均匀半空间模型(各个方向剖分块数不变，X和Y方向剖分尺寸不变)：

模型1：空气剖分尺寸为:8950 m，地下剖分尺寸为:14175 m.

模型2：空气剖分尺寸为:37000 m，地下剖分尺寸为:131850 m

模型3：空气剖分尺寸为:37000 m，地下剖分尺寸为:191750 m

模型4：空气剖分尺寸为:37000 m，地下剖分尺寸为:311750 m

下图为Z向不同剖分大小的四个模型的视电阻率曲线图(相位曲线图和此规律相同，略).

由图 11可以看出，随着Z向剖分尺寸的不断变大，低频段的视电阻率数值模拟结果与解析解的误差也越来越小，从而验证了我的猜想.因此在大地电磁的数值模拟中，要想保证低频时数值模拟结果准确，必须要保证Z向的剖分尺寸足够大.

图 11

Fig. 11

图 11 Z向不同剖分视电阻率对比图 Fig. 11 The comparison chart of resistivity for different partition in z direction

两个三层层状模型的视电阻率和相位的数值模拟结果与解析解差别很小，误差呈现和均匀半空间相类似的情形，但是在低频段(小于0.01 Hz)K模型的数值模拟相对误差比H模型的数值模拟相对误差要小，经过分析认为：对于K模型，基底的电阻率为10 Ω·m，电磁波在低电阻率模型中衰减快，因此也更容易达到底界面边界条件，从而计算出的结果也更加准确.

从图 7~图 10可以看出，对于低阻异常体，视电阻率拟断面图和相位拟断面图均在异常体位置(-75~125 m)出现了闭合的低阻异常区和闭合的低相位值异常区；对于高阻异常体，视电阻率拟断面图和相位拟断面图均在异常体位置(-75 m~125 m)出现了高电阻率和高相位异常.视电阻率和相位对高阻异常体的异常响应没有低阻异常体的异常响应明显，另外低阻相位的异常响应比低阻视电阻率的异常响应效果明显.

图 7

Fig. 7

图 7 低阻异常体视电阻率拟断面图 Fig. 7 Pseudosection map of resistivity for low-resistant anomalous body

图 8

Fig. 8

图 8 低阻异常体相位拟断面图 Fig. 8 Pseudosection map of phase for low-resistant anomalous body

图 9

Fig. 9

图 9 高阻异常体视电阻率拟断面图 Fig. 9 Pseudosection map of resistivity for high-resistant anomalous body

图 10

Fig. 10

图 10 高阻异常体相位拟断面图 Fig. 10 Pseudosection map of phase for high-resistant anomalous body

我们用伽辽金方法推导了含边界条件在内的加权余量方程，并且用矢量有限元的方法对三维的大地电磁场进行了求解，通过建立均匀半空间模型、三层层状模型和三维异常体模型进行正演数值模拟，并且对模拟结果进行了分析，证明了该方法以及程序的正确性.通过对以上数值模拟结果的计算分析得到如下结论：

(1)在利用矢量有限元方法进行三维大地电磁正演时，假设初始电场的偏振方向延X方向，在传统的三维大地电磁边界条件基础上，还需要给定四个垂直侧面的边界条件，前后面：E_z=0；左右面：E_y=0，E_z=0；底界面：E_y=0.

(2)剖分对大地电磁数值模拟的结果影响很大.各个方向剖分的块数越多，数值模拟的结果越精确.频率较大时，Z向剖分区域不大，数值模拟的结果也很精确；但是在低频时，Z向剖分区域不大，将会导致数值模拟的结果误差很大，这在实际计算模拟中要引起注意.要想保证最低频率f₀数值模拟结果的准确性，根据趋肤深度公式δ=503，要将Z向区域增大至1.5~2倍的趋肤深度δ.

(3)从三维异常体的大地电磁正演数值模拟结果中可以看出，低阻体的视电阻率和相位异常响应比高阻体的视电阻率和相位异常响应更加明显；对于低阻异常体，相位的闭合异常响应相对较好.