大型航天器無控飛行再入時間短期預報的軌道攝動方法研究

高興龍，李志輝，陳 欽，丁 娣，彭傲平

(1. 中國空氣動力研究與發(fā)展中心空氣動力學國家重點實驗室，綿陽621000；2. 中國空氣動力研究與發(fā)展中心設備設計及測試技術研究所，綿陽621000；3. 中國空氣動力研究與發(fā)展中心超高速空氣動力研究所，綿陽621000；4. 北京前沿創(chuàng)新中心國家計算流體力學實驗室，北京100191；5. 中國空氣動力研究與發(fā)展中心計算空氣動力研究所，綿陽621000)

1 引言

大型航天器壽命末期離軌再入大氣層時一旦功能失效，只能由受控再入轉為依靠無控飛行軌道衰降再入，無控再入的大型航天器通常由于自身系統(tǒng)如對地通信等故障造成失控，最終自然隕落再入大氣層解體，殘骸碎片墜落地球。 歷史上大型航天器失效后再入大氣層的案例較多，如美國的Skylab 空間站[1]。 大型航天器壽命終止后離軌自然隕落，將在低地球軌道環(huán)境氣動力/矩和地球引力[2]等攝動因素影響下，軌道高度逐漸衰降至最終再入大氣層[3]。 當缺少導航信息或信息服務數(shù)據(jù)中斷時，需要獨立依靠軌道動力學模型對航天器的軌道飛行航跡進行計算預報，盡可能地實現(xiàn)航天器空間位置速度狀態(tài)的可行性預測。 預報初值軌道參數(shù)精度和動力學模型精度直接影響航天器星歷預報精度。 數(shù)值積分法利用確定的初值逐步遞推航天器的位置速度或瞬時軌道根數(shù)，可以充分考慮各種攝動因素的影響，對于短期的高精度軌道預報較適用。

李志輝等[4-5]的研究表明，空氣動力對低軌航天器的軌道預報影響顯著，低地球軌道環(huán)境大型復雜結構航天器的氣動力隨高程顯著改變，如300～200 km 變化幅度達8%～50%。 即使考慮5%的阻力偏差，低軌空間目標的24 h 位置計算精度最大預測偏差可達64 km[6]。 而且李志輝等[4，7]的初步研究表明，對低地球軌道環(huán)境無控航天器的軌道衰降模擬，必須充分考慮目標飛行器軌道衰減過程所受到的跨流域空氣動力特性隨飛行姿態(tài)變化影響。 同時大部分非受控再入航天器都處于姿態(tài)失穩(wěn)飛行狀態(tài)，研究壽命末期大型航天器在低軌甚至超低軌道情況下的軌道動力學特性，需要同時考慮飛行器姿態(tài)的影響[8]。 Hart等[9]利用解析的自由分子流氣動模型開展了低軌空間目標的快速軌道外推和不確定參數(shù)敏感性分析，并對目標在軌壽命進行了預測。 Sommer等[10]對中國天宮一號再入過程的姿態(tài)運動和截面積進行了分析。 Nazarenko[11]利用雙行根數(shù)(Two-Line Element，TLE)軌道數(shù)據(jù)考慮隨機大氣環(huán)境擾動對天宮一號的再入時刻進行了預報研究。

軌道預報精度的誤差因素主要來自定軌初值誤差和動力學模型誤差，通常需要依靠濾波參數(shù)辨識方法盡可能消除誤差因素的影響，提高預報精度。 鄭作亞等[12]將無跡卡爾曼濾波(Unscented Kalman Filter，UKF)算法應用在GPS 衛(wèi)星軌道的短期預報中，有效提高了軌道預報精度和穩(wěn)定性。

傳統(tǒng)的基于常值阻力系數(shù)及數(shù)據(jù)擬合的氣動參數(shù)的預報方法在軌道預報問題上沒有充分考慮姿態(tài)變化引起的空氣動力的影響，因此需要結合姿態(tài)和跨流域氣動環(huán)境將更高精度的氣動特性建模方法融合進軌道動力學模型中。 為此，本文基于稀薄跨流域空氣動力數(shù)值模擬驅動設計的大型航天器氣動特性當?shù)鼗焖偎惴?，求解目標航天器不同高度和飛行姿態(tài)對應的氣動力、力矩系數(shù)，建立氣動融合姿態(tài)動力學軌道攝動模型，結合UKF 算法消除外測軌道星歷數(shù)據(jù)和數(shù)值積分誤差影響，對大型航天器進行短期的軌道外推計算預報。

2 軌道-姿態(tài)攝動模型

2.1 參考坐標系

無控航天器的軌道預報需要考慮姿態(tài)變化所引起的氣動力以及氣動力矩系數(shù)的變化，從而充分考慮姿態(tài)變化和氣動力因素對軌道攝動的影響。 如圖1 所示，建立軌道坐標系Obxryrzr，本體坐標系Obxbybzb， 速度坐標系Obxvyvzv和地心慣性坐標系Oxyz， 描述航天器位置速度和姿態(tài)的變化。

圖1 參考坐標系示意圖Fig.1 Sketch of reference frames

2.2 動力學模型及假設

目標航天器的軌道攝動方程由J2000 地心慣性坐標系下的位置和速度矢量描述，姿態(tài)描述為了避免奇異問題，采用四元數(shù)的方式，應用姿態(tài)四元數(shù)與方向余弦的關系，可得到姿態(tài)四元數(shù)的運動方程[13]，則目標航天器六自由度動力學方程具體形式為式(1)。

式 中， r ＝ x，y，z[ ]T為 位 置 矢 量， v ＝[ vx，vy，vz]T為速度矢量，q ＝ [q1，q2，q3，q4]T為軌道系到星體系的四元數(shù)，ω 為航天器的姿態(tài)角速度；μ 為地球引力常數(shù)， I 為航天器的慣量矩陣。攝動加速度項分別是: aN為N 體引力攝動加速度，本文考慮日、月的引力攝動加速度； aNSE為地球非球形攝動加度，是低軌航天器的主要攝動因素，本文計算采用精度較高的JGM-3 地球引力位模型；aNSL為月球非球形攝動加速度； aR為太陽光壓攝動加速度；aD為J2000 地心慣性坐標系下的大氣阻力攝動加速度，可由體坐標系下的加速度aB經(jīng)過坐標轉換得到，如式(2)所示。

式中，Bq為四元數(shù)的方向余弦矩陣， BOJ為軌道坐標系到地心慣性坐標系的轉換矩陣，aB可以根據(jù)氣動力系數(shù)寫成體坐標系下的三分量形式，如式(3)所示。

式中:CA、CN和CZ分別為目標航天器的軸向力系數(shù)、法向力系數(shù)和側向力系數(shù)，Sref為目標航天器阻力無量綱系數(shù)參考面積，ρ 為大氣密度，本文采用NRLMSISE00 大氣密度模型計算，Vr為航天器相對于大氣的運動速度，可以寫成式(4)。

式中，ωE為地球自轉角速度。

式(1)中右端的擾動力矩分別是:LG為重力梯度力矩，LR為太陽光壓力矩，LM為地磁力矩；LD為氣動擾動力矩，若記滾轉力矩系數(shù)為Cl，俯仰力矩系數(shù)為Cm，偏航力矩系數(shù)為Cn，參考面積為Sref，參考長度為Lref， 同樣可根據(jù)氣動力矩系數(shù)將LD寫成體坐標系下的三分量形式，如式(5)所示。

3 航天器跨流域氣動特性一體化快速計算

