CN112649006A - 一种太阳同步圆轨道的轨道规划方法 - Google Patents

一种太阳同步圆轨道的轨道规划方法 Download PDF

Info

Publication number
CN112649006A
CN112649006A CN202011590479.0A CN202011590479A CN112649006A CN 112649006 A CN112649006 A CN 112649006A CN 202011590479 A CN202011590479 A CN 202011590479A CN 112649006 A CN112649006 A CN 112649006A
Authority
CN
China
Prior art keywords
orbit
point
satellite
latitude
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.)
Pending
Application number
CN202011590479.0A
Other languages
English (en)
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.)
63921 Troops of PLA
Original Assignee
63921 Troops of PLA
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 63921 Troops of PLA filed Critical 63921 Troops of PLA
Priority to CN202011590479.0A priority Critical patent/CN112649006A/zh
Publication of CN112649006A publication Critical patent/CN112649006A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
    • G01C21/24—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 specially adapted for cosmonautical navigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Astronomy & Astrophysics (AREA)
  • Automation & Control Theory (AREA)
  • General Physics & Mathematics (AREA)
  • Radio Relay Systems (AREA)

Abstract

本发明公开了一种太阳同步圆轨道的轨道规划方法，能够结合卫星发射轨道模型、轨道回归特性约束、目标区域光照条件约束和精确地球运动模型(本设计考虑了岁差、章动、自转和极移等地球复杂运动模型)进行太阳同步圆轨道方案设计，实现对目标区域的快速有效探测。本卫星轨道规划设计所要解决的问题就是要根据所选取的目标区域信息、轨道回归特性的约束条件、侦察载荷的约束条件、发射部署参数和精确的地球运动模型等数据，完成针对不同类型目标区域探测的太阳同步圆轨道的规划设计。

Description

一种太阳同步圆轨道的轨道规划方法
技术领域
本发明涉及轨道规划技术领域，具体涉及一种太阳同步圆轨道的轨道规划方法。
背景技术
常规对地观测卫星轨道规划与应急卫星发射轨道规划的区别在于任务设计理念不同。常规对地观测任务首先发射对地观测卫星，在给出具体观测目标后，通过轨道机动将卫星送入观测轨道对目标进行观测。应急卫星发射任务的卫星轨道规划受到空间任务需求的牵引,通常根据侦察目标任务需求进行卫星轨道规划,权衡卫星发射轨道和运行轨道的各项参数,使之在满足应急发射的前提下,可以快速对目标区域进行有效探测。
对于应急卫星发射任务而言，需要综合考虑卫星的发射准备、发射等待和发射入轨过程中的约束，通过多种手段减小卫星发射和入轨各阶段的时间，从而实现对用户需求的快速响应。将卫星送入高轨轨道需要采用大型运载，而大型运载的发射准备时间较长。同时，高轨卫星的发射入轨方式比低轨卫星复杂，卫星入轨时间也较长。当前应急卫星发射通常采用低轨轨道。应急卫星发射过程中在还需要根据发射态势情况计算确定发射窗口，发射窗口的时间和大小由具体任务确定。在应急卫星发射任务中，应根据发射、测控和光照等约束条件具体分析计算其发射窗口。
传统应急卫星发射轨道规划主要是针对入轨后的卫星轨道进行规划设计，通常都没有综合考虑发射轨道模型、轨道回归特性、光学载荷的光照条件约束和地球复杂运动模型等多种因素的影响。
因此，如何克服传统应急卫星发射轨道设计的局限性，结合了卫星发射轨道模型、轨道回归特性约束、目标区域光照条件约束和精确地球运动模型(本设计考虑了岁差、章动、自转和极移等地球复杂运动模型)进行轨道方案设计，实现对目标区域的快速有效探测，是目前亟待解决的问题。
发明内容
有鉴于此，本发明提供了一种太阳同步圆轨道的轨道规划方法，能够结合卫星发射轨道模型、轨道回归特性约束、目标区域光照条件约束和精确地球运动模型(本设计考虑了岁差、章动、自转和极移等地球复杂运动模型)进行太阳同步圆轨道方案设计，实现对目标区域的快速有效探测。
为达到上述目的，本发明的技术方案包括如下步骤：
步骤1、太阳同步圆轨道具有两种入轨模式：即逆行下降和逆行上升2种入轨模式；根据轨道半长轴的上下限，获得逆行上升和逆行下降的轨道倾角搜索范围。
如果卫星载荷类型为CCD光学设备，则以发射部署完成时刻为起点，计算目标点区域24小时内太阳高度角满足门限要求的UTC时间观测窗口。
该轨道规划方法针对逆行上升和逆行下降轨道的倾角搜索范围内所有卫星轨道，计算卫星轨道星下点轨迹16次过境目标纬度圈时星下点的经度分布，并且选取满足侦察要求的最优轨道。设定初始状态为，MO为入轨模式，MO的取值为1或2，分别指代逆行下降和逆行上升2种入轨模式，MO初始值设定为1；j为过境次数，j的初值设定为1；BackupOrbitNum为轨道方案个数，BackupOrbitNum的初值设定为0。
步骤2、判断j的值是否满足j<＝16，若是则进入步骤3，否则进入步骤11；
步骤3、判断MO的值是否满足MO<＝2；若是则进入步骤4，否则将MO重新设定为1，j自增1，返回步骤2。
步骤4、对于第MO种入轨模式轨道倾角遍历范围内所有卫星轨道，计算第j次过境目标纬度圈时星下点经度分布,选取轨道倾角遍历范围内过境时刻星下点经度与目标点经度偏差值最小的轨道。
步骤5、判断步骤4中所选取的偏差值最小的轨道对应的偏差值是否满足预设的指标要求，若是则进入步骤6，否则MO自增1，返回步骤3。
步骤6、判断当前卫星载荷类型是否为CCD光学设备，若是则进入步骤7，否则进入步骤10。
步骤7、判断当前轨道卫星过境目标点时刻是否在CCD光学设备的UTC时间观测窗口内，若是则进入步骤10，否则进入步骤8。
步骤8、计算发射延时，对于当前入轨模式轨道倾角遍历范围内的所有卫星轨道，计算当前第j次过境目标纬度圈时星下点经度分布,选取遍历范围内过境时刻星下点经度与目标点经度偏差值最小的轨道。
步骤9、判断步骤8中所选取的偏差值最小的轨道是否满足预设的指标要求，若是则进入步骤10，否则返回步骤3。
步骤10、当前第BackupOrbitNum个轨道规划方案满足初步要求，将当前第BackupOrbitNum个轨道规划方案存入在BackupOrbit数组中,并且BackupOrbitNum自增1。
步骤11、太阳同步圆轨道规划完成，将当前对应的共BackupOrbitNum个轨道方案存入预先构建的BackupOrbit数组中。
进一步地，步骤4中，对于第MO种入轨模式轨道倾角遍历范围内所有卫星轨道，计算第j次过境目标纬度圈时星下点经度分布,选取轨道倾角遍历范围内过境时刻星下点经度与目标点经度偏差值最小的轨道，具体包括如下步骤：
步骤101、遍历范围内太阳同步圆轨道星下点过境目标纬度圈地心经度分布；
步骤102、轨道倾角遍历范围的参数初始化：轨道倾角遍历初始值domainDegreeLower,遍历间隔0.01度,遍历次数indexNum,索引index＝1；
步骤103、判断是否满足index<＝indexNum；若是则进入步骤104；否则进入步骤105。
步骤104、对当前轨道倾角为IiAngleDegree＝domainDegreeLower+0.01×(index-1)的卫星发射轨道和卫星运行轨道，按照如下S1041～S1043计算之后进入步骤105。
S1041、计算对应当前轨道倾角倾角IiAngleDegree的太阳同步圆轨道半长轴。
S1042、根据当前轨道倾角IiAngleDegree和太阳同步圆轨道半长轴，进行太阳同步圆轨道星下点轨迹计算，得到卫星星下点轨迹多次过境目标纬度圈时的星下点地心经度与目标点地心经度的偏差数据。
S1043、对索引值index自增1，返回步骤103。
步骤105、对于遍历范围内indexNum条太阳同步圆轨道星下点轨迹第j次过境目标纬度圈时刻的星下点地心经度与目标点地心经度的偏差分布进行统计分析处理，选取偏差值最小的轨道方案,最小偏差值为minValue。
步骤106、对最小偏差值minValue进行门限threshold判决，若minValue<＝threshold，则进入步骤107；否则返回无效值。
步骤107、minValue对应的轨道方案保存在minTimeOrbit数组中，并返回有效轨道倾角值。
进一步地，S1041、计算对应当前轨道倾角倾角IiAngleDegree的太阳同步圆轨道半长轴，具体为如下步骤：
S10411、根据从界面输入的轨道入轨点高度上限和下限分别计算太阳同步圆轨道半长轴上限aU和下限aL；根据当前模块的输入参数轨道倾角iiAngleDegreeTemp和轨道半长轴上下限，分别计算当前轨道升交点赤经导数与太阳同步轨道标准升交点赤经导数的差值tempValueU和tempValueL。
S10412、判断若(tempValueL×tempValueU)>0.0，则输出无效值，否则进入S10413。
S10413、轨道半长轴迭代计算初始化:index＝0,差值tempValueL和tempValueU
S10414、判断index<30是否成立，若是进入S10415，否则输出轨道半长轴中间值aM。
S10415、进行轨道半长轴迭代计算：首先计算轨道半长轴中间值aM＝(aU+aL)/2.0；然后根据当前轨道倾角iiAngleDegreeTemp和轨道半长轴中间值aM，计算当前轨道升交点赤经导数与太阳同步轨道标准升交点赤经导数的差值tempValue；index自增1。
S10416、判断(tempValue×tempValueL)>0.0是否成立，若是则令tempValueL＝tempValue；同时令aL＝aM，否则令tempValueU＝tempValue，同时令aU＝aM；返回S10414。
进一步地，S10411和S10415中，计算当前轨道升交点赤经导数与太阳同步轨道标准升交点赤经导数的差值，具体为：
根据轨道半长轴aI、轨道倾角i、地球平均半径Re、地球自转速度ωE和卫星平均角速度n，计算当前轨道升交点赤经导数与太阳同步圆轨道标准升交点赤经导数的偏差值：
Figure BDA0002866868880000051
进一步地，S1042中，根据当前轨道倾角IiAngleDegree和太阳同步圆轨道半长轴，进行太阳同步圆轨道星下点轨迹计算，得到卫星星下点轨迹多次过境目标纬度圈时的星下点地心经度与目标点地心经度的偏差数据；具体为：
S10421、卫星发射轨道参数初始化计算部分：
1.发射点地心经纬度为
Figure BDA0002866868880000061
发射轨道地心角βL和发射轨道倾角i为模块输入参数；
2.发射时刻卫星在发射轨道上真近点角fLradian和入轨时刻卫星在发射轨道上真近点角fIradian的计算；
3.发射轨道偏心率eL、发射轨道半通径p和发射轨道半长轴aL的计算
4.发射点偏近点角ELangleRadian和发射点平近点角MLangleRadian的计算
5.根据上述参数计算卫星入轨时刻tI的计算
6.发射时刻GMST格林尼治平恒星时GL的计算
7.发射时刻发射点在ECI地心惯性坐标系中赤经为：
αL＝OnLimitPiRadian(GL+λL)，其中角度约束函数OnLimitPiRadian(x)将角度x约束到[-π,π)范围内。
S10422、若为下降入轨模式，发射时刻卫星纬度辐角计算如下:
Figure BDA0002866868880000062
轨道升交点赤经:Ω＝OnLimit2PiRadian(αL-tan-1(tan(uL)×cos(i))+π)；其中角度约束函数OnLimit2PiRadian(x)将角度x约束到[0，2π)范围内。
若为上升入轨模式，发射时刻卫星纬度辐角计算如下:
Figure BDA0002866868880000063
轨道升交点赤经:Ω＝OnLimit2PiRadian(αL-tan-1(tan(uL)×cos(i))+0)；
S10423、卫星运行轨道参数初始化部分:
1.入轨时刻卫星纬度辐角uI＝uL+βL,目标点地心经纬度
Figure BDA0002866868880000064
卫星平均角速度n,地球自转转速ωE；
2.入轨点地心纬度为
Figure BDA0002866868880000065
S10424、若为下降入轨模式，则分如下4种情况计算卫星入轨后第1次和第2次过境目标纬度圈的卫星纬度辐角值uT1和uT2；目标点地心纬度为
Figure BDA0002866868880000066
入轨点地心维度
Figure BDA0002866868880000071
第1种情况.目标点在北半球,目标点地心纬度小于等于入轨点地心纬度；
第2种情况.目标点在南半球,目标点地心纬度小于等于入轨点地心纬度；
第3种情况.目标点在北半球,目标点地心纬度大于入轨点地心纬度；
第4种情况.目标点在南半球,目标点地心纬度大于入轨点地心纬度；
若
Figure BDA0002866868880000072
即目标点地心纬度小于等于入轨点地心纬度，则对于入轨后第1次过境，
Figure BDA0002866868880000073
包括第1种情况和第2种情况，如果OnLimit2PiRadian(uT1)<π则属于第1种情况，第一中间变量为Δu＝π-uT1和第二中间变量为utemp＝π+2Δu；否则属于第2种情况：Δu＝uT1-π；utemp＝π-2Δu；则入轨后第2次过境uT2＝uT1+utemp；
若
Figure BDA0002866868880000074
即目标点地心纬度大于入轨点地心纬度，则对于入轨后第1次过境，
Figure BDA0002866868880000075
包括第3种情况和第4种情况：
如果OnLimit2PiRadian(uT1)<π则属第3种情况：以指代量tempRadian指代OnLimit2PiRadian(uT1)，Δu＝tempRadian；utemp＝π-2Δu；否则属于第4种情况：Δu＝fabs(tempRadian-2π)；utemp＝π+2Δu；则入轨后第2次过境uT2＝uT1+utemp；
若为上升入轨模式，分如下6种情况计算卫星入轨后第1次和第2次过境目标纬度圈的卫星纬度辐角值uT1和uT2；
第5种情况.目标点在北半球,目标点纬度大于入轨点纬度,入轨点在下降段；
第6种情况.目标点在南半球,目标点纬度大于入轨点纬度,入轨点在下降段；
第7种情况.目标点在北半球,目标点纬度大于入轨点纬度,入轨点在上升段；
第8种情况.目标点在南半球,目标点纬度大于入轨点纬度,入轨点在上升段；
第9种情况.目标点在北半球,目标点纬度小于入轨点纬度；
第10种情况.目标点在南半球,目标点纬度小于入轨点纬度；
若
Figure BDA0002866868880000081
即目标点地心纬度大于等于入轨点地心纬度，且入轨点在下降段
Figure BDA0002866868880000082
则对于入轨后第1次过境，
Figure BDA0002866868880000083
tempRadian＝OnLimit2PiRadian(uT1)；
若tempRadian＜π属于第5种情况，Δu＝tempRadian；utemp＝π-2Δu；否则属于第6种情况，则Δu＝fabs(tempRadian-2π)；utemp＝π+2Δu；
若
Figure BDA0002866868880000084
即目标点地心纬度大于等于入轨点地心纬度，且入轨点在上升段
Figure BDA0002866868880000085
则对于入轨后第1次过境，
Figure BDA0002866868880000086
tempRadian＝OnLimit2PiRadian(uT1)；
若tempRadian＜π属于第7种情况，Δu＝uT1；utemp＝π-2Δu；否则属于第8种情况，第8种情况不存在；
则入轨后第2次过境uT2＝uT1+utemp；
若
Figure BDA0002866868880000087
即目标点地心纬度小于入轨点地心纬度，则对于入轨后第1次过境，
Figure BDA0002866868880000088
tempRadian＝OnLimit2PiRadian(uT1)；
若tempRadian＜π属于第9种情况，
Figure BDA0002866868880000089
utemp＝π+2Δu；否则属于第10种情况，
Figure BDA00028668688800000810
utemp＝π+2Δu；
则入轨后第2次过境uT2＝uT1+utemp；
S10425、卫星运行轨道星下点过境目标纬度圈时刻地心经度计算:
1.计算卫星运行轨道周期tΩ、升交点赤经导数
Figure BDA00028668688800000811
和近地点辐角导数
Figure BDA00028668688800000812
2.计算卫星入轨后16次过境目标纬度圈的时刻:
t1＝(uT1-uI)/n,t2＝(uT2-uI)/n,t2m+1＝t1+mtΩ,t2m+2＝t2+mtΩ,m＝1,2,...,7
3.计算入轨后第1次和第2次过境目标纬度圈时星下点地心经度LonT(1)和LonT(2)；
4.计算入轨后第3次到第16次卫星星下点轨迹过境目标纬度圈时刻的地心经度:
Figure BDA0002866868880000091
Figure BDA0002866868880000092
5.计算卫星星下点轨迹16次过境目标纬度圈时刻地心经度与目标地心经度的偏差。
进一步地，计算入轨后第1次和第2次过境目标纬度圈时星下点地心经度LonT(1)和LonT(2)；包括如下步骤：
SA1.卫星半长轴aI和当前时刻纬度辐角uT,卫星在地心轨道坐标系下的位置矢量r＝[aI cos(uT) aI sin(uT) 0]′；
SA2.根据相对入轨时刻时间差Δt和升交点赤经导数
Figure BDA0002866868880000093
计算当前时刻卫星运行轨道的升交点赤经
Figure BDA0002866868880000094
其中ΩI为卫星入轨时刻的升交点赤经；
SA3.根据相对入轨时刻时间差Δt和近地点辐角导数
Figure BDA0002866868880000095
计算当前时刻卫星运行轨道的近地点辐角
Figure BDA0002866868880000096
其中ωI为卫星入轨时刻的近地点辐角；
SA4.构建地心轨道坐标系到地心惯性坐标系的转换矩阵,其中i为轨道倾角
Figure BDA0002866868880000097
SA5.对地心轨道坐标系下的位置矢量进行坐标转转获取地心惯性ECI坐标系下的位置矢量peci
peci＝Τm×r
SA6.查询地球运动参数EOP文件,将ECI位置矢量peci转换为ECEF位置矢量pecef；
SA7.根据ECEF位置矢量pecef计算当前时刻圆形轨道卫星的星下点轨迹地心经度Lon。
进一步地，SA 6.查询地球运动参数EOP文件,将ECI位置矢量peci转换为ECEF位置矢量pecef，具体为：
SA601、从EOP文件中读取当前时刻t的6个关键地球运动参数:
1.(TAI-UTC)差值时间dat,单位秒,其中UTC为当前时刻t的世界协调时,原子时TAI＝UTC+dat
2.极移X分量xp，单位为弧度；
3.极移Y分量yp，单位为弧度；
4.(UTC1-UTC)差值时间dut,单位秒,其中UTC1为消除了极移影响后得到的世界时；
5.赤经章动ΔΨ修正量δΔΨ，单位为弧度；
6.交角章动Δε修正量δΔε，单位为弧度；
SA602、计算当前时刻t的岁差转换矩阵P(t)：
Figure BDA0002866868880000101
其中角度变量值ζ、θ和z为历元时刻T的平赤道面和平春分点相对于J2000平赤道和平春分点的3个角度值,其计算定义(公式中的角度单位为度)如下：
ζ＝(2306.2181T+0.30188*T2+0.017998T3)/3600.0
θ＝(2004.3109T-0.42665T2-0.041833T3)/3600.0
z＝(2306.2181T+1.09468T2+0.018203T3)/3600.0
定义历元时刻T指的是地球时从J2000的TT时刻起算的儒略世纪数,T＝(JDTT-2451545.0)/36525.0；其中地球时TT＝TAI+32.184,其对应的儒略日时间为JDTT；
SA603、计算当前时刻t的章动转换矩阵N(t)：
Figure BDA0002866868880000111
其中
Figure BDA0002866868880000112
Figure BDA0002866868880000113
其中章动角分量φi和系数Ai,Bi,Ci,Di根据IAU1980章动数据表进行计算；
其中历元T时刻平黄赤交角为
Figure BDA0002866868880000114
历元T时刻真黄赤交角
Figure BDA0002866868880000115
SA604、计算当前时刻t的地球自转转换矩阵R(t)：
Figure BDA0002866868880000116
其中格林尼治真恒星时
Figure BDA0002866868880000117
其中θGMST＝GMST(JDUT1)为当前世界时UTC1时刻儒略日时间的格林尼治平恒星时；
SA605、计算当前时刻t的极移转换矩阵W(t)：
Figure BDA0002866868880000121
其中极移分量xp和yp从EOP数据中获取；
SA606、t时刻ECI坐标下的坐标向量rECI转化为ECEF坐标下的坐标向量rECEF：rECEF＝[W(t)′][R(t)′][N(t)′][P(t)′]rECI。
有益效果：
本卫星轨道规划设计所要解决的问题就是要根据所选取的目标区域信息、轨道回归特性的约束条件、侦察载荷的约束条件、发射部署参数(包括发射点位置、发射部署完成时刻和发射轨道数据等参数)和精确的地球运动模型等数据，完成针对不同类型目标区域探测的太阳同步圆轨道的规划设计。本发明能够结合卫星发射轨道模型、轨道回归特性约束、目标区域光照条件约束和精确地球运动模型(本设计考虑了岁差、章动、自转和极移等地球复杂运动模型)进行太阳同步圆轨道方案设计，实现对目标区域的快速有效探测。
附图说明
图1为卫星发射和入轨过程示意图；
图2态势1：目标点地心纬度小于入轨点地心纬度，目标点在北半球；
图3态势2：目标点地心纬度小于入轨点地心纬度，目标点在南半球；
图4为态势3：目标点地心纬度大于入轨点地心纬度，目标点在北半球；
图5为态势4：目标点地心纬度大于入轨点地心纬度，目标点在南半球；
图6为态势5：目标点地心纬度大于入轨点地心纬度，入轨点在下降段，目标点在北半球；
图7为态势6：目标点地心纬度大于入轨点地心纬度，入轨点在下降段，目标点在南半球；
图8为态势7：目标点地心纬度大于入轨点地心纬度，入轨点在上升段，目标点在北半球；
图9为态势8(无效)：目标点地心纬度大于入轨点地心纬度，入轨点在上升段，目标点在南半球；
图10为态势9示意图：目标点地心纬度小于入轨点地心纬度，目标点在北半球；
图11为态势10：目标点地心纬度小于入轨点地心纬度，目标点在南半球；
图12太阳同步圆轨道设计的轨道倾角遍历象限图；
图13为本发明实施例中的太阳同步圆轨道设计OnComputeSunOrbit模块的流程图；
图14为本发明实施例中的遍历范围内太阳同步圆轨道星下点过境目标纬度圈地心经度分布计算OnComputeOrbitMinValueSunOrbit模块的流程图；
图15为本发明实施例中的太阳同步圆轨道卫星的发射轨道和运行轨道计算OnComputeOrbitValueSunOrbit模块的流程图；
图16为本发明实施例中的圆形轨道卫星的星下点轨迹计算OnComputeValueLonGCcircle模块流程图；
图17为本发明实施例中的ECI坐标系转ECEF坐标系OnConvertECItoECEF模块流程图；
图18为本发明实施例中的太阳同步圆轨道半长轴迭代计算OnComputeSunOrbitAa模块流程图；
图19为本发明实施例中的太阳同步圆轨道升交点赤经导数偏差计算OnComputeOmegaDotSunOrbit模块流程图。
具体实施方式
下面结合附图并举实施例，对本发明进行详细描述。
本卫星轨道规划设计所要解决的问题就是要根据所选取的目标区域信息、轨道回归特性的约束条件、侦察载荷的约束条件、发射部署参数(包括发射点位置、发射部署完成时刻和发射轨道数据等参数)和精确的地球运动模型等数据，完成针对不同类型目标区域探测的太阳同步圆轨道的规划设计。
根据卫星发射轨道参数和运行轨道参数，对上面轨道规划设计所需要解决的问题构建数学模型。
在理想情况下，卫星在发射和入轨过程受到的推力可近似为“两次冲量”的作用，这两次冲量分别位于发射点L和入轨点I，如图1所示，图中，Re为地球半径，βL为发射轨道的地心角，Θ为发射速度倾角，σ为近地快速覆盖轨道，O为地心。近地快速覆盖轨道可分为发射轨道与运行轨道两部分在发射点L和入轨点I之间的部分为卫星发射轨道；卫星经过入轨点I之后为卫星运行轨道。
近地快速覆盖轨道的发射入轨既要满足快速性的要求又要尽量节省能量。卫星发射方式包括共面发射和非共面发射，由于非共面发射比共面发射需要消耗更多的能量，近地快速覆盖轨道采用共面发射。若发射轨道在入轨点与运行轨道相切，该轨道称为“单共切”轨道。当发射速度倾角Θ给定时，单共切轨道为最佳发射轨道。近地快速覆盖轨道通常采用共面发射，单共切轨道入轨。在共面发射条件下，卫星发射轨道和运行轨道的轨道倾角i(定义为用地轴的北极方向与轨道平面的正法线方向之间的夹角度量，其中0≤i≤π，i＝0或π赤道轨道，i＝π/2为极地轨道)相等，发射轨道和运行轨道的升交点赤经Ω(定义为从北极上看从春分点逆时针转到升交点位置的转角，0≤Ω＜2π)也相等。发射轨道为椭圆轨道，通常可以采用卫星发射时刻tL的轨道6根数{aL,eL,i,Ω,ωL,fL}描述；卫星运行轨道为圆或椭圆轨道，可以采用卫星入轨时刻tI的轨道6根数{aI,eI,i,Ω,ωI,fI}描述。
其中，aL为发射轨道半长轴，eL为发射轨道偏心率，ωL为发射轨道近地点辐角，fL为发射轨道真近点角。
其中，aI为运行轨道半长轴，eI为运行轨道偏心率(0≤eI＜1)，ωI为运行轨道近地点辐角(沿着卫星运动方向从升交点转到近地点的角度，0≤ωI＜2π)，fI为运行轨道真近点角(沿着卫星运动方向从近地点转到卫星位置的角度，0≤fI＜2π)。设置发射点L的地心经纬度为
Figure BDA0002866868880000151
入轨点I的地心经纬度为
Figure BDA0002866868880000152
设置地面目标点的地心经纬度为
Figure BDA0002866868880000153
定义3个时刻点：发射时刻tL、入轨时刻tI和卫星过境目标纬度圈时刻tT。卫星纬度辐角定义为近地点辐角和真近点角之和。定义卫星发射时刻tL的卫星纬度辐角uL＝fL+ωL，卫星入轨时刻tI的卫星纬度辐角uI＝fI+ωI，卫星过境目标纬度圈时刻tT的卫星纬度辐角uT＝fT+ωI。在给定运行轨道半长轴aI和发射轨道地心角βL的条件下，可以由发射点和地面目标点的位置求解发射轨道和运行轨道的轨道根数。
卫星采用单共切轨道入轨，发射轨道与地球相交于发射点，与运行轨道相切于远地点，该点为卫星入轨点。发射点在发射轨道上的真近点角fL＝π-βL，入轨点在发射轨道上的真近点角fI＝π。定义发射轨道半通径计算如下：
p＝aI×(1+eL×cos(fI))
发射轨道偏心率eL的计算如下：
Figure BDA0002866868880000154
发射点的偏近点角EL的计算如下：
Figure BDA0002866868880000155
发射点的平近点角ML的计算如下：
ML＝EL-eL sin(EL)
发射轨道半长轴aL计算如下：
Figure BDA0002866868880000161
卫星入轨时刻tI的计算如下：
Figure BDA0002866868880000162
其中μ＝398600.4415km3/s2为引力常数。
对于地心经纬度为
Figure BDA0002866868880000163
的发射点L，根据球面三角公式，则发射点地心纬度
Figure BDA0002866868880000164
轨道倾角i和发射点卫星纬度辐角uL具有如下重要关系：
Figure BDA0002866868880000165
当前本轨道规划设计中，应急发射点只考虑在国内的情况(即在北半球)，则
Figure BDA0002866868880000166
所以可得0＜uL＜π。
根据卫星发射入轨的四种模式(1逆行下降，2逆行上升，3顺行下降，4顺行上升),可以分别计算发射点卫星纬度辐角uL。
对于卫星下降入轨模式1和模式3的发射点卫星纬度辐角uL计算如下：
Figure BDA0002866868880000167
对于卫星上升入轨模式2和模式4的发射点卫星纬度辐角uL计算如下：
Figure BDA0002866868880000168
定义发射时刻tL的格林尼治平恒星时GMST(瞬时平春分点和格林尼治子午圈之间的夹角)为GL，定义发射时刻tL在ECI地心惯性坐标系中地心经纬度为
Figure BDA0002866868880000169
的发射点L赤经αL计算如下：
αL＝OnLimitPiRadian(GL+λL)
其中角度约束函数OnLimitPiRadian(x)将角度x约束到[-π，π)范围内。对于地心经纬度为
Figure BDA0002866868880000171
的发射点L，根据球面三角公式，则发射点赤经αL、轨道倾角i、发射点卫星纬度辐角uL和卫星轨道升交点赤经Ω具有如下重要关系：
tan(αL-Ω)＝tan(uL)cos(i)
对于卫星下降入轨模式1和模式3的卫星轨道升交点赤经Ω计算如下：
Ω＝OnLimit2PiRadian(αL-tan-1(tan(uL)*cos(i))+π)
对于卫星上升入轨模式2和模式4的卫星轨道升交点赤经Ω计算如下：
Ω＝OnLimit2PiRadian(αL-tan-1(tan(uL)*cos(i)))
其中角度约束函数OnLimit2PiRadian(x)将角度x约束到[0，2π)范围内。
入轨时刻tI的卫星纬度辐角uI＝uL+βL。对于地心经纬度为
Figure BDA0002866868880000172
的入轨点I，根据球面三角公式，则入轨点地心纬度
Figure BDA0002866868880000173
轨道倾角i和入轨点卫星纬度辐角uI具有如下重要关系：
Figure BDA0002866868880000174
如果入轨点在北半球，则卫星纬度辐角uI＜π
如果入轨点在南半球，则卫星纬度辐角uI＞π
不管入轨点是在北半球还是在南半球，入轨点的地心纬度计算如下所示：
Figure BDA0002866868880000175
对于卫星下降入轨模式1和模式3，分4种态势计算卫星入轨后第1次和第2次过境目标纬度圈的卫星纬度辐角值uT1和uT2，这4种态势的示意图如图2～图5所示。
对于卫星上升入轨模式2和模式4，分如下6种态势计算卫星过境目标纬度圈的卫星纬度辐角值uT1和uT2，这6种态势(其中态势8为无效态势)的示意图如图6～图11所示：
获取了卫星入轨后星下点第1次和第2次过境目标纬度圈的卫星纬度辐角值uT1和uT2,可以根据圆形轨道和椭圆轨道的卫星运动方程快速计算后续卫星星下点过境目标纬度圈的卫星纬度辐角值(具体计算参考后续轨道规划问题的求解算法研究部分)。本发明根据卫星星下点过境目标纬度圈的纬度辐角值和过境时刻，可以根据卫星轨道方程计算出卫星过境目标纬度圈时刻在ECI惯性坐标系中的卫星空间位置矢量。通过空间坐标变换，本发明可以获取卫星过境目标纬度圈时刻在ECEF惯性坐标系中的卫星空间位置矢量，即可以获取过境目标纬度圈时刻卫星星下点的大地经度。
为了获取精确的卫星星下点大地经度，需要在ECI惯性坐标系空间位置矢量转ECEF惯性坐标系空间位置矢量的计算过程中，就需要考虑岁差(precession)、章动(nutation)、自转和极移(polar motion)等地球复杂运动模型的影响。只有考虑了这些地球复杂运动模型因素，卫星轨道规划计算结果才可以取得和STK软件仿真结果接近的精度。从1962年1月1号开始到2019年4月20号的地球运动模型参数都可以从STK软件提供的EOP文件中获取,该文件中的数据会定时更新。由于该EOP文件中的数据对于高精度坐标系转换计算非常重要，但在公开的文献中很少介绍这些参数的使用。在卫星轨道规划软件研发过程中，通过查阅文献和开源软件，逐渐掌握了EOP文件中数据的使用方法。
岁差、章动和极移都是引起天极(地球自转轴在天球上的投影)运动的原因，其中岁差与章动引起天极在空间的运动，极移引起其在地球本体内的运动。t时刻ECI坐标下的空间坐标向量rECI与ECEF坐标下的空间坐标向量rECEF之间的关系如下:
rECI＝[P(t)][N(t)][R(t)][W(t)]rECEF；
其中t为当前时刻,P(t)为岁差(precession)转换矩阵,N(t)为章动(nutation)转换矩阵，R(t)为与当前时刻格林尼治真恒星时角(GAST)相关的地球自转矩阵，W(t)为极移(polar motion)转换矩阵。
定义旋转矩阵算子如下:
定义旋转矩阵算子如下:
Figure BDA0002866868880000191
Figure BDA0002866868880000192
Figure BDA0002866868880000193
定义原子时TAI＝UTC+dat,其中UTC为当前时刻t的世界协调时,dat为当前时刻对应的(TAI-UTC)差值,该差值数据可以从EOP文件中获取,单位为秒。定义世界时UTC1＝UTC+dut,dut为当前时刻对应的(UTC1-UTC)差值,该差值数据可以从EOP文件中获取,单位为秒。定义地球时TT＝TAI+32.184,地球时TT对应的儒略日时间为JDTT。
定义时间自变量T＝(JDTT-2451545.0)/36525.0指的是地球时TT从J2000时刻起算的儒略世纪数。则4个关键转换矩阵计算如下所示:
(1).岁差转换矩阵P(t)的计算
岁差描述了地球自转轴指向和春分点的长期变化。岁差转换矩阵P(t)的计算如下:
Figure BDA0002866868880000194
其中角度变量值ζ,θ和z为历元时刻T的平赤道面和平春分点相对于J2000平赤道和平春分点的3个角度值,其计算定义(公式中的角度单位为度)如下:
ζ＝(2306.2181T+0.30188*T2+0.017998T3)/3600.0
θ＝(2004.3109T-0.42665T2-0.041833T3)/3600.0
z＝(2306.2181T+1.09468T2+0.018203T3)/3600.0
(2).章动转换矩阵N(t)计算
章动描述了赤道和春分点的短期和周期性变化。除了长期岁差运动，地球自转轴指向还存在小的周期性扰动即章动。这是由月球和太阳扭矩的月变化和年变化引起的。章动影响包括赤经章动ΔΨ(nutation in longitude)和交角章动Δε(nutation inobliquity)。
章动转换矩阵N(t)的计算如下:
Figure BDA0002866868880000201
赤经章动ΔΨ计算如下:
Figure BDA0002866868880000202
交角章动Δε计算如下:
Figure BDA0002866868880000203
其中δΔΨ为赤经章动ΔΨ修正量,δΔε为交角章动Δε修正量，可以从EOP文件中获取当前时刻对应的δΔΨ和δΔε值。
其中章动角分量φi＝pL,iL+pM,iM+pE,iE+pF,iF+pG,iG
其中整数系数pL,i,pM,i,pE,i,pF,i,pG,i和小数系数Ai,Bi,Ci,Di可以查IAU1980章动数据表(参考:卫星轨道:模型方法和应用)获取。其中L为月球平近点角，M为太阳平近点角，E为月球升交距角，F为日月间平经度差，G为月球轨道升交点平经度,这些中间变量(计算的角度单位为度)计算如下:
L＝134.96298139+(1717915922.6330T+31.310T2+0.064T3)/3600.0
M＝357.52772333+(129596581.2240T-0.577T2-0.012T3)/3600.0
E＝93.27191028+(1739527263.1370T-13.257T2+0.011T3)/3600.0
F＝297.85036306+(1602961601.3280T-6.891T2+0.019T3)/3600.0
G＝125.04452222+(-6962890.5390T+7.455T2+0.008T3)/3600.0
定义
Figure BDA0002866868880000211
为时刻t真黄(ecliptic)赤(equator)交角,其中
Figure BDA0002866868880000212
为时刻t平黄赤交角,Δε为交角章动值。
平黄赤交角
Figure BDA0002866868880000213
(计算的角度单位为度)的计算如下:
Figure BDA0002866868880000214
(3).地球自转转换矩阵R(t)计算
自转转换矩阵R(t)是由格林尼治真恒星时θGAST的影响造成，其计算如下:
Figure BDA0002866868880000215
其中格林尼治真恒星时θGAST(单位为弧度)的计算如下:
Figure BDA0002866868880000216
其中θGMST＝GMST(JDUT1)为当前世界时UTC1时刻的格林尼治平恒星时,JDUT1为世界时UTC1时刻对应的儒略日时间，函数GMST(JD)用于计算儒略日时间JD的格林尼治平恒星时。
其中赤经章动ΔΨ和交角章动Δε的计算参见前面公式。
(4).极移转换矩阵W(t)计算
由于地球不是刚体以及其他一些物理因素的影响，造成地球自转轴在地球体内的运动，引起地球极点在地球表面的位置随时间而变化，称为极移。
极移转换矩阵W(t)计算如下:
Figure BDA0002866868880000217
其中xp为极移X分量,yp为极移Y分量，可以从EOP文件中获取当前时刻对应的xp和yp值。
本发明实施例提供的太阳同步圆轨道设计针对有效倾角范围内的角度值以特定间隔进行遍历，计算每个倾角角度值i对应卫星轨道星下点轨迹在多次过境目标点纬度圈上的经度值，并且计算目标点地心经度和星下点轨迹过境标点纬度圈时地心经度值之间的差值。选取地心经度差值满足指标要求的卫星轨道作为备选方案，最后以优选准则对备选轨道方案进行选取。由于太阳同步圆轨道卫星的轨道倾角大于90度，为逆行轨道。所以太阳同步轨道只有两种入轨模式:1逆行下降,2逆行上升。此外在卫星轨道设计流程中考虑了太阳光照参数对于光学载荷卫星轨道设计的实施要求,并根据轨道参数计算在太阳高度角约束条件下的卫星发射窗口。
对于任意给定的发射点，以发射点为中心，以经纬度线为坐标轴构建四个象限，如图12所示。
对于太阳同步圆轨道只能在象限1和象限2这两个象限内根据输入的卫星轨道高度上下限计算选取一定的角度范围作为卫星轨道倾角的遍历范围。由于太阳同步圆轨道的升交点赤经导数为0.985612288度/天，轨道半长轴与轨道倾角必须满足如下条件：
Figure BDA0002866868880000221
其中轨道半长轴aI、轨道倾角i、地球平均半径Re、地球自转速度ωE和卫星平均角速度n。太阳同步圆轨道设计中，本发明对于2种发射入轨模式遍历范围内的每条卫星轨道，计算卫星入轨后16次星下点过境目标纬度圈时刻的ECEF地心地固坐标系下大地经度值。
当前卫星运行轨道为圆形轨道，选取卫星轨道升交点与轨道近地点重合，则运行轨道的近地点辐角ωI＝0。根据J2摄动条件下的卫星圆形运行轨道方程，可以计算出t1、t2到t16时刻16次卫星星下点过境目标纬度圈时刻在ECEF地心地固坐标系下的大地经度LonT(1)、LonT(2)到LonT(16)。
太阳同步圆轨道设计OnComputeSunOrbit模块对于当前2种入轨模式下的轨道倾角遍历范围内所有卫星轨道，调用太阳同步圆轨道卫星的发射轨道和运行轨道计算OnComputeOrbitValueSunOrbit模块计算卫星入轨后星下点16次过境目标纬度圈时的经度分布，选取遍历范围内过境时刻星下点经度与目标点经度偏差值最小的轨道进行门限判决，并保存有效轨道的计算参数到BackupOrbit数组。在卫星轨道设计流程中还需要考虑了太阳光照参数对于光学载荷卫星轨道设计的实施要求,并根据轨道参数计算在太阳高度角约束条件下的卫星发射窗口。
太阳同步圆轨道设计OnComputeSunOrbit模块流程如图13所示：包括如下步骤：
步骤1、太阳同步圆轨道具有两种入轨模式：即逆行下降和逆行上升2种入轨模式；根据轨道半长轴的上下限，调用OnComputeSunOrbitAngle模块计算获得逆行上升和逆行下降的轨道倾角搜索范围[domainDegreeLowerW,domainDegreeUpperW]。
如果卫星载荷类型为CCD光学设备，则以发射部署完成时刻为起点，计算目标点区域24小时内太阳高度角满足门限要求的UTC时间观测窗口为[m_SunStartTime,m_SunEndTime]。
该轨道设计针对逆行上升和逆行下降轨道的倾角搜索范围内所有卫星轨道，计算卫星轨道星下点轨迹16次过境目标纬度圈时星下点的经度分布，并且选取满足侦察要求的最优轨道；设定初始状态为，MO为入轨模式，MO的取值为1或2，分别指代逆行下降和逆行上升2种入轨模式，MO初始值设定为1；j为过境次数，j的初值设定为1；BackupOrbitNum为轨道方案个数，BackupOrbitNum的初值设定为0；
步骤2、判断j的值是否满足j<＝16，若是则进入步骤3，否则进入步骤11；
步骤3、判断MO的值是否满足MO<＝2；若是则进入步骤4，否则将MO重新设定为1，j自增1，返回步骤2；
步骤4、对于第MO种入轨模式轨道倾角遍历范围内所有卫星轨道，计算第j次过境目标纬度圈时星下点经度分布,选取轨道倾角遍历范围内过境时刻星下点经度与目标点经度偏差值最小的轨道；
步骤5、判断步骤4中所选取的偏差值最小的轨道对应的偏差值是否满足预设的指标要求，若是则进入步骤6，否则MO自增1，返回步骤3；
步骤6、判断当前卫星载荷类型是否为CCD光学设备，若是则进入步骤7，否则进入步骤10；
步骤7、判断当前轨道卫星过境目标点时刻是否在CCD光学设备的UTC时间观测窗口内，若是则进入步骤10，否则进入步骤8；
步骤8、计算发射延时，对于当前入轨模式轨道倾角遍历范围内的所有卫星轨道，计算当前第j次过境目标纬度圈时星下点经度分布,选取遍历范围内过境时刻星下点经度与目标点经度偏差值最小的轨道；
步骤9、判断步骤8中所选取的偏差值最小的轨道是否满足预设的指标要求，若是则进入步骤10，否则返回步骤3；
步骤10、当前第BackupOrbitNum个轨道规划方案满足初步要求，将当前第BackupOrbitNum个轨道规划方案存入在BackupOrbit数组中,并且BackupOrbitNum自增1；
步骤11、太阳同步圆轨道规划完成，将当前对应的共BackupOrbitNum个轨道方案存入预先构建的BackupOrbit数组中。
其中OnComputeSunOrbitAngle模块计算太阳同步圆轨道倾角函数。
步骤4中，对于第MO种入轨模式轨道倾角遍历范围内所有卫星轨道，计算第j次过境目标纬度圈时星下点经度分布,选取轨道倾角遍历范围内过境时刻星下点经度与目标点经度偏差值最小的轨道，采用遍历范围内太阳同步圆轨道星下点过境目标纬度圈地心经度分布计算OnComputeOrbitMinValueSunOrbit模块计算当前倾角遍历范围内所有太阳同步圆轨道计算卫星星下点轨迹过境目标纬度圈时的经度分布。遍历范围内太阳同步圆轨道星下点过境目标纬度圈地心经度分布计算OnComputeOrbitMinValueSunOrbit模块的流程如图14所示。包括如下步骤：
步骤101、启动遍历范围内太阳同步圆轨道星下点过境目标纬度圈地心经度分布计算OnComputeOrbitMinValueSunOrbit模块；
步骤102、轨道倾角遍历范围的参数初始化：轨道倾角遍历初始值domainDegreeLower,遍历间隔0.01度,遍历次数indexNum,索引index＝1；
步骤103、判断是否满足index<＝indexNum；若是则进入步骤104；否则进入步骤105；
步骤104、对当前轨道倾角为IiAngleDegree＝domainDegreeLower+0.01×(index-1)的卫星发射轨道和卫星运行轨道，按照如下S1041～S1043计算之后进入步骤105；
S1041、调用太阳同步圆轨道半长轴迭代计算OnComputeSunOrbitAa模块计算对应当前轨道倾角倾角IiAngleDegree的太阳同步圆轨道半长轴；
S1042、根据当前轨道倾角IiAngleDegree和太阳同步圆轨道半长轴，调用太阳同步圆轨道卫星的发射轨道和运行轨道计算OnComputeOrbitValueSunOrbit模块进行太阳同步圆轨道星下点轨迹计算，得到卫星星下点轨迹多次过境目标纬度圈时的星下点地心经度与目标点地心经度的偏差数据；
S1043、对索引值index自增1，返回步骤103；
步骤105、对于遍历范围内indexNum条太阳同步圆轨道星下点轨迹第j次(j从1到16,为模块输入参数)过境目标纬度圈时刻的星下点地心经度与目标点地心经度的偏差分布进行统计分析处理，选取偏差值最小的轨道方案,最小偏差值为minValue；
步骤106、对最小偏差值minValue进行门限threshold判决，若minValue<＝threshold，则进入步骤107；否则返回无效值；
步骤107、minValue对应的轨道方案保存在minTimeOrbit数组中，并返回有效轨道倾角值。
S1042中，根据当前轨道倾角IiAngleDegree和太阳同步圆轨道半长轴，调用太阳同步圆轨道卫星的发射轨道和运行轨道计算OnComputeOrbitValueSunOrbit模块进行太阳同步圆轨道星下点轨迹计算，得到卫星星下点轨迹多次过境目标纬度圈时的星下点地心经度与目标点地心经度的偏差数据；太阳同步圆轨道卫星的发射轨道和运行轨道计算OnComputeOrbitValueSunOrbit核心计算模块流程图如图15所示，包括如下步骤：
S10421、卫星发射轨道参数初始化计算部分：
1.发射点地心经纬度为
Figure BDA0002866868880000261
发射轨道地心角βL和发射轨道倾角i为模块输入参数；
2.发射时刻卫星在发射轨道上真近点角fLradian和入轨时刻卫星在发射轨道上真近点角fIradian的计算；
3.发射轨道偏心率eL、发射轨道半通径p和发射轨道半长轴aL的计算
4.发射点偏近点角ELangleRadian和发射点平近点角MLangleRadian的计算
5.根据上述参数计算卫星入轨时刻tI(相对于发射准备完成时刻)的计算
6.发射时刻GMST(瞬时平春分点和格林尼治子午圈之间的夹角)格林尼治平恒星时GL的计算
7.发射时刻发射点在ECI地心惯性坐标系中赤经为：
αL＝OnLimitPiRadian(GL+λL)，其中角度约束函数OnLimitPiRadian(x)将角度x约束到[-π,π)范围内。
S10422、若为下降入轨模式，发射时刻卫星纬度辐角计算如下:
Figure BDA0002866868880000271
轨道升交点赤经:Ω＝OnLimit2PiRadian(αL-tan-1(tan(uL)×cos(i))+π)。
若为上升入轨模式，发射时刻卫星纬度辐角计算如下:
Figure BDA0002866868880000272
轨道升交点赤经:Ω＝OnLimit2PiRadian(αL-tan-1(tan(uL)×cos(i))+0)；
S10423、卫星运行轨道参数初始化部分:
1.入轨时刻卫星纬度辐角uI＝uL+βL,目标点地心经纬度
Figure BDA0002866868880000273
卫星平均角速度n,地球自转转速ωE；
2.入轨点地心纬度为
Figure BDA0002866868880000274
S10424、若为下降入轨模式，则分如下4种情况计算卫星入轨后第1次和第2次过境目标纬度圈的卫星纬度辐角值uT1和uT2；目标点地心纬度为
Figure BDA0002866868880000275
入轨点地心维度
Figure BDA0002866868880000276
第1种情况.目标点在北半球,目标点地心纬度小于等于入轨点地心纬度；
第2种情况.目标点在南半球,目标点地心纬度小于等于入轨点地心纬度；
第3种情况.目标点在北半球,目标点地心纬度大于入轨点地心纬度；
第4种情况.目标点在南半球,目标点地心纬度大于入轨点地心纬度；
若
Figure BDA0002866868880000277
即目标点地心纬度小于等于入轨点地心纬度，则对于入轨后第1次过境，
Figure BDA0002866868880000278
包括第1种情况和第2种情况，如果OnLimit2PiRadian(uT1)<π则属于第1种情况，第一中间变量为Δu＝π-uT1和第二中间变量为utemp＝π+2Δu；否则属于第2种情况：Δu＝uT1-π；utemp＝π-2Δu；则入轨后第2次过境uT2＝uT1+utemp；
若
Figure BDA0002866868880000281
即目标点地心纬度大于入轨点地心纬度，则对于入轨后第1次过境，
Figure BDA0002866868880000282
包括第3种情况和第4种情况：
如果OnLimit2PiRadian(uT1)<π则属第3种情况：以tempRadian指代OnLimit2PiRadian(uT1)，Δu＝tempRadian；utemp＝π-2Δu；否则属于第4种情况：Δu＝fabs(tempRadian-2π)；utemp＝π+2Δu；则入轨后第2次过境uT2＝uT1+utemp；
若为上升入轨模式，分如下6种情况计算卫星入轨后第1次和第2次过境目标纬度圈的卫星纬度辐角值uT1和uT2；
第5种情况.目标点在北半球,目标点纬度大于入轨点纬度,入轨点在下降段；
第6种情况.目标点在南半球,目标点纬度大于入轨点纬度,入轨点在下降段；
第7种情况.目标点在北半球,目标点纬度大于入轨点纬度,入轨点在上升段；
第8种情况.目标点在南半球,目标点纬度大于入轨点纬度,入轨点在上升段；
第9种情况.目标点在北半球,目标点纬度小于入轨点纬度；
第10种情况.目标点在南半球,目标点纬度小于入轨点纬度；
若
Figure BDA0002866868880000283
即目标点地心纬度大于等于入轨点地心纬度，且入轨点在下降段
Figure BDA0002866868880000284
则对于入轨后第1次过境，
Figure BDA0002866868880000285
tempRadian＝OnLimit2PiRadian(uT1)；
若tempRadian＜π属于第5种情况，Δu＝tempRadian；utemp＝π-2Δu；否则属于第6种情况，则Δu＝fabs(tempRadian-2π)；utemp＝π+2Δu；
若
Figure BDA0002866868880000286
即目标点地心纬度大于等于入轨点地心纬度，且入轨点在上升段
Figure BDA0002866868880000287
则对于入轨后第1次过境，
Figure BDA0002866868880000288
tempRadian＝OnLimit2PiRadian(uT1)；
若tempRadian＜π属于第7种情况，Δu＝uT1；utemp＝π-2Δu；否则属于第8种情况，第8种情况不存在；
则入轨后第2次过境uT2＝uT1+utemp；
若
Figure BDA0002866868880000291
即目标点地心纬度小于入轨点地心纬度，则对于入轨后第1次过境，
Figure BDA0002866868880000292
tempRadian＝OnLimit2PiRadian(uT1)；
若tempRadian＜π(即目标点在北半球)属于第9种情况，
Figure BDA0002866868880000293
utemp＝π+2Δu；否则属于第10种情况，
Figure BDA0002866868880000294
utemp＝π+2Δu；
则入轨后第2次过境uT2＝uT1+utemp；
S10425、卫星运行轨道星下点过境目标纬度圈时刻地心经度计算:
1.计算卫星运行轨道周期tΩ、升交点赤经导数
Figure BDA0002866868880000295
和近地点辐角导数
Figure BDA0002866868880000296
2.计算卫星入轨后16次过境目标纬度圈的时刻:
t1＝(uT1-uI)/n,t2＝(uT2-uI)/n,t2m+1＝t1+mtΩ,t2m+2＝t2+mtΩ,m＝1,2,...,7
3.调用圆形轨道卫星的星下点轨迹计算OnComputeValueLonGCcircle模块计算入轨后第1次和第2次过境目标纬度圈时星下点地心经度LonT(1)和LonT(2)；
4.计算入轨后第3次到第16次卫星星下点轨迹过境目标纬度圈时刻的地心经度:
Figure BDA0002866868880000297
Figure BDA0002866868880000298
5.计算卫星星下点轨迹16次过境目标纬度圈时刻地心经度与目标地心经度的偏差。
S10426、将计算结果和中间变量存放到OrbitVector数组
计算入轨后第1次和第2次过境目标纬度圈时星下点地心经度LonT(1)和LonT(2)；采用圆形轨道卫星的星下点轨迹计算OnComputeValueLonGCcircle模块，其流程如图16所示:
SA1.卫星半长轴aI和当前时刻纬度辐角uT,卫星在地心轨道坐标系下的位置矢量r＝[aI cos(uT) aI sin(uT) 0]′；
SA2.根据相对入轨时刻时间差Δt和升交点赤经导数
Figure BDA0002866868880000301
计算当前时刻卫星运行轨道的升交点赤经
Figure BDA0002866868880000302
其中ΩI为卫星入轨时刻的升交点赤经；
SA3.根据相对入轨时刻时间差Δt和近地点辐角导数
Figure BDA0002866868880000303
计算当前时刻卫星运行轨道的近地点辐角
Figure BDA0002866868880000304
其中ωI为卫星入轨时刻的近地点辐角；
SA4.构建地心轨道坐标系到地心惯性坐标系的转换矩阵,其中i为轨道倾角
Figure BDA0002866868880000305
SA5.对地心轨道坐标系下的位置矢量进行坐标转转获取地心惯性ECI坐标系下的位置矢量peci
peci＝Τm×r
SA6.查询地球运动参数EOP文件,调用ECI坐标系转ECEF坐标系OnConvertECItoECEF模块将ECI位置矢量peci转换为ECEF位置矢量pecef；
SA7.根据ECEF位置矢量pecef计算当前时刻圆形轨道卫星的星下点轨迹地心经度Lon。
SA6.查询地球运动参数EOP文件,调用ECI坐标系转ECEF坐标系OnConvertECItoECEF模块将ECI位置矢量peci转换为ECEF位置矢量pecef，ECI坐标系转ECEF坐标系OnConvertECItoECEF模块流程图如图17所示，包括如下步骤：
SA601、从EOP文件中读取当前时刻t的6个关键地球运动参数:
1.(TAI-UTC)差值时间dat,单位秒,其中UTC为当前时刻t的世界协调时,原子时TAI＝UTC+dat
2.极移X分量xp，单位为弧度；
3.极移Y分量yp，单位为弧度；
4.(UTC1-UTC)差值时间dut,单位秒,其中UTC1为消除了极移影响后得到的世界时；
5.赤经章动ΔΨ修正量δΔΨ，单位为弧度；
6.交角章动Δε修正量δΔε，单位为弧度；
SA602、计算当前时刻t的岁差转换矩阵P(t)：
Figure BDA0002866868880000311
其中角度变量值ζ、θ和z为历元时刻T的平赤道面和平春分点相对于J2000平赤道和平春分点的3个角度值,其计算定义(公式中的角度单位为度)如下：
ζ＝(2306.2181T+0.30188*T2+0.017998T3)/3600.0
θ＝(2004.3109T-0.42665T2-0.041833T3)/3600.0
z＝(2306.2181T+1.09468T2+0.018203T3)/3600.0
定义历元时刻T指的是地球时从J2000的TT时刻起算的儒略世纪数,T＝(JDTT-2451545.0)/36525.0；其中地球时TT＝TAI+32.184,其对应的儒略日时间为JDTT；
SA603、计算当前时刻t的章动转换矩阵N(t)：
Figure BDA0002866868880000312
其中
Figure BDA0002866868880000321
Figure BDA0002866868880000322
其中章动角分量φi和系数Ai,Bi,Ci,Di根据IAU1980章动数据表进行计算；
其中历元T时刻平黄赤交角为
Figure BDA0002866868880000323
历元T时刻真黄赤交角
Figure BDA0002866868880000324
SA604、计算当前时刻t的地球自转转换矩阵R(t)：
Figure BDA0002866868880000325
其中格林尼治真恒星时
Figure BDA0002866868880000326
其中θGMST＝GMST(JDUT1)为当前世界时UTC1时刻儒略日时间的格林尼治平恒星时；
SA605、计算当前时刻t的极移转换矩阵W(t)：
Figure BDA0002866868880000327
其中极移分量xp和yp从EOP数据中获取；
SA606、t时刻ECI坐标下的坐标向量rECI转化为ECEF坐标下的坐标向量rECEF：rECEF＝[W(t)′][R(t)′][N(t)′][P(t)′]rECI。
S1041、调用太阳同步圆轨道半长轴迭代计算OnComputeSunOrbitAa模块计算对应当前轨道倾角倾角IiAngleDegree的太阳同步圆轨道半长轴，太阳同步圆轨道半长轴迭代计算OnComputeSunOrbitAa模块流程图如图18所示，包括如下步骤：
S10411、根据从界面输入的轨道入轨点高度上限和下限分别计算太阳同步圆轨道半长轴上限aU和下限aL；根据当前模块的输入参数轨道倾角iiAngleDegreeTemp和轨道半长轴上下限，调用太阳同步圆轨道升交点赤经导数偏差计算OnComputeOmegaDotSunOrbit模块分别计算当前轨道升交点赤经导数与太阳同步轨道标准升交点赤经导数的差值tempValueU和tempValueL
S10412、判断若(tempValueL×tempValueU)>0.0，则输出无效值，否则进入S10413；
S10413、轨道半长轴迭代计算初始化:index＝0,差值tempValueL和tempValueU
S10414、判断index<30是否成立，若是进入S10415，否则输出轨道半长轴中间值aM；
S10415、进行轨道半长轴迭代计算：首先计算轨道半长轴中间值aM＝(aU+aL)/2.0；然后根据当前轨道倾角iiAngleDegreeTemp和轨道半长轴中间值aM，调用太阳同步圆轨道升交点赤经导数偏差计算OnComputeOmegaDotSunOrbit模块计算当前轨道升交点赤经导数与太阳同步轨道标准升交点赤经导数的差值tempValue；index自增1；
S10416、判断(tempValue×tempValueL)>0.0是否成立，若是则令tempValueL＝tempValue；同时令aL＝aM，否则令tempValueU＝tempValue，同时令aU＝aM；返回S10414。
太阳同步圆轨道升交点赤经导数偏差计算OnComputeOmegaDotSunOrbit模块流程图如图19所示：具体为：
根据轨道半长轴aI、轨道倾角i、地球平均半径Re、地球自转速度ωE和卫星平均角速度n，计算当前轨道升交点赤经导数与太阳同步圆轨道标准升交点赤经导数的偏差值：
Figure BDA0002866868880000341
对于本发明卫星轨道规划软件计算获得的卫星轨道方案，本发明将卫星轨道方案参数和目标区域参数带入到STK软件中进行仿真评估。根据STK软件的仿真评估结果，可以判断本发明卫星轨道规划设计和轨道评估计算的准确性。
综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (7)

1.一种太阳同步圆轨道的轨道规划方法，其特征在于，包括如下步骤：
步骤1、所述太阳同步圆轨道具有两种入轨模式：即逆行下降和逆行上升2种入轨模式；根据轨道半长轴的上下限，获得逆行上升和逆行下降的轨道倾角搜索范围；
如果卫星载荷类型为CCD光学设备，则以发射部署完成时刻为起点，计算目标点区域24小时内太阳高度角满足门限要求的UTC时间观测窗口；
该轨道规划方法针对逆行上升和逆行下降轨道的倾角搜索范围内所有卫星轨道，计算卫星轨道星下点轨迹16次过境目标纬度圈时星下点的经度分布，并且选取满足侦察要求的最优轨道；设定初始状态为，MO为入轨模式，MO的取值为1或2，分别指代逆行下降和逆行上升2种入轨模式，MO初始值设定为1；j为过境次数，j的初值设定为1；BackupOrbitNum为轨道方案个数，BackupOrbitNum的初值设定为0；
步骤2、判断j的值是否满足j<＝16，若是则进入步骤3，否则进入步骤11；
步骤3、判断MO的值是否满足MO<＝2；若是则进入步骤4，否则将MO重新设定为1，j自增1，返回步骤2；
步骤4、对于第MO种入轨模式轨道倾角遍历范围内所有卫星轨道，计算第j次过境目标纬度圈时星下点经度分布,选取轨道倾角遍历范围内过境时刻星下点经度与目标点经度偏差值最小的轨道；
步骤5、判断步骤4中所选取的偏差值最小的轨道对应的偏差值是否满足预设的指标要求，若是则进入步骤6，否则MO自增1，返回步骤3；
步骤6、判断当前卫星载荷类型是否为CCD光学设备，若是则进入步骤7，否则进入步骤10；
步骤7、判断当前轨道卫星过境目标点时刻是否在CCD光学设备的UTC时间观测窗口内，若是则进入步骤10，否则进入步骤8；
步骤8、计算发射延时，对于当前入轨模式轨道倾角遍历范围内的所有卫星轨道，计算当前第j次过境目标纬度圈时星下点经度分布,选取遍历范围内过境时刻星下点经度与目标点经度偏差值最小的轨道；
步骤9、判断步骤8中所选取的偏差值最小的轨道是否满足预设的指标要求，若是则进入步骤10，否则返回步骤3；
步骤10、当前第BackupOrbitNum个轨道规划方案满足初步要求，将当前第BackupOrbitNum个轨道规划方案存入在BackupOrbit数组中,并且BackupOrbitNum自增1；
步骤11、太阳同步圆轨道规划完成，将当前对应的共BackupOrbitNum个轨道方案存入预先构建的BackupOrbit数组中。
2.如权利要求1所述的方法，其特征在于，所述步骤4中，对于第MO种入轨模式轨道倾角遍历范围内所有卫星轨道，计算第j次过境目标纬度圈时星下点经度分布,选取轨道倾角遍历范围内过境时刻星下点经度与目标点经度偏差值最小的轨道，具体包括如下步骤：
步骤101、遍历范围内太阳同步圆轨道星下点过境目标纬度圈地心经度分布；
步骤102、轨道倾角遍历范围的参数初始化：轨道倾角遍历初始值domainDegreeLower,遍历间隔0.01度,遍历次数indexNum,索引index＝1；
步骤103、判断是否满足index<＝indexNum；若是则进入步骤104；否则进入步骤105；
步骤104、对当前轨道倾角为IiAngleDegree＝domainDegreeLower+0.01×(index-1)的卫星发射轨道和卫星运行轨道，按照如下S1041～S1043计算之后进入步骤105；
S1041、计算对应当前轨道倾角倾角IiAngleDegree的太阳同步圆轨道半长轴；
S1042、根据当前轨道倾角IiAngleDegree和太阳同步圆轨道半长轴，进行太阳同步圆轨道星下点轨迹计算，得到卫星星下点轨迹多次过境目标纬度圈时的星下点地心经度与目标点地心经度的偏差数据；
S1043、对索引值index自增1，返回步骤103；
步骤105、对于遍历范围内indexNum条太阳同步圆轨道星下点轨迹第j次过境目标纬度圈时刻的星下点地心经度与目标点地心经度的偏差分布进行统计分析处理，选取偏差值最小的轨道方案,最小偏差值为minValue；
步骤106、对最小偏差值minValue进行门限threshold判决，若minValue<＝threshold，则进入步骤107；否则返回无效值；
步骤107、minValue对应的轨道方案保存在minTimeOrbit数组中，并返回有效轨道倾角值。
3.如权利要求2所述的方法，其特征在于，所述S1041、计算对应当前轨道倾角倾角IiAngleDegree的太阳同步圆轨道半长轴，具体为如下步骤：
S10411、根据从界面输入的轨道入轨点高度上限和下限分别计算太阳同步圆轨道半长轴上限aU和下限aL；根据当前模块的输入参数轨道倾角iiAngleDegreeTemp和轨道半长轴上下限，分别计算当前轨道升交点赤经导数与太阳同步轨道标准升交点赤经导数的差值tempValueU和tempValueL
S10412、判断若(tempValueL×tempValueU)>0.0，则输出无效值，否则进入S10413；
S10413、轨道半长轴迭代计算初始化:index＝0,差值tempValueL和tempValueU
S10414、判断index<30是否成立，若是进入S10415，否则输出轨道半长轴中间值aM；
S10415、进行轨道半长轴迭代计算：首先计算轨道半长轴中间值aM＝(aU+aL)/2.0；然后根据当前轨道倾角iiAngleDegreeTemp和轨道半长轴中间值aM，计算当前轨道升交点赤经导数与太阳同步轨道标准升交点赤经导数的差值tempValue；index自增1；
S10416、判断(tempValue×tempValueL)>0.0是否成立，若是则令tempValueL＝tempValue；同时令aL＝aM，否则令tempValueU＝tempValue，同时令aU＝aM；返回S10414。
4.如权利要求2所述的方法，其特征在于，所述S10411和S10415中，计算当前轨道升交点赤经导数与太阳同步轨道标准升交点赤经导数的差值，具体为：
根据轨道半长轴aI、轨道倾角i、地球平均半径Re、地球自转速度ωE和卫星平均角速度n，计算当前轨道升交点赤经导数与太阳同步圆轨道标准升交点赤经导数的偏差值：
Figure FDA0002866868870000041
5.如权利要求2所述的方法，其特征在于，所述S1042中，根据当前轨道倾角IiAngleDegree和太阳同步圆轨道半长轴，进行太阳同步圆轨道星下点轨迹计算，得到卫星星下点轨迹多次过境目标纬度圈时的星下点地心经度与目标点地心经度的偏差数据；具体为：
S10421、卫星发射轨道参数初始化计算部分：
1.发射点地心经纬度为
Figure FDA0002866868870000042
发射轨道地心角βL和发射轨道倾角i为模块输入参数；
2.发射时刻卫星在发射轨道上真近点角fLradian和入轨时刻卫星在发射轨道上真近点角fIradian的计算；
3.发射轨道偏心率eL、发射轨道半通径p和发射轨道半长轴aL的计算
4.发射点偏近点角ELangleRadian和发射点平近点角MLangleRadian的计算
5.根据上述参数计算卫星入轨时刻tI的计算
6.发射时刻GMST格林尼治平恒星时GL的计算
7.发射时刻发射点在ECI地心惯性坐标系中赤经为：
αL＝OnLimitPiRadian(GL+λL)，其中角度约束函数OnLimitPiRadian(x)将角度x约束到[-π,π)范围内；
S10422、若为下降入轨模式，发射时刻卫星纬度辐角计算如下:
Figure FDA0002866868870000051
轨道升交点赤经:Ω＝OnLimit2PiRadian(αL-tan-1(tan(uL)×cos(i))+π)；其中角度约束函数OnLimit2PiRadian(x)将角度x约束到[0，2π)范围内；
若为上升入轨模式，发射时刻卫星纬度辐角计算如下:
Figure FDA0002866868870000052
轨道升交点赤经:Ω＝OnLimit2PiRadian(αL-tan-1(tan(uL)×cos(i))+0)；
S10423、卫星运行轨道参数初始化部分:
1.入轨时刻卫星纬度辐角uI＝uL+βL,目标点地心经纬度
Figure FDA0002866868870000053
卫星平均角速度n,地球自转转速ωE；
2.入轨点地心纬度为
Figure FDA0002866868870000054
S10424、若为下降入轨模式，则分如下4种情况计算卫星入轨后第1次和第2次过境目标纬度圈的卫星纬度辐角值uT1和uT2；目标点地心纬度为
Figure FDA0002866868870000055
入轨点地心维度
Figure FDA0002866868870000056
第1种情况.目标点在北半球,目标点地心纬度小于等于入轨点地心纬度；
第2种情况.目标点在南半球,目标点地心纬度小于等于入轨点地心纬度；
第3种情况.目标点在北半球,目标点地心纬度大于入轨点地心纬度；
第4种情况.目标点在南半球,目标点地心纬度大于入轨点地心纬度；
若
Figure FDA0002866868870000057
即目标点地心纬度小于等于入轨点地心纬度，则对于入轨后第1次过境，
Figure FDA0002866868870000061
包括第1种情况和第2种情况，如果OnLimit2PiRadian(uT1)<π则属于第1种情况，第一中间变量为Δu＝π-uT1和第二中间变量为utemp＝π+2Δu；否则属于第2种情况：Δu＝uT1-π；utemp＝π-2Δu；则入轨后第2次过境uT2＝uT1+utemp；
若
Figure FDA0002866868870000062
即目标点地心纬度大于入轨点地心纬度，则对于入轨后第1次过境，
Figure FDA0002866868870000063
包括第3种情况和第4种情况：
如果OnLimit2PiRadian(uT1)<π则属第3种情况：以tempRadian指代OnLimit2PiRadian(uT1)，Δu＝tempRadian；utemp＝π-2Δu；否则属于第4种情况：Δu＝fabs(tempRadian-2π)；utemp＝π+2Δu；则入轨后第2次过境uT2＝uT1+utemp；
若为上升入轨模式，分如下6种情况计算卫星入轨后第1次和第2次过境目标纬度圈的卫星纬度辐角值uT1和uT2；
第5种情况.目标点在北半球,目标点纬度大于入轨点纬度,入轨点在下降段；
第6种情况.目标点在南半球,目标点纬度大于入轨点纬度,入轨点在下降段；
第7种情况.目标点在北半球,目标点纬度大于入轨点纬度,入轨点在上升段；
第8种情况.目标点在南半球,目标点纬度大于入轨点纬度,入轨点在上升段；
第9种情况.目标点在北半球,目标点纬度小于入轨点纬度；
第10种情况.目标点在南半球,目标点纬度小于入轨点纬度；
若
Figure FDA0002866868870000064
即目标点地心纬度大于等于入轨点地心纬度，且入轨点在下降段(UI≥π/2)则对于入轨后第1次过境，
Figure FDA0002866868870000065
tempRadian＝OnLimit2PiRadian(uT1)；
若tempRadian＜π属于第5种情况，Δu＝tempRadian；utemp＝π-2Δu；否则属于第6种情况，则Δu＝fabs(tempRadian-2π)；utemp＝π+2Δu；
若
Figure FDA0002866868870000071
即目标点地心纬度大于等于入轨点地心纬度，且入轨点在上升段(UI＜π/2)则对于入轨后第1次过境，
Figure FDA0002866868870000072
tempRadian＝OnLimit2PiRadian(uT1)；
若tempRadian＜π属于第7种情况，Δu＝uT1；utemp＝π-2Δu；否则属于第8种情况，第8种情况不存在；
则入轨后第2次过境uT2＝uT1+utemp；
若
Figure FDA0002866868870000073
即目标点地心纬度小于入轨点地心纬度，则对于入轨后第1次过境，
Figure FDA0002866868870000074
tempRadian＝OnLimit2PiRadian(uT1)；
若tempRadian＜π属于第9种情况，
Figure FDA0002866868870000075
utemp＝π+2Δu；否则属于第10种情况，
Figure FDA0002866868870000076
utemp＝π+2Δu；
则入轨后第2次过境uT2＝uT1+utemp；
S10425、卫星运行轨道星下点过境目标纬度圈时刻地心经度计算:
1.计算卫星运行轨道周期tΩ、升交点赤经导数
Figure FDA0002866868870000077
和近地点辐角导数
Figure FDA0002866868870000078
2.计算卫星入轨后16次过境目标纬度圈的时刻:
t1＝(uT1-uI)/n,t2＝(uT2-uI)/n,t2m+1＝t1+mtΩ,t2m+2＝t2+mtΩ,m＝1,2,...,7
3.计算入轨后第1次和第2次过境目标纬度圈时星下点地心经度LonT(1)和LonT(2)；
4.计算入轨后第3次到第16次卫星星下点轨迹过境目标纬度圈时刻的地心经度:
Figure FDA0002866868870000079
Figure FDA00028668688700000710
5.计算卫星星下点轨迹16次过境目标纬度圈时刻地心经度与目标地心经度的偏差。
6.如权利要求5所述的方法，其特征在于，所述计算入轨后第1次和第2次过境目标纬度圈时星下点地心经度LonT(1)和LonT(2)；包括如下步骤：
SA1.卫星半长轴aI和当前时刻纬度辐角uT,卫星在地心轨道坐标系下的位置矢量r＝[aIcos(uT) aIsin(uT) 0]′；
SA2.根据相对入轨时刻时间差Δt和升交点赤经导数
Figure FDA0002866868870000081
计算当前时刻卫星运行轨道的升交点赤经
Figure FDA0002866868870000082
其中ΩI为卫星入轨时刻的升交点赤经；
SA3.根据相对入轨时刻时间差Δt和近地点辐角导数
Figure FDA0002866868870000083
计算当前时刻卫星运行轨道的近地点辐角
Figure FDA0002866868870000084
其中ωI为卫星入轨时刻的近地点辐角；
SA4.构建地心轨道坐标系到地心惯性坐标系的转换矩阵,其中i为轨道倾角
Figure FDA0002866868870000085
SA5.对地心轨道坐标系下的位置矢量进行坐标转转获取地心惯性ECI坐标系下的位置矢量peci
peci＝Τm×r
SA6.查询地球运动参数EOP文件,将ECI位置矢量peci转换为ECEF位置矢量pecef；
SA7.根据ECEF位置矢量pecef计算当前时刻圆形轨道卫星的星下点轨迹地心经度Lon。
7.如权利要求6所述的方法，其特征在于，所述SA6.查询地球运动参数EOP文件,将ECI位置矢量peci转换为ECEF位置矢量pecef，具体为：
SA601、从EOP文件中读取当前时刻t的6个关键地球运动参数:
1.(TAI-UTC)差值时间dat,单位秒,其中UTC为当前时刻t的世界协调时,原子时TAI＝UTC+dat
2.极移X分量xp，单位为弧度；
3.极移Y分量yp，单位为弧度；
4.(UTC1-UTC)差值时间dut,单位秒,其中UTC1为消除了极移影响后得到的世界时；
5.赤经章动ΔΨ修正量δΔΨ，单位为弧度；
6.交角章动Δε修正量δΔε，单位为弧度；
SA602、计算当前时刻t的岁差转换矩阵P(t)：
Figure FDA0002866868870000091
其中角度变量值ζ、θ和z为历元时刻T的平赤道面和平春分点相对于J2000平赤道和平春分点的3个角度值,其计算定义(公式中的角度单位为度)如下：
ζ＝(2306.2181T+0.30188*T2+0.017998T3)/3600.0
θ＝(2004.3109T-0.42665T2-0.041833T3)/3600.0
z＝(2306.2181T+1.09468T2+0.018203T3)/3600.0
定义历元时刻T指的是地球时从J2000的TT时刻起算的儒略世纪数,T＝(JDTT-2451545.0)/36525.0；其中地球时TT＝TAI+32.184,其对应的儒略日时间为JDTT；
SA603、计算当前时刻t的章动转换矩阵N(t)：
Figure FDA0002866868870000092
其中
Figure FDA0002866868870000101
Figure FDA0002866868870000102
其中章动角分量φi和系数Ai,Bi,Ci,Di根据IAU1980章动数据表进行计算；
其中历元T时刻平黄赤交角为
Figure FDA0002866868870000103
历元T时刻真黄赤交角
Figure FDA0002866868870000104
SA604、计算当前时刻t的地球自转转换矩阵R(t)：
Figure FDA0002866868870000105
其中格林尼治真恒星时
Figure FDA0002866868870000106
其中θGMST＝GMST(JDUT1)为当前世界时UTC1时刻儒略日时间的格林尼治平恒星时；
SA605、计算当前时刻t的极移转换矩阵W(t)：
Figure FDA0002866868870000107
其中极移分量xp和yp从EOP数据中获取；
SA606、t时刻ECI坐标下的坐标向量rECI转化为ECEF坐标下的坐标向量rECEF：rECEF＝[W(t)′][R(t)′][N(t)′][P(t)′]rECI。
CN202011590479.0A 2020-12-29 2020-12-29 一种太阳同步圆轨道的轨道规划方法 Pending CN112649006A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011590479.0A CN112649006A (zh) 2020-12-29 2020-12-29 一种太阳同步圆轨道的轨道规划方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011590479.0A CN112649006A (zh) 2020-12-29 2020-12-29 一种太阳同步圆轨道的轨道规划方法

Publications (1)

Publication Number Publication Date
CN112649006A true CN112649006A (zh) 2021-04-13

Family

ID=75363680

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011590479.0A Pending CN112649006A (zh) 2020-12-29 2020-12-29 一种太阳同步圆轨道的轨道规划方法

Country Status (1)

Country Link
CN (1) CN112649006A (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112800169A (zh) * 2021-04-15 2021-05-14 航天宏图信息技术股份有限公司 同步带卫星的数据匹配方法、装置、设备及存储介质
CN113632090A (zh) * 2021-06-30 2021-11-09 中国科学院微小卫星创新研究院 全球碳盘点卫星的轨道设计系统
CN114509790A (zh) * 2022-02-17 2022-05-17 北京国电高科科技有限公司 一种基于低轨卫星星座的定位方法及定位系统

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110203422A (zh) * 2019-05-31 2019-09-06 中国人民解放军63729部队 针对面目标区域探测的快速响应卫星轨道设计方法

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110203422A (zh) * 2019-05-31 2019-09-06 中国人民解放军63729部队 针对面目标区域探测的快速响应卫星轨道设计方法

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112800169A (zh) * 2021-04-15 2021-05-14 航天宏图信息技术股份有限公司 同步带卫星的数据匹配方法、装置、设备及存储介质
CN113632090A (zh) * 2021-06-30 2021-11-09 中国科学院微小卫星创新研究院 全球碳盘点卫星的轨道设计系统
CN114509790A (zh) * 2022-02-17 2022-05-17 北京国电高科科技有限公司 一种基于低轨卫星星座的定位方法及定位系统

Similar Documents

Publication Publication Date Title
CN110203422B (zh) 针对面目标区域探测的快速响应卫星轨道设计方法
Xie et al. Spacecraft dynamics and control
Di Mauro et al. Survey on guidance navigation and control requirements for spacecraft formation-flying missions
CN112629543A (zh) 一种大椭圆轨道及小倾角圆轨道的轨道规划方法
CN112649006A (zh) 一种太阳同步圆轨道的轨道规划方法
CN112783183B (zh) 一种太阳同步圆回归轨道的轨道规划方法
CN111427002B (zh) 地面测控天线指向卫星的方位角计算方法
Gil-Fernandez et al. Autonomous vision-based navigation for proximity operations around binary asteroids
CN112713922A (zh) 一种多波束通讯卫星的可见性快速预报算法
Shtark et al. Regional positioning using a low Earth orbit satellite constellation
Zeitlhöfler Nominal and observation-based attitude realization for precise orbit determination of the Jason satellites
Fujita et al. Attitude maneuvering sequence design of high-precision ground target tracking control for multispectral Earth observations
CN106777580B (zh) 近地倾斜轨道发射窗口快速设计方法
Golikov THEONA—a numerical-analytical theory of motion of artificial satellites of celestial bodies
Wood The Evolution of Deep Space Navigation: 2004–2006
Mukhayadi Efficient and high precision momentum bias attitude control for small satellite
CN116125503A (zh) 一种高精度卫星轨道确定及预报算法
Wood The evolution of deep space navigation: 1962-1989
Rondão Modeling and simulation of the ecosat-iii attitude determination and control system
Liu et al. Microsatellite autonomous orbit propagation method based on SGP4 model and GPS data
Wie et al. Attitude and orbit control systems
Centinello III et al. Orbit determination of the Dawn spacecraft with radiometric and image Data
Keil Kalman filter implementation to determine orbit and attitude of a satellite in a molniya orbit
Duxbury A spacecraft-based navigation instrument for outer planet missions
Kardashev et al. Orbit design for the Spektr-R spacecraft of the ground-space interferometer

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