基于機(jī)械能變化的超低軌道大氣密度實(shí)時(shí)估計(jì)

溫生林，閆野

(國(guó)防科學(xué)技術(shù)大學(xué) 航天科學(xué)與工程學(xué)院，湖南 長(zhǎng)沙410073)

0 引言

超低軌道是指大氣層以外而又低于一般航天器軌道高度的空間區(qū)域。隨著空間技術(shù)的發(fā)展和空間應(yīng)用的深入，超低軌道飛行技術(shù)將成為航天領(lǐng)域的新興發(fā)展方向。本文將超低軌道界定為120～300 km之間的飛行軌道。在此高度范圍內(nèi)，大氣非常稀薄，但相對(duì)于外層太空而言，又存在大量氣體分子，在衛(wèi)星高超聲速飛行時(shí)，氣體分子撞擊衛(wèi)星的固壁產(chǎn)生作用力。對(duì)于運(yùn)行于超低軌道上的衛(wèi)星，地球非球形引力和大氣阻力是主要的攝動(dòng)力，尤其是大氣密度相對(duì)傳統(tǒng)衛(wèi)星軌道而言顯著增加。研究發(fā)現(xiàn)，在強(qiáng)大氣阻力作用下，超低軌道會(huì)迅速衰減，大氣阻力是制約超低軌道飛行的最主要因素［1］。同時(shí)，這一空間區(qū)域的大氣密度隨太陽(yáng)活動(dòng)、季節(jié)、地磁場(chǎng)、光照條件等因素的改變發(fā)生劇烈變化，大氣密度存在較大的不確定性［2］。如何在軌實(shí)時(shí)獲取大氣密度是準(zhǔn)確計(jì)算大氣阻力、進(jìn)行超低軌道飛行控制的前提。

當(dāng)前，大氣密度的確定方法主要有三種途徑。一是基于經(jīng)驗(yàn)的大氣密度模型，常用的大氣模型有:指數(shù)模型、改進(jìn)的Harris-Priester模型、Jaccchia-71、Jaccchia-77、DTM、MSIS-90、JB2006 和 JB2008 模型［3］等。目前所使用的各種大氣密度模型都屬于半經(jīng)驗(yàn)公式，且模型誤差大約為10% ～20%，不能反映真實(shí)大氣密度隨太陽(yáng)活動(dòng)、地磁活動(dòng)等因素的影響。另一種途徑是利用最優(yōu)軌道確定理論對(duì)現(xiàn)有密度模型進(jìn)行修正［4-5］，其軌道確定是復(fù)雜的數(shù)學(xué)過(guò)程，計(jì)算量大，難以用于大氣密度的實(shí)時(shí)估計(jì)。此外，也可以利用星載加速度數(shù)據(jù)和相關(guān)測(cè)量數(shù)據(jù)來(lái)計(jì)算大氣密度［6］，這種方法除了要求進(jìn)行高精度的加速度測(cè)量外，還要實(shí)現(xiàn)高精度的姿態(tài)測(cè)量與穩(wěn)定控制，能夠較準(zhǔn)確地計(jì)算衛(wèi)星的阻力系數(shù)。

本文對(duì)超低軌道衛(wèi)星的大氣密度估計(jì)問(wèn)題進(jìn)行研究，利用軌道動(dòng)力學(xué)理論，分析了在地球引力、大氣阻力和控制力作用下超低軌道衛(wèi)星機(jī)械能的變化率，研究了如何利用能量變化實(shí)現(xiàn)大氣密度估計(jì)。

1 大氣密度估計(jì)方法

如圖1所示，OXIYIZI表示地心慣性系，oxyz表示軌道坐標(biāo)系，定義如下:

OXIYIZI:坐標(biāo)原點(diǎn)O位于地心，以赤道面為基準(zhǔn)面，OXI軸指向赤道平面上的平春分點(diǎn)，OZI軸垂直于平赤道面指向北極，OYI軸與OXI軸和OZI軸構(gòu)成右手系。

oxyz:坐標(biāo)原點(diǎn)o位于飛行器的質(zhì)心，ox軸由地心指向飛行器，oy軸在軌道平面內(nèi)垂直于ox軸，指向飛行器運(yùn)動(dòng)方向?yàn)檎?，oz軸與ox軸和oy軸構(gòu)成右手系。

圖1 坐標(biāo)系定義Fig.1 Coordinate Systems

根據(jù)軌道動(dòng)力學(xué)基本原理，衛(wèi)星的角動(dòng)量為:

軌道坐標(biāo)系的坐標(biāo)軸矢量可描述為:

軌道角速度ω= [ωxωyωz]T在軌道坐標(biāo)系中的分量為［7］:

軌道坐標(biāo)系的旋轉(zhuǎn)角速度為:

可得:

地心距矢量r在軌道坐標(biāo)系中可表示為:

對(duì)上式求導(dǎo)，可得衛(wèi)星的速度:

衛(wèi)星的加速度為:

另一方面，在地球引力、大氣阻力加速度和控制力加速度作用下，衛(wèi)星動(dòng)力學(xué)方程為:

式中:U為引力勢(shì)函數(shù);fcontrol為控制加速度，fcontrol=為大氣阻力。大氣阻力加速度在軌道坐標(biāo)系中的表達(dá)式為:

式中:ρ為大氣密度;CD為大氣阻力系數(shù);A為迎風(fēng)面積;m為衛(wèi)星質(zhì)量。

考慮中心引力和 J2項(xiàng)攝動(dòng)，U的具體表達(dá)式為:

引力勢(shì)函數(shù)的梯度為:

結(jié)合式 (8)～式(10)和式(13)可得:

角動(dòng)量的變化率為:

衛(wèi)星單位質(zhì)量的動(dòng)能為:

衛(wèi)星的勢(shì)能為U，根據(jù)式 (17)和式(11)，衛(wèi)星的機(jī)械能為:

機(jī)械能的變化率為:

聯(lián)合式(4)和式(15)，可得:

將式(14)、式(16)和式(20)代入式(19)，可得機(jī)械能的變化率為:

式(21)給出了在地球J2攝動(dòng)、大氣阻力和控制力作用下機(jī)械能變化率的計(jì)算公式，通過(guò)該式可知，機(jī)械能變化率與大氣阻力和控制力有關(guān)，與J2攝動(dòng)無(wú)關(guān)。事實(shí)上，地球引力是保守力，若沒(méi)有非保守力的影響，在高階地球引力模型作用下，衛(wèi)星的動(dòng)能和勢(shì)能之和也是不變的。因此，在地球引力、大氣阻力和控制力共同作用下，大氣密度的計(jì)算公式可表示為:

2 仿真分析

本文以軌道高度為200 km的超低軌道衛(wèi)星為例進(jìn)行分析，仿真條件為:

(1)初始時(shí)刻超低軌道衛(wèi)星的軌道根數(shù)為a=6 578.137 km，e=0，i=30°，Ω =0°，ω =0°，f=0°。軌道初始?xì)v元取為2012年6月1日12時(shí)0分0秒，仿真時(shí)間6 h，仿真步長(zhǎng)1 s。

(2)衛(wèi)星軌道計(jì)算采用STK軟件HPOP軌道計(jì)算模型，考慮70×70階EGM96引力場(chǎng)模型和大氣阻力攝動(dòng)，大氣阻力系數(shù)CD=2.2，衛(wèi)星質(zhì)量為500 kg，面質(zhì)比為0.02，采用Jaccchia-71大氣密度模型，并將STK中的大氣密度作為“真值”。

對(duì)超低軌道衛(wèi)星在地球引力和大氣阻力作用下利用衛(wèi)星機(jī)械能的變化對(duì)大氣密度進(jìn)行估計(jì)，仿真結(jié)果見(jiàn)圖2～圖4。

圖2 機(jī)械能變化曲線(xiàn)Fig.2 Mechanical energy changing curve

圖3 機(jī)械能變化率的變化曲線(xiàn)Fig.3 The rate of mechanical energy changing curve

圖4 大氣密度估計(jì)值和“真值”對(duì)比Fig.4 Comparison between estimated and“true”densities

圖2和圖3分別給出了機(jī)械能和機(jī)械能變化率隨時(shí)間的變化。由圖可知，在大氣阻力作用下，機(jī)械能不斷減小，且呈現(xiàn)出波動(dòng)性加速衰減的狀態(tài)。圖4表示大氣密度隨時(shí)間的變化，可以看出，大氣密度的估計(jì)值和“真值”相一致，利用這種方法可以有效地估計(jì)出大氣密度，與密度“真值”相比，估計(jì)得到的大氣密度最大相對(duì)誤差不超過(guò)5.89%。

3 結(jié)束語(yǔ)

本文研究了利用能量耗散率進(jìn)行超低軌道衛(wèi)星大氣密度的估計(jì)方法?；谲壍绖?dòng)力學(xué)理論，分析了在地球引力、大氣阻力和控制力作用下超低軌道衛(wèi)星的機(jī)械能的變化特性，給出了利用能量變化實(shí)現(xiàn)超低軌道大氣密度的估計(jì)公式。對(duì)典型超低軌道的大氣密度估計(jì)進(jìn)行了仿真，結(jié)果表明，這種估計(jì)方法是有效的，最大相對(duì)誤差不超過(guò)8.29%。該方法所需的信息量少，僅需要超低軌道衛(wèi)星的位置和速度信息，實(shí)現(xiàn)簡(jiǎn)單，對(duì)于超低軌道衛(wèi)星進(jìn)行軌道維持與姿態(tài)控制時(shí)需要獲取高精度的大氣阻力和氣動(dòng)力矩具有重要意義。

［1］ Atsushi N，Masanori H，Masayoshi U.The study of a super low altitude satellite［C］//Proceeding of 58th International Astronautical Congress.Fukuoka，Japan，2007.

［2］ George R G，Ronald JP，Paul JC.Requirements for accurate near-real time atmospheric density correction［C］//AIAA/AAS Astrodynamics Specialist Conference.Denver，CO，2000:1-17.

［3］ Vallado D A，F(xiàn)inkleman D.A critical assessment of satellite drag and atmospheric density modeling［C］//AIAA/AAS Astrodynamics Specialist Conference and Exhibit.Honolulu，USA，2008:1-29.

［4］ James R W.Optimal orbit determination［J］.Advances in the Astronautical Sciences，2002，112(2):1123-1134.

［5］ McLaughlin C A，Bieber B S.Neutral density determined from CHAMPprecision orbits［J］.Advances in the Astronautical Sciences，2008，119(7):167-186.

［6］ 陳光明，符養(yǎng)，薛震剛.利用星載加速度計(jì)數(shù)據(jù)反演高層大氣密度的方法［J］.解放軍理工大學(xué)學(xué)報(bào)，2010，11(3):371-376.

［7］ Chen Weiyue，Jing Wuxing.Differential equations of relative motion under the influence of J2perturbation and air drag［C］//AIAA SPACE 2010 Conference and Exposition.Anaheim，USA，2010:1-13.