CN102338870B - 一种采用前向散射雷达的三维目标跟踪方法 - Google Patents

Abstract

本发明涉及一种采用前向散射雷达的三维目标跟踪方法，属于目标跟踪技术领域。跟踪方法包括三个步骤：进行***模型定义，包括前向散射雷达几何构型的定义，目标模型的定义和观测模型的定义；对目标的初始状态进行估计；将目标初始状态作为滤波初始值，采用扩展卡尔曼滤波算法进行递推滤波，实现前向散射雷达对三维运动目标的跟踪。本发明克服了一般观测噪声背景下高斯-牛顿迭代容易出现矩阵奇异计算不准、以及初值估计误差太大容易陷入局部极值的问题，能够在测量精度不高的情况下获得比较精确的滤波初始值，滤波收敛速度快且精度高、数据率高、计算量小。

Description

一种采用前向散射雷达的三维目标跟踪方法

技术领域

本发明涉及一种采用前向散射雷达的三维目标跟踪方法，属于目标跟踪技术领域。

背景技术

前向散射雷达利用双基地雷达截面积(RCS)的前向散射区(双基地角135°～180°)探测目标。在该区域，双基地RCS一般比单基地RCS大20～40dB，从而能够在前向散射区域实现对低速、小RCS目标(包括隐身目标)的有效检测和跟踪。由于其具有的反隐身和抗低空突防能力，利用前向散射雷达进行目标跟踪近年来受到了较多的关注。

研究表明，连续准谐波信号在前向散射目标探测中的应用最有前景。现有的跟踪方法和实验***均采用连续准谐波信号进行探测。在连续准谐波***中，三维目标的坐标估计是通过测量目标的多普勒频移、方位角和仰角获得的，它们与目标的运动参数(位置、速度)之间的关系是非线性的，因此前向散射雷达的目标跟踪问题可以归结为非线性滤波问题，包括目标的初始状态估计和跟踪保持算法。想要利用现有的跟踪方法(经典方法)获得目标的精确初始状态估计，就要求雷达***必须有非常高的角度测量精度和多普勒测量精度，然而在实际的***中，受限于天线尺寸和***复杂度，方位和俯仰角度的测量精度也非常有限，这种情况下经典方法得到的迭代初值离真实值较远，尤其是在基线方向的速度估计极为不准，误差可达基线方向速度的2～3倍量级，容易导致跟踪保持阶段滤波收敛速度很慢甚至滤波发散。由于经典方法采用高斯-牛顿迭代算法进行滤波，在跟踪保持阶段需要利用若干组观测矢量，经过多次迭代才能得到一个时刻的目标运动参数估计，计算量很大，数据率较低。

发明内容

本发明的目的是为了利用前向散射雷达实现对三维运动目标的初始状态估计和快速收敛滤波，同时提高数据率、减小计算量，提出了一种采用前向散射雷达的三维目标跟踪方法。

本发明的目的是通过以下技术方案实现的。

本发明的一种采用前向散射雷达的三维目标跟踪方法，其步骤如下：

1)进行***模型定义，包括前向散射雷达几何构型的定义，目标模型的定义和观测模型的定义，分别为：

1.1对前向散射雷达几何构型进行定义，设定x，y，z为目标的直角坐标，R_e为接收基地，T_r为发射基地，将发射基地和接收基地之间的线段称为基线，T_g为任意k时刻目标的位置，ψ为目标航迹与基线之间的夹角也即航迹角，θ，β分别为k时刻接收基地测得的目标方位角和仰角，h为目标高度，b为基线长度，设定目标的真实运动轨迹为直线AB，目标轨迹在水平面(xoy平面)上雷达探测范围之内的投影为线段CD；

1.2定义目标模型，设定目标沿着直线AB以航迹角ψ匀速穿过基线，三维跟踪的***状态转移方程为

X(k+1)＝ΦX(k)+Gv(k)                            (1)

设定目标在高度为h的水平面里作匀速直线运动，状态向量可以表示为X(k)＝[x_k y_k h V_x V_y]，x_k，y_k，h为k时刻目标的三维直角坐标，V_x，V_y为该时刻目标的x方向和y方向的速度分量；设定T为采样间隔，v(k)是过程噪声且为零均值高斯白噪声，则状态转移矩阵Φ和噪声分布矩阵G分别为

