基于Gibbs抽樣的結構時域載荷識別

王 婷， 萬志敏， 鄭偉光

(1. 南通職業大學 機械工程學院，江蘇 南通 226000； 2. 華中科技大學 機械科學與工程學院，武漢 430074)

在許多實際工程問題中，精確地掌握作用于結構上的外載荷對結構優化設計、響應重構、故障診斷以及健康監測等起著非常重要的作用。然而，很多情況下的載荷根本無法直接測量，而結構響應卻往往較易獲得，因此載荷識別技術得以發展，主要包括頻域法和時域法兩大類[1]。頻域法首先獲得頻域載荷，再通過傅里葉變換獲得時域載荷，故頻域法對沖擊載荷識別效果不佳。時域法[2-5]直接在時域內識別載荷的時間歷程，故不受載荷形式限制，且識別結果具有明確的物理意義，無變換誤差，在工程中具有良好的應用前景。

另一方面，載荷識別是動力學的第二類反問題，其中的結構矩陣求逆過程往往是病態的。正則化技術被國內外學者廣泛用于獲得該病態問題的穩定解，而其中的關鍵在于正則化參數的選取。Law等[6]基于時域法識別了移動載荷，并引入Tikhonov正則化技術尋求近似解，其中正則化參數采用L-曲線法獲得。Mao 等[7]建立了基于精細積分的載荷識別模型，并采用GCV準則獲得了正則化參數。馬超等[8-9]分別采用一種基于局部最小準則及改進的正則化技術識別了載荷。劉杰[10]基于Green核函數建立時域載荷識別模型，并比較了不同正則化參數的選擇方法，包括L-曲線、GCV準則等。

然而，上述傳統方法獲得的正則化參數均為恒定值，往往很難適用于不同類型載荷的反演。Zhang等[11]提出了貝葉斯框架下的載荷識別頻域方法，該法具有本征的自適應正則化性能。本文將該法拓展到載荷識別的時域法，將載荷與測量噪聲看作是隨機變量，建立載荷識別的貝葉斯模型，并采用Gibbs抽樣法獲得了載荷后驗值。對照Tikhonov正則化的數學程式，可知本文提出的載荷識別法具有可變的正則化參數。數值算例分別采用六自由度彈簧質量模型與四邊簡支板模型為對象，驗證了本文方法的有效性。

1 載荷識別的狀態空間模型

對于線性時不變動態系統，其離散狀態空間方程可以表示為

xt+1=Axt+Bft

(1)

yt=Rxt+Dft

(2)

式中：t∈[1,nt]表示離散時刻；xt為狀態空間變量；ft為外載荷向量；A,B,R和D分別為系統矩陣、輸入矩陣、輸出影響矩陣及輸出傳遞矩陣；yt為響應向量。

假設初始狀態x1為零， 同時令H1=D,Ht=R(A)nt-1B， 將式(1)代入式(2)中，則可得矩陣形式方程如下

Y=Hf

(3)

有H=

2 基于貝葉斯理論的動載荷識別

由于實測的響應都是存在噪聲的，故實際載荷識別方程應為

Y=Hf+W

(4)

式中：W為噪聲向量。

本文采用貝葉斯思想來求解上述動載荷識別問題，即通過后驗概率密度分布來定量未知載荷時間歷程。首先定義向量θ為不確定性參數，包括載荷大小、測量噪聲，那么根據貝葉斯原理，聯合概率后驗概率密度函數可以寫成

p(θ|Y)=p(Y|θ)p(θ)/p(Y)

(5)

2.1 貝葉斯框架下的載荷識別

根據中心極限定理，測量噪聲可以近似表示為均值為0的正態分布。例如，第k個響應的噪聲分布為

(6)

式中： 符號～表示“分布為”。那么似然函數可以寫成

(7)

假設時域載荷向量為共軛先驗分布，則概率分布為

f～N(f0,Cf0)

(8)

(9)

(10)

(11)

式中： 超參數{kW,αW,U0,σU0,kf,αf}均假設已知。結合上述超先驗分布，得到了一個三層貝葉斯載荷識別模型：其一為似然函數層；其二為載荷先驗模型層；其三為底部超參數層。由于分層貝葉斯模型對于超參數的取值非常穩定，故對于上述共軛先驗中的超參數可以取任意值。

2.2 Gibbs抽樣法求解不確定參數

根據上述的似然函數以及參數先驗分布，可以得到所有不確定性參數的聯合后驗概率密度分布如下：

(12)

其中先驗分布

(13)

式(13)是由于各個參數之間相互獨立。另外，由于本文采用的先驗分布以及似然函數都為標準分布(正態與伽馬分布)，那么各個參數的條件概率密度函數分別可以推導出為標準分布形式，具體推導過程如下：

(1) 載荷的條件概率密度函數

(14)

對式(14)右邊取對數，則

(15)

對式(15)取載荷f的一階偏導求極值，得到極值點為

(16)

對式(15)取載荷f的二階偏導求極值，得到極值點為

(17)

(18)

