CN103487800B - 基于残差反馈的多模型高速高机动目标跟踪方法 - Google Patents

Abstract

本发明公开一种基于残差反馈的多模型高速高机动目标跟踪方法，主要解决现有技术计算量大的问题。其步骤是：1)初始化雷达观测值，得到初始状态向量和初始协方差矩阵；2)选择若干运动模型，确定运动模型的状态转移矩阵；3)由运动模型的个数假定初始权值；4)利用状态向量和状态转移矩阵作一步预测；5)利用初始权值对一步预测值加权求和，得到状态估计向量和估计协方差矩阵；6)根据观测值和状态估计向量得到预测误差；7)由预测误差和估计协方差矩阵计算增益矩阵；8)利用增益矩阵和状态估计向量更新状态向量和状态协方差矩阵；9)利用预测误差和状态估计向量更新权值。本发明跟踪性能稳定、计算量小，可用于对临近空间目标的跟踪。

Description

基于残差反馈的多模型高速高机动目标跟踪方法

技术领域

本发明属于雷达技术领域，具体涉及一种对强机动目标的跟踪方法，可用于对临近空间目标的跟踪。

背景技术

由于目标跟踪在国防和民用各领域的重要价值，一直以来很多学者致力于研究目标跟踪算法。随着现代航空的迅速发展，各种飞行器的机动性和运动速度越来越高，但是可靠而精确的跟踪目标始终是目标跟踪系统设计的主要目的与难点。目标跟踪实际上就是对目标状态的跟踪滤波问题，而实现目标跟踪首先要使所建立的目标运动模型与实际的目标运动模型匹配，建立目标的运动模型和自适应滤波是机动目标的两个关键部分。

在当前目标机动能力日益增强的情况下，单模型算法很难满足跟踪速度和精度的需求，而多模型算法可以避免采用单模型时由于目标机动而造成模型的不准确，提高机动目标跟踪性能，从而实现对强机动目标的精确跟踪。经过国内外学者几十年来的不断研究，提出了很多目标模型与算法。目前广泛采用的算法是1984-1989年Blom和Bar-Shalom提出的交互式多模型IMM算法，该算法具有Markov转移概率，且算法中有多个模型并行工作，多个滤波器交互作用得到目标状态估计的结果，由于该模型有较好的自适应能力，因而比较理想。但是IMM算法中由于其模型转移概率是人为设置的，因此引入了人为误差，影响了跟踪性能。除此之外，因为IMM算法需要对每一个模型都进行滤波，计算量过大，对目标跟踪会有延迟。

发明内容

本发明的目的在于针对上述已有技术的不足，提出一种基于残差反馈的多模型高速高机动目标跟踪方法，以降低计算复杂度，提高雷达在目标强机动下的跟踪精度。

实现本发明的技术思路是：利用LMS算法，通过信号的预测残差调整滤波器的权值，从而使估计信号逼近期望信号，以提高跟踪性能，其实现步骤包括如下：

(1)雷达对接受到的机动目标的位置检测值进行N次采样，得到一个长度为N的观测值序列{Z(k)}，k=1，2，…，N；

(2)从观测值序列{Z(k)}中取出前三个采样点的观测值Z(1)，Z(2)，Z(3)，并利用差分法初始化机动目标的运动状态，得到机动目标在第3个采样点的状态向量和初始协方差矩阵P(3)；

(3)根据机动目标的运动特性选择M个目标运动模型，得到M个目标运动模型的状态转移矩阵F_i，i=1，2，…M，M≥2；

(4)根据选择的目标运动模型个数M，确定目标运动模型在第3采样点的权值：

W(3)=[1/M 1/M … 1/M]_1×M，

其中，[·]_1×M表示该矩阵为1行M列的矩阵；

(5)利用步骤(2)中得到第k采样点的状态向量和步骤(3)中得到的目标运动模型的状态转移矩阵F_i进行一步预测，得到对应模型的一步预测值 ${\hat{X}}_i(k+1/k):$

${\hat{X}}_i(k+1/k)=F_i\hat{X}\left(k\right)$

其中，表示第i个目标运动模型在第k采样点预测的机动目标在第k+1采样点的状态向量；

(6)利用步骤(4)中得到的权值对步骤(5)中得到的各个目标运动模型的一步预测值和协方差矩阵P(k)进行加权求和，得到基于LMS算法的状态估计向量和估计协方差矩阵P(k+1/k)：

$\hat{X}(k+1/k)=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{w}_i\left(k\right){\hat{X}}_i(k+1/k)$

