CN113257044B - Ads-b数据的滤波方法、装置、计算机设备及存储介质 - Google Patents

Ads-b数据的滤波方法、装置、计算机设备及存储介质 Download PDF

Info

Publication number
CN113257044B
CN113257044B CN202110775367.0A CN202110775367A CN113257044B CN 113257044 B CN113257044 B CN 113257044B CN 202110775367 A CN202110775367 A CN 202110775367A CN 113257044 B CN113257044 B CN 113257044B
Authority
CN
China
Prior art keywords
gaussian distribution
ads
data
observed
predicted
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110775367.0A
Other languages
English (en)
Other versions
CN113257044A (zh
Inventor
卞磊
薄满辉
唐红武
刘宇
姚远
张聪
王殿胜
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China Travelsky Mobile Technology Co Ltd
Original Assignee
China Travelsky Mobile Technology Co Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by China Travelsky Mobile Technology Co Ltd filed Critical China Travelsky Mobile Technology Co Ltd
Priority to CN202110775367.0A priority Critical patent/CN113257044B/zh
Publication of CN113257044A publication Critical patent/CN113257044A/zh
Application granted granted Critical
Publication of CN113257044B publication Critical patent/CN113257044B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G08—SIGNALLING
    • G08G—TRAFFIC CONTROL SYSTEMS
    • G08G5/00—Traffic control systems for aircraft, e.g. air-traffic control [ATC]
    • G08G5/0073—Surveillance aids
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Physics (AREA)
  • Mathematical Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computational Mathematics (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Computing Systems (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Biology (AREA)
  • Operations Research (AREA)
  • Probability & Statistics with Applications (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本申请提供了一种ADS‑B数据的滤波方法、装置、计算机设备及存储介质，涉及航空技术领域，用于对ADS‑B数据进行滤波，去除ADS‑B数据中的噪声。方法主要包括：获取当前时刻的ADS‑B数据的观测高斯分布和ADS‑B数据的预测高斯分布，以及过程误差矩阵Q，所述观测高斯分布包括协方差矩阵R，所述预测高斯分布包括协方差矩阵P；将所述预测高斯分布中的协方差矩阵P加上所述过程误差矩阵Q；对所述观测高斯分布和加上所述过程误差矩阵Q的所述预测高斯分布进行融合，得到所述当前时刻的校准高斯分布；根据所述当前时刻的校准高斯分布更新所述当前时刻的ADS‑B数据的观测高斯分布；其中，观测的经纬度坐标符合所述观测高斯分布，预测的状态向量符合所述预测高斯分布。

Description

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为本申请提供的观测误差
Figure DEST_PATH_353359DEST_PATH_IMAGE001
的构建示意图；
图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表示飞机在平面坐标系中的位置，粗体的
Figure DEST_PATH_DEST_PATH_IMAGE002
表示预测的状态向量。根据卡尔曼滤波的假设，预测的状态向量
Figure DEST_PATH_113243DEST_PATH_IMAGE002
服从预测高斯分布
Figure DEST_PATH_565084DEST_PATH_IMAGE003
，
Figure DEST_PATH_DEST_PATH_IMAGE004
为高斯分布的均值，P为高斯分布协方差矩阵。预测的状态向量
Figure DEST_PATH_163555DEST_PATH_IMAGE002
表示对飞机运动状态的追踪，协方差矩阵
Figure DEST_PATH_420005DEST_PATH_IMAGE005
表示状态的不确定度。每次预测飞机的运动状态时，都需要计算预测高斯分布的均值与协方差矩阵。在本发明实施例中，采用不同的滤波器，其对应的状态向量
Figure DEST_PATH_914572DEST_PATH_IMAGE002
的描述如下所示：
线性滤波器：独立地处理飞机在x和y方向上的运动，预测的状态向量
Figure DEST_PATH_978343DEST_PATH_IMAGE002
包含飞机的位置、速度和加速度。因此计算初始状态时，需要得到飞机的位置，以及速度
Figure DEST_PATH_DEST_PATH_IMAGE006
，加速度a分别在x、y方向的分量。
Figure DEST_PATH_DEST_PATH_IMAGE007
非线性滤波器：采用了恒定转率和加速度模型（Constant Turn Rate andAcceleration，CTRA）。模型的状态向量里包括：飞机的位置、速度
Figure DEST_PATH_692090DEST_PATH_IMAGE006
、加速度a、偏航角
Figure DEST_PATH_DEST_PATH_IMAGE008
，角速度
Figure DEST_PATH_829810DEST_PATH_IMAGE009
。此处偏航角
Figure DEST_PATH_495278DEST_PATH_IMAGE008
指的是飞机机头方向与正北方向的夹角，在x、y空间中就是与y轴的夹角，顺时针方向为正，逆时针方向为负，取值范围
Figure DEST_PATH_DEST_PATH_IMAGE010
。
Figure DEST_PATH_452869DEST_PATH_IMAGE011
交互式多模型滤波器：融合了3种非线性模型，在CTRA模型的基础上，分别令加速度a为零、角速度
Figure DEST_PATH_360520DEST_PATH_IMAGE009
为零、加速度a和角速度
Figure DEST_PATH_556009DEST_PATH_IMAGE009
都为零，便得到了恒定转率和速度模型（Constant Turn Rate and Velocity，CTRV）、常曲率和加速度模型（Constant Curvatureand Acceleration，CCA）、恒定转向角和速度模型（Constant Steering Angle andVelocity，CSAV）。
具体的，获取当前时刻的ADS-B数据的预测高斯分布，包括：将上一时刻的校准状态向量输入到状态转移公式得到当前时刻的ADS-B数据的观测高斯分布。其中，上一时刻和当前时刻的时间间隔为
Figure DEST_PATH_DEST_PATH_IMAGE012
。需要说明的是，上一时刻的校准状态向量是根据上一时刻的预测高斯分布和上一时刻的观测高斯分布融合得到的。
而初始状态
Figure DEST_PATH_454695DEST_PATH_IMAGE002
计算方法为：在ADS-B中选取第一个高度大于0的点作为起飞时的点，构造状态向量
Figure DEST_PATH_368424DEST_PATH_IMAGE002
的初始分布。首先计算初始分布的均值
Figure DEST_PATH_814187DEST_PATH_IMAGE013
，初始的x和y选择该点观测位置，速度v、加速度a使用历史数据统计的平均起飞速度、加速度。初始角速度
Figure DEST_PATH_254396DEST_PATH_IMAGE009
为0。初始偏航角
Figure DEST_PATH_261666DEST_PATH_IMAGE008
可以根据起飞时的点和下一个点来确定，进而计算速度v和加速度a在x和y方向的分量（即
Figure DEST_PATH_662691DEST_PATH_IMAGE015
和
Figure DEST_PATH_DEST_PATH_IMAGE016
，加速度的分量计算类似）。然后计算初始分布的协方差矩阵，在一个可选的方式中，初始分布的协方差矩阵
Figure DEST_PATH_210347DEST_PATH_IMAGE017
，
Figure DEST_PATH_138684DEST_PATH_IMAGE019
等；在另一个可选方式中，初始分布的协方差矩阵也可以是设置的一个经验值，本发明实施例不做具体限定。其中，Q是过程误差矩阵。
若采用线性滤波器，状态转移公式为：
Figure DEST_PATH_RE-DEST_PATH_IMAGE020
Figure DEST_PATH_DEST_PATH_IMAGE021
Figure DEST_PATH_DEST_PATH_IMAGE022
其中，
Figure DEST_PATH_457801DEST_PATH_IMAGE012
是两个ADS-B点时间戳的差值，单位为秒。
Figure DEST_PATH_DEST_PATH_IMAGE023
和
Figure DEST_PATH_DEST_PATH_IMAGE024
是t时刻高斯分布的均值和方差，
Figure DEST_PATH_188865DEST_PATH_IMAGE025
和
Figure DEST_PATH_DEST_PATH_IMAGE026
是
Figure DEST_PATH_710851DEST_PATH_IMAGE027
时刻分布的均值和方差，F为计算
Figure DEST_PATH_735439DEST_PATH_IMAGE025
和
Figure DEST_PATH_474725DEST_PATH_IMAGE026
所需要的中间变量。
若采用非线性滤波器，CTRA模型的状态转移方程如下：
Figure DEST_PATH_DEST_PATH_IMAGE028
Figure DEST_PATH_256867DEST_PATH_IMAGE029
Figure DEST_PATH_DEST_PATH_IMAGE030
Figure DEST_PATH_785806DEST_PATH_IMAGE031
Figure DEST_PATH_DEST_PATH_IMAGE032
Figure DEST_PATH_992797DEST_PATH_IMAGE033
角速度
Figure DEST_PATH_778350DEST_PATH_IMAGE035
时，计算x和y的方法有所不同：
Figure DEST_PATH_DEST_PATH_IMAGE036
Figure DEST_PATH_552183DEST_PATH_IMAGE037
上述方程组是非线性的，需要用特定的线性化技术来处理，本发明实施例通过无损变换的技术（Unscented Transform，UT）。具体过程为：在t时刻，状态向量
Figure DEST_PATH_573229DEST_PATH_IMAGE002
的分布是高斯分布
Figure DEST_PATH_DEST_PATH_IMAGE038
，在这个分布中以一定方式得到几个状态向量
Figure DEST_PATH_713354DEST_PATH_IMAGE002
的采样点，将这些采样点代入状态转移方程中，分别计算出
Figure DEST_PATH_433923DEST_PATH_IMAGE027
时刻对应的几个状态向量
Figure DEST_PATH_643188DEST_PATH_IMAGE039
，利用这几个点估计
Figure DEST_PATH_608870DEST_PATH_IMAGE027
时刻状态向量
Figure DEST_PATH_728136DEST_PATH_IMAGE002
的高斯分布
Figure DEST_PATH_DEST_PATH_IMAGE040
。本发明实施例采样方式选择了Van der Merwe采样法。
交互式多模型滤波器把三个不同的非线性滤波器组合起来，这三个滤波器的状态转移方程组都是CTRA模型方程组的变体。CTRV模型中，始终有加速度
Figure DEST_PATH_760551DEST_PATH_IMAGE041
；CCA模型中，始终有
Figure DEST_PATH_DEST_PATH_IMAGE042
=0；CSAV模型中，以上两个条件同时满足。计算状态转移时，每个滤波器分别用无损变换计算
Figure DEST_PATH_129216DEST_PATH_IMAGE027
时刻状态向量
Figure DEST_PATH_367430DEST_PATH_IMAGE002
的分布，再根据交互式多模型的方法加权平均。
进一步地，在交互式多模型滤波器中，组合了三种不同的非线性滤波器，分别对应飞机的三种运动状态。CTRV（恒定转率和速度模型）对应飞机转弯的情形，CCA（常曲率和加速度模型）对应飞机加速减速的情形，CSAV（恒定转向角和速度模型）对应飞机匀速直线运动的情形。
在本实施例中，每个ADS-B点的经纬度被投影到平面坐标系，分别用x和y表示。按照卡尔曼滤波的假设，观测向量
Figure DEST_PATH_DEST_PATH_IMAGE043
服从二维的高斯分布。实际计算时，每个点依次得到一个高斯分布。对于其中某一个点，取该点观测的数据（即该点经纬度投影后的x和y）作为高斯分布的均值，再计算一个观测误差R作为高斯分布的协方差矩阵R。
观测误差
Figure DEST_PATH_777421DEST_PATH_IMAGE001
是对观测数据不确定度的建模，
Figure DEST_PATH_810099DEST_PATH_IMAGE001
数值越大，高斯分布的峰越宽，说明观测数据的不确定度越大，即真实位置越有可能离观测数据越远。本发明实施例包含两部分观测误差，第一部分是位置x和y的测量误差，这个误差来自于观测飞机位置的设备；第二部分是时间戳t的误差，即认为ADS-B数据中的时间戳也是不准确的。
但由于滤波器在计算时需要用到固定的时间戳，因此把ADS-B中的时间戳作为固定时间戳，再把时间戳的误差转化为位置的误差。令
Figure DEST_PATH_DEST_PATH_IMAGE044
表示位置误差，即观测位置与飞机真实位置的距离，
Figure DEST_PATH_DEST_PATH_IMAGE045
表示时间戳误差，即ADS-B的时间戳与真实时间戳的差别，根据
Figure DEST_PATH_603742DEST_PATH_IMAGE047
，转化的位置误差与飞机的速度大小成正比（此处
Figure DEST_PATH_427340DEST_PATH_IMAGE006
可以用ADS-B中的速度）。这个误差的方向是沿着速度方向的，因此需要引入飞行的偏航角
Figure DEST_PATH_255619DEST_PATH_IMAGE008
，进行旋转变换。此处偏航角
Figure DEST_PATH_724778DEST_PATH_IMAGE008
是飞机飞行方向与正北方向的角度，顺时针为正，可根据当前点和上一个点的位置来近似计算。
观测误差
Figure DEST_PATH_271297DEST_PATH_IMAGE001
的构建方法是先在沿速度方向和垂直速度方向的空间里构建协方差矩阵
Figure DEST_PATH_DEST_PATH_IMAGE048
，再通过旋转变换得到平面坐标系下的协方差矩阵
Figure DEST_PATH_553111DEST_PATH_IMAGE001
。图2说明了观测误差
Figure DEST_PATH_360530DEST_PATH_IMAGE001
的构建过程。观测的位置作为高斯分布的均值，在图中是坐标的中心点。代表不确定度大小的方差
Figure DEST_PATH_590DEST_PATH_IMAGE049
可以直接相加。叠加了沿速度方向的不确定度后，高斯分布的等概率线从圆形变成椭圆（图中等概率线与均值的距离为一倍标准差），随后椭圆旋转了一个角度
Figure DEST_PATH_768826DEST_PATH_IMAGE008
，长轴是沿着速度方向的。图3是几个连续ADS-B点的观测误差示意。图中横轴为x，纵轴为y，单位是米。
因此，在本发明提供的一个实施例中，获取当前时刻的ADS-B数据的观测高斯分布，包括：将当前时刻的观测的经纬度坐标作为观测高斯分布中的高斯分布均值；通过公式计算观测高斯分布中的协方差矩阵R：
Figure DEST_PATH_DEST_PATH_IMAGE050
Figure DEST_PATH_447807DEST_PATH_IMAGE051
Figure DEST_PATH_DEST_PATH_IMAGE052
Figure DEST_PATH_391623DEST_PATH_IMAGE053
其中，
Figure DEST_PATH_DEST_PATH_IMAGE054
为高斯分布垂直速度方向的标准差，代表观测误差中位置x和y的测量误差，
Figure DEST_PATH_999322DEST_PATH_IMAGE055
为上一时刻t的测量误差，即时间戳t的测量误差，乘以速度
Figure DEST_PATH_753389DEST_PATH_IMAGE006
后转化为位置的误差
Figure DEST_PATH_DEST_PATH_IMAGE056
，
Figure DEST_PATH_471946DEST_PATH_IMAGE056
为高斯分布平行速度方向的标准差；
Figure DEST_PATH_129324DEST_PATH_IMAGE006
为飞机速度，
Figure DEST_PATH_111186DEST_PATH_IMAGE008
为飞行偏航角。高斯分布的方差直接相加，代表不确定度的叠加，因此在沿着速度方向，高斯分布的方差变为
Figure DEST_PATH_713069DEST_PATH_IMAGE057
。A是旋转矩阵，
Figure DEST_PATH_208553DEST_PATH_IMAGE008
是偏航角，也就是飞行方向，正北方向为0度，顺时针方向为正。
Figure DEST_PATH_720437DEST_PATH_IMAGE048
和
Figure DEST_PATH_DEST_PATH_IMAGE058
为计算协方差矩阵R所以需要的中间变量。
一种估计
Figure DEST_PATH_545305DEST_PATH_IMAGE054
和
Figure DEST_PATH_DEST_PATH_IMAGE059
的方法是：取多段历史ADS-B数据，抽出其中直线飞行的片段，用直线拟合轨迹，旋转到与x轴平行，分别得到x和y方向的方差，也就是
Figure DEST_PATH_805122DEST_PATH_IMAGE057
和
Figure DEST_PATH_DEST_PATH_IMAGE060
，二者相减得到
Figure DEST_PATH_DEST_PATH_IMAGE061
。
Figure DEST_PATH_944111DEST_PATH_IMAGE056
除以ADS-B中的速度得到
Figure DEST_PATH_435135DEST_PATH_IMAGE059
。多次计算取平均值。
进一步地，一个航班的ADS-B可能来自于多个不同的数据源（例如空管局数据源、航科院数据源等），每个数据源的
Figure DEST_PATH_257335DEST_PATH_IMAGE054
和
Figure DEST_PATH_974756DEST_PATH_IMAGE059
大小不一样。根据上面的分析，
Figure DEST_PATH_42069DEST_PATH_IMAGE054
和
Figure DEST_PATH_387599DEST_PATH_IMAGE059
越大，观测误差
Figure DEST_PATH_147745DEST_PATH_IMAGE001
越大，认为这个源的数据越不可信。因此计算每个点的
Figure DEST_PATH_585417DEST_PATH_IMAGE001
时，需要根据其所属的数据源选择不同的
Figure DEST_PATH_456421DEST_PATH_IMAGE054
和
Figure DEST_PATH_656459DEST_PATH_IMAGE059
值。
步骤S20，将预测高斯分布中的协方差矩阵P加上过程误差矩阵Q。
在本发明实施例中，在计算每一个时刻的预测高斯分布时，还需要把未考虑到的不确定度（比如飞机改变了机动状态）加到协方差矩阵
Figure DEST_PATH_853085DEST_PATH_IMAGE005
中，才得到最终的预测高斯分布。这个额外的不确定度记作过程误差矩阵Q。
Figure DEST_PATH_DEST_PATH_IMAGE062
是与
Figure DEST_PATH_444298DEST_PATH_IMAGE005
维度相同的矩阵。
线性滤波器：假设有未追踪到的加速度
Figure DEST_PATH_DEST_PATH_IMAGE063
和
Figure DEST_PATH_DEST_PATH_IMAGE064
，过程误差矩阵Q使用离散高斯白噪声，矩阵分为两块，每块3阶。计算
Figure DEST_PATH_56676DEST_PATH_IMAGE062
时需要输入时间间隔和未追踪变量的标准差
Figure DEST_PATH_252165DEST_PATH_IMAGE065
，这里
Figure DEST_PATH_478747DEST_PATH_IMAGE065
就是未追踪到的加速度
Figure DEST_PATH_625432DEST_PATH_IMAGE063
和
Figure DEST_PATH_838239DEST_PATH_IMAGE064
的标准差。每个点计算预测分布时，
Figure DEST_PATH_153814DEST_PATH_IMAGE065
大小相同。为了确定
Figure DEST_PATH_285718DEST_PATH_IMAGE065
的大小，可以尝试多个不同的值，并观察对比效果。一般
Figure DEST_PATH_686743DEST_PATH_IMAGE065
与起飞时加速度大小相近。其中，时间间隔
Figure DEST_PATH_936197DEST_PATH_IMAGE012
取ADS-B数据点时间间隔的平均值。
非线性滤波器：计算采用局部线性化的方法，即扩展卡尔曼滤波（ExpandedKalman Filter，EKF）方法计算过程误差矩阵Q。假设有未追踪到的加速度a和角速度
Figure DEST_PATH_106278DEST_PATH_IMAGE009
，这两个变量的分布也是高斯分布，均值为0，协方差矩阵为
Figure DEST_PATH_DEST_PATH_IMAGE066
。通过CTRA模型的状态转移方程可以求雅可比矩阵
Figure DEST_PATH_815608DEST_PATH_IMAGE067
，计算这两项与状态向量
Figure DEST_PATH_DEST_PATH_IMAGE068
中全部六项的线性相关性。最终计算的
Figure DEST_PATH_140148DEST_PATH_IMAGE062
是a和
Figure DEST_PATH_553812DEST_PATH_IMAGE009
在一次状态转移中对总体协方差的影响。公式如下：
Figure DEST_PATH_578399DEST_PATH_IMAGE069
Figure DEST_PATH_DEST_PATH_IMAGE070
Figure DEST_PATH_130735DEST_PATH_IMAGE071
其中，
Figure DEST_PATH_DEST_PATH_IMAGE072
代表线加速度不确定度的大小，
Figure DEST_PATH_807484DEST_PATH_IMAGE073
代表角速度不确定度的大小。每个点计算预测分布时，
Figure DEST_PATH_900205DEST_PATH_IMAGE072
和
Figure DEST_PATH_44879DEST_PATH_IMAGE073
大小不变。计算
Figure DEST_PATH_564853DEST_PATH_IMAGE067
时用到的状态向量
Figure DEST_PATH_286821DEST_PATH_IMAGE002
均来自于CTRA模型的状态转移方程。
在卡尔曼滤波中，这个情景涉及到了
Figure DEST_PATH_681768DEST_PATH_IMAGE005
与
Figure DEST_PATH_680948DEST_PATH_IMAGE001
大小的问题。而预测时要把
Figure DEST_PATH_637403DEST_PATH_IMAGE062
加到
Figure DEST_PATH_846668DEST_PATH_IMAGE005
上，这个问题实际上转化为了
Figure DEST_PATH_281191DEST_PATH_IMAGE062
与
Figure DEST_PATH_898992DEST_PATH_IMAGE001
相对大小的问题。
Figure DEST_PATH_26348DEST_PATH_IMAGE062
越大，说明模型预测分布的不确定度越大，越倾向于相信观测，但融合的分布越容易被观测数据的噪声影响；而
Figure DEST_PATH_598275DEST_PATH_IMAGE062
越小，越倾向于相信预测，但模型越难以捕捉飞机机动状态的变化。问题的核心就是通过选择合适的
Figure DEST_PATH_102068DEST_PATH_IMAGE062
的大小，在预测和观测之间权衡，在稳定和灵敏之间权衡。
步骤S30，对观测高斯分布和加上过程误差矩阵Q的预测高斯分布进行融合，得到当前时刻的校准高斯分布。
如图4所示，位置的预测分布用圆圈A表示，位置的均值为圆圈A的圆心（从预测的状态向量
Figure DEST_PATH_934895DEST_PATH_IMAGE068
得到），预测高斯分布的协方差矩阵的值越大，圆圈A越大；类似地，位置观测值（也就是ADS-B数据中的值）的分布用圆圈B表示；融合后的分布用圆圈C表示。从图中可以看出，融合后分布的均值更接近于原先方差较小的分布。从图4可以看出，两种信息融合后精度高于原始数据（圆圈更小，也就是方差更小），因此用融合后的位置替换原始数据（观测数据）中的位置，完成对位置的修正。位置更新后，再预测下一个时间点的状态向量，然后再利用下一个状态向量进行更新，两个步骤交替进行。高斯分布的融合用到了
Figure DEST_PATH_DEST_PATH_IMAGE074
时刻预测的状态向量
Figure DEST_PATH_934950DEST_PATH_IMAGE068
和协方差矩阵
Figure DEST_PATH_994173DEST_PATH_IMAGE005
，以及ADS-B中该点的观测位置
Figure DEST_PATH_770499DEST_PATH_IMAGE075
，观测协方差矩阵
Figure DEST_PATH_114891DEST_PATH_IMAGE001
，同时状态向量分布的均值
Figure DEST_PATH_708683DEST_PATH_IMAGE004
和协方差矩阵
Figure DEST_PATH_989623DEST_PATH_IMAGE005
都得到更新。
例如，飞机在跑道上，需要记录飞机每个时刻的位置。一开始飞机在0m处，预测飞机的速度是10m/s，于是1s的时候，预测飞机在10m处（预测飞机在10m处，速度10m/s）。而观测的飞机在20m处，将观测的位置和预测的位置进行融合得到飞机在15m处。也就是说，预测的飞机速度慢了，于是将飞机的速度提高到15m/s（飞机在15m处，速度15m/s）。2s的时候，预测飞机在30m处，而观测的飞机在35m处，接下来又是结合预测和观测的飞机位置修正飞机预测的位置和速度。
步骤S40，根据当前时刻的校准高斯分布更新当前时刻的ADS-B数据的观测高斯分布。
在本发明实施例中，根据当前时刻的校准高斯分布更新当前时刻的观测高斯分布。预测完
Figure DEST_PATH_835219DEST_PATH_IMAGE074
时刻状态向量
Figure DEST_PATH_518005DEST_PATH_IMAGE068
的分布后，还需要利用位置x和y的实际观测值更新状态向量
Figure DEST_PATH_282698DEST_PATH_IMAGE068
的分布，即通过对观测高斯分布和加上过程误差矩阵Q的预测高斯分布进行融合得到当前时刻的校准高斯分布。更新的具体方法可以参考卡尔曼滤波的原理，此处只做示意性介绍。
进一步地，在每次预测时引入自适应机制。当下一个点的时间间隔大于一个阈值时，预测分布的协方差矩阵
Figure DEST_PATH_549469DEST_PATH_IMAGE005
扩大一定倍数（乘以一个常数）。这是因为时间间隔太长，飞机的运动状态会更加地不可预测，因此预测的不确定度需要设置得更大。
进一步地，在每次预测时引入异常点判断机制。当观测的位置与预测位置的马氏距离大于一个阈值时，将该观测标记为异常点，跳过步骤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/s2。
4.位置数据插值
对数据进行初步修正。对于上一步发现的位置异常的点，先去掉x和y的值，再利用前后两点进行线性插值。插值时利用ADS-B中前后两点的位置和三个点的时间戳，计算中间点的位置，替换掉原先x和y的值。
在另一个可选实施例中，对观测高斯分布和加上过程误差矩阵Q的预测高斯分布进行融合得到当前时刻的校准高斯分布之后，还可以包括：使用RTS平滑算法（Rauch–Tung–Striebel smoother）对滤波后的数据进行平滑。滤波是按时间从前往后对数据进行处理，而平滑是从后往前对数据进行处理。平滑后，数据的精度更高，可以替换每个ADS-B点位置的值。
滤波与平滑之后，还需计算以下指标来判断滤波是否发散，即滤波器未能正确追踪飞机的运动状态，偏离了正常航线。一是平均修正距离，即滤波与平滑后，计算每个ADS-B点与原先位置的距离，如果平均值大于某个阈值，认为滤波发散，拒绝使用滤波后的数据。二是计算滤波与平滑后的点距离原先ADS-B轨迹的垂直距离，如果存在偏移较大的点，则拒绝使用滤波平滑后的数据。三是计算轨迹中的尖刺数，如果尖刺数升高，则拒绝使用滤波平滑后的数据。计算指标之后，如果接受滤波与平滑的结果，则将平滑后每一步状态向量
Figure DEST_PATH_667598DEST_PATH_IMAGE068
中的速度值提取出来，补充原先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数据的观测高斯分布，所述上一时刻和所述当前时刻的时间间隔为
Figure DEST_PATH_470469DEST_PATH_IMAGE012
。
所述计算模块20，还用于将飞机的第一个高度大于0的点作为起飞时的点，根据起飞时的点的飞机位置坐标和飞机的历史数据计算初始时刻的ADS-B数据的预测高斯分布。
进一步的，所述获取模块10，具体用于：
将所述当前时刻的观测的经纬度坐标作为所述观测高斯分布中的高斯分布均值；
通过下述公式计算所述观测高斯分布中的协方差矩阵R：
Figure DEST_PATH_406064DEST_PATH_IMAGE050
Figure DEST_PATH_DEST_PATH_IMAGE076
Figure DEST_PATH_566656DEST_PATH_IMAGE052
Figure DEST_PATH_754055DEST_PATH_IMAGE077
其中，
Figure DEST_PATH_536066DEST_PATH_IMAGE054
为高斯分布垂直速度方向的标准差，
Figure DEST_PATH_517928DEST_PATH_IMAGE056
为高斯分布平行速度方向的标准差；
Figure DEST_PATH_526336DEST_PATH_IMAGE059
为上一时刻t的测量误差，
Figure DEST_PATH_750381DEST_PATH_IMAGE006
为飞机速度，
Figure DEST_PATH_527844DEST_PATH_IMAGE008
为飞行偏航角。
在本发明提供的一个实施例中，所述获取模块10，还用于获取历史ADS-B数据，从所述历史ADS-B数据中获取直线飞行数据；
还用于通过对所述直线飞行数据进行直线拟合轨迹，旋转到与x轴平行，得到x方向的方差
Figure DEST_PATH_946187DEST_PATH_IMAGE057
，和y方向的方差
Figure DEST_PATH_769787DEST_PATH_IMAGE060
；
所述计算模块20，还用于计算所述x方向的方差
Figure DEST_PATH_564568DEST_PATH_IMAGE057
和y方向的方差
Figure DEST_PATH_435353DEST_PATH_IMAGE060
的差值得到所述
Figure DEST_PATH_24597DEST_PATH_IMAGE061
；
Figure DEST_PATH_7596DEST_PATH_IMAGE056
除以所述历史ADS-B数据中的速度得到所述
Figure DEST_PATH_340489DEST_PATH_IMAGE059
。
在本发明提供的一个实施例中，若采用线性滤波器，所述状态向量包含飞机的位置、速度
Figure DEST_PATH_325500DEST_PATH_IMAGE006
和加速度a；若采用非线性滤波器，所述状态向量包含飞机的位置、速度
Figure DEST_PATH_820066DEST_PATH_IMAGE006
和加速度a，偏航角
Figure DEST_PATH_24783DEST_PATH_IMAGE008
，角速度
Figure DEST_PATH_489262DEST_PATH_IMAGE009
。
进一步的，所述计算模块20，具体用于：
若采用线性滤波器，状态转移公式为：
Figure DEST_PATH_830245DEST_PATH_IMAGE020
Figure DEST_PATH_259827DEST_PATH_IMAGE021
Figure DEST_PATH_DEST_PATH_IMAGE078
其中，
Figure DEST_PATH_482998DEST_PATH_IMAGE079
和
Figure DEST_PATH_157693DEST_PATH_IMAGE024
是上一时刻t的校准高斯分布中的高斯分布均值和协方差矩阵P，
Figure DEST_PATH_353182DEST_PATH_IMAGE025
和
Figure DEST_PATH_48605DEST_PATH_IMAGE026
是当前时刻
Figure DEST_PATH_460870DEST_PATH_IMAGE027
的校准高斯分布中的高斯分布均值和协方差矩阵P。
关于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）等。
所属领域的技术人员可以清楚地了解到，为了描述的方便和简洁，仅以上述各功能单元、模块的划分进行举例说明，实际应用中，可以根据需要而将上述功能分配由不同的功能单元、模块完成，即将所述装置的内部结构划分成不同的功能单元或模块，以完成以上描述的全部或者部分功能。
以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围，均应包含在本发明的保护范围之内。

Claims (8)

1.一种ADS-B数据的滤波方法，其特征在于，所述方法包括：
获取当前时刻的ADS-B数据的观测高斯分布和ADS-B数据的预测高斯分布，以及过程误差矩阵Q，所述观测高斯分布包括协方差矩阵R，所述预测高斯分布包括协方差矩阵P；
将所述预测高斯分布中的协方差矩阵P加上所述过程误差矩阵Q；
对所述观测高斯分布和加上所述过程误差矩阵Q的所述预测高斯分布进行融合，得到所述当前时刻的校准高斯分布；
根据所述当前时刻的校准高斯分布更新所述当前时刻的ADS-B数据的观测高斯分布；
其中，观测的经纬度坐标符合所述观测高斯分布，预测的状态向量符合所述预测高斯分布；
获取当前时刻的ADS-B数据的观测高斯分布，包括：
将所述当前时刻的观测的经纬度坐标作为所述观测高斯分布中的高斯分布均值；
通过下述公式计算所述观测高斯分布中的协方差矩阵R：
Figure FDA0003410184300000011
σ∥＝vσt
Figure FDA0003410184300000012
R＝ARvAT
其中，σ⊥为高斯分布垂直速度方向的标准差，σ∥为高斯分布平行速度方向的标准差；σt为上一时刻t的测量误差，v为飞机速度，θ为飞行偏航角；取多段历史ADS-B数据，抽出其中直线飞行的片段，用直线拟合轨迹，旋转到与x轴平行，分别得到x和y方向的方差，也就是
Figure FDA0003410184300000013
和
Figure FDA0003410184300000014
二者相减得到
Figure FDA0003410184300000015
σ∥除以ADS-B中的速度得到σt。
2.根据权利要求1所述的ADS-B数据的滤波方法，其特征在于，所述获取当前时刻的ADS-B数据的预测高斯分布，包括：
将上一时刻的校准状态向量输入到状态转移公式得到所述当前时刻的ADS-B数据的预测高斯分布，所述上一时刻和所述当前时刻的时间间隔为Δt。
3.根据权利要求2所述的ADS-B数据的滤波方法，其特征在于，所述方法还包括：
将飞机的第一个高度大于0的点作为起飞时的点，根据所述起飞时的点的飞机位置坐标和飞机的历史数据计算初始时刻的ADS-B数据的预测高斯分布。
4.根据权利要求2所述的ADS-B数据的滤波方法，其特征在于，若采用线性滤波器，所述状态向量包含飞机的位置、速度v和加速度a；若采用非线性滤波器，所述状态向量包含飞机的位置、速度v和加速度a，偏航角θ，角速度ω。
5.根据权利要求4所述的ADS-B数据的滤波方法，其特征在于，所述将上一时刻的校准状态向量输入到状态转移公式得到所述当前时刻的ADS-B数据的预测高斯分布，包括：
若采用线性滤波器，状态转移公式为：
Figure FDA0003410184300000021
xt+Δt＝Fxt
Pt+Δt＝FPtFT
其中，xt和Pt是上一时刻t的校准高斯分布中的高斯分布均值和协方差矩阵P，xt+Δt和Pt+Δt是当前时刻t+Δt的校准高斯分布中的高斯分布均值和协方差矩阵P。
6.一种ADS-B数据的滤波装置，其特征在于，所述装置包括：
获取模块，用于获取当前时刻的ADS-B数据的观测高斯分布和ADS-B数据的预测高斯分布，以及过程误差矩阵Q，所述观测高斯分布包括协方差矩阵R，所述预测高斯分布包括协方差矩阵P；其中，观测的经纬度坐标符合所述观测高斯分布，预测的状态向量符合所述预测高斯分布；
计算模块，用于将所述预测高斯分布中的协方差矩阵P加上所述过程误差矩阵Q；
更新模块，用于根据所述当前时刻的校准高斯分布更新所述当前时刻的ADS-B数据的观测高斯分布；
融合模块，用于对所述观测高斯分布和加上所述过程误差矩阵Q的所述预测高斯分布进行融合，得到所述当前时刻的校准高斯分布；
上述获取模块，具体用于：
将所述当前时刻的观测的经纬度坐标作为所述观测高斯分布中的高斯分布均值；
通过下述公式计算所述观测高斯分布中的协方差矩阵R：
Figure FDA0003410184300000022
σ∥＝vσt
Figure FDA0003410184300000023
R＝ARvAT
其中，σ⊥为高斯分布垂直速度方向的标准差，σ∥为高斯分布平行速度方向的标准差；σt为上一时刻t的测量误差，v为飞机速度，θ为飞行偏航角；取多段历史ADS-B数据，抽出其中直线飞行的片段，用直线拟合轨迹，旋转到与x轴平行，分别得到x和y方向的方差，也就是
Figure FDA0003410184300000031
和
Figure FDA0003410184300000032
二者相减得到
Figure FDA0003410184300000033
σ∥除以ADS-B中的速度得到σt。
7.一种计算机设备，包括存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现如权利要求1至5任一项所述的ADS-B数据的滤波方法。
8.一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现如权利要求1至5任一项所述的ADS-B数据的滤波方法。
CN202110775367.0A 2021-07-09 2021-07-09 Ads-b数据的滤波方法、装置、计算机设备及存储介质 Active CN113257044B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110775367.0A CN113257044B (zh) 2021-07-09 2021-07-09 Ads-b数据的滤波方法、装置、计算机设备及存储介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110775367.0A CN113257044B (zh) 2021-07-09 2021-07-09 Ads-b数据的滤波方法、装置、计算机设备及存储介质

Publications (2)

Publication Number Publication Date
CN113257044A CN113257044A (zh) 2021-08-13
CN113257044B true CN113257044B (zh) 2022-02-11

Family

ID=77190992

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110775367.0A Active CN113257044B (zh) 2021-07-09 2021-07-09 Ads-b数据的滤波方法、装置、计算机设备及存储介质

Country Status (1)

Country Link
CN (1) CN113257044B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US11955014B2 (en) * 2022-01-31 2024-04-09 Honeywell International S.R.O. ADS-B traffic filter
CN115064008B (zh) * 2022-04-14 2023-05-12 中国民用航空总局第二研究所 无人机跑道冲突自主预警系统

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2008076096A (ja) * 2006-09-19 2008-04-03 Mitsubishi Electric Corp 追尾方法及びその装置
CN102082560A (zh) * 2011-02-28 2011-06-01 哈尔滨工程大学 一种基于集合卡尔曼滤波的粒子滤波方法
CN104463841A (zh) * 2014-10-21 2015-03-25 深圳大学 衰减系数自适应的滤波方法及滤波系统
CN104849702A (zh) * 2015-04-30 2015-08-19 中国民航大学 利用ads-b数据的gm-ephd滤波雷达系统误差联合估计方法
CN107589430A (zh) * 2017-09-07 2018-01-16 中国民航大学 基于最小色散方法的ads‑b压制式干扰抑制方法
CN110118556A (zh) * 2019-04-12 2019-08-13 浙江工业大学 一种基于协方差交叉融合slam的机器人定位方法及装置
CN111522044A (zh) * 2020-05-06 2020-08-11 扬州哈工科创机器人研究院有限公司 一种车辆定位方法及装置
CN112362052A (zh) * 2020-10-27 2021-02-12 中国科学院计算技术研究所 一种融合定位方法及系统

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2008076096A (ja) * 2006-09-19 2008-04-03 Mitsubishi Electric Corp 追尾方法及びその装置
CN102082560A (zh) * 2011-02-28 2011-06-01 哈尔滨工程大学 一种基于集合卡尔曼滤波的粒子滤波方法
CN104463841A (zh) * 2014-10-21 2015-03-25 深圳大学 衰减系数自适应的滤波方法及滤波系统
CN104849702A (zh) * 2015-04-30 2015-08-19 中国民航大学 利用ads-b数据的gm-ephd滤波雷达系统误差联合估计方法
CN107589430A (zh) * 2017-09-07 2018-01-16 中国民航大学 基于最小色散方法的ads‑b压制式干扰抑制方法
CN110118556A (zh) * 2019-04-12 2019-08-13 浙江工业大学 一种基于协方差交叉融合slam的机器人定位方法及装置
CN111522044A (zh) * 2020-05-06 2020-08-11 扬州哈工科创机器人研究院有限公司 一种车辆定位方法及装置
CN112362052A (zh) * 2020-10-27 2021-02-12 中国科学院计算技术研究所 一种融合定位方法及系统

Also Published As

Publication number Publication date
CN113257044A (zh) 2021-08-13

Similar Documents

Publication Publication Date Title
Rauch et al. Car2x-based perception in a high-level fusion architecture for cooperative perception systems
CN113257044B (zh) Ads-b数据的滤波方法、装置、计算机设备及存储介质
DE102011120497B4 (de) Systeme und Verfahren zur präzisen Fahrzeugpositionsbestimmung innerhalb einer Fahrspur
US20200216076A1 (en) Method for determining the location of an ego-vehicle
CN109903592B (zh) 一种基于误差理论的高精度航空器自动近地防撞系统地形扫描方法
CN104849702B (zh) 利用ads‑b数据的gm‑ephd滤波雷达系统误差联合估计方法
EP3285245B1 (en) Performance-based track variation for aircraft flight management
EP2370966A2 (en) Methods and system for time of arrival control using time of arrival uncertainty
Kaniewski et al. Estimation of UAV position with use of smoothing algorithms
Peng et al. UAV positioning based on multi-sensor fusion
CN111224710B (zh) 基于卫星空间分布检验的虚拟应答器捕获方法及捕获系统
CN105844969B (zh) 根据气象条件来改进飞行器的飞行轨迹的方法
RU2392198C1 (ru) Прицельно-навигационный комплекс оборудования многофункционального самолета
CN111142143A (zh) 一种基于多源信息融合的进场段飞行技术误差估计方法
US20210158128A1 (en) Method and device for determining trajectories of mobile elements
Herrero et al. Use of map information for tracking targets on airport surface
Breil et al. Multi-agent systems for air traffic conflicts resolution by local speed regulation and departure delay
Du et al. An open data platform for traffic parameters measurement via multirotor unmanned aerial vehicles video
US11579628B2 (en) Method for localizing a vehicle
CN107764273B (zh) 一种车辆导航定位方法及系统
Dasanayaka et al. Analysis of vehicle location prediction errors for safety applications in cooperative-intelligent transportation systems
CN112146615B (zh) 基于多架无人机的边坡监测方法
US20140278181A1 (en) Integrated data registration
CN115775473B (zh) 一种ads-b航空监视系统中的航空器定位系统
CN112154355B (zh) 高精度地图定位方法、系统、平台及计算机可读存储介质

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant