基于Bayesian證據推斷與信息增益的參數化有限元修正模型選擇

尹 濤, 王祥宇, 周 越

(武漢大學 土木建筑工程學院，武漢 430072)

經過幾十年的發展，有限元分析方法在工程中的應用日益廣泛，但由于各種理論假設、邊界條件近似性以及幾何材料特性參數等不確定性因素影響，致使結構有限元分析模型計算響應與實測結構響應之間不可避免地存在一定偏差。基于試驗模態分析方法的有限元模型修正技術正是提高有限元模型精度的有效手段，該方法通過實測低階模態參數數據進行結構有限元模型修正，從而獲得與實際結構較好吻合的修正后有限元模型，并以此為基準模型開展結構響應預測、振動控制、性能評估以及健康監測等領域研究，相關問題已成為目前國內外的研究熱點[1-5]。

為求更好地模擬實際結構，開展結構正向分析，人們往往傾向于建立精細的有限元模型，以求細致描述結構各構件局部及結構整體的物理力學特征，這將顯著增加結構有限元模型的復雜程度，直接導致過于龐大的有限元模型規模。事實上，結構有限元模型修正過程可視為帶約束的參數優化過程，其在優化算法框架內，通過迭代計算不斷調模型參數，使有限元模型的計算響應輸出與實測數據之間的誤差達到極小，屬于典型的結構動力學反問題范疇。但在實際應用中，過于復雜的有限元模型會造成至少兩方面不利影響。一方面，過于龐大的模型規模非常不利于動力學反問題優化分析過程中的反復迭代計算；另一方面，受傳感器通道數量、測量噪聲以及高階模態獲取困難等因素制約，實際獲得用于有限元模型修正的有效測量信息相對于模型待修正參數往往不足，使得模型參數修正過程容易出現病態，同時也會顯著增加待識別模型參數相對于模型誤差及測量噪聲等不確定性因素的敏感性，容易出現所謂‘過擬合’問題[6]，而這對含有大量待定模型參數的復雜有限元模型尤其如此。因此，在保持有限元模型計算響應與結構實測響應吻合程度條件下，如何有效降低有限元分析模型的參數化復雜程度，值得進一步研究。

目前常采用敏感性分析方法來進行有限元模型參數選擇，以減少待修正模型參數數量，即通常選擇對結構響應敏感性較大的模型參數作為待修正參數[7-12]。但對于規模較大的結構模型，該敏感性分析方法所選取的模型參數并不一定能很好地用作待修正參數[13]。針對有限元模型修正中的待定參數選擇問題，本文以概率論和信息論為基礎，發展一種基于Bayesian證據推斷與馬爾科夫鏈蒙特卡洛(Markov Chain Monte Carlo，MCMC)法相結合的有限元修正模型選擇分析方法。通過引入信息增益(Information Divergence)或KL散度(Kullback-Leibler Divergence)[14]指標(該指標在本文具體指從實測數據中提取用于未知模型參數修正的信息量多少)，對有限元模型待修正參數的復雜程度進行懲罰，從而在有限元模型建模參數復雜性與其相應信息論復雜性之間獲得平衡，最終獲得滿足模型與實測數據吻合度的相對簡單有限元修正模型，有效避免由于待修正參數過多而導致的模型過擬合問題。通過某兩層螺栓連接框架結構有限元模型修正的數值仿真與實驗室模型試驗研究，對本文方法正確性與有效性進行驗證。

1 基本理論

1.1 Bayesian多模型證據推斷

考慮目標結構存在NC種可能的有限元參數化建模方式(或稱之為模型類別)，分別定義為C1, C2, …, CNC。本文假定各有限元模型類別之間通過不同的建模參數進行區分，即第j個模型類別Cj(θj)通過其待估參數向量θj來描述，θj=θj(Θj)為模型參數空間Θj?RNj中的向量，Nj表示向量θj所含模型參數的數目。基于Bayes定理，在給定測量數據D的條件下，模型類別Cj的后驗概率密度函數可以表示為

(1)

式中：U為使用者對于各模型類別Cj先驗概率密度函數p(Cj|U)(j=1,2,…,NC)的經驗判斷；其滿足

(2)

式(1)中分子項p(D|U)可由全概率公式給出

(3)

