赵博宇，丁 勇，吴 斌

(哈尔滨工业大学 土木工程学院，哈尔滨 150090)

地震模拟振动台实验在地震工程探究领域有着重要的作用。几十年来众多研究机构兴建各种形式的振动台以满足不同要求的抗震实验。据不完全统计，日本已建造29座振动台，美国18座。 国内起步较晚，于70年代起逐渐采用振动台。例如中国建筑科学研究院抗震研究所3×3m单向模拟振动台，国家地震局工程力学研究所5×5m的双向模拟振动台，同济大学4×4m双向模拟振动台，哈尔滨工业大学3×4m单向模拟振动台。国家有关地震模拟振动台的标准也日益完善,基于振动台的结构抗震性能分析也被广泛研究[1]。

振动台的台面支撑系统多种多样，目前比较常见的有气浮式支撑，滚珠轴承式支撑，静压导轨支撑，交叉十字形钢板弹簧铰支撑[2]。实际上，支撑系统同样伴随台面运动，但是支撑系统究竟对于台面运动影响多大很少被研究。不同样式的支撑系统与台面间的运动模式也不同。随着振动台的长期使用，支撑系统的长期时变效应也会引起振动台的模型的变化。定期对振动台进行模型识别能够使实验人员更加了解振动台现状，从而加以维护。同样，获取准确的振动台的模型可促进振动台的改造与发展。某些实验同样需要得到准确的振动台模型，例如当无法布置剪力量测装置时，一个准确的振动台模型对于结构基底剪力的推算显得尤为重要。

系统识别的手段有多种，近年来，卡尔曼估计算法被广泛应用。1960年，卡尔曼首先创立了卡尔曼滤波算法[3]，依据状态转移的分析方法，对系统应用了矢量矩阵的表示方法，建立了卡尔曼滤波方法，其运算过程清晰简单，对于线性体系有着较好的估计效果。为了将卡尔曼估计算法应用于非线性领域，衍生出了等效线性化的近似方法[4-6]，从而发展出扩展卡尔曼滤波公式，解决了非线性问题。扩展的卡尔曼滤波算法将系统参数作为增广状态向量并入系统状态方程与观测方程中，在对系统的动态测量数据进行滤波的同时，给出系统的参数估计。为了改善卡尔曼滤波算法对非线性系统的识别精度，又发展出了无迹卡尔曼滤波算法[7]等方法。基于对卡尔曼滤波技术的深入研究，其应用也越发广泛[8-10]，例如基于EKF的桥梁结构物理参数识别[11]及橡胶隔震支座参数识别[12]等。

随着振动台使用时间增长，系统的刚度和阻尼均有变化，对振动台模型的识别是有必要的。使用时间较长的振动台，其阻尼特性大多难以确定，且常为时变阻尼，利用一般最小二乘法所求得系统参数往往难以与实际测量结果吻合。考虑到地震模拟振动台系统的复杂特性，尤其是使用时间较久的振动台阻尼的不确定性，本文利用扩展卡尔曼估计算法对哈尔滨工业大学的结构抗震实验室的地震模拟振动台进行了识别。通过仿真及3组不同的实验数据最终确定了振动台的模型，对其质量，刚度，阻尼进行了识别。在不同加载方式下，利用识别的结构参数重构的台面响应与实测台面响应的对比结果表明，本文所识别的振动台参数准确可用。

1 哈尔滨工业大学地震模拟振动台

本文的实验对象为哈尔滨工业大学结构抗震实验室的3×4m单向地震模拟振动台，见图1。台面由等截面焊接钢结构组成。台面支撑系统采用十字型钢板弹簧铰支撑，见图2，其优势在于构造简单，可以不设置侧面导向支撑而保证台面足够的侧向刚度和转动刚度。基础采用总重为720t，主体为13×8×3的倒斗形混凝土基础[2]。

图1 哈尔滨工业大学地震模拟振动台俯视图

1.台面；2.上腿立板；3.立板；4.下腿立板；5.钢板弹簧铰

