CN108875211B - 一种基于fetd与fdtd耦合的二维模型的探地雷达正演方法 - Google Patents

Abstract

本发明公开了一种基于FETD与FDTD耦合的二维模型的探地雷达正演方法，包括以下步骤：加载二维模型参数；将探地雷达数值模拟区域分为模型区域和扩展区域；采用非结构化三角形单元对模型区域进行剖分，采用规则网格对扩展区域进行剖分；设置测量点，在与测量点对应的位置加载激励源；设置初始时间步之前所有时间的电场与辅助场分量均为零，根据FDTD辅助场时间步进公式更新辅助场，根据FETD主场时间步进公式更新电场，直到最大时间步为止，得到测量点处各个时间步的电场值；根据各个测量点处所有时间步的电场值得到探地雷达正演剖面图。本发明拟合精度高、计算快速。

Description

一种基于FETD与FDTD耦合的二维模型的探地雷达正演方法

技术领域

本发明涉及一种地球物理领域的探地雷达数值模拟方法，特别适合几何特征不规则、物性参数分布复杂模型的快速高效正演。

背景技术

探地雷达数值模拟方法很多，包括有限差分法(Finite Difference TimeDomain，FDTD)^[1-3]、有限单元法^[4-5]、时域伪谱法^[6]、辛分块龙格库塔方法^[7]、时域间断有限元(DGM)^[8]、样条小波有限元等算法^[9]，这些算法各有特色。其中时域有限差分法具有直接时域计算、编程容易、节约存储空间和计算时间等优点而备受青睐(Irving&Knight，2006^[10]；Cassidy&Millington，2009^[11]；冯德山等，2010^[12])，但FDTD受CFL稳定性条件的限制、不能与非结构化网格结合，对不规则脱空、任意形体的开裂、渗水通道等复杂隧道衬砌病害拟合效果不好。而非结构化网格的时域有限元算法(FETD)能够对典型隧道病害实现贴体剖分，具有较高拟合精度。而在模拟区域边界采用CPML边界条件，它能有效压制来自截断边界处的反射波，同时，对传输波、隐失波、低频波等综合吸收效果较好，能极大提升FETD边界的处理效果。但CPML边界的最初提出是针对FDTD算法(Berenger，1996^[13])，公式体系并不适用于FETD。尽管有学者(Masoud&Abdolali,2007^[14]，2009^[15]；Donderici&Teixeira，2008^[16]；Correia&Jin，2006^[17])将CFS-PML边界条件应用于FETD的模拟中；但计算公式过于复杂，不便于编程实现。

因此，有必要设计一种拟合精度高、计算快速的二维模型的探地雷达正演方法。

发明内容

本发明所解决的技术问题是，针对现有技术的不足，提供一种基于FETD与FDTD耦合的二维模型的探地雷达正演方法，主场更新时采用时域有限元法求解，辅助场更新时采用时域有限差分法显式计算的耦合方法，充分利用了FETD算法剖分精度高，而FDTD显式递推计算辅助场计算快速的特点。

本发明所提供的技术方案为：

一种基于FETD与FDTD耦合的二维模型的探地雷达正演方法，包括以下步骤：

步骤1、加载二维模型参数；

步骤2、将探地雷达数值模拟区域分为模型区域和扩展区域；模型区域即二维模型所处区域；扩展区域是根据边界物性参数对模型区域进行扩展得到，包括过渡区域和PML区域，其中过渡区域为设置于PML区域和模型区域之间的一个有限宽度的区域；采用非结构化三角形单元对模型区域进行剖分，采用规则网格对扩展区域进行剖分，得到多个离散网格节点；

步骤3、自左向右设置M个测量点，即M个探地雷达接收天线位置，以及与M个探地雷达接收天线位置一一对应的M个探地雷达发射天线位置；初始化测量点序号m＝1；M为大于或等于1的整数；M个测量点分别设置于M个离散网格节点处；

步骤4、在第m个探地雷达发射天线位置加载激励源，计算与其对应的第m个测量点处各个时间步的电场值；

4.1、设置n为时间步序号，n的取值范围为[1，N]，N为最大时间步；设置第1个时间步及第1个时间步之前所有时间的电场与辅助场分量均为零；初始化时间步序号n＝1；

