机动目标跟踪与反跟踪(附录).doc

.
参赛密码
（由组委会填写）
全
第十一届华为杯全国研究生数学建模竞赛
学 校
空军工程大学
参赛队号
90045035
队员姓名
1.唐 茂
2.史 密
3.李世杰
参赛密码
（由组委会填写）

第十一届华为杯全国研究生数学建模竞赛
题 目 机动目标的跟踪与反跟踪
摘 要：
为解决机动目标的跟踪与反跟踪问题，本文综合运用航迹拟合、数据关联、滤波原理，建立了目标的运动模型，并基于自适应卡尔曼滤波原理建立了目标跟踪模型。运用所建立的模型进行仿真，估计和预测了目标机动情况，提出了跟踪与反跟踪策略。
针对问题一，首先将雷达量测数据按照坐标转换公式统一到地心地固坐标系中，采用多项式拟合方法估计目标初始状态；其次综合考虑，确定建立Singer模型来跟踪目标运动状况；最后设计了卡尔曼滤波器估计目标的状态信息，得到目标的航迹。
针对问题二，进行降维处理简化模型，通过数据关联算法分离测量信息，确定测量数据属于哪个目标；依据两个目标的测量数据分别进行卡尔曼滤波得到两个目标各自的航迹图；分析了当雷达一段时间只有一个回波点迹时，怎样使得航迹不丢失。
针对问题三，依据Data3中的数据描点做出点迹图；由于目标运动轨迹较为简单，采用多项式拟合得到目标的拟合航迹一。再依据题目要求，采用问题一所建立的跟踪模型进行处理，得到拟合航迹二。分别对两条航迹进行二次差分得到目标加速度随时间变化规律，并作对比。
针对问题四，依据最小二乘估计，对问题三得到的拟合航迹一做预测；目标着落点即预测航迹与地平面的交点，联立方程求此着落点在以雷达为原点的站心切平面坐标系中坐标。最后对比分析两种预测算法的复杂度，发现复杂度随预测步数指数增加。
针对问题五，建立两种机动动作模型；分析了问题一的跟踪模型对采用以上机动动作的目标的跟踪能力；确定了跟踪策略；采用交互式多模型算法对问题一建立的模型进行完善；调整了跟踪策略。
关键词：目标跟踪；数据关联；卡尔曼滤波；Singer模型；数据拟合；实时预测
一、问题重述
目标跟踪[1]是为了维持对目标当前状态的估计，同时也是对传感器接收到的量测进行处理的过程（多源信息融合），它在军事和民用领域都已经得到了广泛的应用。同时，被跟踪目标为了反跟踪往往会进行机动或释放干扰，这对目标跟踪技术提出了更高的要求。
本题介绍了机动目标跟踪的难点以及目标跟踪处理流程，给定了三组机动目标的测量数据以及雷达坐标和测量误差。其中Data1给定了多个雷达站在不完全相同时刻获得的单个机动目标的测量数据；Data2给定了某个雷达站获得的两个机动目标的测量数据；Data3给定了某个雷达站获得的空间目标的测量数据。
根据已知内容，需要解决的问题如下：
问题一：
1）根据Data1.txt中的数据，分析目标机动发生的时间范围，统计目标加速度的大小和方向。
2）建立对该目标的跟踪模型，利用多个雷达的测量数据估计出目标的航迹。
问题二：
1）根据Data2.txt中的数据，完成各目标的数据关联，形成相应的航迹，并阐明所采用或制定的准则。
2）通过处理保证若出现雷达一段时间只有一个回波点迹的状况，可使得航迹不丢失。
问题三：
1)根据Data3.txt中的数据，分析空间目标的机动变化规律(目标加速度随时间变化)。
2)若采用第1问的跟踪模型进行处理，结果会有哪些变化？
问题四：
对第3问的目标轨迹进行实时预测，估计该目标的着落点的坐标，给出详细结果，并分析算法复杂度。
问题五：
1）Data2.txt数据中的两个目标已被雷达锁定跟踪。分析该目标应采用怎样的有利于逃逸的策略与方案来应对之前所建立的跟踪模型。
2）分析为了保持对目标的跟踪，跟踪策略又应该如何相应地变换 。
二、模型假设
1，假设雷达均处在正常工作状态。
2，假设电磁波在传播过程中符合几何光学定律。
3, 假设目标的机动符合现有的技术水平。
3，假设目标在飞行过程中不会发生坠毁等突发状况。
4，忽略多径效应对测量结果的影响。
三、基本符号说明

