基于等效迎角的氣動融合軌道直接積分計算無控航天器軌道衰降研究

黨雷寧，李志輝，2*，唐小偉，梁 杰，彭傲平

（1.中國空氣動力研究與發展中心超高速空氣動力研究所，綿陽621000；2.國家計算流體力學實驗室，北京100191）

1 引言

低地球軌道（1000 km高度以下）運行的航天器的軌道高度因大氣阻力而衰減。尤其當航天器處于300 km高度以下的超低軌時，大氣密度急劇增加，空氣動力引起軌道的衰降更加顯著。考慮航天器超高速飛行過程空氣動力對軌道模擬計算的影響，不僅對航天器壽命末期可能的再入預警、隕落墜毀分析意義重大，而且是低軌航天器軌道設計和維持的重要研究手段。軌道衰降和再入點預報的準確程度，直接影響航天器再入時間、損毀和落區預報的準確程度。目前無控航天器隕落預報依靠測站軌道星歷數據支撐阻力系數的擬合及外推，精度較低，長期預報誤差可達數月，臨近隕落期預報精度約為1 h～10 min［1］。

大氣阻力對軌道高度的影響有2個方面［2-3］：①沿飛行軌道的熱層大氣總質量密度；②航天器的空氣動力特性，特別是阻力系數。自20世紀60年代以來，研究人員發展了Jacchia、MSIS系列、JB-2008等高層大氣模式［2，4］，極大地提高了熱層大氣密度預測的精準度。但對于軌道預報，現有的大氣模型在描述實際的大氣密度變化上還存在較大誤差，并依賴輻射流量和地磁指數的預報值［5］。

在軌航天器空氣動力系數是航天器外形、飛行高度、速度和飛行姿態的函數［3］。當這些影響因素確定后，可通過地面試驗和數值模擬等稀薄氣體動力學研究手段，獲得準確的空氣動力特性［6-8］。在其它因素都確定或者變化小的情況下，飛行姿態對空氣動力系數影響很大。當航天器正常工作、飛行姿態穩定，往往可以預測得到較準確的氣動力系數。如果航天器處于無控狀態，飛行器自旋，飛行姿態不確定性大，空氣動力系數的預測將十分困難。

大氣密度和空氣動力系數預測模型較大的誤差，使得失效航天器軌道衰降無法得到準確的預報，預報時間越長，誤差越大。Tapley［9］、Montenbruck等［10］較多地依賴于數值手段，通過引入各種待估參數來克服大氣密度和阻力系數的不確定度。蒼中亞等［11］采用TLE軌道根數和SGP4/SDP4模型，發展了基于歷史軌道數據的神經網絡最優阻力系數外推方法，并對國際空間站等大型航天器的軌道衰降進行預測。湯靖師等［5］提出把阻力系數、迎風面積、固定高度密度的組合參數作為外推變量，以開展短弧段和長弧段的軌道預報。從研究現狀看，現有的研究以歷史軌道數據為基礎，往往通過數學手段對阻力系數或者某種組合參數進行外推，從而開展預報，物理意義不明晰。

從空氣動力學的角度看，若航天器飛行姿態確定，則空氣動力系數可以準確計算得到，并隨高度呈現出可用空氣動力學解釋的變化規律。基于這點考慮，并結合傳統的基于歷史數據的軌道外推預報策略，本文提出無控航天器軌道衰降遵循等效迎角的概念和與之相適應的空氣動力快速計算建模方法，把稀薄氣體動力學融合進航天器飛行軌道衰降數值預報，在失效航天器軌道預報的眾多不確定因素中增加可計算建模確定性，提高軌道預報的有效性和精度。

2 軌道動力學計算模型與數值方法

在J2000.0坐標系Oxyz內航天器的運動方程為式（1）。

式中，r是航天器的位置矢量；F是航天器所受合外力引起質心運動的加速度矢量；t是飛行時間。

對于低軌航天器，上述軌道方程組中的外力主要是地球引力、空氣動力和日月引力。本文用地球引力位函數計算地球引力，其中略去田諧項，考慮至四階的帶諧項系數J2、J3、J4。地球的引力位函數可寫為式（2）。

式中，U是地球引力位，φ是地心緯度，r是飛行器質心到地球質心的距離，Re是參考橢球體赤道半徑。引力加速度G可通過對引力位求梯度得到，如式（3）所示。

在日、月引力模型中，將太陽、月球視作質點，根據非慣性坐標系中的受力分析，給出日、月引力模型［2］，其中太陽、月球位置由NASA的JPL DE405星歷計算得到。引力加速度計算公式可寫為式（4）。

式中，GS，L為飛行器受到的日（月）的引力加速度矢量；ΔS，L為慣性系中由太陽（月球）位置指向飛行器位置的矢量，ΔS，L是其大小；μ為太陽（月球）引力常數，μS=1.327 124 38×1020m3/s2，μL=1.327 124 38×1012m3/s2。下標S代表太陽，L代表月球。

空氣動力項模擬的準確性，直接決定軌道動力學方程求解的可靠性，其中大氣密度采用MSIS00模型計算［4］，空氣動力系數需要根據航天器當前飛行高度、速度與空間飛行姿態建模確定，詳見第3節。

文獻［2］討論了航天器運動方程的選擇及求解方法，認為對于方程式（1）中以r和·r為自變量的運動方程，右函數形式簡單但變化快，時間步長較小。基于此，并考慮到減小數值誤差的需求，本文以較小的時間步長（0.01 s），利用十階Adams預測-校正格式對軌道方程組進行數值積分。

3 無控航天器軌道衰降模擬的等效迎角與空氣動力系數一體化計算建模方法

3.1 基于空氣動力精細數值模擬驅動的低軌航天器氣動力特性快速算法

為了獲得航天器軌道衰降與再入全過程準確的氣動力特性并滿足沿軌道生產海量數據的工程需要，李志輝等［12-14］建立了由跨流域空氣動力學精細模擬方法驅動的當地化快速計算工具，并計算分析了大型航天器的氣動力特性［12］。

高超聲速空氣動力特性當地化快速算法是一種基于自由分子流與連續流的當地化橋函數理論關聯方法，基本思想是［7，12-13，15］：當地面元上的氣動力系數只依賴于來流和當地性質，如當地迎角、當地表面作用等；在連續流和自由分子流區域，氣動力特性遵從不同的物理定律，當地面元上的氣動力系數分別用連續流和自由分子流的方法計算；對于連續流和自由分子流之間的稀薄過渡流區，當地的氣動力系數用橋函數進行關聯適應。當地氣動力系數又可細分為當地壓力系數和摩擦力系數。

針對航天器軌道衰降過程空氣動力特性當地化快速計算，本文具體采用的方法［12］為：在連續流區域，采用修正牛頓流理論、二次激波膨脹方法以及背風真空效應法估算當地氣動力系數［16］；對自由分子流區域，使用基于不同模型材料修正的Nocilla壁面反射模型進行壓力與摩擦力系數計算［17］；在稀薄過渡流區域，采用修正的Boettcher/Legge非對稱橋函數理論［15］，建立了可分段描述的非對稱壓力與摩擦力系數關聯橋函數，以對連續流區和自由分子流區的氣動力系數進行關聯，其中包含3個壓力系數關聯參數Knn，p、ΔKnp1、ΔKnp2以及3個摩擦力系數關聯參數Knn，τ、ΔKnτ1、ΔKnτ2。基于這些關聯參數，稀薄過渡流區域的非對稱壓力橋函數Fb，p1和Fb，p2可寫為式（5）。

稀薄過渡流區域的摩擦力橋函數Fb，τ1和Fb，τ2可寫為式（6）。

式（5）～（6）中，Kn0，∞是飛行器克努森數，。

本文首先依靠DSMC方法和求解Boltzmann模型方程氣體動理論統一算法（GKUA）等跨流域氣動力精細模擬手段，對高度250～110 km的流場和氣動力系數進行精細模擬。然后基于氣動力系數精細模擬結果，調試選擇與計算修正這6個關聯參數，使當地化快速預測結果最大程度地符合精細數值模擬結果。

圖1給出利用前述的低軌航天器氣動特性一體化快速算法得到的無控隕落大型航天器B在高度350～100 km的氣動力特性。航天器B的外形包含了中間主體功能艙、兩側的太陽能電池帆板以及其他小部件（圖1（a）），氣動力參考面積取中間主體功能艙橫截面積。由圖1（b）中阻力和升力系數隨迎角變化曲線可見，阻力系數和升力系數隨迎角在0～180°范圍變化較大，阻力系數在0°和180°迎角下最小，在90°迎角下最大；升力系數隨迎角在0～180°范圍呈現反對稱正弦分布，且分別在8°和98°迎角改變符號。圖1（b）還比較了當地化快速算法（圖中標注為Local）得到的氣動力特性與DSMC數值模擬結果，可以看到兩者的變化規律吻合較好。由于DSMC方法是稀薄過渡流區域流場和氣動特性最可靠的計算方法，本文通過與DSMC結果比較，驗證了當地化快速算法的有效性和精度。由圖1（c）可見，阻力系數隨高度變化明顯，整體上隨高度的減小而減小，在高度180 km以上隨高度的減小而緩慢下降，而在高度180 km往下劇烈變化。

圖1 大型航天器B在高度350～100 km的氣動力特性Fig.1 Aerodynam ic characteristics of large spacecraft B at 350～100 km altitude

3.2 無控航天器軌道衰降模擬的等效迎角法

在航天器受到良好控制、姿態穩定的情況下，可采用3.1節所述的空氣動力學建模算法獲得氣動力特性并進行軌道衰降數值預報。然而很多失效航天器通常處于無控飛行狀態，不僅存在隨質心的三自由度軌道運動，還存在繞質心的自旋轉動。自旋航天器飛行姿態上的不確定性會嚴重影響空氣動力特性預測結果，進而影響軌道預報的精度。導致航天器自旋的因素眾多，現有的觀測手段很難準確甄別，且作為精確預測手段考慮轉動角速度的六自由度軌道與姿態動力學計算建模有待深化研究。觀測表明，作為本文研究對象的大型航天器B，因突發故障功能失效，處于無控飛行狀態后，其飛行姿態逐漸失穩，處于繞質心自旋狀態。研究分析同樣表明，歐空局ESA的對地觀測衛星ENVISAT在2012年4月失效后，在緩慢衰降的軌道上以1.5～3°/s的角速度自旋［18］。

從3.1節的計算分析可見，由于航天器B的升力隨姿態角的變化而改變符號，總體上正升力的迎角區間長度與負升力的迎角區間長度相等。因此當航天器自旋轉動一周后，升力對軌道的影響基本抵消。圖1（b）還可看出，升力與阻力相比是小量，升阻比小于0.05，從量值上說明升力的影響小。航天器所受大氣阻力，根據空氣動力學理論表達為式（7）。

其中，Cd是阻力系數，S是參考面積，V是航天器相對于大氣的飛行速度，ρ是大氣密度。軌道動力學領域和空氣動力學領域對參考面積S的定義和計算方法不同。軌道動力學通常認為［2］，S是垂直于速度V方向的迎風面積，隨飛行姿態變化；而空氣動力學通常把參考面積S取為飛行器不隨姿態變化的特征面積［3］，如底部面積、翼面面積、橫截面積等。引言部分已經分析了當前低軌航天器軌道衰降預報的國內外現狀。從空氣動力學角度分析，當前的研究沒有充分利用稀薄氣體動力學的研究成果，或者反演參數數量多，或者用于預報的數學外推缺乏物理意義。

