基于雙參數匹配追蹤算法的不同孔深爆炸地震波特性研究

周　輝， 龍　源, 鐘明壽， 謝興博， 郭　濤(解放軍理工大學 野戰工程學院，南京　210007)

?

基于雙參數匹配追蹤算法的不同孔深爆炸地震波特性研究

周輝， 龍源, 鐘明壽， 謝興博， 郭濤(解放軍理工大學 野戰工程學院，南京210007)

利用爆炸地震波信號的非平穩隨機特性，對其進行Hilbert變換轉換成復數信號，獲得爆炸地震波信號瞬時頻率和相位雙參數。使用雙參數匹配追蹤(Matching Pursuits)算法分解爆炸地震波信號，有效提高了匹配追蹤(MP)算法的掃描效率，降低算法復雜度，較傳統過完備匹配追蹤算法運算速率明顯提高;結合維格納瑞利分布(Wigner Ville Distribution,WVD)分析了不同孔深的爆炸地震波時頻特性，有效消除了WVD交叉干擾項的影響。分析結果表明：在一定深度范圍內，孔深增加，爆炸地震波能量增大，但高頻能量衰減更嚴重；當孔深超過一定深度時，爆炸地震波高頻能量衰減趨于穩定。

雙參數；匹配追蹤；維格納瑞利分布；孔深；爆炸地震波

爆炸地震波的傳播[1-2]受到多種因素的影響，由于介質條件的多變性和邊界條件的復雜性，目前尚不能從理論上對包含有爆炸荷載作用下特定介質本構關系的波動方程進行精確求解。因此，對爆炸地震波特性主要通過實驗和實際測試進行研究[3-5]。通過對實測數據進行歸整和合理剔除后，采用各種理論方法[6-7]在幅域，時域和頻域內分析處理數據，得到地震波的特性和基本規律。

MALLAT等[8]提出的匹配追蹤算法(Matching Pursuit,MP)是信號稀疏分解最常用的一類算法。其基本思想就是在滿足一定條件的情況下，用最少的一類信號原子來表示原信號。該方法有很多優良的性質，如較高的時間—頻率分辨率，暫態結構的局部自適應性，信號結構的參數表示更加靈活等[9-10]。其核心思想是將原信號表示為一系列與信號局部結構特征最佳匹配的時頻原子的線性組合，然后求各時頻原子維格納瑞利分布(Wigner Ville Distribution,WVD)并將其疊加，得到信號的時頻分布，并能很好的降低WVD交叉干擾項的影響[11]。由于子波受振幅、頻率、中心時間和相位控制，為了提高匹配追蹤算法的掃描效率，本文利用Hilbert變換，將地震信號轉換為復數地震信號[12-13]，提前提取地震信號瞬時頻率和相位參數信息來確定動態掃描范圍，提高了運算效率，取得了較好的效果。

本文擬通過在野外進行的多組不同孔深土中爆炸試驗測得的地震波信號，運用雙參數匹配追蹤算法來分析孔深對爆炸地震波特性的影響。

1　雙參數匹配追蹤算法的時頻分析

1.1匹配追蹤算法原理

匹配追蹤算法是一種貪婪算法，是將信號在超完備子波庫中進行分解，超完備意味著信號的分解目標在信號所組成的空間中足夠密集，因此無法保持相互正交的基函數特性，因此信號被分解為子波的線性組合[14]，假設D為進行信號分解的超完備子波庫，信號為f，長度為N，D中的元素滿足：

D={gγ:γ∈Γ}‖gγ‖=1

匹配追蹤算法通過把信號f垂直投影到子波庫D的匹配子波上。設gγ0∈D，則f可以表示為：

f=〈f,gγ0〉gγ0+Rf

(1)

式中：Rf為信號f利用匹配子波gγ0進行近似后的差值。為了使差值盡可能小， 就必須使內積項〈f，gγ0〉盡可能大。很顯然,gγ0與Rf是正交的，因此：

(2)

設R0f=f，且進行了n次迭代(n》0)得到差值Rnf，此時再選擇一個匹配子波gγn∈D，使其匹配Rnf，即

Rnf=〈Rnf,gγn〉gγn+Rn+1f

(3)

