CN103048658A - 基于径向加速度的RA-Singer-EKF机动目标跟踪算法 - Google Patents

Abstract

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

Description

基于径向加速度的RA-Singer-EKF机动目标跟踪算法

技术领域

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

背景技术

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

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

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

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

这种方法具有以下缺陷：

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

发明内容

本发明的目的是提出一种基于径向加速度的RA-Singer-EKF机动目标跟踪算法，解决现有Singer模型算法对目标突发机动的自适应跟踪能力不强，以及速度和加速度估计误差较大的问题。

本发明提出的基于径向加速度的RA-Singer-EKF方法的技术方案包括以下步骤：步骤1：将雷达接收机接收到的线性调频信号s(t)通过采样器以采样间隔T_s进行采样，变为离散信号s(nT_s)，其中n表示采样点序号；将s(nT_s)送入雷达信号处理计算机；

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

步骤2：初始化(设置分解参数)

T设为雷达的脉冲宽度；

λ为雷达波长；

f_s为采样频率；

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

k_v，设为LFM信号的调频率；

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

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

步骤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)& ...& 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；

步骤4：正交匹配追踪(OMP)

本发明采用正交匹配追踪(OMP)方法对信号进行参数提取，和标准匹配追踪(MP)方法相比，该方法在信号参数估计时精度更高，并且稀疏分解收敛速度更快；首先将过完备原子库中的原子与信号进行匹配程度比较，选择与信号最匹配的一组基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}{\mathrm{\&Sigma;}}}<R^{i}s,g_{\mathrm{ri}}>g_{\mathrm{ri}}+R^{L}s$

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

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

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

(1)经过上述对信号s(t)的OMP分解，计算出原子库中匹配的原子，得到信号的稀疏解能量图

i∈(0，U×V)；

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

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

α＝kλ/2，v＝f_dλ/2

步骤6：将得到的径向加速度和径向速度估计值送到雷达数据处理计算机，进行极(球)坐标系到直角坐标系的转换

根据速度与状态向量几何关系图可以得到径向速度和直角坐标系下X轴方向速度v_x和Y轴方向速度v_y的关系：

