樊學(xué)平， 屈 廣， 劉月飛

(1.蘭州大學(xué) 土木工程與力學(xué)學(xué)院， 甘肅 蘭州 7300001； 2.蘭州大學(xué)西部災(zāi)害與環(huán)境力學(xué)教育部重點實驗室， 甘肅 蘭州 730000)

目前，結(jié)構(gòu)健康監(jiān)測(structural health monitoring, SHM)領(lǐng)域的研究主要集中在利用傳感器采集數(shù)據(jù)和監(jiān)測數(shù)據(jù)的應(yīng)用兩個方面.前者主要集中在數(shù)據(jù)壓縮、數(shù)據(jù)恢復(fù)、數(shù)據(jù)獲得技術(shù)和系統(tǒng)組裝技術(shù)等方面[1-5]，目前在硬件和技術(shù)上已經(jīng)趨于成熟，對于后者而言，國內(nèi)外學(xué)者的研究主要集中在結(jié)構(gòu)損傷識別、結(jié)構(gòu)模態(tài)參數(shù)識別以及結(jié)構(gòu)模型修正等方面[6-8].基于SHM數(shù)據(jù)的橋梁可靠性預(yù)測及評定，國內(nèi)外學(xué)者多是基于離線監(jiān)測信息展開研究[9-13]，而基于實時監(jiān)測信息的結(jié)構(gòu)性能動態(tài)預(yù)估相對較少，且離線的結(jié)構(gòu)性能預(yù)測及評定模型在實際工程應(yīng)用中存在一定局限性，如：基于一次回歸函數(shù)與常值函數(shù)的橋梁可靠性預(yù)測分析[9-10]，均未考慮監(jiān)測變量的隨機性和非平穩(wěn)性；基于ARMA(autoregressive moving average model)模型的橋梁構(gòu)件荷載效應(yīng)預(yù)測[11]，以及基于單一或組合貝葉斯動態(tài)模型的橋梁可靠性動態(tài)預(yù)測[12-13]，多是基于平穩(wěn)監(jiān)測信息進行可靠性分析，而考慮到橋梁動力響應(yīng)的復(fù)雜性，基于非線性及非平穩(wěn)監(jiān)測信息的結(jié)構(gòu)性能預(yù)測方法還需進一步研究.近年來，粒子濾波器因在非線性問題中出眾的性能而倍受關(guān)注[14].在SHM領(lǐng)域中，粒子濾波器在結(jié)構(gòu)應(yīng)力預(yù)測、結(jié)構(gòu)可靠性動態(tài)預(yù)測、結(jié)構(gòu)系統(tǒng)識別以及結(jié)構(gòu)損傷識別等方面已取得一些研究成果，如：基于高斯混合粒子濾波器、改進高斯混合粒子濾波器以及折扣高斯粒子濾波器的橋梁可靠性預(yù)測[12,15-17].但粒子濾波器假定重要性采樣是指能夠從一個合理的后驗建議密度函數(shù)中得到一組覆蓋真實狀態(tài)的樣本集合.因此，找到一個最優(yōu)建議密度函數(shù)來指導(dǎo)采樣過程，就能大大提高粒子的利用率，進而可以有效降低粒子退化的影響，提高系統(tǒng)的識別精度.橋梁長期運營中，由于其構(gòu)件受到多種荷載耦合作用，SHM系統(tǒng)所獲取的結(jié)構(gòu)動力響應(yīng)數(shù)據(jù)十分復(fù)雜，并且受諸如氣候條件、橋下水位等外在環(huán)境因素影響較大，因此橋梁應(yīng)力(或撓度)狀況時常會出現(xiàn)變化，這就要求所選取的粒子濾波建議分布具有長期穩(wěn)定性、受外界干擾小、精度高等特點，并且具有能隨橋梁狀況變化而不斷更新的良好自適應(yīng)性.因而，如何選取具有以上特性的建議分布就成為利用粒子濾波方法預(yù)測橋梁時變可靠性的難點所在.

