羌塘盆地位于青藏高原北部, 面积约为1.8×10⁵km²。盆地内有北羌塘坳陷、南羌塘坳陷和中央隆起带3个次级构造单元。发育了巨厚的中生代海相沉积地层, 三叠系、侏罗系地层广泛分布。近些年的研究表明, 上三叠统地层为羌塘盆地油气生成和天然气水合物成藏最好的目的层系^[1-5]。此次研究的主要目的层为上三叠统甲丕拉组(T₃j)、波里拉组(T₃b)和巴贡组(T₃bg)地层。在晚三叠世沉积时期, 研究区内整体以浅水三角洲-碳酸盐岩台地相沉积为主, 沉积岩性主要为碎屑岩和碳酸盐岩。其中, 研究区南部抱布德一带沉积环境主要为台地相和三角洲相, 沉积水体较浅, 水动力较强, 物源供应较丰富, 地层厚度大; 研究区北部沉积环境以水体相对较深的陆棚相为主, 物源供应相对减少, 地层厚度减薄, 总体上呈现出基本对称的地层结构^[6-8]。

雀莫错地区位于羌塘盆地的北羌塘坳陷带东部的雀莫错与各拉丹东山之间, 区内地形起伏大, 地形切割强烈, 地貌形态多样, 既有高山区的冰川侵蚀作用形成的地貌, 也有山麓地带发育的冰碛堤以及低山区的谷地、洼地或湖泊、河谷相间地貌。由于地表条件复杂, 研究区内采集得到的地震资料信噪比较低, 采用传统的处理流程得到的地震剖面, 其分辨率与横向连续性均较差。而反褶积作为一种反演反射系数序列的方法, 可应用于雀莫错地区地震资料的优化处理, 从而有效解决上述问题。

由于地震道的带限特性, 反褶积问题通常是病态的。由此得到的反射系数形式通常有无穷多组, 需要施加适当的约束从而获得合适的结果^[9]。在此前的研究中, 柯西准则^[10]和L₁范数^[11]已被广泛应用于稀疏反射系数反演。MALLAT等^[12]于1993年首先介绍了基于过完备字典的信号分解理论, 并提出了匹配追踪算法(MP)。匹配追踪算法将信号扩展为一个冗余字典中的小波(或原子)的线性叠加。通过寻找最符合条件的小波并计算它与原始地震道之间的残差, 反复迭代直至残差的能量低于预设的阈值, 从而得到反射系数序列的最优解^[13-16]。基于多通道的匹配追踪算法(MCMP)则进一步提高了空间的连续性, 并且获得了高分辨率的时频谱^[17]。在地震反射系数反演的方法中, 基追踪算法是一个相对较新的算法^[18-20]。在使用该方法进行反演的过程中, 由奇、偶偶极子组成字典用于分解反射系数对, 通过使用对数障碍法求得系数向量, 并将获得的系数乘以字典从而得到反射系数序列。将该方法用于求解欠定线性方程组时, 获得的解是唯一并且稀疏的, 这一特性说明基追踪算法适用于求解稀疏反射系数反演问题。然而, 以上介绍的传统方法是一个基于单道的过程, 即每一道的反射系数独立进行反演, 因此获得的反射系数剖面通常缺乏横向连续性。

针对雀莫错地区采集到的地震数据信噪比低, 传统的处理流程得到的地震剖面连续性较差的特点, 本文将F-X预测滤波^[21]与传统的基追踪算法结合起来, 对雀莫错地区的二维地震资料进行了有针对性的处理, 将预测算子与反褶积流程相结合, 最终获得了分辨率高、地下结构连续性较好的地震反射剖面。

1 基于F-X预测滤波的基追踪算法

假设地层介质为一系列具有固定速度和密度的水平层状介质, 地震记录道可表示为反射系数和子波的褶积:

$ d\left( t \right) = w\left( t \right)*r\left( t \right) + n\left( t \right) $

(1)

式中:d(t)为实际观测地震记录; w(t)为地震子波; r(t)为待求反射系数; n(t)为噪声。子波一般可通过高阶统计法求得^[22-23]。

