使用 LSTM 神经网络预测局部环流指数*
LI Hui WANG Runze WANG Chi
(中国科学院国家空间科学中心空间天气国家重点实验室,北京 100190)
(中国科学院大学,北京 100049)
摘要
磁暴期间地磁扰动的局部时间依赖性表明了预报局部磁暴指数的必要性。我们首次利用长短期记忆(LSTM)神经网络构建了SMR-LT(SuperMAG局部环流指数)的预报模型,预报时间从1小时增加到12小时。一般来说,预测性能随着提前时间的增加而降低,SMR-06 指数的预测性能优于 SMR-00、SMR-12 和 SMR-18 指数。对于提前 12 小时的预测,相关系数分别为 0.738,0.608,0.6650.738,0.608,0.665 和 0.613 。为了避免地磁静止期海量数据的过度影响,我们只使用磁暴期间的数据来训练和测试我们的模型,预测指标的改进随提前时间的增加而增加。例如,在预测提前 12 小时的风暴时间 SMR-06 指数时,相关系数和预测效率分别从 0.674 和 0.349 提高到 0.691 和 0.455。对暴雨强度预报模式性能的评估表明,强烈暴雨的相对误差通常小于中等暴雨的相对误差。
关键词 地磁暴、部分环流指数(PRCI)、人工神经网络(ANN)
分类索引 P353
0 简介
地磁暴是地球磁层中的全球性剧烈扰动,是太阳风-磁层能量通过磁重联机制 ^([1-2]){ }^{[1-2]} 耦合增强的结果之一。一般来说,高速太阳风与持续南下的行星际磁场(IMF)会导致地球西向环流的增强,从而引发以地磁场水平分量(H 分量)凹陷为特征的地磁暴。
DstD s t 指数 ^([2-6]){ }^{[2-6]} 。Dst 指数基于赤道附近 4 个地面观测站磁强计的测量结果,代表了低纬度地磁场中几个电流系统(如环流、尾流和查普曼-费拉罗电流)的影响。典型的地磁暴由三个阶段组成。首先,在风暴开始时,Dst 会突然或逐渐上升(初始阶段)。风暴开始后,Dst 保持在高于风暴前值的水平上波动,持续时间从几分钟到几小时不等。其次,随着环流的增强,Dst 值会逐渐减小。
在剧烈波动中迅速下降到最小值(主要阶段)。最后,随着干扰的减少, DD st 缓慢回升到平静水平(恢复阶段)。
地磁暴一直是地球物理学界最具挑战性和探索热情的课题之一。这不仅是因为地磁暴对全球地磁场模式有重大影响,还因为地磁暴是太阳-地球能量耦合链中最重要的一环。此外,地磁暴对通信系统、电力系统、石油管道和太空飞行器等技术基础设施也有严重影响。因此,预测地磁暴具有重要的理论和实用价值。
自 20 世纪 70 年代以来,有关 Dst 预报的研究一直在进行。Burton等人 ^([3]){ }^{[3]} 利用太阳风速度和密度以及IMF的南北分量,开发了一种用于预报 DD st的经验算法,其中包含对太阳风动压、等离子体片的粒子注入以及指数衰减等影响的近似描述,并用常微分方程进行约束。在这一模型的基础上,又提出了一些假设来描述这种非线性行为 ^([6-8]){ }^{[6-8]} 。
除了基于经验公式构建的模型外,还有从数据科学角度对 DstD s t 预测进行的研究。其中一种广受欢迎的方法是通过人工神经网络(ANN)构建机器学习模型。Lundstedt 和 Wintoft ^([9]){ }^{[9]} 开发了一种前馈神经网络,利用太阳风的密度和速度以及 IMF 的 zz 分量来预测前方的 Dst1hD s t 1 \mathrm{~h} ,可以预测磁暴的初始阶段和主要阶段。后来,为了进一步提高 ^([10,11]){ }^{[10,11]} 预报效果,又引入了一种具有特定结构的人工神经网络--循环神经网络(RNN)。不同的优化技术也被用于 Dst 预报。Wei 等人 ^([12]){ }^{[12]} 通过改变隐层神经元的激活函数,使用了多尺度径向基网络。Lazzús 等人 ^([13]){ }^{[13]} 采用蜂群优化算法实现了提前 6 小时的 DstD s t 预测。Gruet 等人 ^([14]){ }^{[14]} 进一步取得了 DstD s t 预测 1∼6h1 \sim 6 \mathrm{~h} 的更好结果。
使用长短时记忆(LSTM)网络来提高效率。此外, kk 近邻(KNN)、非线性自回归移动平均与外生输入(NARMAX)和支持向量机(SVM)等算法也被引入到Dst预测和空间物理学的其他课题中 ^([15-17]){ }^{[15-17]} 。
虽然 Dst 指数有明确的物理解释,但其隐含的简化假设,即来自低纬度或中纬度地面站的方位角不变的 HH 扰动,往往与观测结果不一致。在风暴期间, HH 凹陷往往会出现明显的黎明-黄昏不对称现象(见 Li 等人, ^([18]){ }^{[18]} ,及其中的参考文献)。鉴于在中纬度站观测到的 H - 分量在当地时间(LT)的强烈变化,Newell 和 Gjerloev ^([19]){ }^{[19]} 提出了 SuperMAG 环流指数(SMR 指数)和 SuperMAG 部分环流指数(SMR-LT 指数)。SuperMAG是一个全球性的组织和国家机构合作项目,目前运行着300多个地基磁强计。利用 98 个中低纬度磁强计,SMR 指数在概念上与 Dst 指数相同,其当地时间版本 SMR-LT(包括 SMR-00、SMR-06、SMR12 和 SMR-18)包含以各自当地时间中心点命名的当地时间成分,分别代表午夜、黎明、正午和黄昏。SMR全球指数是SMR-LT指数的算术平均值。
在这项工作中,我们首先从几个不同的角度分析了地磁暴期间 SMR-LT 指数的不对称性。研究结果表明了预测 SMR 区域指数的重要性。然后,我们将人工神经网络方法 LSTM 应用于该主题,提前提供 SMR-LT 指数 1∼12h1 \sim 12 \mathrm{~h} 的预测。为了提高风暴时数据的预测性能,我们仅使用风暴时的数据训练了一个新模型,并分析了该模型在预测磁暴强度方面的性能。
1 部分环流指数预测的重要性
每小时部分环流指数 SMR-LT 和相应的平均指数 SMR 由以下公式得出
超级磁暴*。对 1998 年至 2019 年期间的 318 次磁暴进行了分析。
Li等 ^([18]){ }^{[18]} 提出,由于磁暴过程中部分环流的重要贡献, DstD s t 指数在磁暴主阶段和恢复初期存在明显的晨昏不对称现象。图 1 显示了 2018 年 8 月 26 日磁暴期间部分环流指数的演变。不出所料,SMR-LT 的四个指数在
图 1 2018 年 8 月 26 日磁暴期间的部分环流指数
主要阶段和早期恢复阶段。黄昏侧地磁场扰动的 HH 分量(如 SMR-18 和 SMR-12)的衰减一般大于黎明侧的 HH 分量(如 SMR-00 和 SMR-06)。SMR-00、SMR-06、SMR-12 和 SMR-18 的最小值分别为 -156.9nT,-155.3-156.9 \mathrm{nT},-155.3 nT,-204.4nT\mathrm{nT},-204.4 \mathrm{nT} 和 -209.3 nT。此外,SMR-18 指数最早达到最小值,而 SMR-00 指数则在一小时后达到最小值。在随后的恢复阶段,所有指数几乎相同。
图 2 给出了 318 次磁暴中 SMRLT 最小值的统计特征及其与 SMR 最小值的时间差。上面板显示了风暴期间 SMR-LT 指数的相对强度与风暴强度之间的关系,风暴强度分别由 SMR-LT 最小值与 SMR 最小值和 SMR 最小值的比值定义。一般来说,SMR-18 的相对强度大于 1,表示地磁凹陷较大,而 SMR-06 的相对强度小于 1,表示地磁凹陷较小。SMR-00 和 SMR-12 的结果类似,在大多数情况下相对强度大于 1。
图 2 318 场磁暴中 SMR-LT 最小值的统计特征及其与 SMR 最小值的时间差
情况。318个风暴的SMR-00、SMR-06、SMR-12和SMR-18的平均最小值分别为-111.5 nT、 -80.7nT,-109.4nT-80.7 \mathrm{nT},-109.4 \mathrm{nT} 和-133.8 nT,这与之前的结果一致。此外,随着风暴强度的增大,四个指数的相对强度都趋向于 1,这表明环流趋向于更对称,正如 Li 等 ^([18]){ }^{[18]} 所指出的那样。除风暴强度外,四个部分环流指数的风暴最小值的起始时间也各不相同。下图给出了 SMR-LT 最小值到 SMR 最小值的时滞分布。对于 SMR-00、SMR-12 和 SMR-18,时滞集中在 1 小时之内。时滞的中值相同,均为 0 。平均值分别为 1.5,1.61.5,1.6 和 -0.7 h 。对于 SMR-06,很多情况下时滞大于 0,中位值为 1 h,平均值为 2.2 h。
因此,预测 SMR-LT 指数而不是 SMR 指数应该更适合区域空间天气预警或预报。
2 方法和数据集
2.1 LSTM 网络
LSTM 是一种特殊的 RNN。传统的神经网络无法处理时间序列预测。然而,RNN 可以通过在网络中设计循环来解决这一问题,从而允许信息从网络的一个步骤持续传递到下一个步骤。RNN 可被视为同一网络的多个副本,信息可从一个副本传递到下一个副本。链式性质表明,RNN 与序列和列表密切相关。然而,当相关信息与所需信息之间的差距过大时,RNN 就无法学习连接信息。这个问题被称为 "长期依赖性",Bengio 等人 ^([20]){ }^{[20]} 对其进行了深入探讨,发现它源于 RNN 训练阶段的梯度消失问题。
LSTM 是为避免长期依赖问题而明确设计的。它是由 Hochreiter 和 Schmidhuber ^([21]){ }^{[21]} 提出的,后来被许多研究人员改进和推广,使其在大量问题上发挥了巨大作用。LSTM 也具有链状结构,但重复模块有所不同。LSTM 并非只有一个神经网络层,而是有多个重复模块、
在 LSTM 中有四个,它们以一种非常特殊的方式相互作用。
LSTM 的关键在于细胞状态,它沿着整条链直向下运行,只有一些微小的线性交互。链的一端是 C_(t-1)C_{t-1} ,代表旧的细胞状态,另一端是 C_(t)C_{t} ,代表新的细胞状态。LSTM 的设计目的是在细胞状态中移除或添加信息,并通过称为门的结构进行精确调节。门是选择性地让信息通过的一种方式。它们由一个 sigmoid 神经网层和一个点乘法运算组成。西格码神经网络层输出 0 和 1 之间的数字,描述了每种成分应被允许通过的程度。0 表示完全拒绝,1 表示完全接受。
LSTM 由三个门组成,用于保护和控制细胞状态,即遗忘门层、输入门层和输出门层。遗忘门层是一个 sigmoid 层,用于决定从单元状态中遗忘哪些信息。它接收 h_(t-1)h_{t-1} 和 x_(t)x_{t} ,分别表示前体副本的输出和时间步骤 tt 的输入,并输出一个介于 0 和 1 之间的数字,表示细胞状态 C_(t-1)C_{t-1} 被遗忘的程度。在随后的等式中,将保留 w_(**)w_{*} 和 b_(**)b_{*} 分别代表层的相应权重和偏置的符号。
f_(t)=sigmoid(w_(f)*[h_(t-1),x_(t)]+b_(f)).f_{t}=\operatorname{sigmoid}\left(\boldsymbol{w}_{\boldsymbol{f}} \cdot\left[h_{t-1}, x_{t}\right]+b_{f}\right) .
下一步是决定在细胞状态中存储哪些新信息,这包括两个部分。首先,一个称为 "输入门层 "的 sigmoid 层决定哪些值需要更新。接下来,一个 tanh 层会创建一个可添加到状态中的新候选值向量 tilde(C)_(t)\tilde{C}_{t} 。
{:[i_(t)=sigmoid(w_(i)*[h_(t-1),x_(t)]+b_(i))","],[ tilde(C)_(t)=tanh(w_(C)*[h_(t-1),x_(t)]+b_(C))]:}\begin{aligned}
& i_{t}=\operatorname{sigmoid}\left(\boldsymbol{w}_{i} \cdot\left[h_{t-1}, x_{t}\right]+b_{i}\right), \\
& \tilde{C}_{t}=\tanh \left(\boldsymbol{w}_{C} \cdot\left[h_{t-1}, x_{t}\right]+b_{C}\right)
\end{aligned}
在下一步中,将这两者结合起来,将旧的单元状态 C_(t-1)C_{t-1} 更新为新的单元状态 C_(t)C_{t} 。我们将旧状态乘以 f_(t)f_{t} ,再加上 i_(t) tilde(C)_(t)i_{t} \tilde{C}_{t} 。这就是新的候选值,其比例是我们决定更新每个状态值的大小。
C_(t)=f_(t)C_(t-1)+i_(t) tilde(C)_(t)C_{t}=f_{t} C_{t-1}+i_{t} \tilde{C}_{t}
最后,输出门层将生成细胞状态 C_(t)C_{t} 的过滤版本。首先,将运行一个 sigmoid 层,该层决定输出细胞状态的哪些部分--C_(t)C_{t}。
放。然后,我们将单元状态通过 tanh(将数值推至-1 和 1 之间)并乘以 sigmoid 门的输出,这样就只输出了我们决定输出的部分。
{:[o_(t)=sigmoid(w_(o)*[h_(t-1),x_(t)]+b_(o))","],[h_(t)=o_(t)tanh C_(t)]:}\begin{aligned}
& o_{t}=\operatorname{sigmoid}\left(\boldsymbol{w}_{o} \cdot\left[h_{t-1}, x_{t}\right]+b_{o}\right), \\
& h_{t}=o_{t} \tanh C_{t}
\end{aligned}
在本研究中,使用了几个性能指标来评估机器学习预测模型的有效性。均方根误差 (E_(rms))\left(E_{\mathrm{rms}}\right) 用于表示观察值和预测值之间的误差。
E_(rms)=sqrt((1)/(N)sum_(i=1)^(N)(y_(i)- widehat(y)_(i))^(2))E_{\mathrm{rms}}=\sqrt{\frac{1}{N} \sum_{i=1}^{N}\left(y_{i}-\widehat{y}_{i}\right)^{2}}
皮尔逊相关系数 CC 用于表示预测值和观测值之间的线性关系。 1//-11 /-1 表示完全正相关/负相关,而 0 表示没有线性相关。
{:[C=(Cov(y,( widehat(y))))/(sqrt(Var(y)Var( widehat(y))))=],[(sum_(i=1)^(N)(y_(i)- bar(y)_(i))( widehat(y)_(i)- bar(widehat(y))_(i)))/(sqrt(sum_(i=1)^(N)(y_(i)- bar(y)_(i))^(2)sum_(i=1)^(N)( widehat(y)_(i)- bar(widehat(y))_(i))^(2)))]:}\begin{aligned}
C= & \frac{\operatorname{Cov}(y, \widehat{y})}{\sqrt{\operatorname{Var}(y) \operatorname{Var}(\widehat{y})}}= \\
& \frac{\sum_{i=1}^{N}\left(y_{i}-\bar{y}_{i}\right)\left(\widehat{y}_{i}-\overline{\widehat{y}}_{i}\right)}{\sqrt{\sum_{i=1}^{N}\left(y_{i}-\bar{y}_{i}\right)^{2} \sum_{i=1}^{N}\left(\widehat{y}_{i}-\overline{\widehat{y}}_{i}\right)^{2}}}
\end{aligned}
此外,还使用了预测效率 P_(e)P_{\mathrm{e}} ,其计算方法如下。
P_(e)=1-(sum_(i=1)^(N)(y_(i)- widehat(y)_(i))^(2))/(sum_(i=1)^(N)(y_(i)- bar(y)_(i))^(2)).P_{\mathrm{e}}=1-\frac{\sum_{i=1}^{N}\left(y_{i}-\widehat{y}_{i}\right)^{2}}{\sum_{i=1}^{N}\left(y_{i}-\bar{y}_{i}\right)^{2}} .
P_(e)=1P_{\mathrm{e}}=1 表示完美预测, P_(e)=0P_{\mathrm{e}}=0 表示模型的性能等同于观测数据的算术平均值。 P_(e)P_{\mathrm{e}} 可以是负值,在这种情况下表示模型的预测效果并不比取测试数据的算术平均值好。 CC 只反映时间序列趋势的一致性,与之相反, P_(e)P_{\mathrm{e}} 还可以表示模型特征的振幅,从而检验预测的准确性。
2.3 数据预处理
太阳风参数由 OMNI da-
ta*由美国国家航空航天局(NASA)国家空间科学数据中心(NSSDC)维护。SMR 和 SMR-LT 指数来自 SuperMAG 合作者**。
合理选择特征参数可以大大提高方差网络的学习效率。本研究采用太阳风场 E_(y)E_{y} 、太阳风动压 P_(d)P_{d} 、太阳风速度 V_(p)V_{p} 、太阳风密度 N_(p)N_{p} 、太阳风温度 T_(p)T_{p} 、IMF强度 B_(t)B_{t} 及其 zz 分量 B_(z)B_{z} ,用 D_("sw ")D_{\text {sw }} 表示。
D_(sw)=([E_(y),P_(d),V_(p),N_(p),T_(p),B_(t),B_(z)]).D_{\mathrm{sw}}=\left(\begin{array}{lllllll}
E_{y} & P_{d} & V_{p} & N_{p} & T_{p} & B_{t} & B_{z}
\end{array}\right) .
我们从 OMNI 和 SuperMAG 中选择的数据集时间跨度为 1998 年 1 月 1 日至 2019 年 12 月 31 日,时间分辨率为 1 小时。22 年的数据集分为三个部分,即训练集(1998-2009 年)、验证集(2010-2014 年)和测试集(2015-2019 年)。训练集和验证集都考虑了一个以上的完整太阳活动周期。在这两个周期内,共识别出 257 个磁暴。
我们的目标是提前几个小时预测 SMR-LT 指数(记为 pp ),这在数据结构上意味着要开发一个用于多滞后时间步的多变量时间序列预测的 LSTM 模型。在本研究中,考虑的是 p=1,3,6,12hp=1,3,6,12 \mathrm{~h} 。第一步是为 LSTM 准备数据集,这包括将数据集设定为一个监督学习问题,并对输入变量进行归一化处理。我们将监督学习问题定义为:根据之前多个时间步骤的太阳风参数,预测当前时刻的 SMR-LT 指数。OMNI 数据已经从其观测位置时移到地球弓形震荡位置,这是为了更好地支持太阳风-磁层耦合研究。接下来,我们分别对每个特征进行归一化处理,将物理量转化为可由 ANN 处理的数据流。ANN 预测结果后,我们再应用缩放过程来获得真正的 SMR-LT 指数。
RNN 的重要设计模式有很多不同的例子。我们采用的是隐单元之间具有递归连接的神经网络,该网络读取一个
然后产生单一输出,即 SMR 指数或 SMR-LT 指数之一。我们构建了两个隐藏层,第一层包含 50 个神经元,第二层包含 100 个神经元。根据 LSTM 处理数据的结构特征,每个参数都是时间序列上的多维向量。时间序列的长度(用 ss 表示)会显著影响模型的性能。为了评估这种影响,图 3 给出了在不同 ss 条件下得到的结果。一般来说,当 s=3s=3 时,模型性能最佳。因此,在下面的研究中,我们将时间序列的长度设为 3,因此输入为 3xx83 \times 8 矩阵。
输入 =([E_(y),P_(d),V_(p),N_(p),T_(p),B_(t),B_(z)])o+=\left(\begin{array}{lllllll}E_{y} & P_{d} & V_{p} & N_{p} & T_{p} & B_{t} & B_{z}\end{array}\right) \oplus 索引、 N_("LSTM "):Input(t-p)|->Index(t),p in{1,3,6,12}\mathcal{N}_{\text {LSTM }}: \operatorname{Input}(t-p) \mapsto \operatorname{Index}(t), p \in\{1,3,6,12\} 。
3 项成果
鉴于太阳风的 7 个输入特征,共有 127
图 3 使用不同序列长度 (s) 的模型预测 SMR 指数的性能。四个板块表示考虑不同 pp 时的结果
Lazzús等人的 ^([23]){ }^{[23]} 在提前6小时和12小时的预测结果较好,但在提前1小时和3小时的预测结果较差。Lazzús 等人的 ^([23]){ }^{[23]} 得到了提前 1 小时到 6 小时的预测结果, CC 分别为 0.978,0.8950.978,0.895 和 0.788 。总体而言,我们的结果与已发表的模型大致相当,这验证了我们模型的可靠性。
图 4 给出了模型对测试集 SMR-LT 指数预测的散点图。蓝线代表拟合结果,黑色虚线代表精确预测结果。蓝色
表 1 参考文献[14]、[15]和[16]之间的比较[14]、参考文献[22]、参考文献[23] 和我们提出的模型
p//hp / \mathrm{h} |
坚持不懈 |
我们的模式 |
参考文献 [14] |
参考文献 [22] |
参考文献 [23] |
1 |
0.945 |
0.965 |
0.966 |
0.845 |
0.978 |
3 |
0.853 |
0.903 |
0.923 |
0.872 |
0.895 |
6 |
0.755 |
0.824 |
0.865 |
0.864 |
0.788 |
12 |
0.592 |
0.705 |
- |
0.857 |
- |
p//h Persistence Our model Ref. [14] Ref. [22] Ref. [23]
1 0.945 0.965 0.966 0.845 0.978
3 0.853 0.903 0.923 0.872 0.895
6 0.755 0.824 0.865 0.864 0.788
12 0.592 0.705 - 0.857 -| $p / \mathrm{h}$ | Persistence | Our model | Ref. [14] | Ref. [22] | Ref. [23] |
| :---: | :---: | :---: | :---: | :---: | :---: |
| 1 | 0.945 | 0.965 | 0.966 | 0.845 | 0.978 |
| 3 | 0.853 | 0.903 | 0.923 | 0.872 | 0.895 |
| 6 | 0.755 | 0.824 | 0.865 | 0.864 | 0.788 |
| 12 | 0.592 | 0.705 | - | 0.857 | - |
当 p=1p=1 与黑色虚线的距离越近时,预测结果越好。很明显,当 p=1p=1 时,预测结果与观测结果非常接近。随着 pp 的增加,预测结果(主要是低估)与观测结果的偏差也相应增加。
表 2 给出了我们的预测模型在 SMR-00、SMR-06、SMR-12 和 SMR-18 指数上的定量表现。考虑到之前没有关于 SMR-LT 指数预测的研究,我们将模型的预测结果与持久性操作的结果进行比较,以定量表示模型的能力。总体而言,从 E_(rms),CE_{\mathrm{rms}}, C 和 P_(e)P_{\mathrm{e}} 的角度来看,我们的模型可以很好地预测SMR-LT指数。在用 p=1,3,6,12hp=1,3,6,12 \mathrm{~h} 对SMR-06进行预测时, E_(rms)E_{\mathrm{rms}} 分别为5.752、 6.592,7.6396.592,7.639 和9.976 nT,分别比持久运算结果 3.4%3.4 \% 、 17.7%,22.1%17.7 \%, 22.1 \% 和 18.0%18.0 \% 小; CC 分别为 0.925,0.873,0.8250.925,0.873,0.825 和0.738 ,分别 2.5%,6.2%,12.4%2.5 \%, 6.2 \%, 12.4 \% 、 25.1%25.1 \% 优于持久化操作的结果; P_(e)P_{\mathrm{e}} 为
图 4 模型对测试集上 SMR-LT 指数预测的散点图。蓝线表示拟合结果,黑色虚线表示精确预测结果
0.817,0.759,0.677,0.4490.817,0.759,0.677,0.449 ,分别为 1.6%1.6 \% 、 17.7%,44.7%17.7 \%, 44.7 \% 和 149.4%149.4 \% ,优于持续运算的结果。对于SMR-18的预测, p=1,3,6,12hp=1,3,6,12 \mathrm{~h} 的 E_(rms)E_{\mathrm{rms}}