CN107561934B - 基于多参考时域传递率的仅输出结构模态参数辨识方法 - Google Patents

基于多参考时域传递率的仅输出结构模态参数辨识方法 Download PDF

Info

Publication number
CN107561934B
CN107561934B CN201710732229.8A CN201710732229A CN107561934B CN 107561934 B CN107561934 B CN 107561934B CN 201710732229 A CN201710732229 A CN 201710732229A CN 107561934 B CN107561934 B CN 107561934B
Authority
CN
China
Prior art keywords
formula
matrix
estimated
time domain
transfer rate
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.)
Expired - Fee Related
Application number
CN201710732229.8A
Other languages
English (en)
Other versions
CN107561934A (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.)
Beijing Institute of Technology BIT
Original Assignee
Beijing Institute of Technology BIT
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 Beijing Institute of Technology BIT filed Critical Beijing Institute of Technology BIT
Priority to CN201710732229.8A priority Critical patent/CN107561934B/zh
Publication of CN107561934A publication Critical patent/CN107561934A/zh
Application granted granted Critical
Publication of CN107561934B publication Critical patent/CN107561934B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Other Investigation Or Analysis Of Materials By Electrical Means (AREA)
  • Testing Of Devices, Machine Parts, Or Other Structures Thereof (AREA)

Abstract

本发明公开的基于多参考时域传递率的仅输出结构模态参数辨识方法，属于结构动力学技术领域。本发明将采集的工程结构响应信号分成非参考输出信号与参考输出信号两组，分别作为多参考时域传递率模型的输入和输出，能处理非白输入条件下的结构模态辨识问题，结构不受载荷形式限制；本发明将传统在频域定义的传递率函数拓展到时域，避免结构测量输出时域数据向频域的转换，直接采用时域数据对工程结构进行模态分析，因此，避免转换过程中窗函数的选择，消除转化过程中存在的泄露、频率分辨率降低等误差及不良影响，提高模态参数辨识的精度，同时降低对使用者专业知识的要求。本发明辨识的结构模态参数能够应用于工程结构力学领域解决实际工程问题。

Description

基于多参考时域传递率的仅输出结构模态参数辨识方法
技术领域
本发明涉及一种仅输出结构模态参数辨识方法，尤其涉及一种基于多参考时域传递率的仅输出结构模态参数辨识方法，属于结构动力学技术领域。
背景技术
在日常生活和工业生产中，处于工作条件下的工程结构均会受到环境的激励作用，从而产生振动。当工程结构的固有模态受到外界激励时，会使工程结构产生大幅振动甚至遭到破坏。因此，在进行工程结构的初始设计时，应对其进行动力学分析，考虑工程结构对振动的敏感程度。作为结构动力学特性分析的重要方法和途径，结构模态参数辨识研究显得十分重要。
近几十年来，作为结构动力学特性分析的重要方法和途径，结构模态参数辨识研究成为结构动力学研究的重点之一。结构模态参数辨识可以辨识结构的模态频率、模态振型和模态阻尼，这些参数具有重要的物理意义，除指导工程结构的结构设计外，还可为结构健康监测、结构故障诊断、结构振动控制等方面的应用提供有力的支持。
目前，结构模态参数辨识方法已得到了广泛的关注与发展。根据采用数据集合的不同，目前结构模态参数辨识方法主要分为两类：
第一类是基于输入-输出数据的实验模态分析。该类方法通过传感器采集结构的输入-输出数据，并根据输入-输出数据提取出模态参数。
该类方法中的频域方法大多基于输入-输出数据获得频率响应函数(FrequencyResponse Function,FRF)，之后根据FRF获取模态参数。峰值提取法是最简单的频域模态参数辨识方法，通过提取FRF幅值的峰值确定结构的模态频率，并通过半功率点确定结构的阻尼比。圆拟合方法基于频域的结构模态模型，将FRF绘制在Nyquist图上，FRF的峰值附近的点形成一段圆弧，通过拟合该段圆弧获得结构的阻尼和模态频率。1981年，Phillips等提出复模态指示函数法。通过对FRF矩阵进行奇异值分解，绘制其对于频率的对数幅值图，图中峰值对应结构的模态频率。另外，FRF可写成矩阵分式形式，形成了一系列的参数化辨识方法，如Van der Auweraer等在2001年提出的最小二乘复频域法。
上述基于输入-输出数据的结构模态参数辨识最大的缺陷在于输入数据必须已知。在真实工作状态下，工程结构上作用的环境激励往往是随机且不可测量的；对于工程中的大型结构，传感器的布置十分复杂，且激振设备选择困难，如桥梁、建筑物等。因此，在真实工作状态下，基于输入-输出数据的模态参数辨识方法不适用。
第二类方法是基于仅输出数据的工作模态分析方法。该方法中的输入通常假设为某种形式的随机过程，如白噪声，弥补了第一类方法需要已知输入数据的缺陷。该类方法主要包括时间序列法、随机子空间辨识法、矩阵分式法等。
基于仅输出数据的工作模态分析方法最大的缺陷在于其输入必须为随机激励，如白噪声。在真实工作状态下，这一假设并不能严格成立。例如，转动的叶轮、火箭储箱中晃动的液体等均能产生包含谐波的激励。当激励中的谐波频率与结构固有模态频率相距较远时，谐波将被辨识为虚假模态；当激励中的谐波频率与结构固有模态频率相距很近时，谐波可能会严重影响模态参数的辨识精度。
针对带有谐波分量的激励，可采用滤波的方法将谐波滤掉，但该方法只在谐波频率与结构固有频率相隔较远时才有效，且滤波前需已知谐波频率的大致范围。近年来，采用传递率函数进行仅输出结构模态参数辨识得到广泛的关注与发展。传递率定义为非参数输出通道频域输出数据与参考输出通道频域输出数据的比值，类似于FRF的定义，只是将FRF中的输入变成参考通道的输出，FRF中的响应变成非参考通道的输出。2005年，P.Guillaume等人首次提出基于传递率的工作模态分析(TOMA)方法。2013年，C.Devriendt将所有测量的单参考传递率函数写成矩阵形式，并辨识出系统极点，提高了TOMA方法的稳定性。2010年，C.Devriendt将单参考传递率函数拓展为多参考形式，并采用两个不同工况下的传递率函数成功辨识出系统的模态参数。
基于传递率的工作模态分析方法成功地解决了激励为非白噪声时的模态参数辨识问题，且无须已知激励中谐波成分的频率范围，应用比较方便。
目前，基于传递率的工作模态分析方法均依赖于频率数据。在测得工程结构的时域输出数据之后，必须对时域数据进行转换，使其变成频域数据，转换过程中依赖于傅里叶变换。由于工作状态下测量的时域数据一般不是整周期的，傅里叶变换过程中存在泄漏，导致频域数据存在较大误差。目前普遍采用的抑制泄漏的方法为添加窗函数，但窗函数的选择须十分谨慎，因为添加窗函数会降低频率分辨率，且会增加辨识结果对噪声的敏感度。另外，窗函数会使得相邻频域数据之间的相互作用增强。
总之，目前基于传递率的工作模态分析方法采用频域数据，须对测量的时域数据进行傅里叶变换，变换过程中不可避免使用窗函数。由于窗函数选择需要考虑频率分辨率、噪声敏感度、相邻数据间的相互作用等因素，使得窗函数选择十分复杂。在实际使用中，窗函数的选择一般根据经验选取，没有简单易行的确定方法，在实际应用时很复杂，要求使用者具有深入的专业知识背景。
发明内容
针对基于传递率的工作模态分析方法存在的上述技术问题，本发明公开的基于多参考时域传递率的仅输出结构模态参数辨识方法要解决的问题是，将传统的在频域定义的传递率函数拓展到时域，避免结构测量输出时域数据向频域的转换，直接采用时域数据对工程结构进行模态分析，因此，避免转换过程中窗函数的选择，彻底消除转化过程中存在的泄露、频率分辨率降低等误差及不良影响，提高模态参数辨识的精度，同时降低对使用者专业知识背景的要求。
本发明的目的是通过下述技术方案实现的：
本发明公开的基于多参考时域传递率的仅输出结构模态参数辨识方法，首先利用信号采集系统采集得到工程结构的响应信号，将响应信号分成参考输出信号和非参考输出信号两组。根据参考输出信号与非参考输出信号，构造多参考时域传递率模型。根据赤池信息量准则确定多参考时域传递率模型阶数n，完成多参考时域传递率模型的构造。根据构造的多参考时域传递率模型，估计多参考时域传递率模型系数矩阵。选取多个不同的工况，结构受到的激励幅值、位置或形式不同，均当做不同工况。根据构建的多参考时域传递率模型获得所有不同工况下的待估系数矩阵。根据所有不同工况下估计出的模型系数矩阵，构造并求解广义特征值问题，根据得到的特征值与特征向量即求得结构的模态参数。最后将结构模态参数应用于工程结构设计、工程结构状态监测、工程结构振动控制等领域。
所述的仅输出结构模态参数包括模态频率、模态阻尼及最大值归一化的模态向量。
将本发明公开的基于多参考时域传递率的仅输出结构模态参数辨识方法得到的结构模态参数应用于工程结构力学领域解决实际工程问题。针对任意激励下的工程结构进行模态分析，检测是否符合工程结构标准及隔振等要求，提高结构设计精度与可靠度。结构模态参数对结构故障具有极高的敏感度，将其用于结构的健康状态监测，提高结构健康状态监测的实时性与可靠性。结构模态参数用于结构振动控制领域，根据结构模态参数确定结构对环境激励的敏感程度及发生共振的频率范围，降低结构振动幅值，完成结构振动控制。本发明公开的基于多参考时域传递率的仅输出结构模态参数辨识方法具有广泛的应用前景与效益。
本发明公开的基于多参考时域传递率的仅输出结构模态参数辨识方法，包括以下步骤：
步骤1：采集工程结构的响应信号，并将响应信号进行分组。
利用信号采集系统采集得到工程结构的响应信号，将响应信号分成两组，其中一组为参考输出信号，记为xr，另外一组为非参考输出信号，记为xl；下标r表示参考信号，l表示非参考信号。
其中，
Figure BDA0001387368700000041
为参考输出信号，
Figure BDA0001387368700000044
表示实数集，Nr表示参考通道数，Nt为采集数据的长度；
Figure BDA0001387368700000042
为非参考输出信号，No为输出通道总数。Nr的取值大于等于2，小于等于
Figure BDA0001387368700000043
步骤2：根据步骤1中的参考输出信号xr与非参考输出信号xl，构造多参考时域传递率模型。
将参考输出信号xr当做模型输入，非参考输出信号xl当做模型输出，构造含有外部激励的自回归时间序列模型，如式(1)所示：
Figure BDA0001387368700000051
式(1)中，n为模型阶数；x(t)表示第t个时刻的信号，
Figure BDA0001387368700000052
为第t个时刻的残差序列；Ai表示第i阶AR系数，Bj表示第j阶外部激励系数。Ai与Bj均为待估参数，写成如式(2)所示的矩阵形式：
Θ＝[AT BT]T (2)
式(2)中，Θ为由待估系数矩阵Ai与Bj构成的新矩阵，T表示矩阵转置运算，其中A与B分别定义如式(3)所示：
A＝[A0 A1 … An]T，B＝[B0 B1 … Bn]T (3)
式(1)所示的含有外部激励的自回归时间序列模型即为构造的多参考时域传递率模型。
步骤3：确定式(1)中多参考时域传递率模型阶数n。
依据赤池信息量准则(Akaike Information Criterion，AIC)确定式(1)中的多参考时域传递率模型阶数n。当多参考时域传递率模型阶数n达到最优时，AIC的值取最小。因此，赤池信息量准则AIC值取最小时对应的多参考时域传递率模型阶数即为n值。
步骤4：根据步骤2中构造的多参考时域传递率模型及步骤3中确定的模型阶数n，估计其待估系数矩阵Θ。
步骤4中估计其待估系数矩阵Θ优选选择采用最小二乘法或最大似然法。
步骤4.1：根据步骤2中构造的多参考时域传递率模型及步骤3中确定的模型阶数n，采用最小二乘法估计多参考时域传递率模型待估系数矩阵Θ，具体实现方法如下：
步骤4.1.1：根据步骤2中构造的多参考时域传递率模型，写出最小二乘法对应的待估系数矩阵Θ的费用函数。
最小二乘法估计待估系数矩阵Θ的费用函数JLS如式(4)所示：
Figure BDA0001387368700000061
式(4)中，j表示待估系数矩阵Θ的第j列；T表示矩阵的转置运算；Φ为回归矩阵，形式如式(5)所示：
Φ＝[Φl Φr] (5)
式(5)中的Φl与Φr的定义分别如式(6)所示：
Figure BDA0001387368700000062
步骤4.1.2：根据步骤4.1.1中如式(4)所示的费用函数JLS，使费用函数取得最小值的待估参数值即为待估系数矩阵Θ的估计值，求出待估系数矩阵A与B。
根据式(4)所示的费用函数JLS，待估系数矩阵Θ根据式(7)求得：
Figure BDA0001387368700000063
式(7)中，
Figure BDA0001387368700000064
（·）-1表示矩阵求逆运算。A′由式(8)计算：
A′＝（D′）-1b (8)
式(8)中，D′与b分别定义如式(9)所示：
Figure BDA0001387368700000065
式(9)中，[·][:,:]表示取矩阵的一部分，例如[D][a:b,c:d]表示取矩阵D中行号为a到b、列号为c到d的部分；矩阵D定义如式(10)所示：
D＝L-SR-1SH (10)
式(10)中的L定义如式(11)所示：
L＝Φl TΦl (11)
即完成采用最小二乘法估计多参考时域传递率模型待估系数矩阵Θ步骤。
步骤4.2：根据步骤2中构造的多参考时域传递率模型及步骤3中确定的模型阶数n，采用最大似然法估计多参考时域传递率模型待估系数矩阵Θ，具体实现方法如下：
步骤4.2.1：根据步骤2中构造的多参考时域传递率模型，写出最大似然法对应的待估系数矩阵Θ的费用函数。
采用最大似然法时，将式(1)中的系数矩阵A0取为单位矩阵。根据对数似然函数形式，最大似然估计法的费用函数如式(12)所示：
Figure BDA0001387368700000071
式(12)中，
Figure BDA0001387368700000072
为残差序列e(t)的方差矩阵，且对任意时刻t均相同；|·|表示行列式，tr(·)表示矩阵的迹。
Figure BDA0001387368700000073
为回归矩阵，形式如式(13)所示：
Figure BDA0001387368700000074
式(13)中的
Figure BDA0001387368700000075
与
Figure BDA0001387368700000076
的定义分别如式(14)所示：
Figure BDA0001387368700000077
式(12)中的矩阵X如式所示：
Figure BDA0001387368700000078
式(12)中最大似然法对应的待估系数矩阵
Figure BDA0001387368700000079
如式(16)所示：
Figure BDA00013873687000000710
步骤4.2.2：根据步骤4.2.1中如式(12)所示的费用函数
Figure BDA00013873687000000711
使费用函数取得最大值的待估参数值即为待估系数矩阵
Figure BDA00013873687000000712
的估计值，求出待估系数矩阵A与B。
根据如式(12)所示的费用函数
Figure BDA0001387368700000081
避免最大似然法中的非线性迭代求解，将残差序列e(t)的方差矩阵∑看成与待估系数矩阵
Figure BDA0001387368700000082
无关，得到的待估系数矩阵与最小二乘法得到的结果是等价的。待估系数矩阵
Figure BDA0001387368700000083
由式(17)求得：
Figure BDA0001387368700000084
式(17)中，(·)-1表示矩阵求逆运算。求出待估系数矩阵
Figure BDA0001387368700000085
后，结合式(16)并注意
Figure BDA0001387368700000086
表示No-Nr维单位矩阵，即求出式(3)中完整的待估系数矩阵A与B。
即完成采用最大似然法估计多参考时域传递率模型待估系数矩阵Θ步骤。
步骤5：选取多个不同工况，针对每个工况完成多参考时域传递率模型待估系数矩阵Θ的估计。
选取多个不同的工况，对每一个工况，均执行步骤1、步骤2及步骤4，步骤3只须运行一次，对所有工况均取相同的模型阶数n，获得所有不同工况下的待估参数Ak与Bk。
下标k表示第k个工况，k＝1,2,…,Nk，Nk为不同工况总数。结构受到的激励幅值、位置或形式不同，均当做不同工况，且要求Nk≥2。为区别不同工况的待估系数矩阵，式(3)改写成式(18)的形式：
Ak＝[Ak,0 Ak,1 … Ak,n]T，Bk＝[Bk,0 Bk,1 … Bk,n]T (18)
即完成每个工况下多参考时域传递率模型待估系数矩阵Θ的估计。
步骤6：根据步骤5中得到的所有不同工况下待估系数矩阵Ak与Bk，完成仅输出结构模态参数辨识。
所述的仅输出结构模态参数包括模态频率、模态阻尼及最大值归一化的模态向量。
步骤6.1：根据步骤5中得到的所有不同工况下待估系数矩阵Θk，求解如式(19)所示的广义特征值问题：
Figure BDA0001387368700000087
式(19)中，λm为矩阵Γ的第m个特征值，向量
Figure BDA0001387368700000088
表示λm对应的矩阵Γ的特征向量；矩阵Γ、向量
Figure BDA0001387368700000091
分别定义如式(20)所示：
Figure BDA0001387368700000092
式(20)中，φm为结构的第m阶未最大值归一化的模态向量；矩阵C′i(其中i＝0,1,…,n-1)定义如式(21)所示：
Figure BDA0001387368700000093
式(21)中，矩阵
Figure BDA0001387368700000094
(其中j＝0,1,…,n)定义如式(22)所示：
Figure BDA0001387368700000095
式(22)中，矩阵A与B的第一个下标表示工况数，第二个下标表示步骤2式(1)中的阶数。
式(21)中的矩阵求逆运算要求矩阵Cn为方阵，即要求Nk＝No/(No-Nr)。当Cn不为方阵时，可改用伪逆运算，如式(23)所示：
Figure BDA0001387368700000096
上标“+”表示矩阵伪逆运算。
步骤6.2：根据步骤6.1中求得的特征值λm与未最大值归一化的模态向量φm求出结构的模态频率、模态阻尼及最大值归一化的模态向量。
系统的特征值写成λm＝cm+idm的形式，则结构第m阶模态频率ωm、第m阶模态阻尼ξm与第m阶最大值归一化的模态向量
Figure BDA0001387368700000097
分别定义如式(24)所示：
Figure BDA0001387368700000098
式(24)中，ωm的单位为rad/s，abs(·)表示对向量的每个元素取绝对值，max(·)表示取向量中最大的元素，αm与βm分别定义如式(25)所示：
Figure BDA0001387368700000099
根据式(24)完成基于多参考时域传递率的仅输出结构模态参数辨识。
步骤7：将步骤6辨识的结构模态参数应用于工程结构力学领域解决实际工程问题。优选选择应用于结构设计、状态监测、振动控制领域。
步骤7.1：将步骤6辨识的结构模态参数结构设计领域，实现针对任意激励下的工程结构进行模态分析，根据步骤6辨识的结构模态参数，得到工程结构的固有频率，以此检测工程结构对环境激励的敏感程度，检测是否符合工程结构标准及隔振等要求，提高结构设计精度与可靠度。
步骤7.2：将步骤6辨识的结构模态参数结构状态监测领域，步骤6辨识的结构模态参数对结构故障具有极高的敏感度，根据步骤6辨识的结构模态参数对结构健康状态进行监测，并检测结构是否出现故障，提高结构健康状态监测的实时性与可靠性。
步骤7.3：将步骤6辨识的结构模态参数结构振动控制领域，根据步骤6辨识的结构模态参数确定结构对环境激励的敏感程度及发生共振的频率范围，降低结构振动幅值，完成结构振动控制。
有益效果：
1、本发明公开的基于多参考时域传递率的仅输出结构模态参数辨识方法，属于结构动力学中工作模态分析领域，将采集的工程结构响应信号分成非参考输出信号与参考输出信号两组，并分别作为多参考时域传递率模型的输入和输出，能处理非白输入条件下的结构模态辨识问题，结构不受载荷形式的限制，拓展应用范围。
2、本发明公开的基于多参考时域传递率的仅输出结构模态参数辨识方法，在对结构进行模态分析前，不用已知结构所受载荷中谐波成分的频率范围，与已有技术中的滤波等处理载荷中谐波分量的方法相比，更容易操作，限制条件更少，即使在完全缺乏结构所受激励先验知识的情况下也能进行操作。
3、本发明公开的基于多参考时域传递率的仅输出结构模态参数辨识方法，将已有的在频域内定义的传递率函数拓展到时域。在实际应用中，使用者通过加速度传感器等测量的结构时域输出信号可直接使用，不用通过傅里叶变换变换到频域，因此不用选择窗函数，相对于已有技术中的方法使用更方便，即使在缺乏信号处理等专业知识背景的情况下也能进行操作。
4、本发明公开的基于多参考时域传递率的仅输出结构模态参数辨识方法，能够避免测量的时域数据向频域数据的转换，不仅简单易行，而且去除傅里叶变换中出现的泄露、频率分辨率降低、数据间相互作用加强等误差及不良影响，使得对工程结构的模态分析精度更高。
附图说明
图1为本发明公开的一种基于多参考时域传递率模型的仅输出结构模态参数辨识方法的流程图；
图2为具体实施方式中的四自由度弹簧-粘性阻尼器-质量系统；
图3为具体实施方式中采用多参考时域传递率模型得到的AIC随多参考时域传递率模型阶数的变化曲线；
图4为具体实施方式中的四自由度系统的两个不同工况下的载荷作用位置示意图。其中，图4(A)、图4(B)分别表示工况1和工况2下的载荷作用位置；
图5为具体实施方式中作用在四自由度系统上的非白载荷形式示意图；其中，图5(A)表示载荷的时域曲线，图5(B)表示载荷的功率谱曲线；
图6为具体实施方式中工况1下质量块m1的加速度响应的功率谱曲线，其中曲线上的符号“×”所在的峰值表示该峰值是由作用在系统上的载荷中的谐波分量造成的，并非由系统的动力学特性造成的；
图7为具体实施方式中采用本发明中的多参考时域传递率方法与采用已有技术中的频域方法求得的系统模态频率和阻尼结果，共进行200次Monte Carlo仿真。图7(A)表示得到的频率结果，图7(B)表示得到的阻尼结果；其中，三角形表示由本发明中多参考时域传递率方法得到的结果，空心圆圈表示由已有技术中的频域方法得到的结果。黑色的“+”号表示理论值。
图8为具体实施方式中采用本发明中的多参考时域传递率方法与采用已有技术中的频域方法求得的系统最大值归一化模态向量结果，共进行200次Monte Carlo仿真，图8中绘制的是200次模拟的平均结果。图8(A)为采用本发明中的多参考时域传递率方法辨识出的结果，图8(B)为采用已有技术中的频域方法辨识出的结果。其中，黑色实线、黑色虚线、黑色点画线、黑色点线分别表示辨识出的系统第一阶、第二阶、第三阶、第四阶最大值归一化模态向量；空心圆圈、空心正方形、黑色×、黑色+分别表示系统第一阶、第二阶、第三阶、第四阶最大值归一化模态向量的理论值。
具体实施方式
为了更好地说明本发明的目的和优点，下面通过对一个非白激励下的弹簧-粘性阻尼器-质量块系统进行动力学分析，对本发明做出详细解释。
实施例1：
本实施例的四自由度弹簧-粘性阻尼器-质量系统，如图2所示。本实施例的四自由度弹簧-粘性阻尼器-质量系统包含四个质量块m1、m2、m3、m4，四个粘性阻尼器c1、c2、c3、c4，四个弹簧k1、k2、k3、k4。
四自由度系统的各个变量取值如表1所示。
表1四自由度时变系统参数
Figure BDA0001387368700000121
图2所示的四自由度系统的动力学控制方程如式(26)所示：
Figure BDA0001387368700000122
式(26)中，M、C、K分别为质量矩阵、阻尼矩阵和刚度矩阵，x为系统的响应向量，f(t)为外界激励。M、C、K如式(27)所示：
Figure BDA0001387368700000123
式(27)中，diag(·)表示对角矩阵。
该实施例中的系统响应采用Newmark-β法进行数值求解。采用四个质量块的加速度响应为结构模态分析所用的响应信号样本，将Newmark-β法的求解结果以f＝1024Hz重新采样，记录时间为10s(t∈[0,10s]),信号长度Nt＝10240。
本实施例公开的基于多参考时域传递率模型的仅输出结构模态参数辨识方法，包括以下步骤：
步骤1：采集工程结构的响应信号，并将响应信号进行分组。
利用信号采集系统采集得到工程结构的响应信号，将其分成参考输出信号和非参考输出信号两组。该实施例中参考信号选为质量块m1和m2的加速度响应信号，即x1(t)和x2(t)；非参考信号选为质量块m3和m4的加速度响应信号，即x3(t)和x4(t)。因此，
Figure BDA0001387368700000131
Nr＝2，Nt＝10240，No＝4。
步骤2：根据步骤1中的参考输出信号xr与非参考输出信号xl，构造多参考时域传递率模型。
本实施例中的多参考时域传递率模型如式(28)所示：
Figure BDA0001387368700000132
式(28)中，模型阶数n在步骤3中选取；Ai与Bj均为待估参数，可以写成如式(29)所示的矩阵形式：
Θ＝[AT BT]T (29)
式(29)中，A与B分别定义如式(30)所示：
A＝[A0 A1 … An]T，B＝[B0 B1 … Bn]T (30)
步骤3：确定式(28)中多参考时域传递率模型阶数n。
式(28)中多参考时域传递率模型阶数n的确定依据赤池信息量准则(AkaikeInformation Criterion，AIC)。当多参考时域传递率模型阶数n达到最优时，AIC的值取最小。因此，赤池信息量准则AIC值取最小时对应的多参考时域传递率模型阶数即为n值。
本实施例中的AIC值随阶数n的变化曲线如图3所示，因此n＝4。
步骤4：根据步骤2中构造的多参考时域传递率模型及步骤3中确定的模型阶数n，估计其待估系数矩阵Θ。
步骤4中估计其待估系数矩阵Θ优选选择采用最小二乘法或最大似然法。
步骤4.1：根据步骤2中构造的多参考时域传递率模型及步骤3中确定的模型阶数n，采用最小二乘法估计多参考时域传递率模型待估系数矩阵Θ。
根据最小二乘方法，待估系数矩阵Θ中的子矩阵A与B如式(31)所示：
Figure BDA0001387368700000141
式(31)中，
Figure BDA0001387368700000142
A′由式(32)计算：
A′＝(D′)-1b (32)
式(32)中，D′与b分别定义如式(33)所示：
Figure BDA0001387368700000143
矩阵D定义如式(34)所示：
D＝L-SR-1SH (34)
式(34)中的L定义如式(35)所示：
Figure BDA0001387368700000144
计算矩阵S、L和R时，须用到的矩阵Φl和Φr如式(36)所示：
Figure BDA0001387368700000145
步骤4.2：根据步骤2中构造的多参考时域传递率模型及步骤3中确定的模型阶数n，采用最大似然法估计多参考时域传递率模型待估系数矩阵Θ。
根据最大似然法，待估系数矩阵
Figure BDA0001387368700000146
由式(37)求得：
Figure BDA0001387368700000151
式(13)中的
Figure BDA0001387368700000152
与
Figure BDA0001387368700000153
的定义分别如式(38)所示：
Figure BDA0001387368700000154
式(37)中的矩阵X如式(39)所示：
Figure BDA0001387368700000155
根据式(37)求出待估系数矩阵
Figure BDA0001387368700000156
后，结合式(16)并注意A0＝I2×2，即可求出式(3)中完整的待估系数矩阵A与B。
即完成采用最大似然法估计多参考时域传递率模型待估系数矩阵Θ步骤。
步骤5：选取多个不同工况，针对每个工况完成多参考时域传递率模型待估系数矩阵Θ的估计。
选取多个不同的工况，对每一个工况，均执行步骤1、步骤2及步骤4，模型阶数n恒取为4，获得所有不同工况下的待估参数A与B。
本实施例中选取两个不同工况，Nk＝2。第一个工况载荷作用在质量块m1和m3上(记为工况1)，第二个工况载荷作用在质量块m2和m4上(记为工况2)。
工况1和工况2的载荷作用示意图分别如图4(A)和图4(B)所示。作用在系统上的载荷是非白的，其时域信号如图5(A)所示，其功率谱如图5(B)所示，容易看出其功率谱不为定值，存在谐波分量。
工况1下质量块m1的加速度响应的功率谱如图6所示，图中曲线标“×”位置即为载荷中的谐波分量。若采用已有技术中的工作模态分析方法，该分量会被辨识成结构的固有模态，与实际不符，因此已有技术中的工作模态分析方法必须在输入为白噪声的前提下才有效。
在工况1和工况2下，均执行步骤1、步骤2及步骤4，分别得到工况1下的待估系数矩阵A1、B1与工况2下的待估系数矩阵A2、B2。
即完成每个工况下多参考时域传递率模型待估系数矩阵Θ的估计。
步骤6：根据步骤5中得到的所有不同工况下待估系数矩阵Ak与Bk，完成仅输出结构模态参数辨识。
根据步骤5中得到的工况1和工况2下待估系数矩阵A1、B1与A2、B2，获得结构的模态参数，所述模态参数包括模态频率、模态阻尼及最大值归一化的模态向量。
将工况1下得到的参数矩阵A1、B1与工况2下得到的参数矩阵A2、B2组装成如式(40)所示的矩阵：
Figure BDA0001387368700000161
矩阵C′i(i＝0,1,2,3)可由式(41)求得：
Figure BDA0001387368700000162
求解如式(42)所示的广义特征值问题：
Figure BDA0001387368700000163
式(42)中，矩阵Γ、向量
Figure BDA0001387368700000164
分别定义如式(43)所示：
Figure BDA0001387368700000165
根据式(42)求得特征值λm及结构未最大值归一化的模态向量φm之后，可根据式(24)求得结构的模态频率ωm、模态阻尼ξm及最大值归一化的模态向量，即完成基于多参考时域传递率模型的仅输出结构模态参数辨识。
还包括步骤7：将步骤6辨识的结构模态参数应用于工程结构力学领域解决实际工程问题。优选选择应用于结构设计、状态监测、振动控制领域。
由于本实施例步骤4中估计待估系数矩阵Θ优选选择采用最小二乘法或最大似然法，而最小二乘法和最大似然法得到的结果是等价的，
因此在本实施例中仅给出最小二乘法估计待估系数矩阵Θ得到的结果。最终辨识得到的结构模态频率、模态阻尼及最大值归一化的模态向量分别如图7和图8所示。
由图7能够看出，圆点代表的本实施例中的基于多参考时域传递率模型的仅输出结构模态参数辨识方法，相对于三角形代表的已有频域方法辨识出的模态参数精度更高。图7中本实施例中的方法得到的结果紧密聚集在真实值附近，而已有的频域方法得到的结果非常分散，若在实际应用中只做单次辨识，则已有的频域法得到的结果不可信。由图8可以看出，(A)图代表的本实施例中的基于多参考时域传递率模型的仅输出结构模态参数辨识方法，相对于(B)图代表的已有频域方法辨识出的最大值归一化模态向量与真实值吻合地更好。
因此本实施例公开的方法在实际应用时得到的模态参数精度更高、鲁棒性更强、可信度更高，对接下来的结构设计、状态监测、振动控制等工作更有意义。
本实施例公开的基于多参考时域传递率模型的仅输出结构模态参数辨识方法，能处理非白输入条件下的结构模态辨识问题，结构不受载荷形式的限制，拓展应用范围，而且在对结构进行模态分析前，不用已知结构所受载荷中谐波成分的频率范围，与已有技术中的滤波等处理载荷中谐波分量的方法相比，限制条件更少，即使在完全缺乏结构所受激励先验知识的情况下也能进行操作。同时，本实施例将已有的在频域内定义的传递率函数拓展到时域，实际中测量的结构时域输出信号可直接使用，不用通过傅里叶变换变换到频域，因此不用选择窗函数，相对于已有技术中的方法使用更方便，同时能够避免测量的时域数据向频域数据的转换，不仅简单易行，而且去除傅里叶变换中出现的泄露、频率分辨率降低、数据间相互作用加强等误差及不良影响，使得对工程结构的模态分析精度更高。
以上所述的具体描述，对发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例，用于解释本发明，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (5)

