CN101916241B - 一种基于时频分布图的时变结构模态频率辨识方法 - Google Patents

Abstract

本发明涉及一种基于时频分布图的时变结构模态频率辨识方法，包括以下步骤：1获取被辨识结构的结构动力学响应信号并设定采样时间和频率；2对各个响应信号进行时频变换，得到时频分布系数，并绘制时频分布图；3将时频分布系数写成对应的能量分布形式，并重新排列为列向量；4根据各个响应的时频分布图确定将用于辨识的含有各阶时变模态频率的响应对应的时频分布区域；5采用合适的时频窗函数将时频分布图中对应各阶时变模态频率的能量时频分布最高的部分分别提取出来；6采用加权非线性最小二乘方法对各阶时变模态频率进行估计；7对辨识结果进行误差分析。本发明具有明确的物理意义，使用简单方便，具有较强的适用性和抗干扰能力。

Description

一种基于时频分布图的时变结构模态频率辨识方法

技术领域

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

背景技术

严格的说，在真实物理环境中一切结构(系统)都是时变的，只是在时间尺度上被划分成了不同的层次。当前主要研究的时变结构指在工作过程中迅速改变自身的形状或某些重要参数快速变化的结构，这里的快速表示不能忽略惯性力的作用。许多工程结构表现出这样的时变特征，如列车激励中的车桥系统、飞行过程中液体燃料逐渐减少的运载火箭、气动力附加效应下的飞机、柔性可展开的几何可变航天器、旋转机械等。

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

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

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

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

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

这类方法将数据(结构响应)划分成一个个小的时间段，并在每一个时间段内把结构参数看成是时不变的，然后将每一段内识别值用一定的数据处理技术(如曲线拟合)处理得到模态参数随时间变化的规律。它的特点是估计后一段时间的模态参数时没有用到前面各段的数据信息，对参数变化较快的结构为使估计精度提高必须选取很短的数据段。此方法包括现今较为常用的基于状态空间模型的递推的随机子空间辨识法(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，获取被辨识结构的结构动力学响应信号，通过预处理去掉响应信号中的明显不合理因素，并设定采样时间和频率；

步骤2，对各个响应信号进行时频变换，得到时频分布系数，并绘制时频分布图，所述时频变换方法必须具有明确的能量分布物理意义；本发明采用Reassign Gabor Expansion(RGE)和Smooth Pseudo Wigner-Ville Distribution(SPWVD)两种时频变换方法，这两种时频变换方法能量分布物理意义比较明确，算法实施比较简单，计算效率较高。

步骤3，将时频分布系数写成对应的能量分布形式，并重新排列为列向量，如：

Reassign Gabor Expansion的重新排列能量时频分布系数为：

c_i←|c_m，n|²                (1)

其中，c_m，n为Gabor Expansion系数，下标m和n分别表示时间和频率的离散点的标号，c_i为重新排列的能量时频分布系数，i表示列向量元素的标号；

Smooth Pseudo Wigner-Ville Distribution的重新排列能量时频分布系数为：

c_i←|WVD(t_m，f_n)|           (2)

其中，t和f分别表示时间和频率，下标m和n分别表示时间和频率的离散点的标号，WVD(t_m，f_n)为Wigner-Ville Distribution系数，c_i为重新排列的能量时频分布系数，i表示列向量元素的标号；

步骤4，根据各个响应的时频分布图确定将用于辨识的含有各阶时变模态频率的响应对应的时频分布区域；各阶时变模态频率在时频分布图上均表示为一个连续区域，因此从时频分布图中可以直观的判断出时变模态频率的阶数；时频分布图上各阶时变模态频率的颜色深浅表示该阶模态频率对应的响应信号的能量密度，因此通过时频分布图上各阶时变模态频率的颜色深浅可以直观的判断出对应各阶时变模态频率的能量时频分布最高的时频分布区域，也就是对应各阶时变模态频率，选择出颜色最深的图，得到该图上该阶时变模态频率的分布区域；

步骤5，采用合适的时频窗函数将时频分布图中对应各阶时变模态频率的能量时频分布最高的部分分别提取出来；

$c_i^j=h_j(t,f)c_i---\left(3\right)$

其中，i表示列向量元素的标号，标号j表示包含第j阶时变模态频率的能量时频分布区域，h(t，f)表示时频窗函数。

步骤6，以能量时频分布系数

为加权系数，时间t和频率f为坐标，采用加权非线性最小二乘方法对各阶时变模态频率进行估计，最小二乘目标为使下式的函数值最小：

$s=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{c}_i^j{(f_i\left(t_i\right)-g\left(t_i\right))}^2---\left(4\right)$

其中，f_i(t_i)表示时频域内离散的按列向量排开的时间-频率点，g(t_i)表示被辨识的时变模态频率函数；g(t_i)可以设定为任意次数的多项式，对于每种多项式，通过该步骤可以得到使s值最小的g(t_i)，即本发明最终要获得的时变结构模态频率；

步骤7对时变结构模态频率的辨识结果进行误差分析，评估辨识的正确性和准确性，根据误差分析的结果判断选取步骤6中的哪种多项式作为最终结果。

有益效果

本发明具有明确的物理意义，适用于工程应用领域的各种慢变和快变模态参数的时变结构，使用简单方便，并且不受结构规模大小的限制，具有较强的适用性和抗噪声干扰的能力。

附图说明

图1为三自由度弹簧-阻尼器-质量系统；

图2为激励和三个自由度上的加速度响应的时域波形和能量谱密度；

图3为三个自由度上的加速度响应时频分布；

图4为加权非线性最小二乘法辨识结果。

图5为本发明所述方法的流程图。

具体实施方式

下面结合附图，说明本发明的优选实施方式。

本发明提出并实现了一种以能量时频分布系数为加权系数的加权非线性最小二乘时变结构模态频率辨识方法。下面通过一个随机激励的三自由度时变结构实例对本发明进行进一步说明。

图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(5\right)$

其中M(t)为时变的质量矩阵，M(t)＝M₀(1-0.5t)＝diag{0.2，0.1，0.1}(1-0.5t)，M₀为初始时刻的质量矩阵，阻尼矩阵和刚度矩阵为

$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]\×{10}^5$

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

