基于有限测量信息的两步响应重构方法

2022-06-17史鹏程彭珍瑞董康立

振动与冲击 2022年11期

史鹏程，彭珍瑞，董康立

(1. 兰州交通大学机电工程学院，兰州 730070； 2. 浙江大学生物医学工程与仪器科学学院，杭州 310027)

随着我国经济的飞速发展，近年来国内大跨度桥梁以及高层建筑的数量日益增加，结构各处的响应数据是对其进行实时健康监测的重要基础。然而这些结构的服役时间往往长达几十年甚至上百年，远远超过测量其动态响应的传感器的寿命。因此基于有限的测量信息，对结构其余各处响应进行重构的结构响应重构技术具有重要的意义[1-4]。

现有的响应重构技术大致可分为两类：第一类基于模态分析或传递率的思想，建立已知测量响应和需要重构响应之间的传递关系，从而获得所需的各点的动态响应。Kammer[5]提出基于模态分析的响应重构方法，利用可测量的响应数据和结构模态振型矩阵来重构不可测位置处的动态响应。Zhang等[6-7]采用这种方法，针对梁模型研究了双重传感器下的结构应变、位移响应重构。Ribeiro等[8]最早提出传递率的概念，并将之应用于单自由度系统及多自由度系统的位移响应重构。Law等[9]将位移传递率扩展到了加速度传递率，并将这种方法应用到子结构的响应重构中，在响应计算中结合一阶保持输入逼近，提高了动态响应的准确性。Li等[10]又将基于传递率的响应重构方法从频域扩展到了小波域中，采用离散小波变换对结构和子结构系统分别进行响应重构计算，精度要高于离散傅里叶变换响应重构的精度。此类方法需基于精准的有限元模型，且没有考虑传感器测量噪声对重构精度的影响。

第二类基于卡尔曼滤波(Kalman filter, KF)算法的结构响应重构方法。张笑华等[11]基于KF算法对二维桁架的应变和位移响应进行了重构，同时对位移和应变两种传感器的数量和位置进行了优化。Xu等[12]使用KF对大跨度悬索桥进行应变、位移、加速度的响应重构，并结合逐步消去法对多类型传感器进行优化布置。任鹏等[13]分别使用增广卡尔曼滤波法和传递率法在未知外部激励的情况下对一钢桁架进行了应变响应重构，并对比了两种方法的优劣。董康立等[14]基于激励计算卡尔曼滤波算法验证了对外部激励和各处响应进行重构的有效性，并以响应重构精度为目标，使用萤火虫算法和逐步消去法分别进行传感器优化布置。此类方法是一种不确定性方法，可以有效地降低模型误差和测量噪声对于响应重构带来的影响。但使用KF算法必须事先假设测量噪声为高斯白噪声，而实际工程中的噪声却多为有色噪声，若把有色噪声当作白噪声处理，必然影响响应重构精度。粒子滤波算法可以很好地处理有色噪声的影响，但使用传统粒子滤波算法进行响应重构，会出现粒子贫化现象，从而影响响应重构精度。近年来不少学者将群体智能优化算法引入粒子滤波中[15-16]，替代原有的重采样过程，皆取得了较好的效果。

目前，对各类响应重构方法的研究通常也伴随着以提高响应重构精度为目的的传感器优化布置方案的研究[17-20]，但是在结构的服役过程中传感器的损坏往往是随机的，不可能符合优化布置中的方案。

因此，针对以上问题，首先基于已知测点的应变响应和位移响应，通过建立已知测量响应和所需重构响应的关系，重构出所需重构点的应变和位移响应。针对有色噪声问题，引入粒子滤波算法，同时结合萤火虫算法替代传统重采样过程，改善粒子贫化问题。然后结合改进的粒子滤波算法，减小模型误差和测量噪声对于重构结果的影响，得到最终的响应重构结果，实现利用正常测量的响应对因传感器损坏而丢失或其余未布置传感器处的响应的重构。最后以一个二维桁架结构为例，验证所提方法的有效性。

1 结构响应重构

在线性结构中，结构的应变响应和位移响应可以表示为

(1)

式中：ε为结构应变响应；d为结构位移响应；Ψs为所选择的前s阶应变振型；Φs为所选择的前s阶位移振型；Γs为Ψs和Φs组成的矩阵；q为模态坐标向量。

测量位置的动态位移响应和应变响应同样也可表示为

(2)

应变振型矩阵和位移振型矩阵的关系可表示为

Ψ=TΦ

(3)

式中：T为应变位移转换矩阵，与有限元模型的形函数有关。

结合式(1)和式(2)，可重构出所需位置的动态响应

(4)

(5)

(6)

式中：σε和σd分别为应变和位移的测量噪声标准差。重构式(4)可重新写为

(7)

2 基于状态空间的滤波处理

结构的二阶运动方程可以用模态坐标表示为

(8)

式中：q为模态坐标；ξ为阻尼矩阵；ω0为模态频率矩阵；u为外部激励向量；L为激励映射矩阵；Φ为模态位移振型矩阵。

将式(8)转化为状态空间方程并离散化

(9)

式中：xk，yk和uk分别为离散后的状态向量，观测向量和外部激励向量；A和B分别为离散化后的状态矩阵和输入矩阵；C和D分别为输出矩阵和直接传输矩阵；wk和vk分别为由于模型的不确定性造成的系统噪声和测量误差造成的测量噪声。式(9)中的各项具体表示为

(10)

(11)

(12)

(13)

(14)

通过建立以上状态空间模型，可以将式(7)得到的响应作为状态空间中的观测值，并通过滤波处理得到结构真实响应：

(15)

式中：Ce和De分别为结构对应位置输出矩阵和直接传输矩阵，是矩阵C和D的子矩阵。

3 粒子滤波算法

粒子滤波是一种基于状态空间模型，结合贝叶斯估计和抽样理论的滤波算法。其基本思想是利用一组样本对条件后验密度函数进行近似，用样本均值取代原有的积分运算，获得对系统的最小方差估计。系统状态空间模型可表示为

xk=f(xk-1,wk-1)

(16)

yk=h(xk,vk)

(17)

式中：xk为状态向量；wk-1为系统噪声；yk为观测向量；vk为测量噪声。

设状态初始概率密度p(x0|y0)=p(x0)，则预测方程为

(18)

状态更新方程为

(19)

其中

(20)

采样的重要性函数q(x0:k|y1:k)可以改写为

(21)

则权值公式为

(22)

(23)

式中：δ(·)为狄拉克函数。

概率密度更新公式为

(24)

将权值归一化后计算状态向量

(25)

(26)

以上为粒子滤波算法基本流程，在粒子集的不断迭代更新过程中，会使大多数粒子的权值接近0，少数粒子占有较大的权重，即粒子退化现象。使整个过程中的大量计算耗费在无意义的极小权值粒子上，针对这一问题，Gordon等[21]提出重采样来解决粒子退化问题，即复制大权值粒子，删除小权值粒子。但是传统重采样方法由于删除了小权值粒子，同时复制大权值粒子，反复重采样导致后代粒子多为采样前某几个大权值粒子的后代，从而减小了粒子的多样性，产生粒子贫化问题，导致算法整体性能的下降。

4 萤火虫算法

萤火虫算法是模仿萤火虫运动行为的一种智能算法。萤火虫个体之间的移动主要依赖于荧光亮度和吸引度两个参数，荧光亮度反映萤火虫所处位置的优劣并决定移动方向，吸引度与萤光亮度成正比，决定萤火虫的移动距离。萤火虫个体之间通过向周围萤光亮度更高的萤火虫移动实现寻优。

萤光亮度可以表示为

I=I0e-γr

(27)

式中：I0为初始萤光亮度；γ为荧光吸收系数；r为萤火虫个体之间的距离。

萤火虫之间的吸引度可以表示为

β=β0e-γr2

(28)

式中：β0为初始吸引度。

萤火虫个体之间的运动更新可以表示为

(29)

式中：xi和xj分别为两个萤火虫i和j的位置；α为步长因子。

5 结合萤火虫算法的粒子滤波

5.1 改进萤火虫算法

针对传统重采样方法带来的粒子贫化问题，将萤火虫算法引入粒子滤波中，替代传统的重采样方法，改善粒子贫化问题。由于萤火虫算法在寻优过程中单个萤火虫个体需要与其他所有的萤火虫个体之间计算相互之间的吸引度、萤光亮度等参数，并与其产生相对移动。若将其直接引入粒子滤波，将增加计算复杂度，影响滤波的实时性。同时，由于式(29)中添加随机扰动项避免算法陷入局部最优，但由式(10)可知状态向量中的参数并不在一个数量级，直接加入同一数量级的扰动项必然不可行。

针对以上问题，结合粒子滤波以及响应重构的需要，对于传统萤火虫算法作如下改进后再将其引入粒子滤波算法。

修正运动更新公式

(30)

式中：xbest为k时刻的全局最优值。改进后的运动更新方案使得萤火虫个体i在k时刻只需与全局最优值之间进行参数计算并进行相对移动，这大大降低了运算复杂度，提高了滤波实时性。

修正萤光亮度公式

(31)

5.2 算法步骤

(32)

步骤2将萤火虫算法中的萤火虫个体看作粒子滤波中的粒子，引入改进后的萤火虫算法指导粒子移动：根据式(28)计算粒子i与全局最优值xbest之间的吸引度,并依据式(30)指导粒子向全局最优值移动。

步骤3当所有粒子完成位置更新后，根据式(31)计算每个粒子的荧光度值，并更新最优粒子。

步骤4当荧光亮度小于阈值时，证明预测值与观测值差值已较小，粒子分布已接近真实分布，停止迭代，若未达到阈值则转入“步骤2”。或当迭代次数达到设置的最大迭代次数时，亦停止迭代过程。

步骤5依据式(24)～式(26)计算粒子权重并归一化后计算最终状态输出。

