基于貝葉斯理論的非線性結構模型修正及其動力可靠度分析

丁雅杰，王佐才,2,3，辛 宇,3，戈 壁，袁子青

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

結構在遭受強風、地震、海嘯等極端荷載作用時會呈現出明顯的非線性特性以及復雜的非線性動力學行為，非線性廣泛存在于實際工程結構中，傳統的基于線性模型及響應平穩性假設理論的模型修正方法難以實現非線性結構的模型修正，因此，研究非線性結構模型修正就是要對極端荷載作用下的結構劣化行為及安全性能進行有效評估。

在結構的非線性模型修正過程中，不可避免地受測試噪聲、環境擾動和單元離散化等不確定性誤差的影響，利用實測響應數據對非線性模型參數的不確定性進行量化，是提高非線性模型可靠性的重要手段?；谪惾~斯推理的模型修正方法結合了結構的先驗信息，并利用后驗概率分布量化修正后的模型參數的不確定性，因此，在有限元模型修正領域得到了廣泛應用。例如，BECK 和AU[1]首次利用貝葉斯推理方法對線性結構的有限元模型進行修正，并利用Metropolis-Hastings (MH)算法實現了對剪切框架模型層間剛度的不確定性量化；易偉建等[2]考慮到多維參數的后驗積分求解困難，提出利用Markov Chain Monte Carlo (MCMC)方法實現了對鋼筋混凝土框架結構的損傷分析；劉綱等[3]和張建新[4]基于廣義無偏先驗分布，將相關向量機與MCMC 方法相結合，提高了后驗樣本的計算效率，實現了有限元模型修正的快速計算；WAN 和REN[5]基于貝葉斯框架，利用高斯過程模型實現對一座人行天橋的有限元模型修正，并量化由于測量噪聲所引起的模型參數的不確定性；XIN 等[6]利用結構動力響應主分量的瞬時特征參數構建目標函數，并利用Cramér-Rao 邊界理論實現對非線性模型參數的不確定性量化。

針對非線性結構模型修正問題，很大程度上取決于非線性特征指標的選取[7?8]。傳統的非線性指標包括振動響應時程、動力響應峰值的包絡以及時變模態等，難以全面地反映非線性結構的時變物理特性，因此，構造合適的非線性特征指標是進行非線性模型修正的首要任務。此外，有限元模型修正的最終目的是要定量地評價結構的可靠性，即在概率意義上定量地評價結構的安全程度。利用工程結構在服役期間的實測響應對結構的模型參數進行修正，并對其動力可靠度進行重新評估，能夠較為真實地反映結構的實際安全狀態[9]。根據隨機振動的首次超越破壞準則，進行復合隨機振動系統動力可靠度分析，是實現結構安全性評估的主要方法之一[10?12]。然而，基于跨越分析的動力可靠度理論需要同時對結構響應的聯合概率密度函數和跨越過程進行假定(如Poisson假設或Markov 假設)，導致結構動力可靠度計算結果與真實值產生一定偏差。

考慮到結構振動響應主分量的瞬時加速度幅值反映了結構響應信號的幅值信息，瞬時頻率則隱含了振動響應的相位信息，因此，為了實現結構的非線性模型修正及其動力可靠度計算，本文基于貝葉斯推理方法，利用實測結構動力響應主分量的瞬時特征參數作為非線性指標構建似然函數，結合拒絕延緩自適應(Delayed Rejection and Adaptive Metropolis, DRAM)算法和高斯過程模型用于非線性模型參數后驗概率密度函數的快速計算。同時，利用結構動力響應的概率密度演化方法[13?14]，根據首超破壞準則對概率密度演化方程施加邊界約束條件，采用Newmark-Beta 數值積分法和總變差遞減(Total Variation Diminishing, TVD)差分格式進行非線性隨機結構的動力可靠度計算。最后，通過數值模擬對比隨機荷載作用下確定性模型和非線性隨機模型的動力可靠度計算結果，并應用蒙特卡洛隨機抽樣驗證該方法的準確性。

