CN107357977B - 基于二阶盲辨识的线性结构工作模态参数识别方法及装置 - Google Patents

Abstract

本发明涉及一种基于二阶盲辨识的线性时不变结构工作模态参数识别方法及对应的时不变的工作模态参数识别的装置、一种结合二阶盲辨识和最小二乘广义逆方法的三维结构的工作模态参数识别方法及对应一种时不变三维圆柱壳工作模态参数识别的实验装置，一种时变的工作模态参数识别的方法和一种基于滑动窗二阶盲辨识的线性时变结构工作模态参数识别装置。该方法和装置仅由实测的响应信号就能对时不变或者带有时变特性的动态系统进行工作模态参数在线实时识别，识别出结构(系统)的工作模态参数(模态振型，模态频率)，能实时有效监测系统的动态变化特性，可用于振动控制、设备故障诊断、健康监测以及系统结构分析与优化。

Description

基于二阶盲辨识的线性结构工作模态参数识别方法及装置

技术领域

本发明涉及一种基于二阶盲辨识的线性时不变结构工作模态参数识别方法，在一种时不变的工作模态参数识别的装置产生的数据上进行方法验证。一种结合二阶盲辨识和最小二乘广义逆方法的三维结构的工作模态参数识别方法，在一种时不变三维圆柱壳工作模态参数识别的实验装置产生的数据上进行方法验证。一种基于滑动窗二阶盲辨识的线性时变结构工作模态参数识别方法，以及该方法在仿真数据上进行方法验证。

背景技术

工作模态分析是仅利用结构在工作时而不是在实验中利用传感器获得的振动响应数据中提取出模态参数，其中模态参数包括模态阻尼比，模态振型和模态固有频率。目前工作模态分析方法在解决工程问题上是非常重要的分析方法，特别是在航空，土木工程等领域。

然而，对于大型结构的模态参数不仅仅受到外激励和内激励的影响，而且在工作状态中很难保持不变，也就是说，结构的模态参数是随着时间的变化而变化的，不是时不变的，而是时变的。因此，要识别时变结构的模态参数，如果仅仅依靠时不变的方法是不可能实现的，还得借助时变的方法。国内外学者对时变结构的工作模态参数识别也做了大量的研究，比如有提出TARMA模型，利用小波变换和希尔伯特变换来识别时变结构的模态参数。

发明内容

本发明的目的在于克服现有技术的不足，提供一种基于滑动窗二阶盲辨识的线性时变结构工作模态参数识别方法，还提供基于滑动窗二阶盲辨识的线性时变结构工作模态参数识别方法在设备故障诊断与健康状态监测上的应用。

本发明的技术方案如下：

一种基于二阶盲辨识的线性时不变一维结构工作模态参数识别方法，利用二阶盲辨识方法分解线性时不变一维结构的振动响应信号，从而得到时不变结构的各阶模态参数，其中在二阶盲辨识的解混旋转阶段，利用多个时延非零的相关矩阵同时对角化的方法，消除了由于单个时延选择不适而造成的算法结果的不理想，从而有效识别出时不变一维结构工作模态参数。

所述的基于二阶盲辨识的线性时不变一维结构工作模态参数识别方法，具体步骤如下：

步骤1)时不变结构的模态参数随时间的变化而变化，根据结构动力学理论，对于自由度线性时变振动结构系统，它在物理坐标系统中的运动方程为：

其中，M，C和

分别表示质量矩阵，阻尼矩阵和刚度矩阵，与此同时，它们受结构的影响而随时间发生变化；

表示外载荷的激励向量，

和

分别表示加速度响应信号，速度响应信号和位移响应信号；

表示维度为n×T的矩阵；

步骤2)二阶盲辨识是另外一种盲源分离算法，它利用了信号的另外一种特性—时序结构，即时延协方差信息，要求信号具有不同的谱特性或者不同的相关函数，它的处理对象是时间信号，当然也可以处理统计独立的源信号，为了正确地求解出分离矩阵和源信号，二阶盲辨识算法做出了如下的假设：

假设一：混合矩阵列满秩；

假设二：源信号相互不相关且具有不同的自相关函数；

假设三：源信号是平稳信号；

步骤3)根据上述假设，源信号的协方差矩阵R_s(0)满足如下：

R_s(0)＝E[S(t)S^T(t)]＝I

其中，E[*]表示期望，S(t)表示源信号，I表示单位矩阵。

步骤4)，观测信号的协方差矩阵可表示如下：

其中，X_bss(t)表示观测信号，A表示混合矩阵。

步骤5)对观测信号X_bss(t)进行白化预处理，即进行线性变换：

Z(t)＝VX_bss(t)

其中，

是一个白化矩阵，线性变换的目的是使得协方差矩阵

是一个单位矩阵，去除信号之间的相关性；

步骤6)源信号的延时协方差矩阵R_s(τ)类推定义如下：

R_s(τ)＝E[S(t+τ)S^T(t)]

步骤7)与此同时，源信号白化后的延时协方差矩阵可得到：

R_z(τ)＝E[Z(t+τ)Z^T(t)]

步骤8)选择一组不同的值τ₁,τ₂,…,τ_p，可得到一系列延时协方差矩阵R_s(τ_i)，其中i＝1,…,p，可以得出如下：

Z(t)＝VX_bss(t)＝(VA)S(t)＝US(t)

步骤9)其中，U＝VA，因此，可推导出如下：

R_z(τ_i)＝UR_s(τ_i)U^T,i＝1,2,…,p

步骤10)由于源信号被假设成不相关的，而且源信号白化后的矩阵Z(t)中的各向量相互正交归一化，因此，矩阵U也是正交归一化的，而且

R_s(τ_i)＝U^TR_z(τ_i)U,i＝1,2,…,p

步骤11)可以看出，R_s(τ_i)是一个对角化矩阵，可以看出矩阵U是归一化正交的，因此，利用R_z(τ_i)和最优逼近算法来寻找最优的矩阵U，可计算出分离矩阵W如下：

W＝U^TV

步骤12)同时，源信号S(t)也可计算得到。

一种时不变的工作模态参数识别的装置，包括一根一端简支一端固支的金属梁、n个加速度传感器、数据采集卡、激振器或力锤及电脑终端；所述金属梁设置有n个等分点，其中每等分点上放置有一个加速度传感器，用于采集使用所述激振器或力锤对所述金属梁的某一点施加激励所产生的振动响应信号，所述数据采集卡与所述n个加速度传感器相连用于接收所述振动响应信号并发送给所述电脑终端存储，所述电脑终端包括一模态参数识别模块，所述模态参数识别模块基于电脑终端存储的振动响应信号使用基于二阶盲辨识的线性时不变一维结构工作模态参数识别方法识别出金属梁的模态参数，并与梁的理论模态参数作对比以验证基于二阶盲辨识的线性时不变结构工作模态参数识别方法的正确性。

在实际工程应用中，一般的工程结构都是三维结构，因此，从一维结构到三维结构的模态参数分析，是目前的一大挑战，本发明一种结合二阶盲辨识和最小二乘广义逆方法的三维结构的工作模态参数识别方法，为了求解三维连续体复杂结构的工作模态参数识别，先从其中的一个方向入手，由于每个方向的振动响应都不同，以二阶盲辨识算法为基础，选取响应最大的那个方向用二阶盲辨识进行工作模态参数识别，此时的做法跟一维结构的工作模态参数识别是一致的；当分解出响应最大的那个方向的振动响应，可得到模态响应矩阵和一维模态振型，然后将模态响应矩阵使用最小二乘广义伪逆，回代到另外两个方向的位移响应信号中，然后组装各方向的模态振型，形成识别三维模态振型的算法；最后利用单自由度系统参数识别技术快速傅里叶变换，从模态响应矩阵中得到包括模态固有频率的信息。

所述的结合二阶盲辨识和最小二乘广义逆方法的三维结构的工作模态参数识别方法，具体步骤如下：

步骤1)首先，利用二阶盲辨识方法分解X方向上的振动响应，分解得到X方向上的模态振型Ψ_d×l和结构的模态响应矩阵H_l×T(t)：

(X_three)_d×T(t)≈Ψ_d×lH_l×T(t)

其中(X_three)_d×T(t)表示X方向上的振动响应信号，d为三维结构的阶数，T为时间。

步骤2)由于三个方向的模态响应矩阵是相同的，利用基于二阶盲辨识的三维模态参数识别算法对X方向上进行二阶盲辨识分解之后，需要回代到另外两个方向，模态响应是个矩阵，实际中，是对另外两个方向的响应乘以模态响应矩阵的右伪逆来实现，最小二乘广义逆是在误差平方和最小意义下的一致无偏最优估计，因此利用最小二乘广义伪逆的方法求解Y和Z方向上的模态振型：

其中(Y_three)_d×T(t)和(Z_three)_d×T(t)分别表示Y方向上和Z方向上的振动响应信号，d为三维结构的阶数，T为时间。O_d×l和B_d×l分别表示Y和Z方向上的模态振型矩阵。

步骤3)然后组装各方向的模态振型，形成识别三维模态振型的算法；最后利用单自由度系统参数识别技术快速傅里叶变换，可从模态响应矩阵中得到包括模态固有频率的信息。

一种时不变三维圆柱壳工作模态参数识别的实验装置，包括一个一端固支一端自由的三维薄壁金属圆柱壳、n个三向加速度传感器、数据采集卡、激振器或力锤及电脑终端；所述金属圆柱壳表面上均匀放置n个三向加速度传感器，用于采集使用所述力锤或激振器垂直于圆柱壳表面施加激励产生的三个方向的振动响应信号，所述数据采集卡与所述n个三向加速度传感器相连用于接收所述振动响应信号并发送给所述电脑终端存储，所述电脑终端包括一模态参数识别模块，所述模态参数识别模块基于电脑终端存储的振动响应信号使用结合二阶盲辨识和最小二乘广义逆方法的三维结构的工作模态参数识别方法识别三维圆柱壳的模态参数，并与由有限元分析或理论解相比较，验证结合二阶盲辨识和最小二乘广义逆方法的三维工作模态参数识别方法的正确性。

一种基于滑动窗二阶盲辨识的线性时变一维结构工作模态参数识别方法，其分析的结构是线性时变结构，模态参数是随时间的改变而发生变化的，因此主要结合短时时不变理论与二阶盲辨识算法，利用二阶盲辨识算法在各窗内的统计特性，估计出各时刻的工作模态参数，然后各时刻求得的工作模态参数连接起来，从而实现时变线性结构工作模态参数识别；所述工作模态参数包括各阶模态的固有频率和模态振型。

所述的基于滑动窗二阶盲辨识的线性时变一维结构工作模态参数识别方法，具体步骤如下：

步骤1)时变结构的模态参数随时间的变化而变化，根据结构动力学理论，对于自由度线性时变振动结构系统，它在物理坐标系统中的运动方程为:

其中，M(t)，C(t)和

分别表示随时间变化的质量矩阵，阻尼矩阵和刚度矩阵，与此同时，它们受结构的影响而随时间发生变化；

表示外载荷的激励向量，

和

分别表示加速度响应信号，速度响应信号和位移响应信号；

步骤2)根据“短时时不变”理论，时变离散多自由度系统在时间τ∈[t_begin,t_end]内，它的质量，阻尼和刚度看作是时不变的，因此，在物理坐标系中的动力学方程可表示为：

其中，K＝end，表示时变系统中最后时刻，S'(τ)表示为当t＝τ时刻的时不变结构，S'表示一组有限多个线性时不变结构组成时变结构的集合；

步骤3)对于时不变小阻尼结构，响应数据可以被分为有限多个部分，在第τ个部分，选取一定的窗口的长度，线性系统的模态坐标响应为：

其中，Φ(τ)和q(τ,t)分别表示当第τ个窗口的模态振型矩阵和模态响应向量；

步骤4)当时不变结构的每一阶模态固有频率ω_i都不相等时，各阶模态振型之间满足归一化正交，各阶模态响应相互不相关，如下：

其中，E(q(τ,t)q^T(τ,t))表示两模态振型的期望，Λ"_n×n表示阶数为n的对角矩阵；

步骤5)假设在一个很短的时间段内，系统被看成是短时时不变的，也就是说，时间被划分成有限段，在每一个时间段内，系统被认为是短时时不变的，从而利用时不变结构的工作模态参数识别算法，识别出该时间段的工作模态参数，窗口向右边滑动，即计算下一个时间段内工作模态参数，以此类推，最后将每个时间段按照时间顺序排列起来，从而形成时变结构的模态参数；

其中，响应数据的限定记忆长度为L，n表示的是传感器的个数，T表示采样时间；

步骤6)对于采集到的振动结构的位移响应数据，其模态坐标表示如下：

其中，

表示振动结构的位移响应数据，

是模态振型矩阵，

