一种大椭圆轨道及小倾角圆轨道的轨道规划方法与流程

本发明涉及轨道规划技术领域，具体涉及一种大椭圆轨道及小倾角圆轨道的轨道规划方法。

背景技术：

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

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

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

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

技术实现要素：

有鉴于此，本发明提供了一种大椭圆轨道以及小倾角圆轨道的轨道规划方法，能够根据所选取的目标区域信息、轨道回归特性的约束条件、侦察载荷的约束条件、发射部署参数(包括发射点位置、发射部署完成时刻和发射轨道数据等参数)和精确的地球运动模型等数据，完成针对不同类型目标区域探测的大椭圆轨道以及小倾角圆轨道的规划。

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

步骤1、大椭圆轨道和小倾角圆轨道均具有四种入轨模式：即逆行下降、逆行上升、顺行下降、顺行上升共4种入轨模式；根据轨道半长轴的上下限，获得逆行上升和逆行下降的轨道倾角搜索范围、顺行上升和顺行下降的轨道倾角搜索范围。

如果卫星载荷类型为ccd光学设备，则以发射部署完成时刻为起点，计算目标点区域24小时内太阳高度角满足门限要求的utc时间观测窗口。

该轨道规划方法针对4种入轨模式对应轨道的倾角搜索范围内所有卫星轨道，计算卫星轨道星下点轨迹6次过境目标纬度圈时星下点的经度分布，并且选取满足侦察要求的最优轨道。设定初始状态为，mo为入轨模式，mo的取值为1、2、3和4，分别指代逆行下降、逆行上升、顺行下降、顺行上升共4种入轨模式，mo初始值设定为1；j为过境次数，j的初值设定为1，过境总次数为6；backuporbitnum为轨道方案个数，backuporbitnum的初值设定为0。

步骤2、判断j的值是否满足j<＝6，若是则进入步骤3，否则进入步骤11。

步骤3、判断mo的值是否满足mo<＝4；若是则进入步骤4，否则将mo重新设定为1，j自增1，返回步骤2。

步骤4、对于第mo种入轨模式轨道倾角遍历范围内所有卫星轨道，调用遍历范围内卫星大椭圆轨道星下点过境目标纬度圈地心经度分布计算oncomputeorbitminvalueellipse模块计算第j次过境目标纬度圈时星下点经度分布,选取轨道倾角遍历范围内过境时刻星下点经度与目标点经度偏差值最小的轨道。

步骤5、判断步骤4中所选取的偏差值最小的轨道对应的偏差值是否满足预设的指标要求，若是则进入步骤6，否则mo自增1，返回步骤3。

步骤6、判断当前卫星载荷类型是否为ccd光学设备，若是则进入步骤7，否则进入步骤10。

步骤7、判断当前轨道卫星过境目标点时刻是否在ccd光学设备的utc时间观测窗口内，若是则进入步骤10，否则进入步骤8。

步骤8、计算发射延时，对于当前入轨模式轨道倾角遍历范围内的所有卫星轨道，调用遍历范围内卫星大椭圆轨道星下点过境目标纬度圈地心经度分布计算oncomputeorbitminvalueellipse模块计算当前第j次过境目标纬度圈时星下点经度分布,选取遍历范围内过境时刻星下点经度与目标点经度偏差值最小的轨道。

步骤9、判断步骤8中所选取的偏差值最小的轨道是否满足预设的指标要求，若是则进入步骤10，否则返回步骤3。

步骤10、当前第backuporbitnum个轨道规划方案满足初步要求，将当前第backuporbitnum个轨道规划方案存入在backuporbit数组中,并且backuporbitnum自增1；

步骤11、大椭圆或者小倾角圆轨道规划完成，将当前对应的共backuporbitnum个轨道方案存入预先构建的backuporbit数组中。

进一步地，步骤4和步骤8中，遍历范围内卫星大椭圆轨道星下点过境目标纬度圈地心经度分布计算oncomputeorbitminvalueellipse模块，具体包括如下步骤：

步骤101、开始在遍历范围内轨道星下点过境目标纬度圈地心经度分布计算。

步骤102、轨道倾角遍历范围的参数初始化：轨道倾角遍历初始值domaindegreelower,遍历间隔0.01度,遍历次数indexnum,索引index＝1。

步骤103、判断是否满足index<＝indexnum；若是则进入步骤104；否则进入步骤105。

步骤104、对当前轨道倾角为iiangledegree＝domaindegreelower+0.01×(index-1)的卫星发射轨道和卫星运行轨道，调用大椭圆轨道卫星的发射轨道和运行轨道计算oncomputeorbitvalueellipse模块进行轨道星下点轨迹的计算，返回卫星星下点轨迹多次过境目标纬度圈时的星下点地心经度与目标点地心经度的偏差分布数据,并对索引值index进行累积自增1处理，返回步骤103。

步骤105、对于遍历范围内indexnum条卫星轨道星下点轨迹第j次过境目标纬度圈时的星下点地心经度与目标点地心经度的偏差分布进行统计，选取偏差值最小的轨道方案,并对该最小偏差值minvalue进行门限threshold判决；j取值为1～6。

步骤106、若minvalue<＝threshold，则进入步骤107；否则返回无效值。

步骤107、minvalue对应的轨道方案保存在mintimeorbit数组中，并返回有效轨道倾角值。

进一步地，s104中、大椭圆轨道卫星的发射轨道和运行轨道计算oncomputeorbitvalueellipse模块，具体包括如下步骤：

s10421、卫星发射轨道参数初始化计算部分：

1.发射点地心经纬度为发射轨道地心角β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、若为下降入轨模式，包括顺行下降和逆行下降，发射时刻卫星纬度辐角计算如下:

轨道升交点赤经:ω＝onlimit2piradian(αl-tan^-1(tan(ul)×cos(i))+π)；其中角度约束函数onlimit2piradian(x)将角度x约束到[0，2π)范围内；

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

轨道升交点赤经:ω＝onlimit2piradian(αl-tan^-1(tan(ul)×cos(i))+0)；

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

1.入轨时刻卫星纬度辐角ui＝ul+βl,目标点地心经纬度卫星平均角速度n,地球自转转速ωe；

2.入轨点地心纬度为

s10424、若为下降入轨模式，则分如下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种情况.目标点在南半球,目标点纬度小于入轨点纬度；

若即目标点地心纬度大于等于入轨点地心纬度，且入轨点在下降段则对于入轨后第1次过境，tempradian＝onlimit2piradian(ut1)。

若tempradian＜π属于第5种情况，δu＝tempradian；utemp＝π-2δu；否则属于第6种情况，则δu＝fabs(tempradian-2π)；utemp＝π+2δu。

若即目标点地心纬度大于等于入轨点地心纬度，且入轨点在上升段则对于入轨后第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。

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

①.计算卫星运行轨道周期tω、升交点赤经导数和近地点辐角导数和平近点角导数

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

当为大椭圆轨道时，采用如下方式计算卫星入轨后6次过境目标纬度圈的时刻：根据入轨后第1次进入目标纬度圈时刻的真近点角ff1＝ut1-ui先转为偏近点角ee1，再转为平近点角mm1；根据入轨后第2次进入目标纬度圈时刻的真近点角ff2＝ut2-ui，先转为偏近点角ee2，再转为平近点角mm2；卫星入轨后6次过境目标纬度圈的时刻分别为t1～t6。

当为小倾角圆轨道时，采用如下方式计算卫星入轨后6次过境目标纬度圈的时刻：

t1＝ti+(ut1-ui)/n,t2＝ti+(ut2-ui)/n,

t3＝t1+tω,t4＝t2+tω,t5＝t1+2tω,t6＝t2+2tω

③.调用大椭圆轨道卫星的星下点轨迹计算oncomputevaluelongcellipse模块计算入轨后星下点第1次和第2次过境目标纬度圈时地心经度lont1和lont2。

④.计算入轨后星下点第3次、第4次、第5次和第6次过境目标纬度圈时刻的地心经度，分别为lont3、lont4、lont5以及lont6：

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

s10426、将计算结果和中间变量存放到orbitvector数组。

进一步地，大椭圆轨道卫星的星下点轨迹计算oncomputevaluelongcellipse模块，包括如下步骤：

sa1.卫星半长轴ai、偏心率ei和当前时刻t的真近点角ft,计算当前卫星在地心轨道坐标系下的位置矢量r：

当为大椭圆轨道时，

当为小倾角圆轨道时：

r＝[aicos(ut)aisin(ut)0]′

sa2.根据相对入轨时刻时间差δt和升交点赤经导数计算当前时刻卫星运行轨道的升交点赤经其中ωi为卫星入轨时刻的升交点赤经。

sa3.根据相对入轨时刻时间差δt和近地点辐角导数计算当前时刻卫星运行轨道的近地点辐角其中ωi为卫星入轨时刻的近地点辐角。

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

sa5.对地心轨道坐标系下的位置矢量进行坐标转转获取地心惯性eci坐标系下的位置矢量peci

peci＝tm×r

sa6.查询地球运动参数eop文件,调用eci坐标系转ecef坐标系onconvertecitoecef模块将eci位置矢量peci转换为ecef位置矢量pecef。

sa7.根据ecef位置矢量pecef计算当前时刻圆形轨道卫星的星下点轨迹地心经度lon。

sa8.返回当前时刻大椭圆轨道卫星的星下点轨迹地心经度lon。

进一步地，sa6中的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.30188×t²+0.017998t³)/3600.0

θ＝(2004.3109t-0.42665t²-0.041833t³)/3600.0

z＝(2306.2181t+1.09468t²+0.018203t³)/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；

sa607、返回ecef坐标下的坐标向量recef。

有益效果：

本发明提供了一种大椭圆(小倾角圆)轨道的轨道规划方法，根据所选取的目标区域信息、轨道回归特性的约束条件、侦察载荷的约束条件、发射部署参数(包括发射点位置、发射部署完成时刻和发射轨道数据等参数)和精确的地球运动模型等数据，完成针对不同类型目标区域探测的大椭圆轨道的规划设计。大椭圆轨道设计针对有效倾角范围内的角度值以特定间隔(当前设置为0.01度)进行遍历，计算每个倾角角度值i对应卫星轨道星下点轨迹在多次过境目标点纬度圈上的经度值，并且计算目标点地心经度和星下点轨迹过境标点纬度圈时地心经度值之间的差值。选取地心经度差值满足指标要求的卫星轨道作为备选方案，最后以优选准则对备选轨道方案进行选取。此外在卫星轨道设计流程中考虑了太阳光照参数对于光学载荷卫星轨道设计的实施要求,并根据轨道参数计算在太阳高度角约束条件下的卫星发射窗口。本发明克服了传统应急卫星发射轨道设计的局限性，结合了卫星发射轨道模型、轨道回归特性约束、目标区域光照条件约束和精确地球运动模型(本设计考虑了岁差、章动、自转和极移等地球复杂运动模型)针对大椭圆轨道以及小倾角圆轨道进行轨道方案设计，实现对目标区域的快速有效探测。

附图说明

图1为卫星发射和入轨过程示意图；

图2态势1：目标点地心纬度小于入轨点地心纬度，目标点在北半球；

图3态势2：目标点地心纬度小于入轨点地心纬度，目标点在南半球；

图4为态势3：目标点地心纬度大于入轨点地心纬度，目标点在北半球；

图5为态势4：目标点地心纬度大于入轨点地心纬度，目标点在南半球；

图6为态势5：目标点地心纬度大于入轨点地心纬度，入轨点在下降段，目标点在北半球；

图7为态势6：目标点地心纬度大于入轨点地心纬度，入轨点在下降段，目标点在南半球；

图8为态势7：目标点地心纬度大于入轨点地心纬度，入轨点在上升段，目标点在北半球；

图9为态势8(无效)：目标点地心纬度大于入轨点地心纬度，入轨点在上升段，目标点在南半球；

图10为态势9示意图：目标点地心纬度小于入轨点地心纬度，目标点在北半球；

图11为态势10：目标点地心纬度小于入轨点地心纬度，目标点在南半球；

图12为本发明实施例中轨道设计的轨道倾角遍历象限图；

图13为本发明实施例中的大椭圆(小倾角圆)轨道设计oncomputeellipseorbit模块的流程图；

图14为本发明实施例中的遍历范围内卫星大椭圆轨道星下点过境目标纬度圈地心经度分布计算oncomputeorbitminvalueellipse模块的流程图；

图15为本发明实施例中的大椭圆轨道卫星的发射轨道和运行轨道计算oncomputeorbitvalueellipse模块的流程图；

图16为本发明实施例中的大椭圆轨道卫星的星下点轨迹计算oncomputevaluelongcellipse模块的流程图；

图17为本发明实施例中的eci坐标系转ecef坐标系onconvertecitoecef模块的流程图；

