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

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

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

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

1 正交HHT變換

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

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

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

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

圖1給出了EI Centro地震地面記錄、Taft地震地面記錄和Kobe地震地面記錄，分別采用正交化HHT法(OHHT)和文獻(xiàn)［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波的能量泄漏高達(dá)30%。

由此可知，由正交HHT法得到的Hilbert譜比傳統(tǒng)HHT的Hilbert譜具有更高的精度，更能反映非平穩(wěn)隨機(jī)過程的時－頻特性。

2 非平穩(wěn)過程的生成及統(tǒng)計特性研究

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

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

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

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

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

第三，通過完善“互聯(lián)網(wǎng)+”的建設(shè)，使基礎(chǔ)信息數(shù)據(jù)具有整合性的特質(zhì)，從而保證數(shù)據(jù)的方向型傳輸。進(jìn)而使數(shù)據(jù)的收集、數(shù)據(jù)的運(yùn)用、數(shù)據(jù)的儲存能在校園數(shù)據(jù)庫中進(jìn)行整合、分析和處理。特別是需要根據(jù)“互聯(lián)網(wǎng)+”的功能分層做出相應(yīng)調(diào)研，使教學(xué)、管理和信息的傳遞功能不會因數(shù)據(jù)的混亂而造成整合發(fā)生偏差。

它的均值、協(xié)方差和方差可以由式(8)～式(10)求得?？梢姡琀ilbert譜可以方便的表示任一樣本的幅值函數(shù)和瞬時頻率。第p個和第q個樣本的互協(xié)方差函數(shù)為:

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

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

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

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

3 算例驗證

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

表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譜產(chǎn)生的樣本也能夠很好的反映非平穩(wěn)隨機(jī)過程的頻率變化特性。

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

圖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記錄和四個樣本的加速度反應(yīng)譜，反映出模擬樣本的加速度反應(yīng)譜和原記錄的加速度反應(yīng)譜變化趨勢相一致。

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

4 結(jié)論

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

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

［2］李錦華，李春祥，申建紅.非平穩(wěn)脈動風(fēng)速的數(shù)值模擬［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］王國林，胡 蛟，錢金戈，等.路面對汽車非平穩(wěn)激勵的時域仿真及小波分析［J］.振動與沖擊，2010，29(7):28－32.

［6］朱 昱，馮啟民.地震加速度相位差譜分布的數(shù)字特征［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］胡燦陽，陳清軍.非平穩(wěn)地震地面運(yùn)動局部譜密度正交化HHT估計［J］.同濟(jì)大學(xué)學(xué)報(自然科學(xué)版)，2008，36(9):1164－1169.

［11］樓夢麟，黃天立.正交化經(jīng)驗?zāi)Ｊ椒纸夥椒ǎ跩］.同濟(jì)大學(xué)學(xué)報(自然科學(xué)版)，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.