基于改进匹配追踪算法的化爆地震波信号时频特征提取＊

钟明寿，周 辉，刘 影，龙 源，郭 涛

基于改进匹配追踪算法的化爆地震波信号时频特征提取＊

钟明寿，周 辉，刘 影，龙 源，郭 涛

（解放军理工大学野战工程学院，江苏 南京210007）

通过化爆地震波信号的时频分析可以获取丰富的地下岩层信息，为进一步的化爆地震波特性研究提供支撑。本文中利用爆炸地震波信号的非平稳随机特性，提出一种改进的匹配追踪算法（matching pursuits），该方法能更有效地获取化爆地震波信号的时频信息。该算法首先对地震波信号进行Hibert变换，将其转换成复数信号，获得爆炸地震波信号瞬时频率相位参数，再进行子波分解，从而显著提高了匹配追踪算法的运算速率；用分解后的信号计算其魏格纳威利分布（Wigner Ville distribution），有效地消除了交叉干扰项的影响。将该方法用于实测地震波信号的时频分析，获得了分辨率较高的时频分布图，且运算速率大大提高。

爆炸地震波；匹配追踪；复数信号；时频分析

化爆震源地震勘探是目前油气勘探方法中应用最广泛、效果最显著的重要技术之一。由于地面采集到的爆炸地震波信号是地下多种地质体的综合响应，是一种复合谐波。通过对爆炸地震波相关特征参数的分析，可为确定石油和天然气的存在位置提供重要参考，从而为进一步开采提供决策依据［1～3］。

信号的时频分析一直是学术界和工程界实现特征提取的重要工具和手段之一，通过时频分析得到两个主要的地震相特征参数：振幅和频率，它将信号从时域变换到频域，可以将信号中与频率密切相关的特征信息提取出来，为后续的物理现象识别和工程问题的判别、决策提供依据。信号时频分析方法特别多，但都受到Heisenberg测不准原理的限制。S.Mallat等［4］提出的匹配追踪法（MP），是将信号按字典原子逐步分解，得到用字典中的原子表示的分解信号，再计算其魏格纳威利分布（WVD），从而获取分解信号时频信息［5］［6］。该方法属于二次型时频分析，它有效消除了 WVD交叉干扰项［7］的影响，能获取信号更高分辨率的时频信息。但该方法也存在计算量大的不足，使得运算速率慢，制约了其对大数据量的地震信号的处理能力。

因为地震波信号不可避免地含有各种噪声，所以传统的瞬时频率提取方法都需要用到微分运算，而微分运算对噪声敏感。由于Gaussian函数具有良好的时频聚集性，目前常用Gabor子波函数匹配追踪分解地震波信号［8］。本文中通过在传统基于Gabor原子的匹配追踪算法中引入复信号计算，提出一种新的算法对地震波信号进行分解，并对合成数据和实测地震波信号进行时频分析，该方法可为地震勘探地震波信号时频特征提取提供一种新的有效分析方法。

1 MP算法原理

匹配追踪算法是一种将信号按字典原子逐步分解，得到用字典中原子表示的分解信号的方法［9］，属于稀疏分解范畴。该方法与统计学中使用的投影追踪算法和波形增益矢量量化有密切关系，由于数学原理上的相近，匹配追踪算法可以直接使用L.K.Jones证明投影追踪算法收敛的结论。

设D为信号分解的字典，f为待分解信号，其长度为N，则由定义可知：

设gγ0∈D ，则f 可以表示为：

式中：Rf表示差值。为了使得差值最小，则gγ0与Rf正交，因此：

当进行了n次迭代（n≫0）后，得到差值Rnf，再对其进行分解，即

式中：Rn＋1f为进行了第n＋1次迭代后的差值。重复上述步骤，则迭代m 次后，原始信号f被分解为如下形式：

若选取的子波是收敛的，则经过多次迭代后子波信号就能精确表示原始信号。虽然MP算法是属于非线性迭代过程，产生的子波序列并不一定收敛，考虑到原始分解信号为地震波信号，其能量有限的特点保证了子波的收敛性。

2 改进算法的时频分析

传统的选用Gabor基函数的MP算法分解信号时，先将此基函数的时频参数进行离散化，形成过完备的原子库，再将待分解信号按此原子库进行分解，当分解后的差值信号满足一定的条件时结束分解。本文中在此基础上引进复信号表示方法，提前提取信号瞬时频率和相位信息，并进行搜索，避免对过于庞大的原子库进行搜索，从而提高运算速率。

2.1 复数信号的特征提取

设X（t）是输入的实信号，对X（t）作Hilbert变换：

式中：h（t）为希氏变换因子，因此构造复数解析信号z（t）：

其中A（t）为幅值函数：

φ（t）为相位函数：

求出瞬时频率

2.2 复信号子波分解的算法步骤

对于给定信号，利用Hilbert变换可以提前提取其瞬时频率和相位信息，基于此可以有效缩小信号在原子库中的匹配范围，极大地减少计算量。具体的算法步骤可归纳如下：

（1）选定Gabor原子并对其进行扩展，形成过完备原子库Di（i＝1，2，…，I）；

（2）令初始差值信号r0等于待分解的原始信号X（t）；

（3）对差值信号rm（m＝0，1，2，…，M－1）进行 Hilbert变换，获取其瞬时频率信息和瞬时相位信息，并将获得的频率相位信息带入过完备原子库中，选定将进行子波分解的原子库范围；

（4）从限定的原子库中找出与差值信号rm最匹配的原子，即内积最大的原子dmi，并求出匹配系数cmi，在此基础上将差值信号减去匹配原子，得到新的差值信号rm＋1；

（5）重复步骤（3）和步骤（4），直到差值信号满足需要的条件。

至此，完成了对原始信号X（t）的分解：

2.3 与传统方法的比较

对于一给定量纲一函数

其波形如图1所示。

分别用传统方法和改进方法对其进行分解，再给定阈值小于原始波形两个数量级的条件下，分解后的差值如图2所示。分别计算差值均值，可以得到改进后的差值均值要小于传统方法分解后的差值，且在进行波形分解时，改进算法运算速率明显快于传统方法，对大数据量的地震波信号进行处理时有较大优势。

2.4 时频分析

时频分析是将信号从时域变换到频域，可以将信号中与频率密切相关的特征信息提取出来，为后续的物理现象识别和工程问题的判别、决策提供依据［10］。魏格纳威利分布（WVD）［11］是在定性方面不同于频谱图的一些分布的原型，属于二次型时频分析，发现它的长处和短处已经成为这个领域研究的主要动向。

信号X（t）的 WVD如下：

由WVD的表达式可以看出，WVD时频分析不含窗函数，不受Heisenberg测不准原理的限制，理论上其可以获得较其他时频分析方法分辨率更高的信号时频分布。但WVD存在非线性的不足，即两单个信号 WVD与其合成信号的 WVD并不相同。若令X（t）＝X1（t）＋X2（t），则

式中：2Re［Wx1＋x2（t，ω）］是X1（t）和X2（t）相互影响的结果，称之为交叉干扰项。

由公式（12）可以看到：对两个信号进行的WVD，有时这些值是第一个信号与第二个信号的乘积，因此产生了交叉干扰项。这也是二次型时频分布的不足。如何有效地消除交叉干扰项的影响已经成为二次型时频分析领域重要的研究内容。结合匹配追踪算法将信号进行分解后再进行WVD计算，能很好地消除交叉干扰项。为了更直观地结合匹配追踪算法的WVD计算消除交叉干扰项的过程，选用调制信号进行时频分析，调制信号如图3所示。

对调制信号进行直接WVD计算和匹配追踪分解后的WVD计算，得到的时频分布图如图4所示。可以明显地看到交叉干扰项，存在6个交叉项（其中两个发生了重叠）。

