馬天政, 呂　昊, 張義民

(東北大學(xué) 機(jī)械工程與自動化學(xué)院, 遼寧 沈陽 110819)

?

一種基于Bayes方法的隨機(jī)模型修正方法

馬天政, 呂昊, 張義民

(東北大學(xué) 機(jī)械工程與自動化學(xué)院, 遼寧 沈陽 110819)

提出了一種隨機(jī)模型的修正方法用以估計(jì)結(jié)構(gòu)參數(shù)的統(tǒng)計(jì)特性.基于Bayes方法的參數(shù)估計(jì)原理,將需要修正的結(jié)構(gòu)參數(shù)的均值和方差看作符合一定先驗(yàn)概率分布的隨機(jī)變量,根據(jù)核密度估計(jì)原理構(gòu)建得到似然函數(shù),進(jìn)而使用基于差分進(jìn)化的MCMC方法估計(jì)參數(shù)的后驗(yàn)概率密度，并根據(jù)最大后驗(yàn)概率密度準(zhǔn)則估計(jì)結(jié)構(gòu)參數(shù)的均值和方差.同時使用Kriging方法建立了結(jié)構(gòu)輸入和輸出之間的代理模型,保證計(jì)算精度的同時極大地節(jié)約了計(jì)算時間.數(shù)值算例驗(yàn)證了本方法的可行性.

隨機(jī)模型修正； Kriging代理模型； Bayes方法； MCMC抽樣

近些年來,傳統(tǒng)的確定性有限元模型動力學(xué)修正已經(jīng)在汽車、航空航天等領(lǐng)域得到了廣泛的應(yīng)用,它以單個結(jié)構(gòu)為研究對象,認(rèn)為模型參數(shù)是確定性的值,通過結(jié)構(gòu)的動力學(xué)響應(yīng)間接地估計(jì)結(jié)構(gòu)參數(shù)的取值,從而使得有限元計(jì)算結(jié)果和試驗(yàn)結(jié)果相一致.然而實(shí)際工程中的結(jié)構(gòu)存在著大量的不確定性因素,如由于材料參數(shù)的變異性、加工制造過程中的誤差、零部件連接處剛度的不確定性等因素導(dǎo)致結(jié)構(gòu)的參數(shù)及動力學(xué)特性(如固有頻率、固有振型、頻響函數(shù)等)存在不確定性,這種不確定性可以通過概率統(tǒng)計(jì)的方法來描述.

隨機(jī)有限元動力學(xué)修正以相同條件下制造出來的一批結(jié)構(gòu)為研究對象,通過模態(tài)試驗(yàn)得到一批結(jié)構(gòu)的動力學(xué)統(tǒng)計(jì)特性來間接估計(jì)結(jié)構(gòu)參數(shù)的統(tǒng)計(jì)特性,從而使得數(shù)值計(jì)算得到的響應(yīng)概率密度分布與試驗(yàn)結(jié)果的概率密度分布相一致,可以看作是隨機(jī)有限元方法的逆問題.通過隨機(jī)模型修正,可以對一些難以直接測量的結(jié)構(gòu)參數(shù)(如阻尼、結(jié)構(gòu)連接處的剛度)的統(tǒng)計(jì)特性進(jìn)行識別.

在隨機(jī)模型修正的研究中,諸多學(xué)者采用Monte Carlo方法結(jié)合響應(yīng)面模型根據(jù)結(jié)構(gòu)響應(yīng)的前兩階矩來修正結(jié)構(gòu)參數(shù)的均值和方差[1-3],這種方法需要較大的計(jì)算量.文獻(xiàn)[4]在Monte Carlo方法的基礎(chǔ)上,將隨機(jī)模型修正問題轉(zhuǎn)化為一系列的確定性模型,避免計(jì)算結(jié)構(gòu)響應(yīng)的統(tǒng)計(jì)矩,在一定程度上減少了計(jì)算量.文獻(xiàn)[5]研究了攝動法在隨機(jī)模型修正中的應(yīng)用,這種方法提高了計(jì)算效率,但在使用的過程中需要計(jì)算靈敏度矩陣且推導(dǎo)過程較為繁復(fù).文獻(xiàn)[6]在此基礎(chǔ)上進(jìn)一步將參數(shù)的均值和方差的修正分解為2個相對獨(dú)立的部分分別進(jìn)行迭代計(jì)算,簡化了計(jì)算過程.此外,區(qū)間分析[7-8]的方法也在隨機(jī)模型修正中得到了一定的應(yīng)用.

