CN109032176B - 一种基于微分代数的地球同步轨道确定和参数确定方法 - Google Patents

一种基于微分代数的地球同步轨道确定和参数确定方法 Download PDF

Info

Publication number
CN109032176B
CN109032176B CN201810827230.3A CN201810827230A CN109032176B CN 109032176 B CN109032176 B CN 109032176B CN 201810827230 A CN201810827230 A CN 201810827230A CN 109032176 B CN109032176 B CN 109032176B
Authority
CN
China
Prior art keywords
spacecraft
orbit
earth
representing
state
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201810827230.3A
Other languages
English (en)
Other versions
CN109032176A (zh
Inventor
袁建平
陈建林
代洪华
孙冲
崔尧
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Northwestern Polytechnical University
Original Assignee
Northwestern Polytechnical University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Northwestern Polytechnical University filed Critical Northwestern Polytechnical University
Priority to CN201810827230.3A priority Critical patent/CN109032176B/zh
Publication of CN109032176A publication Critical patent/CN109032176A/zh
Application granted granted Critical
Publication of CN109032176B publication Critical patent/CN109032176B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G05—CONTROLLING; REGULATING
    • G05D—SYSTEMS FOR CONTROLLING OR REGULATING NON-ELECTRIC VARIABLES
    • G05D1/00—Control of position, course or altitude of land, water, air, or space vehicles, e.g. automatic pilot
    • G05D1/10—Simultaneous control of position or course in three dimensions

Abstract

本发明公开了一种基于微分代数的地球同步轨道确定和参数确定方法，选用基于地球同步卫星轨道要素描述的动力学模型，避免了在数值积分过程中使用较大的积分步长且避免动力学模型出现奇异点，并在动力学模型中加入摄动力项，进行多项式形式的积分，得到轨道状态和航天器参数，对得到的参数进行高阶预测，同时进行观测量的高阶预测；利用动力学模型和观测模型的非线性信息，提高估计精度，结合航天器的真实观测值，对轨道状态和航天器参数的高阶预测值进行更新并得到作为初值的高阶估计值，重复上面的实施过程，即可完成地球同步轨道确定和参数确定。本发明不仅可以提高轨道估计的精度，实现参数的高精度估计，还能大幅度降低计算的成本。

Description

一种基于微分代数的地球同步轨道确定和参数确定方法
技术领域
本发明涉及航空航天领域，特别涉及一种基于微分代数的地球同步轨道确定和参数确定方法。
背景技术
近年来，随着航天技术的不断发展和对地球同步轨道卫星需求的增加，位于地球同步轨道区域的物体数量不断增加。欧空间最新的空间物体监测报告显示当前大约有1533颗非涉密物体运行在同步轨道区域，其中502颗具有自控制能力，其余的为不可控卫星或空间碎片。为了避免地球同步轨道卫星之间的碰撞，保证它们的安全，监测它们的当前状态，预测并估计它们在未来一段时间内的状态演化具有十分重要的实际价值。
目前，关于地球同步卫星轨道的描述方法主要包括位置和速度矢量描述，经典的轨道六要素描述以及地球同步轨道要素描述。位置和速度矢量描述方法对航天器运动的描述具有明确的物理意义，然而它忽略了同步轨道卫星相对于地球固连坐标系几乎静止的特点。因此，在动力学模型积分过程中需要较小的积分步长，增加了计算成本。传统的开普勒轨道要素考虑了相对地球静止这一特点，解决了小积分步长引起计算成本增加的问题，然而，该轨道描述方法会在轨道倾角和偏心率为0时产生歧义。为了克服上面的缺点，研究人员发展了地球同步轨道要素。
一般而言，航天器轨道的观测数据由于受到观测器本身的精度和观测过程误差的影响，往往不能直接应用于航天器真实状态的计算。另一方面，航天器的动力学建模只考虑了主要的影响因素而忽略了包括高阶地球扁率摄动在内的一部分难以准确建模的因素，单纯的动力学积分也会存在较大偏差。常用的航天器状态估计方法往往融合观测数据和动力学模型预测。经典的方法是线性卡尔曼滤波或者扩展卡尔曼滤波方法，对于一般的动力学系统，它们能够给出一个较为精确的状态估计。然而，当动力学模型具有较强的非线性或者观测数据时间间隔较大时，这些方法的估计精度会降低，计算成本会增加。为了改善上述缺点，无味滤波器(UKF)得到了快速的发展，但是，大量的采样点会快速的增加无味滤波器的计算成本。另一方面，航天器动力学模型中往往存在航天器参数以及相应的不确定性，如航天器的面质比，在估计轨道的同时，较为准确的估计这些参数具有实际应用价值。由于地面观测设备存在误差，无法通过解算观测数据得到航天器的真实状态；其次，初始状态和航天器参数不可预测的偏差以及动力学模型的未建模因素，致使无法根据航天器的动力学模型计算未来某一时刻的航天器的准确状态。
发明内容
本发明的目的在于提供一种基于微分代数的地球同步轨道确定和参数确定方法，以克服现有技术的不足。
为达到上述目的，本发明采用如下技术方案：
一种基于微分代数的地球同步轨道确定和参数确定方法，首先根据地球同步轨道卫星的动力学特征，选用基于地球同步卫星轨道要素描述的动力学模型，并在动力学模型中加入太阳光压、第三引力体对近地航天器的引力摄动和地球非球形引力摄动三个摄动力项，在微分代数框架下，对加入摄动力项的动力学模型以及航天器参数变化模型进行多项式形式的积分，得到ti+1时刻以初始偏差为变量的k阶多项式表示的轨道状态和航天器参数，运用该k阶多项式形式的轨道状态和航天器参数，结合轨道状态和航天器参数的协方差矩阵，进行轨道状态和航天器参数的高阶预测；然后，根据观测方程和k阶多项式形式的状态解，得到观测量在ti+1时刻的多项式形式的解，结合轨道状态和航天器参数的协方差矩阵，进行观测量的高阶预测；然后结合ti+1时刻航天器的真实观测值，对ti+1时刻轨道状态和航天器参数的高阶预测值进行更新，输出ti+1时刻轨道状态和航天器参数的高阶估计值，最后，将该高阶估计值作为初值，重复上面的实施过程，即可完成地球同步轨道确定和参数确定。
进一步的，采用地球同步轨道要素描述的动力学模型，地球同步轨道要素根据经典的轨道六要素进行定义：
Figure BDA0001742769850000031
其中λ表示相对于本初子午线的恒星时角，
Figure BDA0001742769850000032
径向漂移率用标称地球同步轨道半长轴A＝42164.2km进行无量纲化，GA(t)＝GA(t0)+ωe(t-t0)表示格林尼治恒星时角，ωe表示地球自转角速度；ex、ey表示轨道偏心率e在x，y坐标轴上的投影；Q1、Q2是地球同步轨道要素集合的第五第六个要素，定义如上。使用泊松括号，推导出使用上述地球同步轨道要素描述的动力学模型为：
Figure BDA0001742769850000041
其中，a＝(ar,aθ,ah)表示摄动加速度沿着轨道径向，横向和法向的分量；s＝ω+Ω+ν表示航天器的恒星时角，ωe表示地球自转角速度，r表示航天器相对于地球质心的径向距离，p表示轨道的半通径，h表示轨道角动量的大小。
进一步的，通过地球同步轨道要素得到
Figure BDA0001742769850000042
Figure BDA0001742769850000043
Figure BDA0001742769850000044
其中，s＝ω+Ω+ν表示航天器的恒星时角，μ表示地球引力常数，
Figure BDA0001742769850000045
表示地球同步轨道要素。
进一步的，在微分代数框架下，对ti时刻的航天器状态的一个邻域进行积分，均值为0，标准差为σ的多元正态分布；得到ti+1时刻航天器状态的k阶泰勒多项式近似解，该解是ti时刻航天器状态偏差的函数，当航天器参数存在不确定性时，该k阶泰勒多项式近似解也是航天器参数偏差的函数；使用状态和参数共同组成的协方差矩阵，对该k阶泰勒多项式近似解进行高阶状态预测；同时，对控制参数变化的微分方程在时间区间[ti,ti+1]上进行多项式积分，得到ti+1时刻航天器参数的k阶泰勒多项式近似解，对航天器参数进行预测。
进一步的，使用经典轨道六要素作为桥梁，建立笛卡尔坐标
Figure BDA0001742769850000051
与地球同步轨道要素
Figure BDA0001742769850000052
之间的显式关系：
Figure BDA0001742769850000053
其中，s＝ω+Ω+ν表示航天器的恒星时角，ωe表示地球自转角速度，r表示航天器相对于地球质心的径向距离，p表示轨道的半通径，h表示轨道角动量的大小。
进一步的，摄动力建模，首先分析太阳光压加速度，其简化的动力学模型为
Figure BDA0001742769850000054
其中，aSRP,ECI表示太阳光压加速度在以地心为中心的惯性坐标系下的投影，rr表示太阳相对于航天器的相对位置，P表示在距太阳1个天文单位处单位面积上的太阳压力，Cr表示太阳压力系数，A表示航天器受照横截面积,m表示航天器质量，在对方程(2)进行积分过程中，将aSRP,ECI投影到航天器轨道的径向、横向和法向方向，其投影关系为：
Figure BDA0001742769850000061
其中，aSRP,LVLH表示太阳光压加速度在当地水平当地垂直坐标系中的投影，即沿着航天器轨道的径向，横向和法向方向，
Figure BDA0001742769850000062
表示从以地心为中心的惯性坐标系到以航天器质心为中心的当地水平当地垂直坐标系的转换矩阵，通过以下公式计算：
Figure BDA0001742769850000063
其中，I,J,K分别表示沿着惯性坐标系的坐标轴的单位矢量，i,j,k分别表示沿着当地水平当地垂直坐标系的坐标轴的单位矢量，它们的值分别为：
I＝(1,0,0),J＝(0,1,0),K＝(0,0,1)
Figure BDA0001742769850000064
j＝k×i
其中，r＝(x,y,z)表示航天器的位置，
Figure BDA0001742769850000065
表示航天器的速度，可由方程(8)计算得到；
质量为M的第三引力体对近地航天器的引力摄动加速度表示为：
Figure BDA0001742769850000066
其中，rM和r分别表示天体M和航天器相对于地心的位置矢量，μ＝GM表示天体M的引力常数；aM,ECI表示第三引力体(太阳或月球)对航天器产生的摄动加速度在地心惯性坐标系三个坐标轴上的投影；将aM,ECI投影到航天器轨道的径向、横向和法向方向，其投影关系为：
Figure BDA0001742769850000067
其中，aM,LVLH表示日引力摄动、月引力摄动加速度在当地水平当地垂直坐标系中的投影，即沿着航天器轨道的径向、横向和法向方向，
Figure BDA0001742769850000071
表示从以地心为中心的惯性坐标系到以航天器质心为中心的当地水平当地垂直坐标系的转换矩阵；在以地心为中心的地球固连坐标系中，5×5地球非球形摄动引力加速度在固连坐标系三个坐标轴上的投影anonspherical,ECF可表示为：
Figure BDA0001742769850000072
Figure BDA0001742769850000073
Figure BDA0001742769850000074
Figure BDA0001742769850000075
其中，GM表示地球引力常数，系数Cnm,Snm以及标称的地球半径
Figure BDA0001742769850000076
分别由EGM96地球引力场模型给出；Vn,m和Wn,m可通过下面迭代公式求得：
Figure BDA0001742769850000077
Figure BDA0001742769850000078
Figure BDA0001742769850000079
Figure BDA00017427698500000710
且
Figure BDA0001742769850000081
W00＝0,
Figure BDA0001742769850000082
W10＝0
其中，r＝(x,y,z)表示航天器在地球固连坐标系中的位置矢量，r＝|r|表示航天器到地心的距离，标称的地球半径
Figure BDA0001742769850000083
由EGM96地球引力场模型给出；将anonspherical,ECF投影到航天器轨道的径向、横向和法向方向，其投影关系为
Figure BDA0001742769850000084
其中，anonspherical,LVLH表示地球非球形引力摄动加速度在当地水平当地垂直坐标系中的投影，即沿着航天器轨道的径向、横向和法向方向；anonspherical,ECF表示地球非球形摄动引力加速度在固连坐标系三个坐标轴上的投影；
Figure BDA0001742769850000085
表示从以地心为中心的惯性坐标系到以航天器质心为中心的当地水平当地垂直坐标系的转换矩阵，由方程(9)计算得到；
Figure BDA0001742769850000086
表示从以地心为中心的地球固连坐标系到以地心为中心的惯性坐标系的转换矩阵，
Figure BDA0001742769850000087
其中，ωe表示地球自转角速度，t表示转过的时间。
进一步的，设一个n维变量组成的常微分动力学系统表示为：
Figure BDA0001742769850000088
该微分方程的解可表示为x(t)＝Φ(t；t0,x0)；假设t0时刻航天器的标称状态为
Figure BDA0001742769850000089
初始偏差为δx0，对某一数值初值
Figure BDA00017427698500000810
的邻域
Figure BDA00017427698500000811
进行积分，其中[x0]表示多项式，被称为微分代数变量，在区间[t0,t1]上进行积分后可得到微分方程在t1时刻的解Φ(t1；t0,x0+δx0)的高阶展开式
Figure BDA00017427698500000812
其中
Figure BDA00017427698500000813
其k阶多项式近似解可表示为
Figure BDA0001742769850000091
其中
Figure BDA0001742769850000092
表示标称轨迹状态；δx0＝[δx0,1，…,δx0,n]T表示初始偏差，
Figure BDA0001742769850000093
表示泰勒展开式的系数。
进一步的，在t1时刻，使用k阶多项式形式的解，即方程(22)，预测变量的均值
Figure BDA0001742769850000094
和协方差矩阵P1 -如下
Figure BDA0001742769850000095
Figure BDA0001742769850000096
其中，E[]表示期望值，ki和li表示变量偏差分量的指数，
Figure BDA0001742769850000097
表示过程噪声的协方差矩阵；注意方程(24)中的系数
Figure BDA0001742769850000098
与
Figure BDA0001742769850000099
除了
Figure BDA00017427698500000910
和
Figure BDA00017427698500000911
之外，
Figure BDA00017427698500000912
与对应的
Figure BDA00017427698500000913
均相等。
进一步的，设航天器的观测方程为：
zk+1＝h(xk+1,tk+1)+vk+1 (57)
其中，zk+1表示tk+1时刻的观测量，xk+1表示tk+1时刻的状态预测值，vk+1表示tk+1时刻的观测噪声；将方程(22)代入方程(25)得：
Figure BDA00017427698500000914
其中
Figure BDA00017427698500000915
表示对应于标称轨迹状态的观测值，系数
Figure BDA00017427698500000916
可通过将方程(22)代入方程(25)计算得到，m表示观测方程的数量；在t1时刻，使用k阶多项式形式的解，即方程(26)，预测观测值的均值
Figure BDA00017427698500000917
如下：
Figure BDA00017427698500000918
进一步的，当通过观测设备得到t1时刻新的观测值
Figure BDA00017427698500000919
时，将其融合到预测值中，其公式为：
Figure BDA0001742769850000101
Figure BDA0001742769850000102
K1＝P1 xz(P1 zz)-1 (62)
Figure BDA0001742769850000103
Figure BDA0001742769850000104
其中，系数
Figure BDA0001742769850000105
与
Figure BDA0001742769850000106
除了
Figure BDA0001742769850000107
和
Figure BDA0001742769850000108
之外，
Figure BDA0001742769850000109
与对应的
Figure BDA00017427698500001010
均相等，
Figure BDA00017427698500001011
和P1 +表示变量的均值和协方差矩阵的最终估计值，
Figure BDA00017427698500001012
表示观测噪声的协方差矩阵；得到变量在t1时刻的估计值；重复上述过程，得到任意时刻航天器的轨道状态和航天器参数。
与现有技术相比，本发明具有以下有益的技术效果：
本发明一种基于微分代数的地球同步轨道确定和参数确定方法，首先根据地球同步轨道卫星的动力学特征，选用基于地球同步卫星轨道要素描述的动力学模型，避免了在数值积分过程中使用较大的积分步长且避免动力学模型出现奇异点，并在动力学模型中加入太阳光压、第三引力体对近地航天器的引力摄动和地球非球形引力摄动三个摄动力项，在微分代数框架下，对加入摄动力项的动力学模型以及航天器参数变化模型进行多项式形式的积分，得到ti+1时刻以初始偏差为变量的k阶多项式表示的轨道状态和航天器参数，得到轨道状态和航天器参数的高阶预测；然后，根据观测方程和k阶多项式形式的状态解，得到观测量在ti+1时刻的多项式形式的解，结合轨道状态和航天器参数的协方差矩阵，进行观测量的高阶预测；利用动力学模型和观测模型的非线性信息，提高估计精度，然后结合ti+1时刻航天器的真实观测值，对ti+1时刻轨道状态和航天器参数的高阶预测值进行更新，输出ti+1时刻轨道状态和航天器参数的高阶估计值，最后，将该高阶估计值作为初值，重复上面的实施过程，即可完成地球同步轨道确定和参数确定。本发明不仅可以提高轨道估计的精度，实现参数的高精度估计，还能大幅度降低计算的成本。
进一步的，本发明基于微分代数技术的地球同步卫星高阶轨道确定和参数确定方法，与经典的卡尔曼滤波和扩展卡尔曼滤波算法相比，融合了动力学模型和观测模型的高阶信息，提高了航天器状态估计和参数估计的精度，同时当状态偏差较大时，该方法仍能收敛，因此，其在强非线性系统状态和参数估计，观测量时间间隔较长的状态估计等任务中，比经典的卡尔曼滤波算法表现更加完美。
附图说明
图1为基于微分代数方法的轨道不确定演化流程图。
图2为k阶滤波算法估计航天器的状态和面质比时，航天器面质比的估计误差随时间变化示意图。
图3为k阶滤波算法估计航天器的状态和面质比时，航天器位置的估计误差随时间变化示意图。
图4为k阶滤波算法估计航天器的状态和面质比时，航天器速度的估计误差随时间变化示意图。
具体实施方式
下面结合附图对本发明做进一步详细描述：
本发明一种基于微分代数的地球同步轨道确定和参数确定方法，首先根据地球同步轨道卫星的动力学特征，选用基于地球同步卫星轨道要素描述的动力学模型，并在动力学模型中加入太阳光压、日引力摄动、月引力摄动和地球非球形引力摄动四个摄动力项，其次，由于航天器初始状态和航天器参数存在初始偏差，假设初始偏差是均值为0，标准差为σ的多元正态分布，在微分代数框架下，对地球同步卫星动力学模型以及航天器参数变化模型进行多项式形式的积分，得到ti+1时刻以初始偏差为变量的k阶多项式表示的轨道状态和航天器参数，运用该高阶多项式形式的轨道状态和航天器参数，并结合轨道状态和航天器参数的协方差矩阵，进行轨道状态和航天器参数的高阶预测；然后，根据观测方程和上述k阶多项式形式的状态解，得到观测量在ti+1时刻的多项式形式的解，该解是状态偏差和参数偏差的函数，结合状态和参数的协方差矩阵，亦可进行观测量的高阶预测；结合ti+1时刻航天器的真实观测值，对ti+1时刻航天器状态和参数的高阶预测值进行更新，输出ti+1时刻航天器状态和参数的高阶估计值，最后，将该高阶估计值作为初值，重复上面的实施过程，即可完成地球同步轨道确定和参数确定，当k＝1时，该方法即为经典的扩展卡尔曼滤波算法。
具体包括以下步骤：
a、动力学模型选择：地球同步卫星在地球固连坐标系中运动缓慢，地球同步卫星轨道要素考虑了这个特征，因此可增加积分步长而不损失精度。该问题下，为了验证该方法在进行高阶状态估计的同时能进行参数估计的能力，考虑了太阳光压，日引力摄动、月引力摄动和地球扁率摄动等四个主要摄动力；
b、航天器高阶状态和参数预测：在微分代数框架下，对ti时刻的航天器状态的一个邻域(均值为0，标准差为σ的多元正态分布)进行积分，得到ti+1时刻航天器状态的k阶泰勒多项式近似解，一般来说，该解是ti时刻航天器状态偏差的函数，当航天器参数存在不确定性时，该k阶泰勒多项式近似解也是航天器参数偏差的函数；使用状态和参数共同组成的协方差矩阵，对该k阶泰勒多项式近似解进行高阶状态预测；同时，对控制参数变化的微分方程在时间区间[ti,ti+1]上进行多项式积分，得到ti+1时刻航天器参数的k阶泰勒多项式近似解，亦可对航天器参数进行预测；
c、航天器高阶观测量预测：将航天器状态的k阶多项式形式的解代入观测方程，得到ti+1时刻观测量的k阶多项式形式的解，该解是关于航天器ti时刻状态偏差和参数偏差的函数；使用状态和参数共同组成的协方差矩阵，对该k阶泰勒多项式近似解进行高阶状态预测并得到ti+1时刻航天器的观测量预测值；
d、航天器状态和参数更新：融合ti+1时刻航天器观测量的真实值，对航天器状态和参数进行更新，得到最终的航天器状态和参数估计值；
e、轨道状态估计和参数估计性能评估：在相同初始条件情况下，分别使用1阶，2阶，3阶方法对航天器状态和参数进行估计，并与真实的航天器状态和参数进行比较，分析该方法的性能。
基于微分代数技术的地球同步卫星高阶轨道确定和参数确定方法，与经典的卡尔曼滤波和扩展卡尔曼滤波算法相比，融合了动力学模型和观测模型的高阶信息，提高了航天器状态估计和参数估计的精度，同时当状态偏差较大时，该方法仍能收敛。因此，其在强非线性系统状态和参数估计，观测量时间间隔较长的状态估计等任务中，比经典的卡尔曼滤波算法表现更加完美。
选择动力学模型
根据任务要求，为了在数值积分过程中使用较大的积分步长且避免动力学模型出现奇异点，采用地球同步轨道要素描述的动力学模型。地球同步轨道要素可根据经典的轨道六要素进行定义
Figure BDA0001742769850000141
其中λ表示相对于本初子午线的恒星时角，
Figure BDA0001742769850000142
径向漂移率用标称地球同步轨道半长轴A＝42164.2km进行无量纲化，GA(t)＝GA(t0)+ωe(t-t0)表示格林尼治恒星时角，ωe表示地球自转角速度；ex、ey表示轨道偏心率e在x，y坐标轴上的投影；Q1、Q2是地球同步轨道要素集合的第五第六个要素，定义如上。使用泊松括号，推导出使用上述地球同步轨道要素描述的动力学模型为：
Figure BDA0001742769850000143
其中，a＝(ar,aθ,ah)表示摄动加速度沿着轨道径向，横向和法向的分量；s＝ω+Ω+ν表示航天器的恒星时角，ωe表示地球自转角速度，r表示航天器相对于地球质心的径向距离，p表示轨道的半通径，h表示轨道角动量的大小。它们可通过地球同步轨道要素得到
Figure BDA0001742769850000144
Figure BDA0001742769850000145
Figure BDA0001742769850000151
其中，s＝ω+Ω+ν表示航天器的恒星时角，μ表示地球引力常数，
Figure BDA0001742769850000152
表示地球同步轨道要素。
笛卡尔坐标和地球同步轨道要素之间的关系推导：
在对方程(2)进行积分过程中，需要计算航天器在时刻t的摄动力。然而，摄动力模型是航天器位置矢量的函数。因此，推导笛卡尔坐标和地球同步轨道要素之间显得尤为重要。使用经典轨道六要素作为桥梁，可以得到笛卡尔坐标
Figure BDA0001742769850000153
与地球同步轨道要素
Figure BDA0001742769850000154
之间的显式关系为
Figure BDA0001742769850000155
其中，s＝ω+Ω+ν表示航天器的恒星时角，ωe表示地球自转角速度，r表示航天器相对于地球质心的径向距离，p表示轨道的半通径，h表示轨道角动量的大小。
摄动力建模
针对地球同步卫星轨道的演化情形，主要考虑太阳光压，日引力摄动、月引力摄动和地球扁率这四个显著影响航天器轨道的摄动力。
首先分析太阳光压加速度，其简化的动力学模型为
Figure BDA0001742769850000161
其中，aSRP,ECI表示太阳光压加速度在以地心为中心的惯性坐标系下的投影，rr表示太阳相对于航天器的相对位置，AU表示一个天文单位，P表示在距太阳1个天文单位(1AU)处单位面积上的太阳压力，Cr表示太阳压力系数，A表示航天器受照横截面积,m表示航天器质量。值得注意的是，在对方程(2)进行积分过程中，需要将aSRP,ECI投影到航天器轨道的径向，横向和法向方向，其投影关系为
Figure BDA0001742769850000162
其中，aSRP,LVLH表示太阳光压加速度在当地水平当地垂直坐标系中的投影，即沿着航天器轨道的径向，横向和法向方向，
Figure BDA0001742769850000163
表示从以地心为中心的惯性坐标系到以航天器质心为中心的当地水平当地垂直坐标系的转换矩阵，可通过以下公式计算
Figure BDA0001742769850000164
其中，I,J,K分别表示沿着惯性坐标系的坐标轴的单位矢量，i,j,k分别表示沿着当地水平当地垂直坐标系的坐标轴的单位矢量，它们的值分别为
I＝(1,0,0),J＝(0,1,0),K＝(0,0,1)
Figure BDA0001742769850000165
j＝k×i
其中，r＝(x,y,z)表示航天器的位置，
Figure BDA0001742769850000166
表示航天器的速度，可由方程(8)计算得到。
日引力摄动和月引力摄动称为第三引力体对近地航天器的引力摄动，质量为M的第三引力体对近地航天器的引力摄动加速度表示为：
Figure BDA0001742769850000171
其中，rM和r分别表示天体M和航天器相对于地心的位置矢量，μ＝GM表示天体M的引力常数。aM,ECI表示第三引力体(太阳或月球)对航天器产生的摄动加速度在地心惯性坐标系三个坐标轴上的投影。类似于方程(8)，需要将aM,ECI投影到航天器轨道的径向，横向和法向方向，其投影关系为
Figure BDA0001742769850000172
其中，aM,LVLH表示日引力摄动、月引力摄动加速度在当地水平当地垂直坐标系中的投影，即沿着航天器轨道的径向，横向和法向方向，
Figure BDA0001742769850000173
表示从以地心为中心的惯性坐标系到以航天器质心为中心的当地水平当地垂直坐标系的转换矩阵，由方程(9)计算得到。
在以地心为中心的地球固连坐标系中，5×5地球非球形摄动引力加速度在固连坐标系三个坐标轴上的投影anonspherical,ECF可表示为
Figure BDA0001742769850000174
Figure BDA0001742769850000175
Figure BDA0001742769850000176
Figure BDA0001742769850000177
其中，GM表示地球引力常数，系数Cnm,Snm以及标称的地球半径
Figure BDA0001742769850000181
分别由EGM96地球引力场模型给出。Vn,m和Wn,m可通过下面迭代公式求得
Figure BDA0001742769850000182
Figure BDA0001742769850000183
Figure BDA0001742769850000184
Figure BDA0001742769850000185
且
Figure BDA0001742769850000186
W00＝0,
Figure BDA0001742769850000187
W10＝0
其中，r＝(x,y,z)表示航天器在地球固连坐标系中的位置矢量，r＝|r|表示航天器到地心的距离，标称的地球半径
Figure BDA0001742769850000188
由EGM96地球引力场模型给出。类似于方程(8)，需要将anonspherical,ECF投影到航天器轨道的径向，横向和法向方向，其投影关系为
Figure BDA0001742769850000189
其中，anonspherical,LVLH表示地球非球形引力摄动加速度在当地水平当地垂直坐标系中的投影，即沿着航天器轨道的径向，横向和法向方向；anonspherical,ECF表示地球非球形摄动引力加速度在固连坐标系三个坐标轴上的投影；
Figure BDA00017427698500001810
表示从以地心为中心的惯性坐标系到以航天器质心为中心的当地水平当地垂直坐标系的转换矩阵，由方程(9)计算得到；
Figure BDA00017427698500001811
表示从以地心为中心的地球固连坐标系到以地心为中心的惯性坐标系的转换矩阵，
Figure BDA0001742769850000191
其中，ωe表示地球自转角速度，t表示转过的时间。
常微分方程的高阶多项式近似解
基于微分代数技术，能够得到常微分方程解的任意高阶近似多项式解。其具体过程如下：设一个n维变量组成的常微分动力学系统(方程2)可表示为
Figure BDA0001742769850000192
因此，该微分方程的解可表示为x(t)＝Φ(t；t0,x0)。假设t0时刻航天器的标称状态为
Figure BDA0001742769850000193
初始偏差为δx0。不同于传统的数值积分方法只能从某一数值初值
Figure BDA0001742769850000194
开始积分，微分代数技术可以对某一数值初值
Figure BDA0001742769850000195
的邻域
Figure BDA0001742769850000196
进行积分，其中[x0]表示多项式，被称为微分代数变量。在区间[t0,t1]上进行积分后可得到微分方程在t1时刻的解Φ(t1；t0,x0+δx0)的高阶展开式
Figure BDA0001742769850000197
其中
Figure BDA0001742769850000198
其k阶多项式近似解可表示为
Figure BDA0001742769850000199
其中
Figure BDA00017427698500001910
表示标称轨迹状态；δx0＝[δx0,1，…,δx0,n]T表示初始偏差，
Figure BDA00017427698500001911
表示泰勒展开式的系数。值得注意的是，当航天器参数变化的微分方程需要被包含在方程(21)中，即方程(21)中的变量x既包含状态变量也包含参数变量。
航天器高阶状态和参数预测
在t1时刻，使用k阶多项式形式的解，即方程(22)，预测变量的均值
Figure BDA00017427698500001912
和协方差矩阵P1 -如下
Figure BDA0001742769850000201
Figure BDA0001742769850000202
其中，ki和li表示变量偏差分量的指数，
Figure BDA0001742769850000203
表示过程噪声的协方差矩阵。注意方程(24)中的系数
Figure BDA0001742769850000204
与
Figure BDA0001742769850000205
除了
Figure BDA0001742769850000206
和
Figure BDA0001742769850000207
之外，
Figure BDA0001742769850000208
Figure BDA0001742769850000209
与对应的
Figure BDA00017427698500002010
均相等。
航天器观测值预测
假设航天器的观测方程为
zk+1＝h(xk+1,tk+1)+vk+1 (89)
其中，zk+1表示tk+1时刻的观测量，xk+1表示tk+1时刻的状态预测值，vk+1表示tk+1时刻的观测噪声。使用方程(22)，代入方程(25)可得
Figure BDA00017427698500002011
其中
Figure BDA00017427698500002012
表示对应于标称轨迹状态的观测值，系数
Figure BDA00017427698500002013
可通过将方程(22)代入方程(25)计算得到，m表示观测方程的数量。在t1时刻，使用k阶多项式形式的解，即方程(26)，预测观测值的均值
Figure BDA00017427698500002014
如下
Figure BDA00017427698500002015
航天器状态和参数更新
当通过观测设备得到t1时刻新的观测值
Figure BDA00017427698500002016
时，将其融合到上面的预测值中，其公式为
Figure BDA00017427698500002017
Figure BDA00017427698500002018
K1＝P1 xz(P1 zz)-1 (94)
Figure BDA0001742769850000211
Figure BDA0001742769850000212
其中，系数
Figure BDA0001742769850000213
与
Figure BDA0001742769850000214
除了
Figure BDA0001742769850000215
和
Figure BDA0001742769850000216
之外，
Figure BDA0001742769850000217
与对应的
Figure BDA0001742769850000218
均相等，
Figure BDA0001742769850000219
和P1 +表示变量的均值和协方差矩阵的最终估计值，
Figure BDA00017427698500002110
表示观测噪声的协方差矩阵。至此，我们可以得到变量(包括轨道状态和航天器参数)在t1时刻的估计值；重复上述过程，可估计得到任意时刻航天器的轨道状态和航天器参数。
表1.初始条件和航天器参数
Figure BDA00017427698500002111
实施例：地球同步卫星轨道估计和面质比估计
基于微分代数技术的地球同步卫星轨道和参数高精度估计流程如图1所示，首先根据测量设备性能和观测精度，确定地球同步卫星的初始状态偏差，一般来说，该偏差可被假设为均值为0，标准差为σ的正态分布；然后在微分代数框架下，对地球同步卫星动力学模型以及航天器参数变化模型进行多项式形式的积分，得到ti+1时刻以初始偏差为变量的k阶多项式表示的轨道状态和参数值。运用该高阶多项式形式的航天器状态和参数，并结合状态和参数的协方差矩阵，进行状态和参数的高阶预测；然后，根据观测方程和上述k阶多项式形式的状态解，得到观测量在ti+1时刻的多项式形式的解，该解是状态偏差和参数偏差的函数，结合状态和参数的协方差矩阵，亦可进行观测量的高阶预测；结合ti+1时刻航天器的真实观测值，对ti+1时刻航天器状态和参数的高阶预测值进行更新，输出ti+1时刻航天器状态和参数的高阶估计值。最后，将该高阶估计值作为初值，重复上面的实施过程。以地球同步卫星轨道估计和面质比估计为例，在本例中，航天器的初始轨道参数如表1所示，我们假设航天器的位置和速度以及面质比存在初始偏差，位置偏差我们假设为均值为0，标准差为σx＝σy＝σz＝100km；速度偏差我们假设为均值为0，标准差为
Figure BDA0001742769850000221
面质比偏差我们假设为均值为0，标准差为σA/M＝0.01。本例中，基于微分代数框架的k阶滤波算法被应用于估计航天器的状态和面质比，其中k＝1,2,3，值得注意的是当k＝1时，该估计算法为经典的扩展卡尔曼滤波算法。图2给出了不同阶滤波算法估计航天器面质比的误差；图3和图4分别给出了不同阶滤波算法估计估计航天器位置和速度的误差示意图。图2-4共同表明，高阶的滤波算法相比于经典的扩展卡尔曼滤波算法可以显著提高估计精度。

Claims (10)

1.一种基于微分代数的地球同步轨道确定和参数确定方法，其特征在于，首先根据地球同步轨道卫星的动力学特征，选用基于地球同步卫星轨道要素描述的动力学模型，并在动力学模型中加入太阳光压、第三引力体对近地航天器的引力摄动和地球非球形引力摄动三个摄动力项，在微分代数框架下，对加入摄动力项的动力学模型以及航天器参数变化模型进行多项式形式的积分，得到ti+1时刻以初始偏差为变量的k阶多项式表示的轨道状态和航天器参数，得到轨道状态和航天器参数的高阶预测；然后，根据观测方程和k阶多项式形式的状态解，得到观测量在ti+1时刻的多项式形式的解，结合轨道状态和航天器参数的协方差矩阵，进行观测量的高阶预测；然后结合ti+1时刻航天器的真实观测值，对ti+1时刻轨道状态和航天器参数的高阶预测值进行更新，输出ti+1时刻轨道状态和航天器参数的高阶估计值，最后，将该高阶估计值作为初值，重复上面的实施过程，即可完成地球同步轨道确定和参数确定。
2.根据权利要求1所述的一种基于微分代数的地球同步轨道确定和参数确定方法，其特征在于，采用地球同步轨道要素描述的动力学模型，地球同步轨道要素根据经典的轨道六要素进行定义：
Figure FDA0003058542410000011
其中λ表示相对于本初子午线的恒星时角，
Figure FDA0003058542410000012
径向漂移率用标称地球同步轨道半长轴A＝42164.2km进行无量纲化，GA(t)＝GA(t0)+ωe(t-t0)表示格林尼治恒星时角，ωe表示地球自转角速度；ex、ey表示轨道偏心率e在x，y坐标轴上的投影；Q1、Q2是地球同步轨道要素集合的第五第六个要素，使用泊松括号，推导出使用上述地球同步轨道要素描述的动力学模型为：
Figure FDA0003058542410000021
其中，a＝(ar,aθ,ah)表示摄动加速度沿着轨道径向，横向和法向的分量；s＝ω+Ω+ν表示航天器的恒星时角，r表示航天器相对于地球质心的径向距离，p表示轨道的半通径，h表示轨道角动量的大小。
3.根据权利要求2所述的一种基于微分代数的地球同步轨道确定和参数确定方法，其特征在于，通过地球同步轨道要素得到
Figure FDA0003058542410000022
Figure FDA0003058542410000023
Figure FDA0003058542410000024
其中，μ表示地球引力常数，
Figure FDA0003058542410000025
表示地球同步轨道要素。
4.根据权利要求1所述的一种基于微分代数的地球同步轨道确定和参数确定方法，其特征在于，在微分代数框架下，对ti时刻的航天器状态的一个邻域进行积分，得到ti时刻的航天器状态的多元正态分布，其均值为0，标准差为σ；得到ti+1时刻航天器状态的k阶泰勒多项式近似解，该解是ti时刻航天器状态偏差的函数，当航天器参数存在不确定性时，该k阶泰勒多项式近似解也是航天器参数偏差的函数；使用状态和参数共同组成的协方差矩阵，对该k阶泰勒多项式近似解进行高阶状态预测；同时，对控制参数变化的微分方程在时间区间[ti,ti+1]上进行多项式积分，得到ti+1时刻航天器参数的k阶泰勒多项式近似解，对航天器参数进行预测。
5.根据权利要求2所述的一种基于微分代数的地球同步轨道确定和参数确定方法，其特征在于，使用经典轨道六要素作为桥梁，建立笛卡尔坐标
Figure FDA0003058542410000031
与地球同步轨道要素
Figure FDA0003058542410000032
之间的显式关系：
Figure FDA0003058542410000033
其中，μ表示地球引力常数。
6.根据权利要求5所述的一种基于微分代数的地球同步轨道确定和参数确定方法，其特征在于，摄动力建模，首先分析太阳光压加速度，其简化的动力学模型为
Figure FDA0003058542410000034
其中，aSRP,ECI表示太阳光压加速度在以地心为中心的惯性坐标系下的投影，rr表示太阳相对于航天器的相对位置，AU表示一个天文单位，P表示在距太阳1个天文单位处单位面积上的太阳压力，Cr表示太阳压力系数，A表示航天器受照横截面积,m表示航天器质量，在对方程(2)进行积分过程中，将aSRP,ECI投影到航天器轨道的径向、横向和法向方向，其投影关系为：
Figure FDA0003058542410000041
其中，aSRP,LVLH表示太阳光压加速度在当地水平当地垂直坐标系中的投影，即沿着航天器轨道的径向，横向和法向方向，
Figure FDA0003058542410000042
表示从以地心为中心的惯性坐标系到以航天器质心为中心的当地水平当地垂直坐标系的转换矩阵，通过以下公式计算：
Figure FDA0003058542410000043
其中，I,J,K分别表示沿着惯性坐标系的坐标轴的单位矢量，i,j,k分别表示沿着当地水平当地垂直坐标系的坐标轴的单位矢量，它们的值分别为：
I＝(1,0,0)，J＝(0,1,0)，K＝(0,0,1)
Figure FDA0003058542410000044
其中，r＝(x,y,z)表示航天器的位置，
Figure FDA0003058542410000045
表示航天器的速度，可由方程(6)计算得到；
质量为M的第三引力体对近地航天器的引力摄动加速度表示为：
Figure FDA0003058542410000046
其中，rM和r分别表示天体M和航天器相对于地心的位置矢量，μ＝GM表示天体M的引力常数；aM,ECI表示第三引力体对航天器产生的摄动加速度在地心惯性坐标系三个坐标轴上的投影；将aM,ECI投影到航天器轨道的径向、横向和法向方向，其投影关系为：
Figure FDA0003058542410000051
其中，aM,LVLH表示日引力摄动、月引力摄动加速度在当地水平当地垂直坐标系中的投影，即沿着航天器轨道的径向、横向和法向方向；在以地心为中心的地球固连坐标系中，5×5地球非球形摄动引力加速度在固连坐标系三个坐标轴上的投影anonspherical,ECF可表示为：
Figure FDA0003058542410000052
Figure FDA0003058542410000053
Figure FDA0003058542410000054
Figure FDA0003058542410000055
其中，GM表示地球引力常数，系数Cnm,Snm以及标称的地球半径
Figure FDA0003058542410000056
分别由EGM96地球引力场模型给出；Vn,m和Wn,m可通过下面迭代公式求得：
Figure FDA0003058542410000057
Figure FDA0003058542410000058
Figure FDA0003058542410000059
Figure FDA00030585424100000510
且
Figure FDA0003058542410000061
W10＝0
其中，r＝(x,y,z)表示航天器在地球固连坐标系中的位置矢量，r＝|r|表示航天器到地心的距离；将anonspherical,ECF投影到航天器轨道的径向、横向和法向方向，其投影关系为
Figure FDA0003058542410000062
其中，anonspherical,LVLH表示地球非球形引力摄动加速度在当地水平当地垂直坐标系中的投影，即沿着航天器轨道的径向、横向和法向方向；anonspherical,ECF表示地球非球形摄动引力加速度在固连坐标系三个坐标轴上的投影；
Figure FDA0003058542410000063
表示从以地心为中心的惯性坐标系到以航天器质心为中心的当地水平当地垂直坐标系的转换矩阵，由方程(9)计算得到；
Figure FDA0003058542410000064
表示从以地心为中心的地球固连坐标系到以地心为中心的惯性坐标系的转换矩阵，
Figure FDA0003058542410000065
其中，ωe表示地球自转角速度，t表示转过的时间。
7.根据权利要求1所述的一种基于微分代数的地球同步轨道确定和参数确定方法，其特征在于，设一个n维变量组成的常微分动力学系统表示为：
Figure FDA0003058542410000066
该微分方程的解可表示为x(t)＝Φ(t；t0,x0)；假设t0时刻航天器的标称状态为
Figure FDA0003058542410000067
初始偏差为δx0，对某一数值初值
Figure FDA0003058542410000068
的邻域
Figure FDA0003058542410000069
进行积分，其中[x0]表示多项式，被称为微分代数变量，在区间[t0,t1]上进行积分后可得到微分方程在t1时刻的解Φ(t1；t0,x0+δx0)的高阶展开式
Figure FDA00030585424100000610
其中
Figure FDA00030585424100000611
其k阶多项式近似解可表示为
Figure FDA0003058542410000071
其中
Figure FDA0003058542410000072
表示标称轨迹状态；δx0＝[δx0,1，…,δx0,n]T表示初始偏差，
Figure FDA0003058542410000073
表示泰勒展开式的系数。
8.根据权利要求7所述的一种基于微分代数的地球同步轨道确定和参数确定方法，其特征在于，在t1时刻，使用k阶多项式形式的解，即方程(22)，预测变量的均值
Figure FDA0003058542410000074
和协方差矩阵
Figure FDA00030585424100000719
如下
Figure FDA0003058542410000076
Figure FDA0003058542410000077
其中，E[]表示期望值，ki和li表示变量偏差分量的指数，
Figure FDA0003058542410000078
表示过程噪声的协方差矩阵；注意方程(24)中的系数
Figure FDA0003058542410000079
与
Figure FDA00030585424100000710
除了
Figure FDA00030585424100000711
和
Figure FDA00030585424100000712
之外，
Figure FDA00030585424100000713
与对应的
Figure FDA00030585424100000714
均相等。
9.根据权利要求7所述的一种基于微分代数的地球同步轨道确定和参数确定方法，其特征在于，设航天器的观测方程为：
zk+1＝h(xk+1,tk+1)+vk+1 (25)
其中，zk+1表示tk+1时刻的观测量，xk+1表示tk+1时刻的状态预测值，vk+1表示tk+1时刻的观测噪声；将方程(22)代入方程(25)得：
Figure FDA00030585424100000715
其中
Figure FDA00030585424100000716
表示对应于标称轨迹状态的观测值，系数
Figure FDA00030585424100000717
可通过将方程(22)代入方程(25)计算得到，m表示观测方程的数量；在t1时刻，使用k阶多项式形式的解，即方程(26)，预测观测值的均值
Figure FDA00030585424100000718
如下：
Figure FDA0003058542410000081
10.根据权利要求7所述的一种基于微分代数的地球同步轨道确定和参数确定方法，其特征在于，当通过观测设备得到t1时刻新的观测值
Figure FDA0003058542410000082
时，将其融合到预测值中，其公式为：
Figure FDA0003058542410000083
Figure FDA0003058542410000084
K1＝P1 xz(P1 zz)-1 (30)
Figure FDA0003058542410000085
Figure FDA0003058542410000086
其中，系数
Figure FDA0003058542410000087
与
Figure FDA0003058542410000088
除了
Figure FDA0003058542410000089
和
Figure FDA00030585424100000810
之外，
Figure FDA00030585424100000811
与对应的
Figure FDA00030585424100000812
均相等，
Figure FDA00030585424100000813
和P1 +表示变量的均值和协方差矩阵的最终估计值，
Figure FDA00030585424100000814
表示观测噪声的协方差矩阵；得到变量在t1时刻的估计值；重复上述过程，得到任意时刻航天器的轨道状态和航天器参数。
CN201810827230.3A 2018-07-25 2018-07-25 一种基于微分代数的地球同步轨道确定和参数确定方法 Active CN109032176B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810827230.3A CN109032176B (zh) 2018-07-25 2018-07-25 一种基于微分代数的地球同步轨道确定和参数确定方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810827230.3A CN109032176B (zh) 2018-07-25 2018-07-25 一种基于微分代数的地球同步轨道确定和参数确定方法

Publications (2)

Publication Number Publication Date
CN109032176A CN109032176A (zh) 2018-12-18
CN109032176B true CN109032176B (zh) 2021-06-22

Family

ID=64645039

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810827230.3A Active CN109032176B (zh) 2018-07-25 2018-07-25 一种基于微分代数的地球同步轨道确定和参数确定方法

Country Status (1)

Country Link
CN (1) CN109032176B (zh)

Families Citing this family (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110077627B (zh) * 2019-05-07 2020-08-18 北京航空航天大学 一种空间激光干涉引力波探测器轨道修正方法及系统
CN110781599B (zh) * 2019-10-30 2023-03-31 崔文惠 一种基于航天器多层绝热介质参数的温度估计方法
CN111125874B (zh) * 2019-11-18 2023-08-22 中国人民解放军63686部队 一种动平台高精度测轨预报方法
CN111027204B (zh) * 2019-12-05 2023-07-28 中国人民解放军63620部队 航天发射光、雷、遥与导航卫星测量数据融合处理方法
CN112141366B (zh) * 2020-09-23 2022-03-25 西北工业大学 一种地球轨道中航天器的轨迹预测方法及系统
CN112257343B (zh) * 2020-10-22 2023-03-17 上海卫星工程研究所 一种高精度地面轨迹重复轨道优化方法及系统
CN115320891B (zh) * 2022-10-12 2023-01-24 北京航天驭星科技有限公司 一种基于虚拟卫星的近圆标称轨道控制方法

Family Cites Families (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5100084A (en) * 1990-04-16 1992-03-31 Space Systems/Loral, Inc. Method and apparatus for inclined orbit attitude control for momentum bias spacecraft
US6237876B1 (en) * 2000-07-28 2001-05-29 Space Systems/Loral, Inc. Methods for using satellite state vector prediction to provide three-axis satellite attitude control
CN102116630A (zh) * 2009-12-31 2011-07-06 北京控制工程研究所 一种绕火星探测器的星上快速高精度确定方法
CN103675833A (zh) * 2013-11-27 2014-03-26 福建纳威导航科技有限责任公司 导航卫星轨道确定改进的代数技术
CN104006813A (zh) * 2014-04-03 2014-08-27 中国人民解放军国防科学技术大学 一种高轨卫星的脉冲星/星光角距组合导航方法
CN106202640B (zh) * 2016-06-28 2018-03-27 西北工业大学 日‑地三体引力场中的晕轨道航天器偏置轨道设计方法
CN107310752B (zh) * 2017-05-23 2020-09-08 西北工业大学 一种太阳帆航天器行星圆悬浮轨道之间的转移方法
CN107153209B (zh) * 2017-07-06 2019-07-30 武汉大学 一种短弧段低轨导航卫星实时精密定轨方法
CN107402903B (zh) * 2017-07-07 2021-02-26 中国人民解放军国防科学技术大学 基于微分代数与高斯和的非线性系统状态偏差演化方法

Also Published As

Publication number Publication date
CN109032176A (zh) 2018-12-18

Similar Documents

Publication Publication Date Title
CN109032176B (zh) 一种基于微分代数的地球同步轨道确定和参数确定方法
Gross et al. A comparison of extended Kalman filter, sigma-point Kalman filter, and particle filter in GPS/INS sensor fusion
CN112257343B (zh) 一种高精度地面轨迹重复轨道优化方法及系统
US5909381A (en) System of on board prediction of trajectories for autonomous navigation of GPS satellites
CN109255096A (zh) 一种基于微分代数的地球同步卫星轨道不确定演化方法
Karlgaard et al. Hyper-X post-flight trajectory reconstruction
Montenbruck An epoch state filter for use with analytical orbit models of low earth satellites
CN103438892B (zh) 一种改进的基于ekf的天文自主定轨算法
Liu et al. Autonomous orbit determination and timekeeping in lunar distant retrograde orbits by observing X‐ray pulsars
Rankin et al. VTXO-Virtual Telescope for X-Ray Observations
CN111854765B (zh) 一种中轨道导航卫星轨道长期预报方法
Golikov THEONA—a numerical-analytical theory of motion of artificial satellites of celestial bodies
Biggs et al. In-situ tracking of a solar sail’s characteristic acceleration using a robust active disturbance estimator
Da Forno et al. Autonomous navigation of MegSat1: Attitude, sensor bias and scale factor estimation by EKF and magnetometer-only measurement
Roscoe et al. Force modeling and state propagation for navigation and maneuver planning for cubesat rendezvous, proximity operations, and docking
CN112046793A (zh) 复杂扰动下的空间相对状态快速确定方法
CN112327333A (zh) 一种卫星位置计算方法及装置
Olson Sequential estimation methods for small body optical navigation
CN111547274A (zh) 一种航天器高精度自主目标预报方法
Harwood et al. Long-periodic and secular perturbations to the orbits of Explorer 19 and Lageos due to direct solar radiation pressure
Ibrahim Attitude and orbit control of small satellites for autonomous terrestrial target tracking
Gilani et al. Analysis of fidelities of linearized orbital models using least squares
Singh et al. Feasibility study of quasi-frozen, near-polar and extremely low-altitude lunar orbits
Li et al. X-ray pulsar based autonomous navigation for lunar rovers
Iwata Precision on-board orbit model for attitude control of the advanced land observing satellite (ALOS)

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant