基于改進SCAM算法的結構物理參數識別貝葉斯方法

趙一男，公茂盛，左占宣，高艷濱

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

結構健康監測是目前土木工程領域重要研究方向，而結構損傷識別是結構健康監測技術的核心。結構損傷表現為結構物理參數的變化(如剛度降低)，或由此導致的模態參數的變化(如頻率下降)，因此通過對損傷前后的結構物理參數或模態參數進行識別可以判斷損傷位置和損傷程度。損傷識別方法一般分為確定性方法和不確定性方法[1-2]，實際工程中，由于結構本身屬性復雜、環境噪聲[3]等不確定性因素影響，使得采用基于頻率、振型等確定性方法[4-5]識別結構損傷的準確性和可靠性降低。于是學者們發展確定性方法的同時提出了不確定性方法，而基于貝葉斯估計的結構損傷識別方法因為可以很好地解釋損傷識別過程中所存在的不確定性，成為目前廣受關注和應用的不確定性方法之一。

早在1989年，Beck[6]基于貝葉斯統計框架提出了一種系統識別方法，并明確指出了識別過程中存在不確定性?；谪惾~斯估計的結構損傷識別可以給出結構參數后驗聯合概率分布，通過求解參數后驗邊緣概率分布和最優估計值判斷損傷位置和損傷程度，求解后驗聯合概率分布通常需要積分運算，但實際情況中，結構參數的后驗聯合概率分布是復雜的、高維的，很難通過直接積分求解邊緣概率分布和最優估計值[7]。為解決這一問題，Hastings等[8-11]提出和選擇馬爾可夫蒙特卡羅(Markov Chain Monte Carlo，MCMC)抽樣方法，通過在概率空間中隨機抽樣來逼近參數的后驗分布[12]。但實踐發現，MCMC算法仍存在抽樣收斂速度慢，不能求解高維概率分布等問題，為此Haario等[13]提出了自適應Metropolis算法，提高了傳統MCMC算法的收斂速度，并指出合理選擇提議分布是算法收斂的關鍵，之后又對自適應Metropolis算法進行了改進，提出了逐分量自適應Metropolis(single component adaptive Metropolis，SCAM)算法，該算法可以有效求解高維問題[14]。Simoen等[15]利用SCAM算法對7層結構的損傷進行了識別，表明了SCAM算法的實用性，但SCAM算法的提議分布方差序列不滿足馬爾可夫性質，特別是當選取的抽樣參數初始值與要識別的真實值相差較大且提議分布方差較小時，抽樣效率低、易生成重復樣本序列，從而導致結果誤差較大[16]。

盡管如此，SCAM算法仍然是目前解決貝葉斯估計中求解高維邊緣概率分布的常用方法，廣泛地應用于基于貝葉斯估計的結構損傷識別。本文首先針對SCAM算法的不足，重新定義了提議分布方差計算式，給出了改進的SCAM算法，進一步提高了抽樣效率和計算結果準確性。然后利用貝葉斯理論給出結構物理參數的后驗聯合概率分布，再結合改進的SCAM算法求解結構物理參數的后驗邊緣概率密度分布和最優估計值，從而實現對結構損傷進行識別和評估。最后通過對一個4層數值結構模型進行了識別，驗證了本文所提出算法的有效性和可靠性。

1 基本理論與方法改進

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

對于Nθ個自由度的工程結構，根據最大熵原理，結構剛度參數θ服從伽瑪分布，記θ～Γ(α,β)，表達式為

(1)

式中：θMPE為振動測試前結構剛度參數極大先驗估計值；α為形狀因子；β為逆尺度因子。α,β與參數θ的數學期望μ和標準差σ有關，分別由式(2)和式(3)確定

(2)

(3)

式中，cvi為變異系數。

(4)

(5)

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

(6)

式中，C為標準化常量，保證后驗概率密度函數在定義域內積分為1。通過式(6)即可得到結構參數θ的后驗概率分布和最優估計值。

1.2 SCAM抽樣算法及改進

為了求解式(6)后驗聯合概率分布的邊緣分布，Haario等提出的逐分量自適應Metropolis算法是目前較為有效的一種抽樣方法，得到了較為廣泛地應用。但在SCAM算法中，當前狀態下提議分布方差并不只由上一狀態的樣本候選值決定，而是由之前所有狀態的樣本候選值共同決定，這樣生成的隨機樣本序列將不再滿足馬爾可夫鏈無后效性，不能完全保證樣本序列收斂，會生成大量病態樣本，導致結果誤差較大，因此提議分布的方差直接影響算法收斂性。

研究表明，無論提議分布方差過大還是過小都會導致算法抽樣效率變低，所得最終結果較不可靠。所以必須選擇合適的提議分布方差使得馬爾可夫鏈穩定。在提議分布為正態分布的多維問題中，Roberts等[17-20]給出了1～10維情況下的最優接受率取值，且證明當維度逐漸增加時，最優接受率逐漸減小并趨近于0.234。因此，算法改進的核心問題在于如何合適地調整提議分布方差，使得抽樣過程的接受率收斂到最優接受率。

根據上述分析及最優接受率取值，本文對SCAM算法中的提議分布方差進行重新定義，改進了SCAM算法。新的提議分布方差由式(7)確定

νt=νt-1(1+α(Xt-1,z)-αOAR)2

(7)

式中：νt為t+1時刻提議分布方差；αOAR為最優接受率;α(Xt-1,z)為樣本值為Xt-1，候選值為z時的接受率。由式(7)可知，若取αOAR≈0.234，則α(Xt-1,z)-αOAR∈[-0.234,0.766]，下一狀態的提議分布方差取值范圍為νt∈[0.7662νt-1,1.7662νt-1]。當α(Xt-1,z)≤αOAR時,νt會減小到0.7662νt-1～νt-1；當α(Xt-1,z)>αOAR時，νt會增大到νt-1～1.7662νt-1，隨著迭代步數增加不斷調節提議分布方差的數值，使得接受率趨近于最優。此外，為防止提議分布方差過大或過小，可將其限制在合理區間內。這樣對提議分布方差改進后，SCAM算法在實現自適應性同時，也使得抽樣樣本序列滿足馬爾可夫鏈無后效性。改進的SCAM算法計算流程如圖1所示。

圖1 改進的SCAM算法流程圖Fig.1 Flow chart of improved SCAM algorithm

1.3 算法驗證

本文建立了一個抽樣算例，對二維隨機向量{x1,x2}獨立同分布于α=6，β=1的伽瑪分布進行隨機抽樣，峰值點的橫坐標作為標準值，對改進的SCAM算法可靠性進行驗證。

二維隨機向量{x1,x2}聯合概率密度函數為

(8)

初始向量值選擇[40,40]，初始提議分布方差取[0.1,0.1]，迭代步數定為10 000步，“burn-in”階段步數定為1 500步，分別用SCAM算法和本文改進的SCAM算法對式(8)進行隨機抽樣。

隨機變量x1的馬爾可夫鏈歷史路徑圖如圖2所示，可以看出原始SCAM算法在“burn-in”階段的接受率較高，而在“burn-in”階段之后的接受率較低，出現大量重復樣本，馬爾可夫鏈穩定性較差。原因在于SCAM算法“burn-in”階段的提議分布方差為固定值，如果選擇不合理，計算效率將大大降低，“burn-in”階段之后的提議方差不滿足馬爾可夫鏈無后效性，計算結果誤差較大。本文改進后SCAM算法的提議分布方差一直受接受率的影響而不斷調整，使得馬爾可夫鏈具有較好的穩定性，即使選擇的初始值與標準值相差較大也可以很快達到收斂。圖3為上述兩種算法所得隨機變量x1的累計均值圖，可以看出，改進的SCAM算法收斂速度明顯快于SCAM算法，SCAM算法在2 500步左右出現收斂趨勢，而改進的SCAM算法則在500步左右就能達到收斂。圖4給出了運用上述兩種算法進行隨機抽樣所得隨機變量x1樣本邊緣概率密度曲線與理論邊緣概率密度曲線的對比情況，從圖中可以看出SCAM算法得到的樣本概率密度曲線與理論概率密度曲線差別較大，而改進的SCAM算法得到的樣本概率密度曲線與理論概率密度曲線基本吻合。

圖4 隨機變量x1樣本邊緣概率密度曲線與理論邊緣概率密度曲線對比Fig.4 Comparison of the sample and theoretical marginal probability density curves of random variable x1

圖3 隨機變量x1累計均值隨迭代步數變化趨勢Fig.3 Variation trend of the cumulative mean of random variable x1 with the number of iteration steps

圖2 隨機變量x1馬爾可夫鏈路徑圖Fig.2 Path diagram of the Markov chain of random variable x1

綜上所述，本文改進的SCAM算法克服了原始SCAM算法存在的不足，具有較高的抽樣效率和可靠性。

2 結構識別算例

2.1 結構數值模型

考慮到土木工程結構的損傷大多為累積損傷，建立了一個4層剪切型結構數值模型，假定了多個損傷狀態，如圖5所示，對該模型在多次累積損傷狀態下的物理參數進行識別。結構設計參數見表1，損傷狀態S0～S4見表2，損傷程度依次為S0

表2 損傷狀態表Tab.2 Damage cases

表1 結構設計參數Tab.1 Structural parameters

圖5 四層結構模型Fig.5 4-story structure model