(2) 測量噪聲與模型誤差方差的條件概率密度函數

(19)

(3) 載荷先驗均值的條件概率密度函數

(20)

(4) 載荷先驗方差的條件概率密度函數

(21)

學術界現有眾多的MCMC(Markov Chain Monte Carlo)算法[12-14]，由于本文推導的所有條件后驗概率分布均為標準分布，如正態分布、伽馬分布，故滿足Gibbs抽樣法的應用條件。綜上，基于Gibbs抽樣的載荷識別算法總結如下：

⑥返回第②步，重復迭代抽樣，直到抽樣到足夠多的粒子。

(22)

(23)

2.3 貝葉斯方法與正則化方法比較

由上述貝葉斯方法可得結構動載荷的后驗解為

(24)

可以看出式(24)與傳統正則化方法有些許聯系。標準Tikhonov法的實質是求解滿足最小均方差的載荷，

(25)

即動載荷的正則化解形式如下：

(26)

3 數值分析

采用兩種結構模型來驗證本文方法的可行性，并與Tikhonov正則化方法進行對比。兩種仿真結構模型分別為六自由度彈簧質量結構模型與四邊簡支板結構模型。

3.1 算例 1

六自由度的彈簧質量結構，如圖1所示，共包含6個質量單元和11個彈簧單元，其質量和剛度參數如表1。假設結構為Rayleigh阻尼，結構的6階模態頻率分別為9.76 Hz, 53.26 Hz, 72.54 Hz, 83.51 Hz, 85.58 Hz和162.28 Hz。

圖1 彈簧質量結構模型Fig.1 Mass-spring structural model

表1 結構的質量剛度值

載荷作用于質量塊2上，采用如下形式

F1(t)=50(1-cos(20πt))sin(45πt)

(27)

測量響應采用質量塊3、5的加速度響應，并將噪聲水平為5%的測量噪聲加入測量信號中。采樣頻率采用1 000 Hz，仿真時間為0.3 s。此時，Markov矩陣H的條件數為1.284 5E+008，非常大的條件數意味著矩陣呈病態特性。

針對上述載荷反演不適定性問題，采用Tikhonov正則化技術消除矩陣奇異性，獲得載荷的正則化穩定解。首先，采用L-曲線準則獲得了正則化參數為λ=2.508 6E-005。進而把正則化參數值代入Tikhonov方程(26)計算出載荷識別值，如圖2所示。此時雖然載荷識別值不發散，但與真實值的誤差仍然較大，說明正則化參數的選取不合適。

然后，采用另一種常用的正則化參數選擇準則——GCV準則，來選取正則化參數。GCV選擇參數值為λ=0.007 020 7。采用該參數，代入Tikhonov方程(26)計算出載荷識別值，如圖3所示。由圖可以看出，載荷識別曲線基本和真實曲線一致，整體誤差較小，僅曲線拐點處誤差較大，總體正則化效果要比上述L-曲線正則化方法好得多。另外，對比上述兩種準則選取的參數值，發現通過GCV準則選擇的參數值要遠大于通過L-曲線準則選取的參數值，說明L-曲線選取的參數值偏小，導致矩陣仍然病態。

圖2 載荷識別值與真實值(L-曲線)Fig.2 The actual and identified values of the load (L-curve)

圖3 載荷識別值與真實值(GCV準則)Fig.3 The actual and identified values of the load (GCV criterion)

圖4 載荷識別值與真實值(Gibbs抽樣)Fig.4 The actual and identified values of the load (Gibbs sampling)

圖5 三種正則化參數對比圖Fig.5 The comparison results of the regularization parameter based on these three methods

3.2 算例 2

采用如圖6所示的簡支板為驗證對象，力錘沖擊載荷作用于薄板中心(325 mm，225 mm)處，采用兩個加速度信號(位置如圖6，黑色方框代表加速度傳感器)作為測量信號來識別該沖擊載荷。

不采用正則化技術時，沖擊載荷識別曲線如圖7所示，可以看出載荷識別曲線從開始階段到結束階段逐漸發散，表明該反求系統病態，故采用Tikhonov正則化技術來獲得其穩定解。采用了L-曲線準則和GCV準則來獲取正則化參數分別為0.014 07、0.890 1，那么反求的沖擊載荷分別如圖8和圖9所示，可以看出采用了正則化技術后，識別載荷已不再發散，但由于Tikhonov正則化的存在，反求結果已經偏離了真實解，特別是GCV準則下的正則化參數過大，導致了沖擊載荷峰值與真實值的誤差較大。

圖6 激勵位置以及傳感器布置Fig.6 The load and sensors layout

采用本文方法來識別該沖擊載荷，超參數設置為(kW=10,αW=1)、(U0=0,σU0=1)、(kf=10,αf=0.1)，得到載荷識別值如圖10所示，可以看出載荷識別值既穩定又準確，識別的曲線能夠非常接近真實載荷曲線，特別是沖擊峰值比傳統正則化技術的識別精度要高的多。這是因為本文的載荷識別方法具有本征自適應正則化屬性，其變化的正則化參數曲線如圖11所示，可以看出該正則化參數曲線類似沖擊載荷曲線，整體呈沖擊狀，僅在載荷沖擊時段具有較大的值，在沖擊峰值達到參數峰值。本算例再次說明了基于Gibbs抽樣的動載荷識別方法相較于傳統正則化方法具有自適應正則化的優越性，適合各種類型的載荷識別問題。

