基于线性递推滤波算法的地震反应谱计算方法

徐玉聪

1 长江三峡勘测研究院有限公司(武汉), 武汉市创业街99号，430074

地震反应谱是地震工程研究中的一项重要内容，通过研究单自由度震动体系下质点的运动情况来描述地震条件下结构的最大动力反应[1-2]和地震动的频谱特性。反应谱在强震监测计算分析中既可反映局部场地的地震强度特征，又可用来估算地震作用下工程结构的受力情况，大量的反应谱计算结果还可应用于结构抗震设计。因此，采用有效的方法提高反应谱的精度和计算效率具有重要意义。

反应谱的计算方法主要集中在时域和频域两个维度。频域计算方法的基本原理是在频域计算地震动和传递函数傅里叶频谱的乘积[3]，然后通过傅里叶逆变换得到时间域的反应谱，需要已知频域内精确的传递函数[4]以及完整的地震记录来完成计算，其中长周期反应谱的计算精度受地震记录长度影响较大，因此频域不适合进行实时快速反应谱计算。反应谱时域计算方法一直为研究重点，其基本原理是以地震反应量所表述的动力平衡方程来寻找数字滤波器传递函数或递归公式[5]。时域计算方法有Duhamel逐步积分法、Newmark-β方法[6]、精确解法(ASM)[7]和Z变换方法[8]等，这些方法的计算精度主要受假定时间步长内加速度的特定分布形式所制约[9]。

本文基于单自由度体系运动方程推导计算地震反应谱的一般表达形式，并从传递函数的形式和滤波系数入手，推导用于计算反应谱的递推表达式，然后结合已有研究提出的满足高频和低频收敛条件对应的独立积分常数即可计算得到反应谱结果。

1 单自由度震动体系原理

当地面以一定加速度运动时，单自由度震动体系的运动方程为[2]：

(1)

定义时间加权函数W(t)为：

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

(2)

(3)

(4)

(5)

(6)

式中，Λ1、Λ2、Λ3、Λ4、Λ5、Λ6为一系列自由参数。

