热层大气密度变化特征与太阳辐射和地磁指数的相关性分析

苗 娟，刘四清，李志涛，黄文耿，唐歌实

（1中国科学院空间科学与应用研究中心，北京100190；2航天飞行动力学技术重点实验室，北京100094）

1 引言

大气密度与低轨道航天器所受的大气阻力紧密相关，是影响卫星轨道的重要因素。要准确预测航天器的轨道，就要对大气密度做出精确的预报。从20世纪50年代开始，发展了多种半经验的高层大气模型，形成了 CIRA 系列[1]、Jacchia 系列[2]、DTM 系列[3，4]、MSIS 系列[5，6，7]等各种参考大气模型，为空间活动提供服务。但由于高层大气的变化极其复杂，几十年来虽然大气模型在不断地改进和发展，但一般情况下模型仍存在15%～30%左右的误差，在空间环境扰动期间可达100%甚至更高[8，9]。在航天器测定轨中，大气密度模型误差将转化为轨道误差，影响操控计划、精密定轨、碰撞规避及航天器再入等。因此，如何提高大气模型的预报精度是卫星和航天服务一直致力解决的重要问题。

2000年Tobiska[10]提出在模型中直接用E10.7来替代传统的F10.7能提高SEM卫星轨道衰减的预报精度、2006年Bowman[11]等人建立了两种新的太阳辐射指数S10.7和M10.7，建立了改进的大气密度模型JB2006。2008年Bowman再次增加了一个新的太阳辐射指数Y10.7，并且引入了地磁指数Dst，建立了新模型JB2008[12]。

本文利用CHAMP卫星加速度仪资料反演的大气密度数据，分析大气密度变化与五种辐射指数和两种地磁指数的变化关系，及发生磁暴时大气密度的变化特征，并利用Dst指数对密度变化整体态势进行预报和评估，为大气密度模型的改进奠定基础。

2 大气密度与太阳辐射和地磁指数变化分析

CHAMP卫星是德国地球科学研究中心组织研制的一颗科学小卫星，于2000年7月15日发射，初始高度为454km、运行于倾角为87.3°近圆极轨道上，绕地周期约为90min。本文所采用的数据为美国科罗拉多大学CHAMP研究小组发布的2001年－2008年400km标准高度上的大气密度，此段时间包括了第23太阳活动周峰年及整个下降阶段，为开展大气密度与太阳辐射、地磁暴的变化特征分析提供了充足的数据。

由于纬度、地方时的不同都能引起大气密度的变化，为了突出地磁活动对密度的影响，消除地方时和纬度对密度的影响，文中对卫星运行一个轨道周期内的大气密度进行了平均，以便更加清楚地分析大气密度整体变化特征，并采用一个地方时周期为一个密度平均周期。

分别选择2005年8月27日（F10.7为94，地磁ap指数均小于10）和2001年11月5－7日（F10.7最大为230，ap最大为300）作为平静期和磁暴期的代表，图1、图2分别给出地磁平静和磁暴两种状态下经过平均之后的大气密度变化情况。由图1可以看到，在同一天内，卫星运行各个轨道上的平均密度值变化不大，比较稳定，密度值维持在 1.3×10－12kg/m3，可以推算出密度对卫星轨道衰减的平均影响也比较稳定。图2显示，平均密度能够清晰地反映出密度对磁暴响应的整个过程，并能给出密度的最大变化值和时间，可为轨道衰减幅度分析提供参考。因此可以说，一个轨道周期的平均密度，可以清晰反映密度的整体变化特征和趋势，同时对轨道的平均变化也具有一定意义。

图1 地磁平静时平均密度变化

图2 磁暴时平均密度变化

2.1 大气密度随太阳辐射指数的变化

以F10.7作为太阳辐射能量参数的大气模型已经过了几十年的发展，模型的误差始终未能得到很好的改善，研究者开始寻求能够表征太阳辐射能量的新参量，这比改变百万行的原始大气模型代码容易得多[13]。对此，曾有研究者提出 E10.7，S10.7，M10.7，和 Y10.7等新的太阳辐射指数来改善大气模型[10，11，12]。

