微納衛星星座的Kuhn-Munkres匹配部署優化方法

劉思陽,蒙 濤,2,雷家坤,金仲和,2

(1. 浙江大學微小衛星研究中心，杭州 310027; 2. 浙江省微納衛星研究重點實驗室，杭州 310027)

0 引 言

微納衛星星座因成本低、研制周期短、構型和功能靈活等特點獲得了廣泛的應用空間。為降低發射成本，一箭多星發射成為星座發射與部署組網的主要方式。據統計，2019年全球共進行102次航天發射，單發火箭發射航天器的平均數量接近5，包括一箭60星發射的StarLink衛星以及一箭34星發射的OneWeb衛星。一箭多星發射要求衛星通過軌道控制進入目標軌位，對微納衛星而言仍有許多問題亟待解決。首先，微納衛星部署希望盡量減少燃料預算并使各衛星燃耗量均衡，從而提高星座整體的長期軌控能力。其次，微納衛星面臨推力大小、推力器開機時長受限的問題，傳統大衛星的軌控方案并不適用。同時，為使星座快速投入在軌服務，對部署時間的把控也同樣重要。此外，實際中受發射順序、發射時間差、發射速度等不確定因素的影響，或衛星在部署之前已經歷較長時間在軌而產生位置漂移，許多情況下衛星初始位置在不同的軌道面而并非集中到空間一點。因此，一箭多星發射背景下的微納星座部署問題應視為以燃料消耗和衛星燃耗量均衡性為優化目標，以部署時間、推進能力和始末狀態為約束的綜合優化問題。本文研究一箭多星發射背景下一個軌道面的優化部署，將具有軌道高度、升交點赤經、緯度輻角三維位置差異的衛星群部署到同一軌道面不同目標相位。

星座部署可劃分為部署規劃和部署控制兩階段。部署規劃解決衛星目標位置的分配問題，在以往研究中常隨機指定或依靠簡單的經驗規則實現。因為當衛星初始位置簡化為一點、部署精度要求不高或燃料充足時，衛星目標位置的分配對系統的影響微小。部署控制解決已知始末狀態的軌道控制問題，關鍵在于節約燃料，包括異軌部署控制和同軌部署控制。異軌部署控制一般通過多次發射并配合發射窗口實現，或借助環境攝動實現軌道面的緩慢分離[1-2]。例如，六顆微納衛星組成的FORMOSAT-3/COSMIC星座利用地球扁率攝動下不同高度的軌道升交點赤經漂移速率不同的特征，于2006年4月完成軌道面分離[3]。文獻[4]通過半長軸調整間接控制全球星的升交點回歸經度。文獻[5]研究了利用J2攝動和微推力作用實現半長軸、升交點赤經和傾角調整的四種策略及對應燃耗。此類方法以部署時間為代價，降低了多次發射的發射成本。

同軌部署控制的重點是相位控制，最基本的方法是依靠半長軸差實現相位調整。文獻[6]通過調整半長軸間接調整相對相位進而實現軌道面內的星座相位控制。文獻[7]采用兩次切向機動軌道轉移完成相位調整。文獻[8]根據相同原理提出一種基于李雅普諾夫的集群長期相對相位控制策略。以上方法僅關注相位調整效果，忽略升交點赤經漂移和共面控制。依靠氣動阻力差、太陽光壓差及其與軌道機動的配合作用實現相位控制，可極大的減少燃料消耗[9]。文獻[10]以氣動阻力差為軌道調整的唯一手段，研究了CYGNSS星座的部署策略，獲得同軌8星均勻分布的構型。文獻[11]調整展開機構姿態獲得阻力差，實現同軌衛星相位分離，在仿真中將這一方法應用到Flock星座中部署。利用差分阻力實現無推進能力衛星的軌道修正和衛星位置保持[12-14]的案例也可為星座部署提供參考。但此類方法對于姿態或展開機構具有固定指向需求的衛星并不適用。此外，依靠軌道機動可實現快速入軌，此類方法的關鍵是軌跡與燃料的最優化問題[15-16]。例如，文獻[17]對燃料和時間最優路徑規劃與軌道機動策略進行研究。文獻[18]采用狀態依賴Riccati方程方法給出了無徑向推力的交會對接和編隊星重構的最優軌道控制。此類方法適用于轉移距離較短、燃料充足、對控制精度和時效要求較高的衛星，不局限于同軌部署，因而在編隊捕獲和交會對接中更常見。星座的同軌部署涉及相位分離和共面控制，但目前關于同軌相位分離的研究將所有衛星初始位置近似到一個點且忽視共面控制，即部署過程中的升交點赤經漂移。這種簡化無法評估停泊階段的相對升交點赤經漂移，需浪費大量額外燃料進行升交點赤經修正。

