# 卡尔曼滤波的细致讲解从一维到多维

## 卡尔曼滤波的细致讲解从一维到多维

python代码实现：

# vim: expandtab:ts=4:sw=4
import numpy as np
import scipy.linalg

"""
Table for the 0.95 quantile of the chi-square distribution with N degrees of
freedom (contains values for N=1, ..., 9). Taken from MATLAB/Octave's chi2inv
function and used as Mahalanobis gating threshold.
"""
chi2inv95 = {
1: 3.8415,
2: 5.9915,
3: 7.8147,
4: 9.4877,
5: 11.070,
6: 12.592,
7: 14.067,
8: 15.507,
9: 16.919}

# 对卡尔曼的理解解释的最细致的博客： https://www.jianshu.com/p/d3b1c3d307e0
# 卡尔曼的本质就是将预测分布于测量分布融合为一个更加精准的分布，预测和测量分布都是高斯分布，假设两个分布相互独立，
# 则融合后的分布为两个分布的乘积，卡尔曼系数k也由此得到。K=d1^2 / (d1^2 + d2^2),这里d1为预测分布的方差，d2为测量分布的
# 方差，这是一维的情况，上升到多维，Kg(k)= P(k|k-1) H’ / (H P(k|k-1) H’ + R)，H为参数矩阵,R为测量噪声协方差，P(k|k-1)是X(k|k-1)对应的协方差
class KalmanFilter(object):
"""
A simple Kalman filter for tracking bounding boxes in image space.

The 8-dimensional state space

x, y, a, h, vx, vy, va, vh

contains the bounding box center position (x, y), aspect ratio a, height h,
and their respective velocities.

Object motion follows a constant velocity model. The bounding box location
(x, y, a, h) is taken as direct observation of the state space (linear
observation model).

"""

def __init__(self):
ndim, dt = 4, 1.

# Create Kalman filter model matrices.
self._motion_mat = np.eye(2 * ndim, 2 * ndim)
for i in range(ndim):
self._motion_mat[i, ndim + i] = dt
self._update_mat = np.eye(ndim, 2 * ndim)

# Motion and observation uncertainty are chosen relative to the current
# state estimate. These weights control the amount of uncertainty in
# the model. This is a bit hacky.
self._std_weight_position = 1. / 20
self._std_weight_velocity = 1. / 160

def initiate(self, measurement):
"""Create track from unassociated measurement.

Parameters
----------
measurement : ndarray
Bounding box coordinates (x, y, a, h) with center position (x, y),
aspect ratio a, and height h.

Returns
-------
(ndarray, ndarray)
Returns the mean vector (8 dimensional) and covariance matrix (8x8
dimensional) of the new track. Unobserved velocities are initialized
to 0 mean.

"""
mean_pos = measurement # (x, y, a, h)
mean_vel = np.zeros_like(mean_pos)
mean = np.r_[mean_pos, mean_vel] # 在行方向上相连，即左右拼接

std = [
2 * self._std_weight_position * measurement[3],
2 * self._std_weight_position * measurement[3],
1e-2,
2 * self._std_weight_position * measurement[3],
10 * self._std_weight_velocity * measurement[3],
10 * self._std_weight_velocity * measurement[3],
1e-5,
10 * self._std_weight_velocity * measurement[3]]
covariance = np.diag(np.square(std))  # 得到对角矩阵
return mean, covariance

def predict(self, mean, covariance):
"""Run Kalman filter prediction step.

Parameters
----------
mean : ndarray
The 8 dimensional mean vector of the object state at the previous
time step.
covariance : ndarray
The 8x8 dimensional covariance matrix of the object state at the
previous time step.

Returns
-------
(ndarray, ndarray)
Returns the mean vector and covariance matrix of the predicted
state. Unobserved velocities are initialized to 0 mean.

"""
std_pos = [ # 位置
self._std_weight_position * mean[3],
self._std_weight_position * mean[3],
1e-2,
self._std_weight_position * mean[3]]
std_vel = [ # 速度
self._std_weight_velocity * mean[3],
self._std_weight_velocity * mean[3],
1e-5,
self._std_weight_velocity * mean[3]]
motion_cov = np.diag(np.square(np.r_[std_pos, std_vel])) # 协方差矩阵,反映各维度两两之间的相关性

mean = np.dot(mean, self._motion_mat.T) # 8-dim mean 点积 对角矩阵motion_mat， 得到8-dim new mean
covariance = np.linalg.multi_dot(( # ACA^T+Mc
self._motion_mat, covariance, self._motion_mat.T)) + motion_cov

return mean, covariance

def project(self, mean, covariance): # 状态分布投射到测量空间
"""Project state distribution to measurement space.

Parameters
----------
mean : ndarray
The state's mean vector (8 dimensional array).
covariance : ndarray
The state's covariance matrix (8x8 dimensional).

Returns
-------
(ndarray, ndarray)
Returns the projected mean and covariance matrix of the given state
estimate.

"""
std = [
self._std_weight_position * mean[3],
self._std_weight_position * mean[3],
1e-1,
self._std_weight_position * mean[3]]
innovation_cov = np.diag(np.square(std))

