CN104200002B - 一种粘性阻尼振动信号中的模态参数提取方法 - Google Patents

Abstract

本发明公开了一种粘性阻尼振动信号中的模态参数提取方法，其首先对振动信号进行奈奎斯特均匀采样；然后利用采样信号构建一个自回归矩阵方程，计算自回归矩阵方程中的自回归系数向量的最小二乘解，利用最小二乘解中的所有元素作为系数构造一个Prony多项式，解多项式得到变量的解向量；接着对解向量进行修正，根据修正后的解向量计算含有虚假模态的固有频率向量和固有阻尼比向量；之后以修正后的解向量为基础构造一个矩阵，再通过对采样信号在该矩阵上的投影方程进行稀疏优化求解，获取稀疏向量；最后根据稀疏向量剔除虚假模态，而稀疏向量中的所有非零元素为振型系数；本方法的优点是抗噪声能力强，且提取出的模态参数的精确度高、稳定性好。

Description

一种粘性阻尼振动信号中的模态参数提取方法

技术领域

本发明涉及一种振动信号处理技术，尤其是涉及一种粘性阻尼振动信号中的模态参数提取方法。

背景技术

振动信号中的模态分析在结构的动力特性分析、健康监测、结构探伤等领域具有重要的应用价值，其中模态参数的提取处于基础地位，如何精确提取模态参数一直是学者们关心的问题。当前有很多模态参数提取方法，针对粘性阻尼振动信号中的模态参数提取，最为简单方便、应用较多的是单输入单输出(SISO)方法，具有代表性的是复指数(CE)法，很多方法都是以复指数法为基础建立的。复指数法计算简单，需要的数据量少，当振动信号的噪声低时具有很高的提取精度，但随着噪声的增加该方法将不再适用，因此在实际应用中需要事先对振动信号进行消除噪声处理，然后再用复指数法进行模态参数提取，然而如果事先的噪声消除效果不好，则用复指数法提取出的模态参数的误差会很大。

发明内容

本发明所要解决的技术问题是提供一种粘性阻尼振动信号中的模态参数提取方法，其抗噪声能力强，且提取出的模态参数的精确度高、稳定性高。

本发明解决上述技术问题所采用的技术方案为：一种粘性阻尼振动信号中的模态参数提取方法，其特征在于包括以下步骤：

①设采样间隔为T_s秒，采样点数为2N，对连续的粘性阻尼振动信号进行奈奎斯特均匀采样，得到包含2N个采样值的采样信号，记为x，将x以向量形式表示为：x＝(x₁,x₂,…,x_2N-1,x_2N)^T，其中，T_s>2f，f表示粘性阻尼振动信号的最大频率的数值，N远大于目标系统的阶次，x₁、x₂、x_2N-1和x_2N对应表示x中的第1个采样值、第2个采样值、第2N-1个采样值和第2N个采样值，(x₁,x₂,…,x_2N-1,x_2N)^T为(x₁,x₂,…,x_2N-1,x_2N)的转置向量；

②利用x中的所有采样值构建一个自回归矩阵方程：Ψβ＝h，其中，Ψ表示维数为N×N的自回归矩阵，h＝(x_N+1,x_N+2,…,x_2N)^T，x₃、x_N、x_N+1和x_N+2对应表示x中的第3个采样值、第N个采样值、第N+1个采样值和第N+2个采样值，(x_N+1,x_N+2,…,x_2N)^T为(x_N+1,x_N+2,…,x_2N)的转置向量，β表示自回归系数向量；

然后利用最小二乘法计算β的最小二乘解，记为 其中，Ψ^T为Ψ的转置矩阵，(Ψ^T×Ψ)^-1为(Ψ^T×Ψ)的逆矩阵；

再将表示为其中，β₀、β₁、β_2N-2和β_2N-1对应表示中的第1个元素、第2个元素、第2N-1个元素和第2N个元素，(β₀,β₁,…,β_2N-2,β_2N-1)^T为(β₀,β₁,…,β_2N-2,β_2N-1)的转置向量；

③利用中的所有元素作为系数构造一个Prony多项式：V^2N+β_2N-1V^2N-1+β_2N-2V^2N-2+…+β₁V¹+β₀V⁰＝0，然后解该Prony多项式得到V的2N个根，由2N个根组成V的解向量，再对解向量进行修正，最后根据修正后的解向量计算粘性阻尼振动信号中的含有虚假模态的固有频率向量和含有虚假模态的固有阻尼比向量，对应记为ω和ξ，其中，V表示Prony多项式中的变量，V⁰、V¹、V^2N-2、V^2N-1和V^2N对应表示V的0次方、1次方、2N-2次方、2N-1次方和2N次方，ω和ξ的维数均为m×1维，1≤m≤N；

