基于改進(jìn)匹配追蹤算法的化爆地震波信號時頻特征提取＊

鐘明壽，周 輝，劉 影，龍 源，郭 濤

基于改進(jìn)匹配追蹤算法的化爆地震波信號時頻特征提取＊

鐘明壽，周 輝，劉 影，龍 源，郭 濤

（解放軍理工大學(xué)野戰(zhàn)工程學(xué)院，江蘇 南京210007）

通過化爆地震波信號的時頻分析可以獲取豐富的地下巖層信息，為進(jìn)一步的化爆地震波特性研究提供支撐。本文中利用爆炸地震波信號的非平穩(wěn)隨機(jī)特性，提出一種改進(jìn)的匹配追蹤算法（matching pursuits），該方法能更有效地獲取化爆地震波信號的時頻信息。該算法首先對地震波信號進(jìn)行Hibert變換，將其轉(zhuǎn)換成復(fù)數(shù)信號，獲得爆炸地震波信號瞬時頻率相位參數(shù)，再進(jìn)行子波分解，從而顯著提高了匹配追蹤算法的運(yùn)算速率；用分解后的信號計算其魏格納威利分布（Wigner Ville distribution），有效地消除了交叉干擾項的影響。將該方法用于實測地震波信號的時頻分析，獲得了分辨率較高的時頻分布圖，且運(yùn)算速率大大提高。

爆炸地震波；匹配追蹤；復(fù)數(shù)信號；時頻分析

化爆震源地震勘探是目前油氣勘探方法中應(yīng)用最廣泛、效果最顯著的重要技術(shù)之一。由于地面采集到的爆炸地震波信號是地下多種地質(zhì)體的綜合響應(yīng)，是一種復(fù)合諧波。通過對爆炸地震波相關(guān)特征參數(shù)的分析，可為確定石油和天然氣的存在位置提供重要參考，從而為進(jìn)一步開采提供決策依據(jù)［1～3］。

信號的時頻分析一直是學(xué)術(shù)界和工程界實現(xiàn)特征提取的重要工具和手段之一，通過時頻分析得到兩個主要的地震相特征參數(shù)：振幅和頻率，它將信號從時域變換到頻域，可以將信號中與頻率密切相關(guān)的特征信息提取出來，為后續(xù)的物理現(xiàn)象識別和工程問題的判別、決策提供依據(jù)。信號時頻分析方法特別多，但都受到Heisenberg測不準(zhǔn)原理的限制。S.Mallat等［4］提出的匹配追蹤法（MP），是將信號按字典原子逐步分解，得到用字典中的原子表示的分解信號，再計算其魏格納威利分布（WVD），從而獲取分解信號時頻信息［5］［6］。該方法屬于二次型時頻分析，它有效消除了 WVD交叉干擾項［7］的影響，能獲取信號更高分辨率的時頻信息。但該方法也存在計算量大的不足，使得運(yùn)算速率慢，制約了其對大數(shù)據(jù)量的地震信號的處理能力。

因為地震波信號不可避免地含有各種噪聲，所以傳統(tǒng)的瞬時頻率提取方法都需要用到微分運(yùn)算，而微分運(yùn)算對噪聲敏感。由于Gaussian函數(shù)具有良好的時頻聚集性，目前常用Gabor子波函數(shù)匹配追蹤分解地震波信號［8］。本文中通過在傳統(tǒng)基于Gabor原子的匹配追蹤算法中引入復(fù)信號計算，提出一種新的算法對地震波信號進(jìn)行分解，并對合成數(shù)據(jù)和實測地震波信號進(jìn)行時頻分析，該方法可為地震勘探地震波信號時頻特征提取提供一種新的有效分析方法。

1 MP算法原理

匹配追蹤算法是一種將信號按字典原子逐步分解，得到用字典中原子表示的分解信號的方法［9］，屬于稀疏分解范疇。該方法與統(tǒng)計學(xué)中使用的投影追蹤算法和波形增益矢量量化有密切關(guān)系，由于數(shù)學(xué)原理上的相近，匹配追蹤算法可以直接使用L.K.Jones證明投影追蹤算法收斂的結(jié)論。

設(shè)D為信號分解的字典，f為待分解信號，其長度為N，則由定義可知：

設(shè)gγ0∈D ，則f 可以表示為：

式中：Rf表示差值。為了使得差值最小，則gγ0與Rf正交，因此：

當(dāng)進(jìn)行了n次迭代（n?0）后，得到差值Rnf，再對其進(jìn)行分解，即

式中：Rn＋1f為進(jìn)行了第n＋1次迭代后的差值。重復(fù)上述步驟，則迭代m 次后，原始信號f被分解為如下形式：

若選取的子波是收斂的，則經(jīng)過多次迭代后子波信號就能精確表示原始信號。雖然MP算法是屬于非線性迭代過程，產(chǎn)生的子波序列并不一定收斂，考慮到原始分解信號為地震波信號，其能量有限的特點保證了子波的收斂性。

2 改進(jìn)算法的時頻分析

