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

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

0 引言

橋梁健康監(jiān)測(cè)系統(tǒng)在長期運(yùn)營中積累了大量監(jiān)測(cè)數(shù)據(jù)，它包括均勻數(shù)據(jù)和非均勻數(shù)據(jù)，如何合理地處理這些數(shù)據(jù)是橋梁結(jié)構(gòu)健康監(jiān)測(cè)領(lǐng)域的主要困難之一，目前大量研究主要是基于均勻數(shù)據(jù)集中在模態(tài)參數(shù)識(shí)別、損傷識(shí)別、損傷評(píng)估、模型修正等領(lǐng)域[1-3].這些研究仍難以有效地預(yù)測(cè)和評(píng)估橋梁結(jié)構(gòu)的動(dòng)態(tài)響應(yīng)行為，如： 預(yù)測(cè)均勻監(jiān)測(cè)極值應(yīng)力和非均勻監(jiān)測(cè)極值應(yīng)力.

本研究采用貝葉斯動(dòng)態(tài)線性模型(Bayesian dynamic linear model，BDLM)在橋梁承載力和可靠性預(yù)測(cè)方面已取得一些研究成果，如采用動(dòng)態(tài)線性模型實(shí)現(xiàn)橋梁退化承載力的模擬預(yù)測(cè)[4]； 采用應(yīng)力的單一和組合貝葉斯動(dòng)態(tài)線性模型以及一次二階矩可靠度方法實(shí)現(xiàn)橋梁構(gòu)件的可靠性動(dòng)態(tài)預(yù)測(cè)[4-5].但在預(yù)測(cè)分析過程中存在以下問題需要進(jìn)一步解決： ① 最優(yōu)初始狀態(tài)概率分布如何選擇； ② 非均勻極值應(yīng)力的動(dòng)態(tài)建模問題； ③ 貝葉斯動(dòng)態(tài)線性模型的監(jiān)控機(jī)制問題.

鑒于上述問題，本研究基于非均勻監(jiān)測(cè)得到應(yīng)力信息并建立貝葉斯動(dòng)態(tài)線性模型.考慮到初始狀態(tài)信息概率分布的多樣性，采用K-L準(zhǔn)則選擇最優(yōu)初始狀態(tài)概率分布函數(shù)，并利用折扣因子、單一貝葉斯因子以及累積貝葉斯因子提高和監(jiān)控動(dòng)態(tài)模型的預(yù)測(cè)精度，最后采用實(shí)例驗(yàn)證了所提模型和方法的合理性和適用性.

1 監(jiān)測(cè)極值應(yīng)力的動(dòng)態(tài)線性模型

本研究的極值應(yīng)力是從結(jié)構(gòu)健康監(jiān)測(cè)系統(tǒng)監(jiān)測(cè)得到.其為每天所測(cè)應(yīng)力的最大值信息.由這些極值信息回歸得到的二次函數(shù)[6]，如式(1)所示，比線性函數(shù)[7]可以更合理地反映應(yīng)力極值信息的趨勢(shì)項(xiàng)變化規(guī)律，因而可以用來構(gòu)造狀態(tài)水平方程，如式(2)～(4)所示，極值應(yīng)力的動(dòng)態(tài)建模過程如圖1所示.

h(t)=at2+bt+c

(1)

θt+1=h(t+1)=a(t+1)2+b(t+1)+c

(2)

θt=h(t)=at2+bt+c

(3)

θt-1=h(t-1)=a(t-1)2+b(t-1)+c

(4)

圖1 極值應(yīng)力的貝葉斯動(dòng)態(tài)線性建模 Fig.1 Bayesian dynamic linear modeling of extreme stress

式中:h(t)為由極值信息回歸得到的二次函數(shù)；a，b和c為監(jiān)測(cè)極值應(yīng)力的回歸系數(shù)；t為監(jiān)測(cè)時(shí)間，d；θt為第t天的極值應(yīng)力狀態(tài)值.

由式(2)～(4)，進(jìn)一步可得

(5)

式中:Δt+1為t時(shí)刻到t+1時(shí)刻狀態(tài)水平值的變化量.

考慮到監(jiān)測(cè)極值應(yīng)力的隨機(jī)性和不確定性，基于狀態(tài)信息和狀態(tài)變化信息建立的非均勻極值應(yīng)力動(dòng)態(tài)線性模型如式(6)～(8)所示.式中的應(yīng)力狀態(tài)和狀態(tài)的變化信息分別是指監(jiān)測(cè)信息經(jīng)過重采樣方法(如： 五點(diǎn)三次平滑處理方法)得到的極值應(yīng)力數(shù)據(jù)和重采樣極值應(yīng)力的變化數(shù)據(jù).

監(jiān)測(cè)方程：

yt+1=θt+1+νt+1,νt+1～N[0,Vt+1]

(6)

狀態(tài)方程：

