基于貝葉斯理論嵌套抽樣的結構物理參數識別研究

王坤陽， 公茂盛， 左占宣

(中國地震局工程力學研究所，中國地震局地震工程與工程振動重點實驗室，哈爾濱 150080)

結構健康監測是目前土木工程領域較為熱門與重要的研究方向，而結構損傷識別是結構健康監測技術的核心。結構損傷通常表現為剛度的退化，因此對結構物理參數(如剛度)進行識別對于結構損傷識別是尤為重要的。常用的根據振型、頻率等模態參數進行識別的確定性方法因為環境噪聲、儀器及復雜的結構屬性的影響使得觀測數據中存在很大的誤差，大大增加了參數識別中的不確定性[1-2]，因此基于不確定性的結構損傷識別方法研究是十分必要的。

1989年，Beck[3]基于貝葉斯框架提出了系統參數識別中處理不確定性問題的貝葉斯方法；Vanik[4]通過概率關系將貝葉斯方法應用到了結構健康監測之中。貝葉斯理論通過結構輸入與輸出的關系，建立后驗聯合概率密度函數來計算待識別參數的概率密度值并進行參數估計。實際結構的待識別參數是十分復雜的，所以在結構物理參數識別中，后驗聯合概率分布通常是復雜、高維的，常用的物理參數識別的貝葉斯估計漸進逼近方法[5](asymptotic approximation methods)面對高維概率分布問題難以求解。Beck等[6]提出了利用馬爾可夫蒙特卡洛(MCMC)抽樣方法來解決高維概率分布問題。Sohn等[7-12]使用該方法對不同結構模型的結構物理參數進行了識別。

傳統的馬爾可夫蒙特卡洛方法面對高維概率密度分布時會出現不收斂、采樣效率低與識別結果誤差較大等問題。Cheung曾通過使用混合蒙特卡洛模擬(hybrid Monte Carlo simulation)的方法代替傳統的馬爾可夫蒙特卡洛抽樣方法，雖然提高了抽樣效率，但仍存在識別結果誤差較大的問題[13]。

Skilling[14]在2004年提出了針對解決貝葉斯公式中高維概率分布問題的嵌套抽樣方法(nested sampling)，不僅能加快抽樣速度，更能根據自己的目標與需求來控制抽樣進程，后來Speagle[15]對該算法進行了優化，提出了自適應嵌套抽樣方法，并給出了嵌套抽樣與自適應的嵌套抽樣的相關算例與開源程序。曹彤彤等[16]通過嵌套抽樣對地下水模型算例進行了研究,驗證了其在模型選擇中的有效性。嵌套抽樣方法往往用于計算邊緣似然值并進行模型選擇，很少直接使用該方法計算整體邊緣概率值并應用在結構的物理參數識別之中。本文提出基于加速度時程響應建立結構物理參數相關的后驗聯合概率密度函數，并通過插入不影響整體概率密度函數值的均勻分布π(θ)，代替原有的先驗分布并更改似然函數的方法，應用嵌套抽樣對結構物理參數進行識別。采用所提方法分別對10層結構數值算例與一個3層振動臺試驗鋼筋混凝土(RC)框架模型，進行結構物理參數識別，證實了方法在貝葉斯理論結構物理參數識別中不僅能很好的解決高維概率密度分布問題，精確估計待識別結構參數，并且在實際結構參數識別中也存在很好的適用性。

1 基本理論與方法應用

1.1 結構物理參數識別貝葉斯方法

對于Nθ個自由度的工程結構，根據最大熵原理，結構剛度參數θ服從正態分布，記θ～N(θMPE, 1)，表達式如下

(1)

(2)

(3)

(4)

根據貝葉斯公式，可以求得基于實測響應的結構物理參數θ的后驗聯合概率密度函數

(5)

(6)

根據文獻[17]的分析結果，式中c=0.2，通過式(5)、(6)即可通過抽樣得到結構參數θ的最優估計值。

1.2 嵌套抽樣方法理論及改進

直接應用嵌套抽樣，在先驗分布中生成樣本，通過似然函數概率值進行選擇，不僅會導致抽樣效率低，還可能因為離開先驗分布限制使得識別結果遠遠偏離真實值的情況出現。為解決這一問題，本文假設有取值范圍足夠大的均勻分布與原后驗聯合概率密度函數相乘并作為新的先驗分布

π(θ)=1

(7)

由式(5)、(6)得

(8)

令原后驗聯合概率密度函數

(9)

代替成為新的似然函數，則有：

(10)

(11)

式(10)中dX是先驗分布累積的微分。定義先驗分布累積為

(12)

0

(13)

如上，通過將θ轉換為X便實現了將多維問題轉換為一維問題。為方便計算式(10),認為第k次抽樣后

(14)

取式(13)作為抽樣終止條件，其中ε為足夠小得常數，通常令ε=0.01[18]。

ΔlnZk=lnZk-lnZk-1<ε

(15)

隨著抽樣進程即可以得到Lmax(θ)，即：

(16)

通過建立新的先驗分布與似然函數利用嵌套抽樣性質，便可以取最后一步抽樣得到的θ值作為近似最優識別值。具體嵌套抽樣算法流程如圖1所示。

圖1 嵌套抽樣算法流程圖

2 結構識別算例

2.1 結構數值模型

因為土木工程結構復雜，貝葉斯理論結構物理參數識別應該考慮多參數的高維問題。本文建立了一個10層剪切結構模型,如圖2所示，具體的結構設計參數如表1所示。

圖2 10層剪切型結構模型示意圖

表1 數值模型結構設計參數及識別結果

算例中，假定結構的阻尼比ζ不變且為0.05，待識別的物理參數有結構各層剛度θ=[k1,k2, …,k10]。對結構輸入白噪聲激勵，通過Newmark-β法計算每一層的加速度時程響應，并在激勵與響應中增加5%水平的白噪聲模擬真實的輸入與響應信息。令π(θ)中θ～U(0,10)，代入式(9)得到后驗聯合概率密度后進行抽樣。

2.2 識別結果分析

根據式(15)，當ε滿足停止條件時，取最后一次抽樣的樣本點，即后驗概率密度近似最大估計值作為剛度識別值，具體識別結果見表1。通過對數值算例結果分析，識別結果最大誤差不超過5%，可見嵌套抽樣能很好的解決基于貝葉斯的結構物理參數識別中的高維后驗聯合概率密度函數問題，且較為精確地識別結構物理參數。本文曾采用不添加噪聲進行識別，識別結果的誤差更小，說明所提方法有很好的抗噪性能，但由于篇幅有限，本文沒有給出不考慮噪聲的識別結果。

11例術后切口愈合良好，術后2周傷口拆線。治療后疼痛均消失，足弓恢復，無跛行，術后隨訪3～18個月，未出現疼痛及其他并發癥，恢復正常的生活及工作。末次隨訪時，AOFAS評分為平均（83.00±2.03）分（77～90分），優 4例，良5例，可2例。典型病例見圖1。

圖3 嵌套抽樣進程

2.3 與其它方法結果對比

工程領域最為常用抽樣方法為MCMC抽樣方法，但當預估結構參數與實際結構參數相差過大或待抽樣參數過多時，常常會導致抽樣結果不收斂或不穩定。本文采用MCMC抽樣中較先進的SCAM算法，對該10層結構數值模型進行了識別，識別結果如表2和圖4所示。將該方法識別結果誤差與本文嵌套抽樣結果誤差進行對比，結果如圖5所示。

表2 SCAM方法識別結果

圖4 SCAM抽樣結果

圖5 誤差對比

通過對識別過程及結果分析可知，嵌套抽樣方法可以通過控制終止條件使抽樣結果更加的穩定，很好地解決了多參數抽樣問題，且識別結果更精確。而從誤差對比結果來看，MCMC方法的最大誤差為14.4%，而本文方法最大誤差僅為4.74%，兩者相差3倍；從誤差的平均值來看，MCMC方法誤差的平均值為5.27%，而本文方法平均誤差僅為2.57%，兩者相差2倍。可見本文方法比傳統的MCMC方法具有更高的識別精度。

3 振動臺RC框架模型參數識別