傳統(tǒng)的選用Gabor基函數(shù)的MP算法分解信號時，先將此基函數(shù)的時頻參數(shù)進(jìn)行離散化，形成過完備的原子庫，再將待分解信號按此原子庫進(jìn)行分解，當(dāng)分解后的差值信號滿足一定的條件時結(jié)束分解。本文中在此基礎(chǔ)上引進(jìn)復(fù)信號表示方法，提前提取信號瞬時頻率和相位信息，并進(jìn)行搜索，避免對過于龐大的原子庫進(jìn)行搜索，從而提高運(yùn)算速率。

2.1 復(fù)數(shù)信號的特征提取

設(shè)X（t）是輸入的實信號，對X（t）作Hilbert變換：

式中：h（t）為希氏變換因子，因此構(gòu)造復(fù)數(shù)解析信號z（t）：

其中A（t）為幅值函數(shù)：

φ（t）為相位函數(shù)：

求出瞬時頻率

2.2 復(fù)信號子波分解的算法步驟

對于給定信號，利用Hilbert變換可以提前提取其瞬時頻率和相位信息，基于此可以有效縮小信號在原子庫中的匹配范圍，極大地減少計算量。具體的算法步驟可歸納如下：

（1）選定Gabor原子并對其進(jìn)行擴(kuò)展，形成過完備原子庫Di（i＝1，2，…，I）；

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

（3）對差值信號rm（m＝0，1，2，…，M－1）進(jìn)行 Hilbert變換，獲取其瞬時頻率信息和瞬時相位信息，并將獲得的頻率相位信息帶入過完備原子庫中，選定將進(jìn)行子波分解的原子庫范圍；

（4）從限定的原子庫中找出與差值信號rm最匹配的原子，即內(nèi)積最大的原子dmi，并求出匹配系數(shù)cmi，在此基礎(chǔ)上將差值信號減去匹配原子，得到新的差值信號rm＋1；

（5）重復(fù)步驟（3）和步驟（4），直到差值信號滿足需要的條件。

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

2.3 與傳統(tǒng)方法的比較

對于一給定量綱一函數(shù)

其波形如圖1所示。

分別用傳統(tǒng)方法和改進(jìn)方法對其進(jìn)行分解，再給定閾值小于原始波形兩個數(shù)量級的條件下，分解后的差值如圖2所示。分別計算差值均值，可以得到改進(jìn)后的差值均值要小于傳統(tǒng)方法分解后的差值，且在進(jìn)行波形分解時，改進(jìn)算法運(yùn)算速率明顯快于傳統(tǒng)方法，對大數(shù)據(jù)量的地震波信號進(jìn)行處理時有較大優(yōu)勢。

2.4 時頻分析

時頻分析是將信號從時域變換到頻域，可以將信號中與頻率密切相關(guān)的特征信息提取出來，為后續(xù)的物理現(xiàn)象識別和工程問題的判別、決策提供依據(jù)［10］。魏格納威利分布（WVD）［11］是在定性方面不同于頻譜圖的一些分布的原型，屬于二次型時頻分析，發(fā)現(xiàn)它的長處和短處已經(jīng)成為這個領(lǐng)域研究的主要動向。

信號X（t）的 WVD如下：

由WVD的表達(dá)式可以看出，WVD時頻分析不含窗函數(shù)，不受Heisenberg測不準(zhǔn)原理的限制，理論上其可以獲得較其他時頻分析方法分辨率更高的信號時頻分布。但WVD存在非線性的不足，即兩單個信號 WVD與其合成信號的 WVD并不相同。若令X（t）＝X1（t）＋X2（t），則

式中：2Re［Wx1＋x2（t，ω）］是X1（t）和X2（t）相互影響的結(jié)果，稱之為交叉干擾項。

由公式（12）可以看到：對兩個信號進(jìn)行的WVD，有時這些值是第一個信號與第二個信號的乘積，因此產(chǎn)生了交叉干擾項。這也是二次型時頻分布的不足。如何有效地消除交叉干擾項的影響已經(jīng)成為二次型時頻分析領(lǐng)域重要的研究內(nèi)容。結(jié)合匹配追蹤算法將信號進(jìn)行分解后再進(jìn)行WVD計算，能很好地消除交叉干擾項。為了更直觀地結(jié)合匹配追蹤算法的WVD計算消除交叉干擾項的過程，選用調(diào)制信號進(jìn)行時頻分析，調(diào)制信號如圖3所示。

對調(diào)制信號進(jìn)行直接WVD計算和匹配追蹤分解后的WVD計算，得到的時頻分布圖如圖4所示。可以明顯地看到交叉干擾項，存在6個交叉項（其中兩個發(fā)生了重疊）。

3 合成數(shù)據(jù)計算

為了驗證該方法對提取信號時頻特性的有效性，本文中進(jìn)行了合成信號試算。合成信號（由兩種頻率和兩種相位組成）如圖5所示。

分別直接對合成信號進(jìn)行短時傅里葉變換（short time Fourier transform）獲取時頻特性和進(jìn)行匹配追蹤分解后計算其WVD，結(jié)果如下圖6所示。由圖6可知，本文中算法較好地實現(xiàn)了信號瞬時時頻特性提取，兩種頻率區(qū)分較為明顯；由于此方法沒有窗函數(shù)，其分辨率精度明顯優(yōu)于短時傅里葉變換（STFT）。

