基于星載高精度GPS觀測數(shù)據(jù)的大氣密度反演

苗娟， 任廷領(lǐng)， 龔建村*， 劉四清， 李志濤

1 中國科學(xué)院國家空間科學(xué)中心， 北京　100190 2 中國科學(xué)院大學(xué)， 北京　100049 3 宇航動力學(xué)國家重點實驗室， 西安　710043

?

基于星載高精度GPS觀測數(shù)據(jù)的大氣密度反演

苗娟1,2,3， 任廷領(lǐng)1， 龔建村1*， 劉四清1， 李志濤1

1 中國科學(xué)院國家空間科學(xué)中心， 北京100190 2 中國科學(xué)院大學(xué)， 北京100049 3 宇航動力學(xué)國家重點實驗室， 西安710043

星載高精度GPS觀測數(shù)據(jù)可提供衛(wèi)星速度和位置信息，而衛(wèi)星的運行軌跡又與所處位置的大氣密度緊密相關(guān)，因此可通過求解大氣阻力微分方程，由高精度GPS觀測數(shù)據(jù)反演出衛(wèi)星運行軌跡上的熱層大氣密度．本文從星載高精度GPS觀測數(shù)據(jù)出發(fā)，給出大氣密度的反演方法，以及平均平動參數(shù)nM、反彈道系數(shù)B兩個重要參數(shù)的解算過程，并以天宮一號為例，給出反演結(jié)果與天宮一號觀測數(shù)據(jù)的比對．結(jié)果表明，反演結(jié)果與觀測值符合很好，兩者的均方差在2012年1月1日、2月24日分別為8.6%和8.4%，說明利用星載GPS觀測數(shù)據(jù)反演大氣密度是有效、可行的，可成為今后獲取高精度大氣密度的一種方法．

GPS觀測數(shù)據(jù)； 平運動； 反彈道系數(shù)； 熱層大氣密

1　引言

對于低軌道航天器而言，大氣密度產(chǎn)生的阻力攝動是影響衛(wèi)星軌道的重要因素之一．大氣密度預(yù)報模型的準(zhǔn)確程度對低軌道航天器的精密定軌、軌道預(yù)報、變軌控制、壽命預(yù)測、飛船返回控制，及空間碎片的碰撞規(guī)避等方面均具有重要意義．計算表明，若大氣密度模型存在10%的誤差，神舟八號飛船軌道位置20天后的預(yù)報偏差將達800 km以上(劉舒蒔，2014)．隨著中國載人航天工程的深入開展，航天操作越來越復(fù)雜，工程要求越來越高，這對大氣密度預(yù)報精度也提出了更高要求．

從20世紀50年代后期開始，隨著導(dǎo)彈和航天活動對各種高層大氣模型的不斷要求，發(fā)展了多種半經(jīng)驗的高層大氣模型，形成了CIRA系列(Kallmann-Bijl et al., 1961)、Jacchia系列(Jacchia，1970)、DTM系列(Berger et al.，1998；Bruinsma et al.，2003)、MSIS系列(Hedin，1987，1991；Picone et al.，2002)等各種參考大氣模型，為空間活動提供服務(wù)．由于高層大氣的變化極其復(fù)雜，雖然幾十年來大氣模型在不斷地改進和發(fā)展，但一般情況下，模型仍存在15%～30%左右的誤差(Marcos et al.，1994；Rhoden et al.，2000)．因此，如何進一步提高大氣模型的預(yù)報精度成為模型研究者一直致力解決的問題．改進大氣模型預(yù)報精度主要有兩種方式：一是建立新模型，增加新的控制變量或采用新的擬合函數(shù)，如Bowman(2006；2008a, 2008b)采用新的太陽輻射指數(shù)和地磁指數(shù)建立了JB2006和JB2008模型、Lei等(2012)基于新自然正交函數(shù)分析法建立了400km高度大氣密度變化特征模型；二是在原有模型的基礎(chǔ)上進行修正，如美國空軍空間指揮中心的HASDM計劃(Storz and Bowman，2005)、苗娟和劉四清(2001)開展的基于監(jiān)測數(shù)據(jù)的大氣密度模型實時修正．但無論哪種改進方法，其基礎(chǔ)都是要獲取大量準(zhǔn)確的大氣密度數(shù)據(jù)．

