随着油气勘探技术的发展，愈发复杂的地下构造对资料处理方法提出了更高的要求。在地震资料处理时，要得到精确的地下成像结果，高精度的速度建模是必要的。目前存在大量的速度建模方法，如叠加速度分析，旅行时层析反演和全波形反演等。全波形反演(Full waveform inversion, FWI)以求解波动方程为基础，通过最小化观测数据与模拟数据的L2范数来构建地下介质速度模型^[1]。相较于层析反演等速度建模方法，FWI能够重建更高精度的介质参数模型，满足复杂构造成像对高精度速度模型的要求，且具备岩性反演与储层预测的潜力。目前声波FWI被广泛应用于重建地下纵波速度模型，但声学介质假设无法描述横波，导致多分量地震资料中的横波在处理时成为噪音，干扰最终结果。因此，声波全波形反演已经无法满足多分量地震资料的反演要求，开展弹性波全波形反演(Elastic full waveform inversion, EFWI)研究具有重要意义。

基于Born近似框架的全波形反演算法是一个局部优化问题，隐含着初始模型接近真实模型的假设，当初始模型与真实模型差距较大时，全波形反演容易陷入局部极值解，无法收敛至正确值。为了使反演能够向正确方向收敛，初始模型的构建是关键。缺失低频信息的多分量地震数据使得EFWI极易陷入局部极值，导致反演失败，因此对缺失低频的多分量地震数据进行低频补偿十分必要。同时，全波形反演通过模拟记录与观测记录的匹配来更新速度模型，而模拟记录在一定的子波假设条件下通过求解波动方程得到，因此子波的设置同样会影响弹性波全波形反演的精度，错误的子波会增加“周波跳跃”风险，导致反演失败，而从野外单炮记录中直接获取准确子波非常困难。因此，针对低频缺失和子波未知情况下的地震数据处理算法研究十分必要。

全波形反演的目标函数在低频数据反演时具有更少的局部极值^[2]，但受到采集技术的制约，低频地震数据往往难以获得，或由于低频段地震数据具有较低信噪比而在预处理阶段被滤除，这增加了弹性波全波反演失败的几率。因此许多学者针对低频缺失问题展开了低频补偿研究，Shin等^[3]提出Laplace域反演方法，通过对地震信号施加一个随时间变化的指数衰减项来增加输入数据中的低频成分，并利用这些低频成分恢复速度模型中的低波数信息。但由于该方法存在衰减参数选择困难等问题，被进一步发展为Laplace-Fourier域反演^[4]。Bozdag等^[5]提出利用观测记录的包络来反演背景速度场的思想，Wu等^[6]发展了基于包络的全波形反演算法，推导了包络数据残差的全波形反演目标函数，并通过伴随状态法实现了梯度的有效求取，得到了包含大尺度信息的模型。包乾宗等^[7]研究了对数包络目标函数，并将其与多尺度方法相结合，有效缓解了“周波跳跃”问题。刘张聚等^[8]建立了一种基于时域加权的拉普拉斯——频率域弹性波全波形反演方法，在混合域全波形反演方法的基础上引入拉普拉斯衰减因子降低了全波形反演对于低频成分的依赖，并应用地震道差值模的倒数作为加权系数进行时域加权消除阻尼因子的负面影响。梁展源等^[9]基于地震数据的瞬时振幅包络不随频率移动而改变的性质，提出了频移目标函数，使包络反演能够适应较粗的网格分布，在保持精度的同时降低了包络反演的计算量。但低频补偿方法往往隐含了准确子波假设，当子波不准时，上述低频补偿方法需要额外引入不依赖子波算法，增强了反演算法的不稳定性。此外，基于压缩感知和稀疏性假设的地震数据拓频方法被提出^[10-12]，Lu等^[10]利用稀疏脉冲反褶积提高叠后地震剖面的分辨率；韩立国等^[11]利用反射系数的稀疏性，在Fourier域实现地震数据的全带宽拓频，将拓频得到的低频信息与原数据高频成分构成最终的低频补偿地震数据；Zhang等^[12]分析了低频信息在地震资料处理中的重要性，并将基于压缩感知的低频拓展方法用于处理共成像点道集，进而提高叠加剖面的质量。上述研究多用于拓频和处理叠后剖面，但基于压缩感知和稀疏性假设的低频补偿方法在全波形反演中的应用仍有待进一步研究。

