付亮勇, 胡錦昌,2*, 劉一武, 鄧家權(quán)
1. 北京控制工程研究所,北京 100094 2. 空間智能控制技術(shù)國家級重點(diǎn)實(shí)驗室,北京 100094 3. 63798部隊,西昌 615606
近年來,小行星探測已成為主要航天國家深空探索領(lǐng)域的重點(diǎn)發(fā)展目標(biāo)之一[1].由于小行星具有密度分布不均勻、形狀不規(guī)則和運(yùn)動狀態(tài)復(fù)雜等特點(diǎn),導(dǎo)致小行星的引力場極為不規(guī)則和復(fù)雜;此外,探測器在小行星附近受到太陽中心引力、快速自旋的非球形攝動力以及太陽光壓攝動力等作用,使其在近小行星處軌道動力學(xué)復(fù)雜多變[2-3].不規(guī)則小行星引力場確定問題是小行星探測任務(wù)中下降著陸和采樣返回階段順利實(shí)施的關(guān)鍵難點(diǎn)之一.
小行星的引力場需要采用合適的引力模型來描述,常用的引力建模方法主要可分為級數(shù)逼近法和三維模型逼近法[4].級數(shù)逼近法主要通過無窮級數(shù)來逼近引力勢能,可以得到解析形式,主要有球諧函數(shù)模型法[5]和橢球諧函數(shù)模型法[6];球諧系數(shù)可通過地面測定軌數(shù)據(jù)和星上自主導(dǎo)航算法解算得到,對于還未進(jìn)行探測的小行星來說無法事先獲取其準(zhǔn)確的引力信息[7].三維模型逼近法主要采用簡化的三維模型來逼近小行星不規(guī)則形狀,然后通過數(shù)值計算方法得到引力勢能,進(jìn)而完成引力建模,主要有三軸橢球體模型法[8]和多面體模型法[9];三維模型逼近法可利用天文觀測或掠飛拍攝圖像信息以完成模型建立,可以在任務(wù)發(fā)射前獲取目標(biāo)天體的近似引力模型.
由小行星引力場發(fā)展來看,目前對小行星引力場研究以測定軌數(shù)據(jù)反演為主,而利用三維模型正演為輔;如何將兩者數(shù)據(jù)結(jié)合起來,建立精確且便于軌道分析的不規(guī)則引力場模型是將來研究重點(diǎn).文獻(xiàn)[10]提出一種聯(lián)合導(dǎo)航數(shù)據(jù)和小行星形狀模型建立精確引力場模型的方法,有效解決了密度分布不均勻難題;文獻(xiàn)[11-12]提出利用小行星形狀模型和經(jīng)初步定軌數(shù)據(jù)反演解算的引力場模型對小行星進(jìn)行密度分布估計,然后從估計的密度分布中重新構(gòu)建表面引力場,其精度顯著提高.
目前,公開文獻(xiàn)中,對于較大的目標(biāo)小行星,引力場測量大多是采用穩(wěn)定繞飛軌道,通過地面測定軌和星上相對測量數(shù)據(jù)相結(jié)合的方式進(jìn)行反演.對于弱引力的小行星,則一般適合采用近似直線上下自由運(yùn)動的方式,如隼鳥2號的探測目標(biāo)“龍宮”.對于上下運(yùn)動的測量方式,其難點(diǎn)在于自由運(yùn)動弧段一般較短,難以通過積累長時間測量數(shù)據(jù)以提高反演精度,而對于直徑為百米量級的小行星[13]來說,其特殊性困難在于:引力更弱,其準(zhǔn)確測量難度更大;自由運(yùn)動弧段一般只有數(shù)個小時,數(shù)據(jù)積累長度更短;地面定軌精度不足以支撐引力場測量,需依賴星上相對測量實(shí)現(xiàn);在引力場測量過程中,一般需要同時辨識小行星的自旋角速度大小和方向等參數(shù).針對上述問題,提出一種通過視覺導(dǎo)航方式獲取相應(yīng)的量測信息,進(jìn)而依靠星上自主完成引力場測量的方法;還仿真分析了小行星整體質(zhì)量分布差異下的不同飛行軌跡的反演結(jié)果,進(jìn)而提出根據(jù)整體質(zhì)量分布來規(guī)劃引力測量中探測器的飛行軌跡.
通過定軌數(shù)據(jù)進(jìn)行小行星引力場反演,首先需要建立高精度的軌道動力學(xué)模型,以模擬探測器的軌道真值和推導(dǎo)反演過程所需的系統(tǒng)狀態(tài)矩陣.小行星附近軌道動力學(xué)復(fù)雜,受多種攝動力作用,包括太陽中心引力、小行星中心引力及其快速自旋的不規(guī)則形狀攝動力、大行星引力攝動和太陽光壓攝動等作用[14].在模擬探測器的軌道真值時,系統(tǒng)軌道動力學(xué)模型應(yīng)建立在J2000日心黃道慣性坐標(biāo)系下.
探測器軌道動力學(xué)系統(tǒng)如圖1所示.以太陽為中心引力體的探測器軌道動力學(xué)方程為
圖1 探測器軌道動力學(xué)系統(tǒng)Fig.1 Probe orbit dynamics system
(1)
(2)
式中:μ是小行星引力常數(shù);ra是小行星相對太陽的位置矢量;aE是小行星不規(guī)則引力場去除中心引力部分的非球形攝動加速度,需要注意的是,非球形攝動加速度是在小行星固聯(lián)坐標(biāo)系下定義的.
假設(shè)太陽光壓力對小行星的軌道影響可以忽略,以太陽為中心引力體的小行星軌道動力學(xué)方程為
(3)
對于小行星不規(guī)則引力場,通常可以使用下面的引力能勢函數(shù)表示
(4)
式中:RE為小行星參考半徑,參考量的作用只是用來無量綱化球諧系數(shù),并不影響球諧函數(shù)模型的收斂條件;Cnm、Snm為球諧系數(shù),且有C00=1,Sn0=0;n為階數(shù),m為次數(shù);Vnm、Wnm具體表達(dá)式如下:
(5)
式中:r為檢驗點(diǎn)到中心天體質(zhì)心的距離;Pnm(sinφ)是sinφ的締合勒讓德多項式[15];φ、λ分別是地心緯度和地理經(jīng)度,即探測器受到的不規(guī)則引力攝動是在中心天體固聯(lián)坐標(biāo)系下描述的.
為方便計算探測器所受不規(guī)則引力攝動加速度的解析式,下面給出式(5)的遞推公式[16].
初始條件
遞推公式
(6)
式中,x、y、z分別為檢驗點(diǎn)在小行星固聯(lián)坐標(biāo)系下各坐標(biāo)軸的分量,將式(6)代入式(4)便可得到遞推形式的小行星不規(guī)則引力場的球諧模型.
式(4)中球諧系數(shù)是一組用來描述小行星質(zhì)量和形狀分布的參數(shù),C00=1表示一均質(zhì)球體的引力勢,而后面各級數(shù)項表示對均質(zhì)球體引力勢的修正.
為了方便本文對各球諧系數(shù)和引力常數(shù)的估計,需要對式(4)進(jìn)行變形,下列還是采用原來符號,不加以區(qū)分,但需要注意其中區(qū)別;引力勢能模型重寫為
(7)
式中,C00≠1,并且引力常數(shù)μ=C00·RE,式(7)中的球諧系數(shù)跟式(4)中的相差固定的倍數(shù),下文出現(xiàn)的球諧系數(shù)均是按式(7)約定的,這樣將所有待估計的球諧系數(shù)確定好后,引力常數(shù)也就確定好了,上式中球諧系數(shù)即反映了引力場的全信息.
需要注意的是,當(dāng)描述引力場模型的固聯(lián)坐標(biāo)系原點(diǎn)和小行星質(zhì)心重合時,建立在固聯(lián)系下的球諧系數(shù)一階項C10、C11、S11均為0;并且若固聯(lián)坐標(biāo)系的z軸與小行星慣量主軸方向一致時,球諧系數(shù)項C21、S21均為0[17].
為簡化分析反演算法,文章是在小行星形心固聯(lián)坐標(biāo)系和實(shí)際的小行星固聯(lián)坐標(biāo)系重合下研究的,下文統(tǒng)稱為小行星固聯(lián)坐標(biāo)系.
在實(shí)際任務(wù)中,應(yīng)在進(jìn)行引力場測量前,完成對目標(biāo)小行星三維建模以獲取特征點(diǎn)的三維位置,進(jìn)而構(gòu)建導(dǎo)航特征庫[18-19].利用視覺導(dǎo)航獲取的圖像信息來完成對小行星引力場測量,其觀測量的定義為:導(dǎo)航特征點(diǎn)在小行星固聯(lián)坐標(biāo)系下的三維位置;導(dǎo)航特征點(diǎn)在相機(jī)成像相平面上的像素坐標(biāo),即像元p和像線l坐標(biāo)值,其表達(dá)式為
(8)
圖2 導(dǎo)航相機(jī)成像模型Fig.2 Navigation camera imaging model
(9)
設(shè)慣性系到小行星固聯(lián)系的坐標(biāo)系按照313順序進(jìn)行轉(zhuǎn)動,則式(9)中方向余弦矩陣CFI的具體表達(dá)式為
CFI=C3,ωaCFI0=C3,ωaC3,γC1,βC3,α
(10)
其中:ωa為小行星自旋角速度大小;下標(biāo)i=1,3表示旋轉(zhuǎn)軸,α、β、γ分別表示繞各軸的轉(zhuǎn)角,其中α、β代表了自旋軸的方向;CFI0表示初始時刻小行星固聯(lián)系相對于慣性系的方向余弦矩陣.
系統(tǒng)測量方程
(11)
(12)
其中,δij為狄拉克函數(shù).
假設(shè)待估計的初始狀態(tài)為x0,則測量敏感矩陣為
(13)
為推導(dǎo)反演算法所需要的系統(tǒng)狀態(tài)轉(zhuǎn)移矩陣,首先需要求取探測器在小行星慣性系下的軌道動力學(xué)方程;用式(1)減去(3)可得
(14)
考慮到在實(shí)際引力測量階段,探測器離小行星最遠(yuǎn)距離不超過幾百米,而小行星相距太陽在1 AU左右,可算出式(14)中括號項部分所引起的加速度攝動量級遠(yuǎn)小于受小行星不規(guī)則引力攝動量級,故可經(jīng)簡化消去[20];并且將小行星不規(guī)則攝動加速度aAST用式(6)~(7)替換,經(jīng)整理,最終可得到探測器在小行星慣性系下的軌道動力學(xué)方程為
(15)
式中,v為探測器在小行星慣性系下的速度矢量,a為其他未建模誤差項,rF為探測器在小行星固聯(lián)系下的位置矢量,U為包括各階攝動項的引力勢能函數(shù).
系統(tǒng)的狀態(tài)方程通常由探測器的相對軌道動力學(xué)方程、姿態(tài)運(yùn)動學(xué)方程和姿態(tài)動力學(xué)方程組成.而在引力場測量階段,姿態(tài)可以認(rèn)為是直接測量得到的已知量,因此狀態(tài)方程可以不考慮姿態(tài)運(yùn)動學(xué)和動力學(xué);并且假設(shè)太陽光壓是完全已知的,其模型誤差可忽略,由于觀測時間只有數(shù)小時,忽略其它未建模攝動力的影響.此外,因小行星自旋角速度大小和方向難以事先獲得準(zhǔn)確值,需要作為未知狀態(tài)量一同估計.
由此,系統(tǒng)的狀態(tài)變量X定義為
X=[rT,vT,Cnm,…,Snm,α,β,λ,ωa]T
(16)
因此系統(tǒng)的狀態(tài)方程f(X)表達(dá)式為
(17)
將探測器初始時刻的狀態(tài)預(yù)測值(先驗值)當(dāng)作初值,通過積分得到預(yù)估的探測器軌道,由此得到特征點(diǎn)的像元像線預(yù)估值,然后和光學(xué)相機(jī)實(shí)際測量的特征點(diǎn)像元像線值比較,并根據(jù)最小二乘批處理算法得到初始狀態(tài)的修正值;此過程需要經(jīng)過多次迭代,不斷修正初始先驗?zāi)P蛥?shù)和其它待估狀態(tài)初值,以此得到最終符合收斂要求的引力模型.為了實(shí)現(xiàn)對初始狀態(tài)的更新,需要求取狀態(tài)轉(zhuǎn)移矩陣,為此需要先對狀態(tài)方程線性化,即求取系統(tǒng)狀態(tài)方程f(X)對系統(tǒng)狀態(tài)的偏導(dǎo)數(shù)矩陣,如下所示:
(18)
由此可得到離散化后的狀態(tài)轉(zhuǎn)移矩陣
(19)
對式(13)的敏感矩陣進(jìn)行求解,令系統(tǒng)待估計的初始狀態(tài)x0為
(20)
將式(9)代入式(13),并根據(jù)式(19)可得到系統(tǒng)量測模型的第i的特征點(diǎn)的觀測矩陣為
(21)
其中,Φ(t,t0)|3×[(n+1)2+6]為取狀態(tài)轉(zhuǎn)移矩陣的前3行,xa,0=[α,β,λ,ωa],n是待求取的引力模型階數(shù).
考慮到觀測方程的強(qiáng)非線性和過多估計參數(shù),并且待估參數(shù)量級較小,若在最小二乘迭代過程中,狀態(tài)增量過大會容易導(dǎo)致不收斂;對此可以利用待估參數(shù)的先驗信息來將其盡量保持在收斂范圍內(nèi).因此可將最小二乘殘差意義下的代價函數(shù)擴(kuò)展為[10]
(22)
利用最小二乘批處理算法估計初始狀態(tài)量x0的原則是:通過合理的選擇近似修正項δx0,使線性化的預(yù)測殘差的加權(quán)平方和J(x0)最小.
(23)
式中:Jk、Jk+1分別為上一次和下一次的代價函數(shù)迭代值;η為事先指定的小量;對于迭代終止條件還可以采用,當(dāng)δx0小于某個值等.
理論上通過無窮階球諧函數(shù)模型可以逼近任意形狀的小行星在其布里淵球體外的引力勢,但是由于無窮階本身是無法計算的,通常取其重要項,因此存在截斷誤差,在本文主要分析其四階以內(nèi)的性質(zhì).取余弦面球諧函數(shù)Pnm(sinφ)cos(mλ),繪制其表示的引力勢分布情況如圖3所示,對于正弦面球諧函數(shù)Pnm(sinφ)sin(mλ)而言,它與余弦面球諧函數(shù)只是在經(jīng)度上相位相差90°,前者繞極軸旋轉(zhuǎn)90(°)/m即可得后者,情況本質(zhì)上一致.
圖3 四階余弦面球諧函數(shù)引力勢分布示意圖Fig.3 Gravity potential distribution of four-degree cosine surface spherical harmonic function
圖3中展現(xiàn)了四階余弦面球諧函數(shù)的引力勢情況.第一個是均質(zhì)球體的引力勢分布,余下的第一列是帶諧項,對角線上是扇諧項,其余的是田諧項;暗的地方引力勢為負(fù),表示此處質(zhì)量相比于均質(zhì)球體來說不足;亮的地方引力勢為正,表示此處質(zhì)量相對均質(zhì)球體來說是聚集的.
由上分析可知,球諧系數(shù)反映了小行星的整體質(zhì)量分布情況,因此可以根據(jù)其整體的質(zhì)量分布來規(guī)劃探測器相對于小行星的飛行軌跡,以此來反演出球諧模型.對于質(zhì)量主要沿緯度方向變化的采用南北豎直軌跡;對于質(zhì)量主要沿經(jīng)度方向呈扇形分布的采用東西橫向軌跡;但對于質(zhì)量主要沿經(jīng)緯度方向呈凹凸田塊狀分布的無法直接給出探測器的飛行軌跡,從直觀上看,存在豎直、橫向以及傾斜飛行軌跡,實(shí)際中可能需要用不同方向的軌跡或者相同方向但相對于小行星不同位置的軌跡下量測值組合反演,需要進(jìn)一步研究.
雖然不規(guī)則小行星的密度不均勻,無法獲取其準(zhǔn)確的質(zhì)量分布,但是密度在一定范圍內(nèi)變化,其形狀大致可以決定其整體質(zhì)量的分布,因此可根據(jù)三維建模獲取的圖像信息大致了解其整體質(zhì)量分布,進(jìn)而規(guī)劃探測器飛行軌跡以此進(jìn)行引力場測量,該方法并不依賴于目標(biāo)天體精確的質(zhì)量分布信息.
設(shè)直徑為百米量級的某小行星引力常數(shù)為0.115 m3/s2,參考半徑為100 m,均質(zhì)半徑為50 m;其四階以內(nèi)球諧系數(shù)參考433Eros并加以調(diào)整,上述3種典型質(zhì)量分布差異情況下的球諧系數(shù)設(shè)置如表1所示;其對應(yīng)的探測器飛行軌跡如圖4~圖6所示.
圖4 質(zhì)量呈帶狀分布Fig.4 Mass distribution in bands
表1 整體質(zhì)量分布差異下球諧系數(shù)Tab.1 Spherical harmonic coefficients under different mass distributions
后續(xù)將針對這幾種情況進(jìn)行具體的仿真校驗,驗證在引力場測量中以小行星整體質(zhì)量分布為依據(jù)來規(guī)劃探測器飛行軌跡的合理性.
本節(jié)首先通過仿真分析典型質(zhì)量差異下的引力場測量情況,來對上述提出的探測器軌道規(guī)劃方案進(jìn)行校驗;然后在此基礎(chǔ)上對直徑為百米量級的目標(biāo)小行星進(jìn)行引力場測量,由于其實(shí)際球諧系數(shù)未知,本文參考433Eros的球諧系數(shù)來設(shè)定.除上文中已經(jīng)合理假設(shè)參數(shù)外,設(shè)目標(biāo)小行星自旋周期為30 min,實(shí)際誤差為5%;方向角α、β、γ分別為30°、60°、10°,實(shí)際誤差各5°.
考慮相機(jī)視角條件,探測器的速度脈沖和初始位置在小行星瞬時軌道坐標(biāo)系下給定;此外仿真中導(dǎo)航特征庫是根據(jù)模擬球體和三維建模誤差構(gòu)建出來的;導(dǎo)航敏感器參數(shù)如表2所示.
表2 導(dǎo)航敏感器參數(shù)Tab.2 Parameters of navigation sensor
3種質(zhì)量分布下的球諧系數(shù)標(biāo)稱值已在表1給出;其先驗值和標(biāo)準(zhǔn)差如表3所示.在距離小行星表面30 m位置左右施加速度脈沖,使探測器遠(yuǎn)離小行星,在小行星引力場作用下自由運(yùn)動,各質(zhì)量分布下的探測器軌跡采用圖4、圖5和圖6中設(shè)計的軌跡;仿真時間為2 500 s,特征點(diǎn)提取間隔為30 s左右.
圖5 質(zhì)量呈扇形分布Fig.5 Mass distribution in sectorial shape
圖6 質(zhì)量呈田塊分布Fig.6 Mass distribution in non-strict symmetry
表3 球諧系數(shù)先驗值和標(biāo)稱值Tab.3 Prior and nominal values
由圖7可知,對于整體質(zhì)量沿緯度方向呈帶狀分布的小行星,在引力場測量時,采取南北豎直向軌跡,其最大誤差不超過10%,而對于橫向軌跡,最大誤差超過了20%,由此可見探測器的相對飛行軌跡采取南北豎直向較好.
圖7 質(zhì)量呈帶狀分布下引力測量誤差Fig.7 Measurement errors under banded distribution
由圖8可知,對于整體質(zhì)量沿經(jīng)度方向呈扇形分布的小行星,在引力場測量時,采取東西橫向軌跡,其最大誤差不超過10%,而對于豎直向軌跡,最大誤差可達(dá)到20%,由此可見探測器的相對飛行軌跡采取東西橫向較好.
圖8 質(zhì)量呈扇形分布下引力測量誤差Fig.8 Measurement errors under sectorial distribution
由圖9可知,無論哪種飛行軌跡,其大部分區(qū)域測量誤差都在10%以上;由此可見對于整體質(zhì)量沿經(jīng)緯度方向呈凹凸田塊狀分布的小行星,由于其質(zhì)量在經(jīng)度和緯度方向均存在聚集和缺失,導(dǎo)致單一方向下一條探測器飛行軌跡的測量數(shù)據(jù)無法準(zhǔn)確解算其全局引力信息,應(yīng)當(dāng)組合多條不同軌跡下的定軌數(shù)據(jù)進(jìn)行反演.圖10是在上述橫向軌跡反演的基礎(chǔ)上,將其當(dāng)作先驗值,再次利用一條中緯度附近起始的傾斜軌跡下量測數(shù)據(jù)進(jìn)行迭代優(yōu)化結(jié)果.
由圖10可知,經(jīng)不同軌跡下量測數(shù)據(jù)組合反演后,其引力場測量精度明顯提高,絕大部分區(qū)域測量誤差不超過10%;圖9和圖10仿真結(jié)果對比表明,對于質(zhì)量并不只沿某一方向分布的小行星來說,需要組合不同相對位置下的定軌數(shù)據(jù)以充分反演其引力信息.
圖9 質(zhì)量呈田塊分布下引力測量誤差Fig.9 Measurement errors under asymmetric distribution
圖10 組合軌跡下測量誤差Fig.10 Measurement errors under combined trajectory
利用上述仿真校驗后的結(jié)論,來規(guī)劃相應(yīng)的飛行軌跡對一個四階球諧引力場反演.如前所述,目標(biāo)小行星的四階球諧系數(shù)參考433Eros設(shè)定,但引力常數(shù)和其它參數(shù)不變;狀態(tài)參數(shù)的標(biāo)稱值、先驗值和標(biāo)準(zhǔn)差如表4所示;同樣在距離小行星表面30 m左右處給探測器施加速度脈沖,飛行軌跡仿真時間為2 500 s,特征點(diǎn)提取間隔為30 s左右.
根據(jù)上述球諧系數(shù)可知,若通過均質(zhì)球體展示,其整體質(zhì)量應(yīng)是沿經(jīng)度方向變化,對比實(shí)際的433Eros外形,符合這一理想化結(jié)果.以此為標(biāo)稱值的均質(zhì)不規(guī)則球體,其整體質(zhì)量雖沿經(jīng)度方向變化,但緯度方向也存在少量質(zhì)量聚集和缺失;因此可以根據(jù)上述軌跡規(guī)劃方案,選取探測器橫向方向軌跡來進(jìn)行其整體引力測量,此外在高緯度上添加另一條橫向飛行軌跡,以此來優(yōu)化第一條飛行軌跡下的反演結(jié)果.兩條飛行軌跡如圖11所示;四階球諧系數(shù)和引力常數(shù)反演結(jié)果如表4所示;引力場的整體測量誤差如圖12所示.
表4 四階球諧系數(shù)標(biāo)稱值、標(biāo)準(zhǔn)差、先驗值及估計值Tab.4 Nominal values, priori values, estimated values and standard deviation of the four-degree spherical harmonic coefficients
圖11 探測器飛行軌跡設(shè)計Fig.11 Flight trajectory planning for the probe
圖12 引力場測量誤差分布Fig.12 The distribution of gravity field measurement errors
從表4的估計結(jié)果可知,對于小行星引力常數(shù)的估計誤差可達(dá)到0.2%以內(nèi);從圖12可知,其整體引力測量誤差在5%以內(nèi);可見該測量方法對弱引力小行星的不規(guī)則引力場測量是可行的,精度較高.
由實(shí)驗結(jié)果可知,對于無繞飛軌道的小半徑弱引力目標(biāo)天體,根據(jù)其整體質(zhì)量分布來規(guī)劃探測器的飛行軌跡;然后利用自主導(dǎo)航數(shù)據(jù)可以有效解算其引力常數(shù)和全局引力信息.
針對小行星探測任務(wù)中不規(guī)則弱引力場測量難題,提出一種利用視覺導(dǎo)航數(shù)據(jù)來對其引力場球諧模型反演的方法.在探測器飛行過程中,利用光學(xué)導(dǎo)航相機(jī)獲取特征點(diǎn)的像素信息,結(jié)合其從導(dǎo)航特征庫中匹配得到的三維位置坐標(biāo),然后通過加權(quán)最小二乘批處理算法處理飛行軌跡下的積累數(shù)據(jù)來估計小行星引力常數(shù)和球諧系數(shù);此外討論了如何根據(jù)小行星整體質(zhì)量分布來規(guī)劃探測器的飛行軌跡,以利用較少的飛行軌跡下量測數(shù)據(jù)來完成目標(biāo)天體的全局引力信息測量.仿真結(jié)果表明本方法適用于無繞飛軌道的小半徑弱引力目標(biāo)天體,且適用于密度分布非均勻天體;在充分考慮各種誤差條件下,具有較高的測量精度和一定的工程實(shí)用性.