基于改進(jìn)匹配追蹤算法的化爆地震波信號(hào)時(shí)頻特征提?。?/h1>
2017-12-21 11:15鐘明壽
爆炸與沖擊 2017年6期
鐘明壽，周 輝，劉 影，龍 源，郭 濤
基于改進(jìn)匹配追蹤算法的化爆地震波信號(hào)時(shí)頻特征提?。?/p>
鐘明壽，周 輝，劉 影，龍 源，郭 濤
（解放軍理工大學(xué)野戰(zhàn)工程學(xué)院，江蘇 南京210007）
通過(guò)化爆地震波信號(hào)的時(shí)頻分析可以獲取豐富的地下巖層信息，為進(jìn)一步的化爆地震波特性研究提供支撐。本文中利用爆炸地震波信號(hào)的非平穩(wěn)隨機(jī)特性，提出一種改進(jìn)的匹配追蹤算法（matching pursuits），該方法能更有效地獲取化爆地震波信號(hào)的時(shí)頻信息。該算法首先對(duì)地震波信號(hào)進(jìn)行Hibert變換，將其轉(zhuǎn)換成復(fù)數(shù)信號(hào)，獲得爆炸地震波信號(hào)瞬時(shí)頻率相位參數(shù)，再進(jìn)行子波分解，從而顯著提高了匹配追蹤算法的運(yùn)算速率；用分解后的信號(hào)計(jì)算其魏格納威利分布（Wigner Ville distribution），有效地消除了交叉干擾項(xiàng)的影響。將該方法用于實(shí)測(cè)地震波信號(hào)的時(shí)頻分析，獲得了分辨率較高的時(shí)頻分布圖，且運(yùn)算速率大大提高。
爆炸地震波；匹配追蹤；復(fù)數(shù)信號(hào)；時(shí)頻分析
化爆震源地震勘探是目前油氣勘探方法中應(yīng)用最廣泛、效果最顯著的重要技術(shù)之一。由于地面采集到的爆炸地震波信號(hào)是地下多種地質(zhì)體的綜合響應(yīng)，是一種復(fù)合諧波。通過(guò)對(duì)爆炸地震波相關(guān)特征參數(shù)的分析，可為確定石油和天然氣的存在位置提供重要參考，從而為進(jìn)一步開采提供決策依據(jù)［1～3］。
信號(hào)的時(shí)頻分析一直是學(xué)術(shù)界和工程界實(shí)現(xiàn)特征提取的重要工具和手段之一，通過(guò)時(shí)頻分析得到兩個(gè)主要的地震相特征參數(shù)：振幅和頻率，它將信號(hào)從時(shí)域變換到頻域，可以將信號(hào)中與頻率密切相關(guān)的特征信息提取出來(lái)，為后續(xù)的物理現(xiàn)象識(shí)別和工程問(wèn)題的判別、決策提供依據(jù)。信號(hào)時(shí)頻分析方法特別多，但都受到Heisenberg測(cè)不準(zhǔn)原理的限制。S.Mallat等［4］提出的匹配追蹤法（MP），是將信號(hào)按字典原子逐步分解，得到用字典中的原子表示的分解信號(hào)，再計(jì)算其魏格納威利分布（WVD），從而獲取分解信號(hào)時(shí)頻信息［5］［6］。該方法屬于二次型時(shí)頻分析，它有效消除了 WVD交叉干擾項(xiàng)［7］的影響，能獲取信號(hào)更高分辨率的時(shí)頻信息。但該方法也存在計(jì)算量大的不足，使得運(yùn)算速率慢，制約了其對(duì)大數(shù)據(jù)量的地震信號(hào)的處理能力。
因?yàn)榈卣鸩ㄐ盘?hào)不可避免地含有各種噪聲，所以傳統(tǒng)的瞬時(shí)頻率提取方法都需要用到微分運(yùn)算，而微分運(yùn)算對(duì)噪聲敏感。由于Gaussian函數(shù)具有良好的時(shí)頻聚集性，目前常用Gabor子波函數(shù)匹配追蹤分解地震波信號(hào)［8］。本文中通過(guò)在傳統(tǒng)基于Gabor原子的匹配追蹤算法中引入復(fù)信號(hào)計(jì)算，提出一種新的算法對(duì)地震波信號(hào)進(jìn)行分解，并對(duì)合成數(shù)據(jù)和實(shí)測(cè)地震波信號(hào)進(jìn)行時(shí)頻分析，該方法可為地震勘探地震波信號(hào)時(shí)頻特征提取提供一種新的有效分析方法。
1 MP算法原理
匹配追蹤算法是一種將信號(hào)按字典原子逐步分解，得到用字典中原子表示的分解信號(hào)的方法［9］，屬于稀疏分解范疇。該方法與統(tǒng)計(jì)學(xué)中使用的投影追蹤算法和波形增益矢量量化有密切關(guān)系，由于數(shù)學(xué)原理上的相近，匹配追蹤算法可以直接使用L.K.Jones證明投影追蹤算法收斂的結(jié)論。
設(shè)D為信號(hào)分解的字典，f為待分解信號(hào)，其長(zhǎng)度為N，則由定義可知：
設(shè)gγ0∈D ，則f 可以表示為：
式中：Rf表示差值。為了使得差值最小，則gγ0與Rf正交，因此：
當(dāng)進(jìn)行了n次迭代（n?0）后，得到差值Rnf，再對(duì)其進(jìn)行分解，即
式中：Rn＋1f為進(jìn)行了第n＋1次迭代后的差值。重復(fù)上述步驟，則迭代m 次后，原始信號(hào)f被分解為如下形式：
若選取的子波是收斂的，則經(jīng)過(guò)多次迭代后子波信號(hào)就能精確表示原始信號(hào)。雖然MP算法是屬于非線性迭代過(guò)程，產(chǎn)生的子波序列并不一定收斂，考慮到原始分解信號(hào)為地震波信號(hào)，其能量有限的特點(diǎn)保證了子波的收斂性。
2 改進(jìn)算法的時(shí)頻分析
傳統(tǒng)的選用Gabor基函數(shù)的MP算法分解信號(hào)時(shí)，先將此基函數(shù)的時(shí)頻參數(shù)進(jìn)行離散化，形成過(guò)完備的原子庫(kù)，再將待分解信號(hào)按此原子庫(kù)進(jìn)行分解，當(dāng)分解后的差值信號(hào)滿足一定的條件時(shí)結(jié)束分解。本文中在此基礎(chǔ)上引進(jìn)復(fù)信號(hào)表示方法，提前提取信號(hào)瞬時(shí)頻率和相位信息，并進(jìn)行搜索，避免對(duì)過(guò)于龐大的原子庫(kù)進(jìn)行搜索，從而提高運(yùn)算速率。
2.1 復(fù)數(shù)信號(hào)的特征提取
設(shè)X（t）是輸入的實(shí)信號(hào)，對(duì)X（t）作Hilbert變換：
式中：h（t）為希氏變換因子，因此構(gòu)造復(fù)數(shù)解析信號(hào)z（t）：
其中A（t）為幅值函數(shù)：
φ（t）為相位函數(shù)：
求出瞬時(shí)頻率
2.2 復(fù)信號(hào)子波分解的算法步驟
對(duì)于給定信號(hào)，利用Hilbert變換可以提前提取其瞬時(shí)頻率和相位信息，基于此可以有效縮小信號(hào)在原子庫(kù)中的匹配范圍，極大地減少計(jì)算量。具體的算法步驟可歸納如下：
（1）選定Gabor原子并對(duì)其進(jìn)行擴(kuò)展，形成過(guò)完備原子庫(kù)Di（i＝1，2，…，I）；
（2）令初始差值信號(hào)r0等于待分解的原始信號(hào)X（t）；
（3）對(duì)差值信號(hào)rm（m＝0，1，2，…，M－1）進(jìn)行 Hilbert變換，獲取其瞬時(shí)頻率信息和瞬時(shí)相位信息，并將獲得的頻率相位信息帶入過(guò)完備原子庫(kù)中，選定將進(jìn)行子波分解的原子庫(kù)范圍；
（4）從限定的原子庫(kù)中找出與差值信號(hào)rm最匹配的原子，即內(nèi)積最大的原子dmi，并求出匹配系數(shù)cmi，在此基礎(chǔ)上將差值信號(hào)減去匹配原子，得到新的差值信號(hào)rm＋1；
（5）重復(fù)步驟（3）和步驟（4），直到差值信號(hào)滿足需要的條件。
至此，完成了對(duì)原始信號(hào)X（t）的分解：
2.3 與傳統(tǒng)方法的比較
對(duì)于一給定量綱一函數(shù)
其波形如圖1所示。
分別用傳統(tǒng)方法和改進(jìn)方法對(duì)其進(jìn)行分解，再給定閾值小于原始波形兩個(gè)數(shù)量級(jí)的條件下，分解后的差值如圖2所示。分別計(jì)算差值均值，可以得到改進(jìn)后的差值均值要小于傳統(tǒng)方法分解后的差值，且在進(jìn)行波形分解時(shí)，改進(jìn)算法運(yùn)算速率明顯快于傳統(tǒng)方法，對(duì)大數(shù)據(jù)量的地震波信號(hào)進(jìn)行處理時(shí)有較大優(yōu)勢(shì)。
2.4 時(shí)頻分析
時(shí)頻分析是將信號(hào)從時(shí)域變換到頻域，可以將信號(hào)中與頻率密切相關(guān)的特征信息提取出來(lái)，為后續(xù)的物理現(xiàn)象識(shí)別和工程問(wèn)題的判別、決策提供依據(jù)［10］。魏格納威利分布（WVD）［11］是在定性方面不同于頻譜圖的一些分布的原型，屬于二次型時(shí)頻分析，發(fā)現(xiàn)它的長(zhǎng)處和短處已經(jīng)成為這個(gè)領(lǐng)域研究的主要?jiǎng)酉颉?/p>
信號(hào)X（t）的 WVD如下：
由WVD的表達(dá)式可以看出，WVD時(shí)頻分析不含窗函數(shù)，不受Heisenberg測(cè)不準(zhǔn)原理的限制，理論上其可以獲得較其他時(shí)頻分析方法分辨率更高的信號(hào)時(shí)頻分布。但WVD存在非線性的不足，即兩單個(gè)信號(hào) WVD與其合成信號(hào)的 WVD并不相同。若令X（t）＝X1（t）＋X2（t），則
式中：2Re［Wx1＋x2（t，ω）］是X1（t）和X2（t）相互影響的結(jié)果，稱之為交叉干擾項(xiàng)。
由公式（12）可以看到：對(duì)兩個(gè)信號(hào)進(jìn)行的WVD，有時(shí)這些值是第一個(gè)信號(hào)與第二個(gè)信號(hào)的乘積，因此產(chǎn)生了交叉干擾項(xiàng)。這也是二次型時(shí)頻分布的不足。如何有效地消除交叉干擾項(xiàng)的影響已經(jīng)成為二次型時(shí)頻分析領(lǐng)域重要的研究?jī)?nèi)容。結(jié)合匹配追蹤算法將信號(hào)進(jìn)行分解后再進(jìn)行WVD計(jì)算，能很好地消除交叉干擾項(xiàng)。為了更直觀地結(jié)合匹配追蹤算法的WVD計(jì)算消除交叉干擾項(xiàng)的過(guò)程，選用調(diào)制信號(hào)進(jìn)行時(shí)頻分析，調(diào)制信號(hào)如圖3所示。
對(duì)調(diào)制信號(hào)進(jìn)行直接WVD計(jì)算和匹配追蹤分解后的WVD計(jì)算，得到的時(shí)頻分布圖如圖4所示?？梢悦黠@地看到交叉干擾項(xiàng)，存在6個(gè)交叉項(xiàng)（其中兩個(gè)發(fā)生了重疊）。
3 合成數(shù)據(jù)計(jì)算
為了驗(yàn)證該方法對(duì)提取信號(hào)時(shí)頻特性的有效性，本文中進(jìn)行了合成信號(hào)試算。合成信號(hào)（由兩種頻率和兩種相位組成）如圖5所示。
分別直接對(duì)合成信號(hào)進(jìn)行短時(shí)傅里葉變換（short time Fourier transform）獲取時(shí)頻特性和進(jìn)行匹配追蹤分解后計(jì)算其WVD，結(jié)果如下圖6所示。由圖6可知，本文中算法較好地實(shí)現(xiàn)了信號(hào)瞬時(shí)時(shí)頻特性提取，兩種頻率區(qū)分較為明顯；由于此方法沒(méi)有窗函數(shù)，其分辨率精度明顯優(yōu)于短時(shí)傅里葉變換（STFT）。
在此基礎(chǔ)上加入兩個(gè)同頻率的Gabor原子，再進(jìn)行信號(hào)的時(shí)頻特性提取，結(jié)果如圖7所示。由圖7可知，加入Gabor原子后，本文中的方法仍能很好地提取信號(hào)的瞬時(shí)時(shí)頻特性。
4 實(shí)測(cè)地震波信號(hào)譜分析
以南通市的海安縣境內(nèi)三維地震勘探中化爆地震波監(jiān)測(cè)實(shí)例分析該方法在實(shí)際化爆地震波特性分析中的應(yīng)用。由于地震勘探工區(qū)居民設(shè)施密集，在地震勘探前進(jìn)行了化爆地震波震動(dòng)測(cè)試實(shí)驗(yàn)，測(cè)試儀器為成都中科測(cè)控有限公司生產(chǎn)的TC－4850型爆破測(cè)震儀。測(cè)試地點(diǎn)為工地營(yíng)區(qū)，因工區(qū)大部分為村民民房，因此測(cè)點(diǎn)分別布置在大隊(duì)營(yíng)房的一層。現(xiàn)場(chǎng)共進(jìn)行了4炮測(cè)試實(shí)驗(yàn)：為裝藥1kg，分別距離測(cè)點(diǎn)38、78、102、112m。
圖8為現(xiàn)場(chǎng)爆破震動(dòng)監(jiān)測(cè)布置。
以第2炮監(jiān)測(cè)情況為例，利用改進(jìn)匹配追蹤算法分提取其時(shí)頻特性，測(cè)得化爆地震波信號(hào)如圖9所示。對(duì)其進(jìn)行基于復(fù)數(shù)信號(hào)的匹配追蹤分解，分解結(jié)果如圖10所示。
由圖9可知差值信號(hào)已低于原始信號(hào)一個(gè)數(shù)量級(jí)，且在2s后的原始信號(hào)的波動(dòng)應(yīng)為噪聲，重建信號(hào)很好地達(dá)到了降噪的效果，結(jié)果較為理想。分別對(duì)其進(jìn)行短時(shí)傅里葉變換和做WVD計(jì)算，結(jié)果如圖11所示。
由圖11可以看出，短時(shí)傅里葉變換獲取的爆炸地震波信號(hào)分辨率很低，而子波分解能很好地提取地震波信號(hào)時(shí)頻特性，可以看到此地震波信號(hào)能量主要集中在10～30Hz的頻段上，區(qū)分度較高，為進(jìn)一步的化爆地震波分析提供支撐。分別對(duì)不同起爆點(diǎn)起爆的地震波進(jìn)行WVD計(jì)算，得到地震波時(shí)頻特性如圖12所示。
由圖12可知化爆地震波能量多集中在100Hz以內(nèi)的低頻成分，且隨著傳播距離增大，化爆地震波高頻分量衰減迅速。改進(jìn)的匹配追蹤算法能較好地提取較高分辨率的化爆地震波時(shí)頻分布，為進(jìn)一步的化爆地震波特性研究提供支撐。
5 結(jié) 論
（1）結(jié)合化爆地震波特性預(yù)先提取其瞬時(shí)時(shí)頻特征的改進(jìn)匹配追蹤算法在保證分解精度的同時(shí)能有效提高運(yùn)算速率，較傳統(tǒng)匹配追蹤算法更適合處理大數(shù)據(jù)量的化爆地震波信號(hào)；
（2）基于改進(jìn)匹配追蹤算法進(jìn)行的信號(hào)WVD計(jì)算能有效消除交叉干擾項(xiàng)，很好地利用了二次型時(shí)頻分析不受Heisenberg測(cè)不準(zhǔn)原理限制的優(yōu)勢(shì)，使得獲得的時(shí)頻分布具有更高的分辨率和精確度；
（3）對(duì)實(shí)際地震勘探的化爆地震波信號(hào)進(jìn)行時(shí)頻分析可以看到，改進(jìn)匹配追蹤算法能獲得時(shí)頻分辨率較高的化爆地震波信號(hào)時(shí)頻分布，為進(jìn)一步的化爆地震波分析提供支撐。
［1］ 鐘明壽，龍?jiān)?，謝全民，等.裝藥不耦合系數(shù)對(duì)碳酸鹽巖爆炸地震波能量的影響［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］ 鐘明壽，龍?jiān)?，李興華，等.基于炮孔不同耦合介質(zhì)的孔壁爆炸載荷及比能時(shí)間函數(shù)分析［J］.振動(dòng)與沖擊，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］ 鐘明壽，龍?jiān)?，李興華，等.碳酸鹽巖中炮孔不耦合系數(shù)對(duì)地震激發(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算法的快速地震信號(hào)譜分析［J］.計(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］ 趙中華，王文延.基于時(shí)頻分析的去除魏格納交叉干擾項(xiàng)的方法［J］.北京交通大學(xué)學(xué)報(bào)，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］ 王純偉，楊勝利.基于地震信號(hào)的匹配追蹤算法［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ā)宇，楊長(zhǎng)春.基于 MP方法的地震信號(hào)快速分解算法［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é)華.時(shí)頻分布與地震信號(hào)譜分析研究［D］.成都：成都理工大學(xué)，2006.
［11］ 王培茂.地震信號(hào)時(shí)頻特征表示表示方法及應(yīng)用［D］.長(zhǎng)春：吉林大學(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 國(guó)標(biāo)學(xué)科代碼：13035
A
10.11883／1001－1455（2017）06－0931－08
2016－04－15；
2016－07－27
國(guó)家自然科學(xué)基金項(xiàng)目（51304218，51508569）；江蘇省自然科學(xué)基金面上項(xiàng)目（BK20151449）
鐘明壽（1983— ），男，博士，講師；通信作者：周 輝，495994500＠qq.com。
（責(zé)任編輯 曾月蓉）
猜你喜歡
高階時(shí)頻變換理論與應(yīng)用計(jì)算技術(shù)與自動(dòng)化(2022年1期)2022-04-15 少兒科學(xué)周刊·兒童版(2021年22期)2021-12-11 原子可以結(jié)合嗎？少兒科學(xué)周刊·兒童版(2021年22期)2021-12-11 數(shù)字日照計(jì)和暗筒式日照計(jì)資料對(duì)比分析氣象水文海洋儀器(2021年4期)2021-12-11 帶你認(rèn)識(shí)原子少兒科學(xué)周刊·兒童版(2021年22期)2021-12-11 分?jǐn)?shù)階傅里葉變換改進(jìn)算法在時(shí)頻分析中的應(yīng)用強(qiáng)度與環(huán)境(2021年3期)2021-07-27 紅細(xì)胞壓積與白蛋白差值在繼發(fā)性腹腔感染患者病程中的變化昆明醫(yī)科大學(xué)學(xué)報(bào)(2021年5期)2021-07-22 高聚焦時(shí)頻分析算法研究現(xiàn)代電子技術(shù)(2020年13期)2020-08-07 基于稀疏時(shí)頻分解的空中目標(biāo)微動(dòng)特征分析宇航計(jì)測(cè)技術(shù)(2018年3期)2018-09-08 關(guān)注高中時(shí)代(2017年7期)2018-02-24
鐘明壽，周 輝，劉 影，龍 源，郭 濤
基于改進(jìn)匹配追蹤算法的化爆地震波信號(hào)時(shí)頻特征提?。?/p>
鐘明壽，周 輝，劉 影，龍 源，郭 濤
（解放軍理工大學(xué)野戰(zhàn)工程學(xué)院，江蘇 南京210007）
通過(guò)化爆地震波信號(hào)的時(shí)頻分析可以獲取豐富的地下巖層信息，為進(jìn)一步的化爆地震波特性研究提供支撐。本文中利用爆炸地震波信號(hào)的非平穩(wěn)隨機(jī)特性，提出一種改進(jìn)的匹配追蹤算法（matching pursuits），該方法能更有效地獲取化爆地震波信號(hào)的時(shí)頻信息。該算法首先對(duì)地震波信號(hào)進(jìn)行Hibert變換，將其轉(zhuǎn)換成復(fù)數(shù)信號(hào)，獲得爆炸地震波信號(hào)瞬時(shí)頻率相位參數(shù)，再進(jìn)行子波分解，從而顯著提高了匹配追蹤算法的運(yùn)算速率；用分解后的信號(hào)計(jì)算其魏格納威利分布（Wigner Ville distribution），有效地消除了交叉干擾項(xiàng)的影響。將該方法用于實(shí)測(cè)地震波信號(hào)的時(shí)頻分析，獲得了分辨率較高的時(shí)頻分布圖，且運(yùn)算速率大大提高。
爆炸地震波；匹配追蹤；復(fù)數(shù)信號(hào)；時(shí)頻分析
化爆震源地震勘探是目前油氣勘探方法中應(yīng)用最廣泛、效果最顯著的重要技術(shù)之一。由于地面采集到的爆炸地震波信號(hào)是地下多種地質(zhì)體的綜合響應(yīng)，是一種復(fù)合諧波。通過(guò)對(duì)爆炸地震波相關(guān)特征參數(shù)的分析，可為確定石油和天然氣的存在位置提供重要參考，從而為進(jìn)一步開采提供決策依據(jù)［1～3］。
信號(hào)的時(shí)頻分析一直是學(xué)術(shù)界和工程界實(shí)現(xiàn)特征提取的重要工具和手段之一，通過(guò)時(shí)頻分析得到兩個(gè)主要的地震相特征參數(shù)：振幅和頻率，它將信號(hào)從時(shí)域變換到頻域，可以將信號(hào)中與頻率密切相關(guān)的特征信息提取出來(lái)，為后續(xù)的物理現(xiàn)象識(shí)別和工程問(wèn)題的判別、決策提供依據(jù)。信號(hào)時(shí)頻分析方法特別多，但都受到Heisenberg測(cè)不準(zhǔn)原理的限制。S.Mallat等［4］提出的匹配追蹤法（MP），是將信號(hào)按字典原子逐步分解，得到用字典中的原子表示的分解信號(hào)，再計(jì)算其魏格納威利分布（WVD），從而獲取分解信號(hào)時(shí)頻信息［5］［6］。該方法屬于二次型時(shí)頻分析，它有效消除了 WVD交叉干擾項(xiàng)［7］的影響，能獲取信號(hào)更高分辨率的時(shí)頻信息。但該方法也存在計(jì)算量大的不足，使得運(yùn)算速率慢，制約了其對(duì)大數(shù)據(jù)量的地震信號(hào)的處理能力。
因?yàn)榈卣鸩ㄐ盘?hào)不可避免地含有各種噪聲，所以傳統(tǒng)的瞬時(shí)頻率提取方法都需要用到微分運(yùn)算，而微分運(yùn)算對(duì)噪聲敏感。由于Gaussian函數(shù)具有良好的時(shí)頻聚集性，目前常用Gabor子波函數(shù)匹配追蹤分解地震波信號(hào)［8］。本文中通過(guò)在傳統(tǒng)基于Gabor原子的匹配追蹤算法中引入復(fù)信號(hào)計(jì)算，提出一種新的算法對(duì)地震波信號(hào)進(jìn)行分解，并對(duì)合成數(shù)據(jù)和實(shí)測(cè)地震波信號(hào)進(jìn)行時(shí)頻分析，該方法可為地震勘探地震波信號(hào)時(shí)頻特征提取提供一種新的有效分析方法。
1 MP算法原理
匹配追蹤算法是一種將信號(hào)按字典原子逐步分解，得到用字典中原子表示的分解信號(hào)的方法［9］，屬于稀疏分解范疇。該方法與統(tǒng)計(jì)學(xué)中使用的投影追蹤算法和波形增益矢量量化有密切關(guān)系，由于數(shù)學(xué)原理上的相近，匹配追蹤算法可以直接使用L.K.Jones證明投影追蹤算法收斂的結(jié)論。
設(shè)D為信號(hào)分解的字典，f為待分解信號(hào)，其長(zhǎng)度為N，則由定義可知：
設(shè)gγ0∈D ，則f 可以表示為：
式中：Rf表示差值。為了使得差值最小，則gγ0與Rf正交，因此：
當(dāng)進(jìn)行了n次迭代（n?0）后，得到差值Rnf，再對(duì)其進(jìn)行分解，即
式中：Rn＋1f為進(jìn)行了第n＋1次迭代后的差值。重復(fù)上述步驟，則迭代m 次后，原始信號(hào)f被分解為如下形式：
若選取的子波是收斂的，則經(jīng)過(guò)多次迭代后子波信號(hào)就能精確表示原始信號(hào)。雖然MP算法是屬于非線性迭代過(guò)程，產(chǎn)生的子波序列并不一定收斂，考慮到原始分解信號(hào)為地震波信號(hào)，其能量有限的特點(diǎn)保證了子波的收斂性。
2 改進(jìn)算法的時(shí)頻分析
傳統(tǒng)的選用Gabor基函數(shù)的MP算法分解信號(hào)時(shí)，先將此基函數(shù)的時(shí)頻參數(shù)進(jìn)行離散化，形成過(guò)完備的原子庫(kù)，再將待分解信號(hào)按此原子庫(kù)進(jìn)行分解，當(dāng)分解后的差值信號(hào)滿足一定的條件時(shí)結(jié)束分解。本文中在此基礎(chǔ)上引進(jìn)復(fù)信號(hào)表示方法，提前提取信號(hào)瞬時(shí)頻率和相位信息，并進(jìn)行搜索，避免對(duì)過(guò)于龐大的原子庫(kù)進(jìn)行搜索，從而提高運(yùn)算速率。
2.1 復(fù)數(shù)信號(hào)的特征提取
設(shè)X（t）是輸入的實(shí)信號(hào)，對(duì)X（t）作Hilbert變換：
式中：h（t）為希氏變換因子，因此構(gòu)造復(fù)數(shù)解析信號(hào)z（t）：
其中A（t）為幅值函數(shù)：
φ（t）為相位函數(shù)：
求出瞬時(shí)頻率
2.2 復(fù)信號(hào)子波分解的算法步驟
對(duì)于給定信號(hào)，利用Hilbert變換可以提前提取其瞬時(shí)頻率和相位信息，基于此可以有效縮小信號(hào)在原子庫(kù)中的匹配范圍，極大地減少計(jì)算量。具體的算法步驟可歸納如下：
（1）選定Gabor原子并對(duì)其進(jìn)行擴(kuò)展，形成過(guò)完備原子庫(kù)Di（i＝1，2，…，I）；
（2）令初始差值信號(hào)r0等于待分解的原始信號(hào)X（t）；
（3）對(duì)差值信號(hào)rm（m＝0，1，2，…，M－1）進(jìn)行 Hilbert變換，獲取其瞬時(shí)頻率信息和瞬時(shí)相位信息，并將獲得的頻率相位信息帶入過(guò)完備原子庫(kù)中，選定將進(jìn)行子波分解的原子庫(kù)范圍；
（4）從限定的原子庫(kù)中找出與差值信號(hào)rm最匹配的原子，即內(nèi)積最大的原子dmi，并求出匹配系數(shù)cmi，在此基礎(chǔ)上將差值信號(hào)減去匹配原子，得到新的差值信號(hào)rm＋1；
（5）重復(fù)步驟（3）和步驟（4），直到差值信號(hào)滿足需要的條件。
至此，完成了對(duì)原始信號(hào)X（t）的分解：
2.3 與傳統(tǒng)方法的比較
對(duì)于一給定量綱一函數(shù)
其波形如圖1所示。
分別用傳統(tǒng)方法和改進(jìn)方法對(duì)其進(jìn)行分解，再給定閾值小于原始波形兩個(gè)數(shù)量級(jí)的條件下，分解后的差值如圖2所示。分別計(jì)算差值均值，可以得到改進(jìn)后的差值均值要小于傳統(tǒng)方法分解后的差值，且在進(jìn)行波形分解時(shí)，改進(jìn)算法運(yùn)算速率明顯快于傳統(tǒng)方法，對(duì)大數(shù)據(jù)量的地震波信號(hào)進(jìn)行處理時(shí)有較大優(yōu)勢(shì)。
2.4 時(shí)頻分析
時(shí)頻分析是將信號(hào)從時(shí)域變換到頻域，可以將信號(hào)中與頻率密切相關(guān)的特征信息提取出來(lái)，為后續(xù)的物理現(xiàn)象識(shí)別和工程問(wèn)題的判別、決策提供依據(jù)［10］。魏格納威利分布（WVD）［11］是在定性方面不同于頻譜圖的一些分布的原型，屬于二次型時(shí)頻分析，發(fā)現(xiàn)它的長(zhǎng)處和短處已經(jīng)成為這個(gè)領(lǐng)域研究的主要?jiǎng)酉颉?/p>
信號(hào)X（t）的 WVD如下：
由WVD的表達(dá)式可以看出，WVD時(shí)頻分析不含窗函數(shù)，不受Heisenberg測(cè)不準(zhǔn)原理的限制，理論上其可以獲得較其他時(shí)頻分析方法分辨率更高的信號(hào)時(shí)頻分布。但WVD存在非線性的不足，即兩單個(gè)信號(hào) WVD與其合成信號(hào)的 WVD并不相同。若令X（t）＝X1（t）＋X2（t），則
式中：2Re［Wx1＋x2（t，ω）］是X1（t）和X2（t）相互影響的結(jié)果，稱之為交叉干擾項(xiàng)。
由公式（12）可以看到：對(duì)兩個(gè)信號(hào)進(jìn)行的WVD，有時(shí)這些值是第一個(gè)信號(hào)與第二個(gè)信號(hào)的乘積，因此產(chǎn)生了交叉干擾項(xiàng)。這也是二次型時(shí)頻分布的不足。如何有效地消除交叉干擾項(xiàng)的影響已經(jīng)成為二次型時(shí)頻分析領(lǐng)域重要的研究?jī)?nèi)容。結(jié)合匹配追蹤算法將信號(hào)進(jìn)行分解后再進(jìn)行WVD計(jì)算，能很好地消除交叉干擾項(xiàng)。為了更直觀地結(jié)合匹配追蹤算法的WVD計(jì)算消除交叉干擾項(xiàng)的過(guò)程，選用調(diào)制信號(hào)進(jìn)行時(shí)頻分析，調(diào)制信號(hào)如圖3所示。
對(duì)調(diào)制信號(hào)進(jìn)行直接WVD計(jì)算和匹配追蹤分解后的WVD計(jì)算，得到的時(shí)頻分布圖如圖4所示?？梢悦黠@地看到交叉干擾項(xiàng)，存在6個(gè)交叉項(xiàng)（其中兩個(gè)發(fā)生了重疊）。
3 合成數(shù)據(jù)計(jì)算
為了驗(yàn)證該方法對(duì)提取信號(hào)時(shí)頻特性的有效性，本文中進(jìn)行了合成信號(hào)試算。合成信號(hào)（由兩種頻率和兩種相位組成）如圖5所示。
分別直接對(duì)合成信號(hào)進(jìn)行短時(shí)傅里葉變換（short time Fourier transform）獲取時(shí)頻特性和進(jìn)行匹配追蹤分解后計(jì)算其WVD，結(jié)果如下圖6所示。由圖6可知，本文中算法較好地實(shí)現(xiàn)了信號(hào)瞬時(shí)時(shí)頻特性提取，兩種頻率區(qū)分較為明顯；由于此方法沒(méi)有窗函數(shù)，其分辨率精度明顯優(yōu)于短時(shí)傅里葉變換（STFT）。
在此基礎(chǔ)上加入兩個(gè)同頻率的Gabor原子，再進(jìn)行信號(hào)的時(shí)頻特性提取，結(jié)果如圖7所示。由圖7可知，加入Gabor原子后，本文中的方法仍能很好地提取信號(hào)的瞬時(shí)時(shí)頻特性。
4 實(shí)測(cè)地震波信號(hào)譜分析
以南通市的海安縣境內(nèi)三維地震勘探中化爆地震波監(jiān)測(cè)實(shí)例分析該方法在實(shí)際化爆地震波特性分析中的應(yīng)用。由于地震勘探工區(qū)居民設(shè)施密集，在地震勘探前進(jìn)行了化爆地震波震動(dòng)測(cè)試實(shí)驗(yàn)，測(cè)試儀器為成都中科測(cè)控有限公司生產(chǎn)的TC－4850型爆破測(cè)震儀。測(cè)試地點(diǎn)為工地營(yíng)區(qū)，因工區(qū)大部分為村民民房，因此測(cè)點(diǎn)分別布置在大隊(duì)營(yíng)房的一層。現(xiàn)場(chǎng)共進(jìn)行了4炮測(cè)試實(shí)驗(yàn)：為裝藥1kg，分別距離測(cè)點(diǎn)38、78、102、112m。
圖8為現(xiàn)場(chǎng)爆破震動(dòng)監(jiān)測(cè)布置。
以第2炮監(jiān)測(cè)情況為例，利用改進(jìn)匹配追蹤算法分提取其時(shí)頻特性，測(cè)得化爆地震波信號(hào)如圖9所示。對(duì)其進(jìn)行基于復(fù)數(shù)信號(hào)的匹配追蹤分解，分解結(jié)果如圖10所示。
由圖9可知差值信號(hào)已低于原始信號(hào)一個(gè)數(shù)量級(jí)，且在2s后的原始信號(hào)的波動(dòng)應(yīng)為噪聲，重建信號(hào)很好地達(dá)到了降噪的效果，結(jié)果較為理想。分別對(duì)其進(jìn)行短時(shí)傅里葉變換和做WVD計(jì)算，結(jié)果如圖11所示。
由圖11可以看出，短時(shí)傅里葉變換獲取的爆炸地震波信號(hào)分辨率很低，而子波分解能很好地提取地震波信號(hào)時(shí)頻特性，可以看到此地震波信號(hào)能量主要集中在10～30Hz的頻段上，區(qū)分度較高，為進(jìn)一步的化爆地震波分析提供支撐。分別對(duì)不同起爆點(diǎn)起爆的地震波進(jìn)行WVD計(jì)算，得到地震波時(shí)頻特性如圖12所示。
由圖12可知化爆地震波能量多集中在100Hz以內(nèi)的低頻成分，且隨著傳播距離增大，化爆地震波高頻分量衰減迅速。改進(jìn)的匹配追蹤算法能較好地提取較高分辨率的化爆地震波時(shí)頻分布，為進(jìn)一步的化爆地震波特性研究提供支撐。
5 結(jié) 論
（1）結(jié)合化爆地震波特性預(yù)先提取其瞬時(shí)時(shí)頻特征的改進(jìn)匹配追蹤算法在保證分解精度的同時(shí)能有效提高運(yùn)算速率，較傳統(tǒng)匹配追蹤算法更適合處理大數(shù)據(jù)量的化爆地震波信號(hào)；
（2）基于改進(jìn)匹配追蹤算法進(jìn)行的信號(hào)WVD計(jì)算能有效消除交叉干擾項(xiàng)，很好地利用了二次型時(shí)頻分析不受Heisenberg測(cè)不準(zhǔn)原理限制的優(yōu)勢(shì)，使得獲得的時(shí)頻分布具有更高的分辨率和精確度；
（3）對(duì)實(shí)際地震勘探的化爆地震波信號(hào)進(jìn)行時(shí)頻分析可以看到，改進(jìn)匹配追蹤算法能獲得時(shí)頻分辨率較高的化爆地震波信號(hào)時(shí)頻分布，為進(jìn)一步的化爆地震波分析提供支撐。
［1］ 鐘明壽，龍?jiān)?，謝全民，等.裝藥不耦合系數(shù)對(duì)碳酸鹽巖爆炸地震波能量的影響［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］ 鐘明壽，龍?jiān)?，李興華，等.基于炮孔不同耦合介質(zhì)的孔壁爆炸載荷及比能時(shí)間函數(shù)分析［J］.振動(dòng)與沖擊，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］ 鐘明壽，龍?jiān)?，李興華，等.碳酸鹽巖中炮孔不耦合系數(shù)對(duì)地震激發(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算法的快速地震信號(hào)譜分析［J］.計(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］ 趙中華，王文延.基于時(shí)頻分析的去除魏格納交叉干擾項(xiàng)的方法［J］.北京交通大學(xué)學(xué)報(bào)，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］ 王純偉，楊勝利.基于地震信號(hào)的匹配追蹤算法［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ā)宇，楊長(zhǎng)春.基于 MP方法的地震信號(hào)快速分解算法［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é)華.時(shí)頻分布與地震信號(hào)譜分析研究［D］.成都：成都理工大學(xué)，2006.
［11］ 王培茂.地震信號(hào)時(shí)頻特征表示表示方法及應(yīng)用［D］.長(zhǎng)春：吉林大學(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 國(guó)標(biāo)學(xué)科代碼：13035
A
10.11883／1001－1455（2017）06－0931－08
2016－04－15；
2016－07－27
國(guó)家自然科學(xué)基金項(xiàng)目（51304218，51508569）；江蘇省自然科學(xué)基金面上項(xiàng)目（BK20151449）
鐘明壽（1983— ），男，博士，講師；通信作者：周 輝，495994500＠qq.com。
（責(zé)任編輯 曾月蓉）
猜你喜歡
計(jì)算技術(shù)與自動(dòng)化(2022年1期)2022-04-15
少兒科學(xué)周刊·兒童版(2021年22期)2021-12-11
氣象水文海洋儀器(2021年4期)2021-12-11
少兒科學(xué)周刊·兒童版(2021年22期)2021-12-11
強(qiáng)度與環(huán)境(2021年3期)2021-07-27
昆明醫(yī)科大學(xué)學(xué)報(bào)(2021年5期)2021-07-22
現(xiàn)代電子技術(shù)(2020年13期)2020-08-07
宇航計(jì)測(cè)技術(shù)(2018年3期)2018-09-08
高中時(shí)代(2017年7期)2018-02-24
