随着油气勘探开发的不断深入，小断层、薄储层和岩性圈闭的识别与解释成为地震资料精细解释的重点(Peyton et al., 1998).谱分解技术不仅可提高地震资料对薄储层的解释与预测能力，而且可以从地震数据中挖掘丰富的油气及地质信息，一经推出便引起业界广泛关注，并得到了快速发展(Feng Shen，2003；Burnett et al., 2003；Hall，2006；Castagna and Sun, 2006；撒利明等，2015)，已在地震资料偏移、地层或沉积相解释、油气检测等方面取得广泛应用(陈林, 2009, 韩海英, 2014, 杨继东，2016;王旭谦，2017).谱分解技术的基本原理是选取一组在时间域和频率域均有限的基函数，通过观察待分析信号在基函数上的投影，得到信号在时间域和频率域的联合分布(Partyka et al., 1999).早期的谱分解技术通常采用短时傅里叶变换(STFT)，由于采用固定时窗，变换后的时频分辨率也即固定，不利于低频和高频信号的检测，为了改善傅里叶分析方法对非线性、非平稳信号分析的局限性，许多学者进行了深入研究，发展了许多方法，如小波变换、S变换等.如在时频分析中用长度变化的时窗来变换信号，就可更好地分析局部频率特征，这就是小波变换(CWT)(Mallat，1999；Chakraborty and Okaya, 1995；Marfurt and Kirlin, 2001).S变换是以Morlet小波为基本小波的连续小波变换的延伸，S变换可使用非正交小波将信号从时间域分解到时间-频率域，分解的效果得到了进一步提高(Antoine et al., 2002).但由于S变换中基本小波固定，且无法根据地震信号本身选择基本小波，在此基础上产生了广义S变换，用带参数的调幅简谐波代替S变换中的基本小波，然后利用不同参数的基本小波作线性组合形成最终的基本小波，因此其分频效果良好(Stockwell et al., 1996).除上述方法外，近年来发展起来的WVD(Wigner-Ville分布)、MP(匹配追踪法)以及HHT(希尔伯特-黄变换)等方法.WVD方法不需要时窗，可同时获得较高的时间域分辨率和频率域分辨率，但WVD方法会产生交叉项，这是一种很强的震荡噪音(Li and Zheng, 2007).MP方法遵循自适应信号分解原则，所求得的信号模型或时频特性能自适应地最佳匹配于信号的内部结构，可最有效地表征待分析信号；MP方法是一种信号的稀疏表示方法，它通过寻找数据在一个过完备库上尽可能好的匹配，从而达到稀疏表示的目的(Almoghrabi and Lange, 1986；Çostain and oruh, 2004；Aharon et al., 2006；Mallat and Zhang, 1993；Rubinstein et al., 2010, 张繁昌, 2012).该方法首先将信号表示为一系列时频原子的线性组合，然后通过迭代的方法求出时频原子模型的参数，缺点是时频原子的频率不随时间而变化(Liu and Marfurt, 2005).匹配追踪不需要事先确定一组变换的正交基函数，能够灵活地根据局部的信号特征确定一组非正交的基函数(也称为原子).由于这些原子的集合(即非正交基函数库)是过完备的，所以也称为过完备库(Chen et al., 1998).匹配追踪可以理解为是一种在由原子组成的过完备库上的信号的稀疏表示方法.相比小波变换，匹配跟踪是一种更简单、直接和有效的信号时频分解方法(武国宁，2012).匹配追踪算法的出现使得利用过完备库来实现信号的稀疏表示在大多数情况下成为可能(韩海英，2014).匹配追踪算法的关键在于过完备库的选择和求取, 但由于匹配追踪方法是一种贪婪算法, 且过完备库一般非常庞大，寻找真正的最佳匹配在计算上是非常昂贵和不现实的，目前大多数方法都是采用一组固有的过完备库，其缺点是不能适应不同特征的地震信号, 而且计算效率和分解精度不高.因此，作为贪婪算法, 如何快速、准确地实现从过完备字典中选择最匹配的原子是匹配追踪类算法备受关注的问题.为此，本文在常规MP算法的基础上，通过对常规MP算法进一步分析和研究，提出了一种新的完备库的构建方法和原子选择算法，不仅加快了MP算法收敛速度，提高了信号稀疏表示的质量，而且算法的自适应性和分解精度也大大提高.

匹配追踪算法的关键在于选取过完备库G，选取过完备库G的方法很多，目前大多数方法都是采用一组固有的过完备库G，其缺点是不能适应有不同特征的地震信号，而且计算效率和分解精度不高，为了克服常规MP方法的缺点，本文在对地震信号特别是地震子波研究的基础上，基于词典学习算法，提出了一种改进的过完备库G₀的构建方法，并在此基础上，给出了原子中心位置、原子中心频率的计算方法及原子间距选择策略等，形成了一种新的匹配追踪方法MPSD.下面介绍其算法原理及实现流程.

匹配跟踪利用过完备库来表示原始信号，如果把过完备库表示为G={g_n(t)}，n=0, 1, 2, …, ∞，匹配跟踪可以表示为

(1)

这里g_n(t)为原子即非正交基函数，而{a_n}就是信号s(t)在每个原子上的分解向量.

(2)

理想的过完备库，{a_n}的值在大多数原子上都几乎为零，所以它是一种信号的稀疏表示，能够更有效地代表原始信号的特征(蔡盛盛等，2014；邓承志和曹汉强，2009；李海山等，2014；王燕霞和张弓，2012).

做时频分解，通常选取已知形式的小波信号作为原子来构造过完备库，于是原子g_n(t)可以表示为

(3)

不同于小波变换，这些原子不需要具有正交性，如果母小波则原子

(4)

令位于t_n处的原子的频谱为

(5)

考虑地震信号的时频特点，本文选用Ricker子波来构建过完备库，这样就可以得到当前的原子，令Ricker子波为

(6)

则其对应的频谱为

(7)

所以，对于利用Ricker子波得到的原子，有

(8)

令ω_n=ω₀/L_n为小波的中心频率，得

(9)

其频谱为

(10)

其中ω_n为小波的中心频率，t_n为小波的中心位置.可见，对于此类原子的频谱，其振幅在中心频率处最大，随着远离中心频率其振幅不断衰减，且相位由子波的中心位置决定.确定所有原子以后，原始信号可表示为

(11)

可见，原始信号的频谱可以看成是所有原子频谱的叠加.由于每个原子都是一种中心位置确定的小波，这样就得到了原始信号的时频分解(Gribonval，2001；Liu et al., 2004；Wang，2007).

大量研究表明如果能够获得过完备库G中原子的中心位置，就可以大大减少原子的数量，并可能获得表征实际地震信号的特殊过完备库G，当前原子的中心位置可以通过中心频率来估算，中心频率利用最大相关性从原始地震数据中估计得到，具体可通过计算当前信号的解析信号，得到瞬时能量最高的位置，作为当前原子的中心位置.对于信号c(t)，当前原子的中心位置可定义为

(12)

其中H[·]表示Hilbert变换.

将当前信号在t_n处的瞬时频率ω_n定义为当前原子的中心频率：

(13)

仅仅使用原始地震数据来估计过完备库中原子的中心位置、中心频率，并不能完全解决常规MP方法分解质量和计算效率等关键问题，常规分解算法并不考虑原子间距对分解结果的影响.研究表明，原子间距是影响计算效率和分解质量的关键，据此本文对匹配跟踪算法还进行了以下改进，为避免原子间距过小的问题，引入了最小原子间距的概念，从理论上界定了原子间可允许的最小间距，这样，在原子产生的过程中，每个新的原子位置都要与以前的原子位置进行比较，当所有间距都大于或等于原子最小间距时，新的原子才可以使用.最小原子间距可根据地震资料情况由用户调节，测试中发现取原子最小间距等于10比较理想.

为了减少地震数据噪声、信噪比等对分解结果的影响，除了对地震数据进行精细的保证处理之外，本文方法在分解之前首先对输入的原始地震数据进行了平滑处理，分解后再一次对分解结果进行平滑等处理，这样既能够减轻噪声在产生解析信号时带来的影响，也可以得到较高分解质量的分频结果.

经典的匹配追踪算法的一般流程可表示为

(1) 将输入信号s(t)赋值给当前信号c(t):

(2) 找出过完备库G中和当前信号c(t)最相关(内积最大)的原子g_n(t)：

