航天器可达域的一种确定方法及装置与流程

本发明涉及航天器空间任务规划与设计技术领域，尤其是涉及航天器轨道机动的可达域的一种确定方法及装置。

背景技术：

所谓航天器的可达域，简单来说，就是指在自身燃料限制条件下，航天器通过轨道机动能够到达的所有空间位置点的集合。有的文献也将可达域称为“可达区域”、“可达范围”等，这些名称在本质上是相同的。航天器的可达域既可作为航天器执行各种空间任务的初步可行性判据，又可为空间交会、在轨拦截、空间规避等空间任务的设计提供参考。对可达域进行研究，不仅可以得到航天器轨道机动可达范围的工程判据，还可能得出空间两点边值问题解的存在性的相关数学结论。此外，相关研究成果还可用于行星际轨道任务的研究。

目前，国内外对航天器可达域的研究还比较少，尚未形成一套成熟的理论或方法。现有技术主要是基于变轨脉冲大小固定、方向任意，或大小任意、方向固定，或方向沿轨道切向等简化假设条件来进行研究，适用范围存在较大的局限性，且计算方法主要是计算量巨大的数值遍历算法或精度较差的估算方法。因此，现有技术在实际工程应用中具有较大局限性，且由于计算量和计算精度的问题难以应用于航天器的在线任务规划与设计。

技术实现要素：

有鉴于此，本发明的目的在于提供航天器可达域的一种确定方法及装置，提高航天器可达域的实用性。

第一方面，本发明实施例提供了航天器可达域的一种确定方法，其包括以下步骤：

对于给定的转移轨道面，求取变轨脉冲的异面变轨能量的变化范围；

根据异面变轨能量的变化范围，求取可达域的多个特征椭圆，其中，多个特征椭圆至少包括两个第一边界特征椭圆和两个第二边界特征椭圆；

通过多个特征椭圆中的至少两个第一边界特征椭圆所围区域的并集的边界，逼近并作为可达域的外边界；

通过多个特征椭圆中的至少两个第二边界特征椭圆所围区域的交集的边界，逼近并作为可达域的内边界。

第二方面，本发明实施例提供了航天器可达域的一种确定装置，包括：

异面变轨能量变化范围求取模块，对于给定的转移轨道面，求取变轨脉冲的异面变轨能量的变化范围；

特征椭圆求取模块，用于根据异面变轨能量的变化范围，求取可达域的多个特征椭圆，其中，多个特征椭圆至少包括两个第一边界特征椭圆和两个第二边界特征椭圆；

可达域外边界求取模块，用于通过多个特征椭圆中的至少两个第一边界特征椭圆所围区域的并集的边界，逼近并作为可达域的外边界；

可达域内边界求取模块，用于通过多个特征椭圆中的至少两个第二边界特征椭圆所围区域的交集的边界，逼近并作为可达域的内边界。

第三方面，本发明实施例提供了一种服务器，其包括存储器和处理器，存储器中存储有可在所述处理器上运行的计算机程序，处理器执行所述计算机程序时实现上述第一方面所述方法的步骤。

第四方面，本发明实施例提供了一种具有处理器可执行的非易失的程序代码的计算机可读介质，所述程序代码使所述处理器执行第一方面所述的方法。

本发明实施例带来了以下有益效果：

本发明实施例提供的航天器可达域的一种确定方法，对于给定的转移轨道面，求取变轨脉冲的异面变轨能量的变化范围；根据异面变轨能量的变化范围，求取可达域的多个特征椭圆，其中，多个特征椭圆至少包括两个第一边界特征椭圆和两个第二边界特征椭圆；通过多个特征椭圆中的至少两个第一边界特征椭圆所围区域的并集的边界，逼近并作为可达域的外边界；通过多个特征椭圆中的至少两个第二边界特征椭圆所围区域的交集的边界，逼近并作为所述可达域的内边界，该方法计算量小，计算精度高，提高了可达域的实用性，且非常适合于航天器的在线任务规划与设计。

本发明的其他特征和优点将在随后的说明书中阐述，并且，部分地从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点在说明书、权利要求书以及附图中所特别指出的结构来实现和获得。

为使本发明的上述目的、特征和优点更加清晰易懂，下文特举较佳实施例，并配合所附附图，作详细说明。

附图说明

为了更清楚地说明本发明具体实施方式或现有技术中的技术方案，下面将对具体实施方式或现有技术描述中所需使用的附图作简要介绍。显而易见，下面描述中的附图是本发明的一些实施方式，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例提供的航天器各运动参数在轨道坐标系中的表示；

图2为本发明实施例1提供的航天器可达域的一种确定方法的流程图；

图3为本发明实施例1提供的集合v1确定的轨道簇的分布情况；

图4为本发明实施例1提供的异面变轨能量为给定值时的可达域；

图5为本发明实施例1提供的一值可达域随异面变轨能量变化的典型趋势示意图；

图6为本发明实施例1提供的可达域的内边界和外边界示意图；

