最优系数有限单元法色散介质GPR频率域正演模拟

王珣, 朱磊, 冯德山*, 许德如, 丁思元, 刘硕

1 中南大学 地球科学与信息物理学院, 长沙 4100832 东华理工大学 地球科学学院, 南昌 330013

0 引言

探地雷达(GPR)是一种利用高频电磁波探测地下内部结构及地质体物性变化的地球物理方法,它被广泛应用于城市道路病害、管线探测等浅地表勘探领域(Daniels,2004).土壤、混凝土、岩石等地下浅层介质大多具有色散特性,其介电常数是关于频率的复函数,导致电磁波的相速度和衰减系数不断变化,波形产生衰减与畸变(Wang and Oristaglio,2000;Oden et al.,2007).在GPR数值模拟中,简单地将地下介质设置为非色散介质,会影响GPR信号的精确解译.因此,需开展色散介质的GPR正演研究(Liu et al.,2007),以便更有效地指导雷达工程探测实践(刘四新和曾昭发,2007;朱尉强和黄清华,2016).

目前,GPR针对色散介质的数值模拟方法可分为时间域、频率域两大类.时间域数值模拟方法如时域有限差分法(Teixeira,2008;冯德山等,2014;Giannakis and Giannopoulos,2014;李静等,2010,2016)、时域有限单元法(冯德山等,2013;王洪华和戴前伟,2014;冯德山和王珣,2017;王洪华等,2018;Feng et al.,2018;Liu et al.,2019)、时域伪谱法(Liu and Fan,1999;李展辉等,2009;Huang et al.,2010;Fang and Lin,2012)、时域辛算法(方宏远和林皋,2013;Fang et al.,2013,2019;雷建伟等,2020;Lei et al.,2022)等被众多学者应用到Debye、Drude和Cole-Cole等色散介质的研究中.时域正演方法应用较为广泛,但对色散特性的处理会涉及卷积项,需要存储前一时刻的波场,且针对不同的色散模型需调整时间离散格式(Pratt and Worthington,1990).频率域数值模拟中电位移矢量可表示为电场强度与复相对介电常数的乘积,对色散项的处理更为简单直接,同时不受时间推进算法的限制,从而避免了求解所有时刻波场记录.因此GPR色散介质的频率域正演方法天然地具有原理清晰、易实现的优点(Bitri and Grandjean,1998).频率域正演方法主要有频域有限差分法(Shin and Fan,2012)与频域有限单元法(FEFD)(Jin,2015;Feng et al.,2019).频域有限差分法大都直接离散一阶麦克斯韦方程,求解变量多,计算内存需求较大.FEFD方法基于变分原理和加权余量法,求解变量较少、占用内存小,适应求解复杂电磁问题,因此本文采用FEFD进行色散介质的GPR正演.

利用FEFD方法对介质进行空间离散时,传播波场的采样点数量不足会导致GPR波动方程的数值解与真实解之间误差越来越大,产生数值频散(Ihlenburg and Babuška,1995;Guddati and Yue,2004).通常采用更细的网格尺寸或更高阶的插值基函数去压制数值频散,但这两种策略的计算成本较大,降低了正演效率.此外网格尺寸的细化伴随着一定的数值伪影,而高阶的方法在高频带会出现龙格现象(Marfurt,1984).为避免网格加密与高阶近似所带来的问题,优化低阶有限元的质量、刚度矩阵常被用于压制数值频散:传统方法是将权函数等同于形函数,该方法下就能得到一致质量矩阵(Mullen and Belytschko,1982);或是将单元的质量集中在单元节点上(王月英,2007),得到集中质量矩阵(Marfurt,1984).为了获取更精确的数值解,Marfurt(1984)以及Christon(1999)提出折衷质量矩阵,改善了积分单元的频散特性.此后,Min等(2003)在FEFD中针对不同剖分单元进行加权平均得到优化后的总体刚度、质量矩阵;Guddati和Yue(2004)通过改变双线性矩形单元积分点位置,构造新的质量矩阵与刚度矩阵.以上方法大都单独针对一致矩阵、集中矩阵或二者组合系数进行研究,或在牺牲计算资源的情况下对有限元的刚度矩阵、质量矩阵同时优化,难以获得较高的正演效率与精度.Ogunbo和Shin(2022a,2022b)受优化的九点格式有限差分法启发,采用最速下降法求解出五个参数组合下的优化系数,验证了优化矩阵在不同传播角度及宽频带下,比集中、一致和折衷矩阵具有更优越的压制数值频散的性能.