雖然當(dāng)前已發(fā)射運行的衛(wèi)星眾多，但直接用于大氣密度探測的還比較匱乏．除了國際公開的CHAMP、GRACE(A，B)熱層密度數(shù)據(jù)外，目前公開渠道可獲取的熱層大氣密度數(shù)據(jù)較少．在大量衛(wèi)星兩行軌道根數(shù)(TLEs)的基礎(chǔ)上，Picone等(2006)提出了一種從TLEs數(shù)據(jù)快速得到大氣密度的方法，Lean等(2006)將該方法應(yīng)用于Starshine衛(wèi)星，Emmert等(2004，2008，2009)及任廷領(lǐng)等(2014)使用該方法研究大氣密度的長期變化．但從北美防控司令部(NORAD)公開的TLEs監(jiān)測數(shù)據(jù)來看，數(shù)據(jù)點比較稀疏，每天約1～4個監(jiān)測時間點，這在一定程度上制約了所提取密度的精度，而星載高精度GPS觀測數(shù)據(jù)可以彌補這一缺點．本文以天宮一號上搭載的GPS觀測數(shù)據(jù)為例，開展利用高精度GPS觀測數(shù)據(jù)的大氣密度反演，為后期載人軌道大氣密度模型預(yù)報精度的改善提供數(shù)據(jù)基礎(chǔ)，為今后獲取高精度大氣密度提供一種新方法．

2　基本原理和方法

熱層大氣雖然很稀薄，但相對于高速運行的航天器而言，產(chǎn)生的阻力不可小視，阻力會使航天器動能減少，導(dǎo)致高度下降，軌道收縮，發(fā)生軌道衰變．

衛(wèi)星在大氣中運動，所受阻力F為

(1)

(2)

其中，B=CDA/m，為衛(wèi)星的反彈道系數(shù)，m為衛(wèi)星質(zhì)量．組成該系數(shù)的三個量CD，A，m在衛(wèi)星運行過程中都是變化的，尤其是阻力系數(shù)CD，它不僅與衛(wèi)星的運行姿態(tài)、衛(wèi)星表面的物理性質(zhì)有關(guān)，還與大氣組成成分等因素有關(guān)(Cook，1965)，從而導(dǎo)致B非常難以確定．

大氣阻力對衛(wèi)星的影響主要表現(xiàn)在平均平動參數(shù)nM的變化上．衛(wèi)星軌道在大氣阻力的影響下逐漸衰減，使軌道半長軸a逐漸減小，平均平動nM逐漸增大．在消去地球的非球形引力攝動影響后，衛(wèi)星的半長軸在大氣阻力下的變化規(guī)律可用公式(3)表示(King-Hele，1987)：

(3)

其中，μ=GME，為萬有引力常數(shù)與地球質(zhì)量的乘積．

(4)

其中，fw為無量綱參數(shù)，表示風(fēng)因素的影響，可近似為

(5)

其中，r為衛(wèi)星的地心距，ω=7.292×10-5rad·s-1為地球自轉(zhuǎn)角速度，i為衛(wèi)星軌道傾角．與利用水平風(fēng)場模型HWM-93的計算值相比，該近似所產(chǎn)生的誤差小于3%(Hedin，1996)，而且在對(5)式積分過程中，二者之間的偏差會減小，因此(5)式是對fw很好的近似．

(6)

(7)

(8)

得到λ后，即可得衛(wèi)星相應(yīng)位置上的大氣密度．

3　大氣密度反演過程與結(jié)果

3.1平均平動nM的計算

在利用衛(wèi)星軌道參數(shù)反演大氣密度過程中，平均平動nM是至關(guān)重要的參數(shù)，它消除了地球非球形引力和大氣阻力的周期項攝動及其他攝動力的影響，反映了衛(wèi)星軌道在大氣阻力攝動影響下的長期變化．在沒有人為進行軌道制動及其他特殊情況下，nM隨時間遞增，這就保證了公式(4)式兩側(cè)均為正值，從而確保所求的大氣密度為正值．

需要注意的是，由星載高精度GPS觀測數(shù)據(jù)通過初軌計算可以很容易地求出軌道的瞬時半長軸a，然后再求得瞬時平動n，但該值并不能用來替代nM．雖然二者相差很小，但表示的物理含義卻大不相同．瞬時平動n包含了地球非球形引力和大氣阻力的周期項攝動影響以及其他攝動力的影響，無法直接體現(xiàn)出大氣阻力的影響，并且在積分區(qū)間內(nèi)可能產(chǎn)生隨時間減小的情況，因此需要把瞬時值n換算到平均平動nM，具體方法如下：

