基于非结构化网格的大地电磁2D自适应有限元正演模拟

郭家松， 秦 策， 乃国茹, 谢卓良， 冯 凯

(成都理工大学 地球物理学院，成都 610059; 2.河南理工大学 物理与电子信息学院，焦作 454000)

0 引言

20世纪50年代，苏联学者吉洪诺夫(Tikhonov)[1]以及法国学者卡尼尔(L.Cagniard)[2]提出了大地电磁测深法的理论基础。该方法以天然交变电磁场为场源，基于电磁感应原理，研究地下介质在地表的电磁场分布规律，以此解决相关的地质问题。经过几十年的发展，大地电磁测深法因其成本低廉，对低阻反应敏感、勘探深度大等优点被广泛用于深部构造研究，油气勘探、地下水勘探与环境保护等方面[3-5]。

大地电磁测深法的最终目的是得到与实际地下介质分布情况最接近的地电模型。这个过程需对研究区的实测资料进行一系列正反演计算，使模型响应与实际资料达到最佳拟合。正演是反演的基础，大地电磁测深法的正演建立在Maxwell方程组之上。目前用于大地电磁测深法正演数值模拟的方法主要是有限差分法、有限单元法和积分方程法[6-7]。

这里着重研究有限单元法求解大地电磁正演问题。有限单元法求解变分问题的核心思想是将计算域划分为有限个互不重叠的单元从而在单元上进行插值求解[8]，求解域的剖分结果对计算结果有很大影响。在大地电磁测深法的正演问题中，采用规则化网格进行研究区域解的剖分局限性很大，剖分结构对复杂构造的地电模型拟合程度差，易产生较大的几何离散误差。笔者采用非结构化的四边形网格进行剖分，从而使剖分结果与实际模型能达到最大程度的拟合。

大地电磁具有趋肤效应，高频电磁波主要受浅部介质影响，低频主要受深部介质影响[9]。因而大地电磁正演问题过程中不同频率对网格的稀疏程度的要求不同，利用同一套网格求解不同频率下的大地电磁响应必然有误差。因此，在非结构化四边形网格剖分基础上引入自适应算法，从求解域的粗网格出发，在不同频点下利用网格单元的后验误差估计指导网格的自动加密，优化网格，以此提高大地电磁正演计算精度。

1 正演算法

1.1 边值问题

Maxwell方程组描述了电场、磁场以及电流密度、电荷密度之间的关系，是电磁理论的核心。大地电磁测深正演理论基于Maxwell方程组，从方程组出发，有谐频率为ω的定态电磁场方程为：

▽×E=iωμH

▽×H=(σ-iωε)E

(1)

式中:μ和σ分别表示介质的磁导率电导率；ε为介电常数。

这里主要针对二维地电结构，即介质的电性参数在有明显走向的倾斜岩层、背向斜等地质构造中只沿倾向和深度方向变化。因此，取x轴平行构造走向，y轴垂直于x轴，z轴垂直水平面建立笛卡尔坐标系，将电场和磁场沿坐标轴进行分解。根据上述定态电磁场方程，可得二维介质电场分量与磁场分量所满足的偏微分方程：

▽·(τ▽u)+λu=0

(2)

其中，对于TE极化模式，u=Ex，τ=1/iωμ，λ=σ-iωμ；对于TM极化模式，u=Hx，τ=1/iωμ，λ=σ-iωμ。

与一般求解偏微分方程一样，上述偏微分方程需要在特定的边界条件下求解。考虑如图1所示的研究区域，其中TE极化模式下的研究区域为ABCD，TM极化模式下的研究区域为A′B′CD(2)式引入边界条件后可形成如下边值问题：

图1 研究区域示意图Fig.1 Schematic diagram of the study area

(3)

图2 非结构化四边形网格剖分示意图Fig.2 Schematic diagram of unstructured quadrilateral mesh

1.2 自适应有限单元法

有限元是将边值问题转化为等价的变分问题进而求解。变分问题如下：

(4)

