An Ionospheric Total Electron Content Model with a Storm Option over Japan Based on a Multi-Layer Perceptron Neural Network

Abstract

Ionospheric delay has a severe effect on reducing the accuracy of positioning and navigation of single-frequency receivers. Therefore, it is necessary to construct a precise regional ionospheric model for real-time Global Navigation Satellite System (GNSS) applications. The total electron contents (TECs) of 839 GNSS stations affiliated with the GPS Earth Observation Network were used to build a Japanese ionospheric model (JIM) based on a multi-layer perceptron neural network. During quiet space conditions, the correlation coefficient between the targets and the predictions of the JIM was about 0.98, and the root-mean square error (RMSE) of TEC residuals was ~1.5TECU, while under severe space events, the correlation coefficient increased to 0.99, and the corresponding RMSE dropped to 0.96 TECU. Moreover, the JIM model successfully reconstructed the two-dimensional (time vs latitude) TEC maps, and the TEC maps had evident hourly and seasonal variations. Most of TEC residuals accumulated between universal time 01–06 with an averaged magnitude of 1-2TECU. Furthermore, the JIM model had a perfect prediction performance under various kinds of complex space environments. In the quiet days, the prediction accuracy of the JIM was nearly equal to the global ionosphere map (GIM), and in some moments, the JIM was more competitive than the GIM. In the disturbed days, the RMSEs of TEC residuals were proportional to the solar wind speed and were inversely proportional to the geomagnetic Dst value. The maximum RMSE of the JIM was lower than 2TECU, while the corresponding RMSEs for the IRI and TIE-GCM exceeded 5TECU.

1. Introduction

Ionospheric delay is a major contributor to the reduction of positioning and navigation accuracy, significantly impacting the quality of Global Navigation Satellite System (GNSS) applications. The total electron content (TEC) serves as a primary parameter for representing the intensity of ionospheric activity, and accurate TEC prediction plays a vital role in understanding Earth’s space environment and other associated scientific applications. Currently, TEC is typically estimated using the differential technique of multiple frequency observations with high precision [1]. For the real-time TEC prediction, the NeQuick model [2,3], Klobuchar model [4], Bent model [5], International Reference Ionosphere (IRI) model [6,7], BeiDou global broadcast ionospheric delay correction model (BDGIM) [8] and other physical models [9] are usually recommended. These models are generally successful in reflecting the global ionospheric variability and have good applicability in large-scale regions. However, these models have some shortcomings in small-scale areas, such as low accuracy and low spatial resolution, etc. Therefore, they cannot fully meet the demands of regional GNSS users.

Due to the actual requirement of regional GNSS application, many researchers have attempted to develop regional ionospheric models based on the local GNSS network. Opperman et al. [10] used GPS observations to develop a regional model over South Africa for near real-time TEC mapping, and the model was based on the adjusted spherical harmonic method. Liu et al. [11] built a regional TEC model over China by the spherical cap harmonic analysis, and the spherical cap harmonic model had a better TEC mapping accuracy with smoother residual distributions in both temporal and spatial domains; the prediction accuracies in the following three days ranged from 2.5 to 4.5TECU. Continually, the spherical cap harmonic analysis was used to construct the Australian TEC model successfully [12]. The Australian TEC model not only simulates the daily and seasonal variabilities of electron content, but also captures the ionosphere dynamic evolution under severe geomagnetic storms. Razin et al. [13] used a support vector machine (SVM) to develop a spatial-temporal model for predicting TEC over Iran; the model based on the observations of 37 GPS stations had a good application under intense solar activity.