式中：D為從目標結構獲得的測量數據集合，本文具體指Ns組重復測量得到的前Nm階固有頻率與模態振型，即

(4)

此外，式(1)中p(D|Cj,U)表示模型類別Cj在給定數據集D條件下的證據，忽略用戶經驗判斷U，并結合此前模型類別定義，該證據項可以通過全概率積分形式表示為

(5)

式中：p(θj|Cj)為對于指定模型類別Cj，定義在其模型參數空間Θj中參數向量θj的先驗概率密度函數；對于本文考慮的有限元模型修正問題，假定該先驗概率服從高斯分布。p(D|θj,Cj)為給定數據集D和模型類別Cj條件下的似然函數，即

(6)

式中：J(θj;D,Cj)為有限元模型類別Cj在給定模型參數θj條件下的固有頻率、模態振型計算值與Ns組測量值之間的誤差函數，且假定該誤差服從均值0、方差為σ2的獨立高斯分布，即

(7)

依據Bayes定理，在給定測量數據集D與模型類Cj條件下，模型參數θj的后驗概率密度函數可以表示為

(8)

式中：等式右邊分母項p(D|Cj)即為式(5)給出的證據因子。

為方便，取證據因子的對數形式進行計算，即對式(5)兩邊同時取對數，可得

(9)

考慮式(8)所給出的后驗概率分布，式(9)等號右邊可以進一步表示為

(10)

式中：等號右邊兩個數學期望式分別為

(11)

與

(12)

式(11)反映有限元模型類別Cj的計算數據與實際測量數據的平均吻合程度，即該值越大，表示有限元模型與實測數據越符合。式(12)稱為KL散度或信息增益，其表示對含模型參數θj的有限元模型類Cj進行修正時所需要從測量數據集D中提取的信息量多少，該值恒為非負，且其取較大值時表明對Cj類模型進行修正時需要提取的信息量多。事實上，有限元模型修正過程為特征值反問題求解過程，與函數擬合過程類似，并非模型越復雜越好，過于復雜的模型有時會造成對實測數據的‘過擬合’問題，過度擬合了實測數據中的噪聲，而式(12)中給出的信息增益量度能夠較合理地對模型復雜程度進行懲罰，這在后面的數值仿真與模型實驗部分中將獲得驗證。

1.2 基于后驗分布中MH抽樣的數學期望值估計

Metropolis Hasting(MH)是蒙特卡羅馬爾科夫鏈(MCMC)方法中一類重要的抽樣方法[15]。本文采用該方法分別計算式(11)與式(12)在后驗分布p(θj|D,Cj)中定義的多重數值積分，獲得以數學期望值表示的模型數據吻合度與信息增益，從而通過式(10)得到模型類別Cj以對數形式表達的證據因子。

(13)

(14)

(15)

(16)

此時，式(9)中對數形式的證據因子ln[p(D|Cj) ]可以通過對以上兩個數學期望值求差獲得，因而，式(1)中給出的有限元模型類別Cj的后驗分布則可以通過下式計算為[16]

(17)

式中，M為ln[p(D|Cj) ]在所有模型類別Cj(j=1,2,…,NC)中的最大值。模型后驗概率p(Cj|D)從概率論角度對有限元修正模型的合理選擇進行解釋，即某類有限元模型Cj的后驗概率值越大，表明該類候選模型選作有限元修正模型的可能性越高。

2 數值仿真研究

2.1 框架模型介紹

為驗證本文方法，對某兩層平面螺栓鋼框架進行數值仿真研究，如圖1(a)所示。其中，梁-柱與柱-基礎均采用螺栓連接，模型共4個梁-柱連接部位及2個柱-基礎連接部位。該框架梁柱均采用10#工字鋼，其幾何與材料特性參數分別為：楊氏模量E=2.01×1011N/m2，質量密度ρ=7.85×103kg/m3，截面面積A=1.57×10-3m2及截面慣性矩I=2.21×10-6m4。該兩層框架模型中梁-柱及柱-基礎節點均視為半剛性連接，其轉動剛度采用端部短梁的抗彎剛度來模擬，并通過調整端部短梁的抗彎剛度來模擬螺栓連接松緊程度。

此框架有限元離散模型示意圖及相應測量節點與方向編號，如圖1(b)所示。共有30個單元、30個節點以及84個自由度，6個半剛性節點，且在其中8個自由度方向上布置模態拾振點。取該框架結構前4階模態頻率與振型進仿真計算，對計算頻率與振型分別施加標

