CN104180800A - 基于ads-b***航迹点的修正方法和*** - Google Patents

Abstract

一种基于ADS-B***航迹点的修正方法和***，其方法包括步骤：对航迹点按接收时间进行排列；根据航迹点的状态信息预测当前时刻航迹点的状态预测值以及计算当前时刻航迹点的状态预测协方差矩阵；根据航迹点的状态信息预测当前时刻航迹点的测量预测值以及计算当前时刻航迹点的残差向量协方差矩阵；识别当前时刻航迹点是否为野值，并根据当前时刻航迹点是否为野值确定当前时刻航迹点的衰减因子；根据当前时刻航迹点的衰减因子计算当前时刻滤波器增益；修正当前时刻航迹点的状态预测值以及当前航迹点的状态预测协方差矩阵。根据本发明方案，能够有效识别航迹点中的野值，实现对航迹点中的野值进行实时动态的修正，提高航迹点数据的准确度。

Description

基于ADS-B***航迹点的修正方法和***

技术领域

本发明涉及基于ADS-B的民用航空空中监管领域，特别是涉及一种基于ADS-B***航迹点的修正方法和***。

背景技术

ADS-B(Automatic Dependent Surveillance Broadcast，广播式自动相关监视)监控设备在民用航空空中监管的使用过程中，由于噪声干扰、窜扰、混扰、多径干扰等因素的影响，得到的航迹点数据往往包含了较大的随机误差，造成航迹点数据严重偏离目标真实值，即成为野值。

目前，传统的航迹点野值处理方法是通过模糊控制函数判断新的航迹值是否为野值，如果不是野值，则不对航迹估计值进行调整；如果野值较小，则调整增益实现对野值的修正；如果野值较大，则对野值进行剔除。然而，这种野值的处理方法，存在野值判断不灵敏，无法实现对野值实时动态的修正；对于连续成片的野值，还有可能存在野值的判定和处理误差较大的情况，且对于较大的野值，直接将该野值剔除，无法实现航迹点野值的修正，而且降低了航迹点数据的准确度。

发明内容

基于此，有必要针对背景技术中，现有的航迹点野值处理方法出现野值判断不灵敏，无法实现对野值实时动态的修正，降低了航迹点数据的完整性的问题，提供一种基于ADS-B***航迹点的修正方法，能够有效识别航迹点中的野值，实现对航迹点中的野值进行实时动态的修正，提高航迹点数据的准确度。

为实现上述目的，本发明实施例采用的技术方案如下：

一种基于ADS-B***航迹点的修正方法，包括步骤：

对航迹点按接收时间进行排列；

根据航迹点的状态信息预测当前时刻航迹点的状态预测值以及计算当前时刻航迹点的状态预测协方差矩阵；

根据航迹点的状态信息预测当前时刻航迹点的测量预测值以及计算当前时刻航迹点的残差向量协方差矩阵；

识别当前时刻航迹点是否为野值，并根据当前时刻航迹点是否为野值确定当前时刻航迹点的衰减因子；

根据当前时刻航迹点的衰减因子计算当前时刻滤波器增益；

根据当前时刻航迹点的测量预测值、当前航迹点的残差向量协方差矩阵以及当前时刻滤波器增益修正当前时刻航迹点的状态预测值以及当前航迹点的状态预测协方差矩阵。

根据上述的一种基于ADS-B***航迹点的修正方法，本发明提供一种基于ADS-B***航迹点的修正***，包括预处理单元、计算修正单元；

所述预处理单元对航迹点按接收时间进行排列；

所述计算修正单元根据航迹点的状态信息预测当前时刻航迹点的状态预测值以及计算当前时刻航迹点的状态预测协方差矩阵；根据航迹点的状态信息预测当前时刻航迹点的测量预测值以及计算当前时刻航迹点的残差向量协方差矩阵；识别当前时刻航迹点是否为野值，并根据当前时刻航迹点是否为野值确定当前时刻航迹点的衰减因子；根据当前时刻航迹点的衰减因子计算当前时刻滤波器增益；根据当前时刻航迹点的测量预测值、当前航迹点的残差向量协方差矩阵以及当前时刻滤波器增益修正当前时刻航迹点的状态预测值以及当前航迹点的状态预测协方差矩阵。

根据本发明方案，先对航迹点按接收时间进行排列；然后根据航迹点的状态信息预测当前时刻航迹点的状态预测值以及计算当前时刻航迹点的状态预测协方差矩阵；并预测当前时刻航迹点的测量预测值以及计算当前时刻航迹点的残差向量协方差矩阵；进而识别当前时刻航迹点是否为野值，并根据当前时刻航迹点是否为野值确定当前时刻航迹点的衰减因子；再进一步根据当前时刻航迹点的衰减因子计算当前时刻滤波器增益；最后修正当前时刻航迹点的状态预测值以及当前航迹点的状态预测协方差矩阵。进而达到有效识别航迹点中的野值，实现对航迹点中的野值进行实时动态的修正，提高航迹点数据的准确度。

附图说明

图1为本发明基于ADS-B***航迹点的修正方法的实施例流程图；

图2为本发明基于ADS-B***航迹点的修正方法的一维仿真模型图；

图3为本发明基于ADS-B***航迹点的修正方法的一维仿真修正结果图；

图4为本发明基于ADS-B***航迹点的修正***结构图。

具体实施方式

为使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步的详细说明。应当理解，此处所描述的具体实施方式仅仅用以解释本发明，并不限定本发明的保护范围。

请参阅图1，为本发明基于ADS-B***航迹点的修正方法的实施例流程图：

步骤S101：对航迹点按接收时间进行排列；

具体地，所述按接收时间排列航迹点是指将接收到的航迹点由大地坐标系转化为直角坐标系，再将航迹点的数据按照接收的时间是先后顺序进行排列。

步骤S102：根据航迹点的状态信息预测当前时刻航迹点的状态预测值以及计算当前时刻航迹点的状态预测协方差矩阵；

所述当前时刻航迹点的状态预测值是指当前时刻飞行目标的状态信息，包括飞行目标的位置、速度、加速度等；所述当前时刻航迹点的状态预测协方差矩阵主要是体现当前时刻航迹点的状态预测值的误差，通过调整当前时刻航迹点的状态预测协方差矩阵能够减少当前时刻航迹点的状态预测值的误差。

具体地，在一个实施例中，所述根据航迹点的状态信息预测当前时刻航迹点的状态预测值以及计算当前航迹点的状态预测协方差矩阵通过下述方程式计算得到：

$\hat{X}(k/k-1)=\mathrm{\Φ}(k/k-1)\hat{X}(k-1/k-1)+U(k-1)\stackrel{\‾}{a}(k-1)$

P(k/k-1)＝Φ(k-1)P(k-1/k-1)Φ^T(k-1)+Q(k-1)

其中，为当前时刻航迹点的状态预测值，P(k/k-1)为当前时刻航迹点的状态预测协方差矩阵，为前一时刻航迹点的状态预测值，P(k-1/k-1)为前一时刻航迹点的状态预测协方差矩阵；

所述航迹点状态信息包括前一时刻状态转移矩阵Φ(k-1)、前一时刻控制矩阵U(k-1)、前一时刻机动加速度均值前一时刻航迹点的状态预测协方差矩阵P(k/k-1)、前一时刻状态转移矩阵的转置矩阵Φ^T(k-1)、前一时刻噪声方差矩阵Q(k-1)。

进一步地：

前一时刻状态转移矩阵 $\mathrm{\Φ}(k-1)=\left[\begin{array}{ccc}1& T& \frac{1}{a^2}(-1+\mathrm{aT}+e^{-\mathrm{aT}})\\ 0& 1& \frac{1}{a}(1-e^{\mathrm{aT}})\\ 0& 0& e^{-\mathrm{aT}}\end{array}\right]$

α为机动频率，转弯机动频率α＝1/60，逃避机动频率α＝1/20，大气扰动机动频率α＝1；T为采样时间；e为欧拉数；

前一时刻控制矩阵 $U(k-1)=\left[\begin{array}{c}\frac{1}{a}(-T+\frac{{\mathrm{aT}}^2}{2}+\frac{1-e^{-\mathrm{aT}}}{a})\\ T-\frac{1-e^{-\mathrm{aT}}}{a}\\ 1-e^{-\mathrm{aT}}\end{array}\right]$

当前时刻噪声方差矩阵

$q=\left[\begin{array}{ccc}{q}_{11}& q_{12}& q_{13}\\ q_{12}& q_{22}& q_{23}\\ q_{13}& q_{23}& q_{33}\end{array}\right]$

式中

$q_{11}=\frac{1}{{2a}^5}[1-e^{-2\mathrm{aT}}+2\mathrm{aT}+\frac{{2a}^3{T}^3}{3}-{2a}^2{T}^2-{4\mathrm{aTe}}^{-\mathrm{aT}}]$

