采用隱性卡爾曼濾波器的自適應(yīng)子結(jié)構(gòu)試驗方法

王 濤，張 健，吳 斌

（1．哈爾濱工業(yè)大學(xué)土木工程學(xué)院，黑龍江 哈爾濱 150090；2．黑龍江科技學(xué)院建筑工程學(xué)院，黑龍江 哈爾濱 150027）

引 言

傳統(tǒng)的結(jié)構(gòu)抗震試驗可分為擬靜力試驗、振動臺試驗和擬動力試驗。擬靜力試驗技術(shù)簡單、穩(wěn)定，但是其不考慮地震作用對結(jié)構(gòu)的影響。振動臺試驗?zāi)軌蛑噩F(xiàn)地震動，可以真實地反映地震對結(jié)構(gòu)的作用，但是由于振動臺承載能力及尺寸的限制，往往無法進行大尺度結(jié)構(gòu)試驗。由于小尺度結(jié)構(gòu)模型的動力相似律很難得到滿足，尤其是在彈塑性范圍內(nèi)，試驗結(jié)果往往很難推廣到原型結(jié)構(gòu)中去［1］。

擬動力試驗通過理論計算來考慮結(jié)構(gòu)慣性力和阻尼力的影響，對結(jié)構(gòu)進行慢速加載試驗測量結(jié)構(gòu)的恢復(fù)力。擬動力試驗解決了理論分析計算中恢復(fù)力模型及參數(shù)難以確定的困難，可以獲得結(jié)構(gòu)體系的真實反應(yīng)特征，同時又比振動臺試驗更經(jīng)濟，對試驗裝置要求較低，更易實現(xiàn)。子結(jié)構(gòu)擬動力試驗技術(shù)對結(jié)構(gòu)的關(guān)鍵部位進行足尺實驗或大尺度實驗，其他部位用計算機進行模擬，降低了試驗費用。

傳統(tǒng)的子結(jié)構(gòu)試驗中數(shù)值子結(jié)構(gòu)采用線彈性假定或事先假定的恢復(fù)力模型。假定恢復(fù)力模型與真實恢復(fù)力特性可能存在較大差異，這種差異可能引入較大的額外試驗誤差。當數(shù)值子結(jié)構(gòu)中包含與試驗子結(jié)構(gòu)相同的結(jié)構(gòu)或構(gòu)件時，可以利用試驗子結(jié)構(gòu)試驗觀測數(shù)據(jù)在線預(yù)測或修正相應(yīng)的數(shù)值子結(jié)構(gòu)模型，將這種改進的試驗方法稱為自適應(yīng)子結(jié)構(gòu)試驗方法。Yang等人提出根據(jù)試驗子結(jié)構(gòu)的數(shù)據(jù)在線訓(xùn)練恢復(fù)力神經(jīng)網(wǎng)絡(luò)，并用預(yù)測數(shù)值子結(jié)構(gòu)的恢復(fù)力［2］。

研究發(fā)現(xiàn)，采用神經(jīng)網(wǎng)絡(luò)預(yù)測結(jié)構(gòu)恢復(fù)力時，可能存在恢復(fù)力誤差均值的偏差。而非線性結(jié)構(gòu)對于恢復(fù)力誤差均值的偏差比誤差標準差更敏感，微小的誤差均值偏差可能導(dǎo)致非線性結(jié)構(gòu)位移反應(yīng)很大的差異［3］。

除了可以通過神經(jīng)網(wǎng)絡(luò)方法預(yù)測數(shù)值子結(jié)構(gòu)恢復(fù)力外，還可以通過參數(shù)識別方法（比如隱性卡爾曼濾波器），在線識別恢復(fù)力模型參數(shù)。本文提出基于隱性卡爾曼濾波的自適應(yīng)子結(jié)構(gòu)試驗方法，該方法核心思想是在子結(jié)構(gòu)試驗過程中，通過觀測試驗子結(jié)構(gòu)恢復(fù)力，采用隱性卡爾曼濾波器在線識別試驗子結(jié)構(gòu)模型參數(shù)，實時更新與試驗子結(jié)構(gòu)相同的數(shù)值子結(jié)構(gòu)恢復(fù)力模型參數(shù)。通過數(shù)值仿真研究該方法的精度。

