基于谱元法的低频平面波电磁场数值模拟

万文武, 李长伟, 熊 彬, 高 磊, 宋卫文

(桂林理工大学 a.地球科学学院; b.广西隐伏金属矿产勘查重点实验室, 广西 桂林 541006)

0 引 言

在过去的几十年里, 地球物理电磁法取得了快速发展[1], 出现了大量的研究成果。传统电磁正演主要包含了积分方程法(IEM)、 有限差分法(FDM)和有限元法(FEM)[2-6]。与FDM和FEM相比较, IEM的效率非常高, 主要是其离散的区域和方程组都非常小, IEM解的精度与对称矩阵的稠密性息息相关, 难以达到高精度。FDM使用交错有限差分法离散控制方程, 应用于电磁场中, 计算简单、 编程容易实现, 被广泛应用于电磁正演模拟中，其缺点在于无法使用非结构化网格, 难以模拟复杂地质体。FEM具有较好的灵活性, 能够使用结构化网格和非结构化网格。李长伟等[7]利用仿射坐标变换技术, 应用于平行多面体单元的单元分析, 极大地提高计算效率, 成为有效的电磁正演模拟方法, 特别是在复杂地形的正演研究中。熊彬等[8]借助COMSOL 生成完全非结构化三维网格, 引入Coulomb 规范和适当的边界条件保证矢量位唯一。在有限元方法中, 为提高求解精度需要对网格加密, 会导致计算量增加；另一种提高精度的策略是提高插值多项式的阶数, 即高阶有限元法。Schwarzbach等[9]和Grayver等[10]运用高阶有限元结合面向目标的局部网格细化, 可以达到较高的正演精度。

谱元法(spectral element method, SEM)最早是Patera在1984年提出的求解流体Navier-Stokes方程的一种数值方法[11]。 近年来, SEM在计算地球物理中得到广泛应用, 如Seriani等 、 Komatitsch等提出基于Chebyshev和Legendre多项式的谱元法模拟地震波[12-13], 能够模拟物性不连续介质中的体波和面波的传播特征。随着开源包SPECFEM的出现, 在地震研究中SEM成为解决局部和区域问题的一种流行方法[14-16]。如Magnoni等[17]采用了SEM对2009年意大利中部发生的6.3级地震产生的复杂波场进行模拟；王秀明等[18]引入基于预条件下共轭梯度法的元算法；刘有山等[19]和李琳等[20]基于Fekete和Cohen节点实现了谱元法采用三角网格剖分的地震模拟。在地球物理电磁法领域, 一些研究者将SEM应用于航空电磁和可控源电磁法的正演模拟， 如Lee等[21-22]提出了基于Gauss-Lobatto-Legendre(GLL)多项式矢量电磁波方程和三维瞬变电磁问题， Ren等[23]提出了一种新的时间域伽辽金SEM, 减少了未知数个数； Liu等[24]利用高阶混合阶SEM分析了各向异性和波导, 在边界上使用了旋度一致的GLL多项式和节点基于GLL的标量基函数来计算正向和切向分量得到场值； Huang等[25]利用基于GLL插值多项式基函数的谱元法实现了频域航空正演模拟； Zhou等[26]将基于GLL多项式的谱元法应用到井地电磁中解决了低频电磁场问题； 刘玲等[27]实现了基于GLC插值基函数的谱元法在频域三维海洋可控源电磁正演模拟研究。

与有限元相似, SEM也是将区域剖分成互不重叠的子单元, 在每个子单元内基于谱方法进行插值。 SEM能够用结构化或非结构化网格模拟复杂的几何形状, 但与传统的有限元相比, 谱元法具有指数收敛特性。SEM结合了有限元的灵活性和谱方法的高精度, 减少了内存和CPU的需求。为提高大地电磁法正演计算的效率及精度, 本文将高精度的谱元法应用到低频平面波电磁场数值计算中。 在研究区域使用稀疏网格, 采用GLL多项式构建插值基函数, 利用GLL多项式的完全正交性, 形成对称的质量矩阵, 减少了计算时间, 并通过高阶插值基函数来提高计算精度。与国际标准模型COMMEMI2D-0、 COMMEMI2D-1计算比对, 验证了本文算法的可行性和有效性, 还设计了几个含有构造特征的典型地电模型, 分析其大地电磁响应特征, 为实际野外地质工作提供参考。

1 谱元法基础理论

1.1 控制方程

采用时变因子e-iwt, Maxwell方程可写为

