擴(kuò)展卡爾曼濾波對(duì)時(shí)變參數(shù)追蹤性能的影響研究

唐宏志，姜金輝

（南京航空航天大學(xué)機(jī)械結(jié)構(gòu)力學(xué)及控制國家重點(diǎn)實(shí)驗(yàn)室，南京 210016）

參數(shù)識(shí)別是動(dòng)力學(xué)的第一類逆問題。參數(shù)識(shí)別問題按識(shí)別參數(shù)的種類可以分為模態(tài)參數(shù)識(shí)別與物理參數(shù)識(shí)別；而如果按識(shí)別參數(shù)的方法，可以分為時(shí)域法參數(shù)識(shí)別和頻域法［1?2］參數(shù)識(shí)別。在時(shí)域內(nèi)識(shí)別結(jié)構(gòu)的物理參數(shù)，主要方法有擴(kuò)展卡爾曼濾波（Extended Kalman filter，EKF）、H∞濾波［3］、蒙特卡洛濾波［4］和最小二乘估計(jì)（Least square esti?mate，LSE）［5?6］。EKF 是在經(jīng)典卡爾曼濾波（Kal?man filter，KF）［7?8］的基礎(chǔ)上改進(jìn)得到的，文獻(xiàn)［9?11］使用擴(kuò)展卡爾曼濾波方法對(duì)多自由度線性結(jié)構(gòu)的參數(shù)進(jìn)行識(shí)別。文獻(xiàn)［12?13］為了實(shí)現(xiàn)對(duì)結(jié)構(gòu)時(shí)變參數(shù)的跟蹤，采用自適應(yīng)技術(shù)對(duì)EKF 方法進(jìn)行了改進(jìn)。文獻(xiàn)［14?15］在載荷未知的情況下，將最小二乘法與EKF 方法結(jié)合，識(shí)別時(shí)變的結(jié)構(gòu)參數(shù)，同時(shí)求出未知載荷。

這些文獻(xiàn)都對(duì)擴(kuò)展卡爾曼濾波方法進(jìn)行了改進(jìn)，但并沒有提及濾波參數(shù)的設(shè)置對(duì)識(shí)別效果的影響，尤其是結(jié)構(gòu)中參數(shù)會(huì)隨時(shí)間變化的情況。事實(shí)上，隨著算法的迭代計(jì)算，EKF 算法對(duì)測點(diǎn)響應(yīng)的利用程度會(huì)越來越小，此時(shí)識(shí)別結(jié)果無法追蹤到時(shí)變參數(shù)的變化情況。因此，必須選擇恰當(dāng)?shù)臑V波參數(shù)（模型噪聲的協(xié)方差矩陣Q，測量噪聲的協(xié)方差矩陣R）。本文首先介紹了一種結(jié)構(gòu)參數(shù)識(shí)別方法，適用于外界載荷已知的情況，并試圖通過仿真算例說明如何選擇最優(yōu)的濾波參數(shù)，來使得EKF算法對(duì)時(shí)變結(jié)構(gòu)參數(shù)的跟蹤性能達(dá)到最佳，即算法能快速、準(zhǔn)確地識(shí)別出結(jié)構(gòu)參數(shù)的變化情況。

1 結(jié)構(gòu)參數(shù)識(shí)別方法

本文使用擴(kuò)展卡爾曼濾波方法作為結(jié)構(gòu)參數(shù)的識(shí)別方法。對(duì)于1 個(gè)n自由度的動(dòng)力學(xué)系統(tǒng)，當(dāng)結(jié)構(gòu)中的未知參數(shù)，例如質(zhì)量矩陣M，剛度矩陣K或阻尼矩陣C中的某些參數(shù)待確定時(shí)，參考卡爾曼濾波理論，將這些未知參數(shù)與狀態(tài)向量x一起組成擴(kuò)展?fàn)顟B(tài)向量z

