CN101464205A - 基于减基法的试验模态分析方法 - Google Patents

Abstract

一种基于减基法的试验模态分析方法，其分析方法步骤为：(1)试验模态分析的准备过程；(2)动态信号采集，并采用LabVIEW图形化软件对信号预处理并得出频率响应输入输出传递函数矩阵；(3)在嵌入LabVIEW的Matlab程序实现减基法，得出包括极点和模态参与因子在内的模态参数；(4)建立稳态图，求解模态振型；(5)模态验证和与其它算法进行结果比较，如果不符合要求将修改减基法的参数重新计算模态参数；如果符合要求则进行下一步；(6)模态动画绘制，系统定型。本发明为缩减自由度，将大规模矩阵投影到小规模矩阵，使运算速度大幅度加快，减少了试验强度和时间，大幅提高试验效率，大大提高了试验的精确性。

Description

基于减基法的试验模态分析方法

技术领域

本发明主要涉及到试验模态分析方法领域，特指一种基于减基法的试验模态分析方法，其主要应用于对汽车、飞机、潜艇、大型建筑桥梁等大型或特大型被测物体的试验模态分析。

背景技术

试验模态分析的过程由试验准备过程、数字信号采集与处理和模态参数识别三个部分组成。实验准备过程包括实验平台的安装、实验品的吊装、实验品几何尺寸的绘制、传感器布点的选取、传感器的安装与调试等工作。此阶段的工作是如何利用有限元理论模态分析的结果(如果没有理论分析结果需要凭借经验或反复测试选择安装点)建立准确的系统模型。数字信号采集与处理的研究目标是如何从传感器得到准确而信噪比高的时域信号并将其转变为频域信号或功率谱信号(时域方法不需要转变)。模态参数识别根据运动方程能确定唯一的数学模型这一原则，研究怎样快速而精确地利用输入和输出数据在模态坐标下对传递函数进行参数辨识，并根据这些参数得到稳态图，从而计算出模态振型并绘制振型动画完成实验模态分析的全过程。

对于实际被测对象，在物理上可以看作无限自由度的系统。现有技术集中在利用自由度缩减方法的实验模态分析，即利用有限元模态分析的低阶振型少数关键位置点建立整体模型，每两点之间的区域用插值方法确定振型。由于每阶振型节点、反节点和中间状态点都不同，此种方法容易使插值方法得到的振型和实际振型差别很大。

对于大型系统而言，要解决插值和实际值的差别，必须在关键点以外更多的布置测点。假设系统时域信号经过离散傅立叶变换后，输入参数矩阵为F(ω)，输出参数矩阵为X(ω)，传递函数矩阵(频率响应函数或功率谱)为H(ω)，则可以得到以下关系：[H(ω)]＝[X(ω)][F(ω)]^-1。这会带来一个新的问题，即输入输出传递函数矩阵的维数会随着被测点的增加而增加，从而传递函数矩阵维数也相应增加。由于需要对矩阵求逆，所以带来计算量增大和求解的困难。怎样在保证一定精度的情况下，求解大型传递函数矩阵成为现阶段研究的热门。现阶段许多学者对待这一难点进行了深入的研究。提出了各种各样的模型降阶方法，一般传统的降阶方法可以分为四类：第一类是动力缩聚法，第二类是基于Krylov子空间技术的方法，第三类是采用Karhunen-Loeve展开(或叫做正常正交分解法)技术的方法，还有一类就是基于Hankel范数近似和平衡截断技术的方法。

但以上方法存在三个缺点：其一是，迭代收敛速度较低，尤其是当降阶模型的特征对逼近精确值时，其收敛速度极低；其二是，迭代格式的收敛性证明非常困难，因此到目前为止还没有文献报道这方面的研究；其三是计算量还是较大，尤其是当主自由度数较大时。

发明内容

本发明要解决的问题就在于：针对现有技术存在的技术问题，本发明提供一种能够实现快速计算、精确度高、具有较好的误差控制、能够减少了试验强度和时间，大幅提高试验效率的基于减基法的试验模态分析方法。

为解决上述技术问题，本发明提出的解决方案为：一种基于减基法的试验模态分析方法，其特征在于步骤为：

①、试验模态分析的准备：安装试验平台并放置试验对象，对试验品的几何尺寸进行绘制，在试验品上布置传感器；

