航空器的高度保持性能是航空器对所指定飞行的高度层的符合性能力，反映了由于测高系统误差 (Altimetry System Error，ASE) 和指定高度偏差 (Assigned Attitude Deviation，AAD) 导致的垂向飞行误差。由于不同航空器具有不同的ASE，并且ASE会随机龄的增长或航空器的加改装产生变化，因此国际民航组织 (International Civil Aviation Organization，ICAO) 要求各国民航必须对本国民用航空器的ASE进行定期评估，避免由于航空器ASE过大导致实际飞行高度严重偏离规定高度，从而给航路飞行带来安全隐患[1]。
目前，国际民航通行的航空器高度保持性能评价技术包括机载监控和地基监控2种。机载监控基于与航空器机载设备独立的GPS接收机采集的航空器真实几何位置信息，需要采用跟飞的方式采集数据，由于其自身采集限制，使用机载监控，无法实现批量采集数据，难以实现对整个空域内航空器ASE的批量分析和性能评估。从而，目前欧洲、美国、日本等国家和地区都使用地基监控系统对其空域内的航空器进行高度保持性能监控。传统的地基监控实质上是由4~5个地面站组成，可以覆盖一定空域范围的广域多点定位系统，利用地基监控系统，可以明显改善机载监控系统数量严重不足的问题，但是传统的地基监控系统也存在成本高、建设周期长等不足[2]。
Lebedev[3-4]于2002年提出将广播式自动相关监视 (ADS-B) 作为解算航空器高度保持性能的数据源，进而提出使用临近多航空器高度比较的方式求解航空器ASE值的方法。ADS-B通过航空器自主广播自身位置实现地空和空空监视。利用ADS-B信息的航空器几何高度位置对航空器ASE进行分析，可有效解决传统机载监控和地基监控方式的不足。Aldis和Nixon[5]对Lebedev[3-4]提出的方法进行了进一步分析，提出了其方法在实际使用过程中面临的挑战，说明了其解算精度无法满足ICAO对于航空器高度保持性能监控的要求。基于这一情况，美国FAA的Falk和Martin[6-8]进行了进一步论证，证明了ADS-B数据与机载监控数据源所采集的GPS数据为同源数据，即使用ADS-B下发的航空器几何高度信息作为航空器空域运行的真实高度，并开始与澳大利亚合作使用这一方法验证其实际应用效果。基于实际数据分析，Butcher等[9]发现在澳大利亚空域内，不同地理位置的ADS-B站点测定的ASE值有较大差异。Falk等[10-11]又进一步论证了航空器基于ADS-B的几何高度信息与航空器的大地水准面基准有关，即对于基于椭球高 (Height Above Ellipsoid，HAE) 和平均海平面 (Mean Sea Level，MSL) 高度下发的航空器几何高度信息之间存在差异，并可通过高程异常值加以修正，但是对于航空器的基准面信息在ADS-B下发的数据中并不包含。基于以上结论，Aldis等[12]提出了在基准面信息缺失的情况下，使用回归分析对基于累积数据对航空器基准面进行确定的方法，并在澳大利亚推广使用。
目前，FAA和澳大利亚直接使用ADS-B数据中包含的航空器几何高度信息作为ASE数据源，并使用相同的机载监控流程进行解算。为了保证数据源可以准确反映真实的航迹点高度位置特征，使数据精度满足航空器高度保持性能评价的要求，航空器高度保持性能解算的关键步骤是用恰当的平滑方法对航空器高度信息进行处理，并对航迹点高度保持性能结果数据进行统计分析，而相关方法国外鲜有报道。从而为了解决中国基于ADS-B数据的航空器高度保持性能监控的需要，提出一种适用于ADS-B的平滑方法，同时为了更准确地分析航空器的高度保持性能，在得到其航迹点的ASE值以后，需要利用累积数据对单架航空器的统计分布规律进行分析。
因此，本文在现有研究成果的基础上，针对ADS-B数据精度特征提出基于核平滑的高度数据平滑方法，并结合实际数据比较所提出方法与现行方法的优劣。同时基于航迹点的ASE值提出基于混合正态分布模型的航空器ASE分布拟合分析方法，并给出这一分析方法的实例验证。
1 基于ADS-B的ASE评估 ASE是航空器真实高度和测高系统显示高度之间的高度差，主要由机载测高系统测量不准确导致，是高度保持性能评估要解决的主要问题。
航空器下发的ADS-B数据中包含来自测高系统的气压高度和来自GPS接收机的GPS高度。其中，GPS高度的精度远高于气压高度，可以用于表征航空器的真实气压高度，从而得到航空器ASE。图 1给出了基于ADS-B的ASE评估的一般步骤。图中：TVE为总垂直偏差。
 |
| 图 1 ASE计算的一般步骤 Fig. 1 General steps of ASE calculation |
| 图选项 |
为了对航空器的高度保持性能进行评价，对于获取到的ADS-B数据，应包含航空器的四维位置及高度信息，即样本数据中应包含时间向量Ti，纬度向量Lat, i和经度向量Long, i，GPS几何高度向量Hg, i和模式C高度向量Hmc, i，i=1, 2, …, N，N为原始样本个数。ASE求解的一般步骤如下：
1) 根据样本集{Ti, Hmc, i}推断航空器运行的标准高度层FL，并截取平飞片段{Tm, Hmc, m}，m=1, 2, …, NL，NL为平飞片段的样本个数。将Lat, i、Long, i和Hg, i按{Tm, Hmc, m}进行同步，得到Lat, m、Long, m和Hg, m。
2) 根据样本集{Lat, m, Long, m}计算航空器平飞片断的经纬度相关系数C，以C≥95%为条件截取航空器运行的平直飞片段{Lat, k, Long, k}，k=1, 2, …, NLS，NLS为平直飞片段的样本个数。将Tm、Hg, m和Hmc, m按{Lat, k, Long, k}进行同步，得到Tk、Hg, k和Hmc, k。
3) 对样本集{Tk, Hg, k}、{Tk, Hmc, k}分别进行数据平滑预处理，得到{Tk, Hsg, k}和{Tk, Hsmc, k}。
4) 根据{Tk, Hsmc, k}和FL计算向量AAD, k。
5) 引入气象数据，根据{Tk, Lat, k, Long, k}计算航空器实际运行时给定的FL对应的基于MSL的几何高度向量HgMSLfl, k。
6) 引入高程异常数据，根据{Tk, Lat, k, Long, k, HgMSLfl, k}计算得到基于HAE的几何高度数据HgHAEfl, k。
7) 利用HgMSLfl, k、HgHAEfl, k和Hsg, k计算向量TVEMSL, k和TVEHAE, k。
8) 利用TVEMSL, k、TVEHAE, k和AAD, k计算向量ASEMSL, k和ASEHAE, k。
9) 对ASEMSL, k和ASEHAE, k进行分布模型拟合，得到航空器高度保持性能的最终评价结果。
在整个评价过程中有2个关键步骤，即步骤3) 对航迹高度数据的平滑处理和步骤9) 对航迹点ASE结果的分布拟合分析。
2 基于核回归的ADS-B高度数据平滑 根据航空无线电技术委员会 (Radio Technical Commission for Aeronautics，RTCA) 制定的标准，航空器下发的ADS-B数据中的GPS高度和气压高度的精度均为25 ft (1 ft=0.304 8 m)。因此，需要基于所采集的ADS-B数据对航空器的GPS高度和气压高度进行准确回归，减小航空器高度时间序列的随机噪声，使最终结果更好地反映航空器的运行趋势。
在ADS-B高度数据回归中，由于自变量和因变量之间具有未知的相关关系，因此使用非参数回归平滑方法。本文使用核回归平滑方法实现GPS高度和气压高度的准确计算。
回归曲线拟合的目标是发现变量
 | (1) |
