我是一个 dsp 人,我在大学时只做了最低限度的控制理论。在尝试了解状态空间分析和(离散时间)常规卡尔曼滤波器时,我遇到了一些谷歌/维基百科/我的控制理论书无法启发我的问题。如果你原谅我有些冗长的帖子,试图用我自己的话解释卡尔曼在做什么,试图获得关于我误解了哪些部分的反馈(甚至不了解它为什么这样做)。
这个系统可以用差分方程来描述:
x(k) = A*x(k-1) + B*u(k-1) + q(k-1)
y(k) = H*x(k) + D*u(k) + r(k)
在哪里:
- 大写字母代表矩阵
- 副大写代表向量,
- q~高斯(0,Q)。r〜高斯(0,R),
- u(k) 是外部(已知)刺激,
- x(k) 是内部(未知)状态,
- y(k) 是(已知的)系统输出,
- k 是(标量)离散时间,
类似于底部的参考。始终使用伪 Matlab 表示法。所有随机变量都假定为独立高斯。
试图以一种有点手动的方式将卡尔曼事物超级强加于顶部,我们正在对每个时间 (k) 增量的隐藏状态进行两阶段估计:
所以基本上状态和输出被建模为具有缓慢变化的平均值的高斯分布(待估计),顶部有附加的零均值高斯噪声(待取消)?
“最困难的任务是弄清楚如何以状态空间形式制定估计问题。” {1} “在大多数应用中,内部状态比测量的少数“可观察”参数大得多(更多自由度)。但是,通过组合一系列测量,卡尔曼滤波器可以估计整个内部状态。 " {2}
A 矩阵(“状态转换”或“动态”模型)以及扩展到一些隐藏状态维度的真正作用是什么?我想对于无法观察到的变量,它是关于在 x 中强制执行平滑度(连续性),它的导数等等?对于物理系统,人们可能会假设太空火箭不会立即从位置#1 跳到位置#2(无限速度)或从速度#1 跳到速度#2(无限加速度 = 无限力)。是这样吗?对于像股票市场这样的东西,应用类似的规则似乎更抽象,但也许这只是我不是经济学家。
或者更多的是关于能够接受不同的感官输入,这些输入可能会测量(或者)x,或者它的 n 次导数具有一定的准确性,并且(通常)在推导“自然”感官输出或漂移时具有高频噪声你整合了吗?如果传感器融合是卡尔曼滤波器的“事情”,似乎还有其他(更简单的)方法可以完成类似的事情(我的控制理论书提到了互补滤波器,其中一个简单的滤波器组可以让高频信息从一个传感器,来自另一个的低频)。
“预测”阶段包括对下一个状态的预测,仅基于“旧”信息:上一个状态以及动态系统描述(在时间 k-1 以 100 公里/小时的速度沿某个方向行驶的物体具有在时间 k 的预期位置,除非外力、u() 改变其路线)或噪声在实际状态或其测量中造成一些误差:
x(k) = A * x(k-1) + B * u(k-1);
此外,递归预测状态协方差,烘焙过程噪声:
P(k) = A * P(k-1) * A' + Q;
预计用户将呈现状态 P 和状态噪声 Q 的初始协方差估计。P 将在整个过程中更新,而 Q 是静态的。那么应该如何考虑假设的固定加性噪声源与观察到的(样本)状态协方差之间的差异?
在“更新”阶段,将当前测量值集成到模型中并与预测结果进行比较。测量和模型协方差决定了在新测量中放入多少相对权重。
项目状态到输出,创新/测量残差(预测误差?)
e = y(k)-H*x(k);
创新/测量预测协方差
IS = H*P*H' + R;
%% 卡尔曼增益
Kg = P*H'/IS;
%% 更新状态估计,在当前测量残差中加权
x(k) = x(k) + Kg * e;
%% 更新状态协方差
P = P - Kg*IS*Kg';
“卡尔曼滤波器产生系统状态的估计,作为系统预测状态的平均值和使用加权平均值的新测量值。权重的目的是具有更好(即更小)估计不确定性的值是”值得信赖” 更多... 高增益时,滤波器将更多权重放在最近的测量值上,因此响应性更强。低增益时,滤波器更接近模型预测。在极端情况下,高增益接近 1 会导致估计轨迹更加跳跃,而接近 0 的低增益会消除噪声但会降低响应能力。” {2}。有点像指数平滑{4},只有截止频率是自适应的并且基于系统假设而不是临时调整?
但是,即使对于一阶平滑器,也能够根据与预先计算的平均估计的距离来选择每个输入样本的权重,这意味着存在一些“异常值抑制”。即非线性滤波操作,类似于(实物)中值滤波器?
编辑:在@peter-k 的出色回应的帮助下,我对这个问题有了更多的思考。
我仍然觉得我没有很好地掌握这个“A”矩阵的真正作用。它应该关联从时间 k-1 到时间 k 的状态向量(在这种情况下为 2 元素)。但是它在哪里明确编码 x1 是“平均”而 x2 是导数?为什么非对角线“泄漏”项有任何特殊价值?
实验设置:
%adopted from https://github.com/EEA-sensors/ekfukf
%% Create step signal, then corrupt it with noise
N2 = 150;
X(1:N2) = -1;
X((N2+1):2*N2) = 1;
sd = 0.1;%measurement noise standard deviation
rng('default')
Y = X + sd*randn(size(X));
M = [0; 0];%Mean state estimate
P = diag([0.1 2]);%NxN state initial covariance. Jump start by large diagonal
R = sd^2;%Measurement noise covariance (scalar).
H = [1 0];%observe only x
A = eye(2)+[0 0.1; 0 0];%??
Q = diag([1e-4 1e-3]);%tune parameter
执行循环:
for k=1:size(Y,2)
%% Predict
M = A * M;
P = A * P * A' + Q;
%% Update
e = (Y(k)-H*M);
IS = (R + H*P*H');
Kg(:,k) = P*H'/IS;
M = M + Kg(:,k) * e;
P = P - Kg(:,k)*IS*Kg(:,k)';
%% stash variables for plotting
MM(:,k) = M;
PP(:,:,k) = P;
end
讨论:
我曾以某种方式认为卡尔曼增益会因异常噪声(或无噪声信号中的大步长)而波动,因此有些非线性/时变。我在这里看到它收敛到一个常数。观察到一些平滑,并且观察到一些拖尾/过冲。调整 P、Q 和 A 似乎是这里的关键,而且似乎很重要,参见例如 {6}
参考:
{1} https://users.aalto.fi/~ssarkka/course_k2011/pdf/handout1.pdf
{2} https://en.wikipedia.org/wiki/State-space_representation
{3} https://en.wikipedia.org/wiki/Kalman_filter
{4} https://en.wikipedia.org/wiki/Exponential_smoothing
{5} https://stanford.edu/class/ee363/lectures/kf.pdf
{6} https://minds.wisconsin.edu/bitstream/handle/1793/10890/file_1.pdf?sequence=1&isAllowed=y