首先利用高精度GPS觀測數(shù)據(jù)計算出慣性系下衛(wèi)星的位置矢量r和速度矢量v，通過初軌計算得到軌道的瞬時六根數(shù)，包括傾角i，升交點赤經(jīng)Ω，離心率e，近地點角距ω，平近點角M，瞬時平動n；其次再將六根數(shù)轉(zhuǎn)換為可迭代計算的參數(shù)形式，轉(zhuǎn)換過程不考慮大氣模型的影響；最后采用公式(9)，經(jīng)過反復(fù)迭代，將軌道瞬時值轉(zhuǎn)換為平均值：

(9)

其中，y表示衛(wèi)星態(tài)函數(shù)y(r,v)；x表示所求的6個軌道根數(shù)；Cij=Z-1，Z滿足矩陣(10)；f表示SGP4模型．

(10)

將收斂精度設(shè)置為10-10，即可利用迭代法由態(tài)函數(shù)y(r,v)求出衛(wèi)星的平均軌道根數(shù)．以北美防空聯(lián)合司令部給出的Explore 8衛(wèi)星為例，選取2000年歷元為1.19565055時刻，將該時刻的軌道數(shù)據(jù)代入SGP4模型中，可得態(tài)矢量，然后將該態(tài)矢量用上述方法進行迭代，可得平均軌道根數(shù)，具體參數(shù)見表1，由表可知，采用迭代法計算所得的6個軌道根數(shù)與NORAD公布結(jié)果一致，滿足本文的精度需求，本文用此方法來計算平均軌道根數(shù)．

表1　NORAD所給軌道參數(shù)與迭代法得到參數(shù)誤差比較Table 1　Comparison of errors between iteration method and orbit parameters from NORAD

3.2反彈道系數(shù)的確定

除了平均平動nM之外，在利用上述方法反演大氣密度的過程中，反彈道系數(shù)B(B=CDA/m)是另一個非常重要的參數(shù)，它的精度決定著所求大氣密度的精度．CD和A是隨衛(wèi)星的運動變化的，CD不僅與衛(wèi)星的運行姿態(tài)和表面的物理性質(zhì)有關(guān)，還受大氣組成及狀態(tài)的影響.因此很難確定，但是可以通過計算求出一段時間內(nèi)的平均值，該值可近似認為是真實值．

Bowman等(2002)通過計算31年的衛(wèi)星數(shù)據(jù)，將所得反彈道系數(shù)進行平均，得到一個非常精確的觀測值．但是該方法需要一個高精度的經(jīng)驗大氣模型，而實際上，大氣密度隨時間在變化，阻力系數(shù)CD也會隨近地點高度的變化而變化，這就使所得B值偏離真實值.在本文中，采用Emmert等(2006)的公式，并使用另一種計算方式求解B值．

(11)

(12)

3.3大氣密度反演結(jié)果

天宮一號是我國第一個目標(biāo)飛行器和空間實驗室，于2011年9月29日發(fā)射，分別于2011年11月3日、2012年6月18日、2013年6月13日與神舟八號、神舟九號和神舟十號成功對接．在飛行任務(wù)期間，經(jīng)過多次姿態(tài)和軌道調(diào)控．圖3給出了天宮一號在2012年軌道高度的變化情況，可以看到在任務(wù)過程中進行了4次調(diào)軌．在密度反演過程中，為避免軌道高度突然變化所導(dǎo)致的密度反演異常情況，需選擇一段軌道相對平穩(wěn)的時間段．本文選擇2012年1—3月自主飛行時間段，開展高精度GPS觀測數(shù)據(jù)反演大氣密度研究，其中平均平動nM采用本文計算方法獲得，反彈道系數(shù)B值采用3.2節(jié)中計算得到的0.00939 m2·kg-1．

