张旭帮，王灵梅，孟恩隆，贾成真，原升耀，王 平

(1.山西大学自动化与软件学院，山西 太原 030013)(2.煤炭工业太原设计研究院，山西 太原 030001)

风电机组通常在复杂多变的恶劣环境中运行，传动链中的齿轮、轴承等部件容易出现故障，影响正常的能量传递，造成发电损失，严重时会损坏齿轮箱、发电机等关键部件，形成巨大的运维成本。因此，行之有效的故障诊断方法对于确保风电机组的可靠和安全运行至关重要。由于机组运行受环境噪声影响，信号表现为非平稳、非线性，而目前的时域、频域分析方法效果不佳，因此时频分析法被提出。Huang等[1]提出了经验模态分解(empirical mode decomposition，EMD)方法，在非平稳信号处理方面得到了广泛应用，但它依赖基函数的选取，在递归过程中误差会不断叠加，存在端点效应、虚假分量、模态混叠、缺乏理论支撑等问题；Dragomiretskiy等[2]提出了变分模态分解(variational mode decomposition，VMD)方法，其本质上是通过变分问题使分解结果稳定，虽然克服了模态混叠、端点效应等缺点，但VMD方法受预设参数模态分量数K与惩罚因子α的制约。

近年来，王朝阁等[3]提出利用粒子群算法对最优最小反褶积算法参数进行优化，确定滤波器长度与故障周期搜索范围，对轴承故障进行识别与诊断；夏俊等[4]利用松鼠优化算法优化共振稀疏分解的品质因子与权重系数，得到高低共振分量，经Hilbert包络分析后用于提取轴承微弱故障信息有一定的效果；王奉涛等[5]提出在VMD算法的不同K值下进行信号分解，通过计算各分量与信号总能量，选取最佳本征模态分量(intrinsic mode function，IMF)，避免了信号分解不足或过分解问题；李华等[6]为使VMD算法取得较好的分解效果，对比不同K值下分解所得到的IMF峭度值，确定最优K值，最后进行滤波、包络解调分析判断轴承故障。在上述研究中，仅考虑了模态分量数K值对VMD算法的影响，而忽略了惩罚因子的影响。李华等[7]提出了基于信息熵优化VMD参数的方法，但优化过程是顺序优化，计算不同K值下所得分量的信息熵，确定最小信息熵对应K值，通过既定K值去确定惩罚因子，因此所得参数组合不一定是全局最优解。

针对VMD优化参数难以确定的问题，本文提出了基于改进参数优化VMD的故障特征提取方法。考虑到旋转机械发生故障时所具有的强噪声、周期性特征，用对周期性冲击较敏感的峭度指标和衡量信号时间序列复杂度的样本熵构造综合性能指标，作为优化算法的适应度函数。将振动信号经参数优化的变分模态分解后，筛选综合性能指标较小的模态分量进行重构，此时所得到的重构信号既保留了原始信号的特征，也突出了故障周期冲击特征，最后进行Hilbert包络分析，判断故障部位。通过对山西某风电场的齿轮箱故障数据与风电机组传动系统故障试验台采集的实测故障数据进行对比验证，表明该方法能够在强噪声背景下准确提取故障特征信息。

1 基本原理

1.1 变分模态分解

VMD算法假设实际信号分解所得IMF是具有中心频率的有限带宽，通过循环迭代确定最优解，寻找K个模态分量，经过计算信号梯度平方的L-2范数来估计各分量的带宽，使每个模态的估计带宽之和最小，实现对信号的自适应分解。其对应的约束变分模型表达式为：

(1)

引入二次惩罚因子α和Lagrange乘子λ，利用增广拉格朗日乘子法将式(1)转化为无约束优化模型，表达式如下：

(2)

式中：L(·)表示利用增广拉格朗日乘数法构建的无约束方程；y(t)为式(1)中的所有模态之和；λ,λ(t)为拉格朗日乘子。

(3)

(4)

(5)

具体实现过程如下：

3)根据式(5)更新λ。

4)根据式(6)判断是否收敛，其中ε为任意大于0的正数，代表分解精度。若满足收敛条件则信号分解过程结束；否则迭代次数n=n+1,返回步骤2)继续分解。

(6)

1.2 样本熵