Rn+1f就是進行了n+1次迭代得到的差值。可見匹配追蹤就是利用式(3)描述的一個重復迭代過程，若迭代m次，則可將f表示為如下形式：

(4)

因此經過m次分解計算后，原始信號可表示為m個匹配子波的線性組合，其誤差為第m-1次迭代計算后的差值。雖然MP算法是屬于非線性迭代過程，但能量有限的特點可以保證它的收斂性。

由于Gaussian 函數具有良好的時頻聚集性，一般選用Gabor原子構建原子庫。先將此基函數的時頻參數進行離散化，形成過完備的原子庫，原子庫中Gabor原子可表示為：

(5)

在通過將待分解的信號逐次從原子庫中進行比對迭代，直至差值滿足一定的條件。該方法的缺點是計算量大，運算速度慢。

1.2雙參數匹配追蹤算法

匹配追蹤的最佳基函數求解實質涉及到多極值優化求解，傳統的匹配追蹤算法采用的是貪心算法，雖能獲得全局最優解但存在計算量龐大的不足。在保證對信號分解精度的基礎上，不少人對如何提高匹配追蹤運算效率進行了大量研究[14-16]，考慮到Gabor原子受振幅、頻率、中心時間和相位四個參數控制，對過完備原子庫的掃描實質上是對四個參數的整體尋優過程，若能根據實際情況提前獲取某些參數的取值范圍，便能降低基函數的掃描范圍，提高算法運算速率。

利用Hilbert變換可提前提取其瞬時頻率和相位，但Hilbert變換要求原始信號為單頻率信號，爆炸地震波信號明顯不符合這一條件。由于可將爆炸地震波信號可看成由一定數量的子波組成，這也是進行匹配追蹤的前提。再此基礎上，可對爆炸地震波信號進行Hilbert變換，獲取優勢頻率范圍，再建立動態子波庫[17]，再基于此在有限的動態子波庫范圍內進行匹配搜索，可以有效降低運算量，提高運算速率。

具體算法步驟如下：①離散化Gabor基函數的時頻參數，形成原子庫Di(i=1，2，…，I)；②將原始信號X(t)賦值給初始差值信號r0；③用Hilbert 變換計算差值信號rm(m=0，1，2，…，M-1)；④計算信號的瞬時包絡、瞬時相位和瞬時頻率；⑤找到包絡最大值及其對應的時間位置，計算相應時間處的瞬時頻率和瞬時相位，得到Gabor函數相位角和主頻的估計值；⑥將差值rm從限定的原子庫中找出與原始信號最匹配的原子dmi，求出匹配細數cmi，并將差值減去匹配原子，得到新差值信號rm+1；⑦重復步驟③～⑥步，直到差值小于一定的閾值，完成了對信號X(t)的分解：

(6)

1.3算法復雜度比較

假定對一長度為N的信號進行匹配追蹤分解，為了保證原子庫的過完備性，需對四個參數進行過完備掃描，對中心時間u和頻率v掃描的計算復雜度均為O(N)，振幅掃描計算復雜度為O(log2(N))，相位參數在0～2π之間掃描，故傳統過完備匹配追蹤算法復雜度為(N2.log2(N))。而雙參數匹配追蹤算法已獲取頻率相位參數信息，其算法復雜度為O(N.log2(N))，運算速率明顯提高，且分解信號長度越長，雙參數匹配追蹤算法運算速率優勢越大，更適合處理數據量大的地震波信號。

1.4基于匹配追蹤算法的時頻分析

時頻分析是地震信號數據處理中一種非常重要的信號處理方法，把地震數據變換到頻率域是一系列重要地震資料處理算法和解釋技術的基礎。WVD是由瑞利(Ville)引進信號分析的并對魏格納分布給出了一個似乎合理的論證，并根據特征函數方法推導得出[18]。

信號X(t)的WVD如下：

(7)

由于式中不含有任何的窗函數，因此避免短時傅里葉變換時間分辨率與頻率分辨率相互牽制的矛盾，它的時間-帶寬達到了測不準原理給出的下界。但是維格納-瑞利分布本質不是線性的，即兩信號和的 WVD并不等于每一個信號的 WVD之和。令X(t)=X1(t)+X2(t)，則

Wx1(t,ω)+Wx2(t,ω)+2Re[Wx1+x2(t,ω)]

(8)

