基于代理模型的衛(wèi)星編隊重構(gòu)最短距離建模*

李 鵬，陳琪鋒

(中南大學(xué) 航空航天學(xué)院·長沙·410083)

0 引 言

隨著微電子技術(shù)的快速發(fā)展以及微納小衛(wèi)星技術(shù)的不斷成熟，多顆小衛(wèi)星組成衛(wèi)星編隊共同完成復(fù)雜太空探索任務(wù)已成為當(dāng)今國際航天領(lǐng)域的一個研究熱點和發(fā)展方向。衛(wèi)星編隊在天基合成孔徑雷達、氣象衛(wèi)星三維立體成像和大型空間望遠鏡等領(lǐng)域有著廣泛的應(yīng)用前景[1-2]。

衛(wèi)星編隊構(gòu)形重構(gòu)是其區(qū)別于傳統(tǒng)單顆衛(wèi)星的重要特征和顯著優(yōu)勢之一。構(gòu)形重構(gòu)是指根據(jù)不同的任務(wù)需求，使編隊衛(wèi)星在規(guī)定的有限時間內(nèi)從初始構(gòu)形轉(zhuǎn)變?yōu)槠谕麡?gòu)形的機動過程，涉及軌位分配、最優(yōu)路徑規(guī)劃、星間碰撞避免、節(jié)省燃料和燃料均衡等問題[3]。編隊衛(wèi)星重構(gòu)機動的控制方式主要可以分為兩類，一是基于沖量控制，使用雙脈沖或多脈沖軌道轉(zhuǎn)移策略在多種約束下實現(xiàn)燃料消耗最優(yōu)。W.E.Wiesel[4]使用雙脈沖控制完成隊形重構(gòu)，通過優(yōu)化的方式尋找最佳的控制脈沖實施時刻。張玉錕[5]利用最優(yōu)二沖量控制研究編隊構(gòu)形重構(gòu)的底層控制問題。孟云鶴[6]提出了一種四沖量控制方法，實現(xiàn)了近地軌道的隊形重構(gòu)問題。王有亮等[7]采用遺傳算法和序列二次規(guī)劃結(jié)合的混合算法，研究了J2攝動下的衛(wèi)星編隊重構(gòu)的多脈沖軌跡優(yōu)化。二是基于連續(xù)推力控制，將構(gòu)形重構(gòu)軌跡優(yōu)化問題看作一個帶有狀態(tài)和控制約束的最優(yōu)控制問題，然后使用數(shù)值方法求解。M.Tillerson等[3]將隊形重構(gòu)問題簡化，使用線性規(guī)劃的方式優(yōu)化各衛(wèi)星的軌跡。G.T.Huntington等[8]使用高斯偽譜法，將衛(wèi)星編隊的隊形重構(gòu)問題轉(zhuǎn)化為非線性規(guī)劃問題，因其重構(gòu)過程在多個軌道周期內(nèi)完成，不需要考慮碰撞規(guī)避。Wu B.L.等[9]考慮了在近地軌道環(huán)境下應(yīng)用勒讓德偽譜法處理隊形重構(gòu)問題，能夠處理碰撞規(guī)避及J2攝動的影響，但其計算量較大。黃海濱[10]考慮了重構(gòu)過程中的碰撞避免約束，利用勒讓德偽譜法將衛(wèi)星重構(gòu)問題轉(zhuǎn)換為非線性規(guī)劃問題求解。

以上的研究或者在重構(gòu)過程中沒有考慮碰撞避免，或者考慮了碰撞避免約束，但算法只能處理編隊衛(wèi)星成員數(shù)量較少(如2～3顆衛(wèi)星)的情況。當(dāng)衛(wèi)星數(shù)量增多時，各衛(wèi)星重構(gòu)過程中發(fā)生碰撞的概率大大增加，此時必須考慮碰撞避免約束，需要計算編隊內(nèi)任意兩成員星在整個重構(gòu)過程中的最短距離，使其滿足安全距離約束。這會帶來極大的計算開銷，并增加優(yōu)化求解的困難。因此，有必要尋求一種能快速預(yù)測兩編隊衛(wèi)星重構(gòu)過程中最短距離的方法。

代理模型，又稱為元模型或者響應(yīng)面模型，可以作為解決上述問題的一種途徑。所謂代理模型，是指在原模型過于復(fù)雜或者不清楚輸入和輸出之間的具體數(shù)學(xué)關(guān)系(黑箱函數(shù))時，所建立的一種近似數(shù)學(xué)模型，屬于機器學(xué)習(xí)領(lǐng)域中的回歸類監(jiān)督學(xué)習(xí)[11-13]。它實際上是一種“模型的模型”，無需了解設(shè)計變量與響應(yīng)之間的物理意義，而是通過試驗設(shè)計得到樣本點的響應(yīng)信息，接著利用這些信息通過插值或者擬合的方法建立數(shù)學(xué)模型，從而對未知點的響應(yīng)值進行預(yù)測。代理模型經(jīng)過多年的發(fā)展，已被廣泛應(yīng)用于機械、流體、結(jié)構(gòu)和航空航天等領(lǐng)域。常用的代理模型技術(shù)包括多項式回歸響應(yīng)面(Polynomial Regression Surface，PRS)[14]、克里金(Kriging，KRG)模型[15]、徑向基函數(shù)(Radial Basis Function，RBF)[16]和人工神經(jīng)網(wǎng)絡(luò)(Artificial Neural Network，ANN)[17]等。

本文以近地空間中繞參考星為中心飛行的兩繞飛星為例，基于CW方程狀態(tài)轉(zhuǎn)移矩陣和雙脈沖機動變軌方式，使2顆繞飛星在一個軌道周期內(nèi)從自由繞飛狀態(tài)機動到空間圓編隊狀態(tài)。以編隊成員衛(wèi)星任意的初末狀態(tài)和首末脈沖施加時刻為輸入，重構(gòu)過程間的最短距離為輸出，基于三種不同大小的訓(xùn)練集，分別構(gòu)建了PRS、KRG、RBF、ANN代理模型，從預(yù)測精度和建模效率兩方面進行對比，從而選出最佳的編隊重構(gòu)星間最短距離代理模型。

1 基于CW方程的重構(gòu)最短距離模型

本文考慮了參考星軌道為圓軌道、繞飛星和參考星距離較近、重構(gòu)時間較短且不考慮攝動的情形，因此滿足使用CW方程的3個基本假設(shè)條件(圓軌道、近距離、短時間)[18]。

假設(shè)衛(wèi)星初始時刻t0的狀態(tài)為[r(t0)，v(t0)]T，衛(wèi)星在t時刻的狀態(tài)為[r(t)，v(t)]T，繞飛星相對參考星的相對運動可用CW狀態(tài)轉(zhuǎn)移方程來描述

(1)

根據(jù)CW方程的解析解，繞飛星以參考星為中心自由環(huán)繞飛行的初始條件為[5]

(2)

式中，x0、y0分別為初始時刻t0繞飛星相對參考星的位置矢量在參考星軌道坐標(biāo)系中的徑向和橫向分量；n為參考星環(huán)繞地球運動的軌道角速度。

在式(2)條件下，CW方程的解析解可簡化為三角函數(shù)形式

(3)

由式(3)可知，繞飛星環(huán)繞參考星的運動可分解為軌道平面內(nèi)和軌道面法向的周期性運動，軌道面內(nèi)的運動可由振幅b和相位角φ這2個參數(shù)決定，軌道面法向的運動則由振幅c和相位角ψ這2個參數(shù)決定。因此，可以用b、φ、c、ψ這4個參數(shù)來描述繞飛星相對參考星自由繞飛的運動狀態(tài)。

繞飛星形成空間圓編隊的初始條件為[5]

(4)

其中，r和θ分別表示空間圓編隊的半徑和相位角。

本文中編隊衛(wèi)星重構(gòu)的過程可描述為，初始時刻伴隨衛(wèi)星以參考星為中心自由環(huán)繞飛行，重構(gòu)開始后，通過在任意2個時刻施加速度沖量，使其在一個軌道周期后形成圍繞參考星運動的空間圓編隊。因此，要想得到重構(gòu)過程繞飛星的運動軌跡，就需求出施加的2次沖量大小。

雙脈沖變軌狀態(tài)方程的沖量響應(yīng)式為

(5)

式中，u=[03×1Δvi]T，Δvi表示每次機動施加的速度脈沖矢量；t1和t2分別表示施加2次脈沖的時刻。

由式(5)可知，若已知編隊成員衛(wèi)星的初始狀態(tài)和期望狀態(tài)以及2次首末脈沖的時刻，則可求出雙脈沖機動重構(gòu)所需的2次速度沖量為

(6)

得到2次速度沖量后，將其代入式(5)就可以求出繞飛星在整個重構(gòu)過程的運動狀態(tài)和軌跡。圖1所示為2顆繞飛星從初始自由繞飛狀態(tài)，經(jīng)過2次速度沖量作用后，變軌到空間圓編隊狀態(tài)時，其在軌道平面和三維空間的雙脈沖機動軌跡。圖中圓圈表示施加速度沖量的時刻。

