時頻非平穩地震動模型及結構隨機振動應用*

何浩祥， 范少勇， 閆維明

(北京工業大學工程抗震與結構診治北京市重點實驗室 北京,100124)

引 言

地震動是由震源釋放的地震波引起的地面運動，具有復雜而強烈的時間-空間隨機性和非平穩性[1-2]。即使在同一次地震過程中，不同區域地震動的時頻特性也有所差別。對地震動時頻兩域的非平穩性進行深入研究能夠充分了解地震動的演變機制及隨機性，并可為地震動模擬提供重要的評價依據。此外，關于地震動特性的研究是土木工程抗震設計與動力分析的基礎，建立能夠全面反映地震動非平穩性的模型可以更準確地計算結構隨機響應，為結構精細化隨機振動分析和可靠性評價提供堅實的基礎。

由于通常無法在一次地震中同一區域取得多次地震動記錄，傳統的隨機理論通常假定每一條地震動記錄均符合各態歷經假定，即認為可用時間的平均來代替集合的平均[2-3]，這導致具體地震動的確定性特征和隨機成分沒有被準確分離，且地震動時頻域非平穩性被顯著弱化。如果隨之進行的結構隨機振動分析對非平穩性考慮不充分，則結構隨機響應將更不準確。在傳統的隨機振動研究中通常采用均勻調制方法反映地震動的時域非平穩性，即將真實地震動等效為一個零均值的平穩隨機過程(一般為白噪聲)和一個隨時間漸變的確定性強度包絡函數的Hadamard乘積(各元素點乘)[1-3]。常用的強度包絡函數主要有連續型和分段型兩類[2,4]。采用均勻調制方法得到的地震動時程包絡一般均為單峰型曲線，并不能準確刻畫地震動幅值變化的豐富性和復雜性，其頻域分布具有明顯的平穩性，不能反映地震動頻域的強非平穩特征和能量分布特點。對于多階自振周期與地震動卓越周期接近的結構，采用均勻調制地震動進行隨機振動分析有可能嚴重影響結果的準確性。

為了解決上述不足，Saragoni等[5]提出通過時間段內的Gamma函數調制過濾高斯白噪聲的方法，并提出了一種時域和頻域完全非平穩隨機地震動模型(包含強度和頻率非平穩)。Lin等[6]提出用散粒噪聲、過濾白噪聲或三角級數疊加信號等替換均勻調制方法中的白噪聲。Sanaz[7]提出采用強度包絡函數和時變過濾白噪聲模型共同模擬地震動時頻域非平穩特性。在具體應用以上方法時存在難以精確確定頻帶能量分布的局限，雖然在人工地震動生成中易于實現，而在真實地震動隨機分析中難以應用。梁建文[8]在演化譜理論的基礎上，采用一個譜表示方法來模擬頻率非平穩地震動過程，但此方法依賴目標譜和權重函數。Conte等[9]構建了基于sigma振子過程的完全非平穩模型，并可由真實地震動演變功率譜經自適應最小二乘法擬合確定模型參數，雖然其合成精度相對較高但計算復雜，且對時頻兩域的非平穩細節模擬仍不充分。此外，部分研究者建議在計算結構響應的統計特征時直接將相應的輸入白噪聲功率譜改為傳統的地震動功率譜工程學模型[10]，如Kanai-Tajimi模型及其系列修正模型等，但由于相關的功率譜模型對頻率非平穩性的表征依然不夠精細[11-12]，且該處理方法在理論上也不嚴謹，因此不宜作為有效方法。

因此，針對上述問題開展新型時頻域非平穩地震動模型和相關隨機振動分析方法具有重要的理論與工程意義。筆者從地震動的傳播和演變機制以及隨機特性上進行探究，提出能夠更準確全面反映地震動時頻域非平穩特性的強度包絡模型，并改進傳統的均勻調制方法，進而建立反映時頻非平穩關聯特征的結構隨機振動分析方法。

1 非均勻調制地震動模型