3 合成数据计算

为了验证该方法对提取信号时频特性的有效性，本文中进行了合成信号试算。合成信号（由两种频率和两种相位组成）如图5所示。

分别直接对合成信号进行短时傅里叶变换（short time Fourier transform）获取时频特性和进行匹配追踪分解后计算其WVD，结果如下图6所示。由图6可知，本文中算法较好地实现了信号瞬时时频特性提取，两种频率区分较为明显；由于此方法没有窗函数，其分辨率精度明显优于短时傅里叶变换（STFT）。

在此基础上加入两个同频率的Gabor原子，再进行信号的时频特性提取，结果如图7所示。由图7可知，加入Gabor原子后，本文中的方法仍能很好地提取信号的瞬时时频特性。

4 实测地震波信号谱分析

以南通市的海安县境内三维地震勘探中化爆地震波监测实例分析该方法在实际化爆地震波特性分析中的应用。由于地震勘探工区居民设施密集，在地震勘探前进行了化爆地震波震动测试实验，测试仪器为成都中科测控有限公司生产的TC－4850型爆破测震仪。测试地点为工地营区，因工区大部分为村民民房，因此测点分别布置在大队营房的一层。现场共进行了4炮测试实验：为装药1kg，分别距离测点38、78、102、112m。

图8为现场爆破震动监测布置。

以第2炮监测情况为例，利用改进匹配追踪算法分提取其时频特性，测得化爆地震波信号如图9所示。对其进行基于复数信号的匹配追踪分解，分解结果如图10所示。

由图9可知差值信号已低于原始信号一个数量级，且在2s后的原始信号的波动应为噪声，重建信号很好地达到了降噪的效果，结果较为理想。分别对其进行短时傅里叶变换和做WVD计算，结果如图11所示。

由图11可以看出，短时傅里叶变换获取的爆炸地震波信号分辨率很低，而子波分解能很好地提取地震波信号时频特性，可以看到此地震波信号能量主要集中在10～30Hz的频段上，区分度较高，为进一步的化爆地震波分析提供支撑。分别对不同起爆点起爆的地震波进行WVD计算，得到地震波时频特性如图12所示。

由图12可知化爆地震波能量多集中在100Hz以内的低频成分，且随着传播距离增大，化爆地震波高频分量衰减迅速。改进的匹配追踪算法能较好地提取较高分辨率的化爆地震波时频分布，为进一步的化爆地震波特性研究提供支撑。

5 结 论

（1）结合化爆地震波特性预先提取其瞬时时频特征的改进匹配追踪算法在保证分解精度的同时能有效提高运算速率，较传统匹配追踪算法更适合处理大数据量的化爆地震波信号；

（2）基于改进匹配追踪算法进行的信号WVD计算能有效消除交叉干扰项，很好地利用了二次型时频分析不受Heisenberg测不准原理限制的优势，使得获得的时频分布具有更高的分辨率和精确度；

（3）对实际地震勘探的化爆地震波信号进行时频分析可以看到，改进匹配追踪算法能获得时频分辨率较高的化爆地震波信号时频分布，为进一步的化爆地震波分析提供支撑。

［1］ 钟明寿，龙源，谢全民，等.装药不耦合系数对碳酸盐岩爆炸地震波能量的影响［J］.爆炸与冲击，2011，31（6）：612－618.Zhong Mingshou，Long Yuan，Xie Quanmin，et al.Effects of non－coupling charge coefficients on explosion seismic wave energy in carbonate rocks［J］.Explosion and Shocks Waves，2011，31（6）：612－618.

［2］ 钟明寿，龙源，李兴华，等.基于炮孔不同耦合介质的孔壁爆炸载荷及比能时间函数分析［J］.振动与冲击，2011，30（7）：116－119.Zhong Mingshou，Long Yuan，Li Xinghua，et al.Time function for borehole explosive loading and specific energy based on different coupling mediums［J］.Journal of Vibration and Shock，2011，30（7）：116－119.

