小波包變換中地震信號的結點序號到頻帶序號的轉換算法

白 泉, 韓晶晶, 康玉梅, 邊晶梅

(1.沈陽工業大學建筑與土木工程學院,遼寧 沈陽 110870; 2.賀州學院建筑工程學院,廣西 賀州 542899;3.東北大學資源與土木工程學院,遼寧 沈陽 10004; 4.沈陽化工大學經濟與管理學院,遼寧 沈陽 110142)

技術交流

小波包變換中地震信號的結點序號到頻帶序號的轉換算法

白 泉1, 韓晶晶2, 康玉梅3, 邊晶梅4

(1.沈陽工業大學建筑與土木工程學院,遼寧 沈陽 110870; 2.賀州學院建筑工程學院,廣西 賀州 542899;3.東北大學資源與土木工程學院,遼寧 沈陽 10004; 4.沈陽化工大學經濟與管理學院,遼寧 沈陽 110142)

利用小波包變換對地震非平穩信號進行處理時,結點與頻帶之間存在“跳頻”現象,基于小波包變換算法以及異或運算,對小波包樹頻帶序號和結點序號的關系進行研究,明確信號子空間頻帶與小波包樹結點的排列規律,同時提出一種從結點序號到頻帶序號的轉換算法。基于MATLAB平臺,以唐山(南北向)波為例,驗證該算法的正確性。

小波包; 樹結點; 頻帶排列; 異或運算

0 引言

地震動分析是研究地震危害控制的基礎,也是控制地震危害的前提。小波包變換可將地震信號均勻分解到不同的頻帶上,從而能夠準確分析各個頻帶上地震信號的詳細特征。但在MATLAB平臺中,由于算法的限制,小波包變換后得到的頻帶順序與小波包樹結點順序并不一致,反而會出現“跳頻”現象[1-4]。若忽視這一問題,就不能準確定位頻帶,

進而影響后續時頻分析的正確性[5],因此明確頻帶與小波包樹結點的對應關系對時頻分析意義重大。

本文基于MATLAB平臺,借助小波包算法和異或運算定義,主要對小波包變換后的頻帶序號和樹結點序號的關系進行研究,指出信號子空間頻帶與小波包分解樹結點的排列規律,提出一種從結點序號到頻帶序號的轉換算法,并舉例驗證該轉換算法的正確性,這對提取地震信號局部特征、有層次地展現地震信號等具有重要意義。

1 小波包理論

(1)

利用濾波器組實現小波包變換的關鍵是隔點下采樣,采樣后分信號與原信號的關系為[2]:

(2)

X(N/2+f)=∑x(2n)e-j(N/2+f)·2π·2n/N+ ∑x(2n+1)e-j(N/2+f)·2π·(2n+1)/N

(3)

式中:X(f)表示x(n)的離散傅里葉變換;N為采樣點數。由式(2)和(3)可得:

X′(f)=0.5[X(f)+X(N/2+f)]

(4)

式中:X′(f)為下采樣后得到的序列的傅里葉變換。由離散傅里葉變換的對稱性可知:

X(N/2+f)=X(N/2-f)

(5)

將式(5)代入式(4)中可得:

X′(f)=0.5[X(f)+X(N/2-f)]

(6)

根據采樣定理,當X(f)的采樣頻率為f0時,其有效頻段為(0,f0/2),故X′(f)的采樣頻率為f0/2,有效頻段為(0,f0/4),且有:

(7)

由此可以清楚地分析小波包分解信號的過程:當信號采樣頻率為2f時,其有效頻段為(0,f),在小波包變換中,經過低通和高通濾波器進行第一層分解后,整個頻帶被劃分為低頻頻帶L(0,f/2)和高頻頻帶H(f/2,f);對這兩個頻帶的信號進行下采樣后再分別進行第二次分解,每個頻帶又被劃分為兩個頻帶,其中低頻頻帶劃分得到(0,f/4)和(f/4,f/2)頻帶,高頻頻帶劃分得到(f/2,3f/4)和(3f/4,f)頻帶;如此不斷重復低頻頻帶和高頻頻帶的二進劃分,就可以分離出不同頻帶,從而實現了對頻帶的均勻細分[1,3,7]。

