Tracking a Cruise Missile with a Kalman Filter When Prediction and Measurement Rates Differ

A cruise missile is assumed to fly in a straight line toward a target, while a surveillance radar located at the target measures only the missile’s range.

The setup is as follows:

• Initial distance: 100 km
• Nominal speed: 300 m/s
• Motion: approximately constant velocity, with aerodynamic disturbance modeled as zero-mean white acceleration noise
• Acceleration noise intensity: q = 0.05 m^2/s^3
• Radar measurement rate: 2 Hz
• Radar range error: zero-mean white noise with standard deviation 50 m
• Time update rate: 10 Hz

The goal is to use a Kalman filter to estimate the missile trajectory, and then compare that result with a second case in which both the time update and measurement update are performed at 2 Hz.

In this problem, the key point is that the filter does not receive measurements as often as it propagates the state. With a 10 Hz prediction step and a 2 Hz radar sampling rate, the filter performs five time updates for every one measurement update.

State model

The state contains:

• range
• velocity

The missile is modeled with nearly constant velocity, while random acceleration enters as process noise. Radar measures only range, so velocity must be inferred by the filter from the dynamic model and repeated observations.

Using the notation in the code:

• Continuous-time state matrix:

[ A = \begin{bmatrix}0 & -1 \ 0 & 0\end{bmatrix} ]

• Measurement matrix:

[ H = [1,\;0] ]

• Noise input matrix:

[ \Gamma = \begin{bmatrix}0 \ 1\end{bmatrix} ]

The discrete transition matrix is obtained with a simple first-order discretization:

[ \Phi = I + A T_{update} ]

Simulation idea

The simulation first generates a reference trajectory from the motion model, then adds radar noise to create the measured range sequence. After that, two Kalman filtering schemes are run:

1. 10 Hz time update + 2 Hz measurement update
2. 2 Hz time update + 2 Hz measurement update

The first scheme uses the model more frequently between radar observations. When no new measurement is available, the filter only performs the prediction step and carries the predicted state forward. Once the measurement time arrives, the standard Kalman correction is applied.

MATLAB code

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 功能描述：导弹运动Kalman滤波程序
% 课次：卡尔曼滤波与组合导航 第二次课程
% 时间：2021/5/4
% 作者：Lei Lie
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;clc;close all;
%% 01 初始化参数
T = 300;                   % 仿真时长

T_mea = 1/2;               % 量测采样时间
T_update = 1/10;           % 时间更新间隔
Q = 0.05*T_update;         % 过程噪声
R = 50;                    % 量测噪声
W = sqrt(Q)*randn(1,T);
V = sqrt(R)*randn(1,T);
P0 = diag([10^2,1^2]);

% 系统矩阵
A = [0 -1;0 0];            % 状态矩阵
I = eye(2);
Phi = I + A*T_update;      % 离散化
H = [1,0];                 % 量测矩阵
Gamma = [0;1];

% 初始化
nS = 2; nZ = 1;
count = 1;                % 用于计数
xState = zeros(nS,T);
zMea = zeros(nZ,T);
xKF_10HZ = zeros(nS,T);
xKF_2HZ = zeros(nS,T);
xKF_pre = zeros(nS,T);
xState(:,1) = [100000;300];
zMea(:,1) = H*xState(:,1);
xKF_10HZ(:,1) = xState(:,1);
xKF_2HZ(:,1) = xState(:,1);
%% 02 用模型模拟真实状态
for t = 2:T
    xState(:,t) = Phi*xState(:,t-1) + Gamma*W(:,t);
    zMea(:,t) = H*xState(:,t) + V(t);
