0 引言

跨孔探地雷达层析成像技术不断发展, 已经在水文、环境和工程物探中得到了验证与有效应用。在数据采集过程中, 跨孔雷达将两个天线放置在相隔一定距离的两个井孔中, 一个作为发射源, 另一个作为接收器。发射天线发出的雷达波穿过井孔之间的介质被接收天线采集, 该雷达波包含井孔之间介质的物性参数(如介电常数和电导率)空间分布信息。传统的走时层析成像和衰减层析成像技术虽然可以分别有效地反演出地下介质的介电常数和电导率分布情况, 但由于上述方法基于射线理论, 仅仅使用了雷达波的部分信息, 因而所得到的结果分辨率有限。跨孔雷达全波形反演技术在反演过程中使用雷达波的全部信息, 所得到的地下介质物性参数分布结果的分辨率可以达到亚波长级别。

近些年来, 国际地球物理探测学界对全波形反演(full waveform inversion, FWI)方法持续关注, 波形反演技术在很多领域都取得了显著进展, 尤其是在地震方面, 已经得到了相对成熟的理论架构。早期, Tarantola^[1]将广义最小二乘反演理论应用到了全波形反演方法中, 通过合成数据与观测数据计算出残差并进行反传获得反传波场, 利用反传波场与正传波场的互相关来构建梯度, 形成了一套完整的时间域波形反演方法理论架构。随后, Pratt等^[2-4]将全波形反演方法扩展到了频率域, 提出频率域全波形反演方法; 冯晅等^[5]提出基于逐减随机震源采样法的频率域波形反演方法。与地震全波形反演方法相比, 探地雷达全波形反演方法的研究起步较晚, 且目前仍然处于发展初期。Kuroda等^[6]采用全波形反演方法来解释跨孔雷达数据, 获得了介电常数成像结果。Ernst^[7]等使用全波形反演方法处理跨孔雷达数据, 采用级联更新的方式同时反演了介电常数和电导率, 并取得了较好的效果。Meles等^[8]利用矢量全波形反演技术处理探地雷达数据, 进行了同时迭代介电常数和电导率的全波形反演。在国内, 波形反演也取得了一定的进展。吴俊军等^[9]利用摄动法求解全波形反演的目标函数并详细推导了基础理论, 采用最速下降法在迭代过程中同时更新介电常数和电导率; 孟旭等^[10]进一步将跨孔雷达全波形反演方法发展到了频率域, 推导了基于对数目标函数的频率域全波形反演的基础理论, 使用共轭下降法同时反演介电常数和电导率。在处理实际数据方面, Ernst等^[11]提出一种二维全波形反演策略, 采用级联方式反演两种物性参数, 反演跨孔雷达实际数据。但是, 到目前为止, 由于实际数据信噪比较低, 国内外使用全波形反演方法处理探地雷达实际数据并没有得到良好的应用。

FWI方法广泛采用局部优化的方式进行反演, 但由于其本身对初始模型高度依赖, 面对先验信息缺失的地下介质时, 无法从准确的初始模型出发; 且在实际采集到的数据中包含噪声以及存在低频信息缺失的情况下, 局部优化方法很容易陷入局部极值^[12]。在目前的波形反演研究领域中, 普遍选择为其构建较为准确的初始模型, 来降低非线性程度以获取更好的结果。而对于构建大尺度的初始模型, 低频信息的使用尤为关键。对于实际数据, 由于带宽或数据采集滤波器设计等原因, 低频信息往往难以采集; 虽然在实验室条件下能够采集到一定量的低频数据, 但是难度极高, 成本巨大, 所采集到的数据并不完整且应用难度较大, 无法推广到探地雷达野外数据的处理中。以中心频率为100 MHz的钻孔雷达为例, 难以获得30 MHz以下的有效信息。因此, 重新构建探地雷达数据中的有效低频信息, 成为了利用全波形反演方法处理跨孔雷达数据的重点。

雷达数据的包络能够重构地下模型的大尺度信息, 用其反演长波长分量很合适。最早于2011年, Bozdağ等^[13]进行了包络目标函数的研究, 但是并没有利用包络目标函数进行全波形反演得到良好的结果, 也没有对其应用进行良好的解释与评价; 随后, Chi等^[14-15]以及Wu等^[16]分别使用基于声波方程全波形反演方法处理地震数据, 采用包络目标函数有效反演了地下介质的速度信息, 并使用该方法作为常规波形反演的初始模型, 得到了高分辨率的准确结果。但上述关于包络全波形反演方面的研究目前都集中在地震数据的波形反演。

