中国振动工程学会模态分析高级研修班讲课资料第五章.pptx

Institute of Vibration Engineering 振动工程研究所振动工程研究所 1多输入多输出系统的模态参数识别多输入多输出系统的模态参数识别张永强 高级工程师靖江泰斯特电子有限公司西北工业大学 振动工程研究所Institute of Vibration Engineering 振动工程研究所振动工程研究所 2 概述 单点激励的不足 激励能量不够，且传递过程中损耗过大； 离激励点较远的地方响应信号较弱，信噪比低； 较大激励会造成局部响应过大，产生非线性现象 若激励处于节点位置，系统变成不可控和不可观的； 模态密集时辨识能力较弱； 多输入多输出方法 时域和频域两种方法；Institute of Vibration Engineering 振动工程研究所振动工程研究所 3 多输入多输出频响函数估计 系统的多输入多输出模型 输入与输出Tpfff)()()(21FTLxxx)()()(21XTpmmm)()()(21MTpnnn)()()(21NTpuuu)()()(21UTpvvv)()()(21V实测输入实测输出输入噪声输出噪声真实输入真实输出Institute of Vibration Engineering 振动工程研究所振动工程研究所 4 系统模型 干扰影响 无干扰时 有干扰时 误差（总体误差） 测量误差 信号处理误差 非线性误差)()()()(1jPjijiifHxv)()()( )()()()( )()()(11ijpjijijjPjijiiiEfHnmfHnvxInstitute of Vibration Engineering 振动工程研究所振动工程研究所 5 估计输出噪声估计模型 估计模型 无输入噪声 输出噪声与输入信号无关 估计 右乘F的共轭转置 再求数学期望，得 输出噪声与输入信号不相关时 频响函数估计1HNHFXF输入 向量 H频响 函数 矩阵X输出 向量 N系统 噪声 向量1HHFNFFFXFGHGG0NFG11FFXFGGHGXF输入输出互功率谱密度矩阵 GFF输入自功率谱密度矩阵 GNF输出误差和输入的互功率谱矩阵 Institute of Vibration Engineering 振动工程研究所振动工程研究所 6 特点 是一种欠估计 对输入噪声比较敏感 输入噪声较大时，精度受影响，在共振点附近更是如此 要对输入自谱矩阵求逆，计算量大，且矩阵奇异易导致求逆失败 GFF奇异的几方面原因 某个输入谱为零时 两个或更多输入信号完全相关时 数值计算中的问题:矩阵病态等01HInstitute of Vibration Engineering 振动工程研究所振动工程研究所 7 检验判别工具),(),(),(22iijjjiXXFFXFijGGG若某一输入信号与若某一输入信号与输出信号常相干函输出信号常相干函数等于数等于1，则表示该，则表示该输入信号与该输出输入信号与该输出信号完全相关，该输信号完全相关，该输出完全由该输入产生出完全由该输入产生若两输入信号常相干函若两输入信号常相干函数等于数等于1，则表示两个，则表示两个输入信号完全相关输入信号完全相关常相干函数（表明两信号的因果关系）Institute of Vibration Engineering 振动工程研究所振动工程研究所 8偏相干函数 消除其它输入信号的潜在贡献后，输入与输入、输入与输出、输出与输出之间的相干函数； 如果输入偏相干函数为1，表明两个输入力是相关的。重相干函数 描述某个输出信号与所有已知输入信号之间因果关系； 重相干函数等于1，表示输出xi全部由输入信号f1、 f1、 fp引起； 重相干函数等于0，表示输出xi全部由未知噪声引起的。iiiiiijiixxnnxxPjfxijFxGGGGHmin12)(1Institute of Vibration Engineering 振动工程研究所振动工程研究所 9 估计输入噪声估计模型 估计模型 只有输入噪声 假设输入噪声与输出信号不相关 估计 右乘X的共轭转置 再求数学期望 输入噪声与输出信号不相关时 频响函数估计2H)(MFHXF输入向量 H频响函数矩阵X输出向量 M系统噪声向量2HHX0XMG12FXXXGGHGXF输入输出互功率谱密度矩阵 GFF输入自功率谱密度矩阵 GNF输出误差和输入的互功率谱矩阵 Institute of Vibration Engineering 振动工程研究所振动工程研究所 10 特点 有唯一解的条件是GFX的逆矩阵存在 当激励力数P比响应测点数L小时， GFX的逆不存在 此时可利用最小二乘解，利用GFX的伪逆矩阵，求解 只考虑输入噪声的影响，对输出噪声比较敏感 是一种过估计，即有 在共振点附近， 估计较 估计有较高的精度 在反共振点附近，情况相反2H102HH2H1HInstitute of Vibration Engineering 振动工程研究所振动工程研究所 11 输入输出噪声估计模型 同时考虑输入输出噪声 估计 取 和 的算术平均值 或取其加权平均 估计 取 和 的几何平均值3H2H1H)(21213HHH1213HHH4H2H1H214HHH Institute of Vibration Engineering 振动工程研究所振动工程研究所 12 估计 利用最小二乘法原理，极小化误差矩阵的方法 圆盘结构三种估计对比试验 圆板放置在泡沫塑料衬垫上 采用随机激励 故意造成一些泄漏 人为施加一些噪声vHInstitute of Vibration Engineering 振动工程研究所振动工程研究所 13特征系统实现算法（ERA）最小二乘复频域法 （ PolyMax ）Institute of Vibration Engineering 振动工程研究所振动工程研究所 14 特征系统实现算法（ERA）首先由美国国家航空与宇航局 ( NASA ）所属的Langley 研究中心于1984 年提出。它是一种属于多输入多输出的时域整体模态参数辨识方法。它移植了自动控制理论中的最小实现理论，利用脉冲响应数据，采用奇异值分解的方法，求得系统的特征值与特征向量，从而求得模态参数。该方法于1984年提出后，当年即在美国伽利略航天器的模态分析中应用，次年又在航天飞机机载巨型太阳能帆板的太空模态试验中应用，均取得良好的效果，该方法有最佳的精度，因此是目前比较先进的一种时域参数辨识方法。特征系统实现算法（ERA 法）Institute of Vibration Engineering 振动工程研究所振动工程研究所 15 特征系统实现算法是利用实测的自由响应（脉冲响应函数，相关函数），运用奇异值分解方法，确定系统的阶次和状态方程中的系统矩阵A 、输人矩阵B 和输出矩阵C ，进而求解系统矩阵A 的特征值问题，求得极点与留数，从而确定系统的模态参数。当矩阵A 、B 、C 的阶次为最小时，即为最小实现。此时系统是可控的，又是可观的。基本思路Institute of Vibration Engineering 振动工程研究所振动工程研究所 16 （1）系统的状态方程描述 对一个N 自由度的线性系统，若在P个点激励，在L 个点上测量响应，可用下列状态方程描述： 式中：K 为采样点序号；X ( K ）是在K时刻系统的状态向量，( 2N xl ) ； 为采样间隔时间；Y( K ）是在K时刻的实测响应向量，L x1;F ( K ）是在K时刻系统的输人向量，Pxl ; A 为系统矩阵，2Nx2N ; B为输入矩阵，又称控制矩阵，2Nx P；C 为输出矩阵，又称观测矩阵，L x2N 。 对一线性定常系统，自由响应可用脉冲响应来代替。因此自由响应的最小实现问题常用脉冲响应的最小实现问题来代替。Institute of Vibration Engineering 振动工程研究所振动工程研究所 17（2）脉冲响应矩阵的建立 系统的脉冲响应可由实侧传递函数的拉氏逆变换求得。对各点的脉冲响应函数h ( t ）进行离散采样后，便可得离散的脉冲响应函数序列h ( K ) , K = 1 , 2 , 。设采样的时间间隔为 。在K 时刻，各测量点的脉冲响应可构成下列脉冲响应矩阵： 式中：hij。为j 点激励、i 点测量的脉冲响应函数；L 、P分别为测量点与激励点的数目。 脉冲响应的最小实现问题是已知 及求矩阵A 、B 、C ，并使三重矩阵 A ，B，C 的阶次最小。在求得系统矩阵A 后，再由其特征值与特征向量确定模态参数。Institute of Vibration Engineering 振动工程研究所振动工程研究所 18（3）构成Hankel矩阵 脉冲响应的最小实现一般是从生成Hankel 分块矩阵开始。Hankel 矩阵有如下形式：Institute of Vibration Engineering 振动工程研究所振动工程研究所 19对线性定常系统脉冲响应与矩阵A 、B、C 之间有如下关系：（4）、脉冲响应与三重矩阵【A 、B 、C 】之间的关系 对上式递推，可得Institute of Vibration Engineering 振动工程研究所振动工程研究所 20 继续递推并代入式，可得 式中；P 矩阵称为可观性矩阵；Q 矩阵称为可控性矩阵；、则称为可观、可控性指数，且有 由此即可导出特征系统实现算法的主要计算公式。 Institute of Vibration Engineering 振动工程研究所振动工程研究所 21（5）、特征系统实现算法 由式，当K=1时，有 显然亦有 对 进行奇异值分解， 式中：U 为左奇异向量；V 为右奇异向量；为奇异值矩阵, , U、V是正交归一化矩阵， i称为奇异值，并且有123r。 矩阵的秩即为 系统的阶次。可由不为零的奇异值的个数来确定。NNR22Institute of Vibration Engineering 振动工程研究所振动工程研究所 22U ，V ，和A,B,C的关系其中Institute of Vibration Engineering 振动工程研究所振动工程研究所 23 因此系统的状态方程规范型可写为 由式可见，A矩阵的阶次取决于的阶次，而矩阵 。因此，尽管的阶次很高（Lx P),经过奇异值分解后，属于2N阶。因为 由此可见，系统矩阵A为2N阶方阵。相应的状态矢量X(K)的阶次必为2N阶，它是描述2N阶系统的最小阶次，因此是最小实现。NNR22Institute of Vibration Engineering 振动工程研究所振动工程研究所 24 经上述推导，可以对 做奇异值分解的含意，理解如下； 1）从逼近理论来看， 是 所在子空间的最佳逼近。 2）从信号处理角度来看，用 代替 相当于对数据进行一次维纳滤波。被滤掉的是对应于奇异值为零的那些与输入、输出无关的随机噪声。因此状态方程无需再为噪声提供出口，无需再进行扩阶。3）以最少的参数、最小的阶次来描述系统的特征和进行运算，从而减小了运算量。4）提高了算法的抗噪声干扰能力，避免了在模态转化过程中产生计算误差，即出现虚假模态。)0(HVUT)0(HVUT)0(HInstitute of Vibration Engineering 振动工程研究所振动工程研究所 25（6）、模态参数辨识系统的模态参数可由系统矩阵A 的特征值及特征向量来确定。系统矩系统矩阵A可由下式确定，对矩阵A进行特征值分解，求出特征值矩阵，然后求出特征向量矩阵 式中：z 为特征值矩阵；为特征向量矩阵， 由此便可确定系统的模态振型， Institute of Vibration Engineering 振动工程研究所振动工程研究所 26 然后由下列公式可求的系统的模态频率r及模态阻尼r;矩阵A的特征值与系统特征值之间有如下关系；式中Z r及r均为复数，可写成Institute of Vibration Engineering 振动工程研究所振动工程研究所 27 PolyMax模态识别方法,属于多自由度时域识别法,也称作多参考点最小二乘复频域法( Polyreference least squares complex frequency domain method)， 是最小二乘复频域法(LSCF)的多输入形式，是一种对极点和模态参预因子进行整体估计的多自由度法,一般首先通过实验建立稳态图，以判定真实的模态频率、阻尼和参预因子；建立可以线性化的直交矩阵分式模型，然后基于正则方程缩减最小二乘问题,得到压缩正则方程,于是模态参数可以通过求解最小二乘问题得到。该方法集合了多参考点法和LSCF方法的优点,可以得出非常清晰的稳态图,并且密集空间可以被分离出来,尤其在模态较密集的系统(动力总成系统)，或者FRF数据受到严重噪声污染的情况下仍可以建立清晰的稳态图，识别出高度密集的模态，对每一个模态的频率、阻尼和振型都有很好的识别精度,是国际最新发展并流行的基于传递函数的模态分析方法。 Institute of Vibration Engineering 振动工程研究所振动工程研究所 28 （1）建立频率响应函数模型 多参考点最小二乘复频域识别技术（PRLSCF或PolyMAX）要以频响函数矩阵作为识别的初始数据,其数学模型采用右矩阵分式模型来描述。.在频域中，系统输出 （ ， 其中 为输出点数）和全部输入的关系可用右矩阵分式模型（RMFD）来描述，右矩阵分式模型的表达式为 （1） 式中： 理论频响函数的第行，是输入点数，即激励 数； 分子多项式行向量； 分母多项式矩阵。 且 和 可以表示成如下形式：0,2 , 1Noo0N 1DUHoo iNloCH iNloCU iiNNoCD oU oDInstitute of Vibration Engineering 振动工程研究所振动工程研究所 29 （ ） （2） （3）式中： N多项式阶次 其中分母系数矩阵 和分子系数行向量 是待估计的参数。所有这些系数合并为一个矩阵。 （4） 其中 （5） NrorroBZU00,2 , 1No NrrroAZD0iiNNrRAiNlorRBTTTNT1iNNoNoooRBBB110iiNNNNoRAAA110Institute of Vibration Engineering 振动工程研究所振动工程研究所 30式（2）和式（3）中出现的多项式基函数，一般地，有以下两种选择：.对于连续时域模型，可取为 （6）式中： 缩放因子，用来提高方程的数值状况。 .对于离散时域模型，可取为 （7）式中： 采样周期。 通常采用离散时域模型。 jsriZ21iNs srTjreZsTInstitute of Vibration Engineering 振动工程研究所振动工程研究所 31（2）参数的线性化 通过试验测量出的频率响应函数矩阵 ，用 表示实测频响矩阵的第 行， ，那么关于参数矩阵 的非线性最小二乘（NLS）目标函数可表示为 （8） 式中： 矩阵的复共扼转置； 矩阵的迹，即矩阵的主对角元素之和。 通过对式（8）求极小值，便可以得到频率响应函数矩阵的右分式矩阵模型各系数的估计值，即 矩阵的估计值。 式（8）中的加权非线性最小二乘误差函数被定义： （9 ） ioNNfCHiNfoCH1ofoNfNo, 2 , 1, 2 , 1 ofNoNffNLSoHfNLSonlstr11,H trfofofofofofofofNLSoHDUWHHW,1Institute of Vibration Engineering 振动工程研究所振动工程研究所 32 上式中 是一个加权函数。一般地，为了提高估计的质量，我们采用 （10） 式中： 方差，可用相关函数求取。也可使用公式 （11） 来做加权函数的。这两种加权函数都考虑了测量频响函数数据的好坏：测得频响的方差越小，对目标函数的贡献越大。 非线性误差函数可以经过一个近似的处理为一个线性的问题。实际上，通过对 右乘 ，则可以得到一个关于参数为线性的方程，此加权线性最小二乘（LS）方程误差 为ofWvarofofofHWH var 1varofofWH,NLSof,fD,LSofInstitute of Vibration Engineering 振动工程研究所振动工程研究所 33 （12） 这样式（12）关于参数为线性，将所有频率列， ，它可 用矩阵形式来表示 （13） 10,LSNLSofoffofofofoffNofrforrforrDWUDHDWzBzHA1,2,ffN 1,fLSoooLSooooLSoNXYJ Institute of Vibration Engineering 振动工程研究所振动工程研究所 34其中： （14） （15）式中， Kronecker积。1111111,fffffooNNNooNoNNNNWzzzXCWzzz11111111,fifffffooNoNNNooNoNNNNoNWzzzHYCWzzzHInstitute of Vibration Engineering 振动工程研究所振动工程研究所 35（3）缩减标准方程 加权线性最小二乘估计表达式为 （16）式中：同时，目标函数（16）等价于 （17） 11ooNHLSLSLSoooNoooTToToootrRStrST 11ReNNHoooRXXR11ReiNNNHoooSX YR 11ReiiNNNNHoooTY YR ReTHLStrJ JInstitute of Vibration Engineering 振动工程研究所振动工程研究所 36式中， 是Jacobian矩阵，被如下定义 （18）为使 值最小，将 对系数矩阵 和 求导，并令其为零 （19） （20）J1122100000ofoiooN NNNNNNXYXYJCXY LS LSo 20LSooooRS1,2,ooN 120oNTLSooooSTInstitute of Vibration Engineering 振动工程研究所振动工程研究所 37由式（19）得到 ，把它代入式（20）得 （21）其中，由式（19）和（20）得到标准方程，经过整理，此标准方程的表达式为 （22）1oooR S 1120oNToooooTS R SM11112oiiNNNNNToooooMTS R SR 111222121000022Re000oooooHNNNNTTTNooRSRSJJRSSSST Institute of Vibration Engineering 振动工程研究所振动工程研究所 38 式（21）即为“缩减”标准方程，其中矩阵 维数为 比标准方程式（22）中的 的维数 要小的多。 （4）求解缩减标准方程 通过求解“缩减”标准方程，便可得到分母系数矩阵。根据线性方程组的求解理论，先对系数矩阵施加一个约束。假如，设定系数矩阵中的一个系数矩阵块等于正则常数矩阵（例如设系数矩阵的最后一个矩阵块）,在这种前提下，缩减标准方程变为 （23） 其中M11iiNNNNReHJ J 11oioiNNNNNNA XB1:,1:iiAMN NN N1:1:1iiiBMN N N NNN Institute of Vibration Engineering 振动工程研究所振动工程研究所 39 系数矩阵 的最小二乘估计为 （24） 一旦求得了 ，那么通过 就可得到所有的分子系数 这种方法考虑了标准方程的结构特性，比直接求解方程（22）要快得多。确定了分母系数矩阵 后，通过求解 的伴随矩阵的特征值和特征向量，这样就可以得到了系统的极点和相应的模态参与因子。方程如下 （25）1iLSNXXA BILS1oooR S ,LS o01210000000.000TTTTNNIVVIAAAA Institute of Vibration Engineering 振动工程研究所振动工程研究所 40 上式中， ，矩阵 的最后 行就是模态参与因子；对角阵 的角元记录为 由不稳定的数学极点和稳定的物理结构点两部分组成。记稳定的物理结构极点为 ，通过对这些物理结构极进行转换，便可得出结构的固有频率 和模态阻尼比 ；关系式如下 或 （26） （5）计算频率点和阻尼比点 根据信号与系统基本理论中对系统稳定性的描述：系统的全部极点落于 域左半平面（不包括虚轴），且满足有界输入有界输出原则，系统是稳定的。复特征矩阵 中的复特征值总是以共轭对的形式出现，同时也包含实数（虚轴上），在求解频率点 和阻尼比点 时，对于每个共轭对只取其中一个进行分析，且不考虑实数。,ooN N N NVCViN(1,2,)ioiN Nr sTre rr*,rrrri *2,1rrrrrri siiInstitute of Vibration Engineering 振动工程研究所振动工程研究所 41 复特征矩阵 中的对角元 ，由式（26）， 用 描述，则 （27） （28） （29） 所以 （30） （31）i sTie iReImiii ReImcossiniisi si si siTTiTTiiisisieeeeTiT ReImi sTiiie ImarctanReiisiT1lniisT Im1arctanReiisiTInstitute of Vibration Engineering 振动工程研究所振动工程研究所 42 由此可求得频率 和阻尼比 （32） 在求得的频率 和阻尼比 包含有结构的固有频率 和模态阻尼比 ，因此，必须对所有求得的 和 进行有效的分析和选取，以确定系统真实的固有频率和阻尼比。建立稳态图就是一种行之有效的方法。ii22iiiiii iirriiInstitute of Vibration Engineering 振动工程研究所振动工程研究所 43 （6）建立稳态图 在模态分析中，稳态图是帮助实验者分离结构物理极点和数学极点的一个有力工具，如图所示。通过逐渐增大多项式的阶次 ，且进行相应的重复性分析计算可以建立起稳态图。 模态分析的稳态图 N1.002.05e+3LinearHz79.0e-616.7e-3Amplitude(g/N)ovovoosvoovovovosvoovovooossoovvvvooosvoovo osvooo svoosovosooov ssooovoosoooovssoooovvosoooov ssoooosvoso vovos vso ooo vsosvovvos vsooooo vvoosooosvsvsooooovoovooovsvv vsoooooovsvo voov vvss vv sooovvooovvovvvsv v vv soooosvooosovsv vsv sv vo ooov vvoov ovvvv vsv vv sooooovovovooovvoov vvvoo o oo voovv sovvovvvv vvvooo oo sovovo oovvvvvv vsvoo ooosoovo vvovvvvvvs sv vsoooovooovv vvo vvvvsvs vs svooooooss vos soovovvvvvs vv voov oovvsovvovvvs ovvvso vv ss vs1011121314151617181920212223242526272829303132Institute of Vibration Engineering 振动工程研究所振动工程研究所 44 对于不同的 值，由式（25）得到相应的极点 ，极点再通过转换得出相应的频率 和阻尼比 。稳态图的建立方法为：横坐标是频率点 ，纵坐标是对应的多项式的阶次值 。一般而言，多项式的阶次 取值越大，在稳态图上所反映的极点信息越全面，但与此同时，也带入了许多虚假的极点信息，这样会给极点分离带来困难，经验表明， 值一般取为50左右比较理想。Ni sTie iiNNNInstitute of Vibration Engineering 振动工程研究所振动工程研究所 45