低频缺失下跨孔雷达包络波形反演

刘新彤,刘四新, 孟 旭,傅 磊

1.吉林大学地球探测科学与技术学院，长春 130026 2.澳门科技大学太空科学研究所，澳门 999078

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以下的有效信息。因此，重新构建探地雷达数据中的有效低频信息，成为了利用全波形反演方法处理跨孔雷达数据的重点。

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

1 正演问题表达形式

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

(1)

式中：ε(x)和σ(x)分别表示介电常数和电导率的空间分布；s表示某一个特定的源；μ0表示磁导率，这里假定为常数；Es和Hs分别表示由电流密度源矢量Js产生的电场和磁场。式(1)可简写成如下形式：

(2)

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

MEs=Js或Es=M-1Js。

(3)

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

Es=g*Js。

(4)

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

M-1=g。

(5)

2 目标函数及梯度计算

2.1 传统全波形反演

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

(6)

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

(7)

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

2.2 包络波形反演

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

(8)

(9)

进一步推导获得

(10)

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

(11)

其中，

(12)

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

(13)

其中，

(14)

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

(15)

(16)

(17)

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

[ε(x)k+1]=[ε(x)k]-ζε,k·[Cε(x)k]，

(18)

[σ(x)k+1]=[σ(x)k]-ζσ,k·[Cσ(x)k]。

(19)

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

Cε(x)k=Sε(x)k+

(20)

Cσ(x)k=Sσ(x)k+

(21)

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

[Eenv(ε+κεCε,k,σ)-Eenv(ε,σ)]，

(22)

[Eenv(ε,σ+κσCσ,k)-Eenv(ε,σ)]。

(23)

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

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方法，尤其对于中间两个管线异常体，传统波形反演电导率结果优势明显。

图1 Ricker子波和包络波形(a)及其频谱(b)Fig.1 Ricker wavelet and its envelope (a) and spectrum (b)

图2 去除30 MHz以下低频成分Ricker子波及其包络(a)和频谱(b)Fig.2 Ricker wavelet without low frequency data under 30 MHz and its envelope (a) and spectrum (b)

a.相对介电常数；b.电导率。白色圆圈表示发射源，叉号表示接收器，下同。图3 原始模型Fig.3 Original model

a、b. FWI；c、d. EWI。图4 完整频率信息波形反演结果Fig.4 Waveform inversion results of complete frequency information

接着，我们提取了原始数据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的结果更好，且目标函数收敛，能够为地下地层与异常体从位置和物性参数方面提供定量解释。

a、b. FWI；c、d. EWI。图5 30 MHz以下低频信息反演结果Fig.5 Inversion result of low frequency data under 30 MHz

a、b. FWI；c、d. EWI。图6 去除30 MHz以下低频信息反演结果Fig.6 Inversion result without low frequency under 30 MHz

4 结论

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

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

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

[1]Tarantola A. Inversion of Seismic-Reflection Data in the Acoustic Approximation [J]. Geophysics, 1984, 49: 1259- 1266.

[2] Pratt R G, Shin C, Hicks G J. Gauss-Newton and Full Newton Methods in Frequency Domain Seismic Waveform Inversion [J]. Geophysical Journal International, 1998, 133: 341-362.

[3] Pratt R G. Seismic Waveform Inversion in the Fre-quency Domai: Part 1: Theory and Verification in a Physical Scale Model [J]. Geophysics, 1999, 64: 888-901.

[4] Pratt R G, Shipp R M. Seismic Waveform Inversion in the Frequency Domain: Part 2: Fault Delineation in Sediments Using Crosshole Data [J]. Geophysics, 1999, 64: 902-914.

[5] 冯晅, 鲁晓满, 刘财,等. 基于逐减随机震源采样法的频率域二维黏滞声波方程全波形反演[J]. 吉林大学学报(地球科学版), 2016, 46(6):1865-1873.

Feng Xuan, Lu Xiaoman, Liu Cai, et al. Frequency-Domain Full Waveform Inversion of 2D Viscous Acoustic Wave Equation Using Decreasing Random Shot Subsampling Method[J]. Journal of Jilin University (Earth Science Edition), 2016, 46(6): 1865- 1873.

[6] Kuroda S, Takeuchi M, Kim H J. Full Waveform In-version Algorithm for Interpreting Cross-Borehole GPR Data [J]. Seg Technical Program Expanded Abstracts, 1949, 24(1): 2668.

[7] Ernst J R, Maurer H, Green A G, et al. Full-Wave-form Inversion of Crosshole Radar Data Based on 2-D Finite-Difference Time-Domain Solutions of Maxwell’s Equations [J]. IEEE Transactions on Geoscience and Remote Sensing, 2007, 45(9): 2807-2828.

[8] Meles G A, Kruk J V D, Greenhalgh S A, et al. A New Vector Waveform Inversion Algorithm for Simultaneous Updating of Conductivity and Permittivity Parameters from Combination Crosshole/Borehole-to-Surface GPR Data [J]. IEEE Transactions on Geoscience and Remote Sensing, 2010, 48(9): 3391-3407.

[9] 吴俊军,刘四新,李彦鹏,等．跨孔雷达全波形反演成像方法的研究[J]. 地球物理学报, 2014, 57(5):1623-1635.

Wu Junjun, Liu Sixin, Li Yanpeng, et al. Study of Cross-Hole Radar Tomography Using Full-Waveform Inversion [J]. Chinese Journal of Geophysics, 2014, 57(5): 1623-1635.

[10] 孟旭, 刘四新, 傅磊,等. 基于对数目标函数的跨孔雷达频域波形反演[J]. 地球物理学报,2016,59(5):1875-1887.

Meng Xu, Liu Sixin, Fu Lei, et al. Frequency Domain Waveform Inversion of Cross-Hole GPR Data Based on a Logarithmic Objective Function [J]. Chinese Journal of Geophysics, 2016, 59(5):1875-1887.

[11] Ernst J R, Green A G, Maurer H, et al. Application of a New 2D Time-Domain Full-Waveform Inversion Scheme to Crosshole Radar Data [J]. Geophysics, 2007, 72(5): J53-J64.

[12] Virieux J, Operto S. An Overview of Full-Waveform Inversion in Exploration Geophysics [J]. Geophysics, 2009, 74 (6), WCC1-WCC26.

[14]Chi B X, Dong L G, Liu Y Z. Full Waveform Inversion Based on Envelope Objective Function [C]// 75th EAGE Conference & Exhibition Incorporating Spe Europec 2013. London: EAGE, 2013.

[15]Chi B X, Dong L G, Liu Y Z. Full Waveform Inversion Method Using Envelope Objective Function Without Low Frequency Data [J]. Journal of Applied Geophysics, 2014, 109: 36-46.

[16]Wu R S, Luo J R, Wu B Y. Seismic Envelope Inversion and Modulation Signal Model [J]. Geophysics, 2014, 79(3): WA13-WA21.