基于線性遞推濾波算法的地震反應(yīng)譜計(jì)算方法

徐玉聰

1 長(zhǎng)江三峽勘測(cè)研究院有限公司(武漢), 武漢市創(chuàng)業(yè)街99號(hào)，430074

地震反應(yīng)譜是地震工程研究中的一項(xiàng)重要內(nèi)容，通過(guò)研究單自由度震動(dòng)體系下質(zhì)點(diǎn)的運(yùn)動(dòng)情況來(lái)描述地震條件下結(jié)構(gòu)的最大動(dòng)力反應(yīng)[1-2]和地震動(dòng)的頻譜特性。反應(yīng)譜在強(qiáng)震監(jiān)測(cè)計(jì)算分析中既可反映局部場(chǎng)地的地震強(qiáng)度特征，又可用來(lái)估算地震作用下工程結(jié)構(gòu)的受力情況，大量的反應(yīng)譜計(jì)算結(jié)果還可應(yīng)用于結(jié)構(gòu)抗震設(shè)計(jì)。因此，采用有效的方法提高反應(yīng)譜的精度和計(jì)算效率具有重要意義。

反應(yīng)譜的計(jì)算方法主要集中在時(shí)域和頻域兩個(gè)維度。頻域計(jì)算方法的基本原理是在頻域計(jì)算地震動(dòng)和傳遞函數(shù)傅里葉頻譜的乘積[3]，然后通過(guò)傅里葉逆變換得到時(shí)間域的反應(yīng)譜，需要已知頻域內(nèi)精確的傳遞函數(shù)[4]以及完整的地震記錄來(lái)完成計(jì)算，其中長(zhǎng)周期反應(yīng)譜的計(jì)算精度受地震記錄長(zhǎng)度影響較大，因此頻域不適合進(jìn)行實(shí)時(shí)快速反應(yīng)譜計(jì)算。反應(yīng)譜時(shí)域計(jì)算方法一直為研究重點(diǎn)，其基本原理是以地震反應(yīng)量所表述的動(dòng)力平衡方程來(lái)尋找數(shù)字濾波器傳遞函數(shù)或遞歸公式[5]。時(shí)域計(jì)算方法有Duhamel逐步積分法、Newmark-β方法[6]、精確解法(ASM)[7]和Z變換方法[8]等，這些方法的計(jì)算精度主要受假定時(shí)間步長(zhǎng)內(nèi)加速度的特定分布形式所制約[9]。

本文基于單自由度體系運(yùn)動(dòng)方程推導(dǎo)計(jì)算地震反應(yīng)譜的一般表達(dá)形式，并從傳遞函數(shù)的形式和濾波系數(shù)入手，推導(dǎo)用于計(jì)算反應(yīng)譜的遞推表達(dá)式，然后結(jié)合已有研究提出的滿足高頻和低頻收斂條件對(duì)應(yīng)的獨(dú)立積分常數(shù)即可計(jì)算得到反應(yīng)譜結(jié)果。

1 單自由度震動(dòng)體系原理

當(dāng)?shù)孛嬉砸欢铀俣冗\(yùn)動(dòng)時(shí)，單自由度震動(dòng)體系的運(yùn)動(dòng)方程為[2]：

(1)

定義時(shí)間加權(quán)函數(shù)W(t)為：

W(t)=1+w1Γ+w2Γ2+w3Γ3

(2)

(3)

(4)

(5)

(6)

式中，Λ1、Λ2、Λ3、Λ4、Λ5、Λ6為一系列自由參數(shù)。

將式(1)代入式(2)可得：

(7)

(8)

(9)

(10)

將式(9)和式(10)代入式(8)可得：

dn+1=Adn+Ln

(11)

A和Ln中變量可整理為：

A31=Ω2/D，A32=(2ξΩ+Λ1W1Ω2)/D

A33=1-(1+2Λ4W1ξΩ+Λ2W2Ω2)/D

D=Λ6W1+2Λ5W2ξΩ+Λ3W3Ω2，Ω=ωΔt

2 計(jì)算原理

在Z變換域中輸入、輸出和傳遞函數(shù)的形式為：

Y(Z)=

(12)

式中，na為反饋濾波器階數(shù)；nb為前饋濾波器階數(shù)?？蓪⑹?12)表示為差分方程：

y(m)=B1x(m)+B2x(m-1)+…+

Bnb+1x(m-nb)-A1y(m-1)+

…+Bnay(m-na)

(13)

