CN107561934A - 基于多参考时域传递率的仅输出结构模态参数辨识方法 - Google Patents
基于多参考时域传递率的仅输出结构模态参数辨识方法 Download PDFInfo
- Publication number
- CN107561934A CN107561934A CN201710732229.8A CN201710732229A CN107561934A CN 107561934 A CN107561934 A CN 107561934A CN 201710732229 A CN201710732229 A CN 201710732229A CN 107561934 A CN107561934 A CN 107561934A
- Authority
- CN
- China
- Prior art keywords
- mrow
- msub
- formula
- mtd
- time domain
- 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
Links
Landscapes
- Testing Of Devices, Machine Parts, Or Other Structures Thereof (AREA)
- Other Investigation Or Analysis Of Materials By Electrical Means (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表示非参考信号。
其中,为参考输出信号,表示实数集,Nr表示参考通道数,Nt为采集数据的长度;为非参考输出信号,No为输出通道总数。Nr的取值大于等于2,小于等于
步骤2:根据步骤1中的参考输出信号xr与非参考输出信号xl,构造多参考时域传递率模型。
将参考输出信号xr当做模型输入,非参考输出信号xl当做模型输出,构造含有外部激励的自回归时间序列模型,如式(1)所示:
式(1)中,n为模型阶数;x(t)表示第t个时刻的信号,为第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)所示:
式(4)中,j表示待估系数矩阵Θ的第j列;T表示矩阵的转置运算;Φ为回归矩阵,形式如式(5)所示:
Φ=[Φl Φr] (5)
式(5)中的Φl与Φr的定义分别如式(6)所示:
步骤4.1.2:根据步骤4.1.1中如式(4)所示的费用函数JLS,使费用函数取得最小值的待估参数值即为待估系数矩阵Θ的估计值,求出待估系数矩阵A与B。
根据式(4)所示的费用函数JLS,待估系数矩阵Θ根据式(7)求得:
式(7)中,(·)-1表示矩阵求逆运算。A′由式(8)计算:
A′=(D′)-1b (8)
式(8)中,D′与b分别定义如式(9)所示:
式(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)所示:
式(12)中,为残差序列e(t)的方差矩阵,且对任意时刻t均相同;|·|表示行列式,tr(·)表示矩阵的迹。为回归矩阵,形式如式(13)所示:
式(13)中的与的定义分别如式(14)所示:
式(12)中的矩阵X如式所示:
式(12)中最大似然法对应的待估系数矩阵如式(16)所示:
步骤4.2.2:根据步骤4.2.1中如式(12)所示的费用函数使费用函数取得最大值的待估参数值即为待估系数矩阵的估计值,求出待估系数矩阵A与B。
根据如式(12)所示的费用函数避免最大似然法中的非线性迭代求解,将残差序列e(t)的方差矩阵∑看成与待估系数矩阵无关,得到的待估系数矩阵与最小二乘法得到的结果是等价的。待估系数矩阵由式(17)求得:
式(17)中,(·)-1表示矩阵求逆运算。求出待估系数矩阵后,结合式(16)并注意表示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)所示的广义特征值问题:
式(19)中,λm为矩阵Γ的第m个特征值,向量表示λm对应的矩阵Γ的特征向量;矩阵Γ、向量分别定义如式(20)所示:
式(20)中,φm为结构的第m阶未最大值归一化的模态向量;矩阵C′i(其中i=0,1,…,n-1)定义如式(21)所示:
式(21)中,矩阵(其中j=0,1,…,n)定义如式(22)所示:
式(22)中,矩阵A与B的第一个下标表示工况数,第二个下标表示步骤2式(1)中的阶数。
式(21)中的矩阵求逆运算要求矩阵Cn为方阵,即要求Nk=No/(No-Nr)。当Cn不为方阵时,可改用伪逆运算,如式(23)所示:
上标“+”表示矩阵伪逆运算。
步骤6.2:根据步骤6.1中求得的特征值λm与未最大值归一化的模态向量φm求出结构的模态频率、模态阻尼及最大值归一化的模态向量。
系统的特征值写成λm=cm+idm的形式,则结构第m阶模态频率ωm、第m阶模态阻尼ξm与第m阶最大值归一化的模态向量分别定义如式(24)所示:
式(24)中,ωm的单位为rad/s,abs(·)表示对向量的每个元素取绝对值,max(·)表示取向量中最大的元素,αm与βm分别定义如式(25)所示:
根据式(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四自由度时变系统参数
图2所示的四自由度系统的动力学控制方程如式(26)所示:
式(26)中,M、C、K分别为质量矩阵、阻尼矩阵和刚度矩阵,x为系统的响应向量,f(t)为外界激励。M、C、K如式(27)所示:
式(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)。因此,Nr=2,Nt=10240,No=4。
步骤2:根据步骤1中的参考输出信号xr与非参考输出信号xl,构造多参考时域传递率模型。
本实施例中的多参考时域传递率模型如式(28)所示:
式(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)所示:
式(31)中,A′由式(32)计算:
A′=(D′)-1b (32)
式(32)中,D′与b分别定义如式(33)所示:
矩阵D定义如式(34)所示:
D=L-SR-1SH (34)
式(34)中的L定义如式(35)所示:
计算矩阵S、L和R时,须用到的矩阵Φl和Φr如式(36)所示:
步骤4.2:根据步骤2中构造的多参考时域传递率模型及步骤3中确定的模型阶数n,采用最大似然法估计多参考时域传递率模型待估系数矩阵Θ。
根据最大似然法,待估系数矩阵由式(37)求得:
式(13)中的与的定义分别如式(38)所示:
式(37)中的矩阵X如式(39)所示:
根据式(37)求出待估系数矩阵后,结合式(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)所示的矩阵:
矩阵C′i(i=0,1,2,3)可由式(41)求得:
求解如式(42)所示的广义特征值问题:
式(42)中,矩阵Γ、向量分别定义如式(43)所示:
根据式(42)求得特征值λm及结构未最大值归一化的模态向量φm之后,可根据式(24)求得结构的模态频率ωm、模态阻尼ξm及最大值归一化的模态向量,即完成基于多参考时域传递率模型的仅输出结构模态参数辨识。
还包括步骤7:将步骤6辨识的结构模态参数应用于工程结构力学领域解决实际工程问题。优选选择应用于结构设计、状态监测、振动控制领域。
由于本实施例步骤4中估计待估系数矩阵Θ优选选择采用最小二乘法或最大似然法,而最小二乘法和最大似然法得到的结果是等价的,
因此在本实施例中仅给出最小二乘法估计待估系数矩阵Θ得到的结果。最终辨识得到的结构模态频率、模态阻尼及最大值归一化的模态向量分别如图7和图8所示。
由图7能够看出,圆点代表的本实施例中的基于多参考时域传递率模型的仅输出结构模态参数辨识方法,相对于三角形代表的已有频域方法辨识出的模态参数精度更高。图7中本实施例中的方法得到的结果紧密聚集在真实值附近,而已有的频域方法得到的结果非常分散,若在实际应用中只做单次辨识,则已有的频域法得到的结果不可信。由图8可以看出,(A)图代表的本实施例中的基于多参考时域传递率模型的仅输出结构模态参数辨识方法,相对于(B)图代表的已有频域方法辨识出的最大值归一化模态向量与真实值吻合地更好。
因此本实施例公开的方法在实际应用时得到的模态参数精度更高、鲁棒性更强、可信度更高,对接下来的结构设计、状态监测、振动控制等工作更有意义。
本实施例公开的基于多参考时域传递率模型的仅输出结构模态参数辨识方法,能处理非白输入条件下的结构模态辨识问题,结构不受载荷形式的限制,拓展应用范围,而且在对结构进行模态分析前,不用已知结构所受载荷中谐波成分的频率范围,与已有技术中的滤波等处理载荷中谐波分量的方法相比,限制条件更少,即使在完全缺乏结构所受激励先验知识的情况下也能进行操作。同时,本实施例将已有的在频域内定义的传递率函数拓展到时域,实际中测量的结构时域输出信号可直接使用,不用通过傅里叶变换变换到频域,因此不用选择窗函数,相对于已有技术中的方法使用更方便,同时能够避免测量的时域数据向频域数据的转换,不仅简单易行,而且去除傅里叶变换中出现的泄露、频率分辨率降低、数据间相互作用加强等误差及不良影响,使得对工程结构的模态分析精度更高。
以上所述的具体描述,对发明的目的、技术方案和有益效果进行了进一步详细说明,所应理解的是,以上所述仅为本发明的具体实施例,用于解释本发明,并不用于限定本发明的保护范围,凡在本发明的精神和原则之内,所做的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。
Claims (9)
1.基于多参考时域传递率的仅输出结构模态参数辨识方法,其特征在于:包括以下步骤,
步骤1:采集工程结构的响应信号,并将响应信号进行分组;
利用信号采集系统采集得到工程结构的响应信号,将响应信号分成两组,其中一组为参考输出信号,记为xr,另外一组为非参考输出信号,记为xl;下标r表示参考信号,l表示非参考信号;
其中,为参考输出信号,表示实数集,Nr表示参考通道数,Nt为采集数据的长度;为非参考输出信号,No为输出通道总数;Nr的取值大于等于2,小于等于
步骤2:根据步骤1中的参考输出信号xr与非参考输出信号xl,构造多参考时域传递率模型;
将参考输出信号xr当做模型输入,非参考输出信号xl当做模型输出,构造含有外部激励的自回归时间序列模型,如式(1)所示:
<mrow>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>0</mn>
</mrow>
<mi>n</mi>
</munderover>
<msub>
<mi>A</mi>
<mi>i</mi>
</msub>
<msub>
<mi>x</mi>
<mi>l</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>-</mo>
<mi>i</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>j</mi>
<mo>=</mo>
<mn>0</mn>
</mrow>
<mi>n</mi>
</munderover>
<msub>
<mi>B</mi>
<mi>j</mi>
</msub>
<msub>
<mi>x</mi>
<mi>r</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>-</mo>
<mi>j</mi>
<mo>)</mo>
</mrow>
<mo>+</mo>
<mi>e</mi>
<mrow>
<mo>(</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
</mrow>
式(1)中,n为模型阶数;x(t)表示第t个时刻的信号,为第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,估计其待估系数矩阵Θ;
步骤5:选取多个不同工况,针对每个工况完成多参考时域传递率模型待估系数矩阵Θ的估计;
选取多个不同的工况,对每一个工况,均执行步骤1、步骤2及步骤4,步骤3只须运行一次,对所有工况均取相同的模型阶数n,获得所有不同工况下的待估参数Ak与Bk;
下标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与BV,完成仅输出结构模态参数辨识;
所述的仅输出结构模态参数包括模态频率、模态阻尼及最大值归一化的模态向量;
步骤6.1:根据步骤5中得到的所有不同工况下待估系数矩阵Θk,求解如式(5)所示的广义特征值问题:
式(5)中,λm为矩阵Γ的第m个特征值,向量表示λm对应的矩阵Γ的特征向量;矩阵Γ、向量分别定义如式(6)所示:
式(6)中,φm为结构的第m阶未最大值归一化的模态向量;矩阵C′i(其中i=0,1,…,n-1)定义如式(7)所示:
<mrow>
<msubsup>
<mi>C</mi>
<mrow>
<mi>i</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
<mo>&prime;</mo>
</msubsup>
<mo>=</mo>
<msubsup>
<mi>C</mi>
<mi>n</mi>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msubsup>
<msub>
<mi>C</mi>
<mrow>
<mi>i</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
<mo>,</mo>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
<mo>,</mo>
<mn>2</mn>
<mo>,</mo>
<mo>...</mo>
<mo>,</mo>
<mi>n</mi>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>7</mn>
<mo>)</mo>
</mrow>
</mrow>
式(7)中,矩阵(其中j=0,1,…,n)定义如式(8)所示:
<mrow>
<msub>
<mi>C</mi>
<mi>j</mi>
</msub>
<mo>=</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<msub>
<mi>B</mi>
<mrow>
<mn>1</mn>
<mo>,</mo>
<mi>j</mi>
</mrow>
</msub>
</mtd>
<mtd>
<mrow>
<mo>-</mo>
<msub>
<mi>A</mi>
<mrow>
<mn>1</mn>
<mo>,</mo>
<mi>j</mi>
</mrow>
</msub>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>B</mi>
<mrow>
<mn>2</mn>
<mo>,</mo>
<mi>j</mi>
</mrow>
</msub>
</mtd>
<mtd>
<mrow>
<mo>-</mo>
<msub>
<mi>A</mi>
<mrow>
<mn>2</mn>
<mo>,</mo>
<mi>j</mi>
</mrow>
</msub>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<mo>.</mo>
</mtd>
<mtd>
<mo>.</mo>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>B</mi>
<mrow>
<msub>
<mi>N</mi>
<mi>k</mi>
</msub>
<mo>,</mo>
<mi>j</mi>
</mrow>
</msub>
</mtd>
<mtd>
<mrow>
<mo>-</mo>
<msub>
<mi>A</mi>
<mrow>
<msub>
<mi>N</mi>
<mi>k</mi>
</msub>
<mo>,</mo>
<mi>j</mi>
</mrow>
</msub>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>8</mn>
<mo>)</mo>
</mrow>
</mrow>
式(8)中,矩阵A与B的第一个下标表示工况数,第二个下标表示步骤2式(1)中的阶数;
式(7)中的矩阵求逆运算要求矩阵Cn为方阵,即要求Nk=No/(No-Nr);当Cn不为方阵时,可改用伪逆运算,如式(9)所示:
<mrow>
<msubsup>
<mi>C</mi>
<mrow>
<mi>i</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
<mo>&prime;</mo>
</msubsup>
<mo>=</mo>
<msubsup>
<mi>C</mi>
<mi>n</mi>
<mo>+</mo>
</msubsup>
<msub>
<mi>C</mi>
<mrow>
<mi>i</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msub>
<mo>,</mo>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
<mo>,</mo>
<mn>2</mn>
<mo>,</mo>
<mo>...</mo>
<mo>,</mo>
<mi>n</mi>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>9</mn>
<mo>)</mo>
</mrow>
</mrow>
上标“+”表示矩阵伪逆运算;
步骤6.2:根据步骤6.1中求得的特征值λm与未最大值归一化的模态向量φm求出结构的模态频率、模态阻尼及最大值归一化的模态向量;
系统的特征值写成λm=cm+idm的形式,则结构第m阶模态频率ωm、第m阶模态阻尼ξm与第m阶最大值归一化的模态向量分别定义如式(10)所示:
<mrow>
<msub>
<mi>&omega;</mi>
<mi>m</mi>
</msub>
<mo>=</mo>
<msqrt>
<mrow>
<mo>(</mo>
<msubsup>
<mi>&alpha;</mi>
<mi>m</mi>
<mn>2</mn>
</msubsup>
<mo>+</mo>
<msubsup>
<mi>&beta;</mi>
<mi>m</mi>
<mn>2</mn>
</msubsup>
<mo>)</mo>
</mrow>
</msqrt>
<mo>,</mo>
<msub>
<mi>&xi;</mi>
<mi>m</mi>
</msub>
<mo>=</mo>
<mfrac>
<msub>
<mi>&alpha;</mi>
<mi>m</mi>
</msub>
<msub>
<mi>&omega;</mi>
<mi>m</mi>
</msub>
</mfrac>
<mo>,</mo>
<msub>
<mover>
<mi>&phi;</mi>
<mo>&OverBar;</mo>
</mover>
<mi>m</mi>
</msub>
<mo>=</mo>
<mfrac>
<msub>
<mi>&phi;</mi>
<mi>m</mi>
</msub>
<mrow>
<mi>m</mi>
<mi>a</mi>
<mi>x</mi>
<mrow>
<mo>(</mo>
<mi>a</mi>
<mi>b</mi>
<mi>s</mi>
<mo>(</mo>
<msub>
<mi>&phi;</mi>
<mi>m</mi>
</msub>
<mo>)</mo>
<mo>)</mo>
</mrow>
</mrow>
</mfrac>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>10</mn>
<mo>)</mo>
</mrow>
</mrow>
式(10)中,ωm的单位为rad/s,abs(·)表示对向量的每个元素取绝对值,max(·)表示取向量中最大的元素,αm与βm分别定义如式(11)所示:
<mrow>
<msub>
<mi>&alpha;</mi>
<mi>m</mi>
</msub>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<mrow>
<mn>2</mn>
<mi>&Delta;</mi>
<mi>t</mi>
</mrow>
</mfrac>
<mi>l</mi>
<mi>n</mi>
<mrow>
<mo>(</mo>
<msubsup>
<mi>c</mi>
<mi>m</mi>
<mn>2</mn>
</msubsup>
<mo>+</mo>
<msubsup>
<mi>d</mi>
<mi>m</mi>
<mn>2</mn>
</msubsup>
<mo>)</mo>
</mrow>
<mo>,</mo>
<msub>
<mi>&beta;</mi>
<mi>m</mi>
</msub>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<mrow>
<mi>&Delta;</mi>
<mi>t</mi>
</mrow>
</mfrac>
<msup>
<mi>tan</mi>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
<mfrac>
<msub>
<mi>d</mi>
<mi>m</mi>
</msub>
<msub>
<mi>c</mi>
<mi>m</mi>
</msub>
</mfrac>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>11</mn>
<mo>)</mo>
</mrow>
</mrow>
根据式(10)完成基于多参考时域传递率的仅输出结构模态参数辨识。
2.如权利要求1所述的基于多参考时域传递率的仅输出结构模态参数辨识方法,其特征在于:步骤4中估计其待估系数矩阵Θ选择采用最小二乘法或最大似然法。
3.如权利要求2所述的基于多参考时域传递率的仅输出结构模态参数辨识方法,其特征在于:
步骤4.1:根据步骤2中构造的多参考时域传递率模型及步骤3中确定的模型阶数n,采用最小二乘法估计多参考时域传递率模型待估系数矩阵Θ,具体实现方法如下:
步骤4.1.1:根据步骤2中构造的多参考时域传递率模型,写出最小二乘法对应的待估系数矩阵Θ的费用函数;
最小二乘法估计待估系数矩阵Θ的费用函数JLS如式(12)所示:
<mrow>
<msup>
<mi>J</mi>
<mrow>
<mi>L</mi>
<mi>S</mi>
</mrow>
</msup>
<mrow>
<mo>(</mo>
<mi>&Theta;</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>j</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mrow>
<msub>
<mi>N</mi>
<mi>o</mi>
</msub>
<mo>-</mo>
<msub>
<mi>N</mi>
<mi>r</mi>
</msub>
</mrow>
</munderover>
<msubsup>
<mi>&Theta;</mi>
<mi>j</mi>
<mi>T</mi>
</msubsup>
<msup>
<mi>&Phi;</mi>
<mi>T</mi>
</msup>
<msub>
<mi>&Phi;&Theta;</mi>
<mi>j</mi>
</msub>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>12</mn>
<mo>)</mo>
</mrow>
</mrow>
式(12)中,j表示待估系数矩阵Θ的第j列;T表示矩阵的转置运算;Φ为回归矩阵,形式如式(13)所示:
Φ=[Φl Φr] (13)
式(13)中的Φl与Φr的定义分别如式(14)所示:
步骤4.1.2:根据步骤4.1.1中如式(12)所示的费用函数JLS,使费用函数取得最小值的待估参数值即为待估系数矩阵Θ的估计值,求出待估系数矩阵A与B;
根据式(12)所示的费用函数JLS,待估系数矩阵Θ根据式(15)求得:
<mrow>
<mtable>
<mtr>
<mtd>
<mrow>
<mi>A</mi>
<mo>=</mo>
<mo>&lsqb;</mo>
<msup>
<mi>A</mi>
<mo>&prime;</mo>
</msup>
<mo>;</mo>
<msub>
<mi>I</mi>
<mrow>
<mo>(</mo>
<msub>
<mi>N</mi>
<mi>o</mi>
</msub>
<mo>-</mo>
<msub>
<mi>N</mi>
<mi>r</mi>
</msub>
<mo>)</mo>
<mo>&times;</mo>
<mo>(</mo>
<msub>
<mi>N</mi>
<mi>o</mi>
</msub>
<mo>-</mo>
<msub>
<mi>N</mi>
<mi>r</mi>
</msub>
<mo>)</mo>
</mrow>
</msub>
<mo>&rsqb;</mo>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mi>B</mi>
<mo>=</mo>
<mo>-</mo>
<msup>
<mi>R</mi>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
<msup>
<mi>S</mi>
<mi>H</mi>
</msup>
<mi>A</mi>
</mrow>
</mtd>
</mtr>
</mtable>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>15</mn>
<mo>)</mo>
</mrow>
</mrow>
式(15)中,(·)-1表示矩阵求逆运算;A′由式(16)计算:
A′=(D′)-1b (16)
式(16)中,D′与b分别定义如式(17)所示:
<mrow>
<mtable>
<mtr>
<mtd>
<mrow>
<msup>
<mi>D</mi>
<mo>&prime;</mo>
</msup>
<mo>=</mo>
<msub>
<mrow>
<mo>&lsqb;</mo>
<mi>D</mi>
<mo>&rsqb;</mo>
</mrow>
<mrow>
<mo>&lsqb;</mo>
<mn>1</mn>
<mo>:</mo>
<mi>n</mi>
<mo>&times;</mo>
<mrow>
<mo>(</mo>
<msub>
<mi>N</mi>
<mi>o</mi>
</msub>
<mo>-</mo>
<msub>
<mi>N</mi>
<mi>r</mi>
</msub>
<mo>)</mo>
</mrow>
<mo>,</mo>
<mn>1</mn>
<mo>:</mo>
<mi>n</mi>
<mo>&times;</mo>
<mrow>
<mo>(</mo>
<msub>
<mi>N</mi>
<mi>o</mi>
</msub>
<mo>-</mo>
<msub>
<mi>N</mi>
<mi>r</mi>
</msub>
<mo>)</mo>
</mrow>
<mo>&rsqb;</mo>
</mrow>
</msub>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mi>b</mi>
<mo>=</mo>
<msub>
<mrow>
<mo>&lsqb;</mo>
<mi>D</mi>
<mo>&rsqb;</mo>
</mrow>
<mrow>
<mo>&lsqb;</mo>
<mn>1</mn>
<mo>:</mo>
<mi>n</mi>
<mo>&times;</mo>
<mrow>
<mo>(</mo>
<msub>
<mi>N</mi>
<mi>o</mi>
</msub>
<mo>-</mo>
<msub>
<mi>N</mi>
<mi>r</mi>
</msub>
<mo>)</mo>
</mrow>
<mo>,</mo>
<mrow>
<mo>(</mo>
<mi>n</mi>
<mo>&times;</mo>
<mo>(</mo>
<mrow>
<msub>
<mi>N</mi>
<mi>o</mi>
</msub>
<mo>-</mo>
<msub>
<mi>N</mi>
<mi>r</mi>
</msub>
</mrow>
<mo>)</mo>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
<mo>:</mo>
<mrow>
<mo>(</mo>
<mi>n</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
<mo>&times;</mo>
<mrow>
<mo>(</mo>
<msub>
<mi>N</mi>
<mi>o</mi>
</msub>
<mo>-</mo>
<msub>
<mi>N</mi>
<mi>r</mi>
</msub>
<mo>)</mo>
</mrow>
<mo>&rsqb;</mo>
</mrow>
</msub>
</mrow>
</mtd>
</mtr>
</mtable>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>17</mn>
<mo>)</mo>
</mrow>
</mrow>
式(17)中,[·][:,:]表示取矩阵的一部分,例如[D][a:b,c:d]表示取矩阵D中行号为a到b、列号为c到d的部分;矩阵D定义如式(18)所示:
D=L-SR-1SH (18)
式(18)中的L定义如式(19)所示:
L=Φl TΦl (19)
即完成采用最小二乘法估计多参考时域传递率模型待估系数矩阵Θ步骤。
4.如权利要求2所述的基于多参考时域传递率的仅输出结构模态参数辨识方法,其特征在于:
步骤4.2:根据步骤2中构造的多参考时域传递率模型及步骤3中确定的模型阶数n,采用最大似然法估计多参考时域传递率模型待估系数矩阵Θ,具体实现方法如下:
步骤4.2.1:根据步骤2中构造的多参考时域传递率模型,写出最大似然法对应的待估系数矩阵Θ的费用函数;
采用最大似然法时,将式(1)中的系数矩阵A0取为单位矩阵;根据对数似然函数形式,最大似然估计法的费用函数如式(20)所示:
<mrow>
<msup>
<mi>J</mi>
<mrow>
<mi>M</mi>
<mi>L</mi>
</mrow>
</msup>
<mrow>
<mo>(</mo>
<mover>
<mi>&Theta;</mi>
<mo>&OverBar;</mo>
</mover>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mo>-</mo>
<mfrac>
<mrow>
<msub>
<mi>N</mi>
<mi>t</mi>
</msub>
<mo>-</mo>
<mi>n</mi>
</mrow>
<mn>2</mn>
</mfrac>
<mi>l</mi>
<mi>n</mi>
<mrow>
<mo>(</mo>
<mo>|</mo>
<mo>&Sigma;</mo>
<mo>|</mo>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mfrac>
<mn>1</mn>
<mn>2</mn>
</mfrac>
<mi>t</mi>
<mi>r</mi>
<mrow>
<mo>(</mo>
<msup>
<mo>&Sigma;</mo>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
<msup>
<mrow>
<mo>(</mo>
<mrow>
<mi>X</mi>
<mo>-</mo>
<mover>
<mi>&Phi;</mi>
<mo>&OverBar;</mo>
</mover>
<mover>
<mi>&Theta;</mi>
<mo>&OverBar;</mo>
</mover>
</mrow>
<mo>)</mo>
</mrow>
<mi>T</mi>
</msup>
<mo>(</mo>
<mrow>
<mi>X</mi>
<mo>-</mo>
<mover>
<mi>&Phi;</mi>
<mo>&OverBar;</mo>
</mover>
<mover>
<mi>&Theta;</mi>
<mo>&OverBar;</mo>
</mover>
</mrow>
<mo>)</mo>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>20</mn>
<mo>)</mo>
</mrow>
</mrow>
式(20)中,为残差序列e(t)的方差矩阵,且对任意时刻t均相同;|·|表示行列式,tr(·)表示矩阵的迹;为回归矩阵,形式如式(21)所示:
<mrow>
<mover>
<mi>&Phi;</mi>
<mo>&OverBar;</mo>
</mover>
<mo>=</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<msub>
<mover>
<mi>&Phi;</mi>
<mo>&OverBar;</mo>
</mover>
<mi>l</mi>
</msub>
</mtd>
<mtd>
<msub>
<mover>
<mi>&Phi;</mi>
<mo>&OverBar;</mo>
</mover>
<mi>r</mi>
</msub>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>21</mn>
<mo>)</mo>
</mrow>
</mrow>
式(21)中的与的定义分别如式(22)所示:
式(20)中的矩阵X如式所示:
<mrow>
<mi>X</mi>
<mo>=</mo>
<msubsup>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mrow>
<msub>
<mi>x</mi>
<mi>l</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>n</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
<mtd>
<mrow>
<msub>
<mi>x</mi>
<mi>l</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>n</mi>
<mo>+</mo>
<mn>2</mn>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
<mtd>
<mn>...</mn>
</mtd>
<mtd>
<mrow>
<msub>
<mi>x</mi>
<mi>l</mi>
</msub>
<mrow>
<mo>(</mo>
<msub>
<mi>N</mi>
<mi>t</mi>
</msub>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
<mrow>
<mo>(</mo>
<msub>
<mi>N</mi>
<mi>t</mi>
</msub>
<mo>-</mo>
<mi>n</mi>
<mo>)</mo>
<mo>&times;</mo>
<mo>(</mo>
<msub>
<mi>N</mi>
<mi>o</mi>
</msub>
<mo>-</mo>
<msub>
<mi>N</mi>
<mi>r</mi>
</msub>
<mo>)</mo>
</mrow>
<mi>T</mi>
</msubsup>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>23</mn>
<mo>)</mo>
</mrow>
</mrow>
式(20)中最大似然法对应的待估系数矩阵如式(24)所示:
<mrow>
<mover>
<mi>&Theta;</mi>
<mo>&OverBar;</mo>
</mover>
<mo>=</mo>
<msup>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<msub>
<mi>A</mi>
<mn>1</mn>
</msub>
</mtd>
<mtd>
<mn>...</mn>
</mtd>
<mtd>
<msub>
<mi>A</mi>
<mi>n</mi>
</msub>
</mtd>
<mtd>
<msub>
<mi>B</mi>
<mn>0</mn>
</msub>
</mtd>
<mtd>
<msub>
<mi>B</mi>
<mn>1</mn>
</msub>
</mtd>
<mtd>
<mn>...</mn>
</mtd>
<mtd>
<msub>
<mi>B</mi>
<mi>n</mi>
</msub>
</mtd>
</mtr>
</mtable>
</mfenced>
<mi>T</mi>
</msup>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>24</mn>
<mo>)</mo>
</mrow>
</mrow>
步骤4.2.2:根据步骤4.2.1中如式(20)所示的费用函数使费用函数取得最大值的待估参数值即为待估系数矩阵的估计值,求出待估系数矩阵A与B;
根据如式(20)所示的费用函数避免最大似然法中的非线性迭代求解,将残差序列e(t)的方差矩阵∑看成与待估系数矩阵无关,得到的待估系数矩阵与最小二乘法得到的结果是等价的;待估系数矩阵由式(25)求得:
<mrow>
<mover>
<mi>&Theta;</mi>
<mo>&OverBar;</mo>
</mover>
<mo>=</mo>
<msup>
<mrow>
<mo>(</mo>
<msup>
<mover>
<mi>&Phi;</mi>
<mo>&OverBar;</mo>
</mover>
<mi>T</mi>
</msup>
<mover>
<mi>&Phi;</mi>
<mo>&OverBar;</mo>
</mover>
<mo>)</mo>
</mrow>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
<msup>
<mover>
<mi>&Phi;</mi>
<mo>&OverBar;</mo>
</mover>
<mi>T</mi>
</msup>
<mi>X</mi>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>25</mn>
<mo>)</mo>
</mrow>
</mrow>
式(25)中,(·)-1表示矩阵求逆运算;求出待估系数矩阵后,结合式(24)并注意表示No-Nr维单位矩阵,即求出式(3)中完整的待估系数矩阵A与B;
即完成采用最大似然法估计多参考时域传递率模型待估系数矩阵Θ步骤。
5.如权利要求1、2、3或4所述的基于多参考时域传递率的仅输出结构模态参数辨识方法,其特征在于:还包括步骤7:将步骤6辨识的结构模态参数应用于工程结构力学领域解决实际工程问题。
6.如权利要求5所述的基于多参考时域传递率的仅输出结构模态参数辨识方法,其特征在于:将步骤6辨识的结构模态参数应用于工程结构力学领域解决实际工程问题包括应用于结构设计、状态监测、振动控制领域;
步骤7.1:将步骤6辨识的结构模态参数结构设计领域,实现针对任意激励下的工程结构进行模态分析,根据步骤6辨识的结构模态参数,得到工程结构的固有频率,以此检测工程结构对环境激励的敏感程度,检测是否符合工程结构标准及隔振等要求,提高结构设计精度与可靠度;
步骤7.2:将步骤6辨识的结构模态参数结构状态监测领域,步骤6辨识的结构模态参数对结构故障具有极高的敏感度,根据步骤6辨识的结构模态参数对结构健康状态进行监测,并检测结构是否出现故障,提高结构健康状态监测的实时性与可靠性;
步骤7.3:将步骤6辨识的结构模态参数结构振动控制领域,根据步骤6辨识的结构模态参数确定结构对环境激励的敏感程度及发生共振的频率范围,降低结构振动幅值,完成结构振动控制。
7.基于多参考时域传递率的仅输出结构模态参数辨识方法,其特征在于:首先利用信号采集系统采集得到工程结构的响应信号,将响应信号分成参考输出信号和非参考输出信号两组;根据参考输出信号与非参考输出信号,构造多参考时域传递率模型;根据赤池信息量准则确定多参考时域传递率模型阶数n,完成多参考时域传递率模型的构造;根据构造的多参考时域传递率模型,估计多参考时域传递率模型系数矩阵;选取多个不同的工况,结构受到的激励幅值、位置或形式不同,均当做不同工况;根据构建的多参考时域传递率模型获得所有不同工况下的待估系数矩阵;根据所有不同工况下估计出的模型系数矩阵,构造并求解广义特征值问题,根据得到的特征值与特征向量即求得结构的模态参数;最后将结构模态参数应用于工程结构设计、工程结构状态监测、工程结构振动控制等领域。
8.如权利要求7所述的,其特征在于:所述的仅输出结构模态参数包括模态频率、模态阻尼及最大值归一化的模态向量。
9.如权利要求8所述的,其特征在于:将得到的结构模态参数应用于工程结构力学领域解决实际工程问题;针对任意激励下的工程结构进行模态分析,检测是否符合工程结构标准及隔振等要求,提高结构设计精度与可靠度;结构模态参数对结构故障具有极高的敏感度,将其用于结构的健康状态监测,提高结构健康状态监测的实时性与可靠性;结构模态参数用于结构振动控制领域,根据结构模态参数确定结构对环境激励的敏感程度及发生共振的频率范围,降低结构振动幅值,完成结构振动控制。
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 true CN107561934A (zh) | 2018-01-09 |
CN107561934B CN107561934B (zh) | 2020-04-21 |
Family
ID=60976743
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201710732229.8A Active CN107561934B (zh) | 2017-08-24 | 2017-08-24 | 基于多参考时域传递率的仅输出结构模态参数辨识方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN107561934B (zh) |
Cited By (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN110688729A (zh) * | 2019-08-26 | 2020-01-14 | 南京航空航天大学 | 基于自适应卡尔曼滤波的lstm-idm跟驰特性融合方法、存储介质及设备 |
WO2021077467A1 (en) * | 2019-10-24 | 2021-04-29 | Dalian University Of Technology | Method of complex modal identification for structure with proportional damping |
CN117669211A (zh) * | 2023-12-06 | 2024-03-08 | 南京航空航天大学 | 基于参数化时域传递率的结构参数辨识及方差计算方法 |
CN117669211B (zh) * | 2023-12-06 | 2024-06-25 | 南京航空航天大学 | 基于参数化时域传递率的结构参数辨识及方差计算方法 |
Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN1804612A (zh) * | 2005-12-29 | 2006-07-19 | 西北工业大学 | 基于互相关函数幅值向量的随机振动结构损伤的检测方法 |
CN103207931A (zh) * | 2013-03-11 | 2013-07-17 | 武汉大学 | 基于mm和arma算法的次同步振荡模态辨识方法 |
CN103926538A (zh) * | 2014-05-05 | 2014-07-16 | 山东大学 | 基于aic准则的变阶数rc等效电路模型及实现方法 |
CN104166804A (zh) * | 2014-08-20 | 2014-11-26 | 中国科学技术大学 | 一种基于时频域单源点稀疏成分分析的工作模态辨识方法 |
CN106326530A (zh) * | 2016-08-10 | 2017-01-11 | 北京理工大学 | 一种基于右矩阵分式模型的时变结构模态参数辨识方法 |
CN106354695A (zh) * | 2016-08-22 | 2017-01-25 | 北京理工大学 | 一种仅输出线性时变结构模态参数辨识方法 |
-
2017
- 2017-08-24 CN CN201710732229.8A patent/CN107561934B/zh active Active
Patent Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN1804612A (zh) * | 2005-12-29 | 2006-07-19 | 西北工业大学 | 基于互相关函数幅值向量的随机振动结构损伤的检测方法 |
CN103207931A (zh) * | 2013-03-11 | 2013-07-17 | 武汉大学 | 基于mm和arma算法的次同步振荡模态辨识方法 |
CN103926538A (zh) * | 2014-05-05 | 2014-07-16 | 山东大学 | 基于aic准则的变阶数rc等效电路模型及实现方法 |
CN104166804A (zh) * | 2014-08-20 | 2014-11-26 | 中国科学技术大学 | 一种基于时频域单源点稀疏成分分析的工作模态辨识方法 |
CN106326530A (zh) * | 2016-08-10 | 2017-01-11 | 北京理工大学 | 一种基于右矩阵分式模型的时变结构模态参数辨识方法 |
CN106354695A (zh) * | 2016-08-22 | 2017-01-25 | 北京理工大学 | 一种仅输出线性时变结构模态参数辨识方法 |
Non-Patent Citations (3)
Title |
---|
CHRISTOF DEVRIENDT,等: "The use of transmissibility measurements in output-only modal analysis", 《MECHANICAL SYSTEMS AND SIGNAL PROCESSING》 * |
SI-DA ZHOU,等: "Output-only modal parameter estimator of linear time-varying structural systems based on vector TAR model and least squares support vector machine", 《MECHANICAL SYSTEMS AND SIGNAL PROCESSING》 * |
周思达,等: "基于响应传递率的非白随机激励仅输出结构模态参数辨识", 《振动与冲击》 * |
Cited By (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN110688729A (zh) * | 2019-08-26 | 2020-01-14 | 南京航空航天大学 | 基于自适应卡尔曼滤波的lstm-idm跟驰特性融合方法、存储介质及设备 |
CN110688729B (zh) * | 2019-08-26 | 2023-07-14 | 南京航空航天大学 | 基于自适应卡尔曼滤波的lstm-idm跟驰特性融合方法、存储介质及设备 |
WO2021077467A1 (en) * | 2019-10-24 | 2021-04-29 | Dalian University Of Technology | Method of complex modal identification for structure with proportional damping |
CN117669211A (zh) * | 2023-12-06 | 2024-03-08 | 南京航空航天大学 | 基于参数化时域传递率的结构参数辨识及方差计算方法 |
CN117669211B (zh) * | 2023-12-06 | 2024-06-25 | 南京航空航天大学 | 基于参数化时域传递率的结构参数辨识及方差计算方法 |
Also Published As
Publication number | Publication date |
---|---|
CN107561934B (zh) | 2020-04-21 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN104133950B (zh) | 一种悬臂梁运行模态分析实验方法及装置 | |
CN106897717A (zh) | 基于环境激励数据的多次测试下贝叶斯模型修正方法 | |
CN104993480B (zh) | 基于递推随机子空间的电力系统低频振荡在线辨识方法 | |
CN101706355A (zh) | 基于NExT/ARMA的结构响应分析方法 | |
CN102998118B (zh) | 一种基于形态学滤波和复杂度测度的轴承定量诊断方法 | |
CN112033463B (zh) | 一种核动力设备状态评估与预测一体化方法和系统 | |
CN101916241B (zh) | 一种基于时频分布图的时变结构模态频率辨识方法 | |
CN104132791A (zh) | 一种基于脉冲激励的运行模态分析实验方法及装置 | |
CN108594660B (zh) | 一种时不变结构的工作模态参数识别方法和系统 | |
CN102902981B (zh) | 基于慢特征分析的暴力视频检测方法 | |
CN104165742A (zh) | 一种基于互谱函数的运行模态分析实验方法及装置 | |
CN111307452B (zh) | 一种时变转速下旋转机械智能故障诊断方法 | |
CN102562239A (zh) | 一种航空发动机排气温度的监测方法 | |
CN107271127A (zh) | 基于自迭代主元抽取的工作模态参数识别方法及装置 | |
Mao et al. | The construction and comparison of damage detection index based on the nonlinear output frequency response function and experimental analysis | |
CN107622801A (zh) | 疾病概率的检测方法和装置 | |
EP2342603A2 (en) | Method and apparatus for creating state estimation models in machine condition monitoring | |
CN110866448A (zh) | 基于卷积神经网络和短时傅里叶变换的颤振信号分析方法 | |
CN108919116B (zh) | 基于MCCKAF-FFT-Softmax的海流发电机不平衡定子电流故障诊断方法 | |
CN102829967A (zh) | 一种基于回归模型系数变化的时域故障识别方法 | |
CN107561934A (zh) | 基于多参考时域传递率的仅输出结构模态参数辨识方法 | |
CN107066736A (zh) | 一种基于压缩采样的模态分析及结构冲击监测方法 | |
CN109472288A (zh) | 一种抽水蓄能机组振动混合特征提取与分类方法 | |
CN114925526B (zh) | 一种结合多工况响应的结构模态参数识别方法 | |
CN105844055A (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 |