0 引言

随着全球油气勘探程度地不断深入，薄层和薄互层油藏的重要性日益增加，提高薄储层识别能力对油气勘探具有重要意义^[1]。对于储层精细勘探开发，大规模水平井部署已成为必要，这对地震薄储层预测又提出了更高要求^[2]。由于薄储层顶底界面的地震反射波相互干涉，随着储层厚度减小，基于地震资料开展薄层识别的难度将急剧增加。

薄储层通常具有横向变化快、分布规律复杂等特征。为搞清薄储层的空间分布规律，国内外学者作了大量相关研究，其中刘化清等^[3]立足砂泥岩薄互层地质背景，提出了一套基于地震沉积学的薄砂体预测方法，利用薄层沉积体的平面分布尺度远大于其纵向厚度的普遍规律开展薄储层平面分布预测；杨占龙等^[4]则提出了一种联合井-震时频匹配分析与地震全反射追踪的地震隐性层序界面识别及高频层序格架建立方法，大幅增强了薄油气层系的预测能力；刘恭利等^[5]针对薄互层火成岩厚度不均匀、横向变化大等特征，开展了薄互层形式下火成岩地震响应特征多因素影响分析及厚度预测；李亚哲等^[6]则针对砂砾岩开展了基于地震波形指示反演的薄储层预测研究，实现了砂砾岩储层的定量化预测。目前薄储层预测依然受到诸多因素地影响，如地震分辨率不足、模型依赖严重、横向随机性强、井资料要求高等^[7-10]。

高精度地层波阻抗反演对预测薄储层分布具有重要意义^[1]。波阻抗反演通过将层间地震反射信息转换为地层波阻抗信息，从而进一步突出加强地震资料的地质意义，有利于储层精细分析；特别是在无井或少井区域，通过波阻抗反演预测储层与围岩的波阻抗差异对于油气勘探与开发尤为重要。常用的基于模型的波阻抗反演方法需在前期根据地震、测井等多源信息构建目的层低频波阻抗模型，然后在此模型基础上开展扰动进而逼近实际地层波阻抗^[11]，这一过程对模型数据较为依赖，其反演结果受初始模型精度影响较大。对于地层波阻抗信息，亦可基于地震资料开展地层反射系数估值，进而在此基础上计算地层相对波阻抗信息^[12-15]，该过程无需构建低频波阻抗模型，其最终结果不受模型精度影响，可在无井或少井区域开展应用。由于叠后地震数据通常可视为地层反射系数和带限地震子波的褶积结果，通过地震资料获取地层反射系数可大幅降低带限子波的影响从而提高资料薄层识别能力^[16]。地震反射系数反演主要基于时域褶积和频谱分解技术拓宽资料有效频带，进而提高调谐厚度之下的薄地层成像精度。研究表明^[17-25]，通过压制带限子波对地震分辨率的影响，地震反射系数反演能够突破常规的地震分辨极限（1/4波长），可识别远小于地震调谐厚度的薄储层^[2]，尤其是新的地震反射系数反演理论的引入，对当前中国陆相盆地大量存在的厚度3 m左右的薄储层勘探具有重要意义。该方法近年来已被广泛应用于薄层探测和地层界面识别^[26]。

由于薄层地震响应关系复杂，地震反射系数反演须附加约束条件以降低反演的多解性。其中，最小二乘约束作为降低地震反射系数反演多解性的有效工具已被广泛应用^[19]。作为一种光滑约束，基于最小二乘约束的地震反射系数反演在假设反射系数频带为白谱的条件下采用L2范数约束反演结果^[27]。由于地震频谱的带限性，最小二乘约束下的地震反射系数反演仍存在较强的多解性。为此，研究人员^{[20-21, 28]}引入了稀疏约束（L1范数、Geman-McClure范数、熵范数等）以进一步降低反演多解性。其中，L1范数在地震稀疏反射系数反演中应用最为广泛，其可通过基追踪理论^[29]求解。此外，相应的稀疏约束还包括预设稀疏脉冲数量^[22]和Lq（0 < q < 1）范数^[23]等，然而，目前没有理论证明上述任何一种约束条件对地震反射系数反演的约束效果最佳^{[19, 23]}。