新的辐射指数是否比传统的F10.7能更好地表征太阳辐射对大气密度的影响，本节对五种辐射指数和密度的变化关系进行统计分析。图3给出了F10.7，E10.7，M10.7，S10.7，Y10.7五种太阳辐射指数的月均值比较，为了图形能更清楚地展示各曲线的关系，本图只给出了部分数据（2001－2004年）。明显看出，五种指数的变化规律基本一致，E10.7整体略高于其它几个指数，特别是在太阳活动高年期间更加明显。图4～图8分别给出了大气密度与五种辐射指数分布关系，可以看到密度的变化与五种指数变化趋势一致，相关系统Ro0.97以上。同时也可以看到，在150sfu以下的相关性符合更好，这可能是由于在太阳活动低年，太阳上的突然爆发事件较少，大气密度和太阳辐射相对都比较稳定。

从几种太阳辐射指数与密度都具有很好相关性也可以得到，若不考虑每一种辐射指数的物理含义，单从数据统计意义上来表征辐射对密度的影响的时，任何一种辐射指数均能很好地反映密度的变化，这也可能是为什么长期以来大量模式采用易于获取的F10.7来表征太阳辐射大小的原因之一。

图3 五种辐射指数的比较

图4 大气密度与F10.7

图5 大气密度与E10.7

图6 大气密度与M10.7

图7 大气密度与S10.7

图8 大气密度和Y10.7

2.2 大气密度随地磁指数的变化

除太阳辐射指数外，一直以来地磁指数被作为大气模式的输入来表征地磁活动状态对高层大气的影响。在JB2008模式之前，各类参考大气模式中主要采用Ap/ap/Kp指数作为模式输入参数。2008年，Bowman在新建的JB2008模式中用Dst代替了ap指数[12]。

大气密度与这两种地磁指数的变化关系如何？本文对2001年－2008年400km高度的CHAMP卫星大气密度与ap、Dst两者的变化关系分别进行统计，图9为密度和ap的变化关系，图10为密度和Dst的变化关系。从图可以看到，ap和Dst与密度的变化关系都分为两个阶段。对于ap指数，当ap≥50nT时，ap与密度呈线性相关，相关系数为0.84，而当ap<50nT时，两者的关系更符合指数函数，其相关性为0.99；而对于Dst指数，当Dst≤0时，密度与其有很好的负相关性，线性相关系数达-0.95，而对于Dst＞0，两者关系很难用一个线性或者指数函数表示。通过对Dst大于0的数据进行分析，除了在地磁平静期Dst＞0之外，Dst的突然增加往往都是出现在磁暴前的急始阶段。

图9 大气密度与ap指数变化

图10 大气密度和Dst指数变化

从与密度的线性关系来看，Dst与密度变化更一致。这可能是由于ap指数是由Kp指数等效得到的，属于分级指数，其值并不连续，在最小值0与最大值400之间分为28个等级，导致最多只能有28个格点对其进行统计，而Dst指数属于无级指数，没有人为设定的上下限和间断，这也可能是导致ap指数与密度的线性相关性不如Dst的原因。

2.3 磁暴期间大气密度变化特征

为了进一步分析磁暴对高层大气的影响，本文对2001年－2008年Dst≤-50的60个磁暴事件中的密度及地磁指数ap、Dst进行了统计分析。主要分析磁暴的发生时间、大小和密度变化之间的关系。

事实上，一次典型磁暴的过程往往会伴随太阳辐射的变化，而大气密度的变化中必然包含了这部分影响。要分析不同地磁暴对密度的程度影响，就应剔除密度中太阳辐射变化带来的影响，在同一辐射标准下比较才有意义。这就需要对60个磁暴的进行统一辐射水平“标准化”，即在其它条件不变的情况下把密度都标准化到同一个F10.7情况下，具体做法如下：

其中，ρobs（F）为实测密度，ρmodel（F）为实测太阳辐射指数下的模式计算密度，本文中采用的是NRLMISISE-00国际参考大气模式，ρ′model（F0）为太阳辐射指数标准值情况下的模式计算密度，ρ′obs（F0）为标准后的大气密度。

