CN106096562B - 基于振动信号盲源分离和稀疏成分分析的风电机组齿轮箱故障诊断方法 - Google Patents

Abstract

本发明公开了一种基于振动信号盲源分离和稀疏成分分析的风电机组齿轮箱故障诊断方法。本发明主要包括三部分的算法：一是基于经验模态分解(EMD)、奇异值分解(SVD)和K均值聚类(K‑Means)的源信号数目估计算法；二是基于模糊C均值聚类(FCM)的混叠矩阵估计算法；三是基于最小化l₁范数的源信号估计及诊断算法。整个算法流程引入了盲源分离(BSS)及稀疏成分分析(SCA)的思想。本发明以仿真信号和实际振动信号为测试对象，给出了详细的算法描述，并通过一系列的实验验证了算法在信号处理及故障诊断方面的有效性。

Description

基于振动信号盲源分离和稀疏成分分析的风电机组齿轮箱故 障诊断方法

技术领域

本发明属于振动信号处理和故障诊断方法领域，尤其涉及一种针对风力发电机组齿轮箱的振动信号处理步骤及故障诊断方法。

背景技术

风力发电机组昂贵的运行维护费用是阻碍风电产业快速发展的重要因素之一。随着风电单机容量不断增大，风力发电机组的体积和轮毂高度也不断增大，内部的传动系统受力情况更加复杂，由于风力发电机组故障而引起的事故常有发生，造成巨大的经济损失。齿轮箱位于风力发电机机舱内，是风力发电机传动动力的主要部件和连接主轴和发电机的重要枢纽，它具有结构紧凑、传递精度高、传递力矩大等特点，是风力发电机正常高效运行的保障。齿轮箱内部结构和受力情况复杂，并且常在变工况、变载荷等复杂环境下工作，因此齿轮箱部件在长时间运行过程中极易老化损伤，产生各类故障。因此，研究一套可行的齿轮箱故障诊断方法具有重要意义。

盲源分离技术是目前齿轮箱故障诊断领域较常采用的方法。齿轮箱的故障诊断的实质就是从齿轮箱各部位的振动数据中推断出是否发生故障以及发生故障情况下故障的具体类型。在上述诊断过程中，源信号的数目以及和源信号的信号特征都是未知的，即满足盲源特点，因此将盲源分离技术引入故障诊断是十分适合的并已存在相当多的成功应用案例。然而，现有的盲源分离技术仍然存在一定的局限性：一方面，目前的盲源分离算法大多针对超定或者正定情况，即要求观测信号数目大于或等于源信号数目，对欠定情况(观测信号数目小于源信号数目)下的盲源分离研究。但在很多情况下，这一条件并不满足，从而造成分离出的源信号仍然存在一定的混叠，影响最终的诊断结果。另一方面，真实的齿轮箱故障信号具有非平稳、非线性等特点，现有方法对处理这种信号的表现通常较差。

发明内容

为了改进现有方法，本发明的目的在于提供一种适用于欠定非线性盲源情况下的齿轮箱振动信号处理及故障诊断方法，提高故障诊断的准确性。

本发明的目的通过下述技术方案实现，包括以下步骤：

步骤(1)，从安装在齿轮箱内的m个传感器上获取振动信号，并设振动信号矩阵为x＝[x₁,x₂,...x_m]，其中任意一路传感器信号x_i都为T维向量，即x_i∈R^T，表示每一路信号有T个采样点，并记采样频率为Fs。

步骤(2)选取任意一路传感器信号x_i，对其进行EMD分解，得到IMF矩阵C，EMD分解过程可采用如下方法：

(2-1)初始化r₀＝x_i，q＝1(q为IMF下标，指示IMF数目)；

(2-2)初始化h₀＝r_q-1，k＝1；

(2-3)求取h_k-1的局部极大值和局部极小值点，并以局部极大值点为节点作三次样条插值计算h_k-1上包络u_k-1，以局部极小值点为节点作三次样条插值计算h_k-1下包络l_k-1；

