鄭紅星，周思雨，李偉杰

（1. 哈爾濱工業(yè)大學航天學院，哈爾濱150001；2. 北京空間飛行器總體設(shè)計部，北京100094）

1 引 言

目前各航天大國紛紛提出以“在軌服務(wù)”的方式提高航天器在軌機動能力和工作壽命，升級擴大航天器功能，創(chuàng)造經(jīng)濟效益[1]。美國已經(jīng)開展了XSS[2]、DART[3]、軌道快車[4]等飛行器試驗項目，進行貨物運輸、模塊更換、地面遙操作等技術(shù)驗證和預(yù)先研究[5]。在“在軌服務(wù)”任務(wù)中，當一顆服務(wù)星對多顆在軌運行的目標星進行服務(wù)時稱為“一對多”服務(wù)方式。國外研究表明，“一對多”服務(wù)方式優(yōu)于“一對一”服務(wù)方式，并將成為未來在軌服務(wù)的主流方式[6]。

國內(nèi)外已經(jīng)針對在軌服務(wù)的軌道機動方式和自主決策規(guī)劃等展開了深入的研究。文獻[7]針對未來服務(wù)型航天器總結(jié)了漂移式、組合式、直接彈道式3 種共軌式機動模型，并給出了機動時間、能量、相位的關(guān)系式與仿真曲線。文獻[8]針對圓軌道航天器在軌燃料加注任務(wù)，提出了一種基于聚類分析的在軌加注任務(wù)調(diào)度及優(yōu)化算法。文獻[9]對運行在地球同步軌道上且傾角較小的衛(wèi)星群的“一對多”在軌加注任務(wù)開展了相關(guān)研究，文中將共面轉(zhuǎn)移與變軌道面機動所消耗的燃料進行了比較，并在小角度近似下將該任務(wù)規(guī)劃問題簡化成旅行商（TSP）問題進行求解。文獻[10]針對在軌服務(wù)任務(wù)規(guī)劃問題，建立了一套系統(tǒng)的建模與優(yōu)化方法，深入研究了在軌服務(wù)任務(wù)規(guī)劃模型建模方法。尤其針對在軌加注任務(wù)，提出了一套可應(yīng)用于P2P、“一對多”和混合模式的在軌加注任務(wù)規(guī)劃問題中的系統(tǒng)建模與優(yōu)化方法。許多文獻均將研究重點放在對地球靜止軌道（GEO）衛(wèi)星的服務(wù)上，這是由于GEO 軌道高度較高，發(fā)射衛(wèi)星的代價較大，如果能夠利用發(fā)射到低軌的服務(wù)星為其提供在軌服務(wù)，延長其壽命，將更有利于減小成本[11]。

本文也將重點研究針對GEO軌道上“一對多”的在軌服務(wù)任務(wù)，分別考慮服務(wù)星與工作星軌道共面與異面情況，服務(wù)星的軌道機動采用脈沖變軌的方式，重點通過求解Lambert 問題得到有限轉(zhuǎn)移時間約束下的機動速度增量，以此作為優(yōu)化指標與遺傳算法序列規(guī)劃相結(jié)合，得出使得總脈沖速度增量最小的服務(wù)策略。

2 數(shù)學模型

2.1 軌道轉(zhuǎn)移任務(wù)

為了簡化問題我們忽略攝動等因素的影響，將軌道轉(zhuǎn)移簡化成只考慮地球一個中心天體的二體問題。一個服務(wù)星可以為同一軌道上的多顆工作星提供服務(wù)，即“一對多”。首先介紹三種軌道轉(zhuǎn)移的方式。

2.1.1 霍曼轉(zhuǎn)移

在同一軌道平面內(nèi)，兩不相交圓軌道之間的過渡霍曼轉(zhuǎn)移為能量最省的方式?；袈^渡軌道是一條外切于小圓軌道，內(nèi)切于大圓軌道的橢圓軌道，由兩次脈沖機動實現(xiàn)[1]。

兩次機動的脈沖大小分別為：

2.1.2 異面圓軌道轉(zhuǎn)移

對于沒有時間約束的兩個軌道半徑不同的異面圓軌道之間的轉(zhuǎn)移，本文采用了先在初始軌道平面內(nèi)進行霍曼轉(zhuǎn)移將軌道轉(zhuǎn)移到與新軌道半徑相同的異面軌道上，這里進行兩次機動。然后在兩個異面軌道的交點再施加一次脈沖，用于改變軌道平面，完成轉(zhuǎn)移。

相交的同軌道半徑異面軌道機動速度：

2.1.3 Lambert轉(zhuǎn)移