根據(jù)Zhou等[10]提出的思路，構(gòu)造3階反饋濾波器系數(shù)和4階前饋濾波器系數(shù)，令[1,-A1,A2,-A3]、[B1,B2,B3,B4]為2組濾波器系數(shù)，其中A1、A2、A3、B1、B2、B3、B4分別由式(11)中矩陣A中元素得到：

B1=Δt2λ3W1/D

B2=Δt2[λ3(1-W1)-(A22+A33)λ3W1+

A12λ5W1+A13W1]/D

B3=Δt2λ[-(A22+A33)(1-W1)λ3+

A12λ5(1-W1)+A13(1-W1)+

(A22A33-A23A32)λ3W1-(A12A33-

A13A32)λ5W1+(A12A23-A13A22)W1]/D

B4=Δt2λ[(A22A33-A23A32)(1-W1)λ3-

(A12A33-A13A32)λ5-(1-W1)+

(A12A23-A13A22)(1-W1)]/D

y(m)=B1x(m)+z1(m-1)

(14)

z1(m)=B2x(m)+z2(m-1)-A1y(m)

(15)

z2(m)=B3x(m)+z3(m-1)-A2y(m)

(16)

z3(m)=B4x(m)-A3y(m)

(17)

3 合成數(shù)據(jù)分析計(jì)算

3.1 計(jì)算精度分析

為驗(yàn)證本文方法的計(jì)算精度，采用由2個(gè)頻率合成的正弦簡(jiǎn)諧波f(t)作為單自由度震動(dòng)體系的系統(tǒng)輸入，如式(18)和式(19)所示，其中采樣頻率為200 Hz，時(shí)長(zhǎng)為30 s，f1=0.5 Hz，f2=2.0 Hz。設(shè)置反應(yīng)譜阻尼比為5%，周期為0～4 s，時(shí)間間隔為0.02 s，分別利用LRFM方法、Duhamel逐步積分法、Newmark-β方法和精確解法(ASM)[12]計(jì)算得到反應(yīng)譜結(jié)果，并以ASM方法作為對(duì)比驗(yàn)證另外3種方法的相對(duì)誤差：

f(t)=sin(2π×f1×t)+sin(2π×f2×t)，

t∈[8,15]

(18)

f(t)=0，t∈[0,8)或t∈(15,30)

(19)

采用4種方法計(jì)算得到的位移反應(yīng)譜如圖1(a)所示，從圖中可以看出，4種方法計(jì)算得到的曲線形態(tài)一致，均反映出輸入周期在0.5 s和2 s處的模擬地震動(dòng)信息。位移反應(yīng)譜的相對(duì)誤差(圖1(b))結(jié)果顯示，LRFM方法、Duhamel逐步積分法和Newmark-β方法在周期0～1 s內(nèi)均存在一定誤差，在周期2～4 s內(nèi)誤差較小，在0.5 s處具有明顯突變；LRFM方法在周期0～1 s內(nèi)的相對(duì)誤差結(jié)果優(yōu)于Duhamel逐步積分法和Newmark-β方法，且整體誤差波動(dòng)較小。

圖1 位移反應(yīng)譜計(jì)算結(jié)果及誤差

利用合成數(shù)據(jù)計(jì)算得到的速度反應(yīng)譜和誤差結(jié)果如圖2所示，從圖中可以看出，4種方法計(jì)算得到的速度反應(yīng)譜曲線幾乎完全重合，其中Duhamel逐步積分法的相對(duì)誤差在短周期0～1 s內(nèi)跳動(dòng)較大，最大達(dá)到18%；而Newmark-β方法和LRFM方法的整體誤差較小，基本可達(dá)到ASM方法的計(jì)算效果。

圖2 速度反應(yīng)譜計(jì)算結(jié)果及誤差

圖3為4種方法計(jì)算得到的加速度反應(yīng)譜和誤差結(jié)果，從圖中可以看出，4種方法計(jì)算的加速度反應(yīng)譜曲線僅在個(gè)別峰值點(diǎn)存在突變，而LRFM方法、Duhamel逐步積分法及Newmark-β方法的相對(duì)誤差均小于2%，3種方法的相對(duì)誤差均集中在短周期部分。LRFM方法的相對(duì)誤差基本呈線性正相關(guān)變化，而Duhamel逐步積分法和Newmark-β方法的相對(duì)誤差在正負(fù)方向來(lái)回震蕩，且在0.5 s處存在突變。

圖3 加速度反應(yīng)譜計(jì)算結(jié)果及誤差

通過(guò)對(duì)合成數(shù)據(jù)的計(jì)算結(jié)果進(jìn)行分析可知，LRFM方法的計(jì)算精度和穩(wěn)定性優(yōu)于Duhamel逐步積分法和Newmark-β方法，在位移最大值處的誤差較小，而加速度反應(yīng)譜的計(jì)算誤差整體略大于另外2種方法，但相對(duì)誤差小于2%，計(jì)算結(jié)果滿足精度要求。

