CN105548985A - 基于RAV-Jerk模型的机动目标跟踪方法 - Google Patents

Abstract

本发明公开了一种基于RAV-Jerk模型机动目标跟踪方法，属于雷达机动目标跟踪领域。本发明的方法可快速准确地提供机动目标的径向加速度和径向速度信息，有效提高雷达对机动目标的跟踪性能。本发明的方法主要包括以下步骤：（一）对雷达接收信号进行采样，利用正交匹配追踪（OMP）方法得到目标径向加速度和径向速度；（二）在数据处理阶段将径向加速度和径向速度引入到Jerk机动模型的量测方程和状态方程；（三）采用扩展卡尔曼滤波（EKF）算法实现机动目标跟踪。本发明能够精确实时地反映目标的机动情况，提高了目标跟踪精度，改善了速度和加速度估计精度，工程实现容易，具有较强的工程应用价值和推广前景。

Description

基于RAV-Jerk模型的机动目标跟踪方法

技术领域

本发明隶属于雷达机动目标跟踪领域，适用于高/中脉冲重复频率雷达(如机载脉冲多普勒雷达等)对机动目标的精确跟踪。

背景技术

随着科学技术的发展，现代军用飞机和导弹的机动能力大大增强，对雷达的检测和跟踪性能提出了新的挑战，对高机动目标跟踪也成为雷达目标跟踪领域的难点和重点。现有的机动目标跟踪技术是在雷达测量信息(位置、多普勒速度)的基础上，利用各种机动模型和滤波算法实现对机动目标的跟踪，其中，Singer模型方法是一种常用的有效机动目标跟踪方法，主要由以下3个步骤实现：

(1)将雷达接收机输出的目标回波信号进行A/D变换，送雷达数据处理计算机执行以下步骤；

(2)将机动加速度作为零均值时间相关色噪声建模。设机动加速度的概率密度函数在[-A_max,A_max]上近似服从零均值均匀分布，并假设其时间相关函数为指数衰减形式；

(3)采用卡尔曼滤波算法对机动目标进行自适应跟踪，估计出目标的位置、速度和加速度状态。

这种方法具有以下缺陷：

由于现有雷达在目标跟踪时只能利用目标距离和方位角测量信息，只有多普勒雷达可以利用径向速度和距离、方位角测量信息，因此当目标发生机动时，加速度发生剧变，雷达由于缺少加速度测量信息，会出现跟踪精度低甚至跟踪发散的问题。

发明内容

本发明的目的是提出一种基于RAV-Jerk模型机动目标跟踪方法，将信号处理阶段得到的径向加速度和径向加速度信息引入到机动目标跟踪模型中，解决现有机动模型方法对目标突发机动的自适应跟踪能力不强，以及滤波后速度和加速度估计误差较大的问题。

本发明提出的基于RAV-Jerk模型机动目标跟踪方法的技术方案包括以下步骤：

步骤1：将雷达接收机接收到的线性调频信号s(t)通过采样器以采样间隔T_s进行采样，变为离散信号s(nT_s)，其中n表示采样点序号；将s(nT_s)送入雷达信号处理计算机，初始化分解参数：

T设为雷达的脉冲宽度；

λ为雷达波长；

f_s为采样频率；

f_u设为LFM信号的初始频率；

k_v设为LFM信号的调频率；

G((U×V)×N))设为过完备原子库，U×V为原子库中原子的个数，N＝T/f_s；

λ²设为判断稀疏分解是否完成的相干比阈值；

在雷达信号处理计算机中执行以下步骤：

步骤2：根据机动目标信号回波s(t)的特点，建立与之相匹配的原子g_r，构造出过完备原子库G

(1)根据LFM信号回波s(t)的特点，建立原子g_r＝exp[j2π(f_un+k_vn²/2N)]，r＝1,2,…,N；

(2)设定搜索精度和范围，假设搜索范围f_u的取值为f_u∈[0,U]△f_u，u＝1,2,…,U，U为起始频率的搜索个数，△f_u为多普勒单元，搜索范围k_v的取值为k_v∈[0,V]△k_v，△k_v为调频率单元，v＝1,2,…,V，V为调频率的搜索个数，构造的过完备原子库G为(U×V)×N的矩阵：

$G\left(g_r\right)=\left[\begin{array}{cccc}{g}_r\left(f_{1,}{k}_1\right)& g_r\left(f_{1,}{k}_2\right)& \mathrm{...}& g_r\left(f_{1,}{k}_v\right)\\ g_r\left(f_{2,}{k}_1\right)& g_r\left(f_{2,}{k}_2\right)& \mathrm{...}& g_r\left(f_{2,}{k}_v\right)\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ g_r\left(f_{u,}{k}_1\right)& g_r\left(f_{u,}{k}_2\right)& \mathrm{...}& g_r\left(f_{u,}{k}_v\right)\end{array}\right]$

G＝[g₁,g₂,…,g_N]^T；

步骤3：利用正交匹配追踪方法估计目标径向加速度和径向速度

本发明采用稀疏分解中的正交匹配追踪方法对信号进行参数提取，首先将过完备原子库中的原子与信号进行匹配程度比较，选择与信号最匹配的一组基g_r，将所选的原子利用Grant-Schmidt正交化方法进行正交化处理，然后将信号在这些正交原子构成的空间下进行分解，使其满足：

|<s,g_r(f_u,k_v)>|＝sup|<s,g_r>|

因此，信号可以分解为在最佳原子上的分量和残余两部分，即：

s＝<s,g_r>g_r+R¹s

其中，R¹s是用最佳原子对原信号进行最佳匹配后的残余；然后从原子库中将最匹配的这组基删掉，接下来对最佳匹配后的残余可以不断进行上面同样的分解过程，即：

Rⁱs＝<Rⁱs,g_ri>g_ri+Rⁱ⁺¹s

其中g_ri满足：

<Rⁱs,g_ri>＝sup<Rⁱs,g_r>

