李玉刚 叶庆卫 周 宇 方 宁宁波大学信息科学与工程学院，宁波,315211

基于模态参数提取的随机子空间辨识算法改进

李玉刚 叶庆卫 周 宇 方 宁
宁波大学信息科学与工程学院，宁波,315211

随机子空间辨识(SSI)算法在大型结构的振动检测、损伤识别中有着重要的作用。引入稀疏优化取代最小二乘法来获得尽可能稀疏的状态矩阵，引入K-means算法从众多模态参数中选出真实模态，以避免虚假模态的产生。实验结果表明，所构建的稀疏改进SSI算法能准确提取模态参数，对工程应用具有较大的参考价值。

随机子空间辨识算法；稀疏优化；最小二乘法：模态参数：K-means算法

0 引言

利用结构的动态响应识别结构损伤是近年来发展起来的结构损伤诊断新方法，而参数识别是结构健康监测领域中的重点。目前国内外模态参数的提取方法有频域分解法[1]、PolyMax[2]、NexT[3]等，但是这些算法需要大量工作来确定具有最小量参数的模型，并进行迭代运算，往往存在发散或收敛缓慢等问题。随机子空间识别算法 (stochastic subspace identification，SSI)[4]是基于随机状态空间模型的，只需要确定系统的阶次，且由于算法中的奇异值分解，故该算法不存在收敛问题。

随机子空间识别最关键的问题是确定系统阶次，目前主要的两种方法是奇异值跳跃法和稳定图法。上述两种方法的系统定阶过程都需要人工参与，奇异值跳跃法需要人工找出奇异值跳跃点，受环境噪声干扰较大；稳定图法需要人工找出系统阶次的范围，受人的主观判断影响较大，很容易出现虚假模态现象。

稀疏表示[5-7]是一种高效的搜索算法，在压缩感知中被广泛地应用。在压缩感知中，稀疏表示可以在不丢失原始信号信息的情况下，显著减少观测信号的次数，能有效实现信号的降维处理，具有很好的噪声抗干扰能力。稀疏优化采用稀疏逼近来取代原始信号，从而显著降低处理成本，特别对混合信号的分离有着显著的效果。现阶段的稀疏优化算法主要有梯度投影法[8]、匹配追踪法[9]、正交匹配追踪法(OMP)[10]等。

目前SSI算法能够在系统低阶情况下获取较高精度的模态主频参数和阻尼比，但是模态振型参数的提取还存在些许问题。另外SSI算法需要准确定阶，且受噪声的影响较大。本文在基于协方差的SSI方法(SSI-cov)[11-12]的基础上用稀疏优化方法对其求状态矩阵、输出矩阵环节作了改进，用稀疏优化技术代替最小二乘法，并通过K-means算法来选择系统模态，从而剔除系统虚假模态，减小系统阶次对模态参数的影响，使模态参数更加精确。

1 基于协方差驱动的随机子空间算法识别

在SSI-cov中，可以由采样数据得到可观矩阵Γ如下：

(1)

比较式(1)中两式可得：

Γ2=Γ1Θ

(2)

其中，Θ是对角矩阵;C为系统输出矩阵。然后通过最小二乘法求取状态矩阵A。

现有SSI-cov算法首先需要确定系统阶次，其系统定阶的方法主要是奇异值(SVD)分解法[13]，可以得到较为准确的系统阶次，但是当信号数据过大，外界环境产生的噪声比较复杂时，在确定系统阶次的时候往往会产生较大的偏差，很容易产生虚假模态。另外对输出矩阵C的求解只计算了可观矩阵Γ的第一行，可能会因为系统的复杂性而产生不可避免的偏差，导致结构模态振型出现紊乱。

2 SSI算法改进的稀疏求解

在现有的SSI-cov识别中，式(2)用最小二乘法求解状态矩阵A，然后求解模态频率ω*、阻尼比ξ*。这种方法对模型阶次的限制较大，一旦模型阶次确定不准确，虚假模态就会伴随而生；在以往求取输出矩阵C中，都是把可观矩阵Γ的第一行作为C的固定值，然而当噪声影响较大时，得到的C就会产生偏差，进而导致模态振型Ω*的求取精度出现较大的下降。

首先我们为模型阶次nΔ赋一个较大的值，把可观矩阵Γ转化为

可得出状态矩阵的等式方程：

ΓT=QTPT

(3)

对式(2)我们用高斯随机测量矩阵H1进行观测：

H1Γ2=H1Γ1Θ