With the rapid development of computer performance and deep learning algorithms, the artificial neural network (ANN) has proven to be a powerful tool in modeling the earth dynamic system as well as building an ionospheric prediction system over a single station, small region or even at a global scale [14,15,16,17]. The ANN technique was first utilized to construct a critical frequency of the F2 layer (foF2) over the Grahamstown station (33° S, 26° E) based on the dataset during 1973–1983 [18]; the root-mean-square error (RMSE) of the daily foF2 prediction reached 0.95 MHz. Afterward, the advanced ANN technique was widely applied in modeling TEC variation of the regional ionosphere over low latitudes [19], China [20], Iran [21], Japan [22] and the global region [23]. Many studies have demonstrated that the ANN-based TEC models have a better application in small-scale areas than the IRI model. Tulunay et al. [24] developed a data-driven neural network model of one hidden layer and several neurons for predicting ionospheric TEC over Europe, namely the Middle East Technical University Neural Network (METU-NN). Habarulema et al. [25] used the observations of 10 GPS receivers in South Africa to train a feed-forward neural network, and the result showed that the prediction accuracy of a national model led to improvement compared to the single station model. Ruwali et al. [26] utilized GPS observations from Bengaluru, Lucknow and Guntur to develop a hybrid forecasting model that combines convolutional neural network (CNN) and long-short-term memory (LSTM). The LSTM-CNN model outperformed other ionospheric deep learning algorithms, achieving a minimum RMSE of 1.5TECU and a high correlation coefficient (CC) of 0.99. In addition, Reddybattula et al. [27] also established a TEC model over Bangalore (77.57° E, 13.03° N) based on LSTM that demonstrated a minimum RMSE of 1.6149 and a CC of 0.992. Shi et al. [19] utilized the measurements from 187 GNSS stations and the Constellation Observing System for Meteorology, Ionosphere and Climate (COSMIC-2) to forecast TEC variations in low latitudes; the new model may be a good reference for the structure and variation of the ionosphere over low latitudes. Dabbakuti et al. built an ionospheric model based on the ANN using singular spectrum analysis (SSA) for forecasting the TEC variations in low-middle latitudes under storm conditions, and the results indicated that the precision of the proposed SSA-ANN model improved 43.82% compared to the standard principal component analysis-ANN model [28].

Multiple regional TEC prediction models have concluded that the performance of the empirical model based on the ANN technique is heavily dependent on observations. In some regions with a dense GNSS network, the prediction model was competitive with the global ionospheric map (GIM). The GNSS network in Japan is one of the densest networks in the world; the Geospatial Information Authority of Japan (GSI) operates GNSS-based control stations that cover Japanese archipelago with over 1300 stations at an average interval of about 20 km for crustal deformation monitoring and GNSS surveys. Maruyama built a regional TEC prediction system over Japan using the surface harmonic functional expansion, and the coefficients of the expansion were modeled by the ANN technique. It is well known that the geomagnetic activity has a serious effect on the TEC variations [29], while the TEC variation associated with geomagnetic storm was not considered by the model. Therefore, in the study, we develop a new ionospheric TEC model over Japan by a multi-layer perceptron neural network based on the dual-frequency observations from hundreds of GNSS stations, especially focusing on improving the predicted capability under a severe space environment. In addition, the performance of the new model is evaluated during various space environments compared to the IRI-2016 and the thermosphere–ionosphere electrodynamics general circulation model (TIE-GCM), respectively.

2. Materials and Methods

2.1. Materials

The GPS Earth Observation Network (GEONET) operated by the Geographical Survey Institute (GSI) is a nationwide GPS array of Japan. The GEONET that consists of 1200 GPS permanent stations is one of the largest GPS networks in the world. In this study, the dual-frequency measurements derived from 839 GPS receivers were utilized to construct the Japanese TEC model by a multi-layer perceptron neural network; the geographic location of 839 GPS receivers is shown in Figure 1.

The GPS observations can be obtained from the website https://www.gsi.go.jp/ENGLISH/geonet_english.html (accessed on 15 March 2022). In this study, the original GPS observations were first checked for quality control by the translation, editing and quality checking (TEQC) software. The TEQC software, released by the university-governed consortium, facilitates geoscience research and education using geodesy (UNAVCO). It is a powerful tool in solving pre-processing problems of observations of multiple GNSS systems, including GPS, GLONASS, Galileo, Beidou, etc. [30]. It is provided by the link https://www.unavco.org/software/data-processing/teqc/teqc.html (accessed on 15 March 2022). Second, the slant total electron content (STEC) along the path from receiver to satellite was estimated as follows:

$$I_i = \frac{40.3}{f_i^2}{\displaystyle {\int }_{s}Ne}\mathrm{ds}= \frac{40.3}{f_i^2}STEC$$

(1)

where I_i is the ionospheric delay for the frequency f_i of the signal L_i. Then, the STEC was calculated by the geometry-free linear combinations of pseudo-range and carrier-phase measurements [31]. Finally, the STEC was converted to the vertical TEC (VTEC) through the Equation (2):

$$VTEC=STEC\times \cos (\arcsin (\frac{R\times \sin z}{R+H}))$$

(2)

where R is the Earth radius (6371 km), and z is the zenith angle of satellite. H is the selected altitude of the thin shell, and it was selected as 350 km in this study. It is noted that TEC signifies VTEC throughout the paper. For reducing the TEC error associated with the differential Code Biases (DCB) of multiple GNSS systems, only the observations derived from the GPS system were used in the estimation.

