CN102707279B - 序列ud分解的多目标跟踪方法 - Google Patents

Abstract

本发明公开了一种序列UD分解的多目标跟踪方法，用于解决现有的目标跟踪方法数值结构不稳定导致雷达跟踪过程目标失跟的技术问题。技术方案是对估计误差方差阵进行序列UD分解，建立数值稳定结构模型，没有误差估计的方差阵中有两个半正定矩阵相减，在有限字长的处理系统中能够保证不会产生含有负特征值的对称矩阵。通过对估计误差方差阵的序列UD分解，建立了数值稳定的多目标跟踪结构模型，避免了误差估计的方差阵中两个半正定矩阵相减，在有限字长的处理系统中不会出现数值发散，从而保证了目标跟踪系统的可靠性、避免了雷达跟踪过程目标失跟和整个雷达系统性错误。

Description

序列UD分解的多目标跟踪方法

技术领域

本发明涉及一种雷达多目标跟踪方法，特别涉及一种序列UD分解的多目标跟踪方法，属于信息技术领域。

背景技术

多目标跟踪技术在军用及民用领域均有广泛的应用，可用于空中目标检测、跟踪与攻击，空中导弹防御，空中交通管制，港口和海洋监视等。近年来，随着战场环境的改变，对抗和反对抗技术的发展，产生了背景强杂波、低信噪比、低检测概率和高虚警率等一系列问题，对多目标跟踪方法的精度和准确性提出了更高的要求。

多目标跟踪的目的是将探测器所接收到的量测对应不同的信息源，形成不同观测集合或轨迹，根据轨迹估计被跟踪目标的数目以及每一目标的运动参数，实现对多个目标的跟踪。用于多目标状态估计的基本滤波方法有α-β滤波、α-β-γ滤波、卡尔曼滤波、扩展卡尔曼滤波、高斯和近似、最优非线性滤波、粒子滤波和自适应滤波等。α-β和α-β-γ滤波器由于结构简单，计算量小，在早期计算机资源短缺时应用很广。卡尔曼滤波是多目标跟踪的一种基本方法，但是需要知道系统的精确数学模型，并且只适用于线性系统，限制了算法的应用。扩展卡尔曼滤波将卡尔曼滤波理论扩展到非线性领域，用一个高斯分布来近似状态的条件概率分布；而当近似条件不满足时，高斯和滤波器则用一个高斯分布的加权和来近似状态的条件概率分布。最优非线性滤波使用Makov转移概率来描述目标的动力学过程，具有很好的特性，但是计算量较大，因此一直没有得到广泛应用。粒子滤波采用随机采样，由于计算量太大和粒子退化问题，不适合实际应用。为了改进粒子滤波，无迹卡尔曼滤波采用确定性采样，使得采样的粒子点个数减少，避免了粒子滤波中的粒子点退化问题，因此其应用领域很广。自适应滤波方法通过对目标机动的检测，实时调整滤波器参数或增加滤波器的状态，使滤波器实时适应目标运动，特别适合对机动目标的跟踪；目前，在实际雷达跟踪系统最常用的仍然为JPDA(Joint Probabilistic Data Association，联合概率数据关联)方法(JamesA.Roecker，A Class of Near Optimal JPDA Algorithms，IEEE TRANSACTIONS ONAEROSPACE AND ELECTRONIC SYSTEMS，1994，VOL.30(2)：504-51O)，其它方法大多数是对JPDA方法的简化等。然而，JPDA等方法误差估计的方差阵中有两个半正定矩阵相减，在有限字长的处理系统中会产生含有正负特征值的对称矩阵，导致雷达跟踪过程目标失跟和整个雷达系统性错误。

发明内容

为了解决现有目标跟踪方法数值结构不稳定导致雷达跟踪过程目标失跟的技术缺陷，本发明提供一种序列UD分解的多目标跟踪方法，该方法在在多目标跟踪的时间和测量更新中，通过对估计误差方差阵进行序列UD分解，建立数值稳定结构模型，没有误差估计的方差阵中有两个半正定矩阵相减，在有限字长的处理系统中能够保证不会产生含有负特征值的对称矩阵，可以避免雷达跟踪过程目标失跟和整个雷达系统性错误。

本发明解决其技术问题采用的技术方案是，一种序列UD分解的多目标跟踪方法，其特征包括以下步骤：

1、定义N个目标跟踪中第i个目标的离散化模型为

x_i(k+1)＝Φ(k+1，k)x_i(k)+Λω_i(k)，

其中：

为状态向量，(x，y，z)为目标在地面参考直角坐标系下的位置坐标，ω_i(k)表示方差为对角阵Q_i(k)的过程噪声向量，Q_i(k-1)＝diag[q_i1(k-1)，q_i2(k-1)，q_i3(k-1)]，q_il(k-1)∈R^1×1(l＝1，2，3)为标量；Φ(k+1，k)＝Φ＝diag[Φ₁，Φ₁，Φ₁]为状态转移矩阵， $\mathrm{\Λ}={\∫}_{\mathrm{kT}}^{(k+1)T}\mathrm{\Φ}(k+1,\mathrm{\τ})\mathrm{\Γ}\left(\mathrm{\τ}\right)\mathrm{d\τ}=\left[\begin{array}{ccc}{\mathrm{\Λ}}_1& 0& 0\\ 0& {\mathrm{\Λ}}_1& 0\\ 0& 0& {\mathrm{\Λ}}_1\end{array}\right],$ Γ(t)为系数矩阵， $\mathrm{\Γ}=\left[\begin{array}{ccc}{\mathrm{\Γ}}_1& 0& 0\\ 0& {\mathrm{\Γ}}_1& 0\\ 0& 0& {\mathrm{\Γ}}_1\end{array}\right],$ Γ₁＝[0 0 1]^T， ${\mathrm{\Φ}}_1=\left[\begin{array}{ccc}1& T& \frac{1}{2}{T}^2\\ 0& 1& T\\ 0& 0& 1\\ & & \end{array}\right],$ ${\mathrm{\Λ}}_1={\left[\begin{array}{ccc}\frac{1}{6}{T}^3& \frac{1}{2}{T}^2& T\end{array}\right]}^T,$ T为采样周期；

第i个目标的时间更新为：

x_i(k/k-1)＝Φx_i(k-1/k-1)

${\stackrel{\‾}{U}}_i(k/k-1)={\mathrm{\ΦU}}_i(k-1/k-1)$

${\stackrel{\‾}{D}}_i(k/k-1)=D_i(k-1/k-1)$

U_i(k/k-1)＝[y₄ y₅ … y₁₂]，D_i(k/k-1)＝diag[d₄，d₅，…，d₁₂]

其中：x_i(k/k-1)为第i个目标在kT时刻的一步预测值，为对应的一步预测误差的方差阵，初始的y_j为Y^T(k/k-1)的第j列向量，d_j为

的第j行第j列对角元素； $Y(k/k-1)=\left[\begin{array}{cc}{\stackrel{\‾}{U}}_i(k/k-1)& \mathrm{\Λ}\end{array}\right],$ $\tilde{D}(k/k-1)=\left[\begin{array}{cc}{D}_i(k-1/k-1)& 0\\ 0& Q_i(k-1)\end{array}\right]$ U_i(k/k-1)为单位上三角矩阵，D_i(k/k-1)为对角矩阵；←表示改写，例如a←b为用b值改写a；初始条件为x_i(0/0)和 $U_i(0/0)D_i(0/0)U_i^T(0/0)=P_i(0/0);$

2、第i个目标观测方程为：z_i(k)＝g_i[x_i(k)]+v_i(k)

其中：z_i(k)为对第i个目标的r维观测向量，g_i[x_i(k)]为对应的输出，v_i(k)表示方差为对角阵R_i(k)＝diag[R_i1，R_i2，…，R_ir]测量噪声；计算

$x_i(k/k)=x_i(k/k-1)+G_i\left(k\right)\{\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\λ}}_{\mathrm{ij}}\left(k\right)z_{\mathrm{ij}}\left(k\right)-g_i[x_i(k/k-1)\left]\right\}$

