基于GRACE星載加速度計(jì)數(shù)據(jù)的熱層密度反演

李文文, 李敏,2*, 施闖,2, 趙齊樂(lè),2

1 武漢大學(xué)衛(wèi)星導(dǎo)航定位技術(shù)研究中心， 武漢 430079 2 地球空間信息技術(shù)協(xié)同創(chuàng)新中心， 武漢 430079

?

基于GRACE星載加速度計(jì)數(shù)據(jù)的熱層密度反演

李文文1, 李敏1,2*, 施闖1,2, 趙齊樂(lè)1,2

1 武漢大學(xué)衛(wèi)星導(dǎo)航定位技術(shù)研究中心， 武漢 430079 2 地球空間信息技術(shù)協(xié)同創(chuàng)新中心， 武漢 430079

本文主要研究了利用GRACE星載加速度計(jì)數(shù)據(jù)反演熱層密度.首先聯(lián)合采用GRACE衛(wèi)星2007—2009年星載加速度計(jì)數(shù)據(jù)和星載GPS數(shù)據(jù)進(jìn)行動(dòng)力學(xué)定軌并同時(shí)估計(jì)加速度計(jì)校正參數(shù)，依此對(duì)加速度計(jì)數(shù)據(jù)進(jìn)行了校正.根據(jù)Sentman稀薄空氣動(dòng)力學(xué)方程計(jì)算衛(wèi)星空氣動(dòng)力系數(shù)，對(duì)校正后的加速度計(jì)數(shù)據(jù)進(jìn)行處理，反演得到了該時(shí)期沿軌熱層大氣密度.為分析反演密度的精度，將本文反演得到的GRACE-A衛(wèi)星沿軌密度與Doornbos的解算結(jié)果，以及經(jīng)驗(yàn)密度模型NRLMSISE00，HASDM模型進(jìn)行比較分析.統(tǒng)計(jì)結(jié)果表明，本文反演結(jié)果比Doornbos系統(tǒng)性偏大約5%～8%，二者間的標(biāo)準(zhǔn)差(STD, standard deviation)在10%以內(nèi)，具有較好的符合性.其差異主要是由于采用了不同的加速度標(biāo)校以及空氣動(dòng)力系數(shù)計(jì)算策略.本文反演得到的熱層密度較HASDM模型表現(xiàn)為正的系統(tǒng)性偏差且幅度在4%以內(nèi)，而Doornbos的結(jié)果較HASDM模型約系統(tǒng)性偏小4%～7%，二者與HASDM模型的標(biāo)準(zhǔn)差均為30%左右.另外，本文反演密度與NRLMSISE00模型之間存在約30%～40%的系統(tǒng)性bias，其STD也在30%左右.

熱層密度； GRACE； 加速度計(jì)； HASDM； NRLMSISE00； 空氣動(dòng)力

1 引言

星載加速度計(jì)數(shù)據(jù)包含了衛(wèi)星在軌所受非保守力信息，通過(guò)分離太陽(yáng)輻射光壓(Solar Radiation Pressure, SRP)、地球反照輻射等其他非保守力后可以獲得衛(wèi)星所受的空氣動(dòng)力，進(jìn)而可以用于熱層密度反演.隨著CHAMP、GRACE、GOCE等重力衛(wèi)星計(jì)劃的順利實(shí)施，星載加速度計(jì)資料跨度可達(dá)一個(gè)太陽(yáng)活動(dòng)周期，這為研究熱層密度提供了大量高精度、高時(shí)空分辨率的數(shù)據(jù).Bruinsma和Biancale(2003)，Bruinsma等(2004), Sutton等(2007), Doornbos等(2010), Doornbos(2012), 汪宏波和趙長(zhǎng)印(2009)研究了采用CHAMP、GRACE星載加速度計(jì)實(shí)測(cè)數(shù)據(jù)計(jì)算熱層大氣密度.這些成果極大地促進(jìn)了熱層密度對(duì)太陽(yáng)耀斑、太陽(yáng)風(fēng)暴和地磁暴的響應(yīng)(Bruinsma et al., 2006, 2010; Lei et al., 2011; Liu et al., 2011; Sutton et al., 2005, 2006)、熱層大尺度行擾(Bruinsma et al., 2007a, Bruinsma and Forbes, 2008, 2009; Sutton et al., 2006)、太陽(yáng)輻射變化等對(duì)熱層密度的影響(Guo et al., 2007)、熱層密度年內(nèi)周期性變化(Guo et al., 2008)等方面的研究.

利用CHAMP、GRACE星載加速度計(jì)數(shù)據(jù)反演熱層密度時(shí)首先要對(duì)加速度計(jì)數(shù)據(jù)進(jìn)行校正.通常認(rèn)為GFZ發(fā)布的CHAMP、GRACE衛(wèi)星的加速度計(jì)數(shù)據(jù)在每個(gè)方向均存在尺度誤差和偏差誤差(Case et al., 2010).Gerlach等(2003)，徐天河和楊元喜(2004) 研究了采用能量守恒法進(jìn)行加速度計(jì)校準(zhǔn)，該方法通過(guò)能量守恒方程并采用現(xiàn)有重力場(chǎng)模型估計(jì)標(biāo)校參數(shù)，但是嚴(yán)重依賴于先驗(yàn)重力場(chǎng)模型精度.Helleputte等(2008, 2009)在精密定軌(Precise Orbit Determination, POD)過(guò)程中以加速度計(jì)數(shù)據(jù)代替非保守力，通過(guò)估計(jì)加速度計(jì)標(biāo)校參數(shù)實(shí)現(xiàn)對(duì)CHAMP和GRACE衛(wèi)星加速度計(jì)數(shù)據(jù)的校正.其認(rèn)為如果加速度計(jì)校正參數(shù)精確確定，則聯(lián)合高精度的保守力模型以及校正后的加速度計(jì)數(shù)據(jù)進(jìn)行軌道積分后得到的軌道應(yīng)接近于“真實(shí)軌道”，因此加速度計(jì)校正的好壞可以通過(guò)軌道精度來(lái)體現(xiàn).由于采用該校正方法時(shí)尺度參數(shù)與偏差參數(shù)強(qiáng)反相關(guān)，van Helleputte等(2009)對(duì)尺度和偏差因子引入了強(qiáng)約束，另外為了補(bǔ)償力學(xué)模型誤差還估計(jì)了經(jīng)驗(yàn)力參數(shù).結(jié)果表明該方法定軌結(jié)果與JPL精密科學(xué)軌道(Precision Scientific Orbit, PSO)相比三維RMS(Root Mean Square)為3.5 cm，SLR檢核RMS為2.3 cm，且將尺度因子固定為常數(shù)時(shí)可估計(jì)獲得更為平滑的偏差序列.熊永清等(2011a)提出了利用CHAMP衛(wèi)星精密軌道數(shù)據(jù)作為偽觀測(cè)數(shù)據(jù)，通過(guò)軌道擬合的方式估計(jì)衛(wèi)星初始狀態(tài)向量和加速度計(jì)參數(shù).該方法能夠避免處理GPS數(shù)據(jù)的復(fù)雜性.為消弱尺度因子和線性偏差的相關(guān)性，熊永清等(2011b)在以上方法基礎(chǔ)上聯(lián)合多個(gè)單天法方程進(jìn)行二次校標(biāo)，其多天只解算一組共同的尺度因子，而線性偏差每個(gè)單天解算一組.

衛(wèi)星空氣動(dòng)力系數(shù)CA也是限制反演精度的重要因素之一.CA通常與衛(wèi)星的幾何形狀、表面材質(zhì)物理性質(zhì)(對(duì)大氣粒子的吸附、反射或其他交互作用性質(zhì))以及大氣狀態(tài)，如：大氣溫度，大氣中各種成分的濃度等因素有關(guān)，其可通過(guò)Cook模型(Cook 1965)、Sentman稀薄空氣動(dòng)力學(xué)方程(Sentman, 1961)等模型進(jìn)行計(jì)算.Cook模型中忽略了粒子隨機(jī)熱運(yùn)動(dòng)對(duì)衛(wèi)星的作用.Koppenwallner(2008)，Sutton (2009)的分析表明，對(duì)于CHAMP、GRACE這類衛(wèi)星，其與氣流方向平行的側(cè)面面板較為細(xì)長(zhǎng)，大氣粒子的熱運(yùn)動(dòng)會(huì)導(dǎo)致CA顯著增大.因此這種情況下采用Cook模型計(jì)算的空氣動(dòng)力系數(shù)偏低，會(huì)導(dǎo)致反演得到的熱層密度系統(tǒng)性偏大.Sentman方程(Sentman, 1961)考慮了入射粒子的不規(guī)則熱運(yùn)動(dòng)，并假設(shè)反射粒子全漫反射分布，被認(rèn)為是目前適合用于計(jì)算CHAMP、GRACE衛(wèi)星CA系數(shù)的模型(Sutton 2009).由于低軌衛(wèi)星在軌飛行中，衛(wèi)星表面會(huì)吸附一層氧原子，從而改變了衛(wèi)星表面材料對(duì)大氣粒子吸附、反射等性質(zhì)，Moe等(2004, 2005)在Sentman方程基礎(chǔ)上引入能量調(diào)節(jié)系數(shù)α.能量調(diào)節(jié)系數(shù)決定了粒子是否保持其平均熱能(α=0)或者獲得衛(wèi)星的表面溫度(α=1)，對(duì)CA有顯著影響.考慮到大氣中各種成分對(duì)CA的貢獻(xiàn)是非線性的，Doornbos (2012)在其研究中分別計(jì)算各種大氣成分對(duì)CA的影響并最終取其濃度加權(quán)平均值.