4.2、根据FDTD辅助场时间步进公式，使用第n个时间步的电场和第n-1/2个时间步的辅助场，更新第n+1/2个时间步的辅助场；其中，第n个时间步的电场，是由探地雷达数值模拟区域中所有离散网格节点第n个时间步的电场z方向分量组成的列向量；第n-1/2个时间步的辅助场，是由探地雷达数值模拟区域中所有离散网格节点第n-1/2个时间步的辅助场v方向分量组成的列向量，v＝x,y；具体地，首先由第n个时间步的电场得到探地雷达数值模拟区域中每一个离散网格节点第n个时间步的电场z方向分量，由第n-1/2个时间步的辅助场得到探地雷达数值模拟区域中每一个离散网格节点第n-1/2个时间步的辅助场v方向分量；再根据FDTD辅助场时间步进公式，分别针对探地雷达数值模拟区域中的每一个离散网格节点，使用其第n个时间步的电场z方向分量和第n-1/2个时间步的辅助场v方向分量，更新其第n+1/2个时间步的辅助场v方向分量；最后由探地雷达数值模拟区域中所有离散网格节点第n+1/2个时间步的辅助场v方向分量组成一个列向量，即第n+1/2个时间步的辅助场；

4.3、根据FETD主场时间步进公式，使用第n-1及第n个时间步的电场和第n+1/2个时间步的辅助场更新第n+1个时间步的电场；

4.4、令n＝n+1，返回步骤4.2，直到达到最大时间步为止；

4.5、分别从各个时间步的电场中提取出第m个测量点处各个时间步的电场值；

步骤5、判断当前测量点(第m个测量点)是否为最后一个测量点，即是否有m＝M，若是则转步骤6，否则令m＝m+1，返回步骤4；

步骤6、根据各个测量点处所有时间步的电场值得到探地雷达正演剖面图。

进一步地，所述步骤2中，所述模型区域采用非结构化Delaunay三角形网格进行剖分。

进一步地，所述步骤2中，在模型区域的剖分过程中，设置h_local<v_loacl/(10f_max)，其中v_local表示局部速度，h_loal表示三角形单元的边长尺度，f_max表示子波的最大频率。

进一步地，所述步骤4.2中，FDTD辅助场时间步进公式为：

其中，对于探地雷达数值模拟区域中任一离散网格节点，

和

表示其第n+1/2个时间步的辅助场v方向分量，v＝x,y，6个变量均为辅助变量，

表示其第n个时间步的电场z方向分量；

为差分系数，其中△t为时间步长；μ为磁导率，单位为H/m，σ_v为PML区域内v方向电导率参数，用于控制波的衰减，为大于零的正实数；α_v值的引入则是为了改善PML对低频分量的吸收特性，κ_v值的引入是为了改善PML对倏逝波的吸收特性，在PML区域内，α_v为大于零的正实数，κ_v为大于或等于1的正实数；在PML区域之外，α_v＝0，κ_v＝1；κ_v′，σ_v′，a_v′分别是κ_v，σ_v，a_v关于变量v的导数。

进一步地，所述步骤4.3中，FETD主场时间步进公式为：

其中，Eⁿ表示第n个时间步的电场，

为第n+1/2个时间步的辅助场；f为激励源源向量，fⁿ表示第n个时间步激励源的值，K＝K_x+K_y，T＝T_x+T_y；t表示时间，采用Newmark-β法进行时间离散，t＝n·△t，△t为时间步长；β为Newmark-β法的参数；

其中，Ω代表整个积分区域，即整个探地雷达数值模拟区域；σ为电导率，单位为S/m，ε为介电常数，单位为F/m；N为形函数，上标T为转置符号。

本发明中FETD+FDTD耦合方法原理及公式体系：

由电磁波传播理论可知(Taflove&Brodwin，1975^[18])，GPR遵循的Maxwell方程的二维TM波的电场波动其时域可表示为：

其中，E为电场强度(V/m)，J为电流密度(A/m²)，E_z为电场强度z方向的分量，J_z为电流密度z方向的分量，t表示时间，σ为电导率(S/m)，μ为磁导率(H/m)，ε为介电常数(F/m)。

考虑在PML区域，将该时域方程应用Fourier公式变换到频率域，并引入复数伸展坐标(Kruzouglu&Mittra,1996^[19])，则该波动方程可表示为：

式(2)中，j为虚数单位，ω是角频率，上标＾都表示变量的傅里叶变换，

为E_z、J_z的傅里叶变换；s_v为坐标拉伸因子，计算方法为：

s_v＝κ_v+σ_v/(a_v+jωε₀) v＝x,y. (3)