其中：

为单位上三角矩阵，

为对角矩阵，D_ill(k/k-1)为D_i(k/k-1)的第l行第l列元素，U_ijl(k/k-1)为为U_i(k/k-1)的第j行第l列元素；z_ij(k)为雷达对第i个目标的第j(j＝1，2，…，η)个回波，x_i(k/k)为第i个目标kT时刻的滤波值，λ_ij(k)为权系数，且：

h_ij(k)为

的j列(j＝1，2，…，r)向量，←表示改写，例如a←b为用b值改写a；

3、第i个跟踪估计方法为：

U_i(k/k)＝[y_1+η y_2+η … y_9+η]，D_i(k/k)＝diag[d_1+η，d_2+η，…，d_9+η]

其中：←表示改写，例如c←b为用b值改写c；U_i(k/k)为单位上三角矩阵，D_i(k/k)为对角矩阵；初始的y_j为A^T(k)的第j列向量， $A\left(k\right)=\left[\begin{array}{cc}\stackrel{\‾}{U}\left(k\right)& G_i\left(k\right)d^T(I-{\mathrm{\Ωuu}}^T)\end{array}\right],$ $\tilde{D}\left(k\right)=\left[\begin{array}{cc}\stackrel{\‾}{D}\left(k\right)& 0\\ 0& \mathrm{\Ω}\end{array}\right],$

$u=\left[\begin{array}{c}1\\ 1\\ .\\ .\\ .\\ 1\end{array}\right],$ $d=\left[\begin{array}{c}{\mathrm{\Δ}}_{i,1}^T\left(k\right)\\ {\mathrm{\Δ}}_{i,2}^T\left(k\right)\\ .\\ .\\ .\\ {\mathrm{\Δ}}_{i,m}^T\left(k\right)\end{array}\right];$ Δ_i，j(k)为第j个候选回波信息向量，

Δ_i，j(k)＝z_i，j(k)-g_i[x_i(k/k-1)]。

本发明的有益结果是：该方法在在多目标跟踪的时间和测量更新中，通过对估计误差方差阵的三次序列UD分解，建立了数值稳定的多目标跟踪结构模型，避免了误差估计的方差阵中两个半正定矩阵相减，在有限字长的处理系统中不会出现数值发散，从而保证了多目标跟踪方法的可靠性、避免了雷达跟踪过程目标失跟和整个雷达系统性错误。下面结合实例对本发明作详细说明。

具体实施方式

1、定义N个目标跟踪中第i个目标的离散化模型为

x_i(k+1)＝Φ(k+1，k)x_i(k)+Λω_i(k)，

其中：

为状态向量，(x，y，z)为目标在地面参考直角坐标系下的位置坐标，ω_i(k)表示方差为对角阵Q_i(k)的过程噪声向量，Q_i(k-1)＝diag[q_i1(k-1)，q_i2(k-1)，q_i3(k-1)]，q_il(k-1)∈R^1×1(l＝1，2，3)为标量；Φ(k+1，k)＝Φ＝diag[Φ₁，Φ₁，Φ₁]为状态转移矩阵， $\mathrm{\Λ}={\∫}_{\mathrm{kT}}^{(k+1)T}\mathrm{\Φ}(k+1,\mathrm{\τ})\mathrm{\Γ}\left(\mathrm{\τ}\right)\mathrm{d\τ}=\left[\begin{array}{ccc}{\mathrm{\Λ}}_1& 0& 0\\ 0& {\mathrm{\Λ}}_1& 0\\ 0& 0& {\mathrm{\Λ}}_1\end{array}\right],$ Γ(t)为系数矩阵， $\mathrm{\Γ}=\left[\begin{array}{ccc}{\mathrm{\Γ}}_1& 0& 0\\ 0& {\mathrm{\Γ}}_1& 0\\ 0& 0& {\mathrm{\Γ}}_1\end{array}\right],$ Γ₁＝[0 0 1]^T， ${\mathrm{\Φ}}_1=\left[\begin{array}{ccc}1& T& \frac{1}{2}{T}^2\\ 0& 1& T\\ 0& 0& 1\\ & & \end{array}\right],$ ${\mathrm{\Λ}}_1={\left[\begin{array}{ccc}\frac{1}{6}{T}^3& \frac{1}{2}{T}^2& T\end{array}\right]}^T,$ T为采样周期；

