CN106446503A - 遗忘自协方差矩阵递推主元的时变工作模态识别方法 - Google Patents

Abstract

本发明公开了一种遗忘自协方差矩阵递推主元的时变工作模态识别方法，包括：获取线性时变结构在环境激励下多个振动响应传感器从初始时刻0到时刻k的非平稳信号数据矩阵归一化后求其自协方差矩阵并进行形式的特征向量分解，储存V^(k)。获取下一时刻的时域振动响应信号数据，对新的自协方差矩阵进行递归推导时，加入遗忘因子，分配给新旧数据不同的权重，得到V^(k+1)；循环上述推导步骤，能够得到任意时刻的V^(k)，V^(k)对应k时刻该结构的瞬态工作模态振型矩阵，利用单自由度识别技术对矩阵V^(k)TX_k进行处理，得到k时刻该结构的瞬时工作固有模态频率。该方法能够有效监测线性工程结构工作模态参数的时变结构特性，可被用于设备故障诊断、健康监测以及系统结构分析与优化。

Description

遗忘自协方差矩阵递推主元的时变工作模态识别方法

技术领域

本发明涉及工作模态参数识别领域，特别是涉及一种遗忘自协方差矩阵递推主元的时变工作模态识别方法。

背景技术

随着科学技术的发展与进步，航空航天、建筑、桥梁、海洋、机械等领域的工程结构逐渐向大型化、复杂化、智能化方向发展，结构所承受激励载荷难以测量，若要建立其动力学模型，传统的基于测量输入输出求取系统的模态参数方法难以适用。仅利用环境激励下的振动响应数据就能识别结构模态参数的工作模态参数识别方法越来越受欢迎。并且与此同时，很多工程结构在实际工作运行中，由于各种内外激励作用，其结构参数(质量、刚度、阻尼等)会随时间发生变化，常伴随有时变特性，变成线性时变结构。例如，火箭或导弹在飞行发射过程中，燃料消耗导致质量随时间减小；高速列车会引起桥梁的剧烈振动，由于列车自身带有质量，从而使车与桥耦合系统的动力特性随列车位置的移动而不断变化。除此之外，还有诸如高超声速飞行器结构气动加热、热固耦合振动时变特性等。当前现代工业生产对设备的安全性、稳定性要求逐渐加强，迫切希望对设备实施监测，提前预报故障，否则一旦发生故障，将造成严重的损失。因此，大量的线性时变结构问题亟待解决，研究线性时变结构的参数辨识问题也显得越来越重要，研究设备的状态监测和故障诊断也为现代企业的重要组成模块，具有着重要的意义。

发明内容

本发明的目的在于实现对线性时变结构的时变瞬变工作模态参数的识别，尤其是时变瞬变工作模态振型的识别，为线性时变结构的研究提供准确的参考，提供了一种基于带遗忘因子加权和自协方差矩阵递推主元分析算法的线性时变结构的时变瞬态工作模态参数识别方法及系统。

为实现上述目的，本发明提供了如下方案：

一种遗忘自协方差矩阵递推主元的时变工作模态识别方法，包括：

S101：获取线性时变结构在环境激励下多个振动传感器从初始时刻0到时刻k的非平稳时域振动响应信号数据矩阵

其中，表示维度为m×k的矩阵，m表示在所述线性时变结构上布置的振动传感器检测点个数，k表示时域采样点个数；1≤j≤m；1≤i≤k；

S102：将所述数据矩阵归一化，得到归一化矩阵X_k，求的均值向量与的标准差矩阵Σ_k＝diag{δ_k(1),δ_k(2),…,δ_k(j),…,δ_k(m)}；

S103：求取所述归一化矩阵X_k的自协方差矩阵 所述自协方差矩阵为实对称方阵；

S104：将所述自协方差矩阵分解为的特征值特征向量形式，其中， 为的特征值由大到小排列形成的对角矩阵， 依次为特征值的特征向量，V^(k)T为V^(k)的转置，V^(k)TV^(k)＝I_m×m；

S105：存储矩阵V^(k)，对应k时刻该结构的瞬态正则化工作模态振型矩阵根据矩阵V^(k)与X_k，主成分V^(k)TX_k对应k时刻各阶模态响应矩阵利用单自由度模态识别技术对矩阵进行处理，得到时刻k的瞬时工作模态固有频率；

S106：获取下一时刻(k+1时刻)的振动数据

S107：获取新的时域振动响应数据矩阵

S108：将所述数据矩阵归一化，得到归一化矩阵X_k+1，引入新样本数据后，原始数据矩阵变为 的均值向量的标准差矩阵Σ_k+1＝diag{δ_k+1(1),δ_k+1(2),…,δ_k+1(j),…,δ_k+1(m)}和X_k+1的协方差矩阵按如下公式进行更新：

其中，Σ_k+1＝diag{δ_k+1(1),δ_k+1(2),…,δ_k+1(j),…,δ_k+1(m)}为的标准差矩阵；δ_k+1(j)为的标准差矩阵中的第j个元素，1≤j≤m；为的均值向量；为的均值向量；为均值向量的变化量；为k+1时刻的振动数据的归一化后的结果；

S109：引入遗忘因子γ，将新旧数据分以不同权重，则的均值向量 的标准差Σ_k+1＝diag{δ_k+1(1),δ_k+1(2),…,δ_k+1(j),…,δ_k+1(m)}的更新公式如下所示：

