CN106326530A - 一种基于右矩阵分式模型的时变结构模态参数辨识方法 - Google Patents

Abstract

一种基于右矩阵分式模型的时变结构模态参数辨识方法简称本方法，核心思想为：首先利用测量得到的时变结构在工作状态下的加速度信号进行时频分析，得到非参数化的时变结构的时间相关功率谱函数；之后建立时变结构的右矩阵分式模型，利用最小二乘方法对模型进行估计，得到右矩阵分式模型的待估参数；最后将右矩阵分式模型的待估参数转化为时变结构的模态参数(频率和阻尼比)。本方法所依托的系统为信号采集系统，功能为采集时变结构的加速度信号。本方法提高了时频域时变结构模态参数辨识方法对密集模态和弱模态的辨识效果。

Description

一种基于右矩阵分式模型的时变结构模态参数辨识方法

技术领域

本发明涉及一种基于右矩阵分式模型的时变结构模态参数辨识方法，属于结构动力学技术领域。

背景技术

工程结构在航空航天、汽车、建筑、船舶、机械等众多领域内得到了广泛的应用。工程结构具有承受和传递外载荷，维持一定外形的作用，其动力学特性可以为工程应用中系统的设计与分析提供一定的参考。结构的模态参数(频率和阻尼比等)是描述结构动力学特性的一种手段，模态参数辨识可获取结构的模态参数，在工程领域具有多方面的意义和应用价值。

随着工程结构的大型化、复杂化、高速化，越来越多的结构的时变特性渐渐变得不可忽略。例如，超高声速飞行器在飞行过程中，气动加热会结构的刚度和阻尼发生变化；运载火箭飞行过程中，燃料的消耗将使其总质量和质心位置不断变化；汽车通过桥梁时，桥梁和汽车系统的质量分布将发生变化。对于以上结构特性随时间变化较大的系统，传统的时不变的模态参数辨识方法将很难适用。因此，发展时变结构的模态辨识方法是非常必要的。

目前时变结构的模态参数辨识方法按作用域可分为时域方法和时频域方法。其中，时域方法包括基于时间序列模型的模态参数辨识方法和基于状态空间模型的模态参数辨识方法，时域方法直接利用时域信号并且直接建立时间相关的参数化模型，减少了中间环节所带来的误差，但此类方法对模型阶数较为敏感，并且阶数的选择也较为复杂和耗时；时频域方法包括基于时频分析的非参数方法，基于Hilbert-Huang变换的非参数化方法以及参数化的时频域时变结构模态参数辨识方法。上述两类非参数化方法由于没有建立结构的参数化模型，无法完整的辨识出系统的所有特性和参数，辨识结果的质量依赖于使用者的经验和对物理模型的了解程度。同时，该类方法对于阻尼辨识和模态振型的辨识具有一定困难。

申请号：201210424594.X，标题为“基于时变公分母模型的时频域时变结构模态参数辨识方法”提出了一种参数化的时频域模态参数辨识方法，由于该方法对模型阶数的选择较不敏感，减少了由于模型阶数选择所带来的计算量的增大。同时该方法在高阻尼结构和振型辨识方面具有潜在优势，并且可以通过指定带宽提高辨识精度。然而，该方法对于密集模态的辨识具有一定困难，对弱模态的辨识能力也具有提高的空间。

上述方法虽然减小了计算量，在高阻尼结构和振型辨识方面具有潜在优势，并且可以通过指定带宽提高辨识精度，但仍具有改进的空间。本发明的目的是提高上述基于时变公分母模型的时频域时变结构模态参数辨识方法在具有密集模态和弱模态的时变结构的辨识能力，所提出的一种基于右矩阵分式模型的时变结构模态参数辨识方法是一种参数化的时频域方法，该方法继承了基于时变公分母模型的时频域时变结构模态参数辨识方法在减小计算量、高阻尼结构辨识、振型辨识和指定带宽等方面的优势。同时通过引入多个参考点，可提高对具有密集模态和弱模态时变结构的辨识能力。

发明内容