Bruinsma等(2003,2004)最早采用CHAMP衛(wèi)星加速度計(jì)數(shù)據(jù)研究熱層密度反演，其將校正后的加速度計(jì)數(shù)據(jù)扣除SRP等其他非保守力后得到“大氣阻力觀測(cè)值”，以該阻力觀測(cè)值和由經(jīng)驗(yàn)?zāi)Ｐ兔芏扔?jì)算的阻力的比值作為尺度因子，對(duì)模型密度進(jìn)行改正從而得到了“觀測(cè)密度”.在此基礎(chǔ)上，Bruinsma等(2004)分析了由太陽(yáng)、地磁活動(dòng)以及潮汐引起的熱層密度變化.Sutton等(2007)則系統(tǒng)研究了利用CHAMP、GRACE加速度數(shù)據(jù)反演熱層密度的方法，利用校正后的星固系X軸和Y軸加速度計(jì)數(shù)據(jù)解算熱層密度和大氣風(fēng)速，并分析了加速度計(jì)校準(zhǔn)誤差、大氣阻力/升力系數(shù)誤差、先驗(yàn)風(fēng)速模型誤差、SRP誤差等對(duì)反演結(jié)果精度的影響，其結(jié)果表明反演密度誤差一般小于15%，但是在地磁活動(dòng)劇烈期間會(huì)由于大氣風(fēng)速誤差增加而增至30%.Doornbos(2012)在以上研究基礎(chǔ)上，利用Angara軟件對(duì)衛(wèi)星表面力進(jìn)行精化建模，并分析了多種密度反演方式.國(guó)內(nèi)汪宏波和趙長(zhǎng)印(2008)處理了CHAMP衛(wèi)星2001年至2002年星載加速度計(jì)數(shù)據(jù)，并反演得到了沿軌熱層密度.其利用反演密度對(duì)經(jīng)驗(yàn)大氣模型DTM94和MSIS90的精度進(jìn)行了檢核，結(jié)果表明在太陽(yáng)活動(dòng)峰年DTM94模型誤差為30%～35%，MSIS90模型為25%～30%(汪宏波和趙長(zhǎng)印, 2008).

上述研究中，Bruinsma等(2003,2004)，Sutton等(2007), 汪宏波和趙長(zhǎng)印(2008)采用單天估計(jì)一組偏差和尺度因子的策略進(jìn)行加速度計(jì)校正.然而，根據(jù)van Helleputte等(2008,2009)的結(jié)果，利用該策略校正的加速度計(jì)數(shù)據(jù)進(jìn)行軌道積分的精度僅在10 cm水平，表明其標(biāo)校精度相對(duì)較低.Doornbos(2012)根據(jù)van Helleputte等(2009)的研究結(jié)果進(jìn)行加速度計(jì)校正，其校正時(shí)為了補(bǔ)償保守力模型誤差每600 s估計(jì)一組切向、法向和徑向經(jīng)驗(yàn)力參數(shù).然而，頻繁調(diào)節(jié)的經(jīng)驗(yàn)力參數(shù)同時(shí)還會(huì)吸收部分非保守力，從而導(dǎo)致校正參數(shù)可能出現(xiàn)偏差.Bruinsma等(2003,2004)，汪宏波和趙長(zhǎng)印(2008)在其研究中均采用Cook模型計(jì)算CHAMP衛(wèi)星空氣動(dòng)力系數(shù)，這可能導(dǎo)致其反演的密度偏大.

本文在以上基礎(chǔ)上研究了利用GRACE星載加速度計(jì)數(shù)據(jù)進(jìn)行熱層密度反演.在反演中，采用了與Van Helleputte等 (2009)相同的方式對(duì)GRACE衛(wèi)星加速度計(jì)進(jìn)行標(biāo)校，但是本文在不引入經(jīng)驗(yàn)力參數(shù)的前提下，通過(guò)調(diào)整定軌中的動(dòng)力學(xué)策略，達(dá)到了與前者一致的軌道精度水平，同時(shí)也保證了標(biāo)校參數(shù)的完整性.通過(guò)與已有結(jié)果和模型的比較，對(duì)本文反演密度的精度進(jìn)行了評(píng)價(jià).

2 數(shù)據(jù)源與模型

本文選擇GRACE衛(wèi)星2007—2009年L1b級(jí)的數(shù)據(jù)進(jìn)行處理，其包括星載GPS數(shù)據(jù)(GPS)，加速度計(jì)數(shù)據(jù)(ACC)，姿態(tài)四元數(shù)數(shù)據(jù)(SCA)，衛(wèi)星質(zhì)量數(shù)據(jù)(MAS)，衛(wèi)星機(jī)動(dòng)數(shù)據(jù)(THR)等.數(shù)據(jù)的具體說(shuō)明可參見(jiàn)文獻(xiàn)(Case et al.,2010).在進(jìn)行加速度標(biāo)校時(shí)，本文采用IGS提供的最終GPS精密星歷和30 s采樣率鐘差產(chǎn)品作為數(shù)據(jù)處理的時(shí)空基準(zhǔn)，其軌道、鐘差精度分別在2.5 cm(1D RMS)和75 ps(RMS)量級(jí)(http:∥igs.org/products).反演過(guò)程中，衛(wèi)星姿態(tài)通過(guò)恒星照相儀SCA數(shù)據(jù)計(jì)算，衛(wèi)星在軌質(zhì)量通過(guò)MAS數(shù)據(jù)獲取.為維持KBR儀器相互對(duì)準(zhǔn)，GRACE衛(wèi)星在軌運(yùn)行時(shí)經(jīng)常需要進(jìn)行機(jī)動(dòng)調(diào)整，由此引入的機(jī)動(dòng)力需在反演時(shí)扣除，否則會(huì)引起密度反演誤差.本文在處理時(shí)則依據(jù)THR數(shù)據(jù)剔除了機(jī)動(dòng)時(shí)期的數(shù)據(jù).

為了對(duì)反演密度的精度進(jìn)行評(píng)價(jià)，本文將處理結(jié)果與Doornbos、Sutton博士反演的結(jié)果進(jìn)行了比較.Doornbos、Sutton博士的數(shù)據(jù)成果均可以在線公開獲取，其具體的處理方法與策略可以參考Sutton等(2007)與Doornbos(2012).同時(shí)還對(duì)比了經(jīng)驗(yàn)密度模型，包括：NRLMSISE00模型(Picone et al., 2002)和HASDM模型(Bowman and Storz, 2003).其中，NRLMSISE00模型由美國(guó)海軍實(shí)驗(yàn)室在MSISE系列經(jīng)驗(yàn)密度模型的基礎(chǔ)上進(jìn)行改進(jìn)得到，是目前低軌空間物體軌道確定和預(yù)報(bào)，以及熱層-電離層耦合等空間天氣研究中常用的經(jīng)驗(yàn)?zāi)Ｐ?其誤差一般在10%～30%且隨高度增加而增大，在400～500 km高度處約在25%水平(Macros et al., 2006).HASDM模型則由美國(guó)空軍研究中心開發(fā)，其使用75～80個(gè)低軌空間物體的雷達(dá)跟蹤數(shù)據(jù)在確定軌道的同時(shí)估計(jì)大氣密度的改正值，并以球諧函數(shù)的形式每3 h估計(jì)一次Jacchia71模型中的拐點(diǎn)溫度和夜間散逸層最低溫度(Bowman and Storz, 2003).HASDM模型在200～800 km高度范圍的精度可達(dá)6%～8%(Storz et al., 2005).

3 加速度計(jì)校正

3.1 聯(lián)合GPS與加速度計(jì)數(shù)據(jù)的校正方法

GRACE衛(wèi)星上搭載的SuperSTAR星載加速度計(jì)能以極高的精度測(cè)量衛(wèi)星所受到的非保守力.SuperSTAR的X、Y軸(指向星固系下的Z、X軸)是高敏感軸，而Z軸(指向星固系Y軸)為低敏感軸，三軸的測(cè)量精度都認(rèn)為可以達(dá)到10-9m·s-2(Frommknecht et al., 2003).然而由于星上溫度變化，儀器老化以及儀器開關(guān)切換等原因，會(huì)導(dǎo)致加速度計(jì)測(cè)量的加速度數(shù)據(jù)包含了偏差和尺度變化.因此在使用加速度計(jì)數(shù)據(jù)時(shí)，首先要對(duì)其進(jìn)行校準(zhǔn).由于大氣阻力的作用，切向的加速度計(jì)校正參數(shù)比較容易準(zhǔn)確確定.法向與徑向的尺度與偏差參數(shù)與衛(wèi)星的初始狀態(tài)向量強(qiáng)相關(guān)，且法向是約束最弱的方向(van Helleputte et al., 2008).

GRACE衛(wèi)星搭載的GPS接收機(jī)能夠提供大量連續(xù)的高精度觀測(cè)數(shù)據(jù)，本文聯(lián)合星載GPS和加速度計(jì)數(shù)據(jù)進(jìn)行動(dòng)力學(xué)定軌，并同時(shí)估計(jì)加速度計(jì)校正參數(shù).定軌中動(dòng)力學(xué)模型策略主要參照GFZ最新發(fā)布的Level 2數(shù)據(jù)處理標(biāo)準(zhǔn)0005版本(Dahle et al., 2012)并進(jìn)行了一定調(diào)整.保守力方面選擇精度盡量高的模型，其中地球非球形引力攝動(dòng)采用EIGEN-6C模型計(jì)算，且靜態(tài)重力場(chǎng)截至130階，時(shí)變部分則截至50階，海潮模型采用FES2004模型并截至40階，固體潮、極潮、相對(duì)論等其他保守力則均采用IERS Conventions(2003)的推薦模型(McCarthy and Petit, 2004).非保守力模型則以加速度計(jì)數(shù)據(jù)替代.在使用加速度計(jì)數(shù)據(jù)時(shí)，先由GFZ發(fā)布的初始校正參數(shù)(Bettadpur, 2004)進(jìn)行預(yù)校正，然后在定軌過(guò)程中同時(shí)估計(jì)X、Y、Z三方向的偏差參數(shù)，尺度參數(shù)則不予估計(jì).這樣處理實(shí)際上是將尺度因子固定為常數(shù)，而僅估計(jì)偏差因子.其原因是一方面GRACE衛(wèi)星星上溫度控制較為穩(wěn)定，其尺度參數(shù)可以認(rèn)為是常數(shù)而不隨時(shí)間變化(Bruinsma et al., 2007)；另一方面，由于偏差參數(shù)和尺度參數(shù)強(qiáng)相關(guān)，估計(jì)的偏差參數(shù)會(huì)補(bǔ)償尺度參數(shù)先驗(yàn)值的誤差.Bezděk(2010)，熊永清等(2011b)的研究中也表明了這一現(xiàn)象.偏差參數(shù)一般采用分段常數(shù)模型估計(jì).本文測(cè)試了多組不同分段間隔情況下定軌的精度，發(fā)現(xiàn)估計(jì)間隔為60 min時(shí)軌道精度最佳，因此采用每60 min估計(jì)一組偏差參數(shù)的策略.這也與GFZ的數(shù)據(jù)處理標(biāo)準(zhǔn)(Dahle et al., 2012)中的策略一致.對(duì)于GPS觀測(cè)值，采用無(wú)電離層組合偽距(PC)和無(wú)電離層組合相位(LC)作為基本觀測(cè)模型，觀測(cè)誤差主要考慮GRACE和GPS衛(wèi)星天線相位偏心(phase center offset, PCO)以及GPS衛(wèi)星天線相位中心變化(phase center variation, PCV)、相對(duì)論效應(yīng)、相位纏繞誤差等.模糊度參數(shù)未進(jìn)行固定，而是直接采用浮點(diǎn)解.具體的數(shù)據(jù)處理策略見(jiàn)表1.3.2 校正結(jié)果評(píng)價(jià)

3.2.1 GRACE衛(wèi)星軌道精度

如前所述，利用校正后的加速度計(jì)數(shù)據(jù)結(jié)合保守力模型進(jìn)行軌道積分，所得軌道的精度可以反映加速度計(jì)校正的精度.依據(jù)表1中的數(shù)據(jù)處理策略，對(duì)GRACE衛(wèi)星2007—2009年數(shù)據(jù)進(jìn)行精密定軌處理，并將解算軌道與JPL PSO比較，統(tǒng)計(jì)每個(gè)定軌弧段的軌道差值RMS，其時(shí)序圖見(jiàn)圖1.由圖中可見(jiàn)，以JPL PSO為參考軌道，本文解算的GRACE-A/B衛(wèi)星軌道在切向(Along-track) 與法向(Cross-track) 的RMS均主要在1.0～2.0 cm范圍變化，而徑向(Radial) RMS序列變化相對(duì)較為穩(wěn)定，對(duì)于GRACE-A、B星均在1.5 cm左右變化.表2為三方向軌道差值的平均RMS，其中徑向平均RMS較法向與切向均略差，這可能是由于JPL PSO對(duì)GRACE衛(wèi)星的天線PCV進(jìn)行了改正，而本文并未考慮該影響.GRACE-A，B衛(wèi)星三年軌道差值的平均3D RMS分別為2.6 cm和2.3 cm.

SLR技術(shù)常作為外部檢核的方式用于軌道精度評(píng)價(jià).采用NPT格式SLR數(shù)據(jù)對(duì)本文解算的GRACE衛(wèi)星軌道進(jìn)行檢核，統(tǒng)計(jì)每個(gè)弧段SLR觀測(cè)值殘差(Observed Minus Computed, OMC)的平均值與標(biāo)準(zhǔn)差.統(tǒng)計(jì)時(shí)剔除高度角小于10°以及OMC絕對(duì)值大于50 cm的觀測(cè)值，最終所有剔除的觀測(cè)值約占總觀測(cè)數(shù)量的2.5%，所有弧段的SLR OMC平均值與標(biāo)準(zhǔn)差時(shí)間序列見(jiàn)圖2.可見(jiàn)OMC平均值主要在±3 cm以內(nèi)變化，未見(jiàn)明顯系統(tǒng)偏差.標(biāo)準(zhǔn)差(Standard Deviation, STD)則均主要在5 cm以內(nèi).表3為所有弧段SLR OMC平均值，STD和RMS的均值.同時(shí)對(duì)這三年JPL PSO軌道進(jìn)行了SLR檢核和精度統(tǒng)計(jì).對(duì)比可見(jiàn)，SLR檢核本文軌道的殘差平均值均小于2 mm，JPL PSO則存在約1 cm的系統(tǒng)偏差.二者SLR檢核RMS均在2.5 cm水平.

表1 利用GRACE星載GPS和加速度計(jì)數(shù)據(jù)進(jìn)行精密定軌 以及加速度計(jì)校正的數(shù)據(jù)處理策略Table 1 Data processing strategy for POD and accelerometer data calibration using GPS observations and accelerometer data

表2 GRACE-A/B衛(wèi)星與JPL PSO軌道 比較的平均RMS統(tǒng)計(jì)Table 2 RMS statistics of orbit differences of the GRACE-A/B satellite orbits w.r.t. JPL PSO orbits

表3 本文解算的GRACE-A/B衛(wèi)星軌道 以及JPL PSO軌道的SLR檢核精度統(tǒng)計(jì)Table 3 Mean, STD, RMS statistics of the calculated GRACE-A/B orbits from this paper and the JPL PSO

3.2.2 加速度計(jì)校正參數(shù)

圖3為2007—2009年間GRACE-A和GRACE-B衛(wèi)星加速度計(jì)偏差參數(shù)在星固系X,Y,Z方向的偏差時(shí)間序列.圖中個(gè)別弧段估計(jì)得到的偏差出現(xiàn)明顯跳變，其與衛(wèi)星控制器切換以及衛(wèi)星溫度突變等原因相關(guān)(Bruinsma et al., 2007; Doornbos, 2012).對(duì)于GRACE雙星而言三個(gè)方向的偏差參數(shù)變化都非常平緩，且存在隨時(shí)間的緩慢飄移，其主要由加速度計(jì)的老化導(dǎo)致，這也與Bruinsma等(2007)的結(jié)果一致.根據(jù)估計(jì)的偏差時(shí)序，對(duì)于GRACE-A/B衛(wèi)星其在X方向偏差參數(shù)的年均變化分別僅為1×10-8m·s-2和5×10-9m·s-2.而在Y方向的偏差參數(shù)年均變化相比X方向明顯較大，分別為9×10-8m·s-2和1.84×10-7m·s-2.

圖1 2007—2009年GRACE-A/B衛(wèi)星與JPL PSO軌道比較的RMS序列Fig.1 RMS series of orbit differences with respect to the JPL PSO during 2007—2009

圖2 2007—2009年GRACE-A/B衛(wèi)星SLR殘差平均值(a)與標(biāo)準(zhǔn)差(b)序列Fig.2 The mean, STD series of the SLR OMC of GRACE-A/B satellites during 2007—2009

圖3 GRACE-A/B衛(wèi)星2007—2009年加速度計(jì)偏差序列(A) GRACE-A衛(wèi)星； (B) GRACE-B衛(wèi)星.Fig.3 Time series of estimated bias parameters of GRACE-A/B satellites(A) GRACE-A satellite; (B) GRACE-B satellite.

4 基于加速度計(jì)數(shù)據(jù)的熱層密度反演方法與策略