3.2 計(jì)算效率分析

在計(jì)算地震反應(yīng)譜時(shí)，為獲得更充分的反應(yīng)譜信息，一般選取輸入時(shí)長(zhǎng)較長(zhǎng)的強(qiáng)震動(dòng)記錄，特別是在計(jì)算強(qiáng)震動(dòng)條件下的加速度反應(yīng)譜時(shí)，強(qiáng)震動(dòng)記錄時(shí)長(zhǎng)可達(dá)300 s以上，因此會(huì)產(chǎn)生大量強(qiáng)震動(dòng)記錄，如果反應(yīng)譜計(jì)算耗時(shí)太長(zhǎng)，將不利于反應(yīng)譜的進(jìn)一步分析。為對(duì)比分析LRFM方法、Duhamel逐步積分法和Newmark-β方法的計(jì)算效率，利用式(18)作為系統(tǒng)輸入，合成數(shù)據(jù)計(jì)算時(shí)長(zhǎng)分別設(shè)置為30 s、60 s、90 s、180 s和300 s，采樣頻率為200 Hz，各方法的計(jì)算耗時(shí)如圖4所示。由圖可知，Newmark-β方法的計(jì)算耗時(shí)明顯大于另外2種方法，且隨著時(shí)長(zhǎng)的增加，計(jì)算耗時(shí)呈現(xiàn)線性變化，在300 s處計(jì)算耗時(shí)接近10 s；LRFM方法和Duhamel逐步積分法的整體耗時(shí)在1 s以內(nèi)，兩者的計(jì)算效率基本相當(dāng)，且顯著優(yōu)于Newmark-β方法。

圖4 反應(yīng)譜計(jì)算耗時(shí)對(duì)比

隨著我國(guó)強(qiáng)震動(dòng)監(jiān)測(cè)臺(tái)站的增加，強(qiáng)震及大地震發(fā)生后可獲取數(shù)量多且持續(xù)時(shí)間長(zhǎng)的強(qiáng)震動(dòng)記錄，如2008年汶川地震發(fā)生后，僅四川省數(shù)字強(qiáng)震臺(tái)網(wǎng)獲得的主震記錄就接近400條，大部分紀(jì)錄的持續(xù)時(shí)間為150～160 s[13]，最長(zhǎng)超過(guò)300 s，因此選取合適的方法快速計(jì)算得到反應(yīng)譜以評(píng)估場(chǎng)地的地震強(qiáng)度特征或工程結(jié)構(gòu)的受力情況尤為重要。本文提出的LRFM方法的計(jì)算效率顯著優(yōu)于目前普遍采用的Newmark-β方法，可快速高效地得到計(jì)算結(jié)果。

4 實(shí)際數(shù)據(jù)分析計(jì)算

為驗(yàn)證LRFM方法計(jì)算實(shí)際強(qiáng)震動(dòng)記錄反應(yīng)譜的效果，選取ESM強(qiáng)震數(shù)據(jù)庫(kù)收錄的SNO、PCB、FOC、TLN、AMT、MSCT、CLF、FOS共8個(gè)臺(tái)站記錄的2016-08-24意大利ML6.0地震強(qiáng)震波形數(shù)據(jù)，具體參數(shù)見(jiàn)表1，采樣頻率均為200 Hz。選取的強(qiáng)震動(dòng)記錄震中距為8.5～57.50 km，PGA為47.13～850.80 Gal，卓越周期范圍為0.06～0.56 s，這些記錄可基本反映該地震的強(qiáng)震動(dòng)特征，計(jì)算得到的反應(yīng)譜能夠代表在該地震作用下不同位置的地面運(yùn)動(dòng)，可以很好地驗(yàn)證LRFM方法對(duì)實(shí)際強(qiáng)震動(dòng)記錄反應(yīng)譜的計(jì)算效果。計(jì)算反應(yīng)譜的阻尼系數(shù)選取5%、10%、20%，自振周期為0～4 s，時(shí)間間隔為0.02 s，選取的強(qiáng)震記錄時(shí)程曲線如圖5所示。本文選取的強(qiáng)震記錄具有一般性，主頻可基本覆蓋從低頻到高頻不同范圍的信息。

圖5 地震動(dòng)加速度時(shí)程

表1 強(qiáng)震動(dòng)記錄參數(shù)

