基于神经网络方法的Kp预报模型

刘 杨，罗冰显，刘四清，龚建村

（1中国科学院空间科学与应用研究中心，北京100190；2中国科学院研究生院，北京100049）

1 引言

Kp指数即“行星际3h磁情指数”，通过计算磁纬在48°和63°之间的全球13个地磁台站K指数的加权平均值得到。它与3h时段内地磁扰动有近似对数的关系，是衡量近地空间全球磁扰强度的重要指标之一。Kp指数每天8个值，从0到9共分28级：0，0+，1-，1，1+，2-，2，… ，9-，9。1949 年 Bartels利用早期的地磁观测资料最先引入Kp指数概念，现在能获得1932年之后长达70年无间断Kp数据，对研究日地空间环境极具价值。与Dst和AE指数不同，Kp并没有一个明显的电流体系与之相对应。

统计结果表明，Kp指数与太阳风参数之间有很好的的相关性[1，2]，如太阳风速度、密度和行星际磁场的南向分量等。Kp指数同样与近地空间的一些现象相关，如亚暴注入的位置、磁尾场线拉伸、等离子体片地向边界和等离子体层顶的位置等[3，4]。在空间天气应用中Kp也具有重要的作用，很多磁层和电离层模型需要Kp作为输入参数，如Tsyganenko地磁场模型，大气密度模型，环电流辐射带模型和磁层顶规范预报模型等[5]，OVATION模型需要输入Kp来确定椭圆极光带的赤道向边界[6]；另外可以根据估计Kp的增强对磁层电离层的粒子和电磁环境变化发布预警，能有效保护卫星和地面电力系统免受灾害性的破坏。

作为全球地磁活动指数，Kp有自身的缺陷，通常最终公布的Kp会有数周的延迟，且Kp也不适用于研究较小时间尺度的问题。为了使Kp适用于实时应用的需求，一些机构发展了现报和短期预报模型。如 Gehred 等和 Takahashi等的 Kp 现报模型[7，8]，利用计算Kp的方法处理中磁纬地区几个地磁台站的实时数据得到Kp的估计值，虽然与最终的Kp不完全一致，但具有很好的实时性。基于人工神经网络发展的Kp短期预报模型，主要利用上游太阳风数据和现报Kp值作为输入参数，最早能提前几个小时预报Kp值，主要模型有：Costello模型、Wing等的APL模型、Boberg等的和Bala等的模型[9-12]。

CME、CIR等到达地球时会引起强烈的地磁扰动，使得空间环境异常恶劣，严重影响人类活动和各种技术系统。虽然能实时监测太阳的爆发活动，大致预测扰动传播到地球的时间，但太阳风与磁层电离层的相互作用是十分复杂的非线性动态过程，很难准确地预测地磁暴（Kp>5）的发生。太阳风作为驱动磁层电离层系统的主要能量来源，能量输入的多少决定着系统的行为。利用统计方法得到的等效刻画能量输入的耦合函数主要有Perreault-Akasofu的ε函数 ε＝vB2l2sin4（θ/2）、Stamper等的太阳风功率函数Pα、Newell等的开磁通生成速率函数dΦMP/dt=v4/32sin8/3（θ/2）和 Borovsky 的向日面重联率函数 εB等[13]。最近Lu.JY等用数值模拟的方法，首次定量描述了在不同太阳风条件下电磁能和动能的输入[14]。磁层电离层对太阳风能量输入的响应主要分为三种：直接驱动过程；加载卸载过程；两种过程同时存在[13]。统计分析表明磁层对行星际扰动响应的时间约为1h，也存在长达几小时至数天的响应时间，较长的响应滞后时间意味着输入能量在磁层有存储过程[13，15]。Johnson and Wing指出，在太阳活动高年，磁层主要受行星际条件的控制，在太阳活动低年，主要受磁层内部的动力学过程的控制，并且Kp时间序列在太阳活动低年比在高年有更显著的非线性特性[10，16]。行星际向磁层电离层能量输入的复杂性，磁层电离层对扰动响应时间的不确定性以及磁层内部的动力学过程，使得预报Kp指数的难度更大。

从预报地磁活动指数（Kp、Dst、AE）的准确度方面来衡量能量耦合函数，开磁通生成速率dΦMP/dt的效果最好[17]。Newell等又考虑了太阳风与磁层顶之间的粘滞作用项n1/2v2对能量输入的贡献，综合开磁通生成率和粘滞作用项线性预报地磁活动指数，得到了更好的结果[18]。