表示模态坐标响应，当系统的各阶的模态固有频率不相等是，各阶模态振型相互归一化正交，各阶模态坐标响应相互独立，如下：

其中，

表示两模态坐标响应的期望，Λ"_n×n表示阶数为n的对角矩阵；

步骤7)对于时变结构系统，

是由模态振型向量

组成的，并代表了线性时变结构在L时间段内的统计平均模态振型，可近似估计求得时变结构在(i+(L-1)/2)时刻的瞬间模态振型

同时，

由模态响应函数

构成，并代表了线性时变系统在L时段内的统计平均模态响应，利用单自由度识别技术，可近似估计求得在(i+(L-1)/2)时刻的瞬时模态频率ω_j(i+(L-1)/2)；

步骤8)所识别的模态固有频率是通过快速傅里叶变换得到的，模态固有频率的识别精度取决于频率分辨率△f，当快速傅里叶变换的长度L越长，模态固有频率的识别精度越高，而且频率分辨率△f与采样频率f呈现正比例的关系，频率分辨率△f可定义为如下：

△f＝f/L

步骤9)当线性时变结构的模态参数发生变化时和振动响应信号是非平稳的，如果线性时变结构的模态参数变化非常快，这是一个非平稳程度非常高的时变系统，此时，应该减少滑动窗窗口的长度；定义在一个数据窗内第i阶模态的平均频率变化量为△f_L(i)：

其中，变量f_end(i),f_begin(i),t_end,t_begin分别表示第i阶模态的结束频率，开始频率，整个数据的结束时间和开始时间；

步骤10)滑动窗的窗口的长度L是振动响应数据的采样频率和频率分辨率相除得到的，△f的取值不能取太小，因为不能反映出频率的变化，△f的取值也不能取太大，因为平均频率变化量△f_L(i)不能被识别出来，与此同时，△f_L(i)不能比△f的取值大太多，否则，线性时变结构在一个滑动窗的窗口的长度L内不能被视为时不变的了；

步骤11)模态固有频率的计算，是由单自由度技术计算得到的，利用的是快速傅里叶变换，而快速傅里叶变换的计算复杂度跟滑动窗窗口的长度有直接的关系，根据快速傅里叶变换的原理，滑动窗窗口的长度应该满足如下：

L＝2^α,α＝1,2,…。

本发明的有益效果如下：

本发明所述的一种基于二阶盲辨识的线性时不变一维结构工作模态参数识别方法，利用二阶盲辨识方法分解线性时不变一维结构的振动响应信号，从而得到时不变结构的各阶模态参数，其中在二阶盲辨识的解混旋转阶段，利用多个时延非零的相关矩阵同时对角化的方法，消除了由于单个时延选择不适而造成的算法结果的不理想，从而有效识别出时不变一维结构工作模态参数。

本发明所述的一种结合二阶盲辨识和最小二乘广义逆方法的三维结构的工作模态参数识别方法，以二阶盲辨识算法为基础，选取响应最大的那个方向用二阶盲辨识进行工作模态参数识别，可得到模态响应矩阵和一维模态振型，然后将模态响应矩阵使用最小二乘广义伪逆，回代到另外两个方向的位移响应信号中，然后组装各方向的模态振型，形成识别三维模态振型的算法。最后利用单自由度系统参数识别技术比如快速傅里叶变换等，可从模态响应矩阵中得到模态固有频率等信息。

本发明所述的一种基于滑动窗二阶盲辨识的线性时变结构工作模态参数识别方法，能对带有时变特性的结构进行实时在线的参数识别，识别出系统的工作模态参数(模态振型，模态频率)，实时有效监测系统的动态变化特性，可被用于设备故障诊断、健康监测以及系统结构分析与优化。且该方法是一种工作模态参数识别方法(仅由实测响应信号即可识别出系统的特性)，并从数学理论分析及实验上给予证明，赋予了该方法以物理解释，较之于传统的需要同时测量激励与响应信号的试验模态参数识别技术具有较大的优势。该方法主要思想是，结合短时时不变理论与二阶盲辨识算法，利用二阶盲辨识算法在各窗内的统计特性，估计出各时刻的工作模态参数(包括各阶模态的固有频率和模态振型)，然后各时刻求得的工作模态参数连接起来，从而实现时变线性结构工作模态参数识别。

附图说明

图1是时不变结构工作模态分析装置框图；

图2是时不变简支梁的结构；

图3是时不变结构工作模态分析实验现场布置；

图4是采样频率为2k时SOBI方法识别的固有频率；

图5是采样频率为2k时SOBI方法识别的模态振型；

图6是采样频率为5k时SOBI方法识别的固有频率；

图7是采样频率为5k时SOBI方法识别的模态振型；

图8是采样频率为10k时SOBI方法识别的固有频率；

图9是采样频率为10k时SOBI方法识别的模态振型；

图10是第1118个观测点上三个方向的位移响应数据；

图11是圆柱壳的实际模态振型；

图12是模态阻尼比为0.03时，SOBI和最小二乘广义伪逆方法识别的固有频率；

图13是模态阻尼比为0.03且加5％噪声，SOBI和最小二乘广义伪逆方法识别的固有频率；

图14为当模态阻尼比η＝0.03时三种方法识别的模态振型的形状；

图15是模态阻尼比为0.1时，SOBI和最小二乘广义伪逆方法识别的固有频率；

图16是模态阻尼比为0.1且加5％噪声，SOBI和最小二乘广义伪逆方法识别的固有频率；

图17为当模态阻尼比η＝0.1时三种方法识别的模态振型的形状；

图18是时变三自由度的弹簧振子模型；

图19是高斯白噪声激励和三个响应数据；

图20是三阶固有频率的变化；

图21是在4个不同时刻的真实模态振型；

图22是利用SOBI方法识别的三阶模态固有频率；

图23是实际模态振型和识别的三阶模态振型的比较；

图24是将实际频率和识别频率进行比较；

图25是在四个时刻的模态振型；

图26是三阶模态振型的MAC值；

图27是将实际频率和识别频率进行比较；

图28是四个时刻的三阶模态振型；

图29是将实际频率和识别频率进行比较；

图30是刚性悬臂梁结构有限元模型；

图31是密度变化率为0.005，滑动窗窗口长度为2048时的识别结果；

图32是密度变化率为0.005，滑动窗窗口长度为4096时的识别结果；

图33是密度变化率为0.005，滑动窗窗口长度为8192时的识别结果；

图34是密度变化率为0.08，滑动窗窗口长度为2048时的识别结果；

图35是密度变化率为0.08，滑动窗窗口长度为4096时的识别结果；