end
%% 03-1 Kalman滤波（时间更新为10Hz,量测频率为2Hz）
for t = 2:T
    % 时间更新（时间更新为10Hz）
    xKF_pre(:,t) = Phi*xKF_10HZ(:,t-1);
    P_pre = Phi*P0*Phi' + Gamma*Q*Gamma';
    % 量测更新（量测频率为2Hz）
    if (mod(count,T_mea/T_update)==0)       % 时间更新5次，量测1次
        K = P_pre*H'*pinv(H*P_pre*H'+R);
        xKF_10HZ(:,t) = xKF_pre(:,t) + K*(zMea(:,t)-H*xKF_pre(:,t));
        P0 = (I-K*H)*P_pre;
        count = 1;                          % 计数归0
    else
        xKF_10HZ(:,t) = xKF_pre(:,t);       % 将上一拍的值传给滤波器
        P0 = P_pre;
        count = count + 1;                  % 计数加1
    end
end

%% 03-2 Kalman滤波（时间更新为2Hz,量测频率为2Hz）
for t = 2:T
    % 时间更新（时间更新为2Hz）
    xKF_pre(:,t) = Phi*xKF_2HZ(:,t-1);
    P_pre = Phi*P0*Phi' + Gamma*Q*Gamma';
    % 量测更新（量测频率为2Hz）
    K = P_pre*H'*pinv(H*P_pre*H'+R);
    xKF_2HZ(:,t) = xKF_pre(:,t) + K*(zMea(:,t)-H*xKF_pre(:,t));
    P0 = (I-K*H)*P_pre;
end

%% 04 画图
tPlot = 1:T;
FigWin1=figure('position',[300 300 550 450],'Color',[0.8 0.8 0.8],...
   'Name','01-量测距离误差与估计距离误差的比较','NumberTitle','off');hold on;box on;
plot(tPlot,xState(1,:)-zMea(1,:),'-b','LineWidth',1.5);hold on;
plot(tPlot,xState(1,:)-xKF_10HZ(1,:),'-r','LineWidth',1.5);hold on;
plot(tPlot,xState(1,:)-xKF_2HZ(1,:),'-g','LineWidth',1.5);hold on;
xlabel('时间 t/s');ylabel('误差距离 m');
legend('量测距离误差','估计距离误差(10HZ)','估计距离误差(2HZ)');
title('量测距离误差与估计距离误差的比较');
% 保存图片
saveas(gcf,'01-量测距离误差与估计距离误差的比较.png');

FigWin2=figure('position',[850 300 550 450],'Color',[0.8 0.8 0.8],...
   'Name','02-真实速度与估计速度的比较','NumberTitle','off');hold on;box on;
plot(tPlot,xState(2,:)-xKF_10HZ(2,:),'-r','LineWidth',1.5);hold on;
plot(tPlot,xState(2,:)-xKF_2HZ(2,:),'-g','LineWidth',1.5);hold on;
xlabel('时间 t/s');ylabel('误差速度 m/s');
legend('估计速度(10HZ)','估计速度(2HZ)');
title('估计速度的比较');
% 保存图片
saveas(gcf,'02-真实速度与估计速度的比较.png');

FigWin3=figure('position',[300 200 550 450],'Color',[0.8 0.8 0.8],...
   'Name','03-导弹真实距离、雷达量测距离与估计距离','NumberTitle','off');hold on;box on;
plot(tPlot,xState(1,:),'-k','LineWidth',1.5);hold on;
plot(tPlot,zMea(1,:),'-b','LineWidth',1.5);hold on;
plot(tPlot,xKF_10HZ(1,:),'-r','LineWidth',1.5);hold on;
plot(tPlot,xKF_2HZ(1,:),'-g','LineWidth',1.5);hold on;
xlabel('时间 t/s');ylabel('距离 m');
legend('导弹真实距离','雷达量测轨迹','导弹估计距离(10HZ)','导弹估计距离(2HZ)');
title('导弹真实距离、雷达量测距离与估计距离的比较');
axes('position',[0.25,0.25,0.25,0.25]);
hold on;
plot(tPlot,xState(1,:),'-k','LineWidth',1.5);hold on;
plot(tPlot,zMea(1,:),'-b','LineWidth',1.5);hold on;
plot(tPlot,xKF_10HZ(1,:),'-r','LineWidth',1.5);hold on;
plot(tPlot,xKF_2HZ(1,:),'-g','LineWidth',1.5);hold on;
xlim([T/2,T/2+0.5]);
% 保存图片
saveas(gcf,'03-导弹真实距离、雷达量测距离与估计距离.png');

