目前, 卡尔曼滤波、粒子滤波、基于机器学习的滤波及其改进形式已被广泛应用在目标跟踪、参数估计和状态预测等领域^[1-4]。对于基于单模型的滤波方法而言, 当目标机动十分复杂时, 会出现滤波模型与目标机动模型不匹配的情况, 从而导致滤波精度大幅降低^[5]。

针对这一问题, 交互式多模型(interacting multiple model, IMM)^[6-8]展现出更加优越的性能。近些年, IMM算法得到了很多学者的关注, 很多学者在结构和参数等层面对IMM算法进行了改进。文献[9]基于恒速和当前统计模型设计了IMM算法, 首先利用最小二乘法估计当前统计模型的平均速度, 然后将当前统计模型应用于IMM算法。该方法提高了模型精度, 从而提高了滤波精度。文献[10]从多个方面对IMM算法进行了改进, 包括采用改进的卡尔曼滤波器作为子滤波器、不同模型之间的非对称状态估计和基于熵的模型概率更新公式。文献[7, 11]提出了IMM的一种替代方法, 模型集中的模型由匀加速模型构成, 降低了模型集的复杂程度。在此基础上, 文献[12]提出了一种自适应IMM算法, 首先利用滤波器对目标的加速度进行估计, 然后选取估计加速度附近的值构建模型。该方法可以减少模型集中模型的数量, 在降低计算量的同时提高了模型精度。文献[13]基于二阶马尔可夫链提出了一种二阶IMM算法, 该算法利用了更多的先验信息, 提高了滤波精度。由于上述IMM模型集中的模型种类和数量不变, 所以又称其为固定结构IMM(fixed structure IMM, FSIMM)。

为了避免因模型不匹配而造成的精度误差问题, 在使用IMM算法时, 应使用尽可能多的模型覆盖目标机动模型。但值得注意的是, 单个模型集中模型过多, 同样会降低滤波精度^[14]。在此背景下, 变结构IMM(variable structure IMM, VSIMM)应运而生。VSIMM经过不断的发展和改进, 大致可以分为4类: 模型组切换(model group switching, MGS)、可能模型集(likely mode set, LMS)、期望模型增强(expected mode augmentation, EMA)和自适应网格(adaptive grid, AG)^[15]。其中, MGS将模型集分为模型子集, 一个时间步长只选择一个模型子集进行估计, 模型子集之间根据模型子集转移概率进行切换^[16]。LMS将模型分为3种类型: 不可能的、重要的和主要的模型, 在每个时间步长里, 用于估计的模型子集由主要的模型和接近主要的模型构成^[17-18]。与MGS类似, EMA将一个大的模型集分成小的模型子集, 然后计算下一个时间步长所有模型子集的概率, 选择概率最大的模型子集用于估计^[19]。AG算法结合图理论, 将所有的模型构成一个网格, 利用先验信息和当前数据得到一个局部细化网格, 细化网格中的模型构成候选模型子集, 然后根据规则选择模型, 构成下一时刻用于滤波的模型子集^[20]。

值得注意的是, 文献[19]所提到的EMA方法在计算似然函数和模型子集概率时, 对公式进行了近似。针对这一问题, 本文首先设计了一种新的VSIMM(novel VSIMM, NVSIMM)滤波算法, 给出了精确的数学模型。另外, 在前向NVSIMM滤波的基础上, 对数据进行平滑处理, 设计了VSIMM平滑(VSIMM smoothing, VSIMMS)算法, 前向NVSIMM和后向VSIMMS相结合即为VSIMM滤波和平滑(VSIMM filtering and smoothing, VSIMMFS)算法。

本文安排如下: 第1节对滤波和平滑问题进行了描述, 并给出了IMM滤波算法和IMM平滑算法的数学模型; 第2节设计了VSIMMFS算法, 包括前向NVSIMM滤波和后向VSIMM平滑两个部分; 第3节对本文所提算法的有效性进行仿真验证和讨论; 最后, 第4节给出了本文的结论。

假设目标可能有r个运动模型, 模型集记为Ω={M¹, M², …, M^r}, 模型之间的转换概率矩阵为