图36是密度变化率为0.08，滑动窗窗口长度为8192时的识别结果。

具体实施方式

以下结合附图及实施例对本发明进行进一步的详细说明。

一种基于二阶盲辨识的线性时不变一维结构工作模态参数识别方法，利用二阶盲辨识方法分解线性时不变一维结构的振动响应信号，从而得到时不变结构的各阶模态参数，其中在二阶盲辨识的解混旋转阶段，利用多个时延非零的相关矩阵同时对角化的方法，消除了由于单个时延选择不适而造成的算法结果的不理想，从而有效识别出时不变一维结构工作模态参数。具体包括：

步骤1)时不变结构的模态参数随时间的变化而变化，根据结构动力学理论，对于自由度线性时变振动结构系统，它在物理坐标系统中的运动方程为：

其中，M，C和

分别表示质量矩阵，阻尼矩阵和刚度矩阵，与此同时，它们受结构的影响而随时间发生变化；

表示外载荷的激励向量，

和

分别表示加速度响应信号，速度响应信号和位移响应信号；

表示维度为n×T的矩阵；

步骤2)二阶盲辨识是另外一种盲源分离算法，它利用了信号的另外一种特性—时序结构，即时延协方差信息，要求信号具有不同的谱特性或者不同的相关函数，它的处理对象是时间信号，当然也可以处理统计独立的源信号，为了正确地求解出分离矩阵和源信号，二阶盲辨识算法做出了如下的假设：

假设一：混合矩阵列满秩；

假设二：源信号相互不相关且具有不同的自相关函数；

假设三：源信号是平稳信号；

步骤3)根据上述假设，源信号的协方差矩阵R_s(0)满足如下：

R_s(0)＝E[S(t)S^T(t)]＝I

其中，E[*]表示期望，S(t)表示源信号，I表示单位矩阵。

步骤4)，观测信号的协方差矩阵可表示如下：

其中，X_bss(t)表示观测信号，A表示混合矩阵。

步骤5)对观测信号X_bss(t)进行白化预处理，即进行线性变换：

Z(t)＝VX_bss(t)

其中，

是一个白化矩阵，线性变换的目的是使得协方差矩阵

是一个单位矩阵，去除信号之间的相关性；

步骤6)源信号的延时协方差矩阵R_s(τ)类推定义如下：

R_s(τ)＝E[S(t+τ)S^T(t)]

步骤7)与此同时，源信号白化后的延时协方差矩阵可得到：

R_z(τ)＝E[Z(t+τ)Z^T(t)]

步骤8)选择一组不同的值τ₁,τ₂,…,τ_p，可得到一系列延时协方差矩阵R_s(τ_i)，其中i＝1,…,p，可以得出如下：

Z(t)＝VX_bss(t)＝(VA)S(t)＝US(t)

步骤9)其中，U＝VA，因此，可推导出如下：

R_z(τ_i)＝UR_s(τ_i)U^T,i＝1,2,…,p

步骤10)由于源信号被假设成不相关的，而且源信号白化后的矩阵Z(t)中的各向量相互正交归一化，因此，矩阵U也是正交归一化的，而且

R_s(τ_i)＝U^TR_z(τ_i)U,i＝1,2,…,p

步骤11)可以看出，R_s(τ_i)是一个对角化矩阵，可以看出矩阵U是归一化正交的，因此，利用R_z(τ_i)和最优逼近算法来寻找最优的矩阵U，可计算出分离矩阵W如下：

W＝U^TV

步骤12)同时，源信号S(t)也可计算得到。

一种结合二阶盲辨识和最小二乘广义逆方法的三维结构的工作模态参数识别方法，以二阶盲辨识算法为基础，选取响应最大的那个方向用二阶盲辨识进行工作模态参数识别，可得到模态响应矩阵和一维模态振型，然后将模态响应矩阵使用最小二乘广义伪逆，回代到另外两个方向的位移响应信号中，然后组装各方向的模态振型，形成识别三维模态振型的算法。最后利用单自由度系统参数识别技术比如快速傅里叶变换等，可从模态响应矩阵中得到模态固有频率等信息，具体包括：

步骤1)首先，利用二阶盲辨识方法分解X方向上的振动响应，分解得到X方向上的模态振型Ψ_d×l和结构的模态响应矩阵H_l×T(t)：

(X_three)_d×T(t)≈Ψ_d×lH_l×T(t)

其中(X_three)_d×T(t)表示X方向上的振动响应信号，d为三维结构的阶数，T为时间。

步骤2)由于三个方向的模态响应矩阵是相同的，利用基于二阶盲辨识的三维模态参数识别算法对X方向上进行二阶盲辨识分解之后，需要回代到另外两个方向，模态响应是个矩阵，实际中，是对另外两个方向的响应乘以模态响应矩阵的右伪逆来实现，最小二乘广义逆是在误差平方和最小意义下的一致无偏最优估计，因此利用最小二乘广义伪逆的方法求解Y和Z方向上的模态振型：

其中(Y_three)_d×T(t)和(Z_three)_d×T(t)分别表示Y方向上和Z方向上的振动响应信号，d为三维结构的阶数，T为时间。O_d×l和B_d×l分别表示Y和Z方向上的模态振型矩阵。

步骤3)然后组装各方向的模态振型，形成识别三维模态振型的算法；最后利用单自由度系统参数识别技术快速傅里叶变换，可从模态响应矩阵中得到包括模态固有频率的信息。

一种基于滑动窗二阶盲辨识的线性时变结构工作模态参数识别方法，结合短时时不变理论与二阶盲辨识算法，利用二阶盲辨识算法在各窗内的统计特性，估计出各时刻的工作模态参数(包括各阶模态的固有频率和模态振型)，然后各时刻求得的工作模态参数连接起来，从而实现时变线性结构工作模态参数识别。

具体步骤如下：

步骤1)时变结构的模态参数随时间的变化而变化，根据结构动力学理论，对于自由度线性时变振动结构系统，它在物理坐标系统中的运动方程为:

其中，M(t)，C(t)和

分别表示随时间变化的质量矩阵，阻尼矩阵和刚度矩阵，与此同时，它们受结构的影响而随时间发生变化；

表示外载荷的激励向量，

和

分别表示加速度响应信号，速度响应信号和位移响应信号；

步骤2)根据“短时时不变”理论，时变离散多自由度系统在时间τ∈[t_begin,t_end]内，它的质量，阻尼和刚度看作是时不变的，因此，在物理坐标系中的动力学方程可表示为：

