普渝皓，程實，陳楊*

(1.航天科工空間工程發(fā)展有限公司，北京100854；2.中國航天科工集團二院研究生院，北京100854)

1 引言

近地小行星對地球的威脅持續(xù)存在，特別是大量未發(fā)現(xiàn)編目的近地小行星對地球持續(xù)造成嚴重突發(fā)撞擊威脅。從歷史情況看，直徑1km以上小行星撞擊地球概率較低且短期內(nèi)難以有效防御，直徑10m以下小行星危害較小無需針對性進行防御，因此監(jiān)測預(yù)警的重點目標為直徑10m至1km的近地小行星。進一步考慮工程可行性、應(yīng)對必要性等因素，防御任務(wù)重點為直徑十米級至百米級的小行星。近年來，小行星防御技術(shù)成為研究熱點，其中任務(wù)末段的相對導(dǎo)航技術(shù)與制導(dǎo)技術(shù)是關(guān)鍵技術(shù)[1]。

“雙小行星重定向測試”(Double Asteroid Redirection Test, DART)任務(wù)于北京時間2022年9月27日，對Didymos雙小行星系統(tǒng)中直徑160m的次星Didymos B進行動能防御，仿真命中點距瞄準點約15m[2]，實際命中點距離圖像中心17m。根據(jù)公開報道，DART任務(wù)主飛行器經(jīng)過深空機動與多次中制導(dǎo)修正，滿足中末交接班條件后，進入末制導(dǎo)階段，在此階段使用源于比例引導(dǎo)的制導(dǎo)律，制導(dǎo)算法使用擴展卡爾曼濾波(Extended Kalman Filter, EKF)對交會時刻的投影脫靶量進行估計[3]。濾波器的觀測量是由圖像解算出的慣性坐標下的視線角，并且假設(shè)任務(wù)飛行器與目標小行星之間無相對加速度。根據(jù)投影脫靶量大小進行修正，在交會前5min進行連續(xù)修正，并在交會前2min?？兀瑴蚀_命中目標小行星。面對目標尺寸更小、目標軌道預(yù)報精度更低、交會精度要求更高的任務(wù)場景，上述方法還有改進的空間。

首先，考慮提升濾波算法的精度。以視線角為觀測量的相對導(dǎo)航系統(tǒng)是典型的非線性系統(tǒng)。EKF算法在線性化過程中忽略了非線性函數(shù)泰勒級數(shù)展開的二階及以上的高階項。在初始誤差較大的情況下，特別是在目標軌道預(yù)報精度較低的情況下，有可能引入較大的線性化誤差。無跡卡爾曼(Unscented Kalman Filter, UKF)算法不需要對狀態(tài)方程和觀測方程線性化，沒有對模型高階項的截斷誤差，可以獲得比EKF更好的估計精度和濾波穩(wěn)定性[4]。EKF和UKF算法均需要已知系統(tǒng)噪聲和觀測噪聲的統(tǒng)計特性。噪聲統(tǒng)計特性的偏差將導(dǎo)致濾波精度下降甚至導(dǎo)致濾波發(fā)散。為了克服這一問題，需要引入在線估計噪聲統(tǒng)計特性和修正濾波增益的自適應(yīng)濾波技術(shù)，如基于協(xié)方差匹配的自適應(yīng)濾波[5]、基于極大似然方法的自適應(yīng)濾波[6,7]和基于模糊邏輯的自適應(yīng)濾波[8,9]等。

其次，考慮提升制導(dǎo)指令的計算精度。在公開文獻中，DART任務(wù)在進行投影脫靶量估計時，忽略了相對加速度對制導(dǎo)精度的影響[2]。根據(jù)Lorenzo F的研究表明，DART任務(wù)在四體動力學(xué)環(huán)境下引起的軌道誤差均值為238m，目標軌道不確定性引起的平均軌道誤差為10m[10]。因此，為應(yīng)對更小尺寸小行星的更高精度的防御場景，在投影脫靶量的估計時，需要考慮復(fù)雜引力環(huán)境產(chǎn)生的相對加速度。

