探地雷达(Ground Penetrating Radar, GPR)根据地下介质的电磁性质差异来查找目标异常，是一种高分辨率、无损的浅部地球物理探测方法，在工程质量检测^[1]、公路质量检测^[2]、地质灾害研究^[3]等领域均有广泛应用。

在道路病害研究中，识别并提取GPR剖面中的道路病害特征一直是GPR能否成功应用的关键。在实际工作中，往往存在不同程度的环境干扰，导致检测精度低、测量效果变差。对此，差值检测是一种有效提高数据质量的方法^[4]。建立地下病害属性划分标准，则主要依据回波波组形态、振幅和相位特性、吸收衰减特性等。在理论研究方面，主要通过正演研究电磁信号的传播特性，模拟得到典型病害体的特征图谱来对实际情况进行指导^[5-6]；在工程上，研究探测的重点区域及其测线布置原则，采集参数设置及数据采集质量的评判方法。同时，对正常道路、典型干扰源和典型道路隐性病害进行标记，并对应其地球物理特征和GPR波组特征，可以很好地对相似特征下的病害进行解释，也可为路面养护维修提供决策支持^[7-9]。除病害外，道路中沥青层覆盖的混凝土裂缝也是GPR应用于道路无损检测的重点方向^[10]。同时，对道路采集得到的GPR数据进行特定的处理，能够有效突出病害特征，再将正演或标记的病害特征与处理后的数据进行综合对比分析^[11]，提高路基中隐藏病害的探明精度和后续解释的准确度^[12]。

反演是一种对所研究区域进行成像分析的有效手段，可以更加直观地展示病害类型，为病害解释提供理论依据。在众多GPR反演技术中，基于非线性理论基础的全波形反演(Full Waveform Inversion，FWI)是目前最常用的一种方法。该方法有效利用了全波场信息，通过波动方程约束条件来解决非线性优化问题，经过多次反演迭代得到贴合地下实际情况的物性参数分布，如几何信息特征和介质的介电特性^[13]；在FWI的应用中，激励源子波对反演的结果具有很大影响，使用不同的源子波进行FWI，可以验证不同饱和度的土壤–含水层系统^[14]；稀疏盲反卷积技术是一种优化子波的可靠技术手段，建立地下反射率序列的稀疏表示，提高了反演结果的准确率^[15]；通过分析测量和模型轨迹的拟合以及最终模型的剩余梯度来对独立测量的测井实测数据进行验证分析得出，子波估计主要通过影响初始速度模型进而对FWI结果造成影响，初始速度模型与目标函数联系紧密^[16]；在实际工程实践中仪器所产生的激励源子波估计往往并不准确，若能在建立目标函数的过程中，消除子波影响，可以有效提高FWI的准确性和实用性^[17-18]。

传统目标函数FWI算法都是基于激励源子波已知的假设条件，然而在实际勘探中非常难以获取准确的激励源子波形态，激励源子波估计不准确会降低反演结果的可靠性和准确性。因此，如何降低或消除子波对反演结果的影响至关重要。本文在前人的基础上，将褶积型目标函数的FWI表达形式应用于GPR反演中，在建立目标函数过程中，消除激励源子波估计不准确的影响，提高反演精度，以期为实际工程应用提供理论依据和数据基础。

1   方法原理

1.1   探地雷达正演

二维横磁(Transverse Magnetic, TM)波Maxwell方程如下：

式中：H为磁场强度，A/m；E为电场强度，V/m；ε为介质的介电常数，F/m；μ为介质的磁导率，H/m；σ为电导率，S/m；t为时间，s；J为激励源；x，y，z为方向。本文选用时域有限差分法(Finite Difference Time Domain, FDTD)算法对GPR数据进行正演模拟^[19]，选择卷积完全匹配层(Convolutional Perfectly Matched Layer, CPML)吸收边界用来控制有限模拟区域^[20]。

1.2   传统目标函数全波形反演

FWI就是建立一个目标函数，通过最优化策略找到某个参数值，使得目标函数最小。对于合成数据，由于子波的形态是已知的，因此，可以假设实际采集到的探地雷达数据为$ {E_{\textit{z}}} $分量，对其建立如下的极小化目标泛函^[21]：

式中：E_z为模拟的探地雷达数据(采用猜测参数的正演模拟数据)；${t_{\rm{T}}}$为记录总时长；上标obs为实际的观测数据(实际测量数据、合成剖面数据)；模拟波动方程必须满足Maxwell方程，采用拉格朗日乘子法将无约束优化问题转化为约束优化问题，目标泛函变为：

式中：$[{\psi _x},{\psi _y},\phi ]$为拉格朗日乘子函数；$[{e_1},{e_2},{e_3}]$为伴随方程，其表达式为：

利用分部积分法，注意到初始条件、终止条件为零，且满足自由边界条件。因此，得到：

式中：δ为变量的增量。

式中：x_r为x方向源点位置。

令式(5)、式(6)、式(7)等于0，得到伴随方程表达式为：

伴随方程式(8)的形式与式(4)不同。最后得到相应的梯度公式为：

由于本文使用拟牛顿算法，因此，不需要直接求取Hessian矩阵。

1.3   褶积型目标函数全波形反演

本文采用建立褶积型目标函数的方法，以消除激励源子波的影响，与传统目标函数不同，褶积型目标泛函为：

式中：下标ref为相应的参考道，使用Green函数表示方法，波动方程可以简写为：

式中：d为波动方程表达式，G与f为拆分的矩阵表达式。

将式(11)代入式(10)，得：

式中：上标cal为相应的计算值。

令${{\boldsymbol{G}}_1} = {{\boldsymbol{G}}^{{{\rm{obs}}}}}*{\boldsymbol{G}}_{{\rm{ref}}}^{{{\rm{cal}}}}$，${{\boldsymbol{G}}_2} = {{\boldsymbol{G}}^{{{\rm{cal}}}}}*{\boldsymbol{G}}_{{\rm{ref}}}^{{{\rm{obs}}}}$， ${\boldsymbol{f}} = {{\boldsymbol{f}}^{{{\rm{cal}}}}}*{{\boldsymbol{f}}^{{{\rm{obs}}}}}$，代入式(12)可得：

从化简后的式(12)中可以分析得到，构成目标函数残差项的2个数据项可以看作由相同激励源信号激发产生的，消除了子波估计不准确对于反演的影响。对式(10)求取梯度，有：

其中：$r = E_{\textit{z}}^{{{\rm{obs}}}}*{E_{\textit{z}}}_{,{\rm{ref}}} - {E_{\textit{z}}}*E_{{\textit{z}},{\rm{ref}}}^{{{\rm{obs}}}}$；m为求导参数，代表$\mu 、\sigma、 \varepsilon $；根据拉格朗日乘子法并忽略高次项，得到此时的伴随方程为：

式中：$ \otimes $为相关运算。此外，得到的梯度与式(9)相同。

2   合成数据实验

2.1   路面塌陷异常褶积型目标函数与传统目标函数FWI对比

为了对比2种不同目标函数对于FWI结果的影响并说明本文算法的正确性，建立如图1所示的地面塌陷模型。塌陷模型分为三层，分别为空气层(ε_r =1)，混凝土层(ε_r =9)及基岩层(ε_r =6)，路面塌陷分别为无填充(ε_r =1)及沙土填充(ε_r =3)，并且在地下具有一个异常空洞体(ε_r =1)。

图  1  路面塌陷异常相对介电常数分布

Figure  1.  Relative dielectric constant distribution in case of road collapse exceptions

下载: 全尺寸图片 幻灯片

实验所用激励源为400 MHz的雷克子波，采用屏蔽天线B-Scan的观测方式，得到如图2所示的合成数据正演剖面。

图  2  路面塌陷异常正演剖面

Figure  2.  Forward section of road collapse exceptions

下载: 全尺寸图片 幻灯片

由图2可以看出，路面的2个塌陷异常反映在直达波上，无填充的塌陷部分在直达波部分有明显的下移，并且不存在反射波，但在角点存在绕射异常；沙土填充层的直达波连续，但是由于反射系数的改变导致直达波能量降低，并且具有很明显的反射回波，且角点也存在绕射波；路面内部的空洞位置出现了一个明显的双曲线反射波。约2 ns处的路面模型内部具有十分明显的分层特征，同时路面内部结构的塌陷更为清晰。然而对于约10 ns处空洞下方的层却出现了同相轴的错断，原因是使用地面观测系统时，由于该处上方具有空洞，导致波传播到层界面的时间与其他位置不同。

