基于电磁波衰减补偿的探地雷达三维逆时偏移成像

王洪华， 龚俊波， 梁值欢， 张智， 徐涛

1 桂林理工大学， 地球科学学院， 广西 桂林 541004 2 中国科学院地质与地球物理研究所， 岩石圈演化国家重点实验室， 北京 100029 3 中国科学院地球科学研究院， 北京 100029

0 引言

探地雷达(Ground Penetrating Radar, GPR)作为一种高频脉冲电磁探测技术，以其效率高、分辨率高、无损探测、实时成像等优点，而被广泛应用于工程检测和浅层勘探领域(卢成明等, 2007; Mochales et al., 2008; 刘澜波和钱荣毅, 2015; Atef et al., 2016; 廖红建等, 2016; 郭士礼等, 2019).近年来，工程检测和浅层勘探的日益细化，给GPR数据的快速高精度处理与成像带来极大挑战(程久龙等, 2010; 苏茂鑫等, 2010; 黄忠来和张建中, 2013; 王敏玲等, 2019a).为此，许多学者根据电磁波与弹性波传播规律的相似性和GPR与地震数据接收方式的相似性，将一些成熟的地震数据处理与成像方法如动校正叠加(Ebihara et al., 2000; Perroud and Tygel, 2004)、速度谱分析(Grandjean et al., 2000; Booth et al., 2011)、层析成像(Johnson et al., 2007; Chang and Alumbaugh, 2011)、绕射叠加(Feng and Sato, 2004; Aitken and Stewart, 2004)、克希霍夫偏移(Moran et al., 2000; Porsani and Sauck, 2007)、逆时偏移(Fisher et al., 1992; Leuschen and Plumb, 2001; Bradford, et al., 2018)等引入到GPR数据处理与成像中.其中，逆时偏移因其具有原理简单、计算效率快、成像精度高等优点，而被广泛应用于浅部精细结构的GPR高精度成像中，取得了良好的成像效果(傅磊等, 2014; Liu et al., 2014,2016,2017; Bradford, 2015; Lu et al., 2016;王敏玲等, 2019b; Zhang et al., 2019).

然而，电磁波与弹性波在地下介质传播规律存在显著差异：二阶弹性波动方程只涉及位移关于时间的二次导数(波动项)，然而地下介质电导率的存在，使得二阶电磁波方程还涉及电磁场关于时间的一次导数(衰减项)(Di and Wang, 2004).与弹性波在地下介质中传播时只表现波动特性相比，电磁波既有波动特性，也有表现地下介质吸收的衰减特性.高频电磁波在地下介质中传播时，传播速度和衰减系数是关于介质电导率的函数，电导率越高，传播速度越小，衰减系数越大，能量更易衰减(Bergmann et al., 1998; Neto and Mediros, 2006; 张先武等, 2014; 王洪华等, 2018).因此，在GPR逆时偏移中，考虑地下高电导率介质对电磁波的强吸收衰减作用，对提高高电导率区域的成像质量显得尤为必要.常规GPR逆时偏移在计算反传电磁波场时大都沿用正传电磁波方程，正传电磁波场在高电导率介质中衰减的同时，反传电磁波场会再次衰减，能量非常微弱，难以实现高电导率区域的清晰准确成像(朱尉强和黄清华, 2016; 王敏玲等, 2019a).如何在计算反传电磁波场的同时，对电磁波正传时衰减的能量进行精确补偿，以提高高电导率区域的成像质量是目前该领域的研究热点.Sena等(2006)在裂步-傅里叶偏移算法的基础上，通过在频率域电磁波场外推的同时进行反滤波处理，补偿衰减的电磁波场，有效提高了高衰减区域的成像质量.其后，Oden等(2007)将上述方法应用于GPR频率-波数偏移算法中，数值试验论证了该算法的有效性.然而，反滤波方法大都基于一维衰减模型，难以适用于复杂地质结构的GPR逆时偏移成像(Zhu et al., 2016; 朱尉强和黄清华, 2016).近年来，针对黏弹性介质高精度成像提出的衰减补偿逆时偏移算法为高电导率介质区域的GPR高精度成像提供了一种可行有效的方案(Zhu, 2014; Zhu and Harris, 2014).该算法通过在弹性波场反传过程中，人为改变黏弹性波动方程中衰减项的正负号，以保持逆时外推的时间对称性和反传不变性，精确补偿正传时衰减的弹性波场能量，提高黏弹性介质的成像质量(Zhu et al., 2014; Zhu and Harris, 2015; Zhu, 2016).目前，该算法在黏弹性介质高精度逆时偏移中得到广泛应用(李振春等, 2014; Sun et al., 2016; 田坤等, 2017; 刘财等, 2018; 豆辉和徐逸鹤, 2019)，取得了良好的成像效果.从数学上看，二阶电磁波方程与二阶黏弹性波动方程形式类似，都涉及衰减项，既表现了波动特性，也表现了介质的吸收衰减特性.为此，Zhu等(2016)、朱尉强和黄清华(2016)分别根据两者的相似性，成功将黏弹性波衰减补偿的逆时偏移算法应用于二维高电导率介质结构的GPR成像中，详细推导补偿电磁波正传时衰减能量的反传电磁波方程，并用数值试验论证了该算法应用于提高高电导率介质区域的成像效果的可行性和有效性.

