ADS-B数据的滤波方法、装置、计算机设备及存储介质

技术领域

本申请涉及航空技术领域，尤其涉及一种ADS-B数据的滤波方法、装置、计算机设备及存储介质。

背景技术

ADS-B即广播式自动相关监视系统，是一个集通信与监视于一体的信息系统，支持数据共享和空中监视等优点，是下一代航空运输系统的关键组成部分。机载飞机上的ADS-B系统可以定期播放基于卫星的飞机位置和飞行识别信息，ADS-B接收机与空管系统、其它飞机的机载ADS-B结合起来提供精确、实时的冲突信息。为新航行系统增强和扩展了非常丰富的功能，同时也带来了潜在的经济效益和社会效益。

但由于ADS-B数据通过1090兆赫的数字信息链以每半秒一次的频率发射，同雷达类似，它受“视线”的限制。地面站接收信号的性能取决于高度、离地面站的距离和阻隔信号的地形情况，其监视范围取决于地面站的数量以及分布情况。由于以上提到的接收限制问题，在实际运行环境中，接收到的ADS-B数据会遇到数据异常，比如位置的漂移、时间戳的误差，数据缺失和数据中断的诸多问题。

目前，对ADS-B数据处理方法主要都是直接依据飞机飞行的数据标准，对超过飞机飞行标准的数据直接进行过滤，比如民航飞机的飞行速度不会超过1000千米/小时，对于数据中超出这个范围的异常数据直接进行过滤，这种基于飞机飞行标准的数据清理过滤方法虽然简单，但是会过滤掉处于异常飞行状态下的飞行数据，这对飞机飞行中异常状态的识别和捕捉是不利的。

发明内容

本申请实施例提供一种ADS-B数据的滤波方法、装置、计算机设备及存储介质，用于对ADS-B数据进行滤波，去除ADS-B数据中的噪声。

本发明实施例提供一种ADS-B数据的滤波方法，所述方法包括：

获取当前时刻的ADS-B数据的观测高斯分布和ADS-B数据的预测高斯分布，以及过程误差矩阵Q，所述观测高斯分布包括协方差矩阵R，所述预测高斯分布包括协方差矩阵P；

将所述预测高斯分布中的协方差矩阵P加上所述过程误差矩阵Q；

对所述观测高斯分布和加上所述过程误差矩阵Q的所述预测高斯分布进行融合，得到所述当前时刻的校准高斯分布；

根据所述当前时刻的校准高斯分布更新所述当前时刻的ADS-B数据的观测高斯分布；

其中，观测的经纬度坐标符合所述观测高斯分布，预测的状态向量符合所述预测高斯分布。

本发明实施例提供一种ADS-B数据的滤波装置，所述装置包括：

获取模块，用于获取当前时刻的ADS-B数据的观测高斯分布和ADS-B数据的预测高斯分布，以及过程误差矩阵Q，所述观测高斯分布包括协方差矩阵R，所述预测高斯分布包括协方差矩阵P；其中，观测的经纬度坐标符合所述观测高斯分布，预测的状态向量符合所述预测高斯分布；

计算模块，用于将所述预测高斯分布中的协方差矩阵P加上所述过程误差矩阵Q；

融合模块，用于对所述观测高斯分布和加上所述过程误差矩阵Q的所述预测高斯分布进行融合，得到所述当前时刻的校准高斯分布；

更新模块，用于根据所述当前时刻的校准高斯分布更新所述当前时刻的ADS-B数据的观测高斯分布。

一种计算机设备，包括存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现上述ADS-B数据的滤波方法。

一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序被处理器执行时实现上述ADS-B数据的滤波方法。

本发明提供一种ADS-B数据的滤波方法、装置、计算机设备及存储介质，首先获取当前时刻的ADS-B数据的观测高斯分布和ADS-B数据的预测高斯分布，以及过程误差矩阵Q，观测高斯分布包括协方差矩阵R，预测高斯分布包括协方差矩阵P；然后将预测高斯分布中的协方差矩阵P加上过程误差矩阵Q；对观测高斯分布和加上过程误差矩阵Q的预测高斯分布进行融合，得到当前时刻的校准高斯分布，跟据当前时刻的校准高斯分布更新当前时刻的ADS-B数据的观测高斯分布。从而通过本发明可实现对ADS-B中的经纬度和速度数据进行滤波，去除数据中的噪声，修正位置轨迹和速度数据中的异常点。