将式(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变换域中输入、输出和传递函数的形式为：

Y(Z)=

(12)

式中，na为反馈滤波器阶数；nb为前馈滤波器阶数。可将式(12)表示为差分方程：

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

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

…+Bnay(m-na)

(13)

根据Zhou等[10]提出的思路，构造3阶反馈滤波器系数和4阶前馈滤波器系数，令[1,-A1,A2,-A3]、[B1,B2,B3,B4]为2组滤波器系数，其中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 合成数据分析计算

3.1 计算精度分析

为验证本文方法的计算精度，采用由2个频率合成的正弦简谐波f(t)作为单自由度震动体系的系统输入，如式(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]计算得到反应谱结果，并以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种方法计算得到的曲线形态一致，均反映出输入周期在0.5 s和2 s处的模拟地震动信息。位移反应谱的相对误差(图1(b))结果显示，LRFM方法、Duhamel逐步积分法和Newmark-β方法在周期0～1 s内均存在一定误差，在周期2～4 s内误差较小，在0.5 s处具有明显突变；LRFM方法在周期0～1 s内的相对误差结果优于Duhamel逐步积分法和Newmark-β方法，且整体误差波动较小。

图1 位移反应谱计算结果及误差

利用合成数据计算得到的速度反应谱和误差结果如图2所示，从图中可以看出，4种方法计算得到的速度反应谱曲线几乎完全重合，其中Duhamel逐步积分法的相对误差在短周期0～1 s内跳动较大，最大达到18%；而Newmark-β方法和LRFM方法的整体误差较小，基本可达到ASM方法的计算效果。

图2 速度反应谱计算结果及误差

图3为4种方法计算得到的加速度反应谱和误差结果，从图中可以看出，4种方法计算的加速度反应谱曲线仅在个别峰值点存在突变，而LRFM方法、Duhamel逐步积分法及Newmark-β方法的相对误差均小于2%，3种方法的相对误差均集中在短周期部分。LRFM方法的相对误差基本呈线性正相关变化，而Duhamel逐步积分法和Newmark-β方法的相对误差在正负方向来回震荡，且在0.5 s处存在突变。

图3 加速度反应谱计算结果及误差

通过对合成数据的计算结果进行分析可知，LRFM方法的计算精度和稳定性优于Duhamel逐步积分法和Newmark-β方法，在位移最大值处的误差较小，而加速度反应谱的计算误差整体略大于另外2种方法，但相对误差小于2%，计算结果满足精度要求。

3.2 计算效率分析

在计算地震反应谱时，为获得更充分的反应谱信息，一般选取输入时长较长的强震动记录，特别是在计算强震动条件下的加速度反应谱时，强震动记录时长可达300 s以上，因此会产生大量强震动记录，如果反应谱计算耗时太长，将不利于反应谱的进一步分析。为对比分析LRFM方法、Duhamel逐步积分法和Newmark-β方法的计算效率，利用式(18)作为系统输入，合成数据计算时长分别设置为30 s、60 s、90 s、180 s和300 s，采样频率为200 Hz，各方法的计算耗时如图4所示。由图可知，Newmark-β方法的计算耗时明显大于另外2种方法，且随着时长的增加，计算耗时呈现线性变化，在300 s处计算耗时接近10 s；LRFM方法和Duhamel逐步积分法的整体耗时在1 s以内，两者的计算效率基本相当，且显著优于Newmark-β方法。

图4 反应谱计算耗时对比

随着我国强震动监测台站的增加，强震及大地震发生后可获取数量多且持续时间长的强震动记录，如2008年汶川地震发生后，仅四川省数字强震台网获得的主震记录就接近400条，大部分纪录的持续时间为150～160 s[13]，最长超过300 s，因此选取合适的方法快速计算得到反应谱以评估场地的地震强度特征或工程结构的受力情况尤为重要。本文提出的LRFM方法的计算效率显著优于目前普遍采用的Newmark-β方法，可快速高效地得到计算结果。

4 实际数据分析计算

为验证LRFM方法计算实际强震动记录反应谱的效果，选取ESM强震数据库收录的SNO、PCB、FOC、TLN、AMT、MSCT、CLF、FOS共8个台站记录的2016-08-24意大利ML6.0地震强震波形数据，具体参数见表1，采样频率均为200 Hz。选取的强震动记录震中距为8.5～57.50 km，PGA为47.13～850.80 Gal，卓越周期范围为0.06～0.56 s，这些记录可基本反映该地震的强震动特征，计算得到的反应谱能够代表在该地震作用下不同位置的地面运动，可以很好地验证LRFM方法对实际强震动记录反应谱的计算效果。计算反应谱的阻尼系数选取5%、10%、20%，自振周期为0～4 s，时间间隔为0.02 s，选取的强震记录时程曲线如图5所示。本文选取的强震记录具有一般性，主频可基本覆盖从低频到高频不同范围的信息。

图5 地震动加速度时程

表1 强震动记录参数

以ASM方法的计算结果作为参照，验证LRFM方法计算实际强震动数据反应谱的效果，另外由于不同方法计算得到的位移反应谱误差相差较小，因此本文仅分析速度反应谱和加速度反应谱的计算结果。图6为采用LRFM方法和ASM方法计算9条强震记录得到的速度反应谱结果，除个别台站在周期0～1.0 s之间局部极大值处与ASM方法存在极小偏差外，2种方法的速度反应谱曲线形态几乎一致，均反映出不同阻尼条件下速度反应谱的变化趋势。由此说明，LRFM方法能够达到ASM方法的计算精度，与§3.1得出的结论一致。

图6 不同台站速度反应谱计算结果

图7为采用LRFM方法和ASM方法计算不同强震记录的加速度反应谱，从图中可以看出，与ASM方法相比，LRFM方法的计算误差主要集中在周期0～1.0 s之间，且随着阻尼的增大，误差逐渐增大；5%和10%阻尼计算结果的曲线形态与ASM方法基本重合，在极值点处存在极小误差，基本可忽略；20%阻尼的计算结果在曲线局部与ASM方法的计算结果偏差较大，其中最大误差均出现在曲线极大值或局部极大值处，特别是计算得到的加速度反应谱在卓越周期对应的曲线相邻位置存在局部极大值时，相邻区域内LRFM方法的计算结果与ASM方法出现一定分离，如TLN台EW方向的加速度反应谱计算结果。但从整体计算精度来看，LRFM方法对于不同台站强震动记录加速度反应谱的计算结果可以满足计算精度的要求。

图7 不同台站加速度反应谱计算结果

5 结 语

本文将线性递推滤波算法(LRFM)应用于地震反应谱计算中，对比分析不同方法计算合成数据和实际强震动记录的结果，得出以下结论：

1)LRFM方法计算合成数据得到的反应谱精度和稳定性显著优于Duhamel逐步积分法，略优于Newmark-β方法。与ASM方法相比，LRFM方法计算位移、速度和加速度反应谱的误差均小于2%，满足精度要求。

2)LRFM方法的计算效率与Duhamel逐步积分法相当，且显著优于目前普遍采用的Newmark-β方法，因此对于快速处理强震及大地震产生的数量多且持续时间长的强震动记录具有重要意义。

3)对于实际强震动记录数据的处理效果，LRFM方法与ASM方法计算的速度反应谱结果一致。

4)与ASM方法相比，LRFM方法计算得到实际强震动加速度反应谱结果产生的误差主要集中在周期0～1.0 s之间，且随着阻尼的增大误差逐渐增大，但计算结果总体满足精度要求。

5)针对不同卓越周期和频谱范围的实际强震动记录，采用LRFM方法均可得到可靠的计算结果，稳定性和适用性较好。

综上所述，本文提出的LRFM方法计算耗时短，结果精度高，优于目前普遍采用的Newmark-β方法，可作为一种快速、高效的反应谱计算方法。