目标到雷达距离

滤波增益

俯仰角

均方误差

方位角

误差向量

大地纬度

机动频率

大地经度

统计距离

轴坐标

跟踪门限

轴坐标

多项式系数矩阵

轴坐标

飞行坡度

当前时刻状态

状态转移矩阵

状态估计

残差

量测量

雷达波束宽
四、问题一模型的建立与求解
4.1 问题一分析
Data1给定了多个雷达站在不完全相同时刻测量获得的单个机动目标的距离、方位、俯仰和时间数据信息。由于雷达量测值为目标距离方位和俯仰，为使得观测数据更直观，并简化下文所建立模型的状态方程和量测方程，需要进行坐标转换以统一坐标系：
首先将距离，方位和仰角的极坐标信息转换到站心切平面坐标系中。题目中站心切平面坐标系为：原点设为雷达中心，传感器中心点与当地纬度切线方向指向东为x轴，传感器中心点与当地经度切线方向指向北为y轴，地心与传感器中心连线指向天向的为z轴，目标方位指北向顺时针夹角（从y轴正向向x轴正向的夹角，范围为0~360°），目标俯仰指传感器中心点与目标连线和地平面的夹角，据此，极坐标转换为站心坐标转换公式[2]为：
(4.1)
式中是目标在站心切平面坐标系中的坐标，r为目标量测距离，H为目标俯仰角，A为目标方位角。
依据坐标转换得到的数据，在Matlab中编程（见附件），直接呈现该机动目标的量测数据，如下图所示：

图4.1 雷达一量测值

图4.2 雷达二量测值

图4.3 雷达三量测值
第二步将目标的站心切平面坐标转换到地心地固坐标系中。公式为：
(4.2)
式中为目标在地心地固坐标系下的坐标，为第i个雷达在地心地固坐标系下的坐标，B为大地纬度，L为大地经度。所得结果如下所示：

图4.4 雷达数据合成图
4.2问题一跟踪建模建立与求解
所谓估计就是根据测量得到与状态有关的数据解算出的计算值。对目标的跟踪过程本质上是根据雷达的量测数据估算目标真实位置的过程。其中，最小二乘估计、最小方差估计、极大后验估计、极大似然估计、线性最小方差估计和维纳滤波都是几种最优估计算法[3]。考虑到滤波算法的复杂程度和目标运动特征未知的情况，本组决定采用卡尔曼滤波的方法。下面简述卡尔曼滤波模型的构建步骤：
设系统k时刻的状态,个分量分别表示目标三维坐标下的位置，速度和加速度。
状态方程为:
(4.3)
量测方程为：
(4.4)
其中为量测噪声，为系统噪声。设量测噪声方差阵，滤波步骤如下：
状态一步预测：
(4.5)
状态估计：
(4.6)
滤波增益为：
(4.7)
一步预测均方误差：
(4.8)
估计均方误差：
(4.9)
只要给定初始值和，根据k时刻的量测就可以推算k时刻的状态。
通过以上分析，为了完整卡尔曼滤波方程，需要获得目标初始的状态信息和目标机动的状态方程，既确定目标运动模型。
4.2.1目标初始值的估计
卡尔曼滤波是一种递推算法，启动时必须先给定初值和。如果选取，则滤波过程中估计始终是无偏的[4]，即。如果滤波的起始时刻有量测量，根据卡尔曼状态估计和滤波增益，为：
(4.10)
事实上，如果系统是一致完全随机可控和一致完全随机可观测的，则卡尔曼滤波器一定是一致逐渐稳定的。随着滤波步数的增加，初值和对滤波值的逐渐减小甚至消失，估计逐渐趋向无偏。
然而在实际求解中，很难准确知道目标在初始时刻的状态信息（包括目标位置、速度和加速度）。为使滤波顺利进行，采用拟合[5]量测信息的方法估计目标初始状态。滤波的稳定性以及拟合算法的可行性将在下节做具体阐述。
由于只需要知道目标运动初始状态，这里拟合雷达1的量测值，得到如下图所示的运动轨迹：

