基于滑动窗变步长等变自适应源分离的线性慢时变结构工作模态参数识别

黄海阳，王 成+，李海波，赖雄鸣，缑 锦，张惠臻

(1.华侨大学 计算机科学与技术学院，福建 厦门 361021;2.华侨大学 机电及自动化学院，福建 厦门 361021)

0 引言

模态参数(模态固有频率、模态振型、模态阻尼比)能够很好地反映结构的动力学特性，因此提取模态参数成为结构系统辨识问题的核心[1]。然而，许多大规模的工程结构输入无法测量或代价昂贵，此时工作模态分析方法(又称仅输出辨识方法)能够在工作环境中仅从输出响应中提取结构的模态参数，因此得到广泛关注[2]。目前，工作模态分析已经广泛应用在故障诊断[3]、结构优化设计[4-5]等方面。

现实中工程结构大多是时变的，物理特性(质量、刚度、阻尼等)随时间而变化，其模态参数也会随着时间而变化[6]，例如，火车过桥，火箭发射，旋转机械等都呈现出时变特性[7]。目前，线性时变结构的模态参数识别方法主要分为时域法和频域法(时频分析法)，马志赛等[8]和Zhou等[9]对这些方法进行了综述，并总结了每种方法的优缺点;振动响应信号也不能够提前一次性获得，而是随着时间慢慢采样得到，因此需要在线识别方法[10];滑动窗方法是一种在线识别慢时变系统模态参数的方法，已经应用在一些算法上[11];滑动窗方法相比于递推的方法更能跟踪时变特性，因为其添加了新数据同时舍弃了旧数据，减轻了历史数据的影响[12]。

盲源分离(Blind Source Separation，BSS)是在源信号和混合方式未知的情况下，仅从观测信号中恢复源信号的方法[13]，包括独立成分分析、稀疏成分分析、二阶盲辨识等[14]。Sadhu等[15]对BSS方法在工作模态参数识别上的应用进行了综述，并指出几乎所有基于BSS的模态识别方法都是在离线或批处理模式下进行的，不适用于控制应用和参数的自动化识别。但识别时变结构的模态参数，需要在线的自适应的BSS方法,Luo等[16]对无线自适应处理BSS进行了综合研究;Cardoso等[17]提出了基于连续更新思想的等变自适应源分离(Equivariant Adaptive Separation via Independence, EASI)方法，该方法采用相对梯度，具有等变化性，迭代的在线更新每个时刻的分离矩阵;Amini等[18]将EASI应用到线性时不变结构的工作模态在线辨识中，并指出EASI具有收敛性依赖于源统计量和分析复杂的缺点。EASI是一种在线方法，能够节省存储空间，在数据输入的同时就进行处理，但是其使用的固定步长不能兼顾收敛性和跟踪能力[19]。为此，研究者们提出了大量的变步长EASI算法改进EASI算法的收敛性和稳定性。陆建涛等[20]提出基于分离指标的变步长EASI算法，利用一个非线性单调递增函数自适应调节步长;邬金松等[21]探索了分离程度的概念,并根据信号的分离程度自适应调整EASI算法的步长。但这些变步长方法都只是在时不变系统中提高算法性能，时变系统需要更复杂的变步长方法。目前，最新改进的变步长EASI算法在系统突变后能够自适应追踪并且重新收敛于新的平衡点[22-23],但当系统参数和源信号随时间持续变化时，变步长EASI算法也难以有效跟踪其特性。

综上所述，本文将滑动窗方法[24]与变步长EASI算法相结合，提出一种基于滑动窗变步长等变自适应源分离(Moving Window Variable step-size EASI, MWVEASI)的线性慢时变结构工作模态参数在线识别方法。利用滑动窗的在线跟踪特性和变步长EASI方法的收敛性，识别线性时变结构的瞬态模态固有频率和振型。

1 相关工作

1.1 传统盲源分离问题描述

BSS方法的目的是从观测信号X(t)∈N×T中恢复源信号S(t)∈m×T，混合过程A∈N×m和源信号都是未知的[25]。其中：N是传感器个数，T是采样点数，m是源信号个数。不考虑测量噪声，可以表示为

X(t)=AS(t)。

(1)

为了从观测信号中分离出独立源信号，需要寻找分离矩阵B∈m×N，使得

Y(t)=BX(t)=BAS(t)。

(2)

式中Y(t)=[y1(t),…,ym(t)]T∈m×T是源信号S(t)∈m×T的估计。如图1所示为盲源混合和分离的过程。

