0 引言

由于地下地质构造的复杂性和岩性的多样性，地震子波会随空间和时间变化。储层中流体对子波的影响主要表现为地震记录的振幅、频率、相位的变化。常规的地震信号处理中的褶积、反褶积和地震道反演等都是基于单一地震子波的假设，与实际情况有很大的距离，导致现有的储层预测技术难以满足油气藏勘探的精度要求。匹配追踪地震道分解技术就是将以往只能从宏观上认识的地震道分解为可认识和可控制的不同形状、不同频率的地震子波，分解后不同频段子波数据体对储层的响应特征不同，通过合理的筛选可组合出精细表征储层物性和流体性质的全新的数据体，将其作为波场分析与储层预测的桥梁，可实现由地震反射信息向储层岩性及流体性质信息的转变，从而形成全新的储层及油气预测方法。

目前，将数据向量建模为基向量的线性组合，已大量应用于机器学习、神经科学、信号处理及统计学等领域。匹配追踪则是这类算法的典型代表之一。匹配追踪算法的基本思想就是可以把任何信号分解成一系列匹配子波的线性组合，这些匹配子波来源于一个过完备的匹配子波库。匹配子波库由同一基子波经过伸缩、平移和调制生成。Mallat等^[1]基于Gabor匹配子波库，提出匹配追踪算法，实现了信号自适应稀疏分解。然而，Gabor匹配子波波形与地震波形特征相差较大，不利于地震信号的分解。为此，Liu等^[2-3]先后以Ricker子波和Morlet小波构建超完备匹配子波库，取得了一定的应用效果。李海山等^[4]对匹配追踪算法中的过完备字典进行优化，通过匹配煤层强反射信息消除了其屏蔽作用，但是仍以理论子波作为字典库，分离信号中存在残余分量。张繁昌等^[5]针对这一问题提出利用井标定波形内插获得各道强反射匹配子波库，然而当煤层强反射横向变化较大时，该方法也会遇到类似的问题。何峰等^[6]从井旁地震道提取匹配子波，利用改进的匹配追踪算法去除煤层强反射、凸显目标砂体，但是当煤层强反射横向变化较大时，单纯从井旁道提取的子波作为子波库无法实现地震信号的线性精确表征。杨子鹏等^[7]提出多道联合约束的匹配追踪强反射压制方法，尽管可以考虑匹配过程中强反射的能量变化，但是匹配所用的子波库还是传统的理论子波，对应用效果造成了一定的影响。张在金等^[8]提出层位与子波约束下的匹配追踪去强轴技术，由于采用理论子波作为匹配子波库，剥离强轴后在储层下方发现明显的伴影异常。许璐等^[9]采用基于局部频率约束的动态匹配追踪开展强反射的识别与分离，取得了一定的效果。

由前人的研究成果可知，以理论子波作为匹配子波库不适用于复杂地震信号，虽然一些文献中提出了一些改进方法，但是都有一定的局限性^[10-16]。为此，本文提出了基于字典学习的子波库构建算法，从地震数据中学习获得字典的有限匹配子波而不是预先定义的字典，极大提高了信号分解的匹配度。针对目标区煤层反射不稳定的问题，选择按单个同相轴开展字典学习，从而构建了振幅、相位、主频和极性变化的单峰字典库。从实际应用效果看，基于叠前道集开展煤层强反射衰减处理效果更好，在此基础上进行叠前偏移，不但可以消除煤层强反射的影响，还能使得煤层附近有效信号的成像效果显著改善。在消除煤层强反射影响后，基于Winger-Ville(WV)分布的高分辨率时频属性，可以精确地刻画储层的变化特征^[17-18]。实际资料处理结果表明，与常规匹配追踪强反射压制方法相比，本文方法适用性更强，对于横向变化剧烈的强反射，效果更好。

1 K-SVD字典学习方法 1.1 字典学习的基本原理

在基于稀疏表示的方法中，字典的优劣对稀疏编码的稀疏度有着直接的影响，对信号选择性重构也起着非常关键的作用。一般情况下，超完备字典增加了匹配子波数量，匹配子波个数要大于信号维数，所以在超完备字典的空间域中，信号的表示并不是唯一的，但优点是可以通过学习和训练，自适应地搜索字典中匹配子波的最佳组合，从而更稀疏表示目标信号。字典学习是目前最常用的一种超完备字典构造方法，它能针对不同类型的信号，用训练数据自适应地学习与信号特性相匹配的冗余字典，再从训练好的超完备字典中寻找最能反映信号特征的若干个匹配子波，通过线性组合对信号进行最优稀疏近似。

