基于移动窗卡尔曼滤波算法的结构响应重构

2021-11-17张笑华吴志彪吴圣斌黄梅萍

振动与冲击 2021年21期

张笑华，吴志彪，吴圣斌，黄梅萍

(福州大学土木工程学院，福州 350108)

结构健康监测的首要关键工作就是通过传感器系统正确地测量结构的动力响应。然而由于经济条件、仪器设备和现场环境等因素的制约，安装在结构上的传感器数目始终有限，许多位置并未安装传感器，从而无法为结构健康监测提供完整有效的信息。因此，利用有限位置上的测量响应来重构未测量位置的结构响应具有重要意义。

目前，已有的响应重构方法主要是基于传递的概念，即通过在时域、频域或者小波域里构造传递矩阵和传递函数来实现响应重构[1-3]。这些方法主要考虑了单种响应的重构。已有的研究结果表明：用多类型的测量信息来识别结构参数，其精度高于用单一信息识别的结果，更有利于结构分析、监测和安全评估[4]。因此，Zhang等[5]基于模态理论研究了位移和应变响应重构的方法。Zhu等[6]扩展了该方法，用经典卡尔曼滤波(Kalman filter, KF)算法，考虑传感器位置的优化，重构未测点加速度、位移和应变三种响应。

KF算法是一种以最小均方误差为最佳准则，利用测量响应对系统状态向量进行最优估计的迭代算法。因其优越的性能，该算法也被广泛应用于土木工程领域中，包括线性和非线性结构参数识别[7-8]，结构损伤识别[9]，有限元模型修正[10]等。近年来，也有不少学者将KF算法应用于结构响应的重构。Papadimitriou等[11]利用KF算法估计应力功力谱，用于识别因疲劳而引起的结构损伤累积。Jo等[12]融合加速度和应变测量响应，利用KF算法估计结构的位移响应。Zhang等[13]基于扩展KF算法，通过有限的传感器测量响应，首先估计外部激励，然后重构未测量位置的结构响应。万志敏[14]基于贝叶斯理论，研究了类卡尔曼滤波算法的外激励识别和结构响应重构方法。董康立等[15]分别将逐步消去法和萤火虫算法与类卡尔曼滤波算法相结合，提出了面向多类型传感器优化布置的结构响应重构方法。

目前大多数基于KF算法的结构响应重构方法，需预先假定测量与过程噪声方差已知且为常数。然而，大部分情况下，它们是未知且是时变的。不精确的噪声方差可能导致较大误差的状态向量估计，甚至可能导致滤波发散。Fuad等[16-17]研究了过程和测量噪声方差未知情况下的KF算法。Lai等[18]利用移动窗技术，通过测量响应和系统实际输出之间的残差，实时估计过程噪声方差和测量噪声方差，并用于KF算法以识别结构刚度的变化。Zhang等[19]改进了Lai等的测量噪声估计方法，研究了测量噪声方差未知下的结构响应重构，但该方法需先假定过程噪声方差已知。

基于以上背景，考虑测量和过程噪声方差均未知的情况，研究有限测点下基于移动窗卡尔曼滤波(moving-window Kalman filter, MWKF)算法的结构响应重构。利用有限测点的位移和应变测量响应，先采用移动窗技术估计结构的测量与过程噪声方差，然后基于KF算法重构结构中未安装传感器位置处位移与应变响应，最后以一个平面单跨框架结构为例进行数值模拟和试验分析，验证所提方法的有效性和可行性。

1 MWKF算法

线性系统的结构动力学方程可写成连续的状态空间方程。由于实测数据都是离散的，且实际工程的信号测量中，总是存在系统的不确定性，若将这些不确定性分为测量噪声和过程噪声，则连续的状态空间方程经过离散后，可得到下述的离散时间随机状态空间方程

zk+1=Azk+Buk+wk

(1)

yk=Czk+Duk+vk

(2)

其中：

(3)

(4)

(5)

式中：下标k表示第k时刻；z为离散状态向量；q是模态坐标向量，y是离散观测向量；A、B分别是离散状态矩阵和外部激励的输入矩阵；C，D分别代表测量输出矩阵和传递矩阵。Δt是采样时间间隔；I是单位矩阵；u是外部激励向量；w表示过程噪声向量；v表示测量噪声向量；w和v假定为均值为零且互不相关的高斯白噪声，它们的协方差矩阵分别为Q和R。Φh是质量归一化后的位移振型矩阵；Bu表示激励位置矩阵；ξ是模态阻尼系数矩阵；ω0是模态频率矩阵。

KF算法是无偏递归算法，可最优估计未知状态向量。经典KF算法为

(6)

(7)

先验误差方差Pk|k-1

Pk|k-1=APk-1AT+Q

(8)

状态向量估计均方差Pk

Pk=[I-KkC]Pk|k-1