(1)

其中:E为电场强度(V/m)；H为磁场强度(A/m)；B是磁感应强度矢量(T)；ω为角频率;σ为介质的电导率;ε为介电常数;μ为介质的磁导率;ε、μ分别取真空中的介电常数和介质磁导率, 即ε=ε0=1/36π×10-9(F/m)、μ=μ0=4π×10-7(H/m), 对Maxwell方程组(1)中的第一个方程取旋度, 得到E满足的矢量Helmholtz方程

(2)

1.2 SEM基函数

采用高斯-洛巴托-勒让德(Gauss-Lobatto-Legendre)多项式作为插值基函数, 一维参考域的N阶GLL插值多项式[21]的表达式

(3)

其中:LN(ξj)为Legendre多项式;ξj为GLL插值节点;j=0, 1,…,N。

对于二维等参单元, 其N阶GLL插值基函数为

(4)

(5)

1.3 GLL插值基函数节点和阶数的选择

在SEM中, 选择基函数和用于物理量插值的配置点是决定SEM法能否成功地关键。在标准单元中, 选择不同函数展开形式对计算效率和数值精度都有很大影响。由于Legendre正交多项式的SEM, 在单元数值积分节点和函数插值节点选择的是同一点, 在数值计算过程中产生的质量矩阵是对角阵, 具有简化数值计算和提高计算效率等优点, 从而本文选择GLL积分。根据刘玲等[27]和Yin等[28]的研究, 当插值阶数小于4阶时, 必须通过增加单元个数才能保证SEM的高精度求解, 避免数值计算误差较大问题。理论上, 当多项式插值阶数无限增大时, SEM计算的结果越逼近精确解, 但计算机内存和计算时间相应增加。通常采用4～8阶多项式插值基函数, 考虑到精度和计算量的平衡, 本文选择了4阶的GLL插值多项式。

提高SEM的精度可以通过网格单元数量和增加每个参考单元的插值节点或提高多项式基函数的阶数, 也就是说, 网格大小和插值基函数的阶数直接决定了SEM的精度和耗时。在选取网格大小和基函数的阶数时, 网格尺寸越小和基函数阶数越高, 随之精度越高, 越接近数值解,但需要计算机的内存越大和计算时间越长。网格剖分在SEM计算中是一个极其关键的步骤, 需要不断重复调试适合的网格单元(包括网格形状、 网格数量、 网格尺寸)。

1.4 变分方程离散化

将计算区域Ω离散成一组互不重叠的子区域Ωe。对于任意四面体单元Ωe, 利用一个等参变换将其转换成标准参考单元, 即设物理单元坐标和标准参考单元坐标之间的映射函数为Fe:Λ↔Ω, 对刚度矩阵和质量矩阵使用同样的映射函数转化到参考域中计算。物理单元与参考单元之间的转换如图1所示, 右边为物理坐标(x,y)下的四边形单元, 左边是标准参考(ξ,η)下的四边形单元,ξ,η∈[-1, 1]。

使用协变量映射[29]表达式计算参考域下基函数

(7)

(8)

(9)

将微分方程(2)转化为积分问题, 采用伽辽金

图1 物理单元与标准参考单元的转换

加权残差法, 令加权残差为0, 即

∬Ωevi·(×μ-1(×E)-iω(σ-iωε)E)ds=0,

(10)

其中,vi为加权函数。

基于格林第一定理和狄利克雷边界条件n×E|∂L=0, 把式(10)改写成

∬Ωe(μ-1×Vi)T·(×E)ds-

(11)

(12)

其中,Na是基函数的个数。对单元进行集成, 得到最终的线性方程组

(A-B)E=0,

(13)

其中:A为刚度矩阵;B为质量矩阵。

(14)

(15)

其中,k是单元总数。对于磁场H计算与上述类似的方法对Maxwell方程组中第二个方程取旋度得到, 使用了稳定双共轭梯度算法求解上述线性方程组, 可得到电磁场的分布。

2 算例分析

为了验证算法的正确性和有效性, 本文设计了从简单到复杂的地电模型进行测试, 采用4阶GLL多项式作为插值基函数, 在如下的微型计算机平台上进行运算:处理器 Intel(R) Core(TM)i5-4210 CPU@1.70 GHz 2.40 GHz, 内存8 GB, Windows 10 64位操作系统。

2.1 算法验证

算例1采用国际标准模型库中的COMMEMI 2D-0和COMMEMI 2D-1模型[30]对本文算法进行验证, 模型如图2所示。