式中：θ為待確定的未知參數(shù)向量，長度設(shè)為m；x(t)表示狀態(tài)向量，其長度為2n。令θ?=0，此時(shí)，卡爾曼濾波理論中的狀態(tài)更新方程可以重新寫為

而測量更新方程的形式與觀測量的類別有關(guān)，本文選擇的測量類別是加速度信號(hào)。此時(shí)，測量更新方程可以寫為

式中：M表示大小為n×n的質(zhì)量矩陣；C表示大小為n×n的阻尼矩陣；K表示大小為的n×n的剛度矩陣；p?(t)，p?(t)，p(t)分別代表結(jié)構(gòu)的加速度向量，速度向量，以及位移向量，向量的長度為n；Bu代表激勵(lì)的影響矩陣，與激勵(lì)在結(jié)構(gòu)上的作用位置有關(guān)，由0 和1 組成，它的矩陣大小為s×n；f(t)表示激勵(lì)向量，長度為s；H0為由0 和1 組成的位置矩陣。由于質(zhì)量矩陣M，剛度矩陣K或阻尼矩陣C包含了未知參數(shù)θ，這也就意味著關(guān)于擴(kuò)展?fàn)顟B(tài)向量z(t)的狀態(tài)更新方程fc(z(t)，f(t))是一個(gè)非線性方程，同樣測量更新方程h(z(t)，f(t))也是一個(gè)非線性方程。

利用泰勒展開，將上述兩個(gè)非線性方程（1，2）轉(zhuǎn)化為多項(xiàng)式之和的形式，并略去高階項(xiàng)，僅保留一階多項(xiàng)式，由此可將非線性方程近似看成線性方程

式中：zk-1|k-1為在tk-1時(shí)刻的擴(kuò)展?fàn)顟B(tài)向量；zk|k-1為在tk時(shí)刻的擴(kuò)展?fàn)顟B(tài)向量；?z fc k-1為狀態(tài)更新方程對(duì)擴(kuò)展?fàn)顟B(tài)向量z(t)的一階偏導(dǎo)；?zhk為測量更新方程對(duì)擴(kuò)展?fàn)顟B(tài)向量z(t)的一階偏導(dǎo)，又稱為雅各比矩陣。由泰勒展開定理可知，只要t時(shí)刻的擴(kuò)展?fàn)顟B(tài)向量z(t)與zk-1|k-1，zk|k-1之間大小差距很小，式（4，5）就能近似成立。而要使z(t)與zk-1|k-1，zk|k-1的差距很小，只要時(shí)刻t與時(shí)刻tk-1，tk差距很小即可。

得到離散化的近似線性等式（4，5）后，與卡爾曼濾波理論的方法類似。首先是狀態(tài)更新過程，其先驗(yàn)估計(jì)的均值和估計(jì)誤差的方差可以寫為

由式（7）可知，擴(kuò)展?fàn)顟B(tài)向量中的參數(shù)θ在狀態(tài)更新過程中并不會(huì)變化，僅憑該過程并不能實(shí)現(xiàn)時(shí)變參數(shù)的追蹤。然后對(duì)于測量更新方程，后驗(yàn)估計(jì)的均值和估計(jì)誤差的方差可以寫為

式中：Kk為卡爾曼濾波增益，它表示對(duì)測量響應(yīng)的利用程度。通過方程（8，9）可以使得擴(kuò)展?fàn)顟B(tài)向量zk|k中的參數(shù)θ發(fā)生變化。因此，可以認(rèn)為擴(kuò)展卡爾曼濾波算法是通過測量更新過程實(shí)現(xiàn)了對(duì)參數(shù)的追蹤，只要加大對(duì)測量響應(yīng)的利用程度，就能使zk|k較快地收斂到真實(shí)參數(shù)值。Kk其求法與卡爾曼理論中的相同，可以寫為

由此，可以得到擴(kuò)展卡爾曼濾波算法的流程如表1 所示。

表1 擴(kuò)展卡爾曼濾波法算法流程Table 1 Algorithm flow of extended Kalman filter

2 EKF 算法追蹤結(jié)構(gòu)時(shí)變參數(shù)的影響因素

在前文狀態(tài)更新方程的推導(dǎo)中，假定未知參數(shù)θ不會(huì)隨時(shí)間變化，即θ?=0。這也就意味著僅憑狀態(tài)更新方程，只能識(shí)別出時(shí)不變的參數(shù)，對(duì)時(shí)變參數(shù)或緩變參數(shù)是無法識(shí)別出來的。對(duì)于變化參數(shù)的追蹤，是通過測量更新方程來實(shí)現(xiàn)的?？柭鼮V波增益Kk表示了對(duì)測量更新方程的利用程度。若希望能實(shí)現(xiàn)對(duì)時(shí)變參數(shù)的追蹤，則對(duì)測量更新方程的利用程度必須加大，總的宗旨是能盡量利用測量更新方程，即令Kk在迭代過程中仍保持較大值。

由卡爾曼濾波增益的計(jì)算公式可知，卡爾曼增益Kk由測量噪聲方差R和Pk|k-1共同決定，初始時(shí)刻Kk偏大，對(duì)狀態(tài)的估計(jì)主要依賴于測量更新方程。隨著迭代過程的進(jìn)行，測量對(duì)狀態(tài)的修正作用不斷減少，系統(tǒng)反而依賴于狀態(tài)更新方程，此時(shí)Kk會(huì)逐漸變小。

為了使得濾波過程中，一直保持對(duì)測量更新方程的利用，需要選擇合適的濾波參數(shù)。比如：模型噪聲的協(xié)方差矩陣Q，對(duì)應(yīng)狀態(tài)向量部分的取值較小。對(duì)應(yīng)參數(shù)部分的取值較大，這相當(dāng)于告知濾波系統(tǒng)對(duì)結(jié)構(gòu)未知參數(shù)的估計(jì)存在較大偏差，促使算法識(shí)別的結(jié)構(gòu)參數(shù)會(huì)向真實(shí)值靠攏；測量噪聲的協(xié)方差矩陣R的取值偏小，這相當(dāng)于承認(rèn)濾波系統(tǒng)所測量的響應(yīng)具有很高的置信度，加大系統(tǒng)對(duì)測量更新方程的利用程度，從而加強(qiáng)算法對(duì)結(jié)構(gòu)時(shí)變參數(shù)的追蹤能力；然后是外界的測量噪聲對(duì)參數(shù)識(shí)別效果的影響，算法能在響應(yīng)存在噪聲干擾的情況下識(shí)別出結(jié)構(gòu)未知參數(shù)，這是通過調(diào)整測量噪聲協(xié)方差矩陣R的方式來實(shí)現(xiàn)的；最后是觀測點(diǎn)數(shù)的影響，選擇懸臂梁作為仿真算例，由于懸臂梁既有轉(zhuǎn)角自由度，又有豎向位移自由度，選擇僅僅觀測豎向位移自由度可以看出算法仍能較好地識(shí)別出結(jié)構(gòu)參數(shù)變化的情況。

3 仿真算例

本文選擇一個(gè)三自由度系統(tǒng)，如圖1 所示，討論模型噪聲的協(xié)方差矩陣Q，測量噪聲的協(xié)方差矩陣R的影響，與觀測點(diǎn)數(shù)對(duì)算法跟蹤參數(shù)性能的影響。

圖1 三自由度系統(tǒng)Fig.1 A three degree of freedom system