求解上述变分问题，需要对研究区域进行剖分。本次研究中，借助三维有限元网格生成器-Gmsh将二维研究区域剖分为非结构化的四边形网格(图2)，节点和四边形单元遍及整个研究区域。在单元内采用二次插值，将连续函数的求解转变为节点上的离散函数值的求解，进而将求解域离散为所有单元的线性组合，并利用网格单元内的二次插值函数可使变分表达式转变为向量与矩阵相乘的形式为式(2)。

(5)

其中:ue是网格单元节点场值组成的向量；K1e、K2e以及K3e是单元节点处的形函数组成矩阵；u是研究域所有节点场值形成的向量；Ke是将所有单元上求得的K1e、K2e以及K3e扩展为由全体节点所组成的矩阵。

对式(5)求变分并令其等于零，并加入变分问题的边界条件，可得如下线性方程组

Κu=P

(6)

求解上述线性方程采用LU分解法求得上述线性方程组，从而得到各节点的电磁场值，利用数值算法求得场值沿垂向的偏导数∂u/∂z得到辅助场，将辅助场代入式(7)、式(8)求得视电阻率和阻抗相位响应。

TE模式：

(7)

TM模式：

(8)

为满足网格在不同频率下的计算精度，在有限单元法求解大地电磁正演问题的过程中引入自适应原则。在每个频点下通过单元网格的后验误差估计值去指导网格的加密。本文的误差计算方法基于Kelly等[10]提出的有限单元后验误差估计方法，通过计算单元边界上场值梯度跳变值沿着单元边界的积分来近似作为每个单元的误差得到如下误差计算式：

(9)

式中:cF=hf/2，hf表示是单元K中最长的对角线；〔·〕表示单元K边界上场值梯度的跳变值。

2 模型试算

2.1 层状模型

为检验本文方法的有效性，构建有三层层状模型，模型的地电结构为K型，参数如下：ρ1=10 Ω·m，h1=1 000 m，ρ2=100 Ω·m，h2=1 200 m，ρ3=10 Ω·m。采用本文算法对该模型进行模拟计算，得到数值解与一维解析解的结果如图3所示。

表1 解析解计算误差Tab.1 Calculation error of analytical solution

从表1计算的解析解误差值中可以看出，网格剖分的稀疏程度影响着大地电磁有限元正演模拟的精度，细化次数越多，误差越小，有限元正演模拟的精度越高。图3为TE极化模式下求得的数值解和一维解析解的对比，结果显示网格细化程度越高，解析解和数值解逐渐吻合。但网格剖分较粗糙时，高频段计算出的大地电磁响应与解析解有很大误差，说明粗网格不能满足高频时的计算精度要求。大地电磁波具有趋肤效应，对研究域剖分结果的纵向网格的尺寸与电磁波的最小趋肤深度有关，因此在高频时需要采用较小网格进行计算[11-12]。图4即为100 Hz时的网格细化结果。结果显示在该频点下，网格细化主要集中在浅部，随着细化次数的增加，靠近地表的网格越密，尺寸越小，计算精确度有了明显提高。通过以上分析，利用本文的自适应有限元算法对网格进行局部加密，可以满足不同频率下的网格的疏密程度，且在对初始网格细化两次，数值解和解析解基本吻合，且误差较低，证实了本文算法的有效性。

图3 解析解和数值解的对比Fig.3 Comparison between analytic solution and numerical solution(a)视电阻率；(b)阻抗相位

图4 时不同细化次数下的网格细化结果Fig.4 Mesh refinement results for different refinement times at f=100 Hz(a)N=0；(b)N=1；(c)N=2

图5 断层破碎带模型示意图Fig.5 Schematic diagram of the fault fracture zone

图6 初始网格剖分Fig.6 Initial mesh generation

图7 f=100 Hz网格自适应结果Fig.7 Meshes adaptive result at 100 Hz

2.2 隐伏断层破碎带模型

构建图5所示逆断层破碎带模型，电性参数如图中所示。在水平方向3 000 m～4 000 m范围内有一电阻率为10 Ω.m的低阻破碎带，破碎带宽度为150 m，倾角为60°。

图6和图7分别显示的模型的初始网格剖分结果以及频率f=100 Hz时的网格自适应加密结果。图7显示在断层破碎带以及层位分界处网格被明显加密。