本文將軌道高度、升交點赤經、緯度輻角均不相同的衛星群部署到同一軌道面不同目標相位，這一場景可看作星座的同軌部署，也可看作“空間-平面”的構型變換任務，具有更強的普適性。部署規劃方面，利用衛星初始位置差異和停泊階段相對漂移可在J2攝動的作用下相互補償的基本思想，基于衛星位置可互換的假設，首次嘗試將Kuhn-Munkres(KM)算法引入星座部署并作適應性改進實現衛星與目標位置的優化匹配，節約部分燃料同時提高各星燃耗量的均衡性。對于傾角較低、目標相位分布較廣、初始升交點赤經處于±0.7°范圍內的星群優化效果突出。部署控制方面，利用軌道高度差和J2攝動實現相位分離的同時盡量修正升交點赤經。采用切向N次有限常值推力機動實現軌道機動，適用于推力大小受限，開機時長受限的微納衛星推進系統。

1 星座部署優化問題建模

初始時刻，各星半長軸、升交點赤經、緯度輻角各不相同。為把各星部署至同一軌道面，需要進行緯度輻角、升交點赤經和半長軸調整，如圖1所示。由于J2攝動對偏心率和傾角無長期影響，無法利用J2攝動對二者進行修正，在本文的討論中假設各星傾角和偏心率相同。低軌微納衛星星座一般運行在圓軌道上，本文偏心率取0。

圖1 衛星初始位置與目標位置示意圖

J2攝動下的平均軌道要素攝動方程[19]為：

(1)

式中：a為軌道半長軸;e為偏心率;i為傾角;Ω為升交點赤經;w為近地點輻角;M為平近點角;n為平均角速度;J2為攝動常數;Re為地球半徑。

1.1 緯度輻角調整

(2)

緯度輻角調整所需時間為:

(3)

1.2 升交點赤經調整

圓軌道中，J2攝動引起的升交點赤經漂移速率與平均角速度、半長軸和傾角的關系為:

(4)

由此可利用兩星半長軸差間接調整相對升交點赤經[4]。此外，升交點赤經還可由法向推力直接修正(以下稱直接修正法)。修正ΔΩ所需速度增量為[20]:

(5)

(6)

法向推力也會改變傾角，傾角改變量Δi與法向推力ΔVi的關系為:

(7)

(8)

(9)

式中：ΔVΩ+i為聯合修正燃耗;i0為初始傾角;u1,u2為兩次機動的緯度輻角。

當小推力衛星無法一次性輸出全部速度增量時，可以在最佳修正位置附近適當延長推力器開機時間(總燃耗將大于ΔVmin)，或者在最佳修正位置附近開機較短時間，每軌開機兩次，直至速度增量輸出完畢(總燃耗約為ΔVmin)。為方便計算和比較，下文在升交點赤經修正燃耗的理論計算中均取最小燃耗。在數值仿真中，取修正時間和燃耗的折中方案，即在最佳修正位置附近延長開機時間。

1.3 半長軸調整與推力模型

切向N次有限常值推力機動可實現半長軸調整[20]，對推力大小和開機時長無特殊要求。期間衛星沿特定弧線飛行至目標位置，飛行軌跡由N段推力作用區間和N-1段自由飛行區間組成，可配合推力器開機-關機時間執行。設徑向轉移距離為Δz，則每段推力作用區間的速度增量為:

(10)

總速度增量為:

(11)

式中：n為目標軌道平均角速度。轉移總時間為:

ΔT=TonN+Toff(N-1)

(12)

軌道轉移期間衛星航向轉移距離和轉移角分別為:

(13)

(14)

根據半長軸與偏心率聯合調整方法[20]，在緯度輻角為u和u+π時執行兩次速度增量相等的軌道機動，半長軸調整不改變偏心率。實際執行時可在兩個機動點附近延長開機弧段，每軌執行兩次對稱的軌道機動，配合切向N次有限常值推力機動直至ΔVa輸出完畢。

推力器連續開機時間為Ton，推力間隔時間為Toff，可輸出的最大推力為tmax，則T時刻實際推力t為:

(15)

若推力器持續工作，可將推力間隔時間Toff取0，對應連續常值推力機動。

1.4 星座部署優化模型

半長軸、緯度輻角、升交點赤經三者之間存在耦合關系[21]，緯度輻角和升交點赤經可由軌道高度間接控制。因此可通過優化停泊軌道高度以及衛星與目標緯度輻角之間的匹配關系，盡量使各星的升交點赤經和緯度輻角得到同步控制。當二者沖突時，優先確保緯度輻角的控制精度，最后由直接修正法修正升交點赤經余量。

(16)

(17)

sk星軌道高度轉移所需速度增量為:

(18)

sk星總速度增量(總燃耗)為:

(19)

K個衛星的平均燃耗為:

(20)

燃耗量均衡性可由各星總速度增量的方差評估:

(21)

單星部署時間由停泊軌道運行時間和兩次軌道轉移時間組成:

(22)

(23)

TonNm2+Toff(Nm2-1)

(24)

T=maxTk

(25)

(26)

2 星座部署優化問題求解

鑒于星座部署問題可利用許多先驗知識來提高求解質量和效率，本文通過理論分析逐步確定優化變量或盡量縮小優化變量的取值范圍，將部署優化問題簡化為衛星與目標緯度輻角(目標位置)的匹配。將其看作帶權值的二分圖匹配問題，采用KM算法求最優解。該算法簡單易實現，無需迭代，能夠快速給出優化結果。

2.1 目標軌道與停泊軌道的確定

從節約燃料的角度考慮，初始軌道較高的衛星期望較高的停泊軌道(與目標軌道相比)；初始升交點赤經較小的衛星也期望較高的停泊軌道以便借助軌道高度差修正更多的升交點赤經余量；由于目標緯度輻角未確定，相位對停泊軌道高度的期望不明確。因此，當初始軌道高度和升交點赤經對停泊軌道高度期望一致時，可縮小停泊軌道高度的取值范圍；當二者沖突時則需遍歷全部備選停泊軌道。在此基礎上挑選目標相位，使得三者盡量依靠自然運動和最小燃耗得到同步調整。引入期望因子F來描述衛星對停泊軌道高度的偏好：

(27)

適合衛星sk的停泊軌道半長軸取值范圍為:

(28)

式中：amax為軌道調整量約束。據此可先計算每個衛星對于每個目標位置的最佳停泊軌道高度及最小燃耗，再經過綜合優化匹配為各星分配最合適的位置，停泊軌道高度也隨之確定。步驟如下：

1)衛星sk計算從初始緯度輻角運動至各目標緯度輻角所需走過的弧度值：

(29)

2)衛星sk對于每個目標緯度輻角，遍歷自身全部備選的停泊軌道，根據式(19)、(22)計算到達該緯度輻角的燃料成本和時間成本，并找出自身對于每個緯度輻角的最優停泊軌道及其部署成本。對固定的初始位置和目標位置，部署時間與燃耗呈單調性關系，對一個參數有明確期望時最優解很容易選擇。例如，某星初始高度為495 km，目標高度為498.5 km，計劃走過200°或160°到達目標位置。圖2顯示了較高或較低兩種停泊軌道(與目標軌道相比)對應的部署燃耗與部署時間關系。若要求部署任務在30天內完成，則可選擇圖中陰影區域內的解，但滿足時間要求的最佳解為三角形標注的解。

