一种基于Bayes方法的随机模型修正方法

马天政, 吕　昊, 张义民

(东北大学 机械工程与自动化学院, 辽宁 沈阳 110819)



一种基于Bayes方法的随机模型修正方法

马天政, 吕昊, 张义民

(东北大学 机械工程与自动化学院, 辽宁 沈阳 110819)

提出了一种随机模型的修正方法用以估计结构参数的统计特性.基于Bayes方法的参数估计原理,将需要修正的结构参数的均值和方差看作符合一定先验概率分布的随机变量,根据核密度估计原理构建得到似然函数,进而使用基于差分进化的MCMC方法估计参数的后验概率密度，并根据最大后验概率密度准则估计结构参数的均值和方差.同时使用Kriging方法建立了结构输入和输出之间的代理模型,保证计算精度的同时极大地节约了计算时间.数值算例验证了本方法的可行性.

随机模型修正； Kriging代理模型； Bayes方法； MCMC抽样

近些年来,传统的确定性有限元模型动力学修正已经在汽车、航空航天等领域得到了广泛的应用,它以单个结构为研究对象,认为模型参数是确定性的值,通过结构的动力学响应间接地估计结构参数的取值,从而使得有限元计算结果和试验结果相一致.然而实际工程中的结构存在着大量的不确定性因素,如由于材料参数的变异性、加工制造过程中的误差、零部件连接处刚度的不确定性等因素导致结构的参数及动力学特性(如固有频率、固有振型、频响函数等)存在不确定性,这种不确定性可以通过概率统计的方法来描述.

随机有限元动力学修正以相同条件下制造出来的一批结构为研究对象,通过模态试验得到一批结构的动力学统计特性来间接估计结构参数的统计特性,从而使得数值计算得到的响应概率密度分布与试验结果的概率密度分布相一致,可以看作是随机有限元方法的逆问题.通过随机模型修正,可以对一些难以直接测量的结构参数(如阻尼、结构连接处的刚度)的统计特性进行识别.

在随机模型修正的研究中,诸多学者采用Monte Carlo方法结合响应面模型根据结构响应的前两阶矩来修正结构参数的均值和方差[1-3],这种方法需要较大的计算量.文献[4]在Monte Carlo方法的基础上,将随机模型修正问题转化为一系列的确定性模型,避免计算结构响应的统计矩,在一定程度上减少了计算量.文献[5]研究了摄动法在随机模型修正中的应用,这种方法提高了计算效率,但在使用的过程中需要计算灵敏度矩阵且推导过程较为繁复.文献[6]在此基础上进一步将参数的均值和方差的修正分解为2个相对独立的部分分别进行迭代计算,简化了计算过程.此外,区间分析[7-8]的方法也在随机模型修正中得到了一定的应用.

本文针对工程实际中常见的结构参数服从Gauss分布的情况,提出了一种基于Bayes原理的随机模型修正方法来估计结构参数的均值和方差.在Bayes参数估计过程中,根据核密度估计原理计算结构的似然函数,并使用基于差分进化的MCMC方法计算参数的后验概率密度分布以进行统计推断.为了节约计算时间,建立了Kriging代理模型用来拟合结构输入参数和输出之间的关系.数值算例表明本方法具有良好的适用性.

1　Bayes方法

Bayes方法在进行参数估计时,将待估计的参数看作符合一定先验分布的随机变量,根据系统实际测量得到的动态响应通过Bayes公式计算得到结构参数的后验分布,进而对结构参数进行估计.Bayes方法不仅能给出参数的最优估计,而且能得到参数取值的概率密度分布,此外Bayes方法能够有效地将以往的工程经验(通过先验分布)和结构的实际观测结果融合在一起进行参数估计.

(1)

计算在参数的每个样本点上结构的前Nfreq阶固有频率：

其中：

因此有

(2)

将式(2)代入式(1)得到参数的后验概率密度分布表达式为

(3)

获得了参数的后验概率密度分布后,即可根据最大后验概率密度准则(maximumaposterioriprobabilityestimate,MAP)进行统计推断,即选择概率密度取值最大的点作为参数的估计值.

2　基于差分进化的MCMC抽样

使用Bayes方法时一个很重要的环节在于如何计算参数的后验分布,通常情况下难以获得后验概率密度分布的显示表达式,马尔可夫蒙特卡洛模拟(MCMC)方法是解决这一问题的有效途径.MCMC方法通过随机抽样构建一条极限状态分布为后验分布(3)的马尔可夫链,当马尔可夫链达到平稳状态时从中抽取一定数量的样本,样本的分布即为后验分布(3)的近似.Gibbs抽样和Metropolis-Hastings(MH)抽样式是常用的MCMC抽样方法,然而当参数的后验分布极其复杂时(如结构动力学问题),传统的MCMC抽样方法容易陷入局部最优,难以准确地计算后验概率密度,为此基于进化计算(如遗传算法、差分进化等)的MCMC抽样方法近些年来得到了广泛的研究,本文采用基于差分进化的MCMC方法来计算后验概率[9-10].

其中Vs称为第s条马尔可夫链的“温度”.可以看出,第1条链所对应的极限状态分布即为需要计算的参数后验概率密度分布.由于Vs是介于0和1之间的常数,因此其余Nc-1条链所对应的平稳分布比参数的后验概率密度分布要光滑,Vs越小，就越光滑,同时也意味着更容易进行抽样.在整个抽样计算过程中,一方面Nc条马氏链各自进行MH抽样,另外一方面不同的链条之间通过定义的变异、交换以及差分进化操作进行着信息交换,从而可以有效地计算参数的后验概率密度分布，避免陷入局部最优.

基于差分进化的MCMC抽样的计算过程为：

3)变异：选取一条马氏链进行另外一起MH抽样.

②对第k条马氏链进行MH抽样得到新的样本值Θ+;

5)交换：交换2条马尔可夫链的当前样本值.

③产生均匀分布的随机数u～U(0,1),若u

6)差分进化：来自于智能优化算法中的差分进化方法.

①等概率随机选取3条马氏链j,k,l(1≤j,k,l≤Nc),对第j条马氏链进行差分进化操作.

②产生变异的新样本Θ+：

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

其中γ是权重系数.

⑤产生均匀分布的随机数u～U(0,1),若u

3　Kriging代理模型

在Bayes方法中,使用基于差分进化的MCMC计算结构参数的后验分布时,需要大量计算结构在不同参数下的固有频率,若直接采用结构的有限元模型进行计算显然是不可行的,为了减少计算量、节约计算时间，需要建立原结构模型的代理模型(Metamodel).本文采用Kriging方法拟合输入参数和输出响应之间的关系,和其它方法(如响应面模型、支持向量机、神经网络等)相比,Kriging方法有更好的精度和计算效率,对于非线性模型和具有局部响应突变的模型有着较好的拟合效果[11].

Kriging模型由一个回归模型和一个随机误差模型构成：

其中，xj和ωj分别为变量x和ω的第j个分量,d为设计变量空间的维数.

(4)

(5)

其中，F为系数矩阵,Y为结构响应矩阵,R为相关矩阵.相关参数ρ的最优值可以通过极大似然法估计得到,即

根据式(4)和式(5)求得的Kriging模型参数,对于在未知设计点x*的结构响应y*为

在拟合Kriging模型时,通过拉丁超立方抽样获得样本点集X,其对应的结构响应Y通过结构的有限元模型计算得到.

4　数值算例

三自由度弹簧质量块系统如图1所示.结构的确定性参数为质量m1,m2,m3和弹簧刚度k3,k4,k6，结构的不确定性参数为弹簧刚度k1,k2,k5.

图1　弹簧-质量块系统Fig.1　Spring-mass system

m1=m2=m3=1.0 kg，

k3=k4=1.0 N/m，

k6=3.0 N/m.

结构的不确定性参数为正态分布的随机变量,其分布特性为:

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

使用拉丁超立方抽样抽取100个样本点拟合Kriging模型.为了检验代理模型精度,随机抽取20个样本，分别使用有限元方法和代理模型计算结构的固有频率、最大误差、最小误差及平均误差，使用代理模型的计算结果如表1所示.

表1　Kriging代理模型精度检验

在模型修正过程中,选取参数取值区间上的均匀分布作为参数的无信息先验分布.使用基于差分进化的MCMC方法抽取一条长度为20 000的马尔可夫链,舍掉前5 000个样本作为burn-in阶段,剩下的15 000个样本用来估计参数的后验概率密度,如图2、图3所示.参数均值和标准差的初始估计值和修正值如表2所示,从中可以看出,参数的均值在修正后能够较好地收敛到真实值(最大误差为3.2%),但参数的标准差在修正后与真实值的误差相对较大(最大误差为17.7%),这主要是由于样本数量较少造成的,增加样本的数量会提高标准差的修正值误差.需要说明的是,目前的随机模型修正方法对方差的修正都存在着缺陷[12].

表2　模型修正结果

使用Monte Carlo方法进行500次抽样绘制了修正后的频率分布图,如图4所示,在修正之后结构的频率分布和测量的样本频率分布能够很好地吻合.

图2　MCMC模拟样本Fig.2　Samples of MCMC simulation

图3　结构参数的后验概率密度分布Fig.3　Posterior probability density distribution of the structural parameters

图4　修正后的频率分布云图Fig.4　Updated frequency distribution clouds

5　结　论

本文结合Bayes参数估计原理、核密度估计、基于差分进化的MCMC抽样方法和Kriging代理模型提出了一种随机模型的修正方法,数值算例验证了方法的有效性和可行性.同时,本文提出的方法对于参数的变异性没有要求,可以适用于参数变异性较大的情况；而且不需要计算结构参数的灵敏度矩阵和结构响应的均值、方差等统计矩,易于推广到接触、冲击等具有强非线性的问题中.另外一方面,在实际工程应用中,使用Bayes方法能够充分地将以往的实际经验和结构的观测数据融合在一起.

[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模型的涡轮盘优化设计方法[J].计算机集成制造系统,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] 方圣恩,林友勤,夏樟华.考虑结构参数不确定性的随机模型修正方法[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)；国家重点基础研究发展计划(973计划)项目(2014CB046303)；中央高校基本科研业务费资助项目(02090022115014);“高档数控机床与基础制造装备”科技重大专项课题(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

本刊网址·在线期刊：http://www.journals.zju.edu.cn/gcsjxb