附图说明

图1为本申请提供的ADS-B数据的滤波方法流程图；

图2为本申请提供的观测误差

图3为本申请提供的连续ADS-B点的观测误差示意图；

图4为本申请提供的高斯分布融合示意图；

图5为本申请提供的ADS-B数据的滤波装置的结构示意图。

具体实施方式

为了更好的理解上述技术方案，下面通过附图以及具体实施例对本申请实施例的技术方案做详细的说明，应当理解本申请实施例以及实施例中的具体特征是对本申请实施例技术方案的详细的说明，而不是对本申请技术方案的限定，在不冲突的情况下，本申请实施例以及实施例中的技术特征可以相互组合。

需要说明的是，本发明应用于三种不同的滤波器：线性滤波器、非线性滤波器和交互式多模型滤波器（Interacting Multiple Models，IMM），可根据精度需求和计算资源选择使用。这三种滤波器中，线性滤波器计算速度最快，但精度最低；交互式多模型滤波器精度最高，但计算速度最慢；非线性滤波器介于两者之间。

请参阅图1，为本发明当中的ADS-B数据的滤波方法，该方法用于对ADS-B中的经纬度和速度数据进行滤波，去除数据中的噪声，修正位置轨迹和速度数据中的异常点。该方法具体包括步骤S10-步骤S40：

步骤S10，获取当前时刻的ADS-B数据的观测高斯分布和ADS-B数据的预测高斯分布，以及过程误差矩阵Q。

其中，ADS-B数据用于表示飞机的飞行状态，每条ADS-B数据至少对应有获取时间、飞行高度和飞行速度，具体可以包括：航班号、计划日期、出发机场、到达机场、飞机注册号、飞行经度、飞行纬度、飞行速度、飞行高度。具体的，本发明实施例可与ADS-B数据数据库相连，用于实时获取目标飞机在高空中广播式自动相关监视（Automatic dependentsurveillance–broadcast，ADS-B）数据。广播式自动相关监视，指飞机启动后，无需人工操作或者询问，可以自动地按照一定频率从相关机载设备获取参数向其他飞机或地面站广播飞机的各种状态信息，以供管制员对飞机状态进行监控。

在本发明实施例中，观测高斯分布包括协方差矩阵R，预测高斯分布包括协方差矩阵P，观测的经纬度坐标符合观测高斯分布，预测的状态向量符合预测高斯分布。

其中，x和y表示飞机在平面坐标系中的位置，粗体的

线性滤波器：独立地处理飞机在x和y方向上的运动，预测的状态向量

非线性滤波器：采用了恒定转率和加速度模型（Constant Turn Rate andAcceleration，CTRA）。模型的状态向量里包括：飞机的位置、速度

交互式多模型滤波器：融合了3种非线性模型，在CTRA模型的基础上，分别令加速度a为零、角速度

具体的，获取当前时刻的ADS-B数据的预测高斯分布，包括：将上一时刻的校准状态向量输入到状态转移公式得到当前时刻的ADS-B数据的观测高斯分布。其中，上一时刻和当前时刻的时间间隔为

而初始状态

若采用线性滤波器，状态转移公式为：

其中，

若采用非线性滤波器，CTRA模型的状态转移方程如下：

角速度

上述方程组是非线性的，需要用特定的线性化技术来处理，本发明实施例通过无损变换的技术（Unscented Transform，UT）。具体过程为：在t时刻，状态向量

交互式多模型滤波器把三个不同的非线性滤波器组合起来，这三个滤波器的状态转移方程组都是CTRA模型方程组的变体。CTRV模型中，始终有加速度

