基于正交HHT方法的一種高效地震波仿真研究

胡燦陽(yáng)，陳清軍，徐慶陽(yáng)

（1.南京審計(jì)學(xué)院 江蘇省公共工程審計(jì)重點(diǎn)實(shí)驗(yàn)室，南京 210029；2.同濟(jì)大學(xué) 土木工程防災(zāi)國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海 200092）

0 引言

地震動(dòng)加速度［1］是聯(lián)系地震和結(jié)構(gòu)抗震的橋梁。在進(jìn)行結(jié)構(gòu)非線性分析和模型試驗(yàn)時(shí)都需要一組具有同一統(tǒng)計(jì)特性的地震動(dòng)過(guò)程，然而地震的不確定性使我們無(wú)法得到滿足這一要求的一組地震記錄。所以根據(jù)特定目標(biāo)譜來(lái)模擬地震動(dòng)過(guò)程是地震工程中的一個(gè)重要研究方向。最常見(jiàn)的目標(biāo)譜是符合某實(shí)際地震動(dòng)記錄的反應(yīng)譜或功率譜。文獻(xiàn)［2］以反應(yīng)譜作為目標(biāo)譜，通過(guò)反應(yīng)譜和目標(biāo)譜的轉(zhuǎn)換關(guān)系來(lái)合成人工地震動(dòng)。然而反應(yīng)譜只能間接地反映地震動(dòng)過(guò)程，無(wú)法準(zhǔn)確地表現(xiàn)出地震動(dòng)的時(shí)－頻非平穩(wěn)特性，它只適合于強(qiáng)度非平穩(wěn)情況。為了能真實(shí)反映地震動(dòng)的非平穩(wěn)性質(zhì)，用時(shí)變功率譜作為目標(biāo)譜可以得到更加符合原地震記錄時(shí)頻特性的地震動(dòng)樣本，問(wèn)題的關(guān)鍵是如何獲得能夠準(zhǔn)確反映實(shí)際地震動(dòng)時(shí)－頻特性的時(shí)變譜。本文擬通過(guò)正交HHT 變換估計(jì)地震波的時(shí)變功率譜來(lái)進(jìn)行非平穩(wěn)地震地面運(yùn)動(dòng)模擬的研究。

1 基于正交HHT 變換估計(jì)地震動(dòng)的時(shí)變功率譜

Huang等［3］提出了一種特別適合非平穩(wěn)信號(hào)分析和處理的理論和方法，稱為HHT。該方法是基于信號(hào)局部特征的，因此適用于非平穩(wěn)信號(hào)的處理。相對(duì)于其它時(shí)頻分解方法（如小波分析），HHT可以獲得更為清晰的時(shí)頻譜；相對(duì)于STFT 而言，HHT 可以獲得更多頻率信息，并有效表示出信號(hào)中的瞬態(tài)成分。但是Huang 等所提出的EMD 分解方法在理論上并不能保證嚴(yán)格的正交性，使得直接利用HHT 方法進(jìn)行非平穩(wěn)地震地面運(yùn)動(dòng)局部譜密度估計(jì)時(shí)，會(huì)存在能量泄漏，由它估計(jì)的局部譜密度不能作為原地震信號(hào)的時(shí)變譜密度。而通過(guò)正交HHT［4］變換估計(jì)地震動(dòng)的時(shí)變功率譜可以使能量泄漏的問(wèn)題得到很好的解決，可以作為時(shí)變功率譜的高效估計(jì)方法。

通過(guò)對(duì)EMD 分解得到的各階IMF 分量進(jìn)行正交化處理，得到了完全正交的各階IMF 分量，其步驟如下：

（1）先將原信號(hào)進(jìn)行EMD 分解，得到n個(gè)初始IMF分量c1′（t），c2′（t）…cn′（t）和殘差rn（t）。把EMD分解得到的時(shí)程信號(hào)X（t）中除了殘差的最后一階IMF分量cn′（t），記作c1′（t）；即令c1（t）＝cn′（t），稱為時(shí)程信號(hào)X（t）的第1階正交化IMF分量；

（2）為了得到時(shí)程信號(hào)X（t）的第2 階正交化IMF分量，應(yīng)從cn－1′（t）中消除所含的c1（t）分量，即

式中，β21稱為cn－1′（t）與c1（t）之間的正交化系數(shù)，c2（t）稱為時(shí)程信號(hào)X（t）的第2階正交化IMF 分量。為了得到β21，可將式（1）的兩邊同乘c1（t）并對(duì)時(shí)間t進(jìn)行積分，同時(shí)利用c2（t）與c1（t）的正交性，得到β21，并表示成離散形式，

