CN102043906A - 一种机动目标跟踪的对角型融合方法 - Google Patents

Abstract

本发明提出了一种利用最小方差意义下最优的多传感器信息融合准则联合多个模型输出的对角型交互式多模型融合方法，以实现对带Markov切换跳变参数的机动性目标的精确跟踪。该方法避免了传统方法的概率密度和概率质量的混合，而且有效的区别了状态的不同维度产生的不同影响，经实验证明具有自适应能力强、估计精度高的特点。能迅速准确的跟踪上具有机动运动特性的目标，具有很高的实际价值。

Description

一种机动目标跟踪的对角型融合方法

技术领域

本发明涉及一类带Markov随机跳变参数的机动性目标的跟踪问题，具体地说是利用对角矩阵交互式的融合多个模型输出的方法实现对一类具有模型跳变性目标的精确跟踪。

背景技术

机动目标跟踪在军事和民用领域里有广阔的应用前景，在军事领域中，如反弹道导弹的防御，空防预警，地对空、空对空、舰对空、舰对舰的超视距多目标控测、跟踪与攻击，战场监视空地(海)多目标，精确制导和低空突防，火力控制、威胁估计、态势评估等。在民用领域中，如空中交通管制，机器人的道路规划和障碍躲避，电子医学及视频监控等。

在现今科学技术条件下，目标的速度、加速度与过去相比均有较大不同。高转弯率、超强加减速能力在众多打击利防御武器中已经出现。同时，非常规机动也迅速发展起来，目标能够完成更为复杂的机动运动，很多目标具有主动改变运动模式的能力，在不同的运动模式之间进行瞬时切换。因此，提出适用于具有模型跳变性目标的精确跟踪的方法是目标跟踪问题当前研究的一个热点。

传统的机动目标跟踪方法是自适应的滤波方法，这种方法是根据量侧值和滤波值之间的残差来检测机动是否发生，如发生则需要调整滤波器的参数。这类方法的缺点是既没有充分利用对过去状态估计值的修正信息，也没有充分利用过去量测值之间的关联信息。近年来较为流行的机动目标跟踪方法是多模型方法，这类方法利用半马尔科夫链来描述机动目标跟踪中结构的跳变。具有显著影响的方法有广义伪贝叶斯方法(GPB)和交互式多模型方法(IMM)，其中IMM方法被认为是一种最为有效的混合估计方案，且已经成为该领域的主流方法，然而，机动目标的跟踪涉及两个问题，一个是连续型值的参数估计，如目标状态，另一个是离散型假设的检测，如目标的运动模型。对于连续性值的随机过程通常由概率密度函数(PDF)来描述，而离散型的则通常用概率质量来表示，这就使得基于贝叶斯公式计方的IMM方法以及其它的MM方法不可避免的卷入了概率密度函数值与概率质量的混合计算。这种混合在一些情况下并不会引起很大的误差，但是在运用贝叶斯公式时，由两种不同量级的数值进行了混淆导致最后所获得的并不是一个真正的模型概率质量，而仅仅是一个近似。而且，对于目标跟踪而言，成功跟踪上目标的关键是使得跟踪误差达到最小，但IMM方法并不是以此为宗旨。

基于以上的分析，我们认为如果使用一些新的方法来联合多个模型，避免使用模型概率的近似值，跟踪精度将很可能超过经典的IMM方法以及其他各种多模型方法。这就是形成本发明的动机。

本发明利用最小方差意义下最优的多传感器信息融合准则而非传统上常用的贝叶斯理论用于联合多个模型的输出，提出了实现对一类具有模型跳变性目标精确跟踪的对角型交互式多模型(DIMM)融合方法。

发明内容

本发明的目的在于对带Markov切换跳变参数的离散时间随机系统，设计一种多模型融合方法，使得被跟踪目标的跟踪精度较已有的方法有较明显的提高。

本发明提出的对角型融合方法利用最小方差意义下最优的多传感器信息融合准则，将机动的模型集中获得的各种状态估计进行对角型交互式联合，以使目标的跟踪精度在最小方差意义下达到最优。技术方案具体表述如下：