具体实施方式

下面结合附图并举实施例，对本发明进行详细描述。

本卫星轨道规划设计所要解决的问题就是要根据所选取的目标区域信息、轨道回归特性的约束条件、侦察载荷的约束条件、发射部署参数(包括发射点位置、发射部署完成时刻和发射轨道数据等参数)和精确的地球运动模型等数据，完成针对不同类型目标区域探测的小倾角圆轨道和大椭圆轨道的规划设计。为了验证本卫星轨道规划设计方法的正确性，本发明还需要解决卫星运行轨道方案对目标区域的侦察窗口评估问题和地面接收站对卫星运行轨道方案的跟踪弧段评估问题。

根据卫星发射轨道参数和运行轨道参数，对上面轨道规划设计所需要解决的问题构建数学模型。

在理想情况下，卫星在发射和入轨过程受到的推力可近似为“两次冲量”的作用，这两次冲量分别位于发射点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-elsin(el)

发射轨道半长轴al计算如下：

卫星入轨时刻ti的计算如下：

其中μ＝398600.4415km³/s²为引力常数。

对于地心经纬度为的发射点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)、自转和极移(polarmotion)等地球复杂运动模型的影响。只有考虑了这些地球复杂运动模型因素，卫星轨道规划计算结果才可以取得和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)为极移(polarmotion)转换矩阵。

定义旋转矩阵算子如下:

定义旋转矩阵算子如下:

定义原子时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*t²+0.017998t³)/3600.0

θ＝(2004.3109t-0.42665t²-0.041833t³)/3600.0

z＝(2306.2181t+1.09468t²+0.018203t³)/3600.0

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

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

章动转换矩阵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.310t²+0.064t³)/3600.0

m＝357.52772333+(129596581.2240t-0.577t²-0.012t³)/3600.0

e＝93.27191028+(1739527263.1370t-13.257t²+0.011t³)/3600.0

f＝297.85036306+(1602961601.3280t-6.891t²+0.019t³)/3600.0

g＝125.04452222+(-6962890.5390t+7.455t²+0.008t³)/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值。

本发明实施例提供的大椭圆轨道设计针对有效倾角范围内的角度值以特定间隔(当前设置为0.01度)进行遍历，计算每个倾角角度值i对应卫星轨道星下点轨迹在多次过境目标点纬度圈上的经度值，并且计算目标点地心经度和星下点轨迹过境标点纬度圈时地心经度值之间的差值。选取地心经度差值满足指标要求的卫星轨道作为备选方案，最后以优选准则对备选轨道方案进行选取。此外在卫星轨道设计流程中考虑了太阳光照参数对于光学载荷卫星轨道设计的实施要求,并根据轨道参数计算在太阳高度角约束条件下的卫星发射窗口。

对于任意给定的发射点l，以发射点为中心，以经纬度线为坐标轴构建四个象限，如图12所示。在每个象限内选取一定的角度范围作为卫星轨道倾角的遍历范围。

图12中的四个象限分别对应卫星的四种发射入轨模式：1逆行下降，2逆行上升，3顺行下降，4顺行上升。大椭圆轨道设计的倾角遍历范围选取与小倾角轨道设计的倾角遍历范围选取一致，选取卫星逆行轨道倾角遍历范围[domaindegreelowerw，domaindegreeupperw]度，选取卫星顺行轨道倾角遍历范围[domaindegreelowere，domaindegreeuppere]度。

大椭圆轨道设计中，本发明对于4种发射入轨模式遍历范围内的每条卫星轨道，计算卫星入轨后6次星下点过境目标纬度圈时刻的经度值。

当前卫星运行轨道为椭圆轨道，选取卫星轨道入轨点为大椭圆轨道近地点，则入轨时刻的真近点角fi＝0。根据j2摄动条件下的卫星椭圆运行轨道方程，可以计算出t1、t2、t3、t4、t5和t6时刻卫星星下点过境目标纬度圈时在ecef地心地固坐标系下的大地经度lont1、lont2、lont3、lont4、lont5和lont6，其计算如下所示：

其中大椭圆轨道卫星的星下点轨迹计算oncomputevaluelongcellipse为椭圆轨道卫星星下点轨迹计算模块，m_launchreadytimejd为发射部署准备完成时刻的儒略日时间，ωe为地球自转角速度。为轨道升交点赤经导数，为近地点辐角导数，tω为轨道周期。

大椭圆(小倾角圆)轨道设计oncomputeellipseorbit模块对于当前4种发射入轨模式下的轨道倾角遍历范围内所有卫星轨道，调用遍历范围内卫星大椭圆轨道星下点过境目标纬度圈地心经度分布计算oncomputeorbitminvalueellipse模块计算卫星入轨后卫星星下点6次过境目标纬度圈时刻的大地经度分布，选取遍历范围内过境目标纬度圈时刻星下点经度与目标点经度偏差值最小的轨道进行门限判决，并保存有效轨道的计算参数到backuporbit数组。在卫星轨道设计流程中还需要考虑太阳光照参数对于光学载荷卫星轨道设计的实施要求,并根据轨道参数计算在太阳高度角约束条件下的卫星发射窗口。

小倾角圆轨道设计中，本发明对于4种发射入轨模式遍历范围内的每条卫星轨道，计算卫星入轨后6次星下点过境目标纬度圈时刻的经度值。定义卫星入轨后6次过境目标纬度圈的时刻分别为t1、t2、t3、t4、t5和t6，具体计算如下所示：

t1＝(ut1-ui)/n,t2＝(ut2-ui)/n,t3＝t1+tω,

t4＝t2+tω,t5＝t1+2tω,t6＝t2+2tω

其中为卫星平均角速度，μ＝398600.4415km³/s²为引力常数，轨道周期tω的计算如下所示：

其中地球扁率二阶项j2＝1082.6267*10^-6。

当前卫星运行轨道为圆形轨道，选取卫星轨道升交点与轨道近地点重合，则运行轨道的近地点辐角ωi＝0。根据j2摄动条件下的卫星圆形运行轨道方程，可以计算出t1、t2、t3、t4、t5和t6时刻卫星星下点过境目标纬度圈时在ecef地心地固坐标系下的大地经度lont1、lont2、lont3、lont4、lont5和lont6，其计算如下所示：

其中oncomputevaluelongccircle为圆形轨道卫星星下点轨迹计算模块，m_launchreadytimejd为发射部署准备完成时刻的儒略日时间，ωe为地球自转角速度，轨道升交点赤经导数和近地点辐角导数的计算如下所示：

其中re为地球平均半径，j2为地球扁率二阶项。oncomputesmallorbit模块对于当前4种发射入轨模式下的轨道倾角遍历范围内所有卫星轨道，调用oncomputeorbitminvaluesmall模块计算卫星入轨后卫星星下点6次过境目标纬度圈时刻的大地经度分布，选取遍历范围内过境目标纬度圈时刻星下点经度与目标点经度偏差值最小的轨道进行门限判决，并保存有效轨道的计算参数到backuporbit数组。在卫星轨道设计流程中还需要考虑太阳光照参数对于光学载荷卫星轨道设计的实施要求,并根据轨道参数计算在太阳高度角约束条件下的卫星发射窗口。

本发明提供的针对大椭圆轨道以及小倾角圆轨道的轨道规划方法，大椭圆轨道以及小倾角圆轨道设计oncomputeellipseorbit模块流程如图13所示，包括如下步骤：

步骤1、大椭圆轨道和小倾角圆轨道均具有四种入轨模式：即逆行下降、逆行上升、顺行下降、顺行上升共4种入轨模式；根据轨道半长轴的上下限，获得逆行上升和逆行下降的轨道倾角搜索范围、顺行上升和顺行下降的轨道倾角搜索范围。

如果卫星载荷类型为ccd光学设备，则以发射部署完成时刻为起点，计算目标点区域24小时内太阳高度角满足门限要求的utc时间观测窗口。

该轨道规划方法针对4种入轨模式对应轨道的倾角搜索范围内所有卫星轨道，计算卫星轨道星下点轨迹6次过境目标纬度圈时星下点的经度分布，并且选取满足侦察要求的最优轨道；设定初始状态为，mo为入轨模式，mo的取值为1、2、3和4，分别指代逆行下降、逆行上升、顺行下降、顺行上升共4种入轨模式，mo初始值设定为1；j为过境次数，j的初值设定为1，过境总次数为6；backuporbitnum为轨道方案个数，backuporbitnum的初值设定为0。

步骤2、判断j的值是否满足j<＝6，若是则进入步骤3，否则进入步骤11。

步骤3、判断mo的值是否满足mo<＝4；若是则进入步骤4，否则将mo重新设定为1，j自增1，返回步骤2。

步骤4、对于第mo种入轨模式轨道倾角遍历范围内所有卫星轨道，调用遍历范围内卫星大椭圆轨道星下点过境目标纬度圈地心经度分布计算oncomputeorbitminvalueellipse模块计算第j次过境目标纬度圈时星下点经度分布,选取轨道倾角遍历范围内过境时刻星下点经度与目标点经度偏差值最小的轨道。

步骤5、判断步骤4中所选取的偏差值最小的轨道对应的偏差值是否满足预设的指标要求，若是则进入步骤6，否则mo自增1，返回步骤3。

步骤6、判断当前卫星载荷类型是否为ccd光学设备，若是则进入步骤7，否则进入步骤10。

步骤7、判断当前轨道卫星过境目标点时刻是否在ccd光学设备的utc时间观测窗口内，若是则进入步骤10，否则进入步骤8。

步骤8、计算发射延时，对于当前入轨模式轨道倾角遍历范围内的所有卫星轨道，调用遍历范围内卫星大椭圆轨道星下点过境目标纬度圈地心经度分布计算oncomputeorbitminvalueellipse模块计算当前第j次过境目标纬度圈时星下点经度分布,选取遍历范围内过境时刻星下点经度与目标点经度偏差值最小的轨道。

步骤9、判断步骤8中所选取的偏差值最小的轨道是否满足预设的指标要求，若是则进入步骤10，否则返回步骤3；

步骤10、当前第backuporbitnum个轨道规划方案满足初步要求，将当前第backuporbitnum个轨道规划方案存入在backuporbit数组中,并且backuporbitnum自增1。

步骤11、大椭圆或者小倾角圆轨道规划完成，将当前对应的共backuporbitnum个轨道方案存入预先构建的backuporbit数组中。

本发明实施例针对遍历范围内卫星大椭圆轨道星下点过境目标纬度圈地心经度分布计算oncomputeorbitminvalueellipse模块如图14所示，具体包括如下步骤：

步骤101、开始在遍历范围内轨道星下点过境目标纬度圈地心经度分布；

步骤102、轨道倾角遍历范围的参数初始化：轨道倾角遍历初始值domaindegreelower,遍历间隔0.01度,遍历次数indexnum,索引index＝1；

步骤103、判断是否满足index<＝indexnum；若是则进入步骤104；否则进入步骤105；

步骤104、对当前轨道倾角为iiangledegree＝domaindegreelower+0.01×(index-1)的卫星发射轨道和卫星运行轨道，调用大椭圆轨道卫星的发射轨道和运行轨道计算oncomputeorbitvalueellipse模块进行大椭圆轨道星下点轨迹的计算，返回卫星星下点轨迹多次过境目标纬度圈时的星下点地心经度与目标点地心经度的偏差分布数据,并对索引值index进行累积自增1处理，返回步骤103；

步骤105、对于遍历范围内indexnum条卫星轨道星下点轨迹第j次过境目标纬度圈时的星下点地心经度与目标点地心经度的偏差分布进行统计，选取偏差值最小的轨道方案,并对该最小偏差值minvalue进行门限threshold判决；j取值为1～6。

步骤106、若minvalue<＝threshold，则进入步骤107；否则返回无效值；

步骤107、minvalue对应的轨道方案保存在mintimeorbit数组中，并返回有效轨道倾角值。

s104中、大椭圆轨道卫星的发射轨道和运行轨道计算oncomputeorbitvalueellipse模块如图15所示，具体包括如下步骤：

s10421、卫星发射轨道参数初始化计算部分：

1.发射点地心经纬度为发射轨道地心角βl和发射轨道倾角i为模块输入参数；

2.发射时刻卫星在发射轨道上真近点角flradian和入轨时刻卫星在发射轨道上真近点角firadian的计算；

3.发射轨道偏心率el、发射轨道半通径p和发射轨道半长轴al的计算

4.发射点偏近点角elangleradian和发射点平近点角mlangleradian的计算

5.根据上述参数计算卫星入轨时刻ti，选取卫星入轨点为近地点；

6.进行发射时刻gmst格林尼治平恒星时gl的计算

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

αl＝onlimitpiradian(gl+λl)

