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

原文出處：拓端數(shù)據(jù)公眾號

本文展示了如何使用擴展卡爾曼濾波器進行故障檢測。本文使用擴展的卡爾曼濾波器對一個簡單的直流電機的摩擦力進行在線估計。估計的摩擦力的重大變化被檢測出來，并表明存在故障。

電機模型

電機被模擬成具有阻尼系數(shù)c，轉(zhuǎn)動慣量J，由一個扭矩u驅(qū)動。電機的角速度w和加速度˙w是測量輸出。

為了使用擴展卡爾曼濾波器估計阻尼系數(shù)c，為阻尼系數(shù)引入一個輔助狀態(tài)，并將其導(dǎo)數(shù)設(shè)為零。

因此，模型狀態(tài)，x = [w;c]，和測量，y，方程式為：

使用近似值˙x=xn+1-xnTs將連續(xù)時間方程轉(zhuǎn)換為離散時間，其中Ts為離散采樣周期。這就得到了離散時間模型方程。

指定電機參數(shù)

1. J ?= 10; ? ?% 慣性

2. Ts = 0.01; ?% 采樣時間

指定初始狀態(tài)。

1. x = [...

2. 0; ... % 角速度

3. 1]; % 摩擦力

1. 
   

2. % 以摩擦為狀態(tài)的電機的狀態(tài)更新方程

3. % 狀態(tài)更新方程式

4. x1 = [...

5. x0(1)+Ts/J*(u-x0(1)*x0(2)); ...

6. x0(2)];

1. 
   

2. % 以摩擦為狀態(tài)的電機的測量方程

3. 
   

4. % 輸出。

5. % y - 電機測量元素[角速度；角加速度]。

6. 
   

7. y = [...

8. x(1); ...

9. (u-x(1)*x(2))/J];

電機經(jīng)歷狀態(tài)（過程）噪聲干擾，q，和測量噪聲干擾，r。噪聲項是相加的。

過程和測量噪聲的平均值為零，E[q]=E[r]=0，協(xié)方差Q=E[qq']和R=E[rr']。摩擦狀態(tài)有很高的過程噪聲干擾。這反映了我們希望摩擦力在電機正常運行期間會有變化，并且濾波器能跟蹤這種變化。加速和速度狀態(tài)的噪聲很低，但速度和加速度測量的噪聲相對較大。

指定過程噪聲協(xié)方差。

1. [...

2. 1e-6 0; ... ? % 角速度

3. 


4. 0 1e-2]; ? ? ?% 摩擦力

指定測量噪聲協(xié)方差。

1. [...

2. 1e-4 0; ... ?% 速度測量

3. 0 1e-4]; ? ? % 加速度測量

創(chuàng)建一個擴展的卡爾曼過濾器

創(chuàng)建一個擴展的卡爾曼濾波器來估計模型的狀態(tài)。我們特別關(guān)注阻尼狀態(tài)，因為這個狀態(tài)值的急劇變化表明存在故障事件。
創(chuàng)建一個擴展卡爾曼濾波器對象，并指定狀態(tài)轉(zhuǎn)換和測量函數(shù)的雅各布系數(shù)。

擴展卡爾曼濾波器的輸入?yún)?shù)是之前定義的狀態(tài)轉(zhuǎn)換和測量函數(shù)。初始狀態(tài)值x0、初始狀態(tài)協(xié)方差、過程和測量噪聲協(xié)方差也是擴展卡爾曼濾波器的輸入。在這個例子中，精確的雅各布函數(shù)可以從狀態(tài)轉(zhuǎn)換函數(shù)f和測量函數(shù)h中得到。

1. 


2. % 輸出

3. Jac - 在x處計算出的狀態(tài)雅各布系數(shù)

4. 


5. 


6. % 雅各布系數(shù)

7. Jac = [

8. 1-Ts/J*x(2) -Ts/J*x(1); ...

9. 0 1];

10. 


11. % 電機模型測量方程的雅各布系數(shù)

12. 


13. % 輸出

14. Jac - 在 x 處計算的測量雅各布系數(shù)

15. 


16. 


17. % 雅各布系數(shù)

18. J = [ ...

19. 1 0;

20. -x(2)/J -x(1)/J];

Simulation仿真

為了模擬工廠，創(chuàng)建一個環(huán)路，在電機中引入一個故障（虛構(gòu)的電機劇烈變化）。在模擬回路中，使用擴展的卡爾曼濾波器來估計電機狀態(tài)，并特別跟蹤摩擦狀態(tài)，檢測摩擦力何時發(fā)生統(tǒng)計意義上的變化。

電機被模擬成一個脈沖序列，反復(fù)加速和減速。這種類型的電機操作對于生產(chǎn)線上的采摘機器人來說是典型的。

在模擬電機時，加入與構(gòu)建擴展卡爾曼濾波器時使用的Q和R噪聲協(xié)方差值相似的過程和測量噪聲。對于摩擦，使用一個小得多的噪聲值，因為除了故障發(fā)生時，摩擦大多是恒定的。在模擬過程中人為地誘發(fā)故障。

1. Qv = chol(Q); % 過程噪聲的標(biāo)準(zhǔn)偏差

2. Qv(end) = 1e-2; % 較小的摩擦噪聲

3. Rv = chol(R); % 測量噪聲的標(biāo)準(zhǔn)偏差

使用狀態(tài)更新方程對模型進行仿真，并在模型狀態(tài)中加入過程噪聲。仿真十秒鐘后，強制改變電機的摩擦力。使用模型測量功能來模擬電機傳感器，并在模型輸出中加入測量噪聲。

1. for ct = 1:numel(t)

2. % 模型輸出更新

3. y = y+Rv*randn(2,1); % 添加測量噪聲

4. % 模型狀態(tài)更新

