一种结合多重插补的ADS-B数据筛选算法*

邹文华,左 谛,金开研

(民航数据通信有限责任公司 空管数据部，北京 100191)

0 引 言

广播式自动相关监视技术(Automatic Dependent Surveillance-Broadcast，ADS-B)，是利用机载星基导航和定位系统获得定位信息，通过空-空、地-空数据链通信以广播形式发送数据的航空器运行监视技术[1]。ADS-B地面站接收其覆盖范围内的航空器下发的空地DF17格式的ADS-B报文[2]，地面站对此类报文进行整合、组装后送出CAT21格式的ADS-B报文，报文里包含位置、时间、速度、航空器身份等信息[3]。虽然原则上地面站的处理过程遵循ED129系列协议[4]，但是各家厂商对协议内容有不同解读，因此处理过程存在差异，并且不同地面站在硬件配置、运行状态和输送数据的网络延迟等方面也均存在不同。上述原因最终会导致不同地面站送出的报文在数据项完整性、数据质量和数据实时性等方面存在差异，如果不做处理便直接送出，则会出现大量的重复报文、少量的回跳情况，以及各种数据项不全或低质量报文的数据，因此必须对数据进行处理。

随着国内ADS-B的不断发展[5]，在西部的无雷达地区，ADS-B逐渐作为主用且唯一的监视手段给管制提供服务，地面系统需对新监视技术的监视数据进行融合处理。对于传统雷达数据的处理，地面处理系统对雷达目标的融合一般在时间和空间统一后进行加权融合算法，同时结合各种滤波方式，如卡尔曼滤波，消除雷达系统的随机误差，平滑航迹[6]。在有雷达数据参与融合的算法一般均采用与多雷达融合类似的思想，即加权和滤波结合的方法[7-9]。而在无雷达监视且新监视技术作为主用且唯一监视手段的地区，使用雷达区域的处理方法处理新监视数据并不合适，特别是对于ADS-B技术。ADS-B与雷达相比，有以下两个本质的不同：一是ADS-B报文中最重要的位置信息由机载GPS获得后下发，任何一个地面站获得报文后解出的位置误差理论上应该都一样，与航空器距离地面站的距离无关，不需要对位置进行加权计算；二是机载下发的GPS位置信息本身已经在机载端进行了滤波处理[10-11]，ADS-B地面系统解析输出的位置信息平滑度较高。因此，多地面站覆盖时对多ADS-B数据的融合不用套用传统的多雷达融合算法。

为此，本文利用ADS-B报文中的位置信息具有唯一确定且相对较平滑的特点，采用选优的方式从多条报文中挑选质量最优的ADS-B报文作为筛选后的数据，对于筛选后不合格的数据，根据实际情况，在数据严重缺失或条件极度不理想的条件下，可采用多重插补法提升其完整性。通过对实际多个航班ADS-B数据的算法验证，证实了算法能有效去除重复和回跳报文，并优化提高数据整体质量。

1 航迹建立

本算法为周期型输出的筛选算法，间隔以1 s为例，在筛选程序启动后，收集每一秒内的所有报文，查看是否已经建立航迹，如果未建立航迹，则先建立航迹；如果已经完成航迹的建立，则进入筛选程序。对于尚未建立航迹的报文，应用以下方法确立可信赖的航迹。