本文基于空氣動力學專業角度，把式（7）中的參考面積S作為固定不變的無量綱參考量，提出另一種方案，即結合前述的空氣動力特性當地化快速算法與軌道計算模型，擬合確定航天器的姿態角，即等效迎角。等效迎角不是航天器實際飛行的姿態角，而是一種阻力姿態角，飛行器阻力系數的計算依賴于等效迎角、飛行速度和飛行高度。當采用等效迎角預報時，阻力系數隨飛行高度的變化可根據稀薄氣體動力學獲得的氣動力特性自動獲取，具有實際的物理意義。利用等效迎角預報無控航天器飛行軌道的步驟為：

1）對于某航天器外形，通過3.1節所述由空氣動力學精細模擬方法驅動的氣動力當地化快速算法，得到以速度、高度和迎角為自變量的氣動阻力系數數據表。

2）利用某一弧段的外測軌道星歷數據，擬合得到等效迎角。

3）把歷史弧段的等效迎角直接作為預報弧段的等效迎角；或通過連續幾個歷史弧段的等效迎角，優化擬合得到預報弧段的等效迎角。

4）基于預報弧段的等效迎角計算獲取航天器氣動力特性與飛行航跡，開展軌道預報。

等效迎角的擬合確定，作為初步研究，設計了如下迭代過程：

1）估算等效迎角 αd的范圍： α1≤αd≤α2，在此范圍內取n個離散點 αi（i=1，2，…，n）。

2）選擇外測軌道星歷數據中n個點的位置數據ro，j( j=1，2，…，n)作為比較數據。

3）選擇該歷史弧段的起點作為軌道計算初始點，對于每個迎角的離散值αi，都進行一次軌道計算，計算相同時刻模擬的軌道與n個要比較的歷史位置數據距離的平均值Ed，i，如式（8）所示。

式中，rc，j是與比較數據ro，j在相同時刻的模擬數據。為了加快計算速度，可在計算集群上采用并行計算技術，用n個計算核心分別計算n條軌道。

4）獲得（αi，Ed，i）序列后，則等效迎角 αd為Ed最小值對應的迎角，即式（9）。

4 計算結果及分析

以大型航天器A受控再入和大型航天器B無控飛行軌道衰降為研究對象，開展計算分析，航天器B的外形如圖1（a）所示，航天器A的外形與之相似。大型航天器A受控再入大氣層，在最后一次變軌制動后以固定姿態穩定飛行，迎角保持不變，星上GPS裝置測量到軌道衰降過程中的位置數據。大型航天器B失控后，姿態失穩并處于自旋狀態，通過地面雷達等設備得到軌道衰降過程中的位置數據。

4.1 大型航天器A受控隕落軌道衰降計算分析

大型航天器A受控隕落姿態穩定不變。本文依據軌道衰降過程中固定的飛行姿態，通過3.1節所述由空氣動力學精細模擬方法驅動的低軌航天器跨流域氣動力特性一體化快速算法得到阻力系數，使用3.2節所述的直接積分高精度軌道計算方法進行軌道計算。大型航天器A最后一次變軌后，以固定迎角飛行220～100 km過程的飛行航跡計算結果如圖2所示。由圖可見，在大型航天器A飛行高度從220 km下降到100 km的過程中，計算得到的空間位置（地固坐標系）結果與GPS實測數據吻合較好。詳細比較表明，計算的位置結果和GPS數據偏差不超過10 m。這說明本文發展的軌道計算以及低軌航天器氣動力特性計算模型及方法具有較好的可靠性與精度。

4.2 大型航天器B無控隕落軌道衰降計算分析

4.2.1 長弧段預報

利用3.2節介紹的無控航天器軌道衰降等效迎角方法，對大型航天器B無控飛行350～300 km軌道衰降過程進行分析。圖3（a）給出經過式（8）計算獲得的 （αi，Ed，i）序列曲線，即計算位置平均誤差隨迎角變化曲線，其中所選擇的歷史觀測數據弧長為1個月。根據式（9），圖3（a）曲線最低點所在的迎角就是此弧段擬合確定的等效迎角，量值為42°。預報弧段位于等效迎角擬合弧段之后1個月。把42°等效迎角作為預報弧段的等效迎角，開展預報弧段的軌道預報，得到預報弧段同一時刻預報位置相對測量位置的誤差δr，如圖3（b）所示。可看出，預報的空間位置誤差隨預報時間的增大而增大，預報5天的位置誤差約為42 km，預報10天的位置誤差約為316 km，預報30天的位置誤差約為1137 km。利用軌道速度估算，預報30天的位置誤差相當于在時間上偏差約150 s，占軌道周期約2.8%。這說明了本文發展的等效迎角氣動融合軌道直接積分高精度計算模型，可以得到較好的長弧段預報結果。

