融合航空器监视数据的航空气象参数反演及验证

2022-07-22卢晓光安晨晖许忠睿

信号处理 2022年6期

卢晓光安晨晖张喆许忠睿韩萍

（1.中国民航大学天津市智能信号与图像处理重点实验室，天津 300300；2.中国商飞上海飞机制造有限公司，上海 201324）

1 引言

航空气象参数与飞机的飞行状态参数密切相关，进而严重影响飞行安全和航空公司的正常运营［1-3］，其中主要包括风、温度、压强和密度等［4］。目前常规的航空气象探测手段有气象气球、气象雷达、气象卫星和飞机气象数据中继系统（Aircraft Meteorological Data Relay，AMDAR）等［5-7］，这些探测手段具有大尺度、依赖设备构成和准确度依赖布网密度等特点，其依赖气象气球的发射密度、雷达地面站和卫星的设备建设情况等。如何获得准确实时的气象信息一直是航空气象人员关注的焦点。近年来，由于Mode S 二次监视雷达（Secondary Surveillance Radar，SSR）和广播式自动相关监视（Automatic Dependent Surveillance-Broadcast，ADSB）的快速部署和广泛应用［8-10］，欧洲有大量学者研究了基于这些航空器监视数据的气象参数获取方法。利用航空器监视数据获取气象参数具有很多优势，其直接来自飞机现有的传感器和系统，不会产生额外的费用，数据采集范围遍布全球，且由于数据更新率高因而时间分辨率较高，因此利用航空器监视数据获取气象参数是经济可行的方案。2012 年，Haan 等人使用荷兰航空交通管制局提供的Mode S 数据进行气象参数的反演［11］。2013 年，Leege 等人提出一种基于ADS-B 数据反演风矢量的方法，该方法只在飞机进行大角度转弯时才具有较高的精度［12］。2014 年，Hurter 等人利用飞机的飞行轨迹提取风参数，该方法反演风矢量需要从多架在不同方向飞行的飞机轨迹中提取，飞机转弯、爬升和下降时风矢量的反演性能会下降［13］。2017 年，Sun Jun-zi 等人利用ADS-B 和Mode S 消息反演气象参数，其没有进行数据质量分析，采用欧洲中期天气预报中心（European Centre for Medium-Range Weather Forecasts，ECMWF）气象信息对比验证［14］。2020 年，刘涛等人提出一种基于标准偏差的风矢量估计方法，该方法提高了飞机中小角度转弯时风矢量的估计精度，但没有利用其他气象信息验证反演的风矢量［15］。综上所述，基于单一数据来源或飞行轨迹的气象反演方法一般只在飞机转弯时可用或不能用于所有的飞行阶段，且绝大多数Mode S 数据是空管部门提供的，这些数据不会开放给个人或企业使用，不利于开展更加深入的研究。采用航空器监视数据进行气象反演需要解决的问题包括数据源、数据质量控制、高精度反演、反演结果的验证等方面。

基于此，本文使用OpenSky 网络搜集的开源航空器监视数据，给出一种基于ADS-B 和Mode S 融合数据的气象参数反演方法，并将其应用于单航线和区域性场景。最后，将反演的气象参数与AMDAR和ECMWF 气象信息对比验证，证明了方法的准确性。文中提到的数据质量分析与融合、气象参数的反演等内容，具有十分重要的实际应用价值。同时，本方法能够提供实测雷达无法布置区域的风速、风向、温度等信息。

2 数据采集与预处理

2.1 数据获取

为了获取开源的ADS-B 和Mode S 消息，并将其用于气象参数的反演，本文利用OpenSky 网络中的航空器监视数据监测气象信息，系统构架如图1 所示。OpenSky 网络是一个航空器监视数据开源计划，是基于社区的非营利性接收者网络，旨在收集航空器监视数据，并将收集到的数据免费提供给第三方［16］。这些数据来源于具有S模式接收机的众多志愿者，自2013 年以来不断收集ADS-B 和Mode S等数据，实际接收范围如图2所示。

截止到目前，OpenSky 网络收集了来自全球3000多个传感器的超过20万亿条ADS-B、Mode S和空中防撞系统（Traffic Collision Avoidance System，TCAS）消息，拥有庞大的空中交通监视数据集。OpenSky网络的信号接收机有时会提供错误信息或不能获取飞机发送的所有信息，因此，应用OpenSky网络数据前需要进行数据质量分析并剔除异常值。