$\left\{\begin{array}{c}{v}_{\mathrm{\&rho;}}=v_{\mathrm{\&rho;x}}+v_{\mathrm{\&rho;y}}\\ v_{\mathrm{\&rho;x}}=v_x(x/\sqrt{x^2+y^2})\\ v_{\mathrm{\&rho;y}}=v_y(y/\sqrt{x^2+y^2})\end{array}\right.---\left(1\right)$

同理可得径向加速度和直角坐标系下X轴方向速度α_x和Y轴方向速度α_y的关系：

$\left\{\begin{array}{c}{a}_{\mathrm{\&rho;}}=a_{\mathrm{\&rho;x}}+a_{\mathrm{\&rho;y}}\\ a_{\mathrm{\&rho;x}}=a_x(x/\sqrt{x^2+y^2})\\ a_{\mathrm{\&rho;y}}=a_y(y/\sqrt{x^2+y^2})\end{array}\right.---\left(2\right)$

步骤7：建立Singer模型

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

$R\left(\mathrm{\&tau;}\right)=E\left[a\left(t\right)a(t+\mathrm{\&tau;})\right]={\mathrm{\&sigma;}}_a^2{e}^{-\mathrm{\&alpha;}\left|\mathrm{\&tau;}\right|}$

式中，

α是在区间(t，t+τ)内决定目标机动特性参数的待定参数；α为机动频率，即机动时间常数的倒数，其取值与具体机动状态有关；

为目标的加速度方差，可由α(t)概率密度模型求取。

对时间相关函数R(τ)应用Wiener-Kolmogorov白化程序后，机动加速度可用输入为白噪声的一阶时间相关模型来表征，即

$\stackrel{\&CenterDot;}{a}\left(t\right)=-\mathrm{\&alpha;a}\left(t\right)+\sqrt{2\mathrm{\&alpha;}{\mathrm{\&sigma;}}_a^2}u\left(t\right)=-\mathrm{\&alpha;a}\left(t\right)+{\mathrm{\&omega;}}_1\left(t\right)$

式中u(t)表示均值为零、方差为1的高斯白噪声；ω₁(t)表示均值为零、方差为

的高斯白噪声。

步骤8：确定量测方程：

设状态向量X表示为

X＝[x v_x α_x y v_y α_y]^T    (3)

式中x、v_x、α_x和y、v_y、α_y和分别表示直角坐标系下X轴和Y轴的位置、速度和加速度。

设量测向量Z表示为：

$Z={\left[\begin{array}{cccc}{z}_x^d& z_y^d& z_{\mathrm{rv}}& z_{\mathrm{ra}}\end{array}\right]}^T---\left(4\right)$

式中

和分别表示极-直坐标的无偏转换后直角坐标系下X轴和Y轴的位置；z_rv、z_ra表示信号处理阶段得到的径向速度和径向加速度估计值，设其量测误差

均为零均值高斯白噪声且其协方差为

由(1)、(2)、(3)、(4)式可得状态向量到量测向量的方程：

Z(k+1|k)H[X(k+1|k)]+V(K+1)                         (5)

式中H[X(k+1|k)]表示量测矩阵，V(K+1)表示方差为R(k+1)的量测噪声。

其中

$H(K+1)={\left[\begin{array}{c}x\\ y\\ v_x(x/\sqrt{x^2+y^2})+v_y(y/\sqrt{x^2+y^2})\\ a_x(x/\sqrt{x^2+y^2})+a_y(y/\sqrt{x^2+y^2})\end{array}\right]}_{X=X(k+1|k)}$

$R(k+1)=\mathrm{blkdig}(R^d(k+1),{\mathrm{\&sigma;}}_{\mathrm{rv}}^2,{\mathrm{\&sigma;}}_{\mathrm{ra}}^2)$

雅克比矩阵为：

$H_X(K+1)=\left[{\&dtri;}_X{h}^{\&prime;}\left(X\right)\right]={\left[{\left[\begin{array}{cccccc}\&PartialD;/\&PartialD;x& \&PartialD;/\&PartialD;v_x& \&PartialD;/\&PartialD;a_x& \&PartialD;/\&PartialD;y& \&PartialD;/\&PartialD;v_y& \&PartialD;/\&PartialD;a_y\end{array}\right]}^T\begin{array}{c}\end{array}\left[\begin{array}{c}x\\ y\\ f_1\left(X\right)\\ f_2\left(X\right)\end{array}\right]\right]}_{X=X(k+1|k)}$

式中

$f_1\left(X\right)=v_x(x/\sqrt{x^2+y^2})+v_y{(y/\sqrt{x^2+y^2})}_{X=X(k+1|k)},$

$f_2\left(X\right)=a_x(x/\sqrt{x^2+y^2})+a_y{(y/\sqrt{x^2+y^2})}_{X=X(k+1|k)};$

步骤9：确定状态方程

假设连续时间状态方程可以表示为

$\frac{d}{\mathrm{dt}}X\left(t\right)=\mathrm{AX}\left(t\right)+B\stackrel{\&OverBar;}{a}\left(t\right)+\mathrm{c\&omega;}\left(t\right)$

式中状态向量X(t)＝[x(t) v_x(t) v_y(t) y(t) α_x(t) α_y(t)]^T分别表示目标X轴、Y轴方向位置、速度加速度； $A=\left[\begin{array}{cccccc}0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& -\mathrm{\&alpha;}& 0& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& -\mathrm{\&alpha;}\end{array}\right],$ $B=\left[\begin{array}{cc}0& 0\\ 0& 0\\ \mathrm{\&alpha;}& 0\\ 0& 0\\ 0& 0\\ 0& \mathrm{\&alpha;}\end{array}\right].$

设采样周期为T′_s，对连续时刻的状态方程式进行离散化处理，可得到下列离散状态方程：

$X(k+1)=\mathrm{\&Phi;}(k+1,k)X\left(k\right|k)+U\left(k\right)\stackrel{\&OverBar;}{a}\left(k\right)+w\left(k\right)---\left(6\right)$

式中X(k|k)＝[x(k|k) v_x(k|k) α_x(k|k) y(k|k) v_y(k|k) α_y(k|k)]^T，

$\mathrm{\&Phi;}(k+1,k)=e^{\mathrm{AT}}=\left[\begin{array}{cccccc}0& T& (-1+\mathrm{\&alpha;T}+e^{-\mathrm{\&alpha;T}})/{\mathrm{\&alpha;}}^2& 0& 0& 0\\ 0& 1& (1-e^{-\mathrm{\&alpha;T}})/\mathrm{\&alpha;}& 0& 0& 0\\ 0& 0& e^{-\mathrm{\&alpha;T}}& 0& 0& 0\\ 0& 0& 0& 0& T& (-1+\mathrm{\&alpha;T}+e^{-\mathrm{\&alpha;T}})/{\mathrm{\&alpha;}}^2\\ 0& 0& 0& 0& 1& (1-e^{-\mathrm{\&alpha;T}})/\mathrm{\&alpha;}\\ 0& 0& 0& 0& 0& e^{-\mathrm{\&alpha;T}}\end{array}\right],$

$U\left(k\right)={\&Integral;}_{\mathrm{kT}}^{(k+1)T}{e}^{A[(k+1)T-\mathrm{\&tau;}]}\mathrm{Bd\&tau;}={\&Integral;}_0^T{e}^{\mathrm{At}}\mathrm{Bdt}=\left[\begin{array}{cc}(-T+\mathrm{\&alpha;}{T}^2/2+(1-e^{-\mathrm{\&alpha;T}})/\mathrm{\&alpha;})/\mathrm{\&alpha;}& 0\\ T-\frac{1-e^{-\mathrm{\&alpha;T}}}{\mathrm{\&alpha;}}& 0\\ 1-e^{-\mathrm{\&alpha;T}}& 0\\ 0& (-T+\mathrm{\&alpha;}{T}^2/2+(1-e^{-\mathrm{\&alpha;T}})/\mathrm{\&alpha;})/\mathrm{\&alpha;}\\ 0& T-\frac{1-e^{-\mathrm{\&alpha;T}}}{\mathrm{\&alpha;}}\\ 0& 1-e^{-\mathrm{\&alpha;T}}\end{array}\right],$

$\stackrel{\&OverBar;}{a}\left(k\right)=\left[\frac{4-\mathrm{\&pi;}}{\mathrm{\&pi;}}{(a}_{\max }-a_x\left(k\right|k),\frac{4-\mathrm{\&pi;}}{\mathrm{\&pi;}}(a_{\max }-a_y\left(k\right|k)\right],$

$w\left(k\right)={\&Integral;}_{\mathrm{kT}}^{(k+1)T}{e}^{A[(k+1)T-\mathrm{\&tau;}]}\mathrm{C\&omega;}\left(\mathrm{\&tau;}\right)\mathrm{d\&tau;}={\&Integral;}_0^T\left[\begin{array}{c}\frac{1-e^{-\mathrm{\&alpha;T}}}{\mathrm{\&alpha;}}\\ e^{-\mathrm{\&alpha;T}}\end{array}\right]\mathrm{\&omega;}(t_{k+1}-t)\mathrm{dt},$

Q(k)＝E[w_j(k)·w_j ^T(k)]；

步骤10：EKF算法：

采用量测方程(5)和状态方程(6)的EKF算法流程为：

$X(k+1|k)=A\left(k\right)X\left(k\right|k)+C\left(k\right)\stackrel{\&OverBar;}{a}\left(k\right)$

式中 $\stackrel{\&OverBar;}{a}\left(k\right)=\mathrm{diag}(a_x\left(k\right|k),a_y\left(k\right|k)).$

P(k+1|k)＝A(k)P(k)A^T(k)+Q(k)

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

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

$\mathrm{Kg}(k+1)=P(k+1|k)H_X^T(k+1){[H_X(k+1)P(k+1|k)H_X^T(k+1)+R(K+1)]}^{-1}$

X(k+1|k+1)＝X(k+1|k)+Kg(k+1)In(k+1)

P(k+1|k+1)＝[I-Kg(k+1)H_X(k+1)]P(k+1|k)[I+Kg(k+1)H_X(k+1)]^T

-Kg(k+1)R(k+1)Kg^T(k+1)

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

附图说明

附图1是径向速度与状态向量几何关系图，s表示t时刻目标位置，向量

分别表示直角坐标系下X轴、Y轴速度方向，v_x、v_y分别表示X轴、Y轴方向速度大小；向量

表示极坐标系下径向速度、切向速度方向，v_ρ表示t时刻径向速度大小，v_ρx、v_ρy分别表式v_x、v_y在径向速度方向的分量大小；

附图2是本发明的基于径向加速度的RA-Singer-EKF方法的整体流程图；

附图3是现有的基于Singer模型方法和本发明的基于径向加速度的RA-Singer-EKF方法的位置误差曲线对比图；

附图4是现有的基于Singer模型方法和本发明的基于径向加速度的RA-Singer-EKF方法的速度误差曲线对比图；

附图5是现有的基于Singer模型方法和本发明的基于径向加速度的RA-Singer-EKF方法的加速度误差曲线对比图。

具体实施方式

下面结合附图对本发明的基于径向加速度的RA-Singer-EKF方法进行详细描述。(参照说明书附图2)。

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

其中

    为每次仿真得到的估计值，x为真实值。

表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)

其中， $f_d=\frac{{2\mathrm{xv}}_x+{2\mathrm{xv}}_y}{\mathrm{\&lambda;}\sqrt{x^2+y^2}},$ $k=\frac{{2\mathrm{xa}}_x+{2\mathrm{xa}}_y}{\mathrm{\&lambda;}\sqrt{x^2+y^2}},$ 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)& ...& 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；

步骤4：正交匹配追踪(OMP)

本发明采用正交匹配追踪(OMP)方法对信号进行参数提取，和标准匹配追踪(MP)方法相比，该方法在信号参数估计时精度更高，并且稀疏分解收敛速度更快；首先将过完备原子库中的原子与信号进行匹配程度比较，选择与信号最匹配的一组基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^is，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}{\mathrm{\&Sigma;}}}<R^{i}s,g_{\mathrm{ri}}>g_{\mathrm{ri}}+R^{L}s$

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

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

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