(4)

将其转化成稀疏等价模型：

(5)

这个优化问题是个NP难问题，DONOHO指出，当H1Γ1满足RIP条件时，0范数的优化就可以转化成1范数的优化问题，即

(6)

式中，σ1为误差限度，是一个很小的正值。

对式(3)进行转置并用高斯随机测量矩阵H2观测：

H2ΓT=H2QTPT

(7)

构建稀疏优化模型：

(8)

同样，这也是一个NP难问题，当H2QT满足RIP条件时，0范数的优化就可以转化成1范数的优化问题：

(9)

式中，σ2为误差限度，是一个很小的值。

经验证，H1Γ1和H2QT能够以很高的概率满足RIP条件，符合稀疏求解的条件。但是我们通过式(6)和式(9)求取的稀疏解中有多于系统真实阶次的非零解，这些较小的非零解属于噪声，不过稀疏解中的噪声很小且极不稳定，通过多次运算可以很容易地统计找出，然后运用聚类的方法，对应得出真实模态参数，具体步骤如下。

(1)给系统的模型阶次nΔ赋一个较大的值，选取合适的稀疏度k和误差限度σ1、σ2。

(2)构建随机高斯测量矩阵H1和H2，对式(6)用OMP方法进行稀疏求解，得到矩阵Θ，进而得到矩阵A，再对A进行特征值分解，得到特征值D；得到A之后，构建QT矩阵，然后对式(9)用OMP进行稀疏求解，得到矩阵PT，观察PT的特点，就可以很容易地得到C。

(3) 由求得的A和C计算含有虚假模态的频率向量ω*、阻尼比向量ξ*和模态振型Ω*。

(4)重复步骤N次，统计结果，得到最终的D*和C*，统计步骤如下：①对于步骤(2)，我们得到N个D和C，对每一个D和C求其均值，然后把D与C的元素的绝对值与其各自的均值比较，小于均值的强制置零，然后对照N个D在相同位置的非零个数，记为数组dD，同理，对照C在相同位置的非零个数，记为数组dC；②把数组dD和dC的元素与其均值相减，把差值数据的绝对值记为fD和fC，由K-means算法，我们把数组fD和fC任选两个元素作为聚类中心，然后把剩下的元素与聚类中心比较，把离聚类中心较近的对象归为一类，然后计算所获得的新的聚类，不断重复，直至标准测度函数开始收敛，然后我们统计出均值小的类的下标索引gD和gC；③把下标gD和gC对应的D和C值强制置零，然后得到最终的D*和C*。

(5)根据D*中非零元素的位置，从ω*、ξ*中选择出系统模态ω、ξ；根据C*中非零元素的位置，从Ω*中选择出模态振型Ω。

通过上述方法，如图1所示，我们能够把模型阶次对模态参数提取的影响忽略不计，可以有效剔除虚假模态，并通过聚类的方法，很好地把噪声剔除，较大程度地减小噪声对结果的影响，从而提高算法的消噪能力和识别精度。

图1 SSI算法改进稀疏求解Fig.1 Improvement and solution of SSI algorithm based on sparse representation

3 仿真实验与分析

3.1 仿真实验

本文以MATLAB 作为仿真工具，这里假设一个4次系统的仿真信号：

x=0.4e-0.012tcos(5.3t+0.25)+1.1e-0.057t·cos(3.9t+0.1)-1.3e-0.082tcos(6.6t+0.34)-0.8e-0.034tcos(4.4t+0.18)

(10)

由理论计算可得到系统的固有频率ω、阻尼比ξ和模态振型系数Ω，如表1所示。

表1 系统固有模态参数
Tab.1 Inherent modal parameters of system

ω(Hz)ξ(N·S/m)Ω0.8440.0140.4000.6210.0921.1001.0510.078-1.3000.7010.049-0.800

对式(10)的仿真信号添加信噪比为5 dB的高斯白噪声，以30 Hz的采样频率采集1000个采样点，然后我们为系统阶次n赋值80，然后计算频率向量ω*、阻尼比向量ξ*、模态振型Ω*以及确立等式方程：

(Γ2)200×80=(Γ1)200×80Θ80×80

(11)

(12)

按照如下过程分步计算，最终可得到各阶振型系数及模态。

(2)同样地，构建高斯随机测量矩阵(H2)24×200，对式(12)进行观测，然后由式(9)得到矩阵C1×80。

