CN105891820A - 基于ukf和iufir的机动目标跟踪方法 - Google Patents

Abstract

本发明公开了一种基于UKF和IUFIR的机动目标跟踪方法，它包括以下步骤：建立系统模型，设定目标的状态方程和测量方程；对目标初始状态和误差协方差进行初始化；使用UKF滤波方法进行线性转换；使用IUFIR滤波方法对线性化的数据进行滤波初始化，得到s时刻的广义噪声功率增益分别为和G_s；更新IUFIR滤波方法中的广义噪声功率增益G；对IUFIR滤波增益K更新；使用IUFIR滤波方法进行滤波估计，本发明解决在非线性雷达目标跟踪系统中，状态估计方面因为过程噪声的不确定性而增加预测估计的误差甚至发散的问题。

Description

基于UKF和IUFIR的机动目标跟踪方法

技术领域

本发明涉及空间机动目标跟踪技术领域，具体涉及基于UKF和IUFIR的机动目标跟踪方法。

背景技术

空间机动目标跟踪技术很早就被提出来了，但是它一直是雷达目标跟踪中的难点，随着数据关联等技术在多目标跟踪理论方面取得了开创性的突破，由数据关联技术和滤波技术相结合提出了很多新的目标跟踪算法。但是由于机动目标的机动情况时未知的，这就增加了状态方程的过程噪声难以确定，一般模型的过程噪声均为高斯白噪声，后来也有一些方法来解决这个问题，例如可调白噪声方法、变维滤波方法等，这些方法虽然通过各种方法来降低过程噪声的不确定性来提高跟踪的准确度，但是没有从根本上解决这个问题。

发明内容

本发明克服了现有技术的不足，提供一种基于UKF和IUFIR的机动目标跟踪方法。

为实现上述目的，本发明采用以下技术方案：

一种基于UKF和IUFIR的机动目标跟踪方法，它包括以下步骤：

步骤1、建立系统模型，设定目标的状态方程和测量方程；

步骤2、对目标初始状态和误差协方差进行初始化；

步骤3、使用UKF滤波方法进行线性转换；

步骤4、使用IUFIR滤波方法对线性化的数据进行滤波初始化，得到s时刻的广义噪声功率增益分别为和G_s；

步骤5、更新IUFIR滤波方法中的广义噪声功率增益G；

步骤6、对IUFIR滤波增益K更新；

步骤7、使用IUFIR滤波方法进行滤波估计。

更进一步的技术方案是，所述步骤1中建立的非线性跟踪系统模型中状态方程和测量方程分别为：

x_n＝f(x_n-1,θ_n-1)+w_n-1

z_n＝h_n(x_n,θ_n)+v_n

其中，x_n为目标状态向量，f为状态转移非线性函数，w_n-1为过程噪声，其协方差为Q_n，h为量测函数，v_n为测量噪声，z_n为系统测量向量。

更进一步的技术方案是，所述步骤3使用UKF滤波中的无迹转换进行线性转换得到离散的系统模型为：

x_n＝F_nx_n-1+B_nw_n

z_n＝H_nx_n+v_n

x_n∈R^K×1为n时刻目标状态向量，F_n为n时刻非线性状态函数，w_n∈R^p×1为n时刻过程噪声，v_n∈R^p×1为n时刻量测噪声，z_n∈R^m×1为n时刻系统的量测向量，H_n为n时刻非线性量测函数。

更进一步的技术方案是，所述步骤4使用IUFIR滤波方法对线性化的数据进行滤波初始化，得到s时刻的广义噪声功率增益分别为和G_s，得到：

${\hat{x}}_s={\mathrm{\Ξ}}_{s,0}^{m+1}{\mathrm{\Ψ}}_{s,m}^{-1}{Z}_{s,m}$

$G_s={\mathrm{\Ξ}}_{s,0}^{m+1}{\left({\mathrm{\Ψ}}_{s,m}^T{\mathrm{\Ψ}}_{s,m}\right)}^{-1}{\left({\mathrm{\Ξ}}_{s,0}^{m+1}\right)}^T$

其中，

${\mathrm{\Ξ}}_{s,0}^{m+1}=\underset{i=0}{\overset{s-(m+1)}{\Π}}{F}_{s-i}=F_s{F}_{s-1}...F_{m+1}$

${\mathrm{\Ψ}}_{s,m}^{-1}={\left(H_{s,m}^T{H}_{s,m}\right)}^{-1}{H}_{s,m}^T$

$Z_{s,m}=\left[\begin{array}{cccc}{z}_n^T& z_{n-1}^T& ...& z_m^T\end{array}\right].$