描述地震動的物理和數學模型眾多[4,13]，從應用性的角度可分為地震動工程應用模型和地震動隨機過程模型兩種，前者強調地震動模型及其合成信號的特性應滿足指定的反應譜或功率譜的需求，而后者更注重于地震動時頻非平穩性的描述，并主要應用于結構隨機振動分析領域。筆者的研究主要針對地震動隨機過程模型開展。

在傳統的均勻調制地震動隨機模型中，并不考慮頻率非平穩性，并通常將地震動表示為

A(t)=G(t)n(t)

(1)

其中：A(t)為非平穩地面加速度過程；G(t)為強度包絡函數；n(t)為一零均值的平穩隨機過程。

在上述模型中，強度包絡函數通常是單峰正值函數，無法準確全面地刻畫地震動多次波動的多峰現象，也不能體現地震動正負幅值的差異，因而不能充分反映地震動在時域上的隨機性和非平穩性。更重要的是，強度包絡函數并不能明顯改變平穩隨機過程在頻域上的能量分布，因此上述模型不能反映頻域非平穩特性，更無法直接轉化為有效的速度和位移時程。增強上述模型對時頻域非平穩性的表征能力需要同時對G(t)和n(t)的內涵和形式進行改進，相關的改進可以從地震動的傳播和演變機制的研究成果中得到借鑒。

在地震發生時，震源所產生的尖脈沖經地層介質的吸收后，其波形將會被拉長，這種具有時延的單個反射波稱為地震子波。在地震學和地震勘探中，褶積模型理論認為地震波形是地震子波動態疊加的結果[14-15]。相關研究表明：從地下許多反射界面發生反射時形成的地震子波，其振幅取決于反射界面反射系數的絕對值，極性的正負決定于反射系數的正負, 到達時間的先后取決于界面深度和覆蓋層的波速。實際的地震記錄A(t)可視為由地震子波w(t)和反射系數γ(t)褶積組成，相關公式如下

A(t)=w(t)γ(t)

(2)

比較式(1)和式(2)可發現，二者的數學形式一致但物理意義不同。在褶積模型中，地震子波是振蕩衰減型低頻隨機信號，具有多峰性，幅值也呈現正負交替，這與均勻調制模型中的強度包絡函數有明顯區別。此外，褶積模型中的反射系數是零均值頻域非平穩噪聲信號，這與均勻調制模型中的白噪聲也有所區別。地震子波和反射系數的非平穩特性及其褶積特征能較充分反映地震動的時頻域非平穩性，這為均勻調制模型的改進提供了良好的借鑒。

在諸如褶積模型的地震動分解處理中，地震子波和反射系數依然是復雜隨機信號，可以采用信號處理技術對其進行更深入的分析，而濾波方法是判斷信號頻域非平穩特征的主要手段之一。大量地震動頻譜分析結果表明，地震動的頻域能量主要分布在0～10Hz。由于地震動信號采樣率較高，頻譜分析中最高頻率通常超過25Hz，因此可以通過濾波實現地震信號中重要成分和次要成分的分離。然而，采用諸如低通濾波器的理想濾波技術對信號的處理過于嚴格和刻板，缺乏靈活性和調控能力。相反，采用移動平均算法對信號進行平滑處理，相當于對信號進行了低通濾波，且相應的濾波器具有從通帶到阻帶平緩過渡的特點，因此既能分離局部隨機波動的特征又能凸顯整體變化規律。文獻[16]提出地震動時域信號可分解為多峰平滑信號和非平穩噪聲信號兩部分。多峰平滑信號主要反映了地震產生機理和演變規律的基本特征，雖然由于知識不完備性難以闡釋其細節變化，而非平穩噪聲信號則反映了地震動傳播過程中受到的擾動和強隨機性。在此基礎上，文獻[16]建立由多峰平滑時域信號As(t)和非高斯時域噪聲信號An(t)構成的地震動時域模型A(t)

(3)

非高斯時域噪聲信號An(t)具有如下表達形式

An(t)=As(t)|αNn(t)|

(4)

