基于正交化HHT和隨機相位模擬非平穩隨機過程

胡燦陽，陳清軍，祁 冰，徐慶陽

(1.南京審計學院 江蘇省公共工程審計重點實驗室，南京 210029;2.同濟大學 土木工程防災國家重點實驗室，上海 200092;3.河南省建筑科學研究院有限公司，鄭州 450053)

現實生活中的很多非平穩隨機過程，不僅強度是非平穩的，而且頻率含量也是非平穩的。在地震波分析、地鐵和車輛振動分析、海浪資料分析、風的分析等很多領域里都需要模擬非平穩過程。羅雄等［1］利用觀測的風速記錄，通過快速傅氏變換和希爾伯特變換，建立了非平穩隨機過程模型。李錦華等［2］利用推導出的調制函數，運用諧波合成法模擬的非平穩脈動風速，具有強度和頻率的雙非平穩性，但是模擬計算量非常大，不具有普遍適用性。Yeh等［3］通過頻率調制的隨機過程乘以確定的時變函數來表達非平穩地震動過程。但是它使用穿零率方法來進行頻率調制函數的參數估計，并且頻率調制過程獨立于強度調制，這就使這個模型存在某些限制，只適用于一些特定類型的隨機過程。Liang等［4］基于Priestley演變譜理論，推導了模擬一維單變量非平穩隨機過程公式，并模擬了非平穩地震動過程。王國林等［5］對四輪相關車輛非平穩路面激勵的時域模型進行了建模和仿真，從而可以作為路面輸入為進一步研究提供參考。

從上面可以看到，在很多研究領域都需要非平穩過程的仿真以作為研究的輸入激勵，但在以往的研究中，人們大都把注意力放在幅值譜上，對相位譜的研究卻很少。這主要是因為相位譜要比幅值譜不規則得多［6］。時程所包含的相位特性對這兩個非平穩性都有很大影響。大崎順彥等［7］較早地強調了地震動相位信息在模擬強地震動的重要性。本文即通過正交化HHT方法和考慮隨機過程的相位［8］來進行非平穩隨機過程的模擬。

1 正交HHT變換

Huang等［9］提出了 Hilbert-Huang 變換(HHT)，這種方法特別適合非平穩信號的分析和處理。然而在實際應用中發現，文獻［9］給出的EMD分解得到的IMF分量之間不是完全正交的。由圖1可見，通過傳統的HHT變換進行非平穩地震地面運動譜密度估計時，會存在能量泄漏，EI Centro波的能量泄漏高達40%［10］。樓夢麟等［11］給出了正交化的EMD分解方法，這樣就很好的解決了EMD分解的問題，使各IMF分量之間是嚴格正交的。

通過正交化處理后，從EMD分解得到時程信號X(t)有如下分解形式:

aj(t)和 θj(t)分別為j階信號的瞬時幅值和瞬時相位，由下式計算:

把信號振幅顯示在頻率－時間平面上，就可以得到Hilbert幅值譜 H(ω，t)，稱為 Hilbert譜，記作:

圖1給出了EI Centro地震地面記錄、Taft地震地面記錄和Kobe地震地面記錄，分別采用正交化HHT法(OHHT)和文獻［9］的HHT法所估計的地震波局部譜密度的能量歸一化Husid圖［12］。

圖1 OHHT和HHT估計的地震波能量歸一化Husid圖Fig.1 Unitary Husid curve of seismic wave energy estimated by OHHT and HHT

由圖1可見，無論能量變化過程還是總能量，OHHT的Husid圖都和真實地震記錄完全吻合，說明利用正交化HHT法來估計地震波的局部譜密度具有極高的精度。圖1中的三個地震波的OHHT法不存在能量泄漏，而用HHT法時，三個地震波的能量均有泄漏情況。其中，EI Centro波的能量泄漏高達30%。

由此可知，由正交HHT法得到的Hilbert譜比傳統HHT的Hilbert譜具有更高的精度，更能反映非平穩隨機過程的時－頻特性。

2 非平穩過程的生成及統計特性研究

在工程中，一些實際物理過程不能通過數學方法來精確處理，這就要對它建立一個理論模型進行分析。主要是通過觀測的數據來建立隨機過程模型，觀測的記錄可作為隨機過程的一個樣本?；谡籋HT方法的Hilbert譜能夠精確描述這一樣本幅值和頻率成分隨時間的變化。不同樣本的Hilbert譜是不同的，因此，隨機過程的Hilbert譜可以認為是所有樣本的Hilbert譜的期望值:

然而，實際情況是有些隨機過程，比如地震這樣的隨機過程往往只能記錄到一個樣本，這時，對于隨機過程統計特性的分析和隨機過程的模擬都是困難的。在這種情況下，只能基于工程實際來作出一些假定才能對隨機過程進行分析。在建立非平穩隨機過程的模型時，參數估計是一個非常重要卻又容易出問題的環節。本文假定，已知樣本基于正交HHT方法得到的Hilbert譜即為隨機過程的Hilbert譜，這樣雖和實際情況可能還有一些偏差，但它不再需要假設和進一步的近似分析，是一種便于實現的方法。

從式(1)、式(5)看出，基于Hilbert譜可以通過引入一個隨機變量φj來構建隨機過程:

其中，φj是獨立的相位角，它服從于［0，2π］均勻分布，這樣X(t)就成為一個隨機過程，其均值、協方差和方差分別是:

根據中心極限定理，當n足夠大時，X(t)趨于高斯分布。式(7)是通過Hilbert變換后生成的隨機過程，它的m個樣本函數組成樣本空間［X1(t)，X2(t)，…，Xm(t)］。其中第k個樣本是:

第三，通過完善“互聯網+”的建設，使基礎信息數據具有整合性的特質，從而保證數據的方向型傳輸。進而使數據的收集、數據的運用、數據的儲存能在校園數據庫中進行整合、分析和處理。特別是需要根據“互聯網+”的功能分層做出相應調研，使教學、管理和信息的傳遞功能不會因數據的混亂而造成整合發生偏差。

它的均值、協方差和方差可以由式(8)～式(10)求得。可見，Hilbert譜可以方便的表示任一樣本的幅值函數和瞬時頻率。第p個和第q個樣本的互協方差函數為:

從上式可以看出，協方差函數依賴于各樣本正交IMF分量的幅值、相位函數和隨機相位偏移。比如，若兩個樣本的隨機相位角在統計意義上獨立于所有的φjp和φkq，那么互協方差函數為零。如果兩個樣本之間的相位偏移具有相關性，則可以建立它們的互協方差函數?？梢允箖蓚€樣本的相位偏移具有某種函數關系。比如可以給出兩個樣本同階IMF分量的相位偏移關系:

其中，(j=1，2，…，n，q=2，3，…，m)，ψjq是常數。這樣，不同樣本的同階IMF分量的隨機相位偏移之間是相關的，但不同階IMF分量的隨機相位偏移是獨立的。當 φj1為［0，2π］間的均勻分布時，φjq則服從于［－ψjq，2π－ψjq］上均勻分布，它的分布范圍仍然在2π之內。均值、協方差和方差仍由式(8)～式(10)給出。兩個樣本信號間的互協方差為:

由上式可見，互協方差依賴于不同樣本的同階IMF分量幅值的乘積，也和這些IMF分量的相位差與隨機相位偏移差之和相關。當這兩個差值之和非常小時，余弦函數接近1，這時，兩個樣本的互協方差只與IMF分量的參數有關。當總差值接近于π/2時，互協方差為零。在實測的基礎上，調整常數項 ψjq可以給建立互協方差函數帶來一定的靈活性。在確定的物理機制下，可以假定所有樣本之間的隨機相位偏移一致，那么樣本之間的關系符合式(14)，并且 ψjp=ψjq=0。

上面給出了描述隨機過程各樣本之間相關性的一種方法，通過正交IMF分量的Hilbert譜來估計協方差函數。由式(7)可以看出，通過Hilbert譜分析模型可以很方便地進行非平穩過程的仿真，在式(7)中，改變隨機相位角可以得到隨機過程的樣本函數。模擬樣本的Hilbert譜均值和目標Hilbet譜是一致的。樣本平均值也和目標值相一致，因此，符合模擬標準。同時，也推導出用已知樣本Hilbert譜參數表示的模擬非平穩隨機過程的統計特性函數，便于實際分析和應用。

3 算例驗證

基于EI Centro地震記錄和地鐵振動記錄的樣本Hilbert譜進行數字仿真，來驗證上述方法的正確性。表1～表4分別給出了EI Centro南北向地震地表記錄和上海地鐵1號線常熟路站的一個實測地表振動記錄時程和9個樣本時程的均值和方差?？梢钥闯瞿M樣本的均值和方差都與原記錄很接近，特別是地鐵振動樣本與記錄的方差完全一致，說明模擬樣本對均值的分散程度很好。

表1 EI Centro記錄和模擬樣本的均值Tab.1 Recorded and simlulated El Centro arithmetic average

表2 EI Centro記錄和模擬樣本的方差Tab.2 Recorded and simlulated El Centro variance

表3 地鐵振動記錄和模擬樣本的均值Tab.3 Recorded and simlulated subway vibration arithmetic average

表4 地鐵振動記錄和模擬樣本的方差Tab.4 Recorded and simlulated subway vibration variance

圖2 EI Centro記錄和模擬樣本Fig.2 Recorded and simlulated EI Centro

圖3 地鐵振動記錄和模擬樣本Fig.3 Recorded and simlulated subway vibration

圖4 El Centro記錄和模擬樣本的傅里葉譜Fig.4 Fourier spectra of EI Centro record and its samples

圖5 地鐵振動記錄和模擬樣本的傅里葉譜Fig.5 Fourier spectra of subway vibration record and its samples

圖4、圖5是原記錄和樣本的傅里葉譜。可以看出，無論是以低頻為主的地震動還是高頻的地鐵振動，每種情況的三個樣本頻率分布情況和原記錄都非常相符。因此，基于正交Hilbert譜產生的樣本也能夠很好的反映非平穩隨機過程的頻率變化特性。

為了更加清楚地對比模擬樣本與原記錄的時頻特性，圖6、圖7分別給出了地震記錄與地鐵振動記錄及樣本的小波時頻分布圖。利用Morlet小波對原記錄和相應樣本進行小波分解，得到能準確反映它們時頻非平穩特性的局部譜。由圖可見，利用本文方法模擬生成的非平穩隨機過程樣本的時頻特性都和原記錄非常符合。

圖6 EI Centro記錄和模擬樣本的小波時頻分布Fig.6 Wavelet time-frequency distribution of EI Centro record and its samples

圖7 地鐵振動記錄和模擬樣本的小波時頻分布Fig.7 Wavelet time-frequency distribution of subway vibration record and its samples

圖8給出了El Centro記錄和四個樣本的加速度反應譜，反映出模擬樣本的加速度反應譜和原記錄的加速度反應譜變化趨勢相一致。

圖8 El Centro記錄和模擬樣本的加速度反應譜Fig.8 Acceleration response spectra of El Centro ground motion record and simulations

4 結論

由于正交HHT變換克服了傳統HHT變換存在能量泄漏的缺點，它可以作為一種對非平穩信號進行時頻分析的精確方法?？梢酝ㄟ^正交HHT變換得到的Hilbert譜進行非平穩隨機過程的模擬。首先通過正交化EMD分解，進而進行Hilbert變換，得到瞬時頻率和幅值，在此基礎上進行非平穩隨機過程的模擬。這樣就克服了目前常用的頻率調制或者頻率、幅值雙調制存在的困難與不足。并且給出了非平穩隨機過程的均值、方差和協方差的Hilbert譜表示方法。本文的方法提取了非平穩隨機過程自身的物理特性，在沒有改變它的基礎上進行了非平穩隨機過程的仿真。既不需要假定一個函數作為目標譜，也不用進行分段平穩假設。文中給出了由樣本Hilbert譜參數表示的模擬隨機過程統計特性函數，具有很方便的使用性。并且通過算例驗證了本文方法進行非平穩隨機過程模擬在強度、頻率和反應譜等方面都具有很高的精度，既適用于低頻隨機過程，也適用于高頻隨機過程的模擬。算例還表明:模擬計算速度很快，可以為進一步的研究提供大量非平穩隨機過程樣本，具有很好的適用性和使用性。

［1］羅 雄，李永樂.按非平穩隨機過程模擬自然風速的一種數值方法［J］.西南交通大學學報，2002，37(4):367－370.

［2］李錦華，李春祥，申建紅.非平穩脈動風速的數值模擬［J］.振動與沖擊，2009，28(1):18－23.

［3］Yeh C H，Wen Y K.Model of nonstationary earthquake ground motion and application［C］.Proceedings of Fourth U.S.National Conference on Earthquake Engineering，1990，1:505－ 514.

［4］Liang J，Chaudhuri S R，Shinozuka M.Simulation of nonstationary stochastic processes by spectral representation［J］.J.Engrg Mech，ASCE，2007，133(6):616 －627.

［5］王國林，胡 蛟，錢金戈，等.路面對汽車非平穩激勵的時域仿真及小波分析［J］.振動與沖擊，2010，29(7):28－32.

［6］朱 昱，馮啟民.地震加速度相位差譜分布的數字特征［J］.地震工程與工程振動，1993，13(2):30 －37.

［7］大崎順彥著，呂敏申，謝禮立，譯.地震動的譜分析入門［M］.北京:地震出版社，1980.

［8］金 星，廖振鵬.地震動相位特性的研究［J］.地震工程與工程振動，1993，13(1):7 －13.

［9］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.

［10］胡燦陽，陳清軍.非平穩地震地面運動局部譜密度正交化HHT估計［J］.同濟大學學報(自然科學版)，2008，36(9):1164－1169.

［11］樓夢麟，黃天立.正交化經驗模式分解方法［J］.同濟大學學報(自然科學版)，2007，35(3):293 －298.

［12］Husid R L.Analisis de terremotos:analisis general［J］.Revista del ID1EM，Santiago Chile，1969，8(1):21－42.