$\mathrm{\Φ}=\left[\begin{array}{ccccc}1& 0& 0& T& 0\\ 0& 1& 0& 0& T\\ 0& 0& 1& 0& 0\\ 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 1\end{array}\right]$ $G=\left[\begin{array}{ccc}\frac{T^2}{2}& 0& 0\\ 0& \frac{T^2}{2}& 0\\ 0& 0& \frac{T^2}{2}\\ T& 0& 0\\ 0& T& 0\end{array}\right]---\left(2\right)$

1.3定义观测模型，假设1～k时刻观测矢量可记为：

[(f_d1，θ₁，β₁)，L，(f_dk，θ_k，β_k)]^T

其中，f_dk，θ_k，β_k分别为k时刻多普勒频移、回波方位角和仰角的观测值，则观测方程为：

其中

为k时刻观测矢量(f_dk，θ_k，β_k)，为k时刻状态矢量[x_k y_k h V_x V_y]，观测噪声

为零均值高斯白噪声，噪声方差分别为

2)对目标的初始状态进行估计，具体步骤为：

2.1将雷达对目标取得的观测数据采用最小二乘法进行拟合，设定初始状态估计中采用[(f_d1，θ₁，β₁)，L，(f_dN，θ_N，β_N)]^T共N(N≥24)组观测值进行计算，将这N组观测值分为三组，分别为(f_d1，f_d2，L，f_dN)、(θ₁，θ₂，L θ_N)和(β₁，β₂，L，β_N)，并分别用最小二乘法拟合到一阶多项式，即对测量值进行平滑，得到平滑后的观测值(f_d1LS，f_d2LS，L，f_dNLS)，(θ_1LS，θ_2LS，L，θ_NLS)，(β_1LS，β_2LS，L，β_NLS)；

2.2取步骤2.1平滑后的观测值，即令

f_d1，f_d2，L f_dN＝f_d1LS，f_d2LS，L，f_dNLS，θ₁，θ₂，L，θ_N＝θ_1LS，θ_2LS，L，θ_NLS

β₁，β₂，L β_N＝β_1LS，β_2LS，L，β_NLS

进行初始状态估计，具体步骤为：

设定n(1＜n≤N)时刻目标的状态矢量表示为[x_n，y_n，h，V_x，V_y]，给出x_n，V_x，V_y三者之间的关系式

$\left[\begin{array}{cc}{a}_1& a_3\\ b_1& b_3\end{array}\right]\left[\begin{array}{c}{x}_k\\ V_x\end{array}\right]=\left[\begin{array}{c}{-a}_2\\ {-b}_2\end{array}\right]V_y---\left(5\right)$

其中系数a₁，a₂，a₃，b₁，b₂，b₃的确定方法为

