李璟璟, 張玉兔, 王文從, 胡慧君, 邵思霈
(山東航天電子技術(shù)研究所, 山東 煙臺(tái) 264000)
自20世紀(jì)60年代X射線源天體被發(fā)現(xiàn)以來,基于X射線源的天文導(dǎo)航技術(shù)成為實(shí)現(xiàn)航天器自主導(dǎo)航的又一研究熱點(diǎn)。2016年11月,我國首次發(fā)射了以X射線導(dǎo)航為目的的X射線脈沖星探測衛(wèi)星,獲取了寶貴的在軌數(shù)據(jù)。2017年6月,長四乙發(fā)射的“慧眼”衛(wèi)星,也搭載了X射線探測載荷,通過在軌獲取X射線數(shù)據(jù),同樣能夠?yàn)閄射線導(dǎo)航提供有力支持。X射線脈沖星導(dǎo)航方法主要采用脈沖到達(dá)時(shí)間(time of arrival, TOA)觀測的手段,該方法具有完全自主、全信息、高精度等優(yōu)點(diǎn),但是也存在脈沖折疊時(shí)間長、探測器面積大、TOA估計(jì)精度受限等缺點(diǎn)。此外,要實(shí)現(xiàn)完全的自主導(dǎo)航,還需要通過相當(dāng)長時(shí)間的脈沖星觀測,建立脈沖星數(shù)據(jù)庫[1-4]。相比于TOA觀測技術(shù),X射線成像技術(shù)不需要光子累積時(shí)間,探測器面積小。此外,與傳統(tǒng)光學(xué)波段敏感器相比,X射線成像探測器成本低、重量輕,星圖識(shí)別速度快,既可以作為現(xiàn)有導(dǎo)航手段的補(bǔ)充,也可獨(dú)立提供高精度的導(dǎo)航信息[5-7]。
結(jié)合X射線新型導(dǎo)航手段的優(yōu)點(diǎn),將X射線探測器應(yīng)用于自主天文導(dǎo)航[8]。提出基于X射線源/地心矢量觀測的航天器自主天文導(dǎo)航方法。建立了自主導(dǎo)航觀測模型,并針對系統(tǒng)觀測殘差進(jìn)行了殘差分析和建模,利用UKF濾波器,實(shí)現(xiàn)了航天器自主導(dǎo)航信息的解算。
基于X射線源/地心矢量觀測的自主導(dǎo)航原理如圖1所示。依據(jù)衛(wèi)星、X射線導(dǎo)航星以及地球三者的幾何關(guān)系,可以得到的角度觀測量有星光角距α、地球視角β以及星光仰角γ。星光角距是從衛(wèi)星上觀測到的導(dǎo)航X射線源的星光矢量方向與地心矢量方向之間的夾角。地球視角是從衛(wèi)星上觀測地球邊緣構(gòu)成的圓錐角。星光仰角是從衛(wèi)星上觀測的導(dǎo)航X射線源與地球邊緣的切線方向之間的夾角。結(jié)合衛(wèi)星軌道動(dòng)力學(xué)模型,并運(yùn)用最優(yōu)狀態(tài)估計(jì)方法,可以得到航天器自主導(dǎo)航信息。
圖1 X射線源/地心矢量觀測自主導(dǎo)航原理Fig.1 Autonomous navigation based on X-ray/geocentricvector observation
基于X射線源/地心矢量觀測的航天器自主導(dǎo)航系統(tǒng)主要由探測器單元和導(dǎo)航計(jì)算單元構(gòu)成。其中探測器單元包括X射線探測器和地球敏感器,用于獲取二者矢量信息。星載導(dǎo)航計(jì)算單元包括觀測量合成模塊、星歷數(shù)據(jù)庫和最優(yōu)濾波模塊。系統(tǒng)基本方案如圖2所示。
圖2 自主導(dǎo)航系統(tǒng)方案Fig.2 Autonomous navigation system scheme
兩種探測器方案可以實(shí)現(xiàn)X射線源矢量方向的測量。第1種是編碼孔徑盤像素探測器。通過對天空X射線源進(jìn)行成像,獲取視野中的X射線源特征、輪廓等信息,進(jìn)而獲取X射線源在慣性坐標(biāo)系下的矢量方向。研究表明,利用有效面積為650 cm2的X射線探測器觀測SC0-X1射線源,只需1 ms便可以獲取可識(shí)別的圖像[5]。第2種方案類似于傳統(tǒng)的星掃描儀。通過準(zhǔn)直器和機(jī)械掃描的方法掃描已知X射線源,通過一定的搜索算法確定X射線源單位矢量方向[9]。
以X射線源為觀測對象,通過星圖識(shí)別,獲得X射線源的單位方向矢量。其方位信息來自像平面坐標(biāo)系中的X射線源矢量Sp。設(shè)(xc,yc)為像平面坐標(biāo)系中的像點(diǎn)坐標(biāo),fc為成像焦距,則
(1)
(2)
通過地球敏感器獲取地平信息,并結(jié)合地球物理參數(shù)信息,計(jì)算得到地心在地球敏感器基準(zhǔn)坐標(biāo)系中的矢量。利用紅外地球敏感器輸出的地平掃入脈沖、地平掃出脈沖以及磁基準(zhǔn)脈沖。通過幾何關(guān)系,可以得到地球敏感器坐標(biāo)系中地心方向的矢量為
(3)
衛(wèi)星通過X射線源探測器觀測得到X射線源方向矢量,通過地球敏感器敏感地平并計(jì)算獲得地心矢量。從而進(jìn)一步合成得到關(guān)于兩者夾角的實(shí)際觀測量。為了進(jìn)行狀態(tài)估計(jì),還需要對該實(shí)際觀測量進(jìn)行數(shù)學(xué)上的表達(dá)。由圖1所示的幾何關(guān)系可知,X射線源方向矢量和地心矢量的夾角α為位置的函數(shù),可建立觀測方程為
(4)
式中,rc為J2000地心慣性系下衛(wèi)星的位置矢量;s為衛(wèi)星到X射線源的單位矢量;va為測量噪聲,通常為高斯白噪聲。
如果以星光仰角為觀測矢量,則可直接由地球敏感器獲得地平方向矢量,與星光方向矢量共同構(gòu)建觀測方程
(5)
式中,Re為地球半徑。
將系統(tǒng)觀測殘差建模為加性誤差,則式(4)~式(5)表示的觀測方程可以統(tǒng)一寫為以下形式,即
Yob=Y+ΔY+vob
(6)
式中,Yob為實(shí)際觀測量;Y為真實(shí)量;ΔY為X射線源方向矢量觀測誤差和地球敏感器方向矢量觀測誤差引起的觀測量殘差的和;vob為測量噪聲,服從零均值高斯分布。
根據(jù)式(4),令
M=-(rc·s)/rc=-r·s
(7)
式中,r為探測器在地心慣性坐標(biāo)系中的單位位置矢量。由X射線源方向矢量觀測誤差和地球敏感器方向矢量觀測誤差引起的點(diǎn)積偏差可表示為ΔM,有
ΔM=ΔMs+ΔMe
(8)
式中,ΔMs為X射線源方向矢量觀測誤差引起的點(diǎn)積偏差,ΔMe為地球敏感器方向矢量觀測誤差引起的點(diǎn)積偏差。由ΔM引起的觀測量偏差為ΔY,令觀測量Y′為
Y′=Y+ΔY=arccos(M+ΔM)
(9)
真實(shí)量為
Y=arccos(M)
(10)
可得
ΔM=cos(Y+ΔY)-cos(Y)
(11)
式(11)進(jìn)行麥克勞林展開,有
(12)
略去高階小項(xiàng),整理得
(13)
根據(jù)式(7),對s求偏導(dǎo)數(shù),可得到X射線源觀測誤差引起的觀測量偏差ΔMs為
ΔMs=-r·Δs
(14)
即由X射線源方向矢量角位置偏移Δs引起的誤差傳播函數(shù)。當(dāng)Δs為角秒量級時(shí),ΔMs極小。上式對時(shí)間求導(dǎo),并取絕對值,得
(15)
近地圓軌道衛(wèi)星速率不超過8 km/s,軌道高度一般為200~800 km,X射線源天體角位置誤差一般為角秒量級。由此估算,|ΔMs|一小時(shí)變化率約為10-5量級??梢?X射線源角位置誤差引起的ΔMs極小且變化極為緩慢。
同理,式(7)兩端對r求偏導(dǎo)數(shù),得
ΔMe=-Δr·s
(16)
即由地心矢量角位置偏移Δr引起的誤差傳播函數(shù)。上式對時(shí)間求導(dǎo),并在等式兩邊取絕對值,可得地心方向矢量觀測誤差引起的單位方向矢量點(diǎn)積偏差ΔMe為
(17)
單圓錐掃描地球敏感器在軌測量誤差約為0.25°~1°,雙圓錐掃描地球敏感器在軌測量誤差約為0.05°~0.07°,面陣地球敏感器在軌測量誤差約為0.02°~0.05°,當(dāng)測量誤差為常值誤差或緩變量時(shí)[10],|ΔMe|的變化率極小。對于近地圓軌道衛(wèi)星,當(dāng)采用敏感器光軸在軌道平面內(nèi)與地平切線垂直的最佳UKF濾波安裝方式時(shí)[11,12],真實(shí)觀測量Y約為90°+arcsin(Re/rc),由于ΔM為緩變量,根據(jù)式(13),可將ΔY建立為緩變模型,即
(18)
式中,ws為高斯白噪聲。
衛(wèi)星軌道運(yùn)動(dòng)學(xué)采用牛頓二體引力模型[13],同時(shí)考慮地球非球形引力及其他攝動(dòng)力的作用,衛(wèi)星運(yùn)動(dòng)學(xué)模型建立為
(19)
式中,v為J2000地心慣性系下的衛(wèi)星速度矢量;μ為地球引力常數(shù);P(rc,v)為受攝運(yùn)動(dòng)項(xiàng)。將上式展開,并利用勒讓德多項(xiàng)式化簡,整理后即可得到衛(wèi)星運(yùn)動(dòng)狀態(tài)方程
(20)
式中,狀態(tài)變量取為x=[x,y,z,vx,vy,vz]T,分別為衛(wèi)星在J2000地心慣性系中3個(gè)坐標(biāo)軸上的位置和速度,J2為地球引力二階帶諧項(xiàng)系數(shù)。在這里,考慮到J2攝動(dòng)項(xiàng)相對其他攝動(dòng)項(xiàng)為大項(xiàng),因而在方程中單獨(dú)提出,ΔFx、ΔFy、ΔFz為地球非球形攝動(dòng)的高階攝動(dòng)項(xiàng)和日、月攝動(dòng)以及太陽光壓攝動(dòng)和大氣攝動(dòng)等其他攝動(dòng)力影響。wo為過程噪聲,假定為高斯白噪聲。
為了補(bǔ)償觀測量偏差引起的系統(tǒng)估計(jì)誤差,提高導(dǎo)航參數(shù)估計(jì)精度,根據(jù)系統(tǒng)測量殘差的分析,將X射線源矢量觀測誤差和地心矢量觀測誤差引起的系統(tǒng)殘差擴(kuò)展補(bǔ)充到系統(tǒng)狀態(tài)方程中,選取狀態(tài)量X=[xΔY]T,建立擴(kuò)展?fàn)顟B(tài)方程為
(21)
式中,w=[wows]T為擴(kuò)展?fàn)顟B(tài)噪聲,設(shè)為零均值高斯白噪聲。
式(4)、式(5)和式(21)構(gòu)成的自主導(dǎo)航系統(tǒng)是復(fù)雜的非線性系統(tǒng)。針對這一類系統(tǒng)的最優(yōu)濾波問題,Julier等提出了UKF濾波算法。該算法通過設(shè)計(jì)合理的采樣Sigma點(diǎn)逼近非線性系統(tǒng)狀態(tài)概率密度分布函數(shù),在提高濾波精度的同時(shí)規(guī)避了雅克比矩陣的計(jì)算[14]。經(jīng)過Unscented變換,UKF濾波可以獲得三階準(zhǔn)確度的非線性系統(tǒng)函數(shù)統(tǒng)計(jì)量,在非線性動(dòng)態(tài)濾波中得到越來越多的應(yīng)用[15-16]。因此,本文應(yīng)用UKF算法對自主導(dǎo)航系統(tǒng)狀態(tài)量進(jìn)行濾波計(jì)算,非線性系統(tǒng)的狀態(tài)方程式(21)和觀測方程式(6)經(jīng)離散化后,可整理為以下形式,即
Xk+1=F(Xk)+Wk
(22)
Yob,k=G(Xk)+Vk
(23)
式中,k代表離散化采樣點(diǎn)序列,狀態(tài)向量為Xk∈RL,測量向量為Yob,k∈RM,L和M分別為狀態(tài)向量維數(shù)和觀測向量維數(shù),Wk和Vk分別為狀態(tài)噪聲和測量噪聲,滿足
(24)
式中,Qk≥0為狀態(tài)噪聲方差陣;Rk>0為測量噪聲方差陣;δkj為Kronecker符號。
,i=0
(25)
(26)
(27)
Sigma點(diǎn)在非線性系統(tǒng)中傳播,并對其結(jié)果進(jìn)行加權(quán)計(jì)算,得到狀態(tài)變量和估計(jì)均方誤差陣的一步預(yù)測,即
χi,k/k-1=F(χi,k-1)
(28)
(29)
(30)
重新生成Sigma點(diǎn)并預(yù)測觀測向量為
(31)
(32)
(33)
計(jì)算增益矩陣為
(34)
(35)
(36)
對狀態(tài)變量和估計(jì)均方誤差陣進(jìn)行更新,可得到k時(shí)刻導(dǎo)航信息為
(37)
(38)
(39)
式中,β是狀態(tài)分布參數(shù),β=2 對于高斯分布是最優(yōu)的。
由于觀測量獲取時(shí)間間隔與導(dǎo)航計(jì)算周期往往是不同步的,因而需要基于導(dǎo)航周期最近估計(jì)點(diǎn)進(jìn)行導(dǎo)航時(shí)刻的位速估計(jì),從而獲得滿足導(dǎo)航與控制需求的實(shí)時(shí)位速信息。因此在導(dǎo)航周期中,進(jìn)行兩步獨(dú)立的UKF計(jì)算,一方面根據(jù)觀測量到達(dá)周期進(jìn)行UKF導(dǎo)航估計(jì),從而為下一周期的導(dǎo)航濾波提供狀態(tài)估值;另一方面根據(jù)當(dāng)前時(shí)刻t進(jìn)行UKF導(dǎo)航估計(jì),用于飛行器制導(dǎo)與控制決策。導(dǎo)航周期計(jì)算流程如圖3所示。其中tm 圖3 導(dǎo)航周期計(jì)算流程Fig.3 Calculation flow of navigation periods 為了驗(yàn)證本文提出方法在航天器自主導(dǎo)航應(yīng)用中的有效性,以實(shí)現(xiàn)近地圓軌道衛(wèi)星自主導(dǎo)航為目標(biāo)進(jìn)行仿真。 衛(wèi)星的初始軌道六根數(shù)為 半長軸a=7 208.72 km,偏心率e=0.011,軌道傾角i=65°,升交點(diǎn)赤經(jīng)Ω=30°,近地點(diǎn)角距ω=32°。 J2000地心慣性坐標(biāo)系下初始三軸位置估計(jì)誤差為1 km,初始三軸速度估計(jì)誤差為10 m/s,初始系統(tǒng)殘差估計(jì)誤差為0.5°,濾波器初始位置、速度估計(jì)方差陣為 仿真中所采用的系統(tǒng)過程噪聲矩陣為 Qk=diag[1 m 1 m 1 m 2×10-3m/s 2×10-3m/s 2×10-3m/s]2 測量殘差初始估計(jì)方差為(7°)2,測量殘差狀態(tài)噪聲為(0.001°)2,觀測量噪聲取為(4×10-6)°。 仿真中導(dǎo)航周期取為5 s。 參照實(shí)際敏感器能力,取X射線源矢量方向觀測精度3″,地心矢量方向觀測精度為0.02°,采用標(biāo)準(zhǔn)UKF方法進(jìn)行導(dǎo)航計(jì)算。衛(wèi)星在J2000地心慣性坐標(biāo)系下的位置、速度估計(jì)誤差隨時(shí)間變化如圖4和圖5所示。 圖4 位置估計(jì)誤差Fig.4 Position estimation error 圖5 速度估計(jì)誤差Fig.5 Velocity estimation error 在經(jīng)過一定時(shí)間的收斂過程后,濾波結(jié)果最終達(dá)到穩(wěn)定,圖6為測量殘差的估計(jì)誤差。 圖6 測量殘差估計(jì)誤差Fig.6 Observation residual estimation error 從局部放大圖可見,測量殘差估計(jì)誤差可達(dá)到10-4量級,這表明濾波器能夠?qū)崿F(xiàn)對測量殘差的有效估計(jì)。考慮到濾波收斂時(shí)間,需要在導(dǎo)航方案設(shè)計(jì)中,合理設(shè)置自主導(dǎo)航濾波起始點(diǎn),從而保證任務(wù)時(shí)間點(diǎn)的濾波精度。 按照以上條件進(jìn)行蒙特卡羅仿真,并對穩(wěn)定后的濾波結(jié)果求取均方根誤差均值,三軸位置和速度估計(jì)誤差如表1所示??梢?三軸位置估計(jì)誤差約為100~200 m,三軸速度估計(jì)誤差約為0.1~0.15 m/s。該精度水平能夠滿足航天器在軌自主導(dǎo)航的需求。 表1 均方根誤差均值 Table 1 Mean of RMSE 本文以近地航天器自主天文導(dǎo)航為應(yīng)用背景,提出一種基于X射線源/地心矢量觀測的航天器自主導(dǎo)航方法,將空間X射線探測器拓展應(yīng)用于航天器自主天文導(dǎo)航。進(jìn)行了系統(tǒng)建模與分析,將X射線源矢量觀測誤差和地心矢量觀測誤差引起的系統(tǒng)觀測殘差建模為緩變過程,并擴(kuò)展補(bǔ)充到狀態(tài)方程中。利用UKF算法進(jìn)行了導(dǎo)航信息解算。仿真表明,該方法能夠滿足航天器自主天文導(dǎo)航需求,具有一定的工程應(yīng)用參考價(jià)值。5 仿真分析
6 結(jié) 論