1 隱性卡爾曼濾波器的基本原理

1995年，Julier和Uhlmann基于經(jīng)典卡爾曼濾波器提出了隱性卡爾曼濾波器（Unscented Kalman Filter）［4］，并將它應(yīng)用于非線性系統(tǒng)狀態(tài)估計中［5］。隱性卡爾曼濾波器假定系統(tǒng)的狀態(tài)量及噪聲符合高斯分布，采用確定性采樣方法產(chǎn)生第k步狀態(tài)量的樣本點（σ點），然后由系統(tǒng)狀態(tài)方程得到k＋1步狀態(tài)量的樣本點，通過統(tǒng)計方法計算得到第k＋1步狀態(tài)量的均值、方差和協(xié)方差。這個計算過程就稱之為隱性變換（Unscented Transformation）［5］。確定k＋1步狀態(tài)量的均值、協(xié)方差后，通過經(jīng)典卡爾曼濾波器更新方程就可求得狀態(tài)量的后驗條件概率分布均值及協(xié)方差，其中條件均值即k＋1步狀態(tài)量的估計值［6］，具體過程如下。

設(shè)所考慮的問題可寫成如下離散形式的狀態(tài)方程和觀測方程

式中xk∈Rn為n維隨機狀態(tài)向量；uk為系統(tǒng)輸入；yk＋1∈Rm為m維隨機觀測向量；vk～N（0，Qk）為過程高斯白噪聲，wk＋1～N（0，Rk＋1）為觀測高斯白噪聲，Qk和Rk＋1分別為過程噪聲和觀測噪聲協(xié)方差矩陣。隱性卡爾曼濾波算法分為預(yù)測步和更新步，具體步驟如下［6］：

（1）預(yù)測步：即在不考慮第k＋1步觀測值的情況下，根據(jù)狀態(tài)方程通過UT變換得到第k＋1步狀態(tài)xk＋1估計的均值mk－和協(xié)方差Pk－。首先按下式確定第k步系統(tǒng)狀態(tài)xk的2n＋1個樣本點

然后按式（5），（6）計算預(yù)測狀態(tài)的均值和協(xié)方差［7］

式（5），（6）中的權(quán)重向量Wm及權(quán)重矩陣W由下式確定

常數(shù)α，κ為該方法中與狀態(tài)概率分布類型有關(guān)的參數(shù)，α決定了樣本點在均值周圍的分布［8］。UT變換可使自變量和因變量的均值和方差至少達到泰勒展式的二階精度［9］。對于不同的概率分布，權(quán)值和λ合理的取值，均值和方差的精度可以達到更高。α，β和κ的具體定義及意義參見文獻［8］。當協(xié)方差權(quán)值均取非負數(shù)時，可以保證狀態(tài)量的協(xié)方差矩陣為半正定矩陣。

（2）更新步：根據(jù)觀測方程通過UT變換得到觀測量的均值、協(xié)方差矩陣，以及觀測量與狀態(tài)量的互協(xié)方差矩陣；然后利用Kalman更新方程計算k＋1步狀態(tài)估計均值和協(xié)方差。

計算濾波增益Kk＋1以及更新的狀態(tài)估計均值mk＋1和協(xié)方差Pk＋1

式中yk＋1為第k＋1步的觀測。

2 雙折線恢復(fù)力模型參數(shù)在線識別

雙折線恢復(fù)力模型增量形式可表示成下式

式中Rk為恢復(fù)力；f（θk，uk，tk）為結(jié)構(gòu)恢復(fù)力增量函數(shù)，與模型參數(shù)和加載歷程有關(guān)；θk為恢復(fù)力模型的 參 數(shù) 向 量，θk＝ ［k1，kk2，kRy，k］T，式 中k1，k，k2，k和Ry，k分別為結(jié)構(gòu)第一、第二剛度及 屈 服力；f（θk，uk，tk）＝kkuk＝kk［x（tk＋Δt）－x（tk）］，式中kk應(yīng)根據(jù)滯回階段取k1，k或k2，k，而滯回階段不僅與位移有關(guān)，還與屈服力有關(guān)。