进一步地，在交互式多模型滤波器中，组合了三种不同的非线性滤波器，分别对应飞机的三种运动状态。CTRV（恒定转率和速度模型）对应飞机转弯的情形，CCA（常曲率和加速度模型）对应飞机加速减速的情形，CSAV（恒定转向角和速度模型）对应飞机匀速直线运动的情形。

在本实施例中，每个ADS-B点的经纬度被投影到平面坐标系，分别用x和y表示。按照卡尔曼滤波的假设，观测向量

观测误差

但由于滤波器在计算时需要用到固定的时间戳，因此把ADS-B中的时间戳作为固定时间戳，再把时间戳的误差转化为位置的误差。令

观测误差

因此，在本发明提供的一个实施例中，获取当前时刻的ADS-B数据的观测高斯分布，包括：将当前时刻的观测的经纬度坐标作为观测高斯分布中的高斯分布均值；通过公式计算观测高斯分布中的协方差矩阵R：

其中，

一种估计

进一步地，一个航班的ADS-B可能来自于多个不同的数据源（例如空管局数据源、航科院数据源等），每个数据源的

步骤S20，将预测高斯分布中的协方差矩阵P加上过程误差矩阵Q。

在本发明实施例中，在计算每一个时刻的预测高斯分布时，还需要把未考虑到的不确定度（比如飞机改变了机动状态）加到协方差矩阵

线性滤波器：假设有未追踪到的加速度

非线性滤波器：计算采用局部线性化的方法，即扩展卡尔曼滤波（ExpandedKalman Filter，EKF）方法计算过程误差矩阵Q。假设有未追踪到的加速度a和角速度

其中，

在卡尔曼滤波中，这个情景涉及到了

步骤S30，对观测高斯分布和加上过程误差矩阵Q的预测高斯分布进行融合，得到当前时刻的校准高斯分布。

如图4所示，位置的预测分布用圆圈A表示，位置的均值为圆圈A的圆心（从预测的状态向量

例如，飞机在跑道上，需要记录飞机每个时刻的位置。一开始飞机在0m处，预测飞机的速度是10m/s，于是1s的时候，预测飞机在10m处（预测飞机在10m处，速度10m/s）。而观测的飞机在20m处，将观测的位置和预测的位置进行融合得到飞机在15m处。也就是说，预测的飞机速度慢了，于是将飞机的速度提高到15m/s（飞机在15m处，速度15m/s）。2s的时候，预测飞机在30m处，而观测的飞机在35m处，接下来又是结合预测和观测的飞机位置修正飞机预测的位置和速度。

步骤S40，根据当前时刻的校准高斯分布更新当前时刻的ADS-B数据的观测高斯分布。

在本发明实施例中，根据当前时刻的校准高斯分布更新当前时刻的观测高斯分布。预测完

进一步地，在每次预测时引入自适应机制。当下一个点的时间间隔大于一个阈值时，预测分布的协方差矩阵

进一步地，在每次预测时引入异常点判断机制。当观测的位置与预测位置的马氏距离大于一个阈值时，将该观测标记为异常点，跳过步骤S30，用预测的位置直接替换原始数据（观测的数据）中的位置。

在一个可选实施例中，在获取当前时刻的ADS-B数据的观测高斯分布和ADS-B数据的预测高斯分布之前，还可以对ADS-B数据进行初步的修正，具体包括如下步骤：

1.数据读取

从Kafka的队列中接收实时信息，获取正在飞行的航班，进而在Redis中读取该航班所有的ADS-B数据。

2.数据规整

将ADS-B数据中的经纬度进行投影变换，同时做初步的数据筛选。具体地，使用Lambert投影将球面坐标系下的经纬度数据投影到平面坐标系（x，y），单位为米。投影时，取该航班所有ADS-B经度的中位数作为投影时经度的基准点，纬度的25%和75%分位数作为纬度的基准点。同时，如果该航班开头和结尾的几个点速度、高度均为0，则将这些点去除。

3.初始指标计算