圖1 兩層螺栓框架的模型圖Fig.1 Schematic model of two-storey bolt-connected frame

準差為1%與10%的高斯白噪聲以模擬測量噪聲情況，且假定通過重復測量共獲得30組模態參數數據，即Nm=4，Ns=30。

本算例將此框架結構中所有螺栓半剛性連接的轉動剛度k1,k2,…,k6均視為待修正參數，如圖2所示。NC=6類有限元模型。從C1～C6，模型所含建模參數逐漸增多，即模型復雜程度相應逐漸增加。具體地，C1為相對最簡單模型，其表示所有半剛性連接均采用相同的單一模型參數θ1進行修正。C2含兩個建模參數θ1與θ2，其分別修正所有柱-基礎連接剛度與梁-柱連接剛度。C3通過3個建模參數θ1，θ2與θ3分別修正所有柱-基礎連接剛度與各層梁-柱連接剛度。在C3基礎上，C4進一步區分左、右柱的柱-基礎連接剛度，繼續增加了模型的復雜程度。C5在C3基礎上，進一步將各層左、右柱的梁-柱連接剛度區別開來，共含有θ1～θ5此5個建模參數。C6在C5基礎上進一步對左、右柱-基礎連接剛度k1和k2進行區分，共包含6個模型參數，且該類模型為本算例所考慮的相對最復雜有限元參數化模型類別。

圖2 所考慮的各有限元模型類Fig.2 Different classesof FE models considered

此外，本模型算例共考慮5種仿真工況，各工況具體設置情況，如表1所示。其中kb與kbc分別表示柱-基礎與梁-柱半剛性連接剛度值，本例分別取為0.03EI與0.06EI。工況1為基準狀態情況，即所有梁-柱及柱-基礎半剛性連接剛度均分別取為各自基準值；工況2基于基準狀態，將首層左側梁-柱連接剛度值k3降低為其基準狀態值的0.7倍；工況3在基準工況下，將左側柱-基礎連接剛度k1取為其基準狀態值的0.7倍；工況4進一步考慮兩處梁-柱連接松動情況，即在工況2基礎上，將第二層右側梁-柱連接剛度k6值進行降低；工況5將工況2與工況3中的半剛性節點松動情況進行組合，即研究左側柱-基礎連接剛度k1與第一層左側梁-柱連接剛度k3同時發生降低的情況。

表1 數值仿真工況設置

MCMC算法模擬次數設為10 000，且先驗分布中各模型參數均服從以基準狀態值為均值、標準差為0.3的正態分布，使模型參數空間定義在給定的(0,2]范圍內。以工況2為例，圖3為此工況下模型類C2中模型參數θ1與θ2的MCMC采樣歷程。從圖3可知，在采樣歷程中，模型參數θ1與θ2的值圍繞各自均值在較小范圍內波動。同時，梁-柱連接剛度的整體模型參數θ2值略小于1，這與工況2中所定義的框架第一層左側梁-柱連接剛度值k3降低的事實相符(參見表1與圖2)。此外，模型參數θ1偏離均值的離散度較θ2大，表明該兩層螺栓框架結構中柱-基礎連接剛度的不確定性程度較梁-柱連接剛度大，其與文獻[17]中的有關研究結果規律相符(見圖3)。

圖3 仿真工況2下模型類C2中建模參數的MCMC采樣歷程Fig.3 MCMC sampling history of modeling parameters in C2for case 2

圖4給出工況2下，模型類C2中模型參數θ1與θ2的聯合先驗分布與后驗分布情況對比。從圖4可知，后驗分布范圍僅占先驗分布的一部分，經過Bayesian更新，模型參數θ1與θ2的不確定性程度明顯降低，或者說在有限元模型修正過程中，算法從先驗分布中提取了部分有用信息用于對于模型參數θ1與θ2的更新，依式(12)之定義，其反映從數據D中提取的關于模型類C2的信息增益大小。此外，模型參數θ1的不確定性較θ2大，這與圖3中的結果相符(見圖4)。

圖4 工況2中模型類C2中建模參數的聯合先驗與后驗分布Fig.4 Jointed prior and posterior distributions of modeling parameters in C2 for case 2

