标题

Size: px

Start display at page:

Download "标题"

澄要森
5 years ago
Views:

1 2016 年 11 月 November 2016 第 35 卷第 11 期 Mechanical Science and Technology for Aerospace Engineering Vol.35 No.11 DOI: / j.cnki 一种 AR 模型的稀疏模态参数估计方法李仙果, 周宇, 叶庆卫, 王晓东 ( 宁波大学信息科学与工程学院, 浙江宁波 ) 摘要 : 基于自回归模型 (autoregressive,ar) 的模态参数估计方法需要进行阶次估计, 阶次确定错误会使结果中包含虚假模态或者遗漏模态针对这种情形, 通过将模态选择的问题转化为带有稀疏约束的振型系数问题, 提出了一种基于 AR 模型的稀疏模态参数估计方法先假设一个较大的系统阶次, 采用基于 AR 模型的估计方法计算模态参数, 然后构造振型系数的稀疏优化模型, 多次稀疏求解后用统计和聚类的方法对结果进行校正, 最后根据振型系数鉴定对应模态的真假这种无需提前确定阶次的估计方法既避免了定阶次引起的误差, 又能有效地将噪声模态分离, 适合复杂环境下的参数估计在斜拉索桥模型上验证结果表明, 该方法具有良好的识别效果关键词 : 自回归模型估计优化统计振型系数中图分类号 :TP391 4 文献标识码 :A 文章编号 : (2016) A Sparse Modal Parameter Estimation Method Based on AR Model Li Xianguo 1, Zhou Yu 1, Ye Qingwei 1, Wang Xiaodong 1 ( Faculty of Electrical Engineering and Computer Science, Zhejiang Ningbo , China) Abstract: The modal parameter estimation method based on the autoregressive ( AR ) modal may produce inaccurate modal parameters as well as a fake modal or miss some modals because of wrong orders. To solve this problem, this paper comes up with a new method by solving the coefficients of vibration mode with a sparse constraint so as to produce accurate modal parameters. It first uses the AR model to arrive at the modal parameters by assuming a larger order for the vibration mode and then constructs the sparse optimization model for the coefficient of vibration mode. After that, the statistical results obtained with the multiple solutions for the sparse optimization model are used to adjust the final results. Finally, the paper identifies modal parameters according to the coefficient of vibration mode. This method not only avoids the errors caused by assuming the large order but also enhances noise resistance. The test results on the cable demonstrate the stability and accuracy of this method. Keywords: autoregressive model modal parameter estimation optimization statistical method coefficient of vibration mode 参数识别是结构健康监测领域研究的重点大部分模态参数估计方法都需提前确定系统阶次, 阶次确定的准确与否对参数识别结果有很大影响曲 [1] 线拟合法奇异值方法稳定图方法等方法在大收稿日期 : 基金项目 : 国家自然科学基金项目 ( ) 资助作者简介 : 李仙果 (1987-), 硕士研究生, 研究方向为振动信号模态参数提取,lixianguo2014@ foxmail.com 周宇 ( 联系人 ), 教授, 硕士,zhouyu@ nbu.edu.cn 多数情况下都能准确估计出系统的阶次, 但是在噪声较复杂时可能出现估计错误, 导致估计的模态参数出现遗漏模态或者包含虚假模态的现象, 因此定阶增加了算法的不稳定性稳定图法是一种既能够确定系统阶次, 又能从诸多模态中鉴别出真假模态的算法, 具有诸多优点, 可在强噪声条件下进行模态识别但是该方法在应用中发现获得的一些极少数的极点并不是系统本身的特性, 这导致稳定图法识别的参数中包含了虚假模态, 不少学者对此做了改进与应用 [2] 另外一种重要的模态参数估计方法

