无迹卡尔曼滤波 (Unscented Kalman Filter)

无迹卡尔曼滤波是在卡尔曼滤波的基础上，进行了一系列的改进，得到的一种新的滤波算法。无迹卡尔曼滤波中不需要假设状态变量的变化是线性的，也不用假设观测变量和状态变量的关系是线性的。所以需要我们指定变化的函数。但是，无迹卡尔曼滤波也不是没有限制条件。卡尔曼滤波中观测变量和测量变量都需要服从正态分布，这个条件在这里也是需要的。关于卡尔曼滤波，可以参考之前的博客，这里不再赘述。

无迹变换（Unscented Transform）

根据卡尔曼滤波的计算公式，我们知道，卡尔曼滤波会根据之前的状态变量$x_{t-1}$来预测当前的状态变量$x_t$。然后根据观测变量$y_t$来预测当前的状态变量$x_t$。然后根据两者的协方差矩阵来确定卡尔曼增益$K$，用来判断当前的预测中，根据状态变量的预测的可信度大小以及根据测量变量预测的可信度大小。然后找到一个最优解，即卡尔曼滤波的预测结果。但是对于无迹卡尔曼滤波来讲，协方差矩阵的计算并不容易。因为状态变量的变化以及观测变量和状态变量的关系不是线性的。无法通过矩阵计算直接得到协方差矩阵。那么如何去计算当前的协方差矩阵。目前有两种解决方案：

首先是利用蒙特卡罗方法模拟大量数据点，得到转换后数据的分布均值和方差，然后计算协方差矩阵。这种方法计算的精确程度和模拟数据点数量有关。但是大量的模拟数据点会增加计算量。所以不是一个好的选择。然后就是无迹变换，这种方法能通过很少的采样点数据，计算出转换后数据的分布均值和方差，从而计算协方差矩阵。

对于一组数据 x = { x 1 , x 2 , . . . , x n } x=\{x_1, x_2, ..., x_n\} x={x1​,x2​,...,xn​}，经过一个变换$f$后，得到了另外一组数据 y = { y 1 , y 2 , . . . , y n } y=\{y_1, y_2, ..., y_n\} y={y1​,y2​,...,yn​}。$x\sim N(\overline{x},Q)$，且$y\sim N(\overline{y},R)$。

UT变换需要用到以下几个参数：

首先是采样点的权值分配参数：
$$\begin{gathered} w_m^0=\frac{\lambda }{n+\lambda } \\ {w}_c^0=\frac{\lambda }{n+\lambda }+(1-{\alpha }^2+\beta ) \\ {w}_m^i=w_c^i=\frac{n}{2(n+\lambda )},i=1,...,2n \end{gathered}$$

其中$n$表示 $x$的维度，$\lambda$是一个常数，一般可以表示为$\lambda ={\alpha }^2(n+k)-n$。其中$k\ge 0$，且一般取值为$k=3-n$，$\alpha \in (0,1]$，$\alpha$主要用来调节高阶项对模型的影响，一般取值$\alpha =0.01$。对于高斯分布， $\beta =2$时效果是最优的。

然后是无迹变换的采样，一般无迹变换的采样点个数为$2n+1$。采样过程如下：
$$\begin{gathered} X_0={\overline{x}}_t \\ cho=chol((n+k)P_t) \\ {X}_i=\overline{x}+cho[:,i],i=1,...,n \\ {X}_i=\overline{x}-cho[:,i],i=n+1,...,2n \end{gathered}$$
这里的$chol$表示矩阵的cholesky分解。这里取分解后的下三角矩阵值。然后对于没一个采样点，分别加上或者减去cholesky分解矩阵的第$i$列。

