非平稳地震信号匹配追踪时频分析
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
2 基于匹配追踪的时频表示方法
2． 1 基于 Wigner-Ville 分布的匹配追踪时频表征 通过匹配追踪，地震信号 s（ t） 被分解成一系列
匹配子波 wγn（ n = 0，1，…，N － 1） 之和，通过对各匹 配子波的 Wigner-Ville 分布（ WVD） 求和，便得到原 信号的时频分布［2］
图 5 合成地震记录时频分布
时频分布
时相位等时频信息。式（ 6） 虽然实现了变窗口短时
∞
∑ TFR2 （ t，f） =
gn （ τ － t） wγn （ τ） e －i2πfτ dτ ，（ 6）
－∞
其中，gn（ t） 为与子波 wγn（ t） 对应的时窗。虽然从整
体上看短时傅氏变换的窗口是变化的，但对于某一
交叉项和窗口截断效应。通过与短时傅氏变换和 S 变换时频特征的对比发现，匹配追踪时频表征和瞬时谱参数具
有更高的分辨率。实际数据的应用也表明，匹配追踪分解非常适用于非平稳特征的地震信号的时频分析。
关键词：匹配追踪； 时频分析； 匹配子波； 非平稳信号
中图分类号： P631． 4
文献标识码： A
文章编号： 1000 － 8918（2011）04 － 0546 － 07
第 35 卷第 4 期 2011 年 8 月
物探与化探
GEOPHYSICAL ＆ GEOCHEMICAL EXPLORATION
Vol． 35，No． 4 Aug． ，2011
非平稳地震信号匹配追踪时频分析
张繁昌，李传辉
（ 中国石油大学 地球科学与技术学院，山东 青岛 266555）
摘 要： 根据三步法匹配追踪原理实现了基于雷克子波的地震信号自适应分解，在此基础上，讨论了利用匹配子波
进行地震信号时频表征的方法。由于常规匹配追踪时频是以 Wigner-Ville 分布为基础，得到的时频信息有限，因此
给出了一种基于可调窗口的短时傅氏变换时频表示方法，进而又利用各匹配子波的复谱，引入一种新的时频表示
方法，不仅与 Wigner-Ville 方法具有同等的分辨率和能量聚集特性，而且保留了原信号的最基本时频特征，不存在
由于 Mallat 提出的匹配追踪时频分析方法是用 各匹配子波的 Wigner-Ville 分布表示的，这种时频 表征方式只能提供原信号的时频能量（ 振幅） 分布， 而没有相位信息，而且每个匹配子波的 Wigner-Ville 分布均存在交叉项，在对所有匹配子波的 WignerVille 分布求和时，交叉项能量由于不能实现同相叠 加而受到压制。笔者根据匹配追踪的算法特点，提 出了一种基于可调窗口的 STFT 时频表示方法，通 过逐步递减的窗口长度克服了常规短时傅氏变换单 一分辨率的缺点，而且此方法本身就不产生交叉项，
a—利用式（ 10） 得到的振幅时频分布； b—利用式（ 10） 得到的相位时频分布
图 7 合成地震记录时频表征
3 匹配追踪时频表征与 S 变换的对比
相位雷克子波组成。 对图 8a 的合成地震记录分别进行 S 变换和用
波与原信号越相似，则分解效果越好，由于雷克子波
具有理想的地震子波波形，在对地震信号进行匹配
追踪分解时，可以利用雷克子波作为基本子波来创 建子波库［3］。
设地震信号为 s（ t） ，由式（ 1） 可知，匹配追踪分
解成雷克子波的组合为
∑ s（ t） = ai ·r（ t － ti ，fi ，φi ）
（ 4）
·549·
图 6 左侧曲线是 50 Hz 的雷克子波包络，上部 曲线是其频谱，二者根据式 （ 7 ） 得到其时 频 分 布。 从图 6 看出，匹配子波的包络和频谱分别控制了其 在时域和频域的分布范围，具有较高的能量聚焦和 时频分辨率。与式（ 5） 不同，式（ 7） 给出的是最基本 的信号时频表征，可以得到时频分布的实部、虚部、 振幅、相位等多种信息。而且与式（ 6） 相比，式（ 7）
图 1 不同主频雷克子波及其合成记录
图 2 地震记录与匹配追踪重构结果的比较
·548·
物探与化探
35 卷
图 3 图 1 地震记录的匹配追踪振幅时频分布
图 4 图 1 地震记录的 WVD 振幅时频分布
a—由式（ 6） 得到图 1 地震信号的时频分布； b—直接对图 1 地震信号进行 STFT 得到的时频分布
没有引入任何窗口，因而不存在窗口截断效应。 仍然以图 1 的地震信号来验证式（ 7） 给出的时
频分布。图 7a、b 分别给出了由式（ 7） 计算出的地 震信号振幅和相位的时频分布。对比图 7a 与图 3 可以发现，利用式（ 7） 得到的时频分布，无论是在时 频分辨率还是在能量聚集性上可以与式（ 5） 的计算 结果相媲美。
选择 wγn ，使 |〈Rn s，wγn 〉| 最大，则残差‖Rn + 1 s‖2 最
小。因此，方程（ 2） 中匹配子波的振幅满足如下条件
｛ an ，γn ｝
= arg max ‖Rn+1 s‖2 。 a n，γ n∈Γ
（ 3）
1． 2 基于雷克子波的匹配追踪
匹配追踪算 法 中，创 建 子 波 库 所 使 用 的 基 本 子
主频的匹配子波来说，其短时傅氏变换窗口是不变
的，因此该时频表征方法对原地震信号的正反变换
都不会产生畸变。这样，该方法将匹配追踪得到的
匹配子波与短时傅氏变换结合，形成了变窗口短时
傅氏变换时频表示，而且不存在交叉项干扰。
傅氏变换，但仍然难免窗口的截断效应。
为了避免窗口截断效应，并且得到与相位有关
的信息，根据匹配追踪的算法特点，通过对每个匹配
进而笔者引入一种无窗口截断效应的时频表示方 法，通过匹配子波沿时间轴的能量包络和沿频率轴 的振幅包络双重控制，不仅能够保持 Mallat 方法时 频表示的分辨率，而且不产生交叉项。
1 匹配追踪基本原理
匹配追踪是一种基于子波库扫描的信号自适应
分解方法。设 D 为进行信号分解的超完备子波库，
信号为 s（ t） 。D 中的元素满足 D = ｛ wγ ∶ γ ∈ Γ｝ ， wγ 为其中的一个匹配子波且‖wγ ‖ = 1。匹配追踪 算法每 次 迭 代 都 自 适 应 地 提 取 出 一 个 最 优 子 波
N －1
∑ TFR1 （ t，f） = an 槡WVD［wγn（ t，f） ］ （ 5） n =0
其中，WVD［wγn（ t，f） ］表示匹配子波 wγn （ t） 的 Wigner-Ville 分布。
图 3 为利用式（ 5） 对图 1 最左边的合成地震记 录计算得到的匹配追踪时频分布。从图中可以看 出，合成该地震记录的各子波分量在时频图中均有 体现，而且在 1． 5 s 处相隔较近的两个子波也可以 分开，几乎没有交叉项能量的干扰。图 4 是直接对 图 1 最左边的合成地震信号求 Wigner-Ville 分布的 结果，可见存在很强的交叉项能量，严重影响信号的 时频表示结果。 2． 2 基于短时傅氏变换的匹配追踪时频表征
的中心时间、主频和相位，称为子波的控制参数。每
次迭代中，最优子波 wγn 的控制参数的确定是整个 匹配追踪算法的核心。
1． 1 匹配子波控制参数的优选
对于一道地震信号，利用 Hilbert 变换将其构造
成一个复信号，将复地震信号包络最大值处的时间
作为 tn 的初值，最大值处对应的瞬时频率作为 fn 的 初值，最大值处的瞬时相位作为 φn 的初值。设 tn 的初值为 t0n，采样间隔为 Δtn，则其扫描范围为
收稿日期：2010 － 06 － 19 基金项目：国家 863（ 2006AA09A102） 资助项目
ቤተ መጻሕፍቲ ባይዱ期
张繁昌等： 非平稳地震信号匹配追踪时频分析
·547·
［t0n － m·Δtn ，t0n + m·Δtn ］， m 为正整数，用来控制扫描范围。用同样的方法可
以设置 fn 和 φn 的扫描范围。通过对每个控制参数 在其初始值附近进行扫描，可以确定出每个控制参
wγn，n 为迭代次数。经 N 次迭代之后，
N －1
∑ s（ t） = anwγn（ t） + RNs ，
（ 1）
n =0
其中，an 是第 n 个匹配子波 wγn 的振幅，是迭代 N
次之后的残差，且 R（ 0） s = s。匹配子波 wγn 可以由 γn
= ｛ tn，fn，φn｝ 来刻画，其中 tn、fn 和 φn 分别为子波
数的最优值。tn，fn，φn 三个控制参数一旦确定，匹
配子波 wγn也就随之确定，然后利用
an = | 〈Rn s，wγn 〉|
（ 2）
来确定其振 幅 其 中，〈，〉表 示 内 积 运 算。 在 迭 代 分
解过程中，由于 wγn与 Rn + 1 s 正交，因此有
‖Rn s‖ = | 〈Rn s，wγn 〉| 2 + ‖Rn s‖2 ，
i
其中，ai 、ti 、fi 、φi 分别表示第 i 个雷克子波的振幅、
中心时间、主频和相位。
为了验证基于雷克子波的匹配追踪分解效果，
合成了图 1 所示的地震记录（ 图 1 左边第 1 道） ，该
记录由不同主频和相位的雷克子波叠加在一起组
成，其中包括 2 个零相位的 10 Hz 子波（ 图 1 左边第
短时傅氏变换不产生交叉项，但因其窗口不可 调而使得其时频表示的分辨率受到限制。考虑到各 匹配子波的独立性，笔者提出如下的匹配追踪时频 表征方法： 对各匹配子波进行短时傅氏变换时，对低 频匹配子波采用较大时窗，而对高频匹配子波则采 用较小时窗，这样便可以克服常规短时傅氏变换窗 口固定、分辨率单一的缺点，得到具有较高分辨率的
子波的复谱求和，可以得到如下时频分布表达式
∞
∑ TFR3（ t，f） = anWγn（ f） env［wγn（ t － tn） ］ei2πf（ t－tn） ，（ 7） －∞
其中，Wγn（ f） 为匹配子波 wγn （ t） 的傅氏谱，env［wγn （ t － tn） ］为匹配子波的瞬时包络（ 图 6） 。
基于这种思路，对图 1 所示的合成地震记录采
用长度逐步递减的时窗（ 100 ～ 10 ms） ，分别对从低
频到高频对匹配子波依照式（ 6） 求其时频分布，得
到图 5a 所示的结果。与直接对图 1 的地震信号进
行 100 ms 固定窗口短时傅氏变换相比（ 图 5b） ，时
频分辨率明显改善。但将图 5a 与图 3 对比可见，本
方法的分辨率不如利用式（ 5） 的分辨率高。
2． 3 无窗口截断效应的匹配追踪时频表征
匹配追踪时频分布公式 （ 5 ） 可以保持 Wigner-
Ville 分布的高分辨率，并且压制交叉项的影响，但
它只能给出信号的振幅分布特征，无法从中得到瞬
图 6 50 Hz 雷克子波的时频分布
4期
张繁昌等： 非平稳地震信号匹配追踪时频分析
匹配追踪（ MP） 时频分析方法作为新一代谱分 解技术，在地震解释、储层识别及烃类检测等方面具 有很大的应用潜力［1］。匹配追踪方法由 Mallat 等人 首次提出，其原始算法是基于 Gabor 函数组成的超 完备子波库贪婪算法［2］，并以子波库中各匹配子波 的 Wigner-Ville 分 布 之 和 来 对 原 信 号 进 行 时 频 表 征。Liu J 和 Marfurt 在 2004 和 2005 年先后提出了 基于雷 克 子 波 库 和 Morlet 子 波 库 的 匹 配 追 踪 算 法［3 － 4］，并将地震信号的瞬时特征引入到匹配追踪 算法当中，提高了算法的分解效率； Wang Y 在 Liu J 算法的基础上，提出了三步法匹配追踪方法，其时频 表示方法仍然是 Wigner-Ville 分布［5］。笔者在三步 法原理的基础上，实现了基于雷克子波的匹配追踪 分解，并与短时傅氏变换（ STFT） ［6］和目前流行的 S 变换（ ST） ［7］进行对比，证明匹配追踪时频表征是一 种非常理想的地震信号时频分析工具。
2 道） ，4 个零相位的 30 Hz 子波（ 图 1 左边第 3 道） ，
2 个 90°相位的 50 Hz 子波（ 图 1 右边第 2 道） ，1 个
45°相位的 30 Hz 子波（ 图 1 右边第 1 道） ，各个子波
的振幅相同，记录长度为 1． 6 s，2 ms 采样。将此合
成地震记录进行基于雷克子波的匹配追踪分解，并 利用方程（ 4） 进行重构，其结果如图 2 所示，其中最 左边 1 道是图 1 最左边的原始地震记录，中间的 1 道是匹配追踪重构结果，最右边是前二者的差。通 过图 2 的对比可见，该地震道可以完全由匹配子波 重构，残差曲线接近零值。
2． 1 基于 Wigner-Ville 分布的匹配追踪时频表征 通过匹配追踪，地震信号 s（ t） 被分解成一系列
匹配子波 wγn（ n = 0，1，…，N － 1） 之和，通过对各匹 配子波的 Wigner-Ville 分布（ WVD） 求和，便得到原 信号的时频分布［2］
图 5 合成地震记录时频分布
时频分布
时相位等时频信息。式（ 6） 虽然实现了变窗口短时
∞
∑ TFR2 （ t，f） =
gn （ τ － t） wγn （ τ） e －i2πfτ dτ ，（ 6）
－∞
其中，gn（ t） 为与子波 wγn（ t） 对应的时窗。虽然从整
体上看短时傅氏变换的窗口是变化的，但对于某一
交叉项和窗口截断效应。通过与短时傅氏变换和 S 变换时频特征的对比发现，匹配追踪时频表征和瞬时谱参数具
有更高的分辨率。实际数据的应用也表明，匹配追踪分解非常适用于非平稳特征的地震信号的时频分析。
关键词：匹配追踪； 时频分析； 匹配子波； 非平稳信号
中图分类号： P631． 4
文献标识码： A
文章编号： 1000 － 8918（2011）04 － 0546 － 07
第 35 卷第 4 期 2011 年 8 月
物探与化探
GEOPHYSICAL ＆ GEOCHEMICAL EXPLORATION
Vol． 35，No． 4 Aug． ，2011
非平稳地震信号匹配追踪时频分析
张繁昌，李传辉
（ 中国石油大学 地球科学与技术学院，山东 青岛 266555）
摘 要： 根据三步法匹配追踪原理实现了基于雷克子波的地震信号自适应分解，在此基础上，讨论了利用匹配子波
进行地震信号时频表征的方法。由于常规匹配追踪时频是以 Wigner-Ville 分布为基础，得到的时频信息有限，因此
给出了一种基于可调窗口的短时傅氏变换时频表示方法，进而又利用各匹配子波的复谱，引入一种新的时频表示
方法，不仅与 Wigner-Ville 方法具有同等的分辨率和能量聚集特性，而且保留了原信号的最基本时频特征，不存在
由于 Mallat 提出的匹配追踪时频分析方法是用 各匹配子波的 Wigner-Ville 分布表示的，这种时频 表征方式只能提供原信号的时频能量（ 振幅） 分布， 而没有相位信息，而且每个匹配子波的 Wigner-Ville 分布均存在交叉项，在对所有匹配子波的 WignerVille 分布求和时，交叉项能量由于不能实现同相叠 加而受到压制。笔者根据匹配追踪的算法特点，提 出了一种基于可调窗口的 STFT 时频表示方法，通 过逐步递减的窗口长度克服了常规短时傅氏变换单 一分辨率的缺点，而且此方法本身就不产生交叉项，
a—利用式（ 10） 得到的振幅时频分布； b—利用式（ 10） 得到的相位时频分布
图 7 合成地震记录时频表征
3 匹配追踪时频表征与 S 变换的对比
相位雷克子波组成。 对图 8a 的合成地震记录分别进行 S 变换和用
波与原信号越相似，则分解效果越好，由于雷克子波
具有理想的地震子波波形，在对地震信号进行匹配
追踪分解时，可以利用雷克子波作为基本子波来创 建子波库［3］。
设地震信号为 s（ t） ，由式（ 1） 可知，匹配追踪分
解成雷克子波的组合为
∑ s（ t） = ai ·r（ t － ti ，fi ，φi ）
（ 4）
·549·
图 6 左侧曲线是 50 Hz 的雷克子波包络，上部 曲线是其频谱，二者根据式 （ 7 ） 得到其时 频 分 布。 从图 6 看出，匹配子波的包络和频谱分别控制了其 在时域和频域的分布范围，具有较高的能量聚焦和 时频分辨率。与式（ 5） 不同，式（ 7） 给出的是最基本 的信号时频表征，可以得到时频分布的实部、虚部、 振幅、相位等多种信息。而且与式（ 6） 相比，式（ 7）
图 1 不同主频雷克子波及其合成记录
图 2 地震记录与匹配追踪重构结果的比较
·548·
物探与化探
35 卷
图 3 图 1 地震记录的匹配追踪振幅时频分布
图 4 图 1 地震记录的 WVD 振幅时频分布
a—由式（ 6） 得到图 1 地震信号的时频分布； b—直接对图 1 地震信号进行 STFT 得到的时频分布
没有引入任何窗口，因而不存在窗口截断效应。 仍然以图 1 的地震信号来验证式（ 7） 给出的时
频分布。图 7a、b 分别给出了由式（ 7） 计算出的地 震信号振幅和相位的时频分布。对比图 7a 与图 3 可以发现，利用式（ 7） 得到的时频分布，无论是在时 频分辨率还是在能量聚集性上可以与式（ 5） 的计算 结果相媲美。
选择 wγn ，使 |〈Rn s，wγn 〉| 最大，则残差‖Rn + 1 s‖2 最
小。因此，方程（ 2） 中匹配子波的振幅满足如下条件
｛ an ，γn ｝
= arg max ‖Rn+1 s‖2 。 a n，γ n∈Γ
（ 3）
1． 2 基于雷克子波的匹配追踪
匹配追踪算 法 中，创 建 子 波 库 所 使 用 的 基 本 子
主频的匹配子波来说，其短时傅氏变换窗口是不变
的，因此该时频表征方法对原地震信号的正反变换
都不会产生畸变。这样，该方法将匹配追踪得到的
匹配子波与短时傅氏变换结合，形成了变窗口短时
傅氏变换时频表示，而且不存在交叉项干扰。
傅氏变换，但仍然难免窗口的截断效应。
为了避免窗口截断效应，并且得到与相位有关
的信息，根据匹配追踪的算法特点，通过对每个匹配
进而笔者引入一种无窗口截断效应的时频表示方 法，通过匹配子波沿时间轴的能量包络和沿频率轴 的振幅包络双重控制，不仅能够保持 Mallat 方法时 频表示的分辨率，而且不产生交叉项。
1 匹配追踪基本原理
匹配追踪是一种基于子波库扫描的信号自适应
分解方法。设 D 为进行信号分解的超完备子波库，
信号为 s（ t） 。D 中的元素满足 D = ｛ wγ ∶ γ ∈ Γ｝ ， wγ 为其中的一个匹配子波且‖wγ ‖ = 1。匹配追踪 算法每 次 迭 代 都 自 适 应 地 提 取 出 一 个 最 优 子 波
N －1
∑ TFR1 （ t，f） = an 槡WVD［wγn（ t，f） ］ （ 5） n =0
其中，WVD［wγn（ t，f） ］表示匹配子波 wγn （ t） 的 Wigner-Ville 分布。
图 3 为利用式（ 5） 对图 1 最左边的合成地震记 录计算得到的匹配追踪时频分布。从图中可以看 出，合成该地震记录的各子波分量在时频图中均有 体现，而且在 1． 5 s 处相隔较近的两个子波也可以 分开，几乎没有交叉项能量的干扰。图 4 是直接对 图 1 最左边的合成地震信号求 Wigner-Ville 分布的 结果，可见存在很强的交叉项能量，严重影响信号的 时频表示结果。 2． 2 基于短时傅氏变换的匹配追踪时频表征
的中心时间、主频和相位，称为子波的控制参数。每
次迭代中，最优子波 wγn 的控制参数的确定是整个 匹配追踪算法的核心。
1． 1 匹配子波控制参数的优选
对于一道地震信号，利用 Hilbert 变换将其构造
成一个复信号，将复地震信号包络最大值处的时间
作为 tn 的初值，最大值处对应的瞬时频率作为 fn 的 初值，最大值处的瞬时相位作为 φn 的初值。设 tn 的初值为 t0n，采样间隔为 Δtn，则其扫描范围为
收稿日期：2010 － 06 － 19 基金项目：国家 863（ 2006AA09A102） 资助项目
ቤተ መጻሕፍቲ ባይዱ期
张繁昌等： 非平稳地震信号匹配追踪时频分析
·547·
［t0n － m·Δtn ，t0n + m·Δtn ］， m 为正整数，用来控制扫描范围。用同样的方法可
以设置 fn 和 φn 的扫描范围。通过对每个控制参数 在其初始值附近进行扫描，可以确定出每个控制参
wγn，n 为迭代次数。经 N 次迭代之后，
N －1
∑ s（ t） = anwγn（ t） + RNs ，
（ 1）
n =0
其中，an 是第 n 个匹配子波 wγn 的振幅，是迭代 N
次之后的残差，且 R（ 0） s = s。匹配子波 wγn 可以由 γn
= ｛ tn，fn，φn｝ 来刻画，其中 tn、fn 和 φn 分别为子波
数的最优值。tn，fn，φn 三个控制参数一旦确定，匹
配子波 wγn也就随之确定，然后利用
an = | 〈Rn s，wγn 〉|
（ 2）
来确定其振 幅 其 中，〈，〉表 示 内 积 运 算。 在 迭 代 分
解过程中，由于 wγn与 Rn + 1 s 正交，因此有
‖Rn s‖ = | 〈Rn s，wγn 〉| 2 + ‖Rn s‖2 ，
i
其中，ai 、ti 、fi 、φi 分别表示第 i 个雷克子波的振幅、
中心时间、主频和相位。
为了验证基于雷克子波的匹配追踪分解效果，
合成了图 1 所示的地震记录（ 图 1 左边第 1 道） ，该
记录由不同主频和相位的雷克子波叠加在一起组
成，其中包括 2 个零相位的 10 Hz 子波（ 图 1 左边第
短时傅氏变换不产生交叉项，但因其窗口不可 调而使得其时频表示的分辨率受到限制。考虑到各 匹配子波的独立性，笔者提出如下的匹配追踪时频 表征方法： 对各匹配子波进行短时傅氏变换时，对低 频匹配子波采用较大时窗，而对高频匹配子波则采 用较小时窗，这样便可以克服常规短时傅氏变换窗 口固定、分辨率单一的缺点，得到具有较高分辨率的
子波的复谱求和，可以得到如下时频分布表达式
∞
∑ TFR3（ t，f） = anWγn（ f） env［wγn（ t － tn） ］ei2πf（ t－tn） ，（ 7） －∞
其中，Wγn（ f） 为匹配子波 wγn （ t） 的傅氏谱，env［wγn （ t － tn） ］为匹配子波的瞬时包络（ 图 6） 。
基于这种思路，对图 1 所示的合成地震记录采
用长度逐步递减的时窗（ 100 ～ 10 ms） ，分别对从低
频到高频对匹配子波依照式（ 6） 求其时频分布，得
到图 5a 所示的结果。与直接对图 1 的地震信号进
行 100 ms 固定窗口短时傅氏变换相比（ 图 5b） ，时
频分辨率明显改善。但将图 5a 与图 3 对比可见，本
方法的分辨率不如利用式（ 5） 的分辨率高。
2． 3 无窗口截断效应的匹配追踪时频表征
匹配追踪时频分布公式 （ 5 ） 可以保持 Wigner-
Ville 分布的高分辨率，并且压制交叉项的影响，但
它只能给出信号的振幅分布特征，无法从中得到瞬
图 6 50 Hz 雷克子波的时频分布
4期
张繁昌等： 非平稳地震信号匹配追踪时频分析
匹配追踪（ MP） 时频分析方法作为新一代谱分 解技术，在地震解释、储层识别及烃类检测等方面具 有很大的应用潜力［1］。匹配追踪方法由 Mallat 等人 首次提出，其原始算法是基于 Gabor 函数组成的超 完备子波库贪婪算法［2］，并以子波库中各匹配子波 的 Wigner-Ville 分 布 之 和 来 对 原 信 号 进 行 时 频 表 征。Liu J 和 Marfurt 在 2004 和 2005 年先后提出了 基于雷 克 子 波 库 和 Morlet 子 波 库 的 匹 配 追 踪 算 法［3 － 4］，并将地震信号的瞬时特征引入到匹配追踪 算法当中，提高了算法的分解效率； Wang Y 在 Liu J 算法的基础上，提出了三步法匹配追踪方法，其时频 表示方法仍然是 Wigner-Ville 分布［5］。笔者在三步 法原理的基础上，实现了基于雷克子波的匹配追踪 分解，并与短时傅氏变换（ STFT） ［6］和目前流行的 S 变换（ ST） ［7］进行对比，证明匹配追踪时频表征是一 种非常理想的地震信号时频分析工具。
2 道） ，4 个零相位的 30 Hz 子波（ 图 1 左边第 3 道） ，
2 个 90°相位的 50 Hz 子波（ 图 1 右边第 2 道） ，1 个
45°相位的 30 Hz 子波（ 图 1 右边第 1 道） ，各个子波
的振幅相同，记录长度为 1． 6 s，2 ms 采样。将此合
成地震记录进行基于雷克子波的匹配追踪分解，并 利用方程（ 4） 进行重构，其结果如图 2 所示，其中最 左边 1 道是图 1 最左边的原始地震记录，中间的 1 道是匹配追踪重构结果，最右边是前二者的差。通 过图 2 的对比可见，该地震道可以完全由匹配子波 重构，残差曲线接近零值。