基于响应传递率的非白随机激励仅输出结构模态参数辨识

周思达，刘 莉，杨 武，马志赛

（北京理工大学宇航学院飞行器动力学与控制教育部重点实验室，北京 100081）

基于响应传递率的非白随机激励仅输出结构模态参数辨识

周思达，刘 莉，杨 武，马志赛

（北京理工大学宇航学院飞行器动力学与控制教育部重点实验室，北京 100081）

在对工作状态的航天器等结构进行仅输出模态参数辨识时，激励非白造成了无法剔除谐波引起的虚假模态。为解决此问题，该研究对最新的基于结构响应传递率的模态参数辨识方法进行了改进，提出了一种改进的基于多参考点响应传递率的仅输出模态参数辨识方法。首先，推导了多参考点的结构响应传递率表达式，建立了其左矩阵分式多项式的参数化模型，进而给出了辨识问题的最小二乘估计，利用正则方程Jacobi矩阵的分块性质对最小二乘问题矩阵形式完成了缩减，降低了计算量；然后，通过高维伴随矩阵方法解决了矩阵多项式的特征值求解问题，即多参考多输出的模态参数求解问题，以及通过矩阵伪逆解决了现有方法中载荷工况数与响应点数和参考点数的约束问题；最后，通过两个数值算例对所提方法进行了验证，辨识结果表明：辨识方法能够很好的辨识出结构的模态频率、阻尼比和模态振型，且能够很好的避免激励中含有的谐波分量对辨识结果造成的影响，解决了传统仅输出模态参数辨识中激励非白对辨识结果造成的影响问题。

模态参数辨识；仅输出；响应传递率；非白随机激励

随着航天科技的快速发展，以及相关领域的迫切需要，对大型空间结构的需求变得愈加迫切。由于航天器对结构重量的限制，其结构往往具有大柔性。这样的大柔性造成的结构振动，会对整器姿态控制、关键仪器性能与力学环境、结构健康等方面产生不可忽视的影响。

通常，对于分布参数的结构振动控制，从整器姿态控制的观测量中滤除柔性振动量，以及结构健康监测等方面，都涉及到结构的模态参数，甚至基于模态参数。辨识得到的模态为以上应用提供了有效的低阶模型。获取模态参数的途径有很多，然而不同于理论分析和数值计算，以及地面的模态实验手段，通过测量在轨运行航天器结构动力学响应，从而通过辨识的手段获取航天器结构模态参数，可以很好的考虑航天器在轨运行工作状态中受环境因素的真实影响。

然而，对在轨航天器进行模态参数辨识也面临巨大的挑战。在轨航天器承受来自环境的随机激励，以及自身机构运动和控制元件的驱动力，这些激励是无法准确测量的，因此在轨航天器结构模态参数辨识是仅输出的。这样针对运行结构的仅输出的模态分析被称为工作模态分析，已经在工程［1］、建筑工程［2－3］、航空工程［3］、机械工程［4］等领域得到了一定应用。现有方法多基于白噪声假设，即假设作用在结构上的激励为白噪声，那么结构响应的功率谱函数就只是模态参数（除模态参与因子）和激励幅值的函数，而后者由于白噪声假设，其功率函数为关于频率的常函数，从而可仅由结构响应辨识出结构的模态参数。

由于航天器工作的特殊环境，以及动量轮等姿控器件、液体晃动等产生的谐波激励，造成了作用在航天器结构上的激励为非白噪声，因此现有基于白噪声假设的仅输出模态参数辨识方法无法满足要求。

旋转机械的工作模态分析也遇到了相似的非白噪声激励问题。过去，研究者利用预置陷波滤波器的方法滤除谐波影响，取得了一定的效果［5］。此类方法需要预先知道激励中谐波激励的频率和水平，但是在很多情况下是无法满足的。最近，有学者提出了基于响应传递率的技术处模态参数辨识方法，以克服工作模态分析中激励非白对辨识结果的影响。Devriendt等［6］第一次提出了基于响应传递率模态参数辨识的概念，之后，他们又在多输出、辨识算法、求解方法方面做了一些改进［7－9］。

为解决激励非白对结构模态参数辨识的影响问题，本文对现有基于响应传递率的模态参数辨识方法进行改进，提出了一种利用矩阵伪逆的频域仅输出模态参数辨识方法。

