0 引言

在油气地震勘探中，经常会遇到储层被强反射屏蔽，从而造成储层地震响应特征不明显，振幅属性难以表征储层的发育情况。目前去屏蔽技术多集中在子波分解与重构技术，主要是利用不同频率、振幅、相位的各类原子与地震道匹配^[1-3]，分解成1个原子库，认为最大能量原子对应强背景，在重构中舍去，达到去屏蔽效果。这类技术存在2个问题：一是地震道子波分解在单道进行，无法考虑地层横向关系，容易造成匹配过度和重构剖面连续性差；二是无法确定舍去哪些原子，其地质和地球物理意义不明确^[4]。另外，该技术对计算时窗要求较高，一般是在1个反射轴内进行分解和重构^[5-7]。如1个波峰或波谷、波谷到波谷(波峰到波峰)等，只适用于屏蔽层厚度横向变化比较稳定的地质情况，如松辽盆地T₂页岩层^[3]、鄂尔多斯盆地山西组发育的煤层等地质条件。当屏蔽层厚度横向变化大时会造成欠匹配或匹配过度的情况，影响预测精度。

渤海湾盆地济阳坳陷渤南地区砂砾岩致密储层受膏岩高速体和自身复杂的叠置关系以及非均质性等影响，致使该类储层具有强屏蔽弱反射的响应特征^[8-9]，并且受构造和沉积因素影响，该区膏岩层表现为局部分布，且厚度变化较大^[10-11]。基于子波分解技术的去强反射方法不适用。因此需要开展基于提高分辨率方法的去强反射技术。

针对渤南地区实际资料的强反射屏蔽问题，系统开展基于“钉型”子波+压缩感知技术的方法研究，首先对原始地震数据进行“钉型”子波谱整形处理，实现对子波旁瓣的压制，在此基础上，再进行压缩感知处理，达到去除膏岩层强屏蔽，突出目的层段砂岩的地震响应特征的效果，以期提高砂岩的预测能力，为下一步勘探部署提供支撑。

1 原理与方法 1.1 “钉型”子波谱整形的原理

地震记录为某一主频不变的子波与地层反射系数的褶积，子波的形状对地震分辨率起着关键的作用，雷克(Ricker)子波[图 1(a)]是零相位子波，有1个主峰值和2个小的旁瓣。雷克子波仅由信号参数f确定，从子波频谱上可以看到，f是其峰值频率。带通子波[图 1(b)]也是零相位子波，但其实际上定义的是1个滤波器，带通子波有多个旁瓣。从子波频谱上可以看到，带通子波形态需由4个频率确定，f₁为低截频率，f₂为低通频率，f₃为高通频率，f₄为高截频率。“钉型”子波为改进的宽带Butterworth子波^[12-13]，“钉型”子波具有旁瓣小、周期短、宽频的特点(图 2)

改进的“钉型”子波的频率域表达式为

