基于Kriging模型和改进MCMC算法的随机有限元模型修正

张雪萍, 彭珍瑞, 张亚峰

(兰州交通大学 机电工程学院，兰州 730070)

1 引 言

在工程领域中，一个能够准确表征结构行为的精确有限元模型非常重要[1]。有限元模型修正已成为振动领域研究的热点之一[2]。但在实际工程应用中，由于结构几何尺寸、材料变异性和模型误差等因素的影响，不确定性广泛存在于实际结构的参数和响应中[3]。

为了使动力学模型能够尽可能地反映实际情况，国内外开始运用随机系统理论建立不确定性模型，在有限元模型修正过程中结合概率方法进行分析是十分必要的[4]。不确定性模型修正方法中基于贝叶斯统计理论的研究是一个热点。Beck等[5]将贝叶斯理论引入到模型修正中，并明确了修正的基本思路。韩芳等[6]提出了一种基于信息融合和贝叶斯理论的模型修正方法。但是，贝叶斯方法存在无法直接积分求取参数的似然函数和后验概率公式的积分常数等缺点。因此，贝叶斯方法的参数后验概率密度常用马尔科夫链蒙特卡罗MCMC(Markov Chain Monte Carlo)方法来推断。Beck等[7]又提出了一种基于MH(Metropolis -Hastings)算法的自适应MCMC方法。Green[8]提出了融合改进退火算法的MCMC方法，避免了目标函数陷入局部最优。Lam等[9]使用多个CPU生成多个马尔可夫链(而不是单个链)的并行MCMC方法来提高MCMC方法的计算效率。蒋伟等[10]在标准MCMC方法的基础上，引入差分进化算法，为提高不确定性模型修正中的计算精度提供了一种新手段。然而，在参数维度较高的情况下，虽然MCMC算法采样所得的样本具有一定的收敛性，但马尔可夫链出现大量的停滞阶段，样本的接受率很低[11]。彭珍瑞等[12]引入了布谷鸟算法中新鸟巢更新的思想改进MCMC算法。以上文献大多采用结构的响应为频率和位移等进行模型修正，而在工程实践中会出现选择的响应对局部状态参数不敏感的情况，从而在修正过程中不能得到准确、合理的应力分布。因此，在有限元模型修正方法中，选择一种能够同时反映结构局部状态与全局特征的响应至关重要。同时，在有限元模型修正中，由于所用代理模型精度不高也会使修正结果不理想。

基于上述背景，本文引入花朵授粉算法FPA(Flower Pollination Algorithm)与标准MCMC算法融合，解决在有效提高新样本接受率的同时降低修正参数的相对误差；且由于应变模态包含结构全局的频率信息和能够表征结构局部状态的应变振型信息，因此选用应变模态作为结构响应进行模型修正，并使用蝙蝠算法BA(Bat Algorithm)确定Kriging模型相关系数，建立精确的代理模型，提高结构响应的计算效率。最后，利用三自由度质量-弹簧系统和三维桁架结构两个数值算例验证所提方法的有效性。

2 基于应变模态的贝叶斯方法

2.1 应变模态

应变是由位移经过一次求导得来的，对于结构整体而言，应变与位移之间的关系可表示为[13]

ε=BTφs

(1)

式中ε为结构所有点的应变,B为结构整体几何矩阵,T为局部坐标和总体坐标之间的转换矩阵,φs=φYφTFej ω t为总体坐标的节点位移向量。

根据有限元分析模型，由模态叠加法可得

ε=BTφYφTFejω t=ΨεYφTFejω t=

(2)

2.2 贝叶斯方法

基于贝叶斯统计理论实现有限元模型修正，依据参数θ的先验信息结合观测数据D(模态信息)的似然函数来估计参数的后验概率分布。其简化公式可表示为[5]

P(θ|D)=a·P(D|θ)P(θ)

(3)

假设模型误差向量为e，则在模型修正中，对于待修正参数向量θ得到以下关系，

(4，5)

(6，7)

P(D|θ)=P(ω|θ)·P(Ψε|θ)

(8)

(9)

(10)

为减少可行解主观成分，在贝叶斯方法中引入熵。对模型修正这个反问题，应充分利用最大熵原理作为先验信息和约束条件。引入熵H表达式[7]

(11)

(12)

(13)

3 Kriging代理模型的构造

3.1 Kriging模型

Kriging模型由一个随机过程和一个线性回归模型构成，其模型可表示为[14]

f(xi)=fT(xi)β+z(xi)

(14)

式中β=[β1β2…βp]T是回归模型系数向量,f(xi)=[f1(xi)f2(xi)…fp(xi)]T是变量x的多项式函数,z(x)是协方差非零的随机分布，服从正态分布N(0,σ2)。利用最小二乘估计可得β和σ2的估计值，分别为

(15)

(16)

式中F为样本点向量矩阵，Y为样本点响应列向量，R为空间矩阵，其中元素Ri j=R(xi,xj) (i,j=1,2,…,d)，d为样本数。β和σ2均为θk的函数，则在Kriging模型中唯一决定预测响应精度的未知数是θk。

