发明名称 基于径向加速度的RA-Singer-EKF机动目标跟踪算法
摘要 本发明公开了一种基于径向加速度的RA-Singer-EKF机动目标跟踪算法,属于雷达机动目标跟踪领域。本发明的方法可快速准确地提供机动目标的径向加速度和径向速度信息,有效提高雷达对机动目标的跟踪性能。本发明的方法包括以下步骤:(一)对雷达接收信号进行采样,利用匹配追踪(OMP)方法得到目标径向加速度和径向速度;(二)在数据处理阶段将径向加速度和径向速度进行坐标转换并引入量测方程和状态方程;(三)采用Singer模型和扩展卡尔曼滤波(EKF)算法实现机动目标跟踪。本发明能够精确实时地反映目标的机动情况,提高了目标跟踪精度,改善了速度和加速度估计精度,工程实现容易,具有较强的工程应用价值和推广前景。
申请公布号 CN103048658A 申请公布日期 2013.04.17
申请号 CN201210519350.X 申请日期 2012.11.10
申请人 中国人民解放军海军航空工程学院 发明人 王国宏;贾舒宜;张磊;谭顺成;于洪波
分类号 G01S13/66(2006.01)I 主分类号 G01S13/66(2006.01)I
代理机构 代理人
主权项 1.一种基于径向加速度的RA-Singer-EKF机动目标跟踪方法,是指在信号处理阶段采用正交匹配追踪(OMP)思想对目标径向加速度和径向速度进行提取,在数据处理阶段将径向加速度和径向速度估计值通过坐标转换引入量测方程和和状态方程中,并采用Singer模型和扩展卡尔曼滤波(EKF)方法实现机动目标的精确跟踪,其特征在于包括以下步骤: 步骤1:将雷达接收机接收到的线性调频信号s(t)通过采样器以采样间隔T<sub>s</sub>进行采样,变为离散信号s(nT<sub>s</sub>),其中n表示采样点序号;将s(nT<sub>s</sub>)送入雷达信号处理计算机; 在雷达信号处理计算机中执行以下步骤: 步骤2:初始化(设置分解参数) T设为雷达的脉冲宽度; λ为雷达波长; f<sub>s</sub>为采样频率; T<sub>s</sub>=1/f<sub>s</sub>为采样间隔; f<sub>u</sub>设为LFM信号的初始频率; k<sub>v</sub>设为LFM信号的调频率; G((U×V)×N))设为过完备原子库,U×V为原子库中原子的个数,N=T/f<sub>s</sub>; λ<sup>2</sup>设为判断OMP分解是否完成的相干比阈值; 步骤3:形成过完备原子库 (1)根据LFM信号回波s(t)的特点,建立原子g<sub>r</sub>=ex<sub>p</sub>[j2π(f<sub>u</sub>n+k<sub>v</sub>n<sup>2</sup>/2N)],r=1,2,…,N; (2)设定搜索精度和范围,假设搜索范围f<sub>u</sub>的取值为f<sub>u</sub>∈[0,U]△f<sub>u</sub>,u=1,2,…,U,U为起始频率的搜索个数,△f<sub>u</sub>为多普勒单元,搜索范围Δk<sub>v</sub>的取值为k<sub>v</sub>∈[0,V]Δk<sub>v</sub>,Δk<sub>v</sub>为调频率单元,v=1,2,…,V,V为调频率的搜索个数,构造的过完备原子库G为(U×V)×N的矩阵: <img file="FSA00000818564800011.GIF" wi="945" he="292" />G=[g<sub>1</sub>,g<sub>2</sub>,…,g<sub>N</sub>]<sup>T</sup>; 步骤4:正交匹配追踪(OMP) 本发明采用正交匹配追踪(OMP)方法对信号进行参数提取,和标准匹配追踪(MP)方法相比,该方法在信号参数估计时精度更高,并且稀疏分解收敛速度更快;首先将过完备原子库中的原子与信号进行匹配程度比较,选择与信号最匹配的一组基g<sub>r</sub>,将所选的原子利用 Grant-Schmidt正交化方法进行正交化处理,然后将信号在这些正交原子构成的空间下进行分解,使其满足 |&lt;s,g<sub>r</sub>(f<sub>u</sub>,k<sub>v</sub>)&gt;|=sup|&lt;s,g<sub>r</sub>&gt;| 因此,信号可以分解为在最佳原子上的分量和残余两部分,即 s=&lt;s,g<sub>r</sub>&gt;g<sub>r</sub>+R<sup>1</sup>s 其中,R<sup>1</sup>s是用最佳原子对原信号进行最佳匹配后的残余;然后从原子库中将最匹配的这组基删掉,接下来对最佳匹配后的残余可以不断进行上面同样的分解过程,即: R<sup>i</sup>s=&lt;R<sup>i</sup>s,g<sub>ri</sub>&gt;g<sub>ri</sub>+R<sup>i+1</sup>s 其中g<sub>ri</sub>满足: &lt;R<sup>i</sup>s,g<sub>ri</sub>&gt;=sup&lt;R<sup>i</sup>s,g<sub>r</sub>&gt; 经过i次迭代后能量衰减程度△R<sup>i+1</sup>为:ΔR<sup>i+1</sup>=‖R<sup>i</sup>s‖<sup>2</sup>-‖R<sup>i+1</sup>s‖<sup>2</sup>=ω<sup>2</sup>‖R<sup>i</sup>s‖<sup>2</sup>,其中,‖R<sup>i</sup>s‖<sup>2</sup>和‖R<sup>i+1</sup>s‖<sup>2</sup>分别是信号s的i阶和i+1阶残余能量;ω<sup>2</sup>=λ<sup>2</sup>(R<sup>i</sup>s)表示信号衰减的速率;当λ<sup>2</sup>(R<sup>i</sup>s)≥λ<sup>2</sup>时,停止分解;假设经过L步分解后停止分解,信号被分解为: <img file="FSA00000818564800021.GIF" wi="545" he="118" />用少量的原子L(相对于信号长度N而言,L《N)表示信号的主要成分,即: <img file="FSA00000818564800022.GIF" wi="415" he="118" />步骤5:径向加速度和速度的确定 (1)经过上述对信号s(t)的OMP分解,计算出原子库中匹配的原子,得到信号的稀疏解能量图<img file="FSA00000818564800023.GIF" wi="93" he="76" />i∈(0,U×V);(2)在i∈(0,U×V)范围搜索最大的峰值; (3)找到最大峰值的坐标Am(i),并在原子库G中查找此坐标的位置的g<sub>r</sub>(f<sub>u</sub>,k<sub>v</sub>),就可以得到目标的初始频率f<sub>u</sub>和调频率k<sub>v</sub>,i=1,2,…,N,最后根据以下公式得到机动目标的径向加速度和速度: α=kλ/2,v=f<sub>d</sub>λ/2 步骤6:将得到的径向加速度和径向速度估计值送到雷达数据处理计算机,进行极(球)坐标系到直角坐标系的转换 根据速度与状态向量几何关系图可以得到径向速度和直角坐标系下X轴方向速度v<sub>x</sub>和Y轴方向速度v<sub>y</sub>的关系: <img file="FSA00000818564800031.GIF" wi="1239" he="247" />同理可得径向加速度和直角坐标系下X轴方向速度α<sub>x</sub>和Y轴方向速度α<sub>y</sub>的关系: <img file="FSA00000818564800032.GIF" wi="1221" he="236" />步骤7:建立Singer模型 将机动加速度作为零均值时间相关色噪声建模。该模型中设机动加速度的概率密度函数在[-A<sub>max</sub>,A<sub>max</sub>]上近似服从零均值均匀分布,其时间相关函数为指数衰减形式: <img file="FSA00000818564800033.GIF" wi="538" he="55" />式中,<img file="FSA00000818564800034.GIF" wi="68" he="46" />α是在区间(t,t+τ)内决定目标机动特性参数的待定参数;α为机动频率,即机动时间常数的倒数,其取值与具体机动状态有关;<img file="FSA00000818564800035.GIF" wi="40" he="46" />为目标的加速度方差,可由α(t)概率密度模型求取。对时间相关函数R(τ)应用Wiener-Kolmogorov白化程序后,机动加速度可用输入为白噪声的一阶时间相关模型来表征,即 <img file="FSA00000818564800036.GIF" wi="883" he="77" />式中u(t)表示均值为零、方差为1的高斯白噪声;ω<sub>1</sub>(t)表示均值为零、方差为<img file="FSA00000818564800037.GIF" wi="99" he="51" />的高斯白噪声。步骤8:确定量测方程: 设状态向量X表示为 X=[x v<sub>x</sub> α<sub>x</sub> y v<sub>y</sub> α<sub>y</sub>]<sup>T</sup>    (3) 式中x、v<sub>x</sub>、α<sub>x</sub>和<sub>y</sub>、v<sub>y</sub>、α<sub>y</sub>和分别表示直角坐标系下X轴和Y轴的位置、速度和加速度。 设量测向量Z表示为: <img file="FSA00000818564800038.GIF" wi="1283" he="70" />式中<img file="FSA00000818564800039.GIF" wi="35" he="50" />和<img file="FSA000008185648000310.GIF" wi="35" he="54" />分别表示极-直坐标的无偏转换后直角坐标系下X轴和Y轴的位置;z<sub>rv</sub>、z<sub>ea</sub>表示信号处理阶段得到的径向速度和径向加速度估计值,设其量测误差<img file="FSA000008185648000311.GIF" wi="61" he="43" /><img file="FSA000008185648000312.GIF" wi="39" he="42" />均为零均值高斯白噪声且其协方差为<img file="FSA000008185648000313.GIF" wi="73" he="46" /><img file="FSA000008185648000314.GIF" wi="77" he="47" />由(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)的量测噪声。 其中 <img file="FSA00000818564800041.GIF" wi="1021" he="302" /><img file="FSA00000818564800042.GIF" wi="627" he="50" />雅克比矩阵为: <img file="FSA00000818564800043.GIF" wi="1580" he="259" />式中 <img file="FSA00000818564800044.GIF" wi="931" he="72" /><img file="FSA00000818564800045.GIF" wi="943" he="72" />步骤9:确定状态方程 假设连续时间状态方程可以表示为 <img file="FSA00000818564800046.GIF" wi="539" he="89" />式中状态向量X(t)=[x(t)v<sub>x</sub>(t)v<sub>y</sub>(t)<sub>y</sub>(t)α<sub>x</sub>(t)α<sub>y</sub>(t)]<sup>T</sup>分别表示目标X轴、Y轴方向位置、速度加速度;<img file="FSA00000818564800047.GIF" wi="516" he="366" /><img file="FSA00000818564800048.GIF" wi="234" he="365" />设采样周期为T′<sub>s</sub>,对连续时刻的状态方程式进行离散化处理,可得到下列离散状态方程: <img file="FSA00000818564800049.GIF" wi="1439" he="57" />式中X(k|k)=[x(k|k)v<sub>x</sub>(k|k)α<sub>x</sub>(k|k)y(k|k)v<sub>y</sub>(k|k)α<sub>y</sub>(k|k)]<sup>T</sup>, <img file="FSA000008185648000410.GIF" wi="1368" he="375" /><img file="FSA00000818564800051.GIF" wi="1909" he="478" /><img file="FSA00000818564800052.GIF" wi="885" he="105" /><img file="FSA00000818564800053.GIF" wi="1024" he="168" />Q(k)=E[w<sub>j</sub>(k)·w<sub>j</sub><sup>T</sup>(k)]; 步骤10:EKF算法: 采用量测方程(5)和状态方程(6)的EKF算法流程为: <img file="FSA00000818564800054.GIF" wi="742" he="68" />式中<img file="FSA00000818564800055.GIF" wi="618" he="67" />P(k+1|k)=A(k)P(k)A<sup>T</sup>(k)+Q(k) Z(k+1|k)=H(k+1)X(k+1|k) In(k+1)=Z(k+1)-Z(k+1|k) <img file="FSA00000818564800056.GIF" wi="1543" he="81" />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<sub>X</sub>(k+1)]P(k+1|k)[I+Kg(k+1)H<sub>X</sub>(k+1)]<sup>T</sup>-Kg(k+1)R(k+1)Kg<sup>T</sup>(k+1) 。
地址 264001 山东省烟台二马路188号科研部
您可能感兴趣的专利