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

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

Info

Publication number
CN112783183B
CN112783183B CN202011592823.XA CN202011592823A CN112783183B CN 112783183 B CN112783183 B CN 112783183B CN 202011592823 A CN202011592823 A CN 202011592823A CN 112783183 B CN112783183 B CN 112783183B
Authority
CN
China
Prior art keywords
orbit
track
point
satellite
entering
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
CN202011592823.XA
Other languages
English (en)
Other versions
CN112783183A (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.)
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 CN202011592823.XA priority Critical patent/CN112783183B/zh
Publication of CN112783183A publication Critical patent/CN112783183A/zh
Application granted granted Critical
Publication of CN112783183B publication Critical patent/CN112783183B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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/08—Control of attitude, i.e. control of roll, pitch, or yaw
    • G05D1/0808—Control of attitude, i.e. control of roll, pitch, or yaw specially adapted for aircraft
    • G05D1/0816—Control of attitude, i.e. control of roll, pitch, or yaw specially adapted for aircraft to ensure stability
    • G05D1/0833—Control of attitude, i.e. control of roll, pitch, or yaw specially adapted for aircraft to ensure stability using limited authority control

Landscapes

  • Engineering & Computer Science (AREA)
  • Computer Security & Cryptography (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)
  • 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种入轨模式轨道倾角遍历范围内所有卫星轨道，调用太阳同步圆回归轨道卫星的发射轨道和运行轨道迭代计算OnIterativeSunCompute模块计算第j次过境目标纬度圈时星下点经度分布,选取轨道倾角遍历范围内过境时刻星下点经度与目标点经度偏差值最小的轨道。