经过i次迭代后能量衰减程度△Rⁱ⁺¹为：△Rⁱ⁺¹＝||Rⁱs||²-||Rⁱ⁺¹s||²＝ω²||Rⁱs||²，其中，||Rⁱs||²和||Rⁱ⁺¹s||²分别是信号s的i阶和i+1阶残余能量；ω²＝λ²(Rⁱs)表示信号衰减的速率；当λ²(Rⁱs)≥λ²时，停止分解；假设经过L步分解后停止分解，信号被分解为：

$s=\underset{i=0}{\overset{L-1}{\&Sigma;}}+R^{i}s,g_{ri}>g_{ri}+R^{L}s$

用少量的原子L(相对于信号长度N而言，L＜＜N)表示信号的主要成分，即：

$s=\underset{i=0}{\overset{L-1}{\&Sigma;}}<R^{i}s,g_{ri}>g_{ri}$

经过上述对信号s(t)的稀疏分解，计算出原子库中匹配的原子，得到信号的稀疏解能量图i∈(0,U×V)，在i∈(0,U×V)范围搜索最大的峰值，找到最大峰值的坐标Am(i)，并在原子库G中查找此坐标的位置的g_r(f_u,k_v)，就可以得到目标的初始频率f_u和调频率k_v，i＝1,2,…,N，最后根据以下公式得到机动目标的径向加速度和径向速度估计值：

$a_{rk}^m=k_v\mathrm{\&lambda;}/2,v_{rk}^m=f_u\mathrm{\&lambda;}/2$

将得到的径向加速度和径向速度估计值送到雷达数据处理计算机执行以下步骤：

步骤4：建立目标运动和量测方程，并进行EKF滤波

考虑二坐标雷达数据处理问题，直角坐标系下目标运动模型可表示为：

X(k)＝Φ(k,k-1)X(k-1)+w(k-1)

其中，状态向量为：

$X\left(k\right)={\left[\begin{array}{cccccccc}{x}_k& {\stackrel{\&CenterDot;}{x}}_k& {\stackrel{\&CenterDot;\&CenterDot;}{x}}_k& {\stackrel{\&CenterDot;\&CenterDot;\&CenterDot;}{x}}_k& y_k& {\stackrel{\&CenterDot;}{y}}_k& {\stackrel{\&CenterDot;\&CenterDot;}{y}}_k& {\stackrel{\&CenterDot;\&CenterDot;\&CenterDot;}{y}}_k\end{array}\right]}^T$

式中，和分别表示k时刻目标x方向和y方向的位置、速度、加速度、加加速度。Φ为状态转移矩阵，w(k-1)为零均值方差为Q的高斯白噪声。

状态转移矩阵Φ为：