本文針對工程實(shí)際中常見的結(jié)構(gòu)參數(shù)服從Gauss分布的情況,提出了一種基于Bayes原理的隨機(jī)模型修正方法來估計(jì)結(jié)構(gòu)參數(shù)的均值和方差.在Bayes參數(shù)估計(jì)過程中,根據(jù)核密度估計(jì)原理計(jì)算結(jié)構(gòu)的似然函數(shù),并使用基于差分進(jìn)化的MCMC方法計(jì)算參數(shù)的后驗(yàn)概率密度分布以進(jìn)行統(tǒng)計(jì)推斷.為了節(jié)約計(jì)算時間,建立了Kriging代理模型用來擬合結(jié)構(gòu)輸入?yún)?shù)和輸出之間的關(guān)系.數(shù)值算例表明本方法具有良好的適用性.

1　Bayes方法

Bayes方法在進(jìn)行參數(shù)估計(jì)時,將待估計(jì)的參數(shù)看作符合一定先驗(yàn)分布的隨機(jī)變量,根據(jù)系統(tǒng)實(shí)際測量得到的動態(tài)響應(yīng)通過Bayes公式計(jì)算得到結(jié)構(gòu)參數(shù)的后驗(yàn)分布,進(jìn)而對結(jié)構(gòu)參數(shù)進(jìn)行估計(jì).Bayes方法不僅能給出參數(shù)的最優(yōu)估計(jì),而且能得到參數(shù)取值的概率密度分布,此外Bayes方法能夠有效地將以往的工程經(jīng)驗(yàn)(通過先驗(yàn)分布)和結(jié)構(gòu)的實(shí)際觀測結(jié)果融合在一起進(jìn)行參數(shù)估計(jì).

(1)

計(jì)算在參數(shù)的每個樣本點(diǎn)上結(jié)構(gòu)的前Nfreq階固有頻率：

其中：

因此有

(2)

將式(2)代入式(1)得到參數(shù)的后驗(yàn)概率密度分布表達(dá)式為

(3)

獲得了參數(shù)的后驗(yàn)概率密度分布后,即可根據(jù)最大后驗(yàn)概率密度準(zhǔn)則(maximumaposterioriprobabilityestimate,MAP)進(jìn)行統(tǒng)計(jì)推斷,即選擇概率密度取值最大的點(diǎn)作為參數(shù)的估計(jì)值.

2　基于差分進(jìn)化的MCMC抽樣

使用Bayes方法時一個很重要的環(huán)節(jié)在于如何計(jì)算參數(shù)的后驗(yàn)分布,通常情況下難以獲得后驗(yàn)概率密度分布的顯示表達(dá)式,馬爾可夫蒙特卡洛模擬(MCMC)方法是解決這一問題的有效途徑.MCMC方法通過隨機(jī)抽樣構(gòu)建一條極限狀態(tài)分布為后驗(yàn)分布(3)的馬爾可夫鏈,當(dāng)馬爾可夫鏈達(dá)到平穩(wěn)狀態(tài)時從中抽取一定數(shù)量的樣本,樣本的分布即為后驗(yàn)分布(3)的近似.Gibbs抽樣和Metropolis-Hastings(MH)抽樣式是常用的MCMC抽樣方法,然而當(dāng)參數(shù)的后驗(yàn)分布極其復(fù)雜時(如結(jié)構(gòu)動力學(xué)問題),傳統(tǒng)的MCMC抽樣方法容易陷入局部最優(yōu),難以準(zhǔn)確地計(jì)算后驗(yàn)概率密度,為此基于進(jìn)化計(jì)算(如遺傳算法、差分進(jìn)化等)的MCMC抽樣方法近些年來得到了廣泛的研究,本文采用基于差分進(jìn)化的MCMC方法來計(jì)算后驗(yàn)概率[9-10].

