基于改進粒子濾波算法的Bouc-Wen模型參數在線識別方法

王 濤，劉 媛，潘 毅，孟麗巖，包韻雷

(1.黑龍江科技大學 建筑工程學院，哈爾濱 150022；2. 中國地震局a.工程力學研究所;b. 地震工程與工程振動重點實驗室， 哈爾濱 150080；3. 西南交通大學 土木工程學院， 成都 610031；4. 抗震工程技術四川省重點實驗室， 成都 610031)

土木工程結構及其耗能減震控制元件在大震及超大震作用下表現出強非線性特征[1-2]。受到試驗設備加載能力及試驗費用的制約，傳統結構抗震試驗方法很難準確檢驗大型、復雜結構動力災變過程。結構混合試驗將物理試驗與數值模擬相結合，從而得到整體結構動力響應。目前，結構混合試驗受到了研究者們的廣泛關注[3-4]。受到試驗室作動器數量的限制，數值子結構將不可避免會進入非線性階段，數值模型誤差將會降低混合試驗精度。Bouc-Wen模型采用微分形式的數學表達，非線性模擬精度高和應用普適性強，可以用于模擬結構混合試驗中的數值子結構[5]。然而，Bouc-Wen模型參數本身物理意義并不明確，模型參數值對模型輸出精度具有較大影響，在混合試驗中如何準確給出Bouc-Wen模型參數值是亟需解決的問題。研究者們提出了在線模型更新方法，即利用試驗觀測數據在線識別試驗子結構模型參數，并實時更新數值子結構中Bouc-Wen模型參數[6-7]。

以卡爾曼濾波器為框架的一系列時域在線識別算法，如無跡卡爾曼濾波(unscented Kalman filter, UKF)[8-9]、約束無跡卡爾曼濾波(constrained unscented Kalman filter, CUKF)[10-11]被應用于結構模型更新混合試驗，以實現在線數值模型更新。然而，受到試驗和監(jiān)測過程中可能存在的非高斯噪聲和Bouc-Wen模型自身強非線性的影響，上述參數識別算法的精度和穩(wěn)定性會明顯降低，從而影響整體混合試驗結果精度。因此，找到一種精度更高的非線性參數在線識別算法仍是亟需解決的問題。

粒子濾波(particle filter, PF)算法[12]是一種基于貝葉斯估計和蒙特卡洛方法的在線非線性識別算法，其本質是通過尋找一組在狀態(tài)空間內傳播的隨機樣本來近似狀態(tài)概率密度函數，采用離散樣本模擬連續(xù)函數，以樣本均值代替積分運算，從而獲得狀態(tài)最小方差分布的過程。PF算法從理論上具有比擴展卡爾曼濾波器 (extended Kalman filter, EKF)和UKF算法更高的識別精度。Merwe等[13]對傳統的PF算法進行改進，在重要性采樣過程中利用了最新的觀測信息，更精確地逼近了后驗概率密度函數。Gordon等[14-16]提出了重采樣算法，解決粒子退化問題。李燁等[17]和唐和生等[18]采用PF算法解決結構系統的損傷識別問題，與EKF算法相比, 粒子濾波在非高斯噪聲條件下具有更高的結構模型參數識別精度。樊學平等[19]采用混合高斯粒子濾波器(mixed gaussian particle filter, MGPF)，對監(jiān)測信息狀態(tài)變量的后驗分布參數和監(jiān)測值一步向前預測分布參數進行了預測分析。目前，粒子濾波器在土木工程中的研究及應用仍非常有限，如何進一步提高PF算法的重要性采樣精度并同時削弱粒子匱乏仍是提高算法精度的關鍵問題。

文中首先在標準的PF算法基礎上擬建立一種改進的輔助無跡粒子濾波算法(auxiliary unscented particle filter, AUPF)，給出算法實現步驟。然后，針對單自由度Bouc-Wen模型進行參數在線識別，并與傳統的PF算法、EPF算法、UPF算法識別結果進行對比分析，驗證改進算法的精度和計算效率。最后，通過隔震支座擬靜力試驗驗證了采用AUPF算法在線識別Bouc-Wen模型參數方法的有效性。

1 輔助無跡粒子濾波算法的建立思路

文中提出的AUPF算法在標準粒子濾波算法的基礎上主要進行了兩方面改進：1)采用UKF算法進行重要性采樣，提高非線性系統粒子估計更新精度；2)在重采樣過程中引入輔助因子，修改粒子權值，通過增加粒子多樣性來削弱粒子退化現象。AUPF算法繼承了PF算法原理，理論上可以應用于任意非線性模型參數識別。AUPF算法具體流程如圖1所示，算法主要包括重要性采樣、權值計算和輔助重采樣3個主要環(huán)節(jié)。