图7为本发明实施例2提供的航天器可达域的一种确定装置的示意图；

图8为本发明实施例3提供的仿真实验1中，一系可达域的四个特征椭圆在一系可达域内的分布情况仿真结果；

图9为本发明实施例3提供的仿真实验1中，一系可达域的近似内边界和近似外边界仿真结果；

图10为本发明实施例3提供的仿真实验1中，二系可达域的四个特征椭圆在二系可达域内的分布情况仿真结果；

图11为本发明实施例3提供的仿真实验1中，二系可达域的近似内边界和近似外边界仿真结果；

图12为本发明实施例3提供的仿真实验1中，整个可达域的近似内边界和近似外边界仿真结果；

图13为本发明实施例3提供的仿真实验2中，一系可达域的四个特征椭圆在一系可达域内的分布情况仿真结果；

图14为本发明实施例3提供的仿真实验2中，一系可达域的近似内边界和近似外边界仿真结果；

图15为本发明实施例3提供的仿真实验2中，二系可达域的四个特征椭圆在二系可达域内的分布情况仿真结果；

图16为本发明实施例3提供的仿真实验2中，二系可达域的近似内边界和近似外边界仿真结果；

图17为本发明实施例3提供的仿真实验2中，整个可达域的近似内边界和近似外边界仿真结果。

图标：1-异面变轨能量变化范围求取模块；2-特征椭圆求取模块；3-可达域外边界求取模块；4-可达域内边界求取模块。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合附图对本发明的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

一般来说，航天器的可达域是一个非常复杂的不规则区域，难以用解析方法准确地描述这个区域。但是，在实际应用中，为了判断给定的转移轨道面上的一个目标点是否可达，需要能够快速计算出航天器的可达域。基于此，本发明实施例提供的航天器可达域的一种确定方法及装置计算结果较为准确，且计算量很小，便于工程应用。

如图1所示，定义轨道坐标系ogxgygzg，其坐标原点og位于地心，ogxg轴沿当前时刻航天器地心矢径方向，ogzg轴沿当前时刻航天器动量矩矢量方向，ogyg轴由右手定则确定。设在轨道坐标系ogxgygzg中，航天器初始轨道的动量矩矢量hs＝[0,0,hs]^t，航天器所能施加的最大变轨脉冲为δvmax，在航天器变轨点p1处，航天器的位矢r1＝[r1,0,0]^t，速度v1＝[v1xg,v1yg,0]^t，航天器施加的变轨脉冲δv＝[δvxg,δvyg,δvzg]^t，施加变轨脉冲后的瞬间速度为vk，航天器转移轨道的动量矩矢量为hk＝[hkxg,hkyg,hkzg]^t。

设轨道坐标系ogxgygzg三个坐标轴方向的单位矢量分别为ig、jg和kg，则转移轨道的动量矩矢量：

hk＝r1×vk＝r1×(v1+δv)＝-r1δvzgjg+(hs+r1δvyg)kg1.1

可见，hk是由变轨脉冲δv的δvyg和δvzg两个分量决定的，而与δvxg无关，且hk位于ogygzg平面内，转移轨道面与ogygzg平面垂直。

航天器转移轨道面的空间方位可以用转移轨道面相对于初始轨道面的倾角即相对轨道倾角ir来表示。根据式1.1，相对轨道倾角ir由下式定义

可见，ir是由δvyg和δvzg决定的，与δvxg无关。由于航天器的变轨能力有限，从而相对轨道倾角的大小|ir|一般为一个较小的锐角，即ir∈(-π/2,π/2)。显然，ir与航天器的转移轨道面之间是一一对应的。

设航天器初始轨道远地点处的地心距为ra，速度大小为va，一般情况下，有va＞δv，从而hs/r1≥hs/ra＝va＞δv≥|δvyg|。于是，有hs/r1+δvyg＞0。显然，若δvzg＞0，即航天器向ogzg轴正向变轨，则ir∈(-π/2,0)；若δvzg＝0，即航天器进行同轨道面变轨，则ir＝0；若δvzg＜0，即航天器向ogzg轴负向变轨，则ir∈(0,π/2)。

由于航天器向ogzg轴正向和负向变轨的能力是相同的，从而下面假设δvzg＞0，以确定此时相对轨道倾角ir的变化范围。

令则有δvyg＝δvyzcosλ，δvzg＝δvyzsinλ，其中，λ∈(0,π)。代入式1.2，有

其中，

当时，f(λ)取得最小值：

式中，irm是在给定的δvyz条件下相对轨道倾角ir的最小值，由式1.4可得，可见，当δvyz＝δvmax时，irm取得最小值：

irmin即为航天器助推器所能施加的最大变轨脉冲为δvmax时相对轨道倾角ir的最小值。将ir取最小值irmin时的转移轨道面称为负极轨道面。

又由δvzg＞0，可得ir＜0。从而，当δvzg＞0时，ir∈[irmin,0)。

由于航天器向ogzg轴正向和负向变轨的能力是相同的，从而当δvzg≤0时，ir的最大值为：

