2018, Vol. 61
2. 有色金属成矿预测教育部重点实验室, 长沙 410083
2. Key Laboratory of Metallogenic Prediction of Nonferrous Metals, Ministry of Education, Changsha 410083, China
探地雷达(Ground penetrating radar，GPR)是一种对地下异常体、结构和特性或物体内部不可见目标体进行定位的电磁无损探测技术，它在介质中传播时，其传播路径、波场强度与波形随所通过介质的电性及几何形态而变化(李大心，1994).实际GPR探测中仅仅依靠获取的雷达剖面，即可大致推断探测目标体的位置、大小与尺寸，但多局限于粗略估计，尚不能给出异常体的介电常数、电导率等参数，从而精准界定目标体的属性.而时域全波形GPR反演是新兴起的重要成像手段，它起源于时域地震成像(Lailly，1983；Tarantola，1984)及近代地震数据处理领域(胡光辉等，2014；刘文革和贺振华，2015).它直接利用振幅、走时和相位等雷达波场信息，将正演的雷达记录与实际观测的雷达波形进行匹配，通过使两者的数据之差最小建立目标函数，从而寻找最佳模型参数.显然，全波形反演能够获取介质的介电常数、电导率等物性参数，有助于岩性的精确描述和异常体预测，提高GPR解释精度.目前全波形反演已在钻孔雷达成像中得到广泛关注与研究，分为时间域全波形钻孔雷达数据反演(Ernst et al., 2007；Gloaguen et al., 2005, 2007；Zhou et al., 2007；Meles et al., 2010, 2012a, 2012b；Klotzsche et al., 2010, 2012；Cordua et al., 2012；Belina et al., 2012；吴俊军，2012；吴俊军等，2014)及频率域全波形钻孔雷达数据反演(El Bouajaji et al., 2011；Ellefsen et al., 2011；Yang et al., 2012).
与钻孔GPR相比，地面GPR全波形反演受测量方式制约，反演的不适定性更显著(Meles et al., 2012a)，计算量大，反演效率过低、易受初始模型约束，成像效果不甚理想.为了进一步提高地面GPR全波形反演精度，国内外众多学者开展了深入研究，取得了一系列成果：Moghaddam等(1991)开展小目标体介电常数波形反演成像研究；Cai等(1996)利用二维Maxwell方程进行电磁波速度反演，其结果优于射线追踪成像，但计算量增加若干数量级；Liu等(2004)应用了小波Galerkin算法解积分方程进行GPR成像；Cui等(2004)利用高阶Born近似研究了地下目标体反演技术；Rucker和Ferré(2005)用模拟退火法进行GPR资料的非线性反演，但算法的效率过低；Lavoué等(2014，2015)采用频率域的拟牛顿法对多偏移距GPR数据进行了二维介电常数和电导率的全波形反演.国内，毛立峰(2007)利用改进的非线性共轭梯度反演方法和改进的有限内存拟牛顿约束反演方法，对超宽带雷达数据进行了双参数反演；王兆磊等(2007)用共轭梯度法迭代反演得到了粗糙地表情况下模型的介电常数分布；成艳和张建中(2008)应用GPR线性逆散射进行了地下目标体重建；李壮等(2009)将大范围收敛的同伦方法引入到地下介质参数识别的反演过程中，构造了基于同伦反演的系列算法；丁亮等(2012)将小波多尺度方法结合正则-高斯-牛顿迭代法，得到了Maxwell方程反演的最优解；黄忠来和张建中(2013)提出了一种利用GPR频谱属性同时估算地下层状介质的几何参数与电性参数的全局分步优化反演方法，提高了反演效率和反演结果的可靠性；周辉等(2014)提出了一种不需提取激发脉冲的GPR全波形反演方法，从而避免从实际雷达资料提取激发脉冲的难题，给出了一种GPR反演新思路；李静(2014)采用有限带宽阻抗反演和Monte Carlo反演实现地面雷达和井中雷达随机介质目标反演和参数估计；王洪华(2015)实现了层状介质的GPR阻尼最小二乘法全波形反演，算法能够有效用于实测数据的反演；孟旭(2016)对时间域和频率的钻孔雷达波形反演进行了比较研究.
上述文献中介绍的雷达全波形反演，大多属于介电常数或电导率的单参数反演，在反映目标体的物性参数方面具有一定的局限性.为了使雷达反演剖面具有更好的分辨率，进一步提高GPR资料解释准确度，本文提出了一种时间域全波形优化共轭梯度双参数反演方法.该方法基于地面GPR数据，通过采用GPU并行计算及维度提升策略，大幅提高反演效率.由于该方法同时反演介电常数与电导率双参数，有助于更充分地挖掘雷达剖面中蕴含的信息，达到联合评估、多参数约束、相互验证，提高GPR解释精度的最终目的.
1 时间域全波形优化共轭梯度反演 1.1 正演问题的时域有限差分法由电磁场与电磁波的理论可知(Taflove and Brodwin, 1975)，二维TM(横磁)模式(Hx, Hz, Ey)GPR波满足的Maxwell方程可表示为：
|
(1) |
其中，E为电场强度(V·m-1)，H为磁场强度(A·m-1)，ε为介电常数(F·m-1)，μ为磁导率(H·m-1)，σ为电导率(S·m-1).
由于制约目前雷达剖面法反演的主要因素为反演速度，综合考虑计算速度与并行性加速的便利性，本文采用交错网格的时域有限差分法(FDTD)求解方程(1)，其离散方式与计算步骤如图 1所示，截断边界处采用非分裂的卷积完全匹配层(Roden and Gedney, 2000)对非物理的反射波进行吸收.
|
图 1 FDTD离散方式与更新策略 Fig. 1 The discrete mode and update strategy of FDTD |
如果将式(1)改写为G(m, u)=0的行列式形式，重写方程(1)，则可以得到
|
(2) |
其中L=A∂x+B∂z+C∂t+D为正演算子，u=(HxHzEy)T为雷达波场向量，系数矩阵A，B，C，D为：
|
(3) |
探地雷达全波形反演实质是利用已知的实测数据来重构介质中的模型参数，如求取合理的相对介电常数εr与电导率σ的空间分布，使得正演模拟数据与实测数据之间的拟合最优，其目标函数可以定义为：
|
(4) |
式(4)中N为源的个数.m=(εr, σ)T为模型参数，H0为观测空间的接收区域；dn(x, z, t)为模拟数据，dnobs(x, z, t)为实测数据.如果以算子R表述观测空间中接收器的位置，则模拟数据可以表示为dn(x, z, t)=Ru.
为此，雷达全波形反演求解问题可转化为满足方程(2)的约束条件下，确定方程(4)目标函数的最小值点.利用Lagrange乘数法，该反演问题可转化为无约束优化问题，
|
(5) |
式(5)中，S为新的目标函数，H表示所研究的空间区域(H：0≤x≤X, 0≤z≤Z)，λ为Lagrange乘数.式(5)的解是关于u和λ的Lagrange函数的稳定点，根据新目标函数S关于u的稳定点即可导出伴随方程，因此λ也称伴随波场(刘文革和贺振华，2015)，表示为λ=(Hx*, Hz*, Ey*)T.设伴随场λ满足最终条件
|
(6) |
将式(2)代入式(5)第二项可得
|
(7) |
对式(7)右端项x，z及时间项t分别应用分部积分公式可得：
|
(8) |
根据边值条件u|x=0=0，u|x=X=0，u|z=0=0，u|z=Z=0及初值条件u|t=0=0, 并结合式(6), 式(8)可以简化为：
|
(9) |
其中L*=-A∂x-B∂z-C∂t+D为算子L的伴随算子.
将式(4)与式(9)代入式(5)可得
|
(10) |
然后计算∂S/∂u=0，可以得到第n个的源伴随方程：
|
(11) |
分析式(11)可知，伴随方程也是一种波动方程，以残差场为震源在接收位置进行激发，时间上反向传播.基于一阶方程推导的伴随方程与原始方程(2)的系数矩阵D存在差别，它在原始方程的基础上进行了衰减的补偿.
采用共轭梯度法进行雷达全波形反演，需要进行待反演参数的梯度求取运算，为此，对式(5)求导可得：
|
(12) |
将系数矩阵式(3)的具体参数代入式(12)，可得到目标函数对相对介电常数εr与电导率σ导数的具体表达式：
|
(13) |
|
(14) |
式(13)与式(14)中的波场值Ey及伴随波场值Ey*可分别通过对式(1)及式(11)求解得到.由此可见，共轭梯度法的梯度可以通过残差的伴随波场与激发源正传波场的积分获取，只用到了目标函数及其梯度值，仅需两次FDTD正演，就可以获得目标函数的梯度.在导出目标函数的梯度之后，即可直接应用梯度类寻优算法进行反演.该方法避免了二阶导数(Hessian矩阵)的计算，从而降低了计算量和存储量.
为了使共轭梯度法的迭代结果更好地收敛，常在每一迭代步中，利用当前点处的最速下降方向来生成关于凸二次函数的共轭方向.根据不同的共轭梯度方向修正因子β的构造方法，可得到不同的梯度公式.本文采用Fletcher-Reeves公式(Fletcher and Reeves, 1964)，其算法流程如下：
(1) 给定迭代精度tol和初始点模型m0，令迭代次数k=0，计算目标函数梯度g0；
(2) 如果‖gk‖≤tol，停止，输出m=mk；
(3) 计算搜索方向dk，修正因子β
|
(15) |
|
(16) |
(4) 利用非精确搜索方法确定搜索步长αk；
(5) 更新模型，mk+1=mk+αkdk，通过式(13)与(14)计算梯度gk+1；
(6) 令k=k+1，转到第(2)步，进行循环迭代.
迭代步长是全波形反演中非常关键的参数.如果步长过小，则反演收敛较慢，反演耗时较长；如果步长太大，易导致反演不收敛.而步长的选取方法很多，本文采取基于Wolfe准则的非精确线搜索方法，它需要满足以下两个条件：
|
(17) |
|
(18) |
式(17)与(18)中0＜c1＜c2＜1.设定的初始试探步长选取为：
|
(19) |
其中mean()表示取平均.当进行双参数反演时，将两个参数的梯度乘以不同的因子(其中介电常数的因子为1，电导率的因子由两个参数梯度的量级确定)，然后将两者的梯度合并成一个向量，采取Wolfe约束条件进行步长选取，两个参数的梯度乘以不同因子，因此两个参数梯度的步长是不同的.可以保证每次迭代过程中同时更新两个参数，不增加寻找步长的计算量.
1.3 梯度优化策略在导出目标函数的梯度之后，即可应用常规共轭梯度法进行雷达全波形反演，但是，由于电磁波传播过程中存在几何扩散效应、聚散现象，导致反演迭代过程中，距脉冲源点或接收点较远的区域，其模型参数的反演精度亟待提高，难以较好地重构.为此，本文引入梯度优化方法，采用正演波场和伴随场的能量来计算预条件算子，从而解决距脉冲源点或接收点较远的区域模型的重构和收敛效率问题.该优化预处理的能量梯度公式如下所示：
|
(20) |
其中
|
(21) |
Ws、Wr分别为源点正演波场的预条件算子以及对接收点数据残差的伴随场的预条件算子.gw为优化预处理之后的能量梯度.通过式(20)优化后，可以有效地弱化电磁波传播过程中的几何扩散效应、聚散现象对反演的影响，得到更精确的反演模型参数，并加快收敛效率.
综上所述，应用优化共轭梯度法全波形雷达反演的流程如图 2所示.
|
图 2 GPR全波形反演计算流程 Fig. 2 Flow chart of GPR full waveform inversion |
探地雷达全波形反演计算量大，对计算机要求较高，计算过程中需要多次调用雷达正演计算，过去常局限于在大型计算机、工作站和个人计算机集群上计算，难以在个人微型计算机上实现.因此，怎样提高GPR正演计算效率和降低内存占用，提高计算速度成为关注焦点.为此，作者提出将GPU并行计算与维度提升策略结合起来，以提高反演速度，使雷达全波形反演在微机上计算成为可能.
2.1 基于GPU并行加速正演近年来，随着现代图形处理器单元(GPU)的通用计算能力以及可编程性的提高，基于GPU的并行得到了广大地球物理研究人员的青睐，可以极大地提高微机计算速度.GPU和CPU设备的架构是不同的，GPU适宜计算大量比较简单的任务，CPU的设计是用来运行少量比较复杂的任务.CPU是典型的双核或四核设备，而GPU则拥有多个流多处理器(SM)，每个SM可看作CPU的一个核.
二维GPR波场的单源正演过程中，首先给定波场初始值，通过电场和磁场之间的交替更新计算，将电场和磁场的分量值保存在二维矩阵中，使用矩阵索引来导出电场和磁场的分量值，从而模拟电磁波波场随时间的传播情况.使用GPU并行的二维FDTD的算法流程如图 3所示，图中可见，首先由CPU将模拟所需的场值变量数据初始化，分配给GPU板上的全局内存，用全局内存访问；初始化之后，每个场值分量的计算模块则通过调用GPU内核(Kernel)程序实现，在Kernel中所有节点的计算采用并行方式执行；反复循环该过程，到给定的时间为止.
|
图 3 GPU-FDTD算法 Fig. 3 GPU-FDTD algorithm |
为了更加直观明了地对比GPU并行算法与传统的CPU串行算法特点.本文选取一个线源位于计算域中心的二维简单模型(不含吸收边界)，对比了GPU并行与CPU串行条件下，线源电磁波传播的FDTD算法模拟效率.测试计算在NVIDIA GeForce GTX 960 GPU和Inter(R) Core(TM) i7-6700K CPU上.表 1显示了CPU和GPU并行300个时间步耗时比较.图 4a为CPU和GPU并行计算时间对比图，图中纵坐标采用对数坐标，图 4b为GPU并行对CPU串行计算的加速比图.
|
表 1 CPU、GPU的计算时间(300步) Table 1 Computation time for the CPU only and that of the GPU (300 steps) |
|
图 4 GPU并行与CPU串行计算效率对比 (a) CPU、GPU的计算时间(300步); (b)加速比. Fig. 4 Comparison of GPU parallel and CPU serial computing efficiency (a) Computation time for the CPU only and that of the GPU (300 steps); (b) Speedup factor. |
分析图 4和表 1可知，对于256及其以下小网格规模，CPU实际上更快，此时数据传输时间占更大比重，降低了程序的整体性能.当网格规模大于512时，GPU并行算法相对CPU串行计算有较大的加速比，尤其是当网格规模达到2560时，计算时间从CPU上的2 min降低到GPU的10 s，加速比达到了11.当网格规模为3072及4096时，由于数据所占空间接近GPU内存限制，加速比降低.因此，对于大规模的二维/三维GPR正演采用GPU并行加速计算，能极大地提高模拟效率.
2.2 维度提升策略无论是实际的雷达探测还是正演模拟过程中，由于探地雷达剖面法的工作方式，每采集一道雷达数据，发射天线(脉冲源)移动一固定距离，属于典型的多源问题.若仅在FDTD正演计算时单纯的依靠GPU并行加速计算，尽管可取得较高的加速比，但全波形反演仍需较长的时间.为了进一步加快反演速度，考虑到每个雷达正演剖面需要多次计算移动的脉冲雷达波源，而改变模型参数后又要进行同样的重复性工作，为此作者提出维度提升策略，将传统的多源二维问题提升为图 5所示的三维问题，而每个雷达剖面中多次移动的脉冲雷达波场变为了在三维剖面中的有限的1次或2次激发的三维问题，且每个二维切片之间不进行数据交换，场值更新过程通过GPU-FDTD算法实现，每个三维问题的网格规模均位于高加速比区间，从而实现GPU加速的快速正演.
|
图 5 二维多源GPR正演问题维度提升 Fig. 5 2D multi-source GPR forward dimension promotion |
表 2列出了三种不同算法：CPU串行计算(CPU1)、CPU四核并行计算(CPU2)、维度提升的GPU-FDTD算法(GPU)，模拟不同尺寸离散模型(采用CPML吸收边界，PML层数为10)在50个源，500个时间步所花费的时间.由表 2中可以发现，在测试的几个不同尺寸模型中，GPU算法的耗时相对CPU算法明显减少，加速比大于4.5，特别是420×420的网格，加速比达到了8.因此对于一定规模的多源二维GPR正演问题，维度提升的GPU-FDTD算法可以实现快速正演.
|
|
表 2 三种算法计算时间(50个源，500时间步) Table 2 Computation time for three algorithms (50 sources, 500 steps) |
根据上述方法原理编制了时间域GPR全波形反演程序，首先以图 6a所示的简单地电模型为例，该模型的设置主要为了比较多偏移距与自激自收两种观测方式对反演精度的影响.进行相对介电常数的全波形反演.模拟区域为1.0 m×1.0 m，模型背景相对介电常数5.0，模型区域内设置有大小相等的两个长方形异常体，异常体的尺寸为0.2 m×0.05 m，顶界面埋深为0.3 m，其中左侧异常体的水平中心位置为0.3 m，相对介电常数为1，右侧异常体的水平中心位置为0.7 m，相对介电常数为10，模型区域电导率为0.
|
图 6 (a) 真实模型；(b)初始模型；(c)自激自收反演结果；(d)多偏移距反演结果 Fig. 6 (a) Real model; (b) Initial model; (c) Self-reflection inversion results; (d) Multi-offset inversion results |
利用FDTD算法进行该模型的正演及伴随场计算时，采用的离散步长为0.01 m，网格数为100×100.边界条件为CPML，该完全匹配层设为10层，厚度为0.1 m.分别采用自激自收和多偏移距两种观测方式进行模拟，激励源为500 MHz的雷克子波，采样时间间隔为0.02 ns，时窗长度为10 ns.在地表布设51个反射天线，发射天线的间距为0.02 m，多偏移距观测方式时布设101个接收天线，接收点距为0.01 m，水平方向上覆盖1 m.
图 6为两种接收方式下，迭代50次的反演结果，其中自激自收耗时1902.8 s，多偏移距模式耗时1211.4 s.两种模式下，均能得到较好的反演结果，其中多偏移距由于其接收的数据量要大于自激自收观测方式，且含有较大角度的反射信息，反问题的不适定性相对较小，因此反演效果较好.由于激发点与接收点位于异常体的一侧，这两种观测方式对异常体垂向的分辨效果更好.图 7为两种模式下的目标函数收敛曲线以及异常体位置的介电常数切线图，在相同的反演迭代次数下，相对于自激自收方式，多偏移距观测方式的反演结果与真实模型更为接近，具有更好的收敛性、稳定性.
|
图 7 (a) 目标函数；(b)深度0.33 m介电常数切线；(c)水平0.3 m介电常数切线；(d)水平0.7 m介电常数切线 Fig. 7 (a) Objective function; (b) The permittivity in depth of 0.33 m; (c) The permittivity in distance of 0.3 m; (d) The permittivity in distance of 0.7 m |
起伏界面模型如图 8a所示，该模型的设置主要为了对比优化前后的共轭梯度算法的收敛性、计算效率.模型大小为1.0 m×1.0 m，埋深0.3 m处存在一条起伏界面，起伏界面上方介质的相对介电常数为6.0；下方介质的相对介电常数为5.0.起伏界面下方设置两个直径为0.08 m的圆形异常体，左侧异常体的相对介电常数为10，其圆心位置为(0.25 m，0.50 m)；右侧空洞异常体的相对介电常数为1.0，其圆心位置为(0.75 m，0.50 m).该模型的电导率为0.001 S·m-1.
|
图 8 (a) 真实模型；(b)初始模型；(c)原始梯度g；(d)优化梯度gw Fig. 8 (a) Real model; (b) Initial model; (c) The original gradient g; (d) Optimization gradient gw |
采用的离散步长为0.01 m，网格数100×100，PML边界为10层，厚度为0.1 m.采用多偏移距观测方式对该模型进行模拟，发射源信号采用中心频率为500 MHz的雷克子波，采样时间间隔为0.02 ns，时窗长度为16 ns.在地表布设51个反射天线，发射天线的间距为0.02 m，布设101个接收天线，接收点距为0.01 m，水平方向上覆盖1 m.
图 8为不同梯度公式，相同初始模型，迭代100次的反演结果，原始梯度公式耗时4576.4 s，优化梯度公式耗时4729.0 s；两者均能得到较好的反演结果，起伏界面与真实模型拟合基本一致，两个圆形异常体的位置与真实模型相同，高相对介电常数异常体的形态拟合较好，低相对介电常数异常体的形态有一定的差异.图 9为反演的目标函数收敛曲线以及异常体相应位置的介电常数切线图，图 9a中明显可见，优化之后的梯度相对于原始梯度公式可以在一定程度上提高收敛速度.因此，尽管在同样迭代100次前提下，优化的共轭梯度法耗时更多，但是由于该优化方法具有更好的收敛性，因此，针对同样的反演问题，优化的共轭梯度法可以采用更少迭代次数收敛到需要的反演精度.
|
图 9 (a) 目标函数；(b)深度0.5 m介电常数切线；(c)水平0.25 m介电常数切线；(d)水平0.75 m介电常数切线 Fig. 9 (a) Objective function; (b) The permittivity in depth of 0.5 m; (c) The permittivity in distance of 0.25 m; (d) The permittivity in distance of 0.75 m |
复杂界面模型大小为1.0 m×1.0 m，模型相对介电常数和电导率分布如图 10所示.计算区域离散间隔为0.01 m，网格数100×100.PML边界为10层，厚度为0.1 m.采用多偏移距观测方式对该模型进行模拟，雷克子波的中心频率分别为400 MHz、500 MHz和600 MHz.采样时间间隔为0.0333 ns，时窗长度为20 ns.在地表布设51个反射天线，发射天线的间距为0.02 m，布设101个接收天线，接收点距为0.01 m，水平方向上覆盖1 m.
|
图 10 (a) 介电常数模型；(b)电导率模型 Fig. 10 (a) Permittivity model; (b) Conductivity model |
对复杂界面模型进行不同天线频率下的介电常数与电导率的双参数全波形反演，反演迭代次数均为200次.图 11和图 12分别为不同天线主频下，相同相对介电常数与电导率初始模型，介电常数与电导率的反演结果.图 13为反演的目标函数收敛曲线以及异常体相应位置的介电常数切线图，图 14为异常体相应位置的电导率切线图.
|
图 11 (a) 介电常数初始模型；(b)主频400 MHz介电常数反演结果；(c)主频500 MHz介电常数反演结果；(d)主频600 MHz介电常数反演结果 Fig. 11 (a) Initial model of permittivity; (b) Inversion results of permittivity with 400 MHz; (c) Inversion results of permittivity with 500 MHz; (d) Inversion results of permittivity with 600 MHz |
|
图 12 (a) 电导率初始模型；(b)主频400 MHz电导率反演结果；(c)主频500 MHz电导率反演结果；(d)主频600 MHz电导率反演结果 Fig. 12 (a) Initial model of conductivity; (b) Inversion results of conductivity with 400 MHz; (c) Inversion results of conductivity with 500 MHz; (d) Inversion results of conductivity with 600 MHz |
|
图 13 (a) 目标函数；(b)深度0.5 m介电常数切线；(c)水平0.35 m介电常数切线；(d)水平0.6 m介电常数切线 Fig. 13 (a) Objective function; (b) The permittivity in depth of 0.5 m; (c) The permittivity in distance of 0.35 m; (d) The permittivity in distance of 0.6 m |
|
图 14 (a) 深度0.5 m电导率切线；(b)水平0.35 m电导率切线；(c)水平0.6 m电导率切线 Fig. 14 (a) The conductivity in depth of 0.5 m; (b) The conductivity in distance of 0.35 m; (c) The conductivity in distance of 0.6 m |
由图 11中可以发现，对于介电常数反演结果而言，三个主频均能得到较好的反演结果，起伏界面与真实模型拟合基本一致，两个不规则的异常体的位置与真实模型相同，相对介电常数异常体的形态拟合较好，仅幅值与真实模型有一定的差异.由于真实模型的电导率分布为1~8 mS/m，差异较小，同时电导率对于反射系数的影响较小，由电导率差异引起的异常表现在多偏移距接收数据的比重较小.因此在图 12中的电导率反演结果与真实模型的趋势大体相同，但是两个异常体的形态与幅值较真实模型而言，仍有较大的差距.
图 13a中的目标函数曲线，随着频率的增大收敛速度变缓；由介电常数和电导率的切线图也可以发现，400 MHz的反演结果整体上更为接近真实模型，其反演结果的振荡也相对较小.由梯度的求取公式可知，梯度是由正传播场与残差的反传波场的乘积求和得到，因此雷达波长对于梯度的形态有重要的影响，当异常体尺寸与雷达波长大致相当的条件下，所得的梯度方向可以更快地收敛.在该模型参数下，400 MHz频率的波长能更好地拟合异常形态，因此会有更快的收敛速度.
4 结论(1) 将优化的非线性共轭梯度法应用于时间域全波形探地雷达反演，避免了Hessian矩阵的直接求解，有效降低了计算量和存储量.且反演过程中，每个电场、磁场分量节点的计算，以模块方式调用GPU内核程序采用并行方式执行，当网格数目较大的情况下，相对CPU串行计算有较大的加速比.结合维度提升策略，极大地提高微机计算速度，在普通微机上实现了时间域全波形二维GPR快速反演.
(2) 开展了3个实例模型的反演实验.第一个含高低、阻异常体的规则模型反演结果表明，多偏移距较自激自收的接收方式，含有较大角度的反射信息，信息量更大，FWI的反演效果更好；第二个起伏界面模型，对比了优化共轭梯度前后反演算法的收敛效果和收敛速度，说明优化后的共轭梯度法弱化了源点、接收点的几何扩散效应、聚散现象，具有更高的收敛速度；最后一个复杂起伏界面模型开展了不同主频下的时间域全波形双参数反演说明，当异常体尺寸与雷达波长大致相当的条件下，具有较好的反演收敛速度及反演结果，而双参数反演能给单一的介电常数反演提供更丰富的信息约束，介电常数与电导率反演结果通过相互佐证，互为补充，有助于更为准确地判定原模型.
(3) 由同构模型双参数反演实验表明：由于反演当中不同参数对数据的敏感性不同，导致介电常数与电导率反演精度不一致，双参数反演时的电导率的反演结果稍差；若模型中两个参数的界面存在不一致时，两个参数之间的串扰噪声会更加突出，需要对模型中不同参数进行敏感性分析，并根据该分析结果对反演参数进行尺度变换，以进一步改善同步多参数的反演效果，提高反演精度.
Belin a F A, Irving J, Ernst J R, et al. 2012. Waveform inversion of crosshole georadar data:influence of source wavelet variabilityand the suitability of a single wavelet assumption. IEEE Transactions on Geoscience and Remote Sensing, 50(11): 4610-4625. DOI:10.1109/TGRS.2012.2194154 |
Cai W Y, Qin F H, Schuster G T. 1996. Electromagnetic velocity inversion using 2-D Maxwell's equations. Geophysics, 61(4): 1007-1021. DOI:10.1190/1.1444023 |
Cheng Y, Zhang J Z. 2008. A linear object reconstruction algorithm for ground penetrating radar. Journal of Xiamen University (Natural Science) (in Chinese), 47(1): 31-35. |
Cordua K S, Hansen T M, Mosegaard K. 2012. Monte Carlo full-waveform inversion of crosshole GPR data using multiple-pointgeostatistical a priori information. Geophysics, 77(2): H19-H31. DOI:10.1190/geo2011-0170.1 |
Cui T J, Qin Y, Wang G L, et al. 2004. Low-frequency detection of two-dimensional buried objects using high-order extended Born approximations. Inverse Problems, 20(6): S41-S62. DOI:10.1088/0266-5611/20/6/S04 |
Ding L, Han B, Liu R Z, et al. 2012. Inversion imaging method forconcrete non-destructive testing based on GPR. Chinese Journal of Geophysics (in Chinese), 55(1): 317-326. DOI:10.6038/j.issn.0001-5733.2012.01.032 |
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. DOI:10.1190/1.3554412 |
Ernst J R, Maurer H, Green A G, et al. 2007. Full-waveform inversion of crosshole radar data based on 2-D finite-differencetime-domain solutions of Maxwell's equations. IEEE Transactions on Geoscience and Remote Sensing, 45(9): 2807-2828. DOI:10.1109/TGRS.2007.901048 |
Fletcher R, Reeves C M. 1964. Function minimization by conjugate gradients. The Computer Journal, 7(2): 149-154. DOI:10.1093/comjnl/7.2.149 |
Gloaguen E, Marcotte D, Chouteau M, et al. 2005. Borehole radar velocity inversion using cokriging and cosimulation. Journal of Applied Geophysics, 57(4): 242-259. DOI:10.1016/j.jappgeo.2005.01.001 |
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 |
Hu G H, Wang L X, Fang W B, et al. 2014. Full WaveformInversion Method and Application (in Chinese). Beijing: Petroleum Industry Press.
|
Huang Z L, Zhang J Z. 2013. An inversion method for geometric and electric parameters of layered media using spectrum of GPR signal. Chinese Journal of Geophysics (in Chinese), 56(4): 1381-1391. DOI:10.6038/cjg20130432 |
Klotzsche A, Van Der Kruk J, Meles G A, et al. 2010. Full-waveform inversion of crosshole ground penetrating radar data to characterize a gravel aquifer close to the Thur River, Switzerland.//Proceedings of the XⅢ Internarional Conference on Ground Penetrating Radar. Lecce, Italy: IEEE. http://ieeexplore.ieee.org/document/5550238
|
Klotzsche A, Van Der Kruk J, Meles G, et al. 2012. CrossholeGPR full-waveform inversion of waveguides acting as preferential flow paths within aquifer systems. Geophysics, 77(4): H57-H62. DOI:10.1190/geo2011-0458.1 |
Lailly P. 1983. The seismic inverse problem as a sequence of before stack migrations.//Bednar J B ed. Conference on Inverse Scattering: Theory and Applications. Philadelphia: Society for Industrial and Applied Mathematics, 206-220. http://citeseerx.ist.psu.edu/showciting?cid=284398
|
Lavoué F, Brossier R, Métivier L, et al. 2014. Two-dimensionalpermittivity and conductivity imaging by full waveform inversion of multioffset GPR data:A frequency-domain quasi-Newton approach. Geophysical Journal International, 197(1): 248-268. DOI:10.1093/gji/ggt528 |
Lavoué F, Brossier R, Métivier L, et al. 2015. Frequency-domainmodelling and inversion of electromagnetic data for 2D permittivity andconductivity imaging:An application to the Institut Fresnelexperimental dataset. Near Surface Geophysics, 13(3): 227-241. |
Li D X. 1994. Method and Application of Ground Penetrating Radar (in Chinese). Beijing: Geological Publishing House.
|
Li J. 2014. Ground penetrating radar detection and parameter inversion in stochastic effective medium[Ph. D. thesis] (in Chinese). Changchun: Jilin University.
|
Li Z, Han B, Chen Y. 2009. Wavelet-homotopy hybrid algorithmfor inverse problem of ground penetrating radar. Chinese Journal of Computational Physics (in Chinese), 26(1): 87-93. |
Liu F S, Ling Y, Xia X G, et al. 2004. Wavelet methods for ground penetrating radar imaging. Journal of Computational and Applied Mathematics, 169(2): 459-474. DOI:10.1016/j.cam.2003.12.036 |
Liu W G, He Z H. 2015. Seismic Waveform Inversion (in Chinese). Beijing: Science Press.
|
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. 2012b. 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 |
Meles G A, Greenhalgh S A, Green A G, et al. 2012a. GPR full-waveform sensitivity and resolution analysis using an FDTD adjoint method. IEEE Transactions on Geoscience and Remote Sensing, 50(5): 1881-1896. DOI:10.1109/TGRS.2011.2170078 |
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.
|
Moghaddam M, Chew W C, Oristaglio M. 1991. Comparison of theBorn iterative method and Tarantola's method for an electromagnetic time-domain inverse problem. International Journal of Imaging Systems and Technology, 3(4): 318-333. DOI:10.1002/(ISSN)1098-1098 |
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/(ISSN)1098-2760 |
Rucker D F, Ferré T P A. 2005. Automated water content reconstruction of zero-offset borehole ground penetrating radar data using simulated annealing. Journal of Hydrology, 309(1-4): 1-16. DOI:10.1016/j.jhydrol.2004.11.008 |
Taflove A, Brodwin M E. 1975. Numerical solution of steady-state Electromagnetic scattering Problems using the time-dependentMaxwell'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 H H. 2015. Finite element forward modeling of GPR in special media and its interpretation techniques[Ph. D. thesis] (in Chinese). Changsha: Central South University.
|
Wang Z L, Zhou H, Li G F. 2007. Inversion of ground-penetrating radar data for 2D electric parameters. Chinese Journal of Geophysics (in Chinese), 50(3): 897-904. |
Wu J J. 2012. Research on cross-borehole radar tomography using full-waveform inversion method[Ph. D. thesis] (in Chinese). Changchun: Jilin 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 |
Yang X, Van Der Kruk J, Bikowski J, et al. 2012. Full-waveform inversion of GPR data in frequency-domain.//2012 14th International Conference on Ground Penetrating Radar (GPR). Shanghai, China: IEEE, 324-328. http://www.researchgate.net/publication/259740666_Full-waveform_inversion_of_GPR_data_in_frequency-domain
|
Zhou H, Sato M, Takenaka T, et al. 2007. Reconstruction fromantenna-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 |
成艳, 张建中. 2008. 一种线性探地雷达目标重建算法. 厦门大学学报(自然科学版), 47(1): 31-35. DOI:10.3321/j.issn:0438-0479.2008.01.007 |
丁亮, 韩波, 刘润泽, 等. 2012. 基于探地雷达的混凝土无损检测反演成像方法. 地球物理学报, 55(1): 317-326. DOI:10.6038/j.issn.0001-5733.2012.01.032 |
胡光辉, 王立歆, 方伍宝, 等. 2014. 全波形反演方法及应用. 北京: 石油工业出版社.
|
黄忠来, 张建中. 2013. 利用探地雷达频谱反演层状介质几何与电性参数. 地球物理学报, 56(4): 1381-1391. DOI:10.6038/cjg20130432 |
李大心. 1994. 探地雷达方法与应用. 北京: 地质出版社.
|
李静. 2014.随机等效介质探地雷达探测技术和参数反演[博士论文].长春: 吉林大学. http://cdmd.cnki.com.cn/Article/CDMD-10183-1014268027.htm
|
刘文革, 贺振华. 2015. 地震波形反演. 北京: 科学出版社.
|
李壮, 韩波, 陈勇. 2009. 探地雷达反问题的小波-同伦混合反演算法. 计算物理, 26(1): 87-93. DOI:10.3969/j.issn.1001-246X.2009.01.012 |
毛立峰. 2007.超宽带电磁法正演模拟及反演成像[博士论文].成都: 成都理工大学. http://cdmd.cnki.com.cn/Article/CDMD-10616-2007138668.htm
|
孟旭. 2016.时间域和频率域跨孔雷达波形反演及比较研究[博士论文].长春: 吉林大学. http://cdmd.cnki.com.cn/Article/CDMD-10183-1016090168.htm
|
王洪华. 2015.特殊介质的GPR有限元正演及解释技术研究[博士论文].长沙: 中南大学.
|
王兆磊, 周辉, 李国发. 2007. 用地质雷达数据资料反演二维地下介质的方法. 地球物理学报, 50(3): 897-904. DOI:10.3321/j.issn:0001-5733.2007.03.031 |
吴俊军. 2012.跨孔雷达全波形层析成像反演方法的研究[博士论文].长春: 吉林大学. http://cdmd.cnki.com.cn/Article/CDMD-10183-1012365348.htm
|
吴俊军, 刘四新, 李彦鹏, 等. 2014. 跨孔雷达全波形反演成像方法的研究. 地球物理学报, 57(5): 1623-1635. DOI:10.6038/cjg20140525 |
周辉, 陈汉明, 李卿卿, 等. 2014. 不需提取激发脉冲的探地雷达波形反演方法. 地球物理学报, 57(6): 1968-1976. DOI:10.6038/cjg20140627 |