本節(jié)介紹了利用X軸加速度計(jì)數(shù)據(jù)提取熱層密度的基本方法，并分析了Cook模型和Sentman模型用于計(jì)算空氣動(dòng)力系數(shù)的差異，以及其他非保守力對(duì)反演的影響.基于這些分析，本節(jié)確定了密度反演的具體數(shù)據(jù)處理策略.

4.1 利用X軸加速度計(jì)數(shù)據(jù)反演熱層密度

4.1.1 反演方法

(1)

式中,ρ為大氣密度，vr和vr分別為衛(wèi)星與大氣相對(duì)速度及其單位矢量，m為衛(wèi)星質(zhì)量，CA為整星總的空氣動(dòng)力系數(shù)，其在vr方向的分量一般稱為阻力系數(shù)CD，垂直于vr方向的分量則稱為升力系數(shù)CL，Aref這里是參考面積可以為任一給定正值，uL方向?yàn)樯Ψ较騿挝皇噶?，CD,i，CL,i系數(shù)分別為為面板i的大氣阻力和升力系數(shù).在計(jì)算時(shí)通常將每塊面板CD，CL系數(shù)除以Aref以使其無(wú)量綱.

GRACE衛(wèi)星姿態(tài)相對(duì)穩(wěn)定，其三軸旋轉(zhuǎn)歐拉角較小一般在1°以內(nèi).由于大氣阻力的作用，加速度計(jì)X軸分量信號(hào)強(qiáng)度大且校正精度較高(星固系X軸指向衛(wèi)星沿軌方向)，非常利于提取大氣密度.利用校正后的加速度計(jì)數(shù)據(jù)扣除SRP等非保守力并投影至X軸方向，得到大氣動(dòng)力在X軸的分量aobs,X：

(2)

(3)

其中，acal為校正后的加速度計(jì)數(shù)據(jù)(慣性系下)，asrp為太陽(yáng)輻射壓力，aalb為地球反照輻射壓力，aIR為地球紅外輻射壓力，eX為慣性系轉(zhuǎn)星固系旋轉(zhuǎn)矩陣中X方向的單位向量.根據(jù)大氣動(dòng)力、熱層密度和CA之間的函數(shù)關(guān)系(式1)反解熱層密度：

(4)

由式(2)和(4)可知，利用校正后的加速度計(jì)數(shù)據(jù)反演密度時(shí)，需要對(duì)除空氣動(dòng)力外的其他非保守力進(jìn)行準(zhǔn)確建模和扣除，并準(zhǔn)確計(jì)算衛(wèi)星的空氣動(dòng)力系數(shù).

4.1.2 空氣動(dòng)力系數(shù)模型

CD,CL系數(shù)反映了衛(wèi)星外表面與大氣粒子之間的交互作用的強(qiáng)弱，并與衛(wèi)星外表面材質(zhì)物理特性以及大氣粒子的成分等有關(guān).精確確定CD，CL系數(shù)是目前反演熱層密度的難點(diǎn).由于CHAMP、GRACE這類衛(wèi)星側(cè)面面板較長(zhǎng)，在Cook模型中被忽略的分子熱運(yùn)動(dòng)會(huì)對(duì)其空氣動(dòng)力系數(shù)有較大影響.Sutton(2009) 研究表明對(duì)于CHAMP衛(wèi)星，相比采用Cook模型，采用Sentman方程后反演的熱層密度與經(jīng)驗(yàn)密度之間的標(biāo)準(zhǔn)差降低了約36%，偏差也由約56.4%～61.3%降低至21.2%～23.3%.

根據(jù)Bettadpur(2012)給出的GRACE衛(wèi)星星體幾何模型，圖4顯示了GRACE-A衛(wèi)星面板與星固系示意圖.其中1, 2號(hào)面板為衛(wèi)星前后面，3, 4, 5和6號(hào)面板是衛(wèi)星的側(cè)面，7, 8號(hào)分別為衛(wèi)星底面和頂面，9號(hào)面板為衛(wèi)星圓柱狀底部的平面投影.對(duì)于9號(hào)面板因其法向方向不唯一(可沿衛(wèi)星X-Y平面任意方向)且面積較小，故而未參與整星空氣動(dòng)力系數(shù)的計(jì)算.圖5為分別利用Sentman方程和Cook模型計(jì)算的GRACE-A衛(wèi)星各個(gè)面板的CD系數(shù).由圖4中可見(jiàn)，由于1，8號(hào)面板外法向與相對(duì)速度方向夾角θ>90°，因此在Cook模型中其CD系數(shù)始終為0.2號(hào)面板外法向正對(duì)大氣氣流方向(θ接近0°)，其對(duì)衛(wèi)星的空氣動(dòng)力系數(shù)貢獻(xiàn)最大，Cook模型計(jì)算的結(jié)果與Sentman模型相差約為0.05.相差最大的面板為7號(hào)面板，Sentman模型計(jì)算的CD系數(shù)約高出Cook模型0.4～0.6，這是由于7號(hào)面板是GRACE衛(wèi)星面積最大的面板，兩種模型之間的差異通過(guò)面板面積被放大.對(duì)所有面板的空氣動(dòng)力系數(shù)求和，采用Cook模型和Sentman方程計(jì)算的整星空氣動(dòng)力系數(shù)分別為2.37和3.50，其差異可達(dá)32%(以Sentman模型結(jié)果為參考).

圖4 GRACE-A衛(wèi)星面板與星固系示意圖(Doornbos, 2012)Fig.4 Spacecraft body-fixed directions and panels of GRACE-A satellite

4.1.3 GRACE衛(wèi)星非保守力分析

本文采用Box-wing模型計(jì)算SRP.地球反照輻射壓力采用Bhanderi and Bak(2005)提出的模型計(jì)算，所需的地表反照輻射系數(shù)采用CERES網(wǎng)站提供的1°×1°格網(wǎng)分辨率的月平均地表反射系數(shù)，而地球長(zhǎng)波輻射則參考Rim(1992)給出的經(jīng)驗(yàn)?zāi)Ｐ?由于熱層密度與太陽(yáng)輻射強(qiáng)度相關(guān)，圖6分別給出了在太陽(yáng)活動(dòng)高年(2003年1月3日)和低年時(shí)期(2009年1月2日)GRACE-A衛(wèi)星所受的太陽(yáng)輻射壓力(SRP)，地球輻射壓力(ALB+IR)以及空氣動(dòng)力(AERO)的比較.在熱層反演中，最重要的是大氣阻力在X軸方向的信號(hào)強(qiáng)度.圖中可見(jiàn)，在太陽(yáng)活動(dòng)高年，由于熱層密度較大，其所受大氣阻力遠(yuǎn)大于太陽(yáng)輻射壓力，在X軸方向可以達(dá)到(0.8～3.2)×10-7m·s-2，比太陽(yáng)輻射壓力高一個(gè)量級(jí).而在太陽(yáng)活動(dòng)低年，大氣阻力在X方向稍大于太陽(yáng)輻射壓力，二者均在10-8m·s-2量級(jí).并且圖中可見(jiàn)在大氣阻力X分量極小處，SRP的X分量與其大小幾乎相等.因此，衛(wèi)星所受SRP需要精確計(jì)算.而對(duì)于地球輻射壓力，其量級(jí)通常較小在(1-2)×10-8m·s-2水平且主要作用于衛(wèi)星徑向(星固系Z軸方向)，因此在實(shí)際反演計(jì)算中可以忽略.

4.2 策略分析

根據(jù)4.1節(jié)可知，加速度計(jì)校正，地球輻射壓力，空氣動(dòng)力系數(shù)，衛(wèi)星質(zhì)量，大氣風(fēng)速等均能影響反演結(jié)果.大氣阻力主要作用于衛(wèi)星沿軌方向(即星固系X軸)，從而使得加速度計(jì)數(shù)據(jù)X軸分量信號(hào)強(qiáng)度大得多也更容易校正，這為反演熱層密度提供了有利條件.第3節(jié)中軌道精度的分析結(jié)果也驗(yàn)證了本文校正的結(jié)果具有較高精度.而地球輻射壓力主要作用于衛(wèi)星徑向，在X軸的分量要比大氣阻力低2～3量級(jí)，因此在反演中可以忽略.GRACE衛(wèi)星質(zhì)量年均變化1 kg左右，以GFZ公布的480 kg為參考質(zhì)量，1 kg質(zhì)量誤差引起的密度偏差約為0.2%.在反演時(shí)也可采用GFZ發(fā)布的精確的衛(wèi)星質(zhì)量數(shù)據(jù)(MAS數(shù)據(jù)).大氣風(fēng)速誤差(相對(duì)于旋轉(zhuǎn)大氣)對(duì)密度反演的影響比較復(fù)雜，其導(dǎo)致衛(wèi)星與大氣的相對(duì)速度在量級(jí)與方向上均會(huì)產(chǎn)生誤差，從而引起衛(wèi)星空氣動(dòng)力系數(shù)以及密度計(jì)算誤差.Bruinsma等(2004) 通過(guò)仿真分析表明大氣風(fēng)速與反演誤差幾乎呈線性關(guān)系，緯度分量在赤道地區(qū)每100 m·s-1誤差約引起4%的誤差，在兩極地區(qū)影響可以忽略，而經(jīng)線分量相反，在兩極地區(qū)每100 m·s-1誤差約引起5%的誤差，在赤道地區(qū)則較小為2%.大氣風(fēng)速一般在幾百米每秒量級(jí)，但是在空間環(huán)境擾動(dòng)等時(shí)期，大氣風(fēng)速可以增大至千米每秒量級(jí)，因此在反演時(shí)應(yīng)盡量減小其影響.