根据自协方差矩阵的递推公式得到所述归一化矩阵X_k+1的自协方差矩阵 为X_k+1的转置；

S110：将分解为的特征值特征向量形式，其中， 为的特征值由大到小排列形成的对角矩阵， 依次为特征值的特征向量，V^(k+1)T为V^(k+1)的转置，V^(k+1)TV^(k+1)＝I_m×m；

S111：存储矩阵V^(k+1)，对应k+1时刻该结构的瞬态正则化工作模态振型矩阵

S112：根据矩阵V^(k+1)与X_k+1，主成分V^(k+1)TX_k+1对应k+1时刻各阶模态响应矩阵利用单自由度模态识别技术对进行处理，得到时刻k+1的瞬时工作模态固有频率；

S113：当k小于等于时域采样点总个数n时(k≤n)，将k+1的值代替k的值(k←k+1)，跳转至S106，否则，线性时变结构瞬态工作模态参数识别结束。

一种基于带遗忘因子加权和自协方差矩阵递推主元分析算法的线性时变结构的时变瞬态工作模态参数识别系统，其特征在于，所述系统包括：

201非平稳时域振动响应数据获取单元：用于获取线性时变结构在环境激励下多个振动传感器从初始时刻0到时刻k的非平稳时域振动响应信号数据矩阵

其中，表示维度为m×k的矩阵，m表示在所述线性时变结构上布置的振动传感器检测点个数，k表示时域采样点个数；

202矩阵归一化单元：用于将所述数据矩阵归一化，得到归一化矩阵X_k，求的均值向量与的标准差矩阵Σ_k＝diag{δ_k(1),δ_k(2),…,δ_k(j),…,δ_k(m)}；

203自协方差矩阵求取单元：用于求取所述归一化矩阵X_k的自协方差矩阵 所述自协方差矩阵为实对称方阵；

204自协方差矩阵分解单元：用于将所述自协方差矩阵分解为的特征值特征向量形式，其中， 为的特征值由大到小排列形成的对角矩阵， 依次为特征值的特征向量，V^(k)T为V^(k)的转置，V^(k)TV^(k)＝I_m×m；

205矩阵存储单元：用于存储矩阵V^(k)，对应k时刻该结构的瞬态工作模态振型矩阵；根据矩阵V^(k)与X_k，利用单自由度模态识别技术对矩阵V^(k)TX_k进行处理，得到时刻k的瞬时工作模态固有频率；

206下一时刻数据获取单元：用于获取下一时刻(k+1时刻)的振动响应数据

207新时域振动响应数据获取单元：用于获取新的时域振动响应数据矩阵

208新矩阵递归归一化单元：用于将所述数据矩阵归一化，得到归一化矩阵X_k+1，引入新样本数据后，原始数据矩阵变为 的均值向量的标准差矩阵Σ_k+1＝diag{δ_k+1(1),δ_k+1(2),…,δ_k+1(j),…,δ_k+1(m)}和X_k+1的协方差矩阵按如下公式进行更新：

其中，Σ_k+1＝diag{δ_k+1(1),δ_k+1(2),…,δ_k+1(j),…,δ_k+1(m)}为的标准差矩阵；δ_k+1(j)为的标准差矩阵中的第j个元素，1≤j≤m；为的均值向量；为的均值向量；为均值向量的变化量；为k+1时刻的振动数据的归一化后的结果；

209新自协方差矩阵递归求取单元：引入遗忘因子γ，将新旧数据分以不同权重，则的均值向量的标准差Σ_k+1＝diag{δ_k+1(1),δ_k+1(2),…,δ_k+1(j),…,δ_k+1(m)}的更新公式如下所示：

根据自协方差矩阵的递推公式递推得到所述归一化矩阵X_k+1的自协方差矩阵 为X_k+1的转置；

210：新自协方差矩阵分解单元，用于将分解为的特征值特征向量形式，其中， 为的特征值由大到小排列形成的对角矩阵， 依次为特征值的特征向量，V^(k+1)T为V^(k+1)的转置，V^(k+1)TV^(k+1)＝I_m×m；

211：新矩阵存储单元，用于存储矩阵V^(k+1)，对应k+1时刻该结构的瞬态工作模态振型矩阵；

212：瞬时工作模态参数计算单元，根据矩阵V^(k+1)与X_k+1，利用单自由度模态识别技术对矩阵V^(k+1)TX_k+1进行处理，得到时刻k+1的瞬时工作模态固有频率。

213：跳转单元，用于当k小于等于时域采样点总个数n时(k≤n)，将k+1的值代替k的值(k←k+1)，跳转至下一时刻数据获取单元，否则，线性时变结构瞬时工作模态参数识别结束。

本发明所述的基于带老化遗忘因子和自协方差矩阵递推主元分析算法的线性时变结构工作模态参数在线实时识别方法，能对带有时变特性的线性结构进行在线实时的参数识别，识别出线性时变结构的工作模态参数：瞬时模态振型和瞬时模态频率，能够有效监测结构工作模态参数的动态变化特性，可被用于设备故障诊断、健康监测以及系统结构分析与优化。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例线性时变结构工作模态参数识别方法的流程示意图；

图2为本发明实施例线性时变结构工作模态参数识别系统的示意图；

图3是基于带遗忘因子与自协方差矩阵递推主元分析算法的递推数据模型；