3.1 振動臺試驗RC框架模型

為了驗證基于貝葉斯理論應用嵌套抽樣對于實際結構物理參數識別的適用性，本文對振動臺試驗中3層RC框架結構模型的層間剛度及阻尼比進行了識別。結構模型及結構尺寸簡圖如圖6所示，具體的結構模型設計及試驗工況等情況可以參考文獻[19]。高艷濱曾通過Opensees有限元建模并對該結構進行Pushover分析，得到了該結構的初始狀態物理參數估計值，具體結果如表3所示。

(a) 振動臺試驗模型

表3 初始狀態結構參數估計值

3.2 振動臺試驗RC框架物理參數識別

通過振動臺試驗中結構加速度響應時程數據對結構物理參數及模態參數進行識別。當結構模型遭受低幅值白噪聲激勵時，一般不會引起結構的非線性反應，其響應主要與質量、剛度和阻尼比有關，假設結構質量保持不變，待識別的結構參數只有物理參數剛度以及模態參數阻尼比便只有剛度及阻尼比，即θ=[k1,k2,k3,ζ1,ζ2]。

該結構振動臺模型試驗中，根據輸入地震動峰值加速度大小，可以分為S0、S1、S2、S4等四種損傷狀態，其中S0是未輸入地震動時模型初始狀態，代表結構模型完好狀態；S1為輸入峰值加速度為0.1g地震動輸入后損傷狀態，結構模型發生輕微損傷；S2為輸入峰值加速度為0.2g地震動輸入后損傷狀態，結構1層損傷嚴重，2層輕微破壞；S3為輸入峰值加速度為0.4g地震動輸入后損傷狀態，柱底混凝土壓碎，2層柱端開始出現裂縫；S4為輸入峰值加速度為0.5g地震動后損傷狀態，底層，2層均破壞嚴重，3層輕微破壞。分別用該結構S0～S4工況的白噪聲響應進行識別，其抽樣識別進程如圖7所示，各個工況識別結果如表4所示。

表4 各工況結構參數識別值

3.3 識別結果分析

從識別結果來看，S0狀態識別結果與表2給出的初始值基本相同，隨著加載工況的進行，結構發生了明顯的剛度退化現象，具體如圖7所示，由圖可以推斷結構在振動臺試驗中首先是第1層發生了破壞，其次是第2層發生破壞，最后是第3層，完全符合觀察得到的試驗現象。另外，本文結果與文獻[20]對該結構試驗模型的識別結果相比差別不大，也驗證了本文識別結果的正確性。

圖8 不同工況剛度識別值

為了查看識別效果，本文給出了識別結果和原始記錄對比情況。因篇幅有限，本文只給出了S0和S4狀態記錄對比結果，如圖9和圖10所示，其余結果與此類似。通過圖中計算加速度與觀測加速度的對比可知，模擬識別的加速度時程與觀測記錄的加速度吻合度及識別精度還是較高的。

圖9 S0工況模擬加速度與觀測加速度對比

圖10 S4工況模擬加速度與觀測加速度對比

4 結 論

本文通過修改先驗分布與似然函數的方法對貝葉斯理論嵌套抽樣進行改進，給出了基于加速度時程響應的結構物理參數識別方法，并分別對數值模型和實際振動臺試驗結構模型進行物理參數識別。得到結論如下：

(1) 相較于MCMC抽樣方法，嵌套抽樣能很好地解決結構物理參數識別中的高維后驗聯合概率密度函數問題，即使噪聲水平較高，也能精確的識別出結構物理參數。

(2) 通過對結構數值模型結構物理參數抽樣來看，從抽樣進程可以看出該方法十分高效與可靠，從識別結果來看也具有很好的抗噪性。

(3) 通過對實際振動臺試驗結構模型的識別結果顯示，方法可以有效識別實際結構的物理參數，表明方法可以用在實際工程的結構損傷識別之中。

本文只是對低幅值白噪聲激勵下的結構物理參數剛度及模態參數阻尼比進行了識別，實際結構在較大荷載(如地震)作用下，往往會發生非線性反應，后續工作中將針對這一情況開展進一步研究。