以ASM方法的計(jì)算結(jié)果作為參照，驗(yàn)證LRFM方法計(jì)算實(shí)際強(qiáng)震動(dòng)數(shù)據(jù)反應(yīng)譜的效果，另外由于不同方法計(jì)算得到的位移反應(yīng)譜誤差相差較小，因此本文僅分析速度反應(yīng)譜和加速度反應(yīng)譜的計(jì)算結(jié)果。圖6為采用LRFM方法和ASM方法計(jì)算9條強(qiáng)震記錄得到的速度反應(yīng)譜結(jié)果，除個(gè)別臺(tái)站在周期0～1.0 s之間局部極大值處與ASM方法存在極小偏差外，2種方法的速度反應(yīng)譜曲線形態(tài)幾乎一致，均反映出不同阻尼條件下速度反應(yīng)譜的變化趨勢(shì)。由此說(shuō)明，LRFM方法能夠達(dá)到ASM方法的計(jì)算精度，與§3.1得出的結(jié)論一致。

圖6 不同臺(tái)站速度反應(yīng)譜計(jì)算結(jié)果

圖7為采用LRFM方法和ASM方法計(jì)算不同強(qiáng)震記錄的加速度反應(yīng)譜，從圖中可以看出，與ASM方法相比，LRFM方法的計(jì)算誤差主要集中在周期0～1.0 s之間，且隨著阻尼的增大，誤差逐漸增大；5%和10%阻尼計(jì)算結(jié)果的曲線形態(tài)與ASM方法基本重合，在極值點(diǎn)處存在極小誤差，基本可忽略；20%阻尼的計(jì)算結(jié)果在曲線局部與ASM方法的計(jì)算結(jié)果偏差較大，其中最大誤差均出現(xiàn)在曲線極大值或局部極大值處，特別是計(jì)算得到的加速度反應(yīng)譜在卓越周期對(duì)應(yīng)的曲線相鄰位置存在局部極大值時(shí)，相鄰區(qū)域內(nèi)LRFM方法的計(jì)算結(jié)果與ASM方法出現(xiàn)一定分離，如TLN臺(tái)EW方向的加速度反應(yīng)譜計(jì)算結(jié)果。但從整體計(jì)算精度來(lái)看，LRFM方法對(duì)于不同臺(tái)站強(qiáng)震動(dòng)記錄加速度反應(yīng)譜的計(jì)算結(jié)果可以滿足計(jì)算精度的要求。

圖7 不同臺(tái)站加速度反應(yīng)譜計(jì)算結(jié)果

5 結(jié) 語(yǔ)

本文將線性遞推濾波算法(LRFM)應(yīng)用于地震反應(yīng)譜計(jì)算中，對(duì)比分析不同方法計(jì)算合成數(shù)據(jù)和實(shí)際強(qiáng)震動(dòng)記錄的結(jié)果，得出以下結(jié)論：

1)LRFM方法計(jì)算合成數(shù)據(jù)得到的反應(yīng)譜精度和穩(wěn)定性顯著優(yōu)于Duhamel逐步積分法，略優(yōu)于Newmark-β方法。與ASM方法相比，LRFM方法計(jì)算位移、速度和加速度反應(yīng)譜的誤差均小于2%，滿足精度要求。

2)LRFM方法的計(jì)算效率與Duhamel逐步積分法相當(dāng)，且顯著優(yōu)于目前普遍采用的Newmark-β方法，因此對(duì)于快速處理強(qiáng)震及大地震產(chǎn)生的數(shù)量多且持續(xù)時(shí)間長(zhǎng)的強(qiáng)震動(dòng)記錄具有重要意義。

3)對(duì)于實(shí)際強(qiáng)震動(dòng)記錄數(shù)據(jù)的處理效果，LRFM方法與ASM方法計(jì)算的速度反應(yīng)譜結(jié)果一致。

4)與ASM方法相比，LRFM方法計(jì)算得到實(shí)際強(qiáng)震動(dòng)加速度反應(yīng)譜結(jié)果產(chǎn)生的誤差主要集中在周期0～1.0 s之間，且隨著阻尼的增大誤差逐漸增大，但計(jì)算結(jié)果總體滿足精度要求。

5)針對(duì)不同卓越周期和頻譜范圍的實(shí)際強(qiáng)震動(dòng)記錄，采用LRFM方法均可得到可靠的計(jì)算結(jié)果，穩(wěn)定性和適用性較好。

綜上所述，本文提出的LRFM方法計(jì)算耗時(shí)短，結(jié)果精度高，優(yōu)于目前普遍采用的Newmark-β方法，可作為一種快速、高效的反應(yīng)譜計(jì)算方法。