1. Introduction
The ionosphere is a region of the partially ionized upper atmosphere [
1]. The ionosphere consists of several layers, which are characterized by enhanced levels of ionization. The innermost D-layer is located at heights from 50 to 90 km, the E-layer at 90–150 km, and the outermost F-layer extends from above 150 km. Radio waves propagate through the ionosphere at different group velocities and split into different wave modes, depending on frequency, the electron density, the magnetic field, etc. Various experimental techniques are used for measuring the profiles of electron density (it is also called the electron content or EC) in the ionospheric layers, including ground-based ionosondes, incoherent scatter radars and the spaceborne radio occultation technique.
The GPS radio occultation (at present, it is known as GNSS RO) is a remote sensing technique for vertical profile information on EC in the ionosphere from the height of a low-Earth’s orbit (LEO) satellite down to the bottomside, e.g., as shown in [
2,
3,
4]. The technique of radio occultation is based on the detection of a change in radio signals as they pass from Global Positioning System (GPS) or Global Navigation Satellite System (GNSS) satellites, located at altitudes of about 20,000 km, through the ionosphere to LEO satellites, located at altitudes of several hundred km [
5]. Electromagnetic waves of different frequencies pass the ionized medium at different speeds. The dispersion depends on the amount of ionization. Hence, the electron density can be related to the phase delay for radio signals of different frequencies by using an Abel inversion [
6]. The Abel transform contains a bias term (constant along each GNSS–LEO ray trace). In the real RO measurements, the ray traces vary substantially such that the bias is not constant, which makes the problem even more complicated. Various techniques are used to eliminate the bias [
7,
8]. In routine GNSS RO measurements, the effect of bias term cannot be always eliminated. This results in unexpected negative values in the resultant EC profiles in the bottomside ionosphere [
9,
10].
Recently, two constellations of six satellites FORMOSAT-3/COSMIC (FS-3/C-1) and FORMOSAT-7/COSMIC-2 (FS-7/C-2) have been launched to low-height orbits (below 800 km). Using a GPS occultation experiment (GOX), FS-3/C-1 measures EC profiles from 60 to 800 km at high inclination orbits [
11]. The comparison of the FS-3/C-1 RO technique with the standard techniques of ground-based digisonde and incoherent scatter radars at middle latitudes demonstrated very good agreement between the EC profiles obtained by different techniques, both in the bottom and topside ionosphere [
12]. The mission FS-7/C-2 has a low inclination. It is equipped with a TriGNSS radio occultation system (TGRS) that receives the signals from numerous GNSS satellites [
13]. That makes it possible to scan the low-latitude ionosphere in detail and provides up to 4000 RO profiles per day.
The large amount of EC profiles requires automatic identification of the ionospheric layers. However, the bias problem is a serious obstacle. Namely, it is difficult to establish a firm threshold of EC for each ionospheric layer because the zero level is different for different EC profiles. Alternative techniques are required in order to eliminate the bias problem. One of them is a technique of deep learning for the classification and identification of the ionospheric layers. Recently, this technique was applied to the ionograms scanned by the ground-based ionosondes [
14]. It was shown that a deep-learning method for ionogram automatic scaling (DIAS) can distinguish the signals of ionospheric layers with very high performance. The DIAS method is based on a U-shape structure, which consists of an encoder network and a decoder network. The technique provides the F-score of ~0.91, in contrast to F~0.81, which is provided by the standard method of Automatic Real-Time Ionogram Scaling (ARTIST, e.g., [
15]).
In the present paper, we propose a new approach for the identification and classification of the ionospheric layers. The approach is based on the application of machine learning for the analysis and modeling of EC height profiles, acquired from the spaceborne RO GPS technique in the FS-3/C-1 mission. This technique makes it possible to automatically treat the great amount of EC data for the determination of the key ionospheric layers and avoid the bias problem. It can also be easily applied to the data acquired from the FS-7/C-2 mission.
The study is organized as follows:
- -
Section 2 is a description of the experimental data. Six classes of the EC height profiles are introduced for the modeling process.
- -
Section 3 describes the two convolutional neural network (CNN, [
16,
17]) models used in this study for classification of the EC profiles.
- -
Section 4 demonstrates the result of the CNN model applications for classification of the EC profiles and determination of the key ionospheric layers, such as E (including sporadic) and F1, F2 and F3 layers.
- -
Section 5 is a discussion and comparison of the results obtained from different CNN models.
- -
2. Experimental Data
We use the data of the EC height profiles acquired from the GPS RO experiment onboard the FS-3/C-1 mission in the range from 60 to ~800 km [
18]. The profiles are measured with 1 sec time resolution, such that the height step varies around 1 to 2 km. For the modeling process, the profiles are interpolated with omni-distant steps of 2 km in the range of heights from 80 to 600 km. As a result, the interpolated EC profile consists of 261 points. We selected low-latitude profiles in a 20-degree vicinity of the Taiwan region located under the bulge of the equatorial ionization anomaly (EIA). This region is characterized by strong variability in the ionospheric layers, often including the appearance of the sporadic E-layer (Es) and high-altitude F3 layer.
Examples of the height profiles of interpolated electron density (IEC) are shown in
Figure 1. We distinguish six mutually exclusive kinds of the profiles, representing different ionospheric layers and features. Each panel in
Figure 1 shows the IEC height profile (left) and the trajectory of perigee of RO trace (right).
Figure 1a represents a class of purely the F2 layer (F2) with its maximum value located in the altitude range from 200 to 400 km. The second class (
Figure 1b) corresponds to profiles with a narrow sporadic E-layer (F2E), located at heights around 100 km. The third class (
Figure 1c) contains profiles with a large-scale enhancement of EC below 200 km that corresponds to the intense E-layer (F2ED). The fourth class (
Figure 1d) is an intense F1 layer, located below the F2 layer (F2F1) at heights around 200 km. The fifth class (
Figure 1e) represents profiles with the F3 layer developed above the F2 layer (F2F3). The maximum value of the F3 layer is weaker than that of the F2 layer. The sixth class (
Figure 1f) contains profiles distorted by strong noise fluctuations (Sc). The fluctuations might relate to ionospheric scintillations or other kinds of disturbances.
We analyzed around seven thousand RO profiles obtained during the years from 2011 to 2013 (solar maximum) and manually selected only those containing the prominent features described above (in total, 712 profiles). The number of samples selected for each class is listed in
Table 1. All the samples of each class are presented in
Supplementary Materials. They are used as a ground truth for the classification models (as shown in the next section). It should be noted that we excluded from the consideration bad or corrupted RO profiles. We also did not consider profiles that contained simultaneously several prominent layers, such as Es, F1 and F2 layers, in order to avoid misclassifications during the modeling process.
As one can observe in
Figure 1, the values of EC vary across a very wide range up to 20 (∗10
5/cm
3). In addition, the EC height profiles can have negative values due to the bias problem (as shown in panels (c) and (f) in
Figure 1), which is still unresolved in the GOX experiment onboard the FS-3/C-1 mission [
9,
10]. These features make it difficult to automatically identify the ionospheric layers in several thousand profiles per day, as well as cause issues in the identification of the noisy/disturbed profiles for further detailed analysis. To resolve this problem, we determine the type of each profile by recognition of its pattern in relative units, such that the exact values of EC in the data are not important in this classification task. Therefore, the bias problem is avoided. Considering this characteristic, we normalized each profile, dividing the EC by its maximum value for the model training. This makes it possible to reduce the difficulty of convoluting the diverse values.
We evenly split the six classes into training and test sets by the percentage 80% and 20%, respectively. Due to the insufficient size of the data, we did not further split the validation set from the training set. For the same reason, we applied cross-validation to reduce the bias of the model evaluation. The exact number of samples in the training and the test set, with respect to each class, is shown in
Table 2.
The procedure of cross-validation is described as follows. First, the whole dataset is randomly shuffled, and divided into five subsets. One of the subsets is used as the test set to evaluate the model, and the rest are used to train the model. There are five models trained, and their average performance (specified in
Section 4) are taken as the performance of the model.
3. Classification Models
For the classification of EC height profiles, we used two different CNN models. The models are constructed using the Keras library under the Tensorflow framework [
19]. The input layer of the models consists of 261 nodes, corresponding to the number of points in the EC height profiles interpolated from 80 to 600 km with the 2-km step. The output layer consists of six nodes, corresponding to the number of EC profile classes.
The first model is a simple CNN model based on the work of Kiranyaz et al. [
16]. This model is a basic and simple one-dimensional-CNN, with two types of different layers. It can be distinguished by two parts. As shown in
Figure 2, the first part is composed of 1D-convolutional layers and each of the convolutional layers is followed by an activation and max-pooling layer. The convolution kernel can extract features from the RO profile. The max-pooling layer down-samples the profile, which allows the next convolution kernel to extract more rough features. We choose to repeat it four times to extract the feature by convolution kernels. Then, the data of the feature map are sent into the second part, which is composed of fully connected layers, known as a multi-layer perceptron (MLP). In this part, the output array from the convolutional layers is firstly flattened into a one-dimensional-array, which is then imported to three fully connected (FC) layers and two dropout layers, which randomly sets the input to zero to reduce the overfitting [
20]. In the end of the MLP, the output is set to be six, corresponding to our six classes in the data. Combining these two parts, the CNN model contains 27,110 training parameters.
The second model is a fully convolutional network (FCN) model [
17], composed of 1D convolution kernels, batch normalization layers, max-pooling layers, a dropout layer and a fully connected layer. The batch normalization layer regularizes the input and can accelerate the convergence [
21]. In the model, 1D convolutions of kernel size
k are repeated by
r times. Each time, there are
n features extracted by the kernels. The convolution process is followed by a batch normalization layer and a max-pooling layer to form a pooling process. The pooling process is repeated by
d times, and is then flattened with a dropout rate of
p. The flattened neurons are then fully connected, and output to a vector of six elements using the soft-max function, representing the probabilities of each class. The class with the highest probability is used as the classification result. An illustration of the architecture of the FCN model is shown in
Figure 3. A hyper-parameter search based on the Bayesian optimization technique [
22] is performed to find the most suitable set of the hyper-parameters (
k,
n,
r,
d,
p). The search grid of
k is 7, 9, 11,
n is 8, 16, 24, 32,
r is 1, 2, 3,
d is 1, 2, 3, and
p is 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, respectively.
During each training epoch, the model is fed with a mini-batch from the training set. A categorical cross-entropy function is used to measure the loss between the ground truths and the inference of the model. The Adam optimizer [
23] is used to obtain the gradients, which update the parameters of the model, and reduce the loss. The learning rate to the optimizer is set to 0.001 initially, and is reduced by a factor of 0.1 if no improvement is found within 10 epochs. The training process is stopped early if the performance of the test set does not improve within 20 epochs. After the training, the five models are evaluated using the corresponding test set, and their performances are averaged. The maximum number of trials on the hyperparameters is set to 30. The trial with the best average performance among 30 trials is selected as the best performance of the FCN model.
4. Results
To measure the performance of the model, the following confusion matrix, accuracies, and average accuracy are used.
where
is the total number of samples for the i-th class in the k-th test set,
is the total number of samples classified as the j-th class in the k-th test set,
is the average of five confusion matrices for the models evaluated by the corresponding test sets,
is the average accuracy of the i-th class over the five test sets, and
is the average accuracy over the five test sets and the six classes.
Figure 4a presents a confusion matrix
for the predictions of our CNN model. By extracting the diagonal of the confusion matrix, we obtain the average accuracy of 0.76, 0.75, 0.85, 0.70, 0.33, and 0.22, respectively, for the classes F2, F2E, F2ED, F2F1, F2F3, Sc. The average accuracy for all classes is 0.6. For this model, the F2ED class has the best prediction accuracy and high accuracies (>0.7) can be found for the classes F2, F2E and F2F1. The model performs the worst when identifying the F2F3 and Sc classes, with an accuracy lower than 0.5. In
Figure 4a, we can obviously observe that many profiles of F2F3 and Sc types are predicted as F2E, F2ED and F2F1. In short, the model predicts mostly accurately the F2 layer and prominent bottomside layers D, E and F1. Very poor prediction is found for the F3 layer and disturbed profiles Sc.
Of the 30 trials on the hyper-parameter searching for the FCN model, we found that the combination (k, n, r, d, p) = (7, 24, 3, 3, 0.2) and a batch-size of 4 produces the highest
. The resulting FCN model has 229,446 parameters in total, and its performance on
is 0.69 across the six classes. The corresponding confusion matrix
is shown in
Figure 4b. The average accuracy is 0.86, 0.62, 0.80, 0.76, 0.39, 0.72, respectively, for each class. The five test sets show that the models can produce accuracies above 0.67, and up to 0.96 for F2, F2ED, and F2F1, and that the accuracy can range from 0.52 to 0.77 for F2E and Sc, 0.22 to 0.74 for F2F3.
From
Figure 4b, one can observe that the model has tendencies to predict the F2E class as Sc and to predict F2F3 as F2ED. There are no specific tendencies for the other classes. On average, the model has the best performance for classifying F2, and has poor performance for classifying F2F3. To provide additional information for interested readers, we also calculate the confusion matrix in terms of the exact numbers of samples for CNN and FCN, as shown in
Figure 4c and
Figure 4d respectively. The full list of the classification for the two models is presented in the
Supplementary Materials.
The performances of the CNN and FCN models are compared in
Table 3. The CNN model demonstrates better performances for the prediction of the F2E and F2ED classes. The classes F2, F2F1, F2F3 and Sc are predicted better by the FCN model. It should be noted that the FCN and CNN models have similar performances in the prediction of F2ED and F2F1 classes. However, the FCN model has much better performance in the classification of the F2 and Sc classes. Both models have low performance in the prediction of the F2F3 class.
In order to estimate the overall performances of the models, we also calculate the standards scores, such as recall rate, precision, binary accuracy, and F1-score [
24]. They are defined as follows.
where
are the true positive, false positive, true negative, and false positive predictions for the j-th class in the evaluation of the k-th subset. The scores are presented in
Table 3.
5. Discussion
We have found that both models demonstrate the best classification performance (around 0.8) for the prediction of the F2 and F2ED classes. The former one corresponds to a simple F2 layer, which is mostly prominent at most height profiles. Hence, the recognition of this feature is mostly easy. Another class with high performance is F2ED. This class is characterized by a prominent smooth enhancement of EC at altitudes below 200 km, in the region of the E and D layers. This large-scale structure is also very easy to recognize.
Both models demonstrate moderate accuracy (>0.6) in the identification of the sporadic Es layer (F2E class) and F1 layer (F2F1 class). The numbers in the confusion tables indicate that the two models either have a preference to classify F2E as Sc classes or to classify F2F3 as F2ED. This could be due to the loss of position information in the max-pooling [
25] or noisy RO structures. The lowest performance (<0.4) is found for the prediction of the F2F3 class. This class corresponds to the prominent F3 layer. The statistics of this class are relatively low (see
Table 1) because the F3 layer is relatively rare, which might be a reason for the low performance of the prediction. One must note that although the performance for the F2F3 class is poor, it is still much higher than the estimation of 0.17. On the other hand, the Sc class has the lowest statistics (see
Table 1). This class of small-scale fluctuations (scintillations) is predicted by the FCN model, with accuracy (0.72) much higher than that of the CNN model (0.22). The latter one, however, has the highest performance for the prediction of the large-scale structures of the F2E and F2ED classes. Hence, models with different architectures have different performances for the prediction of different scale structures.
Figure 5 demonstrates examples of misclassifications, which occur very often for each class. The F2 class can be mixed up with the F2F1 class, due to enhancements of EC at heights below 200 km (see
Figure 5a). As shown in
Figure 5b, the prominent E layer can often be misclassified as the Sc class. A large-scale increase in EC at low altitudes in the F2ED class can be interpreted as the F1 layer from the F2F1 class (
Figure 5c). Similarly, an increase in EC in the region of the F1 layer at heights above 100 km can be misclassified as the E layer from the F2ED class (
Figure 5d). Distorted EC profiles (Sc class) can be often mixed up with other classes in the case of prominent features in the EC height profiles, such as an enhanced E layer (see
Figure 5f). The prediction of a relatively weak F3 layer is difficult for both models and it is often misclassified as the F2ED class. In addition to poor statistics, this shortcoming can also be the result of a problematic input dataset. As one can observe in
Figure 5e, the prominent F3 layer at heights >400 km is accompanied by a weak EC enhancement in the bottomside ionosphere at heights <200 km. The latter is the distinct feature of the F2ED class. This situation is often observed, meaning that it is very difficult to select samples with the pure F3 layer.
Although adjusting the weighting in the loss function can improve the poor performance of the imbalanced class [
26], we found that the performance on the unaffected classes also decreased accordingly. As mentioned earlier, due to the size of the dataset, we apply the cross-validation technique to reduce the bias in the model evaluation. However, it is still computationally expensive to analyze the change in error in the model evaluation, with respect to the size of cross-validation [
26,
27,
28]. As a result, further adjustments on the weightings of the loss function, and a more detailed analysis of the splitting size of cross-validation are considered future works.
The average level of accuracies obtained in this study is lower than that obtained in the DIAS ground-based ionosonde technique, which gives the average score of ~0.91 [
11]. For the F2, F1 and E layers, our classification provides accuracies that are comparable with the ARTIST technique [
15]. The DIAS technique resolves the 2D image segmentation problem for 19,000 ionograms. Hence, this very deep learning technique is based on a large number (at least 17 million) of training parameters. While it is interesting to test whether a 1D version of DIAS will benefit the classification of RO profiles, the number of parameters may be too large to be trained for our small dataset. Moreover, it should be noted that the space-borne RO technique of height profiles is essentially different from the ground-based technique of ionograms. The former one might give a global covering of the whole globe, while the latter one is restricted by an observation site, which is strictly localized. Hence, direct comparison of classification accuracies may not be appropriate in the present case.
From the confusion matrices shown in
Figure 4, it can be observed that the number of misclassifications for most of classes is relatively small. In the future, it will be interesting to consider the possibility of applying the present technique for the classification of ionospheric height profiles acquired from the ground-based ionosondes, in addition to the spaceborne RO data. This technique, with some modification, might be also useful for ionospheric modeling, including global TEC models [
29].