采用本文方法對表1中所列各仿真工況進行有限元修正模型選擇研究，分別計算各工況下各類有限元模型的數據吻合度(式(11))、信息增益(式(12))、對數形式的證據因子ln[p(D|Cj) ](式(10))以及后驗概率p(Cj|D)(式(17))，計算結果列于表2，其中，最優模型相關數據用粗體表示。對于工況1，本文方法給出的最優有限元模型類別為C1，其為所有考慮模型類別中最簡單的一類。這表明在基準狀態情況下，對本例僅單個參數就足夠修正該有限元模型，而并非模型越復雜越好。

表2 各數值工況下模型選擇結果

在工況2下，由于其在基準狀態基礎上對兩層框架第一層左側梁-柱連接剛度值k3進行了局部調整，相比于工況1，工況2理論上需要較多的模型參數才能更好地對該局部剛度改變進行刻畫，因而本文方法確定C3為最優模型類的結論符合上述推斷。同時，從表2給出的工況2計算結果中還可知，隨著模型復雜程度增加，即模型參數逐漸增多，有限元模型修正過程就需要從數據D中提取更多的信息用于修正這些額外增加的模型參數，表現在信息增益指標變大。該信息增益指標又反過來對模型的復雜程度起到一定的懲罰作用。例如，工況2中C4與C5兩類模型復雜程度均較C3高，因而相應的信息增益也逐漸增加，表明需要提取更多信息來修正C4與C5中較多的模型參數。但信息增益也并非一定隨著模型復雜程度增加而單調增大，例如，此工況下的C6類模型，其復雜程度相對最高，但對其進行模型修正所需提取的信息量卻較相對簡單的C5類模型少。此外，還應指出，在工況2下，雖然模型類C2較C3簡單，其也將梁-柱連接剛度與柱-基礎連接剛度區分開(見圖2)，但由于其并未進一步區分框架第一層與第二層梁-柱連接剛度，難以細致刻畫僅由k3引起的第一層局部剛度改變，因而，造成模型類C2的數據吻合度相對C3差很多(對比表2中工況2的第2、第3行結果)，因此，盡管修正C2類模型所需提取的信息量較C3類模型相對要少，但綜合來看，模型類C3的證據因子較C2大很多，其導致C3類模型以34.63%的后驗概率p(Cj|D)被選中，而C2類模型被選中的概率僅為0.34%。

表2中所列其余3種工況的計算結果規律與工況1、工況2類似，且對于不同工況，在數據吻合度基本相同條件下，有限元模型修正過程均傾向于選擇較為簡單的參數化模型類別。此外，圖5綜合給出了所有計算工況下各類模型的后驗概率p(Cj|D)，從圖5可知，在綜合考慮信息論表達復雜性條件下，數據吻合程度接近情況下，復雜程度相對較低的參數化模型類別對于有限元模型修正更為適合。

圖5 各仿真工況下 p(Cj |D)的結果對比Fig.5 Comparison of simulation results of p(Cj |D) for all cases

3 實驗驗證

通過實驗室螺栓框架結構模型進一步驗證本文方法，該實驗鋼框架幾何尺寸與材料特性均與此前數值仿真所采用的框架模型對應，同樣由4根10#工字型鋼組成，梁-柱通過角鋼與螺栓進行連接，而柱-基礎則通過柱底端部焊接鋼板與基礎進行螺栓連接，通過扳手控制螺栓松緊程度，本實驗框架整體與連接細部如圖6所示。

模態試驗中傳感器布設位置及方向亦與此前數值仿真保持一致，通過沖擊力錘施加框架平面內的水平向瞬態激勵(見圖6(e))。利用MPS-140801-IEPE信號采集卡將采集到的結構加速度自由響應信號通過USB端口傳輸到PC機，采用MPS-140801 Date Acquisition軟件動態顯示并存儲信號，再通過自編NExT-ERA模態識別程序[18-19]獲得框架結構的前4階固有頻率與模態振型，其中，采樣頻率設為1 000 Hz，采樣時長為10 s。保持激勵位置與方向不變，基準狀態下重復進行多次模態試驗，并獲得Ns=15組模態參數，其中，前4階固有頻率的平均值分別為：17.59 Hz，77.76 Hz，148.85 Hz及159.80 Hz。通過扳手松動框架第一層左側梁-柱連接上部角鋼的一對螺栓，如圖6(d)所示。如圖7所示，并保持此前激勵方式，通過重復模態試驗獲得此狀態下的15組模態參數，且前4階模態頻率均值分別為：16.72 Hz，77.41 Hz，148.74 Hz及157.02 Hz。

