A machine-learning-based electron density (MLED) model in the inner magnetosphere

1.   Introduction

The background electron density can affect the propagation and instability of electromagnetic waves in the magnetosphere (e.g. Chen LJ et al., 2009; Xiao FL et al., 2013; Guo MY et al., 2020; Guan CY et al., 2020; Sauer K et al., 2020). The plasmasphere is composed of low-energy particles, forming a sphere-like reservoir of very cold (~1 eV), fairly dense plasma (~50−10⁴ cm⁻³) that co-rotates with the Earth. Recently, the method for Extreme Ultraviolet (EUV) image reconstruction of the plasmasphere was improved by Huang Y et al. (2021). The low-density region outside the plasmasphere is called the trough. The boundary between the plasmasphere and trough regions is the “plasmapause”. The plasmapause varies with magnetic local time (MLT) and geomagnetic activity (Chappell, 1972; Carpenter and Anderson, 1992; Larsen et al., 2007; Fu HS et al., 2010a). Observations (e.g. Goldstein et al., 2004; Darrouzet et al., 2008; Goldstein et al., 2014) have shown that the enhanced convection electric field during geomagnetic activity leads to erosion of the plasmasphere and formation of a high-density plume in the afternoon sector. In the early stage of its formation, the plasmaspheric plume is generally broad in MLT. It co-rotates with the Earth and becomes narrower as time goes on. In addition, the existence of density troughs inside the main body of the plasmasphere has been confirmed by observation (Fu HS et al., 2010b). Whistler-mode chorus mainly appears in the plasmatrough (Meredith et al., 2012); whistler-mode hiss appears mainly in the plasmasphere (Thorne et al., 1973; Wang JZ et al., 2020). Electrostatic electron cyclotron harmonic waves tend to be generated in the lower (upper) half of harmonic bands in the lower (higher) density region (Zhou QH et al., 2017). In addition, the diffusion coefficients controlling wave-particle interactions are closely related to the background density. Acceleration (scattering loss) of particles by waves generally occurs in the lower (higher) density region (e.g. Xiao FL et al., 2009, 2010, 2015; He JB et al., 2021; Yang C et al., 2021).

Direct measurement of electron density is a challenging task. It can be derived from spacecraft potentials (Escoubet et al., 1997) or, alternativelty, from the upper hybrid resonance frequency f_uh (Kurth et al., 2015). Van Allen Probes running in the inner magnetosphere have provided measurements of electron density (inferred from f_uh) that are more reliable than previous approaches. However, these in-situ measurements detect only the local electron density where the spacecraft is located. A number of empirical models have been developed to obtain a global distribution of the electron density in the inner magnetosphere under various magnetospheric conditions. The model by Carpenter and Anderson (1992) is derived from ISEE-1 data and valid for L shells between 2.25 and 8. The model presents separate empirical functions for the equatorial densities of the plasmaspheric and the trough regions. In this model the trough density depends on L shell; the plasmaspheric density depends on L shell, the long-term solar cycle, and seasonal effects. The plasmapause location is determined by the maximum geomagnetic activity index (Kp) in the preceding 24 hours. The Global Core Plasma Model (GCPM) (Gallagher et al., 2000) is an improvement of several previous models (Carpenter and Anderson, 1992; Gallagher et al., 1995). It incorporates densities of the plasmasphere, the trough, and the polar cap, and considers the influences of the Kp index and MLT on the plasmapause location. The model by Sheeley et al. (2001) is developed from density data from the Combined Release and Radiation Effects Satellite (CRRES). The plasmasphere density is a function of L shell; the trough density is a function of L shell and MLT, both of which are independent of geomagnetic activity. The model is valid for L shells between 3 and 7. Bortnik et al. (2016) developed neural network models for predicting the inner magnetospheric state. They then have published 2D and 3D neural network models for electron density based on data from Time History of Events and Macroscale Interactions during Substorms (THEMIS) Probes (Chu XN et al., 2017a, b; Bortnik et al., 2018). Zhelavskaya et al. (2017) have presented a neural network model for electron density based on EMFISIS data from the Van Allen Probes mission. However, these neural network models have not provided explicit expressions and parameters.

We construct a physics-based nonlinear network to describe the dynamic and nonlinear relationships between electron density and geomagnetic activity levels, solar cycle effects, and seasonal effects. Using data from Van Allen Probes mission during the period from September 25, 2012 to August 30, 2019, we perform machine learning with gradient descent and backpropagation algorithms to find the optimal parameters of the density model. We shall present explicit expressions with few parameters, which can be used straightforwardly by the radiation belt community.

