唐青原，王曉磊

北京控制工程研究所，北京 100190

目前以火星為主的一系列深空大天體探測中的導(dǎo)航任務(wù)大多依然依賴地面測控技術(shù)。隨著深空探測中任務(wù)復(fù)雜度的提高，地面測控由于存在大時延、難以全天候觀測、測量精度有限等問題，越來越難以滿足深空探測的導(dǎo)航要求[1]。未來的深空探測導(dǎo)航關(guān)鍵技術(shù)之一，就是在導(dǎo)航需求超過地面測控能力的情況下能夠?qū)崿F(xiàn)星上自主完成較高精度的軌道確定。

深空自主導(dǎo)航常用的測量信息包括測距、測速及測角，其中深空測距的觀測信息是脈沖星脈沖間隔時間，深空測速的觀測信息是天體光譜。這兩種觀測方式精度較高[2]，但測量技術(shù)尚不成熟，還處于理論研究階段。測角導(dǎo)航，即深空光學(xué)自主導(dǎo)航技術(shù)，是一項(xiàng)已在數(shù)個深空探測項(xiàng)目中得到成功應(yīng)用的成熟技術(shù)。其中最早的應(yīng)用出現(xiàn)在阿波羅系列任務(wù)中[3]，利用觀測月面地標(biāo)獲得位置估計(jì)，這一自主導(dǎo)航技術(shù)適用于接近目標(biāo)天體表面的飛行階段。利用導(dǎo)航相機(jī)觀測天體視線信息的自主導(dǎo)航方法，其首次利用見于1986年的旅行者2號天王星飛掠任務(wù)[4]，此后的伽利略探測器在飛掠一些木星衛(wèi)星時也利用了該技術(shù)[5]，但依然是將自主導(dǎo)航作為地面觀測的補(bǔ)充。近些年，小天體探測中光學(xué)自主導(dǎo)航已成為熱門技術(shù)，2004年歐空局發(fā)射的“羅塞塔”探測器[6]，2005年美國發(fā)射的“深度撞擊”探測器[7]，先后在對彗星小行星等小天體的探測活動中，驗(yàn)證或?qū)崿F(xiàn)了自主導(dǎo)航與控制。2018年抵達(dá)小行星101955Benu的OSIRIS-Rex采樣返回式探測器也采用了小天體接近段自主導(dǎo)航技術(shù)[8]。大行星接近段雙曲線軌道由于動力學(xué)特性不同，自主導(dǎo)航技術(shù)目前還未在軌實(shí)現(xiàn)。文獻(xiàn)[9]中提出了一種利用大視角導(dǎo)航相機(jī)和擴(kuò)展卡爾曼濾波的雙曲線軌道光學(xué)自主導(dǎo)航方法，對金星飛越任務(wù)的軌道進(jìn)行了軌道預(yù)報研究。文獻(xiàn)[10]中研究了利用小型測量儀器進(jìn)行火星接近段自主導(dǎo)航的情形。然而目前大多數(shù)大天體雙曲線軌道預(yù)報都需要利用中心天體視半徑信息，但火星大氣稀薄并可能爆發(fā)全球性風(fēng)暴[11]，此時會出現(xiàn)不可忽視的視半徑錯誤測量，嚴(yán)重影響光學(xué)導(dǎo)航精度。去除視半徑信息的天文自主導(dǎo)航算法此前也多有提出，具有代表性的是文獻(xiàn)[12]中提出的一種觀測火星本身或其衛(wèi)星中心視線矢量相對于參考恒星的視線矢量夾角的方法，以此為基礎(chǔ)水手9號[13]和火星勘測號對這種方式進(jìn)行了在軌驗(yàn)證。這種方法需要預(yù)先確定雙曲線軌道漸近線方向，并且對速度信息的估計(jì)能力弱，多作為地面多普勒導(dǎo)航的補(bǔ)充。

因此，本文提出一種基于純天體視線方向測量的火星接近段自主軌道確定算法，觀測信息為一組連續(xù)測量的火衛(wèi)二中心相對探測器的方向矢量，利用最小二乘法，每測量一組火衛(wèi)二視線矢量得到其中某一時刻探測器完整的軌道速度、位置狀態(tài)估計(jì)量。這種方法不依賴于地面先驗(yàn)信息，可獨(dú)立完成軌道確定。又考慮到觀測過程中火衛(wèi)二和火星時而同時出現(xiàn)在同一視場下，此時可以結(jié)合火星中心視線矢量方向以及火衛(wèi)二的星歷，得到精度較高的飛行器位置估計(jì)量，和第一種方法聯(lián)合可以期望提高導(dǎo)航精度。

