一种基于奇异值分解技术的模型定阶方法

王树青， 林裕裕，孟元栋，高志强

(1.中国海洋大学 山东省海洋工程重点实验室，青岛 266100;2.海洋石油工程股份有限公司，天津 300451)

目前，基于动力特性变化的结构健康检测正受到国内外学者的重视，而模态参数识别非常关键。如果识别的模态参数存在较大的误差或存在虚假模态，这势必会影响到结构健康检测结果。模态参数识别方法一般都是基于一定的系统模型的，如ARMA是基于时间序列模型，SSI和ERA是基于状态空间模型。因此要准确识别系统的模态参数，就需要确定系统的阶次(或参与的模态数)。当前，模型阶次的确定已经成了模态参数识别中非常关键的环节。此外由于测试数据不可避免地包含大量噪声，这给模态参数识别中模型阶次的确定造成了很大的困难。合理地确定模型阶次已成为模态参数识别领域的最为重要的问题。

目前，对模态参数识别中模型阶次的确定问题，已经出现了一些研究方法，如稳定图法［1－2］。该方法通过在频谱图上标示出满足一定条件的稳定极点，并被认为是系统的真实极点。但这种方法不能完全排除噪声模态，特别是随着模型阶数的升高，一些拟合模型的拟合模态往往容易趋于稳定，用稳定图很难完全并正确确定模型阶次［3］。近来，基于奇异值分解的模型定阶与降噪技术［4－6］又引起了国内外有关专家的重视。该类方法一般将奇异值由大到小按降序排列，并将奇异值以最大值归一化，通过画出奇异值归一化曲线，在曲线上找到突降的位置，该处对应的奇异值个数即为模型阶次，也即信号中包含模态数目的两倍。对受噪声影响的数据，奇异值曲线突降不明显，而是趋向于一条水平渐近线。一般认为奇异值曲线开始变为水平的点对应模型的阶次［4－6］。在如何有效确定奇异值曲线突变点上，又发展了奇异谱技术、奇异熵定阶技术等［7－9］。易伟建［10］提出了一种动力系统模型定阶方法一残差期望法。该方法利用系统的脉冲响应函数建立Hankel矩阵，对其进行奇异值分解后，构造残差期望比。在模型阶次处，残差期望会由偏离零值较大变化到接近于零，因此会在此处出现较大的残差期望比。因此在残差期望变化处或残差期望比的最后一个尖峰处则为可判定模型的真实阶次。赵学智［11－12］提出了利用奇异值差分谱、奇异值曲率谱的概念来描述含噪信号的奇异值曲线的转折点情况，进而利用奇异值曲线的最大转折点来确定有效奇异值的个数。

本文提出利用奇异值的相对变化率来确定模型的阶次。该方法利用结构的量测脉冲响应信号构造Hankel矩阵，对其进行奇异值分解后计算奇异值的相对变化率，变化率最大的地方对应着模型的阶次。通过数值算例研究了量测噪声对模型阶次确定的影响，并利用某悬臂梁实验数据验证了本方法的有效性。

1 奇异值分解技术与模型阶次确定

某离散时间系统的状态空间方程为:

其中:A、B、C分别为系统的特征矩阵、输入矩阵和输出矩阵。当输入为脉冲激励时，输出脉冲响应数据Yk可表示为:

利用脉冲响应序列构建Hankel矩阵，得到:

其中:m、n分别为Hankel矩阵的行数和列数。对H(0)进行奇异值分解，得到:

式(4)中，上标“T”表示矩阵转置;矩阵U和V为正交矩阵;∑为对角矩阵，其对角元素为降序排列的奇异值。理论上，超出矩阵秩的奇异值应当为零，即矩阵H的秩恰好等于模型的阶数r，此时:

式(5)中，R=min{m，n}。可以看出，奇异值在模型的阶次r处会产生突降。然而当量测信号受到噪声污染时，非零奇异值的个数远大于模型的阶数r，但此时大于模型阶次的奇异值都是噪声的贡献，其值一般都很小，如果选择一个合适的临界值ε，使得公式(6)成立，即可以确定模型的阶次。

由于噪声的影响，奇异值在模型的阶次处突降变得不是很明显，一般通过观察奇异值曲线，将开始趋向于一条水平渐近线的地方定义为模型的阶次。这种通过肉眼观察的方法在噪声较小的情况下还是适用的;但当噪声变大时，通过肉眼观察也很难确定模型的阶次。能否提出一个量化的模型阶次指标来确定模型的阶次呢?本文提出利用奇异值的相对变化率来作为模型阶次指标，即定义模型阶次指标为:

由于奇异值是降序排列的，在突降点大的地方，模型阶次指标MOCi也是一个大值，即定义模型阶次指标MOCi最大值对应的阶次为模型的阶次。其优点为利用一个量化的指标来确定模型阶次，避免了肉眼观察带来的失误。

2 五自由度质量－弹簧－阻尼系统

