确定地层速度模型是地震勘探的核心目标之一。全波形反演(full waveform inversion, FWI)是在充分利用叠前地震资料中的动力学和运动学信息的基础上反演地层速度模型的主要方法之一^[1], 是目前地球物理勘探领域的研究热点。在该领域, TARANTOLA^[2]创造性地提出了利用反传算法来计算目标函数的梯度, 避免了Jacobian矩阵的求取, 从而实现了时间域全波形反演。GAUTHIER等^[3]和MORA^[4]实现了实际二维地震资料的全波形反演, 从实践上证明了全波形反演具有精细刻画地下构造及岩性特征的能力, 是一种高精度的速度建模手段; 同时, 他们指出了全波形反演理论的固有缺陷, 如计算量巨大、解易陷入局部最小值、低频信息缺失等, 限制了全波形反演技术在实际资料中的应用。

近年来, 随着计算机科学与技术和优化算法的发展, 全波形反演理论和方法的研究有了新的发展, 尤其是GPU并行编码方式的应用使得全波形反演的计算效率大大提升。冯海新等^[5]提出了基于GPU/CPU和震源随机编码技术的混合域全波形反演, 该方法提高了反演精度, 减少了炮间串扰噪声, 在GPU的加速下, 计算效率显著提升。魏哲枫等^[6]使用随机边界技术反传、重建震源波场, 充分发挥了GPU的计算能力, 提高了全波形反演在三维情况下的计算效率和稳健性。即便如此, 基于局部寻优算法求解强非线性问题所伴随的解易陷入局部极小值的问题仍然是全波形反演理论进一步发展的瓶颈。而这种强非线性问题却对初始模型有着较高的要求, 若初始模型与真实模型偏差过大, 计算波场与观测波场相位差大于半个周期, 会产生“周波跳跃”现象, 导致反演陷入局部极小值而无法得到准确的结果。如何为全波形反演提供一个良好的初始速度模型以缓解反演过程中的“周波跳跃”问题是个难题。

前人发展了很多速度建模方法, 如偏移速度分析、叠加速度分析以及速度层析成像等方法均可以建立较为准确的初始速度模型, 并能在一定程度上避免全波形反演过程中的局部极小和“周波跳跃”问题^[7]。此外, VAN LEEUWEN等^[8]提出的波场重构反演(wavefield reconstruction inversion, WRI)是近年来较为新颖的反演理论。他们在WRI方法的目标函数中引入了波动方程, 使目标解的空间扩大到数据空间和模型空间, 然后通过重构真实波场来计算模型梯度, 故而不需要计算反传的残差波场, 提高了计算效率, 同时解更为稳定, 不易陷入局部极小值。段超然等^[9]利用波场重构反演在优化过程中拥有较大自由度的优势, 模拟出低频部分, 并以波场重构反演的结果作为初始模型进行全波形反演, 得到了较好的结果。李辉等^[10]提出了基于高斯束的速度层析方法, 将波动层析与射线层析的优点综合利用, 在Rytov近似和Born近似的基础上阐述了高斯束层析理论, 并应用于角度域成像道集偏移速度分析, 得到了理想的层析结果。

在地震数据中, 不同频率数据对应的周期不同, 相比于低频数据, 高频数据频率更高、周期更短, 在进行计算波场与观测波场匹配时更容易出现“周波跳跃”现象。BUNKS^[11]将多尺度反演方法应用于时间域全波形反演, 引入滤波技术对低频数据进行反演以得到大尺度背景结构, 再用所得结果作为新的初始模型, 对高频数据进行反演得到更为精细的模型结构。PRATT^[12]提出了一种频率域全波形反演方法, 该方法直接在频率域实现, 因而, 很容易实现从低频数据到高频数据的多尺度反演。多尺度反演方法可以有效降低FWI对初始模型的依赖, 可以缓解计算波场与观测波场间的“周波跳跃”现象, 最终获得较好的反演结果。但是, 当数据的低频成分不足时, 利用低通滤波技术难以有效反演出速度模型的长波长结构, 从而导致后续反演结果不佳。研究人员发现利用数学变换提取低频信息, 再与全波形反演结合的方法能取得更好的结果。SHIN等^[13]提出的拉普拉斯域全波形反演方法利用阻尼波场中的零频率分量成功恢复出速度模型的长波长结构, 但该方法的反演深度严重依赖于偏移距和拉普拉斯阻尼系数的选择。CHI^[14]认为, 虽然地震数据自身的低频信息不足, 但在地震数据的包络中却含有丰富的低频成分。在此基础上, 构造地震数据的包络目标函数并进行反演, 用以恢复地下速度模型的长波长结构, 再以此作为常规全波形反演的初始模型, 能有效降低反演的非线性程度, 并显著地提高反演结果的精度。胡勇等^[15]提出了利用解调包络来重构地震记录中缺失的低频数据, 该方法是利用三次样条插值来提取地震包络, 解决了希尔伯特包络在非对称信号情况下出现的振荡现象。但是, 当数据主频过高时, 单纯地采用包络反演与常规全波形反演相结合的方式, 仍会受“周波跳跃”的干扰, 得不到高精度的反演结果。

