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

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

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

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

（解放軍理工大學野戰工程學院，江蘇 南京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。

（責任編輯 曾月蓉）