0 引言

探地雷达(Ground Penetrating Radar, GPR)作为一种采用高频脉冲电磁波来探测地下地质体结构及物性参数的重要浅部地球物理方法(刘澜波和钱荣毅, 2015)，具有无损性、高分辨率、高效率等特点，在工程勘查与环境检测等领域中得到广泛应用(González-Drigo et al., 2008; Chang et al., 2009; Solla et al., 2011; Alani et al., 2013).实验研究表明：GPR的探测对象如土壤(Atteia and Hussein, 2010)、混泥土(Bikowski et al., 2012)、岩石(Deparis and Garambois, 2008)等在GPR工作频率范围内表现出强烈的频散特性，其相对介电常数不是固定常数，而是关于频率的复函数.电磁波在上述高损耗频散介质中传播时能量衰减很快，并常伴有色散与畸变，严重影响了GPR的探测效果和资料解释.因此，深入研究频散介质的GPR正演模拟算法，通过模拟分析介质频散特性对电磁波传播的影响规律，探讨GPR在频散介质中的探测能力，有助于提高实测资料的解释精度(Mangel et al., 2015).

目前描述频散介质中复相对介电常数随频率变化关系的模型主要有德拜(Debye)模型(王洪华和戴前伟, 2014)、洛伦兹(Lorentz)模型(Liu and Fan, 1999)、德鲁(Drude)(Mustafa et al., 2014)和柯尔-柯尔(Cole-Cole)模型(Giannakis et al., 2016).前三种频散介质的复相对介电常数是jω的整数次幂函数，利用傅里叶变换可以将频率域电磁波方程变为相应的时间域整数阶导数方程，借助于递归卷积方法(Giannakis and Giannopoulos, 2014)、辅助方程法或定义位移算子法(Lu et al., 2005)，可快速高精度地实现频散介质的GPR正演模拟.Atteia和Hussein (2010)，Teixeira (2008)，曾昭发和刘四新(2007)、冯德山(2014)、李展辉(2009)等都利用时域有限差分法(Finite Difference Time Domain，FDTD)开发了基于Debye频散介质的GPR正演模拟算法；Liu和Fan(1999)提出并实现了计算Debye频散介质的时域伪谱法(Pseudo spectral Time Domain, PSTD)；Shibayama等(2009)基于梯形递归卷积技术，提出了一种计算Drude-Lorentz频散模型的FDTD正演算法；Shahmansouri和Rashidian(2012)借助于辅助微分方程法，提出了一种计算Drude-Lorentz频散模型的FDTD正演算法.然而实验研究表明，相比Debye，Lorentz，Drude等频散介质模型，Cole-Cole模型更适合描述浅部实际地下介质的频散特性(Giannakis et al., 2012; Chung et al., 2013).不同于前三种频散介质模型，Cole-Cole模型的复相对介电常数是关于jω的分数次幂函数，将其转换到时间域会出现时间分数阶导数，应用数值模拟方法开展Cole-Cole频散介质计算时，现有的复介电常数离散求解技术如递归卷积法、辅助方程法和定义位移算子法都将电磁波方程表示成时间分数阶微分方程，传统的时间域差分法如：中心差分、隐式差分和Nermark差分法等都只能离散时间域整数阶导数，而不能离散时间分数阶导数.因此，如何在时间域离散分数阶导数是开展Cole-Cole频散介质GPR正演计算急需解决的关键问题.为此，国内外学者先后提出了Debye近似法(Giannakis et al., 2012)和Pade近似法两种计算时间域分数阶导数策略(Rekanos, 2012).其中Debye近似法，利用不同阶次的Debye模型的线性组合构建Cole-Cole模型，需要求解非线性优化问题，不仅耗时，而且容易陷入局部最优；而Rekanos和Papadopoulos(2010)提出的Pade近似法是一种对幂级函数进行有理函数逼近的高效算法，具有比截断的泰勒级数逼近更精确的优点.将其应用于Cole-Cole频散介质GPR电磁波方程离散时，可将含频率域分数阶导数的复介电常数展开为一组有理函数求和表达式，进而可将含有时间域分数阶导数Cole-Cole频散介质电磁波方程导出为一组整数阶时间域辅助微分方程，从而方便利用现有的差分法实现电磁波在时间域迭代，具有理论简单、精度高的优点.目前，该算法被广泛应用于求解Cole-Cole频散介质的电磁波传播问题.