圖5 GRACE-A衛(wèi)星各面板大氣阻力系數(shù)(α=0.93, 2003年1月3日)Fig.5 Drag coefficients of each panel of GRACE-A satellite based on Cook and Sentman models on 2003 Jan-3rd (fixed α=0.93)

圖6 不同太陽(yáng)活動(dòng)強(qiáng)度時(shí)期GRACE衛(wèi)星空氣動(dòng)力與太陽(yáng)輻射壓力(A)、地球輻射壓力(B)比較Fig.6 Comparison between aerodynamic, SRP and Earth radiation during different solar activity periods

根據(jù)以上分析，本文按照表4中的策略反演熱層密度.其中，加速度計(jì)數(shù)據(jù)采用3.1節(jié)估計(jì)得到的每60 min一組的標(biāo)校參數(shù)進(jìn)行校正，SRP采用Box-wing模型計(jì)算，Albedo與IR均不予考慮.空氣動(dòng)力系數(shù)采用Sentman模型計(jì)算，因其對(duì)大氣溫度與衛(wèi)星面板溫度變化并不敏感，本文直接采用NRLMSISE00模型計(jì)算大氣溫度，衛(wèi)星表面溫度則采用Doornbos(2012)的推薦值273 K.調(diào)節(jié)系數(shù)是影響空氣動(dòng)力系數(shù)的重要參量，本文研究中采用Sutton博士的推薦值0.93.大氣風(fēng)速采用HWM93模型計(jì)算.此外，為避免衛(wèi)星姿態(tài)調(diào)整機(jī)動(dòng)時(shí)期引入的機(jī)動(dòng)力的影響，本文根據(jù)THR產(chǎn)品(包含衛(wèi)星機(jī)動(dòng)記錄數(shù)據(jù))剔除了衛(wèi)星點(diǎn)火機(jī)動(dòng)時(shí)期的數(shù)據(jù).

表4 熱層密度反演數(shù)據(jù)處理策略Table 4 Data processing strategy of thermosphere density derivation

5 熱層密度反演結(jié)果與統(tǒng)計(jì)分析

采用第4節(jié)提出的策略處理得到了GRACE-A和GRACE-B衛(wèi)星2007—2009年間的沿軌熱層密度.為了對(duì)反演結(jié)果的精度進(jìn)行評(píng)價(jià)，本節(jié)將反演得到的GRACE-A衛(wèi)星2007—2009年沿軌密度與Doornbos(2012)等的反演結(jié)果，以及NRLMSISE00，HASDM等經(jīng)驗(yàn)?zāi)Ｐ瓦M(jìn)行比較分析.

5.1 與Doornbos反演結(jié)果比較

圖7為Doornbos, Sutton博士反演的2007年1月28日GRACE-A衛(wèi)星沿軌密度與本文反演的結(jié)果進(jìn)行對(duì)比，其中ρobs為本文反演結(jié)果，Density_Direct為Doornbos采用單軸法反演結(jié)果，Density_Sutton為Sutton博士反演得到的結(jié)果.Density_Direct的反演策略與表4類似，其差異主要為：(1) 加速度計(jì)校正參數(shù)采用van Helleputte等(2009)中的結(jié)果；(2) 計(jì)算空氣動(dòng)力學(xué)系數(shù)時(shí)分別考慮各種氣體粒子單獨(dú)的貢獻(xiàn).Sutton博士的反演策略也與本文類似，但是其將反演結(jié)果根據(jù)衛(wèi)星軌道緯度進(jìn)行了合并(Sutton et al., 2007)，因此其數(shù)據(jù)點(diǎn)時(shí)刻并不規(guī)則.由圖中可見(jiàn)三者反演密度的變化趨勢(shì)以及幅度均基本一致.

將本文反演結(jié)果ρobs與Density_Direct進(jìn)行比較，對(duì)二者密度比值按年份進(jìn)行對(duì)數(shù)正態(tài)分布擬合并統(tǒng)計(jì)bias和標(biāo)準(zhǔn)差；為了便于對(duì)照分析，將Density_Direct與HASDM模型也進(jìn)行了比較，結(jié)果分別見(jiàn)圖8和表5.其中年均比較數(shù)據(jù)點(diǎn)數(shù)約為2400000，統(tǒng)計(jì)密度比值的采樣間隔設(shè)置為0.002.由圖8可見(jiàn)兩組密度比值擬合的相關(guān)系數(shù)均達(dá)到了0.95以上，且Density_Direct與HASDM密度比值的擬合相關(guān)系數(shù)稍大一些.Density_Direct與ρobs，Density_Direct與HASDM的密度比值均主要分布在0.5～2.0區(qū)間內(nèi)，但是Density_Direct與ρobs的比值分布更為集中，表明二者結(jié)果更為接近.由表5可見(jiàn)，Density_Direct較HASDM模型、ρobs均偏小約4%～8%.Density_Direct與HASDM模型密度比值的STD在30%左右，與ρobs之間STD則均小于10%，表明本文反演結(jié)果與Doornbos的結(jié)果之間符合性較好.這主要是由于二者采用了相同的數(shù)據(jù)進(jìn)行反演，且反演中采用的模型類似.

ρobs與Density_Direct之間的差異主要由計(jì)算空氣動(dòng)力系數(shù)與加速度計(jì)標(biāo)校的策略不一致所造成.Doornbos在計(jì)算CA時(shí)分別考慮了每種氣體粒子的貢獻(xiàn)，而本文在計(jì)算時(shí)則是直接采用大氣粒子的平均分子質(zhì)量.圖9為2007年1月1—15日半個(gè)月內(nèi)二者CA比值的頻數(shù)統(tǒng)計(jì)，其中統(tǒng)計(jì)的空氣動(dòng)力系數(shù)數(shù)量為129600，采樣間隔為0.005.對(duì)該頻數(shù)分布進(jìn)行對(duì)數(shù)正態(tài)擬合的平均值與標(biāo)準(zhǔn)差分別為0.977和1.016，擬合的相關(guān)系數(shù)為0.915.值得注意的是計(jì)算空氣動(dòng)力系數(shù)時(shí)所需要輸入的數(shù)據(jù)或模型參數(shù)均可以很容易保持一致，因此空氣動(dòng)力系數(shù)差異的變化不大，其造成的反演密度差異也應(yīng)在此水平.Doornbos(2012)通過(guò)仿真分析了加速度計(jì)校正誤差對(duì)反演密度的影響，表明當(dāng)X軸加速度計(jì)數(shù)據(jù)存在10 nm·s-2的偏差時(shí)，其造成的密度反演偏差和STD分別為-2.7%和1.5%.本文研究中將尺度參數(shù)固定為常數(shù)并每60 min估計(jì)一組偏差參數(shù)，而van Helleputte等(2009)則單天估計(jì)一組尺度與偏差參數(shù)且在切向、徑向與法向均引入了經(jīng)驗(yàn)力參數(shù)，二者標(biāo)校參數(shù)的差異對(duì)密度反演造成的影響值得進(jìn)一步研究.

圖7 Doornbos，Sutton等研究與本文反演結(jié)果對(duì)比Fig.7 Densities comparison of this paper with respect to results from Doornbos and Sutton

圖8 2007(A), 2008(B)和2009(C)年Density_Direct與HASDM模型以及ρ obs的密度比值的頻數(shù)統(tǒng)計(jì)以及對(duì)數(shù)正態(tài)分布擬合Fig.8 Frequency number count and log-normal fit of ratios of Density_Direct over ρobs and HASDM model densities during year 2007 (A), 2008 (B) and 2009 (C), respectively

圖9 Doornbos與本文計(jì)算的空氣動(dòng)力系數(shù)差異的頻數(shù)統(tǒng)計(jì)和對(duì)數(shù)正態(tài)分布擬合Fig.9 Frequency count and log-normal fit of the differences of the aerodynamic coefficients between Doornbos and this paper

圖10 GRACE-A衛(wèi)星2007年1月28日沿軌反演密度，NRLMSISE00和HASDM模型密度以及密度比值(a) 沿軌ρobs，NRLMSISE00和HASDM密度；(b) ρobs和NRLMSISE00分別與HASDM密度比值；(c) 星下點(diǎn)緯度; (d) 軌道高度； (e) 星下點(diǎn)地方時(shí).Fig.10 The accelerometer-derived, NRLMSISE00 and HASDM model densities of GRACE-A satellite on January 16, 2007, as well as their density ratios(a) Accelerometer-derived, NRLMSISE00 and HASDM model densities; (b) Ratio series of accelerometer-derived and NRLMSISE00 densities over HASDM model densities; (c) Orbit latitude; (d) Orbit altitude; (e) Local time of ground tracks.