1 系統(tǒng)模型

1.1 火星接近段軌道動力學(xué)模型

火星接近段軌道為從火星引力球邊沿開始，經(jīng)歷雙曲線軌道運(yùn)動，直到接近到火星表面200 km左右的軌道段。這一階段飛行器動力學(xué)模型建立在火星赤道慣性坐標(biāo)系O-xiyizi下，系統(tǒng)軌道導(dǎo)航的狀態(tài)向量為：

(1)

(2)

式中：μm為火星引力常數(shù)；m為飛行器質(zhì)量；FS為太陽第三體引力攝動加速度；Fε為所有其他攝動加速度，主要包括火星非球形攝動和太陽光壓攝動。

利用EKF對導(dǎo)航方法進(jìn)行濾波時，該系統(tǒng)的敏感矩陣為：

(3)

式中：aS和aε為太陽第三體引力攝動加速度和其他攝動加速度。

1.2 視線矢量的測量

(1)觀測對象的選擇

火星擁有兩顆形狀不規(guī)則的衛(wèi)星，其中火衛(wèi)二Deimos是距離火星較遠(yuǎn)、質(zhì)量較小的一顆，表1列出了火星兩顆衛(wèi)星的基本參數(shù)[4-6]。

表1 火星衛(wèi)星基本參數(shù)Table 1 Physical parameters of Mars moons

可以看出火衛(wèi)二的軌道半徑較大，適合作為測量的基線；火衛(wèi)二和火衛(wèi)一都是不規(guī)則天體，但火衛(wèi)二有平均直徑，而火衛(wèi)一的不規(guī)則程度過高，不能用直徑來衡量；火衛(wèi)二比火衛(wèi)一體積小，是太陽系內(nèi)最小的衛(wèi)星，對其形心的判斷更精確；且火衛(wèi)二、火衛(wèi)一的軌道偏心率和軌道傾角都非常小，可以近似為沒有傾角的圓軌道。這些特性都使得火衛(wèi)二適合作為輔助觀測對象。

(2)視線矢量測量原理

導(dǎo)航相機(jī)對被觀測天體拍照，經(jīng)過圖像處理可提取目標(biāo)形心，結(jié)合星敏感器數(shù)據(jù)測量出該星光在星敏感器測量坐標(biāo)系下的方向矢量s，再通過星敏感器安裝矩陣Cs進(jìn)行坐標(biāo)轉(zhuǎn)化，得到該星光在飛行器本體坐標(biāo)系下的方向矢量sb，即本體坐標(biāo)系下星光相對飛行器的視線矢量。進(jìn)一步通過姿態(tài)矩陣Cib可得到慣性系下星光相對飛行器的視線矢量。

sb=Cs·s

(4)

si=Cib·sb

(5)

通過測量一個或多個天體的視線矢量，并結(jié)合探測器和天體間的幾何關(guān)系以及最優(yōu)估計(jì)方法，可估計(jì)出探測器在慣性系下的位置速度信息，從而實(shí)現(xiàn)軌道確定，具體方法將在下節(jié)詳細(xì)介紹。

2 觀測方法

2.1 連續(xù)觀測火衛(wèi)二

這里介紹的連續(xù)觀測火衛(wèi)二自主導(dǎo)航方法是在單純角度觀測的高斯初軌確定方法[17]基礎(chǔ)上，將觀測視角的位置由地面轉(zhuǎn)移到火星探測器上，火星是引力中心，火衛(wèi)二圍繞火星轉(zhuǎn)動，因只觀測火衛(wèi)二視線信息，故后文簡稱為“單目標(biāo)導(dǎo)航方法”。探測器進(jìn)入火星影響球后沿雙曲線接近軌道飛行，此時探測器、火星、火衛(wèi)二間的位置幾何關(guān)系可抽象成圖1。

圖1 火衛(wèi)二觀測幾何原理Fig.1 Geometric relation of Deimos observation

(6)

這3個矢量方程也是9個標(biāo)量方程，包含12個未知數(shù)，又由于3個矢量r1、r2和r3屬于同一軌道，根據(jù)角動量守恒原理有：

r2=c1r1+c3r3

(7)

根據(jù)二體運(yùn)動方程，利用拉格朗日系數(shù)，可知探測器軌道某一時刻的狀態(tài)矢量可由任意時刻(t2)的狀態(tài)矢量(r2,v2)來表示，故有：

(8)