圖1 AUPF算法示意圖

2 輔助無跡粒子濾波算法的實現步驟

對于任意一個非線性、非高斯動態(tài)系統，系統狀態(tài)方程、觀測方程分別定義為式(1)和式(2)：

xk=f(xk-1,uk-1)+wk-1,

(1)

yk=h(xk)+vk,

(2)

式中：k為遞推步數；xk為系統的狀態(tài)值，假設xk具有一階馬爾可夫性質。函數f(·)為系統的狀態(tài)方程，uk為系統的輸入；wk-1為系統的過程噪聲；yk為系統的觀測值；h(·)為系統的觀測方程；vk為系統的觀測噪聲；wk-1和vk為2組相互獨立、互不相關的噪聲序列，假定已知其概率密度函數；Qk、Rk分別為過程噪聲協方差矩陣和觀測噪聲協方差矩陣。

2.1 重要性函數的選取

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

新息協方差Py ik|k-1為

(11)

交叉協方差Pxy ik|k-1為

(12)

(13)

利用最新的觀測yk，計算第k步第i個粒子估計的均值和協方差為

(14)

(15)

(16)

在整個過程中，AUPF算法通過UKF算法對非線性模型進行直接處理，得到算法的重要性函數，避免了繁瑣的雅克比矩陣的求解，降低了計算復雜度，同時使得AUPF算法的重要性函數中包含最新的系統觀測信息。

2.2 重要性權值的計算

當由公式(14)計算得到第k步粒子的估計值后，需要通過調整每一個粒子重要性權值，并將每個粒子權值進行歸一化，以更好地逼近狀態(tài)后驗概率密度函數[21]。歸一化后的重要性權值為

(17)

(18)

(19)

2.3 輔助重采樣

(20)

在重采樣過程中，為增加粒子多樣性、減小其中權值較小粒子有效信息的喪失，文中在標準PF算法基礎上引入的輔助因子λ，重新計算粒子的觀測似然函數為

(21)

圖2 AUPF算法流程

3 基于輔助無跡粒子濾波算法的Bouc-Wen模型參數識別驗證

AUPF算法繼承了PF算法原理，是一種完全的非線性估計器，可以識別任意非線性模型參數。Bouc-Wen模型是一種具有代表性的強非線性模型，被廣泛用來模擬結構和構件恢復力特性。文中以一單自由度Bouc-Wen模型為對象，給出應用AUPF算法進行非線性模型參數在線識別的具體實現方法，檢驗算法識別精度。結合結構運動方程，Bouc-Wen模型如式(22)～式(24)所示：

(22)

(23)

F=k0z，

(24)

設Bouc-Wen模型參數真實取值為k0=40 kN/m，β=20，γ=20，n=1.1；對模型進行位移控制加載，輸入位移激勵選用1940年5月19日Imperial Valley地震El Centro(1940，NS)臺站測得的地面運動位移記錄，位移峰值調整為10 cm，如圖3所示，圖中縱坐標d為位移。采用4階Runge-Kutta數值積分方法計算Bouc-Wen系統恢復力，積分步長Δt=0.01 s，積分步數為4 000步。

圖3 位移加載時程曲線

在結構加載過程中，基于當前及之前步結構真實反應和系統輸入，采用AUPF算法在線識別Bouc-Wen模型參數，并與PF算法、EPF算法和UPF算法的識別結果進行對比，用以驗證AUPF算法參數識別精度。為了能夠體現出改進算法在重要性采樣及重采樣的影響，以上4種算法均采用相同的狀態(tài)及參數設置。設系統狀態(tài)為x=[x1,x2,x3,x4,x5]T=[z,k0,β,γ,n]T；系統的觀測為yk=Rk，其中，下標k表示步數，Rk為第k步的結構恢復力；假定系統過程噪聲vk和觀測噪聲wk都服從高斯分布，即均值均為0，協方差分別為Q、W。AUPF算法中輔助因子λ取1.1。

假定狀態(tài)估計初始值為x0=[0， 50， 15， 15， 2]T，過程噪聲協方差為Q=diag(10-8，0.0452，0.000 092，0.000 012，0.0072)，觀測噪聲協方差為W=0.015 kN2；狀態(tài)估計誤差初始協方差為P0=diag(10-6，113.9，15.6，12.7，0.65)。

