一种基于导航数据的应急航天器星历确定方法及系统与流程

1.本发明涉及航天航空领域，具体涉及一种基于导航数据的应急航天器星历确定方法及系统。


背景技术：

2.随着航天与空间技术的飞速发展，越来越多的航天器，包括卫星、载人飞船、空间站、深空探测器等，被发射到太空中，从而使得太空中的可用空间大大减小。据统计，截止2019年10月4日，地球上约有10779个直径超过10cm的空间物体在太空中，其中包括5181个航天器和14598个空间碎片。此外，空间中还存在着大量直径小于10cm的空间碎片。这些都对正常运行的航天器的安全构成了极大的威胁。因此，规避空间碎片已成为一项非常重要的应急航天任务。该任务的前提是确定航天器的星历，并预测其轨道，以确定它是否会与轨道附近的空间碎片相撞。如果有发生碰撞的危险，有必要及时进行轨道控制以避免危险。此外，当某些自然灾害发生时，例如火山爆发、地震和森林火灾等，无法在地面上观察到灾害现场。但是，诸如侦察卫星和遥感卫星之类的航天器却非常适合在这种应急空间侦察任务中对灾害现场进行空间观测，但前提也是要先确定其星历，以便计算这些卫星何时可以运行到灾害发生地方的上空从而进行侦察和探测。除了这些应急空间任务外，航天器还可能发生其他一些紧急情况。例如，导航设备无法提供足够的导航数据、数据包含较大的测量噪声或测量时间间隔较大等，进而导致导航数据不可用，难以解算确定航天器的星历。
3.关于应急航天器的星历确定，美国空军太空司令部(afspc)发布的两行轨道根数(two line elements，tle)可能是最准确的描述。tle不是瞬时轨道根数，而是通过特殊方式消除周期性变化得到的平均轨道根数。另一方面，tle具有标准化的格式，并且只能与sgp(simplified general perturbations)轨道模型配套使用。此外，应急航天任务还强调自主与实时性，以满足时效性的要求。一般地，传统的轨道确定方法将星上导航数据传输到地面站进行计算，然后再将结果传输到航天器。尽管其结果具有较高的精确度，但它占用大量计算资源，并且消耗大量的时间进行数据传输。因此，当考虑到实时性和星载计算能力时，传统的轨道确定方法对于应急航天任务是不可接受的。
4.实际上，大多数轨道模型，例如二体模型、j2摄动模型和高精度轨道预测(highprecision orbit prediction，hpop)模型，都使用瞬时轨道根数进行轨道预报。但是， tle的参数是“平均”轨道根数，与瞬时轨道根数不同，尽管这两种轨道根数非常相似，但是如果将tle替换为瞬时轨道根数代入sgp模型进行轨道预报，将会造成很大的预报误差。因此，将瞬时轨道根数转换为tle是非常必要的。但是，afspc并没有公布这两种轨道根数相互转换的方法。为了解决这个问题，必须建立模型解算可以用于sgp的两行轨道根数。
5.ken cranford在1970年建立了sgp4模型，并用于近地卫星轨道预报。将lane 和cranford的解析理论进行简化，将重力模型和空气动力学模型作为brouwer和功率密度函数的解，便可以得到sgp4模型。该模型考虑了几个长期和长周期扰动项，以确保轨道预报的精度。例如，j2，j3和j4项以及和空气阻力等。但是由于未考虑短周期扰动项，sgp4模型的位
置预报精度只能达到2km。vallado根据tle和sgp4 模型开发了一个轨道预报软件包。胡敏和李卫平等，提出了瞬时轨道根数和平均轨道根数之间的一种转换方法。为了提高太空飞行的安全性，韩蕾、史纪鑫等人，建立了基于sgp4和tle的空间碎片轨道预测方法。但是，几乎所有这些研究都未考虑航天器的自主性和实时性。李梦奇提出了两种方法来修改tle的阻力系数，并将 sgp4的预测精度提高了20％以上。陈俊宇等，对低地球轨道(low earth orbit，leo) 碎片的阻力系数进行研究，并估计了2000多个空间碎片的阻力系数来检验tle的精度。wesam等，运用批量最小二乘技术和测距率来校正历元时刻的tle角度误差，提高了圆形轨道确定的准确性。由于这些研究仅研究了星历确定中tle的部分元素，因此需要对其他元素的解算进行更多的研究。此外，这些方法也可能不适用于应急航天任务。因此，建立一种高效自主的航天器星历确定方法来解算应急航天器的星历则十分重要。


技术实现要素：