②、动态加速度信号采集：采集试验对象动态下加速度信号，并得出频率响应传递函数，假设采样率为bK/s，相邻两采样点的时间间隔为T_s，经过采样后输入信号和输出信号分别变为离散时间序列a_d(n)和b_d(n)，d＝1...l，其中d为时间采样点计数整数，l为总采样点数；

对a(n)和b(n)进行数字信号处理得到输入输出传递函数矩阵，假设系统输入为n维，输出为m维，则输入函数为A^n×n(k)，输出函数为B^n×m(k)，系统自由度为n×m，传递函数方程为A^n×n(k)X^n×m(k)＝B^n×m(k)(1)，直接求解方程可得传递函数X^n×m(k)＝[A^n×n(k)]^-1B^n×m(k)(2)；

③、通过减基法得到极点p_r和模态参与因子矩阵{L_r}，n×np；

考虑所有采样点l，(1)式可写为A_d ^n×n(k)X_d ^n×m(k)＝B_d ^n×m(k)，d＝1…l，在所有l个采样点中均匀取h个样点，(1)式变为A_d ^n×n(k)X_d ^n×m(k)＝B_d ^n×m(k)，d＝1…h；利用(2)式可得出X_d ^n×m(k)，d＝1…h(3)，并写成参数项和基础项乘积的形式为

设

对于每个采样点而言X^n×m(k)＝Z^m×nα(4)，Z^m×n为基础项；对于每个采样点而言，将(4)式代入(1)式得到：A^n×n(k)Z^n×mα＝B^n×m(k)(5)；

假设输出自由度为1，(5)式变为A^n×n(k)Z^n×1α＝B^n×1(k)(6)，将(6)式两端乘以Z^n×1的转置得到： $Z^{T(p\×n)}A{\left(k\right)}^{n\×n}{Z}^{n\×p}\mathrm{\α}=Z^{T(p\×n)}B{\left(k\right)}^{n\×p}\⇒A\′{\left(k\right)}^{p\×p}\mathrm{\α}=B\′{\left(k\right)}^{p\×p}$ ，p＝1(7)，其中p表示矩阵的维数，第一次计算时为1；对于h个样点中每个采样点d而言，(7)式应写成：A＇(k)_d ^pα_d＝B＇(k)_d ^p，p＝1，d＝1…h(8)；此时的Z看作Z_d，求解(8)式得到d个α值；将每个α代入(4)式可以求得传递函数的变种X_d ^n×m(k)＇，d＝1…h，将每个对应采样点的X_d ^n×m(k)和X_d ^n×m(k)＇相减，得到误差为ε_d，d＝1…h，其中误差最大的计为ε_g；如果ε_g小于预设误差ε，则运算停止，现有X_d ^n×m(k)＇，d＝1…h已满足要求，α就是所需传递函数参数；如果ε_g>ε，则将g点对应的解X_g ^n×m(k)相应的Z_g作为新的列加入到Z_d中形成新的

。返回(7)式，用

替换Z进行计算，此时Z自由度变为n×p，p＝p+1。直到误差ε_f小于预设误差ε，α为最终传递函数参数；将α分成零点参数α_A(B的参数)和极点参数α_B(A的参数)，则(2)式变为 $X^{n\×m}\left(k\right)=\frac{{\mathrm{\α}}_A{Z_A}^{n\×m}\left(k\right)}{{\mathrm{\α}}_B{Z_B}^{n\×m}\left(k\right)}$ (9)，将一般传递函数写成模态参数表达式为：

$\left[X\left(\mathrm{\ω}\right)\right]=\underset{r1}{\overset{N_m}{\mathrm{\Σ}}}(\frac{\left\{{\mathrm{\ψ}}_r\right\}{\left\{L_r\right\}}^T}{\mathrm{j\ω}-p_r}+\frac{{\left\{{\mathrm{\ψ}}_r\right\}}^{*}{\left\{L_r\right\}}^H}{\mathrm{j\ω}-{p_r}^{*}})-\frac{\left[\mathrm{LR}\right]}{{\mathrm{\ω}}^2}+\left[\mathrm{UR}\right]---\left(10\right)$

其中{ψ_r}是模态振型，{L_r}是模态参与因子，p_r是极点；