航迹建立和处理是将同一目标的点迹连成航迹的处理过程，高效且可靠的航迹建立是保障航迹稳定更新的前提。一般直接起始法的基本原理是，对于连续收到的目标位置信息Pi(i=1,2，…，N)，如果其中有M(M

图1 航迹建立流程图

Step1 收集N条ADS-B报文，记为{REPORTi}(i=1,2，…，N),N为报文数量，为自定义参数，一般取大于3的数，本文中取5。每条报文中含位置(即经纬度)、时间和速度。

Step2 采用倒叙遍历方式，即时间值从大到小顺序，从{REPORTi}中挑选时间差大于0.4 s且小于4 s的3条报文，记为{RPi}(i=1,2，…，N)。0.4 s和4 s这个两个时间差选择的标准主要基于业务考虑，ADS-B位置更新的频率为0.4～0.6 s[2]，因此0.4 s ADS-B位置更新的最小时间间隔，如果前后报文时间小于0.4 s，则有理由相信这两条报文中的位置可能是相同的，不适宜用来做航迹建立判断。雷达更新的频率为4～12 s。如果前后两个位置的时间差达到了ADS-B自身更新频率的10倍，接近雷达的更新频率，则认为此报文串连续性有欠缺，不适宜用来做航迹建立判断。

Step3 对{RPi}(i=1,2,3)中的3条报文进行互相验证。建立阶段要满足3条连续报文都是可靠的，方法是位置和速度相匹配，且不能超出最大位置差，即判定3条报文是否满足下述公式：

(1)

(P1与P2点的计算方向差-P2与P3点的计算方向差)<120°。

(2)

式中:Vmax为航空器最大可能速度。

如果3条报文满足式(1)和式(2)，则航迹建立成功，输出RP1作为航迹起始点；如果3条报文不满足式(1)和式(2)，则航迹建立不成功，重新遍历寻找满足Step 2的3条报文，找到后执行Step 3，如果遍历完毕没有找到满足Step 2的3条报文，则说明当前{REPORTi}队列中N个报文无法建立航迹，需要释放清空此队列，重新收集报文进行航迹建立过程，即回到Step 1。

2 ADS-B数据筛选算法

2.1 多约束条件下的数据筛选

航迹建立成功后，按第1节中的两个筛选步骤进行融合。首先挑选进入融合的报文。满足公式(3)和(4)的报文可进入2.2节进一步筛选：

Tn+1<([Tn]+1),

(3)

Tn+1>(Tn+0.4)。

(4)

式(3)和(4)中的Tn表示第n个融合周期优选输出的点的报文时间(Time of Message Reception for position，TMR)，[Tn]为第n个融合周期融合后的点的报文中TMR值的整数，Tn+1表示第n+1个融合周期收到的报文中的TMR时间。

公式(3)中的1 s，是基于本文融合算法的周期为1 s；公式(4)中的0.4 s，即进入融合的报文时间必须比上一条报文时间大0.4 s，原因见1.1节中关于0.4 s的说明。

按照上面挑选标准，如果只挑选出了一条报文，且质量达到要求，则作为第n+1个融合周期融合点输出；如果质量未达到要求可采用多重插补法(Multiple Imputation,MI)对缺失数据进行补全达到质量要求后输出。如果有多条报文，则进入2.2节，进行优选融合。

用图2举例说明。U为报文接收时间轴，T为报文TMR时间轴，假设第一个周期输出的融合点为U1/T1点，下个融合周期的时间范围为[[U1]+1,[U1]+2)，[U1]为U1值的整数，那么在这个融合周期时间段内，这架飞机共收到4条报文，这4个点的接收时间和TMR时间(格式记为U2/T2)相对上个融合输出点分别是(U1+0.8)/(T1+0.1)、(U1+1.0)/(T1+1.5)、(U1+1.2)/(T1+1.6)和(U1+1.4)/(T1+0.6)，再根据式(3)和式(4)进行挑选。

图2 进入融合程序的报文挑选示例图

假设上一个点的融合点U1/T1的U1值为0.4，T1值为0.1，那么可以得出(U1+0.8)/(T1+0.1)点不满足式(3)和式(4)，(U1+1.4)/(T1+0.6)不满足式(3)，(U1+1.0)/(T1+1.5)和(U1+1.2)/(T1+1.6)满足条件，可按照2.2节对这两个点进行进一步质量评估，挑选出质量高点，作为此融合周期的融合点输出。

2.2 根据数据质量指标的筛选

2.1节中，第n+1个融合周期挑选出来多条报文，如何从这些报文中挑选出最优的报文，则需要从多方面进行评估。本文选取两个方面评估，分别是报文中数据项的完整性和报文中位置的质量。数据项的完整性是报文中用户关心的数据项的存在情况，用户关心的数据项越丰富，说明报文的数据项完整性越好。报文位置的质量由航空器下发报文中携带的位置数据项不确定性(Navigation Uncertain Category，NUC)进行评价。NUC表征的是ADS-B报文中下发的位置的精度和完好性，取值范围为0～9，数值越大，位置的精度和完好性即数据质量越高。通过式(5)计算参与融合的点的质量水平参数值Qi，挑选Qi值最大的点作为融合后的点:

Qi=Witm×Qitm+Wnuc×Qnuc。

(5)

式中：Qitm为数据项完整性标识参数；Qnuc为报文位置质量标识参数；Witm和Wnuc分别是数据项完整性标识参数和报文位置质量标识参数的权重，两个权重值的和为1。

通过计算所有报文的Qi，选取Qi值最大的报文作为此融合周期的优选融合点输出。

2.2.1Qitm的计算公式

在[Ti,Ti+1)时间段内，参与优选融合的某个报文k的数据项完整性参数Qitm(i,k)由报文中重要数据项的存在性确定。根据民航局空管局在2018年颁布的《民航空管ADS-B及S模式雷达数据应用实施计划(修订版)》[12]中确定的第一阶段中对数据的要求，并结合航空公司在应用过程中的实际需求，本文选择气压高度、航班号、二次代码、地速和时间5个数据项的有无表征完整性分项参数，各个参数根据专家评级法取合适的权重，但各分项权重之和需等于Witm，即

θ1+θ2+θ3+θ4+θ5=Witm。

(6)

式(6)中的θ1、θ2、θ3、θ4和θ5分别是气压高度、航班号、二次代码、地速和时间5个数据项的完整性分项标识参数。

2.2.2Qnuc的计算公式

在[Ti,Ti+1)时间段内，参与融合的某个报文k的NUC参数Qnuc(i,k)的计算公式如下：

Qnuc(i,k)=NUCk/9。

(7)

由于NUC是一个0～9的整数，需要做归一化处理后用于表征ADS-B报文质量参数。

3 数据极度不理想情况下的补全

对于ADS-B报文中的各项指标，根据用途要求综合专家评级法会有不同的权重需求，根据不同的权重要求，可以算出不同要求下的Qi，根据应用场景以及引用要求的差异，结合历史数据确定阈值M，阈值M的确定需结合航空器飞行空域内的环境及流量共同确定。在筛选后，若所筛选Qi>M，则判定报文数据合格；若所筛选Qi

3.1 基于MI算法对缺失数据进行补全

多重填补方法分为三个步骤:

Step1 为每个空值产生一套可能的填补值，这些值反映了无响应模型的不确定性；每个值都被用来填补数据集中的缺失值，产生若干个完整数据集合。

Step2 每个填补数据集合都用针对完整数据集的统计方法进行统计分析。

Step3 对来自各个填补数据集的结果进行综合，产生最终的统计推断，这一推断考虑到了由于数据填补而产生的不确定性。该方法将空缺值视为随机样本。

3.2 基于DA方法的MI插补方法

对于上述筛选过程，如果数据质量达不到要求，则需要对数据进行插补处理。考虑到ADS-B数据类型特点，即报文之间彼此独立，每条报文之间的特征值也是彼此独立，所以采用基于DA的MI算法对数据进行插补。数据扩张(Data Augmentation，DA)算法是马尔科夫链蒙特卡洛(Markov Chain Monte Carlo，MCMC)算法之一，特点在于可以处理任意缺失模式[14]。

DA算法可以通过以下步骤实现:

从Pr(Ymis|Yobs,θ(t))中抽取Ymis(t+1), 从Pr(θ|Yobs,Ymis(t+1))中抽取θ(t+1)，重复该过程多次，这样就建立了一条马尔科夫链：Ymis(1)，θ(1)，Ymis(2)，θ(2)……而且收敛于P(Ymis,θ|Yobs)。其中，θ表示条件概率的参数，Yobs代表观测值，Ymis代表缺失值，P(Ymis,θ|Yobs)则表示在观测值已知的情况下对缺失值进行预测的条件分布。DA算法估计的目的是从收敛的分布中随机抽取Ymis值，替代缺失数据。

3.3 基于DA方法的MI插补方法有效性论证

在通过文献查询和理论分析论证了插补算法在ADS-B数据补全方面的适用性之后，本节通过使用人工生成模拟数据的方法验证插补方法的有效性。选取地面站系统中2019年6月20日航班号为CES5810的1 000个完整ADS-B数据，人工产生缺失率分别为5%、10%、20%、30%、40%、50%，缺失类型为完全随机缺失的缺失数据，再使用基于DA方法的MI插补方法对缺失数据进行插补。而在上述文本中选择的气压高度、航班号、二次代码、地速和时间5项参数中，航班号、二次代码属于定性数据，由飞行员输入和设置，无法通过已有数据建立预测分布，而其中的时间项由于对数据精度要求过高，不适用于用插补法进行插补，故仅对气压高度、地速进行插补。

表1总结比较了原数据与插补后数据的均值、标准差以及皮尔森相似度[15]，用以衡量插补前后数据的相似度。在R语言中完成上述过程[16-17]。由表1可知，插补后的数据的均值与方差与原数据相比均没有太多差异，而插补后的数据与原数据的皮尔森相似度随着缺失率的增长而缓慢下降，但仍然保持着较高相似度，证明了基于DA方法的MI插补方法的有效性。

表1 插补前后对照表

4 筛选算法的仿真过程

4.1 仿真流程

本文取本公司地面站系统中2019年6月19日3个航班的所有地面站数据进行融合，验证融合及选优效果。所选航班信息、跨时、接收到报文的地面站数量见表2。

表2 模拟数据信息表

取气压高度、航班号、二次代码、地速和报文时间5个数据项的有无表征报文完整性的分项参数，取NUC值的归一化后表征报文质量参数，各参数的权重设定如表3所示。

表3 权重参数表

表3中权重值在满足局方文件中对数据项的存在性要求的条件下，兼顾了用户的实际使用需求。NUC、气压高度、航班号、二次代码的权重较高，NUC和气压高度权重高是对位置精度、可靠性和高度存在性的要求，体现了监视技术的核心是四维航迹信息。航班号和二次代码是表征航空器身份的信息，在管制系统中是ADS-B航迹与飞行计划进行相关的关键信息，重要性不言而喻。报文时间项和地速项存在性相对较弱，但也是体现了ADS-B数据自身的优势。ADS-B数据更新率高、及时性强并自身携带高稳定性的速度信息，无需地面系统计算，优势明显。

4.2 仿真结果

根据第2节的参数和权重对这3架航空器进行组合优选融合处理，对融合前后报文数量、报文NUC质量平均值、数据完整性指标、整体质量进行评估，结果如表4所示，可以确定算法去除了大量重复的报文，并且整体质量参数有提高。

表4 仿真结果

结合该运行场景对数据的要求，选取阈值M为0.85，由表4可知，经过筛选后，除地址码为780E81的航班外，另外两个航班的Qi均大于M，符合质量指标要求，而对于地址码为780E81的航班，则需使用上述MI插补方法进行补全。使用上述DA方法的MI插补方法对极端不理想数据地址码为780E81的航班进行插补。

4.3 插补过程模拟

基于DA算法的算法原理，需要基于已有数据生成预测分布，如3.3节中所述，在上述文本中选择的气压高度、航班号、二次代码、地速和时间5项参数中，航班号、二次代码属于定性数据，由飞行员输入和设置，无法通过已有数据建立预测分布，可由人工填写，故依旧仅对气压高度、地速进行插补。地址码为780E81的航班的ADS-B数据(气压高度、地速)的数据缺失情况如图3所示。

图3 航班 “780E81”ADS-B数据缺失情况

由图3可知，地址码为780E81的航班的ADS-B数据的缺失类型为完全随机缺失，且由飞机航行以及气压高度测量原理可知，气压高度与地速之间存在线性关系，所以在DA算法中使用线性回归预测模型，基于已有值，建立预测分布进行插补，插补结果如图4所示，可见经过插补气压高度及地速的缺失率下降为0，Qi值显著增加至0.890 9，且大于阈值M，符合数据质量要求。

图4 插补前后缺失情况对比图

5 结束语

本文所提算法可有效解决不理想条件下ADS-B数据中的重复及回跳现象，确定数据质量指标以快速对报文数据进行筛选，在恶劣环境下亦有良好应用。对于缺失数据，使用多重插补法进行补全，极大地提高了数据完整性。

但同时，在上述流程中，无论是筛选还是插补依旧存在信息的浪费以及误差，想要从根本上解决数据的质量问题，仍应该不断发展航空通信导航技术，提高地面站及航空器下发数据质量。