将ir取最大值irmax时的转移轨道面称为正极轨道面。又由δvzg≤0，可得ir≥0。从而，当δvzg≤0时，ir∈[0,irmax]。

因此，当航天器助推器所能施加的最大变轨脉冲为δvmax时，航天器转移轨道面的相对轨道倾角ir的变化范围为[-irmax,irmax]。显然，航天器的所有转移轨道面都处于正极轨道面和负极轨道面之间。

由式1.1可见，只要δvyg、δvzg确定，转移轨道的动量矩矢量hk便确定了。因此，具有不同δvxg，相同δvyg、δvzg的变轨脉冲δv确定的转移轨道具有相同的动量矩矢量，这些转移轨道构成了位于同一个轨道面内的一个轨道簇，且这个轨道簇经过航天器的变轨点p1。

另外，由于相对轨道倾角ir是由δvyg和δvzg决定的，与δvxg无关，从而最大相对轨道倾角本质上是由变轨时分配给ogyg轴和ogzg轴方向上的变轨能量之和决定的。由式1.4可知，在δvyz确定的条件下，ir的最大值为因此，就轨道倾角的改变量而言，δvyz决定了航天器异面变轨的能力，从而将其称为异面变轨能量。

在本发明中，假设地球为一个质点，忽略各种摄动因素对航天器运动的影响。本发明所讨论的航天器的可达域具体定义如下：

定义：在航天器助推器所能施加的最大变轨脉冲为δvmax的条件下，给定转移轨道面即给定相对轨道倾角ir，航天器从当前位置点p1处变轨所能到达的该转移轨道面上的所有位置点的集合称为航天器在给定转移轨道面上的可达域，简称为可达域。

为便于对本实施例进行理解，首先对本发明实施例所公开的航天器可达域确定方法进行详细介绍。

实施例1：

如图2所示，本发明实施例提供了航天器可达域的一种确定方法，其包括以下步骤：

s1.对于给定的转移轨道面，求取变轨脉冲的异面变轨能量的变化范围。

s2.根据异面变轨能量的变化范围，求取可达域的多个特征椭圆，其中，多个特征椭圆至少包括两个第一边界特征椭圆和两个第二边界特征椭圆。

s3.通过多个特征椭圆中的至少两个第一边界特征椭圆所围区域的并集的边界，逼近并作为可达域的外边界。

s4.通过多个特征椭圆中的至少两个第二边界特征椭圆所围区域的交集的边界，逼近并作为可达域的内边界。

需要说明的是，本实施例不对可达域的外边界和内边界的求取顺序进行限制，同时也不对各个特征椭圆的求取顺序进行限制。

步骤s1中，给定航天器的转移轨道面，即给定相对轨道倾角ir∈[-irmax,irmax]，令η＝tanir，则有

设变轨时分配给ogyg轴和ogzg轴方向上的变轨能量之和为δvyz，则将式1.7代入可得

要使该方程有解，则要求

即

又由于δvyz≤δvmax，从而δvyz∈[δvyzmin,δvmax]，这就是异面变轨能量δvyz的变化范围。

本实施例中步骤s2包括以下步骤：

s21.根据异面变轨能量的变化范围，将可达域划分为一系可达域与二系可达域两个分支；

s22.求取一系可达域的四个特征椭圆，具体包括：

求取异面变轨能量在区间[δvyzmin,δvmax]内变化时，半长轴最大的一值可达域正极椭圆和半长轴最大的一值可达域负极椭圆，并将其作为两个第一边界特征椭圆；

求取异面变轨能量取最小值时对应的一值可达域的正极椭圆与负极椭圆，并将其作为两个第三边界特征椭圆；

s23.求取二系可达域的两个特征椭圆，即异面变轨能量在区间[δvyzmin,δvmax]内变化时，半长轴最小的二值可达域正极椭圆和半长轴最小的二值可达域负极椭圆，并将其作为两个第二边界特征椭圆。

对于步骤s21，当异面变轨能量固定时，即给定δvyz∈[δvyzmin,δvmax]，方程式1.8有两个实根，这两个实根及相应的δvzg值分别如下：

显然，在异面变轨能量δvyz可变时，上述两组解都是δvyz的函数。

定义变轨脉冲δv的两个集合：

则集合v1和v2在同一个轨道面内确定了两个轨道簇。

以集合v1确定的轨道簇为例，研究该轨道簇中各条转移轨道之间的位置关系，从而确定该轨道簇的分布区域。

该轨道簇中所有转移轨道的动量矩矢量：

hk1＝-r1δvzg1jg+(hs+r1δvyg1)kg1.13

对于转移轨道的偏心率矢量：

显然，

与δvxg无关，为一个常量，其中，为转移轨道的半通径，也与δvxg无关。ek1yg、ek1zg与δvxg有关。

偏心率矢量的大小

可见，若δvxg＜-v1xg，则ek1随着δvxg单调递减；若δvxg≥-v1xg，则ek1随着δvxg单调递增。