考虑到实际GPR高频电磁波是在地下三维空间辐射传播，二维逆时偏移难以实现反射波的准确归位和绕射波的完全收敛，成像精度降低(Liu et al., 2017, 2018; 张崇明等, 2019; Zhu et al., 2020).本文在Zhu等(2016)、朱尉强和黄清华(2016)的基础上，开展基于电磁波衰减补偿的三维GPR逆时偏移成像研究.其中，三维时域有限差分法用于计算正传和反传电磁波场，并通过改变反传电磁波方程中衰减项的正负号，以补偿电磁波正传时衰减的能量；零时刻成像条件用于获得三维逆时偏移结果.数值试验论证了本文构建的基于电磁波衰减补偿的三维GPR逆时偏移算法在高电导率区域成像分辨率和抗干扰能力方面的优势.

1 方法原理

1.1 电磁波的衰减特征

根据电磁波场理论，忽略激励源的影响，三维GPR电磁波方程可表示为(冯德山等，2017)：

(1)

假定电场为时谐场，式(1)两边都进行傅里叶变换并整理，可推导电磁波复传播速度为(Carcione, 2014)：

(2)

式中，i为虚数单位，将式(2)展开，可推导电磁波在地下介质中传播的速度和衰减系数表达式为(Zhu et al., 2016)：

(3)

(4)

其中，sgn(x)为符号函数，当x>0时，sgn(x)=1；当x<0时sgn(x)=-1.

由式(3)和(4)可知，电磁波速度和衰减系数是关于介质电导率的函数.图1为均匀介质(相对介电常数为6，频率为400 MHz)中电磁波速度和衰减系数随电导率变化曲线，其中实线和虚线分别是电导率取正值和负值所得.由图可知，与电导率为0.0001 S·m-1时的电磁波速度相比，电导率为0.01 S·m-1时的电磁波速度变化约为0.1%，受电导率变化影响较小；当电导率从0.0001 S·m-1增大到0.01 S·m-1时，衰减系数从0.008增大到0.8，受电导率变化影响较大.由此可见：电导率是影响电磁波能量衰减的关键参数，特别是在高电导率区域中电磁波能量衰减更强.因此，对在高电导率区域采集的GPR数据进行逆时偏移时，补偿电磁波衰减的能量显得尤为必要.

1.2 逆时偏移方法

GPR逆时偏移原理是在构建偏移速度模型的基础上，将实测GPR信号作为边界条件在时间轴上进行逆时外推，当逆推至零时刻时应用相关成像条件获取成像结果，从而实现地下精细结构的高精度成像(王敏玲等, 2019b).根据时间反转原理(Fink, 1992; Zhu et al.,2016; 朱尉强和黄清华, 2016)，电磁波场进行逆时外推满足方程为

图1 电磁波速度(a)和衰减系数(b)随电导率变化曲线Fig.1 Curves of velocity (a) and attenuation coefficient (b) of electromagnetic waves varying with conductivity

(5)

(6)

为避免电磁波逆时外推时的衰减，补偿正传过程中电磁波衰减的能量，Zhu等(2016)提出了一种改变式(5)中衰减项前的正负号方法，人为保持电磁波场逆时外推的反转不变性，即：

(7)

式(7)与式(1)的形式一致，可精确补偿电磁波正传中衰减的能量(朱尉强和黄清华, 2016).