$ {A^2}(f) = \left\{ {\begin{array}{*{20}{l}} {\frac{{\frac{{{f^{2N}}}}{{{f_{\rm{a}}}}}}}{{{{\left( {1 + \frac{f}{{{f_{\rm{a}}}}}} \right)}^{2N}}{{\left( {1 + \frac{f}{{{f_{\rm{b}}}}}} \right)}^{2M}}}}, }&{f \le {f_0}}\\ {\frac{{\frac{{f_0^{2N}}}{{{f_{\rm{a}}}}}}}{{{{\left( {1 + \frac{f}{{{f_{\rm{a}}}}}} \right)}^{2N}}{{\left( {1 + \frac{{{f_0}}}{{{f_{\rm{b}}}}}} \right)}^{2M}}}}, }&{{f_0} < f \le {f_{\rm{b}}} - {A_{\rm{w}}}{\rm{d}}f}\\ {\frac{{\frac{{f_0^{2N}}}{{{f_{\rm{a}}}}}}}{{{{\left( {1 + \frac{{{f_0}}}{{{f_{\rm{a}}}}}} \right)}^{2N}}{{\left( {1 + \frac{{{f_0}}}{{{f_{\rm{b}}}}}} \right)}^{2M}}}} \cdot \cos \left[ {{\rm{ \mathsf{ π} }}\frac{{f - \left( {{f_{\rm{b}}} - {A_{\rm{w}}}{\rm{d}}f} \right)}}{{{A_{\rm{w}}}{\rm{d}}f}}} \right], }&{{f_{\rm{b}}} - {A_{\rm{w}}}{\rm{d}}f < f \le {f_{\rm{b}}}} \end{array}} \right. $

(1)

式中：N = SL/6，M = SH/6，SL和SH分别对应原始地震数据频谱低频部分和高频部分倍频衰减的斜率；f₀为子波主频，Hz；f_a为低截频率，Hz；f_b为高截频率，Hz；df为频率间隔，Hz；A_w为频率域“钉型”子波载波调制因子长度，即原始地震数据频谱高通频率到高频截止频率的长度，其值大小决定高频率成分的多少。

采用雷克子波和“钉型”子波的正演合成记录(图 3)，下伏层砂岩储层受膏岩层强屏蔽的影响，地震响应微弱，“钉型”子波比雷克子波有更强的地震反射，有利于突出强背景下的弱反射。

1.2 压缩感知原理

压缩感知技术利用信号的稀疏特征，可以把高维空间的信号通过测量矩阵投影到一个低维的空间中，通过非线性重构来完美重建信号。首先通过理论推导构造基于压缩感知原理的求解地震反射系数的目标函数，然后采用快速迭代软阈值算法进行求解^[14]，通过去除反演的强反射系数界面，最终得到去强轴后的地震记录。

根据地震褶积模型，y为地震数据，可以描述为地震子波w和地下反射系数x的褶积

$y=w \otimes x+n$

(2)

式中：n为噪音。

式(2)在频率域可以表示为

$Y=W X+N=W F x+N$

(3)

式中：Y为地震记录的Fourier变换；W为地震子波的Fourier变换；X为地下反射系数的Fourier变换；N为噪声n的Fourier变换；F是Fourier变换矩阵。

根据地震勘探基本假设，地下反射系数是随机序列，则其频谱X应是全带宽的，但是由于地震子波的滤波作用使其成为有限带宽的频谱Y，从而损失了许多有用信息。因此，首先需要通过有限带宽的地震数据频谱Y来恢复整个带宽数据X。

压缩感知利用信号的稀疏特性，通过非线性重构来完美重建信号。由于地震反射系数x是稀疏的，所以求解式(3)的过程满足压缩感知理论。与传统常用的压缩感知方法不同的是，这里的采样矩阵W不是完全随机函数。所以只能在一定程度上恢复全带宽能量Y，而不能完全恢复整个频带的能量。

由式(2)和式(3)可知，对于地震记录，A = W F。根据压缩感知原理，式(3)求解可构造如下成本函数并使其为最小

$ \frac{1}{2}\|y-A x\|_{2}^{2}+\lambda\|x\|_{1} \Rightarrow \min $

(4)

式中：||·||₂²和||·||₁分别为L2范数和L1范数；λ是L2范数和L1范数权重的调节因子，λ越大，L1范数占的权重就越大。

式(4)的前半部分采用L2范数约束，是为了在最小平方意义下匹配地震记录频谱。但因为方程(3)是欠定的，最小平方意义下存在多解性；为了克服多解性，引入成本函数第二部分的L1范数约束条件，这也是利用压缩感知原理恢复信号的关键之一。

对式(4)求偏导数并令其为零，可以得到目标函数

$ {\mathit{\boldsymbol{A}}^H}\mathit{\boldsymbol{A}}x - {\mathit{\boldsymbol{A}}^H}y + \lambda {\mathop{\rm sgn}} (x) = 0 $

(5)

令z=A^Hy，B=A^HA，上式简化为

$ z=\boldsymbol{B} x+\lambda \operatorname{sgn}(x) $

(6)

式(6)可以采用快速迭代软软阈值算法(FISTA)求解，具体流程如下：

步骤1：初始化z₁ = x₀，t₁ = 1；

步骤2：$ x_{k}=\operatorname{soft}\left[z_{k}+\frac{1}{\alpha(W F)^{H}\left(y-W F z_{k}\right)}, \frac{\lambda}{\alpha}\right]$；

步骤3：$ t_{(k+1)}=\frac{1+\sqrt{1+4 t_{k}^{2}}}{2}$；

步骤4：$ z_{k+1}=x_{k}+\frac{t_{k-1}}{t_{k+1}}\left(x_{k}-x_{k-1}\right)$。

用步骤2计算x_{k + 1}，并计算其与x_k的误差，重复步骤2—4直到满足精度要求。

其中：α须大于矩阵B的最大特征值；soft为迭代软阈值函数，定义为

$ \operatorname{soft}(u, \alpha)=\operatorname{sgn}(u) \max (|u|-\alpha, 0) $

(7)

1.3 基于钉型子波谱整形和压缩感知去除强反射界面

为保证去强轴反射界面的横向连续性，首先根据前文所述原理方法，对地震数据整道进行压缩感知反演，恢复得到全频带反射系数x；然后设定目标强反射界面层位处2~3个地震波形反射长度的时窗，沿强屏蔽层上下各半个时窗长度，采用余弦窗函数对整道反射系数进行截取。采用余弦窗函数的目的一是保证强屏蔽层反射系数的准确截取；二是避免由于时窗截取导致时窗外地震数据的突变导致的横向连续性变差；其次从整道反演的反射系数中减去截取的强反射系数s，得到去除强反射系数后的反射系数剖面；最后将设计的“钉型”子波褶积去除强反射系数后的反射系数，得到最终去强轴的地震记录。余弦窗函数截取反射系数为

$ s = \cos (n\;win\;{\rm{ \mathsf{ π} }})x $

(8)

式中：s为截取的强反射界面反射系数；x为反演的全频带反射系数；n win为选取时窗。

从采用雷克子波、“钉型”子波和压缩感知去强轴技术组合后的正演记录(图 4)，可以看出，在分辨率足够高的前提下采用“钉型”子波(压制旁瓣) +压缩感知去强轴技术组合可以有效地识别屏蔽储层反射。

2 实际资料应用效果 2.1 研究区地震地质条件

研究区位于渤海湾盆地济阳坳陷渤南地区，主力含油层系为古近系沙河街组沙四上亚段^[15-16]。其地层厚约300 m，分为4个砂层组，1与2砂组发育膏岩，地震表现为强反射特征；3与4砂组为灰色粉细砂岩，受上部膏岩层的屏蔽影响，砂岩储层地震响应特征不清[图 5(a)]。

2.2 去屏蔽处理效果

对原始地震数据进行“钉型”子波谱整形处理，通过对子波旁瓣的压制后得到的地震频谱[图 5(d)]，与原始地震频谱[图 5(b)]相比，地震分辨率得到了一定的提高，其中时间分辨率从30.0 ms缩短为24.8 ms，周期更短，时间分辨率更高，地震主频从14.5 Hz提高到17.5 Hz，绝对带宽从13.0 Hz提高到19.0 Hz，相对带宽从1.5 dB提高到1.8 dB。在此基础上，再进行压缩感知处理，地震分辨率得到进一步的提高，其中时间分辨率、地震主频和绝对带宽都有不同程度的提高[图 5(f)]，目的层段的砂岩地震响应特征更加明显，达到了去强屏蔽的目的。

对比分析去屏蔽前后的地震剖面，并选取研究区2口实钻井A和B的测井岩性资料，分别进行岩性地震标定[图 5(c)，(e)]，可以看出，基于“钉型”子波+压缩感知技术处理后，能有效突出沙四上亚段3与4砂组的地震反射特征，提高了储层的分辨能力，且与测井岩性标定结果吻合，说明基于“钉型”子波+压缩感知技术能够达到有效去除强屏蔽的目的，提高地震勘探对强屏蔽弱反射砂体的识别能力。

如图 6所示，原始振幅切片[图 6(a)]受膏岩层屏蔽影响，导致下伏地层的振幅切片难以反映沉积相和砂体的展布特征。“钉型”子波谱整形后的振幅切片[图 6(b)]得到了一定的改善，但还难以达到去屏蔽的效果，而基于“钉型”子波谱整形基础之上的压缩感知去屏蔽处理后的地层切片[图 6(c)]能够有效反映砂体的展布特征，其显示出的地震地貌特征符合井间沉积特征(扇三角洲前缘沉积形态及北东向和由北向东的物源特征) ^[17-19]。与测井砂地比模式图[图 6(d)]展布特征一致，表现为3条从北东向和由北向东展开的扇三角洲前缘沉积体^[20-22]。

3 结论

(1) 常规基于子波分解与重构技术的去屏蔽方法，由于对分解原子重构时窗的要求严格，只适用于横向稳定分布且厚度变化不大的储层类型，对于横向变化较大的砂砾岩等储层类型无法考虑地层横向关系，容易造成匹配过度和重构剖面连续性差等问题。

(2) “钉型”子波+压缩感知技术是一套综合分析流程，通过提高分辨率达到去除强反射的目的，优于常规子波分解去强屏蔽方法。

(3) 对渤海湾盆地济阳坳陷渤南地区实际资料应用“钉型”子波+压缩感知技术处理，有效地去除了强反射背景，突出了目的层段砂砾岩的地震响应特征，提高了砂砾岩的预测能力，为下一步勘探部署提供了有力支撑。