(3) 重复步骤(1)8次，得到这8次计算的解集{A1,A2,…,A8}，计算特征解集{D1,D2,…,D8},对每一个特征值求取均值，然后将其元素的绝对值与其均值比较，小于均值的强制置零。然后对照8个特征值在相同位置的非零个数，记为数组dD。

(4) 重复步骤(2)8次，得到这8次计算的解集{C1,C2,…,C8},对每一个输出矩阵求取均值，然后将其元素的绝对值与其均值比较，小于均值的强制置零。然后对照8个特征值在相同位置的非零个数，记为数组dC。

(5)将数组dD元素与其均值相减，差值数据的绝对值记为fD，由K-means算法统计出均值小的类的下标索引gD，排除下标索引gD之后，得出余下的下标2、7、24、37，从而求取特征值D*，通过D*就可以找出ω*和ξ*对应的位置。

(6) 将数组dC元素与其均值相减，差值数据的绝对值记为fC，由K-means算法我们统计出均值小的类的下标索引gC，排除下标索引gC之后，得出余下的下标3、13、25、42，从而求输出矩阵C*，通过C*就可以找出对应的模态振型Ω*。

(7)步骤(5)与步骤(6)得出的D*、ω*、ξ*、C*、Ω*如表2所示。

表2 4阶次系统的稀疏改进的模态参数
Tab.2 Modal parameters of sparse improvement based on 4 model order

D*ω*(Hz)ξ*(N·S/m)C*Ω*0.9950.8430.0131.5340.4320.9170.6240.0945.4621.1231.4841.0540.0800.948-1.3050.9340.7030.0510.859-0.812

通过表2，我们可以看出稀疏改进算法在模型阶次过大和高斯白噪声干扰的情况下，仍旧保证了很高的计算精度，很好地证明了稀疏改进算法的优越性。

3.2 模型阶次对系统的影响

3.1节的稀疏求解具有一定的片面性，模型阶次和稀疏度的选取较为主观，所以为了使验证结果更加准确，我们需要再次验证。首先选用两种差距比较明显的较大模型阶次，然后对每个模型阶次下的稀疏度选取不同的值，以信噪比5 dB为干扰信号，各计算80次求均值，得到固有频率见表3，阻尼比见表4，模态振型得到结果见表5。

表3 不同阶次稀疏计算固有频率Tab.3 Calculation of modal frequency of different order based on sparse algorithm Hz

表4 不同阶次稀疏计算阻尼比Tab.4 Calculation of damping ratio of different order based on sparse algorithm N·S/m

表5 不同阶次稀疏计算模态振型
Tab.5 Calculation of vibration pattem of different order based on sparse algorithm

阶次稀疏度模态振型Ω1模态振型Ω2模态振型Ω3模态振型Ω450240.3981.098-1.295-0.79350480.4111.107-1.308-0.811100240.4031.104-1.302-0.803100480.4091.106-1.305-0.809

表3～表5显示，两种相差较大的模型阶次计算结果基本一样，且对稀疏度的敏感程度很低，这对我们稀疏度的选取提供了一个较大的范围，因此，我们可以很容易地得出，稀疏改进算法对模型阶次要求不高，可以很好地在剔除虚假模态的基础上保持计算精度。

3.3 不同噪声下的模态识别精度

为了更好地验证稀疏改进算法的优越性，我们把稀疏改进算法与SSI-cov算法进行对比。

(13)

分别对信号添加不同的加性噪声(脉冲噪声和高斯白噪声以及两种噪声的混合噪声)，为了使结果更加精确，我们计算80次求其平均值，比较两种算法在相同噪声下对信噪比的敏感程度。因篇幅有限，这里只列举4个比较典型的误差对比。如图2～图5所示。

图2 基于高斯白噪声下频率误差对比Fig.2 Comparison of frequency error based on gauss white noise

图3 基于混合噪声下频率误差对比Fig.3 Comparison of frequency error based on mixed noise

图4 基于高斯白噪声下阻尼比误差对比Fig.4 Comparison of damping ratio error based on gauss white noise

图5 基于高斯白噪声下模态振型误差对比Fig.5 Comparison of vibration pattem error based on gauss white noise