本文采用三维时域有限差分法进行接收点电磁波场的逆时外推，零时刻成像条件用于获取成像结果.当电磁波逆时外推至零时刻时，零时刻的电磁波场即为地下结构的成像结果(王敏玲等, 2019b)，可表示为：

(8)

2 数值计算

2.1 方法的可行性和有效性验证

为验证本文提出的三维衰减补偿电磁波场逆时外推方法的可行性和有效性，建立了一个1.5 m×1.5 m×1.5 m三维均匀模型，其相对介电常数εr=8.三维FDTD用于模拟计算时的空间步长均为0.01 m，时间步长为0.015 ns，时间长度为24 ns；激励源是中心频率为400 MHz的雷克子波.首先，将激励源放置于模型的正中心(0.75 m, 0.75 m,0.75 m)，分别将均匀模型的电导率σ设置为 0 S·m-1(无损)、0.001 S·m-1、0.015 S·m-1、-0.015 S·m-1，获得的7 ns时刻Ey分量的波场快照，如图2所示.由图可见，介质电导率越大，电磁波能量衰减更强、能量越弱.当电导率σ为-0.015 S·m-1时，即利用式(7)进行模拟计算时，保持了时间反转不变性和时间对称性，电磁波衰减得到有效补偿，与图2c中的电磁波能量相比，能量被有效恢复；且与图2a无损情况下的电磁波能量相当，说明通过改变电磁波方程中衰减项前的正负号，可有效补偿电磁波在高电导率介质中衰减的能量.

图2 不同电导率均匀模型中7 ns时刻Ey分量的波场快照(a) 0 S·m-1； (b) 0.001 S·m-1； (c) 0.015 S·m-1； (d) -0.015 S·m-1.Fig.2 Snapshots of the Ey wave field of homogenous model at 7ns with different conductivity values

图3为不同电磁波场逆时外推方法在均匀模型正中心位置处接收到的波形对比，其中灰实线为三维FDTD正演在模型正中心接收到的波形；黑点虚线、黑实线和黑虚线分别是模型电导率为0 S·m-1、0.015 S·m-1、-0.015 S·m-1时电磁波场逆时外推接收到的波形，即将模型最外层所有网格点作为接收点接收到的GPR信号进行逆时外推后在模型中心位置处接收的波形.由图3可知：常规不考虑电导率的逆时偏移无法对电磁波衰减进行补偿，如黑点虚线所示；常规考虑电导率的逆时偏移比不考虑电导率的电磁波能量衰减更强，成像结果更差；而本文构建的基于电磁波衰减补偿的逆时偏移结果与正演波形相比，能较好地补偿由电导率引起的电磁波衰减如黑色虚线所示，验证了本文构建的三维衰减补偿电磁波场逆时外推方法的可行性和有效性.

图3 不同逆时外推电磁波场重构方法在均匀模型正中心位置处接收到的波形对比灰实线为正演接收到的波形，黑点虚线、黑实线和黑虚线分别为电导率为0 S·m-1 (无损)、0.015 S·m-1、-0.015 S·m-1时逆时外推接收到的波形.Fig.3 Comparison of reconstructed waveforms at center position of homogenous model by using different reverse time extrapolation methodsThe grey line is the simulated waveform, black dot-dashed line, black line and black dotted line are the reconstructed waves of reverse time extrapolation by using the homogenous model with conductivity of 0 S·m-1, 0.015 S·m-1, and -0.015 S·m-1, respectively.

2.2 空洞模型

为验证本文构建的基于电磁波衰减补偿的三维GPR逆时偏移方法的成像效果，建立了一个大小为0.8 m×2 m×1.1 m的空洞模型，如图4所示.模型被埋深为0.5 m水平界面分为上下两层，其相对介电常数分别为6和8；下层介质的左右两边分别埋有一个大小为0.1 m×0.1 m×0.1 m的正方体空洞，其中心分别位于(0.4 m,0.5 m,0.8 m)、(0.4 m,1.5 m,0.8 m)，如图4a所示.模型的背景电导率为0.001 S·m-1，右侧设置了一个大小为0.4 m×0.4 m×0.8 m高电导率区域，其电导率为0.015 S·m-1，中心位置为(0.4 m,1.5 m,0.7 m)，如图4b所示.利用三维FDTD进行模拟计算时的参数与均匀介质模型相同，平行Y方向X=0～0.8 m之间等距布设9条测线，测线间距为0.1 m，收发天线间距为0.06 m；平行X方向Y=0～2 m之间等距布设11条测线，测线间距为0.2 m，收发天线间距为0.06 m.