圖11 2007(A), 2008(B)和2009(C)年ρ obs和 NRLMSISE00模型與HASDM模型密度比值的頻數(shù)統(tǒng)計(jì)以及對(duì)數(shù)正態(tài)分布擬合Fig.11 Frequency number count and log-normal fit of ratios of ρobs and NRLMSISE00 over HASDM model densities during year 2007 (A), 2008 (B) and 2009 (C) respectively

RatioofDensity_DirectoverHASDMRatioofDensity_Directoverρobs200720082009200720082009Bias-0.070-0.043-0.068-0.057-0.072-0.074STD1.2811.3281.3311.0741.0841.089

5.2 與NRLMSISE00和HASDM經(jīng)驗(yàn)?zāi)Ｐ捅容^

圖10為2007年1月28日GRACE-A衛(wèi)星沿軌的反演密度ρobs，NRLMSISE00和HASDM模型密度以及密度比值.該天為地磁平靜日，其全天Ap指數(shù)均小于10，平均Ap指數(shù)為5，F(xiàn)10.7指數(shù)和81天平均F10.7指數(shù)分別為79.0和84.5 sfu，表明太陽(yáng)活動(dòng)強(qiáng)度較低.為了便于比較，圖中還給出了衛(wèi)星軌道緯度，高度和當(dāng)?shù)貢r(shí).圖中可見(jiàn)軌道高度約在475 km左右.總體上ρobs與經(jīng)驗(yàn)?zāi)Ｐ兔芏鹊淖兓厔?shì)基本一致，但與HASDM更為接近，二者在密度峰值處大小非常接近.而NRLMSISE00模型預(yù)報(bào)密度較ρobs與HASDM均系統(tǒng)性偏大.ρobs與HASDM比值主要分布在0.5～1.5范圍內(nèi)，且在兩極地區(qū)時(shí)出現(xiàn)較大幅度的跳變，其原因可能是在極地地區(qū)大氣風(fēng)速很容易增大至千米每秒量級(jí)，采用HWM93模型計(jì)算的大氣風(fēng)速可能會(huì)存在較大的誤差，從而導(dǎo)致反演結(jié)果的誤差增大.圖 10中同時(shí)可見(jiàn)ρobs和NRLMSISE00模型密度在午夜時(shí)期均出現(xiàn)一個(gè)甚至多個(gè)峰值，這反映的是午夜密度極大(Midnight Density Maximum, MDM)現(xiàn)象.HASDM模型則未能反映出該變化.因此可以發(fā)現(xiàn)在午夜期間，ρobs、NRLMSISE00模型與HASDM模型的密度比值序列均存在峰值現(xiàn)象.

將本文反演的密度ρobs與NRLMSISE00模型預(yù)報(bào)密度以及HASDM模型比較并進(jìn)行統(tǒng)計(jì).圖11為2007—2009年密度比值的頻數(shù)分布統(tǒng)計(jì)，以及對(duì)數(shù)正態(tài)分布擬合.圖中可見(jiàn)，ρobs和NRLMSISE00模型以及HASDM模型之間擬合的相關(guān)系數(shù)均達(dá)到了0.99以上.表6為密度比值的對(duì)數(shù)正態(tài)分布擬合的bias與STD.可見(jiàn)，ρobs與HASDM模型之間bias較小，均在4%以內(nèi)，而與NRLMSISE00模型之間存在約30%～40%的年平均bias，其原因目前認(rèn)為主要是由于2007—2009年太陽(yáng)活動(dòng)異常寧?kù)o從而導(dǎo)致NRLMSISE00模型產(chǎn)生了尺度性bias(Emmert et al., 2010; Solomon et al., 2010; Liu et al.,2011).

表6 2007, 2008和2009年ρobs和NRLMSISE00與HASDM模型密度比值的對(duì)數(shù)正態(tài)分布擬合統(tǒng)計(jì)Table 6 Statistics of log-normal fit of ratios of ρobsand NRLMSISE00 over HASDM model densities during year 2007, 2008 and 2009 respectively

STD統(tǒng)計(jì)結(jié)果來(lái)看，ρobs與NRLMSISE00以及HASDM的STD均在30%左右，但是其與HASDM模型之間的STD均要稍小于NRLMSISE00模型.值得注意的是，ρobs較HASDM模型系統(tǒng)性偏大，而Doornbos反演結(jié)果系統(tǒng)性偏小約4%～7%，二者與HASDM模型之間的STD均在30%左右.

以上的統(tǒng)計(jì)分析表明本文采用加速度計(jì)反演的熱層密度與HASDM模型之間bias較小，而NRLMSISE00模型與本文反演密度以及HASDM均存在較大的bias.一般而言，不同密度模型或數(shù)據(jù)之間的bias主要來(lái)自于空氣動(dòng)力系數(shù)誤差，以及對(duì)熱層密度變化物理機(jī)制的錯(cuò)誤模型化等，因此其常常與高度、緯度、地方時(shí)、太陽(yáng)輻射和地磁活動(dòng)強(qiáng)度等因素有關(guān)(Doornbos et al., 2008; 陳旭杏等, 2013; Shi et al., 2015).本文利用加速度計(jì)數(shù)據(jù)反演密度時(shí)，分段加速度計(jì)校正參數(shù)的誤差雖然會(huì)導(dǎo)致密度反演出現(xiàn)bias，但是從長(zhǎng)時(shí)間來(lái)看加速度計(jì)標(biāo)校誤差是隨機(jī)的(Sutton et al., 2007)，而衛(wèi)星空氣動(dòng)力系數(shù)誤差則會(huì)導(dǎo)致顯著的bias.一般而言，能量調(diào)節(jié)系數(shù)對(duì)空氣動(dòng)力系數(shù)CA的影響最為明顯，調(diào)節(jié)系數(shù)越大CA越小.對(duì)于GRACE衛(wèi)星而言，其調(diào)節(jié)系數(shù)一般認(rèn)為在0.8以上，但是精確確定其調(diào)節(jié)系數(shù)目前仍然是一項(xiàng)挑戰(zhàn)(Sutton, 2009).類似的，對(duì)于NRLMSISE00，HASDM等經(jīng)驗(yàn)?zāi)Ｐ停涫褂玫男l(wèi)星阻力觀測(cè)數(shù)據(jù)等同樣存在由彈道系數(shù)誤差所導(dǎo)致的bias.

另一方面，ρobs與目前公認(rèn)精度最高的模型HASDM之間的標(biāo)準(zhǔn)差可達(dá)30%.其很大部分可能是由于HASDM采用的低階球諧模型僅能在相對(duì)較大時(shí)空尺度上對(duì)熱層密度的平均誤差進(jìn)行校正，而ρobs能夠反映熱層密度更高時(shí)空分辨率上的變化.Lechtenberg 等(2013) 將CHAMP和GRACE加速度計(jì)反演得到的密度與HASDM模型進(jìn)行比較時(shí)發(fā)現(xiàn)二者總體上非常符合，然而前者能夠反應(yīng)熱層密度變化的諸多高頻變化現(xiàn)象，如：午夜密度極大，地磁極尖區(qū)密度增加，大氣行擾等，而HASDM模型無(wú)法探測(cè)到以上變化.本文圖10中也顯示出了類似的現(xiàn)象.其次，本文所研究的數(shù)據(jù)時(shí)段為2007—2009年，屬于太陽(yáng)活動(dòng)低年時(shí)期，而Bowman和Storz(2003) 對(duì)HASDM模型精度進(jìn)行評(píng)價(jià)時(shí)采用的是2001—2002年的數(shù)據(jù)，屬于太陽(yáng)活動(dòng)高年.在太陽(yáng)活動(dòng)低年時(shí)期由于太陽(yáng)輻射較小導(dǎo)致熱層密度為整個(gè)太陽(yáng)活動(dòng)周期內(nèi)最低，因此大氣阻力信號(hào)強(qiáng)度較弱.在考慮其他非保守力如SRP等誤差不變的情況下，大氣阻力減小從而會(huì)使得反演熱層密度的相對(duì)誤差增大.而另一方面，HASDM模型采用低軌空間物體阻力觀測(cè)值來(lái)校正Jacchia71模型，其在太陽(yáng)活動(dòng)低年也同樣存在因阻力信號(hào)減弱所致的阻力觀測(cè)值相對(duì)誤差增大的問(wèn)題.這也是導(dǎo)致本文統(tǒng)計(jì)結(jié)果中ρobs與HASDM模型之間STD較大的原因之一.因此，采用整個(gè)太陽(yáng)活動(dòng)周期內(nèi)的反演密度數(shù)據(jù)將有利于更準(zhǔn)確地評(píng)價(jià)反演精度，特別是在不同太陽(yáng)輻射強(qiáng)度條件下的反演密度精度.另外，對(duì)于反演得到的結(jié)果需要進(jìn)行嚴(yán)格的質(zhì)量控制和粗差剔除.現(xiàn)階段本文僅僅根據(jù)THR產(chǎn)品對(duì)衛(wèi)星點(diǎn)火時(shí)期的數(shù)據(jù)進(jìn)行了刪除，然而對(duì)于極區(qū)大氣風(fēng)速誤差等導(dǎo)致的反演粗差仍然參與了比較，這也可能會(huì)對(duì)統(tǒng)計(jì)結(jié)果產(chǎn)生干擾.