式中：fk和gk為拉格朗日系數(shù)在tk時刻的值，若3次觀測的時間間隔足夠小，fk和gk僅與r2有關(guān)。這樣得到18個方程解18個未知數(shù)，方程組可解。傳統(tǒng)測定軌所用的高斯解析解法對測量誤差容忍度極低，因此本文利用一次完成最小二乘法，放棄求取以上方程組式(6)～(8)的唯一解，轉(zhuǎn)而增加觀測次數(shù)。當(dāng)解不唯一時，求取最小二乘意義下的最優(yōu)解。

最小二乘格式的系統(tǒng)輸入/輸出關(guān)系可以用下式表示：

z(k)=hT(k)θ+n(k)

(9)

式中：z(k)是待辨識系統(tǒng)(過程)的輸出；hT(k)為觀測向量；n(k)為隨機(jī)噪聲，均值為零；θ為待辨識的參數(shù)。則有以下準(zhǔn)則函數(shù)：

(10)

zL=HLθ+nL

(11)

其中：

(12)

(13)

(14)

式中：na為輸出向量的階次。用nb表示輸入向量的階次，則上述線性方程組式(14)有解的條件是L≥na+nb。而L=na+nb。只有nL=[0,0,…,0]T時，θ有唯一解。但只有當(dāng)L>na+nb，nL≠[0,0,…,0]T時，問題才有意義。方程組有不止一個解，可求取最優(yōu)解，所以為了保證辨識精度，數(shù)據(jù)長度L可以取得充分大。

最小二乘算法分為一次完成算法和遞推算法，其中一次完成算法的解為：

(15)

具體到本節(jié)利用火衛(wèi)二中心視線矢量作為測量量的軌道參數(shù)預(yù)報問題，將式(8)代入式(6)中，有：

(16)

待辨識的參數(shù)為t2時刻位置矢量r2和速度矢量v2以及3個測量時刻火衛(wèi)二到探測器的距離ρ1、ρ2、ρ3。將拉格朗日系數(shù)在t2作泰勒展開，截至Δt的4次方：

(17)

若測量間隔時間Δt足夠短可略去2次以上的項(xiàng)，這樣拉格朗日系數(shù)中只出現(xiàn)未知量r2，而更高次項(xiàng)會含有v2，故這一步計(jì)算存在近似處理。

將(16)式寫成最小二乘格式：

(18)

式中：I3為3階單位陣，tri為第i列為全1列，其余列為全0列的3×3階矩陣，此時對應(yīng)最小二乘各個量有：

(19)

從式(19)可以看出至少3個時刻的測量量才能解出一個θ，數(shù)據(jù)長度為3。為了求出“最優(yōu)”辨識參數(shù)，需要取更大的數(shù)據(jù)長度，即更多的測量量。而增加測量量意味著HL矩陣維數(shù)增加，式(15)中的求逆運(yùn)算計(jì)算量急劇增加，所以數(shù)據(jù)長度亦不能取過大?？紤]到計(jì)算的時間成本，數(shù)據(jù)長度取6～7在可接受范圍。又如式(17)所示，若按截至Δt的2次方求拉格朗日系數(shù)，要求Δt足夠小，若Δt過大會導(dǎo)致系統(tǒng)誤差過大，因此待辨識的速度位置矢量應(yīng)避免選擇處于一組測量值兩頭時刻的值。本文選擇t2時刻的軌道半徑作為待估計(jì)量，則有：

(20)

為了彌補(bǔ)數(shù)據(jù)長度的不足，可以進(jìn)一步求取辨識參數(shù)的均值，具體作法可將以下兩種方法結(jié)合：

1)測量參數(shù)不變，隨機(jī)誤差變化，反復(fù)求取同一組辨識參數(shù)，并取均值。

2)除待辨識參數(shù)所在時刻的所有參數(shù)不變，其他測量時刻每次加1，再次測量，即更新測量量，再辨識同一組參數(shù)。如此進(jìn)行n步，再對辨識出的n組參數(shù)求均值。

圖2 準(zhǔn)則函數(shù)變化規(guī)律曲線Fig.2 Criterion function curve

這一變化規(guī)律與軌道預(yù)報誤差的變化規(guī)律一致，故選取在準(zhǔn)則函數(shù)“拐點(diǎn)”處停止觀測，即可表征估計(jì)誤差最小處。此后只進(jìn)行軌道遞推，此時距近點(diǎn)有足夠的時間制定后續(xù)制動計(jì)劃。

2.2 同時觀測火星和火衛(wèi)二

從火星影響球邊界到軌道確定完成的過程中，火衛(wèi)二和火星有時會出現(xiàn)在同一視場下，尤其是當(dāng)探測器距離火星較遠(yuǎn)時，火衛(wèi)二和火星長時間保持在同一視場內(nèi)。因此這種情況下可以考慮加入火星中心相對探測器的視線矢量作為補(bǔ)充觀測量，以期提高軌道預(yù)報的精度。

同時對火衛(wèi)二和火星拍照可以同時得到方向矢量ρ和r的測量值，方便直接計(jì)算出探測器到火星的距離r，因同時觀測火衛(wèi)二和火星視線信息，故后文簡稱為“雙目標(biāo)導(dǎo)航方法”。根據(jù)火星、火衛(wèi)二以及探測器間的幾何關(guān)系有：

rr=Rd+ρρ

(21)

(22)

根據(jù)以上原理，觀測方程可寫為：

(23)

觀測方程的敏感矩陣為：

(24)

3 聯(lián)合導(dǎo)航策略

以上兩種導(dǎo)航方法在適用范圍、預(yù)報精度等方面存在差異。只觀測火衛(wèi)二的單目標(biāo)導(dǎo)航方法可以得到完整的觀測信息，其導(dǎo)航精度隨著探測器沿雙曲線軌道接近火星呈先提高再降低的趨勢，需在軌計(jì)算停止估計(jì)的時刻；同時觀測火衛(wèi)二和火星的方法位置估計(jì)精度很高，沒有速度的直接觀測信息，速度估計(jì)精度較低。因此，將兩種導(dǎo)航方法結(jié)合起來，根據(jù)馬爾可夫最優(yōu)估計(jì)理論可期得到完整的觀測信息和更高精度的軌道估計(jì)結(jié)果。

具體聯(lián)合導(dǎo)航的流程如圖3所示。

圖3中初始軌道參數(shù)指進(jìn)入火星影響球以后雙曲線軌道的初始值，這個值可由上一階段巡航段軌道的導(dǎo)航數(shù)據(jù)給出，也可由地面測控給出一個誤差較大的初始測定軌道。T1代表火衛(wèi)二觀測時刻，T2代表同時觀測火衛(wèi)二及火星時刻。當(dāng)未達(dá)到T2或者僅能在視場內(nèi)觀測到火衛(wèi)二時,只執(zhí)行連續(xù)觀測火衛(wèi)二的自主導(dǎo)航方法;當(dāng)同時觀測火衛(wèi)二及火星時,結(jié)合該方法預(yù)報的軌道半徑模值和只觀測火衛(wèi)二方法預(yù)報的速度矢量及位置單位矢量，得到更高精度的三自由度軌道速度位置預(yù)報值。

圖3 聯(lián)合導(dǎo)航方法計(jì)算流程Fig.3 Flow chart of combined navigation algorithm

4 仿真校驗(yàn)

以常見火星接近段雙曲線軌道為例，從進(jìn)入火星影響球不久，距離火星約5.8×105km處開始，分別對單目標(biāo)和雙目標(biāo)視線導(dǎo)航以及聯(lián)合導(dǎo)航方法進(jìn)行仿真分析。

4.1 仿真條件

仿真條件如表2所示。

表2 仿真條件Table 2 Simulation condition

需要說明，表中給出的軌道初始位置速度信息因篇幅原因保留小數(shù)點(diǎn)后4位，實(shí)際仿真中保留14位。

此外，單目標(biāo)測量間隔時間為2 400 s，雙目標(biāo)測量時間間隔100 s，為節(jié)省計(jì)算量，聯(lián)合導(dǎo)航算法中雙目標(biāo)測量時間間隔放寬至1 000 s。

4.2 仿真結(jié)果

圖4、圖5是連續(xù)單視線矢量測量導(dǎo)航方法的位置、速度預(yù)報誤差曲線，仿真中根據(jù)準(zhǔn)則函數(shù)變化規(guī)律大概于近點(diǎn)前19 h停止觀測，僅進(jìn)行軌道遞推。

圖4 單目標(biāo)觀測軌道位置預(yù)報誤差曲線Fig.4 Altitude error estimation of Deimos-only method

可以看出，隨著探測器接近火星，該方法位置預(yù)報誤差最初階段經(jīng)歷了大幅度震蕩，甚至大于初始軌道誤差，而后迅速收斂，在距近點(diǎn)33 h時位置估計(jì)誤差大小已明顯小于初始誤差，最終位置預(yù)報誤差逐漸收斂至20 km以內(nèi)，在近點(diǎn)處誤差稍有增大。該方法的速度預(yù)報誤差逐漸收斂，除近點(diǎn)跳變區(qū)域外，速度預(yù)報精度可以達(dá)到0.2 m/s以內(nèi)，近點(diǎn)處預(yù)報精度在5 m/s以內(nèi)。

圖6、圖7分別是雙視線矢量測量導(dǎo)航的軌道位置、速度預(yù)報誤差曲線。

圖6 雙目標(biāo)觀測軌道位置預(yù)報誤差曲線Fig.6 Altitude error estimation of Two-target observation method

圖7 雙目標(biāo)觀測軌道速度預(yù)報誤差曲線Fig.7 Velocity error estimation of Two-target observation method

圖6中呈現(xiàn)脈沖較大的區(qū)域是火衛(wèi)二和火星在導(dǎo)航相機(jī)中成像較接近的位置，二者間角度很小，測量誤差影響較大，可設(shè)計(jì)導(dǎo)航流程忽略這部分測量信息，只進(jìn)行軌道遞推，此處保留該區(qū)域是為說明這一特性。其余區(qū)域的位置預(yù)報精度良好，誤差不到10 km，并隨接近火星而逐漸收斂。

根據(jù)第2.2小節(jié)的分析，可以看出由于缺乏速度觀測信息，速度估計(jì)誤差先是緩慢累積，在近點(diǎn)位置誤差累積速度急劇增大，速度預(yù)報精度不理想。

加入火星視線矢量這一補(bǔ)充觀測量后聯(lián)合導(dǎo)航位置預(yù)報誤差曲線如圖8所示。

圖8 聯(lián)合導(dǎo)航方法軌道位置速度預(yù)報誤差曲線Fig.8 Altitude & velocity error estimation of combined navigation algorithm

聯(lián)合后軌道位置預(yù)報精度進(jìn)一步提高，即便在近點(diǎn)脈沖區(qū)域內(nèi)，位置導(dǎo)航精度也達(dá)到了3 km左右；速度預(yù)報誤差也在單目標(biāo)導(dǎo)航方法的導(dǎo)航結(jié)果上進(jìn)一步提高，近點(diǎn)跳變區(qū)域精度達(dá)到1 m/s以內(nèi)，平緩區(qū)域誤差小于0.2 m/s。

進(jìn)一步對火星近點(diǎn)軌道要素信息進(jìn)行預(yù)報，主要仿真了和后續(xù)軌道機(jī)動任務(wù)息息相關(guān)的近點(diǎn)軌道半徑預(yù)報誤差和近點(diǎn)軌道速度預(yù)報誤差，如圖9、圖10所示。

圖9 聯(lián)合導(dǎo)航方法近點(diǎn)半徑預(yù)報誤差曲線Fig.9 Perigee altitude error estimation of combined navigation algorithm

圖10 聯(lián)合導(dǎo)航方法近點(diǎn)速度預(yù)報誤差曲線Fig.10 Perigee velocity error estimation of combined navigation algorithm

近點(diǎn)軌道半徑預(yù)報誤差等價于近點(diǎn)軌道高度預(yù)報誤差，可以看出利用初始軌道信息估計(jì)的近點(diǎn)半徑誤差接近1 000 km，不能滿足任何軌道捕獲任務(wù)需求。經(jīng)過導(dǎo)航預(yù)報后，近點(diǎn)高度的預(yù)報精度達(dá)到了2 km以內(nèi)，近點(diǎn)速度預(yù)報精度收斂至0.8 m/s以內(nèi)，滿足軌道捕獲任務(wù)需求。

5 結(jié)束語

通過理論分析與仿真校驗(yàn)可知，文中提出的基于火衛(wèi)二視線矢量測量的火星接近段自主導(dǎo)航方法能夠通過單純依賴視線矢量測量，在缺少速度直接測量信息的情況下達(dá)到接近段導(dǎo)航后續(xù)任務(wù)所需的導(dǎo)航精度。其中火衛(wèi)二單視線測量法可以實(shí)現(xiàn)通過一組視線矢量測量得到探測器軌道速度信息，而雙視線法可以輔助提高狀態(tài)量的估計(jì)精度。聯(lián)合導(dǎo)航的最終結(jié)果顯示，近點(diǎn)軌道高度預(yù)報精度達(dá)到了2 km，近點(diǎn)軌道速度預(yù)報精度達(dá)到了0.8 m/s。

大天體接近段軌道前后速度特性變化劇烈，非線性較強(qiáng)，難以通過單一觀測手段實(shí)現(xiàn)高精度的自主導(dǎo)航目標(biāo)。本文所討論的基于單純視線矢量測量的自主導(dǎo)航方法在接近火星后因軌道動力學(xué)非線性程度增加，在近點(diǎn)附近預(yù)報精度有所下降。在后續(xù)的研究中，可嘗試對導(dǎo)航相機(jī)拍攝圖像進(jìn)行更多地信息挖掘，利用如火星日照面積等信息進(jìn)一步提高軌道預(yù)報精度。