$q_{12}=\frac{1}{{2a}^4}[e^{-2\mathrm{aT}}+1-{2e}^{-\mathrm{aT}}+{2\mathrm{aTe}}^{-\mathrm{aT}}-2\mathrm{aT}+a^2{T}^2]$

$q_{13}=\frac{1}{{2a}^3}[1-e^{-2\mathrm{aT}}-{2\mathrm{aTe}}^{-\mathrm{aT}}]$

$q_{22}=\frac{1}{{2a}^3}[{4e}^{-\mathrm{aT}}-3-e^{-2\mathrm{aT}}+2\mathrm{aT}]$

$q_{23}=\frac{1}{{2a}^2}[e^{-2\mathrm{aT}}+1-{2e}^{-\mathrm{aT}}]$

$q_{33}=\frac{1}{2a}[1-e^{-2\mathrm{aT}}]$

${\mathrm{\&sigma;}}_a^2(k-1)=\left\{\begin{array}{cc}\frac{4-\mathrm{\&pi;}}{\mathrm{\&pi;}}{[a_{\max }-\stackrel{\&OverBar;}{a}(k-1)]}^2& \stackrel{\&OverBar;}{a}(k-1)\&GreaterEqual;0\\ \frac{4-\mathrm{\&pi;}}{\mathrm{\&pi;}}{[a_{-\max }+\stackrel{\&OverBar;}{a}(k-1)]}^2& \stackrel{\&OverBar;}{a}(k-1)<0\end{array}\right.$

$\stackrel{\&OverBar;}{a}(k-1)=\widehat{\stackrel{\&CenterDot;\&CenterDot;}{x}}(k-1|k-2$

为k-2时刻的加速度预测值。

步骤S103：根据航迹点的状态信息预测当前时刻航迹点的测量预测值以及计算当前时刻航迹点的残差向量协方差矩阵；

所述当前时刻航迹点的测量预测值指的是当前时刻航迹点的位置信息，所述当前时刻航迹点的残差向量协方差矩阵主要是体现当前时刻航迹点的测量预测值的误差，通过调整当前时刻航迹点的残差向量协方差矩阵能够减少当前时刻航迹点的测量预测值的误差。

具体地，在一个实施例中，所述根据航迹点的状态信息预测当前时刻航迹点的测量预测值以及计算当前时刻航迹点的残差向量协方差矩阵通过下述方程式计算得到：

$\hat{Y}(k/k-1)=H\left(k\right)\hat{X}(k/k-1)$

S(k)＝H(k)P(k/k-1)H^T(k)+R(k)

其中，为当前时刻航迹点的测量预测值，S(k)为当前时刻航迹点的残差向量协方差矩阵；

所述航迹点状态信息还包括当前时刻测量矩阵H(k)、当前时刻测量矩阵的转置矩阵H^T(k)、当前时刻观测噪声R(k)。

较优地，在一个实施例中，当前时刻观测噪声R(k)的值可以为30，当前时刻测量矩阵H(k)＝[100]。

步骤S104：识别当前时刻航迹点是否为野值，并根据当前时刻航迹点是否为野值确定当前时刻航迹点的衰减因子；

所述野值是指该航迹点数据严重偏离目标真实值，即为野值；由于要对野值进行修正，所以先判定该航迹点是否为野值，再确定该航迹点的衰减因子。

具体地，在一个实施例中，所述识别当前时刻航迹点是否为野值，并根据当前时刻航迹点是否为野值确定当前时刻航迹点的衰减因子通过下述方程式计算得到：

$\tilde{Y}\left(k\right)=|Y\left(k\right)-\hat{Y}(k/k-1)|\&le;d\sqrt{S\left(k\right)}$

$r\left(k\right)=\left\{\begin{array}{ccc}1& \left|\tilde{Y}\left(k\right)\right|\&le;& d\sqrt{S\left(k\right)}\\ {\left\{\frac{3\sqrt{S\left(k\right)}}{\left|\tilde{Y}\left(k\right)\right|}\right\}}^{1/2}& \left|\tilde{Y}\left(k\right)\right|>& d\sqrt{S\left(k\right)}\end{array}\right.$

其中,为当前时刻航迹点的新息，当时，当前时刻航迹点不是野值，当时，当前时刻航迹点为野值；Y(k)为当前时刻航迹点的测量值，为当前时刻航迹点的测量预测值，d设置为3，r(k)为当前时刻航迹点的衰减因子。

应当指出的是，d的值也可以设置成其他大于零的自然数；在本发明方案中，d设置为3是一个识别当前时刻航迹点是否为野值较优的实施例。

步骤S105：根据当前时刻航迹点的衰减因子计算当前时刻滤波器增益；

所述当前时刻滤波器增益是指用于修正当前时刻航迹点的状态预测值以及当前航迹点的状态预测协方差矩阵的增益值；

具体地，所述根据当前时刻航迹点的衰减因子计算当前时刻滤波器增益通过下述方程式计算得到：

K(k)＝r(k)P(k)H^T(k)S^-1(k)

其中，K(k)为当前时刻滤波器增益，)S^-1(k)为当前时刻航迹点的残差向量协方差矩阵的逆矩阵。

通过调整滤波器衰减因子的大小(即r(k))，改变滤波器增益，实现了野值自适应修正。由上述方程式可以得到，野值的修正是通过衰减因子r(k)对增益进行调整实现的，调整范围为[0,1]。即是野值与预测偏离得越大，限制作用也越大，当非常大时，增益几乎接近于零，当处于正常值时，增益的值为1。

步骤106：根据当前时刻航迹点的测量预测值、当前航迹点的残差向量协方差矩阵以及当前时刻滤波器增益修正当前时刻航迹点的状态预测值以及当前航迹点的状态预测协方差矩阵；

具体地，所述根据当前时刻航迹点的测量预测值、当前航迹点的残差向量协方差矩阵以及当前时刻滤波器增益修正当前时刻航迹点的状态预测值以及当前航迹点的状态预测协方差矩阵通过下述方程式计算得到：

$\hat{X}(k/k)=\hat{X}(k/k-1)+K\left(k\right)[Y\left(k\right)-\hat{Y}(k/k-1)]$

P(k/k)＝P(k/k-1)-K(k)S(k)K^T(k)

其中，为修正后的当前时刻航迹点的状态预测值，P(k/k)为修正后的当前航迹点的状态预测协方差矩阵，K^T(k)为当前时刻滤波器增益转置矩阵。

可以看到，本发明的基于ADS-B***航迹点的修正方法，能够有效识别航迹点中的野值，实现对航迹点中的野值进行实时动态的修正，提高航迹点数据的准确度。

较优地，在一个实施例中，在对航迹点按接收时间进行排列之后，预测当前时刻航迹点的状态预测值以及计算当前时刻航迹点的状态预测协方差矩阵之前还可以包括步骤：

将每两个航迹点的时间间隔与滤波器的采样时间间隔进行比较；

若时间间隔不一致，则判断该两个航迹点之间存在航迹点丢失；

用野值对丢失的航迹点进行填充。

具体地，对丢失的航迹点可以数值为0的野值进行填充，通过找出丢失的航迹点，用野值进行填充，能够保证航迹点数据的完整，在野值被修正后可以获得更准确的航迹点数据。

为了进一步说明本发明的基于ADS-B***航迹点的修正方法，下面建立一个一维仿真模型：

首先，在MATLAB2008(matrix laboratory，矩阵实验室)建立飞机运动方程，然后计算出一维坐标下的航迹点，并对该些航迹点进行随机噪声和随机缺失处理，形成模拟的一维飞机航迹；具体地，请参阅图2，为本发明基于ADS-B***航迹点的修正方法的一维仿真模型图，模拟飞机的航迹如下：

设飞机的初始位置为10000，[1，20]秒做速度为600m/s匀速运动，[20，40]秒做加速度为-2g的减速运动，[40，80]秒做加速度为g(自然加速度)的加速运动，[80，100]秒做匀速运动。在[10，30]中20％的航迹点加入200倍随机噪声，同时在[30，100]航迹点随机引发20％航迹缺失，航迹缺失的点用极大野值填充(例子中用0填充)。

接着，在MATLAB2008中利用本发明的基于ADS-B***航迹点的修正方法对上述模拟飞机的航迹进行修正，得到仿真结果；请参阅图3，为本发明基于ADS-B***航迹点的修正方法的一维仿真修正结果图；仿真结果显示了本发明基于ADS-B***航迹点的修正方法能够有效识别航迹点中的野值，实现对航迹点中的野值进行实时动态的修正，提高航迹点数据的准确度。

根据上述的一种基于ADS-B***航迹点的修正方法，本发明还提供一种基于ADS-B***航迹点的修正***，请参阅图4，为本发明基于ADS-B***航迹点的修正***结构图：

一种基于ADS-B***航迹点的修正***，包括预处理单元10、计算修正单元20；

所述预处理单元10对航迹点按接收时间进行排列；

所述计算修正单元20根据航迹点的状态信息预测当前时刻航迹点的状态预测值以及计算当前时刻航迹点的状态预测协方差矩阵；根据航迹点的状态信息预测当前时刻航迹点的测量预测值以及计算当前时刻航迹点的残差向量协方差矩阵；识别当前时刻航迹点是否为野值，并根据当前时刻航迹点是否为野值确定当前时刻航迹点的衰减因子；根据当前时刻航迹点的衰减因子计算当前时刻滤波器增益；根据当前时刻航迹点的测量预测值、当前航迹点的残差向量协方差矩阵以及当前时刻滤波器增益修正当前时刻航迹点的状态预测值以及当前航迹点的状态预测协方差矩阵。

具体地，所述按接收时间排列航迹点是指将接收到的航迹点由大地坐标系转化为直角坐标系，再将航迹点的数据按照接收的时间是先后顺序进行排列。

具体地，所述计算修正单元20通过下述方程式进行计算：

$\hat{X}(k/k-1)=\mathrm{\&Phi;}(k/k-1)\hat{X}(k-1/k-1)+U(k-1)\stackrel{\&OverBar;}{a}(k-1)$

P(k/k-1)＝Φ(k-1)P(k-1/k-1)Φ^T(k-1)+Q(k-1)

$\hat{Y}(k/k-1)=H\left(k\right)\hat{X}(k/k-1)$

S(k)＝H(k)P(k/k-1)H^T(k)+R(k)

$\tilde{Y}\left(k\right)=|Y\left(k\right)-\hat{Y}(k/k-1)|\&le;d\sqrt{S\left(k\right)}$

$r\left(k\right)=\left\{\begin{array}{ccc}1& \left|\tilde{Y}\left(k\right)\right|\&le;& d\sqrt{S\left(k\right)}\\ {\left\{\frac{3\sqrt{S\left(k\right)}}{\left|\tilde{Y}\left(k\right)\right|}\right\}}^{1/2}& \left|\tilde{Y}\left(k\right)\right|>& d\sqrt{S\left(k\right)}\end{array}\right.$

K(k)＝r(k)P(k)H^T(k)S^-1(k)

$\hat{X}(k/k)=\hat{X}(k/k-1)+K\left(k\right)[Y\left(k\right)-\hat{Y}(k/k-1)]$

P(k/k)＝P(k/k-1)-K(k)S(k)K^T(k)

其中，为当前时刻航迹点的状态预测值，Φ(k-1)为前一时刻状态转移矩阵，U(k-1)为前一时刻控制矩阵，为前一时刻机动加速度均值，为前一时刻航迹点的状态预测值；P(k/k-1)为当前时刻航迹点的状态预测协方差矩阵，Q(k-1)为噪声方差矩阵，P(k-1/k-1)为前一时刻航迹点的状态预测协方差矩阵，Φ^T(k-1)为前一时刻状态转移矩阵的转置矩阵；为当前时刻航迹点的测量预测值，H(k)为当前时刻测量矩阵；S(k)为当前时刻航迹点的残差向量协方差矩阵，R(k)当前时刻观测噪声，H^T(k)为当前时刻测量矩阵的转置矩阵；为当前时刻航迹点的新息，Y(k)为当前时刻航迹点的测量值，为当前时刻航迹点的测量预测值，d设置为3，r(k)为当前时刻航迹点的衰减因子；K(k)为当前时刻滤波器增益，S^-1(k)为当前时刻航迹点的残差向量协方差矩阵的逆矩阵；为修正后的当前时刻航迹点的状态预测值，P(k/k)为修正后的当前航迹点的状态预测协方差矩阵，K^T(k)为当前时刻滤波器增益转置矩阵。

进一步地：

前一时刻状态转移矩阵 $\mathrm{\&Phi;}(k-1)=\left[\begin{array}{ccc}1& T& \frac{1}{a^2}(-1+\mathrm{aT}+e^{-\mathrm{aT}})\\ 0& 1& \frac{1}{a}(1-e^{\mathrm{aT}})\\ 0& 0& e^{-\mathrm{aT}}\end{array}\right]$

α为机动频率，转弯机动频率α＝1/60，逃避机动频率α＝1/20，大气扰动机动频率α＝1；T为采样时间；e为欧拉数；

前一时刻控制矩阵 $U(k-1)=\left[\begin{array}{c}\frac{1}{a}(-T+\frac{{\mathrm{aT}}^2}{2}+\frac{1-e^{-\mathrm{aT}}}{a})\\ T-\frac{1-e^{-\mathrm{aT}}}{a}\\ 1-e^{-\mathrm{aT}}\end{array}\right]$

当前时刻噪声方差矩阵

式中

$q_{11}=\frac{1}{{2a}^5}[1-e^{-2\mathrm{aT}}+2\mathrm{aT}+\frac{{2a}^3{T}^3}{3}-{2a}^2{T}^2-{4\mathrm{aTe}}^{-\mathrm{aT}}]$

$q_{12}=\frac{1}{{2a}^4}[e^{-2\mathrm{aT}}+1-{2e}^{-\mathrm{aT}}+{2\mathrm{aTe}}^{-\mathrm{aT}}-2\mathrm{aT}+a^2{T}^2]$

$q_{13}=\frac{1}{{2a}^3}[1-e^{-2\mathrm{aT}}-{2\mathrm{aTe}}^{-\mathrm{aT}}]$

$q_{22}=\frac{1}{{2a}^3}[{4e}^{-\mathrm{aT}}-3-e^{-2\mathrm{aT}}+2\mathrm{aT}]$

$q_{23}=\frac{1}{{2a}^2}[e^{-2\mathrm{aT}}+1-{2e}^{-\mathrm{aT}}]$

$q_{33}=\frac{1}{2a}[1-e^{-2\mathrm{aT}}]$

${\mathrm{\&sigma;}}_a^2(k-1)=\left\{\begin{array}{cc}\frac{4-\mathrm{\&pi;}}{\mathrm{\&pi;}}{[a_{\max }-\stackrel{\&OverBar;}{a}(k-1)]}^2& \stackrel{\&OverBar;}{a}(k-1)\&GreaterEqual;0\\ \frac{4-\mathrm{\&pi;}}{\mathrm{\&pi;}}{[a_{-\max }+\stackrel{\&OverBar;}{a}(k-1)]}^2& \stackrel{\&OverBar;}{a}(k-1)<0\end{array}\right.$

$\stackrel{\&OverBar;}{a}(k-1)=\widehat{\stackrel{\&CenterDot;\&CenterDot;}{x}}(k-1|k-2$

应当指出的是，d的值也可以设置成其他大于零的自然数；在本发明方案中，d设置为3是一个识别当前时刻航迹点是否为野值较优的实施例。

较优地，在一个实施例中，当前时刻观测噪声R(k)的值可以为30，当前时刻测量矩阵H(k)＝[100]。

通过调整滤波器衰减因子的大小(即r(k))，改变滤波器增益，实现了野值自适应修正。由上述方程式可以得到，野值的修正是通过衰减因子r(k)对增益进行调整实现的，调整范围为[0,1]。即是野值与预测偏离得越大，限制作用也越大，当非常大时，增益几乎接近于零，当处于正常值时，增益的值为1。

可以看到，本发明的基于ADS-B***航迹点的修正***，能够有效识别航迹点中的野值，实现对航迹点中的野值进行实时动态的修正，提高航迹点数据的准确度。

较优地，在一个实施例中，本发明的基于ADS-B***航迹点的修正***，还包括填充单元30；

所述填充单元将每两个航迹点的时间间隔与滤波器的采样时间间隔进行比较；若时间间隔不一致，则判断该两个航迹点之间存在航迹点丢失；用野值对丢失的航迹点进行填充。

具体地，所述填充单元对丢失的航迹点可以数值为0的野值进行填充，通过找出丢失的航迹点，用野值进行填充，能够保证航迹点数据的完整，在野值被修正后可以获得更准确的航迹点数据。

以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。

Claims (10)

1.一种基于ADS-B***航迹点的修正方法，其特征在于，包括步骤:

对航迹点按接收时间进行排列；

根据航迹点的状态信息预测当前时刻航迹点的状态预测值以及计算当前时刻航迹点的状态预测协方差矩阵；

根据航迹点的状态信息预测当前时刻航迹点的测量预测值以及计算当前时刻航迹点的残差向量协方差矩阵；

识别当前时刻航迹点是否为野值，并根据当前时刻航迹点是否为野值确定当前时刻航迹点的衰减因子；

根据当前时刻航迹点的衰减因子计算当前时刻滤波器增益；

根据当前时刻航迹点的测量预测值、当前航迹点的残差向量协方差矩阵以及当前时刻滤波器增益修正当前时刻航迹点的状态预测值以及当前航迹点的状态预测协方差矩阵。

2.根据权利要求1所述的基于ADS-B***航迹点的修正方法，其特征在于，所述根据航迹点的状态信息预测当前时刻航迹点的状态预测值以及计算当前航迹点的状态预测协方差矩阵通过下述方程式计算得到：

P(k/k-1)＝Φ(k-1)P(k-1/k-1)Φ^T(k-1)+Q(k-1)

其中，为当前时刻航迹点的状态预测值，P(k/k-1)为当前时刻航迹点的状态预测协方差矩阵，为前一时刻航迹点的状态预测值，P(k-1/k-1)为前一时刻航迹点的状态预测协方差矩阵；

所述航迹点状态信息包括前一时刻状态转移矩阵Φ(k-1)、前一时刻控制矩阵U(k-1)、前一时刻机动加速度均值前一时刻航迹点的状态预测协方差矩阵P(k/k-1)、前一时刻状态转移矩阵的转置矩阵Φ^T(k-1)、前一时刻噪声方差矩阵Q(k-1)。

3.根据权利要求2所述的基于ADS-B***航迹点的修正方法，其特征在于，所述根据航迹点的状态信息预测当前时刻航迹点的测量预测值以及计算当前时刻航迹点的残差向量协方差矩阵通过下述方程式计算得到：

$\hat{Y}(k/k-1)=H\left(k\right)\hat{X}(k/k-1)$

S(k)＝H(k)P(k/k-1)H^T(k)+R(k)

其中，为当前时刻航迹点的测量预测值，S(k)为当前时刻航迹点的残差向量协方差矩阵；

所述航迹点状态信息还包括当前时刻测量矩阵H(k)、当前时刻测量矩阵的转置矩阵H^T(k)、当前时刻观测噪声R(k)。

4.根据权利要求3所述的基于ADS-B***航迹点的修正方法，其特征在于，所述识别当前时刻航迹点是否为野值，并根据当前时刻航迹点是否为野值确定当前时刻航迹点的衰减因子通过下述方程式计算得到：

$\tilde{Y}\left(k\right)=|Y\left(k\right)-\hat{Y}(k/k-1)|\&le;d\sqrt{S\left(k\right)}$

$r\left(k\right)=\left\{\begin{array}{ccc}1& \left|\tilde{Y}\left(k\right)\right|\&le;& d\sqrt{S\left(k\right)}\\ {\left\{\frac{3\sqrt{S\left(k\right)}}{\left|\tilde{Y}\left(k\right)\right|}\right\}}^{1/2}& \left|\tilde{Y}\left(k\right)\right|>& d\sqrt{S\left(k\right)}\end{array}\right.$

其中,为当前时刻航迹点的新息，当时，当前时刻航迹点不是野值，当时，当前时刻航迹点为野值；Y(k)为当前时刻航迹点的测量值，为当前时刻航迹点的测量预测值，d设置为3，r(k)为当前时刻航迹点的衰减因子。

5.根据权利要求4所述的基于ADS-B***航迹点的修正方法，其特征在于，所述根据当前时刻航迹点的衰减因子计算当前时刻滤波器增益通过下述方程式计算得到：

K(k)＝r(k)P(k)H^T(k)S^-1(k)

其中，K(k)为当前时刻滤波器增益，S^-1(k)为当前时刻航迹点的残差向量协方差矩阵的逆矩阵。

6.根据权利要求5所述的基于ADS-B***航迹点的修正方法，其特征在于，所述根据当前时刻航迹点的测量预测值、当前航迹点的残差向量协方差矩阵以及当前时刻滤波器增益修正当前时刻航迹点的状态预测值以及当前航迹点的状态预测协方差矩阵通过下述方程式计算得到：

$\hat{X}(k/k)=\hat{X}(k/k-1)+K\left(k\right)[Y\left(k\right)-\hat{Y}(k/k-1)]$

P(k/k)＝P(k/k-1)-K(k)S(k)K^T(k)

其中，为修正后的当前时刻航迹点的状态预测值，P(k/k)为修正后的当前航迹点的状态预测协方差矩阵，K^T(k)为当前时刻滤波器增益转置矩阵。

7.根据权利要求1至6任意一项所述的基于ADS-B***航迹点的修正方法，其特征在于，在对航迹点按接收时间进行排列之后，预测当前时刻航迹点的状态预测值以及计算当前时刻航迹点的状态预测协方差矩阵之前还包括步骤：

将每两个航迹点的时间间隔与滤波器的采样时间间隔进行比较；

若时间间隔不一致，则判断该两个航迹点之间存在航迹点丢失；

用野值对丢失的航迹点进行填充。

8.一种基于ADS-B***航迹点的修正***，其特征在于，包括预处理单元、计算修正单元；

所述预处理单元对航迹点按接收时间进行排列；

所述计算修正单元根据航迹点的状态信息预测当前时刻航迹点的状态预测值以及计算当前时刻航迹点的状态预测协方差矩阵；根据航迹点的状态信息预测当前时刻航迹点的测量预测值以及计算当前时刻航迹点的残差向量协方差矩阵；识别当前时刻航迹点是否为野值，并根据当前时刻航迹点是否为野值确定当前时刻航迹点的衰减因子；根据当前时刻航迹点的衰减因子计算当前时刻滤波器增益；根据当前时刻航迹点的测量预测值、当前航迹点的残差向量协方差矩阵以及当前时刻滤波器增益修正当前时刻航迹点的状态预测值以及当前航迹点的状态预测协方差矩阵。

9.根据权利要求8所述的基于ADS-B***航迹点的修正***，其特征在于，所述计算修正单元通过下述方程式进行计算：

$\hat{X}(k/k-1)=\mathrm{\&Phi;}(k/k-1)\hat{X}(k-1/k-1)+U(k-1)\stackrel{\&OverBar;}{a}(k-1)$

P(k/k-1)＝Φ(k-1)P(k-1/k-1)Φ^T(k-1)+Q(k-1)

$\hat{Y}(k/k-1)=H\left(k\right)\hat{X}(k/k-1)$

S(k)＝H(k)P(k/k-1)H^T(k)+R(k)

$\tilde{Y}\left(k\right)=|Y\left(k\right)-\hat{Y}(k/k-1)|\&le;d\sqrt{S\left(k\right)}$

$r\left(k\right)=\left\{\begin{array}{ccc}1& \left|\tilde{Y}\left(k\right)\right|\&le;& d\sqrt{S\left(k\right)}\\ {\left\{\frac{3\sqrt{S\left(k\right)}}{\left|\tilde{Y}\left(k\right)\right|}\right\}}^{1/2}& \left|\tilde{Y}\left(k\right)\right|>& d\sqrt{S\left(k\right)}\end{array}\right.$

K(k)＝r(k)P(k)H^T(k)S^-1(k)

$\hat{X}(k/k)=\hat{X}(k/k-1)+K\left(k\right)[Y\left(k\right)-\hat{Y}(k/k-1)]$

P(k/k)＝P(k/k-1)-K(k)S(k)K^T(k)

其中，为当前时刻航迹点的状态预测值，Φ(k-1)为前一时刻状态转移矩阵，U(k-1)为前一时刻控制矩阵，为前一时刻机动加速度均值，为前一时刻航迹点的状态预测值；P(k/k-1)为当前时刻航迹点的状态预测协方差矩阵，Q(k-1)为噪声方差矩阵，P(k-1/k-1)为前一时刻航迹点的状态预测协方差矩阵，Φ^T(k-1)为前一时刻状态转移矩阵的转置矩阵；为当前时刻航迹点的测量预测值，H(k)为当前时刻测量矩阵；S(k)为当前时刻航迹点的残差向量协方差矩阵，R(k)当前时刻观测噪声，H^T(k)为当前时刻测量矩阵的转置矩阵；为当前时刻航迹点的新息，Y(k)为当前时刻航迹点的测量值，为当前时刻航迹点的测量预测值，d设置为3，r(k)为当前时刻航迹点的衰减因子；K(k)为当前时刻滤波器增益，S^-1(k)为当前时刻航迹点的残差向量协方差矩阵的逆矩阵；为修正后的当前时刻航迹点的状态预测值，P(k/k)为修正后的当前航迹点的状态预测协方差矩阵，K^T(k)为当前时刻滤波器增益转置矩阵。

10.根据权利要求8所述的基于ADS-B***航迹点的修正***，其特征在于，还包括填充单元；

所述填充单元将每两个航迹点的时间间隔与滤波器的采样时间间隔进行比较；若时间间隔不一致，则判断该两个航迹点之间存在航迹点丢失；用野值对丢失的航迹点进行填充。