式(3)中，σ_v为PML区域内v方向电导率参数，控制波的衰减；κ_v值的引入是为了改善PML对倏逝波的吸收特性，α_v值的引入则是为了改善PML对低频分量的吸收特性，σ_v和α_v为大于零的正实数，κ_v为大于或等于1的正实数；在PML区域之外，α_v＝0，κ_v＝1；由于式(2)中右边两项完全类似，仅仅是坐标方向分量的不同，可用类似的方法求取；以式(2)右边第1式为例，应用复合函数求导公式，并将式(3)代入，可得：

式(4)中κ_x′，σ_x′，a_x′分别是κ_x，σ_x，a_x关于变量x的导数。若直接将式(4)转换至时间域非常复杂，为了方便推导，引入变量：

p＝jωε₀+a_x+σ_x/κ_x (5)

将式(5)代入到式(4)，可以得到：

引入辅助变量e_1x，令

对式(7)两边同时乘以p可以得到：

其中

式(9)两边同时乘以p，可以得到：

其中

通过上式可以得到：

将式(5)带入式(8)(10)和(12)，可以得到如下公式：

对式(13)各方程分别作拉普拉斯逆变换可得：

对于y方向可以采用同样的离散方式。

采用上述推导过程(Ma Y,et al.,2014^[20])，则可推导出角点的CFS-PML区域的时域波动方程：

其中式(16)-(18)为辅助场方程，式(15)为主场方程。式(15)-(18)中的参数设置如下所示：

其中ε_r为相对介电常数，

为真空中的波阻抗，l_v为v方向匹配层到计算区域边界的距离，d为完全匹配层的厚度，m为控制衰减快慢的常数，一般取2-4，R是理论反射系数，通常取R＝1.0×10^-2-10^-6，一般取a_max＝πf₀ε₀。

为了能够实现复杂隧道衬砌病害模型的高精度正演，对式(15)采用非结构化有限元进行空间离散，时间离散采用Newmark-β^[21]隐式算法。对(16)-(18)式采用了时域有限差分法进行时空离散，达到辅助场显式递推的目的，可避免方程组的求解。该算法结合了FETD和FDTD两种算法的优点。应用FDTD法求解(16)-(18)时，常采用4阶差分近似空间偏导数，以x方向差分为例：

式中，n，i和j分别为时间步数，x方向和y方向步数。

注意到辅助场计算中，4阶差分导数的求取需要相邻两个点的场值，为此在PML区域和模型区域(内部计算区域)的非规则网格之间添加两个网格宽度的规则网格，该过渡区域为图7(b)所示的浅灰色区域。

根据(20)-(21)中的4阶差分格式，可将PML区域的辅助微分方程(16)-(18)离散为：

其中

为差分系数；对于探地雷达数值模拟区域中任一离散网格节点，

和

表示其第n+1/2个时间步的辅助场v方向分量，v＝x,y，6个变量均为辅助变量。

主场求解采用FETD离散方程(15)，首先得到其弱解形式：

其中

为试函数，考虑到式(23)的边界项积分为0，当试函数

与形函数N相同时，方程(23)的矩阵表达式可为：

其中，E为电场，是由探地雷达数值模拟区域中所有离散网格节点的电场z方向分量组成的列向量，

分别表示E的一阶时间导数和二阶时间导数，f为源向量，K＝K_x+K_y，T＝T_x+T_y。每一项的具体表达式如下所示：

时间离散采用Newmark-β法，有：

将式(25)代入式(24)，可得：

将式(26)进一步简化可以得

其中，Eⁿ表示第n个时间步的电场，

为第n+1/2个时间步的辅助场，是由所有离散网格节点第n+1/2个时间步的辅助场v方向分量组成的列向量；

当β＝0时，Newmark-β法退化为中心差分法，式(27)为有条件稳定；当β≥0.25时，方程(27)为无条件稳定的，随着β的增大，解的精度会降低，计算中通常取β＝0.25。

由此可以得到FETD+FDTD耦合算法的时间推进步骤为：

(1)根据FDTD辅助场时间步进公式(22)，使用Eⁿ和辅助场

计算

其中n表示时间步；

(2)根据FETD主场时间步进公式(27)，使用E^n-1、Eⁿ和辅助场

计算Eⁿ⁺¹；

(3)回到步骤(1)，直到达到最大时间步(最大时间步由模拟时窗长度以及时间步长确定)为止。

有益效果：