(2-4)计算上下包络的均值

(2-5)令h_k＝h_k-1-m_k-1；

(2-6)若h_k满足本征模态函数条件，则令第q个本征模态函数c_q＝h_k；否则，令k＝k+1并转入(2-3)；

(2-7)令残量r_q＝r_q-1-c_q，若r_j至少包含两个极值点，则令q＝q+1并转入(2-2)；否则分解结束，r_q为余量。

(2-8)将上述步骤分离出的q个本征模态函数及余量r_q组合成IMF矩阵C＝[c₁,c₂,...c_q,r_q]。

步骤(3)计算IMF矩阵C的自相关矩阵R_C＝CC^T。

步骤(4)将协方差矩阵进行奇异值分解，由于协方差矩阵为q×q维矩阵，因此得到协方差矩阵q个奇异值，奇异值构成的集合记为{λ₁,λ₂,...λ_q}。

步骤(5)取奇异值的自然对数{logλ₁,logλ₂,...logλ_q}，并对该集合进行K均值聚类，分类数为2，聚类过程可以如下：

(5-1)取两个随机数作为聚类中心，记为{a₁,a₂}；

(5-2)建立目标函数

(5-3)最小化目标函数W，在最小化的过程中，聚类中心{a₁,a₂}会不断被更新，直至收敛；

(5-4)聚类中心收敛后，待分类集合{logλ₁,logλ₂,...logλ_q}会被分为两类，将含数值较大的一类记为G₁，数值较小的一类记为G₂

步骤(6)统计G₁中的元素个数n，将其作为源数目的估计。步骤(1)到步骤(6)的最终目的就是为了得到源数目估计n，用于后续的混叠矩阵及源信号估计。

步骤(7)选取短时傅立叶变换长度L≥64，窗口重叠长度为将各传感器振动信号{x₁,x₂,...x_m}分段，并将各段信号依次进行短时傅立叶变换，将原始信号转化到时频域，得到系数矩阵为X(L,K,m)，即X为L×K×m的三维矩阵，其中表示分段数。

步骤(8)根据系数矩阵X(L,K,m)计算各频率点的能量值，能量值矩阵记为E(L)，E(L)中的元素

步骤(9)在E(L)中选取n个极大值，记录下相应的频率点，频率点集合记为{f₁,f₂,...f_n}。注意，这里的频率点实际上是极大值在E(L)中的坐标，与真实的频率成比例关系，并非真实的频率值。

步骤(10)在频率点集合{f₁,f₂,...f_n}中任取一个频率点f_i，并将其从频率点集合中剔除，从系数矩阵X(L,K,m)提取此频率点的系数矩阵(只取实部)并将其归一化，公式为注意Y为K×m维矩阵。

步骤(11)利用模糊C均值算法对Y进行聚类，得到聚类中心{v₁,v₂}。矩阵Y中的每一行都可以看作一个m维的样本点，由于Y有K行，样本点集合可记为{y₁,y₂,...y_K}。聚类过程可以如下：

(11-1)设定迭代终止误差为ε，设定权重系数为w,将初始化分类数c设定为2，随机初始化隶属度矩阵U＝[u_ij]，隶属度矩阵中的元素u_ij代表样本y_j在第i类中的隶属度，并满足

(11-2)根据公式计算聚类中心{v₁,v₂}；

(11-3)根据公式计算新的隶属度矩阵，记为U′；

(11-4)判断是否满足迭代终止准则，即若||U-U′||＜ε，则停止计算；若||U-U′||≥ε，则令U＝U′并转入(11-2)，继续迭代。

步骤(12)将聚类中心v₁作为混叠矩阵A的一列。

步骤(13)判断频率点集合{f₁,f₂,...f_n}是否为空，假如集合不空，则转入步骤(10)；假如集合为空集，表明已经得到混叠矩阵A的估计，此时混叠矩阵A恰好为m×n维的矩阵。步骤(7)到步骤(13)的最终目的就是得到混叠矩阵的估计。

