基于时变公分母模型的时频域时变结构模态参数辨识方法

技术领域

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

背景技术

在现实的生产和生活中，许多工程结构表现出这样的时变特征，如列车激 励中的车桥系统、飞行过程中液体燃料逐渐减少的运载火箭、气动力附加效应 下的飞机、柔性可展开的几何可变航天器、旋转机械等。

在国内航天领域，大型空间站、新一代运载火箭、大柔性展开式卫星等新 一代的航天器已被列入我国最新的航天发展规划中，成为未来几十年中国航天 器发展的主要方向。大型空间站、新一代运载火箭、大柔度展开式卫星的结构 在运行中无一例外的存在着较强的时变因素，如未来大型空间站的空间对接问 题，现役和未来新一代运载火箭的燃料质量快速消耗，以及大柔度可展开式卫 星的空间展开等。因此，作为时变结构动力学特性分析的重要方法和途径，时 变结构模态参数辨识研究将成为未来航天器结构动力学研究的重点之一。时变 结构模态参数辨识可以辨识时变结构的模态频率、模态振型和模态阻尼，这些 参数具有重要的物理意义，可以为时变结构的结构设计、结构健康监测、结构 故障诊断、结构振动控制等方面的应用提供有力的支持。

按照采用的数学模型的差异区分，现有的时变结构模态参数辨识的方法主 要有四类：

第一类是从传统的时不变结构模态参数辨识发展而来的基于在线递推技术 的时变模态参数辨识方法。

这类方法的基础是传统的时不变结构模态参数辨识方法，不同之处为在每 一个时刻数据序贯地被考虑，老的数据逐渐被遗忘，新的数据不断地加进来， 模态参数的估计值在每一个时刻被修正。这类方法存在两方面的缺陷：第一， 存在着观测数据及遗忘因子(算法)的选取问题，需要在识别精度和跟踪能力 二者之间做折中，并且对于不同结构的相关选取的适应性也很难解决；第二， 这类方法来自传统的模态参数辨识方法，需要结构的输入和输出两方面的响应 信息，因此很难运用于如在飞航天器等只能得到输出响应信号的结构模态参数 辨识。

第二类是基于短时不变假设的模态参数辨识方法。

这类方法将数据(结构响应)划分成一个个小的时间段，并在每一个时间 段内把结构参数看成是时不变的，然后将每一段内识别值用一定的数据处理技 术(如曲线拟合)处理得到模态参数随时间变化的规律。它的特点是估计后一 段时间的模态参数时没有用到前面各段的数据信息，对参数变化较快的结构为 使估计精度提高必须选取很短的数据段。此方法包括现今较为常用的基于状态 空间模型的递推的随机子空间辨识法(N4SID)和时间相关自回归滑动平均模型 (Time-dependent ARMA，TARMA)方法。这类方法的时变结构模态参数辨识 方法发展时间最长，发展的也最为完善。但是一些固有的问题限制了其进一步 发展和应用：第一，短时不变假设限制了此类方法对于快变、突变参数辨识方 面的应用；第二，此类方法需要形式固定、明确的数学模型，如状态空间模型、 时间序列的自回归滑动平均模型等，因此，在辨识中模型的定阶问题十分突出， 模型阶数的不确定将引入无物理意义的虚假模态，造成辨识结果不可用，模型 阶次合理选取、虚假模态的判断等问题更需要进一步深入研究；第三，作为两 种主流的基于短时不变假设的模态参数辨识方法——递推的随机子空间辨识法 和时间相关自回归滑动平均模型各自存在着一些其它的问题：基于状态空间模 型的堆积子空间方法不可避免地要使用QR分解、特征值分解(EVD)或者奇 异值分解(SVD)技术，这必然会带来方法数值实现上的复杂性，对于大型工程结 构，尤其对有在线以及快速辨识要求的问题，这还需要进一步进行研究；基于时间 序列模型的辨识方法研究都不能回避参数跟踪算法的设计，虽然各种改进的最小 二乘法、各种滤波方法不断提出，但是当相同模型使用不同跟踪算法，以及不同模 型应用相同算法结果差异非常大。

第三类是人工神经网络的时变模态参数辨识方法。