考虑如下Markov跳变线性系统

X(k+1)＝F_jX(k)+G_jω(k)    (1)

Z(k)＝H_jX(k)+υ(k)        (2)

其中，X(k)为系统的n维状态变量，包括了目标的位置、速度、加速度等信息；Z(k)是观察目标时收到的量测值；j为目标机动运动的模型标示；ω(k)为系统的过程噪声，且满足：

E[ω_j(k)]＝0；E[ω_j(k)ω_j(k)^T]＝Q_j

这里E[·]表示的是随机变量的期望。υ(k)为系统的量测噪声，且满足

E[υ_j(k)]＝0；E[υ_j(k)υ_j(k)^T]＝R_j

目标由k-1时刻的运动模型i转变到k时刻的运动模型j的发生概率为：

ρ_ij＝Prob[M_j(k)|M_i(k-1)]

其中M_j(k)表示k时刻目标的运动模型为j这一事件，Prob表示概率。

为使目标的跟踪精度在最小方差意义下达到最优，给出如下的多传感器信息融合准则，即

命题1：如果

是状态变量X的一个无偏估计(为表示方便，这里省去了时间指标k)，估计误差记为融合估计记为

状态和估计的元素表示如下：

$X=\left[\begin{array}{c}{x}_1\\ .\\ .\\ .\\ x_n\end{array}\right],{\hat{X}}_j=\left[\begin{array}{c}{\hat{x}}_{j1}\\ .\\ .\\ .\\ {\hat{x}}_{\mathrm{jn}}\end{array}\right],{\hat{X}}_D=\left[\begin{array}{c}{\hat{x}}_{D1}\\ .\\ .\\ .\\ {\hat{x}}_{\mathrm{Dn}}\end{array}\right]$

则通过对角矩阵融合的最优估计为：

${\hat{X}}_D=\underset{j=1}{\overset{s}{\mathrm{\Σ}}}{B}_j{\hat{X}}_j$

这个估计是使得融合估计误差的协方差达到最小，即指标

达到最小。

其中：

B_j＝diag(b_j1，…，b_jn)，j＝1，…，s

b_ji来源于向量β_i(i＝1，…，n)

${\mathrm{\β}}_i=[b_{1i},b_{2i},...,b_{\mathrm{si}}]$

$=\frac{e^T{\left(P^i\right)}^{-1}}{e^T{\left(P^i\right)}^{-1}e}$

式中

其中

是矩阵的第i个对角线上的元素。

基于上述命题给出跟踪具体目标所获得的结果。图1给出的是被跟踪目标的位置航迹。图2显示的是目标在整个跟踪过程中的速度变化。考虑式(1)-(2)所述的Markov跳变线性系统，采用命题1中如下的多传感器信息融合准则，对各滤波器输出的状态估计进行对角型交互式联合，并使用如下的模型集对目标的机动性进行描述：

$F_i=\left[\begin{array}{cccc}1& T& 0& 0\\ 0& 1& 0& 0\\ 0& 1& 0& 0\end{array}\right];G_i=\left[\begin{array}{cc}\frac{T^2}{2}& 0\\ T& 0\\ 0& \frac{T^2}{2}\\ 0& T\end{array}\right];H_i=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\end{array}\right]$

噪声的协方差矩阵为：

$Q_i=\left[\begin{array}{cc}{q}_i& 0\\ 0& q_i\end{array}\right];R_i=\left[\begin{array}{cc}r& r/20\\ r/20& r\end{array}\right]$

其中q₁＝0.01，q₂＝50。r是x或y方向上的量测误差的协方差，采用间隔T＝10秒。目标的初始位置x(0)＝[2100；0；10000；-15]，目标的初始误差协方差矩阵定义为：

$P_i\left(0\right)=\left[\begin{array}{cccc}r& r/T& 0& 0\\ r/T& 2\×r/T^2& 0& 0\\ 0& 0& r& r/T\\ 0& 0& r/T& 2\×r/T^2\end{array}\right]$

