基于機械能變化的超低軌道大氣密度實時估計

溫生林，閆野

(國防科學技術大學 航天科學與工程學院，湖南 長沙410073)

0 引言

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

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

本文對超低軌道衛星的大氣密度估計問題進行研究，利用軌道動力學理論，分析了在地球引力、大氣阻力和控制力作用下超低軌道衛星機械能的變化率，研究了如何利用能量變化實現大氣密度估計。

1 大氣密度估計方法

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

OXIYIZI:坐標原點O位于地心，以赤道面為基準面，OXI軸指向赤道平面上的平春分點，OZI軸垂直于平赤道面指向北極，OYI軸與OXI軸和OZI軸構成右手系。

oxyz:坐標原點o位于飛行器的質心，ox軸由地心指向飛行器，oy軸在軌道平面內垂直于ox軸，指向飛行器運動方向為正，oz軸與ox軸和oy軸構成右手系。

圖1 坐標系定義Fig.1 Coordinate Systems

根據軌道動力學基本原理，衛星的角動量為:

軌道坐標系的坐標軸矢量可描述為:

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

軌道坐標系的旋轉角速度為:

可得:

地心距矢量r在軌道坐標系中可表示為:

對上式求導，可得衛星的速度:

衛星的加速度為:

另一方面，在地球引力、大氣阻力加速度和控制力加速度作用下，衛星動力學方程為:

式中:U為引力勢函數;fcontrol為控制加速度，fcontrol=為大氣阻力。大氣阻力加速度在軌道坐標系中的表達式為:

式中:ρ為大氣密度;CD為大氣阻力系數;A為迎風面積;m為衛星質量。

考慮中心引力和 J2項攝動，U的具體表達式為:

引力勢函數的梯度為:

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

角動量的變化率為:

衛星單位質量的動能為:

衛星的勢能為U，根據式 (17)和式(11)，衛星的機械能為:

機械能的變化率為:

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

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

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

2 仿真分析

本文以軌道高度為200 km的超低軌道衛星為例進行分析，仿真條件為:

(1)初始時刻超低軌道衛星的軌道根數為a=6 578.137 km，e=0，i=30°，Ω =0°，ω =0°，f=0°。軌道初始歷元取為2012年6月1日12時0分0秒，仿真時間6 h，仿真步長1 s。

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

對超低軌道衛星在地球引力和大氣阻力作用下利用衛星機械能的變化對大氣密度進行估計，仿真結果見圖2～圖4。

圖2 機械能變化曲線Fig.2 Mechanical energy changing curve

圖3 機械能變化率的變化曲線Fig.3 The rate of mechanical energy changing curve

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

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

3 結束語

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

［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，Finkleman 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］ 陳光明，符養，薛震剛.利用星載加速度計數據反演高層大氣密度的方法［J］.解放軍理工大學學報，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.