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

Abstract

本发明公开了一种粘性阻尼振动信号中的模态参数提取方法，其通过选用较大的系统阶次的初始值，利用稀疏优化OMP算法计算状态矩阵，之后对可观矩阵中的每一行都运用OMP算法计算输出矩阵，并对所有的输出矩阵求取均值，较大程度上减少了噪声对结果的影响，从而提高了本发明方法的消噪能力和识别精度；然后计算固有频率、固有阻尼比和固有模态振型系数，并运用K‑means算法从众多模态参数中选出有效模态参数，来消除虚假模态，从而大大消弱了系统阶次对模态参数提取精度的影响。

Description

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

技术领域

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

背景技术

随着超高层建筑和大跨桥梁等的兴建，健康监测系统越来越多地应用于此类大型土木工程结构中。这类大型土木工程结构的损伤诊断和受损结构评定、实时监测和安全预警、有限元模型修正等成为了健康监测系统的重点，而粘性阻尼振动信号中的模态参数提取是需要首先解决的问题。现有的模态参数提取方法有很多种，如应用广泛的有随机子空间法。

随机子空间法是一种基于白噪声假定的时域方法。随机子空间法最关键的问题就是需要精准的确定系统阶次，目前系统阶次确定方法主要有两种，分别是奇异值跳跃法和稳定图法，但是，一方面，这两种系统阶次确定方法都需要人工参与观察波形的变化，而人工参与必然会很大程度上影响确定的系统阶次的精确度，从而会导致提取得到的模态参数中存在虚假模态，提取的模态参数不精确；另一方面，随着科学技术的快速发展，模态参数提取除了需要满足高精度和鲁棒性等要求外，还对模态参数提取的自动化程度提出了更高的要求，而这种需要人工参与的系统阶次确定方法显然不能满足发展的需求。此外，随机子空间法受外界环境噪声影响大，也会导致提取得到的模态参数中存在虚假模态，提取的模态参数不精确。

目前，研究人员着力于研究如何提高系统阶次的精确度，来避免产生虚假模态，但通过提高系统阶次的精确度仍不能完全避免虚假模态的产生，因此，有必要研究一种能够消除虚假模态的模态参数提取方法。

发明内容

本发明所要解决的技术问题是提供一种粘性阻尼振动信号中的模态参数提取方法，其削弱了系统阶次对模态参数提取精度的影响，通过消除虚假模态有效地提高了提取的模态参数的精度。

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

①对连续的粘性阻尼振动信号进行奈奎斯特均匀采样，采样间隔为T_S秒，采样点数为2N点，得到包含2N个采样值的采样信号，记为x，将x以列向量形式表示为x＝(x₁,x₂,…,x_2N-1,x_2N)^T；其中，T_S的取值满足奈奎斯特采样定律，N的取值要求远大于目标系统的系统阶次的预估值，(x₁,x₂,…,x_2N-1,x_2N)^T为(x₁,x₂,…,x_2N-1,x_2N)的转置，x₁,x₂,…,x_2N-1,x_2N对应表示x中的第1个采样值、第2个采样值、…、第2N-1个采样值、第2N个采样值；

②利用x中的所有采样值构建一个Hankel矩阵，记为H，H中的每一条副对角线上的元素都相等，然后将H分成两个子矩阵，记为Y_p和Y_f，其中，H的维数为2i×j，Y_p和Y_f的维数均为i×j，i与j满足关系：2i+j-1＝2N，i的取值为人为确定的常数，且n<i<N，n表示目标系统的系统阶次的预估值，x₃、x_j、x_j+1、x_i、x_i+1、x_i+2、x_i+3、x_i+j-1、x_i+j、x_i+j+1、x_2i、x_2i+1、x_2i+j-1对应表示x中的第3个采样值、第j个采样值、第j+1个采样值、第i个采样值、第i+1个采样值、第i+2个采样值、第i+3个采样值、第i+j-1个采样值、第i+j个采样值、第i+j+1个采样值、第2i个采样值、第2i+1个采样值、第2i+j-1个采样值；

③利用Y_p和Y_f构建一个维数为i×i的Toeplitz矩阵，记为T，然后选取T的第1行至第i-1行构成维数为(i-1)×i的第一子矩阵，记为T₁；并选取T的第2行至第i行构成维数为(i-1)×i的第二子矩阵，记为T₂；其中，为Y_p的转置；

④根据状态空间方程的定义对T₁进行分解，得到T₁的可观矩阵和可控反转矩阵，对应记为Γ₁和Δ₁，T₁＝Γ₁Δ₁，Γ₁＝(C_,CA,…,CA^i-2)^T；并对T₁进行SVD分解，得到T₁的第一正交矩阵、第二正交矩阵和正奇异矩阵组成的对角矩阵，对应记为U₁、V₁和S₁，T₁＝U₁S₁V₁ ^T；然后根据T₁＝Γ₁Δ₁和T₁＝U₁S₁V₁ ^T，令

同样，根据状态空间方程的定义对T₂进行分解，得到T₂的可观矩阵和可控反转矩阵，对应记为Γ₂和Δ₂，T₂＝Γ₂Δ₂，Γ₂＝(CA,CA²,…,CA^i-1)^T；并对T₂进行SVD分解，得到T₂的第一正交矩阵、第二正交矩阵和正奇异矩阵组成的对角矩阵，对应记为U₂、V₂和S₂，然后根据T₂＝Γ₂Δ₂和令

其中，Γ₁的维数为(i-1)×1，Δ₁的维数为1×i，C表示维数为1×(i-1)的输出矩阵，A表示维数为(i-1)×(i-1)的状态矩阵，(C,CA,…,CA^i-2)^T为(C,CA,…,CA^i-2)的转置，A^{i -2}为A的i-2次方，V₁ ^T为V₁的转置，为S₁的次方，Γ₂的维数为(i-1)×1，Δ₂的维数为1×i，(CA,CA²,…,CA^i-1)^T为(CA,CA²,…,CA^i-1)的转置，A²为A的2次方，A^i-1为A的i-1次方，为V₂的转置，为S₂的次方；

⑤根据Γ₁＝(C_,CA,…,CA^i-2)^T和Γ₂＝(CA,CA²,…,CA^i-1)^T，确定Γ₁与Γ₂之间的关系为：Γ₂ ^T＝Γ₁ ^TΛ；然后将Γ₁ ^T分解为其中，Γ₂ ^T为Γ₂的转置，Γ₁ ^T为Γ₁的转置，Λ表示维数为(i-1)×(i-1)的对角矩阵，P₁表示维数为1×(i-1)的行向量，P₁＝(C,C,…,C)，Q₁表示维数为(i-1)×(i-1)的对角矩阵，I表示i-1阶的单位矩阵；

⑥根据稀疏优化的原理，确定稀疏度小于或等于i；并根据目标系统的系统阶次的预估值，确定稀疏度大于2n；然后确定稀疏度的取值范围为大于2n且小于或等于i；接着在稀疏度的取值范围内选择一个正整数k作为稀疏度的确定值，k∈(2n,i]；再将大于2n的值i赋值给目标系统的系统阶次作为目标系统的系统阶次的初始值；

⑦令q表示执行的次数，令Q表示重复执行的总次数；其中，q的初始值为1，Q≥2；

⑧利用第q次执行时随机产生的第一高斯随机测量矩阵H₁对Γ₂ ^T＝Γ₁ ^TΛ进行观测，构建得到第q次执行时的第一稀疏优化模型，描述为：并利用第q次执行时随机产生的第二高斯随机测量矩阵H₂对进行观测，构建得到第q次执行时的第二稀疏优化模型，描述为：其中，H₁的维数为M×1，H₂的维数为M×(i-1)，M∈[k-10,k+10]，min为取最小值函数，符号“|| ||₁”为求矩阵的1-范数符号，s.t.表示“受约束于…”，符号“|| ||”为求欧氏距离符号，P₁ ^T为P₁的转置，Q₁ ^T为Q₁的转置，σ₁和σ₂均为值很小的常数；

⑨根据目标系统的系统阶次的初始值和稀疏度的确定值，利用OMP方法对进行稀疏求解，得到对角矩阵Λ；然后根据和求解得到的Λ，计算得到A，并令A_q＝A；之后对A_q进行特征值分解，得到A_q的特征值向量，记为D_q；接着根据和得到的A，计算得到Q₁；同样，根据目标系统的系统阶次的初始值和稀疏度的确定值，利用OMP方法对进行稀疏求解，得到P₁ ^T；最后根据P₁＝(C,C,…,C)得到C，并令C_q＝C；其中，A_q和C_q的初始值均为0；

⑩根据A_q和C_q计算第q次执行得到的连续的粘性阻尼振动信号的固有频率向量、固有阻尼比向量和固有模态振型系数向量，对应记为和其中，和的维数均为1×(i-1)；

判断q≥Q是否成立，如果成立，则结束重复执行的过程，得到D₁,D₂,…,D_q,…,D_Q、C₁,C₂,…,C_q,…,C_Q、和然后执行步骤否则，令q＝q+1，然后返回步骤⑧继续执行；其中，q＝q+1中的“＝”为赋值符号；

对D₁,D₂,…,D_q,…,D_Q各自中的部分元素强制置零，再从D₁,D₂,…,D_q,…,D_Q中任意选择一个作为状态矩阵的最终特征值向量，记为D^*；同样，对C₁,C₂,…,C_q,…,C_Q各自中的部分元素强制置零，再从C₁,C₂,…,C_q,…,C_Q中任意选择一个作为输出矩阵的最终值，记为C^*；

并计算Q次执行共得到的Q个固有频率向量的均值向量Q个固有阻尼比向量的均值向量和Q个固有模态振型系数向量的均值向量

根据D^*中的每个零元素的索引，将中对应索引的模态参数确定为虚假模态；而根据D^*中的每个非零元素的索引，将中对应索引的模态参数确定为有效模态参数，并提取出；

根据D^*中的每个零元素的索引，将中对应索引的模态参数确定为虚假模态；而根据D^*中的每个非零元素的索引，将中对应索引的模态参数确定为有效模态参数，并提取出；

根据C^*中的每个零元素的索引，将中对应索引的模态参数确定为虚假模态；而根据C^*中的每个非零元素的索引，将中对应索引的模态参数确定为有效模态参数，并提取出。

所述的步骤⑧中的H₁和H₂中的每个元素独立的服从均值为0且方差为的高斯分布。

所述的步骤中的D^*和C^*的获取过程为：

_1、在步骤的基础上，计算D₁,D₂,…,D_q,…,D_Q各自中的所有元素的均值，将D_q中的所有元素的均值记为同样，计算C₁,C₂,…,C_q,…,C_Q各自中的所有元素的均值，将C_q中的所有元素的均值记为

_2、比较D₁,D₂,…,D_q,…,D_Q各自中的每个元素与对应的均值的大小，将D₁,D₂,…,D_q,…,D_Q各自中小于对应的均值的元素强制置零；对于D_q，比较D_q中的每个元素与的大小，将D_q中小于的元素强制置零；

同样，比较C₁,C₂,…,C_q,…,C_Q各自中的每个元素与对应的均值的大小，将C₁,C₂,…,C_q,…,C_Q各自中小于对应的均值的元素强制置零；对于C_q，比较C_q中的每个元素与的大小，将C_q中小于的元素强制置零；

_3、在步骤_2的基础上，统计D₁,D₂,…,D_q,…,D_Q中相同索引值的Q个元素中的非零元素的总个数，将D₁,D₂,…,D_q,…,D_Q中索引值为g的Q个元素中的非零元素的总个数记为d_Dg；然后将统计得到的所有数据按序构成一个集合，记为d_D，d_D＝{d_D1,d_D2,…,d_Dg,…,d_DG}；

同样，统计C₁,C₂,…,C_q,…,C_Q中相同索引值的Q个元素中的非零元素的总个数，将C₁,C₂,…,C_q,…,C_Q中索引值为g'的Q个元素中的非零元素的总个数记为d_Cg'；然后将统计得到的所有数据按序构成一个集合，记为d_C，d_C＝{d_C1,d_C2,…,d_Cg',…,d_CG'}；

其中，g的初始值为1，1≤g≤G，G表示D₁,D₂,…,D_q,…,D_Q各自中包含的元素的总个数，d_D1表示D₁,D₂,…,D_q,…,D_Q中索引值为1的Q个元素中的非零元素的总个数，d_D2表示D₁,D₂,…,D_q,…,D_Q中索引值为2的Q个元素中的非零元素的总个数，d_DG表示D₁,D₂,…,D_q,…,D_Q中索引值为G的Q个元素中的非零元素的总个数，g'的初始值为1，1≤g'≤G'，G'表示C₁,C₂,…,C_q,…,C_Q各自中包含的元素的总个数，d_C1表示C₁,C₂,…,C_q,…,C_Q中索引值为1的Q个元素中的非零元素的总个数，d_C2表示C₁,C₂,…,C_q,…,C_Q中索引值为2的Q个元素中的非零元素的总个数，d_CG'表示C₁,C₂,…,C_q,…,C_Q中索引值为G'的Q个元素中的非零元素的总个数；

_4、计算d_D中的所有元素的均值，记为然后计算d_D中的每个元素与的差值的绝对值，将d_Dg与的差值的绝对值记为再将d_D中的所有元素与的差值的绝对值按序构成一个集合，记为f_D，

