基于特征提取和多层感知器的轴承故障检测

摘要:本文是对DC竞赛网站上的轴承故障检测题目的解答。利用此题目中的轴承振动数据,我们使用MATLAB编程从时域和频域提取特征,完成故障分类。利用此种方法,我们实现了100%的检测正确率。这种思路是传统的模式识别方法,相较于卷积神经网络具有较快的训练速度和在线计算速度。

特征提取

与如今广泛应用的深度学习方法不同,使用简单反向传播(BP)算法的多层感知器不能处理过高的输入数据维度。因此,在讨论神经网络分类器的训练之前,我们首先谈谈如何进行特征提取。容易想到,特征提取的两种比较直接的思路是在频域和在时域内进行提取,下面将分别进行讨论。此外,还有小波分析和时频域功率谱等分析方法,本文不予讨论。

时域特征提取

如题目所言,轴承的故障分为内圈、外圈和滚珠故障。在作出轴承的振动图形后,我们发现故障轴承的振动特征是:每间隔一段时间,振动幅值周期性地增大。正常轴承则不存在此现象。待检测轴承还分为三种不同直径。同样作出图形进行比较,可以注意到直径越大,异常振动的幅值也就越大。因此,通过时域特征进行故障识别和分类是可行的。右图绘制出了不同工作状态的轴承的振动信号。

在题目要求中已经指出,数据中各个采样点不是对齐的,不能将其作为特征进行输入。(经笔者测试,使用卷积神经网络直接输入时域序列,可以得到99%的正确率。这是因为卷积层本身具有特征提取的功能。)因此,我们选择时域数据的统计特征,包括:

  • 平均值 \mu_a=\frac{1}{n}\sum_{k=1}^na_k
  • 方均根值,又称为有效值,表征振动能量的强度. a_{RMS}=\sqrt{\frac{1}{n}\sum_{k=1}^{n}a_k^2}
  • 最大值/方均根值。可表征前文指出的周期性振动增强。 \frac{max\{a_k\}}{a_{RMS}}
  • 2-10阶中心矩。部分阶次的中心矩具有明确的物理意义。阶次k=2,即方差。k=3,称为偏度,即随机变量相对于平均值分布的不对称程度。k=4,称为峰度,即变量分布集中的程度(就是概率密度函数的“峰”)。这些变量中有一些应该是无用的,笔者并未逐个分析或尝试,读者可以自行移除其中一部分后在训练网络,观察现象。 \gamma_i=\frac{1}{n}\sum_{k=1}{n}(a_k-\mu_a)^i, i=2,3,\dots ,9,10

下面我们使用MATLAB编程求取上述特征。假设我们已将振动数据导入到792*6000大小的data变量中,并将对应标签存储值792*1的label变量中。需要指出,以下所有数据都必须在归一化之后才能输入神经网络

num_samples = 792;
ave = zeros([1 num_samples]);
for i=1:num_samples %求平均值
    sum = 0;
    for j=1:6000
        sum = sum + data(i,j);
    end
    ave(i) = sum / 6000;
end
central_moments = zeros([10,num_samples]);
for k=2:10  %求k阶中心矩
    for i=1:num_samples
        sum = 0;
        for j=1:6000
            sum = sum + (data(i,j)-ave(i))^k;
        end
        central_moments(k,i) = sum / 6000;
    end
end
mean_square_root = zeros([1 num_samples]);
for i=1:num_samples %下面计算方均根值
    sum = 0;
    for j=1:6000
        sum = sum + data(i,j)^2;
    end
    mean_square_root(i) = sqrt(sum/6000);
end
max_ = zeros([1 num_samples]); %下面计算最大绝对值与有效值之比
for i=1:num_samples
    max_ = max(abs(data(i,:)));
end
peak_d_msr = max_ ./ mean_square_root;

频域特征提取

从频域入手分析信号是常见的思路。实际上,对振动信号的分析与语音数据的分析是类似的。我们采用快速傅里叶变换(Fast Fourier Transform, FFT)求取振动信号的频谱。MATLAB的FFT函数求得的是复数域的双侧频谱。我们对其取模求得实数域上的单侧频谱(即幅频特性,舍弃掉相位信息)。不同故障类型的振动频谱如下图所示。我们并不了解振动数据的采样频率,因此只能使用“归一化频率”的概念。时域共有6000个采样点,则频谱共有3001个点,这是由奈奎斯特采样定理决定的。频谱第1个点为0频率分量,第3001个点对应频率称为1倍归一化频率(即时域采样频率的1/2倍)。左侧图形是完整的频谱图象,右侧图形则是0-0.1倍归一化频率的低频段图象。

频谱图象
低频段对应频谱图象

我们采用以下方式提取频域特征:在频谱的前100个点中,每隔10个点求一个峰值;之后每隔200个点求一个峰值,作为频域特征。事实上,这种做法比较粗糙,但是已经能满足此处问题的需要。更加合理的做法是利用一个“滑动窗口”求峰值,比如1-100,50-150,100-200分别取峰值。这种做法与卷积神经网络中利用卷积核做自动特征提取的思路是类似的。

至此,我们在时域和频域中提取了合计33个特征。

fourier_complex = zeros([num_samples 6000]);
for i = 1:num_samples %求双侧频谱(结果为复数)
    fourier_complex(i,:) = fft(data(i,1:6000));
end
nor = zeros([num_samples 6000]);
fourier = zeros([num_samples,3001]);
for i = 1:num_samples %求单侧频谱(幅频特性)
    nor(i,:) = abs (fourier_complex(i,:) / 6000);
    fourier(i,:) = nor(i,1:3001);
    fourier(i,2:end-1) = 2*fourier(i,2:end-1); 
    %除了0频率分量都要乘以2,因为前面只取了双侧频谱的一半
end

fourier_peaks = zeros(12,num_samples);
for i = 1:10 %前100个点每隔10个点求一个最大值
    for j = 1:num_samples
    fourier_peaks(i,j) = max(fourier(j,(i*10-9):(i*10)));
end
    
end
for i = 2:12
    for j = 1:num_samples
        start_ = 100+i*200-199;
        end_ = 200*i+100;
        fourier_peaks(i+9,j) = max(fourier(j,start_:end_));
    end
end

归一化处理

我们使用MATLAB中的normalize函数将所有样本的各个特征归一化,归一化后的均值为0,方差为1。

%合并为一个特征
final_res = [central_moments;normalized_mean_square_root;normalized_peak_d_msr;fourier_peaks];
for i=1:33
    normalize(final_res(i,:))
end

多层感知器的建构与训练

MATLAB中有一款Neural Net Pattern Recognition工具,能提供一个简单的两层感知器网络。此网络的结构如下图所示。

此网络的输入为前文中提取的33个特征,隐含层神经元激活函数为Sigmoid函数,输出层神经元函数为Softmax函数,输出采用one-hot形式编码。我们将在下文对它们进行简要解释。

待训练参数:权重(Weight)和偏置(Bias)

在网络结构图中看到的W和b即为这两组参数的缩写。权重不在赘述,偏置就是与输入无关,直接加在刺激强度(即激活函数的自变量)上的一个值。这样做的目的是使得神经元可以对远离原点的输入值作出响应。在有些教材和教程的神经网络中是没有Bias这一个参数的,这种神经网络的输入值必须被严格归一化至[-1,1]的区间上,因为对于Sigmoid等激活函数而言,200和2000的刺激强度没有什么区别,其输出都非常接近于1。

Sigmoid激活函数

Sigmoid函数的公式是。选取这一激活函数的目标就是因其可以向网络引入非线性(形如y=kx的线性函数则不行),并且可以使得输出对隐含层输入权重有连续偏导数,可以应用BP(反向传播)算法。这一点在各种教程教材中均反复讲述,因此我们不再赘述。

Softmax激活函数

首先,从这一函数的名称可以直观地理解其功能:输出反应输入的最大值,但又不是“硬”的,一者为1其余为0的形式。其次,此激活函数对于输出层的所有神经元是一个整体,即每个神经元的输出与此层各个神经元上输入的刺激强度均有关。Softmax函数的表达式为:

\sigma(\textbf{z})_j=\frac{e^{Z_j}}{\sum_{k=1}^Ke^{z_k}}

其中自变量z是一个向量,即各神经元的刺激强度;j是输出量的下标。

性能指标:交叉熵(Cross Entropy)

L=-\frac{1}{N} \sum\limits^{N}_{i=1}(yln(\hat{y})+(1-y)ln(1-\hat{y}))

其中:N为待分类的类数,y是预期的输出(1或0),\hat{y}是神经网络实际产生的输出。采用交叉熵的主要原因是可以克服方差代价函数更新权重过慢的问题。

训练方法:变步长共轭梯度法(Scaled Conjugate Gradient Method)

这里不打算对此种算法作数学推导,而仅仅指出变步长共轭梯度法的几点性质:

  • 与梯度下降法相比,此方法的搜索方向不止与梯度有关,还与上一步搜索的方向有关。准确地讲,二者间为共轭关系。
  • 使用Levenberg-Marquardt方法改变搜索步长。
  • 此种方法可以显著提高网络的收敛速度。

训练结果

我们选用50个的隐含层神经元数量。数据集按MATLAB默认设置,取70%为训练集,15%为测试集,15%为验证集。训练结果如下图所示。此多层感知器网络对所有验证样本分类正确,并且经验证,可在题目给出的测试集中取得100%的正确率。

发表评论

电子邮件地址不会被公开。 必填项已用*标注