圖6 兩層螺栓框架實驗模型Fig.6 Experimental two-storey bolt-connected steel frame

圖7 兩層實驗框架螺栓松動部位圖Fig.7 Bolt-loosening location of the two-storey laboratory frame

基于模態試驗測得的兩種狀態下Ns組模態參數數據，采用本文方法進行實驗框架結構的有限元模型選擇研究，其中MCMC算法的控制參數設置與數值仿真一致。對該實驗框架考慮兩種工況，即基準狀態工況與梁-柱連接部位局部螺栓松動工況，其分別與數值仿真中的前兩種工況相對應，計算結果列于表3。與表2所列項目相同，表3分別給出兩種工況下各類有限元模型的數據吻合度、信息增益、證據因子及后驗概率。對于實驗工況1，通過本文方法確定的最優有限元參數化模型類別為所考慮的最簡單模型C1，其與數值仿真工況1結果相符，表明基準狀態條件下，有限元模型修正傾向于選擇相對最簡單的參數化分析模型。同時，該實驗工況1的結果還可以看出，對相對最復雜模型類別C6(含6個模型參數)的修正所需從實測數據D中提取的信息量，比相對簡單的含3參數的C3類模型更少；因而，在兩者模型數據吻合程度接近的情況下，C6類模型的證據因子(或后驗概率)較C3更大，表明此狀態下C6類模型相對C3較優。

表3 實驗工況下模型選擇結果

實驗工況2考慮第一層左側梁-柱半剛性連接局部螺栓松動情況(見圖7)，其計算結果也表現出與仿真工況2相似的規律。從表3可知，模型類別C3的后驗概率最大，被確定為此工況下的最優模型類別。同時，也可以看出，在此局部剛度變化情況下，有限元模型修正傾向于選擇相對較復雜的參數化模型類別，以較好地反映該局部剛度變化特征。但從信息論復雜度表征角度來看，在數據吻合程度類似條件，過于復雜的參數化模型反而不利于模型修正，因為此時需要從實測數據D中提取的信息量最多，其對模型復雜程度的增加起到懲罰作用。例如，表3工況2中，模型類C3與C4數據吻合程度很相近，但由于C4在修正過程中需要提取更多信息，導致其最終的證據因子或后驗概率比C3小。

此外，為更直觀地比較實驗與數值仿真結果，圖8分別比較了實驗與仿真條件下工況1與工況2中各類模型的后驗概率p(Cj|D)，從圖8可知，兩者吻合程度較好，反映出實驗結果與數值仿真的一致性，進一步驗證了本文方法的合理性與正確性。

圖8 模型后驗概率p(Cj |D)的數值仿真與試驗結果對比Fig. 8 Comparison between simulation and experiment results of model posterior PDF p(Cj |D)

4 結 論

本文基于Bayesian證據推斷理論與信息增益指標，并結合MCMC中的Metropolis Hastings抽樣法，對有限元參數化模型修正中待定模型參數的選擇問題進行研究，通過對某兩層螺栓連接鋼框架有限元模型修正進行數值仿真與實驗室模型試驗研究，結果表明：本文方法通過信息增益指標定量表征有限元模型修正過程從測量數據中提取用于修正待定模型參數的信息量大小，以懲罰有限元模型待修正參數的參數化復雜程度，能有效權衡有限元模型建模參數復雜程度與其相應的信息論表征復雜度，以獲得模型參數化相對較簡單的待修正有限元模型。此外，對于本文螺栓連接框架結構而言，柱-基礎半剛性連接模型參數的不確定性較梁-柱半剛性連接的不確定性更大，對結構整體特性影響也相應更顯著，這在該類型結構的實際工程設計過程中需要予以注意。應注意到，待修正模型參數先驗分布形式會對信息增益計算結果產生影響，其選取主要依賴使用者經驗，本文選取滿足最大熵原理和中心極限定律的正態分布作為有限元待修正模型參數先驗分布形式，有關模型參數先驗分布影響規律將在后面的研究中予以進一步探討。

