徐玉聰

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

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

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

本文基于單自由度體系運動方程推導計算地震反應譜的一般表達形式，并從傳遞函數(shù)的形式和濾波系數(shù)入手，推導用于計算反應譜的遞推表達式，然后結(jié)合已有研究提出的滿足高頻和低頻收斂條件對應的獨立積分常數(shù)即可計算得到反應譜結(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。設置反應譜阻尼比為5%，周期為0～4 s，時間間隔為0.02 s，分別利用LRFM方法、Duhamel逐步積分法、Newmark-β方法和精確解法(ASM)[12]計算得到反應譜結(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種方法計算得到的位移反應譜如圖1(a)所示，從圖中可以看出，4種方法計算得到的曲線形態(tài)一致，均反映出輸入周期在0.5 s和2 s處的模擬地震動信息。位移反應譜的相對誤差(圖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 位移反應譜計算結(jié)果及誤差

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

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

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

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

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

3.2 計算效率分析

在計算地震反應譜時，為獲得更充分的反應譜信息，一般選取輸入時長較長的強震動記錄，特別是在計算強震動條件下的加速度反應譜時，強震動記錄時長可達300 s以上，因此會產(chǎn)生大量強震動記錄，如果反應譜計算耗時太長，將不利于反應譜的進一步分析。為對比分析LRFM方法、Duhamel逐步積分法和Newmark-β方法的計算效率，利用式(18)作為系統(tǒng)輸入，合成數(shù)據(jù)計算時長分別設置為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 反應譜計算耗時對比

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

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

為驗證LRFM方法計算實際強震動記錄反應譜的效果，選取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，這些記錄可基本反映該地震的強震動特征，計算得到的反應譜能夠代表在該地震作用下不同位置的地面運動，可以很好地驗證LRFM方法對實際強震動記錄反應譜的計算效果。計算反應譜的阻尼系數(shù)選取5%、10%、20%，自振周期為0～4 s，時間間隔為0.02 s，選取的強震記錄時程曲線如圖5所示。本文選取的強震記錄具有一般性，主頻可基本覆蓋從低頻到高頻不同范圍的信息。

圖5 地震動加速度時程

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

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

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

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

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

5 結(jié) 語

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

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

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

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

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

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

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