更进一步的技术方案是，所述滤波增益和状态估计，具体公式表示为：

$K_l=G_l{H_l}^T$

本技术方案中，无迹卡尔曼滤波方法中的无迹转换过程具体如下，无迹变换是对状态概率密度函数(Probability Density Function,PDF)做近似。它首先选择有限个近似高斯分布离散点(称为σ点)，对每个σ点施以非线性变换，得到一簇UT变换后的线性点，再使用UFIR滤波方法来对转换的线性离散系统进行线性滤波，实现目标的定位和跟踪。

(1)UT变换步骤为：

设为一个随机变量，g代表任意一个非线性函数，实现令y＝g(x),Cov(x，x)＝P_x。UT变换的具体过程如下：

先按照以下方法得到(2n_x+1)个带权值的粒子，m_i为粒子权值，取(i＝1…2n_x),选取原则如下：

$\left\{\begin{array}{cc}{\mathrm{\χ}}_0=\stackrel{\‾}{x}& i=0\\ x_i=\stackrel{\‾}{x}+{\left(\sqrt{(n_x+\mathrm{\λ})P_x}\right)}_i& i=1,2...,n_x\\ x_i=\stackrel{\‾}{x}+{\left(\sqrt{(n_x+\mathrm{\λ})P_x}\right)}_i& i=n_x+1,...,2n_x\end{array}\right.---\left(1\right)$

$\left\{\begin{array}{cc}{m}_0=\mathrm{\λ}/(n_x+\mathrm{\λ}),& i=0\\ m_i=0.5/(n_x+\mathrm{\λ}),& i=1,...,n_x\\ m_i=0.5/(n_x+\mathrm{\λ}),& i=n_x+1,...,2n_x\end{array}\right.---\left(2\right)$

其中，λ为尺度参数，用于调节δ点和先验均值的距离远近，一般取0或3-n_x；是指矩阵平方根的第i行或第i列。

这些δ点通过非线性函数得到，公式如下：

ξ_i＝g(χ_i) i＝0，1，…2n_x (3)

y的均值和方差可通过这些的线性化的δ点得出：

$\stackrel{\‾}{y}=\underset{i=0}{\overset{2n_x}{\Σ}}{m}_i{\mathrm{\ξ}}_i---\left(4\right)$

$P_y=\underset{i=0}{\overset{2n_x}{\Σ}}{m}_i({\mathrm{\ξ}}_i-\stackrel{\‾}{y}){({\mathrm{\ξ}}_i-\stackrel{\‾}{y})}^T---\left(5\right)$

(2)IUFIR滤波方法

IUFIR滤波方法主要针对线性滤波模型，假设离散线性模型系统的状态方程和观测方程如下：

x_n＝F_nx_n-1+B_nw_n (6)

z_n＝H_nx_n+v_n (7)

x_n∈R^K×1为n时刻目标状态向量，F_n为n时刻非线性状态函数，w_n∈R^p×1为n时刻过程噪声，v_n∈R^p×1为n时刻量测噪声，z_n∈R^m×1为n时刻系统的量测向量，H_n为n时刻非线性量测函数。

设在系统中的v_n、w_n为零均值噪声序列，它们互不相关，其协方差分别为：Q_n和R_n。对所有的i、j时刻，噪声的基本统计性质为：E{v_n}＝0，E{w_n}＝0，

当利用IUFIR滤波方法进行滤波时，对多个时刻抽样值取平均值进行滤波(即批处理IUFIR滤波)，此特性能够在忽视噪声统计特性的情况下对信号进行滤波。具体步骤如下：

初始化：通过批处理IUFIR滤波得到初始的s时刻的目标状态初值和广义噪声功率增益分别为和G_s，此时的批处理滤波器窗长为N。

${\hat{x}}_s={\mathrm{\Ξ}}_{s,0}^{m+1}{\mathrm{\Ψ}}_{s,m}^{-1}{Z}_{s,m}$

$G_s={\mathrm{\Ξ}}_{s,0}^{m+1}{\left({\mathrm{\Ψ}}_{s,m}^T{\mathrm{\Ψ}}_{s,m}\right)}^{-1}{\left({\mathrm{\Ξ}}_{s,0}^{m+1}\right)}^T---\left(8\right)$

其中：

${\mathrm{\Ξ}}_{s,0}^{m+1}=\underset{i=0}{\overset{s-(m+1)}{\Π}}{F}_{s-i}=F_s{F}_{s-1}...F_{m+1}---\left(9\right)$

${\mathrm{\Ψ}}_{s,m}^{-1}={\left(H_{s,m}^T{H}_{s,m}\right)}^{-1}{H}_{s,m}^T---\left(10\right)$

$Z_{s,m}=\left[\begin{array}{cccc}{z}_n^T& z_{n-1}^T& ...& z_m^T\end{array}\right]---\left(11\right)$

IUFIR滤波的迭代更新思想为：此时批处理的窗长发生变化为M，且M□N，其中s＝m+M-1，当迭代变量从m+M到n，当l＝n时，就得到n时刻的状态预测估计

${\hat{x}}_l=F{\hat{x}}_{l-1}+K_l(z_l-H_l{F}_l{\hat{x}}_{l-1})---\left(12\right)$

广义噪声功率增益的更新：在机动目标跟踪中，机动会导致量测数据和滤波的数据产生很大的偏差，此偏差能够目标反映机动的情况。但是广义噪声功率增益G只与状态转移矩阵和观测矩阵有关系，并且不同时刻的G之间是相互独立的，所以本发明IUFIR滤波方法通过调整系数γ来自适应调整G的值，提高滤波精度。

在i时刻，选取量测与滤波结果的偏差的平方根η_i来描述机动的基本信息，用公式表示为：

${\mathrm{\η}}_i=\sqrt{\frac{(z_i-H_i{\hat{x}}_i){(z_i-H_i{\hat{x}}_i)}^T}{\mathrm{\κ}}}---\left(13\right)$

其中，z_i为第i时刻的量测数据，为第i时刻的状态预测结果，H_i为第i时刻的状态观测矩阵，κ为目标运动的维数。

由于目标的机动情况只能通过不同时刻的偏差的比值才能观察出来，所以本发明将此时刻的偏差平方根与上一时刻的偏差的平方根的比值来反映机动情况，第i时刻与第i-1时刻的偏差的平方根比值用数学公式表示如下：

${\mathrm{\λ}}_i=\sqrt[\mathrm{\κ}]{{\mathrm{\η}}_i/{\mathrm{\η}}_{i-1}}---\left(14\right)$

同时考虑到野值可能会造成滤波发散的原因，以及一般机动会持续几秒到几十秒，所以G的调整系数γ取比值λ的均值的方法来消除野值影响，同时也不影响对目标机动的描述，在第i时刻得到G得调整系数γ为：

其中，由于必须满足l₀≥m+M+1，所以得到：说明对广义噪声功率增益的调整只能从时刻开始，第l时刻的广义噪声功率增益的计算表达式为：

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

IUFIR滤波方法的滤波增益为：

$K_l=G_l{H}_l^T---\left(17\right)$

IUFIR滤波方法的l时刻的状态预测估计为：

${\hat{x}}_l=F{\hat{x}}_{l-1}+K_l(z_l-H_l{F}_l{\hat{x}}_{l-1})---\left(18\right)$

当l＝n时，得到第n时刻的滤波状态估计为

与现有技术相比，本发明的有益效果是：

本发明解决在非线性雷达目标跟踪系统中，状态估计方面因为过程噪声的不确定性而增加预测估计的误差甚至发散的问题。

附图说明

图1为本发明一种实施例的基于UKF和IUFIR的机动目标跟踪方法的流程图。

图2为本发明一种实施例中跟踪路径和目标实际运行的路径的对比图。

具体实施方式

下面结合附图对本发明作进一步阐述。

如图1所示，一种基于UKF和IUFIR的机动目标跟踪方法，首选引入涉及本实施例的相关参数，并详细描述如下：

g代表任意一个非线性函数，y＝g(x)；

E(x)为x的期望值，Cov(x,x)＝P_x,P_x为协方差；

m_i为粒子权值，λ为尺度参数，和P_y为y的均值和方差：

在k-1时刻的线性化状态和均方差为P_k-1|k-1；

w^m、w^c为权值；

x_n∈R^K×1为n时刻目标状态向量；

F_n为n时刻非线性状态函数；

w_n∈R^p×1为n时刻过程噪声，v_n∈R^p×1为n时刻量测噪声；

z_n∈R^m×1为n时刻系统的量测向量；

H_n为n时刻非线性量测函数；

Q_n和R_n分别为噪声序列v_n、w_n的协方差分别；

IUFIR滤波方法中s时刻的目标状态初值和广义噪声功率增益初始值分别为和G_s；

批处理的窗长M，且M□N，其中s＝m+M-1；

n时刻的状态预测估计

在i时刻的量测与滤波结果的偏差的平方根η_i；

z_i为第i时刻的量测数据，为第i时刻的状态预测结果，H_i为第i时刻的状态观测矩阵；κ为目标运动的维数；

λ_i为第i时刻与第i-1时刻的偏差的平方根比值；

γ为广义噪声功率增益G的调整系数；

G_l为第l时刻的广义噪声功率增益；

K_l为IUFIR滤波方法的滤波增益；

为IUFIR滤波方法的l时刻的状态预测估计。

本发明模拟的环境参数为实验仿真参数：目标在雷达直角坐标系中初始位置为[10km,10km,10km]，初始速度为[100m/s,100m/s,100m/s]，加入随机噪声，仿真帧数为100，采样周期为T＝2s，蒙特卡罗采样次数位100。

步骤1、建立系统模型，设定本方法的目标的状态方程和测量方程，雷达测量系统为非线性的，本发明中假设非线性模型为：

x_n＝f(x_n-1,θ_n-1)+w_n-1 (19)

z_n＝h_n(x_n,θ_n)+v_n (20)

步骤2、对目标初始状态和误差协方差进行初始化，如图1的步骤101对无迹卡尔曼滤波初始化：k＝0，得到如下公式：

${\hat{x}}_{0|0}={\stackrel{\‾}{x}}_0=E\left(x_0\right)---\left(21\right)$

cov(x₀)＝P₀ (22)

步骤3、使用UKF滤波方法进行线性转换，如图1的步骤102、103，假设在k-1时刻的状态为p_k-1|k-1，首先计算在k-1时刻采样的样本点集σ点(i＝1，…2n)：

${\mathrm{\χ}}_{k-1|k-1}^0={\hat{x}}_{k-1|k-1}---\left(23\right)$

${\mathrm{\χ}}_{k-1|k-1}^i={\hat{x}}_{k-1|k-1}+{\left(\sqrt{(n+\mathrm{\λ})p_{k-1|k-1}}\right)}_i,i=1,...n---\left(24\right)$

${\mathrm{\χ}}_{k-1|k-1}^i={\hat{x}}_{k-1|k-1}-{\left(\sqrt{(n+\mathrm{\λ})p_{k-1|k-1}}\right)}_i,i=n+1,...2n---\left(25\right)$

通过计算样本点集σ点中的i＝1，2…，2n来计算p_k|k-1，计算如下：

${\mathrm{\χ}}_k^i=f_k({\mathrm{\χ}}_{k-1|k-1}^i,{\mathrm{\θ}}_{k-1|k-1}^i),i=1,2...,2n---\left(26\right)$

${\hat{x}}_{k|k-1}=\underset{i=0}{\overset{2n}{\Σ}}{w}^m{\mathrm{\χ}}_k^i---\left(27\right)$

$p_{k|k-1}=\underset{i=0}{\overset{2n}{\Σ}}{w}^c({\mathrm{\χ}}_k^i-{\hat{x}}_{k|k-1}){({\mathrm{\χ}}_k^i-{\hat{x}}_{k|k-1})}^T+Q_{k-1}---\left(28\right)$

其中w^m、w^c表达式如下：

w^m＝λ/(n+λ) (29)

w^c＝λ(n+λ)+(1-α²+β) (30)

等式(29)中的β表示状态的分布，对于高斯分布β＝2是最优的，如果状态变量是单变量，则最佳的选择是β＝0。

步骤4，使用改进的IUFIR滤波方法对线性化的数据进行滤波初始化，如图1的步骤104，得到s时刻的广义噪声功率增益分别为和G_s，经UKF的UT转换线性化的系统模型如式(6)和式(7)，此时的批处理滤波器窗长为N＝55。

${\hat{x}}_s={\mathrm{\Ξ}}_{s,0}^{m+1}{\mathrm{\Ψ}}_{s,m}^{-1}{Z}_{s,m}$

$G_s={\mathrm{\Ξ}}_{s,0}^{m+1}{\left({\mathrm{\Ψ}}_{s,m}^T{\mathrm{\Ψ}}_{s,m}\right)}^{-1}{\left({\mathrm{\Ξ}}_{s,0}^{m+1}\right)}^T---\left(31\right)$

$F=\left[\begin{array}{cccccc}1& 0& 0& T& 0& 0\\ 0& 1& 0& 0& T& 0\\ 0& 0& 1& 0& 0& T\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right],Q=\left[\begin{array}{cccccc}0.25T^4& 0& 0& 0.5T^3& 0& 0\\ 0& 0.25T^4& 0& 0& 0.5T^3& 0\\ 0& 0& 0.25T^4& 0& 0& 0.5T^3\\ 0.5T^3& 0& 0& T^2& 0& 0\\ 0& 0.5T^3& 0& 0& T^2& 0\\ 0& 0& 0.5T^3& 0& 0& T^2\end{array}\right],$

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

步骤5、更新IUFIR滤波方法中的广义噪声功率增益G，如图1的步骤105，此时批处理的窗长发生变化为M＝5，且M□N，其中s＝m-M-1，当迭代变量从m+M到n，当l＝n时，得到G_l为

步骤6、如图1的步骤106对IUFIR滤波增益K更新，得到：

步骤7、如图1的步骤107，使用IUFIR滤波方法进行滤波估计，得到l时刻的滤波估计为：

通过上述给定的参数，利用本发明去跟踪目标，虽然计算量增大了，但是滤波精度大大提高，大大降低了漏跟目标的概率，跟踪路径和目标实际运行的路径对比结果如图2所示，表明本发明能稳定精确的跟踪目标。

以上具体实施方式对本发明的实质进行详细说明，但并不能对本发明的保护范围进行限制，显而易见地，在本发明的启示下，本技术领域普通技术人员还可以进行许多改进和修饰，需要注意的是，这些改进和修饰都落在本发明的权利要求保护范围之内。

Claims (6)

1.一种基于UKF和IUFIR的机动目标跟踪方法，其特征在于，它包括以下步骤：

步骤1、建立系统模型，设定目标的状态方程和测量方程；

步骤2、对目标初始状态和误差协方差进行初始化；

步骤3、使用UKF滤波方法进行线性转换；

步骤4、使用IUFIR滤波方法对线性化的数据进行滤波初始化，得到s时刻的广义噪声功率增益分别为和G_s；

步骤5、更新IUFIR滤波方法中的广义噪声功率增益G；

步骤6、对IUFIR滤波增益K更新；

步骤7、使用IUFIR滤波方法进行滤波估计。

2.根据权利要求1所述的基于UKF和IUFIR的机动目标跟踪方法，其特征在于，所述步骤1中建立的非线性跟踪系统模型中状态方程和测量方程分别为：

x_n＝f(x_n-1,θ_n-1)+w_n-1

z_n＝h_n(x_n,θ_n)+v_n

其中，x_n为目标状态向量，f为状态转移非线性函数，w_n-1为过程噪声，其协方差为Q_n，h为量测函数，v_n为测量噪声，z_n为系统测量向量。

3.根据权利要求1所述的基于UKF和IUFIR的机动目标跟踪方法，其特征在于，所述步骤3使用UKF滤波中的无迹转换进行线性转换得到离散的系统模型为：

x_n＝F_nx_n-1+B_nw_n

z_n＝H_nx_n+v_n

x_n∈R^K×1为n时刻目标状态向量，F_n为n时刻非线性状态函数，w_n∈R^p×1为n时刻过程噪声，v_n∈R^p×1为n时刻量测噪声，z_n∈R^m×1为n时刻系统的量测向量，H_n为n时刻非线性量测函数。

4.根据权利要求1所述的基于UKF和IUFIR的机动目标跟踪方法，其特征在于，所述步骤4使用IUFIR滤波方法对线性化的数据进行滤波初始化，得到s时刻的广义噪声功率增益分别为和G_s，得到：

${\hat{x}}_s={\mathrm{\Ξ}}_{s,0}^{m+1}{\mathrm{\Ψ}}_{s,m}^{-1}{Z}_{s,m}$

$G_s={\mathrm{\Ξ}}_{s,0}^{m+1}{\left({\mathrm{\Ψ}}_{s,m}^T{\mathrm{\Ψ}}_{s,m}\right)}^{-1}{\left({\mathrm{\Ξ}}_{s,0}^{m+1}\right)}^T$

其中，

${\mathrm{\Ξ}}_{s,0}^{m+1}=\underset{i=0}{\overset{s-(m+1)}{\Π}}{F}_{s-i}=F_s{F}_{s-1}...F_{m+1}$

${\mathrm{\Ψ}}_{s,m}^{-1}={\left(H_{s,m}^T{H}_{s,m}\right)}^{-1}{H}_{s,m}^T$

$Z_{s,m}=\left[\begin{array}{cccc}{z}_n^T& z_{n-1}^T& ...& z_m^T\end{array}\right].$

5.根据权利要求1所述的基于UKF和IUFIR的机动目标跟踪方法，其特征在于，所述步骤5引进调整系数来更新IUFIR滤波方法中的广义噪声功率增益G：

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

6.根据权利要求1所述的基于UKF和IUFIR的机动目标跟踪方法，其特征在于，所述滤波增益和状态估计，具体公式表示为：

$K_l=G_l{H_l}^T$