图4.5目标拟合轨迹投影
通过差分运算可以得到目标初始速度和加速度的估算值：

图4.6 目标拟合速度投影

图4.7 目标拟合加速度投影
由上图可以得到滤波算法的初始值：

题目中给出了雷达的误差参数：
雷达标号
测距误差(m)
方位角误差( °)
俯仰角误差(°)
1
50
0.4
0.4
2
40
0.3
0.3
3
60
0.5
0.5
下面根据误差理论[6]确定卡尔曼滤波器的估计均方误差
由雷达极坐标测量误差转换为站心切平面误差时
(4.11)
化简得:
(4.12)
同理可得:
(4.13)
(4.14)
令
(4.15)
(4.16)
误差向量为：
(4.17)
(4.18)
其中
(4.19)
可取 作为初始状态，即:
(4.20)
4.2.2目标运动模型的确定
通过对机动目标跟踪算法的学习，本组决定采用自适应跟踪算法。其中Singer模型是自适应跟踪算法的基本模型。Singer模型[7]是一种全局统计模型，考虑了任何类型目标的机动和目标发生各种机动的可能性。因此Singer模型算法适应于各种机动类型的目标跟踪问题。此外,无需进行机动检测,因而在进行目标跟踪时无时间滞后。
Singer模型算法认为机动模型是相关模型[8]，对目标加速度a（t）作为具有指数自相关的零均值随机过程建模:
(4.21)
式中，是在区间内决定目标机动特征的待定参数；是目标的加速度方差，是机动时间的倒数，即机动频率，通常的经验取值范围是：目标是飞机慢速转弯，取60s，对于逃避机动取20s。
Singer对应的离散时间状态方程为：
(4.22)
(4.23)
其离散时间过程噪声V具有协方差，且Q为对称阵。
其中：
(4.24)
4.2.3仿真结果
在确定了滤波初始值和运动模型后，即可估计目标的运动轨迹。滤波模型中的核心过程[9]如下：
for i=1:l-1
Xkk1(:,i)=Fk*Xkk(:,i); %一步预测
Pkk1=Fk*Pkk*Fk'; %一步预测均方误差
Kk=Pkk1*Hk'*inv(Hk*Pkk1*Hk'+R); %滤波增益
Xkk(:,i+1)=Xkk1(:,i)+Kk*(Zkcl(:,i)-Hk*Xkk1(:,i)); %状态估计
Pkk=Pkk1-Pkk1*Hk'*inv(Hk*Pkk1*Hk'+R)*Hk*Pkk1; %估计均方误差
end
滤波结果如下图所示，

图4.8 滤波结果
由此可知目标在被跟踪期间发生的较为明显的机动时间，具体信息由表1给出
表 1 机动发生时间表
机动序号
发生时间(s)
1
36724.4—36744.4
2
36961.4—37199.4
4.3跟踪模型稳定性研究
上文指出，当时，滤波器的估计从开始就是无偏的，且估计均方误差是最小的。但是在本题的具体实践中，很难得到初始状态的验前统计值。如果滤波器随时间的增加，和各自都逐渐不受其初值的影响，则滤波器是稳定的。下面分析本题中所建立的卡尔曼滤波器的稳定性。
滤波稳定条件最早由卡尔曼提出，参考文献[10]提出了离散系统滤波稳定的充分条件。
离散系统状态方程和量测方程为:
(4.25)
(4.26)
式中
(4.27)
如果系统一致完全随机可控和一致完全随机可观测，且 和正定，则卡尔曼滤波器是一直逐渐稳定的。
完全随机可控指：
(4.27)
完全随机可观测指：
(4.28)
结合本题情况，从物理意义可见系统满足随机可控和随机可观测的条件。从的过程来看：
(4.29)
可见，始终有值，使每一步计算都能利用量测中的最新信息，修正旧估计，得到新的实时估计。一般在去较大值得情况下，随估计过程的进行，是逐渐下降的，下降程度去Q和R有关。而却比大，增大的值与Q有关。因此当下降值与增大值相当时，趋于稳定值，滤波呈稳定状态。
实际在问题求算过程中，因为观测值有限，为了防止滤波稳定时间超过量测时间，所以本文采用拟合算法选取目标状态初始值是有必要的和可行的。
五、问题二模型的建立与求解
5.1问题二分析
Data2给定了雷达获得的两个机动目标的测量数据，为了得到两个机动目标各自的航迹，首先要判断信息属于哪个目标，即进行数据关联[11]。
本文先根据Data2中的数据直接得到目标点迹来粗略观测目标的轨迹，通过坐标变换得到目标在站心切平面坐标系中的坐标，描点可得下图：

