CN108490966A - 基于微分代数的静止轨道摄动相对轨迹高阶制导方法 - Google Patents
基于微分代数的静止轨道摄动相对轨迹高阶制导方法 Download PDFInfo
- Publication number
- CN108490966A CN108490966A CN201810093742.1A CN201810093742A CN108490966A CN 108490966 A CN108490966 A CN 108490966A CN 201810093742 A CN201810093742 A CN 201810093742A CN 108490966 A CN108490966 A CN 108490966A
- Authority
- CN
- China
- Prior art keywords
- guidance
- formula
- order
- initial
- track
- 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.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 45
- 230000003068 static effect Effects 0.000 title abstract description 4
- 238000013507 mapping Methods 0.000 claims abstract description 48
- 230000005855 radiation Effects 0.000 claims description 46
- 238000012546 transfer Methods 0.000 claims description 35
- 230000001133 acceleration Effects 0.000 claims description 11
- 239000011159 matrix material Substances 0.000 claims description 7
- 238000005312 nonlinear dynamic Methods 0.000 claims description 6
- 238000006467 substitution reaction Methods 0.000 claims description 6
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims description 3
- 238000000205 computational method Methods 0.000 claims description 3
- 238000000605 extraction Methods 0.000 claims description 3
- 230000000717 retained effect Effects 0.000 claims description 3
- 230000000977 initiatory effect Effects 0.000 claims description 2
- 230000007704 transition Effects 0.000 claims description 2
- 230000017105 transposition Effects 0.000 claims description 2
- 238000004364 calculation method Methods 0.000 abstract description 5
- 238000004422 calculation algorithm Methods 0.000 description 11
- 238000010586 diagram Methods 0.000 description 7
- 238000004458 analytical method Methods 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 3
- 230000000052 comparative effect Effects 0.000 description 2
- 239000012634 fragment Substances 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 230000003287 optical effect Effects 0.000 description 2
- 239000002245 particle Substances 0.000 description 2
- 230000008569 process Effects 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000004891 communication Methods 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000004069 differentiation Effects 0.000 description 1
- 239000000446 fuel Substances 0.000 description 1
- 230000003094 perturbing effect Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 230000004044 response Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G05—CONTROLLING; REGULATING
- G05D—SYSTEMS FOR CONTROLLING OR REGULATING NON-ELECTRIC VARIABLES
- G05D1/00—Control of position, course, altitude or attitude of land, water, air or space vehicles, e.g. using automatic pilots
- G05D1/10—Simultaneous control of position or course in three dimensions
- G05D1/101—Simultaneous control of position or course in three dimensions specially adapted for aircraft
Landscapes
- Engineering & Computer Science (AREA)
- Aviation & Aerospace Engineering (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Automation & Control Theory (AREA)
- Navigation (AREA)
Abstract
本发明涉及一种基于微分代数的静止轨道摄动相对轨迹高阶制导方法，属于航天器轨道动力学技术领域，考虑光压摄动地球静止轨道非线性相对运动轨迹的两脉冲制导问题，首先依据线性化的考虑光压摄动相对运动方程解析求解线性制导脉冲，然后利用微分代数方法，以线性解为参考轨迹求解制导所需速度增量关于初始和终端位置偏差的高阶泰勒多项式，最后将参考轨迹终端脱靶量代入该泰勒多项式，求得所需的制导脉冲，无需多次迭代积分，能够以任意高阶精度计算所需制导脉冲，计算精度高，应用标称问题的制导映射，仅需多项式运算即可求解标称问题附近的偏置制导问题，无需重新积分计算，计算精度和效率高。
Description
技术领域
本发明涉及航天器轨道动力学技术领域，具体涉及一种地球静止轨道航天器制导方法。
背景技术
地球静止轨道卫星在通信、气象、导航、预警等领域发挥着重要作用，面向静止轨道目标在轨服务和碎片清除等任务是当前航天技术发展的热点。实现航天器相对运动轨迹的精确制导，提高制导精度，有助于改善交会任务的燃料消耗、轨迹安全性、鲁棒性等性能指标，保障在轨服务、碎片清除等任务的顺利实施。
地球静止轨道卫星在运行时除了受到地球的中心天体引力外，还受到地球非球形、日月三体引力、太阳光压等摄动因素的影响。对于静止轨道相对运动轨迹，太阳光压摄动的影响远大于其他摄动项，在航天器相距几公里范围内时甚至大于中心引力差分的影响。分析表明，一个轨道周期内，太阳光压可导致静止轨道航天器相对轨迹产生数百米的偏差。因此，太阳光压摄动是进行静止轨道精确相对制导所必须考虑的因素，也是同近地轨道相对制导的主要区别之一。
交会任务中，传统相对轨迹制导算法主要基于线性化的相对动力学模型(CW方程)求解，忽略了高阶引力差分项和摄动因素的影响，终端误差较大。根据线性制导误差修正终端瞄准点，多次迭代求解制导脉冲，是一种有效消除终端脱靶量的方法。但是，该算法每迭代一次需要进行一次高精度轨道积分预报，而且对不同的边值条件需要重新迭代积分，计算量大。
发明内容
为了解决上述技术问题，本发明针对现有技术中的缺陷，提供了一种考虑太阳光压摄动的静止轨道航天器制导算法。
本发明是通过以下技术方案实现的：
一种基于微分代数的静止轨道摄动相对轨迹高阶制导方法，包括以下步骤：
S1：建立静止轨道线性相对运动方程，解析求解线性制导脉冲；
本步骤在基于CW方程的传统相对轨迹制导算法的基础上，将太阳光压摄动因素考虑进动力学模型中，重新解析求解相应的线性制导脉冲，以降低在地球静止轨道的制导误差；
具体的，所述S1包括以下步骤：
S101，建立目标航天器的相对轨道坐标系LVLH(Local Vertical LocalHorizontal)，给定初始参数：
该坐标系原点位于目标航天器的质心o，ox轴沿目标航天器的径向，oz轴沿目标航天器轨道面的法向，oy轴与ox、oz轴构成右手坐标系；给定航天器初始时刻t0的相对状态X0包括相对位置r0和相对速度v0，终端时刻tf的瞄准状态Xf包括相对位置rf和相对速度vf；
S102，建立考虑光压摄动的线性相对运动方程(1)：
ω为目标航天器平均轨道角速度，优选的，对标称地球静止轨道，ω＝7.2921×10-5rad/s；γ＝[γx,γy,γz]T为追踪航天器与目标航天器在惯性空间所受太阳光压摄动加速度矢量之差，即追踪航天器在LVLH坐标系内运动时所受的光压摄动加速度，上标T表示矩阵的转置，光压摄动加速度大小的计算方法如式(2)所示：
p为指向太阳单位面积上所受的光压摄动力，在1天文单位(平均日地距离)处p＝4.56×10-6N/m2；C为航天器表面的平均反射率，当光被全部吸收时C＝1，光被全部反射时C＝2；为航天器面向太阳方向的面质比，其中A为航天器面向太阳方向的横截面积，m为航天器的质量；下标c表示追踪航天器，下标t表示目标航天器；光压摄动方向的分析如附图2所示，LVLH坐标系各方向的光压摄动加速度分量如式(3)，
式(3)中，α为ox轴与光压摄动在目标航天器轨道平面内投影的夹角，β为光压摄动与目标航天器轨道平面的夹角，α称作太阳相位角，当地时间0:00时α＝0，当地时间24:00时α＝2π，设初始时刻t0时太阳相位角为α0，则任意时刻t的太阳相位角如式(4)所示，
α＝α0+ω(t-t0) (4)
由于标称地球静止轨道的轨道倾角为0，因此β即为太阳赤纬，一年中在区间[-23.5°,23.5°]内周期性变化，优选的，实际中静止轨道接近轨迹和制导任务的持续时间一般在一天或几天内，为方便解析求解，作近似处理将β视为常数，且不会引入太大的误差；
S103，推导线性相对运动方程(1)线性动力学模型的解析解：
X(t)＝Φ(t,t0)X(t0)+XSRP(t) (5)
其中，X(t)表示追踪航天器在t时刻LVLH坐标系中的相对状态，X(t0)即为初始时刻t0的相对状态X0，Φ(t,t0)为从初始时刻t0到时刻t的齐次方程的状态转移矩阵，XSRP(t)为时刻t的考虑光压摄动非齐次项的特解；
式(6)中，τ＝t-t0,s*＝sin(ωτ),c*＝cos(ωτ)；
式(7)中，rSRP为考虑光压摄动非齐次项的特解XSRP(t)的位置分量，vSRP为考虑光压摄动非齐次项的特解XSRP(t)的速度分量；
S104，求解线性制导脉冲和
将在终端时刻tf如式(6)的6×6状态转移矩阵Φ(tf,t0)分解为4个3×3的子矩阵，如式(8)所示，
则根据解析解如式(5)，从初始时刻t0到终端时刻tf的线性制导转移轨迹可表示为
式(9)中，线性制导转移轨迹的起始和末端状态如式(10)所示，
由式(9)，可求出初始时刻的线性制导脉冲为
终端时刻的线性制导脉冲为
为进一步提高相对制导精度，考虑相对动力学模型中的高阶引力差分项，下述步骤S2～S4/S4’建立了考虑光压摄动的静止轨道非线性相对运动方程，并以线性制导脉冲为初解，采用微分代数方法求出更为精确的非线性制导脉冲结果；
S2：利用非线性模型预报线性制导轨迹；
具体的，依据考虑光压摄动的非线性相对运动方程，代入线性制导脉冲，采用微分代数方法预报线性制导轨迹，得到预报映射和线性制导终端脱靶量；
进一步的，所述S2步骤包括：
S201，建立考虑光压摄动的静止轨道非线性相对运动方程：
式(13)中，μ为地球引力常数，μ＝3.986×1014m3/s2；a为目标航天器轨道半长轴，优选的，对标称地球静止轨道，a＝42164200m；
S202，以线性制导转移轨迹起始状态为标称值，考虑状态偏差，转移轨迹起始状态如式(14)所示
式(14)中，r1和v1为非线性动力学模型下转移轨迹的初始相对位置和相对速度，和为r1和v1的标称值，即式(10)所示线性转移轨迹初始状态，δr1和δv1为r1和v1相对标称值的偏差；
S203，将转移轨迹起始状态式(14)代入考虑光压摄动的非线性相对运动方程式(13)，应用微分代数方法，求出由转移轨迹起始状态偏差的高阶泰勒多项式描述的转移轨迹末端状态，如式(15)所示，微分代数方法中保留的高阶项阶数为N，N为正整数，综合考虑计算精度和计算量的要求，优选的，N＝3～8；
式(15)中，r2和v2为非线性动力学模型下转移轨迹的末端相对位置和相对速度，和为r2和v2的标称值，δr2和δv2为r2和v2相对标称值的偏差；和为高阶泰勒多项式形式的从起始位置和速度偏差到转移轨迹末端位置和速度偏差的映射，如式(16)所示，
微分代数方法是一种自动微分方法，能够在数值环境下计算非线性映射关于自变量的任意高阶偏导数，利用这些高阶偏导数可以对非线性系统进行任意阶次的泰勒展开，进而使用高阶多项式预报系统状态变量，该方法最初由美国密歇根大学粒子物理学教授Martin Berz提出，基本原理可见参考文献：Berz M.Differential algebraicdescription ofbeam dynamics to very high orders[J].Part.Accel.,1988,24(SSC-152):109-124，目前国内外尚未有应用该方法进行静止轨道航天器摄动高阶制导的研究或应用成果；
S204，求出线性制导的终端位置脱靶量如式(17)，
S3：构建和求解制导映射；
具体的，所述S3包括以下步骤：
S301，依据式(16)构建如下映射：
式(18)中，是转移轨迹起始位置偏差δr1的单位映射；
S302，采用微分代数中对高阶泰勒多项式求逆的方法(具体求逆方法见文献BerzM.Modern map methods in particle beam physics,Academic Press,London,1999.)，对式(18)求逆可得
提取式(19)中计算起始速度偏差δv1的项，得到式(20)由转移轨迹起始和终端位置偏差(δr2,δr1)到起始速度偏差δv1的高阶泰勒多项式制导映射
S4，求解得到初始和终端制导脉冲；
具体的，所述S4包括以下步骤：
S401，将线性制导终端脱靶量代入制导映射，得到高阶初始制导脉冲；
转移轨迹起始位置偏差转移轨迹终端位置偏差为式(17)的线性制导脱靶量代入制导映射式(20)，得到转移轨迹起始速度偏差
计算高阶初始制导脉冲如式(22)
S402，将转移轨迹起始位置偏差和速度偏差代入转移轨迹预报映射式(15)，得到高阶转移轨迹终端状态
计算高阶终端制导脉冲如式(24)：
上述方案可以很好地解决给定初始条件确定的标称情况下相对轨迹高阶制导问题。
进一步的，针对初始和终端瞄准位置发生变化的偏置制导问题，基于步骤S1-S3针对标称问题构建的预报映射和制导映射，本发明还公开另一种制导方法，可以直接求出偏置问题的解，无需重新计算，具体技术方案如下：
一种基于微分代数的静止轨道摄动相对轨迹高阶制导方法，包括以下步骤：
S1～S3同上；
S401’，将线性制导终端脱靶量与偏置问题的位置偏差和代入制导映射，得到高阶偏置初始制导脉冲；
转移轨迹起始位置偏差为偏置问题相对标称问题的起始位置偏差Δr1，转移轨迹终端位置偏差为线性制导脱靶量与偏置问题终端位置偏差Δr2之和，代入制导映射式(20)，得到偏置问题转移轨迹起始速度偏差
计算高阶偏置初始制导脉冲如式(26)
S402’，高阶偏置初始制导脉冲代入预报映射，得到高阶偏置终端制导脉冲；
将转移轨迹起始位置偏差和速度偏差代入转移轨迹预报映射式(15)，得到高阶转移轨迹终端状态
计算高阶偏置终端制导脉冲如式(28)
与现有技术相比，本发明的有益效果有：
1、本发明基于微分代数方法，能够以任意指定阶精度预报非线性系统的状态，因此高阶制导映射可以有效修正线性制导结果，消除终端误差，无需迭代步骤，即具有很高的制导精度；
2、本发明推导的高阶制导映射不仅能够精确求解标称的制导问题，而且对于包含初始和终端位置偏差的偏置制导问题，仅通过一些高效的多项式计算也能够以任意高阶精度求解，无需重新进行高精度轨道积分和迭代运算，兼具很高的计算精度和效率，具备应用于航天器在线自主计算的潜力；
3、本发明提出的算法虽然面向考虑光压摄动的地球静止轨道非线性相对轨迹制导问题，但是具有较强的扩展性，替换步骤S1-S3中的动力学模型和相应的线性化解析制导算法，本发明的算法即可拓展至其他非线性动力学模型下的高阶制导问题。
附图说明
图1为本发明的流程示意图；
图2为本发明步骤S102的光压摄动方向分析示意图；
图3为本发明末端制导轨迹对比图，验证了求解标称制导问题时，本发明的高精度和有效性；
图4为本发明不同初始位置制导轨迹对比图，验证了求解偏置制导问题时，本发明兼具高计算精度和高计算效率；
图5为本发明不同终端位置制导轨迹对比图，验证了求解偏置制导问题时，本发明兼具高计算精度和高计算效率。
具体实施方式
下面将结合本发明的附图，对本发明的技术方案进行清楚、完整地描述，显然，所描述的仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
实施例1
S1：建立静止轨道线性相对运动方程，解析求解线性制导脉冲；
S101，建立目标航天器的相对轨道坐标系LVLH(Local Vertical LocalHorizontal)，给定初始参数：
该坐标系原点位于目标航天器的质心o，ox轴沿目标航天器的径向，oz轴沿目标航天器轨道面的法向，oy轴与ox、oz轴构成右手坐标系；给定航天器初始时刻t0＝0的相对状态X0＝[0,-10km,0,0,0,0]，终端时刻tf＝86400s的瞄准状态Xf＝[0,-1km,0,0,0,0]；
S102，建立考虑光压摄动的线性相对运动方程(1)：
ω为目标航天器平均轨道角速度，对标称地球静止轨道，ω＝7.2921×10-5rad/s；γ＝[γx,γy,γz]T为追踪航天器与目标航天器在惯性空间所受太阳光压摄动加速度矢量之差，即追踪航天器在LVLH坐标系内运动时所受的光压摄动加速度，上标T表示矩阵的转置，光压摄动加速度大小的计算方法如式(2)所示：
下标c表示追踪航天器，下标t表示目标航天器；p为指向太阳单位面积上所受的光压摄动力，给定p＝4.56×10-6N/m2；C为航天器表面的平均反射率，给定Ct＝Cc＝1.2；为航天器面向太阳方向的面质比，给定光压摄动方向的分析如附图2所示，LVLH坐标系各方向的光压摄动加速度分量如式(3)，
式(3)中，α为ox轴与光压摄动在目标航天器轨道平面内投影的夹角，β为光压摄动与目标航天器轨道平面的夹角；α称作太阳相位角，设初始时刻t0时太阳相位角为α0，给定α0＝0，则任意时刻t的太阳相位角如式(4)所示，
α＝α0+ω(t-t0) (4)
由于标称地球静止轨道的轨道倾角为0，因此β即为太阳赤纬，一年中在区间[-23.5°,23.5°]内周期性变化，本实施例中，为方便解析求解，作近似处理将β视为常数，且不会引入太大的误差，给定β＝23.5°；
S103，推导线性相对运动方程(1)线性动力学模型的解析解：
X(t)＝Φ(t,t0)X(t0)+XSRP(t) (5)
其中，X(t)表示追踪航天器在t时刻LVLH坐标系中的相对状态，X(t0)即为初始时刻t0的相对状态X0，Φ(t,t0)为从初始时刻t0到时刻t的齐次方程的状态转移矩阵，XSRP(t)为时刻t的考虑光压摄动非齐次项的特解；
式(6)中，τ＝t-t0,s*＝sin(ωτ),c*＝cos(ωτ)；
式(7)中，rSRP为考虑光压摄动非齐次项的特解XSRP(t)的位置分量，vSRP为考虑光压摄动非齐次项的特解XSRP(t)的速度分量；
S104，求解线性制导脉冲和
将在终端时刻tf如式(6)的6×6状态转移矩阵Φ(tf,t0)分解为4个3×3的子矩阵，如式(8)所示，
则根据解析解如式(5)，从初始时刻t0到终端时刻tf的线性制导转移轨迹可表示为
式(9)中，线性制导转移轨迹的起始和末端状态如式(10)所示，
由式(9)，可求出初始时刻的线性制导脉冲为
终端时刻的线性制导脉冲为
S2：利用非线性模型预报线性制导轨迹：
S201：建立考虑光压摄动的静止轨道非线性相对运动方程：
式(13)中，μ为地球引力常数，μ＝3.986×1014m3/s2；a为目标航天器轨道半长轴，对标称地球静止轨道，a＝42164200m；
S202，以线性制导转移轨迹起始状态为标称值，考虑状态偏差，转移轨迹起始状态如式(14)所示
S203，将转移轨迹起始状态式(14)代入考虑光压摄动的非线性相对运动方程式(13)，应用微分代数方法，求出由转移轨迹起始状态偏差的高阶泰勒多项式描述的转移轨迹末端状态，如式(15)所示，微分代数方法中保留的高阶项阶数为N，给定N＝3
式(15)中，r2和v2为非线性动力学模型下转移轨迹的末端相对位置和相对速度，和为r2和v2的标称值，δr2和δv2为r2和v2相对标称值的偏差；和为高阶泰勒多项式形式的从起始位置和速度偏差到转移轨迹末端位置和速度偏差的映射，如式(16)所示，
S204，求出线性制导的终端位置脱靶量如式(17)，
S3：构建和求解制导映射；
S301，依据式(16)构建如下映射：
式(18)中，Ir1是转移轨迹起始位置偏差δr1的单位映射；
S302，采用微分代数中对高阶泰勒多项式求逆的方法(具体求逆方法见文献BerzM.Modern map methods in particle beam physics,Academic Press,London,1999.)，对式(18)求逆可得
提取式(19)中计算起始速度偏差δv1的项，得到式(20)由转移轨迹起始和终端位置偏差(δr2,δr1)到起始速度偏差δv1的高阶泰勒多项式制导映射
S4，求解得到初始和终端制导脉冲：
S401，将线性制导终端脱靶量代入制导映射，得到高阶初始制导脉冲；
转移轨迹起始位置偏差转移轨迹终端位置偏差为式(17)的线性制导脱靶量代入制导映射式(20)，得到转移轨迹起始速度偏差
计算高阶初始制导脉冲如式(22)
S402，将转移轨迹起始位置偏差和速度偏差代入转移轨迹预报映射式(15)，得到高阶转移轨迹终端状态
计算高阶终端制导脉冲如式(24)：
表1所示为本发明与传统制导方法的结果对比。参见表1和图3可知，线性制导结果与另两种算法的结果存在较大差别，终端制导误差也较大；迭代算法与本发明的高阶制导算法结果相近，均能将航天器准确制导至终端瞄准位置，有效消除线性制导误差。
本实施例验证了求解标称制导问题时，本发明基于微分代数的静止轨道摄动相对轨迹高阶制导方法的高精度和有效性。
表1标称问题各算法制导结果对比
实施例二：
本实施例的过程与实施例一基本相同，其主要不同点为：
给定初始相对位置分别在切向上-20km，-15km，-5km和5km处，终端瞄准位置分别在切向上-1km，1km和径向上-1km，1km处；
本实施例求解初始和终端瞄准位置发生变化的偏置制导问题，即执行步骤S1～S3后，执行S4’；
S401’，将线性制导终端脱靶量与偏置问题的位置偏差和代入制导映射，得到高阶偏置初始制导脉冲；
转移轨迹起始位置偏差为偏置问题相对标称问题的起始位置偏差Δr1，转移轨迹终端位置偏差为线性制导脱靶量与偏置问题终端位置偏差Δr2之和，代入制导映射式(20)，得到偏置问题转移轨迹起始速度偏差
计算高阶偏置初始制导脉冲如式(26)
S402’，高阶偏置初始制导脉冲代入预报映射，得到高阶偏置终端制导脉冲；
将转移轨迹起始位置偏差和速度偏差代入转移轨迹预报映射式(15)，得到高阶转移轨迹终端状态
计算高阶偏置终端制导脉冲如式(28)
相比实施例一而言，本实施例采用与实施例一相同的标称轨迹和制导映射求解位于不同初始和终端位置的偏置制导问题，在实施例一的步骤S3求出制导映射后，在实施例二执行步骤S4’将偏置制导问题的初始和终端条件代入制导映射通过多项式运算求解，无需重新进行高精度轨道积分和迭代修正等运算，能够显著提高计算效率。
参见图4所示不同初始位置制导轨迹对比图中，各条轨迹的终端瞄准点均在切向上-1km处，而标称轨迹的起始点在切向上-10km处，偏置轨迹的起始点分别在切向上-20km、-15km、-5km和5km处。参见图5所示不同终端位置制导轨迹对比图中，各条轨迹的起始点均在切向上-10km处，而标称轨迹的终端瞄准点在切向上-1km处，偏置轨迹的终端瞄准点分别在切向上1km处和径向上-1km、1km处。参见图4和图5可知，对偏置制导问题，本发明基于微分代数的静止轨道摄动相对轨迹高阶制导方法都能够精确求解，将航天器从不同的起始点制导到相应的终端瞄准点；而且该算法的计算过程中，步骤S3采用微分代数方法得到制导映射后，不管是对标称制导问题执行步骤S4，还是对偏置制导问题执行步骤S4’，仅涉及高效的多项式运算，不涉及高精度轨道积分和迭代修正运算，具有很高的求解效率。
本实施例验证了求解偏置制导问题时，本发明基于微分代数的静止轨道摄动相对轨迹高阶制导方法兼具高计算精度和高计算效率，有应用于航天器在线自主制导的潜在价值。
以上所述仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。
Claims (7)
1.一种基于微分代数的静止轨道摄动相对轨迹高阶制导方法，其特征在于，包括以下步骤：
S1：建立考虑光压摄动的静止轨道线性相对运动方程，解析求解线性制导脉冲；
S2：利用非线性模型预报线性制导轨迹；
所述S2，依据考虑光压摄动的非线性相对运动方程，代入线性制导脉冲，采用微分代数方法预报线性制导轨迹，得到预报映射和线性制导终端脱靶量；
S3：构建和求解制导映射；
S4，求解得到初始和终端制导脉冲。
2.根据权利要求1所述的一种基于微分代数的静止轨道摄动相对轨迹高阶制导方法，其特征在于，所述S1具体包括如下步骤：
S101，建立目标航天器的相对轨道坐标系LVLH(Local Vertical Local Horizontal)，给定初始参数：
该坐标系原点位于目标航天器的质心o，ox轴沿目标航天器的径向，oz轴沿目标航天器轨道面的法向，oy轴与ox、oz轴构成右手坐标系；给定航天器初始时刻t0的相对状态X0包括相对位置r0和相对速度v0，终端时刻tf的瞄准状态Xf包括相对位置rf和相对速度vf；
S102，建立考虑光压摄动的线性相对运动方程(1)：
ω为目标航天器平均轨道角速度，对标称地球静止轨道，ω＝7.2921×10-5rad/s；γ＝[γx,γy,γz]T为追踪航天器与目标航天器在惯性空间所受太阳光压摄动加速度矢量之差，即追踪航天器在LVLH坐标系内运动时所受的光压摄动加速度，上标T表示矩阵的转置，光压摄动加速度大小的计算方法如式(2)所示：
p为指向太阳单位面积上所受的光压摄动力，在1天文单位(平均日地距离)处p＝4.56×10-6N/m2；C为航天器表面的平均反射率，当光被全部吸收时C＝1，光被全部反射时C＝2；为航天器面向太阳方向的面质比，其中A为航天器面向太阳方向的横截面积，m为航天器的质量；下标c表示追踪航天器，下标t表示目标航天器；LVLH坐标系各方向的光压摄动加速度分量如式(3)，
式(3)中，α为ox轴与光压摄动在目标航天器轨道平面内投影的夹角，β为光压摄动与目标航天器轨道平面的夹角，α称作太阳相位角，当地时间0:00时α＝0，当地时间24:00时α＝2π，设初始时刻t0时太阳相位角为α0，则任意时刻t的太阳相位角如式(4)所示，
α＝α0+ω(t-t0) (4)
β即为太阳赤纬，一年中在区间[-23.5°,23.5°]内周期性变化；
S103，推导线性相对运动方程(1)线性动力学模型的解析解：
X(t)＝Φ(t,t0)X(t0)+XSRP(t) (5)
其中，X(t)表示追踪航天器在t时刻LVLH坐标系中的相对状态，X(t0)即为初始时刻t0的相对状态X0，Φ(t,t0)为从初始时刻t0到时刻t的齐次方程的状态转移矩阵，XSRP(t)为时刻t的考虑光压摄动非齐次项的特解；
式(6)中，τ＝t-t0,s*＝sin(ωτ),c*＝cos(ωτ)；
式(7)中，rSRP为考虑光压摄动非齐次项的特解XSRP(t)的位置分量，vSRP为考虑光压摄动非齐次项的特解XSRP(t)的速度分量；
S104，求解线性制导脉冲和
将在终端时刻tf如式(6)的6×6状态转移矩阵Φ(tf,t0)分解为4个3×3的子矩阵，如式(8)所示，
则根据解析解如式(5)，从初始时刻t0到终端时刻tf的线性制导转移轨迹可表示为
式(9)中，线性制导转移轨迹的起始和末端状态如式(10)所示，
由式(9)，可求出初始时刻的线性制导脉冲为
终端时刻的线性制导脉冲为
3.根据权利要求2所述的一种基于微分代数的静止轨道摄动相对轨迹高阶制导方法，其特征在于，所述S2具体包括如下步骤：
S201，建立考虑光压摄动的静止轨道非线性相对运动方程：
式(13)中，μ为地球引力常数，μ＝3.986×1014m3/s2；a为目标航天器轨道半长轴，对标称地球静止轨道，a＝42164200m；
S202，以线性制导转移轨迹起始状态为标称值，考虑状态偏差，转移轨迹起始状态如式(14)所示
式(14)中，r1和v1为非线性动力学模型下转移轨迹的初始相对位置和相对速度，和为r1和v1的标称值，即式(10)所示线性转移轨迹初始状态，δr1和δv1为r1和v1相对标称值的偏差；
S203，将转移轨迹起始状态式(14)代入考虑光压摄动的非线性相对运动方程式(13)，应用微分代数方法，求出由转移轨迹起始状态偏差的高阶泰勒多项式描述的转移轨迹末端状态，如式(15)所示，微分代数方法中保留的高阶项阶数为N，N为正整数；
式(15)中，r2和v2为非线性动力学模型下转移轨迹的末端相对位置和相对速度，和为r2和v2的标称值，δr2和δv2为r2和v2相对标称值的偏差；和为高阶泰勒多项式形式的从起始位置和速度偏差到转移轨迹末端位置和速度偏差的映射，如式(16)所示，
S204，求出线性制导的终端位置脱靶量如式(17)，
4.根据权利要求3所述一种基于微分代数的静止轨道摄动相对轨迹高阶制导方法，其特征在于，所述S3包括如下步骤：
S301，依据式(16)构建如下映射：
式(18)中，是转移轨迹起始位置偏差δr1的单位映射；
S302，采用微分代数中对高阶泰勒多项式求逆的方法，对式(18)求逆可得
提取式(19)中计算起始速度偏差δv1的项，得到式(20)由转移轨迹起始和终端位置偏差(δr2,δr1)到起始速度偏差δv1的高阶泰勒多项式制导映射
5.根据权利要求4所述的基于微分代数的静止轨道摄动相对轨迹高阶制导方法，其特征在于，所述S4包括如下步骤：
S401，将线性制导终端脱靶量代入制导映射，得到高阶初始制导脉冲；
转移轨迹起始位置偏差转移轨迹终端位置偏差为式(17)的线性制导脱靶量代入制导映射式(20)，得到转移轨迹起始速度偏差
计算高阶初始制导脉冲如式(22)
S402，将转移轨迹起始位置偏差和速度偏差代入转移轨迹预报映射式(15)，得到高阶转移轨迹终端状态
计算高阶终端制导脉冲如式(24)：
6.根据权利要求4所述的任意一种基于微分代数的静止轨道摄动相对轨迹高阶制导方法，其特征在于，所述S4包括以下步骤：
S401’，将线性制导终端脱靶量与偏置问题的位置偏差和代入制导映射，得到高阶偏置初始制导脉冲；
转移轨迹起始位置偏差为偏置问题相对标称问题的起始位置偏差Δr1，转移轨迹终端位置偏差为线性制导脱靶量与偏置问题终端位置偏差Δr2之和，代入制导映射式(20)，得到偏置问题转移轨迹起始速度偏差
计算高阶偏置初始制导脉冲如式(26)
S402’，高阶偏置初始制导脉冲代入预报映射，得到高阶偏置终端制导脉冲；
将转移轨迹起始位置偏差和速度偏差代入转移轨迹预报映射式(15)，得到高阶转移轨迹终端状态
计算高阶偏置终端制导脉冲如式(28)
7.根据权利要求3-6任一项所述一种基于微分代数的静止轨道摄动相对轨迹高阶制导方法，其特征在于，所述S203步中，微分代数方法中保留的高阶项阶数N＝3～8。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810093742.1A CN108490966B (zh) | 2018-01-31 | 2018-01-31 | 基于微分代数的静止轨道摄动相对轨迹高阶制导方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810093742.1A CN108490966B (zh) | 2018-01-31 | 2018-01-31 | 基于微分代数的静止轨道摄动相对轨迹高阶制导方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN108490966A true CN108490966A (zh) | 2018-09-04 |
| CN108490966B CN108490966B (zh) | 2021-02-05 |
Family
ID=63343934
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810093742.1A Active CN108490966B (zh) | 2018-01-31 | 2018-01-31 | 基于微分代数的静止轨道摄动相对轨迹高阶制导方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108490966B (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111731513A (zh) * | 2020-06-15 | 2020-10-02 | 航天东方红卫星有限公司 | 一种基于单脉冲轨控的高精度引力场中回归轨道维持方法 |
Citations (24)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP0785132A1 (en) * | 1996-01-16 | 1997-07-23 | Globalstar L.P. | Dynamic bias for orbital yaw steering |
| US6305646B1 (en) * | 1999-12-21 | 2001-10-23 | Hughes Electronics Corporation | Eccentricity control strategy for inclined geosynchronous orbits |
| US20040046691A1 (en) * | 2002-08-22 | 2004-03-11 | Barker Lee A. | Method for optimizing spacecraft yaw pointing to minimize geometric pointing error in antenna systems |
| CN101122780A (zh) * | 2007-09-04 | 2008-02-13 | 北京控制工程研究所 | 月球软着陆制导、导航与控制半物理仿真试验系统 |
| CN101226561A (zh) * | 2007-12-28 | 2008-07-23 | 南京航空航天大学 | 用于微型航天器姿态轨道控制系统的微型仿真支持系统及工作方法 |
| CN101762272A (zh) * | 2010-01-18 | 2010-06-30 | 哈尔滨工业大学 | 一种基于可观测度分析的深空自主导航方法 |
| CN106507769B (zh) * | 2009-05-08 | 2013-01-23 | 中国人民解放军国防科学技术大学 | 轨道机动受限的空间交会地面导引段轨道控制方法 |
| CN103093096A (zh) * | 2013-01-15 | 2013-05-08 | 北京航空航天大学 | 卫星轨道的确定方法和装置 |
| CN103729570A (zh) * | 2014-01-21 | 2014-04-16 | 山东大学 | 基于矩阵摄动理论的电力系统振荡模式的匹配方法 |
| CN103728980A (zh) * | 2014-01-08 | 2014-04-16 | 哈尔滨工业大学 | 航天器相对轨道的控制方法 |
| WO2014097263A1 (en) * | 2012-12-20 | 2014-06-26 | Thales Alenia Space Italia S.P.A. Con Unico Socio | Innovative orbit design for earth observation space missions |
| CN106508023B (zh) * | 2012-05-02 | 2014-08-20 | 中国人民解放军国防科学技术大学 | 一种交会航天器地面导引轨道机动故障判别方法 |
| CN104021241A (zh) * | 2014-05-13 | 2014-09-03 | 北京空间飞行器总体设计部 | 一种导航卫星复杂模型的地球反照光压摄动建模方法 |
| CN104309822A (zh) * | 2014-11-04 | 2015-01-28 | 哈尔滨工业大学 | 一种基于参数优化的航天器单脉冲水滴形绕飞轨迹悬停控制方法 |
| CN105183948A (zh) * | 2015-08-13 | 2015-12-23 | 北京空间飞行器总体设计部 | 一种基于二次反射的高精度卫星太阳光压摄动力建模方法 |
| CN105701283A (zh) * | 2016-01-08 | 2016-06-22 | 中国人民解放军国防科学技术大学 | 地球非球形摄动作用下自由段弹道误差传播的分析方法 |
| JP2016124538A (ja) * | 2015-01-07 | 2016-07-11 | 三菱電機株式会社 | 宇宙船の動作を制御する方法および制御システム、並びに宇宙船 |
| CN105865459A (zh) * | 2016-03-31 | 2016-08-17 | 北京理工大学 | 一种考虑视线角约束的小天体接近段制导方法 |
| CN106570285A (zh) * | 2016-11-09 | 2017-04-19 | 中国人民解放军装备学院 | 一种基于状态传递矩阵解析解的J2摄动Lambert问题求解方法 |
| CN106628257A (zh) * | 2016-09-28 | 2017-05-10 | 西北工业大学 | 地球摄动引力场中近地航天器相对运动轨道的保持方法 |
| CN106697333A (zh) * | 2017-01-12 | 2017-05-24 | 北京理工大学 | 一种航天器轨道控制策略的鲁棒性分析方法 |
| CN106970643A (zh) * | 2017-04-27 | 2017-07-21 | 中国人民解放军国防科学技术大学 | 一种解析的卫星非线性相对运动偏差传播分析方法 |
| CN107402903A (zh) * | 2017-07-07 | 2017-11-28 | 中国人民解放军国防科学技术大学 | 基于微分代数与高斯和的非线性系统状态偏差演化方法 |
| CN107506893A (zh) * | 2017-07-17 | 2017-12-22 | 中国人民解放军装备学院 | 一种太阳同步轨道航天器安全管理策略 |
-
2018
- 2018-01-31 CN CN201810093742.1A patent/CN108490966B/zh active Active
Patent Citations (24)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP0785132A1 (en) * | 1996-01-16 | 1997-07-23 | Globalstar L.P. | Dynamic bias for orbital yaw steering |
| US6305646B1 (en) * | 1999-12-21 | 2001-10-23 | Hughes Electronics Corporation | Eccentricity control strategy for inclined geosynchronous orbits |
| US20040046691A1 (en) * | 2002-08-22 | 2004-03-11 | Barker Lee A. | Method for optimizing spacecraft yaw pointing to minimize geometric pointing error in antenna systems |
| CN101122780A (zh) * | 2007-09-04 | 2008-02-13 | 北京控制工程研究所 | 月球软着陆制导、导航与控制半物理仿真试验系统 |
| CN101226561A (zh) * | 2007-12-28 | 2008-07-23 | 南京航空航天大学 | 用于微型航天器姿态轨道控制系统的微型仿真支持系统及工作方法 |
| CN106507769B (zh) * | 2009-05-08 | 2013-01-23 | 中国人民解放军国防科学技术大学 | 轨道机动受限的空间交会地面导引段轨道控制方法 |
| CN101762272A (zh) * | 2010-01-18 | 2010-06-30 | 哈尔滨工业大学 | 一种基于可观测度分析的深空自主导航方法 |
| CN106508023B (zh) * | 2012-05-02 | 2014-08-20 | 中国人民解放军国防科学技术大学 | 一种交会航天器地面导引轨道机动故障判别方法 |
| WO2014097263A1 (en) * | 2012-12-20 | 2014-06-26 | Thales Alenia Space Italia S.P.A. Con Unico Socio | Innovative orbit design for earth observation space missions |
| CN103093096A (zh) * | 2013-01-15 | 2013-05-08 | 北京航空航天大学 | 卫星轨道的确定方法和装置 |
| CN103728980A (zh) * | 2014-01-08 | 2014-04-16 | 哈尔滨工业大学 | 航天器相对轨道的控制方法 |
| CN103729570A (zh) * | 2014-01-21 | 2014-04-16 | 山东大学 | 基于矩阵摄动理论的电力系统振荡模式的匹配方法 |
| CN104021241A (zh) * | 2014-05-13 | 2014-09-03 | 北京空间飞行器总体设计部 | 一种导航卫星复杂模型的地球反照光压摄动建模方法 |
| CN104309822A (zh) * | 2014-11-04 | 2015-01-28 | 哈尔滨工业大学 | 一种基于参数优化的航天器单脉冲水滴形绕飞轨迹悬停控制方法 |
| JP2016124538A (ja) * | 2015-01-07 | 2016-07-11 | 三菱電機株式会社 | 宇宙船の動作を制御する方法および制御システム、並びに宇宙船 |
| CN105183948A (zh) * | 2015-08-13 | 2015-12-23 | 北京空间飞行器总体设计部 | 一种基于二次反射的高精度卫星太阳光压摄动力建模方法 |
| CN105701283A (zh) * | 2016-01-08 | 2016-06-22 | 中国人民解放军国防科学技术大学 | 地球非球形摄动作用下自由段弹道误差传播的分析方法 |
| CN105865459A (zh) * | 2016-03-31 | 2016-08-17 | 北京理工大学 | 一种考虑视线角约束的小天体接近段制导方法 |
| CN106628257A (zh) * | 2016-09-28 | 2017-05-10 | 西北工业大学 | 地球摄动引力场中近地航天器相对运动轨道的保持方法 |
| CN106570285A (zh) * | 2016-11-09 | 2017-04-19 | 中国人民解放军装备学院 | 一种基于状态传递矩阵解析解的J2摄动Lambert问题求解方法 |
| CN106697333A (zh) * | 2017-01-12 | 2017-05-24 | 北京理工大学 | 一种航天器轨道控制策略的鲁棒性分析方法 |
| CN106970643A (zh) * | 2017-04-27 | 2017-07-21 | 中国人民解放军国防科学技术大学 | 一种解析的卫星非线性相对运动偏差传播分析方法 |
| CN107402903A (zh) * | 2017-07-07 | 2017-11-28 | 中国人民解放军国防科学技术大学 | 基于微分代数与高斯和的非线性系统状态偏差演化方法 |
| CN107506893A (zh) * | 2017-07-17 | 2017-12-22 | 中国人民解放军装备学院 | 一种太阳同步轨道航天器安全管理策略 |
Non-Patent Citations (6)
| Title |
|---|
| CHUN-GUANG LI ETAL: "A VERSION OF ISOMAP WITH EXPLCIT MAPPING", 《PROCEEDINGS OF THE FIFTH INTERNATIONAL CONFERENCE ON MACHINE LEARNING AND CYBERNETICS》 * |
| GUANGYAN XU ETAL: "Fuel Optimal Maneuver of Passive and Periodic Satellite Formation:Linear Programming Approach", 《2013 FIFTH INTERNATIONAL CONFERENCE ON INTELLIGENT HUMAN-MACHINE SYSTEMS AND CYBERNETICS》 * |
| YA-ZHONG LUO ETAL: "A review of uncertainty propagation in orbital mechanics", 《PROGRESS IN AEROSPACE SCIENCES》 * |
| ZHEN YANG ETAL: "Two-level optimization approach for Mars orbital long-duration，large non-coplanar rendezvous phasing maneuvers", 《ADVANCE IN SPACE RESEARCH》 * |
| 刘光明等: "基于微分改正的Lambert拦截摄动修正方法研究", 《弹箭与制导学报》 * |
| 张进: "空间交会任务解析摄动分析与混合整数多目标规划方法", 《中国博士学位论文全文数据库 工程科技II辑》 * |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111731513A (zh) * | 2020-06-15 | 2020-10-02 | 航天东方红卫星有限公司 | 一种基于单脉冲轨控的高精度引力场中回归轨道维持方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN108490966B (zh) | 2021-02-05 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106697333B (zh) | 一种航天器轨道控制策略的鲁棒性分析方法 | |
| CN107797130B (zh) | 低轨航天器多点多参数轨道上行数据计算方法 | |
| CN106682274B (zh) | 考虑振幅约束的一种Halo轨道在轨保持方法 | |
| CN102819266B (zh) | 一种拟周期j2不变相对轨道编队飞行控制方法 | |
| CN109656133B (zh) | 一种针对空间走廊跟踪观测的分布式卫星群优化设计方法 | |
| CN110068845B (zh) | 一种基于平根数理论确定卫星理论轨道的方法 | |
| CN110989644A (zh) | 一种考虑目标点多终端约束的飞行器轨迹规划方法 | |
| CN102735231A (zh) | 一种提高光纤陀螺寻北仪精度的方法 | |
| CN107270933A (zh) | 一种基于多星协同的空间碎片运动状态联合确定方法 | |
| CN104501804A (zh) | 一种基于gps测量数据的卫星在轨轨道预报方法 | |
| CN102645223A (zh) | 一种基于比力观测的捷联惯导真空滤波修正方法 | |
| CN111680455B (zh) | 基于搭载形式的撞击探测轨道设计方法和系统 | |
| CN105069311A (zh) | 一种远程火箭发射初态误差传播估计方法 | |
| CN104121930B (zh) | 一种基于加表耦合的mems陀螺漂移误差的补偿方法 | |
| CN108490966A (zh) | 基于微分代数的静止轨道摄动相对轨迹高阶制导方法 | |
| CN104252548A (zh) | 一种燃料最优的火星探测器入轨目标点设计方法 | |
| CN108875174A (zh) | 一种基于多段打靶法的不变拟周期轨道确定方法 | |
| CN107506505B (zh) | 高精度地月自由返回轨道设计方法 | |
| CN103853047A (zh) | 一种基于状态量反馈的小推力跟踪制导方法 | |
| CN105354380A (zh) | 面向摄动因素影响补偿的滑翔弹道快速修正方法 | |
| CN111367305A (zh) | 一种高轨光压作用下导引伴飞稳定性控制方法及系统 | |
| CN103019251A (zh) | 一种强迫绕飞控制方法 | |
| CN106295218B (zh) | 一种快速确定能量最优拦截预测命中点的数值优化方法 | |
| CN114357853A (zh) | 基于遗传算法修正的多圈J2-Lambert转移轨道求解方法 | |
| CN115392540A (zh) | 一种用于月球轨道交会制导的快速预报方法 |
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 |