基于目標分配策略的空間碎片清除方案設計

楊保臻, 錢霙婧

北京工業(yè)大學, 北京 100124

0 引 言

近年來,人類探索宇宙奧秘的航天活動日益頻繁,航天技術飛速發(fā)展的同時也帶來新問題,空間碎片的個數不斷增長,對在軌運行的航天器造成巨大的威脅,發(fā)生碰撞的風險大大增加.空間碎片在軌道上與航天器的平均相對運行速度大約10 km/s,在如此高的速度下毫米級的碎片撞擊都有可能將在軌航天器的外殼擊穿.1978年,美國宇航局科學家KESSLER[1]首次提出一項理論假設,認為當在近地軌道運轉的航天器等的密度達到一定程度時,航天器相互碰撞后產生的碎片能夠形成更多的新撞擊,形成級聯效應,此效應稱為Kessler效應,意味著即使從此不再發(fā)射新的航天器進入太空,近地軌道同樣將覆蓋危險的太空垃圾.由于失去能夠安全運行的軌道,在之后的數百年內太空探索和人造衛(wèi)星的運用將變得無法實施.因此,主動清除空間碎片是解除空間碎片危機的根本方法.

空間碎片的主動清除目前尚處于研究階段[2],設想大致可分為接觸式[3](通過捕獲機構進行抓捕碎片)和非接觸式[4](強激光束照射碎片降低碎片軌道、靜電力增阻對碎片充電使其減速離軌).其中,非接觸式主動清除方法雖然可以避免與空間碎片產生直接接觸,但是同時也增加了相應的難度和不可控制性.接觸式捕獲離軌的清除方式不但可以用來執(zhí)行其他持續(xù)清除碎片任務,同時還可以適用較廣的軌道高度范圍以及碎片清除范圍[5-6].

在執(zhí)行空間碎片清除任務時,首先應考慮數量眾多的空間碎片的清除對象選擇.作為非合作目標的空間碎片,通過掌握其分布的總體情況與碎片的物理特性確定清除對象的大致選擇范圍[7-8].對于一對多的碎片除任務包括清除碎片數量、清除序列的設定等問題,于錫錚等[9]研究了單次任務中在多脈沖推力作用下軌道機動清除12顆集中分布的碎片的方法.BONNAL等[10]研究了單個飛行器清除多碎片的問題,并利用權重系數法求解了同時考慮燃耗和時間的單目標優(yōu)化問題.BéREND和CERF等[11-12]利用分支定界法研究了多碎片清除任務序列優(yōu)化和軌道轉移優(yōu)化問題.MISSEL等[13]采用遺傳算法求解系統中的碎片清除任務規(guī)劃問題.先前的文獻研究從主動清除角度考慮了單個飛行器一對多持續(xù)清除碎片任務,且需清除的多塊碎片有一定的位置分布要求,無法短時間內協同清除大量威脅在軌衛(wèi)星或者飛行器的多塊空間碎片.

不同于以往使用單個飛行器執(zhí)行任務,本文以利用衛(wèi)星星座構型中多顆在軌衛(wèi)星同時釋放帶有抓捕機構的飛行器執(zhí)行碎片清除任務為背景,針對碎片對象選擇及清除序列采用單轉移飛行器進行一對多碎片清除任務.以燃料消耗為指標,利用降低碎片分布密度的DBSACN聚類算法對初步選取的空間碎片群進行碎片分類,并采用基于匈牙利算法的任務分配選擇具體在軌衛(wèi)星下放飛行器執(zhí)行相應碎片持續(xù)清除任務.相對于以往研究可實現快速降低相應區(qū)域碎片分布密度,保障在軌航天器的安全.

1 問題描述

1.1 空間飛行器與衛(wèi)星運動模型

慣性坐標系是指相對于宇宙的其他部分而言沒有加速度和轉動的坐標系.如圖1所示,其原點位于地球質心Oe處,OeX軸在赤道平面內沿著地心Oe與平春分點的連線,指向平春分點;OeZ軸沿地球自轉軸方向從地心指向北極點.OeY軸與OeX、OeZ軸垂直并滿足右手定則.

圖1 地心慣性坐標系(ECI)Fig.1 Geocentric inertial coordinate system (ECI)

衛(wèi)星執(zhí)行碎片清除任務時釋放的空間飛行器在飛行過程中受到地球的引力,二體動力學方程[14]在地心慣性坐標系下表示為

(1)

1.2 空間飛行器交會制導方法

