基于参数降阶模态的分层贝叶斯在线裂纹检测

2024-03-13魏保立冯雅珊罗坤付伟

机床与液压 2024年4期

魏保立，冯雅珊，罗坤，付伟

(1.郑州铁路职业技术学院机电工程学院，河南郑州 451460；2.河南理工大学能源科学与工程学院，河南焦作 454003)

0 前言

在基于状态的维修领域，面向模型的方法已成为解决退化和损伤识别问题的一种常用方法，不仅在检测和定位层面产生了重要作用，而且在更精细的量化阶段也能够发挥其优势[1]。与强烈依赖密集传感信息的纯数据驱动方案不同，系统模型的可用性存在一定程度的冗余，允许提取监测量和模型预测量之间的相互依赖关系[2]。如何设计一种精确的裂纹或损伤检测方法成为健康管理领域的关键。

目前，研究者已经在传统的裂纹或损伤检测的基础上提出了许多具有一定功能的替代方法。在各种方法中，递归贝叶斯估计器是适应实时性能要求的可靠手段，并且已经作为成熟的模型应用在航空航天工程、电信、金融和控制工程等领域[3-4]。在递归贝叶斯框架内，损伤识别本质上被归类为一个状态参数估计问题，这是一个非线性问题，即使对于实际动力学处于线性状态的系统也是如此。文献[5]引入贝叶斯分类方法，优化选择裂纹形核的信息状态变量预测因子，建立了一个简单的标量裂纹形核指示器。文献[6]针对当前疲劳裂纹扩展预测研究中较少考虑不确定因素而导致预测结果偏差大的问题，提出基于动态贝叶斯网络的疲劳裂纹扩展预测方法，以变幅载荷作用下的疲劳裂纹扩展为研究对象。利用统一疲劳寿命预测模型构建疲劳裂纹扩展的物理状态方程。虽然上述方法取得了一定效果，但是与这些滤波器相关的数值问题之一是样本贫化。

为此，研究者提出了一系列基于状态增强的参数估计方法。文献[7]提出了一种基于自然梯度法的变分贝叶斯卡尔曼滤波的变分下界最大化技术，得到了状态估计的变分超参数和相应的误差协方差的估计。文献[8]提出了一种新的综合预测模型贝叶斯优化随机森林-卡尔曼滤波，其中卡尔曼滤波和随机森林算法分别用于预测裂纹的趋势和周期。上述方法中的一个关键挑战是过程和测量噪声项的调整，因为在这些问题中，推理过程是基于初始参数值和互协方差项实现的，从本质上调节了观测状态的信息转化为未观测参数的修正过程。当参数本身表现出很强的相互关联时，问题变得更具挑战性。此外，过程噪声协方差交叉项的静态值可能不足以解决状态参数问题，其中系统的动态特性取决于参数估计，因此需要进行调整，以避免不稳定性。

为解决上述方法中存在的缺点，本文作者提出一种基于参数降阶模态的分层贝叶斯在线裂纹检测方法。在实时性能的约束下，系统动力学由一个参数化的降阶模型表示，该模型与所寻求的特征有关，随后与分层贝叶斯方法融合，解决参数空间采样后生成的多个模型配置解的输入状态估计问题。最后通过仿真验证所提方法的有效性。

1 问题描述

1.1 参数化降阶模型

文中的研究目的是仅输出振动测量值，在线检测结构或机械系统在运行条件下的裂纹，并进行输入状态估计。基本上可通过假设实际振动响应测量的可用性来实现，这些测量是从高保真(High Fidelity，HF)有限元模型和所考虑的系统中综合获得的。为了适应不同的裂纹结构，后者根据裂纹特征即位置、长度和方向进一步参数化。因此，连续时间系统动力学由以下参数相关的二阶微分方程表示：

(1)

