连续推力机动轨道优化设计的贝塞尔曲线法

姚 玮，罗建军，谢剑锋，袁建平

(1.西北工业大学航天学院，西安 710072；2. 航天飞行动力学技术重点实验室，西安 710072；3. 北京航天飞行控制中心，北京 100094)

0 引 言

随着电推进技术的应用以及全电航天器和太阳帆航天器技术的发展，连续推力轨道在航天任务中的应用越来越广泛[1-8]。连续推力轨道机动的核心问题是轨道设计。对于连续推力轨道设计，目前已有的正方法难以满足轨道状态约束，而反方法由于转移轨道方程局限性太强，结果难以令人满意。反方法又叫形状法，学者们针对形状法做了大量研究。Petropoulos和Longuski[9]提出了一种利用指数正弦曲线的形状方法来处理平面内连续推力重力辅助轨迹。该方法只能满足位置一致，而无法实现速度一致。Izzo[10]利用基于指数正弦曲线的形状方法解决多圈兰伯特问题，以地火转移问题为例进行了研究。对于连续推力轨道拦截问题，Petropoulos[11]研究了很多其他的可行方法，包括对数螺旋线，笛卡尔卵形线和卡西尼卵形线。崔平远等[12]提出了一种深空探测任务发射机会的快速搜索算法。尚海滨等[13]针对行星际连续推力Lambert轨道问题，提出了一种用线性化策略简化问题的标称轨道设计方法。Wall和Conway[14]采用5阶和6阶逆多项式形状方法(IP)为各种航天器轨迹问题提供近似最优解。然而，逆多项式方法的最大缺陷是6阶多项式，最多只有5个极值。因此，当转移轨道周期大于3时，逆多项式法不能对椭圆轨道机动问题中转移轨道远地点和近地点的振荡特性进行建模。尚海滨等对Wall的6阶多项式逼近方法开展进一步研究，证明其结果可以为轨道优化过程提供有效的优化初值[15]，并把阶数提高到了15阶和20阶[16]。Wall等[17]又研究了一种余弦逆多项式方法，该方法设计下的轨道矢径可以有更多的极值。另外，Abdelkhalik和Taheri[18-19]还提出了利用有限傅里叶级数来解决未知机动的轨道重构问题的新方法。近几年，Xie等[20]提出了一种解决共面椭圆轨道间连续推力机动问题的简单形状近似方法，该方法利用多项式与椭圆轨道形状方程结合的复合函数完成轨道设计。本文研究的贝塞尔曲线是一种可以通过控制点灵活设计的曲线。学者们通过大量研究，应用贝塞尔曲线解决了机械手抓捕区域问题[21]，无人机路径规划问题[22]，轨迹生成问题[23]，弹道设计问题[24]，智能车轨迹跟踪[25]，车道检测[26]，工业无人运输轨迹规划[27]和多导弹协同航迹规划[28]等问题中的研究。然而，目前还没有研究将贝塞尔曲线用于航天器机动轨道设计和描述中。

本文首先给出了贝塞尔曲线的基本方程，并利用复合函数对机动轨道进行描述，给出了相应的约束条件和控制点设计，作为设计变量；然后，给出了平面切向机动的燃耗计算，作为优化指标；接着，针对推力峰值过大的问题，提出了多段式贝塞尔曲线法，并通过自由控制点的梯度下降修正解决了平面交会问题；最后相应的案例仿真表明新方法在连续推力机动轨道设计中优势明显。

1 贝塞尔曲线

轨道设计方法的核心是设计一个合理的函数作为机动轨道的轨道方程。而贝塞尔曲线法的轨道方程中包含有贝塞尔曲线方程，一个包含n+1个控制点的n阶贝塞尔曲线方程可表示为n阶Bernstein多项式的线性组合[22]，如式(1)所示。

(1)

式中：Pi=[ri，θi]为贝塞尔曲线的控制点参数，(1-s)n-isi为多项式部分，s为贝塞尔曲线的比例参数，当s从0到1改变时，贝塞尔曲线上的点对应的从起点移动至终点。

贝塞尔曲线表达式对s的各阶导数可表述为式(2)～(4)，其中C和A分别代表组合数和排列数。

(2)

(3)

(4)

由式(2)～(4)可得贝塞尔曲线矢径rB对相位方向θB的各阶导数如式(5)～(7)所示。

(5)

(6)

(7)

2 贝塞尔曲线法

本文通过将初始轨道方程、目标轨道方程与贝塞尔曲线结合构成的复合函数对机动轨道的表达式进行设计，并给出具体的转移任务起止点约束条件，然后完成控制点的设计以及燃耗计算。

2.1 复合函数设计

本文设计中复合函数设计为初始轨道方程H1、目标轨道方程H2和两个贝塞尔曲线B1，B2加权相加形式[20]，即转移轨道方程如式(8)所示：

r=rB1rH1+rB2rH2

(8)

机动轨道设计的物理特性实际上包含两部分：首先是需要满足轨道高度的逐渐提升或降低；其次是机动轨道必须满足转移过程中椭圆机动轨道的物理特性，即近地点和远地点的交替出现。而复合函数的设计中，贝塞尔曲线B1，B2的设计使得转移轨道可以灵活设计与优化，并为轨道高度提供了一个全局变化趋势；而H1，H2作为开普勒椭圆轨道方程则为转移轨道引入了椭圆轨道极坐标系下固有的震荡特性。当贝塞尔曲线参数s从0到1缓慢变化时，航天器相应的就从转移轨道的起点逐渐移动到终点，在轨道高度逐渐提升或降低的过程中伴随有椭圆轨道的周期震荡，从而更优地完成轨道转移任务。通过上述的复合函数矢径表达式可得r对θ的一阶导数和二阶导数如式(9)和式(10)所示：

(9)

(10)

2.2 约束条件

轨道转移任务要求转移轨道的位置参数r、位置对相位角一阶导数∂r/∂θ和位置对相位角二阶导数∂2r/∂θ2与初始轨道H1在初始位置θ0，目标轨道H2在终止位置θf匹配，则

(11)

将式(11)代入式(8)～ (10)，轨道转移任务的约束条件可改写为式(12)的形式。

(12)

2.3 控制点设计

在起止点处，贝塞尔曲线参数的n阶导数只受相邻的n个控制点影响，而轨道转移任务要求1阶、2阶导数匹配，即至少需要与起止点相邻的两个控制点被约束才能完成转移任务，因此贝塞尔曲线一共需要至少6个控制点。为了使设计结果更为灵活精确，本文额外引入一个自由控制点来提高贝塞尔曲线阶数，尝试获得更优的轨道设计结果。

首先给出6阶贝塞尔曲线B1方程为

(13)

6阶贝塞尔曲线r和θ关于s的各阶导数为

(14)

(15)

(16)

由式(12)可知，6阶贝塞尔曲线的第一个控制点和最后一个控制点参数如式(17)～(18)所示。

P11=[rB1|s=0,θB1|s=0]=[1,0]

(17)

P17=[rB1|s=1,θB1|s=1]=[0,θf]

(18)

由式(5)～(7)和式(14)～(15)可得，在s=0和s=1处，一阶导数如式(19)、(21)所示，二阶导数如式(20)、(22)所示。

(19)

(20)

(21)

(22)

引入关于控制点P12，P13，P15，P16的四个中间变量k12，k13，k15，k16，分别为

k12=θ12-θ11

(23)

k16=θ17-θ16

(24)

k13=θ11-2θ12+θ13

(25)

k15=θ15-2θ16+θ17

(26)

考虑如式(12)所示的约束条件，P12，P13，P15，P16的控制点参数可表示为

P12=[r12,θ12]=

(27)

P13=[r13,θ13]=

(28)

P15=[r15,θ15]=

(29)

P16=[r16,θ16]=

(30)

P14作为一个自由移动控制点，其参数为

P14=[r14,θ14]

(31)

同理，第二条6阶贝塞尔曲线B2的7个控制点P21，P22，P23，P24，P25，P26，P27的参数分别为

P21=[r21,θ21]=[0,θ0]

(32)

P22=[r22,θ22]=[0,θ21+k22]

(33)

P23=[r23,θ23]=[0,k23+2θ22-θ21]

(34)

P24=[r24,θ24]

(35)

P25=[r25,θ25]=[1,k25+2θ26-θ27]

(36)

P26=[r26,θ26]=[1,θ27-k26]

(37)

P27=[r27,θ27]=[1,θf]

(38)

至此，贝塞尔曲线法的控制点设计就已完成。这样，控制点的参数k12,k13,r14,θ14,k15,k16和k22,k23,r24,θ24,k25,k26将会作为后续轨道优化的优化变量。

2.4 相位角归一化

得到控制点的12个参数之后，最终的转移轨道可以通过式(8)获得。需要注意的是，当计算复合函数时，θB1必须与θB2一致才能做简单的叠加处理。然而，前面的两条贝塞尔曲线是独立设计的，r和θ都是由s求得，对于两条贝塞尔曲线同一个s值将会求出不同的θ。这里就需要做相位归一化处理，即贝塞尔曲线B2的s应由θB1反解获得，进而求得第二条贝塞尔曲线的参数rB2，这样就可以实现θ=θB1=θB2的统一。其中，由θB1反解s的方程可通过式(13)获得，具体表达式如式(39)所示。

