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

（1.中國(guó)空氣動(dòng)力研究與發(fā)展中心超高速空氣動(dòng)力研究所，綿陽(yáng)621000；2.國(guó)家計(jì)算流體力學(xué)實(shí)驗(yàn)室，北京100191）

1 引言

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

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

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

大氣密度和空氣動(dòng)力系數(shù)預(yù)測(cè)模型較大的誤差，使得失效航天器軌道衰降無(wú)法得到準(zhǔn)確的預(yù)報(bào)，預(yù)報(bào)時(shí)間越長(zhǎng)，誤差越大。Tapley［9］、Montenbruck等［10］較多地依賴于數(shù)值手段，通過(guò)引入各種待估參數(shù)來(lái)克服大氣密度和阻力系數(shù)的不確定度。蒼中亞等［11］采用TLE軌道根數(shù)和SGP4/SDP4模型，發(fā)展了基于歷史軌道數(shù)據(jù)的神經(jīng)網(wǎng)絡(luò)最優(yōu)阻力系數(shù)外推方法，并對(duì)國(guó)際空間站等大型航天器的軌道衰降進(jìn)行預(yù)測(cè)。湯靖師等［5］提出把阻力系數(shù)、迎風(fēng)面積、固定高度密度的組合參數(shù)作為外推變量，以開(kāi)展短弧段和長(zhǎng)弧段的軌道預(yù)報(bào)。從研究現(xiàn)狀看，現(xiàn)有的研究以歷史軌道數(shù)據(jù)為基礎(chǔ)，往往通過(guò)數(shù)學(xué)手段對(duì)阻力系數(shù)或者某種組合參數(shù)進(jìn)行外推，從而開(kāi)展預(yù)報(bào)，物理意義不明晰。

從空氣動(dòng)力學(xué)的角度看，若航天器飛行姿態(tài)確定，則空氣動(dòng)力系數(shù)可以準(zhǔn)確計(jì)算得到，并隨高度呈現(xiàn)出可用空氣動(dòng)力學(xué)解釋的變化規(guī)律?；谶@點(diǎn)考慮，并結(jié)合傳統(tǒng)的基于歷史數(shù)據(jù)的軌道外推預(yù)報(bào)策略，本文提出無(wú)控航天器軌道衰降遵循等效迎角的概念和與之相適應(yīng)的空氣動(dòng)力快速計(jì)算建模方法，把稀薄氣體動(dòng)力學(xué)融合進(jìn)航天器飛行軌道衰降數(shù)值預(yù)報(bào)，在失效航天器軌道預(yù)報(bào)的眾多不確定因素中增加可計(jì)算建模確定性，提高軌道預(yù)報(bào)的有效性和精度。

2 軌道動(dòng)力學(xué)計(jì)算模型與數(shù)值方法

在J2000.0坐標(biāo)系Oxyz內(nèi)航天器的運(yùn)動(dòng)方程為式（1）。

式中，r是航天器的位置矢量；F是航天器所受合外力引起質(zhì)心運(yùn)動(dòng)的加速度矢量；t是飛行時(shí)間。

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

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

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

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

空氣動(dòng)力項(xiàng)模擬的準(zhǔn)確性，直接決定軌道動(dòng)力學(xué)方程求解的可靠性，其中大氣密度采用MSIS00模型計(jì)算［4］，空氣動(dòng)力系數(shù)需要根據(jù)航天器當(dāng)前飛行高度、速度與空間飛行姿態(tài)建模確定，詳見(jiàn)第3節(jié)。

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

3 無(wú)控航天器軌道衰降模擬的等效迎角與空氣動(dòng)力系數(shù)一體化計(jì)算建模方法

3.1 基于空氣動(dòng)力精細(xì)數(shù)值模擬驅(qū)動(dòng)的低軌航天器氣動(dòng)力特性快速算法

為了獲得航天器軌道衰降與再入全過(guò)程準(zhǔn)確的氣動(dòng)力特性并滿足沿軌道生產(chǎn)海量數(shù)據(jù)的工程需要，李志輝等［12-14］建立了由跨流域空氣動(dòng)力學(xué)精細(xì)模擬方法驅(qū)動(dòng)的當(dāng)?shù)鼗焖儆?jì)算工具，并計(jì)算分析了大型航天器的氣動(dòng)力特性［12］。