參 考 文 獻

[ 1 ] FRISWELL M I, MOTTERSHEAD J E. Finite element model updating in structural dynamics[M]. Dordrecht: Kluwer Academic Publishers, 1995.

[ 2 ] 張德文, 魏阜旋. 模型修正與破損診斷[M]. 北京：科學出版社, 1999.

[ 3 ] 魏來生. 結構有限元動態模型修正方法綜述[J]. 振動與沖擊, 1998, 17(3): 43-46.

WEI Laisheng.A review on the updating methods of finite element model[J]. Journal of Vibration and Shock, 1998, 17(3): 43-46.

[ 4 ] 尹濤, 余嶺, 朱宏平. 一種基于模型修正的結構損傷識別方法[J]. 振動與沖擊, 2007, 26(6): 59-62.

YIN Tao, YU Ling, ZHU Hongping.Model updating based approach for structural damage detection[J]. Journal of Vibration and Shock, 2007, 26(6): 59-62.

[ 5 ] YU L，YIN T. Damage identification in frame structures based on FE model updating[J]. Journal of Vibration and Acoustics, 2010, 132(5): 1741-1757.

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

[ 7 ] FRISWELL M I, MOTTERSHEAD J E, AHMADIAN H. Combining subset selection and parameter constraints in model updating[J]. Journal of Vibration and Acoustics, 1998, 120(4): 854-859.

[ 8 ] 宋漢文, 王麗煒, 王文亮. 有限元模型修正中若干重要問題[J]. 振動與沖擊, 2003, 22(4): 68-71.

SONG Hanwen, WANG Liwei, WANG Wenliang.Several important problems for updating finite element model[J]. Journal of Vibration and Shock, 2003, 22(4): 68-71.

[ 9 ] 費慶國, 張令彌, 李愛群, 等. 基于統計分析技術的有限元模型修正研究[J]. 振動與沖擊, 2005, 24(3): 23-26.

FEI Qingguo, ZHANG Lingmi, LI Aiqun, et al. Finite element model updating using statistics analysis[J]. Journal of Vibration and Shock, 2005, 24(3): 23-26.

[10] 李延強, 杜彥良. 基于最敏感設計參數的斜拉索有限元動力模型修正[J]. 振動與沖擊, 2009, 28(3): 141-143.

LI Yanqiang, DU Yanliang. Dynamic finite element model updating of stay-cable based on the most sensitivity design variable[J]. Journal of Vibration and Shock, 2009, 28(3): 141-143.

[11] MOTTERSHEAD J E, LINK M, FRISWELL M I. The sensitivity method in finite element model updating: a tutorial[J]. Mechanical Systems and Signal Processing, 2011, 25(7): 2275-2296.

[12] WAN H P, REN W X. Parameter selection in finite-element-model updating by global sensitivity analysis using Gaussian process metamodel[J]. Journal of Structural Engineering, 2014, 141(6): 04014164.

[13] 姜東, 丁繼鋒, 費慶國，等. 一種有限元模型修正中的參數選擇方法[C]∥2010全國固體力學大會. 武漢：固體力學學報, 2011.

[14] GULL S F. Bayesian inductive inference and maximum entropy[C]∥ In Maximum Entropy and Bayesian Methods. Dordrecht: Kluwer Academic Publishers, 1989.

[15] HAARIO H, LAINE M, MIRA A, et al. DRAM: efficient adaptive MCMC[J]. Statistics and Computing, 2006, 16(4): 339-354.

[16] YUEN K V. Bayesian methods for structural dynamics and civil engineering[M]. Hoboken: John Wiley & Sons, 2010.

[17] YIN T, JIANG Q H, YUEN K V. Vibration-based damage detection for structural connections using incomplete modal data by Bayesian approach and model reduction technique[J]. Engineering Structures, 2017, 132(1): 260-277.

[18] YIN T, LAM H F, CHOW H M, et al. Dynamic reduction based structural damage detection of transmission tower utilizing ambient vibration data[J]. Engineering Structures, 2009, 31(9): 2009-2019.

[19] LAM H F, YIN T. Dynamic reduction based structural damage detection of transmission towers: practical issues and experimental verification[J]. Engineering Structures, 2011, 33(5): 1459-1478.