趙建磊, 李海陽,*

(1. 國防科技大學(xué)空天科學(xué)學(xué)院, 湖南 長沙 410073;2. 空天任務(wù)智能規(guī)劃與仿真湖南省重點實驗室, 湖南 長沙 410073)

0 引 言

隨著人類探索地球外層空間活動的不斷深入,空間飛行器尤其是非合作空間飛行器的數(shù)量日益增長。出于對有限空間資源的爭奪,惡意抵近服務(wù)如干擾、偵照、抓捕等行為逐漸增加,飛行器在軌安全受到嚴(yán)重威脅。通常情況下,飛行器根據(jù)實際任務(wù)需要對自身進(jìn)行的軌道機動、軌道維持、抵近觀測、姿態(tài)控制等操作的信息是不公開的。在此背景下,非合作目標(biāo)的軌道機動信息的感知、檢測和識別對及時發(fā)現(xiàn)非合作飛行器的惡意抵近服務(wù)意圖,維護在軌飛行器的安全顯得尤為重要。

國內(nèi)外學(xué)者從不同角度對軌道機動檢測和機動信息識別等問題進(jìn)行研究。Kelecy等人對脈沖小推力作用下的定軌與機動檢測進(jìn)行研究,對比分析了最小二乘估計方法與卡爾曼濾波方法的適用性。Lemmens等人分別提出了基于兩行軌道要素(two line elements, TLE)一致性檢驗的低軌衛(wèi)星機動檢測方法和支持近實時機動檢測的TLE時間序列分析方法。Pittelkau和Acquesta等人分別提出了機動檢測中避免歷史數(shù)據(jù)中噪聲干擾的新方法,其中前者采用中值濾波器,而后者則采用被稱為Box具有噪聲的基于密度的聚類方法(density-based spatial clustering of applications with noise, DBSCAN)。Bai等人利用衛(wèi)星的歷史軌道信息采用聚類的方法進(jìn)行了機動檢測,分析了均值聚類、層次聚類和模糊聚類在機動檢測中的應(yīng)用。Yu等人通過分析機動前后的攝動,利用相對動力學(xué)方程推導(dǎo)出兩種機動模型,提出了基于軌道動力學(xué)模型的機動檢測方法。Daniel等人提出了一種新的基于馬氏距離的二元假設(shè)檢測方法,用于表征和檢測航天器異常。Bao等人提出了一種基于貝葉斯理論和多模型的濾波器,使用最小均方誤差方法來適應(yīng)機動目標(biāo)可能進(jìn)行的多種運動,以此進(jìn)行聯(lián)合檢測和弱目標(biāo)跟蹤。Clark等人通過將TLE分析和狀態(tài)傳播這兩種眾所周知的檢測方法融合成混合檢測方法來提高檢測的計算速度和精度提高,提高了空間態(tài)勢感知(space situational awareness, SSA)算法的性能。Justin等人利用光聲特征檢測方案實現(xiàn)了航天器近實時機動檢測和參數(shù)估計,增強了對小機動事件的識別能力。Jiang等人探討了跟蹤執(zhí)行未知連續(xù)機動的衛(wèi)星問題,提出一種估計衛(wèi)星狀態(tài)和機動加速度的魯棒擴展卡爾曼濾波器,與無偏最小方差輸入和狀態(tài)估計方法相比,該方法在具有相同的位置估計精度的同時機動期間速度估計誤差降低了5倍左右。Christopher從目標(biāo)歷史TLE記錄中獲得的數(shù)據(jù)訓(xùn)練神經(jīng)網(wǎng)絡(luò),用于發(fā)現(xiàn)新的TLE記錄中的異常軌道機動、錯誤標(biāo)記或其他原因引起的意外偏差,該方法能夠避免依賴于個別目標(biāo)軌道狀態(tài)或特定的傳播模型,提高了檢測的準(zhǔn)確性。Li等人利用目標(biāo)的歷史TLE數(shù)據(jù),分別用統(tǒng)計的方法和自組織映射聚類的方法對LEO軌道衛(wèi)星的機動檢測問題進(jìn)行了研究。崔紅正等人提出針對不同推力非合作目標(biāo)的機動檢測方法,結(jié)合地基與天基觀測數(shù)據(jù),可以滿足多數(shù)非合作目標(biāo)軌道機動檢測需求。王慶瑞等人將距離變化率作為特征量,提出一種基于概率判決模型的自適應(yīng)軌道機動檢測方法,該方法基于Neyman-Pearson準(zhǔn)則求解序列數(shù)據(jù)的變軌判決門限。