1.2 EASI算法理论描述

传统BSS方法为离线和批处理的算法，且只适用于平稳环境。因此，研究者们提出了自适应BSS方法，其中，Cardoso等[17]在1996年提出EASI算法，该算法基于连续更新的想法，采用相对梯度，通过迭代得到分离矩阵B(t+1)。

(3)

式中：I∈m×m为单位矩阵;λ(t)为第t时刻的步长;B(t+1)∈m×N取经过多次迭代收敛的t+1时刻分离矩阵。

(4)

式中：V为特征X(t)的特征向量;U为对角矩阵，其对角线由X(t)的特征值组成。

传统的EASI算法式(3)中的λ(t)是固定值。近年来，学者们提出了一系列的变步长方法，本文选取指数衰减型[26]和时间递减型[27]。步长随着时间t增加而减小，分别如式(4)和式(5)所示，使得EASI算法能够兼顾稳定性和收敛性：

(5)

(6)

1.3 基于EASI算法的线性时不变结构工作模态参数识别

线性时不变结构的质量M、阻尼C和刚度矩阵K不随时间改变，其动力学方程为：

(7)

X(t)=ΦQ(t)。

(8)

2 基于MWVEASI算法的线性慢时变结构工作模态参数在线识别

2.1 线性慢时变结构在模态坐标下的表示

根据动力学理论，物理坐标下线性时变结构的动力学方程为：

=F(t),t∈[tbegin,tend]。

(9)

式中：M(t),C(t)和K(t)∈N×N分别表示时变的质量、阻尼和刚度矩阵。基于“时间冻结”和短时时不变理论，线性时变结构在时间[tbegin,tend]内，质量、刚度和阻尼可视为是“冻结”(时不变)的，因而，式(9)可表示为一系列时不变结构S′(ε)构成的集合S′：

(10)

(11)

(12)

(13)

在每个窗内，该时段内的模态振型和固有频率可作为最中间时刻的瞬时模态参数。

2.2 MWVEASI算法

滑窗(moving window)技术利用短时时不变理论，使得在每个窗内非平稳振动响应信号可以看作是平稳的,线性时变结构可以看成是时不变的。随着时间的推移，窗口向右滑动，增加新数据同时删除旧数据，以此类推，如图3所示描述了该过程。

将滑动窗和自适应BSS的变步长EASI算法相结合，提出滑动窗变步长EASI(MWVEASI)的时变BSS算法。在每个窗内，应用变步长EASI方法进行自适应盲源分离。

2.3 基于变步长MWVEASI算法的线性时变结构工作模态参数在线识别方法

在图4的流程中，首先将响应数据分成一系列数据窗。然后，在每个窗内，每个时刻采用不同步长的EASI方法进行迭代。将窗内最后时刻的混合矩阵和源信号估计的结果看成是窗内中间时刻的结果。接着，对EASI识别的源信号使用快速傅里叶变换(Fast Fourier Transform, FFT)或单自由度技术(如自然激励法(Natural Excitation Technique, NExT)等)得到模态固有频率和阻尼比。

2.4 方法的适用范围和理论分析比较

基于滑动窗变步长EASI的线性慢时变结构工作模态参数识别方法的适用范围如下：

(1)系统是线性慢时变的弱阻尼结构；

(2)振动响应传感器个数大于或等于可识别工作模态参数阶数；

(3)振动响应信号的采样频率大于或等于可识别模态固有频率的2倍；

(4)系统未知环境激励应该近似为平稳白噪声。

EASI方法只能识别线性时不变结构的模态参数[19]。变步长的EASI方法通过可变的步长使得算法可以迅速收敛[20]。当系统突变时，变步长EASI方法依旧可以收敛到平衡位置[22-23]。然而，一方面变步长EASI方法重新收敛需要时间，另一方面当时变结构参数随时间持续变化时，EASI方法此时的步长不仅与时间相关，还应与时变特性有关。只随时间变化的步长将导致变步长EASI跟踪时变的能力弱、识别精度低。滑动窗方法将当前时刻数据和前后一定时延时刻的数据组成数据窗，在窗内，基于“时间冻结理论”，数据非平稳程度低，使得窗内变步长EASI方法有较高的精度。滑动窗过程加入了新数据并剔除旧数据，窗内数据包含了当前的时变信息，因而滑动窗能够跟踪时变结构的变化[24]。滑动窗并不改变EASI算法本身的收敛性，因此滑窗变步长EASI收敛速度比滑动窗固定步长EASI快。综上所述，表1描述了EASI方法和滑动窗方法相结合后的特点。