$\left[\begin{array}{ccccc}\left[o\right]& \left[I\right]& \·& \left[o\right]& \left[o\right]\\ \left[o\right]& \left[o\right]& \·& \left[o\right]& \left[o\right]\\ \·& \·& \·& \·& \·\\ \left[o\right]& \left[o\right]& \·& \left[o\right]& \left[I\right]\\ -{\left[{\mathrm{\α}}_0\right]}^T& -{\left[{\mathrm{\α}}_1\right]}^T& \·& -{\left[{\mathrm{\α}}_{p-2}\right]}^T& -{\left[{\mathrm{\α}}_{p-1}\right]}^T\end{array}\right]\left[V\right]=\left[V\right]\left[\mathrm{\Λ}\right]---\left(11\right)$

在(10)式中，方程左边为“友”矩阵，阶数为np×np，特征值矩阵[Λ]的对角线为它的特征值 ${\mathrm{\λ}}_r=e^{p_r\mathrm{\Δt}}(r=1,\·\·\·\mathrm{np})$ ，即求出极点p_r；而特征向量矩阵[V]最后m行即模态参与因子矩阵{L_r}，n×np；

④、建立稳态图，求解模态振型：根据(9)式绘制求和的FRF图，在取不同的模型阶次p时，分别计算出相应的极点和模态参与因子。由于极点p_r和模态参与因子矩阵{L_r}，n×np已求出，利用频率参数直接识别法可由不同的频率对应的k列出上式方程，用线性最小二乘法求出未知模态振型{ψ_r}及上下残余项[UR]和[LR]；

⑤、模态验证和结果对比与分析：主要完成针对结果的评判工作。利用模态判定准则 $\mathrm{MAC}=\frac{{\left|{\left\{{\mathrm{\ψ}}^{*}\right\}}_r^T{\left\{\mathrm{\ψ}\right\}}_s\right|}^2}{\left({\left\{{\mathrm{\ψ}}^{*}\right\}}_r^T{\left\{\mathrm{\ψ}\right\}}_r\right)\left({\left\{{\mathrm{\ψ}}^{*}\right\}}_s^T{\left\{\mathrm{\ψ}\right\}}_s\right)}$ (12)对是否真实模态进行判断。若{ψ}_r和{ψ}_s本质上是同一模态则有MAC≈1，若{ψ}_r和{ψ}_s本质上是不同模态则有MAC≈0。用此方法既用来验证此方法解出的不同阶次模态振型是否满足MAC≈0，也用来比较基于减基法解出的模态振型和其他方法的模态振型的关系，从而判断此方法的有效性。再根据统计学的原理将此方法得到的模态振型和其他方法进行一致性比较，从而验证此方法的有效性。经过模态判定准则的判断和一致性比较，如果符合误差要求则结束整个运算过程，确定减基法的参数h和p值；如果不符合误差要求则调整参数h和p，从步骤③重新开始计算直至符合要求为止。这样确定了减基法的参数h和p值，整体基于减基法的试验模态分析核心计算过程结束。

⑥、模态动画绘制：得出各点每个方向的模态振型，与布点图对号入座，就得到模态振型动画，从而完成整个试验模态分析全过程。

所述步骤②中，采用激振器激振法或力锤法对试验对象进行试验。

所述步骤③中，采用自由度缩减的方法将大型矩阵投影到小型矩阵进行运算。

所述步骤③中，采用误差控制的方法保证计算结果的精度。

所述步骤⑤中，利用模态判定准则进行误差分析和比较。

与现有技术相比，本发明的优点就在于：

1、本发明是一种能够实现快速计算、精确度高、具有较好的误差控制、能够减少了试验强度和时间，大幅提高试验效率的基于减基法的试验模态分析方法；

2、现在实验模态分析正向大系统高精度方向发展。特别对于汽车、飞机、潜艇、大型建筑桥梁等大型或特大型被测物体，不但要考虑缩减自由度下的低阶模态，更要考虑大自由度下的高阶模态。通过本发明的方法，可以看到外界复杂高频激励下被测物体的实际相应状态，使其能相应的改进结构，对降低噪音(对于汽车、飞机、潜艇非常重要)、预判损伤(对于大型建筑桥梁)起决定性的作用。经过实践检验，本发明的方法在这些方面有很好的作用，比现在通用的建立在缩减自由度下的方法得到的模态振型更加精确可靠，同时运算速度也更高。

附图说明

图1是本发明的流程示意图；

图2是本发明中减基法的流程示意图；

图3是具体实施例中传感器布置的示意图；

图4是具体实施例中根据本发明最终得到模态振型动画的示意图。

具体实施方式

以下将结合附图和具体实施例对本发明做进一步详细说明。

如图1和图2所示，本发明公开了一种基于减基法的试验模态分析方法，其步骤为：

(1)、试验模态分析的准备过程，包括检查被测对象是否满足试验的必要条件即线性时不变性、可观测性和Maxwell互易性原理，试验平台的安装，试验品的吊装，试验品几何尺寸的绘制，试验方法的确定包括使用激振器激振法还是力锤法(我们这里采用力锤法)，传感器布点的选取，传感器的安装和初始化等工作。

(2)、动态加速度信号采集，并采用LabVIEW图形化软件对信号预处理并得出频率响应传递函数。

假设输入信号为a(t)，输出信号为b(t)。动态信号采集的关键在于：一、硬件平台采用基于PXI标准的业界领先的美国国家仪器(NI)公司的数据采集与分析平台，使不同板卡和不同机箱之间达到10MHz的同步，能同时对1000通道的动态信号进行采集并保证其精确性和同步性，这使得复杂结构的试验模态分析得以实现。二、动态信号采集卡采用NI公司8通道、独立A/D、24位、102.4K/s采样率的PXI-4472，它的信号分辨率高达110DB；采样方式为Delta-sigma方式，能很好的解决混叠失真问题。这两点优势使动态信号能高精度无失真的用于下一步。三、根据不同的激振方法选择不同的窗函数。激振器法选用矩形窗函数，而锤击法采用指数窗函数。

经过以上三步，假设采样率为bK/s，相邻两采样点的时间间隔为T_s。经过采样后输入输出信号变为离散时间序列a_d(n)和b_d(n)，d＝1...l，其中d为时间采样点计数整数，l为总采样点数。

(3)、对a(n)和b(n)进行数字信号处理得到输入输出传递函数矩阵，关键在于以下三步：

一、业界流行的NI公司的LabVIEW图形化软件和PXI4472能实现无缝连接，使用LabVIEW软件编写的程序能实时调用采集得到的动态离散时域信号。

二、为了输入输出的时域信号转变为频域信号，并使计算机高效的信号处理成为可能，需要进行以下四步：第一步，将离散化的时域信号进行离散时间傅里叶变换(DTFT)，频谱被周期化，于是得到： $A\left(e^{\mathrm{j\ω}}\right)=\underset{n=-\∞}{\overset{+\∞}{\mathrm{\Σ}}}a\left(n\right)e^{-\mathrm{j\ωn}}$ 和 $B\left(e^{\mathrm{j\ω}}\right)=\underset{n=-\∞}{\overset{+\∞}{\mathrm{\Σ}}}b\left(n\right)e^{-\mathrm{j\ωn}},$ 其中ω为圆频率；第二步，再将频域离散化，得到离散周期傅里叶级数(DFS) $A\left(k\right)=\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}\tilde{a}\left(n\right)e^{-j\frac{2\mathrm{\π}}{N}\mathrm{kn}}$ 和 $B\left(k\right)=\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}\tilde{b}\left(n\right)e^{-j\frac{2\mathrm{\π}}{N}\mathrm{kn}},$ 其中N为时间采样点总数或谱线总数，k表示某次谐波。这样时域信号进一步被周期化。第三步，考虑到周期离散化的时域和频域，只取一个周期研究，也就是众所周知的离散傅里叶变换(DFT)，分别得到 $A\left(k\right)=\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}a\left(n\right)W_N^{\mathrm{kN}}$ 和 $B\left(k\right)=\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}b\left(n\right)W_N^{\mathrm{kN}}.$ 第四步，选择快速傅立叶变换(FFT)作为DFT的一种快速算法，A(k)和B(k)结果不变但计算速度大幅度提高。

三、假设系统输入为n维，输出为m维，所以输入函数为A^n×n(k)，输出函数为B^n×m(k)。系统自由度为n×m，传递函数方程为A^n×n(k)X^n×m(k)＝B^n×m(k)(1)，直接求解方程可得传递函数X^n×m(k)＝[A^n×n(k)]^-1B^n×m(k)(2)。

(4)、减基法的实现：

使用减基法进行模态参数拟合是本专利的核心，其采用Matlab进行编程并嵌入LabVIEW中来实现以下算法：