其中，K＝end，表示时变系统中最后时刻，S'(τ)表示为当t＝τ时刻的时不变结构，S'表示一组有限多个线性时不变结构组成时变结构的集合；

步骤3)对于时不变小阻尼结构，响应数据可以被分为有限多个部分，在第τ个部分，选取一定的窗口的长度，线性系统的模态坐标响应为：

其中，Φ(τ)和q(τ,t)分别表示当第τ个窗口的模态振型矩阵和模态响应向量；

步骤4)当时不变结构的每一阶模态固有频率ω_i都不相等时，各阶模态振型之间满足归一化正交，各阶模态响应相互不相关，如下：

其中，E(q(τ,t)q^T(τ,t))表示两模态振型的期望，Λ"_n×n表示阶数为n的对角矩阵；

步骤5)假设在一个很短的时间段内，系统被看成是短时时不变的，也就是说，时间被划分成有限段，在每一个时间段内，系统被认为是短时时不变的，从而利用时不变结构的工作模态参数识别算法，识别出该时间段的工作模态参数，窗口向右边滑动，即计算下一个时间段内工作模态参数，以此类推，最后将每个时间段按照时间顺序排列起来，从而形成时变结构的模态参数；

其中，响应数据的限定记忆长度为L，n表示的是传感器的个数，T表示采样时间；

步骤6)对于采集到的振动结构的位移响应数据，其模态坐标表示如下：

其中，

表示振动结构的位移响应数据，

是模态振型矩阵，

表示模态坐标响应，当系统的各阶的模态固有频率不相等是，各阶模态振型相互归一化正交，各阶模态坐标响应相互独立，如下：

其中，

表示两模态坐标响应的期望，Λ"_n×n表示阶数为n的对角矩阵；

步骤7)对于时变结构系统，

是由模态振型向量

组成的，并代表了线性时变结构在L时间段内的统计平均模态振型，可近似估计求得时变结构在(i+(L-1)/2)时刻的瞬间模态振型

同时，

由模态响应函数

构成，并代表了线性时变系统在L时段内的统计平均模态响应，利用单自由度识别技术，可近似估计求得在(i+(L-1)/2)时刻的瞬时模态频率ω_j(i+(L-1)/2)；

步骤8)所识别的模态固有频率是通过快速傅里叶变换得到的，模态固有频率的识别精度取决于频率分辨率△f，当快速傅里叶变换的长度L越长，模态固有频率的识别精度越高，而且频率分辨率△f与采样频率f呈现正比例的关系，频率分辨率△f可定义为如下：

△f＝f/L

步骤9)当线性时变结构的模态参数发生变化时和振动响应信号是非平稳的，如果线性时变结构的模态参数变化非常快，这是一个非平稳程度非常高的时变系统，此时，应该减少滑动窗窗口的长度；定义在一个数据窗内第i阶模态的平均频率变化量为△f_L(i)：

其中，变量f_end(i),f_begin(i),t_end,t_begin分别表示第i阶模态的结束频率，开始频率，整个数据的结束时间和开始时间；

步骤10)滑动窗的窗口的长度L是振动响应数据的采样频率和频率分辨率相除得到的，△f的取值不能取太小，因为不能反映出频率的变化，△f的取值也不能取太大，因为平均频率变化量△f_L(i)不能被识别出来，与此同时，△f_L(i)不能比△f的取值大太多，否则，线性时变结构在一个滑动窗的窗口的长度L内不能被视为时不变的了；

步骤11)模态固有频率的计算，是由单自由度技术计算得到的，利用的是快速傅里叶变换，而快速傅里叶变换的计算复杂度跟滑动窗窗口的长度有直接的关系，根据快速傅里叶变换的原理，滑动窗窗口的长度应该满足如下：

L＝2^α,α＝1,2,…。

实施例1

本实施例中，所述的基于二阶盲辨识的线性时不变结构工作模态参数识别方法，采用的对象是简支梁，简支梁的特性为一段固定，另一端简支，如图1所示，测量简支梁各个测点的时域信号，再通过模态分析方法来识别出简支梁的模态参数。其中，如图2所示，简支梁的尺寸设置为:长670mm,宽56mm,高8mm，材料是45钢，密度为7.85，质量为2.33，泊松比为0.269。

将简支梁进行13等份，并在简支梁上均匀布置12个加速度传感器，在简支梁右边第二个和第三个加速度传感器之间施加由DH40020控制激振器产生的随机激励信号，实验现场布置如图3所示。利用动态采集分析系统的数据采集分析仪采集实验数据，分别收集了采样频率分别为2k,5k和10k的实验数据，采样时间分别为1分钟。实验验证中，利用DHDAS软件可进行实验模态分析识别出简支梁的模态参数作为参考，与SOBI方法识别出的模态参数进行比较。

当采样频率为2k时，利用SOBI方法对实验数据进行模态识别，识别出来了前三阶模态，如图4所示，为SOBI方法分离的成分的快速傅里叶变换图，峰值处的横坐标表示的是模态的固有频率，表1为SOBI方法识别的固有频率的误差，此误差的计算是将识别的固有频率和DHDAS软件识别的固有频率进行比较。

表1基于SOBI识别的固有频率的误差

当简支梁的固有频率被识别出来后，与此同时，模态振型也被识别出来，如图5所示，表2为识别的模态振型的MAC值。

表2模态振型的MAC值

当采样频率为5k时，利用SOBI方法对实验数据进行模态识别，识别出来了前三阶模态，如图6所示，为SOBI方法分离的成分的快速傅里叶变换图，峰值处的横坐标表示的是模态的固有频率，表3为SOBI方法识别的固有频率的误差，此误差的计算是将识别的固有频率和DHDAS软件识别的的固有频率进行比较。

表3基于SOBI识别的固有频率的误差

当简支梁的固有频率被识别出来后，与此同时，模态振型也被识别出来，如图7所示，表4为识别的模态振型的MAC值。

表4模态振型的MAC值

当采样频率为10k时，利用SOBI方法对实验数据进行模态识别，识别出来了前三阶模态，如图8所示，为SOBI方法分离的成分的快速傅里叶变换图，峰值处的横坐标表示的是模态的固有频率，表5为SOBI方法识别的固有频率的误差，此误差的计算是将识别的固有频率和DHDAS软件识别的固有频率进行比较。

表5基于SOBI识别的固有频率的误差

当简支梁的固有频率被识别出来后，与此同时，模态振型也被识别出来，如图9所示，表6为识别的模态振型的MAC值。

表6模态振型的MAC值