1.基于多参考时域传递率的仅输出结构模态参数辨识方法，其特征在于：包括以下步骤，
步骤1：采集工程结构的响应信号，并将响应信号进行分组；
利用信号采集系统采集得到工程结构的响应信号，将响应信号分成两组，其中一组为参考输出信号，记为xr，另外一组为非参考输出信号，记为xl；下标r表示参考信号，l表示非参考信号；
其中，
Figure FDA0002322423750000011
为参考输出信号，
Figure FDA0002322423750000012
表示实数集，Nr表示参考通道数，Nt为采集数据的长度；
Figure FDA0002322423750000013
为非参考输出信号，No为输出通道总数；Nr的取值大于等于2，小于等于
Figure FDA0002322423750000014
步骤2：根据步骤1中的参考输出信号xr与非参考输出信号xl，构造多参考时域传递率模型；
将参考输出信号xr当做模型输入，非参考输出信号xl当做模型输出，构造含有外部激励的自回归时间序列模型，如式(1)所示：
Figure FDA0002322423750000015
式(1)中，n为模型阶数；x(t)表示第t个时刻的信号，
Figure FDA0002322423750000016
为第t个时刻的残差序列；Ai表示第i阶AR系数，Bj表示第j阶外部激励系数；Ai与Bj均为待估参数，写成如式(2)所示的矩阵形式：
Θ＝[AT BT]T (2)
式(2)中，Θ为由待估系数矩阵Ai与Bj构成的新矩阵，T表示矩阵转置运算，其中A与B分别定义如式(3)所示：
A＝[A0 A1…An]T，B＝[B0 B1…Bn]T (3)
式(1)所示的含有外部激励的自回归时间序列模型即为构造的多参考时域传递率模型；
步骤3：确定式(1)中多参考时域传递率模型阶数n；
依据赤池信息量准则AIC(Akaike Information Criterion)确定式(1)中的多参考时域传递率模型阶数n；当多参考时域传递率模型阶数n达到最优时，AIC的值取最小；因此，赤池信息量准则AIC值取最小时对应的多参考时域传递率模型阶数即为n值；
步骤4：根据步骤2中构造的多参考时域传递率模型及步骤3中确定的模型阶数n，估计其待估系数矩阵Θ；
步骤5：选取多个不同工况，针对每个工况完成多参考时域传递率模型待估系数矩阵Θ的估计；
选取多个不同的工况，对每一个工况，均执行步骤1、步骤2及步骤4，步骤3只须运行一次，对所有工况均取相同的模型阶数n，获得所有不同工况下的待估参数Ak与Βk；
下标k表示第k个工况，k＝1,2,…,Nk，Nk为不同工况总数；结构受到的激励幅值、位置或形式不同，均当做不同工况，且要求Nk≥2；为区别不同工况的待估系数矩阵，式(3)改写成式(4)的形式：
Ak＝[Ak,0 Ak,1…Ak,n]T，Bk＝[Bk,0 Bk,1…Bk,n]T (4)
即完成每个工况下多参考时域传递率模型待估系数矩阵Θ的估计；
步骤6：根据步骤5中得到的所有不同工况下待估系数矩阵Ak与Βk，完成仅输出结构模态参数辨识；
所述的仅输出结构模态参数包括模态频率、模态阻尼及最大值归一化的模态向量；
步骤6.1：根据步骤5中得到的所有不同工况下待估系数矩阵Θk，求解如式(5)所示的广义特征值问题：
Figure FDA0002322423750000021
式(5)中，λm为矩阵Γ的第m个特征值，向量
Figure FDA0002322423750000022
表示λm对应的矩阵Γ的特征向量；矩阵Γ、向量
Figure FDA0002322423750000023
分别定义如式(6)所示：
Figure FDA0002322423750000024
式(6)中，φm为结构的第m阶未最大值归一化的模态向量；矩阵Ci′，定义如式(7)所示：
Figure FDA0002322423750000025
式(7)中，矩阵
Figure FDA0002322423750000026
其中j＝0,1,…,n，定义如式(8)所示：
Figure FDA0002322423750000031
式(8)中，矩阵A与Β的第一个下标表示工况数，第二个下标表示步骤2式(1)中的阶数；
式(7)中的矩阵求逆运算要求矩阵Cn为方阵，即要求Nk＝No/(No-Nr)；当Cn不为方阵时，可改用伪逆运算，如式(9)所示：
Figure FDA0002322423750000032
上标“+”表示矩阵伪逆运算；
步骤6.2：根据步骤6.1中求得的特征值λm与未最大值归一化的模态向量φm求出结构的模态频率、模态阻尼及最大值归一化的模态向量；
系统的特征值写成λm＝cm+idm的形式，则结构第m阶模态频率ωm、第m阶模态阻尼ξm与第m阶最大值归一化的模态向量
Figure FDA0002322423750000033
分别定义如式(10)所示：
Figure FDA0002322423750000034
式(10)中，ωm的单位为rad/s，abs(·)表示对向量的每个元素取绝对值，max(·)表示取向量中最大的元素，αm与βm分别定义如式(11)所示：
Figure FDA0002322423750000035
根据式(10)完成基于多参考时域传递率的仅输出结构模态参数辨识；
将步骤6辨识的结构模态参数应用于工程结构力学领域解决实际工程问题包括应用于结构设计、状态监测、振动控制领域；
步骤7.1：将步骤6辨识的结构模态参数应用于结构设计领域，实现针对任意激励下的工程结构进行模态分析，根据步骤6辨识的结构模态参数，得到工程结构的固有频率，以此检测工程结构对环境激励的敏感程度，检测是否符合工程结构标准及隔振的要求，提高结构设计精度与可靠度；
步骤7.2：将步骤6辨识的结构模态参数应用于结构状态监测领域，步骤6辨识的结构模态参数对结构故障具有极高的敏感度，根据步骤6辨识的结构模态参数对结构健康状态进行监测，并检测结构是否出现故障，提高结构健康状态监测的实时性与可靠性；
步骤7.3：将步骤6辨识的结构模态参数应用于结构振动控制领域，根据步骤6辨识的结构模态参数确定结构对环境激励的敏感程度及发生共振的频率范围，降低结构振动幅值，完成结构振动控制。
2.如权利要求1所述的基于多参考时域传递率的仅输出结构模态参数辨识方法，其特征在于：步骤4中估计其待估系数矩阵Θ选择采用最小二乘法或最大似然法。
3.如权利要求2所述的基于多参考时域传递率的仅输出结构模态参数辨识方法，其特征在于：
步骤4.1：根据步骤2中构造的多参考时域传递率模型及步骤3中确定的模型阶数n，采用最小二乘法估计多参考时域传递率模型待估系数矩阵Θ，具体实现方法如下：
步骤4.1.1：根据步骤2中构造的多参考时域传递率模型，写出最小二乘法对应的待估系数矩阵Θ的费用函数；
最小二乘法估计待估系数矩阵Θ的费用函数JLS如式(12)所示：
Figure FDA0002322423750000041
式(12)中，j表示待估系数矩阵Θ的第j列；T表示矩阵的转置运算；Φ为回归矩阵，形式如式(13)所示：
Φ＝[Φl Φr] (13)
式(13)中的Φl与Φr的定义分别如式(14)所示：
Figure FDA0002322423750000042
步骤4.1.2：根据步骤4.1.1中如式(12)所示的费用函数JLS，使费用函数取得最小值的待估参数值即为待估系数矩阵Θ的估计值，求出待估系数矩阵Α与Β；
根据式(12)所示的费用函数JLS，待估系数矩阵Θ根据式(15)求得：
Figure FDA0002322423750000051
式(15)中，
Figure FDA0002322423750000052
(·)-1表示矩阵求逆运算；A′由式(16)计算：
A′＝(D′)-1b (16)
式(16)中，D′与b分别定义如式(17)所示：
Figure FDA0002322423750000053
式(17)中，[·][:,:]表示取矩阵的一部分，矩阵D定义如式(18)所示：
D＝L-SR-1SH (18)
式(18)中的L定义如式(19)所示：
Figure FDA0002322423750000054
即完成采用最小二乘法估计多参考时域传递率模型待估系数矩阵Θ步骤。
4.如权利要求2所述的基于多参考时域传递率的仅输出结构模态参数辨识方法，其特征在于：
步骤4.2：根据步骤2中构造的多参考时域传递率模型及步骤3中确定的模型阶数n，采用最大似然法估计多参考时域传递率模型待估系数矩阵Θ，具体实现方法如下：
步骤4.2.1：根据步骤2中构造的多参考时域传递率模型，写出最大似然法对应的待估系数矩阵Θ的费用函数；
采用最大似然法时，将式(1)中的系数矩阵A0取为单位矩阵；根据对数似然函数形式，最大似然估计法的费用函数如式(20)所示：
Figure FDA0002322423750000055
式(20)中，
Figure FDA0002322423750000056
为残差序列e(t)的方差矩阵，且对任意时刻t均相同；|·|表示行列式，tr(·)表示矩阵的迹；
Figure FDA0002322423750000057
为回归矩阵，形式如式(21)所示：
Figure FDA0002322423750000058
式(21)中的
Figure FDA0002322423750000059
与
Figure FDA00023224237500000510
的定义分别如式(22)所示：
Figure FDA0002322423750000061
式(20)中的矩阵X如式所示：
Figure FDA0002322423750000062
式(20)中最大似然法对应的待估系数矩阵
Figure FDA0002322423750000063
如式(24)所示：
Figure FDA0002322423750000064
步骤4.2.2：根据步骤4.2.1中如式(20)所示的费用函数
Figure FDA0002322423750000065
使费用函数取得最大值的待估参数值即为待估系数矩阵
Figure FDA0002322423750000066
的估计值，求出待估系数矩阵Α与Β；
根据如式(20)所示的费用函数
Figure FDA0002322423750000067
避免最大似然法中的非线性迭代求解，将残差序列e(t)的方差矩阵Σ看成与待估系数矩阵
Figure FDA0002322423750000068
无关，得到的待估系数矩阵与最小二乘法得到的结果是等价的；待估系数矩阵
Figure FDA0002322423750000069
由式(25)求得：
Figure FDA00023224237500000610
式(25)中，(·)-1表示矩阵求逆运算；求出待估系数矩阵
Figure FDA00023224237500000611
后，结合式(24)并注意
Figure FDA00023224237500000612
Figure FDA00023224237500000613
表示No-Nr维单位矩阵，即求出式(3)中完整的待估系数矩阵A与B；
即完成采用最大似然法估计多参考时域传递率模型待估系数矩阵Θ步骤。
5.如权利要求1、2、3或4所述的基于多参考时域传递率的仅输出结构模态参数辨识方法，其特征在于：还包括步骤7：将步骤6辨识的结构模态参数应用于工程结构力学领域解决实际工程问题。
CN201710732229.8A 2017-08-24 2017-08-24 基于多参考时域传递率的仅输出结构模态参数辨识方法 Expired - Fee Related CN107561934B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710732229.8A CN107561934B (zh) 2017-08-24 2017-08-24 基于多参考时域传递率的仅输出结构模态参数辨识方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710732229.8A CN107561934B (zh) 2017-08-24 2017-08-24 基于多参考时域传递率的仅输出结构模态参数辨识方法