系统响应采用Newmark-β积分格式得到(Newmark-β的积分参数为γ＝0.5，β＝0.1)。激励和三个自由度上的加速度响应的时域波形和能量谱密度如图2所示。采样频率为2000Hz。

图2中(a)为激励的时域波形，(b)(c)(d)为三个自由度上的加速度的时域波形，(e)为激励的能量谱密度，(f)(g)(h)为三个自由度上的加速度的能量谱密度。

本发明所述的基于加权最小二乘时频分布图的时变结构模态频率辨识的具体实现步骤如下：

步骤1，选取如图2(b)(c)(d)所示的三个自由度加速度为辨识所用的响应信号，采样时间为1s，采样频率为2000Hz。

步骤2，对各个响应信号进行时频变换得到时频分布系数，并绘制时频分布图。Reassign Gabor Expansion和Smooth Pseudo Wigner-Ville Distribution表示的时频分布如图3所示。图3中，(a)为第一自由度加速度响应的Reassign GaborExpansion，(b)为第二自由度加速度响应的Reassign Gabor Expansion，(c)为第三自由度加速度响应的Reassign Gabor Expansion，(d)为第一自由度加速度响应的Smooth Pseudo Wigner-Ville Distribution。

步骤3，将时频分布系数写成对应的能量分布形式，并重新排列为列向量。

步骤4，根据各个响应的时频分布图确定将用于辨识的含有各阶时变模态频率的响应和时频分布区域。如图3所示，时频分布图中存在三个黑色的带状区域，从下到上分别表示第一阶、第二阶和第三阶时变模态频率成分的能量时频分布。选择各个自由度响应中能量时频分布最高的时变模态频率，也就是对应各阶时变模态频率，选择出颜色最深的图，即在第三自由度响应的时频分布图中选出第一阶时变模态频率对应的能量时频分布区域；在第1自由度响应的时频分布图中选出第二阶时变模态频率对应的能量时频分布区域；在第一自由度响应的时频分布图中选出第三阶时变模态频率对应的能量时频分布区域。

步骤5，采用矩形时频窗函数将时频分布图中对应各阶时变模态频率的能量时频分布最高的部分分别提取出来。

步骤6，以能量时频分布系数为加权系数，时间t和频率f为坐标，采用加权非线性最小二乘方法对各阶时变模态频率进行辨识。时变模态频率的辨识结果如图4所示。图4中红色圆形图标表示理论计算的结果，绿色菱形图标表示使用本发明的方法进行辨识的结果。

图4中，(a)为以二次多项式为被辨识函数基于Reassigned Gabor Expansion的第一阶时变频率辨识结果，(b)为以二次多项式为被辨识函数基于Reassigned Gabor Expansion的第二阶时变频率辨识结果，(c)为以二次多项式为被辨识函数基于Reassigned Gabor Expansion的第三阶时变频率辨识结果，(d)为以三次多项式为被辨识函数基于Reassigned Gabor Expansion的第一阶时变频率辨识结果，(e)为以三次多项式为被辨识函数基于Reassigned Gabor Expansion的第二阶时变频率辨识结果，(f)为以三次多项式为被辨识函数基于Reassigned GaborExpansion的第三阶时变频率辨识结果，(g)为以二次多项式为被辨识函数基于Smooth Pseudo Wigner-Ville Distribution的第一阶时变频率辨识结果，(h)为以二次多项式为被辨识函数基于Smooth Pseudo Wigner-Ville Distribution的第二阶时变频率辨识结果，(i)为以二次多项式为被辨识函数基于Smooth PseudoWigner-Ville Distribution的第三阶时变频率辨识结果，(j)为以三次多项式为被辨识函数基于Smooth Pseudo Wigner-Ville Distribution的第一阶时变频率辨识结果，(k)为以三次多项式为被辨识函数基于Smooth Pseudo Wigner-VilleDistribution的第二阶时变频率辨识结果，(l)为以三次多项式为被辨识函数基于Smooth Pseudo Wigner-Ville Distribution的第三阶时变频率辨识结果。