圖2 基于燃耗與部署時間的最優停泊軌道選擇

3)遍歷全部衛星，獲得每個衛星對于每個目標緯度輻角的最優停泊軌道高度及其最小燃耗。

2.2 基于KM算法的目標位置優化匹配

衛星需要從K個目標緯度輻角中選擇一個作為自身目標位置。這一過程實際上是各星帶著自身最小到達成本(燃耗)挑選目標位置的過程，使自身總燃耗盡量小的同時還需綜合考慮星座燃耗的優化與均衡，實現星座部署成本的整體優化。

2.2.1Kuhn-Munkres算法

Kuhn-Munkres算法由Kuhn[22]和Munkres[23]二人分別在1955年和1957年獨立提出，簡稱KM算法，被廣泛運用于帶權完全二分圖的最優匹配問題。若將衛星和目標位置看作一個無向圖的點集合，衛星與位置之間的匹配關系看作無向圖的邊，衛星到達各位置所需的最小燃耗看作邊的權值，這一問題即為典型的帶權值的二分圖最優匹配問題，可用KM算法獲得最優解。由于每個衛星必定對應一個目標位置，因此在約束條件設置合理時必定能夠找到該二分圖的完美匹配。二分圖是否存在完美匹配可由Hall匹配定理及其推論得出。若無法得到此二分圖的完美匹配，則表明約束條件設置不合理，應調整部署時間約束或軌道調整量約束。

設帶權完全二分圖G=(S,P,E,W),頂點集S={s1,s2,…,sk}代表衛星，頂點集P={p1,p2,…,pk}代表目標位置。邊集E代表衛星與位置的對應關系，若在規定時間和規定軌道調整量內，某星sk能夠達到某位置j，則邊skpj存在。邊skpj的權值用wkj=w(skpj)描述，表示衛星sk對于位置pj的最小達到成本(最小燃耗)。由于KM算法默認尋找最大權匹配，應先對權值wkj取倒數。KM算法的詳細研究可參考文獻[24]，本文僅簡要說明其過程。Kuhn-Munkres算法步驟如下：

1)為頂點集S、P取平凡標號l,確定l相等子圖Gl，并從Gl中選取任意匹配M。

2)由M出發，用匈牙利算法求出Gl的最大匹配M′。如M′是完美匹配，則M′是最優匹配，計算結束,否則轉步驟3)。

3)計算更新量αl更新標號l并返回步驟1)。

(30)

(31)

KM算法中用到匈牙利算法求解二分圖的最大匹配，匈牙利算法步驟如下：

1)取G的任意匹配M。

2)若M包括S中全部頂點，計算結束,否則轉步驟3)。

3)在S中尋找M匹配的非飽和點s，令Q={s},B=?。其中：符號Q代表圖G中任意頂點集合。符號B代表頂點集P中，在匹配M下與集合A={s|s∈S-u}中頂點配對的頂點集合，其中u為S中關于匹配M的非飽和點。

4)若點集Q的鄰集|NG(Q)|=|B|,計算結束。否則，在NG(Q)-B中任選一個頂點p。

5)如p為M匹配的飽和點，轉步驟6)。否則找到從s到p的M-可擴路R，作對稱差得M=MΔE(R),轉步驟2)。

6)找到與p匹配的點s，令Q=Q∪{s},B=B∪{y}，轉步驟4)。

2.2.2KM算法實現衛星位置匹配

圖3 KM算法實現3星優化匹配過程

1)從s1出發進行匹配，s1與p3頂點標號之和與邊權重相同，匹配成功(圖3(b))。

2)從s2出發進行匹配，只有s2與p3頂點標號之和與邊權重相同但s1與s2產生沖突，匹配失敗。

3)為參與沖突的頂點s1、s2、s3更新標號，左頂點標號減1，右頂點標號加1。可匹配邊變為s1-p3，s2-p3，s1-p1(新增)。即發生沖突時，有n個左頂點和n-1個右頂點參與，通過更新標號，整體能效下降(成本增加)了αl(n-(n-1))。

4)從s2出發進行匹配，s2-p1，匹配成功(圖3(c))。

5)從s3出發進行匹配，沒有符合與s3標號之和等于邊權重的頂點，匹配失敗。

6)為s3更新標號，并從s3出發進行匹配，發生沖突，匹配失敗(圖3(d))。

7)為參與沖突的頂點s1、s2、s3、p1、p3更新標號，可匹配邊變為s1-p1、s1-p3、s2-p1、s2-p2、s2-p3、s3-p3(圖3(e))。用匈牙利算法取目前相等子圖的最大匹配，該最大匹配已經是完美匹配，算法結束(圖3(f))。

可見，每發生一次沖突整體能效將降低一次，但由于每次降低的能效值都是最低值，整體仍是最優匹配。KM算法的核心思想是首先為頂點匹配最優邊，發生沖突時以降低整體能效為代價解決沖突，每次降低的能效均是最小值，再利用匈牙利算法完成最大匹配，最終獲得整體能效最高的匹配方案。

2.3 星座部署算法流程

星座部署算法可由地面或主星執行，解算最佳部署方式后回傳部署指令至各星。以主星集中控制為例，星座啟動部署后，主星收集從星軌道要素并確定參考軌道與參考相位。隨后主星依據參考軌道與當前從星位置信息確定各從星停泊軌道的選擇范圍，并解算各從星到達各目標位置所需的燃料成本。最后利用KM優化算法對各從星與目標位置進行優化匹配，同時確定停泊軌道高度并將信息回傳至從星。從星收到指令后執行軌道機動，達到停泊軌道后自由飛行一段時間并返回目標軌道。最后，主星收集當前各從星升交點赤經信息并以均值為目標升交點赤經，解算從星升交點赤經修正量并回傳，從星執行升交點赤經修正后，星座部署完成。優化后的星座部署算法流程見圖4。

圖4 星座優化部署算法流程圖

3 仿真校驗

初始時刻，10個衛星在軌道高度500 km±5 km，升交點赤經0°±0.5°，緯度輻角0°±10°的空間內隨機分布，軌道偏心率為0，傾角為30°。隨機生成各星初始軌道要素，見表1。其中以緯度輻角u代替近地點輻角和平近點角來表示圓軌道衛星相位。部署任務的約束如下：

表1 衛星初始軌道要素

1)構型約束：各星同平面且相位均勻分布。

2)時間約束：軌道高度部署和相位部署任務在20天內完成。

3)軌道調整量約束：各星原始軌道高度與停泊軌道高度差小于10 km。

3.1 傳統部署方法計算結果

傳統部署過程利用軌道高度差實現相位分離時，通常忽略衛星初始位置差異[19][25]。為公平對比，本文保留衛星高度差異。按逆時針順序對目標位置依次命名為1～10號，如圖5(a)所示，按構型約束相鄰目標位置相位相差36°。傳統部署過程如下：

1)根據軌道高度將衛星分為兩組，軌道較低的一半衛星分配1～5號位置，將執行降軌機動；軌道較高一半衛星分配6～10號位置，將執行抬軌機動，如圖5(b)所示。以軌道高度居中的兩星(S5、S6)的軌道高度均值作為目標軌道高度。

圖5 衛星與目標位置的未優化匹配結果

2)計算各星到達各自目標緯度輻角需要走過的相位角。根據時間約束和軌道高度調整量約束計算各星對應的停泊軌道高度，執行軌道機動并停泊。

3)全部衛星返回目標軌道后，用升交點赤經直接修正法將各星升交點赤經修正至當前升交點赤經中心位置。