其中：Nn(t)為歸一化的高斯白噪聲；α為其調幅系數。

由式(3)和式(4)可得到地震動時域模型的綜合表達式為

(5)

對比式(2)和式(5)可發現，地震動時域模型A(t)與褶積模型在表達形式和物理內涵上均類似，因此可以改進均勻調制模型的不足，為精確細致反映地震動非平穩特性提供了必要的處理方法和分析工具。

上述地震動時域模型是由信號能量平衡條件來約束，雖然精度較高但在結構隨機振動分析中不易處理，因此筆者建議將式(3)和式(5)改為如下的非均勻調制形式

(6)

其中：Acn(t)表示均值為1的白噪聲。

雖然式(6)相當于忽略了原地震動時域模型A(t)的能量表達式中的交叉項，但由于As(t)和An(t)均為隨機性較強的信號，其乘積相對較小，且可通過再次的平滑處理和調整調幅系數α來實現進一步的修正，因此依然能較準確地反映原有地震動信號的隨機特性。

此外，多峰平滑時域信號As(t)可以通過以下兩式計算得到

其中：A0(t)為原始地震動信號；As0(t)為初步平滑信號；S函數為平滑函數，在Matlab程序中可采用smooth命令實現；n和m均為滑動窗寬度或需要考慮的數據點數；sgn為符號函數。

由于As(t)需要反映地震動本質特性，因此試算時可不斷增加滑動窗寬直到其幅值穩定為止。在具體分析中，n和m建議取值區間分別為[150,200]和[5,30]，同時需要選取適當的調幅系數α,使As(t)的最大峰值與原波A(t)最大峰值相同。

為了驗證筆者提出的非均勻調制地震動模型的精確性，應用相關方法對Taft(N21E)地震加速度信號進行分解和合成，并與傳統的均勻調制地震動模型進行比較，相應的結果如圖1所示。對原信號進行平滑處理時，式(7)和式(8)滑動窗寬度n和m分別取為200和5，優化后的調幅系數α為2.25。此外，由于已經知道Taft波的具體信息，因此在均勻調制合成方法中沒有采用常用的強度包絡函數模型，而是通過直接對地震波絕對值進行窗寬為500的平滑處理并歸一化獲得。

圖1 不同地震動模型加速度合成效果對比

由圖1可看出，多峰平滑包絡信號As(t)具有往復振蕩的特點，且最大值不必等于1，這與傳統的時域強度包絡信號有明顯的區別。此外，噪聲信號Acn(t)的幅值相對較大。采用非均勻調制模型合成的地震動與真實地震動在時域上十分接近，而按照均勻調制模型合成的地震動則與真實值有明顯差別。為了對比不同模型在頻域上的精度，對不同信號進行了功率譜分析，如圖2和圖3所示。

圖2 不同地震動模型的加速度功率譜對比

圖3 非均勻調制模型各部分頻譜分布

結果表明，由非均勻調制模型合成的地震動的頻譜分布也與真實地震動接近，而均勻調制模型的頻譜分布較均勻且不能反映卓越頻率特征。從圖3也可看出，多峰平滑包絡信號充分保留了真實地震動的主要頻譜成分，是非均勻調制地震動模型精確性的根本保障。為進一步驗證非均勻調制地震動在低頻信號合成方面的能力和精度，對Taft(N21E)地震位移信號進行合成，滑動窗寬度n和m分別取為150和26，優化后的調幅系數α為4.0，時域和頻域結果如圖4所示。由圖可見，非均勻調制模型依然具有良好的精度，而均勻調制模型對低頻能量的表征精度進一步下降，該結論與前人研究結果吻合[2]。

綜上所述，由于平滑處理后的多峰平滑包絡信號具有正負振蕩特性以及足夠的細節刻畫能力，地震動的主要頻譜成分被充分提取，加之噪聲信號的強隨機性，由二者合成的地震信號能夠精確表征原信號在時域和頻域的非平穩特征。因此，非均勻調制地震動模型既可以根據原波產生大量的地震動隨機信號，從而為基于Monte Carlo方法的結構隨機分析提供可靠的隨機樣本，也可以為基于功率譜的結構時域隨機振動分析提供更精確的激勵模式。