(1) $\boldsymbol{P}=\left(\begin{array}{ccc}p^{11} & \cdots & p^{1 r} \\\vdots & \ddots & \vdots \\p^{r 1} & \cdots & p^{r r}\end{array}\right)$

式中: p^ij(1≤i≤r, 1≤j≤r)为模型i到模型j的转移概率。

系统离散化状态方程如下:

(2) $\boldsymbol{x}_k^j=\boldsymbol{F}^j \boldsymbol{x}_{k-1}^j+\boldsymbol{w}_{k-1}^j, j=1, 2, \cdots, r$

式中: x_k^j为系统的状态向量; F^j为模型j的状态转移矩阵; w_k－1^j是均值为0的高斯白噪声, 其协方差矩阵为Q_k－1^j。

模型j的量测方程为

(3) $\boldsymbol{z}_k=\boldsymbol{H}^j \boldsymbol{x}_k+v_k^j$

式中: z_k为量测向量; H^j为模型j的测量矩阵; v_k^j是均值为0的高斯白噪声, 其协方差矩阵为R_k^j。

基于上述状态方程、量测方程和贝叶斯理论, IMM滤波算法通过利用量测集合Z^k={z₁, z₂, …, z_k}和模型集Ω估计系统在k时刻的状态x_k。而平滑的目的是根据当前所有的观测值Z^k={z₁, z₂, …, z_k}, 估计t时刻的系统状态$\hat{\boldsymbol{x}}_{t \mid k}$, 其中1≤t≤k－1。

IMM算法框图如图 1所示。

IMM算法一般可以分为以下4个步骤^[6-7]。

步骤1  前向输入交互

(4) $\hat{\boldsymbol{x}}_{k-1 \mid k-1}^{o j}=\sum\limits_{i=1}^r \hat{\boldsymbol{x}}_{k-1 \mid k-1}^i \boldsymbol{u}_{k-1 \mid k-1}^{i j}$

(5) $\begin{gathered}\boldsymbol{P}_{k-1 \mid k-1}^{o j}=\sum\limits_{i=1}^r \boldsymbol{u}_{k-1 \mid k-1}^{i j}\left\{\boldsymbol{P}_{k-1 \mid k-1}^i+\right. \\\left.\left[\hat{\boldsymbol{x}}_{k-1 \mid k-1}^i-\hat{\boldsymbol{x}}_{k-1 \mid k-1}^{o j}\right]\left[\hat{\boldsymbol{x}}_{k-1 \mid k-1}^i-\hat{\boldsymbol{x}}_{k-1 \mid k-1}^{o j}\right]^{\mathrm{T}}\right\}\end{gathered}$

式中: $\hat{\boldsymbol{x}}_{k-1 \mid k-1}^{o j}$和$\boldsymbol{P}_{k-1 \mid k-1}^{o j}$分别为模型j的混合状态估计和对应的协方差矩阵；$\hat{\boldsymbol{x}}_k^i-1 \mid k-1$和$\boldsymbol{P}_{k-1 \mid k-1}^i$为模型i的卡尔曼估计和对应的协方差矩阵；u_k－1|k－1^ij为混合概率, 其计算公式为

(6) $u_{k-1 \mid k-1}^{i j}=\frac{p^{i j} u_{k-1}^i}{c_k^j}$

式中: u_k－1ⁱ为模型i的概率, 预测概率$c_k^j=\sum_{i=1}^r p^{i j} u_{k-1}^i=$$p\left(M_k^j \mid \boldsymbol{Z}^{k-1}\right)$为归一化因子。

步骤2  滤波

分别基于每一个模型进行卡尔曼滤波:

步骤2.1  预测状态

(7) $\hat{\boldsymbol{x}}_{k \mid k-1}^{o j}=\boldsymbol{F}_{k \mid k-1}^j \hat{\boldsymbol{x}}_{k-1 \mid k-1}^{o j}$

步骤2.2  预测协方差矩阵

