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

孔 凡，李 杰，2

（1．同濟(jì)大學(xué)土木工程學(xué)院，上海 200092；2．同濟(jì)大學(xué)土木工程防災(zāi)國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海 200092）

引 言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 數(shù)值算例

3．1 模擬多變量地震動(dòng)自／互功率譜估計(jì)

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

圖6 同一地震中的3個(gè)不同場(chǎng)地條件Fig．6 Three different site conditions in a same earthquake event

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

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

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

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

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

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

圖9 軟粘性土場(chǎng)地上歸一化估計(jì)均方值與目標(biāo)均方值的對(duì)比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í)地震加速度時(shí)程的能量時(shí)間－頻率分布

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

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

圖10 汶川地震中不同場(chǎng)地地震動(dòng)加速度記錄Fig．10 Acceleration records of ground motion in Wenchuan earthquake

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

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

4 結(jié) 論

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

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

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

致謝：感謝國(guó)家留學(xué)基金委（CSC）對(duì)本文第一作者在美國(guó)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．