步骤5、判断步骤4中所选取的偏差值最小的轨道对应的偏差值是否满足预设的指标要求，若是则进入步骤6，否则MO自增1，返回步骤3。
步骤6、判断当前卫星载荷类型是否为CCD光学设备，若是则进入步骤7，否则进入步骤10。
步骤7、判断当前轨道卫星过境目标点时刻是否在CCD光学设备的UTC时间观测窗口内，若是则进入步骤10，否则进入步骤8。
步骤8、计算发射延时，对于当前入轨模式轨道倾角遍历范围内的所有卫星轨道，调用太阳同步圆回归轨道卫星的发射轨道和运行轨道迭代计算OnIterativeSunCompute模块计算当前第j次过境目标纬度圈时星下点经度分布,选取遍历范围内过境时刻星下点经度与目标点经度偏差值最小的轨道。
步骤9、判断步骤8中所选取的偏差值最小的轨道是否满足预设的指标要求，若是则进入步骤10，否则返回步骤3。
步骤10、当前第BackupOrbitNum个轨道规划方案满足初步要求，将当前第BackupOrbitNum个轨道规划方案存入在BackupOrbit数组中,并且BackupOrbitNum自增1。
步骤11、太阳同步圆回归轨道规划完成，将当前对应的共BackupOrbitNum个轨道方案存入预先构建的BackupOrbit数组中。
进一步地，步骤4中，太阳同步圆回归轨道卫星的发射轨道和运行轨道迭代计算OnIterativeSunCompute模块，具体包括如下步骤：
步骤201、根据界面输入的轨道高度上下限和轨道回归系数，调用太阳同步圆回归轨道半长轴迭代计算函数OnComputeAa模块计算判断当前输入的轨道高度上下限和轨道回归系数是否有效；如果无效,则返回值Flag＝false；如果有效，则返回值Flag＝true，根据轨道回归系数迭代计算轨道半长轴和轨道倾角，结果保存在SunOrbit数组中；
步骤202、判断返回值Flag是否有效，若是则进入步骤203，否则输出false；
步骤203、根据SunOrbit数组中保存的轨道半长轴和轨道倾角，以及指定的发射区域中心位置参数，调用发射点位置遍历选取进行太阳同步圆回归轨道设计OnComputeOrbitMinValueSunRepeatOrbit模块在发射区域中心位置附近遍历选取发射点位置，计算对应不同发射点位置太阳同步圆回归轨道星下点轨迹第j次过境目标纬度圈时的星下点地心经度与目标点地心经度的偏差分布，选取偏差值最小的轨道设计进行门限判决,并返回有效轨道设计的轨道倾角IiAngle。
进一步地，步骤203中，发射点位置遍历选取进行太阳同步圆回归轨道设计OnComputeOrbitMinValueSunRepeatOrbit模块，具体包括如下步骤：
步骤301、太阳同步圆回归轨道设计循环计算初始化:，从SunOrbit数组中读取当前的轨道半长轴aa和轨道倾角IiAngleDegree；设定当前轨道索引SunOrbitInitialNum初值为0；在发射区域中心位置附近均匀遍历选取发射点位置；
步骤302、调用太阳同步圆轨道卫星的发射轨道和运行轨道计算OnComputeOrbitValueSunOrbit模块计算对应不同发射点位置太阳同步圆回归轨道卫星入轨后星下点16次过境目标纬度圈时星下点地心经度与目标点地心经度的偏差分布；当前轨道索引SunOrbitInitialNum自增1；
步骤303、判断判断当前发射点位置遍历是否完成，若是则进入步骤304，否则选取下一发射点位置，返回步骤302；
步骤304、对于遍历范围内SunOrbitInitialNum条卫星轨道星下点轨迹第j次过境目标纬度圈时的星下点地心经度与目标点地心经度的偏差分布进行统计，选取偏差值最小的轨道方案,并对该最小偏差值minValue进行门限threshold判决；
步骤305、对最小偏差值minValue进行门限threshold判决，若minValue<＝threshold，则进入步骤306；否则返回无效值；
步骤306、minValue对应的轨道方案保存在minTimeOrbit数组中，并返回有效轨道倾角值。
4、如权利要求1～3任意的方法，其特征在于，太阳同步圆轨道卫星的发射轨道和运行轨道计算OnComputeOrbitValueSunOrbit模块，具体包括如下步骤：
S401、卫星发射轨道参数初始化计算部分：
①.发射点地心经纬度为发射轨道地心角βL和发射轨道倾角i为模块输入参数；
②.发射时刻卫星在发射轨道上真近点角fLradian和入轨时刻卫星在发射轨道上真近点角fIradian的计算；
③.发射轨道偏心率eL、发射轨道半通径p和发射轨道半长轴aL的计算
④.发射点偏近点角ELangleRadian和发射点平近点角MLangleRadian的计算
⑤.根据上述参数计算卫星入轨时刻tI的计算
⑥.发射时刻GMST格林尼治平恒星时GL的计算
⑦.发射时刻发射点在ECI地心惯性坐标系中赤经为：
αL＝OnLimitPiRadian(GL+λL)，其中角度约束函数OnLimitPiRadian(x)将角度x约束到[-π,π)范围内。
S402、若为下降入轨模式，发射时刻卫星纬度辐角计算如下:
轨道升交点赤经:Ω＝OnLimit2PiRadian(αL-tan-1(tan(uL)×cos(i))+π)；其中角度约束函数OnLimit2PiRadian(x)将角度x约束到[0，2π)范围内。
若为上升入轨模式，发射时刻卫星纬度辐角计算如下:
轨道升交点赤经:Ω＝OnLimit2PiRadian(αL-tan-1(tan(uL)×cos(i))+0)；
S403、卫星运行轨道参数初始化部分:
1.入轨时刻卫星纬度辐角uI＝uL+βL,目标点地心经纬度卫星平均角速度n,地球自转转速ωE；
2.入轨点地心纬度为
S404、若为下降入轨模式，则分如下4种情况计算卫星入轨后第1次和第2次过境目标纬度圈的卫星纬度辐角值uT1和uT2；目标点地心纬度为入轨点地心维度/>
第1种情况.目标点在北半球,目标点地心纬度小于等于入轨点地心纬度；
第2种情况.目标点在南半球,目标点地心纬度小于等于入轨点地心纬度；
第3种情况.目标点在北半球,目标点地心纬度大于入轨点地心纬度；
第4种情况.目标点在南半球,目标点地心纬度大于入轨点地心纬度；
若即目标点地心纬度小于等于入轨点地心纬度，则对于入轨后第1次过境，/>包括第1种情况和第2种情况，如果OnLimit2PiRadian(uT1)<π则属于第1种情况：，第一中间变量为Δu＝π-uT1；第一中间变量为utemp＝π+2Δu；否则属于第2种情况：Δu＝uT1-π；utemp＝π-2Δu；则入轨后第2次过境uT2＝uT1+utemp；
若即目标点地心纬度大于入轨点地心纬度，则对于入轨后第1次过境，包括第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种情况.目标点在南半球,目标点纬度小于入轨点纬度；
若即目标点地心纬度大于等于入轨点地心纬度，且入轨点在下降段(UI≥π/2)则对于入轨后第1次过境，/>tempRadian＝OnLimit2PiRadian(uT1)；
若tempRadian＜π属于第5种情况，Δu＝tempRadian；utemp＝π-2Δu；否则属于第6种情况，则Δu＝fabs(tempRadian-2π)；utemp＝π+2Δu；
若即目标点地心纬度大于等于入轨点地心纬度，且入轨点在上升段(UI＜π/2)则对于入轨后第1次过境，/>tempRadian＝OnLimit2PiRadian(uT1)；
若tempRadian＜π属于第7种情况，Δu＝uT1；utemp＝π-2Δu；否则属于第8种情况，第8种情况不存在；
则入轨后第2次过境uT2＝uT1+utemp。
若即目标点地心纬度小于入轨点地心纬度，则对于入轨后第1次过境，tempRadian＝OnLimit2PiRadian(uT1)；
若tempRadian＜π属于第9种情况，utemp＝π+2Δu；否则属于第10种情况，/>utemp＝π+2Δu；
则入轨后第2次过境uT2＝uT1+utemp；
S405、卫星运行轨道星下点过境目标纬度圈时刻地心经度计算：
1).计算卫星运行轨道周期tΩ、升交点赤经导数和近地点辐角导数/>
2).计算卫星入轨后16次过境目标纬度圈的时刻:
t1＝(uT1-uI)/n,t2＝(uT2-uI)/n,t2m+1＝t1+mtΩ,t2m+2＝t2+mtΩ,m＝1,2,...,7，m为次数；
3).调用圆形轨道卫星的星下点轨迹计算OnComputeValueLonGCcircle模块计算入轨后第1次和第2次过境目标纬度圈时星下点地心经度LonT(1)和LonT(2)。
4).计算入轨后第3次到第16次卫星星下点轨迹过境目标纬度圈时刻的地心经度:
5).计算卫星星下点轨迹16次过境目标纬度圈时刻地心经度与目标地心经度的偏差。
进一步地，调用圆形轨道卫星的星下点轨迹计算OnComputeValueLonGCcircle模块，具体包括如下步骤：
SA1.卫星半长轴aI和当前时刻纬度辐角uT,卫星在地心轨道坐标系下的位置矢量r＝[aI cos(uT) aI sin(uT) 0]′；
SA2.根据相对入轨时刻时间差Δt和升交点赤经导数计算当前时刻卫星运行轨道的升交点赤经/>其中ΩI为卫星入轨时刻的升交点赤经；
SA3.根据相对入轨时刻时间差Δt和近地点辐角导数计算当前时刻卫星运行轨道的近地点辐角/>其中ωI为卫星入轨时刻的近地点辐角；
SA4.构建地心轨道坐标系到地心惯性坐标系的转换矩阵,其中i为轨道倾角
SA 5.对地心轨道坐标系下的位置矢量进行坐标转转获取地心惯性ECI坐标系下的位置矢量peci
peci＝Tm×r
SA 6.查询地球运动参数EOP文件,调用ECI坐标系转ECEF坐标系OnConvertECItoECEF模块将ECI位置矢量peci转换为ECEF位置矢量pecef；
SA7.根据ECEF位置矢量pecef计算当前时刻圆形轨道卫星的星下点轨迹地心经度Lon。
进一步地，SA 6.ECI坐标系转ECEF坐标系OnConvertECItoECEF模块，具体为：
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)：
其中角度变量值ζ、θ和z为历元时刻T的平赤道面和平春分点相对于J2000平赤道和平春分点的3个角度值,其计算定义(公式中的角度单位为度)如下：
ζ＝(2306.2181T+0.30188T2+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)：
其中
其中章动角分量φi和系数Ai,Bi,Ci,Di根据IAU1980章动数据表进行计算；
其中历元T时刻平黄赤交角为
历元T时刻真黄赤交角
SA604、计算当前时刻t的地球自转转换矩阵R(t)：
其中格林尼治真恒星时其中θGMST＝GMST(JDUT1)为当前世界时UTC1时刻儒略日时间的格林尼治平恒星时；
SA605、计算当前时刻t的极移转换矩阵W(t)：
其中极移分量xp和yp从EOP数据中获取；
SA606、t时刻ECI坐标下的坐标向量rECI转化为ECEF坐标下的坐标向量rECEF：rECEF＝[W(t)′][R(t)′][N(t)′][P(t)′]rECI。
进一步地，步骤201中，太阳同步圆回归轨道半长轴迭代计算函数OnComputeAa模块，具体包括如下步骤：
S501、根据界面输入的轨道高度下限和上限，调用太阳同步圆回归轨道回归系数计算OnComputeSunOrbitBackCoeff模块计算太阳同步圆回归轨道回归系数的上限Qmax和下限Qmin,当前根据界面输入得到的回归系数为m_SatBackCoeff；
S502、判断若Qmin>10000.0)||(Qmax>10000.0成立，则输出无效，否则进入S503；
S503、判断若m_SatBackCoeff<Qmin成立，则输出无效，否则进入S504；
S504、判断若m_SatBackCoeff>Qmax成立，则输出无效，否则进入S505；
S505、根据回归系数m_SatBackCoeff调用太阳同步圆回归轨道半长轴迭代计算OnSunOrbitAa模块迭代计算太阳同步圆回归轨道的半长轴Aa；
S506、判断Aa<0.0是否成立，若是输出无效，否则进入S507；
S507、根据计算获得的太阳同步圆回归轨道半长轴Aa调用根据轨道半长轴计算太阳同步圆轨道倾角函数OnComputeSunOrbitAngle模块迭代计算太阳同步圆回归轨道的轨道倾角IiAngle,并将Aa和IiAngle参数存放在SunOrbit数组中，输出有效。
进一步地，S505中，太阳同步圆回归轨道半长轴迭代计算OnSunOrbitAa模块，具体为：
S801、根据从界面输入的轨道入轨点高度上限和下限分别计算太阳同步圆轨道半长轴上限aU和下限aL；
根据轨道半长轴下限和上限，调用太阳同步圆回归轨道回归系数计算OnComputeSunOrbitBackCoeff模块分别计算太阳同步圆回归轨道回归系数的上限Qmax＝OnComputeSunOrbitBackCoeff(aL)和下限Qmin＝OnComputeSunOrbitBackCoeff(aU)；
当前根据界面输入得到的回归系数为m_SatBackCoeff；
S802、判断若(Qmax>100000.0)||Qmin>100000.0)成立，则返回无效值，否则进入S803；
S803、计算第一差值tempValue1＝Qmax-SatBackCoeff；以及第二差值tempValue2＝Qmin-SatBackCoeff；
S804、判断(tempValue1×tempValue2)>0.0是否成立，若是则返回无效值，否则进入S805；
S805、轨道半长轴迭代计算初始化:迭代次数index初始化为0,差值tempValue1和tempValue2；
S806、判断index<30是否成立，若是则进入S807，否则输出有效值aM；
S807、轨道半长轴迭代计算主体部分包括如下步骤：
计算轨道半长轴中间值aM＝(aU+aL)/2.0,index自增1；
根据轨道半长轴中间值aM，调用太阳同步圆回归轨道回归系数计算OnComputeSunOrbitBackCoeff模块计算当前轨道升交点赤经导数tempValue；
S808、判断tempValue>10000.0是否成立，若是则输出无效值；否则进入S809；
S809、判断tempValue×tempValue1>0.0是否成立，若是则令tempValue1＝tempValue,aL＝aM；否则令tempValue2＝tempValue,aU＝aM；
返回S806。
进一步地，太阳同步圆回归轨道回归系数计算OnComputeSunOrbitBackCoeff模块具体为：
S901、根据输入参数轨道半长轴aa,调用根据轨道半长轴计算太阳同步圆轨道倾角函数OnComputeSunOrbitAngle模块计算太阳同步圆回归轨道的轨道倾角IiAngle；
S902、判断IiAngle>10000.0是否成立，若是则返回无效值，否则进入S903；
S903、根据输入参数轨道半长轴aa和前面获得的有效轨道倾角IiAngle，调用圆轨道回归系数计算OnComputeBackCoeffCircle模块计算太阳同步圆回归轨道回归系数Q，返回太阳同步圆回归轨道回归系数Q。
进一步地，圆轨道回归系数计算OnComputeBackCoeffCircle模块，具体为：
根据轨道半长轴aI、轨道倾角i、地球平均半径Re、地球自转速度ωE和卫星平均角速度n，计算圆轨道回归系数:
返回圆轨道回归系数fa；
根据轨道半长轴计算太阳同步圆轨道倾角函数OnComputeSunOrbitAngle模块，具体包括如下步骤：
S1101.设置太阳同步圆轨道倾角上限angleU＝120度和下限angleL＝96度
根据当前模块的输入参数轨道半长轴aa和轨道倾角上下限，调用太阳同步圆轨道升交点赤经导数偏差计算OnComputeOmegaDotSunOrbit模块分别计算当前轨道升交点赤经导数与太阳同步轨道标准升交点赤经导数的差值tempValueU和tempValueL
S1102.判断(tempValueL×tempValueU)>0.0是否成立，若是则返回无效值，否则进入S1103；
S1103.轨道半长轴迭代计算初始化:迭代次数index＝0,差值tempValueL和tempValueU；
S1104.判断index<30是否成立，若是则进入S1105，否则返回angleM；
S1105.轨道倾角迭代计算主体部分为：
计算轨道倾角中间值angleM＝(angleU+angleL)/2.0
根据当前输入参数轨道半长轴aa和轨道倾角中间值中间值angleM，调用太阳同步圆轨道升交点赤经导数偏差计算OnComputeOmegaDotSunOrbit模块分别计算当前轨道升交点赤经导数与太阳同步轨道标准升交点赤经导数的差值tempValue；
index自增1；
S1106.判断(tempValue×tempValueL)>0.0是否成立，若是则令tempValueL＝tempValue,angleL＝angleM，否则，令tempValueU＝tempValue,angleU＝angleM；
返回S1104；
太阳同步圆轨道升交点赤经导数偏差计算OnComputeOmegaDotSunOrbit模块，具体为：
根据轨道半长轴aI、轨道倾角i、地球平均半径Re、地球自转速度ωE和卫星平均角速度n，计算当前轨道升交点赤经导数与太阳同步圆轨道标准升交点赤经导数的偏差值：
有益效果：
本发明提供了一种太阳同步圆回归轨道的轨道规划方法，太阳同步圆回归轨道设计要满足太阳同步轨道升交点赤经导数约束和回归系数约束这两个指标要求。根据这两个指标要求，可以通过迭代计算获得太阳同步圆回归轨道的轨道半长轴和轨道倾角这两个重要参数。根据计算获得的轨道半长轴和轨道倾角，以及指定发射区域中心位置等参数，在发射区域中心位置附近遍历选取发射点位置，计算对应不同发射点位置太阳同步圆回归轨道星下点多次过境目标纬度圈时的星下点地心经度与目标点地心经度的偏差分布，选取偏差值最小的轨道设计进行门限判决，并保存有效轨道设计。此外在卫星轨道设计流程中考虑了太阳光照参数对于光学载荷卫星轨道设计的实施要求,并根据轨道参数计算在太阳高度角约束条件下的卫星发射窗口。因此本发明能够结合卫星发射轨道模型、轨道回归特性约束、目标区域光照条件约束和精确地球运动模型进行轨道规划，实现了对目标区域的快速有效探测。
附图说明
图1为卫星发射和入轨过程示意图；
图2态势1：目标点地心纬度小于入轨点地心纬度，目标点在北半球；
图3态势2：目标点地心纬度小于入轨点地心纬度，目标点在南半球；
图4为态势3：目标点地心纬度大于入轨点地心纬度，目标点在北半球；
图5为态势4：目标点地心纬度大于入轨点地心纬度，目标点在南半球；
图6为态势5：目标点地心纬度大于入轨点地心纬度，入轨点在下降段，目标点在北半球；
图7为态势6：目标点地心纬度大于入轨点地心纬度，入轨点在下降段，目标点在南半球；
图8为态势7：目标点地心纬度大于入轨点地心纬度，入轨点在上升段，目标点在北半球；
图9为态势8(无效)：目标点地心纬度大于入轨点地心纬度，入轨点在上升段，目标点在南半球；
图10为态势9示意图：目标点地心纬度小于入轨点地心纬度，目标点在北半球；
图11为态势10：目标点地心纬度小于入轨点地心纬度，目标点在南半球；
图12为本发明实施例中的太阳同步圆回归轨道设计OnComputeSunRepeatOrbit模块的流程图；
图13为本发明实施例中的太阳同步圆回归轨道卫星的发射轨道和运行轨道迭代计算OnIterativeSunCompute模块的流程图；
图14为本发明实施例中的发射点位置遍历选取进行太阳同步圆回归轨道设计OnComputeOrbitMinValueSunRepeatOrbit模块的流程图；
图15为本发明实施例中的太阳同步圆轨道卫星的发射轨道和运行轨道计算OnComputeOrbitValueSunOrbit模块的流程图；
图16为本发明实施例中的圆形轨道卫星的星下点轨迹计算OnComputeValueLonGCcircle模块流程图；
图17为本发明实施例中的ECI坐标系转ECEF坐标系OnConvertECItoECEF模块流程图；
图18为本发明实施例中的太阳同步圆回归轨道半长轴迭代计算函数OnComputeAa模块的流程图；
图19为本发明实施例中的太阳同步圆回归轨道半长轴迭代计算OnSunOrbitAa模块的流程图；
图20为本发明实施例中的太阳同步圆回归轨道回归系数计算OnComputeSunOrbitBackCoeff模块的流程图
图21为本发明实施例中的圆轨道回归系数计算OnComputeBackCoeffCircle模块流程图；
图22为本发明实施例中的根据轨道半长轴计算太阳同步圆轨道倾角函数OnComputeSunOrbitAngle模块流程图；
图23为本发明实施例中的太阳同步圆轨道升交点赤经导数偏差计算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的地心经纬度为入轨点I的地心经纬度为/>设置地面目标点的地心经纬度为/>定义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的计算如下：
发射点的偏近点角EL的计算如下：
发射点的平近点角ML的计算如下：
ML＝EL-eL sin(EL)
发射轨道半长轴aL计算如下：
卫星入轨时刻tI的计算如下：
其中μ＝398600.4415km3/s2为引力常数。
对于地心经纬度为的发射点L，根据球面三角公式，则发射点地心纬度/>轨道倾角i和发射点卫星纬度辐角uL具有如下重要关系：
当前本轨道规划设计中，应急发射点只考虑在国内的情况(即在北半球)，则所以可得0＜uL＜π。
根据卫星发射入轨的四种模式(1逆行下降，2逆行上升，3顺行下降，4顺行上升),可以分别计算发射点卫星纬度辐角uL。
对于卫星下降入轨模式1和模式3的发射点卫星纬度辐角uL计算如下：
对于卫星上升入轨模式2和模式4的发射点卫星纬度辐角uL计算如下：
定义发射时刻tL的格林尼治平恒星时GMST(瞬时平春分点和格林尼治子午圈之间的夹角)为GL，定义发射时刻tL在ECI地心惯性坐标系中地心经纬度为的发射点L赤经αL计算如下：
αL＝OnLimitPiRadian(GL+λL)
其中角度约束函数OnLimitPiRadian(x)将角度x约束到[-π，π)范围内。对于地心经纬度为的发射点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。对于地心经纬度为的入轨点I，根据球面三角公式，则入轨点地心纬度/>轨道倾角i和入轨点卫星纬度辐角uI具有如下重要关系：/>
如果入轨点在北半球，则卫星纬度辐角uI＜π
如果入轨点在南半球，则卫星纬度辐角uI＞π
不管入轨点是在北半球还是在南半球，入轨点的地心纬度计算如下所示：
对于卫星下降入轨模式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)转换矩阵。
定义旋转矩阵算子如下:
定义旋转矩阵算子如下:
定义原子时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)的计算如下:
其中角度变量值ζ,θ和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)的计算如下:
赤经章动ΔΨ计算如下:
交角章动Δε计算如下:
其中δΔΨ为赤经章动ΔΨ修正量,δΔε为交角章动Δε修正量，可以从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
定义为时刻t真黄(ecliptic)赤(equator)交角,其中/>为时刻t平黄赤交角,Δε为交角章动值。
平黄赤交角(计算的角度单位为度)的计算如下:
(3).地球自转转换矩阵R(t)计算
自转转换矩阵R(t)是由格林尼治真恒星时θGAST的影响造成，其计算如下:
其中格林尼治真恒星时θGAST(单位为弧度)的计算如下:
其中θGMST＝GMST(JDUT1)为当前世界时UTC1时刻的格林尼治平恒星时,JDUT1为世界时UTC1时刻对应的儒略日时间，函数GMST(JD)用于计算儒略日时间JD的格林尼治平恒星时。
其中赤经章动ΔΨ和交角章动Δε的计算参见前面公式。
(4).极移转换矩阵W(t)计算
由于地球不是刚体以及其他一些物理因素的影响，造成地球自转轴在地球体内的运动，引起地球极点在地球表面的位置随时间而变化，称为极移。
极移转换矩阵W(t)计算如下:
其中xp为极移X分量,yp为极移Y分量，可以从EOP文件中获取当前时刻对应的xp和yp值。
太阳同步圆回归轨道设计要满足太阳同步轨道升交点赤经导数约束和回归系数约束这两个指标要求。根据这两个指标要求，可以通过迭代计算获得太阳同步圆回归轨道的轨道半长轴和轨道倾角这两个重要参数。根据计算获得的轨道半长轴和轨道倾角，以及指定发射区域中心位置等参数，在发射区域中心位置附近遍历选取发射点位置，计算对应不同发射点位置太阳同步圆回归轨道星下点多次过境目标纬度圈时的星下点地心经度与目标点地心经度的偏差分布，选取偏差值最小的轨道设计进行门限判决，并保存有效轨道设计。此外在卫星轨道设计流程中考虑了太阳光照参数对于光学载荷卫星轨道设计的实施要求,并根据轨道参数计算在太阳高度角约束条件下的卫星发射窗口。
圆形轨道回归系数fa的计算如下所示：
太阳同步圆轨道的升交点赤经导数约束如下：
其中轨道半长轴aI、轨道倾角i、地球平均半径Re、地球自转速度ωE和卫星平均角速度n。太阳同步圆回归轨道设计中，本发明对于2种发射入轨模式遍历范围内的每条卫星轨道，计算卫星入轨后16次卫星星下点过境目标纬度圈时刻的ECEF地心地固坐标系下大地经度值。当前卫星运行轨道为圆形轨道，选取卫星轨道升交点与轨道近地点重合，则运行轨道的近地点辐角ωI＝0。根据J2摄动条件下的卫星圆形运行轨道方程，可以计算出t1、t2到t16时刻16次卫星星下点过境目标纬度圈时刻在ECEF地心地固坐标系下的大地经度LonT(1)、LonT(2)到LonT(16)。
太阳同步圆回归轨道设计OnComputeSunRepeatOrbit模块对于当前2种入轨模式下的轨道倾角遍历范围内所有卫星轨道，调用太阳同步圆回归轨道卫星的发射轨道和运行轨道迭代计算OnIterativeSunCompute模块计算卫星入轨后星下点16次过境(根据经验设定过境次数)目标纬度圈时的经度分布，选取遍历范围内过境时刻星下点经度与目标点经度偏差值最小的轨道进行门限判决，并保存有效轨道的计算参数到BackupOrbit数组。在卫星轨道设计流程中还需要考虑太阳光照参数对于光学载荷卫星轨道设计的实施要求,并根据轨道参数计算在太阳高度角约束条件下的卫星发射窗口。
如图12所示为本发明实施例提供的太阳同步圆回归轨道设计OnComputeSunRepeatOrbit模块流程图，即一种太阳同步圆回归轨道的轨道规划方法，包括如下步骤：
步骤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种入轨模式轨道倾角遍历范围内所有卫星轨道，调用太阳同步圆回归轨道卫星的发射轨道和运行轨道迭代计算OnIterativeSunCompute模块计算第j次过境目标纬度圈时星下点经度分布,选取轨道倾角遍历范围内过境时刻星下点经度与目标点经度偏差值最小的轨道。
步骤5、判断步骤4中所选取的偏差值最小的轨道对应的偏差值是否满足预设的指标要求，若是则进入步骤6，否则MO自增1，返回步骤3。
步骤6、判断当前卫星载荷类型是否为CCD光学设备，若是则进入步骤7，否则进入步骤10。
步骤7、判断当前轨道卫星过境目标点时刻是否在CCD光学设备的UTC时间观测窗口内，若是则进入步骤10，否则进入步骤8。
步骤8、计算发射延时，对于当前入轨模式轨道倾角遍历范围内的所有卫星轨道，调用太阳同步圆回归轨道卫星的发射轨道和运行轨道迭代计算OnIterativeSunCompute模块计算当前第j次过境目标纬度圈时星下点经度分布,选取遍历范围内过境时刻星下点经度与目标点经度偏差值最小的轨道。
步骤9、判断步骤8中所选取的偏差值最小的轨道是否满足预设的指标要求，若是则进入步骤10，否则返回步骤3。
步骤10、当前第BackupOrbitNum个轨道规划方案满足初步要求，将当前第BackupOrbitNum个轨道规划方案存入在BackupOrbit数组中,并且BackupOrbitNum自增1。
步骤11、太阳同步圆轨道规划完成，将当前对应的共BackupOrbitNum个轨道方案存入预先构建的BackupOrbit数组中。
图13所示为太阳同步圆回归轨道卫星的发射轨道和运行轨道迭代计算OnIterativeSunCompute模块的流程图，具体包括如下步骤：
步骤201、根据界面输入的轨道高度上下限和轨道回归系数，调用太阳同步圆回归轨道半长轴迭代计算函数OnComputeAa模块计算判断当前输入的轨道高度上下限和轨道回归系数是否有效；如果无效,则返回值Flag＝false；如果有效，则返回值Flag＝true，根据轨道回归系数迭代计算轨道半长轴和轨道倾角，结果保存在SunOrbit数组中。
步骤202、判断返回值Flag是否有效，若是则进入步骤203，否则输出false。
步骤203、根据SunOrbit数组中保存的轨道半长轴和轨道倾角，以及指定的发射区域中心位置参数，调用发射点位置遍历选取进行太阳同步圆回归轨道设计OnComputeOrbitMinValueSunRepeatOrbit模块在发射区域中心位置附近遍历选取发射点位置，计算对应不同发射点位置太阳同步圆回归轨道星下点轨迹第j次(j从1到16,为模块输入参数)过境目标纬度圈时的星下点地心经度与目标点地心经度的偏差分布，选取偏差值最小的轨道设计进行门限判决,并返回有效轨道设计的轨道倾角IiAngle。
步骤203中，其中在发射区域中心位置附近均匀遍历运载火箭发射点用于进行发射点位置遍历选取进行太阳同步圆回归轨道设计OnComputeOrbitMinValueSunRepeatOrbit模块流程图如图14所示，具体包括如下步骤：
步骤301、太阳同步圆回归轨道设计循环计算初始化:，从SunOrbit数组中读取当前的轨道半长轴aa和轨道倾角IiAngleDegree；设定当前轨道索引SunOrbitInitialNum初值为0；在发射区域中心位置附近均匀遍历选取发射点位置。
步骤302、调用太阳同步圆轨道卫星的发射轨道和运行轨道计算OnComputeOrbitValueSunOrbit模块计算对应不同发射点位置太阳同步圆回归轨道卫星入轨后星下点16次过境目标纬度圈时星下点地心经度与目标点地心经度的偏差分布；当前轨道索引SunOrbitInitialNum自增1。
步骤303、判断判断当前发射点位置遍历是否完成，若是则进入步骤304，否则选取下一发射点位置，返回步骤302。
步骤304、对于遍历范围内SunOrbitInitialNum条卫星轨道星下点轨迹第j次(j从1到16,为模块输入参数)过境目标纬度圈时的星下点地心经度与目标点地心经度的偏差分布进行统计，选取偏差值最小的轨道方案,并对该最小偏差值minValue进行门限threshold判决。
步骤305、对最小偏差值minValue进行门限threshold判决，若minValue<＝threshold，则进入步骤306；否则返回无效值。
步骤306、minValue对应的轨道方案保存在minTimeOrbit数组中，并返回有效轨道倾角值。
太阳同步圆轨道卫星的发射轨道和运行轨道计算OnComputeOrbitValueSunOrbit核心计算模块流程图如图15所示，具体包括如下步骤：
S401、卫星发射轨道参数初始化计算部分：
①.发射点地心经纬度为发射轨道地心角βL和发射轨道倾角i为模块输入参数。
②.发射时刻卫星在发射轨道上真近点角fLradian和入轨时刻卫星在发射轨道上真近点角fIradian的计算。
③.发射轨道偏心率eL、发射轨道半通径p和发射轨道半长轴aL的计算。
④.发射点偏近点角ELangleRadian和发射点平近点角MLangleRadian的计算。
⑤.根据上述参数计算卫星入轨时刻tI的计算。
⑥.发射时刻GMST格林尼治平恒星时GL的计算。
⑦.发射时刻发射点在ECI地心惯性坐标系中赤经为：
αL＝OnLimitPiRadian(GL+λL)。
S402、若为下降入轨模式，发射时刻卫星纬度辐角计算如下:
轨道升交点赤经:Ω＝OnLimit2PiRadian(αL-tan-1(tan(uL)×cos(i))+π)。
若为上升入轨模式，发射时刻卫星纬度辐角计算如下:
轨道升交点赤经:Ω＝OnLimit2PiRadian(αL-tan-1(tan(uL)×cos(i))+0)。
S403、卫星运行轨道参数初始化部分:
1.入轨时刻卫星纬度辐角uI＝uL+βL,目标点地心经纬度卫星平均角速度n,地球自转转速ωE。
2.入轨点地心纬度为
S404、若为下降入轨模式，则分如下4种情况计算卫星入轨后第1次和第2次过境目标纬度圈的卫星纬度辐角值uT1和uT2；目标点地心纬度为入轨点地心维度/>
第1种情况.目标点在北半球,目标点地心纬度小于等于入轨点地心纬度；
第2种情况.目标点在南半球,目标点地心纬度小于等于入轨点地心纬度；
第3种情况.目标点在北半球,目标点地心纬度大于入轨点地心纬度；
第4种情况.目标点在南半球,目标点地心纬度大于入轨点地心纬度；
若即目标点地心纬度小于等于入轨点地心纬度，则对于入轨后第1次过境，/>包括第1种情况和第2种情况，如果OnLimit2PiRadian(uT1)<π则属于第1种情况：，第一中间变量为Δu＝π-uT1；第一中间变量为utemp＝π+2Δu；否则属于第2种情况：Δu＝uT1-π；utemp＝π-2Δu；则入轨后第2次过境uT2＝uT1+utemp。
若即目标点地心纬度大于入轨点地心纬度，则对于入轨后第1次过境，/>包括第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种情况.目标点在南半球,目标点纬度小于入轨点纬度；
若即目标点地心纬度大于等于入轨点地心纬度，且入轨点在下降段(UI≥π/2)则对于入轨后第1次过境，/>tempRadian＝OnLimit2PiRadian(uT1)。
若tempRadian＜π属于第5种情况，Δu＝tempRadian；utemp＝π-2Δu；否则属于第6种情况，则Δu＝fabs(tempRadian-2π)；utemp＝π+2Δu。
若即目标点地心纬度大于等于入轨点地心纬度，且入轨点在上升段(UI＜π/2)则对于入轨后第1次过境，/>tempRadian＝OnLimit2PiRadian(uT1)。
若tempRadian＜π属于第7种情况，Δu＝uT1；utemp＝π-2Δu；否则属于第8种情况，第8种情况不存在。
则入轨后第2次过境uT2＝uT1+utemp。
若即目标点地心纬度小于入轨点地心纬度，则对于入轨后第1次过境，tempRadian＝OnLimit2PiRadian(uT1)。
若tempRadian＜π属于第9种情况，utemp＝π+2Δu；否则属于第10种情况，/>utemp＝π+2Δu。
则入轨后第2次过境uT2＝uT1+utemp。
S405、卫星运行轨道星下点过境目标纬度圈时刻地心经度计算：
1).计算卫星运行轨道周期tΩ、升交点赤经导数和近地点辐角导数/>
2).计算卫星入轨后16次过境目标纬度圈的时刻:
t1＝(uT1-uI)/n,t2＝(uT2-uI)/n,t2m+1＝t1+mtΩ,t2m+2＝t2+mtΩ,m＝1,2,...,7
m为次数。
3).调用圆形轨道卫星的星下点轨迹计算OnComputeValueLonGCcircle模块计算入轨后第1次和第2次过境目标纬度圈时星下点地心经度LonT(1)和LonT(2)。
4).计算入轨后第3次到第16次卫星星下点轨迹过境目标纬度圈时刻的地心经度:
5).计算卫星星下点轨迹16次过境目标纬度圈时刻地心经度与目标地心经度的偏差。
圆形轨道卫星的星下点轨迹计算OnComputeValueLonGCcircle模块流程如图16所示，调用的圆形轨道卫星的星下点轨迹计算OnComputeValueLonGCcircle模块，具体包括如下步骤：
SA1.卫星半长轴aI和当前时刻纬度辐角uT,卫星在地心轨道坐标系下的位置矢量r＝[aI cos(uT) aI sin(uT) 0]′；
SA2.根据相对入轨时刻时间差Δt和升交点赤经导数计算当前时刻卫星运行轨道的升交点赤经/>其中ΩI为卫星入轨时刻的升交点赤经；
SA3.根据相对入轨时刻时间差Δt和近地点辐角导数计算当前时刻卫星运行轨道的近地点辐角/>其中ωI为卫星入轨时刻的近地点辐角；
SA4.构建地心轨道坐标系到地心惯性坐标系的转换矩阵,其中i为轨道倾角
SA 5.对地心轨道坐标系下的位置矢量进行坐标转转获取地心惯性ECI坐标系下的位置矢量peci
peci＝Tm×r
SA6.查询地球运动参数EOP文件,调用ECI坐标系转ECEF坐标系OnConvertECItoECEF模块将ECI位置矢量peci转换为ECEF位置矢量pecef；
SA7.根据ECEF位置矢量pecef计算当前时刻圆形轨道卫星的星下点轨迹地心经度Lon。
SA 6中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)：
其中角度变量值ζ、θ和z为历元时刻T的平赤道面和平春分点相对于J2000平赤道和平春分点的3个角度值,其计算定义(公式中的角度单位为度)如下：
ζ＝(2306.2181T+0.30188T2+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)：
其中
其中章动角分量φi和系数Ai,Bi,Ci,Di根据IAU1980章动数据表进行计算。
其中历元T时刻平黄赤交角为
历元T时刻真黄赤交角
SA604、计算当前时刻t的地球自转转换矩阵R(t)：
其中格林尼治真恒星时其中θGMST＝GMST(JDUT1)为当前世界时UTC1时刻儒略日时间的格林尼治平恒星时。/>
SA605、计算当前时刻t的极移转换矩阵W(t)：
其中极移分量xp和yp从EOP数据中获取。
SA606、t时刻ECI坐标下的坐标向量rECI转化为ECEF坐标下的坐标向量rECEF：rECEF＝[W(t)′][R(t)′][N(t)′][P(t)′]rECI。
步骤201中，太阳同步圆回归轨道半长轴迭代计算函数OnComputeAa模块流程图如图18所示，具体包括如下步骤：
S501、根据界面输入的轨道高度下限和上限，调用太阳同步圆回归轨道回归系数计算OnComputeSunOrbitBackCoeff模块计算太阳同步圆回归轨道回归系数的上限Qmax和下限Qmin,当前根据界面输入得到的回归系数为m_SatBackCoeff。
S502、判断若Qmin>10000.0)||(Qmax>10000.0成立，则输出无效，否则进入S503。
S503、判断若m_SatBackCoeff<Qmin成立，则输出无效，否则进入S504。
S504、判断若m_SatBackCoeff>Qmax成立，则输出无效，否则进入S505。
S505、根据回归系数m_SatBackCoeff调用太阳同步圆回归轨道半长轴迭代计算OnSunOrbitAa模块迭代计算太阳同步圆回归轨道的半长轴Aa。
S506、判断Aa<0.0是否成立，若是输出无效，否则进入S507。
S507、根据计算获得的太阳同步圆回归轨道半长轴Aa调用根据轨道半长轴计算太阳同步圆轨道倾角函数OnComputeSunOrbitAngle模块迭代计算太阳同步圆回归轨道的轨道倾角IiAngle,并将Aa和IiAngle参数存放在SunOrbit数组中，输出有效。
S505中，其中太阳同步圆回归轨道半长轴迭代计算OnSunOrbitAa模块流程图如图19所示，具体为：
S801、根据从界面输入的轨道入轨点高度上限和下限分别计算太阳同步圆轨道半长轴上限aU和下限aL。
根据轨道半长轴下限和上限，调用太阳同步圆回归轨道回归系数计算OnComputeSunOrbitBackCoeff模块分别计算太阳同步圆回归轨道回归系数的上限Qmax＝OnComputeSunOrbitBackCoeff(aL)和下限Qmin＝OnComputeSunOrbitBackCoeff(aU)。
当前根据界面输入得到的回归系数为m_SatBackCoeff。
S802、判断若(Qmax>100000.0)||Qmin>100000.0)成立，则返回无效值，否则进入S803。
S803、计算第一差值tempValue1＝Qmax-SatBackCoeff；以及第二差值tempValue2＝Qmin-SatBackCoeff。
S804、判断(tempValue1×tempValue2)>0.0是否成立，若是则返回无效值，否则进入S805。
S805、轨道半长轴迭代计算初始化:迭代次数index初始化为0,差值tempValue1和tempValue2。
S806、判断index<30是否成立，若是则进入S807，否则输出有效值aM。
S807、轨道半长轴迭代计算主体部分包括如下步骤：
计算轨道半长轴中间值aM＝(aU+aL)/2.0,index自增1。
根据轨道半长轴中间值aM，调用太阳同步圆回归轨道回归系数计算OnComputeSunOrbitBackCoeff模块计算当前轨道升交点赤经导数tempValue。
S808、判断tempValue>10000.0是否成立，若是则输出无效值；否则进入S809。
S809、判断tempValue×tempValue1>0.0是否成立，若是则令tempValue1＝tempValue,aL＝aM；否则令tempValue2＝tempValue,aU＝aM。
返回S806。
其中太阳同步圆回归轨道回归系数计算OnComputeSunOrbitBackCoeff模块流程如图20所示，具体为：
S901、根据输入参数轨道半长轴aa,调用根据轨道半长轴计算太阳同步圆轨道倾角函数OnComputeSunOrbitAngle模块计算太阳同步圆回归轨道的轨道倾角IiAngle。
S902、判断IiAngle>10000.0是否成立，若是则返回无效值，否则进入S903。
S903、根据输入参数轨道半长轴aa和前面获得的有效轨道倾角IiAngle，调用圆轨道回归系数计算OnComputeBackCoeffCircle模块计算太阳同步圆回归轨道回归系数Q，返回太阳同步圆回归轨道回归系数Q。
圆轨道回归系数计算OnComputeBackCoeffCircle模块流程如图21所示，具体为：
根据轨道半长轴aI、轨道倾角i、地球平均半径Re、地球自转速度ωE和卫星平均角速度n，计算圆轨道回归系数:
返回圆轨道回归系数fa。
根据轨道半长轴计算太阳同步圆轨道倾角函数OnComputeSunOrbitAngle模块流程如图22所示，具体包括如下步骤：
S1101.设置太阳同步圆轨道倾角上限angleU＝120度和下限angleL＝96度
根据当前模块的输入参数轨道半长轴aa和轨道倾角上下限，调用太阳同步圆轨道升交点赤经导数偏差计算OnComputeOmegaDotSunOrbit模块分别计算当前轨道升交点赤经导数与太阳同步轨道标准升交点赤经导数的差值tempValueU和tempValueL
S1102.判断(tempValueL×tempValueU)>0.0是否成立，若是则返回无效值，否则进入S1103。
S1103.轨道半长轴迭代计算初始化:迭代次数index＝0,差值tempValueL和tempValueU。
S1104.判断index<30是否成立，若是则进入S1105，否则返回angleM。
S1105.轨道倾角迭代计算主体部分为：
计算轨道倾角中间值angleM＝(angleU+angleL)/2.0
根据当前输入参数轨道半长轴aa和轨道倾角中间值中间值angleM，调用太阳同步圆轨道升交点赤经导数偏差计算OnComputeOmegaDotSunOrbit模块分别计算当前轨道升交点赤经导数与太阳同步轨道标准升交点赤经导数的差值tempValue。
index自增1。
S1106.判断(tempValue×tempValueL)>0.0是否成立，若是则令tempValueL＝tempValue,angleL＝angleM，否则，令tempValueU＝tempValue,angleU＝angleM。
返回S1104。
太阳同步圆轨道升交点赤经导数偏差计算OnComputeOmegaDotSunOrbit模块流程如图23所示，具体为：
根据轨道半长轴aI、轨道倾角i、地球平均半径Re、地球自转速度ωE和卫星平均角速度n，计算当前轨道升交点赤经导数与太阳同步圆轨道标准升交点赤经导数的偏差值：
对于本发明卫星轨道规划软件计算获得的卫星轨道方案，本发明将卫星轨道方案参数和目标区域参数带入到STK软件中进行仿真评估。根据STK软件的仿真评估结果，可以判断本发明卫星轨道规划设计和轨道评估计算的准确性。
综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

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种入轨模式轨道倾角遍历范围内所有卫星轨道，调用太阳同步圆回归轨道卫星的发射轨道和运行轨道迭代计算OnIterativeSunCompute模块计算第j次过境目标纬度圈时星下点经度分布,选取轨道倾角遍历范围内过境时刻星下点经度与目标点经度偏差值最小的轨道；
步骤5、判断步骤4中所选取的偏差值最小的轨道对应的偏差值是否满足预设的指标要求，若是则进入步骤6，否则MO自增1，返回步骤3；
步骤6、判断当前卫星载荷类型是否为CCD光学设备，若是则进入步骤7，否则进入步骤10；
步骤7、判断当前轨道卫星过境目标点时刻是否在CCD光学设备的UTC时间观测窗口内，若是则进入步骤10，否则进入步骤8；
步骤8、计算发射延时，对于当前入轨模式轨道倾角遍历范围内的所有卫星轨道，调用太阳同步圆回归轨道卫星的发射轨道和运行轨道迭代计算OnIterativeSunCompute模块计算当前第j次过境目标纬度圈时星下点经度分布,选取遍历范围内过境时刻星下点经度与目标点经度偏差值最小的轨道；
步骤9、判断步骤8中所选取的偏差值最小的轨道是否满足预设的指标要求，若是则进入步骤10，否则返回步骤3；
步骤10、当前第BackupOrbitNum个轨道规划方案满足初步要求，将当前第BackupOrbitNum个轨道规划方案存入在BackupOrbit数组中,并且BackupOrbitNum自增1；
步骤11、太阳同步圆回归轨道规划完成，将当前对应的共BackupOrbitNum个轨道方案存入预先构建的BackupOrbit数组中。
2.如权利要求1所述的方法，其特征在于，所述步骤4中，所述太阳同步圆回归轨道卫星的发射轨道和运行轨道迭代计算OnIterativeSunCompute模块，具体包括如下步骤：
步骤201、根据界面输入的轨道高度上下限和轨道回归系数，调用OnComputeAa模块计算判断当前输入的轨道高度上下限和轨道回归系数是否有效；如果无效,则返回值Flag＝false；如果有效，则返回值Flag＝true，根据轨道回归系数迭代计算轨道半长轴和轨道倾角，结果保存在SunOrbit数组中；
步骤202、判断返回值Flag是否有效，若是则进入步骤203，否则输出false；
步骤203、根据SunOrbit数组中保存的轨道半长轴和轨道倾角，以及指定的发射区域中心位置参数，调用发射点位置遍历选取进行太阳同步圆回归轨道设计OnComputeOrbitMinValueSunRepeatOrbit模块在发射区域中心位置附近遍历选取发射点位置，计算对应不同发射点位置太阳同步圆回归轨道星下点轨迹第j次过境目标纬度圈时的星下点地心经度与目标点地心经度的偏差分布，选取偏差值最小的轨道设计进行门限判决,并返回有效轨道设计的轨道倾角IiAngle，其中，j从1到16,为模块输入参数。
3.如权利要求2所述的方法，其特征在于，所述步骤203中，发射点位置遍历选取进行太阳同步圆回归轨道设计OnComputeOrbitMinValueSunRepeatOrbit模块，具体包括如下步骤：
步骤301、太阳同步圆回归轨道设计循环计算初始化，从SunOrbit数组中读取当前的轨道半长轴aa和轨道倾角IiAngleDegree；设定当前轨道索引SunOrbitInitialNum初值为0；在发射区域中心位置附近均匀遍历选取发射点位置；
步骤302、调用太阳同步圆轨道卫星的发射轨道和运行轨道计算OnComputeOrbitValueSunOrbit模块计算对应不同发射点位置太阳同步圆回归轨道卫星入轨后星下点16次过境目标纬度圈时星下点地心经度与目标点地心经度的偏差分布；当前轨道索引SunOrbitInitialNum自增1；
步骤303、判断判断当前发射点位置遍历是否完成，若是则进入步骤304，否则选取下一发射点位置，返回步骤302；
步骤304、对于遍历范围内SunOrbitInitialNum条卫星轨道星下点轨迹第j次过境目标纬度圈时的星下点地心经度与目标点地心经度的偏差分布进行统计，选取偏差值最小的轨道方案,并对该最小偏差值minValue进行门限threshold判决，其中，j从1到16,为模块输入参数；
步骤305、对最小偏差值minValue进行门限threshold判决，若minValue<＝threshold，则进入步骤306；否则返回无效值；
步骤306、minValue对应的轨道方案保存在minTimeOrbit数组中，并返回有效轨道倾角值。
4.如权利要求1～3中任意一项所述的方法，其特征在于，所述太阳同步圆轨道卫星的发射轨道和运行轨道计算OnComputeOrbitValueSunOrbit模块，具体包括如下步骤：
S401、卫星发射轨道参数初始化计算部分：
①.发射点地心经纬度为发射轨道地心角βL和发射轨道倾角i为模块输入参数；
②.发射时刻卫星在发射轨道上真近点角fLradian和入轨时刻卫星在发射轨道上真近点角fIradian的计算；
③.发射轨道偏心率eL、发射轨道半通径p和发射轨道半长轴aL的计算
④.发射点偏近点角ELangleRadian和发射点平近点角MLangleRadian的计算
⑤.根据上述参数计算卫星入轨时刻tI的计算
⑥.发射时刻GMST格林尼治平恒星时GL的计算
⑦.发射时刻发射点在ECI地心惯性坐标系中赤经为：
αL＝OnLimitPiRadian(GL+λL)，其中角度约束函数OnLimitPiRadian(x)将角度x约束到[-π,π)范围内；
S402、若为下降入轨模式，发射时刻卫星纬度辐角计算如下:
轨道升交点赤经:Ω＝OnLimit2PiRadian(αL-tan-1(tan(uL)×cos(i))+π)；其中角度约束函数OnLimit2PiRadian(x)将角度x约束到[0，2π)范围内；
若为上升入轨模式，发射时刻卫星纬度辐角计算如下:
轨道升交点赤经:Ω＝OnLimit2PiRadian(αL-tan-1(tan(uL)×cos(i))+0)；
S403、卫星运行轨道参数初始化部分:
入轨时刻卫星纬度辐角uI＝uL+βL,目标点地心经纬度卫星平均角速度n,地球自转转速ωE；
入轨点地心纬度为
S404、若为下降入轨模式，则分如下4种情况计算卫星入轨后第1次和第2次过境目标纬度圈的卫星纬度辐角值uT1和uT2；目标点地心纬度为入轨点地心维度/>
第1种情况.目标点在北半球,目标点地心纬度小于等于入轨点地心纬度；
第2种情况.目标点在南半球,目标点地心纬度小于等于入轨点地心纬度；
第3种情况.目标点在北半球,目标点地心纬度大于入轨点地心纬度；
第4种情况.目标点在南半球,目标点地心纬度大于入轨点地心纬度；
若即目标点地心纬度小于等于入轨点地心纬度，则对于入轨后第1次过境，包括第1种情况和第2种情况，如果OnLimit2PiRadian(uT1)<π则属于第1种情况，第一中间变量为Δu＝π-uT1；第一中间变量为utemp＝π+2Δu；否则属于第2种情况：Δu＝uT1-π；utemp＝π-2Δu；则入轨后第2次过境uT2＝uT1+utemp；
若即目标点地心纬度大于入轨点地心纬度，则对于入轨后第1次过境，包括第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种情况.目标点在南半球,目标点纬度小于入轨点纬度；
若即目标点地心纬度大于等于入轨点地心纬度，且入轨点在下降段，即UI≥π/2，则对于入轨后第1次过境，/>tempRadian＝OnLimit2PiRadian(uT1)；
若tempRadian＜π属于第5种情况，Δu＝tempRadian；utemp＝π-2Δu；否则属于第6种情况，则Δu＝fabs(tempRadian-2π)；utemp＝π+2Δu；
若即目标点地心纬度大于等于入轨点地心纬度，且入轨点在上升段，即UI＜π/2，则对于入轨后第1次过境，/>tempRadian＝OnLimit2PiRadian(uT1)；
若tempRadian＜π属于第7种情况，Δu＝uT1；utemp＝π-2Δu；否则属于第8种情况，第8种情况不存在；
则入轨后第2次过境uT2＝uT1+utemp；
若即目标点地心纬度小于入轨点地心纬度，则对于入轨后第1次过境，tempRadian＝OnLimit2PiRadian(uT1)；
若tempRadian＜π属于第9种情况，utemp＝π+2Δu；否则属于第10种情况，/>utemp＝π+2Δu；
则入轨后第2次过境uT2＝uT1+utemp；
S405、卫星运行轨道星下点过境目标纬度圈时刻地心经度计算：
1).计算卫星运行轨道周期tΩ、升交点赤经导数和近地点辐角导数/>
2).计算卫星入轨后16次过境目标纬度圈的时刻:
t1＝(uT1-uI)/n,t2＝(uT2-uI)/n,t2m+1＝t1+mtΩ,t2m+2＝t2+mtΩ,m＝1,2,...,7
m为次数；
3).调用圆形轨道卫星的星下点轨迹计算OnComputeValueLonGCcircle模块计算入轨后第1次和第2次过境目标纬度圈时星下点地心经度LonT(1)和LonT(2)；
4).计算入轨后第3次到第16次卫星星下点轨迹过境目标纬度圈时刻的地心经度:
5).计算卫星星下点轨迹16次过境目标纬度圈时刻地心经度与目标地心经度的偏差。
5.如权利要求4所述的方法，其特征在于，所述调用圆形轨道卫星的星下点轨迹计算OnComputeValueLonGCcircle模块，具体包括如下步骤：
SA1.卫星半长轴aI和当前时刻纬度辐角uT,卫星在地心轨道坐标系下的位置矢量r＝[aI cos(uT) aI sin(uT) 0]′；
SA2.根据相对入轨时刻时间差Δt和升交点赤经导数计算当前时刻卫星运行轨道的升交点赤经/>其中ΩI为卫星入轨时刻的升交点赤经；
SA3.根据相对入轨时刻时间差Δt和近地点辐角导数计算当前时刻卫星运行轨道的近地点辐角/>其中ωI为卫星入轨时刻的近地点辐角；
SA4.构建地心轨道坐标系到地心惯性坐标系的转换矩阵,其中i为轨道倾角
SA5.对地心轨道坐标系下的位置矢量进行坐标转转获取地心惯性ECI坐标系下的位置矢量peci
peci＝Τm×r
SA6.查询地球运动参数EOP文件,调用ECI坐标系转ECEF坐标系OnConvertECItoECEF模块将ECI位置矢量peci转换为ECEF位置矢量pecef；
SA7.根据ECEF位置矢量pecef计算当前时刻圆形轨道卫星的星下点轨迹地心经度Lon。
6.如权利要求5所述的方法，其特征在于，所述SA 6.ECI坐标系转ECEF坐标系OnConvertECItoECEF模块，具体为：
SA601、从EOP文件中读取当前时刻t的6个关键地球运动参数:
(TAI-UTC)差值时间dat,单位秒,其中UTC为当前时刻t的世界协调时,原子时TAI＝UTC+dat
极移X分量xp，单位为弧度；
极移Y分量yp，单位为弧度；
(UTC1-UTC)差值时间dut,单位秒,其中UTC1为消除了极移影响后得到的世界时；
赤经章动ΔΨ修正量δΔΨ，单位为弧度；
交角章动Δε修正量δΔε，单位为弧度；
SA602、计算当前时刻t的岁差转换矩阵P(t)：
其中角度变量值ζ、θ和z为历元时刻T的平赤道面和平春分点相对于J2000平赤道和平春分点的3个角度值,其计算定义如下，公式中的角度单位为度：
ζ＝(2306.2181T+0.30188T2+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)：
其中
其中章动角分量φi和系数Ai,Bi,Ci,Di根据IAU1980章动数据表进行计算；
其中历元T时刻平黄赤交角为
历元T时刻真黄赤交角
SA604、计算当前时刻t的地球自转转换矩阵R(t)：
其中格林尼治真恒星时其中θGMST＝GMST(JDUT1)为当前世界时UTC1时刻儒略日时间的格林尼治平恒星时；
SA605、计算当前时刻t的极移转换矩阵W(t)：
其中极移分量xp和yp从EOP数据中获取；
SA606、t时刻ECI坐标下的坐标向量rECI转化为ECEF坐标下的坐标向量rECEF：rECEF＝[W(t)′][R(t)′][N(t)′][P(t)′]rECI。
7.如权利要求6所述的方法，其特征在于，所述步骤201中，所述OnComputeAa模块，具体包括如下步骤：
S501、根据界面输入的轨道高度下限和上限，调用太阳同步圆回归轨道回归系数计算OnComputeSunOrbitBackCoeff模块计算太阳同步圆回归轨道回归系数的上限Qmax和下限Qmin,当前根据界面输入得到的回归系数为m_SatBackCoeff；
S502、判断若(Qmin>10000.0)||(Qmax>10000.0)成立，则输出无效，否则进入S503；
S503、判断若m_SatBackCoeff<Qmin成立，则输出无效，否则进入S504；
S504、判断若m_SatBackCoeff>Qmax成立，则输出无效，否则进入S505；
S505、根据回归系数m_SatBackCoeff调用OnSunOrbitAa模块迭代计算太阳同步圆回归轨道的半长轴Aa；
S506、判断Aa<0.0是否成立，若是输出无效，否则进入S507；
S507、根据计算获得的太阳同步圆回归轨道半长轴Aa调用根据轨道半长轴计算太阳同步圆轨道倾角函数OnComputeSunOrbitAngle模块迭代计算太阳同步圆回归轨道的轨道倾角IiAngle,并将Aa和IiAngle参数存放在SunOrbit数组中，输出有效。
8.如权利要求7所述的方法，其特征在于，S505中，所述太阳同步圆回归轨道半长轴迭代计算OnSunOrbitAa模块，具体为：
S801、根据从界面输入的轨道入轨点高度上限和下限分别计算太阳同步圆轨道半长轴上限aU和下限aL；
根据轨道半长轴下限和上限，调用太阳同步圆回归轨道回归系数计算OnComputeSunOrbitBackCoeff模块分别计算太阳同步圆回归轨道回归系数的上限Qmax＝OnComputeSunOrbitBackCoeff(aL)和下限Qmin＝OnComputeSunOrbitBackCoeff(aU)；
当前根据界面输入得到的回归系数为m_SatBackCoeff；
S802、判断若(Qmax>100000.0)||(Qmin>100000.0)成立，则返回无效值，否则进入S803；
S803、计算第一差值tempValue1＝Qmax-SatBackCoeff；以及第二差值tempValue2＝Qmin-SatBackCoeff；
S804、判断(tempValue1×tempValue2)>0.0是否成立，若是则返回无效值，否则进入S805；
S805、轨道半长轴迭代计算初始化:迭代次数index初始化为0,差值tempValue1和tempValue2；
S806、判断index<30是否成立，若是则进入S807，否则输出有效值aM；
S807、轨道半长轴迭代计算主体部分包括如下步骤：
计算轨道半长轴中间值aM＝(aU+aL)/2.0,index自增1；
根据轨道半长轴中间值aM，调用太阳同步圆回归轨道回归系数计算OnComputeSunOrbitBackCoeff模块计算当前轨道升交点赤经导数tempValue；
S808、判断tempValue>10000.0是否成立，若是则输出无效值；否则进入S809；
S809、判断tempValue×tempValue1>0.0是否成立，若是则令tempValue1＝tempValue,aL＝aM；否则令tempValue2＝tempValue,aU＝aM；
返回S806。
9.如权利要求8所述的方法，其特征在于，所述OnComputeSunOrbitBackCoeff模块具体为：
S901、根据输入参数轨道半长轴aa,调用根据轨道半长轴计算太阳同步圆轨道倾角函数OnComputeSunOrbitAngle模块计算太阳同步圆回归轨道的轨道倾角IiAngle；
S902、判断IiAngle>10000.0是否成立，若是则返回无效值，否则进入S903；
S903、根据输入参数轨道半长轴aa和前面获得的有效轨道倾角IiAngle，调用太阳同步圆回归轨道回归系数计算OnComputeBackCoeffCircle模块计算太阳同步圆回归轨道回归系数Q，返回太阳同步圆回归轨道回归系数Q。
10.如权利要求9所述的方法，其特征在于，所述太阳同步圆回归轨道回归系数计算OnComputeBackCoeffCircle模块，具体为：
根据轨道半长轴aI、轨道倾角i、地球平均半径Re、地球自转速度ωE和卫星平均角速度n，计算圆轨道回归系数:
返回圆轨道回归系数fa；
所述根据轨道半长轴计算太阳同步圆轨道倾角函数OnComputeSunOrbitAngle模块，具体包括如下步骤：
S1101.设置太阳同步圆轨道倾角上限angleU＝120度和下限angleL＝96度
根据当前模块的输入参数轨道半长轴aa和轨道倾角上下限，调用太阳同步圆轨道升交点赤经导数偏差计算OnComputeOmegaDotSunOrbit模块分别计算当前轨道升交点赤经导数与太阳同步轨道标准升交点赤经导数的差值tempValueU和tempValueL
S1102.判断(tempValueL×tempValueU)>0.0是否成立，若是则返回无效值，否则进入S1103；
S1103.轨道半长轴迭代计算初始化:迭代次数index＝0,差值tempValueL和tempValueU；
S1104.判断index<30是否成立，若是则进入S1105，否则返回angleM；
S1105.轨道倾角迭代计算主体部分为：
计算轨道倾角中间值angleM＝(angleU+angleL)/2.0
根据当前输入参数轨道半长轴aa和轨道倾角中间值中间值angleM，调用太阳同步圆轨道升交点赤经导数偏差计算OnComputeOmegaDotSunOrbit模块分别计算当前轨道升交点赤经导数与太阳同步轨道标准升交点赤经导数的差值tempValue；
index自增1；
S1106.判断(tempValue×tempValueL)>0.0是否成立，若是则令tempValueL＝tempValue,angleL＝angleM，否则，令tempValueU＝tempValue,angleU＝angleM；
返回S1104；
所述太阳同步圆轨道升交点赤经导数偏差计算OnComputeOmegaDotSunOrbit模块，具体为：
根据轨道半长轴aI、轨道倾角i、地球平均半径Re、地球自转速度ωE和卫星平均角速度n，计算当前轨道升交点赤经导数与太阳同步圆轨道标准升交点赤经导数的偏差值：
CN202011592823.XA 2020-12-29 2020-12-29 一种太阳同步圆回归轨道的轨道规划方法 Active CN112783183B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011592823.XA CN112783183B (zh) 2020-12-29 2020-12-29 一种太阳同步圆回归轨道的轨道规划方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011592823.XA CN112783183B (zh) 2020-12-29 2020-12-29 一种太阳同步圆回归轨道的轨道规划方法