本发明的目的是针对现有基于时变公分母模型的时频域时变结构模态参数辨识方法的对密集模态和弱模态辨识结果较差的技术缺陷，提出了一种基于右矩阵分式模型的时变结构模态参数辨识方法。

一种基于右矩阵分式模型的时变结构模态参数辨识方法的核心思想为：首先利用测量得到的时变结构在工作状态下的加速度信号进行时频分析，得到非参数化的时变结构的时间相关功率谱函数；之后建立时变结构的右矩阵分式模型，利用最小二乘方法对模型进行估计，得到右矩阵分式模型的待估参数；最后将右矩阵分式模型的待估参数转化为时变结构的模态参数(频率和阻尼比)。

一种基于右矩阵分式模型的时变结构模态参数辨识方法简称本方法，本方法所依托的系统为信号采集系统，功能为采集时变结构的加速度信号。

一种基于右矩阵分式模型的时变结构模态参数辨识方法，具体步骤如下：

步骤1：利用信号采集系统采集得到的时变结构的加速度信号，采用已有的非参数化估计方法进行时频分析，得到时间相关的功率谱函数；

其中，已有的非参数化估计方法主要包括Wigner-Ville分布(WVD)以及平滑伪Wigner-Ville分布(SPWVD)；

步骤2：建立时变结构的右矩阵分式模型；

时变结构的右矩阵分式模型可以表示为公式(1)的形式：

${\hat{G}}_k(t_{\mathrm{\τ}},{\mathrm{\ω}}_f,\mathrm{\θ})=B_k(t_{\mathrm{\τ}},{\mathrm{\ω}}_f)A^{-1}(t_{\mathrm{\τ}},{\mathrm{\ω}}_f)---\left(1\right)$

其中，为时变结构的右矩阵分式模型，t_τ表示时间变量，ω_f表示频率变量，θ代表待估参数，τ表示时间采样点，f表示频率采样点；B_k为1×N_i的分子系数矩阵，A为N_i×N_i的分母系数矩阵，N_i为参考点的个数；k表示第k个输出点；A^-1(t_τ,ω_f)表示对A(t_τ,ω_f)求逆矩阵；B_k和A可由式(2)和式(3)确定：

其中，∑表示求和符号；表示时频基函数，可以写为时间基函数p_i(t_τ)和频率基函数z_j(ω_f)的乘积，即下标i＝0,1,...,n_t表示时间多项式的阶数，n_t表示时间多项式的最大阶数；j＝0,1,2,...,n_ω表示频率多项式的阶数，n_ω表示频率多项式的最大阶数；b_k,i,j和a_i,j分别为右矩阵分式模型的分子待估参数和分母待估参数，且分子待估参数和分母待估参数统称待估参数；

待估参数可以表示成统一的向量形式如下：

其中，[]表示由向量或矩阵组成的新矩阵，表示复数域，T表示转置，

其中，

步骤3：采用时频域时变结构最小二乘方法对步骤2中建立的右矩阵分式模型进行估计，具体为：

步骤3.1：确定时频域时变结构最小二乘方法的费用函数，所述费用函数为估计右矩阵分式模型待估参数的目标函数；

时频域时变结构最小二乘方法的费用函数可以写为如公式(8)所示：

$l_{LS}=\underset{k=1}{\overset{N_o}{\Σ}}\left[\begin{array}{cc}{{\mathrm{\β}}_k}^H& {\mathrm{\α}}^H\end{array}\right]\left[\begin{array}{cc}{P}_k& Q_k\\ Q_k^T& R_k\end{array}\right]\left[\begin{array}{c}{\mathrm{\β}}_k\\ \mathrm{\α}\end{array}\right]---\left(8\right)$

其中，H表示共轭转置；

步骤3.2：计算右矩阵分式模型的待估参数；

参数α和β_k可以由下式得到：

$\begin{array}{c}\mathrm{\α}=\[{\left(D_p\right)}^{-1}{b}_p;I_{N_i}\]\\ {\mathrm{\β}}_k=-P_k^{-1}{Q}_k\mathrm{\α},k=1,2,...,N_o\end{array}---\left(9\right)$

式中，矩阵D_p＝D{1:((n_t+1)(n_ω+1)-1)N_i,1:((n_t+1)(n_ω+1)-1)N_i}，b_p＝-D{1:((n_t+1)(n_ω+1)-1)N_i,((n_t+1)(n_ω+1)-1)N_i+1:(n_t+1)(n_ω+1)N_i}，{:,:}表示选择矩阵的一部分，如D{1:((n_t+1)(n_ω+1)-1)N_i,1:((n_t+1)(n_ω+1)-1)N_i}表示选择矩阵D从第1行到第((n_t+1)(n_ω+1)-1)N_i行，从第1列到((n_t+1)(n_ω+1)-1)N_i列；为单位矩阵，矩阵Re表示复数的实部，X_k为N_τN_f×(n_t+1)(n_ω+1)的矩阵，Y_k为N_τN_f×(n_t+1)(n_ω+1)N_i的矩阵，可以表示为如(10)和(11)所示：

其中，为步骤1中得到的时间相关的功率谱函数，符号表示kronecker积；

步骤4：将步骤3中估计得到的右矩阵分式模型的待估参数转化为模态参数；

其中，模态参数包括模态频率和阻尼比，分别记为f_r和ξ_r，并通过(12)计算：

$f_r=\frac{Im\left({\mathrm{\λ}}_r\right)}{2\mathrm{\π}},{\mathrm{\ξ}}_r=\frac{\mathrm{Re}\left({\mathrm{\λ}}_r\right)}{\left|{\mathrm{\λ}}_r\right|}---\left(12\right)$

其中，符号Im表示复数的虚部，λ_r为广义伴随矩阵A_c(t_τ)的特征值，即λ_r＝eig(A_c(t_τ))，广义伴随矩阵A_c(t_τ)如(13)所示：

其中，A_i＝P_tα{(i-1)(n_t+1)N_i+1:t_τN_i,i(n_t+1)N_i,:}，i＝1,2,...,n_ω；

其中，

至此，从步骤1到步骤4，完成了一种基于右矩阵分式模型的时变结构模态参数辨识方法。

有益效果

本发明一种基于右矩阵分式模型的时变结构模态参数辨识方法，与基于时变公分母模型的时频域时变结构模态参数辨识方法相比，具有如下有益效果：

本发明所提方法提高了时频域时变结构模态参数辨识方法对密集模态和弱模态的辨识效果。

附图说明

图1为本发明一种基于右矩阵分式模型的时变结构模态参数辨识方法的流程图；

图2为本发明一种基于右矩阵分式模型的时变结构模态参数辨识方法及实施例中的四自由度时变质量-弹簧-阻尼系统；

图3为本发明一种基于右矩阵分式模型的时变结构模态参数辨识方法实施例中的四自由度的时变质量-弹簧-阻尼系统的模态参数的理论值(a)模态频率(b)模态阻尼比；

图4为本发明一种基于右矩阵分式模型的时变结构模态参数辨识方法实施例中的平滑伪Wigner-Ville分布(SPWVD)非参数化方法估计的时间相关功率谱(a)G12(b)G22(c)G32(d)G42；

图5为本发明一种基于右矩阵分式模型的时变结构模态参数辨识方法实施例中的对四自由度时变质量-弹簧-阻尼系统的模态参数的估计值(a)模态频率(b)模态阻尼比；

图6为本发明一种基于右矩阵分式模型的时变结构模态参数辨识方法实施例中的基于公分母模型的时频域时变结构模态参数辨识方法对四自由度时变质量-弹簧-阻尼系统的模态参数的估计值(a)模态频率(b)模态阻尼比。

具体实施方式

本发明以一个四自由度时变质量-弹簧-阻尼系统模拟工作状态下的时变结构，并对本发明进行具体说明。

图1是本发明一种基于右矩阵分式模型的时变结构模态参数辨识方法的流程图。