算例中，假定各損傷狀態下的阻尼比ζ均為0.05，各層質量m不變，待識別的物理參數只有結構各層剛度θ(k)=[k1,k2,k3,k4]。各層剛度參數變異系數取0.2，以保證先驗分布平滑且分布較寬，通過式(6)確定剛度參數的后驗聯合概率分布。對結構基底施加峰值加速度為0.05g的白噪聲激勵，采用Newmark-β法計算結構響應，將各層加速度響應作為觀測值，并分別對無噪聲、5%噪聲和10%噪聲三種噪聲水平下的各工況層間剛度進行識別。初始剛度參數極大先驗估計值取S0損傷狀態剛度參數理論值，后續各狀態剛度參數的極大先驗估計值取上一狀態的極大后驗估計值，通過式(6)求得剛度參數后驗分布，最后利用圖1的計算流程，結合改進的SCAM算法，對各層剛度參數進行逐分量抽樣，求得剛度參數后驗邊緣概率分布，并取剛度參數后驗邊緣概率密度曲線峰值點的橫坐標作為最優估計值。

2.2 識別結果分析

本文采用改進的SCAM算法計算剛度參數的后驗邊緣概率分布，各損傷狀態下各層剛度參數后驗邊緣概率分布如圖6所示。本文只給出了5%噪聲水平下的剛度參數邊緣概率分布，無噪聲和10%水平噪聲下的識別結果與之類似，不再給出。圖中可以清晰地看出結構各損傷狀態下各層剛度參數的后驗概率分布和最優估計值隨損傷程度增加的變化情況。

圖6 各損傷狀態下各層剛度參數邊緣概率分布(5%噪聲)Fig.6 Marginal probability distribution of stiffness parameters of each story ineach damage cases (noise level 5%)

各層剛度參數理論值和最優估計值變化趨勢如圖7所示，從圖中可以看出，各層剛度參數識別值變化趨勢與理論值變化趨勢基本一致，說明本文改進的SCAM算法具有較高的計算準確性。圖8給出了不同噪聲水平下各損傷狀態的各層剛度參數極大后驗估計值與理論值的相對誤差。由圖可知，無噪聲下相對誤差在±4%以內；5%噪聲下相對誤差在±5%以內,10%噪聲下相對誤差在±6%以內。本文所提方法不僅對中等損傷和大損傷(如損傷30%，50%)識別效果較好，而且對小損傷(損傷5%)的識別效果同樣可以接受。從整體上看，隨著噪聲水平的增加，識別結果相對誤差存在增大的趨勢，但在土木工程中仍在可接受范圍內。

圖8 各損傷狀態下各層剛度參數相對誤差Fig.8 Relative errors of the stiffness parameters of each story foreach damage case

圖7 各層剛度參數理論值與最優估計值變化趨勢對比Fig.7 Comparison of the variation trends between the theoretical value and the optimal value of the stiffness parameters of each story

圖9給出了各損傷狀態下各層剛度后驗變異系數變化趨勢，變異系數越小表示識別結果不確定性越小，其可信度越高。在三種噪聲水平下，底層剛度后驗變異系數普遍小于其他三層的剛度后驗變異系數，但各層剛度后驗變異系數均小于先驗變異系數0.2，表明改進方法有效性的同時，也表明基于貝葉斯估計結構物理參數識別方法確實可以降低參數的不確定性。

圖9 各損傷狀態下各層剛度參數變異系數變化趨勢Fig.9 Variation trend of stiffness parameters variation coefficient of each story foreach damage case

該算例表明，基于改進SCAM算法結構物理參數識別貝葉斯方法具有較高的識別精度，能夠準確地識別結構物理參數及其變化,且具有較好的抗噪能力，驗證了改進的SCAM算法有效性和可靠性。

3 結 論

本文通過對SCAM算法中提議分布方差計算方法進行重新定義，提出了改進的SCAM算法，并通過理論分析和數值模擬算例對本文方法的有效性和可靠性進行了驗證。所得結果和結論如下：

(1)通過重新定義提議分布方差所改進的SCAM算法，抽樣樣本序列所構成的馬爾可夫鏈達到了平穩分布，且沒有明顯的趨勢性和周期性，在提高計算準確性的同時提高了計算效率。

(2)基于改進SCAM算法的結構物理參數識別貝葉斯方法可以較好地識別結構物理參數，得到物理參數的后驗邊緣概率分布和最優估計值，識別精度及效果良好，且具有較好的抗噪能力，該方法可以應用于實際工程結構物理參數識別及損傷評估。

本文改進的SCAM算法僅在基本理論、數值模擬中進行了驗證，作者下一步將對實際工程結構及結構振動臺試驗模型進行分析，進一步驗證方法的可靠性、有效性及實用性。