$\left\{\begin{array}{c}\mathrm{\&Phi;}(k,k-1)=diag({\mathrm{\&Phi;}}_{11}\left(k,k-1\right),{\mathrm{\&Phi;}}_{22}\left(k,k-1\right))\\ {\mathrm{\&Phi;}}_{11}(k,k-1)={\mathrm{\&Phi;}}_{22}(k,k-1)=\left[\begin{array}{cccc}1& T& \frac{T^2}{2}& p_1\\ 0& 1& T& q_1\\ 0& 0& 1& r_1\\ 0& 0& 0& s_1\end{array}\right]\end{array}\right.$

式中，T为采样周期，α为机动时间常数；

p₁＝(2-2αT+α²T²-2e^-αT)/(2α³)

q₁＝(e^-αT-1+αT)/α²

r₁＝(1-e^-αT)/α

s₁＝e^-αT

其离散时间过程噪声w具有协方差：

$Q=2{\mathrm{\&alpha;\&sigma;}}_a^2\left[\begin{array}{cccc}{q}_{11}& q_{12}& q_{13}& q_{14}\\ q_{21}& q_{22}& q_{23}& q_{24}\\ q_{31}& q_{23}& q_{33}& q_{34}\\ q_{41}& q_{42}& q_{43}& q_{44}\end{array}\right]$

式中，Q精确表达式为：

$q_{11}=\frac{1}{2{\mathrm{\&alpha;}}^7}\&lsqb;\frac{{\mathrm{\&alpha;}}^5{T}^5}{10}-\frac{{\mathrm{\&alpha;}}^4{T}^4}{2}+\frac{4{\mathrm{\&alpha;}}^3{T}^3}{3}-2{\mathrm{\&alpha;}}^2{T}^2+2\mathrm{\&alpha;}T-3+4e^{-\mathrm{\&alpha;}T}+2{\mathrm{\&alpha;}}^2{T}^2{e}^{-\mathrm{\&alpha;}T}-e^{-2\mathrm{\&alpha;}T}\&rsqb;$

$q_{12}=\frac{1}{2{\mathrm{\&alpha;}}^6}\&lsqb;1-2\mathrm{\&alpha;}T+2{\mathrm{\&alpha;}}^2{T}^2-{\mathrm{\&alpha;}}^3{T}^3+\frac{{\mathrm{\&alpha;}}^4{T}^4}{4}+e^{-2\mathrm{\&alpha;}T}+2\mathrm{\&alpha;}T-2e^{-\mathrm{\&alpha;}T}-{\mathrm{\&alpha;}}^2{T}^2{e}^{-\mathrm{\&alpha;}T}\&rsqb;$

$q_{13}=\frac{1}{2{\mathrm{\&alpha;}}^5}\&lsqb;2\mathrm{\&alpha;}T-{\mathrm{\&alpha;}}^2{T}^2-\frac{1}{3}{\mathrm{\&alpha;}}^3{T}^3-3-2e^{-2\mathrm{\&alpha;}T}+4e^{-\mathrm{\&alpha;}T}+{\mathrm{\&alpha;}}^2{T}^2{e}^{-\mathrm{\&alpha;}T}\&rsqb;$

$q_{14}=\frac{1}{2{\mathrm{\&alpha;}}^4}\&lsqb;1+e^{-2\mathrm{\&alpha;}T}-2e^{-\mathrm{\&alpha;}T}+4e^{-\mathrm{\&alpha;}T}-{\mathrm{\&alpha;}}^2{T}^2{e}^{-\mathrm{\&alpha;}T}\&rsqb;$

$q_{22}=\frac{1}{2{\mathrm{\&alpha;}}^5}\&lsqb;1-e^{-2\mathrm{\&alpha;}T}+2\mathrm{\&alpha;}T+\frac{2{\mathrm{\&alpha;}}^3{T}^3}{3}-2{\mathrm{\&alpha;}}^2{T}^2-4{\mathrm{\&alpha;Te}}^{-\mathrm{\&alpha;}T}\&rsqb;$

$q_{23}=\frac{1}{2{\mathrm{\&alpha;}}^4}\&lsqb;e^{-2\mathrm{\&alpha;}T}+1-2e^{-\mathrm{\&alpha;}T}+2{\mathrm{\&alpha;Te}}^{-\mathrm{\&alpha;}T}-2\mathrm{\&alpha;}T+{\mathrm{\&alpha;}}^2{T}^2\&rsqb;$

$q_{24}=\frac{1}{2{\mathrm{\&alpha;}}^3}\&lsqb;1-e^{-2\mathrm{\&alpha;}T}-2{\mathrm{\&alpha;Te}}^{-\mathrm{\&alpha;}T}\&rsqb;$

$q_{33}=\frac{1}{2{\mathrm{\&alpha;}}^3}\&lsqb;4e^{-\mathrm{\&alpha;}T}-3-e^{-2\mathrm{\&alpha;}T}+2\mathrm{\&alpha;}T\&rsqb;$

$q_{34}=\frac{1}{2{\mathrm{\&alpha;}}^2}\&lsqb;e^{-2\mathrm{\&alpha;}T}+1-2e^{-\mathrm{\&alpha;}T}\&rsqb;$

$q_{44}=\frac{1}{2\mathrm{\&alpha;}}\&lsqb;1-e^{-2\mathrm{\&alpha;}T}\&rsqb;$

普通雷达在极坐标系下可获得位置量测(距离和方位角)，利用正交匹配追踪方法得到的径向速度和径向加速度估计值，测量向量表示为：

$\begin{array}{c}{Z}_k^m={\left[\begin{array}{cccc}{r}_k^m& {\mathrm{\&theta;}}_k^m& v_{rk}^m& a_{rk}^m\end{array}\right]}^T\\ ={\left[\begin{array}{cccc}{r}_k& {\mathrm{\&theta;}}_k& v_{rk}& a_{rk}\end{array}\right]}^T+{\left[\begin{array}{cccc}{\tilde{r}}_k& {\widetilde{\mathrm{\&theta;}}}_k& {\tilde{v}}_{rk}& {\tilde{a}}_{rk}\end{array}\right]}^T\end{array}$

其中

$\left\{\begin{array}{c}{r}_k=\sqrt{x_k^2+y_k^2}\\ {\mathrm{\&theta;}}_k=arctan(y_k/x_k)\end{array}\right.$

式中，为k时刻从雷达信号处理阶段获得的测量向量，其中包括距离量测方位角量测径向速度估计值和径向加速度估计值r_k、θ_k、v_rk和a_rk分别为距离、方位角、径向速度和径向加速度的真值；和为相应的量测误差，假设它们都是均值为零的高斯白噪声序列且相互统计独立，方差分别为和

对r_k分别求一阶和二阶导数，得到径向速度和径向加速度表达式：

$\left\{\begin{array}{c}{v}_{rk}={\stackrel{\&CenterDot;}{r}}_k={\stackrel{\&CenterDot;}{x}}_k\frac{x_k}{\sqrt{x_k^2+y_k^2}}+{\stackrel{\&CenterDot;}{y}}_k\frac{y_k}{\sqrt{x_k^2+y_k^2}}\\ a_{rk}={\stackrel{\&CenterDot;\&CenterDot;}{r}}_k={\stackrel{\&CenterDot;\&CenterDot;}{x}}_k\frac{x_k}{\sqrt{x_k^2+y_k^2}}+{\stackrel{\&CenterDot;\&CenterDot;}{y}}_k\frac{y_k}{\sqrt{x_k^2+y_k^2}}+\frac{{({\stackrel{\&CenterDot;}{x}}_k{y}_k-{\stackrel{\&CenterDot;}{y}}_k{x}_k)}^2}{(x_k^2+y_k^2)\sqrt{x_k^2+y_k^2}}\end{array}\right.$

引入径向加速度和径向速度量测后，量测向量一步预测为：

Z(k|k-1)＝H(k)X(k|k-1)

其中，量测雅克比矩阵为：

$\begin{array}{c}H\left(k\right)=\&lsqb;\&dtri;h^{\&prime;}\left(X\right)\&rsqb;\\ ={\&lsqb;{\left[\begin{array}{cccccc}\&part;/\&part;x& \&part;/\&part;\stackrel{\&CenterDot;}{x}& \&part;/\&part;\stackrel{\&CenterDot;\&CenterDot;}{x}& \&part;/\&part;y& \&part;/\&part;\stackrel{\&CenterDot;}{y}& \&part;/\&part;\stackrel{\&CenterDot;\&CenterDot;}{y}\end{array}\right]}^T\left[\begin{array}{c}{r}_k\\ {\mathrm{\&theta;}}_k\\ v_{vk}\\ a_{rk}\end{array}\right]\&rsqb;}_{X=X\left(k\right|k-1)}\end{array}$

其中

$f_1\left(X\right){|}_{X=X\left(k\right)}=\sqrt{{x_k}^2+y_k^2}$

$f_2\left(X\right){|}_{X=X\left(k\right)}=arctan(y_k/x_k)$

$f_3\left(X\right){|}_{X=X\left(k\right)}={\stackrel{\&CenterDot;}{x}}_k\frac{x_k}{\sqrt{x_k^2+y_k^2}}+{\stackrel{\&CenterDot;}{y}}_k\frac{y_k}{\sqrt{x_k^2+y_k^2}}$

$f_4\left(X\right){|}_{X=X\left(k\right)}={\stackrel{\&CenterDot;\&CenterDot;}{x}}_k\frac{k_k}{\sqrt{x_k^2+y_k^2}}+{\stackrel{\&CenterDot;\&CenterDot;}{y}}_k\frac{x_k}{\sqrt{x_k^2+y_k^2}}+\frac{{({\stackrel{\&CenterDot;}{x}}_k{y}_k-{\stackrel{\&CenterDot;}{y}}_k{x}_k)}^2}{(x_k^2+y_k^2)\sqrt{x_k^2+y_k^2}}$

$H\left(k\right)={\left[\begin{array}{cccccc}\frac{\&part;f_1\left(X\right)}{\&part;x}& 0& 0& \frac{\&part;f_1\left(X\right)}{\&part;y}& 0& 0\\ \frac{\&part;f_2\left(X\right)}{\&part;x}& 0& 0& \frac{\&part;f_2\left(X\right)}{\&part;y}& 0& 0\\ \frac{\&part;f_3\left(X\right)}{\&part;x}& \frac{\&part;f_3\left(X\right)}{\&part;\stackrel{\&CenterDot;}{x}}& 0& \frac{\&part;f_3\left(X\right)}{\&part;y}& \frac{\&part;f_3\left(X\right)}{\&part;\stackrel{\&CenterDot;}{y}}& 0\\ \frac{\&part;f_4\left(X\right)}{\&part;x}& 0& \frac{\&part;f_4\left(X\right)}{\&part;\stackrel{\&CenterDot;\&CenterDot;}{x}}& \frac{\&part;f_4\left(X\right)}{\&part;y}& 0& \frac{\&part;f_4\left(X\right)}{\&part;\stackrel{\&CenterDot;\&CenterDot;}{y}}\end{array}\right]}_{X=X\left(k\right|k-1)}$

其转换测量协方差为：

$R_k^c=\left[\begin{array}{cccc}{\mathrm{\&sigma;}}_r^2& 0& 0& 0\\ 0& {\mathrm{\&sigma;}}_{\mathrm{\&theta;}}^2& 0& 0\\ 0& 0& {\mathrm{\&sigma;}}_v^2& 0\\ 0& 0& 0& {\mathrm{\&sigma;}}_a^2\end{array}\right]$

最后对目标进行EKF滤波：

(1)状态向量一步预测为：

X(k)＝Φ(k|k-1)X(k-1)+W(k-1)

(2)状态向量预测协方差为：

P(k|k-1)＝f×Φ(k|k-1)P(k-1|k-1)+Q(k-1)

(3)量测向量一步预测为：

Z(k|k-1)＝H(k)X(k|k-1)

(4)新息为：

In(k)＝Z(k)-Z(k|k-1)

(5)新息协方差为：

S(k)＝H(k)P(k|k-1)H′(k)+R(k)

(6)卡尔曼增益为：

K(k)＝P(k|k-1)H(k)′S^-1(k)

(7)状态向量一步更新为：

X(k|k)＝X(k|k-1)+K(k)In(k)

(8)状态向量协方差更新为：

P(k|k)＝P(k|k-1)-K(k)S(k)K′(k)

和背景技术相比，本发明的有益效果说明：本发明方法将信号处理阶段得到的径向加速度和径向速度估计值引入量测方程和状态方程中，提高了机动目标跟踪精度，并且改善了滤波后加速度和速度估计精度。

附图说明

附图1是本发明的基于RAV-Jerk模型的机动目标跟踪方法的整体流程图；

附图2是现有的基于Singer模型方法和本发明的基于RAV-Jerk模型的机动目标跟踪方法的位置误差曲线对比图；

附图3是现有的基于Singer模型方法和本发明的基于RAV-Jerk模型的机动目标跟踪方法的速度误差曲线对比图；

附图4是现有的基于Singer模型方法和本发明的基于RAV-Jerk模型的机动目标跟踪方法的加速度误差曲线对比图。

具体实施方式

下面结合附图对本发明的基于RAV-Jerk模型的机动目标跟踪方法进行详细描述。(参照说明书附图1)。

实施例条件：假设雷达波长λ＝8mm，雷达脉冲积累时间T＝10ms，雷达信号处理采样间隔T_s＝5×10^-5秒，目标起始状态：X(0)＝[120000m,-426m/s,0m/s²,2000m,0m/s,0m/s²]^T，目标运动过程历时90s，发生机动时刻及加速度大小如表1所示。雷达数据处理采样间隔为距离量测误差角度量测误差径向速度量测误差径向加速度量测误差进行100次蒙特卡洛仿真，均方根误差为：其中为每次仿真得到的估计值，d为真实值。

表1目标机动运动情况

目标发生机动时刻(s)	t＝31	t＝38	t＝49	t＝61	t＝65	t＝66	t＝81
								X方向加速度值(m/s2)	5	-8	10	0	-10	-5	5
Y方向加速度值(m/s2)	-10	18	-20	30	-8	0	-10

根据上述条件由建立目标信号模型为：

s(t)＝exp(j2π×f_dt+j2π×kt²)+w(t)

其中，w(t)为均值为0，方差为1的高斯白噪声；

将上述模拟信号送入雷达处理计算机中执行以下步骤：

步骤1：将雷达接收机接收到的信号s(t)通过采样器以采样间隔T_s进行采样，变为离散信号s(nT_s)，其中n表示采样点序号，将s(nT_s)送入雷达信号处理计算机；

在雷达信号处理计算机中执行以下步骤：

步骤2：初始化

建立零矩阵G((U×V)×N)；

起始频率搜索个数U设为200；

调频率的搜索个数V设为200；

多普勒单元△f_u设为0.005；

调频率单元△k_v设为0.0025；

设置迭代停止的条件为相干比阈值λ²(Rⁿs)为0.03；

步骤3：形成过完备原子库

(1)根据LFM信号回波s(t)的特点，建立原子g_r＝exp[j2π(f_un+k_vn²/2N)]，r＝1,2,…,N；

(2)设定搜索精度和范围，假设搜索范围f_u的取值为f_u∈[0,U]△f_u，u＝1,2,…,U，U为起始频率的搜索个数，△f_u为多普勒单元，搜索范围△k_v的取值为k_v∈[0,V]△k_v，△k_v为调频率单元，v＝1,2,…,V，V为调频率的搜索个数，构造的过完备原子库G为(U×V)×N的矩阵：

$G\left(g_r\right)=\left[\begin{array}{cccc}{g}_r\left(f_{1,}{k}_1\right)& g_r\left(f_{1,}{k}_2\right)& \mathrm{...}& g_r\left(f_{1,}{k}_v\right)\\ g_r\left(f_{2,}{k}_1\right)& g_r\left(f_{2,}{k}_2\right)& \mathrm{...}& g_r\left(f_{2,}{k}_v\right)\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ g_r\left(f_{u,}{k}_1\right)& g_r\left(f_{u,}{k}_2\right)& \mathrm{...}& g_r\left(f_{u,}{k}_v\right)\end{array}\right]$

G＝[g₁,g₂,…,g_N]^T；

步骤4：利用正交匹配追踪方法估计目标径向加速度和径向速度

本发明采用稀疏分解中的正交匹配追踪方法对信号进行参数提取，首先将过完备原子库中的原子与信号进行匹配程度比较，选择与信号最匹配的一组基g_r，将所选的原子利用Grant-Schmidt正交化方法进行正交化处理，然后将信号在这些正交原子构成的空间下进行分解，使其满足：

|<s,g_r(f_u,k_v)>|＝sup|<s,g_r>|

因此，信号可以分解为在最佳原子上的分量和残余两部分，即：

s＝<s,g_r>g_r+R¹s

其中，R¹s是用最佳原子对原信号进行最佳匹配后的残余；然后从原子库中将最匹配的这组基删掉，接下来对最佳匹配后的残余可以不断进行上面同样的分解过程，即：

Rⁱs＝<Rⁱs,g_ri>g_ri+Rⁱ⁺¹s

其中g_ri满足：

<Rⁱs,g_ri>＝sup<Rⁱs,g_r>

经过i次迭代后能量衰减程度△Rⁱ⁺¹为：△Rⁱ⁺¹＝||Rⁱs||²-||Rⁱ⁺¹s||²＝ω²||Rⁱs||²，其中，||Rⁱs||²和||Rⁱ⁺¹s||²分别是信号s的i阶和i+1阶残余能量；ω²＝λ²(Rⁱs)表示信号衰减的速率；当λ²(Rⁱs)≥λ²时，停止分解；假设经过L步分解后停止分解，信号被分解为：

$s=\underset{i=0}{\overset{L-1}{\&Sigma;}}<R^{L}s,g_{ri}>g_{ri}+R^{L}s$

用少量的原子L(相对于信号长度N而言，L＜＜N)表示信号的主要成分，即：

$s=\underset{i=0}{\overset{L-1}{\&Sigma;}}<R^{i}s,g_{ri}>g_{ri}$

步骤5：径向加速度和径向速度的确定

(1)经过上述对信号s(t)的稀疏分解，计算出原子库中匹配的原子，得到信号的稀疏解能量图i∈(0,U×V)；

(2)在i∈(0,U×V)范围搜索最大的峰值；

(3)找到最大峰值的坐标Am(i)，并在原子库G中查找此坐标的位置的g_r(f_u,k_v)，就可以得到目标的初始频率f_u和调频率k_v，i＝1,2,…,N，最后根据以下公式得到机动目标的径向加速度和径向速度估计值：

$a_{rk}^m=k_v\mathrm{\&lambda;}/2,v_{rk}^m=f_u\mathrm{\&lambda;}/2$

将得到的径向加速度和径向速度估计值送到雷达数据处理计算机执行以下步骤：

步骤6：建立目标运动方程

考虑二坐标雷达数据处理问题，直角坐标系下目标运动模型可表示为：

X(k)＝Φ(k,k-1)X(k-1)+w(k-1)

其中，状态向量为：

$X\left(k\right)={\left[\begin{array}{cccccccc}{x}_k& {\stackrel{\&CenterDot;}{x}}_k& {\stackrel{\&CenterDot;\&CenterDot;}{x}}_k& {\stackrel{\&CenterDot;\&CenterDot;\&CenterDot;}{x}}_k& y_k& {\stackrel{\&CenterDot;}{y}}_k& {\stackrel{\&CenterDot;\&CenterDot;}{y}}_k& {\stackrel{\&CenterDot;\&CenterDot;\&CenterDot;}{y}}_k\end{array}\right]}^T$

式中，x_k、和分别表示k时刻目标x方向和y方向的位置、速度、加速度、加加速度。Φ为状态转移矩阵，w(k-1)为零均值方差为Q的高斯白噪声。

状态转移矩阵Φ为：

$\left\{\begin{array}{c}\mathrm{\&Phi;}(k,k-1)=diag({\mathrm{\&Phi;}}_{11}\left(k,k-1\right),{\mathrm{\&Phi;}}_{22}\left(k,k-1\right))\\ {\mathrm{\&Phi;}}_{11}(k,k-1)={\mathrm{\&Phi;}}_{22}(k,k-1)=\left[\begin{array}{cccc}1& T& \frac{T^2}{2}& p_1\\ 0& 1& T& q_1\\ 0& 0& 1& r_1\\ 0& 0& 0& s_1\end{array}\right]\end{array}\right.$

式中，T为采样周期，α为机动时间常数；

p₁＝(2-2αT+α²T²-2e^-αT)/(2α³)

q₁＝(e^-αT-1+αT)/α²

r₁＝(1-e^-αT)/α

s₁＝e^-αT

其离散时间过程噪声w具有协方差：

$Q=2{\mathrm{\&alpha;\&sigma;}}_a^2\left[\begin{array}{cccc}{q}_{11}& q_{12}& q_{13}& q_{14}\\ q_{21}& q_{22}& q_{23}& q_{24}\\ q_{31}& q_{23}& q_{33}& q_{34}\\ q_{41}& q_{42}& q_{43}& q_{44}\end{array}\right]$

式中，Q精确表达式为：

$q_{11}=\frac{1}{2{\mathrm{\&alpha;}}^7}\&lsqb;\frac{{\mathrm{\&alpha;}}^5{T}^5}{10}-\frac{{\mathrm{\&alpha;}}^4{T}^4}{2}+\frac{4{\mathrm{\&alpha;}}^3{T}^3}{3}-2{\mathrm{\&alpha;}}^2{T}^2+2\mathrm{\&alpha;}T-3+4e^{-\mathrm{\&alpha;}T}+2{\mathrm{\&alpha;}}^2{T}^2{e}^{-\mathrm{\&alpha;}T}-e^{-2\mathrm{\&alpha;}T}\&rsqb;$

$q_{12}=\frac{1}{2{\mathrm{\&alpha;}}^6}\&lsqb;1-2\mathrm{\&alpha;}T+2{\mathrm{\&alpha;}}^2{T}^2-{\mathrm{\&alpha;}}^3{T}^3+\frac{{\mathrm{\&alpha;}}^4{T}^4}{4}+e^{-2\mathrm{\&alpha;}T}+2\mathrm{\&alpha;}T-2e^{-\mathrm{\&alpha;}T}-{\mathrm{\&alpha;}}^2{T}^2{e}^{-\mathrm{\&alpha;}T}\&rsqb;$

$q_{13}=\frac{1}{2{\mathrm{\&alpha;}}^5}\&lsqb;2\mathrm{\&alpha;}T-{\mathrm{\&alpha;}}^2{T}^2-\frac{1}{3}{\mathrm{\&alpha;}}^3{T}^3-3-2e^{-2\mathrm{\&alpha;}T}+4e^{-\mathrm{\&alpha;}T}+{\mathrm{\&alpha;}}^2{T}^2{e}^{-\mathrm{\&alpha;}T}\&rsqb;$

$q_{14}=\frac{1}{2{\mathrm{\&alpha;}}^4}\&lsqb;1+e^{-2\mathrm{\&alpha;}T}-2e^{-\mathrm{\&alpha;}T}+4e^{-\mathrm{\&alpha;}T}-{\mathrm{\&alpha;}}^2{T}^2{e}^{-\mathrm{\&alpha;}T}\&rsqb;$

$q_{22}=\frac{1}{2{\mathrm{\&alpha;}}^5}\&lsqb;1-e^{-2\mathrm{\&alpha;}T}+2\mathrm{\&alpha;}T+\frac{2{\mathrm{\&alpha;}}^3{T}^3}{3}-2{\mathrm{\&alpha;}}^2{T}^2-4{\mathrm{\&alpha;Te}}^{-\mathrm{\&alpha;}T}\&rsqb;$

$q_{23}=\frac{1}{2{\mathrm{\&alpha;}}^4}\&lsqb;e^{-2\mathrm{\&alpha;}T}+1-2e^{-\mathrm{\&alpha;}T}+2{\mathrm{\&alpha;Te}}^{-\mathrm{\&alpha;}T}-2\mathrm{\&alpha;}T+{\mathrm{\&alpha;}}^2{T}^2\&rsqb;$

$q_{24}=\frac{1}{2{\mathrm{\&alpha;}}^3}\&lsqb;1-e^{-2\mathrm{\&alpha;}T}-2{\mathrm{\&alpha;Te}}^{-\mathrm{\&alpha;}T}\&rsqb;$

$q_{33}=\frac{1}{2{\mathrm{\&alpha;}}^3}\&lsqb;4e^{-\mathrm{\&alpha;}T}-3-e^{-2\mathrm{\&alpha;}T}+2\mathrm{\&alpha;}T\&rsqb;$

$q_{34}=\frac{1}{2{\mathrm{\&alpha;}}^2}\&lsqb;e^{-2\mathrm{\&alpha;}T}+1-2e^{-\mathrm{\&alpha;}T}\&rsqb;$

$q_{44}=\frac{1}{2\mathrm{\&alpha;}}\&lsqb;1-e^{-2\mathrm{\&alpha;}T}\&rsqb;$

步骤7：确定量测方程：

雷达在极坐标系下可获得位置量测(距离和方位角)，利用正交匹配追踪方法得到的径向速度和径向加速度估计值，测量向量可表示为：

$\begin{array}{c}{Z}_k^m={\left[\begin{array}{cccc}{r}_k^m& {\mathrm{\&theta;}}_k^m& v_{rk}^m& a_{rk}^m\end{array}\right]}^T\\ ={\left[\begin{array}{cccc}{r}_k& {\mathrm{\&theta;}}_k& v_{rk}& a_{rk}\end{array}\right]}^T+{\left[\begin{array}{cccc}{\tilde{r}}_k& {\widetilde{\mathrm{\&theta;}}}_k& {\tilde{v}}_{rk}& {\tilde{a}}_{rk}\end{array}\right]}^T\end{array}$

其中

$\left\{\begin{array}{c}{r}_k=\sqrt{x_k^2+y_k^2}\\ {\mathrm{\&theta;}}_k=arctan(y_k/x_k)\end{array}\right.$

式中，为k时刻从雷达信号处理阶段获得的测量向量，其中包括距离量测方位角量测径向速度估计值和径向加速度估计值r_k、θ_k、v_rk和a_rk分别为距离、方位角、径向速度和径向加速度的真值。和为相应的量测误差，假设它们都是均值为零的高斯白噪声序列且相互统计独立，方差分别为和

对r_k分别求一阶和二阶导数，得到径向速度和径向加速度表达式：

$\left\{\begin{array}{c}{v}_{rk}={\stackrel{\&CenterDot;}{r}}_k={\stackrel{\&CenterDot;}{x}}_k\frac{x_k}{\sqrt{x_k^2+y_k^2}}+{\stackrel{\&CenterDot;}{y}}_k\frac{y_k}{\sqrt{x_k^2+y_k^2}}\\ a_{rk}={\stackrel{\&CenterDot;\&CenterDot;}{r}}_k={\stackrel{\&CenterDot;\&CenterDot;}{x}}_k\frac{x_k}{\sqrt{x_k^2+y_k^2}}+{\stackrel{\&CenterDot;\&CenterDot;}{y}}_k\frac{y_k}{\sqrt{x_k^2+y_k^2}}+\frac{{({\stackrel{\&CenterDot;}{x}}_k{y}_k-{\stackrel{\&CenterDot;}{y}}_k{x}_k)}^2}{(x_k^2+y_k^2)\sqrt{x_k^2+y_k^2}}\end{array}\right.$

引入径向加速度和径向速度估计值后，量测向量一步预测为：

Z(k|k-1)＝H(k)X(k|k-1)

其中，量测雅克比矩阵为：

$\begin{array}{c}H\left(k\right)=\&lsqb;\&dtri;h^{\&prime;}\left(X\right)\&rsqb;\\ ={\&lsqb;{\left[\begin{array}{cccccc}\&part;/\&part;x& \&part;/\&part;\stackrel{\&CenterDot;}{x}& \&part;/\&part;\stackrel{\&CenterDot;\&CenterDot;}{x}& \&part;/\&part;y& \&part;/\&part;\stackrel{\&CenterDot;}{y}& \&part;/\&part;\stackrel{\&CenterDot;\&CenterDot;}{y}\end{array}\right]}^T\left[\begin{array}{c}{r}_k\\ {\mathrm{\&theta;}}_k\\ v_{vk}\\ a_{rk}\end{array}\right]\&rsqb;}_{X=X\left(k\right|k-1)}\end{array}$

其中

$f_1\left(X\right){|}_{X=X\left(k\right)}=\sqrt{{x_k}^2+y_k^2}$

f₂(X)|_X＝X(k)＝arctan(y_k/x_k)

$f_3\left(X\right){|}_{X=X\left(k\right)}={\stackrel{\&CenterDot;}{x}}_k\frac{x_k}{\sqrt{x_k^2+y_k^2}}+{\stackrel{\&CenterDot;}{y}}_k\frac{y_k}{\sqrt{x_k^2+y_k^2}}$

$f_4\left(X\right){|}_{X=X\left(k\right)}={\stackrel{\&CenterDot;\&CenterDot;}{x}}_k\frac{x_k}{\sqrt{x_k^2+y_k^2}}+{\stackrel{\&CenterDot;\&CenterDot;}{y}}_k\frac{x_k}{\sqrt{x_k^2+y_k^2}}+\frac{{({\stackrel{\&CenterDot;}{x}}_k{y}_k-{\stackrel{\&CenterDot;}{y}}_k{x}_k)}^2}{(x_k^2+y_k^2)\sqrt{x_k^2+y_k^2}}$

$H\left(k\right)={\left[\begin{array}{cccccc}\frac{\&part;f_1\left(X\right)}{\&part;x}& 0& 0& \frac{\&part;f_1\left(X\right)}{\&part;y}& 0& 0\\ \frac{\&part;f_2\left(X\right)}{\&part;x}& 0& 0& \frac{\&part;f_2\left(X\right)}{\&part;y}& 0& 0\\ \frac{\&part;f_3\left(X\right)}{\&part;x}& \frac{\&part;f_3\left(X\right)}{\&part;\stackrel{\&CenterDot;}{x}}& 0& \frac{\&part;f_3\left(X\right)}{\&part;y}& \frac{\&part;f_3\left(X\right)}{\&part;\stackrel{\&CenterDot;}{y}}& 0\\ \frac{\&part;f_4\left(X\right)}{\&part;x}& 0& \frac{\&part;f_4\left(X\right)}{\&part;\stackrel{\&CenterDot;\&CenterDot;}{x}}& \frac{\&part;f_4\left(X\right)}{\&part;y}& 0& \frac{\&part;f_4\left(X\right)}{\&part;\stackrel{\&CenterDot;\&CenterDot;}{y}}\end{array}\right]}_{X=X\left(k\right|k-1)}$

其转换测量协方差为：

$R_k^c=\left[\begin{array}{cccc}{\mathrm{\&sigma;}}_r^2& 0& 0& 0\\ 0& {\mathrm{\&sigma;}}_{\mathrm{\&theta;}}^2& 0& 0\\ 0& 0& {\mathrm{\&sigma;}}_v^2& 0\\ 0& 0& 0& {\mathrm{\&sigma;}}_a^2\end{array}\right]$

步骤8：算法流程：

1)状态向量一步预测为：

X(k)＝Φ(k|k-1)X(k-1)+W(k-1)

2)状态向量预测协方差为：

P(k|k-1)＝f×Φ(k|k-1)P(k-1|k-1)+Q(k-1)

3)量测向量一步预测为：

Z(k|k-1)＝H(k)X(k|k-1)

4)新息为：

In(k)＝Z(k)-Z(k|k-1)

5)新息协方差为：

S(k)＝H(k)P(k|k-1)H′(k)+R(k)

6)卡尔曼增益为：