3.2 BA算法

BA算法从蝙蝠的行为和回声定位能力出发，利用脉冲发射率r和脉冲强度A的多样性来探测、捕食猎物和寻找栖息地，与其他算法相比，BA在准确性和有效性方面远优于其他算法，且没有许多参数要进行调整。因此，采用蝙蝠算法对其进行寻优以提高模型精度。式(17～19)给出了BA算法的主要方程[15]为

fi=fmin+τ(fmax-fmin)

(17)

(18)

(19)

(20)

(21)

3.3 Kriging模型的构造

在建立Kriging模型时，需要先设定相关系数θk，其直接影响Kriging模型的预测精度。采用拉丁超立方抽样方法抽取初始待修正参数区间内的样本，将其分为训练集和测试集(训练集和测试集数据互不相同)，并利用有限元方法计算所对应的响应值。然后，以训练集作为Kriging模型的输入，对应的结构应变模态频率和振型作为响应来构建Kriging模型；以测试集对应的试验响应值与Kriging模型预测响应值之间的均方根误差RMSE作为BA的目标函数即式(22)，寻求最优相关系数。最后，用寻得的最优相关系数构建Kriging模型。

(22)

4 改进的MCMC算法

4.1 MH算法

MH算法是一种基于模拟的MCMC技术，它的一个重要应用是从给定概率分布中抽样，使马尔可夫链稳态于该概率密度。即构造一个Markov链，通过抽样获取后验分布P(θ|D)的样本数据，根据大数定理可获得修正后参数θ的均值和方差。在模型修正中MH算法的步骤如下[16,17]。

(1) 选择具有物理意义的初始值θ0。

(2) 利用当前θt的值，依据建议分布q(·|θt)产生一个候选值θ′。

(3) 根据接受概率函数[16](23)计算接受率

(23)

(4) 从均匀分布u(0,1)中产生随机变量u。

(5) 当α(θt,θ′)>u时接受θ′，使θt +1=θ′，否则取θt +1=θt。

(6) 重复步骤(2～5)，直至收敛。

4.2 基于FPA的MCMC方法

目前，MCMC方法已成为处理复杂高维积分运算的有效工具，但随着参数维度增加，其后验分布更复杂，标准MH算法的样本接受率会随之降低，致使采样出现停滞现象。

FPA具有异花授粉(全局寻优)和自花授粉(局部寻优)两大特性。异花授粉过程可以逃离任何局部景观，进而探索更大的空间；自花授粉过程使相似的解更频繁地选择，从而保证更快地收敛，其收敛速度本质上是指数的，因此寻优效率更高[18]。本文采用FPA改进的MCMC抽样来提高采样接受率及样本多样性。即在N条满足条件的马尔可夫链中,采用FPA的异花授粉和自花授粉再一次寻优来更新候选参数。用待修正参数θ的目标函数J值来度量个体所处位置的优劣，将个体的优胜劣汰过程类比为搜索和优化过程中用较好的可行解取代较差可行解的迭代过程，因此，MCMC 算法的接受率有了很大程度提高的同时增加了样本多样性，提高了寻优效率，且保证样本收敛在合理的范围。

(24)

(25)

(26)

5 模型修正过程

根据上述理论，结合先验信息和已构造的Kriging代理模型产生的预测响应值，由贝叶斯统计理论推导后验概率密度函数。首先，用MCMC算法迭代求解N条马尔可夫链和每次迭代所得的最优值；然后，用FPA-MCMC算法迭代求解得到l条马尔可夫链，即可得到待修正参数的后验概率密度；最后，由后验分布统计特征获得修正后的参数值，完成模型修正，修正流程如图1所示。

6 数值算例

6.1 三自由度弹簧-质量系统

使用如图2所示的三自由度弹簧-质量系统对所提算法的有效性进行验证。该系统的确定性参数为m1=m2=m3=1.0 kg，k3=k4=10 N/m，k6=30 N/m。假设不确定参数的试验真实均值为k1=10 N/m，k2=14 N/m和k5=12 N/m，且服从独立正态分布。结构的固有频率为ω1=3.351 Hz，ω2=6.759 Hz，ω3=9.005 Hz，假定不确定参数的初始值是k1=6 N/m，k2=16 N/m，k5=10 N/m。选用k1，k2和k5为待修正参数，表示为θ={θ1,θ2,θ3}。使用贝叶斯公式计算其后验概率分布，用所提算法对该分布进行采样，迭代次数n=5000，待修正参数θi的取值范围为[5,20],初始值θ0={6,16,10}。

图1 模型修正流程

图2 三自由度弹簧-质量系统