2.   Data and Methods

2.1   Data Set

The high-frequency receiver (HFR) of the Electric and Magnetic Field Instrument Suite and Integrated Science (EMFISIS) instrumentation suite on board the Van Allen Probes mission measures electric components of plasma waves in the frequency range of 10 to 500 kHz, and provides electron density (N) data inferred from f_uh. The density data are available from September 25, 2012 to August 30, 2019. We combined them with MagEphem data from Van Allen Probes and from OMNI (Operating Missions as Nodes on the Internet). Each data sample contains N, Kp index, SYM-H index, the 13-month-average sunspot number ($\overline{R}$), the day of the year (DoY), and the position of the probe (L Shell, MLT and magnetic latitude). The minimum time interval between each sample is 1 minute. We focus on the electron density near the magnetic equatorial plane and have chosen those samples within ±3° of magnetic

latitude to form a data set (1,017,586 observations in total). Samples in the data set are binned as a function of L in steps of 0.5 L and MLT in an interval of 1 h. We plot the corresponding global distribution of the samples in Figure 1. The number of samples increases as L increases. In the region of 1.1 < L < 2.5, there are fewer than 5000 samples in each bin. In the region of 5.0 < L < 6.0, more than 8000 samples are obtained in each bin. Samples occur most frequently from the dusk-to-midnight side. The bin with the most samples (20,000 samples) is located in the region of L = 5.5−6.0 and 15 < MLT < 17. O’Brien and Moldwin (2003) investigated the role of various geomagnetic indices in modeling the plasmapause location. They found that the recent minimum Dst index provides a better model than the maximum Kp index, which is used in GCPM. Therefore, we use the minimum SYM-H in the preceding 24 hours (SYM-H_m24) to represent the geomagnetic activity level. Figures 1b−1d show that the global distribution of samples under various geomagnetic activities is basically the same. The number of samples during SYM-H_m24 > −30 nT, −50 nT < SYM-H_m24 ≤ −30 nT and SYM-H_m24 ≤ −50 nT are 791,340, 147,517 and 78,729, respectively. Under these geomagnetic activity conditions, the orbits of Van Allen Probes have good statistical coverage for 1.1 < L < 6.0 and all magnetic local times near the equator.

Figure  1.  (a) The global distribution of electron density samples within ±3° of the equator during the period from September 25, 2012 to August 30, 2019 (b, c, d) for different levels of geomagnetic activity.

DownLoad: Full-Size Img PowerPoint

2.2   Model Description and Training Algorithm

For the plasmaspheric density, we adopt equations similar to the GCPM (Gallagher et al., 2000):

where L_pp is the location of the midpoint of the plasmapause and a₉ controls the plasmapause gradient. Both terms are associated with Q = log₁₀|SYM-H_m24 |:

where

For the trough, the variation of density with L shell and MLT is described following the form of the trough model proposed by (Sheeley et al., 2001):

By integrating plasmaspheric and trough densities, the electron density is represented by:

where u₁ and u₂ are functions that connect one regional description to another:

It should be mentioned that the electron density varies with season by a cycle of 365 days and with MLT by a cycle of 24 hours. Since both cycles are symmetrically distributed, previous studies (Gallagher et al., 2000; Sheeley et al., 2001) have used cosine functions to describe the density variation with DoY and MLT. Therefore, we choose similar expressions of density variation with DoY and MLT. Meanwhile, some parameter values in Equations (1)−(15) can reasonably describe the density variation with season and MLT. To reduce the computational time and the number of parameters that need to be fitted, the values of these parameters are adopted directly from those previous studies (Gallagher et al., 2000; Sheeley et al., 2001).

We obtain optimal parameters of the density model (k₁–k₁₉) using standard methods for optimizing classical neural networks, i.e., gradient descent and backpropagation algorithms (Rumelhart et al., 1986; Hecht-Nielsen, 1989). We use the mean square error (MSE) as the cost function:

where n is the sample size of the training data set,and ${\widehat{N}}^{\left(i\right)}$ and ${N}^{\left(i\right)}$ are the calculated and measured values of N respectively. MSE is a metric of the difference between calculated and measured values.