θt+1=θt+Δt+1+ωt+1, 1,Δt+1=(ri+1/ri)(Δt+2a)+ωt+1, 2,

(7)

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

(8)

式中：

yt+1為t+1時(shí)刻的應(yīng)力監(jiān)測(cè)極值；θt+1為t+1時(shí)刻極值應(yīng)力的狀態(tài)水平值；N[0·Wt+1]為正態(tài)分布函數(shù)；νt+1為監(jiān)測(cè)誤差項(xiàng)；ωt+1為狀態(tài)噪聲向量/狀態(tài)誤差向量；Dt為t時(shí)刻以及t時(shí)刻以前關(guān)于系統(tǒng)的信息集合，即：Dt={yt,Dt-1},Dt-1為t-1時(shí)刻的信息集，包括mt-1(平均值向量)，Ct-1(方差矩陣)等.另外，假設(shè)vt+1和wt+1相互獨(dú)立，且都與θt+1獨(dú)立，ri+1=ti+1-ti為采樣時(shí)間間隔.

2 貝葉斯動(dòng)態(tài)線性模型主要概率參數(shù)的確定方法

貝葉斯動(dòng)態(tài)線性模型需要確定的主要參數(shù)包括Vt+1，Wt+1，mt和Ct.Vt+1可以通過監(jiān)測(cè)極值應(yīng)力和重采樣極值應(yīng)力的數(shù)據(jù)差近似估計(jì).Wt+1可以通過折扣因子近似得到[6-7]，如式(9)所示.

(9)

對(duì)于狀態(tài)變量θt+1，有K-L信息距離：

(10)

式中:pRT, i(θt+1)為第i個(gè)關(guān)于θt+1的常規(guī)概率分布；pR(θt+1)為核密度估計(jì)得到的θt+1理論概率分布.但是由于核密度得到的理論概率度函數(shù)外延性不好, 所以要選擇合適的常規(guī)概率密度函數(shù)來分析.

對(duì)于狀態(tài)變化變量Δt+1，亦可以通過式(10)選擇關(guān)于Δt+1最合理的常規(guī)概率分布函數(shù).

3 動(dòng)態(tài)線性模型的貝葉斯遞推

圖2 貝葉斯動(dòng)態(tài)線性模型的概率遞推過程Fig.2 Probability recursion processes of BDLM

貝葉斯動(dòng)態(tài)預(yù)測(cè)模型[8-9]適合于對(duì)未來狀態(tài)參數(shù)的預(yù)測(cè)和推斷，其狀態(tài)參數(shù)遞推如圖2所示.

貝葉斯動(dòng)態(tài)線性模型的詳細(xì)遞推步驟：

1)t時(shí)刻狀態(tài)變量的后驗(yàn)分布.

對(duì)于均值mt和方差矩陣Ct，有：

(11)

2)t+1時(shí)刻狀態(tài)變量的先驗(yàn)分布.

(μt+1|Dt)～N[at+1,Rt+1]

(12)

3)t+1時(shí)刻監(jiān)測(cè)變量的一步預(yù)測(cè)分布.

(yt+1|Dt)～N[ft+1,Qt+1]

(13)

式中:ft+1=E(yt+1|Dt)=mθt+(ri+1/ri)(mΔt+2a)，Qt+1=Var(yt+1|y1:t)=Rt+1, 1+Vt+1，pt+1=1/Qt+1為BDLM的預(yù)測(cè)精度.

4)t+1時(shí)刻狀態(tài)變量的后驗(yàn)分布.

(14)

根據(jù)HPD區(qū)域的定義[1-2, 5]，對(duì)于監(jiān)測(cè)值的預(yù)測(cè)區(qū)間(95%的保證率)為：

(15)

4 非均勻貝葉斯動(dòng)態(tài)線性模型的監(jiān)控機(jī)制

本研究主要通過建立貝葉斯因子對(duì)動(dòng)態(tài)模型進(jìn)行監(jiān)控[4-9]，包括單一貝葉斯因子和累積貝葉斯因子.單一貝葉斯因子H(t)為：

(16)

式中:p0(yt|Dt-1)為一步預(yù)測(cè)監(jiān)測(cè)變量的概率密度函數(shù)；p1(yt|Dt-1)為備擇模型的PDF.

基于式(16)，可得累積貝葉斯因子為：

(17)

式中:Ht(k)是基于監(jiān)測(cè)信息yt,yt-1, …,yt-k+1的p0(·)相對(duì)于p1(·)的k步累積貝葉斯因子.

式(16)進(jìn)一步簡(jiǎn)化可得單一貝葉斯因子的具體表達(dá)式為：

(18)