大型航天器的氣動力系數(shù)隨軌道高度和姿態(tài)的變化顯著，在軌道預報過程中，需要根據(jù)航天器再入過程的快速響應，考慮氣動特性計算的高效性。 如天宮一號目標飛行器，低軌道飛行過程是一個自數(shù)百到數(shù)十量級Knudsen 數(shù)高稀薄自由分子流動狀態(tài)多物理場復雜構型極高超聲速流動問題[14-16]。 為了捕捉空氣動力特性對軌道衰降過程影響，本文使用文獻[5]與文獻[17]發(fā)展的當?shù)鼗焖偎惴ǎ?分別在高度區(qū)域[250 km，200 km]、[200 km，120 km]、[120 km，100 km]每間隔5～3 km、2.5 km、2 km 各典型繞流狀態(tài)，使用跨流域空氣動力學數(shù)值方法計算，結果作為各子區(qū)域邊界值[4-5，18]。 而對于連接各子區(qū)域邊界的中間過渡區(qū)，采用修正的Boettcher/Legge 非對稱橋函數(shù)理論[17]，發(fā)展可分段描述的非對稱壓力與摩擦力系數(shù)關聯(lián)橋函數(shù)，在給定的物面角θb，當?shù)孛嬖系膲毫湍Σ亮ο禂?shù)分別為式(6)、(7)。

式中，F(xiàn)b，p和Fb，τ分別是壓力和摩擦力橋函數(shù)，依賴于獨立參數(shù)Kn、Tw/T∞和θb。

對于類天宮一號目標航天器可用三角形非結構網(wǎng)格表征其飛行器物形，采用通用近似的面元法對航天器進行物形處理，當面元劃分足夠細小時，引起的誤差就更小[19]，由此形成低軌道跨流域氣動力特性一體化快速計算技術。 在高稀薄自由分子流區(qū)，可采用基于不同模型材料修正的Nocilla 壁面反射模型進行壓力系數(shù)和摩擦力系數(shù)計算，在Maxwell 平衡態(tài)氣體分子速度分布的假設下，每個面元的壓力和摩擦力系數(shù)為式(8)、(9)。

式中，Tw、T∞分別為物面和來流溫度，σN和σT分別為法向和切向物面動量適應系數(shù)， erf()為誤差函數(shù)，s 為速度比，定義為式(10)。

式中，R 為通用氣體常數(shù)。

由此，充分利用高性能計算機開展DSMC 和求解Boltzmann 模型方程統(tǒng)一算法對相關關聯(lián)參數(shù)進行調(diào)試確定和計算修正，最終獲取適用類天宮一號目標飛行器覆蓋跨流域不同飛行高度和典型姿態(tài)的氣動力和力矩。 圖2 所示為計算得到的在350～100 km 高度下，航天器氣動阻力系數(shù)隨高度變化曲線。 從圖中的阻力系數(shù)變化趨勢可以看出，350～180 km 軌道高度范圍的阻力系數(shù)緩慢減小，之后隨著軌道高度的降低，阻力系數(shù)劇烈下降，這種異?，F(xiàn)象主要是由于跨過該高度，氣流密度、動壓迅速呈量級增加，致大氣阻力系數(shù)迅速減小，雖然大氣阻力有量綱量隨飛行高度降低呈現(xiàn)非線性增大。 圖3 為目標飛行器軸向力、法向力和俯仰力矩系數(shù)隨迎角變化曲線。

圖2 類天宮一號目標飛行器軌降過程氣動阻力變化Fig.2 Changes of aerodynamic drag characteristics of Tiangong-1-type spacecraft with altitude during the orbital decline process

圖3 類天宮一號目標飛行器軌降過程軸向力、法向力和俯仰力矩系數(shù)隨迎角變化曲線Fig.3 Changes of axial， normal force and pitch moment coefficients of Tiangong-1-type spacecraft with angle of attack during the orbital decline process

4 UKF 算法與狀態(tài)估計