图4是模拟环境激励的线性三自由度弹簧振子振动系统；

图5是白噪声激励信号；

图6是第一个振子单元的时域位移响应信号；

图7是第二个振子单元的时域位移响应信号；

图8是第三个振子单元的时域位移响应信号；

图9是三自由度时变结构固有频率变化曲线；

图10是带遗忘因子的递推主元分析算法识别的三自由度时变结构的频率变换曲线；

图11是三自由度时变结构固有频率变化曲线与带遗忘因子的递推主元分析算法识别的三自由度时变结构的频率变换曲线的比较曲线；

图12是50.025s时刻的三自由度时变系统第一阶固有振型与带遗忘因子递推主元分析算法识别的第一阶瞬时振型的比较

图13是50.025s时刻的三自由度时变系统第二阶固有振型与带遗忘因子递推主元分析算法识别的第二阶瞬时振型的比较

图14是50.025s时刻的三自由度时变系统第三阶固有振型与带遗忘因子递推主元分析算法识别的第三阶瞬时振型的比较

图15是50.025s时刻的三自由度时变系统带遗忘因子递推主元分析算法识别的各主成份所占贡献量

图16是500s时刻的三自由度时变系统第一阶固有振型与带遗忘因子递推主元分析算法识别的第一阶瞬时振型的比较

图17是500s时刻的三自由度时变系统第二阶固有振型与带遗忘因子递推主元分析算法识别的第二阶瞬时振型的比较

图18是500s时刻的三自由度时变系统第三阶固有振型与带遗忘因子递推主元分析算法识别的第三阶瞬时振型的比较

图19是500s时刻的三自由度时变系统带遗忘因子递推主元分析算法识别的各主成份所占贡献量

图20是950s时刻的三自由度时变系统第一阶固有振型与带遗忘因子递推主元分析算法识别的第一阶瞬时振型的比较

图21是950s时刻的三自由度时变系统第二阶固有振型与带遗忘因子递推主元分析算法识别的第二阶瞬时振型的比较

图22是950s时刻的三自由度时变系统第三阶固有振型与带遗忘因子递推主元分析算法识别的第三阶瞬时振型的比较

图23是950s时刻的三自由度时变系统带遗忘因子递推主元分析算法识别的各主成份所占贡献量

图24是1400s时刻的三自由度时变系统第一阶固有振型与带遗忘因子递推主元分析算法识别的第一阶瞬时振型的比较

图25是1400s时刻的三自由度时变系统第二阶固有振型与带遗忘因子递推主元分析算法识别的第二阶瞬时振型的比较

图26是1400s时刻的三自由度时变系统第三阶固有振型与带遗忘因子递推主元分析算法识别的第三阶瞬时振型的比较

图27是1400s时刻的三自由度时变系统带遗忘因子递推主元分析算法识别的各主成份所占贡献量

图28是1987.225s时刻的三自由度时变系统第一阶固有振型与带遗忘因子递推主元分析算法识别的第一阶瞬时振型的比较

图29是1987.225s时刻的三自由度时变系统第二阶固有振型与带遗忘因子递推主元分析算法识别的第二阶瞬时振型的比较

图30是1987.225s时刻的三自由度时变系统第三阶固有振型与带遗忘因子递推主元分析算法识别的第三阶瞬时振型的比较

图31是1987.225s时刻的三自由度时变系统带遗忘因子递推主元分析算法识别的各主成份所占贡献量

图32是三自由度时变系统的第一阶振型MAC值比较变化曲线。

图33是三自由度时变系统的第二阶振型MAC值比较变化曲线。

图34是三自由度时变系统的第二阶振型MAC值比较变化曲线。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明的目的是实现对线性时变结构的时变瞬变参数，尤其是时变瞬变模态振型的测量，为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

本发明采用带遗忘因子老化与自协方差矩阵递推主元分析法相结合，形成能够识别线性时变结构时变的瞬态工作模态参数的方法，该方法利用统计平均的方式识别出系统各时刻的瞬时模态参数(瞬时模态频率和瞬时模态振型)，然后曲线拟合各时刻求得的工作模态参数，以达到在线实时监测结构的动力学特性变化的目的。

图1为本发明实施例线性时变结构工作模态参数识别方法的流程示意图，如图1所示，本发明提供的时变结构工作模态参数识别方法具体包括：

S101：获取线性时变结构在环境激励下多个振动传感器从初始时刻0到时刻k的非平稳时域振动响应信号数据矩阵

其中，表示维度为m×k的矩阵，m表示在所述线性时变结构上布置的振动传感器检测点个数，k表示时域采样点个数；1≤j≤m；1≤i≤k；

S102：将所述数据矩阵归一化，得到归一化矩阵X_k，求的均值向量与的标准差矩阵Σ_k＝diag{δ_k(1),δ_k(2),…,δ_k(j),…,δ_k(m)}；

S103：求取所述归一化矩阵X_k的自协方差矩阵 所述自协方差矩阵为实对称方阵；

S104：将所述自协方差矩阵分解为的特征值特征向量形式，其中， 为的特征值由大到小排列形成的对角矩阵， 依次为特征值的特征向量，V^(k)T为V^(k)的转置，V^(k)TV^(k)＝I_m×m；

