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

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

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

1 引言

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

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

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

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

2 軌道-姿態攝動模型

2.1 參考坐標系

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

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

2.2 動力學模型及假設

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

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

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

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

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

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

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

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

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

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

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

式中，R 為通用氣體常數。

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

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

圖3 類天宮一號目標飛行器軌降過程軸向力、法向力和俯仰力矩系數隨迎角變化曲線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 算法與狀態估計

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

6 軌道預報精度分析

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

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

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

7 結論

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

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

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

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