2 1734 第 35 卷是基于 AR 模型的模态参数估计方法, 该方法需要提前定阶, 对噪声敏感, 但它需要的数据量小, 计算速度快, 因此成为很多新方法的基础利用奇异值进行定阶和去噪的方法 [3], 该方法是通过对 AR 基矩阵奇异值分解, 将特征值截断后恢复出新的基矩阵, 然后根据基矩阵本身的特点进行反对角线平均, 如此反复迭代直到特征值不再变化该方法具有较强的定阶和去噪能力, 大大增加了基于 AR 模型的模态参数估计方法的应用范围 [4 6] 稀疏表示往往能够简化一个信号的表示, 便于信号的传输, 或者反映出信号的一些稀疏特征, 为研究人员提供一些隐含的信息信号的稀疏表示 [7] 总是伴随着稀疏优化问题, 通过稀疏优化获得稀疏系数是稀疏表示的核心内容 1 稀疏优化算法当一个信号中大部分元素为零, 只有少数元素为非零元素时, 这个信号被称为稀疏信号文献 [8] 对稀疏信号自身的特殊性进行了详细研究, 其中稀疏优化就是针对这一类信号的一种优化算法通常直接采集得到的稀疏信号较少, 但是大部分非稀疏的信号在其他变换域中是稀疏的因此可以将这类信号在其变换域中进行稀疏表示, 然后再用稀疏优化算法来达到某种优化目的假设一带噪声的信号 y 可以表示为 y = D α + w (1) 式中 : α 为干净信号 s = D α 的稀疏分解系数向量, 具有稀疏的特征 D 为投影字典 w 为噪声如果信号 y 携带的高斯白噪声 w 满足 w 2 σ, 则式 (1) 的数学模型可以改写成稀疏优化模型, 即 min α { 0 (2) subject to y - D α 2 ε 当 σ ε 时上式存在稀疏解上式的求解方法 [9] 主要有贪婪算法和凸优化方法 [10] 相比凸优化方法, 贪婪算法的计算复杂较低, 计算速度快, 在实际中得到了更多的应用, 而在贪婪算法中以正交匹配追踪 [11] (OMP) 算法应用最为广泛用 OMP 算法求解信号稀疏优化问题时通常要知道信号的稀疏度 [12], 当稀疏度未知时 OMP 算法可以采用如下步骤 : Step1 典 D, 误差限度 δ Step2 设定算法输入 : 带噪声信号 y, 过完备字初始化参数 : 残差 r 0 = y, 分解系数 α^ = 0, 索引集 t 0 = [], 子字典 T 0 = [], 迭代因子 f = 1 Step3 1), 运行以下步骤 : 在迭代过程中, 在第 f 次循环时 ( f 1) 相关性最大计算, 寻找最佳原子索引 i f = arg max d i,r f -1,i = 1,2,,N i 式中 : arg 表示取元素的索引 d i 表示字典 D 中的第 i 列向量 2) 更新原子索引集合, t f = t f -1 Ui f 3) 更新子字典,T f = T f -1 Ud if 4) 更新系数估计,α^ [t f ] = T + f y, 其中 T + f = (T H f T f ) T H f 5) 更新残差,r f = y - T f (T + f y) 6) 判断终止条件, 若 r f 2 δ, 迭代结束 Step4 输出稀疏的系数 : α^ = α^ 按照以上步骤的 OMP 算法得到的稀疏解并不唯一 [13], 其随着误差限度 δ 的不同会有所差异 2 基于 AR 模型的模态参数估计已知系统阶次为 n, 脉冲响应序列在第 k 点的值为 p k (k = 0, 2,, m + 2n), 其时间间隔为 dt, 且 m + 1 > 2n, 构造得到 (m + 1) 个样本的自回归方程, 组成关于 a 的线性方程组, 即 C a = b (3) 式中 : C 是 (m + 1) 2n 阶实矩阵, 称为 Hankel 矩阵 b 是 (m+1) 阶实列阵 é p 0 p 1 p 2n ù p 1 p 2 p 2n+1 C = (4) ëp m p m+1 p m+2n-1 û b = [ p 2n+1 p 2n+2 p m+2n ] T (5) 求解式子 (4) 的最小二乘解得 a = (C T C) -1 C T b (6) 此时自回归系数 a 与信号的模态频率 ω 和模态阻尼比 ξ 以及模态振型 γ 存在一定关系 [14], 按照如下步骤可以求得模态频率及模态阻尼比 1) 将自回归系数中的元素表示为 a = (a 0, a 1,, a 2n-1 ), 然后构造 prony 多项式 x 2n + a 2n-1 x 2n a 0 = 0, 求出多项式的根, 尼比 ξ, x^ = (x 0,x 1,,x 2n-1 ) (7) 2) 根据多项式的根, 计算模态频率 ω 和模态阻 ìu = ln(x^) ω = U / dt í 1 ξ = î 1 + (Im(U) / Re(U) 2 (8)

3 第 11 期李仙果等 : 一种 AR 模型的稀疏模态参数估计方法 1735 式中 :Im() 表示取括号内的虚部,Re( ) 表示取括号内的实部 3) 构造关于模态振型 γ 的方程 X = X γ = g (9) é x 0 0 x 0 1 x 0 2n-1 ù x 1 0 x 1 1 x 1 2n ë x j 0 x j 1 x j 2n-1 û (10) g = [p 1 p 2 p 2n-1 ] T (11) 在式 (10) 中 j n, 计算式 (9) 的最小二乘解得到模态振型, 即 γ = (X T X) X T g (12) 这种基于 AR 模型的模态参数估计方法对噪声敏感 [15], 并且需要知道系统的阶次, 当系统阶次估算不准确时会出现遗漏模态或包含虚假模态组成一个矩阵 Ψ d = [γ d1, γ d2,, γ dp ] 4) 对 Ψ d 的行分别进行 K means 聚类 [16], 每行的元素被分为两类, 将元素少的一类置零, 得到新的矩阵 Ψ s 5) 统计矩阵 Ψ s 中的每行中非零元素的数量, 得到统计数组 η, 对 η 聚类分为两类, 将均值绝对值小的一类对应的 Ψ s 中的行置零, 得到 Ψ k 6) 计算 Ψ k 中每行非零元素的均值 ( 元素全部为 0 的行的均值为 0), 得到稀疏的 γ d 根据 γ d 中非零元素的位置从 ω d 中找出系统的模态,γ d 中非零元素为其对应的振型系数系统模态选择流程图如图 1 所示 3 稀疏的模态参数估计 3 1 模态参数的稀疏选择方法首先用 L (L n) 次的自回归方程计算带虚假模态的模态数组 ω d 和振型方程组, 即 X γ d = g (13) 式中 : X 是 h L 的实矩阵 (h > L)g 为实列阵由于 X 的列数远大于 n, 它的很多列对信号 g 并没有贡献, 属于噪声列, 其在 γ d 中对应的系数为 0 因此 γ d 的大多数系数为 0, 可看作是一个稀疏的系数因此可将式 ( 13) 转化为带有稀疏约束的优化模型, 即 : min γ d { 0 (14) X γ d - g δ 式中 δ 为误差限度将 γ d 是信号 g 在冗余字典 X 上的稀疏投影系数, 该稀疏优化模型可以用 OMP 算法求出稀疏的 γ d, 但是在不知道 γ d 的稀疏度的情况下, 用上文 OMP 算法计算的结果并不唯一, 而且精度不高下面采用多次观测求解, 并结合聚类方法对求得的稀疏系数做进一步的优化, 提高其稳定性和精度优化稀疏系数的求解步骤如下 : 1) 构造 L h 阶的高斯随机观测矩阵 G, 用 G 对式子 (13) 进行观测 : G X γ d = G g (15) 2) 确定误差限度 δ, 将式 (15) 转化为带稀疏的优化模型, 并用 OMP 算法估计稀疏的列向量 γ d 3) 重复步骤 1) 2) p 次, 得到 γ d 的 p 个估计, 图 1 系统模态选择流程图这里通过求稀疏的振型系数的方法来判断真假模态, 避免了系统阶次的确定同时稀疏优化算法的引入也提高了算法的抗噪声能力 3 2 稀疏模态参数估计方法仿真假设某一粘性阻尼振动信号的模型为 y 1 = 0 85e t cos(1 885t) + 1 1e -0 04t cos(6 346t) + 0 6e -0 03t cos(3 1416t) (16) 计算出这条信号的模态固有频率为模态阻尼比为振型系数为现在按照离散时间间隔 0 5 s 取 200 个点得到信号 y 2, 对信号 y 2 按照 8 db 添加高斯白噪声, 得到信号 y n 用信号 y n 建立一个 100 阶的自回归方程组, 计算得到含虚假模态的模态数组 ω d, 它们各包含 100 个元素, 同时得到对应的振型方程组为 X γ = g (17) 通过用高斯随机测量矩阵对该方程组进行 8 次观测并分别用 OMP 算法进行稀疏求解, 得到这 8 次计算的稀疏解组成的矩阵 Ψ, 由于数据量较大, 现将矩阵 Ψ 前 5 列的部分结果列出如表 1

１７３６ Ψ １表１稀疏矩阵 Ψ 的部分结果 Ψ ２ Ψ ３ Ψ ４ Ψ ５７６２５８５３８７５９８１２３第３５卷８３１５３７５１２３９７９９５４１１１８８８２５５３６１２５８３６中方法１是稀疏改进的模态参数估计方法使用前未对信号做任何处理方法２是基于ＡＲ模型的模态参数估计方法使用前用奇异值方法进行定阶和去噪声处理

4 １７３６ Ψ １表１稀疏矩阵 Ψ 的部分结果 Ψ ２ Ψ ３ Ψ ４ Ψ ５７６２５８５３８７５９８１２３第３５卷８３１５３７５１２３９７９９５４１１１８８８２５５３６１２５８３６中方法１是稀疏改进的模态参数估计方法使用前未对信号做任何处理方法２是基于ＡＲ模型的模态参数估计方法使用前用奇异值方法进行定阶和去噪声处理５５９６８１２１２通过行聚类将每一行元素中差异明显的且属于少数的元素置零例如表中的第３行中的第３个元素第５行中第２第５个元素统统被置零得到新矩阵 Ψｄ然后统计 Ψｄ每行中非零元素的数量得到统计数组 η 对 η 用Ｋｍｅａｎｓ聚类方法分为两类表１中第２４６行对应的 η 中的元素被分为类图３１其他的被分为类２因为类２均值的绝对值小于类１所以类２在 Ψｄ中对应的行置零最终保留的就是我们要的稀疏解矩阵 Ψｓ计算 Ψｓ中每一行的均值得到稀疏的振型系数数组 γ 其中非零元素分别为８１５５７８１１２４根据数组 γ 中非零元素的位置从模态数组 ω ｄ ξ ｄ中筛选出固有模态频率 ω 为３５９８９固有阻尼比 ξ 为３１５３４ γ 中非零元素则为其对应的振型系数通过与信号真实参数的比较可以看出在高斯白噪声信噪比为８ｄｂ的条件下稀疏改进的参数识别算法可以得到准确的模态参数为了更好地表现本算法的效果现在分别采用基于ＡＲ模型的模态参数估计和稀疏改进的模态参数估计两种方法估计参数然后分别用所估计的参数恢复出响应信号通过比较恢复信号的信噪比来评判算法的优劣原始响应信号与加噪声后的信号对比图如图２所示添加噪声为高斯白噪声与脉冲噪声的混合噪声其信噪比为８８ｄｂ图２恢复信号对比图图３显示方法１恢复信号的信噪比为１９２ｄｂ方法２恢复信号的信噪比为１５３ｄｂ可见在高斯白噪声与随机脉冲噪声的混合情况下无定阶的参数识别方法能够恢复出信号具有更高的信噪比因此具有更好的抗噪声性能４斜拉索桥模型的实例研究利用斜拉索实验模型验证该稀疏改进的模态参数估计方法的有效性其中斜拉索模型如图４所示该模型的理论固有频率为７６２６３８９９固有阻尼比为２８２５４利用加速图４斜拉索桥模型图度传感器对该模型采样２个点采样频率为５Ｈｚ得到的采样信号去均值后如图５所示混合噪声下的信号波形图图３为２种方法识出的参数恢复出的信号其ｈｔｔｐｗｗｗｎｗｐｕｊｏｕｒｎａｌｓｏｒｇｃｎ图５斜拉索桥模型的采样信号