圖4 地震動位移信號合成效果對比

2 全非平穩地震動模型

目前計算演變譜的方法雖然較豐富，但精度通常不高。小波包分解方法采用正交小波將信號分解成各尺度上的正交小波包分量，再對各分量用其相應的小波包基函數為窗函數進行時頻變換，其分解機理明確，由此得到的時變功率譜從理論表達上是準確的，且分辨率高。文獻[16]的研究表明，將小波包分解產生的Paley序子頻帶序列按升頻進行重新排列，并采用正交性和緊支撐性的離散Meyer(dmey)小波作為小波包基函數可以更準確地描述模擬地震記錄信號時頻變化特征。

基于上述研究并參考非均勻調制地震模型的構建模式，筆者提出采用二維平滑技術的全非平穩地震動演變譜模型A(ω,t)，該模型由多峰平滑時頻域譜As(t)和時頻域非高斯噪聲譜An(t)構成，表達式為

A(ω,t)=As(ω,t)+An(ω,t)=

As(ω,t)[1+αNn(t)]=As(ω,t)Acn(t)

(9)

其中：As(t)和Acn(t)的維數均為p×q，p為頻率向量的維數，q為時域向量的維數。

為了驗證上述全非平穩地震動演變譜模型的精確性，采用離散Meyer小波函數對Taft(N21E)地震加速度信號進行9層小波包分解，從而獲得地震動演變譜，如圖5(a)所示；利用穩健補償最小二乘計算方法[19]對該演變譜進行二維平滑處理，總滑動窗寬度為2，獲得的多峰平滑演變譜如圖5(b)所示；生成調幅系數α為2.50的隨機噪聲譜，如圖5(c)所示；最后按式(9)合成地震動演變譜，如圖5(d)所示。為了便于觀察，上述圖采用譜絕對值表示原波頻譜和多峰平滑演變譜等結果。

從結果可看出：平滑演變譜與原波演變譜在總體變化規律上類似，但幅值較小，合成后的演變譜與原波演變譜十分接近，能同時反映時頻域的非平穩變化特征。分別計算真實演變譜和合成演變譜的時域和頻域的邊緣條件，即分別求演變譜在時域和頻域的疊加值，結果如圖6和圖7所示。由圖可知，全非平穩合成地震動的頻域合成精度較高，而時域合成精度略低，主要原因在于目前的二維平滑算法不易實現時頻兩域同精度的協調平滑，該問題可以通過改進算法得到進一步解決。

圖5 全非平穩地震譜分解與合成

圖6 全非平穩地震動時域信號

此外，還可基于小波包分解分別獲得由均勻調制模型和非均勻調制模型生成的地震波的演變譜，如圖8所示。結果表明，非均勻調制模型地震波演變譜與原波演變譜吻合較好，其頻域上的非平穩性也較準確細致。均勻調制模型地震波的時頻分布與真實波有較大差異，在地震動模擬和結構分析時需謹慎使用。為了進一步比較以上不同地震動模型的合成精度，筆者選取了Ⅰ～Ⅳ場地的10條地震波進行相關分析，并以合成偏差率作為精度評判指標。時域合成偏差率定義為原波與合成波各點幅值平方差的累積和與原波幅值平方和的比值；頻域合成偏差率定義為原波與合成波的功率譜幅值差的累積和與原波功率譜的比值。均勻調制模型地震動在時域和頻域上的平均合成偏差率分別為0.292和0.304；非均勻調制地震動的平均合成偏差率分別為0.076和0.069；全非平穩地震動的平均合成偏差率分別為0.231和0.091?？梢?，非均勻調制地震動模型和全非平穩地震動模型的精度均比傳統的均勻調制模型高，在地震動合成或結構隨機振動分析中采用這兩種模型將有效提升計算精度。

圖7 全非平穩地震動功率譜