小波包分解過程可以用分解樹來描述,如圖1所示。圖中Uj,n表示小波包分解的子空間編號,亦可稱作結點編號;j表示分解層數,其對應的分解尺度為2j,第j層一共有2j個正交基;n表示第j分解層的子空間編號,第j層的第n個子空間由集合{uj,n(t-k)}構成。

由小波包變換理論知識可知,同一層的各個子空間之間相互正交且每一層上所有子空間的直和構成這個信號空間。

2 地震信號頻帶劃分參數

對地震動非平穩信號進行小波包變換時,主要參數有分解層數和分解路徑。

圖1 小波包分解Fig.1 Wavelet packet decomposition

2.1 分解層數

在進行小波包分解時,除了確定小波基函數,還需要確定分解層數。這里選用字母j表示小波包分解層數,通常分解層數按下式選取[2,7]:

(8)

式中:Ls為輸入信號的長度。對于地震非平穩信號,采樣時間間隔一般為0.01 s或0.02 s,即采樣頻率為100 Hz或50 Hz。地震動信號一般持續時間為10～20 s,若采樣頻率為50 Hz,則信號長度約為29～210。綜合考慮頻帶寬度及分辨率,本文認為分解層數取7一般能滿足地震信號的分析要求。

2.2 分解路徑

如上所述的小波包分解樹中,樹結點從左至右依次編號為0,1,2,…2j-1,稱作結點序號。地震信號的小波包分解實際上是把信號通過一組濾波器組合(低通或高通)后,到達底端的樹結點。設定0表示低通濾波,1表示高通濾波,則各濾波器組合可由一個二進制數表示[8],稱該二進制數為實際路徑。例如,實際路徑101表示信號經過一次高通濾波器,一次低通濾波器,最后又經過一次高通濾波器,如圖2中虛線所示。該二進制數對應的十進制數就是分解終端樹結點的頻帶序號,因此實際路徑又稱頻帶路徑。與頻帶路徑相對應,把結點序號的二進制數稱作結點路徑或理想路徑。

圖2 101分解路徑Fig.2 Decomposition path 101

除了頻帶路徑和結點路徑,每個結點還對應一個通頻帶。通頻帶是指通過該結點的頻帶范圍,與頻帶序號有關。若信號的有效頻帶為f,分解層數為j,則頻帶序號m對應的通頻帶為m*(f/2j)-(m+1)*(f/2j),其中m=0,1,2,…,2j-1。

需要說明的是,分解層數和分解路徑之間存在的關系:分解路徑的每個二進制數的位數即為獲得該頻帶所需要的分解層數,這樣可以準確獲取各個頻帶分解路徑,尤其是以0開頭的分解路徑,同時為后續準確得到頻帶序號及對應的通頻帶做鋪墊。

3 頻帶序號與結點序號的變換

3.1 小波包分解的頻帶排列規律

地震信號經過小波包分解后,可看作由每層所有子空間組合構成。小波包分解采用隔點采樣算法,得到的頻帶信號分解分量按照自然序號排列,而不是頻率由小到大排列,因此各個頻帶信號之間存在排序紊亂現象。若想對各頻帶分信號進行準確分析,就需要找到其實際的頻帶序號排列規則。

通過對多個信號進行不同層次的小波包分解總結,得到的規律如下:由高通濾波器得到的頻段,進行下一階段濾波時高低頻段要發生位置互換。這有助于準確定位所需頻段對應的小波包樹結點,將其用小波包樹的形式表示出來,如圖3所示。

圖3 對應頻帶的小波包分解樹Fig.3 Wavelet packet decomposition tree of corresponding frequency band

