最优过程噪声协方差Q估计下的结构响应重构

Structural Response Reconstruction under Optimal Process Noise Covariance Q Estimation

摘要

为解决卡尔曼滤波器在结构响应重构时过程噪声协方差难以获取的问题，提出两阶段最优过程噪声协方差Q估计（two-stage optimal Q estimation， TOQE）算法. 在固定窗长下建立关于Q的最小化误差函数，基于新息理论推导搜索区间，利用一维两阶段搜索策略加快过程噪声协方差Q的搜索，并使用SG（Savitzky-Golay， SG）平滑算法对窗口之间估计的过程噪声协方差进行平滑处理，实现最优过程噪声协方差的估计，使用卡尔曼滤波器实现结构的响应重构. 通过垂直轴风力发电机仿真和外伸梁实验验证算法的有效性，并应用重构的加速度响应在L1/L2正则化方法下实现外伸梁结构的激励识别. 结果表明，相较应用传统最优过程噪声协方差Q估计（optimal Q estimation， OQE）算法进行响应重构，TOQE算法在5%、10%、15%测量噪声下均能确定一维搜索的区间范围，避免了OQE算法搜索区间的主观选取，实现了最优过程噪声协方差的估计，且搜索效率提高约50%~60%. 不同精度评价标准下的TOQE算法响应重构误差有所降低，应用重构的加速度响应能较好识别外伸梁结构的外部激励. TOQE算法表现出较好的实时性和鲁棒性.

Abstract

To address the challenge of obtaining process noise covariance in Kalman filters for structural response reconstruction， a two-stage optimal Q estimation （TOQE） algorithm was proposed. Under a fixed window length， a minimization error function concerning Q was established， and the search interval was derived based on the innovation theory. A one-dimensional two-stage search strategy was employed to accelerate the search for the process noise covariance Q. The Savitzky-Golay （SG） smoothing algorithm was used to smooth the estimated process noise covariance between windows， achieving the optimal process noise covariance estimation. Finally， the Kalman filters was used to reconstruct the structural response. The effectiveness of the algorithm was validated through simulations of a vertical axis wind turbine and experiments on a cantilever beam. The reconstructed acceleration response was applied to identify the excitation of the cantilever beam structure under L1/L2 regularization. The results show that compared to the traditional optimal Q estimation （OQE） algorithm for response reconstruction， the TOQE algorithm can determine the search interval for one-dimensional search under 5%， 10%， and 15% measurement noise， respectively， avoiding the subjective selection of OQE search interval and achieving optimal process noise covariance estimation. The search efficiency was improved by 50%~60%. Under different accuracy evaluation standards， the response reconstruction error of the TOQE algorithm is reduced， and the reconstructed acceleration response can better identify the external excitation of the cantilever beam structure. The TOQE algorithm demonstrates good real-time performance and robustness.

Graphical abstract

关键词

过程噪声协方差 / 响应重构 / Savitzky-Golay / 卡尔曼滤波 / 激励识别

Key words

process noise covariance / response reconstruction / Savitzky-Golay / Kalman filters / excitation identification

结构的振动响应是结构健康监测的重要信息，大多数结构健康监测配置中，结构振动响应的数据样本容量取决于有限的传感器数量^［1-5］. 在实际工程中，传感器的安装数量和安装位置受现实因素的制约以及布置的传感器在服役过程中失效都会导致数据样本容量的不足和结构关键位置响应数据的缺失^［6］. 因此，使用响应重构技术对结构关键位置、传感器失效位置的响应进行重构至关重要^［7］.

确定性响应重构方法依赖于精确的有限元模型和测量响应，从而建立已知响应与未知响应之间的关系，实现对未知响应的重构^［8-9］. 在实际工程中，精确的有限元模型难以获取，结构的过程噪声（模态分析试验产生的模型误差）严重影响响应重构的精度. 近年来，卡尔曼滤波（Kalman filters， KF）算法被引入响应重构中，用来处理结构的过程噪声和传感器测量噪声，以提高响应重构的精度. 然而，KF算法需要获取过程噪声和测量噪声协方差矩阵的先验知识^［10-11］.噪声的估计技术目前可分为三类：协方差匹配技术^［12］、相关性技术^［13］和贝叶斯技术^［14］. 协方差匹配技术利用移动窗内的状态样本来计算协方差矩阵. Zhang等^［15］改进移动窗并将其应用到响应重构中，但是主观选取的加权因子严重影响噪声方差的估计精度. 相关性技术和贝叶斯技术的参数设置严重影响噪声协方差的估计，计算成本过高^［16-17］. 这些方法描述了一般的自适应滤波方案，并基于寻找过程和测量噪声协方差矩阵的稳态估计，而不是当前的最优值. Vettori等^［18］根据以上问题提出了最优过程噪声协方差Q估计（optimal Q estimation， OQE）算法，较好地估计出风机叶片的过程噪声协方差，并对作用在风机叶片的激励进行识别. 但是，OQE算法在进行一维搜索时，搜索范围通过主观经验设置，搜索时间过长. 对于设置的窗口之间数据误差未进行处理，只考虑了传感器精度对于过程噪声协方差估计的影响.

目前，大多数响应重构方法主要对测量噪声进行处理^［15］，测量噪声通过传感器采集混叠在测量信号当中，测量噪声方差的估计可通过对采集信号的频域、时域处理获得^［11］. 而过程噪声表征模态分析试验过程中的模型误差，误差大小严重影响响应重构的精度，且误差的估计无法用信号处理的方式实现，响应重构对难以获取的过程噪声通常按经验设置为已知常数^［19-20］，这在实际工程中难以实现.