为了说明激励源子波估计不准确对于FWI结果的影响，将激励源调整为450 MHz的雷克子波。采用混凝土层介电常数构成的均一介质当作反演的初始模型，传统目标函数FWI的结果如图3a所示。改为褶积型目标函数使用同样的参数对该模型进行反演，得到如图3b所示结果。

图  3  路面塌陷异常传统目标函数和褶积型目标函数FWI结果对比

Figure  3.  Comparison between results from conventional objective function and convolution type objective function

下载: 全尺寸图片 幻灯片

对比图3可知，当激励源子波估计不准确时，使用传统目标函数进行反演得到的结果仅仅可以看出模型的大体趋势，直达波连续无错断，也没有扭曲变形，无法对塌陷和空洞的位置和形态信息进行准确反映，同时，0.6 ns处的分层也并不清晰，且反演的具体介电常数数值是不准确的，极大地影响了对于下地介质分布的人工判读，容易造成错误的解释。而使用褶积型目标函数进行反演得到的结果，准确性得到了很大的提升，2个地面塌陷、空洞和层位信息都得到了较好的恢复，且介电常数数值对应良好，能够从介电常数的分布上清楚地确定该处的异常体类型。

特别需要注意的是图3b中对于路面内部的空洞异常体，反演得到的异常体大小并不对应，并且在空洞异常体上方出现了虚假的高介电常数异常。该现象的原因是地面观测这种观测方式导致了对于深层和靠近模型边界的异常得到的回波信号较少，导致该位置在剖面上的波形并不完整，且反演的多解性也在一定程度上深化了这一现象。若想规避这一现象则需要获取更多的地下介质回波信息的约束，如钻孔雷达、声波结果等。

为了能够量化说明FWI算法的精确性，建立如下式所示的重构误差函数:

式中：p为评价的参数，这里为相对介电常数；p₀为初始模型参数；p_true为真实模型参数；p_k为第k次迭代后的模型参数。计算$\Delta (p)$得到如图4所示的模型重构误差曲线。

图  4  路面塌陷异常FWI重构误差曲线

Figure  4.  Error curve of ground surface collapse model reconstruction

下载: 全尺寸图片 幻灯片

由图4可知，随着反演迭代次数的增加，模型重构误差逐渐减小，最后趋于稳定，同时由于探地雷达观测方式的限制，信息量不足以将地下介质信息完全恢复，因此，重构误差会逐渐稳定且不再下降。

在实际工作中通常还需要考虑算法的效率问题，本文基于Ubuntu 20.04操作系统，使用Python3.8编程语言进行全波形反演实验，CPU配置为Intel Xeon(R) CPU E5-2690 3.00 GHz × 20，反演过程中使用20个进程并行的策略，得到FWI效率对比(表1)。

表  1  合成数据全波形反演效率对比

Table  1.  Comparison of full-waveform inversion efficiency for synthetic data

算法	网格大小	使用道数	记录时长/ns	迭代次数	单次迭代耗时/s	总耗时/s
传统算法	120×300	300	30	20	45.234 3	1 885.363 8
本文算法	120×300	300	30	20	47.234 6	1 269.668 5

下载: 导出CSV 

| 显示表格

由表1可知，由于本文算法需要进行卷积运算，因此，单次迭代耗时略大于传统算法，FWI总耗时与单次迭代耗时并不是线性关系，所以传统目标函数FWI算法的总耗时反而大于本文提出的算法，这是由于错误的激励源子波估计，造成了整体算法的不稳定，需要多次进行迭代步长的选取，而本文算法由于消除了激励源子波估计对于反演的影响，因此，整体算法更为稳定，从而获得了更高的算法效率。

2.2   地下空洞异常褶积型目标函数FWI实验

为了说明基于褶积型目标函数FWI对于地下复杂介质的反演效果和算法的正确性，建立如图5所示的地下模型。该模型分为两层，分别为地表层(ε_r =4)及基岩层ε_r =9)，并且在地表层设置了6个异常体，分别用于模拟地下空洞(ε_r =1)和地下管线(ε_r =6, ε_r =8, ε_r =9, ε_r =10)。