由图1可以看出，本发明首先对时变结构加速度信号进行1、时频分析，得出时间相关功率谱，其次2、建立时变结构的右矩阵分式模型，之后3、采用最小二乘方法对右矩阵分式模型进行估计，得到右矩阵分式模型的待估参数，最后4、将右矩阵分式模型的待估参数转化为模态参数。

具体到本实施例，阐述的是本发明所提的方法用于四自由度时变质量-弹簧-阻尼系统，简称四自由度系统，如图2所示。

图2中m₁、m₂、m₃和m₄分别为四自由度系统的四个质量块的质量，c₁、c₂、c₃和c₄分别为四自由度系统的四个阻尼器的阻尼系数，k₁、k₂、k₃和k₄分别为四自由度系统的四个弹簧的刚度系数，四自由度系统的结构动力学方程可以表示为：

$M\left(t\right)\stackrel{\·\·}{x}\left(t\right)+C\stackrel{\·}{x}\left(t\right)+Kx\left(t\right)=f\left(t\right)---\left(15\right)$

式中和x(t)分别表示为加速度、速度和位移向量；f(t)表示施加在第二个自由度上的高斯白噪声激励；M(t)、C和K分别表示结构的质量、阻尼和刚度矩阵，其中质量矩阵随时间变化，可分别表示如下：

$M=(1-0.5t)\left[\begin{array}{cccc}{m}_1& & & \\ & m_2& & \\ & & m_3& \\ & & & m_4\end{array}\right],C=\left[\begin{array}{cccc}{c}_1+c_2& -c_2& 0& 0\\ -c_2& c_2+c_3& -c_3& 0\\ 0& -c_3& c_3+c_4& -c_4\\ 0& 0& -c_4& c_4\end{array}\right],$

$K=\left[\begin{array}{cccc}{k}_1+k_2& -k_2& 0& 0\\ -k_2& k_2+k_3& -k_3& 0\\ 0& -k_3& k_3+k_4& -k_4\\ 0& 0& -k_4& k_4\end{array}\right]---\left(16\right)$

式中，m₁＝m₂＝m₄＝0.1kg，m₃＝0.5kg，k₁＝k₂＝k₃＝k₄＝0.98×10⁵N·m-¹，c₄＝15N·m-¹·s，c₁＝c₁＝c₃＝2.5N·m^-1·s。

本实施例中的四自由度时变质量-弹簧-阻尼系统的模态参数的理论值如图3所示；其中，图3(a)为频率的理论值，横坐标Time(s)代表时间，纵坐标Frequency(Hz)代表频率，╳、○、□和△分别代表第1、2、3和4阶频率；图3(b)为阻尼比的理论值，横坐标Time(s)代表时间，纵坐标Damping Ratio(％)代表阻尼比，╳、○、□和△分别代表第1、2、3和4阶阻尼比。

从图3(a)中可以看出，该系统的第二阶模态频率和第三阶模态频率较为接近。

本实施例中，一种基于右矩阵分式模型的时变结构模态频率辨识方法的具体实现步骤如下：

步骤A：利用数值算法计算得到四自由度系统的加速度信号，信号采样频率为1024Hz，采样时间为1s；利用非参数化估计方法中的平滑伪Wigner-Ville分布(SPWVD)方法，得到时间相关的功率谱函数如图4所示，Time(s)坐标表示时间，Frequency(Hz)坐标表示频率，Amplitude坐标表示时间相关的功率谱的幅值。

图4中颜色较深的且突起的位置为所辨识系统的模态频率位置；

步骤B：建立时变结构的右矩阵分式模型，给定右矩阵分式模型中的时间采样数N_τ＝32，频率采样数N_f＝256，时间多项式阶数n_t＝5，频率多项式阶数n_ω＝12，参考点数N_i＝2，参考点为第2个和第3个自由度。

步骤C：采用时频域时变结构最小二乘方法对步骤2中建立的右矩阵分式模型进行估计，根据本发明所提方法编写程序，计算得到右矩阵分式模型的待估参数θ；

