CN102508238B - 一种基于坐标旋转变换的雷达跟踪方法 - Google Patents

Abstract

本发明公开了一种基于坐标旋转变换的雷达跟踪方法，该方法基于随机变量相关系数和坐标旋转变换的原理，以达到在定量度量雷达滤波系统模型中的量测方程的非线性程度的基础上，降低雷达滤波系统模型中的量测方程非线性程度，从而提高雷达跟踪效果的目的；具体步骤为：定义二维雷达量测极-直坐标转换线性度ρ；当r₀、σ_r和σ_a为确定的任意数值时，通过坐标旋转变换

使得ρ在

时取得最大值；设置目标运动状态，进行滤波初始化；分析由k时刻到k+1时刻的目标运动状态并建立坐标旋转变换滤波模型下的状态方程和量测方程；在坐标旋转变换滤波模型下，选取滤波算法进行滤波，得到k+1时刻目标状态的估计值；直到达到所设置的雷达跟踪时间长度，跟踪结束。

Description

一种基于坐标旋转变换的雷达跟踪方法

技术领域

本发明涉及雷达跟踪技术，属于雷达信号处理领域，具体涉及一种基于坐标旋转变换的雷达跟踪方法。 

背景技术

对目标运动状态进行估计时，学者Kalman将状态变量法引入到滤波理论中来，使状态空间描述与离散时间更新联系起来，对状态进行线性最小均方根误差估计，应用最为广泛，并称之为卡尔曼滤波方法。在雷达进行目标跟踪的过程中，由于在直角坐标系中易于对目标的运动状态进行描述，所以，目标状态方程通常是在直角坐标系中建立的。然而，对目标位置的量测通常是在极坐标系中得到的，即在极坐标系中，进行目标位置相对于雷达的距离、方位角或俯仰角（包括3D雷达的俯仰角）的量测；这就使得目标的运动状态参量（即目标动态参数）与雷达量测值之间的关系是非线性的，所以，雷达跟踪系统必然是非线性系统，这样就不能采用经典的Kalman滤波算法对目标进行跟踪。 

目前，在雷达跟踪系统中，为了解决非线性滤波问题，主要有三类解决方法：一是近似非线性方程的方法，例如基于泰勒展开式的扩展卡尔曼滤波（EKF）方法；二是估计目标状态一阶矩和二阶矩的方法，例如基于不敏变换，并沿用卡尔曼滤波框架的不敏卡尔曼滤波（UKF）算法；三是直接估计状态后验概率密度函数的方法，例如基于重要性采样原理的粒子滤波（PF）算法。然而，目前所有的滤波算法都是在事先选定的坐标系下进行目标状态估计的，这样，滤波系统的非线性程度是固定的，最终只能通过改进滤波算法来提高目标跟踪效果。 

在雷达跟踪系统中，滤波系统模型包括状态方程和量测方程。通常雷达量测值是在极坐标系下得到的，为了方便后续描述，将该极坐标系记为雷达极坐标系，将与该极坐标系对应的直角坐标系记为雷达直角坐标系XOY。 

雷达二维极坐标量测产生的基本原理为： 

假设目标是在二维平面内运动，目标的观测值是在雷达极坐标下得到的，雷达极坐标的原点与雷达直角坐标系的原点重合，则目标在雷达直角坐标系下位置的真实值为： 

