童科偉,邱亮,高朝輝,唐慶博,李揚,焉寧

(1.中國運載火箭技術(shù)研究院研究發(fā)展中心,北京100076;2.中國空間技術(shù)研究院,北京100081)

1 前言

航天器壽命末期離軌再入時會在大氣層內(nèi)燒蝕,但個別大型航天器如火箭末級、上面級、大型衛(wèi)星、空間站等無法完全燒毀,在無控再入地球時有可能對人類生命及財產(chǎn)安全造成巨大災(zāi)難。隨著未來高密度發(fā)射,在軌的航天器越來越多,有必要計算航天器再入預(yù)報與落點散布問題。

針對目前航天器壽命末期離軌再入預(yù)報及落點散布問題,計劃設(shè)計一種能夠考慮軌道攝動的高精度軌道/再入彈道計算分析軟件,并考慮基于蒙特卡羅仿真技術(shù)對標(biāo)準(zhǔn)再入彈道落點偏差進(jìn)行仿真計算,解決航天器再入軌跡計算以及再入落點散布分析問題。

本文的研究成果可應(yīng)用于在軌安全領(lǐng)域,對火箭末級、上面級離軌再入軌跡進(jìn)行預(yù)報,對人類生命財產(chǎn)安全防護(hù)具有重要意義。

2 算法設(shè)計

航天器軌道計算方法通常有解析法 (如擬平均根數(shù)法[1])和數(shù)值積分法 (如Cowell法[2,3])。隨著計算機技術(shù)的高速發(fā)展,計算精度更高的數(shù)值積分法得到了廣泛應(yīng)用。本文針對數(shù)值積分法用于標(biāo)準(zhǔn)再入軌道預(yù)報以及離軌再入落點偏差分析進(jìn)行討論。

再入軌道預(yù)報標(biāo)準(zhǔn)算法需要對軌道動力學(xué)模型進(jìn)行建模,還需要描述坐標(biāo)系統(tǒng)轉(zhuǎn)換、時間系統(tǒng)等,坐標(biāo)系統(tǒng)轉(zhuǎn)換又需要對地球自轉(zhuǎn)與公轉(zhuǎn)運動(涉及到地心固連系到慣性系轉(zhuǎn)換)進(jìn)行建模。

再入軌道預(yù)報標(biāo)準(zhǔn)算法是偏差計算的基礎(chǔ),基于再入軌道預(yù)報標(biāo)準(zhǔn)算法,根據(jù)再入初始軌道參數(shù) (位置、速度),以及模型偏差分布,可以計算再入體的標(biāo)準(zhǔn)再入軌跡以及偏差再入軌跡,實現(xiàn)再入體的再入軌跡預(yù)報以及落點散布。

2.1 標(biāo)準(zhǔn)軌道預(yù)報器設(shè)計

2.1.1 地球運動

地球相對慣性空間運動的概貌如下:地球繞其自轉(zhuǎn)軸自轉(zhuǎn),大約一天自轉(zhuǎn)一周,自轉(zhuǎn)軸在慣性空間的指向表現(xiàn)為進(jìn)動和章動運動,且自轉(zhuǎn)軸還表現(xiàn)為極移運動。同時地球還繞太陽公轉(zhuǎn),公轉(zhuǎn)周期略大于一年,公轉(zhuǎn)軌道平面與赤道平面不重合,兩個平面的交線在慣性空間的指向幾乎保持不變,兩個平面的夾角約為23.4°,這個夾角稱為黃赤交角。圖1描述了地球自轉(zhuǎn)的歲差和章動影響因素。此外地球自轉(zhuǎn)軸并不總是與地球地理北極保持一致,稱為極移運動。歲差、章動和極移描述了地球的自轉(zhuǎn)運動[2]。

地球?qū)嶋H自轉(zhuǎn)軸并不總是與地球地理北極保持一致,天文學(xué)稱地球?qū)嶋H自轉(zhuǎn)軸為極軸。由于日月等三體引力對地球運動的攝動影響,極軸偏離地理北極稱為極移運動。表現(xiàn)為繞地球地理北極、周期為一年的螺旋線運動。

圖1 地球歲差和章動Fig.1 Earth precession and nutation

地球作為在慣性空間繞極軸自旋的自旋體,由于環(huán)境力矩的存在,主要是太陽和月球引力梯度力矩,使自旋軸表現(xiàn)為進(jìn)動和章動運動。進(jìn)動使得春分點不斷向西移動,周期為26000年,即黃道面和赤道面的交線每26000年在天球上自西向東旋轉(zhuǎn)一周。地球進(jìn)動運動通常稱為歲差。

由于地球赤道附近隆起,地球質(zhì)量分布不均勻,月球等其他三體引力導(dǎo)致地球自轉(zhuǎn)角動量軸與地球自轉(zhuǎn)軸不重合。地球自轉(zhuǎn)軸繞地球自轉(zhuǎn)角動量軸的運動,稱為地球自轉(zhuǎn)軸章動運動。相對于地球自轉(zhuǎn)軸的進(jìn)動,章動運動為短周期項運動,章動周期為18.6年。此外章動還具有半月的周期相。章動角最大不超過0.006°。地球自轉(zhuǎn)章動運動與月球運動密切相關(guān)。

2.1.2 坐標(biāo)系轉(zhuǎn)換

常用的坐標(biāo)系有地固坐標(biāo)系和地心慣性坐標(biāo)系 (J2000坐標(biāo)系)。

一般軌道積分采用的方法如Cowell法等,需要在地心慣性系 (ECI)下直接積分,常用的參考系為J2000平春分點平赤道慣性系 (天文學(xué)中常稱為FK5,目前的國際天球坐標(biāo)系ICRS定義與之重合[4])。而GPS等的觀測數(shù)據(jù)取自WGS84坐標(biāo)系 (據(jù)ITRF的介紹,基于GPS觀測數(shù)據(jù)定義的新版WGS84坐標(biāo)系,與國際大地參考坐標(biāo)系ITRS或地心固連坐標(biāo)系ECEF的誤差僅為分米級,認(rèn)為二者屬同一坐標(biāo)系)[4],而且地球引力場模型也是建立在ECEF坐標(biāo)系下的。因此如何從地心固連坐標(biāo)系到地心慣性系的轉(zhuǎn)化是航天器動力學(xué)的一個基本問題。目前常用的理論是國際天文聯(lián)合會 (IAU)提出的IAU1980模型、更精密的IAU2000和IAU2006模型等[4]。

IAU1980、IAU2000和 IAU2006的具體實現(xiàn)可參考SOFA[5](Standards of Fundamental Astronomy,基礎(chǔ)天文標(biāo)準(zhǔn)庫)。其中 IAU2000和IAU2006的精度更高。

對于章動模型,IAU1980理論采用具有106項的三角系數(shù)和來考慮其主要項,另外還有一個4個系數(shù)的Seidelmann修正項,以及85項系數(shù)和的行星作用修正項,其精度約在10-4角秒[5]。IAU2000模型基于天文學(xué)理論的發(fā)展和大量的觀測數(shù)據(jù),目前的精度已達(dá)10-7角秒,其公式包括1365項附加項,其中678項為日月章動系數(shù)項,687項為行星章動系數(shù)項[5]。

J2000坐標(biāo)轉(zhuǎn)換矩陣涉及到歲差、章動、格林尼治恒星時和極移。

在某歷元UTC時刻,地固系到J2000的轉(zhuǎn)換矩陣可寫成[4,5]:

(1)極移矩陣

極移量 (xp,yp)的求解為:(xp,yp)=(x,y)IERS+(Δx, Δy)tidal+(Δx, Δy)nutation

極移量主要是由IERS根據(jù)天文觀測給出的,每周都有新的觀測數(shù)據(jù)。此外,由于地球潮汐和章動的影響,會對極移有微小的修正 (Δx,Δy)tidal和 (Δx,Δy)nutation。(x,y)IERS由IERS給出的觀測數(shù)據(jù)計算求得, (Δx, Δy)tidal和 (Δx, Δy)nutation可由公式計算得到。

(2)自轉(zhuǎn)矩陣

地球自轉(zhuǎn)角θ的求解根據(jù)轉(zhuǎn)換方法的不同有不同的求解方式。

(3)歲差章動矩陣

其中,常值偏差矩陣B,歲差矩陣P(t)和章動矩陣N(t)如下:

章動量(Δψ, Δε)為:

上式中,各模型的計算方法可以參考SOFA[5]。

2.1.3 時間系統(tǒng)

一個平太陽日的長度為24小時,即86400s。實際上由于地球自轉(zhuǎn)角速度存在非常微小的變化,導(dǎo)致地球自轉(zhuǎn)360.985647°時經(jīng)歷時間存在非常微小的變化。因此平太陽日長度不僅存在短周期波動,而且存在長期項波動。盡管平太陽日差別很小,但每天的逐步積累將使得平太陽日與地球不再匹配,失去了協(xié)調(diào)地球公轉(zhuǎn)和自轉(zhuǎn)運動計時,指導(dǎo)人們?nèi)粘Ｉ畹哪康?。天文學(xué)家引入“跳秒”(Leap Second)來克服平太陽日誤差的長期項積累,以協(xié)調(diào)日常計時和地球自轉(zhuǎn)運動。