(9)

卡尔曼滤波增益Kk

Kk=Pk|k-1CT[CPk|k-1CT+R]-1

(10)

算法中的R和Q通常需要通过大量实际实验的统计特性或在算法实际操作之前凭经验评估。在土木工程领域，它们常被假定为常数且为恒定。然而，事实上它们经常是时变的，尤其是在结构长期的运营服役期间。因此，本文基于移动窗技术，首先估计噪声方差阵R和Q。

定义测量噪声σk为实际响应与测量响应之间的差值。由于实际响应难以获得，这里用测量响应的加权平均来代替实际响应，具体计算如下

yk，smooth=λ1yk-1，m+λ2yk，m+λ3yk+1，m

(11)

式中：λ1,λ2,λ3为加权系数；下标‘m’表示测量的。用yk，smooth代替实际响应计算测量噪声

σk=yk，m-yk，smooth

(12)

(13)

其中

(14)

式中，Nr表示移动窗长度。

(15)

其中γi为加权系数，计算如下

(16)

式中，tr(·)表示计算矩阵的迹。测量噪声方差的估计仅需要通过测量数据来获得而不涉及任何结构特性。

测量残值即新息(innovation)ek定义为测量估计值与实际测量值之差，它表示存在于测量响应中而滤波器无法预测的那部分信息，即

(17)

其对应的自相关系数Hτ定义如下

(18)

式中，Nq为移动窗长度。测量残差的自相关系Hτ与先验误差方差之间有如下关系

当τ=0时

Hτ=CPk|k-1CT+R

(19)

当τ>0时

Hτ=

CA(I-KkC)τ-1A[Pk|k-1CT-Kk(CPk|k-1CT+R)]

(20)

(21)

式中：a为控制收敛速率的系数矩阵；H0为测量残值的方差矩阵。

2 结构位移和应变响应重构

考虑应变与位移两种响应，观测向量为

(22)

因此，式(2)中的C和D矩阵可写为

(23)

式中:Ψh=BsdΦh为应变模态矩阵;Bsd是应变-位移转换矩阵，与有限元划分的单元类型形函数有关。这里将输出矩阵C分为Cm和Ce两部分，分别对应于测量的位置和需响应重构的位置。

由于同时考虑了应变和位移两种响应，二者在数值上有量级差异，引起对应输出矩阵C的高度病态，导致卡尔曼增益矩阵计算结果不准确。可采用应变和位移对应的测量噪声标准差来规则化输出矩阵Cm

(24)

其中，

(25)

(26)

未测点的响应时程按下式计算

(27)

实际上，并非所有的振型对结构的整体动力响应具有同等重要的贡献，因此在使用基于MWKF算法重构未测点结构响应时，可根据实际情况仅考虑对动态响应贡献最大的前几模态振型，忽略高阶模态的影响。

3 数值模拟分析

采用三层单跨框架为数值算例验证所提方法的有效性和可靠性。该框架每层高度为0.6 m，跨度为0.6 m，梁与柱的横截面尺寸分别为50 mm×8.8 mm和50 mm×5 mm，弹性模量为206 GPa，密度为7 850 kg/m3，泊松比为0.3。同时也建立了相应的有限元模型，模型共有53个节点和54个单元，其中节点1和节点20为固定端，模型共有153个自由度，如图1所示。根据模态振型对减小响应重构误差的贡献情况，选取了前6阶模态为目标模态进行未测点的响应重构[20]。前6阶的模态频率分别为3.78 Hz、11.21 Hz、17.17 Hz、46.06 Hz、54.48 Hz和66.17 Hz。根据文献[21]的方法得到传感器位置，如图1所示。

外荷载以随机激励为例，利用状态空间方程并考虑过程噪声的影响计算结构动力时程响应。该随机激励的频带宽度为0.5～50 Hz，最大幅值为100 N，施加于节点38-x向。假设应变对应的过程噪声方差阵为Qs=1×10-15·I，位移对应的过程噪声方差阵Qd=1×100·I。提取对应有安装传感器位置处的动力响应并人为加入3%白噪声后作为测量响应。

图1 框架结构传感器布置Fig.1 Sensor placement of the frame model

图2是测量噪声方差时程曲线。图2(a)对应于节点19处位移测量响应，图2(b)对应于单元7处应变测量响应。由图2可知，估计的测量噪声随着时间的推移在整个过程中不断变化且会在某一值上下波动。由于用测量响应的加权平均值代替了实际响应，因此实际的测量噪声方差和估计的测量噪声方差均值之间稍有差异。图3为节点19处位移响应对应的过程噪声方差时程与单元7处的应变响应对应的过程噪声方差时程。由图3可知，估计的过程噪声方差也随着时间的变化而变化，应变对应的过程噪声方差趋于1×10-15，位移对应的过程噪声方差趋于1×100。不论是应变还是位移，过程噪声方差都趋于预设值。

