楊建平，帥曉勇，陶黃林

?

被淹沒地震信號(hào)的小波熵檢測(cè)與自動(dòng)識(shí)別方法

*楊建平1，帥曉勇1，陶黃林2

（1. 井岡山大學(xué)電子與信息工程學(xué)院，江西，吉安 343009；2. 井岡山大學(xué)數(shù)理學(xué)院，江西，吉安 343009）

為探測(cè)大震前的微震，保護(hù)大型煤礦、油田和礦山等重要設(shè)施，急需地震信號(hào)的實(shí)時(shí)處理、自動(dòng)識(shí)別和提取地震初至點(diǎn)等地震數(shù)據(jù)處理技術(shù)。采用了小波變換和信息熵理論相結(jié)合的一種具有多分辨率的復(fù)雜度參數(shù)——小波熵，該參數(shù)能夠從被淹沒環(huán)境中清晰地顯示出勘探數(shù)據(jù)中地震波到達(dá)所帶來的變化。結(jié)合實(shí)測(cè)數(shù)據(jù)進(jìn)行了仿真，并對(duì)比了單一的小波變換、數(shù)字帶通濾波器的監(jiān)測(cè)效果，結(jié)果表明小波熵參數(shù)能夠更好地自動(dòng)識(shí)別微震初至點(diǎn)。

小波變換；信息熵；FIR帶通濾波器；地震波；初至點(diǎn)

0 引言

通過地震觀測(cè)儀記錄、采集地震波并設(shè)法從波形中獲取地震的關(guān)鍵參數(shù)是地震勘探中一項(xiàng)最基本、也是最重要的內(nèi)容之一。分析和處理地震勘探數(shù)據(jù)，能夠及時(shí)檢測(cè)和識(shí)別大地震前的微震以及制定對(duì)策、及時(shí)保護(hù)大型煤礦、油田和礦山等重要設(shè)施帶來幫助。然而，通常在地震勘探數(shù)據(jù)的采集過程中，由于受到探測(cè)區(qū)內(nèi)以及周邊各種機(jī)械設(shè)備的影響，采集數(shù)據(jù)中不可避免地包含著各種噪聲的干擾，有時(shí)干擾的強(qiáng)度甚至?xí)艽螅ǖ卣鹦盘?hào)完全被淹沒），因此如何從地震勘探數(shù)據(jù)中自動(dòng)檢測(cè)出地震波初至點(diǎn)、發(fā)震時(shí)間段等參數(shù)將非常有利于地震勘探的研究，減少地震災(zāi)害帶來的影響。

通常，所勘探的地震波數(shù)據(jù)為非平穩(wěn)時(shí)間信號(hào)，變化強(qiáng)度、變化頻率都具有典型的非平穩(wěn)信號(hào)特征，因而采用傳統(tǒng)的、基于傅里葉分析的、具有全局性意義的功率譜方法無法表征非平穩(wěn)信號(hào)中最根本的特征——時(shí)頻的局域化特性。目前多用時(shí)變譜來刻畫時(shí)頻的非平穩(wěn)特性，如Iyama、Liu X以及Basu等人[1-6]分別利用在各尺度下能量無泄漏的正交小波基（或雙正交小波基）對(duì)地震波作小波變換，得到地震波的小波時(shí)變譜來探討具有非平穩(wěn)特性的地震波，分析小波變換系數(shù)與地震波演變之間的關(guān)系，取得了一定的研究成果，這些成果也表明小波變換能夠較好地處理地震波信號(hào)。另一方面，由于地震波數(shù)據(jù)是地震沿地球表層的復(fù)雜振動(dòng)過程中的一個(gè)觀測(cè)資料，是復(fù)雜地質(zhì)系統(tǒng)中的復(fù)雜現(xiàn)象，目前各種復(fù)雜現(xiàn)象、復(fù)雜系統(tǒng)的研究已成為非線性科學(xué)的前沿課題，如各種研究中經(jīng)常使用的關(guān)聯(lián)維、最大李雅普諾夫指數(shù)[7]、Lempel-Ziv方法的算法復(fù)雜度[8]、近似熵復(fù)雜度[9]等一系列描述非線性特征的參數(shù)來定量評(píng)價(jià)復(fù)雜系統(tǒng)的動(dòng)力學(xué)規(guī)律。本文欲用小波變換的時(shí)頻特性，結(jié)合描述系統(tǒng)復(fù)雜度的熵特征，首先提出并使用小波熵[10-11](wavelet entropy)分析方法對(duì)地震波的發(fā)生時(shí)間進(jìn)行提取，目前尚沒有發(fā)現(xiàn)其他文獻(xiàn)已采用該方法進(jìn)行地震淹沒信息的識(shí)別，該方法有利于計(jì)算機(jī)系統(tǒng)自動(dòng)識(shí)別、實(shí)時(shí)處理等操作的執(zhí)行。

