基于探地雷达的路基地下异常体全波形反演

2022-12-12胡磊磊黄德军

煤田地质与勘探 2022年11期

胡磊磊，陈康，黄德军，杨荣

(1.广西壮族自治区地球物理勘察院，广西柳州 545005；2.广西壮族自治区第七地质队，广西柳州 545005)

探地雷达(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 就是建立一个目标函数，通过最优化策略找到某个参数值，使得目标函数最小。对于合成数据，由于子波的形态是已知的，因此，可以假设实际采集到的探地雷达数据为Ez分量，对其建立如下的极小化目标泛函[21]：

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

式中：[ ψx,ψy,ϕ]为拉格朗日乘子函数；[e1,e2,e3]为伴随方程，其表达式为：

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

式中：δ为变量的增量。

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

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

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

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

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

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

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

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

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

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

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

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

2 合成数据实验

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

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

图1 路面塌陷异常相对介电常数分布Fig.1 Relative dielectric constant distribution in case of road collapse exceptions

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

图2 路面塌陷异常正演剖面Fig.2 Forward section of road collapse exceptions

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

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

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

图3 路面塌陷异常传统目标函数和褶积型目标函数FWI结果对比Fig.3 Comparison between results from conventional objective function and convolution type objective function

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

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

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

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

图4 路面塌陷异常FWI 重构误差曲线Fig.4 Error curve of ground surface collapse model reconstruction

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

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

表1 合成数据全波形反演效率对比Table 1 Comparison of full-waveform inversion efficiency for synthetic data

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

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

图5 地下空洞异常相对介电常数分布Fig.5 Relative dielectric constant distribution in case of underground void exceptions

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

图6 路面塌陷异常正演剖面Fig.6 Forward section of road collapse exceptions

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

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

图7 褶积型目标函数FWI 结果Fig.7 Convolution type objective function FWI result

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

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

图8 地下复杂介质模型重构误差曲线Fig.8 Error curve of reconstruction for underground complex medium models

表2 合成数据全波形反演效率Table 2 Full-waveform inversion efficiency for synthetic data

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

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

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

3.1 地面塌陷实测数据

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

图9 实测路面现场工作记录Fig.9 Field work record for measured pavement

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

图10 实测路面的GPR 剖面Fig.10 GPR profile of measured pavement

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

图11 实测数据全波形反演介电常数分布Fig.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 剖面Fig.12 GPR profile of measured pavement

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

图13 实测数据全波形反演介电常数分布Fig.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

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