基于擴展卡爾曼估計算法的地震模擬振動臺模型識別

趙博宇，丁 勇，吳 斌

(哈爾濱工業(yè)大學 土木工程學院，哈爾濱 150090)

地震模擬振動臺實驗在地震工程探究領(lǐng)域有著重要的作用。幾十年來眾多研究機構(gòu)興建各種形式的振動臺以滿足不同要求的抗震實驗。據(jù)不完全統(tǒng)計，日本已建造29座振動臺，美國18座。 國內(nèi)起步較晚，于70年代起逐漸采用振動臺。例如中國建筑科學研究院抗震研究所3×3m單向模擬振動臺，國家地震局工程力學研究所5×5m的雙向模擬振動臺，同濟大學4×4m雙向模擬振動臺，哈爾濱工業(yè)大學3×4m單向模擬振動臺。國家有關(guān)地震模擬振動臺的標準也日益完善,基于振動臺的結(jié)構(gòu)抗震性能分析也被廣泛研究[1]。

振動臺的臺面支撐系統(tǒng)多種多樣，目前比較常見的有氣浮式支撐，滾珠軸承式支撐，靜壓導軌支撐，交叉十字形鋼板彈簧鉸支撐[2]。實際上，支撐系統(tǒng)同樣伴隨臺面運動，但是支撐系統(tǒng)究竟對于臺面運動影響多大很少被研究。不同樣式的支撐系統(tǒng)與臺面間的運動模式也不同。隨著振動臺的長期使用，支撐系統(tǒng)的長期時變效應也會引起振動臺的模型的變化。定期對振動臺進行模型識別能夠使實驗人員更加了解振動臺現(xiàn)狀，從而加以維護。同樣，獲取準確的振動臺的模型可促進振動臺的改造與發(fā)展。某些實驗同樣需要得到準確的振動臺模型，例如當無法布置剪力量測裝置時，一個準確的振動臺模型對于結(jié)構(gòu)基底剪力的推算顯得尤為重要。

系統(tǒng)識別的手段有多種，近年來，卡爾曼估計算法被廣泛應用。1960年，卡爾曼首先創(chuàng)立了卡爾曼濾波算法[3]，依據(jù)狀態(tài)轉(zhuǎn)移的分析方法，對系統(tǒng)應用了矢量矩陣的表示方法，建立了卡爾曼濾波方法，其運算過程清晰簡單，對于線性體系有著較好的估計效果。為了將卡爾曼估計算法應用于非線性領(lǐng)域，衍生出了等效線性化的近似方法[4-6]，從而發(fā)展出擴展卡爾曼濾波公式，解決了非線性問題。擴展的卡爾曼濾波算法將系統(tǒng)參數(shù)作為增廣狀態(tài)向量并入系統(tǒng)狀態(tài)方程與觀測方程中，在對系統(tǒng)的動態(tài)測量數(shù)據(jù)進行濾波的同時，給出系統(tǒng)的參數(shù)估計。為了改善卡爾曼濾波算法對非線性系統(tǒng)的識別精度，又發(fā)展出了無跡卡爾曼濾波算法[7]等方法?；趯柭鼮V波技術(shù)的深入研究，其應用也越發(fā)廣泛[8-10]，例如基于EKF的橋梁結(jié)構(gòu)物理參數(shù)識別[11]及橡膠隔震支座參數(shù)識別[12]等。

隨著振動臺使用時間增長，系統(tǒng)的剛度和阻尼均有變化，對振動臺模型的識別是有必要的。使用時間較長的振動臺，其阻尼特性大多難以確定，且常為時變阻尼，利用一般最小二乘法所求得系統(tǒng)參數(shù)往往難以與實際測量結(jié)果吻合?？紤]到地震模擬振動臺系統(tǒng)的復雜特性，尤其是使用時間較久的振動臺阻尼的不確定性，本文利用擴展卡爾曼估計算法對哈爾濱工業(yè)大學的結(jié)構(gòu)抗震實驗室的地震模擬振動臺進行了識別。通過仿真及3組不同的實驗數(shù)據(jù)最終確定了振動臺的模型，對其質(zhì)量，剛度，阻尼進行了識別。在不同加載方式下，利用識別的結(jié)構(gòu)參數(shù)重構(gòu)的臺面響應與實測臺面響應的對比結(jié)果表明，本文所識別的振動臺參數(shù)準確可用。

