馬天政, 呂　昊, 張義民

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

?

一種基于Bayes方法的隨機模型修正方法

馬天政, 呂昊, 張義民

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

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

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

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

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

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

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

1　Bayes方法

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

(1)

計算在參數(shù)的每個樣本點上結構的前Nfreq階固有頻率：

其中：

因此有

(2)

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

(3)

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

2　基于差分進化的MCMC抽樣

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

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

基于差分進化的MCMC抽樣的計算過程為：

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

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

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

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

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

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

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

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

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

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

3　Kriging代理模型

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

Kriging模型由一個回歸模型和一個隨機誤差模型構成：

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

(4)

(5)

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

根據(jù)式(4)和式(5)求得的Kriging模型參數(shù),對于在未知設計點x*的結構響應y*為

在擬合Kriging模型時,通過拉丁超立方抽樣獲得樣本點集X,其對應的結構響應Y通過結構的有限元模型計算得到.

4　數(shù)值算例

三自由度彈簧質(zhì)量塊系統(tǒng)如圖1所示.結構的確定性參數(shù)為質(zhì)量m1,m2,m3和彈簧剛度k3,k4,k6，結構的不確定性參數(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.

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

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

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

表1　Kriging代理模型精度檢驗

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

表2　模型修正結果

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

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

圖3　結構參數(shù)的后驗概率密度分布Fig.3　Posterior probability density distribution of the structural parameters

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

5　結　論

本文結合Bayes參數(shù)估計原理、核密度估計、基于差分進化的MCMC抽樣方法和Kriging代理模型提出了一種隨機模型的修正方法,數(shù)值算例驗證了方法的有效性和可行性.同時,本文提出的方法對于參數(shù)的變異性沒有要求,可以適用于參數(shù)變異性較大的情況；而且不需要計算結構參數(shù)的靈敏度矩陣和結構響應的均值、方差等統(tǒng)計矩,易于推廣到接觸、沖擊等具有強非線性的問題中.另外一方面,在實際工程應用中,使用Bayes方法能夠充分地將以往的實際經(jīng)驗和結構的觀測數(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法的結構動力學模型確認[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)化設計方法[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] 方圣恩,林友勤,夏樟華.考慮結構參數(shù)不確定性的隨機模型修正方法[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.

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

馬天政(1987—),男,遼寧鞍山人,博士,從事隨機模型修正及動力學可靠性研究,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