圖8 均勻與非均勻調制模型時頻譜對比

3 基于全非平穩譜模型隨機振動分析

由于筆者提出的非均勻調制模型和全非平穩模型的形式已與均勻調制模型有明顯差別，因此需要在傳統隨機振動分析理論的基礎上進行適當推廣。

若采用非平穩隨機激勵f(t)表示為前述多峰平滑時域信號和均值為1的白噪聲的點乘形式，則在該改進均勻調制非平穩隨機激勵f(t)的作用下，初始靜止的線性體系響應可表示為

(10)

其中：h為脈沖響應函數；τ為0～t的任一時刻。

由于多峰平滑包絡的均值為零，故激勵和體系響應均具有零均值，且其協方差函數為

Ky(t1,t2)=E[y(t1)y(t2)]=

Ry(t1,t2)-E[y(t1)]E[y(t2)]=Ry(t1,t2)

(11)

響應y(t)的協方差函數和自相關函數Ry(t)相等，且有

(12)

根據Wiener-Khinchin定理，對于寬平穩隨機過程Acn(t)有

E[Acn(τ1)Acn(τ2)]=Rx(τ1-τ2)=

(13)

其中：Scn為非均勻調制模型中噪聲Acn(t)的自功率譜密度函數。

因此，式(11)可表示為

(14)

其中

(15)

令t1=t2=t，則可得到響應y(t)的方差為

(16)

其中：σy為響應的標準差。

響應y(t)的自演變譜密度為

Sy(ω,t)=I*(ω,t)I(ω,t)Scn(ω)=

(17)

對于多自由度線性結構，在非均勻調制地震激勵A(t)=As(t)Acn(t)的作用下，其隨機動力方程為

(18)

令Y(t)=Φq(t)，其中Φ為結構振型向量，將其帶入上式并在方程兩側同乘ΦT，則有

-γ*M*As(t)Acn(t)

(19)

其中：M*=ΦTMΦ;C*=ΦTCΦ;K*=ΦTKΦ;γ*=ΦTMI/M*。

考慮到結構參數矩陣與振型的正交性，上式可表示為n個單自由度方程

(20)

該系統的響應Y(t)的方差為

(21)

若由于結構自由度數較多或計算效率偏低，可考慮采用虛擬激勵法計算結構隨機響應[20-21]?；谔摂M激勵法構造一個虛擬的外部激勵，則式(18)變為

(22)

仍可采用振型分解法得到該系統的響應Y(t)的方差，有

(23)

(24)

上述分析方法是針對非均勻調制激勵的，當采用完全非平穩地震模型時，需要按式(9)計算全非平穩地震動演變譜模型A(ω,t)，并按下式計算Ii

(25)

其他計算公式類同。由于本研究提出的非平穩地震動隨機模型的精度主要取決于多峰平滑信號的時頻域特性，加之實際結構的動力分析需要考慮地震動傳播條件和場地特征的影響，因此對于具體結構的隨機分析應選取相應場地的典型地震動或滿足工程要求的人工生成地震波進行計算，并確保其多峰平滑信號具有代表性。

4 算例與分析