2.2. Multi-Layer Perceptron Neural Network

The ANN technique is a powerful method in learning the nonlinear relationships between various kinds of variables in depth, and it is widely applied in modeling earth systems, image recognition, data classification, etc. [32]. The ANN technique contains many types of algorithms, including Back Propagation Network (BP), Convolutional Neural Networks (CNN), Radial Basis Function (RBF), Recurrent Neural Network (RNN), etc. [33]. Among them, the BP network is one of the most classical and successful algorithms [34]. Usually, a BP network has two stages: feed-forward and back-forward. In a feed-forward process, a BP network learns the nonlinear relationship between various inputs and predicts a corresponding value, then the predicted residual is estimated by the difference between target and prediction. In a back-forward process, the predicted residual is transferred from the output layer to the input layer, and then the input weights are adjusted.

A simple BP network consists of one Input layer, one hidden layer and one output layer. Compared to the BP network with one hidden layer, multiple hidden layers could accelerate the learning efficiency and increasing the prediction accuracy partly. Therefore, a BP network with multiple hidden layers is utilized to construct an ionospheric model for predicting the TEC variations over Japan, namely Japanese Ionospheric Model (JIM). Previous studies have found the ionospheric TEC had a distinct daily period that was dominated by the Earth’s rotation [35]. Therefore, the day of year (DOY) and universal time (UT) were selected as the time parameters of the JIM. The parameters DOY and UT were normalized as sin(2pi*DOY/365.25), cos(2pi*DOY/365.25), sin(2pi*UT/24) and cos(2pi*UT/24). Moreover, the ionospheric TEC is inversely proportional to the geographic latitude that is associated with the intensity of solar irradiation. Thus, the geographic longitude (Lon) and geographic latitude (Lat) were chosen as the location parameters of the JIM. Furthermore, the solar parameters, including solar wind speed (Vsw) and 10.7 cm solar flux (F10.7), and the geomagnetic indices that contain Dst and Kp were also selected. The solar-geomagnetic parameters were obtained from the Space Physics Data Facility, Goddard Space Flight Center (GSFC) via the website (https://omniweb.gsfc.nasa.gov/form/dx1.html (accessed on 15 March 2022).

The training of the ANN model was stopped when the RMSE reached its minimum value. The optimal structure of the ANN model consisted of three hidden layers with corresponding neurons of 30, 28 and 24, respectively. The structure of the JIM model is shown in Figure 2. The structure of the JIM was adjusted for thousands of times, and the statistical result of prediction errors with different structures is shown in Figure 3. The result indicates that the minimum RMSE of the JIM model reached 1.4TECU.

3. Results

3.1. Evaluation of the JIM’s Accuracy for Different Sample Datasets

This study utilized GNSS observations from the GEONET between 2009 and 2021 to develop the JIM model. Specifically, days when the solar wind speed (Vsw) exceeded 600 km/s or the Dst value drops below −50 nT were defined as “space disturbing periods.” In summary, the quiet-time JIM model was constructed using GNSS data spanning more than 11 years, whereas the storm-time JIM model utilized data spanning approximately 560 days. In the training stage, the samples were divided into three parts: the training set, validation set and test set with the corresponding 70%, 15% and 15% percent of the total samples, respectively. In addition, the observations under the severe solar or geomagnetic activities were also trained as a “Storm” option. In the prediction stage, if the space-geomagnetic indices exceeded the threshold value (Vsw > 600km/s or Dst < −50nT), the “storm” option would be activated automatic. The prediction performances of the JIM model under quiet or “Storm” conditions were evaluated, and the accuracies in each sample dataset are shown in Figure 4.

The left panels describe the prediction performances of the JIM model in training set, validation set, test set and whole set under quiet space environment. The color bar signifies the correlation coefficient between prediction values and targets, while the lighter color signifies a higher correlation coefficient. The result indicates that the JIM had a good performance in the simulations. In four datasets, all of the correlation coefficients almost reached 0.98. In the training stage, Figure 4a shows that the accuracy of the JIM model for the training dataset was highest with a RMSE of 1.4974TECU. In the other three sample datasets, the RMSEs of the JIM were 1.4985, 1.5021 and 1.4982TECU, respectively.

Similar to the left panels, the right panels signify the prediction performances of the JIM model under “Storm” condition. The result shows that the JIM model had a better performance under the severe space environment, and the predictions and targets surrounded the regression line closely. All of the correlation coefficients in Figure 4e–h exceeded 0.99, and the minimum RMSE reached 0.96TECU. Figure 4 demonstrates that the JIM model based on the multi-layer perceptron neural network remarkably enhanced the prediction capability of TEC variation over Japan under the complex space environment. Especially under disturbed space events, the prediction performance of the JIM model remained at a stable level.

3.2. Evaluation of the JIM’s Prediction Ability in Reconstructing TEC Map

In order to investigate the capability of the JIM model in reproducing the spatial maps, the averaged latitudinal variations of TEC during the equinoxes and solstices of 2016 (solar maximum year) were simulated, as shown in Figure 5. The latitudinal TECs were calculated as the means of the electron contents over 839 GNSS stations with an interval of 1° latitude. It was found that the TEC had a remarkable hourly variation, which expresses that the TECs during local time (LT) 10–14 were larger than that in the evening. The amount of electron contents in the equinoxes was obviously larger than that in the solstices. In the spring equinox, the TEC remained in a low level (<10TECU) at dawn. From LT08, the TEC abruptly enhanced and reached a maximum value of 35TECU at noon. The TEC gradually decreased from LT14 and dropped to a quiet level (<10TECU) after LT19. Moreover, the scope of high TEC had a remarkable extension, and we found that the high TEC extended to 45° N in the spring equinox. Figure 5c shows that the TEC variation in the autumn equinox agreed well with the spring equinox, except for the amount of electron content. The TEC peak was centered at LT13 with a maximum value of ~23TECU, and the high TEC reached 37° N. The hourly variations of TEC in the solstices were consistent with the equinoxes, but the duration of high TEC lasted a bit longer or shorter. For example, in the summer solstice, the high TEC increased from LT06 and dropped after LT21, which is associated with the longer sunlit illumination time. The condition was inversed in the winter solstice, as shown in Figure 5d. At that time, the subsolar point moves to the Tropic of Capricorn, which reduces the sunlit illumination time over Japan. The hourly TEC variations simulated by the JIM model in four seasons agree well with [21], which demonstrates the JIM has a strong ability in reconstructing the spatial TEC maps over Japan.

3.3. Temporal-Spatial Variations of the JIM-TEC

In order to validate the prediction accuracy of JIM model in a long-term scale, the daily variations of TEC over station STK2, 0203 and TSKB during DOY18-47, 2016 were estimated, see the left column in Figure 6. Meanwhile, the corresponding TEC timeseries were also simulated by the JIM model, see the middle column in Figure 6. Finally, the prediction residual was estimated by TEC(resi) = TEC(pred)—TEC(obs), as shown in the right column of Figure 6. The panels in the left column show that most of high TECs accumulated in the period of UT01-05 during DOY18-47, 2021, and the values were within 12-16TECU. However, two diurnal peaks of TEC occurred in DOY24-28; the first peak appeared in UT01-03, and the second peak happened in UT05-07. The left column shows that the JIM model not only reconstructed the daily change of TEC, but also successfully simulated the double-peak structure of TEC during DOY24-28. The TEC residuals were proportional to the amount of electron content, which shows that most of significant TEC residuals accumulated in UT01-06. Figure 6c shows that both positive and negative residuals were observed over station STK2, and the absolute amplitude of TEC residual was within 4TECU. The maximum prediction residual appeared in UT00-02, DOY24-26, and in other periods, the absolute TEC residuals drop to 2TECU. The TEC residuals over station 0203 and TSKB were smaller, and the magnitude of TEC residuals during UT00-06 ranged from −2 to 2TECU. In other hours, the prediction residuals were lower than 0.5TECU.

To further evaluate the spatial prediction performance of the JIM over the whole Japan region, the TECs over 839 GNSS stations at UT06, UT12 and UT18 on 22 March 2021 were estimated by the GEONET observations. In addition, the TEC maps were also simulated by the JIM, GIM and IRI-2016 model, respectively. In the study, the GIM products were provided by the Center for Orbit Determination in Europe (CODE) via the link (http://ftp.aiub.unibe.ch/CODE/ (accessed on 15 March 2022). The accuracy of CODE-GIM was evaluated using observations from the Crustal Movement Observation Network of China (CMONOC), which is located at latitudes similar to those of Japan. The results indicate that the RMSE of CODE-GIM was about 2-4TECU under the magnetic quiet circumstances, and for the magnetic disturbances circumstances, the RMSE increased to 5.7TECU [36]. The IRI-2016 model has demonstrated a high level of proficiency in predicting ionospheric variability and accurately depicting the spatial characteristics of global ionospheric convection, including the equatorial ionization anomaly (EIA), Weddell Sea anomaly and winter anomaly, among others. This assertion is supported by previous research [7,37]. For the purpose of this study, the observed solar-terrestrial indices of F10.7 and Dst were utilized as inputs for the IRI-2016 model to generate TEC maps. The F10.7 and Dst indices were obtained from GSFC (https://omniweb.gsfc.nasa.gov/form/dx1.html (accessed on 15 March 2022).

The comparative results of three kinds of prediction maps are shown in Figure 7. Figure 7 shows distinct TEC latitudinal variations according to the amount of electron content, which expresses that the TEC was larger than 16TECU in 26–36° N latitude, following the 12-14TECU within 36–40° N latitude, and the TEC dropped to ~10TECU in 40–45° N. Figure 7b,c indicate that the TEC maps predicted by JIM and GIM agreed well with the actual TEC, while the IRI-2016 model significantly underestimated the TEC in the southern region of Japan, such as Kyushu, Shikoku, etc. At UT12 and UT18, the phenomenon was more evident. Both JIM and GIM succeeded to reconstruct high-accuracy TEC maps compared to GEONET observations; among them, the accuracy of the GIM was better, and the JIM underestimated the TEC with a magnitude of 1-2 TECU. However, the prediction accuracy of the IRI-2016 abruptly decreased, and the maximum prediction residual reached 4-5TECU. For example, the TEC derived from the GEONET observations was about 12TECU at UT12, while the corresponding TEC simulated by the IRI-2016 dropped to ~7TECU.

3.4. Prediction Performances of the JIM under Various Complex Space Environments

Prediction accuracy is the most important indicator in evaluating the performance of a new model, which has been verified partly in the Section 3.1–3.3. Because the GIM is estimated by the observations of hundreds of GNSS receivers globally, the GIM file has a delayed release by 1–2 days. Therefore, the GIM is useless in the application of real-time ionospheric delay correction. In conclusion, the real-time ionospheric delay is only corrected by the means of differential technique or ionospheric models. The physics-based TIE-GCM model released by the High-Altitude Observatory (HAO) of the National Center for Atmospheric Research has a distinct advantage of seasonal variations and spatial resolutions over empirical models. The model solves the three-dimensional momentum, energy and continuity equations for neutral and ion species at each time step using a semi-implicit, fourth-order, centered finite difference scheme on each pressure surface, and plays an important role in numerical simulation and physical explanation of ionospheric dynamic convection. In addition, as we know, severe space events usually trigger remarkable ionospheric perturbations, such as travelling ionospheric disturbances (TIDs), ionospheric scintillation, plasma bubbles, etc. [16,31,38]. Under this condition, the disturbed ionosphere affects the radio propagation and reduces the accuracy of positioning and navigation. Therefore, it is necessary to evaluate the prediction capability of the JIM model under complex space environments compared to the TIE-GCM and IRI-2016 model, and the results are shown in Figure 8 and Figure 9.

In Figure 8, the TECs over 839 GNSS stations during 22 March 2020 are estimated as reference values. During the period, the space environment was quiet, the solar wind velocity was lower than 470 km/s and the Dst value was larger than −20nT. The left panels show the Japanese TEC maps at UT06, UT12 and UT18, respectively. Figure 8a shows that the TEC gradually depleted from 16TECU in 27° N to 12TECU in 45° N at UT06. Both the JIM and IRI-2016 overestimated the TEC over the southern region of Japan, while the TIE-GCM model underestimated it compared to references. For example, the reference TECs over the Kyushu, Shikoku and Chugoku regions were about 12-13TECU, the corresponding JIM-TECs ranged from 13-15TECU and the maximum IRI-TECs exceeded 16TECU. Meanwhile, the TEC maps predicted by the TIE-GCM were within 8-10TECU. The comparative result indicates that the prediction performance of the JIM model is better than other two models.

The prediction results at UT12 and UT18 are consistent with that at UT06. Figure 8f shows that the JIM model estimated the TEC with a magnitude of 2TECU at UT12, and the main prediction residuals were located at Kyushu and Shikoku. However, both the TIE-GCM and IRI underestimated the TECs over the whole Japan region with a magnitude of 3-4TECU. Among them, the prediction performance of the TIE-GCM model was worse. The maximum prediction residuals appeared in the northern region (Tohoku and Hokkaido) of Japan with a magnitude of 4TECU, as shown in Figure 8g. The prediction residuals between the three models were more distinct at UT18. The prediction performance of the JIM was stable, and the residual was within 2TECU. However, the prediction performances of the TIE-GCM and IRI were significantly reduced, and the prediction accuracies abruptly dropped to 4-5TECU, especially in Figure 8k.

The X9.3 solar flare, the strongest eruption since 2005, erupted at UT11:53, 6 September 2017, during which large amounts of enhanced energy and momentum emitted from coronal mass ejections (CMEs) triggered severe thermospheric-ionospheric (TI) perturbations. The TI perturbations lasted several hours, followed by two geomagnetic storms occurring on 7–8 September 2017, when the solar wind speed exceeded 800 km/s and the Dst dropped to −142 nT [39]. Therefore, the TEC maps estimated by the observations of GEONET on 8 September 2017 were selected as references, and the corresponding TEC maps simulated by the JIM, TIE-GCM and IRI-2016 are shown in Figure 9.

The top panels show that the JIM and IRI underestimated the TEC at UT06. The JIM underestimated about 2TECU in Kyushu, and for IRI-2016, the prediction residual increased to 3–4TECU. However, the prediction performance of the TIE-GCM model was not satisfied. Figure 9c shows the TECs over the southern region of Japan enhanced to 24TECU, and the maximum prediction residual reached 6–8TECU. In addition, the prediction TEC over the Hokkaido was also 4–6TECU larger than the references. At UT12, the JIM-TECs were almost consistent with the actual values, and the prediction residual was smaller than 2TECU. The TIE-GCM overestimated the TECs over Kyushu, Shikoku and Chugoku with a magnitude of 3–4TECU, but underestimated by about 2TECU over Hokkaido. On the contrary, the IRI-2016 model remarkably underestimated the TECs over the whole Japan region, and the averaged prediction residual was about 4TECU. Similar to Figure 8, the JIM-TECs agreed well with the reference at UT18 under severe geomagnetic storms, and both TIE-GCM and IRI failed to reconstruct the precise TEC maps.

The RMSEs of the JIM, GIM, TIE-GCM and IRI-2016 under quiet and severe geomagnetic activities were calculated, and the statistical result is shown in

Table 1. RMSEs of the JIM, GIM, TIE-GCM and IRI-2016 at UT06, UT12 and UT18 under quiet and severe geomagnetic days. Reference: GNSS-TEC.

	Quiet Day (22 March 2020)			Storm Day (8 September 2017)
	UT06	UT12	UT18	UT06	UT12	UT18
JIM	2.01	0.75	0.60	1.51	0.88	0.89
GIM	0.93	0.37	0.65	1.58	0.62	0.74
TIE-GCM	3.15	3.28	4.87	8.83	2.79	5.34
IRI-2016	4.62	3.12	3.32	3.80	5.32	3.74

. We can find that the prediction accuracy of the JIM was nearly equal to the GIM; in some conditions, the JIM had a more competitive edge. For example, in the quiet space condition, the prediction RMSEs of the JIM at UT06, UT12 and UT18 on 22 March 2020 were 2.01, 0.75 and 0.60TECU, respectively. The prediction accuracies of the GIM were higher at UT06 and UT12 with the amplitudes of 0.93 and 0.37TECU, respectively. However, at UT18, the RMSE of the GIM reached 0.65TECU, which was larger than that of the JIM with an amplitude of 0.05TECU. This phenomenon was also observed at UT06 on 8 September 2017 when a strong geomagnetic storm occurred. At this moment, the prediction RMSE of the GIM was 1.58TECU, while the corresponding RMSE of the JIM was only 1.51TECU. In other moments, the prediction RMSEs of the JIM were larger than the GIM, with amplitudes ranging from 0.26 to 0.15TECU. Compared to the JIM, the prediction performances of the TIE-GCM and IRI-2016 were worse under any space environment. The RMSEs of the TIE-GCM and IRI-2016 were 3–4-times larger than that of the GIM, which is consistent with the results of Figure 7, Figure 8 and Figure 9. During the quiet day, the RMSEs of the TIE-GCM and IRI-2016 ranged from 3 to 5TECU. The RMSEs suddenly increased under severe geomagnetic activity. For example, the RMSE of the TIE-GCM increased to 8.83TECU at UT06, 8 September 2017, and the RMSE of the IRI-2016 increased to 5.32TECU at UT12. The result of Table 1 demonstrates that the prediction capability of the JIM is significantly superior to the TIE-GCM and IRI-2016. The JIM model not only predicted the TEC variations during quiet days precisely, but also successfully simulated the ionospheric dynamic transportation under geomagnetic activities.

The results in Figure 8 and Figure 9 show that the JIM has a powerful predictability not only in quiet space condition, but also under severe geomagnetic storms. To further to investigate the predicted accuracy of the JIM depending on the season and solar-terrestrial environments, the RMSEs of the TECs simulated by the JIM, IRI-2016 and TIE-GCM during different months, as well as the solar wind speeds and geomagnetic Dst values, were calculated, as shown in Figure 10. Figure 10a shows the monthly TEC RMSEs based on JIM, IRI-2016 and TIE-GCM from January to December; the prediction RMSEs of the JIM-TECs were much smaller than that of the IRI and TIE-GCM in each month. Generally, the RMSEs of the JIM-TECs are smaller than 2TECU in each month. Among them, the maximum RMSE occurred in March with a magnitude of 1.86TECU, and the minimum RMSE was about 1.26TECU in January. Both the corresponding prediction RMSEs of the IRI and TIE-GCM exceeded 2TECU, and the performance of the IRI-2016 was better than the TIE-GCM. In January, the RMSEs of the IRI and TIE-GCM were 2.7 and 3.6TECU, respectively. Then, the prediction residuals of two models gradually increased and reached a maximum level in March with amplitudes of 4.05 and 4.56TECU, respectively.

Figure 10b,c show the prediction performances of the JIM, IRI and TIE-GCM under various intensities of solar-geomagnetic activities. It was found that the prediction RMSEs of three models were proportional to the intensity of space condition, and the JIM performed much better than the IRI-2016 and TIE-GCM when the “Storm” option was activated. For example, the RMSE of JIM-TECs was about 0.8TECU when the solar wind speed was lower than 500 km/s, and the corresponding RMSEs of the IRI and TIE-GCM were 2.6 and 2.75TECU, respectively. With the increment of solar wind speed, the prediction residuals of the IRI and TIE-GCM abruptly enhanced. The maximum RMSEs of the IRI and TIE-GCM reached 4.85 and 5.69TECU when the solar wind speed exceeded 800 km/s. However, the prediction accuracy of the JIM aided by the “Storm” option under severe solar flares was controlled at a satisfied level. Even the solar wind speed reached 800 km/s, the RMSE of the JIM-TEC was only 1.95TECU.

Similar to the solar wind speed, the RMSEs of the JIM, IRI and TIE-GCM were also proportional to the intensity of geomagnetic activity. Figure 10c shows that the RMSEs of TECs predicted by the JIM, IRI and TIE-GCM were 1.86, 4.7 and 5.36TECU when the geomagnetic Dst index dropped to −210 nT, respectively. As the geomagnetic activity gradually recovered to a normal background level, the predicted accuracy of three models enhanced significantly. For example, the RMSEs of three models decreased to 0.85, 2.1 and 3.23 TECU under the quiet geomagnetic condition (Dst > −30 nT).

4. Discussion

The comparative results in Section 3 demonstrate that the JIM model based on artificial neural network has a powerful capability in capturing ionospheric spatial-temporal features over Japan. The ANN model with the “Storm” option not only predicted Japanese ionospheric variation with higher accuracy under quiet solar-geomagnetic activities but also had a better prediction performance during severe space events compared to traditional empirical and physical models, including IRI-2016 and TIE-GCM. Figure 10 shows that the minimum RMSE of the JIM during the quiet space environment (Vsw < 500 km/s and Dst > −30 nT) was about 0.8TECU, and even in disturbed space events, the maximum predicted RMSE was smaller than 2TECU. Ansari et al. [40] used long-term GNSS measurements from 01 January 2011 to 31 December 2018 to build an ionospheric TEC model over Southwest Japan through a modified linear timeseries model (LTM). Similar to the JIM, the LTM also included the geomagnetic indices (Northern Polar Cap (PCN) and Disturbances (Dst)) and the solar-planetary geomagnetic disturbance indices (F10.7 and Ap). The comparison results show that the RMSEs between the observed and the TEC predicted by the LTM varied from 1.57 to 1.66TECU, while for the auto-regressive moving average (ARMA) model and the global ionospheric map (GIM), the predicted RMSEs were 1.84–1.93TECU and 9.67–9.84TECU, respectively. In addition, Maruyama [22] constructed a regional reference model of TEC for predicting two-dimensional (time vs. latitude) TEC maps over Japan. The time-latitude variation average was first determined by the surface harmonic functional expansion, and then the coefficients of the expansion were modeled by a neural network with input parameters of the season (day of year) and solar activity (F10.7 and sunspot number). Though the model succussed to capture the time-latitude distribution of TEC map, the predicted accuracy was relatively low (3.29TECU) and the storm-induced TEC variations were not considered. Furthermore, Li et al. [12] built an Australian Ionospheric Model (AIM) based on a combination method of artificial neural network and spherical cap harmonic analysis (SCHA). Due to the sparse GNSS stations located at the central vast land of Australia, the work first built an AIM model for predicting TEC values over the area that lacks observations; then, the observations of the Australian regional GPS network (ARGN) and the prediction TECs were utilized to construct the AIM. During the equinoxes and solstices of 2017 (solar minimum year), the RMSEs of TEC residuals were 2.16, 1.57, 1.68 and 1.98TECU, respectively. The prediction accuracy of the AIM was slightly worse than the JIM in this study. Compared to the geographic location of the GNSS stations in ARGN, the location of the GNSS stations in GEONET is significantly denser. Therefore, we believe that adding sufficient observations is a solution to replace other mathematical methods to improve the prediction performance of the ANN model. The above results demonstrate that the JIM model developed by a multi-layer perceptron neural network not only had a higher predicted accuracy than other ANN model, but also performed better under any space conditions compared to other methods, including LTM, ARMA, the surface harmonic expansion method, etc.

5. Conclusions

In this study, the long-term observations of 839 GNSS stations were used to constructure a regional reference model (JIM) for predicting TEC variations over Japan. Large amounts of GNSS observations were trained by a multi-layer perceptron neural network, as well as a “Storm” option. The JIM model had a good prediction performance under complex space environments, not only in quiet solar-geomagnetic activities but also in severe space events when the “Storm” option was activated. This study draws the following main conclusions:

(1)

The whole samples were divided into the training set, validation set and test set for learning the JIM. Under the quiet space condition, the correlation coefficients between targets and predictions for three datasets were 0.98, 0.97993 and 0.97994, and the corresponding RMSEs of prediction residuals were 1.4974, 1.4985 and 1.5021TECU, respectively. The performance of the JIM was better during severe space events; the correlation coefficients for three parts all exceeded 0.99, and the corresponding RMSEs were 0.96027, 0.95356 and 0.95767TECU, respectively.

(2)

The JIM had a strong capability in reconstructing the two-dimensional (time vs latitude) TEC maps over Japan. The JIM was successful in reproducing the spatial TEC maps during equinoxes and solstices, and the TEC maps had evident hourly and seasonal variations. The maximum TEC appeared in the spring equinox, following the autumn equinox, and the minimum values occurred in solstices. Moreover, the TEC timeseries simulated by the JIM were nearly consistent with the targets over GNSS stations STK2, 0203 and TSKB. Most of TEC residuals accumulated in UT01-06 with a maximum magnitude of 4TECU, while in other moments, the averaged magnitude of TEC residuals was lower than 1TECU.

(3)

The JIM had a perfect prediction performance under various kinds of complex space environments. During the 2021 spring equinox (Vsw < 600 km/s and Dst > −30 nT), both the predictions of JIM and GIM agreed well with the target TECs. The JIM usually tended to underestimate the TEC with a magnitude of 1-2TECU, and in some moments, the JIM had a more competitive edge than the GIM. Under severe geomagnetic storm on 8 September 2017, the performance of the JIM remained at a stable level. The RMSEs of the TEC residuals of the JIM at UT06, UT12 and UT18 were 1.51, 0.88 and 0.89TECU, while the corresponding RMSEs of the TEC residuals simulated by the IRI-2016 and TIE-GCM were 3–4 times larger than that of the JIM. Moreover, the TEC residuals had an evident monthly variation; the maximum residual occurred in March and April, and the minimum residual appeared in December. Furthermore, the magnitude of TEC residual was proportional to the solar wind speed and was inversely proportional to the geomagnetic Dst value. Even in a severe disturbed space environment, the TEC residual of the JIM was still lower than 2TECU, while the corresponding residuals for the IRI-2016 and TIE-GCM exceeded 5TECU.