橋梁耦合極值應力的貝葉斯動態預測

戈 壁，王佐才，2，辛 宇，3，李 舒，丁雅杰，孫曉彤

（1.合肥工業大學土木與水利工程學院，安徽合肥 230009；2.土木工程防災減災安徽省工程技術研究中心，安徽合肥 230009；3.安徽省基礎設施安全檢測與監測工程實驗室，安徽合肥 230009；4.清華大學合肥公共安全研究院，安徽合肥 230601）

引言

隨著橋梁健康監測系統在實際工程中的普及，系統能夠實時監測橋梁結構在服役期內的響應及變形［1-3］。其中，在環境及荷載作用下結構內部產生的應力反映了結構局部的受力狀態，同時也是結構安全運營狀態的重要指標之一。因此，應力監測在橋梁健康監測系統中具有十分重要的地位。應力監測可分為監測數據的收集和應用兩個部分。經過幾十年的研發，數據收集技術都已相對成熟［4-6］。目前國內外學者的研究主要集中在基于歷史監測數據的應用方面，并在橋梁結構損傷識別［7-9］、荷載效應分析建模［10］、模態參數識別［11］等方面取得了一定的成果。然而，在基于橋梁極值應力監測數據的結構安全性能評估與預測方面還有待發展。橋梁極值應力蘊含著結構安全性能的重要信息，對極值應力進行分析以及動態預估，不但能夠動態、直觀地反映在役橋梁運營狀況下的極值應力狀態，為橋梁結構的維護與養護提供重要依據，而且能夠為橋梁的動態安全性能評估提供理論依據。因此，如何合理地利用海量歷史監測數據，實現結構響應的動態預估是結構健康監測目前亟需解決的關鍵問題。

目前，基于實際監測數據對結構響應進行動態預估已經有了一些研究。例如，馮海月等［12］基于廣義Pareto 分布提出了一種改進的獨立風暴法，可較好地預測車輛荷載效應。陽霞等［13-14］改進了廣義Pareto 分布模型的閾值選取方法，并可利用實際監測數據建立具有一定規律的概率統計模型。研究表明改進后的概率統計模型能夠更好地擬合出由車輛荷載引起的橋梁應變右尾分布，結合極值統計理論能夠對橋梁剩余服役期內的應變極值進行估計。Fan 等［15］提出可利用極值應力監測數據以及高斯混合粒子濾波器實現橋梁極值應力的動態預測。但是，由于極值應力是多種荷載效應的耦合，直接對耦合監測數據進行擬合以及動態預估，其效果并不理想。近年來，由于貝葉斯動態模型可結合監測數據的客觀信息和主觀信息，能夠對橋梁極值應力進行有效預測，使該方法得到了長足的發展［16-17］。Ni 等認為在役橋梁由于承受了多種荷載才導致應變/應力響應數據結構的多模態性。為了更好地預測耦合極值應力，分別提出了一種具有較強非平穩動態過程建模能力的貝葉斯動態線性模型［18］和一種參數貝葉斯混合模型［19］。通過試驗驗證，改進后的模型對極值應力的預測精度具有一定程度地提高。Liu等［20］考慮到監測得到的耦合極值應力是非平穩數據與平穩數據的耦合，提出了一種利用局部多項式理論建立的貝葉斯動態線性趨勢模型。該模型能有效地減小監測數據中非平穩趨勢項對預測結果的影響，并隨著預測模型參數的更新，橋梁耦合極值應力的預測結果越來越合理。樊學平等［21］則將橋梁耦合極值應力簡單分為具有趨勢性的溫度荷載效應和具有隨機性的車輛荷載效應并提出了基于解耦極值應力的動態耦合線性模型。模型首先使用簡單的移動平均法數據解耦，再分別建立貝葉斯動態耦合線性模型并進行極值應力預測分析。樊學平等［22］認為由于監測數據具有動態性、隨機性以及周期性等特點，因此可利用Fourier 函數對初始應力狀態數據進行回歸分析，進而運用Fourier 動態線性模型并對橋梁極值應力進行預測。從上述研究可以看出，基于貝葉斯方法的動態模型能夠對橋梁極值應力進行動態預測，并且利用解耦極值應力數據，可使預測模型精度進一步提升。