$P(k+1/k)=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{w}_i\left(k\right)(F_{i}P\left(k\right)F_i^T+Q_i\left(k\right))$

其中，w_i(k)表示在第k采样点权值W(k)的第i个元素；[·]^T表示矩阵的转置；P(k)表示机动目标在第k采样点的协方差矩阵；Q_i(k)表示第i个目标运动模型在第k采样点的过程噪声协方差；

(7)根据步骤(1)中得到的雷达对机动目标位置的观测值Z(k+1)和步骤(6)中得到的基于LMS算法的状态估计向量通过下式求解机动目标在第k+1采样点的预测误差

$\tilde{Z}(k+1)=Z(k+1)-\hat{Z}(k+1/k)$

其中，其表示第k采样点预测机动目标在第k+1采样点的观测值；Z(k+1)表示雷达在第k+1采样点对机动目标位置的观测值；H表示机动目标的量测矩阵；

(8)由步骤(7)中得到的预测误差和步骤(6)中得到的估计协方差矩阵P(k+1/k)确定预测误差的协方差矩阵S(k+1)，继而得到增益矩阵G(k+1)：

S(k+1)=HP(k+1/k)H^T+R(k+1)

G(k+1)=P(k+1/k)H^TS^-1(k+1)，

其中，S(k+1)是预测误差在第k个采样点的协方差矩阵；R(k+1)是观测噪声在第k+1采样点的协方差矩阵；G(k+1)表示第k+1采样点的滤波增益矩阵；[·]^-1表示矩阵的逆；

(9)利用步骤(8)中得到的增益矩阵G(k+1)和步骤(6)中得到的状态估计向量完成对状态向量和协方差矩阵P(k)的更新，即：

$\hat{X}(k+1)=\hat{X}(k+1/k)+G(k+1)\tilde{Z}(k+1)$

P(k+1)=P(k+1/k)-G(k+1)S(k+1)G^T(k+1)

其中，表示机动目标更新后在第k+1采样点的状态向量；P(k+1)表示机动目标在第k+1采样点的协方差矩阵；[·]^T表示矩阵的转置；

(10)通过步骤(7)中得到的预测误差和步骤(5)中得到的目标运动模型的一步预测值对权值W(k)进行更新，得到更新后的权值：

$W(k+1)=W\left(k\right)+2\mathrm{\μ}{\tilde{Z}}^T(k+1)H\left[\begin{array}{cccc}{\hat{X}}_1(k+1/k)& {\hat{X}}_2(k+1/k)& \·\·\·& {\hat{X}}_M(k+1/k)\end{array}\right]$

其中，W(k+1)表示第k+1采样点目标运动模型的权值，μ表示步长因子，M为采用的目标运动模型的个数；

(11)判断跟踪是否完成，即判断k是否满足k≤N-1，若条件满足，则k递增，并继续执行步骤(5)；若条件不满足，则跟踪过程结束。

本发明与现有技术相比的优点在于：

1.跟踪性能高

本发明中引入了类似IMM算法中转移概率的权值，但与传统IMM算法中转移概率固定不同的是，本发明中滤波器可以自适应调节权值，从而提高了跟踪精度，减小了误差，同时提高了对机动目标的跟踪性能。

2.计算简单

现有IMM算法中每次滤波时对每个目标运动模型都要进行卡尔曼滤波，即每次滤波过程需要进行M次卡尔曼滤波；本发明中每个目标运动模型只需要进行一步预测，每次滤波过程都只需要进行一次卡尔曼滤波，大幅度降低了计算复杂度。

附图说明

图1为本发明的流程图；

图2为本发明对机动目标跟踪的效果示意图；

图3为本发明与IMM方法对机动目标跟踪三维位置误差局部比较图；

具体实施方式

参照图1，本发明实现步骤如下：

步骤1，对运动状态和模型参数进行初始化。

1a)雷达对其接收到的机动目标的位置检测值进行N次采样，得到一个长度为N的观测值序列{Z(k)}，k=1，2，…，N；

1b)从观测值序列中取出前三个采样点的观测值Z(1)，Z(2)，Z(3)，其中：

Z(1)=[xx(1) yy(1) zz(1)]^T

Z(2)=[xx(2) yy(2) zz(2)]^T

Z(3)=[xx(3) yy(3) zz(3)]^T

式中xx(·)，yy(·)，zz(·)分别表示机动目标在x轴、y轴、z轴的位置，[·]^T表示矩阵的转置；