(8) $\boldsymbol{P}_{k \mid k-1}^{o j}=\boldsymbol{F}_{k \mid k-1}^j \boldsymbol{P}_{k-1 \mid k-1}^{o j} \boldsymbol{F}_{k \mid k-1}^j{ }^{\mathrm{T}}+\boldsymbol{Q}_{k-1}^j$

步骤2.3  计算卡尔曼增益

(9) $\boldsymbol{K}_k^{o j}=\boldsymbol{P}_{k \mid k-1}^{o j} \boldsymbol{H}_k^{j \mathrm{~T}}\left[\boldsymbol{H}_k^j \boldsymbol{P}_{k \mid k-1}^{o j} \boldsymbol{H}_k^j+\boldsymbol{R}_k^j\right]^{-1}$

步骤2.4  滤波估计

(10) $\hat{\boldsymbol{x}}_{k \mid k}^j=\hat{\boldsymbol{x}}_{k \mid k-1}^{o j}+\boldsymbol{K}_k^{o j}\left[z_k-\boldsymbol{H}_k^j \hat{\boldsymbol{x}}_{k \mid k-1}^{o j}\right]$

步骤2.5  滤波值的协方差

(11) $\boldsymbol{P}_{k \mid k}^j=\left[\boldsymbol{I}-\boldsymbol{K}_k^{o j} \boldsymbol{H}_k^j\right] \boldsymbol{P}_{k \mid k-1}^{o j}$

步骤3  模型概率更新

模型概率更新公式为

(12) $u_k^j=p\left\{M_k^j \mid \boldsymbol{Z}^k\right\}=\frac{p\left\{z_k \mid M_k^j, \boldsymbol{Z}^{k-1}\right\} p\left\{M_k^j \mid \boldsymbol{Z}^{k-1}\right\}}{p\left(z_k \mid \boldsymbol{Z}^{k-1}\right)}$

式中: 分母p(z_k|Z^k－1)为归一化因子, 对分子积分即可得到, p{M_k^j|Z^k－1}=c_k^j在步骤1中已经求出, p{z_k|M_k^j, Z^k－1}服从高斯分布, 其具体形式为

(13) $z_k \sim N\left(\boldsymbol{H}_k^j \hat{\boldsymbol{x}}_{k \mid k-1}^{o j}, \boldsymbol{S}_k^j\right)$

式中: $\boldsymbol{H}_k^j \hat{\boldsymbol{x}}_{k \mid k-1}^{o j}$为均值; $\boldsymbol{S}_k^j=\boldsymbol{H}_k^j \boldsymbol{P}_{k \mid k-1}^{o j} \boldsymbol{H}_k^{j \mathrm{~T}}+\boldsymbol{R}_k^j$为方差。p{z_k|M_k^j, Z^k－1}的概率密度为

(14) $\begin{gathered}\mathit{\Lambda}_k^j=p\left\{z_k \mid M_k^j, \boldsymbol{Z}^{k-1}\right\}= \\\frac{1}{(2 \pi)^{n / 2}\left|\boldsymbol{S}_k^j\right|^{1 / 2}} \exp \left\{-\frac{1}{2}\left(\gamma_k^j\right)^{\mathrm{T}}\left(\boldsymbol{S}_k^j\right)^{-1} \gamma_k^j\right\}\end{gathered}$

式中: $\gamma_k^j=z_k-\boldsymbol{H}_k^j \hat{\boldsymbol{x}}_{k-1 \mid k-1}^{o j}$。

步骤4  输出

IMM滤波的状态估计和协方差估计分别为

(15) $\hat{\boldsymbol{x}}_{k \mid k}=\sum\limits_{j=1}^r \hat{\boldsymbol{x}}_{k \mid k}^j u_k^j$

(16) $\boldsymbol{P}_{k \mid k}=\sum\limits_{j=1}^r u_k^j\left\{\boldsymbol{P}_{k \mid k}^j+\left[\hat{\boldsymbol{x}}_{k \mid k}^j-\hat{\boldsymbol{x}}_{k \mid k}\right]\left[\hat{\boldsymbol{x}}_{k \mid k}^j-\hat{\boldsymbol{x}}_{k \mid k}\right]^{\mathrm{T}}\right\}$