$\left\{\begin{array}{c}{x}_0=r_0\cos a_0\\ y_0=r_0\sin a_0\end{array}\right.---\left(1\right)$

其中，r₀为雷达极坐标下的目标真实距离，a₀为雷达极坐标下的目标真实方位角，x₀为目标在雷达直角坐标系下的横坐标，y₀为目标在雷达直角坐标系下的纵坐标。 

而在实际中，由于雷达测量精度的限制，雷达并不能获得目标真实的距离和方位角，而是包含量测误差的目标观测值，那么，目标在雷达极坐标下的量测为： 

$\left[\begin{array}{c}r\\ a\end{array}\right]=\left[\begin{array}{c}{r}_0\\ a_0\end{array}\right]+\left[\begin{array}{c}{n}_r\\ n_a\end{array}\right]---\left(2\right)$

其中，r为雷达极坐标下的量测距离，a为雷达极坐标下的量测方位角，n_r为雷达测距噪声，n_a为雷达测角噪声。 

将雷达极坐标系下的目标量测值转换到雷达直角坐标系，则二维雷达量测极-直坐标转换方程为： 

$\left[\begin{array}{c}x=r\cos a\\ y=r\sin a\end{array}\right]---\left(3\right)$

其中，x为目标在雷达直角坐标系下的量测值的横坐标，y为目标在雷达直角坐标系下的量测值的纵坐标。 

在雷达直角坐标系XOY下，滤波系统的状态方程为： 

X(k+1)＝f(X(k))+V(k)              （4） 

其中，X(k+1)为k+1时刻目标运动的状态矢量，f(·)为k时刻目标运动的状态转移函数，V(k)为k时刻目标运动的过程噪声，并假定V(k)是零均值的高斯白噪声，V(k)的方差为： 

E[V(k)V^T(j)]＝Q(k)δ_kj           （5） 

其中，V(j)为j时刻目标运动的过程噪声，Q(k)为k时刻V(k)的协方差矩阵，δ_kj为Kronecker Delta函数，其数学表示如下： 

${\mathrm{\δ}}_{\mathrm{kj}}=\left\{\begin{array}{cc}1& k=j\\ 0& k\≠\end{array}---\left(6\right)\right.$

在雷达极坐标系下，滤波系统的量测方程为： 

Z(k)＝h[x(k)]+W(k)                 （7） 

其中，k时刻的量测值Z(k)＝[r(k)a(k)^]T，k时刻的量测函数  $h\left[x\left(k\right)\right]=\left[\begin{array}{c}\sqrt{x^2\left(k\right)+y^2\left(k\right)}\\ \arctan [y\left(k\right)/x\left(k\right)]\end{array}\right],$ 该量测函数的非线性使得滤波系统的量测方程是非线性的；x(k)为k时刻在XOY下目标位置分解到x轴的值，y(k)为k时刻在XOY下的目标位置分解到y轴的值，r(k)为k时刻目标在雷达极坐标下的量测距离，a(k)为k时刻目标在雷达极坐标下的量测方位角，W(k)为k时刻的量测噪声，并假定W(k)为零均值的高斯白噪声， $W\left(k\right)=\left[\begin{array}{c}{n}_r\left(k\right)\\ n_a\left(k\right)\end{array}\right],$ W(k)的方差为： 

E(W(k)W^T(j))＝R(k)δ_kj           （8） 

其中，W(j)为j时刻的量测噪声，R(k)为k时刻W(k)的协方差矩阵，n_r(k)为k时刻雷达直角坐标系下的雷达测距噪声，n_a(k)为k时刻雷达直角坐标系下的雷达测角噪声。 

在实际的许多雷达跟踪系统中，雷达的目标动态参数与雷达量测值之间的 关系是非线性的，由雷达极坐标量测所带来的量测方程的非线性，会对滤波及相应的目标跟踪效果产生影响。目前，大部分滤波算法都是建立在公式（4）和公式（7）的基础上对目标进行跟踪的。当描述雷达目标运动状态和雷达测量值的坐标系确定后，相应的量测方程的非线性也就确定了。所以需要建立一种坐标系，在该坐标系下，在保证目标运动特性没有发生变化的情况下，来达到降低滤波系统模型中的量测方程的非线性程度的目的。 

发明内容

有鉴于此，本发明提供了一种基于坐标旋转变换的雷达跟踪方法，该方法基于随机变量相关系数和坐标旋转变换的原理，以达到在定量度量雷达滤波系统模型中的量测方程的非线性程度的基础上，降低雷达滤波系统模型中的量测方程非线性程度，从而提高雷达跟踪效果的目的。 

本发明所提供的方法的具体设计步骤如下： 

步骤S00：在雷达极坐标下的目标真实距离r₀、雷达测距噪声标准差σ_r和雷达测角噪声标准差σ_a为确定的任意数值的情况下，得到当二维雷达量测极-直坐标转换的非线性程度最小，即坐标转换的线性度最大时，雷达极坐标下的目标真实方位角a₀的取值范围。 

1）对二维雷达量测极-直坐标转换方程中的x、y分别在（r₀，a₀）处进行二元泰勒展开，并保留至一阶项，将x的一阶泰勒展开式记为随机变量g，将y的一阶泰勒展开式记为随机变量k。 

$x\≈r_0\cos a_0+(r-r_0)\cos a_0-(a-a_0)r_0\sin a_0\stackrel{\mathrm{\Δ}}{=}g---\left(9\right)$

$y\≈r_0\sin a_0+(r-r_0)\sin a_0-(a-a_0)r_0\cos a_0\stackrel{\mathrm{\Δ}}{=}k---\left(10\right)$

通常地，雷达测距噪声n_r和雷达测角噪声n_a是统计独立的，且： 

$n_r~N(0,{\mathrm{\σ}}_r^2)---\left(11\right)$

$n_a~N(0,{\mathrm{\σ}}_a^2)---\left(12\right)$

其中，σ_r为雷达测距噪声标准差，σ_a为雷达测角噪声标准差。 

由公式（2）、（9）和（11），则随机变量g服从如下正态分布： 

$g~N(r_0\cos a_0,{\mathrm{\σ}}_r^2{\left(\cos a_0\right)}^2+{\mathrm{\σ}}_a^2{r_0}^2{\left(\sin a_0\right)}^2)---\left(13\right)$

由公式（2）、（10）和（12），则随机变量k服从如下正态分布： 

$k~N(r_0\sin a_0,{\mathrm{\σ}}_r^2{\left(\sin a_0\right)}^2+{\mathrm{\σ}}_a^2{r_0}^2{\left(\cos a_0\right)}^2)---\left(14\right)$

2）由随机变量的相关系数的定义可知，随机变量x与随机变量g的相关系数ρ_xg为： 

${\mathrm{\ρ}}_{\mathrm{xg}}=\frac{E\left\{\right[x-E\left(x\right)\left]\right[g-E\left(g\right)\left]\right\}}{\sqrt{E\left\{\right[x-E\left(x\right){]}^2\}\sqrt{E\left\{\right[g-E\left(g\right){]}^2\}}}}---\left(15\right)$

结合协方差的性质，则有： 

${\mathrm{\ρ}}_{\mathrm{xg}}=\frac{E\left(\mathrm{xg}\right)-E\left(x\right)E\left(g\right)}{\sqrt{E\left(x^2\right)-[E\left(x\right){]}^2\sqrt{E\left\{\right[g-E\left(g\right){]}^2\}}}}---\left(16\right)$

根据文献M.Miller and0.Dmmmond.Coordinate Transformation Bias in Target Tracking.In Proceedings of SPIE Conference on Signal and Data Processing of Small Targets1999,pages409424,1999.SPIE Vol.3809.中的记载，可知： 

$E\left[\cos a\right]=e^{-\frac{1}{2}{\mathrm{\σ}}_a^2}\cos a_0---\left(17\right)$

$E\left[\sin a\right]=e^{-\frac{1}{2}{\mathrm{\σ}}_a^2}\sin a_0---\left(18\right)$

$E\left[n_a\cos a\right]={-e}^{-\frac{1}{2}{\mathrm{\σ}}_a^2}{\mathrm{\σ}}_a^2\sin a_0---\left(19\right)$

$E\left[n_a\sin a\right]=e^{-\frac{1}{2}{\mathrm{\σ}}_a^2}{\mathrm{\σ}}_a^2\cos a_0---\left(20\right)$

根据公式（17）可得： 

${\cos }^{2}a]=E[\frac{1}{2}1+\cos 2a]=\frac{1}{2}[1+e^{-2{\mathrm{\σ}}_a^2}\cos 2a_0]---\left(21\right)$

根据公式（18）可得： 

$E\left[{\sin }^{2}a\right]=E[\frac{1}{2}1-\cos 2a]=\frac{1}{2}[1-e^{-2{\mathrm{\σ}}_a^2}\cos 2a_0]---\left(22\right)$

由公式（2）和（11），则： 

$r~N(r0,{\mathrm{\σ}}_r^2)---\left(23\right)$

由公式（2）和（12），则： 

$a~N(a0,{\mathrm{\σ}}_a^2)---\left(24\right)$

由公式（2）、（3）、（11）、（17）、（19）、（21）和（23），并结合公式（16），则ρ_xg的解析表达式如下： 

${\mathrm{\ρ}}_{\mathrm{xg}}=\sqrt{\frac{{\mathrm{\σ}}_r^2{\left(\cos a_0\right)}^2+{r_0}^2{\mathrm{\σ}}_a^2{\left(\sin a_0\right)}^2}{({r_0}^2+{\mathrm{\σ}}_r^2)0.5[e^{{\mathrm{\σ}}_a^2}+e^{-{\mathrm{\σ}}_a^2}\cos 2a_0]-{r_0}^2{\left(\cos a_0\right)}^2}}---\left(25\right)$

同理，随机变量y与随机变量k的相关系数ρ_yk为： 

${\mathrm{\ρ}}_{\mathrm{xg}}=\sqrt{\frac{{\mathrm{\σ}}_r^2{\left(\sin a_0\right)}^2+{r_0}^2{\mathrm{\σ}}_a^2{\left(\cos a_0\right)}^2}{({r_0}^2+{\mathrm{\σ}}_r^2)0.5[e^{{\mathrm{\σ}}_a^2}+e^{-{\mathrm{\σ}}_a^2}\cos 2a_0]-{r_0}^2{\left(\sin a_0\right)}^2}}---\left(26\right)$

3）取经步骤S00的第2）步得到的ρ_xg和ρ_yk的最小值，定义二维雷达量测极-直坐标转换线性度ρ为： 

ρ＝min(ρ_xg,ρ_yk)             （27） 

根据随机变量坐标转换相关系数ρ的定义可知，ρ是一个无量纲的量，并且0≤ρ≤1。 

由式（25）、（26）和（27）可知，二维雷达量测极-直坐标转换的非线性与雷达极坐标下的目标真实距离r₀、雷达极坐标下的目标真实方位角a₀、雷达测距噪声标准差σ_r和雷达测角噪声标准差σ_a有关。 

当r₀、σ_r和σ_a为确定的任意数值时，二维雷达量测极-直坐标转换线性度ρ在 

时取得最大值。 

步骤S01：当雷达极坐标下的目标真实距离r₀、雷达极坐标下的目标真实方位角a₀、雷达测距噪声标准差σ_r和雷达测角噪声标准差σ_a为确定的任意数值时，通过将雷达直角坐标系XOY顺时针旋转角度使得旋转后的直角坐标系，即旋转雷达直角坐标系X^cOY^c，与X^cOY^c对应的极坐标系，即旋转雷达极坐标系，之间的极-直坐标转换线性度最大。 

由步骤S00可知，当r₀、σ_r和σ_a固定时，二维雷达量测极-直坐标转换线性度ρ在

时取得最大值。然而在雷达直角坐标系XOY下，当r₀、σ_r和σ_a固定时，雷达极坐标下的目标真实方位角a₀不一定满足

使二维雷达量测极-直坐标转换线性度ρ取最大值，故需要将雷达直角坐标系XOY顺时针旋转角度

变为旋转雷达直角坐标系X^cOY^c，相应地，雷达直角坐标量测变为： 

其中，x^c为目标在旋转雷达直角坐标系X^cOY^c下的量测值的横坐标；y^c为目标在旋转雷达直角坐标系X^cOY^c下的量测值的纵坐标；坐标旋转矩阵为： 

由公式（3）、（28）和（29），则在旋转雷达直角坐标系X^cOY^c下的雷达直角坐标量测值为： 

由公式（3）和（30）可知，雷达直角坐标系XOY顺时针旋转角度后，目标在雷达极坐标下的量测[ra]^T转换为旋转雷达极坐标系X^cOY^c下的

因此，二维雷达量测极-直坐标转换线性度ρ(a₀,r₀,σ_r,σ_a)转换为

 相应地，在步骤S00中得到的：当r₀、σ_r和σ_a固定时，二维雷达量测极-直坐标转换线性度ρ(a₀,r₀,σ_r,σ_a)在

时取得最大值，转换为：当r₀、a₀、σ_r和σ_a固定时，二维旋转雷达量测极-直坐标转换线性度

在满足公式（31）时取得最大值，即当r₀、a₀、σ_r和σ_a固定时，通过选择旋转角度

满足公式（31），使得二维旋转雷达量测极-直坐标转换线性度

最大，进而达到当r₀、σ_r和σ_a固定时，通过坐标旋转变换，使得旋转雷达极坐标系与旋转雷达直角坐标系X^cOY^c之间的极-直坐标转换线性取得最大值。 

由于直-极坐标转换，即滤波系统的量测函数是极-直坐标转换的反函数，因而当极-直坐标转换的线性程度大时，则直-极坐标转换的线性程度亦大。因而，可以通过坐标旋转变换使得极-直坐标转换的线性度最大，来降低滤波系统的量测方程的非线性程度。 

步骤S02：设置目标运动状态，包括目标运动的初始位置和速度、雷达跟踪时间长度和采样时间间隔，进行滤波初始化。 

步骤S03：目标运动状态由k时刻递推得到k+1时刻的实现步骤。 

1）由于实际中无法得到雷达极坐标系下的目标真实方位角a₀，因而将k+1时刻的雷达极坐标系下的目标方位角的量测值a(k+1)代入式（31）求得k+1时刻的坐标旋转变换的角度值

2）根据所得到的

将雷达直角坐标系XOY顺时针旋转角度值

得到k+1时刻的旋转雷达直角坐标系X^cOY^c和相应的旋转雷达极坐标系。 

3）根据所得到的k+1时刻的旋转雷达直角坐标系X^cOY^c和相应的旋转雷达极坐标系，通过分析坐标旋转变换前后坐标系之间的关系，即雷达直角坐标系和旋转雷达直角坐标系、雷达极坐标系和旋转雷达极坐标系之间的关系，相应 地得到旋转雷达直角坐标系下的状态方程和旋转雷达极坐标系下的量测方程，并将二者记为坐标旋转变换滤波系统模型。 

a、坐标旋转变换滤波系统模型的状态方程 

①由公式（28）可知：在k+1时刻目标的位置在雷达直角坐标系XOY下与其在旋转雷达直角坐标系X^cOY^c下的关系为： 

其中，x^c(k+1)为在X^cOY^c下目标位置分解到x轴的值；y^c(k+1)为在X^cOY^c下目标位置分解到y轴的值；x(k+1)为在XOY下目标位置分解到x轴的值；y(k+1)为在XOY下目标位置分解到y轴的值；

为对应于坐标旋转变换的角度值

的坐标旋转矩阵，即： 

在k+1时刻目标的速度在XOY下与其在X^cOY^c下的关系为： 

其中，为在X^cOY^c下目标速度分解到x轴的值；

为在X^cOY^c下目标速度分解到y轴的值；v_x(k+1)为在XOY下目标速度分解到x轴的值；v_y(k+1)为在XOY下目标速度分解到y轴的值。 

在k+1时刻目标的加速度在XOY下与其在X^cOY^c下的关系为： 

其中，

为在X^cOY^c下目标加速度分解到x轴的值；

为在X^cOY^c下目标加速度分解到y轴的值；a_x(k+1)为在XOY下目标加速度分解到x轴的值；a_y(k+1)为在XOY下目标加速度分解到y轴的值。 

②若在雷达直角坐标系XOY下，仅考虑目标位置和目标速度，即d＝2，则目标运动状态向量X(k+1)为： 

X(k+1)＝[x(k+1)v_x(k+1)y(k+1)v_y(k+1)]^T   （36） 

由公式（32）和（33）可知旋转雷达直角坐标系X^cOY^c下的目标运动状态向量X^c(k+1)与X(k+1)的关系为： 

X^c(k+1)＝AX(k+1)                （37） 

其中， 

为Kronecker积，I_d为d×d维单位矩阵，由d＝2，则取I₂，进一步地，由公式（38），则有： 

若在雷达直角坐标系XOY下，考虑目标位置、目标速度和目标加速度，则d＝3，在公式（38）中，I_d取I₃。 

③经过步骤②的坐标旋转变换，即由雷达直角坐标系XOY变为旋转雷达直角坐标系X^cOY^c，相应地，由公式（4），坐标旋转变换滤波系统模型的状态方程为： 

X^c(k+1)＝f^c(X^c(k))+V^c(k)          （40） 

其中，X^c(k+1)为在k+1时刻X^cOY^c下的目标运动的状态矢量，f^c(·)为在k时刻X^cOY^c下的目标运动的状态转移函数，V^c(k)为在k时刻X^cOY^c下的目标运动的过程噪声。 

由公式（37）可知：在k时刻，旋转雷达直角坐标系X^cOY^c下的目标运动的状态转移函数f^c(X^c(k))与其在雷达直角坐标系XOY下的关系为： 

f^c(X^c(k))＝Af(A^-1X^c(k))        （41） 同理可得，在k时刻，旋转雷达直角坐标系X^cOY^c下的目标运动的过程噪声V^c(k)与其在雷达直角坐标系XOY下的关系为： 

V^c(k)＝AV(k)                 （42） 

V^c(k)的方差为： 

E[V^c(k)(V^c(j))^T]＝AQ(k)A^Tδ_kj   （43） 

其中，V^c(j)为在j时刻X^cOY^c下的目标运动的过程噪声。 

由于坐标旋转并没有改变目标的运动特性，因此，若目标在雷达直角坐标系下做匀速运动，那么目标在旋转雷达直角坐标系下仍然做匀速运动，只是在不同的雷达直角坐标系下，目标的运动状态分解到相应的各个坐标轴的值不一样。 

b、坐标旋转变换滤波系统模型的量测方程 

经过步骤S03的步骤2）的坐标旋转变换，由公式（7），在k+1时刻，坐标旋转变换滤波系统模型的量测方程为： 

Z^c(k+1)＝h^c[X^c(k+1)]+W^c(k+1)          （44） 

其中，k+1时刻的旋转量测值Z^c(k+1)＝[r^c(k+1)a^c(k+1)]^T，r^c(k+1)为k+1时刻目标在旋转雷达极坐标下的量测距离，a^c(k+1)为k+1时刻目标在旋转雷达极坐标系下的量测方位角；W^c(k1)为k+1时刻旋转雷达直角坐标系X^cOY^c下的量测噪声，由于雷达直角坐标转换到雷达极坐标的坐标转换关系与旋转雷达直角坐标转换到旋转雷达极坐标的坐标转换关系是相同的，则有： 

h^c(·)＝h(·)                （45）即： 

$h^c\left[X^c(k+1)\right]=\left[\begin{array}{c}\sqrt{{\left[x^c(k+1)\right]}^2+{\left[y^c(k+1)\right]}^2}\\ \arctan [y^c(k+1)/x^c(k+1)]\end{array}\right]---\left(46\right)$

由于旋转雷达极坐标系与雷达极坐标系的坐标原点是相同，故坐标旋转不改变目标量测距离，即r^c(k+1)＝r(k+1)；而目标在旋转雷达极坐标系中的量测方位角a^c(k+1)相比于雷达极坐标系下的量测方位角发生了变化，变化的角度值为坐标旋转角度值

并且测距噪声和测角噪声并不随着坐标旋转的变化而变化；对比公式（7），则有： 

Z^c(k+1)＝h^cX^c(k+1)]+W^c(k)        （47） 