表1 基于滑动窗变步长EASI与其他线性慢时变结构工作模态参数识别的比较

3 仿真验证

3.1 仿真系统和数据集介绍

为验证所提出的MWVEASI算法的正确性和精度，在MATLAB/Simulink中，设计了一个线性时变的三自由度弹簧振子，如图5所示。

图5所示的三自由度弹簧结构在慢时变条件下的动力学方程为：

(14)

设置3个物块的初始位移为零，阻尼比都为c1=c2=c3=0.01 N.s/m，刚度都是k1=k2=k3=1 000 N/m。质量m2=m3=1 kg，m1随时间缓慢变化如下：

(15)

m1(t)所受到的激励F1(t)为高斯白噪声，采样间隔为0.025 s，采样频率为40 Hz，仿真时长为2 000 s。取50 s后的时变数据进行研究，三自由度结构在MATLAB/Simulink的建模如图6所示。在Simulink模块中利用Runge-Kutta算法求解得到位移响应信号X(t)数据集，如图7所示。更具体的仿真细节可见文献[24]。最终得到的位移响应信号数据集是一个3行80 001列的矩阵，即N=3，T=80 001。

3.2 评价指标

理论的模态固有频率、振型、阻尼比通过有限元仿真得到。采用模态置信度准则(MAC)作为评价指标衡量模态振型识别的精度，定义如下：

(16)

其中：φi为第i阶理论值，φj为第j阶识别值。0≤MAC≤1，MAC越接近1说明精度越高。

3.3 参数设置

滑动窗的窗长L取1 024。变步长EASI方法的参数取值，即式(4)和式(5)分别为λ0=0.014,t0=200,α=0.0025,β=2。固定步长EASI方法取步长为0.001 4。B的初值设为：B0=0.6×rands(3,3)，其中rands()是MATLAB的一个函数，表示生成3行3列的每个元素是0～1之间随机数的矩阵。ε=10-4，l_max=5000，t=l=i=1。

收敛和停机准则：

(3)当t=L时，i=i+1，进入下一个窗，重置EASI方法初值，t=l=1，返回步骤(1)，直到i+L>T时停止。

3.4 工作模态参数识别结果

本文对两种变步长EASI算法、相应的两种加滑动窗变步长方法以及滑窗固定步长EASI的结果进行了对比，图8和图9分别表示5种算法的频率、MAC值识别结果与理论值的对比。表2统计了各种算法未识别的模态参数个数。

表2 未识别的模态参数百分比

3.5 识别结果分析

(1)从图8b、图8d、图9b、图9d可以看出，基于滑动窗EASI方法能够识别线性慢时变结构的模态参数。而图8a、图8c、图9a、图9c则表明，不加窗的变步长EASI算法在识别时变快的第三阶模态参数时效果不佳。对比图8e与8a、图8c，滑窗EASI方法时变快时识别的值仍有较高的精度，因而跟踪性比变步长EASI好。

(2)对比图8e与8b、图8d，可以看出，滑窗固定步长EASI方法识别模态参数精度没有滑窗变步长EASI方法高。从图8b、图8d、图9b、图9d可以看出，基于滑动窗指数衰减变步长EASI算法识别精度比基于滑动窗迭代递减变步长EASI高。

(3)从图8、图9和表2可以看出，MWVEASI算法未能识别的模态参数时刻点的个数明显少于MWEASI和VEASI，说明MWVEASI算法收敛速度快，更容易收敛、稳定性更强。解决这些未识别的时刻点的办法是将上一时刻点的分离矩阵收敛值作为该时刻点分离矩阵的迭代计算初值，以及改进变步长格式和步长参数，加快算法收敛速度。

4 结束语

本文将滑动窗方法与变步长EASI方法相结合，应用于识别线性慢时变结构的工作模态参数。在质量慢时变的三自由度弹簧振子结构的仿真数据集上的识别结果表明，两种不同的MWVEASI方法都能够很好地识别线性慢时变结构的模态振型、固有频率，与滑窗固定步长EASI、变步长EASI算法相比，MWVEASI能够更好地跟踪时变特性、识别精度高、收敛速度快。

但是，变步长也存在方式和参数选择的问题；滑动窗方法也带来了计算量的增加和窗长设置的问题；如何实现窗与窗之间的递推，减少时空复杂度也是进一步的研究方向。与此同时，通过具体的工程实例的实验数据来验证所提出的新方法也是进一步的工作方向。