1 非線性結構模型修正及其動力可靠度分析理論

1.1 基于貝葉斯理論的非線性結構模型修正

隨機有限元模型修正實質上是一個逆向不確定性量化過程，根據實測結構輸出響應的不確定性逆推參數的不確定性，此過程用貝葉斯公式可表示為[15]：

式中：a和b分 別為參數 θ的下、上界限值。假設Y?(θ)為由瞬時加速度幅值和瞬時頻率組成的實測結構動力響應，Y(θ)為有限元模型計算的理論值，基于貝葉斯方法可建立如下模型修正表達式：

式中： θ為待修正模型參數組成的向量，θ=(θ1,θ2,···,θm)T∈Rm，m為 參數個數； ε為考慮模型誤差、測量噪聲等不確定性因素產生的誤差向量，通?？杉俣榉牧憔档母咚狗植?，記為ε ～N(0, Σε) ， Σε為協方差矩陣。對于進行s次獨立重復試驗的實測結構動力響應樣本，其似然函數可表示為：

1.1.1 DRAM 算法

基于貝葉斯理論的后驗概率密度函數形式復雜且與實測響應數據有關，對于具有高維參數空間的非線性結構模型修正問題，其解析解難以直接利用積分方法獲得，因而，通常采用MCMC 抽樣方法來近似估計待修正參數的后驗分布[17]。對于傳統的MCMC 方法，如MH 算法，當建議分布函數產生的候選樣本參數 θ?被拒絕后，則后驗樣本將保持不變，采樣出現“停滯”。此外，后驗樣本的收斂性與建議分布的協方差函數有關，對于偏離模型參數真實分布的建議分布函數，會導致后驗樣本燃燒期過長、不易收斂，降低了采樣效率。為了提高后驗樣本的接受率、加速Markov鏈的收斂，HAARIO 等[18]在AM 算法的基礎上結合DR 策略，提出了具有自適應策略的DRAM 算法，主要包括以下兩方面內容：

1)設置多重建議分布函數。根據當前參數θt利用建議分布函數q1(θt,θ?) 生 成候選樣本 θ?，并計算樣本接受率：

2)自適應更新建議分布的協方差函數。根據已生成的后驗樣本改變建議分布的協方差矩陣Ct，通過不斷更新建議分布函數來逼近參數的真實概率分布，從而提高采樣效率、縮短燃燒期，此過程可表示為：

式中：C0為建議分布初始協方差；n0為非自適應長度；sm為僅與待修正參數維度m有關的比例因子； cov(θ0,θ1,···,θt?1)為前(t?1)個已生成的后驗樣本協方差矩陣； ε 為 常數；Im為m維單位矩陣。DRAM 算法流程如圖1 所示，其中建議分布函數個數設為q=2。

圖1 DRAM 算法流程圖Fig. 1 Flow chart of DRAM algorithm

1.1.2 高斯過程模型

DRAM 算法通過全局自適應策略提高了后驗樣本的遍歷性，但仍存在計算量大的缺點，每組待修正參數后驗樣本的迭代均需調用有限元模型進行結構動力分析，因而計算效率十分低下。為了實現后驗概率密度函數的快速計算，本文采用高斯過程模型替代有限元模型進行結構響應計算，從而大幅提高后驗樣本的計算效率。

高斯過程模型是基于貝葉斯框架理論發展起來的一種非參數模型，具有計算速度快、擬合精度高以及能夠定量預測不確定性等特點[19]。高斯過程模型完全由均值函數和協方差函數所確定，記為：

式中：m(x) 為高斯過程均值函數；k(x,x′)為協方差函數。本文采用零均值的高斯過程模型，協方差函數包含自相關項、線性回歸項、常數項以及噪聲項，即：

由式(13)可知，基于訓練樣本可得到測試點x?處的預測值f?的均值m??及其方差k??。因此，對于DRAM 算法生成的候選樣本 θ?，根據式(14)可計算得到相應的動力響應預測值Y(θ?)，利用高斯過程模型預測非線性結構動力響應，避免了耗時的非線性有限元模型時程分析，因而，能夠大大提高非線性結構模型修正效率。

1.2 基于概率密度演化的復合隨機結構動力可靠度分析

有限元模型修正的最終目的是要定量地評價結構的可靠性，即在概率意義上定量地評價結構的安全程度。利用實測動力響應數據對結構有限元模型進行修正，并對其動力可靠度進行重新評估，能夠較為真實地反映結構的實際安全狀態。

對于考慮激勵和結構參數隨機性的非線性隨機動力系統，其隨機變量記為 Θ=(θ1,θ2,···,θs)，其中s為隨機變量個數，則多自由度結構運動方程可記為：

式中：X¨(t)、X˙(t)、X(t)分別為結構響應加速度、速度和位移向量；M為質量矩陣；f為包含非線性阻尼力和恢復力的結構內力向量；F(Θ,t)為隨機激勵向量。對于一般的結構振動系統，式(15)有依賴于隨機變量 Θ的唯一解，設系統的狀態量Z(t)記為：

式中：pZΘ(z,θ,t) 為 向量 [Z(t),Θ]的聯合概率密度函數； ?t、 ?Θ分別為振動系統的時間域和 Θ的參數空間。進一步求偏導可得廣義概率密度演化方程[20 ? 21]：

式中，Z˙(θ,t)=?H(θ,t)/?t為Z(t) 在 { Θ=θ}條件下的速度；式(18)的初始條件為：

對于式(18)和式(19)構成的偏微分方程，可利用TVD 格式的差分方法進行求解[22]，通過積分得到隨時間變化的結構動力響應演化概率密度函數pZΘ(z,θ,t)。

基于首超破壞準則的結構動力可靠度可記為[23]：

式中， ?s為安全域。當發生事件 {Xlim(τ)∈?f}時，則結構失效， ?f為失效域。首超破壞準則條件下一旦結構響應超出安全域邊界 ??s進入失效域 ?f，則該事件的概率不再返回安全域，其留在安全域內的總概率即為結構的動力可靠度。因此，對式(18)施加位移邊界約束條件：

結合邊界條件式(21)、初始條件式(19)可求解廣義概率密度演化方程式(18)，進而得到結構的動力可靠度：

2 非線性結構模型修正

2.1 模型概況

本文以具有非線性Bouc-Wen 模型的層間剪切型框架結構為例，根據貝葉斯推理方法實現結構非線性模型修正及其參數不確定性量化。如圖2所示，該結構為五層剪切框架，層高3m，結構質量自頂層至底層依次為7.1×104kg、3.1×105kg、4.8×105kg、5.1×105kg、5.3×105kg。柱截面尺寸為500 mm×500 mm，初始彈性模量E=3.2×104MPa，其單元恢復力曲線為圖3 所示的Bouc-Wen 模型，數學表達式為：

圖2 五層非線性剪切模型 /mFig. 2 Five-layer nonlinear shear model

圖3 Bouc-Wen 模型Fig. 3 Bouc-Wen model

式中：k為結構的彈性剛度；z為滯變位移；u為相對位移； α為水平剛度屈服比；A、 β 、 γ、μ為Bouc-Wen 模型滯回環控制參數，非線性Bouc-Wen模型參數的初始設計名義值如表1 所示。

表1 非線性Bouc-Wen 模型參數的名義值Table 1 Nominal value of nonlinear Bouc-Wen model parameters

2.2 高斯過程替代模型建立

取非線性Bouc-Wen 模型參數作為待修正參數，其名義值如表1 所示。利用Opensees 有限元軟件計算該剪切框架在峰值加速度PGA=0.20g的Northridge 地震動作用下的結構動力響應，其頂層加速度時程曲線如圖4 所示。利用離散解析模式 分 解(Discrete Analytical Mode Decomposition,DAMD)方法[24]和Hilbert 變換對頂層加速度響應進行信號分解，并提取結構振動響應主分量的瞬時加速度幅值和瞬時頻率慢變成分。由于主分量的瞬時特征參數相比于加速度響應隨時間變化較為緩慢，因此，分別選取瞬時加速度幅值和瞬時頻率慢變成分的30 個局部峰值點作為非線性指標來構建似然函數，如圖5 所示。為了提高非線性結構模型修正計算效率，首先利用高斯過程模型建立非線性模型參數與結構瞬時特征參數的回歸關系，定義非線性模型參數在名義值的80%范圍內波動，故待修正參數[ α ,A,μ, β , γ]的上下界限值分別為[0.72, 1.8, 5.4, 0.9, 0.9]和[0.08, 0.2, 0.6,0.1, 0.1]。本算例中，根據Sobol 序列采樣法抽取200 組Bouc-Wen 模型參數訓練樣本，并利用Opensees 有限元模型計算其對應的結構主分量瞬時特征參數。采用留一法交叉驗證高斯過程模型的擬合精度，以第10 s 處結構瞬時加速度幅值為例，其相應的交叉驗證殘差的正態分位圖如圖6所示。由圖6 可知，留一法交叉驗證殘差的正態分位圖趨于直線，表明殘差的正態特性好，因而可判斷建立的高斯過程模型能夠很好地預測結構的瞬時特征參數。

圖4 結構頂層的加速度響應Fig. 4 Acceleration response of top floor

圖5 結構主分量瞬時特征參數Fig. 5 Instantaneous characteristics of principal component response

圖6 殘差正態分位圖Fig. 6 Normal quantile-quantile plot

2.3 考慮不確定性的結構非線性模型修正

2.3.1 考慮測量不確定性的非線性模型修正

為了考慮測量不確定性的影響，對非線性模型參數名義值下的頂層加速度響應，施加30 組噪信比為5%的高斯白噪聲，得到30 組響應數據，并以此作為實際結構的測量值。分別采用DRAM算法和MH 算法進行貝葉斯模型修正，設置迭代終止次數T=8000，非自適應長度n0=200，建議分布函數個數q=2。同時結合已建立的高斯過程模型來提高非線性模型參數后驗概率分布的計算效率。隨機選取待修正參數的初始值為[ α ,A,μ, β , γ]ini=[0.2, 0.4, 2.5, 0.8, 0.7]，分別利用DRAM 和MH 算法計算得到參數 γ的Markov 鏈如圖7(a)和圖7(b)所示，非線性Bouc-Wen 模型修正結果列出如表2，從圖7 和表2 可以看出，DRAM 算法通過自適應調節建議分布函數，提高了后驗樣本的遍歷性，與MH 算法計算結果相比，DRAM 計算得到的后驗樣本均值與名義值誤差較小，樣本接受率高，且待修正參數后驗樣本變異系數均小于0.04，因此，利用DRAM 進行參數后驗概率估計可以較為快速、準確地量化模型參數的不確定性。同時，基于高斯過程模型和DRAM 算法的非線性模型修正過程運行時間為42 s，與傳統基于有限元模型計算的蒙特卡洛抽樣法相比計算效率得到了大幅提高。

圖7 參數γ 的馬爾科夫鏈Fig. 7 Markov chain of γ

表2 非線性Bouc-Wen 模型參數修正結果Table 2 Updated results of nonlinear Bouc-Wen model parameters

2.3.2 考慮結構參數不確定性的非線性模型修正