计算原始ADS-B数据中的尖刺数、速度越界和加速度越界的点的数量。其中，尖刺是指三个ADS-B点的位置形成了小于20度的夹角；速度越界是指原始ADS-B中速度这一项的值，或者相邻两点根据位置距离和时间差计算的速度超过340m/s；加速度越界是指用相邻两点的速度差和时间差计算加速度，或是用连续三点的距离和时间差计算加速度，大于5m/s

4.位置数据插值

对数据进行初步修正。对于上一步发现的位置异常的点，先去掉x和y的值，再利用前后两点进行线性插值。插值时利用ADS-B中前后两点的位置和三个点的时间戳，计算中间点的位置，替换掉原先x和y的值。

在另一个可选实施例中，对观测高斯分布和加上过程误差矩阵Q的预测高斯分布进行融合得到当前时刻的校准高斯分布之后，还可以包括：使用RTS平滑算法（Rauch–Tung–Striebel smoother）对滤波后的数据进行平滑。滤波是按时间从前往后对数据进行处理，而平滑是从后往前对数据进行处理。平滑后，数据的精度更高，可以替换每个ADS-B点位置的值。

滤波与平滑之后，还需计算以下指标来判断滤波是否发散，即滤波器未能正确追踪飞机的运动状态，偏离了正常航线。一是平均修正距离，即滤波与平滑后，计算每个ADS-B点与原先位置的距离，如果平均值大于某个阈值，认为滤波发散，拒绝使用滤波后的数据。二是计算滤波与平滑后的点距离原先ADS-B轨迹的垂直距离，如果存在偏移较大的点，则拒绝使用滤波平滑后的数据。三是计算轨迹中的尖刺数，如果尖刺数升高，则拒绝使用滤波平滑后的数据。计算指标之后，如果接受滤波与平滑的结果，则将平滑后每一步状态向量

应理解，上述实施例中各步骤的序号的大小并不意味着执行顺序的先后，各过程的执行顺序应以其功能和内在逻辑确定，而不应对本发明实施例的实施过程构成任何限定。

在一实施例中，提供一种ADS-B数据的滤波装置，该ADS-B数据的滤波装置与上述实施例中ADS-B数据的滤波方法一一对应。如图5所示，所述ADS-B数据的滤波装置各功能模块详细说明如下：

获取模块10，用于获取当前时刻的ADS-B数据的观测高斯分布和ADS-B数据的预测高斯分布，以及过程误差矩阵Q，所述观测高斯分布包括协方差矩阵R，所述预测高斯分布包括协方差矩阵P；其中，观测的经纬度坐标符合所述观测高斯分布，预测的状态向量符合所述预测高斯分布。

计算模块20，用于将所述预测高斯分布中的协方差矩阵P加上所述过程误差矩阵Q。

融合模块30，用于对所述观测高斯分布和加上所述过程误差矩阵Q的所述预测高斯分布进行融合，得到所述当前时刻的校准高斯分布。

更新模块40，用于根据所述当前时刻的校准高斯分布更新所述当前时刻的ADS-B数据的观测高斯分布。

所述获取模块10，具体用于将上一时刻的校准状态向量输入到状态转移公式得到所述当前时刻的ADS-B数据的观测高斯分布，所述上一时刻和所述当前时刻的时间间隔为

所述计算模块20，还用于将飞机的第一个高度大于0的点作为起飞时的点，根据起飞时的点的飞机位置坐标和飞机的历史数据计算初始时刻的ADS-B数据的预测高斯分布。

进一步的，所述获取模块10，具体用于：

将所述当前时刻的观测的经纬度坐标作为所述观测高斯分布中的高斯分布均值；

通过下述公式计算所述观测高斯分布中的协方差矩阵R：

其中，

在本发明提供的一个实施例中，所述获取模块10，还用于获取历史ADS-B数据，从所述历史ADS-B数据中获取直线飞行数据；

还用于通过对所述直线飞行数据进行直线拟合轨迹，旋转到与x轴平行，得到x方向的方差

所述计算模块20，还用于计算所述x方向的方差

在本发明提供的一个实施例中，若采用线性滤波器，所述状态向量包含飞机的位置、速度

进一步的，所述计算模块20，具体用于：