实施例2

主要研究具有复杂三维结构的圆柱壳，其中圆柱壳的边界条件是两端简支。在圆柱壳表面布置一定数量的振动传感器，利用比利时LMS公司一套数采前端LMS SCADS-XIII，Test lab 9B系统，DELL M65，记录三个方向的振动响应。其中，圆柱壳的参数设置为：圆柱壳的厚度为0.005米，长度为0.37米，半径为0.1825米，弹性模量为205Gpa，材料的泊松比为0.3，还有材料的密度为7850kg/m³。在仿真中，模态阻尼比η共有两种情况，分别是0.03和0.1。然后，圆柱壳的观测点布置为，绕着圆柱壳的轴向平均分成38圈，其中在每圈中均匀布置115个观测点，那么共有d＝38×115＝4370个观测点，对圆柱壳施加高斯白噪声激励。采样频率为5120Hz，采样时间为1秒，因此T＝5120。利用LMS Virtual.lab有限元方法进行计算，从每个观测点获取到3中不同阻尼比下的X，Y和Z这3个方向的结构位移响应数据，形成3个方向的响应数据集合，如图10所示，是第1118个观测点三个方向的位移响应数据。

对于三种不同阻尼比下的结构位移响应数据，利用SOBI与最小二乘广义逆结合的三维结构工作模态参数识别方法来计算模态振型和模态频率，并且为了评价识别的效果，与将3个方向的直接矩阵组装后利用SOBI方法求解三维结构的工作模态参数，将LMSVitual.lab系统中去阻尼情况下利用有限元方法计算得到的模态振型和固有频率表示真实的模态振型和固有频率，如图11所示，为圆柱壳真实的模态振型。

根据图10所示，在X和Y方向上的位移响应数据的响应量级为10^-12，而Z方向的位移响应数据的响应量级为10^-13，比X和Y方向上的位移响应小，因此，首先利用SOBI方法分解X方向上的位移响应数据，分解得到X方向上的模态振型和结构的模态响应矩阵，利用最小二乘广义伪逆的方法求解Y和Z方向上的模态振型。然后组装三个方向的模态振型，形成识别三维模态振型的算法。最后利用单自由度系统参数识别技术比如快速傅里叶变换(FFT)等，从模态响应矩阵中得到固有频率。

当模态阻尼比η＝0.03时，利用SOBI和最小二乘广义伪逆的三维工作模态参数识别方法所识别出来的固有频率如图12所示。

当模态阻尼比η＝0.03时，表7比较了不同方法下识别的固有频率的识别精度，分别为去阻尼后的有限元方法，所提出的利用SOBI和最小二乘广义伪逆相结合的三维工作模态参数识别方法以及将3个方向的直接矩阵组装后利用SOBI方法求解三维结构的工作模态参数。表8比较了后两种方法下识别的模态振型的MAC值。

表7比较不同方法下的固有频率

表8当阻尼比为0.03时不同方法识别的模态振型的MAC值

根据所提出的SOBI和最小二乘广义伪逆相结合的三维工作模态参数识别方法识别的模态振型，计算每阶模态振型的自MAC值和与其他阶数的模态振型的互MAC值，可得出表9所示的内容。

表9每阶模态振型的自MAC值互MAC值

如表9可知，每阶模态振型的自MAC均为1，互MAC值均接近于0。

当模态阻尼比η＝0.03并加5％的高斯白噪声时，利用SOBI和最小二乘广义伪逆结合的三维工作模态参数识别方法所识别出来的固有频率如图13所示。

当模态阻尼比η＝0.03，加5％高斯白噪声时，表10比较了不同方法下识别的固有频率的识别精度，分别为去阻尼后的有限元方法，所提出的利用SOBI和最小二乘广义伪逆相结合的三维工作模态参数识别方法以及将3个方向的直接矩阵组装后利用SOBI方法求解三维结构的工作模态参数。表11比较了后两种方法下识别的模态振型的MAC值。

表10比较不同方法下的固有频率当加5％高斯白噪声

表11当阻尼比为0.03时不同方法识别的模态振型的MAC值当加5％高斯白噪声

当模态阻尼比η＝0.03，如图14所示为三种方法识别的模态振型的形状，分别为利用有限元求取无阻尼的结构的实际模态振型的形状，提出的方法SOBI和最小二乘广义伪逆相结合的三维工作模态参数识别方法在无噪声情况下的识别的模态振型的形状和在加5％高斯白噪声情况下的识别的模态振型的形状。

当模态阻尼比η＝0.1时，利用SOBI和最小二乘广义伪逆的三维结构工作模态参数识别方法所识别出来的固有频率如图15所示。

当模态阻尼比η＝0.1时，表12比较了不同方法下识别的固有频率的识别精度，分别为去阻尼后的有限元方法，所提出的利用SOBI和最小二乘广义伪逆相结合的三维工作模态参数识别方法以及将3个方向的直接矩阵组装后利用SOBI方法求解三维结构的工作模态参数。表13比较了后两种方法下识别的模态振型的MAC值。

表12比较不同方法下的固有频率

表13当阻尼比为0.1时不同方法识别的模态振型的MAC值

当模态阻尼比η＝0.1并加5％的高斯白噪声时，利用SOBI和最小二乘广义伪逆结合的三维工作模态参数识别方法所识别出来的固有频率如图16所示。

当模态阻尼比η＝0.1，加5％高斯白噪声时，表14比较了不同方法下识别的固有频率的识别精度，分别为去阻尼后的有限元方法，所提出的利用SOBI和最小二乘广义伪逆相结合的三维工作模态参数识别方法以及将3个方向的直接矩阵组装后利用SOBI方法求解三维结构的工作模态参数。表15比较了后两种方法下识别的模态振型的MAC值。

表14比较不同方法下的固有频率当加5％高斯白噪声

表15当阻尼比为0.1时不同方法识别的模态振型的MAC值当加5％高斯白噪声

当模态阻尼比η＝0.1时，如图17所示三种方法识别的模态振型的形状，分别为利用有限元求取无阻尼的结构的实际模态振型的形状，提出的方法SOBI和最小二乘广义伪逆相结合的三维工作模态参数识别方法在无噪声情况下的识别的模态振型的形状和在加5％高斯白噪声情况下的识别的模态振型的形状。

实施例3

本实施例中，所述的基于滑动窗二阶盲辨识的线性时变结构工作模态参数识别方法采用三自由度弹簧振子模拟时变结构，其中，

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，限定记忆长度为L＝1024,傅里叶变换的频率分辨率△f＝0.039Hz，平均频率变化量△f_L(1)＝9.25×10^-4Hz，△f_L(2)＝0.0130Hz以及△f_L(3)＝0.0426Hz。