在初始狀態(tài)下，k1=k4= 200 N / m，k2=k3= 100 N / m；m1=m2=m3= 1 kg。并假設(shè)這個(gè)三自由度系統(tǒng)的阻尼是比例阻尼，且有C=αM+βK，α=0.05，β=0.02。在1.5 s 時(shí)該結(jié)構(gòu)中的剛度k1從200 N / m 突然下降到150 N / m。假設(shè)激勵(lì)作用的位置和形式已知，質(zhì)量塊m1受到激勵(lì)f=sin(5π·t)+2·sin(2π·t)，由此可計(jì)算得到結(jié)構(gòu)各自由度的響應(yīng)。為了便于討論，選取加速度響應(yīng)作為測量向量，且不施加噪聲到測量向量中去。使用本文提出的擴(kuò)展卡爾曼濾波算法，此時(shí)激勵(lì)已知（作用位置，幅值大?。粗獏?shù)僅選擇k1，則擴(kuò)展?fàn)顟B(tài)向量z=[p p?k1]T是一個(gè)7×1列向量。

由于參數(shù)是未知的，所以初始參數(shù)的選擇必然會(huì)有較大誤差，而結(jié)構(gòu)的初始位移和速度（即狀態(tài)向量）卻可以事先測得，所以可以設(shè)定成真實(shí)值。相應(yīng)的初始估計(jì)誤差協(xié)方差P0|0關(guān)于狀態(tài)向量的部分可以設(shè)定為一個(gè)較小的值，初始估計(jì)誤差協(xié)方差P0|0關(guān)于參數(shù)的部分必須設(shè)定為一個(gè)較大的值。令初始擴(kuò)展?fàn)顟B(tài)向量z0，初始誤差方差P0|0分別寫為

3.1 模型噪聲協(xié)方差矩陣Q 的影響

討論模型噪聲協(xié)方差矩陣Q的取值對(duì)時(shí)變參數(shù)識(shí)別的影響。Q為一個(gè)7 階對(duì)角矩陣，前6 階對(duì)應(yīng)的是狀態(tài)更新方程對(duì)結(jié)構(gòu)的位移和速度所產(chǎn)生的誤差，第7 階代表的是狀態(tài)更新方程對(duì)參數(shù)k1所產(chǎn)生的誤差。兩者是不一樣的，因此必須分開討論。令測量噪聲的協(xié)方差矩陣R為對(duì)角元素為10-4的三階對(duì)角矩陣。首先考慮模型噪聲的協(xié)方差矩陣Q參數(shù)部分取值的影響，假定Q關(guān)于狀態(tài)向量部分的取值為10-10，關(guān)于參數(shù)部分的取值為x，即

x分別取值為1，10-1，10-2，10-3時(shí)，討論其參數(shù)部分的取值對(duì)參數(shù)識(shí)別效果的影響。同樣，為了討論狀態(tài)向量部分取值的影響，令Q關(guān)于參數(shù)部分的取值為10-1，關(guān)于狀態(tài)向量部分的取值為x。EKF 算法識(shí)別參數(shù)的效果如圖2 所示。

從圖2 可以看出，模型誤差Q關(guān)于參數(shù)部分傾向于取一個(gè)偏大的值，模型誤差Q關(guān)于狀態(tài)向量的部分傾向于取一個(gè)偏小的值。它本質(zhì)上是使Kk保持一個(gè)較大的值，算法對(duì)參數(shù)的追蹤性能也就會(huì)比較好；但Kk過大，會(huì)使得識(shí)別參數(shù)出現(xiàn)較大的波動(dòng)。在本算例中，狀態(tài)向量部分取值10-10是合適的；參數(shù)部分取值10-1是合適的。這樣的參數(shù)設(shè)置能很好地追蹤到參數(shù)的變化情況，且參數(shù)識(shí)別圖也不會(huì)有太大的波動(dòng)。

圖2 模型噪聲的協(xié)方差矩陣Q 對(duì)參數(shù)識(shí)別效果的影響Fig.2 Influence of covariance matrix Q of model noise on parameter identification effect

3.2 測量噪聲協(xié)方差矩陣R 的影響