為合理應(yīng)用監(jiān)測信息對橋梁動態(tài)可靠度指標進行實時在線預(yù)測分析，基于橋梁監(jiān)測極值應(yīng)力(或撓度)信息，利用貝葉斯動態(tài)線性模型和粒子濾波器相融合的改進粒子濾波預(yù)測算法，結(jié)合FOSM(first order second moment)方法對橋梁可靠度指標進行了動態(tài)預(yù)測，并進行了實例和設(shè)計試驗驗證分析.改進粒子濾波方法利用貝葉斯動態(tài)線性預(yù)測模型為動態(tài)粒子提供建議分布，并充分利用動態(tài)健康監(jiān)測數(shù)據(jù)對建議分布進行遞推更新，有效提高了粒子的利用率，并且引入貝葉斯動態(tài)模型折扣因子[12]，使得粒子濾波在遞推過程中具有較好的自適應(yīng)性.

1 監(jiān)測數(shù)據(jù)的動態(tài)模型

SHM系統(tǒng)在長期運營過程中積累了大量監(jiān)測應(yīng)力(或撓度)信息，定義每小時(或每天)監(jiān)測應(yīng)力(或撓度)的極大值為監(jiān)測極值應(yīng)力(或撓度)信息.定義{Xt,t=1,2,…}為極值應(yīng)力(或撓度)狀態(tài)時間序列，Xt表示t時刻的極值應(yīng)力(或撓度)狀態(tài).所建動態(tài)方程基于以下兩點假設(shè)：

(1) {Xt,t=1,2,…}具有馬爾科夫性.

(2) 監(jiān)測信息{Zt,t=1,2,…}之間相互獨立，且Zt僅和狀態(tài)變量Xt相關(guān)，Zt與Xt呈線性關(guān)系.

橋梁極值應(yīng)力(或撓度)的動態(tài)模型如下：

狀態(tài)方程為

Xt=f(Xt-1)+wt,wt～N(0,Wt)

(1)

監(jiān)測方程為

Zt=Xt+vt,vt～N(0,Vt)

(2)

初始狀態(tài)信息為

p(Xt-1|Dt-1)=N(mt-1,Ct-1)

(3)

式(1)～(3)中：Xt為狀態(tài)變量；f(·)為狀態(tài)轉(zhuǎn)移函數(shù)，可以通過極值應(yīng)力(或撓度)監(jiān)測數(shù)據(jù)利用五點三次平滑處理之后的時間序列近似回歸擬合得到[12]；wt為狀態(tài)白噪聲;Wt為狀態(tài)誤差的方差；Zt為t時刻的極值應(yīng)力(撓度)監(jiān)測值；vt為監(jiān)測白噪聲；Vt為監(jiān)測誤差的方差；初始狀態(tài)變量的概率密度函數(shù)(probability density function, PDF)p(Xt-1|Dt-1)可根據(jù)既有極值應(yīng)力(或撓度)信息平滑得到的時間序列樣本統(tǒng)計得到；N(·)為正態(tài)PDF.

基于先驗信息集Dt和當前監(jiān)測數(shù)據(jù)Zt+1估計狀態(tài)Xt+1的后驗PDF,即p(Xt+1|Dt+1)，求解過程包含以下兩個階段：

(1) 預(yù)測.已知t時刻狀態(tài)的后驗PDFp(Xt|Dt)，利用式(1)可得t+1時刻狀態(tài)Xt+1的先驗PDF為

(4)

(2) 測量更新.在t+1時刻，利用當前監(jiān)測極值Zt+1，由Bayes公式更新式(4)，可得后驗PDF為

(5)

其中

p(Ztt+1|Xt+1)dXt+1

(6)

為極值的一步預(yù)測PDF.

2 改進的粒子濾波算法

結(jié)合重要性采樣與重采樣方法，改進的粒子濾波器(IPF)通過蒙特卡羅方法以參數(shù)化的分布逼近式(4)和式(5)中的預(yù)測PDF與后驗狀態(tài)PDF，是式(4)和(5)近似求解的一種策略.

2.1 基于BDLM(Bayesian dynamic linear model)的近似重要性采樣分布