④首先以修正后的解向量为基础构造一个矩阵Φ，然后对x在Φ上的投影方程Φθ＝x进行多次稀疏求解，得到θ的多个稀疏的解向量，再通过优化并统计稀疏的解向量获取一个最终的稀疏向量最后根据稀疏向量剔除ω和ξ中的虚假模态，得到粘性阻尼振动信号中的固有频率、固有阻尼比，而中的所有非零元素为粘性阻尼振动信号中的振型系数，其中，θ表示维数为m×1的投影系数向量，1≤m≤N。

所述的步骤③的具体过程为：

③-1、利用中的所有元素作为系数构造一个Prony多项式：V^2N+β_2N-1V^2N-1+β_2N-2V^2N-2+…+β₁V¹+β₀V⁰＝0，然后解该Prony多项式得到V的2N个根，由2N个根组成V的解向量，记为 其中，V表示Prony多项式中的变量，V⁰、V¹、V^2N-2、V^2N-1和V^2N对应表示V的0次方、1次方、2N-2次方、2N-1次方和2N次方，V₁、V₂和V_2N对应表示V的第1个根、第2个根和第2N个根；

③-2、剔除中2范数值不在区间(0,1)内的所有元素及在中没有自身的共轭元素的所有元素，假设共留下m对共轭元素，则再剔除每对共轭元素中的一个元素，留下m个非共轭元素，实现对的修正，得到修正后的解向量，记为 其中，1≤m≤N，和对应表示中的第1个元素、第2个元素和第m个元素；

③-3、令得到R＝(R₁,R₂,…,R_m)，其中，ln()为自然对数函数，R₁、R₂和R_m对应表示R中的第1个元素、第2个元素和第m个元素，1≤r≤m，表示中的第r个元素；

③-4、根据R计算粘性阻尼振动信号中的含有虚假模态的固有频率向量和含有虚假模态的固有阻尼比向量，对应记为ω和ξ，ξ＝(ξ₁，ξ₂，…，ξ_m)^T，其中，符号“||”为取绝对值符号，为的转置向量，(ξ₁，ξ₂，…，ξ_m)^T为(ξ₁，ξ₂，…，ξ_m)的转置向量，1≤r≤m，R_r表示R中的第r个元素，Im(R_r)表示求取R_r的虚部，Re(R_r)表示求取R_r的实部。

所述的步骤④的具体过程为：

④-1、以为基础构造一个维数为2N×m的矩阵，记为Φ，然后将x在Φ上的投影方程表示为：Φθ＝x，其中，θ表示维数为m×1的投影系数向量；

④-2、对Ψ进行奇异值分解，得到Ψ的所有奇异值，然后计算Ψ的所有奇异值的平均值，再统计值大于平均值的奇异值的数量，记为P，最后根据P的值估计得到目标系统的阶次为其中，P≥2；

④-3、构造一个维数为2P×2N～4P×2N的高斯随机测量矩阵，记为G，然后对G中的每个列向量进行归一化操作，得到归一化后的矩阵，记为

④-4、利用对Φθ＝x进行测量，得到然后采用基于内点法的l₁优化算法函数l1_ls对进行计算，得到θ的解向量；

④-5、重复执行步骤④-3至步骤④-4，共执行K次，得到K个θ的解向量，将第k个θ的解向量记为 其中，1≤k≤K，K≥1，θ_1,k、θ_2,k和θ_m,k对应表示中的第1个元素、第2个元素和第m个元素，(θ_1,k,θ_2,k,…,θ_m,k)^T为(θ_1,k,θ_2,k,…,θ_m,k)的转置向量；

④-6、根据K个θ的解向量构造一个维数为m×K的矩阵，记为Γ，然后将Γ中绝对值小于设定的阈值的所有元素的值均置为零，得到一个新的矩阵，记为再统计中的每行中非零元素的数量，得到维数为m×1的统计结果向量，记为Q，其中，θ_1,1、θ_2,1和θ_m,1对应表示第1个θ的解向量中的第1个元素、第2个元素和第m个元素，θ_1,2、θ_2,2和θ_m,2对应表示第2个θ的解向量中的第1个元素、第2个元素和第m个元素，θ_1,K、θ_2,K和θ_m,K对应表示第K个θ的解向量中的第1个元素、第2个元素和第m个元素；

④-7、计算Q中值最大的非零元素与值最小的非零元素之间的相对误差，即最大值减去最小值的差除以最大值，如果相对误差大于设定的误差阈值，则用k-means聚类方法对Q中的所有非零元素进行聚类，得到两类，然后计算每类中的所有元素的平均值，接着提取出平均值小的一类中的每个元素在Q中的索引，再根据提取出的所有索引将中对应的行中的所有元素置零，得到一个新的矩阵，记为最后计算中的每行中的非零元素的平均值，由m个平均值构成稀疏向量，记为 如果相对误差小于或等于设定的误差阈值，则直接计算中的每行中的非零元素的平均值，由m个平均值构成稀疏向量，记为 其中，上述如果或中的一行中的所有元素均为零元素，则认为该行中的非零元素的平均值为0，和对应表示或中的第1行中的非零元素的平均值、第2行中的非零元素的平均值和第m行中的非零元素的平均值；