如图18所示，为基于三自由度时变弹簧振子(图19)测量的位移响应信号以及给时变系统施加的白噪声激励。

如图20所示，为通过理论计算的三阶时变的固有频率变化曲线，由于时变结构中模态振型时刻变化，难以列举出所有振型，基于此，在2000s仿真时间内，随机选取50.025s、675s、1350s、1942.475s(避免随机振动的影响，选取50s时刻之后的数据进行计算)，如图21所示，分别为50.025s、675s、1350s、1942.475s各时刻通过理论计算的各阶模态振型。在t＝50s之前，结构是时不变的系统，因此，可以用SOBI方法识别出此时间段的模态振型和模态固有频率，在图22中，快速傅里叶变换的峰值出即为识别的模态固有频率，而且表16为实际模态固有频率和识别的模态固有频率的相对误差，在图23中比较了实际模态振型和识别的模态振型，它的MAC值在表17中。

表16实际模态固有频率和识别的模态固有频率的相对误差

表17模态振型的MAC值

对于每一个滑动窗的窗口，SOBI方法只运行一次，图24为基于滑动窗SOBI方法识别的模态固有频率随时间的变化。当模态固有频率被识别出来时，模态振型也可被识别出来，图25为在时刻t＝50.025s,t＝675.00s,t＝1350.00s和t＝1942.475s时识别的模态振型。比较实际模态振型和利用滑动窗SOBI方法的时变结构的工作模态参数识别的模态振型，时变的模态振型的MAC值可如图26所示。

对于基于滑动窗ICA方法的时变结构的工作模态参数参数识别，ICA方法选用的目标函数为负熵，选用的优化方法为拟牛顿迭代法。对于基于滑动窗ICA方法的时变结构的工作模态参数参数识别方法可分为两种情况，第一种情况，对于每一个窗口内，ICA方法的迭代过程迭代一次，如图27所示，为时变结构的模态固有频率，在某些窗口内，利用ICA方法没法识别出工作模态参数，因为在计算分离矩阵时陷入了局部最优，未识别的模态参数的百分比如表18所示。

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

利用滑动窗ICA方法识别的时变结构的模态振型，我们选取四个时刻t＝50.025s,t＝675.00s,t＝1350.00s和t＝1942.475s，识别的时变结构的模态振型如图28所示。表19为上述选取的四个时刻用滑动窗ICA方法的线性时变结构工作模态参数识别方法识别的模态振型的MAC值。

表19四个时刻的模态振型的MAC值

对于基于滑动窗ICA方法的时变结构的工作模态参数参数识别方法的另外一种情况，对于每一个窗口内，ICA方法的迭代过程迭代60次，图29为识别的时变结构的模态固有频率，它未识别的模态参数的百分比如表20所示。

表20未识别的模态固有频率是百分比

实施例4

本实施例中，所述的基于滑动窗二阶盲辨识的线性时变结构工作模态参数识别方法采用一维悬臂梁结构模拟时变结构，对于一维悬臂梁结构，在不考虑剪切变形的情况下，通过有限元建模，将一维悬臂梁结构均匀分成40个单元，如图30所示，并且只考虑梁的横向位移和转角，不考虑轴向位移。一维连续体悬臂梁的参数设置为：悬臂梁的长度Len＝1m,0.02m宽and 0.02m高，横截面积为Area＝Wide×High＝4×10^-4m²,惯性矩为I＝Wide×(High)³/12,杨氏模量为E＝2.1×10¹¹N/m²,泊松比为u＝0.3,密度为ρ₀＝7860kg/m³。F表示为对悬臂梁结构施加的高斯白噪声激励。

在有限元方法中，通常连续的悬臂梁被离散为有限的多自由度单元后建立二阶常微分方程组形式的运动控制方程，其中单个单元质量矩阵M^e、刚度矩阵K^e和阻尼矩阵C^e(假设梁的阻尼为比例阻尼)可分别表示为：

C^e＝β_MM^e+β_KK^e

其中L表示的是滑动窗的窗口的长度，β_M和β_K都表示比例系数。然后将单元质量矩阵M^e、刚度矩阵K^e和阻尼矩阵C^e组装成系统的总的质量矩阵M^total，总的刚度矩阵C^total和总的阻尼矩阵K^total，如下所示：

因此，通过有限元方法计算悬臂梁结构的工作模态参数，它的模态固有频率，模态振型和模态阻尼比可表示如下：

其中，M^total _r,K^total _r和C^total _r分别表达第r阶的模态质量矩阵，模态刚度矩阵和模态阻尼矩阵。

在仿真验证中，悬臂梁的密度是随时间而变化的，它的变化情况包括两种速率，一个是0.005，一个是0.08，如下所示：

其中，仿真时间为4s，系统的采样频率为f_s＝10000Hz。

为了避免振动系统在起振阶段受到随机激励的影响，因此在0s到0.5s，将系统设置为恒定不变的，在0.5s之后，系统再实时地发生变化。在实验中，假设系统的初始条件为零，在得到各时刻的集总系统质量矩阵、刚度矩阵和阻尼矩阵后，在悬臂梁的自由端施加白噪声激励，采用Newmark-β求解梁上每个节点的位移(或加速度)响应信号，参数设置如下：Newmark-β积分时间步长1/10000s，阻尼系数β_M＝4×10^-4，β_K＝1×10^-7。

在悬臂梁的自由端施加白噪声激励后得到响应数据是，我们利用基于滑动窗SOBI方法的时变结构的工作模态参数识别方法进行识别，识别的模态参数与利用有限元方法计算得到的模态参数以及真实的模态参数进行对比。

当密度的变化速率为0.005时，设置滑动窗的窗口的长度为L＝2048，此时的快速傅里叶变换的频率分辨率为△f＝4.88Hz，平均频率变化率△f_L(1)＝0.157Hz,△f_L(2)＝1.0932Hz,△f_L(3)＝3.0609Hz。图31表示的是基于滑动窗SOBI方法的时变结构的工作模态参数识别方法识别的模态固有频率和各阶模态振型的MAC值。由于在0s到0.5s，将系统设置为恒定不变的，在0.5s之后，系统再实时地发生变化，因此对于模态振型的MAC值变化是在时间0.625s到3.8977s内变化。