人工神经网络已经被广泛地应用于非线性系统辨识问题，并取得良好的效果 但大部分研究工作还仅局限于时不变系统，只是近几年来被推广到时变系统。将 人工神经网络用于时变模态参数辨识领域研究的公开发表的文献很少，其主要 集中在针对简单结构(系统)的机理性研究。对于真实的复杂结构还存在算法 复杂、计算效率低和辨识精度差等问题。

第四类方法是基于时频分析的非参数化时频域的时变结构模态参数辨识方 法。

从信号分析的角度来看，时变结构在工作环境下的结构动力学响应信号是 非平稳随机信号。

经典傅立叶变换经过一个世纪的发展，已成为信号处理领域最强有力的分 析方法和工具，这主要是由它的正交性和鲜明的物理意义以及快速简洁的计算 方法所决定的。但是，由于傅立叶变换是对时间求积，去掉了非平稳信号中的 时变信号，因而要求信号是平稳的，对时变非平稳信号难以充分刻画。为了满 足对突变信号、非平稳信号分析的要求，1946年，Gabor提出了加窗傅立叶变 换分析方法，亦称短时傅立叶变换(short-time Fourier transform，STFT)，通过 适当窗函数的选取，就可以实现一定程度上的时频分析，但是由于时间分辨率 与频率分辨率要受到窗函数宽度的限制，总是不能同时到达最佳。1948年，Ville 提出了著名的维格纳-威尔分布(Wigner-Ville distribution，WVD)。它作为一 种能量型时频联合分布，与其他时频分布相比有许多优良性质，如真边缘性、 弱支撑性、平移不变性等，是一个非常有用的非平稳信号分析工具。由于多信 号的维格纳-威尔分布出现交叉项，在不少场合会限制其应用效果，所以后来研 究人员在此基础上，提出了多种改进形式，如指数分布、广义指数分布、广义 双线性时频分布等，其中广义双线性时频分布又称为科恩类能量型时频分布。 后来在此基础上，人们又提出了科恩类时频分布等方法，这些时频分析方法在 非平稳随机信号分析中得到了广泛的应用并取得了许多令人满意的结果。

近十年，由于时频分析在非平稳随机信号分析方面的优势，越来越多的研 究者运用时频分析来进行时变和非线性系统辨识的研究。时频分析方法对时变 和非线性结构模态参数进行辨识也渐渐成为模态参数辨识研究领域的热点之 一。2000年Ghanem将结构动力学控制微分方程在一系列小波基上展开，用小 波系数来代替原来的物理响应，并采用了求解展开方程的方法辨识了系统的模 态参数；2003年Zhang和Xu通过对一个简单的时变结构响应的Gabor变换辨 识了结构的模态频率；2007年Roshan-Ghias采用解析推导的方式对一个单自由 度系统和一个三自由度系统自由振动下的响应进行了WVD和SPWVD变换，并 根据变换结果估计了系统的模态频率和阻尼比。

现有的基于时频分析的时频域时变结构模态参数辨识方法都是非参数化 的，虽然有的方法能够很好的辨识出时变结构的模态频率，但是非参数化的方 法不同程度的依赖于使用者的主观意识和经验，并且对于随机激励下的时变结 构模态阻尼比辨识还没有很好的办法。

发明内容

本发明针对航空器和航天器时变结构模态参数辨识问题，提出了一种时变 公分母模型的时频域时变结构模态参数辨识方法，其基本思路为：首先对具有 时变特性的航空器或航天器结构在工作环境下测量得到的结构动力学响应信号 进行时频分析，得到非参数化估计的对应时变结构的时间相关功率谱函数，然 后以时变公分母模型为时变结构动力学的参数化模型，通过时频域的最小二乘 方法估计出时变公分母模型的待估参数，最后利用估计出来的时变公分母模型 的待估参数计算出对应时变结构的模态频率和模态阻尼比。

本发明的具体实现步骤如下：

步骤1，根据被辨识航空器或航天器时变结构工作状态、感兴趣的时间范围 以及被辨识航空器或航天器时变结构的主要频率范围，如运载火箭约0-20Hz、 一般卫星约0-100Hz、大型可展开卫星约0-10Hz、飞机机翼约0-50Hz等，设定 辨识所需的采样时间和采样频率，并对辨识结构的结构动力学响应信号进行采 集。