两个常速模型之间有如下的Markov转移关系：

$\left[{\mathrm{\π}}_{\mathrm{ij}}\right]=\left[\begin{array}{cc}0.95& 0.05\\ 0.05& 0.95\end{array}\right]$

Markov链的初始分布为：

我们基于根均方误差(RMSE)来比较方法的性能。RMSE定义为：

${\mathrm{RMSE}}_k=\sqrt{\frac{1}{M}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{[{\hat{X}}^i\left(k\right)-X^i\left(k\right)]}^2}$

其中

和Xⁱ(k)分别表示第i次Monto Carlo运行获得的k时刻状态估计和目标更新的真实状态，M表示独立的Monto Calo运行的次数。

图3和图4显示了当量测噪声中参数r取为1000时，本发明提出的DIMM方法和经典的IMM方法分别在x和y方向上位置估计精度的比较。图5和图6显示的是两个方法分别在x和y方向上速度估计精度的比较。从图中可以看出，DIMM方法的速度估计明显优于经典的IMM方法，特别的速度估计上的优势，虽然DIMM方法第二次机动发生时位置估计误差有稍大的摆动，但是在很短的时间内，误差摆动迅速恢复稳定，出现这一现象是由于DIMM方法在高速机动发生时，对速度的迅速改变给予了更多的关注。

本发明利用对角矩阵作为各滤波器输出的状态估计的联合权重，避免了传统方法的概率密度和概率质量的混合，而且有效的区别了状态的不同维度产生的不同影响。经实验证明具有自适应能力强、估计精度高的特点。能迅速准确的跟踪上具有机动运动特性的目标。

本发明有以下一些技术特征：

(1)使用多个模型的联合来描述以及跟踪具有机动性运动特点的目标；

(2)使用对角矩阵作为联合多个模型状态输出的权重。

附图说明

图1为被跟踪目标的位置航迹；

图2为被跟踪目标的速度变化曲线；

图3为DIMM方法与IMM方法在x方向上位置估计RMSE变化曲线；

图4为DIMM方法与IMM方法在y方向上位置估计RMSE变化曲线；

图5为DIMM方法与IMM方法在x方向上速度估计RMSE变化曲线；

图6为DIMM方法与IMM方法在y方向上速度估计RMSE变化曲线；

具体实施方式

对角型融合方法的执行步骤。

以前述的命题为基础对具有图1和图2运动特性的目标进行跟踪，可给出对角型融合方法的执行步骤如下：

步骤1：使用如下的模型集对目标的机动性进行描述：

$F_i=\left[\begin{array}{cccc}1& T& 0& 0\\ 0& 1& 0& 0\\ 0& 1& 0& 0\end{array}\right];G_i=\left[\begin{array}{cc}\frac{T^2}{2}& 0\\ T& 0\\ 0& \frac{T^2}{2}\\ 0& T\end{array}\right];H_i=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\end{array}\right]$

噪声的协方差矩阵为：

$Q_i=\left[\begin{array}{cc}{q}_i& 0\\ 0& q_i\end{array}\right];R_i=\left[\begin{array}{cc}r& r/20\\ r/20& r\end{array}\right]$

目标的初始位置为X(0)＝[2100；0；10000；-15]，初始误差协方差矩阵为：

$P_i\left(0\right)=\left[\begin{array}{cccc}r& r/T& 0& 0\\ r/T& 2\×r/T^2& 0& 0\\ 0& 0& r& r/T\\ 0& 0& r/T& 2\×r/T^2\end{array}\right]$

其中q₁＝0.01，q₂＝50，r＝1000。

步骤2：初始化Markov链的分布为：

且取两个模型之间的Markov转移关系为如下矩阵：

$\left[{\mathrm{\π}}_{\mathrm{ij}}\right]=\left[\begin{array}{cc}0.95& 0.05\\ 0.05& 0.95\end{array}\right]$