6 總結(jié)

本文介紹了利用加速度計(jì)數(shù)據(jù)反演熱層密度的模型與方法，并對(duì)GRACE衛(wèi)星2007—2009年的數(shù)據(jù)進(jìn)行了處理與分析.首先利用星載GPS數(shù)據(jù)進(jìn)行動(dòng)力學(xué)精密定軌，并以加速度計(jì)數(shù)據(jù)代替非保守力，估計(jì)得到了該時(shí)段內(nèi)GRACE衛(wèi)星的加速度計(jì)標(biāo)校參數(shù).在此基礎(chǔ)上利用校正后的加速度計(jì)數(shù)據(jù)扣除太陽(yáng)輻射光壓，并根據(jù)Sentman方程計(jì)算衛(wèi)星空氣動(dòng)力系數(shù)，反演得到了2007—2009年沿GRACE衛(wèi)星軌道的熱層大氣密度.本文研究的主要結(jié)論如下：

(1) 對(duì)不同太陽(yáng)活動(dòng)強(qiáng)度時(shí)期的衛(wèi)星所受空氣動(dòng)力、太陽(yáng)輻射壓力，地球反照輻射等非保守力進(jìn)行比較，結(jié)果表明在太陽(yáng)活動(dòng)高年衛(wèi)星阻力信號(hào)較強(qiáng)，其在X軸分量比SRP高出一個(gè)量級(jí)；而太陽(yáng)活動(dòng)低年時(shí)期由于衛(wèi)星阻力信號(hào)變?nèi)?，可能?dǎo)致反演密度誤差增加.

(2) 利用GRACE衛(wèi)星2007—2009年的GPS與加速度計(jì)數(shù)據(jù)進(jìn)行動(dòng)力學(xué)定軌處理，通過(guò)軌道精度分析對(duì)加速度計(jì)校正精度進(jìn)行了評(píng)價(jià).結(jié)果表明，采用本文處理策略，利用GPS和加速度計(jì)數(shù)據(jù)定軌結(jié)果與JPL PSO在同一精度量級(jí)，GRACE-A/B星與JPL PSO軌道差值3D RMS平均值分別為2.6 cm和2.3 cm，SLR軌道檢核的RMS平均值分別為2.4 cm和2.6 cm.通過(guò)該方法估計(jì)得到的三個(gè)方向的偏差參數(shù)變化都非常平緩，X方向偏差參數(shù)的年均變化在10 nm·s-2量級(jí).

(3) 分析了利用GRACE衛(wèi)星加速度計(jì)進(jìn)行密度反演的數(shù)據(jù)處理策略，對(duì)反演結(jié)果進(jìn)行分析與評(píng)估.對(duì)GRACE衛(wèi)星2007—2009年的加速度計(jì)數(shù)據(jù)進(jìn)行處理得到了該時(shí)期內(nèi)沿GRACE衛(wèi)星軌道的熱層大氣密度.將本文反演的結(jié)果與Doornbos解算的結(jié)果，經(jīng)驗(yàn)密度模型NRLMSISE00和HASDM進(jìn)行對(duì)比分析.結(jié)果表明本文反演結(jié)果與Doornbos結(jié)果之間存小于8%的系統(tǒng)性bias，二者STD統(tǒng)計(jì)值在10%以內(nèi)，具有較好的符合性.二者之間的差異主要是由不同的加速度標(biāo)校和空氣動(dòng)力系數(shù)計(jì)算策略造成.本文反演結(jié)果較HASDM偏大，bias在4%以內(nèi)，而Doornbos結(jié)果較HASDM模型值偏小，bias在7%以內(nèi)，二者與HASDM模型比較的STD均在30%左右.本文結(jié)果與NRLMSISE00模型之間存在約30%～40%的系統(tǒng)bias，這主要是由于NRLMSISE00模型在2007—2009年太陽(yáng)活動(dòng)異常寧?kù)o期間存在顯著模型偏差所致.

致謝 感謝ISDC提供GRACE衛(wèi)星相關(guān)數(shù)據(jù)，以及IGS提供GPS精密軌道和鐘差產(chǎn)品.感謝Doornbos和Sutton博士提供基于GRACE衛(wèi)星加速度計(jì)反演的密度數(shù)據(jù).感謝各基金項(xiàng)目對(duì)論文研究的支持.

Bettadpur S V. 2004. Recommendation for a-priori bias & scale parameters for level-1B ACC data (Release 00). GRACE TN-04-02.

Bettadpur S V. 2012. GRACE Product Specification Document, GRACE 327-720. The University of Texas at Austin.

Bezděk A. 2010. Calibration of accelerometers aboard GRACE satellites by comparison with POD-based nongravitational accelerations.JournalofGeodynamics, 50(5): 410-423.

Bhanderi D D V, Bak T. 2005. Modeling Earth Albedo for Satellites in Earth Orbit.∥ AIAA Guidance, Navigation, and Control Conference and Exhibit. San Francisco, California: AIAA.

Bowman B R, Storz M F. 2003. High Accuracy Satellite Drag Model (HASDM) review.AdvancesintheAstronauticalSciences, 116: 1943-1952. Bruinsma S, Biancale R. 2003. Total densities derived from accelerometer data.JournalofSpacecraftandRockets, 40(2): 230-236. Bruinsma S, Tamagnan D, Biancale R. 2004. Atmospheric densities derived from CHAMP/STAR accelerometer observations.PlanetaryandSpaceScience, 52(4): 297-312. Bruinsma S, Forbes J M, Nerem R S, et al. 2006. Thermosphere density response to the 20-21 November 2003 solar and geomagnetic storm from CHAMP and GRACE accelerometer data.JournalofGeophysicalResearch:SpacePhysics, 111(A6): A06303.

Bruinsma S, Biancale R, Perosanz F. 2007. Calibration parameters of the CHAMP and GRACE accelerometers. Poster presented at the GRACE Science Team Meeting.

Bruinsma S L, Forbes J M. 2007. Global observation of traveling atmospheric disturbances (TADs) in the thermosphere.GeophysicalResearchLetters, 34(14): L14103. Bruinsma S L, Forbes J M. 2008. Medium- to large-scale density variability as observed by CHAMP.SpaceWeather, 6(8): S08002. Bruinsma S L, Forbes J M. 2009. Properties of traveling atmospheric disturbances (TADs) inferred from CHAMP accelerometer observations.AdvancesinSpaceResearch, 43(3): 369-376.

Bruinsma S L, Forbes J M. 2010. Anomalous behavior of the thermosphere during solar minimum observed by CHAMP and GRACE.JournalofGeophysicalResearch:SpacePhysics, 115(A11): A11323.

Case K, Kruizinga G L H, Wu S C. 2010. GRACE Level 1B Data Product User Handbook, JPL Publication D-22027.

Chen X X, Hu X, Xiao C Y, et al. 2013. Correction method of the low earth orbital neutral density prediction based on the satellites data and NRLMSISE-00 model.ChineseJ.Geophys. (in Chinese), 56(10): 3246-3254, doi: 10.6038/cjg20131003.

Cook G E. 1965. Satellite drag coefficients.PlanetaryandSpaceScience, 13(10): 929-946.

Dahle C, Flechtner F, Gruber C, et al. 2012. GFZ GRACE Level-2 Processing Standards Document for Level-2 Product Release 0005. Deutsches GeoForschungsZentrum GFZ.

Doornbos E, Klinkrad H, Visser P. 2008. Use of two-line element data for thermosphere neutral density model calibration.AdvancesinSpaceResearch, 41(7): 1115-1122.

Doornbos E, van Den Ijssel J, Luehr H, et al. 2010. Neutral density and crosswind determination from arbitrarily oriented multiaxis accelerometers on satellites.JournalofSpacecraftandRockets, 47(4): 580-589.

Doornbos E. 2012. Thermospheric Density and Wind Determination from Satellite Dynamics. Berlin Heidelberg: Springer.

Emmert J, Lean J, Picone J. 2010. Record-low thermospheric density during the 2008 solar minimum.GeophysicalResearchLetters, 37(12): 245-269.Frommknecht B, Oberndorfer H, Flechtner F, et al. 2003. Integrated sensor analysis for GRACE-development and validation.AdvancesinGeosciences, 1: 57-63. Gerlach C, Sneeuw N, Visser P, et al. 2003. CHAMP gravity field recovery with the energy balance approach: first results.∥ First CHAMP Mission Results for Gravity, Magnetic and Atmospheric Studies. Berlin Heidelberg: Springer, 134-139.

Guo J P, Wan W X, Forbes J M, et al. 2007. Effects of solar variability on thermosphere density from CHAMP accelerometer data.JournalofGeophysicalResearch:SpacePhysics, 112(A10): A10308.

Guo J P, Wan W X, Forbes J M, et al. 2008. Interannual and latitudinal variability of the thermosphere density annual harmonics.JournalofGeophysicalResearch:SpacePhysics, 113(A8): A08301.

