Statistical Characteristics of Nighttime Medium‐Scale Traveling Ionospheric Disturbances From 10‐Years of Airglow Observation by the Machine Learning Method

For the first time, we used the machine learning method to analyze the statistical occurrence and propagation characteristics of nighttime medium‐scale traveling ionospheric disturbances (MSTIDs) from October 2011 to December 2021 observed by the all‐sky airglow imager deployed at Xinglong (40.4°N, 117.6°E, 30.5° MLAT), China. We developed a program code using the algorithms to identify and extract the propagation and morphological features of MSTIDs in 630 nm airglow images automatically. The classification model and detection model have accuracies of 96.9% and 70%–85%, respectively. We identified 611 MSTID events from 749,888 airglow images, and obtained the following statistical results: (a) the MSTIDs occurrence peaked at 2200–2300 local time in summer and 2300–2400 in winter; (b) the annual average of horizontal wavelength and velocity are 160–311 km and 98–133 m/s, respectively; (c) among 611 events, 589 MSTIDs propagated southwestward. Fifteen events are northeastward and all of them are periodic MSTIDs, most of which occurred between April and August; (d) the annual trend of relative intensity perturbation (%) shows a negative correlation with the horizontal phase speed; (e) horizontal wavelengths of MSTIDs are independent of the solar activity. Further analyses found those southwestward propagating MSTIDs are consistent with the Es‐Perkins coupling theory, while those non‐southwestward ones could be related to the atmospheric gravity waves and other possible sources. The northeastward events exhibit morphological and seasonal characteristics, which cannot be explained by the Perkins instability, more simultaneous observations (GPS‐TEC, OH airglow, etc.) are required to reveal the mechanism behind these characteristics.

satellite system (GNSS) and the global positioning system (GPS) receivers, Yu et al. (2016) found the occurrence of nighttime MSTIDs is consistent with mid-latitude spread F (MSF). Compared with the GNSS/GPS-TEC, an all-sky airglow imager (ASAI) can record the two-dimensional structure of nighttime MSTID through optical observation with higher spatial resolution . Therefore, airglow images are often used to investigate morphological features, such as wavelength, speed, direction, etc. of nighttime MSTIDs (Garcia et al., 2000;Taylor et al., 1998). According to OI 630 nm airglow images from 2006 to 2012 in the Central Pacific and South American sectors, Duly et al. (2013) presented the occurrence rate of MSTID in the Southern Hemisphere peaks at the solstice months (June and December) and has an inverse relationship with the F10.7 cm solar flux. Huang et al. (2016) utilized the airglow images and GPS data to statistically analyze the MSTID over central China during 2013China during -2015 showing the northwest-southeast alignment and southwestward propagation. Huang et al. (2016) combined analyses of the airglow images and GPS data to investigate the characteristic parameters of nighttime MSTIDs over central China during 2013China during -2015 showing the northwest-southeast alignment and southwestward propagation. Also based on an all-sky imager, Xu et al. (2021) reported the morphological features of nighttime MSTIDs during 2012-2014 over Xinglong (40.4°N,117.6°E,30.5° MLAT), China, and obtained consistent results with Huang et al. (2016). In the opposite hemisphere, Candido et al. (2008) first obtained the optically statistical features of nighttime MSTIDs over Brazil and found an antidependence of MSTID on solar activity.
The unique northwestern-southeastern (northeastern-southwestern) alignment and southwestern (northwestern) propagation of nighttime MSTIDs in the Northern (Southern) Hemisphere have been widely reported (e.g., Martinis et al., 2019;Rajesh et al., 2016;Saito et al., 1998). The Perkins instability (Hamza, 1999;Perkins, 1973) is a possible candidate for explaining the wavefront and propagating preferences of nighttime MSTIDs (Garcia et al., 2000;Kelley & Fukao, 1991). Numerical simulated studies conducted by Yokoyama et al. (2008) first successfully reproduced this northwestern-southeastern aligned wavefront of nighttime MSTIDs in Japan. However, the very slow growth rate and the northeastward wave vectors of nighttime MSTIDs derived from the Perkins instability in the Northern Hemisphere are not consistent with the observations (Garcia et al., 2000). Investigators (e.g., Haldoupis et al., 2003;Kelley & Fukao, 1991) suggested that other mechanisms are also involved in the initiation or generation of nighttime MSTIDs. It is now generally accepted that coupling between E and F regions (Cosgrove et al., 2004;Kelley et al., 2003;Otsuka et al., 2004;Yokoyama et al., 2009) and polarization electrical field (Kelley & Fukao, 1991; are supposed to be the alternative mechanisms to accelerate the development of Perkins instability. Statistical investigations are also carried out to explore the potential mechanisms of nighttime MSTIDs. Based on 942 nights of observation at Arecibo from 2002 to 2007 years, Martinis et al. (2010) found the occurrence of nighttime MSTIDs peaked at both solstices and explained the two peaks with E layer/F layer and inter-hemispheric couplings. Yu et al. (2016) compared the TEC retrieved from GNSS at two stations in China during 2001-2012, finding a strong connection between MSTIDs and MSF. They suggested that nighttime MSTIDs could have the same generation mechanism as the mid-latitude MSF. Paulino et al. (2016) studied 98 events from September 2000 to November 2010 and found other important mechanisms besides the Perkins instability should be also involved in the generation of nighttime MSTIDs.
Most of the previous statistical researches, especially the airglow data, were based on manual statistics, which are laborious and may introduce subjective errors. Compared with the time-consuming manual statistic, automatic detection by program is more effective and can avoid subjective errors. To identify MSTID from TEC data, Cheng et al. (2021) developed an autonomous algorithm utilizing both the three-dimensional Fast Fourier Transform (3D FFT) and support vector machine methods. But 3D FFT requires a time series of MSTID-containing images and cannot identify a single airglow image. It is difficult to build a system based on prior knowledge to recognize MSTIDs from the complex environment of airglow images automatically. Machine learning is thus imperative for extracting effective information from massive airglow images (Lai et al., 2019). The ASAI deployed at Xinglong, China has accumulated 749,888 airglow images by taking pictures each night since 2011. We developed an automatic identification program code based on the algorithms of convolutional neural network (CNN) and faster regional convolutional neural network to detect nighttime MSTIDs in airglow images and extract the characteristic parameters of MSTIDs for statistical analysis. We statistically analyzed the morphological features of nighttime MSTIDs over Xinglong, including the occurrence rate, horizontal speed, propagation direction, and wavelength based on the airglow data near a solar cycle (from October 2011 to December 2021). The solar activity dependences of the annual occurrence rate, relative intensity perturbation (%), and horizontal wavelength were also investigated. We further discussed the possible reasons for these results. To our best knowledge, this is the first comprehensively and systematically statistical analysis of long-term continuous optical observation (630 nm) of nighttime MSTIDs by using a machine learning method. This paper is organized as follows: The data processing, including observation, learning of classification model and detection model, is introduced with details in Section 2. The statistical results of MSTIDs, related discussions, and comparisons with previous studies are presented in Section 3. Summaries and conclusions are in Section 4.

Data
All airglow images used in this study were captured by the ASAI located at Xinglong (40.4°N, 117.6°E; 30.5° MLAT), China from October 2011 to December 2021. The ASAI consists of a Mamiya 24 mm/f4.0 fisheye lens with a 180° field of view (FOV), a filter to improve the OI emission at 630 nm, a charge coupled device with 1,024 × 1,024 pixels, and other supporting devices. During the runtime, three different exposure times of 5, 3, and 1.5 min were used in different year.
Four steps of image processing are shown in Figure 1. First, we categorized raw airglow images by using a classification model to distinguish the images on clear nights from those on unclear nights (overexposure, overcast, twilight, etc). The main purpose is to filter banded clouds and light bands out to reduce interference to later MSTID detection. After a preprocessing operation, including azimuth correction, star-removing, and deleting background noise by subtracting the mean value of 1 hr, images of clear nights were further projected to the geographical coordinate. The details of projection method of airglow images can be found in the article by Garcia et al. (1997). The size of a projected image is 1,600 × 1,600 pixels, corresponding to the actual sky of 1,600 × 1,600 km 2 . In the third step, the detection model was applied to identify MSTID wavefronts in the projected airglow images. In the last step, the propagation parameters, including wavelength, horizontal speed, and propagation direction were calculated based on the automatic detection code.

Classification Mode
In this part, we trained a classification mode based on CNN to select images of clear nights from massive original airglow images. We adapt the network architecture of the CNN algorithm (Lai et al., 2019). Previous statistical study on morphological features of atmospheric gravity waves (AGWs) at the height of mesopause (∼87 km) observed by the Xinglong ASAI has proved effective and reliable for our code. Our CNN consists of input layer, two convolution sets, flatten layer and full connection layer in turn. Both convolution sets are composed of convolutional layer, max pooling layer and dropout layer, but with different parameters. Images in the training set are loaded by the input layer and then processed by different layers. The output result of the CNN is five percentage values, normalized by the softmax function, representing the confidence degree of the loaded airglow images in the five categories (starring, brighten, etc. shown in Figure 2) of the loaded airglow images. The loss of CNN is measured by the cross-entropy loss function (see more details in Lai et al. (2019)). In the training process, the learnable parameters were optimized to minimize the loss until the loss-epoch curves became flat. After 150 epochs of training, a classification model with high prediction accuracy (96.9%) and low loss was generated (shown in Figure 3). We only need to distinguish the starring class, which is ideal for observing MSTID, from the other four classes in practice, hence the misjudgments between the other four classes can be ignored.
1909 original airglow images from the observations in 2012 were manually classified into five categories: starring, brightening, moonlight, light band and cloudy (shown in Figure 2). The number (percentage) of selected image in each category is 562(29.4%), 70 (3.7%), 517 (27.1%), 160 (8.4%), and 600 (31.4%), respectively. Among the manually classified images, 1,554 were used for training sets and 355 for validation sets. Both sets were quadrupled by rotating the images for 90°, 180°, and 270°. An area out of the square of 512 × 512 pixels at the center of an image was cut off when training and classifying, in order to avoid the interference of possible artificial light at the edge of FOV.
The feature maps of five categories were visualized to reveal the criteria of classification (shown in the second line of Figure 2). The brightness in the feature maps indicates the scoring weights in classification. The 90° symmetry was caused by image rotation. The airglow images containing structures similar or partially similar to a feature map will score higher in the corresponding classification and are sorted into that category. Comparing the feature maps with the corresponding classified image, we confirmed that the classification model focused on the correct features for classification.

MSTID Detection
The detection model for locating MSTID wavefronts was generated through the training of Faster RCNN based on a data set consisting of 336 representative real projected airglow images and 250 simulated images. The Milky Ways (if any) and MSTID wavefronts were labeled manually in the total 586 images, respectively.
Considering the labeled MSTID wavefronts are similar, we simulated some images containing wavefronts to prevent overfitting (shown in Figure 4a). A wavefront is defined by the function: where "fix" is a rounding function. d and p are random integers in the ranges from 2 to 6 and 20 to 200, respectively. n is an array of natural numbers from 1 to b. b is rounded at random in range of [20,40], determining the width of the simulated wave surface. When generating a wave surface, array x is rounded continuously in the interval [c, 1600-c], where c is a random integer in range of [50,100]. In different simulated images, the gray scale of generated wave surface varies randomly from 0.5 to 1.5 times of the average background.
The training of MSTID detection model was initiated with the faster_rcnn_ resnet101 _voc07 pre-trained model (Nealwu & Pkulzc, 2016). We obtained the MSTID detection model after one hundred thousand rounds of training. The detection model can locate on MSTID wavefronts and the Milky Way in  a projected airglow image precisely (shown in Figure 4b). The precision is measured with intersection over union (IoU): where GT and DR denote the ground truth and detection results, respectively. For most of images, the IoU of our detection model reaches 70%-85% (shown in Figures 4c and 4d), which is much higher than the general qualification value of 50% (Rosebrock, 2016). Even in a complex case of two wavefronts, the model can still accurately locate the two targets separately.

Result and Discussion of MSTID Detection
Through the classification model and detection model, MSTID wavefronts can be identified in airglow images automatically. The classification model was applied to all data from 2011 to 2021 to sort out the airglow images of clear nights. And then the detection model was utilized to locate MSTIDs in the images taken in clear nights with rectangular boxes. The region in the rectangle was binarized to enhance the edge of an MSTID wavefront. Then the wavefront was fitted to a line to extract parameters, since the MSTID wavefront can be viewed as an arc with a large radius of curvature. Regarding objects with a confidence greater than 90% as MSTID, we detected 611 MSTID events from October 2011 to December 2021.

Occurrence Rate
The occurrence rate is defined as the ratio of MSTID duration to observable duration in all clear nights. The observable duration is the sum of the duration of photos classified as clear night. The MSTID occurrence rate was analyzed with respects of years, seasons, and nocturnal hours. Figure 5 shows the monthly and yearly occurrence rates of MSTIDs from 2011 to 2021. The annual occurrence rates were low in 2014 and high in 2019, indicating an inverse correlation with solar activity. A major peak in the summer months (black bars) and a second peak in the winter months (gray bars) can also be seen. The negative correlation of the seasonal occurrence rate with the annual solar activity is obvious in spring and winter, but not significant in summer. The occurrence rate of autumn varied greatly in different years and shows no correlation with solar activity may due to the low occurrences of clear nights and MSTID events in autumn. Figure 6 illustrates the monthly average occurrence rate of MSTIDs, duration time of clear night and MSTID during 2012-2021. The monthly occurrence rate has a minimum in March and increases significantly from May to August. The seasonal occurrence rate in spring (March-April), summer (May-August), autumn (September-October), and winter (November-February) reach 3.5%, 16%, 4.1%, and 5.1%, respectively. The maximum monthly occurrence rate (29.4%) appears in summer, the next is in winter. Figure 7 illustrates the hourly occurrence rate of nighttime MSTIDs. The occurrence rate has a peak from 2200 to 2400 local time (LT). The peaked hour varies in different seasons: 2200-2300 (LT) in summer and 2300-2400 (LT) in spring and winter.
Our statistics in near a solar cycle (2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020)(2021) show that the annual occurrence rate had a maximum (minimum) value in the minimum (maximum) year of solar activity, confirming the anticorrelation between the MSTID occurrence and solar activity reported by the previous studies (e.g., Candido et al., 2008;Narayanan et al., 2014;Pimenta et al., 2008). The opposite tendency could be explained by the Perkins instability. According to Perkins' theory (Perkins, 1973), the growth rate of Perkins instability is inversely proportional to the collision frequency of neutral-ions. And empirical models indicate the collision frequency of neutral-ions increases with solar F10.7 cm solar flux index (Hedin, 1991). Hence the annual occurrence rate is decreasing (increasing) at ascending (descending) phase of solar activity (Hazeyama et al., 2022). . The linear growth rate of Perkins instability, which is positively correlated with the occurrence rate, can be also used to explain the seasonal vibration (Bowman, 1992;Kotake et al., 2007). The linear growth rate of Perkins instability is inversely proportional to the neutral density in the upper atmosphere. The neutral density trends to minima at solstices and maxima at equinoxes (Martinis et al., 2010), leading to a higher (lower) linear growth in summer and winter. In addition, the difference between summer and winter may be related to the E and F region coupling, which was introduced to explain the formation of nighttime MSTIDs Otsuka et al., 2008). Therefore, the Figure 6. Monthly average medium-scale traveling ionospheric disturbance occurrence rate, duration, and observable time. peak E-region process in summer (Arras et al., 2008) may be the reason for the major peak of occurrence rate in summer. Since the E-region process is at its peak in winter of the Southern Hemisphere, the minor peak in winter can be explained by the inter-hemispheric coupling (Otsuka et al., 2004).
The peaked hourly occurrence rate of those MSTIDs is between 2100 and 0100 LT. This observed trend resembles those in Kotake et al. (2007) We notice that except the similar trend, the value of occurrence is lower than that in the previous literature. It is because that when calculating the clear sky duration, we counted the corresponding duration of each image including the isolated images. However, in order to avoid false positive, we applied a strict threshold when determining the MSTID-containing images so that several pictures at the beginning and end of each event were not included.

Horizontal Propagation Direction and Phase Speed
The horizontal phase speed and direction of MSTIDs were counted based on the average phase speed and direction of each single event. Here 611 MSTID events were counted in Figure 9 by season. Propagation directions of southwestward MSTID events are classified into 5° bins (shown in Figure 8). The propagation directions are distributed from 200° to 260° and peak at 235°-240°. Five hundred eighty-nine MSTID events propagated southwestward and 15 MSTID events propagated northeastward. 12 of 15 northeastward events were captured in April and August. All northwestward events are periodic MSTIDs. Two northwestward and five southeastward MSTIDs were observed. Annual average propagation directions varied between 231° and 228° in azimuth (clockwise from north), without discernible solar activity dependence. The major southwest ward propagation and a minor northeastward peak were consistent with the previous literature Tsuchiya et al., 2018).
The monthly average phase speed varies from 102 to 126 m/s, without an obvious seasonal difference. The annual horizontal phase speed ranges from 98 to 133 m/s (shown in Table 1), showing no significant solar activity dependence.
The southwestward propagation of MSTIDs in the North Hemisphere is widely reported (Martinis et al., 2010;Miller et al., 1997;Shiokawa et al., 2003). The characteristic of southwestward propagation suggests the electrodynamical forces make a key effect on the generation of MSTIDs (Garcia et al., 2000;Kelley & Fukao, 1991;Miller et al., 1997). Kelley and Makela (2001) proposed a dipole polarization electric field model and invoked polarization of the Perkins structures in the direction parallel to their long axis, which produce a northwestern electric field as the cause of the southwestward propagation. An alternative explanation is independent of electric field. Based on numerical simulation, Yokoyama et al. (2009) attributed the southwestward propagation of TIDs in both E and F region to the southward component of neutral wind in the E s -layer.
It is notable that all of 15 northeastward MSTID events are periodic MSTIDs with multiple waveforms, 12 of which occurred between April and August. For the cause of northeastward MSTID, there are two possibilities. One possible situation is that MSTIDs propagating southwestward reversed their propagation directions (Shiokawa  Shiokawa et al. (2008) reported an MSTID event that turning the propagation direction from southwest to northeast, and suggested the reverse was caused by northeastward thermospheric neutral wind or the simultaneous F layer height decrease. According to the model simulation of Kelley and Makela (2001), the change of neutral wind direction from southwest to northeast would lead to the reverse of MSTID propagation direction. Wu et al. (2021) observed a northeastward MSTID event at Fuke (19.5°N, 109.1°E) and suggested that the reverse of MSTID was one possible reason of this unusual propagation; they attributed the reversion to the enhancement of poleward neutral wind or the northeastward movement of background plasma. Another possibility is that these MSTIDs were presumably triggered by AGW, which can enhance the Perkins instability naturally, leading to periodic MSTIDs. When approaching the F layer, AGWs can disturb the plasma density periodically and then this perturbation is amplified by the Perkins instability to form MSTIDs (Kelley & Fukao, 1991;Rathi et al., 2021). The second evidence for AGW is most of these MSTID events occurred between April and August, which is consistent with the seasonal propagation trend of AGWs captured by the OH band ASAI at Xinglong station. The northeastward propagation trend of AGWs over Xinglong is only spotted in summer . The filtering effect of westerly winds prevailing at the critical height (20-70 km) in summer is the reason for this seasonal difference. Statistical study of OH and OI nightglow images captured at Japan also reported the northeastward propagation direction in summer (Ejiri et al., 2003). Third, AGWs due to the deep convections which have large scales are believed to generate disturbances in the ionosphere, leading to MSTIDs (Vadas, 2007). The northeast propagating MSTID mainly occurs between April and August, which is the period of frequent convective weather. Consequently, AGWs should be another possible explanation for the northeastward MSTIDs. Katamzi-Joseph et al. (2022) recently reported the eastward MSTIDs based on observation over South Africa and suggested that besides the gravity waves, the sporadic E and neutral are also candidates for the mechanism behind the unusual propagations. More simultaneous ionograms and OH airglow data are needed before we draw conclusion of the mechanism.

Relative Intensity
The relative intensity I r of an MSTID wavefront is defined as: where I max and I min indicate the maximum and minimum gray values of an MSTID wavefront, respectively. I avg denotes the average gray values of the background in the local region located by the detection model. The anticorrelation of the relative intensity with the solar activity reported by Otsuka et al. (2021) and Tsuchiya et al. (2018) is not found in our statistical result. Instead, the annual relative intensity (shown in Figure 10a) seems negatively correlated with the phase speed (shown in Figure 10b). Since observation based on the OI emission of 630 nm mainly occurs in the range of 250 ± 5 km altitude, more observations are required to discover whether the negative correlation between the relative intensity and horizontal phase speed is caused by accident or other factors (e.g., the development of economy lightens the environment of the dark night).

Horizontal Wavelength
Horizontal wavelength was measured for events with multiple wavefronts in FOV. No obvious relationship with solar activity is found in the annual wavelength trend (shown in Figure 11). Seasonal average wavelengths in spring, summer, autumn, and winter are 282, 246, 239, and 265 km, respectively, showing no obvious seasonal change.

Periodic and Single Dark Band MSTIDs
According to the number of wavefronts in the FOV, MSTIDs can be categorized into two kinds: single dark band MSTID (SDBM) with one wavefront in the FOV and periodic MSTID (PMSTID) with multiple wavefronts in the FOV (e.g., Figueiredo et al., 2018). Table 2 (Figueiredo et al., 2018) and Southern (Amorim et al., 2011;Pimenta et al., 2008) Hemispheres. Propagation of PMSTIDs is mainly in the southwest and northeast directions. The horizontal phase speed of two kinds of MSTID is comparable to the values in the report of Sau et al. (2018). The SDBMs are believed to be created by E-F coupling that seeded the Perkins instability (Rathi et al., 2021), so the southwestward propagation of SDBMs can be explained by the reason we discussed in Section 3.2. Regarding PMSTIDs, the northeastward events are possibly due to AGWs propagated to F layer. E-F coupling generates both SDBMs and PMSTIDs (Rathi et al., 2021) so the southwestward PMSTIDs may be partly caused by E-F couplings, and partly triggered by AGWs.

Summaries and Conclusions
We developed an MSTID automatic extraction program code, including a classification model based on CNN and an MSTID detection model based on the Fast RCNN algorithm. MSTID-containing airglow images were simulated to enlarge the data set for training. The classification model has an accuracy of 96.9% and the detection model can mark the MSTID waveforms with accurate localization rate above 70%. With the program code, we processed 749,888 airglow images captured by an all-sky imager deployed at Xinglong, China near a solar cycle (from October 2011 to December 2021) and detected 611 MSTID events. The characteristic parameters of these events are extracted for comprehensive statistical analysis. The mainly propagation direction, speed, relative intensity, and horizontal wavelength are comparable to manual statistics during 2012-2014 at Xinglong station   (shown in Table 3), which further proves the reliability of our auto detection. The results are summarized as follow:  1. The annual occurrence rate of MSTID is negatively correlated with solar activity, reaching a peak (valley) value in 2019 (2014), the year of solar activity minimum (maximum). The seasonal occurrence rate has a major peak in summer and a minor peak in winter, which can be explained by the E-and F-region coupling and inter-hemispheric couplings, respectively. The negative correlation between the annual change of summer occurrence rate and solar activity is not obvious. The peaked values of the hourly occurrence rate in 2019 and 2014 occurred at 2200-2300 LT in summer, and 2300-2400 LT in winter, respectively. 2. Among 611 MSTID events, 589, 15, 5, and 2 events propagated to the southwest, northeast, northwest, and southeast, respectively. The southwestward events have a wavefront azimuth of 200°-260° with respect with the geographic meridian, with the highest proportion between 235° and 240°. All 15 northeastward events are periodic MSTIDs, in which 12 events occurred between April and August. The northeastward MSTIDs are possibly caused by the reversed southwestward MSTIDs or the AGW propagated in the thermospheric height. The annual and monthly average phase speeds for southwestward and northeastward MSTIDs vary between 98-133 and 102-126 m/s, respectively. Neither the propagation direction nor the phase speed shows the dependence of solar activity. In addition, there is no obvious difference in the phase speed for all seasons. 3. The annual trend of relative intensity shows negatively correlated with the phase speed. 4. The horizontal wavelength does not show an obvious relationship with either season or solar activity.
5. The ratio of periodic MSTIDs to the total events reaches a lowest value of 18.9% in 2014 and a highest value of 45.2% in 2018. Most SDBMs propagate to the southwest. In 22 non-southwest ward events, there are only two SDBMs but 20 periodic MSTIDs.
The morphological and seasonal characteristics of northeastward MSTIDs are reported for the first time, which needs more investigation (e.g., GPS-TEC, OH airglow) to explain. Our investigation generated a database containing all characteristics of nighttime MSTTIDs from October 2011 to December 2021, which will benefit further research in this field.