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

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

Info

Publication number
CN107357977A
CN107357977A CN201710500228.0A CN201710500228A CN107357977A CN 107357977 A CN107357977 A CN 107357977A CN 201710500228 A CN201710500228 A CN 201710500228A CN 107357977 A CN107357977 A CN 107357977A
Authority
CN
China
Prior art keywords
mrow
modal
time
msub
matrix
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201710500228.0A
Other languages
English (en)
Other versions
CN107357977B (zh
Inventor
王成
王建英
张天舒
赖雄鸣
张忆文
陈叶旺
张惠臻
何霆
皮慧娟
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Huaqiao University
Original Assignee
Huaqiao University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Huaqiao University filed Critical Huaqiao University
Priority to CN201710500228.0A priority Critical patent/CN107357977B/zh
Publication of CN107357977A publication Critical patent/CN107357977A/zh
Application granted granted Critical
Publication of CN107357977B publication Critical patent/CN107357977B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/13—Architectural design, e.g. computer-aided architectural design [CAAD] related to design of buildings, bridges, landscapes, production plants or roads
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Geometry (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • General Physics & Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Architecture (AREA)
  • Civil Engineering (AREA)
  • Structural Engineering (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Complex Calculations (AREA)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)

Abstract

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

Description

基于二阶盲辨识的线性结构工作模态参数识别方法及装置
技术领域
本发明涉及一种基于二阶盲辨识的线性时不变结构工作模态参数识别方法，在一种时不变的工作模态参数识别的装置产生的数据上进行方法验证。一种结合二阶盲辨识和最小二乘广义逆方法的三维结构的工作模态参数识别方法，在一种时不变三维圆柱壳工作模态参数识别的实验装置产生的数据上进行方法验证。一种基于滑动窗二阶盲辨识的线性时变结构工作模态参数识别方法，以及该方法在仿真数据上进行方法验证。
背景技术
工作模态分析是仅利用结构在工作时而不是在实验中利用传感器获得的振动响应数据中提取出模态参数，其中模态参数包括模态阻尼比，模态振型和模态固有频率。目前工作模态分析方法在解决工程问题上是非常重要的分析方法，特别是在航空，土木工程等领域。
然而，对于大型结构的模态参数不仅仅受到外激励和内激励的影响，而且在工作状态中很难保持不变，也就是说，结构的模态参数是随着时间的变化而变化的，不是时不变的，而是时变的。因此，要识别时变结构的模态参数，如果仅仅依靠时不变的方法是不可能实现的，还得借助时变的方法。国内外学者对时变结构的工作模态参数识别也做了大量的研究，比如有提出TARMA模型，利用小波变换和希尔伯特变换来识别时变结构的模态参数。
发明内容
本发明的目的在于克服现有技术的不足，提供一种基于滑动窗二阶盲辨识的线性时变结构工作模态参数识别方法，还提供基于滑动窗二阶盲辨识的线性时变结构工作模态参数识别方法在设备故障诊断与健康状态监测上的应用。
本发明的技术方案如下：
一种基于二阶盲辨识的线性时不变一维结构工作模态参数识别方法，利用二阶盲辨识方法分解线性时不变一维结构的振动响应信号，从而得到时不变结构的各阶模态参数，其中在二阶盲辨识的解混旋转阶段，利用多个时延非零的相关矩阵同时对角化的方法，消除了由于单个时延选择不适而造成的算法结果的不理想，从而有效识别出时不变一维结构工作模态参数。
所述的基于二阶盲辨识的线性时不变一维结构工作模态参数识别方法，具体步骤如下：
步骤1)时不变结构的模态参数随时间的变化而变化，根据结构动力学理论，对于自由度线性时变振动结构系统，它在物理坐标系统中的运动方程为：
其中，M，C和分别表示质量矩阵，阻尼矩阵和刚度矩阵，与此同时，它们受结构的影响而随时间发生变化；表示外载荷的激励向量，和分别表示加速度响应信号，速度响应信号和位移响应信号；表示维度为n×T的矩阵；
步骤2)二阶盲辨识是另外一种盲源分离算法，它利用了信号的另外一种特性—时序结构，即时延协方差信息，要求信号具有不同的谱特性或者不同的相关函数，它的处理对象是时间信号，当然也可以处理统计独立的源信号，为了正确地求解出分离矩阵和源信号，二阶盲辨识算法做出了如下的假设：
假设一：混合矩阵列满秩；
假设二：源信号相互不相关且具有不同的自相关函数；
假设三：源信号是平稳信号；
步骤3)根据上述假设，源信号的协方差矩阵Rs(0)满足如下：
Rs(0)＝E[S(t)ST(t)]＝I
其中，E[*]表示期望，S(t)表示源信号，I表示单位矩阵。
步骤4)，观测信号的协方差矩阵可表示如下：
其中，Xbss(t)表示观测信号，A表示混合矩阵。
步骤5)对观测信号Xbss(t)进行白化预处理，即进行线性变换：
Z(t)＝VXbss(t)
其中，是一个白化矩阵，线性变换的目的是使得协方差矩阵是一个单位矩阵，去除信号之间的相关性；
步骤6)源信号的延时协方差矩阵Rs(τ)类推定义如下：
Rs(τ)＝E[S(t+τ)ST(t)]
步骤7)与此同时，源信号白化后的延时协方差矩阵可得到：
Rz(τ)＝E[Z(t+τ)ZT(t)]
步骤8)选择一组不同的值τ1,τ2,…,τp，可得到一系列延时协方差矩阵Rs(τi)，其中i＝1,…,p，可以得出如下：
Z(t)＝VXbss(t)＝(VA)S(t)＝US(t)
步骤9)其中，U＝VA，因此，可推导出如下：
Rz(τi)＝URs(τi)UT,i＝1,2,…,p
步骤10)由于源信号被假设成不相关的，而且源信号白化后的矩阵Z(t)中的各向量相互正交归一化，因此，矩阵U也是正交归一化的，而且
Rs(τi)＝UTRz(τi)U,i＝1,2,…,p
步骤11)可以看出，Rs(τi)是一个对角化矩阵，可以看出矩阵U是归一化正交的，因此，利用Rz(τi)和最优逼近算法来寻找最优的矩阵U，可计算出分离矩阵W如下：
W＝UTV
步骤12)同时，源信号S(t)也可计算得到。
一种时不变的工作模态参数识别的装置，包括一根一端简支一端固支的金属梁、n个加速度传感器、数据采集卡、激振器或力锤及电脑终端；所述金属梁设置有n个等分点，其中每等分点上放置有一个加速度传感器，用于采集使用所述激振器或力锤对所述金属梁的某一点施加激励所产生的振动响应信号，所述数据采集卡与所述n个加速度传感器相连用于接收所述振动响应信号并发送给所述电脑终端存储，所述电脑终端包括一模态参数识别模块，所述模态参数识别模块基于电脑终端存储的振动响应信号使用基于二阶盲辨识的线性时不变一维结构工作模态参数识别方法识别出金属梁的模态参数，并与梁的理论模态参数作对比以验证基于二阶盲辨识的线性时不变结构工作模态参数识别方法的正确性。
在实际工程应用中，一般的工程结构都是三维结构，因此，从一维结构到三维结构的模态参数分析，是目前的一大挑战，本发明一种结合二阶盲辨识和最小二乘广义逆方法的三维结构的工作模态参数识别方法，为了求解三维连续体复杂结构的工作模态参数识别，先从其中的一个方向入手，由于每个方向的振动响应都不同，以二阶盲辨识算法为基础，选取响应最大的那个方向用二阶盲辨识进行工作模态参数识别，此时的做法跟一维结构的工作模态参数识别是一致的；当分解出响应最大的那个方向的振动响应，可得到模态响应矩阵和一维模态振型，然后将模态响应矩阵使用最小二乘广义伪逆，回代到另外两个方向的位移响应信号中，然后组装各方向的模态振型，形成识别三维模态振型的算法；最后利用单自由度系统参数识别技术快速傅里叶变换，从模态响应矩阵中得到包括模态固有频率的信息。
所述的结合二阶盲辨识和最小二乘广义逆方法的三维结构的工作模态参数识别方法，具体步骤如下：
步骤1)首先，利用二阶盲辨识方法分解X方向上的振动响应，分解得到X方向上的模态振型Ψd×l和结构的模态响应矩阵Hl×T(t)：
(Xthree)d×T(t)≈Ψd×lHl×T(t)
其中(Xthree)d×T(t)表示X方向上的振动响应信号，d为三维结构的阶数，T为时间。
步骤2)由于三个方向的模态响应矩阵是相同的，利用基于二阶盲辨识的三维模态参数识别算法对X方向上进行二阶盲辨识分解之后，需要回代到另外两个方向，模态响应是个矩阵，实际中，是对另外两个方向的响应乘以模态响应矩阵的右伪逆来实现，最小二乘广义逆是在误差平方和最小意义下的一致无偏最优估计，因此利用最小二乘广义伪逆的方法求解Y 和Z方向上的模态振型：
其中(Ythree)d×T(t)和(Zthree)d×T(t)分别表示Y方向上和Z方向上的振动响应信号，d为三维结构的阶数，T为时间。Od×l和Bd×l分别表示Y和Z方向上的模态振型矩阵。
步骤3)然后组装各方向的模态振型，形成识别三维模态振型的算法；最后利用单自由度系统参数识别技术快速傅里叶变换，可从模态响应矩阵中得到包括模态固有频率的信息。
一种时不变三维圆柱壳工作模态参数识别的实验装置，包括一个一端固支一端自由的三维薄壁金属圆柱壳、n个三向加速度传感器、数据采集卡、激振器或力锤及电脑终端；所述金属圆柱壳表面上均匀放置n个三向加速度传感器，用于采集使用所述力锤或激振器垂直于圆柱壳表面施加激励产生的三个方向的振动响应信号，所述数据采集卡与所述n个三向加速度传感器相连用于接收所述振动响应信号并发送给所述电脑终端存储，所述电脑终端包括一模态参数识别模块，所述模态参数识别模块基于电脑终端存储的振动响应信号使用结合二阶盲辨识和最小二乘广义逆方法的三维结构的工作模态参数识别方法识别三维圆柱壳的模态参数，并与由有限元分析或理论解相比较，验证结合二阶盲辨识和最小二乘广义逆方法的三维工作模态参数识别方法的正确性。
一种基于滑动窗二阶盲辨识的线性时变一维结构工作模态参数识别方法，其分析的结构是线性时变结构，模态参数是随时间的改变而发生变化的，因此主要结合短时时不变理论与二阶盲辨识算法，利用二阶盲辨识算法在各窗内的统计特性，估计出各时刻的工作模态参数，然后各时刻求得的工作模态参数连接起来，从而实现时变线性结构工作模态参数识别；所述工作模态参数包括各阶模态的固有频率和模态振型。
所述的基于滑动窗二阶盲辨识的线性时变一维结构工作模态参数识别方法，具体步骤如下：
步骤1)时变结构的模态参数随时间的变化而变化，根据结构动力学理论，对于自由度线性时变振动结构系统，它在物理坐标系统中的运动方程为:
其中，M(t)，C(t)和分别表示随时间变化的质量矩阵，阻尼矩阵和刚度矩阵，与此同时，它们受结构的影响而随时间发生变化；表示外载荷的激励向量，和分别表示加速度响应信号，速度响应信号和位移响应信号；
步骤2)根据“短时时不变”理论，时变离散多自由度系统在时间τ∈[tbegin,tend]内，它的质量，阻尼和刚度看作是时不变的，因此，在物理坐标系中的动力学方程可表示为：
其中，K＝end，表示时变系统中最后时刻，S'(τ)表示为当t＝τ时刻的时不变结构，S'表示一组有限多个线性时不变结构组成时变结构的集合；
步骤3)对于时不变小阻尼结构，响应数据可以被分为有限多个部分，在第τ个部分，选取一定的窗口的长度，线性系统的模态坐标响应为：
其中，Φ(τ)和q(τ,t)分别表示当第τ个窗口的模态振型矩阵和模态响应向量；
步骤4)当时不变结构的每一阶模态固有频率ωi都不相等时，各阶模态振型之间满足归一化正交，各阶模态响应相互不相关，如下：
其中，E(q(τ,t)qT(τ,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阶模态的平均频率变化量为△fL(i)：
其中，变量fend(i),fbegin(i),tend,tbegin分别表示第i阶模态的结束频率，开始频率，整个数据的结束时间和开始时间；
步骤10)滑动窗的窗口的长度L是振动响应数据的采样频率和频率分辨率相除得到的，△f的取值不能取太小，因为不能反映出频率的变化，△f的取值也不能取太大，因为平均频率变化量△fL(i)不能被识别出来，与此同时，△fL(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)根据上述假设，源信号的协方差矩阵Rs(0)满足如下：
Rs(0)＝E[S(t)ST(t)]＝I
其中，E[*]表示期望，S(t)表示源信号，I表示单位矩阵。
步骤4)，观测信号的协方差矩阵可表示如下：
其中，Xbss(t)表示观测信号，A表示混合矩阵。
步骤5)对观测信号Xbss(t)进行白化预处理，即进行线性变换：
Z(t)＝VXbss(t)
其中，是一个白化矩阵，线性变换的目的是使得协方差矩阵是一个单位矩阵，去除信号之间的相关性；
步骤6)源信号的延时协方差矩阵Rs(τ)类推定义如下：
Rs(τ)＝E[S(t+τ)ST(t)]
步骤7)与此同时，源信号白化后的延时协方差矩阵可得到：
Rz(τ)＝E[Z(t+τ)ZT(t)]
步骤8)选择一组不同的值τ1,τ2,…,τp，可得到一系列延时协方差矩阵Rs(τi)，其中i＝1,…,p，可以得出如下：
Z(t)＝VXbss(t)＝(VA)S(t)＝US(t)
步骤9)其中，U＝VA，因此，可推导出如下：
Rz(τi)＝URs(τi)UT,i＝1,2,…,p
步骤10)由于源信号被假设成不相关的，而且源信号白化后的矩阵Z(t)中的各向量相互正交归一化，因此，矩阵U也是正交归一化的，而且
Rs(τi)＝UTRz(τi)U,i＝1,2,…,p
步骤11)可以看出，Rs(τi)是一个对角化矩阵，可以看出矩阵U是归一化正交的，因此，利用Rz(τi)和最优逼近算法来寻找最优的矩阵U，可计算出分离矩阵W如下：
W＝UTV
步骤12)同时，源信号S(t)也可计算得到。
一种结合二阶盲辨识和最小二乘广义逆方法的三维结构的工作模态参数识别方法，以二阶盲辨识算法为基础，选取响应最大的那个方向用二阶盲辨识进行工作模态参数识别，可得到模态响应矩阵和一维模态振型，然后将模态响应矩阵使用最小二乘广义伪逆，回代到另外两个方向的位移响应信号中，然后组装各方向的模态振型，形成识别三维模态振型的算法。最后利用单自由度系统参数识别技术比如快速傅里叶变换等，可从模态响应矩阵中得到模态固有频率等信息，具体包括：
步骤1)首先，利用二阶盲辨识方法分解X方向上的振动响应，分解得到X方向上的模态振型Ψd×l和结构的模态响应矩阵Hl×T(t)：
(Xthree)d×T(t)≈Ψd×lHl×T(t)
其中(Xthree)d×T(t)表示X方向上的振动响应信号，d为三维结构的阶数，T为时间。
步骤2)由于三个方向的模态响应矩阵是相同的，利用基于二阶盲辨识的三维模态参数识别算法对X方向上进行二阶盲辨识分解之后，需要回代到另外两个方向，模态响应是个矩阵，实际中，是对另外两个方向的响应乘以模态响应矩阵的右伪逆来实现，最小二乘广义逆是在误差平方和最小意义下的一致无偏最优估计，因此利用最小二乘广义伪逆的方法求解Y 和Z方向上的模态振型：
其中(Ythree)d×T(t)和(Zthree)d×T(t)分别表示Y方向上和Z方向上的振动响应信号，d为三维结构的阶数，T为时间。Od×l和Bd×l分别表示Y和Z方向上的模态振型矩阵。
步骤3)然后组装各方向的模态振型，形成识别三维模态振型的算法；最后利用单自由度系统参数识别技术快速傅里叶变换，可从模态响应矩阵中得到包括模态固有频率的信息。
一种基于滑动窗二阶盲辨识的线性时变结构工作模态参数识别方法，结合短时时不变理论与二阶盲辨识算法，利用二阶盲辨识算法在各窗内的统计特性，估计出各时刻的工作模态参数(包括各阶模态的固有频率和模态振型)，然后各时刻求得的工作模态参数连接起来，从而实现时变线性结构工作模态参数识别。
具体步骤如下：
步骤1)时变结构的模态参数随时间的变化而变化，根据结构动力学理论，对于自由度线性时变振动结构系统，它在物理坐标系统中的运动方程为:
其中，M(t)，C(t)和分别表示随时间变化的质量矩阵，阻尼矩阵和刚度矩阵，与此同时，它们受结构的影响而随时间发生变化；表示外载荷的激励向量，和分别表示加速度响应信号，速度响应信号和位移响应信号；
步骤2)根据“短时时不变”理论，时变离散多自由度系统在时间τ∈[tbegin,tend]内，它的质量，阻尼和刚度看作是时不变的，因此，在物理坐标系中的动力学方程可表示为：
其中，K＝end，表示时变系统中最后时刻，S'(τ)表示为当t＝τ时刻的时不变结构，S'表示一组有限多个线性时不变结构组成时变结构的集合；
步骤3)对于时不变小阻尼结构，响应数据可以被分为有限多个部分，在第τ个部分，选取一定的窗口的长度，线性系统的模态坐标响应为：
其中，Φ(τ)和q(τ,t)分别表示当第τ个窗口的模态振型矩阵和模态响应向量；
步骤4)当时不变结构的每一阶模态固有频率ωi都不相等时，各阶模态振型之间满足归一化正交，各阶模态响应相互不相关，如下：
其中，E(q(τ,t)qT(τ,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阶模态的平均频率变化量为△fL(i)：
其中，变量fend(i),fbegin(i),tend,tbegin分别表示第i阶模态的结束频率，开始频率，整个数据的结束时间和开始时间；
步骤10)滑动窗的窗口的长度L是振动响应数据的采样频率和频率分辨率相除得到的，△f的取值不能取太小，因为不能反映出频率的变化，△f的取值也不能取太大，因为平均频率变化量△fL(i)不能被识别出来，与此同时，△fL(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/m3。在仿真中，模态阻尼比η共有两种情况，分别是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
本实施例中，所述的基于滑动窗二阶盲辨识的线性时变结构工作模态参数识别方法采用三自由度弹簧振子模拟时变结构，其中，m2＝1kg，m3＝1kg； k1＝1000N/m，k2＝1000N/m，k3＝1000N/m；c1＝0.01N.s/m，c2＝0.01N.s/m， c3＝0.01N.s/m。初始位移及速度为0。物块1受到均值为0，方差为1的高斯白噪声激励F1。基于Matlab/Simulink仿真，采样间隔为0.025s，采样频率为40Hz，仿真时间为2000s，限定记忆长度为L＝1024,傅里叶变换的频率分辨率△f＝0.039Hz，平均频率变化量△fL(1)＝9.25×10-4Hz，△fL(2)＝0.0130Hz以及△fL(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-4m2,惯性矩为 I＝Wide×(High)3/12,杨氏模量为E＝2.1×1011N/m2,泊松比为u＝0.3,密度为ρ0＝7860kg/m3。F表示为对悬臂梁结构施加的高斯白噪声激励。
在有限元方法中，通常连续的悬臂梁被离散为有限的多自由度单元后建立二阶常微分方程组形式的运动控制方程，其中单个单元质量矩阵Me、刚度矩阵Ke和阻尼矩阵Ce(假设梁的阻尼为比例阻尼)可分别表示为：
Ce＝βMMe+βKKe
其中L表示的是滑动窗的窗口的长度，βM和βK都表示比例系数。然后将单元质量矩阵Me、刚度矩阵Ke和阻尼矩阵Ce组装成系统的总的质量矩阵Mtotal，总的刚度矩阵Ctotal和总的阻尼矩阵Ktotal，如下所示：
因此，通过有限元方法计算悬臂梁结构的工作模态参数，它的模态固有频率，模态振型和模态阻尼比可表示如下：
其中，Mtotal r,Ktotal r和Ctotal r分别表达第r阶的模态质量矩阵，模态刚度矩阵和模态阻尼矩阵。
在仿真验证中，悬臂梁的密度是随时间而变化的，它的变化情况包括两种速率，一个是 0.005，一个是0.08，如下所示：
其中，仿真时间为4s，系统的采样频率为fs＝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，平均频率变化率△fL(1)＝0.157Hz,△fL(2)＝1.0932Hz, △fL(3)＝3.0609Hz。图31表示的是基于滑动窗SOBI方法的时变结构的工作模态参数识别方法识别的模态固有频率和各阶模态振型的MAC值。由于在0s到0.5s，将系统设置为恒定不变的，在0.5s之后，系统再实时地发生变化，因此对于模态振型的MAC值变化是在时间0.625s 到3.8977s内变化。
当密度的变化速率为0.005时，设置滑动窗的窗口的长度为L＝4096，此时的快速傅里叶变换的频率分辨率为△f＝2.44Hz，平均频率变化率△fL(1)＝0.3138Hz,△fL(2)＝2.1864Hz, △fL(3)＝6.1218Hz。图32表示的是基于滑动窗SOBI方法的时变结构的工作模态参数识别方法识别的模态固有频率和各阶模态振型的MAC值。由于在0s到0.5s，将系统设置为恒定不变的，在0.5s之后，系统再实时地发生变化，因此对于模态振型的MAC值变化是在时间 0.7049s到3.7953s内变化。
当密度的变化速率为0.005时，设置滑动窗的窗口的长度为L＝8192，此时的快速傅里叶变换的频率分辨率为△f＝1.22Hz，平均频率变化率△fL(1)＝0.6275Hz,△fL(2)＝4.3727Hz, △fL(3)＝12.2436Hz。图33表示的是基于滑动窗SOBI方法的时变结构的工作模态参数识别方法识别的模态固有频率和各阶模态振型的MAC值。由于在0s到0.5s，将系统设置为恒定不变的，在0.5s之后，系统再实时地发生变化，因此对于模态振型的MAC值变化是在时间 0.9097s到3.5905s内变化。
当密度的变化速率为0.08时，设置滑动窗的窗口的长度为L＝2048，此时的快速傅里叶变换的频率分辨率为△f＝4.88Hz，平均频率变化率△fL(1)＝0.157Hz,△fL(2)＝1.0932Hz, △fL(3)＝3.0609Hz。图34表示的是基于滑动窗SOBI方法的时变结构的工作模态参数识别方法识别的模态固有频率和各阶模态振型的MAC值。由于在0s到0.5s，将系统设置为恒定不变的，在0.5s之后，系统再实时地发生变化，因此对于模态振型的MAC值变化是在时间0.625s 到3.8977s内变化。
当密度的变化速率为0.08时，设置滑动窗的窗口的长度为L＝4096，此时的快速傅里叶变换的频率分辨率为△f＝2.44Hz，平均频率变化率△fL(1)＝0.3138Hz，△fL(2)＝2.1864Hz，△fL(3)＝6.1218Hz。图35表示的是基于滑动窗SOBI方法的时变结构的工作模态参数识别方法识别的模态固有频率和各阶模态振型的MAC值。由于在0s到0.5s，将系统设置为恒定不变的，在0.5s之后，系统再实时地发生变化，因此对于模态振型的MAC值变化是在时间 0.7049s到3.7953s内变化。
当密度的变化速率为0.08时，设置滑动窗的窗口的长度为L＝8192，此时的快速傅里叶变换的频率分辨率为△f＝1.22Hz，平均频率变化率△fL(1)＝0.6275Hz,△fL(2)＝4.3727Hz, △fL(3)＝12.2436Hz。图36表示的是基于滑动窗SOBI方法的时变结构的工作模态参数识别方法识别的模态固有频率和各阶模态振型的MAC值。由于在0s到0.5s，将系统设置为恒定不变的，在0.5s之后，系统再实时地发生变化，因此对于模态振型的MAC值变化是在时间 0.9097s到3.5905s内变化。
上述实施例仅是用来说明本发明，而并非用作对本发明的限定。只要是依据本发明的技术实质，对上述实施例进行变化、变型等都将落在本发明的权利要求的范围内。

Claims (8)

1.一种基于二阶盲辨识的线性时不变一维结构工作模态参数识别方法，其特征在于：利用二阶盲辨识方法分解线性时不变一维结构的振动响应信号，从而得到时不变结构的各阶模态参数，其中在二阶盲辨识的解混旋转阶段，利用多个时延非零的相关矩阵同时对角化的方法，消除了由于单个时延选择不适而造成的算法结果的不理想，从而有效识别出时不变一维结构工作模态参数。
2.根据权利要求1所述的基于二阶盲辨识的线性时不变一维结构工作模态参数识别方法，其特征在于，具体步骤如下：
步骤1)时不变结构的模态参数随时间的变化而变化，根据结构动力学理论，对于自由度线性时变振动结构系统，它在物理坐标系统中的运动方程为：
<mrow> <mi>S</mi> <mo>:</mo> <mi>M</mi> <mover> <mi>X</mi> <mo>&amp;CenterDot;&amp;CenterDot;</mo> </mover> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>+</mo> <mi>C</mi> <mover> <mi>X</mi> <mo>&amp;CenterDot;</mo> </mover> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>+</mo> <mi>K</mi> <mi>X</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>F</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> </mrow>
其中，M，C和分别表示质量矩阵，阻尼矩阵和刚度矩阵，与此同时，它们受结构的影响而随时间发生变化；表示外载荷的激励向量，和分别表示加速度响应信号，速度响应信号和位移响应信号；表示维度为n×T的矩阵；
步骤2)二阶盲辨识是另外一种盲源分离算法，它利用了信号的另外一种特性—时序结构，即时延协方差信息，要求信号具有不同的谱特性或者不同的相关函数，它的处理对象是时间信号，当然也可以处理统计独立的源信号，为了正确地求解出分离矩阵和源信号，二阶盲辨识算法做出了如下的假设：
假设一：混合矩阵列满秩；
假设二：源信号相互不相关且具有不同的自相关函数；
假设三：源信号是平稳信号；
步骤3)根据上述假设，源信号的协方差矩阵Rs(0)满足如下：
Rs(0)＝E[S(t)ST(t)]＝I
其中，E[*]表示期望，S(t)表示源信号，I表示单位矩阵。
步骤4)，观测信号的协方差矩阵可表示如下：
<mrow> <msub> <mi>R</mi> <mi>x</mi> </msub> <mrow> <mo>(</mo> <mn>0</mn> <mo>)</mo> </mrow> <mo>=</mo> <mi>E</mi> <mo>&amp;lsqb;</mo> <msub> <mi>X</mi> <mrow> <mi>b</mi> <mi>s</mi> <mi>s</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <msubsup> <mi>X</mi> <mrow> <mi>b</mi> <mi>s</mi> <mi>s</mi> </mrow> <mi>T</mi> </msubsup> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> <mo>=</mo> <msub> <mi>AR</mi> <mi>s</mi> </msub> <mrow> <mo>(</mo> <mn>0</mn> <mo>)</mo> </mrow> <msup> <mi>A</mi> <mi>T</mi> </msup> </mrow>
其中，Xbss(t)表示观测信号，A表示混合矩阵。
步骤5)对观测信号Xbss(t)进行白化预处理，即进行线性变换：
Z(t)＝VXbss(t)
其中，是一个白化矩阵，线性变换的目的是使得协方差矩阵是一个单位矩阵，去除信号之间的相关性；
步骤6)源信号的延时协方差矩阵Rs(τ)类推定义如下：
Rs(τ)＝E[S(t+τ)ST(t)]
步骤7)与此同时，源信号白化后的延时协方差矩阵可得到：
Rz(τ)＝E[Z(t+τ)ZT(t)]
步骤8)选择一组不同的值τ1,τ2,…,τp，可得到一系列延时协方差矩阵Rs(τi)，其中i＝1,…,p，可以得出如下：
Z(t)＝VXbss(t)＝(VA)S(t)＝US(t)
步骤9)其中，U＝VA，因此，可推导出如下：
Rz(τi)＝URs(τi)UT,i＝1,2,…,p
步骤10)由于源信号被假设成不相关的，而且源信号白化后的矩阵Z(t)中的各向量相互正交归一化，因此，矩阵U也是正交归一化的，而且
Rs(τi)＝UTRz(τi)U,i＝1,2,…,p
步骤11)可以看出，Rs(τi)是一个对角化矩阵，可以看出矩阵U是归一化正交的，因此，利用Rz(τi)和最优逼近算法来寻找最优的矩阵U，可计算出分离矩阵W如下：
W＝UTV
步骤12)同时，源信号S(t)也可计算得到。
3.一种时不变的工作模态参数识别的装置，其特征在于，包括一根一端简支一端固支的金属梁、n个加速度传感器、数据采集卡、激振器或力锤及电脑终端；所述金属梁设置有n个等分点，其中每等分点上放置有一个加速度传感器，用于采集使用所述激振器或力锤对所述金属梁的某一点施加激励所产生的振动响应信号，所述数据采集卡与所述n个加速度传感器相连用于接收所述振动响应信号并发送给所述电脑终端存储，所述电脑终端包括一模态参数识别模块，所述模态参数识别模块基于电脑终端存储的振动响应信号使用基于二阶盲辨识的线性时不变一维结构工作模态参数识别方法识别出金属梁的模态参数，并与梁的理论模态参数作对比以验证基于二阶盲辨识的线性时不变结构工作模态参数识别方法的正确性。
4.一种结合二阶盲辨识和最小二乘广义逆方法的三维结构的工作模态参数识别方法，其特征在于：在实际工程应用中，一般的工程结构都是三维结构，因此，从一维结构到三维结构的模态参数分析，是目前的一大挑战，为了求解三维连续体复杂结构的工作模态参数识别，先从其中的一个方向入手，由于每个方向的振动响应都不同，以二阶盲辨识算法为基础，选取响应最大的那个方向用二阶盲辨识进行工作模态参数识别，此时的做法跟一维结构的工作模态参数识别是一致的；当分解出响应最大的那个方向的振动响应，可得到模态响应矩阵和一维模态振型，然后将模态响应矩阵使用最小二乘广义伪逆，回代到另外两个方向的位移响应信号中，然后组装各方向的模态振型，形成识别三维模态振型的算法；最后利用单自由度系统参数识别技术快速傅里叶变换，从模态响应矩阵中得到包括模态固有频率的信息。
5.根据权利要求4所述的结合二阶盲辨识和最小二乘广义逆方法的三维结构的工作模态参数识别方法，其特征在于，具体步骤如下：
步骤1)首先，利用二阶盲辨识方法分解X方向上的振动响应，分解得到X方向上的模态振型Ψd×l和结构的模态响应矩阵Hl×T(t)：
(Xthree)d×T(t)≈Ψd×lHl×T(t)
其中(Xthree)d×T(t)表示X方向上的振动响应信号，d为三维结构的阶数，T为时间。
步骤2)由于三个方向的模态响应矩阵是相同的，利用基于二阶盲辨识的三维模态参数识别算法对X方向上进行二阶盲辨识分解之后，需要回代到另外两个方向，模态响应是个矩阵，实际中，是对另外两个方向的响应乘以模态响应矩阵的右伪逆来实现，最小二乘广义逆是在误差平方和最小意义下的一致无偏最优估计，因此利用最小二乘广义伪逆的方法求解Y和Z方向上的模态振型：
<mrow> <msub> <mi>O</mi> <mrow> <mi>d</mi> <mo>&amp;times;</mo> <mi>l</mi> </mrow> </msub> <mo>=</mo> <msub> <mrow> <mo>(</mo> <msub> <mi>Y</mi> <mrow> <mi>t</mi> <mi>h</mi> <mi>r</mi> <mi>e</mi> <mi>e</mi> </mrow> </msub> <mo>)</mo> </mrow> <mrow> <mi>d</mi> <mo>&amp;times;</mo> <mi>T</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <msubsup> <mi>H</mi> <mrow> <mi>T</mi> <mo>&amp;times;</mo> <mi>l</mi> </mrow> <mi>T</mi> </msubsup> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <msup> <mrow> <mo>&amp;lsqb;</mo> <msub> <mi>H</mi> <mrow> <mi>l</mi> <mo>&amp;times;</mo> <mi>T</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <msubsup> <mi>H</mi> <mrow> <mi>T</mi> <mo>&amp;times;</mo> <mi>l</mi> </mrow> <mi>T</mi> </msubsup> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> </mrow> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> </mrow> 2
<mrow> <msub> <mi>B</mi> <mrow> <mi>d</mi> <mo>&amp;times;</mo> <mi>l</mi> </mrow> </msub> <mo>=</mo> <msub> <mrow> <mo>(</mo> <msub> <mi>Z</mi> <mrow> <mi>t</mi> <mi>h</mi> <mi>r</mi> <mi>e</mi> <mi>e</mi> </mrow> </msub> <mo>)</mo> </mrow> <mrow> <mi>d</mi> <mo>&amp;times;</mo> <mi>T</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <msubsup> <mi>H</mi> <mrow> <mi>T</mi> <mo>&amp;times;</mo> <mi>l</mi> </mrow> <mi>T</mi> </msubsup> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <msup> <mrow> <mo>&amp;lsqb;</mo> <msub> <mi>H</mi> <mrow> <mi>l</mi> <mo>&amp;times;</mo> <mi>T</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <msubsup> <mi>H</mi> <mrow> <mi>T</mi> <mo>&amp;times;</mo> <mi>l</mi> </mrow> <mi>T</mi> </msubsup> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> </mrow> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> </mrow>
其中(Ythree)d×T(t)和(Zthree)d×T(t)分别表示Y方向上和Z方向上的振动响应信号，d为三维结构的阶数，T为时间。Od×l和Bd×l分别表示Y和Z方向上的模态振型矩阵。
步骤3)然后组装各方向的模态振型，形成识别三维模态振型的算法；最后利用单自由度系统参数识别技术快速傅里叶变换，可从模态响应矩阵中得到包括模态固有频率的信息。
6.一种时不变三维圆柱壳工作模态参数识别的实验装置，包括一个一端固支一端自由的三维薄壁金属圆柱壳、n个三向加速度传感器、数据采集卡、激振器或力锤及电脑终端；所述金属圆柱壳表面上均匀放置n个三向加速度传感器，用于采集使用所述力锤或激振器垂直于圆柱壳表面施加激励产生的三个方向的振动响应信号，所述数据采集卡与所述n个三向加速度传感器相连用于接收所述振动响应信号并发送给所述电脑终端存储，所述电脑终端包括一模态参数识别模块，所述模态参数识别模块基于电脑终端存储的振动响应信号使用结合二阶盲辨识和最小二乘广义逆方法的三维结构的工作模态参数识别方法识别三维圆柱壳的模态参数，并与由有限元分析或理论解相比较，验证结合二阶盲辨识和最小二乘广义逆方法的三维工作模态参数识别方法的正确性。
7.一种基于滑动窗二阶盲辨识的线性时变一维结构工作模态参数识别方法，其特征在于：分析的结构是线性时变结构，模态参数是随时间的改变而发生变化的，因此主要结合短时时不变理论与二阶盲辨识算法，利用二阶盲辨识算法在各窗内的统计特性，估计出各时刻的工作模态参数，然后各时刻求得的工作模态参数连接起来，从而实现时变线性结构工作模态参数识别；所述工作模态参数包括各阶模态的固有频率和模态振型。
8.根据权利要求7所述的基于滑动窗二阶盲辨识的线性时变一维结构工作模态参数识别方法，其特征在于，具体步骤如下：
步骤1)时变结构的模态参数随时间的变化而变化，根据结构动力学理论，对于自由度线性时变振动结构系统，它在物理坐标系统中的运动方程为:
<mrow> <mi>S</mi> <mo>:</mo> <mi>M</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mover> <mi>X</mi> <mo>&amp;CenterDot;&amp;CenterDot;</mo> </mover> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>+</mo> <mi>C</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mover> <mi>X</mi> <mo>&amp;CenterDot;</mo> </mover> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>+</mo> <mi>K</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mi>X</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>F</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>,</mo> <mi>t</mi> <mo>&amp;Element;</mo> <mo>&amp;lsqb;</mo> <msub> <mi>t</mi> <mrow> <mi>b</mi> <mi>e</mi> <mi>g</mi> <mi>i</mi> <mi>n</mi> </mrow> </msub> <mo>,</mo> <msub> <mi>t</mi> <mrow> <mi>e</mi> <mi>n</mi> <mi>d</mi> </mrow> </msub> <mo>&amp;rsqb;</mo> </mrow>
其中，M(t)，C(t)和分别表示随时间变化的质量矩阵，阻尼矩阵和刚度矩阵，与此同时，它们受结构的影响而随时间发生变化；表示外载荷的激励向量，和分别表示加速度响应信号，速度响应信号和位移响应信号；
步骤2)根据“短时时不变”理论，时变离散多自由度系统在时间τ∈[tbegin,tend]内，它的质量，阻尼和刚度看作是时不变的，因此，在物理坐标系中的动力学方程可表示为：
<mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msup> <mi>S</mi> <mo>&amp;prime;</mo> </msup> <mover> <mo>=</mo> <mi>&amp;Delta;</mi> </mover> <mo>{</mo> <msup> <mi>S</mi> <mo>&amp;prime;</mo> </msup> <mrow> <mo>(</mo> <mi>&amp;tau;</mi> <mo>)</mo> </mrow> <mo>:</mo> <mi>M</mi> <mrow> <mo>(</mo> <mi>&amp;tau;</mi> <mo>)</mo> </mrow> <mover> <mi>X</mi> <mo>&amp;CenterDot;&amp;CenterDot;</mo> </mover> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>+</mo> <mi>C</mi> <mrow> <mo>(</mo> <mi>&amp;tau;</mi> <mo>)</mo> </mrow> <mover> <mi>X</mi> <mo>&amp;CenterDot;</mo> </mover> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>+</mo> <mi>K</mi> <mrow> <mo>(</mo> <mi>&amp;tau;</mi> <mo>)</mo> </mrow> <mi>X</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>F</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>,</mo> <mi>t</mi> <mo>&amp;Element;</mo> <mo>&amp;lsqb;</mo> <msub> <mi>t</mi> <mrow> <mi>b</mi> <mi>e</mi> <mi>g</mi> <mi>i</mi> <mi>n</mi> </mrow> </msub> <mo>,</mo> <msub> <mi>t</mi> <mrow> <mi>e</mi> <mi>n</mi> <mi>d</mi> </mrow> </msub> <mo>&amp;rsqb;</mo> <mo>}</mo> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>&amp;tau;</mi> <mo>=</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mi>k</mi> </msub> <mo>+</mo> <msub> <mi>t</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>)</mo> </mrow> <mo>,</mo> <mi>k</mi> <mo>=</mo> <mn>0</mn> <mo>,</mo> <mn>1</mn> <mo>,</mo> <mo>...</mo> <mo>,</mo> <mi>K</mi> <mo>,</mo> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>t</mi> <mn>0</mn> </msub> <mo>=</mo> <msub> <mi>t</mi> <mrow> <mi>b</mi> <mi>e</mi> <mi>g</mi> <mi>i</mi> <mi>n</mi> </mrow> </msub> <mo>,</mo> <msub> <mi>t</mi> <mi>K</mi> </msub> <mo>=</mo> <msub> <mi>t</mi> <mrow> <mi>e</mi> <mi>n</mi> <mi>d</mi> </mrow> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced>
其中，K＝end，表示时变系统中最后时刻，S'(τ)表示为当t＝τ时刻的时不变结构，S'表示一组有限多个线性时不变结构组成时变结构的集合；
步骤3)对于时不变小阻尼结构，响应数据可以被分为有限多个部分，在第τ个部分，选取一定的窗口的长度，线性系统的模态坐标响应为：
其中，Φ(τ)和q(τ,t)分别表示当第τ个窗口的模态振型矩阵和模态响应向量；
步骤4)当时不变结构的每一阶模态固有频率ωi都不相等时，各阶模态振型之间满足归一化正交，各阶模态响应相互不相关，如下：
<mrow> <msubsup> <mover> <mi>&amp;phi;</mi> <mo>&amp;RightArrow;</mo> </mover> <mi>i</mi> <mi>T</mi> </msubsup> <mrow> <mo>(</mo> <mi>&amp;tau;</mi> <mo>)</mo> </mrow> <msub> <mover> <mi>&amp;phi;</mi> <mo>&amp;RightArrow;</mo> </mover> <mi>j</mi> </msub> <mrow> <mo>(</mo> <mi>&amp;tau;</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <mn>0</mn> <mo>,</mo> <mi>i</mi> <mo>&amp;NotEqual;</mo> <mi>j</mi> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mn>1</mn> <mo>,</mo> <mi>i</mi> <mo>=</mo> <mi>j</mi> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>,</mo> <mi>&amp;tau;</mi> <mo>=</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <mrow> <mo>(</mo> <msub> <mi>t</mi> <mi>k</mi> </msub> <mo>+</mo> <msub> <mi>t</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>)</mo> </mrow> </mrow>
其中，E(q(τ,t)qT(τ,t))表示两模态振型的期望，Λ"n×n表示阶数为n的对角矩阵；
步骤5)假设在一个很短的时间段内，系统被看成是短时时不变的，也就是说，时间被划分成有限段，在每一个时间段内，系统被认为是短时时不变的，从而利用时不变结构的工作模态参数识别算法，识别出该时间段的工作模态参数，窗口向右边滑动，即计算下一个时间段内工作模态参数，以此类推，最后将每个时间段按照时间顺序排列起来，从而形成时变结构的模态参数；
其中，响应数据的限定记忆长度为L，n表示的是传感器的个数，T表示采样时间；
步骤6)对于采集到的振动结构的位移响应数据，其模态坐标表示如下：
<mrow> <msubsup> <mi>X</mi> <mi>L</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </msubsup> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>&amp;ap;</mo> <msubsup> <mi>&amp;Phi;</mi> <mi>L</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </msubsup> <msubsup> <mi>Q</mi> <mi>L</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </msubsup> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>,</mo> <mi>i</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mn>2</mn> <mo>,</mo> <mo>...</mo> <mo>,</mo> <mi>T</mi> <mo>+</mo> <mn>1</mn> <mo>-</mo> <mi>L</mi> </mrow>
其中，表示振动结构的位移响应数据，是模态振型矩阵，表示模态坐标响应，当系统的各阶的模态固有频率不相等是，各阶模态振型相互归一化正交，各阶模态坐标响应相互独立，如下：
<mrow> <msup> <mrow> <mo>(</mo> <msubsup> <mi>&amp;Phi;</mi> <mi>L</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </msubsup> <mo>)</mo> </mrow> <mi>T</mi> </msup> <msubsup> <mi>&amp;Phi;</mi> <mi>L</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </msubsup> <mo>=</mo> <msub> <mi>I</mi> <mrow> <mi>n</mi> <mo>&amp;times;</mo> <mi>n</mi> </mrow> </msub> </mrow>
<mrow> <mi>E</mi> <mo>&amp;lsqb;</mo> <msubsup> <mi>Q</mi> <mi>L</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </msubsup> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <msup> <mrow> <mo>(</mo> <msubsup> <mi>Q</mi> <mi>L</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </msubsup> <mo>(</mo> <mi>t</mi> <mo>)</mo> <mo>)</mo> </mrow> <mi>T</mi> </msup> <mo>&amp;rsqb;</mo> <mo>=</mo> <msubsup> <mi>&amp;Lambda;</mi> <mrow> <mi>n</mi> <mo>&amp;times;</mo> <mi>n</mi> </mrow> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </msubsup> </mrow>
其中，表示两模态坐标响应的期望，Λ"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阶模态的平均频率变化量为△fL(i)：
<mrow> <msub> <mi>&amp;Delta;f</mi> <mi>L</mi> </msub> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mi>L</mi> <mi>f</mi> </mfrac> <mo>&amp;times;</mo> <mfrac> <mrow> <msub> <mi>f</mi> <mrow> <mi>e</mi> <mi>n</mi> <mi>d</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>f</mi> <mrow> <mi>b</mi> <mi>e</mi> <mi>g</mi> <mi>i</mi> <mi>n</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo> </mrow> </mrow> <mrow> <msub> <mi>t</mi> <mrow> <mi>e</mi> <mi>n</mi> <mi>d</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>t</mi> <mrow> <mi>b</mi> <mi>e</mi> <mi>g</mi> <mi>i</mi> <mi>n</mi> </mrow> </msub> </mrow> </mfrac> </mrow>
其中，变量fend(i),fbegin(i),tend,tbegin分别表示第i阶模态的结束频率，开始频率，整个数据的结束时间和开始时间；
步骤10)滑动窗的窗口的长度L是振动响应数据的采样频率和频率分辨率相除得到的，△f的取值不能取太小，因为不能反映出频率的变化，△f的取值也不能取太大，因为平均频率变化量△fL(i)不能被识别出来，与此同时，△fL(i)不能比△f的取值大太多，否则，线性时变结构在一个滑动窗的窗口的长度L内不能被视为时不变的了；
步骤11)模态固有频率的计算，是由单自由度技术计算得到的，利用的是快速傅里叶变换，而快速傅里叶变换的计算复杂度跟滑动窗窗口的长度有直接的关系，根据快速傅里叶变换的原理，滑动窗窗口的长度应该满足如下：
L＝2α,α＝1,2,…。
CN201710500228.0A 2017-06-27 2017-06-27 基于二阶盲辨识的线性结构工作模态参数识别方法及装置 Active CN107357977B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710500228.0A CN107357977B (zh) 2017-06-27 2017-06-27 基于二阶盲辨识的线性结构工作模态参数识别方法及装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710500228.0A CN107357977B (zh) 2017-06-27 2017-06-27 基于二阶盲辨识的线性结构工作模态参数识别方法及装置