图5.1 观测数据描点图
由于目标的不规则运动以及误差的影响，得到的点迹图杂乱无章，很难提取有效信息。因此我们进行了降维处理，得到目标运动轨迹在水平面的投影如下图：

图5.2 目标水平面运动轨迹
从图中可以清楚的看到两条轨迹，即目标飞行包线在水平面的投影。再以这两条粗略的轨迹为基础，依据5.2中的数据关联算法分别得到两条不同的投影，从而得到两个目标的航迹。
5.2问题二模型的建立与求解
1）主成分线性回归得到初始航迹投影
从已给数据可知雷达的扫描周期较小（1s左右），且目标据雷达距离相对较远。可以假定在一小段时间内雷达扫描两个目标的先后顺序没有发生变化，由于目标的机动要受到诸多限制，所以这种假定可以认为是合理的。考查前200个点，分别取奇数点和偶数点得到以下两图，把它们作为两个目标的初始航迹投影粗略值：

图5.3 初始航迹估计
通过主成分线性回归[12]得：

图5.4 初始航迹
由此得到粗略的初始航迹，可以为下面的数据关联算法和卡尔曼滤波提供所需的初始值。
2)数据关联算法得到两个目标的航迹投影
计算目标的预测位置和有效回波之间的距离：

(5.1)
其中为卡尔曼滤波的残差，表示的协方差，为定义的统计距离。当满足时，则保留该点击，不满足的则剔除。其中为跟踪门限（采用矩形跟踪门），如此可以得到与航迹配对的所有可能的点迹。若此时只有一个点迹满足要求,就取该点迹为目标的新点迹。反之，当跟踪门内的点迹大于1个时,则将最小的点迹作为新的点迹。
运用该数据关联算法得到航迹投影如图：

图5.5 水平面航迹图
3）卡尔曼滤波得精确航迹投影
根据卡尔曼滤波精确估计目标的真实状态:

图5.6 卡尔曼滤波航迹图
4）对Z轴数据进行分析，得到三维航迹图
分别对两个目标z轴方向的数据进行分析，得到下图：

图5.7 图5.8
观察两图可知两个目标的飞行高度基本没有变化，可以认为它们只是在水平面内做了机动。综合以上分析可以最终得到目标的航迹图如下：

图5.9 目标三维航迹图
5）雷达一段时间只有一个回波点迹状况下的处理方法
卡尔曼滤波[13]可以根据先验估计值和观测值得到较准确的估计值。当雷达在一段时间内提供不了新的观测值时，可以将下一时刻的估计值作为量测值[14]进行滤波，如此循环得到目标航迹状态估计。显然，这种方法只适合于短时间内的状态估计，如果雷达长时间无法提供新的观测数据就很有可能造成滤波发散，使滤失去作用。
六、问题三模型的建立与求解
6.1问题三分析
Data3给定了某雷达对空间目标的测量数据，通过坐标转换将测量数据转换到站心切平面坐标系，描点可得下图：

图6.1 目标点迹图
该目标的运动轨迹比较简单，分别对x-t，y-t，z-t进行多项式拟合（x、y、z为目标点坐标，t为时间），得到目标拟合后的位置数据，并画出航迹图。再通过差分运算即可得到目标加速度随时间变化情况。根据题目要求，再采用第一问的模型对数据进行处理，分析对结果的影响。
6.2问题三模型的建立与求解
1）采用多项式拟合算法建立目标运动轨迹模型，结果如下图：

图6.2 航迹拟合图
2）分析目标在运动过程中，x、y、z坐标值相对时间t的变化情况，对其进行差分得到目标在三个方向的速度变化情况，再进行一次差分得到目标加速度在三维方向的变化情况，结果如下图：