1 地震勘探數(shù)據(jù)的小波分解思路

小波分析中對(duì)地震波信號(hào)的分解是基于小波變換的多分辨率分解算法[12]，算法采用正交小波基（或者是雙正交小波基）將地震波信號(hào)分解為多個(gè)不同頻帶上的分信號(hào)之和，這個(gè)分解過程本質(zhì)上是地震波信號(hào)的濾波過程，即將地震波信號(hào)重復(fù)通過一組高通濾波器和一組低通濾波器，逐個(gè)進(jìn)行濾波以達(dá)到將地震波信號(hào)逐級(jí)分解，由高通濾波器輸出地震波信號(hào)的高頻細(xì)節(jié)分量，低通濾波器輸出地震波信號(hào)的低頻粗略分量。經(jīng)過第一次濾波后得到的兩個(gè)地震波分信號(hào)所占頻帶寬度相等，即分別占原地震波信號(hào)總頻帶的二分之一，分解后地震波信號(hào)的采樣頻率降低了一倍；后面的每一步分解都是對(duì)其上一步的低頻段地震波信號(hào)重復(fù)使用一對(duì)高通和低通濾波器，進(jìn)而得到下一低頻層次上的兩個(gè)分解分量（低頻中的高頻和低頻中的低頻分量）。

2 勘探信號(hào)的小波熵

為了將地震勘探信號(hào)中地震時(shí)間段內(nèi)、外的小波系數(shù)的差異體現(xiàn)出來，而背景信號(hào)和地震波信號(hào)的各頻帶信息也不盡詳細(xì)，為表征這種差異熵?zé)o疑是一個(gè)最好的參數(shù)。信息熵是在一定的狀態(tài)下定位系統(tǒng)的一種信息測(cè)度，可以對(duì)勘探信號(hào)的未知程度進(jìn)行度量，可以用來估計(jì)這類隨機(jī)信號(hào)的復(fù)雜性。對(duì)于地震這個(gè)不確定性系統(tǒng)，若用一個(gè)取有限個(gè)值的隨機(jī)變量表示其狀態(tài)特征，取值為的概率為，且，則的某一結(jié)果所得到的信息用表示，的信息熵為：。

地震勘探信號(hào)小波熵具有鮮明的物理意義，是通過計(jì)算地震勘探信號(hào)的小波分解系數(shù)的信息熵所計(jì)算出來的一種復(fù)雜度，其值的大小與小波系數(shù)大小的分布特征有關(guān)，小波系數(shù)的大小又與地震勘探信號(hào)各頻率成分有關(guān)，因而其小波熵值的大小與地震波各頻率的分布有關(guān)，但與幅值的絕對(duì)大小無關(guān)，能反映地震勘探信號(hào)的不確定性（變化），小波熵大，則表明地震勘探信號(hào)中出現(xiàn)的變化多，地震勘探信號(hào)變化的速率快，即大的小波熵值對(duì)應(yīng)地震勘探信號(hào)變化是無序而復(fù)雜的，其復(fù)雜度也大；反之，復(fù)雜度越小，則說明發(fā)生變化的速率慢，變化是規(guī)則的，周期性較強(qiáng)。

3 仿真研究