$\left\{\begin{array}{cc}{a}_1=({\tan \mathrm{\θ}}_n-{\tan \mathrm{\θ}}_1)& b_1=({\tan \mathrm{\θ}}_n-{\tan \mathrm{\θ}}_2)\\ a_2=-(n-1)T& b_2=-(n-2)T\\ a_3={\tan \mathrm{\θ}}_1(n-1)T& b_3={\tan \mathrm{\θ}}_2(n-2)T\end{array}\right.---\left(6\right)$

解方程组(5)，得到x_n，V_x，y_n与V_y之间的关系表达式分别为

$x_n=\frac{a_3{b}_2-a_2{b}_3}{a_1{b}_3-a_3{b}_1}\·V_y=d_{n1}{V}_y---\left(7\right)$

$V_x=\frac{a_2{b}_1-a_1{b}_2}{a_1{b}_3-a_3{b}_1}\·V_y=d_{n2}{V}_y---\left(8\right)$

y_n＝tanθ_kgd_n1V_y＝d_n3V_y          (9)

给出n时刻目标状态分量h，x_n，y_n之间的关系式

$h={\tan \mathrm{\β}}_n\sqrt{{x_n}^2+{y_n}^2}---\left(10\right)$

将式(7)和(9)代入(10)可得h与V_y之间的关系式

$h={\tan \mathrm{\β}}_n\sqrt{d_{n1}^2+d_{n3}^2}\left|V_y\right|=d_{n4}\left|V_y\right|---\left(11\right)$

给出f_dn与V_y之间的关系式

$f_{\mathrm{dn}}+\frac{1}{\mathrm{\λ}}[\frac{(d_{n1}{d}_{n2}+d_{n3})V_y^2}{\sqrt{(d_{n1}^2+d_{n3}^2+d_{n4}^2)V_y^2}}+\frac{d_{n3}{V}_y^2-(b-d_{n1}{V}_y)d_{n2}{V}_y}{\sqrt{{(b-d_{n1}{V}_y)}^2+d_{n3}^2{V}_y^2+d_{n4}^2{V}_y^2}}]=0---\left(12\right)$

将d_n1，d_n2，d_n3，d_n4代入(12)，得到关于V_y的一元非线性方程，将该一元非线性方程中变量n依次取n＝(N/2+1):N，计算每个n对应的d_n1，d_n2，d_n3，d_n4，得到如下非线性方程组

$\left\{\begin{array}{c}{f}_{d2}+\frac{1}{\mathrm{\λ}}[\frac{(d_{21}{d}_{22}+d_{23})V_y^2}{\sqrt{(d_{21}^2+d_{23}^2+d_{24}^2)V_y^2}}+\frac{d_{23}{V}_y^2-(b-d_{21}{V}_y)d_{22}{V}_y}{\sqrt{{(b-d_{21}{V}_y)}^2+d_{23}^2{V}_y^2+d_{24}^2{V}_y^2}}]=0\\ f_{d3}+\frac{1}{\mathrm{\λ}}[\frac{(d_{31}{d}_{32}+d_{33})V_y^2}{\sqrt{(d_{31}^2+d_{33}^2+d_{34}^2)V_y^2}}+\frac{d_{33}{V}_y^2-(b-d_{31}{V}_y)d_{32}{V}_y}{\sqrt{{(b-d_{31}{V}_y)}^2+d_{33}^2{V}_y^2+d_{34}^2{V}_y^2}}]=0\\ M\\ f_{\mathrm{dN}}+\frac{1}{\mathrm{\λ}}[\frac{(d_{N1}+d_{N2}+d_{N3})V_y^2}{\sqrt{(d_{N1}^2+d_{N3}^2+d_{N4}^2)V_y^2}}+\frac{d_{N3}{V}_y^2-(b-d_{N1}{V}_y)d_{N3}{V}_y}{\sqrt{{(b-d_{N1}{V}_y)}^2+d_{N3}^2{V}_y^2+d_{N4}^2{V}_y^2}}]=0\end{array}\right.---\left(13\right)$

采用Levenberg-Marquardt算法解非线性方程组(13)，得到V_y的最优解，根据该V_y的最优解，基于公式(7)、(8)、(9)、(11)，求解得到x_N，V_x，y_N，h的状态值，也即得到N时刻目标状态(x_N，y_N，h，V_x，V_y)的数值解；

3)将步骤2)得到的目标初始状态也即N时刻目标状态(x_N，y_N，h，V_x，V_y)的数值解作为滤波初始值，采用扩展卡尔曼滤波(EKF)算法进行递推滤波，得到N及N以后各个时刻目标的状态估计，即实现了前向散射雷达对三维运动目标的跟踪。

有益效果

本发明与现有技术相比，其优势在于：首先，在目标初始状态估计中，基于解析推导获得非线性方程组，并采用Levenberg-Marquardt算法解方程组，克服了一般观测噪声背景下高斯-牛顿迭代容易出现矩阵奇异计算不准、以及初值估计误差太大容易陷入局部极值的问题，能够在测量精度不高的情况下获得比较精确的滤波初始值；其次，利用得到的精确初始状态估计值，与扩展卡尔曼滤波(EKF)算法相结合进行递推滤波，滤波收敛速度快且精度高、数据率高、计算量小。

附图说明

图1为本发明中三维前向散射雷达的几何结构示意图；

图2为本发明实施例中采用的目标初始状态估计方法与经典方法估计的初值在相同条件下进行EKF滤波的位置估计均方误差对比示意图；

图3为本发明实施例中采用的目标初始状态估计方法与经典方法估计的初值在相同条件下进行EKF滤波的速度估计均方误差对比示意图；

图4为本发明实施例中采用的三维目标跟踪方法与经典方法、不敏卡尔曼滤波(UKF)算法进行目标跟踪的仿真对比效果图；

图5为本发明实施例中采用的三维目标跟踪方法与经典方法的基线方向位置估计均方误差对比示意图；

图6为本发明实施例中采用的三维目标跟踪方法与经典方法的基线方向速度V_x估计均方误差对比示意图；