(a)軌道平面

由式(3)可知，2顆繞飛星的初始狀態(tài)由b、c、φ、ψ這4個參數(shù)共同決定，期望狀態(tài)由空間圓編隊的半徑r和相位角θ決定。兩繞飛星的重構(gòu)軌跡則由施加沖量時刻t11、t12、t21、t22共同決定，其中t11、t12為繞飛星1施加2次脈沖的時刻，t21、t22為繞飛星2施加2次脈沖的時刻。另外需要注意的是，本文所建立的最短距離代理模型并不針對給定初末狀態(tài)的2顆編隊成員衛(wèi)星重構(gòu)過程，而是建立在滿足編隊約束時任意2顆編隊成員衛(wèi)星重構(gòu)過程之上的，即2顆成員衛(wèi)星初始自由繞飛時的狀態(tài)參數(shù)是任意的，期望空間圓編隊中的狀態(tài)參數(shù)也是完全任意的。因此，滿足從自由繞飛狀態(tài)重構(gòu)到空間圓編隊約束的任意2顆編隊衛(wèi)星重構(gòu)過程中的最短距離dmin可表示為

(7)

式中，T表示參考星的軌道周期；下標(biāo)1和2分別表示衛(wèi)星編隊的2顆繞飛衛(wèi)星；bmin、bmax分別表示編隊衛(wèi)星繞飛時在軌道平面內(nèi)所允許的最小和最大振幅；cmin、cmax分別表示編隊衛(wèi)星繞飛時在軌道面法向所允許的最小和最大振幅；rmin、rmax分別表示期望空間圓編隊的最小和最大半徑。

這個黑箱函數(shù)有16維輸入，1維輸出。在設(shè)計空間進行拉丁超立方抽樣得到訓(xùn)練集后，經(jīng)過訓(xùn)練即可得到重構(gòu)過程最短距離的代理模型。

2 代理模型簡介和模型性能指標(biāo)

2.1 PRS

PRS是一種基于最小二乘法的回歸分析方法，其中二次多項式響應(yīng)面模型在工程和實際應(yīng)用中最為常見和廣泛。本文中的多項式響應(yīng)面也選用二次多項式，假設(shè)輸入變量與預(yù)測響應(yīng)值具有如下關(guān)系[11]

(8)

式中，β為待估計的參數(shù)，可通過已知樣本點信息使用最小二乘法估計得到

(9)

2.2 RBF

RBF是一種常用的多元散點插值預(yù)測算法。該方法通過核函數(shù)將復(fù)雜高維問題轉(zhuǎn)化為簡單線性問題，RBF的輸入變量和預(yù)測響應(yīng)具有如下關(guān)系[16]

(10)

2.3 克里金法(KRG)

KRG模型是一種基于隨機過程的估計方差最小的無偏估計模型，也是一種插值模型。KRG模型區(qū)別于其他模型的最大特點是，它不僅能提供未知函數(shù)預(yù)測點的響應(yīng)，還能給出預(yù)估值的方差估計。KRG模型是一種插值模型，其插值結(jié)果為已知樣本點函數(shù)值的線性加權(quán)。為求解加權(quán)系數(shù)，KRG模型引入統(tǒng)計學(xué)假設(shè)，將未知函數(shù)看成是某個高斯靜態(tài)隨機過程的具體實現(xiàn)，即未知函數(shù)在設(shè)計空間中每一點處的響應(yīng)都是隨機變量，服從均值為fTβ、方差均為σ2的正態(tài)分布。KRG模型的輸入變量和預(yù)測響應(yīng)具有如下關(guān)系[15]

(11)

R=(R(xi,xj))ij∈Rm×m
r=[R(x1,x),R(x2,x),…,R(xm,x)]

(12)

式中，R(xi，xj)稱為2個樣本點的相關(guān)函數(shù)，是2個樣本點歐氏距離的函數(shù)，代表不同位置處隨機變量的相關(guān)性，且R的值隨著距離的增加而減小，即距離為零時R=1；距離無窮大時R=0。KRG模型的性能由其基函數(shù)和相關(guān)函數(shù)共同確定。常見的基函數(shù)有常數(shù)項、一次和二次多項式。在未知函數(shù)的整體趨勢未知的情況下，常數(shù)項通常被認(rèn)為是一個合適的基函數(shù)選擇，因此本文選用常數(shù)項作為KRG模型的基函數(shù)，此時亦稱為普通KRG模型。常用的相關(guān)函數(shù)模型有高斯函數(shù)、各項同性高斯指數(shù)函數(shù)、各項異性高斯指數(shù)函數(shù)、三次樣條函數(shù)和Matern函數(shù)。本文選用綜合性能相對穩(wěn)定的Matern函數(shù)作為KRG模型的相關(guān)函數(shù)。

2.4 ANN

ANN作為機器學(xué)習(xí)領(lǐng)域中最為常見的算法，由于其對非線性函數(shù)具有很強的擬合能力而得到了廣泛的應(yīng)用。神經(jīng)網(wǎng)絡(luò)最基本的單位是神經(jīng)元，將多個神經(jīng)元按照一定的規(guī)則組織在一起，就構(gòu)成了最基本的ANN。最常用的神經(jīng)網(wǎng)絡(luò)為BP神經(jīng)網(wǎng)絡(luò)。BP神經(jīng)網(wǎng)絡(luò)的結(jié)構(gòu)主要包括三層：輸入層、隱含層和輸出層[17]。有關(guān)BP神經(jīng)網(wǎng)絡(luò)的具體原理本文不再贅述，本文所采用的BP神經(jīng)網(wǎng)絡(luò)來自MATLAB中的神經(jīng)網(wǎng)絡(luò)擬合工具箱，隱含層的神經(jīng)元數(shù)量設(shè)置為100。

2.5 模型性能評價指標(biāo)

本文分別從模型精度和時間效率兩方面來評價各代理模型的優(yōu)劣。其中模型精度主要關(guān)注代理模型的全局誤差準(zhǔn)則，因此選擇量綱為1的決定系數(shù)R2、相對平均誤差ERA和有量綱的均方根誤差ERMS[11]。R2和ERA的值均在0～1之間變化，且R2的值越高表示代理模型性能越好，其預(yù)測精度越高；相反，ERA和ERMS的值越低表示模型性能越好，其預(yù)測精度越高。

(13)

(14)

(15)

3 仿真算例

3.1 仿真條件

參考星初始時刻軌道參數(shù)：軌道高度h=622km，偏心率e=0，軌道傾角i=45°，升交點赤經(jīng)Ω=30°，近地點角距ω=45°，真近點角f=60°。繞飛星的初始狀態(tài)的振幅和期望狀態(tài)的空間圓半徑范圍為：bmin=0.2km，bmax=1km，cmin=0.2km，cmax=1km，rmin=0.5km，rmax=1km。

根據(jù)式(7)，在設(shè)計空間用拉丁超立方設(shè)計方法分別抽樣5000、10000、20000個樣本點，組成3個不同大小的訓(xùn)練集，分別構(gòu)建PRS、KRG、RBF、ANN代理模型。接著采用蒙特卡羅法抽樣得到2000個測試點組成測試集，用于測試和驗證各代理模型的精度。最后分別對比各代理模型的精度和訓(xùn)練集大小對代理模型性能的影響。

3.2 結(jié)果和對比

3.2.1 精度

(1)訓(xùn)練集大小為5000

圖2給出了訓(xùn)練集大小為5000時各代理模型的精度對比，如前所述，R2越接近1，表示模型精度越高；ERA和ERMS越小，表示模型精度越高。圖2(a)顯示，在訓(xùn)練集大小為5000時KRG模型精度最高，R2為0.6760，ERA為0.3178；RBF模型的精度稍遜于KRG模型；ANN模型的精度最差，R2僅有0.4117。原因可能是本問題的維數(shù)較高，神經(jīng)網(wǎng)絡(luò)訓(xùn)練得到的模型若想取得比較理想的性能，通常需要大量的訓(xùn)練點。此外，盡管KRG模型的精度最高，但是其模型性能仍然比較差，圖2(b)顯示，其ERMS為0.0968km，接近100m，難以滿足實際需要。因此，需要增加訓(xùn)練集的樣本點，以期望改善模型性能。

(2)訓(xùn)練集大小為10000

圖3給出了訓(xùn)練集大小為10000時各代理模型的精度對比，圖3(a)顯示，此時KRG模型的精度仍然最高，R2為0.7802，ERA為0.2159；RBF模型的精度仍然次于KRG模型；但此時ANN的精度有所提升，R2為0.6313，ERA為0.2827，接近于RBF模型。此外，相較于訓(xùn)練集大小為5000時，各代理模型的精度都得到了一定程度的提升。圖3(b)顯示，表現(xiàn)最好的KRG模型的ERMS為0.0735km，比訓(xùn)練集大小為5000時提升了24%左右，其他模型的均方根誤差也都有所下降。