利用初始外測目標飛行器軌道信息，直接積分軌道姿態(tài)動力學模型，同時根據(jù)預測狀態(tài)信息可以確定對應高度和姿態(tài)的氣動系數(shù)。 依次順序迭代計算。 同時結合UKF 思想，開展數(shù)據(jù)濾波，得到對航天器位置速度和姿態(tài)的跟蹤結果。 首先將位置、速度矢量的軌道參數(shù)和四元數(shù)、角速度矢量的姿態(tài)參數(shù)定義狀態(tài)量x 為式(11)。

同時將該問題以狀態(tài)方程、輸出方程和測量方程的形式表示如(12)所示[20]。

式中， x 為狀態(tài)向量，對應目標航天器的位置、速度、姿態(tài)角和角速度，y 為觀測向量，u 為系統(tǒng)的輸入，對于無控的目標航天器u ＝O；z 為外測得到的離散時間序列，包含測量誤差。 F 和G為過程和測量噪聲附加矩陣。 通過測量結果z 來對x 進行跟蹤和估計，采用UKF 算法對上述方程進行處理。 在UKF 算法中，對于非線性系統(tǒng)的狀態(tài)估計問題，其狀態(tài)隨機變量需要定義成原狀態(tài)量和過程、測量噪聲隨機變量的組合，即式(13)。

這樣系統(tǒng)擴展成L 維空間，相應的狀態(tài)協(xié)方差矩陣可寫成[21]式(14)。

式中， Px為狀態(tài)量的協(xié)方差矩陣， Rw和Rv為過程和測量噪聲協(xié)方差矩陣。 與EKF 算法相比，UKF 具有二階精度，且不需要計算Jacobian 和Hessian 矩陣，使得算法更易實現(xiàn)，并可用于不可微函數(shù)估計。

無跡變換(Unscented Transformation，UT)是無跡卡爾曼濾波的核心[21]，它通過構造有限樣本點(sigma points)來近似隨機變量在非線性變換下的統(tǒng)計特性。 比例無跡變換(Scaled Unscented Transformation，SUT)是在UT 的基礎上加入比例系數(shù)，用于在協(xié)方差矩陣保持正定的狀態(tài)下對樣本點進行調(diào)節(jié)。 基于SUT 變換的標準UKF 算法可分為預測和矯正兩步進行，算法流程如下:

1)初始化。 由先驗信息給出狀態(tài)量估計和協(xié)方差，假設噪聲信號為零均值白噪聲。

2)構造采樣點。 由SUT 變換構造2L＋1 個采樣點，即式(15)。

引入系數(shù)γ，其表達式為式(16)。

3)預測。 由k 狀態(tài)預估第k＋1 步的狀態(tài)量及均方差如式(17)所示。

4)校正。 由觀測量對預測量進行校正，首先計算增益矩陣K，由狀態(tài)量估計值可給出輸出量估計及協(xié)方差矩陣，然后由協(xié)方差矩陣計算增益矩陣，再由增益矩陣對系統(tǒng)狀態(tài)量和協(xié)方差矩陣進行校正。 若噪聲為累積白噪聲，則可去掉過程和測量噪聲。

與EKF 算法相比，UKF 具有二階精度，且不需要計算Jacobian 和Hessian 矩陣，使得算法更易實現(xiàn)，并可用于不可微函數(shù)估計。

該方法可以克服直接積分六自由度動力學模型長時間所帶來的誤差積累、長時間預報精度降低的問題，同時可以顯著提高計算效率，便于對軌道和姿態(tài)動力學模型分別進行誤差分析。

5 無控航天器軌道攝動模型仿真結果

利用前文所建立的數(shù)學模型和濾波方法對類天宮一號目標飛行器的無控飛行軌道變化進行計算驗證，圖4～圖7 為預測的該目標飛行器軌道根數(shù)和姿態(tài)變化。 計算時間取4.5 天，觀察其變化趨勢。 圖4 為類天宮一號目標飛行器軌道高度衰降變化曲線的仿真結果與TLE 觀測數(shù)據(jù)的對比結果，從初步的計算結果可以看出，軌道預測的仿真結果與觀測結果相符。 圖5 為軌道根數(shù)變化曲線，可以進一步看出目標飛行器臨近再入前幾天的軌道變化趨勢。