Lambert 問題是航天動力學中的兩點邊界值問題，主要研究如何根據(jù)兩個位置矢量和轉(zhuǎn)移時間確定航天器轉(zhuǎn)移軌道[1]。根據(jù)初始時刻服務(wù)星的軌道六根數(shù)可以得到其位置和速度(r10、v10)，工作星的軌道六根數(shù)可以得到同一時刻工作星的(r20、v20)。給定轉(zhuǎn)移時間t，迭代求解出工作星在t時刻后的位置速度(r2f、v2f)，經(jīng)過t時間后服務(wù)星也要到達同一位置。已知服務(wù)星初始位置、轉(zhuǎn)移時間和最終位置后，就可以采用轉(zhuǎn)化為Lambert 兩點邊值問題求解出初始速度v1t與最終速度v1f。由此可以得到兩次的脈沖速度：Δv1=v1t-v10和Δv2=v2t-v1f。

2.2 任務(wù)規(guī)劃

該服務(wù)任務(wù)規(guī)劃可以理解為服務(wù)星選擇合適的路徑來對工作星進行在軌服務(wù)，它將每顆工作星服務(wù)一次后回到自身軌道，使得任務(wù)完成后消耗的總推進劑最少。該問題實際上就是不對稱的TSP。從圖論的角度來看，該問題實質(zhì)是在一個帶權(quán)完全無向圖中，找一個權(quán)值最小的回路。由于該問題的可行解是所有頂點的全排列，隨著頂點數(shù)的增加，會產(chǎn)生組合爆炸，它是一個非確定性多項式完全（Nondeterministic Polynomial Complete，NPC）問題。

3 分析與討論

3.1 同面軌道服務(wù)

考慮一個服務(wù)星為三個工作星提供在軌服務(wù)，三顆工作星間隔120°均布GEO 軌道平面上。服務(wù)星工作在10000km軌道上。如圖1分別為工作星軌道與服務(wù)星軌道共面和異面兩種情況。

圖1 同面軌道轉(zhuǎn)移（左）與異面軌道轉(zhuǎn)移（右）Fig.1 Same-plane orbit transfer(left)and different-plane orbit transfer(right)

將發(fā)出在軌服務(wù)的指令的時刻記為初始時刻，假設(shè)此時三顆工作星的軌道六根數(shù)分別為：

服務(wù)衛(wèi)星的軌道六根數(shù)為：

假設(shè)轉(zhuǎn)移時間tf= 7200s。為方便后面的任務(wù)，采用整數(shù)編碼，服務(wù)衛(wèi)星編號0，工作星分別編號1、2、3。

服務(wù)星依次軌道轉(zhuǎn)移到第一顆工作星，第二顆工作星再到第三顆工作星，結(jié)束服務(wù)后再返回原軌道。前幾次服務(wù)通常由于存在時間限制采用Lambert軌道轉(zhuǎn)移，最后一次返回原軌道不存在時間限制，所以采用霍曼轉(zhuǎn)移即可。仿真結(jié)果計算見表1。

MATLAB仿真得到的部分結(jié)果如圖2。藍色為工作星運行軌道，綠色為服務(wù)星運行軌道，紅色虛線為轉(zhuǎn)移軌道。轉(zhuǎn)移初始時刻服務(wù)星位置為綠點，工作星初始位置藍點，轉(zhuǎn)移結(jié)束后同時到達紅點處。

表1 同面機動速度Table1 Same-face maneuvering speed

3.2 異面軌道服務(wù)

異面時的情況如圖1右圖所示。服務(wù)指令發(fā)出時刻假設(shè)三顆工作星的軌道六根數(shù)與3.1 節(jié)中相同，服務(wù)衛(wèi)星軌道六根數(shù)為：

圖2 同面軌道轉(zhuǎn)移仿真結(jié)果圖Fig.2 Simulation results of same-plane orbit transfer

同樣假設(shè)服務(wù)星依次軌道轉(zhuǎn)移到第一顆工作星，第二顆工作星再到第三顆工作星，結(jié)束服務(wù)后再返回原軌道。前幾次服務(wù)采用Lambert 軌道轉(zhuǎn)移，最后一次返回無時間限制，所以采用先霍曼再異面的轉(zhuǎn)移方式。仿真結(jié)果計算如表2。

表2 異面機動速度Table2 Different-plane maneuvering speed

圖3 異面軌道轉(zhuǎn)移仿真結(jié)果圖Fig.3 Simulation results of different-plane orbit transfer

3.3 基于遺傳算法的規(guī)劃