本文中把60个磁暴都标准化到同一太阳辐射水平150sfu后，分析密度、Dst、ap指数之间的变化关系。图11（a）—（d）给出了4个较大磁暴事件的例子，图中的蓝色竖直虚线分别表示密度开始增长和恢复到磁暴前水平。通过比较密度、Dst和ap之间变化，可以看到：1）密度的变化和两种地磁指数的变化时间都比较一致，但相比于ap，密度的开始上升时间与Dst的开始时间更加吻合；2）密度的恢复时间要比ap的恢复时间慢；但要比Dst恢复的快，对60个磁暴事件统计表明，大部分磁暴在Dst恢复至60%～80%时，密度就已恢复到磁暴前水平。3）图11（b）给出的是 ap 指数同为 300，Dst分别为－373、－289的连续两个磁暴过程，可以看到，两个磁暴中的密度变化有较大差别，从变化幅度来看，密度的变化与Dst更加一致。

为了更进一步分析磁暴时密度的变化，对辐射指数标准化之后的60个磁暴事件进行统计，具体规定与做法如下：1）磁暴之前的水平：密度出现连续上升之前的几点平均值密度值作为磁暴前的密度水平，对应的几个点的平均Dst/ap为磁暴前的地磁水平，分别表示为 ρbegin，Dstbegin，apbegin；2）对一个磁暴事件，密度的最大值作为这个磁暴时的最大密度，Dst/ap的极值表征此次磁暴强度，分别表示为ρmax，Dstmax，apmax；3）当密度再次回到磁暴前的水平，表示磁暴对大气密度的影响结束。对60个磁暴事例分别进行了 ρbegin，Dstbegin，apbegin，ρmax，Dstmax，apmax统计，图12、13给出分别用Dst和ap表征最大磁暴与密度最大值之间的关系。

图11 （a）－（d）大气密度与 Dst、ap 的变化比较

图12 Dst与密度变化关系

图13 ap与密度变化关系

比较可以看到，无论用Dst还是用ap来表征磁暴程度，从总体趋势来看，磁暴的大小与密度的变化是一致的，磁暴越大，密度增加幅度越大，这与前面总的统计结果一致；同时从个例来看，并非一个大磁暴对应的大气密度一定大于一个小磁暴的密度，这个结论与陈光明的研究结果一致。分析其原因可能是由于磁暴对大气密度的影响是一个复杂的过程，用一个地磁指数可能无法全面地描述它的变化；如果只用地磁指数表征磁暴期大气密度的变化，两个指数相比，Dst与密度之间的线性相关系数为-0.86，ap与密度之间的线性相关系为0.71，Dst与密度的相关性好于ap。

2.4 大气密度和Dst的一元线性拟合

上面分析可以看到，无论是一般情况下，还是在磁暴事件中，密度的整体变化与ap、Dst指数都有很好的相关性，并且密度与Dst变化更一致。基于Dst和密度变化关系，可以建立一个密度、Dst之间的一元线性回归方程，验证是否可以通过Dst指数的变化就可以较好地预测、特别是磁暴时大气密度的整体态势变化，一元回归方程如下：

其中，ρ密度，Dst为地磁指数，P0，P1为待定系数。

预留出2003年的数据作为验证之外，利用2001－2008年的其它所有CHAMP大气密度数据和Dst进行P0、P1系数拟合，再用所拟合得到的系数及实际的Dst指数对2003年的密度进行计算，并与CHAMP实测结果及NRLMSISE-00模式计算结果进行比对。图14是为2003年11月18－23日的一次磁暴过程，比较看出，模式对磁暴的反映很不足，而直接采用Dst值拟合的密度能很好地反映密度的整体态势变化。为了进一步验证拟合结果的适应性，用CHAMP拟合的系数对GRACE-A进行验证，图15为GRACE-A卫星2003年8月16－22日的一次磁暴过程，同样可以看到，拟合结果比模式计算更接近于实测结果。

图14 CHAMP验证

图15 GRACE-A验证

3 结论与讨论