步骤3：初始化对角矩阵权重B_j(0)：

计算 $P^i\left(0\right)=\left[\begin{array}{cc}{P}_1^{\left(\mathrm{ii}\right)}\left(0\right)& 0\\ 0& P_2^{\left(\mathrm{ii}\right)}\left(0\right)\end{array}\right]$

其中

是k时刻矩阵

的第i个对角线上的元素，

计算 ${\mathrm{\β}}_i\left(0\right)=\frac{e^T{\left(P^i\left(0\right)\right)}^{-1}}{e^T{\left(P^i\left(0\right)\right)}^{-1}e}$

$=[b_{1i},b_{2i},...,b_{\mathrm{si}}]$

B_j(0)＝diag(b_j1，…，b_jn)，j＝1，2

其中 $e={\left[\begin{array}{c}1\\ .\\ .\\ .\\ 1\end{array}\right]}_{s\×1}$

步骤4：k＝1；

步骤5：计算与每个模型匹配的混合初始对角矩阵权重B_i|j(k|k)(i，j＝1，2)：

$B_{i|j}\left(k\right|k)=\frac{{\mathrm{\π}}_{\mathrm{ij}}{B}_i(k-1)}{\underset{i=1}{\overset{2}{\mathrm{\Σ}}}{\mathrm{\π}}_{\mathrm{ij}}{B}_i(k-1)}$

$=\left[\begin{array}{cc}\frac{{\mathrm{\π}}_{\mathrm{ij}}{b}_{i1}}{\underset{i=1}{\overset{2}{\mathrm{\Σ}}}{\mathrm{\π}}_{\mathrm{ij}}{b}_{i1}}& 0\\ 0& \frac{{\mathrm{\π}}_{\mathrm{ij}}{b}_{\mathrm{in}}}{\underset{i=1}{\overset{2}{\mathrm{\Σ}}}{\mathrm{\π}}_{\mathrm{ij}}{b}_{\mathrm{in}}}\end{array}\right]$

步骤6：计算与每个模型匹配的滤波器初始状态

和相应的协方差矩阵P_0j(k|k)：

${\hat{X}}_{0j}\left(k\right|k)=\underset{i=1}{\overset{2}{\mathrm{\Σ}}}{B}_{i|j}\left(k\right|k){\hat{X}}_i(k-1)$

$P_{0j}\left(k\right|k)=\underset{i=1}{\overset{2}{\mathrm{\Σ}}}{B}_{i|j}\left(k\right|k)\{P_i(k-1)+{[{\hat{X}}_i(k-1)-{\hat{X}}_{0j}\left(k\right|k)][{\hat{X}}_i(k-1)-{\hat{X}}_{0j}\left(k\right|k)]}^T\}$

步骤7：用Kalman滤波获得每个模型在k时刻的状态估计

以及相应的协方差矩阵P_j(k)：

${\hat{X}}_j\left(k\right|k-1)=F_j{\hat{X}}_{0j}\left(k\right|k)$

P_j(k|k-1)＝F_jP_0j(k|k)F_j ^T+G_jQ_jG_j ^T

$r_j\left(k\right)=Z\left(k\right)-H_j{\hat{X}}_j\left(k\right|k-1)$

S_j(k)＝H_jP_j(k|k-1)H_j ^T+R_j(k)

K_j(k)＝P_j(k|k-1)H_j ^TS_j(k)^-1

${\hat{X}}_j\left(k\right)={\hat{X}}_j\left(k\right|k-1)+K_j\left(k\right)r_j\left(k\right)$

P_j(k)＝(I-K_j(k)H_j)P_j(k|k-1)

步骤8：计算更新的对角矩阵权重B_j(k)：

$P^i\left(k\right)=\left[\begin{array}{cc}{P}_1^{\left(\mathrm{ii}\right)}\left(k\right)& 0\\ 0& P_2^{\left(\mathrm{ii}\right)}\left(k\right)\end{array}\right]$