④-8、根据中的所有非零元素的索引，从ω中选出对应的所有元素，这些元素为粘性阻尼振动信号中的固有频率，并从ξ中选出对应的所有元素，这些元素为粘性阻尼振动信号中的固有阻尼比，而中的所有非零元素为粘性阻尼振动信号中的振型系数。

所述的步骤④-6中设定的阈值为0.05。

所述的步骤④-7中设定的误差阈值为20％。

在执行完所述的步骤①后执行所述的步骤②之前，对x进行预处理，预处理的过程为：对x进行去均值处理，以剔除x中的直流分量，使去均值处理后得到的信号的均值为0，将去均值处理后得到的信号记为x₀；然后令x＝x₀，其中，x＝x₀中的“＝”为赋值符号。

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

1)本发明方法利用采样信号中的所有采样值构建一个高阶次的自回归矩阵方程，由于自回归矩阵方程的阶次远高于实际值，因此获得的模态数量多于实际的模态数量，而噪声模态可以作为虚假模态被剔除掉，即在优化求解时给噪声提供了出口，可以更好地降低噪声影响，有效保留信号的原始信息，这样即使事先对振动信号进行消除噪声处理的效果不好，利用本发明方法也能提取得到精确稳定的模态参数，即本发明方法的抗噪声能力强，且提取出的模态参数的精确度高、稳定性好。

2)本发明方法通过解投影系数向量的最稀疏解来选择固有频率和固有阻尼比，避免了人工取阈值来选择模态的弊端，没有受到外界因素的影响，选择出的模态具有更高的可信度，降低了结果中包含虚假模态或者遗漏系统模态这两种情况出现的可能性。

3)本发明方法通过求的稀疏解来确定振型系数，所使用的稀疏求解方法属于搜索方法，当矩阵病态时，用其他算法求解，例如最小二乘算法，其结果会因为矩阵的病态性而误差较大，而且会因为解不稀疏而无法确定正确的振型系数，而搜索算法则受病态矩阵的影响很小，解中非零元素即为振型系数，因此稀疏求解算法的使用使本发明方法求解得到的振型系数具有高准确性和高稳定性的特点。

4)本发明方法将虚假模态剔除问题转化为一个稀疏优化问题，将选择真实模态的过程当作求最稀疏解的过程，在构造高斯随机测量矩阵时其稀疏度的确定不需要很精确，即高斯随机测量矩阵的行数可以取2P～4P行，具有较大的灵活性，使得本发明方法的实施更加简单，即使确定目标系统的阶次时有误差对结果影响也很小。

附图说明

图1为本发明方法的流程框图；

图2为对理论信号采样后得到的采样信号的波形图；

图3为理论信号与利用本发明方法提取得到的模态参数重构的信号的对比示意图。

具体实施方式

以下结合附图实施例对本发明作进一步详细描述。

实施例一：

本实施例提出的一种粘性阻尼振动信号中的模态参数提取方法，其在复指数法的基础上引入稀疏优化算法来消除虚假模态，其处理过程为：对连续的粘性阻尼振动信号进行奈奎斯特均匀采样得到采样信号，利用采样信号中的采样值构建一个高阶次的自回归矩阵，求出自回归系数；根据复指数法的方法理论，在自回归系数的基础上求出包含虚假模态的固有频率向量和固有阻尼比向量及信号在新基上的投影方程，其投影系数即为振型系数；对信号在新基上的投影方程进行稀疏优化求解，得到振型系数的稀疏解，根据稀疏解中非零元素筛选出模态参数。本发明方法的流程框图如图1所示，其包括以下步骤：

①设采样间隔为T_s秒，采样点数为2N，对连续的粘性阻尼振动信号进行奈奎斯特均匀采样，得到包含2N个采样值的采样信号，记为x，该信号x反映了目标系统(如桥梁等)的系统结构特性，这些系统结构特性综合映射到了模态空间中被称为模态参数，通过求模态参数最终可以得到目标系统的固有特性，将x以向量形式表示为：x＝(x₁,x₂,…,x_2N-1,x_2N)^T，其中，T_s>2f，f表示粘性阻尼振动信号的最大频率的数值，f的值可通过现有技术估计得到，要求N远大于目标系统的阶次，这为后续进行稀疏求解提供了可能，目标系统比如为桥梁等，在具体实施时在区间[50,400]内取一个正整数作为N的具体值，一般情况下N取值越大则计算复杂度越高，x₁、x₂、x_2N-1和x_2N对应表示x中的第1个采样值、第2个采样值、第2N-1个采样值和第2N个采样值，(x₁,x₂,…,x_2N-1,x_2N)^T为(x₁,x₂,…,x_2N-1,x_2N)的转置向量。