6.本发明的目的在于克服上述技术缺陷，提出了两种基于卫星导航数据的方法来确定航天器星历。当导航数据有限时，则可利用迭代方法解算，因为仅依靠某一时刻的航天器飞行状态矢量即可执行迭代方法。如果导航数据足够多时，则可以利用滤波方法解算航天器星历。两种方法都具有很高的效率，并且可以满足tle格式的精度要求。
7.为实现上述目的，本发明的实施例1提出了一种基于导航数据的应急航天器星历确定方法，所述方法包括：
8.从航天器的导航设备获取导航数据，所述导航数据包括：航天器飞行状态矢量；
9.当导航数据不足，则基于导航数据使用航天器星历迭代解算方法计算得到航天器的星历；否则，基于导航数据使用ukf方法计算得到航天器的星历；所述导航数据不足是指基于导航数据使用ukf方法无法收敛到所设定的精度。
10.作为上述方法的一种改进，基于导航数据使用航天器星历迭代解算方法计算得到航天器的星历，具体包括：
11.步骤s1)从导航设备获取航天器当前飞行状态矢量；将其转换到temed坐标系得到该坐标系下的飞行状态矢量ρ0＝[r
0 v0]
t
，r0为位置矢量，v0为速度矢量；迭代次数的初始值k＝1；
[0012]
步骤s2)根据轨道动力学理论，基于ρ
k-1
计算航天器在第k次迭代下的瞬时轨道根数ik为轨道倾角，ωk为升交点赤经，ek为偏心率，ωk为近地点角距，mk为平近点角，nk为平均运动角速度，为大气阻力系数，为常数；
[0013]
步骤s3)将计算得到的航天器瞬时轨道根数视为tle并代入sgp4模型进行轨道预报，得到当前迭代步数下的航天器飞行状态矢量ρk＝[r
k vk]
t
；
[0014]
步骤s4)利用误差参数δρ＝ρ
k-ρ0修正航天器的飞行状态，得到ρ
k 1
＝ρk δρ；
[0015]
步骤s5)若误差参数|δρ|小于给定值；则αk为所求的平均轨道根数α0，即航天器的星历；否则，k加1后转入步骤s2)。
[0016]
作为上述方法的一种改进，所述基于导航数据使用ukf方法计算得到航天器的星历，具体包括：
[0017]
步骤t1)获取星历状态初值x0和状态协方差初值p0，令p
0,0
＝p0，滤波递推
次数的初始值k＝1；
[0018]
步骤t2)获得一组sigma点集和相应的权值：
[0019][0020][0021]
λ＝α2(l κ)-l
[0022]
其中，ω(i)为权重，λ是缩放比例参数；l代表状态矢量的维数；α代表主要比例因子，该因子可控制选取的sigma点的分布状态，其取值范围是10-3
≤α≤1；β代表次要比例因子，β≥0；κ代表三次比例因子，其值为零；
[0023]
步骤t3)计算上述sigma点集的一步预测值
[0024][0025]
步骤t4)通过对sigma点集的一步预测值加权求和得到状态量的一步预测及协方差矩阵p
k 1,k
：
[0026][0027][0028]
其中，q为过程噪声协方差阵；
[0029]
步骤t5)根据得到的一步预测值，再次根据ut变换，得到新的sigma点集
[0030][0031]
步骤t6)将新的sigma点集代入离散测量方程，得到每个点对应的预测的测量量
[0032][0033]
其中，表示预测位置矢量，v(k)表示测量噪声矢量；
[0034]
步骤t7)通过对步骤t6)的sigma点集加权求和得到预测的测量量
[0035]
[0036][0037][0038]
其中，为测量量的方差，为状态和测量量的协方差；r为测量噪声协方差阵；
[0039]
步骤t8)计算得到滤波增益矩阵k
k 1
：
[0040][0041]
步骤t9)更新状态矢量和协方差矩阵p
k 1,k 1
：
[0042][0043][0044]
其中，z
k 1
为第k 1的真实测量值；
[0045]
步骤t10)判断是否小于阈值，若判断结果为肯定的，则为所求的平均轨道根数α0，即航天器的星历；否则，k加1后转入步骤t2)。
[0046]
作为上述方法的一种改进，所述获取星历状态初值x0和状态协方差初值p0，具体包括：
[0047]
设航天器的当前飞行状态矢量为ρ0＝[r
0 v0]
t
，r0＝[r
01 r
02 r
03
]为位置矢量， r
01
,r
02
,r
03
分别为三个方向的位置坐标值；v0＝[v
01 v
02 v
03
]为速度矢量，v
01
,v
02
,v
03
分别为三个方向的速度值；则动量矩矢量h为：
[0048]
h＝[h
1 h
2 h3]
t
＝r0×
v0[0049]
h1,h2,h3为动量矩矢量h三个方向的标量；
[0050]
动量矩大小h为：
[0051][0052]
则轨道倾角i为：
[0053][0054]
偏心率矢量e为：
[0055][0056]
其中，μ＝3.986004356
×
10
14
为地球引力常数；e1,e2,e3为偏心率矢量e三个方向的标量；
[0057]
则偏心率大小e为：
[0058][0059]
由此计算得到轨道半长轴a：
[0060][0061]
轨道平均运动角速度n为：
[0062][0063]
升交点矢量n为：
[0064]
n＝[n
1 n
2 n3]
t
＝[0 0 1]
t
×h[0065]
其中，n1,n2,n3为升交点矢量n三个方向的标量；
[0066]
升交点矢量的大小n为：
[0067][0068]
考虑到升交点赤经ω的范围为0至360度，所以分两种情况计算，当n2≥0时，
[0069][0070]
当n2＜0时，
[0071][0072]
近地点角距ω分两种情况，当e3≥0时，
[0073][0074]
当e3＜0时，
[0075][0076]
当r
03
≥0时，升交点角距u为：
[0077][0078]
其中，距离r0为
[0079]
当r
03
＜0时，升交点角距u为：
[0080][0081]
则真近点角f为：
[0082]
f＝u-ω
[0083]
当f＜180
°
时，偏近点角e为：
[0084]
[0085]
当f≥180
°
时，偏近点角e为：
[0086][0087]
则平近点角m为：
[0088]
m＝e-esin(e)
[0089]
则星历状态初值x0为：
[0090]
x0＝[i ω e ω m n b
*
]
t
[0091]b*
为大气阻力系数，根据经验值给定；
[0092]
状态协方差初值p0为：
[0093]
p0＝diag([10-7 10-7 10-7 10-7 10-7 10-7 10-7
])。
[0094]
本发明的实施例2提供了一种基于导航数据的应急航天器星历确定系统，所述系统包括：数据获取模块和航天器星历计算模块；
[0095]
数据获取模块，用于从航天器的导航设备获取导航数据，所述导航数据包括：航天器飞行状态矢量；
[0096]
航天器星历计算模块，用于当导航数据不足，则基于导航数据使用航天器星历迭代解算方法计算得到航天器的星历；否则，基于导航数据使用ukf方法计算得到航天器的星历；所述导航数据不足是指基于导航数据使用ukf方法无法收敛到所设定的精度。
[0097]
本发明的优势在于：
[0098]
1、本发明的方法提出了快速有效的迭代和滤波解算方法，解算精度达到了两行轨道根数标准格式精度要求；
[0099]
2、本发明的方法可操作性强，步骤简单，且针对观测数据不足和充足两种情况均可解算航天器星历；
[0100]
3、本发明的方法适用于应急航天器星历解算场景，占用星上计算资源少，自主性和实时性强。
附图说明
[0101]
图1为本发明的迭代方法的流程图；
[0102]
图2为本发明的ukf方法的流程图；
[0103]
图3(a)为轨道倾角和迭代次数的示意图；
[0104]
图3(b)为偏心率和迭代次数的示意图；
[0105]
图4(a)为迭代算法(没有测量噪声)解算的星历进行轨道预报的位置误差的示意图；
[0106]
图4(b)为迭代算法(没有测量噪声)解算的星历进行轨道预报的速度误差的示意图；
[0107]
图5为测量噪声对迭代方法的影响的示意图；
[0108]
图6(a)为不同测量噪声下的迭代方法的轨道预报的位置误差示意图；
[0109]
图6(b)为不同测量噪声下的迭代方法的轨道预报的速度误差示意图；
[0110]
图7为ekf和ukf结果的精度对比示意图；
[0111]
图8(a)为ekf和ukf结果的轨道倾角滤波过程对比图；
[0112]
图8(b)为ekf和ukf结果的偏心率滤波过程对比图；
[0113]
图8(c)为ekf和ukf结果的平均运动角速度滤波过程对比图；
[0114]
图8(d)为ekf和ukf结果的滤波残差滤波过程对比图；
[0115]
图9(a)为ekf和ukf结果进行轨道预报的位置误差示意图；
[0116]
图9(b)为ekf和ukf结果进行轨道预报的速度误差示意图；
[0117]
图10(a)为3σ测量噪声对滤波算法的影响示意图；
[0118]
图10(b)为10σ测量噪声对滤波算法的影响示意图；
[0119]
图11(a)为3σ测量噪声下ekf和ukf进行轨道预报的位置误差示意图；
[0120]
图11(b)为3σ测量噪声下ekf和ukf进行轨道预报的速度误差示意图；
[0121]
图11(c)为10σ测量噪声下ekf和ukf进行轨道预报的位置误差示意图；
[0122]
图11(d)为10σ测量噪声下ekf和ukf进行轨道预报的速度误差示意图。
具体实施方式
[0123]
下面结合附图对本发明的技术方案进行详细说明。
[0124]
本发明的实施例1提出了一种基于导航数据的应急航天器星历确定方法，包括：
[0125]
1建模与坐标转换
[0126]
sgp4模型的变量是两行轨道根数。两行轨道根数的元素为sgp类型的平均运动角速度n0，轨道倾角i0，偏心率e0，平近点角m0，近地点角距ω0，升交点赤经ω0，平均运动角速度的一阶变化率平均运动角速度的二阶变化率和sgp类型的大气阻力系数所有变量都是通过特定方式消除周期性变化而获得的“平均”值。另外，它们是被定义在航天器进入其轨道的历元时刻。将tle代入为sgp4模型可以预报飞行状态，包括航天器的位置和速度，可以描述如下：
[0127][0128]
其中s(t)、r(t)和v(t)都是行向量，s(
·
)代表sgp4模型，α0代表tle的元素，t代表从历元时刻t0开始的时间间隔：
[0129][0130]
已知sgp4模型和tle的参考坐标系为temed(true equator and mean equinoxof date)参考系统。但是，导航数据的参考坐标系是wgs84(world geodetic system1984)。因此，需要将导航数据经过两次坐标转换，从wgs84转换为j2000.0，再从j2000.0转换为temed。
[0131]
2航天器星历迭代解算方法
[0132]
导航数据是迭代方法和滤波方法解算航天器星历的输入。两种方法是基于导航数据不足和充足两种情况建立的。在正常情况下，导航数据是与时间匹配的连续数据集。当导航数据充足时，滤波方法和迭代方法都可以解算航天器星历。但是，当发生特殊情况(例如导航设备故障)时，导航数据可能不可用。因此，在这种情况下只有有限数量的导航数据可以利用，滤波方法可能会由于测量数据不足而无法收敛，导致无法解算出满足精度的航天器星历。相反，迭代方法却可以很好地工作，因为一组有效的导航数据即可执行迭代过程。
因此，迭代方法在两种情况下均能有效地解算航天器星历，而滤波方法只有在导航数据充足时才能有效地解算航天器星历。
[0133]
根据轨道动力学理论，利用包含航天器位置和速度的gps导航数据，可以计算出航天器的瞬时轨道根数。但是由于消除周期性变化的特定过程未知，使得确定航天器的平均轨道根数也即星历变得非常困难。为了解决这个问题，将该过程视为“黑匣子”，利用“输出”来求解“输入”。航天器的星历为“输入”，航天器的飞行状态为“输出”，基于“黑匣子”的迭代方法过程如图1所示。
[0134]
首先，从导航设备获取航天器当前飞行状态矢量然后，将其转换到 temed坐标系得到ρ0＝[r
0 v0]
t
。接着，根据轨道动力学理论计算航天器在当前迭代步数下的瞬时轨道根数由于迭代方法无法解算大气阻力系数，因此根据经验值将设置为常数。然后，将计算得到的航天器瞬时轨道根数视为tle并代入sgp4模型进行轨道预报，得到当前迭代步数下的航天器飞行状态矢量ρk＝[r
k vk]
t
,(k＝1,2,3...)。接着，利用误差参数δρ＝ρ
k-ρ0修正航天器的飞行状态，得到ρ
k 1
＝ρk δρ。重复上述过程直到误差参数|δρ|小于给定值。则最后一次迭代的结果即为所求的平均轨道根数，也即航天器的星历。
[0135]
3航天器星历滤波解算方法
[0136]
3.1ekf解算方法
[0137]
通过迭代方法，可以根据航天器任意时刻的飞行状态矢量解算得到星历。但是，迭代方法的缺点也不容忽视。首先，它不能充分利用所有导航数据。其次，由于在整个迭代过程中未更新用来对比的标准飞行状态矢量，因此无法消除测量噪声的影响。另外，它不能解算大气阻力系数。最后，一旦无法获得导航数据中的飞行状态矢量，其功能将无法实现。为了充分利用所有导航数据并消除测量噪声的影响，建立了扩展卡尔曼滤波方法来解算航天器星历。实际上，它是基于状态估计理论建立的，其中星历是状态变量，而飞行状态是测量变量。
[0138]
类似地，假设航天器的星历为x＝[i ω e ω m n b
*
]
t
，位置测量数据为根据定义可知，星历在初始历元时刻是常值。因此，没有状态更新，也即状态转移矩阵是单位矩阵。测量方程如下：
[0139]
y＝hx v
[0140]
其中y＝[x y z]
t
是测量变量，h是测量矩阵，v是服从高斯噪声理论的测量噪声。
[0141]
由于sgp4模型的复杂性与非线性，因此无法使用来计算测量矩阵。但是可通过数值差商来计算测量矩阵h。首先，使用sgp4模型计算t
k-1
时刻的 y
k-1
＝[x
k-1 y
k-1 z
k-1
]
t
。然后，分别给星历的元素增加一个小量。在这里，以轨道倾角为例：
[0142]
x
k-1/i
＝[i δi ω e ω m n b
*
]
t
[0143]
然后利用sgp4模型预报航天器当前时刻的飞行状态y
k-1/i
＝[x
k-1/i y
k-1/i z
k-1/i
]
t
。进而可得到如下的数值差商：
[0144][0145]
其他六个元素的差商通过相同的方式获得，则测量矩阵为：
[0146][0147]
最后，即可建立完整的ekf(extended kalman filter)滤波方程。
[0148]
3.2ukf解算方法
[0149]
尽管ekf方法可以有效消除测量噪声并充分利用导航数据，但仍存在一些不足。首先，它对非线性测量方程进行了线性化，导致存在线性化误差，这可能使滤波发散。此外，用数值差商代替偏导数会增加计算的复杂度，占用过多计算资源。因此，建立无迹卡尔曼滤波方法来弥补扩展卡尔曼滤波方法的不足。
[0150]
ukf(unscented kalman filter)最特殊优势的是利用无迹变换(unscentedtransform，ut)而不是线性化来处理均值和协方差的传递。对于ut，为了获取原始状态矢量的均值和协方差，选择了一组选定的sigma点。然后，通过非线性模型进行映射。接着，计算sigma点的均值和协方差。这样，对于任意的非线性函数，均值和协方差可以达到二阶精度。
[0151]
类似地，状态变量由x表示，状态转换矩阵也为单位矩阵。测量方程为：
[0152][0153]
其中，是状态矢量。
[0154]
由于导航数据是离散的，并且仅使用了位置矢量，因此可以将上述公式改为：
[0155][0156]
其中，s(x,k)r表示预测位置矢量，v表示测量噪声矢量。
[0157]
sigma点及其权重通过以下方程计算：
[0158][0159][0160]
其中，是协方差矩阵p的平方根的第j列；
[0161]
λ＝α2(l κ)-l是缩放比例参数；l代表状态矢量的维数；α代表主要比例因子，该因子可控制选取的sigma点的分布状态，其经验取值范围是10-3
≤α≤1；β代表次要比例因子，通常β≥0；κ代表三次比例因子，其值通常为零。
[0162]
显然，ukf没有线性化。本质上，它实现了无迹变换以使sigma点的均值和协方差与原始特征相匹配。然后，利用sigma点的非线性映射得到状态量的概率密度函数。
[0163]
具体地，ukf方法解算星历的步骤包括：
[0164]
步骤t1)获取星历状态初值x0和状态协方差初值p0，令p
0,0
＝p0，迭代次数的初始值k＝1；
[0165]
所述获取星历状态初值x0和状态协方差初值p0，具体包括：
[0166]
设航天器的当前飞行状态矢量为ρ0＝[r
0 v0]
t
，r0＝[r
01 r
02 r
03
]为位置矢量， r
01
,r
02
,r
03
分别为三个方向的位置坐标值；v0＝[v
01 v
02 v
03
]为速度矢量，v
01
,v
02
,v
03
分别为三个方向的速度值；则动量矩矢量h为：
[0167]
h＝[h
1 h
2 h3]
t
＝r0×
v0[0168]
h1,h2,h3为动量矩矢量h三个方向的标量；
[0169]
动量矩大小h为：
[0170][0171]
则轨道倾角i为：
[0172][0173]
偏心率矢量e为：
[0174][0175]
其中，μ＝3.986004356
×
10
14
为地球引力常数；e1,e2,e3为偏心率矢量e三个方向的标量；
[0176]
则偏心率大小e为：
[0177][0178]
由此计算得到轨道半长轴a：
[0179][0180]
轨道平均运动角速度n为：
[0181][0182]
升交点矢量n为：
[0183]
n＝[n
1 n
2 n3]
t
＝[0 0 1]
t
×h[0184]
其中，n1,n2,n3为升交点矢量n三个方向的标量；
[0185]
升交点矢量的大小n为：
[0186][0187]
考虑到升交点赤经ω的范围为0至360度，所以分两种情况计算，当n2≥0时，
[0188]
[0189]
当n2＜0时，
[0190][0191]
近地点角距ω分两种情况，当e3≥0时，
[0192][0193]
当e3＜0时，
[0194][0195]
当r
03
≥0时，升交点角距u为：
[0196][0197]
其中，距离r0为
[0198]
当r
03
＜0时，升交点角距u为：
[0199][0200]
则真近点角f为：
[0201]
f＝u-ω
[0202]
当f＜180
°
时，偏近点角e为：
[0203][0204]
当f≥180
°
时，偏近点角e为：
[0205][0206]
则平近点角m为：
[0207]
m＝e-esin(e)
[0208]
则星历状态初值x0为：
[0209]
x0＝[i ω e ω m n b
*
]
t
[0210]b*
为大气阻力系数，根据经验值给定；
[0211]
状态协方差初值p0为：
[0212]
p0＝diag([10-7 10-7 10-7 10-7 10-7 10-7 10-7
])。
[0213]
步骤t2)获得一组sigma点集和相应的权值：
[0214]
[0215][0216]
λ＝α2(l κ)-l
[0217]
其中，ω(i)为权重，λ是缩放比例参数；l代表状态矢量的维数；α代表主要比例因子，该因子可控制选取的sigma点的分布状态，其取值范围是10-3
≤α≤1；β代表次要比例因子，β≥0；κ代表三次比例因子，其值为零；
[0218]
步骤t3)计算上述sigma点集的一步预测值
[0219][0220]
步骤t4)通过对sigma点集的一步预测值加权求和得到状态量的一步预测及协方差矩阵p
k 1,k
：
[0221][0222][0223]
其中，q为过程噪声协方差阵；
[0224]
步骤t5)根据得到的一步预测值，再次根据ut变换，得到新的sigma点集
[0225][0226]
步骤t6)将新的sigma点集代入离散测量方程，得到每个点对应的预测的测量量
[0227][0228]
其中，表示预测位置矢量，v(k)表示测量噪声矢量；
[0229]
步骤t7)通过对步骤t6)的sigma点集加权求和得到预测的测量量
[0230][0231][0232][0233]
其中，为测量量的方差，为状态和测量量的协方差；r为测量噪声协方差阵；
[0234]
步骤t8)计算得到滤波增益矩阵k
k 1
：
[0235][0236]
步骤t9)更新状态矢量和协方差矩阵p
k 1,k 1
：
[0237][0238][0239]
其中，z
k 1
为第k 1的真实测量值；
[0240]
步骤t10)判断是否小于阈值，若判断结果为肯定的，则为所求的平均轨道根数α0，即航天器的星历；否则，k加1后转入步骤t2)。
[0241]
4数值仿真
[0242]
方法和模型已经建立，下面将通过数值仿真来验证有效性和准确性。航天器的参数是从stk(systems tool kit)的tle数据库中选择的，而导航数据是由stk仿真模拟得到的。选择tle-00005航天器来实施仿真，它是stk软件包的tle数据库中的一个航天器，是美国海军在1958年3月17日发射的，其名称为vanguard1。它的平均运动为10.839419606圈/天，其偏心率为0.1852049。经过换算可知，它的周期约为133分钟，其轨道为椭圆轨道。它的这些特点与一些可执行应急航天任务的 leo航天器相似。此外，当周期小于255min时，sgp4的性能更好。因此，选择了该航天器，其两行轨道根数如表1所示：
[0243]
表1：tle-00005航天器两行轨道根数
[0244][0245]
4.1迭代方法仿真
[0246]
利用stk仿真，得到航天器在wgs84系中的飞行状态。任意选择一组状态量，文中选取2007年5月18日12：00：19.609时刻的位置和速度矢量。表2列出了不含测量噪声的该航天器的位置和速度矢量。对位置和速度矢量进行坐标转换，得到 temed系下的飞行状态矢量，代入迭代方法进行仿真。结果如表3所示，以轨道倾角和偏心率的绝对误差来说明迭代过程，如图3(a)和图3(b)所示。
[0247]
表2：tle-00005航天器的飞行状态参数(wgs84)
[0248][0249]
表3：迭代方法结果与stk值对比
[0250][0251]
根据表3中的比较，迭代方法具有很高的精度。当没有测量噪声时，四个角度的精度为至少10-6
度，偏心精度为10-10
，平均运动角速度精度为10-8
圈/天。根据图2，大约3次迭代即可收敛，9次迭代可达到给定的精度要求。利用由迭代方法解算的tle预报该卫星的位置和速度，并将其与stk仿真值进行比较，结果如图4(a) 和图4(b)所示。
[0252]
在图4(a)和图4(b)中，我们可以看到一天内，三个方向上的位置预报误差都在1m以内，而速度误差都在10-3
m/s以内。这表明在没有测量噪声的情况下，迭代方法是有效的，且解算精度很高。
[0253]
如前所述，迭代方法容易受到测量噪声的影响。为了说明这种影响，在导航数据中叠加服从高斯噪声理论的随机误差。噪声的均值设置为零，而噪声的方差分别设置为10m，20m和30m。解算得到的tle绝对误差如图5所示。根据图5，绝对误差随测量噪声增大而增大。然后，利用tle预报航天器的位置和速度，并与stk 值进行比较。结果如图6(a)和图6(b)所示。
[0254]
当测量噪声增大时，位置和速度预报误差也会增大。如果测量噪声的方差大于 10m，则位置和速度预报误差将分别达到1000m和1m/s以上。由此可见，迭代方法确实容易受到测量噪声的影响。还可以看到位置和速度预报误差在400分钟左右最小，这是因为仿真使用的是2007年5月18日12：00：19.609时刻(仿真时刻)的数据来解算星历。根据图2中的迭代方法流程，当预报值与仿真时刻的值之间的误差小于给定参数时，迭代终止并输出结果。因此，通过迭代方法解算的星历与仿真时刻密切相关，且预报误差在仿真时刻时最小。该航天器tle的初始时刻是 07138.20994918，也即2007年5月18日05：02：19.609，因此，初始时刻和仿真时刻的时间间隔约400分钟(418分钟)。
[0255]
4.2滤波方法仿真
[0256]
滤波参数的设置决定了滤波器的性能，因此，首先要仔细设置滤波器的参数。计算ekf测量矩阵的小量设置如下：
[0257]
[δi δω δe δω δm δn δb
*
]
t
＝[10-7 10-7 10-7 10-7 10-7 10-7 10-7
]
t
[0258]
测量噪声参数设置如下：
[0259]
r＝diag(10,10,10)
[0260]
通过设置过程噪声矩阵来调节滤波过程，避免滤波器发散。同样，我们根据这个原则设置ekf的过程噪声。以一天为预报时长，分别利用stk软件和本文的程序进行轨道预报，将stk预报值作为标准值，通过叠加小量到tle上，使本文的程序预报值与标准值之间的位置误差约为1m(σ10)，基于此，可以估计过程噪声矩阵。初始状态协方差的设置要避免滤波器发散，通过参考过程噪声来设置它。进而，给出过程噪声和初始状态协方差矩阵的值如下：
[0261]
q＝diag(10-7
,10-7
,10-7
,10-7
,10-7
,10-10
,10-7
)
[0262]
p0＝diag(10-7
,10-7
,10-7
,10-7
,10-7
,10-7
,10-7
)
[0263]
ukf的参数设置如下公式所示：
[0264]
[n α κ β]＝[7 0.05 0 1]
[0265]
特别地，ekf和ukf方法的过程噪声、测量噪声和收敛条件均相同。测量数据也由stk仿真得到，且测量时间间隔设置为1分钟。
[0266]
根据无测量噪声的导航数据，利用二体轨道预报模型计算滤波的初值，也即航天器星历的初始状态值，如下：
[0267]
x0＝[34.2502 178.4066 0.1855494 103.8038 277.4321 10.830012742
ꢀ-
1.1
×
10-4
]
t
[0268]
然后进行仿真，结果如表4和图7所示。选择轨道倾角、偏心率、平均运动角速度和均方根误差(root mean square error， rmse)的绝对误差来显示ekf和ukf的滤波过程，如图8(a)和图8(b)所示。
[0269]
表4：ekf和ukf方法结果
[0270][0271]
从表4、图8(a)和图8(b)可以知道，ekf和ukf都能有效解算星历。当测量噪声为10m时，ekf结果中，四个角度的精度至少为10-5
度，偏心率精度为10-9
，平均运动角速度精度为10-9
rev/day，大气阻力系数精度为10-7
；ukf结果中，四个角度的精度至少为10-5
度，偏心率精度为10-7
，平均运动角速度精度为10-10
rev/day，大气阻力系数精度为10-6
。
[0272]
此外，使用滤波方法解算的结果进行轨道预报并与stk值进行比较，结果如图 9(a)和图9(b)所示，位置误差均在5m以内，速度误差均在0.005m/s以内。可以推断，当测量噪声较小时，ekf和ukf的性能都很好，但ekf略优于ukf。
[0273]
为了说明大噪声对滤波方法的影响，下面叠加噪声到导航数据中然后再进行仿真。类似地，噪声的均值设置为零，而方差分别设置为3σ(σ＝10)和10σ。解算结果的绝对误差如图10(a)和图10(b)所示。将解算结果代入sgp4模型进行轨道预报并与stk值相比较，得到的位置和速度预报误差如图11(a)、图11(b)、图11 (c)和图11(d)所示。
[0274]
可以知道，tle的绝对误差以及位置和速度预报误差均随着噪声方差的增大而增大。但是在两个仿真算例中，位置预报误差均小于测量噪声的方差，且在大噪声条件下，ukf的性能优于ekf。
[0275]
对于其他特殊情况，如较大的初始误差，同样进行了数值仿真。通过之前仿真，可以知道初始状态值是由二体模型在导航数据不包含测量噪声条件下计算得到的。在此，对用来计算初始状态的导航数据叠加噪声，然后再利用二体模型计算初始状态，从而进行较大初始误差情况的数值仿真。
[0276]
测量噪声的均值设置为零，方差分别设置为10σ(σ＝10)和100σ，其他参数设置与
正常情况一致，仿真结果如表5和表6所示。
[0277]
表5：大初始误差条件下的ekf和ukf结果(10σ)
[0278][0279]
表6：大初始误差条件下的ekf和ukf结果(100σ)
[0280][0281]
从表5和表6中可以看出，在较大初始误差条件下，ekf和ukf仍可以解算得到满足tle格式精度的航天器星历。
[0282]
本文研究了基于导航数据的航天器星历自主确定方法。通过建立迭代和滤波方法，实现了单点和多点导航数据情况下的航天器星历解算。对于迭代方法，在没有测量噪声的情况下，3次迭代后开始收敛，9次迭代则能收敛到给定精度。使用迭代方法解算的tle进行轨道预报，在一天时间内，三个方向上的位置预报误差均在1m 之内，速度预报误差均在10-3
m/s之内。但是，迭代方法容易受到测量噪声的影响，且测量噪声越大，绝对误差越大。对于滤波方法，在较小测量噪声条件下，ekf的性能要优于ukf。当测量噪声方差为10m时，ekf和ukf的解算结果，四个角度的精度为10-5
度，偏心度精度为10-7
，平均运动角速度精度为10-9
rev/day，大气阻力系数精度为10-6
。利用ekf和ukf解算的tle进行轨道预报，在一天时间内，三个方向上的位置预报误差均在5m以内，速度预报误差均在0.005m/s以内。在某些特殊情况下，如较大初始误差和较大测量噪声，ekf和ukf仍然有效。ukf在大噪声条件下的性能优于ekf，而ekf在较大初始误差下的性能优于ukf。
[0283]
对于实际应用，可以根据本文设置的仿真情况进行。当导航设备正常工作并可以获得连续有效的导航数据时，可以使用滤波方法或迭代方法解算星历，并根据需要实时更新以预报航天器的飞行状态。但是，迭代方法易受测量噪声的影响，而滤波方法在大噪声和大初始误差情况下仍能有效执行。因此，在正常、大噪声和大初始误差条件下，首先推荐滤波方法。当检测到某些严重的特殊情况，如导航设备出现故障，导航数据不连续，或者大量导航数据不可用时，建议使用迭代方法来解算航天器星历。
[0284]
本发明的实施例2提供了一种基于导航数据的应急航天器星历确定系统，所述系统包括：数据获取模块和航天器星历计算模块；
[0285]
数据获取模块，用于从航天器的导航设备获取导航数据，所述导航数据包括：航天器飞行状态矢量；
[0286]
航天器星历计算模块，用于当导航数据不足，则基于导航数据使用航天器星历迭代解算方法计算得到航天器的星历；否则，基于导航数据使用ukf方法计算得到航天器的星历；所述导航数据不足是指基于导航数据使用ukf方法无法收敛到所设定的精度。
[0287]
最后所应说明的是，以上实施例仅用以说明本发明的技术方案而非限制。尽管参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明的技术方案进行修改或者等同替换，都不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。