非平穩(wěn)隨機過程功率譜密度估計的小波方法

孔 凡，李 杰，2

（1．同濟大學(xué)土木工程學(xué)院，上海 200092；2．同濟大學(xué)土木工程防災(zāi)國家重點實驗室，上海 200092）

引 言

工程結(jié)構(gòu)在役期間，會受到各種非平穩(wěn)隨機動力的作用，如地震、強風(fēng)以及海浪等。其非平穩(wěn)性不僅表現(xiàn)在依賴于時間的幅值上，而且其頻譜特性也是隨時間變化的［1］。由平穩(wěn)隨機過程經(jīng)過確定性函數(shù)調(diào)制所得的非平穩(wěn)過程并不能合理反映頻域非平穩(wěn)性質(zhì)［2］。獲取非平穩(wěn)隨機過程的時間－頻率特性，成為描述隨機動力激勵、計算結(jié)構(gòu)隨機動力響應(yīng)的重要一環(huán)。

最早的頻率分析工具為Fourier變換：它在頻域內(nèi)分辨率可視為δ函數(shù)，沒有時間分辨率［3，4］，或者精確地講，其時間信息湮滅在相位信息之中［3］。為了改進(jìn)Fourier變換的不足，人們提出了Gabor展 開／短 時 Fourier 變 換 （Short Time Fourier Transform，STFT）以及 Wigner－Ville分布（Wigner－Ville Distribution，WVD）方法［5～9］。這些方法雖然能體現(xiàn)非平穩(wěn)過程的某些時間－頻率特征，但它們都有各自的缺點和不足［4］。近30年來發(fā)展起來的小波分析［3，10～13］方法，因其具有良好的時－頻分辨率、豐富的小波基以及逐漸完備的數(shù)學(xué)背景，使其在工程方面得到了廣泛的應(yīng)用。

在利用小波分析估計非平穩(wěn)功率譜密度方面，Basu等人建立了小波變換系數(shù)的均方值與功率譜密度之間的關(guān)系［14～17］；Liang等人在文獻(xiàn)［18］中使用了這一方法估計由非平穩(wěn)過程譜表達(dá)方法合成的非平穩(wěn)隨機過程的功率譜密度。值得指出的是，Basu關(guān)系式只適用于頻域內(nèi)互不重疊的小波族［19］。Spanos等人發(fā)展了利用正交且頻域內(nèi)不重疊的諧和小波（Harmonic Wavelet，HW）以及廣義諧和小波（Generalized Harmonic Wavelet，GHW）估計非平穩(wěn)隨機過程功率譜密度的方法（簡稱Spanos－Tratskas方 法）［10，20～26］。隨 后，結(jié) 合 Wold－Cramer非平穩(wěn)隨機過程［27，28］，Spanos和Failla提出了適合一般非正交小波（簡稱為“一般小波”）估計非平穩(wěn)隨機過程演變功率譜密度的方法（簡稱Spanos－Failla方法）［29］；Huang等人把這種方法推廣于多維隨機過程以及強風(fēng)的演變功率譜估計［19］；Chakraborty等人采用這種方法計算了大跨橋梁的隨機非平穩(wěn)響應(yīng)［30］；基 于 Spanos－Failla方 法，F(xiàn)ailla等 人 發(fā) 展 了非平穩(wěn)隨機過程小波譜的思想［31］，并且證明此方法能推廣應(yīng)用于一般非Cramer－Wold模型的非平穩(wěn)隨機過程。