1 多参考的结构响应传递率

由于传递率定义为结构不同位置力或响应量频谱的比值，表达的是结构不同位置力或响应量的传递关系，通常用来评价结构隔振的效果［10－11］。然而结构响应传递率参数化表达式中并不显含系统极点信息，因此，长期以来并未如结构系统传递函数、功率谱函数等显含系统极点信息的函数一样，被广泛的应用于结构模态参数辨识。

下面通过结构系统的传递函数、不同工况下的结构响应来说明结构响应传递率的性质和伪传递率的定义。

若定义多输入多输出（MIMO）的结构系统传递函数为H（s）∈CNo，Ni（其中No和Ni分别为输出和输入数），那么结构响应为：

式中，X（s）∈CNo×1为响应（输出），F（s）∈CNi×1为激励（输入）。

为方便定义结构响应传递率，将响应分为参考响应和非参考响应，即X（s）＝［Xr（s） Xp（s）］T，其中r 和q分别表示参考点集合和非参考点集合。那么，式（1）可写为：

式中，Xr（s）∈CNr×1，Xq（s）∈C（No－Nr）×1，Nr为参考点数，Hr（s）∈CNr×Ni，Hq（s）∈C（No－Nr）×Ni。

根据上式，文献［12］定义了多参考点的结构响应传递率：

将式（2）代入式（3）得：

若激励F（s）不为0，则有：

如果Nr＝Ni，即输入点个数等于参考点个数，式（5）可以写为：

如果Nr＞Ni，则可以用伪逆来代替矩阵Hr（s）的逆，则式（5）可以写为：

如果Nr＜Ni，则无法利用如式（6）或式（7）中Hr（s）的逆或伪逆来定义多参考点的结构响应传递率。

类似于单参考点的结构响应传递率，除在极点附近，不同工况的结构响应传递率是不同的，多参考点的结构响应传递率也与不同激励点上作用的激励相对幅值有关，即与工况有关。

假设在Ni个输入点上作用种工况，用l指出不同的工况Fl（s）。对于Ni个输入点，存在Ni组独立（线性无关）的激励分布组合出工况Fl（s），如：

式中，Fl，j为第l个工况的第j组独立的激励分布，βl，j（s）为第l个工况的第j组独立的激励在s域上的相对幅值。将式（8）写为矩阵形式为：

式中，Fl∈CNi×Ni为激励分布矩阵，βl（s）∈CNi×1为激励在s域上的相对幅值。

根据式（2），式（9）可写为：

多参考点的系统传递函数也可写为公分母形式的分式多项式参数化模型，如H（s）＝N（s）／D（s），因此，式（11）可改写成：

如式（12）所示，多参考点传递函数矩阵的逆与非参考点传递函数矩阵相乘导致结构响应传递率中的极点信息被消去。因此，要用响应传递率进行模态参数辨识，还需要利用关于传递函数和传递率的更多性质，具体方法将在下文给出。

2 基于多参考伪传递率的模态参数辨识

2.1 传递率的参数化模型

如式（3），传递率与系统传递函数的形式一致，表示频域的两个向量的传递关系，因此可采用矩阵分式模型的形式对多参考点的传递率进行参数化建模。又如式（3）定义，传递率矩阵T（s）为（No－Nr）×Nr的矩阵，因此采用左矩阵分式多项式［13］对第l个工况下的多参考点的传递率建立参数化模型：

式中，

式中，Al（Ω）为（No－Nr）×（No－Nr）矩阵，Bl（Ω）为（No－Nr）×Nr矩阵，nA和nB为对应多项式阶数为频率基函数。为提高参数估计的矩阵数值条件，以z域多项式为频率基函数，即＝e－jωkiΔT。

2.2 从传递率的参数化模型到模态参数

为建立参数化的传递率模型与传递率性质的关系，即如何利用参数化模型中的待估参数计算模态参数，下面利用了多输入多输出传递函数矩阵秩为1的性质完成推导。

由于多输入多输出传递函数矩阵在极点处秩为1，且在极点附近所有列与模态振型成正比［11］，因此，式（18）在极点附近有：

利用Nl种工况，结合式（13）和（19），极点λm和模态振型φm可由以下的矩阵多项式的广义特征值求解得到。