相比FDTD和PSTD法，时域有限元法(Finite element time domain，FETD)采用非结构化网格离散复杂介质与结构的几何模型，并能自然满足自由边界条件，具有计算精度高的优点.因此，一些学者将FETD引入到GPR正演模拟中：底青云和王妙月(1999)推导了含衰减项的电磁波有限元方程，实现了复杂地质体的FETD正演模拟；戴前伟等(2012)采用二次插值基函数的FETD开展了GPR有限元正演，提高了正演精度；戴前伟等(2012)采用三角网格的FETD实现了复杂结构的GPR正演；冯德山等(2012)结合透射边界条件和Sarma边界条件各自的优势，提出了一种混合吸收边界条件，在一定程度上改善了GPR有限元正演的精度；冯德山和王珣等(2016)基于四边形网格的FETD实现了GPR正演；王洪华和戴前伟(2014)等开发了Debye频散介质的FETD正演模拟算法；Lu等(2005)利用离散时域Galerkin有限元法实现了Debye频散介质的GPR正演.上述文献主要集中于非频散介质、Debye频散介质的GPR正演模拟，将FETD应用到Cole-Cole频散介质的GPR正演模拟中还未见报道.此外，在人工边界处理上，由于GPR有限元正演模拟基于二阶电磁波方程，基于一阶电磁波方程提出的各种完全匹配层吸收边界条件(Perfectly matched layers boundary condition, PML)还主要集中于FDTD、PSTD算法中，难以直接应用于GPR有限元正演中.因此，传统的Sarma边界条件(Wang et al., 2014)、傍轴近似边界条件(底青云和王妙月, 1999)及它们的混合边界条件(冯德山等, 2012)被广泛应用于GPR有限元正演，而基于PML边界条件的FETD在声波、弹性波、面波(Collino和Tsogka, 2001; Komatitsch和Tromp, 2003; Drossaert和Giannopoulos, 2007; 李义丰等, 2010; 刘有山等, 2013; Zhao和Shi, 2013)正演模拟中得到成功应用，模拟结果论证了在人工截断边界处理上，以变分形式将PML边界条件引入到有限元方程组中可取得良好的吸收效果.

为实现Cole-Cole频散介质的GPR有限元高精度模拟，本文基于Pade逼近算法将含有时间域分数阶导数的电磁波方程变为一组整数阶辅助微分方程，开发了一种适用于Cole-Cole频散介质的GPR有限元正演模拟算法.采用吸收效果较好的PML边界条件处理来自截断边界处的超强反射波；应用集中质量矩阵技术以提高计算效率；时间域离散采用保能量的Newmark差分算法以在保证计算精度的同时提高计算效率.数值模拟结果表明，本文构建的GPR有限元正演算法可用于复杂Cole-Cole频散介质结构模拟，且具有较高的计算精度.

1 Pade逼近算法求解Cole-Cole频散介质中的电磁波方程组

根据电磁波理论，GPR在频散介质中的传播规律满足Maxwell方程组，考虑二维地电模型的走向为z轴，则TM模式的电磁波方程组为

(1a)

(1b)

(1c)

式中H_x，H_y分别为x, y方向的磁场强度(A·m^-1)；E_z为z方向的电场强度(V·m^-1)；ε₀为真空的介电常数(F·m^-1)；ε_r为介质的相对介电常数(F·m^-1)，它是关于频率ω的函数；μ₀表示真空的磁导率(H·m^-1).

考虑介质电导率的影响，在GPR工作频带范围内，Cole-Cole频散介质模型的复相对介电常数可表示为

(2)

式中σ表示介质的电导率(S·m^-1)，ε_r⁰和ε_r^∞分别为频率趋近于0和无穷大的相对介电常数，τ为驰豫时间，β(0≤β≤1)为衰减因子，当β=1时，式(2)退化为Debye频散介质.

将式(1b)、(1c)和(2)代入到式(1a)，整理可得

(3)

式(3)中辅助变量J的表达式为

(4)

对式(4)进行反傅里叶变换，并将其转变到时间域可得

(5)