（3）采用與上述相同的方法，從EMD 分解得到時(shí)程信號(hào)X（t）的第n－j階初始IMF 分量中消除所含的從第1階到第j階正交化IMF 分量，則可得到時(shí)程信號(hào)X（t）的第j＋1 階正交化IMF 分量cj＋1（t）（j＝1，…，n－1）。具體的正交化計(jì)算為：

為了得到βj＋1，i，將式（3）的兩邊同乘ck（t），（k≤j）并對(duì)時(shí)間t進(jìn)行積分，利用ck（t）與ci（t），（i≠k）以及cj＋1（t）間的正交性，令i＝k時(shí)即可得到βj＋1，i，并將βj＋1，i表示成離散形式為，

經(jīng)過(guò)上述步驟以后，X（t）被分解成了如下形式，

由此，分量cj（t）（j＝1，2，…，n）之間是完全正交的，通過(guò)線性變換得到的c＊j（t）（j＝1，2，…，n）之間也是完全正交的。這樣，信號(hào)X（t）被分解成為n個(gè)正交的IMF 分量c＊j（t）（j＝1，2，…，n）及余量rn（t）的和。上述提取正交IMF分量過(guò)程并沒(méi)有改變EMD 分解的原有過(guò)程，只是在提取了IMF 后對(duì)它們進(jìn)行了正交化，即在數(shù)學(xué)意義上的重組。因此，這個(gè)正交化過(guò)程，能在嚴(yán)格保證IMF 屬性的基礎(chǔ)上，使之具有準(zhǔn)確的正交性。

對(duì)于正交后的任一IMF 分量c＊（t）作Hilbert變換即得到（t）：

式中P代表柯西主值，則對(duì)于c＊（t）的解析信號(hào)z（t）為：

a（t）和θ（t）分別為信號(hào)的瞬時(shí)頻率和瞬時(shí)相位。把信號(hào)振幅顯示在頻率－時(shí)間平面上，就可以得到Hilbert幅值譜H（ω，t），稱為Hilbert譜，記作：

H（ω，t）精確地描述了信號(hào)的幅值隨時(shí)間和頻率的變化規(guī)律。由此可以得到原隨機(jī)信號(hào)的局部功率譜密度［5］，即

由此可見(jiàn)，通過(guò)正交HHT 變換（OHHT）的時(shí)頻特性，將隨機(jī)過(guò)程或樣本的能量進(jìn)行了時(shí)頻局部化，即其能量在時(shí)間和頻率構(gòu)成的平面上被展開(kāi)，如（9）式，從而得到了隨機(jī)過(guò)程或樣本的時(shí)頻局部特性——局部譜密度。信號(hào)通過(guò)正交HHT 變換后得到的Hilbert譜與局部譜密度的這種簡(jiǎn)單的關(guān)系，使利用正交HHT 變換進(jìn)行局部譜密度的估計(jì)非常方便。

利用Husid［6］圖可以定量的檢驗(yàn)估計(jì)局部譜密度對(duì)地震波的時(shí)間局部特性反映的準(zhǔn)確性。利用（9）式可得到地震波能量在時(shí)間上的分布，把每一時(shí)刻和此前所有時(shí)刻的能量累加可以得到地震波的Husid圖。通過(guò)對(duì)比它和實(shí)際地震波的能量Husid圖就可以定量的檢驗(yàn)正交HHT 變換所估計(jì)的地震波局部譜密度對(duì)強(qiáng)度非平穩(wěn)估計(jì)的準(zhǔn)確性。圖1為對(duì)三條著名地震波分別采用正交HHT 法和常規(guī)HHT 法所估計(jì)的地震波局部譜密度的能量歸一化Husid圖。

圖1 OHHT和HHT估計(jì)的地震波能量歸一化Husid圖

地震波Husid圖是能量隨時(shí)間積累的曲線，反映了地震波強(qiáng)度隨時(shí)間變化的局部特性。從圖1可以看出，用正交HHT 法估計(jì)地震地面運(yùn)動(dòng)的局部譜密度得到的Husid圖模擬值與地震波真實(shí)值吻合得很好，而且模擬值總能量和真實(shí)值一致，說(shuō)明正交HHT 法在估計(jì)局部譜密度時(shí)沒(méi)有能量泄漏，表明了它在反映地震波能量在時(shí)間上的分布方面有著很高的精度。而直接用HHT 法對(duì)這三條波估計(jì)時(shí)都有較大的能量泄露，EI Centro波的能量泄露高達(dá)30%。因此，在對(duì)原信號(hào)局部功率譜密度進(jìn)行估計(jì)時(shí)，各IMF的正交性是必需的。

2 基于正交HHT 變換的非平穩(wěn)地震地面運(yùn)動(dòng)模擬

Scanlan和Sachs［7］建議了利用三角級(jí)數(shù)法和演變譜模擬地震波的表達(dá)式：

