结构损伤识别方法可分为基于数据与基于模型两大类，其中基于数据的方法有神经网络与支持向量机等，基于模型的方法有贝叶斯方法、优化算法等．基于模型的损伤识别方法须建立一个能反映实际结构的基准有限元模型，因此无法避免有限元模型建立过程中带来的误差．这些误差属于损伤识别的不确定性(包括测量噪声、方法误差和建模误差)，采用优化算法不能兼顾到这些不确定性．采用贝叶斯损伤识别方法时可以考虑这些不确定性．文献[1]提出了一种用于多层框架结构损伤检测的贝叶斯概率方法．文献[2]提出了一种有限长周期支撑结构损伤概率识别方法，结合贝叶斯理论与马尔科夫蒙特卡罗(MCMC)模拟算法，实现基于测量自振频率变化的有限长周期支撑结构单元损伤的概率识别．文献[3]提出了一种新型混合贝叶斯优化算法，该方法可以有效地加快算法的收敛速度和收敛精度．文献[4]实现了有限元分析程序和复合链抽样技术的协同，得到了模型参数后验分布的统计特征．文献[5]提出了一种自适应采样与延迟拒绝算法相结合的参数估计新方法，有效提高了抽样效率．对于贝叶斯损伤识别中的后验概率密度函数的求解，最常用的方法为MCMC方法．文献[6]将MCMC方法引入到贝叶斯损伤识别领域，避免直接从修正参数的后验概率密度函数抽样，通过较简单的中间条件概率逐步抽样，但此方法用于高维问题存在的抽样困难．当标准MH算法在高维参数空间采样时，建议分布选取不当会对采样结果造成影响[7]，导致损伤识别结果不精确．文献[8]将延缓拒绝自适应MH算法引入到模型修正中，通过自适应算法实现自主调整采样步距；通过延缓拒绝算法提高新样本接受概率，从而有效克服了标准MH算法对高维待修正参数收敛较慢或无法收敛的问题．文献[9]将变异算子和PSO算法应用到每个粒子上作为MCMC方法的转移核，粒子有可能移动到一个更好位置，直到收敛，但无法实现自适应．针对标准MH算法易受建议分布影响的问题，本研究提出了一种基于粒子群算法中粒子位置更新方式的自适应MH算法(自适应PSO-MH算法)，该方法无须选择建议分布，可以有效避免因建议分布选取不当而导致马尔科夫链收敛不精确的问题．1 贝叶斯损伤识别方法1.1　结构损伤参数结构损坏后，结构的质量改变通常可以忽略不计，其损伤程度可通过结构刚度折减因子θ=(θ1,θ2,⋯,θn)T来反映，对于完整结构，θ=(0,0,⋯,0)T．结构单元刚度矩阵与整体刚度矩阵之间的关系可以表示为[10]ked=(1-θe)keu    (0≤θe≤1)；(1)Ked=TrTkedTr；(2)Kd=∑i=1nK¯ed，(3)式中：ked和keu分别为损伤情况下与完整情况下的单元刚度矩阵；Tr为转换矩阵；Ked为坐标转换后的单元刚度矩阵；K¯ed为整体坐标系下的单元刚度矩阵；Kd为整体刚度矩阵；θe为第e个单元的刚度折减因子，e为单元编号，e=1，2，…，n，n为结构单元数．1.2　贝叶斯理论基本公式假设θ的先验分布为π(θ)，则样本X的联合概率密度分布可表示为p(X|θ)=∏l=1Nvp(xl|θ)，(4)式中：θ为后验分布的自变量；l为样本数，l=1，2，…，Nv，Nv为样本总数．后验分布可表示为π(θ|X)=p(X|θ)π(θ)∫θp(X|θ)π(θ)dθ=cp(X|θ)π(θ), (5)式中∫θp(X|θ)π(θ)dθ为一个常数，c为其倒数．1.3　基于贝叶斯方法的损伤识别目标函数基于模态参数推导的后验概率密度函数(目标函数)已经在文献[11-12]中有研究，本研究采用均匀先验分布[8]推导的后验概率密度函数式[13]，即π(θ|D)=cexp(-J(θ)/2); (6)      J(θ)=∑j=1Nm∑r=1m(f̂r,j-fr(θ))2cov fr+(ψ̂r,j-ψr(θ))T(ψ̂r,j-ψr(θ))cov ψr,式中：D为实测自振频率与振型；f̂r,j和ψ̂r,j为分别为第j次实测第r阶的自振频率和振型；fr(θ)和ψr(θ)分别为第r阶理论自振频率与振型，须要注意的是频率残差远大于振型残差；r为模态阶数，r=1，2，…，m，m为最大模态阶数；j为模态测试次数，j=1，2，…，Nm，Nm为总测试次数；cov f和cov ψ分别为自振频率和振型的协方差，cov fr和cov ψr为其第r个对角线元素，cov f和cov ψ的计算方法为cov f=diag(σ12,σ22,⋯,σm2)，σr2=1Nm∑j=1Nm(f̂r,j-f̂¯r)2；cov ψ=diag(λ12,λ22,⋯,λm2)，λr2=1Nm∑j=1Nm||ψ̂r,j-ψ̂¯r,j||2，式中：f̂¯r为第r阶实测自振频率的均值；ψ̂r,j为第j次实测第r阶振型的均值．结构损伤参数通常是一个高维向量，其后验分布复杂，须采用MCMC方法来近似求解后验分布．2 MCMC采样方法MCMC方法建立一个与目标问题相似的概率模型，通过大量抽样获取一个平稳分布的马尔科夫链，样本统计特征值可作为目标问题近似解．2.1　标准MH算法MH算法须要使用两个分布：一个是平稳分布(后验概率密度分布)；另一个是建议分布，用于产生新的抽样样本．MH抽样算法的具体步骤如下．步骤1 初始化马尔科夫链的状态θi=θ1，令i=1,i=1,2,…,N，N为总迭代次数；步骤2 从建议分布q中抽取一个候选样本值θ*；步骤3 计算候选样本值θ*的接受概率α(θi,θ*)=min1,π(θ*|D)q(θi|θ*)π(θi|D)q(θ*|θi)；步骤4 从均匀分布U(0，1)中生成一个随机数u，若u≤α(θi,θ*)，则接受候选样本值θ*，令θi+1=θ*，否则令θi+1=θi；步骤5 令i=i+1，重复步骤2~步骤4，直到得到一个平稳分布的马尔科夫链．2.2　自适应PSO-MH算法a. PSO算法中粒子位置更新方式为[14]Vi=c1r1(θiibest-θi)+c2r2β(i)(θigbest-θi)；(7)θi+1=θi+Vi，(8)式中：Vi为第i次迭代的速度；θigbest和θiibest分别为种群最优值和个体最优值；θi为第i次迭代后损伤参数；c1和c2为加速度因子，均取1.494 45；r1和r2为从区间[0，1]内随机抽取的随机数；β(i)为步长缩放因子，且β(i)=μexp[-0.5(i/N)2]，其中μ为缩放系数，取μ=4．θigbest的更新方式为θigbest=θibest    (i≤N0或D(θi,θi+1,⋯,θi-N0)≥vt);      E(θi,θi+1,⋯,θi-N0)  (iN0且D(θi,θi+1,⋯,θi-N0)vt),式中：θibest为自适应样本值；E(X)与D(X)分别为X的期望与方差；N0为初始适应阶段长度，本研究取400；vt为方差容差，本研究取1×10-5．θiibest的更新方式为：θ1ibest=θ1； D1=|θiibest-θigbset|； D2=|θi+1-θigbest|．若D2D1，则θi+1ibest=θi+1，否则θi+1ibest=θiibest．b. 自适应PSO-MH算法的采样步骤如下：步骤1 初始化马尔科夫链；步骤2 设置初始适应阶段长度N0与容差vt；步骤3 从自适应均匀分布U(θibest-Diadp,θibest+Diadp)中随机抽取θibest，其中Diadp=exp(θi-5)；步骤4 采用式(8)产生马尔科夫链候选值θ*并计算θ*与θi的后验概率密度函数；步骤5 计算样本接受率r=min1,π(θ*|D)π(θi|D)；步骤6 从均匀分布U(0，1)中生成一个随机数u，若u≤r，则接受候选样本值θ*，令θi+1=θ*，否则令θi+1=θi；步骤7 重复步骤3~6，直到得到一个平稳分布的马尔科夫链．3 数值算例采用如图1所示的混凝土横梁数值算例来验证自适应PSO-MH算法的统计效率．将横梁划分为16个单元，每个单元的材料属性为：弹性模量E=3.45×1010 Pa；密度ρ=2 500 kg/m3；梁长8 m；宽0.3 m；高0.1 m．10.13245/j.hust.220818.F001图1横梁示意图损伤参数的真实分布以实际损伤程度为均值，方差为0.01的正太分布来模拟．图2和3分别为自适应PSO-MH算法与标准MH算法生成的马尔科夫链．在MCMC过程中抽取样本意味着连续样本之间可能存在相关性．为了评估序列的相关性，本研究绘制了每个参数在不同滞后次数下的样本ACF函数图(见图4和5，图中：F为自相关函数值；L为滞后次数)，均取后1 000次迭代数据．表1为自适应PSO-MH算法与标准MH算法的马尔科夫链估计值与蒙特卡罗标准差(MCSE)．10.13245/j.hust.220818.F002图2自适应PSO-MH算法生成的马尔科夫链10.13245/j.hust.220818.F003图3标准MH算法生成的马尔科夫链10.13245/j.hust.220818.F004图4自适应PSO-MH算法ACF函数图10.13245/j.hust.220818.F005图5标准MH算法ACF函数图10.13245/j.hust.220818.T001表1马尔科夫链估计值与蒙特卡罗标准差损伤参数损伤单元实际损伤自适应PSO-MH算法标准MH算法估计值误差/%标准差/10-5估计值误差/%标准差/10-4θ130.050.0511.317.080.0523.234.31θ250.100.1000.493.580.0981.841.81θ370.300.2980.544.610.2971.118.42θ4100.500.5010.197.630.5010.235.36从图2和3可看出：自适应PSO-MH算法和标准MH算法相比，可以更快收敛到实际损伤值，且生成的马尔科夫链更平稳．从图4和5可看出：自适应PSO-MH算法生成的马尔科夫链ACF值呈现出截尾状态，表明马尔科夫链的连续样本间相关性低，同时进一步证明生成的马尔科夫链为平稳的．而标准MH算法的ACF值始终呈现出拖尾状态，表明马尔科夫链高度自相关．具有高度相关性的马尔可夫链在参数空间中缓慢移动，须要更多迭代和更长计算时间才能以接近目标分布的概率访问参数空间不同区域．从表1可看出：两种算法生成的马尔科夫链都能收敛到实际损伤值，但是自适应PSO-MH算法生成的马尔科夫链的估计值比标准MH算法的估计值误差更小，同时蒙特卡罗标准差也更小，表明自适应PSO-MH算法生成的马尔科夫链统计效果优于标准MH算法．4 试验研究采用LANL实验室的一个八自由度结构[15]作为试验实例来验证自适应PSO-MH算法的有效性．该结构由八个弹簧连接的质量块组成，通过改变连接两个质量块的弹簧刚度特性来模拟损伤．每个质量块由一个带孔的圆形铝盘与聚四氟乙烯衬套组成，厚25.4 mm，直径76.2 mm．质量块在一个高度抛光的钢杆上滑动，钢杆支撑质量块并约束它们沿钢杆移动．靠近激振器的质量块质量为559.3 g，其他质量块的质量为419.4 g．用刚度系数小于原弹簧的另一个线性弹簧替换原弹簧，模拟模型中的线性损伤．替换的弹簧可位于任何相邻质量之间，从而模拟不同的损伤位置．替换的弹簧可以有不同程度的刚度折减，以模拟不同程度的损伤．原弹簧的刚度系数为56.7 kN/m．通过将第五个弹簧替换为另一个刚度较低的弹簧(刚度降低55%)来模拟损伤[16-17]．该结构由一个电动激振器随机激发，测量每个质量块的水平加速度，得到八组加速度数据，对加速度数据进行快速傅里叶变换，得到结构的自振频率与振型．4.1　模型修正采用PSO算法对有限元模型进行修正，实测结构的频率和振型与修正后模型的频率和振型及模态置信准则(MAC)如表2和图6所示，图6中：A为相对振幅；O为测点号．一阶模态为结构的刚体模态，已忽略．10.13245/j.hust.220818.T002表2试验与数值模态参数模态阶数频率/Hz误差/%MAC实验数值222.8823.000.530.998 7344.7544.690.130.998 4465.3865.600.340.998 610.13245/j.hust.220818.F006图6无损状态下结构实测和数值振型从表2可看出：实测频率与数值模型频率非常接近，误差较小．同时，从图6中可发现二者的模态振型也非常相似，对比表2中的MAC值可以看出，各阶振型的MAC值基本接近于1，表明二者振型匹配效果良好．以上结果表明模型修正效果显著，修正后模型与实际结构的模态高度匹配，可用于后续损伤识别．4.2　损伤识别首先，基于修正的结构有限元模型构建关于损伤程度与模态参数的函数，即输入结构损伤单元位置及损伤程度，输出结构自振频率与振型；然后，利用输出的模态参数构建贝叶斯损伤识别的后验概率密度函数，采用自适应PSO-MH算法与标准MH算法求解该后验概率密度函数，生成实际损伤程度的马尔科夫链(见图7)；最后，根据生成的马尔科夫链估计值得到损伤程度，马尔科夫链的估计值与蒙特卡罗标准差如表3所示．为了评估马尔科夫链的自相关性，绘制了其ACF函数图(见图8)．10.13245/j.hust.220818.F007图7实际损伤程度的马尔科夫链10.13245/j.hust.220818.T003表3马尔科夫链估计值与蒙特卡罗标准差采样方法估计值误差/%标准差/10-4PSO-MH0.5283.946.04MH0.5166.2411.6010.13245/j.hust.220818.F008图8ACF函数图从图7和8可看出：自适应PSO-MH算法在实际结构中生成的马尔科夫链依然平稳，ACF值呈现出截尾状态，表明马尔科夫链的连续样本间相关性低，而对于标准MH算法，ACF值则呈现拖尾状态，表明马尔科夫链高度自相关．通过表3可看出：相比于标准MH算法，自适应PSO-MH算法的马尔科夫链估计值与实际损伤程度更接近，且蒙特卡罗标准差更小，统计效果更佳．5 结论本研究将PSO算法中粒子位置更新方式融入到标准MH算法中，并引入自适应均分分布抽样来更新粒子最优值，构建了自适应PSO-MH算法．通过一个横梁数值算例对比了该算法与标准MH算法之间的差异，并采用LANL实验室的一个八自由度结构试验验证了该算法的有效性．a. 自适应PSO-MH算法无须选取建议分布，避免了标准MH算法因建议分布选取不当而造成的损伤识别结果不精确．通过数值和试验结果可以看出，自适应PSO-MH算法损伤识别精度高于标准MH算法．b. 自适应PSO-MH算法生成的马尔科夫链平稳，蒙特卡罗标准差小，统计效果优于标准MH算法．c. 自适应PSO-MH算法生成的马尔科夫链ACF值呈截尾状态，马尔科夫链的连续样本间相关性低，所需计算时间更少，收敛速度更快．