电导率分块线性变化的二维高频大地电磁法有限元正演模拟

2012-01-11戴前伟张富强杨坤平李芸芸

物探化探计算技术 2012年5期

关键词：分块电导率电阻率

戴前伟，张富强，杨坤平，李芸芸

（中南大学地球科学与信息物理学院，湖南长沙 410083）

电导率分块线性变化的二维高频大地电磁法有限元正演模拟

戴前伟，张富强，杨坤平，李芸芸

（中南大学地球科学与信息物理学院，湖南长沙 410083）

为了适应电性参数在水平方向和垂直方向是连续变化的实际需要，这里采用矩形剖分，线性插值的有限元法研究了电导率分块线性变化的高频大地电磁高精度，快速正演模拟。首先给出了二维高频大地电磁测深的变分问题以及电导率分块线性变化时的有限元数值解法，编制了相应的有限元程序。利用该程序对层状模型进行了计算，并与解析解的结果做了对比分析，证明了程序的正确性和有效性，然后对典型电导率随深度线性变化的模型和典型地堑模型进行了正演模拟，结果表明能够有效地解决电导率在水平和垂直方向分块线性变化的高频大地电磁正演问题。

高频大地电磁法；电导率分块线性变化；正演模拟

0 前言

随着中深度地球物理勘探的迫切需要，以美国EH－4电导率成像系统为代表的高频大地电磁法，在我国地球物理勘探行业应用越来越广泛［1、2］，其勘探深度在地下1 000m 以内，采集的信号频率范围在10Hz～64KHz之间，开展高频大地电磁场的数值模拟有着重要的理论意义与实际意义［3］。国内、外众多学者，都对大地电磁法的正演模拟进行了研究：Wannamaker［4］等采用矩形对角线二次划分的剖分方式研究了二维带地形问题的大地电磁的有限元正演模拟；陈乐寿［5］等采用有限元法对二维大地电磁正演计算并提出了一些改进的方法；徐世浙［6、7］等采用矩形剖分的方式对电导率分块线性变化的大地电磁法进行了正演计算；陈小斌［8］研究了大地电磁正演计算中地形的影响；刘云［9］等采用自适应地形四边形网格剖分进行了大地电磁场的模拟；汤井田［10］等采用矩形剖分的有限单元法对高频大地电磁进行了有限元正演计算；阮百尧等［11、12］采用矩形剖分方式的有限单元法，对电导率分块线性变化地电断面电阻率测深进行了有限元数值模拟；欧东新［13］等对二维电导率线性变化测井进行了正演计算，取得了良好的效果。

大地电磁侧重研究区域地质构造划分，而高频大地电磁（EH－4）侧重于研究局部区域的电性变化，相比于大地电磁正演模拟的区域大小，高频大地电磁的计算区域相比要小，网格的大小要比较小的反应局部地质体电性的有微小变化。此外，高频大地电磁正演模拟所采用的频点是固定的，而大地电磁采用的频率是天然电磁场的随机频点，受野外电磁场干扰比较大。

前人在研究大地电磁或高频大地电磁有限元模拟时，大都假设电导率分块均匀，但实际工作中的岩矿石，由于组成、温度、湿度及孔隙度等的变化，其电性参数在水平方向和垂直方向是连续变化的。因此，研究采用矩形剖分，线性插值的有限元法进行高频大地电磁的数值模拟，有着重要的理论和实际意义。

1 变分问题

假定地电结构是二维的，选取右旋直角坐标系的原点在地面上，y轴正方向垂直向上，z轴方向平行于走向方向（见图1），即介质的电性参数仅是x、y两个坐标的函数。当平面电磁波以任意角度向地面入射时，地下介质中的电磁波总以平面波形式，几乎垂直地向下传播。这就使电磁场的各分量沿z方向都没有变化，即z方向的偏导数为零。根据Maxwell方程，并考虑到∂／∂z＝0，得到两个独立的方程组，并以z分量为准，分别称为TE模式和TM模式。

图1 大地电磁二维正演研究区域Fig.1 Studied area of 2Dforward modeling in magnetotelluric

（1）TE模式：

式（3）和式（4）可以统一表示为：

对于TE模式：

对于TM模式：

相应的边值问题归结为：

方程组（6）与下列变分问题等价［14］：

2 高频大地电磁法有限元求解

用有限单元法解变分问题（7）的步骤如下：

（1）网格剖分。首先将求解的二维区域剖分成矩形单元，如图2所示（在图2中，1、2、…代表节点；①、②、…代表单元）。

图2 网格剖分及节点编号示意图Fig.2 Sketch map of mesh division and node number

（2）线性插值。图3（a）是母单元示意图，图3（b）是子单元示意图，两个单元间的坐标变换关系为：

其中 x0、y0是子单元中点的坐标；a、b是子单元的两个边长，微分关系为：

图3 母单元、子单元示意图Fig.3 Sketch map of parent element and sub－element