Z^c(k+1)＝Z(k+1)+B(k+1)              （48） 

W^c(k)＝W(k)                     （50） 

其中，W^c(k)为k时刻旋转雷达直角坐标系X^cOY^c下的量测噪声。 

步骤S04:在坐标旋转变换滤波模型下，选取滤波算法进行滤波，得到k+1时刻目标状态的估计值。 

步骤S05:重复步骤S03和步骤S04，直到达到步骤2）设置的雷达跟踪时间长度，跟踪结束。 

有益效果： 

本发明针对雷达极坐标量测方程的非线性问题，采用二维雷达量测极-直坐标转换后的变量与二维雷达量测极-直坐标转换方程的一阶泰勒展开式之间的相关系数进行定量度量，并且定义二维雷达极-直坐标转换线性度为两个相关系数的最小值，将该最小值作为度量二维雷达量测极-直坐标转换方程的非线性程度的数值；在给出的计算二维雷达量测坐标转换非线性程度的公式中，得到影响雷达极坐标下目标的量测值与雷达直角坐标下目标动态参数间的非线性程度的因素，包括雷达极坐标系下的真实距离、雷达极坐标系下的目标真实方位角、雷达测距噪声标准差以及雷达测角噪声标准差，其中，雷达极坐标系下的目标 真实方位角对坐标转换非线性的的影响处于主导地位。雷达极坐标下目标的量测值与雷达直角坐标下目标动态参数间的非线性程度随着目标真实方位角（记为a₀）的正弦的平方sin²(a₀)变化而变化，当0≤sin²(a₀)≤12时，坐标转换的非线性程度随着sin²(a₀)的增大而降低。当12≤sin²(a₀)≤1时，坐标转换的非线性程度随着sin²(a₀)的增大而增大，其中，坐标转换的非线性程度在sin²(a₀)＝12时最小。 

