基于ADSB系统航迹点的预测方法和系统.pdf

一种基于ADS-B系统航迹点的预测方法和系统，其方法包括步骤：将接收到的大地坐标系航迹点转化为直角坐标系航迹点并按接收时间进行排列；判断接收到的航迹点是否为初始航迹点；若是，根据航迹点的状态信息正向预测下一时刻航迹点的状态预测值、下一时刻航迹点的测量预测值以及计算下一时刻航迹点的状态预测协方差矩阵、下一时刻航迹点的残差向量协方差矩阵；计算下一时刻航迹点的衰减因子并由其计算下一时刻滤波器增益；对下一时刻航迹点的状态预测值以及下一时刻航迹点的状态预测协方差矩阵进行修正，得到最终正向预测的下一时刻航迹点的状态预测值以及最终下一时刻航迹点的状态预测协方差矩阵。根据本发明方案，能够精确及时地对航迹点进行预测。

1.  一种基于ADS-B系统航迹点的预测方法，其特征在于，包括步骤：
将接收到的大地坐标系航迹点转化为直角坐标系航迹点并按接收时间对收到的航迹点进行排列；
判断接收到的航迹点是否为初始航迹点；
若是，根据航迹点的状态信息正向预测下一时刻航迹点的状态预测值以及计算下一时刻航迹点的状态预测协方差矩阵，并正向预测下一时刻航迹点的测量预测值以及计算下一时刻航迹点的残差向量协方差矩阵；
根据下一时刻航迹点的测量预测值以及下一时刻航迹点的残差向量协方差矩阵计算下一时刻航迹点的衰减因子以及根据该衰减因子计算下一时刻滤波器增益；
以下一时刻航迹点的测量预测值、下一时刻航迹点的残差向量协方差矩阵以及下一时刻滤波器增益对下一时刻航迹点的状态预测值以及下一时刻航迹点的状态预测协方差矩阵进行修正，得到最终正向预测的下一时刻航迹点的状态预测值以及最终下一时刻航迹点的状态预测协方差矩阵。

2.  根据权利要求1所述的基于ADS-B系统航迹点的预测方法，其特征在于，所述根据航迹点的状态信息正向预测下一时刻航迹点的状态预测值以及计算下一时刻航迹点的状态预测协方差矩阵，并正向预测下一时刻航迹点的测量预测值以及计算下一时刻航迹点的残差向量协方差矩阵通过下述方程式计算得到：
X^(k+1/k)=Φ(k)X^(k/k)+U(k)a‾(k)]]>
P(k+1/k)＝Φ(k)P(k/k)Φ^T(k)+Q(k)
Y^(k+1/k)=H(k)X^(k+1/k)]]>
S(k+1)＝H(k)P(k+1/k)H^T(k)+R(k)
其中，为下一时刻航迹点的状态预测值，P(k+1/k)为下一时刻航迹点的状态预测协方差矩阵，为下一时刻航迹点的测量预测值，S(k+1)为下一时刻航迹点的残差向量协方差矩阵；
所述航迹点状态信息包括当前时刻航迹点的状态预测值当前时刻航迹点的状态预测协方差矩阵P(k/k)、当前时刻正向状态转移矩阵Φ(k)、当前时刻控制矩阵U(k)、当前时刻机动加速度均值当前时刻正向噪声方差矩阵 Q(k)、当前时刻正向状态转移矩阵的转置矩阵Φ^T(k)、当前时刻测量矩阵H(k)、当前时刻测量矩阵的转置矩阵H^T(k)、当前时刻观测噪声R(k)。

3.  根据权利要求2所述的基于ADS-B系统航迹点的预测方法，其特征在于，所述根据下一时刻航迹点的测量预测值以及下一时刻航迹点的残差向量协方差矩阵计算下一时刻航迹点的衰减因子以及根据该衰减因子计算下一时刻滤波器增益通过下述方程式计算得到：
Y~(k+1)=|Y(k+1)-Y^(k+1/k)|]]>
r(k+1)=1|Y~(k+1)|≤3S(k+1){3S(k+1)|Y~(k+1)|}1/2|Y~(k+1)|>3S(k+1)]]>
K(k+1)＝r(k+1)P(k+1)H^T(k)S^-1(k+1)
其中，为新息，Y(k+1)为下一时刻航迹点的测量值，r(k+1)为下一时刻航迹点的衰减因子，K(k+1)为下一时刻滤波器增益。