FigWin4=figure('position',[850 200 550 450],'Color',[0.8 0.8 0.8],...
   'Name','04-导弹真实速度与估计速度的比较','NumberTitle','off');hold on;box on;
plot(tPlot,xState(2,:),'-b','LineWidth',1.5);hold on;
plot(tPlot,xKF_10HZ(2,:),'-r','LineWidth',1.5);hold on;
plot(tPlot,xKF_2HZ(2,:),'-g','LineWidth',1.5);hold on;
xlabel('时间 t/s');ylabel('速度 m/s');
legend('真实速度','估计速度(10HZ)','估计速度(2HZ)');
title('真实速度与估计速度的比较');
% 保存图片
saveas(gcf,'04-真实速度与估计速度的比较.png');

FigWin5=figure('position',[300 100 550 450],'Color',[0.8 0.8 0.8],...
   'Name','05-估计距离与估计速度','NumberTitle','off');hold on;box on;
subplot(211);
plot(tPlot,xKF_10HZ(1,:),'-r','LineWidth',1.5);hold on;
plot(tPlot,xKF_2HZ(1,:),'-b','LineWidth',1.5);hold on;
legend('10HZ估计值','2HZ估计值');
xlabel('时间 t/s');ylabel('估计距离 m');
title('距离估计值与速度估计值');
subplot(212);
plot(tPlot,xKF_10HZ(2,:),'-r','LineWidth',1.5);hold on;
plot(tPlot,xKF_2HZ(2,:),'-b','LineWidth',1.5);hold on;
xlabel('时间 t/s');ylabel('估计速度 m/s');
% 保存图片
saveas(gcf,'05-估计距离与估计速度.png');

What the code is doing

A few details in the implementation are worth noting.

The simulation length is set to 300 steps. The radar sampling interval is 0.5 s and the time update interval is 0.1 s. Process noise is scaled by the update step as Q = 0.05*T_update, while measurement noise is set as R = 50 in the code.

The true state is initialized as:

• range = 100000 m
• velocity = 300 m/s

During simulation, the state evolves through the discrete model, and each radar observation is generated by adding measurement noise to the true range.

For the 10 Hz / 2 Hz filter, a counter is used to decide when to perform the measurement update. Most of the time, the algorithm only predicts:

• propagate state with Phi
• propagate covariance with Phi*P0*Phi' + Gamma*Q*Gamma'

Every fifth prediction step, it performs the Kalman correction using the latest range measurement.

For the 2 Hz / 2 Hz case, the filter updates with a measurement every step, so prediction and correction occur at the same rate.

Simulation plots

The following figures compare the raw radar measurements, the true missile trajectory, and the two filter outputs.

Comparison of the two update strategies

The comparison focuses on whether increasing the time update frequency improves tracking when the measurement frequency stays the same.

With 10 Hz prediction and 2 Hz measurement, the filter uses the motion model more often between radar samples. That gives it more frequent state propagation and typically produces a smoother estimate, especially for velocity, which is not directly measured.

With both updates running at 2 Hz, the filter still works, but it relies on a coarser propagation cycle. Since the system is evolving continuously while measurements arrive only every 0.5 s, a lower prediction rate means less frequent use of the dynamic model between observations.

From the plotted results, the filtered estimates are clearly better than the noisy range measurements, and the 10 Hz time update case shows an advantage over the 2 Hz time update case, particularly in the quality of the state estimation between measurement instants.

This example highlights a practical point in Kalman filtering: even when sensors report slowly, raising the prediction rate can still improve tracking performance, as long as the system model is reasonable and the process noise is handled properly.