基于改進(jìn)卡爾曼濾波的結(jié)構(gòu)損傷識別方法研究

沈 皓，常 軍

（蘇州科技大學(xué) 土木工程學(xué)院，江蘇 蘇州 215011）

卡爾曼濾波算法是近幾十年發(fā)展起來的一種行之有效的結(jié)構(gòu)損傷識別方法。卡爾曼濾波算法是一種基于最小方差準(zhǔn)則的線性無偏估計(jì)方法，能夠?qū)崿F(xiàn)有限觀測量下的結(jié)構(gòu)狀態(tài)的實(shí)時(shí)估計(jì)[1]。對于土木工程結(jié)構(gòu)的卡爾曼濾波算法損傷識別，國內(nèi)外學(xué)者已經(jīng)做過很多研究。王堅(jiān)等[2]提出了擴(kuò)展卡爾曼濾波，解決了一般的線性問題，但原始狀態(tài)誤差較大或非線性程度較高時(shí)，會(huì)嚴(yán)重影響濾波精度；寶志雯等[3-4]提出了基于推廣卡爾曼濾波算法的剪切型結(jié)構(gòu)模型的參數(shù)識別方法，實(shí)現(xiàn)了結(jié)構(gòu)剛度和阻尼的識別。周麗等[5]提出了基于自適應(yīng)卡爾曼濾波的識別方法，并對剪切線性模型實(shí)現(xiàn)了剛度和阻尼的識別。

在系統(tǒng)識別方面，已經(jīng)有很多時(shí)域分析方法在相應(yīng)文獻(xiàn)中提出，比如：最小二乘估計(jì)[6]，卡爾曼濾波[7]，序列非線性最小二乘法等[8]。Julier[9]提出了無跡卡爾曼濾波，采用Sigma點(diǎn)近似其均值與方差，簡化了計(jì)算，但帶來的問題是對高于3維狀態(tài)向量估計(jì)時(shí)，其精度達(dá)不到要求；王路[10]提出了自適應(yīng)迭代無跡卡爾曼濾波，其實(shí)現(xiàn)了自動(dòng)迭代的功能，不需要人工設(shè)置迭代，可還是不能解決高維的非線性問題，精度不夠理想；Arasaratnam等[11]提出的容積卡爾曼濾波（Cubature Kalman Filter，CKF），其核心思想是對于非線性高斯系統(tǒng)，利用數(shù)值積分理論，保證任意非線性高斯?fàn)顟B(tài)的后驗(yàn)均值和方差可以近似為三階多項(xiàng)式，但由于初始誤差較大，降低了CKF對其估計(jì)的精度；為了提高CKF運(yùn)行速度，穆靜[12]提出了迭代容積卡爾曼濾波（Iterated CKF，CKF），但其不具有自動(dòng)迭代的功能，需要人工設(shè)定迭代的次數(shù)。為此，本文借鑒遺忘因子的濾波算法和AIUKF的思想，把遺忘因子與自適應(yīng)迭代容積卡爾曼濾波（Adaptive Iteration CKF，AICKF）相結(jié)合，這樣既可以起到遺忘因子的作用，減少歷史數(shù)據(jù)對濾波結(jié)果的影響，又可以提高濾波算法本身的準(zhǔn)確性和處理非線性問題的能力。

1 KF算法的結(jié)構(gòu)損傷識別

1.1 離散狀態(tài)空間方程和KF算法

線性系統(tǒng)的結(jié)構(gòu)動(dòng)力學(xué)方程可寫成連續(xù)的狀態(tài)空間方程。如式（1）和（2）所示。

KF算法是無偏移遞歸算法，可最優(yōu)估計(jì)未知狀態(tài)向量。經(jīng)典KF算法為第k步時(shí)估計(jì)先驗(yàn)狀態(tài)向量

第k步時(shí)刻的狀態(tài)向量

先驗(yàn)誤差方差

狀態(tài)向量估計(jì)均方差

卡爾曼濾波增益

在經(jīng)典KF算法中，非線性動(dòng)態(tài)系統(tǒng)中過程噪聲協(xié)方差矩陣Q，在實(shí)際中往往很難準(zhǔn)確獲得，而導(dǎo)致濾波器發(fā)散；濾波算法本身精度不高，會(huì)降低處理問題的能力。因此本文借鑒遺忘算法因子的濾波算法與AICKF相結(jié)合，這樣既可以發(fā)揮遺忘因子的作用，減少歷史數(shù)據(jù)對濾波結(jié)果的影響避免發(fā)散，又可以提高濾波精度和處理非線性問題的能力。

2 改進(jìn)后KF算法結(jié)構(gòu)損傷識別

2.1 AICKF算法

借鑒AIUKF算法的思想，在CKF基礎(chǔ)上加入自適應(yīng)迭代，在濾波算法上CKF精度就高于UKF，在加入自適應(yīng)迭代后一方面迭代再次提高了濾波的精度，另一方面自適應(yīng)省去了人工去設(shè)置迭代的問題。

考慮如下離散時(shí)間非線性動(dòng)態(tài)系統(tǒng)

式中，xk為系統(tǒng)狀態(tài)向量；zk為量測值。假設(shè)過程噪聲wk-1和量測噪聲vk相互獨(dú)立，是均值為零，協(xié)方差分別為Qk-1和Rk的高斯白噪聲。

容積卡爾曼濾波算法首先計(jì)算加權(quán)函數(shù)為標(biāo)準(zhǔn)正態(tài)分布密度的積分的基本容積點(diǎn)和對應(yīng)的權(quán)值

m為系統(tǒng)的狀態(tài)維數(shù)；ζj為第j個(gè)容積點(diǎn)；基本容積點(diǎn)按照下列方式產(chǎn)生，記n維單位向量為e=[1,0,…,0]T，使用[1]表示對的元素進(jìn)行全排列和改變元素符號產(chǎn)生的點(diǎn)集，稱為完整全對稱點(diǎn)集，[1]j表示點(diǎn)集中[1]的第j個(gè)點(diǎn)；wj為對應(yīng)點(diǎn)的權(quán)值。

2.2 遺忘因子的計(jì)算

遺忘因子[13-14]是為了限制濾波器的長度，讓濾波值中的歷史數(shù)據(jù)占的權(quán)重變小，新數(shù)據(jù)權(quán)重變大。

2.3 遺忘因子和AICKF結(jié)合的濾波算法

該算法的不同之處在于方差預(yù)測加入遺忘因子，即

這就意味著隨著過程噪聲協(xié)方差矩陣不準(zhǔn)確時(shí)，加入遺忘因子不會(huì)對預(yù)測協(xié)方差產(chǎn)生較大影響，遺忘因子降低了歷史數(shù)據(jù)對濾波結(jié)果的影響，增加了當(dāng)前測量值在估計(jì)中的比例，提高了濾波精度[15-16]。具體算法步驟如下。

（1）初始化。初始化公式

（2）時(shí)間更新。首先，計(jì)算容積點(diǎn)及經(jīng)非線性狀態(tài)方程傳播的容積點(diǎn)

式中：Sk-1=chol（Pk-1），chol（）表示矩陣的喬列斯基分解。

然后，計(jì)算狀態(tài)和方差預(yù)測，見下式

（3）計(jì)算容積點(diǎn)經(jīng)量測方程的傳遞值和測量值。見下式。

（4）判斷是否迭代。計(jì)算出預(yù)測值與實(shí)際觀測值的適應(yīng)度函數(shù)，體積點(diǎn)的傳遞值與實(shí)際觀測值的適應(yīng)度函數(shù)，并根據(jù)適應(yīng)度函數(shù)的比值確定采樣點(diǎn)與目標(biāo)真實(shí)估計(jì)值的偏差，從而自適應(yīng)確定是否進(jìn)行迭代重采樣。具體步驟，首先是定義適應(yīng)度函數(shù)，預(yù)測測量值與實(shí)際觀測值的適應(yīng)度函數(shù)為

容積點(diǎn)傳遞值與實(shí)際觀測值的適度函數(shù)為

適應(yīng)度函數(shù)比

其中：zk為實(shí)際觀測值；z^k為預(yù)測測量值；Zj，k為容積點(diǎn)傳遞值；Rk為觀測噪聲方差。

然后，判斷是否迭代數(shù)據(jù)。

（5）若ρ＜1，表示采樣點(diǎn)有效逼近真實(shí)估計(jì)，則不進(jìn)行迭代；接著往下計(jì)算量測更新。首先，計(jì)算新息方差和協(xié)方差

然后，計(jì)算增益。見下式

最后，計(jì)算系統(tǒng)的狀態(tài)更新和協(xié)方差更新。見下式

（6）若ρ≥1，表示采樣點(diǎn)與真實(shí)估計(jì)偏差較大，則進(jìn)行迭代，步驟如下。

第二，重新時(shí)間更新。

第四，重新量測更新。

設(shè)迭代終止時(shí)的迭代次數(shù)為N，則k時(shí)刻的狀態(tài)估計(jì)與協(xié)方差估計(jì)為

3 數(shù)值模擬分析

