CN105806613A - 一种基于阶比复杂度的行星齿轮箱故障诊断方法 - Google Patents

Abstract

本发明属于旋转机械故障诊断技术领域，尤其涉及一种基于阶比复杂度的行星齿轮箱故障诊断方法，特别适用于变工况下的风电行星齿轮箱的故障诊断方法。包括：振动信号预处理；将行星轮系的故障分为两类；故障特征参数；对行星行星齿轮箱变工况下的振动信号进行分析；实例验证发现。本发明采用阶比分析方法可以避免变工况的影响，经过阶比分析处理后的角域信号的Lempel‑Ziv复杂度较时域信号更加平稳。行星齿轮箱振动角域信号的Lempel‑Ziv复杂度能够对行星齿轮局部故障和分布故障进行识别，通过仿真和实际案例的分析，验证了本发明的有效性，同时在缺少故障数据的情况下也可以作为故障发生与否的特征参数。

Description

一种基于阶比复杂度的行星齿轮箱故障诊断方法

技术领域

本发明属于旋转机械故障诊断技术领域，尤其涉及一种基于阶比复杂度的行星齿轮箱故障诊断方法，特别适用于变工况下的风电行星齿轮箱的故障诊断方法。

背景技术

行星齿轮箱具有传动比大、承载能力强、传动效率高等特点，广泛地应用于风力发电、直升机等大型复杂机械设备中。恶劣的工作环境，常常引发齿轮箱中关键部件的故障发生。针对定轴齿轮箱的故障诊断研究已取得了初步的成效，但行星齿轮箱不同于固定中心轴旋转的传统齿轮箱，行星齿轮箱中的齿轮运动是典型的复合运动，其振动响应比传统的齿轮箱更为复杂，故障的诊断与识别难度相较传统齿轮箱更大，传统齿轮箱故障诊断方法很难直接照搬应用于行星齿轮箱的故障诊断。

齿轮箱的故障诊断主要对振动信号进行分析，振动信号能够很好地描述齿轮振动的故障信息，但由于齿轮箱处于变工况的运行状态下，直接对振动信号进行分析存在一定的误差。在传统齿轮箱和轴承的振动信号的复杂度分析的基础上，针对行星齿轮箱变工况的运行条件，本发明采用阶比分析与EMD经验模态分解方法对振动信号进行预处理，分析行星齿轮箱不同运行状态下的复杂度，进而区分其故障状态。通过对风电场故障机组的振动信号进行分析，验证了本发明的有效性。

发明内容

针对现有技术的不足，本发明提供一种变工况行星齿轮箱故障诊断分析方法。

本发明的技术方案是这样实现的：

一种基于阶比复杂度的行星齿轮箱故障诊断方法，包括如下步骤：

(1)振动信号预处理：根据阶比跟踪技术，将处于变工况的风电行星齿轮箱传感器所采集的振动信号进行预处理，将非线性、非平稳的时域信号转化为具有平稳性的角域信号，避免了使用硬件方式实现等角度采样的昂贵成本；基于线性插值方法的非平稳振动时域信号的阶比重构技术，将等时间间隔采样的非平稳振动时域信号转化为具有平稳特性的角域振动信号，保证了行星齿轮箱振动角域信号的整周期性，避免了变工况运行环境的影响；

(2)将行星轮系的故障分为两类：分布故障和局部故障；对行星轮系的分布故障和局部故障的特征频率的计算进行分析；以华锐SL1500风电行星齿轮箱为研究对象，对两级行星轮系、一级平行级行星齿轮箱的振动信号进行仿真；

(3)故障特征参数：Lempel-Ziv复杂度在滚动轴承和简单平行级齿轮的故障诊断中已经表明，可以作为轴承故障诊断的故障特征参数，对Lempel-Ziv复杂度的概念进行分析，进一步确认Lempel-Ziv复杂度可以作为故障诊断的特征参数；

(4)对行星行星齿轮箱变工况下的振动信号进行分析，发现非平稳振动信号的Lempel-Ziv复杂度会有所波动，但经过阶比重采样处理后的角域信号可以有效避免了变工况对Lempel-Ziv复杂度的影响；由此可以验证，可以作为变工况下故障诊断特征参数；Lempel-Ziv复杂度可以作为恶劣工况环境下的行星齿轮箱的状态识别的特征参数；在步骤(1)、(2)的基础上，采用阶比复杂度对复杂结构的行星齿轮箱故障进行诊断，计算不同故障模式下的阶比复杂度，通过观察各故障所对应的阶比复杂度的变化情况，发现阶比复杂度可以作为行星齿轮箱故障模式识别的特征参数；该步骤利用阶比复度，通过对风电行星齿轮箱处于不同状态下的故障特征阶比复杂度进行计算，实现行星齿轮箱故障模式的识别；

(5)实例验证发现，行星齿轮箱在不同故障模式下的阶比复杂度相差较大，可以作为行星齿轮箱状态的评价，在缺少充足故障数据的情况下，也可以作为齿轮箱健康状态的评估；经该步骤实例验证，利用阶比复度，可以对风电行星齿轮箱故障状态进行识别，在缺少故障数据的情况下，也可以作为行星齿轮箱健康状态的评估参数。

步骤(1)所述的振动信号预处理：是对齿行星轮箱振动信号预处理；行星齿轮箱处于变工况的工作环境下，其采集的振动信号为非平稳信号，这对齿轮箱的故障诊断产生很大的困难，采用阶比重采样的方法对振动信号进行角域重采样，将非平稳的时域振动信号转化为平稳或准平稳特性的角域信号，便于对振动信号的分析；阶比分析技术的核心在于获得相对参考轴的恒定角增量采样数据，因此需要能准确获得阶次采样的时刻及相应的基准转速，即实现阶次跟踪；常见的阶比跟踪方法有硬件阶次跟踪法、计算阶比跟踪法和基于瞬时频率估计的阶比跟踪法；采用计算阶比跟踪法实现振动信号的重采样计算，通过重采样可将机械变速过程中产生的与转速有关的振动信号有效的分离出来，同时对与转速无关的信号起到一定的抑制作用，按转角位置重新分配信号的采样间隔，可以剔除转速变化对信号随机性的影响；实际的齿轮箱振动信号中一般情况下都含有多种干扰成分，这就使得其故障特征的提取变得比较困难；EMD方法可以根据信号的局部时变特性，自适应的将任意一个复杂信号分解为一系列分量，通过相关系数法则对信号进行重构，剔除原始信号中的干扰成分；采用阶比分析与EMD分解方法可以有效的对振动信号进行预处理，对下一步的信号分析做好了准备；

由于运行工况的复杂性和部件间的耦合性，风电机组的振动数据具有明显的非线性和非平稳特性，因此必须采用振动数据预处理技术对机组振动数据进行前置处理；本发明基于阶比重采样技术的风电机组振动数据处理方法，该方法可以有效地解决风电机组振动信号的变转速问题，可以将具有非平稳性特性的振动时域信号转化为具有平稳特性的振动角域信号，其基本原理是将采用软件算法将等时间间隔采样的时域信号转化为等角度采样模式下的角域振动信号，转化后的角域振动信号即保证了原始时域信号的真实性，有具有了正周期特性，在角域振动信号的基础上提取的振动特征提取能够准确有效地反映风电机组变工况运行条件下的故障状态；

其阶比重采样技术的主要步骤如下所示：

①假设振动测点获得的振动数据为{x₁,x₂,x₃,…x_n-2,x_n-1,x_n}，首先将该等时间间隔采样获得的时域信号，转化为等时间间隔采样的角域信号，该计算过程中假设转速短时间内按照等角加速度的方式变速，采用等时间间隔采样时振动频率最低的振动周期所包含的采样点数越多，为保证重构信号不会失真，以最低振动频率周期包含的点数作为插值的标准值参考值，将该参考值转化为等角度采样时旋转一周所包含的点数；该点数即为等角度采样时转子每转应采集额点数，转速越快每转采集的点数越少，因此转速参考应以最高转速为参考；假设振动测点的最低频率成分为f₀，齿轮箱传动比为n(三个传动级的传动比分别为n₁,n₂,n₃)，按照发电机侧的最高转速要求，正常发电时主轴的转速范围为n_min～n_max，则以高速轴转速作为参考转速，其值为n_ck＝n₁×n₂×n₃×n_max，计算得到阶比重采样的角域信号重构的插值参考点数为n＝n_ck*f_o/60，即表示采用等角度采样方式时每转所采集的点数应为n个，假设振动信号分析周期内包含的总点数为N，则重新定义的角域序列为：Theta(N)＝0:2π/n:2π(N-1)/n；采用线性插值方法将Theta(N)插入等时间间隔采样下的角域信号中，获得等角度采样方式下的振动角域信号；