式中：mn(x) 为m(x) 的核估计；K(·) 为核函数；hn为核估计的窗宽。
对于核函数，通常使用的有：均匀核K(u)=
图 2给出了3种不同的核函数对于平滑效果的比较。
 |
| 图 2 不同核函数的比较 Fig. 2 Comparison of different kernel functions |
| 图选项 |
然而研究表明，相对于核函数的选择，平滑结果的精确度更应该说是窗宽选择的问题[13]，核估计的窗宽对平滑程度有较大影响。窗宽较大的情况下，每个观测值的影响区间较大，则结果更为平滑，过大的窗宽往往会造成“过平滑”的情况，反之窗宽过小的情况会造成平滑程度的缺乏，导致“平滑不足”[1, 14]。从而解决ADS-B数据平滑问题的关键点是为核平滑方法找到一个适合的窗宽。图 3给出了不同窗宽平滑效果的比较。图中：h为窗宽。
 |
| 图 3 不同窗宽的比较 Fig. 3 Comparison of different bandwidths |
| 图选项 |
本文使用Bowman和Azzalini[15]在1997年建议的最优窗宽，即
 | (2) |
式中：median表示中位数；n为总样本量。
图 4给出了本文方法与FAA平滑方法的比较情况。可以看出，较FAA的平滑结果，使用核平滑方法具有更强的稳定性，平滑结果更好。
 |
| 图 4 2种平滑方法的结果比较 Fig. 4 Results comparison of two smoothing method |
| 图选项 |
基于核回归平滑，可以对ADS-B数据中的高度数据进行适当的平滑，进而得到航空器航迹近似真实的GPS高度和气压高度。
3 航空器ASE分布拟合 由于GPS高度的精度远高于气压高度，因此可将GPS高度作为气压高度的真实值，二者之差为航空器测高系统测量气压高度的误差，即ASE。通过对每架航空器航迹点ASE的分布密度估计，可以评估航空器高度保持性能。
对基于航空器ADS-B数据计算的ASE，采用核密度估计方法估计其概率分布。对于航空器的ASE样本数据，可以认为其为取自连续分布p(x) 的数据
 | (3) |
本文算例中的ASE核密度估计结果如图 5所示。可以看出，航空器相对于ASEMSL，ASEHAE的分布模型较为分散，从而可以进一步确定这一航空器的基准面为MSL。
 |
| 图 5 航空器ASE的核密度估计 Fig. 5 Kernel density estimation of aircraft ASE |
| 图选项 |
除此之外，航空器ASE的核密度分布拟合的一项重要功能是可以对航空器ASE的部分其他特征进行描述。在此以另外一个实例对这一问题进行分析，图 6给出了这一分析的实例。
 |
| 图 6 某航空器基准面与日均值分析 Fig. 6 Height reference and daily mean value analysis of an aircraft |
| 图选项 |
图 6(a)为航空器的日均值分析，横轴为航空器的运行时间。可以看出，日均值回归线斜率约为0.06，并没有明显的趋势性变化，应该用总体的均值作为对该航空器性能的评价，但是可以发现图 6中的ASE日均值点明显较为分散。
从航空器实际运行的角度来看，每一架航空器都具有2套独立的测高学设备，在飞行过程中由飞行员自行确定使用哪套测高学设备进行航空器高度控制。从而可以初步预计，发生这种ASE日均值点离散的原因是该航空器的2套测高学设备具有不同均值的高度保持性能。显然，对于具有此类特性的航空器，简单地用均值对其性能进行评价无法真实反映其高度保持性能，所以，在对其性能进行评价时，应考虑其他方法，更好地反映其离散特征。
本节需要解决的实际问题为：如何从核密度估计的结果中发现离散情况的真实原因，并对真实均值进行分析。
从统计学角度来看，航空器2套独立的测高学设备应具备不同的高度保持性能分布特征，反映在总体的分布上则为2个分布的混合。从而这一实际问题对应的数学问题为对ASE进行基于混合分布的模型拟合，其难点则在于对混合分布中各参数的估计。
利用2种正态分布的混合分布对航空器ASE进行分析，即
 | (4) |
式中：μ1和μ2为2个正态分布的均值；σ1和σ2为2个正态分布的标准差；α1和α2为2个正态分布的权重系数。进而通过对该分布的拟合，估计其参数，得到的μ1和μ2分别代表2个独立的测高学设备所对应的ASE值。
在此，对2种正态分布的混合分布的分布参数进行求解，分布拟合的初值根据表 1确定。表中：sx为样本标准差。
表 1 个体航空器分布拟合初值 Table 1 Initial values of individual aircraft distribution fitting
| 参数 | 拟合初值 |
| α1 | 0.5 |
| α2 | 0.5 |
| μ1 |  |
| μ2 |  |
| σ1 |  |
| σ2 |  |
表选项
图 7给出了混合模型式 (4) 的拟合结果。图 7(a)中，紫色和绿色密度曲线为混合分布的2个拟合结果，图 7(b)中的红色密度曲线为2个正态分布概率密度的混合。从结果中可以看出，拟合结果很好，2个正态分布的均值相差100 ft左右。
 |
| 图 7 航空器ASE混合分布拟合 Fig. 7 Mixed distribution fitting of aircraft ASE |
| 图选项 |
图 8给出了该航空器的核密度估计情况。
 |
| 图 8 航空器ASE核密度估计的双峰结果 Fig. 8 Double peak results for kernel density estimation of aircraft ASE |
| 图选项 |
可以看出，经混合分布拟合后，该航空器的ASEMSL核密度曲线具有明显的双峰特征，体现了2套独立的测高学系统之间航空器高度保持性能的差异。
经过实际分析可以发现，利用基于混合分布的模型拟合和参数估计，可以有效地对ASE离散的问题进行判别，从而准确估计航空器的ASE性能。
对于具有ASE双峰特性的航空器的评价，不能简单使用航空器ASE的均值作为其评价方法，根据混合分布模型拟合的结果，应使用偏差绝对值较大的分布均值作为航空器最终的高度保持性能评价结果。
4 实际应用 本文基于ADS-B数据的航空器高度保持性能分析方法已经在中国RVSM航空器高度保持性能的实际分析中正式使用。
利用这一方法，中国地区监控组织已经对飞经中国空域同时具有ADS-B out能力的航空器进行高度保持性能分析。截至2016年1月，共累计结算航迹点ASE两亿八千多万条，对累计近7 000架航空器的高度保持性能进行了分析。分析结果满足了ICAO的相关要求，成为全球范围内认可的航空器高度保持性能分析结果。
在目前正根据中国民航ADS-B总体规划建设的“ADS-B航空器高度保持性能地基监控系统”中，以本文给出的分析方法作为主要核心算法加以运用，图 9给出了系统的初步运行界面。
 |
| 图 9 ADS-B航空器高度保持性能地基监控系统界面 Fig. 9 Graphical user interface of ground-based ADS-B aircraft height keeping performance monitoring system |
| 图选项 |
5 结论 本文在介绍了航空器高度保持性能监控的背景下，给出了如何利用ADS-B数据进行航迹点的ASE求解，同时给出了利用基于ADS-B数据的航空器航迹点ASE值对航空器总体高度保持性能进行评价的不同方面。
由于利用ADS-B技术对航空器高度保持性能进行分析属于一个新兴的方法，所以本文中所给出的方法和分析流程为首次提出，结果表明：
1) 基于核回归的ADS-B高度数据平滑方法可以较好地对精度为25 ft的ADS-B高度数据进行平滑，进而得到航空器航迹近似真实的GPS高度和气压高度。
2) 利用基于混合分布的模型拟合和参数估计，可以有效地对ASE的双峰离散的问题进行判别，从而准确估计航空器的ASE性能。
3) 利用本文所提出的方法，可以对基于ADS-B数据计算的航空器高度保持性能进行评价，给出个体航空器性能分析结论。
参考文献
| [1] | ICAO.Manual on implementation of a 300 m (1000 ft) vertical separation minimum between FL 290 and FL 410 inclusive:Doc 9574-AN/934[S].2nd ed.Montreal:ICAO, 2002:9-12. |
| [2] | ICAO.Operating procedures and practices for regional monitoring agencies in relation to the use of a 300 m (1000 ft) vertical separation minimum between FL 290 and FL 410 inclusive:Doc 9937-AN/477[S].Montreal:ICAO, 2010:74-76. |
| [3] | LEBEDEV B.Method of ASE monitoring:ICAO-SASP-WG/WHL/1-WP/07[R].Canberra:ICAO SASP, 2002:2-4. |
| [4] | LEBEDEV B.Another method of ASE monitoring:ICAO-SASP-WG/WHL/2-WP/10[R].Montreal:ICAO SASP, 2002:1-5. |
| [5] | ALDIS G, NIXON D.Use of ADS-B for height-keeping monitoring:ICAO-SASP-WG/WHL/12-WP/22[R].Santiago:ICAO SASP, 2007:12-13. |
| [6] | FAA.Investigation into the use of automatic dependent surveillance-broadcast data for monitoring aircraft altimetry system error:ICAO-SASP-WG/WHL/13-IP/13[R].Montreal:ICAO SASP, 2008:3-4. |
| [7] | MARTIN L, FALK C.An update to the investigation into the use of automatic dependent surveillance-broadcast data for monitoring aircraft altimetry system error:ICAO-SASP-WG/WHL/14-WP/17[R].Paris:ICAO SASP, 2008:6-12. |
| [8] | FALK C.An update to the investigation into the use of automatic dependent surveillance-broadcast data for monitoring aircraft altimetry system error:ICAO-SASP-WG/WHL/15-IP/06[R].Montreal:ICAO SASP, 2009:4-5. |
| [9] | BUTCHER R, FALK C, ALDIS G.Trial-use of the automatic dependent surveillance-broadcast data for monitoring aircraft altimetry system error:ICAO-SASP-WG/WHL/16-WP/24[R].Auckland:ICAO SASP, 2009:8-9. |
| [10] | FALK C, ALDIS G, BUTCHER R.An update to the investigation into the use of automatic dependent surveillance-broadcast data for monitoring aircraft altimetry system error:ICAO-SASP-WG/WHL/17-WP/21[R].Montreal:ICAO SASP, 2010:19-21. |
| [11] | FALK C, ALDIS G, BUTCHER R.Progress on the research conducted to determine the use of automatic dependent surveillance-broadcast data for monitoring aircraft altimetry system error:ICAO-SASP-WG/WHL/18-WP/12[R].Brussels:ICAO SASP, 2010:2-5. |
| [12] | ALDIS G, BARRY S, FALK C, et al.Large-scale study of the use of automatic dependent surveillance-broadcast data for monitoring aircraft altimetry system error:ICAO-SASP-WG/WHL/19-WP/16[R].Montreal:ICAO SASP, 2011:23-26. |
| [13] | 薛留根. 现代非参数统计[M].北京: 科学出版社, 2015: 173. XUE L G. Modern nonparametric statistics[M].Beijing: Science Press, 2015: 173.(in Chinese) |
| [14] | JEFFREY S S. Smoothing methods in statistics[M].Berlin: Springer, 1996: 54-55. |
| [15] | BOWMAN A W, AZZALINI A. Applied smoothing techniques for data analysis[M].New York: Oxford University Press, 1997: 31. |