由于橋梁耦合極值應力主要受到溫度荷載效應、車輛荷載效應和風荷載效應等的影響，因此橋梁極值應力監測數據的解耦方法通常按照所受荷載類別進行解耦并分類，這就使得解耦數據仍可能是多階模態響應的耦合。這樣，利用解耦數據建立的貝葉斯動態模型依然存在著模型相對復雜，不能很好地體現橋梁極值應力周期性、隨機性等數據特點的問題。為了解決上述問題，本文引入了Hilbert 信號分解以及解調技術，對數據進行解耦。該方法不但能夠對具有密集模態與頻率疊混的非平穩結構響應信號進行有效分解，而且能夠很好地保存信號幅度、頻率和相位的信息［23-25］。同時，在得到單分量極值應力數據后，即可利用Hilbert 信號解調技術，回歸得到狀態變量函數的參數，建立單分量極值應力數據的Hilbert 動態線性模型（Hilbert Dynamic Linear Model，HDLM）。隨后，利用貝葉斯方法對模型進行概率遞推，實現單分量極值應力的動態預測。將各階單分量極值應力預測值進行相加即可得到橋梁耦合極值應力的動態預測值。最后，利用在役橋梁一年的監測數據對本文所提方法的有效性進行了驗證。

1 橋梁耦合極值應力的動態預測

橋梁健康監測系統在長期運營期間可以監測大量的應力信息，而結構的應力響應則反映了諸多荷載（如溫度荷載，車輛荷載，風荷載等）作用下結構的內部特征。因此，橋梁結構的應力響應具有復雜的耦合性。為了更好地動態預測橋梁耦合極值應力，本文首先建立了單分量極值應力的HDLM，進而提出了BHDLM。其中，HDLM 由監測方程與狀態方程組成，狀態變量表示了監測數據的變化趨勢，而狀態方程則描述了狀態變量如何隨時間變化，監測方程描述了監測變量與狀態變量之間的關系。其中，預測模型的主要流程如圖1所示。

從圖1 流程圖中可以看出，在數據解耦過程中，定義時間區間（t-1，t］內，監測的最大應力為橋梁耦合極值應力Yt，并將監測得到的極值應力數據視為一個時間序列。在得到Yt后，首先利用了Hilbert信號分解技術對數據進行解耦［23-26］，得到有限多個具有單一頻率成分的單分量極值應力數據yi，t。隨后，基于Hilbert 平方解調算法（Hilbert Square Demodulation，HSD）回歸得到各階單分量極值應力數據yi，t的HDLM。一旦建立了單分量極值應力的HDLM，即可通過貝葉斯方法對模型中的參數實現概率遞推，實現對極值應力的動態預測，進而可得到橋梁耦合極值應力的動態預測結果。

圖1 橋梁耦合極值應力預測模型流程圖Fig.1 The flowchart of the prediction model for bridge coupled extreme stress

2 單分量極值應力的HDLM 及參數確定

2.1 橋梁極值應力解耦

直接利用耦合極值應力進行建模與預測，不但增加了極值應力以及預測模型的復雜性，而且不利于應力預測的精度［20-21］。因此，為了減小監測數據的復雜性，滿足建立預測模型的需要，對極值應力進行解耦應是HDLM 建模過程中的必要步驟。

在以往的研究中，橋梁耦合極值應力被認為是多種荷載效應的耦合，因此解耦過程中通常按荷載效應類型進行分類。一般利用簡單移動平均、二次移動平均等擬合算法，即可得到結構監測應力的趨勢項。趨勢項的均值可看作結構自重荷載效應，而除去均值的趨勢項因其具有較為明顯的周期性，可認為是環境溫度作用下的荷載效應。監測數據與趨勢項的差值則被認為與車輛荷載效應和風荷載效應相關。可以看出上述解耦過程雖然步驟簡單、便于實現，但過程中存在著經驗成分，因此解耦過程并不十分充分。

