CN110210129B - 自适应有限元gpr频率域正演方法 - Google Patents

Abstract

本发明公开了一种自适应有限元GPR频率域正演方法，包括设定正演方法的参数；计算当前网格的数值解和局部误差指示值；标记需细化的单元集合；对标记网格进行细化；进行数值求解得到最终的自适应有限元GPR正演结果。本发明提供的这种自适应有限元GPR正演方法，在后验误差估计策略的基础上，采用自适应有限元方法开展GPR频率域正演，将EPML边界条件引入到频率域正演计算当中，简化了参数优化过程，使得吸收层厚度较小情况下便能达到较好的吸收效果，能够降低节点自由度，将很大程度上节约计算成本。

Description

自适应有限元GPR频率域正演方法

技术领域

本发明具体涉及一种自适应有限元GPR频率域正演方法。

背景技术

探地雷达(Ground Penetrating Radar，GPR)是利用高频电磁波对地下结构或者物体内部不可见目标体进行定位的一种无损探测技术(Daniels，2004；曾昭发，2010)。通过GPR正演，可以加深对GPR传播规律与反射剖面的认识，提高雷达数据解释水平(Giannopoulos,2005；Irving&Knight，2006)。GPR的数值模拟可分为时间域与频率域两大类。其中时间域模拟方法应用较为广泛，包括时域有限差分(Finite Difference TimeDomain，FDTD)(Chen&Huang,1998；Cassidy&Millington,2009；Diamanti&Giannopoulos，2009；Lin et al.，2012；Zhang et al.，2016；Warren et al.，2016)、时域有限单元(Finite Element Time Domain，FETD)(Feng&Wang,2017；Feng et al.,2018；Wang&Dai，2013)、伪谱法(Pseudo-spectral Method，PSM)(Liu&Fan,1999；Huang et al.,2010)、时域间断有限元(Discontinuous Galerkin Finite Element Time Domain，DGFETD)(Lu,etal.,2005)、龙格库塔法(Fang&Lin,2012)等算法。而频率域正演是频率域全波形反演的基础(Keskinen et al.，2017；Ren et al.，2018；Nilot et al.，2018)，通过合理选取少数频率分量数据的反演，即可达到与时间域反演一致的反演结果(Pratt&Worthington，1990)，且频率域正演的各个频率分量相对独立，易于并行运算，无时间频散，频带选取灵活，因此，开展快速、有效的频率域正演对提高反演的计算效率十分重要。

目前，频率域正演方法主要有频域有限差分(Finite Difference FrequencyDomain，FDFD)(Shin&Fan，2012)与频域有限单元(Finite Element Frequency Domain，FEFD)(Jin，1993)两大类。其中，FDFD具有思路清晰、编程简单的优点，但不能与非结构化网格结合，对不规则目标体、复杂的物性分界面的拟合效果不好。而FEFD能够与非结构化网格结合，可以较好地拟合目标体的形状，但针对一些复杂隧道衬砌裂隙病害、在地质应力源点、断层断点等局部大梯度问题，若网格剖分较粗则达不到精度要求，若网格剖分过细则易导致过采样，增加了计算内存，降低了计算效率。为了兼顾运算效率与计算精度，讲究网格节点与单元密度合理分配的自适应有限元法(Adaptive Finite Element Method，AFEM)成为最优方案(Key&Weiss，2006；Liu et al.，2018)，它能根据模拟对象特点，采用后验误差估计算法，由粗及细动态调节网格大小及形函数的阶数，实现网格疏密、形函数阶次自动调控，从而有利于捕捉到解的局部突变特征，实现GPR正问题的高效求解。自适应有限元目前应用非常广泛，Key等学者使用非结构化网格实现了海洋电磁法(Key&Weiss,2006；Frankeet al.,2007；Li&Pek,2008；Li&Key,2007；Key&Ovall,2011)的自适应有限元正演模拟。自适应有限元目前主要分为3类，仅调整单元分布密度的h型、仅改变形函数阶次的p型、同时调整单元分布密度及改变形函数阶次的h-p型。本文选取易于实现的h型自适应策略。

将自适应有限元应用于频率域GPR正演时，截断边界的处理也是必不可少的关键研究内容。在所有的边界条件当中，完美匹配层(PML)(Berenger，1994，1996)无疑为目前应用最广泛、效果最好的边界条件。但传统Berenger PML边界条件提出时主要针对FDTD，且为了取得最佳吸收效果，需要对吸收边界厚度、层数、反射系数等多个参数进行多次调试取优(Gedney，1996)，计算成本较大，并且对不同的计算条件或不同模型，需要重新进行最优吸收参数的选取，此过程相对繁琐。

发明内容

本发明的目的在于提供一种计算速度更快、参数优化计算成本更低且吸收效果更好的自适应有限元GPR频率域正演方法。

本发明提供的这种自适应有限元GPR频率域正演方法，包括如下步骤：

S1.设定正演方法的参数；

S2.根据步骤S1设定的参数，计算当前网格的数值解；

S3.计算局部误差指示值；

S4.标记需细化的单元集合；

S5.对步骤S4标记的网格进行细化；

S6.在步骤S5进行细化后的网格基础上，进行数值求解，从而得到最终的自适应有限元GPR正演结果。

步骤S1所述的参数包括正演频率、模型参数和初始网格。

所述的初始网格，具体为采用Delaunay三角形网格，并对网格进行剖分从而得到非结构化网格剖分结果。

步骤S2所述的计算当前网格的数值解，具体为采用如下步骤进行计算：

A.采用FEFD算法进行计算网格的数值解；

B.采用精确MPL(Exact PML，EPML)处理边界条件。

步骤S3所述的计算局部误差指示值，具体为采用后验误差估计策略计算所有节点处的局部误差指示值。

步骤S4所述的标记需细化的单元集合，具体为将所有的局部误差指示值按照数值从大到小的顺序进行排列，确定并标记局部误差指示值大于设定阈值的三角单元集合，从而得到需细化的单元集合。

步骤S5所述的对步骤S4标记的网格进行细化，具体为按照标记的单元集合中局部误差指示值从大至小遍历所有三角单元，并进行网格细化；同时，若细化后的网格满足给定的精度要求，则执行步骤S6，否则返回步骤S2。

本发明提供的这种自适应有限元GPR频率域正演方法，在后验误差估计策略的基础上，采用自适应有限元方法开展GPR频率域正演，将EPML边界条件引入到正演计算当中，简化了参数优化过程，使得吸收层厚度较小情况下便能达到较好的吸收效果，能够降低节点自由度，将很大程度上节约计算成本。

附图说明

图1为本发明方法的方法流程示意图。

图2为本发明方法的超收敛恢复技术计算局部误差的示意图。

图3为本发明方法的不同PML方法的误差对比图。

图4为本发明方法的100MHz频率下解析解、加载常规Berenger PML及加载EPML的数值解波场对比图。

图5为本发明方法的100MHz频率下不同位置处理论解与解析解的对比图。

图6为本发明方法的100MHz频率下数值解与解析解的波场值绝对误差分布图。

图7为本发明方法的100MHz频率下加载EPML的数值解与解析解之间的误差对比图。

图8为本发明方法的复杂隧道衬砌裂隙病害模型及时间域正演剖面图。

图9为本发明方法的x＝0.2m处时间域单道波形对比图。

图10为本发明方法的激发点在0.5m处不同频率下的自适应网格示意图。

图11为本发明方法的激发点在0.5m处不同频率下的频率域波场分布图。

图12为本发明方法的激发点在0.5m处使用AFEFD+EPML正演出的频率域波场转换至时间域后，不同时刻的时间域波场快照示意图。

具体实施方式

如图1所示为本发明方法的方法流程示意图：本发明提供的这种自适应有限元GPR频率域正演方法，包括如下步骤：

S1.设定正演方法的参数，包括正演频率、模型参数和初始网格；

初始网格采用Delaunay三角形网格，并对网格进行剖分从而得到非结构化网格剖分结果；

S2.根据步骤S1设定的参数，计算当前网格的数值解；具体为采用如下步骤进行计算：

A.采用FEFD算法进行计算网格的数值解；

根据电磁波传播理论(Taflove&Brodwin，1975)，雷达波在有耗媒质中的传播遵循Maxwell方程组，其二维TM波的频率域表示形式为：

其中E_z为频域率电场强度(V/m)，J_z为频率域电流密度(A/m²)，μ为磁导率(H/m)，ε为介电常数(F/m)，σ为电导率(S/m)，ω＝2πf为角频率(rad/s)，e^iωt为时谐因子，i为虚数单位。

假定介电常数ε和电导率σ是实数，再引入表示波的传播和衰减的复介电常数，进一步简化上式可得

考虑到GPR探测的目标对象多为非磁性介质，取其磁导率μ＝μ₀μ_r＝μ₀＝4π×10^-1H/m(μ_r＝1)。已知自由空间的介电常数为ε₀＝8.854×10^-12F/m，相对介电常数表达式为ε_r＝ε/ε_r，则

设

为自由空间波数，将频率域波数的定义为：

可以进一步改写成为如下形式：

其中，μ_r＝1根据电磁场理论，将时谐场的各向异性介质的频率域波动方程表示为：

其中，α_x＝α_y＝1/μ_r，β＝μ_r；

对上式进行积分处理，根据分部积分原理，该方程相应的弱解形式为：

其中，

为权函数，而由于计算边界为完美电导，二维TM波的频率域表示式的右边为0，因此，弱解形式可简化为：

计算区域划分单元后，积分可以写成各个单元积分之和：

其中，e表示单元标号，N_e为单元总个数，Ω^e表示在第e个单元上的积分。若采用三角形单元进行网格剖分离散，并将场值函数u在各单元上应用线性插值：

根据Galerkin法，令权函数

i＝1,2,3，对于任意三角形单元，有如下单元积分：

若假设单元足够小，则介质参数在单元内近似为均匀的，与积分无关，则可改写为：

其中

为单元刚度矩阵，

为单元质量矩阵，f_i ^e为单元源向量：

再将单元场向量、激励源及单元系数矩阵扩展成整体矩阵，有如下矩阵表达：

-Ku+Mu＝f

其中系统矩阵K，M，和源向量f分别由它们的单元级别矩阵K^e，M^e，和源向量f^e组合而成，u是模拟波场向量；

B.采用EPML处理边界条件；

在有限区域进行正演时，截断边界的处理也是必不可少的关键研究内容。在所有的边界条件当中，完美匹配层(Perfect Matching Layer,PML)(Berenger,1994)无疑为目前应用最广泛、效果最好的边界条件。

回顾有限元控制方程式，若引入PML边界条件，在PML吸收边界区域，式中a_x＝s_y/s_x，a_y＝s_x/s_y，β＝s_xs_y，s_x和s_y为坐标拉伸变量，则上式可具体表示为：

传统Berenger PML边界条件提出时主要针对FDTD，为了取得最佳吸收效果，对不同的计算条件或不同模型，均需要重新对吸收边界厚度、层数、反射系数等多个参数进行最优化的选取，此过程相对繁琐，计算成本较大。本发明使用Bermúdez等人2004年提出的一种基于无限吸收函数的精确完全匹配层(Exact Perfectly Matched Layer，EPML)，其吸收函数在垂直于PML厚度增加方向的积分无穷大，使得最终得到的PML边界条件下的解与无界区域下的真实解一致(Bermúdez et al.，2004；Bermúdez et al.，2007；Rabinovich et al.，2010；Cimpeanu et al.,2015)：

∫σ_i(ξ)dξ＝+∞i＝x,y

坐标拉伸变量定义为s_x＝1-iσ_x/ω，s_y＝1-iσ_y/ω，其中k为GPR波动方程频率域波数，以x方向为例，电导率σ_x具体计算公式为：

其中L_PML为PML层的厚度，l_x表示PML内部单元中心到模拟区域边界的距离，γ为自由变量；

EPML提供了一种简洁有效的PML方案，将这种精确PML方法引入到GPR频率域正演当中，将极大简化吸收参数的最优化过程，理论上能够提高正演效率，对于正反演的开展有着重要的研究意义；

S3.计算局部误差指示值；

采用Zienkiewicz-Zhu后验误差估计策略(Zienkiewicz&Zhu，1992a；Zienkiewicz&Zhu，1992b)，在当前网格数值解的基础上，该策略被用来计算出所有节点处的局部误差指示值η；

已知三角网格剖分节点坐标及波场值，设网格中任一节点P为Zienkiewicz-Zhu后验误差估计的超级收敛点，对于围绕P点的所有三角单元而言，都可以得到一个线性近似梯度。即对于任意三角单元，内部的场值分布满足线性插值关系，可以求得每个三角单元的一阶偏导数：

其中Δ表示任一三角单元的面积，已知围绕P点的所有三角单元的一阶偏导数，则超收敛点P点处的一阶偏导数定义为：

其中n_p表示包含节点P的三角单元个数，Δ_i表示三角单元i的面积，u_i表示三角形单元i的场值。由此可求得区域所有超级收敛点处的一阶偏导数(梯度值)u_x和u_y。再对每个节点的梯度值分别求x和y偏导数，可得到4个二阶偏导数：

至此，可定义节点P的局部误差指示值为

η＝(u_xx+u_xy+u_yy+u_yx)²

对三角剖分的所有节点进行上述处理，则完成了有限元计算区域的后验误差估计；

S4.标记需细化的单元集合；

将所有的局部误差指示值按照数值从大到小的顺序进行排列，引入

提出的标记策略(

1996)，确定局部误差指示值大于设定阈值的三角单元集合；将满足细化条件的所有节点进行一次性标记，为下一步网格细化做准备；

将所有的局部误差指示值按照数值从大到小的顺序进行排列，进而确定局部误差指示值大于设定阈值的三角单元集合。若定义区域内部所有误差指示值的总和为：

其中η_τ为局部误差指示值，T表示全体三角单元，τ表示局部三角单元。为了在所有的三角单元集合中选取大于设定阈值的单元进行标记，引入

提出的标记策略(

1996)，并定义标记单元集合为M：

其中θ为标记控制参数，若θ较大，则选中的区域更大，网格更精细化，若θ较小，则选中的区域较小，往往更容易圈定异常体，易获得最佳网格。参考Cheng&Zhzng(2006)的细化方案，θ的取值范围控制在0.2至0.5之间，将满足细化条件的所有节点进行一次性标记，为下一步网格细化做准备；

S5.对步骤S4标记的网格进行细化；

如果某一节点被标记，则包含该节点的所有单元均需要进行细化，为方便执行，具体实施时，按照标记的单元集合中局部误差指示值从大至小来遍历所有三角单元进行网格细化。网格细化方法参照Chen&Zhang在2006年提出的细化方案(Chen&Zhang，2006)。当遍历完所有三角形，则可实现标记集合中的一次网格细化。若细化后的网格满足给定的精度要求，执行步骤S6，否则返回步骤S2，重复此过程，并根据需要调节网格细化的次数及细化范围，即可达到网格疏密自动调控的自适应网格剖分的目的；

S6.在步骤S5进行细化后的网格基础上，进行数值求解，从而得到最终的自适应有限元GPR正演结果。

以下结合两个具体实施例对本发明方法进行进一步说明：

a.均匀介质模型EPML数值实验

为分析EPML在GPR频率域正演中的吸收效果，选取100MHz、400MHz、900MHz三个典型频率，在均匀介质区域进行GPR二维FEFD正演模拟，PML吸收边界分别使用常规BerengerPML和EPML，比较两种PML边界条件下均匀介质区域线电流辐射场的数值解与解析解之间的误差。已知均匀介质区域线电流辐射场的解析式为(Harrington，1968)：

其中，

(·)为第二类零阶Hankel函数，k为频率域波数，ρ为距线电流源的距离。定义误差：

其中，N为模拟区域内节点数目，u_TRUE为与解析解相对应的值，u_FEFD为与数值解相对应的值，由于频率域解为复数，u_TRUE和u_FEFD可对应取解的实部、虚部和绝对值；

模拟在规则三角网格剖分情况下进行，计算出PML区域中不同节点个数(即不同PML厚度)条件下的误差。激励源在模拟区域中心，模拟区域背景介质的相对介电常数为3，电导率为5.0×10^-5S/m。其中100MHz频率下区域大小为8.0m×8.0m，离散网格间距取0.05m；400MHz频率下区域大小为4.0m×4.0m，离散网格间距取0.01m；900MHz频率下区域大小为2.0m×2.0m，离散网格间距取0.005m。