步骤2，从步骤1采集的响应信号中任意选取参考信号，并结合参考信号分 别对各个响应信号进行时频分析，得到响应信号的时间相关功率谱函数 G_k(t_τ，ω_f)。其中，t_τ表示时间采样点，ω_f表示频率采样点，下标τ=1，2，...，N_τ， f=1，2，...，N_f，N_τ为时间采样数，N_f为频率采样数，k=1，2，...，N_sN_r，N_s为输 出响应数，N_r为响应信号参考点数。

本发明采用平滑伪Wigner-Ville分布(Smooth Pseudo平滑伪Wigner-Ville  Distribution，SPWVD)来计算被辨识结构响应信号的功率谱函数，实施简单， 计算效率高，且能够较好的抑制分布中交叉项。

步骤3，根据采样时间和采样频率以及被辨识航空器或航天器的时变特点， 建立时变公分母模型如下：

${\hat{G}}_k(t_{\tau },{\omega }_f)=\frac{B_k(t_{\tau },{\omega }_f)}{A(t_{\tau },{\omega }_f)}---\left(1\right)$

其中，分子多项式和分母多项式分别为：

其中，为时频基函数(时间多项式和频率多项式分别为i和j阶， i=0，1，2，...，n_t，j=0，1，2，...，n_ω，n_t为时间多项式阶数，n_ω为频率多项式阶数)， 分子多项式系数b_k，i，j和公分母a_i，j写为向量形式：

$B_{k,j}={[b_{k,0,j},b_{k,1,j},...,b_{k,n_t,j}]}^T,$$A_j={[a_{0,j},a_{1,j},...,a_{n_t,j}]}^T---\left(4\right)$

令：

${\beta }_k={[B_{k,0}^T,B_{k,1}^T,...,B_{k,n_{\omega }}^T]}^T;$$\alpha ={[A_0^T,A_1^T,...,A_{n_{\omega }}^T]}^T---\left(5\right)$

则有：

$\theta ={[{\beta }_1^T,...,{\beta }_k^T,...,{\beta }_{N_o{N}_i}^T,{\alpha }^T]}^T---\left(6\right)$

其中，θ为时变公分母模型中的待估参数向量。

步骤4，采用最小二乘参数方法得到时变公分母模型中的待估参数向量θ。

具体过程如下：

定义最小二乘的花费函数：

$l_{\mathrm{LS}}\left(\theta \right)=\underset{k=1}{\overset{N_r{N}_s}{\Sigma }}\left(\begin{array}{cc}{\beta }_k^T& {\alpha }^T\end{array}\right)\left(\begin{array}{cc}{P}_k& Q_k\\ Q_k^T& R_k\end{array}\right)\left(\begin{array}{c}{\beta }_k\\ \alpha \end{array}\right)---\left(7\right)$

其中，$P_k=\mathrm{Re}\left({\Theta }_k^H{\Theta }_k\right),$$Q_k=\mathrm{Re}\left({\Theta }_k^H{\Xi }_k\right)$和$R_k=\mathrm{Re}\left({\Xi }_k^H{\Xi }_k\right).$

其中W_k(t_τ，ω_f)为权函数。

利用下式求得约束后的待估公分母参数向量α′：

D′α′=b′      (10)

其中，$D^{\prime }=D_{[1:(n_t+1)(n_{\omega }+1)-\mathrm{1,1}:(n_t+1)(n_{\omega }+1)-1]},$$b^{\prime }=-D_{[1:(n_t+1)(n_{\omega }+1)-1,(n_t+1)(n_{\omega }+1)]},$$D=2{\Sigma }_{k=1}^{N_r{N}_s}(R_k-Q_k^H{P}_k^{-1}{Q}_k).$

令$\alpha =\left(\begin{array}{c}{\alpha }^{\prime }\\ 1\end{array}\right),$并通过下式计算出待估分子参数向量β_k：

${\beta }_k=P_k^{-1}{Q}_k\alpha ---\left(11\right)$

从而得到时变公分母模型中的待估参数向量θ，进而确定时变公分母模型。

步骤5，根据被辨识航空器或航天器模态参数的应用需要，如运载火箭振动 控制系统的控制频率、卫星结构健康监测时的时间间隔要求等，给定需要计算 模态参数的时刻点t_τ′，并利用步骤4得到的公分母参数向量α以及时变公分母 模型，计算出给定时刻t_τ′下辨识得到的模态频率f_r和模态阻尼比ξ_r：

