劉俊堯，賈建芳，李文強(qiáng)，宋錦武

(1. 中北大學(xué) 電氣與控制工程學(xué)院，山西 太原 030051； 2. 中北大學(xué) 信息與通信工程學(xué)院，山西 太原 030051；3. 廣州廣電計(jì)量檢測股份有限公司，廣東 廣州 510656)

0 引 言

航天器軌道轉(zhuǎn)移優(yōu)化技術(shù)是空間技術(shù)的重要組成部分. 目前對于軌道轉(zhuǎn)移優(yōu)化的研究大多基于脈沖變軌[1]. 在真實(shí)工作狀態(tài)下進(jìn)行軌道轉(zhuǎn)移時(shí)，由于發(fā)動(dòng)機(jī)推力制約，脈沖式變軌并不符合軌道轉(zhuǎn)移的實(shí)際情況[2]，特別是發(fā)動(dòng)機(jī)推力較小的航天器在遠(yuǎn)距離軌道轉(zhuǎn)移時(shí)，所需時(shí)間較長，存在很長的推力弧段，此時(shí)脈沖式變軌假設(shè)不再成立[3]. 與脈沖式軌道轉(zhuǎn)移技術(shù)相比，基于有限推力的軌道轉(zhuǎn)移技術(shù)能夠真實(shí)反映航天器的工作狀態(tài)，在提高研究結(jié)果真實(shí)性、 精確性方面有很大的優(yōu)越性.

基于有限推力的航天器軌道轉(zhuǎn)移優(yōu)化問題的本質(zhì)是泛函空間最優(yōu)控制問題[4]，根據(jù)優(yōu)化過程中對動(dòng)力學(xué)模型有限維近似處理方式不同，大致分為間接法、 直接法和混合法三種. 間接法[5-7]的本質(zhì)是以極小值原理以及非線性規(guī)劃作為理論基礎(chǔ)，通過相應(yīng)的邊界條件將軌道轉(zhuǎn)移優(yōu)化問題轉(zhuǎn)化為兩點(diǎn)邊值求解問題，雖能得到很光滑的解, 但由于對初始猜測異常要求很高，求解時(shí)困難較大； 直接法[8-10]的一般思路是對待優(yōu)化軌道進(jìn)行離散處理，利用優(yōu)化算法在每一個(gè)軌道段內(nèi)求解若干個(gè)待優(yōu)化參數(shù)，轉(zhuǎn)化后待優(yōu)化參數(shù)較多，需要在很大范圍內(nèi)進(jìn)行參數(shù)搜索，計(jì)算性能隨著節(jié)點(diǎn)數(shù)增加會有較大程度下降，算法魯棒性較差； 混合法[11-13]研究以優(yōu)化算法的實(shí)現(xiàn)為主要方向，對問題本身具有很大依賴性，求解時(shí)需要針對具體問題對算法做大量修改.

智能優(yōu)化算法[14-17]是一種啟發(fā)式優(yōu)化算法，一般是針對具體問題設(shè)計(jì)相關(guān)的算法，智能優(yōu)化算法技術(shù)性、 應(yīng)用性強(qiáng)，計(jì)算速度快且理論要求不高，得到國內(nèi)外學(xué)者廣泛關(guān)注并且發(fā)展迅速.

針對有限推力下遠(yuǎn)程導(dǎo)引變軌問題，本文采用改進(jìn)遺傳算法求解遠(yuǎn)程軌道能量最優(yōu)的多脈沖變軌策略，將求解得到的多脈沖變軌策略在變軌點(diǎn)轉(zhuǎn)換為推力弧段上的有限推力，通過發(fā)動(dòng)機(jī)推力計(jì)算出航天器加速度，對加速度進(jìn)行積分求得航天器瞬時(shí)速度增量，將速度增量不斷累加完成變軌，實(shí)現(xiàn)有限推力下遠(yuǎn)程導(dǎo)引變軌，并對有限推力下遠(yuǎn)程導(dǎo)引變軌策略速度增量方向進(jìn)行修正，修正后變軌策略精度有較大程度提高，與目標(biāo)軌道位置誤差大幅減小，能夠反映航天器在真實(shí)工作狀態(tài)下的軌道轉(zhuǎn)移情況.

1 遠(yuǎn)程導(dǎo)引變軌原理

1.1 特殊點(diǎn)變軌

根據(jù)特殊點(diǎn)變軌特點(diǎn)，變軌方案如圖1 所示.

圖1 特殊點(diǎn)變軌方法Fig.1 Special point change method

特殊點(diǎn)變軌大致分為四步：

1) 調(diào)整近地點(diǎn)高度： 初始軌道遠(yuǎn)地點(diǎn)處施加跡向脈沖量Δv1；

2) 調(diào)整軌道傾角： 緯度幅角為0或π處施加法向脈沖量Δv2，將異面軌道調(diào)整為共面軌道；

3) 調(diào)整遠(yuǎn)地點(diǎn)高度： 轉(zhuǎn)移軌道近地點(diǎn)處施加跡向脈沖量Δv3；

4) 調(diào)整軌道偏心率： 在目標(biāo)軌道遠(yuǎn)地點(diǎn)處施加跡向脈沖量Δv4，實(shí)現(xiàn)橢圓軌道向圓軌道轉(zhuǎn)換.

1.2 脈沖量計(jì)算

二體條件下，基于軌道動(dòng)力學(xué)方程，根據(jù)變軌前的軌道參數(shù)和變軌后軌道參數(shù)，可以求得每次變軌所需要脈沖量，每次施加脈沖量計(jì)算如下：

第1次變軌

(1)

第2次變軌

(2)

(3)

第3次變軌

(4)

第4次變軌

(5)

式中：n表示軌道角速度；e表示軌道偏心率；hp表示近地點(diǎn)高度；ha表示遠(yuǎn)地點(diǎn)高度； Δi表示軌道傾角修正量； ΔΩ表示升交點(diǎn)赤經(jīng)修正量；r表示航天器地心距；a表示軌道半長軸.

1.3 脈沖量轉(zhuǎn)換為有限推力方法

根據(jù)初始軌道參數(shù)和目標(biāo)軌道參數(shù)求解得到脈沖量，將脈沖量根據(jù)航天器發(fā)動(dòng)機(jī)推力、 比沖在發(fā)動(dòng)機(jī)工作時(shí)間段內(nèi)轉(zhuǎn)換為有限推力，將長時(shí)間變軌過程轉(zhuǎn)換為短時(shí)間內(nèi)的連續(xù)變軌，轉(zhuǎn)換方法如圖2 所示.

圖2 脈沖量向有限推力轉(zhuǎn)換Fig.2 Pulse thrust to finite thrust conversion

圖2 中，ti表示連續(xù)變軌過程中的時(shí)間節(jié)點(diǎn)，Δvi表示速度增量.

在軌道轉(zhuǎn)換過程中，需要確定發(fā)動(dòng)機(jī)開始工作時(shí)間、 終止工作時(shí)間以及發(fā)動(dòng)機(jī)產(chǎn)生的加速度.

發(fā)動(dòng)機(jī)工作時(shí)間為

(6)

式中：tb為發(fā)動(dòng)機(jī)開始工作時(shí)刻；tf為發(fā)動(dòng)機(jī)工作終止時(shí)刻；M為發(fā)動(dòng)機(jī)消耗燃料的質(zhì)量；F為發(fā)動(dòng)機(jī)提供的推力. 發(fā)動(dòng)機(jī)燃料消耗量M由航天器質(zhì)量、 發(fā)動(dòng)機(jī)比沖以及速度增量決定，計(jì)算公式為

(7)

式中：m為航天器質(zhì)量； Δv為速度增量；isp為發(fā)動(dòng)機(jī)比沖.

1.4 有限推力轉(zhuǎn)換修正

在脈沖變軌向有限推力變軌轉(zhuǎn)換過程中，在軌道轉(zhuǎn)移弧段對發(fā)動(dòng)機(jī)推力產(chǎn)生的航天器加速度進(jìn)行積分，求得節(jié)點(diǎn)處的速度增量. 在qsw本地軌道坐標(biāo)系中，X′軸是從地球中心指向衛(wèi)星的單位矢量；Z′軸是垂直衛(wèi)星軌道平面且方向向上的單位矢量；Y′軸是由(x,y,z)形成規(guī)則確定的單位矢量，如圖3 所示.

圖3 本地軌道坐標(biāo)系Fig.3 Local orbital coordinate system

圖3 中，O-XYZ是ECI坐標(biāo)系，O′-X′Y′Z′是qsw本地軌道坐標(biāo)系. 在qsw本地軌道坐標(biāo)系中，軌道為圓軌時(shí)，Y′軸方向?yàn)檐壍狼芯€方向，與航天器速度方向一致； 軌道為近圓軌道時(shí)，Y′軸方向與橢圓軌道切線方向存在少量偏差，與航天器速度方向基本保持一致.