式(5)中，Δε_r=ε_r⁰－ε_r^∞，由于存在时间分数阶导数项，利用中心差分法在时间域离散求解非常困难.为此，采用Pade近似方法，将(jω)^β用有理函数离散化，令x=j(ω/ω₀)，式(4)可写为

(6)

式(6)中：ω₀为GPR脉冲激励源的中心频率.由于x=0是函数f(x)=x^β的导数的奇点，函数f(x)在x=0处的泰勒展开不存在.因此，在半径为|x－w|的区域内由于包含奇点x=0，f(x)的泰勒展开不存在.此外，f(x)在其他点处的泰勒展开式中包含分数阶导数，利用时间域差分法计算非常困难.为此，采用Pade逼近对f(x)在w处进行有理函数展开，为了计算简便，取w=1，其Pade逼近表达式为(Rekanos, 2012)

(7)

式(7)中M, N分别为分母和分子多项式的阶数，q₀=1，a_{k, k－l}=(k！)/(l!(k－l)!)(－1)^k－l，有理函数的系数p_n和q_m的计算方法见参考文献(Rekanos and Papadopoulos, 2010).

将式(7)代入式(6)中可得

(8)

式(8)中系数A_m和B_n的计算公式为

(9)

对比式(2)、(4)和(8)可知，利用Pade逼近算法计算的复相对介电常数为

(10)

将式(8)进行反傅里叶变换，并转换到时间域可得

(11)

式(3)和式(11)即为电磁波在Cole-Cole频散介质中传播满足的电磁波方程.

2 Cole-Cole频散介质GPR有限元方程及PML边界条件

根据PML边界条件理论(Chew and Liu, 1996)，PML边界条件表述为在频率域引入复数伸展坐标系.为构建电磁波在Cole-Cole频散介质中传播满足的方程式(3)和(11)的PML边界条件，建立复数伸展坐标系为

(12)

式(12)中，为复数空间坐标，j为虚数单位，d_p为吸收边界衰减函数，是随坐标p衰减的实函数，d_p的经验值为d_p=－(n+1)log(R)V_p/(2δ)·(d/δ)ⁿ，δ为PML层的厚度；V_p为电磁波传播速度；d为吸收层到计算区域边界的距离；n为控制衰减系数快慢的常数，一般取2~3；R为理论反射系数，一般取1.0×10^-2~10^-6.

由式(12)可知，函数关于p和的一阶偏导数有如下关系：

(13)

将式(13)代入式(3)可得

(14)

令E_z=E_{1, z}+E_{2, z}，式(14)可拆分成如下方程组

(15)

对式(15)进行反傅里叶变换，可得到加载了PML边界条件的时间域方程为

(16)

式(16)中，P_x, P_y为引入的中间变量，d′_x=∂d_x(x)/dx, d′_y=∂d_y(y)/dy.在计算区域内d_x(x)=0，d_y(y)=0，在吸收层内d_x(x)≠0，d_y(y)≠0用以消除人工边界处的虚假反射波.加载PML边界条件的方程组为

(17)

方程(17)相应的弱形式为

(18)

式(18)中，Ω, Γ分别是计算区域及边界；n=(n_x, n_y)^T为计算区域边界的外法线向量；φ=(φ_x, φ_y)^T为试函数.在PML边界上利用电磁场为0的dirichlet边界条件，上边界为自由边界条件.

根据Galerkin有限元理论(王洪华和戴前伟, 2014)，式(18)形成的有限元方程组为

(19)

式(19)只在PML区域进行计算，分裂变量E_{1, z}, E_{2, z}只在PML区域内进行存储，相比于在整个计算区域内存储，节约了存储量.式(19)中系数矩阵的计算表达式为

(20)

(21)

(22)

式(20)、(21)、(22)中，N是计算区域内的形函数，N_x和N_y分别是计算区域内形函数关于x和y的导数.

3 Newmark差分法求解GPR有限元方程

频散介质GPR有限方程组式(19)中，前四式的一般形式为

(23)

其中，M为质量矩阵，C为阻尼矩阵，K₁, K₂为刚度矩阵，S为源向量；和分别表示E_z关于时间t的一阶、二阶导数.上述有限元形成的矩阵是一个大型的稀疏矩阵，全部存储需要很大的内存空间.本文对刚度矩阵采用压缩存储行(CSR)格式，只需存储其非零元素，可极大减少所需存储空间(刘有山等, 2013).