其中：u∈Rn是位移向量，n表示有限元模型空间离散化产生的自由度数；M(θ)、C(θ)、K(θ)∈Rn×n，分别是质量矩阵、阻尼矩阵和刚度矩阵；θ∈Rm是参数向量，此处表示裂纹的几何特征。等式(1)的右侧由力矢量组成p(t)∈Rnp和相应的选择矩阵Sp∈Rn×np，它指定p(t)选择相应的自由度。应注意的是，上述方程通常是有限元离散化的结果，其应足够密集，以详细地捕捉系统的响应。因此，在在线设置中使用这些方程的积分可能会消耗大量资源。

为了使系统动力学预测的计算效率更高，从而允许在线解决输入、状态和参数估计问题，将式(1)的解投影到一个子空间上，该子空间的维数为r，比全订单模型小得多(r≪n)，由于系统动力学依赖于参数向量θ，减少的基础也是θ，这意味着投影步长可以表示为

u(t)=V(θ)q(t)

(2)

其中：V(θ)∈Rn×r包含约化空间的基向量；q(t)∈Rr是对应广义坐标的向量。将式(2)代入式(1)并结合uT得到降阶运动方程：

V(θ)TSpp(t)

(3)

通过本征正交分解(POD)可以获得全局约化基，从而求解多个参数样本的运动方程θj(j=1,2,…,ns)，ns表示样本的数量。然后，通过收集所有参数样本位移向量的迭代次数构建Snapshot矩阵，如式(4)所示：

S=[U(θ1)U(θ2) …U(θns)]

(4)

其中：U(θj)∈Rn×nt指定从等式(1)的解中获得的Snapshot矩阵θ=θj，nt表示解决方案步骤的数量，即Snapshot的数量。最后，可以通过最初提取的Snapshot矩阵的奇异值分解(Singular Value Decomposition，SVD)来执行缩减步骤，然后仅保留第r个显著性来截断相应的基向量，这些基向量用于形成全局约简基V∈Rn×r。对于具有大量时间步长的动态问题，可以增加一个额外的步骤来减少Snapshot矩阵的存储需求，该步骤涉及在Snapshot级别通过SVD进行缩减。

由于参数相关问题的全局POD方法通常会导致一个超大的基础V，因此，为了消除计算增益，文中采用了一种基于聚类的POD方案，其中包含局部基和网格变形，作为表示不同裂纹形态的一种手段。根据前述的过程，对参数空间的某些区域提取约简基V，这些区域通过k-均值聚类算法进行提取，并使用网格变形对未经训练的裂纹配置的解进行插值。

1.2 状态空间表示

为了将建模的系统动力学与实际振动测量值融合，需要将参数化的降阶运动方程即公式(3)转换为连续的时间-状态-空间表示，其表示为

(5)

y(t)=G(θ)x(t)+J(θ)p(t)

(6)

(7)

(8)

输出矩阵和馈通矩阵的结构G(θ)和J(θ)分别由观测到的响应量确定。对于结构系统，这些量通常是应变和加速度信号，因此得出以下表达式：

G(θ)=

(9)

(10)

其中：Sa∈Rna×n是一个布尔矩阵，用于选择加速度测量和；Sd，ε∈Rnd×n是与应变测量相关的选择矩阵。在处理有限元模型时，应变是在单元级别计算的，因此Sd，ε选择相应单元自由度的位移，然后使用矩阵将它转换为应变De(θ)∈Rnε×nd。后者具有块对角形式，每个块条目包含元素的变形矩阵。

2 递阶状态输入参数估计

2.1 状态输入估计

式(5)(6)在稀疏和噪声输出的基础上，仅对测量值进行离散化。为了适应离散的自然响应测量，将式(5)、(6)转移到离散时间域，如下所示：

xk+1=Ad(θ)xk+Bd(θ)pk+vk

(11)

yk=Gd(θ)xk+Jd(θ)pk+wk

(12)