下面以一算例說明隱性卡爾曼濾波方法的識別效果。真實結(jié)構(gòu)恢復(fù)力模型參數(shù)分別為：第一剛度40 000kN／m，第二剛度1 000kN／m，屈服力400 kN。對結(jié)構(gòu)進行位移控制加載，作用于結(jié)構(gòu)的位移命令選用El Centro（1940，NS）的地面運動位移記錄，位移峰值為5cm。濾波器參數(shù)設(shè)置為α＝1，β＝0，κ＝0。利用隱性卡爾曼濾波在線識別雙折線恢復(fù)力模型參數(shù)時，設(shè)系統(tǒng)狀態(tài)為xk＝［Rkk1，kk2，k Ry，k］T，結(jié)構(gòu)的狀態(tài)方程和觀測方程分別為

噪聲協(xié)方差矩陣，結(jié)構(gòu)狀態(tài)量初始狀態(tài)的均值m0和協(xié)方差矩陣P0如下：

采用Intel雙核CPU（2．33GHz）的普通臺式計算機，數(shù)值仿真耗時6．20s，結(jié)果如圖1～5所示。從圖1和2的仿真結(jié)果可以看出，在線識別預(yù)測得到的結(jié)構(gòu)恢復(fù)力與真實值十分接近，兩者幾乎重合。從圖3～5的結(jié)構(gòu)恢復(fù)力模型參數(shù)識別時程曲線可以看出，隱性卡爾曼濾波可以快速地識別出結(jié)構(gòu)的恢復(fù)力模型參數(shù)，所識別出的結(jié)果與真實值相差均不超過1%，具有很高的精度。另外，在仿真時觀測量恢復(fù)力中加入了觀測噪聲，證明該方法對折線型恢復(fù)力模型參數(shù)在線識別也具有很好的魯棒性。經(jīng)仿真分析發(fā)現(xiàn)，隱性卡爾曼濾波對于狀態(tài)方程初始狀態(tài)量的協(xié)方差矩陣值比較敏感。當初始狀態(tài)量的協(xié)方差矩陣取為主對角矩陣時，隱性卡爾曼濾波相應(yīng)參數(shù)識別效果主要與初始狀態(tài)量協(xié)方差矩陣中對應(yīng)主元素絕對值有關(guān)，與其他值關(guān)聯(lián)不大。

圖1 結(jié)構(gòu)滯回曲線Fig．1 Structure hysteretic curves

圖2 恢復(fù)力時程曲線Fig．2 Time history curves of restoring force

圖3 第一剛度識別時程曲線Fig．3 Time history curves of the first stiffness identification

圖4 第二剛度識別時程曲線Fig．4 Time history curves of the second stiffness identification

圖5 屈服力識別時程曲線Fig．5 Time history curves of the yield force identification

3 基于隱性卡爾曼濾波的自適應(yīng)子結(jié)構(gòu)擬動力試驗方法

從第2節(jié)算例可以看出，隱性卡爾曼濾波能快速準確地識別出結(jié)構(gòu)的恢復(fù)力模型參數(shù)。因此本文將隱性卡爾曼濾波應(yīng)用于子結(jié)構(gòu)試驗中，提出了基于隱性卡爾曼濾波器自適應(yīng)子結(jié)構(gòu)試驗方法。該方法的具體步驟如下所示：

（1）初始化。建立試驗子結(jié)構(gòu)恢復(fù)力參數(shù)識別的狀態(tài)方程及觀測方程，根據(jù)經(jīng)驗確定系統(tǒng)噪聲、觀測噪聲協(xié)方差矩陣，結(jié)構(gòu)系統(tǒng)初始狀態(tài)量的均值及協(xié)方差矩陣。

（2）利用數(shù)值積分方法計算出數(shù)值子結(jié)構(gòu)的位移，把位移命令發(fā)送給作動器。