图7为本发明实施例中采用的三维目标跟踪方法与经典方法的垂直基线方向速度V_y估计均方误差对比示意图。

具体实施方式

下面结合附图和实施例对本发明做进一步说明。

实施例

采用俄罗斯A.B.Blyakhman等人研制的双基地前向散射雷达实验***中的参数进行三维目标跟踪仿真。取工作波长为0.77m，基线长度为40km，数据率为1Hz，多普勒频移标准差为1，方位角测量标准差为0.5°，仰角测量标准差为0.5°，过程噪声标准差为1，目标速度为150m/s，航迹角为15°，目标初始位置的三维直角坐标为(20，-6，2)km。

一种采用前向散射雷达的三维目标跟踪方法，其步骤如下：

本发明的一种采用前向散射雷达的三维目标跟踪方法，其步骤如下：

1)进行***模型定义，包括前向散射雷达几何构型的定义，目标模型的定义和观测模型的定义，分别为：

1.1对前向散射雷达几何构型进行定义，如图1所示，设定x，y，z为目标的直角坐标，R_e为接收基地，T_r为发射基地，将发射基地和接收基地之间的线段称为基线，T_g为任意k时刻目标的位置，ψ为目标航迹与基线之间的夹角，也即航迹角，θ，β分别为k时刻接收基地测得的目标方位角和仰角，h为目标高度，b为基线长度，设定目标的真实运动轨迹为直线AB，目标轨迹在水平面(xoy平面)上雷达探测范围之内的投影为线段CD；

1.2定义目标模型，设定目标沿着直线AB以航迹角ψ匀速穿过基线，三维跟踪的***状态转移方程为

X(k+1)＝ΦX(k)+Gv(k)                     (14)

设定目标在高度为h的水平面里作匀速直线运动，状态向量可以表示为X(k)＝[x_k y_k h V_x V_y]，x_k，y_k，h为k时刻目标的三维直角坐标，V_x，V_y为该时刻目标的x方向和y方向的速度分量；设定T为采样间隔，v(k)是过程噪声，为零均值高斯白噪声，则状态转移矩阵Φ和噪声分布矩阵G分别为

$\mathrm{\Φ}=\left[\begin{array}{ccccc}1& 0& 0& T& 0\\ 0& 1& 0& 0& T\\ 0& 0& 1& 0& 0\\ 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 1\end{array}\right]$ $G=\left[\begin{array}{ccc}\frac{T^2}{2}& 0& 0\\ 0& \frac{T^2}{2}& 0\\ 0& 0& \frac{T^2}{2}\\ T& 0& 0\\ 0& T& 0\end{array}\right]---\left(15\right)$

1.3定义观测模型，假设1:k时刻观测矢量可记为：

[(f_d1，θ₁，β₁)，L，(f_dk，θ_k，β_k)]^T

其中，f_dk，θ_k，β_k分别为k时刻多普勒频移、回波方位角和仰角的观测值，则观测方程为：

其中

为k时刻观测矢量(f_dk，θ_k，β_k)，

为k时刻状态矢量[x_k y_k h V_x V_y]，观测噪声

为零均值高斯白噪声，噪声方差分别为

2)对目标的初始状态进行估计，具体步骤为：

2.1将雷达对目标取得的观测数据采用最小二乘法进行拟合，设定初始状态估计中采用[(f_d1，θ₁，β₁)，L，(f_dN，θ_N，β_N)]^T共N(N≥24)组观测值进行计算，将这N组观测值分为三组，分别为(f_d1，f_d2，L，f_dN)、(θ₁，θ₂，L θ_N)和(β₁，β₂，L，β_N)，并分别用最小二乘法拟合到一阶多项式，即对测量值进行平滑，得到平滑后的观测值(f_d1LS，f_d2LS，L，f_dNLS)，(θ_1LS，θ_2LS，L，θ_NLS)，(β_1LS，β_2LS，L，β_NLS)；

2.2取步骤2.1平滑后的观测值，即令

f_d1，f_d2，L f_dN＝f_d1LS，f_d2LS，L，f_dNLS，θ₁，θ₂，L，θ_N＝θ_1LS，θ_2LS，L，θ_NLS

β₁，β₂，L β_N＝β_1LS，β_2LS，L，β_NLS

进行初始状态估计，具体步骤为：

设定n(1＜n≤N)时刻目标的状态矢量表示为[x_n，y_n，h，V_x，V_y]，由于x_n，V_x，V_y满足下面关系

