【MATLAB】卡尔曼滤波器的原理及仿真(初学者专用)
文章目錄
- 0.引言
- 1.場景預設
- 2.卡爾曼濾波器
- 3.仿真及效果
0.引言
\qquad 本文參考了Matlab對卡爾曼濾波器的官方教程及幫助文檔(Kalman Filter)。官方教程的B站鏈接如下,在此對分享資源的Up主表示感謝。(如不能正常播放或需要看中文字幕,請點擊此處B站鏈接)
【官方教程】卡爾曼濾波器教程與MATLAB仿真(全)(中英字幕)
另外提供友情鏈接如下:
\qquad本文不會完全照搬視頻中的所有內容,只會介紹有關卡爾曼濾波器關于定位方面的知識。卡爾曼濾波器除最原始的版本(KF)外,其延伸版本主要有三種——擴展卡爾曼濾波器(EKF)、無跡卡爾曼濾波器(UKF)、粒子濾波器(PF)。它們的運算復雜度依次遞增,其中KF、EKF、UKF是建立在隨機噪聲是高斯分布的基礎上的;PF理論則沒有此預設。它們的關系如下圖所示:
\qquad首先,KF是一種狀態觀測器。狀態觀測器是針對可觀系統,根據輸出yyy和輸入uuu對系統內部狀態uuu進行觀測的結構單元。其構圖如下(參考視頻鏈接):
原系統方程為:
{x˙=Ax+Buy=Cx\begin{cases} \dot{x}=Ax+Bu \\ y=Cx \end{cases} {x˙=Ax+Buy=Cx?
采用狀態觀測器的觀測系統方程為:
{x^˙=Ax^+Bu+K(y?y^)y^=Cx^\begin{cases} \dot{\hat{x}}=A\hat{x}+Bu+K(y-\hat{y}) \\ \hat{y}=C\hat{x} \end{cases} {x^˙=Ax^+Bu+K(y?y^?)y^?=Cx^?
為保證觀測器的lim?t→∞(x?x^)T(x?x^)=0\lim_{t\rightarrow\infty}(x-\hat{x})^T(x-\hat{x})=0t→∞lim?(x?x^)T(x?x^)=0
需要A?KCA-KCA?KC是負定的。而卡爾曼濾波器就是一種狀態觀測器,只不過它是隨機系統的狀態觀測器,其結構框圖如下:
在輸入uuu和動態系統Plant中間會引入過程噪聲www,而在輸出yyy(即測量)和實際測量yvy_vyv?中間會引入傳感器噪聲vvv,而卡爾曼濾波器則是根據u,yvu,y_vu,yv?求得測量的最優估計yey_eye?。
1.場景預設
\qquad 假設一輛汽車在做勻速直線運動(一維場景),駕駛員通過油門控制了汽車的加速度恒定。忽略汽車所受的阻力、質量變化,并假設駕駛員操作會給汽車的牽引力造成一定的過程噪聲。選擇汽車的位移和速度為狀態變量x=[p,v]Tx=[p,v]^Tx=[p,v]T,則狀態變量導數為汽車的速度和加速度即x˙=[v,a]T\dot{x}=[v,a]^Tx˙=[v,a]T,選取控制變量u=au=au=a,測量量為汽車的位移。則狀態方程如下:
{x˙=[0100]x+[01]uy=[01]x\begin{cases} \dot{x}=\begin{bmatrix} 0 &1\\0 & 0 \end{bmatrix}x+\begin{bmatrix} 0 \\ 1 \end{bmatrix}u\\[2ex] y=\begin{bmatrix} 0 & 1 \end{bmatrix}x \end{cases} ??????x˙=[00?10?]x+[01?]uy=[0?1?]x?
選取采用間隔dtdtdt,則狀態方程離散化后變為:
{x(n)=[1dt01]x(n?1)+[0dt]u(n?1)y(n)=[01]x(n)\begin{cases} x(n)=\begin{bmatrix} 1 &dt \\ 0 & 1 \end{bmatrix}x(n-1)+\begin{bmatrix} 0 \\ dt \end{bmatrix}u(n-1)\\[2ex] y(n)=\begin{bmatrix} 0 & 1 \end{bmatrix}x(n) \end{cases} ??????x(n)=[10?dt1?]x(n?1)+[0dt?]u(n?1)y(n)=[0?1?]x(n)?
按照慣例定義
A=[1dt01],B=[0dt],C=[01],I=[1001]A=\begin{bmatrix} 1 &dt \\ 0 & 1 \end{bmatrix},B=\begin{bmatrix} 0 \\ dt \end{bmatrix},C=\begin{bmatrix} 0 & 1 \end{bmatrix}, I=\begin{bmatrix} 1 &0 \\ 0 & 1 \end{bmatrix} A=[10?dt1?],B=[0dt?],C=[0?1?],I=[10?01?]
由于系統具有一定的過程噪聲www和測量噪聲vvv。引入噪聲的離散系統狀態方程為:
{x(n)=Ax(n?1)+B(u(n?1)+w′)=Bu(n)+wyv(n)=Cx(n)+v\begin{cases} x(n)=Ax(n-1)+B(u(n-1)+w')=Bu(n)+w\\ y_v(n)=Cx(n)+v \end{cases} {x(n)=Ax(n?1)+B(u(n?1)+w′)=Bu(n)+wyv?(n)=Cx(n)+v?
為加以區別,使用yvy_vyv?在本文中表示含噪聲的真實測量。定義隨機噪聲的方差:
E((w?w ̄)T(w?w ̄))=Q,E((v?v ̄)T(v?v ̄))=RE((w-\overline{w})^T(w-\overline{w}))=Q,E((v-\overline{v})^T(v-\overline{v}))=RE((w?w)T(w?w))=Q,E((v?v)T(v?v))=R
汽車的初始真實位移為0.2,真實速度為0即x0=[0.2,0]Tx_0=[0.2,0]^Tx0?=[0.2,0]T.
2.卡爾曼濾波器
\qquad為方便初學者入門,本文不會從貝葉斯估計的角度證明KF的公式,但會將其應用步驟以盡可能簡單的手段列出。
\qquad嚴格地來說,如果步驟4要成立,需要E((w?w ̄)T(x?x ̄))=0E((w-\overline{w})^T(x-\overline{x}))=0E((w?w)T(x?x))=0即www和vvv下相互獨立,且系統方程中直接傳輸矩陣D=0D=0D=0,但考慮到本文是寫給初學者的,所以在此默認了兩個條件是成立的。
\qquad狀態協方差即為P,測量協方差為CPCTCPC^TCPCT。卡爾曼濾波器算法收斂一般是指的是測量協方差收斂。
3.仿真及效果
在此附上Matlab的仿真代碼
% 模擬要求汽車在一維空間的加速和減速過程 % 控制變量u是汽車的加速度 % 狀態變量x=[p,v],x^hat=[v,a] % w為控制變量的隨機擾動,v為測量的隨機擾動 % Q為w的方差,R為v的方差,假設w與v相互獨立 clear dt = 0.1; % 采樣間隔 N = 100; % 仿真數 Q = diag([0,0.05]); R = 1; A = [1,dt;0,1]; B = [0;dt]; C = [1,0]; P = Q; I = eye(2); x0 = [0.2;0]; % 初始位置和速度 xh0 = [0;0]; % x0的估計 u = ones(1,N); % 加速度恒定 w = sqrt(Q)*randn(2,N); % 控制變量的誤差2*N v = sqrt(R)*randn(1,N); % 測量誤差1*N ye_list = zeros(size(u)); % 估計值 yv_list = zeros(size(u)); % 測量值 y_list = zeros(size(u)); % 實際值 cov_list = zeros(size(u)); % 測量方差 for i = 1:numel(u)xreal = A*x0 + B*u(i); % 真實的狀態變量yreal = C*x0; % 真實的測量x1 = xreal + w(:,i); % 含噪聲的狀態變量yv = yreal + v(i); % 含噪聲的測量xfe = A*xh0 + B*u(i); % 先驗的狀態變量Pfe = A*P*A'+ Q; % 先驗的狀態變量協方差K = Pfe*C'/(C*Pfe*C'+R); % 卡爾曼最優增益xh1 = xfe + K*(yv-C*xfe); % 當前的狀態估計ye = C*xh1;P = (I-K*C)*Pfe;x0 = x1;xh0 = xh1;y_list(i) = yreal;yv_list(i) = yv;ye_list(i) = ye;cov_list(i) = C*P*C'; end ax = (1:N).*dt; figure(1); subplot(2,2,1) plot(ax,y_list,ax,yv_list,ax,ye_list) legend('實際','測量','估計','Location','best') title('汽車的位移') ylabel('位移/m') xlabel('時間/s') subplot(2,2,2) plot(ax,yv_list-y_list,ax,ye_list-y_list) legend('測量','估計','Location','best') title('汽車的位移誤差') ylabel('位移/m') xlabel('時間/s') subplot(2,2,[3,4]) plot(ax,cov_list) legend('測量方差','Location','best') title('測量方差') ylabel('方差/m^2') xlabel('時間/s')本文設定的采樣間隔dt=0.1,x0=[0.2,0]T,x^0=[0,0]Tdt=0.1,x_0=[0.2,0]^T,\hat{x}_0=[0,0]^Tdt=0.1,x0?=[0.2,0]T,x^0?=[0,0]T,注意由于w=Bw′w=Bw'w=Bw′,而B的第一行為0,故QQQ的對角線第一元素必定為0.仿真的效果圖如下:
實際值即狀態方程的輸出y=Cxy=Cxy=Cx,測量值即含噪聲的輸出yv=Cx+vy_v=Cx+vyv?=Cx+v,估計值為卡爾曼濾波器對測量的最優估計值。
\qquad希望本文對您有所幫助,謝謝閱讀!如有異議,歡迎評論區指正!
另外對本話題感興趣的讀者可以繼續閱讀有關適用于非線性系統的擴展卡爾曼濾波算法的介紹。
創作挑戰賽新人創作獎勵來咯,堅持創作打卡瓜分現金大獎總結
以上是生活随笔為你收集整理的【MATLAB】卡尔曼滤波器的原理及仿真(初学者专用)的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: Android实现点击两次返回键退出
- 下一篇: Android开发:利用Activity