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

技术领域

本发明涉及轨道规划技术领域，具体涉及一种太阳同步圆回归轨道的轨道规划方法。

背景技术

常规对地观测卫星轨道规划与应急卫星发射轨道规划的区别在于任务设计理念不同。常规对地观测任务首先发射对地观测卫星，在给出具体观测目标后，通过轨道机动将卫星送入观测轨道对目标进行观测。应急卫星发射任务的卫星轨道规划受到空间任务需求的牵引,通常根据侦察目标任务需求进行卫星轨道规划,权衡卫星发射轨道和运行轨道的各项参数,使之在满足应急发射的前提下,可以快速对目标区域进行有效探测。

对于应急卫星发射任务而言，需要综合考虑卫星的发射准备、发射等待和发射入轨过程中的约束，通过多种手段减小卫星发射和入轨各阶段的时间，从而实现对用户需求的快速响应。将卫星送入高轨轨道需要采用大型运载，而大型运载的发射准备时间较长。同时，高轨卫星的发射入轨方式比低轨卫星复杂，卫星入轨时间也较长。当前应急卫星发射通常采用低轨轨道。应急卫星发射过程中在还需要根据发射态势情况计算确定发射窗口，发射窗口的时间和大小由具体任务确定。在应急卫星发射任务中，应根据发射、测控和光照等约束条件具体分析计算其发射窗口。

传统应急卫星发射轨道规划主要是针对入轨后的卫星轨道进行规划设计，通常都没有综合考虑发射轨道模型、轨道回归特性、光学载荷的光照条件约束和地球复杂运动模型等多种因素的影响。

因此，如何克服传统应急卫星发射轨道设计的局限性，结合卫星发射轨道模型、轨道回归特性约束、目标区域光照条件约束和精确地球运动模型(本设计考虑了岁差、章动、自转和极移等地球复杂运动模型)进行轨道方案设计，实现对目标区域的快速有效探测，是目前亟待解决的问题。

发明内容

有鉴于此，本发明提供了一种太阳同步圆回归轨道的轨道规划方法，能够结合卫星发射轨道模型、轨道回归特性约束、目标区域光照条件约束和精确地球运动模型进行轨道规划，实现了对目标区域的快速有效探测。

为达到上述目的，本发明的技术方案包括如下步骤：

步骤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、卫星发射轨道参数初始化计算部分：

①.发射点地心经纬度为

②.发射时刻卫星在发射轨道上真近点角fLradian和入轨时刻卫星在发射轨道上真近点角fIradian的计算；

③.发射轨道偏心率e

④.发射点偏近点角ELangleRadian和发射点平近点角MLangleRadian的计算

⑤.根据上述参数计算卫星入轨时刻tI的计算

⑥.发射时刻GMST格林尼治平恒星时G

⑦.发射时刻发射点在ECI地心惯性坐标系中赤经为：

α

S402、若为下降入轨模式，发射时刻卫星纬度辐角计算如下:

轨道升交点赤经:Ω＝OnLimit2PiRadian(α

若为上升入轨模式，发射时刻卫星纬度辐角计算如下:

轨道升交点赤经:Ω＝OnLimit2PiRadian(α

S403、卫星运行轨道参数初始化部分:

1.入轨时刻卫星纬度辐角u

2.入轨点地心纬度为

S404、若为下降入轨模式，则分如下4种情况计算卫星入轨后第1次和第2次过境目标纬度圈的卫星纬度辐角值u

第1种情况.目标点在北半球,目标点地心纬度小于等于入轨点地心纬度；

第2种情况.目标点在南半球,目标点地心纬度小于等于入轨点地心纬度；

第3种情况.目标点在北半球,目标点地心纬度大于入轨点地心纬度；

第4种情况.目标点在南半球,目标点地心纬度大于入轨点地心纬度；

若

若

如果OnLimit2PiRadian(u

若为上升入轨模式，分如下6种情况计算卫星入轨后第1次和第2次过境目标纬度圈的卫星纬度辐角值u

第5种情况.目标点在北半球,目标点纬度大于入轨点纬度,入轨点在下降段；

第6种情况.目标点在南半球,目标点纬度大于入轨点纬度,入轨点在下降段；

第7种情况.目标点在北半球,目标点纬度大于入轨点纬度,入轨点在上升段；

第8种情况.目标点在南半球,目标点纬度大于入轨点纬度,入轨点在上升段；

第9种情况.目标点在北半球,目标点纬度小于入轨点纬度；

第10种情况.目标点在南半球,目标点纬度小于入轨点纬度；

若

若tempRadian＜π属于第5种情况，Δu＝tempRadian；u

若

若tempRadian＜π属于第7种情况，Δu＝u

则入轨后第2次过境u

若

若tempRadian＜π属于第9种情况，

则入轨后第2次过境u

S405、卫星运行轨道星下点过境目标纬度圈时刻地心经度计算：

1).计算卫星运行轨道周期t

2).计算卫星入轨后16次过境目标纬度圈的时刻:

t

3).调用圆形轨道卫星的星下点轨迹计算OnComputeValueLonGCcircle模块计算入轨后第1次和第2次过境目标纬度圈时星下点地心经度Lon

4).计算入轨后第3次到第16次卫星星下点轨迹过境目标纬度圈时刻的地心经度:

5).计算卫星星下点轨迹16次过境目标纬度圈时刻地心经度与目标地心经度的偏差。

进一步地，调用圆形轨道卫星的星下点轨迹计算OnComputeValueLonGCcircle模块，具体包括如下步骤：

SA1.卫星半长轴a

SA2.根据相对入轨时刻时间差Δt和升交点赤经导数

SA3.根据相对入轨时刻时间差Δt和近地点辐角导数

SA4.构建地心轨道坐标系到地心惯性坐标系的转换矩阵,其中i为轨道倾角

SA 5.对地心轨道坐标系下的位置矢量进行坐标转转获取地心惯性ECI坐标系下的位置矢量p

p

SA 6.查询地球运动参数EOP文件,调用ECI坐标系转ECEF坐标系OnConvertECItoECEF模块将ECI位置矢量p

SA7.根据ECEF位置矢量p

进一步地，SA 6.ECI坐标系转ECEF坐标系OnConvertECItoECEF模块，具体为：

SA601、从EOP文件中读取当前时刻t的6个关键地球运动参数:

1.(TAI-UTC)差值时间dat,单位秒,其中UTC为当前时刻t的世界协调时,原子时TAI＝UTC+dat

2.极移X分量x

3.极移Y分量y

4.(UTC1-UTC)差值时间dut,单位秒,其中UTC1为消除了极移影响后得到的世界时；

5.赤经章动ΔΨ修正量δΔΨ，单位为弧度；

6.交角章动Δε修正量δΔε，单位为弧度；

SA602、计算当前时刻t的岁差转换矩阵P(t)：

其中角度变量值ζ、θ和z为历元时刻T的平赤道面和平春分点相对于J2000平赤道和平春分点的3个角度值,其计算定义(公式中的角度单位为度)如下：