人工神经网络是模仿人类大脑的神经网络结构和功能建立的一种信息处理系统，能够解决高度复杂的非线性问题。目前神经网络方法广泛应用于空间环境预报，除了上文中提到的Kp短期预报模型外，也建立了对Dst指数、AE指数、同步轨道相对论电子通量和太阳质子事件的神经网络预报模型[19-23]。

Kp为3h精度指数，其衡量的是这3h时段内地磁H分量变化的最大幅值，为了更好的寻找行星际扰动源与Kp之间的非线性关系，提高Kp的预报效果，把开磁通生成速率dΦMP/dt和粘滞作用项n1/2v2加入到神经网络中来。为此，建立了三个模型以不同的方式预报未来3h内的Kp值：①模型1输入位于L1点的ACE卫星提供的当前太阳风速度和密度，行星际磁场强度，以及开磁通生成率dΦMP/dt和粘滞作用项n1/2v2；②模型2的输入是在模型1的基础上加入现报Kp；③模型3输入9h延迟dΦMP/dt的和n1/2v2，以及当前的太阳风和行星际磁场参数。

2 数据选取与处理

构建从1998年到2010年长达13年的太阳风和Kp数据集，覆盖了几乎整个第23太阳活动周以及24活动周的前几年。位于L1点的ACE卫星提供了地球上游的太阳风速度密度和行星际磁场参数，这些数据来自 NASA/CDAWEB（http：//cdaweb.gsfc.nasa.gov/），开磁通生成率dΦMP/dt和粘滞作用项n1/2v2；可利用太阳风参数计算得到，最终的Kp实测值来自德国的GeoForschungsZentrum（ftp：//ftp.gfz-potsdam.de/pub/home/obs/kp-ap/），现报 Kp 来自 NOAA/SWPC（http：//www.swpc.noaa.gov/）。

本文始终是对3h时段Kp的预报，把Kp固定在每一时段结束的时间点上。将ACE卫星SWEPAM 1min精度和MAG 4min精度的数据进行30min平均，作为模型1和模型2的输入；模型3的输入数据的精度为3h，选取的是ACE卫星1h精度的数据，然后进行3h平均。开磁通量输入率和粘滞作用项也相应地处理成30min和3h精度。Kp实测值几乎是连续的，但太阳风数据存在缺失，尤其是模型3考虑了时间延迟效应，会使得有效样本数量大大减少。

3 神经网络的建立

BP神经网络是一种使用误差反向传播学习算法的前向网络。采用了包含输入层、隐层和输出层的3层BP神经网络来预报Kp，输入层与隐层之间采用带偏差的对数S型激活函数，隐层和输出层之间采用线性激活函数。将13年的数据分为两个集合，1999-2001、2003-2005和2007-2010共10年的数据为训练验证集，主要用于神经网络的训练，寻找最优的网络权值；1998、2002和2006年3年的数据为测试集，用于测试训练好的网络，评估其预报效果。测试集包含了23活动周不同阶段的大量数据，以此来检验模型的预报能力随太阳活动周的变化。在大量实验的基础上，根据不同模型的输入样本量，确定隐层神经元节点数量为12到20之间，训练算法采用了具有较快收敛速度的Levenberg-Marquardt法。

对测试结果的分析主要用到以下的几个统计参数：线性相关系数R、均方根误差RMSE、平均相对误差ARV，具体定义为：

其中T为目标值，即Kp实测值，O为网络的输出值，即Kp预报值，N代表测试样本的总数。

4 Kp预报模型

4.1 输入太阳风参数、dΦMP/dt和n1/2v2提前1～3.5h预报Kp：模型1

模型1的输入参数为太阳风速度v，密度n，行星际磁场总强度B，By分量，Bz分量，开磁通生成率dΦMP/dt和粘滞作用项n1/2v2。神经网络训练样本的构造过程中，把实测的3h时段Kp固定在该时段结束的时刻上，每天8个值对应的时刻分别为世界时3点，6点，……，24点。考虑到ACE卫星观测到的太阳风传到地球的时间以及磁层对太阳风的响应时间，假设对实测Kp值产生影响的行星际条件的时间为该时段结束前1h至4h，若t为对应的时刻（t=3UT，6UT，……，24UT），Kpt为网络输出的目标值，所选输入参数的时间为t-4时刻至t-1时刻。把t-4到t-1时刻之间的3个小时平均分为六段，每段为30min，那么每一个Kp值都对应6个输入条件组合，Kp与输入参数的非线性关系写成函数形式为：