本文首先證明，基于諧和小波和廣義諧和小波估計隨機過程功率譜密度的Spanos－Tratskas方法為Spanos－Failla方法的特殊情況；比較了文獻(xiàn)［14～16］所提之修正的 Littlewood－Paley（Modified Littlewood－Paley，MLP）小波與GHW在估計非平穩(wěn)隨機過程EPSD方面的特點。研究表明，當(dāng)使用Spanos－Failla方法估計功率譜密度時，小波尺度的選取、小波邊緣效應(yīng)以及非平穩(wěn)過程調(diào)制函數(shù)的變化快慢均對功率譜的估計有較大影響。數(shù)值計算表明，補零延長（Zeros Padding）對于邊緣效應(yīng)有較好的消除作用。隨后，在最近提出的局部平穩(wěn)小波過程（Locally Stationary Wavelet process，LSW）的基礎(chǔ)上［32，33］，提出了一種新的估計非平穩(wěn)隨機過程的方法，它不僅適用于正交小波，而且適用于一般非正交小波。有趣的是，本文建議方法不僅可與Spanos－Failla方法在形式上統(tǒng)一起來，而且廣義諧和小波的應(yīng)用，使二者退化為同一形式。簡言之，本文從隨機過程LSW模型的角度，證明了Spanos－Tratskas方法不僅適用于基于 Wold－Cramer模型的Priestley演變功率譜密度，而且可以估計任意非平穩(wěn)隨機過程的時變功率譜密度。這一點，與作者剛注意到的文獻(xiàn)［31］有類似的結(jié)論：Failla從另一個角度證明了Failla－Spanos方法能適用于任意非平穩(wěn)隨機過程。為了驗證所提之方法的有效性，文章最后給出了基于GHW且在考慮場地土與行波效應(yīng)情況下的多維地震動互／自功率譜的數(shù)值算例。在確定性時程時－頻譜分析方面，以汶川地震中離斷層最近的綿竹清平波和Ⅲ類遠(yuǎn)場地的西安波為例，利用GHW對二者進(jìn)行了時頻譜分析，討論了導(dǎo)致近場地和遠(yuǎn)場地時－頻譜迥異的物理原因。

1 非平穩(wěn)隨機過程的EPSD估計

考慮形非平穩(wěn)隨機過程的Wold－Cramer模型

的雙邊功率譜密度。Priestley為此類非平穩(wěn)隨機過程定義了（雙邊）演變功率譜密度［28］

可以證明，非平穩(wěn)隨機過程f（t）的演變功率譜密度可以表達(dá)為

式中W（ar，b）為非平穩(wěn)隨機過程f（t）在尺度ar和時間b處的小波變換，且

式（5）可視為非平穩(wěn)隨機過程小波譜估計值［31］：在特定時間點b處，隨機過程的瞬時功率譜密度由各階小波的Fourier變換之模的平方構(gòu)造而成，cj（b）為瞬時各階加權(quán)值。

另一方面，當(dāng)使用正交且在頻域內(nèi)緊支無重疊的GHW估計功率譜時，Spanos等人給出了如下隨機過程功率譜與其小波變換系數(shù)之間的關(guān)系式［23～26］

式中 （mj，nj）為GHW的尺度因子；k為其時間平移因子；Wψ［（mj，nj），k］為非平穩(wěn)隨機過程在相應(yīng)尺度和時間點上的 GHW 變換且mjΔω＜ωj≤njΔω。

可證明當(dāng)式（5）中一般小波退化為GHW時，與式（8）是統(tǒng)一的。事實上，當(dāng)小波為GHW時，由于其在頻域內(nèi)的特殊性，僅依賴于小波形式的矩陣Q為對角陣，且

結(jié)合式（9）與式（6），易證基于 GHW 的 Spanos－Failla方法退化為如式（8）所示的Spanos－Tratskas方法。具體而言，由式（9）可知，頻域內(nèi)緊支且不重疊導(dǎo)致了矩陣Q的對角化，進(jìn)而導(dǎo)致了系數(shù)cj（b）與小波變換均方值E［｜W（aj，b）｜2］具有正比關(guān)系。非平穩(wěn)隨機過程小波譜思想和GHW的引入，更清楚地表明了非平穩(wěn)隨機過程瞬時功率譜密度與各階GHW的Fourier變換均方值之間的關(guān)系，如圖1所示。顯然，GHW對瞬時功率譜的估計類似于函數(shù)的分段常數(shù)插值。值得指出的是，類似于式（8），文獻(xiàn)［14～16］利用MLP小波從不同的角度也得到了類似于Spanos－Tratskas方法的非平穩(wěn)隨機過程瞬時PSD與小波系數(shù)模的均方值之間的關(guān)系