粒子濾波的有效性依賴于重要性采樣分布的合理性.較為合理的重要性采樣分布能夠有效減少具有偏差的樣本，提升粒子濾波性能.因此，引入BDLM[12]在重要性采樣階段近似產(chǎn)生重要性采樣PDF.假設(shè)采樣過程中抽取了n個粒子，則更新粒子集合的步驟如下：

(1)t時刻每個粒子的狀態(tài)后驗PDF為

(7)

(2)t+1時刻的狀態(tài)先驗PDF為

(8)

(3)t+1時刻的一步預(yù)測PDF為

(9)

(4)t+1時刻的狀態(tài)后驗PDF為

(10)

式(10) 用來近似模擬重要性采樣分布.

2.2 改進的粒子濾波算法

粒子濾波器能夠為監(jiān)測值與狀態(tài)值的后驗概率分布進行高精度表達[18]，在非線性系統(tǒng)中具有出色的表現(xiàn).然而，傳統(tǒng)粒子濾波器以給定的高斯分布作為建議分布，并沒有充分利用最新監(jiān)測信息對建議分布進行動態(tài)更新[19]，因此通常需要大量粒子樣本來近似接近系統(tǒng)的后驗PDF，并且在狀態(tài)更新時計算量較大，具有一定盲目性.特別是當橋梁的應(yīng)力(或撓度)環(huán)境發(fā)生變化時，由粒子濾波產(chǎn)生的樣本就會出現(xiàn)較大偏差.因而有效減少樣本數(shù)量的自適應(yīng)采樣策略是該算法的重點所在.因此，將自適應(yīng)性及魯棒性良好的BDLM應(yīng)用到粒子濾波中，為其提供重要性采樣，不僅可以獲取更為準確的建議分布，而且能夠提升粒子濾波算法的魯棒性與自適應(yīng)性，同時大大減少了粒子濾波算法的樣本需求量，降低了算法復(fù)雜度.

改進后的算法步驟如下：

(2) 由式(7)～(10)遞推得到每個粒子的建議分布，亦即在t時刻用BDLM以及最新監(jiān)測信息Zt+1更新粒子，則第i個粒子更新為

(11)

利用式(2)、(8)和(11)可得粒子權(quán)重為

(12)

將權(quán)重進行歸一化，可得

(13)

(3) 利用重采樣算法[20]，根據(jù)歸一化權(quán)值的大小對粒子集合進行復(fù)制和淘汰.

(14)

(6) 重復(fù)(2)～(5)步，可實現(xiàn)監(jiān)測變量的動態(tài)預(yù)測.

一步預(yù)測精度可由均方誤差(EMSE)來衡量，即

(15)

顯然，均方誤差越小，預(yù)測精度越高.

改進算法流程如圖1所示.

圖1 改進的粒子算法流程

3 可靠度指標預(yù)測計算

3.1 工程實例1預(yù)測公式

基于天津富民橋主纜的應(yīng)力監(jiān)測數(shù)據(jù)，結(jié)合FOSM方法對主纜可靠性進行動態(tài)預(yù)測.主纜鋼絲的抗拉強度按照[σ]=1 670 MPa進行計算，主纜截面失效模式對應(yīng)的功能函數(shù)為

g(σp，t)=[σ]-σp,t

(16)

式中:σp,t為t時刻預(yù)測的極值應(yīng)力，且隨機變量之間相互獨立.

主纜鋼絲抗拉強度的平均值為μσ=1 670 MPa，變異系數(shù)δσ=0.15，則利用FOSM方法可得主纜的動態(tài)監(jiān)測可靠度指標βt為

(17)

考慮到預(yù)測過程的不確定性與隨機性，則動態(tài)可靠度指標預(yù)測公式為

(18)

式中:μσp,t和Vσp,t分別為t時刻主纜極值應(yīng)力基于改進粒子濾波器的一步預(yù)測均值和方差.

本算例中的極值應(yīng)力為每小時監(jiān)測應(yīng)力的極大值.

3.2 工程實例2預(yù)測公式

參考已有文獻[12]，可得美國I-39北橋第二跨橫梁可靠度指標預(yù)測公式為

(19)