作为一种重要的回归和分类算法，稀疏贝叶斯学习理论^[30]近年来亦被引入地震稀疏反射系数反演^{[17, 19, 24-25]}。与传统稀疏表征理论不同，稀疏贝叶斯学习在贝叶斯框架下通过边缘概率最大化获取稀疏表征结果。当前主流的稀疏贝叶斯学习理论主要包括基于序贯算法的稀疏贝叶斯学习（SBLSA）^[31]和基于最大期望算法（EM）的稀疏贝叶斯学习（SBL-EM）^[32]，二者差异主要体现在边缘概率最大化算法的不同。在基于SBL-SA理论的地震反射系数反演中，通过在每次迭代中增加、删除或重估相应的稀疏脉冲以提高边缘概率最大化程度^{[17, 31]}，进而获得地震数据对应的稀疏反射系数。研究表明^[31]，SBL-SA理论存在局部收敛等问题，会降低地震稀疏反射系数反演的精度。为此，Yuan等^[19]提出了一种基于SBL-EM的地震稀疏反射系数反演方法，其通过经典的EM算法而非序贯算法开展边缘概率最大化计算，大幅度提高了地层稀疏反射系数的表征能力。与传统SBL-SA理论相比，SBL-EM能够有效提高薄储层预测精度，在同等资料品质条件下薄层识别能力较传统方法可提高约一倍。本次研究在基于SBL-EM地震反射系数反演的基础上联合设计了线性FIR滤波器来开展目的层相对波阻抗预测，以期能够大幅度提高薄层识别精度，并能为当前我国广泛存在的3 m左右的薄储层勘探奠定基础。

1 地震反射系数反演 1.1 稀疏贝叶斯学习理论

稀疏贝叶斯学习是解决回归和分类问题的有效方法^[30]，目前已在地震反射系数反演领域得到广泛应用^{[17, 19, 24-25]}。设地层反射系数序列为r (t)，且满足零均值高斯分布，方差为ε (t)，其中t表示时间，则r对应的高斯先验分布可写为

$ p(r \mid \varepsilon)=(2 \pi)^{-\frac{L}{2}} \prod\limits_{i=1}^{L} \varepsilon_{i}^{-\frac{1}{2}} \exp \left(-\frac{r\left(t_{i}\right)^{2}}{2 \varepsilon_{i}}\right) $

(1)

式中：p(.|.)为概率分布，L为r (t)对应的样点数。

设数据噪声e(t)的概率分布特征与位置t无关且满足零均值高斯分布，高斯方差为常数χ。则对应的似然函数可写为

$ $

(2)

式中：ψ = Dr + e，D为正演矩阵；H为向量ψ的长度；T为矩阵转置；Φ = χI，I为单位阵。

根据贝叶斯理论^[33]，稀疏反射系数的后验分布可表示为

$ \begin{aligned} p(r \mid \psi, \varepsilon, \chi) &=\frac{p(\psi \mid r, \chi) p(r \mid \varepsilon)}{\int p(\psi \mid r, \chi) p(r \mid \varepsilon) \mathrm{d} r} \\ &=(2 \pi)^{-\frac{L}{2}}\left|\sum\right|^{-\frac{1}{2}} \exp \left[-\frac{1}{2}(r-\mu)^{T} \Sigma^{-1}(r-\mu)\right] \end{aligned} $

(3)

其高斯期望μ和方差∑分别为