根據模態疊加理論，結構的位移響應可表示為各階模態貢獻之和。如果只考慮結構某一方向的正應變，則應變與位移之間存在如下關系式：

式中ε為正應變；φi為第i階的位移模態振型；qi為模態坐標。而應力與應變在一定的比例極限范圍內可視為線性關系。則Yt同樣可以表示為具有n個有效頻率ωi（i=1，2，…，n）成分耦合的實數時間序列。為了能夠充分解耦，得到具有單一頻率成分的單分量極值應力數據。本文提出利用Hilbert 信號分解算法對橋梁耦合極值應力進行分解［23-26］。一旦確定了Yt中的頻率，即可利用Hilbert 變換的信號分解算法對Yt進行解耦，分解成n個具有單一頻率的單分量成分yi，t（i=1，2，…，n）。則極值應力Yt可進一步表示為：

2.2 單分量極值應力解調

在利用Hilbert 信號分解算法將橋梁耦合極值應力Yt分解為有限多個單分量極值應力yi，t后，可進一步利用HSD 算法對yi，t進行解調，回歸得到狀態變量θi，t的函數。HSD 算法是一種針對具有單一頻率成分數據的解調算法，其原理可簡單概述如下［27］：

假設單分量極值應力數據yi，t的回歸函數θi，t可簡化為：

式中Ai，Bi，Di，ωγ，i以及ωi均為回歸系數。首先，構建yi，t的解析信號Zi，t：

式中 H［·］為Hilbert 變換算子。則回歸函數θi，t中的li，t可表示為：

聯立公式（3）和（5），可進一步得到hi，t：

其中，為了減小單分量極值應力數據中噪聲以及解調誤差對回歸函數θi，t的影響，可對公式（5）和（6）作進一步優化［27］。則回歸系數Ai可認為是解調結果li，t的均值。li，t中的相位函數ωγ，it+Di可通過下式得到：

其中，系數ωγ，i可通過對相位函數求導得到，進而可回歸得到Di。結合公式（5）和（7），回歸系數Bi可表示為：

而hi，t中的相位函數ωit可表示為：

再對相位函數進行求導可得到回歸系數ωi。

2.3 單分量極值應力的HDLM

根據公式（3）～（9），可得到各階yi，t的單分量極值應力回歸函數θi，t。則θi，t+1的表達式可由下式得到：

式中g（·）表示非線性函數。可以看出，狀態變量隨時間呈非線性關系。利用公式（10）可得到的歷史監測應力數據的動態非線性模型。然而，由于貝葉斯方法更適用于線性模型，因此在得到動態非線性模型后往往需要進行線性化處理［28］。本文利用Taylor 級數展開技術，將g（θi，t）在初始狀態數據的期望值M1，t附近處展開，即可得到公式（10）對應的線性函數：

其中，初始狀態數據可利用Hilbert 信號分解算法得到，則初始狀態數據期望值M1，t可由數據的均值近似得到；為已知項，且與θi，t無關；о（θi，t-M1，t）為g（θt）在M1，t處展開的高階項。公式（11）可進一步改寫為：

式中θi，t+1-θi，t為狀態數據的一階差分；βi，t=g（M1，t）-G1·M1，t+（G1-1）θi，t。通過公式（11）和（12）可以看出，狀態t+1 時刻的狀態變量θi，t+1不僅和前一時刻的θi，t有關，也可能與狀態數據的一階差分變量βi，t有關。考慮βi，t只與前一時刻的一階差分變量βi，t-1有關，并忽略公式（11）和（12）中的高階項，橋梁單分量極值應力的Hilbert 動態線性狀態方程可表示為：

模型誤差方差Wt+1無法直接從狀態數據中得到，可根據初始狀態數據的方差Ct近似計算得到［21］：