$\left[\begin{array}{cc}{a}_1& a_3\\ b_1& b_3\end{array}\right]\left[\begin{array}{c}{x}_k\\ V_x\end{array}\right]=\left[\begin{array}{c}{-a}_2\\ {-b}_2\end{array}\right]V_y---\left(18\right)$

其中系数a₁，a₂，a₃，b₁，b₂，b₃的确定方法为将平滑后的前n组观测值

代入下式

$\left\{\begin{array}{cc}{a}_1=({\tan \mathrm{\θ}}_n-{\tan \mathrm{\θ}}_1)& b_1=({\tan \mathrm{\θ}}_n-{\tan \mathrm{\θ}}_2)\\ a_2=-(n-1)T& b_2=-(n-2)T\\ a_3={\tan \mathrm{\θ}}_1(n-1)T& b_3={\tan \mathrm{\θ}}_2(n-2)T\end{array}\right.---\left(19\right)$

解方程组(5)，得到x_n，V_x，y_n与V_y之间的关系表达式为

$x_n=\frac{a_3{b}_2-a_2{b}_3}{a_1{b}_3-a_3{b}_1}\·V_y=d_{n1}{V}_y---\left(20\right)$

$V_x=\frac{a_2{b}_1-a_1{b}_2}{a_1{b}_3-a_3{b}_1}\·V_y=d_{n2}{V}_y---\left(21\right)$

y_n＝tanθ_kgd_n1V_y＝d_n3V_y            (22)

根据n时刻目标状态分量h，x_n，y_n满足下面关系

$h={\tan \mathrm{\β}}_n\sqrt{{x_n}^2+{y_n}^2}---\left(23\right)$

将式(7)和(9)代入(10)可得h与V_y之间的关系式

$h={\tan \mathrm{\β}}_n\sqrt{d_{n1}^2+d_{n3}^2}\left|V_y\right|=d_{n4}\left|V_y\right|---\left(24\right)$

根据f_dn与V_y满足下面关系

$f_{\mathrm{dn}}+\frac{1}{\mathrm{\λ}}[\frac{(d_{n1}{d}_{n2}+d_{n3})V_y^2}{\sqrt{(d_{n1}^2+d_{n3}^2+d_{n4}^2)V_y^2}}+\frac{d_{n3}{V}_y^2-(b-d_{n1}{V}_y)d_{n2}{V}_y}{\sqrt{{(b-d_{n1}{V}_y)}^2+d_{n3}^2{V}_y^2+d_{n4}^2{V}_y^2}}]=0---\left(25\right)$

将d_n1，d_n2，d_n3，d_n4代入(12)，得到关于V_y的一元非线性方程，将该一元非线性方程中变量n依次取n＝(N/2+1):N，计算每个n对应的d_n1，d_n2，d_n3，d_n4，得到如下非线性方程组

$\left\{\begin{array}{c}{f}_{d2}+\frac{1}{\mathrm{\λ}}[\frac{(d_{21}{d}_{22}+d_{23})V_y^2}{\sqrt{(d_{21}^2+d_{23}^2+d_{24}^2)V_y^2}}+\frac{d_{23}{V}_y^2-(b-d_{21}{V}_y)d_{22}{V}_y}{\sqrt{{(b-d_{21}{V}_y)}^2+d_{23}^2{V}_y^2+d_{24}^2{V}_y^2}}]=0\\ f_{d3}+\frac{1}{\mathrm{\λ}}[\frac{(d_{31}{d}_{32}+d_{33})V_y^2}{\sqrt{(d_{31}^2+d_{33}^2+d_{34}^2)V_y^2}}+\frac{d_{33}{V}_y^2-(b-d_{31}{V}_y)d_{32}{V}_y}{\sqrt{{(b-d_{31}{V}_y)}^2+d_{33}^2{V}_y^2+d_{34}^2{V}_y^2}}]=0\\ M\\ f_{\mathrm{dN}}+\frac{1}{\mathrm{\λ}}[\frac{(d_{N1}+d_{N2}+d_{N3})V_y^2}{\sqrt{(d_{N1}^2+d_{N3}^2+d_{N4}^2)V_y^2}}+\frac{d_{N3}{V}_y^2-(b-d_{N1}{V}_y)d_{N3}{V}_y}{\sqrt{{(b-d_{N1}{V}_y)}^2+d_{N3}^2{V}_y^2+d_{N4}^2{V}_y^2}}]=0\end{array}\right.---\left(26\right)$

采用Levenberg-Marquardt算法解非线性方程组(13)，得到V_y的最优解；

根据上述V_y的最优解，基于公式(7)、(8)、(9)、(11)，求解得到x_N，V_x，V_y，h的状态值，也即得到N时刻目标状态(x_N，y_N，h，V_x，V_y)的数值解；

步骤2.2中x_n，V_x，V_y满足的关系，h，x_n，y_n满足的关系，以及f_dn与V_y满足的关系推导过程如下：

根据***观测方程，对于n时刻的观测值(f_dn，θ_n，β_n)，有

$f_{\mathrm{dn}}=-\frac{1}{\mathrm{\λ}}[\frac{x_n{V}_x+y_n{V}_y}{\sqrt{x_n^2+y_n^2+h^2}}+\frac{y_n{V}_y-(b-x_n)V_x}{\sqrt{{(b-x_n)}^2+y_n^2+h^2}}]---\left(27\right)$

${\mathrm{\θ}}_n=\arctan \frac{y_n}{x_n}---\left(28\right)$

${\mathrm{\β}}_n=\arctan \frac{h}{\sqrt{{x_n}^2+{y_n}^2}}---\left(29\right)$

由方程(28)可得

$\frac{y_1}{x_1}={\tan \mathrm{\θ}}_{1}L\frac{y_n}{x_n}={\tan \mathrm{\θ}}_n---\left(30\right)$

由目标运动模型可知y₁与y_n之间存在关系

$\left\{\begin{array}{c}{y}_1=y_n-(n-1){\mathrm{TV}}_y\\ x_1=x_n-(n-1){\mathrm{TV}}_x\end{array}\right.---\left(31\right)$

将公式(31)代入(30)可得到

$\frac{y_n-(k-1){\mathrm{TV}}_y}{x_n-(k-1){\mathrm{TV}}_x}={\tan \mathrm{\θ}}_1---\left(32\right)$

同时结合(30)和(32)，有

(tanθ_n-tanθ₁)x_n-(n-1)TV_y+tanθ₁(n-1)TV_x＝0(33)

同理可得

(tanθ_n-tanθ₂)x_n-(n-2)TV_y+tanθ₂(n-2)TV_x＝0(34)

联立(33)和(34)可列出如下方程组

令方程组(35)中x_n，V_y，V_x的系数分别为(a₁，a₂，a₃)，(b₁，b₂，b₃)，即

$\left\{\begin{array}{cc}{a}_1=({\tan \mathrm{\θ}}_n-{\tan \mathrm{\θ}}_1)& b_1=({\tan \mathrm{\θ}}_n-{\tan \mathrm{\θ}}_2)\\ a_2=-(n-1)T& b_2=-(n-2)T\\ a_3={\tan \mathrm{\θ}}_1(n-1)T& b_3={\tan \mathrm{\θ}}_2(n-2)T\end{array}\right.---\left(36\right)$

则方程组(35)可写成如下形式

$\left[\begin{array}{cc}{a}_1& a_3\\ b_1& b_3\end{array}\right]\left[\begin{array}{c}{x}_k\\ V_x\end{array}\right]=\left[\begin{array}{c}{-a}_2\\ {-b}_2\end{array}\right]V_y---\left(37\right)$

即x_n，V_x，V_y满足的关系；

解方程(37)可得

$x_n=\frac{a_3{b}_2-a_2{b}_3}{a_1{b}_3-a_3{b}_1}\·V_y=d_{n1}{V}_y---\left(38\right)$

$V_x=\frac{a_2{b}_1-a_1{b}_2}{a_1{b}_3-a_3{b}_1}\·V_y=d_{n2}{V}_y---\left(39\right)$

y_n＝tanθ_kgd_n1V_y＝d_n3V_y      (40)

根据方程(29)

$h={\tan \mathrm{\β}}_{n}g\sqrt{{x_n}^2+{y_n}^2}---\left(41\right)$

即h，x_n，y_n满足的关系；

将(38)，(40)代入(41)可得

$h={\tan \mathrm{\β}}_n\sqrt{d_{n1}^2+d_{n3}^2}\left|V_y\right|=d_{n4}\left|V_y\right|---\left(42\right)$

将式(38)，(39)，(40)，(42)代入方程(27)，得到关于V_y的一元非线性方程

$f_{\mathrm{dn}}+\frac{1}{\mathrm{\λ}}[\frac{(d_{n1}{d}_{n2}+d_{n3})V_y^2}{\sqrt{(d_{n1}^2+d_{n3}^2+d_{n4}^2)V_y^2}}+\frac{d_{n3}{V}_y^2-(b-d_{n1}{V}_y)d_{n2}{V}_y}{\sqrt{{(b-d_{n1}{V}_y)}^2+d_{n3}^2{V}_y^2+d_{n4}^2{V}_y^2}}]=0---\left(43\right)$

即f_dn与V_y的关系；

3)将步骤2)得到的目标初始状态也即N时刻目标状态(x_N，y_N，h，V_x，V_y)的数值解作为滤波初始值采用扩展卡尔曼滤波(EKF)算法进行递推滤波，得到≥N的各个时刻目标的状态估计，即实现了前向散射雷达对三维运动目标的跟踪。

仿真对比实验目标初始状态估计精度对滤波收敛速度和精度的影响如图2、图3所示，分别采用经典方法和本发明的方法估计的目标初始状态作为初始值进行EKF滤波，完成100次蒙特卡洛仿真并取平均；图2中，实线是采用经典方法得到的初始状态估计值作为滤波初值，然后进行EKF滤波的X(基线)方向位置估计均方误差曲线；虚线是采用新方法得到的初始状态估计值作为滤波初值，进行EKF滤波的X(基线)方向位置估计均方误差曲线；由两条曲线对比可知，本文提出的初值估计方法在相同测量精度下可以得到比较精确的滤波初值，从而使EKF算法快速收敛，跟踪精度也比较高；而经典方法估计的初始值使得EKF滤波的收敛速度很慢，在目标穿越基线时刻还没有收敛到最终的跟踪精度范围内；图3图例中的“新方法基线方向速度”和“新方法垂直基线方向速度”分别表示采用本发明的方法得到的初值进行EKF滤波的速度V_x，V_y估计误差曲线，“经典方法基线方向速度”和“经典方法垂直基线方向速度”分别表示采用经典方法得到的初值进行EKF滤波的速度V_x，V_y估计误差曲线，对比图3中四条曲线可以发现，采用经典方法估计的初始状态作为EKF滤波初值时，X方向的速度估计误差很大，导致滤波收敛慢、精度低；本文提出的新方法估计的初值代入EKF滤波时，V_x和V_y速度估计误差都很小，因此滤波收敛快，位置估计精度比较高；仿真对比实验不同跟踪方法的性能比较如图4～图7所示，在观测噪声为(0.5Hz，0.5°，0.5°)时，经典方法得到的初始状态估计误差比较大，而且在后面的滤波过程中并没有明显的减小初始估计误差，滤波性能不稳定，滤波精度较差；本文发明的跟踪方法得到的初始状态估计误差较小(尤其是基线方向的位置和速度)，在后面的滤波过程中可以有效的减小误差并收敛，UKF算法滤波精度也比较高，但是性能不及本发明的方法，且滤波误差有增大的趋势；仿真实验对比不同跟踪方法的计算量如下，进行100次蒙特卡洛仿真并取平均，三种方法完成一次目标跟踪的平均运行时间分别为：经典方法需要8.564927s，本发明的方法需要0.14912884s，UKF算法需要0.32414323s。

上述结果说明，在前向散射雷达的目标跟踪中采用本发明提供的技术方案，对于测量误差的容忍度较大，在一般测量精度下，该方法得到的初始状态估计误差比经典方法的估计误差小一个数量级；采用该方法估计的初始状态作为滤波初值进行扩展卡尔曼滤波(EKF)，可以有效提高滤波算法的收敛速度和估计精度，从而能够在目标穿越基线之前形成稳定的航迹估计，计算量小，数据率高；可以满足前向散射雷达对三维匀速运动目标进行跟踪的应用要求。

以上所述为本发明的较佳实施例而已，本发明不应该局限于该实施例和附图所公开的内容。凡是不脱离本发明所公开的精神下完成的等效或修改，都落入本发明保护的范围。

Claims (1)

1.一种采用前向散射雷达的三维目标跟踪方法，其特征在于步骤如下：

1)进行***模型定义，包括前向散射雷达几何构型的定义，目标模型的定义和观测模型的定义，分别为：

1.1对前向散射雷达几何构型进行定义，设定x，y，z为目标的直角坐标，R_e为接收基地，T_r为发射基地，将发射基地和接收基地之间的线段称为基线，T_g为任意k时刻目标的位置，ψ为目标航迹与基线之间的夹角也即航迹角，θ，β分别为k时刻接收基地测得的目标方位角和仰角，h为目标高度，b为基线长度，设定目标的真实运动轨迹为直线AB，目标轨迹在水平面也即xoy平面上雷达探测范围之内的投影为线段CD；

1.2定义目标模型，设定目标沿着直线AB以航迹角ψ匀速穿过基线，三维跟踪的***状态转移方程为

X(k+1)＝ΦX(k)+Gv(k)                     (1)

设定目标在高度为h的水平面里作匀速直线运动，状态向量可以表示为X(k)＝[x_k y_k h V_x V_y]，x_k，y_k，h为k时刻目标的三维直角坐标，V_x，V_y为该时刻目标的x方向和y方向的速度分量；设定T为采样间隔，v(k)是过程噪声且为零均值高斯白噪声，则状态转移矩阵Ф和噪声分布矩阵G分别为 

，

1.3定义观测模型，假设1～k时刻观测矢量可记为：

[(f_d1，θ₁，β₁)，...，(f_dk，θ_k，β_k)]^T

其中，f_dk，θ_k，β_k分别为k时刻多普勒频移、回波方位角和仰角的观测值，则观测方程为：

其中

为k时刻观测矢量(f_dk，θ_k，β_k)， 为k时刻状态矢量[x_k y_k h V_x V_y]，观测噪声 

为零均值高斯白噪声，噪声方差分别为 

2)对目标的初始状态进行估计，具体步骤为：

2.1将雷达对目标取得的观测数据采用最小二乘法进行拟合，设定初始状态估计中采用[(f_d1，θ₁，β₁)，...，(f_dN，θ_N，β_N)]^T共N组观测值进行 计算，将这N组观测值分为三组，分别为(f_d1，f_d2，...，f_dN)、(θ₁，θ₂，...θ_N)和(β₁，β₂，...，β_N)，并分别用最小二乘法拟合到一阶多项式，即对测量值进行平滑，得到平滑后的观测值(f_d1LS，f_d2LS，...，f_dNLS)，(θ_1LS，θ_2LS，...，θ_NLS)，(β_1LS，β_2LS，...，β_NLS)；

2.2取步骤2.1平滑后的观测值，即令

f_d1，f_d2，... f_dN＝f_d1LS，f_d2LS，...，f_dNLS，θ₁，θ₂，...，θ_N＝θ_1LS，θ_2LS，...，θ_NLS，β₁，β₂，... β_N＝β_1LS，β_2LS，...，β_NLS

进行初始状态估计，具体步骤为：

设定n(1＜n≤N)时刻目标的状态矢量表示为[x_n，y_n，h，V_x，V_y]，给出x_n，V_x，V_y三者关系式

其中系数a₁，a₂，a₃，b₁，b₂，b₃的确定方法为