$f_r=\frac{\mathrm{Im}\left({\lambda }_r\right)}{2\pi },$${\xi }_r=\frac{\mathrm{Re}\left({\lambda }_r\right)}{\left|{\lambda }_r\right|}---\left(12\right)$

其中，Im和Re分别为取括号中值的虚部和实部，λ_r为以α为系数的t_τ′时刻下的 时变公分母模型中分母多项式A(t_τ′，ω)的第r个根，ω为被辨识航空器或航天器 的频率变量。

有益效果

本发明从参数化的时频域的角度出发，给出一种基于时频分析的参数化时 频域的时变结构模态参数辨识方法，物理意义清晰；适用于航空和航天工程应 用领域的时变结构的模态参数辨识，并且所需使用者的参与度较低，具有使用 简单和方便的特点。

附图说明

图1为具体实施方式中的三自由度弹簧-阻尼器-质量系统；

图2为具体实施方式中的三自由度时变结构随机响应的SPWVD；

图3为具体实施方式中辨识得到的模态参数。

具体实施方式

为了更好的说明本发明的目的和优点，下面通过一个随机激励下的三自由 度时变结构实例来模拟工作状态下的航空器和航天器时变结构，对本发明进行 进一步说明。

本实施例的三自由度弹簧-阻尼器-质量系统，如图1所示。

三自由度系统的参数为k₁=k₂=k₃=10⁵，c₁=1.0，c₂=0.5，c₃=0.5，初始质量 为m₁(0)=0.2，m₂(0)=0.1，m₃(0)=0.1。系统的动力学控制方程为：

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

其中M(t)为时变的质量矩阵，M(t)=M₀(1-0.5t)，M₀为初始时刻的质量矩阵， 阻尼矩阵和刚度矩阵为

$C=\left(\begin{array}{ccc}0.2& 0& 0\\ 0& 0.2& 0\\ 0& 0& 0.1\end{array}\right),$$C=\left(\begin{array}{ccc}1.5& -0.5& 0\\ -0.5& 1.0& -0.5\\ 0& -0.5& 0.5\end{array}\right),$$K=\left(\begin{array}{ccc}2& -1& 0\\ -1& 2& -1\\ 0& -1& 1\end{array}\right)\times {10}^5$

f(t)为施加在三自由度上的幅值为100个单位的白噪声激励。

本实施例中的激励为作用在系统第一个自由度上的仿真Gauss白噪声。系 统响应采用Newmark-β数值积分方法计算得到(γ=0.5，β＝0.1)，其中积分步长 为1/4096s。信号采样频率为2048Hz，采样时间为1s。

本实施例中，基于时变公分母模型的时频域时变结构模态频率辨识的具体 实现步骤如下：

步骤1，三个自由度加速度为辨识所用的响应信号，采样频率为1024Hz， 采样时间为1s。

步骤2，将第一个响应点的加速度响应信号作为进行时频分析时的参考信 号，并通过平滑伪Wigner-Ville分布得到各响应信号与参考信号的时间相关功率 谱函数。40次平均的平滑伪Wigner-Ville分布如图2所示。

步骤3，建立时变公分母模型，然后给定时变公分母模型中的时间采样数 N_τ=32、频率采样数N_f=256、时间多项式阶数n_t=5和频率多项式阶数 n_ω=32，并初始化待估参数向量θ。

步骤4，采用最小二乘方法估计出参数向量θ。其中，W_k(t_τ，ω_f)为权函数， 本实例设为1。

步骤5，根据参数向量θ中的公分母参数向量α，计算出每一时刻的结构模 态频率和模态阻尼比。

采用本实例的参数辨识出来的时变结构模态频率和模态阻尼比如图3所示， 其中实线为辨识值，黑圈为理论值。

由此可见本发明能够很好的辨识出时变结构的模态频率，较好的辨识出时 变结构的模态阻尼比。由于其仅需要结构的响应信号作为输入，因此，适合对 工作状态的时变结构进行模态频率辨识。另一方面，由于在整个过程中使用者 仅需设置一个初步的参数即可得到最终的模态频率和模态阻尼比的辨识结果， 用户参与度较低，使用十分简单方便。

以上所述的具体描述，对发明的目的、技术方案和有益效果进行了进一步 详细说明，所应理解的是，以上所述仅为本发明的具体实施例，用于解释本发 明，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的 任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。