图3为本发明方法的不同PML方法的误差对比图：Ue表示加载EPML，Ug表示加载常规Berenger PML，Re、Im、Abs分别表示解的实部、虚部、绝对值。(a)天线频率100MHz时的误差对比图；(b)天线频率400MHz时的误差对比图；(c)天线频率900MHz时的误差对比图。图3(a)、(b)和(c)展示了天线频率分别为100MHz、400MHz和900MHz时，常规Berenger PML吸收边界和无参数PML吸收边界条件下的数值解与解析解的实部、虚部和绝对值的误差随着PML中节点数增加的变化情况。误差Error的计算公式已给出。从图3中可看出，在误差较小时，EPML的节点个数比常规Berenger PML的节点个数少；即使PML区域节点个数小于4，EPML的误差也几乎接近等于最小误差；相反地，常规Berenger PML的误差更大且受PML厚度的影响较大。详细来分析，以图3(a)中线Abs(Ue)和Abs(Ug)为例，当误差接近最小值(大约0.002)时，EPML只需要6个节点的PML厚度，而常规Berenger PML则需要20个节点的PML厚度，并且类似的特征我们也可以从图3(b)和图3(c)中总结得到。此例说明了在正演中使用EPML，可以减少PML区域的节点数，在更小的PML厚度条件下达到较好的吸收效果，从而减少计算量，提高正演效率。

以天线频率为100MHz、PML包含4个节点(即PML厚度为0.2m)的正演模拟为例，直观比较常规Berenger PML和EPML边界条件下的正演结果，区域大小为8.0m×8.0m，离散网格间距取0.05m。

图4为本发明方法的100MHz频率下解析解、加载常规Berenger PML及加载EPML的数值解波场对比图：(a)解析解的实部；(b)解析解的虚部；(c)解析解的绝对值；(d)加载常规Berenger PML的数值解的实部；(e)加载常规Berenger PML的数值解的虚部；(f)加载常规Berenger PML的数值解的绝对值；(g)加载EPML的数值解的实部；(h)加载EPML的数值解的虚部；(i)加载EPML的数值解的绝对值。图4从上至下分别为解析解、常规Berenger PML条件下的数值解以及EPML条件下的数值解的波场图，从左至右分别为解的实部、虚部和绝对值。分析图4(a)-(c)中解析解的波场图可知，频率域GPR二维波场呈圆环状分布，能量随着距离源的位置的增加而逐渐降低。图4(g)-(i)的EPML正演波场分布情况与解析解相似，而图4(d)-(f)的常规Berenger PML数值解则与解析解之间存在较大差异，计算区域和吸收层之间的虚假反射使得波场的圆环状分布变得紊乱，这一点从图4(f)解的绝对值更容易看出。

图5为本发明方法的100MHz频率下不同位置处理论解与解析解的对比图：其中实线表示理论解，点虚线表示加载常规Berenger PML的数值解，虚线表示加载EPML的数值解。(a)y＝0m处解的实部；(b)y＝2m处解的实部；(c)y＝4m处解的实部；(d)y＝0m处解的虚部；(e)y＝2m处解的虚部；(f)y＝4m处解的虚部；(g)y＝0m处解的绝对值；(h)y＝2m处解的绝对值；(i)y＝4m处解的绝对值。图5为y＝0m，y＝2m，y＝4m处解析解及加载两种PML边界条件下的频率域单道波形对比图。图中可见，除激励源加载处，EPML条件下的数值解与解析解在计算区域几乎保持着高度吻合，且雷达波在进入吸收层后几乎呈直线下降为0值；在靠近激励源的地方解的数值较大，常规Berenger PML条件下的数值解与解析解的差异很难分辨，而在离激励源较远的区域，其与解析解的误差较大，波场存在着变形、震荡，仔细观察可发现雷达波进入吸收层后并没有被完全吸收，与解析解的偏差比EPML大。

图6为本发明方法的100MHz频率下数值解与解析解的波场值绝对误差分布图：(a)加载常规Berenger PML的数值解的实部减解析解的实部；(b)加载常规Berenger PML的数值解的虚部减解析解的虚部；(c)加载常规Berenger PML的数值解的绝对值减解析解的绝对值；(d)加载EPML的数值解的实部减解析解的实部；(e)加载EPML的数值解的虚部减解析解的虚部；(f)加载EPML的数值解的绝对值减解析解的绝对值。图6给出了加载两种不同PML边界条件下的数值解与解析解的残差分布图。图6(d)-(f)为加载EPML边界条件的数值解与解析解的残差波场分布图，它显然比图6(a)-(c)中的常规Berenger PML方法的残差的大小、残差分布均呈现出优越性，吸收效果更好，这也与图4和图5中得到的结论相吻合。

除了与PML厚度紧密相关的PML中节点数量，在规则网格中，影响PML方法所产生的误差大小的另一个重要因素就是离散网格间距。仍以天线频率为100MHz为例，设置不同的PML厚度及离散网格间距，模拟区域大小为8.0m×8.0m，进行基于EPML方法的GPR二维FEFD正演。