字典学习的主要思想是利用包含K个匹配子波的字典矩阵D∈R^m×K稀疏线性表示原始样本Y∈R^m×N(其中N为样本个数，m是单个样本或匹配子波的维度)，有

$ \mathit{\boldsymbol{Y}} \approx \mathit{\boldsymbol{D X}} $

(1)

式中X∈R^K×N，为稀疏矩阵。需满足约束条件：X尽可能稀疏；D的每一列是一个归一化向量。可以将上述问题描述为优化问题

$ \min\limits_{\boldsymbol{D}, \boldsymbol{X}}\left[\|\boldsymbol{Y}-\boldsymbol{D} \boldsymbol{X}\|_{\mathrm{F}}^2+\lambda \varphi(\boldsymbol{X})\right] $

(2)

式中：||·||_F²表示Frobenius范数，定义为矩阵各项元素的绝对值平方的总和；λ为稀疏度控制参数，通常0 < λ < 1；φ(X)表示对X进行稀疏性约束，可用零阶、一阶或其他阶范数。

对于有两个未知的优化变量D、X的双准则问题，需要在稀疏表示误差和表示系数的稀疏性之间进行权衡，可以通过依次固定字典D和稀疏系数X，交替迭代求解，以降低求解难度。具体可分为以下两个主要步骤。

(1) 求解稀疏系数。假定字典D已知，则上述问题变为稀疏分解问题，可用稀疏分解算法，如LASSO(Least Absolute Shrinkage and Selection Operator)、OMP(Orthogonal Matching Pursuit)等方法求解优化问题

$ \min\limits_{\boldsymbol{X}}\left[\|\boldsymbol{Y}-\boldsymbol{D} \boldsymbol{X}\|_{\mathrm{F}}^2+\lambda \varphi(\boldsymbol{X})\right] $

(3)

(2) 字典训练。假定稀疏表示系数X已知，求解优化问题

$ \min\limits_{\boldsymbol{D}}\|\boldsymbol{Y}-\boldsymbol{D} \boldsymbol{X}\|_{\mathrm{F}}^2 $

(4)

训练出自适应的字典D。

两个步骤不断交替更新迭代，循环多次后，就可以得到与信号Y特征相匹配的超完备字典D以及对应的稀疏系数X。

1.2 K-SVD字典学习算法

K-SVD是一种关于稀疏表示的字典学习算法，该算法采用迭代交替学习方式，使用SVD(Singular Value Decomposition)优化并更新字典以更好地拟合数据。算法步骤如下。

(1) 初始化。从原始样本Y中随机取K个列向量{d₁, d₂, …, d_K}作为初始字典的匹配子波，得到初始字典D⁽⁰⁾∈R^m×K。令J=1，其中J为迭代次数。

(2) 稀疏编码。基于上一步得到的字典D^(J－1)，使用复数域快速匹配追踪得到每一个信号样本对应的稀疏系数向量的近似逼近解，进行稀疏编码，得到X^(J)。

(3) 逐列更新字典D^(J－1)：

1) 当更新d_k(k=1, 2, …, K)时，计算误差矩阵

$ \boldsymbol{E}_k=\boldsymbol{Y}-\sum\limits_{j \neq k}^K \boldsymbol{d}_j \boldsymbol{x}_j^{\mathrm{T}} $

(5)

式中x_j^T为X^T的第j个行向量；

2) 求取稀疏矩阵X^T第k个行向量x_k^T元素不为0的索引的集合ω_k={i|1≤i≤N, X_{k, i}≠0}；

3) 从E_k取出ω_k对应的列，组成新矩阵E′_k；

4) 对E′_k做奇异值分解E′_k=UΣV^T，其中U、V是正交矩阵，Σ是对角矩阵；

5) 取U的第1列更新字典中的d_k，令(x′_k)^T=Σ_{1, 1}V₁^T，其中V₁为V的第一个列向量。用(x′_k)^T更新X中的x_k^T；

6) 重复步骤1)~5)，更新所有子波。