圖1為中國數(shù)字地震臺(tái)網(wǎng)的長(zhǎng)春（ChangChun）地震信號(hào)波形圖，采樣間隔為0.02 s，采樣頻率為50 Hz，從圖中可以看出地震信號(hào)淹沒在振動(dòng)很強(qiáng)、長(zhǎng)周期的主干擾成分之中，另外還包含有隨機(jī)噪聲，很難識(shí)別地震波信號(hào)的初至點(diǎn)及其結(jié)束時(shí)刻。

圖1 淹沒在長(zhǎng)周期干擾中的地震信號(hào)

3.1 采用數(shù)字濾波器的波形識(shí)別

長(zhǎng)春地震信號(hào)主要受簡(jiǎn)單的單一長(zhǎng)周期成分干擾，其所含隨機(jī)噪聲與地震波相似（全部頻率成分的分布相對(duì)均衡），因而若能消除此長(zhǎng)周期干擾，則應(yīng)能更清晰地觀察出地震波的初至點(diǎn)。

這里設(shè)計(jì)一個(gè)有限長(zhǎng)脈沖響應(yīng)數(shù)字濾波器——FIR帶通濾波器，其阻帶的通帶邊界頻率（歸一化頻率）為wp=[0.032 0.2]，以濾除狹窄帶寬的具有單一頻率的周期成分（接近正弦信號(hào)的波形），獲取消除長(zhǎng)周期干擾成分的地震波信號(hào)，帶通濾波器的幅度響應(yīng)曲線和相位響應(yīng)曲線如圖2中的（a）、(b)所示。將選用的長(zhǎng)春臺(tái)地震勘探信號(hào)輸入此帶通濾波器以濾除長(zhǎng)周期干擾成分，獲得圖3所示的輸出波形（消除主干擾成分的地震勘探信號(hào)），觀察圖形可以粗略看出地震的初至點(diǎn)——圖中60 s附近波形幅度明顯變大的地方。但是，這種變化信息所反映的地震初至點(diǎn)、終止時(shí)刻，用人眼僅可粗略判別，若用自動(dòng)控制系統(tǒng)進(jìn)行自動(dòng)識(shí)別則非常困難，其原因有兩個(gè)：（1）如在濾波輸出波形的29 s附近處，隨機(jī)成分的振動(dòng)也比較強(qiáng)，自動(dòng)識(shí)別時(shí)很難設(shè)計(jì)出合適的算法將此種情況排除在外。（2）在60 s附近地震波形的隨機(jī)波動(dòng)性很強(qiáng)，借助計(jì)算機(jī)進(jìn)行自動(dòng)識(shí)別也將非常困難。

圖2 帶通濾波器的幅度和相頻響應(yīng)曲線：(a)幅度響應(yīng)曲線 (b)相位響應(yīng)曲線

圖3 帶通濾波器的輸出結(jié)果

3.2 地震信號(hào)的小波分解與識(shí)別

3.2.1 地震信號(hào)的小波系數(shù)

所勘探的地震波信號(hào)實(shí)際上可看作由兩部分疊加而成，一是背景信號(hào)（無地震波成分），另一是真正的地震波信號(hào)。對(duì)勘探信號(hào)進(jìn)行小波分解可分兩步來理解：背景信號(hào)的小波分解，其各頻帶成分體現(xiàn)了背景信號(hào)中各種頻率的分布特征；地震波信號(hào)的小波分解包含了地震波的各種頻率信息，因?yàn)榈卣鸩ㄓ傻竭_(dá)到結(jié)束只經(jīng)歷了采集信號(hào)的部分時(shí)間段，也就是地震波信號(hào)在震前、震后的振幅趨于0。對(duì)背景信號(hào)與地震波信號(hào)的合成信號(hào)（勘探信號(hào)）進(jìn)行小波分解，地震時(shí)間段上的小波分解系數(shù)（相比于地震時(shí)間段外）在各頻帶上變大。由于小波變換能突出信號(hào)的變化細(xì)節(jié)，這種變化比勘探信號(hào)中地震時(shí)間段內(nèi)原信號(hào)的變化更明顯，而地震時(shí)間段外的小波分解系數(shù)與背景信號(hào)的分解情況一致，因而可以利用地震時(shí)間段內(nèi)和地震時(shí)間段外的小波系數(shù)的差異來更有效地提取地震時(shí)間段的信息。

