附加列插值法在起伏地形条件下二维大地电磁正演中的研究与应用

吴 頔，辛慧翠，卜传新

(1.湖南省有色地质勘查局二四七队，湖南 长沙 410129；2.承德石油高等专科学校 计算机与信息工程系，河北 承德 067000)

起伏地形对大地电磁场影响很大，研究表明二维带地形条件下TM 和TE模式电阻率曲线具有不同的畸变特征，这种畸变既可能是电阻率数值发生畸变也可能是发生静态位移, 或者是二者的综合影响，畸变一般与相对波长的地形几何尺度有关，大地电磁测深地形影响的复杂性使得对它的有效校正较为困难, 本文利用双二次插值有限元法分别推导出了TM和TE模式在复杂地形条件下的二维正演公式，采用整体法和附加列法两种方法进行计算，并对计算结果进行了分析和评价。

1 变分问题

1.1 基本理论

根据麦克斯韦方程，假定二维地下电性结构，各坐标轴关系如图1所示，在均匀半空间里当平面电磁波以任何角度入射地面时，地下介质中的电磁波总以平面波形式近似垂直地向下传播，将(1)、(2)式按分量展开，并考虑到∂/∂z=0，得两个独立的方程组，以垂直分量为准，分别命名为TE和TH极化模式[1]。

TE极化模式：

(1)

TM极化模式：

(2)

式(1)、(2)可以统一写成下式：

(3)

其中：μ是介质的磁导率，σ是电导率，ε是介电常数，ω是角频率， e-iωt是时间因子，是二维哈密顿算符。

1.2 边值问题和变分问题

求解微分方程必须给出边界条件，建立以AB、CD为边界的模型(见图1和图2)，边值问题可以归纳为:

(4)

与边值问题等价的变分问题可以归纳为：

(5)

2 有限元法

采用不均匀矩形单元对区域进行剖分，每个矩形网格单元上取4个顶点和4条边的中点, 每个矩形网格的编号如下图所示[2-3]，图3为母单元, 图4为子单元, 两者的坐标关系为:

(6)

其中x0、y0为子单元的中点,a、b为子单元的两个边长。

将(5)式分解为单元积分之和，当单元中的τ和λ是常数时，各项单元积分可分别写为：

(7)

(8)

对(8)式求变分，由于δu具有任意性，因此可得：

Ku=0

(9)

加载式(4)中AB的边界条件，形成右端项进而可得到大型稀疏方程组，利用迭代法求解线性方程组。视电阻率计算公式如下：

TE型极化：

(10)

TM型极化：

(11)

在求得地表各节点的u后还要计算其对y方向的导数，取近地表的4个等距节点的u及节点1与节点4间的距离来计算导数，将计算结果代入公式(10)、(11)即可计算电阻率[4-6]。

3 程序实现及验证

3.1 程序实现

一般利用有限元模拟计算大地电磁2D正演分为以下几个步骤：1)网格剖分，单元格内电导率连续变化；2)单元内插值并进行单元分析，一般以整个模型为研究对象，对所有网格单元进行编号；3)总体合成；4)代入边界条件，解线性方程组；5)解出各节点u值，求导后计算视电阻率，计算过程中对模型整体进行插值因此这种方法也被称为整体法。在此基础上作者提出了一种改进的插值方法，以第N列(N=1…L，L为网格列数)网格为对象，增加一列与其完全一致的网格，以这两列为对象进行编号和计算，然后取这两列网格地表三个节点中的第二个节点及其下面的三个节点的U值来计算视电阻率值，所得结果即为整个模型中第N列网格表层中间点的视电阻率值，以此类推，从第一列循环计算到最后一列，从本质上讲该方法是对若干个层状模型进行组合，可有效改善一般算法中因地形起伏引起的分界面上的电阻率畸变，该方法称之为附加列法，计算思路见图5。

3.2 程序验证

为验证程序可靠性，设计三层一维模型，第一层厚度为300 m，电阻率值500 Ω·m，第二层厚度为500 m，电阻率值50 Ω·m，第三层厚度为900 m，电阻率值1 500 Ω·m，经正演计算出不同频率的电阻率计算(图6)。

一维模型电阻率只在深度方向发生变化，理论上无论TE模拟或者TM模式，电阻率计算结果与解析解应近似，由上图可以看出整体法和列分解法在TE和TM模式下都能较好的模拟电阻率变化，这也验证了方法的正确性。

4 二维模型计算

设计二维模型长度4 km，第一层中间厚度250 m，最大厚度为450 m，最小厚度为100 m，分别模拟山脊和沟谷；第二层厚度为300 m；第三层厚度为450 m(图7)。

图8～图11是应用整体法和列分解法计算的视电阻率拟断面图。其中图8是TE模式下整体法电阻率计算结果，该结果能较清晰地反映地形引起的电阻率变化且与地形构造大致吻合，但在地形变化分界面上存在由地下电性结构变化引起的干扰，表层高阻电性层未能较好体现。图9是TM模式整体法电阻率计算结果，正演结果在模拟山脊和沟谷处都出现了贯穿了整个频率带的低阻假异常，且两个低阻假异常间还出现了高阻假异常。图10和图11是TE和TM模式下附加列法电阻率计算结果，由图可知两种极化模式下的正演效果都得到了改善，其中TE模式下表层高阻得到了连续且清晰地反映，也很好的表现出了其随地形变化的特征，在TM模式下正演结果取得了很大的提高，消除了浅部(高频部分)高阻假异常以及模拟山脊和沟谷处所对应的低阻假异常，清楚地反映了模型的地形特征。

5 结语

本文利用有限元数值模拟技术在整体法基础上提出了更为精细的附加列插值法，对起伏地形进行了正演模拟，通过不同方法在TM和TE模式下计算结果的对比分析得出结论：1)整体法模拟中，TE模式正演结果比TM模式正演结果更加接近理论模型，TE模式正演结果基本反映出了地形模型的地电结构特征，模型表层存在些许干扰形成的高阻假异常，但影响较轻微。2)整体法TM模式下电阻率计算结果发生严重的畸变和位移，在地形发生变化的分界面上出现了严重干扰，假异常贯穿了整个频带且几乎掩盖了地形。3)附加列法模拟中，TE和TM模式下的正演结果基本相同，能够很好地反映表层地电特征，基本消除了由地形引起的干扰及假异常。4)通过对模型正演计算结果的分析和研究, 总结了起伏地形电阻率畸变规律，为正确识别和选取极化模式及大地电磁数据校正与处理提供了参考。