5 第 11 期李仙果等 : 一种 AR 模型的稀疏模态参数估计方法 1737 取图 5 信号的前 100 个点进行稀疏改进的模态参数估计, 估计结果如表 2 所示表 2 斜拉索桥模型的模态参数估计结果序号理论值估计值相对误差频率 % 阻尼比 % 频率 % 阻尼比 % 频率 % 阻尼比 % 通过表 2 数据对比可以明显的看出, 稀疏改进的模态参数估计方法估计的参数与桥梁的理论参数非常接近, 其中对于频率的识别误差在 1% 以内, 阻尼比的识别误差相对较大但不超过 10% 因此该稀疏改进的参数估计方法是有效的 5 结论软件仿真表明, 稀疏改进方法能够在复杂噪声条件下准确估计出系统参数, 而且恢复出的信号比传统的基于 AR 模型的模态参数估计方法恢复出的信号有更高的信噪比斜拉索桥模型的模态估计试验再次表明, 改进算法可以在没有定阶的情况下识别出系统的模态参数该稀疏改进方法对真假模态的鉴别是基于振型系数的稀疏解, 而稀疏解的不唯一性使结果稳定性不高尽管采样聚类优化的方法进行了处理, 但在极少数情况下仍然有出现噪声模态或者遗漏系统模态的可能如果能同时将模态参数本身的估计问题转化为稀疏求解问题, 那么结果中出现噪声模态或者遗漏系统模态的可能性就会进一步下降 [ 参考文献 ] [1] Fan J L, Zhang Z Y, Hua H X. Data processing in subspace identification and modal parameter identification of an arch bridge[ J]. Mechanical Systems and Signal Processing, 2007,21(4): [2] 常军, 张启伟, 孙利民. 稳定图方法在随机子空间识别模态参数中的应用 [J]. 工程力学,2007,24(2):39 44 Chang J, Zhang Q W, Sun L M. Application of stabilization diagram for modal parameter identification using stochastic subspace method [ J ]. Engineering Mechanics, 2007,24(2):39 44 ( in Chinese) [3] 曾小军, 黄宜坚. 基于 AR 模型和支持向量机的故障诊断法 [ J]., 2010, 29 ( 7): ,980 Zeng X J, Huang Y J. Fault diagnosis based on AR model and support vector machine [ J ]. Mechanical Science and Technology for Aerospace Engineering, 2010,29(7): ,980 ( in Chinese) [4] Donoho D L, Elad M. Optimally sparse representation in general ( non orthogonal ) dictionaries via l 1 minimization[ J]. Proceedings of the National Academy of Sciences of the United States of America, 2003,100 (5): [5] 戴琼海, 付长军, 季向阳. 压缩感知研究 [ J]. 计算机学报,2011,34(3): Dai Q H, Fu C J, Ji X Y. Research on compressed sensing[ J]. Chinese Journal of Computers, 2011, 34 (3): ( in Chinese) [6] Tsaig Y, Donoho D L. Extensions of compressed sensing[j]. Signal Processing, 2006,86(3): [7] Ye Q W, Sun Y, Wang X D, et al. An improved LLE algorithm with sparse constraint [ J ]. Journal of Computational and Theoretical Nanoscience, 2013, 10 (12): [8] Elad M, Bruckstein A M. A generalized uncertainty principle and sparse representation in pairs of bases [ J]. IEEE Transactions on Information Theory, 2002, 48 ( 9 ): [9] Mallat S G, Zhang Z F. Matching Pursuit with time frequency dictionaries[ J]. IEEE Transactions on Signal Processing, 1993,41(12): [10] Shanno D F, Vanderbei R J. Interior point methods for nonconvex nonlinear programming: orderings and higher order methods [ J]. Mathematical Programming, 2000, 87(2): [11] Chen S S, Donoho D L, Saunders M A. Atomic decomposition by basis pursuit [ J ]. SIAM Review, 2001,43(1): [12] Donoho D L. Compressed sensing[j]. IEEE Transactions on Information Theory, 2006,52(4): [13] Candes E J, Romberg J, Tao T. Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information[ J]. IEEE Transactions on Information Theory, 2006,52(2): [14] 曹树谦, 张文德, 萧龙翔. 振动结构模态分析 : 理论实验与应用 [M]. 天津 : 天津大学出版社,2001: Cao S Q, Zhang W D, Xiao L X. The modal analysis theory, experiment and application of structure vibration [ M ]. Tianjin: Tianjin University Press, 2001: ( in Chinese) [15] Huang C S, Su W C. Identification of modal parameters of a time invariant linear system by continuous wavelet transformation [ J ]. Mechanical Systems and Signal Processing, 2007,21(4): [16] Brohée S, Van Helden J. Evaluation of clustering algorithms for protein protein interaction networks [ J]. BMC Bioinformatics, 2006,7:488

Vol. 36 ( 2016 ) No. 6 J. of Math. (PRC) HS, (, ) :. HS,. HS. : ; HS ; ; Nesterov MR(2010) : 90C05; 65K05 : O221.1 : A : (2016)

Vol. 36 ( 2016 ) No. 6 J. of Math. (PRC) HS, (, ) :. HS,. HS. : ; HS ; ; Nesterov MR(2010) : 90C05; 65K05 : O221.1 : A : (2016) Vol. 36 ( 6 ) No. 6 J. of Math. (PRC) HS, (, 454) :. HS,. HS. : ; HS ; ; Nesterov MR() : 9C5; 65K5 : O. : A : 55-7797(6)6-9-8 ū R n, A R m n (m n), b R m, b = Aū. ū,,., ( ), l ū min u s.t. Au = b, (.)