圖7 沖擊載荷識別值與真實值(無正則化)Fig.7 The actual and identified values of the impact load (No regularization)

圖8 沖擊載荷識別值與真實值(L-曲線)Fig.8 The actual and identified values of the impact load (L-curve)

圖9 沖擊載荷識別值與真實值(GCV準則)Fig.9 The actual and identified values of the impact load (GCV criterion)

圖10 沖擊載荷識別值與真實值(Gibbs抽樣)以及局部放大圖Fig.10 The actual and identified values of the impact load (Gibbs sampling)

圖11 正則化參數對比圖Fig.11 The comparison results of the regularization parameter based on these two methods

4 結 論

本文在貝葉斯框架下建立了動載荷識別的Gibbs抽樣方法。該法具有本征的自適應正則化性能，正則化參數即為噪聲方差與載荷方差的比值。數值算例驗證了本文方法的有效性，并與Tikhonov正則化方法對比，結果表明本文方法的精度要高于Tikhonov正則化法。當前的載荷識別方法大都以準確的結構模型為基礎，而本文的方法可以進一步考慮結構的模型誤差，此為作者下一步的研究內容。

[ 1 ] INOUE H, HARRIGAN J J, REID S R. Review of inverse analysis for indirect measurement of impact force [J]. Applied Mechanics Reviews, 2001, 54(6): 503-524.

[ 2 ] 孫興盛，劉杰，丁飛，等. 基于矩陣攝動的隨機結構動態載荷識別技術[J]. 機械工程學報, 2014, 50(13): 148-156.

SUN Xingsheng, LIU jie, DING Fei, et al. Identification method of dynamic loads for stochastic structures based on matrix perturbation theory [J]. Journal of Mechanical Engineering, 2014, 50(13): 148-156.

[ 3 ] 楊智春，賈有. 動載荷的識別方法[J]. 力學進展, 2015: 29-54.

YANG Zhichun, JIA You. The identification of dynamic loads [J]. Advances in Mechanics, 2015: 29-54.

[ 4 ] 朱濤，肖守訥，陽光武. 一種新的時域動態載荷識別方法 [J]. 西南交通大學學報, 2012, 47(6): 968-973.

ZHU Tao, XIAO Shoune, YANG Guangwu. A new time domain method for force identification [J]. Journal of Southwest Jiaotong University, 2012, 47(6): 968-973.

[ 5 ] 徐菁，張方，姜金輝，等. 基于擬靜態初值的載荷識別數值修正算法[J]. 振動與沖擊, 2016, 35(2): 39-44.

XU Jing, ZHANG Fang, JIANG Jinhui, et al. Numerical correcting algorithm for load identification based on quasi-static initial value [J]. Journal of Vibration and Shock, 2016, 35(2): 39-44.

[ 6 ] LAW S S, FANG Y L. Moving force identification: optimal state estimation approach [J]. Journal of Sound and Vibration, 2001, 239(2): 233-254.

[ 7 ] MAO Y M, GUO X L, ZHAO Y. A state space force identification method based on Markov parameters precise computation and regularization technique [J]. Journal of Sound and Vibration, 2010, 329(15): 3008-3019.

[ 8 ] 馬超，華宏星. 正則化技術在狀態空間載荷識別中的應用 [J]. 振動. 測試與診斷, 2014, 34(6): 1154-1158.

MA Chao, HUA Hongxing. State space force identification based on regularization technique [J]. Journal of Vibration, Measurement & Diagnosis, 2014, 34(6): 1154-1158.

[ 9 ] 馬超，華宏星. 基于改進正則化方法的狀態空間載荷識別技術[J]. 振動與沖擊, 2015, 34(11): 146-149.

MA Chao, HUA Hongxing. State space load identification technique based on an improved regularized method [J]. Journal of Vibration and Shock, 2015, 34(11): 146-149.

[10] 劉杰. 動態載荷識別的計算反求技術研究[D]. 長沙：湖南大學, 2011.

[11] ZHANG E, ANTONI J, FEISSEL P. Bayesian force reconstruction with an uncertain model [J]. Journal of Sound and Vibration, 2012, 331(4): 798-814.

[12] HASTINGS W K. Monte Carlo sampling methods using Markov chains and their applications [J]. Biometrika, 1970, 57(1): 97-109.

[13] GEMAN S, GEMAN D. Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images [J]. IEEE Transactions on Pattern Analysis and Machine Intelligence ,1984(6): 721-741.

[14] GILKS W R. Introducing markov chain monte carlo, markov chain monte carlo in practice [M](eds. WR Gilks, S. Richardson and DJ Spiegelhalter), 1—19. Chapman & Hall, London, 1996.