本文提出一种基于盲稀疏脉冲反褶积的地震数据低频补偿方法。在地震信号稀疏性假设的前提下，采用L₁范数约束^[13]构造目标函数，将地震信号压缩为时域尖脉冲；然后利用富含低频成分的子波与时域尖脉冲褶积的方法实现宽频数据的低通滤波，得到低频重构地震记录，再利用重构记录进行EFWI得到初始纵横波速度模型。本文算法用于低频重构的子波是已知的，因此将重构数据作为全波形反演的输入数据时无需额外引入不依赖子波算法^[14]，保证了EFWI的稳定进行。

1 常规弹性波全波形反演基本原理

二维各向同性介质下的位移-应力弹性波方程为

$ \left\{\begin{array}{l} \rho \frac{\partial^2 u_x}{\partial t^2}=\frac{\partial \tau_{x x}}{\partial x}+\frac{\partial \tau_{x z}}{\partial z} \\ \rho \frac{\partial^2 u_z}{\partial t^2}=\frac{\partial \tau_{z z}}{\partial z}+\frac{\partial \tau_{x z}}{\partial x} \\ \frac{1}{\rho} \tau_{x x}=v_{\mathrm{p}}^2 \frac{\partial u_x}{\partial x}+\left(v_{\mathrm{p}}^2-2 v_{\mathrm{s}}^2\right) \frac{\partial u_z}{\partial z} \\ \frac{1}{\rho} \tau_{z z}=v_{\mathrm{p}}^2 \frac{\partial u_z}{\partial z}+\left(v_{\mathrm{p}}^2-2 v_{\mathrm{s}}^2\right) \frac{\partial u_x}{\partial x} \\ \frac{1}{\rho} \tau_{x z}=v_{\mathrm{s}}^2\left(\frac{\partial u_x}{\partial z}+\frac{\partial u_z}{\partial x}\right) \end{array}\right.。$

(1)

式中：ρ为介质密度，u_x、u_z分别表示质点在x和z方向的位移分量。τ_xx、τ_zz分别表示x和z方向的正应力，τ_xz表示切应力，v_p和v_s分别表示纵波速度和横波速度，t表示时间，x和z表示笛卡尔坐标系下的x和z方向。

利用模拟记录与观测记录的L₂范数构建时间域弹性波全波形反演目标函数为

$ J(\boldsymbol{m})=\frac{1}{2} \sum\limits_{n s} \sum\limits_{n r} \int_0^{\mathrm{T}}(\boldsymbol{u}-\boldsymbol{d})^2 \mathrm{~d} t。$

(2)

式中：m表示模型参数；ns和nr别表示炮数和道数；u和d分别表示模拟数据和观测数据。

利用伴随状态法求得目标泛函的梯度为

$ \boldsymbol{g}=\frac{\partial J}{\partial \boldsymbol{m}}=\int_0^{\mathrm{T}} \frac{\partial^2 \boldsymbol{u}}{\partial_t{ }^2} \hat{\boldsymbol{u}} \mathrm{~d} t 。$

(3)

式中u为正传波场；$\hat{\boldsymbol{u}}$为伴随波场，由伴随震源反传产生。

公式(2)所示目标函数对应的伴随源为R=u－d。最终的模型迭代更新公式为

$ \boldsymbol{m}_{j+1}=\boldsymbol{m}_j-\alpha_j \boldsymbol{g}_j 。$

(4)

式中，j表示迭代次数。相较于声波全波形反演，在多个参数相互耦合的情况下实现更多参数的反演使得EFWI的非线性问题更加突出。前人通过分步反演策略来降低EFWI的非线性程度^[15]，在这类策略中，低频数据对于重建速度的低波数背景成分至关重要，本文采用盲稀疏脉冲反褶积算法，对缺失低频信息的多分量地震数据进行低频补偿并应用于EFWI中，有效恢复了速度模型的低波数成分。

2 盲稀疏脉冲反褶积方法原理

本文提出基于盲稀疏脉冲反褶积的全波形反演方法，根据褶积模型^[16]，地震数据可以表征为有限频带的地震子波与宽频反射系数序列褶积的结果：

$ \boldsymbol{d}(t)=\boldsymbol{w}(t) * \boldsymbol{r}(t)。$

(5)

式中：d(t)表示地表接收到的地震信号；w(t)是子波；r(t)是地下反射系数序列；*表示褶积运算。

反褶积的目标就是通过反问题求解方法得到宽频反射系数序列，实现拓频。在单炮地震记录中，r(t)具有更广泛的含义，它不仅包含非自激自收情况下某一特定入射角的反射系数信息，而且包含了该入射角情况下的入射波能量等信息，如果将地震信号d(t)视为输入的震源信号w(t)经过大地滤波器作用后的输出信号，则r(t)相当于大地滤波器的固有响应，我们称之为传播与反射响应(Propagation and reflective response, PRR)，为便于区分，我们用i(t)表示，则大地滤波器的输入输出关系为

$ \boldsymbol{d}(t)=\boldsymbol{w}(t) * \boldsymbol{i}(t)。$

(6)

式(6)表明，地震记录可以视为震源子波w(t)对传播与反射响应i(t)进行带通滤波的结果。如果能够从地震记录中反褶积得到i(t)，然后将i(t)与低频震源子波褶积，就可以得到低频重构记录。因此，关键在于准确提取i(t)。

将式(6)改写为矩阵形式

$ \boldsymbol{d}=\boldsymbol{W} \boldsymbol{i} 。$

(7)

式中：$\boldsymbol{d}=\left[\begin{array}{ccc} d_{1, 1} & \cdots & d_{n r, 1} \\ \vdots & \ddots & \vdots \\ d_{1, T} & \cdots & d_{n r, T} \end{array}\right], \boldsymbol{i}=\left[\begin{array}{ccc} i_{1, 1} & \cdots & i_{n r, 1} \\ \vdots & \ddots & \vdots \\ i_{1, T} & \cdots & i_{n r, T} \end{array}\right]$；nr表示地震道数；W表示托普利兹矩阵，本文采用近偏移距直达波来构建该矩阵，表达式为

$ \boldsymbol{W}=\left[\begin{array}{cccc} w_1 & 0 & \cdots & 0 \\ \vdots & w_1 & \ddots & \vdots \\ w_m & \vdots & \ddots & 0 \\ 0 & w_m & \vdots & w_1 \\ \vdots & \ddots & \ddots & \vdots \\ 0 & \cdots & \ddots & w_m \end{array}\right]。$

(8)

则从d中提取i(t)等价于求解如下目标函数

$ \min\limits _i\|\boldsymbol{d}-\boldsymbol{W} \boldsymbol{i}\| 。$

(9)

我们利用压缩感知^[17-18]理论求解式(9)，压缩感知理论利用信号的稀疏性，通过使稀疏域中非零值个数最少来重建信号。由于地震反射系数往往具有稀疏性^[11]，所以式(9)的求解满足压缩感知的假设条件。对于式(9)的求解，L₀范数虽然具有最好的稀疏性，但基于L₀范数稀疏约束项的地震反演数学模型是一个非指定性多项式难题，该数学模型无法直接判断是否有解，只能通过假设来猜测它的解并通过模型验证正确性。L₁范数是L₀范数的凸逼近，既保证了反射系数的稀疏性，又解决了基于稀疏约束的地震反演数学模型无法求解的问题。此外，受接收仪器限制和采集环境等的影响，观测数据中往往存在噪音干扰，导致提取的地下传播与反射响应存在误差，因此，我们引入基于最小二乘原理的保真项$\frac{1}{2}\|\boldsymbol{d}-\boldsymbol{W} \boldsymbol{i}\|_2^2$和L₁范数正则化项对式(9)进一步约束，构建出如下目标函数

$ J(\boldsymbol{i})=\frac{1}{2}\|\boldsymbol{d}-\boldsymbol{W} \boldsymbol{i}\|_2^2+\lambda\|\boldsymbol{i}\|_1 。$

(10)

式中，λ是惩罚参数，用于平衡保真项和正则化项，通过调节λ的值可以控制解的稀疏程度。

采用盲稀疏脉冲反褶积算法将多分量地震记录压缩为一系列尖脉冲，二维情况下的目标函数如下：

$ J_x\left(\boldsymbol{i}_x\right)=\frac{1}{2}\left\|\boldsymbol{d}_x-\boldsymbol{W i}_x\right\|_2^2+\lambda\left\|\boldsymbol{i}_x\right\|_1 \text {, } $

(11-a)

$ J_z\left(\boldsymbol{i}_z\right)=\frac{1}{2}\left\|\boldsymbol{d}_z-\boldsymbol{W i}_z\right\|_2^2+\lambda\left\|\boldsymbol{i}_z\right\|_1 。$

(11b)

式中：i_x、i_z分别表示x分量传播与反射响应和z分量传播与反射响应；d_x、d_z分别表示x和z分量地震记录。观测记录中的纵横波可以看作是震源子波与i(t)的褶积，我们采用快速迭代收缩阈值算法(FISTA)^[19]进行求解，提取i(t)。

3 盲稀疏脉冲反褶积目标函数求解

对式(10)所示目标函数中的保真项在i^(k)处做二阶泰勒展开

$ \begin{gathered} \underset{i}{\operatorname{argmin}} f(\boldsymbol{i})=\underset{i}{\operatorname{argmin}} f\left(\boldsymbol{i}^{(k)}\right)+\nabla f\left(\boldsymbol{i}^{(k)}\right)^{\mathrm{T}}\left(\boldsymbol{i}-\boldsymbol{i}^{(k)}\right)+ \\ \frac{1}{2}\left(\boldsymbol{i}-\boldsymbol{i}^{(k)}\right)^{\mathrm{T}} \nabla^2 f\left(\boldsymbol{i}^{(k)}\right)\left(\boldsymbol{i}-\boldsymbol{i}^{(k)}\right) 。\end{gathered} $

(12)

式(12)中，$\nabla^2 f(\boldsymbol{i})$与i无关，因此$\nabla^2 f\left(\boldsymbol{i}^{(k)}\right)$可表示为常数项$\frac{1}{t}$；常数项对于求解最小值没有影响，因此去除$f\left(\boldsymbol{i}^{(k)}\right)$，通过“配方法”得到

$ \underset{i}{\operatorname{argmin}} f(\boldsymbol{i})=\underset{i}{\operatorname{argmin}} \frac{1}{2 t}\left\|\boldsymbol{i}-\left[\boldsymbol{i}^{(k)}-t \nabla f\left(\boldsymbol{i}^{(k)}\right)\right]\right\|_2^2 。$

(13)

式(10)可以表示为

$ \boldsymbol{i}^{k+1}=\underset{i}{\operatorname{argmin}} \sum\limits_j \frac{1}{2 t}\left(\boldsymbol{i}_j-\left[\boldsymbol{i}_j^{(k)}-t \nabla f\left(\boldsymbol{i}_j^{(k)}\right)\right]\right)^2+\lambda\left|\boldsymbol{i}_j\right| 。$

(14)

令$\boldsymbol{p}\left(i_j\right)=\frac{1}{2 t}\left(\boldsymbol{i}_j-\left[\boldsymbol{i}_j^{(k)}-t \nabla f\left(\boldsymbol{i}_j^{(k)}\right)\right]\right)^2+ \lambda\left|\boldsymbol{i}_j\right|$，当$\frac{\partial \boldsymbol{p}}{\partial i_j}=0$时式(13)取得最小值，移项后有

$ \boldsymbol{i}^{(k+1)}=D\left(\boldsymbol{i}^{(k)}-t \nabla\left\|\boldsymbol{d}-\boldsymbol{W i}^{(k)}\right\|_2^2, \lambda * { step }\right) 。$

(15)

式中：D(·)为软阈值函数；$\text { step }=\frac{1}{2 \phi}$为FISTA算法每次迭代的步长；$\phi=\operatorname{eig}_{\max }\left(\boldsymbol{W}^{\mathrm{T}} * \boldsymbol{W}\right) ; \operatorname{eig}(\bullet)$表示矩阵的特征值。

4 数值算例 4.1 盲稀疏脉冲反褶积算法测试

随机生成如图 1(a)所示的稀疏序列作为PRR，与图 1(b)所示的主频为20 Hz雷克子波褶积得到合成记录如图 1(c)。

以图 1(c)作为输入，假定图 1(b)未知，测试本文算法求取i(t)的精度。将构建的子波褶积矩阵带入式(7)的目标函数中，即可利用FISTA算法提取得到PRR，图 2(a)为利用上述思路得到i(t)与准确解的对比。图 2(a)中求取i(t)所用的参数为λ=0.3，迭代次数为300次。由图 2(a)可见，本文算法求取的PRR与真实的反射系数相比较，相位信息准确，不存在时移误差，具有较好的保幅性。

利用求得的PRR与一个富含低频信息的子波褶积就可以得到低频信息丰富的地震信号。相较于雷克子波，高斯子波的频率分布较窄，具有良好的频率集中性，且高斯子波在低频部分的能量较高，与PRR褶积得到的地震信号能够更好的保留低频成分。图 2(b)为利用求得的PRR与一个主频20 Hz高斯子波褶积后的结果。图 2(c)展示了原地震信号与低频重构地震信号的频谱，显然图 2(b)中重构地震信号包含了更为丰富的低频信息，以这种信号作为EFWI的输入更有利于EFWI建立初始纵横波速度模型。

4.2 模型试算

纵横波速度模型如图 3(a)、3(b)所示，空间离散后的网格大小为10 m×10 m。模拟所用的观测系统如下：地表纵波源激发、炮间距50 m，炮点起始位置0 m，震源为主频8 Hz的Ricker子波(见图 3(c))，图 3(d)为其振幅谱，共58炮，全排列291个三分量检波器接收，采样间隔0.001 s，记录长度4 s。

图 3(e)、3(f)为正演得到炮点位于650 m处的x分量和z分量记录。对该记录进行高通滤波，滤除3 Hz以下的频率成分，得到缺失低频信息的两分量记录(见图 4(a)、4(b))，采用图 3(c)所示的子波对低频缺失两分量记录进行不依赖子波的弹性波全波形反演^[20-23]，反演所用的初始模型如图 4(c)、4(d)所示，反演结果如图 4(e)、4(f)所示，反演过程中选择最小偏移距道作为参考道^[24]，共迭代120次。可以看出，由于初始模型缺少先验信息，且合成记录中缺乏低频成分，不依赖子波的全波形反演出现速度错估，反演陷入局部极值，无法建立较高精度的纵横波速度模型。

基于相同的初始速度模型，进行盲稀疏脉冲反褶积的弹性波全波形反演。首先，利用盲稀疏脉冲反褶积算法，逐道提取合成记录的传播与反射响应，采用低频成分丰富的震源子波与PRR褶积得到低频丰富的重构地震记录，实现低频波场重建。

设置惩罚参数λ=1。从低频缺失的两分量合成记录中计算得到的传播和反射响应如图 5(a)、图 5(b)所示；为了验证传播与反射响应提取的正确性，将其与主频8 Hz高斯子波(见图 5(c)，图 5(d)为高斯子波的归一化振幅谱)褶积得到两分量低频重构记录如图 5(e)、图 5(f)所示，将重构记录与主频8 Hz高斯子波正演得到的合成记录对比，部分道的对比结果如图 5(g)、5(h)所示。图 5(i)展示了两分量低频缺失的合成记录和重构记录的频谱，由图可知，相较于原始合成记录，重构地震记录的低频成分更丰富。

采用巴特沃斯滤波器对图 5(e)、图 5(f)所示低频重构记录进行低通滤波, 截止频率为5 Hz, 为了缓解重构记录中振幅误差的影响, 强调信号相位的作用, 采用互相关目标函数计算伴随震源^[25]; 将滤波后的记录作为输入, 迭代20次得到低频重构记录的时间域全波形反演结果(见图 6), 纵横波速度模型中大尺度信息被正确重建, 说明本文方法可以在低频缺失情况下建立良好的初始模型。

将低频缺失的合成记录作为输入，以图 6(a)、6(b)为初始纵横波速度模型，并选择每一炮的最小偏移距道作为参考道，以图 3(c)所示子波作为震源进行基于褶积法的不依赖子波弹性波全波形反演，共迭代120次，得到最终反演结果如图 7(a)、7(b)所示，图 7(c)和图 7(d)为模型位置1 300 m处初始速度、反演速度与真实速度的对比图。

对比图 7(a)、图 7(b)和图 4(e)、图 4(f)可知，相较于地震记录低频缺失情况下不依赖子波的全波形反演算法，基于盲稀疏脉冲反褶积的弹性波全波形反演结合不依赖子波的全波形反演没有出现速度错估。这说明基于盲稀疏脉冲反褶积的全波形反演提供的初始模型可以有效缓解“周波跳跃”现象，有利于反演的稳定迭代。

4.3 联合低频重建与多尺度反演方法

全波形反演是一个强非线性问题，Born近似的散射项是阻碍反演正确收敛的重要原因，但这种负面影响随着震源频率的降低而减弱，此外，低频信号具有更长的波长，能够有效避免“周波跳跃”现象。为了缓解全波形反演的强非线性问题，分频带的多尺度反演策略被提出^[2]，该方法要求在反演前对地震信号进行频带划分，由低频段向高频段逐级反演。为了分析本文方法得到的初始模型对多尺度反演的贡献，以图 3(a)、3(b)为准确模型，以图 4(c)、4(d)所示的线性速度模型作为初始模型，反演参数和观测系统的布置同4.2节所述，震源采用主频为12 Hz的雷克子波，对比分析如下三种反演方法：(1)使用线性初始模型进行常规全波形反演；(2)以低频重构记录反演得到的速度模型代替线性初始模型，进行后续反演；(3)以低频重构记录反演得到的速度模型代替线性速度模型，进行多尺度全波形反演。

图 8展示了直接以线性速度模型为初始模型进行常规全波形反演的结果，反演结果与真实速度模型差距较大，这主要是由于线性初始速度模型中不含有构造信息，震源子波主频较高，常规时间域全波形反演陷入局部极值造成的。以线性速度模型为初始模型，以本文方法得到的低频重构记录作为输入，反演结果如图 9所示，相比常规全波形反演结果，速度模型中大尺度构造信息得到正确恢复，但精细地质构造没有被修正。采用巴特沃斯滤波器将合成记录分为三个频段(0~5、5~9、9~12 Hz)进行多尺度全波形反演，每个频段迭代30次。以线性速度模型为初始模型的多尺度全波形反演结果如图 10所示，以图 9所示速度模型为初始模型进行多尺度全波形反演的反演结果如图 11。对比图 10和图 11，由于原地震数据中缺乏低频信息，常规多尺度全波形反演缺失长波长模型更新信息，不能缓解FWI对初始模型的依赖，反演精度较差，而以图 9所示速度模型为初始模型的多尺度全波形反演实现了较高精度的纵横波速度模型重建，本文方法能够从低频缺失的地震数据中获取地下介质模型的长波长信息，联合多尺度全波形反演有效降低了全波形反演的非线性，避免落入局部极值，实现了纵横波速度模型的高精度重建。

5 结语

本文从压缩感知理论和稀疏反演方法出发，利用地震反射系数的稀疏性构造稀疏约束反演的目标函数，实现了L₁范数约束下的盲反褶积算法。该算法将地震信号压缩为宽频反射脉冲，并通过与已知的富含低频成分的子波褶积得到低频重构波场，实现了低频缺失地震记录的低频补偿，由于褶积的低频子波已知，可直接用于全波形反演算法中，实现低频缺失和子波不准情况下全波形反演初始模型的有效构建。模型测试结果表明，在子波不准和地震记录低频缺失情况下，基于盲稀疏脉冲反褶积的全波形反演算法能够建立良好的低频初始模型，并且避免了子波差异对反演的影响，有效缓解了“周波跳跃”现象，降低弹性波全波形反演陷入局部极值的可能。