K(k)＝P(k|k-1)H(k)′S^-1(k)

7)状态向量一步更新为：

X(k|k)＝X(k|k-1)+K(k)In(k)

8)状态向量协方差更新为：

P(k|k)＝P(k|k-1)-K(k)S(k)K′(k)

为了比较对现有技术的基于Singer模型方法在相同条件下进行了仿真，附图2是现有的基于Singer模型方法和本发明的基于RAV-Jerk模型的机动目标跟踪方法的位置误差曲线对比图，附图3是现有的基于Singer模型方法和本发明的基于RAV-Jerk模型的机动目标跟踪方法的速度误差曲线对比图，附图4是现有的基于Singer模型方法和本发明的基于RAV-Jerk模型的机动目标跟踪方法的加速度误差曲线对比图。从仿真结果可以看出：本发明和背景技术相比，目标跟踪精度有了明显的提高，同时滤波后速度和加速度估计精度也有了很大改善。

Claims (4)

1.一种基于RAV-Jerk模型的机动目标跟踪方法，其特征在于包括以下步骤：

步骤1，将雷达接收机接收到的线性调频信号s(t)进行采样，送入雷达信号处理计算机；

步骤2，根据机动目标信号回波s(t)的特点，建立与之相匹配的原子g_r，构造出过完备原子库G；