步骤D：将右矩阵分式模型参数转化为结构的模态参数(模态频率和阻尼比)，如图5所示，其中虚线为理论值，实线为估计值；图5(a)为所辨识结构的频率，横坐标Time(s)代表时间，纵坐标Frequency(Hz)代表频率；图5(b)为所辨识结构的阻尼比，其中，Mode 1、Mode 2、Mode 3和Mode 4分别对应第1、2、3和4阶阻尼比；横坐标Time(s)代表时间，纵坐标表示阻尼比；

基于公分母模型的时频域时变结构模态参数辨识方法对四自由度时变质量-弹簧-阻尼系统的模态参数的估计结果如图6所示，其中虚线为理论值，实线为估计值；图6(a)为所辨识结构的频率，横坐标Time(s)代表时间，纵坐标Frequency(Hz)代表频率；图6(b)为所辨识结构的阻尼比，其中，Mode 1、Mode 2、Mode 3和Mode 4分别对应第1、2、3和4阶模态的阻尼比；横坐标Time(s)代表时间，纵坐标表示阻尼比；

对比基于公分母模型的时频域时变结构模态参数辨识方法对四自由度时变质量-弹簧-阻尼系统的模态参数的估计值，可以看出本发明提出的一种基于右矩阵分式模型的时变结构模态参数辨识方法可以估计全部四阶模态；而基于公分母模型的时频域时变结构模态参数辨识方法只能估计出两个频率相近模态中的一个，即：第2阶模态。

同时，对于较弱的第1阶模态，本发明提出的一种基于右矩阵分式模型的时变结构模态参数辨识方法可以较好的估计出来，而基于公分母模型的时频域时变结构模态参数辨识方法未能估计出来。说明本发明所提出的方法可以提高时频域时变结构模态参数辨识方法对密集模态和弱模态的辨识效果，对于真实工作状态下具有密集模态的时变结构的模态参数辨识具有潜在优势。

以上所述为本发明的较佳实施例而已，本发明不应该局限于该实施例和附图所公开的内容。凡是不脱离本发明所公开的精神下完成的等效或修改，都落入本发明保护的范围。

Claims (6)

1.一种基于右矩阵分式模型的时变结构模态参数辨识方法，其特征在于：

一种基于右矩阵分式模型的时变结构模态参数辨识方法简称本方法，本方法所依托的系统为信号采集系统，功能为采集时变结构的加速度信号；

本方法的核心思想为：首先利用测量得到的时变结构在工作状态下的加速度信号进行时频分析，得到非参数化的时变结构的时间相关功率谱函数；之后建立时变结构的右矩阵分式模型，利用最小二乘方法对模型进行估计，得到右矩阵分式模型的待估参数；最后将右矩阵分式模型的待估参数转化为时变结构的模态参数(频率和阻尼比)。

2.如权利要求1所述的一种基于右矩阵分式模型的时变结构模态参数辨识方法，其特征还在于：

本方法，具体步骤如下：

步骤1：利用信号采集系统采集得到的时变结构的加速度信号，采用已有的非参数化估计方法进行时频分析，得到时间相关的功率谱函数；

步骤2：建立时变结构的右矩阵分式模型；

步骤3：采用时频域时变结构最小二乘方法对步骤2中建立的右矩阵分式模型进行估计；

步骤4：将步骤3中估计得到的右矩阵分式模型的待估参数转化为模态参数；

至此，从步骤1到步骤4，完成了一种基于右矩阵分式模型的时变结构模态参数辨识方法。

3.如权利要求2所述的一种基于右矩阵分式模型的时变结构模态参数辨识方法，其特征还在于：

步骤1中，已有的非参数化估计方法主要包括Wigner-Ville分布(WVD)以及平滑伪Wigner-Ville分布(SPWVD)。

4.如权利要求2所述的一种基于右矩阵分式模型的时变结构模态参数辨识方法，其特征还在于：

步骤2中，时变结构的右矩阵分式模型可以表示为公式(1)的形式：

${\hat{G}}_k(t_{\mathrm{\τ}},{\mathrm{\ω}}_f,\mathrm{\θ})=B_k(t_{\mathrm{\τ}},{\mathrm{\ω}}_f)A^{-1}(t_{\mathrm{\τ}},{\mathrm{\ω}}_f)---\left(1\right)$