图6.3 x轴方向加速度变化情况

图6.4 y轴方向加速度变化情况

图6.5 z轴方向加速度变化情况
3）采用第一问的跟踪模型对数据进行处理，得到数据的航迹图如下：

图6.6 卡尔曼滤波跟踪轨迹图
上方的轨迹为采用第一问的跟踪模型得到的结果，下方的轨迹为采用拟合算法得到的轨迹图。
4） 在第一问跟踪模型这一前提下，分析目标在运动过程中x、y、z坐标值相对时间t的变化情况，对其进行两次差分得到目标在三个方向的加速度变化情况结果如下图：

图6.7 x轴方向加速度变化情况

图6.8 y轴方向加速度变化情况

图6.9 z轴方向加速度变化情况
七、问题四的求解
4.1多项式拟合预测原理
以x轴方向为例，为得到x轴方向上x坐标与量测时间的关系，多项式拟合即根据已知的k个量测点，i时刻量测数据位，i=1,2，..K，求出一条m次多项式曲线P(t)，使其在i处的取值尽量接近。设所求多项式
(7.1)
令
(7.2)
求解过程就是选取适合的多项式系数，使得最小。
上式的解可以用矩阵的解来表示，其中
，， (7.3)
记，则
(7.4)
同理可得y轴和z轴上的多项式，三个多项式组合可以得到目标三维的参数方程。
4.2着落点求解
为估计目标的着落点，需要联立大地平面方程和目标运动轨迹方程，方程的解即是目标的预测着落点。
问题三中目标运动参数方程为

当z=0时，求得，
着落点在站心切平面坐标系下为
4.3方法复杂度分析
由多项式拟合原理可以看到，已知k个量测点的情况下，要预测k+1时刻的点需要之前所有k个数据做最小二乘。而最小二乘求解的运算量与矩阵阶数的三次方成正比。当预测第k+2时刻时，算法需要把预测的第k+1个量作为量测值，故需要之前k+1个数据参与运算。分析得知该方法的运算复杂度是随预测点数的增加而成三次方的速度增加。
若采用问题一所建立的卡尔曼滤波模型进行预测，当量测值用完时，为了预测下一点的状态，把下一时刻的一步估计值作为量测值处理。这样做虽然引入了一定的误差，但可以使得卡尔曼能继续向下递推，从而得到期望的预测值。
卡尔曼的特点是采用了地推算法，不同时刻的量测值不必储存下来，而是经实时处理提炼成被估计状态的信息，随着滤波步数的增加，提取出的信息浓度逐渐增加[15]。对于对状态的预测，算法的复杂度不会随预测步数的增加而增加。
八、问题五模型的建立与求解
8.1问题五分析
当目标发现自身已经被雷达锁定后，可以采取一定的机动动作来逃避雷达的跟踪。为了确定目标的逃逸策略，我们建立了两种飞机典型机动动作模型[16]，即盘旋机动模型和蛇形机动模型。以问题一中提出的跟踪模型为基础分析了这两种机动动作对雷达跟踪能力的影响，从而确定目标反跟踪的策略。
为了保持对目标的跟踪，通过交互式多模算法对问题一中的模型进行完善，并仿真分析新模型的跟踪效果，调整跟踪策略。
8.2反跟踪模型的建立与求解
1）建立目标盘旋机动模型
当目标进行盘旋机动时，目标的盘旋运动方程应满足：
(8.1)
式中:为垂直于迎面空气流的升力;为与推力相反的阻力;为飞行坡度;为飞机质量;为目标盘旋机动速度;为目标盘旋半径;为重力加速度。根据目标盘旋运动方程，可以得到目标盘旋半径为：
(8.2)
盘旋是经典的机动方式之一，如果将目标假设为飞机，做动作设计时还要考虑实际飞行中飞行员所能承受的最大载荷。飞机作盘旋时，比飞机平飞时的速度扩大倍，比平飞时所需的推力扩大倍。
2）建立目标蛇形机动模型
目标在进行蛇形机动时，其运动方程满足如下方程组：
(8.3)
其中和为目标机动参数。
3）采用问题一中建立的跟踪模型分别对进行盘旋机动和蛇形机动的目标进行跟踪，得到跟踪结果如下图：