针对以上问题，本文提出两阶段最优过程噪声协方差Q估计（two-stage optimal Q estimation， TOQE）算法进行卡尔曼滤波响应重构. 推导新息和过程噪声协方差之间的传递关系，获取搜索区间；在设置的固定窗长内，获取KF处理后的观测向量误差最小化函数，并通过一维两阶段搜索策略进行过程噪声协方差估计，减少搜索时长；对窗口间的噪声协方差进行平滑处理，实现最优过程噪声协方差估计，并用KF进行结构的响应重构，获取未布置传感器位置的加速度响应数据. 通过仿真和实验验证算法的有效性，并将实验重构的加速度响应用于激励识别，以验证重构响应的准确性.

1 卡尔曼滤波响应重构

线性离散状态空间方程为：

\left\{\begin{array}{c}{\mathit{x}}_{k+\mathrm{1}}=\mathit{A}{\mathit{x}}_k+\mathit{B}{\mathit{u}}_k+{\mathit{w}}_k\\ {\mathit{y}}_k=\mathit{C}{\mathit{x}}_k+\mathit{D}{\mathit{u}}_k+{\mathit{v}}_k\end{array}\right.（1）

式中：下标k表示第k时刻；\mathit{x}和\mathit{y}分别为离散后系统的状态向量和观测向量；\mathit{u}为外部激励向量；\mathit{A}、\mathit{B}、\mathit{C}和\mathit{D}分别为离散化后的系统状态矩阵、输入矩阵、输出矩阵和直接传递矩阵；\mathit{w}和\mathit{v}分别为系统过程噪声和测量噪声，分别假设均值为0、协方差矩阵为{\mathit{Q}}_k=E\left[{\mathit{w}}_k{\mathit{w}}_k^{\mathrm{T}}\right]和{\mathit{R}}_k=E\left[{\mathit{v}}_k{\mathit{v}}_k^{\mathrm{T}}\right]的高斯白噪声，其中E\left[·\right]表示均值.

为处理过程噪声和传感器测量噪声，采用KF进行响应重构. KF算法基于状态空间方程，针对线性系统的状态进行最优估计，分为预测和更新两个步骤.

预测步：

{\mathit{x}}_{k|k-\mathrm{1}}=\mathit{A}{\mathit{x}}_{k-\mathrm{1}|k-\mathrm{1}}+\mathit{B}{\mathit{u}}_{k-\mathrm{1}}（2）

{\mathit{P}}_{k|k-\mathrm{1}}=\mathit{A}{\mathit{P}}_{k-\mathrm{1}|k-\mathrm{1}}{\mathit{A}}^{\mathrm{T}}+{\mathit{Q}}_k（3）

更新步：

{\mathit{K}}_k={\mathit{P}}_{k|k-\mathrm{1}}{\mathit{C}}^{\mathrm{T}}{\left(\mathit{C}{\mathit{P}}_{k|k-\mathrm{1}}{\mathit{C}}^{\mathrm{T}}+{\mathit{R}}_k\right)}^{-\mathrm{1}}（4）

{\mathit{y}}_{k|k-\mathrm{1}}=\mathit{C}{\mathit{x}}_{k|k-\mathrm{1}}+\mathit{D}{\mathit{u}}_k（5）

{\mathit{x}}_{k|k}={\mathit{x}}_{k|k-\mathrm{1}}+{\mathit{K}}_k\left({\mathit{y}}_k-{\mathit{y}}_{k|k-\mathrm{1}}\right)（6）

{\mathit{P}}_{k|k}={\mathit{P}}_{k|k-\mathrm{1}}-{\mathit{K}}_k\mathit{C}{\mathit{P}}_{k|k-\mathrm{1}}（7）

式中：{\mathit{x}}_{k|k-\mathrm{1}}为先验状态估计；{\mathit{P}}_{k|k-\mathrm{1}}为先验估计误差协方差矩阵；{\mathit{y}}_{k|k-\mathrm{1}}为先验状态估计对应的先验观测向量；{\mathit{K}}_k为卡尔曼增益矩阵；{\mathit{x}}_{k|k}为后验状态估计；{\mathit{P}}_{k|k}为后验估计误差协方差矩阵.

由式（1）可得测量方程和重构方程分别为：

\left\{\begin{array}{l}{\mathit{y}}_k^{\mathrm{m}}={\mathit{C}}^{\mathrm{m}}{\mathit{x}}_k+{\mathit{D}}^{\mathrm{m}}{\mathit{u}}_k\\ {\mathit{y}}_k^{\mathrm{s}}={\mathit{C}}^{\mathrm{s}}{\mathit{x}}_k+{\mathit{D}}^{\mathrm{s}}{\mathit{u}}_k\end{array}\right.（8）

式中：上标\mathrm{m}和\mathrm{s}分别表示测量位置和关心位置；{\mathit{C}}^{\mathrm{m}}和{\mathit{C}}^{\mathrm{s}}分别表示测量位置和关心位置的输出矩阵，是\mathit{C}的子矩阵；{\mathit{D}}^{\mathrm{m}}和{\mathit{D}}^{\mathrm{s}}分别表示测量位置和关心位置的直接传递矩阵，是\mathit{D}的子矩阵；{\mathit{y}}_k^{\mathrm{m}}和{\mathit{y}}_k^{\mathrm{s}}分别表示传感器位置的测量响应和关心位置的重构响应，其中包括加速度响应\mathit{a}、速度响应、位移响应和应变响应.

将KF滤波后的{\mathit{x}}_{k|k}代入重构方程获得未知测量位置的重构响应数据{\mathit{y}}_k^{\mathrm{s}}，代入测量方程获得观测向量滤波估计后的响应{\mathit{y}}_k^{\mathrm{m}\mathrm{s}}.

响应重构精度评价标准采用平均绝对百分比误差E_MAPE和相对百分比误差E_RPE ^［21］. 两个评价指标的计算式分别如下：

E_{\mathrm{M}\mathrm{A}\mathrm{P}\mathrm{E}}=\frac{\mathrm{1}}{L}{\displaystyle \sum _{k=\mathrm{1}}^L}\left|\frac{{\mathit{y}}_k^{\mathrm{s}}-{\mathit{y}}_k^{\mathrm{r}}}{{\mathit{y}}_k^{\mathrm{r}}}\right|\times \mathrm{100}\mathrm{\%}（9）

E_{\mathrm{R}\mathrm{P}\mathrm{E}}=\frac{‖{\mathit{y}}_k^{\mathrm{s}}-{\mathit{y}}_k^{\mathrm{r}}‖}{‖{\mathit{y}}_k^{\mathrm{r}}‖}\times \mathrm{100}\mathrm{\%}（10）

式中：L表示测量响应数据长度；{\mathit{y}}_k^{\mathrm{r}}为不含测量噪声的真实响应；\left|·\right|表示绝对值；‖·‖表示2范数.

2 两阶段最优过程噪声协方差Q估计

建立KF后由后验状态估计得到的响应{\mathit{y}}_k^{\mathrm{m}\mathrm{s}}和传感器测量响应{\mathit{y}}_k^{\mathrm{m}}的误差E^［18］为：

E_i=\sqrt[]{{\displaystyle \sum _{k=i\times N}^{(i+\mathrm{1})\times N}}\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}{\mathrm{n}}_Q{‖{\mathit{y}}_k^{\mathrm{m}}-{\mathit{y}}_k^{\mathrm{m}\mathrm{s}}‖}^{\mathrm{2}}}（11）

式中：i=\mathrm{1,2},\cdots ,L/N表示窗口顺序；N表示窗口长度；误差E_i的维度取决于传感器的数量.

通过一维搜索最优的过程噪声协方差Q_k来获取卡尔曼滤波器在窗口长度N内所需的过程噪声协方差矩阵{\mathit{Q}}_k：

{\mathit{Q}}_k={\left[\begin{array}{cccc}\mathrm{1}& \mathrm{0}& \cdots & \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \cdots & \mathrm{0}\\ ⋮& ⋮& & ⋮\\ \mathrm{0}& \mathrm{0}& \cdots & \mathrm{1}\end{array}\right]}_{n\times n}\times Q_k（12）

式中：n为模态截断数；Q_k为过程噪声协方差.

OQE算法依赖主观选取的一维搜索区间，使用小步长逐次搜索，导致区间范围不确定，搜索时间长.

考虑传感器测量响应带噪，定义新息{\mathit{\eta }}_k^［15］为：

\begin{array}{l}{\mathit{\eta }}_k={\mathit{y}}_k^{\mathrm{m}}+\mathit{v}-{\mathit{y}}_{k|k-\mathrm{1}}=\\ [{\mathit{C}}^{\mathrm{m}}(\mathit{A}{\mathit{x}}_{k-\mathrm{1}}+\mathit{B}{\mathit{u}}_{k-\mathrm{1}})+{\mathit{D}}^{\mathrm{m}}{\mathit{u}}_k]+{\mathit{v}}_k-\\ [{\mathit{C}}^{\mathrm{m}}(\mathit{A}{\mathit{x}}_{k-\mathrm{1}|k-\mathrm{1}}+\mathit{B}{\mathit{u}}_{k-\mathrm{1}}+{\mathit{w}}_k)+{\mathit{D}}^{\mathrm{m}}{\mathit{u}}_k]=\\ {\mathit{C}}^{\mathrm{m}}\mathit{A}({\mathit{x}}_{k-\mathrm{1}}-{\mathit{x}}_{k-\mathrm{1}|k-\mathrm{1}})+{\mathit{v}}_k-{\mathit{C}}^{\mathrm{m}}{\mathit{w}}_k=\\ {\mathit{C}}^{\mathrm{m}}\mathit{A}({\mathit{x}}_{k-\mathrm{1}}-{\mathit{x}}_{k-\mathrm{1}|k-\mathrm{1}})+\sqrt[]{{\mathit{R}}_k}-{\mathit{C}}^{\mathrm{m}}{\mathit{w}}_k\end{array}（13）

式中：{\mathit{x}}_{k-\mathrm{1}}表示第k-\mathrm{1}时刻系统的状态向量.由于{\mathit{x}}_{k-\mathrm{1}|k-\mathrm{1}}是{\mathit{x}}_{k-\mathrm{1}}的无偏估计， 式（13）可以改写为：

{\mathit{w}}_k={{\mathit{C}}^{\mathrm{m}}}^{+}({\mathit{\eta }}_k-\sqrt[]{R_k})（14）

式中： +表示伪逆；{\mathit{C}}^{\mathrm{m}}表示非方阵. 引入最小二乘法求解伪逆矩阵，式（14）改写为：

{\mathit{w}}_k=({{\mathit{C}}^{\mathrm{m}}}^{\mathrm{T}}{\mathit{C}}^{\mathrm{m}}{)}^{-\mathrm{1}}{{\mathit{C}}^{\mathrm{m}}}^{\mathrm{T}}({\mathit{\eta }}_k-\sqrt[]{{\mathit{R}}_k})（15）

通过式（15）计算Q_k，但由于伪逆矩阵和无偏估计项的影响，计算的Q_k会分布在一定区间范围，取其数量级最大最小值作为一维两阶段搜索区间.

\left[\mathrm{1}\times {\mathrm{10}}^{-q_{\mathrm{1}}},\mathrm{9}\times {\mathrm{10}}^{-q_{\mathrm{2}}}\right]（16）

式中：q_{\mathrm{1}}和q_{\mathrm{2}}分别为计算的最小数量级和最大数量级数值.

在以上一维两阶段搜索区间内，设置大步长\mathrm{1}\times {\mathrm{10}}^{-q_{\mathrm{1}}\le -a_{\mathrm{1}}\le -q_{\mathrm{2}}}，a_{\mathrm{1}}每次增加1个数量级，当获取到大步长下最优a_{\mathrm{1}\mathrm{m}\mathrm{a}\mathrm{x}}时，设置小步长，通过获取大步长下最优a_{\mathrm{1}\mathrm{m}\mathrm{a}\mathrm{x}}生成小步长搜索区间，b\times {\mathrm{10}}^{-(a_{\mathrm{1}\mathrm{m}\mathrm{a}\mathrm{x}}+\mathrm{1})\le -a_{\mathrm{2}}\le -(a_{\mathrm{1}\mathrm{m}\mathrm{a}\mathrm{x}}-\mathrm{1})}，a_{\mathrm{2}}每次增加1个数量级，b=\mathrm{1,2},\cdots ,\mathrm{9}依次循环取值，获得估计的过程噪声协方差b_{\mathrm{m}\mathrm{a}\mathrm{x}}\times {\mathrm{10}}^{-a_{\mathrm{2}\mathrm{m}\mathrm{a}\mathrm{x}}}. 通过这种两阶段搜索策略来加快Q_k的搜索速度.

Q_k的获取在窗口长度N内进行，在整个数据维度上产生的L/N个窗口之间会因为卡尔曼滤波器在每个窗口初始值的变化产生误差，进而造成Q_k估计的波动.

SG滤波是一种基于信号局部特性的多项式拟合方法，广泛应用于信号的平滑和去噪. SG滤波器使用多项式p(l)拟合信号^［22］：

p(l)={\displaystyle \sum _{m=\mathrm{0}}^g}{h}_m{l}^m（17）

式中：l为窗口大小；h为多项式系数；g为拟合阶数.

因此，采用SG滤波器对估计的Q_k进行平滑处理，可以减小窗口间的数据衔接误差，从而获得最优的Q_k.

通过不同数量传感器产生的误差E，得到不同数量的最优Q_k，将这些最优的Q_k进行平均处理，提供给KF进行响应重构.

响应重构流程如图1所示. 首先，每个窗口下，利用新息计算一维两阶段搜索区间，通过建立的滤波误差最小化函数E逐步搜索过程噪声协方差. 在一维两阶段搜索区间的限制下，过程噪声协方差的搜索通过大步长和小步长搜索完成，当所有窗口内的过程噪声协方差计算结束后，通过SG滤波对窗口之间的数据进行平滑处理，获取最优过程噪声协方差. 最后，利用KF实现未测量位置的加速度响应.

3 算 例

3.1 数值算例

垂直轴风力发电机易受垂向疲劳载荷的作用，导致风机疲劳断裂^［23］，采用某型号3 kW垂直轴海上风力发电机的叶片和风轮为仿真研究对象，结构如图2所示. 叶片厚度0.06 m，风轮宽度0.2 m、厚度 0.1 m，结构采用6061型铝合金材料，弹性模量为 71 GPa，密度为2 700 kg/m³，泊松比为0.33. 结构共 4 230个节点，只考虑每个节点的垂向自由度. 通过截断模态法使用结构前6阶模态振形建立状态空间方程进行响应重构^［24］，结构频率分别为7.3 Hz、 9.6 Hz、10.8 Hz、21.8 Hz、25.2 Hz、44.3 Hz，将建模后的模型作为真实模型. 为模拟风力引起的结构振动，在风机叶片顶部垂向施加周期激励，激励的频带宽度为20~75 Hz，激励最大值为20 N.

为模拟加速度传感器带噪测量响应，利用 式（18）在真实响应{\mathit{y}}_k^{\mathrm{r}}数据中添加测量噪声^［19］.

{\mathit{y}}_k=\lambda \cdot \mathrm{s}\mathrm{t}\mathrm{d}\left({\mathit{y}}_k^{\mathrm{r}}\right)\cdot \mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{n}（18）

式中：\lambda为噪声因子，用于控制测量噪声等级；std（·）为标准差；randn为正态分布的随机数.

采用NA、NB、NC、ND、NE、NF节点处的垂向加速度响应用于结构响应重构，验证所提TOQE算法的有效性，仿真运行环境为Intel（R） Core（TM） i7-6500U处理器，通过MATLAB 2018b进行算法编写. 为了进行方案对比和算法鲁棒性验证，设置两种方案.

方案1：通过NA节点和NB节点的加速度响应作为测量响应，分别使用OQE和TOQE算法估计最优Q_k进行KF响应重构，重构NC节点和ND节点加速度响应.

针对风力发电机模型设置过程噪声协方差，使用TOQE算法在不同等级的测量噪声下进行Q_k的估计. 为了验证TOQE算法一维两阶段搜索区间的正确性，设置3种过程噪声协方差，即Q^{\mathrm{1}}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}、Q^{\mathrm{2}}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{6}}和Q^{\mathrm{3}}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{8}}；设置λ分别为5%、10%和15% 3个等级的测量噪声，对应的测量噪声协方差分别为R^{\mathrm{1}}、R^{\mathrm{2}}和R^{\mathrm{3}}. 使用TOQE算法进行一维两阶段搜索区间计算，结果如表1所示.

