4.1 仿真事例介绍

当风电机组传动系统出现故障时，故障振动信号往往表现为调制形式，因此解调分析成为风电机组传动系统故障诊断的一种重要的信号处理方法。本文采用近年来被广泛使用的基于EMD的Hilbert解调法，即通过EMD分解自适应地将多分量的调制故障信号分解为若干个IMF之和，每一个IMF分量都是单分量的调幅和调频信号（AM-FM信号）：

对表征工频背景信号的IMF分量之前的表征故障调制特征的高频IMF分量进行Hilbert变换解调就得到故障信号的瞬时幅值：

式中 c _i (t)为对各IMF分量进行Hilbert变换后的信号， 以上各IMF分量的瞬时幅值包含丰富的故障信息，为提取其中的故障特征测度，本文提出对以上IMF分量瞬时幅值构建Shannon熵：

各类故障样本的归一化瞬时幅值Shannon熵是信息熵的一种，是衡量整个信号源在总体上的平均不确定性的量度，Shannon熵越大，信息越确定，信息的无序程度越小，反之越大。正常设备信号的能量是随机分布的，它携带信息的不确定度大，Shannon熵值就小；异常设备的信号中增加了确定性故障成分，且这种成分随故障程度的增加而增加，Shannon熵值便增大，因此在设备故障诊断领域可利用Shannon熵这种对信号整体能量分布特性的概括能力，即用Shannon熵值检测信号概率分布的均匀性来反映设备所处的状态。把训练样本的归一化Shannon熵矢量T 输入优化LS-LPWSVM对其进行训练，以训练好的优化LS-LPWSVM来判别待测样本的故障位置、类型或程度。

4.2仿真过程

取三个仿真信号 x(t)、y(t)和 z(t)，其信号表达式分别为：

x(t)=5cos(20πt)+10cos(40πt)+15cos(60πt)+20cos(80πt)

y(t)=5[1+cos(4πt)]cos(20πt)+10[1+cos(4πt)]cos(40πt)+15[1+cos(4πt)]cos(60πt)+20[1+cos(4πt)]cos(80πt)

z(t)=10+5cos(10πt)cos2πt+2cos(5πt)+sin(25πt)sin(30πt)+sin(40πt)

其中x(t)是一些谐波信号的合成， y (t)为调幅信号的合成，z (t)包含 3个分量，1个单分量 AM-FM 信号，1个调幅信号和 1个调频信号。对这些信号进行离散，采样频率为 1024 Hz，长度为60 ×1s。它们的时域波形如下图所示。

将每个仿真信号分成60组样本，显然每个样本时间长度为1s，含有1024个数据点。从中随机抽取10组用于样本训练，其余50组作为测试样本，由于需要区分三类信号，因此需要2层LS-SVM,第1层LS-SVM将信号x(t)从三类信号中分离出来，第2层LS-SVM将y(t)和z(t)区分开来。将作为训练样本的三类信号的Shannon熵矢量及其模式编号（x(t)→1、y(t)→2、z(t)→3）输入到两层均设置为Poly核，γ=10，d=1，c=1的LS-SVM中进行训练，其中第1层LS-SVM用来训练三种故障样本，第2层LS-SVM用来训练滚动体故障和外圈故障样本。再将作为测试样本的Shannon熵矢量输入到训练好的两层LS-SVM中去，

4.3 Matlab软件源程序代码

clear;

for     a=1:1:60;

       t=a-1:1/1024:a; %      data(a,:)=5*cos(20*pi*t)+10*cos(40*pi*t)+15*cos(60*pi*t)+20*cos(80*pi*t);

data(60+a,:)=5*(1+cos(4*pi*t)).*cos(20*pi*t)+10*(1+cos(4*pi*t)).*cos(40*pi*t)+15*(1+cos(4*pi*t)).*cos(60*pi*t)+20*(1+cos(4*pi*t)).*cos(80*pi*t);        data(120+a,:)=(10+5*cos(10*pi*t)).*cos(2*pi*t+2*cos(5*pi*t))+4*sin(25*pi*t).*sin(30*pi*t)+sin(40*pi*t.*t);

end

for     i=1:1:10;

        r=fix(120*rand(1))

        traindata(i,:)=data(r,:)%选取十组作为训练样本     

end

for     i=1:1:10;

        imf=emd(traindata(i,:));

%经验模态分解:训练数据里有10组，每个得到2个imf信号

        [A,f,tt] = hhspectrum(imf(1:end-1,:));

%希尔伯特变换得到幅值A和频率f

        AA=A(1:2,:)

%取imf的前两层

%以下构建熵的特征向量T

%以下构建x的每个训练样本的熵特征向量sxT(:,i)

        E=(AA.^2).*log(AA.^2)  

        for z=1:1:2

          for j=1:1:1023

             EX(z,i)=sum(E(z,j)); 

          end

        end

%利用幅值A求第i层imf分量的熵值E（i）

end

%以下对shang进行归一化处理

      STx =mapminmax(EX,0,1); 

%得到归一化瞬时幅值熵值矢量

gam = 10; %正则参数

sig2 =1;%RBF核函数的参数

type = 'classification';%用于分类

sx=STx

sx=sx'

sy=sign(sx(:,1)-sx(:,2));

[alpha,b] = trainlssvm({sx,sy,type,gam,sig2,'RBF_kernel'});

[alpha,b] = trainlssvm({sx,sy,type,gam,sig2,'RBF_kernel','original'});

[alpha,b] = trainlssvm({sx,sy,type,gam,sig2,'RBF_kernel','preprocess'});

%以上构造检验样本。Ytest是计算出来的检验值%

plotlssvm({sx,sy,type,gam,sig2,'RBF_kernel'},{alpha,b});%画图%