图8.1 盘旋机动时的跟踪结果

图8.2 蛇形机动时的跟踪结果
分析可知，当目标进行持续盘旋机动时，采用问题一的模型跟踪轨迹一段时间后就偏离了真实航迹，当目标的真实航迹与预测值的偏差超过跟踪门限时，不能保持对目标的持续有效跟踪。所以目标可以采取持续进行盘旋机动的方式来进行反跟踪。
跟踪门限的确定采用如下方法：
(8.4)
其中为目标距离雷达距离，为雷达波束宽度。
8.3 改进的跟踪模型
当目标作出反跟踪机动时，为了保持对目标的跟踪，本组采用交互式多模算法[17]对已有模型进行改进。
通过上文分析不难看出，当滤波器中只使用单一的运动模型，可以对特定的机动进行跟踪。但是当目标改变机动模式，单一的模型很难适应采用多种模式机动的目标。故考虑使用交互式多模型（IMM）解决问题。IMM采用多个模型并行工作，模型间通过马尔科夫链进行切换，使用不同的模型匹配目标不同的运动状态，从而完成对机动目标的跟踪。
设计交互式多模型的模型集一般要包含一个非机动模型以及一个或多个机动模型,非机动模型通常采用二阶匀速CV模型或三阶匀加速CA模型,机动模型一般采用Singer模型,国内一般采用“当前”统计模型, 涉及转弯机动一般采用协调转弯CT模型 。
结合本题具体情况，综合滤波准确性和模型算法的复杂度，本组决定使用CV，CA，CT三个模型的交互多模型构建模型。该模型中包含了三个滤波器、一个模型概率估计器、一个交互式作用器和一个估计混合器。模型总体框架如下图所示。

8.3 IMM模型系统框图
图中是基于三个模型基础上的状态估计，为模型j的状态估计。为模型的可能性向量，为模型的概率向量。为k-1时刻j滤波器的输出。为交互作用的结果，它作为k时刻滤波器j的输入。为k时刻的量测量。
该IMM模型包含了三个模型，采样时刻的模型为（j=1,2,3）。、分别为CV、VA模型， 分别为目标的位置、速度、加速度分量；为均值为零，方差为的高斯白噪声。
状态方程：
(8.5)
CV模型坐标状态向量：
(8.6)
一步转移阵：
(8.7)
CA模型模型坐标状态向量：
(8.8)
一步转移阵：
(8.9)
为CT模型[18]，可模拟目标机动转弯时的运动过程。状态变量为，， 分别表示直角坐标系下x轴和y轴的分量。
离散的转弯模型系统状态方程表示为：
(8.10)
在IMM模型中，假设目标在任意时刻的运动状态可以用以上三个模型描述；假设模型间的切换在采样时刻发生，并用马尔科夫链来描述模型切换。
具体的流程如下：
1、参数设置
仿真过程中，如果将状态变量设为，x可以很好描述目标三维空间下三个维度的位置，速度和加速度。但是这样在模型滤波中,部分变量就会表示成9x9的矩阵参与运算,势必需要很大的计算量。考虑算法的复杂度，实际情况中，对空间三维分量分别滤波。这样做可能会因为没有考虑到在坐标转换后各方向之间的协方差影响,从而滤波后不是最优的结果。但是由于分开考虑能够降低计算复杂度,而且如果各个状态都能够得到比测量值更优的值,那么总体来看,目标的位置还是比测量值更加接近真值。这样就实现了目标跟踪过程中滤波的目的。
设X轴方向上状态变量： ；Y轴方向上状态变量；垂直方向上的状态变量为。
下面步骤针对X轴方向状态变量进行滤波。
2、状态估计的交互式作用
设从模型i转移到模型j的转移矩阵为。一般由先验知识设定，根据文献[17]提供的参数，本次模拟设
(8.11)
根据L.Campo提出的与状态驻留时间相关的交互式多模型算法[19]：
令为k-1时刻j滤波器的状态估计，为相应的状态协方差阵，为k-1时刻模型j的概率，且。交互后结果为：
(8.12)
式中
(8.13)