對(duì)于測量噪聲的協(xié)方差矩陣R，其最優(yōu)取值與測量向量中的噪聲有關(guān)?，F(xiàn)討論對(duì)參數(shù)識(shí)別效果的影響。令模型噪聲的協(xié)方差矩陣Q狀態(tài)向量部分取值10-10，參數(shù)部分取值10-1。令測量噪聲的協(xié)方差矩陣R為對(duì)角元素為x的三階對(duì)角矩陣，當(dāng)x取不同值時(shí)，參數(shù)的識(shí)別結(jié)果如圖3 所示。

圖3 測量噪聲的協(xié)方差矩陣R 對(duì)參數(shù)識(shí)別效果的影響Fig.3 Influence of covariance matrix R of measurement noise on parameter identification effect

測量噪聲協(xié)方差矩陣R取一個(gè)較小值時(shí)，此時(shí)卡爾曼濾波增益保持一個(gè)較大的值，算法的追蹤參數(shù)性能較好。但值得注意的是，R的取值不能太小，它需要與觀測向量中的噪聲相適應(yīng)，如果選擇過小的值，導(dǎo)致這些估計(jì)的誤差超過了實(shí)際誤差，識(shí)別的參數(shù)必然不會(huì)與實(shí)際情況相同，嚴(yán)重時(shí)會(huì)發(fā)生發(fā)散。因此，需要為測量噪聲協(xié)方差矩陣R取一個(gè)合適的值，確保算法能較快地識(shí)別出結(jié)構(gòu)參數(shù)值，識(shí)別結(jié)果不會(huì)發(fā)散。

3.3 噪聲影響

當(dāng)響應(yīng)中存在噪聲干擾時(shí)，本文提出的算法仍能正確識(shí)別出參數(shù)變化情況。在此分別討論測量響應(yīng)中不含噪聲與含有5%噪聲的情況。令模型噪聲的協(xié)方差矩陣Q狀態(tài)向量部分取值10-10，參數(shù)部分取值10-1。由前文分析可知，測量噪聲的協(xié)方差矩陣R的最優(yōu)取值與噪聲水平的大小有關(guān)。對(duì)于響應(yīng)中不含噪聲時(shí)，設(shè)定測量噪聲的協(xié)方差矩陣R是對(duì)角元素為10-4的三階對(duì)角矩陣，可以看出參數(shù)識(shí)別效果很好；而對(duì)于響應(yīng)含有5%噪聲時(shí)，此時(shí)R的取值必須大于10-4，與觀測噪聲水平相適應(yīng)。若選取R是對(duì)角元素為10-2或10-3的三階對(duì)角矩陣，可以看出算法能較好地識(shí)別出參數(shù)變化情況。

由圖4 可以看出，無論測量向量中是否包含噪聲，算法實(shí)踐對(duì)時(shí)變參數(shù)的追蹤是通過調(diào)整測量噪聲協(xié)方差矩陣R的方式來實(shí)現(xiàn)的。

圖4 測量噪聲對(duì)參數(shù)識(shí)別效果的影響Fig.4 Effect of measurement noise on parameter identification

3.4 觀測點(diǎn)數(shù)影響

為說明觀測點(diǎn)數(shù)的影響，取一個(gè)懸臂梁作為仿真對(duì)象，梁長度L=1 m，矩形橫截面尺寸為寬b=0.1 m，高h(yuǎn)=0.01 m，材料密度ρ= 2 770 kg/m3，彈性模量E=71e+9 Pa，設(shè)懸臂梁的阻尼特性為瑞利阻尼，且有α=1，β=8e-7。再由公式I=b·h3/12，A=b·h可以推得該懸臂梁的截面慣性矩I與面積A。該懸臂梁的示意如圖5 所示。

圖5 懸臂梁結(jié)構(gòu)Fig.5 Cantilever beam structure