4.  根据权利要求3所述的基于ADS-B系统航迹点的预测方法，其特征在于，所述以下一时刻航迹点的测量预测值、下一时刻航迹点的残差向量协方差矩阵以及下一时刻滤波器增益对下一时刻航迹点的状态预测值以及下一时刻航迹点的状态预测协方差矩阵进行修正，得到最终正向预测的下一时刻航迹点的状态预测值以及最终下一时刻航迹点的状态预测协方差矩阵通过下述方程式计算得到：
X^(k+1/k+1)=X^(k+1/k)+K(k+1)[Y(k+1)-H(k)X^(k+1/k)]]]>
P(k+1/k+1)＝P(k+1/k)-K(k+1)S(k+1)K^T(k+1)
其中，为修正后最终预测的下一时刻航迹点的状态预测值，P(k+1/k+1)为修正后最终的下一时刻航迹点的状态预测协方差矩阵，K^T(k+1)为下一时刻滤波器增益转置矩阵。

5.  根据权利要求1所述的基于ADS-B系统航迹点的预测方法，其特征在于，在判断接收到的航迹点不是初始航迹点时，包括步骤：
根据航迹点的状态信息反向预测前一时刻航迹点的状态预测值以及计算前 一时刻航迹点的状态预测协方差矩阵，并计算前一时刻增益矩阵；
根据前一时刻增益矩阵对前一时刻航迹点的状态预测值以及计算前一时刻航迹点的状态预测协方差矩阵进行修正，得到修正后最终反向预测的前一时刻航迹点的状态预测值以及前一时刻航迹点的状态预测协方差矩阵。

6.  根据权利要求5所述的基于ADS-B系统航迹点的预测方法，其特征在于，所述根据航迹点的状态信息反向预测前一时刻航迹点的状态预测值以及计算前一时刻航迹点的状态预测协方差矩阵，并计算前一时刻增益矩阵通过下述方程式计算得到：
X^(k-1/k)=FX^(k/k)]]>
P(k-1/k)＝FP(k/k)F^T+Q
K＝P(k-1/k)H^T/(HP(k-1/k)H^T+R)
其中，为前一时刻航迹点的状态预测值，P(k-1/k)为前一时刻航迹点的状态预测协方差矩阵，为当前时刻航迹点的状态预测值，P(k/k)为当前时刻航迹点的状态预测协方差矩阵，K为前一时刻增益矩阵，H为测量矩阵，H^T为测量矩阵的转置矩阵，R为反向观测噪声均值；所述航迹点的状态信息包括反向预测状态转移矩阵F、反向随机噪声方差Q。

7.  根据权利要求6所述的基于ADS-B系统航迹点的预测方法，其特征在于，所述根据前一时刻增益矩阵对前一时刻航迹点的状态预测值以及计算前一时刻航迹点的状态预测协方差矩阵进行修正，得到修正后最终反向预测的前一时刻航迹点的状态预测值以及前一时刻航迹点的状态预测协方差矩阵通过下述方程式计算得到：
X^(k-1/k-1)=X^(k-1/k)+K[Z(k)-HX^(k-1/k)]]]>
P(k-1/k-1)＝(I-KH)P(k-1/k)
其中，为修正后最终预测的前一时刻航迹点的状态预测值，P(k-1/k-1)为修正后最终的前一时刻航迹点的状态预测协方差矩阵，Z(k)为航迹测量值，I为3阶单位矩阵。

8.  一种基于ADS-B系统航迹点的预测系统，其特征在于，包括预处理单元、判断单元、计算预测单元；
所述预处理单元将接收到的大地坐标系航迹点转化为直角坐标系航迹点并按接收时间对收到的航迹点进行排列；
所述判断单元判断接收到的航迹点是否为初始航迹点；
所述计算预测单元在判断单元判断接收到的航迹点为初始航迹点后，根据航迹点的状态信息正向预测下一时刻航迹点的状态预测值以及计算下一时刻航迹点的状态预测协方差矩阵，并正向预测下一时刻航迹点的测量预测值以及计算下一时刻航迹点的残差向量协方差矩阵；根据下一时刻航迹点的测量预测值以及下一时刻航迹点的残差向量协方差矩阵计算下一时刻航迹点的衰减因子以及根据该衰减因子计算下一时刻滤波器增益；以下一时刻航迹点的测量预测值、下一时刻航迹点的残差向量协方差矩阵以及下一时刻滤波器增益对下一时刻航迹点的状态预测值以及下一时刻航迹点的状态预测协方差矩阵进行修正，得到最终正向预测的下一时刻航迹点的状态预测值以及最终下一时刻航迹点的状态预测协方差矩阵。

9.  根据权利要求8所述的基于ADS-B系统航迹点的预测系统，其特征在于，所述计算预测单元在判断单元判断接收到的航迹点不是初始航迹点后，根据航迹点的状态信息反向预测前一时刻航迹点的状态预测值以及计算前一时刻航迹点的状态预测协方差矩阵，并计算前一时刻增益矩阵；根据前一时刻增益矩阵对前一时刻航迹点的状态预测值以及计算前一时刻航迹点的状态预测协方差矩阵进行修正，得到修正后最终反向预测的前一时刻航迹点的状态预测值以及前一时刻航迹点的状态预测协方差矩阵。