为了减小环境噪声对提取故障特征的影响，本文引入样本熵对所采集的实际信号进行评价。样本熵由非线性动力学研究人员在20世纪末提出[8]，可以衡量实际信号的随机性和复杂性。

假设实际信号是由N个数据组成的时间序列{x(n)}=x(1),x(2),…,x(N)，则序列样本熵SampEn(m,r)为：

(7)

式中：r为相似容限阈值；m为初始分段的维数；Bm(r)和Am+1(r)分别是由时间序列组成的m维与(m+1)维向量序列，经过相空间重构、重构向量距离计算，满足相似容限r下匹配m与(m+1)个数据点的概率。

当N为有限值时，实际信号经VMD算法分解后模态分量IMFi的样本熵SampEn(m,r,IMFi)可用式(8)估计：

(8)

式中：i为信号分解得到的模态分量序号；r一般取输入时间序列IMFi标准差的0.2倍。若模态分量的自相似性越高，信号复杂度越低，包含的周期性故障信息越多，则样本熵值越小；反之，模态分量的自相似性较低，信号复杂度越高，包含周期性特征越少，则样本熵值越大。

2 故障特征提取流程

由于故障信息往往包含在由故障冲击成分引起的幅值调制信号中，为使所分解的模态分量包含较多的故障特征信息，本文所选目标函数在样本熵的基础上引入了对早期冲击故障较为敏感的峭度指标，各模态分量的峭度qdi定义为：

(9)

式中:E(·)表示括号中变量值的数学期望；μi，σi分别为IMFi的均值与方差；IMFi(j)为第i个分量中的第j个元素。故综合性能指标C作为遗传算法的目标函数可表示为：

C=SampEn(m,r,IMFi)+1/abs(qdi-3)

(10)

具体优化步骤如下：

步骤1，种群初始化。种群中的个体由模态分解个数K和惩罚因子α组成，不同个体的[K,α]组合不同。

步骤2，对初始种群进行选择、交叉和变异。按照目标函数式(10)对种群中每个个体的目标函数值进行求解，得到目标函数最小值Cmin所对应的个体，并记录每一代个体进化的最佳适应度与平均适应度。

步骤3，根据[K,α]组合的最优值，由VMD算法对振动信号进行分解，选取各模态分量与原始信号相关性较大且综合性能指标较小的分量进行Hilbert包络谱分析，得出诊断结果。

基于以上理论，本文所提出的基于改进参数优化VMD的故障特征提取流程如图1所示。

图1 改进参数优化VMD故障特征提取流程

3 仿真信号验证

为验证本文所提方法的有效性，对齿轮箱行星轮故障仿真信号进行特征提取，不失一般性，本文不考虑谐波因素，只考虑载波频率和调制频率的基频与振动传递路径对振动信号的影响[9]，故齿轮箱行星轮故障可用如式(11)的调幅调频信号x(t)表示[10]。

x(t)=[1-cos(2πfrt)][1+Acos(2πfst)]·cos[2πfmt+Bsin(2πfst+φ)]+n(t)

(11)

式中：fr=30 Hz，为行星架的旋转频率；fs=75 Hz，为行星轮的故障特征频率；fm=1 000 Hz，为齿轮啮合频率；A,B为信号的调制系数，均为1；初始相位φ=0；n(t)为高斯白噪声。

为使仿真信号更接近实际现场工况，仿真信号中加入信噪比为-12 dB的高斯白噪声n(t)，采样频率为12 kHz，分析数据点数为8 196个。齿轮箱行星轮故障仿真信号分析结果如图2所示。

图2 仿真信号分析

图2(a)为Fourier频谱图，由图可以看出，啮合频率fm、啮合频率与行星轮故障频率的组合(fm+fs)以及在其两侧出现的以行星架转频的边频带占据主导地位，符合行星轮发生故障时具有的特征；图2(b)为包络谱图，图中并未出现任何特征频率成分，特征信息已被噪声淹没；图2(c)为对VMD算法的参数寻优的过程，设置寻优范围K∈[3,10]，α∈[500,2 000]，得到模态分量与惩罚因子的最优组合为[8,1 397]。对VMD所得到的各模态分量，分别计算加入噪声信号的相关性系数与综合性能指标C，结果见表1。从表中可以看出，IMF2分量与原始信号的相关性系数最大，并且综合性能指标最小，故将IMF2分量作为重构信号进行分析，结果如图3所示。