由圖3可知,小波包空間分解的子頻帶按照結點序號順序排列時,并不是完全按照頻率遞增順序排列的。比如,對于分解的第三層,結點序號依次為0,1,2,3,4,5,6,7,是自然順序,而頻帶序號依次為0,1,3,2,7,6,4,5。這種排列順序驗證了上述規律的正確性,即由高通濾波器得到的頻段進行下一階段濾波時,高低頻段要發生位置互換。

3.2 結點序號到頻帶序號的轉換關系

當分解層數較少時,可以根據上述規律手動依次得到每個結點的頻帶路徑和通頻帶,進而得到完整的小波包樹。但是當分解層數較多時,手動獲取比較困難,因此需要找到結點序號到頻帶序號之間的換算關系。

按照圖3的3層小波包分解樹可以得到每個結點的結點分解路徑和頻帶分解路徑。比如,第3層小波包分解結點序號為7的結點,其結點路徑為111,對應的頻帶濾波路徑為HLH,可得到頻帶路徑為101,對應的頻帶序號為5,通頻帶為(5*f/23,6*f/23)。通過多個結點的分析,得到結點序號到頻帶序號的轉換過程:

對結點序號進行二進制轉換得到結點分解路徑(一個二進制數),從左至右對該路徑的每位二進制數與其左側所有位數進行異或運算,得到頻帶分解路徑(一個新的二進制數),對新的二進制數進行十進制轉換,即可得到頻帶序號。比如,經過3層小波包分解得到的小波包樹,結點序號為a,a的結點分解路徑的二進制數為ABC,異或運算后的頻帶分解路徑二進制數為XYZ,則X=A、Y=A?B、Z=A?B?C,對二進制數XYZ轉化為十進制數b,則b即為結點a對應的頻帶序號。

上述的“異或”是一種數學邏輯運算符,數學符號為“?”,計算機符號為“XOR”,其運算法則為:

(9)

根據上述換算過程和運算法則,在MATLAB平臺上編制程序,可以實現小波包分解到任意分解層數時,地震信號的結點序號與頻帶序號的轉換。表1為7層小波包變換下的結點序號與經過上述轉換得到的頻帶序號。

4 實例驗證

選取唐山波南北向(北京觀測)的強震記錄進行實例分析,驗證本文所述地震信號頻帶特性分析中頻帶序號與結點序號之間轉換算法的正確性。地震信號采樣點數為2 000,采樣的時間間隔為0.01 s,采用sym5小波函數,在MATLAB平臺上編制程序對唐山波進行7層小波包變換,得到各樹結點的頻譜圖。本例中,采樣頻率為100 Hz,根據采樣定理,有效頻率為0～50 Hz,進行7層分解后每個頻帶寬度為f/27,即50/27Hz。限于篇幅,給出結點序號分別為11、12、34、36、53、55的頻譜圖,如圖4所示。

表1 結點序號與頻帶序號

續表1

