基于谱元法的广域电磁法三维正演模拟

2022-04-08徐锦通汤井田

地球物理学报 2022年4期

徐锦通，汤井田,2,3,4*

1 中南大学地球科学与信息物理学院，长沙 410083 2 中南大学有色金属成矿预测与地质环境监测教育部重点实验室，长沙 410083 3 有色资源与地质灾害探查湖南省重点实验室，长沙 410083 4 自然资源部覆盖区深部资源勘查工程技术创新中心，合肥 230001

0 引言

可控源音频大地电磁法(CSAMT)是一种重要的地球物理勘探方法.具有观测范围广、工作效率高、测量精度高、适应性强等特点，被广泛应用于矿产资源勘查、地下水及工程勘查等领域.但也存在一定的缺点，可控源音频大地电磁法沿用了大地电磁法的观测方式以及卡尼亚视电阻率计算公式，为满足远区测量要求，通常需要在距离发射源数千米至数十千米的区域内接收正交的水平电磁场信号来获得卡尼亚视电阻率和相位信息，背离了采用人工源使信号增强的初衷，并且在野外实际施工中，由于测区场地条件的限制，有时难以满足远区测量的要求.采用的卡尼亚视电阻率计算公式是在假设远区条件下，将电磁场精确表达式中的高次项近似舍去，然后基于阻抗定义推导得到的，引入了一定的人为误差(何继善，2010).

广域电磁法是一类人工源频率域电磁勘探方法(汤井田和何继善，1994，2005；何继善，2010)，广域电磁法继承了CSAMT 使用人工场源克服场源随机性的优点，采用较短的收发距进行测量，摒弃了CSAMT远区信号微弱的劣势，用适合于全域的不进行简化的公式进行视电阻率计算，扩展了观测适用的范围.同时，广域电磁法发射源一次发送的是含多个主频成分伪随机电流信号，而非变频发送，并且在测量时只需要观测电磁场的一个分量，而不是观测一对正交的电磁场分量，大大提高了观测速度、精度和野外工作效率.

广域电磁法解释目前主要采用的是一维和二维反演技术.然而在地形和地质条件复杂地区，地下目标体实际是典型的三维异常体.此时，若仍然采用一维和二维反演，很难得到可靠的的反演结果，必须发展三维电磁反演技术.三维反演的一个重要难点是如何实现快速高精度的三维正演，目前主流的三维电磁数值模拟方法主要包括积分方程法(Xiong et al.,1992; Kuvshinov et al.2002; Zhdanov et al.,2006)、有限体积法(Jahandari and Farquharson, 2015;Lu and Farquharson, 2020)、有限差分法(Haber et al.,2000; Weiss and Newman,2003；卢杰和李予国, 2019)、有限元法(Ren et al., 2013, 2014; Grayver,2015,秦策等, 2020).广域电磁法三维正演方面，李帝铨等(2013)采用积分方程法实现了E-Ex广域电磁法三维正演；彭荣华等(2018)实现了基于二次耦合势的有限体积法广域电磁法三维正演计算；周印明等(2021)实现了基于矢量位和标量位的广域电磁法三维有限元数值模拟.

谱元法是将有限元和谱方法结合发展而来的一种高精度数值模拟方法.同时兼具有限元的灵活性和谱方法的高精度.与传统高阶有限元方法相比，由于采用了高斯正交基函数，插值节点为高斯分布型节点，可避免高阶等距插值的Runge现象.Patera(1984)最早提出基于Gauss-Lobatto-Chebyshev 基函数的谱元法，将其应用流体力学数值模拟.Rønquist和Patera(1987)提出了另一种基函数Gauss-Lobatto-Legendre基函数的谱元法.基于GLL基函数的谱元法结合降阶的高斯积分可得到完全对角的质量矩阵，基于GLC多项式的谱元法可以解析计算单元矩阵，理论上可获得更高的精度.在地球物理领域里，谱元法最早被应用于地震勘探中(Seriani et al.,1997;Komatitsch and Tromp, 1999，2002a, 2002b).在电磁勘探领域，Lee等(2006)将谱元法应用于电阻抗层析成像的正演求解器开发，证明了谱元法的高精度及快速收敛特性.Lee等(2006)将基于GLL多项式的混合阶谱元法用于求解三维矢量电磁波动方程和三维瞬变电磁问题中，与高阶有限元相比计算时间大幅减小.Zhou等(2016)将谱元法与区域分解结合应用于大规模航空瞬变电磁正演计算，有效的节省了计算时间，克服了一定程度的较低频率下由于计算机精度有限导致的系统矩阵严重病态的低频崩溃问题.Zhou和Shi等(2017)提出一种混合谱元法并且将散度约束条件加入到谱元法中，来克服低频崩溃问题，在井下电磁勘探正演中取得了理想的效果.Huang 等(2017，2019)实现了基于GLL基函数谱元法频率域航空电磁法正演及各向异性正演.刘玲等(2018)实现了基于GLC基函数谱元法的海洋电磁正演模拟.陈汉波等(2019)实现了基于混合阶矢量基函数的海洋电磁法各项异性正演.Yin等(2019)等实现了基于GLC多项式的航空电磁正演模拟.