从图2中可以看出，在高斯白噪声下，SSI-cov识别的频率误差起伏较大，在γSN为0～10 dB范围内,波动特别明显，可看出SSI-cov算法对噪声较为敏感，而稀疏改进算法起伏波动比较平缓，有很高的稳定性；从图3中可以看出，在混合噪声下，SSI-cov识别的主频精度比在高斯白噪声下相对较差，对噪声更加敏感，在信噪比15 dB以内，相对误差变化起伏更大，而稀疏改进算法相对误差变化平缓；从图4中可以看出，SSI-cov识别的阻尼比在信噪比5 dB以内，已经产生了明显的失真，当信噪比越小时，识别精度越差，而稀疏改进算法识别精度较高，表现出了很高的稳定性；从图5中可以看出，稀疏改进算法的相对误差始终在SSI-cov之下，特别当信噪比较小时，误差值更为明显。

4 实测信号分析

本文研究对象的工程数据采自宁波某斜拉索大桥。大桥全长67 m，由102根直径为0.15 m的拉索构成拉索支撑系统。在采集过程中，采用WS-ZHT2振动设备和双传感器采集振动信号，双传感器安装在拉索和梁端的铰支部位，能有效感应索-梁耦合的拉索振动，如图6所示。

图6 传感器安装图Fig.6 Sensors installation diagram

传感器采集到的信号都是含噪信号，但是我们不知道噪声的类型是否规则，不清楚无噪声原始信号的波形，只清楚采集到的信号波形。这里我们以第6号到第12号斜拉索的振动采集信号为例，因篇幅有限，只显示第10号斜拉索的振动波形图，如图7所示，信号具有较强的环境干扰，从而对模态参数提取的精度和噪声抗干扰性要求较高。用本文算法对采集到的信号提取基频参数，与SSI-cov算法作比较，来进一步验证稀疏改进SSI算法的优越性。

图7 第10号斜拉索振动波形图Fig.7 Vibration waveform of No.10 stay cable

在图7所示波形图的基础上，两种算法的测定结果(测定10次，求其均值)如表6所示。

表6 两种识别算法对斜拉索基频识别情况
Tab.6 Frequency identification of two algorithms for the stay cable

斜拉索号参考基频(Hz)识别方法计算结构基频(Hz)相对误差A60.7528稀疏改进SSI-cov0.75800.76520.00520.0124A80.7863稀疏改进SSI-cov0.79660.80380.01030.0175A100.8011稀疏改进SSI-cov0.80920.81530.00810.0142A120.7932稀疏改进SSI-cov0.80060.80580.00740.0126

从表6中数据可以看出，稀疏改进算法识别具有更好的精确度，相对误差普遍较小，更加接近参考基频(大桥管理服务有限公司提供)。稀疏改进算法解决了SSI-cov当阶次确定有偏差的情况下，识别精度不准确的问题，另外稀疏改进算法对噪声不敏感，本身在稀疏优化的过程中，就已经具有了很好的消噪能力。因此，我们得出如下结论：稀疏改进算法对系统阶次的要求不高，即使出现定阶错误对结果影响也可忽略，以及具有很好的抗噪性，在较大的噪声干扰下，仍然具有良好的识别精度。

5 结论

综上所述，稀疏改进SSI算法具有工程可行性。相较于经典的SSI算法的最小二乘法，稀疏优化的OMP算法在有噪声或噪声较大时，能够较好地提取状态矩阵和输出矩阵，准确提取模态参数，表现出更好的抗噪声性能，鲁棒性有明显的提升；相较于经典的SSI算法的系统定阶，稀疏改进SSI的聚类算法在系统阶次方面几乎不需要定阶，有效地解决了因模型阶次过高而产生的虚假模态问题，扩大了模型阶次的选择范围。

在仿真实验中，通过与SSI-cov算法的对比，很好地验证了稀疏改进SSI算法的抗噪性以及解决模型阶次过高产生的虚假模态的能力；对宁波某大桥的斜拉索的基频测定，更加有力地证明了稀疏改进算法的优越性。因此，稀疏改进SSI算法具有一定的工程应用价值。

[1] BRINCKER R, ZHANG L, ANDERSEN P. Modal Identification from Ambient Responses Using Frequency Domain Decomposition [C]//18thIMAC. San Antorio. Texas,2000:65-630.[2] GUILLAUME P, VERBOVEN P, VANLANDUIT S. Frequency-domain Maximum Llikelihood Identification of Modal Parameters with Confidence Intervals [C]//In Proceedings of ISMA 23, the International Confere-nce on Noise and Vibration Engineering. Leuven, Belgium,1998:16-18.

[3] JAMES G H, CARNE T G, LAUFFER J P. The Natural Excitation Technique for Modal Parameter Extraction from Operating Structures [J]. Modal Analysis: Int. J. Analytical and Experimental Modal Analysis,1995,10(4):260-227.