圖4 類天宮一號目標飛行器軌道高度變化Fig.4 Altitude variation of Tiangong-1-type space-

圖5 類天宮一號目標飛行器軌道根數(shù)變化Fig.5 Orbital elements variation of Tiangong-1-type spacecraft

從圖6 和圖7 的結果可以看出，目標飛行器姿態(tài)角和角速度的參數(shù)變化具有一定的隨機性，當前的測量手段還無法對無控大型航天器的姿態(tài)變化參數(shù)給出較精確的估計結果。 本文通過算例初步驗證了考慮姿態(tài)變化關聯(lián)氣動參數(shù)方法的可行性，未來還需要結合隨機理論和統(tǒng)計學方法進一步分析大型無控航天器姿態(tài)變化對軌道預報精度的影響，完善氣動融合軌道數(shù)值預報方法的技術手段。 總之，本文計算得到的軌道和姿態(tài)變化的預測結果均與測站偶爾觀測理論預期吻合相容，尤其是姿態(tài)動力學積分計算收斂、結果合理。

圖6 類天宮一號目標飛行器姿態(tài)角速度變化Fig.6 Attitude angle velocity variation of Tiangong-1-type spacecraft

圖7 類天宮一號目標飛行器姿態(tài)角變化Fig.7 Attitude angle variation of Tiangong-1-type spacecraft

6 軌道預報精度分析

將經(jīng)氣動融合參數(shù)辨識攝動模型計算得到的軌道預報數(shù)據(jù)與目標航天器實際觀測數(shù)據(jù)所獲取的軌道星歷信息進行對比，分析不同時刻的計算誤差，得到地心慣性坐標系下的位置矢量偏差和平均軌道高度偏差隨時間推進的變化趨勢，如表1 所示。

表1 無控航天器軌道預報精度Table 1 Orbit prediction accuracy of uncontrolled spacecraft

氣動力作為非保守力是引起目標航天器軌道形狀和大小改變的主要攝動因素之一。 從誤差分析結果可看出，24 h 內(nèi)的位置矢量偏差最大約為30 m，平均軌道高度偏差最大約為16 m，實際目標航天器的軌道高度在平均軌道高度上下做周期性運動；從計算精度變化趨勢來看，誤差都是隨著時間推進先增大后逐漸收斂，計算偏差的收斂有助于提高模型預報的精度，證實本文提出研究建立經(jīng)卡爾曼濾波誤差分析與參數(shù)辨識的氣動融合軌道攝動模型方案，對無控航天器軌道衰降模擬正確性有較強模擬能力。

7 結論

1) 本文基于J2000 地心慣性坐標系軌道攝動方程，發(fā)展UKF 和跨流域氣動特性一體化快速計算方法，建立氣動融合軌道攝動預報模型。

2) 針對大型航天器無控飛行軌道變化進行仿真與誤差精度對比分析，驗證了基于UKF 跨流域氣動參數(shù)解算軌道攝動模型可靠性。

3) 該方法可以為壽命末期的大型無控航天器結束軌道衰降轉為再入的時刻和彈道信息預報提供參考，并為壽命末期的目標航天器無控再入解體過程(桁架結構熱力響應變形軟化、復合材料熱解燒蝕、解體)數(shù)值預報的計算分析提供具有一定精度的軌道衰降再入輸入條件。

4) 進一步發(fā)展空氣動力融合軌道衰降的濾波誤差分析與參數(shù)辨識的軌道攝動模型，用于大型航天器無控飛行中長期軌道預報，充分發(fā)揮其快速高效是有價值意義之路。