表2顯示了采用傳統部署方法得到的各星停泊軌道高度、返回目標軌道后需要修正的升交點赤經余量和總燃耗。當部署時間約束(不包括升交點赤經修正)為20天時，采用傳統方法部署得到的各星平均燃耗為21.620 m/s，燃耗方差為367.709 m2/s2。

表2 衛星目標位置、停泊軌道高度與燃料預算-未優化

3.2 基于KM的部署優化算法優化結果

運用本文提出的部署優化算法得到各星與目標位置的匹配關系如圖6(b)所示。各星停泊軌道高度優化結果、返回目標軌道后需要修正的升交點赤經余以及總燃耗量見表3。當部署時間約束(不包括升交點赤經修正)為20天時，采用優化部署算法得到的各星平均燃耗為9.034 m/s，燃耗方差為5.594 m2/s2。對比傳統部署方法結果(平均燃耗21.620 m/s，方差367.709 m2/s2)可見，當構型約束、部署時間約束與軌道調整量約束相同的情況下，本文提出的優化算法顯著優于未經優化的傳統部署方案，優化后的結果在降低各星燃耗的同時，大大提高各星燃耗量的均衡性。

圖6 衛星與目標位置的優化匹配結果

3.3 KM的部署優化數值仿真

表3 衛星目標位置、停泊軌道高度優化結果與燃料預算

圖7 半長軸部署效果

圖8 緯度輻角部署效果

圖9 升交點赤經部署效果

表4 衛星燃料預算仿真結果

圖10 S1星推力施加情況

3.4 部署優化算法適用范圍討論

顯然，依靠停泊階段J2攝動能夠修正的升交點赤經量是有限的，且與停泊軌道高度、傾角、偏心率、停泊時間和緯度輻角調整量有關，即：

(32)

圖11 升交點赤經修正量與傾角、緯度輻角調整量的關系

當依靠停泊階段J2攝動差修正的升交點赤經量恰好等于升交點赤經全部待調量時，衛星返回目標軌道時的升交點赤經恰為目標升交點赤經，無需花費額外燃料用于升交點赤經修正，此時算法效益實現最大化。然而實際中難以遇到幾個參數完美匹配的情況，只能借此修正部分升交點赤經。用升交點赤經直接修正燃料ΔVΩ在衛星總燃耗ΔV中的占比定義優化算法效率：

(33)

即借助J2攝動修正的升交點赤經量越多，則升交點赤經余量越少，ΔVΩ在ΔV中占比越小，優化算法效率k越高。顯然k與軌道參數和緯度輻角調整量直接相關。500 km圓軌道，停泊時間為20天時，依據式(14)、(16)、(33)可繪制不同傾角下優化算法效率與初始升交點赤經分布范圍的關系，如圖12所示。可見該優化算法更適用于低傾角軌道以及緯度輻角調整量較大的情形。若以0.5為算法效率臨界值，則該算法針對初始升交點赤經在±0.7°范圍內分布的星群有較為突出的表現。具體場景下仍需根據軌道參數和緯度輻角調整量，參考式(32～33)評估適用性。

圖12 優化算法效率與傾角、初始升交點赤經分布的關系

4 結 論

針對異面至同面的低軌微納衛星星座部署問題，本文提出一種基于KM匹配的星座部署優化方法。利用J2攝動下的停泊階段相對漂移來補償初始位置差異，引入KM算法實現停泊軌道和目標位置的優化，在相位部署的同時盡量修正升交點赤經。采用N次有限常值推力實現軌道機動，適用于推力受限的微納衛星。在相同時間、目標構型和軌道高度調整量約束下，優化部署算法可將最小燃耗均值和方差由21.62 m/s、367.71 m2/s2降至9.03 m/s、5.59 m2/s2。最后對優化算法適用范圍做出討論。該方法彌補了傳統部署方法對衛星初始位置差異的忽略問題，算法簡單易實現，能夠快速給出優化結果，適用于在軌實時規劃，也可以擴展到星座變構或重構等應用場景中。