同时进一步补充了零均值高斯噪声项vk∈R2r和wk∈Rny，其协方差矩阵分别为Qv∈R2r×2r和Qw∈Rny×ny。矩阵Ad(θ)、Bd(θ)、Gd(θ)和Jd(θ)是离散时间对应的A(θ)、B(θ)、G(θ)和J(θ)，其零阶保持表达式如下所示：Ad(θ)=eA(θ)dt、Bd(θ)=(A(θ)-I)·A-1(θ)B(θ)、Gd(θ)=G(θ)和Jd(θ)=J(θ)。系统参数θ在估计过程中发生变化时，可以使用泰勒级数展开，以牺牲精度为代价，更有效地逼近状态方程矩阵，从而得到Ad(θ)≈I+A(θ)dt+0.5A2(θ)dt2+O(dt3)及Bd(θ)≈B(θ)dt。

(13)

式(13)使用了过程模型的马尔可夫特性：p(xk|xk-1,Yk-1)=p(xk|Yk-1)。在计算先验知识后，可以使用贝叶斯定理进行更新：

(14)

其递推公式如下：

(15)

用由状态空间表示的似然函数p(yk|xk)的测量方程定义。所寻求的过滤概率密度函数的估计为p(xk|Yk)，在工程应用中，状态的统计矩通常是有意义的，可以计算如下：

(16)

(17)

E[·]表示期望算子和上述表达式的积分极限(-∞，+∞)。在线性系统方程组和加性高斯噪声的特殊情况下，由著名的卡尔曼滤波器获得式(11)(12)的解。虽然该系统由线性表达式(9)(10)描述，并且假设加性高斯噪声，即额外估计未知参数的问题θ和输入pk在状态xk下变为非线性，因此需要替代解决方案策略。序列蒙特卡罗(Sequential Monte Carlo，SMC)采样器构成了这一大类方法，也将在文中使用。

尽管状态增强是解决状态参数估计问题的一种直接技术，但该方法的实现强烈依赖于初始状态变量和追求的参数之间的耦合，以及初始状态和参数值的均值和协方差。由于参数本身不可观测，因此其估计取决于状态变量(可能部分可观测)。这意味着，参数动态由状态参数协方差的耦合项控制，它本质上调节了状态更新转化为未观测参数的修正。因此，状态和参数估计的精度与系统的实际误差统计密切相关，根据定义可知，这些误差统计是未知的，只能得到近似值。

为了避免使用未知参数增加状态向量，文中使用了一组滤波器，使它适应所寻求的系统参数。在这种情况下，未知的裂纹几何特征或通常要估计的系统参数为一组随机变量，状态和参数估计问题可分解为

(18)

(19)

图1 提出算法的结构框图

式(19)中所描述的分解本质上是Rao-Blackwellization操作，然而，由于某些原因，文中方案与标准Rao-Blackwellization粒子过滤器(RBPF)有所不同。也就是说，所提方法基于输入状态、参数更新和校正步骤的不同时间尺度。此外，参数动力学不受虚构模型的控制，在此类问题中通常假设随机游动。相反，通过进化策略(Evolutionary Strategies，ES)探索参数空间，即协方差矩阵自适应进化策略(Covariance Matrix Adaptation Evolutionary Strategies，CMA-ES)中使用的策略[9]。因此，收敛参数估计的关键要素由CMA-ES的步长保证，该步长强制参数核收缩。

(20)

联合输入状态估计可以使用AKF递归获得，AKF本质上是一个标准的卡尔曼滤波器(KF)，在增广状态空间模型上操作。

(21)

(22)

随后通过如下遗忘因子更新测量噪声的实际协方差矩阵：

(23)

其中：δQw表示用户定义的遗忘因子，用于调节预期更新的级别。δQw会减弱括号中的术语，得到非自适应方案，而如果δQw的值较小，则更新是突然进行的。文中所采用的协方差更新方案优于更为成熟的噪声识别方法，例如自协方差最小二乘(Autocovarianoe Least Squares，ALS)方法，由于可以通过递归实现该方法，因此能够在线检测。

