Envelope-waveform Inversion Based on Multi-offset Ground Penetrating Radar Data

1.   引言

探地雷达(Ground Penetrating Radar, GPR)是一种用于研究浅地表电磁特性的无损检测技术，主要基于电磁波传播原理，通过激发并记录高频电磁波(20～1000 MHz)在地下介质中(地下障碍物或不连续地质体)的反射、绕射等传播特性来探测电性结构分布。目前被广泛应用于土木工程场地勘查^[1]、地雷探测^[2]、环境与水文地质研究^[3,4]、农业评估^[5,6]、考古调查^[7,8]、极地考察^[9]以及航空航天^[10,11]等前沿工程与科学研究领域。

探地雷达的分辨率和探测深度极易受到脉冲源频率和地下电导率分布的限制，其中，高频电磁脉冲具有较高的分辨率，但其在地下介质中的衰减更快，因此穿透深度较小。根据测量方式的差异，探地雷达可进一步被划分为共偏移距和多偏移距两种类型。相较于共偏移距的测量方式，多偏移距测量可以灵活调整发射和接收天线之间的距离，捕获到的绕射波信息较为可靠，有助于实现精确的速度分析以及对地下电性参数的定量估计。因此开展基于多偏移距探地雷达数据的研究是十分必要的。

目前，探地雷达数据的解释工作已经从单纯推断目标体空间位置、形态、尺寸以及层界面等信息，上升到对介质电性参数(介电常数、电导率)的反演估计和对近地表结构的精细刻画。现行的反演方法主要可分为基于射线理论的层析成像方法和基于全波场理论的全波形反演方法。其中，全波形反演(Full-Waveform Inversion, FWI)最初起源于时域地震成像，后被推广到频率域，其同时兼顾了波场的运动学和动力学特征(振幅和相位)，是一种高分辨率、高精度的成像方法。全波形反演已在地震勘探领域获得广泛应用，但由于地面雷达测量方式对地下介质的照明能力有限，基于地面雷达数据的全波形反演具有挑战性。Lavoué等人^[12]采用拟牛顿算法针对多偏移距地面雷达数据开展频率域FWI，同时重构了介电常数和电导率的2维分布。Feng等人^[13,14]将井间地震FWI和GPR FWI集成到一个联合程序中，用于对不同地球物理参数(P波速度、介电常数和电导率)模型的定量成像；Feng等人^[15]将多尺度反演策略和吉洪诺夫正则化(Tikhonov regularization, TV)约束方法结合应用于时域GPR FWI过程，实现了介电常数和电导率的同步反演；王珣等人^[16]在时间域采用改进的全变差正则化策略，开展了GPR双参数多尺度同步反演研究；Huai等人^[17]在时间域提出了一种基于模型的层剥离FWI，有效增强了模型深部的重构质量和反演分辨率。

然而，当利用梯度类最优化算法求解基于Born近似框架下的全波形反演问题时，一个严重偏离正解的初始模型会在波形匹配过程中产生“周波跳跃”问题，使得目标函数倾向陷入到局部极值^[18]。因此，全波形反演高度依赖精确的初始模型来确保反演过程的正确收敛。低频成分对于恢复长波背景速度结构进而构建大尺度初始模型至关重要，但是实际采集的数据中所携带的低频信息十分有限，这就使得研究低频缺失情况下的全波形反演方法显得尤为重要。最初在地震勘探领域提出并发展起来的包络反演是用于重构低频成分比较知名的一类方法。Bozdağ等人^[19]最早构建了基于瞬时相位和包络的全波形反演目标函数，该种方法旨在利用Hilbert变换来实现波形相位和振幅信息的有效分离，从而降低了反问题的非线性程度。Chi等人^[20,21]和Wu等人^[22]构建了基于包络残差的全波形反演目标函数，并利用伴随状态法求取了梯度算子，利用包络反演成功重建了地下介质的大尺度背景信息。Wu等人^[23]提出了调制-褶积信号模型并在研究中赋予了包络数据以明确的物理意义。刘新彤等人^[24]研究了低频缺失下基于包络目标函数的全波形反演，在有效还原低频信息的同时提供了对地下介质的定量解释，然而这一研究仅针对跨孔雷达的观测模式。与跨孔雷达相比，地面雷达对地下介质的照明角度有限，这本身就增加了反问题的非线性和不确定性；当地面雷达缺失低频成分的情况下，也会极大地制约全波形反演的稳定性和精确度。因此，研究基于地面雷达数据的包络全波形反演具有重要意义。