式中:2Re[Wx1+x2(t,ω)]是X1(t)和X2(t)互WVD，稱之為交叉項。

由式(8)可以看到：有時魏格納分布在時間和頻率上把這些值置于兩個信號的中間；有時這些值又處在時-頻平面和所期望的東西爭奪位置。因此產生了交叉項。交叉項極大的干擾時頻分布，同時也抑制了二次型時頻分布的推廣。若能結合匹配追蹤算法，將信號分解成基本原子的線性疊加，分別對每個信號做WVD，在將其疊加，能達到很好的時頻分析效果。例如，圖1是由四種頻率成分的調制信號。

圖1　四頻率成分信號Fig.1　Signals with four frequencies

分別直接計算WVD和對其進行子波分解后計算WVD，結果如圖2所示。

圖2　Wigner-Ville 時頻圖Fig.2　Time-frequency of Wigner-Ville distribution

在圖2(a)中，處理后的信號細節部分主要集中在三個時間段和三個頻率段上，而原始四頻率成分信號主要是由兩兩頻率或相位相近的信號組成，這表明處理后的信號并不能很好的還原原始合成信號的細節部分。這是因為由式(8)可以看到WVD是在時間域和頻率域將這些值置于兩個信號之間，產生了干擾項。四頻率成分信號兩兩相互形成干擾項，共產生 6 個交叉項(其中中間陰影部分的交叉項為重疊交叉項，為由左上與右下交叉項和左下與右上交叉項的重疊)。

若將WVD對t進行積分，可以得到：

(9)

式(9)表明，若將WVD在某頻率處對t積分等于該頻率處的能量譜。

2　不同孔深土中爆炸試驗

2.1試驗方法

本次試驗的方法是進行多組不同孔深的爆炸試驗，裝藥選用的是75 gTNT藥柱，采用電起爆法。分組進行試驗，孔深依次為0 m、0.2 m、0.4 m、0.5 m、0.6 m、0.8 m。炸藥豎直放置在孔底部，不回填。試驗方案(見圖3)。

圖3　土中爆炸試驗方案Fig.3　Testing program of the explosion in earth

分別在距離炮孔3 m、4 m、5 m處布置傳感器。傳感器采用中國地震局工程力學研究所研制的891-Ⅱ型磁電式拾振傳感器，主要用于感應地震波、機械振動和各種沖擊信號。爆破測震儀采用成都中科測控有限公司生產的TC-4850型爆破測震儀，主要用于對地震波、機械振動和各種沖擊信號進行記錄。儀器為自適應量程，一般為35 cm/s。試驗現場布置(見圖4)，設置儀器采樣時間為10 s，采樣率8 kHz。

各組試驗具體的爆破參數見表1

圖4　現場試驗布置Fig.4　Field test arrangment

表1　孔深試驗爆破參數

2.2試驗結果

根據試驗所測地震波數據得到各測點不同孔深土中爆炸地震波質點振動速度峰值主頻情況見圖5。

由圖5(a)可知隨著孔深增加，地震波質點振動速度增大，表明炸藥的能量更多的作用于土壤內部，較少的耗散到空氣中，使得地震波的能量增大；當在地面爆炸時，產生的地震波明顯較土中爆炸產生的地震波弱。

圖5　各測點處振動峰值與主頻Fig.5　The vibration peak and main frequency of each measuring point

由圖5(b)可知，地面爆炸產生的地震波主頻較高，而土中爆炸產生的爆炸地震波主頻低，為幾十赫茲(Hz)。表明高頻分量的地震波在土中衰減較快。

3　不同孔深爆炸地震波特性分析

3.1地震波信號分析

利用雙參數匹配追蹤算法分解爆炸地震波，以孔深60 cm，5 m處測點為列，分解結果見圖6。

圖6　爆炸地震波信號分解Fig.6　Decomposition of explosion seismic wave

圖6為在使得殘差信號比原始信號低一個數量級的標準下進行的分解，因其算法復雜度為O(N.log2N))，較過完備匹配追蹤少一個信號長度，使得其運算速率明顯提高。選定的殘差信號標準保證了分解結果的有效性和準確性，并達到了去噪的效果，更能反應實際地震波特性。對其進行時頻分析，得到其WVD和頻譜能量圖見圖7。

