基于反應數據的貝葉斯概率方法在恢復力模型選擇中的應用

劉佩，袁泉，魏慶朝

(北京交通大學 土木建筑工程學院，北京，100044)

模型選擇的標準是考察包含識別參數的模型對待識別結構進行描述的精確程度。數學模型是建立在對原始數據規律性挖掘的基礎上的，但是，不同的建模方法對于原信息的規律把握不同，目前沒有一個絕對優秀的建模方法，只能根據原始數據的特點選擇相對合適的模型。對于模型選擇在系統識別中的應用，需要根據具體識別問題的特點嘗試得到，但由于對待識別結構系統的認識不足以及計算效率的差異，通常會產生識別參數的不當選取問題，或者得到多個不同的模型，導致模型對待識別結構的描述總是存在不完善之處。為了從中找出合理的模型，有必要對模型選擇方法進行研究，來判斷模型及其參數選擇的優劣。過度復雜的模型可能使估計或者預測的方差偏大，而過度簡單的模型又可能導致較大的估計或預測偏差，對于統計學中傳統的模型選擇方法，如AIC 和BIC等[1-2]，這些模型選擇方法忽視了模型選擇過程所帶來的不確定性。傳統的數據分析方法分成模型選擇和基于選定的模型進行統計推斷這2 個步驟。在進行第2步推斷時，人們通常把第1 步選定的模型當作真實的數據產生過程而忽略了模型選擇過程中的不確定性，導致低估實際的方差或均方誤差。在系統識別中，通常采用的是在一系列模型中找到最優模型，但是，因為系統識別是反演問題，所以，很可能存在多個模型都可以對結構特性進行很好的模擬，因此，需要發展更具一般性的模型選擇方法。張立濤等[3]結合結構識別問題的求解可采用正則化技術來改善其不適定的特點，引入不同識別結果所對應目標函數與識別結果的范數之間的對應規律作為先驗條件，提出了一種基于L 曲線的模型確認方法。Goulet 等[4]提出了一種多模型方法來考慮不確定性及模型假定，該方法首先生成大量的候選模型，對反應進行預測，其次對結構進行靜力試驗，將模型不確定性及測量不確定性之和與設定的界限值對比，來判斷接受或拒絕該模型。考慮到模型選擇過程中的不確定性，基于貝葉斯理論的模型選擇方法[5]可以同時考慮多個模型對結構反應預測的影響，還可以自動對復雜模型加以限制，定量計算得到選擇各模型的概率。貝葉斯統計推斷一般包括參數識別和模型選擇：基于觀測數據，通過貝葉斯理論或者結合一些抽樣方法進行模型參數識別可以得到更新的模型參數概率密度函數；而通過貝葉斯理論進行模型選擇可以評估各模型用來估計特定觀測數據的概率能力[6-11]。本文作者通過推導用于計算各模型選擇概率的證據表達式，建立基于貝葉斯理論的模型選擇方法計算框架；當已知數據點足夠多且為全局可識別情況時，可通過拉普拉斯漸近估計解法利用模型參數識別結果對證據進行計算。隨后將基于貝葉斯理論的模型選擇方法用于恢復力模型的選擇中，首先是對雙線性單自由度體系在地震波作用下的恢復力模型的選擇，重點將基于貝葉斯理論的模型選擇方法用于實測數據的模型選擇，即采用實測滯回曲線數據的密肋復合墻試件的恢復力模型的選擇，可利用各模型參數識別結果計算選擇各模型的概率。

1 基于貝葉斯理論的模型選擇方法

1.1 計算框架

令y 表示模擬的或實測的結構反應數據，模型選擇 的 目 的 是 利 用 y 從 一 系 列 模 型μ={μj:j=1,2, …,NM}中選擇最有可能的模型。根據貝葉斯定理，可以得到基于反應數據y 選擇各模型的概率為

根據全概率公式，已知 y 時模型μj的證據p(y| μj)可表示為

其中：j=1, 2, …, NM； θj表示模型 μj的參數矢量；p(θj| μj)為模型參數的先驗概率密度函數，可根據經驗確定； p(y | θj, μj)為模型參數的似然函數。

對于全局可識別情況，當數據y 足夠多時， θj的后驗概率密度函數可通過高斯分布進行估計[5]，則根據拉普拉斯漸近估計解法 p(y| μj)可近似表示為