3.2.2 基于小波變換的地震信息檢測(cè)

對(duì)地震勘探信號(hào)進(jìn)行5層小波分解，圖4為低頻重構(gòu)信號(hào)及1～5層高頻重構(gòu)信號(hào)，低頻信號(hào)主要由長(zhǎng)周期干擾成分構(gòu)成，各高頻為噪聲及地震波的分解成分，若對(duì)不同頻帶取閾值濾除噪聲部分，如文獻(xiàn)[6]那樣可以達(dá)到很好的消噪效果。

圖4 地震勘探信號(hào)的5層小波分解，橫軸為時(shí)間，縱軸為幅度

小波變換將地震波分解成各頻帶的不同分量，根據(jù)小波變換的特性，地震波到達(dá)時(shí)間段上呈現(xiàn)更大的變化細(xì)節(jié)。從2～5層高頻重構(gòu)信號(hào)中均可觀察出地震波時(shí)刻在60 s附近，但在每一層信號(hào)變化的時(shí)間點(diǎn)上都有差異。這樣差異對(duì)于地震勘探儀來說無論采用哪一層高頻信號(hào)都很難借助自動(dòng)識(shí)別系統(tǒng)進(jìn)行自動(dòng)檢測(cè)并給出初至點(diǎn)及地震時(shí)間段，因?yàn)楦鞲哳l信號(hào)分量包含了地震波信號(hào)中的不同頻率成分，頻率變化信息所對(duì)應(yīng)的時(shí)間位置也各不相同。若用任一高頻段去分析地震波信號(hào)都將不夠全面，得不到反映地震情況的準(zhǔn)確信息，因而不能采用單一的小波變換方法進(jìn)行較為準(zhǔn)確的檢測(cè)、識(shí)別及提取相關(guān)參數(shù)。

3.3 小波熵算法自動(dòng)檢測(cè)地震波的初至點(diǎn)

由于采集地震勘探信號(hào)時(shí)已將周圍各種環(huán)境信息一并收集，這其中包含有周圍的環(huán)境信號(hào)，想直接從中檢測(cè)出地震的初至點(diǎn)尤其是當(dāng)?shù)卣鹦盘?hào)微弱時(shí)是非常困難的。當(dāng)環(huán)境信息相對(duì)穩(wěn)定時(shí)，隨著地震波到達(dá)使得整個(gè)勘探信號(hào)的成分增多，小波變換的高頻系數(shù)的幅度將增大，因而信號(hào)的小波熵復(fù)雜度變大；當(dāng)?shù)卣鸩ㄍＶ箷r(shí)，不含地震波的勘探信號(hào)小波熵復(fù)雜度將變小。因此，可以通過對(duì)勘探信號(hào)進(jìn)行小波變換求取各分量系數(shù)（或重構(gòu)分信號(hào)），再根據(jù)各分量系數(shù)并利用信息熵理論的計(jì)算公式求得小波熵值，該小波熵值的大小情況能夠反映地震波的到達(dá)信息。若采用一個(gè)移動(dòng)的窗口并依次計(jì)算各移動(dòng)窗口內(nèi)的小波熵，可獲得一小波熵的變化曲線，該曲線上不同點(diǎn)的小波熵值大小能夠反映不同窗口內(nèi)的勘探信號(hào)是否含有地震波。于是，為了探明地震波初至點(diǎn)、地震波時(shí)間段，可設(shè)計(jì)算法對(duì)小波熵曲線上的不同點(diǎn)進(jìn)行識(shí)別，以檢測(cè)出地震波的初至點(diǎn)等參數(shù)。通過上述的分析，設(shè)計(jì)識(shí)別算法步驟如下：

