0 引言

随着油气勘探、开发的不断深入，地震资料分辨率与储层预测精度的矛盾日益突出。迫切需要提高地震资料的分辨率和成像精度，依托高分辨率地震资料对储层的结构特征和沉积层序进行综合解释和精细描述。

业界广泛应用的提高地震资料分辨率的传统方法，如地表一致性反褶积、谱模拟反褶积和谱白化反褶积等，是基于线性维纳滤波器，存在反射系数白噪或地震子波最小相位的假设条件。由于地震资料的频带是带限的，这些反褶积方法虽然在一定程度上能拓宽地震记录的频谱，但并未突破频带限制，不能恢复反射系数的全频带特征。

非线性反褶积方法摆脱了子波最小相位和反射系数白噪的假设，一定程度上突破有效频带的限制，更大幅度地提高地震记录分辨率。Wiggins^[1]提出了最小熵反褶积，在“稀疏度”准则下实现了混合相位子波反褶积处理；Sacchi^[2]和Ulrych等^[3]提出利用柯西准则进行非线性稀疏反褶积，通过重加权迭代方法恢复稀疏的反射系数序列；Velis^[4]、云美厚等^[5]提出随机稀疏脉冲反褶积，通过检测稀疏反射系数的位置和大小提高地震记录的分辨率；潘树林等^[6]提出了自适应步长(FISTA)算法稀疏脉冲反褶积，以提高目标函数的收敛性。然而这些非线性反褶积方法都是基于单道运算的，每道的结果只与本道有关。由于忽视了道与道之间的空间结构联系，单道反演结果容易产生横向不连续性问题。为了克服反褶积结果的横向连续性差的问题，许多学者进行了持续研究和探索，近年还提出了通过引入横向约束项的多道反褶积或反演方法。

多道反褶积或反演方法可分为两类。一类是基于模型驱动的，如Yuan等^[7]在反演过程中通过对地震资料求取纵向和横向一阶导数作为稀疏约束项，以压制随机噪声和孤立噪声，得到较好去噪效果；Zhang等^[8]提出具有横向约束的多道基追踪反演法，在横向约束上采用一个“Z”形的空间导数算子，以提高反演结果的横向一致性。这类多道反演法比传统单道反演方法更多地利用了空间结构信息，因此能得到更好的反演结果。然而模型驱动的多道反演方法的横向约束算子是固定的，反演结果易受模型结构的影响，不同的模型得到的反演结果差别很大。

另一类是基于数据驱动的多道反演方法。如Clapp等^[9]从偏移结果中估算倾角信息得到一个非稳态的各向异性约束算子，并将其引入反射波层析成像中，可得到更连续的层析成像结果；Kazemi等^[10]提出稀疏多道盲反褶积方法，通过消除子波的影响，得到道与道之间的关系，再进行多道反演；Yuan等^[11]提出变换域波阻抗多道反演方法，通过在降维地震数据上强加一个约束项实现同时多道波阻抗反演，解决了单道反演方法的横向不连续性问题；Hamid等^[12]提出基于地质结构约束的波阻抗反演方法，在模型或数据空间内估算各位置的地质结构倾角作为横向约束算子，提高了反演的连续性；Yuan等^[13]提出一种稳定的基于反演的多道补偿方法，通过求取f-k域反射系数的目标函数的最小值展开多道同时反演，相对于常规单道方法能更好地保护同相轴的横向连续性；Nose-Filho等^[14]提出基于最小熵反褶积的改进SMBD方法的快速算法，通过设计反褶积滤波器使归一化的混合L₁和L₂范数的目标函数取得最小值，该方法计算量较小，效率高。

数据驱动的多道反演方法是从地震数据本身出发，利用地震道与道之间的关系，求取一个横向约束算子或将地震数据变换到某个域求取横向约束算子，该算子具有自适应性，不会随地震数据而变化，因此反演结果更稳定、可靠。