双线性插值的形函数为：

在单元e上对电磁场u和参数τ、λ、k进行双线性插值。用ui、τi、λi（i＝1、2、3、4）表示各节点的u、τ、λ，单元中的u、τ、λ可表示为式（11）。

（3）单元分析。根据有限单元法的原理，对于区域的积分，可以划分为各个矩形单元的积分之和。

式（7）第一项：

其中 K1e＝（k1ij）；k1ij＝k1ji；ue＝（ui）T；i、j＝1、2、3、4

式（7）第二项：

其中 K2e＝（k2ij）；k2ij＝k2ji；ue＝1、2、3、4

式（7）第三项：

式（2）右侧最后一项线积分只对CD边界进行，但单元的23边落在CD边界上时，线积分为式（14）。

其中 K3e＝（k3ij）；k3ij＝k3ji；ue＝（ui）T；i、j＝1、2、3、4

（4）总体合成。将以上各单元积分相加，地单元（i，j）的Fe（u），将各单元的Fe（u）相加，便得到总体系数矩阵的F（u）为式（15）。

（5）视电阻率和阻抗相位的计算。令式（15）为零，得线性方程组，加入上边界条件u｜AB＝1，解线性代数方程组即可得各节点的u值，它代表各节点的（对TE模式）或（对TM模式）。在求得地表各节点的u后，再利用数值方法求出场值沿地表的偏导数∂u／∂y，代入式（16）、式（17），便可计算视电阻率和相位。

TE模式：

3 模型的计算

3.1 模型1－水平层状模型

水平均匀层状模型分四层：第一层电阻率为100Ω·m，厚度为 50m；第二层电阻率为1 000Ω·m，厚度为100m；第三层电阻率为500Ω·m，厚度为400m；基底电阻率为10Ω·m。正演所采用的频率为10Hz～10 000Hz。

用基于Fortran语言程序编写的电导率分块线性变化的高频大地电磁进行模拟计算，得出视电阻率和相位值并绘制出图4。如图4所示，在所研究的频点上，程序计算的结果与解析解都很接近，误差在3%以内。这说明程序在TE、TM二种模式下，计算均匀层状模型方面的正确性和有效性。在TM模式下，视电阻率值稍大于解析解，相位值稍小于解析解；在TE模式下，电阻率稍小于解析解，相位值稍大于解析解。纵观两种模式下的误差，高频段相对明显，高频段的误差比低频段的误差要稍大些，如果网格间距能够更小，更易于接近解析解。

图4 模型1曲线图Fig.4 Curves of model 1

3.2 模型2－电导率线性变化的层状模型

图5 为三层地电模型示意图，第一层的电阻率ρ1＝100Ω·m，厚度为h1＝30m；第二层的厚度h2＝150m；第三层的电阻率ρ3＝20Ω·m。其中第二层电阻率变化过渡层，其电阻率是深度线性关系为：

图6是用所编制的程序计算的TE及TM模式的视电阻率和相位拟断面图。从图6上可以看出：电阻率值随着频率的减小，从约100Ω·m逐渐线性减小到20Ω·m，与计算模型的电阻率变化相一致。此外，无论是TE模式还是TM模式，都能清楚地反映出中间层电阻率线性减小的变化和地堑构造形成的异常。尤其在TE模式的视电阻率和相位的拟断面图上，更能反映中间层电阻率线性变化的趋势。在TM模式下的视电阻率断面图上，中间层能明显地反应出来。但是在深度方向上，形态被明显向下拉伸，比较而言，TM模式的相位断面图对模型的反应不如TE模式的效果好。

图5 模型2Fig.5 The parameters of model 2

图6 模型2视电阻和相位拟断面图Fig.6 The apparent resistivity and phase pseudo－section map of model 2

3.3 模型3－地堑模型

图7 为地堑模型，区域大小为1 200m×500m。第一层的电阻率为10Ω·m，厚度为50m；在厚度为200m的中间存在一个地堑，长为150m，宽为150m；其中第二层电阻率变化过渡层，其电阻率是深度线性关系为：

具体参数详见图7。

图8是用所编制的程序计算的TE，TM模式的视电阻率和相位拟断面图。从图8上可以看出：电阻率值随着频率的减小从约10Ω·m逐渐线性减小到100Ω·m。在碰到地堑构造时，电阻率曲线明显发生变化，形成一个凹陷区域，与计算模型的电阻率变化相一致。无论是TE模式还是TM模式，都能清楚地反映出第二层电阻率线性减小的变化和地堑结构形成的异常。相比之下TE和TM模式的视电阻率拟断面图，都能够清楚地反映出地堑构造模型的电阻率实际分布情况。在TM模式下的相位断面图上，第二层能明显地反应出来，在地堑位置形成了向深部延伸的低阻异常。比较而言，TE模式的相位断面图对模型的反应不如TM模式的效果好。