ζ＝(2306.2181T+0.30188T

θ＝(2004.3109T-0.42665T

z＝(2306.2181T+1.09468T

定义历元时刻T指的是地球时从J2000的TT时刻起算的儒略世纪数,T＝(JD

SA603、计算当前时刻t的章动转换矩阵N(t)：

其中

其中章动角分量φ

其中历元T时刻平黄赤交角为

历元T时刻真黄赤交角

SA604、计算当前时刻t的地球自转转换矩阵R(t)：

其中格林尼治真恒星时

SA605、计算当前时刻t的极移转换矩阵W(t)：

其中极移分量x

SA606、t时刻ECI坐标下的坐标向量r

进一步地，步骤201中，太阳同步圆回归轨道半长轴迭代计算函数OnComputeAa模块，具体包括如下步骤：

S501、根据界面输入的轨道高度下限和上限，调用太阳同步圆回归轨道回归系数计算OnComputeSunOrbitBackCoeff模块计算太阳同步圆回归轨道回归系数的上限Qmax和下限Qmin,当前根据界面输入得到的回归系数为m_SatBackCoeff；

S502、判断若Qmin>10000.0)||(Qmax>10000.0成立，则输出无效，否则进入S503；

S503、判断若m_SatBackCoeff

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模块，具体为：

根据轨道半长轴a

返回圆轨道回归系数f

根据轨道半长轴计算太阳同步圆轨道倾角函数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模块，具体为：

根据轨道半长轴a

有益效果：

本发明提供了一种太阳同步圆回归轨道的轨道规划方法，太阳同步圆回归轨道设计要满足太阳同步轨道升交点赤经导数约束和回归系数约束这两个指标要求。根据这两个指标要求，可以通过迭代计算获得太阳同步圆回归轨道的轨道半长轴和轨道倾角这两个重要参数。根据计算获得的轨道半长轴和轨道倾角，以及指定发射区域中心位置等参数，在发射区域中心位置附近遍历选取发射点位置，计算对应不同发射点位置太阳同步圆回归轨道星下点多次过境目标纬度圈时的星下点地心经度与目标点地心经度的偏差分布，选取偏差值最小的轨道设计进行门限判决，并保存有效轨道设计。此外在卫星轨道设计流程中考虑了太阳光照参数对于光学载荷卫星轨道设计的实施要求,并根据轨道参数计算在太阳高度角约束条件下的卫星发射窗口。因此本发明能够结合卫星发射轨道模型、轨道回归特性约束、目标区域光照条件约束和精确地球运动模型进行轨道规划，实现了对目标区域的快速有效探测。

附图说明

图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所示，图中，R

近地快速覆盖轨道的发射入轨既要满足快速性的要求又要尽量节省能量。卫星发射方式包括共面发射和非共面发射，由于非共面发射比共面发射需要消耗更多的能量，近地快速覆盖轨道采用共面发射。若发射轨道在入轨点与运行轨道相切，该轨道称为“单共切”轨道。当发射速度倾角Θ给定时，单共切轨道为最佳发射轨道。近地快速覆盖轨道通常采用共面发射，单共切轨道入轨。在共面发射条件下，卫星发射轨道和运行轨道的轨道倾角i(定义为用地轴的北极方向与轨道平面的正法线方向之间的夹角度量，其中0≤i≤π，i＝0或π赤道轨道，i＝π/2为极地轨道)相等，发射轨道和运行轨道的升交点赤经Ω(定义为从北极上看从春分点逆时针转到升交点位置的转角，0≤Ω＜2π)也相等。发射轨道为椭圆轨道，通常可以采用卫星发射时刻t

其中，a

其中，a

卫星采用单共切轨道入轨，发射轨道与地球相交于发射点，与运行轨道相切于远地点，该点为卫星入轨点。发射点在发射轨道上的真近点角f

p＝a

发射轨道偏心率e

发射点的偏近点角E

发射点的平近点角M

M

发射轨道半长轴a

卫星入轨时刻t

其中μ＝398600.4415km

对于地心经纬度为

当前本轨道规划设计中，应急发射点只考虑在国内的情况(即在北半球)，则

根据卫星发射入轨的四种模式(1逆行下降，2逆行上升，3顺行下降，4顺行上升),可以分别计算发射点卫星纬度辐角u

对于卫星下降入轨模式1和模式3的发射点卫星纬度辐角u

对于卫星上升入轨模式2和模式4的发射点卫星纬度辐角u

定义发射时刻t

α

其中角度约束函数OnLimitPiRadian(x)将角度x约束到[-π，π)范围内。对于地心经纬度为

tan(α

对于卫星下降入轨模式1和模式3的卫星轨道升交点赤经Ω计算如下：

Ω＝OnLimit2PiRadian(α

对于卫星上升入轨模式2和模式4的卫星轨道升交点赤经Ω计算如下：

Ω＝OnLimit2PiRadian(α

其中角度约束函数OnLimit2PiRadian(x)将角度x约束到[0，2π)范围内。

入轨时刻t

如果入轨点在北半球，则卫星纬度辐角u

如果入轨点在南半球，则卫星纬度辐角u

不管入轨点是在北半球还是在南半球，入轨点的地心纬度计算如下所示：

对于卫星下降入轨模式1和模式3，分4种态势计算卫星入轨后第1次和第2次过境目标纬度圈的卫星纬度辐角值u

对于卫星上升入轨模式2和模式4，分如下6种态势计算卫星过境目标纬度圈的卫星纬度辐角值u

获取了卫星入轨后星下点第1次和第2次过境目标纬度圈的卫星纬度辐角值u

为了获取精确的卫星星下点大地经度，需要在ECI惯性坐标系空间位置矢量转ECEF惯性坐标系空间位置矢量的计算过程中，就需要考虑岁差(precession)、章动(nutation)、自转和极移(polar motion)等地球复杂运动模型的影响。只有考虑了这些地球复杂运动模型因素，卫星轨道规划计算结果才可以取得和STK软件仿真结果接近的精度。从1962年1月1号开始到2019年4月20号的地球运动模型参数都可以从STK软件提供的EOP文件中获取,该文件中的数据会定时更新。由于该EOP文件中的数据对于高精度坐标系转换计算非常重要，但在公开的文献中很少介绍这些参数的使用。在卫星轨道规划软件研发过程中，通过查阅文献和开源软件，逐渐掌握了EOP文件中数据的使用方法。

岁差、章动和极移都是引起天极(地球自转轴在天球上的投影)运动的原因，其中岁差与章动引起天极在空间的运动，极移引起其在地球本体内的运动。t时刻ECI坐标下的空间坐标向量r

r

其中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对应的儒略日时间为JD

定义时间自变量T＝(JD

(1).岁差转换矩阵P(t)的计算

岁差描述了地球自转轴指向和春分点的长期变化。岁差转换矩阵P(t)的计算如下:

其中角度变量值ζ,θ和z为历元时刻T的平赤道面和平春分点相对于J2000平赤道和平春分点的3个角度值,其计算定义(公式中的角度单位为度)如下:

ζ＝(2306.2181T+0.30188*T

θ＝(2004.3109T-0.42665T

z＝(2306.2181T+1.09468T

(2).章动转换矩阵N(t)计算

章动描述了赤道和春分点的短期和周期性变化。除了长期岁差运动，地球自转轴指向还存在小的周期性扰动即章动。这是由月球和太阳扭矩的月变化和年变化引起的。章动影响包括赤经章动ΔΨ(nutation in longitude)和交角章动Δε(nutation inobliquity)。

章动转换矩阵N(t)的计算如下:

赤经章动ΔΨ计算如下:

交角章动Δε计算如下:

其中δΔΨ为赤经章动ΔΨ修正量,δΔε为交角章动Δε修正量，可以从EOP文件中获取当前时刻对应的δΔΨ和δΔε值。

其中章动角分量φ

其中整数系数p

L＝134.96298139+(1717915922.6330T+31.310T

M＝357.52772333+(129596581.2240T-0.577T

E＝93.27191028+(1739527263.1370T-13.257T

F＝297.85036306+(1602961601.3280T-6.891T

G＝125.04452222+(-6962890.5390T+7.455T

定义

平黄赤交角

(3).地球自转转换矩阵R(t)计算

自转转换矩阵R(t)是由格林尼治真恒星时θ

其中格林尼治真恒星时θ

其中θ

其中赤经章动ΔΨ和交角章动Δε的计算参见前面公式。

(4).极移转换矩阵W(t)计算

由于地球不是刚体以及其他一些物理因素的影响，造成地球自转轴在地球体内的运动，引起地球极点在地球表面的位置随时间而变化，称为极移。

极移转换矩阵W(t)计算如下:

其中x

太阳同步圆回归轨道设计要满足太阳同步轨道升交点赤经导数约束和回归系数约束这两个指标要求。根据这两个指标要求，可以通过迭代计算获得太阳同步圆回归轨道的轨道半长轴和轨道倾角这两个重要参数。根据计算获得的轨道半长轴和轨道倾角，以及指定发射区域中心位置等参数，在发射区域中心位置附近遍历选取发射点位置，计算对应不同发射点位置太阳同步圆回归轨道星下点多次过境目标纬度圈时的星下点地心经度与目标点地心经度的偏差分布，选取偏差值最小的轨道设计进行门限判决，并保存有效轨道设计。此外在卫星轨道设计流程中考虑了太阳光照参数对于光学载荷卫星轨道设计的实施要求,并根据轨道参数计算在太阳高度角约束条件下的卫星发射窗口。

圆形轨道回归系数f

太阳同步圆轨道的升交点赤经导数约束如下：

其中轨道半长轴a

太阳同步圆回归轨道设计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、卫星发射轨道参数初始化计算部分：

①.发射点地心经纬度为

②.发射时刻卫星在发射轨道上真近点角fLradian和入轨时刻卫星在发射轨道上真近点角fIradian的计算。

③.发射轨道偏心率e

④.发射点偏近点角ELangleRadian和发射点平近点角MLangleRadian的计算。

⑤.根据上述参数计算卫星入轨时刻tI的计算。

⑥.发射时刻GMST格林尼治平恒星时G

⑦.发射时刻发射点在ECI地心惯性坐标系中赤经为：

α

S402、若为下降入轨模式，发射时刻卫星纬度辐角计算如下:

轨道升交点赤经:Ω＝OnLimit2PiRadian(α

若为上升入轨模式，发射时刻卫星纬度辐角计算如下:

轨道升交点赤经:Ω＝OnLimit2PiRadian(α

S403、卫星运行轨道参数初始化部分:

1.入轨时刻卫星纬度辐角u

2.入轨点地心纬度为

S404、若为下降入轨模式，则分如下4种情况计算卫星入轨后第1次和第2次过境目标纬度圈的卫星纬度辐角值u

第1种情况.目标点在北半球,目标点地心纬度小于等于入轨点地心纬度；

第2种情况.目标点在南半球,目标点地心纬度小于等于入轨点地心纬度；

第3种情况.目标点在北半球,目标点地心纬度大于入轨点地心纬度；

第4种情况.目标点在南半球,目标点地心纬度大于入轨点地心纬度；

若

若

如果OnLimit2PiRadian(u

若为上升入轨模式，分如下6种情况计算卫星入轨后第1次和第2次过境目标纬度圈的卫星纬度辐角值u

第5种情况.目标点在北半球,目标点纬度大于入轨点纬度,入轨点在下降段；

第6种情况.目标点在南半球,目标点纬度大于入轨点纬度,入轨点在下降段；

第7种情况.目标点在北半球,目标点纬度大于入轨点纬度,入轨点在上升段；

第8种情况.目标点在南半球,目标点纬度大于入轨点纬度,入轨点在上升段；

第9种情况.目标点在北半球,目标点纬度小于入轨点纬度；

第10种情况.目标点在南半球,目标点纬度小于入轨点纬度；

若

若tempRadian＜π属于第5种情况，Δu＝tempRadian；u

若

若tempRadian＜π属于第7种情况，Δu＝u

则入轨后第2次过境u

若

若tempRadian＜π属于第9种情况，

则入轨后第2次过境u

S405、卫星运行轨道星下点过境目标纬度圈时刻地心经度计算：

1).计算卫星运行轨道周期t

2).计算卫星入轨后16次过境目标纬度圈的时刻:

t

m为次数。

3).调用圆形轨道卫星的星下点轨迹计算OnComputeValueLonGCcircle模块计算入轨后第1次和第2次过境目标纬度圈时星下点地心经度Lon

4).计算入轨后第3次到第16次卫星星下点轨迹过境目标纬度圈时刻的地心经度:

5).计算卫星星下点轨迹16次过境目标纬度圈时刻地心经度与目标地心经度的偏差。

圆形轨道卫星的星下点轨迹计算OnComputeValueLonGCcircle模块流程如图16所示，调用的圆形轨道卫星的星下点轨迹计算OnComputeValueLonGCcircle模块，具体包括如下步骤：

SA1.卫星半长轴a

SA2.根据相对入轨时刻时间差Δt和升交点赤经导数

SA3.根据相对入轨时刻时间差Δt和近地点辐角导数

SA4.构建地心轨道坐标系到地心惯性坐标系的转换矩阵,其中i为轨道倾角

SA 5.对地心轨道坐标系下的位置矢量进行坐标转转获取地心惯性ECI坐标系下的位置矢量p

p

SA6.查询地球运动参数EOP文件,调用ECI坐标系转ECEF坐标系OnConvertECItoECEF模块将ECI位置矢量p

SA7.根据ECEF位置矢量p

SA 6中ECI坐标系转ECEF坐标系OnConvertECItoECEF模块流程如图17所示，具体为：

SA601、从EOP文件中读取当前时刻t的6个关键地球运动参数:

1.(TAI-UTC)差值时间dat,单位秒,其中UTC为当前时刻t的世界协调时,原子时TAI＝UTC+dat

2.极移X分量x

3.极移Y分量y

4.(UTC1-UTC)差值时间dut,单位秒,其中UTC1为消除了极移影响后得到的世界时；

5.赤经章动ΔΨ修正量δΔΨ，单位为弧度；

6.交角章动Δε修正量δΔε，单位为弧度；

SA602、计算当前时刻t的岁差转换矩阵P(t)：

其中角度变量值ζ、θ和z为历元时刻T的平赤道面和平春分点相对于J2000平赤道和平春分点的3个角度值,其计算定义(公式中的角度单位为度)如下：

ζ＝(2306.2181T+0.30188T

θ＝(2004.3109T-0.42665T

z＝(2306.2181T+1.09468T

定义历元时刻T指的是地球时从J2000的TT时刻起算的儒略世纪数,T＝(JD

SA603、计算当前时刻t的章动转换矩阵N(t)：

其中

其中章动角分量φ

其中历元T时刻平黄赤交角为

历元T时刻真黄赤交角

SA604、计算当前时刻t的地球自转转换矩阵R(t)：

其中格林尼治真恒星时

SA605、计算当前时刻t的极移转换矩阵W(t)：

其中极移分量x

SA606、t时刻ECI坐标下的坐标向量r

步骤201中，太阳同步圆回归轨道半长轴迭代计算函数OnComputeAa模块流程图如图18所示，具体包括如下步骤：

S501、根据界面输入的轨道高度下限和上限，调用太阳同步圆回归轨道回归系数计算OnComputeSunOrbitBackCoeff模块计算太阳同步圆回归轨道回归系数的上限Qmax和下限Qmin,当前根据界面输入得到的回归系数为m_SatBackCoeff。

S502、判断若Qmin>10000.0)||(Qmax>10000.0成立，则输出无效，否则进入S503。

S503、判断若m_SatBackCoeff

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所示，具体为：

根据轨道半长轴a

返回圆轨道回归系数f

根据轨道半长轴计算太阳同步圆轨道倾角函数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所示，具体为：

根据轨道半长轴a

对于本发明卫星轨道规划软件计算获得的卫星轨道方案，本发明将卫星轨道方案参数和目标区域参数带入到STK软件中进行仿真评估。根据STK软件的仿真评估结果，可以判断本发明卫星轨道规划设计和轨道评估计算的准确性。

综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。