由表1可得，在3个等级5%、10%、15%测量噪声下，设置的3种不同过程噪声协方差均处在一维两阶段搜索区间的数量级范围内.

设置过程噪声协方差为Q^{\mathrm{2}}，根据搜索区间［1×10^-12，9×10^-2］，在3个等级5%、10%、15%测量噪声下，使用OQE和TOQE算法进行过程噪声协方差搜索，估计最优Q_k. 传感器采集数据长度L=2 000，设置搜索窗口长度N=\mathrm{100}，在初始窗口长度和5%测量噪声下，过程噪声协方差搜索情况如图3所示.

由图3（a）可知，在5%的测量噪声下，在NA节点和NB节点估计的加速度测量响应误差曲线下，使用OQE算法搜索的最优Q_k均为6×10^-7. 由图3（b）可知，在NA节点和NB节点估计的加速度测量响应误差曲线下，使用TOQE算法进行一维两阶段大步长搜索，确定的Q_k均为1×10^-6，由此生成的TOQE算法小步长搜索区间为1×10^-7至9×10^-5. 由图3（c）可得，在生成的TOQE算法小步长搜索区间内，使用TOQE算法搜索的Q_k均为6×10^-7，且响应误差曲线均在一维两阶段搜索区间内出现最小值，从而验证了区间的正确性.

TOQE与OQE算法在两节点下搜索的Q_k相同，由图4可知，相比OQE算法，TOQE算法对各窗口之间的数据进行了平滑处理，使得估计的最优Q_k能够在窗口之间平滑过渡，减小数据衔接误差.

图5、图6分别为OQE和TOQE算法在10%和15%的测量噪声下最优过程噪声协方差估计.与图4对比可得，随着测量噪声的增大，窗口之间数据的衔接能力变差，估计的Q_k波动增加，但相对于OQE算法，所提TOQE算法在窗口之间平滑过渡，曲线波动平缓，更贴近真实噪声协方差.

方案1中过程噪声协方差搜索时长如由表2所示，由表2可知，在一维两阶段搜索区间内，TOQE算法能够快速搜索Q_k，相对于OQE算法，效率提高约50%~60%.

利用OQE和TOQE算法估计的最优Q_k进行KF响应重构，NC节点和ND节点加速度重构响应分别如图7和图8所示. 不同测量噪声下，应用OQE和TOQE算法估计的最优Q_k进行响应重构，NA节点和NB节点加速度响应都能够较好地重构NC节点和ND节点加速度响应. 加速度响应重构精度如表3所示，由表3可知，TOQE和OQE算法响应重构精度E_MAPE和E_RPE均较高，但TOQE算法响应重构误差均有所降低.

方案2：设置NE节点和NF节点的加速度响应作为测量响应，分别使用OQE和TOQE算法估计最优Q_k进行KF响应重构，重构NC节点和ND节点的加速度响应. 在与方案1相同的噪声条件下，使用TOQE算法进行一维两阶段搜索区间计算，在3个等级5%、10%、15%的测量噪声下，设置的3种不同过程噪声协方差均处在一维两阶段搜索区间的数量级范围内，且与方案1所得区间一致.

设置过程噪声协方差为Q^{\mathrm{2}}，在3个等级5%、10%、15%测量噪声下，使用OQE和TOQE算法进行Q_k搜索，估计最优Q_k. 结果分别如图9、图10、图11所示.

由图9~图11可知，OQE和TOQE算法在3个等级5%、10%、15%测量噪声下进行最优Q_k的估计，结论与方案1相同，随着噪声的增大，估计的Q_k波动增加，但相对于OQE算法，TOQE算法在窗口之间平滑过渡，曲线波动平缓，更贴近真实噪声协方差.

过程噪声协方差搜索时长如表4所示，由表4可知，在一维两阶段搜索区间，相对于OQE算法，TOQE算法在方案2下效率同样提高约50%~60%.

利用OQE和TOQE算法估计最优Q_k进行KF响应重构，重构NC节点和ND节点的加速度响应. TOQE算法响应重构的E_MAPE在1.0%~1.1%，E_RPE在1.5%~1.7%；OQE算法响应重构的E_MAPE在1.1%~1.2%，E_RPE在1.6%~1.8%. TOQE算法响应重构误差均有所降低.

两种方案下，TOQE算法计算的一维两阶段搜索区间均相同，搜索速度快，估计的最优Q_k均在真实值1×10^-6附近波动.使用TOQE算法估计的最优Q_k进行KF响应重构，所得加速度响应重构精度均较高. 相对OQE算法，所提TOQE算法均能够快速、准确地得到最优过程噪声协方差，并实现结构在未测量位置的加速度响应重构.

3.2 实验分析

在工程实际应用中，许多工程结构都可以简化为外伸梁形式. 为了进一步验证所提方法的有效性，以外伸梁作为实验对象，进行实验分析，外伸梁结构如图12所示.

外伸梁尺寸为2 000 mm\times100 mm\times10 mm（长\times宽\times高），将该外伸梁划分为20个单元，21个节点，只考虑每个节点的垂向自由度.

采用INV3062-C2（L）型网络分布式采集仪系统DASP-V11信号分析软件进行实验. 设置采样频率为2 000 Hz，通过MSC-3弹性力锤对外伸梁施加垂向激励，完成对外伸梁的加速度信号采集和模态分析工作. 通过模态截断法^［24］，使用实验测得结构前 5阶模态振型建立状态空间模型进行响应重构，结构振型频率分别为7.7 Hz、28.4 Hz、62.1 Hz、110.7 Hz、166.3 Hz，其中建模产生的模型误差（即过程噪声的真实值）无法获取. 传感器布置如图13所示，在第11节点施加垂向锤击激励，测量第4、7、10和15节点处的垂向加速度信号用于模态分析后的模型误差估计和结构的响应重构，以验证所提方法的可行性.

在封闭的室内使用传感器进行数据采集作为节点真实响应，为了评估测量噪声的影响，在真实响应中添加\lambda分别为5%、10%和15%的测量噪声，其对应的测量噪声协方差为R^{\mathrm{1}}、R^{\mathrm{2}}和R^{\mathrm{3}}^［19-20］.

方案3：将4节点和15节点的加速度响应作为测量响应，分别用OQE和TOQE算法估计的最优Q_k进行KF响应重构，重构7节点和10节点加速度响应.

TOQE算法分别在5%、10%、15%等级的测量噪声下，进行一维两阶段搜索区间计算，结果如表5所示.

从表5可知，分别在5%、10%、15%等级的测量噪声下，一维两阶段搜索区间均为［1×10^-20，9×10^-7］.

根据搜索区间［1×10^-20，9×10^-7］，分别使用OQE和TOQE算法进行Q_k搜索，估计最优Q_k. 传感器采集数据长度L=\mathrm{2}\mathrm{000}，设置搜索窗口长度N=\mathrm{100}，在初始窗口长度和5%测量噪声下，过程噪声协方差的搜索情况如图14所示.

在5%的测量噪声下，由图14（a）（b）可得，在TOQE算法计算的一维两阶段搜索区间内，4节点和15节点估计的加速度测量响应误差曲线均出现了最小值点，验证了计算的一维两阶段搜索区间的正确性. 由图14（a）可得，在4节点和15节点估计的加速度测量响应误差曲线下，使用OQE算法搜索的最优Q_k分别为\mathrm{4}\times {\mathrm{10}}^{-\mathrm{13}}和\mathrm{1}\times {\mathrm{10}}^{-\mathrm{12}}. 由图14（b）可得，在4节点和15节点估计的加速度测量响应误差曲线下，使用TOQE算法进行一维两阶段大步长搜索，确定的Q_k均为\mathrm{1}\times {\mathrm{10}}^{-\mathrm{12}}，由此生成的TOQE算法小步长搜索区间为\mathrm{1}\times {\mathrm{10}}^{-\mathrm{13}}至\mathrm{9}\times {\mathrm{10}}^{-\mathrm{11}}. 由图14（c）可得，在生成的TOQE算法小步长搜索区间内，使用TOQE算法搜索的Q_k分别为\mathrm{4}\times {\mathrm{10}}^{-\mathrm{13}}和\mathrm{1}\times {\mathrm{10}}^{-\mathrm{12}}. 与仿真不同，受不同节点传感器数据采集误差的影响，不同节点所得协方差表现出一定的差异性.

同一窗口长度内，TOQE与OQE算法在同一节点下搜索Q_k相同. 方案3中5%测量噪声-最优过程噪声协方差估计如图15所示. 由图15可得，相比OQE算法，TOQE算法对各窗口之间的数据进行了平滑处理，使得估计的最优Q_k能够在窗口之间平滑过渡，数据波动明显减弱.

图16、图17分别为方案3中OQE和TOQE算法在10%和15%测量噪声下最优的Q_k估计. 对比图15可得，3种测量噪声下，所得最优Q_k均在\mathrm{1}\times {\mathrm{10}}^{-\mathrm{12}}附近波动，与仿真算例得出的结论相同，随着测量噪声的增大，窗口之间数据的衔接能力变差，估计的最优Q_k波动增加，但相对于OQE算法，所提TOQE算法估计的最优Q_k在窗口之间平滑过渡，曲线波动平缓，其估计的Q_k符合经验认知.

方案3中过程噪声协方差搜索时长如表6所示.由表6可得，与仿真结论相同，在相同的一维两阶段搜索区间内，所提TOQE算法能够快速搜索Q_k，相对于OQE算法，效率提高约50%~60%.

利用OQE和TOQE算法估计的最优Q_k进行KF响应重构，7节点和10节点的加速度重构响应分别如图18和图19所示，不同测量噪声下，应用OQE和TOQE算法估计的最优Q_k进行响应重构，4节点和15节点加速度响应都能够较好地重构7节点和10节点加速度响应. 方案3中加速度响应重构精度如表7所示，由表7可知，在E_MAPE和E_RPE 两种精度评价标准下，TOQE和OQE算法响应重构精度均较高，相比OQE算法，TOQE算法精度均有所提高.

方案4：使用7节点和10节点的加速度响应作为测量响应，分别使用OQE和TOQE算法估计最优Q_k进行KF响应重构，重构4节点和15节点的加速度响应.

在相同噪声条件下，使用TOQE算法进行一维两阶段搜索区间计算. 在5%、10%、15%的测量噪声下，TOQE算法计算的一维搜索区间为［1×10^-20， 9×10^-7］，与方案3所得的区间相同.

根据搜索区间［1×10^-20，9×10^-7］，在5%、10%、15%的测量噪声下，分别使用OQE和TOQE算法进行Q_k搜索，估计最优Q_k. 结果分别如图20~图22所示.

由图20~图22可得，OQE和TOQE算法分别在5%、10%、15%的测量噪声下进行最优Q_k的估计，所得最优Q_k均在\mathrm{1}\times {\mathrm{10}}^{-\mathrm{12}}附近波动，且TOQE能够快速搜索Q_k，相对于OQE算法，效率提高约50%~60%，且与方案3得出的结论完全一致.

利用OQE和TOQE算法估计最优Q_k进行KF响应重构，重构4节点和15节点的加速度响应. TOQE算法响应重构的E_MAPE为4.5%~6.8%，E_RPE为6.8%~9.9%；OQE算法响应重构的E_MAPE为4.8%~7.2%，E_RPE为7.2%~10.4%. TOQE算法响应重构误差均有所降低.

由于模态试验所产生的模型误差真实值无法获得，无法将估计的最优过程噪声协方差与真实值进行对比，但是，在两种方案下，TOQE算法计算的一维两阶段搜索区间均相同，估计的最优过程噪声协方差均在\mathrm{1}\times {\mathrm{10}}^{-\mathrm{12}}附近波动，符合过程噪声协方差的经验认知. 且通过前后对照，在计算的搜索区间内，测量响应误差曲线均出现了最小值点，可验证搜索区间的正确性和算法估计的合理性，说明TOQE算法能够有效估计模态分析试验所产生的过程噪声协方差（模型误差）.

因此，在不同重构方案下，所提TOQE算法均能够快速、准确地得到最优过程噪声协方差，并实现结构未测量位置的加速度响应重构.

4 激励识别

为了进一步评价TOQE算法进行卡尔曼滤波响应重构的效果，应用外伸梁实验重构的加速度响应进行激励识别. 通过构建传递矩阵\mathit{H}，获得激励识别方程：

\mathit{u}=({\mathit{H}}^{\mathrm{T}}{\mathit{H})}^{-\mathrm{1}}{\mathit{H}}^{\mathrm{T}}{\mathit{y}}^{\mathrm{m}}（19）

将15%的噪声下重构的7节点和10节点加速度响应数据用于激励识别. 使用传统L1正则化方法和L2正则化方法进行激励识别^［25］，激励识别效果分别如图23和图24所示.

对识别激励进行E_RPE计算，使用L1正则化方法进行激励识别，重构响应识别的激励与真实激励的误差为20.7%，传感器采集的加速度响应识别的激励与真实激励的误差为20.5%.使用L2正则化方法进行激励识别，重构响应识别的激励与真实激励的误差为18.6%，传感器采集的加速度响应识别的激励与真实激励的误差为18.2%. 在15%的噪声数据下，TOQE算法重构的7节点和10节点加速度响应与传感器采集响应在激励识别效果上几乎一致，产生误差的原因是传递矩阵\mathit{H}存在的不确定性. 由此可得，所提方法TOQE算法重构的响应数据可进一步用于冲击载荷识别，具有一定的工程适用性.

5 结 论

1）建立了两阶段最优过程噪声协方差Q估计

的卡尔曼滤波响应重构方法. 仿真和实验结果表明，不同方案下，所提方法响应重构精度均较高，鲁棒性好.

2）在不同等级的测量噪声下，通过TOQE算法获得了不同过程噪声协方差的一维搜索区间，避免了OQE算法搜索区间的主观选取. 使用SG滤波对估计的过程噪声协方差进行平滑处理，获得最优过程噪声协方差估计，有效减缓了窗口之间的数据波动.

3）相比OQE算法，TOQE算法提出的两阶段搜索策略在多种方案下均较大限度提高一维搜索速度，实现过程噪声协方差的快速估计.

4）对TOQE算法重构的加速度响应进行应用，外伸梁锤击激励识别结果表明，重构数据和传感器数据具有相同的效力，响应重构数据能够较好地应用于激励识别.

与一般的自适应滤波方案无法通过自身对照来确定估计的噪声协方差的合理性不同，TOQE算法通过不同方案和算法的前后对照有效确定响应重构的过程噪声协方差，且没有初始参数设置的影响. 因此，TOQE算法能够有效地获取过程噪声协方差表征的模型误差，这对于难以获取的模型误差以及在实验室环境下模态试验的后续分析具有重要意义.

参考文献

[1] 周云，危俊杰，李剑，等 ．基于InSAR技术的大跨桥梁温度变形监测研究［J］．湖南大学学报（自然科学版），2024，51（3）：39-50．

[2] ZHOU Y， WEI J J， LI J，et al ．Research on temperature deformation monitoring of long-span bridges based on InSAR technology ［J］．Journal of Hunan University （Natural Sciences），2024，51（3）：39-50．（in Chinese）

[3] THANKACHAN P， FIDA A， MADHAVAN PILLAI T M ．Health monitoring of steel structures using Cuckoo Search algorithm-based ANN［J］．Structures，2024，61：105933．

[4] 封周权， 王文赞， 华旭刚， 等 ．基于模态参数与改进萤火虫算法的结构模型修正［J］．湖南大学学报（自然科学版）， 2022， 49（11）：252-259．

[5] FEN Z Q， WANG W Z， HUA X G，et al. Structural model updating based on modal parameters and modified firefly algorithm［J］．Journal of Hunan University （Natural Sciences）， 2022，49（11）：252-259．（in Chinese）

[6] GOMASA R， TALAKOKULA V， KALYANA RAMA JYOSYULA S，et al ．A review on health monitoring of concrete structures using embedded piezoelectric sensor［J］．Construction and Building Materials，2023，405：133179．

[7] XU C J， LEI H G， WANG G Q ．Fatigue life and fatigue reliability assessment for long-span spatial structure based on long-term health monitoring data［J］．Structures，2023，47：586-594．

[8] NI P H， LI Q， HAN Q，et al ．Substructure approach for Bayesian probabilistic model updating using response reconstruction technique ［J］．Mechanical Systems and Signal Processing，2023，183：109624．

[9] YANG J S， FU Z Y， ZOU Y F，et al. A response reconstruction method based on empirical mode decomposition and modal synthesis method［J］．Mechanical Systems and Signal Processing，2023，184：109716．

[10] ZHOU K， WANG Y N， QIAO B J，et al ．Single-sensor-based dynamic response reconstruction of blades under base excitation ［J］．Mechanical Systems and Signal Processing，2023，192：110217．

[11] ZHANG X H， ZHU Z M， YUAN G K，et al ．Adaptive mode selection integrating Kalman filter for dynamic response reconstruction ［J］．Journal of Sound and Vibration，2021，515：116497．

[12] LIU Y R， WANG L， QIU Z P，et al ．A dynamic force reconstruction method based on modified Kalman filter using acceleration responses under multi-source uncertain samples［J］．Mechanical Systems and Signal Processing，2021，159：107761．

[13] HE J J， GUAN X F， LIU Y M ．Structural response reconstruction based on empirical mode decomposition in time domain ［J］．Mechanical Systems and Signal Processing，2012，28：348-366．

[14] LAI Z L， LEI Y， ZHU S Y，et al ．Moving-window extended Kalman filter for structural damage detection with unknown process and measurement noises［J］. Measurement，2016，88：428-440．

[15] WERNITZ S， CHATZI E， HOFMEISTER B，et al ．On noise covariance estimation for Kalman filter-based damage localization ［J］．Mechanical Systems and Signal Processing，2022，170：108808．

[16] TEYMOURI D， SEDEHI O， KATAFYGIOTIS L S，et al ．Input-state-parameter-noise identification and virtual sensing in dynamical systems：a Bayesian expectation-maximization （BEM） perspective［J］. Mechanical Systems and Signal Processing，2023，185：109758．

[17] ZHANG X H， WU Z B ．Dual-type structural response reconstruction based on moving-window Kalman filter with unknown measurement noise［J］. Journal of Aerospace Engineering， 2019， 32（4）： 04019029．

[18] TEYMOURI D， SEDEHI O， KATAFYGIOTIS L S，et al ．A Bayesian expectation-maximization （BEM） methodology for joint input-state estimation and virtual sensing of structures［J］. Mechanical Systems and Signal Processing，2022，169：108602．

[19] YUEN K V， KUOK S C ．Online updating and uncertainty quantification using nonstationary output-only measurement ［J］．Mechanical Systems and Signal Processing，2016，66：62-77．

[20] VETTORI S， DI LORENZO E， PEETERS B，et al ．An adaptive-noise augmented Kalman filter approach for input-state estimation in structural dynamics［J］. Mechanical Systems and Signal Processing， 2023， 184： 109654．

[21] ZHU Z M， ZHU S Y ．Asynchronous Kalman filtering for dynamic response reconstruction by fusing multi-type sensor data with arbitrary sampling frequencies［J］．Mechanical Systems and Signal Processing，2024，215：111395．

[22] ZHU Z M， LU J B， ZHU S Y ．Multi-rate Kalman filtering for structural dynamic response reconstruction by fusing multi-type sensor data with different sampling frequencies［J］. Engineering Structures，2023，293：116573．

[23] NING Z S， DENG Z W， LI J W，et al. Co-estimation of state of charge and state of health for 48 V battery system based on cubature Kalman filter and H-infinity［J］. Journal of Energy Storage，2022，56：106052．

[24] JI J L， HUANG Y J， PI M Q，et al ．Performance improvement of on-chip mid-infrared waveguide methane sensor using wavelet denoising and Savitzky-Golay filtering［J］. Infrared Physics & Technology，2022，127：104469．

[25] CHANG H， LI D Y， ZHANG R Y，et al ．Effect of discontinuous biomimetic leading-edge protuberances on the performance of vertical axis wind turbines［J］. Applied Energy，2024，364：123117．

[26] 郑国峰，陈柏先，陈文，等 ．基于应变响应的结构动态载荷识别方法［J］．振动·测试与诊断，2023，43（4）：779-786．

[27] ZHENG G F， CHEN B X， CHEN W，et al ．Study on dynamic excitation identification method of components based on strain measurements［J］. Journal of Vibration，Measurement & Diagnosis，2023，43（4）：779-786．（in Chinese）

[28] HOU Z L， YU L， LIANG Y，et al ．Integrating L1 and weighted L2 regularization for moving force identification from combined response measurements［J］．Measurement，2024，228：114337．

基金资助

国家自然科学基金资助项目(62161018)

National Natural Science