式中:μM和σM分別為基于改進粒子濾波器預(yù)測得到的極值應(yīng)力平均值和標準差；μR和σR分別為按照規(guī)范計算的抗力的平均值和標準差；μS和σS分別為由鋼板恒載所引起的應(yīng)力的平均值和標準差；μC和σC分別為由混凝土恒載所引起的應(yīng)力的平均值和標準差；γM=1.15是傳感器的修正系數(shù).

由文獻[12]可知：μR=380 MPa,σR=26.6 MPa，結(jié)合由結(jié)構(gòu)自重產(chǎn)生的應(yīng)力分布參數(shù),μS=116.3 MPa,σS=116.3 MPa×0.04=4.65 MPa,μC=108.8 MPa以及σC=108.8 MPa×0.04=4.35 MPa，可得第二跨橫梁的可靠度指標預(yù)測公式為

(20)

參考式(17)，可得第二跨橫梁的可靠度指標監(jiān)測值βt為

(21)

式中:Mt為t時刻的極值應(yīng)力監(jiān)測值.

本算例中的極值應(yīng)力為每天監(jiān)測應(yīng)力的極大值.

3.3 工程實例3預(yù)測公式

參考已有文獻[12]，可得伊通河橋上游主梁動態(tài)可靠度指標預(yù)測公式為

(22)

參考式(17)和(21)，可得主梁動態(tài)可靠度指標監(jiān)測值βt為

(23)

式中:f(t)為t時刻的極值撓度監(jiān)測值.

本算例中的極值撓度為每小時監(jiān)測撓度的極大值.

4 算例分析

4.1 工程實例1計算分析

天津富民橋的結(jié)構(gòu)及其示意圖分別如圖2和圖3所示.富民橋總長340.3 m，橋?qū)?0 m.橋梁動力響應(yīng)復(fù)雜，結(jié)構(gòu)受溫度荷載影響較大.基于此懸索橋主纜的監(jiān)測信息，對主纜可靠性進行動態(tài)預(yù)測.

圖2 天津富民橋

圖3 富民橋結(jié)構(gòu)示意圖

此主纜共安裝了2個傳感器，取其中歷史監(jiān)測極值應(yīng)力絕對值較大的FBG01192傳感器作為實時監(jiān)控對象.定義每小時監(jiān)測應(yīng)力的最大值為監(jiān)測極值應(yīng)力.從2009年8月1日到2009年9月3日，對此橋進行了810 h的健康監(jiān)測.其中前300 h的監(jiān)測極值應(yīng)力如圖4所示.對其采用五點三次平滑方法[9]處理后的數(shù)據(jù)如圖5所示，處理后的數(shù)據(jù)近似作為初始狀態(tài)信息，此初始狀態(tài)信息用正態(tài)PDF模擬.

圖4 極值應(yīng)力監(jiān)測值(工程實例1)

圖5 極值應(yīng)力初始值與監(jiān)測值(工程實例1)

采用前300 h的極值應(yīng)力信息，并結(jié)合式(1)～(3)，建立的動態(tài)模型如下：

狀態(tài)方程為

Xt=Xt-1-0.006 6+wt,wt～N(0,Wt)

(24)

監(jiān)測方程為

Zt=Xt+vt,vt～N(0,0.103 2)

(25)

初始狀態(tài)信息為

Xt-1|Dt-1～N(147.634 1,6.261 7)

(26)

利用改進的粒子濾波算法，結(jié)合式(7)～(15)與式(24)～(26)，基于第300小時到第809小時的監(jiān)測極值應(yīng)力數(shù)據(jù)，對第301小時到第810小時的極值應(yīng)力進行實時動態(tài)一步預(yù)測，結(jié)果如圖6所示.一步預(yù)測均方誤差可由式(15)計算，連續(xù)50次濾波結(jié)果的EMSE見圖7.可以看出,一步預(yù)測值與實測值幾乎重合，并且均方誤差較小(0.068左右波動)且平穩(wěn)可觀，驗證了模型對于監(jiān)測應(yīng)力預(yù)測的有效性.將預(yù)測結(jié)果與基本粒子濾波(PF)[21]及差分自回歸移動平均模型(ARIMA)[22]進行比較，預(yù)測結(jié)果對比如圖8所示.經(jīng)計算得到粒子濾波的均方誤差為0.146 8(平均50次濾波EMSE的均值)，ARIMA模型的均方誤差為0.475 8.由此可見改進粒子濾波器的預(yù)測精度十分可觀，進一步驗證了所提模型的準確性.

