劉佩，袁泉，魏慶朝

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

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

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

1.1 計算框架

令y 表示模擬的或?qū)崪y的結(jié)構(gòu)反應(yīng)數(shù)據(jù)，模型選擇 的 目 的 是 利 用 y 從 一 系 列 模 型μ={μj:j=1,2, …,NM}中選擇最有可能的模型。根據(jù)貝葉斯定理，可以得到基于反應(yīng)數(shù)據(jù)y 選擇各模型的概率為

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

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

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

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

1.2 模型參數(shù)最有可能值的確定

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

假定模型計算反應(yīng)為q，則

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

則

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

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

2 所提方法在恢復(fù)力模型選擇中的應(yīng)用

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

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

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

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

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

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

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

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

假定各模型參數(shù)的先驗分布相互獨立且為均勻分布，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 所示為各模型在地震作用下根據(jù)式(8)識別得到的最有可能值。

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

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

2.2 低周反復(fù)加載下的密肋復(fù)合墻試件

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

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

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

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

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

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

這2 種模型都可以用來模擬構(gòu)件的反應(yīng)，但要對這2 種模型識別結(jié)果的差別定量化，則需要應(yīng)用基于貝葉斯理論的模型選擇方法確定，以便決定模擬構(gòu)件反應(yīng)時每種模型的權(quán)重。

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

根據(jù)統(tǒng)計學(xué)理論中先驗分布[15-16]的特點，其主要來源于經(jīng)驗和歷史資料。本算例采用無信息先驗分布，并根據(jù)試驗數(shù)據(jù)確定模型參數(shù)先驗分布的取值區(qū)間。設(shè)各模型參數(shù)的先驗分布相互獨立且為均勻分布：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。經(jīng)計算發(fā)現(xiàn)，模型參數(shù)的先驗分布對模型參數(shù)的最有可能值及Hessian 矩陣的結(jié)果沒有影響。

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

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

選擇μm,2的概率為

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

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

3 結(jié)論

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

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

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

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

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

[2] 張新雨, 鄒國華. 模型平均方法及其在預(yù)測中的應(yīng)用[J]. 統(tǒng)計研究, 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] 張立濤, 李兆霞, 張宇峰, 等. 結(jié)構(gòu)識別計算中基于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] 劉佩, 袁泉, 魏慶朝. 概率方法在恢復(fù)力模型參數(shù)識別中的應(yīng)用分析[J]. 中南大學(xué)學(xué)報(自然科學(xué)版), 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] 茆師松, 王靜龍, 濮曉龍. 高等數(shù)理統(tǒng)計[M]. 2 版. 北京: 高等教育出版社, 2006: 368-372.MAO Shisong, WANG Jinglong, PU Xiaolong. Advanced mathematical statistics[M]. 2nd ed. Beijing: Higher Education Press, 2006: 368-372.