使用国际民用航空组织（International Civil Aviation Organization，ICAO）地址、时间和高度为基准，融合ADS-B 和Mode S 数据。实际使用中，Mode S信号会出现错误，错误的信号可能会使ICAO 地址错误或产生一个新的ICAO 地址。例如，如果一位发生移位，则奇偶校验可能会导致ICAO 地址略有不同，由于国际民航组织的地址通常是按顺序分配的，因此该错误地址可能与其他飞机的ICAO 地址相同，容易产生混淆。为了避免这些问题发生，利用ADS-B 信息验证Mode S 消息，剔除异常信息。对于气象反演，需要ADS-B 和Mode S 融合后的数据，单独发送ADS-B 或Mode S 消息的飞机会被剔除，所以ADS-B 消息验证Mode S 不会剔除有用信息。将具有相同ICAO 地址和时间的ADS-B 和Mode S信息组合在一起，如果两信息中的高度、地速等参数的差异低于设定的阈值，则认为该条Mode S 消息是正确的。

2.2 数据质量分析与异常值处理

OpenSky 网络数据主要集中在欧洲和美国等地区，且欧洲地区的数据质量较好，因此主要采用欧洲地区的数据进行气象参数的反演与验证，如图2所示。通过OpenSky 网络，获取了2020 年6 月5 日至2020 年6 月14 日共10 天的ADS-B 和Mode S 消息，数据接收范围主要集中在德国、荷兰和法国等地区。因数据量过大，每天接收的范围有些许不同，但都集中在欧洲地区，图3 显示了10 天内途径数据接收设置范围内的飞机数量，该数量以飞机的ICAO地址码为基准统计。

为了进一步了解飞机发送的ADS-B 和Mode S消息的数量，统计了10天内每天的ADS-B和Mode S消息数量，如图4 所示。由图4 可知，Mode S 消息数据量为ADS-B 数据的15 倍左右，同时接收到的Mode S 消息中Mode S 增强监视（Enhanced Surveil⁃lance，EHS）占比约为70%，这表明欧洲Mode S EHS的使用较广泛。Mode S EHS 消息中包含真空速、指示空速、马赫数和磁航向等，这些数据被用于监视飞机的速度、高度、位置等，同时还是反演气象参数的关键信息。

基于OpenSky 网络的低成本信号接收机有时会提供错误信息或不能获取飞机发送的所有消息，这为信息的解码和气象参数的估计带来困难。马赫数和温度等参数不会短时间内剧烈变化，对于飞机的单个参数，随时间推移应用滑动窗口滤波算法。首先剔除明显无效数据，对于窗口的每个位置，检查窗口的中心点是不是离群点，计算窗口中点的中值，并计算该中值与数据中心点的差值，当这个差值大于量化步长的设定倍数时，将其中心值认定为异常值，如图5 所示。大多数异常值是由错误的格式引起的，此时该消息中的其他字段很可能也是错误的，因此直接丢弃异常值的整条消息。ADS-B 和Mode S 中包含不同的数据项，对解码后的明码信息进行数据质量分析，具有以下条件的任何消息都将被丢弃：（1）缺少关键项：经度、纬度、气压高度、地速、目标身份码、真空速和磁航向等；（2）异常数据项：气压高度异常、几何高度异常、几何和气压高度差值异常、航迹点丢失和航班号异常等。

对2020 年6 月11 日一天的ADS-B 数据进行分析，共得到458 万条有效消息，251 万条缺少关键项消息，27 万条异常数据项消息，如图6 所示。其中约30%的数据缺失了数据项，进一步统计表明大部分数据缺失了气压高度、速度和磁航向等，而这些数据项是反演气象参数的关键，因此这些数据将被丢弃。异常数据中大部分为气压高度数值异常和几何和气压高度差值异常，由文献［17］可知几何和气压高度的差值大部分在650 ft 之内，若几何和气压高度差值超过650 ft，认定该数据为异常数据。

3 气象参数反演算法

3.1 基本思想