矩阵多项式的广义特征值需展成如下形式按照普通矩阵特征值问题进行求解：

式中，n＝max（nA，nB），Γ为方阵，称为式（20）的高维伴随矩阵，其中的子矩阵须为方阵。文献［9］为满足子矩阵为方阵，给出了约束Nl＝No／（No－Nr）。然而，这个约束在实际应用中是很难满足的。因此，本文利用矩阵伪逆，实现广义矩阵的方阵化。即令：

2.3 模型参数的估计

如前述，只要得到不同工况下的待估参数矩阵Al，i和Bl，i，就能够利用式（22）计算出结构的模态参数。本节将给出待估参数的最小二乘估计方法。

根据传递率的定义和传递率的参数化模型，此参数估计问题的最小二乘花费函数可表示为：

式中，al，m和bl，m分别为Al（Ω）和Bl（Ω）第m行的行向量，Nf为频带内的频率采样点个数，最小二乘的方程误差为：

将式（14）代入式（25），可得：

那么，最小二乘的花费函数也可以写成相应的矩阵形式：

那么，最小二乘的正则方程为：

为减小最小二乘求解的计算量，将原最小二乘问题转化为减缩形式。对原正则方程（式（29））进行处理：

由于若传递率参数化模型中矩阵Al（Ω）和Bl（Ω）同乘一个常数，而待估的传递率T＾l（Ω）并不会变化，因此如式（25）的传递率参数化模型是参数冗余的，当它被用来做最小二乘估计时导致了正则方程（式（32））仅有平凡解或无解。那么，为了消除参数化模型中的参数冗余，需对待估参数矩阵Θl中的某一个子块进行约束。本文将Ψl中Al，nA约束为单位矩阵，那么，约束后的正则方程为：

解方程式（33）得到Ψl，进而通过（31）计算出Ξl，至此待估参数Al，j和Bl，j就通过最小二乘估计得到了。

3 算例一

本节利用一个4自由度的弹簧－阻尼－质量系统对本文研究的方法进行验证。4自由度的弹簧－阻尼－质量系统如图1所示。

图1 四自由度的弹簧－阻尼－质量系统Fig．1 4－DOF system

作用在如图1所示的激励为理想白噪声，即激励的频谱是频率的常函数。工况数Nl＝2，两种工况的相对幅值如图2所示。

图2 两种工况激励的相对幅值Fig．2 Relative Amplitude of excitation for the two load cases

图3 稳定图（算例一）Fig．3 Stabilization diagram（Example 1）

响应参考点数Nr＝2，分别为响应点1和2。令nA＝nB，并从4取值到20构造稳定图，如图3所示（o为不稳定，f、d和v分别为频率、阻尼比或振型稳定，s为频率、阻尼比和振型都稳定）。

由如图3所示的稳定图可以看出4自由度系统的4阶模态（连续的“s”表示），辨识得到的模态频率和阻尼比与理论值对比如表1所示。

表1 辨识得到的模态频率和阻尼比与理论值对比（算例一）Tab.1 Comparison between estimated and theoretical modal frequency and damping ratio（Example 1）

辨识得到的模态振型与理论值对比如表2所示。（由于仅输出模态参数辨识无法得到模态缩放因子，因此为方便对比，辨识值和理论值都进行最大值归一化）。

表2 辨识得到的模态振型与理论值对比Tab.2 Comparison between estimated and theoretical mode shapes

由表1和表2可以看出，在没有噪声干扰的情况下，本文给出的辨识方法能够十分精确的辨识出结构的模态频率和阻尼比。

4 算例二

为验证方法对实际结构模态参数辨识的有效性和正确性，本节给出一个有限元数值算例。

4.1 模型

有限元模型模拟一个在空间飞行的火箭结构，边界条件为自由－自由。有限元模型如图4所示，结构采用梁单元模拟，锥形整流罩部分采用变截面梁单元实现。结构总长20 m，直径1．4 m，分为20个单元，共21个节点。模型阻尼矩阵由Rayleigh阻尼公式得到，其中α＝1．0，β＝0．000 25。

方便问题描述和结果分析，本算例仅研究横向一个面内的响应和模态特性，即仅考虑如图4所示的y方向的激励以及y方向的平动和绕z方向的转动自由度。

图4 火箭有限元模型（三维显示）Fig．4 Finite element model of a rocket（3D display）