步驟1 采用db6小波基（選取原則是正交性好，信號(hào)分解后的細(xì)節(jié)更明顯）對(duì)地震波信號(hào)進(jìn)行三層小波分解，信號(hào)分解為低頻系數(shù)和一、二、三層高頻系數(shù)，進(jìn)行單支重構(gòu)得a3、d3、d2、d1等四個(gè)分信號(hào)，即；

步驟2 滑動(dòng)窗口法計(jì)算各窗口內(nèi)的小波熵：用移動(dòng)窗口逐步計(jì)算每一個(gè)窗口內(nèi)的小波熵，滑動(dòng)步長(zhǎng)為1，即用一個(gè)寬度為的時(shí)間窗從四個(gè)分信號(hào)的第1點(diǎn)開始截取長(zhǎng)為的一段，根據(jù)小波熵定義求出此時(shí)間窗的小波熵；

步驟3 將時(shí)間窗向右移動(dòng)一個(gè)點(diǎn)，得到四個(gè)分信號(hào)在第二窗口的小波熵；如此下去，可得到整個(gè)地震勘探數(shù)據(jù)各個(gè)時(shí)間點(diǎn)上的小波熵，即獲得勘探數(shù)據(jù)的小波熵曲線；

步驟4 計(jì)算小波熵曲線的差分獲取波形變化率的正、負(fù)值，正值為上升沿，負(fù)值對(duì)應(yīng)下降沿；以（最大值+平均值）/2為閾值搜索小波熵曲線中的閾值點(diǎn)，閾值點(diǎn)和上升沿重合時(shí)為地震波初至點(diǎn)，緊接著的閾值點(diǎn)和下降沿重合點(diǎn)為地震波終止點(diǎn)，并計(jì)算兩閾值點(diǎn)間的時(shí)間寬度；繼續(xù)按此方式搜索，尋找閾值點(diǎn)和上升沿（或下降沿），獲取下一個(gè)地震初至點(diǎn)，直到所有勘探信號(hào)結(jié)束。

圖5 地震勘探信號(hào)的小波熵曲線

圖5為地震勘探信號(hào)的小波熵曲線，其波形中小波熵“山峰”清晰可見，這段時(shí)間內(nèi)的小波熵大，包含有地震波到達(dá)所致的變化信息，即這個(gè)過程對(duì)應(yīng)地震波時(shí)間段。根據(jù)設(shè)定的閾值和上升沿（或下降沿）可檢測(cè)出地震到達(dá)時(shí)間點(diǎn)即初至點(diǎn)為56.82 s，終止點(diǎn)為68 s，地震時(shí)間為11.18 s。

4 討論與工作展望

本文雖然僅以淹沒在長(zhǎng)周期為主干擾成分的地震勘探信號(hào)為例，但文中的多分辨率小波熵算法能夠廣泛地適用于檢測(cè)地震波初至點(diǎn)的各種不同干擾情形，這個(gè)結(jié)論可以從小波熵處理信號(hào)的理論特點(diǎn)進(jìn)行推斷：只要地震勘探信號(hào)中的各種干擾成分在地震波到達(dá)前后性質(zhì)基本相同，則由于地震波到達(dá)時(shí)刻信號(hào)的成分增多、復(fù)雜度變大，此時(shí)的小波熵必定變大；而當(dāng)?shù)卣鸩ㄏr(shí)，信號(hào)的成分減少、復(fù)雜度變小，因而地震波時(shí)間段的復(fù)雜度大，其小波熵值也大，能夠從中識(shí)別和自動(dòng)提取地震波的初至點(diǎn)。該算法在實(shí)時(shí)監(jiān)控、自動(dòng)識(shí)別中能夠準(zhǔn)確地判斷微震初至點(diǎn)，但程序在計(jì)算中存在計(jì)算量較大的特點(diǎn)，今后的研究將致力于減少運(yùn)算量。

[1] Iyama J, Kuwamura H. Application of wavelets to analysis and simulation of earthquake motions [J]. Earthquake Engineering and Structural Dynamics, 1999, 28(3)：255-272.