1c)利用差分法对Z(1)，Z(2)，Z(3)进行运算，获得初始状态向量如下：

$\hat{X}\left(3\right)=\left[\begin{array}{c}\mathrm{xx}\left(3\right)\\ (\mathrm{xx}\left(3\right)-\mathrm{xx}\left(2\right))/T\\ ((\mathrm{xx}\left(3\right)-\mathrm{xx}\left(2\right))/T-(\mathrm{xx}\left(2\right)-\mathrm{xx}\left(1\right))/T)/T\\ \mathrm{yy}\left(3\right)\\ (\mathrm{yy}\left(3\right)-\mathrm{yy}\left(2\right))/T\\ ((\mathrm{yy}\left(3\right)-\mathrm{yy}\left(2\right))/T-(\mathrm{yy}\left(2\right)-\mathrm{yy}\left(1\right))/T)/T\\ \mathrm{zz}\left(3\right)\\ (\mathrm{zz}\left(3\right)-\mathrm{zz}\left(2\right))/T\\ ((\mathrm{zz}\left(3\right)-\mathrm{zz}\left(2\right)/T)-(\mathrm{zz}\left(3\right)-\mathrm{zz}\left(2\right)/T))/T\end{array}\right],$

其中，(xx(3)-xx(2))/T、(yy(3)-yy(2))/T、(zz(3)-zz(2))/T分别表示机动目标在第3采样点处的x轴、y轴、z轴方向速度；((xx(3)-xx(2))/T-(xx(2)-xx(1))/T)/T、((yy(3)-yy(2))/T-(yy(2)-yy(1))/T)/T、((zz(3)-zz(2)/T)-(zz(3)-zz(2)/T))/T分别表示机动目标在第3采样点的x轴、y轴、z轴方向加速度；T为采样间隔；

1d)根据初始状态向量获得初始协方差矩阵P(3)：

$P\left(3\right)=E\left[\hat{X}\left(3\right){\hat{X}}^T\left(3\right)\right],$

其中，E[·]表示数学期望；[·]^T表示矩阵的转置。

步骤2，选取目标运动模型，确定状态转移矩阵。

常用的目标运动模型包括二阶匀速CV模型，三阶匀加速CA模型，协调转弯CT模型，singer模型和“当前”统计模型，根据机动目标的运动特性，在常用的目标运动模型中选取M个目标运动模型，得到M个目标运动模型的状态转移矩阵F_i，i＝1，2，…M，M≥2。

步骤3，由目标运动模型的个数，确定目标运动模型在第3采样点的权值W(3)：

W(3)＝[1/M 1/M … 1/M]_1×M，

其中，[·]_1×M表示该矩阵为1行M列的矩阵。

步骤4，利用状态向量和状态转移矩阵进行一步预测，得到对应模型的一步预测值

${\hat{X}}_i(k+1/k)=F_i\hat{X}\left(k\right),$

其中，表示第i个目标运动模型在第k个采样点预测的机动目标在第k+1采样点的状态向量；表示机动目标在第k个采样点的状态向量；F_i表示第i个目标运动模型的状态转移矩阵。

步骤5，计算基于LMS算法的最优状态估计向量和估计协方差矩阵。

5a)对各个目标运动模型的一步预测值进行加权求和，得到基于LMS算的最优状态估计向量

$\hat{X}(k+1/k)=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{w}_i\left(k\right){\hat{X}}_i(k+1/k),$

其中，w_i(k)表示在第k采样点权值W(k)的第i个元素；

5b)对协方差矩阵P(k)进行加权求和，得到基于LMS算法的估计协方差矩阵P(k+1/k)：

$P(k+1/k)=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{w}_i\left(k\right)(F_{i}P\left(k\right){F_i}^T+Q_i\left(k\right)),$

其中，P(k)表示机动目标在第k采样点的协方差矩阵；Q_i(k)表示第i个目标运动模型在第k采样点的过程噪声协方差：

Q_i(k)＝C_i(k)qC_i(k)^T，

其中，C_i(k)表示第i个目标运动模型在第k采样点对应的控制矩阵，由选取的运动模型不同而取不同的矩阵；q为过程噪声方差，由机动目标位置的数量级确定。

步骤6，计算机动目标观测值的新息。

根据雷达对机动目标位置的观测值Z(k+1)和基于LMS算法的最优状态估计向量用下式求解机动目标预测值在第k+1采样点的预测误差 $\tilde{Z}(k+1):$

$\tilde{Z}(k+1)=Z(k+1)-\hat{Z}(k+1/k)$