本文首先以B平面坐標系作為相對導(dǎo)航坐標系，并給出該坐標系下的自主相對導(dǎo)航系統(tǒng)模型。其次將適用于線性系統(tǒng)的Sage-Husa噪聲估計器推廣到非線性系統(tǒng)，結(jié)合UKF算法得到一種自適應(yīng)無跡卡爾曼(Adaptive Unscented Kalman Filter, AUKF)濾波器，并引入?yún)f(xié)方差匹配的思想，進一步提高濾波的快速性和穩(wěn)定性。在投影脫靶量的估計計算中，考慮復(fù)雜引力環(huán)境下的影響，提出一種利用投影脫靶量直接計算制導(dǎo)指令的制導(dǎo)方法。最后以具有潛在威脅的小行星Bennu為目標，設(shè)計典型任務(wù)場景并生成飛行數(shù)據(jù)，對本文提出的模型及算法進行仿真驗證。

2 自主相對導(dǎo)航方法

2.1 系統(tǒng)模型

2.1.1相對導(dǎo)航坐標系

本文考慮單目僅測角相對導(dǎo)航系統(tǒng)，以目標天體的俯仰角和方位角為觀測量。任務(wù)末段飛行器相對目標天體的位置矢量和速度矢量夾角很小，甚至重合，沿視線方向的信息無法通過連續(xù)觀測獲得。若在J2000日心黃道坐標系中描述狀態(tài)方程，視線方向距離無法估計的特性將反映到三個慣性坐標軸上，導(dǎo)致三軸的位置估計均不可信，無法保證三軸的位置信息的估計精度。因此，需要將沿視線方向的距離獨立出來，提高垂直于視線方向的估計置信度。

B平面坐標系是根據(jù)任務(wù)飛行器相對目標天體進入漸近線方向定義的參考坐標系，如圖1所示，其數(shù)學(xué)定義為：以目標天體為中心原點，S軸為進入軌道雙曲線漸近線的方向矢量，某參考方向矢量記為N，S與N叉乘作為T軸，R軸由S軸和T軸按右手螺旋法則確定。在本文中，S矢量方向可選為初始時刻任務(wù)飛行器相對目標天體的速度V∞的方向。取N為任務(wù)飛行器軌道平面的法線。即：

圖1 B平面坐標系示意圖

(1)

式中：r和v分別表示任務(wù)飛行器在慣性系下的位置和速度矢量。則B平面坐標系慣性坐標系的轉(zhuǎn)換矩陣為：

(2)

視線方向與相對位置方向重合，相對位置方向和相對速度方向幾乎重合。因此B平面坐標系可以將視線方向的距離獨立出來。由于S軸近似表示飛行器進入B平面的方向，T軸和R軸方向的位置信息可以描述任務(wù)飛行器與小行星在B平面內(nèi)的相對位置。故選取B平面坐標系作為相對導(dǎo)航坐標系。

2.1.2 系統(tǒng)狀態(tài)和觀測方程

該任務(wù)的中心天體為太陽，同時飛行器與小行星還要受其他大行星引力及太陽光壓力等攝動力的影響。飛行器在接近目標天體的末段，其與目標天體所處的位置相差不大，二者所受的太陽引力、攝動力相差很小。在不進行機動時，任務(wù)飛行器相對目標小行星近似做勻速直線運動。由于近似所忽略的加速度可視為系統(tǒng)狀態(tài)噪聲。任務(wù)飛行器在B平面坐標系下的離散系統(tǒng)狀態(tài)方程為：

(3)

式中：X=[xyzvxvyvz]T表示任務(wù)飛行器在B平面坐標系下的相對位置和速度；t表示濾波步長，Wk-1表示系統(tǒng)狀態(tài)噪聲。

假設(shè)光學(xué)相對導(dǎo)航相機捷聯(lián)安裝在任務(wù)飛行器本體，輸出在飛行器本體坐標系下描述的目標視線角，包括俯仰角和方位角。因此觀測模型如下：

(4)

式中：qγ和qλ分別表示目標的俯仰角和方位角；Vk表示觀測噪聲；xb、yb、zb表示r在本體坐標系b系下的投影。xb、yb、zb與系統(tǒng)狀態(tài)x、y、z的轉(zhuǎn)換關(guān)系為：

(5)

2.2 AUKF濾波算法

假定非線性系統(tǒng)的狀態(tài)方程和觀測方程具有離散形式，可表示為：

(6)

式中：X表示被估計的隨機向量，Z表示對X的觀測為隨機向量，f和h表示非線性向量函數(shù)，Wk和Vk表示不相關(guān)的噪聲序列，且Wk和Vk滿足：