实际上，奈奎斯特采样的采样频率与采样点数应根据具体环境确定，例如大桥斜拉索的基频大约在1Hz左右，以10Hz为采样频率，采样500个点就足以满足需求，过高的采样频率(即过短的采样间隔)反而会影响模态参数的提取效果；另外，对连续的粘性阻尼振动信号进行奈奎斯特采样应该在目标系统稳定时进行，这样可获取有效的采样信号。

②利用x中的所有采样值构建一个自回归矩阵方程：Ψβ＝h，其中，Ψ表示维数为N×N的自回归矩阵，h＝(x_N+1,x_N+2,…,x_2N)^T，x₃、x_N、x_N+1和x_N+2对应表示x中的第3个采样值、第N个采样值、第N+1个采样值和第N+2个采样值，(x_N+1,x_N+2,…,x_2N)^T为(x_N+1,x_N+2,…,x_2N)的转置向量，β表示自回归系数向量。在实际处理过程中，为了避免自回归矩阵Ψ的病态性，要求N的取值为目标系统的阶次的30～50倍最好，即：如果目标系统的阶次为3，则N的取值最好为90～150。

然后利用最小二乘法计算β的最小二乘解，记为 其中，Ψ^T为Ψ的转置矩阵，(Ψ^T×Ψ)^-1为(Ψ^T×Ψ)的逆矩阵。

再将表示为其中，β₀、β₁、β_2N-2和β_2N-1对应表示中的第1个元素、第2个元素、第2N-1个元素和第2N个元素，(β₀,β₁,…,β_2N-2,β_2N-1)^T为(β₀,β₁,…,β_2N-2,β_2N-1)的转置向量。

③利用中的所有元素作为系数构造一个Prony多项式：V^2N+β_2N-1V^2N-1+β_2N-2V^2N-2+…+β₁V¹+β₀V⁰＝0，然后解该Prony多项式得到V的2N个根，由2N个根组成V的解向量，再对解向量进行修正，最后根据修正后的解向量计算粘性阻尼振动信号中的含有虚假模态的固有频率向量和含有虚假模态的固有阻尼比向量，对应记为ω和ξ，其中，V表示Prony多项式中的变量，V⁰、V¹、V^2N-2、V^2N-1和V^2N对应表示V的0次方、1次方、2N-2次方、2N-1次方和2N次方，ω和ξ的维数均为m×1维，1≤m≤N。

在此具体实施例中，步骤③的具体过程为：

③-1、利用中的所有元素作为系数构造一个Prony多项式：V^2N+β_2N-1V^2N-1+β_2N-2V^2N-2+…+β₁V¹+β₀V⁰＝0，然后解该Prony多项式得到V的2N个根，由2N个根组成V的解向量，记为 其中，V表示Prony多项式中的变量，V⁰、V¹、V^2N-2、V^2N-1和V^2N对应表示V的0次方、1次方、2N-2次方、2N-1次方和2N次方，V₁、V₂和V_2N对应表示V的第1个根、第2个根和第2N个根。

③-2、为了缓解后续剔除虚假模态的压力，此处对明显不符合条件的根进行剔除以减少虚假模态的产生，具体过程为：剔除中2范数值不在区间(0,1)内的所有元素及在中没有自身的共轭元素的所有元素，假设共留下m对共轭元素，则再剔除每对共轭元素中的一个元素，留下m个非共轭元素，实现对的修正，得到修正后的解向量，记为 其中，1≤m≤N，和对应表示中的第1个元素、第2个元素和第m个元素。

③-3、令得到R＝(R₁,R₂,…,R_m)，其中，ln()为自然对数函数，R₁、R₂和R_m对应表示R中的第1个元素、第2个元素和第m个元素，1≤r≤m，表示中的第r个元素。

③-4、根据R计算粘性阻尼振动信号中的含有虚假模态的固有频率向量和含有虚假模态的固有阻尼比向量，对应记为ω和ξ，ξ＝(ξ₁，ξ₂，…，ξ_m)^T，其中，符号“||”为取绝对值符号，为的转置向量，(ξ₁，ξ₂，…，ξ_m)^T为(ξ₁，ξ₂，…，ξ_m)的转置向量，1≤r≤m，R_r表示R中的第r个元素，Im(R_r)表示求取R_r的虚部，Re(R_r)表示求取R_r的实部。由于步骤②中构造的自回归矩阵Ψ使用的阶次远高于目标系统的阶次，因此计算得到的固有频率向量和固有阻尼比向量中包含了大量的虚假模态。

④首先以修正后的解向量为基础构造一个矩阵Φ，然后对x在Φ上的投影方程Φθ＝x进行多次稀疏求解，得到θ的多个稀疏的解向量，再通过优化并统计稀疏的解向量获取一个最终的稀疏向量最后根据稀疏向量剔除ω和ξ中的虚假模态，得到粘性阻尼振动信号中的固有频率、固有阻尼比，而中的所有非零元素为粘性阻尼振动信号中的振型系数，其中，θ表示维数为m×1的投影系数向量，1≤m≤N。