在前向IMM滤波的基础上, 对数据进一步进行平滑处理, 得到后向IMM平滑算法, 其算法框图如图 2所示, 包括以下几个步骤^[21]。

步骤1  后向状态交互

后向混合状态和对应的协方差为

(17) $\hat{\boldsymbol{x}}_{t+1 \mid k}^{o j}=\sum\limits_{i=1}^N u_{t+1 \mid k}^{i \mid j} \hat{\boldsymbol{x}}_{t+1 \mid k}^i$

(18) $\boldsymbol{P}_{t+1 \mid k}^{o j}=\sum\limits_{i=1}^N u_{t+1 \mid k}^{i \mid j}\left\{\boldsymbol{P}_{t+1 \mid k}^i+\left[\hat{\boldsymbol{x}}_{t+1 \mid k}^i-\hat{\boldsymbol{x}}_{t+1 \mid k}^{o j}\right]\left[\hat{\boldsymbol{x}}_{t+1 \mid k}^i-\hat{\boldsymbol{x}}_{t+1 \mid k}^{o j}\right]^{\mathrm{T}}\right\}$

式中: t+1=k时刻$\hat{\boldsymbol{x}}_{t+1 \mid k}^i$ 和$\boldsymbol{P}_{t+1 \mid k}^i$的值是已知的, t+1 < k的值可以在步骤2中不断地递推得到; u_t+1|k^i|j为后向混合概率, 其计算公式为

(19) $u_{t+1 \mid k}^{i \mid j}=\frac{1}{d_j} b_{i j} u_{t+1 \mid k}^i$

式中: $d_j=\sum_{l=1}^N b_{l j} u_{t+1 \mid k}^l$为归一化因子, 同样t+1=k时刻u_t+1|k^l的值是已知的, t+1 < k的值可以在步骤3中不断地递推得到; b_ij为后向模型转移概率, 其计算公式为

(20) $b_{i j}=\frac{1}{e_i} p_{j i} u_{t \mid t}^j$

式中: $e_i=\sum_{l=1}^N p_{l i} u_{t \mid t}^l$为归一化因子。

步骤2  卡尔曼平滑^[22]

平滑值和对应的协方差为

(21) $\hat{\boldsymbol{x}}_{t \mid k}^j=\hat{\boldsymbol{x}}_{t \mid t}^j+\boldsymbol{A}_{t \mid k}^j\left(\hat{\boldsymbol{x}}_{t+1 \mid k}^j-\hat{\boldsymbol{x}}_{t+1 \mid k}^{o j}\right)$

(22) $\boldsymbol{P}_{t \mid k}^j=\boldsymbol{P}_{t \mid t}^j+\boldsymbol{A}_{t \mid k}^j\left(\boldsymbol{P}_{t+1 \mid k}^j-\boldsymbol{P}_{t+1 \mid k}^{o j}\right) \boldsymbol{A}_{t \mid k}^j{ }^{\mathrm{T}}$

式中: $\boldsymbol{A}_{t \mid k}^j=\boldsymbol{P}_{t \mid t}^j \boldsymbol{F}_t^{j \mathrm{~T}}\left(\boldsymbol{P}_{t+1 \mid k}^{o j}\right)^{-1}$。

步骤3  平滑模型概率更新

平滑后的模型概率计算公式为

(23) $u_{t \mid k}^j=\sum\limits_{i=1}^N b_{i j} u_{t+1 \mid k}^i$

步骤4  IMM平滑状态交互

后向IMM平滑的状态输出和对应的协方差为

(24) $\hat{\boldsymbol{x}}_{t \mid k}=\sum\limits_{j=1}^N u_{t \mid k}^j \hat{\boldsymbol{x}}_{t \mid k}^j$

(25) $\boldsymbol{P}_{t \mid k}=\sum\limits_{j=1}^N u_{t \mid k}^j\left\{\boldsymbol{P}_{t \mid k}^j+\left[\hat{\boldsymbol{x}}_{t \mid k}^j-\hat{\boldsymbol{x}}_{t \mid k}\right]\left[\hat{\boldsymbol{x}}_{t \mid k}^j-\hat{\boldsymbol{x}}_{t \mid k}\right]^{\mathrm{T}}\right\}$