本文在Ogunbo和Shin(2022a,2022b)的研究基础上,结合有限元质量、刚度矩阵自身的约束条件,利用最小二乘方法得到仅需三个参数组合下的精确优化系数,所需的优化参数更少,避免了对源项进行归一化处理,具有更好的普适性.同时引入基于无限吸收函数的精确完全匹配层(EPML)(Feng et al.,2019)进一步提高了正演效率.通过频散曲线、算法收敛性、吸收边界以及复杂色散模型正演等对比实验,验证了本文算法相较于一致、集中和折衷矩阵等常规FEFD算法在压制数值频散方面更具优势,且能有效提高吸收边界性能.

1 基于EPML的色散介质GPR二维频率域控制方程

1.1 色散介质模型GPR频率域波动方程

根据Maxwell方程组,可由时间域标量波动方程推导出GPR二维TM波频率域传播方程.假设电磁波在z方向没有变化,引入时谐因子ejωt得到:

(1)

(2)

代入(1)式得到Helmholtz方程:

(3)

在GPR中常见的色散模型有Cole-Cole模型和Debye模型,许多土壤和岩石可以使用Cole-Cole函数去模拟(Giannakis et al.,2016).若不考虑电导率的影响,相对介电常数可表示为(Cole and Cole,1941,1942):

(4)

式中ε∞为介质的光频介电常数,ε0为静态介电常数,τ为偶极子的弛豫时间,参数β(0≤β≤1)表示弛豫时间分散的程度,参数β=0时,式(4)简化为Debye色散介质(Atteia and Hussein,2010).

1.2 EPML吸收边界条件

FEFD法应用于GPR正演时,需要对模拟区域进行人为截断.目前在截断边界处对反射波吸收效果较好的边界条件为完全匹配层(PML)(Berenger,1994;王洪华等,2019),结合GPR控制方程:

(5)

其中,sx与sy为坐标拉伸变量.常规PML取得最佳吸收效果的过程本质上是求解一个多参数最优化问题,不仅计算繁琐,且得到的参数不具有普适性(Gedney,1996).EPML吸收边界的提出,简化了边界条件的加载过程(Bermúdez et al.,2004;Cimpeanu et al.,2015).参考Feng等(2019)有关EPML边界条件的讨论,其sx与sy计算公式为

sx=1-iσx/k,sy=1-iσy/k,

(6)

电导率σx和σy表示EPML中不同方向的吸收函数,σx具体计算公式为

(7)

其中,LPML为PML层厚,dx表示PML内部单元中心到模拟区域边界的距离.式(7)所表达的吸收函数形式仅包含最基本的PML厚度及节点位置,可以极大简化参数优化过程,提高正演效率.

2 最优系数的频率域有限元法

2.1 有限元法

根据Galerkin加权余量法,并代入完美导体边界条件,得到式(5)的弱解形式:

=∬wjωμ0JzdΩ,

(8)

其中w为权函数,利用双线性插值基函数对场值Ez进行展开:Ez=∑NiEe,Ni为形函数.单元场向量Ee可表示为

(9)

(10)

其中Ke为单元刚度矩阵,Me为单元质量矩阵,fe为单元源向量.其具体表达式:

Me=∬Ωesxsyk2NiNjdΩ,

fe=∬Ωejωμ0JzNidΩ.

将单元场向量、源向量及单元刚度、质量矩阵扩展成整体矩阵,形成多源系统的线性方程组:

AE=f,

(11)

式中A是与频率相关的所有源共用的系数矩阵,f是离散化激励源项,E是对应节点的离散电场值.

2.2 最优参数选取