4.2 工况与响应仿真分析

假设作用在结构上的有两种工况：一种是作用在整流罩部分节点上的随机载荷，模拟气动力的随机波动部分；另一种是作用在火箭尾部上的发动机摆动在y方向的随机分量。通常液体火箭的液体燃料、发动机摆动等导致了周期载荷，为考虑这样的周期载荷，本算例模拟其极端情况，即考虑两个单频率分量的谐波作用在结构上，分别是6 Hz和46．5 Hz，幅值为随机载荷方差的103（仅为考验本文的辨识算法，实际情况不会达到此量级）。利用数值积分得到结构在两种工况下的加速度响应，从而获取结构响应的频谱信息X（Ω）。两个工况下1号点频率响应如图5所示。

图5 两个工况下1号点频率响应（上：工况1；下：工况2）Fig．5 Frequency responses of Point 1 （upper：load case 1；lower：load case 2）

如图5所示，结构的频率响应中在6 Hz和46．5 Hz处存在由于谐波引起的峰（黑色虚线处）。若采用传统的基于功率谱函数的仅输出辨识方法进行模态参数辨识将无可避免的在这两个峰值处产生虚假模态。

4.3 辨识结果

通过本文给出的基于结构响应传递率的模态参数辨识方法，仅利用结构的响应信息进行模态参数辨识，其中，辨识带宽选为0－120 Hz，21个点的响应皆用于辨识（No＝21），响应参考点数Nr＝2，分别为响应点3和响应点4。令nA＝nB，并从8取值到20构造稳定图（如图6所示）。

如图6所示，稳定图中的4列连续的稳定极点指出了带宽内的四阶模态。模态频率和阻尼比辨识结果如表3所示。

图6 稳定图（算例二）Fig．6 Stabilization diagram（Example 2）

表3 辨识得到的模态频率和阻尼比与理论值对比（算例二）Tab.3 Comparison between estimated and theoretical modal frequency and damping ratio（Example 2）

由于仅输出的模态参数辨识无法得到模态缩放因子，因此为方便对比，辨识值和理论值都对其最大值进行归一化。模态频率和阻尼比辨识结果如图7所示（线为理论值，点为辨识结果）。

图7 辨识得到的模态振型与理论值对比Fig．7 Comparison between estimated and theoretical mode shapes

由表3和图7可以看出，在实际激励的情况下，即使激励中含有明显的谐波分量，本文提出的辨识方法能够很好的辨识出结构的模态频率、阻尼比和模态振型，且辨识结果不受谐波影响，即辨识结果中不含有谐波造成的虚假模态。

5 结 论

本文推导了多参考点的结构响应传递率表达式，建立了其左矩阵分式多项式的参数化模型，进而给出了辨识问题的最小二乘估计格式，利用Jacobi矩阵的分块性质对最小二乘问题矩阵形式完成了缩减，降低了计算量，通过高维伴随矩阵方法解决了矩阵多项式的特征值求解问题，即多参考多输出的模态参数求解问题，以及通过矩阵伪逆解决了现有方法中载荷工况数与响应点数和参考点数的约束问题。最后，本文通过两个数值算例对所提方法进行了验证，辨识结果表明：本文给出的辨识方法能够很好的辨识出结构的模态频率、阻尼比和模态振型，且能够很好的避免激励中含有的谐波分量对辨识结果造成的影响，解决了传统仅输出模态参数辨识中激励非白对辨识结果造成的影响问题。

［1］李火坤，练继建．高拱坝泄流激励下基于频域法的工作模态参数识别［J］．振动与冲击，2008，27（7）：149－53．

LI Huo-kun，LIAN Ji-jian．Operational modal parameter identification of high arch dam under flood discharge excitation based on frequency domain decomposition method［J］．Journal of Vibration and Shock，2008，27（7）：149－53．

［2］Peeters B，De Roeck G．Reference-based stochastic subspace identification for output-only modal analysis［J］．Mechanical Systems and Signal Processing，1999，13（6）：855－878．

［3］Pickrel C R，White P J．Flight flutter testing of transport aircraft：in-flight modal analysis［C］／／proceedings of the IMAC 21，the International Modal Analysis Conference，Kissimmee（FL），USA，F，2003．