(7)

式中，δkj表示Kronecker算子。

針對線性系統(tǒng)，Sage和Husa提出了基于極大后驗估計的噪聲統(tǒng)計估值器[11]。將此噪聲統(tǒng)計估值器推廣到非線性系統(tǒng)，可以得到式(6)對應(yīng)的非線性時變噪聲統(tǒng)計無偏遞推估值器為：

(8)

(9)

(10)

(11)

(12)

(13)

基于無跡卡爾曼濾波算法，結(jié)合式(8)～(11)所示的Sage-Husa噪聲統(tǒng)計估值器以及協(xié)方差匹配思想，可以得到自適應(yīng)無跡卡爾曼濾波算法。針對式(6)所示的非線性系統(tǒng)，本文提出的自適應(yīng)無跡卡爾曼濾波算法實現(xiàn)步驟如下：

步驟1，選定濾波初值

(14)

(15)

(16)

(17)

對k=1,2,3…，執(zhí)行：

(18)

i=1,2,…,n

(19)

i=n+1,n+2,…,2n

(20)

步驟3，時間更新

(21)

(22)

(23)

步驟4，觀測更新

(24)

(25)

(26)

(27)

步驟5，狀態(tài)更新

(28)

(29)

(30)

式(28)表示狀態(tài)增益矩陣更新；式(29)表示狀態(tài)變量值更新；式(30)表示狀態(tài)誤差方差陣更新。

步驟6，噪聲更新

需要特別指出的是，當(dāng)系統(tǒng)噪聲方差陣Q和觀測噪聲方差陣R均未知時，由于無法從輸出中區(qū)分系統(tǒng)噪聲和觀測噪聲造成的影響，Sage-Husa噪聲估計器無法同時對二者進行估計[12]。因此，需要在一種噪聲已知的條件下估計另一種噪聲。本文對飛行器相對目標小行星的運動采取了近似處理，忽略的相對加速度為系統(tǒng)噪聲，且系統(tǒng)噪聲為時變噪聲，而觀測噪聲與敏感器的精度和安裝有關(guān)，可以預(yù)先確定。故本文對系統(tǒng)噪聲進行估計：式(8)表示對系統(tǒng)噪聲均值的估計；式(10)表示對系統(tǒng)噪聲方差的估計。

步驟7，濾波發(fā)散抑制

協(xié)方差匹配的基本原理是：檢驗殘差與其理論統(tǒng)計特征的相容性，即在濾波過程中，判斷殘差的統(tǒng)計分布在原假設(shè)條件的相容性，當(dāng)殘差在原假設(shè)下不相容，則需要根據(jù)殘差與理論統(tǒng)計特征一致的原則，對協(xié)方差陣進行修正。

根據(jù)協(xié)方差匹配判據(jù)，濾波器收斂時應(yīng)滿足：

(31)

式(31)中，εk表示殘差序列；γ為事先設(shè)定的可調(diào)系數(shù)，取值范圍為γ≥1;tr[·]表示求跡運算。若式(31)不成立，則對Pk∣k-1做如下修正：

(32)

式(32)中，λk的確定規(guī)則如下：

(33)

(34)

(35)

式(35)中，μ表示衰減系數(shù)，取值范圍為0<μ≤1，一般取值為0.95。該系數(shù)能進一步提高濾波器的快速跟蹤能力。當(dāng)增大μ時，k時刻之前信息所占的比重減小，當(dāng)前殘差的影響增大。由于沿視線方向的分量是不可觀測的，因此在進行修正時，沿視線方向的分量不進行修正。該方法有很強的關(guān)于突變狀態(tài)的跟蹤能力，并且在濾波達到穩(wěn)態(tài)時，仍保持對于緩變以及突變狀態(tài)的跟蹤能力。

綜上所述，AUKF濾波器計算流程圖如圖2。

圖2 濾波器計算流程圖

3 投影脫靶量制導(dǎo)

3.1 投影脫靶量常值模型估計與計算

根據(jù)任務(wù)飛行器深空自主導(dǎo)航和相對導(dǎo)航結(jié)果，能夠確定當(dāng)前時刻任務(wù)飛行器與目標小行星在慣性空間的位置與速度，進而可以確定兩者的引力加速度GT、GI和飛行器與目標小行星之間的引力差GT/I。針對交會過程的末制導(dǎo)段，交會時刻的引力差很小，即GT/I=0。因此，當(dāng)前時刻引力差常值模型用式(36)表示：