表1 各模态分量指标计算

图3 重构信号分析

图3(a)与图2(a)对比可以看出，除了有图2(a)中的主导频率以外，出现了额外的特征组合频率(fm+2fs+fr)和(fm+2fs-fr)，并对信号中的噪声有一定抑制作用，提高了信号信噪比。图3(b)与图2(b)对比可以看出，由于噪声影响，图2(b)中并未出现任何故障特征信息，而图3(b)中可以清晰地看出行星架的旋转频率fr、行星轮的故障特征频率fs以及(fs+fr)和(fs-fr)的组合频率，验证了本文所提方法对故障特征提取的有效性。

4 实测数据验证

4.1 试验台验证

如图4所示，采用实验室风电机组传动系统故障试验台验证本文所提方法的有效性。在齿轮箱箱体及所在轴安装美国PCB公司生产的加速度传感器采集振动信号，在高、低速轴分别配备光电编码器采集转速信号，并由变频器控制输出转速，利用自主研制的数据采集仪采集高速轴振动信号，其中试验台参数见表2。

本文在转速为1 200 rad/min、采样频率为12.8 kHz下采集1 s高速轴小齿轮断齿的振动加速度，由试验台参数可计算得到高速轴小齿轮转频fr=1 200/60=20 Hz，齿轮啮合频率fm=20×18=360 Hz。采用遗传算法寻优得到模态分量与惩罚因子的最优组合为[8,1 668]，对VMD分解算法得到的模态分量进行选取、重构，最后进行包络谱分析，结果如图5所示。

图4 风电机组传动系统故障试验台

表2 试验台参数

图5 齿轮断齿数据分析

图5(a)是齿轮断齿状态下未经本文方法优化重构的信号包络谱，从图中可以看出，谱图中并未出现断齿齿轮所在轴的转频与啮合频率。经本文方法处理的包络谱如图5(b)所示，从图中可以清晰地看到断齿齿轮的2倍转频、6倍转频占主导地位，除此之外还有很多的倍转频，因此可以据此判断高速轴齿轮发生了故障。

4.2 工程实例验证

为验证本文所提方法在实际工程中的有效性，利用笔者所在团队自主研发并已实际应用于山西某风电场的状态监测与故障诊断系统，在齿轮箱表面安装加速度传感器，采集风电机组传动链振动数据。已知该风电场2#风机在2017年3月初因齿轮箱故障停机，本文选取的风电场监视控制与数据采集(supervisory control and data acquisition，SCADA)系统实测数据为故障前4个月中每月某时刻的数据，见表3，对其进行分析，根据齿轮箱结构参数和高速轴转速，计算得到齿轮箱二级平行级齿轮的啮合频率为fm=Z2fz/60，其中Z2和fz分别为二级平行级高速轴小齿轮齿数和发电机转速。机组额定转速为1 800 rad/min，齿轮箱由行星级、一级平行级、二级平行级组成，其参数见表4。

表3 风电场实测数据

对该风电场4个时间点的齿轮箱高速轴齿轮振动加速度信号进行分析，结果如图6所示。经过分析发现，在图6(a)中振动信号包络谱较为整齐，且幅值较小，并未出现额定转速下的高速轴啮合频率级转频；而图6(b)、图6(c)、图6(d)均出现对应时间点数据下齿轮啮合频率的1、2、3倍频，且占据主导地位，包络谱幅值较2016年11月22日剧增，说明在2016年11月22日至12月1日之间二级平行级齿轮处出现故障，直至停机，与现场的实际情况相符。

表4 齿轮箱参数

图6 风电场实测故障数据分析

5 结束语

针对风电机组传动链常处于强噪声环境中，而现有VMD算法仍依靠人为经验设定，导致过分解、欠分解、特征提取效果不佳等问题，本文引入样本熵与峭度对VMD算法预设参数进行自适应优化，有效提高输入信号的信噪比，实现故障处的特征频率幅值增强。实验结果表明，本文所提方法对强噪声环境下提取故障特征具有较好的效果，在实际风电机组传动链故障诊断领域有较高的实用价值。