气象参数反演平台的目的就是通过OpenSky网络收集ADS-B 和Mode S 消息，将网络获取的消息进行解码、融合、质量分析等处理，并利用融合后的信息反演出风矢量、温度、压强等气象参数。为了确定风矢量，需要有真空速矢量和地速矢量；为了估计空气温度，需要有马赫数和真空速；为了估计压强，需要有飞机的气压高度，这些所需的数据项及其来源如表1 所示。Mode-S 信号中的各种飞机状态和意图信息由56 位的Comm-B数据选择器（Binary Data Selector，BDS）决定，其包含在信号的消息字段中，由两个十六进制字符表示。对于气象参数的估计，若ADS-B 和Mode S 信息包含相同的数据项，此时优先使用ADS-B 信息。因为对于Mode S 解码具有很大的难度，我们采用概率解码器进行解码，ADS-B 信息具有更好的稳定性和确定性，且具有更高精度的高度信息，这对获取高精度的温度和压强等参数至关重要。

表1 ADS-B和Mode S的可用数据项Tab.1 Available data items of ADS-B and Mode S

3.2 压强和密度反演

飞机上通常使用气压高度表示飞行高度，来自飞机上的气压高度计。根据预先设定的参考压力值计算垂直距离，工作原理是测量飞机周围的大气压数值，然后将数值转换为飞机气压高度。压强和高度的关系可以表示为：

其中，g为重力加速度，hb、Pb和Tb分别为大气层的基本高度、气压和温度，λ为固定常数，如表2所示。

表2 国际标准大气层参数Tab.2 International standard atmospheric parameters

对式（1）进行反写，即可得到大气压强的计算公式：

其中，R为特定气体常数，R=287.05 J/(kg ⋅K)，其他参数的具体数值由飞机的实际飞行高度决定，如表2所示。

空气密度可以由理想气体定律得出，此时密度计算公式可以表示为：

其中，ρ为密度，单位为kg/m3，T为局部空气温度，单位为K。

3.3 空气温度反演

可以使用真空速、指示空速和马赫数计算温度，马赫数是真空速与空气中声速的比值：

其中，c为声速，大小为340 m/s，Vt为真空速，单位为m/s。

声速不是恒定的，主要取决于温度，还受到湿度和其他因素的影响，但它们的影响要小得多。得到马赫数、真空速和指示空速后，温度可以表示为：

其中，T为温度，单位为K，Vi为指示空速，单位为m/s，M为马赫数，没有单位，a0为海平面声速，大小为340.3 m/s。

除了上述方法外，还可以通过几何高度、气压高度和高度差值来推算温度值。此方法需要两架飞机或一架处于上升、下降的飞机来确定某一高度层的大气温度。ADS-B 数据中包含气压高度和几何高度，文献［17］给出了几何和气压高度的差值，可以利用高度差值确定大气层温度，此时温度可以表示为：

其中，hp是距国际标准大气的压力高度，单位为m，a是几何和气压高度的绝对差值，单位为m，p0为海平面气压，数值为101325 Pa，T0为海平面温度，数值为288.15 K。

3.4 风矢量反演

空速矢量是地速与风速矢量的差值，可通过地速矢量和真空速矢量组成的速度矢量模型间接推导出风速。当飞机遇到风时，爬升和下降阶段的飞行轨迹角很小，飞机的速度分量在水平方向最大，垂直分量可以忽略不计。因此速度模型可以简化为仅包含真空速、地速和水平风矢量三个参数的矢量关系模型，此时地速矢量是风矢量和真空速矢量的总和，如图7所示，数学模型可表示为：

其中，W为风矢量，Vg为地速，Vt为真空速。χa为磁航向，是真空速与磁北方向的夹角；χg为轨迹角，是地速与地理正北方向的夹角；χω为风矢量与地理正北方向的夹角。真空速可以从BDS 50中直接获得，也可以从BDS 60中的指示空速和马赫数计算得出：

其中，T为温度，单位为K，Vt和Vi分别为真空速和指示空速，单位为m/s，M为马赫数，没有单位，γ为比热容比，γ=cp/cv=1.397774。

根据航空学基本原理，任何矢量0°都指向地理正北方向，角度沿顺时针方向增加。当没有风时，飞机将精确地沿着预期的航向角飞行，真空速等于地速；当存在风时，飞机会偏离预期或计划的轨迹，此时实际的飞行方向偏离预期航向角。在数学上，地速矢量减去真空速矢量即可得出风矢量。因此，只需获得地速、轨迹角、真空速和磁航向就可以直接获取风速和风向参数。为了化简式（7）的计算，矢量W、Vg、Vt分解为东西分量和南北分量，此时风矢量可以表示为：