(36)

式中：ti表示當(dāng)前時刻，tf表示交會時刻。由式(36)可知，雖然投影脫靶量的計算是一個常數(shù)，但是在每次制導(dǎo)計算過程中，通過C0更新來反映當(dāng)前的引力狀態(tài)。將式(36)帶入的式(37)中且令控制力F=0，并對其進行積分：

(37)

VT/I(t)=C0(ti)(t-ti)+VT/I(ti)

(38)

RT/I(ti)

(39)

式中，aT/I、VT/I(t)和RT/I(t)分別表示飛行器與小行星的相對加速度、相對速度和相對位置。

飛行器與小行星的相對距離到達最小值時，應(yīng)滿足的條件是:

VT/I(t)·RT/I(t)=0

(40)

將式(38)和式(39)帶入式(40)后得到以交會時間為未知量的方程如下所示：

(tf-ti)+RT/I(ti)

(41)

式(41)中，Rmiss(ti)表示ti時刻的投影脫靶量(Projected Miss Distance)。根據(jù)式(36)-式(41)，當(dāng)C0為0的時候常值模型下的投影脫靶量表達式與不考慮相對引力影響的表達式完全一致。另外，式(41)中通過C0引入了用于表達不同引力的高階項，使得表達式能夠更加準確的近似表達引力差的二次項。

3.2 制導(dǎo)律

圖3 比例導(dǎo)引與投影脫靶量的關(guān)系

從圖3中可以看出，投影脫靶量還可以表示為：

(42)

(43)

比例導(dǎo)引法是根據(jù)飛行器的旋轉(zhuǎn)角速度與目標視線角速度成正比得出的一種導(dǎo)引方法。廣義比例導(dǎo)引是針對比例導(dǎo)引的一種改善形式，其一般形式為

(44)

將式(43)帶入式(44)，可以得到基于投影脫靶量的制導(dǎo)指令計算形式：

(45)

令tgo=tf-ti，式(45)可以簡化表示為：

(46)

上述給出的過載指令均為連續(xù)變大小控制量，考慮飛行器無法產(chǎn)生連續(xù)可變的制導(dǎo)過載，需要設(shè)計高精度的控制輸出變換方法，通過脈寬調(diào)制等手段生成離散開關(guān)控制量，驅(qū)動推力器完成實際的推力輸出。采用脈沖調(diào)寬調(diào)頻方法(PWPF)，對推力指令進行分配。對于PWPF中的4個參數(shù)：放大系數(shù)、時間常數(shù)、推力器的開關(guān)閾值，在實際應(yīng)用中，通過試湊的方法經(jīng)仿真試驗來確定。

4 仿真實驗及分析

4.1 仿真場景

基于潛在威脅小行星Bennu構(gòu)建小行星防御典型場景，進行仿真驗證。Bennu開普勒軌道信息如表1所示[13]。以軌道信息標稱值設(shè)計飛行任務(wù)，任務(wù)于2036年9月18日0時開始，交會時間為2037年1月16日0時，交會前4h進入末段導(dǎo)航與制導(dǎo)階段。小行星實際軌道根據(jù)星歷誤差范圍隨機選擇，小行星直徑選擇為30m。為驗證算法在實際應(yīng)用中的有效性，除星歷誤差外引入多種誤差源，進行蒙特卡洛仿真分析，誤差項及數(shù)值如表2所示。

表1 Bennu開普勒軌道信息(時間：59000.0 MJD)

表2 誤差源及誤差值

4.2 仿真結(jié)果與分析