步骤(14)分别从系数矩阵X(L,K,m)中提取各元素的实部和虚部并各自压缩重组为一个新的二维矩阵，实部矩阵记为RealX(m,L×K)，虚部矩阵记为ImgX(m,L×K)。

步骤(15)建立目标函数且满足As^t＝x^t，其中A为步骤(13)估计得到的混叠矩阵，s^t为源信号在时频域中的一列，x^t为输入矩阵的一列，这里输入矩阵可以是RealX(m,L×K)或者ImgX(m,L×K)。该步骤建立的目标函数的意义在于最小化分段频域源信号矩阵的l₁范数，从而得到分段时频域源信号的估计。

步骤(16)最小化目标函数，得到分段时频域源信号的估计。为了使源信号更加直观和方便诊断，对分段时频域源信号执行逆短时傅立叶变化，得到时域源信号{s₁,s₂,...s_n}。

步骤(17)对源信号{s₁,s₂,...s_n}进行傅立叶变换，得到n个频域波形，判断各个频域波形的故障频率点是否存在明显波形，假如存在，则说明存在该类型故障，否则反之。

本发明引入了盲源分离技术来解决齿轮箱故障诊断问题。齿轮箱振动信号具有时变非线性的特点，利用盲源分离算法直接分离时域信号通常存在比较严重的混叠现象。实际上，齿轮箱振动信号的频率成分与故障类型息息相关，说明信号在频域上能更直观反映故障诊断的本质；另一方面，结合振动信号时变的特点，频率组成比例会随时间产生波动，因此需要结合时域和频域综合考虑。同时，本发明算法引入了稀疏成分分析，利用齿轮箱振动信号在时频域上稀疏的特性，将传感器时域振动信号通过短时傅立叶变换转化时频域求解，随后将计算结果转换回时域或者频域分析，充分利用了振动信号的特性。

本发明与现有振动信号处理方法相比，有益效果包括：本发明将时域振动信号转化到时频域进行求解，更加适用于非线性、非平稳信号的处理；本发明对传感器数目没有要求，适用于超定、正定以及欠定盲源分离问题，弥补了以独立成分分析为基础的盲源分离算法无法处理欠定盲源问题的不足之处；方法的鲁棒性较强，不容易受到噪音干扰，同时减弱了信号间的混叠现象，更加突出一些较微弱的故障信号，避免一些较微弱的故障信号被无用信号覆盖从而造成错判，利于早期故障的发现与维护，从而降低了风场的维护成本。

附图说明

图1算法框架流程图；

图2源信号数目估计算法示意图；

图3不同信噪比下源信号数目估计算法的准确度；

图4五路模拟齿轮箱振动信号的仿真源信号；

图5两路线性混叠情况下的仿真传感器信号；

图6线性混叠情况下第一路传感器信号经验模态分解得到的本征模态函数及余量；

图7五路线性混叠情况下的源信号估计图；

图8两路非线性混叠情况下的仿真传感器信号；

图9五路非线性混叠情况下的源信号估计图；

图10两路真实传感器振动信号时域波形；

图11运用本发明算法恢复得到的三路源信号频域图。

具体实施方式

下面结合附图对本发明做进一步详述：

图1展示了本发明的整体算法框架，从图中可以看到，本发明主要包括三部分算法。在接收到安装在齿轮箱内的m个传感器振动信号之后，此时未知量有两个，即源信号的数量n及源信号本身，因此第一部分算法的目的就是估计源信号的数量n，这里采用的算法有经验模态分解(EMD)、奇异值分解(SVD)以及K均值聚类(K-Means)。图2具体展示了第一部分算法：任意选择一路传感器信号进行EMD分解，分解结束后会得到一个由本征模态函数组成的IMF矩阵，求取IMF自相关矩阵的奇异值并利用K均值聚类算法二分类，数值较大一类的奇异值个数就是源信号数量n的估计。确定了源信号数量n实质上就确定了混叠矩阵的维数，因而可以执行第二部分估计混叠矩阵的算法，这里用到的主要方法是模糊C均值聚类(FCM)，该部分的算法实际上还利用了振动信号在时频域上的稀疏性。混叠矩阵估计完毕后，第三部分算法开始执行，本发明运用最小化l₁范数的方法还原原始信号。之后通过寻找故障频率点来确定故障类型。