結點序號1617181920212223結點路徑00100000010001001001000100110010100001010100101100010111頻帶路徑00111110011110001110000111010011000001100100110110011010頻帶序號3130282924252726結點序號2425262728293031結點路徑00110000011001001101000110110011100001110100111100011111頻帶路徑00100000010001001001100100100010111001011000101000010101頻帶序號1617191823222021結點序號3233343536373839結點路徑01000000100001010001001000110100100010010101001100100111頻帶路徑01111110111110011110001111010111000011100101110110111010頻帶序號6362606156575958結點序號4041424344454647結點路徑01010000101001010101001010110101100010110101011100101111頻帶路徑01100000110001011001101100100110111011011001101000110101頻帶序號4849515055545253結點序號4849505152535455結點路徑01100000110001011001001100110110100011010101101100110111頻帶路徑01000000100001010001101000100100111010011001001000100101頻帶序號3233353439383637結點序號5657585960616263結點路徑01110000111001011101001110110111100011110101111100111111頻帶路徑01011111011100010110001011010101000010100101010110101010頻帶序號4746444540414342結點序號6465666768697071結點路徑10000001000001100001010000111000100100010110001101000111頻帶路徑11111111111110111110011111011111000111100111110111111010頻帶序號127126124125120121123122結點序號7273747576777879結點路徑10010001001001100101010010111001100100110110011101001111頻帶路徑11100001110001111001111100101110111111011011101001110101頻帶序號112113115114119118116117結點序號8081828384858687結點路徑10100001010001101001010100111010100101010110101101010111頻帶路徑11000001100001110001111000101100111110011011001001100101頻帶序號96979998103102100101結點序號8889909192939495結點路徑10110001011001101101010110111011100101110110111101011111頻帶路徑11011111101110110110011011011101000110100111010111101010頻帶序號111110108109104105107106結點序號96979899100101102103結點路徑11000001100001110001011000111100100110010111001101100111頻帶路徑10000001000001100001110000101000111100011010001001000101頻帶序號6465676671706869結點序號104105106107108109110111結點路徑11010001101001110101011010111101100110110111011101101111頻帶路徑10011111001110100110010011011001000100100110010111001010頻帶序號7978767772737574結點序號112113114115116117118119結點路徑11100001110001111001011100111110100111010111101101110111頻帶路徑10111111011110101110010111011011000101100110110111011010頻帶序號9594929388899190結點序號120121122123124125126127結點路徑11110001111001111101011110111111100111110111111101111111頻帶路徑10100001010001101001110100101010111101011010101001010101頻帶序號8081838287868485

按照本文給出的變換算法,參照表1的頻帶序號與結點序號的對應關系,取圖4中的各點如表2所列。

圖4 樹結點頻譜Fig.4 Spectra of tree nodes

結點序號頻帶序號通頻帶結點序號頻帶序號通頻帶11135．078～5．468365621．875～22．2661283．125～3．516533814．844～15．234346023．437～23．828553714．453-14．844注：頻帶序號m對應的通頻帶為m?(f/2j)-(m+1)?(f/2j)，單位為Hz。

5 結論

小波包變換可將地震信號均勻分解到不同的頻帶上,但由于算法限制,變換后得到的樹結點順序與頻帶順序并不一致,即“跳頻”。本文針對這一問題,基于MATLAB平臺對小波包變換中小波包分解樹結點序號與頻帶序號所存在的錯亂現象進行了簡單介紹,并從小波包算法和異或運算的角度,通過濾波器組的基本原理來實現小波包變換的算法,明確并分析了頻帶均勻劃分理論。借助分解路徑,總結了小波包變換的頻帶順序排列規律,并由此得到了小波包分解樹的結點序號轉換到頻帶序號的算法,為確定頻帶和小波包分解樹各結點的對應關系提供了一種有效的方法。若不考慮能量泄漏,圖中給出的各結點通頻帶范圍與本文所提供的轉換算法結果基本一致,由此驗證了本文所述地震信號從結點序號轉換到頻帶序號的算法的正確性,對有層次地展現地震信號等工程應用具有重要意義。

References)

[1] 曾宇清,王衛東,賀啟庸.按頻帶順序排列的小波包新算法及應用[J].力學學報,1998,30(2):186-192. ZENG Yu-qing,WANG Wei-dong,HE Qi-yong.Theory and Application of New Wavelet Packets Algorithm with Results in Order of Frequency-bands[J].Acta Mechanica Sinica,1998,30(2):58-64.(in Chinese)

[2] 薛蕙,楊仁剛,郭永芳.小波包變換(WPT)頻帶劃分特性的分析[J].電力系統及其自動化學報,2003,15(2):5-8. XUE Hui,YANG Ren-gang,GUO Yong-fang.Frequency Division Character of Wavelet Packet Transform[J].Proceedings of the EPSA,2003,15(2):5-8.(in Chinese)

[3] 曾憲偉,趙衛明,盛菊琴.小波包分解樹結點與信號子空間頻帶的對應關系及其應用[J].地震學報,2008,30(1):90-96. ZENG Xian-wei,ZHAO Wei-ming,SHENG Ju-qin.Corresponding Relationships between Nodes of Decomposition Tree of Wavelet Packet and Frequency Bands of Signal Subspace[J].Acta Seismologica Sinica,2008,30(1):90-96.(in Chinese)