空間飛行器交會對接過程一般分為4個階段:遠距離導引段、近距離導引段、平移靠攏段和對接段[15-16].本文重點在于多碎片清除任務設計,制導方法采用Lambert制導來估算遠距離導引段算法中衛(wèi)星釋放空間飛行器捕捉碎片以及空間飛行器持續(xù)清除碎片所需的燃料消耗.

空間飛行器在大氣層外飛行受到地球引力作用與其構成二體模型.二體運動相對慣性空間的運動方程是二階常微分方程[17].具體到當前的二體問題,若初始位置矢量和速度矢量已知,可以通過二體運動方程計算空間飛行器相應時間的位置和速度.此問題稱為初值問題.若初始時刻和終點的位置矢量已知,則求得空間飛行器初始位置處的速度矢量.此問題稱為兩點邊值問題.

Lambert問題是在固定時間約束下的兩點邊值問題,是航天器軌道動力學的經典和基本問題.通過求解Lambert問題計算衛(wèi)星釋放的帶有捕獲機構的飛行器飛行到目標碎片附近所需的速度增量,如圖2所示,僅考慮中心天體C的引力,飛行器Lambert轉移軌道為從起始位置P1點飛行到終點位置P2點,因此相對于中心天體C,確定起始位置矢量與終點位置矢量,飛行時間為Δt.起點P1的地心距為r1,終點P2的地心距為r2,飛行路徑所對應的地心角為θ.此兩點軌道轉移所需要的時間Δt僅與始末位置的幾何構型以及軌道半長徑有關,而與其他軌道參數無關,通過給出初始時刻與終點時刻的位置矢量r1、r2、轉移軌道的角度θ、轉移時間Δt并通過Lambert求解算法則可得到空間飛行器變軌機動后初始時刻與終點時刻的速度矢量[18].

圖2 Lambert轉移軌道幾何構型Fig.2 Lambert transfer orbital geometry

2 目標碎片選擇與分類

2.1 目標碎片初步選擇

地球周圍的空間碎片不計其數,按照尺寸大小可分為3大類:(1)大型空間碎片:直徑超過10 cm,占據碎片總質量的99%;(2)中型碎片:直徑在1～10 cm之間,也成為了危險碎片;(3)小型碎片:直徑小于1 cm,數量巨大且難以探測[19-20].空間碎片的分布密集區(qū)域主要集中在3部分:500～2 000 km的LEO區(qū)域、20 000 km的中軌區(qū)域以及36 000 km的GEO區(qū)域.在空間碎片分布的3大密集區(qū)域中,LEO相對運動速度較大并且運動方向存在交叉,碰撞事件的發(fā)生概率較高,因此本文選定LEO區(qū)域的大型碎片與火箭箭體為初步清除目標.

根據Space-Track網站公布的數據,如圖3～5所示,LEO區(qū)域中在軌的中小型碎片分別有11 326個、1 643個,在軌大型碎片有116個,在軌火箭箭體934個.中小型碎片的絕大部分是由于大型碎片或者火箭箭體碰撞導致,為抑制Kessler效應,應優(yōu)先考慮清除箭體等大型碎片,并且箭體等結構堅固,對捕獲機構要求不高.根據圖5中空間碎片分布,大型空間碎片和火箭箭體在軌道傾角97°～100°間分布較為密集.考慮到減少燃料消耗,除減少較大傾角的軌道機動外,還應選擇升交點赤經變動在一定范圍內的碎片.

圖3 LEO在軌中小型碎片Fig.3 Small and medium fragments in LEO

圖4 LEO在軌大型碎片Fig.4 Large debris in LEO

圖5 LEO在軌火箭箭體Fig.5 Rocket bodies in LEO

2.2 基于密度的聚類算法-DBSCAN算法

通過確定部分軌道要素與碎片大小的方法已經初步選定一組目標碎片,希望通過空間飛行器實現一對多的碎片清除任務要根據碎片特性進一步分類.聚類是無監(jiān)督學習任務之一,聚類是數據挖掘中的概念[21],就是按照某個特定標準或特征(如距離)把一個數據集分割成不同的類或簇,需要通過分析將數據對象中潛在的特征信息發(fā)掘并設定合適的評判標準,使得同一個簇內的數據對象的相似性盡可能大,同時不在同一個簇中的數據對象的差異性也盡可能大.聚類后同一類的數據盡可能聚集到一起,不同類數據盡量分離.DBSCAN算法[22]的聚類定義很簡單:由密度可達關系導出的最大密度相連的樣本集合,即為本文最終聚類的一個類別,或者說一個簇.