图2 模型COMMEMI 2D-0(a)及COMMEMI 2D-1(b)示意图

在图2a中，x∈[-50, -10]区域上电阻率为10 Ωm、x∈[-10, 10]区域上电阻率为1 Ωm、x∈[10, 50]上电阻率为2 Ωm。频率取1/300 Hz, 计算地面上x=-25、 -15、 -7、 7、 15、 30 km处的视电阻率, 本文算法与积分方程法的结果对比如图3所示, 两种方法都能够反映出垂直地层电阻率的变化趋势, 结果几乎一致, 说明了SEM算法用于2DMT正演的正确性。

图3 IE与SEM的视电阻率

图2b是在电阻率为100 Ωm的均匀半空间中嵌入一个电阻率为0.5 Ωm的异常体, 异常体在x方向上的位置是-0.5～0.5 km, 长为2 km, 宽为1 km, 上顶面离地面 0.25 km。 计算地面上x=0.5、 1、 2、 4、 8、 16 km处的视电阻率。频率取0.1和10 Hz时的SEM与IE的计算结果对比如图4所示, 两种方法在不同频域时的模拟结果几乎一致, 误差小于5%, 能够反映出地下低阻异常体的特征, 说明了SEM算法用于2DMT正演计算的准确性。

2.2 二维模型谱元法效率分析

分别采用本文算法(SEM)和双二次插值的有限元算法(FEM)计算模型COMMEMI 2D-0, 对二者的计算结果进行了对比。实验中设计了由粗到细4套不同的网格, 在粗网格1上应用本文算法的4阶谱元法, 在逐渐细化的网格上应用FEM进行计算。

图4 IE与SEM在频率0.1(a)和10 Hz(b)的视电阻率

表1给出了SEM和FEM的自由度和计算时间比较, 图5是FEM在不同尺寸的网格(Grid2、Grid4)和SEM、 IEM的计算结果比对。结果表明, 相同网格尺寸下SEM比FEM(Grid2)精度高、 耗时长, 但在同等精度下(FEM(Grid4))SEM所用计算时间更少。因此与FEM相比, SEM在粗网格下能够保证精度, 减少计算时间。

表1 谱元法与有限元法的计算时间对比

图5 IE与SEM(Grid1)、 FEM(Grid2/Grid4)的视电阻率

2.3 水平地层中的异常体模拟

图6 层状地层低阻异常体(a)及高阻异常体(b)示意图

图7 低阻异常体(a)及高阻异常体(b)视电阻率等值线图

2.4 复杂地形模拟

2.4.1 地垒模拟 为研究起伏地形条件下平面波电磁响应, 算例3设计了一个地垒模型, 如图8所示。在地垒模型中, 隆起部位上端宽为10 km, 下端宽为20 km, 凸起垂直高度为5 km, 地下均匀半空间的电阻率为100 Ωm, 空气层电阻率取1010Ωm, 频率为f=0.1、 1和100 Hz, 计算地面上的视电阻率和相位, 结果如图9所示。起伏地形对视电阻率和相位均发生了畸变, 视电阻率随着地形的凸起而增大, 且频率越大响应越明显, 与视电阻率相比, 相位受地形的影响更为明显。

图8 地垒模型示意图

2.4.2 覆盖层下的垂直断层模拟 算例4设计了覆盖层下有不同岩性的垂直断层带, 如图10所示。上层覆盖层电阻率为50 Ωm； 第二层是直立的断层带, 电阻率依次为30、 1、 100、 30 Ωm； 第三层为沉积层，其电阻率为100 Ωm； 第四层为基岩, 电阻率为0.1 Ωm。在频率0.01～10 Hz范围内进行正演计算, 得到如图11所示的视电阻率分布。

可看出SEM能够较好反映覆盖层下的断层, 在频率为1 Hz的时候反映了覆盖层的电阻率情况, 探测深度随着频率的减小而增大; 频率在0.15 Hz时的视电阻率响应接近垂直断层, 在x方向上-60～-20 km和30～60 km显示了30 Ωm的异常, 在-20～10 km能够准确得到低阻异常(1 Ωm), 在10～30 km清晰反映了高阻异常(100 Ωm); 再次减小频率, 得到沉积层和基岩电阻率, 但在-15～10 km处受到低阻断层异常的影响, 视电阻较低, 说明在层状地层深部受到低阻断层的影响较大, 而在覆盖层中受到高阻断层的影响更加明显, 说明在层状地层中浅部受到高阻异常的影响较大。