(4) J+1⇒J，重复步骤(2)、步骤(3)，直到达到指定的迭代步数或收敛到指定的误差。

2 复数域快速匹配追踪

匹配追踪算法在每次匹配时通常要扫描全部过完备匹配子波库，其计算量较大、计算速度慢。因此，使用复数域快速匹配追踪(Complex Domain Fast Matching Pursuit Decomposition，CFMP)法提高计算效率。CFMP方法需要根据实地震记录计算复地震记录，从复信号中获取先验信息(如振幅包络最大值处的时间、瞬时频率以及瞬时相位等)，然后在先验信息的约束下扫描部分匹配子波得到最佳匹配，完成复信号的重构后再将结果返回实数域，从而在一定程度减小匹配中的计算量，提高算法效率。

实信号s(t)的复信号可表示为

$ S(t)=s(t)+\mathrm{i} s^*(t) $

(6)

式中s^*(t)为s(t)的Hilbert变换，其实质是90°相移。复信号的瞬时振幅、瞬时频率和瞬时相位分别为

$ A(t)=\sqrt{s^2(t)+\left[s^*(t)\right]^2} $

(7)

$ f_l=\left.\frac{1}{2 \pi} \frac{\mathrm{d}}{\mathrm{d} t} \tan ^{-1} \frac{s^*(t)}{s(t)}\right|_{t=u_l} $

(8)

$ \theta_l=\left.\tan ^{-1} \frac{s^*(t)}{s(t)}\right|_{t=u_l} $

(9)

式中：u_l是第l个波包络中心对应的时间；f_l和θ_l是相应的瞬时频率和瞬时相位。基于式(7)~式(9)，可选取K-SVD字典学习获得的子波峰值点的瞬时振幅、瞬时频率和瞬时相位作为这个子波的参数γ，构成学习子波库的一个匹配子波，这与传统的匹配追踪方法基于理论子波形成的子波库类似。基于信号的瞬时信息可以快速在字典库中寻找到相应的匹配子波。

假设K-SVD字典学习的子波库D={d_γ(t)| ||d_γ(t)=1}是Hilbert空间H中的向量集，每一个子波d_γ(t)都经过标准化，则复匹配子波可表示为

$ \boldsymbol{g}_\mathit{\boldsymbol{\gamma}}(t)=\boldsymbol{d}_\mathit{\boldsymbol{\gamma}}(t)+\mathrm{i} \boldsymbol{d}_\gamma^*(t) $

(10)

式中d_γ^*(t)为d_γ(t)的Hilbert变换。所有的g_γ(t)构成了复匹配子波库G。

再假设复信号F∈H，CFMP方法的目的是将F表示成从G中选择的复匹配子波的线性组合。

匹配追踪从初始残量R⁽⁰⁾=F开始。首先求取残量的瞬时振幅、频率和相位，得到残量的参数γ′。从G中选取参数γ与γ′相近的复匹配子波g_γ(t)计算新的残量

$ \boldsymbol{R}^{(n+1)}=\boldsymbol{R}^{(n)}-\left\langle\boldsymbol{R}^{(n)}, \boldsymbol{g}_\mathit{\boldsymbol{\gamma}}^{(n)}\right\rangle \boldsymbol{g}_\mathit{\boldsymbol{\gamma}}^{(n)} $

(11)

式中：n为匹配次数；R⁽ⁿ⁾为第n次匹配后的复残量；〈R⁽ⁿ⁾, g_γ⁽ⁿ⁾〉表示残差与子波的内积。其中g_γ⁽ⁿ⁾正交于R⁽ⁿ⁺¹⁾，因此有

$ \left\|\boldsymbol{R}^{(n+1)}\right\|_2^2=\left\|\boldsymbol{R}^{(n)}\right\|_2^2-\left|\left\langle\boldsymbol{R}^{(n)}, \boldsymbol{g}_\mathit{\boldsymbol{\gamma}}^{(n)}\right\rangle\right|^2 $

(12)

为了使||R⁽ⁿ⁺¹⁾||₂最小，则需要所选择的复子波g_γ⁽ⁿ⁾∈G使〈R⁽ⁿ⁾, g_γ⁽ⁿ⁾〉最大，即最佳匹配子波的选择条件为