VSIMMFS算法包括前向NVSIMM滤波和后向VSIMM平滑两个部分, 这两个部分分别在后续两个小节中展开。

为了避免模型集中模型数量过多, 减少模型误差, 本节设计了NVSIMM算法。NVSIMM滤波算法包括多个模型子集, 每个模型子集独立运行IMM, 选取概率最高的模型子集的估计结果作为最终的估计状态输出。NVSIMM算法框图如图 3所示。

算法包括以下4个步骤。

步骤1  并行独立IMM状态估计

针对不同的模型子集, 分别运行IMM算法得到状态估计值$\hat{\boldsymbol{x}}_{k \mid k}(n)$、对应的协方差P_k|k(n), 以及各个模型子集中模型概率u_kⁱ(n)、似然函数Λ_kⁱ(n)和预测概率c_kⁱ(n), 这里n是模型子集的编号, 1≤n≤N, i是模型子集中模型的编号。

步骤2  计算每个模型子集的似然函数

模型子集的似然函数为

(26) $\begin{gathered}p\left(z_k \mid \boldsymbol{Z}^{k-1}, \Pi_k(n)\right)= \\\sum\limits_i p\left(z_k \mid \boldsymbol{Z}^{k-1}, M_k^i, \Pi_k(n)\right) p\left(M_k^i \mid \boldsymbol{Z}^{k-1}, \Pi_k(n)\right)= \\\sum\limits_i \mathit{\Lambda}_k^i(n) c_k^i(n)\end{gathered}$

式中: Π_k(n)表示第n个模型子集, Λ_kⁱ(n)和c_kⁱ(n)分别在步骤3和步骤1中求出。记p(z_k|Z^k－1, Π_k(n))为C_k(n), 则有:

(27) $C_k(n)=\sum\limits_i \mathit{\Lambda}_k^i(n) c_k^i(n)$

步骤3  计算模型子集概率

记模型子集概率p(Π_k(m)|Z^k)=κ_k(n), 则有:

(28) $\begin{gathered}\kappa_k(n)=\frac{p\left(z_k \mid \boldsymbol{Z}^{k-1}, \Pi_k(n)\right) p\left(\Pi_k(n) \mid \boldsymbol{Z}^{k-1}\right)}{\sum\limits_n p\left(z_k \mid \boldsymbol{Z}^{k-1}, \Pi_k(n)\right) p\left(\Pi_k(n) \mid \boldsymbol{Z}^{k-1}\right)}= \\\frac{p\left(z_k \mid \boldsymbol{Z}^{k-1}, \Pi_k(n)\right) p\left(\Pi_k(n) \mid \boldsymbol{Z}^{k-1}\right)}{p\left(z_k \mid \boldsymbol{Z}^{k-1}\right)}\end{gathered}$

式中: p(z_k|Z^k－1, Π_k(n))在步骤2中已经求出, 分母为归一化因子, p(Π_k(n)|Z^k－1)根据全概率公式有:

(29) $\begin{gathered}p\left(\Pi_k(n) \mid \boldsymbol{Z}^{k-1}\right)= \\\sum\limits_m p\left(\Pi_k(n) \mid \Pi_{k-1}(m), \boldsymbol{Z}^{k-1}\right) p\left(\Pi_{k-1}(m) \mid \boldsymbol{Z}^{k-1}\right)\end{gathered}$

式中: κ_k－1(m)=p(Π_k－1(m)|Z^k－1)为上一时刻的模型子集概率, 应用马尔可夫性质可知模型子集转移概率与观测值无关, 即:

(30) $\begin{gathered}p\left(\Pi_k(n) \mid \Pi_{k-1}(m), \boldsymbol{Z}^{k-1}\right)= \\p\left(\Pi_{k+1}(n) \mid \Pi_k(m)\right)=\sum\limits_j \prod\limits_i p_{i j}\end{gathered}$

步骤4  NVSIMM滤波状态估计

选择概率最大的模型子集的IMM估计结果作为最终的状态估计输出, 首先求出概率最高的模型子集编号:

(31) $n_m=\max\limits _n\left\{\kappa_k(n)\right\}$

从而可以得到最终的状态估计和对应的协方差为

(32) $\hat{\boldsymbol{x}}_{k \mid k}=\hat{\boldsymbol{x}}_{k \mid k}\left(n_m\right)$

(33) $\boldsymbol{P}_{k \mid k}=\boldsymbol{P}_{k \mid k}\left(n_m\right)$

本节在NVSIMM的基础上, 对数据进行进一步平滑, 得到VSIMMS算法。VSIMMS算法框图如图 4所示。

算法包括以下4个步骤。

步骤1  并行独立IMM状态平滑

针对不同的模型子集, 分别利用后向IMM平滑算法, 得到每个模型的平滑估计$\hat{\boldsymbol{x}}_{t \mid k}(n)$、对应的协方差矩阵P_t|k(n)。另外, 对模型全集Ω运行IMM滤波算法, 得到状态估计值$\hat{\boldsymbol{x}}_{t \mid k})$(Ω)和协方差估计值P_t|k(Ω)。

步骤2  计算平滑后的模型子集的似然函数

应用全概率公式和马尔可夫性质, 可以得到模型子集的似然函数为

(34) $\begin{gathered}p\left(z_k \mid \boldsymbol{Z}^{k-1}, \Pi_t(n)\right)= \\\sum\limits_m p\left(z_k \mid \boldsymbol{Z}^{k-1}, \Pi_k(m), \Pi_t(n)\right) p\left(\Pi_k(m) \mid \boldsymbol{Z}^{k-1}, \Pi_t(n)\right)= \\\sum\limits_m p\left(z_k \mid \boldsymbol{Z}^{k-1}, \Pi_k(m)\right) p\left(\Pi_k(m) \mid \boldsymbol{Z}^{k-1}, \Pi_t(n)\right)= \\\sum\limits_m p\left(z_k \mid \boldsymbol{Z}^{k-1}, \Pi_k(m)\right) p\left(\Pi_k(m) \mid \Pi_t(n)\right)\end{gathered}$

式中: p(z_k|Z^k－1, Π_k(m))为前向NVSIMM滤波过程中模型子集的似然函数, 应用全概率公式和马尔可夫性质, p(Π_k(m)|Π_t(n))可以转化为

(35) $\begin{gathered}p\left(\Pi_k(m) \mid \Pi_t(n)\right)= \\\sum\limits_{n_{t+1}} p\left(\Pi_k(m) \mid \Pi_{t+1}\left(n_{t+1}\right), \Pi_t(n)\right) p\left(\Pi_{t+1}\left(n_{t+1}\right) \mid \Pi_t(n)\right)= \\\sum\limits_{n_{t+1}} p\left(\Pi_k(m) \mid \Pi_{t+1}\left(n_{t+1}\right)\right) p\left(\Pi_{t+1}\left(n_{t+1}\right) \mid \Pi_t(n)\right)= \\\sum\limits_{n_{t+1}}\left[\sum\limits_{n_{t+2}} p\left(\Pi_k(m) \mid \Pi_{t+2}\left(n_{t+2}\right)\right) p\left(\Pi_{t+2}\left(n_{t+2}\right) \mid \Pi_{t+1}\left(n_{t+1}\right)\right)\right] \cdot \\p\left(\Pi_{t+1}\left(n_{t+1}\right) \mid \Pi_t(n)\right)=\cdots \cdots= \\\sum\limits_{n_{t+1}}\left\{\sum _ { n _ { t + 2 } } \left[\cdots \sum\limits_{n_{k-1}} p\left(\Pi_k(m) \mid \Pi_{k-1}\left(n_{k-1}\right)\right) \cdot\right.\right. \\\left.\left.p\left(\Pi_{k-1}\left(n_{k-1}\right) \mid \Pi_{k-2}\left(n_{k-2}\right)\right) \cdots\right] p\left(\Pi_{t+2}\left(n_{t+2}\right) \mid \Pi_{t+1}\left(n_{t+1}\right)\right)\right\} \cdot \\p\left(\Pi_{t+1}\left(n_{t+1}\right) \mid \Pi_t(n)\right)\end{gathered}$