设航天器的变轨点p1在转移轨道上的真近点角为θk1，则由

由此可知：

1.若v1xg+δvxg＞0，则ek1×r1与hk1同向，这等价于θk1∈(0,π)。

2.若v1xg+δvxg＝0，则ek1×r1＝0，这等价于θk1＝0或π(具体值由ek1确定)。

3.若v1xg+δvxg＜0，则ek1×r1与hk1反向，这等价于θk1∈(π,2π)。

设航天器从变轨点p1处经过转移角θ(范围为(0,2π))后到达点p2，点p2处地心矢径的大小为r2，则由r1＝pk1/(1+ek1cosθk1)可得

从而，

1.若0＜θ＜π，则

(1)当δvxg＜-v1xg时，θk1∈(π,2π)。此时，ek1随δvxg单调递减。由于ek1cosθk1＝ek1xg为一个定值，且ek1sinθk1＜0，从而ek1sinθk1随δvxg单调递增。由于sinθ＞0，从而ek1sinθk1sinθ随δvxg单调递增。因此，r2随δvxg单调递增。

(2)当δvxg＝-v1xg时，θk1＝0或π。此时，

(3)当δvxg＞-v1xg时，θk1∈(0,π)。此时ek1随δvxg单调递增。由于ek1cosθk1＝ek1xg为一个定值，且ek1sinθk1＞0，从而ek1sinθk1随δvxg单调递增。由于sinθ＞0，从而ek1sinθk1sinθ随δvxg单调递增。因此，r2随δvxg单调递增。

由上述(1)、(2)和(3)可知，当0＜θ＜π时，r2随δvxg单调递增。

2.若θ＝π，则

可见，r2与δvxg无关，为一个常量。也就是说，轨道簇中的所有转移轨道都经过ogxg轴负半轴上的同一点[-rc1,0,0]^t，记为点pc1。

3.若π＜θ＜2π，则通过与1中类似的讨论可知，r2随δvxg单调递减。

综合1、2和3可得到如下结论：

(a)轨道簇中的所有转移轨道有两个共同的交点p1和pc1，这两个交点分别位于ogxg轴的正半轴和负半轴上，而且轨道簇中的任意两个转移轨道有且只有这两个交点。

(b)当转移角θ∈(0,π)时，对于相同的转移角θ，转移终点p2的地心距r2随着δvxg的增大逐渐增大；当转移角为θ∈(π,2π)时，对于相同的转移角θ，转移终点p2的地心距r2随着δvxg的增大逐渐减小。

根据以上结论，可以得到集合v1确定的轨道簇的分布情况如图3所示。

在集合v1中，将δvxg取最小值时的变轨脉冲确定的转移轨道称为负极椭圆，将δvxg取最大值时的变轨脉冲确定的转移轨道称为正极椭圆。显然，该轨道簇分布区域的边界完全是由正极椭圆和负极椭圆组成的。在给定的转移轨道面上，设正极椭圆所围区域为a⁺，负极椭圆所围区域为a^-，则该轨道簇的分布区域a1可以表示为该区域称为一值可达域。

类似地，集合v2确定的轨道簇的分布情况与集合v1确定的轨道簇的分布情况完全类似。只不过，根据式1.7可得：

由于hs+r1δvygi＞0(i＝1,2)，且δvyg1＞δvyg2，从而，相应的两个轨道簇的动量矩矢量的大小hk1＞hk2，两个轨道簇的半通径pk1＞pk2，两个轨道簇与ogxg轴负半轴交点的地心距rc1＞rc2。因此，两个轨道簇与轨道坐标系的ogxg轴的负半轴的交点pc1和pc2不重合。

将集合v2确定的轨道簇的分布区域称为二值可达域，记为a2。在δvyz确定的条件下，航天器的可达域a＝a1∪a2，如图4中的整个阴影区域所示。

以上在异面变轨能量δvyz确定的条件下研究了航天器的一值可达域和二值可达域。如前所述，对于一个确定的转移轨道面(即ir确定)，异面变轨能量δvyz是可变的，且δvyz∈[δvyzmin,δvmax]。当δvyz变化时，一值可达域和二值可达域也将随之变化，且这两个区域的变化趋势是不同的。当δvyz取集合[δvyzmin,δvmax]中所有的确定值时，所有相应的一值可达域的并集称为一系可达域，所有相应的二值可达域的并集称为二系可达域。一系可达域和二系可达域为整个可达域的两个分支，两者的并集即为整个可达域。

步骤s22中，两个第一边界特征椭圆和两个第三边界特征椭圆可以通过分析一系可达域的分布情况获得，而一系可达域的分布可以通过分析异面变轨能量δvyz变化时相应的一值可达域的变化规律获得。由于一值可达域是由正极椭圆和负极椭圆围成的，从而要获取δvyz变化时一值可达域的变化规律，只需得到δvyz变化时一值可达域的正极椭圆和负极椭圆的变化趋势即可。根据其变化趋势，选取合适的正极椭圆和负极椭圆，作为两个第一边界特征椭圆，以逼近一系可达域的外边界，同时，选取合适的正极椭圆和负极椭圆，作为两个第三边界特征椭圆，以逼近一系可达域的内边界。