2 扩展卡尔曼滤波算法识别结构参数

2.1 扩展卡尔曼滤波算法

N个自由度的结构运动微分方程为

(1)

式中M，C，K分别为n×n维质量矩阵，阻尼矩阵和刚度矩阵；x(t)=[x1,x2,x3,…xn]T为n阶位移向量；J为激励位置矩阵；F为n阶外荷载向量。取2n+m阶Z向量为状态向量

(2)

(3)

式中w(t)为零均值过程噪声。

系统的离散的观测方程可以表达为

Y(i)=h[Z(i),i]+v(i)

(4)

式中Y(i)为在第i×Δt时刻的l阶观测向量，v(i)为观测噪声向量，其为0均值的高斯白噪声。

使用扩展卡尔曼估计算法，需将系统状态方程与观测方程线性化，可将式 (3)，(4)写为

f[Z(t),t]≈f[Z(i/i),t]+G(i)(Z-Z(i/i))

(5)

h(Z(i+1))≈h(Z(i+1/i))+

H(i+1)(Z(i+1)-Z(i+1/i))

(6)

式中Z(i/i)为状态量在第t=i×Δt时刻的最优估计值，Z(i+1/i)为基于第t=i×Δt时刻观测，对于第t=(i+1)×Δt时刻状态量的估计值。式(5)，(6)中

G(i)=[∂f(Z,t)/∂Z]Z=Z(i/i)

(7)

H(i+1)=

[∂h(Z(i+1)/∂Z(i+1))]Z(i+1)=Zi+1/i

(8)

由扩展卡尔曼滤波估计得最优估计为

Z(i+1/i+1)=Z(i+1/i)+

K[Y(i+1)-H(i+1)Z(i+1/i)]

(9)

式(6)，(9)中Z(i+1/i)由下式计算

(10)

式(9)中K为卡尔曼滤波增益

K=P(i+1/i)HT(i+1)

[H(i+1)P(i+1/i)HT(i+1)+R(i+1)])-1

(11)

式(11)中H和h(Z(i+1))由式(6)和(8)给出，Z(i+1/i)的协方差矩阵

P(i+1/i)=Φ(i+1/i)P(i/i)

ΦT(i+1/i)+Qk

(12)

式(12)中状态转移矩阵

(13)

最优估计协方差矩阵可表示为

P(i+1/i+1)=E[(Z(i+1)-Z(i+1/i+1))

(Z(i+1)-Z(i+1/i+1))T]=

[I-KH(i+1)]P(i+1/i)[I-KH(i+1)]T+

KR(i+1)KT

(14)

综上所述，扩展的卡尔曼滤波估计可由下述步骤进行，

(15a)

P(i+1/i)=

Φ(i+1/i)P(i/i)ΦT(i+1/i)+Qk

(15b)

Z(i+1/i+1)=Z(i+1/i)+

K[Y(i+1)-H(i+1)Z(i+1/i)]

(15c)

K=P(i+1/i)HT(i+1)

[H(i+1)P(i+1/i)

HT(i+1)+R(i+1)]-1

(15d)

P(i+1/i+1)=E[(Z(i+1)-Z(i+1/i+1))

(Z(i+1)-Z(i+1/i+1))T]=

[I-KH(i+1)]P(i+1/i)[I-KH(i+1)]T+

KR(i+1)KT

(15e)

H(i+1)=

[∂h(Z(i+1)/∂Z(i+1))]Z(i+1)=Z(i+1/i)

(15f)

(15g)

2.2 扩展卡尔曼滤波算法识别地震模拟振动台模型

一般单向地震模拟振动台可以简化为平面内单自由度结构，有运动方程

(16)

式中m,c，k分别为地震模拟振动台的质量，阻尼和刚度；x为台面位移，F为台面受力。认为地震模拟振动台的质量，刚度，阻尼皆为未知，通过测量台面位移，对地震模拟振动台模型进行修正。

取状态向量

Z(t)=[Xi,X2,X3,X4,X5]T=[Xd,Xv,θ]T=

(17)

(18)

若取台面位移为观测量，则系统观测方程为

Y(i+1)=[1 0 0 0 0]×Z(i+1/i)

(19)

由(15a)-(15g)即可计算得到最优状态量，从而更新修正地震模拟振动台模型

3 仿真算例

仿真模型：单自由度结构，结构待识别参数有质量M=4t，阻尼C=0.5 kN/(m/s)，刚度K=443 kN/m。输入的地震动选取峰值加速度为350 gal的El Centro(1940NS)波，噪声水平10%。采用扩展卡尔曼滤波估计算法。通过观测位移响应，识别振动台的质量，刚度，阻尼。计算时长20 s，计算步长0.01 s。

仿真结果如图3所示。

图3(e)中计算加速度是由估计的结构参数，基于结构运动微分方程方程(16)重构计算而来。测量加速度由真实结构参数，基于运动微分方程计算得来。由式(19)计算得到计算加速度时程与测量加速度时程的相对误差为6.77%，由此认为参数识别比较准确。

(19)

4 地震模拟振动台的模型识别与验证

4.1 模型识别

针对哈尔滨工业大学结构抗震实验室的地震模拟振动台进行了三组实验，通过首节对哈尔滨工业大学地震模拟振动台的介绍，将其计算模型简化为单自由度结构。采用扩展卡尔曼估计算法识别地震模拟振动台质量，刚度与阻尼。

三组实验加载都依据1940年5月19日Imperial Valley地震El Centro(1940NS)台站测得的地面运动位移记录，振动台施加的命令位移峰值相应按比例变为为10.3mm，14.7mm，18.2mm。第三组实验识别效果如图4所示。

图5中加速度计算值是由各组估计的结构参数，基于结构运动微分方程方程(15)重构计算而来。加速度测量值是由布置于台面的加速度传感器测得。三组实验的识别结果见表2。各组加速度测量值时程与加速度计算值时程的相对误差见表3。

图4 第三组实验参数识别结果

图5 各组实验加速度计算值与测量值对比

表2 三组实验识别结果

表3 计算与测量加速度时程相对误差Tab.3 Relative errors between calculatedand measured accelerations

由表2识别结果可见，三组实验识别结果相差不大，因此可估计地震模拟振动台质量为5.7t，刚度约为46 kN/m，阻尼约为1.5 kN/(m/s)。由表3可见各组重构计算加速度响应时程皆与实测台面加速度响应时程存在一定误差，误差主要由时程相位差导致，然而由图5可见,各组加速度计算值与测量值峰值和形状吻合良好，由此验证了参数识别的准确性。

4.2 模型验证

利用4.1节所识别的振动台模型，通过另外两种不同加载方式进一步验证所识别模型的准确性。振动台参数选取为：质量5.7t，刚度46 kN/m，阻尼1.5 kN/(m/s)。加载由位移控制。

工况1：首先进行时长3 s，周期为1秒的正弦加载，位移峰值为10 mm。正弦加载后空载1 s进行位移峰值为23.4mm的El Centro(1940NS)波6秒加载，其对应加速度峰值为220 gal。

工况2：进行持时为25 s，位移峰值为14.2 mm的Taft波加载，其对应加速度峰值为150 gal。

图6(a)为工况一下力传感器所测台面荷载时程。图6(b)、(c)分别为利用识别的振动台模型参数和实测台面受力，基于运动微分方程计算得到的台面加速度响应、位移响应与实测响应的对比图。图7(a)为工况二下力传感器所测台面荷载时程。图7(b)～(c)为响应计算值与实测值时程对比，计算方法与工况一相同。

从图6和图8可以看出，4.1节所识别的振动台模型参数可以准确的计算出不同地震波作用下振动台的响应。工况1和工况2下的位移计算时程与实测时程的相对误差分别为11.7%、13.3%，加速度计算时程与实测时程的相对误差分别为17.1%，18.4%。相对误差主要来源于时程曲线的相位差，相位差的出现与本次试验中假定阻尼为时不变阻尼有关。然而，加速度与位移时程的峰值大小和形状皆吻合良好，由此可判断本文所提出方法有效的建立了地震模拟振动台模型。