[4] V L Pham,K P Wong.Antidistortion Method for Wavelet Transform Filter Banks and Nonstationary Power System Waveform Harmonic Analysis[J].IET Proceedings-Generation,Transmission and Distribution,2001,148(2):117-122.

[5] 許康生,李英,李秋紅.近地震波的小波相對能量分布特征分析[J].地震工程學報,2013,35(1):166-170. XU Kang-sheng,LI Ying,LI Qiu-hong.Distribution Characteristics of Wavelet Relative Energy on Near-earthquake Wave[J].China Earthquake Engineering Journal,2013,35(1):166-170.(in Chinese)

[6] 彭玉華.小波變換與工程應用[M].北京:科學出版社,1999:21-69. PENG Yu-hua.Wavelet Transform and Engineering Application[M].Beijing:Science Press,1999:21-69.(in Chinese)

[7] 楊世錫,項文娟,葉紅仙.基于Lab-VIEW的機械故障信號小波包分解和重構[J].機電工程,2007,24(7):14-16. YANG Shi-xi,XIANG Wen-juan,YE Hong-xian.Decomposition and Reconstruction with Wavelet Packet of Mechanical Fault Signals Based on Lab-VIEW.[J].Mechanical & Electrical Engnering Magazine,2007,24(7):14-16.(in Chinese)

[8] 吝伶艷,宋建成,謝特列.小波包頻帶檢索改進算法及其在電機故障診斷中的應用[J].太原理工大學學報,2012,43(3):391-395. LIN Ling-yan,SONG Jian-cheng,XIE Te-lie.The Improved Algorithm of Wavelet Packet Frequency Band Retrieval and Its Application in Motor Fault Diagnosis[J].Journal of Taiyuan University of Technology,2012,43(3):391-395.(in Chinese)

Conversion Algorithm from Order Number of Nodes to that of Frequency Bands for Seismic Signals Using Wavelet Packet Transform

BAI Quan1, HAN Jing-jing2, KANG Yu-mei3, BIAN Jing-mei4

(1.SchoolofArchitectureandCivilEngineering,ShenyangUniversityofTechnology,Shenyang110870,Liaoning,China; 2.SchoolofArchitecturalEngineering,HezhouUniversity,Hezhou542899,Guangxi,China; 3.CollegeofResourcesandCivilEngineering,NortheasternUniversity,Shenyang110004,Liaoning,China; 4.SchoolofEconomicsandManagement,ShenyangUniversityofChemicalTechnology,Shenyang110142,Liaoning,China)

The "frequency hopping" phenomenon often occurs when using the wavelet packet transform to process non-stationary seismic signals. To solve this problem, we examine the relationship between the order number of the wavelet packet tree bands and the order number of the nodes, based on the wavelet packet transform and exclusive or (XOR) operation algorithm. We clarify the arrangement rule between the frequency bands of the signal subspace and the tree nodes of the wavelet packet. Moreover, we propose a conversion algorithm for seismic signals from the order number of the nodes to that of the frequency bands. Based on the MATLAB platform, we use the Tangshan wave (NS direction) as an example and verify the correctness of the algorithm.

wavelet packet; tree node; frequency bands arrangement; XOR algorithm

2015-12-25 基金項目:國家自然科學基金(51204029,51308348);遼寧省教育廳基金項目(L2013050);賀州學院自然科學基金資助項目(2016ZZZK06)

白 泉,男,副教授,主要從事隨機荷載模擬及結構動力反應分析等方面的研究。E-mail:baiquan@163.com。

韓晶晶，女，助教，主要從事地震動特性分析及地震荷載模擬等方面的研究。E-mail:hanjingjing_stu@163.com。

P315.3+1

A

1000-0844(2016)06-0991-06

10.3969/j.issn.1000-0844.2016.06.0991