面向小推力變軌的天文組合自主導航方法

馬廣富，王 偉,，張 偉，黃慶龍，彭玉明，張 曉

(1. 哈爾濱工業大學控制科學與工程系，哈爾濱 150001；2. 上海衛星工程研究所，上海 200240； 3. 上海市深空探測技術重點實驗室，上海 200240)

0 引 言

小行星探測活動是當前國內外航天活動的重點和熱點，包括小行星探測在內的多次深空任務也已納入國家深空重大專項。小行星探測任務由于能量需求大，所以高比沖、長壽命、高可靠的電推進系統是工程任務實施的優選動力形式。采用電推進方式能有效增強探測器變軌能力，提高任務靈活性，降低任務成本等優勢，但電推進方式給探測器帶來一些新問題。電推進方式推力較小，需要長期連續開機工作，期間為確保推力方向正確且帆板對日定向，探測器與地面的通信無法保證連續可用，地面無線電導航方法無法為探測器提供可靠的導航信息。此外，電推進方式對連續小推力建模的不準確性，只依靠探測器自主軌道遞推等方法應用受限，因此，采用小推力變軌的小行星探測器亟需一種能夠提供實時連續、高精度導航信息的自主導航方法。

當前小行星探測工程任務主要依賴于天文光學圖像測角方法進行自主導航，該方法對觀測的目標天體約束限制較多，導致小行星探測器的導航信息不能滿足實時連續高精度的需求[1-3]。本文在近年來深空導航領域新提出的直接獲得探測器速度信息方法的基礎上[4-6]，結合測角導航方法，給出了一種適用于小推力變軌小行星探測器的組合自主導航方法。

基于電推進方式的小行星探測器，由于存在太空環境未知因素多、動力學環境變化大等因素，導航系統的狀態具有不確定性強、狀態模型變化較大等特點，這些因素對導航系統狀態估計性能會有很大影響[7]。目前，國內外學者針對導航系統狀態模型多變的情況，提出了采用多模型的狀態估計方法[8-10]。交互式多模型算法(Interacting multiple model，IMM)是在二十世紀九十年代由Blom和Bar-Shalom提出，同時他們也對IMM算法做了大量改進研究[11-12]。IMM算法計算消耗適中，通過對各子濾波器估計結果融合輸出作為濾波算法的初值，同時利用輸出結果的交互概率實現各子濾波器的初始化，達到基于不同模型組合跟蹤實際導航系統變化的效果。

在上述研究基礎上，學者們對IMM算法提出了很多的改進[13-15]。Johnston和Krishnamurthy通過重新計算權值給出了新的IMM算法[16]。Li等[17]通過對IMM算法結構進行改進，提出了變結構的交互式多模型(VSIMM)算法，但是由于結構復雜，計算量加大，在實際應用中有一定的局限性。近年來，徐田來等融合參數自適應與IMM算法得到了新的IMM算法[18]，但由于加入了參數自適應，算法的計算量較大。有些學者通過引入神經網絡等智能算法計算IMM算法的模型轉移概率[19-20]，或者對各個子模型的參數做修正，進而更準確地對系統狀態進行估計[21-22]，上述這些針對IMM的改進方法計算效率不高，不利于在實際工程中應用。

當前針對基于電推進方式變軌的小行星探測器，在推力矢量不確定性較大、動力學環境未知因素多等情況下的導航狀態估計問題研究尚少。本文針對持續小推力變軌的小行星探測器對導航系統自主性強、實時性高的需求，研究了一種小推力變軌的天文測角測速組合導航方法。首先，根據工程實踐分析并建立了電推進變軌過程中的動力學模型，給出了天文光學圖像測角融合天文光譜測速的組合導航方案；為了克服系統模型集合小、先驗信息的不準確性對導航系統估計精度的影響，提出結合Sage-Husa自適應估計器的AIMM算法，以較少的的模型個數實現對實際狀態的覆蓋；同時結合UKF算法給出AIMM-UKF算法，以提高組合導航系統的魯棒性和抗干擾能力。最后，通過數學仿真對天文測角測速組合導航系統的性能進行對比驗證，分析本文所提出的導航方法的有效性。

1 天文組合導航系統模型

1.1 導航系統狀態方程

依據軌道動力學分析，建立小行星探測器天文測角測速組合自主導航系統的狀態方程。本文以小行星探測任務的巡航段為對象，在動力學建模時主要考慮的因素有探測器在飛行過程中的太陽、地球、小行星等天體引力，太陽光壓攝動，探測器的修正推力等。