（3）對試驗子結(jié)構(gòu)進行加載，將測得的觀測量（恢復(fù)力）傳遞給結(jié)構(gòu)模型識別模塊。

（4）利用隱性卡爾曼濾波及測得的觀測量對試驗子結(jié)構(gòu)恢復(fù)力模型參數(shù)進行在線識別，根據(jù)識別結(jié)果更新數(shù)值子結(jié)構(gòu)相應(yīng)恢復(fù)力模型參數(shù)。

（5）轉(zhuǎn)至第（2）步，直至試驗結(jié)束。

下面對圖6所示兩層帶支撐框架結(jié)構(gòu)進行自適應(yīng)子結(jié)構(gòu)擬動力試驗的數(shù)值仿真，檢驗基于隱性卡爾曼濾波子結(jié)構(gòu)試驗方法的性能。

圖中下層層間支撐作為試驗子結(jié)構(gòu)，其他部位作為數(shù)值子結(jié)構(gòu)?？蚣芟聦蛹吧蠈咏Y(jié)構(gòu)的質(zhì)量分別為MN1＝MN2＝2 000t，剛度分別為KN1＝KN2＝80 000kN／m，阻尼系數(shù)分別為CN1＝CN2＝1 550 kN／（m／s）。上下層結(jié)構(gòu)層間支撐均采用 Bouc－Wen模型，其表達式如下

假定下層試驗子結(jié)構(gòu)支撐模型參數(shù)取值為kb＝40 000kN／m，β＝60，γ＝40，n＝1．1。支撐與樓面的夾角均為28．81°。地震動選取 El Centro（1940，NS）地震記錄，地震動峰值加速度為4m／s2。數(shù)值積分算法采用4階Runge－Kutta方法，積分時間步長為0．01s。

在自適應(yīng)試驗中，利用UKF識別試驗子結(jié)構(gòu)支撐的剛度kb以及滯回參數(shù)β，γ和n，并實時更新數(shù)值子結(jié)構(gòu)支撐模型參數(shù)。觀測量選用試驗子結(jié)構(gòu)的恢復(fù)力。試驗子結(jié)構(gòu)的狀態(tài)方程和觀測方程分別為下式所示

式中˙x′為試驗支撐加載速度，第k步加載速度˙x′k可以根據(jù)試驗子結(jié)構(gòu)位移加載命令差分得到。試驗子結(jié)構(gòu)狀態(tài)量的初始均值m0、協(xié)方差矩陣P0和過程噪聲的方差矩陣Qk如下式所示。

圖6 支撐結(jié)構(gòu)子結(jié)構(gòu)試驗示意圖Fig．6 Schematic of the substructure testing for the frame－brace structure

在觀測向量中加入?yún)f(xié)方差為恢復(fù)力真實值標準差1%的觀測高斯白噪聲。由于位移相關(guān)型子結(jié)構(gòu)試驗中的試驗子結(jié)構(gòu)的狀態(tài)方程中不存在阻尼，且狀態(tài)量之間的聯(lián)系相對較差，因此抗噪能力相對較差。當觀測量加入較大噪聲時有時會出現(xiàn)計算發(fā)散。自適應(yīng)子結(jié)構(gòu)試驗仿真耗時7．01s，試驗仿真結(jié)果如圖7～14所示。

圖7 下層結(jié)構(gòu)層間位移時程曲線Fig．7 Time history curves of the lower structure interstory displacement

圖8 上層結(jié)構(gòu)層間位移時程曲線Fig．8 Time history curves of the upper structure interstory displacement