1 哈爾濱工業(yè)大學地震模擬振動臺

本文的實驗對象為哈爾濱工業(yè)大學結(jié)構(gòu)抗震實驗室的3×4m單向地震模擬振動臺，見圖1。臺面由等截面焊接鋼結(jié)構(gòu)組成。臺面支撐系統(tǒng)采用十字型鋼板彈簧鉸支撐，見圖2，其優(yōu)勢在于構(gòu)造簡單，可以不設(shè)置側(cè)面導向支撐而保證臺面足夠的側(cè)向剛度和轉(zhuǎn)動剛度?；A(chǔ)采用總重為720t，主體為13×8×3的倒斗形混凝土基礎(chǔ)[2]。

圖1 哈爾濱工業(yè)大學地震模擬振動臺俯視圖

1.臺面；2.上腿立板；3.立板；4.下腿立板；5.鋼板彈簧鉸

2 擴展卡爾曼濾波算法識別結(jié)構(gòu)參數(shù)

2.1 擴展卡爾曼濾波算法

N個自由度的結(jié)構(gòu)運動微分方程為

(1)

式中M，C，K分別為n×n維質(zhì)量矩陣，阻尼矩陣和剛度矩陣；x(t)=[x1,x2,x3,…xn]T為n階位移向量；J為激勵位置矩陣；F為n階外荷載向量。取2n+m階Z向量為狀態(tài)向量

(2)

(3)

式中w(t)為零均值過程噪聲。

系統(tǒng)的離散的觀測方程可以表達為

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

(4)

式中Y(i)為在第i×Δt時刻的l階觀測向量，v(i)為觀測噪聲向量，其為0均值的高斯白噪聲。

使用擴展卡爾曼估計算法，需將系統(tǒng)狀態(tài)方程與觀測方程線性化，可將式 (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=i×Δt時刻的最優(yōu)估計值，Z(i+1/i)為基于第t=i×Δt時刻觀測，對于第t=(i+1)×Δt時刻狀態(tài)量的估計值。式(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)

由擴展卡爾曼濾波估計得最優(yōu)估計為

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)的協(xié)方差矩陣

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

ΦT(i+1/i)+Qk

(12)

式(12)中狀態(tài)轉(zhuǎn)移矩陣

(13)

最優(yōu)估計協(xié)方差矩陣可表示為

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 擴展卡爾曼濾波算法識別地震模擬振動臺模型

一般單向地震模擬振動臺可以簡化為平面內(nèi)單自由度結(jié)構(gòu)，有運動方程

(16)

式中m,c，k分別為地震模擬振動臺的質(zhì)量，阻尼和剛度；x為臺面位移，F(xiàn)為臺面受力。認為地震模擬振動臺的質(zhì)量，剛度，阻尼皆為未知，通過測量臺面位移，對地震模擬振動臺模型進行修正。

取狀態(tài)向量

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

(17)

(18)

若取臺面位移為觀測量，則系統(tǒng)觀測方程為

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

(19)

由(15a)-(15g)即可計算得到最優(yōu)狀態(tài)量，從而更新修正地震模擬振動臺模型

3 仿真算例

仿真模型：單自由度結(jié)構(gòu)，結(jié)構(gòu)待識別參數(shù)有質(zhì)量M=4t，阻尼C=0.5 kN/(m/s)，剛度K=443 kN/m。輸入的地震動選取峰值加速度為350 gal的El Centro(1940NS)波，噪聲水平10%。采用擴展卡爾曼濾波估計算法。通過觀測位移響應，識別振動臺的質(zhì)量，剛度，阻尼。計算時長20 s，計算步長0.01 s。

仿真結(jié)果如圖3所示。

圖3(e)中計算加速度是由估計的結(jié)構(gòu)參數(shù)，基于結(jié)構(gòu)運動微分方程方程(16)重構(gòu)計算而來。測量加速度由真實結(jié)構(gòu)參數(shù)，基于運動微分方程計算得來。由式(19)計算得到計算加速度時程與測量加速度時程的相對誤差為6.77%，由此認為參數(shù)識別比較準確。

(19)

4 地震模擬振動臺的模型識別與驗證

4.1 模型識別

針對哈爾濱工業(yè)大學結(jié)構(gòu)抗震實驗室的地震模擬振動臺進行了三組實驗，通過首節(jié)對哈爾濱工業(yè)大學地震模擬振動臺的介紹，將其計算模型簡化為單自由度結(jié)構(gòu)。采用擴展卡爾曼估計算法識別地震模擬振動臺質(zhì)量，剛度與阻尼。

三組實驗加載都依據(jù)1940年5月19日Imperial Valley地震El Centro(1940NS)臺站測得的地面運動位移記錄，振動臺施加的命令位移峰值相應按比例變?yōu)闉?0.3mm，14.7mm，18.2mm。第三組實驗識別效果如圖4所示。

圖5中加速度計算值是由各組估計的結(jié)構(gòu)參數(shù)，基于結(jié)構(gòu)運動微分方程方程(15)重構(gòu)計算而來。加速度測量值是由布置于臺面的加速度傳感器測得。三組實驗的識別結(jié)果見表2。各組加速度測量值時程與加速度計算值時程的相對誤差見表3。

圖4 第三組實驗參數(shù)識別結(jié)果

圖5 各組實驗加速度計算值與測量值對比

表2 三組實驗識別結(jié)果

表3 計算與測量加速度時程相對誤差Tab.3 Relative errors between calculatedand measured accelerations

由表2識別結(jié)果可見，三組實驗識別結(jié)果相差不大，因此可估計地震模擬振動臺質(zhì)量為5.7t，剛度約為46 kN/m，阻尼約為1.5 kN/(m/s)。由表3可見各組重構(gòu)計算加速度響應時程皆與實測臺面加速度響應時程存在一定誤差，誤差主要由時程相位差導致，然而由圖5可見,各組加速度計算值與測量值峰值和形狀吻合良好，由此驗證了參數(shù)識別的準確性。

4.2 模型驗證

利用4.1節(jié)所識別的振動臺模型，通過另外兩種不同加載方式進一步驗證所識別模型的準確性。振動臺參數(shù)選取為：質(zhì)量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)分別為利用識別的振動臺模型參數(shù)和實測臺面受力，基于運動微分方程計算得到的臺面加速度響應、位移響應與實測響應的對比圖。圖7(a)為工況二下力傳感器所測臺面荷載時程。圖7(b)～(c)為響應計算值與實測值時程對比，計算方法與工況一相同。

從圖6和圖8可以看出，4.1節(jié)所識別的振動臺模型參數(shù)可以準確的計算出不同地震波作用下振動臺的響應。工況1和工況2下的位移計算時程與實測時程的相對誤差分別為11.7%、13.3%，加速度計算時程與實測時程的相對誤差分別為17.1%，18.4%。相對誤差主要來源于時程曲線的相位差，相位差的出現(xiàn)與本次試驗中假定阻尼為時不變阻尼有關(guān)。然而，加速度與位移時程的峰值大小和形狀皆吻合良好，由此可判斷本文所提出方法有效的建立了地震模擬振動臺模型。

圖6 工況一驗證

圖7 工況二驗證

5 結(jié) 論

本文建立了哈爾濱工業(yè)大學地震模擬振動臺的模型。首先進行了單自由度結(jié)構(gòu)仿真，利用擴展的卡爾曼估計算法，將結(jié)構(gòu)質(zhì)量，阻尼，剛度同時引入狀態(tài)量，仿真驗證了算法可以基于位移響應準確的識別結(jié)構(gòu)參數(shù)。隨后，通過3組振動臺實驗成功的識別了哈爾濱工業(yè)大學振動臺的模型。最后利用其它不同加載方式下的臺面計算響應與實測響應的時程對比可知，所識別的振動臺參數(shù)準確可信。由此本文解決了振動臺模型難以準確建立的問題。此研究不但成功的完成了擴展卡爾曼濾波算法在噪聲較大的環(huán)境下的多參數(shù)識別，而且更新了長期不確定的振動臺模型。本文研究方法同樣適用于其他類型地震模擬振動臺模型的建立與修正，并可擴展于某些難以測定參數(shù)的復雜結(jié)構(gòu)模型識別。

[1] 公茂盛,謝禮立,歐進萍.結(jié)構(gòu)振動臺模型模態(tài)參數(shù)識別新方法研究[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] 潘景龍.單向模擬地震振動臺設(shè)計中的若干問題討論[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] 吳子燕,丁蘭,劉書奎.基于廣義卡爾曼濾波的橋梁結(jié)構(gòu)物理參數(shù)識別的子結(jié)構(gòu)法[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.