Publications (2)

Publication Number Publication Date
CN112783183A CN112783183A (zh) 2021-05-11
CN112783183B true CN112783183B (zh) 2023-12-19

Family

ID=75753146

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011592823.XA Active CN112783183B (zh) 2020-12-29 2020-12-29 一种太阳同步圆回归轨道的轨道规划方法

Country Status (1)

Country Link
CN (1) CN112783183B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114415716B (zh) * 2021-12-17 2023-02-28 哈尔滨工业大学 一种维持星座构型的方法、装置及介质
CN114021068B (zh) * 2022-01-05 2022-03-11 成都国星宇航科技有限公司 太阳同步圆轨道卫星受晒因子计算方法、装置及电子设备

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110203422A (zh) * 2019-05-31 2019-09-06 中国人民解放军63729部队 针对面目标区域探测的快速响应卫星轨道设计方法
CN111443722A (zh) * 2020-03-23 2020-07-24 上海航天控制技术研究所 一种编队卫星定时段自主保持方法
CN111731513A (zh) * 2020-06-15 2020-10-02 航天东方红卫星有限公司 一种基于单脉冲轨控的高精度引力场中回归轨道维持方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
FR3047813B1 (fr) * 2016-02-16 2019-11-08 Airbus Defence And Space Sas Procede de commande de guidage d'attitude d'un satellite, satellite, pluralites de satellites et programme d'ordinateur associe
US11643225B2 (en) * 2017-07-21 2023-05-09 The Aerospace Corporation Interlocking, reconfigurable, reconstitutable, reformable cell-based space system

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110203422A (zh) * 2019-05-31 2019-09-06 中国人民解放军63729部队 针对面目标区域探测的快速响应卫星轨道设计方法
CN111443722A (zh) * 2020-03-23 2020-07-24 上海航天控制技术研究所 一种编队卫星定时段自主保持方法
CN111731513A (zh) * 2020-06-15 2020-10-02 航天东方红卫星有限公司 一种基于单脉冲轨控的高精度引力场中回归轨道维持方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
太阳同步卫星的轨道设计;陈洁, 汤国建;上海航天(第03期);全文 *