圖6 基于改進粒子濾波器的極值應(yīng)力預(yù)測值與實測值(工程實例1)

Fig.6Monitoring and predicted extreme stresses based on IPF

圖7 連續(xù)50次濾波的預(yù)測均方誤差(工程實例1)

圖8 3種算法極值應(yīng)力預(yù)測結(jié)果比較(工程實例1)

結(jié)合式(17)和式(18)對主纜可靠性進行實時預(yù)測分析，得到結(jié)果如圖9所示，可以看出可靠度指標預(yù)測值與實測值幾乎重合，很好地預(yù)測了可靠度指標的變化趨勢，而考慮了不確定性與隨機性的可靠度指標預(yù)測值略低于實測值，但由于考慮了不確定性，預(yù)測結(jié)果更符合工程實際.

圖9 可靠度指標的實測值與一步預(yù)測值(工程實例1)

4.2 工程實例2計算分析

采用美國I-39北橋第二跨橫梁83 d的極值應(yīng)力監(jiān)測數(shù)據(jù)對所提算法進行驗證.極值應(yīng)力數(shù)據(jù)如圖10所示.

圖10 極值應(yīng)力監(jiān)測值(工程實例2)

首先基于前50 d的監(jiān)測極值應(yīng)力數(shù)據(jù)，結(jié)合式(1)～(3)建立動態(tài)模型如下:

狀態(tài)方程為

Xt=Xt-1-0.056 3+wt,wt～N(0,Wt)

(27)

監(jiān)測方程為

Zt=Xt+vt,vt～N(0,8.814)

(28)

初始狀態(tài)信息為

Xt-1|Dt-1～N(24.335 1,7.138 7)

(29)

然后采用所提IPF預(yù)測算法，結(jié)合式(7)～(15)與式(27)～(29)，基于第50天到第82天的監(jiān)測極值應(yīng)力信息，對第51天到第83天的極值應(yīng)力進行動態(tài)預(yù)測，預(yù)測結(jié)果如圖11所示.并與PF[21]及 ARIMA[22]進行比較，預(yù)測結(jié)果如圖12所示.利用式(15)得到IPF、PF以及ARIMA(0,0,0)的EMSE分別為18.897 7、37.612 7、27.238 3，由此可見IPF準確性較好.

結(jié)合式(20)～(21)，可得可靠度指標預(yù)測值和監(jiān)測值如圖13所示.可以看出,可靠度指標預(yù)測值較好反映了監(jiān)測可靠度指標的變化趨勢，而考慮了不確定性與隨機性的可靠度指標預(yù)測值略低于實測值，但由于考慮了不確定性，預(yù)測值更符合工程實際.

圖11 基于改進粒子濾波器的極值應(yīng)力預(yù)測值與實測值(工程實例2)

Fig.11 Monitored and predicted extreme stresses

圖12 3種算法極值應(yīng)力預(yù)測結(jié)果比較(工程實例2)

圖13 可靠度指標的實測值與一步預(yù)測值(工程實例2)

4.3 工程實例3計算分析

采用吉林長春伊通河橋上游主梁70 h的極值撓度監(jiān)測數(shù)據(jù)對所提算法進行驗證.極值撓度數(shù)據(jù)如圖14所示.

首先基于前50 h的監(jiān)測極值撓度數(shù)據(jù)，結(jié)合式(1)～(3)建立動態(tài)模型如下:

狀態(tài)方程為

Xt=Xt-1-0.049 5+wt,wt～N(0,Wt)

(30)

監(jiān)測方程為

Zt=Xt+vt,vt～N(0,1.421 3)

(31)

圖14 極值撓度監(jiān)測值(工程實例3)

初始狀態(tài)信息為

Xt-1|Dt-1～N(45.869 2,0.953 9)

(32)