其中，为时变结构的右矩阵分式模型，t_τ表示时间变量，ω_f表示频率变量，θ代表待估参数，τ表示时间采样点，f表示频率采样点；B_k为1×N_i的分子系数矩阵，A为N_i×N_i的分母系数矩阵，N_i为参考点的个数；k表示第k个输出点；A^-1(t_τ,ω_f)表示对A(t_τ,ω_f)求逆矩阵；B_k和A可由式(2)和式(3)确定：

其中，∑表示求和符号；表示时频基函数，可以写为时间基函数p_i(t_τ)和频率基函数z_j(ω_f)的乘积，即下标i＝0,1,...,n_t表示时间多项式的阶数，n_t表示时间多项式的最大阶数；j＝0,1,2,...,n_ω表示频率多项式的阶数，n_ω表示频率多项式的最大阶数；b_k,i,j和a_i,j分别为右矩阵分式模型的分子待估参数和分母待估参数，且分子待估参数和分母待估参数统称待估参数；

待估参数可以表示成统一的向量形式如下：

其中，[]表示由向量或矩阵组成的新矩阵，表示复数域，T表示转置，

其中，

5.如权利要求2所述的一种基于右矩阵分式模型的时变结构模态参数辨识方法，其特征还在于：

步骤3，具体为：

步骤3.1：确定时频域时变结构最小二乘方法的费用函数，所述费用函数为估计右矩阵分式模型待估参数的目标函数；

时频域时变结构最小二乘方法的费用函数可以写为如公式(8)所示：

其中，H表示共轭转置；

步骤3.2：计算右矩阵分式模型的待估参数；

参数α和β_k可以由下式得到：

$\begin{array}{c}\mathrm{\α}=\[{\left(D_p\right)}^{-1}{b}_p;I_{N_i}\]\\ {\mathrm{\β}}_k=-P_k^{-1}{Q}_k\mathrm{\α},k=1,2,...,N_o\end{array}---\left(9\right)$

其中，矩阵D_p＝D{1:((n_t+1)(n_ω+1)-1)N_i,1:((n_t+1)(n_ω+1)-1)N_i}，b_p＝-D{1:((n_t+1)(n_ω+1)-1)N_i,((n_t+1)(n_ω+1)-1)N_i+1:(n_t+1)(n_ω+1)N_i}，{:,:}表示选择矩阵的一部分，如D{1:((n_t+1)(n_ω+1)-1)N_i,1:((n_t+1)(n_ω+1)-1)N_i}表示选择矩阵D从第1行到第((n_t+1)(n_ω+1)-1)N_i行，从第1列到((n_t+1)(n_ω+1)-1)N_i列；为单位矩阵，矩阵 Re表示复数的实部，X_k为N_τN_f×(n_t+1)(n_ω+1)的矩阵，Y_k为N_τN_f×(n_t+1)(n_ω+1)N_i的矩阵，可以表示为如(10)和(11)所示：

其中，为步骤1中得到的时间相关的功率谱函数，符号表示kronecker积。

6.如权利要求2所述的一种基于右矩阵分式模型的时变结构模态参数辨识方法，其特征还在于：

步骤4中，模态参数包括模态频率和阻尼比，分别记为f_r和ξ_r，并通过(12)计算：

$f_r=\frac{\mathrm{Im}\left({\mathrm{\λ}}_r\right)}{2\mathrm{\π}},{\mathrm{\ξ}}_r=\frac{\mathrm{Re}\left({\mathrm{\λ}}_r\right)}{\left|{\mathrm{\λ}}_r\right|}---\left(12\right)$

其中，符号Im表示复数的虚部，λ_r为广义伴随矩阵A_c(t_τ)的特征值，即λ_r＝eig(A_c(t_τ))，广义伴随矩阵A_c(t_τ)如(13)所示：

其中，A_i＝P_tα{(i-1)(n_t+1)N_i+1:t_τN_i,i(n_t+1)N_i,:}，i＝1,2,...,n_ω；

其中，