②将等时间间隔采集的振动时域信号转化为相应的角域振动信号，由于风电机组变桨系统的功能作用，风电机组转速变动相对平缓，即可假设风电机组转速为均角加速度变速，其中机组转速信号可从发电机转速编码器获取，因此假设每秒钟的时间内的机组转速为均角加速度变化的；设振动测点某十分钟内获得的主轴转速数据为n_a，下十分钟获得的主轴转速数据为n_b，对应的角速度分别为ω_b和ω_a，因此该时间段内的机组传动链振动测点参考转速变化曲线表示为：

(ω_b-ω_a)·t+ω_a(1)

对该角速度曲线在时间尺度上进行积分获得角度公式：

$\mathrm{\θ}={\∫}_{t_1}^{t_2}\[\left({\mathrm{\ω}}_b-{\mathrm{\ω}}_a\right)\·t+{\mathrm{\ω}}_a\]dt---\left(2\right)$

通过式(1)将等时间间隔采样获得时域振动信号转化为对应的角域信号；

③采用线性插值的方法对等时间间隔的角域信号进行插值，其每转动一圈所采集的点数在步骤①中已经计算获得；对等时间间隔采集角域信号的角度按照插值点数进行划分，求得需要进行插值的点的角度坐标{θ₁,θ₂,…,θ_n-1,θ_n}，对这些点在等时间间隔采样角域信号上进行线性插值，获得{x_θ1,x_θ2,…x_θn-1,x_θn}，通过公式[x_t(i-1)-x_θ(i)]×[x_t(i)-x_θ(i)]≤0搜索位于插值点两侧的真实值，通过线性插值公式

$x_{\mathrm{\θ}\left(i\right)}=\frac{x_{t\left(i\right)}\·\left(\mathrm{\θ}\right(i)-t(i-1))}{t\left(i\right)-t\left(i-1\right)}---\left(3\right)$

计算线性插值后插入坐标点的纵坐标，获得等角度角域信号{x_θ1,x_θ2,…x_θn-1,x_θn}；至此我们通过步骤2中的1-3步实现了等时间间隔采样时域信号向等角度采样角域信号的转化，转化后的角域信号具有平稳特性和整周期特性；

至此已经通过阶比重采样技术实现了变工况条件下的非平稳时域振动信号向平稳角域振动信号的转化，此时的振动信号已经具有整周期特性，对具有整周期特性的角域振动信号提取故障预警和诊断指标能够准确、有效地反映故障发生与否或发展程度；但要注意的是，角域采样中，与等时间采样频率S_f对应的为阶次采样频率S_o(SamplingfrequencyofOrdertracking)，即参考轴每旋转一周所采集的等角度数据点数，为了保证原始信号中信息的完整性，阶次采样与时域采样一样，都需要满足采样定理：阶次采样频率必须大于最高阶次成分O_max的两倍，即

S_o＞O_max(4)；

对重采样的等角度信号作傅里叶变换后得到的谱图为阶次谱；利用等角度采样信号和阶次谱进行角域信号的分析和信号的时域、频域分析中的各参量具有一一对应的关系；

重采样后的振动角域信号常常包含噪声，为了降低此成分的干扰，本文选用基于EMD方法进行降噪处理，此方法可以剔除噪声，降低低频干扰有效突出高频固有振动，易于在复杂信号中提取故障特征频率；EMD分解是一种自适应的信号分解方法，其优点表现为：(1)基函数自动产生：从EMD分解过程可以得出，该方法根据原始信号的特点，自适应的选取最优的基函数，而用小波分解需要预先选择基函数，首先基函数的选择比较麻烦，分解过程中基函数一旦选定也不能改变；(2)具有自适应的滤波特性；(3)自适应的多分辨率；EMD分解后的IMF包含了从高到低的频率成分；由于EMD分解过程中存在端点效应，插值误差，过度分解等情况，造成分解结果可能存在伪分量，和原始信号无关；伪分量中可能含有和故障频率重合的频率，应该想办法将这些造成影响的伪分量剔除掉；通过互相关系数准则可以有效识别伪分量，实现方法是计算IMF分量与原始信号的互相关系数，通过系数的大小可识别伪分量；

分量与原始信号S的互相关系数为：

${\mathrm{\ρ}}_{s,{\hat{c}}_j}=max\left(R_{s,{\hat{c}}_j}\right(\mathrm{\τ}))/max\left(R_s\right(\mathrm{\τ}\left)\right)---(5);$

其中，R_s(τ)为原信号的自相关系数；

通过互相关系数准则可剔除伪分量，然后进行重构信号，降低噪声的干扰；经过阶比重采样以及EMD降噪处理后的信号为角域平稳信号。

步骤(2)所述的对行星轮系的分布故障和局部故障的特征频率的计算进行分析：

风电机组齿轮增速箱结构多样，传动比大，为减小齿轮箱的尺寸，一般为行星齿轮结构，这里对某一风电机组行星齿轮箱进行分析，其风电行星齿轮箱结构为两级行星轮、一级平行齿轮结构；

风电机组行星齿轮箱的两级行星轮系和平行级齿轮结构；

行星齿轮箱不同于传统的齿轮箱，其结构由太阳轮、行星轮、齿圈和行星架组成；一般情况下，齿圈固定不动，太阳轮绕自身的中心轴线旋转，行星轮不仅自转，还围绕太阳轮公转；行星轮既和太阳轮啮合，又和齿圈啮合；行星齿轮箱中多个齿轮的复合运动造成了振动信号的复杂性；在振动的监测中，传感器一般安装在齿圈或与之相连的箱体上采集振动信号，太阳轮-行星轮以及行星轮-齿圈啮合副的啮合点对传感器的位置随行星架旋转变化，使得啮合点至传感器之间的振动传递路径发生变化，这种时变的传递路径对振动测试信号产生调幅调制效应；

齿轮故障一般分为两类：分布式故障和局部故障；行星轮系发生不同故障类型和故障位置时，其振动信号模型和故障特征频率会有所不同，对故障发生时的振动信号进行分析可以实现对故障的诊断；

风电行星齿轮箱一般分为行星轮系、平行级齿轮，对于行星轮系和平行级齿轮，其故障可分为局部故障和分布式故障；在单级行星齿轮箱中，太阳轮–行星轮和行星轮–齿圈两种啮合副的啮合频率相同；通常齿圈固定不动，太阳轮、行星轮和行星架旋转，在这种情况下，啮合频率：

f_m＝f_cZ_r＝(f_s ^(r)-f_c)Z_s(6)；

式中：Z_r和Z_s分别为齿圈和太阳轮的齿数；f_m为啮合频率；f_c为行星架的旋转频率；f_s ^(r)为太阳轮的绝对旋转频率；

太星轮局部故障特征频率为：

$f_s=\frac{f_m}{Z_s}N---\left(7\right);$

式中：f_m为啮合频率；Z_s为太阳轮齿数；N为行星轮数量；

行星轮局部故障特征频率为：

$f_p=2\frac{f_m}{Z_p}---\left(8\right);$

式中：f_m为啮合频率；Z_p为行星轮齿数；齿圈局部故障特征频率为：

$f_r=\frac{f_m}{Z_r}N---\left(9\right);$

式中：f_m为啮合频率；Z_r为齿圈齿数；

行星齿轮箱中各种齿轮的分布式故障特征频率等于齿轮相对行星架(太阳轮和齿圈故障)或齿圈(行星轮故障)的旋转频率；已知行星齿轮箱的啮合频率f_m和某个齿轮的齿数Z_g，则该齿轮相对行星架(太阳轮和齿圈故障)或齿圈(行星轮故障)的旋转频率：

f_g＝f_g/Z_g(10)；

则太阳轮、行星轮和齿圈分布式故障的特征频率分别为：

f_s'＝f_m/Z_s(11)；

f_p'＝f_m/Z_p(12)；

f_r'＝f_m/Z_r(13)；

式中：f_m为啮合频率；f_s'、f_p'、f_r'太阳轮、行星轮和齿圈分布式故障的特征频率；Z_s为太阳轮齿数；Z_p为行星轮的齿数；Z_r为齿圈齿数；

以与主轴相连接的一级行星轮系的行星架为参考转速，对行星齿轮箱中各级的各个齿轮的局部故障和分布故障的特征阶比进行计算；

这里的信号为matlab仿真信号，信号中忽略齿轮箱中轴承以及各齿轮之间的影响，假设各齿轮箱传动级之间的振动影响不存在，对该方法进行仿真验证针对行星轮系处于正常、分布故障、局部故障时的振动信号模型参考文献，这里就不再赘述；

由此，仿真行星齿轮箱各级处于正常工况时，其振动信号模型为:

x₁(t)＝A·[1-cos(2π·3f_c1·t)]cos(2π·f_m1·t+θ₁)

+B·[1-cos(2π·3f_c2·t)]·cos(2π·f_m2·t+θ₂)+C·cos(2π·f_m3·t+θ₃)(14)

行星齿轮箱一级行星轮系太阳轮发生局部故障时，其振动信号模型为:

行星齿轮箱一级行星轮系太阳轮发生分布故障时，其振动信号模型为:

式(9)～(11)中：x₁(t)、x₂(t)、x₃(t)为行星齿轮箱处于正常一级行星轮系太阳轮发生局部故障、分布故障时的振动信号序列；t为时间序列；θ₁、θ₂、θ₃、φ、为初始相位；f_m1、f_m2、f_m3为各级啮合频率；f_c1、f_c2为一、二级行星架的旋转频率；为一级太阳轮的绝对旋转频率；f_s1、f_s1'为一级太阳轮发生局部故障和分布故障时的特征频率；A、B、C为无量纲常数，行星齿轮箱各个状态时值会有所不同，这里就不再详述；各个振动信号采用频率为8192HZ。

步骤(3)所述的故障特征参数：在行星齿轮箱发生故障时往往使其频率组成及各成分幅、频特性发生变化，从而使时域波形发生不同程度的畸变，这种畸变也是行星齿轮箱齿轮故障在时域信号中的直接反映；信号的波形、频谱结构发生变化时，其信号的复杂程度也发生变化；姜建东用Lempel-Ziv复杂性测度定量评估大机组运行状态；Hong研究了Lempel-Ziv复杂度指标用于评估故障轴承损伤状态；实践表明，这一指标是衡量有限时间序列复杂程度的高效工具；

序列的Lempel-Ziv复杂度可按经过次循环得到，计算过程如下：

(1)初始化S_v,0＝{}，Q₀＝{}，C_N(0)＝0，r＝1；令由于Q_r不属于S_v,r-1,则C_N(r)＝C_N(r-1)，Q_r＝{}，r＝r+1；

(2)令判断Q_r是否属于若属于，则C_N(r)＝C_N(r-1)，r＝r+1，重复过程(2)；

(3)若不属于，则C_N(r)＝C_N(r-1)+1，Q_r＝{}，r＝r+1，重复过程(2)；

例如序列{0000···}的分段形式为{0·000···}，C_N＝2；符号序列{010101···}的分段形式为{0·1·0101···}，C_N就为3；出于对数据严格的可比性考虑，Lempel和Ziv进一步提出了一种归一化基准公式，将值界定在[0-1]之间，详细算法为：

$0\≤C_{nN}=\frac{C_N\left(N\right)}{C_{UL,N}}\≤1---\left(12\right)$

$C_{UL,N}=\underset{N\→\mathrm{\∞}}{\lim }{C}_N\left(N\right)=\underset{N\→\mathrm{\∞}}{\lim }\frac{N}{(1-\mathrm{\β}){\log }_{k}N}\≈\frac{N}{{\log }_{k}N}---\left(13\right)$

其中：k为中S_N元素的个数(对于二进制S_N序列，k＝2)；根据Lempel-Ziv复杂度算法和文献可知，一个序列的复杂度越大，说明它的添加操作越多，序列的内涵分量也越多，描述给定符号序列所需最少的、互不相同的分段形式越多，给定的序列周期性越弱，出现新内涵模式的速度也就越快；Lempel-Ziv复杂度的物理意义在于它反映了一个时间序列随着长度的增加出现新模式的速率；复杂度越大，说明数据在窗口长度时期内随时间出现的新变化越多，发生新变化的速度越快，表明这一时期的数据变化是无序而复杂的；反之复杂度越小，则说明发生新变化的速度越慢，数据变化是规则的，周期性越强；因此，振动信号的Lempel-Ziv复杂度指标能够客观反映出系统状态发生变化的情况。

步骤(4)所述的实现行星齿轮箱故障模式的识别：是对行星齿轮箱正常状态下的变转速振动信号进行分析，振动信号长度为20s，其振动信号为变转速信号，振动信号的幅值与频率随转速的减小而有所降低；根据变转速振动信号经过阶比重采样和EMD分解方法处理后的角域平稳信号，得出，角域信号相较于时域信号更加平稳，信号的采样间隔更加规律；

对行星齿轮箱的振动时域信号与角域信号的复杂度进行计算；时域信号进行计算的长度分别为3s和1.5s的信号长度，对角域信号进行计算的长度分别为和的信号长度，通过计算不同信号长度的Lempel-Ziv复杂度数值可以看出，时域信号的复杂度数值波动较大，不够平稳，而角域信号的复杂度数值趋于一致，更加稳定，说明，角域信号的复杂度相较于时域信号更能反映行星齿轮箱的运行状态；

借助阶比分析方法对行星齿轮箱不同运行状态下的平稳角域信号进行Lempel-Ziv复杂度的计算，从左往右依次对应行星齿轮箱正常、齿圈分布故障、行星轮分布故障、太阳轮分布故障、齿圈局部故障、行星轮局部故障、太阳轮局部故障的Lempel-Ziv复杂度变化趋势，得出具体数值，

从结果中可以看出，行星齿轮箱正常状态的振动最为简单，局部故障发生时其非线性相较更大，相较于分布故障具有更高的Lempel-Ziv值；通过Lempel-Ziv复杂度数值，能够对行星齿轮箱的运行状态进行有效的识别。

步骤(5)所述的实例验证发现：对某风电场华锐SL1500风电齿轮箱振动数据进行分析，振动信号的采用频率为5120Hz，信号长度为6s，三段信号分别为齿轮箱正常状态、齿轮箱高速级齿轮磨损状态和齿轮箱高速级齿轮断齿状态下的振动信号，各段振动信号处于不同转速，时域信号不够平稳，通过对振动时域信号进行阶比重采样和EMD分解方法的预处理，计算三个状态下的信号的Lempel-Ziv复杂度；

由于缺少齿轮箱各个运行状态下的运行数据，难以对齿轮箱各个故障状态进行识别，但从上面三个运行状态下的Lempel-Ziv数值仍可以看出，Lempel-Ziv复杂度可以作为表征齿轮箱运行状态的特征参数；同时在行星齿轮箱故障数据充足的情况下，可以对行星齿轮的故障进行划分，实现对每一故障模式下的Lempel-Ziv复杂度的区间的划分。

本发明的优点及有益效果如下：

本发明采用阶比分析方法可以避免变工况的影响，经过阶比分析处理后的角域信号的Lempel-Ziv复杂度较时域信号更加平稳，Lempel-Ziv复杂度指标用于评价不同故障状态动力系统的变化，它是一种表征时间序列里出现新模式的速率的方法，这个方法最先由Lempel和Ziv提出，因此取名为Lempel-Ziv复杂度，其受变工况的影响较小，可以避免变工况的影响。行星齿轮箱振动角域信号的Lempel-Ziv复杂度能够对行星齿轮局部故障和分布故障进行识别，通过仿真和实际案例的分析，验证了本发明的有效性，同时在缺少故障数据的情况下也可以作为故障发生与否的特征参数。

附图说明

图1振动数据阶比重采样流程图

图2为本发明中行星齿轮箱结构简略图；

图3为本发明中行星齿轮箱转速信号和振动时域信号；

图4为本发明中行星齿轮箱振动角域信号；

图5为本发明中正常状态下时域信号与角域信号的Lempel-Ziv复杂度数值；

图6为本发明中不同状态的Lempel-Ziv数值变化趋势；

图7为本发明中行星齿轮箱三个运行状态的振动时域信号；

图8为本发明中基于Lempel-Ziv复杂度的行星齿轮箱故障诊断流程。

具体实施方式

本发明是一种基于阶比复杂度的行星齿轮箱故障诊断方法，包括如下步骤：

(1)振动信号预处理。根据阶比跟踪技术，将处于变工况的风电行星齿轮箱传感器所采集的振动信号进行预处理，将非线性、非平稳的时域信号转化为具有平稳性的角域信号，避免了使用硬件方式实现等角度采样的昂贵成本。基于线性插值方法的非平稳振动时域信号的阶比重构技术，将等时间间隔采样的非平稳振动时域信号转化为具有平稳特性的角域振动信号，保证了行星齿轮箱振动角域信号的整周期性，避免了变工况运行环境的影响。