Publications (2)

Publication Number Publication Date
CN107561934A CN107561934A (zh) 2018-01-09
CN107561934B true CN107561934B (zh) 2020-04-21

Family

ID=60976743

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710732229.8A Expired - Fee Related CN107561934B (zh) 2017-08-24 2017-08-24 基于多参考时域传递率的仅输出结构模态参数辨识方法

Country Status (1)

Country Link
CN (1) CN107561934B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110688729B (zh) * 2019-08-26 2023-07-14 南京航空航天大学 基于自适应卡尔曼滤波的lstm-idm跟驰特性融合方法、存储介质及设备
CN110749655B (zh) * 2019-10-24 2021-05-07 大连理工大学 一种针对比例阻尼结构的复模态辨识方法
CN117669211B (zh) * 2023-12-06 2024-06-25 南京航空航天大学 基于参数化时域传递率的结构参数辨识及方差计算方法

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN100420944C (zh) * 2005-12-29 2008-09-24 西北工业大学 基于互相关函数幅值向量的随机振动结构损伤的检测方法
CN103207931A (zh) * 2013-03-11 2013-07-17 武汉大学 基于mm和arma算法的次同步振荡模态辨识方法
CN103926538B (zh) * 2014-05-05 2016-10-05 山东大学 基于aic准则的变阶数rc等效电路模型及实现方法
CN104166804B (zh) * 2014-08-20 2018-01-30 中国科学技术大学 一种基于时频域单源点稀疏成分分析的工作模态辨识方法
CN106326530A (zh) * 2016-08-10 2017-01-11 北京理工大学 一种基于右矩阵分式模型的时变结构模态参数辨识方法
CN106354695B (zh) * 2016-08-22 2019-09-17 北京理工大学 一种仅输出线性时变结构模态参数辨识方法