c1s6+c2s5+c3s4+c4s3+c5s2+c6s1+c7s0=θB1

(39)

其中，系数矩阵可表示为

(40)

2.5 速度增量计算

前文给出了平面机动轨道设计优化问题中具体的优化变量，而本节将给出优化问题的指标函数，即累积速度增量ΔV。其具体推导和计算流程如下:

在小推力机动轨道形状近似的方法中，方程通常定义在极坐标系下。航天器在机动轨道起止点的速度和位置必须分别与初始轨道和目标轨道一致。

轨道形状方程为

(41)

由式(41)可求得rH关于相位角θ的各阶导数如下

(42)

(43)

(44)

航天器的速度可分解成两个分量：切向速度Vθ和径向速度Vr，如式(45)、(46)所示。

(45)

(46)

飞行航迹角γ和速度大小V由式(47)～(48)给出:

(47)

(48)

(49)

(50)

(51)

(52)

极坐标系下轨道动力学方程为

(53)

(54)

(55)

将式(55)代入式(53)并假设δ=γ，即推力方向与速度方向保持一致，δ为发动机安装角，则θ对t的一阶导数为

(56)

从式(47)～(49)和(56)可以看出，转移任务要求航天器在机动轨道起止点的位置r，位置对相位角的一阶偏导∂r/∂θ和位置对相位角的二阶偏导∂2r/∂θ2分别与初始轨道和目标轨道相匹配，这就是机动任务的约束条件。

相位角θ对时间t的二阶导数可由式(56)推导获得，表达式如下:

(57)

将式(56)～(57)代入式(47)和式(54)可得推力表达式[20]为

(58)

则相应地做积分可知，平面轨道转移时间Δt和累积速度增量ΔV分别为

(59)

(60)

以上得到的累积速度增量就是优化问题的指标函数Q=ΔV，再结合前面机动轨道设计给出的具体优化变量，接下来就可以利用优化算法通过迭代寻找最优机动轨道。

综上，平面转移轨道设计与优化的基本流程如图1所示，首先通过复合函数设计给出轨道形状方程的通式；接着，分别利用复合函数和控制点设计给出具体轨道方程和优化变量；最后，以累积速度增量为指标，通过优化给出最优结果。

图1 平面转移轨道设计与优化流程图Fig.1 The flow chart of coplanar transfer problem

3 多段贝塞尔曲线法

为了降低推力峰值，本文进一步将整个过程用多段贝塞尔曲线进行描述，分段点参数通过优化获得，完成综合最优机动轨道设计。

多段贝塞尔曲线法将在每条贝塞尔曲线中引入一个自由设计的分段点将曲线拆分为两条，且在分段点处二阶连续。机动轨道优化中指标函数计算与优化算法与前文相同，但在优化过程中额外添加Tmax<0.018的约束(该值为文献[20]设计结果中的最大推力峰值)。如表1所示，多段贝塞尔曲线法中包含四条贝塞尔曲线。其中两条曲线B1和B2分别由Ba，Bb和Bc，Bd组成，这样多段贝塞尔曲线法中的优化参数共计有32个。

4 贝塞尔曲线法交会轨道设计

前文利用贝塞尔曲线复合函数实现了机动轨道的设计与优化，而这里的机动轨道实际上是转移轨道，只要求航天器速度矢量和位置矢量满足要求，然而对机动时间没有任何约束。当考虑固定机动时间为tfixed，轨道转移问题就变成了轨道交会问题，为了满足这个额外的机动时间约束，本节就将在图1的基础上，通过在轨道设计和轨道优化过程之间引入梯度下降修正算法来实施修正，满足时间匹配约束，从而实现交会轨道的设计与优化。轨道机动时间与轨道周期直接相关，周期受轨道高度影响，而在本文中轨道高度又受自由控制点的矢径值影响，因此需要对自由控制点的矢径值进行修正以满足交会任务。

表1 多段贝塞尔曲线法控制点属性表Table 1 List of the control points for multi-segment method

图2 平面交会轨道设计与优化流程图Fig.2 The flow chart of coplanar rendezvous problem

如图2所示，交会轨道优化设计的具体流程与转移轨道设计十分相似，唯一的区别是额外的梯度下降修正算法部分。在转移任务中，机动轨道设计部分给出的是可行转移轨道，机动轨道优化部分则是在这些可行轨道中优化出最优转移轨道；而对于交会任务，首先要把机动轨道设计部分给出的可行转移轨道，通过梯度下降修正算法变成可行交会轨道，即实现机动时间匹配，然后再通过机动轨道优化部分优化出燃料最优的交会轨道。

当通过梯度下降对贝塞尔曲线法转移轨道设计的设计结果中两条贝塞尔曲线的自由控制点的矢径值r14和r24进行修正时，如图3所示，两个自由控制点的位置从P14和P24移动到P′14和P′24，对应的两条贝塞尔曲线也相应改变。通过这样的修正处理使得机动时间达到任务需要的值，从而实现交会任务，具体设计过程如下(1)～(4)所示。

(1)在转移任务每个个体确定之后，根据式(61)～(63)给出梯度下降算法的第一代个体，其中下标r和t分别代表交会任务和转移任务。

Δt(r)(1)=Δt(t)

(61)

r14(r)(1)=r14(t)

(62)

r24(r)(1)=r24(t)

(63)

(2)接着，通过式(64)～(65)，利用较小的数值偏差σ生成第二代个体，通常σ取1即可。

r14(r)(2)=r14(r)(1)+σ

(64)

r24(r)(2)=r24(r)(1)+σ

(65)

(3)根据r14(r)(2)和r24(r)(2)以及转移任务优化结果的其他参数计算可得Δt(r)(2)，进而利用如式(66)～(67)进行迭代生成之后的每一代个体，直到满足如式(68)所示的终止条件后停止，并将最后一代的ΔV(r)作为优化算法的指标函数。

r14(r)(n)=r14(r)(n-1)+

(66)

r24(r)(n)=r24(r)(n-1)+

(67)

|Δtfixed-Δt(r)(n)|<10-5

(68)

(4)最后，通过SA-DE-RM[29]算法的参数优化即可获得满足时间约束的最优交会轨道。

5 数值仿真

5.1 贝塞尔曲线法转移轨道设计仿真

优化过程中，每一代个体数设为10，迭代代数为100。贝塞尔曲线法的目标函数为minQ=ΔV，这是为了寻找燃耗最优的转移轨道。参数N代表转移所需轨道周期数，T代表每一时刻所需的推力，是一个没有任何约束的参数，完全由设计结果确定。贝塞尔曲线的前三个控制点的矢径参数r1，r2，r3均为1，而后三个控制点的矢径参数r4，r5，r6均为0，这是因为这6个控制点矢径参数均受式(12)所示的约束条件限制。设计结果中两条贝塞尔曲线的具体参数和变化趋势如表3和图4所示，其中k12，k13，r14，θ14，k15，k16和k22，k23，r24，θ24，k25，k26均为SA-DE-RM[29]算法的优化变量。在这些变量中，r14，θ14，r24，θ24是不受任何约束的自由变量。从图4可以看出，控制点P24的矢径参数r24为负数，这是因为它是一个完全自由的优化变量，不受任何约束限制。图5和图6分别给出了贝塞尔曲线法设计结果与文献[20]方法在极坐标系和直角坐标系下的轨道对比图，尽管两种方法都获得了一个可行的转移轨道，但设计结果差异明显。图7给出了贝塞尔曲线法推力变化曲线图，设计结果的累积速度增量为0.2044VU，而文献[20]方法中的设计结果为0.3231VU，可以看出贝塞尔曲线法的设计结果明显更优，设计结果相比于参考文献降低了36.7%。需要注意的是文献[30]方法在文献[20]方法的基础上通过轨道数值优化后的结果为0.2995VU，该结果虽略有提升，但仍远高于本文贝塞尔曲线法的设计结果。

表2 初始轨道和目标轨道参数Table 2 Parameters of initial orbit and final orbit

表3 贝塞尔曲线法控制点参数Table 3 Parameters of the control points for Bezier method

当转移周期数N从0到5改变时，多案例仿真下，贝塞尔曲线法总是可以获得一个优于文献设计结果的转移轨迹，累积速度增量降低了15%～40%，具体结果对比见表4。极坐标系和直角坐标系下轨道图如图8和图9所示。可以看出，转移过程中轨道总是保持在初始轨道和目标轨道之间，转移的起止点处，航天器的位置、速度都能与初始轨道和目标轨道完美匹配。图9中任意两个终止点间距都为2π，对应一个轨道周期。

图4 贝塞尔曲线法设计结果图Fig.4 The details of the Bezier curves

图7 推力变化曲线对比图Fig.7 The comparison of the thrust profiles

推力变化情况如图10所示。在N=0时，贝塞尔曲线法在累积速度增量ΔV和最大推力峰值Tmax方面均优于文献。而当N>0时，贝塞尔曲线法总是得到一个更小的累积速度增量ΔV，但同时推力峰值Tmax也更大，然而从图10可以看出，贝塞尔曲线法仅仅在推力峰值附近，推力更高，其他时刻推力均更小，因此最终的累积速度增量优于文献[20]。另外，需要注意，当N>3时，IP方法[14]无解，而本文方法和文献[20]方法均可获取有效的解。

图8 多案例仿真中极坐标系下轨道图Fig.8 Radius profiles when N=0～5

图9 多案例仿真中直角坐标下轨道图Fig.9 The transfer orbits when N=0～5

图10 多案例仿真中推力曲线图Fig.10 The thrust profiles when N=0～5

N贝塞尔曲线法文献[20]方法ΔVΔtTmaxΔVΔtTmax00.655312.89380.42590.884413.66870.980110.193768.76250.04770.288555.26030.029320.199194.75160.03900.3076110.28790.016730.2044170.47600.04030.3231162.88760.018040.2287191.71620.03660.3590215.35930.015850.2417263.46450.03150.4031267.82110.0132

5.2 多段贝塞尔曲线法转移轨道设计仿真

多段贝塞尔曲线法仿真中，参数设置均与表2相同，唯一区别是添加了最大推力峰值约束(仿真中取值为0.018DU/TU2)。多段贝塞尔曲线法设计结果参数如表5所示，图11给出了具体的四条贝塞尔曲线形式，仿真结果中，累积速度增量为0.1535VU，机动持续时间为169.8TU，最大推力峰值为0.013DU/TU2。与表4结果对比可以看出多段贝塞尔曲线法无论累积速度增量和最大推力峰值都远小于文献[20]方法，其中累积速度增量仅为文献的47.5%，这是因为多段法在降低最大推力峰值的同时，由于分段造成的控制点增多，实际相当于阶数也得以提升，因此设计过程也更加灵活，从而可以获得更优的累积速度增量。

图12、图13给出多段贝塞尔曲线法的轨道图。当本文在仿真中将进化代数增加至1000时，累积速度增量可达到0.135VU(文献[20]中为0.3231VU)，机动时间和最大推力峰值分别为178.75TU，0.0136DU/TU2。可以看到进一步优化的累积速度增量仅为文献[20]的41%，同时最大推力需求和机动时间也都更小，设计结果全面优于文献[20]方法。

表5 多段贝塞尔曲线法控制点参数Table 5 Parameters of the control points for multi-segment Bezier method

图11 多段贝塞尔曲线法设计结果图Fig.11 The details of the multi-segment Bezier curves

5.3 贝塞尔曲线法交会轨道设计仿真

当任务约束要求机动时间Δtfixed=180TU时，贝塞尔曲线法交会轨道设计与贝塞尔曲线法转移轨道设计的仿真结果对比如表6所示，可以看出交会任务结果满足机动时间约束。图14为最优交会轨道与最优转移轨道对比图，可以看出经过梯度下降修正之后，轨道发生了改变以满足交会任务时间需求。

图12 多段贝塞尔曲线法设计结果极坐标系下轨道图Fig.12 Radius profiles of multi-segment Bezier method

图13 多段贝塞尔曲线法设计结果直角坐标系下轨道图Fig.13 The transfer orbits of multi-segment Bezier method

图14 转移轨道和交会轨道直角坐标系下轨道对比图Fig.14 The comparison of transfer and rendezvous orbits

参数转移任务交会任务是否满足时间约束否是累积速度增量ΔV0.2044VU0.2059VU任务时间Δt170.4760TU180.0000TU最大推力需求Tmax0.0403VU/TU20.0401VU/TU2

6 结 论

本文提出了基于贝塞尔曲线的机动轨道设计方法。首先，将贝塞尔曲线的基本特点与轨道机动问题相结合，提出用贝塞尔曲线与轨道形状方程的复合函数描述和设计机动轨道的方法，并给出相应的约束条件，控制点设计和累计速度增量计算，所设计的控制点和计算得到的累计速度增量即为后续优化过程的优化变量与指标函数。仿真结果表明贝塞尔曲线法的设计结果远优于其他方法，但仍存在推力峰值过大的问题。进而，论文通过自由分段点的引入提出了多段式贝塞尔曲线法，在解决推力峰值过大问题的基础上，进一步优化了最优指标，最终燃耗仅为对比方法的41%。在此基础上，通过梯度下降算法对自由控制点进行修正来调整机动时间，解决了固定机动时间约束下的轨道交会问题。仿真结果证明本文提出的贝塞尔曲线法能更优地解决平面连续推力轨道转移和交会问题。