为了评估测量和过程噪声未知条件下，基于KF算法重构响应的精度，定义了相对百分比误差，计算如下

(a) 节点19

(b) 单元7图2 测量噪声方差Fig.2 Measurement noise covariance

(a) 节点19

(b) 单元7图3 过程噪声方差Fig.3 Process noise covariance

(28)

式中:std指标准差;下标j表示的是第j个单元或节点。

图4显示了使用KF算法(R与Q均已知)与使用MWKF算法(R与Q均未知)重构未测点响应得到的相对百分比误差对比。其中，R与Q均已知时，位移响应重构的相对百分比误差大部分都在2%以内，平均相对百分比误差0.81%；应变响应重构的相对百分比误差大部分都在8%以内，平均相对百分比误差为2.08%；R与Q均未知时，位移响应重构的相对百分比误差在3%以内，平均相对百分比误差为1.49%；应变重构的相对百分比误差都在10%以内，平均相对百分比误差为2.93%。两种情况得到的响应重构结果精度相差不大。图5为使用MWKF算法得到的单元4应变重构响应时程、节点19-x向位移重构响应时程与相应计算响应时程的对比。由图5可知，重构响应时程与有限元软件计算得到的响应时程吻合良好。数值算例的结果表明，在测量和过程噪声方差未知的情况下，基于MWKF算法，可对测量和过程噪声方差进行实时估计，并用于结构未测点的响应重构，能取得较好的重构精度。

(a) 应变

(b) 位移图4 相对百分误差Fig.4 Relative percentage errors

(a) 单元4应变响应

(b) 节点19-x向位移响应图5 计算响应与重构响应时程对比Fig.5 Comparison of reconstructed and real time history responses

4 实验室框架试验验证

根据数值算例中框架结构的几何尺寸和物理参数，在实验室制作了相应的单跨三层钢框架，模型如图6所示。框架的材料选用Q235钢，框架底座通过两个L型角钢锚固在固定于地面的工字型钢梁底座上。本次试验中应变片的型号为BX120-5AA；位移传感器采样sensorpart FT 25 RA型激光位移计和KEYENCE IL-300型激光位移传感器；数据采集系统的型号为JM5959采集系统。

图6 实验室框架结构Fig.6 Experimental model

文章提出的基于MWKF算法的结构响应重构方法是基于有限元模型计算的振型向量。因此，首先须对有限元模型进行修正，确保有限元模型与实际结构的动力特性基本一致。为此，本次试验首先进行了结构模态试验，以获得结构的模态频率和振型。表1列出了前6阶的频率值和振型的MAC(modal assurance criterion)值。从表中可以看到，频率的测量值与计算值对比，最大误差为1.59%，其余误差基本都低于或者接近1%。前6阶模态振型的MAC值均与1.0非常接近。这些对比结果表明有限元模型与实验室模型的动力特性有很高的相似度，因此可以用有限元模型来预测实际结构的动态特性，不需要进行有限元模型修正。

表1 前6阶频率值及MAC值

用于响应重构试验的传感器位置如图7所示。除了12个测量点外，在节点11和26额外安装了2个激光位移计，在单元4、18、28、39和47增加贴了5个应变片，用于检验重构响应的精度。本文提出的响应重构方法与外荷载的类型无关。本试验施加的外荷载以锤击激励为例，施加在节点35-x位置上。为了验证不同的采样频率对重构方法的影响，本试验考虑了两种情况下的响应重构：

工况1：采样频率200 Hz；

工况2：采样频率500 Hz。

响应重构时考虑前6阶模态，加权因子、移动窗口大小与数值模拟分析中的取值相同。采集的动力响应首先用70 Hz的截止频率进行滤波，去除高频响应的影响，然后结合测量的位移和应变响应，对测量和过程噪声进行实时估计，并用于结构未测点的位移和应变响应重构。

图7 框架结构传感器布置Fig.7 Sensor placement of frame structure

将测量响应替代式(28)中的计算响应y，计算相对百分比误差。表2列出了两个工况的相对百分比误差。工况1中，应变和位移响应重构最大相对百分比误差分别为14.72%和9.75%；工况2中，应变和位移响应重构最大相对百分比误差分别为15.33%和11.15%。两个工况的相对百分比误差大部分都低于15%，采样频率不论是为500 Hz还是200 Hz，二者重构结果精度相差不大。图8和图9分别为工况1估计的测量噪声方差时程和估计的过程噪声方差时程，对应于单元7的应变和节点32的位移。图10和图11分别为工况1在时域和频域内单元4的应变重构响应、节点11-x位移重构与测量响应之间的对比，无论是时域内还是频域内，重构响应与测量响应之间的误差都很小，吻合良好。