高超聲速空氣動(dòng)力特性當(dāng)?shù)鼗焖偎惴ㄊ且环N基于自由分子流與連續(xù)流的當(dāng)?shù)鼗瘶蚝瘮?shù)理論關(guān)聯(lián)方法，基本思想是［7，12-13，15］：當(dāng)?shù)孛嬖系臍鈩?dòng)力系數(shù)只依賴于來(lái)流和當(dāng)?shù)匦再|(zhì)，如當(dāng)?shù)赜恰?dāng)?shù)乇砻孀饔玫?；在連續(xù)流和自由分子流區(qū)域，氣動(dòng)力特性遵從不同的物理定律，當(dāng)?shù)孛嬖系臍鈩?dòng)力系數(shù)分別用連續(xù)流和自由分子流的方法計(jì)算；對(duì)于連續(xù)流和自由分子流之間的稀薄過(guò)渡流區(qū)，當(dāng)?shù)氐臍鈩?dòng)力系數(shù)用橋函數(shù)進(jìn)行關(guān)聯(lián)適應(yīng)。當(dāng)?shù)貧鈩?dòng)力系數(shù)又可細(xì)分為當(dāng)?shù)貕毫ο禂?shù)和摩擦力系數(shù)。

針對(duì)航天器軌道衰降過(guò)程空氣動(dòng)力特性當(dāng)?shù)鼗焖儆?jì)算，本文具體采用的方法［12］為：在連續(xù)流區(qū)域，采用修正牛頓流理論、二次激波膨脹方法以及背風(fēng)真空效應(yīng)法估算當(dāng)?shù)貧鈩?dòng)力系數(shù)［16］；對(duì)自由分子流區(qū)域，使用基于不同模型材料修正的Nocilla壁面反射模型進(jìn)行壓力與摩擦力系數(shù)計(jì)算［17］；在稀薄過(guò)渡流區(qū)域，采用修正的Boettcher/Legge非對(duì)稱橋函數(shù)理論［15］，建立了可分段描述的非對(duì)稱壓力與摩擦力系數(shù)關(guān)聯(lián)橋函數(shù)，以對(duì)連續(xù)流區(qū)和自由分子流區(qū)的氣動(dòng)力系數(shù)進(jìn)行關(guān)聯(lián)，其中包含3個(gè)壓力系數(shù)關(guān)聯(lián)參數(shù)Knn，p、ΔKnp1、ΔKnp2以及3個(gè)摩擦力系數(shù)關(guān)聯(lián)參數(shù)Knn，τ、ΔKnτ1、ΔKnτ2?；谶@些關(guān)聯(lián)參數(shù)，稀薄過(guò)渡流區(qū)域的非對(duì)稱壓力橋函數(shù)Fb，p1和Fb，p2可寫為式（5）。

稀薄過(guò)渡流區(qū)域的摩擦力橋函數(shù)Fb，τ1和Fb，τ2可寫為式（6）。

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

本文首先依靠DSMC方法和求解Boltzmann模型方程氣體動(dòng)理論統(tǒng)一算法（GKUA）等跨流域氣動(dòng)力精細(xì)模擬手段，對(duì)高度250～110 km的流場(chǎng)和氣動(dòng)力系數(shù)進(jìn)行精細(xì)模擬。然后基于氣動(dòng)力系數(shù)精細(xì)模擬結(jié)果，調(diào)試選擇與計(jì)算修正這6個(gè)關(guān)聯(lián)參數(shù)，使當(dāng)?shù)鼗焖兕A(yù)測(cè)結(jié)果最大程度地符合精細(xì)數(shù)值模擬結(jié)果。

圖1給出利用前述的低軌航天器氣動(dòng)特性一體化快速算法得到的無(wú)控隕落大型航天器B在高度350～100 km的氣動(dòng)力特性。航天器B的外形包含了中間主體功能艙、兩側(cè)的太陽(yáng)能電池帆板以及其他小部件（圖1（a）），氣動(dòng)力參考面積取中間主體功能艙橫截面積。由圖1（b）中阻力和升力系數(shù)隨迎角變化曲線可見(jiàn)，阻力系數(shù)和升力系數(shù)隨迎角在0～180°范圍變化較大，阻力系數(shù)在0°和180°迎角下最小，在90°迎角下最大；升力系數(shù)隨迎角在0～180°范圍呈現(xiàn)反對(duì)稱正弦分布，且分別在8°和98°迎角改變符號(hào)。圖1（b）還比較了當(dāng)?shù)鼗焖偎惴ǎ▓D中標(biāo)注為L(zhǎng)ocal）得到的氣動(dòng)力特性與DSMC數(shù)值模擬結(jié)果，可以看到兩者的變化規(guī)律吻合較好。由于DSMC方法是稀薄過(guò)渡流區(qū)域流場(chǎng)和氣動(dòng)特性最可靠的計(jì)算方法，本文通過(guò)與DSMC結(jié)果比較，驗(yàn)證了當(dāng)?shù)鼗焖偎惴ǖ挠行院途?。由圖1（c）可見(jiàn)，阻力系數(shù)隨高度變化明顯，整體上隨高度的減小而減小，在高度180 km以上隨高度的減小而緩慢下降，而在高度180 km往下劇烈變化。

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

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

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

從3.1節(jié)的計(jì)算分析可見(jiàn)，由于航天器B的升力隨姿態(tài)角的變化而改變符號(hào)，總體上正升力的迎角區(qū)間長(zhǎng)度與負(fù)升力的迎角區(qū)間長(zhǎng)度相等。因此當(dāng)航天器自旋轉(zhuǎn)動(dòng)一周后，升力對(duì)軌道的影響基本抵消。圖1（b）還可看出，升力與阻力相比是小量，升阻比小于0.05，從量值上說(shuō)明升力的影響小。航天器所受大氣阻力，根據(jù)空氣動(dòng)力學(xué)理論表達(dá)為式（7）。

其中，Cd是阻力系數(shù)，S是參考面積，V是航天器相對(duì)于大氣的飛行速度，ρ是大氣密度。軌道動(dòng)力學(xué)領(lǐng)域和空氣動(dòng)力學(xué)領(lǐng)域?qū)⒖济娣eS的定義和計(jì)算方法不同。軌道動(dòng)力學(xué)通常認(rèn)為［2］，S是垂直于速度V方向的迎風(fēng)面積，隨飛行姿態(tài)變化；而空氣動(dòng)力學(xué)通常把參考面積S取為飛行器不隨姿態(tài)變化的特征面積［3］，如底部面積、翼面面積、橫截面積等。引言部分已經(jīng)分析了當(dāng)前低軌航天器軌道衰降預(yù)報(bào)的國(guó)內(nèi)外現(xiàn)狀。從空氣動(dòng)力學(xué)角度分析，當(dāng)前的研究沒(méi)有充分利用稀薄氣體動(dòng)力學(xué)的研究成果，或者反演參數(shù)數(shù)量多，或者用于預(yù)報(bào)的數(shù)學(xué)外推缺乏物理意義。

本文基于空氣動(dòng)力學(xué)專業(yè)角度，把式（7）中的參考面積S作為固定不變的無(wú)量綱參考量，提出另一種方案，即結(jié)合前述的空氣動(dòng)力特性當(dāng)?shù)鼗焖偎惴ㄅc軌道計(jì)算模型，擬合確定航天器的姿態(tài)角，即等效迎角。等效迎角不是航天器實(shí)際飛行的姿態(tài)角，而是一種阻力姿態(tài)角，飛行器阻力系數(shù)的計(jì)算依賴于等效迎角、飛行速度和飛行高度。當(dāng)采用等效迎角預(yù)報(bào)時(shí)，阻力系數(shù)隨飛行高度的變化可根據(jù)稀薄氣體動(dòng)力學(xué)獲得的氣動(dòng)力特性自動(dòng)獲取，具有實(shí)際的物理意義。利用等效迎角預(yù)報(bào)無(wú)控航天器飛行軌道的步驟為：

1）對(duì)于某航天器外形，通過(guò)3.1節(jié)所述由空氣動(dòng)力學(xué)精細(xì)模擬方法驅(qū)動(dòng)的氣動(dòng)力當(dāng)?shù)鼗焖偎惴?，得到以速度、高度和迎角為自變量的氣?dòng)阻力系數(shù)數(shù)據(jù)表。

2）利用某一弧段的外測(cè)軌道星歷數(shù)據(jù)，擬合得到等效迎角。

3）把歷史弧段的等效迎角直接作為預(yù)報(bào)弧段的等效迎角；或通過(guò)連續(xù)幾個(gè)歷史弧段的等效迎角，優(yōu)化擬合得到預(yù)報(bào)弧段的等效迎角。

4）基于預(yù)報(bào)弧段的等效迎角計(jì)算獲取航天器氣動(dòng)力特性與飛行航跡，開(kāi)展軌道預(yù)報(bào)。

等效迎角的擬合確定，作為初步研究，設(shè)計(jì)了如下迭代過(guò)程：

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

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

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

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

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

4 計(jì)算結(jié)果及分析

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

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

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

4.2 大型航天器B無(wú)控隕落軌道衰降計(jì)算分析

4.2.1 長(zhǎng)弧段預(yù)報(bào)

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

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

圖3 大型航天器B無(wú)控飛行軌道衰降長(zhǎng)弧段預(yù)報(bào)結(jié)果（預(yù)報(bào)弧段為1個(gè)月，等效迎角擬合弧段為預(yù)報(bào)弧段的前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）

等效迎角反映了失穩(wěn)自旋航天器在擬合弧段的平均阻力姿態(tài)角大小。當(dāng)直接應(yīng)用到預(yù)報(bào)弧段時(shí)，隱含的假設(shè)是預(yù)報(bào)弧段的航天器無(wú)控飛行姿態(tài)與擬合弧段的飛行姿態(tài)具有繼承連續(xù)性，彼此相差不大。圖3中軌道高度范圍（350～300 km）變化不大，外界擾動(dòng)因素有可能沒(méi)有發(fā)生明顯變化，使得此時(shí)的無(wú)控航天器自旋運(yùn)動(dòng)模式和飛行姿態(tài)沒(méi)有發(fā)生明顯改變。

4.2.2 短弧段預(yù)報(bào)

當(dāng)利用等效迎角氣動(dòng)融合軌道直接積分建模方法開(kāi)展短弧段預(yù)報(bào)時(shí)，擬合弧段需要保證較合適的長(zhǎng)度，以反映等效迎角的變化。利用大型航天器B再入前5天的外測(cè)軌道星歷數(shù)據(jù)與3.2節(jié)介紹的方法，得到等效迎角的變化規(guī)律，如圖4所示，其中擬合弧長(zhǎng)為1天?？煽闯?，再入前5天等效迎角變化較大，從再入前倒數(shù)第5天的32.3°減小到再入當(dāng)天的20.9°。

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

圖5給出再入前倒數(shù)第3天軌道位置預(yù)報(bào)偏差，預(yù)報(bào)中采用的等效迎角通過(guò)再入前倒數(shù)第4天外測(cè)軌道星歷觀測(cè)數(shù)據(jù)擬合得到?？煽闯?，預(yù)報(bào)誤差隨著預(yù)報(bào)時(shí)間的增大而增大，預(yù)報(bào)5 h軌道位置誤差為1.3 km，預(yù)報(bào)24 h軌道位置誤差為19.6 km。利用軌道速度估算，預(yù)報(bào)24 h的時(shí)間偏差約為2.6 s。這說(shuō)明本文等效迎角氣動(dòng)融合軌道直接積分模型方法可以得到準(zhǔn)確性高的短弧段預(yù)報(bào)結(jié)果。

圖5 大型航天器B無(wú)控再入大氣層前倒數(shù)第3天軌道預(yù)報(bào)結(jié)果誤差（按倒數(shù)第4天外測(cè)軌道數(shù)據(jù)擬合等效迎角計(jì)算）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 結(jié)論

1）對(duì)大型航天器A受控隕落220～100 km過(guò)程飛行航跡的計(jì)算分析，驗(yàn)證了本文由空氣動(dòng)力精細(xì)數(shù)值模擬驅(qū)動(dòng)的低軌航天器氣動(dòng)特性快速算法以及軌道動(dòng)力學(xué)預(yù)報(bào)模型及方法的有效性。

2）在大型航天器B無(wú)控失穩(wěn)自旋飛行軌道衰降計(jì)算分析中，得到了與外測(cè)軌道星歷數(shù)據(jù)吻合一致的長(zhǎng)弧段和短弧段預(yù)報(bào)結(jié)果，在停止外測(cè)軌道數(shù)據(jù)供給5 h內(nèi)，地心慣性系位置預(yù)報(bào)偏差低于1.5 km，驗(yàn)證了基于等效迎角的氣動(dòng)融合軌道直接積分計(jì)算模型對(duì)無(wú)控航天器軌道衰降預(yù)報(bào)的合理性。

3）所建立基于等效迎角的氣動(dòng)融合軌道數(shù)值預(yù)報(bào)模型，是一種初步的方法。下一步可結(jié)合空氣動(dòng)力學(xué)中的氣動(dòng)參數(shù)辨識(shí)理論以及軌道動(dòng)力學(xué)中的反演方法，對(duì)等效迎角進(jìn)行更準(zhǔn)確合理的擬合優(yōu)化確定，為大型航天器再入解體分析提供更準(zhǔn)確的再入點(diǎn)數(shù)據(jù)。