考虑所有采样点l，(1)式可写为A_d ^n×n(k)X_d ^n×m(k)＝B_d ^n×m(k)，d＝1…l，在所有l个采样点中均匀取h个样点，(1)式变为A_d ^n×n(k)X_d ^n×m(k)＝B_d ^n×m(k)，d＝1…h。利用(2)式可得出X_d ^n×m(k)，d＝1…h(3)，并写成参数项和基础项乘积的形式为

设对于每个采样点而言X^n×m(k)＝Z^m×nα(4)，Z^m×n为基础项。

为了起到“减基”的作用，需要将大维空间方程投影到减缩空间中，必须将新的个体样本对应的解向量作为基向量Z＇进行迭代缩减运算。

对于每个采样点而言，将(4)式代入(1)式得到：A^n×n(k)Z^n×mα＝B^n×m(k)(5)。

根据在步骤(1)中的Maxwell互易性原理可知对多输入多输出系统参数矩阵的一行或一列的就可以得到模态参数。为了简化运算，实现多次循环求解模态参数，假设输出自由度为1，(5)式变为A^n×n(k)Z^n×1α＝B^n×1(k)(6)

(6)式两端乘以Z^n×1的转置得到：

$Z^{T(p\×n)}A{\left(k\right)}^{n\×n}{Z}^{n\×p}\mathrm{\α}=Z^{T(p\×n)}B{\left(k\right)}^{n\×p}\⇒A\′{\left(k\right)}^{p\×p}\mathrm{\α}=B\′{\left(k\right)}^{p\×p},p=1---\left(7\right)$

其中p表示矩阵的维数，第一次计算时为1。对于h个样点中每个采样点d而言，(7)式应写成：A＇(k)_d ^pα_d＝B＇(k)_d ^p，p＝1，d＝1…h(8)。此时的Z看作Z_d。

求解(8)式得到d个α值。将每个α代入(4)式可以求得传递函数的变种X_d ^n×m(k)＇，d＝1…h，将每个对应采样点的X_d ^n×m(k)和X_d ^n×m(k)＇相减，得到误差为ε_d，d＝1…h。其中误差最大的计为ε_g。如果ε_g小于预设误差ε，则运算停止，现有X_d ^n×m(k)＇，d＝1…h已满足要求，α就是所需传递函数参数。

如果ε_g>ε，则将g点对应的解X_g ^n×m(k)相应的Z_g作为新的列加入到Z_d中形成新的

。返回(7)式，用

替换Z进行计算，此时Z自由度变为n×p，p＝p+1。直到误差ε_f小于预设误差ε，α为最终传递函数参数。

将α分成零点参数α_A(B的参数)和极点参数α_B(A的参数)，则(2)式变为 $X^{n\×m}\left(k\right)=\frac{{\mathrm{\α}}_A{Z_A}^{n\×m}\left(k\right)}{{\mathrm{\α}}_B{Z_B}^{n\×m}\left(k\right)}\left(9\right).$ 将一般传递函数写成模态参数表达式为：

$\left[X\left(\mathrm{\ω}\right)\right]=\underset{r1}{\overset{N_m}{\mathrm{\Σ}}}(\frac{\left\{{\mathrm{\ψ}}_r\right\}{\left\{L_r\right\}}^T}{\mathrm{j\ω}-p_r}+\frac{{\left\{{\mathrm{\ψ}}_r\right\}}^{*}{\left\{L_r\right\}}^H}{\mathrm{j\ω}-{p_r}^{*}})-\frac{\left[\mathrm{LR}\right]}{{\mathrm{\ω}}^2}+\left[\mathrm{UR}\right]---\left(10\right)$

其中{ψ_r}是模态振型，{L_r}是模态参与因子，p_r是极点。我们利用α_A扩展的“友”矩阵(Companion Matrix)的特征值分解得出极点和模态参与因子。

$\left[\begin{array}{ccccc}\left[o\right]& \left[I\right]& \·& \left[o\right]& \left[o\right]\\ \left[o\right]& \left[o\right]& \·& \left[o\right]& \left[o\right]\\ \·& \·& \·& \·& \·\\ \left[o\right]& \left[o\right]& \·& \left[o\right]& \left[I\right]\\ -{\left[{\mathrm{\α}}_0\right]}^T& -{\left[{\mathrm{\α}}_1\right]}^T& \·& -{\left[{\mathrm{\α}}_{p-2}\right]}^T& -{\left[{\mathrm{\α}}_{p-1}\right]}^T\end{array}\right]\left[V\right]=\left[V\right]\left[\mathrm{\Λ}\right]---\left(11\right)$

在(10)式中，方程左边为“友”矩阵，阶数为np×np，特征值矩阵[A]的对角线为它的特征值 ${\mathrm{\λ}}_r=e^{p_r\mathrm{\Δt}}(r=1,\·\·\·\mathrm{np})$ ，即求出极点p_r。而特征向量矩阵[V]最后m行即模态参与因子矩阵{L_r}，n×np。

(5)、稳态图的建立和模态阶次的确定：

根据(9)式绘制求和的FRF图。在取不同的模型阶次p时，分别计算出相应的极点和模态参与因子。如果模态频率、阻尼比和模态参与因子在规定的容差范围内不随p的取值不同而变化，就在图上注明符号“S”，并认为此时模态频率为系统某阶频率，从而确定在分析频带分为内物理模态的阶次N。

(6)、LSFD法求解模态振型

由于极点p_r和模态参与因子矩阵{L_r}，n×np已求出，利用LSFD(频率参数直接识别法)可由不同的频率对应的k列出上式方程，用线性最小二乘法求出未知模态振型{ψ_r}及上下残余项[UR]和[LR]。

(7)、模态验证和与其他主流方法的比较

根据模态判定准则(MAC，见(11)式)对采用不同方法计算的同一物理振型进行一致性的比较。我们将减基法和工业标准LSCE方法对某白车身进行MAC值的比较。其中{ψ_r}为减基法求得的模态振型，

为它的共轭转置矩阵；{ψ}_s为LSCE法求得的模态振型，

为它的共轭转置矩阵。比较结果见表1。

$\mathrm{MAC}=\frac{{\left|{\left\{{\mathrm{\ψ}}^{*}\right\}}_r^T{\left\{\mathrm{\ψ}\right\}}_s\right|}^2}{\left({\left\{{\mathrm{\ψ}}^{*}\right\}}_r^T{\left\{\mathrm{\ψ}\right\}}_r\right)\left({\left\{{\mathrm{\ψ}}^{*}\right\}}_s^T{\left\{\mathrm{\ψ}\right\}}_s\right)}---\left(12\right)$

表1 减基法与LSCE方法的MAC值比较

从表1的对比结果中可看出，前5阶频率值非常接近，对应模态的对角线对称MAC值也相近，说明两种方法结果有一致性。并且对角线上的MAC值非常接近1，非对角线上的MAC值全部小于2％并接近0说明对同一物理振型的估计接近1，对不同物理振型的估计接近0，两种方法均符合模态分析实际情况。如果对于同一振型的MAC比较值大于2％或不满足预定的要求则说明减基法的参数需要调整。一般采用调整采样点d和预设误差标准ε的方法加以解决。适当增加采样点数一定程度增加计算精度但延长了计算时间，而采用调整取点方式能达到很好的效果，采用类似正态积分函数取点法代替平均取点法，即在某频率段模态密集区密集取点而其它区域稀疏取点。这样同样的采样点数效率大幅度提高。调整预设误差标准ε则要根据具体要求而定。一般情况下，MAC比较值大于预设则降低ε直到比较值小于预设；MAC比较值小于预设则ε还有提高的空间。

以某大货车车头为被测对象，在精度相近并同时符合实际情况的条件下，经过测算减基法所需时间在100×1自由度情况下和LSCE方法接近，但在10000×1自由度情况下减少15％，充分体现了其在大自由度时计算效率的优势。

从表1还可以看出基于减基法解出的模态振型和其他方法的模态振型的关系。对角线上的值接近1而非对角线的值接近0，说明对角线对应的是同一模态而非对角线对应的是不同模态，从而判断减基法的有效性。再根据统计学的原理将此方法得到的模态振型和其他方法进行一致性比较，从而验证此方法的有效性。经过模态判定准则的判断和一致性比较，如果符合误差要求则结束整个运算过程，确定减基法的参数h和p值；如果不符合误差要求则调整参数h和p，从步骤③重新开始计算直至符合要求为止。这样确定了减基法的参数h和p值，整体基于减基法的试验模态分析核心计算过程结束。

