发明名称 一种基于坐标旋转变换的雷达跟踪方法
摘要 本发明公开了一种基于坐标旋转变换的雷达跟踪方法,该方法基于随机变量相关系数和坐标旋转变换的原理,以达到在定量度量雷达滤波系统模型中的量测方程的非线性程度的基础上,降低雷达滤波系统模型中的量测方程非线性程度,从而提高雷达跟踪效果的目的;具体步骤为:定义二维雷达量测极-直坐标转换线性度ρ;当r<sub>0</sub>、σ<sub>r</sub>和σ<sub>a</sub>为确定的任意数值时,通过坐标旋转变换<img file="dda0000098655750000011.GIF" wi="65" he="43" />使得ρ在<img file="dda0000098655750000012.GIF" wi="214" he="107" />时取得最大值;设置目标运动状态,进行滤波初始化;分析由k时刻到k+1时刻的目标运动状态并建立坐标旋转变换滤波模型下的状态方程和量测方程;在坐标旋转变换滤波模型下,选取滤波算法进行滤波,得到k+1时刻目标状态的估计值;直到达到所设置的雷达跟踪时间长度,跟踪结束。
申请公布号 CN102508238A 申请公布日期 2012.06.20
申请号 CN201110311058.4 申请日期 2011.10.14
申请人 北京理工大学 发明人 陈新亮;曾涛;李春霞
分类号 G01S13/72(2006.01)I;G01S7/295(2006.01)I 主分类号 G01S13/72(2006.01)I
代理机构 北京理工大学专利中心 11120 代理人 郭德忠;李爱英
主权项 1.一种基于坐标旋转变换的雷达跟踪方法,其特征在于,包括:步骤S00、二维雷达量测极-直坐标转换方程为:<maths num="0001"><![CDATA[<math><mrow><mfenced open='{' close=''><mtable><mtr><mtd><mi>x</mi><mo>=</mo><mi>r</mi><mi>cos</mi><mi>a</mi></mtd></mtr><mtr><mtd><mi>y</mi><mo>=</mo><mi>r</mi><mi>sin</mi><mi>a</mi></mtd></mtr></mtable></mfenced><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>1</mn><mo>)</mo></mrow></mrow></math>]]></maths>其中,x为目标在雷达直角坐标系下的量测值的横坐标,y为目标在雷达直角坐标系下的量测值的纵坐标;对式(2)中的x、y分别在(r<sub>0</sub>,a<sub>0</sub>)处进行二元泰勒展开,并保留至一阶项,将x的一阶泰勒展开式记为随机变量g,将y的一阶泰勒展开式记为随机变量k;x≈r<sub>0</sub>cosa<sub>0</sub>+(r-r<sub>0</sub>)cosa<sub>0</sub>-(a-a<sub>0</sub>)r<sub>0</sub>sina<sub>0</sub>□g    (2)y≈r<sub>0</sub>sina<sub>0</sub>+(r-r<sub>0</sub>)sina<sub>0</sub>+(a-a<sub>0</sub>)r<sub>0</sub>cosa<sub>0</sub>□k    (3)其中,r<sub>0</sub>为雷达极坐标下的目标真实距离,a<sub>0</sub>为雷达极坐标下的目标真实方位角,定义二维雷达量测极-直坐标转换线性度ρ为ρ<sub>xg</sub>和ρ<sub>yk</sub>二者的最小值:ρ=min(ρ<sub>xg</sub>,ρ<sub>yk</sub>)    (4)其中:ρ<sub>xg</sub>为随机变量x与随机变量g的相关系数,ρ<sub>yk</sub>为随机变量y与随机变量k的相关系数,且<maths num="0002"><![CDATA[<math><mrow><msub><mi>&rho;</mi><mi>xg</mi></msub><mo>=</mo><msqrt><mfrac><mrow><msubsup><mi>&sigma;</mi><mi>r</mi><mn>2</mn></msubsup><msup><mrow><mo>(</mo><mi>cos</mi><msub><mi>a</mi><mn>0</mn></msub><mo>)</mo></mrow><mn>2</mn></msup><mo>+</mo><msup><msub><mi>r</mi><mn>0</mn></msub><mn>2</mn></msup><msubsup><mi>&sigma;</mi><mi>a</mi><mn>2</mn></msubsup><msup><mrow><mo>(</mo><mi>sin</mi><msub><mi>a</mi><mn>0</mn></msub><mo>)</mo></mrow><mn>2</mn></msup></mrow><mrow><mrow><mo>(</mo><msup><msub><mi>r</mi><mn>0</mn></msub><mn>2</mn></msup><mo>+</mo><msubsup><mi>&sigma;</mi><mi>r</mi><mn>2</mn></msubsup><mo>)</mo></mrow><mn>0.5</mn><mo>[</mo><msup><mi>e</mi><msubsup><mi>&sigma;</mi><mi>a</mi><mn>2</mn></msubsup></msup><mo>+</mo><msup><mi>e</mi><mrow><mo>-</mo><msubsup><mi>&sigma;</mi><mi>a</mi><mn>2</mn></msubsup></mrow></msup><mi>cos</mi><mn>2</mn><msub><mi>a</mi><mn>0</mn></msub><mo>]</mo><mo>-</mo><msup><msub><mi>r</mi><mn>0</mn></msub><mn>2</mn></msup><msup><mrow><mo>(</mo><mi>cos</mi><msub><mi>a</mi><mn>0</mn></msub><mo>)</mo></mrow><mn>2</mn></msup></mrow></mfrac></msqrt><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>5</mn><mo>)</mo></mrow></mrow></math>]]></maths><maths num="0003"><![CDATA[<math><mrow><msub><mi>&rho;</mi><mi>yk</mi></msub><mo>=</mo><msqrt><mfrac><mrow><msubsup><mi>&sigma;</mi><mi>r</mi><mn>2</mn></msubsup><msup><mrow><mo>(</mo><mi>sin</mi><msub><mi>a</mi><mn>0</mn></msub><mo>)</mo></mrow><mn>2</mn></msup><mo>+</mo><msup><msub><mi>r</mi><mn>0</mn></msub><mn>2</mn></msup><msubsup><mi>&sigma;</mi><mi>a</mi><mn>2</mn></msubsup><msup><mrow><mo>(</mo><mi>cos</mi><msub><mi>a</mi><mn>0</mn></msub><mo>)</mo></mrow><mn>2</mn></msup></mrow><mrow><mrow><mo>(</mo><msup><msub><mi>r</mi><mn>0</mn></msub><mn>2</mn></msup><mo>+</mo><msubsup><mi>r</mi><mi>r</mi><mi>a</mi></msubsup><mo>)</mo></mrow><mn>0.5</mn><mo>[</mo><msup><mi>e</mi><msubsup><mi>&sigma;</mi><mi>a</mi><mn>2</mn></msubsup></msup><mo>-</mo><msup><mi>e</mi><mrow><mo>-</mo><msubsup><mi>&sigma;</mi><mi>a</mi><mn>2</mn></msubsup></mrow></msup><mi>cos</mi><mn>2</mn><msub><mi>a</mi><mn>0</mn></msub><mo>]</mo><mo>-</mo><msup><msub><mi>r</mi><mn>0</mn></msub><mn>2</mn></msup><msup><mrow><mo>(</mo><mi>sin</mi><msub><mi>a</mi><mn>0</mn></msub><mo>)</mo></mrow><mn>2</mn></msup></mrow></mfrac></msqrt><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>6</mn><mo>)</mo></mrow></mrow></math>]]></maths>σ<sub>r</sub>为雷达测距噪声标准差,σ<sub>a</sub>为雷达测角噪声标准差;步骤S01、由式(4),当r<sub>0</sub>、σ<sub>r</sub>和σ<sub>a</sub>为确定的任意数值时,二维雷达量测极-直坐标转换线性度ρ在<img file="FDA0000098655720000021.GIF" wi="214" he="107" />时取得最大值,此时,将雷达直角坐标系XOY顺时针旋转角度<img file="FDA0000098655720000022.GIF" wi="65" he="30" />变为旋转雷达直角坐标系X<sup>c</sup>OY<sup>c</sup>,相应地,雷达直角坐标量测变为:<img file="FDA0000098655720000023.GIF" wi="1122" he="143" />其中,x<sup>c</sup>为目标在旋转雷达直角坐标系X<sup>c</sup>OY<sup>c</sup>下的量测值的横坐标;y<sup>c</sup>为目标在旋转雷达直角坐标系X<sup>c</sup>OY<sup>c</sup>下的量测值的纵坐标;坐标旋转矩阵<img file="FDA0000098655720000024.GIF" wi="115" he="63" />为:<img file="FDA0000098655720000025.GIF" wi="1261" he="137" />由式(7)和(8),则在旋转雷达直角坐标系X<sup>c</sup>OY<sup>c</sup>下的雷达直角坐标量测值为:<img file="FDA0000098655720000026.GIF" wi="1178" he="146" />由式(9),并结合所述当r<sub>0</sub>、σ<sub>r</sub>和σ<sub>a</sub>为确定的任意数值时,二维雷达量测极-直坐标转换线性度ρ在<img file="FDA0000098655720000027.GIF" wi="214" he="107" />时取得最大值,通过坐标旋转变换,使得:<img file="FDA0000098655720000028.GIF" wi="1273" he="107" />步骤S02、设置目标运动状态,包括目标运动的初始位置、目标运动的速度、雷达跟踪时间长度和采样时间间隔,进行滤波初始化;步骤S03、目标运动状态由k时刻递推得到k+1时刻的具体步骤如下:①将k+1时刻的雷达极坐标系下的目标方位角的量测值a(k+1)代入式(10)求得k+1时刻的坐标旋转变换的角度值<img file="FDA0000098655720000029.GIF" wi="192" he="64" />②将雷达直角坐标系XOY顺时针旋转角度值<img file="FDA00000986557200000210.GIF" wi="196" he="64" />建立k+1时刻的旋转雷达直角坐标系X<sup>c</sup>OY<sup>c</sup>和相应的旋转雷达极坐标系;③在k+1时刻,分别分析雷达直角坐标系XOY和旋转雷达直角坐标系X<sup>c</sup>OY<sup>c</sup>、雷达极坐标系和旋转雷达极坐标系之间的关系,相应地得到旋转雷达直角坐标系X<sup>c</sup>OY<sup>c</sup>下的状态方程和旋转雷达极坐标系下的量测方程,并将二者记为坐标旋转变换滤波系统模型;a、坐标旋转变换滤波系统模型的状态方程结合式(7),则:在k+1时刻目标的位置在雷达直角坐标系XOY下与其在旋转雷达直角坐标系X<sup>c</sup>OY<sup>c</sup>下的关系为:<img file="FDA0000098655720000031.GIF" wi="1361" he="148" />其中,x<sup>c</sup>(k+1)为在旋转雷达直角坐标系X<sup>c</sup>OY<sup>c</sup>下目标位置分解到x轴的值;x<sup>c</sup>(k+1)为在旋转雷达直角坐标系X<sup>c</sup>OY<sup>c</sup>下目标位置分解到y轴的值;x(k+1)为在雷达直角坐标系XOY下目标位置分解到x轴的值;y(k+1)为在雷达直角坐标系XOY下目标位置分解到y轴的值;<img file="FDA0000098655720000032.GIF" wi="115" he="63" />为对应于坐标旋转变换的角度值<img file="FDA0000098655720000033.GIF" wi="170" he="64" />的坐标旋转矩阵,即:<img file="FDA0000098655720000034.GIF" wi="1480" he="163" />在k+1时刻目标的速度在雷达直角坐标系XOY下与其在旋转雷达直角坐标系X<sup>c</sup>OY<sup>c</sup>下的关系为:<img file="FDA0000098655720000035.GIF" wi="1370" he="150" />其中,<img file="FDA0000098655720000036.GIF" wi="175" he="64" />为在旋转雷达直角坐标系X<sup>c</sup>OY<sup>c</sup>下目标速度分解到x轴的值;<img file="FDA0000098655720000037.GIF" wi="175" he="66" />为在旋转雷达直角坐标系X<sup>c</sup>OY<sup>c</sup>下目标速度分解到y轴的值;v<sub>x</sub>(k+1)为在雷达直角坐标系XOY下目标速度分解到x轴的值;v<sub>y</sub>(k+1)为在雷达直角坐标系XOY下目标速度分解到y轴的值;在k+1时刻目标的加速度在XOY下与其在X<sup>c</sup>OY<sup>c</sup>下的关系为:<img file="FDA0000098655720000038.GIF" wi="1354" he="150" />其中,<img file="FDA0000098655720000041.GIF" wi="178" he="64" />为在旋转雷达直角坐标系X<sup>c</sup>OY<sup>c</sup>下目标加速度分解到x轴的值;<img file="FDA0000098655720000042.GIF" wi="181" he="66" />为在旋转雷达直角坐标系X<sup>c</sup>OY<sup>c</sup>下目标加速度分解到y轴的值;a<sub>x</sub>(k+1)为在雷达直角坐标系XOY下目标加速度分解到x轴的值;a<sub>y</sub>(k+1)为在雷达直角坐标系XOY下目标加速度分解到y轴的值;若在雷达直角坐标系XOY下,仅考虑目标位置和目标速度,即d=2,则目标运动状态向量X(k+1)为:X(k+1)=[x(k+1)v<sub>x</sub>(k+1)y(k+1)v<sub>y</sub>(k+1)]<sup>T</sup>    (15)由式(11)和(12),则旋转雷达直角坐标系X<sup>c</sup>OY<sup>c</sup>下的目标运动状态向量X<sup>c</sup>(k+1)与X(k+1)的关系为:X<sup>c</sup>(k+1)=AX(k+1)    (16)其中,<img file="FDA0000098655720000043.GIF" wi="1218" he="79" /><img file="FDA0000098655720000044.GIF" wi="39" he="39" />为Kronecker积,I<sub>d</sub>为d×d维单位矩阵,由d=2,则取I<sub>2</sub>,进一步地,由式(17),则有:<img file="FDA0000098655720000045.GIF" wi="1654" he="163" />若在雷达直角坐标系XOY下,考虑目标位置、目标速度和目标加速度,则d=3,相应地,在式(17)中,I<sub>d</sub>取I<sub>3</sub>;经过坐标旋转变换,即由雷达直角坐标系XOY变为旋转雷达直角坐标系X<sup>c</sup>OY<sup>c</sup>,相应地,坐标旋转变换滤波系统模型的状态方程为:X<sup>c</sup>(k+1)=f<sup>c</sup>(X<sup>c</sup>(k))+V<sup>c</sup>(k)    (19)其中,X<sup>c</sup>(k+1)为在k+1时刻旋转雷达直角坐标系X<sup>c</sup>OY<sup>c</sup>下的目标运动的状态矢量,f<sup>c</sup>(·)为在k时刻旋转雷达直角坐标系X<sup>c</sup>OY<sup>c</sup>下的目标运动的状态转移函数,V<sup>c</sup>(k)为在k时刻旋转雷达直角坐标系X<sup>c</sup>OY<sup>c</sup>下的目标运动的过程噪声;结合式(16),在k时刻,旋转雷达直角坐标系X<sup>c</sup>OY<sup>c</sup>下的目标运动的状态转移函数f<sup>c</sup>(X<sup>c</sup>(k))与其在雷达直角坐标系XOY下的关系为:f<sup>c</sup>(X<sup>c</sup>(k))=Af(A<sup>-1</sup>X<sup>c</sup>(k))    (20)同理,在k时刻,旋转雷达直角坐标系X<sup>c</sup>OY<sup>c</sup>下的目标运动的过程噪声V<sup>c</sup>(k)与其在雷达直角坐标系XOY下的关系为:V<sup>c</sup>(k)=AV(k)    (21)假定V<sup>c</sup>(k)是零均值的高斯白噪声,结合式(21),则V<sup>c</sup>(k)的方差为:E[V<sup>c</sup>(k)(V<sup>c</sup>(j))<sup>T</sup>]=AQ(k)A<sup>T</sup>δ<sub>kj</sub>    (22)其中,V<sup>c</sup>(j)为在j时刻旋转雷达直角坐标系X<sup>c</sup>OY<sup>c</sup>下的目标运动的过程噪声;b、坐标旋转变换滤波系统模型的量测方程经过步骤S03中步骤②的坐标旋转变换,在k+1时刻,坐标旋转变换滤波系统模型的量测方程为:<img file="FDA0000098655720000051.GIF" wi="1290" he="80" />其中,k+1时刻的旋转量测值Z<sup>c</sup>(k+1)=[r<sup>c</sup>(k+1)a<sup>c</sup>(k+1)]<sup>T</sup>,r<sup>c</sup>(k+1)为k+1时刻目标在旋转雷达极坐标下的量测距离,a<sup>c</sup>(k+1)为k+1时刻目标在旋转雷达极坐标系下的量测方位角;W<sup>c</sup>(k+1)为k+1时刻旋转雷达直角坐标系X<sup>c</sup>OY<sup>c</sup>下的量测噪声;由于雷达直角坐标转换到雷达极坐标的坐标转换关系与旋转雷达直角坐标转换到旋转雷达极坐标的坐标转换关系是相同的,则有:h<sup>c</sup>(·)=h(·)    (24)即:<maths num="0004"><![CDATA[<math><mrow><msup><mi>h</mi><mi>c</mi></msup><mo>[</mo><msup><mi>X</mi><mi>c</mi></msup><mrow><mo>(</mo><mi>k</mi><mo>+</mo><mn>1</mn><mo>)</mo></mrow><mo>]</mo><mo>=</mo><mfenced open='[' close=']'><mtable><mtr><mtd><msqrt><msup><mrow><mo>[</mo><msup><mi>x</mi><mi>c</mi></msup><mrow><mo>(</mo><mi>k</mi><mo>+</mo><mn>1</mn><mo>)</mo></mrow><mo>]</mo></mrow><mn>2</mn></msup><mo>+</mo><msup><mrow><mo>[</mo><msup><mi>y</mi><mi>c</mi></msup><mrow><mo>(</mo><mi>k</mi><mo>+</mo><mn>1</mn><mo>)</mo></mrow><mo>]</mo></mrow><mn>2</mn></msup></msqrt></mtd></mtr><mtr><mtd><mi>arctan</mi><mo>[</mo><msup><mi>y</mi><mi>c</mi></msup><mrow><mo>(</mo><mi>k</mi><mo>+</mo><mn>1</mn><mo>)</mo></mrow><mo>/</mo><msup><mi>x</mi><mi>c</mi></msup><mrow><mo>(</mo><mi>k</mi><mo>+</mo><mn>1</mn><mo>)</mo></mrow><mo>]</mo></mtd></mtr></mtable></mfenced><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>25</mn><mo>)</mo></mrow></mrow></math>]]></maths>由于旋转雷达极坐标系与雷达极坐标系的坐标原点是相同,故坐标旋转不改变目标量测距离,即r<sup>c</sup>(k+1)=r(k+1);而目标在旋转雷达极坐标系中的量测方位角a<sup>c</sup>(k+1)相比于雷达极坐标系下的量测方位角发生了变化,变化的角度值为坐标旋转角度值<img file="FDA0000098655720000061.GIF" wi="195" he="64" />并且测距噪声和测角噪声并不随着坐标旋转的变化而变化,则有:Z<sup>c</sup>(k+1)=h<sup>c</sup>[X<sup>c</sup>(k+1)]+W<sup>c</sup>(k)    (26)Z<sup>c</sup>(k+1)=Z(k+1)+B(k+1)        (27)<img file="FDA0000098655720000062.GIF" wi="1261" he="142" />W<sup>c</sup>(k)=W(k)                   (29)其中,W<sup>c</sup>(k)为k时刻旋转雷达直角坐标系X<sup>c</sup>OY<sup>c</sup>下的量测噪声;步骤S04、在坐标旋转变换滤波模型下,选取滤波算法进行滤波,得到k+1时刻目标状态的估计值;步骤S05、重复步骤S03和步骤S04,直到达到步骤S02中设置的雷达跟踪时间长度,跟踪结束。
地址 100081 北京市海淀区中关村南大街5号