根據(jù)3.1和3.2 節(jié)的分析，已得到同面變軌和異面變軌的模型，并計算得出0 ?1、0 ?2、0 ?3、1 ?2、2 ?3、1 ?3 的機動速度，接下來需要求解使得機動速度總和最小的在軌服務(wù)方案。該問題隸屬于整數(shù)規(guī)劃問題，解空間為離散變量，模型的梯度等信息不存在。針對此類問題的特點，遺傳算法具有較強的優(yōu)勢，它通過模仿物種群中優(yōu)勝劣汰的選擇機制以及種群中優(yōu)勢個體的繁衍進化來實現(xiàn)優(yōu)化功能[1]，因此本文采用遺傳算法求解。算法流程如圖4所示。

3.3.1 求解過程

圖4 算法流程圖Fig.4 Algorithm flowchart

隨機產(chǎn)生一個種群作為初始種群，同時計算這個初始種群的適應(yīng)度——取機動速度總和的倒數(shù)作為適應(yīng)度，適應(yīng)度越高，被保留概率越大。值得注意的是，由于整個在軌服務(wù)的模型需從服務(wù)星開始機動，所以種群的首序號不為0時，適應(yīng)度函數(shù)為無窮大。每次完成交叉、變異運算之后也需要重新評價。

交叉運算操作如圖5和圖6所示，選擇兩個個體，在[1，sz]（sz為服務(wù)星和工作星的總數(shù)量）的范圍內(nèi)，隨機產(chǎn)生兩個交叉位min和max（min

圖5 交叉前示意圖Fig.5 Schematic diagram before crossing

圖6 交叉后示意圖Fig.6 Schematic diagram after crossing

其中，A、B是未發(fā)生交叉之前的染色體（個體），A’和B’是A、B 發(fā)生交叉之后產(chǎn)生的新個體。通過沖突檢測，發(fā)現(xiàn)交叉之后會產(chǎn)生同一個基因在同一條染色體上重復(fù)出現(xiàn)，此為交叉的沖突，如圖6粉色標注部分。對圖6消除基因沖突的示意圖如圖7~10所示。

變異運算操作首先在隨機產(chǎn)生兩個變異位min、max，其中min

圖7 將染色體A沖突基因8改為染色體B對應(yīng)位置的5Fig.7 Change chromosome A conflict gene 8 to 5 corresponding to chromosome B

圖8 將染色體A沖突基因5改為染色體B對應(yīng)位置的7Fig.8 Change chromosome A conflict gene 5 to 7 corresponding to chromosome B

圖9 將染色體A沖突基因7改為染色體B對應(yīng)位置的4Fig.9 Change chromosome A conflict gene 7 to 4 corresponding to chromosome B

圖10 對染色體B同樣進行消除沖突Fig.10 Conflict elimination on chromosome B

圖11 變異結(jié)果Fig.11 Variation results

算法的終止條件是達到最大代數(shù)的迭代次數(shù)，每一次迭代結(jié)束后將得到的路徑長度和當前代數(shù)（第幾代）記錄在數(shù)組中，然后在搜索完成之后，將數(shù)組中記錄的最短路徑和對應(yīng)的代數(shù)輸出，作為搜索結(jié)果。

3.3.2 結(jié)果分析

求解TSP 的遺傳算法是利用MATLAB 軟件編程的。由于遺傳算法包含了隨機搜索方法，所求的最優(yōu)解不一定是實際最優(yōu)的。在求解過程中，發(fā)現(xiàn)遺傳算法得到的結(jié)果的精確度除了和交叉算子、變異算子、適應(yīng)度計算方法有關(guān)，還受交叉概率、變異概率、迭代次數(shù)的影響。對于本文算法，在一定范圍內(nèi)，迭代次數(shù)也大，變異概率越小，遺傳算法的精確度越高；執(zhí)行時間隨著迭代次數(shù)的增加而增加。當交叉概率為0.8，變異概率為0.5，最大代數(shù)為10000 時，能得到較理想的結(jié)果，如圖12所示。

圖12 搜索過程圖Fig.12 Search process diagram

遺傳算法得到的最小機動速度總和為44.00km/s，最短路線為0 →1 →2 →3 →0。

4 結(jié) 論

本文忽略軌道攝動等影響因素簡化軌道轉(zhuǎn)移問題，對“一對多”模式的在軌服務(wù)任務(wù)規(guī)劃問題進行了研究。分別考慮了服務(wù)星與工作星軌道共面與異面的情況，假設(shè)軌道機動采用脈沖變軌的方式，針對不同的約束條件通過求解Lambert 問題、霍曼轉(zhuǎn)移或異面圓軌道轉(zhuǎn)移得到相應(yīng)的機動速度增量，并將總機動速度增量作為優(yōu)化指標，采用遺傳算法序列規(guī)劃得出使得總脈沖速度增量最小的服務(wù)策略。