(8)、模态动画绘制：

对于模态动画绘制，我们以某大货车车头为被测对象进行说明。第(1)步试验模态分析的准备过程中有传感器布点的步骤，实际布点图如图3所示。第(7)步得出各点每个方向的模态振型，与布点图对号入座，就得到模态振型动画(图4为动画截图，移去某些非关键点)，从而完成整个试验模态分析全过程。

以上所述仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (5)

1、一种基于减基法的试验模态分析方法，其特征在于步骤为：

①、试验模态分析的准备：安装试验平台并放置试验对象，对试验品的几何尺寸进行绘制，在试验品上布置传感器；

②、动态信号采集：采集试验对象动态下加速度信号，并得出频率响应传递函数，假设采样率为bK/s，相邻两采样点的时间间隔为T_s，经过采样后输入信号和输出信号分别变为离散时间序列a_d(n)和b_d(n)，d＝1...l，其中d为时间采样点计数整数，l为总采样点数；

对a(n)和b(n)进行数字信号处理得到输入输出传递函数矩阵，假设系统输入为n维，输出为m维，则输入函数为A^n×n(k)，输出函数为B^n×m(k)，系统自由度为n×m，传递函数方程为A^n×n(k)X^n×m(k)＝B^n×m(k)(1)，直接求解方程可得传递函数X^n×m(k)＝[A^n×n(k)]^-1B^n×m(k)(2)；

③、通过减基法得到极点p_r和模态参与因子矩阵{L_r}，n×np；

考虑所有采样点l，(1)式可写为 ${A_d}^{n\×n}\left(k\right){X_d}^{n\×m}\left(k\right)={B_d}^{n\×m}\left(k\right),$ d＝1…l，在所有l个采样点中均匀取h个样点，(1)式变为 ${A_d}^{n\×n}\left(k\right){X_d}^{n\×m}\left(k\right)={B_d}^{n\×m}\left(k\right),d=1\·\·\·h;$ 利用(2)式可得出

并写成参数项和基础项乘积的形式为

设

对于每个采样点而言X^n×m(k)＝Z^m×nα(4)，Z^m×n为基础项；对于每个采样点而言，将(4)式代入(1)式得到：A^n×n(k)Z^n×mα＝B^n×m(k)(5)；

假设输出自由度为1，(5)式变为A^n×n(k)Z^n×1α＝B^n×1(k)(6)，将(6)式两端乘以Z^n×1的转置得到： $Z^{T(p\×n)}A{\left(k\right)}^{n\×n}{Z}^{n\×p}\mathrm{\α}=Z^{T(p\×n)}B{\left(k\right)}^{n\×p}\⇒A\′{\left(k\right)}^{p\×p}\mathrm{\α}=B\′{\left(k\right)}^{p\×p},p=1---\left(7\right),$ 其中p表示矩阵的维数，第一次计算时为1；对于所有h个样点中每个采样点d而言，(7)式应写成： $A\′{{\left(k\right)}_d}^p{\mathrm{\α}}_d=B\′{{\left(k\right)}_d}^p,$ p＝1，d＝1…h(8)；此时的Z看作Z_d，求解(8)式得到d个α值；将每个α代入(4)式可以求得传递函数的变种

将每个对应采样点的

和

相减，得到误差为ε_d，d＝1…h，其中误差最大的计为ε_g；如果ε_g小于预设误差ε，则运算停止，现有

已满足要求，α就是所需传递函数参数；如果ε_g>ε，则将g点对应的解相应的Z_g作为新的列加入到Z_d中形成新的

返回(7)式，用

替换Z进行计算，此时Z自由度变为n×p，p＝p+1；直到误差ε_f小于预设误差ε，α为最终传递函数参数；将α分成零点参数α_A和极点参数α_B，则(2)式变为 $X^{n\×m}\left(k\right)=\frac{{\mathrm{\α}}_A{Z_A}^{n\×m}\left(k\right)}{{\mathrm{\α}}_B{Z_B}^{n\×m}\left(k\right)}---\left(9\right),$ 将一般传递函数写成模态参数表达式为：