通过上述分析可知，准确估计源信号的数量n是保证后续算法有效的前提，错误的估计直接影响到源信号的恢复。为了验证本发明算法的有效性，通过Matlab构造出若干个仿真源信号来模拟振动源信号，混叠矩阵随机选取，将仿真源信号混合并加以不同强度的噪声，每改变1dB的信噪比(SNR)便运用本发明算法重复一百次仿真实验(不同仿真实验的混叠矩阵都是重新随机选取的)，统计得到不同信噪比下源信号数目估计算法的准确度，如图3所示。从图中可以看到，当信噪比升高时，算法的估计准确度随之提高，最终收敛为1，表示信噪比达到一定值后预测准确度为100％，同时，当信噪比在3％附近时，预测准确度也能保持在50％左右。

图4为五路模拟齿轮箱振动信号的仿真源信号时域图，各路信号分别表示齿轮箱内某个轴承或者齿轮产生的振动冲击信号，各路信号之间可以互相混叠，振动传感器接收到的信号就是它们混叠后的信号。图5和图8分别是五路模拟源信号在线性混叠情况与非线性混叠情况下得到的传感器信号，分别对它们执行以下步骤：

步骤(1)，由于传感器信号有两路，因此m＝2，并设振动信号矩阵为x＝[x₁,x₂]，其中任意一路传感器信号x_i都为1024维向量，表示每一路信号有1024个采样点，采样频率为1Hz。

步骤(2)选取任意一路传感器信号x_i，这里选取x₁，对其进行EMD分解，得到IMF矩阵C，线性混叠情况下第一路传感器信号经验模态分解得到的本征模态函数及余量见图6，EMD分解过程如下：

(2-1)初始化r₀＝x_i，q＝1(q为IMF下标，指示IMF数目)；

(2-2)初始化h₀＝r_q-1，k＝1；

(2-3)求取h_k-1的局部极大值和局部极小值点，并以局部极大值点为节点作三次样条插值计算h_k-1上包络u_k-1，以局部极小值点为节点作三次样条插值计算h_k-1下包络l_k-1；

(2-4)计算上下包络的均值

(2-5)令h_k＝h_k-1-m_k-1；

(2-6)若h_k满足本征模态函数条件，则令第q个本征模态函数c_q＝h_k；否则，令k＝k+1并转入(2-3)；

(2-7)令残量r_q＝r_q-1-c_q，若r_j至少包含两个极值点，则令q＝q+1并转入(2-2)；否则分解结束，r_q为余量。

(2-8)将上述步骤分离出的q个本征模态函数及余量r_q组合成IMF矩阵C＝[c₁,c₂,...c_q,r_q]。

步骤(3)计算IMF矩阵C的自相关矩阵R_C＝CC^T。

步骤(4)将协方差矩阵进行奇异值分解，由于协方差矩阵为9×9维矩阵，因此得到协方差矩阵9个奇异值，奇异值构成的集合记为{λ₁,λ₂,...λ₉}。

步骤(5)取奇异值的自然对数{logλ₁,logλ₂,...logλ₉}，并对该集合进行K均值聚类，分类数为2，聚类过程如下：

(5-1)取两个随机数作为聚类中心，记为{a₁,a₂}；

(5-2)建立目标函数

(5-3)最小化目标函数W，在最小化的过程中，聚类中心{a₁,a₂}会不断被更新，直至收敛；

(5-4)聚类中心收敛后，待分类集合{logλ₁,logλ₂,...logλ_q}会被分为两类，将含数值较大的一类记为G₁，数值较小的一类记为G₂

步骤(6)统计G₁中的元素个数，得到n＝5，将其作为源数目的估计。

步骤(7)设定短时傅立叶变换长度L＝512，窗口重叠长度Overlap＝504，将各传感器振动信号{x₁,x₂}分段，并将各段信号依次进行短时傅立叶变换，将原始信号转化到时频域，得到系数矩阵为X(512,65,2)。

步骤(8)根据系数矩阵X(512,65,2)计算各频率点的能量值，能量值矩阵记为E(L)，E(L)中的元素

步骤(9)在E(L)中选取5个极大值，记录下相应的频率点，频率点集合记为{f₁,f₂,...f₅}。注意，这里的频率点实际上是极大值在E(L)中的坐标，与真实的频率成比例关系，并非真实的频率值。

步骤(10)在频率点集合{f₁,f₂,...f₅}中任取一个频率点f_i，并将其从频率点集合中剔除，从系数矩阵X(512,65,2)提取此频率点的系数矩阵(只取实部)并将其归一化，公式为注意Y为65×2维矩阵。

步骤(11)利用模糊C均值算法对Y进行聚类，得到聚类中心{v₁,v₂}。矩阵Y中的每一行都可以看作一个2维的样本点，由于Y有K＝65行，样本点集合可记为{y₁,y₂,...y_K}。聚类过程如下：

(11-1)设定迭代终止误差为ε＝1e-5，设定权重系数为w＝2.0,将初始化分类数c设定为2，随机初始化隶属度矩阵U＝[u_ij]，隶属度矩阵中的元素u_ij代表样本y_j在第i类中的隶属度，并满足

(11-2)根据公式计算聚类中心{v₁,v₂}；

(11-3)根据公式计算新的隶属度矩阵，记为U′；

(11-4)判断是否满足迭代终止准则，即若||U-U′||＜ε，则停止计算；若||U-U′||≥ε，则令U＝U′并转入(11-2)，继续迭代。

步骤(12)将聚类中心v₁作为混叠矩阵A的一列。

步骤(13)判断频率点集合{f₁,f₂,...f_n}是否为空，假如集合不空，则转入步骤(10)；假如集合为空集，表明已经得到混叠矩阵A的估计，此时混叠矩阵A恰好为2×5维的矩阵。步骤(7)到步骤(13)的最终目的就是得到混叠矩阵的估计。

步骤(14)分别从系数矩阵X(512,65,2)中提取各元素的实部和虚部并各自压缩重组为一个新的二维矩阵，实部矩阵记为RealX(2,512×65)，虚部矩阵记为ImgX(2,512×65)，表示两者都是2×33280维矩阵。

步骤(15)建立目标函数且满足As^t＝x^t，其中A为步骤(13)估计得到的混叠矩阵，s^t为源信号在时频域中的一列，x^t为输入矩阵的一列，这里输入矩阵可以是RealX(2,512×65)或者ImgX(2,512×65)。该步骤建立的目标函数的意义在于最小化分段频域源信号矩阵的l₁范数，从而得到分段时频域源信号的估计。

步骤(16)最小化目标函数，得到分段时频域源信号的估计。为了使源信号更加直观和方便诊断，对分段时频域源信号执行逆短时傅立叶变化，得到时域源信号{s₁,s₂,...s₅}。

步骤(17)对源信号{s₁,s₂,...s₅}进行傅立叶变换，得到5个频域波形，判断各个频域波形的故障频率点是否存在明显波形，假如存在，则说明存在该类型故障，否则反之。由于本实例中仿真源信号已知，不需要再转换到频域，可以直接观察恢复信号与源信号的时域波形。

图7是将线性混叠的传感器信号代入上述步骤得到五路的源信号估计，对比图4中的原始信号，需要注意的是，本发明提出的算法并不保证恢复信号保持原始信号的次序性，即图7中的e₁、e₂、e₃、e₄和e₅分别是图4中的s₁、s₄、s₂、s₅和s₃的估计。观察对应的信号，可以发现估计信号较好地恢复了原始信号，从时域波形上几乎无法分辨两者的区别。为了定量衡量算法对源信号的恢复程度，可以计算估计信号与对应信号的相关系数，见表I。图9则是将非线性混叠的传感器信号代入上述步骤得到五路的源信号估计，对比图4可以看到非线性情况下的分离效果相对于线性情况效果稍差，说明非线性情况下的恢复信号仍然存在一定程度的混叠，但还是能分辨出大部分信号对应的源信号。

表I线性混叠情况下恢复信号与对应原始信号相关系数

上述实施例中用到的信号都是仿真信号，下面采用真实振动数据对算法进行验证。振动信号来自美国凯斯西储大学轴承数据中心，实施过程与上述过程类似，这里不再赘述。图10为两路真实传感器振动信号时域波形，采样频率为12000Hz，电动机转速为1772rpm，轴承类型为深沟球轴承，具体型号为6205-2RS JEM SKF，故障类型为轴承内圈故障，理论故障频率为159.93Hz。图11为运用本发明算法恢复得到的三路源信号频域图。从图中可以看出，e₁、e₂和e₃分别占据了低频段、中频段和高频段，并且在e₁的159.7Hz频率点附近存在明显的波峰，说明存在内圈故障。

应当理解的是，上述实施例为本发明较佳的实施例子，但本发明的实施方式并不受上述实施例的限制。对本领域普通技术人员来说，可以根据上述说明加以改进或变换，而所有这些改进和变换都应属于本发明所附权利要求的保护范围。

Claims (4)

1.基于振动信号盲源分离和稀疏成分分析的风电机组齿轮箱故障诊断方法,其特征在于：

步骤(1)，从安装在齿轮箱内的m个传感器上获取振动信号，并设振动信号矩阵为x＝[x₁,x₂,...x_m]，其中任意一路传感器信号x_i都为T维向量，即x_i∈R^T，表示每一路信号有T个采样点，并记采样频率为Fs；

步骤(2)选取任意一路传感器信号x_i，对其进行EMD分解，分离出的q个本征模态函数及余量r_q，得到本征模态函数IMF矩阵C，IMF矩阵C＝[c₁,c₂,...c_q,r_q]；

步骤(3)计算IMF矩阵C的自相关矩阵R_C＝CC^T；

步骤(4)将协方差矩阵进行奇异值分解，由于协方差矩阵为q×q维矩阵，因此得到协方差矩阵q个奇异值，奇异值构成的集合记为{λ₁,λ₂,...λ_q}；

步骤(5)取奇异值的自然对数{logλ₁,logλ₂,...logλ_q}，并对该集合进行K均值聚类，分类数为2，将含数值较大的一类记为G₁，数值较小的一类记为G₂；

步骤(6)统计G₁中的元素个数n，将其作为源数目的估计；

步骤(7)选取短时傅立叶变换长度L≥64，窗口重叠长度为将各传感器振动信号{x₁,x₂,...x_m}分段，并将各段信号依次进行短时傅立叶变换，将原始信号转化到时频域，得到系数矩阵为X(L,K,m)，即X为L×K×m的三维矩阵，其中表示分段数；

步骤(8)根据系数矩阵X(L,K,m)计算各频率点的能量值，能量值矩阵记为E(L)，E(L)中的元素

步骤(9)在E(L)中选取n个极大值，记录下相应的频率点，即各极大值在E(L)中的坐标，频率点集合记为{f₁,f₂,...f_n}；

步骤(10)在频率点集合{f₁,f₂,...f_n}中任取一个频率点f_i，并将其从频率点集合中剔除，从系数矩阵X(L,K,m)提取此频率点的系数矩阵的实部并将其归一化，公式为其中，Y为K×m维矩阵；

步骤(11)利用模糊C均值算法对Y进行聚类，得到聚类中心{v₁,v₂}， 矩阵Y中的每一行都可以看作一个m维的样本点，由于Y有K行，样本点集合可记为{y₁,y₂,...y_K}；

步骤(12)将聚类中心v₁作为混叠矩阵A的一列；

步骤(13)判断频率点集合{f₁,f₂,...f_n}是否为空，假如集合不空，则转入步骤(10)；假如集合为空集，表明已经得到混叠矩阵A的估计，此时混叠矩阵A恰好为m×n维的矩阵；

步骤(14)分别从系数矩阵X(L,K,m)中提取各元素的实部和虚部并各自压缩重组为一个新的二维矩阵，实部矩阵记为RealX(m,L×K)，虚部矩阵记为ImgX(m,L×K)；

步骤(15)建立目标函数且满足As^t＝x^t，其中A为步骤(13)估计得到的混叠矩阵，s^t为源信号在时频域中的一列，x^t为输入矩阵的一列，所述的输入矩阵是RealX(m,L×K)或者ImgX(m,L×K)；

步骤(16)最小化目标函数，得到分段时频域源信号的估计，为了使源信号更加直观和方便诊断，对分段时频域源信号执行逆短时傅立叶变化，得到时域源信号{s₁,s₂,...s_n}；

步骤(17)对源信号{s₁,s₂,...s_n}进行傅立叶变换，得到n个频域波形，判断各个频域波形的故障频率点是否存在明显波形，假如存在，则说明存在该类型故障，否则反之。

2.根据权利要求1所述的基于振动信号盲源分离和稀疏成分分析的风电机组齿轮箱故障诊断方法,其特征在于，所述的步骤(2)中EMD分解过程如下：

(2-1)初始化r₀＝x_i，q＝1，q为IMF下标，指示IMF数目；

(2-2)初始化h₀＝r_q-1，k＝1；

(2-3)求取h_k-1的局部极大值和局部极小值点，并以局部极大值点为节点作三次样条插值计算h_k-1上包络u_k-1，以局部极小值点为节点作三次样条插值计算h_k-1下包络l_k-1；

(2-4)计算上下包络的均值

(2-5)令h_k＝h_k-1-m_k-1；

(2-6)若h_k满足本征模态函数条件，则令第q个本征模态函数c_q＝h_k；否则，令k＝k+1并转入(2-3)；

(2-7)令残量r_q＝r_q-1-c_q，若r_j至少包含两个极值点，则令q＝q+1并转入(2-2)；否则分解结束，r_q为余量；

(2-8)将上述步骤分离出的q个本征模态函数及余量r_q组合成IMF矩阵C＝[c₁,c₂,...c_q,r_q]。

3.根据权利要求1所述的基于振动信号盲源分离和稀疏成分分析的风电机组齿轮箱故障诊断方法,其特征在于，所述的步骤(5)中聚类过程如下：

(5-1)取两个随机数作为聚类中心，记为{a₁,a₂}；

(5-2)建立目标函数

(5-3)最小化目标函数W，在最小化的过程中，聚类中心{a₁,a₂}会不断被更新，直至收敛；

(5-4)聚类中心收敛后，待分类集合{logλ₁,logλ₂,...logλ_q}会被分为两类，将含数值较大的一类记为G₁，数值较小的一类记为G₂。

4.根据权利要求1所述的基于振动信号盲源分离和稀疏成分分析的风电机组齿轮箱故障诊断方法,其特征在于，所述的步骤(11)中聚类过程如下：

(11-1)设定迭代终止误差为ε，设定权重系数为w,将初始化分类数c设定为2，随机初始化隶属度矩阵U＝[u_ij]，隶属度矩阵中的元素u_ij代表样本y_j在第i类中的隶属度，并满足

(11-2)根据公式计算聚类中心{v₁,v₂}；

(11-3)根据公式计算新的隶属度矩阵，记为U′；

(11-4)判断是否满足迭代终止准则，即若||U-U′||＜ε，则停止计算；若||U-U′||≥ε，则令U＝U′并转入(11-2)，继续迭代。