在土木工程领域中,有限元分析已成为结构计算不可或缺的手段。但受施工误差、离散误差以及物理参数、边界条件取值的影响[1],有限元模型的计算值与结构的实测值之间必然存在偏差,限制了有限元方法在结构损伤诊断、优化设计方面的应用。为此,有限元模型修正已成为土木工程领域一个重要的研究分支。按是否考虑待修正参数、测试数据等的不确定性,有限元模型修正可分为确定性方法和不确定性方法两类。国内外大量学者针对确定性模型修正方法进行了深入研究,取得了丰硕的成果。而不确定模型修正能综合考虑先验信息、测试样本信息等多种不确定因素,已成为有限元模型修正的前沿研究方向。
目前,不确定性有限元模型修正采用的理论方法有概率理论、模糊理论和区间理论等[2-4],其中最常用的是基于概率理论的贝叶斯方法,最初由Beck教授于1998年提出[5-6],为解决贝叶斯复杂后验分布的计算问题,其又于 2002年引入 MCMC(Markov chain Monte Carlo)算法,并以钢框架结构为例实现了2参数的不确定性有限元模型修正[7]。随后,基于 MCMC的贝叶斯有限元模型修正得到了快速发展,易伟建基于 MCMC算法成功实现了某 4层框架结构的有限元模型修正[8]。房长宇等[9]利用MCMC算法达成了某7根后张无粘结预应力简支梁的模型修正。韩方结合 MCMC算法和最大信息熵原理,对某电站地下厂房结构进行了有限元模型修正[10]。
针对传统 MCMC算法采样效率低下的难题,Ching等[11]提出了 TMCMC(transitional MCMC)算法,其基本思路是启用易于采样的目标概率密度函数代替复杂的目标概率密度函数,从而加快了采样效率,并在2自由度的结构上取得了较好的修正效果。Cheung等[12]建立了 HMCMC(Hamiltonian MCMC)算法,通过使用动力学方法,避免了 MH(metropolis-hastings)算法的局部随机游走行为,完成了 10层建筑的模型修正和不确定因素分析。刘纲等[13]建立了 DRAM-MH (delayed rejection adaptive metropolis)算法,通过引入延缓拒绝算法、自适应策略,解决了MH算法在高维待修正参数下收敛速度较慢以及采样“停滞”的问题,完成了 4层两跨钢框架实验室结构的有限元模型修正。
除以上几种方法外,迄今学术界已提出了较多的改进MCMC算法[14-16],这些算法提高了MCMC的计算效率,但这些算法多采用单条马尔科夫链,在待修正参数维数较高时计算精度不佳,且在混合目标分布下容易陷入局部最优解。针对这两个问题,本文引入差分进化马尔科夫算法,通过多条马氏链之间的差分运算来缩小条件分布探索的区间,自适应引导优化搜索的方向[17];在此基础上,引入子空间采样,通过自主选择维度空间采样以加快计算效率;引入异常链剔除算法以克服误采样,从而建立基于多链差分算法的贝叶斯有限元模型修正方法。最后,利用数值模拟与实验结构模型对所提出方法的有效性进行验证。
贝叶斯模型修正将待修正参数θ以及实测数据x均视为随机变量,设θ的初始(先验)分布为π(θ),则贝叶斯公式可表述为:
式中:p(θ|x)为修正后的参数概率分布,可简称为后验分布;p(x|θ)为θ给定下的条件分布,通常称为似然函数为x的边缘概率密度函数,是一个不依赖于θ的常数。
将实测获得的Nn阶结构响应量组成向量X,并假定X与有限元模型对应的理论计算值组成的向量X( θ)符合线性模型:
其中:ε ∈RNn,为测试误差、环境波动等所有不确定性误差产生的总和向量,通常假定其是均值为零、协方差为cov的正态分布,故似然函数p( x| θ)p(x|θ)采用向量方式可表示为:
其中,Ns为结构响应量的测试次数。为减少人为假定对修正结果的影响,取先验分布 π(θ)为均匀分布[18]:
其中,ai、bi为第i个待修正参数θi的上、下限值。假设结构响应量包括较易测试得到的结构频率ω和振型φ,将式(3)、式(4)代入式(1)可得:
式(5)表明,后验分布可直接通过积分进行求解。但由于实际结构中待修正参数维度较高,后验分布往往无法采用显示表达式,故需借助数学方法来获取近似的后验分布。由于式(5)中c为常数,因此仅选取J( θ)建立目标函数F( θ):
从而将求后验分布p(θ|X)的最大值变成了求J(θ)的最小值。
MCMC算法本质上是一个Monto Carlo综合程序,目的是在给定的初始值下,根据条件分布不断采样,假定马氏链在n时刻已经收敛了,则[θn,…,θT]的样本即为表征p(θ|X)的样本,也称为马氏链的平稳分布 π(x),再采用大数定理估计 π(x)的某些特征值,作为p(θ|X)的参数值,也即待修正参数更新后的值。
马氏链在开始采样阶段存在不稳定期,即样本统计的平均值或方差变化较大,为保证马氏链经过一定时间的采样后能平稳收敛到后验分布,需构造一个满足细致平稳条件的转移矩阵P:
式中,i,j∈[0,T],T为马氏链的迭代终止时刻。而实际中往往难以找到合适的转移矩阵P,故通常在某个已知的转移矩阵Q的基础上,引入接受概率α,将式(8)改写成如下形式:
式中,q( i, j)表示马氏链从状态i转移到状态j的概率,也称为条件分布q(j|i),为矩阵Q的元素。且α(i,j)满足以下条件:
此时P可通过Q乘以α得到,只需利用α(i,j)对马氏链从状态i以q(i,j)的概率跳转到状态j时进行评价,即是否接受转移,就能从概率的角度保证马氏链最终均能收敛到平稳分布。
具体采样步骤如下:
1)计算马氏链的初始状态θ0,即从先验分布中采样获取θ0。
2)在第t个时刻,利用当前θt值,从条件概率分布q(θ|θt)中采样获取候选样本θ*,并根据式(10)计算接受概率α(θt,θ*)。
3)从[0,1]的均匀分布中随机采样获取u,当u<α(θt,θ*)时,θt+1=θ*,否则,θt+1=θt。
4)重复步骤 2)~步骤 3),直到达到迭代收敛条件。
5)从获取的全部样本θ0,θ1,…,θT中去除非平稳样本(即n时刻之前的样本),再估计剩余样本的某些统计特征,作为修正后的θ值。
传统的 MCMC算法只有一条马氏链,且依赖于条件分布的选取,条件分布的方差较大时易出现经常拒绝候选点θ*的“停滞现象”;条件分布方差选取较小时易出现马氏链的“原地踏步现象”,从而使样本不再更新直至迭代结束,导致在高维待修正参数下计算精度与效率不佳。本节通过引入多链差分马氏链算法(DE-MC, differential evolutionmonto carlo)以克服高维参数模型修正过程中的“停滞现象”和“原地踏步现象”。在此基础上,引入子空间采样以及异常链检测的组合算法,通过自适应选择维度空间进行采样,并结合异常链检测算法,以提升高维待修正参数下的计算精度与效率,从而快速逼近目标分布。
DE-MC算法采用多条马氏链,例如N条马氏链Y 1 , Y 2 ,… ,YN 构成种群Y。然后通过各条马氏链之间的交叉、变异和选择操作实现条件分布的更新。在第t次进化(迭代)中,第i条马氏链的条件分布向量为:
式中,上标f和g表示随机选取的马氏链索引,f, g ∈ { 1,2,… ,i- 1 ,i+1,… ,N };γ是从[0.4,1]中随机抽取的常数;e是从正态分布N(0,k)d中抽取的向量组,其中d为待修正参数的维数,k通常取10 -6。
通过计算接受概率判别是否接受
,其中F为目标函数。当r>1时,
马氏链迭代初期所得样本是非平稳的,用来估计待修正参数的统计值会降低修正精度,因此采用Gelman等[19]提出的定量收敛准则(即比例缩小因子S)来区分非平稳和平稳阶段:
式中:C为每条马氏链的进化代数;D为用于评价的马氏链条数,且D≤N;B/C为D条马氏链平均值的方差;W为D条马氏链方差的平均值。当S<1.2时,表明此时样本开始趋于平稳。
通常采用比例缩小因子S来判定样本趋于平稳时停止迭代;也可根据大量试算,根据经验人为设定迭代次数。针对有限元模型修正,本文通过试算发现采用比例缩小因子停止迭代计算所得的平稳样本较少,而当迭代次数T≥300d时,所得平稳样本的数量才能精确获取修正后的参数值,所以推荐采用人为设定迭代次数T≥300d。
当待修正参数维数大幅增加时,如果仍对所有维度进行差分计算(采样),采样效率较低。此时,引入子空间采样以提高计算效率,即第i条链在第t次进化中采用部分维度来确定候选样本
式中,为第i条链的第t代的样本,下标
′为选取的采样维度(?≤d),初始迭代时选取?′=d。Id ′为′阶单位矩阵。ψd′、ξd′是从多元均匀分布U n (- q, q)、 N n (0,q′)中独立采样获取的,通常取
其中δ为用于产生候选样本的配对链对数;
为随机选取的平行链的向量差,
的计算式为:
式中,其 中
。
采样完成后,按式(15)判别是否接受
式中,CR为交叉概率,且CR=m/ nCR,其中m=1,2,… ,nCR,nCR =3;Uj是从[0,1]的均匀分布中抽取的随机数,j ∈[1 ,2,… ,d ′]。
多链差分算法在迭代初期存在马氏链的不稳定期,此时某些链可能会偏离正常轨迹较远(异常链),继续沿用此类异常链进行采样无疑会增加收敛至平稳期的时间。
因此,在非平稳期(即S≥1.2时)采用异常链检测准则(Inter-Quartile-Range(IQR))[20]进行异常链的统计与剔除:
式中,Q1、Q3表示各条马氏链在当前状态的下四分位数、上四分位数。当 IQ R <(Q1-Ω)/2时,记为异常链,其中Ω为各条马氏链在当前状态下至少50%样本后验密度的对数值。
将多链差分、子空间采样与异常链剔除算法融合,基于多链差分的贝叶斯有限元模型修正算法的具体步骤如下:
1)根据待修正参数的先验信息,随机产生N条链的初始样本,并计算其后验概率,用于异常链统计。设定迭代次数T,按式(11)开始进行N条链的平行进化。
2)在第t次进化中,针对第i(i=1,2,…,N)条链,按式(11)计算条件分布向量再按式(13)、式(14)产生候选样本
并按式(15)确定是否接受候选样本分量
以及进行维度更新。
3)按式(17)计算接受概率,并从[0,1]中随机产生u,如果
则接受新样本
否则仍保留
4)根据式(16)统计并剔除异常链。
5)根据式(12)判断马氏链是否平稳。
6)重复步骤2)~步骤5),直至迭代次数达到设定值T。
7)待修正参数θ特征值计算:采用种群中平稳期的所有样本估计θ的某些特征值,获取修正后的参数值。
建立如图1所示的简支梁模型,梁的横截面尺寸为0.4 m×0.6 m,初始弹性模量E0取3.3×1010 Pa,密度ρ取2500 kg/m3。沿梁长方向划分10个单元,记为编号①~⑩。同时设置 9个振型测点,记为编号 1~9。
图1 简支梁模型
Fig.1 Mode of simple supported beam
选取 10个单元的弹性模量作为待修正参数,并假定各单元弹性模量的真实均值Ei相互独立,且服从正态分布Ni(Ei,0.12),如表1所示。通过随机抽取Ei代入简支梁模型中,并采用求特征值和特征向量法求出结构的前6阶频率和前3阶振型,共模拟计算 20次。限于篇幅,文中仅给出第一阶频率与第三阶振型均值,分别如图2和图3所示。
表1 弹性模量的预设均值
Table 1 Default value of elastic modulus
?
图2 简支梁一阶频率
Fig.2 First order frequency of simple supported beam
图3 简支梁20次模态测试中第三阶振型均值
Fig.3 The average value of the third mode shape in 20 mode tests of simple supported beam.
根据 20次频率、振型实测值,采用多链差分进化算法进行模型修正,取N=10,T=3000,π(θ)为[0,1]的均匀分布,各参数的取值区间为[1.0E0,3.0E0]。然后根据第2.4小节的步骤1)~步骤7)得出修正参数的值。
限于篇幅,文中仅给出一条马氏链和三条马氏链在10个维度上的收敛情况,分别如图4、图5所示,图中右边框的“x”标记表示各参数的真实值。从图4和图5可以看出,所有参数最终均能很好的收敛到真实值。
图4 某条马尔科夫链的收敛情况
Fig.4 The convergence of a Markov chain
图5 三条马尔科夫链的收敛情况
Fig.5 The convergence of three Markov chains
为进一步验证算法精度,取各参数的初始值均为1.0 E0进行修正,修正前后的参数值和误差如图6所示。从图6可知,修正前各参数的均值误差基本在9%以上,单元9更是高达到47%左右,而修正后各参数均值误差基本降低至1.2%以下,单元9仅0.5%左右,因此,无论从整体或单个参数角度来看,所提算法的修正效果均较好。
图6 修正前后待修正参数值对比
Fig.6 Value of parameter comparison before and after updating
假定频率、振型的噪声在不同模态阶数和测试点上相互独立,将“实测”频率和振型加上一定的噪声来考虑实测数据的不确定性,噪声幅值分别取为2%、5%和10%。修正后各参数的后验误差如图7所示。从图7可以看出,随着噪声的增大,虽然各参数的误差值也逐渐增大,但基本保持在10%以内,表现出对噪声具有良好的容忍性。
图7 不同噪声下参数误差
Fig.7 Error of different parameters under different noise
文献[13]所提DRAM (delayed rejection adaptive metropolis)算法是对传统MCMC算法的改进,但其仍为单马氏链算法,为对比本文所提方法与单马氏链的计算效果,选取该算法进行对比,计算结果如表2、表3所示。由表2、表3可知,多链差分进化算法修正后的频率误差均低于1%,而DRAM算法各阶频率误差均大于 1%;本文方法所得振型MAC值均为0.999,而DRAM算法修正后的MAC值最高为0.981,表明本文所提算法在精度上明显优于 DRAM 算法。在计算时间方面,采用相同的个人计算机,本文方法所需时间为316 s,DRAM所花费的时间为110 s。
表2 两种算法频率修正的对比
Table.2 Frequency comparison after updating of the two kinds of algorithm
?
表3 两种算法振型修正后的对比
Table.3 Mode shape comparison after updating of the two kinds of algorithm
?
在实验室搭建四层两跨钢框架结构,结构的开间与层高均为 350 mm,梁、柱均由钢板组成,单个钢板尺寸为 350 mm×65 mm×4 mm,梁和柱通过L型节点板和两个螺栓连接,结构底部与刚性底座通过螺栓固定,基础采用12 mm厚钢板模拟,在楼层节点处共布置4个加速度传感器,如图8中实线圆圈所示。
图8 试验框架模型
Fig.8 Experiment frame model
在框架顶层施加水平向单点脉冲激励(力锤提供),测试框架结构加速度时程响应曲线,并利用频域分解法(FDD)获取结构的基本振动特性。为获取多次测试数据,共进行了 20次测试。限于篇幅,文中仅给出结构的前4阶频率以及第1阶、第4阶振型,如表4和图9所示。
表4 前4阶频率测试值
Table.4 The value of first four measured frequency
阶数 均值/Hz 标准差 标准差/均值/(%)第1阶 4.88 0.06 1.13第2阶 15.63 0.17 1.09第3阶 27.34 0.27 1.00第4阶 35.16 0.40 1.14
图9 第1阶和第4阶实测振型(均值)
Fig.9 The first and fourth measured mode shape
在有限元模型中,将各层梁、各跨柱分别划分为1个单元,并采用欧拉梁单元模拟,如图10所示。初始弹性模量E0均取210 GPa,初始质量密度ρ0均取为7800 kg/m3,泊松比取为0.3。将同一层2根梁的弹性模量E划分为一个待修正参数,质量密度ρ划分为一个待修正参数,并按楼层1~楼层4分别记为E1~E4,ρ1~ρ4;同一层 3根柱的E 划分为一个待修正参数,质量密度ρ划分为一个待修正参数,并按楼层 1~楼层 4 分别记为E5~E8,ρ5~ρ8。因此,整个框架共划分为16个待修正参数。
图10 单元号和传感器布设示意图
Fig.10 Unit number and sensor layout diagram
采用灵敏度分析对初始待修正参数进行筛选[21],即将E1~E8和ρ1~ρ8分别降低一定的比例,并计算所得前4阶频率的改变率,结果如图11所示。从图11可知,单元E5~E8和ρ1的灵敏度均小于2%(图11中虚线以下)。E5~E8灵敏度小是因为梁的刚度对框架侧向频率影响较小,且弹性模量与刚度呈正相关,故弹性模量对面内频率的影响也较小;ρ1灵敏度小是由于首层柱与地面相连,质量密度对频率的影响较小。在剔除以上5个灵敏度小于2%的参数后,选取剩余的11个参数作为最终的待修正参数。
图11 各单元物理参数的灵敏度
Fig.11 Sensitivity of the physical parameters of each unit
在多链差分法中,取N=10,并设定迭代次数T=5000,各参数初始区间分布为[0.5E0,1.5E0]。在DRAM 算法中,各参数的初始值均取为E0,迭代次数以及各参数初始区间分布的选取与多链差分进化算法一致。修正后的频率与振型MAC值分别如表5、表6所示。
表5 修正后各阶频率对比
Table 5 Frequency comparison after updating
阶数 测试均值/Hz DRAM算法 多链差分进化算法均值 误差/(%)均值 误差/(%)第1阶 4.88 4.87 0.24 4.90 0.34第2阶 15.63 15.60 0.16 15.73 0.65第3阶 27.34 26.72 2.28 26.82 1.91第4阶 35.16 35.87 2.04 35.34 0.52
表6 修正前后各阶振型的MAC值对比
Table 6 MAC values comparison before and after updating
各阶MAC值 第1阶 第2阶 第3阶 第4阶修正前 0.972 0.986 0.964 0.986 DRAM 0.998 0.996 0.990 0.996多链差分进化算法 0.999 0.999 0.998 0.999
从表5可知,两种算法修正后的前3阶频率的误差基本相当,但对于第4阶频率,多链差分算法的精度明显高于DRAM法;从表6可知,对于结构的各阶振型,多链差分算法的 MAC值均大于DRAM法的MAC值。总的来看,多链差分进化算法较 DRAM 算法在高阶频率和各阶振型的修正精度有较大优势。
从计算时间来看,本文方法耗时 5 h;DRAM算法耗时 3 h。这主要是因为本文方法采用多条马氏链并行运算,其样本数量远大于 DRAM 算法。应指出的是,有限元模型修正往往不需要在线、实时完成,且随着当前云计算等的发展,有望弥补本文方法耗时量大的不足。
针对MCMC算法在高维待修正参数下的计算精度不佳、过于依赖条件分布的选取、在混合目标分布下容易陷入局部最优等问题,本文提出了结合差分进化、子空间采样以及异常链检测的多链差分进化有限元模型修正方法,通过数值算例与试验结构的模型修正表明:
(1)在无测试噪声影响下,简支梁数值模型采用多链差分算法修正后,待修正参数的误差由 9%以上降低至1.2%以下,在加入噪声后,待修正参数的误差仍小于 10%,表明本文方法的修正精度较高,且具有良好的抗噪性。
(2)不同算法修正结果表明,DRAM算法修正后各阶频率误差均在 1%以上,而多链差分进化算法修正后频率误差降至1%以下,同时振型的MAC值接近 1.0,表明多链差分方法修正精度高于DRAM算法。
(3)框架试验模型修正中,多链差分进化算法相比 DRAM 算法在高阶频率与各阶振型上的修正精度更高,这对实际结构中在多参数下的有限元模型修正具有重要的意义。
[1]Mosavi A A, Sedarat H, O’Connor S M, et al.Calibrating a high-fidelity finite element model of a highway bridge using a multi-variable sensitivity-based optimization approach[J].Structure and Infrastructure Engineering,2014, 10(5): 627―642.
[2]Collins J D, Hart G C, Hasselman T K, et al.Statistical identification of structures[J].AIAA Journal, 1974,12(2): 185―190.
[3]Jiang C, Han X, Liu G R.Optimization of structures with uncertain constraints based on convex model and satisfaction degree of interval[J].Computer Methods in Applied Mechanics and Engineering, 2007, 196(49):4791―4800.
[4]Hanss M.Applied fuzzy arithmetic[M].Berlin Heidelberg: Springer-Verlag, 2005: 128.
[5]Beck J L, Katafygiotis L S.Updating models and their uncertainties I: Bayesian statistical framework[J].Journal of Engineering Mechanics, 1998, 124(4): 455―461.
[6]Katafygiotis L S, Beck J L.Updating models and their uncertainties II: Model identify ability[J].Journal of Engineering Mechanics, 1998, 124(4): 463―467.
[7]Beck J L, Au S K.Bayesian updating of structural models and reliability using Markov chain Monte Carlo simulation[J].Journal of Engineering Mechanics-ASCE,2002, 128(4): 380―391.
[8]易伟建, 周云, 李浩,等.基于贝叶斯统计推断的框结构损伤诊断研究[J].工程力学, 2009, 26(5): 121―129.Yi Weijian, Zhou Yun, Li Hao, et al.Damage assessment research on frame structure based on Bayesian statistical inference[J].Engineering Mechanics, 2009, 26(5):121―129.(in Chinese)
[9]房长宇, 张耀庭.基于参数不确定性的预应力混凝土梁模型修正[J].华中科技大学学报(自然科学版),2011, 39(11): 87―91.Fang Changyu, Zhang Yaoting.Model updating of prestressed concrete beams using parameters uncertainty[J].Journal of Huazhong University of Science and Technology (Natural Science Edition), 2011, 39(11):87―91.(in Chinese)
[10]韩芳, 钟冬望, 汪君.基于贝叶斯法的复杂有限元模型修正研究[J].振动与冲击, 2012, 31(1): 39―43.Han Fang, Zhong Dongwang, Wang Jun.Complicated finite element model updating based on Bayesian method[J].Journal of Vibration and Shock, 2012, 31(1): 39―43.(in Chinese)
[11]Ching J Y, Chen Y C.Transitional Markov chain Monte Carlo method for Bayesian model updating, model class selection, and model averaging[J].Journal of Engineering Mechanics, 2007, 133(7): 816―832.
[12]Cheung S H, Beck J L.Bayesian model updating using hybrid Monte Carlo simulation with application to structural dynamic models with many uncertain parameters[J].Journal of Engineering Mechanics, 2009,135(4): 243―255.
[13]刘纲, 罗钧, 秦阳, 等.基于改进 MCMC 方法的有限元模型修正研究[J].工程力学, 2016, 33(6): 138―145.Liu Gang, Luo Jun, Qin Yang, et al.A finite element model updating method based on improved MCMC method[J].Engineering Mechanics, 2016, 33(6): 138―145.(in Chinese)
[14]Wang F S, Li X C, Lu M Y.Adaptive Hamiltonian MCMC sampling for robust visual tracking[J].Multimedia Tools and Applications, 2017, 76(11):13087―13106.
[15]Agostinho N B, Machado K S, Werhli A V.Inference of regulatory networks with a convergence improved MCMC sampler[J].BMC Bioinformatics, 2015,306(16): 1―10.
[16]Kennedy D A, Dukic V, Dwyer G.Combining principal component analysis with parameter line-searches to improve the efficacy of Metropolis-Hastings MCMC[J].Environmental and Ecological Statistics, 2015, 22(2):247―274.
[17]Mengersen K L, Robert C P.IID sampling using self-avoiding population Monte Carlo: The pinball sampler[J].Bayesian Statistics, 2003, 7: 277―292.
[18]茆诗松, 程依明, 濮晓龙, 等.概率论与数理统计教程[M].北京: 高等教育出版社, 2011: 107―108.Mao Shisong, Chen Yiming, Pu Xiaolong, et al.Probability and statistics[M].Beijing: Higher Education Press, 2011: 107―108.(in Chinese)
[19]Gelman A, Rubin D B.Inference from iterative simulation using multiple sequences[J].Statistical Science,1992, 7(4): 457―472.
[20]Vrugt J A, ter Braak C J F, Diks C G H, et al.Accelerating markov chain monte carlo simulation by differential evolution with self-adaptive randomized subspace sampling[J].International Journal of Nonlinear Sciences and Numerical Simulation, 2009, 10(3): 273―290.
[21]Kim G H, Park Y S.An improved updating parameter selection method and finite element model update using multiobjective optimisation technique[J].Mechanical Systems and Signal Processing, 2014, 18(1): 59―78.
BAYESIAN FINITE ELEMENT MODEL UPDATING METHOD BASED ON MULTI-CHAIN DIFFERENTIAL EVOLUTION