一种基于无偏有限冲击响应滤波器(UFIR)的改进机动目标跟踪方法

摘要

本发明提供一种基于无偏有限冲击响应滤波器（UFIR）的改进机动目标跟踪方法，其实现步骤如下：（1）建立离散时变线性系统模型；（2）估计UFIR滤波器最优窗长等参数；（3）每一时刻利用K阶批处理UFIR滤波器生成迭代初始条件和初始广义噪声功率增益矩阵（GNPG）；（4）再利用改进两阶段式类卡尔曼UFIR滤波器得到该时刻的滤波结果。本发明利用量测与滤波结果之间的偏差能够部分反映目标机动这一事实，引入一个广义噪声功率增益调整系数对GNPG进行时时调整，实现了对目标机动的自适应，从而提高了类卡尔曼UFIR滤波器对于机动目标的跟踪效果。

说明书

技术领域

本发明属于目标跟踪领域，涉及一种基于无偏有限冲击响应滤波器（UFIR）的改进机动 目标跟踪方法。

背景技术

机动目标跟踪一直是雷达目标跟踪中的难点问题之一，原因就在于常规的卡尔曼滤波 中，目标的机动被理解为状态空间模型中状态方程的过程噪声。而由于目标的机动情况未知， 就造成状态方程的过程噪声难以确定进而影响了滤波的效果。针对卡尔曼滤波的这一缺点， 之后提出了很多改进算法，诸如可调白噪声算法（CN），变维滤波算法（VD），多模型（MM） 及交互式多模型（IMM）算法。这些算法大都通过设置多个过程噪声或者通过各种方法降低 过程噪声的不确定性以期获得比卡尔曼滤波更好的性能。但从本质上这些优化算法都没有从 根本上解决滤波对于过程噪声统计特性的依赖。

有限冲击响应滤波器（finite impulse response,FIR）在滤波过程中能够完全无视噪声的 统计特性，正好解决了过程噪声难以确定的问题。特别是类卡尔曼无偏有限冲击响应滤波器 （unbiased finite impulse response,UFIR）的提出大大降低了UFIR滤波器的计算负担，使其 在工程上能够得以应用。同时上述许多对卡尔曼滤波器在机动条件下的改进算法都能够被推 广到UFIR滤波器中。

发明内容

本发明要解决的技术问题是：在现有的UFIR滤波器中，随着系统状态空间模型的确定， 广义噪声功率增益（GNPG）也就随之确定。这就表明在整个滤波过程中，无论量测与滤波 结果存在怎样的偏差，新息都将被赋予相同的权重。这显然没有充分利用量测数据所有的信 息来降低对于航迹的不确定性。针对这一情况，提出了一种基于无偏有限冲击响应滤波器 （UFIR）的改进机动目标跟踪方法。该方法通过各个时刻量测与滤波结果之间的偏差来动态 调整GNPG，使算法具有了对于目标机动的自适应能力。

本发明解决其技术问题所采用的技术方案是：一种基于无偏有限冲击响应滤波器（UFIR） 的改进机动目标跟踪方法，实现步骤如下：首先建立离散时变线性系统模型；然后估计UFIR 滤波器最优窗长等参数；接着每一时刻利用K阶批处理UFIR滤波器生成迭代初始条件和初 始GNPG；最后再利用改进的两阶段式类卡尔曼UFIR滤波器得到每一时刻的滤波结果。具 体包括以下步骤：

步骤（1）、建立离散时变线性系统模型；

步骤（2）、估计UFIR滤波器最优窗长等参数；

步骤（3）、每一时刻利用K阶批处理UFIR滤波器生成迭代初始条件和初始广义噪声功 率增益矩阵（GNPG）；

步骤（4）、利用改进的两阶段式类卡尔曼UFIR滤波器得到该时刻的滤波结果。

所述步骤（1）建立离散时变线性系统模型如下：

离散时变线性系统模型用状态空间模型的状态和观测方程描述为：

x_n＝F_nx_n-1+B_nw_n，

z_n＝H_nx_n+v_n，

其中为第n时刻目标的状态，同理为第n-1时刻的状态，为系 统所得目标第n时刻的量测，为第n时刻系统的状态转移矩阵，为第n时 刻系统的量测矩阵，假定过程噪声矢量量测噪声矢量分布均为零均值高 斯白噪声；并假定两个噪声矢量互不相关，为第n时刻过程噪声与系统状态的关联 矩阵，同时定义x_n的p阶估计为表示在已知n时刻的量测时n+p时刻状态的估计值。

所述步骤（2）估计UFIR滤波器最优窗长等参数：