DBSCAN是通過一組鄰域來描述樣本集的緊密程度,參數(∈,Minpts)用來描述鄰域的樣本分布緊密程度.其中,用∈描述某一樣本的鄰域距離閾值,Minpts描述某一樣本的距離為∈的鄰域中樣本個數的閾值,樣本集記為D=(x1,x2,…,xm),DBSCAN算法的關鍵性定義如下:

(1)∈-鄰域:對于樣本集中任一樣本xj∈D,其∈-鄰域包含樣本集D中與xj的距離不大于∈的子樣本集,即N∈(xj)∩{xi∈D|distance(xi,xj)≤∈},這個子樣本集的個數記為|N∈(xj)|.

(2)核心對象:對于任一樣本xj∈D,如果其∈-鄰域對應的N∈(xj)至少包含Minpts個樣本,即如果|N∈(xj)|≥Minpts,則xj是核心對象.

(3)密度直達:若xi位于xj的∈-鄰域中,且xi是核心對象,則稱xi由xj密度直達.但是反之不一定正確,即此時并不能確定xj由xi密度直達, 除非且xi也是核心對象.

(4)密度可達:對于xi和xj,如果存在樣本序列p1,p2,…,pT,滿足p1=xi,pT=xj,且pt+1由pt密度直達,則稱xj由xi密度可達.即密度可達滿足傳遞性.

通過DBSCAN算法劃分的簇中有一個或者多個核心對象.當某個簇里只有一個核心對象時,則簇里其他的非核心對象樣本都在此核心對象的∈-鄰域里;當某個簇里有多個核心對象,則簇里任意一個核心對象的∈-鄰域中必定有一個其他的核心對象從而實現兩個核心對象間的密度可達.這些核心對象的∈-鄰域里所有的樣本的集合組成一個聚類簇.

根據初步選擇出的一組目標碎片,將其所在位置相互間轉移所需的速度增量集合設定為一樣本集,∈-鄰域描述為空間飛行器從某碎片到另一碎片所需的燃料消耗閾值,Minpts描述為到達某一碎片后,周圍存在燃料消耗允許范圍內潛在的清除碎片數量閾值.

如圖6所示,Minpts=5,紅色的點都是核心對象,因為其∈-鄰域至少有5個樣本.黑色的樣本是非核心對象.所有核心對象密度直達的樣本在以紅色核心對象為中心的超球體內,若不在超球體內,則不能密度直達.用綠色箭頭連起來的核心對象組成了密度可達的樣本序列.在上述密度可達的樣本序列的∈-鄰域內所有樣本構成一個簇.

圖6 DBSACN分類概念說明圖Fig.6 DBSACN classification concept diagram

針對碎片清除任務,碎片分類算法設計如下:

算法1:碎片分類算法輸入:樣本集D,鄰域參數∈,樣本距離度量Minpts輸出:碎片簇向量C與噪聲Is與核心對象Vcore1 由D生成n×n成本矩陣M,n維全零向量C、全邏輯零向量Vs、Is, Q=02 for i =1 to n do3 if Vs(i)為0 then4 將Vs(i)改為邏輯1,尋找M中當前行滿足∈要求的對象Neighbors5 if Neighbors數量

輸入:樣本集D=(x1,x2,…,xm),鄰域參數(∈,Minpts),樣本距離度量方式(碎片間轉移所需燃料消耗)

輸出:碎片樣本集劃分出的簇C={C1,C2,…,Ck}與噪聲

1)初始化核心對象集合Ωcur=φ,初始化聚類簇數,初始化未訪問碎片樣本集合Γ=D,簇劃分C=φ;

2)對于j=1,2,…,m,按下述步驟找出所有核心對象:

a)通過距離度量方式,找到樣本xj的∈-鄰域子樣本集N∈(xj);

b)若子樣本集中樣本個數滿足|N∈(xj)|≥Minpts,將樣本xj加入核心對象樣本集合:Ω=Ω∪{xj}.

3)若核心對象集合Ω=φ,跳出步驟2,計算完成,否則轉入步驟4;

4)在核心對象集合中,任選一核心對象o,初始化當前簇核心對象序列Ωcur={o}, 初始化類別序號k=k+1,初始化當前簇樣本集合Ck={o},更新未訪問樣本集合Γ=Γ-{o};

5)如果當前簇核心對象序列Ωcur=φ,則當前聚類簇Ck生成完畢,更新簇劃分C={C1,C2,…,Ck}更新核心對象集合Ω=Ω-Ck,轉入步驟3.否則更新核心對象集合Ω=Ω-Ck;

6)在當前簇核心對象序列Ωcur中選取一核心對象o′,通過鄰域距離閾值∈找出所有的∈鄰域子樣本集N∈(o’),令Δ=N∈(o’)∩Γ,更新當前簇樣本集合Ck=Ck∪Δ,更新未訪問樣本集合Γ=Γ-Δ,更新Ωcur=Ωcur∪(Δ∩Ω)-o’,返回步驟5.

3 基于匈牙利算法的清除碎片任務分配策略

3.1 求解指派問題的匈牙利算法

KUHN提出的匈牙利算法(Hungarian algorithm)是一種關于指派問題的求解方法,其引用了匈牙利數學家康尼格的一個關于矩陣中獨立零元素個數的定理:矩陣中獨立零元素的個數等于能夠覆蓋所有零元素的最少直線數[23-24].結合本文研究,矩陣中的不同行和列確定的位置元素對應不同空間飛行器去抓捕不同碎片所需的燃料消耗.算法基本思想是修改效率矩陣的行或列,使得每行或每列中至少有個零元素,通過修改直至在不同行列中最少存在一個零元素,得到與所有零元素位置相對應的任務分配方案并且為效率矩陣中的最優(yōu)分配,此方案使任務花費的成本最小.

假設一顆衛(wèi)星可以釋放n個空間飛行器,將所有衛(wèi)星能夠釋放的空間飛行器編入執(zhí)行任務集,則衛(wèi)星釋放空間飛行器清除碎片指派問題模型.

(1)決策變量

若xij=1,即指派第i個空間飛行器到達第j碎片簇的第一個核心對象附近并持續(xù)清除碎片簇中相應核心對象.若xij=0,即不指派第i個空間飛行器到達第j碎片簇的第一個核心對象附近并持續(xù)清除碎片簇中相應核心對象.

(2)目標函數

(2)

式中,Cij為第i個空間飛行器到達第j碎片簇的第一個核心對象附近所需的燃料消耗.

(3)約束條件

a)xij只能為0或1;

求解指派任務步驟流程如圖7所示,其中關鍵步驟定義如下:

圖7 求解指派問題的匈牙利算法流程Fig.7 Hungarian algorithm flow for solving assignment problem

(1)行列歸約.尋找每行和每列中最小元素,分別從每行和每列中減去這個最小元素.

(2)指派任務(確定獨立零元素).圈零法:依次尋找只有一個零元素的行或列并圈出,并劃去該零元素所在的列或行中其他零元素,此時可能出現3種情況:

a)每行都有獨立零元素,個數滿足m=n時得到最優(yōu)解,跳出所有計算步驟,算法完成.

b)存在未被標記的零元素,并且其所在行列中未被標記的零元素均至少有兩個,可得到最優(yōu)解.此時從剩余零元素最少的行或列開始,選零元素畫圈,然后劃掉同行同列的其它零元素,反復進行,直到所有零元素均被圈出或劃掉為止.

c)不存在未被標記過的零元素,但圈零個數m

(3)畫蓋零線.利用最少的水平線和垂直線覆蓋所有的零元素:

a)對效率矩陣中所有不含圈零元素的行打√;

b)對打√的行中所有零元素所在列打√;

c)對所有打√的列中圈零元素所在行打√;

d)重復上述第b)和c)步,直到不能繼續(xù)為止;

e)對未打√的每一行畫一直線,對已打√的每一列畫一縱線.

(4)更新矩陣.跟經過畫蓋零線后的矩陣進一步交換增加零元素,在未被直線覆蓋過的元素中找出最小元素,將打√行的各元素減去這個最小元素,同時將打√列的各元素加上這個最小元素.

3.2 整體清除碎片任務分配

每顆衛(wèi)星可以釋放n個空間飛行器,將所有衛(wèi)星能夠釋放的空間飛行器編入執(zhí)行任務集與碎片簇位置速度數據作為輸入,整體任務分配流程如圖8所示.

圖8 整體任務分配實現流程Fig.8 Overall task allocation implementation process

飛行器挑選準則:通過空間飛行器任務集與計算空間飛行器到達碎片簇間所需速度增量生成速度增量集,依據速度增量集將所有空間飛行器到達相應碎片簇首個核心對象附近所需速度增量按照由小到大順序排列,后根據所需執(zhí)行任務的空間飛行器數量首先選擇到達相應碎片簇首個核心對象附近所需速度增量的飛行器,即速度增量集排列后的第一行數據對應的空間飛行器編號,若數量沒有達到執(zhí)行任務的飛行器數量要求,則繼續(xù)從下一行選取直到滿足數量要求.值得注意的是,流程中替換執(zhí)行任務的飛行器時從先前挑選終止處繼續(xù)進行.