[4] PEETERS B, ROECK G D. Reference-based Stochastic Subspace Identification for Output-only Modal Analysis [J].Mechanical Systems and Signal Processing,1999,13(6):855-878.

[5] DONOHO D L. Compressed Sensing [J]. IEEE Transactions on Information Theory,2006,52(4):1289-1306.

[6] 宋欢欢，叶庆卫，王晓东，等. 基于稀疏AR建模信号去噪研究与应用[J].振动与冲击,2015,34(6):127-131. SONG Huanhuan, YE Qingwei, WANG Xiaodong, et al. Study and Application of Signal Denoising Based on Sparse AR Model[J]. Journal of Vibration and Shock,2015,34(6):127-131.

[7] HAUPT J, BAJWA W U, RABBAT M, et al. Compressed Sensing for Networked Data [J]. IEEE Signal Processing Magazine,2008,25(2):92-101.

[8] BLUMENSATH T, DAVIES M E. Gradient Pursuits [J]. IEEE Trans. on Signal Processing,2008,56(6):2370-2382.

[9] MALLAT S G, ZHANG Z. Matching Pursuits with Time-frequency Dictionaries [J]. IEEE Transactions on Signal Processing,1993,41(12):3397-3415.

[10] TROPP J. Greed Is Good: Algorithmic Results for Sparse Approximation [J]. IEEE Trans. on Information Theory,2004,50(10):2231-2242.

[11] BOONYAPINYO V, JANESUPASAEREE T. Data-driven Stochastic Subspace Identification of Flutter Derivatives of Bridge Decks [J].Journal of Wind Engineering and Industrial Aerodynamics, 2010,98(12):784-799.

[12] 李永军，马立元，王天辉，等.协方差驱动子空间模态参数辨识方法改进分析[J].中国机械工程，2012,23(13):1533-1536. LI Yongjun, MA Liyuan, WANG Tianhui, et al. An Improvement on Subspace Modal Parameter Identification Algorithm Driven by Covariance [J]. China Mechanical Engineering,2012,23(13):1533-1536.

[13] 赵学智，叶邦彦.分量形成方式对奇异值分解信号处理效果的影响[J].上海交通大学学报,2011,45(3):368-374. ZHAO Xuezhi, YE Bangyan. The Influence of Formation Manner of Component on Signal Processing Effect of Singular Value Decomposition[J]. Journal of Shanghai Jiaotong University,2011,45(3):368-374.

[14] CANDES E J, TAO T. Decoding by Linear Programming[J].IEEE Transactions on Information Theory,2010,34(4):435-443.

[15] 张成，杨海蓉，韦穗.基于随机间距稀疏Toeplitz测量矩阵的压缩传感[J].自动化学报，2012,38(8):1362-1369. ZHANG Cheng, YANG Hairong, WEI Hui. Compressive Sensing Based on Deterministic Spar-se Toeplitz Measurement Matrices with Random Pitch [J]. Acta Automatica Sinica,2012,38(8):1362-1369.

(编辑 袁兴玲)

Improvement of SSI Algorithm Based on Extraction of Modal Parameters

LI Yugang YE Qingwei ZHOU Yu FANG Ning

Faculty of Information Science and Enginner，Ningbo University， Ningbo，Zhejiang,315211

SSI algorithm played an important role in the large structure vibration detection and damage identification. Sparse optimization solution was introduced to replace the least square method that was used to get sparser state matrix. K-means algorithm was introduced to elect real modal parameters from many modal parameters so as to eliminate the false modals effectively. The experimental results show that optimization solution of SSI algorithm may accurately extract modal parameters. The work herein has reference values in engineering applications.

stochastic subspace identification (SSI) algorithm; sparse optimization; least square method; modal parameter; K-means algorithm

2016-07-20

国家自然科学基金资助项目(51675286，61071198)；浙江省自然科学基金资助项目(LY13F010015)；宁波市自然科学基金资助项目(2012A610019)；浙江省科技创新团队资助项目(2013TD21)

TP391.4

10.3969/j.issn.1004-132X.2017.01.012

李玉刚，男，1991年生。宁波大学信息科学与工程学院硕士研究生。主要研究方向为振动信号处理。E-mail:630818644@qq.com。叶庆卫，男，1970年生。宁波大学信息科学与工程学院副教授。周 宇，男，1962年生。宁波大学信息科学与工程学院教授。方 宁，1992年生。宁波大学信息科学与工程学院硕士研究生。