本文采用的時間系統(tǒng)主要有格林尼治恒星時、協(xié)調(diào)世界時、儒略日 (JD)、相對儒略日(MJD)、地心動力學(xué)時、歷書時等,這些時間系統(tǒng)的定義、相互關(guān)系具體可見 [1]。

2.1.4 軌道攝動模型設(shè)計

由于軌道涉及到的攝動因素很多,而考慮的攝動因素越多,計算越復(fù)雜,需要的計算量也越大,因此在進(jìn)行軌道計算時,需要根據(jù)任務(wù)的等級選擇合適的模型[2、6-7],軌道積分器則選擇ODE45[8]。本文主要考慮以下主要攝動因素:

地球非球形攝動:地球非球形攝動是軌道的主要攝動,應(yīng)根據(jù)特定軌道及不同軌道精度要求而采用不同的地球模型。地球非球形引力中包括帶諧項和田諧項,其中J2帶諧項是影響地球軌道最重要的項。本文采用70×70階次的EGM96模型。

大氣阻力攝動:大氣阻力攝動是近地軌道的另一種主要攝動,在再入軌道計算中必須考慮大氣阻力的影響。目前廣泛應(yīng)用的高精度大氣模型當(dāng)屬Jacchia系列和MSIS系列。Akins等[9]比較了MSIS大氣模型系列和Jacchia大氣模型系列對于定軌精度的影響,結(jié)果表明,最新的NRLMSISE 2000模型略好于Jacchia大氣模型,因此本文選用NRLMSISE 2000作為參照標(biāo)準(zhǔn)。

太陽光壓攝動和日月三體攝動:在精確軌道預(yù)報時一般需要考慮太陽光壓和日月三體攝動的影響,特別是長時間在軌飛行時。再入段時間較短,可忽略這兩種攝動力的影響。本文中計算太陽光壓攝動和日月三體攝動時太陽和月球的位置基于NASA的DE405模型獲得。

2.2 偏差軌道預(yù)報器設(shè)計

航天器軌道參數(shù)以及自身特征參數(shù)一般存在偏差,因此基于標(biāo)準(zhǔn)軌道計算再入落點無法得到統(tǒng)計學(xué)意義上的落點散布。為了定量分析航天器軌道參數(shù)偏差及自身特征參數(shù)偏差對最終落點散布的綜合影響,需以上文標(biāo)準(zhǔn)軌道預(yù)報算法為基礎(chǔ),借助統(tǒng)計學(xué)工具進(jìn)行分析。

本文以蒙特卡羅仿真方法對航天器再入落點偏差進(jìn)行分析。蒙特卡羅仿真 (Monte Carlo Simulation,MCS)方法,又稱隨機仿真法或統(tǒng)計仿真法[10]。其主要思想是在計算機上仿真實際概率過程,然后加以統(tǒng)計處理,具有思想新穎、直觀性強、簡便易用等優(yōu)點。

蒙特卡羅仿真法最重要的參數(shù)是仿真次數(shù)N。N與蒙特卡羅仿真法的仿真結(jié)果精度直接相關(guān)。一般來說,仿真次數(shù)越多,仿真結(jié)果精度越高。但仿真結(jié)果精度越高,仿真次數(shù)必然越大,相應(yīng)帶來的計算量就越大。因此,需要選擇適當(dāng)?shù)姆抡娲螖?shù),使得仿真結(jié)果滿足精度要求,又不至于造成計算量過大。

偏差軌道預(yù)報器蒙特卡羅初值選取的流程如下[10]:

對每個輸入變量的取值范圍按假定的概率密度函數(shù)以等概率 (取值為1/n)分為n個互不重疊的區(qū)間間隔;對每個輸入變量在每個間隔內(nèi)的取值按各自的概率密度分布隨機抽樣;再對x1的n個取值隨機地與x2的n個取值組成n個配對,然后這n個配對再與x3的n個取值隨機組合,以此類推,可得一組n個抽樣的k維隨機變量組值;對輸入變量的隨機配對進(jìn)行有效篩選,選擇出一個適當(dāng)?shù)拈g隔配合。

把以上初值分別帶入標(biāo)準(zhǔn)軌道預(yù)報器進(jìn)行仿真即可以得到對應(yīng)的偏差,對落點偏差可以進(jìn)行常用的統(tǒng)計分析,得到落點偏差均值、方差等特征參數(shù)。

3 軌道預(yù)報器測試