s10422、若为下降入轨模式，包括顺行下降和逆行下降，发射时刻卫星纬度辐角计算如下:

轨道升交点赤经:ω＝onlimit2piradian(αl-tan^-1(tan(ul)×cos(i))+π)；

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

轨道升交点赤经:ω＝onlimit2piradian(αl-tan^-1(tan(ul)×cos(i))+0)；

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

1.入轨时刻卫星纬度辐角ui＝ul+βl,目标点地心经纬度卫星平均角速度n,地球自转转速ωe；

2.入轨点地心纬度为

s10424、若为下降入轨模式，则分如下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种情况.目标点在南半球,目标点纬度小于入轨点纬度；

若即目标点地心纬度大于等于入轨点地心纬度，且入轨点在下降段则对于入轨后第1次过境，tempradian＝onlimit2piradian(ut1)；

若tempradian＜π属于第5种情况，δu＝tempradian；utemp＝π-2δu；否则属于第6种情况，则δu＝fabs(tempradian-2π)；utemp＝π+2δu；

若即目标点地心纬度大于等于入轨点地心纬度，且入轨点在上升段则对于入轨后第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；

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

①.计算卫星运行轨道周期tω、升交点赤经导数和近地点辐角导数和平近点角导数

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

当为大椭圆轨道时，采用如下方式计算卫星入轨后6次过境目标纬度圈的时刻：根据入轨后第1次进入目标纬度圈时刻的真近点角ff1＝ut1-ui先转为偏近点角ee1，再转为平近点角mm1；根据入轨后第2次进入目标纬度圈时刻的真近点角ff2＝ut2-ui，先转为偏近点角ee2，再转为平近点角mm2；

卫星入轨后6次过境目标纬度圈的时刻分别为t1～t6；

当为小倾角圆轨道时，采用如下方式计算卫星入轨后6次过境目标纬度圈的时刻：

t1＝ti+(ut1-ui)/n,t2＝ti+(ut2-ui)/n,

t3＝t1+tω,t4＝t2+tω,t5＝t1+2tω,t6＝t2+2tω

③.调用大椭圆轨道卫星的星下点轨迹计算oncomputevaluelongcellipse模块计算入轨后星下点第1次和第2次过境目标纬度圈时地心经度lont1和lont2；

④.计算入轨后星下点第3次、第4次、第5次和第6次过境目标纬度圈时刻的地心经度，分别为lont3、lont4、lont5以及lont6：

⑤.计算卫星星下点轨迹6次过境目标纬度圈时刻地心经度与目标地心经度的偏差；

s10426、将计算结果和中间变量存放到orbitvector数组。

大椭圆轨道卫星的星下点轨迹计算oncomputevaluelongcellipse模块，如图16所示，包括如下步骤：

sa1.卫星半长轴ai、偏心率ei和当前时刻t的真近点角ft,计算当前卫星在地心轨道坐标系下的位置矢量：

当为大椭圆轨道时，

当为小倾角圆轨道时：

r＝[aicos(ut)aisin(ut)0]′

sa2.根据相对入轨时刻时间差δt和升交点赤经导数计算当前时刻卫星运行轨道的升交点赤经其中ωi为卫星入轨时刻的升交点赤经；

sa3.根据相对入轨时刻时间差δt和近地点辐角导数计算当前时刻卫星运行轨道的近地点辐角其中ωi为卫星入轨时刻的近地点辐角；

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

sa5.对地心轨道坐标系下的位置矢量进行坐标转转获取地心惯性eci坐标系下的位置矢量peci

peci＝tm×r

sa6.查询地球运动参数eop文件,调用eci坐标系转ecef坐标系onconvertecitoecef模块将eci位置矢量peci转换为ecef位置矢量pecef；

sa7.根据ecef位置矢量pecef计算当前时刻圆形轨道卫星的星下点轨迹地心经度lon；

sa8.返回当前时刻大椭圆轨道卫星的星下点轨迹地心经度lon。

sa6中的eci坐标系转ecef坐标系onconvertecitoecef模块，其流程如图17所示，具体包括如下步骤：

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.30188×t²+0.017998t³)/3600.0

θ＝(2004.3109t-0.42665t²-0.041833t³)/3600.0

z＝(2306.2181t+1.09468t²+0.018203t³)/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；

sa607、返回ecef坐标下的坐标向量recef。

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

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