图8(a)、图8(b)分别是TE、TM两种极化模式下断层破碎带模型的视电阻率-频率拟断面图。从图8中可以看出，K型地电模型的三层电性层可以从结果图中清晰分辨，且电阻率与实际地层电阻基本一致，计算精度高。高频部分断层破碎带处出现等值线错断。上盘的表层厚度反应明显比下盘厚度薄，反应上盘地层有相对下盘地层而言有抬升。低频时，TE模式破碎带下方地层显示为凹陷形态。TM极化模式断面图在破碎带下方电阻率值较两侧低，形成纵向的带状低阻体，基底在破碎带下方被低阻体切断，形成断裂。

图8 视电阻率-频率拟断面图Fig.8 Apparent resistivity-frequency pseudosection(a)TE; (b)TM

图9 阻抗相位-频率拟断面图Fig.9 Impedance phase-frequency pseudosections(a)TE; (b)TM

图10 对称褶皱模型示意图Fig.10 Schematic diagram of symmetry fold model

图9(a)、图9(b)分别是TE、TM两种极化模式断层破碎带模型的阻抗相位-频率拟断面图。从图9中可以看出，高频时 可明显看出浅表地层底界面埋深有差别，左高右低。在频率为范围内，相位剖面显示为低异常，TM极化模式的等值线出现断裂。低频时，破碎带底部相位等值线呈现凹陷形态，TM模式为凸起，但反应不明显。

2.3 对称褶皱模型

建立图10所示的对称褶皱模型，模型划分有四层电性层，为HK型。电性层层厚及电阻率参数及褶皱的大小如图10所示，其中对称褶皱存在于水平方向3 000 m～6 000 m的范围内。以波长和波幅定义对称褶皱构造的大小，其中褶皱波长为2 000 m，波峰的幅度为500 m。

图11和图12分别显示的对称褶皱模型的初始网格剖分结果以及当频率f=100 Hz时，对初始网格细化两次之后的自适应加密结果。通过细化前后的网格剖分结果可以看出，细化两次之后在浅部，网格被明显细化加密。

图13(a)、图13(b)分别显示的是TE、TM两个极化模式对称褶皱模型视电阻率-频率拟断面图。 图13显示高频时， 可清晰地反应地层起伏。频率在100 Hz，视电阻率值等值线在褶皱处断裂，褶皱处电阻率值较两侧地层低。低频时TE模式剖面显示地层并未有明显的起伏变化，而TM模式剖面图清晰显示存在地层起伏。且随着频率降低，褶皱的波谷处凹陷形态反应越加明显，且波谷波幅随频率变低增大。

图11 模型初始网格剖分Fig.11 Initial mesh generation

图12 f=100 Hz网格自适应结果Fig.12 Meshes adaptive result at 100 Hz

图13 视电阻率-频率拟断面图Fig.13 Apparent resistivity-frequency pseudosection(a)TE; (b)TM

图14 阻抗相位-频率拟断面图Fig.14 Impedance phase-frequency pseudosections(a)TE; (b)TM

图14(a)、图14(b)分别显示的是阻抗相位-频率拟断面图。高频部分，地层有明显起伏，形态与模型褶皱形态吻合。但频率越低，起伏幅度越小。在频率范围内，同一频点相位值出现断裂。频率范围为1 Hz～10 Hz，TM模式拟断面图在褶皱波谷处显示为低阻圈闭异常。低频时阻抗相位在波谷处显示为地层的隆起。TE极化模式低频时并未表现起伏变化。

3 结论

采用了三维网格剖分器-Gmsh对研究模型进行非结构化四边形网格剖分，可以使剖分结果与实际模型达到最大程度上的拟合，能够真实的模拟地形起伏，断层构造，褶皱等复杂地质结构模型。采用自适应的有限元方法求解大地电磁正演问题，在每个频点下利用网格单元的后验误差估计指导网格的局部加密，满足了网格对不同频点计算精度的要求，进一步优化了网格质量和数量，避免了人为网格剖分所造成的影响，从而提高了二维正演计算精度。