模型監(jiān)控準(zhǔn)則[4-9]為：k=3時(shí)，H(t)小于0.15，則此檢測(cè)值為異常值，需去除.異常值去除之后，可以通過累計(jì)貝葉斯因子的變化曲線來體現(xiàn)貝葉斯動(dòng)態(tài)模型的預(yù)測(cè)精度，累計(jì)貝葉斯因子越大，說明預(yù)測(cè)精度越好，模型的不確定性越小.

5 算例分析

I-39北橋[10-11]建于1961年，對(duì)此橋橫向第二跨鋼板的跨中梁底極值應(yīng)力進(jìn)行了83 d的監(jiān)測(cè)，其中前60 d為均勻采集，后23 d為非均勻采集.此位置應(yīng)力集中比較明顯，所以監(jiān)測(cè)的應(yīng)力數(shù)據(jù)，可以較好把握結(jié)構(gòu)的安全狀態(tài)[10-11].應(yīng)力是由溫度荷載、車輛荷載以及收縮徐變?cè)斐傻? 橋梁自重造成的應(yīng)力不包括在此應(yīng)力信息中.前60 d均勻采集的數(shù)據(jù)如圖3所示.

圖3 平滑處理后的極值應(yīng)力與監(jiān)測(cè)極值應(yīng)力對(duì)比 圖4 平滑處理后的極值應(yīng)力與對(duì)應(yīng)極值應(yīng)力差值對(duì)比

現(xiàn)采用前60 d均勻采集的極值應(yīng)力數(shù)據(jù)，通過建立非均勻極值應(yīng)力貝葉斯動(dòng)態(tài)線性模型(如式(19)～(21)所示)，預(yù)測(cè)后23 d非均勻的極值應(yīng)力信息.

為了得到初始狀態(tài)信息的分布參數(shù)，采用五點(diǎn)三次平滑處理方法對(duì)前60 d的極值應(yīng)力進(jìn)行重采樣.然后近似得到極值應(yīng)力狀態(tài)水平值和相對(duì)應(yīng)的狀態(tài)水平差值，如圖4所示.

對(duì)水平值和水平差值分別進(jìn)行了K-S檢驗(yàn), 可知水平值和水平差值都分別服從正態(tài)分布和對(duì)數(shù)正態(tài)分布，但通過K-L信息測(cè)度計(jì)算(式(10)所示)，水平值和水平差值的最優(yōu)概率分布均為正態(tài)分布.因此建立的初始狀態(tài)信息的概率分布如式(21)所示.

本算例基于式(1)～(8)，建立的貝葉斯動(dòng)態(tài)線性模型為

監(jiān)測(cè)方程：

yt+1=θt+1+νt+1,νt+1～N[0,Vt+1]

(19)

狀態(tài)方程：

θt+1=θt+Δt+1+ωt+1, 1,Δt+1=(ri+1/ri)(Δt-0.002 294)+ωt+1, 2,

(20)

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

(21)

基于式(1)～(8)、式(11)～(15)、式(19)～(21)和圖2，可得預(yù)測(cè)的非均勻極值應(yīng)力和預(yù)測(cè)精度分別如圖5～6所示，動(dòng)態(tài)模型的預(yù)測(cè)應(yīng)力與監(jiān)測(cè)應(yīng)力基本一致，而且預(yù)測(cè)精度越來越好.

圖5 監(jiān)測(cè)和預(yù)測(cè)得到的極值應(yīng)力數(shù)據(jù) 圖6 動(dòng)態(tài)模型的預(yù)測(cè)精度

基于式(16)～(21)，可知單一和累計(jì)貝葉斯因子如圖7～8所示.從圖7可以看出單一貝葉斯因子均大于0.15，說明監(jiān)測(cè)的極值應(yīng)力均為正常值，即無異常值.從圖8可以看出，隨著監(jiān)測(cè)應(yīng)力的不斷修正，累積貝葉斯因子越來越大，說明動(dòng)態(tài)模型的預(yù)測(cè)精度越來越好，模型的預(yù)測(cè)不確定性越來越小.

圖7 單一貝葉斯因子 圖8 累計(jì)貝葉斯因子

6 結(jié)論

本研究提出了一種預(yù)測(cè)非均勻極值應(yīng)力的貝葉斯動(dòng)態(tài)方法，并給出對(duì)應(yīng)動(dòng)態(tài)模型的監(jiān)控機(jī)制.通過實(shí)例驗(yàn)證得出： 預(yù)測(cè)的非均勻極值應(yīng)力和監(jiān)測(cè)的極值應(yīng)力基本一致，而且預(yù)測(cè)精度越來越好, 并且監(jiān)控機(jī)制對(duì)動(dòng)態(tài)模型進(jìn)行了很好的監(jiān)控分析.本研究只提出橋梁構(gòu)件非均勻極值應(yīng)力的貝葉斯動(dòng)態(tài)預(yù)測(cè)方法，而對(duì)于橋梁體系極值應(yīng)力的動(dòng)態(tài)預(yù)測(cè)方法有待進(jìn)一步研究分析.