S105：存储矩阵V^(k)，对应k时刻该结构的瞬态正则化工作模态振型矩阵根据矩阵V^(k)与X_k，主成分V^(k)TX_k对应k时刻各阶模态响应矩阵利用单自由度模态识别技术对矩阵进行处理，得到时刻k的瞬时工作模态固有频率；

S106：获取下一时刻(k+1时刻)的振动数据

S107：获取新的时域振动响应数据矩阵

S108：将所述数据矩阵归一化，得到归一化矩阵X_k+1，引入新样本数据后，原始数据矩阵变为 的均值向量的标准差矩阵Σ_k+1＝diag{δ_k+1(1),δ_k+1(2),…,δ_k+1(j),…,δ_k+1(m)}和X_k+1的协方差矩阵按如下公式进行更新：

其中，Σ_k+1＝diag{δ_k+1(1),δ_k+1(2),…,δ_k+1(j),…,δ_k+1(m)}为的标准差矩阵；δ_k+1(j)为的标准差矩阵中的第j个元素，1≤j≤m；为的均值向量；为的均值向量；为均值向量的变化量；为k+1时刻的振动数据的归一化后的结果；

S109：引入遗忘因子γ，将新旧数据分以不同权重，则的均值向量 的标准差Σ_k+1＝diag{δ_k+1(1),δ_k+1(2),…,δ_k+1(j),…,δ_k+1(m)}的更新公式如下所示：

根据自协方差矩阵的递推公式得到所述归一化矩阵X_k+1的自协方差矩阵 为X_k+1的转置；

S110：将分解为的特征值特征向量形式，其中， 为的特征值由大到小排列形成的对角矩阵， 依次为特征值的特征向量，V^(k+1)T为V^(k+1)的转置，V^(k+1)TV^(k+1)＝I_m×m；

S111：存储矩阵V^(k+1)，对应k+1时刻该结构的瞬态正则化工作模态振型矩阵

S112：根据矩阵V^(k+1)与X_k+1，主成分V^(k+1)TX_k+1对应k+1时刻各阶模态响应矩阵利用单自由度模态识别技术对进行处理，得到时刻k+1的瞬时工作模态固有频率；

S113：当k小于等于时域采样点总个数n时(k≤n)，将k+1的值代替k的值(k←k+1)，跳转至S106，否则，线性时变结构瞬态工作模态参数识别结束。

在步骤S112中V^(k+1)即为所述线性时变结构在时刻k+1的模态振型矩阵，利用单自由度模态识别技术对矩阵V^(k+1)TX_k+1进行处理，得到时刻k+1的瞬时模态固有频率。

上述步骤的原理为：基于主成分分析原理，X_k可被唯一分解为X_k＝V^(k)(V^(k)TX_k)，X_k为多个振动传感器从初始时刻0到时刻k的非平稳时域振动响应信号数据矩阵，建立PCA初始化模型，其中，是主元分析中的变换阵，V^(k)TX_k是从初始时刻0到时刻k收集的样本长度为k的非平稳时域振动响应信号X_k的主成分，各主成分彼此之间不相关；

对于从初始时刻0到时刻k收集的长度为k的非平稳时域振动响应信号X_k在模态坐标下表示为X_k≈ΦQ_k；其中，正则化模态振型矩阵满足Φ^TΦ＝I_m×m，模态坐标响应矩阵

其中的各阶模态响应相互独立；

可以认为是线性时变结构在k时刻的瞬时模态振型的一个近似估计。

是线性时变结构在k时段内的模态坐标响应，利用单自由度模态识别技术(通过傅立叶变换，最高峰值处对应模态频率)，可以识别在k时刻的瞬时模态固有频率。

因为相互独立必定不相关，所以基于主成分分析原理，正则化模态振型矩阵对应主元分析中的线性混叠矩阵各阶模态响应矩阵为主成分分析中的主成分V^(k)TX_k。

本发明所述的基于带遗忘因子的自协方差矩阵递推主元分析算法的线性时变结构工作模态参数识别方法，能对带有时变特性的线性结构进行在线实时的参数识别，识别出线性时变结构的时变瞬态工作模态参数：时变的瞬时模态振型和瞬时模态频率，能够有效监测结构工作模态参数的动态变化特性，可被用于设备故障诊断、健康监测以及系统结构分析与优化。

本发明还提供了一种基于带遗忘因子加权和自协方差矩阵递推主元分析算法的线性时变结构的时变瞬态工作模态参数识别系统，图2为本发明实施例基于带遗忘因子加权和自协方差矩阵递推主元分析算法的线性时变结构的时变瞬态工作模态参数识别系统的示意图，如图2所示，该系统包括：时域振动响应数据获取单元201，用于获取线性时变结构在环境激励下多个振动传感器从初始时刻0到时刻k的非平稳时域振动响应信号数据矩阵

其中，表示维度为m×k的矩阵，m表示在所述线性时变结构上布置的振动传感器检测点个数，k表示时域采样点个数；矩阵归一化单元202，用于将所述数据矩阵归一化，得到归一化矩阵X_k，求的均值向量与的标准差矩阵Σ_k＝diag{δ_k(1),δ_k(2),…,δ_k(j),…,δ_k(m)}；自协方差矩阵求取单元203，用于求取所述归一化矩阵X_k的自协方差矩阵 所述自协方差矩阵为实对称方阵；自协方差矩阵分解单元204，用于将所述自协方差矩阵分解为的特征值特征向量形式，其中， 为的特征值由大到小排列形成的对角矩阵， 依次为特征值的特征向量，V^(k)T为V^(k)的转置，V^(k)TV^(k)＝I_m×m；矩阵存储单元205，用于存储矩阵V^(k)，对应k时刻该结构的瞬态工作模态振型矩阵；根据矩阵V^(k)与X_k，利用单自由度模态识别技术对矩阵V^(k)TX_k进行处理，得到时刻k的瞬时工作模态固有频率；下一时刻数据获取单元206，用于获取下一时刻(k+1时刻)的振动响应数据新时域振动响应数据获取单元207，用于获取新的时域振动响应数据矩阵 新矩阵递归归一化单元208，用于将所述数据矩阵归一化，得到归一化矩阵X_k+1，引入新样本数据后，原始数据矩阵变为 的均值向量 的标准差矩阵Σ_k+1＝diag{δ_k+1(1),δ_k+1(2),…,δ_k+1(j),…,δ_k+1(m)},

和X_k+1的协方差矩阵

其中，Σ_k+1＝diag{δ_k+1(1),δ_k+1(2),…,δ_k+1(j),…,δ_k+1(m)}为的标准差矩阵；δ_k+1(j)为的标准差矩阵中的第j个元素，1≤j≤m；为的均值向量；为的均值向量；为均值向量的变化量；为k+1时刻的振动数据的归一化后的结果；新自协方差矩阵递归求取单元209，引入遗忘因子γ，将新旧数据分以不同权重，则的均值向量 的标准差Σ_k+1＝diag{δ_k+1(1),δ_k+1(2),…,δ_k+1(j),…,δ_k+1(m)},

根据自协方差矩阵的递推公式递推得到所述归一化矩阵X_k+1的自协方差矩阵 为X_k+1的转置；

新自协方差矩阵分解单元210：用于将分解为的特征值特征向量形式，其中， 为的特征值由大到小排列形成的对角矩阵， 依次为特征值的特征向量，V^(k+1)T为V^(k+1)的转置，V^(k+1)TV^(k+1)＝I_m×m；；新矩阵存储单元211，用于存储矩阵V^(k+1)，对应k+1时刻该结构的瞬态工作模态振型矩阵；瞬时工作模态参数计算单元212，根据矩阵V^(k+1)与X_k+1，利用单自由度模态识别技术对矩阵V^{(k +1)T}X_k+1进行处理，得到时刻k+1的瞬时工作模态固有频率；跳转单元213，用于当k小于等于时域采样点总个数n时(k≤n)，将k+1的值代替k的值(k←k+1)，跳转至下一时刻数据获取单元，否则，线性时变结构瞬时工作模态参数识别结束。

本发明所述的基于带遗忘因子自协方差矩阵递推主元分析算法的时变结构工作模态参数识别系统，能对带有时变特性的线性结构进行在线实时的参数识别，识别出线性时变结构的工作模态参数：时变的瞬时工作模态振型和时变的瞬时工作模态频率，能够有效监测结构工作模态参数的动态变化特性，可被用于设备故障诊断、健康监测以及系统结构分析与优化。

本实施例中，所述的基于带遗忘因子和自协方差矩阵递推主元分析的线性时变结构工作模态参数识别装置采用三自由度弹簧振子模拟时变结构，如图4所示，其中，m₂＝1kg，m₃＝1kg；k₁＝1000N/m，k₂＝1000N/m，k₃＝1000N/m；c₁＝0.01N.s/m，c₂＝0.01N.s/m，c₃＝0.01N.s/m。初始位移及速度为0。物块1受到均值为0，方差为1的高斯白噪声激励F₁。基于Matlab/Simulink仿真，采样间隔为0.025s，采样频率为40Hz，仿真时间为2000s。

图5是白噪声激励信号；

图6是第一个振子单元的时域位移响应信号；

图7是第二个振子单元的时域位移响应信号；

图8是第三个振子单元的时域位移响应信号；

图9是三自由度时变结构固有频率变化曲线；

图10是带遗忘因子的递推主元分析算法识别的三自由度时变结构的频率变换曲线；

图11是三自由度时变结构固有频率变化曲线与带遗忘因子的递推主元分析算法识别的三自由度时变结构的频率变换曲线的比较曲线；

如图9所示，为通过理论计算的三阶时变的固有频率；如图10所示，为通过带遗忘因子递推主元分析算法识别的时变频率变化曲线；如图11所示，为时变频率理论值与该算法识别值比较曲线。

图12是50.025s时刻的三自由度时变系统第一阶固有振型与带遗忘因子递推主元分析算法识别的第一阶瞬时振型的比较

图13是50.025s时刻的三自由度时变系统第二阶固有振型与带遗忘因子递推主元分析算法识别的第二阶瞬时振型的比较

图14是50.025s时刻的三自由度时变系统第三阶固有振型与带遗忘因子递推主元分析算法识别的第三阶瞬时振型的比较

图15是50.025s时刻的三自由度时变系统带遗忘因子递推主元分析算法识别的各主成份所占贡献量

图16是500s时刻的三自由度时变系统第一阶固有振型与带遗忘因子递推主元分析算法识别的第一阶瞬时振型的比较

图17是500s时刻的三自由度时变系统第二阶固有振型与带遗忘因子递推主元分析算法识别的第二阶瞬时振型的比较

图18是500s时刻的三自由度时变系统第三阶固有振型与带遗忘因子递推主元分析算法识别的第三阶瞬时振型的比较

图19是500s时刻的三自由度时变系统带遗忘因子递推主元分析算法识别的各主成份所占贡献量

图20是950s时刻的三自由度时变系统第一阶固有振型与带遗忘因子递推主元分析算法识别的第一阶瞬时振型的比较

图21是950s时刻的三自由度时变系统第二阶固有振型与带遗忘因子递推主元分析算法识别的第二阶瞬时振型的比较

图22是950s时刻的三自由度时变系统第三阶固有振型与带遗忘因子递推主元分析算法识别的第三阶瞬时振型的比较

图23是950s时刻的三自由度时变系统带遗忘因子递推主元分析算法识别的各主成份所占贡献量

图24是1400s时刻的三自由度时变系统第一阶固有振型与带遗忘因子递推主元分析算法识别的第一阶瞬时振型的比较

图25是1400s时刻的三自由度时变系统第二阶固有振型与带遗忘因子递推主元分析算法识别的第二阶瞬时振型的比较

图26是1400s时刻的三自由度时变系统第三阶固有振型与带遗忘因子递推主元分析算法识别的第三阶瞬时振型的比较

图27是1400s时刻的三自由度时变系统带遗忘因子递推主元分析算法识别的各主成份所占贡献量

图28是1987.225s时刻的三自由度时变系统第一阶固有振型与带遗忘因子递推主元分析算法识别的第一阶瞬时振型的比较

图29是1987.225s时刻的三自由度时变系统第二阶固有振型与带遗忘因子递推主元分析算法识别的第二阶瞬时振型的比较

图30是1987.225s时刻的三自由度时变系统第三阶固有振型与带遗忘因子递推主元分析算法识别的第三阶瞬时振型的比较

图31是1987.225s时刻的三自由度时变系统带遗忘因子递推主元分析算法识别的各主成份所占贡献量

由于线性时变结构中模态振型时刻变化，难以列举出所有振型，基于此，在2000s仿真时间内，随机选取50.025s、500s、950s、1400s、1987.225s(避免随机振动的影响，选取50s时刻之后的数据进行计算)，如图12ˉ图15所示为50.025s时刻的各阶振型通过理论计算与遗忘因子递推主元分析算法识别的比较以及各主成份累积贡献率；，如图16ˉ图19所示为500s时刻的各阶振型通过理论计算与遗忘因子递推主元分析算法识别的比较以及各主成份累积贡献率；，如图20ˉ图23所示为950s各时刻的各阶振型通过理论计算与遗忘因子递推主元分析算法识别的比较以及各主成份累积贡献率；，如图24ˉ图27所示为1400s时刻的各阶振型通过理论计算与遗忘因子递推主元分析算法识别的比较以及各主成份累积贡献率；，如图28ˉ图31所示为1987.225s时刻的各阶振型通过理论计算与遗忘因子递推主元分析算法识别的比较以及各主成份累积贡献率；从各图中可以看出，遗忘因子递推主元分析算法能基本很好的识别各振型，且在各时刻第一主成份占主要成分；如图32、33、34所示，为三自由度时变系统的MAC值比较变化曲线，从图中MAC值变化曲线可以看出，在选取合适的遗忘因子后，模态振型的识别效果是较好的。

112.775s、600s、1000s、1500s、1987.2s各时刻模态置信度(MAC)比较如表1ˉ表5所示。

可以发现，各时刻识别效果比较好，振型非常接近；其中表5中的第二、三识别振型由于第二、三主成份贡献率非常接近而发生交换识别。

112.775s与1987.2s时刻的三阶固有频率比较结果，由此发现频率发生了变化，如表6所示。

根据表6，可得知本发明所述的装置为线性时变结构。

表1：112.775s时刻的固有振型与FRPCA识别振型的MCA比较

表2：600s时刻的固有振型与FRPCA识别振型的MCA比较

表3：1000s时刻的固有振型与FRPCA识别振型的MCA比较

表4：1500s时刻的固有振型与FRPCA识别振型的MCA比较

表5：1987.2s时刻的固有振型与FRPCA识别振型的MCA比较

表6：112.775s与1987.2s时刻的固有频率比较

本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。

本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。

Claims (2)

1.一种遗忘自协方差矩阵递推主元的时变工作模态识别方法，其特征在于，包括：

S101：获取线性时变结构在环境激励下多个振动传感器从初始时刻0到时刻k的非平稳时域振动响应信号数据矩阵

其中，表示维度为m×k的矩阵，m表示在所述线性时变结构上布置的振动传感器检测点个数，k表示时域采样点个数；1≤j≤m；1≤i≤k；

S102：将所述数据矩阵归一化，得到归一化矩阵X_k，求的均值向量与的标准差矩阵Σ_k＝diag{δ_k(1),δ_k(2),…,δ_k(j),…,δ_k(m)}；

S103：求取所述归一化矩阵X_k的自协方差矩阵 所述自协方差矩阵为实对称方阵；

S104：将所述自协方差矩阵分解为的特征值特征向量形式，其中， 为的特征值由大到小排列形成的对角矩阵， 依次为特征值的特征向量，V^(k)T为V^(k)的转置，V^(k)TV^(k)＝I_m×m；