(a)決定系數(shù)R2和平均相對誤差ERA

(a)決定系數(shù)R2和平均相對誤差ERA

(b)均方根誤差ERMS圖3 訓(xùn)練集大小為10000時各代理模型的精度Fig.3 Accuracy of each surrogate model when the training set size is 10000

(3)訓(xùn)練集大小為20000

圖4給出了訓(xùn)練集大小為20000時各代理模型的精度對比，圖4顯示，當(dāng)訓(xùn)練集大小增加到20000時，PRS模型的精度并沒有明顯提升，其他三種模型的精度都得到了提升。其中KRG模型的精度最高，R2達到了0.8240，ERA為0.1775，降到了0.2以下，ERMS為0.0639km，比較接近60m。ANN模型的精度此時超過了RBF模型，由此可見，ANN模型對訓(xùn)練點的大小最為敏感，隨著訓(xùn)練點的增加，ANN模型精度得到了較大的提升。

(a)決定系數(shù)R2和平均相對誤差ERA

3.2.2 效率

代理模型的效率即為模型訓(xùn)練和預(yù)測未知點響應(yīng)所花費的時間，時間越短，表示其效率越高。效率同時還取決于訓(xùn)練集的樣本大小和所使用的計算平臺，本文以10000個樣本點的訓(xùn)練集和2000個測試點組成的測試集為例，統(tǒng)計了各代理模型訓(xùn)練和預(yù)測未知點消耗的時間。在一臺CPU主頻為2.80GHz的電腦上使用MATLAB軟件進行所有的計算。表1給出了各模型訓(xùn)練以及預(yù)測2000個測試點所需的大概時間。

表1 訓(xùn)練集大小為10000時各代理模型的效率

由表1可知，PRS模型的效率是最高的，模型訓(xùn)練和預(yù)測的時間都少于1s。同時發(fā)現(xiàn)，KRG模型的訓(xùn)練是非常耗時的，需要花費5～6h。這主要有兩方面的原因，一是因為在構(gòu)建KRG代理模型時，為了獲得最佳的模型精度，需要對相關(guān)函數(shù)模型中的參數(shù)(或超參數(shù))進行優(yōu)化，使得其似然函數(shù)達到最大值，本文研究的問題有16個輸入變量，因此該優(yōu)化問題為一個16維的優(yōu)化問題；二是因為在構(gòu)建KRG模型時需要求相關(guān)矩陣Rm×m的逆，m為樣本點的個數(shù)，本例中為10000，因此當(dāng)問題維數(shù)和樣本量很大時，訓(xùn)練KRG模型所花費的時間較長。而對于模型的預(yù)測時間，同樣發(fā)現(xiàn)KRG模型耗時最長，預(yù)測2000個測試點所需的時間為5～10s，這主要是因為KRG模型預(yù)測時需要計算未知點和各樣本點組成的相關(guān)向量。然而需要指出的是，當(dāng)時間步長為1s時，計算2000個測試點的真實輸出響應(yīng)耗時197s，與之相比，KRG模型預(yù)測所花費的時間仍然非常短，大大提高了計算效率。

4 結(jié) 論

本文考慮編隊衛(wèi)星成員數(shù)量較多時，為了生成燃料最優(yōu)的無碰撞重構(gòu)軌跡，需要計算任意2顆衛(wèi)星重構(gòu)過程的最短距離。為減少優(yōu)化過程中的計算開銷，本文以近地空間中繞參考星飛行的2顆衛(wèi)星從初始自由繞飛狀態(tài)重構(gòu)機動到空間圓編隊狀態(tài)為例，基于CW方程和雙脈沖機動變軌策略，分別構(gòu)建了PRS、KRG、RBF、ANN模型，并從預(yù)測精度和建模效率兩方面進行了對比。仿真結(jié)果表明，PRS模型的效率最高，但其精度較差，且隨著訓(xùn)練集的增大，其精度并無明顯提升。KRG模型的精度最高。隨著訓(xùn)練集的增加，KRG和ANN模型的精度都得到了明顯改善。KRG模型的缺點是由于問題維數(shù)較大和樣本容量大導(dǎo)致模型訓(xùn)練耗時較長，且KRG模型預(yù)測未知點響應(yīng)的時間同樣耗時最長，但與真實模型相比，其計算耗時仍然很短，因此可以使用KRG模型作為衛(wèi)星編隊重構(gòu)過程中星間最短距離建模的代理模型，用于后續(xù)的優(yōu)化求解，以生成燃料最優(yōu)的無碰撞重構(gòu)規(guī)劃軌跡。