雷 博，鄭 輝，2

（1.上海交通大學(xué) 振動、沖擊、噪聲研究所，上海 200240；2.高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心，上海 200240）

載荷信息的獲取是系統(tǒng)動力學(xué)響應(yīng)計算、動態(tài)設(shè)計和故障分析的重要組成部分，在實際工程中具有重要的意義。但由于工程環(huán)境中存在諸多限制條件，使得多數(shù)情況下，載荷無法直接測量[1]。因此，建立正確有效的載荷識別方法是進(jìn)行后續(xù)動力學(xué)分析和動態(tài)設(shè)計的前提。

載荷識別是指利用系統(tǒng)響應(yīng)和系統(tǒng)特性對系統(tǒng)所受載荷進(jìn)行求解，屬于動力學(xué)第二類反問題。Thite等提出利用直接求逆法和最小二乘正則化方法進(jìn)行設(shè)備噪聲傳遞路徑上的載荷識別，后者利用正則化方法解決了傳遞矩陣的病態(tài)性問題[2–3]。楊立娟等利用Tikhonov正則化方法進(jìn)行懸臂殼載荷識別，并通過混合傳遞路徑分析確定了該結(jié)構(gòu)的主要振動傳遞路徑[4]。張方等推導(dǎo)了基于時間有限元的動載荷識別模型，將時域動載荷識別的逆卷積關(guān)系轉(zhuǎn)變?yōu)閺V義正交域的線性算子逆運(yùn)算，大大降低了該類問題的復(fù)雜程度[5]。隨著計算機(jī)技術(shù)的發(fā)展，載荷識別也可通過智能算法進(jìn)行。智能算法可避免傳遞矩陣的病態(tài)性問題，同時可降低運(yùn)算量。如Yan等利用微遺傳算法對復(fù)合加筋板所受沖擊載荷進(jìn)行了識別，結(jié)果表明該算法具有較高的計算效率[6]。沙瑞華利用三種神經(jīng)網(wǎng)絡(luò)算法對水電機(jī)組所受載荷進(jìn)行了識別，發(fā)現(xiàn)基于LM優(yōu)化算法的BP網(wǎng)絡(luò)對各類動載荷識別的效果優(yōu)于附加動量法BP網(wǎng)絡(luò)和RBF網(wǎng)絡(luò)[7]。

由此可見，對動載荷識別的研究已取得很多成果，然而考慮動載荷識別過程中由于測量誤差等不確定因素所帶來識別誤差的研究尚不多見。為了降低不確定性因素帶來的動載荷識別誤差，本文從系統(tǒng)運(yùn)動方程出發(fā)，考慮測量誤差這一不確定性因素，利用貝葉斯估計理論，結(jié)合馬爾可夫蒙特卡羅方法，建立了基于貝葉斯估計的動載荷識別方法，并通過仿真和實驗數(shù)據(jù)對比對該方法的有效性進(jìn)行了驗證。

1 載荷和誤差參數(shù)的概率分布

1.1 貝葉斯估計理論[8–9]

傳統(tǒng)的統(tǒng)計推斷方法利用統(tǒng)計量的先驗信息進(jìn)行推斷，而貝葉斯估計理論同時考慮統(tǒng)計量的先驗信息和觀測數(shù)據(jù)中包含的有效信息。

將系統(tǒng)所受載荷看作系統(tǒng)參數(shù)，對系統(tǒng)參數(shù)進(jìn)行估計的貝葉斯公式為

式中M表示給定的結(jié)構(gòu)模型，D表示觀測數(shù)據(jù)，θ表示系統(tǒng)參數(shù)，p(θ|D,M)為給定模型和觀測數(shù)據(jù)下參數(shù)的后驗概率分布，p(D|θ,M)為給定模型在輸入D下參數(shù)的似然函數(shù)（觀測數(shù)據(jù)信息），p(θ|M)為給定模型下參數(shù)的先驗概率分布（先驗信息），通常由經(jīng)驗獲得，p(D|M)為觀測數(shù)據(jù)的邊緣概率分布，與系統(tǒng)參數(shù)無關(guān)，起正則化作用，其計算公式為

貝葉斯估計的實質(zhì)是在利用觀測數(shù)據(jù)對系統(tǒng)參數(shù)的先驗概率分布進(jìn)行修正，進(jìn)而得到系統(tǒng)參數(shù)的后驗概率分布，以此作為估計推斷的依據(jù)。

1.2 概率分布估計

系統(tǒng)的運(yùn)動方程為

式中Fm=[F1，F(xiàn)2，…Fm]T表示系統(tǒng)所受載荷的頻域列向量，Yn表示系統(tǒng)響應(yīng)的頻域列向量，Hn×m表示系統(tǒng)傳遞矩陣，e表示由不確定因素引起的預(yù)測誤差，服從均值為零、協(xié)方差矩陣為σe2In×n的正態(tài)分布。

假設(shè)系統(tǒng)所受載荷的先驗概率分布分別為均勻分布，即