為了模擬地震記錄的時(shí)頻非平穩(wěn)特性，這里f（t，ωk）采用時(shí)變譜，Δω＝（ωk－ωk－1）為時(shí)變譜的頻率間隔，Φk是［0，2π］區(qū)間均勻分布、相互獨(dú)立的隨機(jī)相位角。用普遍使用的多重過(guò)濾［8］或短時(shí)傅立葉變換（STFT）［9］估計(jì)地震波時(shí)變譜不僅速度慢、精度不高，而且如果在譜估計(jì)時(shí)為了提高精度，需要把Δω取得很小，導(dǎo)致在估計(jì)時(shí)變譜時(shí)非常費(fèi)時(shí)，難以滿足大量地震波模擬的需要。由圖1可見(jiàn)，常規(guī)HHT 功率譜存在能量泄漏，無(wú)法作為目標(biāo)譜來(lái)模擬地震波，本文使用正交HHT 功率譜作為時(shí)變譜進(jìn)行地震波的模擬。利用正交化后HHT 變換可以快速、準(zhǔn)確地估計(jì)出能夠全面反映地震波時(shí)－頻非平穩(wěn)特性的時(shí)變譜，其精度和速度明顯高于STFT和多重過(guò)濾，使得快速進(jìn)行大量非平穩(wěn)地震波模擬成為可能。

以El Centro地震波和1條Landers地震記錄為目標(biāo)進(jìn)行模擬。圖2為El Centro波地震記錄和生成的3條樣本。圖3為L(zhǎng)anders波地震記錄及生成的3條樣本。

由圖2和圖3可以看出，2 條波的樣本曲線都能反映和原地震記錄一樣的強(qiáng)度隨時(shí)間的變化趨勢(shì)，特別是Landers波的2組加速度峰值的大小和分布都得到了很好的再現(xiàn)。這里每條樣本的峰值和原記錄不一定完全一致，但是可以從統(tǒng)計(jì)上來(lái)研究模擬波強(qiáng)度峰值的精度。為此，同時(shí)用多重過(guò)濾功率譜和STFT 功率譜作為目標(biāo)譜來(lái)模擬地震波，同樣各生成1000個(gè)樣本。限于篇幅，這里不做出它們的樣本曲線。圖4和圖5分別是用這3種目標(biāo)功率譜生成1000個(gè)樣本的標(biāo)準(zhǔn)差。

圖2 El Centro波地震記錄和生成的樣本

圖3 Landers波地震記錄和生成的樣本

圖4 El Centro波1000個(gè)樣本標(biāo)準(zhǔn)差

圖5 Landers波1000個(gè)樣本標(biāo)準(zhǔn)差

對(duì)于每個(gè)地震波，不同樣本的峰值大小、峰值個(gè)數(shù)和峰值出現(xiàn)位置都不相同。以樣本的加速度峰值大小為例：對(duì)于El Centro波，正交HHT 功率譜生成3個(gè)樣本的最大峰值分別是相應(yīng)標(biāo)準(zhǔn)差的2.09、1.63、1.46倍；對(duì)應(yīng)于多重過(guò)濾功率譜，3個(gè)樣本的最大峰值分別是相應(yīng)標(biāo)準(zhǔn)差的2.20、3.00、2.69倍；對(duì)應(yīng)于STFT 功率譜，3個(gè)樣本的最大峰值分別是相應(yīng)標(biāo)準(zhǔn)差的2.72、1.93、2.99 倍。對(duì)于Landers波，它有2組峰值。對(duì)應(yīng)于正交HHT 功率譜，3個(gè)樣本的第1 組峰值分別是相應(yīng)標(biāo)準(zhǔn)差的2.38、2.14、1.83 倍；第2 組 峰 值 分 別 是 相 應(yīng) 標(biāo) 準(zhǔn) 差 的2.17、1.84、1.82倍。對(duì)應(yīng)于多重過(guò)濾功率譜，3個(gè)樣本的第1 組峰值分別是相應(yīng)標(biāo)準(zhǔn)差的3.55、2.12、2.52倍；第2 組 峰 值 分 別 是 相 應(yīng) 標(biāo) 準(zhǔn) 差 的3.34、2.08、2.56倍。對(duì)應(yīng)于STFT 功率譜，3個(gè)樣本的第1組峰值分別是相應(yīng)標(biāo)準(zhǔn)差的2.96、2.32、2.28 倍；第2 組峰值分別是相應(yīng)標(biāo)準(zhǔn)差的2.56、2.09、2.67倍。從樣本過(guò)程峰值大小的離散程度來(lái)看，模擬產(chǎn)生地震波的強(qiáng)度大小比較理想。

由上面可以看出，以正交HHT 功率譜作為目標(biāo)譜并通過(guò)三角級(jí)數(shù)法生成的地震波能很好地反映原地震記錄的強(qiáng)度非平穩(wěn)特性。

此外，從圖2和圖3也可以直觀地看出樣本頻率隨時(shí)間的變化趨勢(shì)和原記錄相符合。圖6和圖7分別是2個(gè)波和任意一個(gè)模擬樣本的傅立葉幅值譜。由圖可見(jiàn)，樣本和原記錄的頻譜特性比較接近。為了更好地檢驗(yàn)強(qiáng)震記錄和生成樣本的時(shí)－頻非平穩(wěn)特性，圖8和圖9分別為2個(gè)波和它們的樣本在時(shí)－頻面上的Hilbert譜，由此可以看出能量在時(shí)－頻面上的分布情況。從圖上可以看出，樣本的Hilbert時(shí)－頻譜分布和原記錄吻合的較好，樣本能較準(zhǔn)確地反映原強(qiáng)震記錄的時(shí)－頻非平穩(wěn)特性。

圖6 El Centro波和模擬樣本的傅氏頻譜

圖7 Landers波和模擬樣本的傅氏頻譜

圖8 El Centro波和模擬樣本的Hilbert譜

圖9 Landers波和模擬樣本的Hilbert譜

3 結(jié)論

在總結(jié)常用非平穩(wěn)地震時(shí)變功率譜估計(jì)的基礎(chǔ)上，本文通過(guò)正交HHT 功率譜和三角級(jí)數(shù)法來(lái)模擬非平穩(wěn)地震波。正交HHT 功率譜避免了常規(guī)HHT 功率譜的能量泄漏，在時(shí)間和頻率上都能很好地反映原地震記錄的非平穩(wěn)特性，因此以它作為目標(biāo)譜生成樣本的強(qiáng)度和頻率都是非平穩(wěn)的。任意一個(gè)樣本過(guò)程的時(shí)變譜不一定符合目標(biāo)譜，但是在統(tǒng)計(jì)意義上嚴(yán)格符合。算例表明，本文的方法既能滿足地震波非平穩(wěn)的宏觀特性，又在數(shù)值上具有很好的模擬精度，這樣可以為某一強(qiáng)震記錄補(bǔ)充大量具有相同統(tǒng)計(jì)特性的樣本波，有利于結(jié)構(gòu)隨機(jī)響應(yīng)時(shí)程分析。

正交化HHT 法可以有效地避免EMD 分解時(shí)的能量泄漏，但不能解決端點(diǎn)效應(yīng)［10］等HHT 自身的缺陷，所以正交HHT 功率譜還有可能存在失真現(xiàn)象。如果能夠解決端點(diǎn)效應(yīng)問(wèn)題，本方法將會(huì)很好地應(yīng)用于地震波仿真和其它非平穩(wěn)信號(hào)（如風(fēng)波、水波、鐵路和公路的振動(dòng)信號(hào)等）的模擬，會(huì)具有廣泛地應(yīng)用前景。

［1］ 肖本夫，萬(wàn)永革，祁玉萍.強(qiáng)干擾環(huán)境下有效數(shù)字地震信號(hào)的提?。跩］.華北地震科學(xué)，2011，（1）：6－9.

［2］ 蔣溥，戴麗思.工程地震學(xué)概論［M］.北京：地震出版社，1993：90－94.

［3］ Huang N E，Shen Z，Long S R，et al.The empirical mode decomposition and the Hilbert spectrum for nonlinear and non－stationary time series analysis［J］.Proc.R.Soc.Lond.A，1998，454：903－995.

［4］ 樓夢(mèng)麟，黃天立.正交化經(jīng)驗(yàn)?zāi)Ｊ椒纸夥椒ǎ跩］.同濟(jì)大學(xué)學(xué)報(bào)（自然科學(xué)版），2007，35（3）：293－298.

［5］ 胡燦陽(yáng)，陳清軍.基于HHT 法的地震地面運(yùn)動(dòng)局部譜密度估計(jì)［J］.振動(dòng)與沖擊，2007，26（10）：126－131.

［6］ Husid R L.Analisis de Terremotos：Analisis General［J］.Revista del ID1EM，1969，Santiago Chile，8（1）：21－42.

［7］ Shinozuka J，Jan C M.Simulation of random processes and its applications［J］.J.Sound and Vibration，1972，25（1）：111－128.

［8］ Priestly M B.Power spectral analysis of non－stationary random process［J］.J.Sound and Vibration，1967，6：86－97.

［9］ Mark W D.Spectral analysis of the convolution and filtering of non－stationary stochastic processes［J］.J.Sound and Vibration，1970，11：19－63.

［10］ 羅奇峰，石春香.Hilbert－Huang變換理論及其計(jì)算中的問(wèn)題［J］.同濟(jì)大學(xué)學(xué)報(bào)（自然科學(xué)版），2003，31（6）：637－640.