其中，其表示第k采样点预测机动目标在第k+1采样点的观测值；Z(k+1)表示雷达在第k+1采样点对机动目标位置的观测值；H表示机动目标的量测矩阵：

$H=\left[\begin{array}{ccccccccc}1& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 1& 0& 0\end{array}\right].$

步骤7，计算增益矩阵G(k+1)。

7a)由估计协方差矩阵P(k+1/k)确定预测误差的协方差矩阵S(k+1)

S(k+1)＝HP(k+1/k)H^T+R(k+1)，

其中，S(k+1)是预测误差在第k个采样点的协方差矩阵；R(k+1)是观测噪声在第k+1采样点的协方差矩阵；

7b)由预测误差的协方差矩阵S(k+1)和估计协方差矩阵P(k+1/k)，计算增益矩阵G(k+1)：

G(k+1)＝P(k+1/k)H^TS^-1(k+1)，

其中，G(k+1)表示第k+1采样点的增益矩阵；[·]^-1表示矩阵的逆。

步骤8，对机动目标状态向量和状态协方差矩阵的进行更新。

8a)利用增益矩阵G(k+1)、状态估计向量和预测误差完成对状态向量的更新，即：

$\hat{X}(k+1)=\hat{X}(k+1/k)+G(k+1)\tilde{Z}(k+1),$

其中，表示机动目标更新后在第k+1采样点的状态向量；

8b)利用增益矩阵G(k+1)、估计协方差矩阵P(k+1/k)和预测误差的协方差矩阵S(k+1)完成对协方差矩阵P(k)的更新，即：

P(k+1)＝P(k+1/k)-G(k+1)S(k+1)G^T(k+1)

其中，P(k+1)表示机动目标更新后在第k+1采样点的协方差矩阵。

步骤9，通过预测误差和目标运动模型的一步预测值完成对权值W(k)的更新，即：

$W(k+1)=W\left(k\right)+2\mathrm{\μ}{\tilde{Z}}^T(k+1)H\left[\begin{array}{cccc}{\hat{X}}_1(k+1/k)& {\hat{X}}_2(k+1/k)& \·\·\·& {\hat{X}}_M(k+1/k)\end{array}\right]$

其中，W(k+1)表示更新后得到的第k+1采样点机动目标运动模型的权值，μ表示步长因子，保证算法收敛的μ的取值范围为：0＜μ＜1/λ_max，λ_max为矩阵x_in的自相关矩阵的最大特征值，

$x_{\mathrm{in}}=H\left[\begin{array}{cccc}{\hat{X}}_1(k+1/k)& {\hat{X}}_2(k+1/k)& \·\·\·& {\hat{X}}_M(k+1/k)\end{array}\right],$ M为采用的目标运动模型的个数。

步骤10，判断跟踪是否完成，即判断k是否满足k≤N-1，若条件满足，则k递增，并继续执行步骤(4)；反之，则跟踪过程结束。

本发明的效果可以通过以下仿真进一步说明：

1.仿真条件：

假设目标的起始位置为[50 125 75]km，起始速度为[0.7 0.65 0.6]km/s，x轴方向的加速度a1为0.15sin(0.04πt)km/s²，y轴方向的加速度a2为0.1cos(0.04πt)km/s²，z轴方向的加速度a3为0.12sin(0.02πt)km/s²。距离误差的标准差为方位角误差的标准差为仰角测量误差的标准差为对目标进行30s的观测，采样间隔T为0.25s，Monte Carlo仿真100次。

选择的运动模型是singer模型、“当前”统计模型和匀加速CA模型。IMM算法中的概率转移矩阵p为：

$p=\left[\begin{array}{ccc}0.8& 0.1& 0.1\\ 0.1& 0.8& 0.1\\ 0.1& 0.1& 0.8\end{array}\right].$

2.仿真内容

仿真1，采用本发明方法对仿真条件中所述的目标进行跟踪仿真，结果如图2。从图2可以看出，本发明方法的跟踪轨迹很接近机动目标运动的真实轨迹。

仿真2，对用本发明方法与传统IMM方法的跟踪误差进行比较，结果如图3。其中图3(a)表示本发明方法和传统IMM算法在x轴上的跟踪误差比较；图3(b)表示本发明方法和传统IMM算法在y轴上的跟踪误差比较；图3(c)表示本发明方法和传统IMM算法在z轴上的跟踪误差比较。

图3表明本发明方法的跟踪误差在整体趋势上小于IMM算法的跟踪误差，且跟踪稳定，更适合对高速高机动目标的跟踪。

