拓端数据tecdat|matlab实现扩展卡尔曼滤波(EKF)进行故障检测

原文链接：http://tecdat.cn/?p=22467 

原文出处：拓端数据公众号

本文展示了如何使用扩展卡尔曼滤波器进行故障检测。本文使用扩展的卡尔曼滤波器对一个简单的直流电机的摩擦力进行在线估计。估计的摩擦力的重大变化被检测出来，并表明存在故障。

电机模型

电机被模拟成具有阻尼系数c，转动惯量J，由一个扭矩u驱动。电机的角速度w和加速度˙w是测量输出。

为了使用扩展卡尔曼滤波器估计阻尼系数c，为阻尼系数引入一个辅助状态，并将其导数设为零。

因此，模型状态，x = [w;c]，和测量，y，方程式为：

使用近似值˙x=xn+1-xnTs将连续时间方程转换为离散时间，其中Ts为离散采样周期。这就得到了离散时间模型方程。

指定电机参数

1.  
   J = 10; % 惯性
2.  
   Ts = 0.01; % 采样时间

指定初始状态。

1.  
   x = [...
2.  
   0; ... % 角速度
3.  
   1]; % 摩擦力

1.  
    
2.  
   % 以摩擦为状态的电机的状态更新方程
3.  
   % 状态更新方程式
4.  
   x1 = [...
5.  
   x0(1)+Ts/J*(u-x0(1)*x0(2)); ...
6.  
   x0(2)];

1.  
    
2.  
   % 以摩擦为状态的电机的测量方程
3.  
    
4.  
   % 输出。
5.  
   % y - 电机测量元素[角速度；角加速度]。
6.  
    
7.  
   y = [...
8.  
   x(1); ...
9.  
   (u-x(1)*x(2))/J];

电机经历状态（过程）噪声干扰，q，和测量噪声干扰，r。噪声项是相加的。

过程和测量噪声的平均值为零，E[q]=E[r]=0，协方差Q=E[qq']和R=E[rr']。摩擦状态有很高的过程噪声干扰。这反映了我们希望摩擦力在电机正常运行期间会有变化，并且滤波器能跟踪这种变化。加速和速度状态的噪声很低，但速度和加速度测量的噪声相对较大。

指定过程噪声协方差。

1.  
   [...
2.  
   1e-6 0; ... % 角速度
3.  
    
4.  
   0 1e-2]; % 摩擦力

指定测量噪声协方差。

1.  
   [...
2.  
   1e-4 0; ... % 速度测量
3.  
   0 1e-4]; % 加速度测量

创建一个扩展的卡尔曼过滤器

创建一个扩展的卡尔曼滤波器来估计模型的状态。我们特别关注阻尼状态，因为这个状态值的急剧变化表明存在故障事件。
创建一个扩展卡尔曼滤波器对象，并指定状态转换和测量函数的雅各布系数。

扩展卡尔曼滤波器的输入参数是之前定义的状态转换和测量函数。初始状态值x0、初始状态协方差、过程和测量噪声协方差也是扩展卡尔曼滤波器的输入。在这个例子中，精确的雅各布函数可以从状态转换函数f和测量函数h中得到。

1.  
    
2.  
   % 输出
3.  
   Jac - 在x处计算出的状态雅各布系数
4.  
    
5.  
    
6.  
   % 雅各布系数
7.  
   Jac = [
8.  
   1-Ts/J*x(2) -Ts/J*x(1); ...
9.  
   0 1];
10.  
     
11.  
    % 电机模型测量方程的雅各布系数
12.  
     
13.  
    % 输出
14.  
    Jac - 在 x 处计算的测量雅各布系数
15.  
     
16.  
     
17.  
    % 雅各布系数
18.  
    J = [ ...
19.  
    1 0;
20.  
    -x(2)/J -x(1)/J];

Simulation仿真

为了模拟工厂，创建一个环路，在电机中引入一个故障（虚构的电机剧烈变化）。在模拟回路中，使用扩展的卡尔曼滤波器来估计电机状态，并特别跟踪摩擦状态，检测摩擦力何时发生统计意义上的变化。

电机被模拟成一个脉冲序列，反复加速和减速。这种类型的电机操作对于生产线上的采摘机器人来说是典型的。

在模拟电机时，加入与构建扩展卡尔曼滤波器时使用的Q和R噪声协方差值相似的过程和测量噪声。对于摩擦，使用一个小得多的噪声值，因为除了故障发生时，摩擦大多是恒定的。在模拟过程中人为地诱发故障。

1.  
   Qv = chol(Q); % 过程噪声的标准偏差
2.  
   Qv(end) = 1e-2; % 较小的摩擦噪声
3.  
   Rv = chol(R); % 测量噪声的标准偏差

使用状态更新方程对模型进行仿真，并在模型状态中加入过程噪声。仿真十秒钟后，强制改变电机的摩擦力。使用模型测量功能来模拟电机传感器，并在模型输出中加入测量噪声。

1.  
   for ct = 1:numel(t)
2.  
   % 模型输出更新
3.  
   y = y+Rv*randn(2,1); % 添加测量噪声
4.  
   % 模型状态更新
5.  
   xSig(:,ct) = x0;