(2)将行星轮系的故障分为两类：分布故障和局部故障。对行星轮系的分布故障和局部故障的特征频率的计算进行分析。以华锐SL1500风电行星齿轮箱为研究对象，对两级行星轮系、一级平行级行星齿轮箱的振动信号进行仿真。

(3)故障特征参数。Lempel-Ziv复杂度在滚动轴承和简单平行级齿轮的故障诊断中已经表明，可以作为轴承故障诊断的故障特征参数，对Lempel-Ziv复杂度的概念进行分析，进一步确认Lempel-Ziv复杂度可以作为故障诊断的特征参数。

(4)对行星行星齿轮箱变工况下的振动信号进行分析，发现非平稳振动信号的Lempel-Ziv复杂度会有所波动，但经过阶比重采样处理后的角域信号可以有效避免了变工况对Lempel-Ziv复杂度的影响。由此可以验证，可以作为变工况下故障诊断特征参数。Lempel-Ziv复杂度可以作为恶劣工况环境下的行星齿轮箱的状态识别的特征参数。在步骤(1)、(2)的基础上，采用阶比复杂度对复杂结构的行星齿轮箱故障进行诊断，计算不同故障模式下的阶比复杂度，通过观察各故障所对应的阶比复杂度的变化情况，发现阶比复杂度可以作为行星齿轮箱故障模式识别的特征参数。该步骤利用阶比复度，通过对风电行星齿轮箱处于不同状态下的故障特征阶比复杂度进行计算，可以实现行星齿轮箱故障模式的识别。

(5)实例验证发现，行星齿轮箱在不同故障模式下的阶比复杂度相差较大，可以作为行星齿轮箱状态的评价，在缺少充足故障数据的情况下，也可以作为齿轮箱健康状态的评估。经该步骤实例验证，利用阶比复度，可以对风电行星齿轮箱故障状态进行识别，在缺少故障数据的情况下，也可以作为行星齿轮箱健康状态的评估参数。

下面针对本发明方法作详细的解释如下：

步骤1，齿行星轮箱振动信号预处理。行星齿轮箱处于变工况的工作环境下，其采集的振动信号为非平稳信号，这对齿轮箱的故障诊断产生很大的困难，采用阶比重采样的方法对振动信号进行角域重采样，将非平稳的时域振动信号转化为平稳或准平稳特性的角域信号，便于对振动信号的分析。阶比分析技术的核心在于获得相对参考轴的恒定角增量采样数据，因此需要能准确获得阶次采样的时刻及相应的基准转速，即实现阶次跟踪。常见的阶比跟踪方法有硬件阶次跟踪法、计算阶比跟踪法和基于瞬时频率估计的阶比跟踪法等。采用计算阶比跟踪法实现振动信号的重采样计算，通过重采样可将机械变速过程中产生的与转速有关的振动信号有效的分离出来，同时对与转速无关的信号起到一定的抑制作用，按转角位置重新分配信号的采样间隔，可以剔除转速变化对信号随机性的影响。实际的齿轮箱振动信号中一般情况下都含有多种干扰成分，这就使得其故障特征的提取变得比较困难。EMD方法可以根据信号的局部时变特性，自适应的将任意一个复杂信号分解为一系列分量，通过相关系数法则对信号进行重构，剔除原始信号中的干扰成分。采用阶比分析与EMD分解方法可以有效的对振动信号进行预处理，对下一步的信号分析做好了准备。

由于运行工况的复杂性和部件间的耦合性，风电机组的振动数据具有明显的非线性和非平稳特性，因此必须采用振动数据预处理技术对机组振动数据进行前置处理。本发明基于阶比重采样技术的风电机组振动数据处理方法，该方法可以有效地解决风电机组振动信号的变转速问题，可以将具有非平稳性特性的振动时域信号转化为具有平稳特性的振动角域信号，其基本原理是将采用软件算法将等时间间隔采样的时域信号转化为等角度采样模式下的角域振动信号，转化后的角域振动信号即保证了原始时域信号的真实性，有具有了正周期特性，在角域振动信号的基础上提取的振动特征提取能够准确有效地反映风电机组变工况运行条件下的故障状态。

其阶比重采样技术的主要步骤如下所示：

①假设振动测点获得的振动数据为{x₁,x₂,x₃,…x_n-2,x_n-1,x_n}，首先将该等时间间隔采样获得的时域信号，转化为等时间间隔采样的角域信号，该计算过程中假设转速短时间内按照等角加速度的方式变速，采用等时间间隔采样时振动频率最低的振动周期所包含的采样点数越多，为保证重构信号不会失真，以最低振动频率周期包含的点数作为插值的标准值参考值，将该参考值转化为等角度采样时旋转一周所包含的点数。该点数即为等角度采样时转子每转应采集额点数，转速越快每转采集的点数越少，因此转速参考应以最高转速为参考。假设振动测点的最低频率成分为f₀，齿轮箱传动比为n(三个传动级的传动比分别为n₁,n₂,n₃)，按照发电机侧的最高转速要求，正常发电时主轴的转速范围为n_min～n_max，则以高速轴转速作为参考转速，其值为n_ck＝n₁×n₂×n₃×n_max，计算得到阶比重采样的角域信号重构的插值参考点数为n＝n_ck*f_o/60，即表示采用等角度采样方式时每转所采集的点数应为n个，假设振动信号分析周期内包含的总点数为N，则重新定义的角域序列为：Theta(N)＝0:2π/n:2π(N-1)/n；采用线性插值方法将Theta(N)插入等时间间隔采样下的角域信号中，获得等角度采样方式下的振动角域信号。

②将等时间间隔采集的振动时域信号转化为相应的角域振动信号，由于风电机组变桨系统的功能作用，风电机组转速变动相对平缓，即可假设风电机组转速为均角加速度变速，其中机组转速信号可从发电机转速编码器获取，因此假设每秒钟的时间内的机组转速为均角加速度变化的。设振动测点某十分钟内获得的主轴转速数据为n_a，下十分钟获得的主轴转速数据为n_b，对应的角速度分别为ω_b和ω_a，因此该时间段内的机组传动链振动测点参考转速变化曲线表示为：

(ω_b-ω_a)·t+ω_a(1)

对该角速度曲线在时间尺度上进行积分获得角度公式：

$\mathrm{\θ}={\∫}_{t_1}^{t_2}\[\left({\mathrm{\ω}}_b-{\mathrm{\ω}}_a\right)\·t+{\mathrm{\ω}}_a\]dt---\left(2\right)$

通过式(1)将等时间间隔采样获得时域振动信号转化为对应的角域信号。

③采用线性插值的方法对等时间间隔的角域信号进行插值，其每转动一圈所采集的点数在步骤①中已经计算获得。对等时间间隔采集角域信号的角度按照插值点数进行划分，求得需要进行插值的点的角度坐标{θ₁,θ₂,…,θ_n-1,θ_n}，对这些点在等时间间隔采样角域信号上进行线性插值，获得{x_θ1,x_θ2,…x_θn-1,x_θn}，通过公式[x_t(i-1)-x_θ(i)]×[x_t(i)-x_θ(i)]≤0搜索位于插值点两侧的真实值，通过线性插值公式

$x_{\mathrm{\θ}\left(i\right)}=\frac{x_{t\left(i\right)}\·\left(\mathrm{\θ}\right(i)-t(i-1))}{t\left(i\right)-t\left(i-1\right)}---\left(3\right)$

计算线性插值后插入坐标点的纵坐标，获得等角度角域信号{x_θ1,x_θ2,…x_θn-1,x_θn}。至此我们通过步骤2中的1-3步实现了等时间间隔采样时域信号向等角度采样角域信号的转化，转化后的角域信号具有平稳特性和整周期特性，阶比重采样角域信号重构的流程见图1。

至此已经通过阶比重采样技术实现了变工况条件下的非平稳时域振动信号向平稳角域振动信号的转化，此时的振动信号已经具有整周期特性，对具有整周期特性的角域振动信号提取故障预警和诊断指标能够准确、有效地反映故障发生与否或发展程度。但要注意的是，角域采样中，与等时间采样频率S_f对应的为阶次采样频率S_o(SamplingfrequencyofOrdertracking)，即参考轴每旋转一周所采集的等角度数据点数，为了保证原始信号中信息的完整性，阶次采样与时域采样一样，都需要满足采样定理：阶次采样频率必须大于最高阶次成分O_max的两倍，即

S_o＞O_max(4)

对重采样的等角度信号作傅里叶变换后得到的谱图为阶次谱。利用等角度采样信号和阶次谱进行角域信号的分析和信号的时域、频域分析中的各参量具有一一对应的关系，如表2所示。

重采样后的振动角域信号常常包含噪声，为了降低此成分的干扰，本文选用基于EMD方法进行降噪处理，此方法可以剔除噪声，降低低频干扰有效突出高频固有振动，易于在复杂信号中提取故障特征频率。EMD分解是一种自适应的信号分解方法，其优点表现为：(1)基函数自动产生：从EMD分解过程可以得出，该方法根据原始信号的特点，自适应的选取最优的基函数，而用小波分解需要预先选择基函数，首先基函数的选择比较麻烦，分解过程中基函数一旦选定也不能改变；(2)具有自适应的滤波特性；(3)自适应的多分辨率。EMD分解后的IMF包含了从高到低的频率成分。由于EMD分解过程中存在端点效应，插值误差，过度分解等情况，造成分解结果可能存在伪分量，和原始信号无关。伪分量中可能含有和故障频率重合的频率，应该想办法将这些造成影响的伪分量剔除掉。通过互相关系数准则可以有效识别伪分量，实现方法是计算IMF分量与原始信号的互相关系数，通过系数的大小可识别伪分量。

分量与原始信号S的互相关系数为：

${\mathrm{\ρ}}_{s,{\hat{c}}_j}=max\left(R_{s,{\hat{c}}_j}\right(\mathrm{\τ}))/max\left(R_s\right(\mathrm{\τ}\left)\right)---\left(5\right)$

其中，R_s(τ)为原信号的自相关系数。

通过互相关系数准则可剔除伪分量，然后进行重构信号，降低噪声的干扰。经过阶比重采样以及EMD降噪处理后的信号为角域平稳信号。

步骤2，故障特征频率计算。风电机组齿轮增速箱结构多样，传动比大，为减小齿轮箱的尺寸，一般为行星齿轮结构，本发明对某一风电机组行星齿轮箱进行分析，其风电行星齿轮箱结构为两级行星轮、一级平行齿轮结构，其结构如图2所示。

风电机组行星齿轮箱的两级行星轮系和平行级齿轮结构，其结构参数如表2所示。

行星齿轮箱不同于传统的齿轮箱，其结构由太阳轮、行星轮、齿圈和行星架组成。一般情况下，齿圈固定不动，太阳轮绕自身的中心轴线旋转，行星轮不仅自转，还围绕太阳轮公转。行星轮既和太阳轮啮合，又和齿圈啮合。行星齿轮箱中多个齿轮的复合运动造成了振动信号的复杂性。在振动的监测中，传感器一般安装在齿圈或与之相连的箱体上采集振动信号，太阳轮-行星轮以及行星轮-齿圈啮合副的啮合点对传感器的位置随行星架旋转变化，使得啮合点至传感器之间的振动传递路径发生变化，这种时变的传递路径对振动测试信号产生调幅调制效应。

齿轮故障一般分为两类：分布式故障和局部故障。行星轮系发生不同故障类型和故障位置时，其振动信号模型和故障特征频率会有所不同，对故障发生时的振动信号进行分析可以实现对故障的诊断，行星轮系分布故障和局部故障特征频率与振动信号模型在参考文献中有详细讲述，这里就不再赘述。

风电行星齿轮箱一般分为行星轮系、平行级齿轮，对于行星轮系和平行级齿轮，其故障可分为局部故障和分布式故障。在单级行星齿轮箱中，太阳轮–行星轮和行星轮–齿圈两种啮合副的啮合频率相同。通常齿圈固定不动，太阳轮、行星轮和行星架旋转，在这种情况下，啮合频率：

f_m＝f_cZ_r＝(f_s ^(r)-f_c)Z_s(6)；

式中：Z_r和Z_s分别为齿圈和太阳轮的齿数；f_m为啮合频率；f_c为行星架的旋转频率；f_s ^(r)为太阳轮的绝对旋转频率。

太阳轮局部故障特征频率为：

$f_s=\frac{f_m}{Z_s}N---\left(7\right);$

式中：f_m为啮合频率；Z_s为太阳轮齿数；N为行星轮数量。

行星轮局部故障特征频率为：

$f_p=2\frac{f_m}{Z_p}---\left(8\right);$

式中：f_m为啮合频率；Z_p为行星轮齿数。齿圈局部故障特征频率为：

$f_r=\frac{f_m}{Z_r}N---\left(9\right);$

式中：f_m为啮合频率；Z_r为齿圈齿数。

行星齿轮箱中各种齿轮的分布式故障特征频率等于齿轮相对行星架(太阳轮和齿圈故障)或齿圈(行星轮故障)的旋转频率。已知行星齿轮箱的啮合频率f_m和某个齿轮的齿数Z_g，则该齿轮相对行星架(太阳轮和齿圈故障)或齿圈(行星轮故障)的旋转频率：

f_g＝f_m/Z_g(10)；

则太阳轮、行星轮和齿圈分布式故障的特征频率分别为：

f_s'＝f_m/Z_s(11)；

f_p'＝f_m/Z_p(12)；

f_r'＝f_m/Z_r(13)；

式中：f_m为啮合频率；f_s'、f_p'、f_r'太阳轮、行星轮和齿圈分布式故障的特征频率；Z_s为太阳轮齿数；Z_p为行星轮的齿数；Z_r为齿圈齿数。

以与主轴相连接的一级行星轮系的行星架为参考转速，对行星齿轮箱中各级的各个齿轮的局部故障和分布故障的特征阶比进行计算，如表3所示。

本发明的信号为matlab仿真信号，信号中忽略齿轮箱中轴承以及各齿轮之间的影响，假设各齿轮箱传动级之间的振动影响不存在，对本发明的方法进行仿真验证针对行星轮系处于正常、分布故障、局部故障时的振动信号模型参考文献，这里就不再赘述。

由此，可以仿真行星齿轮箱各级处于正常工况时，其振动信号模型为:

x₁(t)＝A·[1-cos(2π·3f_c1·t)]cos(2π·f_m1·t+θ₁)

+B·[1-cos(2π·3_fc2·t)]·cos(2π·_fm2·t+θ₂)+C·cos(2π·_fm3·t+θ₃)(14)

行星齿轮箱一级行星轮系太阳轮发生局部故障时，其振动信号模型为:

行星齿轮箱一级行星轮系太阳轮发生分布故障时，其振动信号模型为:

式(14)～(16)中：x₁(t)、x₂(t)、x₃(t)为行星齿轮箱处于正常一级行星轮系太阳轮发生局部故障、分布故障时的振动信号序列；t为时间序列；θ₁、θ₂、θ₃、φ、为初始相位；f_m1、f_m2、f_m3为各级啮合频率；f_c1、f_c2为一、二级行星架的旋转频率；为一级太阳轮的绝对旋转频率；f_s1、f_s1'为一级太阳轮发生局部故障和分布故障时的特征频率；A、B、C为无量纲常数，行星齿轮箱各个状态时值会有所不同，这里就不再详述。各个振动信号采用频率为8192HZ。

步骤3，故障特征参数。行星齿轮箱发生故障时往往使其频率组成及各成分幅、频特性发生变化，从而使时域波形发生不同程度的畸变，这种畸变也是行星齿轮箱齿轮故障在时域信号中的直接反映。信号的波形、频谱结构发生变化时，其信号的复杂程度也发生变化。姜建东用Lempel-Ziv复杂性测度定量评估大机组运行状态。Hong研究了Lempel-Ziv复杂度指标用于评估故障轴承损伤状态。实践表明，这一指标是衡量有限时间序列复杂程度的高效工具。

序列的Lempel-Ziv复杂度可按经过次循环得到，计算过程如下：

(1)初始化S_v,0＝{}，Q₀＝{}，C_N(0)＝0，r＝1；令由于Q_r不属于S_v,r-1,则C_N(r)＝C_N(r-1)，Q_r＝{}，r＝r+1；

(2)令判断Q_r是否属于若属于，则C_N(r)＝C_N(r-1)，r＝r+1，重复过程(2)；

(3)若不属于，则C_N(r)＝C_N(r-1)+1，Q_r＝{}，r＝r+1，重复过程(2)。

例如序列{0000···}的分段形式为{0·000···}，C_N＝2。符号序列{010101···}的分段形式为{0·1·0101···}，C_N就为3。出于对数据严格的可比性考虑，Lempel和Ziv进一步提出了一种归一化基准公式，将值界定在[0-1]之间，详细算法为：

$0\≤C_{nN}=\frac{C_N\left(N\right)}{C_{UL,N}}\≤1---\left(17\right)$

$C_{UL,N}=\underset{N\→\mathrm{\∞}}{\lim }{C}_N\left(N\right)=\underset{N\→\mathrm{\∞}}{\lim }\frac{N}{(1-\mathrm{\β}){\log }_{k}N}\≈\frac{N}{{\log }_{k}N}---\left(18\right)$