图7为本发明方法的100MHz频率下加载EPML的数值解与解析解之间的误差对比图。图7为模拟区域中数值解与解析解的绝对值误差示意图。每一条带有不同标记符号的线段描绘了相同离散网格间距条件下解的绝对值误差随着PML厚度的改变的变化情况。不同的离散网格间距的误差虽然不尽相同，但其趋势相近，即随着PML厚度的增加，误差逐渐减小，最终保持不变。比较相同PML厚度条件下，带有相同标记符号的线段表示的不同离散网格间距条件下的误差变化曲线，可发现，随着离散网格间距的减小，误差呈现下降趋势。但与此同时，离散网格间距减小，计算时间相应增加。例如当PML厚度为0.2m时，离散网格间距分别取0.0125m、0.02m、0.025m、0.0333m、0.05m和0.1m时，计算时间分别为74.926s、11.739、5.643s、2.276s、0.690s和0.158s。当网格间距较小时，增加离散网格单元数目对数值解的优化会使得计算成本急剧增大，继续加密网格的意义不大(例如，当离散网格间距足够小为0.0125m时，误差小于0.001，PML厚度的增加对误差的影响很小)，故而在实际计算中应当选择合适的网格大小，兼顾计算精度和速率。

b.复杂衬砌裂隙病害AFEFD正演模拟

在隧道施工过程中，因地质条件复杂、施工环境恶劣、施工工艺等问题，容易导致衬砌中形成不规则裂隙病害。

图8为本发明方法的复杂隧道衬砌裂隙病害模型及时间域正演剖面图。(a)包含三角剖分结果的模型图，白色虚框为图9中的放大区域；(b)GPR时间域正演剖面图，白色虚线表示图8中单道曲线所在位置，椭圆虚线框出的区域A为裂隙反射界面，矩形虚线框出的区域B为下方不规则围岩界面的反射界面。为了更好地认识衬砌裂隙模型的雷达剖面特征，建立图8(a)所示的衬砌裂隙模型，采用非结构化网格对该模型区域进行离散。模型中部区域为混凝土，相对介电常数为9，电导率为1.0×10^-5S/m；不规则起伏面下方为围岩，相对介电常数为7，电导率为1.0×10^-3S/m；模型中部模拟裂隙形状，裂隙内填充空气，相对介电常数为1，电导率为0；模型上方为0.1m厚的空气层。激励源为900MHz的Ricker子波，收发天线位置如图8(a)所示，“·”表示激发天线位置，“×”表示接收天线位置，收发天线距地表0.01m，激发点起始位置为0.01m，收发距为0.02m，两者自左向右同步移动，采样间隔为0.01m，时间步长为0.01ns，时窗长度为15ns，共计录197道雷达数据。图8(b)为多频率正演数据变换至时间域的雷达记录剖面，图中清晰可见直达波和裂隙产生的反射界面(椭圆虚线框出的区域A)，下方由不规则界面引起的反射界面(矩形虚线框出的区域B)也依稀可辨，同时，正演记录的下方也出现了振幅较大的由裂隙引起的多次反射。

参考上例图3(c)中900MHz频率下的误差变化，可推断EPML只需0.03m厚(图3(c)中为6个节点)便能够达到很好的吸收效果，而常规Berenger PML需要0.1m(图3(c)中为20个节点)才能达到与EPML同等的模拟精度。图9为本发明方法的x＝0.2m处时间域单道波形对比图。图9为x＝0.2m处时间域的单道雷达波形对比图，实线表示加载0.03m厚的EPML，虚线与点线分别表示加载0.1m及0.03m两种不同厚度的常规Berenger PML。图中可加载3种不同PML的自适应有限元的模拟结果在原始图中都能较好地吻合；但是，在中部框线突出显示的8-15ns的放大波形中可以看到细微的差别，图中可见实线表示的0.03m厚的EPML与虚线表示的0.1m厚的常规Berenger PML数据较好地拟合，而0.03m厚的常规Berenger PML与前面两者相比，波形曲线存在一定的偏离。为比较PML厚度引起的计算时间的差异，选取单个频率900MHz，加载单个激励源，给出初始网格剖分情况下的正演时间。三种不同厚度的PML的求解时间分别为0.109s(0.03m EPML)，0.173s(0.1m常规Berenger PML)，0.110s(0.03m常规Berenger PML)。PML厚度为0.03m时，两种PML单次求解的时间大致相同，而PML厚度为0.1m时，计算时间却增加了50％。以上对比说明EPML达到较好的吸收效果时，与常规Berenger PML相比，加载PML的厚度较小，同时计算时间较少。

图10为本发明方法的激发点在0.5m处不同频率下的自适应网格示意图：显示区域为图8(a)中白色虚框标出的区域。(a)400MHz频率下的自适应网格；(b)900MHz频率下的自适应网格；(c)1500MHz频率下的自适应网格。图11为本发明方法的激发点在0.5m处不同频率下的频率域波场分布图：虚线绘出了异常体轮廓。(a)400MHz频率下加载EPML的数值解的实部；(b)400MHz频率下加载EPML的数值解的虚部；(c)400MHz频率下加载EPML的数值解的绝对值；(d)900MHz频率下加载EPML的数值解的实部；(e)900MHz频率下加载EPML的数值解的虚部；(f)1500MHz频率下加载EPML的数值解的绝对值；(g)1500MHz频率下加载EPML的数值解的实部；(h)1500MHz频率下加载EPML的数值解的虚部；(i)1500MHz频率下加载EPML的数值解的绝对值。

图10表示不同频率下的自适应加密网格，图11展示了激励源在(0.5m，-0.01m)处加载时频率为400MHz(a，b，c)、900MHz(d，e，f)和1500MHz(g，h，i)时的正演波场分布情况。随着频率的增加，网格逐渐加密(400MHz、900MHz和1500MHz频率下，离散三角单元数分别为47235,73020和96707)，电磁波波长越来越短，且相比较于时间域，频率域波场的异常反射特征更加明显。