Figure 2 shows the constructed network of neuron-like units according to Equations (1)−(15). This is a nonlinear network because these equations are nonlinear. The orange nodes represent the input parameters of the density model. The intermediate variables and the output of the model are represented by the green nodes. We follow the blue arrows to calculate the density and cost function of the model. The derivatives of the cost function with respect to the parameter ${k}_{i}$ (i = 1−19) are calculated using the backpropagation algorithm (Rumelhart et al., 1986; Hecht-Nielsen, 1989) following the red arrows. The Levenberg-Marquardt algorithm is sensitive to the local optimum when training neural networks, i.e., the results obtained from different initial values vary widely. We adopt the mini-batch and Adam algorithms commonly used in non-convex optimization problems to find the global optimum efficiently. The metric used to evaluate the density model performance is the Adjusted R-squared coefficient (${R}_{a}^{2}$):

Figure  2.  The structure of the network to build the density model. The orange nodes represent the input parameters. The green nodes represent intermediate variables and the output of the model. The numbers of the equations corresponding to each node and the parameters to be optimized are shown in the node.

DownLoad: Full-Size Img PowerPoint

where $\overline{N}$ is the mean density of data set and p is the number of fitting parameters in the model (p = 19 for our model). The value range of ${R}_{a}^{2}$ is (−∞, 1]; the larger the ${R}_{a}^{2}$ is, the better the model fits.

In order to avoid the contingency of results caused by the random division of the data set, 5-fold cross-validation (Ojala and Garriga, 2009) is used. The dataset is divided into five blocks randomly and evenly, without overlap between any of them. We perform five training sessions, taking one block in each turn as the test set and the remaining four blocks as the training set. One thousand iterations are performed in each training session; the mini-batch of each iteration is 200 samples randomly selected from the training data set. We obtain five different models from the five training sessions and calculate ${R}_{a}^{2}$ for each model based on the test set of each training session. The larger ${R}_{a}^{2}$ is, the more accurate is the model on the test set. We choose the model with the largest ${R}_{a}^{2}$ as the optimal one.

3.   Results and Discussion

Table 1 lists the optimal values of parameters in Equations (1)−(12). Figure 3 shows the regression of the modeled (N_model) and observed (N_obs) electron densities for various data sets. The classes of each dataset and the corresponding correlation coefficients (R) are shown in the panel headers. The data points are essentially distributed close to the dashed line representing the perfect fit, suggesting that the MLED model fits the observations well. The samples near the perfect fit line in Figure 3 are much larger than those away from the perfect fit line. The correlation of observations with MLED predictions is greater than 0.87 in all three datasets. Thus, the MLED model appears to offer promising prediction performance. The correlation coefficient on the test dataset is 0.90, implying that the MLED model has good generalizability. In machine learning tasks, overfitting may occur when the model is too complex and the metrics of the training set are much better than those of the test set. The concise network used here avoids overfitting. For such concise networks, it is possible that the metrics of the test set are slightly higher than those of the training set, as shown in Figure 3.

Table  1.  The optimal values of parameters in Equations (1)−(12).

Parameter	Value		Parameter	Value
k1	−0.43		k11	2.1726
k2	4.4		k12	−1.5138
k3	0.09867		k13	−7.9149
k4	−0.2555		k14	44.2847
k5	0.0022		k15	14.4557
k6	−1.956		k16	1.2235
k7	0.015		k17	11.5
k8	−1.5226		k18	187.7
k9	0.0319		k19	88.2
k10	5.7689

 | Show Table

DownLoad: CSV

Figure  3.  Regressions between the observed electron density and the simulation results of the MLED model for various data sets. The dashed line represents a perfect model fit.

DownLoad: Full-Size Img PowerPoint

In Table 2 we present the ${R}_{a}^{2}$ coefficient and the root mean square error (RMSE) for the training and test sets under different geomagnetic activities. The ${R}_{a}^{2}$ coefficient and RMSE of the MLED model are better than those of the GCPM model in both the training and test sets. Figure 4 shows the equatorial electron densities derived from the MLED model. Since the magnitudes of k₃−k₅ are very small, the electron density is insensitive to variations in DoY and $\overline{R}$. In Figure 4 we set DoY = 30 and $\overline{R}$ = 90 in the calculation of N. The results show that N is very sensitive to the variation of L. When L varies from 1.2 to 7.0, N drops by 3 orders of magnitude. The electron density changes dramatically at the plasmapause, from greater than 1000 to less than 50 as L increases by 0.1. At fixed L, the electron density is less sensitive to the change in MLT. The electron density in the plume region is roughly one order of magnitude higher than that outside the plume region. The SYM-H index has a significant effect on the locations of the plasmapause and the plume region. During low geomagnetic activity (e.g. SYM-H_m24 = −5 nT), the plasmapause is located at L ≈ 4.8 in the sector from midnight to post-noon and there is a high-density plume on the duskside. During active geomagnetic storms, the plasmapause is compressed toward the Earth and the plume region rotates westward, with a pattern similar to that of the observational results.

