2. 上海大学 人工智能研究院,上海 200444
2. School of Artificial Intelligence, Shanghai University, Shanghai 200444, China
无人水面艇的目标跟踪可被视为一个随机过程,而目标跟踪问题的建模是通过数学表达式描述机动目标当前运动模式,并准确估计机动目标位置、速度等运动状态[1]。常见模型有匀速运动模型(Constant Velocity,CV) 、匀速转弯运动模型(Constant Turn,CT)、匀加速度运动模型(Constant Acceleration,CA) 、Singer模型[2]以及当前统计(Current Statistical,CS)模型。但单一的运动模型在处理复杂的机动目标跟踪问题上,容易出现较大误差,甚至出现目标丢失情况[3]。针对单一模型的不足,Bar-Shalom 提出了交互式多模型 (Interactive Multiple Model,IMM)算法[4]。该算法将目标运动映射到多个已知的模型集中,对每个模型并行滤波。其原理是利用 Markov 概率转移矩阵进行不同子模型之间的交互和切换,并且使得这些模型都相互作用,从各个模型的滤波器得出的状态估计,最后经过加权得到最终状态估计[5]。
IMM通过组合不同模型进行跟踪,在一定程度上能适应目标多种机动情况,相对于单模型来说,提高了对目标跟踪的自适应性。传统的IMM算法采用的是固定模型集和固定参数,这使得处理复杂场景的机动目标跟踪问题时,只能通过增加模型个数来提高测量的精度。一方面这会降低整体算法的运算速度,在另一方面模型之间的竞争也会降低跟踪的精度问题[6]。在实际跟踪过程中,以船舶为主的水面障碍物由于强随机性和高动态性,其运动模型会发生频繁切换,直接将具有固定模型集和固定参数的IMM算法应用于无人水面艇目标跟踪场景时,实验表明会造成算法稳定性和跟踪精度不高的问题。
有学者就此提出了一些改进方案。例如变结构多模型(Variable Structure Multiple Model,VSMM)算法[7],异构模型集协同滤波(Novel Interactive Multiple Model,Novel-IMM)算法[8]以及FAIMM(Fuzzy-logic Adaptive Interactive Multiple Model)算法[9]。尽管VSMM 算法中利用时变模型集来替代IMM中的固定模型集,但是频繁的切换、激活以及初始化模型集会导致滤波估计稳定性下降,影响跟踪效果。Novel-IMM 算法在面对机动目标运动模式时,不需要频繁地切换场景,具有较好的跟踪效果,但是计算量却是成倍增加[10]。而FAIMM算法经过实验表明,对机动目标运动模式的变化比较敏感,容易输出错误模型,影响跟踪效果。因此,针对无人水面艇目标跟踪场景,除了需要计算结果的实时性,还对跟踪准确性提出了较高要求,以上算法均难以同时满足。
1 交互式多模型算法IMM算法引入了马尔可夫状态转移矩阵[11],内部包含多个模型和与之对应的滤波器,算法将每个模型都参与到状态估计的运算中,各模型之间按照马尔可夫转移概率进行切换[12]。假设 IMM 采用 N 个运动模型,则模型集
将上一时刻各运动模型在滤波器中的输出结果进行交互,作为当前迭代过程中的初始状态。
$X_{k-1\mid k-1}^{oj}=\sum _{i=1}^{N}{X}_{k-1\mid k-1}^{i}\cdot {\mu }_{k-1}^{ij},i,j=1,2,\dots ,n,$ | (1) |
$ \begin{split} {P}_{k-1\mid k-1}^{oj}=& \sum _{i=1}^{N}{\mu }_{k-1}^{ij}\cdot [{P}_{k-1\mid k-1}^{i}+ ({X}_{k-1\mid k-1}^{i}-{X}_{k-1\mid k-1}^{oj})\times\\ & \left.{\left({X}_{k-1\mid k-1}^{i}-{X}_{k-1\mid k-1}^{oj}\right)}^{\mathrm{T}}\right]。\end{split} $ | (2) |
其中:
$ {\mu }_{k-1}^{ij}=\frac{1}{\overline{{C}_{j}}}{{\text{π}} }^{ij}{\mu }_{k-1}^{i},$ | (3) |
$ \overline{{C}_{j}}=\sum _{i=1}^{N}{{\text{π}} }^{ij}\cdot {\mu }_{k-1}^{i}。$ | (4) |
式中:
容积卡尔曼滤波器(Cubature Kalman Filter,CKF)[13],相比其他非线性滤波器如EKF、UKF等,其数值精度更高,稳定性更好。本文采用文献[14]提出的SRCKF滤波器进行目标状态估计。其过程如下:
步骤1 初始化。给定初始条件
${\widehat{\boldsymbol{x}}}_{0}=E\left({\boldsymbol{x}}_{0}\right),$ | (5) |
$ {\boldsymbol{P}}_{0}=\mathrm{cov}\left({\boldsymbol{x}}_{0}\right)。$ | (6) |
步骤2 计算容积点
$ \left\{\begin{split}&{\xi }_{j}=\sqrt{\dfrac{m}{2}}[l{]}_{j},\\& {\omega }_{j}=\dfrac{1}{m},\end{split}\right.j=1,2,\cdots m。$ | (7) |
式中:
步骤3 时间更新
1)对k-1时刻的协方差矩阵进行Cholesky 分解
$ {\boldsymbol{P}}_{k-1}={{{\boldsymbol{S}}}}_{k-1}{{{\boldsymbol{S}}}}_{k-1}^{{\rm{T}}}。$ | (8) |
2)计算状态预测容积点和传播容积点
$ {\boldsymbol{\zeta }}_{j,k-1}={\boldsymbol{S}}_{k-1}{\boldsymbol{\xi }}_{j}+{\widehat{\boldsymbol{x}}}_{k-1},$ | (9) |
$ {\boldsymbol{X}}_{j,\frac{k}{k}-1}^{\mathrm{*}}=f\left({\zeta }_{j,k-1}\right)。$ | (10) |
3)状态预测估计和预测估计协方差
$ \left\{\begin{aligned} & {\widehat{\boldsymbol{x}}}_{k|k-1}=\sum _{j=1}^{m}{\omega }_{j}{\boldsymbol{X}}_{j,\frac{k}{k}-1}^{\mathrm{*}},\\ & {\boldsymbol{S}}_{k\mid k-1}=\mathrm{Tria}\left(\left[{\chi }_{k\mid k-1}^{\mathrm{*}},{\boldsymbol{S}}_{Q,k-1}\right]\right)。\end{aligned}\right. $ | (11) |
式中:
步骤4 量测更新
1)计算量测容积点和传播容积点
$ \left\{\begin{aligned} & {\zeta }_{j,k|k-1}={\boldsymbol{S}}_{k|k-1}{\boldsymbol{\xi }}_{j}+{\widehat{\mathbf{x}}}_{k|k-1},\\ & {\boldsymbol{y}}_{j,k|k-1}^{\mathrm{*}}=h\left({\zeta }_{j,k|k-1}\right)。\end{aligned}\right. $ | (12) |
2)量测估计和量测估计协方差
$ \left\{\begin{aligned} & {\widehat{\boldsymbol{y}}}_{k|k-1}=\sum _{j=1}^{m}{\omega }_{j}{\boldsymbol{y}}_{\mathrm{j},\mathrm{k}|\mathrm{k}-1}^{\mathrm{*}},\\ & {\boldsymbol{S}}_{yy}=Tria\left(\right[{\eta }_{k\mid k-1},\sqrt{{\boldsymbol{R}}_{k-1}}\left]\right)。\end{aligned}\right. $ | (13) |
$ \begin{split} {\eta }_{k\mid k-1}=& \frac{1}{\sqrt{2n}}\left[\left({\boldsymbol{y}}_{1,{k}|{k}-1}^{\mathrm{*}}-{\widehat{\boldsymbol{y}}}_{k|k-1}\right),\left({\boldsymbol{y}}_{2,{k}|{k}-1}^{\mathrm{*}}-{\widehat{\boldsymbol{y}}}_{k|k-1}\right),\cdots ,\right.\\ & \left.\left({\boldsymbol{y}}_{2n,k|k-1}^{\mathrm{*}}-{\widehat{\boldsymbol{y}}}_{k|k-1}\right)\right]。\\[-10pt] \end{split} $ | (14) |
3)获得互协方差矩阵
$ {\boldsymbol{P}}_{xy}={\boldsymbol{\gamma }}_{k\mid k-1}{\left({\boldsymbol{\eta }}_{k\mid k-1}\right)}^{{\rm{T}}} ,$ | (15) |
$ \begin{split} {\gamma }_{k\mid k-1}=& \frac{1}{\sqrt{2n}}\left[\left({\zeta }_{j,k|k-1}-{\widehat{\mathbf{x}}}_{k|k-1}\right),\left({\zeta }_{j,k|k-1}-{\widehat{\mathbf{x}}}_{k|k-1}\right),\cdots ,\right.\\ & \left.\left({\zeta }_{j,k|k-1}-{\widehat{\mathbf{x}}}_{k|k-1}\right)\right]。\\[-10pt] \end{split} $ | (16) |
4)增益更新
$ {\boldsymbol{K}}_{k}=({\boldsymbol{P}}_{xy}∕{\boldsymbol{S}}_{yy}^{{\rm{T}}})∕{\boldsymbol{S}}_{yy} 。$ | (17) |
5)计算状态估计值和平方根形式的估计协方差
$ \left\{\begin{aligned} & {\widehat{\boldsymbol{x}}}_{k}={\widehat{\boldsymbol{x}}}_{k|k-1}+{\boldsymbol{K}}_{k}\left({\boldsymbol{y}}_{k}-{\widehat{\boldsymbol{y}}}_{k|k-1}\right),\\ & {\boldsymbol{S}}_{k}=\mathrm{Tria}\left(\left[{\gamma }_{k\mid k-1}-{\boldsymbol{K}}_{k}{\eta }_{k\mid k-1},{\boldsymbol{K}}_{k}\sqrt{{\boldsymbol{R}}_{k-1}}\right]\right)。\end{aligned}\right. $ | (18) |
通常假设模型的似然函数服从正态分布。因此k时刻模型j匹配的似然函数以及更新后的模型概率分别为
$ \left\{\begin{aligned} & {\mathrm{\Lambda }}_{k}^{j}=\frac{1}{\sqrt{\left|2{\text{π}} {\mathbf{S}}_{k}^{j}\right|}}\mathrm{exp}\left[-\frac{1}{2}{\left({\boldsymbol{e}}_{k}^{j}\right)}^{{\rm{T}}}{\left({\boldsymbol{S}}_{k}^{j}\right)}^{-1}{\boldsymbol{e}}_{k}^{j}\right],\\ & {u}_{k}^{j}=\frac{{\mathrm{\Lambda }}_{k}^{j}\overline{{C}_{j}}}{\displaystyle\sum _{i=1}^{n}{\mathrm{\Lambda }}_{k}^{i}\overline{{C}_{i}}}。\end{aligned}\right. $ | (19) |
其中,滤波残差
$ \left\{\begin{aligned} & {e}_{k}^{j}={\boldsymbol{y}}_{k}-{\widehat{\boldsymbol{y}}}_{k|k-1},\\ & {\boldsymbol{S}}_{k}^{j}={\boldsymbol{S}}_{zz,k|k-1}{\boldsymbol{S}}_{zz,k|k-1}^{{\rm{T}}}。\end{aligned}\right.$ | (20) |
将各个滤波器输出的估计值根据模型概率进行加权融合,得到IMM算法的最终状态估计输出
$ \left\{\begin{aligned} & {\widehat{\boldsymbol{x}}}_{k}=\sum _{j=1}^{r}{u}_{k}^{j}{\widehat{\boldsymbol{x}}}_{k}^{j},\\ & {\boldsymbol{P}}_{k}=\sum _{j=1}^{r}{u}_{k}^{j}\left[{\boldsymbol{P}}_{k}^{j}+\left({\widehat{\boldsymbol{x}}}_{k}^{j}-{\boldsymbol{x}}_{k}\right){\left({\widehat{\boldsymbol{x}}}_{k}^{j}-{\boldsymbol{x}}_{k}\right)}^{{\rm{T}}}\right]。\end{aligned}\right.$ | (21) |
假设N=3,则IMM算法结构图如图1所示。
IMM算法综合考虑了各模型之间的切换,相比传统单模型跟踪,提供了更高的跟踪精度和更快的收敛速度。当机动目标包含多种运动模式时,IMM模型需包含尽可能多的运动模型才能更好匹配目标运动状态,而过多的运动模型也意味着非匹配模型的竞争,引入误差,导致跟踪精度下降的同时,增大算法运算量。
2 SWA-IMM为解决IMM算法固定模型集带来的局限性问题,同时满足无人艇在水面机动目标跟踪场景下对跟踪精度的要求,水面机动目标常见且高频的3种运动模型分别是CV、CT、CA模型,令模型数量n=3,
根据模型权重,初始化每个独立滑动窗口模型集的模型概率转移矩阵
滑窗判定方法核心是m/n滑窗逻辑法[15],序列
SWA-IMM算法流程图如图3所示。
由于滑动窗口中,模型集前几个时刻先验信息及后验信息较少,而SWA-IMM算法的实现依赖于窗口内长度为n的序列。因此前n-1时刻
$ \begin{split} & {\overline{\boldsymbol{X}}}_{k}=\frac{1}{3}{\left[\begin{array}{ccc}{\lambda }_{cv}& {\lambda }_{ct}& {\lambda }_{ca}\end{array}\right]}^{{\rm{T}}}\left[\begin{array}{ccc}{\widehat{\boldsymbol{x}}}_{k}\left(1\right)& {\widehat{\boldsymbol{x}}}_{k}\left(2\right)& {\widehat{\boldsymbol{x}}}_{k}\left(3\right)\end{array}\right],\\ & {\overline{\boldsymbol{P}}}_{k}=\frac{1}{3}{\left[\begin{array}{ccc}{\lambda }_{cv}& {\lambda }_{ct}& {\lambda }_{ca}\end{array}\right]}^{{\rm{T}}}\left[\begin{array}{ccc}{\boldsymbol{P}}_{k}\left(1\right)& {\boldsymbol{P}}_{k}\left(2\right)& {\boldsymbol{P}}_{k}\left(3\right)\end{array}\right]。\end{split} $ | (22) |
式中,
当窗口内序列长度大于n时,可进行滑窗判定。输入目标状态以及状态协方差矩,并执行IMM算法得到状态估计
滑窗判定主要是将窗口中,概率值最大值对应的模型
1)获取最佳匹配模型
检测数m和序列n一起构成滑动窗口的逻辑判定。对于序列中的n个时刻,获取每个时刻
$ F\left(w,m\right)=\left\{\begin{array}{cc}\varnothing,& count\left({\boldsymbol{M}}_{\boldsymbol{w}}^{\boldsymbol{*}}\right) < m,\\ {\boldsymbol{M}}^{\boldsymbol{*}},& count\left({\boldsymbol{M}}_{\boldsymbol{w}}^{\boldsymbol{*}}\right)\geqslant m。\end{array}\right. $ | (23) |
2)设定滑窗判定系数
$ {\lambda }_{cv},{\lambda }_{ct},{\lambda }_{ca}=\left\{\begin{array}{c}1,F\left(w,m\right)={\boldsymbol{M}}^{\boldsymbol{*}},\\ 0 ,F\left(w,m\right)\ne {\boldsymbol{M}}^{\boldsymbol{*}}。\end{array}\right. $ | (24) |
式中:
当多个滑动窗口均判定其权重模型有效时,取所有窗口结果的算术平均值作为结果。
$ \begin{split} {\overline{\boldsymbol{X}}}_{k}=& \dfrac{1}{{\lambda }_{cv}+{\lambda }_{ct}+{\lambda }_{ca}}{\left[\begin{array}{ccc}{\lambda }_{cv}& {\lambda }_{ct}& {\lambda }_{ca}\end{array}\right]}^{{\rm{T}}}\times\\ & \left[\begin{array}{ccc}{\widehat{\boldsymbol{x}}}_{k}\left(1\right)& {\widehat{\boldsymbol{x}}}_{k}\left(2\right)& {\widehat{\boldsymbol{x}}}_{k}\left(3\right)\end{array}\right],\\ {\overline{\boldsymbol{P}}}_{k}=& \dfrac{1}{{\lambda }_{cv}+{\lambda }_{ct}+{\lambda }_{ca}}{\left[\begin{array}{ccc}{\lambda }_{cv}& {\lambda }_{ct}& {\lambda }_{ca}\end{array}\right]}^{{\rm{T}}}\times\\ & \left[\begin{array}{ccc}{\boldsymbol{P}}_{k}\left(1\right)& {\boldsymbol{P}}_{k}\left(2\right)& {\boldsymbol{P}}_{k}\left(3\right)\end{array}\right]。\end{split} $ | (25) |
假设目标在x和y方向的初始位置分别为[100 m, 100 m],初始加速度为[0.8 m/s2, 0.8 m/s2]。目标的轨迹由多段不同运动状态下的路径组成:该目标在0~20 s时做加速度[0.8 m/s2, 0.8 m/s2]的匀加速度直线运动;20~150 s 时做角速度为0.05 rad/s的匀速转弯运动;150~180 s时做加速度为[1 m/s2, 1 m/s2]的匀加速度运动。产生的机动目标运动仿真轨迹如图4所示。
仿真中目标运动存在匀速直线运动、匀速转弯运动、匀加速直线运动3种状态。假设采样间隔T=1 s,设定量测向量
$ P=\left[\begin{array}{ccc}0.90& 0.25& 0.25\\ 0.25& 0.90& 0.25\\ 0.25& 0.25& 0.90\end{array}\right]。$ | (26) |
对于SAC-IMM算法,设定滑动窗口中模型概率初值为:
$ \varPi =\left[\begin{array}{ccc}0.6& 0.2& 0.2\\ 0.2& 0.6& 0.2\\ 0.2& 0.2& 0.6\end{array}\right]。$ | (27) |
设定滑动窗口中概率转移矩阵
$ {PP}_{i}=\left[\begin{array}{ccc}0.9& 0.05& 0.05\\ 0.05& 0.9& 0.05\\ 0.05& 0.05& 0.9\end{array}\right],i=1,2,3。$ | (28) |
图5为笛卡尔坐标系下 IMM、SWA-IMM 算法的跟踪轨迹。可知,SWA-IMM算法相比传统IMM算法在机动目标发生运动状态改变时,能够更加有效地结合量测结果,对目标进行轨迹跟踪。
图6所示为目标真实运动轨迹、量测轨迹以及2种滤波算法轨迹的局部对比图。可以发现,无论目标处于匀速转弯运动阶段还是匀加速运动阶段,本文提出的SWA-IMM滤波跟踪算法相比传统IMM算法,在滤波轨迹上更接近于真实轨迹。
对于2种算法性能,使用目标位置、速度、加速度的均方根误差作为评价指标。
$ {\boldsymbol{RMSE}}=\sqrt{\frac{1}{M}\sum _{j=1}^{M}{\left({\boldsymbol{X}}_{i}\left(k\right)-{\overline{\boldsymbol{X}}}_{i}\left(k\right)\right)}^{2}} 。$ | (29) |
式中:
根据表1中的数据分析可知,SWA-IMM算法相较于固定模型集IMM算法,在位置、速度和加速度性能上分别提高了33%、39%和17%。
分析图7~图9可知,目标处于弱机动或强机动状态时,本文提出的算法相比较传统IMM算法无论在位置还是速度亦或是加速度上的跟踪误差都有所提升,在匹配目标运动模型时效果良好。这是因为SWA-IMM算法通过记录历史信息,当滑动窗口判定目标运动状态发生改变时,及时调窗口模型集中失配模型和配准模型的权重比,即提高配准模型的概率,从而比固定模型集得到更加准确的匹配结果。由图10可看出,在目标运动的0~30 s和150~180 s内,权重模型的概率分布基本接近目标实际运动状态;在目标运动的130~150 s内,即目标处于匀加速转弯状态时,由于滤波结果不稳定,多个窗口满足滑窗判定,最终的模型集为多个窗口模型集融合结果,因此最佳模型的匹配度相比于其他运动阶段会有所下降。
本文针对无人水面艇目标跟踪场景,对固定模型集的交互式多模型算法进行改进,提出了SWA-IMM算法。实验表明,相比传统IMM算法,SWA-IMM 算法能更加准确获取运动目标速度、加速度等信息,使得无人艇在跟踪水面机动目标时能更加快速的匹配当前目标的运动状态,并完成目标的跟踪定位以及运动信息获取。
[1] |
沙俊臣. 基于模型集自适应多模交互滤波的机动目标跟踪算法研究及性能分析[D]. 济南: 山东大学, 2013.
|
[2] |
ROBERT A Singer. Estimating optimal tracking filter performance for manned maneuvering targets[J]IEEE Transactions on Aerospace and Electronic Systems, 1970, 4(6): 473–483.
|
[3] |
陈亮. 机动目标跟踪关键技术研究[D]. 哈尔滨: 哈尔滨工程大学, 2012.
|
[4] |
MAZOR E, AVERBUCH A, BAR-SHALOM Y. Interacting multiple model methods in target tracking: a survey[J]. IEEE Transactions on Aerospace and Electronic Systems, 1998, 1(34): 103-123. |
[5] |
盛涛, 夏海宝, 肖冰松. 基于AIMM-SRCKF的机动目标跟踪算法[J]. 电子测量与仪器学报, 2021, 35(1): 159-164. DOI:10.13382/j.jemi.B2003199 |
[6] |
潘媚媚, 曹运合, 王宇, 等. 基于机动判别的变结构交互多模型跟踪算法[J]. 系统工程与电子技术, 2019, 41(4): 730-736. |
[7] |
毛少锋, 冯新喜, 文莎, 等. 具有固定滞后平滑的变结构机动目标跟踪[J]. 重庆邮电大学学报(自然科学版), 2015, 27(2): 189-196. |
[8] |
QU Hongquan, PANG Liping, LI Shaohong. A novel interacting multiple model algorithm[J]. Signal Processing, 2009, 11(89): 2171-2177. |
[9] |
赵楚楚, 王子微, 丁冠华, 等. 基于模糊逻辑的改进自适应IMM跟踪算法[J]. 信号处理, 2021, 37(5): 724-734. DOI:10.16798/j.issn.1003-0530.2021.05.005 |
[10] |
王耀林, 盖梦欧, 周敏. IMM算法在空中复杂机动目标跟踪中的应用[J]. 现代信息科技, 2020, 4(13): 9-11. DOI:10.19850/j.cnki.2096-4706.2020.13.003 |
[11] |
戴定成, 姚敏立, 蔡宗平, 等. 改进的马尔可夫参数自适应IMM算法[J]. 电子学报, 2017, 45(5): 1198-1205. DOI:10.3969/j.issn.0372-2112.2017.05.024 |
[12] |
朱洪峰, 熊伟, 崔亚奇, 等. 基于加速度的马尔可夫参数自适应IMM算法[J]. 火力与指挥控制, 2019, 44(11): 46-50+57. |
[13] |
ARASARATNAM I, HAYKIN S. Cubature kalman filters[J]. IEEE Transactions on Automatic Control, 2009, 6(54): 1254-1269. |
[14] |
郝燕玲, 杨峻巍, 陈亮, 等. 平方根容积卡尔曼滤波器[J]. 弹箭与制导学报, 2012, 32(2): 169-172. DOI:10.3969/j.issn.1673-9728.2012.02.047 |
[15] |
董志荣. 论航迹起始方法[J]. 情报指挥控制系统与仿真技术, 1999(2): 1-7. |