其中，Vx为风的东西分量，Vy为风的南北分量。

在Mode S 信号中，磁航向的值存在偏差，此偏差由设备测量误差或局部磁场引起误差等原因造成，偏差值一般与飞机的型号相关。飞机使用磁偏角查询表将真实航向从惯性参考单元（Inertial Ref⁃erence Unit，IRU）转换为磁航向［18］。使用地面飞机的ADS-B 和Mode S 消息来确定磁偏角，地面的飞机不受风的影响，此时飞机的磁航向等于地面的轨迹，地面飞机为机场中起飞、着陆或滑行的飞机。磁地理偏差是地理真北与磁北之间的差，为时间、经度、纬度和海拔高度的函数，如图8所示。真实的磁航向可表示为：

其中，∠MT为校正后的磁航向，∠M为Mode-S消息解码后的原始磁航向，∠MG为磁北与地理真北的夹角，∠MD为飞机设备的磁航向偏差。

4 实验结果与分析

4.1 单飞机气象参数验证

4.1.1 基于ADS-B数据的风矢量验证

ADS-B 数据中只包含地速矢量，不包含真空速矢量，且地速和真空速一般不同时存在，为此，文献［12］、［15］等采用圆拟合方法来计算风矢量。为了验证该风矢量反演方法的适用范围和准确性，通过三组ADS-B 数据进行验证分析。三组实验分别选取了2020 年6 月11 日的49°大角度转弯飞行、20°中角度转弯飞行和8°小角度转弯飞行三种情况进行实验分析，分别记为实验1、实验2 和实验3。得到的风矢量结果如图9 所示，蓝色圆点为输入的地速分量，黑色的虚线为拟合圆，红色边框的绿色点为圆心，风矢量的计算结果如表3所示。

表3 反演风矢量与ECMWF风矢量对比Tab.3 Comparison between inverse wind vector and ECMWF wind vector

通过以上三组实验可以看出基于ADS-B 数据风矢量反演方法的精度受到飞机转弯角度的影响。当飞机的转弯角超过40°时，反演的风速和风向与ECMWF 风矢量基本相似；当飞机的转弯角小于20°时，反演的风速和风向与ECMWF 风矢量信息相差较大；当飞机的转弯角小于10°时，风矢量的反演结果与实际情况严重不符。而在实际情况下，大部分飞机的转弯角度都小于20°，因此基于ADS-B 数据的风矢量反演方法具有很大的局限性，只能用于极少的航班，且只能用于进行大角度转弯的飞机。

4.1.2 基于融合数据的风矢量验证

为了验证气象反演算法的适用性，单飞机气象参数验证分别选取了不同地区和机型的两架航班。实验4选取途经德国地区的波音B763客机，实验5选取途经比利时地区的空客A319 客机。实验4 为2020 年7 月8 日02 时03 分至25分起飞至巡航阶段共计22 分钟的飞行数据，飞行高度为1875 ft 至36025 ft，将实验4飞机标记为飞机4。为了验证气象反演算法的正确性，选取飞机4 相同时间段的AMDAR 数据，来源于文献［19］，共计41 个数据点。飞机4 气象参数反演结果与AMDAR 数据验证如图10所示。

图10 中红色点为风速、风向和温度的估计值，蓝色点为飞机4 传输的AMDAR 气象数据。由图可知风速和风向的估计曲线与AMDAR 数据体现出较好的一致性，温度的估计曲线与AMDAR 温度完全吻合，估计结果精度较高。实验结果表明气象参数反演算法针对单个飞机的反演结果是可行的。为了更加直观的观察估计风速的变化趋势与AMDAR风速的差异，采用四次多项式拟合算法对风速估计值进行拟合平滑处理，如图11所示。拟合后的风速值曲线更加平滑，降低了风速值瞬间波动对风速估计的影响，提高了风速估计的稳定性。

实验5选取2020年7月9日10时04分至20分共16 分钟的飞行数据，飞行高度为1050 ft 至36050 ft，将实验2 飞机标记为飞机5。同时选取飞机5 相同时间段内的AMDAR 气象数据，共计48 个数据点，飞机5 气象反演结果与AMDAR 数据验证如图12所示。

由图12 可知，风速和风向的估计值与AMDAR数据吻合度较高，没有出现较大的偏差，温度估计值与AMDAR 温度值完全吻合，反演结果真实可信。通过实验4 和实验5 可知，基于ADS-B 和Mode S 融合数据的气象参数反演方法对单架航班具有较高的估计精度，能对飞机飞行提供帮助。同时该方法对温度的估计精度极高，温度反演曲线与AMDAR温度值变化趋势完全一致。

4.2 多飞机气象参数验证

多飞机实验选取指定地区和时间段内的所有航班进行区域性的气象参数反演，并将反演的风速、风向和温度与ECMWF 气象信息对比验证，数据来源于文献［20］。为了验证气象反演算法的适用性，多飞机气象参数验证分别选取了只包含陆地区域和同时包含陆地和海洋区域的两组实验。多飞机气象参数反演具有极大的数据量，为了更好的体现该地区的气象信息，采用箱型图对反演的气象参数进行展示。利用箱型图的统计学特性［21］，可以更好的表示出特定区域内气象信息的分布情况。实验6 选取2020 年6 月11 日途经德国、法国、比利时等陆地区域上午5：30 至6：30 一个小时内所有飞机的ADS-B 和Mode S 数据，经数据解码、质量分析和融合后，共得到12 多万条监视数据，如图13所示。

通过上述气象反演算法得出该地区的气象参数，并将风速、风向和温度信息与该地区上午6时的ECMWF 气象信息对比验证，如图14 所示。图中每一个高度层都绘制了一个箱型图，每个箱型图从左至右分别为最小值、上四分位数、中位数、下四分位数与最小值，黑色加号为异常数据。箱体中的黄色短线为估计气象参数的中位数，绿色短线为平均数；蓝线折线和黄线折线分别为ECMWF 数据中风速、风向和温度的平均值和中值。

由图14 可知，ECMWF 风速和风向曲线大部分处于箱型图的箱内，证明反演的该地区的风矢量信息与ECMWF风矢量基本吻合；ECMWF温度曲线全部在箱型图的箱内，证明反演的温度信息与ECMWF 温度完全吻合。为了更加直观的展示该地区风速和风向的具体数值和分布趋势，绘制了该地区不同高度和经纬度的风杆图，如图15所示。图中横纵坐标分别为纬度和经度，每隔2 km 进行一次风矢量参数的统计，绿色圆圈表示该地风速为零。

为了验证气象参数反演算法各气象要素日变化的敏感性，分别选取了实验6 相同地区2020 年6 月11 日11∶30-12∶30 和17∶30-18∶30 两个小时的航空器监视数据进行验证，气象参数反演结果和ECMWF气象信息对比验证结果，如图16、图17所示。

通过图15、图16、图17 可知，2020 年6 月11 日5∶30-6∶30、11∶30-12∶30 和17∶30-18∶30 早中晚三个时间段的ECMWF 风速、风向和温度曲线绝大多数处于箱型图的箱内，反演的风速、风向和温度信息与ECMWF 气象参数吻合度极高，证明该气象参数反演方法具有较高的日变化敏感性。

实验7选取2020年6月12日途经英国、德国、比利时和荷兰等海洋和陆地区域上午12∶30 至13∶30一个小时内所有飞机的ADS-B 和Mode S 数据，经数据解码、质量分析和融合后，共得到18 多万条监视数据，如图18所示。

将反演的风速、风向和温度等信息与该地区下午13 时的ECMWF 气象信息对比验证，如图19 所示。由图19 可以看出ECMWF 风速和风向曲线大部分处在箱型图的箱内，证明该地区的风速和风向反演结果与ECMWF 风矢量信息基本吻合；ECMWF温度曲线全部在箱型图的箱内，证明反演的温度信息与实际温度吻合度极高。为了更加直观的展示该地区风速和风向的分布趋势，绘制了该地区不同高度和经纬度的风杆图，如图20所示。

通过两组实验可知，风速随着高度的升高逐渐增大，高度达到30000 ft以上时，风速开始逐渐减小；风向随着高度升高逐渐趋于稳定，波动范围更小，这是因为高空风比较稳定。温度随着高度的升高不断降低，达到36000 ft以上时，逐渐趋于稳定，这是因为对流层顶的温度比较稳定。本气象参数反演方法对多架飞机区域性气象参数反演具有较高的估计精度，且不受陆地和海洋等区域的影响，温度估计精度极高，反演的气象参数能对特定区域内的气象信息进行补充，可用于气象预报，进而保障飞机的飞行安全。