為了進一步驗證非平穩地震動模型及其在結構隨機振動分析中的有效性，筆者建立了層數分別為1～20的20個二維剪切型鋼筋混凝土框架結構并進行隨機振動分析。為便于比較不同結構的動力響應規律，各個結構的層高均為3 m，層重均為5.2×108kg。假定各個結構的無阻尼基本周期為0.06n(n為相應總層數)，則可算出各層的等效剪切剛度，進而建立剛度矩陣。假定各結構阻尼比均為5%，則可構建Rayleigh阻尼矩陣，最終按式(16)建立結構動力學方程。將Taft波作為水平向地震動輸入，按照隨機振動方法計算結構隨機地震響應。為了研究地震動合成模型對結構響應的影響，筆者基于均勻調制模型、非均勻調制模型和全非平穩模型分別對以上20個具有不同自振周期結構的標準差進行計算分析?；诰鶆蛘{制模型和非均勻調制模型的5層結構的1層、3層和5層的位移標準差如圖9所示，其他典型結構頂層位移標準差結果依次如圖10所示，5層結構的頂層位移演變譜絕對值如圖11所示。

圖9 結構典型樓層位移標準差

圖10 典型結構頂層位移標準差比較

圖11 結構的頂層位移演變譜絕對值

從以上結果可以看出，基于非均勻調制模型和全非平穩模型的結果較接近，且其波動規律與地震動幅值的變化規律相關，表明其隨機響應能夠較細致地反映地震動時域非平穩性對結構的影響。此外，當基本周期較小時基于上述兩種模型的結構標準差隨時間變化劇烈，隨著周期增大其變化趨緩，反映了地震動頻域非平穩性的影響。與此相反，由于均勻調制模型反映時頻域非平穩性的能力均明顯不足，相應的結果通常變化平緩且最大值偏小，不能充分反映地震動非平穩性對結構的影響。

為了進一步對比3種方法結果的差異并探究其原因，提取基于均勻調制模型和全非平穩模型的所有結構頂層位移標準差最大值，并與基于非均勻調制模型的結果相比獲得最大值的比值，根據結構基本頻率(基本周期倒數)列出上述結果，如圖12所示。

圖12 結構標準差最大值比值規律

可以發現，均勻調制模型最大值均小于非均勻調制模型最大值，且其變化規律與原地震動的功率譜外形相近。這是因為均勻調制模型地震動功率譜的頻率能量分布是均勻的，而非均勻調制模型地震動的頻譜特征能夠充分反映地震動頻域非平穩性，因而與真實功率譜吻合，而二者的差異仍然具有真實功率譜的特征，最終導致按照均勻調制模型計算短周期結構(基本周期小于0.3s)時很可能嚴重低估結構響應。此外，如果地震動的長周期成分較顯著或結構屬于大型復雜結構時，由均勻調制模型得到隨機響應也可能偏低，不利于結構的安全設計。由全非平穩模型得到的結構響應與非均勻調制模型結果較接近，但由于其在時頻兩域的非平穩特征相對不夠精確，因此建議作為參考方法。

綜上，可認為按照傳統均勻調制模型計算結構隨機響應不能充分獲得結構時頻域的非平穩特性，結果偏于不安全，基于筆者提出的非均勻調制模型的結構隨機響應計算方法能夠較好地解決上述問題?；谌瞧椒€模型的結構響應計算方法在原理上更嚴格而全面，但由于現有多維平滑技術的協調性和精度略低，因此需要進行更深入的研究和改進。

5 結 論

1) 地震動在時域和頻域均具有復雜的隨機性和非平穩性，傳統的地震動隨機過程模型通常采用均勻調制模型，并不能全面細致地表征地震動時頻非平穩性。筆者根據地震褶積模型理論以及之前的研究成果，提出由多峰平滑信號和非平穩噪聲信號點乘而成的非均勻調制地震動模型。該模型能夠在時域上充分反映地震波的波動和衰減特性，其頻譜分布也與真實地震動十分接近，因此可以同時表征時頻域非平穩性。在此基礎上，根據演變譜原理和需求提出了以時間和頻率為變量的全非平穩地震動建模方法。分析結果表明，筆者提出的兩種模型在時域和頻域的合成精度上均優于均勻調制模型。

2) 針對非均勻調制模型和全非平穩模型的形式和特點，在傳統隨機振動分析理論的基礎上對結構隨機振動響應求解方法進行了改進。通過算例證明了按照均勻調制模型計算結構隨機響應不能充分獲得結構時頻域的非平穩特性，且結果偏于不安全，而利用筆者提出的兩種模型可以獲得更精確的結構響應概率特征和更細致的非平穩性演變過程。本研究主要針對具體地震動特性開展，在今后研究中可根據不同場地條件和工程需求獲取具有典型非平穩特征和概率統計特征的地震波形，實現結構隨機響應特征的分析和預測。