6.  
   % 诱发摩擦力的变化
7.  
   if t(ct) == 10
8.  
   x1(2) = 10; % 步骤变化
9.  
    
10.  
    x1n = x1+Qv*randn(nx % 加入过程噪声

Significant friction change at 10.450000

为了从电机测量值中估计电机状态，使用扩展卡尔曼滤波器的预测和纠正命令。

1.  
    
2.  
    
3.  
   % 使用扩展卡尔曼滤波器进行状态估计
4.  
   x_corr = correct(ekf,y,u(ct),J,Ts); % 根据当前测量结果修正状态估计。
5.  
   predict(ekf,u(ct),J,Ts); % 根据当前状态和输入预测下一个状态。

为了检测摩擦力的变化，使用4秒的移动窗口来计算估计的摩擦力平均值和标准偏差。在最初的7秒之后，锁定计算的平均值和标准差。这个最初计算出的平均值是摩擦力的预期无故障平均值。7秒后，如果估计的摩擦力与预期的无故障平均值相差超过3个标准差，这就意味着摩擦力有了明显的变化。为了减少噪音和估计摩擦力的影响，在与3个标准差的界限比较时，使用估计摩擦力的平均值。

1.  
   % 计算估计平均值和标准偏差。
2.  
   else
3.  
   % 存储计算出的平均数和标准差，不需要
4.  
   %重新计算。
5.  
   fMean(ct) = fMean(ct-1)
6.  
   % 使用预期的摩擦力平均值和标准偏差来检测
7.  
   %摩擦力变化。
8.  
   estFriction = mean(xSigEst(2,
9.  
   if fChanged(ct) && ~fChanged(ct-1)
10.  
    % 检测摩擦变化信号的上升沿|fChanged|

使用估计的状态来计算估计的输出。计算测量输出和估计输出之间的误差，并计算出误差统计。误差统计可用于检测摩擦力的变化。这一点将在后面详细讨论。

1.  
    
2.  
   kurtosis(ySigEst(1,idx)-ySig(1,idx));
3.  
   kurtosis(ySigEst(2,idx)-ySig(2,idx))];

扩展的卡尔曼滤波器性能

请注意，在10.45秒时检测到了一个摩擦变化。我们现在描述一下这个故障检测规则是如何得出的。首先检查仿真结果和过滤器的性能。

1.  
   figure,
2.  
   plot(t, Sig(1,:) Sig(2,:));

模型的输入输出响应表明，很难直接从测量信号中检测出摩擦力的变化。扩展的卡尔曼滤能够估计状态，特别是摩擦状态。比较真实的模型状态和估计状态。估计的状态显示了对应于3个标准差的置信区间。

1.  
   plot(t, True(1,:), t, Est(1,:), ...
2.  
    

请注意，滤波器的估计值跟踪了真实值，而且置信区间仍然有界。检查估计误差可以更深入地了解滤波器。

plot(t,True(1,:)-Est(1,:)
 

误差图显示，滤波器在10秒的摩擦力变化后进行了调整，并将估计误差降低到了零。然而，误差图不能用于故障检测，因为它们依赖于对真实状态的了解。将测量的状态值与加速度和速度的估计状态值进行比较，可以提供一种检测机制。

plot(t,Sig(1,:-Est(1,:)

加速度误差图显示，在引入故障的10秒左右，平均误差有微小的差异。查看误差统计，看看是否可以从计算的误差中检测出故障。加速度和速度误差预计是正态分布的（噪声模型都是高斯的）。因此，加速度误差的峰度可能有助于识别由于摩擦力的变化和由此产生的误差分布从对称变为不对称的变化情况。

plot(t,Kur(1,:)

忽略估计器仍在收敛和收集数据的前4秒，误差的峰度相对稳定，在3（高斯分布的预期峰度值）附近有微小的变化。因此，在这个应用中，误差统计不能被用来自动检测摩擦力的变化。在这个应用中，使用误差的峰度也是很困难的，因为过滤器正在适应并不断地将误差推向零，只给出了一个误差分布与零不同的短暂时间窗口。

因此在这个应用中，使用估计的摩擦力的变化提供了自动检测电机故障的最好方法。来自已知无故障数据的摩擦力估计值（平均值和标准偏差）提供了摩擦力的预期界限，当这些界限被超过时，很容易检测出来。下面的图强调了这种故障检测方法。

plot(t,x,[nan t],[Mean+3*STD,Mean-3*STD]

摘要

这个例子展示了如何使用扩展的卡尔曼滤波器来估计一个简单的直流电动机的摩擦力，并使用摩擦力估计值进行故障检测。

---

最受欢迎的见解

1.matlab使用经验模式分解emd 对信号进行去噪

2.Matlab使用Hampel滤波去除异常值

3.matlab偏最小二乘回归(PLSR)和主成分回归(PCR)

4.matlab预测ARMA-GARCH 条件均值和方差模型

5.matlab中使用VMD(变分模态分解) 

6.matlab使用贝叶斯优化的深度学习

7.matlab贝叶斯隐马尔可夫hmm模型

8.matlab中的隐马尔可夫模型(HMM)实现

9.matlab实现MCMC的马尔可夫切换ARMA – GARCH模型