天宮一號搭載有大氣密度探測器，天宮一號飛行姿態(tài)的變化，使探測器的探測效果受到影響，在三軸穩(wěn)定對地期間，數(shù)據(jù)觀測質(zhì)量最好．為了進行反演結(jié)果與探測結(jié)果比對，本文選擇2012年1月1—12日和2月23日—3月4日兩個三軸穩(wěn)定時間段，采用文中的密度反演方法，利用天宮一號的高精度GPS觀測數(shù)據(jù)反演大氣密度．圖4、圖5分別給出了2012年1月1日、2月24日的反演結(jié)果，包括反演密度和探測值的比較，反演密度/模型密度，衛(wèi)星高度、緯度和經(jīng)度的變化，數(shù)據(jù)時間間隔為2 s．圖4中，由于2012年1月1日09∶56∶02—11∶37∶00GPS觀測數(shù)據(jù)的缺失，密度值、λ(ρGPS/ρMSIS)、高度、經(jīng)度、緯度出現(xiàn)了跳躍．

圖1　1999年Explore 8 反彈道系數(shù)變化Fig.1　Variation of B value for Explorer 8 in 1999

圖2　2012年天宮一號反彈道系數(shù)變化Fig.2　Variation of B value for Tiangong-1 in 2012

圖3　2012年天宮一號軌道變化Fig.3　Variation of Tiangong-1 orbit in 2012

圖4　2012年1月1日GSP觀測數(shù)據(jù)反演得到的大氣密度與天宮一號密度測量值比較Fig.4　Comparisons between GPS-derived density and observed density on 1 January in 2012

圖5　2012年2月24日GPS觀測數(shù)據(jù)反演得到的大氣密度與天宮一號密度測量值比較Fig.5　Comparisons between GPS-derived density and observed density on 24 February in 2012

從圖4、圖5中的對比結(jié)果可以看到，基于GPS觀測數(shù)據(jù)反演的大氣密度(Gpsden表示)與探測結(jié)果(Obsden表示)均符合很好，通過誤差計算，1月1日、2月24日的反演密度與探測值之間的均方差分別為8.6%和8.4%，遠小于大氣模型的平均誤差．自主運行期間，受大氣密度拖曳的影響，1月1日天宮一號高度處于357～380 km范圍內(nèi)，2月24日天宮一號高度下降至345～370 km．通過圖4、圖5中λ的變化，可以看到基于GPS觀測數(shù)據(jù)反演的大氣密度與模型計算結(jié)果的比值在1上下波動，這也從另一個方面驗證了反演結(jié)果的合理性．λ>1時，說明天宮一號所在位置的當(dāng)?shù)孛芏纫哂谀Ｐ陀嬎阒?，?1則反之．

4　結(jié)論與討論

無論研究大氣密度的變化規(guī)律，還是提高大氣模型的預(yù)報精度，其基礎(chǔ)首先是要有大量準(zhǔn)確的密度數(shù)據(jù)．在目前大氣密度直接探測數(shù)據(jù)相對匱乏的情況下，利用星載高精度GPS觀測數(shù)據(jù)反演大氣密度，可以為今后獲取高精度大氣密度數(shù)據(jù)提供一種新途徑．本文以天宮一號為例，開展利用高精度GPS觀測數(shù)據(jù)反演大氣密度研究，結(jié)論如下：

(1) 利用GPS觀測數(shù)據(jù)反演大氣密度過程中，平運動nM是至關(guān)重要的參數(shù)，本文利用迭代法得到高精度的平均軌道根數(shù)，其中傾角i、升交點赤經(jīng)Ω、離心率e、近地點角距ω、平近點角M與NORAD公布數(shù)據(jù)一致，平均平運動nM的誤差<10-7．

(2) 基于近地點高度相似的兩顆衛(wèi)星其B值修正量一致的假設(shè)，以Starshine I作為參考衛(wèi)星，計算得到1999年7—12月Explorer 8的反彈道系數(shù)與Emmert等的結(jié)果一致，天宮一號2012年1—3月的反彈道系數(shù)為0.00939 m2·kg-1．

(3) 通過2012年1月1日、2月24日的反演大氣密度與觀測結(jié)果的比較，可以看到，利用天宮一號的GPS觀測數(shù)據(jù)反演得到的大氣密度與觀測值符合很好，反演密度與觀測值之間的均方差分別為8.6%和8.4%，說明利用GPS觀測數(shù)據(jù)反演大氣密度是有效、可行的，可成為今后獲取高精度大氣密度的一種新途徑．

本文僅以天宮一號為例開展了大氣密度反演驗證．下一步，將針對多顆衛(wèi)星的GPS觀測數(shù)據(jù)，評估和分析不同太陽活動和地磁條件下的大氣密度反演結(jié)果．