UFIR滤波器最优窗长通过最小化估计误差的均方值的所有状态或第k个状态分量求取， 所述的均方值计算方式如下：

$V_n=E\left\{(z_n-H{\hat{x}}_n){(z_n-H{\hat{x}}_n)}^T\right\}$

其中，为第n时刻滤波器对目标状态的估计值，为系统所得目标第n时 刻的量测，为系统的量测矩阵；

为了避免选取极小值时出现歧义，用于估计的点数应足够多并且提前平滑目标函数；

初始批处理窗长K并不会对滤波效果产生太大影响，为了降低计算量，一般所选取批处 理窗长远小于最优滤波窗长。

所述步骤（3）利用K阶批处理UFIR滤波器生成迭代初始条件和初始广义噪声功率增 益矩阵（GNPG），具体为：

迭代初始条件以及初始广义噪声功率增益矩阵G_s由K阶批处理UFIR滤波器计算得 到：

其中    

$Z_{s,m}={[z_s^T{z}_{s-1}^T···z_m^T]}^T,$

$H_{s,m}={\bar{H}}_{s,m}{F}_{s,m},$

${\bar{H}}_{s,m}=\mathrm{diag}(H_s{H}_{s-1}···H_m),$

式中，Z_s,m为由m时刻到s时刻的量测所构成的矩阵；为从m+1时刻到(s-0)时 刻系统状态转移矩阵的乘积，描述了系统从m时刻到s+1时刻的状态转移关系；描述 了将m时刻到s时刻的量测变换到m时刻状态的关系；为由m时刻到s时刻系统量测 矩阵构成的对角矩阵；F_n,m描述了将m时刻的状态变换到m时刻到s时刻的状态的关系； H_s,m描述了从m时刻到s时刻的状态变换到从m时刻到s时刻的量测的关系，而 s＝m+K-1，m为批处理UFIR滤波器的起始时刻，K为批处理窗长。

所述步骤（4）利用改进的两阶段式类卡尔曼UFIR滤波器得到滤波结果：

改进两阶段式类卡尔曼UFIR滤波器描述为：

${\hat{x}}_l=F_l{\hat{x}}_{l-1}+K_l(z_l-H_l{F}_l{\hat{x}}_{l-1}),$

其中      $K_l=G_l{H}_l^T,$

$G_l={\gamma }_l*{[H_l^T{H}_l+{\left(F_l{G}_{l-1}{F}_l^T\right)}^{-1}]}^{-1},$

式中，为估计第n时刻系统状态的迭代过程中所求得的第l时刻状态的估计值；同理 为估计第n时刻系统状态的迭代过程中所求得的第l-1时刻状态的估计值，F_l、H_l、z_l分别为第l时刻系统的状态转移矩阵、量测矩阵以及量测，其中迭代变量l从m+K变化到 n，当l＝n时，得到滤波器的迭代输出K_l表示第l时刻的新息增益矩阵，G_l表示第l 时刻的广义噪声功率增益矩阵，而γ_l即为广义噪声功率增益调整系数，通过下式计算得到：

其中描述为相邻两时刻偏差的均方根的比值， 表示取最优窗长一半的整数部分。而η_i描述为第i时 刻量测与滤波结果偏差的均方根，其中为机动目标运动的维数，Z_i为系统第i时刻所得量 测，H为系统的量测矩阵，为估计第n时刻系统状态的迭代过程中所求得的第i时刻状态 的估计值。

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

（1）、本发明将UFIR滤波器应用到机动目标跟踪中，解决了由目标机动的不确定而带 来的状态空间模型过程噪声难以确定的问题；

（2）、本发明由于不需要过程噪声的统计特性，表现出对于噪声变化更强的鲁棒性；

（3）、本发明利用量测与滤波结果的偏差来时时调整GNPG，实现对目标机动的自适应。

附图说明

图1为本发明基于UFIR滤波器的改进机动目标跟踪算法的流程图；

图2为机动目标仿真场景；

图3为最优条件下卡尔曼滤波器、原始UFIR及改进UFIR滤波器的机动目标跟踪性能； 其中，图3（a）卡尔曼最优情况下，UFIR滤波器与卡尔曼滤波器对比图（x方向）；图3（b） 卡尔曼最优情况下，UFIR滤波器与卡尔曼滤波器对比图（y方向）；图3（c）卡尔曼最优情 况下，UFIR滤波器与卡尔曼滤波器对比图（z方向）；