$ \boldsymbol{g}_\mathit{\boldsymbol{\gamma}}^{(n)}=\operatorname{argsup}\left|\left\langle\boldsymbol{R}^{(n)}, \boldsymbol{g}_\mathit{\boldsymbol{\gamma}}^{(n)}\right\rangle\right| $

(13)

式中argsup表示上确界。

通过式(13)选择的匹配子波能使每次迭代计算的残量最小。

3 基于匹配子波的WV高分辨率时频分析

由于地下地质构造的复杂性和岩石的多样性，地震子波在传播过程中是随空间和时间变化的，储层中流体对子波的影响表现为地震记录的振幅、频率、相位的变化。相关研究表明WV分布具有极高的时频分辨率，它的时间分辨率和频率分辨率之积达到了测不准原理的极限分辨率^{[1, 19]}。由卷积定理可知，当信号为多分量信号时，其WV分布必然会出现交叉干扰项，出现很多“虚假信号”，这是WV分布的致命缺点。基于匹配追踪子波分解的高分辨率时频谱是地震数据体时频转换的一种新手段，将子波波形分解和高分辨率的WV时频分析方法结合起来，在提高时频分辨率的基础上，可有效消除交叉项的影响^[1]。子波波形分解后，地震信号已经变成了一系列具有不同振幅、频率、相位的地震子波，每个子波代表地震信号的局部特征，通过WV时频分析技术，针对单个子波进行时频分析，最后形成高分辨率的时频谱。

单个最佳匹配子波的WV分布可以表示为WVD[g_γ⁽ⁿ⁾]，多个分解匹配子波的WV分布可重构地震记录的WV分布

$ \operatorname{TFR}(t, f)=\sum\limits_{n=1}^M a_n \sqrt{\mathrm{WVD}\left[\boldsymbol{g}_\mathit{\boldsymbol{\gamma}}^{(n)}(t)\right]} $

(14)

式中：a_n为最佳匹配子波的振幅系数；M表示参与计算时频谱的匹配子波个数。

图 1a是由17个不同频率、不同振幅和不同时刻的子波叠加组成的信号。图 1b是该信号的传统WV时频分布，存在严重交叉干扰项，如图中红色箭头所指的两组信号能量之间多出来的能量分布。

图 1c是该信号基于匹配追踪算法的WV分布，可见17个匹配子波的WV分布的能量非常聚焦，时频分辨率高，没有交叉项的干扰。

4 实际资料应用 4.1 字典学习匹配子波库的构建

以中国西部A三维工区典型煤层强反射为例，说明字典学习子波库的构建过程。图 2是一条典型地震剖面(Inline2882)，煤层表现为强反射特征，呈两峰夹一谷。但该煤层的峰谷反射并不稳定，强煤层反射附近存在一些弱反射(白色箭头所示)，因此用完整的两峰夹一谷的强煤层反射特征构建时频匹配子波存在困难。

为此，采用单极性同相轴子波库构建的思路，利用字典学习，分别构建波峰和波谷的子波库，并分别提取峰、谷强反射特征，然后再进行强反射特征的加权衰减处理。图 3为基于字典学习构建的部分子波。全工区地震数据总共有3476856道，由于煤层反射在大部分区域相对稳定，因此学习生成的子波个数设定为3000，既保证子波库的冗余度，又能保证匹配追踪分解计算的速度。

4.2 叠后煤层强反射分离及加权衰减

通过字典学习，构建好基础子波库，基于层控的选择性重构可以按极性沿着煤层响应的峰谷特征分别提取相应的同相轴。匹配重构中，基于层控信息，按峰谷能量单独提取同相轴。图 4a是提取图 2所示的煤层响应的上部波峰部分，图 4b是提取的煤层响应的波谷部分，图 4c是提取的煤层响应的下部波峰部分，图 4d是合并后的煤层响应。由图 4d可见，如果按常规的带旁瓣的子波进行匹配追踪，在强反射波谷分叉的地方就很难处理。因此，对于横向变化较为复杂的煤层强反射，采用峰、谷分别提取的策略是可行的。

图 5是对图 2煤层强响应加权衰减70%能量后的地震剖面。与图 2对比可以看出，已有效消除了煤层对上、下弱反射的影响，弱反射能量和剩余的煤层响应的能量相当，为后续的储层精细描述和反演工作提供了高质量的数据基础。