Table  2.  The ${R}_{a}^{2}$ coefficients and RMSE of both models under various geomagnetic activities.

		Training			Test
		MLED	GCPM		MLED	GCPM
All samples	${R}_{a}^{2}$	0.7383	−0.2282		0.7702	0.1126
	RMSE	259.0	561.1		242.9	477.3
SYM-Hm24 > −30 nT	${R}_{a}^{2}$	0.7485	−0.1423		0.8085	0.3448
	RMSE	254.7	542.9		227.9	450.8
−50 nT < SYM-Hm24 ≤ −30 nT	${R}_{a}^{2}$	0.7075	−0.4212		0.7661	0.1904
	RMSE	266.2	586.7		242.3	421.6
SYM-Hm24 ≤ −50 nT	${R}_{a}^{2}$	0.5769	−1.5331		0.5666	−1.2918
	RMSE	301.8	738.4		291.7	670.8

 | Show Table

DownLoad: CSV

Figure  4.  Simulated equatorial electron densities under different geomagnetic activities.

DownLoad: Full-Size Img PowerPoint

We extracted samples from the data set at midnight (MLT = 0), dawn (MLT = 6), noon (MLT = 12) and dusk (MLT = 18) during different levels of geomagnetic activity. We use the measured data (L, MLT, $\overline{R}$, DoY, Kp index, and SYM-H index) as input parameters for the GCPM and MLED models, and compare the simulation results with the measured densities. As shown in Figure 5, in the region of L < 2.5 the results of the MLED model are generally consistent with the measured densities while the GCPM model yields significantly higher densities. The MLED model better presents the pattern of decreasing density with increasing L. In contrast, the density calculated by the GCPM model decreases rapidly with increasing L, much lower than the measured density for L > 5 in most cases. Furthermore, the MLED model provides a better picture, than does the GCPM model, of how the plasmapause moves earthward with increasing levels of geomagnetic activity.

Figure  5.  Comparison of the simulation results of MLED (red dots) and GCPM (blue dots) models with the measured densities (black dots) for four different magnetic local times under different geomagnetic activities.

DownLoad: Full-Size Img PowerPoint

The ${R}_{a}^{2}$ coefficient and RMSE of GCPM and MLED models for MLT = 0, 6, 12 and 18 are listed in Tables 3 and 4, respectively. Under all conditions, the ${R}_{a}^{2}$ coefficient (RMSE) of the MLED model is higher (lower) than that of the GCPM model. The ${R}_{a}^{2}$ coefficient of the GCPM model is less than 0 in most cases and its maximum value is 0.3847. The minimum value of the ${R}_{a}^{2}$ coefficient of the MLED model is 0.4437 and the maximum value is 0.8102. The ratio of the RMSE of the MLED model to that of the GCPM model is in the range of 0.26 to 0.69.

Table  3.  Comparison of the ${R}_{a}^{2}$ coefficient of MLED and GCPM models for different MLT.

		MLT = 0	MLT = 6	MLT = 12	MLT = 18
SYM-Hm24 > −30 nT	MLED	0.7267	0.7444	0.7875	0.8102
	GCPM	−0.6418	−0.4475	0.2170	0.3847
−50 nT < SYM-Hm24 ≤ −30 nT	MLED	0.6873	0.7685	0.6784	0.6994
	GCPM	−0.9493	−1.1429	0.1843	0.3661
SYM-Hm24 ≤ −50 nT	MLED	−0.7856	0.4437	0.4910	0.6311
	GCPM	−2.2801	−2.9358	−2.0681	−0.2622

 | Show Table

DownLoad: CSV

Table  4.  Comparison of the RMSE of MLED and GCPM models for different MLT.

		MLT = 0	MLT = 6	MLT = 12	MLT = 18
SYM-Hm24 > −30 nT	MLED	291.6	223.2	231.9	215.4
	GCPM	714.6	531.2	445.2	387.9
−50 nT < SYM-Hm24 ≤ −30 nT	MLED	304.5	305.6	257.8	211.8
	GCPM	760.3	929.8	410.6	307.6
SYM-Hm24 ≤ −50 nT	MLED	213.2	417.1	320.3	231.9
	GCPM	834.0	1109	786.3	428.9

 | Show Table