3、模型修正
将，作为k时刻模型j的输入，得到相应的滤波输出为，。
4、模型可能性计算
模型j滤波残差为，相应的协方差为。则：
(8.14)
(8.15)
那么模型j的可能性为：
(8.16)
5、模型概率更新
模型j的概率更新为：
(8.17)
其中
(8.18)
6、模型输出
，为k时刻交互式的输出，则有：
(8.19)
(8.20)
九、总结与展望
为期四天的建模竞赛紧张忙碌而又充实。选题后本组成员对目标跟踪领域的相关知识做了认真的学习和研究，结合本题的具体要求，我们主要做了如下工作：
1. 为准确描述机动目标的运动，建立了目标运动的Singer模型和卡尔曼滤波器；基于Matlab编写跟踪目标的程序；利用Data1的量测信息估计目标位置；分析了目标机动情况。
2. 分析Data2中的数据，采用数据关联方法分离Data2中两个目标的量测信息；最终拟合出目标航迹，实现对多目标的跟踪。
3. 针对Data3中的量测数据，采用多项式拟合算法得出目标航迹；再采用问题一建立的跟踪模型估计目标机动；通过对航迹的估计，寻找目标着落点。
4. 采用问题一的模型，对两种典型机动情况下的目标进行了跟踪，发现目标采用持续盘旋机动可以有效躲避雷达的跟踪；反之为了使雷达保持对目标跟踪，可以采用交互式多模型算法对已有模型进行完善。
由于时间仓促，对机动目标的跟踪与反跟踪问题研究的不是很透彻。主要存在的问题有：建立模型适用范围较窄；自适应能力较差；算法应用过于死板；缺少精密的误差控制过程。在将来的学习过程中有待于进一步完善。
十、参考文献
[1] 权太范.目标跟踪新理论与技术.北京：国防工业出版社，2009
[2] 边少锋，李文魁.卫星导航系统概论.北京：电子工业出版社， 2005
[3] 秦永元，张洪钺，汪叔华等.卡尔曼滤波与组合导航原理.西安：西北工业大学出版社，2012
[4] 方青，梅晓春，张育平.用于机动目标跟踪的Kalman滤波器的设计.雷达科学与技术.2006,4(1):50-55.
[5] 张尚剑，刘永智.用滑动窗多项式拟合法实时预测运动目标轨迹.光电工程.2003,30(4):24-27.
[6] 杨惠连，张涛.误差理论与数据处理.天津：天津大学出版社，1992
[7] 何友，修建娟，张晶炜，关欣等.雷达数据处理及应用（第二版）.北京：电子工业出版社， 2012
[8] 谭顺成，王国宏，王娜.MM-Singer模型的机动目标跟踪算法.火力与指挥控制.2012,37(2)
[9] 赵书兰.Matlab建模与仿真.北京：清华大学出版社，2013
[10] Sweeting D.Some design aspects of the WG13 rigid rotor helicopter.AGARD-CP-86-71
[11] 冯洋.多目标跟踪的数据关联算法研究.西安电子科技大学硕士学位论文.2008
[12] 韩崇昭，朱红艳，段战胜等.多源信息融合.北京：清华大学出版社， 2010
[13] Greg Welch，Gary Bishop.An Introduction to the Kalman Filter. UNC-Chapel Hill,TR 95-041,April 5,2004.
[14] 朱自谦，胡士强等.机载雷达多目标跟踪技术.北京：国防工业出版社， 2013
[15] 孙福明.机动目标跟踪状态估计与数据关联技术的研究.中国科学技术大学博士学位论文.2007
[16] 刘思源，张蕾，姚佩阳，万路军.战术飞行机动动作仿真与实现.空军工程大学学报自然科学版.2011,12(5):30-34.
[17] 陆晶莹.高速高机动目标IMM跟踪算法研究.南京理工大学硕士学位论文.2010
[18] 王硕，刘丽.基于机动转弯目标的自适应交互式多模型算法.计算机仿真.2013:30（4）：169-172.
[19] L.Campo,P.Mookeriee,Y.Bar-Shalom. State Estimation for System with a Sojoum-Time-Dependent Markov Swithing Model.IEEE Trans on Auto,1991,36(2):238-243.
附录
clear;clc;
%%%%%%%%%%%%%%雷达地心地固坐标计算%%%%%%%%%%%%%%%%%%
[x01,y01,z01]=change(122.1,40.5,0);[x02,y02,z02]=change(122.4,41.5,0);[x03,y03,z03]=change(122.7,40.9,0);
load('Data1.mat');
load('Data4.mat');
% %%%%%%%%%%%%1号雷达球-直角坐标变换%%%%%%%%%
logic1=(Data1(:,5)==1);
x1=Data1(logic1,1).*cos(pi/2-Data1(logic1,2)/180*pi).*cos(Data1(logic1,3)/180*pi);
y1=Data1(logic1,1).*sin(pi/2-Data1(logic1,2)/180*pi).*cos(Data1(logic1,3)/180*pi);
z1=Data1(logic1,1).*sin(Data1(logic1,3)/180*pi);
% %%%%%%%%%%%%多项式拟合%%%%%%%%%%%%%%%
tt1=Data1(logic1,4)-Data1(1,4);
p=polyfit(tt1,y1,4);
YY1=p(1)*tt1.^4+p(2)*tt1.^3+p(3)*tt1.^2+p(4)*tt1+p(5);
p1=polyfit(tt1,x1,4);
XX1=p1(1)*tt1.^4+p1(2)*tt1.^3+p1(3)*tt1.^2+p1(4)*tt1+p1(5);
p2=polyfit(tt1,z1,4);
ZZ1=p2(1)*tt1.^4+p2(2)*tt1.^3+p2(3)*tt1.^2+p2(4)*tt1+p2(5);
%%%%%%%%%%%%2号雷达球-直角坐标变换%%%%%%%%%
logic2=(Data1(:,5)==2);
x2=Data1(logic2,1).*cos(Data1(logic2,2)/180*pi).*cos(Data1(logic2,3)/180*pi);
y2=Data1(logic2,1).*sin(Data1(logic2,2)/180*pi).*cos(Data1(logic2,3)/180*pi);
z2=Data1(logic2,1).*sin(Data1(logic2,3)/180*pi);
% %%%%%%%%%%%%多项式拟合%%%%%%%%%%%%%%%
tt2=Data1(logic2,4)-Data1(1,4);
p=polyfit(tt2,y2,4);
YY2=p(1)*tt2.^4+p(2)*tt2.^3+p(3)*tt2.^2+p(4)*tt2+p(5);
p1=polyfit(tt2,x2,4);
XX2=p1(1)*tt2.^4+p1(2)*tt2.^3+p1(3)*tt2.^2+p1(4)*tt2+p1(5);
p2=polyfit(tt2,z2,4);
ZZ2=p2(1)*tt2.^4+p2(2)*tt2.^3+p2(3)*tt2.^2+p2(4)*tt2+p2(5);
%%%%%%%%%%%%3号雷达球-直角坐标变换%%%%%%%%%
logic3=(Data1(:,5)==3);
x3=Data1(logic3,1).*cos(Data1(logic3,2)/180*pi).*cos(Data1(logic3,3)/180*pi);
y3=Data1(logic3,1).*sin(Data1(logic3,2)/180*pi).*cos(Data1(logic3,3)/180*pi);
z3=Data1(logic3,1).*sin(Data1(logic3,3)/180*pi);
%%%%%%%%%%%%%多项式拟合%%%%%%%%%%%%%%%
tt3=Data1(logic3,4)-Data1(1,4);
p=polyfit(tt3,y3,4);
YY3=p(1)*tt3.^4+p(2)*tt3.^3+p(3)*tt3.^2+p(4)*tt3+p(5);
p1=polyfit(tt3,x3,4);
XX3=p1(1)*tt3.^4+p1(2)*tt3.^3+p1(3)*tt3.^2+p1(4)*tt3+p1(5);
p2=polyfit(tt3,z3,4);
ZZ3=p2(1)*tt3.^4+p2(2)*tt3.^3+p2(3)*tt3.^2+p2(4)*tt3+p2(5);
tt=Data1(logic3,4)-Data1(1,4);
TT=linspace(min(tt),max(tt),max(tt)-min(tt)+1);
%%%%%%%%%%%%%%雷达站心坐标-地心地固坐标转换%%%%%%%%%%%%%%%%%%
[xx1,yy1,zz