進(jìn)一步假設(shè)誤差參數(shù)σe2的先驗概率分布為倒伽馬分布IG(α，β)，即

式中α為形狀參數(shù)，β為尺度參數(shù)。

將系統(tǒng)所受載荷看作系統(tǒng)參數(shù)，由式（2）、式（4）、式（5）、式（6）可得系統(tǒng)所受載荷和誤差參數(shù)的聯(lián)合后驗分布為

對傳遞矩陣進(jìn)行奇異值分解為

式中，上標(biāo)H表示共軛轉(zhuǎn)置。

由式（7）、式（8）可得系統(tǒng)所受載荷和誤差參數(shù)的邊緣概率分布分別為

2 載荷抽樣估計

2.1 馬爾可夫蒙特卡羅方法[10]

馬爾科夫蒙特卡羅方法是根據(jù)統(tǒng)計抽樣理論獲取近似解的一種統(tǒng)計試驗方法，即構(gòu)造平穩(wěn)分布為目標(biāo)概率分布的馬爾科夫鏈，通過隨機(jī)抽樣獲得容量足夠大的樣本。根據(jù)大數(shù)定律，每個參數(shù)的樣本值分布近似于其邊緣概率分布，因此以該樣本的統(tǒng)計特征（期望）作為滿足目標(biāo)概率分布的參數(shù)估計值。在本文建立的方法中，其目標(biāo)概率分布為載荷和誤差參數(shù)的聯(lián)合后驗分布。

2.2 Gibbs抽樣算法[11]

馬爾科夫蒙特卡羅方法常用的抽樣算法為Metropolis-Hastings(MH)抽樣算法和Gibbs 抽樣算法。

MH抽樣算法在每次抽取樣本后需要計算接受概率，根據(jù)接受概率的大小判斷是否接受本次抽取的樣本，若接受，則在本次抽取樣本的基礎(chǔ)上繼續(xù)抽取，若不接受，則需要重新抽取，直至接受，抽樣效率較低。而Gibbs抽樣算法可看作MH抽樣算法中接受概率為1的特例，即每次抽取的樣本均被接受，大大提高了抽樣效率。

構(gòu)造目標(biāo)概率分布為載荷和誤差參數(shù)的聯(lián)合后驗概率分布的馬爾可夫鏈，并通過Gibbs抽樣算法進(jìn)行抽樣的具體流程如圖1所示。

圖1 Gibbs抽樣算法流程圖

3 仿真和實驗驗證

3.1 仿真算例

在基于本文方法進(jìn)行載荷識別時，需要動力學(xué)系統(tǒng)的傳遞矩陣和輸出響應(yīng)。仿真算例以兩端簡支梁為例，通過ANSYS有限元仿真，獲得相應(yīng)的傳遞矩陣和輸出響應(yīng)。該梁長度為1 m，密度為2 700 kg/m3，彈性模量為71 GPa，泊松比為0.3，結(jié)構(gòu)阻尼損耗因子為0.002，前3階固有頻率分別為118.95 Hz、473.8 Hz、1 058.8 Hz，仿真中載荷（F1、F2）和響應(yīng)參考點（1、2、3）的信息見表1，表中位置均表示與簡支梁左端的距離。

表1 載荷幅值、位置及響應(yīng)參考點位置

前兩次在不同位置施加單位激勵力，進(jìn)行諧響應(yīng)分析，獲得傳遞矩陣，而第三次施加一定大小的力，獲得系統(tǒng)加速度響應(yīng)。

考慮不確定性因素的影響，在仿真獲得的傳遞矩陣和加速度響應(yīng)中均加入信噪比為40 dB的高斯白噪聲，噪聲模型為

式中Nr取分布為N(0，A(ω)10-(B/20))的正態(tài)分布隨機(jī)數(shù)，A(ω)表示響應(yīng)或傳遞函數(shù)，B表示信噪比，Nu取0到1上的均勻分布隨機(jī)數(shù)。

利用添加噪聲后的傳遞矩陣和加速度響應(yīng)，假設(shè)誤差參數(shù)的先驗分布為IG（3,1），分別通過Gibbs抽樣算法和基于奇異值分解的載荷識別方法[12]，對第三次仿真施加的載荷進(jìn)行識別，結(jié)果如圖2所示。

以載荷頻率1 200 HZ時為例，使用Gibbs抽樣算法所得載荷的抽樣曲線如圖3所示，樣本直方圖如圖4所示。

對比圖2的結(jié)果可知，在大部分頻段內(nèi)，兩種方法的識別結(jié)果均較好，而在梁的固有頻率附近，通過Gibbs抽樣算法獲得的識別結(jié)果較奇異值分解法的識別結(jié)果更接近真實值，說明該方法能夠降低不確定因素對識別結(jié)果的影響。

將梁模型的邊界條件改為兩端固支，其余參數(shù)均不變，使用Gibbs抽樣算法所得載荷識別結(jié)果如圖5所示。

圖2 簡支梁所受載荷的識別結(jié)果