第i个目标的时间更新为：

x_i(k/k-1)＝Φx_i(k-1/k-1)

${\stackrel{\‾}{U}}_i(k/k-1)={\mathrm{\ΦU}}_i(k-1/k-1)$

${\stackrel{\‾}{D}}_i(k/k-1)=D_i(k-1/k-1)$

U_i(k/k-1)＝[y₄ y₅ … y₁₂]，D_i(k/k-1)＝diag[d₄，d₅，…，d₁₂]

其中：x_i(k/k-1)为第i个目标在kT时刻的一步预测值，为对应的一步预测误差的方差阵，初始的y_j为Y^T(k/k-1)的第j列向量，d_ij为

的第j行第j列对角元素； $Y(k/k-1)=\left[\begin{array}{cc}{\stackrel{\‾}{U}}_i(k/k-1)& \mathrm{\Λ}\end{array}\right],$ $\tilde{D}(k/k-1)=\left[\begin{array}{cc}{D}_i(k-1/k-1)& 0\\ 0& Q_i(k-1)\end{array}\right]$ U_i(k/k-1)为单位上三角矩阵，D_i(k/k-1)为对角矩阵；D_ll(k/k-1)为D_i(k/k-1)的第l行第l列元素，U_jl(k/k-1)为为U_i(k/k-1)的第j行第l列元素；←表示改写，例如a←b为用b值改写a；初始条件为x_i(0/0)和 $U_i(0/0)D_i(0/0)U_i^T(0/0)=P_i(0/0);$

2、第i个目标观测方程为：z_i(k)＝g_i[x_i(k)]+v_i(k)

其中：z_i(k)为对第i个目标的r维观测向量，g_i[x_i(k)]为对应的输出，v_i(k)表示方差为对角阵R_i(k)＝diag[R_i1，R_i2，…，R_ir]测量噪声；计算

$x_i(k/k)=x_i(k/k-1)+G_i\left(k\right)\{\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\λ}}_{\mathrm{ij}}\left(k\right)z_{\mathrm{ij}}\left(k\right)-g_i[x_i(k/k-1)\left]\right\}$

其中：

为单位上三角矩阵，

为对角矩阵，D_ill(k/k-1)为D_i(k/k-1)的第l行第l列元素，U_ijl(k/k-1)为为U_i(k/k-1)的第j行第l列元素；z_ij(k)为雷达对第i个目标的第j(j＝1，2，…，η)个回波，x_i(k/k)为第i个目标kT时刻的滤波值，λ_ij(k)为权系数，且：

取g_i[x_i(k)]＝[r_i(k) α_i(k) β_i(k)]^T，r_i为雷达能测量斜距、α_i为高低角、β_i方位角，且