以J2000日心黃道慣性坐標系為基準，小行星探測器位置矢量設為r=[x,y,z]T，速度矢量設為v=[vx,vy,vz]T，則組合導航系統的狀態變量可寫為X(t)=[x,y,z,vx,vy,vz]T，此時系統狀態方程可表示如下：

(1)

在第二行的方程中，第一項是太陽引力攝動，μS是太陽引力常數；第二項是主要行星的第三體引力攝動，ri表示第i個考慮的行星位置矢量，ui表示該行星的引力常數，M=1,2,3,…是所考慮的行星數量編號，rSi表示行星相對小行星探測器的位置矢量；第三項是探測器的太陽光壓攝動。前三項是深空探測器一般都需要考慮的影響因素，第四項是探測器推力加速度項，一般地對推力加速度項是按某個小時段脈沖變軌的形式處理。但對于持續小推力變軌的小行星探測器，由于推力矢量長時間存在，需要考慮較多長期變化因素影響，有必要對模型做進一步的細化。

對于采用電推進方式的小行星探測器，根據工程實踐分析經驗，其推力加速度項模型可表示為：

(2)

式中：第一項是探測器的修正加速度，U是發動機的推力矢量，m是航天器的瞬時質量；第二項是推力常值偏差；第三項是在推力常值偏差外，呈有周期性緩變的推力偏差；最后一項是其他未建模誤差，可按照高斯白噪聲處理。Isp為發動機比沖；g0為海平面重力加速度。

基于上述分析，小行星探測器巡航段的軌道動力學方程如下：

(3)

式中:W是狀態模型噪聲。

1.2 導航系統量測方程

圖1 測角導航原理示意圖Fig.1 Principle of angle measurement navigation

在小行星探測器飛行過程中，導航目標天體的視星等將隨著探測器逐步靠近小行星而慢慢變小，探測器上可通過高精度光學圖像導航敏感器得到目標天體的圖像，通過提取目標天體相對小行星探測器的視線方向矢量，即可獲得探測器相對目標天體的角度信息。所得到的量測信息作為濾波估計器的輸入即為天文測角自主導航[17]。以小行星和太陽為目標天體的測角導航方法其原理如圖1所示。

基于上述分析，測角導航量測方程如下：

(4)

式中：r表示探測器在J2000日心黃道慣性系的位置矢量；lps表示探測器敏感期測量得到的太陽視線矢量、lpa表示探測器敏感期測量得到的小行星視線方向矢量；Vps、Vpa表示對應量測系統噪聲。

圖2 測速導航原理示意圖Fig.2 Principle of velocity measurement navigation

小行星探測器在飛行過程中，相對空間中的太陽或者恒星會因為相對速度的多普勒效應，探測器上的光譜敏感器接收的光譜波長會產生漂移。依據多普勒原理，探測器敏感器在飛行過程中接收波長的漂移量與其靜止時接收到的波長之比等于探測器相對太陽或者恒星的視向速度與光速之比?；谏鲜龇治?，可得到以小行星探測器相對太陽、2顆恒星的視向速度為量測的導航方案，其原理如圖2所示。

測速導航量測方程如下：

(5)

式中：vSpe1，vSpe2，vSpe3分別表示探測器上光譜敏感器測得的相對太陽及2顆恒星的視向速度數值；r，v表示探測器相對太陽的位置、速度矢量；vStar1，vStar2和lStar1，lStar2分別表示2顆恒星在J2000日心黃道慣性系下的速度和視線方向矢量；V表示量測系統噪聲。

綜上，小行星探測器的導航量測模型可用如下一般形式描述：

Z=h(X,t)+V

(6)

2 組合導航算法設計

基于電推進方式變軌的小行星探測器，由于推力持續時間長、推力影響因素多、探測器的動力學環境變化大等因素，前述的組合導航系統不確定性強、狀態模型變化較大，這些因素對導航系統狀態估計精度會有很大影響。IMM算法利用多個子模型的模型集來描述系統可能的運動狀態，對小行星探測器導航系統狀態具有很好的跟蹤效果。下文基于IMM算法，對其開展研究分析并給出改進的導航濾波算法。

2.1 IMM算法及其改進

交互式多模型算法包括四個部分：輸入交互、濾波估計、模型更新及融合輸出。

IMM算法利用多個子模型描述導航系統的可能運動狀態，通過計算不同模型概率及其轉移概率，完成導航濾波估計過程中的輸入交互。輸入交互后的數據經過導航系統濾波后，更新不同導航模型概率權值，并進而得到多個模型狀態和對應協方差矩陣的融合值，最后通過濾波估計器融合輸出導航系統狀態的估計結果。