图5a、b分别为空洞模型X方向和Y方向上三维GPR正演切片.由图5a可见：8 ns处，Y=1.5 m附近出现较强的反射波，这是由于高电导率区域与背景电导率差异明显，反射系数不为零所致.11 ns附近出现上、下层介质分界面产生的水平反射波，波形能量强、易识别；受高电导率区域(Y=1.5 m附近)的影响，电磁波在传播过程中出现较为明显地衰减，波形能量较弱，如X=0.3 m、0.4 m、0.5 m位置处的正演切片所示.16 ns开始出现空洞产生的双曲线绕射波，空洞正上方测线(X=0.4 m)的正演切片中绕射波能量最强，其他测线上的正演切片中绕射波能量随距离增大而变弱、出现时间变长.低电导率区域(Y=0.5 m附近)中空洞产生的绕射波能量比高电导率区域(Y=1.5 m附近)能量更强、波形更明显，这是由于电磁波能量在高电导率区域衰减更强所致.分析图5b中Y方向上的正演切片可得到类似的结论.

图4 空洞模型的示意图(a) 相对介电常数分布； (b) 电导率分布.Fig.4 Schematic diagrams of void GPR model(a) Relative permittivity； (b) Conductivity.

图5 空洞模型的三维正演剖面(a) X方向； (b) Y方向.Fig.5 3D GPR forward profile of void model(a) X direction； (b) Y direction.

利用本文构建的基于电磁波衰减补偿的三维GPR逆时偏移算法对图5所示的三维正演剖面进行逆时偏移成像，并与常规逆时偏移和介质无损情况下的逆时偏移结果进行对比，获得结果如图6所示.由图6可知：三种逆时偏移成像剖面中水平界面产生的反射波能量得到准确归位，空洞产生的绕射波完全收敛，成像结果清晰准确.但三种逆时偏移方法对高电导率区域的成像分辨率存在明显差别：图6a展示的常规三维GPR逆时偏移结果中，由于未考虑电导率对电磁波能量衰减的影响，高电导率区域处的水平界面与空洞位置处的成像非常模糊、不易被识别；这是由于电磁波在高电导率区域进行逆时外推时能量再次衰减所致.与图6a相比，图6b所示的基于电磁波衰减补偿的GPR三维逆时偏移结果中高电导率区域中衰减的电磁波能量得到较好补偿，水平界面和空洞的成像能量得到较好的恢复，成像结果更清晰、准确；且与介质无损(电导率为0)情况下的三维GPR逆时偏移结果图6c吻合较好.

为更好地分析基于电磁波衰减补偿的三维GPR逆时偏移成像方法对高电导率区成像的优势，提取X=0.4 m、Y=1.5 m的单道波形对比，如图7所示.由图7可见，三种GPR逆时偏移成像结果中，水平界面和空洞成像位置与真实位置相符；相比介质无损情况下三维GPR逆时偏移结果中的波形能量，常规GPR逆时偏移结果波形能量衰减约86%；基于电磁波衰减补偿的三维GPR逆时偏移结果中的能量得到有效恢复，且与介质无损情况下的逆时偏移结果基本吻合.由此可见：基于电磁波衰减补偿的三维GPR逆时偏移可有效补偿电磁波在高电导率介质中传播损失的能量，大大提升了目标体的成像精度和分辨率，其结果更有利于后续雷达资料的解释.

图6 空洞模型三维GPR正演数据的逆时偏移剖面(a) 常规逆时偏移结果； (b) 衰减补偿逆时偏移结果； (c) 介质无损情况下的逆时偏移结果.Fig.6 RTM imaging results of 3D forward GPR profile of void model(a) Conventional RTM； (b) Aattenuation compensated RTM； (c) Conventional RTM in lossless media.

2.3 分层界面模型

图8是一个大小为0.8 m×2 m×1.5 m的分层界面模型，从上至下分为四层，各层的介电参数和几何参数分布如图所示，在第二层介质左边存在一个局部高电导率区域，其电导率为0.012 S·m-1.利用三维FDTD对该模型进行模拟的激励源是中心频率为600 MHz的雷克子波，其余计算参数与空洞模型相同.

