基于期望最大算法的空间事件及异常值探测*

刘劲宏 吴晨韵 徐劲 杜建丽 雷祥旭

1(重庆交通大学智慧城市学院 重庆 400074)

2(中国科学院紫金山天文台 南京 210023)

3(山东理工大学建筑工程学院 淄博 255000)

0 引言

美国空间监测网（Space Surveillance Network，SSN）在全球布设了29 个测站，能够连续跟踪观测空间碎片，监测数据通过轨道处理后以两行根数（TLE）的格式进行轨道编目进入ELSET 数据库，并同时在网站*www.space-track.org发布，是可公开使用的种类最全、数量最多的空间碎片轨道编目数据库[1-3]。TLE 数据在空间活动研究中具有重要作用，例如利用TLE 数据反演热层大气密度、估计弹道系数、再入预报和碰撞预警等[4-8]。然而许多研究表明，TLE精度存在较大波动的现象，轨道根数出现离散程度不一的异常值以及与趋势相违背的连续TLE 数据等[9,10]。由于TLE 精度未提供，直接使用TLE 数据将引起潜在风险或降低后续研究成果的稳定性和可靠性。因此，必须识别ELSET 数据库中的异常TLE数据以及空间事件。

目前，ELSET 数据库中异常TLE 的产生机理尚未被完全理解，并且不同研究者仅关注特定轨道区域的少部分空间碎片，造成异常数据的清理方法和流程不同。针对低轨碎片和大椭圆轨道碎片，Águeda等[11]和Dolado-Perez等[12]将TLE 轨道根数转化为分点根数，使其变化平滑，采用迭代最小二乘多项式拟合轨道参数的变化趋势，将离散明显的异常值剔除。针对存在轨道机动的卫星，Patera[13]提出了滑动窗口-多项式拟合方法，探测卫星轨道机动，通过计算轨道参数与其回归值的离散度确定是否存在变轨，然而离散度越小，越难区分轨道机动和空间环境扰动产生的异常值。针对地球同步轨道的火箭体，Lidtke等[3]利用滑动窗口和多项式提出了多步清理策略，给出多种异常值阈值的定义，以便清理不同的轨道参数。然而以上方法无法统一轨道参数异常值的清理流程和算法，适用的空间碎片对象有限，并且在实际处理过程中存在异常值判定不合理等问题。

为克服以上问题，Liu等[14]提出了基于期望最大（Expectation Maximization,EM）算法的滑动窗口-多项式拟合预报技术，解决了现有算法存在的关键问题，重新定义了基于高斯统计分布特征的阈值，统一了轨道参数的异常值探测流程，适用于所有轨道区域空间碎片异常值清理。探测结果表明：对于相对误差幅度 ≥ 5×10-3的平运动和偏心率异常值，可以实现100%的探测；对于相对误差幅度 ≥0.01的轨道倾角异常值，可实现74%～98%的探测；对于相对误差幅度 ≥5×10-3的轨道倾角异常值，可实现37%～68%的探测。

在实际清理过程中，由于异常值没有真实的参考值，对TLE 中探测的异常值进行严格评估通常是不可能的，这是所有清理算法常见的问题。本文在上述工作基础上，基于EM 算法的滑动窗口-多项式拟合预报技术，识别受复杂空间环境及轨道机动影响的空间碎片，给出该方法的清理步骤和最终清理效果。

1 TLE 异常值清理算法原理

TLE 轨道根数异常值清理算法基于极大验后估计得到多项式的拟合-预报表达式，通过引入EM 算法精确地估计误差的期望和方差，该算法涉及的变量包括阈值、预报值个数和滑动窗口大小（滑动窗口包含的TLE 数据数量）[14]。

1.1 基于极大验后估计的线性回归和预报表达式

假设观测值y可表示为如下多项式：

根据极大验后估计原理[15]，如果给定一系列观测值，对未知参数进行估计，其极大验后估计的期望和方差分别为

相应预报值的期望和方差表达式如下：

1.2 异常值及其空间事件探测原理

图1 给出了基于EM 算法的异常值及其空间事件探测流程[14]。假定滑动窗口内观测值总数始终为NTLE=9，对滑动窗口内的变化趋势进行拟合，并剔除离散度较大的值（窗口内阈值设定）；i为窗口内最后一个TLE 观测值，假定同时预报窗口外k=6 个观测值，i+1,...,i+6 为窗口外观测值时刻的预报值，则根据窗口外阈值设定，预报结果有以下三种情况。

图1 基于EM 算法的异常值及空间事件探测流程图（黑色圆点代表TLE 观测值，红色圆点代表预报失败值）Fig.1 Flowchart of outlier detection based on EM algorithm (Black dot represents the TLE observation value,and the red dot represents the prediction failure value)

情况1预报值与观测值的差全部落入置信区间内，滑动窗口向前移动k=6，继续向前预报。

情况2预报值与观测值的差全部落在置信区间之外，形成新的子序列且滑动窗口以第i+1 个实测值为起始点，重新初始化滑动窗口。

情况3预报值与观测值的差部分落在置信区间之外（i+1 和i+6 ），第i+1为异常值，滑动窗口向前移动至第i+5个观测值，重新初始化滑动窗口。

通过以上三种情况的区分，可以探测异常值和空间事件，且使用滑动窗口外多个预报值同时进行判断，使得空间事件的探测更加合理。经过以上方式的滑动窗口前移，TLE 序列中明显的异常值被清理，而次级别的异常值可能依然存在于TLE 序列中，因而整个清理过程需要重复多次。

2 TLE 异常值探测算法应用

2.1 阈值设定策略对异常值探测的影响

Liu等[14]讨论了异常值及空间事件探测的影响因素，包括滑动窗口大小、滑动窗口和预报值的阈值设定策略、预报个数以及多项式拟合阶数。研究表明，滑动窗口大小在60～100 内，预报个数为8，多项式阶数为3 时较为合适。阈值策略包括以下三种。

策略1滑动窗口内阈值和预报阈值均采用 3σ。

策略2滑动窗口内阈值采用2σ，预报阈值采用3σ。

策略3滑动窗口内阈值和预报阈值均采用 2σ。

碎片NORAD ID 为13025，处理时间范围为1982 年1 月9 日至1988 年11 月20 日，共计6.9 年，对轨道倾角采用三种阈值策略清理，结果如图2 所示[14]。策略1 适合探测异常值少、变化平滑的TLE 子序列，产生的异常值数量最少，而部分变化复杂的子序列仍然存在异常值；策略2 产生的异常值数量适中；策略3 产生的异常值数量最多，对于变化平滑的子序列，策略3 过于严格，其对于变化复杂的子序列却很适合。因此，在实际清理过程中，三种策略需交替使用，首先采用策略1 清理，将变化平滑的子序列分离，然后采用策略2 或者策略3 对复杂的子序列做进一步清理。

图2 采用三种阈值策略清理轨道倾角，碎片NORAD ID13025（Ariane 1 R/B）Fig.2 Inclination outlier detection by three strategies for NORAD ID13025 (Ariane 1 R/B)

2.2 TLE 异常值清理流程

针对轨道参数（平运动、偏心率、轨道倾角和Bstar）中存在的异常值，介绍异常值清理算法的实际步骤，并展示部分碎片清理结果。在无特殊说明时，清理算法涉及的参数预报个数默认为8，多项式阶数为3，其他参数如滑动窗口大小、阈值策略，根据实际情况进行选择。选用卫星HAIYANG 2 A（NORAD ID：37781），该卫星存在轨道机动，且经历的空间环境包括复杂和平滑变化时期，适用于对异常值清理算法的实际清理流程进行说明，处理的TLE 数据时间范围为2018 年1 月1 日至2021 年9 月8 日。清理步骤如下。

第1步清理平运动序列中的异常值，滑动窗口大小为200，阈值采用策略1，清理1 次，序列中存在的显著异常值完全被清理，清理结果如图3 所示。判定的TLE 子序列起始点和终止点，用于辅助清理TLE 中的异常值。

图3 平运动序列中的异常值清理结果Fig.3 Outlier detection in mean motion

第2步清理偏心率序列中的异常值，滑动窗口大小为200，阈值采用策略1，连续清理2 次。清理结果见图4，共计产生16 个子序列（编号1～16），其中第1，3，4，6 和14 子序列存在明显异常值，应进一步对其进行异常值清理。子序列1 包括406 个观测值，采用策略1，滑动窗口大小分别用40 和80 各清理1 次。子序列3 和4 合并为一个子序列，共包括467 个观测值，采用策略2，滑动窗口大小为100，连续清理2 次。子序列6 包括200 个观测值，采用策略1，滑动窗口大小为40，清理1 次。子序列14 包括289 个观测值，采用策略1，滑动窗口大小为40，清理1 次。各个子序列清理的异常值结果见图5，偏心率异常值最终的清理结果见图6。

图4 偏心率序列中的异常值初步清理结果Fig.4 Preliminary results of outlier detection in eccentricity series

图5 各个子序列的异常值清理结果Fig.5 Outlier detection results of each subsequence

图6 偏心率序列中的异常值最终清理结果Fig.6 Final results of outlier detection in eccentricity series

第3步清理轨道倾角序列中的异常值。轨道倾角序列中存在大量多个（8～30）连续观测值近乎水平直线的分布。采用2 阶多项式拟合，窗口大小20 进行清理，异常值清理结果如图7 所示。

图7 轨道倾角序列中的异常值清理结果Fig.7 Results of outlier detection in inclination

第4 步清理Bstar 序列中的异常值，滑动窗口大小为80，阈值采用策略1。清理结果如图8 所示，由于Bstar 变化没有规律，仅清理了3 个异常值明显的观测值。

图8 Bstar 序列中的异常值清理结果Fig.8 Results of outlier detection in Bstar

经过以上步骤，TLE 序列中的异常值得以清理。在清理过程中需注意以下几方面问题。

（1）异常值存在复杂性，清理一次未必能成功，需重复多次。

（2）对轨道根数初次清理时，滑动窗口大小宜选用较大值（100～200），便于拟合长期趋势变化，随后根据判定的子序列进一步清理。

（3）子序列的观测值数量一般较少，滑动窗口大小宜选用较小值（20～40），若子序列中出现异常值明显（子序列6）或部分趋势相违背（子序列3，4）的情况，宜选用较大值（100～200）。

（4）针对变化平滑的子序列，阈值采用策略1，而变化稍微复杂的子序列需要采用策略2。

将爱情与诗歌相捆绑，更加完美地展现诗人对“爱情”的赞美和欣赏是该诗的思维修辞。“爱情”是文学作品永恒的主题，《你眼睛的弧线》更是毫无例外。艾吕雅曾有一首诗集命名为《爱情·诗歌》(L’amour La poésie)，题名“爱情·诗歌”明显蕴含了两层含义：爱情与诗歌；爱情即诗歌。可以看出，在作者心中，爱情与诗歌虽不相同，但却是难解难分。诗人无法离开“爱情”的滋养，他与爱情的关系就如同“白日取决于空气的纯净”，他甚至可以将自己的生命都托付于爱情(我的血液也在你的眼神中流淌)。可以说艾吕雅的整部作品、他的一生都在为爱而生，为爱奋斗。

（5）针对变化特殊的子序列，例如轨道倾角，其多项式阶数和窗口阈值需作调整。

2.3 TLE 异常值清理结果

从ELSET 数据库中选取4 个碎片，展示异常值清理结果（见图9）。从图9 可知，TLE 序列中并非所有轨道根数都包含异常值，碎片27386 是卫星Envisat，存在轨道机动，所探测的异常值数量最多，其次是卫星O3BFM9（NORAD ID：40351）。碎片12446是TITAN 34 B 的火箭体残骸，异常值主要来源于轨道倾角，碎片23358 是USA 40 R/B 的碎片，雷达反射面积等级为SMALL，所探测的异常值数量最少。

图9 空间碎片12446，23358，27386 和40351 轨道倾角异常值清理结果Fig.9 Outlier detection of inclination for debris 12446,23358,27386 and 40351

3 结论

针对TLE 序列中存在的大量异常值，提出一种基于期望最大算法的滑动窗口-多项式拟合预报方法，给出清理步骤及部分碎片清理结果，得到结论如下。

（1）对于长周期TLE 序列的清理，滑动窗口大小应尽量大（100～200），阈值采用策略1，实现异常值的清理和子序列的判定，并进一步对子序列进行逐个清理。

（3）若TLE 序列包含与趋势变化相违背的明显的连续多个异常值，采用滑动窗口100～200 和阈值策略1 将其清理。

（4）对于变化平缓的TLE 序列，阈值采用策略1；对于变化稍微复杂的TLE 序列，阈值采用策略2。

（5）由于异常值的存在影响多项式的拟合与预报，整个过程需要反复清理，降低异常值的影响，一般1～3 次。

由于TLE 异常值产生的原因复杂，在实际清理中需结合TLE 序列的变化特点，合理选择参数组合，实现清理结果最优。为便于清理，在基于EM 算法的滑动窗口-多项式拟合预报技术基础上开发了可视化异常清理软件，可通过重复性清理实验找到最佳清理参数设置和最佳清理效果。作者利用卫星（ID:36985）2011 年的TLE 数据进行轨道机动探测，初步发现异常清理后，能显著降低其他因素的干扰，提升轨道机动探测成功率，更易识别小推力轨道机动。因此，后续将在已有研究基础上对观测值实时异常判定及轨道机动探测开展更深入的研究。