工程实际中，由于航天器所能施加的最大变轨脉冲δvmax比较有限，远小于航天器的速度大小，从而对于任意变轨脉冲δv，转移轨道的偏心率矢量相对于初始轨道的偏心率矢量变化都比较小。同时，δvyz的变化区间[δvyzmin,δvmax]相对来说也是一个比较小的区间，在此区间内，相应一值可达域的正极椭圆和负极椭圆的偏心率矢量变化都很小，也即这两个椭圆的形状和长轴方向变化很小。一值可达域随δvyz变化的典型趋势如图5所示。图5中，i、ii、iii和iv四个区域分别为δvyz＝δv1、δv2、δv3和δv4(δv1＜δv2＜δv3＜δv4)时的一值可达域。可见，随着δvyz的增大，一值可达域的正极椭圆和负极椭圆的形状和长轴方向变化不大，而半长轴却显著变化。因此，在研究正极椭圆和负极椭圆随δvyz的变化趋势时，可以忽略这两个椭圆偏心率矢量的变化，而只需研究这两个椭圆的半长轴随δvyz的变化规律即可。

正极椭圆半长轴随异面变轨能量的变化规律

对于一值可达域，正极椭圆对应的变轨脉冲为其中δvyg1、δvzg1的表达式如式1.11所示。在正极椭圆上，p1点处的速度为：

由活力公式可知，正极椭圆的半长轴

显然，akp随单调递增，且由式1.22可得

将δvyg1对δvyz求导可得

从而，有

由轨道坐标系的定义可知，v1yg＞0。从而，

1.若v1xg≤0则即随δvyz单调递增，从而akp随δvyz单调递增。

2.若v1xg＞0，则令可得

由易知，δvc∈[δvyzmin,δvmax]。从而有：

(1)当δvyz≤δvc时，即随δvyz单调递增，从而akp随δvyz单调递增。

(2)当δvyz＞δvc时，即随δvyz单调递减，从而akp随δvyz单调递减。

因此，当δvyz＝δvc时，正极椭圆的半长轴akp取得最大值，而akp的最小值在区间[δvyzmin,δvmax]的一个端点处取得。

根据式1.24可知：

从而，若则当δvyz＝δvyzmin时，正极椭圆半长轴akp取得最小值；若则当δvyz＝δvmax时，正极椭圆半长轴akp取得最小值。

负极椭圆半长轴随异面变轨能量的变化规律

对于一值可达域，负极椭圆对应的变轨脉冲为其中δvyg1、δvzg1的表达式如式1.11所示。在负极椭圆上，p1点处的速度为：

从而

类似地，可知

1.若v1xg≥0，则即随δvyz单调递增，从而akn随δvyz单调递增。

2.若v1xg＜0，则有：

(1)当δvyz≤δvc时，即随δvyz单调递增，从而akn随δvyz单调递增。

(2)当δvyz＞δvc时，即随δvyz单调递减，从而akn随δvyz单调递减。

因此，当δvyz＝δvc时，负极椭圆的半长轴akn取得最大值而akn的最小值在区间[δvyzmin,δvmax]的一个端点处取得。

根据式1.30可知，

从而，若则当δvyz＝δvyzmin时，负极椭圆半长轴akn取得最小值；若则当δvyz＝δvmax时，负极椭圆半长轴akn取得最小值。

一系可达域的特征椭圆的选取

由一系可达域的定义可知，一系可达域实际上是异面变轨能量δvyz从其最小值δvyzmin逐渐变化到其最大值δvmax时一值可达域所扫过的区域。

由式1.20可知，当δvyz在区间[δvyzmin,δvmax]内变化时，相应的一值可达域与轨道坐标系ogxg轴负半轴的交点pc1的地心距rc1随一值可达域中轨道的半通径pk1单调递增。由式1.21可知，一值可达域中轨道的动量矩矢量的大小hk1随δvyg1单调递增，从而pk1随δvyg1单调递增。又由式1.25可知，δvyg1随δvyz单调递增。从而rc1随δvyz单调递增。因此，当δvyz＝δvyzmin时，rc1取得最小值；当δvyz＝δvmax时，rc1取得最大值。

可见，一系可达域实际上是一个空心区域，这是因为当异面变轨能量δvyz在区间[δvyzmin,δvmax]内变化时，pc1的地心距rc1存在最小值(相应的交点记为pc1min)，在变轨点p1与pc1min之间有一个区域是不可达的，如图5所示。