数值频散的程度取决于剖分单元大小、波的传播角度以及单位波长内的网格点数(梁昊,2016),选择合适的参数可以使数值频散最小化.根据Galerkin法,设四边形单元剖分下x与y方向上单元边长为lx与ly,Ke与Me具体表示为

(12)

(13)

对于传统的有限元方法:a=2,b=1,c=4,d=2,e=1.利用16个四边形单元对二维无源区域进行离散,剖分单元边长为h,如图1所示.

图1 计算区域剖分示意图

二维频率域电磁TM波,在xy平面内电场设其表达式为

E=E0e-jk(xcosθ+ysinθ),

(14)

其中,θ为传播方向与x轴的夹角.设波长为λ,单位波长内的网格点数为G,电磁波的传播速度为v,可知有:λ=2πv/ω,k=ω/v,G=λ/h.

剖分区域内部的节点周围都有八个节点,以图1中编号13的节点为例,节点的离散场值Ei与每个节点对应系数Ai满足关系(Chen et al.,2013):

A7E7+A12E12+A17E17+A8E8+A13E13+A18E18

+A9E9+A14E14+A19E19=0,

(15)

设节点13处的电场值为u0,式(15)可表示为

(16)

根据欧拉公式ejx=cosx+jsinx,式(16)可表示为

(17)

由波的传播规律可知:kh=2π/G.定义:

(18)

将(18)式代入式(17),由于u0不为零,化简得:

(19)

此处k是通过有限元方法得到的数值解波数,数值解相速度为Vph=ω/k.设真实相速度为V,则归一化相速度Vph/V可表示为

Vph/V的值越接近1,表示数值解与真实解的误差越小.Vph/V与a、b、c、d、e等参数相关,可通过选取最优化的参数使其接近1,设置目标函数:

J(a,b,c,d,e,G,θ)=

(21)

通常θ与G的取值范围(Shin and Sohn,1998)分别为:θ∈[0,π/2],G∈[Gmin,Gmax]=[4,400].

观察参数a、b、c、d、e在式(21)中分布特征,同比放大缩小一组参数会存在无数组优化参数,且改变相应的单元系数矩阵,需要对源向量进行归一化处理.因此需要对参数施加约束条件,根据质量守恒定理,令a、b、c、d、e满足:

a+b=3,c+2d+e=9.

(22)

根据上式约束条件,目标函数确定好G与θ后,实际上只有三个独立参数.若单独采用一个约束条件a+b=3或c+2d+e=9,则会有四个独立参数组合.不同参数组合下的目标函数形式不同.以参数组合ade为例,将式(22)代入式(21)中得到:

6G2(1-P-Q+PQ)a+2π2(2-P-Q)d

+2π2(1-PQ)e=18π2+9G2(2PQ-P-Q),

(23)

相应的矩阵形式为

(24)

其中l与r分别为G、θ的离散个数:m=1,2,…,l;n=1,2,…,r.令G、θ的离散值分别为Gn和θm,矩阵中元素S的表达式如下:

本文采用最小二乘算法求解超定方程式(24).结合约束条件(22),可获得优化系数a、b、c、d、e的值.为了方便对比,分别求出了bde组合、abde组合以及abcde组合下的优化系数,结果如表1所示.分析表1可知,ade与bde组合求出的优化系数相同,abde组合下的优化系数值大小与前两组组合基本相同.五个参数的优化参数abcde组合存在多解性,且需要对源项权重进行修改才能得到与真实解相同的场值.因此本文采用约束条件式(22)所得到三参数组合优化系数更加的准确,普适性更好.

表1 不同参数组合下得到的最优化参数取值

选取ade组合下的最优化参数,并得到常规、集中以及折衷类型(Ogunbo and Shin,2022a)下的a,b,c,d,e不同参数取值,如表2所示.

表2 不同类型参数取值

2.3 频散及精度分析