綜上，根據貝葉斯模型選擇方法，可根據p(y | μj) P(μj| μ)對各模型進行排序，最有可能的模型為使 p(y | μj) P(μj| μ)取得最大值的模型。模型證據 p(y| μj)可通過式(3)計算得到。模型先驗分布P (μj| μ)可通過均勻分布來進行描述。

1.2 模型參數最有可能值的確定

其中，k=1/ p(y)為正規化常數。則使后驗分布取最大值，相當于使 -ln[ p(y | θ) p(θ)]取最小值。

假定模型計算反應為q，則

其中，e 為模型誤差，假定e 在不同的測試點處為獨立同分布的均值為0，標準差為σ 的高斯分布，則

則

對式(7)的σ 求一階導數并令其等于0，假定模型參數的先驗分布為均勻分布，可得模型誤差方差的最有可能值為

根據式(8)并通過Matlab 中的fminsearch函數對模型參數進行優化，可得到模型參數的最有可能值，進而可計算得到最有可能值處 -ln[ p(y | θ) p(θ)]的Hessian 矩陣。

2 所提方法在恢復力模型選擇中的應用

本文將基于貝葉斯理論的模型選擇方法用于恢復力模型的選擇：(1) 受地震作用的雙線性單自由度體系恢復力模型采用模擬輸入輸出數據；(2) 低周反復加載下密肋復合墻試件恢復力模型采用實測滯回曲線數據。

2.1 地震作用下的單自由度雙線性體系

考慮單自由度雙線性滯回體系的運動方程：

其中：x 為位移；m 為質量；c 為阻尼系數且假定為線性黏滯阻尼；k1為屈服前剛度；k2為屈服后剛度；xy為屈服位移；fh(x; k1, k2, xy)為恢復力，其特性見圖1。

圖1 單自由度體系雙線性恢復力模型Fig.1 Bilinear restoring forcing model of SDOF system

假定已知反應數據為結構在地震作用下所得位移加上其5%均方差的高斯白噪聲。考慮2 個恢復力模型，各模型均采用零均值的離散高斯白噪聲作為預測誤差模型。對模型1(μs,1)：雙線性滯回體系，屈服前剛度k1＞0，屈服后剛度k2＞0，黏滯阻尼c=0，屈服位移為xy，預測誤差標準差為 ση。對模型2(μs,2)：雙線性滯回體系，屈服前剛度k1＞0，屈服后剛度k2＞0，黏滯阻尼c＞0，屈服位移為xy，預測誤差標準差為 ση，該模型包含產生模擬反應數據的精確模型。

圖2 地震波及其作用下體系的反應Fig.2 Earthquake record and corresponding displacement of system

圖3 地震波作用下體系的滯回環Fig.3 Hysteretic curve of system under earthquake record

假定各模型參數的先驗分布相互獨立且為均勻分布，k1的分布范圍為(0, 2) N/m，k2的分布范圍為(0, 0.5)N/m，c 的分布范圍為(0, 0.1) N·s/m，xy的分布范圍為(0, 0.1) m， ση的分布范圍為(0, 0.01) m。

表1 所示為各模型在地震作用下根據式(8)識別得到的最有可能值。

地震作用下基于反應數據選擇各模型的概率P (μj| y, μ)，j=1, 2 可根據式(1)計算得到，其中各模型的證據可根據式(3)計算得到，各模型的先驗概率取P (μj| μ)=1/2，最終計算結果為P (μs,1| y, μ) ≈0，P (μs,2| y, μ)=1。結果表明，基于反應數據選擇最有可能模型的概率接近于1，選擇另外的模型概率近似為0，表明對該體系反應進行預測時，除最有可能模型外的模型都可以忽略。另外，該算例表明進行系統識別時，對于一個實際結構，最優的模型依賴于已知的系統數據，應該只對與系統識別中系統所受的激勵水平接近的激勵作用下的系統反應進行預測。

表1 單自由度體系各雙線性恢復力模型參數的最有可能值Table 1 Most probable values of model parameters for each bilinear restoring force model of SDOF system

2.2 低周反復加載下的密肋復合墻試件

將基于貝葉斯理論的模型選擇方法用于已知實測數據的恢復力模型選擇中。作為該模型選擇方法的研究基礎，劉佩等[13]進行了基于貝葉斯理論的密肋復合墻試件恢復力模型參數的識別。

文獻[13]中給出了7 塊密肋復合墻試件在低周反復加載下發生剪切破壞的滯回曲線，標準密肋復合墻試件的實測滯回曲線見圖4。從圖4 可知：試件出現了明顯的剛度退化、強度退化、滑移和捏攏現象。

圖4 標準密肋復合墻試件實測滯回曲線Fig.4 Tested hysteretic curve of typical multi-gird composite wall specimen

根據以上現象，文獻[13]提出了一種適用于剪切破壞的密肋復合墻體的恢復力模型，本文記為μm,1，見圖5。模型達到極限荷載之前骨架曲線假定為雙線性。模型的控制參數為：屈服前剛度k1，屈服后剛度k2，屈服荷載fy及其對應位移xy，極限荷載fu及其對應位移xu，考慮極限荷載后強度降低現象的參數Δd，考慮滑移捏攏現象的參數Δp 和Δs。模型的滯回規律假定為：(1) 達到極限荷載前，前一次循環結束之后再加載時和反向加載時，直線指向前一次循環的最大變形點；(2) 達到極限荷載后，前一次循環結束之后再加載時和反向加載時，直線指向由前一次循環的最大變形點與Δd 之和對應的位移及極限荷載確定的點的位置處；(3) 正向卸載及反向卸載直線通過考慮滑移捏攏參數軸正負Δs 位置前剛度為k1，之后指向力軸的考慮滑移捏攏參數負正Δp 位置處。

圖5 密肋復合墻體的恢復力模型μ1Fig.5 Restoring force model μ1 of multi-grid composite walls

對于模型 μ1，若將極限荷載后強度降低現象通過將曲線指向上一循環的最大荷載降低一定的值與最大變形確定的位置表示，同樣可以作為用來模擬密肋復合墻體反應的恢復力模型，即將考慮極限荷載后強度降低現象的參數Δd 變為ΔF[14]，本文將該恢復力模型記為μm,2，與 μ1的滯回規律假定(2)對比，模型 μ2將其變為達到極限荷載后，前一次循環結束之后再加載時和反向加載時，直線指向由前一次循環的最大變形點對應的位移及荷載與ΔF 之差確定的點的位置處。

這2 種模型都可以用來模擬構件的反應，但要對這2 種模型識別結果的差別定量化，則需要應用基于貝葉斯理論的模型選擇方法確定，以便決定模擬構件反應時每種模型的權重。

采用標準密肋復合墻試件的實測滯回曲線作為已知數據，對選擇2 個模型μm,1和μm,2的概率進行計算。

根據統計學理論中先驗分布[15-16]的特點，其主要來源于經驗和歷史資料。本算例采用無信息先驗分布，并根據試驗數據確定模型參數先驗分布的取值區間。設各模型參數的先驗分布相互獨立且為均勻分布：k1的分布范圍為(5, 15) kN/mm，k2的分布范圍為(0, 5)kN/mm，fy的分布范圍為(45, 75) kN，Δp 的分布范圍為(1, 7) kN，Δs 的分布范圍為(5, 35) kN，Δd 的分布范圍為(2, 6) mm，ΔF 的分布范圍為(3, 9) kN，σ 的分布范圍為(1, 3) mm。經計算發現，模型參數的先驗分布對模型參數的最有可能值及Hessian 矩陣的結果沒有影響。

根據貝葉斯理論，由式(8)識別得到的各模型參數的最有可能值見表2。

設模型先驗概率P(μm,1| μ)=P(μm,2| μ)=1/2，則依據已知數據并根據式(1)得選擇μm,1的概率為

選擇μm,2的概率為

由式(10)和(11)可得計算結果為P (μm,1| y, μ)=0.697 7，P (μm,2| y, μ)=0.302 3?？梢姡簩υ嚰磻M行預測時，2 個模型均應考慮，可將依據已知數據選擇各模型的概率作為權重系數，來考慮各模型在反應預測時所占的比例。

表2 密肋復合墻體各恢復力模型參數的最有可能值Table 2 Most probable values of model parameters for each restoring force model of multi-grid composite wall

3 結論

(1) 將基于貝葉斯理論的模型選擇方法用于采用實測數據的恢復力模型的選擇中，其中通過貝葉斯理論對恢復力模型參數進行了識別。

(2) 通過推導各模型的證據表達式，并結合貝葉斯定理得到選擇各模型的概率，建立了基于貝葉斯理論的模型選擇方法計算框架。該方法不僅可以識別得到模型參數的最有可能值，還可以識別得到最有可能的模型，并且考慮了模型及模型參數的不確定性。