與單模型及其他多模型處理方法相比，IMM-UKF算法具有很強的交互性和自適應性，具有很高的效費比，有利于處理小行星探測任務中具有持續推力變軌、動力學環境未知因素多等情況復雜的導航狀態估計問題。但IMM-UKF算法本身也存在一些局限性，由于導航系統模型概率是基于先驗信息給出[23-24]，而先驗信息往往由于環境不確定性、存在量測噪聲等因素，具有一定的不準確性，對導航系統的濾波估計精度會有較大的影響。

對于測角測速組合導航系統，為解決上述不準確性影響，保證交互式多模型中的模型集合覆蓋系統的所有可能狀態，IMM算法通常需要建立較大的模型集合，但是過多的模型會消耗探測器上很多的計算資源。為了減輕導航算法的計算量，匹配小行星探測器的計算能力，通過引入自適應參數估計器，在實際模型參數附近逼近真實的模型參數，實現對IMM算法進行改進并得到AIMM算法。AIMM算法可以利用較小的模型集合實現對實際組合導航狀態模型的覆蓋。

針對天文測角測速組合自主導航系統，本文將針對線性系統的Sage-Husa次優極大后驗噪聲估計器推廣到非線性系統[25-26]，得到非線性時變噪聲無偏遞推估計器，具體計算如下：

(7)

(8)

式中：

(9)

(10)

(11)

(12)

式中：遺忘因子b∈[0.95, 0.99]。

AIMM算法將利用Sage-Husa自適應估計器，對導航系統量測噪聲誤差進行自適應估計，從而增強IMM算法輸出對導航系統狀態模型的逼近，減小由于模型較少導致導航估計精度下降的影響。

2.2 濾波器模型集設計

在小行星探測器電推進持續變軌過程中，由于導航系統動力學環境變化大、推力矢量模型不精準等因素，導航系統的狀態模型變化較大。結合前述AIMM算法，考慮在導航系統狀態模型上引入一組系統噪聲和觀測噪聲，設計濾波器的模型集。將式(3)和式(6)離散化處理，并引入子模型記號mi，則有如下濾波器模型集：

(13)

式中：W(k,mj)和V(k+1,mj)分別是系統狀態模型和量測模型的噪聲；j=1,2,3，…，N表示不同子模型序號，通過選取不同大小的量測噪聲得到不同的子模型。

2.3 AIMM-UKF算法設計

圖3 AIMM-UKF算法原理框圖Fig.3 Principle block diagram of AIMM-UKF

結合AIMM算法及導航系統狀態模型集合的設計，下面給出AIMM-UKF算法的具體過程。

1)輸入交互

首先對各個子模型進行初始化并做交互運算，將前一時刻探測器系統狀態的交互運算結果作為各個子濾波器的輸入，此時導航系統狀態的估計值和協方差陣為：

(14)

(15)

式中:μji(k+1|k)是預測模型轉移概率：

(16)

2)濾波估計

(17)

(18)

同時，結合自適應參數估計，可得：

(19)

濾波增益為：

(20)

更新各個子濾波器的系統狀態及其協方差陣：

(21)

Rj(k)](Kj(k))T

(22)

3)模型更新

結合上述步驟，計算不同子濾波器的新息及其對應協方差陣，并同時更新μji(k+1|k)。

各個子濾波器的新息及其協方差陣如下：

(23)

(24)

更新模型概率時，在k時刻與第j模型匹配的似然函數為：

(25)

(26)

上式中的C為歸一化常數，表達式為：

(27)

4)融合輸出

以更新后的模型概率μji(k+1)作為各個子模型與探測器真實運動狀態逼近程度的指標，對當前各個子模型的輸出進行加權融合，得到的結果作為導航系統最終的輸出：

(28)

(29)

3 仿真分析

3.1 仿真條件及參數選取

本文仿真所需要的基本參數歸結如下。

1)以J2000日心黃道慣性坐標系為參考基準，選取熱點探測目標谷神星作為小行星探測工程任務背景開展數學仿真分析，檢驗本文所提出的一種小推力變軌的天文測角測速組合自主導航方法的有效性。仿真中的真實軌道數據由STK產生，組合導航系統軌道動力學模型包括太陽引力攝動、光壓攝動及探測器長期持續的電推力等影響。谷神星等天體的參數通過DE406星歷庫獲得。小行星探測器的相關參數如表1所示。

表1 小行星探測器巡航段參數Table 1 Parameters of asteroid probe cruising period

2)測速導航方案中的2顆恒星數據通過依巴谷星表獲得，對應的編號分別為35468和60178，相應的星歷數據及誤差如表2所示。

表2 測速導航恒星參數Table 2 Parameters of stars for velocity navigation