${\mathrm{\β}}_i\left(k\right)=\frac{e^T{\left(P^i\left(k\right)\right)}^{-1}}{e^T{\left(P^i\left(k\right)\right)}^{-1}e}$

$=[b_{1i},b_{2i},...,b_{\mathrm{si}}]$

B_j(k)＝diag(b_j1，…，b_jn)

步骤9：依据更新的对角矩阵权重计算融合估计

${\hat{X}}_D\left(k\right)=\underset{j=1}{\overset{2}{\mathrm{\Σ}}}{B}_j\left(k\right){\hat{X}}_j\left(k\right)$

步骤10：若k＝100，则跟踪结束，否则执行步骤11；

步骤11：令

执行步骤5至步骤10。

图3和图4显示了当量测噪声中参数r取为1000时，本发明提出的DIMM方法和经典的IMM方法分别在x和y方向上位置估计精度的比较。图5和图6显示的是两个方法分别在x和y方向上速度估计精度的比较。从图中可以看出，DIMM方法的速度估计明显优于经典的IMM方法，特别是在速度估计上的优势，虽然DIMM方法第二次机动发生时位置估计误差有稍大的摆动，但是在很短的时间内，误差摆动迅速恢复稳定，出现这一现象是由于DIMM方法在高速机动发生时，对速度的迅速改变给予了更多的关注。

本发明利用对角矩阵作为各滤波器输出的状态估计的联合权重，避免了传统方法的概率密度和概率质量的混合，而且有效的区别了状态的不同维度产生的不同影响。经实验证明具有自适应能力强、估计精度高的特点。能迅速准确的跟踪上具有机动运动特性的目标。

Claims (3)

1.一种机动目标跟踪的对角型融合方法，该方法包括：

步骤1：使用如下的模型集对目标的机动性进行描述：

$F_i=\left[\begin{array}{cccc}1& T& 0& 0\\ 0& 1& 0& 0\\ 0& 1& 0& 0\end{array}\right];G_i=\left[\begin{array}{cc}\frac{T^2}{2}& 0\\ T& 0\\ 0& \frac{T^2}{2}\\ 0& T\end{array}\right];H_i=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\end{array}\right]$

噪声的协方差矩阵为

$Q_i=\left[\begin{array}{cc}{q}_i& 0\\ 0& q_i\end{array}\right];R_i=\left[\begin{array}{cc}r& r/20\\ r/20& r\end{array}\right]$

目标的的初始位置X(0)＝[2100；0；10000；-15]，初始误差协方差矩阵为：

$P_i\left(0\right)=\left[\begin{array}{cccc}r& r/T& 0& 0\\ r/T& 2\×r/T^2& 0& 0\\ 0& 0& r& r/T\\ 0& 0& r/T& 2\×r/T^2\end{array}\right]$

其中q₁＝0.01，q₂＝50，r＝1000。

步骤2：初始化Markov链的分布为：

且取两个模型之间的Markov转移关系为如下矩阵：

$\left[{\mathrm{\π}}_{\mathrm{ij}}\right]=\left[\begin{array}{cc}0.95& 0.05\\ 0.05& 0.95\end{array}\right]$

步骤3：初始化对角矩阵权重B_j(0)：

计算 $P^i\left(0\right)=\left[\begin{array}{cc}{P}_1^{\left(\mathrm{ii}\right)}\left(0\right)& 0\\ 0& P_2^{\left(\mathrm{ii}\right)}\left(0\right)\end{array}\right]$

其中

是k时刻矩阵

的第i个对角线上的元素，

计算 ${\mathrm{\β}}_i\left(0\right)=\frac{e^T{\left(P^i\left(0\right)\right)}^{-1}}{e^T{\left(P^i\left(0\right)\right)}^{-1}e}$

$=[b_{1i},b_{2i},...,b_{\mathrm{si}}]$

B_j(0)＝diag(b_j1，…，b_jn)，j＝1，2

其中 $e={\left[\begin{array}{c}1\\ .\\ .\\ .\\ 1\end{array}\right]}_{s\×1}$

步骤4：k＝1