然后根据采样点计算经过转换函数$f$后的数据分布的均值和方差。
$$\begin{gathered} {\overline{y}}_t=\sum _{i=0}^{2n}{w}_i^{m}f(X_i) \\ {P}_t^y=\sum _{i=0}^{2n}(w_i^c(f(Xi)-{\overline{y}}_t)(f(Xi)-{\overline{y}}_t{)}^T \end{gathered}$$

无迹卡尔曼滤波（Unscented Kalman Filter）

根据上述无迹变换的方法，对卡尔曼滤波做一些改进，就得到了无迹卡尔曼滤波，计算方法如下：

1. 给定数据集 D = { X , Y } D=\{X,Y\} D={X,Y}。其中 X = { x 1 , x 2 , . . . , x t } X=\{x_1, x_2, ..., x_t\} X={x1​,x2​,...,xt​}， Y = { y 1 , y 2 , . . . , y t } Y=\{y_1, y_2, ..., y_t\} Y={y1​,y2​,...,yt​}以及状态变量$X$的转换函数$f$和根据测量变量$Y$计算状态变量的转换函数$g$。其中$\delta$和$ϵ$都是服从高斯分布的0均值噪声，这部分也可以放在转换函数$f$和$g$里面，效果是一样的。
   $$\begin{gathered} x_t=f(x_{t-1})+{\delta }_t,\delta \sim N(0,Q) \\ {y}_t=g(x_t)+{ϵ}_t,ϵ\sim N(0,R) \end{gathered}$$

2. 根据UT变换，计算状态转换函数的估计值$x_t^{\prime }$，以及协方差矩阵${P^{\prime }}_t^{xx}$。UT变换中的输入变量的均值用$x_{t-1}$表示。状态变量的协方差在初始化时可以用对角矩阵表示
   $$\begin{gathered} X_0=x_{t-1} \\ cho=chol((n+k)P_{t-1}^{xx}) \\ {X}_i=x_{t-1}+cho[:,i],i=1,...,n \\ {X}_i=x_{t-1}-cho[:,i],i=n+1,...,2n \\ {x}_t^{\prime }=\sum _{i=0}^{2n}{w}_i^{m}f(X_i) \\ {P^{\prime }}_t^{xx}=\sum _{i=0}^{2n}(w_i^c(f(Xi)-x_t^{\prime })(f(Xi)-x_t^{\prime }{)}^T+Q \end{gathered}$$

3. 根据UT变换，计算测量变量的转换函数的估计值$y_t^{\prime }$，以及协方差矩阵$P_t^{yx}$。UT变换中的输入变量的均值用$x_t^{\prime }$表示
   $$\begin{gathered} X_0=x_t^{\prime } \\ cho=chol((n+k){P^{\prime }}_t^{xx}) \\ {X}_i=x_t^{\prime }+cho[:,i],i=1,...,n \\ {X}_i=x_t^{\prime }-cho[:,i],i=n+1,...,2n \\ {y}_t^{\prime }=\sum _{i=0}^{2n}{w}_i^{m}g(X_i) \\ {P}_t^{yy}=\sum _{i=0}^{2n}{w}_i^c(g(Xi)-y_t^{\prime })(g(Xi)-y_t^{\prime }{)}^T+R \end{gathered}$$

4. 根据测量变量预测值$y_t‘$以及观测变量预测值$x_t^{\prime }$，计算测量变量和观测变量的协方差矩阵$P_t^{xy}$。注意，这里的$f(X_i)$和$g(X_i)$中的输入$X_i$表示的时不同的UT变换的输入。可以参考上述步骤
   $$P_t^{xy}=\sum _{i=0}^{2n}{w}_i^c(f(X_i)-x_t^{\prime })(g(X_i)-y_t^{\prime }{)}^T$$

5. 根据协方差矩阵$P_t^{xx}$和$P_t^{yx}$，计算卡尔曼增益$K$
   $$K=P_t^{xy}(P_t^{yy}{)}^{-1}$$

6. 根据卡尔曼增益$K$计算估计值$x_t$
   $$x_t=x_t^{\prime }+K(y_t-y_t^{\prime })$$

7. 更新状态变量的协方差矩阵$P_t^{xx}$
   $$P_t^{xx}={P^{\prime }}_t^{xx}-K{P}_t^{yy}{K}^T$$

实例

用阿里天池数据集中的发动机数据，这组数据主要包括37个发动机的数据和1个target的数据。需要建立模型根据发动机的数据预测target的数据。这里状态变量和测量变量的函数都用线性函数表示。用CC和RMSE作为评价指标，得到的结果如下图所示。可以看出CC=0.9，RMSE=0.4，属于比较好的解码效果了。