步骤3，利用正交匹配追踪方法估计目标径向加速度和径向速度，送入雷达数据处理计算机；

步骤4，建立目标运动和量测方程，并进行EKF滤波。

2.根据权利要求1所述的基于RAV-Jerk模型的机动目标跟踪方法，其特征在于，步骤2具体为：

步骤2-1，根据机动目标回波具有LFM信号的特点，建立原子g_r＝exp[j2π(f_un+k_vn²/2N)]；

步骤2-2，设定搜索f_u和k_v精度和范围，构造的过完备原子库G为(U×V)×N的矩阵：

$G\left(g_r\right)=\left[\begin{array}{cccc}{g}_r\left(f_{1,}{k}_1\right)& g_r\left(f_{1,}{k}_2\right)& ...& g_r\left(f_{1,}{k}_v\right)\\ g_r\left(f_{2,}{k}_1\right)& g_r\left(f_{2,}{k}_2\right)& ...& g_r\left(f_{2,}{k}_v\right)\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ g_r\left(f_{u,}{k}_1\right)& g_r\left(f_{u,}{k}_2\right)& ...& g_r\left(f_{u,}{k}_v\right)\end{array}\right]$

G＝[g₁,g₂,…,g_N]^T。

3.根据权利要求1所述的基于RAV-Jerk模型的机动目标跟踪方法，其特征在于，步骤3具体为：