在此具体实施例中，步骤④的具体过程为：

④-1、以为基础构造一个维数为2N×m的矩阵，记为Φ，然后将x在Φ上的投影方程表示为：Φθ＝x，其中，θ表示维数为m×1的投影系数向量。由于矩阵Φ是通过中的元素按照指数方式构造得到的，这对之间的差别进行了放大，因此当对投影方程进行稀疏求解时便不会因为中的两个元素接近而选择Φ中错误的列，最终得到的稀疏的投影系数向量解决了振型系数选择的问题，避免了人为筛选带来的误差。

④-2、对Ψ进行奇异值分解，得到Ψ的所有奇异值，然后计算Ψ的所有奇异值的平均值，再统计值大于平均值的奇异值的数量，记为P，最后根据P的值估计得到目标系统的阶次为其中，P≥2。

在实际实施时，也可进行以下操作获得目标系统的阶次，具体过程为：对Ψ进行奇异值分解，得到Ψ＝SVD，V为对角矩阵，其对角线元素为奇异值，将绝对值小于小正数阈值(如0.0558)的奇异值置为零，得到V^*；然后得到Ψ^*＝SV^*D，对Ψ^*的反对角线元素用其平均值代替，得到新的矩阵之后对继续进行奇异值分解，重复操作，直到奇异值稳定，得到非零的奇异值的数量，记为P，最后根据P的值估计得到目标系统的阶次为其中，P≥2。

④-3、构造一个维数为2P×2N～4P×2N的高斯随机测量矩阵，记为G，然后对G中的每个列向量进行归一化操作，得到归一化后的矩阵，记为在此，高斯随机测量矩阵的维数可以为2P×2N，也可以为4P×2N，也可以为3P×2N等。

④-4、利用对Φθ＝x进行测量，得到然后采用基于内点法的l₁优化算法函数l1_ls对进行计算，得到θ的稀疏解向量。

④-5、重复执行步骤④-3至步骤④-4，共执行K次，得到K个θ的解向量，将第k个θ的解向量记为 其中，1≤k≤K，K≥1，重复执行次数越多，则最终提取的模态参数的精确度越高，但计算复杂度也会相应增加，因此折衷考虑提取精度与计算复杂度，在具体实施时可取K＝10，θ_1,k、θ_2,k和θ_m,k对应表示中的第1个元素、第2个元素和第m个元素，(θ_1,k,θ_2,k,…,θ_m,k)^T为(θ_1,k,θ_2,k,…,θ_m,k)的转置向量。

④-6、根据K个θ的解向量构造一个维数为m×K的矩阵，记为Γ，然后将Γ中绝对值小于设定的阈值的所有元素的值均置为零，得到一个新的矩阵，记为再统计中的每行中非零元素的数量，得到维数为m×1的统计结果向量，记为Q，其中，θ_1,1、θ_2,1和θ_m,1对应表示第1个θ的解向量中的第1个元素、第2个元素和第m个元素，θ_1,2、θ_2,2和θ_m,2对应表示第2个θ的解向量中的第1个元素、第2个元素和第m个元素，θ_1,K、θ_2,K和θ_m,K对应表示第K个θ的解向量中的第1个元素、第2个元素和第m个元素，在本实施例中取设定的阈值为0.05。

④-7、计算Q中值最大的非零元素与值最小的非零元素之间的相对误差，即Q中值最大的非零元素减去最小的非零元素的差除以最大的非零元素，如果相对误差大于设定的误差阈值，则用k-means聚类方法对Q中的所有非零元素进行聚类，得到两类，然后计算每类中的所有元素的平均值，接着提取出平均值小的一类中的每个元素在Q中的索引，再根据提取出的所有索引将中对应的行中的所有元素置零，得到一个新的矩阵，记为最后计算中的每行中的非零元素的平均值，由m个平均值构成稀疏向量，记为 如果相对误差小于或等于设定的误差阈值，则直接计算中的每行中的非零元素的平均值，由m个平均值构成稀疏向量，记为 其中，上述如果或中的一行中的所有元素均为零元素，则认为该行中的非零元素的平均值为0，和对应表示或中的第1行中的非零元素的平均值、第2行中的非零元素的平均值和第m行中的非零元素的平均值，在本实施例中取设定的误差阈值为20％。

④-8、根据中的所有非零元素的索引，从ω中选出对应的所有元素，这些元素为粘性阻尼振动信号中的固有频率，并从ξ中选出对应的所有元素，这些元素为粘性阻尼振动信号中的固有阻尼比，而中的所有非零元素为粘性阻尼振动信号中的振型系数。

实施例二：

本实施例提出的模态参数提取方法的步骤①至步骤④的过程与实施例一提出的模态参数提取方法的步骤①至步骤④的过程相同，不同之处仅在于在执行完步骤①后执行步骤②之前，对x进行预处理，预处理的过程为：对x进行去均值处理，以剔除x中的直流分量，即用x的各采样值减去x的均值，使去均值处理后得到的信号的均值为0，将去均值处理后得到的信号记为x₀；然后令x＝x₀，其中，x＝x₀中的“＝”为赋值符号。