在本发明给出的坐标旋转变换滤波系统模型下，将与雷达极坐标系对应的雷达直角坐标系旋转一定的角度

使得

从而使旋转后雷达直角坐标系对应的雷达极坐标系中的量测方程的非线性程度降到最低。同时，由于坐标旋转变换是线性的，在进行坐标旋转之后，目标运动特性并没有发生改变，即通过使用本发明所提供的坐标旋转滤波系统模型，可以在不改变目标运动特性的情况下，大大地降低量测方程的非线性程度，这样就降低了滤波系统固有的非线性程度，从而达到了提高雷达跟踪效果的目的。 

此外，在坐标旋转变换滤波系统模型下可以应用所有的非线性滤波算法，在该系统模型下，所有的非线性滤波算法的目标滤波效果相比于不经过坐标旋转变换的滤波系统模型下的滤波效果得到了显著的提高。 

附图说明

图1为本发明实施例中雷达极坐标与雷达直角坐标对应关系示意图； 

图2为本发明所提供的二维坐标旋转示意图； 

图3为滤波位置均方根误差对比仿真图； 

图4为滤波ANEES对比仿真图。 

具体实施方式

下面结合附图并举实施例，对本发明进行详细描述。 

图1为本发明所提供的雷达极坐标与雷达直角坐标对应关系示意图，图2为二维坐标旋转示意图，假设目标在二维平面内做匀速直线运动，具体步骤如下： 

步骤S00：当雷达极坐标下的目标真实距离r₀、雷达测距噪声标准差σ_r和雷达测角噪声标准差σ_a为确定的任意数值时，由下式求解二维雷达量测极-直坐标转换的线性度ρ， 

ρ＝min(ρ_xg,ρ_yk)              （1） 

可得，二维雷达量测极-直坐标转换的线性度ρ在

时取得最大值。 

步骤S01：当雷达极坐标下的目标真实距离r₀、雷达极坐标下的目标真实方位角a₀、雷达测距噪声标准差σ_r和雷达测角噪声标准差σ_a为确定的任意数值时，通过将雷达直角坐标系XOY顺时针旋转角度使得旋转雷达极坐标系与旋转雷达直角坐标系X^cOY^c之间的极-直坐标转换线性取得最大值；其中，

满足下式： 

步骤S02：设置目标的初始位置为（50Km，0.5Km）、目标的初始速度为（100m/s，100m/s），相应地，d＝2；设置雷达跟踪时间长度为40s，采样时间间隔为0.1s。并且采用两点起始法进行滤波初始化。 

步骤S03：目标状态由k时刻递推得到k+1时刻的实现步骤。 

1)将k+1时刻的目标方位角的量测值a(k+1)代入公式（2）求得k+1时刻的坐标旋转变换的角度值

2)根据所得到的

将雷达直角坐标系XOY顺时针旋转角度值

得到k+1时刻的旋转雷达直角坐标系X^cOY^c和相应的旋转雷达极坐标系。 