10.  根据权利要求9所述的基于ADS-B系统航迹点的预测系统，其特征在于，所述计算预测单元通过下述方程式进行正向预测计算：
X^(k+1/k)=Φ(k)X^(k/k)+U(k)a‾(k)]]>
P(k+1/k)＝Φ(k)P(k/k)Φ^T(k)+Q(k)
Y^(k+1/k)=H(k)X^(k+1/k)]]>
S(k+1)＝H(k)P(k+1/k)H^T(k)+R(k)
Y~(k+1)=|Y(k+1)-Y^(k+1/k)|]]>
r(k+1)=1|Y~(k+1)|≤3S(k+1){3S(k+1)|Y~(k+1)|}1/2|Y~(k+1)|>3S(k+1)]]>
K(k+1)＝r(k+1)P(k+1)H^T(k)S^-1(k+1)
X^(k+1/k+1)=X^(k+1/k)+K(k+1)[Y(k+1)-H(k)X^(k+1/k)]]]>
P(k+1/k+1)＝P(k+1/k)-K(k+1)S(k+1)K^T(k+1)
其中，为下一时刻航迹点的状态预测值，P(k+1/k)为下一时刻航迹点的状态预测协方差矩阵；为当前时刻航迹点的状态预测值，P(k/k)为当前时刻航迹点的状态预测协方差矩阵；Φ(k)为当前时刻正向状态转移矩阵、U(k)为当前时刻控制矩阵、为当前时刻机动加速度均值、Q(k)为当前时刻正向噪声方差矩阵、Φ^T(k)为当前时刻正向状态转移矩阵的转置矩阵；为下一时刻航迹点的测量预测值，Y(k+1)为下一时刻航迹点的测量值,S(k+1)为下一时刻航迹点的残差向量协方差矩阵；H(k)为当前时刻测量矩阵、H^T(k)为当前时刻测量矩阵的转置矩阵、R(k)为当前时刻观测噪声；r(k+1)为下一时刻航迹点的衰减因子，K(k+1)为下一时刻滤波器增益；为新息；为修正后最终预测的下一时刻航迹点的状态预测值，P(k+1/k+1)为修正后最终的下一时刻航迹点的状态预测协方差矩阵；
和/或
通过下述方程式进行反向预测计算：
X^(k-1/k)=FX^(k/k)]]>
P(k-1/k)＝FP(k/k)F^T+Q
K＝P(k-1/k)H^T/(HP(k-1/k)H^T+R)
X^(k-1/k-1)=X^(k-1/k)+K[Z(k)-HX^(k-1/k)]]]>
P(k-1/k-1)＝(I-KH)P(k-1/k)
其中，为前一时刻航迹点的状态预测值，P(k-1/k)为前一时刻航迹点的状态预测协方差矩阵；为当前时刻航迹点的状态预测值，P(k/k)为当前时刻航迹点的状态预测协方差矩阵；F为反向预测状态转移矩阵、Q为反 向随机噪声方差；K为前一时刻增益矩阵，H^T为测量矩阵的转置矩阵，H为测量矩阵，R为反向观测噪声均值；为修正后最终预测的前一时刻航迹点的状态预测值，P(k-1/k-1)为修正后最终的前一时刻航迹点的状态预测协方差矩阵，Z(k)为航迹测量值，I为3阶单位矩阵。

基于ADS-B系统航迹点的预测方法和系统
技术领域
本发明涉及基于ADS-B的民用航空空中监管领域，特别是涉及一种基于ADS-B系统航迹点的预测方法和系统。
背景技术
ADS-B(Automatic Dependent Surveillance Broadcast，广播式自动相关监视)系统在现实使用过程中，由于信号干扰等因素的影响，航迹点数据往往包含了较大的随机误差，或者因为监控环境周围地形、飞机转弯等方面的原因，出现了部分航迹甚至连续长时间航迹的缺失。因此提升航迹数据滤波性能和完整性问题成为提高ADS-B系统监控性能的关键点。
目前，常用的滤波算法通过平均机动加速度与相应的加速度方差自适应来提高跟踪精度，但其将目标机动频率和最大加速度设定为固定值，不符合机动目标的性质，且当前加速度与最大加速度差值较大时，对非机动或者弱机动目标跟踪时精度较低，不适合民航飞机非机动情况较多的跟踪要求。还有些滤波算法在上述的算法中加入了机动检测，并通过模糊函数实现机动频率和最大加速度的模糊控制，但这些算法存在实时性差、自适应不理想、预测精度不高等问题。
发明内容
基于此，有必要针对背景技术中现有的航迹点跟踪方法存在实时性差、自适应不理想造成获取的航迹点数据精确度低、丢失严重的问题，提供一种基于ADS-B系统航迹点的预测方法，能够精确地、及时地对航迹点进行预测，并且自适应性高，进而获得完整、准确的航迹点数据。
为实现上述目的，本发明实施例采用的技术方案如下：
一种基于ADS-B系统航迹点的预测方法，包括步骤：
将接收到的大地坐标系航迹点转化为直角坐标系航迹点并按接收时间对收 到的航迹点进行排列；
判断接收到的航迹点是否为初始航迹点；
若是，根据航迹点的状态信息正向预测下一时刻航迹点的状态预测值以及计算下一时刻航迹点的状态预测协方差矩阵，并正向预测下一时刻航迹点的测量预测值以及计算下一时刻航迹点的残差向量协方差矩阵；
根据下一时刻航迹点的测量预测值以及下一时刻航迹点的残差向量协方差矩阵计算下一时刻航迹点的衰减因子以及根据该衰减因子计算下一时刻滤波器增益；
以下一时刻航迹点的测量预测值、下一时刻航迹点的残差向量协方差矩阵以及下一时刻滤波器增益对下一时刻航迹点的状态预测值以及下一时刻航迹点的状态预测协方差矩阵进行修正，得到最终正向预测的下一时刻航迹点的状态预测值以及最终下一时刻航迹点的状态预测协方差矩阵。
根据上述的一种基于ADS-B系统航迹点的预测方法，本发明的方案还提供一种基于ADS-B系统航迹点的预测系统，包括预处理单元、判断单元、计算预测单元；
所述预处理单元将接收到的大地坐标系航迹点转化为直角坐标系航迹点并按接收时间对收到的航迹点进行排列；
所述判断单元判断接收到的航迹点是否为初始航迹点；
所述计算预测单元在判断单元判断接收到的航迹点为初始航迹点后，根据航迹点的状态信息正向预测下一时刻航迹点的状态预测值以及计算下一时刻航迹点的状态预测协方差矩阵，并正向预测下一时刻航迹点的测量预测值以及计算下一时刻航迹点的残差向量协方差矩阵；根据下一时刻航迹点的测量预测值以及下一时刻航迹点的残差向量协方差矩阵计算下一时刻航迹点的衰减因子以及根据该衰减因子计算下一时刻滤波器增益；以下一时刻航迹点的测量预测值、下一时刻航迹点的残差向量协方差矩阵以及下一时刻滤波器增益对下一时刻航迹点的状态预测值以及下一时刻航迹点的状态预测协方差矩阵进行修正，得到最终正向预测的下一时刻航迹点的状态预测值以及最终下一时刻航迹点的状态预测协方差矩阵。
根据本发明方案，先将接收到的大地坐标系航迹点转化为直角坐标系航迹点并按接收时间对收到的航迹点进行排列，接着判断接收到的航迹点是否为初始航迹点，若是，则根据航迹点的状态信息正向预测下一时刻航迹点的状态预测值以及计算下一时刻航迹点的状态预测协方差矩阵，并正向预测下一时刻航迹点的测量预测值以及计算下一时刻航迹点的残差向量协方差矩阵；根据下一时刻航迹点的测量预测值以及下一时刻航迹点的残差向量协方差矩阵计算下一时刻航迹点的衰减因子以及根据该衰减因子计算下一时刻滤波器增益；最后修正航迹点的状态预测值以及下一时刻航迹点的状态预测协方差矩阵，进而精确地、及时地得到最终正向预测的下一时刻航迹点的状态预测值以及最终下一时刻航迹点的状态预测协方差矩阵。
附图说明
图1为本发明基于ADS-B系统航迹点的预测方法实施例流程图；
图2本发明基于ADS-B系统航迹点的预测方法仿真对比图；
图3为图2里面标识1的局部放大图；
图4为图2里面标识2的局部放大图；
图5为本发明基于ADS-B系统航迹点的预测系统结构图。
具体实施方式
为使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步的详细说明。应当理解，此处所描述的具体实施方式仅仅用以解释本发明，并不限定本发明的保护范围。
请参阅图1，为本发明基于ADS-B系统航迹点的预测方法实施例流程图：
步骤S101：将接收到的大地坐标系航迹点转化为直角坐标系航迹点并按接收时间对收到的航迹点进行排列；
具体地，所述大地坐标系航迹点是指采用大地经、纬度和大地高来描述航迹点的空间位置，将大地坐标系航迹点通过现有的投影变换方法转化为直角坐标系航迹点，进而以接收时间对收到的航迹点进行排列。
步骤S102：判断接收到的航迹点是否为初始航迹点；
具体地，所述判断接收到的航迹点是否为初始航迹点是指在接收到3个或3个以上的连续有效的航迹点时，即判断接收到的航迹点为初始航迹点；而且，若在监测的过程中连续有效的航迹点缺失大于3个时，重新判断接收到的航迹点是否为初始航迹点；应当指出，也可以设置在接收到4个连续有效的航迹点时，判断接收到的航迹点为初始航迹点，即设置所述接收连续有效的航迹点的个数可以根据具体应用进行设定。
步骤S103：若是，根据航迹点的状态信息正向预测下一时刻航迹点的状态预测值以及计算下一时刻航迹点的状态预测协方差矩阵，并正向预测下一时刻航迹点的测量预测值以及计算下一时刻航迹点的残差向量协方差矩阵；
所述下一时刻航迹点的状态预测值是指下一时刻飞行目标的状态信息，包括飞行目标的位置、速度、加速度等；所述下一时刻航迹点的状态预测协方差矩阵主要是体现下一时刻航迹点的状态预测值的误差，通过调整下一时刻航迹点的状态预测协方差矩阵能够减少下一时刻航迹点的状态预测值的误差。
所述下一时刻航迹点的测量预测值指的是下一时刻航迹点的位置信息，所述下一时刻航迹点的残差向量协方差矩阵主要是体现下一时刻航迹点的测量预测值的误差，通过调整下一时刻航迹点的残差向量协方差矩阵能够减少下一时刻航迹点的测量预测值的误差。
具体地，在一个实施例中，所述根据航迹点的状态信息正向预测下一时刻航迹点的状态预测值以及计算下一时刻航迹点的状态预测协方差矩阵，并正向预测下一时刻航迹点的测量预测值以及计算下一时刻航迹点的残差向量协方差矩阵可以通过下述方程式计算得到：
X^(k+1/k)=Φ(k)X^(k/k)+U(k)a‾(k)]]>
P(k+1/k)＝Φ(k)P(k/k)Φ^T(k)+Q(k)
Y^(k+1/k)=H(k)X^(k+1/k)]]>
S(k+1)＝H(k)P(k+1/k)H^T(k)+R(k)
其中，为下一时刻航迹点的状态预测值，P(k+1/k)为下一时刻航迹点的状态预测协方差矩阵，为下一时刻航迹点的测量预测值， S(k+1)为下一时刻航迹点的残差向量协方差矩阵；
所述航迹点状态信息包括当前时刻航迹点的状态预测值当前时刻航迹点的状态预测协方差矩阵P(k/k)、当前时刻正向状态转移矩阵Φ(k)、当前时刻控制矩阵U(k)、当前时刻机动加速度均值当前时刻正向噪声方差矩阵Q(k)、当前时刻正向状态转移矩阵的转置矩阵Φ^T(k)、当前时刻测量矩阵H(k)、当前时刻测量矩阵的转置矩阵H^T(k)、当前时刻观测噪声R(k)。
较优地，在一个实施例中，当前时刻观测噪声R(k)的值可以为30，当前时刻测量矩阵H(k)＝[1 0 0]。
进一步地，
当前时刻正向状态转移矩阵Φ(k)=1T1α2(-1+αT+e-αT)011α(1-e-αT)00e-αT]]>
当前时刻控制矩阵U(k)=1α(-T+αT22+1-e-αTα)T-1-e-αTα1-e-αT]]>
其中，α为机动频率，转弯机动频率α＝1/60，逃避机动频率α＝1/20，大气扰动机动频率α＝1；T为采样时间；e为欧拉数；
当前时刻机动加速度均值a‾(k)=x··^(k|k-1)=x··^(k-1|k-1);]]>其中，为k-1时刻的加速度预测值，为k-1时刻的加速度估计值；加速度协方差的表达式为：σa2(k+1)=|x·^(k|k)-x·^(k|k-1)|]]>
其中，为k时刻的速度估计值，为k时刻的速度预测值；
当前时刻噪声方差矩阵
q=q11q12q13q12q22q23q13q23q33]]>
式中
q11=12α5[1-e-2αT+2αT+2α3T33-2α2T2-4αTe-αT]]]>
q12=12α4[e-2αT+1-2e-αT+2αTe-αT-2αT+α2T2]]]>
q13=12α3[1-e-2αT-2αTe-αT]]]>
q22=12α3[4e-αT-3-e-2αT+2αT]]]>
q23=12α2[e-2αT+1-2e-αT]]]>
q33=12α[1-e-2αT]]]>
步骤S104：根据下一时刻航迹点的测量预测值以及下一时刻航迹点的残差向量协方差矩阵计算下一时刻航迹点的衰减因子以及根据该衰减因子计算下一时刻滤波器增益；
具体地，所述根据下一时刻航迹点的测量预测值以及下一时刻航迹点的残差向量协方差矩阵计算下一时刻航迹点的衰减因子以及根据该衰减因子计算下一时刻滤波器增益可以通过下述方程式计算得到：
Y~(k+1)=|Y(k+1)-Y^(k+1/k)|]]>
r(k+1)=1|Y~(k+1)|≤3S(k+1){3S(k+1)|Y~(k+1)|}1/2|Y~(k+1)|>3S(k+1)]]>
K(k+1)＝r(k+1)P(k+1)H^T(k)S^-1(k+1)
其中，为新息，Y(k+1)为下一时刻航迹点的测量值，r(k+1)为下一时刻航迹点的衰减因子，K(k+1)为下一时刻滤波器增益。
步骤S105：以下一时刻航迹点的测量预测值、下一时刻航迹点的残差向量协方差矩阵以及下一时刻滤波器增益对下一时刻航迹点的状态预测值以及下一 时刻航迹点的状态预测协方差矩阵进行修正，得到最终正向预测的下一时刻航迹点的状态预测值以及最终下一时刻航迹点的状态预测协方差矩阵；
具体地，所述以下一时刻航迹点的测量预测值、下一时刻航迹点的残差向量协方差矩阵以及下一时刻滤波器增益对下一时刻航迹点的状态预测值以及下一时刻航迹点的状态预测协方差矩阵进行修正，得到最终正向预测的下一时刻航迹点的状态预测值以及最终下一时刻航迹点的状态预测协方差矩阵可以通过下述方程式计算得到：
X^(k+1/k+1)=X^(k+1/k)+K(k+1)[Y(k+1)-H(k)X^(k+1/k)]]]>
P(k+1/k+1)＝P(k+1/k)-K(k+1)S(k+1)K^T(k+1)
其中，为修正后最终预测的下一时刻航迹点的状态预测值，P(k+1/k+1)为修正后最终的下一时刻航迹点的状态预测协方差矩阵，K^T(k+1)为下一时刻滤波器增益转置矩阵。
可以看到，本发明的基于ADS-B系统航迹点的预测方法，通过加速度和加速度方差的自适应调节能够精确地、及时地对航迹点进行预测，并且自适应性高。
较优地，在一个实施例中，在判断接收到的航迹点不是初始航迹点时，包括步骤：
根据航迹点的状态信息反向预测前一时刻航迹点的状态预测值以及计算前一时刻航迹点的状态预测协方差矩阵，并计算前一时刻增益矩阵；
根据前一时刻增益矩阵对前一时刻航迹点的状态预测值以及计算前一时刻航迹点的状态预测协方差矩阵进行修正，得到修正后最终反向预测的前一时刻航迹点的状态预测值以及前一时刻航迹点的状态预测协方差矩阵。
具体地，所述根据航迹点的状态信息反向预测前一时刻航迹点的状态预测值以及计算前一时刻航迹点的状态预测协方差矩阵，并计算前一时刻增益矩阵通过下述方程式计算得到：
X^(k-1/k)=FX^(k/k)]]>
P(k-1/k)＝FP(k/k)F^T+Q
K＝P(k-1/k)H^T/(HP(k-1/k)H^T+R)
其中，为前一时刻航迹点的状态预测值，P(k-1/k)为前一时刻航迹点的状态预测协方差矩阵，为当前时刻航迹点的状态预测值，P(k/k)为当前时刻航迹点的状态预测协方差矩阵，K为前一时刻增益矩阵，H为测量矩阵，H^T为测量矩阵的转置矩阵，R为反向观测噪声均值；所述航迹点的状态信息包括反向预测状态转移矩阵F、反向随机噪声方差Q。
进一步地，
反向预测状态转移矩阵F=1TT2/201T001]]>
反向随机噪声方差Q＝δ*G，δ为的高斯白噪声方差均值，可以设为10，T为采样时间，G＝[T³/6 T²/2 T]，反向观测噪声均值R＝30，H＝[1 0 0]。
具体地，所述根据前一时刻增益矩阵对前一时刻航迹点的状态预测值以及计算前一时刻航迹点的状态预测协方差矩阵进行修正，得到修正后最终反向预测的前一时刻航迹点的状态预测值以及前一时刻航迹点的状态预测协方差矩阵通过下述方程式计算得到：
X^(k-1/k-1)=X^(k-1/k)+K[Z(k)-HX^(k-1/k)]]]>
P(k-1/k-1)＝(I-KH)P(k-1/k)
其中，为修正后最终预测的前一时刻航迹点的状态预测值，P(k-1/k-1)为修正后最终的前一时刻航迹点的状态预测协方差矩阵，Z(k)为航迹测量值，I为3阶单位矩阵。
通过上述的基于ADS-B系统航迹点的预测方法反向预测前一时刻航迹点，进一步地提高了预测航迹点的精度，而且增大了预测的范围。
为了进一步说明本发明的基于ADS-B系统航迹点的预测方法，建立以下测试模型：
首先，选取国内某一个航线航班的历史测量数据，在MATLAB2008(matrix laboratory，矩阵实验室)建立仿真模型，然后利用上述发明的基于ADS-B系统航迹点的预测方法对该航线航班进行预测。请参阅图2，为本发明基于ADS-B系统航迹点的预测方法仿真对比图；进一步地，请参阅图3，为图2里面标识1 的局部放大图；以及图4，为图2里面标识2的局部放大图。
从图3以及图4中原始航迹测量值与滤波航迹预测值的对比，可以看到本发明基于ADS-B系统航迹点的预测方法能够精确地、及时地对航迹点进行预测，并且自适应性高。
根据上述的一种基于ADS-B系统航迹点的预测方法，本发明还提供一种基于ADS-B系统航迹点的预测系统，请参阅图5，为本发明基于ADS-B系统航迹点的预测系统结构图：
一种基于ADS-B系统航迹点的预测系统，包括预处理单元10、判断单元20、计算预测单元30；
所述预处理单元10将接收到的大地坐标系航迹点转化为直角坐标系航迹点并按接收时间对收到的航迹点进行排列；
所述判断单元20判断接收到的航迹点是否为初始航迹点；
所述计算预测单元30在判断单元20判断接收到的航迹点为初始航迹点后，根据航迹点的状态信息正向预测下一时刻航迹点的状态预测值以及计算下一时刻航迹点的状态预测协方差矩阵，并正向预测下一时刻航迹点的测量预测值以及计算下一时刻航迹点的残差向量协方差矩阵；根据下一时刻航迹点的测量预测值以及下一时刻航迹点的残差向量协方差矩阵计算下一时刻航迹点的衰减因子以及根据该衰减因子计算下一时刻滤波器增益；以下一时刻航迹点的测量预测值、下一时刻航迹点的残差向量协方差矩阵以及下一时刻滤波器增益对下一时刻航迹点的状态预测值以及下一时刻航迹点的状态预测协方差矩阵进行修正，得到最终正向预测的下一时刻航迹点的状态预测值以及最终下一时刻航迹点的状态预测协方差矩阵。
其中，所述大地坐标系航迹点是指采用大地经、纬度和大地高来描述航迹点的空间位置，所述预处理单元10将大地坐标系航迹点通过现有的投影变换方法转化为直角坐标系航迹点，进而以接收时间对收到的航迹点进行排列。
所述判断单元20判断接收到的航迹点是否为初始航迹点是指在接收到3个或3个以上的连续有效的航迹点时，即判断接收到的航迹点为初始航迹点；而且，若在监测的过程中连续有效的航迹点缺失大于3个时，重新判断接收到的 航迹点是否为初始航迹点；应当指出，也可以设置在接收到4个连续有效的航迹点时，判断接收到的航迹点为初始航迹点，即设置所述接收连续有效的航迹点的个数可以根据具体应用进行设定。
所述下一时刻航迹点的状态预测值是指下一时刻飞行目标的状态信息，包括飞行目标的位置、速度、加速度等；所述下一时刻航迹点的状态预测协方差矩阵主要是体现下一时刻航迹点的状态预测值的误差，通过调整下一时刻航迹点的状态预测协方差矩阵能够减少下一时刻航迹点的状态预测值的误差。
所述下一时刻航迹点的测量预测值指的是下一时刻航迹点的位置信息，所述下一时刻航迹点的残差向量协方差矩阵主要是体现下一时刻航迹点的测量预测值的误差，通过调整下一时刻航迹点的残差向量协方差矩阵能够减少下一时刻航迹点的测量预测值的误差。
具体地，所述计算预测单元30通过下述方程式进行正向预测计算：
X^(k+1/k)=Φ(k)X^(k/k)+U(k)a‾(k)]]>
P(k+1/k)＝Φ(k)P(k/k)Φ^T(k)+Q(k)
Y^(k+1/k)=H(k)X^(k+1/k)]]>
S(k+1)＝H(k)P(k+1/k)H^T(k)+R(k)
Y~(k+1)=|Y(k+1)-Y^(k+1/k)|]]>
r(k+1)=1|Y~(k+1)|≤3S(k+1){3S(k+1)|Y~(k+1)|}1/2|Y~(k+1)|>3S(k+1)]]>
K(k+1)＝r(k+1)P(k+1)H^T(k)S^-1(k+1)
X^(k+1/k+1)=X^(k+1/k)+K(k+1)[Y(k+1)-H(k)X^(k+1/k)]]]>
P(k+1/k+1)＝P(k+1/k)-K(k+1)S(k+1)K^T(k+1)
其中，为下一时刻航迹点的状态预测值，P(k+1/k)为下一时刻航迹点的状态预测协方差矩阵；为当前时刻航迹点的状态预测值，P(k/k)为当前时刻航迹点的状态预测协方差矩阵；Φ(k)为当前时刻正向状态转移矩阵、U(k)为当前时刻控制矩阵、为当前时刻机动加速度均值、Q(k)为当前时刻正 向噪声方差矩阵、Φ^T(k)为当前时刻正向状态转移矩阵的转置矩阵；为下一时刻航迹点的测量预测值，Y(k+1)为下一时刻航迹点的测量值,S(k+1)为下一时刻航迹点的残差向量协方差矩阵；H(k)为当前时刻测量矩阵、H^T(k)为当前时刻测量矩阵的转置矩阵、R(k)为当前时刻观测噪声；r(k+1)为下一时刻航迹点的衰减因子，K(k+1)为下一时刻滤波器增益；为新息；为修正后最终预测的下一时刻航迹点的状态预测值，P(k+1/k+1)为修正后最终的下一时刻航迹点的状态预测协方差矩阵；
所述新息是指判断航迹点是否为野值的一个参数；进一步地，若Y(k+1)的值为0，则可以通过野值判定方式判断其Y(k+1)为缺失的航迹点，并对其进行航迹的预测。
较优地，在一个实施例中，当前时刻观测噪声R(k)的值可以为30，当前时刻测量矩阵H(k)＝[1 0 0]。
进一步地，
当前时刻正向状态转移矩阵Φ(k)=1T1α2(-1+αT+e-αT)011α(1-e-αT)00e-αT]]>
当前时刻控制矩阵U(k)=1α(-T+αT22+1-e-αTα)T-1-e-αTα1-e-αT]]>
其中，α为机动频率，转弯机动频率α＝1/60，逃避机动频率α＝1/20，大气扰动机动频率α＝1；T为采样时间；e为欧拉数；
当前时刻机动加速度均值a‾(k)=x··^(k|k-1)=x··^(k-1|k-1);]]>其中，为k-1时刻的加速度预测值，为k-1时刻的加速度估计值；加速度协方差的表达式为：σa2(k+1)=|x·^(k|k)-x·^(k|k-1)|]]>
其中，为k时刻的速度估计值，为k时刻的速度预测值；
当前时刻噪声方差矩阵
q=q11q12q13q12q22q23q13q23q33]]>
式中
q11=12α5[1-e-2αT+2αT+2α3T33-2α2T2-4αTe-αT]]]>
q12=12α4[e-2αT+1-2e-αT+2αTe-αT-2αT+α2T2]]]>
q13=12α3[1-e-2αT-2αTe-αT]]]>
q22=12α3[4e-αT-3-e-2αT+2αT]]]>
q23=12α2[e-2αT+1-2e-αT]]]>
q33=12α[1-e-2αT]]]>
可以看到，本发明的基于ADS-B系统航迹点的预测系统，通过加速度和加速度方差的自适应调节能够精确地、及时地对航迹点进行预测，并且自适应性高。
较优地，在一个实施例中，所述计算预测单元30在判断单元20判断接收到的航迹点不是初始航迹点后，根据航迹点的状态信息反向预测前一时刻航迹点的状态预测值以及计算前一时刻航迹点的状态预测协方差矩阵，并计算前一时刻增益矩阵；根据前一时刻增益矩阵对前一时刻航迹点的状态预测值以及计算前一时刻航迹点的状态预测协方差矩阵进行修正，得到修正后最终反向预测的前一时刻航迹点的状态预测值以及前一时刻航迹点的状态预测协方差矩阵。
具体地，所述计算预测单元30通过下述方程式进行反向预测计算：
X^(k-1/k)=FX^(k/k)]]>
P(k-1/k)＝FP(k/k)F^T+Q
K＝P(k-1/k)H^T/(HP(k-1/k)H^T+R)
X^(k-1/k-1)=X^(k-1/k)+K[Z(k)-HX^(k-1/k)]]]>
P(k-1/k-1)＝(I-KH)P(k-1/k)
其中，为前一时刻航迹点的状态预测值，P(k-1/k)为前一时刻航迹点的状态预测协方差矩阵；为当前时刻航迹点的状态预测值，P(k/k)为当前时刻航迹点的状态预测协方差矩阵；F为反向预测状态转移矩阵、Q为反向随机噪声方差；K为前一时刻增益矩阵，H^T为测量矩阵的转置矩阵，H为测量矩阵，R为反向观测噪声均值；为修正后最终预测的前一时刻航迹点的状态预测值，P(k-1/k-1)为修正后最终的前一时刻航迹点的状态预测协方差矩阵，Z(k)为航迹测量值，I为3阶单位矩阵。
进一步地，
反向预测状态转移矩阵F=1TT2/201T001]]>
反向随机噪声方差Q＝δ*G，δ为的高斯白噪声方差均值，设为10，T为采样时间G＝[T³/6 T²/2 T]，反向观测噪声均值R＝30，H＝[1 0 0]。
通过上述的基于ADS-B系统航迹点的预测系统反向预测前一时刻航迹点，进一步地提高了预测航迹点的精度，而且增大了预测的范围。
以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。