基于改进SCAM算法的结构物理参数识别贝叶斯方法

赵一男，公茂盛，左占宣，高艳滨

(中国地震局工程力学研究所 中国地震局地震工程与工程振动重点实验室，哈尔滨 150080)

结构健康监测是目前土木工程领域重要研究方向，而结构损伤识别是结构健康监测技术的核心。结构损伤表现为结构物理参数的变化(如刚度降低)，或由此导致的模态参数的变化(如频率下降)，因此通过对损伤前后的结构物理参数或模态参数进行识别可以判断损伤位置和损伤程度。损伤识别方法一般分为确定性方法和不确定性方法[1-2]，实际工程中，由于结构本身属性复杂、环境噪声[3]等不确定性因素影响，使得采用基于频率、振型等确定性方法[4-5]识别结构损伤的准确性和可靠性降低。于是学者们发展确定性方法的同时提出了不确定性方法，而基于贝叶斯估计的结构损伤识别方法因为可以很好地解释损伤识别过程中所存在的不确定性，成为目前广受关注和应用的不确定性方法之一。

早在1989年，Beck[6]基于贝叶斯统计框架提出了一种系统识别方法，并明确指出了识别过程中存在不确定性。基于贝叶斯估计的结构损伤识别可以给出结构参数后验联合概率分布，通过求解参数后验边缘概率分布和最优估计值判断损伤位置和损伤程度，求解后验联合概率分布通常需要积分运算，但实际情况中，结构参数的后验联合概率分布是复杂的、高维的，很难通过直接积分求解边缘概率分布和最优估计值[7]。为解决这一问题，Hastings等[8-11]提出和选择马尔可夫蒙特卡罗(Markov Chain Monte Carlo，MCMC)抽样方法，通过在概率空间中随机抽样来逼近参数的后验分布[12]。但实践发现，MCMC算法仍存在抽样收敛速度慢，不能求解高维概率分布等问题，为此Haario等[13]提出了自适应Metropolis算法，提高了传统MCMC算法的收敛速度，并指出合理选择提议分布是算法收敛的关键，之后又对自适应Metropolis算法进行了改进，提出了逐分量自适应Metropolis(single component adaptive Metropolis，SCAM)算法，该算法可以有效求解高维问题[14]。Simoen等[15]利用SCAM算法对7层结构的损伤进行了识别，表明了SCAM算法的实用性，但SCAM算法的提议分布方差序列不满足马尔可夫性质，特别是当选取的抽样参数初始值与要识别的真实值相差较大且提议分布方差较小时，抽样效率低、易生成重复样本序列，从而导致结果误差较大[16]。

尽管如此，SCAM算法仍然是目前解决贝叶斯估计中求解高维边缘概率分布的常用方法，广泛地应用于基于贝叶斯估计的结构损伤识别。本文首先针对SCAM算法的不足，重新定义了提议分布方差计算式，给出了改进的SCAM算法，进一步提高了抽样效率和计算结果准确性。然后利用贝叶斯理论给出结构物理参数的后验联合概率分布，再结合改进的SCAM算法求解结构物理参数的后验边缘概率密度分布和最优估计值，从而实现对结构损伤进行识别和评估。最后通过对一个4层数值结构模型进行了识别，验证了本文所提出算法的有效性和可靠性。

1 基本理论与方法改进

1.1 结构物理参数识别贝叶斯方法

对于Nθ个自由度的工程结构，根据最大熵原理，结构刚度参数θ服从伽玛分布，记θ～Γ(α,β)，表达式为

(1)

式中：θMPE为振动测试前结构刚度参数极大先验估计值；α为形状因子；β为逆尺度因子。α,β与参数θ的数学期望μ和标准差σ有关，分别由式(2)和式(3)确定

(2)

(3)

式中，cvi为变异系数。

(4)

(5)

根据贝叶斯公式，可以求得基于实测响应的结构物理参数θ的后验联合概率密度函数

(6)

式中，C为标准化常量，保证后验概率密度函数在定义域内积分为1。通过式(6)即可得到结构参数θ的后验概率分布和最优估计值。

1.2 SCAM抽样算法及改进

为了求解式(6)后验联合概率分布的边缘分布，Haario等提出的逐分量自适应Metropolis算法是目前较为有效的一种抽样方法，得到了较为广泛地应用。但在SCAM算法中，当前状态下提议分布方差并不只由上一状态的样本候选值决定，而是由之前所有状态的样本候选值共同决定，这样生成的随机样本序列将不再满足马尔可夫链无后效性，不能完全保证样本序列收敛，会生成大量病态样本，导致结果误差较大，因此提议分布的方差直接影响算法收敛性。

研究表明，无论提议分布方差过大还是过小都会导致算法抽样效率变低，所得最终结果较不可靠。所以必须选择合适的提议分布方差使得马尔可夫链稳定。在提议分布为正态分布的多维问题中，Roberts等[17-20]给出了1～10维情况下的最优接受率取值，且证明当维度逐渐增加时，最优接受率逐渐减小并趋近于0.234。因此，算法改进的核心问题在于如何合适地调整提议分布方差，使得抽样过程的接受率收敛到最优接受率。

根据上述分析及最优接受率取值，本文对SCAM算法中的提议分布方差进行重新定义，改进了SCAM算法。新的提议分布方差由式(7)确定

νt=νt-1(1+α(Xt-1,z)-αOAR)2

(7)

式中：νt为t+1时刻提议分布方差；αOAR为最优接受率;α(Xt-1,z)为样本值为Xt-1，候选值为z时的接受率。由式(7)可知，若取αOAR≈0.234，则α(Xt-1,z)-αOAR∈[-0.234,0.766]，下一状态的提议分布方差取值范围为νt∈[0.7662νt-1,1.7662νt-1]。当α(Xt-1,z)≤αOAR时,νt会减小到0.7662νt-1～νt-1；当α(Xt-1,z)>αOAR时，νt会增大到νt-1～1.7662νt-1，随着迭代步数增加不断调节提议分布方差的数值，使得接受率趋近于最优。此外，为防止提议分布方差过大或过小，可将其限制在合理区间内。这样对提议分布方差改进后，SCAM算法在实现自适应性同时，也使得抽样样本序列满足马尔可夫链无后效性。改进的SCAM算法计算流程如图1所示。

图1 改进的SCAM算法流程图Fig.1 Flow chart of improved SCAM algorithm

1.3 算法验证

本文建立了一个抽样算例，对二维随机向量{x1,x2}独立同分布于α=6，β=1的伽玛分布进行随机抽样，峰值点的横坐标作为标准值，对改进的SCAM算法可靠性进行验证。

二维随机向量{x1,x2}联合概率密度函数为

(8)

初始向量值选择[40,40]，初始提议分布方差取[0.1,0.1]，迭代步数定为10 000步，“burn-in”阶段步数定为1 500步，分别用SCAM算法和本文改进的SCAM算法对式(8)进行随机抽样。

随机变量x1的马尔可夫链历史路径图如图2所示，可以看出原始SCAM算法在“burn-in”阶段的接受率较高，而在“burn-in”阶段之后的接受率较低，出现大量重复样本，马尔可夫链稳定性较差。原因在于SCAM算法“burn-in”阶段的提议分布方差为固定值，如果选择不合理，计算效率将大大降低，“burn-in”阶段之后的提议方差不满足马尔可夫链无后效性，计算结果误差较大。本文改进后SCAM算法的提议分布方差一直受接受率的影响而不断调整，使得马尔可夫链具有较好的稳定性，即使选择的初始值与标准值相差较大也可以很快达到收敛。图3为上述两种算法所得随机变量x1的累计均值图，可以看出，改进的SCAM算法收敛速度明显快于SCAM算法，SCAM算法在2 500步左右出现收敛趋势，而改进的SCAM算法则在500步左右就能达到收敛。图4给出了运用上述两种算法进行随机抽样所得随机变量x1样本边缘概率密度曲线与理论边缘概率密度曲线的对比情况，从图中可以看出SCAM算法得到的样本概率密度曲线与理论概率密度曲线差别较大，而改进的SCAM算法得到的样本概率密度曲线与理论概率密度曲线基本吻合。

图4 随机变量x1样本边缘概率密度曲线与理论边缘概率密度曲线对比Fig.4 Comparison of the sample and theoretical marginal probability density curves of random variable x1

图3 随机变量x1累计均值随迭代步数变化趋势Fig.3 Variation trend of the cumulative mean of random variable x1 with the number of iteration steps

图2 随机变量x1马尔可夫链路径图Fig.2 Path diagram of the Markov chain of random variable x1

综上所述，本文改进的SCAM算法克服了原始SCAM算法存在的不足，具有较高的抽样效率和可靠性。

2 结构识别算例

2.1 结构数值模型

考虑到土木工程结构的损伤大多为累积损伤，建立了一个4层剪切型结构数值模型，假定了多个损伤状态，如图5所示，对该模型在多次累积损伤状态下的物理参数进行识别。结构设计参数见表1，损伤状态S0～S4见表2，损伤程度依次为S0

表2 损伤状态表Tab.2 Damage cases

表1 结构设计参数Tab.1 Structural parameters

图5 四层结构模型Fig.5 4-story structure model

算例中，假定各损伤状态下的阻尼比ζ均为0.05，各层质量m不变，待识别的物理参数只有结构各层刚度θ(k)=[k1,k2,k3,k4]。各层刚度参数变异系数取0.2，以保证先验分布平滑且分布较宽，通过式(6)确定刚度参数的后验联合概率分布。对结构基底施加峰值加速度为0.05g的白噪声激励，采用Newmark-β法计算结构响应，将各层加速度响应作为观测值，并分别对无噪声、5%噪声和10%噪声三种噪声水平下的各工况层间刚度进行识别。初始刚度参数极大先验估计值取S0损伤状态刚度参数理论值，后续各状态刚度参数的极大先验估计值取上一状态的极大后验估计值，通过式(6)求得刚度参数后验分布，最后利用图1的计算流程，结合改进的SCAM算法，对各层刚度参数进行逐分量抽样，求得刚度参数后验边缘概率分布，并取刚度参数后验边缘概率密度曲线峰值点的横坐标作为最优估计值。

2.2 识别结果分析

本文采用改进的SCAM算法计算刚度参数的后验边缘概率分布，各损伤状态下各层刚度参数后验边缘概率分布如图6所示。本文只给出了5%噪声水平下的刚度参数边缘概率分布，无噪声和10%水平噪声下的识别结果与之类似，不再给出。图中可以清晰地看出结构各损伤状态下各层刚度参数的后验概率分布和最优估计值随损伤程度增加的变化情况。

图6 各损伤状态下各层刚度参数边缘概率分布(5%噪声)Fig.6 Marginal probability distribution of stiffness parameters of each story ineach damage cases (noise level 5%)

各层刚度参数理论值和最优估计值变化趋势如图7所示，从图中可以看出，各层刚度参数识别值变化趋势与理论值变化趋势基本一致，说明本文改进的SCAM算法具有较高的计算准确性。图8给出了不同噪声水平下各损伤状态的各层刚度参数极大后验估计值与理论值的相对误差。由图可知，无噪声下相对误差在±4%以内；5%噪声下相对误差在±5%以内,10%噪声下相对误差在±6%以内。本文所提方法不仅对中等损伤和大损伤(如损伤30%，50%)识别效果较好，而且对小损伤(损伤5%)的识别效果同样可以接受。从整体上看，随着噪声水平的增加，识别结果相对误差存在增大的趋势，但在土木工程中仍在可接受范围内。

图8 各损伤状态下各层刚度参数相对误差Fig.8 Relative errors of the stiffness parameters of each story foreach damage case

图7 各层刚度参数理论值与最优估计值变化趋势对比Fig.7 Comparison of the variation trends between the theoretical value and the optimal value of the stiffness parameters of each story

图9给出了各损伤状态下各层刚度后验变异系数变化趋势，变异系数越小表示识别结果不确定性越小，其可信度越高。在三种噪声水平下，底层刚度后验变异系数普遍小于其他三层的刚度后验变异系数，但各层刚度后验变异系数均小于先验变异系数0.2，表明改进方法有效性的同时，也表明基于贝叶斯估计结构物理参数识别方法确实可以降低参数的不确定性。

图9 各损伤状态下各层刚度参数变异系数变化趋势Fig.9 Variation trend of stiffness parameters variation coefficient of each story foreach damage case

该算例表明，基于改进SCAM算法结构物理参数识别贝叶斯方法具有较高的识别精度，能够准确地识别结构物理参数及其变化,且具有较好的抗噪能力，验证了改进的SCAM算法有效性和可靠性。

3 结 论

本文通过对SCAM算法中提议分布方差计算方法进行重新定义，提出了改进的SCAM算法，并通过理论分析和数值模拟算例对本文方法的有效性和可靠性进行了验证。所得结果和结论如下：

(1)通过重新定义提议分布方差所改进的SCAM算法，抽样样本序列所构成的马尔可夫链达到了平稳分布，且没有明显的趋势性和周期性，在提高计算准确性的同时提高了计算效率。

(2)基于改进SCAM算法的结构物理参数识别贝叶斯方法可以较好地识别结构物理参数，得到物理参数的后验边缘概率分布和最优估计值，识别精度及效果良好，且具有较好的抗噪能力，该方法可以应用于实际工程结构物理参数识别及损伤评估。

本文改进的SCAM算法仅在基本理论、数值模拟中进行了验证，作者下一步将对实际工程结构及结构振动台试验模型进行分析，进一步验证方法的可靠性、有效性及实用性。