谱元法当前主要应用于航空、海洋电磁法，在陆地可控源电磁法中的应用尚未见到，本文在前人工作的基础上将谱元法引入陆地可控源电磁法的分支——广域电磁法的正演计算.目前针对广域电磁法三维正演研究主要是基于位场方程，需要通过差分间接得到电场，会带来精度上的损失，本文采用电场双旋度方程作为控制方程直接获得高精度的电场解，利用一维模型和三维验证算法的正确性、精度及适用性，设计典型三维地电模型研究广域电磁法响应特征及其探测能力.

1 基本理论

1.1 广域电磁法边值问题

从微分形式的麦克斯韦方程组出发，取时谐因子为e-iω t，可得电性源广域电磁法满足的频率域麦克斯韦方程组

(1)

(2)

(3)

(4)

式中：D、E为电位移矢量和电场强度；B、H为磁感应强度和磁场强度；σ为电导率；ω为角频率；μ为磁导率；ε为介电常数；Js为外部激发的电性源.

人工源电磁法正演需考虑源的加载及源奇异性的影响，为了避免场源的奇异性，采用二次场算法，将总场分解成背景场和二次场，背景场由解析解计算得到，二次场满足的公式(1)、(2)的形式为

(5)

(6)

其中

(7)

式中，σp为背景电导率；Ep为一次电场.

将(6)式代入(5)式，消去磁场可得二次电场满足的矢量赫姆霍兹方程：

(8)

其中，k2=iωμσ+ω2με，为波数.

为保证电场的唯一性，需要考虑边界条件的加载，本文在正演时对计算区域进行扩边处理，当扩边区域足够大时可以近似认为二次电场已衰减为0，因此本文在外边界处选取 Dirichlet 边界条件如下

Es|∂Ω=0,

(9)

(8)式和(9)式共同构成广域电磁法边值问题.

1.2 边值问题的谱元求解

采用谱元法求解上述边值问题，首先将求解区域划分成若干个六面体单元，在每个单元内二次电场可用关于x，y，z方向的高阶矢量插值基函数展开

(10)

利用伽辽金加权残差法处理电场矢量方程(8)，令单元内的加权残差积分为0，并应用第一矢量格林定理与Dirichlet边界条件可得

(11)

式中，Ne为单元总数，Nb为单元棱边总数.

(11)式可写成矩阵方程的形式：

(A-k2M)Es=B,

(12)

式中，A为刚度矩阵，M为质量矩阵，Es为待求二次场列向量，B为右端源向量.

谱元法中可采用基于Gauss-Lobatto-Legendre(GLL)多项式的基函数和Gauss-Lobatto-Chebyshev(GLC)多项式基函数，本文选取的是GLL基函数，GLL基函数是定义在标准参考区间[-1，1]上的，在三维参考单元中ξ,η,ζ∈[-1,1]×[-1,1]×[-1,1]中，N阶矢量基函数可表示为一维基函数的混合阶张量积(Lee et al.,2006)：

(13)

其中一维GLL基函数表达式为

(14)

式中，LN(ξ)、L′N(ξ)为N阶Legendre正交多项式及其导数.ξi为GLL插值节点，即方程(1-ξ2)L′N(ξ)=0的零点.

物理单元的插值基函数与参考单元插值基函数的映射关系如下

(15)

(16)

其中J为坐标映射雅克比矩阵

(17)

由此可将(11)式中的物理坐标系下单元矩阵转换为标准参考坐标系下的表达形式

(18)

(19)

进而通过在参考域中求解线性方程组(12)可以获得各棱边上的场值，并通过式(10)插值获得任意点的场值.方程组的求解本文采用直接解法开源求解器PARDISO进行求解.

1.3 E-Ex广域视电阻率

广域电磁法根据场源形式以及观测方式可以分为多种，考虑到野外实际情况，目前为止采用水平电流源发射信号测量电场的x单分量的E-Ex广域电磁法应用最为广泛，本文研究主要针对这类广域电磁法.

这里以电场水平分量Ex来阐述E-Ex广域电磁法和广域视电阻率的概念.均匀大地表面水平电流源的电场x分量的计算公式为

(20)

式中，I为电流;dL为电偶极源的长度；k为均匀半空间的波数；r为收发距，即观测点距偶极子中心的距离；φ为电偶极源方向和源的中点到接收点矢径之间的夹角.

根据式(21)可以定义广域意义上的视电阻率：

(21)

式中，

ΔVMN=Ex·MN,

(22)

(23)

式中，MN为接收电极距.不难看出，公式(21)为自变量为复数的指数复变函数，采用一般的代数方法无法简单的提取出其中含有的的电阻率信息，广域电磁法提出的方法是采用计算机迭代的方法进行提取.

2 数值试验

为验证本文算法的正确性和精度，设计了典型地电模型进行正演计算与对比分析.计算平台为Linux Fedora 21小型AMAX服务器，内存250 GB，主频为2.3 GHz.

2.1 一维层状模型精度验证

为验证本文算法的正确性及精度，设计如图1所示的三层模型.第一层电阻率为10 Ωm ，厚度为500 m，第二层电阻率为100 Ωm，厚度为200 m，第三层电阻率为500 Ωm.发射源为位于坐标原点的电偶极子，电流为1 A，长度为1 m，接收点位置坐标为(6000 m，0，0)，发射频率为1～8192 Hz对数等间隔的25个频点.整个计算区域为100 km×100 km×100 km.

图1 一维层状模型示意图

分别采用3阶、4阶基函数的谱元法进行正演计算，并与有限单元法及解析解对比.图2给出了不同阶数基函数的Ex分量幅值随频率变化曲线及与解析解的相对误差.从图中可以看出，有限元法、谱元法Ex响应均与解析解拟合较好，高频部分误差相对于低频要小，精度更高，这是因为高频电磁波衰减较快，更加符合扩边区域外边界二次电场衰减为零的边界条件假设. 对于谱元法，随着基函数阶数的提升，相对误差显著减小，表明谱元法可以通过提高基函数的阶数达到提高精度的目的.表1统计了3阶谱元法与有限元两种方法的计算参数，从表中可以看出，当基函数的阶数为3阶时，两者的精度相当，但谱元法的自由度相对于有限元更少，相应的计算时间也更少，说明高阶谱元法可以在稀疏网格下以更少的自由度获得更高精度的解，对于网格依赖性较低，从而可以节省计算资源消耗.

图2 一维层状模型Ex响应及相对误差曲线

表1 谱元法与有限元法计算时间统计表

2.2 三维模型对比测试

为进一步验证本文算法的精度及对三维模型的有效性，设计如图3所示的三维低阻体模型.背景电阻率为50 Ωm，低阻异常体电阻率为5 Ωm，顶部埋深为100 m，几何尺寸为120 m×200 m×400 m，发射源为沿x方向的电偶极子源，长度为100 m，中心坐标为(50 m，0 m，0 m)，距离异常体中心1 km，发射电流为1 A，发射频率为3 Hz.采用背景为50 Ωm均匀半空间模型进行一次场计算.在异常体正上方沿x方向布置一条测线，测线范围为400～1400 m，点距为20 m.整个模拟区域为80 km×80 km×80 km,最小网格尺寸为50 m×50 m×50 m，计算时间为153 s.

图3 三维低阻体模型示意图

图4为本文算法结果与积分方程法(汤井田等，2018)，有限元法(Ansari and Farquharson，2014)等计算结果Ex分量实部与虚部对比曲线.从图中结果可以看出，本文的结果与前人积分方程、有限元法结果均吻合较好，验证了本文算法对三维模型良好的适用性及精度，可以用来进行广域电磁法三维正演分析.

图4 低阻体模型Ex响应曲线

2.3 广域电磁法响应特征分析

2.3.1 低阻异常体

为分析广域电磁法响应特征及探测能力，设计如图5所示的三维低阻异常体模型，电阻率为1 Ωm的低阻体赋存在电阻率为100 Ωm的背景介质中，异常体大小为500 m×500 m×500 m，顶部埋深为300 m.发射采用沿x方向的电偶极子源，长度为1 m，发射电流为1 A，频率为0.5～8192 Hz对数等间隔的15个频点.源的坐标为(0 m，-5000 m，0 m)，取异常体中心在地面投影为坐标原点.在异常体正上方沿x坐标轴布置一条测线，长度为2 km，点距为50 m.整个模拟区域为80 km×80 km×80 km,最小网格尺寸为100 m×100 m×50 m，计算时间为127 s.

图5 三维低阻体模型示意图

利用本文算法对该模型进行正演，并同时计算卡尼亚视电阻率和广域视电阻率.图6为该模型卡尼亚视电阻率与广域视电阻率频率切片图.从图中可以看出无论是卡尼亚视电阻率还是广域视电阻率均能较好的反映出低阻异常体的存在，表现出明显的低阻响应.但在频率低于10 Hz时，卡尼亚视电阻率整体都出现了明显的升高，大多数观测值都高于背景电阻率，这是因为随着频率的降低，电磁波在地下介质传播的趋肤深度变大，此时测线所在位置已经不满足远区(五倍以上趋肤深度)的条件，出现了由场源引起的非平面波效应，卡尼亚视电阻率公式不再适用，产生了畸变.而广域视电阻率在低频段仍能较好的反映出异常体存在及地下电阻率信息，不受远区测量的限制，表明广域电磁法相对于传统可控源电磁法具有更大的观测范围.

图6 低阻体模型卡尼亚视电阻率与广域视电阻率切片图

在异常体上方布置一2 km×2 km的测网，线距为100 m，点距为100 m.发射频率为8 Hz，其他源参数保持不变.图7为卡尼亚视电阻率与广域视电阻率测区平面图.源的位置为(0 m，-5000 m,0 m),图中从上到下观测位置离源的距离依次减小.图中可以看出，卡尼亚视电阻率与广域视电阻率均对低阻异常体有响应，随着测线与源的距离的减小，卡尼亚视电阻率出现了急剧的增大，无法真实客观的反映出地下电阻率信息.而广域视电阻率在整个测区内均能较好的反应出异常体的位置及地下电阻率信息.表明了广域电磁法在较小收发距下也能获得地下真实地电信息，即在同等收发距下相对于CSAMT 具有更深的探测能力.

图7 低阻体模型卡尼亚视电阻率与广域视电阻率平面分布图(8 Hz)

2.3.2 组合异常体

为进一步探究广域电磁法对复杂异常体的分辨能力，设计如图8 所示的组合异常体模型.均匀半空间背景电阻率为100 Ωm，低阻异常体电阻率为1 Ωm，高阻异常体电阻率为1000 Ωm，顶部埋深均为200 m，大小均为400 m×400 m×400 m.在异常体上方布置一2 km×2 km的测网，测线沿x方向布置，线距100 m，点距50 m.发射源的位置为(0 m，-5000 m，0 m)，采用x方向电偶源，长度为1 m，发射电流为1 A，发射频率为8 Hz、32 Hz.整个模拟区域为80 km×80 km×80 km,最小网格尺寸为100 m×100 m×50 m，计算时间为162 s.

图8 组合异常体模型示意图

图9为电场Ex分量和磁场Hy分量平面等值线图.由图可知，电场、磁场均表现出随收发距增大而衰减的规律.电场较好的反映出了地下异常体的存在，低阻体上方呈现低电场值的特征，等值线向源一侧弯曲，高阻体上方呈现高电场值特征，等值线向背离源一侧弯曲.而磁场未出现明显的高低阻异常特征.图10分别为8 Hz、32 Hz卡尼亚视电阻率与广域视电阻率平面分布图.从图中可以看出，在32 Hz时，所有测点基本满足远区条件，卡尼亚视电阻率与广域视电阻率表现出类似的特征，未出现明显的畸变.频率为8 Hz时，部分测点已不满足远区条件，卡尼亚视电阻率受非平面波效应的影响，在靠近源一侧的部分测点视电阻率值出现明显的增大，与地下真实电阻率相差较大，并且对高阻异常体响应有一定的覆盖.而广域视电阻率仍能较好的反映出两个异常体的几何边界位置和地下电阻率信息，体现了广域电磁法对复杂地下异常体良好的分辨能力.此外，我们可以看出，无论是卡尼亚视电阻率还是广域视电阻率对于低阻体的响应相对于高阻体都要更强.这是因为低阻异常体吸引电流，使得异常体上方电场等值线变密，而高阻体排斥电流，使得异常体上方电场等值线变疏，低阻异常体电场响应更强.而磁场对高阻体，低阻体均不敏感，因此无论是基于电场和磁场的比值提取到的卡尼亚视电阻率还是仅用电场值提取到的广域视电阻率对于低阻体都更加敏感.