$\left\{\begin{array}{c}{r}_i=\sqrt{x_i^2+y_i^2+z_i^2}\\ {\mathrm{\α}}_i={\tan }^{-1}\frac{z_i}{\sqrt{x_i^2+y_i^2}}\\ {\mathrm{\β}}_i={\tan }^{-1}\frac{x_i}{y_i}\end{array}\right.$

$H_i\left(k\right)=\frac{{\∂g}_i\left[x_i\left(k\right)\right]}{{\∂x}_i\left(k\right)}{|}_{x_i\left(k\right)=x_i(k/k-1)}$

$={\left[\begin{array}{ccccccccc}\frac{x_i}{\sqrt{x_i^2+y_i^2+z_i^2}}& 0& 0& \frac{y_i}{\sqrt{x_i^2+y_i^2+z_i^2}}& 0& 0& \frac{z_i}{\sqrt{x_i^2+y_i^2+z_i^2}}& 0& 0\\ \frac{{-x}_i{z}_i}{(x_i^2+y_i^2+z_i^2)\sqrt{x_i^2+y_i^2}}& 0& 0& \frac{{-y}_i{z}_i}{(x_i^2+y_i^2+z_i^2)\sqrt{x_i^2+y_i^2}}& 0& 0& \frac{\sqrt{x_i^2+y_i^2}}{(x_i^2+y_i^2+z_i^2)}& 0& 0\\ \frac{y_i}{x_i^2+y_i^2}& 0& 0& \frac{{-x}_i}{x_i^2+y_i^2}& 0& 0& 0& 0& 0\end{array}\right]}_{x_i\left(k\right)=x_i(k/k-1)};$

3、第i个跟踪估计方法为：

U_i(k/k)＝[y_1+η y_2+η … y_9+η]，D_i(k/k)＝diag[d_1+η，d_2+η，…，d_9+η]

其中：←表示改写，例如c←b为用b值改写c；U_i(k/k)为单位上三角矩阵，D_i(k/k)为对角矩阵；初始的y_j为A^T(k)的第j列向量， $A\left(k\right)=\left[\begin{array}{cc}\stackrel{\‾}{U}\left(k\right)& G_i\left(k\right)d^T(I-{\mathrm{\Ωuu}}^T)\end{array}\right],$ $\tilde{D}\left(k\right)=\left[\begin{array}{cc}\stackrel{\‾}{D}\left(k\right)& 0\\ 0& \mathrm{\Ω}\end{array}\right],$

$u=\left[\begin{array}{c}1\\ 1\\ .\\ .\\ .\\ 1\end{array}\right],$ $d=\left[\begin{array}{c}{\mathrm{\Δ}}_{i,1}^T\left(k\right)\\ {\mathrm{\Δ}}_{i,2}^T\left(k\right)\\ .\\ .\\ .\\ {\mathrm{\Δ}}_{i,m}^T\left(k\right)\end{array}\right];$ Δ_i，j(k)为第j个候选回波信息向量，

Δ_i，j(k)＝z_i，j(k)-g_i[x_i(k/k-1)]。

Claims (1)

1.一种序列UD分解的多目标跟踪方法，其特征在于包括以下步骤：

（1）、定义N个目标跟踪中第i个目标的离散化模型为

x_i(k+1)=Φ(k+1,k)x_i(k)+Λω_i(k)，

其中：

为状态向量，（x,y,z）为目标在地面参考直角坐标系下的位置坐标，ω_i(k)表示方差为对角阵Q_i(k)的过程噪声向量，Q_i(k-1)=diag[q_i1(k-1),q_i2(k-1),q_i3(k-1)]，q_il(k-1)∈R^1×1(l=1,2,3)为标量；Φ(k+1,k)=Φ=diag[Φ₁,Φ₁,Φ₁]为状态转移矩阵，

$\mathrm{\Λ}={\∫}_{\mathrm{kT}}^{(k+1)T}\mathrm{\Φ}(k+1,\mathrm{\τ})\mathrm{\Γd\τ}=\left[\begin{array}{ccc}{\mathrm{\Λ}}_1& 0& 0\\ 0& {\mathrm{\Λ}}_1& 0\\ 0& 0& {\mathrm{\Λ}}_1\end{array}\right],$ Γ为系数矩阵， $\mathrm{\Γ}=\left[\begin{array}{ccc}{\mathrm{\Γ}}_1& 0& 0\\ 0& {\mathrm{\Γ}}_1& 0\\ 0& 0& {\mathrm{\Γ}}_1\end{array}\right],$ Γ₁=[0 0 1]^T， ${\mathrm{\Φ}}_1=\left[\begin{array}{ccc}1& T& \frac{1}{2}{T}^2\\ 0& 1& T\\ 0& 0& 1\end{array}\right],$ ${\mathrm{\Λ}}_1={\left[\begin{array}{ccc}\frac{1}{6}{T}^3& \frac{1}{2}{T}^2& T\end{array}\right]}^T,$ T为采样周期；

第i个目标的时间更新为：

$\begin{array}{c}{x}_i(k/k-1)=\mathrm{\Φ}{x}_i(k-1/k-1)\\ {\stackrel{\‾}{U}}_i(k/k-1)=\mathrm{\Φ}{U}_i(k-1/k-1)\\ {\stackrel{\‾}{D}}_i(k/k-1)=D_i(k-1/k-1)\end{array}$

U_i(k/k-1)=[y₄ y₅ … y₁₂],D_i(k/k-1)=diag[d₄,d₅,…,d₁₂]

其中：x_i(k/k-1)为第i个目标在kT时刻的一步预测值，为对应的一步预测误差的方差阵，初始的y_j为Y^T(k/k-1)的第j列向量，d_j为的第j行第j列对角元素； $Y(k/k-1)=\left[\begin{array}{cc}{\stackrel{\‾}{U}}_i(k/k-1)& \mathrm{\Λ}\end{array}\right],$ $\tilde{D}(k/k-1)=\left[\begin{array}{cc}{D}_i(k-1/k-1)& 0\\ 0& Q_i(k-1)\end{array}\right]$

U_i(k/k-1)为单位上三角矩阵，D_i(k/k-1)为对角矩阵；←表示改写，a←b为用b值改写a；初始条件为x_i(0/0)和 $U_i(0/0)D_i(0/0)U_i^T(0/0)=P_i(0/0);$

（2）、第i个目标观测方程为：z_i(k)=g_i[x_i(k)]+v_i(k)

其中：z_i(k)为对第i个目标的r维观测向量，g_i[x_i(k)]为对应的输出，v_i(k)表示方差为对角阵R_i(k)=diag[R_i1,R_i2,…,R_ir]的测量噪声；计算

$x_i(k/k)=x_i(k/k-1)+G_i\left(k\right)\{\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\λ}}_{i,j}\left(k\right)z_{i,j}\left(k\right)-g_i[x_i(k/k-1)\left]\right\}$