由圖5 可知，該懸臂梁被劃分為4 個(gè)單元，它的總體剛度矩陣和總體質(zhì)量矩陣是由各個(gè)單元的剛度矩陣Ki與單元質(zhì)量矩陣Mi組合而成的。因此該懸臂梁既包括轉(zhuǎn)角自由度，又包括豎向位移自由度，即4 個(gè)轉(zhuǎn)角自由度和4 個(gè)豎向位移自由度。其中Ki，Mi可寫為

為了便于討論，把梁單元的剛度特性Ei Ii作為單元發(fā)生剛度變化的指標(biāo)。假定梁單元①的剛度特性E1I1在第1.5～3.5 s 發(fā)生線性下降，下降至原來值的70%。其他單元的參數(shù)并不發(fā)生變化。在響應(yīng)中添加1% 的白噪聲，利用本文所提出的算法，把該懸臂梁前3 個(gè)單元的剛度特性作為未知參數(shù)向量，即擴(kuò)展?fàn)顟B(tài)向量為z=[p，p?，E1I1，E2I2，E3I3]T，依次識(shí)別出各個(gè)剛度系數(shù)。

選擇只觀測第2，3，4，5 節(jié)點(diǎn)上的豎直加速度，不測量轉(zhuǎn)角自由度，因此這也算是部分觀測。設(shè)模型噪聲的協(xié)方差矩陣G，以及測量噪聲的協(xié)方差矩陣R分別為

觀測全部測點(diǎn)與部分測點(diǎn)的區(qū)別僅需改變測量方程中的位置矩陣H0。若選擇只記錄結(jié)構(gòu)的豎向加速度響應(yīng)，因此可寫為

同理可求出既觀測轉(zhuǎn)角自由度，又觀測豎向自由度情況下的參數(shù)識(shí)別效果，由此可得參數(shù)識(shí)別結(jié)果如圖6～8 所示。

圖6 剛度系數(shù)E1I1識(shí)別效果圖Fig.6 Identification effect of stiffness coefficient E1I1

圖7 剛度系數(shù)E2I2識(shí)別效果圖Fig.7 Identification effect of stiffness coefficient E2I2

圖8 剛度系數(shù)E3I3識(shí)別效果圖Fig.8 Identification effect of stiffness coefficient E3I3

由圖6～8 可以看出，在含有噪聲的情況下，僅觀測部分測點(diǎn)也能追蹤到結(jié)構(gòu)參數(shù)的變化情況，且效果比較好。這說明了擴(kuò)展卡爾曼濾波算法能在測點(diǎn)部分觀測的情況下識(shí)別參數(shù)變化。

4 結(jié)論

在載荷已知的情況下，擴(kuò)展卡爾曼濾波算法對(duì)結(jié)構(gòu)參數(shù)的識(shí)別主要是通過測量更新方程實(shí)現(xiàn)的，而卡爾曼濾波增益Kk表示了算法對(duì)測量更新方程的利用程度。為了實(shí)現(xiàn)算法對(duì)時(shí)變參數(shù)的追蹤，本文提出了通過設(shè)置合適的濾波參數(shù)來保持較大的Kk。這其中包括：模型誤差Q關(guān)于參數(shù)部分取較大值，模型誤差Q關(guān)于狀態(tài)向量的部分取較小的值。還討論了測量噪聲對(duì)參數(shù)識(shí)別效果的影響，測量噪聲的協(xié)方差矩陣R必須與測量噪聲相適應(yīng)，并盡量選取一個(gè)較小值。仿真結(jié)果表明，這樣的濾波參數(shù)設(shè)置使得EKF 算法在運(yùn)行過程中能保持對(duì)測量更新方程的利用，從而能追蹤到結(jié)構(gòu)參數(shù)的變化情況，識(shí)別結(jié)構(gòu)的時(shí)變參數(shù)。最后，通過懸臂梁仿真算例，證明了在部分測點(diǎn)觀測的情況下，算法也能較好地識(shí)別出結(jié)構(gòu)的參數(shù)變化情況，體現(xiàn)了EKF 算法的優(yōu)越性。