步骤7对时变结构模态频率的辨识结果进行误差分析，分别考察辨识结果与时频分布图和理论值之间的绝对和相对方均根误差值。

辨识结果的误差分析如表1所示。

表1误差分析

在表中，p-RMSE表示辨识值与能量时频分布系数之间的方均根误差，t-RMSE表示辨识值与理论值之间的方均根误差，t-rRMSE表示辨识值与理论值之间的相对方均根误差；

如表1所示，选用二次多项式和选用三次多项式辨识结果与时频分布图和理论值之间的绝对和相对方均根误差值均在较低的水平，因此认为辨识结果均是准确的。因此对于本实施例来说，选用二次多项式已经可以满足应用需求。

由此可见本发明能够较好的辨识出时变结构的模态频率，并具有明确的物理意义。由于其仅需要结构的响应信号作为输入，因此，适合对工作状态的时变结构进行模态频率辨识。另一方面，本发明采用的加权最小二乘辨识方法有较强的适用性和抗噪声干扰的能力。

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

Claims (3)

1.一种基于时频分布图的时变结构模态频率辨识方法，其特征在于，包括以下步骤：

步骤1，获取被辨识结构的结构动力学响应信号，通过预处理去掉响应信号中的明显不合理因素，并设定采样时间和频率；

步骤2，对各个响应信号进行时频变换，得到时频分布系数，并绘制时频分布图，所述的时频变换方法可为Reassign Gabor Expansion或Smooth PseudoWigner-Ville Distribution任意一种，这两种时频变换方法具有明确的能量分布物理意义；

步骤3，将时频分布系数写成对应的能量分布形式，并重新排列为列向量ci；

步骤4，根据各个响应的时频分布图确定将用于辨识的含有各阶时变模态频率的响应对应的时频分布区域；各阶时变模态频率在时频分布图上均表示为一个连续区域，因此从时频分布图中可以直观的判断出时变模态频率的阶数；时频分布图上各阶时变模态频率的颜色深浅表示该阶模态频率对应的响应信号的能量密度，因此通过时频分布图上各阶时变模态频率的颜色深浅可以直观的判断出对应各阶时变模态频率的能量时频分布最高的时频分布区域，也就是对应各阶时变模态频率，选择出颜色最深的图，得到该图上该阶时变模态频率的分布区域；

步骤5，采用合适的时频窗函数将时频分布图中对应各阶时变模态频率的能量时频分布最高的部分分别提取出来；

$c_i^j=h_j(t,f)c_i---\left(1\right)$

其中，i表示列向量元素的标号，标号j表示包含第j阶时变模态频率的能量时频分布区域，h(t，f)表示时频窗函数。

步骤6，以能量时频分布系数

为加权系数，时间t和频率f为坐标，采用加权非线性最小二乘方法对各阶时变模态频率进行估计，最小二乘目标为使下式的函数值最小：

$s=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{c}_i^j{(f_i\left(t_i\right)-g\left(t_i\right))}^2---\left(2\right)$

其中，f_i(t_i)表示时频域内离散的按列向量排开的时间-频率点，g(t_i)表示被辨识的时变模态频率函数；g(t_i)可以设定为任意次数的多项式，对于每种多项式，通过该步骤可以得到使s值最小的g(t_i)，即本发明最终要获得的时变结构模态频率；

步骤7对时变结构模态频率的辨识结果进行误差分析，评估辨识的正确性和准确性，根据误差分析的结果判断选取步骤6中的哪种多项式作为最终结果。

2.根据权利要求1所述的一种基于时频分布图的时变结构模态频率辨识方法，其特征在于，步骤3中当时频变换方法为Reassign Gabor Expansion方法时所述的重新排列的能量时频分布系数c_i为：

c_i←|c_m，n|²

其中，c_m，n为Gabor Expansion系数，下标m和n分别表示时间和频率的离散点的标号，i表示列向量元素的标号。

3.根据权利要求1所述的一种基于时频分布图的时变结构模态频率辨识方法，其特征在于，步骤3中当时频变换方法为Smooth Pseudo Wigner-VilleDistribution方法时所述的重新排列的能量时频分布系数c_i为：

c_i←|WVD(t_m，f_n)|

其中，t和f分别表示时间和频率，下标m和n分别表示时间和频率的离散点的标号，WVD(t_m，f_n)为Wigner-Ville Distribution系数，i表示列向量元素的标号。

Images