为避免在使用萤火虫算法指导粒子进行移动时迭代次数过多而使得粒子过于收敛，导致粒子多样性降低，但同时又要保证粒子对高似然区域有着较好的覆盖能力，能反应出系统真实的状态。为此，在“步骤4”中设置算法的最大迭代次数和荧光度阈值控制其迭代次数。

6 数值算例

采用二维桁架模拟整个响应重构过程，以验证所提方法的有效性。如图1所示。该桁架结构共有14个节点，25个单元和25个自由度。弹性模量、密度和泊松比分别为210 GPa、7 850 kg/m3和0.3。选取前六阶模态振型为目标模态进行响应重构，模态频率分别为17.72 Hz、58.89 Hz、85.48 Hz、113.67 Hz、169.98 Hz、220.22 Hz。在第4节点的Y方向施加随机白噪声激励，激励的时程曲线如图2所示。提取10、11、12单元的应变响应和4、5、7节点的Y方向的位移响应，并加入3%的有色噪声作为正常工作的传感器的测量响应，并将结构其余单元、节点未加噪声的响应yr作为其余所需重构位置以及传感器失效位置的响应。有色噪声的功率谱如图3所示。由图3可知，有色噪声的功率谱分布是不均匀的。

图1 二维桁架结构Fig.1 Two-dimensional truss structure

图2 随机白噪声激励Fig.2 Excitation of random white noise

图3 有色噪声功率谱Fig.3 Power spectrum of colored noise

将提出的两步响应重构方法与文献[11]中的响应重构方法进行对比，由于文献[11]中的KF算法不适合处理有色噪声问题，因此将改进粒子滤波算法替换原有的KF算法直接进行响应重构。在改进的粒子滤波算法中，设置最大迭代次数为15次，荧光度向量求模后设置阈值为10。

在结合萤火虫算法的粒子滤波中，由于粒子维数较大，不容易直观地看出是否发生粒子贫化或粒子退化现象，因此通过粒子权重，间接反映粒子质量。某一时刻的粒子权重，如图4所示。从图4可知，粒子权重分布均匀，不存在粒子退化或粒子贫化问题。

图4 粒子权重Fig.4 Weight of particles

方法1本文提出的分两步进行响应重构方法，首先对各类所需重构响应进行直接计算，第二步再使用粒子滤波进行滤波处理。

方法2基于粒子滤波算法通过测量响应直接进行响应重构。

设定9单元、19单元的应变传感器以及9节点、11节点Y方向的位移传感器失效，通过现有的10单元、11单元、12单元的应变传感器和4节点、5节点、7节点的Y方向的位移传感器测得的响应，对传感器失效处的响应以及其他所有各点的响应进行重构。

两种方法分别对传感器失效处的应变和位移响应重构的结果，如图5和图6所示。从图可知，两种方法重构的响应曲线均能较好地拟合测量响应曲线，均能对于缺失响应进行重构。

(a) 9单元应变

(b) 19单元应变图5 重构应变响应与测量响应对比Fig.5 Comparison of reconstructed and measured strain responses

(a) 9节点Y向位移

(b) 11节点Y向位移图6 重构位移响应与测量响应对比Fig.6 Comparison of reconstructed and measured displacement responses

为比较两种方法在重构精度上的差异，引入相对百分比误差(relative percentage errors, RPE)来衡量精度，具体公式如下：

(33)

两种重构方法对于传感器缺失处的响应重构结果的RPE如表1所示，对比结果可以看出方法1的重构精度要更高。

表1 相对百分比误差Tab.1 Relative percentage errors

用RPE比较两种方法对于其余所有位置响应的重构结果精度，结果如图7所示。从图7可知，方法1的应变响应重构的RPE基本都在0.4%以内，平均RPE为0.36%；位移响应重构的RPE基本都在0.5%以内，平均RPE为0.44%；方法2的应变响应重构的RPE基本在0.7%以内，平均RPE为0.48%；位移响应重构的RPE基本在0.8%以内，平均RPE为0.61%。方法1与方法2相比，由于方法1第一步已经计算出了含噪声的各类响应，再进行粒子滤波处理时增加了状态空间方程的观测量，因此相比方法2，其重构精度更高。

图7 相对百分比误差Fig.7 Relative percentage errors

综上分析可以看出，两步响应重构方法可以实现对于传感器失效处的响应和其余未安装传感器处的响应进行有效地重构，相比于直接使用测量响应通过滤波算法进行重构的方法重构精度更高。

为了探究在传感器位置发生改变时，两步响应重构方法精度是否会下降，即对测点的随机变化是否具有一定的鲁棒性，将响应测量位置改为9单元、13单元的应变响应和6节点、8节点、10节点、12节点的Y方向的位移响应，失效传感器为10单元的应变传感器和11节点Y方向的位移传感器。再次进行响应重构并对比重构精度，更改测点后此方法仍能较好地对因传感器损坏而丢失的响应进行重构，如图8所示。对于其余所有位置的响应的重构精度并未发生较大变化，验证了所提方法的鲁棒性，如图9所示。