傳統(tǒng)子結(jié)構(gòu)試驗數(shù)值子結(jié)構(gòu)中支撐恢復(fù)力模型采用Bouc－Wen模型，模型參數(shù)采用自適應(yīng)子結(jié)構(gòu)試驗中試驗子結(jié)構(gòu)模型參數(shù)估計初始值，分別為kb＝50 000kN／m，β＝40，γ＝50，n＝2。從圖7～10可以看出，位移相關(guān)型自適應(yīng)子結(jié)構(gòu)試驗結(jié)果與結(jié)構(gòu)真實反應(yīng)基本吻合，較傳統(tǒng)子結(jié)構(gòu)試驗結(jié)果有很大改善。由于數(shù)值子結(jié)構(gòu)中支撐恢復(fù)力模型的初始參數(shù)與實際值差別較大，在傳統(tǒng)子結(jié)構(gòu)試驗仿真中，數(shù)值子結(jié)構(gòu)仍基本處于線性段，未進入非線性。從圖11～14可以看出自適應(yīng)子結(jié)構(gòu)試驗可以很快地準確地識別出結(jié)構(gòu)剛度k和恢復(fù)力模型參數(shù)n，恢復(fù)力模型參數(shù)β和γ也很快趨近于真實值。雖然β和γ最終值與真實值有一定的差異，但對試驗精度并無太大影響。

圖10 上層支撐滯回曲線Fig．10 Hysteretic curves of upper brace

圖11 剛度識別時程曲線Fig．11 Time history curves of the stiffness identification

圖12 β識別時程曲線Fig．12 Time history curves ofβidentification

圖13 γ識別時程曲線Fig．13 Time history curves ofγidentification

圖14 n識別時程曲線Fig．14 Time history curves of nidentification

4 結(jié) 論

本文利用隱性卡爾曼濾波在線識別結(jié)構(gòu)的恢復(fù)力模型參數(shù)，并提出了基于隱性卡爾曼濾波在線識別的自適應(yīng)子結(jié)構(gòu)試驗方法，對兩層框架支撐結(jié)構(gòu)進行了試驗仿真分析。結(jié)果表明基于隱性卡爾曼濾波的恢復(fù)力模型參數(shù)在線識別方法具有較好的精度和收斂速度，且耗時較短?；陔[性卡爾曼濾波的自適應(yīng)子結(jié)構(gòu)試驗方法與傳統(tǒng)子結(jié)構(gòu)試驗方法相比，其結(jié)構(gòu)反應(yīng)更接近真實值。

［1］ 邱法維，錢稼茹，陳志鵬．結(jié)構(gòu)抗震實驗方法［M］．北京：科學(xué)出版社，2000：1—153．

Qiu F W，Qian J R，Chen Z P．Seismic Test Methodology for Structures［M］．Beijing：Science Press，2000：1—153．

［2］ Yang W J，Nakano Y．Substructure online test by using real－time hysteresis modeling with a neural network［A］．Proceedings of the First International Conference on Advances in Experimental Structural Engineering［C］．2005：267—274．

［3］ 張?。赃m應(yīng)子結(jié)構(gòu)擬動力試驗方法［D］．哈爾濱：哈爾濱工業(yè)大學(xué)，2010．

Zhang J．Adaptive substructure pseudo－dynamic testing method［D］．Harbin：School of Civil Engineering，Harbin Institute of Technology，2010．

［4］ Julier S J，Uhlmann J，Durrant－Whyte H F．A new approach for filtering nonlinear systems［A］．Ameri－can Control Conference［C］．Seattle， Washington，1995：1 628—1 632．

［5］ Julier S J，Uhlmann J．Unscented filtering and nonlinear estimation ［J］．Proceedings of The IEEE，2004，92（3）：401—422．

［6］ Julier S J，Uhlmann J．A general method for approximation nonlinear transformation of probability distributions［R］．Technical Report，Robotics Research Group，Department of Engineering Science，University of Oxford，1996：1—27．

［7］ Xiong K，Zhang H Y，Chan C W．Performance evalu－ation of UKF－based nonlinear filtering［J］．Automatica，2006，42（2）：261—270．

［8］ Merwe R．Sigma－point Kalman filters for probabilistic inference in dynamic state－space models［D］．Oregon Health ＆Science University，2004．

［9］ Merwe R，Wan E A．The square－root unscented Kalman filter for state and parameter－estimation［A］．Proc．of the International Conference on Acoustics，Speech and Signal Processing，IEEE［C］．Salt Lake City，2001：3 461—3 464．