DownLoad: CSV

Figure 6 shows the distribution of the relative error as a function of L for MLED and GCPM models under different geomagnetic activities. The mean relative errors and standard deviations are illustrated as bar plots. The relative error is defined as the ratio of the absolute error of the modeled value to the observed value: (N_model − N_obs)/N_obs. The standard deviation is used to estimate the uncertainty of the density model predictions. The absolute values of mean relative errors of the MLED model are less than ~0.3 and significantly smaller than those of the GCPM model at different L shells. Figures 6d−6f further confirm the results of Figure 5 that the GCPM model produces overestimates at lower L shells and underestimates at higher L shells. The MLED model can be incorporated into the previously developed radiation belt models (e.g. Glauert and Horne, 2005; Xiao FL et al., 2009, 2010; Su ZP et al., 2010; Shprits et al., 2015) to better forecast the dynamic evolution of energetic electrons in the radiation belt.

Figure  6.  The distribution of the relative error as a function of L for (a−c) MLED and (d−f) GCPM models under different geomagnetic activities.

DownLoad: Full-Size Img PowerPoint

4.   Results and Discussion

During the period from September 25, 2012 to August 30, 2019, Van Allen Probes covered all magnetic local times in the low latitude region of L = 1.1 − 6.0 and provided a good opportunity to measure electron density on the global scale. We develop a machine-learning-based model of the equatorial electron density (MLED) in the inner magnetosphere based on density data from the Van Allen Probes mission. This MLED model is a physics-based nonlinear network. Gradient descent and backpropagation algorithms are used to find the optimal parameters of the nonlinear network. This approach can efficiently treat nonlinear relationships among physical quantities in space plasma environments. Compared to recently developed neural network density models (Chu XN et al., 2017a, b; Zhelavskaya et al., 2017; Bortnik et al., 2018), this MLED model provides explicit expressions with very few parameters and is therefore easily described and deployed.

The MLED model describes the variation of electron density in the inner magnetosphere as a function of L, MLT, SYM-H index, DoY, and $\overline{R}$. Different from previous neural network models, it employs fundamental physical principles to describe the variation of electron density. It models the electron densities of the plasmasphere and trough separately, and predicts the plasmapause location under different geomagnetic activity levels. Moreover, this MLED model uses cosine functions to describe the density variation with DoY and MLT based on the periodic variation of electron density with season and MLT. Under different geomagnetic conditions, the electron densities calculated by this model agree remarkably well with the Van Allen Probes measurements and present a better picture of the inward movement of the plasmapause with increasing geomagnetic activity. This MLED model should improve the accuracy of forecasting the dynamic evolution of energetic electrons in the radiation belt.

It should be mentioned that the MLED model adopts a concise network structure with the input parameters L, MLT, $\overline{R}$, DoY, and SYM-H_m24 index. We use the minimum SYM-H in the preceding 24 hours (SYM-H_m24) to represent the prevailing geomagnetic activity level. This means that the MLED model can hardly reflect the plasmaspheric evolution on a short time scale, such as the highly dynamic variation of the plasmasphere on a timescale of several hours during the main phase of geomagnetic storms, and the plasmaspheric reformation on a timescale of several minutes during substorms (Su ZP et al., 2018). To reflect variations of the plasmasphere during geomagnetic storms and substorms, and thus make the results more accurate, it may be necessary to employ a more complex neural network (e.g., Recurrent Neural Network with short-term memory capabilities) with higher temporal resolution sequences of SYM-H and AE indices as input parameters. Increasing the number of input parameters and deepening the number of layers can be expected to improve the performance of the model, but doing so may not be cost-effective; for example, increasing R from ~0.9 to ~0.95, especially when solving large-scale problems, may incur high computational costs while yielding limited performance improvement. In the future, we will investigate how additional input parameters and layers of the network would affect the performance of the model, carefully balancing the performance benefits of each modification against its computational cost.

Acknowledgments

This work is supported by the National Natural Science Foundation of China grants 42074198, 41774194, 41974212 and 42004141, Natural Science Foundation of Hunan Province 2021JJ20010, Science and Technology Innovation Program of Hunan Province 2021RC3098, and Foundation of Education Bureau of Hunan Province for Distinguished Young Scientists 20B004. All the Van Allen Probes data are publicly available at https://cdaweb.gsfc.nasa.gov/pub/data/rbsp/. The OMNI data are obtained online (https://spdf.gsfc.nasa.gov/pub/data/omni/).