步骤3-1，利用正交匹配追踪将信号s(t)稀疏分解后，得到信号的稀疏解能量图i∈(0,U×V)，在i∈(0,U×V)范围搜索最大的峰值，找到最大峰值的坐标Am(i)，并在原子库G中查找此坐标的位置的g_r(f_u,k_v)，就可以得到目标的初始频率f_u和调频率k_v，i＝1,2,…,N；

步骤3-2，根据以下公式得到机动目标的径向加速度和径向速度估计值：

$a_{rk}^m=k_v\mathrm{\&lambda;}/2,v_{rk}^m=f_u\mathrm{\&lambda;}/2.$

4.据权利要求1所述的基于RAV-Jerk模型的机动目标跟踪方法，其特征在于，步骤4具体为：步骤4-1，在数据处理阶段引入正交匹配追踪方法得到的径向速度和径向加速度估计值测量向量表示为：

$\begin{array}{c}{Z}_k^m={\left[\begin{array}{cccc}{r}_k^m& {\mathrm{\&theta;}}_k^m& v_{rk}^m& a_{rk}^m\end{array}\right]}^T\\ ={\left[\begin{array}{cccc}{r}_k& {\mathrm{\&theta;}}_k& v_{rk}& a_{rk}\end{array}\right]}^T+{\left[\begin{array}{cccc}{\tilde{r}}_k& {\widetilde{\mathrm{\&theta;}}}_k& {\tilde{v}}_{rk}& {\tilde{a}}_{rk}\end{array}\right]}^T\end{array}$