式中δ為折扣因子，取值范圍為0.95～0.98。折扣因子δ對預測精度影響很大，選擇合適的δ能夠使貝葉斯動態線性模型的預測精度快速收斂［15］。狀態變量以及狀態一階差分量的方差矩陣Ct可表示為：

其中，方差矩陣中對角元素C1，t，C2，t可分別通過計算初始狀態數據以及初始狀態一階差分數據的方差獲得。θi，t與βi，t的協方差C3，t可由下式近似得到：

式中M2，t為初始狀態數據一階差分結果的期望值，可由數據的均值近似得到。

利用監測數據以及狀態變量θi，t+1，可得到單分量極值應力的Hilbert 動態線性模型中的監測方程。假設監測值與狀態變量呈線性關系，并考慮監測誤差的影響，監測方程可表示為：

式中yi，t+1為t+1 時刻單分量極值應力監測值；vt+1為監測誤差并假設監測誤差服從期望為0，方差為Vt+1的正態分布。根據公式（17）可知，對單分量極值應力yi，t與狀態數據θi，t之間的差值進行統計分析，可近似得到監測誤差的方差Vt+1。

結合公式（11），（15）和（16），t時刻以及t時刻之前初始狀態信息的后驗概率分布可整理為：

式中θi，t為狀態變量向量；為狀態變量和狀態一階差分量的期望向量；Dt表示t以及t時刻之前監測數據的信息集合，包括t時刻的監測數值yi，t，t-1 以及t-1 時刻之前的Mt-1和Ct-1。

3 HDLM 的貝葉斯概率遞推

在建立單分量極值應力的HDLM 之后，結合貝葉斯方法可得到BHDLM。根據t時刻狀態變量θi，t的后驗分布可對t+1 時刻的極值應力yi，t+1進行預測。再利用貝葉斯方法以及實際監測數據對模型中t+1 時刻的狀態變量后驗分布θi，t+1概率進行修正，使得下一時刻的預測結果更加合理。這樣，隨著模型中的參數不斷更新，BHDLM 即可實現對單分量極值應力的動態預測并且精度不斷提高。模型的概率遞推過程如圖2所示。

圖2 BHDLM 的概率遞推過程Fig.2 The probability recursive process of BHDLM

BHDLM 遞推過程中參數更新如下所示［19-22，27］：

（1）t+1 時刻的狀態變量先驗概率密度分布

結合公式（18）給出了狀態變量θi，t在t時刻的后驗概率密度分布。根據狀態方程可得到t+1 時刻的θi，t+1和βi，t+1先驗概率密度分布，表示為：

簡化可得：

式中

其中，R1，t+1=C1，t+2C3，t+C2，t+W1，t+1，R2，t+1=C2，t+W2，t+1，R3，t+1=C3，t+C2，t+W3，t+1。

（2）t+1 時刻監測變量的一步向前預測概率密度分布P（yi，t+1|Dt）

結合狀態變量的先驗概率P（θi，t+1|Dt）以及監測方程，即可得到監測變量一步向前預測概率密度：

簡化可得監測變量概率密度：

式中 預測值ft+1=at+1（1，1）；預測方差Qt+1=Rt+1（1，1）+Vt+1。而預測方差的倒數則為BHDLM的預測精度［21］。

（3）t+1 時刻的狀態變量后驗概率

密度分布P（θt+1|Dt+1）

當已知t+1 時刻的監測值Yt+1后，結合公式（20）和（22）以及貝葉斯方法，t+1 時刻的狀態變量θt+1后驗概率密度分布可表示為：

簡化可得：

式中 期望向量Mt+1中M1，t+1=at+1（1，1）+A1，t+1·et+1，M2，t+1=at+1（1，2）+A2，t+1·et+1，其中，適應性系數A1，t+1=Rt+1（1，1）/Qt+1，A2，t+1=Rt+1（1，2）/Qt+1，一步預測誤差et+1=yi，t+1-ft+1；方差矩陣Ct+1中C1，t+1=A1，t+1·Vt+1，C2，t+1=Rt+1（2，2）-A2，t+1·Rt+1（1，2），C3，t+1=A2，t+1·Vt+1。

4 實橋分析

派河大橋主橋上部結構為飛燕式鋼箱拱橋，跨徑組合為54 m+130 m+54 m。拱圈和鋼梁采用Q345D 鋼材，彈性模量為206 GPa。為了全面監測橋梁在運營期間的極值應力應變，全橋共選取5 個關鍵截面布置了32 個應變傳感器，應變監測截面布置圖如圖3所示。其中，鋼縱梁的應變傳感器安裝在鋼縱梁內部，每個截面安裝4 個傳感器。為了了解拱腳處鋼縱梁的應力狀態，定義每小時監測最大應力為橋梁耦合極值應力Yt，并提取了12#墩附近2-2 截面Sx9～Sx12 監測點一年（共8571 h）的歷史監測數據。其中，2-2 截面鋼縱梁應力監測點布置如圖4所示。

圖3 應力監測截面布置圖（單位：mm）Fig.3 The layout diagram of stress monitored section（Unit：mm）

4.1 基于BHDLM 的橋梁耦合極值應力分析

為了驗證本文所提方法的有效性，現選取Sx12測點歷史監測數據中前3000 h 作為建模數據，建立各階單分量極值應力的HDLM。再利用建立的模型以及貝葉斯方法，對剩余5572 h 的極值應力進行動態預測。其中，Sx12 測點采集的歷史監測數據如圖5所示。

對建模數據進行傅里葉變換，結果如圖6所示。從圖6 中可以看出，歷史監測數據是多模態響應耦合的數據。

圖6 Sx12 測點極值應力數據的傅里葉譜圖Fig.6 The Fourier spectrum of the extreme stress data of the Sx12 measure point

結合傅里葉分析結果和Hilbert 信號分解算法對前3000 h 耦合極值應力進行解耦，可得到4 個具有單一頻率的單分量成分yi，t（i=1，2，3，4）。基于公式（3）～（9），分別對yi，t進行解調，即可得到單分量極值應力的狀態數據回歸函數θi，t（i=1，2，3，4）。解調結果為：

首先，基于y1，t以及θ1，t，對橋梁耦合極值應力中的一階單分量應力數據進行預測。根據θ1，t的回歸公式可得θ1，t+1：

則基于公式（25），（26）和（12）可知，θ1，t≈θ1，t+1，即可認為θ1，t+1僅和前一時刻的θ1，t有關，與狀態數據的 一 階 差 分 變 量β1，t無 關，且G1=。其中，監測應力和初始狀態信息如圖7所示。

圖7 第一階解耦應力數據以及初始狀態信息Fig.7 The first order decoupled stress data and initial state information

則初始狀態信息為：

則根據公式（13），狀態方程可簡化為：

其中，折扣因子δ的取值為0.95。監測方程為：

在得到0～2999 h 監測數據的動態線性模型后，結合公式（19）～（24），可對3000～8573 h 的數據進行動態預測，一步預測結果以及預測精度分別如圖8 和9所示。從圖8 的預測結果中可以看出，第一階單分量極值應力的預測值與監測數據近似相等，證明本文所提模型預測效果較好。而結合圖9 的預測精度結果可以看出，預測應力區間可以包含所有的監測應力數據，隨著單分量極值應力數據的遠程更新，精度越來越高并最終趨于穩定，說明本文所提模型的合理性。

圖8 第一階解耦極值應力數據預測結果Fig.8 The predicted results of the first order decoupled extreme stress data

圖9 BHDLM 的預測精度Fig.9 The prediction precision of BHDLM

同樣的，基于y2，t以及θ2，t，對橋梁耦合極值應力中的二階單分量應力數據進行預測。根據θ2，t的回歸公式可得θ2，t+1：

基于公式（25），（30）和（12）可知，監測方程中θ2，t+1與前一時刻的θ2，t和狀態數據的一階差分變量β2，t有關。利用差分計算可得狀態變量轉移矩陣中其中，監測應力和初始狀態信息如圖10所示，初始狀態信息的一階差分值如圖11所示。