本文在时间域跨孔雷达波形反演的基础上, 使用包络目标函数, 采用求导法, 详细推导了包络波形反演(EWI)方法的基础理论, 将EWI成功推广到了探地雷达(ground penetrating radar, GPR)领域, 并采用共轭梯度法同时反演了两种物性参数(介电常数和电导率); 通过模拟原始数据低频信息缺失的情况验证该方法在低频缺失情况下的反演能力, 反演结果在为地下介质做出评价的同时, 可以作为初始模型, 降低非线性程度, 进一步参与到传统波形反演当中, 来获取分辨率更高的结果。该全波反演方法的正演部分采用高阶时间域有限差分算法, 使用了CPML(convolutional perfect matched layer)吸收边界, 并在每一次迭代过程中, 对梯度进行了中值滤波等优化方式。由于反演过程中数据量巨大, 采用了Matlab的分布式计算引擎(matlab distributed computing engine, MDCE)并行运算的方式, 以提高在普通微机上的运行效率。

1 正演问题表达形式

在讨论探地雷达地球物理问题时, 地下介质的物性参数分布情况使用介电常数和电导率来描述, 并假设磁导率是均匀不变的。在时间域, 对于任意时间点t和空间x, Maxwell方程如下^[7]:

式中: ε(x)和σ(x)分别表示介电常数和电导率的空间分布; s表示某一个特定的源; μ₀表示磁导率, 这里假定为常数; E_s和H_s分别表示由电流密度源矢量J_s产生的电场和磁场。式(1)可简写成如下形式:

式中, M(ε, σ)代表Maxwell方程的一个线性算子。忽略磁场项:

在时间域, 空间任意位置的电场可以用源与对应位置格林函数之间的卷积来表示:

式中, g为M算子对应的格林函数。显然, M^-1包含离散近似的格林函数^[2-4]。因此,

2 目标函数及梯度计算 2.1 传统全波形反演

传统全波形反演的主要目的是寻找目标函数最小时对应的介电常数ε与电导率σ的空间分布。传统目标函数表达式为

式中: d表示接收器; τ表示观测时间; E(ε, σ)和E^obs分别为电场正演数据和实际数据。式(6)表示源s在接收点d处、观测时间为τ时全部误差的总和。目标函数对应的梯度由正传合成波场与反传伴随波场的互相关获得:

式中: p代表模型参数ε和σ; v为虚拟源向量, 对于不同的模型参数具有不同的表达式; r(ξ)=E(ξ)-E^obs(ξ), 代表反传波场的后向残场源。

2.2 包络波形反演

本文通过将包络目标函数引入到Meles等^[7]的矢量波形反演技术当中, 重新构建包络全波形反演基础理论。不同于Meles等使用扰动目标函数一阶近似来推导梯度的方法, 本文使用对包络目标函数求导的方式来推导梯度。式(6)为反演中最广泛使用的传统目标函数, 本文给出基于该目标函数基础上的包络目标函数, 即合成包络波场与实际包络波场差的2-范数:

式中: E_env(ε, σ)是E(ε, σ)的包络, 定义为, E_env(ε, σ)=, E_{s, d, τ}(ε, σ)是场源s、接收器d在时间点τ的源场数据, 是E_{s, d, τ}(ε, σ)的希尔伯特变换。通过目标函数(8)对模型参数p进行求导来获得梯度, 则

进一步推导获得

式中, 对合成波场的偏微分可以通过使用模型参数p对公式(3)两端求偏导获得:

其中,

在时间域, 空间任意位置的电场可以认为是源与对应位置格林函数之间的卷积:

其中,

将式(5)代入式(13)并写成积分的形式:

式中, 的时间方向与正传场相反, 表示后向传播场。因此, 式(14)可以认为是虚拟源v与后向传播场之间的零相位卷积。这与Meles等利用扰动法来推导目标函数所获取的梯度结果一致。根据式(14), 我们给出具体的梯度表达式:

式中, 代表残场之和R_s的后向传播过程,

为了实现最大程度的收敛, 本文使用共轭梯度法来更新模型:

式中: ζ_{ε, k}和ζ_{σ, k}分别为第k次迭代中介电常数和电导率的迭代步长; C_ε(x)_k和C_σ(x)_k分别为介电常数和电导率对应的共轭梯度方向, 由当前和上一次的梯度计算得到。

当k=1时, 满足C_ε(x)₁=▽S_ε(x)₁和C_σ(x)₁=▽S_σ(x)₁。通过沿各自C_ε、C_σ方向寻找极值点, 我们能够在同一次迭代中实现介电常数和电导率的同步反演。迭代步长由以下公式获得:

式中, κ_ε、κ_σ为不同的小稳定因子。在反演过程中必须小心地为其选择适当的值并且该值随着迭代不断更新。

3 数值模拟

与传统的FWI方法对比完整波形不同, 我们的方法比较的是合成数据的包络与观测数据的包络。这里我们以Ricker子波为例来分析包络转换的作用。图 1a显示了一道Ricker子波的波形与它对应的包络, 通过对比可以看出, Ricker子波信号本身变化剧烈, 而包络波形变化相对缓慢。图 1b为子波频谱与包络频谱的对比图, 可以很明显看出存在低频信息通过非线性转换被加入到了包络频谱中。

传统的跨孔雷达FWI在雷达数据低频信息充足且初始模型足够准确的情况下可以获得良好的反演结果; 但是当低频信息缺少时, FWI方法往往无法恢复大尺度的背景信息, 甚至错误地更新模型并导致反演最终无法收敛, 非常容易陷入局部最小。雷达数据缺少低频信息这种状况在很多雷达实例中存在, 面对类似情况的实例, FWI方法难以发挥作用; 但EWI方法有能力在数据缺少大量低频信息的情况下更新模型的长波长信息, 其原因是它的伴随波场含有低频信息。我们去除Ricker子波30 MHz以下的低频成分, 图 2a中为去除低频成分后的Ricker子波及其包络, 可以看出去除低频成分后, Ricker子波信号变动更剧烈, 而包络曲线依然平稳。图 2b为其频谱对比图, 与图 1b相比较, 可以非常明显地看出, 在除去低频成分后, Ricker子波频谱低频部分虽然被移除, 但其包络低频部分几乎不受影响, 仍然有明显的低频信息被添加。

为了验证EWI方法的成像能力, 我们建立了如图 3所示的复杂小模型, 模型尺寸为6 m×6 m, 包含3个地层:第一层与第三层物性参数相同, 相对介电常数均为5, 电导率为0.001 S·m^-1; 中间层的相对介电常数为5.5, 电导率为0.002 8 S·m^-1, 并埋藏两个柱状异常体, 位于深度3 m处, 间隔1 m, 相对介电常数为7, 电导率为0.008 S·m^-1。在本次模拟中共13个发射源, 每组发射对应13个接收器。发射器与接收器位置均为深度0~6 m, 并以0.5 m等间隔排列; 天线中心频率为100 MHz。

首先, 我们对比了两种方法基于均匀初始模型反演的结果。均匀初始模型的相对介电常数为5.5, 电导率为0.002 8 S·m^-1。初始模型物性参数与原始模型中间层背景物性参数一致。由于原始数据包含了完整且正确的长波长信息, 所以FWI和EWI两种反演方法均能够获得良好的反演结果。其中传统的全波形反演方法所得到的结果分辨率略高于基于包络的波形反演结果, 这是因为FWI方法使用的信息更加完整, 但也有更高的非线性程度。从图 4反演结果来看, FWI方法在波形信息完整时, 反演结果要优于EWI方法, 尤其对于中间两个管线异常体, 传统波形反演电导率结果优势明显。

接着, 我们提取了原始数据30 MHz以下的低频信息, 以验证FWI和EWI两种方法对低频信息的反演能力。如图 5所示, FWI方法在单独反演低频信息时, 反演能力明显下降, 图 5a、b是迭代20次的结果, 虽然能够反演出一定的地层和异常体信息, 但已经远远不及EWI方法(图 5c、d); 并没有选择更多的迭代次数是因为由于低频缺失, FWI方法容易陷入局部最小, 随着迭代的进行, 结果将会朝着错误的方向发展, 最终无法收敛, 完全无法准确反映地下信息。而包络波形反演能够获得良好的结果, 目标函数收敛, 且最终结果(图 5c、d)与图 4中使用全部波形信息反演的结果相差不大。该组模拟证明, EWI方法对低频信息的利用能力高于FWI方法。对比提取低频信息的反演结果与完整波形信息的反演结果可以看出, 无论使用FWI还是EWI方法, 由于高频信息的缺失, 反演对于中间细节的刻画明显下降, 极易陷入局部最小。

最后, 我们模拟现实中低频缺失的情况, 以了解当低频缺失、FWI方法无法获得理想反演结果的情况下, EWI方法是否仍然可以获得良好结果, 使跨孔雷达波形反演体系更加完善且更符合实际应用。图 6中, 由于相同的原因, FWI结果与图 5一致, 均为20次迭代的结果。由图 6可见:在低频信息缺失的情况下, 无论是介电常数还是电导率, FWI均无法获得与利用完整波形信息进行FWI时一致的结果, 目标函数不收敛; 与FWI相对比, EWI的结果更好, 且目标函数收敛, 能够为地下地层与异常体从位置和物性参数方面提供定量解释。

4 结论

1) 与全波形反演方法相比, 包络波形反演方法是一种对频率成分更加稳定且非线性更低的波形反演方法。跨孔雷达数值模拟证明该方法在跨孔雷达领域有着很高的应用价值, 尤其对于低频信息缺失的数据, 包络波形反演可以在传统波形反演不适用的情况下仍然提供良好的反演结果。

2) 虽然对于完整的波形信息, 全波形反演方法结果要优于包络波形反演的结果, 尤其在电导率的结果上全波形反演方法优势明显; 但在低频信息缺失时, 包络波形反演无论是介电常数还是电导率, 结果明显优于全波形反演方法, 可以为地下物性参数提供解释。

3) 包络波形反演方法是对跨孔雷达波形反演的一种补充方法, 使得跨孔雷达波形反演符合更多应用实例。