致謝感謝審稿專家對本文提出的寶貴建議．

Berger C, Biancale R, Ill M, et al. 1998. Improvement of the empirical thermospheric model DTM: DTM94—a comparative review of various temporal variations and prospects in space geodesy applications.J.Geod., 72(3): 161-178.

Bowman B R. 2002. True satellite ballistic coefficient determination for HASDM.∥AIAA/AAS Astrodynamics Specialist Conference and Exhibit. California, Monterey: AIAA, AAS.

Bowman B R. 2006. A new empirical thermospheric density model JB2006 using new solar indices.∥ AIAA/AAS Astrodynamics Specialist Conference and Exhibit 21-24 August 2006. Keystone, Colorado: AIAA, AAS. Bowman B R, Tobiska W K, Marcos F A, et al. 2008. The JB2006 empirical thermospheric density model.J.Atmos.Solar-Terr.Phys., 70(5): 774-793.

Bowman B R. 2008. A new empirical thermospheric density model JB2008 using new solar and geomagnetic indices.∥ AIAA/AAS Astrodynamics Specialist Conference 18-21 August 2008. Honolulu, Hawaii: AIAA/AAS. Bruinsma S, Thuillier G, Barlier F. 2003. The DTM-2000 empirical thermosphere model with new data assimilation and constraints at lower boundary: Accuracy and properties.J.Atmos.Solar-Terr.Phys., 65(9): 1053-1070.

Cook G E. 1965. Satellite drag coefficients.Planet.SpaceSci., 13(10): 926-964.

Emmert J T, Picone J M, Lean J L, et al. 2004. Global change in the thermosphere: Compelling evidence of a secular decrease in density.J.Geophys.Res., 109(A2): A02301, doi: 10.1029/2003JA010176.

Emmert J T, Meier R R, Picone J M, et al. 2006. Thermospheric density 2002—2004: TIMED/GUVI dayside limb observations and satellite drag.J.Geophys.Res., 111(A10): A10S16, doi: 10.1029/2005JA011495.

Emmert J T, Picone J M, Meier R R. 2008. Thermospheric global average density trends, 1967—2007, derived from orbits of 5000 near-Earth objects.Geophys.Res.Lett., 35(5): L05101, doi: 10.1029/2007GL032809. Emmert J T. 2009. A long-term data set of globally averaged thermospheric total mass density.J.Geophys.Res., 114(A6): A06315, doi: 10.1029/2009JA014102. Hedin A E. 1987. MSIS-86 thermospheric model.J.Geophys.Res., 92(A5): 4649-4662.

Hedin A E. 1991. Extension of the MSIS thermosphere model into the middle and lower atmosphere.J.Geophys.Res., 96(A2): 1159-1172.

Hedin A E, Fleming E L, Manson A H, et al. 1996. Empirical wind model for the upper, middle and lower atmosphere.J.Atmos.Terr.Phys., 58(13): 1421-1447.

Jacchia L G. 1970. New static models of the thermosphere and exosphere with empirical temperature profiles. Technical Report 313, Cambridge, MA: Smithsonian Astrophysical Observatory.

Kallmann-Bijl H, Boyd R L F, Lagow H, et al. 1961. CIRA 1961: COSPAR international reference atmosphere 1961. Amsterdam: North-Holland Publishing Company. King-Hele D. 1987. Satellite Orbits in an Atmosphere: Theory and Application. Glasgow: Blackie and Son Ltd., 291.

Lean J L, Picone J M, Emmert J T, et al. 2006. Thermospheric densities derived from spacecraft orbits: Application to the Starshine satellites.J.Geophys.Res., 111(A4): A04301, doi: 10.1029/2005JA011399. Lei J H, Matsuo T, Dou X K, et al. 2012. Annual and semiannual variations of thermospheric density: EOF analysis of CHAMP and GRACE data.J.Geophys.Res., 117(A1): A01310, doi: 10.1029/2011Ja017324.

Liu S H, 2014. The investigation of thermospheric density change during geomagnetic storm and model correction (in Chinese)[Ph. D. thesis]. Beijing: University of Chinese Academy Sciences. Marcos F A, Bass J N, Baker C R, et al. 1994. Neutral density models for aerospace applications.∥32nd Aerospace Sciences Meeting and Exhibit. Reno, NY: AIAA.

Miao J, Liu S Q, Li Z T, et al. 2001. Atmospheric density calibration using the real-time satellite observation.ChineseJ.SpaceSci. (in Chinese), 31(4): 459-466.

Picone J M, Hedin A E, Drob D P, et al. 2002. NRLMSISE-00 empirical model of the atmosphere: Statistical comparisons and scientific issues.J.Geophys.Res., 107(A12): SIA 15-1-SIA 15-16, doi: 10.1029/2002JA009430.

Picone J M, Emmert J T, Lean J L. 2005. Thermospheric densities derived from spacecraft orbits: Accurate processing of two-line element sets.J.Geophys.Res., 110(A3): A03301, doi: 10.1029/2004JA010585.

Ren T L, Miao J, Liu S Q, et al. 2014. Research on thermospheric densities derived from two-line element sets.ChineseJ.SpaceSci. (in Chinese), 34(4): 426-433.

Rhoden E A, Forbes J M, Marcos F A. 2000. The influence of geomagnetic and solar variabilities on lower thermosphere density.J.Atmos.Solar-Terr.Phys., 62(1): 999-1013.

Storz M F, Bowman B R, Branson M J I, et al. 2005. High accuracy satellite drag model (HASDM).Adv.SpaceRes., 36(12): 2497-2505.

附中文參考文獻

劉舒蒔. 2014. 暴時熱層大氣密度時空分布規(guī)律研究及模式修正[博士論文]. 北京: 中國科學(xué)院大學(xué).

苗娟, 劉四清, 李志濤等. 2001. 基于實時觀測數(shù)據(jù)的大氣密度模式修正. 空間科學(xué)學(xué)報, 31(4): 459-466.

任廷領(lǐng), 苗娟, 劉四清等. 2014. 利用衛(wèi)星兩行軌道根數(shù)反演熱層密度. 空間科學(xué)學(xué)報, 34(4): 426-433.

(本文編輯汪海英)

Thermospheric density derived from onboard GPS observation data

MIAO Juan1,2,3, REN Ting-Ling1, GONG Jian-Cun1*, LIU Si-Qing1, LI Zhi-Tao1

1NationalSpaceScienceCenter,ChineseAcademyofSciences,Beijing100190,China2UniversityofChineseAcademyofSciences,Beijing100049,China3StateKeyLaboratoryofAstronauticDynamics,Xi′an710043,China

Onboard GPS observation data of a satellite includes accurate information of velocity and location, which are closely related to atmospheric density. Therefore, these GPS information can be used to derive thermospheric density through integration of differential equation. This paper presents a new method of deriving atmospheric density with a high temporal resolution from precise orbit data of low earth orbiting (LEO) space objects, and also presents the solution procedure of mean motionnMand the inverse ballistic coefficient B which are the two most important parameters for retrieving density. Tiangong-1 is taken as an example to evaluate the effectiveness of the method. The result shows that the GPS-derived density is in good agreement with observed density with the average error 8.6% and 8.4% respectively on 1st January and 24th February in 2012 for Tiangong-1. This result indicates the method provides an effective and reliable way to obtain extensive and accurate thermospheric density.

GPS parameter; Mean orbit element; Inverse ballistic coefficient; Thermospheric density

10.6038/cjg20161004.

宇航動力學(xué)國家重點實驗室開放基金項目(2014ADL-DW0303，2016ADL-DW304)和國家863計劃(2015AA7033102B)聯(lián)合資助.

苗娟，女，1975年生，副研究員，主要從事中高層大氣和預(yù)報方法研究. E-mail：miaoj@nssc.ac.cn

龔建村，男，1961年生，研究員，博士生導(dǎo)師，主要從事空間物理、空間環(huán)境預(yù)報和效應(yīng)預(yù)測研究. E-mail: gongjc@nssc.ac.cn

10.6038/cjg20161004

P352

2014-10-19，2016-08-31收修定稿

苗娟， 任廷領(lǐng)， 龔建村等. 2016. 基于星載高精度GPS觀測數(shù)據(jù)的大氣密度反演. 地球物理學(xué)報，59(10):3566-3572,

Miao J, Ren T L, Gong J C, et al. 2016. Thermospheric density derived from onboard GPS observation data.ChineseJ.Geophys. (in Chinese),59(10):3566-3572,doi:10.6038/cjg20161004.