图12为本发明方法的激发点在0.5m处使用AFEFD+EPML正演出的频率域波场转换至时间域后，不同时刻的时间域波场快照示意图：虚线绘出了异常体轮廓。(a)t＝1.2ns时的波场分布图；(b)t＝2.8ns时的波场分布图；(c)t＝3.8ns时的波场分布图；(d)t＝4.9ns时的波场分布图；(e)t＝6.3ns时的波场分布图；(f)t＝7.5ns时的波场分布图。图12为激发点位于(0.5m，-0.01m)处由频率域变换至时间域的不同时刻的波场快照。图12(a)中，电磁波由激发点向外传播，在空气和下方介质处发生了分离，在空气层中传播较快，在下方介质中传播较慢。图12(b)中电磁波继续传播，传播到空气层中的波被上方PML吸收，传播到下方介质中的电磁波遇到裂隙开始发生反射。图12(c)中由裂隙产生的反射波向上传播，原电磁波前继续向外传播，下方波前接触到弯曲界面。图12(d)中由裂隙产生的反射波一部分穿过空气层，速度加快，最终被上方PML吸收，一部分在与空气层接触的地方发生反射；原电磁波前继续向外传播，在弯曲界面处产生的反射波较小，几乎不可分辨。图12(e)中由空气层与下方介质产生的二次反射波向下传播，下方波前接触到弯曲界面，原电磁波前继续向外传播，最终被边界处的PML吸收。图12(f)中由空气层与下方介质产生的二次反射波在裂隙处再次发生反射，波场分离，原始电磁波继续向外传播，并被PML吸收。

结合上述的两个实施例，本发明的优点在于：

a.均匀介质模型EPML数值实验

为分析EPML在GPR频率域正演中的吸收效果，选取100MHz、400MHz、900MHz三个典型频率，在均匀介质区域进行GPR二维FEFD正演模拟，PML吸收边界分别使用常规BerengerPML和EPML，比较两种PML边界条件下均匀介质区域线电流辐射场的数值解与解析解之间的误差。

模拟在规则三角网格剖分情况下进行，计算出PML区域中不同节点个数(即不同PML厚度)条件下的误差。激励源在模拟区域中心，模拟区域背景介质的相对介电常数为3，电导率为5.0×10^-5S/m。其中100MHz频率下区域大小为8.0m×8.0m，离散网格间距取0.05m；400MHz频率下区域大小为4.0m×4.0m，离散网格间距取0.01m；900MHz频率下区域大小为2.0m×2.0m，离散网格间距取0.005m。图3(a)、(b)和(c)展示了天线频率分别为100MHz、400MHz和900MHz时，常规Berenger PML吸收边界和无参数PML吸收边界条件下的数值解与解析解的实部、虚部和绝对值的误差随着PML中节点数增加的变化情况。从图3中可看出，在误差较小时，EPML的节点个数比常规Berenger PML的节点个数少；即使PML区域节点个数小于4，EPML的误差也几乎接近等于最小误差；相反地，常规Berenger PML的误差更大且受PML厚度的影响较大。详细来分析，以图3(a)中线Abs(Ue)和Abs(Ug)为例，当误差接近最小值(大约0.002)时，EPML只需要6个节点的PML厚度，而常规Berenger PML则需要20个节点的PML厚度，并且类似的特征我们也可以从图3(b)和图2(c)中总结得到。此例说明了在正演中使用EPML，可以减少PML区域的节点数，在更小的PML厚度条件下达到较好的吸收效果，从而减少计算量，提高正演效率。

以天线频率为100MHz、PML包含4个节点(即PML厚度为0.2m)的正演模拟为例，直观比较常规Berenger PML和EPML边界条件下的正演结果，区域大小为8.0m×8.0m，离散网格间距取0.05m。图4从上至下分别为解析解、常规Berenger PML条件下的数值解以及EPML条件下的数值解的波场图，从左至右分别为解的实部、虚部和绝对值。分析图4(a)-(c)中解析解的波场图可知，频率域GPR二维波场呈圆环状分布，能量随着距离源的位置的增加而逐渐降低。图4(g)-(i)的EPML正演波场分布情况与解析解相似，而图4(d)-(f)的常规BerengerPML数值解则与解析解之间存在较大差异，计算区域和吸收层之间的虚假反射使得波场的圆环状分布变得紊乱，这一点从图4(f)解的绝对值更容易看出。

图5为y＝0m，y＝2m，y＝4m处解析解及加载两种PML边界条件下的频率域单道波形对比图。图中可见，除激励源加载处，EPML条件下的数值解与解析解在计算区域几乎保持着高度吻合，且雷达波在进入吸收层后几乎呈直线下降为0值；在靠近激励源的地方解的数值较大，常规Berenger PML条件下的数值解与解析解的差异很难分辨，而在离激励源较远的区域，其与解析解的误差较大，波场存在着变形、震荡，仔细观察可发现雷达波进入吸收层后并没有被完全吸收，与解析解的偏差比EPML大。

图6给出了加载两种不同PML边界条件下的数值解与解析解的残差分布图。图6(d)-(f)为加载EPML边界条件的数值解与解析解的残差波场分布图，它显然比图6(a)-(c)中的常规Berenger PML方法的残差的大小、残差分布均呈现出优越性，吸收效果更好，这也与图4和图5中得到的结论相吻合。