其中：k为中S_N元素的个数(对于二进制S_N序列，k＝2)。根据Lempel-Ziv复杂度算法和文献可知，一个序列的复杂度越大，说明它的添加操作越多，序列的内涵分量也越多，描述给定符号序列所需最少的、互不相同的分段形式越多，给定的序列周期性越弱，出现新内涵模式的速度也就越快。Lempel-Ziv复杂度的物理意义在于它反映了一个时间序列随着长度的增加出现新模式的速率。复杂度越大，说明数据在窗口长度时期内随时间出现的新变化越多，发生新变化的速度越快，表明这一时期的数据变化是无序而复杂的。反之复杂度越小，则说明发生新变化的速度越慢，数据变化是规则的，周期性越强。因此，振动信号的Lempel-Ziv复杂度指标能够客观反映出系统状态发生变化的情况。

步骤4，故障模式的识别。对行星齿轮箱正常状态下的变转速振动信号进行分析，振动信号长度为20s，其转速信号与振动信号如图3所示，可以看出，振动信号为变转速信号，振动信号的幅值与频率随转速的减小而有所降低。图4为变转速振动信号经过阶比重采样和EMD分解方法处理后的角域平稳信号，从图中可以看出，角域信号相较于时域信号更加平稳，信号的采样间隔更加规律。

对行星齿轮箱的振动时域信号与角域信号的复杂度进行计算。时域信号进行计算的长度分别为3s和1.5s的信号长度，对角域信号进行计算的长度分别为和的信号长度，通过计算不同信号长度的Lempel-Ziv复杂度数值可以看出，时域信号的复杂度数值波动较大，不够平稳，而角域信号的复杂度数值趋于一致，更加稳定，可以说明，角域信号的复杂度相较于时域信号更能反映行星齿轮箱的运行状态，如图5所示。

借助阶比分析方法对行星齿轮箱不同运行状态下的平稳角域信号进行Lempel-Ziv复杂度的计算，图6从左往右依次对应行星齿轮箱正常、齿圈分布故障、行星轮分布故障、太阳轮分布故障、齿圈局部故障、行星轮局部故障、太阳轮局部故障的Lempel-Ziv复杂度变化趋势，表4为其具体数值。

从结果中可以看出，行星齿轮箱正常状态的振动最为简单，局部故障发生时其非线性相较更大，相较于分布故障具有更高的Lempel-Ziv值。通过Lempel-Ziv复杂度数值，可以对行星齿轮箱的运行状态进行有效的识别。

步骤5，实例验证。对某风电场华锐SL1500风电齿轮箱振动数据进行分析，振动信号的采用频率为5120Hz，信号长度为6s，三段信号分别为齿轮箱正常状态、齿轮箱高速级齿轮磨损状态和齿轮箱高速级齿轮断齿状态下的振动信号，各段振动信号处于不同转速，时域信号不够平稳，如图7所示。通过对振动时域信号进行阶比重采样和EMD分解方法的预处理，计算三个状态下的信号的Lempel-Ziv复杂度，其数值如表5所示。

由于缺少齿轮箱各个运行状态下的运行数据，难以对齿轮箱各个故障状态进行识别，但从上面三个运行状态下的Lempel-Ziv数值仍可以看出，Lempel-Ziv复杂度可以作为表征齿轮箱运行状态的特征参数。同时在行星齿轮箱故障数据充足的情况下，可以对行星齿轮的故障进行划分，实现对每一故障模式下的Lempel-Ziv复杂度的区间的划分。根据本发明，基于Lempel-Ziv复杂度的行星齿轮箱故障诊断流程如图8所示。

表1时域、角域对比关系。

表2风电行星齿轮箱结构参数。

表3行星齿轮箱故障特征频率计算。

表4行星齿轮箱不同状态的Lempel-Ziv复杂度。

表5行星齿轮箱三个状态下的Lempel-Ziv复杂度。

正常	高速级齿轮磨损	高速级齿轮断齿
			0.2628	0.4023	0.5625

Claims (6)

1.一种基于阶比复杂度的行星齿轮箱故障诊断方法，其特征在于：包括如下步骤：

(1)振动信号预处理：根据阶比跟踪技术，将处于变工况的风电行星齿轮箱传感器所采集的振动信号进行预处理，将非线性、非平稳的时域信号转化为具有平稳性的角域信号，避免了使用硬件方式实现等角度采样的昂贵成本；基于线性插值方法的非平稳振动时域信号的阶比重构技术，将等时间间隔采样的非平稳振动时域信号转化为具有平稳特性的角域振动信号，保证了行星齿轮箱振动角域信号的整周期性，避免了变工况运行环境的影响；

(2)将行星轮系的故障分为两类：分布故障和局部故障；对行星轮系的分布故障和局部故障的特征频率的计算进行分析；以华锐SL1500风电行星齿轮箱为研究对象，对两级行星轮系、一级平行级行星齿轮箱的振动信号进行仿真；

(3)故障特征参数：Lempel-Ziv复杂度在滚动轴承和简单平行级齿轮的故障诊断中已经表明，可以作为轴承故障诊断的故障特征参数，对Lempel-Ziv复杂度的概念进行分析，进一步确认Lempel-Ziv复杂度可以作为故障诊断的特征参数；

(4)对行星行星齿轮箱变工况下的振动信号进行分析，发现非平稳振动信号的Lempel-Ziv复杂度会有所波动，但经过阶比重采样处理后的角域信号可以有效避免了变工况对Lempel-Ziv复杂度的影响；由此可以验证，可以作为变工况下故障诊断特征参数；Lempel-Ziv复杂度可以作为恶劣工况环境下的行星齿轮箱的状态识别的特征参数；在步骤(1)、(2)的基础上，采用阶比复杂度对复杂结构的行星齿轮箱故障进行诊断，计算不同故障模式下的阶比复杂度，通过观察各故障所对应的阶比复杂度的变化情况，发现阶比复杂度可以作为行星齿轮箱故障模式识别的特征参数；该步骤利用阶比复度，通过对风电行星齿轮箱处于不同状态下的故障特征阶比复杂度进行计算，实现行星齿轮箱故障模式的识别；

(5)实例验证发现，行星齿轮箱在不同故障模式下的阶比复杂度相差较大，可以作为行星齿轮箱状态的评价，在缺少充足故障数据的情况下，也可以作为齿轮箱健康状态的评估；经该步骤实例验证，利用阶比复度，可以对风电行星齿轮箱故障状态进行识别，在缺少故障数据的情况下，也可以作为行星齿轮箱健康状态的评估参数。

2.根据权利要求1所述的一种基于阶比复杂度的行星齿轮箱故障诊断方法，其特征在于：步骤(1)所述的振动信号预处理：

是对齿行星轮箱振动信号预处理；行星齿轮箱处于变工况的工作环境下，其采集的振动信号为非平稳信号，这对齿轮箱的故障诊断产生很大的困难，采用阶比重采样的方法对振动信号进行角域重采样，将非平稳的时域振动信号转化为平稳或准平稳特性的角域信号，便于对振动信号的分析；阶比分析技术的核心在于获得相对参考轴的恒定角增量采样数据，因此需要能准确获得阶次采样的时刻及相应的基准转速，即实现阶次跟踪；常见的阶比跟踪方法有硬件阶次跟踪法、计算阶比跟踪法和基于瞬时频率估计的阶比跟踪法；采用计算阶比跟踪法实现振动信号的重采样计算，通过重采样可将机械变速过程中产生的与转速有关的振动信号有效的分离出来，同时对与转速无关的信号起到一定的抑制作用，按转角位置重新分配信号的采样间隔，可以剔除转速变化对信号随机性的影响；实际的齿轮箱振动信号中一般情况下都含有多种干扰成分，这就使得其故障特征的提取变得比较困难；EMD方法可以根据信号的局部时变特性，自适应的将任意一个复杂信号分解为一系列分量，通过相关系数法则对信号进行重构，剔除原始信号中的干扰成分；采用阶比分析与EMD分解方法可以有效的对振动信号进行预处理，对下一步的信号分析做好了准备；

由于运行工况的复杂性和部件间的耦合性，风电机组的振动数据具有明显的非线性和非平稳特性，因此必须采用振动数据预处理技术对机组振动数据进行前置处理；本发明基于阶比重采样技术的风电机组振动数据处理方法，该方法可以有效地解决风电机组振动信号的变转速问题，可以将具有非平稳性特性的振动时域信号转化为具有平稳特性的振动角域信号，其基本原理是将采用软件算法将等时间间隔采样的时域信号转化为等角度采样模式下的角域振动信号，转化后的角域振动信号即保证了原始时域信号的真实性，有具有了正周期特性，在角域振动信号的基础上提取的振动特征提取能够准确有效地反映风电机组变工况运行条件下的故障状态；