为了更有效地说明本发明方法的有效性和可行性，进行实验验证。

假设存在一条无噪声的连续的粘性阻尼振动信号，其表达式为：其中，y₀表示连续的粘性阻尼振动信号，cos()为求余弦函数，t表示时间变量。忽略y₀的初始相位的影响，可以推导得到y₀的固有频率为：0.84、0.63、1.06，固有阻尼比为：0.012、0.09、0.078，振型系数为：0.6、1.1、-1.4。

(1)以采样间隔为T_s＝0.2秒均匀采样2N＝200个点，得到离散信号，对该离散信号y₀添加高斯随机白噪声得到采样信号x，其信噪比为8db。采样信号x的波形如图2所示。

(2)构建一个维数为100×100的自回归矩阵Ψ，Ψ中的第一行元素为x中的第1个至第100个采样值，Ψ中的第二行元素为x中的第2个至第101个采样值，依次类推，直至Ψ中的第100行，该行元素为x中的第100个至第199个采样值，即为：那么可以构造出自回归矩阵方程：Ψβ＝h，其中h＝(x₁₀₁,x₁₀₂,…,x₂₀₀)^T，利用最小二乘法可以得到

(3)求解多项式V²⁰⁰+β₁₉₉V¹⁹⁹+β₁₉₈V¹⁹⁸+…+β₁V¹+β₀V⁰＝0的根，得到 剔除中2范数值不在区间(0,1)内的所有元素及在中没有自身的共轭元素的所有元素，剩余23对共轭元素，取其中23个非共轭的元素，得到修正后的解向量 根据模态分析理论计算得到含有虚假模态的固有频率向量ω＝(2.0876,0.7498,…,2.3208)^T，含有虚假模态的固有阻尼比向量ξ＝(0.0001，0.0060，…，0.1471)^T，固有频率向量与固有阻尼比向量的长度均为23，包含虚假模态。

(4)以为基础构造一个维数为200×23的矩阵Φ，其第一行为第二行为以此类推，第200行为即然后将x在Φ上的投影方程表示为：Φθ＝x，其中，θ表示维数为23×1的投影系数向量；对Ψ进行奇异值分解，得到Ψ的所有奇异值，然后计算Ψ的所有奇异值的平均值，再统计值大于平均值的奇异值的数量P＝6，得到目标系统的阶次为构造一个维数为12×200的高斯随机测量矩阵G，并对矩阵G中的每一列列向量进行归一化，得到归一化后的矩阵用矩阵对投影方程Φθ＝x进行测量，得到然后采用基于内点法的l₁优化算法函数l1_ls计算得到θ的解向量；重复构造高斯随机测量矩阵并计算θ的解向量10次，得到然后将Γ中绝对值小于设定的阈值0.05的所有元素的值均置为零，得到一个新的矩阵再统计中的每行中非零元素的数量，得到维数为23×1的统计结果向量Q，Q＝(0,0,8,5,0,0,0,4,0,4,0,0,4,0,0,4,0,8,8,0,0,0,0)^T；计算Q中值最大的非零元素与值最小的非零元素之间的相对误差为50％，大于20％，用k-means聚类方法对Q中的所有非零元素进行聚类，得到两类，第一类：8,8,8，第二类：5,4,4,4,4，因为第二类的平均值小于第一类的平均值，提取出第二类中的每个元素在Q中的索引，再根据提取出的所有索引将中对应的行中的所有元素置零，即将中的第4、8、10、13、16行中的所有元素置零，此时得到的新的矩阵中仅第3、18、19行包含非零元素，最后计算中的每行中的非零元素的平均值，由m个平均值构成稀疏向量，记为 第3、18、19行中的非零元素的平均值分别为-1.1242,0.9536,0.4852；中的非零元素的索引分别为3、18、19，从ω中选出第3、18、19个元素，得到粘性阻尼振动信号中的固有频率为1.0540、0.6276、0.8398，从ξ中选出第3、18、19个元素，得到粘性阻尼振动信号中的固有阻尼比为0.0687、0.0906、0.0167，中的第3、18、19个元素为0.6、1.1、-1.4为粘性阻尼振动信号中的振型系数。

比较利用本发明方法提取得到的模态参数(固有频率、固有阻尼比和振型系数)与理论的模态参数，可以看出利用本发明方法在信噪比为8db的情况下提取得到的模态参数仍然具有较高的精确度。图3给出了理论信号y₀与利用本发明方法提取得到的模态参数重构的信号的对比，从图3中可以直观地表现出本发明方法提取的模态参数的精度。

Claims (6)

1.一种粘性阻尼振动信号中的模态参数提取方法，其特征在于包括以下步骤：