系統狀態(tài)方程為

(25)

(26)

設系統觀測方程為

yk=Fk=x1,kx2,k+wk。

(27)

將基于AUPF算法、UPF算法、EPF算法和PF算法得到的Bouc-Wen模型在線參數識別結果進行對比，如圖4所示。可以看出，對于k0和n2個參數的識別效果4種算法大致相同，基本都收斂到了真實值附近，且收斂速度大致相似，其中AUPF算法得到的識別值同真實值吻合度最好。PF算法、EPF算法、UPF算法、AUPF算法得到的參數β識別終值相對誤差分別為22.34%、14.95%、3.33%、2.46%，參數γ識別終值相對誤差分別為20.41%、17.25%、14.03%、5.58%。與其他3種算法結果相比， AUPF算法提高了模型參數的識別精度，同時顯著減小了參數識別值收斂過程波動幅度。由于Bouc-Wen模型參數β和γ為控制滯回環(huán)形狀參數，本身無物理意義。在相同的模型輸入下，不同參數β和γ的組合可以得到相同的模型輸出，因此，導致反問題的參數識別值可能并不唯一，識別誤差增大。AUPF算法在重要性采樣具有更高的非線性變換精度，在重采樣過程中豐富了粒子多樣性，有效削弱粒子退化。算法性能決定了算例識別結果的優(yōu)劣，具有普遍意義。因此，當算例中的Bouc-Wen模型參數的真實值取值發(fā)生變化，在相同的條件下，4種算法識別結果仍會有類似的規(guī)律，由于篇幅有限，沒有給出采取其他參數真實值時的識別結果。

圖4 Bouc-Wen模型參數識別值

為了能定量評價算法識別精度，定義一次獨立仿真的均方根誤差(root-mean-square error，RMSE)為

(28)

PF算法及其改進的EPF算法、UPF算法和AUPF算法本質上均為隨機性參數識別算法，4種算法均基于蒙特卡洛隨機采樣方法，因此，即使在相同的參數初值條件下，每一種算法在每一次仿真得到的參數識別值都是不同的，即參數識別結果具有隨機性。為了檢驗算法識別結果的隨機性，采用4種濾波算法分別進行了10次獨立仿真，統計識別結果來對比分析不同算法的識別精度和收斂性，更具有說服性。在本算例中的10次獨立仿真中，系統輸入、Bouc-Wen模型初始參數和算法初始參數均相同，隨機性主要來自算法生產粒子的隨機性。4種算法參數識別值的均方根誤差與仿真次數的關系曲線如圖5所示，圖中的橫坐標為仿真次數。

圖5 模型參數識別值均方根誤差

可以看出，AUPF算法的參數識別整體誤差明顯低于PF算法、EPF算法和UPF算法，而且誤差波動幅度顯著降低。可見，由于AUPF算法利用最新觀測信息修正粒子，同時通過引入輔助因子增加了粒子多樣性。因此，AUPF算法明顯高于PF算法、EPF算法和UPF算法的識別精度。

統計10次獨立仿真在線參數識別值的均方根誤差RMSE均值、相對誤差RE(relative error)均值如圖6所示。可以明顯看出，文中提出的AUPF算法的參數識別值均方根誤差均值和相對誤差均值整體上都要小于PF算法、EPF算法和UPF算法參數識別誤差。AUPF算法得到的4組參數識別值的均方根誤差整體上比PF算法、EPF算法和UPF算法結果誤差減小了81.5%、37.7%和8.0%，AUPF算法得到的4組參數識別值的相對誤差整體上比PF算法、EPF算法和UPF算法結果誤差減小了87.3%、39.0%和61.8%。由此可見，10次仿真4組參數識別值的平均均方根誤差和平均相對誤差2種評價指標均表明AUPF算法精度均高于其它3種算法。需要注意的是，算法在取得較高參數識別精度同時，也需要付出更多的計算耗時。10次仿真4種算法的單步平均計算耗時如圖7所示。其中，AUPF算法與UPF算法的單步平均計算耗時幾乎相同，均為0.20 s，相當于標準粒子濾波算法耗時的2倍和 EPF算法的6倍左右。其原因主要是由于AUPF算法和UPF算法在重要性采樣時，對于每一個粒子均值均需要采用UKF方法計算，這將顯著增加算法的計算耗時。

圖6 參數識別值均方根誤差及相對誤差均值

圖7 單步平均計算耗時

4 橡膠隔震支座參數識別