此外，引用分频处理的多尺度反演策略也有利于缓解全波形反演的“周波跳跃”问题，提高反演的稳定性。频率多尺度方法可按低频到高频依次反演多个离散频率或频率组，其中，利用低频记录对局部极值的不敏感特性来重建平滑的大尺度背景结构，随后再利用高频记录刻画局部精细结构^[25]。Boonyasiriwat等人^[26]最初在地震勘探领域通过选择最优频带和维纳滤波来实现更高效的时间域多尺度全波形反演；随后，这一策略被成功借鉴到探地雷达的全波形反演中。Meles等人^[27]通过在模型更新过程中由低到高逐渐增加频率成分来扩展数据的频率带宽，从而实现了对跨孔-井地数据的时-频域联合反演。Lavoué等人^[12]开展了与频率采样有关的多参数全波形反演研究，通过数值实验证实了宽频带数据的同时反演可以有效增强双参数同步重建质量。Li等人^[28]提出了一种基于跨孔雷达数据的Laplace域波形反演，并提出了频率跳跃的多尺度波形反演策略，实现了反演过程的快速收敛。

本文在现有的探地雷达时间域全波形反演框架下，回顾了地震勘探中信号包络提取和常规包络反演方法的基本原理与实施流程，并将其成功推广到探地雷达领域，提出了基于地面多偏移距雷达数据的包络-波形反演方法，实现了对地下介电常数模型的有效重构。采用3层结构介电常数模型仿真模拟了数据低频信息缺失的情况，验证了包络-波形反演在增强数据所携带的大尺度构造信息方面的能力与优势，有效克服了常规全波形反演在缺失低频信息时无法正常工作的缺陷。此外，在所构建的包络-波形反演框架下，又有针对性地加入了频率多尺度策略，在重建大尺度宏观背景构造的基础上，也有利于实现对小尺度结构的精细刻画。本文还进行了抗噪性测试并验证了包络-波形反演在噪声条件下仍具有较强的反演能力。

2.   包络-波形反演

2.1   探地雷达信号包络提取

一个解析信号$\tilde y\left( t \right)$是由一个实信号$y\left( t \right)$和该实信号所对应的希尔伯特变换${y_{\rm{H}}}\left( t \right)$组成的

其中，${y_{\rm{H}}}\left( t \right)$可以表示为

其中，P代表柯西主值。

解析信号$\tilde y\left( t \right)$也可以通过瞬时振幅${{E}}\left( t \right)$以及瞬时相位$\phi \left( t \right)$来表示

其中，瞬时振幅${{E}}\left( t \right)$可表示为

其中，${{E}}\left( t \right)$也称为实信号$\tilde y\left( t \right)$的包络。值得注意的是，仅允许在时域而非频率域中将信号分离为上述相位和振幅的形式。

因此，探地雷达信号的包络也可以通过求取解析信号的振幅而获得，即

其中，${{{E}}_{{\rm{syn}}}}$表示探地雷达合成数据的包络，d(t)表示探地雷达合成数据， d_H(t)则表示合成数据的希尔伯特变换。由希尔伯特变换的定义可知：信号包络的提取实际上是一种非线性运算过程，这是因为d_H(t)相当于是将合成数据d(t)的相位旋转90°而获得的结果。

图1(a)和图1(b)分别给出了主频为80 MHz的源子波及其包络的波形与频谱数据；接下来，利用高通滤波器剔除了低于40 MHz的全部以及40～80 MHz的部分低频成分，缺失低频的子波及其包络的波形与频谱数据分别如图1(c)和图1(d)所示。从图1(d)不难看出：尽管源子波缺失了低频成分，但是其包络数据仍包含丰富的低频信息。

