Maneuvering Target Tracking Algorithm Based on the Adaptive Augmented State Interracting Multiple Model
-
摘要: 现有的增广状态-交互式多模型算法存在着依赖于量测噪声协方差矩阵这一先验信息的问题。当先验信息未知或不准确时,算法的跟踪性能将会下降。针对上述问题,该文提出一种自适应的变分贝叶斯增广状态-交互式多模型算法VB-AS-IMM。首先,针对增广状态的跳变马尔科夫系统,该文给出了联合估计增广状态和量测噪声协方差矩阵的变分贝叶斯推断概率模型。其次,通过理论推导证明了该概率模型是非共轭的。最后,通过引入一种“信息反馈+后处理”方案,提出联合后验密度的次优求解方法。所提算法能够在线估计未知的量测噪声协方差矩阵,具有更强的鲁棒性和适应性。仿真结果验证了算法的有效性。Abstract: The existing Augmented State-Interracting Multiple Model (AS-IMM) algorithm suffers from the problem that it relies on the prior information of the covariance matrix of the measurement noise. When the prior information is unavailable or inaccurate, the tracking performance of AS-IMM will be degraded. In order to overcome this problem, a novel adaptive Bayesian Variational Augmented State-Interracting Multiple Model (VB-AS-IMM) algorithm is proposed. Firstly, the variational Bayesian inference probabilistic model of the augmented state and the covariance matrix of the measurement noise for the jump Markovarian system is presented. Secondly, the probabilistic model is proven to be non-conjugated. Finally, by introducing a novel post processing method, the suboptimal solution to calculate the joint posterior distribution is proposed. The proposed algorithm can estimate the unknown covariance matrix of the measurement noise online, thus it is more robust and has higher adaptability. Simulation result verifies good performance of the proposed algorithm.
-
1. 引言
机动目标跟踪一直是目标跟踪领域备受关注的热点问题,其难点在于目标运动的不确定性。针对这一难题,交互式多模型(Interracting Multiple Model, IMM)算法[1]已被证明是目前最有效的手段。通过引入模型交互步骤,IMM具有1阶广义伪贝叶斯(Generalized Pseudo Bayesian, GPB)算法[2]的计算复杂度优势,同时兼备2阶GPB算法的跟踪性能,实现了跟踪精度与算法复杂度的折中。IMM已被广泛应用于各类场景下的机动目标跟踪问题。
固定步长平滑(Fixed Lag Smoothing, FLS)是一种在线的状态估计器[3],其在量测数据与待估状态之间设置一个固定的延迟步长L,利用到当前时刻k为止接收到的量测数据来估计k-L时刻的状态。由于其性能优于滤波估计器,且可以在线工作,因而被广泛应用于信号处理、目标跟踪、导航和计算机视觉等领域。Rauch[4]给出了非跳变单模型的FLS的实现形式。然而,文献[5]指出Rauch提出的FLS的实现形式存在着数值不稳定的问题。文献[6]通过引入状态增广的方式,利用传统的卡尔曼滤波(Kalman Filter, KF)来实现数值稳定的FLS,并得到了广泛的应用。文献[7]给出了跳变线性马尔科夫系统的FLS的实现形式,其采用了假设管理与剪枝的方法,算法实现较为复杂。文献[8]进一步假设增广状态后验密度可用与模型数目相同的高斯混合分布给出,进而提出了增广状态-交互式多模型(Augmented State-IMM, AS-IMM)算法。AS-IMM实现形式简单,且能实现目标状态和模型概率的平滑。文献[9]将AS-IMM扩展到了非线性系统。除了在线的FLS外,跳变马尔科夫系统的离线固定区间平滑(Fixed Interval Smoothing, FIS)也得到广泛研究,可参考文献[10]。需要注意的是,传统的AS-IMM算法假设量测协方差矩阵先验已知。而当量测协方差矩阵先验未知或先验信息不准确时,AS-IMM算法的跟踪性能将下降。
针对量测噪声分布未知条件下的目标跟踪问题,近年来基于变分贝叶斯[11,12](Variational Bayesian, VB)的自适应滤波方法备受关注。文献[13]提出了联合估计目标状态和量测噪声方差的变分贝叶斯-自适应卡尔曼滤波(Variational Bayesian Adaptive Kalman Filter, VB-AKF)。文献[14]改进了VB-AKF,进一步考虑了非对角阵的量测协方差矩阵的VB推断方法。此外,VB方法已扩展至跳变马尔科夫系统[15],非线性自适应滤波[16]和非高斯自适应滤波[17]等问题。然而,对于跳变马尔科夫系统的自适应平滑问题,目前尚未有报导,这也是本文的出发点。
本文针对传统的AS-IMM算法存在着依赖于量测协方差矩阵先验信息的问题,提出了一种自适应变分贝叶斯增广状态-交互式多模型VB-AS-IMM算法。该算法通过引入一种“信息反馈+后处理”的近似策略,利用VB求解量测噪声协方差矩阵的统计分布。所提算法能够联合估计机动目标运动状态和量测噪声分布,具有更强的鲁棒性和适应性,更适用于机动目标的跟踪问题。
2. 问题描述
对于机动目标的跟踪,本文考虑包含M个模型的跳变线性马尔科夫系统,其中第j个模型的状态空间模型为
$$ \qquad\quad {{{x}}_k} = {{F}}_{k - 1}^j{{{x}}_{k - 1}} + {{v}}_{k - 1}^j $$ (1) $$ \qquad\quad {{{z}}_k} = {{H}}_k^j{{{x}}_k} + {{{n}}_k} $$ (2) 其中,
${{{x}}_k}$ 和${{{z}}_k}$ 分别表示状态向量和量测向量。${{v}}_{k - 1}^j \sim {\cal{N}}\left( {0,{{Q}}_{k - 1}^j} \right)$ 和${{{n}}_k} \sim {\cal{N}}\left( {0,{{{R}}_k}} \right)$ 分别表示服从零均值高斯分布的过程噪声和量测噪声,且两者相互独立。模型的先验概率及马尔科夫转移概率分别为$$ P\left( {m_0^j} \right) = \mu _0^j;\;\;P\left( {m_k^j\left| {m_{k - 1}^i} \right.} \right) = {\pi _{ij}} $$ (3) 需要注意的是式(1)和式(2)中的模型包含一个假设:跳变线性马尔科夫系统的量测噪声是一致的,即
${{R}}_k^i = {{R}}_k^j = {{{R}}_k}$ ,$i,j \in \left[ {1,M} \right]$ 。这样假设的缘由如下:(1)机动目标运动的复杂性主要是由于状态方程具有多态性;(2)量测噪声主要由目标信噪比、传感器的探测精度等因素决定,可认为与目标运动状态无关。FLS可以概括为:给出平滑步长
$L$ ,利用量测数据${{{z}}_{1:k}} = \left\{ {{{{z}}_i}} \right\}_{i = 1}^k$ 来估计$k - L$ 时刻的后验密度函数$p\left( {{{{x}}_{k - L}}\left| {{{{z}}_{1:k}}} \right.} \right)$ 。当${{{R}}_k}$ 未知时,需要同时估计${{{R}}_k}$ 和目标状态${{{x}}_{k - L}}$ ,因而需要求解联合后验密度函数$p\left( {{{{x}}_{k - L}},{{{R}}_k}\left| {{{{z}}_{1:k}}} \right.} \right)$ 。3. 传统的增广状态-交互式多模型算法概述
定义增广的状态向量为
${{{X}}_k} = \left[ {{{x}}_k^{\rm{T}}}\ \ {{{x}}_{k - 1}^{\rm{T}}}\ \ ··· \ {{{x}}_{k - L}^{\rm{T}}} \right]^{\rm{T}}$ ,则跳变线性马尔科夫系统的增广状态空间模型可表示为$$ {{{X}}_k} = {{F}}_{k - 1}^{a\left( j \right)}{{{X}}_{k - 1}} + {{\varGamma}} {{{v}}_{k - 1}} $$ (4) $$ {{{z}}_k} = {{H}}_k^{a\left( j \right)}{{{X}}_k} + {{{n}}_k} = {{H}}_k^j\left[ {{{{I}}_{{n_x}}}}\ 0\ ··· \ 0 \right]{{{X}}_k} + {{{n}}_k}\!\!\!\!\! $$ (5) 其中,
$j \in \left[ {\begin{array}{*{20}{c}} 1&M \end{array}} \right]$ ;${n_x}$ 表示状态向量的维数;增广状态转移矩阵${{F}}_{k - 1}^{a\left( j \right)}$ 和状态控制矩阵${{\varGamma}} $ 参考文献[8]定义。文献[8]将$p\left( {{{{X}}_k}\left| {m_k^j,{{{z}}_{1:k}}} \right.} \right)$ 近似为由$M$ 个分量组成的高斯混合分布,提出了AS-IMM算法。为了后文推导方便,本文将AS-IMM的实现步骤概述如下:(1) 模型重初始化:与传统的IMM步骤一致,可给出增广状态的模型初始值
$$ \begin{split} {\hat{\tilde{ X}}}_{k - 1\left| {k - 1} \right.}^j \,& \triangleq {\rm{E}}\left[ {p\left( {{{{X}}_{k - 1}}\left| {m_k^j,{{{z}}_{1:k - 1}}} \right.} \right)} \right] \\ & = \sum\limits_{i = 1}^M {\mu _{k - 1\left| {k - 1} \right.}^{i\left| j \right.}{{\hat X}}_{k - 1\left| {k - 1} \right.}^i} \end{split} $$ (6) $$ \begin{split} {\tilde{ P}}_{k - 1\left| {k - 1} \right.}^j =\,& \sum\limits_{i = 1}^M \mu _{k - 1\left| {k - 1} \right.}^{i\left| j \right.}[ {{P}}_{k - 1\left| {k - 1} \right.}^i \\ & + {{\left( {{{\hat X}}_{k - 1\left| {k - 1} \right.}^i - {\hat{\tilde{ X}}}_{k - 1\left| {k - 1} \right.}^j} \right)}^{\rm{T}}}\\ & \left.\left( {{{\hat X}}_{k - 1\left| {k - 1} \right.}^i - {\hat{\tilde{ X}}}_{k - 1\left| {k - 1} \right.}^j} \right) \right] \end{split} $$ (7) 其中,
${{\hat X}}_{k - 1\left| {k - 1} \right.}^i$ 和${{P}}_{k - 1\left| {k - 1} \right.}^i$ 分别表示第$k - 1$ 时刻,第$i$ 个模型的状态估计值及其协方差矩阵。$\mu _{k - 1\left| {k - 1} \right.}^{i\left| j \right.} = P\left( {m_{k - 1}^i\left| {m_k^j,{{{z}}_{1:k - 1}}} \right.} \right)$ 表示混合概率,满足$\mu _{k - 1\left| {k - 1} \right.}^{i\left| j \right.} = \dfrac{{{\pi _{ij}}\mu _{k - 1}^i}}{{{c_j}}},\;\;{c_j} = \displaystyle\sum\nolimits_{i = 1}^M {{\pi _{ij}}\mu _{k - 1}^i}$ 。需要说明的是,对于${\tilde{ P}}_{k - 1\left| {k - 1} \right.}^j$ 的初始化,只需要计算${\tilde{ P}}_{k - 1\left| {k - 1} \right.}^{j\left( {i,i} \right)}$ 和${\tilde{ P}}_{k - 1\left| {k - 1} \right.}^{j\left( {0,i} \right)}$ ,其中${\tilde{ P}}_{k - 1\left| {k - 1} \right.}^{j\left( {i,i} \right)}$ 和${\tilde{ P}}_{k - 1\left| {k - 1} \right.}^{j\left( {0,i} \right)}$ 分别表示${\tilde{ P}}_{k - 1\left| {k - 1} \right.}^j$ 的分块对角线元素和分块第1行元素,具体定义参考文献[8]。(2) 模型条件滤波:在给定的重初始化的增广状态和协方差矩阵的条件下,模型条件滤波包含3个步骤。
(a) 状态预测:
$$ \begin{split} & {{\hat X}}_{k\left| {k - 1} \right.}^{j\left( 0 \right)} = {{F}}_{k - 1}^j{\hat{\tilde{ X}}}_{k - 1\left| {k - 1} \right.}^{j\left( 0 \right)},\;\;{{\hat X}}_{k\left| {k - 1} \right.}^{j\left( i \right)} = {\hat{\tilde{ X}}}_{k - 1\left| {k - 1} \right.}^{j\left( {i - 1} \right)},\\ & i = 1,2, ··· ,L\\[-10pt] \end{split} $$ (8) $$ \begin{split} & {{P}}_{k\left| {k - 1} \right.}^{j\left( {0,0} \right)} = {{F}}_{k - 1}^j{\tilde{ P}}_{k - 1\left| {k - 1} \right.}^{j\left( {0,0} \right)}{\left( {{{F}}_{k - 1}^j} \right)^{\rm{T}}} + {{Q}}_{k - 1}^j,\\ & {{P}}_{k\left| {k - 1} \right.}^{j\left( {i,i} \right)} = {\tilde{ P}}_{k - 1\left| {k - 1} \right.}^{j\left( {i - 1,i - 1} \right)},{{P}}_{k\left| {k - 1} \right.}^{j\left( {0,i} \right)} = {{F}}_{k - 1}^j{\tilde{ P}}_{k - 1\left| {k - 1} \right.}^{j\left( {0,i - 1} \right)} \end{split} $$ (9) (b) 新息残差及其协方差矩阵求解:
$$ \quad {{v}}_k^j = {{{z}}_k} - {{H}}_k^{a\left( j \right)}{{\hat X}}_{k\left| {k - 1} \right.}^j = {{{z}}_k} - {{H}}_k^j{{\hat X}}_{k\left| {k - 1} \right.}^{j\left( 0 \right)}\!\!\!\!\!\! $$ (10) $$ \quad {{S}}_k^j = {{H}}_k^j{{\hat X}}_{k\left| {k - 1} \right.}^{j\left( 0 \right)}{{P}}_{k\left| {k - 1} \right.}^{j\left( {0,0} \right)}{\left( {{{H}}_k^j} \right)^{\rm{T}}} + {{{R}}_k} $$ (11) 因此,模型
$m_k^j$ 的似然函数为:$\varLambda _k^j = {\cal{N}}\left( {{{v}}_k^j\left| {0,{{S}}_k^j} \right.} \right)$ 。(c) 滤波更新:
滤波增益为
$$ {{K}}_k^{j\left( i \right)} \!=\! {\left( {{{P}}_{k\left| {k - 1} \right.}^{j\left( {0,i} \right)}} \right)^{\rm{T}}}{\left( {{{H}}_k^j} \right)^{\rm{T}}}{\left( {{{S}}_k^j} \right)^{ - 1}}, i = 0,1, ··· ,L $$ (12) 状态估计及其协方差矩阵为
$$ \quad{{\hat X}}_{k\left| k \right.}^{j\left( i \right)} = {{\hat X}}_{k\left| {k - 1} \right.}^{j\left( i \right)} + {{K}}_k^{j\left( i \right)}{{v}}_k^j,\;\;i = 0,1, ··· ,L $$ (13) $$ \quad {{P}}_{k\left| k \right.}^{j\left( {i,i} \right)} = {{P}}_{k\left| {k - 1} \right.}^{j\left( {i,i} \right)} - {{K}}_k^{j\left( i \right)}{{S}}_k^j{\left( {{{K}}_k^{j\left( i \right)}} \right)^{\rm{T}}} $$ (14) $$ \quad {{P}}_{k\left| k \right.}^{j\left( {0,i} \right)} = {{P}}_{k\left| {k - 1} \right.}^{j\left( {0,i} \right)} - {{K}}_k^{j\left( 0 \right)}{{S}}_k^j{\left( {{{K}}_k^{j\left( i \right)}} \right)^{\rm{T}}} $$ (15) (3) 模型概率更新与矩匹配:
$$ \mu _k^j = P\left( {m_k^j\left| {,{{{z}}_{1:k}}} \right.} \right) = {{\mu _{k\left| {k - 1} \right.}^j\varLambda _k^j} / {\sum\nolimits_{i = 1}^M {\mu _{k\left| {k - 1} \right.}^i\varLambda _k^i} }} $$ (16) $$ {{{\hat X}}_{k\left| k \right.}} = \sum\limits_{j = 1}^M {\mu _k^j{{\hat X}}_{k\left| k \right.}^j} $$ (17) $$ {{{P}}_{k\left| k \right.}} \!=\! \sum\limits_{j = 1}^M {\mu _k^j\!\left[ {{{P}}_{k\left| k \right.}^j \!+\! {{\left( {{{\hat X}}_{k\left| k \right.}^j \!-\! {{{{\hat X}}}_{k\left| k \right.}}} \right)}^{\rm{T}}}\!\!\left( {{{\hat X}}_{k\left| k \right.}^j \!-\! {{{{\hat X}}}_{k\left| k \right.}}} \right)} \right]} $$ (18) 因此,平滑的状态
${{{\hat x}}_{k - L\left| k \right.}}$ 和协方差矩阵${{{P}}_{k - L\left| k \right.}}$ 为${{{\hat x}}_{k - L\left| k \right.}} = {{\hat X}}_{k\left| k \right.}^{\left( L \right)}$ ,${{{P}}_{k - L\left| k \right.}} = {{P}}_{k\left| k \right.}^{\left( {L,L} \right)}$ 。(4) 平滑模型概率:
平滑模型概率
$P\left( {m_{k - i}^j\left| {,{{{z}}_{1:k}}} \right.} \right)$ 满足$$ \begin{split} & P\left( {m_{k - i}^j\left| {,{{{z}}_{1:k}}} \right.} \right) \\ & \approx \frac{{{\cal{N}}\left( {{{{{\hat x}}}_{k - i\left| k \right.}}\left| {{{\hat x}}_{k - i\left| {k - i} \right.}^j,{{P}}_{k - i\left| {k - i} \right.}^j} \right.} \!\!\right)\!\!P\!\left( {m_{k - i}^j\left| {,{{{z}}_{1:k - i}}} \right.} \!\right)}}{{\displaystyle\sum\limits_{i = 1}^n {{\cal{N}}\!\!\left( {{{{{\hat x}}}_{k - i\left| k \right.}}\!\left| {{{\hat x}}_{k - i\left| {k - i} \right.}^i,{{P}}_{k - i\left| {k - i} \right.}^i} \right.} \!\!\right)\!\!P\!\left( {m_{k - i}^i\left| {,{{{z}}_{1:k - i}}} \right.} \!\right)} }} \end{split} $$ (19) 4. 自适应的增广状态-交互式多模型算法
4.1 联合后验密度的变分贝叶斯推断概率模型
对于自适应的AS-IMM算法需要估计联合后验密度函数
$p\left( {{{{X}}_k},{{{R}}_k}\left| {{{{z}}_{1:k}}} \right.} \right)$ 。对于跳变线性马尔科夫系统,需要引入隐变量${{{I}}_k}$ ,定义见文献[15]。${{{I}}_k}$ 表示模型的辨识因子,满足$I_k^j \in \left\{ 0\ 1 \right\}$ ,$j \in \left[ 1\ M \right]$ ,且$\displaystyle\sum\nolimits_{j = 1}^M {I_k^j = 1}$ 。因此,算法需要求解后验密度为$p\left( {{{{X}}_k},{{{R}}_k},{{{I}}_k}\left| {{{{z}}_{1:k}}} \right.} \right)$ 。在传统的变分贝叶斯推断当中,当后验密度满足平均场假设,且概率模型满足共轭性时,后验密度可通过固定点迭代求得解析结果[11]。而对于增广的状态空间模型,本文得出定理1。定理1 对于跳变线性马尔科夫系统,给定增广的状态空间模型式(4)和式(5),后验密度
$p\left( {{{{X}}_k},{{{R}}_k},{{{I}}_k}\left| {{{{z}}_{1:k}}} \right.} \right)$ 的变分推断概率模型非共轭。证明 根据平均场假设,后验密度可以近似为
$$p\left( {{{{X}}_k},{{{R}}_k},{{{I}}_k}\left| {{{{z}}_{1:k}}} \right.} \right) \approx q\left( {{{{X}}_k}} \right)q\left( {{{{R}}_k}} \right)q\left( {{{{I}}_k}} \right)$$ (20) 此时,后验密度可通过如下的固定点迭代求解:
$$ \ln q\left( {{{{I}}_k}} \right) = {{\rm{E}}_{{{{X}}_k},{{{R}}_k}}}\left[ {\ln p\left( {{{{X}}_k},{{{R}}_k},{{{z}}_k},{{{I}}_k}\left| {{{{z}}_{1:k - 1}}} \right.} \right)} \right] + {C_{{{{I}}_k}}} $$ (21) $$ \ln q\left( {{{{X}}_k}} \right) = {{\rm{E}}_{{{{I}}_k},{{{R}}_k}}}\left[ {\ln p\left( {{{{X}}_k},{{{R}}_k},{{{z}}_k},{{{I}}_k}\left| {{{{z}}_{1:k - 1}}} \right.} \right)} \right] + {C_{{{{X}}_k}}} $$ (22) $$ \ln q\left( {{{{R}}_k}} \right) = {{\rm{E}}_{{{{X}}_k},{{{I}}_k}}}\left[ {\ln p\left( {{{{X}}_k},{{{R}}_k},{{{z}}_k},{{{I}}_k}\left| {{{{z}}_{1:k - 1}}} \right.} \right)} \right] + {C_{{{{R}}_k}}} $$ (23) 其中,
${C_{{{{I}}_k}}}$ ,${C_{{{{X}}_k}}}$ 和${C_{{{{R}}_k}}}$ 分别表示关于${{{I}}_k}$ ,${{{X}}_k}$ 和${{{R}}_k}$ 的常数。给定对数联合似然函数$$ \begin{split} & \ln p\left( {{{{X}}_k},{{{R}}_k},{{{z}}_k},{{{I}}_k}\left| {{{{z}}_{1:k - 1}}} \right.} \right) = \ln p\left( {{{{z}}_k}\left| {{{{X}}_k},{{{R}}_k},{{{I}}_k}} \right.} \right) \\ & \quad + \ln p\left( {{{{X}}_k}\left| {{{{I}}_k},{{{z}}_{1:k - 1}}} \right.} \right) + \ln p\left( {{{{I}}_k}\left| {{{{z}}_{1:k - 1}}} \right.} \right)\\ & \quad + \ln p\left( {{{{R}}_k}\left| {{{{z}}_{1:k - 1}}} \right.} \right)\\[-10pt] \end{split} $$ (24) 其中,给定
${{{X}}_{k - 1}}$ 的高斯混合先验分布,满足$$ \begin{split} p\left( {{{{X}}_{k - 1}}\left| {{{{z}}_{1:k - 1}}} \right.} \right) =\,& \sum\limits_{j = 1}^M p\left( {{{{X}}_{k - 1}}\left| {m_{k - 1}^j,{{{z}}_{1:k - 1}}} \right.} \right)\\ & \cdot P\left( {m_{k - 1}^j\left| {{{{z}}_{1:k - 1}}} \right.} \right)\\[-15pt] \end{split} $$ (25) $p\left( {{{{X}}_k}\left| {{{{I}}_k},{{{z}}_{1:k - 1}}} \right.} \right)$ 和$p\left( {{{{I}}_k}\left| {{{{z}}_{1:k - 1}}} \right.} \right)$ 可由交互多模型的重初始化和模型条件预测求解得到,即$$p\left( {{{{X}}_k}\left| {{{{I}}_k},{{{z}}_{1:k - 1}}} \right.} \right) = \prod\limits_{j = 1}^M {{\cal{N}}{{\left( {{{{X}}_k}\left| {{{\hat X}}_{k\left| {k - 1} \right.}^j,{{P}}_{k\left| {k - 1} \right.}^j} \right.} \right)}^{I_k^j}}} $$ (26) $$ p\left( {{{{I}}_k}\left| {{{{z}}_{1:k - 1}}} \right.} \right) = \prod\limits_{j = 1}^M {{{\left( {\mu _{k\left| {k - 1} \right.}^j} \right)}^{I_k^j}}} $$ (27) 为了满足共轭性,先验分布
$p\left( {{{{R}}_{k - 1}}\left| {{{{z}}_{1:k - 1}}} \right.} \right)$ 取为逆-威沙特分布[14],即$p\left( {{{{R}}_{k - 1}}\left| {{{{z}}_{1:k - 1}}} \right.} \right) = {\cal{I}}{\cal{W}}\left( {{{\hat v}_{k - 1}},{{{{\hat V}}}_k}} \right)$ 。通过引入遗忘因子,$p\left( {{{{R}}_k}\left| {{{{z}}_{1:k - 1}}} \right.} \right)$ 可由Beta-Bartlett演化模型[14]给出,即$p\left( {{{{R}}_k}\left| {{{{z}}_{1:k - 1}}} \right.} \right) = {\cal{I}}{\cal{W}} \left( {{{{R}}_k}\left| {{{\hat v}_{k\left| {k - 1} \right.}},{{{\hat{ V}}}_{k\left| {k - 1} \right.}}} \right.} \right)$ ,其中$$ \qquad {\hat v_{k\left| {k - 1} \right.}} = \lambda {\hat v_{k - 1}} + \left( {1 - \lambda } \right)\left( {{n_z} - 1} \right) $$ (28) $$ \qquad {{{\hat V}}_{k\left| {k - 1} \right.}} = \lambda {{{\hat V}}_k} $$ (29) 其中,
$\lambda $ 表示遗忘因子,取值在$\left( {0.9,1} \right]$ 范围内。根据式(22),$q\left( {{{{X}}_k}} \right)$ 由式(30)决定$$ \begin{split} \ln q\left( {{{{X}}_k}} \right) =\,& {{\rm{E}}_{{{{I}}_k},{{{R}}_k}}}\left[ {\ln p\left( {{{{X}}_k},{{{R}}_k},{{{z}}_k},{{{I}}_k}\left| {{{{z}}_{1:k - 1}}} \right.} \right)} \right] + {C_{{{{X}}_k}}} \\ =\,& \sum\limits_{j = 1}^M \mu _k^j\left[ \ln \left( {{\cal{N}}\left( {{{{z}}_k}\left| {{{H}}_k^{a\left( j \right)}{{{X}}_k},{\rm{E}}{{\left[ {{{R}}_k^{ - 1}} \right]}^{ - 1}}} \right.} \right)} \right) \right.\\ & \left.+ \ln \left( {{\cal{N}}\left( {{{{X}}_k}\left| {{{\hat X}}_{k\left| {k - 1} \right.}^j,{{P}}_{k\left| {k - 1} \right.}^j} \right.} \right)} \right) \right] \\ =\,& \sum\limits_{j = 1}^M {\mu _k^j\ln {\cal{N}}\left( {{{{X}}_k}\left| {{{\hat X}}_{k\left| k \right.}^j,{{P}}_{k\left| k \right.}^j} \right.} \right)} \\[-21pt] \end{split} $$ (30) 其中
$$ {{P}}_{k\left| k \right.}^j = {\left( {{{\left[ {{{P}}_{k\left| {k - 1} \right.}^j} \right]}^{ - 1}} + {{\left[ {{{H}}_k^{a\left( j \right)}} \right]}^{\rm{T}}}{\rm{E}}\left[ {{{R}}_k^{ - 1}} \right]{{H}}_k^{a\left( j \right)}} \right)^{ - 1}} $$ (31) $$ \begin{split} {{\hat X}}_{k\left| k \right.}^j =\,& {{P}}_{k\left| k \right.}^j\left[ {{\left( {{{P}}_{k\left| {k - 1} \right.}^j} \right)}^{ - 1}}{{\hat X}}_{k\left| {k - 1} \right.}^j \right.\\ & \left.+ {{\left( {{{H}}_k^{a\left( j \right)}} \right)}^{\rm{T}}}{\rm{E}}\left[ {{{R}}_k^{ - 1}} \right]{{{z}}_k} \right] \end{split} $$ (32) 显然,式(31)和式(32)表明,
${{\hat X}}_{k\left| k \right.}^j$ 和${{P}}_{k\left| k \right.}^j$ 是由模型匹配的增广卡尔曼滤波求解与式(12)至式(15)是一致的,其中量测协方差矩阵${{{R}}_k}$ 的取值为${{{R}}_k} = {{\rm{E}}_{q\left( {{{{R}}_k}\left| {{{{z}}_{1:k - 1}}} \right.} \right)}}{\left[ {{{R}}_k^{ - 1}} \right]^{ - 1}}$ 。由式(30)可知,$q\left( {{{{X}}_k}} \right)$ 不是高斯混合分布。同理可证明,$q\left( {{{{I}}_k}} \right)$ 和$q\left( {{{{R}}_k}} \right)$ 也不满足共轭性。因而,后验密度$p\left( {{{{X}}_k},{{{R}}_k},{{{I}}_k}\left| {{{{z}}_{1:k}}} \right.} \right)$ 的变分推断概率模型为非共轭的。证毕4.2 联合后验密度的近似变分求解
定理1表明了完全采用变分推断求解
$q\left( {{{{I}}_k}} \right)$ ,$q\left( {{{{X}}_k}} \right)$ 和$q\left( {{{{R}}_k}} \right)$ 是没有解析结果的,这与文献[15]的结论类似,其区别在于文献[15]并未考虑${{{R}}_k}$ 的估计,且系统模型为非增广的线性跳变马尔科夫系统。文献[15]提出了一种$q\left( {{{{I}}_k}} \right)$ 和$q\left( {{{{X}}_k}} \right)$ 的近似变分求解方法,该方法首先采用传统的IMM求解$q\left( {{{{I}}_k}} \right)$ ;而给定$q\left( {{{{I}}_k}} \right)$ 后,$q\left( {{{{X}}_k}} \right)$ 由变分贝叶斯方法求解。需要说明的是,文献[15]的研究对象为非增广状态的跳变马尔科夫模型,为了便于理解,本文并没有区分增广状态${{{X}}_k}$ 和非增广状态${{{x}}_k}$ 。受文献[15]启发,本文提出一种“信息反馈+后处理”的次优处理方法,所提方法基于以下事实:(1) 给定量测协方差矩阵
${{{R}}_k}$ 后,$q\left( {{{{I}}_k}} \right)$ 和$q\left( {{{{X}}_k}} \right)$ 可由传统的AS-IMM求解[8];(2) 在传统的AS-IMM当中,
$p\left( {{{{X}}_k}\left| {{{{z}}_{1:k}}} \right.} \right)$ 可通过矩匹配近似为高斯分布[8];(3) 在非跳变高斯状态空间模型中,
$q\left( {{{{R}}_k}} \right)$ 的变分贝叶斯推断是解析的[14]。基于上述观点,所提方法的实现步骤如下:
步骤1 求解隐变量和增广状态的后验分布。固定
$q\left( {{{{R}}_k}} \right)$ ,取定${{{R}}_k}$ 为:${{{R}}_k} = {{\rm{E}}_{p\left( {{{{R}}_k}\left| {{{{z}}_{1:k - 1}}} \right.} \right)}}{\left[ {{{R}}_k^{ - 1}} \right]^{ - 1}}$ ,通过AS-IMM求解$q\left( {{{{I}}_k}} \right)$ 和$q\left( {{{{X}}_k}} \right)$ 。其中,${{{I}}_k}$ 服从多项式分布,即$q\left( {{{{I}}_k}} \right) = {\rm{Mul}}\left( {{{{I}}_k}\left| {{\mu _k}} \right.} \right)$ ,$\mu _k^j \propto \mu _{k\left| {k - 1} \right.}^j \varLambda _k^j\left( {{{{z}}_k}} \right)$ ;$q\left( {{{{X}}_k}} \right)$ 的求解见第3节。步骤2 矩匹配近似。根据式(17)和式(18)中的矩匹配步骤,求解近似的后验分布
$q\left( {{{{X}}_k}} \right)$ 。步骤3 求解量测噪声协方差矩阵的后验分布。给定状态后验密度
$q\left( {{{{X}}_k}} \right) \approx {\cal{N}}\left( {{{{X}}_k}\left| {{{{{\hat X}}}_{k\left| k \right.}},{{{P}}_{k\left| k \right.}}} \right.} \right)$ ,$q\left( {{{{R}}_k}} \right)$ 可根据高斯条件下的变分方法求解,即$$ \ln q\left( {{{{R}}_k}} \right) = {{\rm{E}}_{q\left( {{{{X}}_k}} \right)}}\left[ {\ln p\left( {{{{X}}_k},{{{R}}_k},{{{z}}_k}\left| {{{{z}}_{1:k - 1}}} \right.} \right)} \right] + {C_{{{{R}}_k}}} $$ (33) 显然,
${{{R}}_k}$ 的后验密度服从逆-威沙特分布,即${{{R}}_k} \sim {\cal{I}}{\cal{W}}\left( {{v_k},{{{V}}_k}} \right)$ ,其中$$\quad \left. \begin{aligned} & {{{\hat v}_k} = {{\hat v}_{k\left| {k - 1} \right.}} + 1} \\ &{{{{{\hat V}}}_k} = {{{{\hat V}}}_{k\left| {k - 1} \right.}} + {{{A}}_k}} \end{aligned} \right\} $$ (34) $$ \begin{split} \quad{{{A}}_k} =\,& \left( {{{{z}}_k} - {{H}}_k^a{{{{\hat X}}}_{k\left| k \right.}}} \right){\left( {{{{z}}_k} - {{H}}_k^a{{{{\hat X}}}_{k\left| k \right.}}} \right)^{\rm{T}}} \\ & + {\left( {{{H}}_k^a} \right)^{\rm{T}}}{{{P}}_{k\left| k \right.}}{{H}}_k^a \end{split} $$ (35) 其中,
${{H}}_k^a$ 表示模型中任意取定的一个量测矩阵,即${{H}}_k^{a\left( j \right)},\forall j$ 。因而,${{{R}}_k}$ 的估计值为:${{{\hat R}}_k} \!=\! \dfrac{{{{{{\hat V}}}_{k\left| {k - 1} \right.}} + {{{A}}_k}}}{{{{\hat v}_{k\left| {k - 1} \right.}} + 1}}$ 。步骤4 求解平滑模型概率。根据式(19)进行模型概率平滑。
为了行文方便,本文将所提出的自适应的AS-IMM算法简记为VB-AS-IMM。对于VB-AS-IMM,本文有如下说明:(1)步骤1中固定
$q\left( {{{{R}}_k}} \right)$ 时,对${{R}}_k^{ - 1}$ 求期望是关于$p\left( {{{{R}}_k}\left| {{{{z}}_{1:k - 1}}} \right.} \right)$ 操作的,其代表着利用了将第$k - 1$ 时刻估计的协方差矩阵反馈到当第$k$ 时刻,因此步骤1包含了一种信息反馈处理。(2)所提方法可以看作一种后处理。首先,采用传统的AS-IMM进行状态估计;而后,给出矩匹配的估计值${{{\hat X}}_{k\left| k \right.}}$ 和${{{P}}_{k\left| k \right.}}$ 后,所提算法按照高斯条件求解$q\left( {{{{R}}_k}} \right)$ 。(3)传统的VB方法中,需要通过固定点迭代来求解各分量。由于$q\left( {{{{X}}_k}} \right)$ 的求解过程与传统的AS-IMM是一致的,因而对于$q\left( {{{{X}}_k}} \right)$ 的优化可一步实现,因而不用迭代的方式。5. 仿真实验
本节通过仿真实验验证VB-AS-IMM算法对机动目标的跟踪性能,并与现有的IMM[1]算法和AS-IMM[8]算法进行了对比。仿真实验采用匹配IMM(Matched IMM, MIMM)和匹配AS-IMM(Matched AS-IMM, MAS-IMM)分别作为非平滑和平滑多模型算法的最优性能基准,同时采用均方根误差(Root Mean Square Error, RMSE)来衡量状态估计性能。为了衡量算法对量测协方差矩阵的估计性能,本文引入了平均Kullback-Leibler散度(Averaged Kullback-Leibler Divergence, AKLD)指标,定义为
$${\rm{AKLD}} \triangleq \frac{1}{N}\sum\limits_{n = 1}^N {{D_{KL}}\left( {\hat p\left( {{{{n}}_k}} \right)\left\| {p\left( {{{{n}}_k}} \right)} \right.} \right)} $$ (36) 其中,
$ {D_{KL}}\left( \cdot \right)$ 表示求解两个概率分布的KLD;$ {\hat p\left( {{{{n}}_k}} \right)}$ 表示采用估计$ {{\hat{ R}}_k}$ 构成的量测噪声分布;$ {p\left( {{{{n}}_k}} \right)}$ 表示真实的量测噪声分布;N=100表示蒙特卡洛仿真次数。仿真实验考虑雷达对慢速机动目标的跟踪问题。坐标系选择为雷达为原点所构成的2维直角坐标系。系统状态向量为
${{{x}}_k} = {[{x_k}\;\;{\dot x_k}\;\;{y_k}\;\;{\dot y_k}]^{\rm{T}}}$ 。仿真实验考虑3个运动模型,分别为匀速(CV)模型,右转弯机动(CT1)模型和左转弯机动(CT2)模型。模型的非增广状态方程为CV模型:
$${{{x}}_k} = {{{F}}_{{\rm{CV}}}}{{{x}}_{k - 1}} + {{{Q}}_{{\rm{CV}}}},\;\;{{{Q}}_{{\rm{CV}}}} = \sigma _1^2{{{Q}}_1}$$ (37) CT1和CT2模型:
$${{{x}}_k} = {{{F}}_{{\rm{CT}}}}{{{x}}_{k - 1}} + {{{Q}}_{{\rm{CT}}}},\;\;{{{Q}}_{{\rm{CT}}}} = \sigma _2^2{{{Q}}_2}$$ (38) 公式中的
${{{F}}_{{\rm{CV}}}}$ ,${{{Q}}_1}$ ,${{{F}}_{{\rm{CT}}}}$ 和${{{Q}}_2}$ 可参考文献[18]。雷达对目标斜距和方位角进行测量,量测方程为$${{{z}}_k} = \left[ {\begin{array}{*{20}{c}} {\sqrt {{x_k}^2 + {y_k}^2} } \\ {\arctan \left(\dfrac{{{y_k}}}{{{x_k}}}\right)} \end{array}} \right] + {{{w}}_k}$$ (39) 其中,
${x_k}$ 和${y_k}$ 分别为目标在直角坐标系中的$X$ 轴坐标和$Y$ 轴坐标;量测噪声协方差矩阵${{{R}}_k} \!=\! {\rm{diag}}\left( {\left[ {\delta _{k,r}^2\; \delta _{k,\theta }^2} \right]} \right)$ ;$\delta _{k,r}^{}$ 为雷达测距误差标准差,$\delta _{k,\theta }^{}$ 为雷达测角误差标准差。仿真中,$\delta _{k,r}^{}$ 和$\delta _{k,\theta }^{}$ 设置为:$\delta _{k,r}^{} = 60\;{\rm{m}}$ ,$\delta _{k,\theta }^{} = {0.2^{ \circ} }$ 。目标的初始状态设置为${{{x}}_0} = \left[ {100\;{\rm{km}}}\ \ {20\;{\rm{m/s}}} \ \ {100\;{\rm{km}}} \ \ {0\;{\rm{m/s}}} \right]^{\rm{T}}$ 。由于量测方程为非线性方程,因此仿真中采用去偏转换的卡尔曼滤波(Converted Measurement Kalman Filter, CMKF)[19]。雷达扫描时间间隔$T = 5\;{\rm{s}}$ ,观测总时长为400帧。过程噪声设置为:$\sigma _1^2 = {10^{ - 4}}$ ,$\sigma _2^2 = {10^{ - 4}}$ 。目标在第1~80帧、161~240帧和321~400帧作匀速运动;在第81~160帧作右转弯运动,角速度为$w = {{ - {{0.45}^{ \circ} }} / {\rm{s}}}$ ;在第241~320帧作左转弯运动,角速度为$w = {{{{0.45}^{ \circ} }} / {\rm{s}}}$ 。仿真场景如图1所示。仿真实验考虑量测协方差矩阵中
$\delta _{k,r}^{}$ 和$\delta _{k,\theta }^{}$ 的先验信息并不准确,在VB-AS-IMM, IMM和AS-IMM中的量测协方差矩阵为真实去偏转换量测协方差矩阵的10倍。而在MIMM和MAS-IMM中,量测协方差矩阵为真实的去偏转换的量测协方差矩阵。各算法的模型初始概率均设置为$\left[ {{1 / 3}}\ \ {{1 / 3}}\ \ {{1 / 3}} \right]$ ,马尔科夫模型转移概率矩阵为$${{P}} = \left[ {\begin{array}{*{20}{c}} {0.950}&{0.025}&{0.025} \\ {0.025}&{0.950}&{0.025} \\ {0.025}&{0.025}&{0.950} \end{array}} \right]$$ (40) 在算法VB-AS-IMM中,协方差矩阵分布的初始参数设置为
${\hat v_0} = 5$ ,${{{\hat V}}_0} = 10{{R}}_p^0$ ,其中${{R}}_p^0$ 表示初始时刻真实的去偏转换量测协方差矩阵;遗忘因子设置为$\lambda = 0.98$ 。AS-IMM和VB-AS-IMM得平滑步长$L = 10$ 。图2给出了量测协方差矩阵先验信息不准确时算法跟踪性能。图2(a)给出了VB-AS-IMM的AKLD,图中的AKLD随时间收敛,最后约为0。这表明所提出的VB-AS-IMM算法能够在线学习量测协方差矩阵的统计分布,并将这一信息反馈至估计器来提高状态估计性能。图2(b), 2(c)给出了各算法的位置均方根误差和速度均方根误差。图3给出了各算法的模型概率,其中图3(a)为CV的模型概率,图3(b)为CT1的模型概率,图3(c)为CT2的模型概率。由图2(b)和图2(c)可看出,传统的IMM和MIMM算法在模型切换时,会呈现出均方根误差抬升的现象。通过分析IMM和MIMM的模型概率可知,造成RMSE抬升是由于滤波类算法的模型概率估计结果并不精确,这引起了模型竞争。由于采用了平滑处理,AS-IMM的性能优于滤波算法。此外,由图3可知,采用平滑处理能够显著提高模型概率的准确度且减少模型切换的过渡时间。由于存在模型失配,AS-IMM算法估计的模型概率较MAS-IMM来说精度更差,这也导致了算法的RMSE较MAS-IMM明显增大。由于VB-AS-IMM能够在线学习量测协方差矩阵的统计分布,因而并不依赖于先验信息,其性能与最优的MAS-IMM基本一致。
表1统计了在整个跟踪时间段内,各算法的平均RMSE。由表1可知,在量测协方差矩阵先验信息不准确时,IMM的位置RMSE较MIMM提高了60.22 m,速度RMSE提高了1 m/s。AS-IMM的位置RMSE较MAS-IMM提高了12.04 m,速度RMSE提高了0.52 m/s。VB-AS-IMM的RMSE较IMM有明显的改善,与AS-IMM相近。
表 1 算法的平均RMSE算法 位置RMSE(m) 速度RMSE(m/s) MIMM 131.69 2.37 IMM 191.91 3.37 MAS-IMM 78.90 0.95 AS-IMM 90.94 1.47 VB-AS-IMM 79.74 0.96 6. 结论
本文针对现有的增广状态-交互式多模型算法存在着依赖于量测噪声协方差矩阵这一先验信息的问题,提出了一种自适应的增广状态-交互式多模型算法VB-AS-IMM,本文方法引入了一种“信息反馈+后处理”的近似方案,实现了目标状态和未知量测噪声协方差矩阵的联合求解。仿真结果表明,相比于非平滑的IMM算法,采用平滑处理的AS-IMM和VB-AS-IMM能显著改善跟踪精度。这主要表现在两个方面:一是平滑算法利用了更多的量测新息;二是平滑算法能显著提高模型概率的准确性且减少模型切换的过渡时间。本文提出的VB-AS-IMM能够在线学习未知量测噪声协方差矩阵的统计分布,相较于传统的IMM和AS-IMM,算法具有鲁棒性和较高的跟踪精度。因此,VB-AS-IMM更适用于量测协方差矩阵先验未知或不准确条件下的机动目标跟踪问题。
-
表 1 算法的平均RMSE
算法 位置RMSE(m) 速度RMSE(m/s) MIMM 131.69 2.37 IMM 191.91 3.37 MAS-IMM 78.90 0.95 AS-IMM 90.94 1.47 VB-AS-IMM 79.74 0.96 -
BLOM H A P and BAR-SHALOM Y. The interacting multiple model algorithm for systems with Markovian switching coefficients[J]. IEEE Transactions on Automatic Control, 1988, 33(8): 780–783. doi: 10.1109/9.1299 CHANG C B and ATHANS M. State estimation for discrete systems with switching parameters[J]. IEEE Transactions on Aerospace and Electronic Systems, 1978, AES–14(3): 418–425. doi: 10.1109/TAES.1978.308603 ANDERSON B D O and MOORE J B. Optimal Filtering[M]. Englewood Cliffs: Prentice-Hall, 1979: 165–190. RAUCH H. Solutions to the linear smoothing problem[J]. IEEE Transactions on Automatic Control, 1963, 8(4): 371–372. doi: 10.1109/TAC.1963.1105600 KELLY C N and ANDERSON B D O. On the stability of fixed-lag smoothing algorithms[J]. Journal of the Franklin Institute, 1971, 291(4): 271–281. doi: 10.1016/0016-0032(71)90183-9 MOORE J B. Discrete-time fixed-lag smoothing algorithms[J]. Automatica, 1973, 9(2): 163–174. doi: 10.1016/0005-1098(73)90071-X MATHEWS V J and TUGNAIT J K. Detection and estimation with fixed lag for abruptly changing systems[J]. IEEE Transactions on Aerospace and Electronic Systems, 1983, AES–19(5): 730–739. doi: 10.1109/TAES.1983.309374 CHEN Bing and TUGNAIT J K. Interacting multiple model fixed-lag smoothing algorithm for Markovian switching systems[J]. IEEE Transactions on Aerospace and Electronic Systems, 2000, 36(1): 243–250. doi: 10.1109/7.826326 MORELANDE M R and RISTIC B. Smoothed state estimation for nonlinear Markovian switching systems[J]. IEEE Transactions on Aerospace and Electronic Systems, 2008, 44(4): 1309–1325. doi: 10.1109/TAES.2008.4667711 LOPEZ R and DANÈS P. Low-complexity IMM smoothing for jump Markov nonlinear systems[J]. IEEE Transactions on Aerospace and Electronic Systems, 2017, 53(3): 1261–1272. doi: 10.1109/TAES.2017.2669698 BLEI D M, KUCUKELBIR A, and MCAULIFFE J D. Variational inference: A review for statisticians[J]. Journal of The American Statistical Association, 2017, 112(518): 859–877. doi: 10.1080/01621459.2017.1285773 TZIKAS D G, LIKAS A C, and GALATSANOS N P. The variational approximation for bayesian inference[J]. IEEE Signal Processing Magazine, 2008, 25(6): 131–146. doi: 10.1109/MSP.2008.929620 SARKKA S and NUMMENMAA A. Recursive noise adaptive kalman filtering by variational bayesian approximations[J]. IEEE Transactions on Automatic Control, 2009, 54(3): 596–600. doi: 10.1109/TAC.2008.2008348 AGAMENNONI G, NIETO J I, and NEBOT E M. Approximate inference in state-space models with heavy-tailed noise[J]. IEEE Transactions on Signal Processing, 2012, 60(10): 5024–5036. doi: 10.1109/TSP.2012.2208106 MA Yanjun, ZHAO Shunyi, and HUANG Biao. Multiple-model state estimation based on variational Bayesian inference[J]. IEEE Transactions on Automatic Control, 2019, 64(4): 1679–1685. doi: 10.1109/TAC.2018.2854897 DONG Peng, JING Zhongliang, and LEUNG H. Variational bayesian adaptive Cubature information filter based on Wishart distribution[J]. IEEE Transactions on Automatic Control, 2017, 62(11): 6051–6057. doi: 10.1109/TAC.2017.2704442 XU Hong, XIE Wenchong, YUAN Huadong, et al. Fixed-point iteration Gaussian sum filtering estimator with unknown time-varying non-Gaussian measurement noise[J]. Signal Processing, 2018, 153: 132–142. doi: 10.1016/j.sigpro.2018.07.017 LI Xiaorong and JILKOV V P. Survey of maneuvering target tracking. Part I. Dynamic models[J]. IEEE Transactions on Aerospace and Electronic Systems, 2003, 39(4): 1333–1364. doi: 10.1109/TAES.2003.1261132 ZHOU Gongjian, GUO Zhengkun, CHEN Xi, et al. Statically fused converted measurement Kalman filters for phased-array radars[J]. IEEE Transactions on Aerospace and Electronic Systems, 2018, 54(2): 554–568. doi: 10.1109/TAES.2017.2760798 -