①设采样间隔为T_s秒，采样点数为2N，对连续的粘性阻尼振动信号进行奈奎斯特均匀采样，得到包含2N个采样值的采样信号，记为x，将x以向量形式表示为：x＝(x₁,x₂,…,x_2N-1,x_2N)^T，其中，T_s>2f，f表示粘性阻尼振动信号的最大频率的数值，N远大于目标系统的阶次，x₁、x₂、x_2N-1和x_2N对应表示x中的第1个采样值、第2个采样值、第2N-1个采样值和第2N个采样值，(x₁,x₂,…,x_2N-1,x_2N)^T为(x₁,x₂,…,x_2N-1,x_2N)的转置向量；

②利用x中的所有采样值构建一个自回归矩阵方程：Ψβ＝h，其中，Ψ表示维数为N×N的自回归矩阵， $\mathrm{\Ψ}=\left[\begin{array}{cccc}{x}_1& x_2& ...& x_N\\ x_2& x_3& ...& x_{N+1}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ x_N& x_{N+1}& ...& x_{2N-1}\end{array}\right],$ h＝(x_N+1,x_N+2,…,x_2N)^T，x₃、x_N、x_N+1和x_N+2对应表示x中的第3个采样值、第N个采样值、第N+1个采样值和第N+2个采样值，(x_N+1,x_N+2,…,x_2N)^T为(x_N+1,x_N+2,…,x_2N)的转置向量，β表示自回归系数向量；

然后利用最小二乘法计算β的最小二乘解，记为 其中，Ψ^T为Ψ的转置矩阵，(Ψ^T×Ψ)^-1为(Ψ^T×Ψ)的逆矩阵；

再将表示为其中，β₀、β₁、β_2N-2和β_2N-1对应表示中的第1个元素、第2个元素、第2N-1个元素和第2N个元素，(β₀,β₁,…,β_2N-2,β_2N-1)^T为(β₀,β₁,…,β_2N-2,β_2N-1)的转置向量；

③利用中的所有元素作为系数构造一个Prony多项式：V^2N+β_2N-1V^2N-1+β_2N-2V^2N-2+…+β₁V¹+β₀V⁰＝0，然后解该Prony多项式得到V的2N个根，由2N个根组成V的解向量，再对解向量进行修正，最后根据修正后的解向量计算粘性阻尼振动信号中的含有虚假模态的固有频率向量和含有虚假模态的固有阻尼比向量，对应记为ω和ξ，其中，V表示Prony多项式中的变量，V⁰、V¹、V^2N-2、V^2N-1和V^2N对应表示V的0次方、1次方、2N-2次方、2N-1次方和2N次方，ω和ξ的维数均为m×1维，1≤m≤N；

④首先以修正后的解向量为基础构造一个矩阵Φ，然后对x在Φ上的投影方程Φθ＝x进行多次稀疏求解，得到θ的多个稀疏的解向量，再通过优化并统计稀疏的解向量获取一个最终的稀疏向量最后根据稀疏向量剔除ω和ξ中的虚假模态，得到粘性阻尼振动信号中的固有频率、固有阻尼比，而中的所有非零元素为粘性阻尼振动信号中的振型系数，其中，θ表示维数为m×1的投影系数向量，1≤m≤N。

2.根据权利要求1所述的一种粘性阻尼振动信号中的模态参数提取方法，其特征在于所述的步骤③的具体过程为：

③-1、利用中的所有元素作为系数构造一个Prony多项式：V^2N+β_2N-1V^2N-1+β_2N-2V^2N-2+…+β₁V¹+β₀V⁰＝0，然后解该Prony多项式得到V的2N个根，由2N个根组成V的解向量，记为 其中，V表示Prony多项式中的变量，V⁰、V¹、V^2N-2、V^2N-1和V^2N对应表示V的0次方、1次方、2N-2次方、2N-1次方和2N次方，V₁、V₂和V_2N对应表示V的第1个根、第2个根和第2N个根；

③-2、剔除中2范数值不在区间(0,1)内的所有元素及在中没有自身的共轭元素的所有元素，假设共留下m对共轭元素，则再剔除每对共轭元素中的一个元素，留下m个非共轭元素，实现对的修正，得到修正后的解向量，记为 其中，1≤m≤N，和对应表示中的第1个元素、第2个元素和第m个元素；

③-3、令得到R＝(R₁,R₂,…,R_m)，其中，ln()为自然对数函数，R₁、R₂和R_m对应表示R中的第1个元素、第2个元素和第m个元素，1≤r≤m，表示中的第r个元素；

③-4、根据R计算粘性阻尼振动信号中的含有虚假模态的固有频率向量和含有虚假模态的固有阻尼比向量，对应记为ω和ξ，ξ＝(ξ₁，ξ₂，…，ξ_m)^T，其中，符号“||”为取绝对值符号，为的转置向量，(ξ₁，ξ₂，…，ξ_m)^T为(ξ₁，ξ₂，…，ξ_m)的转置向量，1≤r≤m，R_r表示R中的第r个元素，Im(R_r)表示求取R_r的虚部，Re(R_r)表示求取R_r的实部。