圖7　爆炸地震波時頻分析Fig.7　Time-frequency analysis on explosion seismic wave

由圖可知爆炸地震波能量譜有較高的分辨率，地震波能量集中在低頻部分，且主頻處的能量并不一定為爆炸地震波頻譜能量的峰值。

3.2不同孔深地震波特性分析

分別計算爆炸地震波頻譜能量，得到最大能量及其對應的頻率如圖8所示。

圖8　各測點頻譜峰值能量及應頻率Fig.8　The peak energy and frequency of each test points

由圖8可知，在一定深度范圍內，孔深增加，爆炸地震波峰值能量增加，更多的向低頻分量集中，高頻分量衰減更嚴重，但當孔深超過一定深度時，峰值能量和頻率變化不大，表明此時孔深的影響較小。因此，淺孔爆破時，孔深對土中爆炸地震波頻譜能量有較大影響，在進行爆炸地震波危害防護時應予以考慮。

為了進一步量化分析不同孔深爆炸地震波特性，分別計算不同孔深爆炸地震波峰值能量隨傳播距離的衰減斜率，計算結果如圖9所示。

圖9　不同孔深處爆炸地震波峰值能量衰減斜率Fig.9　The slope of peak energy attenuation with different hole depth

計算結果表明，在一定深度范圍內，孔深增加，爆炸地震波的能量衰減更嚴重，但當深度超過一定范圍時，爆炸地震波峰值能量衰減不再增加。

4　結　論

(1) 雙參數匹配追蹤算法能夠極大的提高分解信號的運算速率，符合處理大量地震信號數據的需求；

(2) 利用匹配追蹤分解地震信號后進行WVD計算，能夠充分利用其沒有窗函數的高分辨率的優勢并能很好消除其交叉干擾項的影響；

(3) 在一定范圍內，孔深增加，爆炸地震波峰值能量增加，高頻衰減更嚴重，當孔深超過一定深度時，爆炸地震波峰值能量衰減不再增加；

(4) 淺孔土中爆破作業時，在進行爆炸地震波危害防護時應予以考慮合理的孔深參數。

[1] 張義平，吳桂義. 爆破地震波特性研究[J]. 礦業研究與開發，2007，27(6)：68-72.

ZHANG Yiping, WU Guiyi. Study on characteristics of blasting-caused seismic wave[J]. Mining R & D,2007,27(6):68-72.

[2] 龍源，馮長根，徐全軍，等. 爆破地震波在巖石介質中傳播特性與數值計算研究[J]. 工程爆破，2000，6(3)：1-7.

LONG Yuan, FENG Changgen, XU Quanjun, et al. Study on profagation characteristics of blasting seismic waves in a rock medium and nuerical calculation[J]. Engineering Blasting, 2000, 6(3): 1-7.

[3] 施鵬，鄧國強，楊秀敏，等. 土中爆炸地沖擊能量分布研究[J]. 爆炸與沖擊，2006，26(3)：240-244.

SHI Peng, DENG Guoqiang, YANG Xiumin, et al. Study on ground shock energy distribution of explosion in soil[J]. Explosion and Shock Waves, 2006, 26(3): 240-244.

[4] 趙躍堂，鄭大亮，王明洋，等. 飽和土介質中爆炸波傳播的時間特征試驗研究[J]. 防災減災工程學報，2004，24(2)：162-167.

ZHAO Yuetang, ZHENG Daliang, WANG Mingyang, et al. Experimental investigation on time characteristics of explosive wave propagation[J]. Journal of Disaster Prevention and Mitigation Engineering, 2004, 24(2): 162-167.

[5] 施鵬，辛凱，楊秀敏，等. 土中裝藥不同埋深爆炸試驗研究[J]. 工程力學，2006，23(12)：171-174.

SHI Peng, XIN Kai, YANG Xiumin, et al. Experimental study of explosion with different burial depths soil[J]. Engineering Mechanics, 2006, 23(12): 171-174.

[6] 李興華，龍源，紀沖，等. 基于小波包變換的高程差對爆破震動信號影響分析[J]. 振動與沖擊，2013，32(4)：44-47.

LI Xinghua, LONG Yuan, JI Chong, et al. Influence of height difference based on wavelet packets transformation on blasting vibration signals[J]. Journal of Vibration and Shock, 2013, 32(4): 44-47.