2.1 数据的模拟

采用文献［6］中的数值算例，建立一个五自由度的质量－弹簧－阻尼系统的数值模型如图1所示。单元的质量、刚度和阻尼系数分别为mn=50kg、kn=2.9×107N/m、cn=1000 N·s/m。通过特征值分析，可以得到模态频率的理论值为:34.499Hz、100.700Hz、158.730Hz、203.880Hz、232.520Hz;5 阶模态阻尼比的理论值为:0.0037374、0.010909、0.017197、0.022092、0.025198。

图1 五自由度质量－弹簧－阻尼系统Fig.1 Mass-spring-damper system

在第一自由度输入单位脉冲激励，从而可以计算系统的脉冲响应函数。采用Matlab编制程序，得到系统的脉冲响应函数和频率响应函数，如图2所示，其中采样频率为500Hz。频率响应函数共有5个峰值，分别对应着5阶模态频率;另外也可以看出，第5阶模态(频率为232.520Hz)不明显，即该阶模态对脉冲响应的贡献非常小。

2.2 噪声对模型定阶的影响

在没有噪声影响的情况下，取128个点(27个点)的第1自由度脉冲响应序列构建 Hankel矩阵(H64×64)，对其进行奇异值分解，得到归一化的奇异值曲线如图3(a)所示。从中可以看出奇异值曲线在第10个奇异值处出现陡降，其后面的奇异值接近于零，即判定该模型的阶次为10阶(模态自由度的2倍)。利用归一化的奇异值计算奇异值的相对变化率——模型阶次指标，如图3(b)所示。很明显，模型阶次指标的最大值正好对应10阶。用模型阶次指标可以很好地确定该模型的正确阶次。

实测信号不可避免地会受到各种噪声的影响。为了研究不同程度的噪声对模型阶次的影响，利用Matlab模拟了高斯白噪声，叠加到脉冲响应序列中。噪声水平通过一个百分比来定量描述，该百分比定义为白噪声的标准差和精确信号的标准差之比。

分别模拟1%、2.5%、5%、10%、20%的噪声，叠加到脉冲响应序列中，然后进行模型阶次计算，计算结果如图4所示。图4中，左侧为不同噪声水平对应的归一化奇异值曲线，右侧为对应不同噪声的模型阶次指标。可以看出，随着量测噪声的增强，模型阶次指标的最大值在减小，也就是说奇异值曲线中的突降点变化越来越小，如左图所示。当噪声较小时(如小于10%)，利用模型阶次指标可以很好地确定模型的阶次为10。当量测噪声变大时，可以判定模型阶次由10变成了8，即丢失了1阶模态。这点可以从图2中的频率响应函数得到解释。如前所述，第5阶模态的峰值不明显，对脉冲响应的贡献较小，更容易受到噪声的干扰。在噪声较严重时，第5阶模态已经完全被噪声淹没了，所以模型的阶次由10降到了8。同时也可以看出，随着噪声程度的增加，按照传统的方法(比如按照奇异值趋于一条水平渐近线的地方)来定义模型阶次，越来越困难;因为随着噪声程度的增加，奇异值趋于一条水平渐近线地方并不明显，如图4所示。

图4 不同噪声情况下的奇异值曲线与模型阶次指标(数据长度:27个点)Fig.4 Normalized singular values and model order indicators with different noise level(data length:27points)

2.3 数据长度对模型阶次的影响

2.2 节利用128个点(27个点)的脉冲响应数据进行模型阶次确定。为了研究数据长度对模型阶次确定的影响，分别取64(26个点)、256(28个点)进行计算，研究归一化奇异值曲线与模型阶次指标的变化(篇幅所限，图形未列出)。研究发现随着脉冲响应序列数据长度的增加，模型阶次越容易受到噪声的影响。

3 海洋平台模型

为了检验本文提出的模型定阶方法的有效性，下面考虑一个更加复杂的三维空间结构数值算例。该空间结构为某海洋平台模型，如图 5所示，图中的数字为模型的节点编号。考虑到实际海洋平台的阻尼比一般都不大，取前两阶阻尼比为0.01来确定比例阻尼系数α和β。平台模型的前6阶频率如表1所示。

图5 海洋平台模型Fig.5 The finite element model of an offshore platform

表1 海洋平台模型的前6阶频率(Hz)Tab.1 The first six natural frequencies of the offshore platform model(Hz)

在节点29同时施加x和y向脉冲激励，计算各个节点的脉冲响应，其中采样频率为20Hz，采样点数为1024点。然后叠加不同程度的噪声，研究模型定阶的有效性。

首先研究没有噪声干扰情况下的模型定阶。对三维空间结构，水平振动响应分为x向和y向。分别取128个点(27个点)的x向和y向脉冲响应数据进行模型阶次确定，结果如图6和图7所示，其中(a)为脉冲响应，(b)频率响应函数，(c)为模型阶次。从图中可以看出，利用x向和y向脉冲响应数据确定的模型阶次都是8，即都包含了4阶模态。分析图6(b)可以看出，这4阶模态分别为x向1阶弯曲模态、1阶扭转模态、x向2阶模态与2阶扭转模态。而图7(b)中的4阶模态分别为y向1阶弯曲模态、1阶扭转模态、y向2阶模态与2阶扭转模态。

图6 x向脉冲响应模型定阶Fig.6 Model order using x-dir impulse

图7 y向脉冲响应模型定阶Fig.7 Model order using y-dir impulse

如第2节，分别模拟1%、2.5%、5%、10%、20%的噪声，叠加到x向脉冲响应序列中，然后进行模型阶次计算，计算结果如图8所示。可以看出，随着量测噪声的增强，模型阶次指标的最大值在减小，但始终可以确定模型阶次为8。在不同噪声水平下，利用y向脉冲响应进行模型阶次确定结果同图8类似，由于版面限制没有列出。同2.3节类似，分别采取数据长度为26、27、28等不同时长的脉冲响应进行模态阶次确定，研究发现，对本节的三维空间结构，数据长度对模型阶次确定的影响非常小。

图8 不同噪声影响x向脉冲响应模型定阶Fig.8 Model order determination using x-dir impulse under different noise level

4 实验数据检验

为了验证本文提出方法的有效性，采用某悬臂梁进行冲击实验，布置传感器量测其振动响应，采样频率为200Hz。图9(a)为量测的脉冲响应，图9(b)为其频响函数，可以看出在100Hz以内共有3阶模态，其频率分别为 5.52Hz、34.71Hz、97.2Hz。图10 为其归一化奇异值曲线和模型阶次指标，可以看出，利用本文提出的模型阶次指标可以很容易确定的模型的阶次为6(即包含的3阶模态的2倍)。同时从图10中也可以发现，如果利用奇异值趋于水平渐近线不容易确定模型的正确阶次。

5 结论

(1)利用结构的量测脉冲响应信号构造Hankel矩阵，利用奇异值分解技术得到其奇异值，提出利用奇异值的相对变化率来确定模型的阶次。数值算例和模型实验数据验证了本方法的有效性。

(2)利用一个简单的5自由度体系和一个三维的空间结构数值算例数据研究了噪声程度对模型阶次确定的影响。研究发现随着噪声强度的增加，对响应贡献较小的模态容易受到噪声干扰，此时模型的阶次的确定会受到影响。在5自由度系统中，第5阶模态对脉冲响应的贡献较小，在较大的噪声情况下被湮没，因此模型阶次确定中无法包含此阶模态;而对海洋平台数值算例则没有出现此种情况。

(3)利用2个数值算例进行了数据长度对模型阶次确定的研究。研究发现，对简单5自由度系统，随着数据量的增加，噪声对模型阶次的影响也在加剧。但对复杂的三维空间结构，发现数据长度对模型阶次确定的影响非常小。从应用的角度来说，建议可以采用几种不同的数据长度通过比较来确定模型的阶次。

［1］罗文波，夏松波，王 莉，等.对稳定图识别真假模态的一点改进［J］.应用力学学报，1999，16(2):112－116.

［2］常 军，孙利民，张启伟.基于两阶段稳定图的随机子空间识别结构模态参数［J］.地震工程与工程振动，2008，28(3):47－51.

［3］易伟建，刘 翔.动力系统模型阶次的确定［J］.振动与冲击，2008，27(11):12 －18.

［4］Sanliturk K Y，Cakar O.Noise elimination from measured frequency response functions［J］.Mechanical Systems and Signal Processing，2005，19:615 －631.

［5］张 波，李健君.基于Hankel矩阵与奇异值分解(SVD)的滤波方法以及在飞机颤振试验数据预处理中的应用［J］.振动与冲击，2009，28(2):162－167.

［6］Sau-Lon James Hu，Bao X X，Li H J.Model order determination and noise removalfor modalparameter estimation［J］.Mechanical Systems and Signal Processing，2010，24(6):1605–1620.

［7］吕志民，张武军，徐金梧.基于奇异谱的降噪方法及其在故障诊断技术中的应用［J］.机械工程学报，1999，35(3):85－88.

［8］杨文献，任兴民，姜节胜.基于奇异熵的信号降噪技术研究［J］.西北工业大学学报，2001，19(3):368－371.

［9］练继建，李火坤，张建伟.基于奇异熵定阶降噪的水工结构振动模态ERA识别方法［J］.中国科学E辑:技术科学，2008，38(9):1398 －1413.

［10］赵学智，叶邦彦，陈统坚.奇异值差分谱理论及其在车床主轴箱故障诊断中的应用［J］.机械工程学报，2010，46(1):100－108.

［11］赵学智，叶邦彦，陈统坚.基于奇异值曲率谱的有效奇异值选择［J］.华南理工大学学报，2010，38(6):11 －16.