当密度的变化速率为0.005时，设置滑动窗的窗口的长度为L＝4096，此时的快速傅里叶变换的频率分辨率为△f＝2.44Hz，平均频率变化率△f_L(1)＝0.3138Hz,△f_L(2)＝2.1864Hz,△f_L(3)＝6.1218Hz。图32表示的是基于滑动窗SOBI方法的时变结构的工作模态参数识别方法识别的模态固有频率和各阶模态振型的MAC值。由于在0s到0.5s，将系统设置为恒定不变的，在0.5s之后，系统再实时地发生变化，因此对于模态振型的MAC值变化是在时间0.7049s到3.7953s内变化。

当密度的变化速率为0.005时，设置滑动窗的窗口的长度为L＝8192，此时的快速傅里叶变换的频率分辨率为△f＝1.22Hz，平均频率变化率△f_L(1)＝0.6275Hz,△f_L(2)＝4.3727Hz,△f_L(3)＝12.2436Hz。图33表示的是基于滑动窗SOBI方法的时变结构的工作模态参数识别方法识别的模态固有频率和各阶模态振型的MAC值。由于在0s到0.5s，将系统设置为恒定不变的，在0.5s之后，系统再实时地发生变化，因此对于模态振型的MAC值变化是在时间0.9097s到3.5905s内变化。

当密度的变化速率为0.08时，设置滑动窗的窗口的长度为L＝2048，此时的快速傅里叶变换的频率分辨率为△f＝4.88Hz，平均频率变化率△f_L(1)＝0.157Hz,△f_L(2)＝1.0932Hz,△f_L(3)＝3.0609Hz。图34表示的是基于滑动窗SOBI方法的时变结构的工作模态参数识别方法识别的模态固有频率和各阶模态振型的MAC值。由于在0s到0.5s，将系统设置为恒定不变的，在0.5s之后，系统再实时地发生变化，因此对于模态振型的MAC值变化是在时间0.625s到3.8977s内变化。

当密度的变化速率为0.08时，设置滑动窗的窗口的长度为L＝4096，此时的快速傅里叶变换的频率分辨率为△f＝2.44Hz，平均频率变化率△f_L(1)＝0.3138Hz，△f_L(2)＝2.1864Hz，△f_L(3)＝6.1218Hz。图35表示的是基于滑动窗SOBI方法的时变结构的工作模态参数识别方法识别的模态固有频率和各阶模态振型的MAC值。由于在0s到0.5s，将系统设置为恒定不变的，在0.5s之后，系统再实时地发生变化，因此对于模态振型的MAC值变化是在时间0.7049s到3.7953s内变化。

当密度的变化速率为0.08时，设置滑动窗的窗口的长度为L＝8192，此时的快速傅里叶变换的频率分辨率为△f＝1.22Hz，平均频率变化率△f_L(1)＝0.6275Hz,△f_L(2)＝4.3727Hz,△f_L(3)＝12.2436Hz。图36表示的是基于滑动窗SOBI方法的时变结构的工作模态参数识别方法识别的模态固有频率和各阶模态振型的MAC值。由于在0s到0.5s，将系统设置为恒定不变的，在0.5s之后，系统再实时地发生变化，因此对于模态振型的MAC值变化是在时间0.9097s到3.5905s内变化。

上述实施例仅是用来说明本发明，而并非用作对本发明的限定。只要是依据本发明的技术实质，对上述实施例进行变化、变型等都将落在本发明的权利要求的范围内。

Claims (2)

1.一种结合二阶盲辨识和最小二乘广义逆方法的三维结构的工作模态参数识别方法，其特征在于：为了求解三维连续体复杂结构的工作模态参数识别，先从其中的一个方向入手，由于每个方向的振动响应都不同，以二阶盲辨识算法为基础，选取响应最大的那个方向用二阶盲辨识进行工作模态参数识别，此时的做法跟一维结构的工作模态参数识别是一致的；当分解出响应最大的那个方向的振动响应，可得到模态响应矩阵和一维模态振型，然后将模态响应矩阵使用最小二乘广义伪逆，回代到另外两个方向的位移响应信号中，然后组装各方向的模态振型，形成识别三维模态振型的算法；最后利用单自由度系统参数识别技术快速傅里叶变换，从模态响应矩阵中得到包括模态固有频率的信息；

具体步骤如下：

步骤1)首先，利用二阶盲辨识方法分解X方向上的振动响应，分解得到X方向上的模态振型Ψ_d×l和结构的模态响应矩阵H_l×T(t)：

(X_three)_d×T(t)≈Ψ_d×lH_l×T(t)

其中(X_three)_d×T(t)表示X方向上的振动响应信号，d为三维结构的阶数，T为时间；

步骤2)由于三个方向的模态响应矩阵是相同的，利用基于二阶盲辨识的三维模态参数识别算法对X方向上进行二阶盲辨识分解之后，需要回代到另外两个方向，模态响应是个矩阵，实际中，是对另外两个方向的响应乘以模态响应矩阵的右伪逆来实现，最小二乘广义逆是在误差平方和最小意义下的一致无偏最优估计，因此利用最小二乘广义伪逆的方法求解Y和Z方向上的模态振型：

其中(Y_three)_d×T(t)和(Z_three)_d×T(t)分别表示Y方向上和Z方向上的振动响应信号，d为三维结构的阶数，T为时间；O_d×l和B_d×l分别表示Y和Z方向上的模态振型矩阵；

步骤3)然后组装各方向的模态振型，形成识别三维模态振型的算法；最后利用单自由度系统参数识别技术快速傅里叶变换，可从模态响应矩阵中得到包括模态固有频率的信息。

2.一种用于实现如权利要求1所述的工作模态参数识别方法的时不变三维圆柱壳工作模态参数识别的实验装置，包括一个一端固支一端自由的三维薄壁金属圆柱壳、n个三向加速度传感器、数据采集卡、激振器或力锤及电脑终端；所述金属圆柱壳表面上均匀放置n个三向加速度传感器，用于采集使用所述力锤或激振器垂直于圆柱壳表面施加激励产生的三个方向的振动响应信号，所述数据采集卡与所述n个三向加速度传感器相连用于接收所述振动响应信号并发送给所述电脑终端存储，所述电脑终端包括一模态参数识别模块，所述模态参数识别模块基于电脑终端存储的振动响应信号使用结合二阶盲辨识和最小二乘广义逆方法的三维结构的工作模态参数识别方法识别三维圆柱壳的模态参数，并与由有限元分析或理论解相比较，验证结合二阶盲辨识和最小二乘广义逆方法的三维工作模态参数识别方法的正确性。

Images