$ \left\{\begin{array}{l} \sum=\left(\Gamma^{-1}+\chi^{-1} D^{T} D\right)^{-1} \\ \mu=\chi^{-1} \sum D^{T} \psi \end{array}\right. $

(4)

式中：Γ = diag (ε)。由上式可知，反演期望μ是超参数[χ, ε]的函数，而[χ, ε]可在稀疏贝叶斯学习框架下通过边缘概率最大化的方式求取。设对应的边缘概率密度函数为

$ p(\psi \mid \Upsilon)=(2 \pi)^{-\frac{H}{2}}|\Theta|^{-\frac{1}{2}} \exp \left(-\frac{1}{2} \psi^{T} \Theta^{-1} \psi\right) $

(5)

式中：ϒ = [ε, χ]且Θ = χI + DΓD^T。

当[ε, χ]等于真实值时，p(ψ| ϒ)应为最大，对此可通过迭代更新ϒ使得p(ψ| ϒ)取得极大值；当获取最优ϒ_ML = [ε_ML, χ_ML]后，稀疏反射系数r (t)可表征为

$ r(t)=\left(D^{T} D+\chi_{M L} \Gamma_{M L}^{-1}\right)^{-1} D^{T} \psi $

(6)

为使边缘概率密度函数p(ψ| ϒ)取得最大值，等同于取-log p(ψ| ϒ)的最小值，其中-logp(ψ| ϒ) = log| Θ | + ψ^TΘ^-1Ψ。由于-logp(ψ| ϒ)表达式中不含r (t)，因此引入EM算法设r为隐变量，并计算下式的极大值

$ \begin{array}{l} Q\left( \Upsilon \right) = {E_{r\mid \psi ;{\Upsilon ^{\left( {old} \right)}}}}[\log p(\psi \mid \Upsilon ;\chi )] + {E_{r\mid \psi ;{\Upsilon ^{\left( {old} \right)}}}}[\log p(r;\varepsilon )]\\ \;\;\;\;\;\;\;\; = Q(\chi ) + Q(\varepsilon ) \end{array} $

(7)

观测上式可知，式中第一项仅包含超参数χ，第二项仅包含ε。因此，对上述两项分别求导可得到2个超参数[ε, χ]对应的学习规则（learning rule）。对于ε，其学习规则可写为

$ {\varepsilon _i} \leftarrow \sum\nolimits_{ii} + \mu _i^2 $

(8)

式中：$ \sum\nolimits_{ii} $为方差∑的第i个对角元素；μ_i为μ的第i个元素。

相应地，χ的学习规则可近似表征为

$ \chi \longleftarrow \frac{\|\psi-D r\|_{2}^{2}+\chi^{(o l d)} \sum\limits_{i=1}^{L}\left(\frac{1-\sum\nolimits_{ii}}{\varepsilon_{i}^{(o l d)}}\right)}{H} $

(9)

式中χ^(old)表示χ的上一次迭代结果。然后基于式（8）—（9），对ϒ = [ε, χ]开展迭代寻优确定最佳[ε_ML, χ_ML]，然后将其带入式（6）即可获取地层反射系数的最佳估值。

1.2 地层稀疏反射系数计算

研究区位于中国松辽盆地北部，区内发育有大型背斜，其构造特征从深到浅具有较好的继承性。前期研究结果表明^[34-35]，研究区发育多个含油层系，本次研究目的层所在层系为三角洲前缘亚相，具有单层厚度小、井间连续性差、平面上大面积错叠连片等特征；油藏类型相对复杂、多样，岩性与构造或断层共同控藏，具有较强的油气富集能力。

为满足研究区薄储层精细勘探，实施了高密度三维地震勘探，地震资料采集面元为10 m×10 m（图 1）。Xline和Inline方向分别有860道和840道地震数据，覆盖面积约290 km²，纵向时间为1 100~ 1 500 ms，采样间隔为1 ms；资料主频约55 Hz。采用商业软件提取了上述三维地震对应的地震子波（图 2），并结合研究区测井数据开展了合成记录对比（图 3）。由图 3可看出，基于提取的地震子波和测井曲线的合成地震数据与实际地震道具有较强的相关性，验证了子波提取的准确性。

在上述分析基础上，开展地震反射系数反演以提高薄层分辨能力。地震反射系数反演从数据作用域可分为时间域反演和频率域反演（谱反演）^[36]，本次基于稀疏贝叶斯学习的地震反射系数反演方法在时间域和频率域均可使用。以频率域反演为例，根据三维地震资料的频谱特征，设反演有效频段为5~90 Hz，频率采样间距为1 Hz，并在此基础上开展地震稀疏反射系数反演（图 4）。与原始地震数据体相比，地震反射系数反演抑制了地下薄层的邻层干涉效应，获得的地层反射系数更加清晰地展示了研究区的地质细节，提高了薄储层的分辨能力。

为了准确评价反演效果，并验证本次研究方法的优势，基于多种稀疏理论（基追踪、SBL-SA，SBL-EM）开展地震反射系数反演，然后结合提取的地震子波与不同反演结果开展正演模拟，最终通过对比模拟数据与实际资料的残余方差（Residual Variance，RV）来评价本次研究方法的准确性。由图 5可看出，与传统方法相比，当前方法计算的地层反射系数精度相对较高，为高精度地层相对波阻抗信息地预测奠定了基础。

2 地层相对波阻抗计算

地层反射系数表征地下地层界面对地震波的反射强度，其难以同传统地震振幅等属性一样开展切片分析，因此须将其转换为地层相对波阻抗以开展薄储层的精细分析^[12]。对于任意地震道，其对应的反射系数r (t)和地层波阻抗Z (t)的关系如下

$ r\left(t_{i}\right)=\frac{Z\left(t_{i+1}\right)-Z\left(t_{i}\right)}{Z\left(t_{i+1}\right)+Z\left(t_{i}\right)} $

(10)

由此可知，任意地层波阻抗Z (t_n)(1 ≤ n ≤ L)可通过反射系数r (t)表示为

$ Z\left(t_{n}\right)=Z\left(t_{n-1}\right) \frac{1+r\left(t_{n-1}\right)}{1-r\left(t_{n-1}\right)}=Z\left(t_{1}\right) \prod\limits_{i=1}^{n} \frac{1+r\left(t_{i}\right)}{1-r\left(t_{i}\right)} $

(11)

设最上部地层波阻抗Z (t₁)为1，则Z (t_n)可表示第n层的地层相对波阻抗且满足

$ Z\left(t_{n}\right)=\prod\limits_{i=1}^{n} \frac{1+r\left(t_{i}\right)}{1-r\left(t_{i}\right)} $

(12)

该式建立了地层反射系数和相对波阻抗的转换关系。由上式可知，任意地层的相对波阻抗估值在一定范围内和其上部所有层间反射系数有关，如果其上方任意位置处的地层反射系数计算有误或精度较低，则会累积对下方地层的相对波阻抗估值产生影响，形成低频累积误差。根据上述公式直接计算的地层相对波阻抗体存在较大的低频累计误差（图 6），由于该误差的存在，图 6中地层相对波阻抗对地层叠置关系表征不明确，无法开展薄层精细描述。因此有必要对上述低频累积误差予以压制，以提高其薄层分辨能力。

针对地层相对波阻抗计算中的低频累积误差问题，设计了一种线性FIR滤波器开展地层相对波阻抗滤波处理，以抑制低频累积误差。FIR滤波器通常具有稳定性强、精度高等优点，在信号处理中具有广泛应用。为了提高滤波效率，采用了线性相位FIR滤波器，其比非线性相位FIR滤波器减少了近50%的乘法运算^[37]，大幅度提高了滤波运算效率。对于上述地层相对波阻抗体（图 6），经过多次测试设定滤波阶数为n= 600，频率系数F =[0，0.012，0.032，1]，对应的高通滤波频率门槛值为[6 Hz，16 Hz]，剔除低频累积误差以提高地层相对波阻抗体的薄层描述能力，滤波效果如图 7所示。

由图 7（a）可知，原始相对波阻抗剖面几乎无法表征研究区的地层叠置关系，由于低频累积误差的存在，图 7（b）中相对波阻抗曲线（红色）趋势基本受低频累积误差（黑色）控制，不能有效突出其所蕴含的地层信息。经过滤波处理后，滤除低频误差后的相对波阻抗曲线[图 7（c）]对地下地层的表征能力得到明显提升。基于设计的线性FIR滤波器，通过参数试验[图 7（d）]选择最佳滤波参数开展高通滤波，处理后的地层相对波阻抗剖面[图 7（e）]能够清晰地展现地下地层在垂向上和横向上的展布特征，滤除的低频累积误差如图 7（f）所示。图 1中地震数据体包含近73万道数据，完成相对波阻抗计算和高通滤波处理共计用时15 s，验证了算法的高效率，滤波处理后的相对波阻抗体如图 8（a）所示，图 8（b）为对应的低频累积误差。

对比原始地震数据和滤波后的地层相对波阻抗的频谱特征（图 9）可知，地层相对波阻抗频谱的有效频带较宽，这对薄储层勘探具有重要意义。在此基础上开展地震振幅和相对波阻抗的切片对比，结果如图 10所示，地层相对波阻抗切片对研究区内地质体沉积特征（河道）较地震振幅表征得更加清晰、准确。

3 薄储层识别

如图 11所示，目的层位于松辽盆地北部埋深约1 695 m处，储层厚度约3 m，对应地震时间1 350 ms左右。测井分析可知，薄储层位置处声波时差值相对邻层增加，代表储层速度相对降低，而同时密度值相对降低，因此该薄储层相对波阻抗较邻层表现为低阻特征[图 11（a）]。在上述分析基础上，基于图 11（b）中层位处提取了该位置上地震均方根振幅属性[图 12（a）]和对应的地层相对波阻抗[图 12（b）]切片。此外，本次研究方法亦基于商业软件和研究区12口井资料开展了地层绝对波阻抗预测，并在相同位置提取了波阻抗切片[图 12（c）]以便结果对比。

由于地震分辨率和反射波邻层干涉效应的影响，地震均方根振幅[图 12（a）]对层厚小于地震波长一定比例的薄层的刻画能力通常不足。相比之下，由于地层相对波阻抗在计算过程中降低了地震邻层干涉效应的影响，同时亦拓宽了资料有效频带，因此其对薄储层的刻画能力更强[图 12（b）]。图中地震均方根振幅与研究区12口井数据的匹配率为67%，地层相对波阻抗匹配率为83%，符合率提高了近1/4。图 12（b）-（c）对比了计算的地层相对波阻抗与联合测井资料所预测的地层波阻抗，由图 12（c）可以看出，油井基本位于低阻区域而水井或干井位于高阻区域。将图 12（c）中刻画的波阻抗边界平移至图 12（b）的相对波阻抗切片上，观测可知，该边界亦能基本区分地层相对波阻抗的高低波阻抗区域，验证了本次研究方法的有效性，为薄储层精细勘探提供了重要参考信息。

4 结论

（1）由于地震资料分辨率及薄层干涉效应等因素的影响，基于地震振幅信息预测地下薄储层通常多解性较强。通过开展地震反射系数反演，进而计算地层相对波阻抗，能够有效提高薄储层平面展布形态预测精度。

（2）与基于传统方法的地震反射系数反演方法相比，新的稀疏贝叶斯学习理论（SBL-EM）提高了地层反射系数的反演精度，为高精度地层相对波阻抗预测奠定了基础。

（3）地层相对波阻抗计算过程中会产生低频累积误差，对薄储层识别影响较大。通过设计的线性FIR滤波器开展地层相对波阻抗高通滤波，能够极大地抑制低频累积误差，从而提高薄储层识别精度。引入SBL-EM理论开展地层相对波阻抗直接反演可作为后续研究重点之一，以提高相对波阻抗预测的稳定性。