圖3 1 200 Hz時簡支梁所受載荷抽樣曲線

圖4 1 200 Hz時簡支梁所受載荷樣本直方圖

圖5 固支梁所受載荷識別結(jié)果

結(jié)果表明，識別精度與兩端簡支邊界條件下類似，在大部分頻段內(nèi)，載荷識別結(jié)果均較好，僅在固有頻率734.41 Hz和1 425.6 Hz附近具有較大誤差。因此，在固支邊界條件下，本文方法同樣適用。

3.2 實驗驗證

由于在工程應(yīng)用中，系統(tǒng)結(jié)構(gòu)復(fù)雜，且存在加工誤差、安裝誤差等多種因素的影響，為了驗證該方法在工程應(yīng)用中的適用性，故實驗對象采用結(jié)構(gòu)相對復(fù)雜的鋼制夾層梁，且安裝條件為兩端固支。該梁長0.468 m，寬0.03 m，高0.024 m，實驗臺架見圖6，從左至右，加速度傳感器分別安裝在響應(yīng)參考點1至6。

使用激振器施加一定大小的激勵力，穩(wěn)定后，通過LMS數(shù)采系統(tǒng)采集激勵力信號和加速度響應(yīng)信號。

圖6 實驗臺架

選用參考點1、3、6作為本次載荷辨識的響應(yīng)參考點，建立該結(jié)構(gòu)的有限元模型，通過仿真所得各響應(yīng)參考點的傳遞函數(shù)如圖7所示。

圖7 響應(yīng)參考點傳遞函數(shù)

利用仿真所得的傳遞函數(shù)和測量所得各響應(yīng)參考點的加速度響應(yīng)，通過Gibbs抽樣算法計算結(jié)構(gòu)所受的激勵力，計算所得的激勵力和測量所得的激勵力如圖8所示。

圖8 激勵力實測和計算結(jié)果

從圖8可以看出，基于貝葉斯估計理論，利用Gibbs抽樣算法進(jìn)行載荷識別，識別結(jié)果與實測結(jié)果基本相符，能夠滿足工程實際中的識別精度要求。

4 結(jié)語

本文從系統(tǒng)運(yùn)動方程出發(fā)，考慮不確定性因素的影響，對載荷和誤差參數(shù)服從的先驗概率分布進(jìn)行了假設(shè)，從而推導(dǎo)了載荷和誤差參數(shù)的后驗聯(lián)合概率分布，并分別得到二者的后驗邊緣概率分布。基于二者的后驗邊緣概率分布，通過馬爾科夫蒙特卡羅方法，采用Gibbs抽樣算法對系統(tǒng)所受載荷進(jìn)行了估計，并通過仿真算例和實驗數(shù)據(jù)對本方法進(jìn)行了驗證。結(jié)果表明，該方法滿足工程實際的要求，且在一定程度上降低了不確定性因素對載荷識別的影響，對提高載荷識別精度具有參考意義。

[1]胡寅寅，率志君，李玩幽，等.設(shè)備載荷識別與激勵源特性的研究現(xiàn)狀[J].噪聲與振動控制，2011，31(4)：1-5.

[2]THITE A N,THOMPSON D J.The quantification of structure-borne transmission paths by inverse methods.Part 1:Improved singular value rejection methods[J].Journal of Sound&Vibration,2003,264(2):411-431.

[3]CHOI H G,THITE A N,THOMPSON D J.Comparison of methods for parameter selection in Tikhonov regularization with application to inverse force determination[J].Journal of Sound&Vibration,2007,304(3-5):894-917.

[4]楊立娟，羅新杰，苗慧慧，等.混合傳遞路徑分析方法[J].噪聲與振動控制，2016，36(6)：12-15.

[5]張方，朱德懋.動載荷識別的時間有限元模型理論及其應(yīng)用[J].振動與沖擊，1998，17(2)：1-4.

[6]YAN G,ZHOU L.Impact load identification of composite structure using genetic algorithms[J].Journal of Sound&Vibration,2009,319(3-5):869-884.

[7]沙瑞華.基于神經(jīng)網(wǎng)絡(luò)的水電機(jī)組動載識別研究[D].大連：大連理工大學(xué)，2005.

[8]BOLSTAD W M.Introduction to Bayesian statistics,second edition[M].2007.

[9]BERNARDO J M,SMITH A F M.Bayesian theory[M].2008.

[10]曹晨.基于MCMC方法的統(tǒng)計模型的參數(shù)估計[D].南京：南京航空航天大學(xué)，2007.

[11]劉書奎，吳子燕，張玉兵，等.基于Gibbs抽樣的馬爾科夫蒙特卡羅方法在結(jié)構(gòu)物理參數(shù)識別及損傷定位中的研究[J].振動與沖擊，2011，30(10)：203-207.

[12]郭榮，房懷慶，裘剡，等.基于Tikhonov正則化及奇異值分解的載荷識別方法[J].振動與沖擊，2014，33(6)：53-58.