其中，i=t-3.5，t-3，t-2.5，t-2，t-1.5，t-1。没有把这 6个组合全部应用于神经网络的训练，而是先对其进行评估，找出可能导致地磁场产生最强扰动的一组值，作为与Kpt对应的输入条件.这利用了Newell等提出的Kp与开磁通输入率和粘滞作用项的线性拟合公式[18]：

分别计算6个输入组合对应的Kplinear值，比较这6个值，选取其中Kplinear最大的一组作为最终的输入。通过这种方法，得到10年中用于训练和验证的样本数为22163。

神经网络训练完成之后，用包含1998年、2002年和2006年的数据进行测试。测试集的输入同样为30min精度，每3h时段有6组行星际条件，每一组条件输入到神经网络之后，输出一个对Kpt的预测值，第一次预测在t-3.5时刻，之后每30min有一个输出值，随时间向前推移，最后一次预测在t-1时刻。每一次预测之后更新为网络输出的最大值，预测曲线始终单调上升，最终的预测值为6个预测值中最大的。测试结果表明，Kp观测值与预测值之间的线性相关系数为0.88，均方根误差为0.65，平均相对误差为 0.23。图 1给出了 1998、2002和2006年Kp实测值与预测值之间对比的散点图，拟合曲线为y=0.77x+0.56，图2给出的是对1998年11月和2006年12月两次磁暴事件的预测结果。

从应用角度来看，模型1的优点是在t-3.5时刻就完成了第一次预测，即时间提前量为3.5h，表一给出了提前3.5h的对不同等级磁暴的预报效果的统计，如第一行所示，1998年、2002年和2006年共有288个时段Kp值达到5，提前3.5h能够预报出这一变化的有93次，所占比例为33%；三年中共有10次Kp值达到8，其中有2次提前3.5h能够预报出。

图1 模型1预测Kp与实测Kp的对比（1998、2002、2006年）

图2 模型1对1998年11月7日至10日、2006年12月14日至17日的两次磁暴的预测结果

表1 模型1提前3.5h预报的结果

模型1的输入参数均为30min精度，这样能有效反映行星际条件的瞬时变化，保证不同情况的扰动（尤其是行星际磁场By和Bz）不会因时间的积分效应平均掉，而且大部分行星际扰动从L1点传到地球的时间都超过30min，每隔30min对Kp进行一次预报，也能够对即将到来强烈地磁扰动发出警告。

4.2 输入太阳风参数、现报Kp、dΦMP/dt和n1/2v2提前1～3.5h 预报 Kp：模型 2

现报Kp是对地磁场磁扰状态的实时估计，与最终的Kp有一定的差别，但两者的相关性非常高，基本上能实时地反映磁层的状态。在模型1的基础上把现报Kp加入到神经网络输入中，能有效降低Kp序列的非线性效应，更加准确的预报将来的Kp。NOAA/SWPC现报Kp每3h发布一次，时段与Kp相同，约有几分钟的时间延迟。模型2训练集的构造原则是在模型1选出的训练样本中始终加入最新的现报Kp：如果选出的太阳风条件对应的时刻为t-3.5和t-3，则输入的是t-6时刻的现报Kp；如果选出的太阳风条件对应的时刻为t-2.5到t-1，则输入t-3时刻的现报Kp。模型2同样选取Kp实测值作为神经网络训练和测试的目标值。对测试集的测试结果为，Kp预测值与实测值的相关系数达到了0.90，均方根误差为0.62，平均相对误差为0.20。线性相关系数超过了模型1，并且降低了预测值与实测值的误差。图3给出了模型2测试集的Kp实测值与预测值之间的对比，拟合曲线为y=0.79x+0.55，图4与图2相同，给出了模型2对两次磁暴事件的预测。

图3 模型2预测Kp与实测Kp的对比（1998、2002、2006年）

图4 模型2对1998年11月7日至10日、2006年12月14日至17日两次磁暴的预测结果

4.3 输入太阳风参数、dΦMP/dt和n1/2v2提前3h预报Kp：模型 3

模型3输入当前的太阳风速度v，密度n，行星际磁场总强度B，By分量，Bz分量，以及9h延迟的dΦMP/dt和n1/2v2提前3h预报Kp。对于延迟时间的考虑，通过试验测试了从6h到18h的不同组合，发现延迟时间超过9h，网络的性能没有明显改进。若t为当前时刻，为了预测t+3时刻的Kp，把模型的输入与输出写成函数表达式为：

网络训练完成之后，利用模型3对测试集的数据进行测试，结果表明Kp预测值与实测值的相关系数为0.85，均方根误差为0.72，平均相对误差为0.27。与前两个模型相比，预报时间提前，所付出的代价是各项统计指标均变差。图5给出了预测值与实测值的对比，图6是对1998年11月和2006年12月两次事件的预测结果。

4.4 Kp模型的预报效果

从2002年的数据中随机选取了50d做测试，这段时间包含了三次较强的磁暴过程。图7给出了测试结果，（a）、（b）和（c）表示对应的太阳风条件，（d）、（e）和（f）分别为模型1、模型2和模型3的输出与实测Kp的对比。三个模型的预测曲线与太阳风条件有很好的对应关系，当行星际源发生剧烈扰动时，预测Kp明显上升。模型1和2基本上准确的预报出了磁暴（Kp>5）的发生，模型3因为提前时间更长的原因，结果偏低。

图5 模型3预测Kp与实测Kp的对比（1998、2002、2006年）

图6 模型3对1998年11月7日至10日、2006年12月14日至17日的两次磁暴的预测结果

图7 从2002年随机选取的50d的测试结果，（a）-（c）分别为行星际磁场总强度、By及Bz分量，太阳风速度和密度（d）-（e）中黑色线为Kp实测值，红色线分别表示三个模型的预测值

表2 三个模型在不同年份中的统计特性

表2列出了三个模型对不同年份数据的测试结果，从表2中可以看出，同一个模型在太阳活动低年的表现优于高年；不同模型在2006年的各项指标几乎一致，这与太阳活动低年地磁平静，Kp指数维持在较低水平有关，2002年的测试结果则显示模型1和2明显好于模型3。

5 讨论与结论

为了满足空间环境预报的需要，发展了三个不同的Kp指数短期预报模型，利用开磁通生成率和太阳风磁层粘滞作用项与Kp之间较强的线性相关性，把这两个参数同ACE卫星直接观测到的太阳风参数一起加入到神经网络中来，训练神经网络，寻找行星际条件与地磁扰动之间的非线性关系。

模型1和模型3的输入只与太阳风参数有关，模型2需要输入反映最新磁层状态的现报Kp值。模型1和模型2通过对每一Kp时段的六组行星际条件进行评估，找到可能产生最强磁扰的一组条件，巧妙地构造了神经网络训练集。测试结果表明，模型1和模型2实际预报的平均时间提前量约为2.5h，对Kp>5事件预报的平均提前时间也是2.5h，两者有很好的一致性。

通过对随机选取的磁暴事件进行测试，可以看出模型1在磁暴的起始阶段能准确地预测Kp指数的迅速增加，预测的磁暴强度与实测值符合得很好，但由于模型1输入的太阳风参数的积分时间较短且完全假设磁暴由行星际扰动所驱动，可能会导致磁暴恢复相期间对Kp预报值偏低，持续性较差。总的来说，模型1对响应时间较短、能量释放迅速的磁暴预测较好，这也是大部分仅以较短时间延迟的行星际参数为输入预报Kp模型的特点。

模型2加入了现报Kp，提高了对磁暴恢复相以及太阳活动低年由磁层内部活动引发的磁暴的预报能力，从图3、图4和图7（e）中可以看出现报Kp对模型1的改进。模型3增加了开磁通输入率和粘滞作用项的延迟时间，增大了预报的提前时间量，但图7（f）显示了Kp预测曲线与实测结果有一定的滞后时间。

为了更好地评估站点模型，表3列出了主要的基于神经网络方法的Kp预报模型和具体参数[9-12]。表3中列出的模型根据输入量的不同主要可分为两类：第一类是完全输入上游太阳风参数，第二类同时输入地磁参量（实时Kp估计值）与太阳风参数。通过比较这两类模型的相关系数和均方根误差，可以看到加入Kp估计值能提高Kp的预报精度。根据业务需求和输入量构造方式的不同，所有模型的时间提前量从1h到4h不等，其中APL模型的提前时间是通过计算太阳风从L1点传播到磁层顶的时间得到的，本文的模型能从提前3.5h至1h连续的对某一时段内的Kp做出预测。比较本文模型1，APL 3和Bala的Model 3得到，开磁通生成率和粘粘滞作用项的加入提高了预报精度，直接输入的太阳风参数也起到了重要作用。APL 1由于输入了15min精度的Kp估计值，取得了更好的效果，本文的模型2受3h Kp估计值的限制，结果稍差，但均方根误差明显低于Bala的Model 1和Model 2。总之，受益于ACE卫星提供的几乎覆盖第23活动周的数据，以及太阳风磁层能量耦合函数的引入，本文的三个模型都取得了很好的预报效果。◇