本发明在模拟区域采用Delaunay非结构化网格剖分，应用FETD算法求解，而在模拟区域边界(PML区域)采用FDTD进行辅助场求解，方便卷积完全匹配层(CPML)边界条件的加载与计算，该算法能对几何特征不规则、物性参数分布复杂模型进行快速高效正演，如对任意形状不规则的隧道衬砌病害模型实现精确模拟。与以前的混合算法不同的是应用了非结构化网格剖分，而在模拟区域边界采用CPML边界条件，它能有效压制来自截断边界处的反射波，同时，对传输波、隐失波、低频波等综合吸收效果较好，能极大提升FETD边界的处理效果。本发明拟合精度高、计算快速。

附图说明

图1为隧道衬砌不密实模型图及雷达正演剖面图；图1(a)为模型剖分示意图，图1(b为雷达正演剖面图；

图2为隧道衬砌渗漏水模型图及雷达正演剖面图；图2(a)为模型图，图2(b)为雷达正演剖面图；

图3为无水空洞、含水空洞短时傅里叶分析；3(a)和图3(b)分别为无水空洞与含水脱空的上界面反射信号，图3(c)和图3(d)分别为无水空洞和含水脱空频谱图；

图4为隧道衬砌渗漏水模型图及雷达正演剖面图；图4(a)为模型图，图4(b)为雷达正演剖面图；

图5为无水裂隙、含水裂隙频谱对比；图5(a)和图5(b)分别为无水通道和含水通道上下界面反射波，图5(c)和图5(d)分别为无水通道和含水通道频谱图；

图6为钢筋网干扰下的衬砌裂缝及脱空模型图及雷达正演剖面图；图6(a)为模型图剖分示意图，图6(b)为雷达正演剖面图；

图7为雷达正演剖面图；图7(a)为模型示意图，图7(b)为耦合算法网格生成示意图；

图8为3种方法计算的单道雷达数据；

图9为本发明流程图。

具体实施方式

下面以本发明应用于隧道衬砌典型病害模型为例进行说明。

隧道在经历了多年运营后，由于围岩压力变化，衬砌材料老化、强度降低等因素，导致衬砌发生开裂，衬砌掉块、剥离、脱空、渗漏水等现象，严重威胁隧道运营安全。应用本发明对衬砌不密实模型、隧道渗漏水模型、衬砌脱空模型、钢筋网干扰下的衬砌裂缝及脱空模型进行正演计算，以指导实测衬砌病害模型的探地雷达资料解释。

1.1隧道不密实模型图

在隧道施工过程中，因地质条件复杂、施工环境恶劣、施工工艺等问题，导致衬砌中形成不密实部位。衬砌中的不密实部位会降低衬砌强度，容易造成衬砌结构破坏。衬砌不密实模型如图1(a)所示的。模型的长与宽分别为1.0m×0.5m，采用非结构化网格对该模型内部区域进行离散，离散时的内部的三角单元数为147019个，节点个数为70890个；地表上方设置0.1m的空气层。介质的电性参数如表1所示。激励源采用900MHz的零相位Ricker子波，时间步长为0.01ns，时窗长度为12ns。收发天线距地表0.01m，收发距为0.1m，收发天线位置如图1(a)所示，“·”表示激发天线位置，“×”表示接收天线位置，两者自左向右同步移动，采样间隔为0.005m，共计录181道雷达数据。

图1(a)中的上层介质为混凝土衬砌，下层介质为围岩，黑色点状区域为不密实部位，黑点表示小空洞，内部充填空气。衬砌不密实模型的FETD+FDTD耦合算法雷达正演剖面如图1(b)所示，区域A出现大量不规则抛物线状反射信号叠加，是由于不密实部位存在大量小空洞，电磁波信号发生多次反射叠加形成的，不密实部位的强反射同相轴发生错乱、不连续甚至断开；区域B为雷达探测到的围岩分界面，信号强度较弱。

1.2隧道衬砌脱空模型图

隧道在施工时超挖未回填密实导致隧道围岩与衬砌之间存在脱空区。衬砌中的空洞(脱空)会降低衬砌强度，使衬砌结构受力不均，应力集中，容易造成衬砌结构破坏。为了更好地认识衬砌脱空模型的雷达剖面特征，建立图2(a)所示的衬砌空洞(脱空)模型。模型的长与宽分别为1.0m×0.5m，离散时的内部非结构化三角单元数为210074个，节点个数为105425个，激励源采用900MHz的零相位Ricker子波，时间步长为0.02ns，时窗长度为20ns。天线位置与隧道衬砌不密实模型相同。其中左侧区域为无水空洞，右下方区域为富水衬砌脱空区，相关介质的电性参数如表1所示。

表1介质电性参数

Tab.1 Electrical parameters of media

衬砌脱空模型的雷达正演剖面如图2(b)所示，图中可见，电磁波遇到空洞后产生强反射，区域A电磁波强反射同相轴为一大型抛物线，比较容易识别为一个大型空洞；区域B为雷达探测到的围岩分界面，信号强度较弱；区域C电磁波强反射同相轴为一大型抛物线，结合围岩分界面位置和深度信息，可以推断为衬砌含水脱空区，该衬砌脱空与区域的雷达图像上表现为下界面的反射波更强，且相位发生了反转。区域D为含水脱空的下界面反射波。

为了区分雷达剖面图2(b)中含水脱空与无水空洞，方便对异常进行定性，对无水空洞和含水脱空分别进行短时傅里叶分析，得到频谱对比图3(c)-(d)，分析图3(a)-(b)可知，空洞与含水脱空的上界面反射信号相位相反，上界面反射波主频与直达波主频基本相同。图3(c)中10ns处的反射异常系多反射波叠加，空洞下界面反射波无法分辨，对应主频在1200MHz左右；图3(d)中含水脱空界面反射波(12ns)的主频在位于600MHz左右，含水通道的反射波主频明显降低。

1.3隧道渗漏水模型图

衬砌周围岩石破碎，存在导水通道，水通过衬砌中的缝隙渗漏出，衬砌渗漏水能够使混凝土衬砌风化剥蚀、腐蚀破坏衬砌中的钢筋，对衬砌结构造成破坏；在冬季寒冷地区，渗漏水会造成衬砌挂冰，衬砌中的空洞、裂缝成为储水空间，使衬砌由于冻胀而发生破坏。为了更好地认识衬砌渗漏水模型的雷达剖面特征，建立图4所示的衬砌渗漏水模型。模型的长与宽分别为5.0m×2.5m，采用非结构化网格对该模型内部区域进行离散，离散时的内部三角单元数为613293个，节点个数为307451个，介质的电性参数如表1所示。模拟的时间步长为0.05ns，时窗长度为70ns。激励源采用400MHz的零相位Ricker子波，收发天线距地表0.025m，收发距为0.1m，自左向右同步移动，采样间隔为0.025m，共记录196道雷达数据。

图4(a)隧道衬砌漏水模型图中的灰色介质为混凝土衬砌，黄色介质为围岩，中间黑区域为围岩破碎带，右侧黑分叉区域为围岩裂隙、破碎带和裂隙构成的导水通道，内部充水，左侧白色区域为围岩裂隙、破碎带构成的无水裂隙。相关介质的电性参数如表1所示。

衬砌渗漏水模型的雷达正演剖面如图4(b)所示，图4(b)中标注的区域A为衬砌与围岩的分界面，该起伏不平界面的形状能得到大致体现；区域B为中间破碎带上界面产生的反射波，由于水与围岩的介电常数之差较大，电磁波遇到含水通道后产生很强的反射振幅，电磁波反射信号同相轴的形态能大致反映围岩破碎带、裂缝的形态，但倾角变得更平缓；区域C与区域D分别为含水裂隙与无水裂隙的端点处的绕射波及裂隙通道的反射波，显然由于水的介电常数与围岩的差异较大，则含水裂隙端点处的绕射波及裂隙通道的反射波能量更强，对比两者的电磁波反射信号相位发生了180°反转；区域E为电磁波遇到围岩破碎带时产生比较杂乱的强反射信号。

为了进一步对不同裂隙进行定性，对无水通道和含水通道分别进行短时傅里叶分析，得到图5(c)-(d)所示的频谱对比图，图中可见，两者直达波主频基本相同，都在400MHz左右，图5(a)与图5(b)中无水通道上下界面反射波(40ns处)基本叠加在一起，其对应位置主频在600MHz左右；图5(b)与5(d)中含水通道上界面反射波的主频在400MHz左右(20ns)，含水通道的下界面反射波(36ns)主频明显降低，这是由于水对高频电磁波的衰减较大，高频反射波能量较低之故。

1.4钢筋网干扰下的衬砌裂缝及脱空模型图

衬砌病害常常不是孤立存在的，多与钢筋网共存，给衬砌病害雷达探测造成强干扰，提高了雷达探测与解释的难度。钢筋网干扰下的衬砌裂缝及脱空模型如图6(a)所示。模型图中上层介质为混凝土衬砌，下层介质为围岩，灰色的小圆圈为钢筋网，左边区域为衬砌脱空，右侧白色交叉线条为裂隙，模型的长与宽分别为1.0m×0.5m。离散时的内部非结构化三角单元数为959583个，节点个数为48373个，激励源采用900MHz的零相位Ricker子波，时间步长为0.01ns，时窗长度为12ns。天线位置与隧道衬砌不密实模型相同，介质的电性参数如表1所示。

钢筋网干扰下的衬砌裂缝及脱空模型的雷达正演剖面如图6(b)所示。图中可见，由于钢筋为金属材质，对电磁波信号表现为全反射，钢筋产生了很强的抛物线状绕射信号，由于钢筋的强反射干扰导致钢筋下层区域探测效果较差，但仍可以看到一处强反射信号(如区域A所示)，可以判断为空洞，其位置及尺寸与模型较为一致；区域B可以看到反射信号同相轴连续一段距离，推断为衬砌裂缝，且同相轴形态与衬砌裂缝形态较为一致；由于钢筋产生的电磁波干扰信号很强，导致难以识别围岩分界面。

本申请与以往技术的效果对比

采用如图7(a)所示简单的圆形水洞模型对混合算法的网格生成方式进行介绍。模型区域大小为0.5m×0.6m，模型上方设置0.1m的空气层，背景介质为混凝土，混凝土介质中有一个埋深为0.15m，半径为0.08m圆形水洞，介质的电性参数如表1所示。FETD+FDTD混合算法的网格剖分示意图如图7(b)所示，将模拟区域分为模型区域(A)和扩展区域(包括过渡区域B+PML区域C)，模型区域采用非结构化三角形单元进行剖分，扩展区域采用规则网格进行剖分。对于模型区域而言，需要模拟任意复杂模型，而非结构化Delaunay三角形网格剖分质量好，单元密度易控制，能够实现物性参数分布复杂、几何特征不规则的精确拟合，因此特别适用于模型区域的剖分。在剖分过程为了避免数值频散，将网格的空间步长与网格单元的局部速度关联起来，其满足的条件为h_local<v_loacl/(10f_max)，其中v_local表示局部速度，h_loal表示三角形单元的边长尺度，f_max表示子波的最大频率。

扩展区域是根据边界物性参数对模型区域进行扩展得到的，扩展的宽度与模拟的区域大小以及天线的频率相关，通常取22*h_backroun，其中h_backroun表示背景介质速度计算得到的网格单元尺度。在扩展区域采用结构化规则网格剖分，方便FDTD进行辅助场的计算。为了方便模型区域的非结构化网格与扩展区域结构化网格的过渡，在规则网格与非规则网格之间设计一个有限宽度的过渡区域，该过渡区域为2个(由差分阶次确定)有限差分网格间距，见图7(b)中浅灰色区域(B区域)。编程过程中需要建立从规则网格的有限元总体编码到有限差分局部编码的映射，以方便辅助场的显式递推计算。

整个模拟区域需要进行时域有限元方法的计算，扩展区域的形函数可以选取四边形双线性插值形函数或者采用与模型区域相同的三角形线性插值形函数。时间域离散采用Newmark-β隐式算法。该隐式算法突破CFL稳定性条件的限制，可以根据采样需求选取合适的时间步长。

为了说明本发明中FETD+FDTD算法的计算效率及精度，将该算法与粗网格FDTD、细网格FDTD算法进行对比。在Tx处(0.2,-0.04)加载主频为900MHz的Ricker子波，并在Rx处(0.3,-0.04)接收信号，模拟时窗为20ns。其中，细网格FDTD的空间步长为0.001m和时间步长为0.002ns满足频散条件min(△x,△y)<v_min/(10f_max)和稳定性条件

模拟区域剖分为501×601，模拟时间297.52s；粗网格模拟步长根据背景介质的速度参数进行设置，其空间步长为0.004m，时间步长为0.009ns，模拟区域剖分126×151，模拟时间4.02s。混合算法的网格空间步长根据网格单元的局部速度参数进行设置，内部剖分133748个单元，67150个节点，时间步长取0.02ns，模拟时间为116.52s。计算环境为Intel(R)Core^TMi7-6700K CPU，@4.00GHz，16.0GB的内存物理地址扩展，Window10操作系统的Dell XPS8900台式机。

三种方法模拟得到的单道雷达记录如图8所示。图8中黑色的曲线代表细网格方法，点划线代表粗网格方法，虚线代表FETD+FDTD耦合算法。3种算法结果基本吻合，一致性较好，在直达波的位置，粗网格的场值较另外两种方法稍大，2.6ns与12.3ns附近的正负波峰很好地对应了圆形异常体的上下界面；图8中这两个位置处虚线(FETD+FDTD耦合算法)与黑色的曲线(细网格)吻合得更好，点划线(粗网格)与另外两条曲线稍有偏差。如果认为细网格的计算结果与解析解接近，说明FETD+FDTD耦合算法在较短运行时间的前提下，计算结果具有较高的精度，在计算效率与计算精度方面取得了较好的平衡。

参考文献：

[1]Zhang W Y,Hao T,Chang Y,Zhao Y H.Time-frequency analysis ofenhanced GPR detection of RF tagged buried plastic pipes.NDT&E International,2017,92(6):88-96.

[2]Yee K S.Numerical solution of initial boundary value problemsinvolving Maxwell’s equations in isotropic media.IEEE Transactions onAntennas and Propagation,1996,14(3):302-307.

[3]Chen H W,Huang T M.Finite-difference time-domain simulation of GPRdata[J].Journal of Applied Geophysics.1998,40(1-3)：139-163.

[4]Feng D S,Chen C S.Finite element method GPR forward simulationbased on mixed boundary condition.Chinese J.Geophys.(in Chinese),2012,55(11):3774-3785.

[5]Feng D S,Wang X.Convolution perfectly matched layer for theFinite-Element Time-Domain method modeling of Ground Penetrating Radar,Chinese J.Geophys.(in Chinese),2017,60(1):413-423.

[6]Liu Q H,Fan G-X.Simulations of GPR in dispersive media using afrequency-dependent PSTD algorithm.IEEE Transactions on Geoscience and RemoteSensing,1999,37(5):2317-2324.

[7]Fang H Y,Lin G.Symplectic partitioned Runge–Kutta methods for two-dimensional numerical model of ground penetrating radar.Computers&Geosciences,2012,49:323–329.

[8]Lu T，Cai W，Zhang P W.Discontinuous Galerkin Time-Domain Method forGPR Simulation in Dispersive Media[J].IEEE Transactions on Geoscience andRemote Sensing，2005，43(1)：72-80.

[9]Feng D S,Wang X.The GPR simulation of bi-phase random concretemedium based on finite element method of B-spline wavelet on the interval,Chinese Journal of Geophysics,2016,59(8):3098-3109.

[10]Irving J,Knight R..Numerical modeling of ground-penetrating radarin 2-D using MATLAB.Computers&Geosciences,2006,32(9):1247-1258.

[11]Cassidy N J,Millington T M.The application of finite-differencetime-domain modelling for the assessment of GPR in magnetically lossymaterials.Journal of Applied Geophysics,2009,67(4):296-308.

[12]Feng D S,Chen C S.Finite element method GPR forward simulationbased on mixed boundary condition.Chinese J.Geophys.(in Chinese),2012,55(11):3774-3785.

[13]Berenger J P.Three-dimensional perfectly matched layer for theabsorption of electromagnetic waves[J].Journal of Computational Physics,1996,127(2):363-379.

[14]Masoud M,Abdolali A.Optimization of the perfectly matched layerfor the finite-element time-domain method.IEEE Microwave and wirelesscomponents letters.2007,17(1):10-12.

[15]Masoud M,Abdolali A.Complex frequency shifted-perfectly matchedlayer for the finite-element time-domain method.International Journal ofElectronics and communications.

.2009.,63:72-77.

[16]Donderici B,Teixeira F L.Conformal perfectly matched layer forthe mixed finite element time-domain method.IEEE Transactions on Antennas andPropagation,2008,56(4):1017-1026.

[17]Correia D,Jin J M.2006.Performance of regular PML,CFS-PML,andsecond-order PML for waveguide problems.Microwave and Optical TechnologyLetters,48(10):2121-2126.

[18]Taflove A,Brodwin M E.Numerieal solution of steady-state EMscattering Problems using the time-dependent Maxwell’s equations.IEEETrans.MicrowaveTheory Tech,1975,23:623-630.

[19]Kuzuoglu M,Mittra R.Frequency dependence of the constitutiveparameters of causal perfectly matched anisotropic absorbers.IEEE Microwaveand Guided Wave Letters,1996,6(12):447-449.

[20]Ma Y,Yu J,Wang Y.An efficient complex-frequency shifted-perfectlymatched layer for second-order acoustic wave equation[J].InternationalJournal for Numerical Methods in Engineering,2014,97(2):130-148.

[21]Newmark N M.A Method of Computation for Structural Dynamics.ProcAsce,1959,85(1):67-94.

Claims (3)

1.一种基于FETD与FDTD耦合的二维模型的探地雷达正演方法，其特征在于，包括以下步骤：

步骤1、加载二维模型参数；

步骤2、将探地雷达数值模拟区域分为模型区域和扩展区域；模型区域即二维模型所处区域；扩展区域是根据边界物性参数对模型区域进行扩展得到，包括过渡区域和PML区域，其中过渡区域为设置于PML区域和模型区域之间的一个有限宽度的区域；采用非结构化三角形单元对模型区域进行剖分，采用规则网格对扩展区域进行剖分，得到多个离散网格节点；

步骤3、自左向右设置M个测量点，即M个探地雷达接收天线位置，以及与M个探地雷达接收天线位置一一对应的M个探地雷达发射天线位置；初始化测量点序号m＝1；M为大于或等于1的整数；M个测量点分别设置于M个离散网格节点处；

步骤4、在第m个探地雷达发射天线位置加载激励源，计算与其对应的第m个测量点处各个时间步的电场值；

4.1、设置n为时间步序号，n的取值范围为[1，N]，N为最大时间步；设置第1个时间步及第1个时间步之前所有时间的电场与辅助场分量均为零；初始化时间步序号n＝1；

4.2、根据FDTD辅助场时间步进公式，使用第n个时间步的电场和第n-1/2个时间步的辅助场，更新第n+1/2个时间步的辅助场；其中，第n个时间步的电场，是由探地雷达数值模拟区域中所有离散网格节点第n个时间步的电场z方向分量组成的列向量；第n-1/2个时间步的辅助场，是由探地雷达数值模拟区域中所有离散网格节点第n-1/2个时间步的辅助场v方向分量组成的列向量，v＝x,y；其中，FDTD辅助场时间步进公式为：

其中，对于探地雷达数值模拟区域中任一离散网格节点，

和

表示其第n+1/2个时间步的辅助场v方向分量，v＝x,y，6个变量均为辅助变量，

表示其第n个时间步的电场z方向分量；

为差分系数，其中△t为时间步长；μ为磁导率，单位为H/m，σ_v为PML区域内v方向电导率参数，用于控制波的衰减，为大于零的正实数；α_v值的引入则是为了改善PML对低频分量的吸收特性，κ_v值的引入是为了改善PML对倏逝波的吸收特性，在PML区域内，α_v为大于零的正实数，κ_v为大于或等于1的正实数；在PML区域之外，α_v＝0，κ_v＝1；κ′_v，σ′_v，a′_v分别是κ_v，σ_v，a_v关于变量v的导数；

4.3、根据FETD主场时间步进公式，使用第n-1及第n个时间步的电场和第n+1/2个时间步的辅助场更新第n+1个时间步的电场；其中，FETD主场时间步进公式为：

其中，Eⁿ表示第n个时间步的电场，

为第n+1/2个时间步的辅助场；f为激励源源向量，fⁿ表示第n个时间步激励源的值，K＝K_x+K_y，T＝T_x+T_y；t表示时间，采用Newmark-β法进行时间离散，t＝n·△t，△t为时间步长；β为Newmark-β法的参数；

其中，Ω代表整个积分区域，即整个探地雷达数值模拟区域；σ为电导率，单位为S/m，ε为介电常数，单位为F/m；N为形函数，上标T为转置符号；

4.4、令n＝n+1，返回步骤4.2，直到达到最大时间步为止；

4.5、分别从各个时间步的电场中提取出第m个测量点处各个时间步的电场值；

步骤5、判断当前测量点是否为最后一个测量点，即是否有m＝M，若是则转步骤6，否则令m＝m+1，返回步骤4；

步骤6、根据各个测量点处所有时间步的电场值得到探地雷达正演剖面图。

2.根据权利要求1所述的基于FETD与FDTD耦合的二维模型的探地雷达正演方法，其特征在于，所述步骤2中，所述模型区域采用非结构化Delaunay三角形网格进行剖分。

3.根据权利要求2所述的基于FETD与FDTD耦合的二维模型的探地雷达正演方法，其特征在于，所述步骤2中，在模型区域的剖分过程中，设置h_local<v_loacl/(10f_max)，其中v_local表示局部速度，h_loal表示三角形单元的边长尺度，f_max表示子波的最大频率。

Images