式中: p(Π_t+1(n_t+1)|Π_t(n)), …, p(Π_k(m)|Π_k－1(n_k－1))为模型子集转移概率, 在第2.1节步骤3中已经求出, 记模型子集转移概率$p\left(\Pi_k(m) \mid \Pi_{k-1}\left(n_{k-1}\right)=P_{k \mid k-1}^{m \mid n_{k-1}}\right.$, 则式(35)变为

(36) $P_{k \mid t}^{m \mid n}=\sum\limits_{n_{t+1}}\left\{\sum\limits_{n_{t+2}}\left(\cdots \sum\limits_{n_{k-1}} P_{k \mid k-1}^{m \mid n_{k-1}} P_{k-1 \mid k-2}^{n_{k-1}\mid n_{k-2}} \cdots\right) P_{t+2 \mid t+1}^{n_{t+2}\mid n_{t+1}}\right\} P_{t+1 \mid t}^{\left.n_{t+1}\right|n_t}$

记$C_t^s(n)=p\left(z_k \mid Z^{k-1}, \Pi_t(n)\right)$, 则式(34)变为

(37) $C_t^s(n)=\sum\limits_m C_k(m) P_{k \mid t}^{m \mid n}$

步骤3  计算平滑后的模型子集概率

根据贝叶斯定理可以得到模型子集概率为

(38) $\begin{gathered}p\left(\Pi_t(n) \mid \boldsymbol{Z}^k\right)= \\ \frac{p\left(z_k \mid \boldsymbol{Z}^{k-1}, \Pi_t(n)\right) p\left(\Pi_t(n) \mid \boldsymbol{Z}^{k-1}\right)}{\sum\limits_n p\left(z_k \mid \boldsymbol{Z}^{k-1}, \Pi_t(n)\right) p\left(\Pi_t(n) \mid \boldsymbol{Z}^{k-1}\right)}\end{gathered}$

式中: p(Π_t(n)|Z^k－1)在上一时刻的平滑过程中已经求出, 为已知值。p(z_k|Z^k－1, Π_t(n))在步骤2中求出。记p(Π_t(n)|Z^k)=κ_t|k^s(n), 式(38)可以变换为

(39) $\kappa_{t \mid k}^s(n)=\frac{p\left(z_k \mid \boldsymbol{Z}^{k-1}, \Pi_t(n)\right) \kappa_{t \mid k-1}^s(n)}{\sum\limits_n p\left(z_k \mid \boldsymbol{Z}^{k-1}, \Pi_t(n)\right) \kappa_{t \mid k-1}^s(n)}$

步骤4  VSIMMS状态估计

概率最高的模型子集编号为

(40) $n_m^s=\max \limits_n\left\{\kappa_{t \mid k}^s(n)\right\}$

判断概率最高的模型子集在T_switch时间段内是否发生变化。若发生变化, VSIMMS的状态估计和对应的协方差为

(41) $\hat{\boldsymbol{x}}_{t \mid k}=\hat{\boldsymbol{x}}_{t \mid k}(\mathit{\Omega})$

(42) $\boldsymbol{P}_{t \mid k}=\boldsymbol{P}_{t \mid k}(\mathit{\Omega})$

若不发生变化, VSIMMS的状态估计和对应的协方差为

(43) $\hat{\boldsymbol{x}}_{t \mid k}=\hat{\boldsymbol{x}}_{t \mid k}\left(n_m^s\right)$

(44) $\boldsymbol{P}_{t \mid k}=\boldsymbol{P}_{t \mid k}\left(n_m^s\right)$

为了验证VSIMMFS算法的有效性, 本节首先针对三维空间中的目标运动轨迹展开分析讨论。目标的初始位置、速度和加速度分别为(10 km, 40 km, 30 km)、(300 m/s, 0, 0)和(0, 0, 0)。目标的机动参数如下:

(1) 30 s的匀速直线运动;

(2) 30 s的匀加速运动, 加速度大小为(-10 m/s², -10 m/s², -10 m/s²);

(3) 30 s的匀速转弯运动, 转弯角速率为(0.2 rad/s, 0.2 rad/s, 0);

(4) 30 s的匀速直线运动。

目标运动轨迹如图 5所示。

利用蒙特卡罗法分析对比IMM、NVSIMM和VSIMMFS的性能, 其中IMM算法包括3个模型: 匀速直线运动(constant velocity, CV)模型、匀加速运动(constant acceleration, CA)模型和匀速转弯(constant turn, CT)模型^[23-25]。VSIMMFS算法包括两个模型子集, 第一个模型子集包括CV、CA, 第二个模型子集包括CV、CT; NVSIMM滤波算法包括的模型子集和VSIMM相同。可以注意到, VSIMMFS和NVSIMM仅需两个原始模型子集即可包括所有的模型, 但这两种算法的原始模型子集并不包括模型子集CA、CT。设模型转移概率矩阵为

$\boldsymbol{P}=\left(\begin{array}{ccc}0.95 & 0.025 & 0.025 \\0.025 & 0.95 & 0.025 \\0.025 & 0.025 & 0.95\end{array}\right)$

假设系统过程噪声的标准差为1 m/s², 测量噪声的标准差为50 m, T_switch=20 s, 仿真结果如图 6~图 9所示。图 6给出了IMM算法中模型概率变化曲线, 其中, 概率最大的曲线所对应的模型与目标实际机动模型相同。图 7给出了NVSIMM和VSIMMFS中模型子集概率曲线, 由图 7可知, 平滑后(VSIMMFS)的模型子集概率大于NVSIMM的模型子集概率, 说明平滑后的模型子集与目标机动模型匹配程度更高。另外, 对比图 6和图 7可知, 图 6中概率最大的曲线所对应的模型属于图 7中概率高的模型子集, 也就是说, 概率最高的模型子集始终包括目标实际的机动模型。图 8给出了IMM、NVSIMM和VSIMMFS对目标运动轨迹的跟踪效果, 这3种方法都能够对目标轨迹进行较好的估计, 三者之间的差异具体体现在图 9中。图 9给出了IMM、NVSIMM和VSIMMFS对目标位置估计的均方根误差(root mean square error, RMSE)曲线。从图 9中可以看出, VSIMMFS的RMSE最小, VSIMMFS次之, IMM的RMSE最大。值得注意的是, 概率最大的模型子集发生变化时, NVSIMM算法的性能会短暂下降, 而VSIMMFS算法对位置估计的RMSE与IMM相同, 原因在于在设计VSIMMFS算法时, 当最大概率的模型子集短暂发生时, 利用基于模型全集的IMM估计值作为VSIMMFS的估计值, 从而避免了VSIMMFS算法短暂的性能下降。

为了进一步体现VSIMMFS算法的优越性, 将VSIMMFS与基于深度强化学习的目标跟踪算法(maneuvering target tracking based on deep reinforcement learning, MTTDRL)^[27]、粒子滤波(particle filter, PF)^[28]、基于粒子群优化的PF(particle swarm optimization algorithm-based PF, PSO-PF)^[29]、基于混沌PSO-PF(chaos PSO-PF, CPSO-PF)^{[26, 30]}进行对比。为了保证对比的有效性, 仿真参数参考文献[26]和文献[27]设计, 目标的运动轨迹如图 10所示。图 11分析对比了上述算法对目标位置的跟踪性能。由图 11(a)可知, 相比于PF、PSO-PF和CPSO-PF, VSIMMFS的位置跟踪的RMSE较小; 从图 11(b)可以看出, 利用VSIMMFS对目标跟踪时的位置偏差更小。所以综上可知, VSIMMFS的跟踪性能优于上述其他算法。

为了进一步提高对机动目标的跟踪精度, 本文设计了一种VSIMMFS算法。该方法包括前向VSIMM滤波和后向VSIMMS两个部分, 通过并行独立运行IMM滤波和IMM平滑, 选择概率最大的模型子集的估计结果作为最终的状态估计值。通过与IMM、NVSIMM、PF、PSO-PF、CPSO-PF和MTTDRL算法比较, VSIMMFS具有更好的跟踪性能。