以上方法中基于歷史軌道信息的機動檢測方法只能判斷是否發(fā)生機動,不能對機動參數(shù)進(jìn)行有效的定量識別,而對于機動信息的識別,大多是基于天基或地基平臺對目標(biāo)的連續(xù)測角或測距追蹤,結(jié)合動力學(xué)模型和濾波方法進(jìn)行機動軌道定軌及機動識別。但對于多數(shù)非合作目標(biāo),難以對其進(jìn)行連續(xù)追蹤觀測并且動力學(xué)模型是未知的,只能通過稀疏短弧段觀測數(shù)據(jù)在無機動假設(shè)下定軌得到若干個離散的軌道狀態(tài)。而當(dāng)衛(wèi)星發(fā)生軌道機動時,通過軌道測量得到的衛(wèi)星軌道根數(shù)相比于不發(fā)生軌道機動時的軌道根數(shù),通常含有較大的誤差,且在短時間內(nèi)有限條軌道根數(shù)無法通過濾波等手段處理誤差。

針對以上問題,本文基于稀疏短弧段觀測數(shù)據(jù)定軌得到的有限條軌道測量狀態(tài),對非合作飛行器的機動參數(shù)進(jìn)行識別。首先建立基于稀疏軌道根數(shù)的飛行器軌道反演模型,然后基于軌道反演對飛行器機動參數(shù)進(jìn)行識別,最后通過算例對本文的方法進(jìn)行校驗。

1 非合作飛行器軌道反演模型

1.1 基本假設(shè)

在本文的研究中,做以下假設(shè)。

(1) 當(dāng)變軌推力作用時間遠(yuǎn)遠(yuǎn)小于變軌前后的軌道周期時,可以假設(shè)推力隨時間變化的函數(shù)是脈沖函數(shù),使得脈沖矢量與原推力產(chǎn)生的沖量一樣,因此本文中涉及的變軌沖量均考慮為脈沖沖量。

(2) 軌道測量僅存在隨機誤差,不存在系統(tǒng)誤差。

(3) 以美國戰(zhàn)略司令部發(fā)布的TLE中軌道機動后恢復(fù)時間為例,空間飛行器在發(fā)生機動后,通過測定軌得到穩(wěn)定的TLE軌道通常需要2～7天,本文假設(shè)飛行器在得到穩(wěn)定的軌道之前不會施加下一次軌道機動。

(4) 在非合作飛行器的一個軌道周期內(nèi),至少能得到一段短弧段觀測數(shù)據(jù),結(jié)合假設(shè)(3),即對于軌道周期小于2天的非合作飛行器,相鄰短弧段觀測數(shù)據(jù)之間至多考慮存在一次軌道機動。

1.2 沖量作用下飛行器軌道預(yù)報模型

在二體假設(shè)下,天體為質(zhì)點,空間飛行器只受中心引力作用。而實際飛行任務(wù)中,中心天體一般為非理想球體,因此飛行器還會受中心天體非球形引力、大氣阻力、太陽光壓、第三體引力等各種攝動因素的影響。

在J2000地心慣性坐標(biāo)系中的飛行器動力學(xué)模型如下所示:

(1)

當(dāng)飛行器的初始運動狀態(tài)可知,在脈沖機動作用下,飛行器在任意時刻的軌道狀態(tài)表示如下:

(2)

1.3 飛行器軌道反演

(3)

定義軌道反演的狀態(tài)偏差為在所有測量點時刻(=0,1,…,),反演軌道狀態(tài)與測量狀態(tài)的加權(quán)偏差和,即

(4)

(5)

=DU/VU

(6)

(7)

式中: ,(=,,)為對應(yīng)于沿軌道不同方向的測量誤差的比例系數(shù)。