4 仿真算例

考慮到清除LEO附近的空間碎片,算例采用基于Walker星座隨機生成LEO附近容易受此處碎片威脅的一星座構型中的在軌衛(wèi)星釋放飛行器執(zhí)行碎片清除任務,其中飛行器由抓捕手臂和平臺兩部分構成,平臺包括結構、通訊、制導等系統構成,星座衛(wèi)星具體軌道數據如表1所示,衛(wèi)星軌道半長軸為7 575.3 km,偏心率為0.042 86,軌道傾角為1.435 rad.

表1 衛(wèi)星星座軌道要素數據Tab.1 Orbit element data of satellite constellation

根據先前的空間碎片分析對在軌碎片進行初步選擇,并選擇升交點赤經差值在20°范圍的多組碎片,如圖9所示共計282顆.

圖9 目標碎片分布Fig.9 Target fragments distribution

如圖10所示,DBSCAN算法中鄰域參數(∈,Minpts)在任務背景下的具體表現為(速度增量限制、超球體內碎片密度),6組任務鄰域參數分別為(2 km/s,4塊)、(2 km/s,3塊)、(1.5 km/s,3塊)、(3 km/s,4塊)、(3 km/s,5塊)和(3 km/s,6塊),飛行器在碎片間轉移時間設定為100 s,以此進行分類,分別選出45塊、118塊、58塊、139塊、96塊和67塊碎片,分別分為14個、41個、22個、45個、30個和22個碎片簇.

圖10 目標碎片分類Fig.10 Target fragments classification

圖10(a)～(f)中藍色叉號為噪聲,彩色各點不同顏色代表不同的碎片簇,同一顏色為同一碎片簇.圖11(a)～(f)中為各個碎片簇中的核心對象數量情況.

圖11 各碎片簇中核心對象數量Fig.11 Number of core objects in each debris cluster

假設每顆衛(wèi)星最多能釋放3個空間飛行器,6次整體碎片清除任務如圖12(a)～(f)所示,黃色為執(zhí)行任務的衛(wèi)星釋放空間飛行器,紅色為碎片簇第一個核心對象.

圖12 整體碎片清除任務Fig.12 Overall fragments removal task

假設每顆衛(wèi)星最多能釋放3個空間飛行器,空間飛行器與碎片交會制導時間為100 s,對速度增量的限制在10 km/s,6次整體碎片清除任務中具體所需速度增量情況與任務分配情況如圖13(a)～(f)所示,圖中標記為碎片編號.

圖13 碎片清除任務分配Fig.13 Debris clearing task allocation

仿真結果可知,在碎片分類中,參數約束較小的B組與D組選取了更多碎片進行分類,分成了40個以上的碎片簇,并且有接近10個碎片簇中有多個核心對象;有參數約束較大的A組選取了較少碎片分類,分成了20個以下的碎片簇,并且絕大部分碎片簇中只包含一個核心對象.6次整體碎片清除任務星座構型的75顆衛(wèi)星分別采用了其中10顆、27顆、17顆、27顆、20顆和14顆衛(wèi)星執(zhí)行各類碎片簇的清除任務,其中包含有些衛(wèi)星被多次匹配,釋放1～3個空間飛行器執(zhí)行不同碎片簇的清除任務,并保證了交會制導所需速度增量大部分保持在6 km/s以下,DBSACN聚類算法隨著超球體中碎片密度的降低與所需最大燃料消耗要求的提高,更容易從目標碎片中分出多顆碎片執(zhí)行任務,并且碎片簇數量越多,分類越散,許多碎片簇中只包含一個核心對象,沒有進行一對多的碎片清除任務.

5 結 論

本文通過對空間碎片分布情況進行分析,初步選擇一組目標碎片,并利用碎片清除任務實際情況根據碎片間轉移所需燃料消耗及碎片分布密度對所選碎片進行不同種類的有效區(qū)分,劃分不同的碎片簇.根據求解指派問題的匈牙利算法通過星座中在軌衛(wèi)星設計整體碎片清除分配任務,仿真結果表明通過星座構型中10～30顆在軌衛(wèi)星可以短時間內在較小燃料消耗成本下完成任務.