分别采用标准MCMC算法与FPA-MCMC算法进行5000次迭代后得到待修正参数的马尔可夫链如图3所示。可以看出，两种算法得到的马尔科夫链均能收敛到试验真实值，但改进的FPA-MCMC 算法能够有效克服采样停滞的缺点，其抽样多样性优于标准MCMC算法。由式(23)计算得到改进算法的样本接受率为74.5%，与同维度下标准MCMC算法3%的接受率相比，改进算法的采样效果更好。其中，两种算法迭代所得θ1的马尔可夫链如图4所示，可以看出，所提算法遍历性远优于标准MCMC算法；图5为θ1的概率分布直方图，可以看出，其后验均值基本上与预设的真实值相吻合。

虽然初始值的选取不会影响收敛的最终结果，但对参数的选取误差越大，则燃烧期越长，因此去掉样本集的前300次迭代值，降低所选初始值不佳的影响。然后，检验后验样本的正态概率结果，如 图6 所示，可以看出，改进算法抽样得到的样本点基本与假设的参数分布线重合，因此样本的正态分布假设成立。并对不确定性参数θi的置信区间做出了概率性的估计，如图7所示，在置信椭圆内为含95%样本参数的区域。改进算法经过5000次迭代后，待修正参数基本位于置信椭圆轮廓内。

图3 三条马尔科夫链

图4 θ1两种算法的马尔科夫链

由后验分布统计特征得到修正后的不确定性参数值和响应值，修正结果列入表1和表2，可知待修正参数在修正后均接近假设的真实值，且修正后参数对应响应值的相对误差都小于0.02%，表明所提算法修正效果较好。

图5 θ1概率分布直方图

图6 参数的正态概率检验图

图7 置信椭圆的散点图

表1 修正参数的对比

表2 修正前后频率对比

6.2 三维桁架结构

图8 三维桁架模型

采用拉丁超立方抽样抽取2500组样本参数，分为2000组训练集和500组测试集，选取前6阶应变模态频率值和第60杆单元的前3阶应变模态振型值作为目标响应，通过训练集的样本和其响应值来构造Kriging模型。设置BA的种群大小为20，迭代次数为100，从区间[0.1,50]内寻得最优相关系数θk，其中第5阶频率的迭代过程如图9所示。

图9 第5阶频率迭代曲线

由优化所得θk值构建Kriging模型，并重新在参数取值范围内抽取500组样本，根据式(22)计算得到RMSE为1.808×10-7，其中第5阶频率的拟合情况如图10所示。可以看出，Kriging模型预测值和真实值几乎全部重合，说明所建立的Kriging模型预测精度很高，可以替代有限元模型进行模型修正。

图10 模型预测

按假设试验模型随机产生50组试验样本，将样本代入结构得到应变模态作为仿真试验响应，然后根据试验响应的均值和方差修正有限元模型，两种算法经20000次采样产生的马尔可夫链如图11所示。可以看出，改进算法产生的马尔可夫链一直在试验模型参数值附近变化，其接受率与遍历性均有所提高，改进算法的采样效果较好。为进一步验证所提方法的有效性，假定频率和模态振型的噪声在不同模态阶数上相互独立，对试验频率和模态振型加入5%的随机噪声来考虑试验数据的不确定性，得到如图12所示的马尔科夫链，可以看出，待修正参数均在均值附近收敛，对噪声具有良好的鲁棒性。

图11 两种算法的马尔科夫链

图12 FPA-MCMC算法马氏链

修正后参数的后验正态概率结果如图13所示。可以看出，改进算法抽样所得样本点与试验模型参数值分布线基本重合，能够较好地模拟后验样本。并对修正参数θi的置信区间做出了概率性的估计如图14所示，在置信椭圆内为含95%样本参数的区域，可以看出，待修正参数基本位于置信椭圆轮廓内的高概率区域。修正后参数的正态概率分布直方图如图15所示，可以看出，修正后参数的后验均值基本上与预设的真实值相吻合，表明后验分布的样本符合正态分布。

表3为参数θi的修正结果，可知当初始值选取误差较大时,仍可得到较精确的结果。表4 和表5为修正所得应变模态振型和频率的结果，可知其修正值和对应的响应值的相对误差均较小，因此所提算法的修正结果精度较高，模型修正效果较好。

图13 正态概率检验图

图14 置信椭圆的散点图

图15 参数的概率分布直方图

表3 修正参数对比

表4 修正前后模态振型对比

表5 修正前后频率对比

7 结 论

(1) 在贝叶斯统计理论的基础上，将FPA的局部优化和全局优化特性融入MCMC算法，保证了候选样本的合理性，也提高了寻优效率，且对噪声也具有一定的鲁棒性，实现了较高维待修正参数下的随机有限元模型修正。

(2) 在构造Kriging模型时，利用BA对Kriging模型的相关系数进行寻优，使所构造的Kriging模型具有良好的拟合精度和预测能力，能代替有限元模型进行迭代计算，使修正结果更可靠且提高了模型修正计算效率。

(3) 修正结果显示，所提修正方法增加了马尔可夫链的遍历性和样本合理性，同时大大地提高了样本接受率，克服了传统MH算法在参数维度高的情况下接受率低的缺点。