本文引入构造导向滤波算子作为约束条件，完全基于数据驱动，利用反演思想实现多道反褶积，在提高地震记录分辨率的同时，显著地改善了地震资料的空间连续性。

1 方法原理 1.1 稀疏反褶积方法

地震信号可表示为地震子波与反射系数的褶积

$ s(t)=w(t) * r(t) $

(1)

式中：w(t)为地震子波；r(t)为反射系数序列；s(t)为合成的无噪地震记录。

实际地震资料往往包含很多噪声，可表示为

$ d(t)=s(t)+n(t) $

(2)

式中：d(t)为实际的地震记录；n(t)为随机噪声。

常规稀疏反褶积方法假设反射系数是稀疏的脉冲序列，则稀疏反褶积可看成是一个稀疏序列的恢复问题。通过引入稀疏正则化约束项，在贝叶斯数学框架下，利用最小二乘法求取稀疏的反射系数。根据贝叶斯原理，稀疏反褶积目标函数为

$ J=\|\boldsymbol{W} \boldsymbol{r}-\boldsymbol{d}\|_2^2+\mu_1\|\boldsymbol{r}\|_1 $

(3)

式中：W为子波褶积矩阵；μ₁为正则化参数。第二项为稀疏约束项，可选择柯西约束、Sech分布、Huber分布、修正柯西约束等。本文选用修正柯西约束，在反射系数较小时，修正柯西约束具有较好的识别和保护弱反射系数能力^[15]。

目标函数式(3)的优化问题是非线性的，可通过共轭梯度法、重加权迭代法求解。

1.2 构造导向滤波算子构建

构造导向滤波算子构建的关键是如何计算地震波同相轴的局部倾角属性。为避免噪声影响，应用Fomel等^[16]提出的希尔伯特变换代替导数运算求取局部倾角，具体表示为

$ \sigma=-\frac{H_{\mathrm{HT} x}}{H_{\mathrm{HT} y}+\xi} $

(4)

式中：H_HTx、H_HTy分别为希尔伯特变换时间方向和空间方向的分量；ξ为极小非零常数。在实际计算过程中，通过求取希尔伯特变换时间方向和空间方向的分量，得到地震波同向轴的局部倾角属性，可减小直接求导带来的高频噪声的放大效应，同时该方法运算效率较高。

得到局部倾角σ之后，可计算出每一道的预测因子Q_{i, j}(σ_i)，进而得到构造导向滤波算子T=1/D。其中D可表示为

$ \boldsymbol{D}=\left[\begin{array}{cccc} \boldsymbol{I} & 0 & \cdots & 0 \\ -Q_{1, 2}\left(\sigma_1\right) & & & \\ \cdots & \cdots & \cdots & \cdots \\ 0 & 0 & -Q_{N-1, N}\left(\sigma_{N-1}\right) & \boldsymbol{I} \end{array}\right] $

(5)

式中：I为单位矩阵；Q_{i, j}(σ_i)为第i道~第j道的预测因子；σ_i为局部倾角。

1.3 横向约束多道反褶积

在求得构造导向滤波算子T后，将其作为横向约束算子加入反演的目标函数中，同时纵向约束选用修正柯西约束，得到多道反褶积方法的最终目标函数

$ J=\|\boldsymbol{W r}-\boldsymbol{d}\|_2^2+\mu_1\|\boldsymbol{r}\|_1+\mu_2\|\boldsymbol{T} \boldsymbol{r}\|_2^2 $

(6)

这里μ₁、μ₂均为稀疏约束项正则化参数，分别用于提高纵向分辨率和控制空间横向连续性。通过引入横向约束项，使每一道反演结果都与其附近地震道有关，从而改善反演结果的横向连续性，并提高其横向分辨率。

本文使用分裂Bregman优化算法求解式(6)，该迭代求解过程分为三个步骤。

(1) 利用共轭梯度法求解线性问题