為了驗證AUPF算法對于真實物理試驗中進行參數識別的有效性，采用直徑為300 mm的鉛芯橡膠隔震支座擬靜力試驗所測得的水平剪力和位移數據，進行了在線模型參數識別。參數識別時，假定隔震支座恢復力模型為Bouc-Wen模型。擬靜力試驗模型為LRB300鉛芯橡膠隔震支座，其質量分別為81 kg和82 kg。鉛芯橡膠隔震支座的設計承載力為566 kN，橡膠直徑為300 mm，橡膠總厚度為48 mm，支座高度為100 mm，一次形狀系數為9.375，二次形狀系數為6.250，水平等效剛度為1.017 kN/mm，豎向剛度為608 kN/mm。隔震支座實物如圖8所示。試驗在哈爾濱工業(yè)大學結構與力學實驗中心完成，試驗加載設備是華龍20 MN 動態(tài)壓剪試驗機，豎向最大加載壓力為2 000 t，行程為±100 mm；水平最大加載壓力為200 t，行程為±500 mm。試驗機采樣頻率為0.01 Hz。隔震支座在工作臺上由2塊擋板固定下連接板，當工作臺移動到試驗位置附近時，調節(jié)工作臺使之居中，并固定上連接板。

圖8 LRB300隔震支座實物圖

橡膠隔震支座Bouc-Wen模型參數識別時程曲線如圖9所示。參數的識別終值分別為：k0=0.382 8 kN/mm、β=-0.007 73、γ=0.009 6和n=1.374 37。將AUPF算法識別得到的支座滯回曲線與試驗真實滯回曲線進行對比，如圖10所示。

圖9 Bouc-Wen模型參數識別值

圖10 隔震支座滯回曲線

可以看出，當基于隔震支座真實試驗數據采用AUPF算法進行在線參數識別時，除了Bouc-Wen模型參數n收斂較慢，模型參數k0、β和γ均可以很快地收斂于穩(wěn)定值。參數n識別值收斂較慢主要有兩方面的原因：一是存在模型誤差；二是模型具有強非線性。AUPF算法是一種基于模型的參數識別算法，當假定模型與真實系統之間存在模型誤差，就會降低算法識別精度。由圖4(b)可見，當識別算法中假定的系統模型和真實系統模型均為相同的Bouc-Wen模型時，即算法不存在模型誤差時，模型參數n的識別值收斂較快，基本可以收斂到參數的真實值，可見算法對參數n具有較好的識別效果。由圖9(d)可見，模型參數n的識別值具有一定的波動性，并沒有很快地收斂于穩(wěn)定值，其主要原因是由于識別算法采用Bouc-Wen模型來近似真實的隔震器的滯變關系仍存在一定模型誤差，模型誤差會降低識別算法的識別精度。另外，模型中的時變滯變位移z為參數n的指數函數，具有較高程度的非線性。以上原因導致參數n識別值具有時變性，以補償模型誤差的不利影響。

由于Bouc-Wen模型參數沒有明確的物理意義，很難確定隔震支座所對應的模型參數真實值，因此并不能直接評價參數識別值的精度。為了驗證識別參數準確性，將識別參數在線計算得到的水平恢復力和試驗測量得到的恢復力進行對比，可以看出，滯回曲線識別值和實驗值吻合較好，表明采用AUPF算法在線識別Bouc-Wen模型參數具有較高的識別精度，同時也表明Bouc-Wen模型可以很好地模擬鉛芯橡膠隔震支座力學行為。

5 結 論

采用AUPF算法對一種強非線性Bouc-Wen模型進行了在線參數識別，分析了AUPF算法的精度和計算效率，得到以下結論。

1)基于標準PF算法提出了一種改進的AUPF算法，采用UKF算法進行粒子重要性采樣提高粒子非線性變化精度，同時引入輔助因子提高重采樣粒子多樣性，削弱粒子退化。

2)與標準PF算法、EPF算法以及UPF算法相比，提出的AUPF算法得到的Bouc-Wen模型參數識別值的均方根誤差整體上減小了81.5%、37.7%和8.0%，并能有效降低模型參數識別過程中的波動幅度，提高了強非線性模型在線參數識別精度。

3)AUPF算法平均單步計算耗時為0.20 s，與PF算法、EPF算法相比，計算耗時有所增加。需考慮減少UKF算法重要性采樣計算負荷，提高算法計算效率。

4)采用鉛芯橡膠支座擬靜力試驗數據，驗證了AUPF算法在線識別Bouc-Wen模型參數方法的有效性。