為對比本文提出的算法與其他算法的相對導(dǎo)航效果，在不施加機動的情況下，分別采用EKF算法、UKF算法、不含協(xié)方差匹配的AUKF算法以及含協(xié)方差匹配的AUKF算法(本文提出的完整算法)進行仿真。目標小行星在J2000日心黃道坐標系的初始位置為[-50119639522.77，130393530773.86，13960229459.83]m，初始速度為[-32208.87，-8047.33，-730.88]m/s，任務(wù)飛行器的初始位置為[-50183158838.61，130466720546.96，13974363604.73]m，初始速度為[-27784.43，-13118.19，-1716.38]m/s。在計算濾波器的初值時，目標小行星的初始位置和速度由小行星星歷主值計算得出，任務(wù)飛行器的初始位置和速度根據(jù)真實位置設(shè)置測控誤差，誤差大小如表2所示。本文中采用的濾波初值為[-97946239.90，-241607.05，88444.30，6801.47，0.00，0.00]。濾波器的參數(shù)設(shè)置一致。系統(tǒng)狀態(tài)噪聲方差陣Q0=diag(1×105,1×105,1×105,1×102,1×102,1×102)，觀測噪聲方差陣R0=diag(6×10-8,7×10-8)，初始狀態(tài)誤差方差陣P0=diag(1×108,1×106,1×106,1×102,1×102,1×102)，濾波步長為5s，其他參數(shù)按上文推薦取值給定。由于飛行器沿視線方向不可觀測，在濾波不發(fā)散的情況下，該方向的位置和速度的變化趨勢與按照初始值遞推一致。故僅列出Y軸和Z軸方向的濾波結(jié)果，如圖4、5所示。

圖4 B平面坐標系下位置估計誤差

圖5 B平面坐標系下速度估計誤差

由圖4、5可以看出，EKF算法隨著濾波時間的增加，累計的誤差變大，出現(xiàn)了發(fā)散的情況；UKF算法由于預(yù)先給出的噪聲統(tǒng)計特性存在偏差，估計精度較差；引入噪聲估計器后，濾波精度提高，但是收斂速度變慢；利用協(xié)方差匹配的思想進行修正，濾波收斂時間從約2000s減小至約1000s。交會前300s，y、z方向的位置濾波誤差為6.89m和0.19m，速度濾波誤差為0.0016m/s和0.00028m/s。本文提出的相對導(dǎo)航算法具有快速、高精度的特點。

在上述條件下對本文給出的僅測角導(dǎo)航模型相對導(dǎo)航和投影脫靶量制導(dǎo)算法進行仿真驗證，飛行器質(zhì)量800kg，推力大小為5N，交會前兩分鐘停止控制。制導(dǎo)過程中，B平面坐標系下垂直于視線方向的投影脫靶量的估計值和真實值變化如圖6所示。進行200次蒙特卡羅仿真驗證，最終交會誤差散布圖如圖7所示。最終交會誤差的均方差3σ=5.22m，最大交會誤差為6.3m。

圖6 末制導(dǎo)階段飛行器與導(dǎo)航濾波器表現(xiàn)情況

圖7 B平面內(nèi)命中點散布圖

5 結(jié)論

小行星防御中的目標尺寸更小、目標軌道預(yù)報精度更低、交會精度要求更高等問題對相對導(dǎo)航的估計精度和制導(dǎo)精度提出了更高要求。本文首先采用B平面坐標系作為相對導(dǎo)航坐標系，分離出不可觀的信息，建立僅測角導(dǎo)航的系統(tǒng)狀態(tài)方程及觀測方程。由于觀測方程為非線性方程，本文以UKF算法作為基礎(chǔ)濾波算法，避免了線性化后帶來的截斷誤差?？紤]系統(tǒng)方程含有未知統(tǒng)計特性的噪聲，將Sage-Husa噪聲估計器推廣到非線性形式，對未知噪聲進行在線估計。針對Sage-Husa算法對濾波穩(wěn)定性及快速性帶來的影響，采用協(xié)方差匹配思想和奇異值分解方法，加快濾波收斂，提高濾波的穩(wěn)定性并且保證算法不失效。其次利用相對導(dǎo)航濾波結(jié)果，建立了引力差產(chǎn)生的相對加速度的常值模型，計算考慮相對加速度的投影脫靶量，并設(shè)計了基于投影脫靶量的制導(dǎo)指令計算方法。最后利用典型小行星防御任務(wù)場景進行算法精度驗證。由計算機仿真結(jié)果表明，本文建立的模型及濾波算法可以對可觀測坐標軸方向的狀態(tài)進行快速準確估計?；谕队懊摪辛康闹茖?dǎo)算法，在復(fù)雜引力環(huán)境下，能夠?qū)崿F(xiàn)5.22m(3σ)的交會精度。本文的相關(guān)研究，能夠為我國未來小行星防御任務(wù)相對導(dǎo)航與制導(dǎo)系統(tǒng)的設(shè)計提供一定的參考和理論依據(jù)。