(3) 从当前信号减去找到的原子对应的分量：

(4) 如果‖c(t)‖ < 设定阈值，程序结束，否则回到第(2) 步继续.

基于2.1节描述，经过1) 选用Ricker子波构建原子；2) 利用原始地震数据最大相关性估计过完备库G的原子的中心位置和能量；3) 引进最小原子间距和原始地震数据平滑等改进后，本文MPSD算法流程可表示为

(1) 将输入信号s(t)赋值给当前信号c(t):

(2) 计算当前信号c(t)的解析信号：

其中H[·]表示Hilbert变换；

(3) 原子的中心位置由以下方法估计：

(4) 其中心频率ω_n=ω₀/L_n的估计定义为当前信号c(t)在t_n处的瞬时频率；

(5) 构建当前原子g_n(t)：

并将此原子加入过完备库G；

(6) 计算此原子对应的分量：

(7) 从当前信号减去找到的原子所对应的分量：

(8) 如果‖c(t)‖ < 设定阈值，说明找到了足够的原子；否则回到第(2) 步继续；

(9) 利用所有找到的原子分量，计算输入信号的时频分解：

其中

相对于其他分解算法，匹配追踪算法极其耗时，如何提高计算效率，分解效率和数据量存在密切关系，为了进一步提高计算效率，本文在算法优化的基础上，针对匹配追踪算法的特点，采用CPU多线程策略，极大地提高了计算效率.以Xeon系列E5-2670处理器(主频2.30 Hz)为例.图 1a为单CPU情况下本文算法与短时窗FFT算法计算效率对比，二者同时利用单线程对不同的样点数进行谱分解.可见在CPU主频一定的情况下，随着参与计算样点数增加FFT与本文算法的耗时都在快速增加，但本文算法的计算耗时的速度增大更快，在较多样点数时更加耗时，这是由于本文算法计算时使用大量冗余算子所影响.在CPU主频一定的情况下(图 1b)固定参与计算的样点数，通过CPU多线程进行并行改造，对两种算法进行耗时对比，可以看到随着参与运算的线程的增多，两种算法的计算效率得到大幅度提升，但随着线程的增多耗时性能提升速率相对降低.本文算法与FFT算法都具有较好的并行化特征，通过多线程并行加速能够达5~6倍的加速效果，且运算的样点数越多，并行加速效果越明显.

图 1

Fig. 1

图 1 基于CPU的多线程优化策略 (a)本文算法与FFT方法计算效率对比；(b)多线程加速对比. Fig. 1 Optimization strategy based on multi-thread of CPU (a) Efficiency comparison between MP and FFT; (b) Multi-threaded acceleration comparison.

由于该算法中过完备库是利用输入地震信号得到的，且引入了Ricker子波构建原子、最小原子间距以及地震数据平滑、CPU多线程等优化策略.因此，此算法不仅加快了收敛速度，而且提高了地震信号稀疏表示的质量，分解精度得到了极大提高.

为验证该方法的有效性，本文进行了相应的数值模型试算.设计主频为30 Hz、延续长度为200 ms、子波峰值在70 ms处的Ricker子波，然后分别利用小波变换CWT、S变换、常规匹配追踪算法、改进后的匹配追踪算法MPSD进行时频分析，主要测试分解算法的时-频局部性是否良好，是否具有较高的时间和频率域分辨率.图 2a为主频30 Hz的Ricker子波，图 2b—2e分别为小波变换CWT、S变换ST、原始匹配追踪算法、改进后的匹配追踪算法进行时频分解结果.可以看出，对于分频后的图像，两种匹配追踪算法的结果在时间和频率域的分辨率比小波变换(CWT)和S变换(ST)都要高，时频局部性好.由于此输入只包含一个Ricker子波，稀疏表示的原子数目很少，而且间距都比较大，所以改进后的匹配追踪算法(MPSD)和原始的匹配追踪算法结果类似.测试结果表明，改进后的MPSD算法具有良好的时-频局部性.通过模型试算，不仅验证了此算法可提高信号稀疏表示的质量，加快算法的收敛速度，同时也可提高匹配跟踪算法的自适应性，进而提高分解精度.

图 2

Fig. 2