通过试错法调整遗忘因子δQw，直到滤波器产生次优或接近最优的结果。文中采用了同样的启发式方法，这证明了该方法可以提供足够准确的结果，这可以归因于这样一个事实，即整个方法的关键因素不是每个输入状态估计器本身的最优性，而是滤波器组的公平权重，这也由次优自适应方案保证。然而，对于全自动校准δQw，更新策略可以与优化方案相结合，同时可以在LOO设置中进一步实施，以确保对创新序列的无偏估计，从而确保自适应的最佳性。

2.2 参数估计

(24)

当处于当前步骤k，可以假设yk-1满足p(yk-1|θj)=p(yk-1)=1，同样适用之前的所有观测yk-1，那么p(Yk-1|θj)=p(Yk-1)=1。因此p(θj|yk-1)=p(θj|Yk-1)=p(θj)，根据贝叶斯规则，在式(24)中进行替换后，得出以下递推公式：

(25)

参数向量的离散状态θ应达到参数空间要求的分辨率，在离线应用中，可以通过参数空间密集离散化来解决，原因为计算时间不是一个关键因素。然而，在需要实时性能且只允许有限数量粒子的在线应用中，情况并非如此。

(26)

2.3 更新策略

参数的更新不是由随机游走模型决定的，通常是在增广贝叶斯方案中假设的。相反，参数通过CMA-ES进行更新，超过简并阈值时调用CMA-ES。这意味着参数更新步骤是在与执行输入和状态校正的时间尺度不同的时间尺度下执行的。从这个层面上说，每个重采样事件都会创建新的参数样本，这些样本来自多元正态分布，因此：

(27)

(28)

(29)

算法1：用于参数和输入状态估计仅输出的分层滤波器。

(5)fork=1，2，…，Tdo

(6)forj=1，2，…，nθdo

(12)end for

(16)forj=1，2，…，nθdo

(20)end for

(21)end for

(30)

(31)

(32)

其中：

(33)

c1和cμ分别是更新协方差矩阵的排名第一和第nθ，p的学习率。算法1中记录了实现分层滤波器的详细步骤。

3 仿真案例

3.1 机身面板

通过2个仿真示例验证所提出的方法，利用稀疏传感信息识别裂纹特征。为了避免反演问题，每个示例的合成振动测量值的生成都基于高保真模型(HFM)，这2种应用中的高保真模型与用于反向问题的模型不同。即，HFM使用高度精细的网格构建，该网格独立于ROM使用的网格，而其材料特性的最小空间变化遵循高斯分布，其平均值等于弹性模量的标称值，其标准偏差等于平均值的1%。最后，向模拟的振动响应信号中随机添加高斯白噪声，每个信号的水平等于相应RMS的3%，以便创建一个尽可能接近实际的模拟装置。所有应用程序都在MATLAB中实现，硬件设施为3.80 GHz Intel Xeon E3-1275四核处理器的工作站，内存为32 GB。

下面的仿真研究基于裂纹面之间没有接触的假设。在没有拉伸载荷的情况下，这种假设可能不现实，拉伸载荷会打开裂纹，对于厚度为零的裂纹也不现实。然而，当前工作的主要目的不是准确地表示裂纹附近的物理现象，而是捕捉裂纹导致的局部刚度降低对系统整体响应的影响。此外，在这2种应用中，pROMs的训练空间决定参数搜索空间，同时选择参数样本的数量，以便在不牺牲精度的情况下最小化计算成本。然而，最小参数样本的数量与追求参数的数量有关。因此，参数的增加会相应地扩大参数样本量。

图2 面内应变(a)和面外加速度振动测量(b)