(1)经过上述对信号s(t)的OMP分解，计算出原子库中匹配的原子，得到信号的稀疏解能量图

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，最后根据以下公式得到机动目标的径向加速度和速度：

α＝kλ/2，v＝f_dλ/2

步骤6：将得到的径向加速度和径向速度估计值送到雷达数据处理计算机，进行极(球)坐标系到直角坐标系的转换

根据速度与状态向量几何关系图(参照说明书附图1)可以得到径向速度和直角坐标系下X轴方向速度v_x和Y轴方向速度v_y的关系：

$\left\{\begin{array}{c}{v}_{\mathrm{\&rho;}}=v_{\mathrm{\&rho;x}}+v_{\mathrm{\&rho;y}}\\ v_{\mathrm{\&rho;x}}=v_x(x/\sqrt{x^2+y^2})\\ v_{\mathrm{\&rho;y}}=v_y(y/\sqrt{x^2+y^2})\end{array}\right.---\left(7\right)$

同理可得径向加速度和直角坐标系下X轴方向速度α_x和Y轴方向速度α_y的关系：

$\left\{\begin{array}{c}{a}_{\mathrm{\&rho;}}=a_{\mathrm{\&rho;x}}+a_{\mathrm{\&rho;y}}\\ a_{\mathrm{\&rho;x}}=a_x(x/\sqrt{x^2+y^2})\\ a_{\mathrm{\&rho;y}}=a_y(y/\sqrt{x^2+y^2})\end{array}\right.---\left(8\right)$

步骤7：建立Singer模型

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

$R\left(\mathrm{\&tau;}\right)=E\left[a\left(t\right)a(t+\mathrm{\&tau;})\right]={\mathrm{\&sigma;}}_a^2{e}^{-\mathrm{\&alpha;}\left|\mathrm{\&tau;}\right|}$

式中，

α是在区间(t，t+τ)内决定目标机动特性参数的待定参数；α为机动频率，即机动时间常数的倒数，其取值与具体机动状态有关；为目标的加速度方差，可由α(t)概率密度模型求取。

对时间相关函数R(τ)应用Wiener-Kolmogorov白化程序后，机动加速度可用输入为白噪声的一阶时间相关模型来表征，即

$\stackrel{\&CenterDot;}{a}\left(t\right)=-\mathrm{\&alpha;a}\left(t\right)+\sqrt{2\mathrm{\&alpha;}{\mathrm{\&sigma;}}_a^2}u\left(t\right)=-\mathrm{\&alpha;a}\left(t\right)+{\mathrm{\&omega;}}_1\left(t\right)$

式中u(t)表示均值为零、方差为1的高斯白噪声；ω₁(t)表示均值为零、方差为

的高斯白噪声。

步骤8：确定量测方程：

设状态向量X表示为

X＝[x v_x α_x y v_y α_y]^T                                  (9)

式中x、v_x、α_x和y、v_y、α_y和分别表示直角坐标系下X轴和Y轴的位置、速度和加速度。

设量测向量Z表示为：

$Z={\left[\begin{array}{cccc}{z}_x^d& z_y^d& z_{\mathrm{rv}}& z_{\mathrm{ra}}\end{array}\right]}^T---\left(10\right)$

式中

和

分别表示极-直坐标的无偏转换后直角坐标系下X轴和Y轴的位置；z_rv、z_ra表示信号处理阶段得到的径向速度和径向加速度估计值，设其量测误差

均为零均值高斯白噪声且其协方差为

由(7)、(8)、(9)、(10)式可得状态向量到量测向量的方程：

Z(k+1|k)＝H[X(k+1|k)]+V(K+1)                        (11)

式中H[X(k+1|k)]表示量测矩阵，V(K+1)表示方差为R(k+1)的量测噪声。

其中

$H(K+1)={\left[\begin{array}{c}x\\ y\\ v_x(x/\sqrt{x^2+y^2})+v_y(y/\sqrt{x^2+y^2})\\ a_x(x/\sqrt{x^2+y^2})+a_y(y/\sqrt{x^2+y^2})\end{array}\right]}_{X=X(k+1|k)}$

$R(k+1)=\mathrm{blkdig}(R^d(k+1),{\mathrm{\&sigma;}}_{\mathrm{rv}}^2,{\mathrm{\&sigma;}}_{\mathrm{ra}}^2)$

雅克比矩阵为：

$H_X(K+1)=\left[{\&dtri;}_X{h}^{\&prime;}\left(X\right)\right]={\left[{\left[\begin{array}{cccccc}\&PartialD;/\&PartialD;x& \&PartialD;/\&PartialD;v_x& \&PartialD;/\&PartialD;a_x& \&PartialD;/\&PartialD;y& \&PartialD;/\&PartialD;v_y& \&PartialD;/\&PartialD;a_y\end{array}\right]}^T\begin{array}{c}\end{array}\left[\begin{array}{c}x\\ y\\ f_1\left(X\right)\\ f_2\left(X\right)\end{array}\right]\right]}_{X=X(k+1|k)}$

式中

$f_1\left(X\right)=v_x(x/\sqrt{x^2+y^2})+v_y{(y/\sqrt{x^2+y^2})}_{X=X(k+1|k)},$

$f_2\left(X\right)=a_x(x/\sqrt{x^2+y^2})+a_y{(y/\sqrt{x^2+y^2})}_{X=X(k+1|k)};$

步骤9：确定状态方程

假设连续时间状态方程可以表示为

$\frac{d}{\mathrm{dt}}X\left(t\right)=\mathrm{AX}\left(t\right)+B\stackrel{\&OverBar;}{a}\left(t\right)+\mathrm{c\&omega;}\left(t\right)$

式中状态向量X(t)＝[x(t) v_x(t) v_y(t) y(t) α_x(t) α_y(t)]^T分别表示目标X轴、Y轴方向位置、速度加速度； $A=\left[\begin{array}{cccccc}0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& -\mathrm{\&alpha;}& 0& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& -\mathrm{\&alpha;}\end{array}\right],$ $B=\left[\begin{array}{cc}0& 0\\ 0& 0\\ \mathrm{\&alpha;}& 0\\ 0& 0\\ 0& 0\\ 0& \mathrm{\&alpha;}\end{array}\right].$

设采样周期为T′_s，对连续时刻的状态方程式进行离散化处理，可得到下列离散状态方程：

$X(k+1)=\mathrm{\&Phi;}(k+1,k)X\left(k\right|k)+U\left(k\right)\stackrel{\&OverBar;}{a}\left(k\right)+w\left(k\right)---\left(12\right)$

式中X(k|k)＝[x(k|k)v_x(k|k) α_x(k|k) y(k|k) v_y(k|k) α_y(k|k)]^T，

$\mathrm{\&Phi;}(k+1,k)=e^{\mathrm{AT}}=\left[\begin{array}{cccccc}0& T& (-1+\mathrm{\&alpha;T}+e^{-\mathrm{\&alpha;T}})/{\mathrm{\&alpha;}}^2& 0& 0& 0\\ 0& 1& (1-e^{-\mathrm{\&alpha;T}})/\mathrm{\&alpha;}& 0& 0& 0\\ 0& 0& e^{-\mathrm{\&alpha;T}}& 0& 0& 0\\ 0& 0& 0& 0& T& (-1+\mathrm{\&alpha;T}+e^{-\mathrm{\&alpha;T}})/{\mathrm{\&alpha;}}^2\\ 0& 0& 0& 0& 1& (1-e^{-\mathrm{\&alpha;T}})/\mathrm{\&alpha;}\\ 0& 0& 0& 0& 0& e^{-\mathrm{\&alpha;T}}\end{array}\right],$

$U\left(k\right)={\&Integral;}_{\mathrm{kT}}^{(k+1)T}{e}^{A[(k+1)T-\mathrm{\&tau;}]}\mathrm{Bd\&tau;}={\&Integral;}_0^T{e}^{\mathrm{At}}\mathrm{Bdt}=\left[\begin{array}{cc}(-T+\mathrm{\&alpha;}{T}^2/2+(1-e^{-\mathrm{\&alpha;T}})/\mathrm{\&alpha;})/\mathrm{\&alpha;}& 0\\ T-\frac{1-e^{-\mathrm{\&alpha;T}}}{\mathrm{\&alpha;}}& 0\\ 1-e^{-\mathrm{\&alpha;T}}& 0\\ 0& (-T+\mathrm{\&alpha;}{T}^2/2+(1-e^{-\mathrm{\&alpha;T}})/\mathrm{\&alpha;})/\mathrm{\&alpha;}\\ 0& T-\frac{1-e^{-\mathrm{\&alpha;T}}}{\mathrm{\&alpha;}}\\ 0& 1-e^{-\mathrm{\&alpha;T}}\end{array}\right],$

$\stackrel{\&OverBar;}{a}\left(k\right)=\left[\frac{4-\mathrm{\&pi;}}{\mathrm{\&pi;}}{(a}_{\max }-a_x\left(k\right|k),\frac{4-\mathrm{\&pi;}}{\mathrm{\&pi;}}(a_{\max }-a_y\left(k\right|k)\right],$

$w\left(k\right)={\&Integral;}_{\mathrm{kT}}^{(k+1)T}{e}^{A[(k+1)T-\mathrm{\&tau;}]}\mathrm{C\&omega;}\left(\mathrm{\&tau;}\right)\mathrm{d\&tau;}={\&Integral;}_0^T\left[\begin{array}{c}\frac{1-e^{-\mathrm{\&alpha;T}}}{\mathrm{\&alpha;}}\\ e^{-\mathrm{\&alpha;T}}\end{array}\right]\mathrm{\&omega;}(t_{k+1}-t)\mathrm{dt},$

Q(k)＝E[w_j(k)·w_j ^T(k)]；

步骤10：EKF算法：

采用量测方程(11)和状态方程(12)的EKF算法流程为：

$X(k+1|k)=A\left(k\right)X\left(k\right|k)+C\left(k\right)\stackrel{\&OverBar;}{a}\left(k\right)$

式中 $\stackrel{\&OverBar;}{a}\left(k\right)=\mathrm{diag}(a_x\left(k\right|k),a_y\left(k\right|k)).$

P(k+1|k)＝A(k)P(k)A^T(k)+Q(k)

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

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

$\mathrm{Kg}(k+1)=P(k+1|k)H_X^T(k+1){[H_X(k+1)P(k+1|k)H_X^T(k+1)+R(K+1)]}^{-1}$

X(k+1|k+1)＝X(k+1|k)+Kg(k+1)In(k+1)

P(k+1|k+1)＝[I-Kg(k+1)H_X(k+1)]P(k+1|k)[I+Kg(k+1)H_X(k+1)]^T

-Kg(k+1)R(k+1)Kg^T(k+1)

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

Claims (1)

1.一种基于径向加速度的RA-Singer-EKF机动目标跟踪方法，是指在信号处理阶段采用正交匹配追踪(OMP)思想对目标径向加速度和径向速度进行提取，在数据处理阶段将径向加速度和径向速度估计值通过坐标转换引入量测方程和和状态方程中，并采用Singer模型和扩展卡尔曼滤波(EKF)方法实现机动目标的精确跟踪，其特征在于包括以下步骤： 

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

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

步骤2：初始化(设置分解参数) 

T设为雷达的脉冲宽度； 

λ为雷达波长； 

f_s为采样频率； 

T_s＝1/f_s为采样间隔； 

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

k_v设为LFM信号的调频率； 

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

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

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

(1)根据LFM信号回波s(t)的特点，建立原子g_r＝ex_p[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＝[g₁，g₂，…，g_N]^T； 

步骤4：正交匹配追踪(OMP) 

本发明采用正交匹配追踪(OMP)方法对信号进行参数提取，和标准匹配追踪(MP)方法相比，该方法在信号参数估计时精度更高，并且稀疏分解收敛速度更快；首先将过完备原子库中的原子与信号进行匹配程度比较，选择与信号最匹配的一组基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步分解后停止分解，信号被分解为： 

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

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

(1)经过上述对信号s(t)的OMP分解，计算出原子库中匹配的原子，得到信号的稀疏解能量图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，最后根据以下公式得到机动目标的径向加速度和速度： 

α＝kλ/2，v＝f_dλ/2 

步骤6：将得到的径向加速度和径向速度估计值送到雷达数据处理计算机，进行极(球)坐标系到直角坐标系的转换 

根据速度与状态向量几何关系图可以得到径向速度和直角坐标系下X轴方向速度v_x和Y轴方向速度v_y的关系： 

同理可得径向加速度和直角坐标系下X轴方向速度α_x和Y轴方向速度α_y的关系： 

步骤7：建立Singer模型 

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

式中，

α是在区间(t，t+τ)内决定目标机动特性参数的待定参数；α为机动频率，即机动时间常数的倒数，其取值与具体机动状态有关；

为目标的加速度方差，可由α(t)概率密度模型求取。 

对时间相关函数R(τ)应用Wiener-Kolmogorov白化程序后，机动加速度可用输入为白噪声的一阶时间相关模型来表征，即 

式中u(t)表示均值为零、方差为1的高斯白噪声；ω₁(t)表示均值为零、方差为的高斯白噪声。 

步骤8：确定量测方程： 

设状态向量X表示为 

X＝[x v_x α_x y v_y α_y]^T    (3) 

式中x、v_x、α_x和_y、v_y、α_y和分别表示直角坐标系下X轴和Y轴的位置、速度和加速度。 

设量测向量Z表示为： 

式中

和

分别表示极-直坐标的无偏转换后直角坐标系下X轴和Y轴的位置；z_rv、z_ea表示信号处理阶段得到的径向速度和径向加速度估计值，设其量测误差

均为零均值高斯白噪声且其协方差为

由(1)、(2)、(3)、(4)式可得状态向量到量测向量的方程： 

Z(k+1|k)＝H[X(k+1|k)]+V(K+1)                      (5) 

式中H[X(k+1|k)]表示量测矩阵，V(K+1)表示方差为R(k+1)的量测噪声。 

其中 

雅克比矩阵为： 

式中 

步骤9：确定状态方程 

假设连续时间状态方程可以表示为 

式中状态向量X(t)＝[x(t)v_x(t)v_y(t)_y(t)α_x(t)α_y(t)]^T分别表示目标X轴、Y轴方向位置、速度加速度；

设采样周期为T′_s，对连续时刻的状态方程式进行离散化处理，可得到下列离散状态方程： 

式中X(k|k)＝[x(k|k)v_x(k|k)α_x(k|k)y(k|k)v_y(k|k)α_y(k|k)]^T， 

Q(k)＝E[w_j(k)·w_j ^T(k)]； 

步骤10：EKF算法： 

采用量测方程(5)和状态方程(6)的EKF算法流程为： 

式中

P(k+1|k)＝A(k)P(k)A^T(k)+Q(k) 

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

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

X(k+1|k+1)＝X(k+1|k)+Kg(k+1)In(k+1) 

P(k+1|k+1)＝[I-Kg(k+1)H_X(k+1)]P(k+1|k)[I+Kg(k+1)H_X(k+1)]^T

-Kg(k+1)R(k+1)Kg^T(k+1) 。

Images