取非線性Bouc-Wen 模型參數為待修正參數，為了考慮結構參數的不確定性，設其為符合正態分布的隨機變量，均值為表1 所示參數初始設計名義值，標準差均為0.1。根據結構非線性模型參數的概率分布隨機生成30 組樣本，并代入有限元模型計算得到30 組結構頂層加速度響應，根據DAMD 方法和Hilbert 變換提取每組結構響應主分量的瞬時特征參數作為“測試值”。以第10 s 處主分量響應的瞬時加速度幅值為例，計算得到30 次測試值如圖8 所示，利用表1 所示參數名義值計算得到的瞬時加速度期望值為1.19 m/s2，從圖8 可以看出模擬得到的瞬時加速度幅值響應在期望值附近波動，滿足貝葉斯推斷中誤差的高斯分布假設。

圖8 模擬的30 次瞬時加速度幅值Fig. 8 30 sets of simulated instantaneous acceleration amplitudes

結合高斯過程替代模型，利用DRAM 算法進行貝葉斯模型修正，設置迭代終止次數T=8000，非自適應長度n0=200，建議分布函數個數q=2，待修正參數初始值設為 [ α ,A,μ, β , γ]ini= [0.2, 0.4,2.5, 0.8, 0.7]。以參數 γ為例，其Markov 鏈和剔除燃燒期后的樣本直方分布分別如圖9(a)和圖9(b)所示。非線性結構模型修正結果如表3 所示。從圖9和表3 可知，待修正參數的Markov 鏈燃燒期較短，樣本接受率相比于考慮測量不確定性的計算結果較高，原因歸結于生成測量樣本的非線性模型參數嚴格服從正態分布，DRAM 算法通過自適應策略快速調整建議分布協方差函數逼近待修正參數的真實分布，加速了后驗樣本的收斂。表3中后驗樣本均值與名義值的誤差較小，因此，采用本文提出的基于貝葉斯理論的非線性模型修正方法可以應用于考慮多種不確定性誤差的非線性結構模型修正問題。

表3 非線性結構模型修正結果Table 3 Updated results of nonlinear model

圖9 參數γ 的后驗樣本Fig. 9 Posterior samples of γ

3 結構動力可靠度分析

3.1 非線性動力響應概率密度演化分析

利用實測結構動力響應對結構模型參數進行識別，并對其動力可靠度進行計算，是結構安全性評估的重要手段。以第2 節中地震動作用下具有隨機非線性參數的剪切框架為例，進行結構動力可靠度分析，其中非線性模型參數采用2.3.1 節中表2 所示的DRAM 算法修正結果。結構受Northridge 波形地震動作用，分別以地震動峰值加速度(PGA)和非線性Bouc-Wen 模型參數為隨機變量，其概率分布統計如表4 所示。

表4 非線性模型隨機參數與激勵隨機參數Table 4 Random parameters of nonlinear model and excitation

記框架結構頂層位移為x(t)，同時考慮非線性結構模型參數與激勵參數的隨機性，根據概率密度演化方法計算得到x(t)在不同時刻的概率密度曲線及累積分布函數，如圖10 所示。圖11 給出了不同時段內結構動力響應的概率密度函數曲面及其對應的等值線圖，從圖10 和圖11 可以看出，結構響應的概率密度函數隨著時間不斷變化且趨于復雜，概率密度曲線具有明顯的演化特征不再服從規則的概率分布函數，因此，在傳統結構動力可靠度分析中將結構響應假設為以正態或極值分布為代表的規則分布，可能會導致結構可靠度計算結果出現偏差。

圖10 不同時刻概率分布Fig. 10 Probability distribution of different time

圖11 不同時刻概率密度演化Fig. 11 Probability density evolution of different time

此外，根據概率密度演化方法計算得到的結構頂層位移均值和標準差及其與蒙特卡洛隨機抽樣的比較見圖12。在蒙特卡洛模擬中，首先根據隨機參數的概率分布抽取10 000 組樣本，然后利用抽取的樣本參數進行確定性結構響應計算，并統計結構位移響應的均值和標準差。圖12 對比結果表明，基于概率守恒原理的廣義概率密度演化方法與隨機蒙特卡洛模擬結果相比誤差較小，因而能夠較為真實、準確地反映實際結構動力響應概率分布的演化過程。