为了利用基追踪算法求解反射系数, 可将反射系数剖面看作一系列地层反射系数对的叠加, 而反射系数对可进一步表示为偶反射对(r_e)和奇反射对(r_o)的加权叠加, 如图 1所示。

图 1因此, 反射系数序列可表示为:

$ \begin{array}{l} r\left( t \right) = \sum\limits_{n = 1}^{{n_{{\rm{max}}}}} {} \sum\limits_{l = 1}^{{l_{{\rm{max}}}}} {[{a_{n, l}}{r_{\rm{e}}}\left( {t, l, n} \right)} + \\ \;\;\;\;\;\;{\rm{ }}{b_{n, l}}{r_{\rm{o}}}\left( {t, l, n} \right)] \end{array} $

(2a)

$ {r_{\rm{e}}}\left( {t, l, n} \right) = \delta (t-l\Delta t) + \delta (t-l\Delta t + n\Delta t) $

(2b)

$ {r_{\rm{o}}}\left( {t, l, n} \right) = \delta (t-l\Delta t)-\delta (t-l\Delta t + n\Delta t) $

(2c)

式中:Δt为采样率; nΔt为两个反射层之间的时间厚度; a_{n, l}和b_{n, l}表示相应的反射系数对的权重系数; l用于调整该奇或偶反射对在地震道上的位置, 可上移或下移, 称为地震道时序。考虑到所有可能的地层厚度, n的范围可以从1变化到表示最大地层厚度的常数n_max。

进一步将(2a)式写为矩阵矢量相乘的形式, 可以得到:

$ \mathit{\boldsymbol{r}} = \mathit{\boldsymbol{G'm}} $

(3)

式中:r表示反射系数的向量形式; m是由a_{n, l}和b_{n, l}构成的系数的向量形式; G′是由奇、偶偶极子组成的楔形反射矩阵。因为楔形模型是用于分解反射系数的过完备字典, 系数向量m通常比r更加稀疏, 意味着m更适于在稀疏约束条件下求解。

根据(1)式所描述的褶积模型, (3)式等号两侧都可以与子波进行褶积, 从而得到地震道的矩阵矢量相乘形式:

$ \mathit{\boldsymbol{d}} = \mathit{\boldsymbol{Gm}} $

(4)

式中:d表示地震道的向量形式; G是子波与G′褶积得到的楔形地震响应矩阵。从公式(3)可知待求的未知量为a_{n, l}及b_{n, l}, 与公式(4)中的未知量相同, 因此可以利用基追踪算法基于输入的地震道信号, 已知的子波及楔形反射矩阵求得这些未知量, 并且代入公式(3), 从而获得最终的反射系数。

在反演过程中, 系数向量m的初始值由下式求得:

$ \mathit{\boldsymbol{m}} = {({\mathit{\boldsymbol{G}}^{\rm{T}}}\mathit{\boldsymbol{G}})^{-1}}{\mathit{\boldsymbol{G}}^{\rm{T}}}\mathit{\boldsymbol{d}} $

(5)

基于初始模型, m可以通过最小化目标函数((6)式)求得:

$ J = {\left\| {\mathit{\boldsymbol{d}}-\mathit{\boldsymbol{Gm}}} \right\|_2} + \lambda {\left\| \mathit{\boldsymbol{m}} \right\|_1} $

(6)

式中:λ为控制数据残差及L₁范数模型约束两者平衡的权衡参数。在基追踪算法中, (6)式可以转换为一个线性规划问题, 并可以利用原始对偶对数障碍法加以求解^[24-25]。得到矢量m后, 反射系数可以通过公式(3)求得。

一个含有平层和倾斜结构的合成地震数据(图 2a)被作为输入来验证传统基追踪方法的有效性。图 2b为采用传统基追踪算法得到的反射系数剖面, 可见, 其结果非常干净, 并具有连续的结构。

以上合成地震数据的实例表明了基追踪算法的有效性。然而, 当输入地震道的信噪比较低时, 反演得到的反射系数剖面横向连续性将会变差。图 3a为添加了随机噪声后的合成地震数据(信噪比为2.0), 图 3b为利用基追踪算法得到的反射系数结果。可以发现, 直接使用基追踪算法求解时, 可以去除绝大多数的噪声, 但一些有效的地震信号也会被破坏, 从而影响了剖面的横向连续性。为了进一步提高在低信噪比背景下得到的反射系数的连续性, 将F-X预测滤波与传统基追踪算法相结合用于低信噪比地震数据的反射系数求取。

F-X预测滤波方法基于以下假设, 即地震信号图像可以由频域中一系列线性相干的反射叠加组成。在此假设条件下, 地震信号是严格的正弦形式, 并在频域具有可预测性。因此, 对某一固定频率f, 某一道地震信号可以通过其相邻的地震道预测获得:

$ \begin{array}{l} {D_i}\left( f \right) = {P_1}D{'_{i-1}}\left( f \right) + {P_2}D{'_{i-2}}\left( f \right) + \ldots + \\ \;\;\;\;\;\;\;\;\;\;\;\;{\rm{ }}{P_L}D{'_{i-L}}(f) \end{array} $

(7)

式中:P_i表示长度为L的预测滤波器; D_i为频域中经过平滑后的某道地震信号; D′_i为平滑之前的该道地震信号。

由于输入的地震信号和反演得到的反射系数有相似的空间连续性, 根据其可预测性, 反射系数可表示为:

$ \begin{array}{l} {R_i}\left( f \right) = {P_1}R{\prime _{i-1}}\left( f \right) + {P_2}R{\prime _{i-2}}\left( f \right) + \ldots + \\ \;\;\;\;\;\;\;\;\;\;\;\;{P_L}R{\prime _{i-L}}(f) \end{array} $

(8)

式中:R_i和R′分别表示频域中平滑之后及之前的反射系数。

基于F-X预测滤波的基追踪算法具体步骤描述如下^[26-27]。

步骤1:使用传统的基追踪算法计算反射系数r′。

步骤2:对求得的r′应用F-X预测滤波得到滤波后的反射系数r, 从而有效去除反射系数剖面上的噪声, 增强剖面的横向连续性。需要指出的是, F-X预测滤波在频域中进行, 而基追踪算法在时间域中实现, 因此, 在实际操作中, 必须先利用傅里叶变换将反射信号转换到频域, 待预测滤波完成后再变换回时间域。

步骤3:根据(1)式将滤波后得到的反射系数r与已知的子波褶积, 得到合成地震记录d。

步骤4:利用(5)式重新计算系数向量m的初始值。

步骤5:重复基追踪算法流程, 根据步骤4中得到的最新初始值m, 迭代最小化目标函数直至获得满意的结果。

根据地震道的可预测性, 采用基于F-X预测滤波的基追踪算法可以有效提高最终反射系数剖面的横向连续性。为了验证算法的效果, 将添加噪声的合成地震记录(图 3a)作为输入地震信号来进行效果分析, 结果如图 4所示, 可以看出, 相比于传统的基追踪算法, 采用基于F-X预测滤波的基追踪算法能够获得更干净和更连续的反射系数剖面。

2 实际地震数据处理及效果分析

以雀莫错地区QMC13线为例进行实际地震数据处理及效果分析。如图 5所示, 由于在雀莫错地区采集时使用的是可控震源, 导致单炮记录上初至难以拾取, 静校正问题突出。同时, 由于全区不同位置上的单炮记录的信噪比及目的层品质均较低, 主要存在面波、声波、线性干扰、异常振幅和高低频噪声等(图 6), 原始叠加能量纵向衰减快, 横向变化大, 平面上同一测线炮间能量不均衡, 导致产状无法识别。

在雀莫错地区二维地震资料的处理流程中, 针对初至难以拾取的问题, 选取高程静校正方案解决静校正问题^[28-31], 同时测试了采用多种保真、保幅的叠前去噪方法压制面波及线性干扰等干扰波的效果。通过测试叠前反褶积参数, 有效提高了地震资料的分辨率, 拓宽了信号的频带。采用速度分析和剩余静校正迭代处理技术, 解决了速度建模以及短波长问题。在做好噪声压制和精细速度分析以提高成像精度的基础上, 采用多次速度调整偏移, 同时结合构造及断层的合理性分析, 最终建立精细偏移速度场, 确保偏移归位成像效果。图 7为利用上述处理流程得到的QMC13线偏移剖面, 可以看出, 剖面的分辨率及横向连续性不甚理想。

鉴于此, 本文采用基于F-X预测滤波的基追踪算法来实现地震反射系数反演, 以期提高地震资料的分辨率及连续性。图 8a和图 8b分别为使用传统的基追踪算法及优化后的算法得到的反射系数剖面。虽然直接观察整体剖面无法明确分辨反射系数剖面的改善效果, 但对比图 7和图 8可以看出, 在进行反射系数反演后的剖面上, 噪声被有效地去除, 整体构造形态和地层展布规律更清晰、更易于被刻画出来。

为了进一步展示F-X预测滤波基追踪算法的效果, 选取图 8a和图 8b剖面中时间0.15~0.45s、道号1400~1500的局部剖面进行放大展示(图 9)。对比红圈中的波组特征可以看出, 在使用改进算法后, 波组更加连续, 一些细小的层位也得到了恢复。

图 10为反射系数及原始地震数据的全剖面振幅谱, 其中黑色曲线表示原始输入信号的振幅谱, 红色曲线表示利用基于F-X预测滤波的基追踪算法求得的反射系数的振幅谱。从图 10可以看出, 反演后, 振幅谱的高频部分被很好地补偿, 频带得到了大幅度拓宽, 说明采用F-X预测滤波基追踪算法可以求解得到连续的反射系数剖面, 并且能有效地提高剖面的分辨率。

利用基于F-X预测滤波约束的基追踪算法得到的反射系数剖面对研究区目的层进行构造解释。研究区整体上呈周低中高的构造格局, 背斜主要由中央古隆起带由南向北的挤压作用形成。该区现今的构造形态主要以褶皱为主, 断层断距较小, 背斜形态较完整, 上三叠统地层沉积厚度较大且在区域上分布稳定, 为天然气水合物提供了有利的成藏条件。图 11和图 12分别为采用传统基追踪算法和基于F-X预测滤波的基追踪算法得到的反射系数剖面的构造解释结果。对比图 11和图 12可以看出，使用基于F-X预测滤波的基追踪算法处理后的地震剖面分辨率得到提高，波组特征明显、连续性增强，能够更加清晰地刻画出研究区构造形态及地层展布规律。

3 结论

传统基追踪算法在反演复杂构造区二维地震资料数据时, 存在地震剖面的分辨率和连续性欠佳的问题。本文采用基于F-X预测滤波的基追踪算法, 在原有基追踪算法的迭代中加入了F-X预测算子, 对每次迭代的结果进行预测滤波, 在增强结果的连续性后将其作为下一次迭代的初值, 从而考虑了不同地震道之间的相关性, 在提高分辨率的同时有效提高了反射系数剖面的连续性。

在实现该算法的过程中, 一方面利用基追踪算法控制反射系数的稀疏程度, 另一方面利用F-X预测滤波器控制反射系数的横向连续性。对合成记录进行测试, 发现本方法能够在保证反射系数剖面横向连续性及分辨率的基础上, 有效地压制噪声。

由于羌塘盆地特殊的地理位置及环境条件约束, 地震资料普遍存在品质较差(如低信噪比)等问题, 地震资料的精细处理成为西藏羌塘地区油气及天然气水合物勘探突破的关键点。雀莫错地区的实例应用结果表明, 使用该方法处理后的地震反射系数剖面分辨率得到提高, 波组特征明显、连续性增强, 有利于更加清晰地刻画出研究区构造形态及地层展布特征。说明采用本方法可以获得高分辨率及高连续性的二维地震资料, 在类似的复杂环境及构造地区二维地震资料处理过程中, 具有较好的使用效果与应用前景。