若采用线性滤波器，状态转移公式为：

其中，

关于ADS-B数据的滤波装置的具体限定可以参见上文中对于ADS-B数据的滤波方法的限定，在此不再赘述。上述设备中的各个模块可全部或部分通过软件、硬件及其组合来实现。上述各模块可以硬件形式内嵌于或独立于计算机设备中的处理器中，也可以以软件形式存储于计算机设备中的存储器中，以便于处理器调用执行以上各个模块对应的操作。

在一个实施例中，提供了一种计算机设备，该计算机设备可以是服务器。该计算机设备包括通过系统总线连接的处理器、存储器、网络接口和数据库。其中，该计算机设备的处理器用于提供计算和控制能力。该计算机设备的存储器包括非易失性存储介质、内存储器。该非易失性存储介质存储有操作系统、计算机程序和数据库。该内存储器为非易失性存储介质中的操作系统和计算机程序的运行提供环境。该计算机设备的网络接口用于与外部的终端通过网络连接通信。该计算机程序被处理器执行时以实现一种ADS-B数据的滤波方法。

在一个实施例中，提供了一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，处理器执行计算机程序时实现以下步骤：

获取当前时刻的ADS-B数据的观测高斯分布和ADS-B数据的预测高斯分布，以及过程误差矩阵Q，所述观测高斯分布包括协方差矩阵R，所述预测高斯分布包括协方差矩阵P；

将所述预测高斯分布中的协方差矩阵P加上所述过程误差矩阵Q；

对所述观测高斯分布和加上所述过程误差矩阵Q的所述预测高斯分布进行融合，得到所述当前时刻的校准高斯分布；

根据所述当前时刻的校准高斯分布更新所述当前时刻的ADS-B数据的观测高斯分布；

其中，观测的经纬度坐标符合所述观测高斯分布，预测的状态向量符合所述预测高斯分布。

在一个实施例中，提供了一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现以下步骤：

获取当前时刻的ADS-B数据的观测高斯分布和ADS-B数据的预测高斯分布，以及过程误差矩阵Q，所述观测高斯分布包括协方差矩阵R，所述预测高斯分布包括协方差矩阵P；

将所述预测高斯分布中的协方差矩阵P加上所述过程误差矩阵Q；

对所述观测高斯分布和加上所述过程误差矩阵Q的所述预测高斯分布进行融合，得到所述当前时刻的校准高斯分布；

根据所述当前时刻的校准高斯分布更新所述当前时刻的ADS-B数据的观测高斯分布；

其中，观测的经纬度坐标符合所述观测高斯分布，预测的状态向量符合所述预测高斯分布。

本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的计算机程序可存储于一非易失性计算机可读取存储介质中，该计算机程序在执行时，可包括如上述各方法的实施例的流程。其中，本申请所提供的各实施例中所使用的对存储器、存储、数据库或其它介质的任何引用，均可包括非易失性和/或易失性存储器。非易失性存储器可包括只读存储器（ROM）、可编程ROM（PROM）、电可编程ROM（EPROM）、电可擦除可编程ROM（EEPROM）或闪存。易失性存储器可包括随机存取存储器（RAM）或者外部高速缓冲存储器。作为说明而非局限，RAM以多种形式可得，诸如静态RAM（SRAM）、动态RAM（DRAM）、同步DRAM（SDRAM）、双数据率SDRAM（DDRSDRAM）、增强型SDRAM（ESDRAM）、同步链路（Synchlink） DRAM（SLDRAM）、存储器总线（Rambus）直接RAM（RDRAM）、直接存储器总线动态RAM（DRDRAM）、以及存储器总线动态RAM（RDRAM）等。

所属领域的技术人员可以清楚地了解到，为了描述的方便和简洁，仅以上述各功能单元、模块的划分进行举例说明，实际应用中，可以根据需要而将上述功能分配由不同的功能单元、模块完成，即将所述装置的内部结构划分成不同的功能单元或模块，以完成以上描述的全部或者部分功能。

以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围，均应包含在本发明的保护范围之内。