同样，计算d_C中的所有元素的均值，记为然后计算d_C中的每个元素与的差值的绝对值，将d_Cg'与的差值的绝对值记为再将d_C中的所有元素与的差值的绝对值按序构成一个集合，记为f_C，

其中，表示d_D1与的差值的绝对值，表示d_D2与的差值的绝对值，表示d_DG与的差值的绝对值，表示d_C1与的差值的绝对值，表示d_C2与的差值的绝对值，表示d_CG'与的差值的绝对值；

_5、在f_D中任选两个元素作为聚类中心；然后根据所选的两个聚类中心，利用K-means算法将f_D中的所有元素分为两个类；接着计算f_D的两个类各自中的所有元素的均值；之后提取出小值均值对应的类中的每个元素在f_D中的索引值；再将提取出的所有索引值按序构成一个集合，记为h_D；

同样，在f_C中任选两个元素作为聚类中心；然后根据所选的两个聚类中心，利用K-means算法将f_C中的所有元素分为两个类；接着计算f_C的两个类各自中的所有元素的均值；之后提取出小值均值对应的类中的每个元素在f_C中的索引值；再将提取出的所有索引值按序构成一个集合，记为h_C；

_6、将经步骤_2处理后得到的D₁,D₂,…,D_q,…,D_Q各自中索引值与h_D中的每个索引值相同的元素强制置零；

同样，将经步骤_2处理后得到的C₁,C₂,…,C_q,…,C_Q各自中索引值与h_C中的每个索引值相同的元素强制置零；

_7、在步骤_6的基础上，从D₁,D₂,…,D_q,…,D_Q中任选一个，假设选择D_q，则令D^*＝D_q；同样，从C₁,C₂,…,C_q,…,C_Q中任选一个，假设选择C_q，则令C^*＝C_q。

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

本发明方法通过选用较大的系统阶次的初始值，利用稀疏优化OMP算法代替最小二乘法计算状态矩阵，之后对可观矩阵中的每一行都运用OMP算法计算输出矩阵，并对所有的输出矩阵求取均值，较大程度上减少了噪声对结果的影响，从而提高了本发明方法的消噪能力和识别精度；然后计算固有频率、固有阻尼比和固有模态振型系数，并运用K-means算法从众多模态参数中选出有效模态参数，来消除虚假模态，从而大大消弱了系统阶次对模态参数提取精度的影响。

附图说明

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

图2a为在采样频率为15Hz、采样点数为1000、系统阶次的初始值为80、稀疏度的确定值为24、给采样信号添加高斯白噪声的情况下，高斯白噪声的信噪比逐次递增，利用本发明方法提取的固有频率与理论值之间的相对误差及利用随机子空间算法提取的固有频率与理论值之间的相对误差的对比；

图2b为在采样频率为15Hz、采样点数为1000、系统阶次的初始值为80、稀疏度的确定值为24、给采样信号添加混合噪声(由脉冲噪声和高斯白噪声组成)，混合噪声的信噪比逐次递增，利用本发明方法提取的固有频率与理论值之间的相对误差及利用随机子空间算法提取的固有频率与理论值之间的相对误差的对比；

图2c为在采样频率为15Hz、采样点数为1000、系统阶次的初始值为80、稀疏度的确定值为24、给采样信号添加混合噪声(由脉冲噪声和高斯白噪声组成)，混合噪声的信噪比逐次递增，利用本发明方法提取的固有阻尼比与理论值之间的相对误差及利用随机子空间算法提取的固有阻尼比与理论值之间的相对误差的对比；

图2d为在采样频率为15Hz、采样点数为1000、系统阶次的初始值为80、稀疏度的确定值为24、给采样信号添加混合噪声(由脉冲噪声和高斯白噪声组成)，混合噪声的信噪比逐次递增，利用本发明方法提取的固有模态振型系数与理论值之间的相对误差及利用随机子空间算法提取的固有模态振型系数与理论值之间的相对误差的对比。

具体实施方式

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

本发明提出的一种粘性阻尼振动信号中的模态参数提取方法，其流程框图如图1所示，其包括以下步骤：

①对连续的粘性阻尼振动信号进行奈奎斯特均匀采样，采样间隔为T_S秒，采样点数为2N点，得到包含2N个采样值的采样信号，记为x，将x以列向量形式表示为x＝(x₁,x₂,…,x_2N-1,x_2N)^T，x反映了目标系统(如桥梁、汽车等)的系统结构特性，这些系统结构特性综合映射到了模态空间中被称为模态参数，通过求模态参数最终可以得到目标系统的固有特性；其中，T_S的取值满足奈奎斯特采样定律，即T_S>2f，f表示连续的粘性阻尼振动信号的最大频率，f的值可以通过现有技术估计得到，在本实施例中取N的取值要求远大于目标系统的系统阶次的预估值，为后续稀疏求解做铺垫，在本实施例中在区间[200,500]内取一个正整数作为N的具体值，如取N＝500，一般不建议N的取值大于500，因为N的取值越大计算复杂度越高，而区间[200,500]是通过多次实验确定的，目标系统的系统阶次n的预估值由人为预估得到，如目标系统的系统阶次的预估值为8，在本实施例中目标系统的系统阶次的预估值精确与否无关紧要，(x₁,x₂,…,x_2N-1,x_2N)^T为(x₁,x₂,…,x_2N-1,x_2N)的转置，x₁,x₂,…,x_2N-1,x_2N对应表示x中的第1个采样值、第2个采样值、…、第2N-1个采样值、第2N个采样值。

实际上，奈奎斯特均匀采样的采样间隔与采样点数应根据具体环境确定，例如大桥斜拉索的基频在1Hz左右，以15Hz为采样频率，如以秒为采样间隔，采样1000个点即N＝500就可以满足要求。

②利用x中的所有采样值构建一个Hankel矩阵，记为H，H中的每一条副对角线上的元素都相等，然后将H分成两个子矩阵，记为Y_p和Y_f，其中，H的维数为2i×j，Y_p和Y_f的维数均为i×j，i和j为两个重要的控制参数，i与j满足关系：2i+j-1＝2N，i的取值为人为确定的常数，且n<i<N，n表示目标系统的系统阶次的预估值，在本实施例中取i＝80、j＝841、n＝8，x₃、x_j、x_j+1、x_i、x_i+1、x_i+2、x_i+3、x_i+j-1、x_i+j、x_i+j+1、x_2i、x_2i+1、x_2i+j-1对应表示x中的第3个采样值、第j个采样值、第j+1个采样值、第i个采样值、第i+1个采样值、第i+2个采样值、第i+3个采样值、第i+j-1个采样值、第i+j个采样值、第i+j+1个采样值、第2i个采样值、第2i+1个采样值、第2i+j-1个采样值。