然后采用所提IPF預(yù)測算法，結(jié)合式(7)～(15)與式(30)～(32)，基于第50小時到第69小時的監(jiān)測極值撓度信息，對第51小時到第70小時的極值撓度進行動態(tài)預(yù)測，預(yù)測結(jié)果如圖15所示.并與PF[21]及 ARIMA[22]進行比較，預(yù)測結(jié)果對比如圖16所示.利用式(15)得到IPF、PF以及ARIMA(1,0,2)的EMSE分別為0.907 6、1.628 3、1.427 7，由此可見IPF預(yù)測性能較好.

圖15 基于改進粒子濾波器的極值撓度預(yù)測值與實測值(工程實例3)

Fig.15Monitored and predicted extreme deflections based on IPF

圖16 3種算法極值撓度預(yù)測結(jié)果比較(工程實例3)

結(jié)合式(22)和式(23)，可得可靠度指標分析結(jié)果如圖17所示.可以看出可靠度指標預(yù)測值很好地反映了監(jiān)測可靠度指標的變化趨勢，而考慮了不確定性與隨機性的可靠度指標預(yù)測值略低于實測值，但由于考慮了不確定性，預(yù)測結(jié)果更符合工程實際.

圖17 可靠度指標的實測值與一步預(yù)測值(工程實例3)

4.4 橋梁環(huán)境突變的假設(shè)試驗

由于橋梁的應(yīng)力環(huán)境受外界因素影響較大，當在役橋梁氣候條件、橋下水位等外在環(huán)境因素發(fā)生突變時，橋梁的應(yīng)力狀況會也會隨之出現(xiàn)變化.假設(shè)4.1中富民橋在監(jiān)測到第600小時候橋梁環(huán)境發(fā)生突變，例如受到洪水影響，現(xiàn)假定其實際監(jiān)測信息受到如下變化：

(33)

現(xiàn)利用改進的粒子濾波算法，結(jié)合式(7)～(15)、式(24)～(26)以及式(33)，基于第300小時到第809小時的監(jiān)測極值應(yīng)力數(shù)據(jù)，對第301小時到第810小時的監(jiān)測極值應(yīng)力進行實時動態(tài)一步預(yù)測，同時與PF方法結(jié)果進行對比，預(yù)測結(jié)果如圖18所示.可以明顯地看出，由于IPF的建議分布在更新過程中借助了突變后的應(yīng)力信息，因此仍能產(chǎn)生較為合理的建議分布，有效地對突變后的監(jiān)測信息進行預(yù)測，而PF算法由于粒子集合的分布不能很好地覆蓋真實值，因此當監(jiān)測信息突變時出現(xiàn)了難以避免的的濾波發(fā)散.

圖18 橋梁環(huán)境突變下的應(yīng)力預(yù)測結(jié)果

5 結(jié)論

提出了一種用于橋梁動態(tài)可靠度指標在線預(yù)測的改進粒子濾波算法.基于天津富民橋主纜、美國I-39北橋第二跨橫橋以及伊通河橋上游主梁的監(jiān)測數(shù)據(jù)對所提算法進行了驗證分析，為證實算法的適用性和魯棒性，進行了假設(shè)試驗驗證，通過假設(shè)試驗可以模擬多種情況下橋梁極值應(yīng)力(或撓度)的預(yù)測分析.

基于監(jiān)測極值應(yīng)力(或撓度)信息對極值應(yīng)力(或撓度)進行了預(yù)測分析，預(yù)測結(jié)果與實測值變化趨勢基本保持一致，預(yù)測結(jié)果的均方誤差也較為理想，且通過對比發(fā)現(xiàn)改進的粒子濾波器預(yù)測精度十分可觀；此外，在橋梁極端情況下，IPF的預(yù)測性能良好，具有良好的預(yù)測穩(wěn)定性以及對突變數(shù)據(jù)的靈敏性與自適應(yīng)性.而對橋梁結(jié)構(gòu)可靠性進行預(yù)測分析時，可靠度指標的一步預(yù)測值與實測值標化趨勢十分吻合，而考慮了隨機性影響所得到的預(yù)測可靠度指標從整體上略微低于實測值，但更符合實際情況.