Claims (1)

1.一种基于最小均方误差的强机动目标跟踪方法，包括如下步骤：

(1)雷达对接收到的机动目标的位置检测值进行N次采样，得到一个长度为N的观测值序列{Z(k)}，k＝1,2,",N；

(2)从观测值序列{Z(k)}中取出前三个采样点的观测值Z(1),Z(2),Z(3)，并利用差分法初始化机动目标的运动状态，得到机动目标在第3个采样点的状态向量和初始协方差矩阵P(3)；

(3)根据机动目标的运动特性选择M个目标运动模型，得到M个目标运动模型的状态转移矩阵F_i，i＝1,2,…M，M≥2；

(4)根据选择的目标运动模型个数M，确定目标运动模型在第3采样点的权值：

W(3)＝[1/M 1/M … 1/M]_1×M，

其中，[·]_1×M表示该矩阵为1行M列的矩阵；

(5)利用第k采样点的状态向量和步骤(3)中得到的目标运动模型的状态转移矩阵F_i进行一步预测，得到对应模型的一步预测值

${\hat{X}}_i(k+1/k)=F_i\hat{X}\left(k\right)$

其中，表示第i个目标运动模型在第k采样点预测的机动目标在第k+1采样点的状态向量；

(6)利用权值对步骤(5)中得到的各个目标运动模型的一步预测值和协方差矩阵P(k)进行加权求和，得到基于LMS算法的状态估计向量和估计协方差矩阵P(k+1/k)：

$\hat{X}(k+1/k)=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{w}_i\left(k\right){\hat{X}}_i(k+1/k)$

$P(k+1/k)=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{w}_i\left(k\right)(F_{i}P\left(k\right){F_i}^T+Q_i\left(k\right))$

其中，w_i(k)表示在第k采样点权值W(k)的第i个元素；[·]^T表示矩阵的转置；P(k)表示机动目标在第k采样点的协方差矩阵；Q_i(k)表示第i个目标运动模型在第k采样点的过程噪声协方差；

(7)根据步骤(1)中得到的雷达对机动目标位置的观测值Z(k+1)和步骤(6)中得到的基于LMS算法的状态估计向量通过下式求解机动目标在第k+1采样点的预测误差

$\tilde{Z}(k+1)=Z(k+1)-\hat{Z}(k+1/k)$

其中，其表示第k采样点预测机动目标在第k+1采样点的观测值；Z(k+1)表示雷达在第k+1采样点对机动目标位置的观测值；H表示机动目标的量测矩阵；

(8)由步骤(7)中得到的预测误差和步骤(6)中得到的估计协方差矩阵P(k+1/k)确定预测误差的协方差矩阵S(k+1)，继而得到增益矩阵G(k+1)：

S(k+1)＝HP(k+1/k)H^T+R(k+1)

G(k+1)＝P(k+1/k)H^TS^-1(k+1)，

其中，S(k+1)是预测误差在第k个采样点的协方差矩阵；R(k+1)是观测噪声在第k+1采样点的协方差矩阵；G(k+1)表示第k+1采样点的滤波增益矩阵；[·]^-1表示矩阵的逆；

(9)利用步骤(8)中得到的增益矩阵G(k+1)和步骤(6)中得到的状态估计向量完成对状态向量和协方差矩阵P(k)的更新，即：

$\hat{X}(k+1)=\hat{X}(k+1/k)+G(k+1)\hat{Z}(k+1)$

P(k+1)＝P(k+1/k)-G(k+1)S(k+1)G^T(k+1)

其中，表示机动目标更新后在第k+1采样点的状态向量；P(k+1)表示机动目标在第k+1采样点的协方差矩阵；[·]^T表示矩阵的转置；

(10)通过步骤(7)中得到的预测误差和步骤(5)中得到的目标运动模型的一步预测值对权值W(k)进行更新，得到更新后的权值：

$W(k+)=W\left(k\right)+2\mathrm{\μ}{\tilde{Z}}^T(k+1)H\left[\begin{array}{cccc}{\hat{X}}_1(k+1/k)& {\hat{X}}_2(k+1/k)& ...& {\hat{X}}_M(k+1/k)\end{array}\right]$

其中，W(k+1)表示第k+1采样点目标运动模型的权值，μ表示步长因子，M为采用的目标运动模型的个数；

(11)判断跟踪是否完成，即判断k是否满足k≤N-1，若条件满足，则k递增，并继续执行步骤(5)；若条件不满足，则跟踪过程结束。