本发明涉及航天测量与控制领域，尤其涉及一种基于tle根数的在轨航天器变轨检测方法。
背景技术：
随着人类航天科技的不断进步与发展，人类太空活动的剧烈增加，在轨运行的航天器呈不断大幅增加的趋势。根据美国航天局公布的tle数据统计，目前在轨空间目标总数超过16000个，其中在轨航天器数量超过4000个，而在轨时间最长的vanguard1卫星，至今已超过60年。不同航天器由于其担负的任务不同、轨道类型不同等因素，需要不定期的进行轨道调整以保证完成其任务。对于空间目标管理方，对于航天器变轨的信息预先是未知的，因此需要及时从探测数据中发现航天器轨道机动行为，以及时作出轨道探测策略的调整。由于在轨空间目标数量巨大，每天对单个目标的观测圈次是非常有限的。如果在轨航天器进行了轨道调整，而对应的管理部门不能及时发现并匹配至该目标，则可能失去该目标的编目管理。
近年来，随着太空领域竞争的进一步加剧，太空安全在国家安全体系中的地位日益凸显，快速识别在轨航天器变轨事件，准确分析变轨量及变轨操作对整体太空态势的影响都至关重要，对保障国家太空安全具有非常重要的意义。
由于美国tle根数的具有更新快、误差小和易获得的特点，基于tle根数开展在轨航天器变轨检测分析已成为广大学者研究的热点。目前，有的学者提出了sacm方法，该方法从tle数据中的轨道半长轴变化出发，从样本数据中提取得到半长轴变化阈值，进而利用该阈值对探测数据进行判断，但该方法未能给出具体变轨时间及变轨量的计算方法，且判定阈值的选取方法存在缺陷；有的学者提出了基于tle数据的综合判定方法，但该方法对于一段时间内连续变轨的情况未能很好的检测识别。
鉴于上述原因，本发明针对在轨航天器空间轨道运行规律及轨道控制的特点，提出了一种基于tle根数进行在轨航天器变轨检测的方法，能够有效的解决航天器在轨运行过程中轨道发生变化的时间、轨道变化量以及多次变轨识别的问题。
技术实现要素：
本发明的目的就在于为了解决上述问题而提供一种基于tle根数的在轨航天器变轨检测方法。
本发明通过以下技术方案来实现上述目的：
本发明在给定的时间段内，在轨航天器发生变轨的时间及轨道半长轴变化量的计算方法，包括如下步骤：
步骤一：在给定的时间段[t0,t1](t1-t0≥30天)内，已知在轨航天器共有k条按时间升序排列的tle根数，则根据每条tle根数可计算得到航天器轨道半长轴
步骤二：根据步骤一中的半长轴序列ai，计算相邻两条tle根数的半长轴差，即δai＝ai 1-ai，其中i＝1,2,…k-1；
步骤三：根据步骤二中的δai(i＝1,2,…k-1)，计算其绝对值的平均值为
步骤四：对于低轨航天器(轨道高度小于2000千米)，取σa0＝5米；对于中高轨航天器(轨道高度大于2000千米)，取σa0＝100米；若根据步骤三中计算的
步骤五：重复步骤三至步骤四的流程，直至步骤四中
步骤六：将步骤二中的半长轴差δai(i＝1,2,…k-1)的绝对值|δai|与步骤五中的变轨筛选门限δa0逐一进行比较，若|δai|≥δa，则表示ai 1为航天器变轨后的半长轴，其对应的tle根数为变轨后的根数；若连续有p个半长轴ai 1,ai 2,…,ai p被判定为变轨后的半长轴且对应的半长轴差δai,δai 1,…,δai p-1符号相同，则ai对应航天器此次变轨前更新的最后一条tle根数，ai p对应航天器此次变轨结束后更新的第一条tle根数；
步骤七：将步骤六中ai和ai p对应的tle根数历元分别记为ts,te；在时间段[ts,te]内，基于sgp4/sdp4模型和航天器变轨前后对应的两条tle根数，计算其在j2000惯性坐标系中一分钟一点的位置分别为[x1qy1qz1q]，[x2qy2qz2q]，其中q＝1,2,…,[1440×(tn-ts)]，[1440×(tn-ts)]为(te-ts)×1440的整数部分，根据上述位置分量，则可计算相对距离
步骤八：根据步骤七中的相对距离dq序列，采用拉格朗日插值算法即可得到时间段[ts,te]内相对距离的最小值dqmin，此最小距离对应的时刻记为tmin，则tmin即为航天器的变轨时刻；基于步骤六中航天器变轨前后的两条tle根数，计算其在tmin时刻的平均半长轴，分别记为
本发明的有益效果在于：
本发明是一种基于tle根数的在轨航天器变轨检测方法，与现有技术相比，本发明：针对在轨航天器空间轨道运行规律及轨道控制的特点，提出了一种基于tle根数的航天器变轨时间及轨道半长轴变化量的计算方法，设置了动态变轨门限值，有效的解决航天器在轨运行过程中轨道发生变化的时间、轨道变化量以及多次变轨识别的问题。
附图说明
图1是2019年6月7日至2019年9月13日期间目标29499半长轴变化趋势图；
图2是2019年6月7日至2019年9月13日期间目标29499半长轴差的变化趋势图；
图3是2019年6月7日至2019年9月6日期间目标41469半长轴变化趋势图；
图4是2019年6月7日至2019年9月6日期间目标41469半长轴差的变化趋势图；
具体实施方式
下面结合附图对本发明作进一步说明：
本发明在给定的时间段内，在轨航天器发生变轨的时间及轨道半长轴变化量的计算方法，包括如下步骤：
步骤一：在给定的时间段[t0,t1](t1-t0≥30天)内，已知在轨航天器共有k条按时间升序排列的tle根数，则根据每条tle根数可计算得到航天器轨道半长轴
步骤二：根据步骤一中的半长轴序列ai，计算相邻两条tle根数的半长轴差，即δai＝ai 1-ai，其中i＝1,2,…k-1；
步骤三：根据步骤二中的δai(i＝1,2,…k-1)，计算其绝对值的平均值为
步骤四：对于低轨航天器(轨道高度小于2000千米)，取σa0＝5米；对于中高轨航天器(轨道高度大于2000千米)，取σa0＝100米；若根据步骤三中计算的
步骤五：重复步骤三至步骤四的流程，直至步骤四中
步骤六：将步骤二中的半长轴差δai(i＝1,2,…k-1)的绝对值|δai|与步骤五中的变轨筛选门限δa0逐一进行比较，若|δai|≥δa，则表示ai 1为航天器变轨后的半长轴，其对应的tle根数为变轨后的根数；若连续有p个半长轴ai 1,ai 2,…,ai p被判定为变轨后的半长轴且对应的半长轴差δai,δai 1,…,δai p-1符号相同，则ai对应航天器此次变轨前更新的最后一条tle根数，ai p对应航天器此次变轨结束后更新的第一条tle根数；
步骤七：将步骤六中ai和ai p对应的tle根数历元分别记为ts,te；在时间段[ts,te]内，基于sgp4/sdp4模型和航天器变轨前后对应的两条tle根数，计算其在j2000惯性坐标系中一分钟一点的位置分别为[x1qy1qz1q]，[x2qy2qz2q]，其中q＝1,2,…,[1440×(tn-ts)]，[1440×(tn-ts)]为(te-ts)×1440的整数部分，根据上述位置分量，则可计算相对距离
步骤八：根据步骤七中的相对距离dq序列，采用拉格朗日插值算法即可得到时间段[ts,te]内相对距离的最小值dqmin，此最小距离对应的时刻记为tmin，则tmin即为航天器的变轨时刻；基于步骤六中航天器变轨前后的两条tle根数，计算其在tmin时刻的平均半长轴，分别记为
1、近地航天器变轨检测
从美国公布的tle根数中选择norad编号为29499的航天器作为分析目标，该航天器是欧洲气象卫星开发组织在2006年发射的一颗气象卫星。从2019年6月7日至2019年9月13日公布的tle根数中提取目标29499的根数，其半长轴的变化趋势如图1所示：
从图1可知，在2019年6月7日至2019年9月13日期间，该目标有2次明显的变轨行为。
采用本发明的方法，目标29499半长轴差的绝对值|δa'j|及4次迭代门限示意图如图2所示：
图2中，筛选门限历次迭代过程中计算得到的
表1目标29499变轨筛选门限迭代过程结果
从表1可知，经过4次迭代后，
表2目标29499变轨分析结果
表1和表2的计算结果表明，采用本发明的方法可以准确的检测低轨航天器的变轨时间及轨道半长轴变化量。
2、中高轨航天器变轨检测
从美国公布的tle根数中选择norad编号为41469的航天器作为分析目标，该航天器是印度在2016年发射的一颗区域导航卫星。从2019年6月7日至2019年9月6日公布的tle根数中提取目标41469的根数，其半长轴的变化趋势如图3所示：
从图3可知，在2019年6月7日至2019年9月6日期间，该目标有3次明显的变轨行为。
采用本发明的方法，目标41469半长轴差的绝对值|δa'j|及2次迭代门限示意图如图4所示：
图4中，筛选门限历次迭代过程中计算得到的
表3目标41469变轨筛选门限迭代过程结果
从表3可知，经过2次迭代后，
表4目标41469变轨分析结果
表3和表4的计算结果表明，采用本发明的方法可以准确的检测高轨道航天器的变轨时间及轨道半长轴变化量。
以上显示和描述了本发明的基本原理和主要特征及本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。
1.一种基于tle根数的在轨航天器变轨检测方法，其特征在于：在给定的时间段内，在轨航天器发生变轨的时间及轨道半长轴变化量的计算方法，包括如下步骤：
步骤一：在给定的时间段[t0,t1](t1-t0≥30天)内，已知在轨航天器共有k条按时间升序排列的tle根数，则根据每条tle根数可计算得到航天器轨道半长轴
步骤二：根据步骤一中的半长轴序列ai，计算相邻两条tle根数的半长轴差，即δai＝ai 1-ai，其中i＝1,2,…k-1；
步骤三：根据步骤二中的δai(i＝1,2,…k-1)，计算其绝对值的平均值为
步骤四：对于低轨航天器(轨道高度小于2000千米)，取σa0＝5米；对于中高轨航天器(轨道高度大于2000千米)，取σa0＝100米；若根据步骤三中计算的
步骤五：重复步骤三至步骤四的流程，直至步骤四中
步骤六：将步骤二中的半长轴差δai(i＝1,2,…k-1)的绝对值|δai|与步骤五中的变轨筛选门限δa0逐一进行比较，若|δai|≥δa，则表示ai 1为航天器变轨后的半长轴，其对应的tle根数为变轨后的根数；若连续有p个半长轴ai 1,ai 2,…,ai p被判定为变轨后的半长轴且对应的半长轴差δai,δai 1,…,δai p-1符号相同，则ai对应航天器此次变轨前更新的最后一条tle根数，ai p对应航天器此次变轨结束后更新的第一条tle根数；
步骤七：将步骤六中ai和ai p对应的tle根数历元分别记为ts,te；在时间段[ts,te]内，基于sgp4/sdp4模型和航天器变轨前后对应的两条tle根数，计算其在j2000惯性坐标系中一分钟一点的位置分别为[x1qy1qz1q]，[x2qy2qz2q]，其中q＝1,2,…,[1440×(tn-ts)]，[1440×(tn-ts)]为(te-ts)×1440的整数部分，根据上述位置分量，则可计算相对距离
步骤八：根据步骤七中的相对距离dq序列，采用拉格朗日插值算法即可得到时间段[ts,te]内相对距离的最小值dqmin，此最小距离对应的时刻记为tmin，则tmin即为航天器的变轨时刻；基于步骤六中航天器变轨前后的两条tle根数，计算其在tmin时刻的平均半长轴，分别记为