mean = np.dot(self._update_mat, mean) # update mean: Xk = M Xk-1
covariance = np.linalg.multi_dot(( # update cov: Cov_k = M Cov_k-1 M^T + innovation_cov
self._update_mat, covariance, self._update_mat.T))
return mean, covariance + innovation_cov

def multi_predict(self, mean, covariance):
"""Run Kalman filter prediction step (Vectorized version).

Parameters
----------
mean : ndarray
The Nx8 dimensional mean matrix of the object states at the previous
time step.
covariance : ndarray
The Nx8x8 dimensional covariance matrics of the object states at the
previous time step.

Returns
-------
(ndarray, ndarray)
Returns the mean vector and covariance matrix of the predicted
state. Unobserved velocities are initialized to 0 mean.

"""
std_pos = [
self._std_weight_position * mean[:, 3],
self._std_weight_position * mean[:, 3],
1e-2 * np.ones_like(mean[:, 3]),
self._std_weight_position * mean[:, 3]]
std_vel = [
self._std_weight_velocity * mean[:, 3],
self._std_weight_velocity * mean[:, 3],
1e-5 * np.ones_like(mean[:, 3]),
self._std_weight_velocity * mean[:, 3]]
sqr = np.square(np.r_[std_pos, std_vel]).T

motion_cov = []
for i in range(len(mean)):
motion_cov.append(np.diag(sqr[i]))
motion_cov = np.asarray(motion_cov)

mean = np.dot(mean, self._motion_mat.T) # 预测位置，卡尔曼滤波器5个公式的第一个
left = np.dot(self._motion_mat, covariance).transpose((1,0,2))
covariance = np.dot(left, self._motion_mat.T) + motion_cov  # 预测协方差，卡尔曼滤波器5个公式的第二个

return mean, covariance

def update(self, mean, covariance, measurement):
"""Run Kalman filter correction step.

Parameters
----------
mean : ndarray
The predicted state's mean vector (8 dimensional).
covariance : ndarray
The state's covariance matrix (8x8 dimensional).
measurement : ndarray
The 4 dimensional measurement vector (x, y, a, h), where (x, y)
is the center position, a the aspect ratio, and h the height of the
bounding box.

Returns
-------
(ndarray, ndarray)
Returns the measurement-corrected state distribution.

""" # 卡尔曼解释博客：https://blog.csdn.net/weixin_42905141/article/details/99690849
projected_mean, projected_cov = self.project(mean, covariance) # 将状态分布投射到测量空间

chol_factor, lower = scipy.linalg.cho_factor( # 得到cholesky因子， lower表示是否是下三角函数
projected_cov, lower=True, check_finite=False)
kalman_gain = scipy.linalg.cho_solve(  # 卡尔曼滤波5个公式的第3个公式，计算卡尔曼增益
(chol_factor, lower), np.dot(covariance, self._update_mat.T).T,
check_finite=False).T
innovation = measurement - projected_mean # 测量偏差=测量值-状态值

new_mean = mean + np.dot(innovation, kalman_gain.T)  # 卡尔曼滤波5个公式第4个公式，根据预测和测量值更新位置
new_covariance = covariance - np.linalg.multi_dot((  # 卡尔曼滤波5个公式第5个公式，根据预测和测量值更新协方差
kalman_gain, projected_cov, kalman_gain.T))
return new_mean, new_covariance

def gating_distance(self, mean, covariance, measurements,
only_position=False, metric='maha'):
"""Compute gating distance between state distribution and measurements.

A suitable distance threshold can be obtained from chi2inv95. If
only_position is False, the chi-square distribution has 4 degrees of
freedom, otherwise 2.

Parameters
----------
mean : ndarray
Mean vector over the state distribution (8 dimensional).
covariance : ndarray
Covariance of the state distribution (8x8 dimensional).
measurements : ndarray
An Nx4 dimensional matrix of N measurements, each in
format (x, y, a, h) where (x, y) is the bounding box center
position, a the aspect ratio, and h the height.
only_position : Optional[bool]
If True, distance computation is done with respect to the bounding
box center position only.

Returns
-------
ndarray
Returns an array of length N, where the i-th element contains the
squared Mahalanobis distance between (mean, covariance) and
measurements[i].

"""
mean, covariance = self.project(mean, covariance)
if only_position:
mean, covariance = mean[:2], covariance[:2, :2]
measurements = measurements[:, :2]

d = measurements - mean
if metric == 'gaussian':
return np.sum(d * d, axis=1)
elif metric == 'maha':
cholesky_factor = np.linalg.cholesky(covariance)
z = scipy.linalg.solve_triangular(
cholesky_factor, d.T, lower=True, check_finite=False,
overwrite_b=True)
squared_maha = np.sum(z * z, axis=0)
return squared_maha
else:
raise ValueError('invalid distance metric')


12-18
07-22

04-06 3万+
09-09
08-22
09-15 7万+
03-12 1930
01-16 1万+
07-11 3469