图  1  子波及其包络的波形与频谱信息

下载: 全尺寸图片 幻灯片

2.2   反演框架

常规包络反演目标函数的定义方式借鉴了全波形反演定义目标函数的基本思想，即采用观测数据包络与模拟数据包络的残差的2范数作为目标函数^[23]

其中，m表示模型参数(介电常数$\varepsilon$、电导率$\sigma$)，$\Delta E$表示包络的残差，E_obs和E_syn分别表示观测数据和模拟数据的包络，数据的包络可利用式(5)计算获得，p表示包络的能量阶数。

一般而言，可将反演过程视为解一个无约束优化问题，在反复迭代更新过程中，使得所构造的目标函数最小化。本文采用局部优化算法来实现模型的优化更新。在局部优化算法中，模型的更新方程可表示为

其中，m_k和m_k+1分别表示第k次和第k+1次迭代的模型参数，${\alpha _k}$表示更新步长，$\delta {m_k}$则表示模型的更新方向，为了避免应用复杂的局部优化算法对反演结果的影响，突出包络-波形反演策略的可靠性，本文采用最简单的最速下降算法来求取$\delta {m_k}$。

在包络反演算法框架中，求取梯度算子是至关重要的。类比于常规全波形反演方法，包络反演的梯度也可以通过1阶伴随状态法求得^[29,30]，即通过震源正传波场与伴随源反传波场的零延迟互相关来间接构造梯度。基于上述原理，横电磁波(Transverse Magnetic, TM)模式下包络反演介电常数的梯度公式可表示为

其中， $\left[ {\Delta Ed\left( t \right)E_{{\rm{syn}}}^{{{p - 2}}} - {{H}}\left\{ {\Delta E{d_{\rm{H}}}\left( t \right)} \right\}E_{{\rm{syn}}}^{{{p - 2}}}} \right]$可看作伴随源，而$\partial d\left( t \right)/\partial \varepsilon$则为敏感核函数，即线性Fréchet导数。

2.3   频率多尺度策略

在包络-波形反演框架下，又进一步引入频率多尺度策略，该策略的基本思想在于：借助滤波手段将观测数据和子波分解为所需的频率成分，在反演的初始阶段先利用低频成分构造大尺度背景速度构造，并将其作为后续反演的初始模型；接下来，再利用高频成分精细刻画地下介质中的小尺度弱扰动。这种多尺度策略从 “空间尺度”来细化地下介质模型，可以有效确保反演过程的正确收敛并提高反演结果的精度。

滤波是实现频率多尺度的关键，这里借鉴Boonyasiriwat等人^[26]的滤波思想，即通过维纳滤波器将子波和观测数据分解到不同的频带范围。在包络-波形反演框架下，借助子波包络设计的维纳滤波器公式为

其中，$\omega$表示角频率，f_Wi表示维纳滤波器，E_or表示原始子波去低频后取包络，这里假定预先已知，E_ta为目标频率子波的包络，$\varepsilon$表示控制数值溢出的稳定因子。

2.4   包络-波形反演实施步骤

本文所采用的包络-波形反演可概括为一种分段式反演策略，即先通过包络反演来重建介质的大尺度背景信息，并将其作为后续反演的初始模型，再利用常规全波形反演构造小尺度精细结构。值得注意的是：在低频缺失的情况下，由于常规全波形反演无法构建高频子波到低频子波的维纳滤波器，因此无法采用频率多尺度策略。低频成分缺失情况下包络-波形反演的具体实施步骤可概括如下：

阶段1　包络反演：

步骤1　输入：初始模型、观测记录、源子波；

步骤2　利用高通滤波器，构造缺失低频的观测记录；

步骤3　选取频率多尺度策略所需反演频率，设定迭代次数；

步骤4　计算模拟记录；

步骤5　利用式(5)计算缺失低频观测记录的包络和模拟记录的包络；