在脈沖變軌向有限推力變軌轉(zhuǎn)換過程中，將qsw本地軌道坐標(biāo)系Y′軸方向近似為速度增量方向，通過確定節(jié)點(diǎn)處航天器速度矢量，計(jì)算速度增量方向.

2 遺傳算法

遺傳算法是一種高效、 并行的全局搜索方法. 它的求解過程是基于概率搜索技術(shù)，以目標(biāo)函數(shù)值作為搜索目標(biāo)，同時(shí)在多個(gè)搜索點(diǎn)對目標(biāo)進(jìn)行搜索. 標(biāo)準(zhǔn)遺傳算法全局最優(yōu)搜索能力強(qiáng)，但是在實(shí)踐應(yīng)用中，也存在局部搜索能力差的問題，無法保證算法完全收斂. 針對標(biāo)準(zhǔn)遺傳算法存在的問題，對遺傳算法的選擇和交叉操作進(jìn)行改進(jìn)，在對群體根據(jù)適應(yīng)度值高低進(jìn)行排序的基礎(chǔ)上，用個(gè)體最優(yōu)與其他偶數(shù)位的所有個(gè)體進(jìn)行交叉，每次交叉產(chǎn)生兩個(gè)新的個(gè)體，在交叉過后，對新產(chǎn)生的群體進(jìn)行多點(diǎn)變異產(chǎn)生子群體，再計(jì)算其適應(yīng)度值，然后與父群體合并，并且根據(jù)適應(yīng)度值進(jìn)行排序，取前NP個(gè)個(gè)體為新群體，進(jìn)行下一次遺傳操作. 改進(jìn)遺傳算法運(yùn)算流程如圖4 所示.

二體條件下，式(1)中的n1、e1、hp由初始軌道根數(shù)確定； 式(3)中的n3、e3由初始軌道根數(shù)及Δv1確定，h′a為目標(biāo)軌道遠(yuǎn)地點(diǎn)高度，ha為初始軌道遠(yuǎn)地點(diǎn)高度； 式(5)中的n4、a4、e4由目標(biāo)軌道根數(shù)確定. 選取第1次變軌時(shí)轉(zhuǎn)移軌道近地點(diǎn)高度h′p為唯一優(yōu)化變量，參數(shù)設(shè)置如下： 基因數(shù)目D=1； 染色體數(shù)目NP=100； 最大遺傳代數(shù)G=1 000； 交叉概率Pc=0.8； 變異概率Pm=0.1. 構(gòu)建目標(biāo)函數(shù)

fmin=Δv=∑(|Δv1|+|Δv2|+|Δv3|).

(8)

圖4 改進(jìn)遺傳算法流程圖Fig.4 Improved genetic algorithm flow chart

改進(jìn)遺傳算法求解h′p最優(yōu)解時(shí)，先計(jì)算Δv1、 Δv2、 Δv3的適應(yīng)度并按升序排列，選擇適應(yīng)度最高的Δv1、 Δv2、 Δv3染色體個(gè)體與其他偶數(shù)位的所有染色體個(gè)體進(jìn)行交叉并進(jìn)行多點(diǎn)變異產(chǎn)生子代染色體Δv′1、 Δv′2、 Δv′3，計(jì)算子代染色體適應(yīng)度并按升序排列，將Δv′1、 Δv′2、 Δv′3染色體與Δv1、 Δv2、 Δv3染色體合并，計(jì)算合并后的適應(yīng)度并按升序排列，選取前NP個(gè)適應(yīng)度最高值作為此代最優(yōu)個(gè)體.

3 仿真算例與分析

3.1 軌道根數(shù)設(shè)置

初始軌道根數(shù)和目標(biāo)軌道根數(shù)設(shè)置如表1. 表1 中，a表示軌道半長軸；e表示偏心率；i表示軌道傾角；Ω表示升交點(diǎn)赤經(jīng)；ω表示近地點(diǎn)幅角；anv表示真近點(diǎn)角.

由于初始軌道與目標(biāo)軌道共面，無需改變軌道傾角，僅需要施加3次脈沖量，實(shí)現(xiàn)遠(yuǎn)程導(dǎo)引變軌.

表1 軌道根數(shù)參數(shù)

3.2 仿真結(jié)果與分析

以第一次變軌時(shí)轉(zhuǎn)移軌道近地點(diǎn)高度h′p為優(yōu)化變量，采用改進(jìn)遺傳算法求解總速度增量的最小值. 結(jié)果顯示，當(dāng)h′p為6 982.039 km時(shí)，總速度增量最小值為66.699 m/s，其中第1次變軌的速度增量Δv1為28.855 m/s，第2次變軌的速度增量Δv2為33.122 m/s，第3次變軌的速度增量Δv3為 4.722 m/s，航天器遠(yuǎn)程導(dǎo)引變軌序列如表2 所示.

表2 脈沖變軌方案

將初始得到的脈沖變軌方案在變軌點(diǎn)處轉(zhuǎn)換為推力弧段上的有限推力，設(shè)置航天器初始質(zhì)量為1 000 kg，發(fā)動(dòng)機(jī)比沖為220 s，分別施加50 N, 100 N, 200 N的推力，轉(zhuǎn)換結(jié)果如表3 所示.

表3 50 N推力下變軌方案

表4 100 N推力下變軌方案

由表3～表5 可知，有限推力變軌方案消耗的燃料會多于脈沖變軌方案； 推力越小，軌道上推力弧段越長，變軌結(jié)果誤差越大. 由表5可知，推力為200 N時(shí)，航天器在經(jīng)過3次變軌后，軌道半長軸為6 998.295 km，與目標(biāo)軌道半長軸誤差為1.715 km，相對位置誤差為3.89 km. 轉(zhuǎn)換后的連續(xù)推力變軌方案速度不斷迭加，在不同節(jié)點(diǎn)處不斷變化，但速度增量方向一直保持脈沖變軌處方向不變，對轉(zhuǎn)換結(jié)果造成較大影響，并且變軌持續(xù)時(shí)間，誤差越大. 對200 N推力下變軌方案進(jìn)行修正，結(jié)果如表6 所示.

表5 200 N推力下變軌方案

表6 有限推力變軌方案修正結(jié)果

由表6 可知，有限推力變軌轉(zhuǎn)換修正方案精度更高，與目標(biāo)軌道半長軸誤差為0.957 km，相對位值誤差為1.33 km. 經(jīng)過3次變軌，軌道根數(shù)與目標(biāo)軌道根數(shù)對比如表7 所示.

表7 軌道根數(shù)對比

表7 中，真近點(diǎn)角差異是由于推力弧段的存在，導(dǎo)致變軌結(jié)束時(shí)刻有差異. 由表7 可知，脈沖變軌方案轉(zhuǎn)換為有限推力方案，僅有軌道半長軸、 偏心率、 真近點(diǎn)角存在少量差異，軌道傾角、 升交點(diǎn)赤經(jīng)、 近地點(diǎn)幅角保持一致； 有限推力修正方案精度更高，可以反映航天器在真實(shí)飛行環(huán)境中的變軌情況. 航天器有限推力方案變軌仿真結(jié)果如圖5 所示.

由圖5 可知，變軌后軌道高度提升，軌道仍為圓軌； 由圖6 可知，航天器軌道半長軸在機(jī)動(dòng)時(shí)間段內(nèi)不斷增大，非機(jī)動(dòng)時(shí)間段內(nèi)基本保持不變，經(jīng)過3次變軌，軌道半長軸與目標(biāo)軌道半長軸高度一致.

圖5 有限推力方案3D仿真圖Fig.5 Finite thrust scheme 3D simulation diagram

與脈沖變軌相比，有限推力變軌所需燃料更多，是由于有限推力變軌弧段存在引力損失； 不進(jìn)行推力方向修正的變軌方案的誤差要高于進(jìn)行推力方向修正的變軌方案，是因?yàn)橥ㄟ^修正推力方向可找到時(shí)間節(jié)點(diǎn)處推力的方位角和俯仰角.

4 結(jié) 論

本文對有限推力下遠(yuǎn)程導(dǎo)引變軌進(jìn)行了研究. 通過改進(jìn)遺傳算法的選擇和交叉操作，求解得到燃料最優(yōu)的多脈沖變軌方案，在脈沖點(diǎn)處將脈沖量轉(zhuǎn)換為推力弧段上的有限推力，實(shí)現(xiàn)多脈沖機(jī)動(dòng)遠(yuǎn)程導(dǎo)引變軌向有限推力下遠(yuǎn)程導(dǎo)引變軌的轉(zhuǎn)換，并對有限推力變軌方案速度增量方向進(jìn)行修正，修正結(jié)果精度有較大程度提高，證明方法確實(shí)有效. 基于有限推力的遠(yuǎn)程導(dǎo)引變軌方案可以反映航天器在真實(shí)工作狀態(tài)下環(huán)境中變軌的過程，并取得較高的精度.