根据公式(20)可得到归一化相速度Vph/V与单位波长内网格点数G的频散关系如图2所示.随着G的减少与传播角度的改变,归一化相速度曲线逐渐偏离1,数值频散随之产生.一致、集中矩阵的曲线偏离程度较大,其次是折衷矩阵;而优化矩阵的归一化相速度曲线与1的偏离程度相较于前三种方法极小.归一化相速度与1的偏差在小于1%时,一致矩阵和集中矩阵的G至少大于12.9,折衷矩阵的G至少大于6.2.而优化系数矩阵的G仅需大于4.8的情况下,其归一化相速度与1的最大绝对误差即可小于0.2%.通过频散曲线分析,优化系数矩阵与一致、集中和折衷矩阵相比,在压制数值频散方面的性能极好,可显著降低归一化相速度的误差.

图2 (a) 一致、集中、折衷和优化矩阵的归一化相速度曲线; (b) 优化矩阵归一化相速度曲线

值得注意的是,当G为4.8时,优化矩阵的四条归一化相速度曲线会相交于一点.此点的意义十分重要:首先此点对应的归一化相速度与1的误差较小,说明压制频散的效果较强.其次,此点是不同传播角度下归一化相速度曲线相交的点,意味着满足此点的网格剖分条件下,即使在不同传播角度,优化矩阵都能达到相同的压制数值频散的效果.因此该点对应的G是最优化的单位波长下的网格点数.

3 探地雷达有限元色散介质模型正演模拟

3.1 网格精度对比

3.1.1 不同网格尺寸精度对比试验

模拟区域为10.0 m×10.0 m的Debye色散模型:ε∞=10.0,ε0=15.0,σ=0.001 S·m-1,τ为3 ns,将中心频率为100 MHz的雷达波脉冲激励源置于模拟区域中心.利用四种不同方法对模拟区域进行频率域正演,改变正方形剖分单元边长h,测试在不同网格尺寸下优化矩阵算法的收敛性.四种方法在满足相同剖分条件下的运行时间基本相同,h分别为0.04 m、0.05 m、0.08 m、0.1 m、0.2 m,五种剖分条件下的计算时间分别为0.8 s、0.5 s、0.22 s、0.15 s、0.06 s.

表3为不同方法在不同剖分网格大小下数值解与解析解的误差的二范数.同等自由度下,优化矩阵下的网格精度相较于集中、一致、折衷矩阵分别提高了89.48%、87.91%、53.57%.图3是四种方法误差二范数与剖分单元边长h的收敛曲线,随着h的减小,优化矩阵下误差二范数总是最小,其次是折衷、一致、集中矩阵.该实验证明优化矩阵在保障正演效率的同时,也能提高正演模拟的精度.

表3 不同网格大小下四种方法的误差L2范数

图3 不同方法下数值解误差2范数收敛曲线

3.1.2 最优网格点精度对比实验

模拟区域激励源与介电参数不变,取单元边长h=0.186 m,根据归一化相速度理论曲线,对应优化矩阵最优单位波长网格点数(G=4.8),剖分单元数量200×200.将源置于模拟区域中心,接收点坐标(27.714 m,9.486 m).采样时间为0.1 ns,采样点数为3000,加载10层常规PML,对模拟区域进行频率域正演并利用傅里叶逆变换至时间域,得到相应波场快照及单道波形,四种方法的计算时间分别都约为190 s.

图4为100 ns时刻的波场快照图,黑色箭头表示数值频散发生的区域.其中图4a黑色实心圆表示激励源的位置,黑色倒三角表示接收点位置.为方便对比四种方法的频散特征,本文将四种方法在100 ns时刻的1/4波形拼接为一个波场快照.图4b表示为四种方法与解析解的残差对比图.此时波未到边界处不会产生反射,且图4a中的频散现象属于正演方法所带来的数值频散而非色散介质产生的色散现象.数值频散效应在一致矩阵和集中矩阵中比较明显:一致矩阵方法的频散效应超前,集中矩阵方法的频散现象滞后.折衷矩阵的频散效应既有超前也有滞后,但相比前两种方法,其频散程度较小.而优化矩阵相较于其他方法未见明显的数值频散现象.

图4 (a) 在均匀频散介质下100ns时刻不同方法的波场快照图;(b) 不同方法与解析解的绝对误差