其中Vs稱為第s條馬爾可夫鏈的“溫度”.可以看出,第1條鏈所對應(yīng)的極限狀態(tài)分布即為需要計(jì)算的參數(shù)后驗(yàn)概率密度分布.由于Vs是介于0和1之間的常數(shù),因此其余Nc-1條鏈所對應(yīng)的平穩(wěn)分布比參數(shù)的后驗(yàn)概率密度分布要光滑,Vs越小，就越光滑,同時也意味著更容易進(jìn)行抽樣.在整個抽樣計(jì)算過程中,一方面Nc條馬氏鏈各自進(jìn)行MH抽樣,另外一方面不同的鏈條之間通過定義的變異、交換以及差分進(jìn)化操作進(jìn)行著信息交換,從而可以有效地計(jì)算參數(shù)的后驗(yàn)概率密度分布，避免陷入局部最優(yōu).

基于差分進(jìn)化的MCMC抽樣的計(jì)算過程為：

3)變異：選取一條馬氏鏈進(jìn)行另外一起MH抽樣.

②對第k條馬氏鏈進(jìn)行MH抽樣得到新的樣本值Θ+;

5)交換：交換2條馬爾可夫鏈的當(dāng)前樣本值.

③產(chǎn)生均勻分布的隨機(jī)數(shù)u～U(0,1),若u

6)差分進(jìn)化：來自于智能優(yōu)化算法中的差分進(jìn)化方法.

①等概率隨機(jī)選取3條馬氏鏈j,k,l(1≤j,k,l≤Nc),對第j條馬氏鏈進(jìn)行差分進(jìn)化操作.

②產(chǎn)生變異的新樣本Θ+：

Θ+=Θ(j)(i)+γ(Θ(k)(i)-Θ(l)(i))，

其中γ是權(quán)重系數(shù).

⑤產(chǎn)生均勻分布的隨機(jī)數(shù)u～U(0,1),若u

3　Kriging代理模型

在Bayes方法中,使用基于差分進(jìn)化的MCMC計(jì)算結(jié)構(gòu)參數(shù)的后驗(yàn)分布時,需要大量計(jì)算結(jié)構(gòu)在不同參數(shù)下的固有頻率,若直接采用結(jié)構(gòu)的有限元模型進(jìn)行計(jì)算顯然是不可行的,為了減少計(jì)算量、節(jié)約計(jì)算時間，需要建立原結(jié)構(gòu)模型的代理模型(Metamodel).本文采用Kriging方法擬合輸入?yún)?shù)和輸出響應(yīng)之間的關(guān)系,和其它方法(如響應(yīng)面模型、支持向量機(jī)、神經(jīng)網(wǎng)絡(luò)等)相比,Kriging方法有更好的精度和計(jì)算效率,對于非線性模型和具有局部響應(yīng)突變的模型有著較好的擬合效果[11].

Kriging模型由一個回歸模型和一個隨機(jī)誤差模型構(gòu)成：

其中，xj和ωj分別為變量x和ω的第j個分量,d為設(shè)計(jì)變量空間的維數(shù).

(4)

(5)

其中，F(xiàn)為系數(shù)矩陣,Y為結(jié)構(gòu)響應(yīng)矩陣,R為相關(guān)矩陣.相關(guān)參數(shù)ρ的最優(yōu)值可以通過極大似然法估計(jì)得到,即

根據(jù)式(4)和式(5)求得的Kriging模型參數(shù),對于在未知設(shè)計(jì)點(diǎn)x*的結(jié)構(gòu)響應(yīng)y*為

在擬合Kriging模型時,通過拉丁超立方抽樣獲得樣本點(diǎn)集X,其對應(yīng)的結(jié)構(gòu)響應(yīng)Y通過結(jié)構(gòu)的有限元模型計(jì)算得到.

4　數(shù)值算例

三自由度彈簧質(zhì)量塊系統(tǒng)如圖1所示.結(jié)構(gòu)的確定性參數(shù)為質(zhì)量m1,m2,m3和彈簧剛度k3,k4,k6，結(jié)構(gòu)的不確定性參數(shù)為彈簧剛度k1,k2,k5.