Also Published As

Publication number Publication date
CN107561934A (zh) 2018-01-09

Similar Documents

Publication Publication Date Title
CN107561934B (zh) 基于多参考时域传递率的仅输出结构模态参数辨识方法
CN112329855B (zh) 基于自适应字典的欠定工作模态参数识别方法及检测方法
Chen et al. Simultaneous identification of structural parameters and input time history from output-only measurements
CN112507606B (zh) 基于rbf网络的欠定工作模态参数识别方法及检测方法
CN107844627A (zh) 一种仅输出时变结构模态参数贝叶斯估计方法
Hu et al. Model order determination and noise removal for modal parameter estimation
CN101706355A (zh) 基于NExT/ARMA的结构响应分析方法
CN104133950A (zh) 一种悬臂梁运行模态分析实验方法及装置
CN114925526B (zh) 一种结合多工况响应的结构模态参数识别方法
CN107168063A (zh) 基于集成变量选择型偏最小二乘回归的软测量方法
CN109446552B (zh) 多轴相关随机激励下结构疲劳寿命时域计算方法
CN101587007A (zh) 识别柔性桥梁结构动力参数的惟输出小波基分析方法
CN109753690A (zh) 基于计算流体力学的非线性非定常气动力降阶方法
CN109254321B (zh) 一种地震激励下快速贝叶斯模态参数识别方法
CN110657934A (zh) 一种电动振动台在线修正迭代控制方法
CN111308979A (zh) 基于多率迟延状态空间模型的辨识方法及u控制系统
CN106326530A (zh) 一种基于右矩阵分式模型的时变结构模态参数辨识方法
CN106446502B (zh) 带遗忘因子的特征向量递推的时变工作模态在线识别方法
CN110414150B (zh) 一种桥梁时变系统的张量子空间连续系统识别方法
CN111274630A (zh) 一种用于工程结构柔度识别的物理模态提取方法
CN106446503B (zh) 遗忘自协方差矩阵递推主元的时变工作模态识别方法
CN117892046A (zh) 一种基于动态模式分解的时变模态仅输出递推辨识方法
CN106980722A (zh) 一种脉冲响应中谐波成分的检测和去除方法
CN111428342B (zh) 一种基于频域谱分解的随机动载荷识别方法
CN115630273A (zh) 增量式时变结构工作模态参数实时识别方法及系统

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
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20200421