其阶比重采样技术的主要步骤如下所示：

①假设振动测点获得的振动数据为{x₁,x₂,x₃,…x_n-2,x_n-1,x_n}，首先将该等时间间隔采样获得的时域信号，转化为等时间间隔采样的角域信号，该计算过程中假设转速短时间内按照等角加速度的方式变速，采用等时间间隔采样时振动频率最低的振动周期所包含的采样点数越多，为保证重构信号不会失真，以最低振动频率周期包含的点数作为插值的标准值参考值，将该参考值转化为等角度采样时旋转一周所包含的点数；该点数即为等角度采样时转子每转应采集额点数，转速越快每转采集的点数越少，因此转速参考应以最高转速为参考；假设振动测点的最低频率成分为f₀，齿轮箱传动比为n(三个传动级的传动比分别为n₁,n₂,n₃)，按照发电机侧的最高转速要求，正常发电时主轴的转速范围为n_min～n_max，则以高速轴转速作为参考转速，其值为n_ck＝n₁×n₂×n₃×n_max，计算得到阶比重采样的角域信号重构的插值参考点数为n＝n_ck*f_o/60，即表示采用等角度采样方式时每转所采集的点数应为n个，假设振动信号分析周期内包含的总点数为N，则重新定义的角域序列为：Theta(N)＝0:2π/n:2π(N-1)/n；采用线性插值方法将Theta(N)插入等时间间隔采样下的角域信号中，获得等角度采样方式下的振动角域信号；

②将等时间间隔采集的振动时域信号转化为相应的角域振动信号，由于风电机组变桨系统的功能作用，风电机组转速变动相对平缓，即可假设风电机组转速为均角加速度变速，其中机组转速信号可从发电机转速编码器获取，因此假设每秒钟的时间内的机组转速为均角加速度变化的；设振动测点某十分钟内获得的主轴转速数据为n_a，下十分钟获得的主轴转速数据为n_b，对应的角速度分别为ω_b和ω_a，因此该时间段内的机组传动链振动测点参考转速变化曲线表示为：

(ω_b-ω_a)·t+ω_a(1)

对该角速度曲线在时间尺度上进行积分获得角度公式：

$\mathrm{\θ}={\∫}_{t_1}^{t_2}\[\left({\mathrm{\ω}}_b-{\mathrm{\ω}}_a\right)\·t+{\mathrm{\ω}}_a\]dt---\left(2\right)$

通过式(1)将等时间间隔采样获得时域振动信号转化为对应的角域信号；

③采用线性插值的方法对等时间间隔的角域信号进行插值，其每转动一圈所采集的点数在步骤①中已经计算获得；对等时间间隔采集角域信号的角度按照插值点数进行划分，求得需要进行插值的点的角度坐标{θ₁,θ₂,…,θ_n-1,θ_n}，对这些点在等时间间隔采样角域信号上进行线性插值，获得{x_θ1,x_θ2,…x_θn-1,x_θn}，通过公式[x_t(i-1)-x_θ(i)]×[x_t(i)-x_θ(i)]≤0搜索位于插值点两侧的真实值，通过线性插值公式

$x_{\mathrm{\θ}\left(i\right)}=\frac{x_{t\left(i\right)}\·\left(\mathrm{\θ}\right(i)-t(i-1))}{t\left(i\right)-t(i-1)}---\left(3\right)$

计算线性插值后插入坐标点的纵坐标，获得等角度角域信号{x_θ1,x_θ2,…x_θn-1,x_θn}；至此我们通过步骤2中的1-3步实现了等时间间隔采样时域信号向等角度采样角域信号的转化，转化后的角域信号具有平稳特性和整周期特性；

至此已经通过阶比重采样技术实现了变工况条件下的非平稳时域振动信号向平稳角域振动信号的转化，此时的振动信号已经具有整周期特性，对具有整周期特性的角域振动信号提取故障预警和诊断指标能够准确、有效地反映故障发生与否或发展程度；但要注意的是，角域采样中，与等时间采样频率S_f对应的为阶次采样频率S_o(SamplingfrequencyofOrdertracking)，即参考轴每旋转一周所采集的等角度数据点数，为了保证原始信号中信息的完整性，阶次采样与时域采样一样，都需要满足采样定理：阶次采样频率必须大于最高阶次成分O_max的两倍，即

S_o＞O_max(4)；

对重采样的等角度信号作傅里叶变换后得到的谱图为阶次谱；利用等角度采样信号和阶次谱进行角域信号的分析和信号的时域、频域分析中的各参量具有一一对应的关系；

重采样后的振动角域信号常常包含噪声，为了降低此成分的干扰，本文选用基于EMD方法进行降噪处理，此方法可以剔除噪声，降低低频干扰有效突出高频固有振动，易于在复杂信号中提取故障特征频率；EMD分解是一种自适应的信号分解方法，其优点表现为：(1)基函数自动产生：从EMD分解过程可以得出，该方法根据原始信号的特点，自适应的选取最优的基函数，而用小波分解需要预先选择基函数，首先基函数的选择比较麻烦，分解过程中基函数一旦选定也不能改变；(2)具有自适应的滤波特性；(3)自适应的多分辨率；EMD分解后的IMF包含了从高到低的频率成分；由于EMD分解过程中存在端点效应，插值误差，过度分解等情况，造成分解结果可能存在伪分量，和原始信号无关；伪分量中可能含有和故障频率重合的频率，应该想办法将这些造成影响的伪分量剔除掉；通过互相关系数准则可以有效识别伪分量，实现方法是计算IMF分量与原始信号的互相关系数，通过系数的大小可识别伪分量；

分量与原始信号S的互相关系数为：

${\mathrm{\ρ}}_{s,{\hat{c}}_j}=max\left(R_{s,{\hat{c}}_j}\right(\mathrm{\τ}\left)\right)/max\left(R_s\right(\mathrm{\τ}\left)\right)---\left(5\right);$

其中，R_s(τ)为原信号的自相关系数；

通过互相关系数准则可剔除伪分量，然后进行重构信号，降低噪声的干扰；经过阶比重采样以及EMD降噪处理后的信号为角域平稳信号。

3.根据权利要求1所述的一种基于阶比复杂度的行星齿轮箱故障诊断方法，其特征在于：步骤(2)所述的对行星轮系的分布故障和局部故障的特征频率的计算进行分析：

风电机组齿轮增速箱结构多样，传动比大，为减小齿轮箱的尺寸，一般为行星齿轮结构，这里对某一风电机组行星齿轮箱进行分析，其风电行星齿轮箱结构为两级行星轮、一级平行齿轮结构；

风电机组行星齿轮箱的两级行星轮系和平行级齿轮结构；

行星齿轮箱不同于传统的齿轮箱，其结构由太阳轮、行星轮、齿圈和行星架组成；一般情况下，齿圈固定不动，太阳轮绕自身的中心轴线旋转，行星轮不仅自转，还围绕太阳轮公转；行星轮既和太阳轮啮合，又和齿圈啮合；行星齿轮箱中多个齿轮的复合运动造成了振动信号的复杂性；在振动的监测中，传感器一般安装在齿圈或与之相连的箱体上采集振动信号，太阳轮-行星轮以及行星轮-齿圈啮合副的啮合点对传感器的位置随行星架旋转变化，使得啮合点至传感器之间的振动传递路径发生变化，这种时变的传递路径对振动测试信号产生调幅调制效应；

齿轮故障一般分为两类：分布式故障和局部故障；行星轮系发生不同故障类型和故障位置时，其振动信号模型和故障特征频率会有所不同，对故障发生时的振动信号进行分析可以实现对故障的诊断；

风电行星齿轮箱一般分为行星轮系、平行级齿轮，对于行星轮系和平行级齿轮，其故障可分为局部故障和分布式故障；在单级行星齿轮箱中，太阳轮–行星轮和行星轮–齿圈两种啮合副的啮合频率相同；通常齿圈固定不动，太阳轮、行星轮和行星架旋转，在这种情况下，啮合频率：

$f_m=f_c{Z}_r=(f_s^{\left(r\right)}-f_c)Z_s---\left(6\right);$

式中：Z_r和Z_s分别为齿圈和太阳轮的齿数；f_m为啮合频率；f_c为行星架的旋转频率；为太阳轮的绝对旋转频率；

太星轮局部故障特征频率为：

$f_s=\frac{f_m}{Z_s}N---\left(7\right);$

式中：f_m为啮合频率；Z_s为太阳轮齿数；N为行星轮数量；

行星轮局部故障特征频率为：

$f_p=2\frac{f_m}{Z_p}---\left(8\right);$

式中：f_m为啮合频率；Z_p为行星轮齿数；齿圈局部故障特征频率为：

$f_r=\frac{f_m}{Z_r}N---\left(9\right);$

式中：f_m为啮合频率；Z_r为齿圈齿数；

行星齿轮箱中各种齿轮的分布式故障特征频率等于齿轮相对行星架(太阳轮和齿圈故障)或齿圈(行星轮故障)的旋转频率；已知行星齿轮箱的啮合频率f_m和某个齿轮的齿数Z_g，则该齿轮相对行星架(太阳轮和齿圈故障)或齿圈(行星轮故障)的旋转频率：

f_g＝f_g/Z_g(10)；

则太阳轮、行星轮和齿圈分布式故障的特征频率分别为：

f_s'＝f_m/Z_s(11)；

f_p'＝f_m/Z_p(12)；

f_r'＝f_m/Z_r(13)；

式中：f_m为啮合频率；f_s'、f_p'、f_r'太阳轮、行星轮和齿圈分布式故障的特征频率；Z_s为太阳轮齿数；Z_p为行星轮的齿数；Z_r为齿圈齿数；

以与主轴相连接的一级行星轮系的行星架为参考转速，对行星齿轮箱中各级的各个齿轮的局部故障和分布故障的特征阶比进行计算；

这里的信号为matlab仿真信号，信号中忽略齿轮箱中轴承以及各齿轮之间的影响，假设各齿轮箱传动级之间的振动影响不存在，对该方法进行仿真验证针对行星轮系处于正常、分布故障、局部故障时的振动信号模型参考文献，这里就不再赘述；

由此，仿真行星齿轮箱各级处于正常工况时，其振动信号模型为:

x₁(t)＝A·[1-cos(2π·3f_c1·t)]cos(2π·f_m1·t+θ₁)

+B·[1-cos(2π·3f_c2·t)]·cos(2π·f_m2·t+θ₂)+C·cos(2π·f_m3·t+θ₃)(14)

行星齿轮箱一级行星轮系太阳轮发生局部故障时，其振动信号模型为:

行星齿轮箱一级行星轮系太阳轮发生分布故障时，其振动信号模型为:

式(9)～(11)中：x₁(t)、x₂(t)、x₃(t)为行星齿轮箱处于正常一级行星轮系太阳轮发生局部故障、分布故障时的振动信号序列；t为时间序列；θ₁、θ₂、θ₃、φ、为初始相位；f_m1、f_m2、f_m3为各级啮合频率；f_c1、f_c2为一、二级行星架的旋转频率；为一级太阳轮的绝对旋转频率；f_s1、f_s1'为一级太阳轮发生局部故障和分布故障时的特征频率；A、B、C为无量纲常数，行星齿轮箱各个状态时值会有所不同，这里就不再详述；各个振动信号采用频率为8192HZ。

4.根据权利要求1所述的一种基于阶比复杂度的行星齿轮箱故障诊断方法，其特征在于：步骤(3)所述的故障特征参数：

在行星齿轮箱发生故障时往往使其频率组成及各成分幅、频特性发生变化，从而使时域波形发生不同程度的畸变，这种畸变也是行星齿轮箱齿轮故障在时域信号中的直接反映；信号的波形、频谱结构发生变化时，其信号的复杂程度也发生变化；姜建东用Lempel-Ziv复杂性测度定量评估大机组运行状态；Hong研究了Lempel-Ziv复杂度指标用于评估故障轴承损伤状态；实践表明，这一指标是衡量有限时间序列复杂程度的高效工具；

序列的Lempel-Ziv复杂度可按经过次循环得到，计算过程如下：

(1)初始化S_v,0＝{}，Q₀＝{}，C_N(0)＝0，r＝1；令由于Q_r不属于S_v,r-1,则C_N(r)＝C_N(r-1)，Q_r＝{}，r＝r+1；

(2)令判断Q_r是否属于若属于，则C_N(r)＝C_N(r-1)，r＝r+1，重复过程(2)；

(3)若不属于，则C_N(r)＝C_N(r-1)+1，Q_r＝{}，r＝r+1，重复过程(2)；

例如序列{0000···}的分段形式为{0·000···}，C_N＝2；符号序列{010101···}的分段形式为{0·1·0101···}，C_N就为3；出于对数据严格的可比性考虑，Lempel和Ziv进一步提出了一种归一化基准公式，将值界定在[0-1]之间，详细算法为：

$0\≤C_{nN}=\frac{C_N\left(N\right)}{C_{UL,N}}\≤1---\left(12\right)$

$C_{UL,N}=\underset{N\→\mathrm{\∞}}{\lim }{C}_N\left(N\right)=\underset{N\→\mathrm{\∞}}{\lim }\frac{N}{(1-\mathrm{\β}){\log }_{k}N}\≈\frac{N}{{\log }_{k}N}---\left(13\right)$

其中：k为中S_N元素的个数(对于二进制S_N序列，k＝2)；根据Lempel-Ziv复杂度算法和文献可知，一个序列的复杂度越大，说明它的添加操作越多，序列的内涵分量也越多，描述给定符号序列所需最少的、互不相同的分段形式越多，给定的序列周期性越弱，出现新内涵模式的速度也就越快；Lempel-Ziv复杂度的物理意义在于它反映了一个时间序列随着长度的增加出现新模式的速率；复杂度越大，说明数据在窗口长度时期内随时间出现的新变化越多，发生新变化的速度越快，表明这一时期的数据变化是无序而复杂的；反之复杂度越小，则说明发生新变化的速度越慢，数据变化是规则的，周期性越强；因此，振动信号的Lempel-Ziv复杂度指标能够客观反映出系统状态发生变化的情况。

5.根据权利要求1所述的一种基于阶比复杂度的行星齿轮箱故障诊断方法，其特征在于：步骤(4)所述的实现行星齿轮箱故障模式的识别：

是对行星齿轮箱正常状态下的变转速振动信号进行分析，振动信号长度为20s，其振动信号为变转速信号，振动信号的幅值与频率随转速的减小而有所降低；根据变转速振动信号经过阶比重采样和EMD分解方法处理后的角域平稳信号，得出，角域信号相较于时域信号更加平稳，信号的采样间隔更加规律；

对行星齿轮箱的振动时域信号与角域信号的复杂度进行计算；时域信号进行计算的长度分别为3s和1.5s的信号长度，对角域信号进行计算的长度分别为和的信号长度，通过计算不同信号长度的Lempel-Ziv复杂度数值可以看出，时域信号的复杂度数值波动较大，不够平稳，而角域信号的复杂度数值趋于一致，更加稳定，说明，角域信号的复杂度相较于时域信号更能反映行星齿轮箱的运行状态；

借助阶比分析方法对行星齿轮箱不同运行状态下的平稳角域信号进行Lempel-Ziv复杂度的计算，从左往右依次对应行星齿轮箱正常、齿圈分布故障、行星轮分布故障、太阳轮分布故障、齿圈局部故障、行星轮局部故障、太阳轮局部故障的Lempel-Ziv复杂度变化趋势，得出具体数值，

从结果中可以看出，行星齿轮箱正常状态的振动最为简单，局部故障发生时其非线性相较更大，相较于分布故障具有更高的Lempel-Ziv值；通过Lempel-Ziv复杂度数值，能够对行星齿轮箱的运行状态进行有效的识别。

6.根据权利要求1所述的一种基于阶比复杂度的行星齿轮箱故障诊断方法，其特征在于：步骤(5)所述的实例验证发现：

对某风电场华锐SL1500风电齿轮箱振动数据进行分析，振动信号的采用频率为5120Hz，信号长度为6s，三段信号分别为齿轮箱正常状态、齿轮箱高速级齿轮磨损状态和齿轮箱高速级齿轮断齿状态下的振动信号，各段振动信号处于不同转速，时域信号不够平稳，通过对振动时域信号进行阶比重采样和EMD分解方法的预处理，计算三个状态下的信号的Lempel-Ziv复杂度；

由于缺少齿轮箱各个运行状态下的运行数据，难以对齿轮箱各个故障状态进行识别，但从上面三个运行状态下的Lempel-Ziv数值仍可以看出，Lempel-Ziv复杂度可以作为表征齿轮箱运行状态的特征参数；同时在行星齿轮箱故障数据充足的情况下，可以对行星齿轮的故障进行划分，实现对每一故障模式下的Lempel-Ziv复杂度的区间的划分。