5. xSig(:,ct) = x0;

6. % 誘發(fā)摩擦力的變化

7. if t(ct) == 10

8. x1(2) = 10; % 步驟變化

9. 


10. x1n = x1+Qv*randn(nx ?% 加入過程噪聲

Significant friction change at 10.450000

為了從電機測量值中估計電機狀態(tài)，使用擴展卡爾曼濾波器的預(yù)測和糾正命令。

1. 


2. 


3. % 使用擴展卡爾曼濾波器進行狀態(tài)估計

4. x_corr = correct(ekf,y,u(ct),J,Ts); % 根據(jù)當(dāng)前測量結(jié)果修正狀態(tài)估計。

5. predict(ekf,u(ct),J,Ts); % 根據(jù)當(dāng)前狀態(tài)和輸入預(yù)測下一個狀態(tài)。

為了檢測摩擦力的變化，使用4秒的移動窗口來計算估計的摩擦力平均值和標(biāo)準(zhǔn)偏差。在最初的7秒之后，鎖定計算的平均值和標(biāo)準(zhǔn)差。這個最初計算出的平均值是摩擦力的預(yù)期無故障平均值。7秒后，如果估計的摩擦力與預(yù)期的無故障平均值相差超過3個標(biāo)準(zhǔn)差，這就意味著摩擦力有了明顯的變化。為了減少噪音和估計摩擦力的影響，在與3個標(biāo)準(zhǔn)差的界限比較時，使用估計摩擦力的平均值。

1. % 計算估計平均值和標(biāo)準(zhǔn)偏差。

2. else

3. % 存儲計算出的平均數(shù)和標(biāo)準(zhǔn)差，不需要

4. %重新計算。

5. fMean(ct) = fMean(ct-1)

6. % 使用預(yù)期的摩擦力平均值和標(biāo)準(zhǔn)偏差來檢測

7. %摩擦力變化。

8. estFriction = mean(xSigEst(2,

9. if fChanged(ct) && ~fChanged(ct-1)

10. % 檢測摩擦變化信號的上升沿|fChanged|

使用估計的狀態(tài)來計算估計的輸出。計算測量輸出和估計輸出之間的誤差，并計算出誤差統(tǒng)計。誤差統(tǒng)計可用于檢測摩擦力的變化。這一點將在后面詳細討論。

1. 


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

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

擴展的卡爾曼濾波器性能

請注意，在10.45秒時檢測到了一個摩擦變化。我們現(xiàn)在描述一下這個故障檢測規(guī)則是如何得出的。首先檢查仿真結(jié)果和過濾器的性能。

1. figure,

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

模型的輸入輸出響應(yīng)表明，很難直接從測量信號中檢測出摩擦力的變化。擴展的卡爾曼濾能夠估計狀態(tài)，特別是摩擦狀態(tài)。比較真實的模型狀態(tài)和估計狀態(tài)。估計的狀態(tài)顯示了對應(yīng)于3個標(biāo)準(zhǔn)差的置信區(qū)間。

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

2. 


請注意，濾波器的估計值跟蹤了真實值，而且置信區(qū)間仍然有界。檢查估計誤差可以更深入地了解濾波器。

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

誤差圖顯示，濾波器在10秒的摩擦力變化后進行了調(diào)整，并將估計誤差降低到了零。然而，誤差圖不能用于故障檢測，因為它們依賴于對真實狀態(tài)的了解。將測量的狀態(tài)值與加速度和速度的估計狀態(tài)值進行比較，可以提供一種檢測機制。

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

加速度誤差圖顯示，在引入故障的10秒左右，平均誤差有微小的差異。查看誤差統(tǒng)計，看看是否可以從計算的誤差中檢測出故障。加速度和速度誤差預(yù)計是正態(tài)分布的（噪聲模型都是高斯的）。因此，加速度誤差的峰度可能有助于識別由于摩擦力的變化和由此產(chǎn)生的誤差分布從對稱變?yōu)椴粚ΨQ的變化情況。

plot(t,Kur(1,:)

忽略估計器仍在收斂和收集數(shù)據(jù)的前4秒，誤差的峰度相對穩(wěn)定，在3（高斯分布的預(yù)期峰度值）附近有微小的變化。因此，在這個應(yīng)用中，誤差統(tǒng)計不能被用來自動檢測摩擦力的變化。在這個應(yīng)用中，使用誤差的峰度也是很困難的，因為過濾器正在適應(yīng)并不斷地將誤差推向零，只給出了一個誤差分布與零不同的短暫時間窗口。

因此在這個應(yīng)用中，使用估計的摩擦力的變化提供了自動檢測電機故障的最好方法。來自已知無故障數(shù)據(jù)的摩擦力估計值（平均值和標(biāo)準(zhǔn)偏差）提供了摩擦力的預(yù)期界限，當(dāng)這些界限被超過時，很容易檢測出來。下面的圖強調(diào)了這種故障檢測方法。

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

摘要

這個例子展示了如何使用擴展的卡爾曼濾波器來估計一個簡單的直流電動機的摩擦力，并使用摩擦力估計值進行故障檢測。

最受歡迎的見解

1.matlab使用經(jīng)驗?zāi)Ｊ椒纸鈋md 對信號進行去噪

2.Matlab使用Hampel濾波去除異常值

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

4.matlab預(yù)測ARMA-GARCH 條件均值和方差模型

5.matlab中使用VMD(變分模態(tài)分解)?

6.matlab使用貝葉斯優(yōu)化的深度學(xué)習(xí)

7.matlab貝葉斯隱馬爾可夫hmm模型

8.matlab中的隱馬爾可夫模型(HMM)實現(xiàn)

9.matlab實現(xiàn)MCMC的馬爾可夫切換ARMA – GARCH模型