实际上,上述特征是由归一化相速度与1的关系所决定的:归一化相速度大于1时会出现数值频散超前的现象,而小于1时会出现频散滞后的现象.从图2中可知,一致矩阵和集中矩阵的曲线分别在大于1与小于1的区域分布,折衷矩阵为前两者的线性组合,其归一化相速度曲线在1两侧都有分布;而优化矩阵的归一化相速度曲线与1的误差极小,数值频散现象并不明显.

图5a为接收点处四种方法数值解与解析解的单道波形图,图5b则是四种方法对应的单道波形与解析解的绝对误差.波在135 ns左右到达,一致与集中矩阵方法的频散效应分别在波形到达前和到达后较为明显,且二者的波形与解析解相差较大.折衷矩阵方法的单道波形整体与解析解拟合较好,但在小图以及135 ns前的波形中可以看出,仍然存在一定的频散现象.红色虚线所代表的优化矩阵方法,无论是整体还是局部,其波形都与解析解最为接近,从误差图中可以明显看出其压制数值频散的优势.

图5 (a) G=4.8网格剖分下均匀频散介质下不同方法单道波形图; (b) 不同方法下数值解与解析解绝对误差

值得注意的是,从图5b中可知,一致、集中和折衷矩阵方法的波形在220 ns后还有反射波能量残留.实际上该反射波是由数值频散产生的,可见传统方法加载10层的常规PML的吸收效果并不理想,而优化矩阵方法不仅能够压制数值频散,也能适应常规PML吸收边界.

3.2 吸收边界对比

受模拟区域介电参数以及剖分单元数量的限制,正演模拟时难以将G控制在最优网格点,需要一种更具普适性的吸收边界.本节考虑G不等于最优网格点的剖分下,讨论EPML与常规PML的吸收效果对比.选择中心频率为100 MHz的雷达波脉冲激励源,模拟区域8 m×8 m,剖分单元大小h=0.04 m,模拟区域介电参数如表4所示.左侧一半设置为Debye介质模型,对应G约为22.3;右侧一半介质设置为Cole-Cole模型,对应G约为20.7.将源置于模拟区域中心,左右接收点坐标分别为(2.12 m,2.04 m)和(5.96 m,2.04 m).对模拟区域进行频率域正演并利用傅里叶逆变换至时间域得到相应波场快照及其单道波形.

表4 模型介电参数

图6为PML与EPML吸收边界在58 ns时刻的吸收效果对比图,黑色实心圆表示激励源的位置,黑色倒三角表示接收点位置,黑色矩形框表示模拟区域边界,黑色虚线表示介质分界面.加载常规PML层厚度为0.2 m时,波形在边界处有明显反射波.0.4 m厚度的常规PML与0.2 m的EPML下的波形,其反射波能量明显减弱.从表5可以看出,增加常规PML的厚度,会牺牲正演运行时间以及内存占用;而EPML只需加载0.2 m便可达到0.4 m厚度常规PML边界的吸收效果.

表5 不同吸收边界运行参数

图6 (a)—(c) 不同吸收边界厚度下的波场快照; (d)—(f) 对应吸收边界残差对比图

图7为不同吸收边界条件下在接收点处的单道波形.图中可见加载3种不同PML的直达波都能与参考解较好地吻合.为了分析细微的差别,将图7中加黑色框线67～148 ns的单道雷达波形进行放大得到小图,在不同的色散介质中,0.2 m的EPML与0.4 m的常规PML二者吸收效果相近,且都能与参考解较好地拟合;而0.2 m的常规PML下的反射波能量扰动十分明显.上述结果表明EPML在能够在保证吸收效果的同时,也能提高模拟效率.

图7 R1 (a) 和R2 (b) 接收点处不同吸收边界厚度下的单道波形

3.3 复杂色散介质模型正演