式中，为k时刻从雷达信号处理阶段获得的测量向量，其中包括距离量测方位角量测径向速度估计值和径向加速度估计值r_k、θ_k、v_rk和a_rk分别为距离、方位角、径向速度和径向加速度的真值，和为相应的量测误差；

步骤4-2，对r_k分别求一阶和二阶导数，得到径向速度和径向加速度表达式：

$\left\{\begin{array}{c}{v}_{rk}={\stackrel{\&CenterDot;}{r}}_k={\stackrel{\&CenterDot;}{x}}_k\frac{x_k}{\sqrt{x_k^2+y_k^2}}+{\stackrel{\&CenterDot;}{y}}_k\frac{y_k}{\sqrt{x_k^2+y_k^2}}\\ a_{rk}={\stackrel{\&CenterDot;\&CenterDot;}{r}}_k={\stackrel{\&CenterDot;\&CenterDot;}{x}}_k\frac{x_k}{\sqrt{x_k^2+y_k^2}}+{\stackrel{\&CenterDot;\&CenterDot;}{y}}_k\frac{y_k}{\sqrt{x_k^2+y_k^2}}+\frac{{({\stackrel{\&CenterDot;}{x}}_k{y}_k-{\stackrel{\&CenterDot;}{y}}_k{x}_k)}^2}{(x_k^2+y_k^2)\sqrt{x_k^2+y_k^2}}\end{array}\right.$

步骤4-3，量测向量一步预测表示为：

Z(k|k-1)＝H(k)X(k|k-1)

其中，量测雅克比矩阵为：

$\begin{array}{c}H\left(k\right)=\&lsqb;\&dtri;h^{\&prime;}\left(X\right)\&rsqb;\\ ={\&lsqb;{\left[\begin{array}{cccccc}\&part;/\&part;x& \&part;/\&part;\stackrel{\&CenterDot;}{x}& \&part;/\&part;\stackrel{\&CenterDot;\&CenterDot;}{x}& \&part;/\&part;y& \&part;/\&part;\stackrel{\&CenterDot;}{y}& \&part;/\&part;\stackrel{\&CenterDot;\&CenterDot;}{y}\end{array}\right]}^T\left[\begin{array}{c}{r}_k\\ {\mathrm{\&theta;}}_k\\ v_{rk}\\ a_{rk}\end{array}\right]\&rsqb;}_{X=X\left(k\right|k-1)}\end{array}$

其转换测量协方差为：

$R_k^c=\left[\begin{array}{cccc}{\mathrm{\&sigma;}}_r^2& 0& 0& 0\\ 0& {\mathrm{\&sigma;}}_{\mathrm{\&theta;}}^2& 0& 0\\ 0& 0& {\mathrm{\&sigma;}}_v^2& 0\\ 0& 0& 0& {\mathrm{\&sigma;}}_a^2\end{array}\right]$

和分别为量测误差的方差。