Publications (2)

Publication Number Publication Date
CN107357977A true CN107357977A (zh) 2017-11-17
CN107357977B CN107357977B (zh) 2021-03-09

Family

ID=60272708

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710500228.0A Active CN107357977B (zh) 2017-06-27 2017-06-27 基于二阶盲辨识的线性结构工作模态参数识别方法及装置

Country Status (1)

Country Link
CN (1) CN107357977B (zh)

Cited By (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108090270A (zh) * 2017-12-12 2018-05-29 华南理工大学 一种基于形态学滤波和盲源分离的暂态振荡参数识别方法
CN108491608A (zh) * 2018-03-06 2018-09-04 大连理工大学 传感器数量不完备时结构模态识别的稀疏分量分析方法
CN109238447A (zh) * 2018-09-12 2019-01-18 西北工业大学 一种系绳振动信号的盲源分离方法
CN109992834A (zh) * 2019-03-05 2019-07-09 中国人民解放军海军勤务学院 改进型盲源分离的结构模态识别方法
CN110598173A (zh) * 2019-08-31 2019-12-20 中国人民解放军陆军工程大学 基于运用效能的作战系统操作行为分析方法
CN110705041A (zh) * 2019-09-12 2020-01-17 华侨大学 一种基于easi的线性结构工作模态参数识别方法
CN110901689A (zh) * 2019-11-19 2020-03-24 华东交通大学 一种基于模态识别的轨道结构扣件松脱检测方法
CN112255121A (zh) * 2020-12-23 2021-01-22 天津航天瑞莱科技有限公司 基于Matlab叶片疲劳极限评估方法
CN112506058A (zh) * 2020-12-03 2021-03-16 华侨大学 一种线性时变结构的工作模态参数识别方法及系统
CN112905958A (zh) * 2021-01-27 2021-06-04 南京国电南自电网自动化有限公司 基于测控装置的短时数据窗遥测数据状态辨识方法及系统
CN113076517A (zh) * 2021-04-01 2021-07-06 重庆大学 基于希尔伯特变换的土木工程结构动态监测相位评估方法
CN113686528A (zh) * 2021-07-28 2021-11-23 华南理工大学 一种结构-tld系统的子系统动力特性检测方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103528667A (zh) * 2013-10-23 2014-01-22 东北大学 一种基于激光扫描的圆柱壳模态振型测试装置及方法
CN105159865A (zh) * 2015-07-01 2015-12-16 华侨大学 复杂声振模拟实验环境下进行不相关多源频域载荷识别的装置和方法
CN106202977A (zh) * 2016-08-17 2016-12-07 华南理工大学 一种基于盲源分离算法的低频振荡模式分析方法
CN106446502A (zh) * 2016-07-21 2017-02-22 华侨大学 带遗忘因子的特征向量递推的时变工作模态在线识别方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103528667A (zh) * 2013-10-23 2014-01-22 东北大学 一种基于激光扫描的圆柱壳模态振型测试装置及方法
CN105159865A (zh) * 2015-07-01 2015-12-16 华侨大学 复杂声振模拟实验环境下进行不相关多源频域载荷识别的装置和方法
CN106446502A (zh) * 2016-07-21 2017-02-22 华侨大学 带遗忘因子的特征向量递推的时变工作模态在线识别方法
CN106202977A (zh) * 2016-08-17 2016-12-07 华南理工大学 一种基于盲源分离算法的低频振荡模式分析方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
A.BELOUCHRANI等: "《Second-order blind separation of temporally correlated sources》", 《PROC. INT. CON# DIGITAL SIGNAL PROCESSING. CITESEER》 *
康兴无等: "《振动筛试验模态分析》", 《2012年LMS中国用户大会论文集》 *
张亚飞等: "《最小二乘广义逆求解方法研究及应用》", 《应用科技》 *
王成等: "《Operational modal analysis for slow linear time-varying structures based on moving window second order blind identification》", 《SIGNAL PROCESSING》 *

Cited By (19)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108090270A (zh) * 2017-12-12 2018-05-29 华南理工大学 一种基于形态学滤波和盲源分离的暂态振荡参数识别方法
CN108491608B (zh) * 2018-03-06 2021-06-08 大连理工大学 传感器数量不完备时结构模态识别的稀疏分量分析方法
CN108491608A (zh) * 2018-03-06 2018-09-04 大连理工大学 传感器数量不完备时结构模态识别的稀疏分量分析方法
CN109238447A (zh) * 2018-09-12 2019-01-18 西北工业大学 一种系绳振动信号的盲源分离方法
CN109992834A (zh) * 2019-03-05 2019-07-09 中国人民解放军海军勤务学院 改进型盲源分离的结构模态识别方法
CN110598173A (zh) * 2019-08-31 2019-12-20 中国人民解放军陆军工程大学 基于运用效能的作战系统操作行为分析方法
CN110598173B (zh) * 2019-08-31 2023-06-09 中国人民解放军陆军工程大学 基于运用效能的作战系统操作行为分析方法
CN110705041A (zh) * 2019-09-12 2020-01-17 华侨大学 一种基于easi的线性结构工作模态参数识别方法
CN110705041B (zh) * 2019-09-12 2022-12-23 华侨大学 一种基于easi的线性结构工作模态参数识别方法
CN110901689B (zh) * 2019-11-19 2020-09-29 华东交通大学 一种基于模态识别的轨道结构扣件松脱检测方法
CN110901689A (zh) * 2019-11-19 2020-03-24 华东交通大学 一种基于模态识别的轨道结构扣件松脱检测方法
CN112506058A (zh) * 2020-12-03 2021-03-16 华侨大学 一种线性时变结构的工作模态参数识别方法及系统
CN112255121B (zh) * 2020-12-23 2021-05-14 天津航天瑞莱科技有限公司 基于Matlab叶片疲劳极限评估方法
CN112255121A (zh) * 2020-12-23 2021-01-22 天津航天瑞莱科技有限公司 基于Matlab叶片疲劳极限评估方法
CN112905958A (zh) * 2021-01-27 2021-06-04 南京国电南自电网自动化有限公司 基于测控装置的短时数据窗遥测数据状态辨识方法及系统
CN112905958B (zh) * 2021-01-27 2024-04-19 南京国电南自电网自动化有限公司 基于测控装置的短时数据窗遥测数据状态辨识方法及系统
CN113076517A (zh) * 2021-04-01 2021-07-06 重庆大学 基于希尔伯特变换的土木工程结构动态监测相位评估方法
CN113076517B (zh) * 2021-04-01 2022-09-30 重庆大学 基于希尔伯特变换的土木工程结构动态监测相位评估方法
CN113686528A (zh) * 2021-07-28 2021-11-23 华南理工大学 一种结构-tld系统的子系统动力特性检测方法

Also Published As

Publication number Publication date
CN107357977B (zh) 2021-03-09

Similar Documents

Publication Publication Date Title
CN107357977B (zh) 基于二阶盲辨识的线性结构工作模态参数识别方法及装置
Zhou et al. Blind source separation based vibration mode identification
CN107271127B (zh) 基于自迭代主元抽取的工作模态参数识别方法及装置
CN108594660B (zh) 一种时不变结构的工作模态参数识别方法和系统
Hu et al. Model order determination and noise removal for modal parameter estimation
Poozesh et al. Modal parameter estimation from optically-measured data using a hybrid output-only system identification method
Guan et al. Data-driven methods for operational modal parameters identification: A comparison and application
CN114912547A (zh) 一种线性时变结构的欠定工作模态参数识别方法和系统
CN114925526B (zh) 一种结合多工况响应的结构模态参数识别方法
Brewick et al. On the application of blind source separation for damping estimation of bridges under traffic loading
Prawin et al. Damage detection in nonlinear systems using an improved describing function approach with limited instrumentation
Lin et al. Mechanism of principal component analysis in structural dynamics under ambient excitation
CN116911049B (zh) 单段振动响应数据的结构模态参数不确定性量化方法
Pandiya et al. Direct estimation of residues from rational-fraction polynomials as a single-step modal identification approach
CN106446502A (zh) 带遗忘因子的特征向量递推的时变工作模态在线识别方法
Nejad et al. An investigation on the capability of proper orthogonal modes in determining the natural frequencies and damping ratios of linear structural systems
Kim et al. Comparison study of output-only subspace and frequency-domain methods for system identification of base excited civil engineering structures
Kaljević et al. Stochastic boundary elements in elastostatics
CN112506058A (zh) 一种线性时变结构的工作模态参数识别方法及系统
Ghobadi et al. AOSID: An analytical solution to the output‐only system identification problem to estimate physical parameters and unknown input simultaneously
Adhikari Random eigenvalue problems revisited
Mirzazadeh et al. Uncertainty quantification in polysilicon MEMS through on-chip testing and reduced-order modelling
Cara et al. Estimation of modal parameters in structures using multiple time history records
Panda et al. Mastering Complex Modes: A New Method for Real-Time Modal Identification of Vibrating Systems
Gibanica Experimental-analytical dynamic substructuring: a state-space approach

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant