徐玉聰

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

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

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

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

1 單自由度震動體系原理

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

(1)

定義時間加權(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 計算原理

在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ù)分析計算

3.1 計算精度分析

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

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種方法計算得到的位移反應(yīng)譜如圖1(a)所示，從圖中可以看出，4種方法計算得到的曲線形態(tài)一致，均反映出輸入周期在0.5 s和2 s處的模擬地震動信息。位移反應(yīng)譜的相對誤差(圖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)的相對誤差結(jié)果優(yōu)于Duhamel逐步積分法和Newmark-β方法，且整體誤差波動較小。

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

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

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

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

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

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

3.2 計算效率分析

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

圖4 反應(yīng)譜計算耗時對比

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

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

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

圖5 地震動加速度時程

表1 強震動記錄參數(shù)

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

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

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

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

5 結(jié) 語

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

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

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

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

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

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

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