其中：

为单位上三角矩阵，

为对角矩阵，D_ill(k/k-1)为D_i(k/k-1)的第l行第l列元素，U_ijl(k/k-1)为U_i(k/k-1)的第j行第l列元素；z_i,j(k)为雷达对第i个目标的第j(j=1,2,…,η)个回波，x_i(k/k)为第i个目标kT时刻的滤波值，λ_i,j(k)为权系数，且： $\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\λ}}_{i,j}\left(k\right)=1,$ $H_i\left(k\right)={\frac{{\∂g}_i\left[x_i\left(k\right)\right]}{{\∂x}_i\left(k\right)}|}_{x_i\left(k\right)=x_i(k/k-1)},$ h_ij(k)为

的j列(j=1,2,…,r)向量，←表示改写，a←b则表示为用b值改写a；

（3）、第i个目标跟踪估计方法为：

U_i(k/k)=[y_1+η y_2+η … y_9+η],D_i(k/k)=diag[d_1+η,d_2+η,…,d_9+η]

其中：U_i(k/k)为单位上三角矩阵，D_i(k/k)为对角矩阵；初始的y_j为A^T(k)的第j列向量， $A\left(k\right)=\left[\begin{array}{cc}\stackrel{\‾}{U}\left(k\right)& G_i\left(k\right)d^T(1-\mathrm{\Ω}{\mathrm{uu}}^T)\end{array}\right],$ $\tilde{D}\left(k\right)=\left[\begin{array}{cc}\stackrel{\‾}{D}\left(k\right)& 0\\ 0& \mathrm{\Ω}\end{array}\right],$

$u=\left[\begin{array}{c}1\\ 1\\ .\\ .\\ .\\ 1\end{array}\right],$ $d=\left[\begin{array}{c}{\mathrm{\Δ}}_{i,1}^T\left(k\right)\\ {\mathrm{\Δ}}_{i,2}^T\left(k\right)\\ .\\ .\\ .\\ {\mathrm{\Δ}}_{i,m}^T\left(k\right)\end{array}\right];$

△_i,j(k)为第j个候选回波信息向量，

△_i,j(k)=z_i,j(k)-g_i[x_i(k/k-1)]。