步骤6　确定当前迭代频率，根据式(9)构造维纳滤波器，作用于观测记录和模拟记录包络上，获得指定频率的观测记录和模拟记录的包络；

步骤7　根据梯度式(8)，计算整个模型的梯度；

步骤8　利用不精确线搜索方法计算更新步长；

步骤9　更新模型，迭代次数+1，并重复步骤4—步骤9，直到完成设定的全部迭代运算；

步骤10　输出反演结果，并将其作为第2阶段常规全波形反演的初始模型。

阶段2　常规全波形反演：

步骤1　输入：包络反演的结果作为初始模型、观测记录、源子波；

步骤2　利用高通滤波器，构造缺失低频的观测记录和源子波；

步骤3　设定迭代次数；

步骤4　计算模拟记录；

步骤5　利用梯度公式$\dfrac{{\partial F}}{{\partial \varepsilon }} = \displaystyle\sum\limits_{{\rm{sr}}} \mathop \int \nolimits_0^T E_{\rm{y}}^*\left( {{{s,r,t}}} \right)$$\dfrac{{\partial {E_{\rm{y}}}\left( {{{s,r,t}}} \right)}}{{\partial t}}{{{\rm{d}}t}}$，计算整个模型的梯度；

步骤6　利用不精确线搜索方法计算更新步长；

步骤7　更新模型，迭代次数+1，并重复步骤4—步骤7，直到完成设定的全部迭代运算；

步骤8　输出反演结果，并将其作为包络-波形反演的最终结果。

3.   数值测试

3.1   完整频率信息情况测试

为了验证包络-波形反演策略在地面多偏移距采集模式下的反演效果，采用如图2(a)所示的3层结构介质模型来进行数值测试。本文讨论的是时间域单参数包络-波形反演，因此仅反演单一参数介电常数，这里采用的是相对介电常数${\varepsilon _{{r}}} = \varepsilon /{\varepsilon _0},$${\varepsilon _0} \approx 8.85 \times {10^{ - 12}}{\rm{F/m}}$，而非绝对介电常数$\varepsilon$。图2(a)给出了相对介电常数真实模型，模型大小为10 m×11 m，水平和垂直方向的有限差分网格步长相同，均为h=0.1 m，共计生成101×111=11211个网格点。模型自上而下由3部分组成：模型顶部具有一个0～1.7 m厚的均匀介质层($\varepsilon _{{r}}^{{\rm{top}}} = 5$)；中间层的相对介电常数值为$\varepsilon _r^{{\rm{mid}}} = 5.5$，具有一个起伏底界面，深度介于4.1～5.9 m；第3层嵌入了两个圆柱形夹杂体，该层背景介质的相对介电常数为$\varepsilon _{{r}}^{{\rm{bot}}} = 5$，两个直径为0.5 m的目标体的中心分别位于(3.0 m, 7.7 m)和(6.0 m, 7.7 m)，相对介电常数值均为$\varepsilon _{{r}}^{{\rm{ta}}} = 7$。模型顶部设置了一个厚度为0.3 m的空气层。本文采用地面多偏移距观测方式，即在地表原点每隔2个网格点放置1个发射天线，共计56个发射天线，采用Blackman-Harris窗函数的1阶导数作为激励源，子波主频为80 MHz，记录时长500 ns，采样间隔为0.1 ns；每隔1个网格点放置1个接收天线，共计布设111个接收天线。采用如图2(b)所示的地下均匀半空间介质模型作为反演的初始模型($\varepsilon _{{r}}^{{\rm{init}}} = 5$)，在反演过程中电导率值始终被设置为常数($\sigma = 0.001\;{\rm{S}} \cdot {{\rm{m}}^{{{ - 1}}}}$)。

图  2  3层结构介质模型

下载: 全尺寸图片 幻灯片