除了与PML厚度紧密相关的PML中节点数量，在规则网格中，影响PML方法所产生的误差大小的另一个重要因素就是离散网格间距。仍以天线频率为100MHz为例，设置不同的PML厚度及离散网格间距，模拟区域大小为8.0m×8.0m，进行基于EPML方法的GPR二维FEFD正演。图7为模拟区域中数值解与解析解的绝对值误差示意图。每一条带有不同标记符号的线段描绘了相同离散网格间距条件下解的绝对值误差随着PML厚度的改变的变化情况。不同的离散网格间距的误差虽然不尽相同，但其趋势相近，即随着PML厚度的增加，误差逐渐减小，最终保持不变。比较相同PML厚度条件下，带有相同标记符号的线段表示的不同离散网格间距条件下的误差变化曲线，可发现，随着离散网格间距的减小，误差呈现下降趋势。但与此同时，离散网格间距减小，计算时间相应增加。例如当PML厚度为0.2m时，离散网格间距分别取0.0125m、0.02m、0.025m、0.0333m、0.05m和0.1m时，计算时间分别为74.926s、11.739、5.643s、2.276s、0.690s和0.158s。当网格间距较小时，增加离散网格单元数目对数值解的优化会使得计算成本急剧增大，继续加密网格的意义不大(例如，当离散网格间距足够小为0.0125m时，误差小于0.001，PML厚度的增加对误差的影响很小)，故而在实际计算中应当选择合适的网格大小，兼顾计算精度和速率。

b.复杂衬砌裂隙病害AFEFD正演模拟

在隧道施工过程中，因地质条件复杂、施工环境恶劣、施工工艺等问题，容易导致衬砌中形成不规则裂隙病害。为了更好地认识衬砌裂隙模型的雷达剖面特征，建立图8(a)所示的衬砌裂隙模型，采用非结构化网格对该模型区域进行离散。模型中部区域为混凝土，相对介电常数为9，电导率为1.0×10^-5S/m；不规则起伏面下方为围岩，相对介电常数为7，电导率为1.0×10^-3S/m；模型中部模拟裂隙形状，裂隙内填充空气，相对介电常数为1，电导率为0；模型上方为0.1m厚的空气层。激励源为900MHz的Ricker子波，收发天线位置如图8(a)所示，“·”表示激发天线位置，“×”表示接收天线位置，收发天线距地表0.01m，激发点起始位置为0.01m，收发距为0.02m，两者自左向右同步移动，采样间隔为0.01m，时间步长为0.01ns，时窗长度为15ns，共计录197道雷达数据。图8(b)为多频率正演数据变换至时间域的雷达记录剖面，图中清晰可见直达波和裂隙产生的反射界面(椭圆虚线框出的区域A)，下方由不规则界面引起的反射界面(矩形虚线框出的区域B)也依稀可辨，同时，正演记录的下方也出现了振幅较大的由裂隙引起的多次反射。

参考上例图3(c)中900MHz频率下的误差变化，可推断EPML只需0.03m厚(图3(c)中为6个节点)便能够达到很好的吸收效果，而常规Berenger PML需要0.1m(图3(c)中为20个节点)才能达到与EPML同等的模拟精度。图9为x＝0.2m处时间域的单道雷达波形对比图，实线表示加载0.03m厚的EPML，虚线与点线分别表示加载0.1m及0.03m两种不同厚度的常规Berenger PML。图中可见加载3种不同PML的自适应有限元的模拟结果在原始图中都能较好地吻合；但是，在框线突出显示的8-15ns的放大波形中可以看到细微的差别，图中可见实线表示的0.03m厚的EPML与虚线表示的0.1m厚的常规Berenger PML数据较好地拟合，而0.03m厚的常规Berenger PML与前面两者相比，波形曲线存在一定的偏离。为比较PML厚度引起的计算时间的差异，选取单个频率900MHz，加载单个激励源，给出初始网格剖分情况下的正演时间。三种不同厚度的PML的求解时间分别为0.109s(0.03m EPML)，0.173s(0.1m常规Berenger PML)，0.110s(0.03m常规Berenger PML)。PML厚度为0.03m时，两种PML单次求解的时间大致相同，而PML厚度为0.1m时，计算时间却增加了50％。以上对比说明EPML达到较好的吸收效果时，与常规Berenger PML相比，加载PML的厚度较小，同时计算时间较少。

图10表示不同频率下的自适应加密网格，图11展示了激励源在(0.5m，-0.01m)处加载时频率为400MHz(a，b，c)、900MHz(d，e，f)和1500MHz(g，h，i)时的正演波场分布情况。随着频率的增加，网格逐渐加密：400MHz、900MHz和1500MHz频率下，离散三角单元数分别为47235,73020和96707，电磁波波长越来越短，且相比较于时间域，频率域波场的异常反射特征更加明显。

图12为激发点位于(0.5m，-0.01m)处由频率域变换至时间域的不同时刻的波场快照。图12(a)中，电磁波由激发点向外传播，在空气和下方介质处发生了分离，在空气层中传播较快，在下方介质中传播较慢。图12(b)中电磁波继续传播，传播到空气层中的波被上方PML吸收，传播到下方介质中的电磁波遇到裂隙开始发生反射。图12(c)中由裂隙产生的反射波向上传播，原电磁波前继续向外传播，下方波前接触到弯曲界面。图12(d)中由裂隙产生的反射波一部分穿过空气层，速度加快，最终被上方PML吸收，一部分在与空气层接触的地方发生反射；原电磁波前继续向外传播，在弯曲界面处产生的反射波较小，几乎不可分辨。图12(e)中由空气层与下方介质产生的二次反射波向下传播，下方波前接触到弯曲界面，原电磁波前继续向外传播，最终被边界处的PML吸收。图12(f)中由空气层与下方介质产生的二次反射波在裂隙处再次发生反射，波场分离，原始电磁波继续向外传播，并被PML吸收。

本发明中，基于非结构化网格AFEFD方法和EPML吸收边界条件是主要内容，为了提高GPR的频率域正演效率，重点分析引入的EPML的吸收性能和复杂目标体的自适应模拟结果，一些结论如下：

与基于指数吸收函数的一般PML吸收边界条件相比，本发明中提出的EPML吸收边界条件具有良好的吸收效果，且开发的自适应有限元策略可以自动控制网格密度，这为GPR模拟，特别是针对复杂的地质构造，提供了另一种可行方案。

Claims (6)

1.一种自适应有限元GPR频率域正演方法，包括如下步骤：

S1.设定正演方法的参数；

S2.根据步骤S1设定的参数，计算当前网格的数值解；具体为采用如下步骤进行计算：

A.采用频域有限单元算法进行计算网格的数值解；

根据电磁波传播理论，雷达波在有耗媒质中的传播遵循Maxwell方程组，其二维TM波的频率域表示形式为：

其中E_z为频域率电场强度V/m，J_z为频率域电流密度A/m²，μ为磁导率H/m，ε为介电常数F/m，σ为电导率S/m，ω＝2πf为角频率rad/s，e^iωt为时谐因子，i为虚数单位；

假定介电常数ε和电导率σ是实数，再引入表示波的传播和衰减的复介电常数，进一步简化上式可得

考虑到探地雷达探测的目标对象多为非磁性介质，取其磁导率μ＝μ₀μ_r＝μ₀＝4π×10^{- 1}H/m；μ_r＝1；已知自由空间的介电常数为ε₀＝8.854×10^-12F/m，相对介电常数表达式为ε_r＝ε/ε_r，则

设

为自由空间波数，将频率域波数定义为：

进一步改写成为如下形式：

其中，μ_r＝1根据电磁场理论，将时谐场的各向异性介质的频率域波动方程表示为：

其中，α_x＝α_y＝1/μ_r，β＝μ_r；

对上式进行积分处理，根据分部积分原理，该方程相应的弱解形式为：

其中，

为权函数，而由于计算边界为完美电导，二维TM波的频率域表示式的右边为0，因此，弱解形式可简化为：

计算区域划分单元后，积分可以写成各个单元积分之和：

其中，e表示单元标号，N_e为单元总个数，Ω^e表示在第e个单元上的积分；若采用三角形单元进行网格剖分离散，并将场值函数u在各单元上应用线性插值：

根据Galerkin法，令权函数

i＝1,2,3，对于任意三角形单元，有如下单元积分：

假设介质参数在单元内为均匀的，与积分无关，则可改写为：

其中

为单元刚度矩阵，

为单元质量矩阵，f_i ^e为单元源向量：

再将单元场向量、激励源及单元系数矩阵扩展成整体矩阵，有如下矩阵表达：

-Ku_u+Mu_u＝f

其中系统矩阵K，M，和源向量f分别由它们的单元级别矩阵K^e，M^e，和源向量f^e组合而成，u_u是模拟波场向量；

B.采用精确完美匹配层处理边界条件；

引入完美匹配层边界条件，在完美匹配层吸收边界区域，式中a_x＝s_y/s_x，a_y＝s_x/s_y，β＝s_xs_y，s_x和s_y为坐标拉伸变量，则上式可具体表示为：

本发明使用精确完美匹配层，其吸收函数在垂直于完美匹配层厚度增加方向的积分无穷大，使得最终得到的完美匹配层边界条件下的解与无界区域下的真实解一致：

∫σ_i(ξ)dξ＝+∞ i＝x,y

坐标拉伸变量定义为s_x＝1-iσ_x/ω，s_y＝1-iσ_y/ω，其中k为探地雷达波动方程频率域波数，以x方向为例，电导率σ_x具体计算公式为：

其中L_PML为完美匹配层的厚度，l_x表示完美匹配层内部单元中心到模拟区域边界的距离，γ为自由变量；

S3.计算局部误差指示值；

S4.标记需细化的单元集合；

S5.对步骤S4标记的网格进行细化；

S6.在步骤S5进行细化后的网格基础上，进行数值求解，从而得到最终的自适应有限元探地雷达正演结果。

2.根据权利要求1所述的自适应有限元GPR频率域正演方法，其特征在于步骤S1所述的参数包括正演频率、模型参数和初始网格。

3.根据权利要求2所述的自适应有限元GPR频率域正演方法，其特征在于所述的初始网格，具体为采用Delaunay三角形网格，并对网格进行剖分从而得到非结构化网格剖分结果。

4.根据权利要求3所述的自适应有限元GPR频率域正演方法，其特征在于步骤S3所述的计算局部误差指示值，具体为采用后验误差估计策略计算所有节点处的局部误差指示值。

5.根据权利要求4所述的自适应有限元GPR频率域正演方法，其特征在于步骤S4所述的标记需细化的单元集合，具体为将所有的局部误差指示值按照数值从大到小的顺序进行排列，确定并标记局部误差指示值大于设定阈值的三角单元集合，从而得到需细化的单元集合。

6.根据权利要求5所述的自适应有限元GPR频率域正演方法，其特征在于步骤S5所述的对步骤S4标记的网格进行细化，具体为按照标记的单元集合中局部误差指示值从大至小遍历所有三角单元，并进行网格细化；同时，若细化后的网格满足给定的精度要求，执行步骤S6，否则返回步骤S2。

Images