［3］ 钟明寿，龙源，李兴华，等.碳酸盐岩中炮孔不耦合系数对地震激发效果影响的分析［J］.石油地球物理勘探，2011，46（2）：165－169.Zhong Mingshou，Long Yuan，Li Xinghua，et al.Influence of borehole non－coupling coefficients on detonation effectiveness of seismic source in carbonate rock［J］.Oil Geophysical Prospecting，2011，46（2）：165－169.

［4］ Mallat S，Zhang Z.Matching pursuit with time－frequency dictionaries［J］.IEEE Transactiong on Signal Processing，1993，41（12）：3397－3415.

［5］ Wang Y H.Seismic time－frequency spectral decomposition by matching pursuit［J］.Geophysics，2007，72（1）：13－20.

［6］ 刘继承，富爽.基于 MP算法的快速地震信号谱分析［J］.计算机技术与发展，2010，20（7）：231－234.Liu Jicheng，Fu Shuang.Fast spectral analysis of seismic signal based on matching pursuits algorithm［J］.Computer Technology and Development，2010，20（7）：231－234.

［7］ 赵中华，王文延.基于时频分析的去除魏格纳交叉干扰项的方法［J］.北京交通大学学报，2010，34（6）：81－85.Zhao Zhonghua，Wang Wenyan.Simulation of method to eliminate cross－term in Wigner distribution based on matching pursuits in time－frequency domain［J］.Journal of Beijing Jiaotong University，2010，34（6）：81－85.

［8］ 王纯伟，杨胜利.基于地震信号的匹配追踪算法［J］.科教前沿，2010，7（6）：443；460.Wang Chunwei，Yang Shengli.Based on the matching pursuit algorithm of seismic signal［J］.Science and Education，2010，7（6）：443；460.

［9］ 陈发宇，杨长春.基于 MP方法的地震信号快速分解算法［J］.地球物理学进展，2007，22（6）：1692－1697.Chen Fayu，Yang Chanchun.Seismic signal’s decomposition based on matching pursuit method ［J］.Progress in Geophysics，2007，22（6）：1692－1697.

［10］ 陈学华.时频分布与地震信号谱分析研究［D］.成都：成都理工大学，2006.

［11］ 王培茂.地震信号时频特征表示表示方法及应用［D］.长春：吉林大学，2008.

Time－frequency analysis of explosion seismic signal based on improved matching pursuit

Zhong Mingshou，Zhou Hui，Liu Ying，Long Yuan，Guo Tao
（College of Field Engineering，PLA University of Science ＆ Technology，Nanjing210007，Jiangsu，China）

Abundant underground rock information can be obtained through the analysis of time－frequency of explosive seismic wave signals，which lays a foundation for further research on characteristics of explosive seismic waves.In this study，based on the non－stationary random nature of the explosive seismic signals，we presented an improved matching pursuits（MP）algorithm to extract the timefrequency information of explosive seismic signals，namely using the Hilbert transformation to convert the signal into a plural signal to acquire the instantaneous frequency and phase parameter of the signal and then decomposing the signal to a series of wavelets.Thus the operation rate was significantly improved and the Wigner Ville distributions（WVD）and its time frequency were calculated with the influence of the cross interference term effectively eliminated.This method has been used to take a time－frequency analysis of an actual explosive seismic wave signal and has acquired a time－frequency distribution with a high resolution and operation rate.

explosive seismic wave；matching pursuits；plural signal；time－frequency analysis

O389 国标学科代码：13035

A

10.11883／1001－1455（2017）06－0931－08

2016－04－15；

2016－07－27

国家自然科学基金项目（51304218，51508569）；江苏省自然科学基金面上项目（BK20151449）

钟明寿（1983— ），男，博士，讲师；通信作者：周 辉，495994500＠qq.com。

（责任编辑 曾月蓉）