③利用Y_p和Y_f构建一个维数为i×i的Toeplitz矩阵，记为T，然后选取T的第1行至第i-1行构成维数为(i-1)×i的第一子矩阵，记为T₁；并选取T的第2行至第i行构成维数为(i-1)×i的第二子矩阵，记为T₂；其中，为Y_p的转置。

④根据状态空间方程的定义对T₁进行分解，得到T₁的可观矩阵和可控反转矩阵，对应记为Γ₁和Δ₁，T₁＝Γ₁Δ₁，Γ₁＝(C_,CA,…,CA^i-2)^T；并对T₁进行SVD分解，得到T₁的第一正交矩阵、第二正交矩阵和正奇异矩阵组成的对角矩阵，对应记为U₁、V₁和S₁，T₁＝U₁S₁V₁ ^T；然后根据T₁＝Γ₁Δ₁和T₁＝U₁S₁V₁ ^T，令使Γ₁已知。

同样，根据状态空间方程的定义对T₂进行分解，得到T₂的可观矩阵和可控反转矩阵，对应记为Γ₂和Δ₂，T₂＝Γ₂Δ₂，Γ₂＝(CA_,CA²,…,CA^i-1)^T；并对T₂进行SVD分解，得到T₂的第一正交矩阵、第二正交矩阵和正奇异矩阵组成的对角矩阵，对应记为U₂、V₂和S₂，然后根据T₂＝Γ₂Δ₂和令使Γ₂已知。

其中，Γ₁的维数为(i-1)×1，Δ₁的维数为1×i，C表示维数为1×(i-1)的输出矩阵，C未知，A表示维数为(i-1)×(i-1)的状态矩阵，A未知，(C,CA,…,CA^i-2)^T为(C,CA,…,CA^i-2)的转置，A^i-2为A的i-2次方，V₁ ^T为V₁的转置，U₁的奇异值大于S₁的奇异值，S₁的奇异值大于V₁的奇异值，为S₁的次方，Γ₂的维数为(i-1)×1，Δ₂的维数为1×i，(CA,CA²,…,CA^i-1)^T为(CA,CA²,…,CA^i-1)的转置，A²为A的2次方，A^i-1为A的i-1次方，为V₂的转置，U₂的奇异值大于S₂的奇异值，S₂的奇异值大于V₂的奇异值，为S₂的次方。

⑤根据Γ₁＝(C_,CA,…,CA^i-2)^T和Γ₂＝(CA_,CA²,…,CA^i-1)^T，确定Γ₁与Γ₂之间的关系为：Γ₂ ^T＝Γ₁ ^TΛ；然后将Γ₁ ^T分解为其中，Γ₂ ^T为Γ₂的转置，Γ₁ ^T为Γ₁的转置，Λ表示维数为(i-1)×(i-1)的对角矩阵，P₁表示维数为1×(i-1)的行向量，P₁＝(C,C,…,C)，Q₁表示维数为(i-1)×(i-1)的对角矩阵，I表示i-1阶的单位矩阵。

⑥根据稀疏优化的原理，确定稀疏度小于或等于i；并根据目标系统的系统阶次的预估值，确定稀疏度大于2n，使稀疏度大于n的2倍是为了使稀疏度更加精确；然后确定稀疏度的取值范围为大于2n且小于或等于i；接着在稀疏度的取值范围内选择一个正整数k作为稀疏度的确定值，k∈(2n,i]，在本实施例中取k＝24,i＝80；再将大于2n的值i赋值给目标系统的系统阶次作为目标系统的系统阶次的初始值。

⑦令q表示执行的次数，令Q表示重复执行的总次数；其中，q的初始值为1，Q≥2。

⑧利用第q次执行时随机产生的第一高斯随机测量矩阵H₁对Γ₂ ^T＝Γ₁ ^TΛ进行观测，构建得到第q次执行时的第一稀疏优化模型，描述为：并利用第q次执行时随机产生的第二高斯随机测量矩阵H₂对进行观测，构建得到第q次执行时的第二稀疏优化模型，描述为：其中，H₁的维数为M×1，H₂的维数为M×(i-1)，M∈[k-10,k+10]，在本实施例中取M＝k，min为取最小值函数，符号“|| ||₁”为求矩阵的1-范数符号，s.t.表示“受约束于…”，符号“|| ||”为求欧氏距离符号，P₁ ^T为P₁的转置，Q₁ ^T为Q₁的转置，σ₁和σ₂均为值很小的常数，在本实施例中取σ₁＝σ₂＝0.01。

在此具体实施例中，步骤⑧中的H₁和H₂中的每个元素独立的服从均值为0且方差为的高斯分布。

⑨根据目标系统的系统阶次的初始值和稀疏度的确定值，利用OMP(orthogonalmatching pursuit)方法对进行稀疏求解，得到对角矩阵Λ；然后根据和求解得到的Λ，计算得到A，并令A_q＝A；之后对A_q进行特征值分解，得到A_q的特征值向量，记为D_q；接着根据和得到的A，计算得到Q₁；同样，根据目标系统的系统阶次的初始值和稀疏度的确定值，利用OMP(orthogonalmatching pursuit)方法对进行稀疏求解，得到P₁ ^T；最后根据P₁＝(C,C,…,C)得到C，并令C_q＝C；其中，A_q和C_q的初始值均为0。

⑩根据A_q和C_q计算第q次执行得到的连续的粘性阻尼振动信号的固有频率向量、固有阻尼比向量和固有模态振型系数向量，对应记为和其中，和的维数均为1×(i-1)。

判断q≥Q是否成立，如果成立，则结束重复执行的过程，得到D₁,D₂,…,D_q,…,D_Q、C₁,C₂,…,C_q,…,C_Q、和然后执行步骤否则，令q＝q+1，然后返回步骤⑧继续执行；其中，q＝q+1中的“＝”为赋值符号，和中均含有虚假模态。

对D₁,D₂,…,D_q,…,D_Q各自中的部分元素强制置零，再从D₁,D₂,…,D_q,…,D_Q中任意选择一个作为状态矩阵的最终特征值向量，记为D^*；同样，对C₁,C₂,…,C_q,…,C_Q各自中的部分元素强制置零，再从C₁,C₂,…,C_q,…,C_Q中任意选择一个作为输出矩阵的最终值，记为C^*。

并计算Q次执行共得到的Q个固有频率向量的均值向量Q个固有阻尼比向量的均值向量和Q个固有模态振型系数向量的均值向量

在此具体实施例中，步骤中的D^*和C^*的获取过程为：

_1、在步骤的基础上，计算D₁,D₂,…,D_q,…,D_Q各自中的所有元素的均值，将D_q中的所有元素的均值记为同样，计算C₁,C₂,…,C_q,…,C_Q各自中的所有元素的均值，将C_q中的所有元素的均值记为

_2、比较D₁,D₂,…,D_q,…,D_Q各自中的每个元素与对应的均值的大小，将D₁,D₂,…,D_q,…,D_Q各自中小于对应的均值的元素强制置零；对于D_q，比较D_q中的每个元素与的大小，将D_q中小于的元素强制置零。

同样，比较C₁,C₂,…,C_q,…,C_Q各自中的每个元素与对应的均值的大小，将C₁,C₂,…,C_q,…,C_Q各自中小于对应的均值的元素强制置零；对于C_q，比较C_q中的每个元素与的大小，将C_q中小于的元素强制置零。

_3、在步骤_2的基础上，统计D₁,D₂,…,D_q,…,D_Q中相同索引值的Q个元素中的非零元素的总个数，将D₁,D₂,…,D_q,…,D_Q中索引值为g的Q个元素中的非零元素的总个数记为d_Dg；然后将统计得到的所有数据按序构成一个集合，记为d_D，d_D＝{d_D1,d_D2,…,d_Dg,…,d_DG}。

同样，统计C₁,C₂,…,C_q,…,C_Q中相同索引值的Q个元素中的非零元素的总个数，将C₁,C₂,…,C_q,…,C_Q中索引值为g'的Q个元素中的非零元素的总个数记为d_Cg'；然后将统计得到的所有数据按序构成一个集合，记为d_C，d_C＝{d_C1,d_C2,…,d_Cg',…,d_CG'}。

其中，g的初始值为1，1≤g≤G，G表示D₁,D₂,…,D_q,…,D_Q各自中包含的元素的总个数，d_D1表示D₁,D₂,…,D_q,…,D_Q中索引值为1的Q个元素中的非零元素的总个数，d_D2表示D₁,D₂,…,D_q,…,D_Q中索引值为2的Q个元素中的非零元素的总个数，d_DG表示D₁,D₂,…,D_q,…,D_Q中索引值为G的Q个元素中的非零元素的总个数，g'的初始值为1，1≤g'≤G'，G'表示C₁,C₂,…,C_q,…,C_Q各自中包含的元素的总个数，d_C1表示C₁,C₂,…,C_q,…,C_Q中索引值为1的Q个元素中的非零元素的总个数，d_C2表示C₁,C₂,…,C_q,…,C_Q中索引值为2的Q个元素中的非零元素的总个数，d_CG'表示C₁,C₂,…,C_q,…,C_Q中索引值为G'的Q个元素中的非零元素的总个数。

_4、计算d_D中的所有元素的均值，记为然后计算d_D中的每个元素与的差值的绝对值，将d_Dg与的差值的绝对值记为再将d_D中的所有元素与的差值的绝对值按序构成一个集合，记为f_D，

同样，计算d_C中的所有元素的均值，记为然后计算d_C中的每个元素与的差值的绝对值，将d_Cg'与的差值的绝对值记为再将d_C中的所有元素与的差值的绝对值按序构成一个集合，记为f_C，

其中，表示d_D1与的差值的绝对值，表示d_D2与的差值的绝对值，表示d_DG与的差值的绝对值，表示d_C1与的差值的绝对值，表示d_C2与的差值的绝对值，表示d_CG'与的差值的绝对值。

_5、在f_D中任选两个元素作为聚类中心；然后根据所选的两个聚类中心，利用K-means算法将f_D中的所有元素分为两个类；接着计算f_D的两个类各自中的所有元素的均值；之后提取出小值均值对应的类中的每个元素在f_D中的索引值；再将提取出的所有索引值按序构成一个集合，记为h_D。

同样，在f_C中任选两个元素作为聚类中心；然后根据所选的两个聚类中心，利用K-means算法将f_C中的所有元素分为两个类；接着计算f_C的两个类各自中的所有元素的均值；之后提取出小值均值对应的类中的每个元素在f_C中的索引值；再将提取出的所有索引值按序构成一个集合，记为h_C。

_6、将经步骤_2处理后得到的D₁,D₂,…,D_q,…,D_Q各自中索引值与h_D中的每个索引值相同的元素强制置零。

同样，将经步骤_2处理后得到的C₁,C₂,…,C_q,…,C_Q各自中索引值与h_C中的每个索引值相同的元素强制置零。

_7、在步骤_6的基础上，从D₁,D₂,…,D_q,…,D_Q中任选一个，假设选择D_q，则令D^*＝D_q；同样，从C₁,C₂,…,C_q,…,C_Q中任选一个，假设选择C_q，则令C^*＝C_q。

根据D^*中的每个零元素的索引，将中对应索引的模态参数确定为虚假模态；而根据D^*中的每个非零元素的索引，将中对应索引的模态参数确定为有效模态参数，并提取出。

根据D^*中的每个零元素的索引，将中对应索引的模态参数确定为虚假模态；而根据D^*中的每个非零元素的索引，将中对应索引的模态参数确定为有效模态参数，并提取出。

根据C^*中的每个零元素的索引，将中对应索引的模态参数确定为虚假模态；而根据C^*中的每个非零元素的索引，将中对应索引的模态参数确定为有效模态参数，并提取出。

为了进一步说明本发明方法的可行性和有效性，对本发明方法进行仿真对比试验。

在本实施例中的自由响应仿真信号为

y＝0.4e^-0.0101tcos(5.3t+0.01)+1.1e^-0.0567tcos(3.9t+0.02)+1.3e^-0.0827tcos(6.6t+0.05)，其中，y表示3阶连续的粘性阻尼振动信号，cos()为求余弦函数，t为时间变量，忽略y的初始相的影响，由自由响应仿真信号的仿真原理可得y的固有频率为：ω₁＝0.8439、ω₂＝0.6210、ω₃＝1.0510，固有阻尼比为：ξ₁＝0.012、ξ₂＝0.079、ξ₃＝0.091，模态振型系数为：Ω₁＝0.4、Ω₂＝1.1、Ω₃＝1.3。为了验证本发明方法的优越性，按照上述步骤，取k＝M＝24,i＝80。

1)以采样间隔为T_s＝0.067秒均匀采样2N＝1000个点，得到采样信号x，对x添加信噪比从0～30dB的噪声，由于篇幅过大，因此这里只取信噪比为5dB的高斯白噪声。

2)由于本发明方法中的步骤①至步骤⑦对不同的噪声过程都相同，因此在此不再赘述。令Q＝8，这里构建一个维数为24×1的第一高斯随机测量矩阵H₁，然后利用OMP方法得到Λ，进而得到A，再计算A的特征值向量D；重复构造第一高斯随机测量矩阵H₁共Q＝8次，并计算得到8个D、8个含有虚假模态的固有频率向量8个含有虚假模态的固有阻尼比向量和8个含有虚假模态的固有模态振型系数向量同样，构建一个维数为24×79的第二高斯随机测量矩阵H₂，然后利用OMP方法得到C；重复构造第二高斯随机测量矩阵H₂共Q＝8次，并计算得到8个C；