解方程组(5)，得到x_n，V_x，y_n与V_y之间的关系表达式分别为

y_n＝tanθ_k·d_n1V_y＝d_n3V_y        (9)

给出n时刻目标状态分量h，x_n，y_n的关系式

将式(7)和(9)代入(10)可得h与V_y之间的关系式 

给出f_dn与V_y之间的关系式

将d_n1，d_n2，d_n3，d_n4代入(12)，得到关于V_y的一元非线性方程，将该一元非线性方程中变量n依次取n＝(N/2+1):N，计算每个n对应的d_n1，d_n2，d_n3，d_n4，得到如下非线性方程组

采用Levenberg-Marquardt算法解非线性方程组(13)，得到V_y的最优解，根据该V_y的最优解，基于公式(7)、(8)、(9)、(11)，求解得到x_N，V_x，y_N，h的状态值，也即得到N时刻目标状态[ x_N，y_N，h，V_x，V_y ]的数值解；

3)将步骤2)得到的目标初始状态也即N时刻目标状态[ x_N，y_N，h，V_x，V_y ]的数值解作为滤波初始值，采用扩展卡尔曼滤波算法进行递推滤波，得到N及N以后各个时刻目标的状态估计，即实现了前向散射雷达对三维运动目标的跟踪。 

Images