根据以上分析，选取δvyz＝δvyzmin时对应的一值可达域的正极椭圆与负极椭圆作为两个特征椭圆，这两个特征椭圆称为第三边界特征椭圆，用于逼近一系可达域的内边界，分别记为和这两个特征椭圆各自所围区域的交集的边界可以作为一系可达域的近似内边界，记为li_inner。这是因为当δvyz＝δvyzmin时，相应的一值可达域与轨道坐标系ogxg轴负半轴的交点pc1的地心距rc1取得最小值。当δvyz从δvyzmin开始逐渐增大时，只要δvyz≤δvc，相应的一值可达域的正极椭圆和负极椭圆就在逐渐增大，从而相应的一值可达域不会扫过边界li_inner内的区域。当δvyz＞δvc时，一值可达域的正极椭圆(当v1xg＞0时)或者负极椭圆(当v1xg＜0时)开始减小，相应的一值可达域即便扫过边界li_inner内的区域，所扫过的这部分区域也很小。因此，所选取的这两个特征椭圆可以有效逼近一系可达域的内边界。

对于一系可达域的外边界，由于在异面变轨能量δvyz变化时一值可达域的正极椭圆和负极椭圆的偏心率矢量变化都很小，从而，半长轴最大的正极椭圆所围区域基本上可以覆盖δvyz变化时所有的正极椭圆，半长轴最大的负极椭圆所围区域也基本上可以覆盖δvyz变化时所有的负极椭圆，因此可以选取半长轴最大的正极椭圆和半长轴最大的负极椭圆作为两个特征椭圆，这两个特征椭圆称为第一边界特征椭圆，用于逼近一系可达域的外边界，分别记为和这两个特征椭圆各自所围区域的并集的边界可以作为一系可达域的近似外边界，记为li_outer。

综上所述，一系可达域共有四个特征椭圆，即两个第一边界特征椭圆和两个第三边界特征椭圆其中：和对应的变轨脉冲分别为：

根据一值可达域的正极椭圆和负极椭圆的半长轴随异面变轨能量的变化规律可知，和对应的变轨脉冲分别为：

显然，利用变轨脉冲可以直接求出相应的特征椭圆。

对于步骤s23，二系可达域的分析方法与一系可达域的分析方法类似，在此不予赘述，直接给出相关结论。

与一系可达域类似，二系可达域也是一个空心区域，且有四个特征椭圆。选取异面变轨能量δvyz在区间[δvyzmin,δvmax]内变化时，半长轴最小的二值可达域正极椭圆和半长轴最小的二值可达域负极椭圆作为两个特征椭圆，这两个特征椭圆称为第二边界特征椭圆，用于逼近二系可达域的内边界，分别记为和这两个特征椭圆各自所围区域的交集的边界可以作为二系可达域的近似内边界，记为lii_inner。两个特征椭圆和对应的变轨脉冲分别为：

选取δvyz＝δvyzmin时对应的二值可达域的正极椭圆和负极椭圆作为两个特征椭圆，并将其分别记为和这两个特征椭圆各自所围区域的并集的边界可以作为二系可达域的近似外边界，记为lii_outer。两个特征椭圆和对应的变轨脉冲分别为：

但是，注意到，当δvyz＝δvyzmin时，方程式1.8有两个相等实根，即δvyg1(δvyzmin)＝δvyg2(δvyzmin)和δvzg1(δvyzmin)＝δvzg2(δvyzmin)。从而，根据式1.33、1.34和式1.39、1.40可知，

这说明，特征椭圆和分别与两个第三边界特征椭圆和重合，从而和无需再进行计算。

因此，对于二系可达域，只需计算和这两个第二边界特征椭圆即可。

对于步骤s3和s4，由式1.41可见，当δvyz＝δvyzmin时，相应的一值可达域与二值可达域是重合的，这两者与ogxg轴的交点pc1和pc2也是重合的，即rc1＝rc2。当δvyz由最小值δvyzmin逐渐增大到最大值δvmax时，点pc1的地心距rc1单调增加，而点pc2的地心距rc2单调递减，从而，在ogxg轴负方向上，一值可达域向外扩张，形成一系可达域，二值可达域向内收缩，形成二系可达域。因此，一系可达域和二系可达域全向衔接，形成一个连续的空心区域，即整个可达域。

由于一系可达域和二系可达域的特征椭圆逼近一系可达域和二系可达域的内、外边界，从而可以利用这些特征椭圆逼近整个可达域的内边界和外边界，实现可达域的快速计算。由式1.41可知，一系可达域和二系可达域共有八个特征椭圆，而其中共有六个独立的特征椭圆：两个第一边界特征椭圆和两个第三边界特征椭圆和以及两个第二边界特征椭圆和

实际使用时，优先但不限于采用下列方法之一计算可达域的外边界：

方法1：通过两个第一边界特征椭圆所围区域的并集的边界，逼近并作为可达域的外边界。

方法2：通过两个第一边界特征椭圆以及两个第二边界特征椭圆所围区域的并集的边界，逼近并作为可达域的外边界。

方法3：通过两个第一边界特征椭圆以及两个第三边界特征椭圆所围区域的并集的边界，逼近并作为可达域的外边界。

方法4：通过两个第一边界特征椭圆、两个第二边界特征椭圆，以及两个第三边界特征椭圆所围区域的并集的边界，逼近并作为可达域的外边界。

优先但不限于采用下列方法之一计算可达域的内边界：

方法1：通过两个第二边界特征椭圆所围区域的交集的边界，逼近并作为可达域的内边界。

方法2：通过两个第一边界特征椭圆和两个第二边界特征椭圆所围区域的交集的边界，逼近并作为可达域的内边界。

方法3：通过两个第二边界特征椭圆和两个第三边界特征椭圆所围区域的交集的边界，逼近并作为可达域的内边界。

方法4：通过两个第一边界特征椭圆、两个第二边界特征椭圆，以及两个第三边界特征椭圆所围区域的交集的边界，逼近并作为可达域的内边界。

由于这六个特征椭圆都位于可达域内，且其中用于逼近一系可达域或二系可达域内边界的特征椭圆在有些地方也非常贴近相应可达域分支的外边界，而其中用于逼近一系可达域或二系可达域外边界的特征椭圆在有些地方也非常贴近相应可达域分支的内边界，从而利用计算可达域外边界的第4种方法和计算可达域内边界的第4种方法计算可达域，可以得到较高精度的计算结果。

如图6所示，在转移轨道面内，从地心og引射线l，射线l相对于ogxg轴的转移角为α∈[0,2π)，且该射线与可达域的内边界和外边界分别相交于点qmin和qmax，这两点的地心距离都是转移角α的函数，分别记为rmin(α)和rmax(α)。

设射线l分别与特征椭圆和相交于点q1、q2、q3、q4、q5和q6，且点qi(i＝1,2,...,6)的地心距离为ri(α)，则可达域的内、外边界点qmin和qmax的地心距离rmin(α)和rmax(α)可以用式1.42近似计算：

利用式1.42即可快速获得较高精度的可达域的内边界和外边界，即得到整个可达域。

实施例2：

如图7所示，本实施例提供了航天器轨道机动的可达域的一种确定装置，包括：异面变轨能量变化范围求取模块1、特征椭圆求取模块2、可达域外边界求取模块3和可达域内边界求取模块4，异面变轨能量变化范围求取模块1用于对给定的转移轨道面，求取变轨脉冲的异面变轨能量的变化范围；特征椭圆求取模块2用于根据异面变轨能量的变化范围，求取可达域的多个特征椭圆，其中，多个特征椭圆至少包括两个第一边界特征椭圆和两个第二边界特征椭圆；可达域外边界求取模块3用于通过多个特征椭圆中的至少两个第一边界特征椭圆所围区域的并集的边界，逼近并作为可达域的外边界；可达域内边界求取模块4用于通过多个特征椭圆中的至少两个第二边界特征椭圆所围区域的交集的边界，逼近并作为可达域的内边界。

本发明实施例提供的可达域确定装置与实施例1提供的可达域确定方法具有相同的技术特征，所以也能解决相同的技术问题，达到相同的技术效果。

其中，本实施例中的特征椭圆求取模块，包括：

可达域划分单元，用于根据异面变轨能量的变化范围，将可达域划分为一系可达域与二系可达域，且二系可达域近似位于一系可达域的内部；

两个第一边界特征椭圆求取单元，用于求取异面变轨能量变化过程中半长轴最大的一值可达域正极椭圆和半长轴最大的一值可达域负极椭圆，并将其作为两个第一边界特征椭圆；

两个第二边界特征椭圆求取单元，用于求取异面变轨能量变化过程中半长轴最小的二值可达域正极椭圆和半长轴最小的二值可达域负极椭圆，并将其作为两个第二边界特征椭圆。

进一步，特征椭圆求取模块，还包括：

两个第三边界特征椭圆求取单元，用于求取异面变轨能量取最小值时对应的一值可达域的正极椭圆与负极椭圆，并将其作为两个第三边界特征椭圆。

本发明实施例还提供了一种服务器，其包括存储器和处理器，存储器中存储有可在所述处理器上运行的计算机程序，处理器执行所述计算机程序时实现上述实施例1所述方法的步骤。

根据前述实施例，本发明实施例提供了一种具有处理器可执行的非易失的程序代码的计算机可读介质，所述程序代码使前述实施例所述的处理器执行实施例1所述的方法。

实施例3：

下面通过两组仿真实验来验证实施例1所述的航天器可达域的一种确定方法的准确性和快速性。仿真实验采用计算可达域外边界的第4种方法和计算可达域内边界的第4种方法，即式1.42计算可达域，所采用的计算平台为一个处理器为intel(r)core(tm)i7-3770cpu，安装内存为4.00gb的个人台式计算机。在仿真中，设在轨道坐标系中，航天器在变轨点p1处的位矢r1＝[9250000,0,0]^t(m)，速度v1＝[3500,5900,0]^t(m/s)，航天器所能施加的最大变轨脉冲δvmax＝600m/s。从而，由式1.6可知，航天器转移轨道面的相对轨道倾角ir的变化范围为[-5.838,5.838](°)。

实验1：共面变轨时的可达域计算

