基于貝葉斯估計(jì)的動(dòng)載荷識(shí)別方法

雷 博，鄭 輝，2

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

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

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

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

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

1.1 貝葉斯估計(jì)理論[8–9]

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

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

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

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

1.2 概率分布估計(jì)

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

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

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

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

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

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

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

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

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

2 載荷抽樣估計(jì)

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

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

2.2 Gibbs抽樣算法[11]

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

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

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

圖1 Gibbs抽樣算法流程圖

3 仿真和實(shí)驗(yàn)驗(yàn)證

3.1 仿真算例

在基于本文方法進(jìn)行載荷識(shí)別時(shí)，需要?jiǎng)恿W(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)參考點(diǎn)（1、2、3）的信息見表1，表中位置均表示與簡支梁左端的距離。

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

前兩次在不同位置施加單位激勵(lì)力，進(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ù)的先驗(yàn)分布為IG（3,1），分別通過Gibbs抽樣算法和基于奇異值分解的載荷識(shí)別方法[12]，對(duì)第三次仿真施加的載荷進(jìn)行識(shí)別，結(jié)果如圖2所示。

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

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

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

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

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

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

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

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

3.2 實(shí)驗(yàn)驗(yàn)證

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

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

圖6 實(shí)驗(yàn)臺(tái)架

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

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

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

圖8 激勵(lì)力實(shí)測(cè)和計(jì)算結(jié)果

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

4 結(jié)語

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

[1]胡寅寅，率志君，李玩幽，等.設(shè)備載荷識(shí)別與激勵(lì)源特性的研究現(xiàn)狀[J].噪聲與振動(dòng)控制，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].噪聲與振動(dòng)控制，2016，36(6)：12-15.

[5]張方，朱德懋.動(dòng)載荷識(shí)別的時(shí)間有限元模型理論及其應(yīng)用[J].振動(dòng)與沖擊，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òng)載識(shí)別研究[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)計(jì)模型的參數(shù)估計(jì)[D].南京：南京航空航天大學(xué)，2007.

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

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