为了确保反演过程的正确收敛及反演结果的置信度，本文采用了频率多尺度策略，即按照2次函数的变化趋势，从15～80 MHz中选择了10个反演频率，分别为15, 19, 24, 30, 37, 44, 52, 61, 70, 80 MHz。包络-波形反演的第1阶段包络反演和第2阶段常规全波形反演，分别利用上述的10个频率按由低到高的次序依次进行10次迭代，共计执行200次迭代运算，包络-波形反演结果如图3(a)所示。为了便于进一步对比不同方法的反演效果，这里给出了仅采用常规全波形反演的结果(图3(b))，反演过程中，上述10个反演频率各进行20次迭代，以确保与包络-波形反演的迭代运算次数一致。

图  3  3层结构介质模型的反演结果

下载: 全尺寸图片 幻灯片

两种反演策略都能够清晰地反映模型中异常的整体分布特征。其中，包络-波形反演主要利用低频成分重构大尺度背景信息，因此对于具有起伏底界面的中间层的刻画要明显优于常规全波形反演。图4给出了分别沿图2(a)所示的真实模型和图3所示的反演结果中不同位置(x=1.5 m, x=3.0 m, x=6.0 m)提取的参数纵向分布曲线；可以看出：包络-波形反演结果中间层左下方介质的分布要更接近真实模型。尽管如此，由于常规全波形反演方法本身包含了全部波形信息，因此其在小尺度结构的刻画精度和分辨率则更高。这里还分别计算了包络-波形反演结果和常规全波形反演结果与真实模型的相关系数，分别为0.9868和0.9828，两者与真实模型均极为相似，但包络-波形反演重构能力略高。

图  4  反演结果抽道对比

下载: 全尺寸图片 幻灯片

3.2   缺失低频成分情况测试

在全波形反演中，低频数据对保证反演的稳定性至关重要。研究表明：采用低频波场进行波形匹配时，发生“周波跳跃”的临界值较大，因而低频波场更有利于缓解反演的跳周问题。然而，在数据实际采集过程中，往往难以获得数据有效的低频分量，因此，研究低频信息缺失情况下的反演方法则显得十分重要。本节仍采用3层结构介质模型来测试低频信息缺失下包络-波形方法的反演能力。

首先，采用主频为80 MHz的子波对如图2(a)所示的介电常数真实模型进行正演计算，获得探地雷达合成数据。接下来，利用高通滤波器滤除数据的低频成分从而对低频缺失情况进行仿真模拟。图5给出了高通滤波前后第1个源对应的第1道雷达记录及其相应的频谱。从中可以看出原始雷达合成数据包含了完整的频率信息，频带范围为0～300 MHz(图5(b))；而高通滤波器有效剔除了低于40 MHz的全部以及40～80 MHz的部分低频成分(图5(d))；此外，通过对比滤波前后的单道雷达记录(图5(a)和图5(c))不难发现，缺失了低频成分的雷达记录的波形在振幅上也有所变化。

图  5  第1个源对应的第1道探地雷达合成记录及其频谱

下载: 全尺寸图片 幻灯片

基于2.4节给出的低频成分缺失情况下包络-波形反演的具体实施步骤，仍然采用3.1节所选取的10个反演频率：在包络反演阶段，10个反演频率先各进行10次迭代；由于低频缺失条件下常规全波形反演无法应用频率多尺度策略，因而常规全波形反演阶段仅利用缺失低频的源子波(原始主频为80 MHz)进行100次迭代运算。由此，包络-波形反演共计进行200次迭代更新，反演结果如图6(a)所示。这里仍然给出了仅采用常规全波形反演的结果(图6(b))，利用缺失低频的源子波(原始主频为80 MHz)进行了200次迭代运算。

图  6  低频缺失情况下模型的反演结果

下载: 全尺寸图片 幻灯片