圖1　彈簧-質(zhì)量塊系統(tǒng)Fig.1　Spring-mass system

m1=m2=m3=1.0 kg，

k3=k4=1.0 N/m，

k6=3.0 N/m.

結(jié)構(gòu)的不確定性參數(shù)為正態(tài)分布的隨機(jī)變量,其分布特性為:

σk1=σk2=σk5=0.2 N/m.

使用拉丁超立方抽樣抽取100個樣本點(diǎn)擬合Kriging模型.為了檢驗(yàn)代理模型精度,隨機(jī)抽取20個樣本，分別使用有限元方法和代理模型計(jì)算結(jié)構(gòu)的固有頻率、最大誤差、最小誤差及平均誤差，使用代理模型的計(jì)算結(jié)果如表1所示.

表1　Kriging代理模型精度檢驗(yàn)

在模型修正過程中,選取參數(shù)取值區(qū)間上的均勻分布作為參數(shù)的無信息先驗(yàn)分布.使用基于差分進(jìn)化的MCMC方法抽取一條長度為20 000的馬爾可夫鏈,舍掉前5 000個樣本作為burn-in階段,剩下的15 000個樣本用來估計(jì)參數(shù)的后驗(yàn)概率密度,如圖2、圖3所示.參數(shù)均值和標(biāo)準(zhǔn)差的初始估計(jì)值和修正值如表2所示,從中可以看出,參數(shù)的均值在修正后能夠較好地收斂到真實(shí)值(最大誤差為3.2%),但參數(shù)的標(biāo)準(zhǔn)差在修正后與真實(shí)值的誤差相對較大(最大誤差為17.7%),這主要是由于樣本數(shù)量較少造成的,增加樣本的數(shù)量會提高標(biāo)準(zhǔn)差的修正值誤差.需要說明的是,目前的隨機(jī)模型修正方法對方差的修正都存在著缺陷[12].

表2　模型修正結(jié)果

使用Monte Carlo方法進(jìn)行500次抽樣繪制了修正后的頻率分布圖,如圖4所示,在修正之后結(jié)構(gòu)的頻率分布和測量的樣本頻率分布能夠很好地吻合.

圖2　MCMC模擬樣本Fig.2　Samples of MCMC simulation

圖3　結(jié)構(gòu)參數(shù)的后驗(yàn)概率密度分布Fig.3　Posterior probability density distribution of the structural parameters

圖4　修正后的頻率分布云圖Fig.4　Updated frequency distribution clouds

5　結(jié)　論

本文結(jié)合Bayes參數(shù)估計(jì)原理、核密度估計(jì)、基于差分進(jìn)化的MCMC抽樣方法和Kriging代理模型提出了一種隨機(jī)模型的修正方法,數(shù)值算例驗(yàn)證了方法的有效性和可行性.同時,本文提出的方法對于參數(shù)的變異性沒有要求,可以適用于參數(shù)變異性較大的情況；而且不需要計(jì)算結(jié)構(gòu)參數(shù)的靈敏度矩陣和結(jié)構(gòu)響應(yīng)的均值、方差等統(tǒng)計(jì)矩,易于推廣到接觸、沖擊等具有強(qiáng)非線性的問題中.另外一方面,在實(shí)際工程應(yīng)用中,使用Bayes方法能夠充分地將以往的實(shí)際經(jīng)驗(yàn)和結(jié)構(gòu)的觀測數(shù)據(jù)融合在一起.

[1] BAO Nuo,WANG Chun-jie.A Monte Carlo simulation based inverse propagation method for stochastic model updating[J].Mechanical Systems and Signal Processing,2015(60/61):928-944．

[2] MOTTERSHEAD J E,MARES C,FRISWELL M I.Stochastic model updating:part 1：theory and simulated example[J].Mechanical Systems and Signal Processing,2006,20(7):1674-1695．

[3] 陳志國,鄧忠民,畢司峰.基于Monte Carlo法的結(jié)構(gòu)動力學(xué)模型確認(rèn)[J].振動與沖擊,2013,32(16):76-81．

CHEN Zhi-guo,DENG Zhong-min,BI Si-feng.Structural dynamics model validation based on Monte Carlo method[J].Journal of vibration and shock,2013,32(16):76-81．

[4] FANG Shen-gen,REN Wei-xin,PERERA R.A stochastic model updating method for parameter variability quantification based on response surface models and Monte Carlo simulation[J].Mechanical Systems and Signal Processing,2012,33:83-96．

[5] KHODAPARAST H H,MOTTERSHEAD J E,FRISWELL M I.Perturbation methods for the estimation of parameter variability in stochastic model updating[J].Mechanical Systems and Signal Processing,2008,22(8):1751-1773．

[6] GOVERS Y,LINK M.Stochastic model updating-Covariance matrix adjustment from uncertain experimental modal data[J].Mechanical Systems and Signal Processing,2010,24(3):696-706．

[7] KHODAPARAST H H,MOTTERSHEAD J E,BADCOCK K J.Interval model updating with irreducible uncertainty using the Kriging predictor[J].Mechanical Systems and Signal Processing,2011,25(4):1204-1226．

[8] FANG Shen-gen,ZHANG Qiu-hu,REN Wei-xin.An interval model updating strategy using interval response surface models[J].Mechanical Systems and Signal Processing,2015(60/61):909-927．

[9] JASRA A,STEPHENS D A,HOLMES C C.On population-based simulation for static inference [J].Statistics and Computing,2015(60/61):909-927．

[10] NICHOLS J M,MOORE E Z,MURPHY K D.Bayesian identification of a cracked plate using a population-based Markov chain Monte Carlo method[J].Computers & Structures,2011,89(13/14):1323-1332．

[11] 黃章俊,王成恩.基于Kriging模型的渦輪盤優(yōu)化設(shè)計(jì)方法[J].計(jì)算機(jī)集成制造系統(tǒng),2010,16(5):905-911．

HUANG Zhang-jun,WANG Cheng-en.Turbine discs optimization design based on Kriging model[J].Computer Integrated Manufacturing Systems,2010,16(5):905-911．

[12] 方圣恩,林友勤,夏樟華.考慮結(jié)構(gòu)參數(shù)不確定性的隨機(jī)模型修正方法[J].振動、測試與診斷,2014,34(5)：832-837．

FANG Sheng-en,LIN You-qin,XIA Zhang-hua.Stochastic model updating method considering the uncertainties of strucutral parameters[J].Journal of Vibration,Measurement & Diagnosis,2014,34(5):832-837．

Stochastic model updating based on Bayesian method

MA Tian-zheng, Lü Hao, ZHANG Yi-min

(School of Mechanical Engineering & Automation，Northeastern University, Shenyang 110819, China)

A new method for stochastic model updating was proposed to estimate the statistical information of the structural parameters.According to the principle of Bayesian method,the parameters’ mean value and variance to be estimated were regarded as random variables and the likelihood functions were constructed by using kernel density estimation method.The posterior distribution of the parameters were calculated by utilizing the population-based MCMC simulation method and then the parameters’ mean value and variance could be obtained based on MAP (maximum a posterior) principle.A surrogate model based on Kriging method was established,which greatly saved the computational cost.Numerical example demonstrated the effectiveness of the method．

stochastic model updating; Kriging agent model; Bayesian method; MCMC sampling

2016-01-06.

國家自然科學(xué)基金重點(diǎn)資助項(xiàng)目(51135003)；國家自然科學(xué)基金資助項(xiàng)目(U1234208)；國家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃(973計(jì)劃)項(xiàng)目(2014CB046303)；中央高?；究蒲袠I(yè)務(wù)費(fèi)資助項(xiàng)目(02090022115014);“高檔數(shù)控機(jī)床與基礎(chǔ)制造裝備”科技重大專項(xiàng)課題(2013ZX04011011).

馬天政(1987—),男,遼寧鞍山人,博士,從事隨機(jī)模型修正及動力學(xué)可靠性研究,E-mail:zpaprecv@sohu.com.http://orcid.org//0000-0002-9509-9802

10.3785/j.issn. 1006-754X.2016.03.002

O 327; TU 311

A

1006-754X(2016)03-0206-06

本刊網(wǎng)址·在線期刊：http://www.journals.zju.edu.cn/gcsjxb