利用基于电磁波衰减补偿的三维GPR逆时偏移算法对该模型正演数据进行逆时偏移成像，并与常规GPR逆时偏移结果和介质无损情况下的逆时偏移结果进行对比，如图9所示.由图9可见：常规逆时偏移结果中由于电磁波在高电导率区域传播时衰减较强，附近分界面的成像模糊、能量微弱，特别是最下层的分界面难以被识别；相比常规逆时偏移结果，基于电磁波衰减补偿的逆时偏移结果中，高电导率区域附近的反射界面的成像能量得到有效补偿且与介质无损情况下的逆时偏移结果相当，成像结果清晰可见，且易被识别.

图7 图6中X=0.4 m、 Y=1.5 m的单道波形对比Fig.7 Comparison of single-trace waveforms at position of X=0.4 m, Y=1.5 m in Fig.6

图8 分层界面模型示意图(a) 相对介电常数分布； (b) 电导率分布.Fig.8 Schematic diagrams of layered interface GPR model(a) Relative permittivity； (b) Conductivity.

图9 分层界面模型三维GPR正演数据的逆时偏移剖面(a) 常规逆时偏移结果； (b) 衰减补偿逆时偏移结果； (c) 介质无损情况下的逆时偏移结果.Fig.9 RTM imaging results of the 3D forward GPR profile of layered interface model(a) Conventional RTM algorithm； (b) Attenuation compensated RTM algorithm； (c) The result conventional RTM with lossless media.

为验证基于电磁波衰减补偿的三维GPR逆时偏移算法的抗干扰能力，对X=0.4 m处的正演剖面图10a分别施加50%、70%和90%噪声后的GPR剖面如图10b、c、d所示.由图可知，噪声越强，界面产生的反射波越难以被识别，施加90%噪声后的GPR剖面中，第三层界面的反射波难以被识别.

图11是对施加噪声的GPR剖面进行三维GPR衰减补偿逆时偏移和常规逆时偏移结果对比，其中图11 a、c、e是常规逆时偏移剖面，图11 b、d、f是衰减补偿逆时偏移结果.由图可知：施加的噪声越强，成像结果越模糊，杂波干扰越强；但不同程度噪声干扰下的GPR剖面中的反射波均得到准确归位，且与其真实位置相符，也易被识别.特别是施加90%噪声情况下，部分有效波已经被严重污染，但通过逆时偏移仍然能对其进行较好成像.对比常规GPR逆时偏移和衰减补偿GPR逆时偏移结果可知，本文构建的基于电磁波衰减补偿的三维GPR逆时偏移成像算法具有较强的抗干扰能力.值得一提是，基于电磁波衰减补偿的逆时偏移在对反射波能量进行补偿的同时，也存在对噪声信号进行补偿的现象.

3 结论

(1) 本文从电磁波的衰减特性和逆时偏移原理出发，通过改变三维反传电磁波方程中包含电导率的衰减项的正负号，保持电磁波反传的时间对称性和不变性，以精确补偿电磁波在正传中衰减的能量，构建了一种基于电磁波衰减补偿的三维GPR逆时偏移算法.其中，三维FDTD用于计算正传和反传电磁波场，零时刻成像条件用于获取地下介质的成像结果.

(2) 两个典型三维GPR模型的正演剖面的电磁波衰减补偿逆时偏移和常规逆时偏移结果对比表明：本文构建的基于电磁波衰减补偿的三维GPR逆时偏移算法可精确补偿电磁波在地下介质传播时衰减的能量，高电导率区域的成像分辨率更高，抗干扰能力更强，其结果更有利于指导后续雷达剖面的解译.

图10 施加不同比例噪声的GPR剖面(a) 未施加噪声； (b) 50%； (c) 70%； (d) 90%.Fig.10 GPR profile added with noises of different proportions(a) Without noise; (b) 50%; (c) 70%; (d) 90%.

图11 图10中GPR剖面的衰减补偿逆时偏移和常规逆时偏移结果对比(a)、(c)和(e)是常规逆时偏移结果； (b)、(d)和(f)是衰减补偿逆时偏移结果.Fig.11 Comparison of attenuation compensated RTM and conventional RTM of the GPR profile shown in Fig.10(a), (c) and (e) Conventional RTM；(b), (d) and (f) Attenuation compensated RTM.