为了进一步说明本文方法相对于传统匹配追踪方法消除煤层强反射的优势，特别选择一段反射较为稳定的煤层，对比二者的效果。图 6a是原始地震剖面；应用本文方法分别提取煤层强反射的上波峰、波谷和下波峰，再将三者叠加，结果如图 6b所示；图 6c为将本文方法提取的煤层强反射能量衰减70%后的重构地震剖面，可见，煤层强反射消除得自然、平顺，没有异常噪声引入。图 6d是基于传统匹配追踪方法、以Ricker子波作为基函数提取的煤层响应，存在一些噪声干扰。由于Ricker子波的波峰和旁瓣之间的能量差异较大，因此消除煤层强反射后，强波谷反射出现很多异常，下波峰的强反射消除也不理想，出现明显的信号残余(图 6e)。由此可见：本文提出的算法通过峰谷单独提取和处理，具有较好的自由度，效果更加理想；传统的匹配追踪分解方法，由于子波旁瓣和波峰之间相互影响，很难精确消除煤层强反射的影响，处理后易产生假象。

4.3 叠前煤层强反射分离及加权衰减

叠后煤层强反射衰减处理可显著消除强煤层反射的影响，但不能根本上改变煤层上、下弱反射信号的成像质量，因此有必要将煤层反射衰减处理用于叠前道集，以改善弱反射的成像质量。图 7a是原始动校正后的CMP道集，在道集上可见较强的煤层反射。图 7b是应用本文方法提取煤层强反射后、按能量加权衰减70%后的CMP动校道集，衰减后煤层的反射量与上、下其他反射波相当。

图 8是在叠前道集上提取煤层强反射，通过能量衰减处理，再进行叠前时间偏移获得的地震剖面，与图 5所示的叠后强煤层反射衰减相比，图中黄色、蓝色和紫色箭头所指的煤层上、下原来被强煤层反射掩盖的弱反射的连续性增强，信噪比更高，成像质量得到显著提高。

4.4 强反射衰减前、后WV时频分析

在煤层强反射衰减处理后，基于高分辨率的WV时频分析方法，说明叠前煤层强反射衰减处理对时频属性的影响。图 9为过Z3井的主测线Inline2566煤层强反射叠前衰减前、后WV单频(40Hz)能量属性对比。从图 9a可以看出，衰减前，时频属性主要反映煤层的强能量时频特征，其它反射细节模糊；从图 9b则可明显看出，衰减后，煤层上、下弱反射能量增强且更连续，储层细节清晰可见。从Z3井钻遇的目的层三工河组二段上亚段(J₁s₂¹)(黄色箭头所示)看，煤层强反射衰减后目的层时频属性剖面煤层上、下有效层段反射特征清晰，原来模糊不清的目标反射能量明显加强，在过井点的时频属性特征和伽马曲线的特征匹配度更高，为储层的精细描述奠定了良好的资料基础。

图 10是去煤层强反射叠前衰减前、后单频数据体沿J₁s₂¹层提取的均方根振幅属性，其中，Z3、Z110、Z102、Z106、Z5都是实际见产井，都处在储层相对发育区域，而强煤层衰减前显示Z110、Z102、Z5三口井都落在储层不发育区域，叠前煤层强反射衰减后目的层时频属性平面图与实钻结果吻合度更高，断层边界更清晰，可以更精确地刻画储层的空间展布。

5 结论与认识

(1) 基于字典学习的匹配子波库构建是可行的，由于波形特征来源于实际数据，更加符合实际资料的波形特征，可更精确地实现信号的匹配分解；

(2) 对于横向变化比较复杂的强煤层反射，用带旁瓣的子波进行匹配分解会出现假象，可采用单峰子波的匹配分解和重构，结果更精确、适用性更强；

(3) 基于快速匹配追踪的地震信号分解和重构，可有效提取强煤层反射信息，通过消除或衰减处理，可适度消除叠后数据中煤层强反射的影响；

(4) 在叠前道集上选择性提取煤层反射，通过适当的加权衰减重构，再进行叠前时间偏移，可有效的提高煤层附近目的层的成像质量；

(5) 在叠前煤层强反射衰减处理后，煤层上、下弱反射能量有所提高，再在子波分解基础上开展高分辨率的WV时频分析，则储层反射特征清晰且连续性较强，以较高的时频分辨率为煤层附近或下伏储层的精细刻画奠定了基础。