S105：存储矩阵V^(k)，V^(k)对应k时刻该结构的瞬态正则化工作模态振型矩阵根据矩阵V^(k)与X_k，主成分V^(k)TX_k对应k时刻各阶模态响应矩阵利用单自由度模态识别技术对矩阵进行处理，得到时刻k的瞬时工作模态固有频率；

S106：获取下一时刻(k+1时刻)的振动数据

S107：获取新的时域振动响应数据矩阵

S108：将所述数据矩阵归一化，得到归一化矩阵X_k+1，引入新样本数据后，原始数据矩阵变为 的均值向量的标准差矩阵Σ_k+1＝diag{δ_k+1(1),δ_k+1(2),…,δ_k+1(j),…,δ_k+1(m)}和X_k+1的协方差矩阵按如下公式进行更新：

${\stackrel{\→}{b}}_{k+1}=\frac{k}{k+1}{\stackrel{\→}{b}}_k+\frac{1}{k+1}{\stackrel{\→}{x}}_{k+1}^0$

${\mathrm{\σ}}_{k+1}^2\left(j\right)=\frac{k}{k+1}\left({\mathrm{\σ}}_k^2\right(j)+{\mathrm{\Δb}}_{k+1}^2(j\left)\right)+\frac{1}{k+1}\left|\right|x_{k+1}^0\left(j\right)-b_{k+1}\left(j\right)|{|}^2$

${\stackrel{\→}{x}}_{k+1}={\mathrm{\Σ}}_{k+1}^{-1}({\stackrel{\→}{x}}_{k+1}^0-{\stackrel{\→}{b}}_{k+1})$

$C_{X_{k+1}{X_{k+1}}^T}=\frac{k}{k+1}{\mathrm{\Σ}}_{k+1}^{-1}({\mathrm{\Σ}}_k{C}_{X_k{X_k}^T}{\mathrm{\Σ}}_k+\mathrm{\Δ}{\stackrel{\→}{b}}_{k+1}\mathrm{\Δ}{\stackrel{\→}{b}}_{k+1}^T){\mathrm{\Σ}}_{k+1}^{-1}+\frac{1}{k+1}{\stackrel{\→}{x}}_{k+1}{\stackrel{\→}{x}}_{k+1}^T$

其中，Σ_k+1＝diag{δ_k+1(1),δ_k+1(2),…,δ_k+1(j),…,δ_k+1(m)}为的标准差矩阵；δ_k+1(j)为的标准差矩阵中的第j个元素，1≤j≤m；为的均值向量；为的均值向量；为均值向量的变化量；为k+1时刻的振动数据的归一化后的结果；

S109：引入遗忘因子γ，将新旧数据分以不同权重，则的均值向量 的标准差Σ_k+1＝diag{δ_k+1(1),δ_k+1(2),…,δ_k+1(j),…,δ_k+1(m)}的更新公式如下所示：

${\stackrel{\→}{b}}_{k+1}=\frac{\mathrm{\γ}k}{k+1}{\stackrel{\→}{b}}_k+(1-\frac{k\mathrm{\γ}}{k+1}){\stackrel{\→}{x}}_{k+1}^0$

${\mathrm{\σ}}_{k+1}^2\left(j\right)=\frac{k\mathrm{\γ}}{k+1}\left({\mathrm{\σ}}_k^2\right(j)+{\mathrm{\Δb}}_{k+1}^2(j\left)\right)+(1-\frac{k\mathrm{\γ}}{k+1})\left|\right|x_{k+1}^0\left(j\right)-b_{k+1}\left(j\right)|{|}^2$

${\stackrel{\→}{x}}_{k+1}={\mathrm{\Σ}}_{k+1}^{-1}({\stackrel{\→}{x}}_{k+1}^0-{\stackrel{\→}{b}}_{k+1})$

根据自协方差矩阵的递推公式得到所述归一化矩阵X_k+1的自协方差矩阵 为X_k+1的转置；

S110：将分解为的特征值特征向量形式，其中， 为的特征值由大到小排列形成的对角矩阵， 依次为特征值的特征向量，V^(k+1)T为V^(k+1)的转置，V^(k+1)TV^(k+1)＝I_m×m；

S111：存储矩阵V^(k+1)，V^(k+1)对应k+1时刻该结构的瞬态正则化工作模态振型矩阵

S112：根据矩阵V^(k+1)与X_k+1，主成分V^(k+1)TX_k+1对应k+1时刻各阶模态响应矩阵利用单自由度模态识别技术对进行处理，得到时刻k+1的瞬时工作模态固有频率；

S113：当k小于等于时域采样点总个数n时(k≤n)，将k+1的值代替k的值(k←k+1)，跳转至S106，否则，线性时变结构瞬态工作模态参数识别结束。

2.一种基于带遗忘因子加权和自协方差矩阵递推主元分析算法的线性时变结构的时变瞬态工作模态参数识别系统，其特征在于，所述系统包括：

非平稳时域振动响应数据获取单元：用于获取线性时变结构在环境激励下多个振动传感器从初始时刻0到时刻k的非平稳时域振动响应信号数据矩阵

其中，表示维度为m×k的矩阵，m表示在所述线性时变结构上布置的振动传感器检测点个数，k表示时域采样点个数；

矩阵归一化单元：用于将所述数据矩阵归一化，得到归一化矩阵X_k，求的均值向量与的标准差矩阵Σ_k＝diag{δ_k(1),δ_k(2),…,δ_k(j),…,δ_k(m)}；