根据本发明方法的步骤计算8个固有频率向量的均值向量计算8个固有阻尼比向量的均值向量计算8个固有模态振型系数向量的均值向量 和的维数均为1×79；

8个D经过步骤_1至步骤_4的过程处理后，得到维数为1×79的统计结果向量f_D＝(0,0,5,2,0,0,3,",0,0,3)，然后利用K-means算法对f_D进行聚类操作，得到两个类，第一类为7、7、6，第二类为2、3、3、4、5、5，由于第二类的平均值小于第一类的平均值，因此提取出第二类中的每个元素在f_D中的索引值，再对第一次强制置零后的8个D中的第3、4、7、13、43、79列中的所有元素强制置零，此时8个D中只有第9、22、47列中包含有非零元素。

8个C经过步骤_1至步骤_4的过程处理后，利用K-means算法对f_C进行聚类操作，得到两个类，并提取出小值均值的类中的每个元素在f_C中的索引值，再对第一次强制置零后的8个C中的部分元素强制置零，此时8个C中只有第5、26、53列中包含有非零元素。

3)从中选出第9、22、47个元素，得到粘性阻尼振动信号中的固有频率为0.8501、0.6217、1.0523，从中选出第9、22、47个元素，得到粘性阻尼振动信号中的固有阻尼比为0.015、0.076、0.093，从中选出第5、26、53个元素，得到粘性阻尼振动信号中的固有模态振型系数为0.405、1.106、1.292。比较y的固有频率ω₁＝0.8439、ω₂＝0.6210、ω₃＝1.0510与本实验中得到的固有频率0.8501、0.6217、1.0523，比较y的固有阻尼比ξ₁＝0.012、ξ₂＝0.079、ξ₃＝0.091与本实验中得到的固有阻尼比0.015、0.076、0.093，比较y的模态振型系数Ω₁＝0.4、Ω₂＝1.1、Ω₃＝1.3与本实验中得到的固有模态振型系数0.405、1.106、1.292，可以看出本发明方法得到的参数较系统的理论值是较为精确的。

为了更好的证明本发明方法的优越性，对本发明方法与随机子空间法进行仿真对比。首先，仿真信号可知，仿真信号是一个V＝3阶的系统信号，然后由随机子空间法以3阶的系统阶次算出相应信噪比下的固有频率向量的均值向量固有阻尼比向量的均值向量和固有模态振型系数向量的均值向量而本发明方法将系统阶次的初始值赋值为80，在同样的信噪比的情况下用本发明方法获得的固有频率向量的均值向量固有阻尼比向量的均值向量和固有模态振型系数向量的均值向量

利用公式计算使用本发明方法和随机子空间法提取的固有频率的值与固有频率的理论值之间的相对误差，利用公式计算使用本发明方法和随机子空间法提取的固有阻尼比的值与固有阻尼比的理论值之间的相对误差，利用公式计算使用本发明方法和随机子空间法提取的固有模态振型系数的值与固有模态振型系数的理论值之间的相对误差，其中，符号“||”为取绝对值符号，γ_ω表示本发明方法和随机子空间法提取的固有频率与理论值之间的相对误差，V＝3为目标系统的固有阶次，ω_v为理论计算得到的第v个固有频率，为利用本发明方法和随机子空间法提取的第v个固有频率，γ_ξ表示本发明方法和随机子空间法提取的固有阻尼比与理论值之间的相对误差，ξ_v为理论计算得到的第v个固有阻尼比，为利用本发明方法和随机子空间法提取的第v个固有阻尼比，γΩ表示本发明方法和随机子空间法提取的固有模态振型系数与理论值之间的相对误差，Ω_v为理论计算得到的第v个固有模态振型系数，为利用本发明方法和随机子空间法提取的第v个固有模态振型系数。

图2a给出了在采样频率为15Hz、采样点数为1000、系统阶次的初始值为80、稀疏度的确定值为24、给采样信号添加高斯白噪声的情况下，高斯白噪声的信噪比逐次递增，利用本发明方法提取的固有频率与理论值之间的相对误差及利用随机子空间算法提取的固有频率与理论值之间的相对误差的对比；图2b给出了在采样频率为15Hz、采样点数为1000、系统阶次的初始值为80、稀疏度的确定值为24、给采样信号添加混合噪声(由脉冲噪声和高斯白噪声组成)，混合噪声的信噪比逐次递增，利用本发明方法提取的固有频率与理论值之间的相对误差及利用随机子空间算法提取的固有频率与理论值之间的相对误差的对比；图2c给出了在采样频率为15Hz、采样点数为1000、系统阶次的初始值为80、稀疏度的确定值为24、给采样信号添加混合噪声(由脉冲噪声和高斯白噪声组成)，混合噪声的信噪比逐次递增，利用本发明方法提取的固有阻尼比与理论值之间的相对误差及利用随机子空间算法提取的固有阻尼比与理论值之间的相对误差的对比；图2d给出了在采样频率为15Hz、采样点数为1000、系统阶次的初始值为80、稀疏度的确定值为24、给采样信号添加混合噪声(由脉冲噪声和高斯白噪声组成)，混合噪声的信噪比逐次递增，利用本发明方法提取的固有模态振型系数与理论值之间的相对误差及利用随机子空间算法提取的固有模态振型系数与理论值之间的相对误差的对比。在图2a至图2d中，SSI-cov代表随机子空间算法，稀疏优化代表本发明方法。

从图2a所示的曲线可以看出，在高斯白噪声下，利用随机子空间算法提取的固有频率与理论值之间的相对误差起伏较大，在信噪比10dB以内，波动特别明显，可看出对噪声较为敏感，而利用本发明方法提取的固有频率与理论值之间的相对误差起伏波动比较平缓，有很高的稳定性；从图2b所示的曲线可以看出，在脉冲噪声下，利用随机子空间算法提取的固有频率的精度比在高斯白噪声下利用随机子空间算法提取的固有频率的精度相对较差，对噪声也更加敏感，在信噪比15dB以内，相对误差变化起伏更大，而利用本发明方法提取的固有频率与理论值之间的相对误差变化平缓；从图2c所示的曲线可以看出，利用随机子空间算法提取的固有阻尼比在信噪比5dB以内，已经产生了明显的失真，当信噪比越小时，提取精度越差，而利用本发明方法提取的固有阻尼比的精度较高，表现出了很高的稳定性；从图2d所示的曲线可以看出，利用本发明方法提取的固有模态振型系数与理论值之间的相对误差始终在利用随机子空间算法提取的固有模态振型系数与理论值之间的相对误差之下，特别当信噪比较小时，误差值更为明显。

Claims (3)

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

