摘要
以两行轨道根数为数据来源，提出采用逆向移动滑窗的自适应轨道机动检测方法。建立检测窗口的自适应配置模型以获得不同目标的检测窗口，采用交叉弧段预报思想更精确地估计机动时刻，通过逆向移动滑窗方法实现近实时机动检测。将检测结果与真实机动以及已有文献中的检测结果进行比较，结果表明：新方法能够适应不同轨道空间目标机动检测的窗口配置，且在给定仿真条件下检测成功率可以达到90%以上，机动时刻估计准确度与已有方法对比有显著的提升；仿真结果表明新方法可以应用于近实时的轨道机动检测。
Abstract
Using two line elements as data source, an adaptive method using antidromic window moving for detecting maneuvers of space objects was proposed. Adaptive window configuration model was established to get detection window for different objects. Meanwhile, cross arc prediction was integrated in order to obtain the maneuver time more precisely. Applying antidromic window moving method, the maneuver in the latest epochs can be detected. The results were compared with real maneuvers as well as in the existing literature. It is shown that the proposed method can adapt to the window configuration of maneuvering detection of space objects on different orbits. Besides, the detection success rate can reach more than 90% under the given simulation conditions. Moreover, the accuracy of maneuver time estimation is observably higher than that of existed method. Simulations of the approach proposed have demonstrated good performance in near-real-time maneuver detecting.
在不间断的太空发射活动下，地外空间累积的在轨航天器数量庞大，且有不断增加的趋势。在这样的空间环境下，对卫星的机动行为进行检测是提高己方航天器在轨安全性，对非合作目标进行威胁评估和行为识别的重要手段，能为空间态势感知提供有力支撑。
对于卫星的历史机动检测，分析两行轨道根数（two line elements，TLE）包含的轨道信息变化是有效的策略之一。目前基于TLE进行卫星机动检测的方法主要分为移动窗口曲线拟合法和轨道预报误差拟合法。
移动窗口曲线拟合法方面，刘二江等[1]结合SGP4/SDP4轨道预报模型，采用马氏距离判别法对航天器是否发生轨道机动进行识别。张栩晨等[2]提出基于平经度及其漂移率的地球同步轨道（geosynchronous orbit，GEO）卫星东西机动周期提取和预报方法，实现对赤道带GEO非合作目标的机动检测。张炜等[3]使用双行根数作为数据源，对双行根数进行预处理，检测轨道机动事件并剔除误差较大根数，对天宫一号姿态变化进行判定。崔红正等[4]针对非合作目标，结合地基与天基观测数据，提出不同推力作用下轨道机动检测策略。Li等[5]采用离散小波变换方法对TLE数据进行降噪处理，计算半长轴和轨道倾角的检测阈值进行机动检测。Patera[6]通过在一定长度的时间序列内移动窗口得到轨道误差，从统计意义上设定异常检测门限进行机动检测。Lemmens等[7]分析了同一空间目标任意TLE数据集之间的一致性，结合背景模型进行机动检测。Liu等[8]构造了一种基于期望最大化算法的滤波器，通过滤波器实现对航天器机动的识别。移动窗口曲线拟合法通过对窗口内的编目数据进行多项式拟合，将窗口内最后一个数据点的拟合值与实际值的差作为预报误差。此类方法基本从数据处理入手，无法有效地对卫星轨道变化的物理规律进行分析，且检测精度还不够高。
轨道预报误差拟合法方面，张涛涛等[9]选择低地球轨道（low Earth orbit，LEO）的半长轴或轨道倾角作为特征轨道参数，基于预报偏差进行机动检测。许晓丽等[10]对不同类型轨道预报误差的演化规律和特征进行了分类讨论，对基于TLE的航天器机动检测有潜在实用价值。李涛等[11]对特征轨道参数的预报误差进行高斯和模型拟合，从预报误差概率分布中学习得到轨道的异常检测门限，实现机动检测。Kelecy等[12]在历史TLE数据中选择一定长度的两相邻数据段分别进行多项式拟合，在数据段的外推时间中点计算预报值之差，通过数据段的移动构造偏差数据序列进行机动检测。Bai等[13]通过对轨道预报误差进行聚类分析，实现量级不同的机动检测。实际上，轨道预报误差拟合的思路也被应用于其他数据类型的轨道机动检测中。Xi等[14]提出了一种基于因果推理的方法，利用光学传感器得到的轨道数据，生成时间序列上的残差信息，对机动进行准确的检测。Kaderali等[15]提出一种改进的扩展卡尔曼滤波方法，并建立了轨道机动和观测信息之间的关联性指标。Singerman等[16]提出一种利用符号动力学的方法对雷达跟踪的目标进行机动检测。Aguilar Marsillach等[17]提出了一种新的基于马氏距离的二元假设检测方法对航天器异常行为进行表征。这些方法将轨道预报和移动窗口曲线拟合的数学方法结合，有效利用了轨道变化的物理规律，提高了机动检测的精度。但以上方法主要存在三点不足：一是没有对检测窗口大小的设置进行详细说明，依赖的是经验，而该参数对机动检测结果有显著的影响；二是机动时刻的估计模型不够精细；三是通过正向移动滑窗的方式无法识别末端数据段的机动。
本文在文献[11]的基础上提出自适应的轨道预报误差拟合法进行卫星的历史机动检测。首先，在进行机动检测前加入数据预处理并使用局部加权回归原理的平滑（locally weighted scatterplot smoothing，LOWESS）取代原有的三次样条曲线平滑提高TLE数据的信噪比；其次，提出机动检测窗口大小的自适应方法，采用交叉弧段预报进行机动时刻估计；再次，通过逆向移动滑窗方式，避免数据末端无法识别机动的问题，从而实现最新TLE和前序历史TLE数据间的机动检测；最后，通过仿真算例与文献[11]对比，检验本文改进方法对轨道机动检测的有效性。
1 自适应轨道机动检测模型
1.1 常规轨道预报误差拟合法
文献[11]是应用常规轨道预报误差拟合法对卫星历史机动进行检测的典型文献，流程如图1所示。该方法在获取目标轨道的编目值和预报值之前未进行数据预处理，采用三次样条曲线平滑在降噪的同时没有很好地保留目标轨道的详细特征，目标TLE轨道数据的信噪比还有待提升。
本文相对文献[11]所做出的改进主要体现在以下三方面。
首先，文献[11]在生成偏差数据序列时借鉴了移动窗口的思想，设置检测窗口，窗口向后移动以生成大量与预报时间相关联的偏差数据，但对不同的卫星进行机动检测时依赖经验配置检测窗口，难以先验获取最佳窗口，在对大量轨道特性不同的卫星进行机动检测时，固定的检测窗口往往达不到理想的检测效果。因此需要引入自适应检测窗口配置方法，使得对更多目标进行机动检测时能够自适应计算窗口宽度参数，从而达到较高的检测精度。
其次，文献[11]在进行轨道预报误差的获取和分布拟合时无法获得数据序列末端TLE数据点的误差分布，从而导致文献[11]检测模式下的轨道预报误差拟合法进行机动检测时无法检测到最近历元时刻航天器的机动情况，本文提出的逆向预报滑窗方法能有效解决这一问题，如图2所示。
图2逆向移动滑窗示意图
Fig.2Diagram of reverse window moving
再次，文献[11]对机动时刻的估计模型从数学统计的思想出发，以异常数据段中异常值最多的数据点作为空间目标机动的发生时刻，对机动时刻的估计不够精细，且无法得到关于机动量的任何信息。本文考虑采用动力学含义更完备的基于交叉弧段预报的机动时刻估计方法，对机动参数进行估计。
1.2 基于逆向移动滑窗的轨道预报误差拟合法
文献[11]中轨道预报误差拟合法采用正向移动滑窗的模式，该模式下的机动检测方法无法对最新TLE时刻中产生的机动进行检测，而最新测定轨数据的机动检测对近实时的目标异动告警具有重要意义。针对该问题，本节提出基于逆向移动滑窗的轨道预报误差拟合法，如图2所示。
以TLE数据序列的末端数据为预报起点，采用SGP4预报模型依次向TLE数据序列的起点进行逆向移动滑窗轨道预报，预报至窗口最后的数据点，而后预报起点变更为末端数据的前一个数据点，重复这一过程即可获得通过逆向移动滑窗得到的包含末端信息的轨道预报误差。该方法可以在检测模式上弥补采用正向移动滑窗的传统方法无法对最近历元时刻机动进行检测的缺陷。
1.3 轨道数据的预处理及平滑方法改进
进行机动检测时，空间目标的轨道误差由同一时刻轨道根数的实际值与预报值作差得到。轨道根数的实际值直接从TLE编目数据中得到，其预报值则通过与TLE配套的SGP4/SDP4轨道预报模型预报得到。与文献[11]相比，本文对平滑方法进行了改进。
由于测轨、定轨以及模型精度（如大气模型、地球模型等）等因素的影响，编目值和预报值会在其真值附近抖动，如图3所示。
图3平滑处理前数据特性示意图
Fig.3Diagram of data characteristics before smoothing
抖动现象可能导致偏差数据即使在正常情况下也会变大，从而增加虚警的概率。因此有必要对原始数据进行平滑处理，本文对基于三次样条插值曲线拟合的平滑方法和基于LOWESS方法进行了比较，并采用LOWESS方法取代原有的三次样条曲线平滑方法。
LOWESS方法是一种利用局部加权多项式拟合来平滑数据的非参数方法[18]，对有一定趋势或者季节性的数据有较好的拟合效果。使用该方法对特征参数序列进行平滑处理时，以待平滑的数据点为中心，对邻近一定长度的数据段作加权多项式回归，在时间序列上越靠近待平滑点的数据点在多项式拟合时权重越大。对于所有的数据点均可得到这样的回归线，每条回归线中心值的连线构成该特征参数序列数据的平滑曲线。
选用经典权函数作为多项式拟合时候的权值函数，即：
(1)
式中：ωj为邻近数据点的权值；ti表示待平滑数据点的历元时刻；tj表示邻近数据点的历元时刻且i-l≤j≤i+l，l为邻近数据点与中心数据点之间的数据个数；d（ti）表示邻近数据段中与ti间隔最长的时间距离。
1.4 机动检测窗口自适应改进
1.4.1 检测窗口定义及影响
以TLE数据作为原始数据进行卫星历史机动检测时，检测窗口定义为目标TLE序列中连续的数据点集合，如图4所示。检测窗口的大小即该集合中数据点的数量，虚线框中给出了大小为4的检测窗口，连续的检测窗口构成了窗口数据序列。
图4检测窗口示意图
Fig.4Detection window diagram
设目标TLE的数据点总数为n，检测窗口大小为w，采用轨道预报误差拟合法进行机动检测时，窗口数据序列的大小可表示为
(2)
将检测窗口内的TLE数据点依次向后预报得到预报值序列，每个预报值对应着预报时间，将预报值与预报值所在历元时刻的编目值作差，窗口数据序列内所有窗口执行此过程生成误差数据序列，误差数据序列中误差数据的个数可表示为
(3)
式中，S是关于检测窗口w的离散函数，一般而言，用于历史机动检测的TLE数据点n足够大，可近似作连续处理。误差数据序列存在极值，式（3）中对w求导，当w取（n+2）/2时，得到S的极值。
当检测窗口选择极值点时，误差数据最多，从统计角度看拟合效果最好。而实际应用时要综合考虑SGP4预报精度和检测效率。一颗轨道高度约为655 km的卫星经过SGP4模型预报7 d时，半长轴误差达到了0.42 km，且随着时间的增长而继续增大，预报30 d时达到1.43 km[19]。本文将检测窗口和检测对象的轨道特征及TLE数据特性结合，提出自适应检测窗口设置方法。
1.4.2 自适应检测窗口设计
对任意空间目标，TLE数据序列的时间跨度TD指其最后一个数据点所在历元时刻与第一个数据点所在历元时刻之差，则TLE频率f可定义为数据总数与时间跨度的比值。
(4)
本文对两个概念进行定义：当检测窗口设置为某个值时，机动检测的成功率达到峰值，即使增大检测窗口能满足同样的检测成功率，但效率更低，这样的情况称为“饱和检测”；满足“饱和检测”条件的检测窗口，称为最佳检测窗口。
对未失效的空间目标而言，其TLE频率与轨道高度之间存在某种联系，一般而言，轨道高度越高，TLE频率越低。LEO空间目标的TLE频率在1~5个/d之间，而GEO空间目标的TLE频率在1个/d左右。低轨道、中轨道和高轨道地球卫星三类空间目标的轨道高度与TLE频率关系如图5所示。
图5轨道高度与TLE频率的关系概貌图
Fig.5Overview of the relationship between the orbital altitude and TLE frequency
同时，从大量仿真分析与真实机动的对比来看，基于轨道预报误差拟合法机动检测的检测效果呈现以下特点：
1）对于相同的空间目标，不同时间段的TLE频率也不一样，对应的最佳检测窗口也有区别，例如表1中空间目标编号为28376的最佳检测窗口大小变化。
表1空间目标28376最佳检测窗口
Tab.1 Best detection windows of space object 28376
2）对于TLE频率接近的不同空间目标，即使轨道高度相差较大，对应的最佳检测窗口也不会有太大区别，例如表2中轨道高度相差较大的三个空间目标的最佳检测窗口大小变化。
表2TLE频率接近的不同轨道高度的空间目标最佳检测窗口
Tab.2 Best detection windows of space objects at various orbital altitude with almost identical TLE frequency
从检测效果的规律总结来看，检测窗口作为机动检测效果的最大影响参数，其大小与空间目标的TLE频率间存在关系。通过对10个不同空间目标的仿真分析，检测窗口大小与TLE频率之间的倍数关系近似呈现为图6所示关系。
图6检测窗口大小与TLE频率之间的倍数关系
Fig.6Multiple relationship between detection window size and TLE frequency
考虑到检测窗口一般为整数，则最佳检测窗口大小与TLE频率之间的关系以经验公式表示为
(5)
其中，〈·〉表示对括号内的数值四舍五入取整。
式（5）展示的经验公式以多项式的形式给出了最佳检测窗口大小的设置方式，能适应绝大部分空间目标的机动检测要求。本文仅以10个目标数据拟合式（5）的原因是可公开查询到的、已知机动数据的目标有限，在有更多标定数据的情况下，可基于本文建模思想获得更好的经验公式。
1.4.3 基于交叉弧段的机动估计
文献[11]以异常数据段中异常值数量最多的数据点对应时刻为机动时刻，该方法估计出来的机动时刻与实际的机动时刻理论上存在较大误差。本节在脉冲假设下，建立基于交叉弧段的机动估计模型进行机动时刻与机动量信息的估计。
设异常数据段中异常值数量最多的数据点对应历元时刻为tm，其相邻数据点历元时刻为t0和tf（t0＜tm＜tf），交叉弧段预报时间步数为d，则预报时间步长，通常可将时间步长选为10 s以内。将t0时刻的位置速度预报至tf时刻，每一步对应的预报值位置速度为Si；tf时刻的位置速度预报至t0时刻，每一步对应的预报值位置速度为S′i，如图7所示。
图7基于交叉弧段预报的机动估计示意图
Fig.7Diagram of maneuver estimation based on crossover arc prediction
式（6）定义了预报值之间的距离，通过位置矢量差的模表征。
(6)
假设距离集合dis中第j步对应的距离为最小值，则对应时刻可近似为机动发生的时刻，即
(7)
该时刻对应的速度矢量差即为空间目标在该时刻的机动量。
(8)
2 仿真校验
2.1 平滑方式改进效果
以空间目标33105为例，该空间目标信息如表3所示。
表3卫星轨道信息
Tab.3 Orbit information
截取其中2018年9月8日至12月17日的一段TLE编目数据，预报值通过对前一个TLE进行轨道预报得到，采用文献[11]方法中的三次样条平滑处理得到的结果如图8所示，与之对比的是LOWESS平滑处理后的数据结果，如图9所示。
图8三次样条平滑后的数据特性
Fig.8Data characteristics after cubic spline smoothing
从图8和图9经过平滑后的数据特性与图3未经处理的数据特性对比来看，LOWESS平滑方法在去除大量噪声的同时较好地保留了特征参数的变化性质，而三次样条平滑处理之后的参数序列误差明显较大。考虑到进行轨道机动检测时，过度拟合容易丢失机动引起的异常信息，增大漏检率，本文采用LOWESS方法进行数据的平滑处理，进一步减少野值的同时尽可能地保留含有轨道机动变化特征的信息。
图9LOWESS平滑后的数据特性
Fig.9Data characteristics after LOWESS smoothing
2.2 自适应机动检测仿真校验
选取表4中空间目标作为机动检测对象，值得注意的是，所选取目标均有真实机动数据可进行对比，但由于真实机动数据来源有限，大部分是处于LEO的空间目标，且时间跨度有所区别。其中卫星编号指卫星在美国态势感知网站[20]上的编目号。
表4机动检测对象轨道信息
Tab.4 Orbital information of detection objects
机动检测效果以检测成功个数、漏检个数、虚检个数作为检测效果的评判指标[21]。序列中真实存在的机动没有被检测出来属于漏检。当序列中某处没有发生机动，检测结果却显示有机动，则属于虚检。表4中空间目标的检测结果如表5所示。其中真实机动从部分有公布的数据[22-23]中得到。
表5机动检测结果
Tab.5 Results of maneuver detection
从检测结果来看，6个目标总体检测成功率达到了90.625%，漏检率为9.375%，虚检率几乎可以忽略不计，从检测效果方面验证了本文所提方法的有效性。表5中的漏检信息如表6所示。
表6漏检机动信息
Tab.6 Maneuver information of leak detection
表6中的漏检可能是航天器机动量级较小导致的。值得注意的是，使用者可以自行调整窗口大小观察检测结果的变化，对部分目标可以实现漏检率为0的效果，但是实际上大部分空间目标是没有机动先验信息的，使得调节参数较为困难，本文提供的自适应检测窗口给出参考值能够帮助解决这种困难。
以空间目标41240、23710以及37384为例，对检测效果进行进一步描述。
41240是美国于2016年1月17日发射的遥感卫星JASON3，检测效果如图10所示。根据国际激光测距服务组织发布的真实机动数据[22]，2017年1月1日至2019年1月1日，该卫星共执行6次轨道机动，本文检测结果与其能够较好吻合。
图10JASON3机动检测结果
Fig.10Maneuver detection results of JASON3
23710是加拿大于1995年11月4日发射的一颗遥感卫星，名称是RADARSAT-1。检测结果显示其在2010年共执行6次机动，如图11所示。
图11RADARSAT-1机动检测结果
Fig.11Maneuver detection results of RADARSAT-1
37384是中国于2011年4月9日发射的一颗北斗组网卫星——BeiDou C08。文献[23]显示，该卫星在2015年1月9日13:52:00进行了一次机动；本文检测结果显示其在2015年1月9日21:42:58进行了一次机动，如图12所示。
表7中给出了本文方法与文献[11]方法之间的机动参数估计结果对比。
图12BeiDou C08机动检测结果
Fig.12Maneuver detection results of BeiDou C08
表7机动参数估计结果对比
Tab.7 Comparison of maneuver estimation results
2.3 近实时机动检测仿真校验
以编号为27424的空间目标为例，27424是美国于2002年5月4日发射的一颗遥感卫星，名为AQUA。AUQA卫星在2019年9月11日左右发生过一次机动，选取该卫星TLE序列中2019年9月12日和9月11日的两相邻TLE，将其加入已有历史数据中，应用本文方法得到检测结果如图13所示。
图13本文方法对最近历元时刻中的机动检测结果
Fig.13Detection results of the maneuver in newest epoch time by the proposed method
由图13可知，本文方法在2019年9月11日 19:29:26检测出一次机动，有效识别出最近历元时刻的机动。文献[11]的检测方法无法对新产生的机动进行有效识别，仅能作为历史机动检测方法。原因是，本文提出的逆向移动滑窗方法解决了文献[11]方法无法获得最新历元时刻预报误差信息的问题，从而实现对最新历元时刻的机动检测。本文方法实现了近实时机动检测，大幅减少了对最新历元时刻中的机动进行识别和检测的响应时间。