圖1 非平穩(wěn)隨機過程瞬時EPSD與其GHW估計Fig．1 Instantaneous PSDof non－stationary stochastic process and the GHW estimation

其中，K為與小波容許條件［11］以及離散格式相關(guān)的常數(shù)。

基于上述背景，有如下注記：

（1）MLP小波是一種正交二進(jìn)小波（Dyadic Wavelet），為諧和小波的實部。因此，式（10）顯然適用于MLP小波。而對于其他在頻域內(nèi)非緊支且重疊的小波，則宜使用Spanos－Failla方法估計演變功率譜。導(dǎo)致這種現(xiàn)象發(fā)生的原因為：其一，由于GHW或MLP小波在頻域內(nèi)的特別性，使其具有明確的尺度－頻率關(guān)系，而其他在頻域內(nèi)非緊支且重疊的小波（簡稱為“其他小波”）沒有明顯的這種關(guān)系；其二，頻域內(nèi)一般小波在不同尺度處重疊部分的能量被重復(fù)計算，導(dǎo)致了利用其他小波和式（10）所估計的演變功率譜密度大于目標(biāo)演變功率譜密度［19］。圖2所示為不同尺度Morlet小波在頻域內(nèi)的表達(dá)，清楚地反映了這一點。

圖2 不同尺度的小波在頻域內(nèi)的重疊Fig．2 Mean square value of Fourier transform of Morlet wavelets at different scales

（2）雖然利用實值的MLP小波能較好地估計單變量隨機過程的演變自功率譜密度，卻不能估計多變量隨機過程的演變互功率譜密度。因此，對于多維／多變量非平穩(wěn)隨機過程演變互功率譜密度的估計，復(fù)值的GHW顯然更具優(yōu)勢。除此之外，由于MLP小波是二進(jìn)小波，不同尺度小波在頻域內(nèi)所占頻寬依2的倍數(shù)增長。當(dāng)隨機過程的演變功率譜為慢變、寬頻且峰值出現(xiàn)在低頻處時，MLP小波尚能較好估計；如果演變功率譜在高頻處出現(xiàn)快變峰值，MLP小波不能滿足要求，如圖3所示。因此，雖然文獻(xiàn)［15，16］中所建議的MLP小波，限制了其母小波頻寬，因此對于高頻峰值的EPSD估計起到一定改進(jìn)效果，但小波個數(shù)的增加卻帶來了額外的計算量。

圖3 L－P小波對高頻峰值的瞬時功率譜的估計Fig．3 L－P wavelet based estimation of EPSD with a peak at high frequency

（3）在利用Spanos－Failla方法估計EPSD時，適當(dāng)選擇小波尺度至關(guān)重要。計算表明，尺度因子選擇不當(dāng)會導(dǎo)致低頻處估計功率譜失真或估計功率譜異常波動。圖4（a）所示為在不同尺度選擇下，基于Morlet小波的一致調(diào)制非平穩(wěn)隨機過程的功率譜估計，調(diào)制前的平穩(wěn)PSD選為Kanai－Tajimi譜。從圖4（b）易見不同尺度下Morlet小波的頻域峰值構(gòu)成：由于不當(dāng)?shù)男〔ǔ叨冗x擇，導(dǎo)致圖4（b）中不同尺度小波的Fourier變換異常明顯，從而不能較好地逼近目標(biāo)演變功率譜。計算實踐表明，小波尺度的選擇決定了矩陣Q的條件數(shù)，而后者直接影響了功率譜估計的精確程度。

圖4 不同小波尺度選擇下目標(biāo)與估計功率譜在7s時的對比Fig．4 Comparisons of target and estimated PSDs at 7s，with different wavelet scales

圖5 無數(shù)據(jù)延長和補零延長時非平穩(wěn)隨機過程估計功率譜密度Fig．5 PSD estimations of nonstationary stochastic process with different data－paddings

2 基于局部平穩(wěn)小波過程的功率譜估計

具有嚴(yán)格數(shù)學(xué)基礎(chǔ)、如式（11）所示的 Wold－Cramer非平穩(wěn)隨機過程模型，根本上是由經(jīng)典的平穩(wěn)隨機過程的譜分解理論

擴展而來，其非平穩(wěn)性由時間慢變的調(diào)制函數(shù)A（ω，t）描 述。Dahlhaus 在 非 平 穩(wěn) 隨 機 過 程 的Wold－Cramer模型的同一理論框架下，提出了一種局部平穩(wěn)過程（Local Stationary Process）［35］。直觀地講，當(dāng)非平穩(wěn)隨機過程的局部區(qū)間能以平穩(wěn)過程近似表達(dá)時，此過程可以稱為局部平穩(wěn)隨機過程。由于Dahlhaus模型使用的展開基仍為Fourier基函數(shù)，其整體非平穩(wěn)性由局部平穩(wěn)的時變傳遞函數(shù)（ω）提供［35，36］。實際上，F(xiàn)ourier其函數(shù)由于其頻域分辨率很高，卻無時間分辨率，本質(zhì)上無法作為非平穩(wěn)隨機過程模型的基函數(shù)；且根據(jù)Balian－Low定理，對Fourier基函數(shù)加以光滑緊支窗雖然能改善其分辨率，卻失去了Fourier基函數(shù)的正交性［3］。基于此考慮，Ombao等人在通過對Fourier基函數(shù)施加兩個特殊構(gòu)造的光滑窗后，提出了一種在時－頻域內(nèi)緊支且正交的光滑局部指數(shù)基（Smooth Localized complex EXponential basis，SLEX），這種基函數(shù)可以認(rèn)為是Fourier基的局部形式［36］。通過以SLEX基代替式（11）中的Fourier基，Ombao等人得到了一種類于式（11）離散形式的非平穩(wěn)隨機過程的SLEX模型。與Dahlhaus模型不同的是，SLEX模型提供了一種對時變功率譜在顯式劃分的時－頻域上分段常數(shù)插值的途徑［36］。因此，SLEX模型也描述了一種局部平穩(wěn)的非平穩(wěn)隨機過程。小波基函數(shù)由于同時具有一定的時－頻域分辨率，Nason等人和Eckley等人提出了一種以非抽樣小波為基函數(shù)的局部平穩(wěn)小波過程（Locally Stationary Wavelet process，LSW）［32，33］。實質(zhì)上，此模型可視為以小波基代替Fourier基后式（11）的離散形式，它直接定義了一種反映非平穩(wěn)隨機過程能量在不同時間－尺度上分布的小波譜。因此本節(jié)的出發(fā)點在于，如果在LSW模型中引入具有尺度－頻率顯式關(guān)系的小波基，如諧和小波或廣義諧和小波，可望得到一種簡單的估計非平穩(wěn)隨機過程的時變功率譜密度的方法。

局部平穩(wěn)小波過程表達(dá)為［33］

式中ξj，k為歸一化正交隨機序列；ψj，k為尺度變換和時間平移后的非抽樣小波基函數(shù)；wj，k為相應(yīng)尺度和時間處的幅值系數(shù)；∑j與∑k分別代表對小波和尺度平移的離散時間點求和。

在如式（12）中的平移時間點和尺度下，上述非平穩(wěn)隨機過程的小波變換為

其中上標(biāo)＊表示復(fù)共軛，且

結(jié)合式（12）與（13），當(dāng)LSW 中的小波與式（13）中的小波相同時，則

式中Il，i（j，k）為此小波的重建核（Reproducing Kernel）［37］，即

將式（15）兩邊乘以其復(fù)共軛并求期望，在考慮到ξl，i的正交歸一性后，可得

式中 小波的重建核Il，i（j，k）與式（6）中的Qr，s對應(yīng)。

另一方面，將式（13）兩端進(jìn)行小波內(nèi)積并求期望［37］，

同時在考慮到小波重建核的性質(zhì)

和ξl，i的正交歸一性后，式（18）可以寫為

因此，非平穩(wěn)隨機過程的瞬時功率譜表達(dá)為

不難看出，非平穩(wěn)LSW過程的譜估計式（21）與非平穩(wěn)隨機過程演變功率譜估計式（5）類似：前者的時變系數(shù)由式（17）求出，后者的時變系數(shù)由式（6）求出。因此，二者在形式上可視為是統(tǒng)一的。

當(dāng)LSW模型中的小波函數(shù)為GHW時，考慮到GHW在頻域內(nèi)的性質(zhì)，易證

將式（22）代入式（17），得

綜合式（22）與（23），可證式（21）退化為式（8）。由式（8）可知，兩非平穩(wěn)隨機過程的互功率譜密度可以估計為

因此，本文所建議方法的意義在于，從數(shù)學(xué)上解決了LSW非平穩(wěn)隨機過程的時變譜估計問題。對于以GHW為基函數(shù)的LSW非平穩(wěn)隨機過程，其時變譜估計方法在形式上與Priestley演變功率譜的GHW估計方法相同。換言之，式（8）不僅適用于Priestley演變功率譜估計，而且適用于基于LSW模型的隨機過程時變譜估計。

3 數(shù)值算例

3．1 模擬多變量地震動自／互功率譜估計

前已述基于一般小波的非平穩(wěn)隨機過程（Wold－Cramer和LSW非平穩(wěn)隨機過程模型）功率譜估計，以及它們的GHW退化形式，顯示了GHW在估計非平穩(wěn)隨機過程時變功率譜方面的優(yōu)勢。為了驗證GHW在估計多變量非平穩(wěn)隨機過程互／自功率譜的可靠性，以同一地震中沿地震波主要傳播方向分布的3個不同場地為例（如圖6所示），在考慮行波效應(yīng)和場地土條件不同的情況下，估計了三變量非平穩(wěn)隨機過程的互／自演變功率譜密度。場地條件分別為巖石或硬土、低粘性土以及粘性土或細(xì)砂。以Clough－Penzien譜表征平穩(wěn)地震動的功率譜；非平穩(wěn)隨機過程樣本由平穩(wěn)樣本乘以考慮行波效應(yīng)的一致調(diào)制函數(shù)產(chǎn)生。考慮場地土條件和行波效應(yīng)的三變量平穩(wěn)隨機過程樣本由譜表達(dá)方法給出［38］，樣本個數(shù)為500個，時間步長為Δt＝0．02s，時程取樣點為1 024個，譜表達(dá)方法的諧和項取為512個，截止頻率為 Nyquist頻率ωu＝π／Δtrad／s。為了達(dá)到時間－頻率分辨率的平衡，不同尺度GHW的頻寬選擇為ni－mi＝12。鑒于篇幅限制，考慮行波效應(yīng)和場地土條件等參數(shù)的多變量譜表達(dá)方法在本文中不再重述，其設(shè)置可參見文獻(xiàn)［38］。

圖6 同一地震中的3個不同場地條件Fig．6 Three different site conditions in a same earthquake event

圖7（a）～（c）所示分別為硬土、低粘性土及粘性土場地上在6s處的目標(biāo)與估計自功率譜；圖7（d）～（f）所示分別為軟粘性土－硬土、低粘性土－硬土及軟粘性土－低粘性場地上在6s處的目標(biāo)與估計互功率譜圖。可見，基于GHW能較好地給出多變量非平穩(wěn)隨機過程的自／互功率譜估計。圖8（a）～（d）所示為軟粘性土和硬土場地上，互譜實部與虛部的目標(biāo)和估計譜值。限于篇幅，其他場地上自／互譜的實／虛部的目標(biāo)和估計譜值，茲不一一列出。

圖7 不同場地條件下6s時非平穩(wěn)地震動瞬時自／互功率譜密度Fig．7 Instantaneous auto／cross－PSDs at 6sof non－stationary ground motion on different sites

圖8 軟粘性土和硬土場地上互譜實部與虛部的目標(biāo)和估計譜值Fig．8 Target and estimated cross－PSDs on soft clay－stiff soil site

最后，以粘性土場地上地震動的自功率譜為例，圖9對比了歸一化后功率譜密度的0階矩根，即隨機過程目標(biāo)均方值與估計值對比；表1為各時間點處估計均方值與目標(biāo)均方值的相對誤差。

可見，估計均方相對誤差較大處一般發(fā)生在地震動開始、峰值以及結(jié)尾處，其中以地震動開始處為甚，除此之外，其最大誤差均未超過10%。由于在實際工程中更關(guān)心的是峰值時間段處隨機地震動的統(tǒng)計特性，可見所建議方法精度能滿足工程實際需要。

表1 軟粘性土場地上各時間點處估計均方值與目標(biāo)瞬時均方值的相對誤差Tab．1 Relative errors between the target and estimated mean sqaure value of the process on the soft clay site

圖9 軟粘性土場地上歸一化估計均方值與目標(biāo)均方值的對比Fig．9 Comparison of the normalized target and estimated mean square value of the process on the soft clay site

3．2 汶川8．0級地震加速度時程的能量時間－頻率分布

本文所建議方法不僅可用于隨機過程的時變功率譜估計，而且可以對確定性時程進(jìn)行時間－頻率分析，以得到確定性時程能量的時間－頻率分布信息。

2008年5月12日在四川省汶川境內(nèi)發(fā)生了Ms8．0級特大地震，震中位于北緯31．021°，東經(jīng)103．367°，震源深度14km。本文選取本次地震中離龍門山斷層最近（斷層距0．74km）的四川綿竹清平臺和沿斷層破裂方向且離震中較遠(yuǎn)（震中距637km）的西安臺所測得的東西向強震加速度為例，利用本文所建議方法估計了地震動能量的時間－頻率分布。圖10為兩地測得的強震加速度時程；圖11為基于GHW變換的兩地加速度時程能量的時間－頻率分布。

圖10 汶川地震中不同場地地震動加速度記錄Fig．10 Acceleration records of ground motion in Wenchuan earthquake

圖11 汶川地震動中不同場地上加速度記錄能量的時間－頻率分布Fig．11 Energy distributions of acceleration records in Wenchuan earthquake

由圖10僅能得出兩地地震動加速度時程最大值：其中離斷層較近的綿竹清平臺東西向最大加速度達(dá)到了814．13gal，離震中較遠(yuǎn)的西安臺東西向最大加速度只有52．71gal。二者能量的時間－頻率分布卻有很大不同：前者頻域上能量分布較為廣泛，有較明顯能量的區(qū)間為0～200rad／s，主要集中在0～50rad／s范圍內(nèi)，時域上能量主要分布在35～80 s之間，其能量時－頻譜主軸主要平行于頻率軸，在較短時間內(nèi)包含了大量頻率信息；后者頻域上能量分布較前者集中，其中有較明顯能量的區(qū)間為0～20 rad／s，主要集中在0～5rad／s范圍內(nèi)，時域上能量分布較為廣泛，除起震部分時間段外，其他時間段均有明顯能量分布，其能量時－頻主軸主要平行于時間軸。究其原因，由于前者離震源機制較近，雖然綿竹清平臺處于自由場地之上，卻仍具有直接來自于震源機制且幅值較大的豐富頻率成分，時－頻譜表現(xiàn)為時域較窄而頻域?qū)拵У奶匦裕环粗靼才_處于渭河斷裂帶附近的第Ⅲ類遠(yuǎn)場地，由于在傳播過程中地震波能量不斷被吸收、不同頻率成分地震波的頻散效應(yīng)以及場地的過濾作用，致使其時－頻譜表現(xiàn)為時間寬帶而頻域較窄的特性。

4 結(jié) 論

本文首先對基于小波的非平穩(wěn)隨機過程演變功率譜密度估計作了若干討論。指出了小波系數(shù)的均方值與演變功率譜密度之間關(guān)系式（如式（8）或（10））的應(yīng)用范圍，闡述了此關(guān)系式實為Spanos－Failla方法的特殊形式，理清了MLP小波與GHW之間的關(guān)系以及它們在估計EPSD應(yīng)用之中的差別，因此進(jìn)一步澄清了基于小波估計非平穩(wěn)隨機過程EPSD的物理意義。研究表明，在利用一般小波估計功率譜時，不僅小波的尺度選擇至關(guān)重要，且非平穩(wěn)調(diào)制函數(shù)的慢變特性、小波的邊緣效應(yīng)對估計值也有較大影響。數(shù)值試驗表明，補零延長能較好地消除邊緣效應(yīng)帶來的影響。

根據(jù)一種最近提出的隨機過程的局部平穩(wěn)小波過程模型，推導(dǎo)了一種新的非平穩(wěn)隨機過程時變功率譜估計的方法。指出了所建議的新方法與基于Wold－Cramer隨機過程模型方法的關(guān)聯(lián)性：二者在形式上是統(tǒng)一的。當(dāng)基于LSW模型的新方法中的小波基為非抽樣廣義諧和小波時，此方法退化為Spanos－Tratskas關(guān)系式。因此，可以認(rèn)為Spanos－Tratskas方法不僅可用于 Wold－Cramer模型的非平穩(wěn)隨機過程，而且可應(yīng)用于任意的隨機過程。

為了驗證GHW在估計多維非平穩(wěn)隨機過程的可靠性，給出了同一地震中不同場地條件下，在考慮行波效應(yīng)和不同場地土?xí)r，目標(biāo)與估計自／互功率譜的對比算例。目標(biāo)瞬時均方值與估計瞬時均方值的對比，顯示此方法能滿足工程實際的精度要求。最后，以汶川地震中離斷層最近的綿竹清平臺以及西安臺所測得的地震加速度為例，用文中所建議方法對二者進(jìn)行了加速度能量的時－頻譜分析。顯示了二者在時－頻譜方面的明顯不同：前者表現(xiàn)為能量分布時域較窄而頻域?qū)拵В笳弑憩F(xiàn)為能量分布時間寬帶而頻域較窄的特性。

致謝：感謝國家留學(xué)基金委（CSC）對本文第一作者在美國Rice大學(xué)以聯(lián)合培養(yǎng)博士研究生進(jìn)行訪問期間給予的資助。

［1］ Trifunac M D．Response envelope spectrum and interpretation of strong earthquake ground motion ［J］．Bulletin of the Seismological Society of America，1971，61（2）：343—356．

［2］ Wang J，F(xiàn)an L，Qian S，et al．Simulations of non－stationary frequency content and its importance to seismic assessment of structures［J］．Earthquake Engineering＆Structural Dynamics，2002，31（4）：993—1 005．

［3］ Qian S．Introduction to Time－Frequency and Wavelet Transforms［M］．Pretice Hall，2001．

［4］ Qian S，Chen D．Joint Time－frequency Analysis：Methods and Applications［M］．New Jersey：Prentice Hall PTR，1996．

［5］ Gabor D．Theory of communication ［J］．Journal of the IEEE，1946，93（III）：429—457．

［6］ Wexler J，Raz S．Discrete Gabor Expansions［J］．Signal Processing，1990，21（3）：207—220．

［7］ Qian S，Chen D．Discrete Gabor Transform ［M］．New York，NY，ETATS－UNIS：Institute of Electrical and Electronics Engineers，1993．

［8］ Wigner E P．On the quantum correction for the thermodynamic equilibrium ［J］．Physics Review，1932，40：749—759．

［9］ Ville J．Theorie at applications de la notion de signal analytique［J］．Cables Transm，1948（2）：61—74．

［10］Newland D E．An Introduction to Random Vibrations，Spectral and Wavelet Analysis ［M］．New York：Longman Scientific ＆ Technical，1993．

［11］Daubechies I．Ten Lectures on Wavelets［M］．Philadelphia：Society for Industrial and Applied Mathematics，1992．

［12］Mallat S．Multiresolution approximation and wavelets［J］．Tansation of America Mathematics Society，1989，315：69—88．

［13］Grossmann A，Morlet J．Decomposition of Hardy function into square intergrable wavelets of constant shape［J］．SIAM J．Mathematics Anual，1984，15：723—736．

［14］Basu B．Wavelet－based stochastic seismic response of a duffing oscillator［J］．Journal of Sound and Vibration，2001，245（2）：251—260．

［15］Basu B，Gupta V K．Stochastic seismic response of single－degree－of－freedom systems through wavelets［J］．Engineering Structures，2000，22（12）：1 714—1 722．

［16］Basu B，Gupta V K．Seismic response of SDOF systems by wavelet modeling of nonstationary processes［J］．Journal of Engineering Mechanics，1998，124（10）：1 142—1 150．

［17］Iyama J，Kuwamura H．Application of wavelets to analysis and simulation of earthquake motions ［J］．Earthquake Engineering ＆ Structural Dynamics，1999，28（3）：255—272．

［18］Liang J－W，Chaudhuri S R，Shinozuka，M．Simulation of nonostationary stochastic process by spectral representation ［J］．Journal of Engineering Mechanics，2007，133（6）：616—627．

［19］Huang G，Chen X．Wavelets－based estimation of multivariate evolutionary spectra and its application to nonstationary downburst winds ［J］． Engineering Structures，2009，31（4）：976—989．

［20］Newland D E．Practical signal analysis：Do wavelets make any difference？［A］．Proceedings of the 16th ASME Biennial Conference on Vibration and Noise［C］．Sacramento，1997．

［21］Newland D E．Harmonic and musical wavelets［A］．Proceedings：Mathematical and Physical Sciences［C］．1994a：605—620．

［22］Newland D E．Harmonic wavelet analysis［A］．Proceedings of the Royal Society of London．Series A：Mathematical and Physical Sciences［C］．1993：203—225．

［23］Spanos P，Tezcan J，Tratskas P．Stochastic processes evolutionary spectrum estimation via harmonic wavelets［J］．Computer Methods in Applied Mechanics and Engineering，2005，194（12－16）：1 367—1 383．

［24］Tratskas P，Spanos P D．Linear multi－degree－of－freedomsystem stochastic response by using the harmonic wavelet transform ［J］．Journal of Applied Mechanics，2003，70（5）：724．

［25］Tratskas P． Wavelet－based excitation representation and response determination of linear and nonlinear systems［D］．Houston：Rice University，2001．

［26］Spanos P D，Kougioumtzoglou I A．Harmonic wavelet－based statistical linearlization for response evolutionary power spectrum determination．［J］．Probabilistic Engineering Mechanics，2011，doi：10．1016／j．probengmech．2011．05．008．

［27］Priestley M B．Spectral Analysis and Time Series［M］．Academic Press，1981．

［28］Priestley M B．Evolutionary spectra and non－stationary process［J］．Journal of the Royal Statistical Society，Series B，1965，27：204—237．

［29］Spanos P D，F(xiàn)ailla G．Evolutionary spectra estimation using wavelets［J］．Journal of Engineering Mechanics，2004，130（8）：952—960．

［30］Chakraborty A，Basu B．Nonstationary response analysis of long span bridges under spatially varying differential support motions using continuous wavelet transform ［J］．Journal of Engineering Mechanics，ASCE，2008，134（2）：155—162．

［31］Failla G，Pappatico M，Cundari G A．A wavelet－based spectrum for non－stationary processes［J］．Mechanics Research Communications，2011，38（5）：361—367．

［32］Eckley I A，Nason G P，Treloar R L．Locally stationary wavelet fields with application to the modeling and analysis of image texture ［J］．Applied Statistics，2010，59：595—616．

［33］Nason G P，von Sachs R，Kroisandt G．Wavelet processes and adaptive estimation of the evolutionary wavelet spectrum ［J］．Journal of the Royal Statistical Society：Series B（Statistical Methodology），2000，62（2）：271—295．

［34］Matz G，Hlawatsch F，Kozek W．Generalized evolutionary spectral analysis and the Weyl spectrum of nonstationary random processes ［M］．New York，NY，ETATS－UNIS：Institute of Electrical and Electronics Engineers，1997．

［35］Dahlhaus R．Fitting time series models to non－stationary processes［J］．The Annals of Statistics，1997，25：1—37．

［36］Ombao H，Raz J，Von Sachs R，et al．The SELX model of a non－stationary random process［J］．Annals of the Institute of Statistical Mathematics，2002，54：171—200．

［37］Rao M R，Bopardkar S A．Wavelet Transforms：Introduction to Theory and Applications ［M］．Massachusetts：Addison Wesley Longman，1998．

［38］Deodatis G．Non－stationary stochastic vector processes：seismic ground motion applications［J］．Probabilistic Engineering Mechanics，1996，11（3）：149—167．