将系统动力学投影到一个由25个基向量构成的参数化降阶空间中。模型的参数化是根据裂纹特性进行的，即裂纹位置坐标xc∈[-0.15,0.15]、yc∈[-0.1,0.1]和长度lc∈[0.03，0.12]，其定义在x-y系统的投影如图2所示，参考裂纹与图3的相应网格。文中认为压力是面板表面上的均匀压力，其大小通过估计来确定。但这并不是文中方法的局限性，当使用降阶模型时，它可以通过高斯过程模型[13]或广义力项[14]轻松扩展，以适应非均匀压力的估计。

图3 参考网格和裂纹配置

表1 调整机身面板应用程序的参数值

使用文中方案估计的参数值如图4(a)所示，其中算法在模拟的前几秒钟内收敛到实际参数值。基于搜索空间的随机探索，算法运行5次得到的相应参数估计。图4(a)(b)分别为参数样本在2个不同时刻的演变情况，即t=0 s和t=7 s。可以看出：最初搜索分布在空间的很大一部分，并随着时间的推移而收敛到由黄色三角形标记的实际参数向量区域。t=7 s后，对参数空间重新采样，但是搜索分布已经收敛到非常接近目标值的位置，新粒子对估计的参数的改进不大。

图4 估计的参数值

使用pROM获得的状态估计值与用于生成合成振动数据的模型状态不可直接比较，因为后者是一个全阶模型，其状态与pROM的状态不一致。因此，根据测量和未测量位置的响应估计，可以证明预测状态的准确性。图5展示了点1、3处的预测位移响应，其位置如图2所示，从正向模拟中获得相应的噪声信号。因此，图2中定义的未测量点A、B、C和D处的预测位移响应如图6所示。可以看出：尽管信号的平均值存在一定偏差，但在4个位置都能很好地捕捉到实际的系统动力学。上述曲线图中，删除模拟时间的前5 s，因为响应主要由瞬态动力学控制，其量级比稳态响应大一个数量级，并且当与稳态响应一起绘制时，后者不可见。

图6 预测位移响应

图7 输入估计曲线

图8 图形化描述

3.2 翼箱面板

第二个示例是翼箱面板，如图9所示。面板由一块3 mm厚的板组成，该板由2个C形槽肋支撑，肋位于2个边缘x方向，并由2个角截面纵梁进一步加固，由虚线表示其位置。实际面板在每个肋的两端受到弹性支承的约束，弹性支承基本上代表了互连元件的刚度，但为了简单起见，文中认为面板完全固定在肋的两端。这种固定的设置使系统比实际情况稍硬，其动态特性对局部变化不太敏感，其识别问题比文献[15]中研究的问题更具挑战性。假定面板的材料特性为线弹性，且E=73 GPa、ν=0.3、密度ρ=2 700 kg/m3，在第一和第四振型中进一步引入了2%的比例阻尼。

图9 机翼箱面板几何图形和裂纹参数定义的示意

图9还显示了面板的尺寸(单位：m)和传感器位置，以及第一桁条上的裂缝，该裂缝垂直于构件的纵轴。用于识别问题的pROM由15个基向量组成，从y方向开始，裂纹相对位置xc进行参数化，范围在距纵梁中点0.2 ～2 m，或笛卡尔坐标系下0.15～0.55 m之间，长度为lc∈[0.002 5，0.022 5] m。图10描述了用于ROM构造的面板参考离散化以及参考裂纹。在操作条件下测试面板的振动响应，假设动力学由均匀压力驱动，该压力的大小为高斯白噪声过程，并沿振动方向施加在板上z轴，与之前的应用类似，文中根据量级识别该输入p(t)。

图10 参考离散化(a)以及参考裂纹(b)

表2 调整机翼箱面板应用程序的参数

在计算性能方面，整个算法平均需要472 s，模拟周期为30 s。算法中最耗时的部分是新系统矩阵的生成，大约需要90%的执行时间，即419 s，而剩下的53 s用于依次执行所有粒子的预测、滤波和噪声适应步骤。后者几乎满足实时处理，可以通过并行处理，从而实现数字孪生。在这种情况下，可以持续检测物理系统，并根据测得的响应更新响应状态。另外，降低阶数建模部分的实现在计算时间方面没有进行优化，这也可以在进一步并行化时显著提高速度。

从算法的5次不同运行中获得的参数估计如图11所示，根据模拟的不同时刻的估计历史时间和参数样本。从图11(a)可以观察到：在前5 s内，参数估计非常接近实际值。橙色轨道的参数估计如图11所示，其相应颜色的粒子权重在整个仿真时间内的更新如图12所示。在前2.5 s内，粒子快速探索搜索空间，到达目标附近。在密集重采样事件中，当简并阈值达到时，整个权重集中到单个粒子，从而产生高度接近1的峰值。然后，在目标附近继续搜索，直到t=8.5 s附近，粒子被略微拉开，经过一段密集的重采样后，在t=11.5 s附近返回目标点。此后，算法保持接近目标，只进行有限的调整，对参数估计没有任何显著的改进。需要强调的是，参数值在整个模拟过程中不断更新，因此图12中所示的权值曲线不是恒定的参数值。

图11 参数估计结果(a)和粒子权重(b)

图12 权值曲线

运行系统的预测响应时间如图13所示。将观测到应变和加速度的点1处的预测位移和点2处的加速度与图13(a)中的测量结果进行比较。这些信号再次出现在不包含瞬态动力学的时间窗口中，从2.5 s开始。值得注意的是，在测量点处较好地捕获振动信号，这在很大程度上是由于噪声调谐，下面将进一步讨论。另外，在图13(b)中观察到略有不同的图片，其中以位移和加速度表示未测量点A和B的响应。点A的动力学在位移方面得到了充分的追踪，在响应的峰值处出现了一定的差异。这可以部分归因于预测裂纹参数中的误差以及系统模型的降阶，该降阶是基于实时性能选择的，以牺牲准确性为代价。在实际监控场景中，可以相应地调整系统顺序，以平衡计算需求，提高精度。相反，点B处预测加速度的误差主要与估计输入的误差有关，如图14所示。虽然位移和速度估计仅取决于未观测状态估计的质量，但加速度预测强烈依赖于输入估计。

图13 预测响应时间

图14 输入估计

最后，为了证明测量噪声项的协方差矩阵自适应的效果，并进一步对遗忘因子δQw的选择值进行推理，对δQw的不同值进行实验。通常，噪声协方差矩阵的调整旨在实现滤波器的最佳性能，可以通过使用不同的最佳性能测试及噪声的白度进行验证。在此例中，使用预测残差的频域表示在点1处的定性检查最优性e1(t)，图15(a)中绘制了遗忘因子的3个不同值。图15(b)中显示了相应的时域残差。为了滤除密集重采样事件引入的偏差，重新初始化滤波器并导致残差突然增加，图15所示的结果是基于12.5 s后的样本，由充分收敛的模型生成这些样本，如图11(a)所示。

图15 遗忘因子测量值和预测值之间的差异

遗忘因子通常可以进行更贪婪地探索，文献[15]中提出了遗忘因子，但是，这不在文中的研究范围。文中旨在证明协方差矩阵自适应性对算法性能的影响，并为实际应用中遗忘因子的调整提供进一步的指导。图15中显示的输入、状态和参数结果是令δQw=4，如图15(a)所示，与在试验中获得的光谱相比，产生了白色光谱δQw=3和5。另外，δQw意味着协方差矩阵的适应非常缓慢，测量噪声仍然接近初始值。另一方面，接近1的值会根据计算的残差突然改变协方差矩阵。在此例中，由于预测和测量之间存在较大差异，小于2.5的δQw在某些过滤器中产生不稳定性，即远离目标。相反，当遗忘因子值较大时，这些问题并没有表现出来，但是，其预测精度较低，反映在高度相关的残差上，与δQw=5相似。