步骤5：计算与每个模型匹配的混合初始对角矩阵权重B_i|j(k|k)(i，j＝1，2)：

$B_{i|j}\left(k\right|k)=\frac{{\mathrm{\π}}_{\mathrm{ij}}{B}_i(k-1)}{\underset{i=1}{\overset{2}{\mathrm{\Σ}}}{\mathrm{\π}}_{\mathrm{ij}}{B}_i(k-1)}$

$=\left[\begin{array}{cc}\frac{{\mathrm{\π}}_{\mathrm{ij}}{b}_{i1}}{\underset{i=1}{\overset{2}{\mathrm{\Σ}}}{\mathrm{\π}}_{\mathrm{ij}}{b}_{i1}}& 0\\ 0& \frac{{\mathrm{\π}}_{\mathrm{ij}}{b}_{\mathrm{in}}}{\underset{i=1}{\overset{2}{\mathrm{\Σ}}}{\mathrm{\π}}_{\mathrm{ij}}{b}_{\mathrm{in}}}\end{array}\right]$

步骤6：计算与每个模型匹配的滤波器初始状态

和相应的协方差矩阵P_0j(k|k)：

${\hat{X}}_{0j}\left(k\right|k)=\underset{i=1}{\overset{2}{\mathrm{\Σ}}}{B}_{i|j}\left(k\right|k){\hat{X}}_i(k-1)$

$P_{0j}\left(k\right|k)=\underset{i=1}{\overset{2}{\mathrm{\Σ}}}{B}_{i|j}\left(k\right|k)\{P_i(k-1)+{[{\hat{X}}_i(k-1)-{\hat{X}}_{0j}\left(k\right|k)][{\hat{X}}_i(k-1)-{\hat{X}}_{0j}\left(k\right|k)]}^T\}$

步骤7：用Kalman滤波获得每个模型在k时刻的状态估计

以及相应的协方差矩阵P_j(k)

${\hat{X}}_j\left(k\right|k-1)=F_j{\hat{X}}_{0j}\left(k\right|k)$

P_j(k|k-1)＝F_jP_0j(k|k)F_j ^T+G_jQ_jG_j ^T

$r_j\left(k\right)=Z\left(k\right)-H_j{\hat{X}}_j\left(k\right|k-1)$

S_j(k)＝H_jP_j(k|k-1)H_j ^T+R_j(k)

K_j(k)＝P_j(k|k-1)H_j ^TS_j(k)^-1

${\hat{X}}_j\left(k\right)={\hat{X}}_j\left(k\right|k-1)+K_j\left(k\right)r_j\left(k\right)$

P_j(k)＝(I-K_j(k)H_j)P_j(k|k-1)

步骤8：计算更新的对角矩阵权重B_j(k)：

$P^i\left(k\right)=\left[\begin{array}{cc}{P}_1^{\left(\mathrm{ii}\right)}\left(k\right)& 0\\ 0& P_2^{\left(\mathrm{ii}\right)}\left(k\right)\end{array}\right]$

${\mathrm{\β}}_i\left(k\right)=\frac{e^T{\left(P^i\left(k\right)\right)}^{-1}}{e^T{\left(P^i\left(k\right)\right)}^{-1}e}$

$=[b_{1i},b_{2i},...,b_{\mathrm{si}}]$

B_j(k)＝diag(b_j1，…，b_jn)

步骤9：依据更新的对角矩阵权重计算融合估计：

${\hat{X}}_D\left(k\right)=\underset{j=1}{\overset{2}{\mathrm{\Σ}}}{B}_j\left(k\right){\hat{X}}_j\left(k\right)$

步骤10：若k＝100，则跟踪结束，否则执行步骤11；

步骤11：令

执行步骤5至步骤10。

2.根据权利要求1所述的方法，其特征在于，使用多个模型的联合来描述以及跟踪具有机动性运动特点的目标。

3.根据权利要求1所述的方法，其特征在于，使用对角矩阵作为联合多个模型状态输出的权重。

Images