Newmark差分算法(Newmark, 1959)最初由Newmark提出，是一种时间积分算法，具有保能量的特性(当δ=0, γ=1/2时).该算法将点更新到点，式(23)的Newmark差分算法的实现过程如下：

(24)

式中，n为时间步数，a为加速度，当γ≤1/2是能量消耗的；γ＞1/2是能量增长的；当且仅当δ=0, γ=1/2是保能量的. 表示电场对时间t的一次积分，利用中心差分法，可推导其差分格式为

(25)

式(19)第5式涉及到关于n阶时间导数的计算，根据差分原理，矢量F的n阶中心差分公式为

(26)

其中，l_n=[n/2]表示对n/2取整, Q_mⁿ表示n阶中心差分式中的第m个中心差分参数.当n=0, 1, 2, 3, 4时的中心差分参数如表 1所示.

将式(24)、(25)、(26)代入式(19)进行时间域离散和迭代，则可计算不同时刻的电磁波场值，从而实现Cole-Cole频散介质的GPR有限元正演.

4 数值算例 4.1 Pade逼近算法的精度及阶数选择

在上述理论基础上，编制了相应的Cole-Cole频散介质的复介电常数Pade逼近程序，为验证Pade逼近算法的正确性并选择合适的阶数，设计了一种Cole-Cole频散介质Ⅰ，它的相对介电常数为ε_r^∞=3, ε_r⁰=6，电导率为σ=0.0005 s·m^-1，弛豫时间τ=100 ps，衰减因子β=0.50.利用Pade逼近算法进行计算时的频率范围为10 MHz~10 GHz，分母多项式阶数分别取为M=1, 2, …, 10，相应的分子多项式分别取为N=M，N=M-1，M-2.定义Pade逼近的相对均方根误差为

(27)

式(27)中，, ε_r分别为Pade逼近算法计算的复相对介电常数值和真实相对介电常数值，f_L, f_H分别为计算的上限频率值和下限频率值.

图 1为利用Pade逼近算法对Cole-Cole频散介质Ⅰ的相对复介电常数进行计算的相对均方根误差随分母多项式阶数M的变化曲线，由图可知：随着分母多项式阶数M的增加，Pade逼近的相对均方根误差逐渐减小；当分母多项式阶数M相同时，相比分子多项式阶数N=M-1和N=M-2，当N=M时的相对均方根误差最小.当分母多项式阶数M增加时，由式(19)可知，FETD正演迭代公式增多，计算存储量和耗时也相应增加.因此，我们建议在多项式阶数和精度之间进行折衷选择.当N=M=4时，Cole-Cole频散介质Ⅰ的复相对介电常数的相对均方根误差约为0.0098%，计算精度能够满足一般的工程应用需求，后文的正演实例多项式阶数默认取N=M=4.当N=M=4时，利用Pade逼近算法计算的Cole-Cole频散介质Ⅰ的复相对介电常数与解析解对比，如图 2所示.由图可知：在GPR工作频段(10 MHz~10 GHz)范围内，频散介质的相对介电常数不是一个常数，而是随频率变化的函数；复相对介电常数的实部和虚部曲线与解析解拟合得非常好，误差分别为0.003%和0.002%.计算结果表明，Pade逼近算法能够高精度地计算Cole-Cole频散介质的复相对介电常数，在兼顾计算精度和效率的前提下，选择多项式阶数N=M=4是合理的.

4.2 正演模拟算例 4.2.1 算例1

为验证基于Pade逼近的Cole-Cole频散介质GPR有限元正演算法的正确性和有效性，分析GPR在频散介质中的传播特征及探测能力，建立了一个5.0 m×4.0 m的均匀Cole-Cole频散介质模型，模型的介电参数与Cole-Cole频散介质Ⅰ相同.Cole-Cole频散介质GPR有限元正演时的计算区域和吸收层被剖分成19981个非结构化三角形网格，网格内采用线性插值基函数，计算区域外的PML吸收边界层厚度为0.2 m.GPR脉冲激励源为雷克子波，其中心频率为500 MHz，时窗长度为30 ns，采样时间间隔为0.02 ns.采用宽角法观测方式，发射源位于(2.0 m，2.0 m)处；接收源与发射源的水平距离R分别为0.5 m，1.0 m，1.5 m，2.0 m.