Also Published As

Publication number Publication date
CN112783183A (zh) 2021-05-11

Similar Documents

Publication Publication Date Title
CN110203422B (zh) 针对面目标区域探测的快速响应卫星轨道设计方法
CN112629543A (zh) 一种大椭圆轨道及小倾角圆轨道的轨道规划方法
WO2017113567A1 (zh) 火星探测器自主导航方法
CN105905317B (zh) 一种卫星对日定向控制系统及其控制方法
CN107450582B (zh) 一种基于星上实时规划的相控阵数传引导控制方法
CN112649006A (zh) 一种太阳同步圆轨道的轨道规划方法
CN112783183B (zh) 一种太阳同步圆回归轨道的轨道规划方法
CN108508918B (zh) 一种静轨遥感卫星数传天线高精度实时对地指向控制方法
CN111427002B (zh) 地面测控天线指向卫星的方位角计算方法
CN109946728B (zh) 一种适用于卫星用户站数字跟踪接收机的程序跟踪方法
CN102175241A (zh) 一种火星探测器巡航段自主天文导航方法
CN110146093A (zh) 双体小行星探测自主协同光学导航方法
Anderson et al. Gravity results from Pioneer 10 Doppler data
CN112713922A (zh) 一种多波束通讯卫星的可见性快速预报算法
Fujita et al. Attitude maneuvering sequence design of high-precision ground target tracking control for multispectral Earth observations
CN106777580B (zh) 近地倾斜轨道发射窗口快速设计方法
CN109269508A (zh) 一种卫星相对小行星视觉自主导航方法
CN113310496A (zh) 一种确定月地转移轨道的方法及装置
McAdams et al. MESSENGER mission design and navigation
Abrahamson et al. Dawn orbit determination team: Trajectory modeling and reconstruction processes at Vesta
CN106767824B (zh) 一种计算双探测器在地外天体表面相对位置的方法
CN116125503A (zh) 一种高精度卫星轨道确定及预报算法
Liu et al. Microsatellite autonomous orbit propagation method based on SGP4 model and GPS data
CN111475767B (zh) 考虑地球自转影响的最小能量弹道严格构造方法
Duxbury A spacecraft-based navigation instrument for outer planet missions

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