表3 本文的模型与现有模型的比较

［1］Crooker N U，Gringauz K I.On the low correlation between longterm averages of solar wind speed and geomagnetic event activity after 1976.J Geophys.Res.1993，98.59

［2］Papitashvili V O，Papitashvili N E，King J H.Solar cycle effects in planetary geomagnetic activity：Analysis of 36-year long OMNI dataset.Geophys.Res.Lett.2000，27，2797-2800

［3］Wing S，Newell P T，Sibeck D G，et al.Large statistical study of the entry of interplanetary

magnetic field Y-component into the magnetosphere.Geophys.Res.Lett.1995，22，2086-2086

［4］Sergeev V A，Malkov L，Mirsula K.Testing the isotropic boundary algorithm method to evaluate the magnetic field configuration in the tail.J Geophy.Res.1993，98，7609-7620

［5］Detman T，Joselyn J A.Real-time Kp predictions from ACE real time solar wind.AIP Conf.Proc.1999，471，729-732

［6］Newell P T，Sotirellis T，Carbary J F，et al.OVATION：Oval Variation， Assessment， Tracking， Intensity， and Online Nowcasting.Ann.Geophys.2002，20，1039-1047

［7］Gehred P A，Cliffswallow W，Schroeder J D.A comparison of USAF Ap and Kp indices and Gottingen indices.Tech.Memo.ERLSEL，NOAA，Silver Spring，Md.1995

［8］Takahashi K，Toth B A，Olson J V.An automated procedure for near-real-time Kp estimates.J Geophy.Res.2001，106：21-017

［9］Costello K A.Moving the Rice MSFM into a real-time forecast mode using the solar wind driven forecast mode.PHD dissertation，Huston，Texas，Rice Univ，1997

［10］Wing S，Johnson J R，Jen J，et al.Kp forecast model.J Geophy.Res.2005，110，A04203，doi：10.1029/2004JA010500

［11］Boberg F，Wintoft P，Lundstedt H.Real time Kp prediction from solar wind data using neural networks.Phys.Chem.Earth.2000，25，275-280

［12］Bala R，Reiff P H，Landivar J E.Real-time prediction of magnetosphere activity using the Boyle Index.Space Weather，2009，7，S04003，doi：10.1029/2008SW000407

［13］徐文耀.太阳风-磁层-电离层耦合过程中的能量支出.空间科学学报，2011，31(1)：1～14

［14］Lu J Y，Jing H，Liu Z Q，Kabin K，and Jiang Y，Energy transfer across the magnetopause for northward and southward interplanetary magnetic fields，J Geophy.Res.，doi：10.1002/jgra.50093

［15］Bargatze L F， BarkerD N， McPherron R L， etal.Magnetospheric impulse response for many levels of geomagnetic activity.J Geophy.Res.，1985.90(A7)：6387-6394

［16］Johnson J R，Wing S.A solar cycle dependence of nonlinearity of magnetospheric activity.J Geophy.Res.2005，110，A04211，doi：10.1029/2004JA010638

［17］Newell P T，Sotirelis T，Liou K，et al.A nearly universal solar wind-magnetosphere coupling function inferred from 10 magnetospheric state variables.J Geophy.Res.2007， 112，A01206，doi：10.1029/200JA012015

［18］Newell P T，Sotirellis T，Liou K，et al.Pairs of solar windmagnetosphere coupling function：Combining a merging term with a viscous term work best.J Geophy.Res.2008，113，A04218，doi：10.1029/2007JA012825

［19］Wu J G，Reiff P H.Geomagnetic storm predictions fron solar wind data with the use of neural networks.J Geophy.Res.1997，102，255

［20］Wu J G，Lundstedt H.Geomagnetic storm predictions from solar wind data with the use of dynamic neural networks.1997，102，14255-14268

［21］Gleisner H，Lundstedt H.Response of the auroral electrojets to the solar wind modeled with neural networks.J Geophy.Res.1997，102，14269-14278

［22］龚建村，薛炳森，刘四清等.神经网络方法在太阳质子事件短期预报中的应用.空间科学学报，2003，23(6)：443～451

［23］ PerryK L， GinetG P， LingA G， etal.Comparing geosynchronous relativistic electron prediction models.Space Weather，2010，8，S12002，doi：10.1029/2010SW000581