图 3为利用基于Pade逼近的Cole-Cole频散介质FETD正演算法对均匀频散介质Ⅰ模型进行计算的结果与解析解(Giannakis et al., 2012)对比图，由图可知：FETD的正演结果与解析解拟合得非常好，平均相对误差为0.35%，电磁波传播到模型边界外几乎没有出现明显的反射波，说明本文采用的PML边界条件对边界处的人工反射电磁波具有较好的吸收效果.对比结果验证了本文提出的基于Pade逼近的Cole-Cole频散介质GPR有限元正演算法的正确性，且具有较高的精度.

为进一步分析GPR在频散介质中的传播特征及探测能力，对收发距R分别为0.5 m，1.0 m，1.5 m，2.0 m的观测方式进行正演计算，并与相应非频散介质的正演结果进行对比，如图 4所示.为了定量分析电磁波在Cole-Cole频散介质中的衰减规律和子波延伸程度，对图 4中的信号取包络，并定义半振幅宽度为子波宽度，提取的非频散介质与频散介质中不同收发距的电磁波振幅分别为(1.179, 0.725, 0.576, 0.486)和(1.052, 0.373, 0.153, 0.073).可以看到，随着收发距R的增大，电磁波振幅在频散和非频散介质中都不断变小；相比于非频散介质中电磁波振幅下降趋势，频散介质中的电磁波能量下降更快.收发距R分别为0.5 m, 1.0 m, 1.5 m, 2.0 m时，频散介质中电磁波振幅与非频散介质中电磁波振幅比分别为：0.89, 0.51, 0.27, 0.15，由此可见，电磁波能量在非频散介质中衰减更强，严重降低了GPR的探测深度和能力.

图 4中提取的非频散介质中不同收发距R的子波都约为1.14 ns，而相应的频散介质中的子波宽度分别为(1.28, 1.50, 1.80, 2.19)ns.当收发距为2.0 m时，电磁波在频散介质中的子波宽度比非频散介质中的子波宽度拉伸了近2倍.由此可见，由于介质的频散，电磁波在其中传播时子波宽度不断拉伸，子波持续时间更长，有效带宽不断减小，严重降低了GPR探测的分辨率.

当收发距为1.0 m，衰减因子β分别为0.2，0.4，0.6，0.8，1.0时，利用基于Pade逼近的Cole-Cole频散介质GPR有限元正演算法对均匀频散介质模型进行计算得到的对比结果如图 5所示.由图可知，随着衰减因子β增大，电磁波振幅急剧变小；当衰减因子为1.0时的电磁波振幅约为衰减因子为0.2时的电磁波振幅的1/20.此外，随着衰减因子β增大，电磁波的初至时间不断减小，且子波的持续时间不断增大，有效带宽不断减小.对比结果表明：衰减因子是影响GPR探测能力和分辨率的重要影响因素，在其他参数相同条件下，衰减因子越大，GPR探测分辨率越低，衰减因子越小，GPR探测分辨率越高.

4.2.2 算例2

图 6为一个大小为8.0 m×5.0 m的复杂Cole-Cole频散介质模型示意图.模型被一不规则界面分为上、下两层，介电参数参考文献(Giannakis et al., 2012)进行设置.上层介质为干沙，其相对介电常数ε_{r, 1}^∞=3.7, ε_{r, 1}⁰=5.4，电导率σ₁=0.00045 S·m^-1，驰豫时间τ₁=80 ps，衰减因子β₁=0.3；下层介质为湿沙，其相对介电常数ε_{r, 2}^∞=5.6, ε_{r, 2}⁰=8.9，电导率σ₁=0.002 S·m^-1，驰豫时间τ₁=110 ps，衰减因子β₂=0.25.不规则界面中间为一个由等边三角形的两边组成的凹槽地形，等边三角形的三个顶点位置分别为(3.0 m，1.0 m)，(4.0 m，2.0 m)，(5.0 m，1.0 m).下层介质的左右两边分别放置一个充填空气的圆状异常体，两个圆状异常体的圆心分别位于(1.5 m，2.0 m)和(6.5 m，2.0 m)位置，半径都为0.5 m.基于Pade逼近的Cole-Cole频散介质GPR有限元正演算法的计算区域和吸收层被剖分为293786个非结构化三角形网格，每个三角形网格内采用线性插值基函数；PML吸收边界层厚度占0.5 m.采用等偏移距观测方式，收发天线的距离为0.01 m；GPR波脉冲激励源为雷克子波，其中心频率为500 MHz；时窗长度为30 ns，采样时间间隔为0.01 ns.

图 7a为利用基于Pade逼近的Cole-Cole频散介质GPR有限元正演算法对该模型进行计算得到的GPR剖面；图 7b为模型二中背景介质为非频散介质的GPR有限元正演剖面.由图可见：剖面两边出现了一条非常明显的水平反射波，对应模型中水平界面的反射；由于模型中存在一个三角形凹槽地形，因此在水平反射波中间出现一个“Ⅴ”字形的反射波，“Ⅴ”字形两翼分别得到延伸；在“Ⅴ”字形反射波下方出现一条明显的双曲线反射，这是由于三角形底部顶点产生绕射所致.图中左右两边都出现明显的双曲线反射波，分别对应左右两边的圆状异常体产生的绕射.对比图 7a和图 7b可知：图 7a中的水平反射波，“Ⅴ”字形反射波和圆状异常体的双曲线反射波能量都比图 7b中能量弱，说明电磁波在Cole-Cole频散介质衰减更快.图 7a中频散介质中的水平反射波约在20 ns出现，空气管状异常体顶部产生的双曲线反射波约30 ns开始出现，并且反射波持续时间更长；而图 7b中的水平反射波约在18 ns开始出现，空气管状异常体顶部产生的双曲线反射波在约25 ns开始出现，反射波持续时间更短，由此可见电磁波在频散介质中传播时速度降低.

为了更好地分析介质频散对电磁波传播的影响，分别提取图 7a和7b中水平位置1.5 m处的单道波形进行对比，如图 8所示.由图可见：相比电磁波在非频散介质中的能量，GPR在Cole-Cole频散介质中传播时直达波、水平界面的反射波和圆状异常体的反射波能量都比非频散介质中的要弱；并且随传播时间的增加，电磁波在频散介质与非频散介质的反射波能量差异越大.对比频散与非频散介质中水平界面和圆状异常体的反射波可知：频散介质中的水平界面和圆状异常体的反射波的持续时间更长，电磁波子波被严重拉伸，分辨率降低.

为了对比电磁波信号在Cole-Cole频散介质与非频散介质中传播过程中的频率变化，对图 8中的波形去除直达波后进行傅里叶变换，获得的频谱对比图如9所示.由图可见，在主频带范围内，相同频率处的电磁波在Cole-Cole频散介质中能量比在非频散介质中的能量更低.在主频200 MHz附近，相比非频散介质的电磁波能量，Cole-Cole频散介质的电磁波能量约衰减了10 dB，严重降低了GPR的探测深度和分辨率.通过对比这些频散与非频散介质的雷达剖面图，单道波形图和单道频谱图能更全面准确地认识电磁波在频散介质中的传播规律，可更精确、更可靠地指导GPR资料解释.

5 结论

(1) 从Maxwell方程组出发，在介质频散满足Cole-Cole关系条件下，采用Pade逼近算法将含有时间域分数阶导数的方程导出为一组整数阶辅助微分方程，推导了Cole-Cole频散介质中的整数阶时间域电磁波方程.然后，在复数坐标系下，通过在频率域电磁波方程中引入2个中间变量并将变换到时间域，从而以变分形式将PML边界条件加载到GPR有限元方程组中，并给出了基于Pade逼近的Cole-Cole频散介质GPR有限元正演算法.

(2) 均匀Cole-Cole频散介质模型的数值模拟结果表明：本文提出的基于Pade逼近的Cole-Cole频散介质GPR有限元计算结果与解析解拟合程度好，误差小于0.35%，验证了本文提出的算法的正确性和有效性.

(3) 设计了一个复杂GPR模型，应用基于Pade逼近的Cole-Cole频散介质GPR有限元正演算法对该模型进行正演计算，获得相应的剖面图、单道波形图和单道频谱图，将该模拟结果与相应非频散介质的计算结果进行对比，分析了电磁波在Cole-Cole频散介质中传播具有衰减更强、子波带宽增大，分辨率更低等传播规律，有助于雷达资料更可靠、更准确的解释.