图4为变噪声条件下卡尔曼滤波器、原始UFIR及改进UFIR滤波器的机动目标跟踪性 能；其中，图4（a）噪声未知情况下，UFIR滤波器与卡尔曼滤波器对比图（x方向）；图4 （b）噪声未知情况下，UFIR滤波器与卡尔曼滤波器对比图（y方向）；图4（c）噪声未知 情况下，UFIR滤波器与卡尔曼滤波器对比图（z方向）。

具体实施方式

下面结合附图及具体实施方式详细介绍本发明。

本发明首先建立离散时变线性系统模型；然后估计UFIR滤波器最优窗长等参数；接着 每个时刻利用K阶批处理UFIR滤波器生成迭代初始条件和初始广义噪声功率增益矩阵 （GNPG）；最后再利用改进两阶段式类卡尔曼UFIR滤波器得到每一时刻的滤波结果。其实 施流程如图1所示，具体包含以下4个步骤：

1、建立离散时变线性系统模型

离散时变线性系统模型用状态空间模型的状态和观测方程描述为：

x_n＝F_nx_n-1+B_nw_n，

z_n＝H_nx_n+v_n，

其中为第n时刻目标的状态，同理x_n-1为第n-1时刻的状态，为系统所 得目标第n时刻的量测，为第n时刻系统的状态转移矩阵，为第n时刻系 统的量测矩阵。假定过程噪声矢量量测噪声矢量分布均为零均值高斯白 噪声，即E{w_n}＝0，E{v_n}＝0；并假定两个噪声矢量互不相关，即对于任意的i,j满足 E{w_iv_j^T}＝0。进而假定E[w_kw_j^T]＝Q(k)δ_kj，E[v_kv_j^T]＝R(k)δ_kj，Q(k)为过程噪声协方差 矩阵，R(k)为量测噪声协方差矩阵,其中δ_kj为Kronecker Delta函数。为第n时刻过 程噪声与系统状态的关联矩阵。同时定义x_n的p阶估计为表示在已知n时刻的量测时 n+p时刻状态的估计值。

2、估计UFIR滤波器最优窗长等参数

初始批处理窗长K并不会对滤波效果产生太大影响，但K的取值不易过大，因为这样将 造成UFIR滤波器计算量的急剧增加。特别当初始批处理窗长接近滤波窗长时，类卡尔曼 UFIR滤波器的计算量将接近批处理UFIR滤波器。所以一般批处理窗长只要求远远小于最 优滤波窗长即可。

UFIR滤波器最优窗长N_opt可以通过均方值（mean square value，MSV）：

$V_n=E\left\{(z_n-H{\hat{x}}_n){(z_n-H{\hat{x}}_n)}^T\right\}$

求取。其中为第n时刻对目标状态的估计值，为系统所得目标第n时刻的量 测，为系统的量测矩阵。

通过最小化V_n的所有状态或第k个状态分量的估计误差，能够求得N_opt的估计值或

${\hat{N}}_{\mathrm{opt}}\cong \underset{n}{\arg \min }\left(\frac{\partial }{\partial n}\mathrm{tr}{V}_n\right)+1,$

${\hat{N}}_{\mathrm{kopt}}\cong \underset{n}{\arg \min }\left(\frac{\partial }{\partial n}{V}_{\left(\mathrm{kk}\right)n}\right)+1,$

其中V_(kk)n是V_n的第k行k列的值。与的选取仅仅取决于具体的应用需要。从 K-1时刻开始直到或取得第一个极小值时的时刻，假设为时刻n，那n+1就是最优 窗长的估计值。为了避免选取极小值时出现歧义，必须保证用于估计的点数足够多并且提前 平滑目标函数。

3、利用K阶批处理UFIR滤波器生成迭代初始条件和初始广义噪声功率增益矩阵 （GNPG）

迭代初始条件以及初始广义噪声功率增益矩阵G_s由K阶批处理UFIR滤波器计算得 到：

其中      

$Z_{s,m}={[z_s^T{z}_{s-1}^T···z_m^T]}^T,$

$H_{s,m}={\bar{H}}_{s,m}{F}_{s,m},$

${\bar{H}}_{s,m}=\mathrm{diag}(H_s{H}_{s-1}···H_m),$

式中，Z_s,m为由m时刻到s时刻的量测所构成的矩阵；为从m+1时刻到(s-0)时 刻系统状态转移矩阵的乘积，描述了系统从m时刻到s+1时刻的状态转移关系；描述 了将m时刻到s时刻的量测变换到m时刻状态的关系；为由m时刻到s时刻系统量测 矩阵构成的对角矩阵；F_n,m描述了将m时刻的状态变换到m时刻到s时刻的状态的关系； H_s,m描述了从m时刻到s时刻的状态变换到从m时刻到s时刻的量测的关系，而 s＝m+K-1，m为批处理UFIR滤波器的起始时刻，K为批处理窗长。