在本次实验中，取ir＝0°，即航天器进行共面变轨。由式1.10和式1.27可得，δvyzmin＝360.399m/s，δvc＝547.663m/s，从而δvyz∈[360.399,600](m/s)。由式1.33-1.38可得，六个独立的特征椭圆和对应的变轨脉冲分别如下：

一系可达域的四个特征椭圆和在一系可达域内的分布情况如图8所示。图8中，两个实线椭圆分别为和两个虚线椭圆分别为和灰色区域为实际的一系可达域。可见，和非常贴近一系可达域的外边界，和非常贴近一系可达域的内边界，同时也在部分地方非常贴近一系可达域的外边界。利用这四个特征椭圆得到的一系可达域的近似内边界li_inner和近似外边界li_outer如图9中的两条闭合实线所示。可见，这两个近似边界是比较准确的，尤其是近似内边界。因此，选取这四个椭圆作为特征椭圆以逼近一系可达域的内、外边界是合理有效的。

类似地，二系可达域的四个特征椭圆和(即)在二系可达域内的分布情况如图10所示。图10中，两个实线椭圆分别为和两个虚线椭圆分别为和灰色区域为实际的二系可达域。可见，和非常贴近二系可达域的外边界，同时也在部分地方贴近二系可达域的内边界，和非常贴近二系可达域的内边界。二系可达域的近似内边界lii_inner和近似外边界lii_outer如图11中的两条闭合实线所示。可见，这两个近似边界是比较准确的，尤其是近似外边界。因此，选取这四个椭圆作为二系可达域的特征椭圆是非常有效的。

根据上述六个独立的特征椭圆，利用所采用的可达域计算方法得到的可达域近似内边界和近似外边界如图12中的两条闭合实线所示，图12中，灰色区域为实际的整个可达域。由图12可见，对于共面变轨的情形，所采用的可达域计算方法是非常准确的。

经统计，所采用的可达域计算方法整个计算过程耗时0.015s。

实验2：异面变轨时的可达域计算

在本次实验中，取ir＝4.669°。此时，航天器实施异面变轨。由式1.10和式1.27可得，δvyzmin＝480.299m/s，δvc＝571.116m/s，从而，δvyz∈[480.299,600](m/s)。

为便于描述，在转移轨道面内建立坐标系ogx1y1z1，其坐标原点位于地心og，ogx1轴与ogxg轴重合，ogz1轴沿转移轨道的动量矩矢量hk方向，ogy1轴由右手定则确定。显然，航天器的可达域位于ogx1y1平面内。由式1.33-1.38可得，六个独立的特征椭圆和对应的变轨脉冲分别如下：

类似地，一系可达域的四个特征椭圆和在一系可达域内的分布情况如图13所示。利用这四个特征椭圆得到的一系可达域的近似内边界li_inner和近似外边界li_outer如图14中的两条闭合实线所示，可见，这两个近似边界是比较准确的，这也说明，选取这四个椭圆作为一系可达域的特征椭圆是非常有效的。

二系可达域的四个特征椭圆和(即)在二系可达域内的分布情况如图15所示。二系可达域的近似内边界lii_inner和近似外边界lii_outer如图16中的两条闭合实线所示，可见，这两个近似边界是比较准确的，这也说明，选取这四个椭圆作为二系可达域的特征椭圆是非常有效的。

利用所采用的可达域计算方法得到的可达域近似内边界和近似外边界如图17所示。由图17可见，对于异面变轨的情形，所采用的可达域计算方法也是非常准确的。

经统计，所采用的可达域计算方法整个计算过程耗时0.016s。

综合以上两个仿真实验的结果可见，所采用的可达域计算方法对于共面变轨和异面变轨的情形都是非常准确的，且由于计算过程只涉及到六个特征椭圆的求取，从而计算量非常小，计算时间在十几毫秒量级。因此，该方法非常适合于实际工程应用，尤其是在线应用的场合。

本发明实施例所提供的航天器可达域的一种确定方法及装置的计算机程序产品，包括存储了程序代码的计算机可读存储介质，所述程序代码包括的指令可用于执行前面方法实施例中所述的方法，具体实现可参见方法实施例，在此不再赘述。

所属领域的技术人员可以清楚地了解到，为描述的方便和简洁，上述描述的系统和装置的具体工作过程，可以参考前述方法实施例中的对应过程，在此不再赘述。

所述功能如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：u盘、移动硬盘、只读存储器(rom，read-onlymemory)、随机存取存储器(ram，randomaccessmemory)、磁碟或者光盘等各种可以存储程序代码的介质。

术语“第一”、“第二”、“第三”仅用于描述目的，而不能理解为指示或暗示相对重要性。

最后应说明的是：以上所述实施例，仅为本发明的具体实施方式，用以说明本发明的技术方案，而非对其限制，本发明的保护范围并不局限于此，尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，其依然可以对前述实施例所记载的技术方案进行修改或可轻易想到变化，或者对其中部分技术特征进行等同替换；而这些修改、变化或者替换，并不使相应技术方案的本质脱离本发明实施例技术方案的精神和范围，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。