①对连续的粘性阻尼振动信号进行奈奎斯特均匀采样，采样间隔为T_S秒，采样点数为2N点，得到包含2N个采样值的采样信号，记为x，将x以列向量形式表示为x＝(x₁,x₂,…,x_2N-1,x_2N)^T；其中，T_S的取值满足奈奎斯特采样定律，N的取值要求大于目标系统的系统阶次的预估值，在区间[200,500]内取一个正整数作为N的具体值，(x₁,x₂,…,x_2N-1,x_2N)^T为(x₁,x₂,…,x_2N-1,x_2N)的转置，x₁,x₂,…,x_2N-1,x_2N对应表示x中的第1个采样值、第2个采样值、…、第2N-1个采样值、第2N个采样值；

②利用x中的所有采样值构建一个Hankel矩阵，记为H，H中的每一条副对角线上的元素都相等，然后将H分成两个子矩阵，记为Y_p和Y_f，其中，H的维数为2i×j，Y_p和Y_f的维数均为i×j，i与j满足关系：2i+j-1＝2N，i的取值为常数，且n＜i＜N，n表示目标系统的系统阶次的预估值，x₃、x_j、x_j+1、x_i、x_i+1、x_i+2、x_i+3、x_i+j-1、x_i+j、x_i+j+1、x_2i、x_2i+1、x_2i+j-1对应表示x中的第3个采样值、第j个采样值、第j+1个采样值、第i个采样值、第i+1个采样值、第i+2个采样值、第i+3个采样值、第i+j-1个采样值、第i+j个采样值、第i+j+1个采样值、第2i个采样值、第2i+1个采样值、第2i+j-1个采样值；

③利用Y_p和Y_f构建一个维数为i×i的Toeplitz矩阵，记为T，然后选取T的第1行至第i-1行构成维数为(i-1)×i的第一子矩阵，记为T₁；并选取T的第2行至第i行构成维数为(i-1)×i的第二子矩阵，记为T₂；其中，为Y_p的转置；

④根据状态空间方程的定义对T₁进行分解，得到T₁的可观矩阵和可控反转矩阵，对应记为Γ₁和Δ₁，T₁＝Γ₁Δ₁，Γ₁＝(C,CA,…,CA^i-2)^T；并对T₁进行SVD分解，得到T₁的第一正交矩阵、第二正交矩阵和正奇异矩阵组成的对角矩阵，对应记为U₁、V₁和S₁，T₁＝U₁S₁V₁ ^T；然后根据T₁＝Γ₁Δ₁和T₁＝U₁S₁V₁ ^T，令

同样，根据状态空间方程的定义对T₂进行分解，得到T₂的可观矩阵和可控反转矩阵，对应记为Γ₂和Δ₂，T₂＝Γ₂Δ₂，Γ₂＝(CA,CA²,…,CA^i-1)^T；并对T₂进行SVD分解，得到T₂的第一正交矩阵、第二正交矩阵和正奇异矩阵组成的对角矩阵，对应记为U₂、V₂和S₂，然后根据T₂＝Γ₂Δ₂和令

其中，Γ₁的维数为(i-1)×1，Δ₁的维数为1×i，C表示维数为1×(i-1)的输出矩阵，A表示维数为(i-1)×(i-1)的状态矩阵，(C,CA,…,CA^i-2)^T为(C,CA,…,CA^i-2)的转置，A^i-2为A的i-2次方，V₁ ^T为V₁的转置，为S₁的次方，Γ₂的维数为(i-1)×1，Δ₂的维数为1×i，(CA,CA²,…,CA^i-1)^T为(CA,CA²,…,CA^i-1)的转置，A²为A的2次方，A^i-1为A的i-1次方，为V₂的转置，为S₂的次方；

⑤根据Γ₁＝(C_,CA,…,CA^i-2)^T和Γ₂＝(CA,CA²,…,CA^i-1)^T，确定Γ₁与Γ₂之间的关系为：Γ₂ ^T＝Γ₁ ^TΛ；然后将Γ₁ ^T分解为其中，Γ₂ ^T为Γ₂的转置，Γ₁ ^T为Γ₁的转置，Λ表示维数为(i-1)×(i-1)的对角矩阵，P₁表示维数为1×(i-1)的行向量，P₁＝(C,C,…,C)，Q₁表示维数为(i-1)×(i-1)的对角矩阵，I表示i-1阶的单位矩阵；

⑥根据稀疏优化的原理，确定稀疏度小于或等于i；并根据目标系统的系统阶次的预估值，确定稀疏度大于2n；然后确定稀疏度的取值范围为大于2n且小于或等于i；接着在稀疏度的取值范围内选择一个正整数k作为稀疏度的确定值，k∈(2n,i]；再将大于2n的值i赋值给目标系统的系统阶次作为目标系统的系统阶次的初始值；

⑦令q表示执行的次数，令Q表示重复执行的总次数；其中，q的初始值为1，Q≥2；

⑧利用第q次执行时随机产生的第一高斯随机测量矩阵H₁对Γ₂ ^T＝Γ₁ ^TΛ进行观测，构建得到第q次执行时的第一稀疏优化模型，描述为：并利用第q次执行时随机产生的第二高斯随机测量矩阵H₂对进行观测，构建得到第q次执行时的第二稀疏优化模型，描述为：其中，H₁的维数为M×1，H₂的维数为M×(i-1)，M∈[k-10,k+10]，min为取最小值函数，符号“|| ||₁”为求矩阵的1-范数符号，s.t.表示“受约束于…”，符号“|| ||”为求欧氏距离符号，P₁ ^T为P₁的转置，Q₁ ^T为Q₁的转置，σ₁和σ₂均为常数；

⑨根据目标系统的系统阶次的初始值和稀疏度的确定值，利用OMP方法对进行稀疏求解，得到对角矩阵Λ；然后根据和求解得到的Λ，计算得到A，并令A_q＝A；之后对A_q进行特征值分解，得到A_q的特征值向量，记为D_q；接着根据和得到的A，计算得到Q₁；同样，根据目标系统的系统阶次的初始值和稀疏度的确定值，利用OMP方法对进行稀疏求解，得到P₁ ^T；最后根据P₁＝(C,C,…,C)得到C，并令C_q＝C；其中，A_q和C_q的初始值均为0；

⑩根据A_q和C_q计算第q次执行得到的连续的粘性阻尼振动信号的固有频率向量、固有阻尼比向量和固有模态振型系数向量，对应记为和其中，和的维数均为1×(i-1)；

判断q≥Q是否成立，如果成立，则结束重复执行的过程，得到D₁,D₂,…,D_q,…,D_Q、C₁,C₂,…,C_q,…,C_Q、和然后执行步骤否则，令q＝q+1，然后返回步骤⑧继续执行；其中，q＝q+1中的“＝”为赋值符号；