本文利用CHAMP卫星加速度仪资料反演的大气密度数据，分析了 F10.7，E10.7，Mg10，S10，Y10 五种太阳辐射指数和两种地磁指数Dst、ap与密度之间的变化关系，得到如下结论：

1）密度的变化与五种指数变化一致，相关程度都在0.97以上。若不考虑辐射指数的物理含义，如从数据统计意义上来表征辐射对密度的影响，任何一种辐射指数均能很好地反映密度的变化。

2）密度与地磁指数ap、Dst的变化都具有很好的相关性，通过对60个磁暴事例的统计，可看到在磁暴期，无论从上升时间及增加幅度上，密度的变化与Dst更加一致，这可以为今后模式的建立和改进提供参考。

3）在密度与Dst变化关系的基础上，建立一种简单的利用Dst指数预测大气密度整体态势变化的拟合关系，通过对CHAMP和GRACE-A两个卫星的验证，可以看到大气密度整体态势的变化能够通过Dst指数的变化很好地反映出来，并且比利用国际参考大气NRLMSISE-00计算的结果更接近观测值，这可以为航天器测定轨提供参考。

本文研究的数据是以CHAMP为基础，通过GRACE-A数据进行了适应性验证，但因这两颗卫星同属于极轨卫星，对于其它轨道倾角卫星的适用性还有待于进一步的分析。同时，本文中只构建了密度和Dst的一元线性回归关系，暂时未考虑ap指数，而鉴于ap指数的通用性和对Dst的互补性特征，后期将进一步开展大气密度与两种指数的加权分析。

［1］H.Kallmann-Bijl，R.L.F.Boyd，H.Lagow，S.M.Poloskov，and W.Priester （eds.）.Cira 1961:Cospar International Reference Atmosphere 1961，North-Holland Publishing Company，Amsterdam，1961.

［2］Jacchia，L.G..New static Models of the thermosphere and exosphere with empirical temperature models.Technical Report 313，Smithsonian Astrophysical Observatory，1970.

［3］Berger，C.，Biancale，R.，Ill，M.，Barlier，F.Imporvement of th empirical thermosphere model DTM:DTM-94—a comparative review of various temporal variations and prospects in space geodesy applications.Journal of Geodesy 72（3），161-178，1998.

［4］Bruinsma，S.，Thuillier，G.，Barlier，F..The DTM-2000 empirical thermosphere model with new data assiminlation and constraints at lower boundary:accuracy and properties，Journal of Atmosphere and Solar-terriestrial Physics65，1053-1070，2003.

［5］Hedin，A..MSIS-86 thermosphereic model.Journal of Geophysical Research 92（A5），4649-4662，1987.

［6］Hedin，A.E.Extention of the MSIS thermospheric model into the middle and lower atmosphere.Journal of Geophysical Research 96（A2），1159-1172，1991.

［7］Picone，J.，Hedin，A.，Drob，D.，Aikin，A..NRLMSISE-00 empirical model of the atmophere:statistical comparisons and scientific issues.Journal of Geophysical Research 107（A12），2002.

［8］Marcos，F.，Bass，J.，Baker，C.，Boner，W..Neutral density models for aerospace applications.In:32nd Aerospace Sciences Meeting and Exhibit，January 10－13，1994/Reno，NY，1994

［9］Rhoden，E.，Forbes，J.，Marcos，F..The influence of geomagnetic and solar variabilities on lower thermosphere density.J.Atmos.Sol.-Terr.Phys.62，999－1013，2000

［10］Tobiska W K..Validating the solar EUV proxy E10.7.J Geophys Res，2001，106（A12）:29969—29978

［11］Bowman，B.R..The JB2006 empirical thermospheric density model，Journal of Atmospheric and Solar-Terrestrial Physics 70（2008）774-793.

［12］Bowman，B.R..A New Empirical Thermospheric Density Model JB2008 Using New Solar and Geomagnetic Indices，AIAA/AAS Astrodynamics Specialist Conference 18-21 August 2008，Honolulu，Hawaii

［13］Tobiska，W.K..Forecast E10.7 for Improved LEO Satellite Operations，American Institute of Aeronautics and Astronautics，2002