图  5  地下空洞异常相对介电常数分布

Figure  5.  Relative dielectric constant distribution in case of underground void exceptions

下载: 全尺寸图片 幻灯片

采用400 MHz的雷克子波作为激励源，采用屏蔽天线B Scan的观测方式，得到如图6所示的合成数据正演剖面。

图  6  路面塌陷异常正演剖面

Figure  6.  Forward section of road collapse exceptions

下载: 全尺寸图片 幻灯片

由图6可知，模型正演剖面十分复杂，由于异常体的相对介电常数各不相同，因此，剖面中的反射波能量也有所差异，具体表现为：相对介电常数相差越大，得到的反射波能量也越强。同时需要注意的是，由于管线存在垂直方向的重叠，正演剖面中的反射回波双曲线存在耦合。

为了说明激励源子波估计不准确对于FWI结果的影响，采用600 MHz的雷克子波作为反演所使用的激励源，采用地表层介电常数构成的均一介质当作全波形反演的初始模型，使用褶积型目标函数，得到如图7所示的反演结果。

图  7  褶积型目标函数FWI结果

Figure  7.  Convolution type objective function FWI result

下载: 全尺寸图片 幻灯片

对比图5和图7可知，使用褶积型目标函数FWI得到的结果由于消除了激励源子波估计不准确带来的影响，得到了较好的反演结果，较为明显地区分了地层，6个异常体位置及数值也得到了较好的对应。同时需要注意的是图7中反演得到的异常体大小并不对应，同时空洞异常体上方出现了虚假的高介电常数异常，造成这个现象的原因是观测系统采用地面观测，这种观测方式导致了对于深层和靠近模型边界的异常得到的回波信号较少，导致该位置在剖面上的波形并不完整，并且由于全波形算法本身存在的多解性，从而导致了这一现象。同时，对于距离较近的不同相对介电常数的异常体也能够通过反演算法进行良好的对应，进而说明了基于褶积型目标函数的FWI算法对于复杂地下介质的正确性。

为了能够量化说明算法的精确性和效率，使用式(16)所示的重构误差函数，得到如图8所示的模型重构误差曲线，并且通过记录得到反演效率(表2)。

图  8  地下复杂介质模型重构误差曲线

Figure  8.  Error curve of reconstruction for underground complex medium models

下载: 全尺寸图片 幻灯片

表  2  合成数据全波形反演效率

Table  2.  Full-waveform inversion efficiency for synthetic data

网格大小	使用道数	记录时长/ns	迭代次数	单次迭代耗时/s	总耗时/s
200×400	200	40	20	156.834 0	4 186.338 7

下载: 导出CSV 

| 显示表格

由图8可知，随着反演迭代次数的增加，模型重构误差逐渐减小，且下降幅度逐渐变缓，同时由于探地雷达观测方式的限制，信息量不足以将地下介质信息完全恢复。

由表2可知，由于使用了更细的模型网格，更深的测量深度增加了记录时长，因此，增加了单次迭代耗时，同时算法的总耗时随之增加。通过分析得出，对于更为精细的反演模型需要采用更细的网格，但同时造成了反演算法效率的下降，因此，需要根据实际的反演目标选择合适的网格大小。

3   实测数据全波形反演试验

3.1   地面塌陷实测数据

为了验证本文算法的实用性，采用实测的GPR剖面数据进行算法试验。本次试验路面数据采集自某城市路面塌陷项目，仪器选用MALA探地雷达仪器，测量区域表面为平整的混凝土路面，周围道路暂时封闭，无其他人为干扰。图9为现场实际工作记录图，其中红色方框为实际的未填充的塌陷区域，黑色有向线段表示选取测线的实际走向。此次探测的测量区域深度为2.00 m，测线长度为1.00 m，使用400 MHz的屏蔽天线，时窗长度为40 ns，采用距离测量模式获取589道数据，每道的采样点数设置为512。

图  9  实测路面现场工作记录

Figure  9.  Field work record for measured pavement

下载: 全尺寸图片 幻灯片

由图10可知，直达波同相轴具有明显的断裂，深部存在一个能量不强的异常反射。可以判断该段道路数据存在多个异常，且路面存在塌陷部分。剖面的末端浅部存在大量的异常波形，但由于同相轴没有发生断裂现象，可以推测异常不在地表，可能为地下浅埋管线。