［4］吕中亮，杨昌棋，安培文，等．多点激励模态参数识别方法研究进展［J］．振动与冲击，2011，30（1）：197－203．

LÜ Zhong-liang，YANG Chang-qi，AN Pei-wen，et al．Progress on modal parameter identification with multipleexcitation［J］．Journal of Vibration and Shock，2011，30 （1）：197－203．

［5］Hermans L，Van Der Auweraer H，Hatami A，et al．InflightModal testing and analysis of ahelicopter［C］／／proceedings of the Proceedings of the 17th International Modal Analysis Conference（IMAC17），Kissimmee，FL，USA，F，1999．

［6］Devriendt C，Guillaume P．The use of transmissibility measurements in output-only modal analysis［J］．Mechanical Systems and Signal Processing，2007，21（7）：2689－2696．

［7］Devriendt C，Guillaume P．Identification of modal parameters from transmissibility measurements［J］．Journal of Sound and Vibration，2008，314（12）：343－356．

［8］Devriendt C，De Sitter G，Vanlanduit S，et al．Operational modal analysis in the presence of harmonic excitations by the use of transmissibility measurements［J］．Mechanical Systems and Signal Processing，2009，23（3）：621－35．

［9］Weijtjens W，Sitter G D，Devriendt C，et al．Operational modalanalysisbasedonmultivariabletransmissibility functions：revisited［C］／／proceedings of the IMAC-XXXI，International Modal Analysis Conference，Garden Grove，California，United States，F，2013．

［10］吴川，刘学文，黄醒春．短型浮置板轨道系统隔振性能研究［J］．振动与冲击，2008，27（8）：74－76．

WU Chuan，LIU Xue-wen，HUANG Xing-chun．Study on effectiveness of short floating slab track system［J］．Journal of Vibration and Shock，2008，27（8）：74－76．［11］涂奉臣，张洋，陈照波，等．采用磁流变阻尼器的整星隔振平台试验研究［J］．振动与冲击，2012，31（21）：33 －36．

TU Feng-chen，ZHANG Yang，CHEN Zhao-bo，et al．Tests for a whole-spacecraft vibration isolation platform using a magnetorheological damper［J］．Journal of Vibration and Shock，2012，31（21）：33－36．

［12］Maia N M M，Silva J M M，Ribeiro A M R．The transmissibility concept in multi-degree-of-freedom systems ［J］．Mechanical Systems and Signal Processing，2001，15 （1）：129－137．

［13］Allemang R J，Brown D L．A unified matrix polynomial approach to modal identification［J］．Journal of Sound and Vibration，1998，211（3）：301－322．

Output-only structural modal parameter estimation under non-white excitations based on response transmissibility

ZHOU Si-da，LIU Li，YANG Wu，MA Zhi-sai
（School of Aerospace Engineering，Beijing Institute of Technology，Beijing 100081，China）

In cases of output-only modal parameter estimation for operational engineering structures，such as aerospace structures，the non-white excitations lead to non-physical modes being not able to be identified．For these cases，the recently developed response transmissibility-based methods of modal parameter estimation were studied and an improved poly-reference response transmissibility-based output-only approach of modal parameter estimation was proposed．Firstly，the parametric model of the poly-reference response transmissibility was introduced，it was presented by the left matrix fraction polynomials，the formulations of the least square estimation for unknown parameters were proposed，and the normal equation of the least square estimation was reduced by using the block property of the corresponding Jacobian matrix，it largely decreases the computational costs．Secondly，based on the high-dimensional accompany matrix method，the eigenvalue problem of the matrix polynomials was solved，it led to the solution to the modal parameters．In addition，the original constraints between the number of loading conditions and the number of references were removed by using the pseudo inverse of matrix．Finally，two numerical examples validated the proposed method．The results illustrated that the proposed method can be used to estimate the modal parameters very well and has the advantage of avoiding the influence of harmonics，so the problem due to non-white excitations is solved．

modal parameter estimation；output-only；response transmissibility；non-white noise excitation

N945．14；TH113．1；O324；TN911

A

10．13465／j．cnki．jvs．2014．23．009

北京理工大学基础研究基金（20120142009）

2013－09－12 修改稿收到日期：2013－12－12

周思达男，博士，讲师，1984年7月生

刘莉女，教授，博士生导师，1964年生