以下分別以在軌飛行以及離軌再入兩種情況作為標(biāo)準(zhǔn)工況對本文的算法進(jìn)行測試。

3.1 軌道預(yù)報器測試 (在軌飛行)

以400km近地軌道飛行為例對本文的軌道預(yù)報器算法進(jìn)行測試。軌道計算時間2天。積分結(jié)果與STK/HPOP模型對比,二者位置和速度差別分別如圖2所示 (限于篇幅僅給出X軸誤差):

由計算結(jié)果可見,位置誤差在兩天內(nèi)量級不超過10m,速度誤差在兩天內(nèi)不超過0.01m/s。

經(jīng)分析,發(fā)現(xiàn)二者的誤差主要是由大氣模型和太陽光壓模型不一致引起,其中大氣模型誤差造成的影響更大一些,但是這些誤差不影響最終軌道預(yù)報器的結(jié)果。本文的結(jié)果與文獻(xiàn) [6、7]的結(jié)果也吻合。因此驗證了本文方法的計算可靠性。

3.2 軌道預(yù)報器測試 (再入)

下面以100km軌道再入為例對本文的軌道預(yù)報器算法進(jìn)行測試。軌道積分結(jié)果與STK/HPOP模型對比,二者位置和速度的差別分別如圖3所示 (限于篇幅僅給出最大的Y軸誤差):

圖2 在軌飛行測試X軸位置誤差 在軌飛行測試X軸速度誤差Fig.2 Position and velocity error of X axis for on-orbit flight test

圖3 再入測試Y軸位置誤差 再入測試Y軸速度誤差Fig.3 Position and velocity error of Y axis for reentry flight test

由計算結(jié)果可見,在再入過程中,X、Y、Z向位置誤差在再入末段量級不超過500m,速度誤差在再入末段不超過2.5m/s。經(jīng)分析發(fā)現(xiàn)這主要是由于大氣模型不一致造成?？紤]到在再入過程中,各種誤差因素 (如燒蝕、解體等)的綜合影響遠(yuǎn)大于大氣模型的影響,可以認(rèn)為本軟件的計算結(jié)果與STK基本一致,因此驗證了本軟件的計算可靠性。

4 再入落點偏差分析案例

現(xiàn)在以一個100km軌道的航天器再入實例為例進(jìn)行分析,其中初始軌道根數(shù)、航天器相關(guān)特征參數(shù)以及相關(guān)偏差參數(shù) (3σ值,假定呈正態(tài)分布)如表1所示。

表1 再入落點偏差設(shè)計案例參數(shù)表Tab.1 Input Diagram for Reentry Footprint Dispersion Analysis

軌道開始初始時刻:

2010-09-22 00∶00∶46.0。

軟件計算100次蒙特卡羅仿真,計算的落點經(jīng)緯度散布如圖4所示。偏差分析計算得到的落點平均緯度為55.46deg,標(biāo)準(zhǔn)差為2.36deg,平均經(jīng)度為311.50deg,標(biāo)準(zhǔn)差為1.16deg;而基于標(biāo)準(zhǔn)軌道預(yù)報器計算得到的緯度為55.45deg,經(jīng)度為311.40deg。

由計算結(jié)果可見,給定初始軌道偏差及航天器特征參數(shù)偏差分布,基于軌道偏差散布計算方法,可以得到航天器再入落點散布統(tǒng)計特性。

5 結(jié)論

航天器壽命末期離軌再入時會在大氣層內(nèi)燒蝕,但個別大型火箭末級、上面級、大型衛(wèi)星、空間站等無法完全燒毀,無控再入地球時有可能對人類生命及財產(chǎn)安全造成巨大災(zāi)難。未來高密度的發(fā)射,使在軌的航天器越來越多,因此有必要計算航天器再入預(yù)報與落點散布問題。本文針對目前在軌航天器壽命末期離軌再入預(yù)報及落點散布問題,設(shè)計了一種能夠考慮軌道攝動的高精度軌道/再入彈道計算分析軟件,并考慮基于蒙特卡羅仿真技術(shù)對標(biāo)準(zhǔn)再入彈道落點偏差進(jìn)行仿真計算,克服了STK等商業(yè)軟件無法計算落點散布的缺點,可解決傳統(tǒng)火箭末級、上面級再入軌跡計算以及再入落點散布分析問題。本文的研究成果可用于在軌安全領(lǐng)域,對航天器離軌再入進(jìn)行預(yù)報,對人類生命財產(chǎn)安全防護(hù)具有工程意義。

圖4 再入落點偏差分析落點經(jīng)緯度散布Fig.4 Longitude and latitude dispersion for reentry footprint analysis