3.根据权利要求2所述的一种粘性阻尼振动信号中的模态参数提取方法，其特征在于所述的步骤④的具体过程为：

④-1、以为基础构造一个维数为2N×m的矩阵，记为Φ，然后将x在Φ上的投影方程表示为：Φθ＝x，其中，θ表示维数为m×1的投影系数向量；

④-2、对Ψ进行奇异值分解，得到Ψ的所有奇异值，然后计算Ψ的所有奇异值的平均值，再统计值大于平均值的奇异值的数量，记为P，最后根据P的值估计得到目标系统的阶次为其中，P≥2；

④-3、构造一个维数为2P×2N～4P×2N的高斯随机测量矩阵，记为G，然后对G中的每个列向量进行归一化操作，得到归一化后的矩阵，记为

④-4、利用对Φθ＝x进行测量，得到然后采用基于内点法的l₁优化算法函数l1_ls对进行计算，得到θ的解向量；

④-5、重复执行步骤④-3至步骤④-4，共执行K次，得到K个θ的解向量，将第k个θ的解向量记为 其中，1≤k≤K，K≥1，θ_1,k、θ_2,k和θ_m,k对应表示中的第1个元素、第2个元素和第m个元素，(θ_1,k,θ_2,k,…,θ_m,k)^T为(θ_1,k,θ_2,k,…,θ_m,k)的转置向量；

④-6、根据K个θ的解向量构造一个维数为m×K的矩阵，记为Γ， $\mathrm{\Γ}=\left[\begin{array}{cccccc}{\mathrm{\θ}}_{\mathrm{1,1}}& {\mathrm{\θ}}_{\mathrm{1,2}}& ...& {\mathrm{\θ}}_{1,k}& ...& {\mathrm{\θ}}_{1.K}\\ {\mathrm{\θ}}_{\mathrm{2,1}}& {\mathrm{\θ}}_{\mathrm{2,2}}& ...& {\mathrm{\θ}}_{2,k}& ...& {\mathrm{\θ}}_{2,K}\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ {\mathrm{\θ}}_{m,1}& {\mathrm{\θ}}_{m,2}& ...& {\mathrm{\θ}}_{m,k}& ...& {\mathrm{\θ}}_{m,K}\end{array}\right],$ 然后将Γ中绝对值小于设定的阈值的所有元素的值均置为零，得到一个新的矩阵，记为再统计中的每行中非零元素的数量，得到维数为m×1的统计结果向量，记为Q，其中，θ_1,1、θ_2,1和θ_m,1对应表示第1个θ的解向量中的第1个元素、第2个元素和第m个元素，θ_1,2、θ_2,2和θ_m,2对应表示第2个θ的解向量中的第1个元素、第2个元素和第m个元素，θ_1,K、θ_2,K和θ_m,K对应表示第K个θ的解向量中的第1个元素、第2个元素和第m个元素；

④-7、计算Q中值最大的非零元素与值最小的非零元素之间的相对误差，即最大值减去最小值的差除以最大值，如果相对误差大于设定的误差阈值，则用k-means聚类方法对Q中的所有非零元素进行聚类，得到两类，然后计算每类中的所有元素的平均值，接着提取出平均值小的一类中的每个元素在Q中的索引，再根据提取出的所有索引将中对应的行中的所有元素置零，得到一个新的矩阵，记为最后计算中的每行中的非零元素的平均值，由m个平均值构成稀疏向量，记为 如果相对误差小于或等于设定的误差阈值，则直接计算中的每行中的非零元素的平均值，由m个平均值构成稀疏向量，记为 其中，上述如果或中的一行中的所有元素均为零元素，则认为该行中的非零元素的平均值为0，和对应表示或中的第1行中的非零元素的平均值、第2行中的非零元素的平均值和第m行中的非零元素的平均值；

④-8、根据中的所有非零元素的索引，从ω中选出对应的所有元素，这些元素为粘性阻尼振动信号中的固有频率，并从ξ中选出对应的所有元素，这些元素为粘性阻尼振动信号中的固有阻尼比，而中的所有非零元素为粘性阻尼振动信号中的振型系数。

4.根据权利要求3所述的一种粘性阻尼振动信号中的模态参数提取方法，其特征在于所述的步骤④-6中设定的阈值为0.05。

5.根据权利要求4所述的一种粘性阻尼振动信号中的模态参数提取方法，其特征在于所述的步骤④-7中设定的误差阈值为20％。

6.根据权利要求1至5中任一项所述的一种粘性阻尼振动信号中的模态参数提取方法，其特征在于在执行完所述的步骤①后执行所述的步骤②之前，对x进行预处理，预处理的过程为：对x进行去均值处理，以剔除x中的直流分量，使去均值处理后得到的信号的均值为0，将去均值处理后得到的信号记为x₀；然后令x＝x₀，其中，x＝x₀中的“＝”为赋值符号。