图 2 主频为30 Hz的时频分解 (a) Ricker子波; (b) CWT; (c) ST; (d)原始MPSD; (e)改进的MPSD. Fig. 2 Frequency decomposition of 30 Hz (a) Ricker Wavelet; (b) CWT; (c) ST; (d) Original MPSD; (e) Improved MPSD.

实际地震数据经过时频分解后，结合三维可视化等技术可实现对储层、油气等特殊地质现象的精细解释.以下实例来自于某碳酸盐岩储层(采样率1 ms，主频约35 Hz，频带宽度约7~85 Hz)，图 3中为某一道地震数据的时频分解结果，其中图 3a为原始地震道，图 3b为CWT分解结果，图 3c为常规MP分解结果，图 3d为本文方法分解结果.图中显示出不同的反射率，在时频域可以精确地识别单个反射目标.对比图 3b—3d可以得出，本文方法具有较高的分解精度，反射能量团更加集中，时频分辨率更高.

图 3

Fig. 3

图 3 一个地震道的分频结果对比 (a)原始地震道；(b) CWT分解结果；(c)常规匹配追踪结果；(d)改进的匹配追踪结果. Fig. 3 Comparison of different frequency decomposing results (a) Original seismic signal; (b) CWT; (c) Original MPSD; (d) Improved MPSD.

图 4a—4c是用该方法产生的分频剖面，如在图中1.85 s、2.25 s处反射振幅从10 Hz到50 Hz逐渐增强并随后减弱，30 Hz左右最强.这种现象的出现要归结到储层含气特征，由于储层含气，使共振频率朝较高的一侧迁移，当照亮这些共振频率时，就很容易识别.因此利用不同频率的分解结果可以帮助实现在频率域分辨单个目标反射.图 5为10 Hz、30 Hz、50 Hz沿层切片RGB三原色多谱图象合成结果.左边为本文方法分解，右图为短时FFT方法分解，中等频谱区调谐的特征用绿色表示，低反射区表示为深色，其他用黄色或橙色表示.合成后图中储层得到进一步精细刻画，可以很好地对其局部特征进行精细分析，图中(图 5a)右上区域内地层的非均质性刻画更加清晰，可以精细识别多个溶洞，另外可以看到北西南东方向的早期断裂被后期发育的北东向断裂所切割，断裂系统更加清晰，特征明显，更加有利于精细解释.

图 4

Fig. 4

图 4 本文方法分解不同频率结果对比 (a) 10 Hz；(b) 30 Hz; (c) 50 Hz. Fig. 4 Comparison of different frequency component by improved MPSD

图 5

Fig. 5

图 5 不同频率分量(10 Hz红色、30 Hz绿色、50 Hz蓝色)多谱图象合成结果对比 (a)本文方法; (b)短时窗FFT方法. Fig. 5 Comparison of Fusion results from different frequency decomposed methods (a) Improved MPSD; (b) FFT.

利用改进的MPSD技术，可以对一些特殊的地质目标，如河道、岩性异常体、低频阴影等进行识别，并可进一步利用分解结果，开展基于图像学的多谱图象合成.离散型的频谱分解数据也可以通过动画形式帮助解释人员分析潜在的地质目标，并对这些潜在目标进行精细解释，极大地挖掘有效信息.这可以在部分程度上解决利用单一谱分量不能展示全部信息的缺陷，提高解释精度.

本文在对比分析常规MP方法的基础上，提出了一种新的分解算法MPSD，此算法不仅可加快收敛速度，提高信号稀疏表示的质量，而且分解精度和算法的适应性也大大提高.具体有如下几点认识：1) 提出了一种基于词典学习算法的完备库的构建方法，不仅加快了计算速度，而且提高了信号稀疏表示的质量；2) 利用Ricker子波构建原子，利用原始地震数据最大相关性估计过完备库G的原子的位置和能量，进一步提高了该算法的适应性；3) 为避免原子间距过小问题，引进了最小原子间距的概念，使分解质量和计算效率进一步提高；4) 模型试算和实际资料应用表明此本文方法更适合于油气检测、地震衰减研究和地震信号重构以拓展频带等，具有较好的应用前景.

本文中得到西北分院地球物理所GeoFrac裂缝预测软件研发项目组、GeoSeisQC地震采集质量监控软件项目组的大力支持，对他们的工作表示衷心感谢！