(3) 根據雙線性滯回體系在地震作用下響應的模擬數據對恢復力模型進行選擇，若選擇最有可能模型的概率接近于1，則進行反應預測時可以忽略其余模型的影響。

(4) 根據密肋復合墻試件在低周反復荷載下的實測滯回曲線數據對恢復力模型進行選擇，得到了各模型的定量概率，結果表明：基于反應數據選擇各模型的概率相差不大，此時需同時考慮各模型來進行反應預測，并可將選擇各模型的概率作為反應預測的權重系數。

[1] Beck J L. Bayesian system identification based on probability logic[J]. Structural Control and Health Monitoring, 2010, 17(7):825-847.

[2] 張新雨, 鄒國華. 模型平均方法及其在預測中的應用[J]. 統計研究, 2011, 28(6): 97-102.ZHANG Xinyu, ZOU Guohua. Model averaging method and its application in forecast[J]. Statistical Research, 2011, 28(6):97-102.

[3] 張立濤, 李兆霞, 張宇峰, 等. 結構識別計算中基于L 曲線的模型確認方法研究[J]. 振動與沖擊, 2011, 30(11): 36-41.ZHANG Litao, LI Zhaoxia, ZHANG Yufeng, et al. Methodology of L-curve based model validation in structural identification[J].Journal of Vibration and Shock, 2011, 30(11): 36-41.

[4] Goulet J A, Kripakaran P, Smith F C. Multimodel structural performance monitoring[J]. ASCE Journal of Structural Engineering, 2010, 136(10): 1309-1318.

[5] Beck J L, Yuen K V. Model selection using response measurements: Bayesian probabilistic approach[J]. ASCE Journal of Engineering Mechanics, 2004, 130(2): 192-203.

[6] Mthembu L, Marwala T, Friswell M I, et al. Model selection in finite element model updating using the Bayesian evidence statistic[J]. Mechanical Systems and Signal Processing, 2011,25(7): 2399-2412.

[7] Cheung S H, Beck J L. Calculation of posterior probabilities for Bayesian model class assessment and averaging from posterior samples based on dynamic system data[J]. Computer-Aided Civil and Infrastructure Engineering, 2010, 25(5): 304-321.

[8] Ching J, Chen Y C. Transitional Markov chain Monte Carlo method for Bayesian model updating, model class selection, and model averaging[J]. ASCE Journal of Engineering Mechanics,2007, 133(7): 816-832.

[9] Cheung S H, Beck J L. Bayesian model updating using hybrid Monte Carlo simulation with application to structural dynamic models with many uncertain parameters[J]. ASCE Journal of Engineering Mechanics, 2009, 135(4): 243-255.

[10] Muto M, Beck J L. Bayesian updating and model class selection for hysteretic structural models using stochastic simulation[J].Journal of Vibration and Control, 2008, 14(1/2): 7-34.

[11] Beck J L, Au S K. Bayesian updating of structural models and reliability using Markov Chain Monte Carlo simulation[J].Journal of Engineering Mechanics, 2000, 128(4): 380-391.

[12] Beck J L, Katafygiotis L S. Updating models and their uncertainties. Ⅰ: Bayesian statistical framework[J]. Journal of Engineering Mechanics, 1998, 124(4): 455-461.

[13] 劉佩, 袁泉, 魏慶朝. 概率方法在恢復力模型參數識別中的應用分析[J]. 中南大學學報(自然科學版), 2013, 44(9):3843-3848.LIU Pei, YUAN Quan, WEI Qingchao. Application of a probabilistic approach to restoring force model parameter identification[J]. Journal of Central South University (Science and Technology), 2013, 44(9): 3843-3848.

[14] Qu Z, Ye L P. Strength deterioration models based on effective hysteretic energy dissipation for RC members under cyclic loading[C]// 7th International Conference on Urban Earthquake Engineering & 5th International Conference on Earthquake Engineering. Tokyo: Tokyo Institute of Technology, 2010:851-856.

[15] Siu N O, Kelly D L. Bayesian parameter estimation in probabilistic risk assessment[J]. Reliability Engineering and System Safety, 1998, 62(1/2): 89-116.

[16] 茆師松, 王靜龍, 濮曉龍. 高等數理統計[M]. 2 版. 北京: 高等教育出版社, 2006: 368-372.MAO Shisong, WANG Jinglong, PU Xiaolong. Advanced mathematical statistics[M]. 2nd ed. Beijing: Higher Education Press, 2006: 368-372.