Koppenwallner G. 2008. Comment on special section: New perspectives on the satellite drag environments of earth, mars, and venus.JournalofSpacecraftandRockets, 45(6): 1324-1327. Lechtenberg T F, McLaughlin C A, Locke T, et al. 2013. Thermospheric density variations: observability using precision satellite orbits and effects on orbit propagation.SpaceWeather, 11(1): 34-35.

Lei J H, Thayer J P, Lu G, et al. 2011. Rapid recovery of thermosphere density during the October 2003 geomagnetic storms.JournalofGeophysicalResearch:SpacePhysics, 116(A3): A03306.

Liu H X, Doornbos E, Yamamoto M, et al. 2011. Strong thermospheric cooling during the 2009 major stratosphere warming.GeophysicalResearchLetters, 38(12): L12102. Marcos F A, Bowman B R, Sheehan R E. 2006. Accuracy of Earth′s thermospheric neutral density models.∥ AIAA/AAS Astrodynamics Specialist Conference and Exhibit, Guidance, Navigation, and Control and Co-located Conferences. Keystone, CO.: AIAA.

McCarthy D D, Petit G. 2004. IERS conventions 2003. IERS Technical Note No. 32. Frankfurt am Main: Verlag des Bundesamts für Kartographie und Geod?sie.

Moe K, Moe M M, Rice C J. 2004. Simultaneous analysis of multi-instrument satellite measurements of atmospheric density.JournalofSpacecraftandRockets, 41(5): 849-853.

Moe K, Moe M M. 2005. Gas-surface interactions and satellite drag coefficients.PlanetaryandSpaceScience, 53(8): 793-801.

Montenbruck O, Gill E. 2000. Satellite Orbits: Models, Methods and Applications. Heidelberg: Springer-Verlag.

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

Rim H J. 1992. TOPEX Orbit Determination Using GPS Tracking System [Ph. D. thesis]. Austin: The University of Texas.

Sentman L H. 1961. Free Molecule Flow Theory and its Application to the Determination of Aerodynamic Forces. DTIC Document.

Shi C, Li W W, Li M, et al. 2015. Calibrating the scale of the NRLMSISE00 model during solar maximum using the two line elements dataset.AdvancesinSpaceResearch, 56(1): 1-9. Solomon S C, Woods T N, Didkovsky L V, et al. 2010. Anomalously low

solar extreme-ultraviolet irradiance and thermospheric density during solar minimum.GeophysicalResearchLetters, 37(16): L16103. Storz M F Bowman B R, Branson M J I, et al. 2005. High accuracy satellite drag model (HASDM).AdvancesinSpaceResearch, 36(12): 2497-2505.

Sutton E K, Forbes J M, Nerem R S. 2005. Global thermospheric neutral density and wind response to the severe 2003 geomagnetic storms from CHAMP accelerometer data.JournalofGeophysicalResearch:SpacePhysics, 110(A9): A09S40.

Sutton E K, Forbes J M, Nerem R S, et al. 2006. Neutral density response to the solar flares of October and November, 2003.GeophysicalResearchLetters, 33(22): L22101.

Sutton E K, Nerem R S, Forbes J M. 2007. Density and winds in the thermosphere deduced from accelerometer data.JournalofSpacecraftandRockets, 44(6): 1210-1219.

Sutton E K. 2009. Normalized force coefficients for satellites with elongated shapes.JournalofSpacecraftandRockets, 46(1): 112-116.

Van Helleputte T, Visser P. 2008. GPS based orbit determination using accelerometer data.AerospaceScienceandTechnology, 12(6): 478-484. Van Helleputte T, Doornbos E, Visser P. 2009. CHAMP and GRACE accelerometer calibration by GPS-based orbit determination.AdvancesinSpaceResearch, 43(12): 1890-1896. Wang H B, Zhao C Y. 2008. Use CHAMP/STAR accelerometer data to evaluate atmospheric density models during solar maximum year.ActaAstronomicaSinica(in Chinese), 49(2): 168-178. Xiong Y Q, Wang H B, Zhao C Y, et al. 2011a. A rapid and effective calibration method for CHAMP accelerometer.ProgressinAstronomy(in Chinese), 29(2): 228-237. Xiong Y Q, Wang H B, Xu X L, et al. 2011b. The calibration research for GRACE accelerometer.ScientiaSinicaPhysica,Mechanica&Astronomica(in Chinese), 41(11): 1319-1327.

Xu T H, Yang Y X. 2004. Calibration for CHAMP accelerometry data based on known Earth gravity field model.ActaGeodaeticaetCartographicSinica(in Chinese), 33(3): 200-204.

附中文參考文獻(xiàn)

陳旭杏, 胡雄, 肖存英等. 2013. 基于衛(wèi)星數(shù)據(jù)和NRLMSISE-00模型的低軌道大氣密度預(yù)報(bào)修正方法. 地球物理學(xué)報(bào), 56(10): 3246-3254, doi: 10.6038/cjg20131003.

汪宏波, 趙長(zhǎng)印. 2008. 用CHAMP加速儀數(shù)據(jù)校驗(yàn)太陽(yáng)活動(dòng)峰年的大氣模型精度. 天文學(xué)報(bào), 49(2): 168-178.

熊永清, 汪宏波, 趙長(zhǎng)印等. 2011a. CHAMP加速儀資料的快速校標(biāo)研究. 天文學(xué)進(jìn)展, 29(2): 228-237.

熊永清, 汪宏波, 許曉麗等. 2011b. GRACE加速儀資料定標(biāo)研究. 中國(guó)科學(xué): 物理學(xué)、力學(xué)、天文學(xué), 41(11): 1319-1327.

徐天河, 楊元喜. 2004. 利用現(xiàn)有重力場(chǎng)模型求定CHAMP衛(wèi)星加速度計(jì)修正參數(shù). 測(cè)繪學(xué)報(bào), 33(3): 200-204.

(本文編輯 汪海英)

Thermosphere mass density derivation using on-board accelerometer observations from GRACE satellites

LI Wen-Wen1, LI Min1,2*, SHI Chuang1,2, ZHAO Qi-Le1,2

1ResearchCenterofGNSS,WuhanUniversity,Wuhan430079，China2CollaborativeInnovationCenterforGeospatialTechnology,Wuhan430079,China

This paper mainly discusses the thermosphere mass density derivation using the GRACE on-board accelerometer data. The precision of the accelerometer calibration as well as the satellite aerodynamics model have significant impact on the derivation. In this paper the accelerometer calibration is realized through precise orbit determination (POD) process by using the on-board GPS and accelerometer data, while the aerodynamic coefficients of the GRACE twin satellites are calculated using the Sentman equations. The GRACE accelerometer data during year 2007—2009 are processed and the densities along the satellites′ orbit during this period are obtained. To assess the precision of the accelerometer-derived densities, densities from the GRACE-A satellite during 2007—2009 are compared to the results from Doornbos (2012), as well as the empirical NRLMSISE00 and HASDM models. The results indicate that the densities from this paper are systematically about 8% larger than those from Doornbos (2012), while the STD between these two independent dataset is under 10%. Their differences are mainly attributed to the different strategies adopted for accelerometer calibration and aerodynamic coefficients calculation. The accelerometer-derived densities in this paper are systematically larger than the HASDM model with a bias about 4%, while the densities from Doornbos are 4%～7% lower. The STDs of the derived densities and the results of Doornbos with respect to HASDM model are both around 30%. The bias and STD between the accelerometer-derived densities and the empirical NRLMSISE00 model are 30%～40% and around 30%, respectively.

Thermosphere mass density; GRACE; Accelerometer; HASDM; NRLMSISE00; Aerodynamic

10.6038/cjg20160903.

國(guó)家自然科學(xué)基金項(xiàng)目(41274049，41375041，41274049，41375041，41231174)，國(guó)家高技術(shù)研究發(fā)展計(jì)劃(2014AA123101)，地球空間環(huán)境與大地測(cè)量教育部重點(diǎn)實(shí)驗(yàn)室開放研究基金(15-02-05)資助.

李文文，男，1990年生，在讀博士，主要研究方向包括星載GPS低軌衛(wèi)星精密定軌，利用低軌空間物體數(shù)據(jù)反演高層大氣密度等.

E-mail:cheeselee@whu.edu.cn

10.6038/cjg20160903

P228

2015-08-05，2016-07-19收修定稿

李文文, 李敏, 施闖等. 2016. 基于GRACE星載加速度計(jì)數(shù)據(jù)的熱層密度反演. 地球物理學(xué)報(bào)，59(9):3159-3174,

Li W W, Li M, Shi C, et al. 2016. Thermosphere mass density derivation using on-board accelerometer observations from GRACE satellites.ChineseJ.Geophys. (in Chinese),59(9):3159-3174,doi:10.6038/cjg20160903.

*通訊作者 李敏，男，1983年生，2011年畢業(yè)于武漢大學(xué)，獲博士學(xué)位，副教授，目前研究方向包括利用低軌空間物體數(shù)據(jù)反演高層大氣密度，低軌衛(wèi)星精密軌道確定，GNSS精密定位等.E-mail：limin@whu.edu.cn