在此基礎(chǔ)上加入兩個同頻率的Gabor原子，再進(jìn)行信號的時頻特性提取，結(jié)果如圖7所示。由圖7可知，加入Gabor原子后，本文中的方法仍能很好地提取信號的瞬時時頻特性。

4 實測地震波信號譜分析

以南通市的海安縣境內(nèi)三維地震勘探中化爆地震波監(jiān)測實例分析該方法在實際化爆地震波特性分析中的應(yīng)用。由于地震勘探工區(qū)居民設(shè)施密集，在地震勘探前進(jìn)行了化爆地震波震動測試實驗，測試儀器為成都中科測控有限公司生產(chǎn)的TC－4850型爆破測震儀。測試地點為工地營區(qū)，因工區(qū)大部分為村民民房，因此測點分別布置在大隊營房的一層。現(xiàn)場共進(jìn)行了4炮測試實驗：為裝藥1kg，分別距離測點38、78、102、112m。

圖8為現(xiàn)場爆破震動監(jiān)測布置。

以第2炮監(jiān)測情況為例，利用改進(jìn)匹配追蹤算法分提取其時頻特性，測得化爆地震波信號如圖9所示。對其進(jìn)行基于復(fù)數(shù)信號的匹配追蹤分解，分解結(jié)果如圖10所示。

由圖9可知差值信號已低于原始信號一個數(shù)量級，且在2s后的原始信號的波動應(yīng)為噪聲，重建信號很好地達(dá)到了降噪的效果，結(jié)果較為理想。分別對其進(jìn)行短時傅里葉變換和做WVD計算，結(jié)果如圖11所示。

由圖11可以看出，短時傅里葉變換獲取的爆炸地震波信號分辨率很低，而子波分解能很好地提取地震波信號時頻特性，可以看到此地震波信號能量主要集中在10～30Hz的頻段上，區(qū)分度較高，為進(jìn)一步的化爆地震波分析提供支撐。分別對不同起爆點起爆的地震波進(jìn)行WVD計算，得到地震波時頻特性如圖12所示。

由圖12可知化爆地震波能量多集中在100Hz以內(nèi)的低頻成分，且隨著傳播距離增大，化爆地震波高頻分量衰減迅速。改進(jìn)的匹配追蹤算法能較好地提取較高分辨率的化爆地震波時頻分布，為進(jìn)一步的化爆地震波特性研究提供支撐。

5 結(jié) 論

（1）結(jié)合化爆地震波特性預(yù)先提取其瞬時時頻特征的改進(jìn)匹配追蹤算法在保證分解精度的同時能有效提高運(yùn)算速率，較傳統(tǒng)匹配追蹤算法更適合處理大數(shù)據(jù)量的化爆地震波信號；

（2）基于改進(jìn)匹配追蹤算法進(jìn)行的信號WVD計算能有效消除交叉干擾項，很好地利用了二次型時頻分析不受Heisenberg測不準(zhǔn)原理限制的優(yōu)勢，使得獲得的時頻分布具有更高的分辨率和精確度；

（3）對實際地震勘探的化爆地震波信號進(jìn)行時頻分析可以看到，改進(jìn)匹配追蹤算法能獲得時頻分辨率較高的化爆地震波信號時頻分布，為進(jìn)一步的化爆地震波分析提供支撐。

［1］ 鐘明壽，龍源，謝全民，等.裝藥不耦合系數(shù)對碳酸鹽巖爆炸地震波能量的影響［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］ 鐘明壽，龍源，李興華，等.基于炮孔不同耦合介質(zhì)的孔壁爆炸載荷及比能時間函數(shù)分析［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］ 鐘明壽，龍源，李興華，等.碳酸鹽巖中炮孔不耦合系數(shù)對地震激發(fā)效果影響的分析［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］.計算機(jī)技術(shù)與發(fā)展，2010，20（7）：231－234.Liu Jicheng，F(xiàn)u Shuang.Fast spectral analysis of seismic signal based on matching pursuits algorithm［J］.Computer Technology and Development，2010，20（7）：231－234.

［7］ 趙中華，王文延.基于時頻分析的去除魏格納交叉干擾項的方法［J］.北京交通大學(xué)學(xué)報，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］ 陳發(fā)宇，楊長春.基于 MP方法的地震信號快速分解算法［J］.地球物理學(xué)進(jìn)展，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］ 陳學(xué)華.時頻分布與地震信號譜分析研究［D］.成都：成都理工大學(xué)，2006.

［11］ 王培茂.地震信號時頻特征表示表示方法及應(yīng)用［D］.長春：吉林大學(xué)，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 國標(biāo)學(xué)科代碼：13035

A

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

2016－04－15；

2016－07－27

國家自然科學(xué)基金項目（51304218，51508569）；江蘇省自然科學(xué)基金面上項目（BK20151449）

鐘明壽（1983— ），男，博士，講師；通信作者：周 輝，495994500＠qq.com。

（責(zé)任編輯 曾月蓉）