计算D₁,D₂,…,D_q,…,D_Q各自中的所有元素的均值，将D_q中的所有元素的均值记为比较D₁,D₂,…,D_q,…,D_Q各自中的每个元素与对应的均值的大小，将D₁,D₂,…,D_q,…,D_Q各自中小于对应的均值的元素强制置零；对于D_q，比较D_q中的每个元素与的大小，将D_q中小于的元素强制置零；再从D₁,D₂,…,D_q,…,D_Q中任意选择一个作为状态矩阵的最终特征值向量，记为D^*；同样，计算C₁,C₂,…,C_q,…,C_Q各自中的所有元素的均值，将C_q中的所有元素的均值记为比较C₁,C₂,…,C_q,…,C_Q各自中的每个元素与对应的均值的大小，将C₁,C₂,…,C_q,…,C_Q各自中小于对应的均值的元素强制置零；对于C_q，比较C_q中的每个元素与的大小，将C_q中小于的元素强制置零；再从C₁,C₂,…,C_q,…,C_Q中任意选择一个作为输出矩阵的最终值，记为C^*；

并计算Q次执行共得到的Q个固有频率向量的均值向量Q个固有阻尼比向量的均值向量和Q个固有模态振型系数向量的均值向量

根据D^*中的每个零元素的索引，将中对应索引的模态参数确定为虚假模态；而根据D^*中的每个非零元素的索引，将中对应索引的模态参数确定为有效模态参数，并提取出；

根据D^*中的每个零元素的索引，将中对应索引的模态参数确定为虚假模态；而根据D^*中的每个非零元素的索引，将中对应索引的模态参数确定为有效模态参数，并提取出；

根据C^*中的每个零元素的索引，将中对应索引的模态参数确定为虚假模态；而根据C^*中的每个非零元素的索引，将中对应索引的模态参数确定为有效模态参数，并提取出。

2.根据权利要求1所述的一种粘性阻尼振动信号中的模态参数提取方法，其特征在于所述的步骤⑧中的H₁和H₂中的每个元素独立的服从均值为0且方差为的高斯分布。

3.根据权利要求1所述的一种粘性阻尼振动信号中的模态参数提取方法，其特征在于所述的步骤中的D^*和C^*的获取过程为：

_1、在步骤的基础上，计算D₁,D₂,…,D_q,…,D_Q各自中的所有元素的均值，将D_q中的所有元素的均值记为同样，计算C₁,C₂,…,C_q,…,C_Q各自中的所有元素的均值，将C_q中的所有元素的均值记为

_2、比较D₁,D₂,…,D_q,…,D_Q各自中的每个元素与对应的均值的大小，将D₁,D₂,…,D_q,…,D_Q各自中小于对应的均值的元素强制置零；对于D_q，比较D_q中的每个元素与的大小，将D_q中小于的元素强制置零；

同样，比较C₁,C₂,…,C_q,…,C_Q各自中的每个元素与对应的均值的大小，将C₁,C₂,…,C_q,…,C_Q各自中小于对应的均值的元素强制置零；对于C_q，比较C_q中的每个元素与的大小，将C_q中小于的元素强制置零；

_3、在步骤_2的基础上，统计D₁,D₂,…,D_q,…,D_Q中相同索引值的Q个元素中的非零元素的总个数，将D₁,D₂,…,D_q,…,D_Q中索引值为g的Q个元素中的非零元素的总个数记为d_Dg；然后将统计得到的所有数据按序构成一个集合，记为d_D，d_D＝{d_D1,d_D2,…,d_Dg,…,d_DG}；

同样，统计C₁,C₂,…,C_q,…,C_Q中相同索引值的Q个元素中的非零元素的总个数，将C₁,C₂,…,C_q,…,C_Q中索引值为g'的Q个元素中的非零元素的总个数记为d_Cg'；然后将统计得到的所有数据按序构成一个集合，记为d_C，d_C＝{d_C1,d_C2,…,d_Cg',…,d_CG'}；

其中，g的初始值为1，1≤g≤G，G表示D₁,D₂,…,D_q,…,D_Q各自中包含的元素的总个数，d_D1表示D₁,D₂,…,D_q,…,D_Q中索引值为1的Q个元素中的非零元素的总个数，d_D2表示D₁,D₂,…,D_q,…,D_Q中索引值为2的Q个元素中的非零元素的总个数，d_DG表示D₁,D₂,…,D_q,…,D_Q中索引值为G的Q个元素中的非零元素的总个数，g'的初始值为1，1≤g'≤G'，G'表示C₁,C₂,…,C_q,…,C_Q各自中包含的元素的总个数，d_C1表示C₁,C₂,…,C_q,…,C_Q中索引值为1的Q个元素中的非零元素的总个数，d_C2表示C₁,C₂,…,C_q,…,C_Q中索引值为2的Q个元素中的非零元素的总个数，d_CG'表示C₁,C₂,…,C_q,…,C_Q中索引值为G'的Q个元素中的非零元素的总个数；

_4、计算d_D中的所有元素的均值，记为然后计算d_D中的每个元素与的差值的绝对值，将d_Dg与的差值的绝对值记为再将d_D中的所有元素与的差值的绝对值按序构成一个集合，记为f_D，

同样，计算d_C中的所有元素的均值，记为然后计算d_C中的每个元素与的差值的绝对值，将d_Cg'与的差值的绝对值记为再将d_C中的所有元素与的差值的绝对值按序构成一个集合，记为f_C，

其中，表示d_D1与的差值的绝对值，表示d_D2与的差值的绝对值，表示d_DG与的差值的绝对值，表示d_C1与的差值的绝对值，表示d_C2与的差值的绝对值，表示d_CG'与的差值的绝对值；

_5、在f_D中任选两个元素作为聚类中心；然后根据所选的两个聚类中心，利用K-means算法将f_D中的所有元素分为两个类；接着计算f_D的两个类各自中的所有元素的均值；之后提取出小值均值对应的类中的每个元素在f_D中的索引值；再将提取出的所有索引值按序构成一个集合，记为h_D；

同样，在f_C中任选两个元素作为聚类中心；然后根据所选的两个聚类中心，利用K-means算法将f_C中的所有元素分为两个类；接着计算f_C的两个类各自中的所有元素的均值；之后提取出小值均值对应的类中的每个元素在f_C中的索引值；再将提取出的所有索引值按序构成一个集合，记为h_C；

_6、将经步骤_2处理后得到的D₁,D₂,…,D_q,…,D_Q各自中索引值与h_D中的每个索引值相同的元素强制置零；

同样，将经步骤_2处理后得到的C₁,C₂,…,C_q,…,C_Q各自中索引值与h_C中的每个索引值相同的元素强制置零；

_7、在步骤_6的基础上，从D₁,D₂,…,D_q,…,D_Q中任选一个，若选择D_q，则令D^*＝D_q；同样，从C₁,C₂,…,C_q,…,C_Q中任选一个，若选择C_q，则令C^*＝C_q。