圖2 大型航天器A受控隕落220～100 km高度空間位置計算與GPS測量數據比較Fig.2 Comparison of orbit position between computation and GPS measurement for controlled reentry spacecraft A at 220～100 km altitude

圖3 大型航天器B無控飛行軌道衰降長弧段預報結果（預報弧段為1個月，等效迎角擬合弧段為預報弧段的前1月，軌道高度350～300 km）Fig.3 Com putational results of orbit decay for large spacecraft B（The prediction arc is onemonth and the inversion arc is onemouth before prediction arc.Orbit altitude is 350～300 km）

等效迎角反映了失穩自旋航天器在擬合弧段的平均阻力姿態角大小。當直接應用到預報弧段時，隱含的假設是預報弧段的航天器無控飛行姿態與擬合弧段的飛行姿態具有繼承連續性，彼此相差不大。圖3中軌道高度范圍（350～300 km）變化不大，外界擾動因素有可能沒有發生明顯變化，使得此時的無控航天器自旋運動模式和飛行姿態沒有發生明顯改變。

4.2.2 短弧段預報

當利用等效迎角氣動融合軌道直接積分建模方法開展短弧段預報時，擬合弧段需要保證較合適的長度，以反映等效迎角的變化。利用大型航天器B再入前5天的外測軌道星歷數據與3.2節介紹的方法，得到等效迎角的變化規律，如圖4所示，其中擬合弧長為1天。可看出，再入前5天等效迎角變化較大，從再入前倒數第5天的32.3°減小到再入當天的20.9°。

圖4 大型航天器B無控隕落再入大氣層前5日根據每日觀測數據得到的等效迎角Fig.4 Equivalent-angle-of-attacks of large spacecraft B obtained from dailymeasured data during 5 days before uncontrolled re-entry

圖5給出再入前倒數第3天軌道位置預報偏差，預報中采用的等效迎角通過再入前倒數第4天外測軌道星歷觀測數據擬合得到。可看出，預報誤差隨著預報時間的增大而增大，預報5 h軌道位置誤差為1.3 km，預報24 h軌道位置誤差為19.6 km。利用軌道速度估算，預報24 h的時間偏差約為2.6 s。這說明本文等效迎角氣動融合軌道直接積分模型方法可以得到準確性高的短弧段預報結果。

圖5 大型航天器B無控再入大氣層前倒數第3天軌道預報結果誤差（按倒數第4天外測軌道數據擬合等效迎角計算）Fig.5 Prediction error of orbit position on the 3rd day before uncotrolled re-entry of large spacecraft B（calculated by the fitted equivalent-angle-of-attack from measured orbit data on the 4th day before re-entry）

5 結論

1）對大型航天器A受控隕落220～100 km過程飛行航跡的計算分析，驗證了本文由空氣動力精細數值模擬驅動的低軌航天器氣動特性快速算法以及軌道動力學預報模型及方法的有效性。

2）在大型航天器B無控失穩自旋飛行軌道衰降計算分析中，得到了與外測軌道星歷數據吻合一致的長弧段和短弧段預報結果，在停止外測軌道數據供給5 h內，地心慣性系位置預報偏差低于1.5 km，驗證了基于等效迎角的氣動融合軌道直接積分計算模型對無控航天器軌道衰降預報的合理性。

3）所建立基于等效迎角的氣動融合軌道數值預報模型，是一種初步的方法。下一步可結合空氣動力學中的氣動參數辨識理論以及軌道動力學中的反演方法，對等效迎角進行更準確合理的擬合優化確定，為大型航天器再入解體分析提供更準確的再入點數據。