2020, Vol. 63
2. 有色金属成矿预测教育部重点实验室, 长沙 410083
2. Key Laboratory of Metallogenic Prediction of Nonferrous Metals, Ministry of Education, Changsha 410083, China
探地雷达(GPR)是一种对地下异常体、结构或物体内部不可见目标体进行定位的电磁无损探测技术(Daniels，2004；王敏玲等，2019).全波形反演(FWI)是新兴起的重要成像手段，其最初是为地震勘探中的声波和弹性波方程开发的(Lailly，1983；Tarantola，1984；Pratt et al., 1998；Xu and McMechan, 2014；Dagnino et al., 2014).目前全波形反演已在GPR成像中得到广泛关注与研究.它可以在时间域(Gloaguen et al., 2007；Zhou et al., 2007；Klotzsche et al., 2010；Cordua et al., 2012；Belina et al., 2012；吴俊军等，2014；李尧，2017；胡周文，2018)、频率域(El Bouajaji et al., 2011；Ellefsen et al., 2011；Yang et al., 2013)、拉普拉斯域求解(孟旭，2016).GPR全波形反演方法能够获取介质的介电常数、电导率，有助于岩性的精确描述和异常体预测，是重构近地表复杂地质结构的有力工具.
在常规的探地雷达反演中，通常假设介质的电导率为已知的背景值，仅反演介电常数.然而介电常数和电导率在雷达数据的定量解释和近地表地质解释中都扮演着重要角色，全波形反演中获取可靠的电导率参数已成为迫切需求.随着大规模科学计算能力的提升，多参数FWI逐渐成为可能，并迅速成为FWI的一个重要研究方向.多参数全波形反演会存在参数间的串扰，使反演问题的性态变差，增加多参数全波形反演的非线性程度，严重影响反演结果的精度.为了减轻多参数FWI的参数串扰，提高全波形反演中电导率重建的精度，Ernst等(2007)提出了一种介电常数与电导率交替更新的反演方法.毛立峰(2007)利用改进的非线性共轭梯度反演方法和改进的有限内存BFGS反演方法，对超宽带雷达数据进行了双参数反演.Meles等(2010)通过使用不同的参数步长同步更新介电常数和电导率，提高电导率的成像精度.Busch等(2012)开发了一种通过地面GPR多偏移距数据估计地下分层模型的介电常数和电导率值的FWI方案.Kalogeropoulos等(2013)将介电常数和电导率的GPR全波形反演用于混凝土中氯化物变化规律和分布的评估.Lavoué等(2014)讨论了参考介电常数、参考电导率和参数比例因子的定义范围，在二维频率域实现了相同的步长下同步更新介电常数和电导率.Nilot等(2018)通过模型平均值对参数进行变换，使用归一化因子在无记忆拟牛顿方法的框架中同步更新介电常数和电导率.Jazayeri等(2018)利用GPR全波形反演实现地下圆柱形目标体的直径和填充材料的电性参数估计.
目前全波形多参数反演在数学上常被抽象为一个二次型优化问题.从最优化求解角度可分为梯度引导类与牛顿类解法.梯度引导类算法中应用较多的有最速下降法和共轭梯度法.最速下降法直接用目标函数梯度的反方向作为搜索方向，简单易实施，Ernst等(2007)、Meles等(2010)均采用该算法求解GPR FWI问题，但其收敛速度较慢.共轭梯度法采用共轭方向代替梯度反方向作为反演的搜索方向，每次迭代所需的共轭方向由当前梯度方向与上次迭代的共轭方向组合产生，冯德山和王珣(2018)、俞海龙(2019)采用共轭梯度法求解GPR双参数反演问题，然而该算法对于双参数反演而言其收敛效果仍不理想.牛顿类算法是通过利用目标函数的一阶偏导的梯度算子和二阶偏导的Hessian矩阵进行约束，具有二阶性，在理论上比梯度类算法的收敛速度更快，但是Hessian矩阵的存储和求取的空间复杂度较高.针对牛顿算法中Hessian矩阵求取存在的问题，具有良好的数值效果和快速收敛性的BFGS算法被提出，但该算法的Hessian矩阵的存储量及计算量仍然较大.而L-BFGS(Nocedal and Wright, 2006)不需要显式形成Hessian矩阵及其逆矩阵，仅需要提供Hessian逆矩阵的一个初始近似矩阵，通过保存部分梯度算子和模型修正量的信息求取伪Hessian矩阵，大大减小了存储量和运算量，同时对多参数反演中参数反演精度具有一定的调节功能.
由于GPR全波形反演是一个高度非线性、不适定问题，对初始模型具有很高的依赖性，易陷入局部极小值及周波跳跃现象(Meles et al., 2012).多尺度策略能有效地避免反演结果陷入局部极小值，它将反演问题分解为不同尺度，采用若干个低频带到高频带的逐频带反演，根据不同尺度上的反演目标函数的特征去求解反演问题，从而逐步搜索到全局极值点(Bunks et al., 1995；Boonyasiriwat et al., 2009).为改善这种问题，Meles等(2012)提出了一种基于频域-时域分析的改进全波形新算法；Feng等(2019b)基于多尺度策略采用非结构化网格FETD正演算法对常见衬砌病害模型进行了反演成像研究；刘新彤(2019)提出了一种GPR相位波形反演策略，该反演策略通过对GPR数据和子波积分处理减弱目标函数的非线性，能有效的实现对介电常数和电导率参数的单独反演成像；Zhang等(2019)基于波场积分与子波滤波处理提出了一种改进的全波形反演方法，结合逆时偏移算法能有效改善成像质量.模型正则化同样可以减少全波形反演的不适定性，使多参数反演具有更好的稳定性(Lavoué et al., 2014；Watson，2016)；Lavoué等(2014)将Tikhonov正则化应用于GPR的全波形反演求解中改善了电导率模型中的高振幅振荡现象，然而Tikhonov正则化的反演倾向于产生平滑模型；Watson(2016)使用近似的全变差正则化约束介电常数模型，使其反演结果具有尖锐界面，这种近似全变差方法可能不稳定并且反演的收敛性不能保证；Feng等(2019a)开发了基于近似全变差正则化的GPR双参数反演算法；Ren(2018)采用交叉梯度结构约束交替反演介电常数和电导率模型，修正了参数模型的结构，提高了电导率模型的准确性.
本文采用改进的全变差(MTV)正则化方案(Lin and Huang, 2015)，实现了二维地面GPR FWI的时间域介电常数和电导率参数同步反演.该反演策略通过改进的全变差正则化方案对介电常数和电导率进行模型约束，这种改进的全变差正则化由两个正则化项组成：L2范数项和L1范数全变差项，可以改进全变差反演的鲁棒性，提高反演的精度.采用L-BFGS优化算法同步反演介电常数和电导率，它对双参数反演中介电常数与电导率的参数尺度具有调节功能，同时避免对Hessian矩阵的直接存储和精确求解，大大减小了存储量和运算量.通过加载介电常数与电导率参数调节因子，改善双参数反演中的串扰问题，结合频率多尺度策略和改进全变差正则化使优化问题逐步搜索到全局极值点，避免陷入局部极值，可以达到提高GPR多参数全波形反演中介电常数、电导率的反演精度的目的.
1 时间域正反演问题 1.1 探地雷达FDTD正演由电磁场与电磁波的理论可知(Taflove and Brodwin, 1975)，二维GPR波满足的Maxwell方程可表示为
|
(1) |
式(1)中L为正演算子，u为波场向量，j为场源：
|
(2) |
上角标T表示转置，Hx、Hz为磁场强度分量(A·m-1)，Ey为电场强度分量(V·m-1)，Jy为电流密度分量(A·m-2).系数矩阵A、B、C、D为
|
(3) |
其中，ε为介电常数(F·m-1)，μ为磁导率(H·m-1)，σ为电导率(S·m-1).由于GPR反演中需要多次调用GPR正演，因此，正演的效率及并行加速的便利性非常关键，为此，本文选用基于交错网格的时域有限差分法(FDTD)进行GPR正演，时间空间离散精度均为2阶，截断边界处采用非分裂的卷积完全匹配层(Roden and Gedney, 2000).
1.2 反演数据目标函数构建探地雷达全波形反演实质是利用已知的实测数据来重构地下介质中的模型参数：介电常数ε与电导率σ的空间分布.根据正演模拟数据与实测数据之间的拟合最优，定义数据目标函数为
|
(4) |
式(4)中，M为源的个数，N为每个源的接收器个数，T0表示雷达记录的时窗长度，rj是第j个接收器空间坐标向量，Eiobs(rj, t)是第i个源激发在rj处接收到的观测数据，Ei(m, rj, t)是第i个源对猜测模型正演计算的模拟数据，模型介质参数向量m：
|
(5) |
全波形反演就是寻求满足目标函数S(m)极小值的模型介质参数向量m.由于全波形反演计算量太大，为了减少计算量，本文采用L-BFGS算法(Nocedal and Wright, 2006)对式(4)进行求解，反演迭代过程中需要多次求解目标函数的导数.目标函数(4)的Fréchet导数为
|
(6) |
其中m的变分δm与残差vi(m, rj, t)分别为
|
(7) |
δEi(m, rj, t)是式(2)中δui位于r=rj的第三项，δui是第i个源的场向量ui的变分：
|
(8) |
因为在t=0时刻电磁波尚未传播，初始条件δui(m, r, 0)=0.
下面推导δu(m, r, t)与u(m, r, t)的关系，δu是δm引起u的变化量，根据式(2)有：
|
(9) |
|
(10) |
联立式(9)与式(10)整理可得：
|
(11) |
其中δC和δD分别为C和D变分：
|
(12) |
忽略高阶项δC∂tδu和δDδu，式(12)可化为
|
(13) |
相应的δui与ui满足如下关系：
|
(14) |
为了使目标函数的梯度可以显式表达，引入伴随场w=(Hx* Hz* Ey*)T，定义算子L*为算子L的伴随算子，根据伴随作用：
|
(15) |
其中〈, 〉表示时间和空间内积，根据定义式(15)可以化为
|
(16) |
定义伴随场wi(m, r, t)方程满足如下微分方程和终止条件：
|
(17) |
|
(18) |
其中iy是一个y方向的单位向量，δ(r－rj)为Dirac函数.
将式(14)与(17)代入式(16)可得：
|
(19) |
对式(19)左端的广义函数求体积积分得：
|
(20) |
将式(20)代入式(6)得：
|
(21) |
将式(21)右端项展开，并改写为空间内积形式有：
|
(22) |
其中：
|
(23) |
|
(24) |
地面探地雷达探测时通常只接收与测线方向垂直的电场分量，因此在反演中仅有一个电场分量Ey，无磁场分量Hx与Hz.实际的GPR全波形双参数反演过程中，为了更准确地对模型参数定性，需要对介电常数、电导率同时反演.由于介电常数、电导率在数量级上相差很大，给反演计算带来了诸多不便.因此，如何设计一个能处理不同参数单位和敏感度的多参数反演策略，是GPR全波形双参数反演的关键(Meles et al., 2012).
考虑到相对介电常数εr=ε/ε0可以根据真空介电常数来定义，因此，类似地可以引入相对电导率σr=σ/σ0，取参考介质的电导率σ0=1/η0，其中η0=120πΩ为自由空间波阻抗，可以保证相对介电常数εr和相对电导率σr处于同一个量级，然后采用Lavoué等(2014)的做法，在反演过程中引入无量纲比例因子β，将模型参数m设定为相对介电常数和相对电导率的线性组合形式(εr，σr/β)，重写尺度变化之后的模型向量和梯度向量的明确表达式为
|
(25) |
式(25)中ε0与η0为常量，β为可调整的比例因子.这样，在优化过程中通过控制σr对εr的权重，避免由相对电导率与相对介电常数定义不准确引起反演过程的不稳定性.
2 时间域GPR全波形反演优化策略GPR全波形反演是一个高度非线性、不适定问题，对初始模型具有很高的依赖性，易陷入局部极小值及周波跳跃现象.为了使全波形反演具有更好的稳定性、精度更高，本文引入了多尺度反演策略和改进全变差模型约束.其中多尺度反演策略通过从低频到高频多尺度逐步反演，避免反演陷入局部极小值，有效提高寻优的准确性.而改进全变差模型约束将可靠的先验信息加入到正则化中，通过约束反演解的迭代方向，使反演结果的信息更加完整.
2.1 多尺度串行反演策略雷达数据中的低频分量主要包含地下较大的构造体信息，无法重构与波长相比较小的细节信息.而高频分量中包含较小构造体的细节信息，易发生多次散射，非线性更强.因此，将高/低频分量结合起来进行多尺度反演，是一种较好的策略.本文采用Boonyasiriwat等(2009)文献中提出的多尺度反演策略，它将反演问题分解为不同尺度，并采用Wiener低通滤波器，对观测数据与激励源子波滤波得到低频带信息，采用2~3个低频带到高频带的逐频反演，根据不同尺度上的反演目标函数的特征去求解反演问题，从而逐步搜索到全局极值点，避免陷入局部极值.其中低通滤波器采用Wiener滤波器(Boonyasiriwat et al., 2009)：
|
(26) |
其中fwiener为频率域Wiener低通滤波器，Wtarget为目标激励源子波频谱，Woriginal(ω)为原始激励源子波频谱，上标*表示共轭，τ为一个防止分母为零的常量小数，如果该值选取过大，会导致滤波子波形态与目标子波不同，这里设τ=10－4.可以将数据变换到频率域，低通滤波到目标频段，再反变换回时间域.需要注意的是，激励源函数和观测数据都要进行低通滤波.
2.2 改进全变差模型约束解决GPR全波形反演不适定性最常用的方法是进行模型正则化，通过加入先验模型约束的正则化项，使反问题更加稳定，反演结果更接近物理实际.但由于光滑性的缘故，Tikhonov正则化方法易导致目标区域与背景区域边界模糊.全变差正则化方法能有效改善Tikhonov正则化边界的过度光滑，与背景区域区别更加明显，重构的目标体边缘轮廓更加清晰.引入了一个全变差正则化后，形成新的目标函数：
|
(27) |
其中Φm(m)为模型参数目标函数，可表示为
|
(28) |
式中λ为正则化因子，‖·‖TV表示全变差算子，对于一个二维的n×n的离散模型m1，其可以表示为
|
(29) |
TV算子对于齐次函数是不可微的，不适合基于梯度的优化算法(Watson，2016)，常采用近似的全变差算子进行代替：
|
(30) |
然而，利用具有近似的传统全变差正则化方案的全波形反演可能产生反演伪像，并且由于数据拟合项和全变差正则化项引起的非线性以及反演的收敛对平滑参数τ高度敏感，将导致反演问题十分不稳定，且不能保证反演的收敛.因此本文使用改进全变差方案提高全波形反演的精度和收敛性.引入改进的全变差正则化方法(Lin and Huang, 2015)，目标函数可定义为
|
(31) |
其中λ和γ均为正的正则化参数，v=(vε, vσ)T为辅助变量，vε和vσ表示介电常数和电导率对应的先验模型向量.可将问题进一步分解为2个交替最小化子问题：
|
(32) |
|
(33) |
这两个子问题分别起着不同的作用：第一个子问题使用基于Tikhonov正则化和多参数的先验模型vi(k)的求解m(k)的全波形反演问题；第二个子问题进行2次标准的L2-TV最小化方法求解不同参数的vi(k)，以保持反演结果mi(k)中的分界面清晰.
根据方程(32)和(33)可以发现：参数λ仅存在于方程(32)的子问题中，它的作用是保持数据目标函数和Tikhonov正则化项之间的平衡.方程(32)中的Tikhonov正则化反演能保证反演的稳定性.另一方面，γ仅用于子问题方程(33)中，该方程对于vi(k)是标准的L2-TV最小化问题，它能抑制mi(k)中的假象，并保留清晰的界面.交替迭代求解方程(32)和(33)，这不仅能增强界面的清晰度，同时保证全波形反演更加稳定.
根据文献(陈小斌等，2005), 第k次迭代的正则化因子为
|
(34) |
其中λi(0)为一个无量纲数，本文取λε(0)=0.1，λσ(0)=0.05.式(34)可以实现正则化因子的自适应选取.
改进的全变差问题需要交替求解两个最小化子问题.在第k次迭代中，根据式(32)使用辅助变量vi(k)求解模型m(k)然后根据式(33)使用mi(k)求解辅助变量vi(k).得到的vi(k)将在下一个迭代步骤中方程式中(32)使用.对于第一个迭代步骤，初始模型vi(0)=mi(0).由此可以得到一个反演序列：
|
(35) |
本文使用L-BFGS方法求解方程式mi(k).使得对于第k迭代：
|
(36) |
式中a(k)为搜索步长，p(k)=－B(k)g0(k)表示搜索方向，其中g0(k)表示目标函数Φ0(m)的梯度，B(k)为逆Hessian矩阵的近似，可以用两个循环递归算法计算(Nocedal and Wright, 2006)，步长a(k)采用强Wolfe准则的非精确线搜索方法计算.在每次迭代更新中需要对目标函数梯度进行计算：
|
(37) |
子问题式(33)是一个标准的L2-TV全变差最小化问题，本文采用Split-Bregman方法(SB-TV，Goldstein and Osher, 2009)求解方程式(33)中的vi(k).
3 反演影响因素分析 3.1 多尺度反演策略的影响设置图 1所示的模型，模型背景为均匀介质，其相对介电常数和电导率分别为4和3 mS·m-1，在该均匀介质中包含了6个字母形状异常体.红色异常体的相对介电常数为8，电导率为10 mS·m-1，蓝色异常体的相对介电常数为1，电导率为0.模拟区域的边界设有20个“×”型点源和100个“o”表示的接收器，激励源为200 MHz的Ricker子波，时窗长度为75 ns，离散网格的空间步长为0.05 m，时间步长为0.1 ns.本例反演中采用均匀背景介质作为初始模型.
|
图 1 (a) 模型1的相对介电常数和(b)电导率分布 “×”和“o”符号代表发射机和接收机的位置. Fig. 1 Relative permittivity (a) and conductivity (b) distributions of model 1 Symbols "×"and "o"represent positions of sources and receivers, respectively. |
为了计算目标函数中的残差及多尺度反演中的高低频雷达数据，先用有限差分法对预设模型开展正演.正演时将Ricker子波激励源置于(2 m, 0 m)，接收天线布设100道，其中1~25道、26~50道、51~75道、76~100道的接收天线分别为上边界z=0 m、右边界x=5 m、下边界z=5 m、左边界x=0 m处，均匀等间距自左至右接收雷达信号.图 2a、b分别为模型1的真实模型和初始均匀模型的正演剖面图，图 2c为图 2a、b两者相减的初始残差，它是由模型异常体形成的记录.图 3a蓝色点划线表示中心频率为200 MHz的原始Ricker子波，黑色实线表示频率为80 MHz的目标Ricker子波，红色虚线表示经Wiener滤波器滤波后所得的子波.图 3b为图 2a中的时域共剖点雷达记录经过Wiener滤波后得到80 MHz数据.
|
图 2 模型1的200 MHz中心频率的雷克子波的时域共炮点记录 (a)真实模型；(b)初始均匀模型；(c)残差. Fig. 2 Time-domain shot gathers computed for model 1 at position (2, 0) with the time-derivative of a Ricker wavelet of central frequency 200 MHz (a) True model; (b) Initial homogeneous model; (c) Residuals. |
|
图 3 (a) 蓝色短划线是200 MHz的原始Ricker子波，黑色实线是80 MHz的目标Ricker子波，红色点划线是Wiener滤波所得的波形；(b) Wiener滤波后得到的80 MHz时域共炮点记录 Fig. 3 (a) The blue dotted line is the original 200 MHz Ricker wavelet, the black solid line is the target 80 MHz Ricker wavelet, and the red dash is the low-pass filtered wavelets such that the peak frequency is about 100 MHz; (b) Shot gather low-pass filtered by the Wiener filter to the frequency band of an 80 MHz Ricker wavelet |
为分析多尺度策略的影响，以模型1为例进行单尺度和多尺度反演计算.反演过程中首先令β=1，λ=0(未进行模型约束)，仅仅探讨单尺度与多尺度对反演效果的影响，从而限制其他因素的干扰.单尺度和多尺度反演都在4核CPU的微机上并行运算，初始模型都设为均匀背景模型，反演终止条件有两个，一个是相对目标函数Φ0(k)/Φ0(0)的比值小于1e-5，二是达到最大迭代次数200.单尺度反演直接对200 MHz原始数据进行反演，其介电常数与电导率的反演结果分别如图 4a、d所示；多尺度反演时，使用多个频带的源和数据，分别采用Wiener低通滤波后的80 MHz起始频率和200 MHz的原始频带.数据首先在低频带反演50次；然后再进行高频带反演，直至达到反演停止标准.多尺度的仅低频反演结果如图 4b、e所示，而多尺度综合双参数反演结果如图 4c、d所示.
|
图 4 模型1不同尺度策略的反演结果 (a)和(d)给出了应用原始数据的单尺度反演结果；(b)和(e)给出了应用滤波数据低频带的反演结果；(c)和(f)给出了多尺度串行反演策略的最终反演结果. Fig. 4 Inversion results of model 1 with different scale strategies (a) and (d) Single-scale inversion using the original data; (b) and (e) Multiple-scale inversion using the low-frequency filtering data; (c) and (f) Final inversion results of the multiple-scale inversion. |
在图 4a的介电常数反演剖面中，低介电常数异常体的形态重构较好，而高介电常数异常体畸变严重，无法准确识别异常形态和幅值.图 4d的电导率反演剖面的异常体形态存在严重扭曲，大小、空间位置与真实模型存在较大偏差.多尺度的低频介电常数反演剖面图 4b与图 4e电导率反演剖面都能刻画出异常体的大致轮廓，但异常体形态出现了一定程度的畸变.图 4c与图 4f的多尺度串行反演结果的介电常数与电导率剖面都能很好地重构字母异常体，其形态、大小、空间位置都能准确反映，反演效果最好，但是背景中仍然存在反演伪像.对比单尺度反演、多尺度仅低频反演和多尺度反演可以发现：双参数反演时，介电常数的反演结果较电导率的反演效果更好，而多尺度综合反演由于充分结合了高、低频的信息，对异常体的大小、形态、空间位置重构更准确，分界面更加明显，异常体的轮廓更清晰可辨，反演出的背景及异常体的电性参数与真实模型基本一致.
为了更精确地描述几种反演方式在反演精度方面的异同，在模型图 1中选取两条截线：AA′与BB′，截线分别位于深度z=1.6 m、z=3.25 m处.图 5a、c为深度z=1.6 m AA′处的介电常数和电导率单道反演切线图，图 5b、d为深度z=3.25 m BB′处的介电常数和电导率单道反演切线图.对比介电常数单道切线，黑色虚线表示的单尺度介电常数反演曲线在高介电常数位置与黑色实线表示的真实模型幅值存在较大的偏差，曲线振荡；而蓝色实线表示的低频反演曲线波峰位置与黑色实线表示的真实模型能较好地吻合，说明异常体形态归位比较准确，但幅值存在一定偏差，重构结果有待提高.红色实线表示多尺度反演曲线与黑色实线表示的真实模型曲线能较好地拟合，无论是幅值大小、曲线拐点位置都能较好地对应，存在轻微振荡.对比电导率反演结果，黑色虚线与蓝色实线与黑色实线存在较大的偏差，说明电导率的反演精度不高，而红色实线与黑色实线位置及幅值都能大致拟合，说明多尺度反演的电导率精度能够得到有效保证.两个不同深度截取的单道反演曲线对比表明，相对于单尺度反演，多尺度策略的反演结果与真实模型更为接近，特别是介电常数与真实模型基本吻合，电导率反演结果趋势一致，界面稍显模糊.
|
图 5 模型1不同尺度策略的参数截线 (a) AA′的介电常数截线；(b) BB′的介电常数截线；(c) AA′的电导率截线；(d) BB′的电导率截线. Fig. 5 Parameter sections in model 1 with different scale strategies (a) Section of permittivity at AA′; (b) Section of permittivity at BB′; (c) Section of conductivity at AA′; (d) Section of conductivity at BB′. |
图 6a为归一化的数据目标函数曲线图，图中可见，红色实线表示的多尺度的目标函数比黑色虚线表示的单尺度目标函数低一个数量级，进一步说明了应用多尺度策略的反演具有更高的精度.图 6b为模型重构误差(‖mk－mtrue‖2/‖m0－mtrue‖2)曲线，图中实线表示多尺度重构误差曲线，虚线表示单尺度重构误差曲线，而红色代表电导率，蓝色代表介电常数，结合4条误差曲线可知，相比单尺度反演曲线，多尺度反演结果与真实模型更为接近，可以收敛到全局最小，具有比单尺度波形反演成像更快的收敛速度，拟合效果更好.而且，相比电导率重构误差，介电常数重构误差更少，结果与真实模型更为接近.
|
图 6 模型1不同反演策略的收敛性和误差分析 (a)归一化目标函数；(b)模型重建误差. Fig. 6 Convergence and error analysis of different inversion strategies used for model 1 (a) Normalized objective function; (b) Model reconstruction error. |
双参数反演过程中，由于介电常数与电导率在数值上相差较大，同步反演时需要引入参数调节因子β，合理的β值可以保证反演快速稳定的收敛到全局极小值，因此，最优的β值选取非常重要.为了进一步探讨参数调节因子β对反演效果的影响，下面通过选取一组β值进行反演测试，反演过程中采用与3.1节相同的多尺度策略及终止条件，将反演结果的最终数据目标函数和模型重构误差绘制在图 7中.图中可见，当β∈[1, 1.5]时，最终的数据目标函数以及介电常数和电导率的模型重构误差达到一个合理的区间.而图 7a中的最终数据目标函数随β因子变化趋势与图 7b中的重构误差趋势基本相同.因此在实际反演中，可以开展多个β值的反演测试，根据反演结果的最终数据拟合差选择合适的β值.当β=1.25时，最终的数据目标函数值达到最小.
|
图 7 模型1使用不同参数调整因子β的收敛性和误差分析 (a)归一化的目标函数；(b)模型重构误差. Fig. 7 Convergence and error analysis of based on different parameter adjustment factors β used in model 1 (a) Normalized objective function; (b) Model reconstruction error. |
图 8a—c为β=0.5, 1.25, 2.0时介电常数反演结果，图 8d—f为电导率反演结果.分析图 8a、b、c可知，不同的β取值对介电常数的反演结果仅有细微的影响.观察电导率反演剖面图 8d、e、f，随着β取值的增大，电导率参数在反演中所占权重会逐步增大，反演结果逐步向真实模型接近，当β取值过大时，电导率反演结果会发生畸变，所有反演结果中均存在一定的反演伪像.对比图 6与图 7可以发现，β取1与1.25时的反演结果区别不大，因此β的取值并不能进一步改善反演效果.
|
图 8 模型1使用的不同参数调整因子β的反演结果 (a)和(d)是当β=0.5时的反演结果; (b)和(e)是当β=1.25时的反演结果; (c)和(f)是β=2.0时的反演结果. Fig. 8 Inversion results based on different parameter adjustment factors β used in model 1 (a) and (d) β=0.5; (b) and (e) β=1.0; (c) and (f) β=2.0. |
为了进一步改善反演效果，在反演过程中加载全变差模型约束.仍以3.1节模型为例，设置三组双参数反演实验：(S1)不加载正则化项；(S2)加载传统的全变差正则化项；(S3)加载改进全变差正则化项.反演过程中取β=1.25，其他的设置与3.2节一致.图 9a—c为三组不同正则化方案的介电常数反演剖面图，图 9d—f为三组不同正则化方案的电导率反演剖面图；图 10为三组不同正则化方案的单道切线图.综合对比反演结果图可知，方案S1的图 8a、d反演效果最差，异常体分界面模糊且存在振荡现象；方案S2的图 9b、e反演效果相对较好，异常体分界面更加清晰，振荡现象相对减小；方案S3的图 9c、f反演图像最清晰，分界面最为明显，振荡现象明显减弱，异常体重构效果最好.若不加载正则化，由于目标函数仅强调数据的拟合，导致目标函数的非线性，全变差正则化项的引入可以使数据拟合与模型约束之间达到更好的平衡，而传统的全变差正则化方法虽然在一定程度上压制了背景的非物理振荡，但是整体效果还有待提高，改进的全变差正则化方法更好地压制了背景的非物理振荡，提高反演剖面的分辨率，使得分界面更加明显，同时改进的全变差正则化方法使介质电导率幅值更加接近真实模型，在一定程度上提高了电导率模型的反演质量.
|
图 9 模型1使用的不同正则化方案的反演结果 (a)和(d)是未加载正则化时的反演结果; (b)和(e)是传统全变差正则化方案的反演结果; (c)和(f)是改进全变差方案的反演结果. Fig. 9 Inversion results of different regularization schemes used by model 1 (a) and (d) Without regularization; (b) and (e) Traditional total variation regularization scheme; (c) and (f) Improved total variation regularized inversion. |
|
图 10 模型1不同正则化方案的参数截线 (a) AA′的介电常数截线；(b) BB′的介电常数截线；(c) AA′的电导率截线；(d) BB′的电导率截线. Fig. 10 Parameter sections in model 1 with the different regularization schemes (a) Permittivity at AA′; (b) Permittivity at BB′; (c) Conductivity at AA′; (d) Conductivity at BB′. |
为了比较改进全变差正则化方法与传统全变差正则化方法对于反演收敛性的影响，并定量说明改进全变差正则化方法的重构效果，表 1为三组不同正则化方案反演结果参数，传统全变差正则化方法模型重构误差较不加载正则化有所降低，由于平滑参数的引入，导致模型约束对平滑参数高度敏感，使反演更加不稳定；而改进全变差正则化方法使用先验模型约束，从而保证了反演的稳定性，使反演更加快速稳定收敛，同时显著降低了模型重构误差.
|
表 1 模型1不同正则化方案反演结果参数 Table 1 Parameters of inversion results the different regularization schemes in model 1 |
为了更好的模拟实际情况，设计如图 11所示的模型进行反演测试，模型区域为10 m×5 m，其中图 11a的相对介电常数模型的变化范围设为3~30区间，图 11c的电导率模型的变化范围在0~20 mS·m-1区间.地表之上设置0.5 m的空气层，上层复杂分层区域代表地表相对干燥的砂层回填覆盖区域；深度约4 m的界面代表了地下水位；地下具有变化范围较大的不规则层状结构，模型中左侧浅部存在两个强衰减层，电导率为10 mS·m-1，这可能掩盖下面的结构；模型中部存在两个透镜体.
|
图 11 模型2的真实的相对介电常数(a)和电导率分布(c)，反演初始模型的介电常数(b)和电导率(d) 红色“×”和绿色“o”符号代表发射机和接收机的位置. Fig. 11 True relative permittivity (a) and conductivity (c) distributions of model 2, initial relative permittivity (b) and conductivity (d) distributions Symbols "×"and "o" represent the positions of sources and receivers. |
离散网格步长为0.05 m，将该模型区域离散为111×201网格(不含PML区域)，通过限定反演中空气层中的参数，避免源和接收器位置处的奇异性，待求的地下介电常数和电导率参数为20301个.采用多偏移距共炮点的GPR数据采样方式，设置21个源，水平间隔0.5 m，101个接收器，水平间隔0.1 m，源和接收器均位于地面之上两个网格点-0.1 m深处.反演的介电常数与电导率的初始模型设置如图 11b、d所示，它们仅仅描绘了介质参数变化的大致趋势，是图 11a、c中的真实模型高斯平滑后的结果.
为了验证该算法对噪声的适应性，将白噪声添加到正演数据中，含噪数据的信噪比(SNR)为25 dB.图 12为β=0.4时利用改进全变差正则化多尺度反演含噪数据的重构结果，在图 12a介电常数图像中，该反演曲线没有对噪声信息进行拟合，得到了较为准确的浅部构造信息，尽管浅部层中含有一些振荡，但浅部的分层仍然相当准确地重建，分辨率随着深度增加显著降低，深部的界面无法识别，介电常数幅值与真实模型存在差距.图 13b为电导率反演剖面图，图中可见，地下结构的主要趋势在电阻率剖面图中也得到较好的反映，说明该反演方法具有一定的抗噪性，但反演的电导率深度曲线与真实模型存在偏差，图像重构效果比介电常数差.
|
图 12 左边：β=0.4时利用MTV多尺度策略反演噪声数据得到的介电常数(a)和电导率(b)反演结果；右边：沿图 11中AA′、BB′和CC′处的一维深度切线 黑色曲线表示真实模型，红色曲线表示反演结果，蓝色曲线表示初始模型. Fig. 12 Left: Permittivity (a) and conductivity models (b) obtained by the inversion of noisy data using a scaling factor β=0.4 and multiple-scale with MTV schema. Right: Vertical logs extracted along the white lines indicated on the AA′, BB′ and CC′ of Fig. 11 Black curves denote the true model, red curves denote the inverted model and blue curves denote the initial model. |
|
图 13 用来反演的B-scan雷达实测数据 Fig. 13 B-scan measured radar data used in inversion |
图 13所示的归一化B-Scan剖面为实验室沙坑管线模型的实际测量数据，在填满均匀干石英砂的砂槽中埋设2根塑料空气管，直径分别为0.16 m和0.11 m，埋深分别为0.19 m和0.27 m.采用美国GSSI公司生产的SIR-4000型探地雷达仪器900 MHz天线进行测量，共采集42道数据，道间距为5 cm，每道1024个采样点.从剖面图中我们无法直接准确识别管线的材质、精确埋深和管径大小.
为测试本文算法的实用性，对该数据进行反演测试.反演前，采用周辉等(2014)所提出的方法对图 13所示观测数据进行预处理，并采用双曲线拟合方法计算得背景介质模型的相对介电常数为3.8；反演过程设置内部反演网格为217×80，其中10层为空气层，网格大小为0.01 m，反演初始模型为εr=3.8，σ=1 mS·m-1的均匀背景模型.
图 14分别为不加载模型约束、传统全变差以及改进全变差的双参数反演结果.总体而言，从三者反演结果图中均可以准确判断两个管线的材质、埋深和管径参数，反演结果正确.对比不同正则化策略的反演结果，可以发现改进全变差正则化能有效压制噪声和背景振荡，重构结果清晰，与真实物理模型更为接近，表明改进全变差正则化策略能改善反演的稳定性，提高结果的精度，对实测数据具有更好的适用性.
|
图 14 实测数据使用的不同正则化方案的反演结果 (a)和(d)是未加载正则化时的反演结果; (b)和(e)是传统全变差正则化方案的反演结果; (c)和(f)是改进全变差方案的反演结果. Fig. 14 Inversion results of different regularization schemes used by measured data (a) and (d) Without regularization; (b) and (e) Traditional total variation regularization scheme; (c) and (f) Improved total variation regularized inversion. |
(1) 提出了一种多尺度全变差正则化GPR全波形反演算法，并将该算法成功用于雷达数据的介电常数与电导率双参数同步反演中.通过引入参数调节因子，实现了相同的下降步骤内同时更新介电常数和电导率值，结合并行方式在普通微机上实现了高精度的GPR反演.
(2) 将多尺度反演策略应用到探地雷达全波形反演中，它将反演问题分解为不同尺度，先将雷达数据进行低频带反演，再将低频带的反演结果作为高频带反演的初始模型，通过采用2~3个低频带到高频带的逐频反演，从而逐步搜索到全局极值点，避免陷入局部极值，与单尺度的反演结果对比表明，多尺度反演具有更高的反演精度.
(3) 在双参数全波形雷达反演中引入了参数调节因子β及全变差模型约束，并详细探讨了两者对反演效果的影响.通过选取合适的参数调节因子，能够有效解决介电常数与电导率在数值上相差较大的问题；而改进全变差正则化模型约束的引入可以使数据拟合与模型约束之间达到更好的平衡，使反演更加快速稳定收敛.
(4) 开展了添加SNB=25db噪声的合成数据和实测数据的全波形反演测试，结果表明：本文提出的反演算法对于含噪声的合成数据以及实测数据均具有良好的适应性，能够较好地重构介电常数与电导率的参数分布，反演图像界面清晰，复杂模型中透镜体、断层以及物理模型中的管线都得到较好还原，可在探地雷达反演中推广使用.
附录将式(2)代入式(16)右端项有：
|
(A1) |
对式(A1)右边第一项有：
|
(A2) |
根据电磁波的衰减特性，有
|
(A3) |
同理可以求得式(A2)右边第二项：
|
(A4) |
对式(A2)右边第三项做类似的处理：
|
(A5) |
根据初始条件式δu|t=0=0，设伴随场w满足终止条件w|t=T=0，因此式(A5)可以化为
|
(A6) |
将式(A3)与(A4)代入式(16)可得：
|
(A7) |
因此：
|
(A8) |
Belina F A, Irving J, Ernst J R, et al. 2012. Waveform inversion of crosshole georadar data:influence of source wavelet variability and the suitability of a single wavelet assumption. IEEE Transactions on Geoscience and Remote Sensing, 50(11): 4610-4625. DOI:10.1109/TGRS.2012.2194154 |
Boonyasiriwat C, Valasek P, Routh P, et al. 2009. An efficient multiscale method for time-domain waveform tomography. Geophysics, 74(6): WCC59-WCC68. DOI:10.1190/1.3151869 |
Bunks C, Saleck F M, Zaleski S, et al. 1995. Multiscale seismic waveform inversion. Geophysics, 60(5): 1457-1473. DOI:10.1190/1.1443880 |
Busch S, Van der Kruk J, Bikowski J, et al. 2012. Quantitative conductivity and permittivity estimation using full-waveform inversion of on-ground GPR data. Geophysics, 77(6): H79-H91. DOI:10.1190/geo2012-0045.1 |
Chen X B, Zhao G Z, Tang J, et al. 2005. An adaptive regularized inversion algorithm for magnetotelluric data. Chinese Journal of Geophysics (in Chinese), 48(4): 937-946. DOI:10.3321/j.issn:0001-5733.2005.04.029 |
Cordua K S, Hansen T M, Mosegaard K. 2012. Monte Carlo full-waveform inversion of crosshole GPR data using multiple-point geostatistical a priori information. Geophysics, 77(2): H19-H31. DOI:10.1190/geo2011-0170.1 |
Dagnino D, Sallarès V, Ranero C R. 2014. Scale-and parameter-adaptive model-based gradient pre-conditioner for elastic full-waveform inversion. Geophysical Journal International, 198(2): 1130-1142. DOI:10.1093/gji/ggu175 |
Daniels D J. 2004. Ground Penetrating Radar. 2nd ed. UK: IEE.
|
El Bouajaji M, Lanteri S, Yedlin M. 2011. Discontinuous Galerkin frequency domain forward modelling for the inversion of electric permittivity in the 2D case. Geophysical Prospecting, 59(5): 920-933. |
Ellefsen K J, Mazzella A T, Horton R J, et al. 2011. Phase and amplitude inversion of crosswell radar data. Geophysics, 76(3): J1-J12. |
Ernst J R, Maurer H, Green A G, et al. 2007. Full-waveform inversion of crosshole radar data based on 2-D finite-difference time-domain solutions of Maxwell's equations. IEEE Transactions on Geoscience and Remote Sensing, 45(9): 2807-2828. DOI:10.1109/TGRS.2007.901048 |
Feng D S, Cao C, Wang X. 2019a. Multiscale full-waveform dual-parameter inversion based on total variation regularization to on-ground GPR data. IEEE Transactions on Geoscience and Remote Sensing, 57(11): 9450-9465. DOI:10.1109/TGRS.2019.2926626 |
Feng D S, Wang X. 2018. Fast Ground Penetrating Radar double-parameter inversion based on GPU-parallel by time-domain full waveform optimization conjugate gradient method. Chinese Journal of Geophysics (in Chinese), 61(11): 4647-4659. DOI:10.6038/cjg2018L0531 |
Feng D S, Wang X, Zhang B. 2019b. Improving reconstruction of tunnel lining defects from ground-penetrating radar profiles by multi-scale inversion and bi-parametric full-waveform inversion. Advanced Engineering Informatics, 41: 100931. DOI:10.1016/j.aei.2019.100931 |
Gloaguen E, Giroux B, Marcotte D, et al. 2007. Pseudo-full-waveform inversion of borehole GPR data using stochastic tomography. Geophysics, 72(5): J43-J51. DOI:10.1190/1.2755929 |
Goldstein T, Osher S. 2009. The split Bregman method for L1-regularized problems. SIAM Journal on Imaging Sciences, 2(2): 323-343. DOI:10.1137/080725891 |
Hu Z W. 2018. Research on GPR concrete nondestructive testing technology based on full waveform inversion[Master's thesis] (in Chinese). Beijing: China University of Geosciences.
|
Jazayeri S, Klotzsche A, Kruse S. 2018. Improving estimates of buried pipe diameter and infilling material from ground-penetrating radar profiles with full-waveform inversion. Geophysics, 83(4): H27-H41. DOI:10.1190/geo2017-0617.1 |
Kalogeropoulos A, Van Der Kruk J, Hugenschmidt J, et al. 2013. Full-waveform GPR inversion to assess chloride gradients in concrete. NDT & E International, 57: 74-84. |
Klotzsche A, van der Kruk J, Meles G A, et al. 2010. Full-waveform inversion of cross-hole ground-penetrating radar data to characterize a gravel aquifer close to the Thur River, Switzerland. Near Surface Geophysics, 8(6): 635-649. DOI:10.3997/1873-0604.2010054 |
Lailly P. 1983. The seismic inverse problem as a sequence of before-stack migrations.//Conference on Inverse Scattering: Theory and Applications. Philadelphia, PA: SIAM, 206-220.
|
Lavoué F, Brossier R, Métivier L, et al. 2014. Two-dimensional permittivity and conductivity imaging by full waveform inversion of multi offset GPR data:A frequency-domain quasi-Newton approach. Geophysical Journal International, 197(1): 248-268. DOI:10.1093/gji/ggt528 |
Li Y. 2017. The crosshole GPR ahead prospecting method and its application for adverse geology in tunnel construction[Ph.D. thesis] (in Chinese). Ji'nan: Shandong University.
|
Lin Y Z, Huang L J. 2015. Acoustic-and elastic-waveform inversion using a modified total-variation regularization scheme. Geophysical Journal International, 200(1): 489-502. |
Liu X T. 2019. Research on multi-scale waveform inversion method of ground penetrating radar[Ph.D. thesis] (in Chinese). Changchun: Jilin University.
|
Mao L F. 2007. The forward simulation and inversion imaging of the ultra wide band electromagnetic method[Ph.D. thesis] (in Chinese). Chengdu: Chengdu University of Technology.
|
Meles G, Greenhalgh S, Van der Kruk J, et al. 2012. Taming the non-linearity problem in GPR full-waveform inversion for high contrast media. Journal of Applied Geophysics, 78: 31-43. DOI:10.1016/j.jappgeo.2011.12.001 |
Meles G A, Van der Kruk J, Greenhalgh S A, et al. 2010. A new vector waveform inversion algorithm for simultaneous updating of conductivity and permittivity parameters from combination crosshole/borehole-to-surface GPR data. IEEE Transactions on Geoscience and Remote Sensing, 48(9): 3391-3407. DOI:10.1109/TGRS.2010.2046670 |
Meng X. 2016. Time domain and frequency domain waveform inversion of crosshole radar data and comparative research[Ph.D. thesis] (in Chinese). Changchun: Jilin University.
|
Nilot E, Feng X, Zhang Y, et al. 2018. Multiparameter Full-waveform inversion of on-ground GPR using Memoryless quasi-Newton (MLQN) method.//2018 17th International Conference on Ground Penetrating Radar (GPR). Rapperswil, Switzerland: IEEE, 1-4.
|
Nocedal J, Wright S J. 2006. Numerical Optimization. 2nd ed. New York: Springer.
|
Pratt R G, Shin C, Hicks G J. 1998. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion. Geophysical Journal International, 133(2): 341-362. DOI:10.1046/j.1365-246X.1998.00498.x |
Ren Q C. 2018. Inverts permittivity and conductivity with structural constraint in GPR FWI based on truncated Newton method. Journal of Applied Geophysics, 151: 186-193. DOI:10.1016/j.jappgeo.2018.02.025 |
Roden J A, Gedney S D. 2000. Convolution PML (CPML):An efficient FDTD implementation of the CFS-PML for arbitrary media. Microwave and Optical Technology Letters, 27(5): 334-339. DOI:10.1002/1098-2760(20001205)27:5<334::AID-MOP14>3.0.CO;2-A |
Taflove A, Brodwin M E. 1975. Numerical solution of steady-state electromagnetic scattering Problems using the time-dependent Maxwell's equations. IEEE Transactions on Microwave Theory and Techniques, 23(8): 623-630. DOI:10.1109/TMTT.1975.1128640 |
Tarantola A. 1984. Inversion of seismic reflection data in the acoustic approximation. Geophysics, 49(8): 1259-1266. DOI:10.1190/1.1441754 |
Wang M L, Liang Z H, Wang H H, et al. 2019. Review of reverse time migration in ground penetrating radar. Progress in Geophysics (in Chinese), 34(5): 2087-2096. DOI:10.6038/pg2019CC0169 |
Watson F M. 2016. Better imaging for landmine detection: An exploration of 3D full-wave inversion for ground-penetrating radar[Ph.D. thesis]. Manchester, UK: Manchester University.
|
Wu J J, Liu S X, Li Y P, et al. 2014. Study of cross-hole radar tomography using full-waveform inversion. Chinese Journal of Geophysics (in Chinese), 57(5): 1623-1635. DOI:10.6038/cjg20140525 |
Xu K, McMechan G A. 2014. 2D frequency-domain elastic full-waveform inversion using time-domain modeling and a multistep-length gradient approach. Geophysics, 79(2): R41-R53. DOI:10.1190/geo2013-0134.1 |
Yang X, Klotzsche A, Meles G, et al. 2013. Improvements in crosshole GPR full-waveform inversion and application on data measured at the Boise Hydrogeophysics Research Site. Journal of Applied Geophysics, 99: 114-124. DOI:10.1016/j.jappgeo.2013.08.007 |
Yu H L. 2019. A time domain full waveform inversion of GPR based on Modified PRP conjugate gradient method[Master's thesis] (in Chinese). Changchun: Jilin University.
|
Zhang F K, Liu B, Liu L B, et al. 2019. Application of ground penetrating radar to detect tunnel lining defects based on improved full waveform inversion and reverse time migration. Near Surface Geophysics, 17(2): 127-139. DOI:10.1002/nsg.12032 |
Zhou H, Sato M, Takenaka T, et al. 2007. Reconstruction from antenna-transformed radar data using a time-domain reconstruction method. IEEE Transactions on Geoscience and Remote Sensing, 45(3): 689-696. DOI:10.1109/TGRS.2006.888140 |
Zhou H, Chen H M, Li Q Q, et al. 2014. Waveform inversion of ground-penetrating radar data without needing to extract exciting pulse. Chinese Journal of Geophysics (in Chinese), 57(6): 1968-1976. DOI:10.6038/cjg20140627 |
陈小斌, 赵国泽, 汤吉, 等. 2005. 大地电磁自适应正则化反演算法. 地球物理学报, 48(4): 937-946. DOI:10.3321/j.issn:0001-5733.2005.04.029 |
冯德山, 王珣. 2018. 基于GPU并行的时间域全波形优化共轭梯度法快速GPR双参数反演. 地球物理学报, 61(11): 4647-4659. DOI:10.6038/cjg2018L0531 |
胡周文. 2018.基于全波形反演的GPR混凝土无损检测技术研究[硕士论文].北京: 中国地质大学(北京).
|
李尧. 2017.隧道施工不良地质跨孔雷达超前探测方法与工程应用[博士论文].济南: 山东大学.
|
刘新彤. 2019.探地雷达多尺度波形反演方法研究[博士论文].长春: 吉林大学.
|
毛立峰. 2007.超宽带电磁法正演模拟及反演成像[博士论文].成都: 成都理工大学.
|
孟旭. 2016.时间域和频率域跨孔雷达波形反演及比较研究[博士论文].长春: 吉林大学.
|
王敏玲, 梁值欢, 王洪华, 等. 2019. 探地雷达逆时偏移成像方法研究现状及进展. 地球物理学进展, 34(5): 2087-2096. DOI:10.6038/pg2019CC0169 |
吴俊军, 刘四新, 李彦鹏, 等. 2014. 跨孔雷达全波形反演成像方法的研究. 地球物理学报, 57(5): 1623-1635. DOI:10.6038/cjg20140525 |
俞海龙. 2019.基于修正PRP共轭梯度法的探地雷达时间域全波形反演[硕士论文].长春: 吉林大学.
|
周辉, 陈汉明, 李卿卿, 等. 2014. 不需提取激发脉冲的探地雷达波形反演方法. 地球物理学报, 57(6): 1968-1976. DOI:10.6038/cjg20140627 |