图  10  实测路面的GPR剖面

Figure  10.  GPR profile of measured pavement

下载: 全尺寸图片 幻灯片

使用本文的褶积型目标函数FWI算法对该数据进行反演计算，采用均一模型作为反演的初始模型，激励源选用400 MHz雷克子波，得到的反演结果如图11所示。

图  11  实测数据全波形反演介电常数分布

Figure  11.  Dielectric constant distribution in full-waveform inversion for measured data

下载: 全尺寸图片 幻灯片

由图11可知，传统目标函数的反演结果中没有任何有效信息，这是由于实际工作中的仪器激励源子波未知，在反演的过程中所使用的子波估计所迭代的结果中的错误累计会逐渐增大，最终导致这一结果。而褶积型目标函数的反演结果反映了该路面不仅存在路面塌陷，在路面塌陷下方存在明显的充填不密实，并且在路面深部同样存在高介电常数的方形异常体，根据相对介电常数的分布可推测为规则排列的金属管线。同样，对于剖面末端的浅层异常，根据相对介电常数的分布，可推测为浅埋的金属电缆，并且通过查阅相关资料，最终确定了该处存在2处不同埋深的管线分布。

3.2   地下异常实测数据

选择MALA探地雷达仪器进行GPR实际探测，测量区域表面为干燥的沙土路面，采用测线模式采集数据。已知地下埋藏未知深度的管线，此次探测的测量区域深度为4.00 m，测线长度为8.00 m，使用600 MHz的屏蔽天线，时窗长度为60 ns，采用距离测量模式获取682道数据，每道的采样点数设置为512。采集到的实测数据如图12所示。

图  12  实测路面GPR剖面

Figure  12.  GPR profile of measured pavement

下载: 全尺寸图片 幻灯片

由图12可知，该数据表面存在层异常，并且能够看到3个反射异常。使用本文的褶积型目标函数FWI算法对该数据进行反演计算，采用均一模型作为反演的初始模型，采用600 MHz雷克子波作为反演的激励源，得到的反演结果如图13所示

图  13  实测数据全波形反演介电常数分布

Figure  13.  Dielectric constant distribution in full-waveform inversion for measured data

下载: 全尺寸图片 幻灯片

由图13可知，显示了3个管线异常的相对位置及埋深情况，且层位信息得到了较好的恢复，并且从反演得到的相对介电常数分布，可以推测出该管线为非金属管线，初步判断为PVC材质的小管径管线。

同时，为了说明本文提出的算法针对实测数据的反演效率，记录了模型参数和反演效率(表3)。

表  3  实测数据全波形反演效率

Table  3.  Full-waveform inversion efficiency for measured data

	网格大小	使用道数	记录时长/ns	迭代次数	单次迭代耗时/s	总耗时/s
实测1	200×200	200	40	40	35.590 2	1 502.007 3
实测2	200×400	400	60	40	247.106 4	11 723.320 0

下载: 导出CSV 

| 显示表格

通过对实测数据的FWI试验，证明了本文算法对于实测数据的适用性，同时，由于无法对实测数据激励源子波进行良好的估计，因此，采用传统目标函数FWI将对反演结果造成较大的影响，但是本文提出的褶积型目标函数则能够有效地避免这一问题，并且记录了相应的算法耗时，结果表明该算法效率虽然不能做到实时处理，但针对复杂且重要的未知地下异常分布的分析依然能够进行指导，为生产实际提供有效的依据。

4   结 论

a. 针对城市道路塌陷隐患探测问题，提出探地雷达适用的褶积型目标函数的全波形反演算法。

b. 该算法能够克服常规反演目标函数由于激励源子波估计不准确对全波形反演结果准确性造成的影响，取得良好的反演结果；在计算效率和反演准确度上，均优于常规的全波形反演。

c. 将该算法应用于GPR实测数据反演中，不仅较好识别异常体位置及层位信息，且可从反演的相对介电常数中推测出异常体材质，验证了本文方法的实用性。

d. 该算法可以得到较为准确的相对介电常数分布，与 GPR 剖面相互验证，提高反演精度，为实际工程应用提供理论依据和数据基础。