[7] XING Xiaojun, HE Xilei, PU Yong, et al. High-precision frequency attenuation analysis and its application[J]. Applied Geophysics, 2011, 8(4): 337-343.

[8] MALLAT S, ZHANG Z. Matching pursuit with time-frequency dictionaries[J]. IEEE Trans Signal Processing，1993,41(12): 3397-3415.

[9] WANG Y H. Seismic time-frequency spectral decomposition by matching pursuit[J]. Geophysics,2007,72(1):13-20.

[10] 劉繼承，富爽. 基于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.

[11] 趙中華，王文延. 基于時頻分析的去除魏格納交叉干擾項的方法[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.

[12] 黃玱，金河. 地震資料復數道分析的解釋基礎和應用[J]. 石油物探，1984，22(3)：1-19.

HUANG Qiang, JIN He. The oretical basis for interpretation of the seismic data by complex trace analysis and its application[J]. Geophysical Prospecting for Petroleum, 1984, 22(3): 1-19.

[13] 武國寧，曹思遠，孫娜. 基于復數道地震記錄的匹配追蹤算法及其在儲層預測中的應用[J]. 地球物理學報，2012，55(6)：2027-2034.

WU Guoning, CAO Siyuan, SUN Na. Matching pursuit method based on complex seismic traces and its application of hydrocarbon exploration[J]. Chinese Journal of Geophysics, 2012, 55(6): 2027-2034.

[14] 邵君. 基于MP的信號稀疏分解算法研究[D]. 北京：中國科學院研究生學院(軟件研究所)，2002.

[15] 熊純，何衛鋒. 基于匹配追蹤算法的自適應時頻分解改進及應用[J]. 長沙航空職業技術學院學報，2007，7(4)：44-48.

XIONG Chun, HE Weifeng. Application of adaptive time-frequency decomposition with modified matching pursuit [J]. Changsha Aeronautical Vocational and Technical College Journal, 2007, 7(4): 44-48.

[16] 孟慶豐. 信號特征表示方法與應用[D]. 西安：西安電子科技大學，2006.

[17] 張繁昌，李傳輝，印興耀. 基于動態匹配子波庫的地震數據快速匹配追蹤[J]. 石油地球物理勘探，2010，45(5)：667-673.

ZHANG Fanchang, LI Chuanhui, YIN Xingyao. The fast matching pursuit of seismic data based on dynamic matching library [J]. Oil Geophysical Prospecting, 2010, 45(5): 667-673.

[18] 王培茂. 地震信號時頻特征表示方法[D]. 長春：吉林大學， 2008.

Characteristics analysis of explosion seismic waves with different hole depth based on the method of double parameters matching pursuit

ZHOU Hui, LONG Yuan, ZHONG Mingshou, XIE Xingbo, GUO Tao

(College of Field Engineering, PLA University of Science & Technology, Nanjing 210007, China)

Taking advantage of the non-stationary random nature of explosion seismic signals, and using the Hilbert transformation to covert the signal into a plurality signal, the two parameters, namely the instantaneous frequency and phase parameter of the signal were acquired. The method of double parameter matching pursuit was introduced to decompose explosion seismic wave signals. The scan efficiency of the matching pursuit (MP) algorithm, the complexity of the algorithm, and the operation rate of the algorithm were effectively improved compared with the ordinary complete matching pursuit algorithm. With resorting to the Wigner Ville distribution, in which the influence of its cross interference terms was effectively eliminated, the time-frequency characteristics of the explosion seismic waves with different hole depth were analysed. The results show that: in a certain depth range, with the increase of hole depth, the energy of the seismic wave is increased, while the attenuation of the explosion seismic wave energy of high frequency also becomes more intense. When the hole depth is over a certain value, the attenuation of the explosion seismic wave energy of high frequency remains at a stable level.

double parameter; matching pursuit; WVD; hole depth; explosion seismic wave

國家自然科學基金項目(51304218;51178460)

2015-10-09修改稿收到日期：2016-03-16

周輝 男，博士生，1990年生

龍源 男，博士，教授，1958年生

E-mail:long_yuan@sohu.com

TB236

A DOI:10.13465/j.cnki.jvs.2016.14.013