卡爾曼濾波器(Kalman Filter, KF)是一種最優估計技術,可用於根據角度測量來估算速度。在這裡,我們假設有一個系統,其中角度是可測量的,而速度是我們要估算的隱藏狀態。 --- ## **問題定義** 我們考慮一個一維系統: - 狀態向量:$x = \begin{bmatrix} \theta \\ \dot{\theta} \end{bmatrix}$(角度和角速度) - 狀態轉移模型: $$ x_{k+1} = A x_k + w_k $$ 其中: $$ A = \begin{bmatrix} 1 & T \\ 0 & 1 \end{bmatrix} $$ $T$ 是離散時間間隔,$w_k$ 是過程噪聲。 - 測量方程: $$ z_k = H x_k + v_k $$ 其中: $$ H = \begin{bmatrix} 1 & 0 \end{bmatrix} $$ $z_k$ 是測得的角度,$v_k$ 是測量噪聲。 --- ## **Matlab Code** 以下是使用卡爾曼濾波器來估計角速度的 Matlab 代碼: ```matlab clc; clear; close all; % 參數設定 T = 0.1; % 取樣時間 A = [1 T; 0 1]; % 狀態轉移矩陣 H = [1 0]; % 觀測矩陣 Q = [1e-4 0; 0 1e-2]; % 過程噪聲協方差 R = 1e-2; % 測量噪聲協方差 P = eye(2); % 初始協方差矩陣 x_est = [0; 0]; % 初始狀態 [theta; theta_dot] % 生成模擬數據 N = 100; % 時間步數 theta_true = zeros(1, N); theta_dot_true = zeros(1, N); z = zeros(1, N); for k = 2:N theta_dot_true(k) = 0.1 * sin(0.1 * k); % 真實速度變化 theta_true(k) = theta_true(k-1) + T * theta_dot_true(k); z(k) = theta_true(k) + sqrt(R) * randn; % 加入測量噪聲 end % 卡爾曼濾波估計 x_est_hist = zeros(2, N); for k = 1:N % 預測步驟 x_pred = A * x_est; P_pred = A * P * A' + Q; % 更新步驟 K = P_pred * H' / (H * P_pred * H' + R); x_est = x_pred + K * (z(k) - H * x_pred); P = (eye(2) - K * H) * P_pred; % 儲存結果 x_est_hist(:, k) = x_est; end % 畫圖 figure; subplot(2,1,1); plot(1:N, theta_true, 'g', 1:N, z, 'b.', 1:N, x_est_hist(1,:), 'r'); legend('True Angle', 'Measured Angle', 'Estimated Angle'); xlabel('Time Step'); ylabel('Angle (rad)'); subplot(2,1,2); plot(1:N, theta_dot_true, 'g', 1:N, x_est_hist(2,:), 'r'); legend('True Velocity', 'Estimated Velocity'); xlabel('Time Step'); ylabel('Angular Velocity (rad/s)'); ``` --- ## **代碼解釋** 1. **系統初始化**: - 設定狀態轉移矩陣 $A$、觀測矩陣 $H$ 及協方差矩陣 $Q, R$。 - 初始化狀態 $x_{est}$ 和協方差 $P$。 2. **模擬真實數據**: - 使用 $\theta_{dot}(k) = 0.1 \sin(0.1 k)$ 產生角速度變化。 - 通過數值積分獲得角度,並添加測量噪聲。 3. **卡爾曼濾波迴圈**: - **預測步驟**: - 計算預測狀態 $x_{pred} = A x_{est}$。 - 更新誤差協方差矩陣 $P_{pred}$。 - **更新步驟**: - 計算卡爾曼增益 $K$。 - 更新狀態估計 $x_{est}$。 - 更新誤差協方差矩陣 $P$。 4. **繪圖**: - 比較真實角度、測量角度與估計角度。 - 比較真實角速度與估計角速度。 這樣就能使用卡爾曼濾波器來估算角速度,並在有噪聲的情況下得到較準確的結果! ## 限制式 ### 連續form $$I\ddot{\alpha}<0.6$$$$I\ddot{\beta}<0.6$$ ### Discrete form $$I \frac{\alpha_{k+1} - 2\alpha_k + \alpha_{k-1}}{T^2} < 0.6$$$$I \frac{\beta_{k+1} - 2\beta_k + \beta_{k-1}}{T^2} < 0.6$$