图8为Debye色散介质的城市道路病害模型,模拟区域2 m×0.8 m,采用边长0.0125 m的规则四边形剖分单元,剖分网格160×60.激励源主频900 MHz,放置在空气层中,激发接收位置如图中所示.第一层为随机双相介质混凝土层(冯德山和王珣,2016),背景介质为水泥砂浆,其中骨料散射体占比45%,最大和最小骨料粒径分别为0.02 m和0.0125 m,该层水平位置0.75 m处存在一个宽度为0.0125 m贯通空气裂隙.第二层为色散土壤,此层存在一个大小0.2 m×0.1 m的矩形色散异常体,中心位置为(1.6 m,0.3 m).第三层为埋有两根PVC管的均匀介质,PVC管壁厚度0.0125 m,内径0.1 m,左侧PVC管中心位置为(0.5 m,0.6 m),管内充盈液体;右侧为PVC空管,中心位置为(1.5 m,0.6 m).模型具体参数如表6所示.对模拟区域进行频率域正演,并采用傅里叶逆变换至时间域,可得到相应不同时刻波场快照及剖面图.

表6 模型介电参数表

图8 城市道路病害模型图

将激励源置于空气层中心,不同算法正演结果如图9所示.图9(a、c、e)为加载常规PML下一致矩阵的FEFD算法下波场快照,图9(b、d、f)为加载EPML下优化矩阵的FEFD算法下的波场快照,二者吸收边界厚度相同.电磁波在5.2 ns经过色散矩形异常体,在6.5 ns时到达PVC空气管并开始产生反射波,7.8 ns时刻经过PVC左管产生反射波.如图中黑色箭头所示,图9(a、c、e)中,一致矩阵下的正演结果出现了超前的数值频散,这验证了前文的分析,而优化矩阵下的正演结果未出现明显的数值频散现象.由于加载了EPML吸收边界,在7.8 ns时刻的优化矩阵算法正演结果未出现反射波,而一致矩阵算法波场快照出现了反射波,验证了EPML良好的吸收性能.

图9 (a)、(c)、(e) 常规算法下不同时刻波场快照; (b)、(d)、(e) 优化算法下不同时刻波场快照

将发射天线和接收天线均放在空气层,采用自激自收的方式,偏移距0.0125 m,每次同步移动0.0125 m,共记录160道数据,得到传统算法与优化算法下正演剖面如图10所示.

图10 (a) 常规算法正演剖面图; (b) 优化算法下正演剖面图

从左右剖面图都可看出,混凝土层和色散土壤层的边界是清晰的,而色散土壤层与均匀介质层的边界较为模糊.在混凝土层,由于随机双相介质中骨料的影响,雷达波散射较严重,波形发生扭曲,出现较多的干扰杂波,同时在该层水平位置0.75 m处可以看到由空气裂隙产生的绕射波.色散土壤层中可观察到弧顶在水平位置1.6 m处的矩形异常体的双曲线反射波.均匀介质层中,左右两侧PVC管的双曲线反射波,其弧顶分别在水平位置0.5 m、1.5 m处.右侧的PVC空管的双曲线反射波并不完整,这是由于空管上方的色散矩形体对雷达波的衰减较强,导致波到达空管处的能量较弱.

一致矩阵算法下的剖面图10a中黑色箭头表示数值频散较为明显的地方,且从第二三层界面处可以看出图中的相位、幅值也发生了一定变化.从红色箭头可以看出,在常规PML下可以看到明显的角点反射,而优化矩阵下算法由于加载了EPML,在边缘处未出现明显的反射波.

4 结论

针对常规FEFD算法在色散介质中存在的数值频散问题,提出了一种基于EPML的最优系数FEFD算法.该算法采取归一化相速度与1的误差最小策略,通过最小二乘法求解目标函数,实现了仅需三个优化参数组合下对有限元刚度矩阵与质量矩阵的优化.实验结果表明:优化矩阵的归一化相速度在单位波长仅需4.8个网格点下便可达到误差小于0.2%的精度,而传统方法不仅需要更多的网格点,且精度更低.通过与传统方法的频散特征对比分析,验证了优化矩阵在压制频散方面优良特性.EPML吸收边界条件的引入,简化了吸收参数优化过程,5层EPML即可达到10层常规PML的吸收效果,能够有效节约计算成本,提高色散介质GPR频率域正演模拟精度与效率.

致谢衷心感谢编辑以及审稿专家提出的建设性意见,对文章质量的提高有很大帮助.