3)視線方向矢量量測中測角精度為1″(3σ)，太陽光譜測速儀測速精度為1 m/s (3σ)，恒星光譜測速儀測速精度為2 m/s (3σ)。

4)相關濾波算法的基本參數設置如表3所示。

5)組合導航系統在交互式多模型中選取4個子模型，模型參數設置如下：

表3 濾波算法參數設置Table 3 Parameters of filtering algorithms

(30)

(31)

多模型集合設置如下：

模型1：Q1=0.01Q；R1=100R。

模型2：Q1=100Q；R1=100R。

模型3：Q1=0.01Q；R1=0.01R。

模型4：Q1=100Q；R1=0.01R。

本文對不同濾波算法采用Monte Carlo方法進行仿真，那么每次濾波估計的均方根誤差可以記為：

(32)

在完整的時間序列下，均方根誤差可以如下描述：

ψ=[ψ1,ψ2,…,ψn]T

(33)

為了更準確衡量導航濾波算法的精度，在濾波估計收斂后，選取濾波過程末期的L個點的ψ的均值作為整個導航濾波算法的精度指標(3σ)，即：

(34)

通常算法的復雜度是判斷該算法是否可用于實時處理的依據，在實時性要求高的自主導航系統中作為一項參考指標。算法的復雜度是指其自身在計算中所需的乘法次數，其中不含動力學、量測等產生的時間。在硬件相同時，可采用每次運行的總時間作為濾波算法復雜度的直觀度量。

3.2 仿真結果及分析

基于AIMM-UKF算法的測角測速組合自主導航位置與速度估計結果如圖4和圖5所示，從圖中可以看出，組合導航估計很快收斂。

圖4 AIMM-UKF算法三軸位置估計結果Fig.4 Three axis position estimation results of AIMM-UKF

圖5 AIMM-UKF算法三軸速度估計結果Fig.5 Three axis velocity estimation results of AIMM-UKF

在仿真算例中，AIMM-UKF采用了4個不同的模型，隨著導航濾波逐步遞推，AIMM-UKF算法對模型權值不斷更新調整。在外部動力學環境不確定情況下，本文所給出的AIMM-UKF算法模型轉移概率很快收斂穩定，算法可以很快地匹配到近真實的系統模型，即算法具有很高的自適應性。

為了對比所給出的AIMM-UKF算法性能，下面分別從導航估計結果的精度、計算效率指標上，對比分析AIMM-UKF算法和普通UKF算法、IMM-UKF算法的性能，如圖6和圖7所示。

圖6 不同導航濾波算法的位置估計結果Fig.6 Position estimation results of different filtering algorithms

圖7 不同導航濾波算法的速度估計結果Fig.7 Velocity estimation results of different filtering algorithms

從上述不同導航濾波算法的估計結果可知，在導航濾波估計的初期，AIMM-UKF算法由于受到初始協方差矩陣權重分配、初始概率分布等因素與實際狀態不匹配的影響，濾波估計會有較大超調情況，但隨著濾波過程中多模型概率分布的更新，AIMM-UKF算法能夠更好地逼近真實的動力學模型和量測模型，它會比普通的UKF和IMM-UKF更快地收斂到穩態，最終AIMM-UKF算法的位置和速度估計精度也會更高。

不同導航濾波算法的估計結果如表4所示。分析可知，普通UKF算法由于模型偏差較大，在導航估計后期逐漸震蕩發散；IMM-UKF算法導航估計結果較好；本文所給出的AIMM-UKF算法的位置和速度估計精度比IMM-UKF算法的位置估計精度和速度估計精度分別提升了30.76%和40.55%。

表4 不同導航濾波算法的性能對比Table 4 Performance comparison of different navigation filtering algorithms

4 結 論

針對電推進形式持續小推力變軌的小行星探測任務對導航系統自主性強、實時性高的需求，根據工程實踐分析并建立了電推進變軌過程中的動力學模型，給出了天文光學圖像測角融合天文光譜測速的小行星探測組合導航方法。為克服小行星探測器推力的不確定性，提出了采用自適應交互式多模型無跡卡爾曼濾波(AIMM-UKF)算法，以較少的模型個數實現對導航系統狀態的覆蓋，克服了模型集合先驗信息不準確對導航精度的影響，提高了組合導航系統的魯棒性和抗干擾能力。最后，通過基于小行星探測工程任務巡航段的仿真分析，對測角測速組合導航系統的性能進行了對比驗證，結果表明本文所提出的組合導航方法估計精度更高、計算消耗更小，可滿足小行星探測工程任務對自主導航的高性能需求。