狀態(tài)偏差反映了通過軌道測量數(shù)據(jù)反演的軌道與測量軌道的接近程度,狀態(tài)偏差越小,說明反演的軌道與測量軌道更接近,同時與測量軌道的離散狀態(tài)相比,通過測量狀態(tài)反演的軌道是一條連續(xù)的軌道,減小了測量數(shù)據(jù)中的隨機誤差的影響,反演得到的軌道更加接近真實軌道。因此,可以通過最優(yōu)化方法使?fàn)顟B(tài)偏差最小,以此反演非合作飛行器的軌道。

2 基于軌道反演的機動參數(shù)識別方法

2.1 機動參數(shù)識別優(yōu)化模型

(1) 優(yōu)化變量

由式(3)和式(4)可知,狀態(tài)偏差與反演軌道的初始軌道狀態(tài)和中間機動參數(shù)有關(guān),即

(8)

因此,將優(yōu)化變量選為反演軌道初始狀態(tài)和中間機動,即

(9)

(2) 優(yōu)化變量邊界約束

(10)

式中:,,,,,分別為飛行器軌道的半長軸、偏心率、軌道傾角、升交點赤經(jīng)、近地點角距和真近點角;和分別為根據(jù)飛行器軌道類型設(shè)定的半長軸邊界。

(3) 目標(biāo)函數(shù)

目標(biāo)函數(shù)取為狀態(tài)偏差最小,即

(11)

值得注意的是,由于在設(shè)計變量中引入的中間機動通常是多于實際機動的,多余的機動在優(yōu)化結(jié)果中為小值,因此在優(yōu)化流程的最后進(jìn)行一次小機動值的修正,即若存在設(shè)計變量|Δ|(=1,2,…,)的數(shù)值小于一定門限,認(rèn)為在該時刻未發(fā)生機動。

(12)

2.2 “全局尋優(yōu)+局部修正”的雙重優(yōu)化求解

上述問題可以歸結(jié)為非線性規(guī)劃問題,關(guān)于這種非線性規(guī)劃問題的求解,有很多算法可以使用,如差分進(jìn)化(differential evolution, DE)算法、序列二次規(guī)劃(sequential quadratic programming, SQP)算法等。此處選取DE算法和SQP算法,其中DE算法全局搜索能力較好,但效率不高,尤其當(dāng)軌道外推使用高精度模型進(jìn)行時,會花費大量時間;序列二次規(guī)劃算法局部尋優(yōu)效果好,但對初值比較敏感。因此,在使用高精度軌道外推模型求解時,為了提高搜索效率以及避免初值陷入局部最優(yōu),提出一種“全局尋優(yōu)+局部修正”的雙重優(yōu)化算法進(jìn)行計算,計算流程如圖1所示。

圖1 雙重優(yōu)化求解流程

在二體條件下,通過DE算法尋找滿足約束的全局最優(yōu)解;

以步驟1中最優(yōu)解為初值,利用SQP算在高精度模型條件下,尋找滿足約束的最優(yōu)解;

剔除小機動值,重新進(jìn)行一遍SQP優(yōu)化流程。

3 仿真校驗

本節(jié)對上述的飛行器軌道反演模型及機動參數(shù)識別方法進(jìn)行仿真校驗。仿真在計算機處理器為Intel i7-7700 CPU 3.60 GHz,運行內(nèi)存為16 GB的C++環(huán)境下進(jìn)行。

3.1 仿真參數(shù)設(shè)定

地球同步靜止軌道由于具有獨特的高軌和地球同步特性,成為了通信、氣象、導(dǎo)彈預(yù)警、廣播電視和數(shù)據(jù)中繼等重要衛(wèi)星的首選軌道,但同樣也將面臨嚴(yán)峻的惡意服務(wù)考驗。因此,對地球同步靜止軌道衛(wèi)星機動信息的感知、檢測和識別顯得尤為重要。本節(jié)以地球同步靜止軌道為例,對本文提出的模型和方法進(jìn)行校驗。

表1 仿真初始軌道根數(shù)

表2 標(biāo)稱軌道中間機動參數(shù)

對上述標(biāo)稱軌道每隔7 h選取該時刻的軌道狀態(tài)作為測量樣本,并施加不同的均值為0的高斯隨機定軌測量誤差,以此模擬實際情況下稀疏短弧段定軌得到的有限條離散軌道測量數(shù)據(jù)。不同工況下的定軌誤差如表3所示。

表3 定軌誤差標(biāo)準(zhǔn)差

3.2 算例分析

使用本文描述的機動參數(shù)識別方法,不同測量點時刻對應(yīng)狀態(tài)偏差的權(quán)重因子均取為1,優(yōu)化目標(biāo)函數(shù)取為式(7),目標(biāo)函數(shù)中沿軌道不同方向的測量誤差的比例系數(shù), (=,,)取值如下所示：

(13)

其中，取為015 m/s,DE算法種群規(guī)模為240,代數(shù)為300,SQP算法最大迭代代數(shù)為50。

每種工況分別進(jìn)行50次蒙特卡羅仿真,機動參數(shù)識別結(jié)果與標(biāo)稱機動參數(shù)的偏差如表4所示。從結(jié)果來看,定軌誤差越大,則機動參數(shù)的識別誤差越大。

表4 機動參數(shù)識別結(jié)果

為進(jìn)一步說明雙重優(yōu)化算法的收斂情況及機動參數(shù)識別效果,選取工況2中的一次仿真,初始軌道測量數(shù)據(jù)如表5所示。雙重優(yōu)化算法的迭代收斂情況如圖2所示,仿真顯示雙重優(yōu)化算法中DE算法在100代左右基本收斂,SQP算法在20代左右基本收斂。

圖2 雙重優(yōu)化算法收斂情況

表5 離散軌道測量數(shù)據(jù)

在表5中，定軌結(jié)果以經(jīng)典軌道根數(shù)的形式給出,括號中數(shù)據(jù)依次代表半長軸、偏心率、軌道傾角、升交點赤經(jīng)、近地點角距、真近點角,其中半長軸的單位為m,角度的單位為(°)。

同時,為比較說明本文優(yōu)化算法的優(yōu)勢,分別用DE算法和SQP算法對算例進(jìn)行求解。表6給出了不同算法求解的最大機動識別誤差和計算時間,比較可知,SQP算法機動識別誤差較大,原因是SQP算法對初值較為敏感,隨機給定初值容易使算法陷入局部收斂;DE算法計算時間較長,原因在于直接使用高精度軌道外推模型進(jìn)行全局搜索時,由于計算量驟增,導(dǎo)致計算效率極低。相對而言,本文提出的“全局尋優(yōu)+局部修正”的雙重優(yōu)化算法對本文的模型適應(yīng)性很好。

表6 不同算法機動識別誤差及計算時間

經(jīng)SQP高精度修正后的飛行器機動參數(shù)反演結(jié)果如表7和表8所示,分別給出了剔除小機動前后的機動參數(shù)識別情況?？梢钥闯?剔除小機動后的識別結(jié)果與不剔除小機動相比結(jié)果更加接近標(biāo)稱軌道施加的機動,機動識別的速度誤差在0.08 m/s以內(nèi),時間誤差在3 min以內(nèi)。

表7 飛行器機動參數(shù)識別結(jié)果

表8 飛行器機動參數(shù)識別結(jié)果(剔除小機動)

為進(jìn)一步說明反演軌道與實際軌道的偏差,圖3給出了在仿真時段反演得到的軌道與標(biāo)稱軌道的位置偏差及軌道跡向、法向和徑向偏差。可以看出,在整個仿真時段,反演軌道與標(biāo)稱軌道的最大位置偏差不超過1 200 m。

圖3 反演軌道與標(biāo)稱軌道的偏差

4 結(jié) 論

本文在稀疏軌道信息條件下,針對非合作飛行器軌道機動參數(shù)的識別問題,探討了一種可行的思路。建立了飛行器軌道反演模型及基于軌道反演的機動參數(shù)識別方法,并提出相應(yīng)的“全局尋優(yōu)+局部修正”的雙重優(yōu)化求解算法。仿真結(jié)果顯示,該方法能夠在稀疏軌道信息的條件下,對非合作飛行器的機動參數(shù)進(jìn)行較為有效的識別;空間飛行器定軌誤差越大,則機動參數(shù)的識別誤差越大;工況2中算例機動參數(shù)識別結(jié)果顯示,剔除小機動后的識別結(jié)果與不剔除小機動相比結(jié)果偏差更小,反演軌道與標(biāo)稱軌道的位置偏差不超過1 200 m,機動識別的速度誤差在0.08 m/s以內(nèi),時間誤差在3 min以內(nèi)。