[2] Liu X. Time-Arrival Location of Seismic P-Wave based on Wavelet Transform Modulus Maxima [J]. Journal of Multimedia, 2013, 8(1): 32-39.

[3] Basu B, Gupta V K．Seismic response of SDOF systems by wavelet modeling of non-stationary processes [J]. Journal of Engineering Mechanics, 1998, 124（10）： 1142-1150.

[4] Ao M S, Hu Y J, Zhao B, et al. Extraction of High-Rate GPS Seismic Wave Signals with Wavelet Packets Decomposition [J]. Earth Science/Diqiu Kexue, 2012, 37(5).

[5] Mukherjee S, Gupta V K. Wavelet-based characterization of design ground motions [J]. Earthquake Engineering and Structural Dynamics, 2002,31(5):1173-1190.

[6] En D, Wei N N, Wei H H, et al. Research of the Seismic Wave Data Processing Based on the Wavelet Transform[J]. Applied Mechanics and Materials, 2013, 273: 184-187.

[7] Krystal AD，Zaidman C，Greenside HS，et al．The largest Lyapunov exponent of the EEG during ECT seizures as a measure of ECT seizure adequacy[J]. Electroencephalography and Clinical Neurophysiology，1997,103 (6):599-606．

[8] Shaw FZ，Chen RF，Tsao HW，et al．Algorithmic complexity as an index of cortical function in awake and pentobarbital-anesthetized rats[J].Journal of neuroscience methods，1999，93(2)：101-110．

[9] Pincus S．Approximate entropy as an irregularity measure for financial data [J]．Econometric Reviews, 2008, 27(4)：329-362.

[10] 楊建平. 地磁場(chǎng)變化的小波熵復(fù)雜度分析方法[J]. 地球物理學(xué)進(jìn)展, 2010, 25(5): 1605-1611.

[11] 楊建平,朱平．分析Lorenz系統(tǒng)動(dòng)力學(xué)特征的新方法[J]. 計(jì)算機(jī)工程與應(yīng)用,2012,48(23):230-233．

[12] Mallat SG．A theory for multi-resolution signal decomposition: the wavelet representation [J]. IEEE Transaction on Pattern Analysis and Machine Intelligence, 1989,11:674-693.

METHOD OF DETECTION BY WAVELET ENTROPY AND IDENTIFICATION AUTOMATICALLY FOR SUBMERGED SEISMIC SIGNAL

*YANG Jian-ping1, SHUAI Xiao-yong1, TAO Huang-lin2

(1. School of Electronics and Information Engineering, Jinggangshan University, Ji’an, Jiangxi 343009, China 2. School of Mathematics and Physics, Jinggangshan University, Ji’an, Jiangxi 343009, China)

In order to detect the micro-seismic before large earthquake, protect the important facilities, such asthe large coal, oil, mine and so on. It’s an urgent need for seismic data processing technique, such as real-time process, recognize automatically and extract the submerged seismic onset point. A multi-resolution complexity parameter was acquired based on the wavelet transform and the theory of information entropy, the parameter can clearly shows the change in the exploration data from the arrivals of seismic waves. A simulation was done with the exploration data, Comparison of the monitoring effect of wavelet transform or digital band-pass filter, the results show that the parameter can be very good at the micro seismic onset point for automatic identification.

wavelet transform；information entropy；FIR band-pass filter；seismic wave；onset point

1674-8085(2015)04-0043-06

P318

A

10.3969/j.issn.1674-8085.2015.04.008

2015-03-30；修改日期：2015-06-22

國家自然科學(xué)基金項(xiàng)目（31260238）；江西省自然科學(xué)基金項(xiàng)目（20151BAB207063）

*楊建平(1970-)，男，江西吉安人，副教授，主要從事信號(hào)與信息處理研究(Email：yangjp9273@163.com);

帥曉勇(1977-)，男，江西吉安人，講師，主要從事通信信息處理(Email：jgsxythl@126.com);

陶黃林(1981-)，女，湖北黃岡人，講師，主要從事智能ERP、數(shù)據(jù)挖掘研究(Email：jgsxythl@126.com).