受前人工作的启发, 本文先对高频数据作低通滤波处理, 采用低频数据的对数包络进行反演, 重构出速度模型的长波长结构, 为后续全波形反演提供一个较为准确的初始模型; 再采用巴特沃斯低通滤波器实现从低频到高频的分频多尺度全波形反演。此反演方法降低了反演过程陷入局部极小值的可能, 抑制了“周波跳跃”现象对反演结果的不良影响。数值试验结果证实, 该方法可以更好地反演出地下复杂结构的速度模型。

1 包络反演理论

包络反演目标函数的定义和梯度的求解与常规全波形反演所采取的策略极为相似, 两者都是基于最小二乘准则建立目标函数, 并用反传算法计算梯度。在本节中首先简单回顾了常规全波形反演理论的目标函数及其梯度求解方法, 然后介绍如何利用希尔伯特变换提取地震数据的包络, 最后推导出包络反演中目标函数的建立及其梯度的求解。

1.1 常规全波形反演的目标函数建立及其梯度求解

在常规全波形反演中, 目标函数(J)通常根据计算波场与观测波场之间残差能量最小来定义:

$ J=\sum\limits_{i=1}^{{{n}_{\text{s}}}}{\sum\limits_{j=1}^{{{n}_{\text{r}}}}{\|{{\mathit{\boldsymbol{u}}}_{i, j}}-{{\mathit{\boldsymbol{d}}}_{i, j}}{{\|}^{^{2}}}}} $

(1)

式中:n_s, n_r分别表示激发点数和接收点数; u_{i, j}与d_{i, j}分别表示第i炮、第j道的计算波场和观测波场^[16]。

对目标函数关于模型参数m求导得到其梯度:

$ \frac{\partial J}{\partial m}=\sum\limits_{i=1}^{{{n}_{\text{s}}}}{\sum\limits_{j=1}^{{{n}_{\text{r}}}}{\left[ \frac{\partial {{\mathit{\boldsymbol{u}}}_{i, j}}}{\partial m}\cdot {{\mathit{\boldsymbol{r}}}_{i, j}} \right]}} $

(2)

式中: r_{i, j}= u_{i, j}- d_{i, j}, 即残差波场^[16]。在实际计算中, 一般使用反传算法来求解目标函数的梯度^[2]。其简略步骤为:正向传播源波场, 反向传播残差波场r_{i, j}, 最后计算二者的零延时互相关得到梯度。

1.2 地震数据的包络

对于任意实信号f(t)都有一个对应的解析信号$ \hat{f}\left( t \right) $, 该解析信号是一个不含负频率分量的复信号:

$ \hat{f}\left( t \right)=f\left( t \right)+\text{i}H\left\{ f\left( t \right) \right\} $

(3)

H{f(t)}表示实信号f(t)的希尔伯特变换, 希尔伯特变换定义为:

$ H\left\{ f\left( t \right) \right\}=-\frac{1}{\text{ }\!\!\pi\!\!\text{ }}K\int _{_{-\infty }}^{^{+\infty }}\frac{f\left( \tau \right)}{t-\tau }\text{d}\tau $

(4)

式中:K为柯西主值。解析信号$ \hat{f}\left( t \right)$可表示为由瞬时振幅和瞬时相位构成的关系式:

$ \hat{f}\left( t \right)=E\left( t \right){{\text{e}}^{^{\text{i}\mathit{\Phi }(t)}}} $

(5a)

$ E\text{ }\left( t \right)=\sqrt{{{\left[ f\left( t \right) \right]}^{2}}+{{\left[ H\{f\left( t \right)\} \right]}^{2}}} $

(5b)

$ \Phi \left( t \right)=\text{arctan}\frac{H\left\{ f\left( t \right) \right\}}{f\left( t \right)} $

(5c)

式中:E(t)为瞬时振幅, 即为信号f(t)的包络; Φ(t)为瞬时相位^[17]。图 1显示了单道地震数据及其包络, 可以看出, 数据的包络反映了波形在时间域上的宏观变化。图 2显示了单道地震数据及其包络的频谱, 可以看出, 数据包络内含有丰富的低频信息。