$\left[X\left(\mathrm{\ω}\right)\right]=\underset{r1}{\overset{N_m}{\mathrm{\Σ}}}(\frac{\left\{{\mathrm{\ψ}}_r\right\}{\left\{L_r\right\}}^T}{\mathrm{j\ω}-p_r}+\frac{{\left\{{\mathrm{\ψ}}_r\right\}}^{*}{\left\{L_r\right\}}^H}{\mathrm{j\ω}-{p_r}^{*}})-\frac{\left[\mathrm{LR}\right]}{{\mathrm{\ω}}^2}+\left[\mathrm{UR}\right]---\left(10\right)$

其中{ψ_r}是模态振型，{L_r}是模态参与因子，p_r是极点；

$\left[\begin{array}{ccccc}\left[o\right]& \left[I\right]& \·& \left[o\right]& \left[o\right]\\ \left[o\right]& \left[o\right]& \·& \left[o\right]& \left[o\right]\\ \·& \·& \·& \·& \·\\ \left[o\right]& \left[o\right]& \·& \left[o\right]& \left[I\right]\\ -{\left[{\mathrm{\α}}_0\right]}^T& -{\left[{\mathrm{\α}}_1\right]}^T& \·& -{\left[{\mathrm{\α}}_{p-2}\right]}^T& -{\left[{\mathrm{\α}}_{p-1}\right]}^T\end{array}\right]\left[V\right]=\left[V\right]\left[\mathrm{\Λ}\right]---\left(11\right)$

在(10)式中，方程左边为“友”矩阵，阶数为np×np，特征值矩阵[Λ]的对角线为它的特征值 ${\mathrm{\λ}}_r=e^{p_r\mathrm{\Δt}}(r=1,\·\·\·\mathrm{np})$ ，即求出极点p_r；而特征向量矩阵[V]最后m行即模态参与因子矩阵{L_r}，n×np；

④、建立稳态图，求解模态振型：根据(9)式绘制求和的FRF图，在取不同的模型阶次p时，分别计算出相应的极点和模态参与因子；由于极点p_r和模态参与因子矩阵{L_r}，n×np已求出，利用频率参数直接识别法可由不同的频率对应的k列出上式方程，用线性最小二乘法求出未知模态振型{ψ_r}及上下残余项[UR]和[LR]；

⑤、模态验证和结果对比与分析：主要完成针对结果的评判工作，利用模态判定准则 $\mathrm{MAC}=\frac{{\left|{\left\{{\mathrm{\ψ}}^{*}\right\}}_r^T{\left\{\mathrm{\ψ}\right\}}_s\right|}^2}{\left({\left\{{\mathrm{\ψ}}^{*}\right\}}_r^T{\left\{\mathrm{\ψ}\right\}}_r\right)\left({\left\{{\mathrm{\ψ}}^{*}\right\}}_s^T{\left\{\mathrm{\ψ}\right\}}_s\right)}---\left(12\right)$ 对是否真实模态进行判断；若｛ψ}_r和｛ψ}_s本质上是同一模态则有MAC≈1，若{ψ}_r和｛ψ}_s本质上是不同模态则有MAC≈0；用此方法既用来验证此方法解出的不同阶次模态振型是否满足MAC≈0，也用来比较基于减基法解出的模态振型和其他方法的模态振型的关系，从而判断此方法的有效性；再根据统计学的原理将此方法得到的模态振型和其他方法进行一致性比较，从而验证此方法的有效性；经过模态判定准则的判断和一致性比较，如果符合误差要求则结束整个运算过程，确定减基法的参数h和p值；如果不符合误差要求则调整参数h和p，从步骤③重新开始计算直至符合要求为止；这样确定了减基法的参数h和p值，整体基于减基法的试验模态分析核心计算过程结束；

⑥、模态动画绘制：得出各点每个方向的模态振型，与布点图对号入座，就得到模态振型动画，从而完成整个试验模态分析全过程。

2、根据权利要求1所述的基于减基法的试验模态分析方法，其特征在于：所述步骤②中，采用激振器激振法或力锤法对试验对象进行试验。

3、根据权利要求2所述的基于减基法的试验模态分析方法，其特征在于：所述步骤③中，采用自由度缩减的方法将大型矩阵投影到小型矩阵进行运算。

4、根据权利要求2所述的基于减基法的试验模态分析方法，其特征在于：所述步骤③中，采用误差控制的方法保证计算结果的精度。

5、根据权利要求2所述的基于减基法的试验模态分析方法，其特征在于：所述步骤⑤中，利用模态判定准则进行误差分析和比较。

Images