3)根据所得到的k+1时刻的旋转雷达直角坐标系X^cOY^c和相应的旋转雷达极 坐标系，得到旋转雷达直角坐标系下的状态方程和旋转雷达极坐标系下的量测方程。 

a、坐标旋转变换滤波系统模型的状态方程 

①若在雷达直角坐标系XOY下，仅考虑目标位置和目标速度，即d＝2，则目标的运动状态向量X(k+1)为： 

X(k+1)＝[x(k+1)v_x(k+1)y(k+1)v_y(k+1)]^T     （3） 

旋转雷达直角坐标系X^cOY^c下的目标运动状态向量X^c(k+1)与X(k+1)的关系为： 

X^c(k+1)＝AX(k+1)                  （4） 

由d＝2，则取I₂，进一步地， 

②由于假定目标在二维平面内做匀速直线运动。因而，坐标旋转变换滤波系统模型的状态方程为： 

X^c(k+1)＝F^cX^c(k)+V^c(k)              （6） 

其中， 

$X^c\left(k\right)=[x^c\left(k\right)v_x^c\left(k\right)y^c\left(k\right)v_y^c\left(k\right){]}^T---\left(7\right)$

根据步骤S03的步骤1）得到的坐标旋转的角度值

并由公式（5）求得A。并且公式（6）中的状态转移矩阵为: 

$F^c=\left[\begin{array}{cccc}1& 0& T& 0\\ 0& 1& 0& T\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]---\left(8\right)$

其中，T为采样间隔，取T＝0.1s，假定旋转雷达直角坐标系下各个坐标轴的过程噪声为高斯白噪声，那么，公式（6）中的V^c(k)为， 

$V^c\left(k\right)=\left[\begin{array}{cc}{0.5T}^2& 0\\ T& 0\\ 0& {0.5T}^2\\ 0& T\end{array}\right]\left[\begin{array}{c}{n}_x\left(k\right)\\ n_y\left(k\right)\end{array}\right]---\left(9\right)$

其中，n_x(k)为k时刻旋转雷达直角坐标系下目标运动的过程噪声分解到x轴的值，n_y(k)为k时刻旋转雷达直角坐标系下目标运动的过程噪声分解到y轴的值，n_x和n_y的均值均为0，n_x的方差σ_x＝0.01m²/s，n_y的方差σ_y＝0.01m²/s。 

化简式（6）得到本实施例的坐标旋转变换滤波系统模型的状态方程为： 

$\left[\begin{array}{c}{x}^c(k+1)\\ v_x^c(k+1)\\ y^c(k+1)\\ v_y^c(k+1)\end{array}\right]=\left[\begin{array}{cccc}1& 0& T& 0\\ 0& 1& 0& T\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x}^c\left(k\right)\\ v_x^c\left(k\right)\\ y^c\left(k\right)\\ v_y^c\left(k\right)\end{array}\right]+\left[\begin{array}{cc}{0.5T}^2& 0\\ T& 0\\ 0& {0.5T}^2\\ 0& T\end{array}\right]\left[\begin{array}{c}{n}_x\left(k\right)\\ n_y\left(k\right)\end{array}\right]---\left(10\right)$

b、坐标旋转变换滤波系统模型的量测方程 

假设k时刻旋转雷达直角坐标系X^cOY^c下的量测噪声W^c(k)中的雷达测距噪声n_r(k)为零均值的高斯白噪声，设置其标准差σ_r＝50m，雷达测角噪声n_a(k)为零均值的高斯白噪声，设置其标准差σ_a＝0.5°，其中： 

$W^c\left(k\right)=W\left(k\right)=\left[\begin{array}{c}{n}_r\left(k\right)\\ n_a\left(k\right)\end{array}\right]---\left(11\right)$

将步骤S03的步骤1）得到的坐标旋转的角度值

代入下式： 

将Z(k+1)＝[r(k+1)a(k+1)]^T以及由上式计算得到的B(k+1)，代入下式： 

Z^c(k+1)＝Z(k+1)+B(k+1)              （13） 

将 $h^c{X}^{c}k+1=\left[\begin{array}{c}\sqrt{x^c{(k+1)}^2+y^c{(k+1)}^2}\\ \arctan [y^c(k+1)/x^c(k+1)]\end{array}\right]`$ 公式（13）以及由式（11）计算得到的W^c(k)代入下式： 

Z^c(k+1)＝h^cX^c(k+1)]+W^c(k)         （14） 

化简式（14）得到本实施例的坐标旋转变换滤波系统模型的量测方程为： 

步骤S04：选择一阶扩展卡尔曼滤波算法（FEKF）和不敏卡尔曼滤波算法（UKF）两种非线性滤波算法，采用MATLAB软件分别进行仿真，并在坐标旋转变换滤波模型下进行滤波，得到k+1时刻目标状态的估计值。 

步骤S05：重复步骤S03和步骤S04，直至达到步骤S02所设置的雷达跟踪时间长度40s时，跟踪结束。 

图3为滤波位置均方根误差对比仿真图，图4为滤波ANEES对比仿真图，本实施例的仿真参数如下： 

表1仿真参数 

假设目标在二维平面内做匀速直线运动： 

图3的仿真结果显示了，采用本发明所提供的坐标旋转变换滤波系统模型的滤波算法：FEKF、UKF并分别记为CFEKF和CUKF，以及未采用本发明所提供的坐标旋转变换滤波系统模型的滤波算法：FEKF、UKF和二阶扩展卡尔曼滤波算法（SEKF）各种滤波算法得到的目标状态估计的位置均方根误差曲线分别与目标状态估计的位置的后验克拉美罗界（PCRLB）之间的关系图。 

其中，目标状态估计的位置的PCRLB的表示了各种滤波算法的下界，滤波位置均方根误差越小越好，当某种滤波算法的滤波均方根误差曲线达到PCRLB时，该算法最优。 

图4的仿真结果显示了，采用本发明所提供的坐标旋转变换滤波系统模型的滤波算法：CFEKF和CUKF，以及未采用本发明所提供的坐标旋转变换滤波系统模型的滤波算法：FEKF、UKF和二阶扩展卡尔曼滤波算法（SEKF）各种滤波算法的滤波ANEES仿真结果图。图4表示对各个滤波算法的一致性的检验，当滤波算法落入两条平行于x轴的直线之间，则表示该算法满足一致性。 

由图3和图4可知，UKF的滤波效果优于SEKF，SEKF优于FEKF；CFEKF和CUKF的滤波效果优于FEKF、SEKF和UKF，并且，CFEKF滤波曲线由最初的接近于PCRLB到最终达到PCRLB，CUKF滤波曲线始终与PCRLB满足一致性。因而基于坐标旋转变换滤波系统模型的滤波算法的目标状态估计效果要远远优于不采用坐标旋转变换的滤波系统模型的滤波算法。 

以上所述的具体描述，对发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。 

Claims (1)

1.一种基于坐标旋转变换的雷达跟踪方法，其特征在于，包括：

步骤S00、二维雷达量测极-直坐标转换方程为：

$\left\{\begin{array}{c}x=r\cos a\\ y=r\sin a\end{array}\right.---\left(1\right)$

其中，x为目标在雷达直角坐标系下的量测值的横坐标，y为目标在雷达直角坐标系下的量测值的纵坐标；

对式（2）中的x、y分别在（r₀，a₀）处进行二元泰勒展开，并保留至一阶项，将x的一阶泰勒展开式记为随机变量g，将y的一阶泰勒展开式记为随机变量k；

$x\≈r_0\cos a_0+(r-r_0)\cos a_0-(a-a_0)r_0\sin a_0\stackrel{\mathrm{\Δ}}{=}g---\left(2\right)$

$y\≈r_0\sin a_0+(r-r_0)\sin a_0-(a-a_0)r_0\cos a_0\stackrel{\mathrm{\Δ}}{=}k---\left(3\right)$

其中，r₀为雷达极坐标下的目标真实距离，a₀为雷达极坐标下的目标真实方位角，

定义二维雷达量测极-直坐标转换线性度ρ为ρ_xg和ρ_yk二者的最小值：

ρ＝min(ρ_xg,ρ_yk)               （4）

其中：ρ_xg为随机变量x与随机变量g的相关系数，ρ_yk为随机变量y与随机变量k的相关系数，且

${\mathrm{\ρ}}_{\mathrm{xg}}=\sqrt{\frac{{\mathrm{\σ}}_r^2{\left(\cos a_0\right)}^2+{r_0}^2{\mathrm{\σ}}_a^2{\left(\sin a_0\right)}^2}{({r_0}^2+{\mathrm{\σ}}_r^2)0.5[e^{{\mathrm{\σ}}_a^2}+e^{{-\mathrm{\σ}}_a^2}\cos {2a}_0]-{r_0}^2{\left(\cos a_0\right)}^2}}---\left(5\right)$

${\mathrm{\ρ}}_{\mathrm{yk}}=\sqrt{\frac{{\mathrm{\σ}}_r^2{\left(\sin a_0\right)}^2+{r_0}^2{\mathrm{\σ}}_a^2{\left(\cos a_0\right)}^2}{({r_0}^2+{\mathrm{\σ}}_r^2)0.5[e^{{\mathrm{\σ}}_a^2}+e^{{-\mathrm{\σ}}_a^2}\cos {2a}_0]-{r_0}^2{\left(\sin a_0\right)}^2}}---\left(6\right)$

σ_r为雷达测距噪声标准差，σ_a为雷达测角噪声标准差；

步骤S01、由式（4），当r₀、σ_r和σ_a为确定的任意数值时，二维雷达量测极-直坐标转换线性度ρ在

时取得最大值，此时，将雷达直角坐标系XOY顺时针旋转角度

变为旋转雷达直角坐标系X^cOY^c，相应地，雷达直角坐标量测变为：

其中，x^c为目标在旋转雷达直角坐标系X^cOY^c下的量测值的横坐标；y^c为目标在旋转雷达直角坐标系X^cOY^c下的量测值的纵坐标；坐标旋转矩阵

为：

由式（7）和（8），则在旋转雷达直角坐标系X^cOY^c下的雷达直角坐标量测值为：

由式（9），并结合所述当r₀、σ_r和σ_a为确定的任意数值时，二维雷达量测极-直坐标转换线性度ρ在时取得最大值，通过坐标旋转变换，使得：

步骤S02、设置目标运动状态，包括目标运动的初始位置、目标运动的速度、雷达跟踪时间长度和采样时间间隔，进行滤波初始化；

步骤S03、目标运动状态由^k时刻递推得到k+1时刻的具体步骤如下：

①将k+1时刻的雷达极坐标系下的目标方位角的量测值a(k+1)代入式（10）求得k+1时刻的坐标旋转变换的角度值

②将雷达直角坐标系XOY顺时针旋转角度值

建立k+1时刻的旋转雷达直角坐标系X^cOY^c和相应的旋转雷达极坐标系；

③在k+1时刻，分别分析雷达直角坐标系XOY和旋转雷达直角坐标系X^cOY^c、雷达极坐标系和旋转雷达极坐标系之间的关系，相应地得到旋转雷达直角坐标系X^cOY^c下的状态方程和旋转雷达极坐标系下的量测方程，并将二者记为坐标旋转变换滤波系统模型；

a、坐标旋转变换滤波系统模型的状态方程

结合式（7），则：在k+1时刻目标的位置在雷达直角坐标系XOY下与其在旋转雷达直角坐标系X^cOY^c下的关系为：

其中，x^c(k+1)为在旋转雷达直角坐标系X^cOY^c下目标位置分解到x轴的值；y^c(k+1)为在旋转雷达直角坐标系X^cOY^c下目标位置分解到y轴的值；x(k+1)为在雷达直角坐标系XOY下目标位置分解到x轴的值；y(k+1)为在雷达直角坐标系XOY下目标位置分解到y轴的值；

为对应于坐标旋转变换的角度值

的坐标旋转矩阵，即：

在k+1时刻目标的速度在雷达直角坐标系XOY下与其在旋转雷达直角坐标系X^cOY^c下的关系为：

其中，为在旋转雷达直角坐标系X^cOY^c下目标速度分解到x轴的值；

为在旋转雷达直角坐标系X^cOY^c下目标速度分解到y轴的值；v_x(k+1)为在雷达直角坐标系XOY下目标速度分解到x轴的值；v_y(k+1)为在雷达直角坐标系XOY下目标速度分解到y轴的值；

在k+1时刻目标的加速度在XOY下与其在X^cOY^c下的关系为：

其中，为在旋转雷达直角坐标系X^cOY^c下目标加速度分解到x轴的值；

为在旋转雷达直角坐标系X^cOY^c下目标加速度分解到y轴的值；a_x(k+1)为在雷达直角坐标系XOY下目标加速度分解到x轴的值；a_y(k+1)为在雷达直角坐标系XOY下目标加速度分解到y轴的值；

若在雷达直角坐标系XOY下，仅考虑目标位置和目标速度，即d＝2，则目标运动状态向量X(k+1)为：

X(k+1)＝[x(k+1)v_x(k+1)y(k+1)v_y(k+1)]^T   （15）

由式（11）和（12），则旋转雷达直角坐标系X^cOY^c下的目标运动状态向量X^c(k+1)与X(k+1)的关系为：

X^c(k+1)＝AX(k+1)                （16）

其中，

为Kronecker积，I_d为d×d维单位矩阵，由d＝2，则取I₂，进一步地，由式（17），则有：

若在雷达直角坐标系XOY下，考虑目标位置、目标速度和目标加速度，则d＝3，相应地，在式（17）中，I_d取I₃；

经过坐标旋转变换，即由雷达直角坐标系XOY变为旋转雷达直角坐标系X^cOY^c，相应地，坐标旋转变换滤波系统模型的状态方程为：

X^c(k+1)＝f^c(X^c(k))+V^c(k)         （19）

其中，X^c(k+1)为在k+1时刻旋转雷达直角坐标系X^cOY^c下的目标运动的状态矢量，f^c(·)为在k时刻旋转雷达直角坐标系X^cOY^c下的目标运动的状态转移函数，V^c(k)为在k时刻旋转雷达直角坐标系X^cOY^c下的目标运动的过程噪声；

结合式（16），在k时刻，旋转雷达直角坐标系X^cOY^c下的目标运动的状态转移函数f^c(X^c(k))与其在雷达直角坐标系XOY下的关系为：

f^c(X^c(k))＝Af(A^-1X^c(k))       （20）

同理，在k时刻，旋转雷达直角坐标系X^cOY^c下的目标运动的过程噪声V^c(k)与其在雷达直角坐标系XOY下的关系为：

V^c(k)＝AV(k)                  21）

假定V^c(k)是零均值的高斯白噪声，结合式（21），则V^c(k)的方差为：

E[V^c（k）(V^c(j）)^T]＝AQ(k）A^Tδ_kj       （22）

其中，V^c(j)为在j时刻旋转雷达直角坐标系X^cOY^c下的目标运动的过程噪声；

b、坐标旋转变换滤波系统模型的量测方程

经过步骤S03中步骤②的坐标旋转变换，在k+1时刻，坐标旋转变换滤波系统模型的量测方程为：

Z^c（k+1）＝h^c[X^c（k+1）]+W^c（k+1）               （23）

其中，k+1时刻的旋转量测值Z^c(k+1)＝[r^c(k+1)a^c(k+1)]^T，r^c（k+1）为k+1时刻目标在旋转雷达极坐标下的量测距离，a^c（k+1）为k+1时刻目标在旋转雷达极坐标系下的量测方位角；W^c（k+1）为k+1时刻旋转雷达直角坐标系X^cOY^c下的量测噪声；

由于雷达直角坐标转换到雷达极坐标的坐标转换关系与旋转雷达直角坐标转换到旋转雷达极坐标的坐标转换关系是相同的，则有：

h^c(·)＝h(·)                 （24）即：

$h^c\left[X^c(k+1)\right]=\left[\begin{array}{c}\sqrt{[x^c(k+1){]}^2+[y^c(k+1){]}^2}\\ \arctan [y^c(k+1)/x^c(k+1)]\end{array}\right]---\left(25\right)$

由于旋转雷达极坐标系与雷达极坐标系的坐标原点是相同，故坐标旋转不改变目标量测距离，即r^c(k+1)＝r(k+1)；而目标在旋转雷达极坐标系中的量测方位角a^c(k+1)相比于雷达极坐标系下的量测方位角发生了变化，变化的角度值为坐标旋转角度值

并且测距噪声和测角噪声并不随着坐标旋转的变化而变化，则有：

Z^c(k+1)＝h^cX^c(k+1)]+W^c(k)       （26）

Z^c(k+1)＝Z(k+1)+B(k+1)           （27）

W^c（k）＝W（k）               （29）

其中，W^c(k)为k时刻旋转雷达直角坐标系X^cOY^c下的量测噪声；

步骤S04、在坐标旋转变换滤波模型下，选取滤波算法进行滤波，得到k+1时刻目标状态的估计值；

步骤S05、重复步骤S03和步骤S04，直到达到步骤S02中设置的雷达跟踪时间长度，跟踪结束。

Images