不同类型的包络反映了不同的波形宏观变化, 所提取的低频成分也不同。设E(t)^p表示数据包络的p次方, 理论上讲, p可以取任意正整数, 但在实际反演中, p值一般取1或2, 分别表示数据包络的一次方和二次方。包络的指数形式(图 3)表明:较大的p值(p=2)会给予早期到达的数据(如直达波)更大的权重从而压制后期到达的数据, 导致深层的反演效果不佳; 较小的p值(p=1)则会压制低频而突出后期的反射波, 难以重构出速度模型的长波长结构。设log_aE(t)表示数据包络的对数, 取自然数e为底数(a=e), 得自然对数包络, 简称为对数包络^[18]:

$ \text{lo}{{\text{g}}_{a\text{=e}}}E\left( t \right)=\text{ln}E\left( t \right) $

(6)

从图 4可以看出, 地震数据的对数包络同样能够反映出波形的宏观变化。图 5与图 6分别给出了不同类型的包络及其频谱。对比对数包络和指数包络可以看出, 对数包络在反映波形宏观变化的同时, 没有像包络二次方那样给予早期数据过大的权重, 压制后期数据; 与包络一次方相比, 对数包络更为温和而不尖锐, 有效提取了数据的低频成分。因此, 本文研究采用对数包络。

此外, 取正演模拟的单炮记录(采样间隔0.5 ms, 采样长度1.5 s), 从更大尺度上来分析地震波场及其包络所包含的信息。由Marmousi模型(部分)正演得到的单炮地震记录及其包络(图 7)可以看出, 单炮记录完整地记录了地下结构的复杂变化, 并且分辨率较高; 但包络只含有极少的细节信息, 其主要信息反映的是地下大尺度的背景构造。

1.3 包络反演目标函数的建立及其梯度的求解

与常规全波形反演类似, 包络反演的目标函数J_env是根据计算波场包络与观测波场包络之间残差能量最小来定义:

$ {{J}_{\text{env}}}=\sum\limits_{i=1}^{{{n}_{\text{s}}}}{\sum\limits_{j=1}^{{{n}_{\text{r}}}}{\|\text{ln}E({{\mathit{\boldsymbol{u}}}_{i, j}})-\text{ln}E({{\mathit{\boldsymbol{d}}}_{i, j}}){{\|}^{^{2}}}}} $

(7)

对方程(7)关于模型参数m求导, 可以得到目标函数的梯度:

$ \frac{\partial {{J}_{\text{env}}}}{\partial m}=\sum\limits_{i=1}^{{{n}_{\text{s}}}}{\sum\limits_{j=1}^{{{n}_{\text{r}}}}{\left[ \frac{\partial \text{ln}E({{\mathit{\boldsymbol{u}}}_{i, j}})}{\partial m}\cdot {{\mathit{\boldsymbol{r}}}_{\text{env}}}_{_{i, j}} \right]}} $

(8)

式中: r_{envi, j}=lnE(u_{i, j})-lnE(d_{i, j})表示计算波场与观测波场两者数据包络间的残差。与常规全波形反演类似, 包络反演目标函数的梯度也是通过反传算法进行求解。

1.4 基于地震数据包络的多尺度全波形反演方法

虽然地震数据自身的低频信息不足, 但包络中却含有丰富的低频成分, 地震数据的包络所提取出的低频信息与地震数据的主频可以相互对应。图 8为不同频率下地震数据及其包络的频谱, 可以看出, 当地震数据的主频较高时, 其包络的频率也较高。因此, 较高频率的数据如果不作任何处理, 仍会受“周波跳跃”现象的干扰而影响反演的稳定性。

本文借鉴BUNKS^[11]的分频多尺度反演思想, 对高主频的地震数据首先做低通滤波处理, 再进行包络反演(低频包络), 然后, 以包络反演的结果作为全波形反演的初始模型, 并在后续反演过程中, 引入巴特沃斯低通滤波器, 进行分频多尺度的全波形反演。

2 数值计算

选取Marmousi模型的一部分(图 9a)作为真实模型(图 9b)进行反演数值计算, 采用不含任何构造信息的一维线性模型作为初始模型(图 9c), 模型的主要参数为:网格点数320×120, 空间网格大小5 m。在模型周围还设置了10个网格点宽度的卷积完全匹配层(convolutional perfectly matched layer, CPML)作为吸收层^[19]。采用地表激发地表接收的观测系统, 在地表均匀分布了19个炮点和140个接收点, 时间采样间隔为0.5 ms, 采样时间为1.5 s, 震源采用主频为30 Hz的Ricker子波。

2.1 包络反演构造初始模型