$ \begin{aligned} & \boldsymbol{r}^{k+1}=\underset{m}{\operatorname{argmin}}\|\boldsymbol{d}-\boldsymbol{W} \boldsymbol{r}\|_2^2+\frac{\lambda}{2}\left\|\boldsymbol{x}^k-\boldsymbol{m}-\boldsymbol{b}_x^k\right\|_2^2+ \\ & \mu_x\|\boldsymbol{T} \boldsymbol{r} \boldsymbol{\|}\|_2^2 \end{aligned} $

(7)

(2) 利用快速软阈值算法求解非线性问题

$ \boldsymbol{x}^{k+1}=\underset{x}{\operatorname{argmin}} \mu_z\|x\|_1+\frac{\lambda}{2}\left\|\boldsymbol{x}-\boldsymbol{r}^{k+1}-\boldsymbol{b}_x^k\right\|_2^2 $

(8)

(3) 更新辅助变量

$ \boldsymbol{b}_x^{k+1}=\boldsymbol{b}_x^k+\boldsymbol{r}^{k+1}-\boldsymbol{x}^{k+1} $

(9)

式中：k为迭代次数；变量x^k和b_x^k的初始值为零向量；λ为正则化因子，一般设为1。

在进行横向约束多道反褶积时，时窗可选择包含目的层范围且上、下各取10 ms。在三维情况下，可采取逐线方式进行，并采取多核并行计算以提高计算效率^[17-18]。

2 模型试验

设计图 1a所示的反射系数模型，将其与主频为30 Hz、采样间隔为1 ms的零相位雷克子波褶积，得到合成的无噪地震记录(图 1b)。对该合成记录加入信噪比(SNR)为5的随机噪声(图 2a)，得到该含噪记录的修正柯西约束单道反褶积结果(图 2b)。再对含噪记录求取构造导向算子(图 3a)，并将其作为横向约束加入多道反褶积中，得到反褶积结果(图 3b)。对比可知，在第三薄层及下伏倾斜地层中，多道反褶积不仅纵向分辨率比单道反褶积高，且表现出较强横向连续性；同时在楔形体尖灭地层，多道反褶积的分辨率较单道反褶积有明显改善。

另外，在抗噪性方面，修正柯西约束单道反褶积在含噪情形下会去掉一部分有效信号，并且同相轴出现了抖动且在横向上不连续，而本文方法(多道反褶积)所得结果呈现较干净且横向连续性很好的面貌，且抗噪性更强。

3 实际资料处理

松辽盆地M区块实际地震资料(图 4a)中优势频带范围是15~45 Hz(图 5中黑线)，利用本文基于构造导向滤波的多道反褶积方法处理后得到结果剖面(图 4b)，相应的频谱分析结果(图 5中红线)显示其优势频带范围扩展为15~65 Hz。同时，从图 4中区内两口井的井震标定结果可知，通过本文方法处理后波组关系更准确(图 4左侧圆圈所示)，井震标定更一致(图 4右侧圆圈所示)，相关系数由0.53升至约0.88。

进一步，针对三维叠后数据体通过逐线处理方式得到多道反褶积后的数据体，处理过程中采用GPU提高运算效率。图 6是对图 4地震资料提取的沿层切片。通过对比可看出，在引入构造导向滤波约束后再进行多道反褶积，可在较充分地保持空间连续性的情况下，大幅度地提高地震数据横向分辨率，其小构造更明显，断层更清晰，更有利于后续的精细储层预测。

4 结论

(1) 现有的多道非线性反褶积方法主要分为模型驱动和数据驱动两大类。数据驱动类方法完全从数据本身出发，比模型驱动类方法具有更强适应性，因此反褶积结果常优于模型驱动类方法。

(2) 本文提出的基于构造导向滤波的多道反褶积方法是一种数据驱动的多道非线性反褶积方法，它直接从数据中求出构造导向滤波算子作为横向约束项，不仅可提高地震资料分辨率，还能有效提高地震剖面的横向连续性。

(3) 本文方法通过逐线进行三维地震数据反演，效率有待提高；同时，后期需探究并完善三维反演目标函数约束下的多道反褶积方法，以更显著提高反演的效率。