圖12 結構頂層位移時程均值和標準差對比Fig. 12 Comparison of time history mean and standard deviation of top floor displacement

3.2 結構動力可靠度計算

為了考慮參數不確定性對結構動力可靠度的影響，分別建立了兩種非線性振動模型：1)僅考慮激勵隨機性的確定性非線性模型，結構模型參數為表1 所示的名義值，記為名義模型；2)同時考慮激勵和結構模型參數的不確定性，按照表4隨機參數概率分布建立的復合隨機模型。

根據首次超越破壞準則，利用廣義概率密度演化方法，分別對不同界限條件下的結構動力可 靠 度R(t) 進 行計算，其 中R(t)=P{|x(τ)|≤xlim,τ∈[0,t]}，t=30 s。為了驗證計算結果的準確性，利用蒙特卡洛模擬方法對復合隨機模型的動力可靠度進行隨機抽樣計算。

分別利用頂層位移和速度作為邊界條件指標，計算非線性結構的動力可靠度。不同位移界限條件下結構動力可靠度計算結果如圖13 和表5所示；不同速度界限條件下的計算結果見圖14 和表6。由表5 和表6 可知，基于概率密度演化的結構動力可靠度計算結果與蒙特卡洛抽樣模擬相比誤差較小，具有較高的計算精度，且名義模型在不同類型界限條件下的動力可靠度均大于復合隨機模型，非線性模型參數的不確定性降低了結構的動力可靠度。此外，當名義模型分別在位移和速度界限條件下動力可靠度達到1.000 時，復合隨機模型的動力可靠度分別為0.9808 和0.9475，由此可知，以速度作為界限條件相比于位移使得復合隨機模型的可靠度計算結果偏低，不同類型界限指標對結構動力可靠度的計算結果會產生一定的影響，因此，在結構安全評估中應充分考慮非線性模型參數的不確定性，使得結構動力可靠度計算結果更加合理、準確。

圖13 不同位移界限條件下結構動力可靠度Fig. 13 Dynamic reliability under different displacement boundary conditions

表5 不同位移界限條件下結構動力可靠度對比Table 5 Dynamic reliability for different displacement boundary conditions

圖14 不同速度界限條件下結構動力可靠度Fig. 14 Dynamic reliability under different velocity boundary conditions

表6 不同速度界限條件下結構動力可靠度對比Table 6 Comparison of dynamic reliability for different velocity boundary conditions

4 結論

本文基于實測結構動力響應主分量的瞬時特征參數構建貝葉斯推理的似然函數，并結合DRAM 算法和高斯過程替代模型，實現了結構的非線性模型修正和參數的不確定性量化。同時，根據非線性模型修正結果，采用廣義概率密度演化方法對實測結構進行安全可靠度評估。通過地震激勵下非線性框架結構的數值模擬結果，得到以下結論：

(1)結構動力響應瞬時特征參數包含了振動系統的幅值和相位信息，能夠較為全面地反映時變非線性結構的物理特性。利用有限個局部峰值點的瞬時特征參數建立似然函數，能夠準確識別結構的非線性模型參數并對其不確定性進行量化。

(2)通過將DRAM 算法和高斯過程替代模型相結合，克服了傳統MCMC 方法樣本接受率低、燃燒期長等缺點，大大提高了待修正參數后驗樣本的計算效率，簡化了非線性結構模型修正的過程。

(3)根據首次超越破壞準則，利用廣義概率密度演化方法能夠準確計算結構的動力可靠度。與確定性模型相比，考慮非線性模型參數不確定性的復合隨機模型能夠更加真實地反映結構的安全狀態，其動力可靠度計算結果總體偏低，因此，在結構安全狀態評估中，考慮模型參數的不確定性能夠更好地反映實際結構的抗震性能水平。