圖10 第二階解耦應力數據以及初始狀態信息Fig.10 The second order decoupled stress data and initial state information

圖11 初始狀態信息的一階差分值Fig.11 The first order differential data of initial state information

則根據圖10 和11，初始狀態信息可表示為：

其中

而根據公式（13），狀態方程為：

其中，折扣因子δ的取值范圍為 0.95，。

監測方程為：

利用貝葉斯方法，對3000～8572 h 的數據進行動態預測，預測結果如圖12所示，預測值與監測值接近，預測效果較好。預測精度如圖13所示，隨著數據的更新，預測精度趨于穩定。

圖12 第二階解耦極值應力數據預測結果Fig.12 The predicted results of the second order decoupled extreme stress data

圖13 BHDLM 的預測精度Fig.13 The prediction precision of BHDLM

將第3和4階單分量極值應力數據按上述步驟進行動態預測，則橋梁耦合極值應力即為各階單分量極值應力動態預測值相加，預測結果如圖14所示。

圖14 Sx12 測點橋梁耦合極值應力預測結果Fig.14 The predicted results of bridge coupled extreme stress of the Sx12 measure point

按照上述流程，利用BHDLM 可進一步得到Sx9，Sx10 和Sx11 這3 個測點3000 h 以后的橋梁耦合極值應力動態預測結果，預測結果如圖15～17所示。可以看出，耦合后的極值應力動態預測值與實際監測值接近，說明了本文所提方法能夠較好地預測橋梁耦合極值應力。

圖15 Sx9 測點橋梁耦合極值應力預測結果Fig.15 The predicted results of bridge coupled extreme stress of the Sx9 measure point

4.2 耦合極值應力預測準確度對比分析

圖16 Sx10 測點橋梁耦合極值應力預測結果Fig.16 The predicted results of bridge coupled extreme stress of the Sx10 measure point

圖17 Sx11 測點橋梁耦合極值應力預測結果Fig.17 The predicted results of bridge coupled extreme stress of the Sx11 measure point

為了驗證BHDLM 的在預測橋梁耦合極值應力方面的優越性，本文利用貝葉斯動態線性模型（Bayesian Dynamic Linear Model，BDLM）對Sx12測點處的極值應力進行了動態預測。同時，采用均方根誤差（Root Mean Squared Error，RMSE）作為模型預測準確度指標，用于模型預測準確度的比較分析。其中，預測準確度指標RMSE可表示為：

式中uy，t為橋梁耦合極值應力預測值；Yt為橋梁耦合極值應力監測值。

BDLM 預測結果如圖18所示。根據圖18 的預測結果并結合公式（34），可得BDLM 預測準確度指標RMSE為2.59。而利用本文提出的BHDLM 進行動態預測時，模型預測準確度指標RMSE則為1.60。證明了本文所提方法相比于BDLM 具有更高的預測準確度，更具優越性。

圖18 利用BDLM 得到的橋梁耦合極值應力預測結果Fig.18 The predicted results of bridge coupled extreme stress by using BDLM

5 結論

橋梁耦合極值應力可視為有限多個單分量極值應力的累加。在考慮到單分量極值應力周期性、隨機性等數據特點的情況下，提出了單分量極值應力的BHDLM。通過對單分量極值應力的動態預測，進而可得到橋梁耦合極值應力的預測結果。最后，利用派河大橋一年的實際監測數據進行分析驗證，得到如下結論：

（1）通過Hilbert 信號分解及解調技術可對橋梁耦合極值應力數據進行解耦，并建立各階單分量極值應力數據的HDLM。

（2）基于解耦得到的單分量極值應力數據以及BHDLM，可對單分量極值應力進行合理預測，預測精度隨著模型參數的更新不斷提升。

（3）橋梁耦合極值應力的最終預測結果與監測數據基本一致，具有較好的預測準確度。驗證了本文提出方法在實際工程中的可行性與有效性。