龔秋武 毛永軍 楊舒農(nóng) 陳鼎 白云

(1 中國酒泉衛(wèi)星發(fā)射中心, 甘肅酒泉 732750)(2 中國人民解放軍95983部隊, 甘肅酒泉 732750)

據(jù)統(tǒng)計，截止2019年10月4日，太空中約有19 779個直徑超過10 cm的空間物體，其中包括5181個航天器和14 598個空間碎片[1]，這些都對正常運行的航天器的安全構(gòu)成了極大的威脅。因此，規(guī)避空間碎片已成為一項非常重要的應急航天任務。而該任務的前提是確定航天器的軌道根數(shù)，并預測其軌道，以確定它是否會與軌道附近的空間碎片相撞。此外，當某些自然災害發(fā)生時，例如火山爆發(fā)、地震和森林火災等，無法在地面上觀察到災害現(xiàn)場。但是，偵察衛(wèi)星和遙感衛(wèi)星等航天器卻可用于這種應急空間偵察任務，但前提也是要先確定其軌道根數(shù)，以便計算它們何時飛行到災害發(fā)生地的上空，從而開展空間偵察和探測。

對于航天器的軌道根數(shù)確定，美國空軍太空司令部(Air Force Space Command，AFSPC)發(fā)布的雙行根數(shù)(Two Line Elements，TLEs)可能是最準確的描述[2-3]。TLEs不是瞬時軌道根數(shù)，而是通過特殊方式消除周期性變化得到的“平均軌道根數(shù)”；TLEs具有標準化的格式，并且只能與SGP軌道模型配套使用。此外，應急航天任務還強調(diào)實時與自主性，以滿足時效性的要求。傳統(tǒng)的軌道確定方法將星上導航數(shù)據(jù)傳輸?shù)降孛嬲具M行計算，然后再將結(jié)果傳輸?shù)胶教炱鳌１M管其結(jié)果具有較高的精確度，但占用大量計算資源，且消耗大量時間進行數(shù)據(jù)傳輸。因此，當考慮到實時性和星載計算能力時，傳統(tǒng)的軌道確定方法對于應急航天任務是不可接受的[4-5]。

實際上，絕大多數(shù)軌道模型，例如二體模型、J2擾動和高精度軌道預測(High Precision Orbit Prediction，HPOP)模型等，都使用瞬時軌道根數(shù)進行軌道預報[6]。但是，TLEs是“平均”軌道根數(shù)，與瞬時軌道根數(shù)不同，盡管這兩種軌道根數(shù)非常相似，但是如果將TLEs替換為瞬時軌道根數(shù)代入SGP模型進行軌道預報，將會造成較大的預報誤差。因此，將瞬時軌道根數(shù)轉(zhuǎn)換為TLEs是非常必要的。但是，AFSPC并沒有公布這兩種軌道根數(shù)相互轉(zhuǎn)換的算法。為了解決這個問題，必須建立對應的算法將瞬時軌道根數(shù)轉(zhuǎn)換為可以用于SGP的雙行根數(shù)[7]。

文獻[8]建立了SGP4模型，并用于近地衛(wèi)星軌道預報，將文獻[9]的解析理論進行簡化，將重力模型和空氣動力學模型作為Brouwer和功率密度函數(shù)的解，便可以得到SGP4模型。該模型考慮了幾個長期和長周期擾動項，例如，J2，J3和J4項和空氣阻力等，以確保軌道預報的精度。但是由于未考慮短周期擾動項，SGP4模型的位置預報精度只能達到2 km[10]。文獻[11]根據(jù)TLEs和SGP4模型開發(fā)了一個軌道預報軟件包。文獻[12-13]提出了瞬時軌道根數(shù)和平均軌道根數(shù)之間的一種轉(zhuǎn)換方法。為了提高太空飛行的安全性，文獻[14-15]建立了基于SGP4和TLEs的空間碎片軌道預測方法。文獻[16]提出了兩種方法來修改TLEs的阻力系數(shù)，并將SGP4的預測精度提高了20%以上。文獻[17]對LEO碎片的阻力系數(shù)進行研究，并估計了2000多個空間碎片的阻力系數(shù)來檢驗TLEs的精度。文獻[18]運用批量最小二乘技術(shù)和測距率來校正歷元時刻的TLEs角度誤差，提高了圓形軌道確定的準確性。由于上述研究僅研究了TLEs部分元素的解算和轉(zhuǎn)換方法，因此需要對其他元素的解算和轉(zhuǎn)換開展更多研究。

本文即以航天器雙行根數(shù)為狀態(tài)變量，以飛行位置和速度為測量變量，提出基于無跡卡爾曼濾波的航天器雙行根數(shù)轉(zhuǎn)換方法，有效轉(zhuǎn)換得到TLEs的基本元素[19]，且達到TLEs格式的精度要求。

1 建模與坐標轉(zhuǎn)換

(1)

式中：s(t)、r(t)和v(t)都是矢量，用列向量表示；S(·)代表SGP4模型；α0代表TLEs的元素；t代表從歷元時刻t0開始的時間間隔。

(2)

已知SGP4模型和TLEs的參考坐標系為真春分點和平二分點(TEMED)參考系統(tǒng)。但是，導航數(shù)據(jù)的參考坐標系是世界大地測量系統(tǒng)(World Geodetic System 1984，WGS84)。因此，需要將導航數(shù)據(jù)經(jīng)過兩次坐標轉(zhuǎn)換，從WGS84轉(zhuǎn)換為J2000.0，再從J2000.0轉(zhuǎn)換為TEMED。第一次轉(zhuǎn)換可參考文獻[20]，第二次轉(zhuǎn)換可參考文獻[12]和[21]。

2 瞬時軌道根數(shù)與雙行根數(shù)轉(zhuǎn)換方法

文獻[12]和[13]中，均提出了瞬時軌道根數(shù)轉(zhuǎn)換為雙行根數(shù)的迭代算法，算法根據(jù)航天器任意時刻的飛行狀態(tài)矢量來實現(xiàn)瞬根與雙行根數(shù)的轉(zhuǎn)換。但是，迭代算法不能充分利用航天器的全部飛行狀態(tài)數(shù)據(jù)；其次，由于在整個迭代過程中未更新用來對比的標準飛行狀態(tài)矢量，因此無法消除測量噪聲的影響；另外，迭代算法不能解算出TLEs中的大氣阻力系數(shù)。

為了彌補上述算法的不足，充分利用航天器飛行狀態(tài)數(shù)據(jù)并消除測量噪聲的影響，本節(jié)建立卡爾曼濾波算法將瞬根轉(zhuǎn)換為雙行根數(shù)?？柭鼮V波是基于狀態(tài)估計理論而建立的，最常用的是在此卡爾曼濾波理論上發(fā)展而來的擴展卡爾曼濾波(Extended Kalman Filter，EKF)和無跡卡爾曼濾波(Unscented Kalman Filter，UKF)等。

本文基于UKF建立航天器雙行根數(shù)轉(zhuǎn)換方法，其中航天器雙行根數(shù)是狀態(tài)變量，而飛行狀態(tài)是測量變量。UKF最特殊優(yōu)勢的是利用無跡變換(Unscented Transform，UT)而不是線性化來處理均值和協(xié)方差的傳遞[22]。對于UT，為了獲取原始狀態(tài)矢量的均值和協(xié)方差，選擇了一組選定的sigma點。然后，通過非線性模型進行映射。接著，計算sigma點的均值和協(xié)方差。這樣，對于任意的非線性函數(shù)[23-26]，均值和協(xié)方差可以達到二階精度。

未選擇擴展卡爾曼濾波主要是考慮到EKF需要對非線性測量方程進行線性化，會產(chǎn)生線性化誤差，嚴重時使得濾波發(fā)散，其次由于SGP4模型的復雜性與非線性，無法使用偏微分計算測量矩陣，而是用數(shù)值差商則會大量增加計算量。

X(t)=I7X0

(3)

式中：X(t)表示t歷元時刻雙行根數(shù)；I7表示7維單位陣；X0表示初始歷元時刻雙行根數(shù)。

測量方程類似于式(1)

Z(t)=S(X0,t)+V

(4)

由于測量數(shù)據(jù)是離散的，因此可以將式(4)改為

Z(k)=rk=S(X0,k)r+V(k),(k=1,2,3…)

(5)

式中：S(X0,k)r表示離散狀態(tài)的預測位置矢量，V(k)表示離散狀態(tài)的測量噪聲矢量。

sigma點及其權(quán)重通過式(6)和(7)計算

(6)

(7)

詳細的過程可參考文獻[26]。顯然，UKF沒有線性化。它實現(xiàn)了UT以使sigma點的均值和協(xié)方差與原始特征相匹配，然后將應用sigma點的非線性映射得到狀態(tài)量的概率密度函數(shù)。

3 數(shù)值仿真

前兩節(jié)中，建立了算法和模型，下面將通過數(shù)值仿真來驗證有效性和準確性。航天器的參數(shù)是從STK軟件的TLEs數(shù)據(jù)庫中選擇的，而測量數(shù)據(jù)是由STK仿真模擬得到的。選擇tle-00005航天器來實施仿真(tle-00005航天器是STK軟件包的TLEs數(shù)據(jù)庫中的一個航天器，其名稱為Vanguard1是美國海軍在1958年3月17日發(fā)射的，周期約為133 min)，由SGP4模型的特點可知，當軌道周期小于255 min時，SGP4的預報性能更好[2]。根據(jù)UKF算法仿真需要，挑選該航天器TLEs中的部分參數(shù)，如表1所示。

表1 tle-00005航天器雙行根數(shù)Table 1 TLEs of tle-00005 spacecraft

濾波參數(shù)的設置決定了濾波器的性能，UKF的參數(shù)設置為

(8)

測量數(shù)據(jù)由STK仿真得到，且測量時間間隔設置為1 min。根據(jù)GPS和北斗等導航設備的定位精度可知，目前民用的定位精度約為10 m(σ)量級[27]，因而測量噪聲參數(shù)設置為

R=diag(10,10,10)

(9)

式中：diag(·)表示對角矩陣。

以一天為預報時長，分別利用STK軟件和本文的程序進行軌道預報，將STK預報值作為標準值，通過疊加小量到TLEs上，使本文的程序預報值與標準值之間的位置誤差約為1 m(σ/10)，基于此，可以估計過程噪聲矩陣。初始狀態(tài)協(xié)方差的設置要避免濾波器發(fā)散，通過參考過程噪聲來設置它。進而，給出過程噪聲和初始狀態(tài)協(xié)方差矩陣的值如下：

Q=diag(10-7,10-7,10-7,10-7,10-7,10-10,10-7)

(10)

P0=diag(10-7,10-7,10-7,10-7,10-7,10-7,10-7)

(11)

利用STK仿真，得到航天器在WGS84系中的飛行狀態(tài)。任意選擇一組狀態(tài)量，本文中選取2007年5月18日05：02：19.609時刻(歷元初始時刻)的位置和速度矢量，如表2所示(這些數(shù)據(jù)均不含測量噪聲)。

表2 tle-00005航天器的飛行狀態(tài)參數(shù) Table 2 Flight state of tle-00005 spacecraft

根據(jù)二體運動方程的軌道積分、活力積分等公式，可推導得到軌道根數(shù)與位置/速度的關(guān)系[28]，利用無測量噪聲的位置和速度數(shù)據(jù)，即可計算得到濾波狀態(tài)初值的前6個參數(shù)，但第7個參數(shù)大氣阻力系數(shù)無法通過該方法計算，因此該參數(shù)文中采用經(jīng)驗值，取為-1.000 0×10-4，則航天器雙行根數(shù)的初始狀態(tài)值為

(12)

然后代入UKF算法進行計算，并將得到的雙行根數(shù)與表1的參考值進行比較，結(jié)果見表3和圖1。其中圖1橫坐標的數(shù)字1～7分別代表軌道傾角、升交點赤經(jīng)、偏心率、近地點角距、平近點角、平均運動角速度和大氣阻力系數(shù)等7個TLEs參數(shù)(均轉(zhuǎn)換為誤差大小，不含符號)。選擇偏心率和平均運動角速度的絕對誤差曲線來展示UKF濾波過程，如圖2所示。

表3 UKF算法結(jié)果Table 3 Results of UKF

圖1 UKF精度(誤差大小，不含符號)Fig.1 Accuracy of UKF results (without sign)

圖2 UKF的濾波過程Fig.2 Procedure of UKF

從表3和圖2可知，UKF能有效將瞬根轉(zhuǎn)換為雙行根數(shù)。當測量噪聲位置誤差為10 m時，4個角度的精度至少為10-5(°)，偏心率精度為10-7，平均運動角速度精度為10-9圈/天，大氣阻力系數(shù)精度為10-6。

此外，使用濾波算法的結(jié)果進行軌道預報并與STK值進行比較，結(jié)果如圖3所示，位置誤差在10 m以內(nèi)，速度誤差在0.008 m/s以內(nèi)?？梢酝茢啵敎y量噪聲較小時，UKF精度很高。

圖3 利用UKF結(jié)果進行軌道預報的位置和速度誤差Fig.3 Position and velocity propagation deviations of UKF results

為了說明大噪聲對濾波算法的影響，下面疊加噪聲到導航數(shù)據(jù)中然后再進行仿真。類似地，噪聲的均值設置為零，而方差分別設置為3σ和10σ，其中σ=10。結(jié)果的絕對誤差如圖4所示，其中橫坐標的數(shù)字1～7所表示的參數(shù)與圖1相同。將結(jié)果代入SGP4模型進行軌道預報并與STK值相比較，得到的位置和速度預報誤差如圖5所示。

圖4 測量噪聲對濾波算法的影響(誤差大小，不含符號)Fig.4 Influence of measurement noise to filter algorithm (without sign)

圖5 利用UKF在大噪聲條件下的結(jié)果進行軌道預報的位置和速度誤差Fig.5 Position and velocity propagation deviations of UKF results under big noise condition

從圖4、5可知，TLEs的絕對誤差以及位置和速度預報誤差均隨著噪聲方差的增大而增大。但是在兩個仿真算例中，位置預報誤差均小于測量噪聲的方差。

對于其他特殊情況，如較大的初始誤差，同樣進行了數(shù)值仿真。通過之前仿真，可以知道初始狀態(tài)值是由二體模型在導航數(shù)據(jù)不包含測量噪聲條件下計算得到的。在此，對用來計算初始狀態(tài)的導航數(shù)據(jù)疊加噪聲，然后再利用二體模型計算初始狀態(tài)，從而進行較大初始誤差情況的數(shù)值仿真。

測量噪聲的均值設置為零，方差分別設置為10σ和100σ，其中σ=10[28]，其他參數(shù)設置與正常情況一致，仿真結(jié)果如表4所示。

表4 大初始誤差條件下的結(jié)果Table 4 Results under large initial error condition

從表4中可以看出，在較大初始誤差條件下，UKF算法仍可以得到滿足TLEs格式精度的航天器雙行根數(shù)。

4 結(jié)束語

本文研究了基于無跡卡爾曼濾波的航天器雙行軌道根數(shù)轉(zhuǎn)換方法。在較小測量噪聲條件下，UKF性能很高。當測量噪聲位置誤差為10 m時，UKF的結(jié)果，4個角度的精度為10-5(°)，偏心度精度為10-7，平均運動角速度精度為10-9圈/天，大氣阻力系數(shù)精度為10-6，均達到了TLEs標準格式的精度要求。利用轉(zhuǎn)換得到的TLEs進行軌道預報，在一天時間內(nèi)，位置預報誤差在10 m以內(nèi)，速度預報誤差均在0.008 m/s以內(nèi)。在某些特殊情況下，如較大初始誤差和較大測量噪聲，UKF仍能有效轉(zhuǎn)換得到滿足精度要求的航天器雙行根數(shù)。