图9 地垒视电阻率(a)及相位(b)

图10 覆盖层下的垂直断层示意图

2.4.3 简单断层与异常体模拟 算例5给出了一个更为复杂的地电模型(图12), 其背景电阻率为100 Ωm, 存在电阻率为10 Ωm的4个低阻体, 分别位于x∈[0,14]、z∈[8,12];x∈[8,12]、z∈[26,30];x∈[20, 40]、z∈[8, 12], 在x为20～30 km是倾斜的;x∈[30,34]、z∈[22,30]。在其深度为2～3 km处有一条高阻断层, 电阻率为1 000 Ωm，厚度为0.5 km, 在x方向20 km处向下错动0.4 km。在频率0.01～1 Hz范围内计算得到频率电阻率等值线图(图13)及不同频率的相位曲线(图14)。

图11 频率视电阻分布

图12 断层岩体下的异常体示意图

在图13中电阻率等值线图能够精确反映地电结构, 清晰地显示出高阻断层和低阻体的存在, 右边低阻倾斜体没有完全吻合的原因是由于使用阶梯网格和网格尺寸较大所致。

图14显示了不同频率的相位变化曲线, 相位曲线更加清晰的反映了各个深度的地电结构。在地表1 Hz的相位清晰地反映了高阻断层岩, 0.25 Hz的相位能够明显反映浅部两边低阻体和中间的背景场在0.063 Hz体现了上下低阻体中间位置的电场变化, 在低阻体所在位置均有变化趋势, 在右端低阻体相位变化趋势比左端低阻体相位变化平缓, 是由于受到右端倾斜低阻体的影响。 在0.037 Hz相位(图14b)能够清晰地反映底部低阻体和背景场, 能够体现两个低阻体在x方向的准确位置; 在0.024 Hz仍能够反映低阻体的存在; 在0.01 Hz相位趋于平缓, 与背景场趋于一致。

图13 视电阻率等值线图

2.4.4 深部地电结构简化模拟 为研究深层流体穿过高阻岩石圈形成传导通道时的电磁响应特征[31-33], 设计了算例6，如图15所示。 第一层为水平层低阻体, 层厚1 km, 电阻率为30 Ωm； 第二层为高阻体, 电阻率为2 000 Ωm, 层厚20 km, 含有2个垂直狭窄的通道ρ1(其高度为7 km), 下半部分是一个地垒结构的低阻体, 电阻率为10 Ωm, 模型的类型取决于ρ1的选择。观测点取传导通道中心点O和离中心点25 km的点R, 在频率0.01～10 Hz范围内进行计算，得到存在流体通道(ρ1区域)和不存在流体通道在不同观测点处的视电阻率，如图16所示。

图14 不同频率(1 、 0.25、 0.063、 0.037、 0.024、 0.01 Hz)相位曲线

图15 深部地电结构简化模型示意图

图中实线表示ρ1=2 000 Ωm, 即不存在垂直的流体通道情况下, 在中心点O和远点R处的视电阻曲线, 低频时无论观测点R还是O都准确地反映了深部地电结构, 在高频时, 远点R的视电阻率高于中心点O, 这是因为中心点O处的视电阻率受到底部低阻地垒的影响; 虚线表示ρ1=10 Ωm, 即存在一个双向低阻溶液流体通道, 在中心点O和远点R处的视电阻率曲线, 视电阻率曲线和无导电通道的形态相似, 但在中心点O处, 视电阻率低于无流体通道时的视电阻率值, 在远点R处视电阻率与无流体通道时的值一致。 由分析可知, 大地电磁响应对低阻流体通道有相当高的灵敏度, 这些通道将近地表和深部导体连接起来, 形成闭合的电流电路, 导致地壳和上地幔中的大地电流重新分布。

图16 视电阻率对比

3 结 论

本文实现了基于谱元法的低频平面波电磁场的数值模拟。在数值实验中, 对国际标准模型COMMEMI2D-0、 COMMEMI2D-1进行计算, 验证了算法的有效性和准确性, 与有限元方法对比, 表明了本文算法具有更高的计算精度和效率, 通过设计复杂地电模型, 表明了谱元法应用于大地电磁正演的有效性, 并进一步分析了这些典型模型的大地电磁响应特征, 对实际工作提供了借鉴和指导意义。为进一步提高算法的实用性和计算效率, 将在以后的研究中考虑将有限元和谱元法相结合。