图7给出了分别沿图2(a)所示的真实模型和图6所示的反演结果中不同位置(x=1.5 m, x=3.0 m, x=6.0 m)提取的参数纵向分布曲线。通过两种反演结果对比(图6)以及抽道对比结果(图7)不难发现：在低频信息缺失的情况下，常规全波形反演恢复的介电常数模型图像模糊，几乎无法识别关于地下异常和地层分布的有效反射信号(图6(b))；此外，与真实模型相比，重建结果的背景值被错误估计，推测反演结果可能陷入了局部极小值，由此可见，缺失低频波场对全波形反演的重建质量会造成一定的影响。相比之下，包络-波形反演很好地恢复了模型从浅到深的整体特征，对模型的还原度较高；反演图像中，具有起伏底界面的层状介质清晰可见，并且可以准确辨别模型深部嵌入的两个夹杂体的位置信息，与3.1节具有完整频带信息的反演结果(图3(a))相似度极高，这足以说明：包络-波形反演可以在低频成分缺失且常规全波形反演失效时发挥作用。此外，缺失低频情况下，两种反演结果与真实模型相关系数分别为0.9747和0.9351，相较于常规全波形反演，包络-波形反演有效降低了反演误差。

图  7  反演结果抽道对比

下载: 全尺寸图片 幻灯片

3.3   抗噪性测试

在实际探地雷达勘探中，采集的雷达记录中往往包含严重的干扰噪声，这对后续的反演和成像都将造成极大的影响。因此，本节将对基于探地雷达数据的包络-波形反演方法进行抗噪性测试。真实模型仍然采用如图2(a)所示的3层结构介质模型。发射和接收天线的布设方式、反演参数、反演策略的选取均与3.1节无噪测试相同。在合成的观测记录中加入了经过均值滤波处理的白噪声，这样做是为了降低白噪声的频率以更好地匹配原始观测记录的主频；噪声总能量约为原始合成记录总能量的0.5%。图8(a)和图8(b)分别给出了包络-波形反演以及常规全波形反演的重建结果。两种反演方法都能准确地恢复模型的宏观构造信息，对具有起伏底界面的层状介质的层位刻画都比较准确，但是常规全波形反演图像中，模型第2层左下角(4～6 m)介质的均匀性欠佳；图9给出了分别沿图2(a)所示的真实模型和图8所示的反演结果中不同位置(x=1.5 m, x=3.0 m, x=6.0 m)提取的参数纵向分布曲线；可以看出：模型深度范围在6 m以下的部分包络-波形反演(蓝色实线)要略优于常规全波形反演(绿色实线)。通过分别计算两种反演结果与真实模型的相似度也可有力证明上述结论(相关系数分别为0.9808和0.9784)。

图  8  添加噪声后模型的反演结果

下载: 全尺寸图片 幻灯片

图  9  反演结果抽道对比

下载: 全尺寸图片 幻灯片

4.   结束语

本文依据电磁波与地震波在动力学和运动学特征上的高度相似性，将地震勘探中用于解决低频缺失情况的包络反演方法迁延至探地雷达研究领域，构建了适用于地面多偏移距雷达采集模式的包络-波形反演框架。经过研究分析，本文得到如下几方面的结论：

(1)低频数据的缺乏会影响大尺度扰动介质的恢复，从而导致常规全波形反演对初始模型具有高度依赖性。雷达记录的包络波动和衰减携带了有效的低频信息，可用于估计大尺度(长波长)的背景结构。包络-波形反演既保证了地下大尺度强扰动构造的有效重建，又兼顾了全波形反演可以精细刻画小尺度弱扰动目标体的优势，因此在近地表电磁勘探领域具有很高的应用前景。

(2)本文通过对3层结构介电常数模型进行模拟研究发现：当雷达记录包含了完整的频带信息时，常规全波形反演方法在细节信息的刻画方面要优于包络-波形反演；但在低频信息缺失的情况下，常规全波形反演则无法提供对地下物性参数的准确估计。相反，包络-波形反演在低频成分缺失时仍能给出稳定且质量较高的反演结果。

(3)在包络-波形反演框架下，加入了频率多尺度的策略，实现了组间频率由低到高串行反演，该策略可以有效压制反演假象，提高反演精度。然而，低频缺失的情况下，常规全波形反演方法无法构建高频子波到低频子波的维纳滤波器，因而无法开展频率多尺度策略。通过抗噪性测试，进一步验证了包络-波形反演方法对白噪声具有不错的抵抗性。