采用簡支梁為數(shù)值算例驗(yàn)證所提方法的有效性和可靠性。該簡支梁共20個(gè)單元，每個(gè)單元長1 m，梁截面尺寸0.5 m×1 m，結(jié)構(gòu)的彈性模量2×107Pa，泊松比0.2，密度2 500 kg/m3。利用有限元分析軟件ABAQUS建立了相應(yīng)的有限元模型。通過降低單元?jiǎng)偠葋砟M結(jié)構(gòu)損傷狀況，簡支梁模型如圖1所示。

圖1 簡支梁模型

對各節(jié)點(diǎn)輸入不同的高斯白噪聲，模擬環(huán)境激勵(lì)。激勵(lì)頻率為500 Hz，激勵(lì)時(shí)長20 s。每個(gè)單元的初始剛度為103 136 N/m。首先模擬單元3處發(fā)生損傷，剛度損傷5%，損傷后剛度值為97 980 N/m。考慮噪聲影響，以白噪聲激勵(lì)模擬隨機(jī)荷載，作用各節(jié)點(diǎn)，取噪聲水平5%。圖2是以上述工況為例，得到3號節(jié)點(diǎn)的加速度響應(yīng)信號；圖3是3號單元KF剛度識別曲線；圖4是改進(jìn)后KF剛度識別曲線。

圖2 3號節(jié)點(diǎn)的加速度響應(yīng)

從圖3可看出，簡支梁模型發(fā)生損傷時(shí)，傳統(tǒng)卡爾曼濾波算法對損傷狀況較為敏感，能夠有效識別其損傷。圖4表明改進(jìn)后的卡爾曼濾波識別損傷結(jié)果。通過圖3和圖4的比較能夠很直觀的看出改進(jìn)后的卡爾曼濾波識別值能夠快速且穩(wěn)定地收斂到真實(shí)值。表1是選取簡支梁的部分單元參數(shù)識別結(jié)果對比?？梢钥闯龈倪M(jìn)后的卡爾曼濾波誤差值普遍降低了2%左右，說明改進(jìn)后的卡爾曼濾波更能夠準(zhǔn)確有效的識別損傷。

表1 簡支梁參數(shù)識別結(jié)果對比

圖3 3號單元KF剛度識別曲線

圖4 3號單元改進(jìn)后KF剛度識別曲線

圖5是噪聲下剛度收斂時(shí)程曲線。從圖5可以看出，改進(jìn)后的卡爾曼濾波相比于傳統(tǒng)的卡爾曼濾波的剛度識別曲線能夠更快，更穩(wěn)定的收斂到真實(shí)值。表2可以看出改進(jìn)后的卡爾曼濾波算法在5%噪聲水平時(shí)，剛度參數(shù)識別誤差為1%左右，相比于傳統(tǒng)卡爾曼濾波誤差降低了5%，說明改進(jìn)后的卡爾曼濾波對于結(jié)構(gòu)損傷識別具有較強(qiáng)的抗噪性能。

表2 噪聲下簡支梁參數(shù)識別結(jié)果對比

圖5 噪聲下剛度收斂時(shí)程曲線

為了研究改進(jìn)后的卡爾曼濾波不僅能識別單點(diǎn)損傷，還能識別多點(diǎn)損傷，假定單元2和單元4模擬損傷，模擬剛度下降10%，剛度值和為92 822 N/m，取噪聲水平為5%。

表3是多點(diǎn)損傷簡支梁結(jié)構(gòu)參數(shù)識別結(jié)果。從表3可以不難看出，簡支梁模型發(fā)生兩點(diǎn)損傷時(shí)，改進(jìn)后的卡爾曼濾波識別精度相比于傳統(tǒng)的卡爾曼濾波識別能力更強(qiáng)，它能夠有效地提取結(jié)構(gòu)的損傷信號，進(jìn)而準(zhǔn)確的識別結(jié)構(gòu)損傷狀況，說明改進(jìn)后的卡爾曼濾波具有多點(diǎn)損傷識別的能力。

表3 多點(diǎn)損傷簡支梁結(jié)構(gòu)參數(shù)識別結(jié)果

4 結(jié)論

（1）新改進(jìn)的算法與傳統(tǒng)的卡爾曼濾波算法相比較，通過數(shù)值模擬的簡支梁進(jìn)行參數(shù)識別，結(jié)果表明能夠更加精確的識別結(jié)構(gòu)構(gòu)件的剛度進(jìn)而診斷結(jié)構(gòu)損傷位置和程度效果。

（2）通過上述數(shù)據(jù)結(jié)果表明該方法還可以有效減少誤差，提高濾波精度。數(shù)值模擬結(jié)果表明該方法可以排除噪聲因素的干擾，識別結(jié)構(gòu)損傷發(fā)生的部位和程度，具有較好的抗噪性。