图6 工况一验证

图7 工况二验证

5 结 论

本文建立了哈尔滨工业大学地震模拟振动台的模型。首先进行了单自由度结构仿真，利用扩展的卡尔曼估计算法，将结构质量，阻尼，刚度同时引入状态量，仿真验证了算法可以基于位移响应准确的识别结构参数。随后，通过3组振动台实验成功的识别了哈尔滨工业大学振动台的模型。最后利用其它不同加载方式下的台面计算响应与实测响应的时程对比可知，所识别的振动台参数准确可信。由此本文解决了振动台模型难以准确建立的问题。此研究不但成功的完成了扩展卡尔曼滤波算法在噪声较大的环境下的多参数识别，而且更新了长期不确定的振动台模型。本文研究方法同样适用于其他类型地震模拟振动台模型的建立与修正，并可扩展于某些难以测定参数的复杂结构模型识别。

[1] 公茂盛,谢礼立,欧进萍.结构振动台模型模态参数识别新方法研究[J].振动工程学报,2010,23(2):230-236.

GONG Mao-sheng,XIE Li-li,OU Jin-ping.A method for modal parameter identification of structural shaking table model[J].Journal of Vibration Engineering，2010，23(2): 230-236.

[2] 潘景龙.单向模拟地震振动台设计中的若干问题讨论[J].哈尔滨建筑工程学院学报，1990，23(2):90-99.

PAN Jing-long.Discuss on a number of problem in design of the unidirectional shaking table[J].Journal of Harbin Institute of Technology，1990，23(2):90-99.

[3] Kalman R E.A new approach to linear filtering and prediction problems[J].Transaction of the ASME-Journal of Basic Engineering,1960: 35-45.

[4] Yang J N，Lin S.On-line identification of nonlinear hysteretic structures using an adaptive tracking technique[J].International Journal of Non-linear Mechanics，2004，39:1481-1491.

[5] Hoshiya M，Saito E.Structural identification by extended Kalman filter[J].Journal of Engineering Mechanics，1984，110(12): 1757-1770.

[6] Corigliano A，Mariani S.Parameter identification in explicit structural dynamics: performance of the extended Kalman filter[J].Computer Methods in Applied Mechanics and Engineering，2004，193(36): 3807-3835.

[7] Wan E A，Van Der Merwe R.The unscented Kalman filter for nonlinear estimation[C]//Adaptive Systems for Signal Processing，Communications，and Control Symposium 2000.AS-SPCC.The IEEE 2000.IEEE，2000: 153-158.

[8] Hoshiya M，Saito E,.Structural identiifcation by extended kalman filter[J].Journal of Engineering Mechanics (ASCE),1984,110(12): 1757-1771.

[9] Maruyama O，Hoshiya M.System identification of an experimental model by extended Kalman filter[J].Proceedings of Structural Safety and Reliability，ICOSSA Swet & Zeitinger: Lisse，CD-ROM,2001.

[10] Zhou L，Wu S，Yang J N.Experimental study of an adaptive extended Kalman filter for structural damage identification[J].Journal of Infrastructure Systems，2008，14(1): 42-51.

[11] 吴子燕,丁兰,刘书奎.基于广义卡尔曼滤波的桥梁结构物理参数识别的子结构法[J].计算力学学报,2007,24(4):425-428.

WU Zi-yan,DING Lan,LIU Shu-kui.Identification on physical parameters of bridge structures based on extended Kalman filter[J].Chinese Journal of Computational Mechanics，2007，24(4):425-428.

[12] Yin Q，Zhou L，Yang J N.An experimental study on AEKF method for health monitoring of base-isolated structures[C]//Proc.of SPIE Vol，2010，7647: 764709-1.