自协方差矩阵求取单元：用于求取所述归一化矩阵X_k的自协方差矩阵 所述自协方差矩阵为实对称方阵；

自协方差矩阵分解单元：用于将所述自协方差矩阵分解为的特征值特征向量形式，其中， 为的特征值由大到小排列形成的对角矩阵， 依次为特征值的特征向量，V^(k)T为V^(k)的转置，V^(k)TV^(k)＝I_m×m；

矩阵存储单元：用于存储矩阵V^(k)，对应k时刻该结构的瞬态工作模态振型矩阵；根据矩阵V^(k)与X_k，利用单自由度模态识别技术对矩阵V^(k)TX_k进行处理，得到时刻k的瞬时工作模态固有频率；

下一时刻数据获取单元：用于获取下一时刻(k+1时刻)的振动响应数据

新时域振动响应数据获取单元：用于获取新的时域振动响应数据矩阵

新矩阵递归归一化单元：用于将所述数据矩阵归一化，得到归一化矩阵X_k+1，引入新样本数据后，原始数据矩阵变为 的均值向量的标准差矩阵Σ_k+1＝diag{δ_k+1(1),δ_k+1(2),…,δ_k+1(j),…,δ_k+1(m)}和X_k+1的协方差矩阵按如下公式进行更新：

${\stackrel{\→}{b}}_{k+1}=\frac{k}{k+1}{\stackrel{\→}{b}}_k+\frac{1}{k+1}{\stackrel{\→}{x}}_{k+1}^0$

${\mathrm{\σ}}_{k+1}^2\left(j\right)=\frac{k}{k+1}\left({\mathrm{\σ}}_k^2\right(j)+{\mathrm{\Δb}}_{k+1}^2(j\left)\right)+\frac{1}{k+1}\left|\right|x_{k+1}^0\left(j\right)-b_{k+1}\left(j\right)|{|}^2$

${\stackrel{\→}{x}}_{k+1}={\mathrm{\Σ}}_{k+1}^{-1}({\stackrel{\→}{x}}_{k+1}^0-{\stackrel{\→}{b}}_{k+1})$

$C_{X_{k+1}{X_{k+1}}^T}=\frac{k}{k+1}{\mathrm{\Σ}}_{k+1}^{-1}({\mathrm{\Σ}}_k{C}_{X_k{X_k}^T}{\mathrm{\Σ}}_k+\mathrm{\Δ}{\stackrel{\→}{b}}_{k+1}\mathrm{\Δ}{\stackrel{\→}{b}}_{k+1}^T){\mathrm{\Σ}}_{k+1}^{-1}+\frac{1}{k+1}{\stackrel{\→}{x}}_{k+1}{\stackrel{\→}{x}}_{k+1}^T$

其中，Σ_k+1＝diag{δ_k+1(1),δ_k+1(2),…,δ_k+1(j),…,δ_k+1(m)}为的标准差矩阵；δ_k+1(j)为的标准差矩阵中的第j个元素，1≤j≤m；为的均值向量；为的均值向量；为均值向量的变化量；为k+1时刻的振动数据的归一化后的结果；

新自协方差矩阵递归求取单元：引入遗忘因子γ，将新旧数据分以不同权重，则的均值向量的标准差Σ_k+1＝diag{δ_k+1(1),δ_k+1(2),…,δ_k+1(j),…,δ_k+1(m)}的更新公式如下所示：

${\stackrel{\→}{b}}_{k+1}=\frac{\mathrm{\γ}k}{k+1}{\stackrel{\→}{b}}_k+(1-\frac{k\mathrm{\γ}}{k+1}){\stackrel{\→}{x}}_{k+1}^0$

${\mathrm{\σ}}_{k+1}^2\left(j\right)=\frac{k\mathrm{\γ}}{k+1}\left({\mathrm{\σ}}_k^2\right(j)+\mathrm{\Δ}{\stackrel{\→}{b}}_{k+1}^2(j\left)\right)+(1-\frac{k\mathrm{\γ}}{k+1})\left|\right|x_{k+1}^0\left(j\right)-b_{k+1}\left(j\right)|{|}^2$

${\stackrel{\→}{x}}_{k+1}={\mathrm{\Σ}}_{k+1}^{-1}({\stackrel{\→}{x}}_{k+1}^0-{\stackrel{\→}{b}}_{k+1})$

根据自协方差矩阵的递推公式递推得到所述归一化矩阵X_k+1的自协方差矩阵 为X_k+1的转置；

新自协方差矩阵分解单元，用于将分解为的特征值特征向量形式，其中， 为的特征值由大到小排列形成的对角矩阵， 依次为特征值的特征向量，V^(k+1)T为V^(k+1)的转置，V^(k+1)TV^(k+1)＝I_m×m；

新矩阵存储单元，用于存储矩阵V^(k+1)，对应k+1时刻该结构的瞬态工作模态振型矩阵；

瞬时工作模态参数计算单元，根据矩阵V^(k+1)与X_k+1，利用单自由度模态识别技术对矩阵V^(k+1)TX_k+1进行处理，得到时刻k+1的瞬时工作模态固有频率；

跳转单元，用于当k小于等于时域采样点总个数n时(k≤n)，将k+1的值代替k的值(k←k+1)，跳转至下一时刻数据获取单元，否则，线性时变结构瞬时工作模态参数识别结束。