4、利用改进的两阶段式类卡尔曼UFIR滤波器得到滤波结果

改进两阶段式类卡尔曼UFIR滤波器描述为：

${\hat{x}}_l=F_l{\hat{x}}_{l-1}+K_l(z_l-H_l{F}_l{\hat{x}}_{l-1}),$

其中       $K_l=G_l{H}_l^T,$

$G_l={\gamma }_l*{[H_l^T{H}_l+{\left(F_l{G}_{l-1}{F}_l^T\right)}^{-1}]}^{-1}$

式中，为估计第n时刻系统状态的迭代过程中所求得的第l时刻状态的估计值；同理 为估计第n时刻系统状态的迭代过程中所求得的第l-1时刻状态的估计值，F_l、H_l、z_l分别为第l时刻系统的状态转移矩阵、量测矩阵以及量测，其中迭代变量l从m+K变化到 n，当l＝n时，得到滤波器的迭代输出K_l表示第l时刻的新息增益矩阵，G_l表示第l 时刻的广义噪声功率增益矩阵，而γ_l即为广义噪声功率增益调整系数，通过下式计算得到：

其中描述为相邻两时刻偏差的均方根的比值， 表示取最优窗长一半的整数部分。而ηi描述为第i时 刻量测与滤波结果偏差的均方根，其中为机动目标运动的维数，Z_i为系统第i时刻所得量 测，H为系统的量测矩阵，为估计第n时刻系统状态的迭代过程中所求得的第i时刻状态 的估计值。

下面通过仿真对比改进的UFIR滤波器、卡尔曼滤波器以及原始UFIR滤波器的性能对 本发明进行验证。仿真场景如图2所示，其中时刻1到时刻20为加速区，时刻20到时刻40 为匀速区，时刻40到时刻70为转弯区，时刻70到时刻100为匀速区。仿真系统描述为用 状态空间表示的时不变系统，状态方程和量测方程都在直角坐标系下进行。状态x定义为xyz 方向上的位移、速度和加速度的九维向量，其中：

$F=\left(\begin{array}{ccccccccc}1& 0& 0& T& 0& 0& {0.5T}^2& 0& 0\\ 0& 1& 0& 0& T& 0& 0& {0.5T}^2& 0\\ 0& 0& 1& 0& 0& T& 0& 0& {0.5T}^2\\ 0& 0& 0& 1& 0& 0& T& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& T& 0\\ 0& 0& 0& 0& 0& 1& 0& 0& T\\ 0& 0& 0& 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1\end{array}\right),$

$H=\left(\begin{array}{ccccccccc}1& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0& 0& 0& 0\end{array}\right),$

B＝I,

$Q=\left(\begin{array}{ccccccccc}{0.25T}^4& 0& 0& {0.5T}^3& 0& 0& {0.5T}^2& 0& 0\\ 0& {0.25T}^4& 0& 0& {0.5T}^3& 0& 0& {0.5T}^2& 0\\ 0& 0& {0.25T}^4& 0& 0& {0.5T}^3& 0& 0& {0.5T}^2\\ {0.5T}^3& 0& 0& T^2& 0& 0& T& 0& 0\\ 0& {0.5T}^3& 0& 0& T^2& 0& 0& T& 0\\ 0& 0& {0.5T}^3& 0& 0& T^2& 0& 0& T\\ {0.5T}^2& 0& 0& T& 0& 0& 1& 0& 0\\ 0& {0.5T}^2& 0& 0& T& 0& 0& 1& 0\\ 0& 0& {0.5T}^2& 0& 0& T& 0& 0& 1\end{array}\right),$

$R=\left(\begin{array}{ccc}400& 0& 0\\ 0& 400& 0\\ 0& 0& 400\end{array}\right)$

其中T为采样间隔，仿真中T＝0.1。

首先当量测被附加一个协方差为100的零均值高斯白噪声时，卡尔曼滤波器近似达到最 优。通过计算求出此时UFIR滤波器的最优窗长N_opt＝55，并取批处理长度K=5。进行100 次蒙特卡洛实验得到滤波结果与真实航迹在转弯阶段的均方根误差如图3图所示。接着将附 加噪声的协方差调整为16，观察三种滤波器对于噪声变化的鲁棒性，仿真结果如图4所示。 可以看出本发明所提出的算法对机动目标的跟踪性能最佳。

本发明未详细阐述部分属于本领域技术人员的公知技术。