在进行包络反演构造初始模型的数值试验里, 采用两种方式进行对比, 一种是对数据不作低通滤波处理, 直接进行包络反演; 另一种是先对数据作低通滤波处理, 滤去其高频成分再进行包络反演(低频包络)。图 10为普通包络、低频包络构造的初始速度模型, 可以看出:由普通包络反演得到的速度模型能够在一定程度上重构速度模型的大尺度构造, 但是效果较差, 其根本原因在于在数据主频较高的情况下, 普通包络提取的低频成分不足, 反演过程受“周波跳跃”现象干扰; 低频包络可以提取数据的极低频成分, 反演结果较好地反映了速度模型的长波长结构, 所展示的背景构造更为平滑。

2.2 反演结果对比

在6个尺度(0~5, 0~10, 0~15, 0~20, 0~25, 0~30 Hz)下进行分频多尺度反演, 将由低频数据反演得到的结果作为高频数据反演的初始模型, 逐级反演, 提高反演过程稳定性, 最终得到准确的反演结果。按照设计的数值实验参数, 分别进行了常规全波形反演、普通包络+多尺度反演和低频包络+多尺度反演, 结果如图 11所示。

常规全波形反演结果(图 11a)只能刻画出真实模型的大致形态和轮廓, 由于缺乏准确的初始模型且数据主频较高, 因此反演结果受“周波跳跃”现象的影响较大, 反演结果精度不高, 未能刻画出精细的构造。普通包络+多尺度反演结果(图 11b)显示, 当初始模型含有一定的大尺度构造(图 10a), 并在后续反演过程中引入分频多尺度策略时, 高频数据的反演过程变得稳定, 反演精度得到大幅度提升; 在模型中部区域(400~1 200 m), 复杂结构基本得以重构, 地层界面清晰。低频包络+多尺度反演结果(图 11c)显示, 一个良好的初始模型(图 10b)对反演精度的提升效果极其显著, 此反演方法所重构出的速度模型, 整体构造形态与真实模型基本吻合, 模型的横/纵向分辨率高, 各地质层分界明显, 浅/深部区域的细节构造得以准确刻画。

2.3 反演结果量化分析

在真实模型(图 9b)中抽取地面距离为x=200 m, x=900 m和x=1 400 m处的速度对图 11所示的3个反演结果与真实模型进行速度曲线拟合分析, 结果如图 12所示。

从图 12可以看出:低频包络反演结果的速度曲线拟合程度最高, 其次为普通包络反演, 常规全波形反演结果的速度拟合程度最低; 在x=200 m处, 由于靠近模型边界, 覆盖次数较低, 低频成分难以获得, 具有高非线性性, 常规全波形反演与普通包络反演得到的速度曲线只能得到速度随深度的大体变化且与真实模型相差甚远, 而低频包络反演所得到速度曲线可以准确反映出速度的宏观趋势和微小变化, 并且高度拟合真实模型的速度变化曲线。

表 1给出了x=200 m, x=900 m和x=1 400 m处的反演结果与真实模型的速度扰动比。从表 1可以看出, 相比另两种反演方法, 低频包络+多尺度反演结果的速度扰动比大幅度降低, 凸显了该方法的优越性。

综合图 11、图 12和表 1可以明显看出3种反演方法的优劣差异:相比常规全波形反演, 利用包络提取地震数据的低频成分作初步反演, 恢复出大尺度结构, 为后续反演提供初始模型, 最终反演结果的精度得以提升且对细节的刻画能力更强; 对于高主频地震数据, 低频包络的反演结果比普通包络的反演结果更为平滑, 更能反映出速度模型的长波长分量, 作为后续反演的初始模型有着更好的结果; 包络反演与多尺度反演方法的有机结合, 极大地避免了“周波跳跃”现象对反演过程的影响, 最终提高了反演精度; 与前两种方法相比, 低频包络对覆盖次数较低区域的反演精度提升十分明显, 减少了对大偏移距数据的依赖。

3 结论

1) 利用数据的对数包络提取数据的低频信息可以重构出模型的大尺度背景, 而低频包络反演的结果则能反映出速度模型的大尺度轮廓和背景构造, 可作为全波形反演的初始模型。

2) 频率越高的数据在全波形反演中越容易受“周波跳跃”的影响。对于高频地震数据, 分频多尺度反演以低频数据反演结果作为高频数据的初始模型, 逐级反演, 可以有效地解决数据主频过高、反演不收敛、最终结果精度过低的问题。

3) 本文采用包络反演与多尺度反演有机结合的反演策略, 主要针对低频成分占比不高且数据主频较高的地震数据难以用常规全波形反演成功反演的情况, 该反演方法可操作性高, 能够明显提升反演的精度。