A ROTI-Aided Equatorial Plasma Bubbles Detection Method

: In this study, we present a Rate of Total Electron Content Index (ROTI)-aided equatorial plasma bubbles (EPBs) detection method based on a Global Navigation Satellite System (GNSS) ionospheric Total Electron Content (TEC). This technique seeks the EPBs occurrence time according to the ROTI values and then extracts the detrended ionospheric TEC series, which include EPBs signals using a low-order, partial polynomial ﬁtting strategy. The EPBs over the Hong Kong area during the year of 2014 were detected using this technique. The results show that the temporal distribution and occurrence of EPBs over the Hong Kong area are consistent with that of previous reports, and most of the TEC depletion error is smaller than 1.5 TECU (average is 0.63 TECU), suggesting that the detection method is feasible and highly accurate. Furthermore, this technique can extract the TEC depletion series more effectively, especially for those with a long duration, compared to previous method.


Introduction
Equatorial plasma bubbles (EPBs) are very frequent ionospheric events associated with the depletion in plasma density, which are widely believed to be generated by the Rayleigh-Taylor (R-T) instability in the F-layer [1,2]. The EPBs can cause the scintillations and huge ionospheric gradients of radio signals, leading to severe effects on navigation and communication systems [3]. As a result, the study on climatology and the forecasting of EPBs has remained an open question among the scientific community during the last two decades [4][5][6][7][8][9][10][11][12][13][14][15][16]. Many measurements, such as in situ satellite data [4,5], ionosonde [6], and airglow imaging [9,16] are applied to investigate the EPBs characteristics. In recent years, the Global Navigation Satellite System (GNSS) ionospheric Total Electron Content (TEC) has been an important tool for the study of EPBs due to its high spatial-temporal resolution [10][11][12][13][14][15].
To study the EPBs using GNSS TEC data, the first step is to extract the EPBs signals in the TEC series. One common method is using a polynomial fitting with a specific window length to fit the TEC series; then the detrended TEC series that might include EPBs signals are obtained by subtracting the fitted TEC series from the observed TEC series [11][12][13][14][15]. However, the detrended TEC series achieved by this method not only contain depletion but also enhancement, which is similar to the feature of the Medium-Scale Travelling Ionospheric Disturbances (MSTIDs). This exerts an adverse effect on distinguishing EPBs and MSTIDs. Magdaleno et al. [17] presented an Ionospheric Bubble Seeker (IBS) method to detect the EPBs. The IBS method applies the TEC variation and its variance to acquire the initial and end times of the EPBs; then, a linear fitting strategy is used to calculate the background TEC values during the period of EPBs. This method can obtain a good profile of the EPBs compared to the common method. However, the assumption that the profile of the background TEC is a straight line during the period of EPBs is inaccurate, especially for a long duration.
High-precision EPBs signals can help to understand the fine characteristics of EPBs and evaluate their effect on GNSS positioning accurately. To further improve the accuracy of the extracted EPBs signals, in this study, we attempted to employ the so-called Rate of TEC Index (ROTI) to seek the occurrence times of the EPBs and then use the low-order, polynomial fitting to calculate the background TEC values during the period of EPBs. In the following sections, we present the ROTI-aided EPBs detection method in detail and test it using GNSS data in Hong Kong.

Data
GNSS observations are the main data employed in this study. We downloaded the 30ssampling GNSS observation files of station HKKT (receiver type is LEICA GRX1200+GNSS) in 2014 from the public website of Hong Kong's satellite reference network [18].

Methods
The previous IBS method applied the variance of TEC variation to acquire the initial and end times of the EPBs [17]. According to the definition, the ROTI could be seen as the square root of the variance of TEC variation [19]. However, ROTI values are more widely used in monitoring ionospheric irregularities [11,20,21]. So, the "by-product" of the method in this study, namely the ROTI, is also meaningful.
The definition of ROTI is the standard deviation of the Rate of TEC (ROT) [19]. ROT at epoch t i can be computed using the following equation: where M(θ) is the slant factor depending on satellite elevation (θ) suggested by Sanz et al. [22]. Ionospheric TEC time series can be computed using the GNSS dualfrequency, geometry-free, combined carrier phase observations for each pair of satellite and receiver [23]. The ROTI over 5 min intervals of ROT data can be calculated as follows: where the notation · is the averaging operation. The unit of ROT and ROTI is TECU per min, where TECU is the unit of TEC (1 TECU = 10 16 /m 2 ). The ROTI series are then employed to seek the occurrence time of EPBs. Two thresholds, namely "up threshold" and "down threshold", are adopted. For a period of TEC series (T 1 , T 2 ), if the ROTI values are larger than the down threshold and the maximum values are larger than the up threshold, then this period could have an EPB event. Then, we use the low-order polynomial (third order in this case) to fit the TEC data during (T 1 − 7.5 min, T 1 ) and (T 2 , T 2 + 7.5 min). That is to say, the TEC data during the event (T 1 , T 2 ) are not used (see Figure 1). We call this strategy the partial fitting, as opposed to the general fitting method. Lastly, the detrended ionospheric TEC series, including possible EPB signals, are extracted by subtracting the fitted TEC series from the observed TEC series. If the TEC depletion maximum (called "depth") and duration of the EPB candidate meet the defined thresholds then it can be confirmed as an EPB event or vice versa. Table 1 lists the corresponding parameters of the EPBs detection method in this study. The minimums of depth and duration are set to 4 TECU and 10 min, respectively. According to the experimental results, the value of up ROTI could be set to 0.1-0.3 TECU/min (0.2 TECU/min here). The value of down ROTI is very important as it is related to the duration. Figure 2 presents the percent distribution of the ROTI values calculated for station HKKT during 2014. Considering that EPB generally occurs at night, only ROTI values from 6:00-18:00 LT are counted. As shown in Figure 2, the vast majority of the ROTI values are smaller than 0.04 TECU/min. Here, the value of down ROTI is set to 0.04 TECU/min. Furthermore, to distinguish from other events which lead not only to TEC depletion but also TEC enhancement (such as MSTIDs), the maximum of detrended TEC should be less than one-third of the depth. for station HKKT during 2014. Considering that EPB generally occurs at night, only ROTI values from 6:00-18:00 LT are counted. As shown in Figure 2, the vast majority of the ROTI values are smaller than 0.04 TECU/min. Here, the value of down ROTI is set to 0.04 TECU/min. Furthermore, to distinguish from other events which lead not only to TEC depletion but also TEC enhancement (such as MSTIDs), the maximum of detrended TEC should be less than one-third of the depth.     for station HKKT during 2014. Considering that EPB generally occurs at night, only ROTI values from 6:00-18:00 LT are counted. As shown in Figure 2, the vast majority of the ROTI values are smaller than 0.04 TECU/min. Here, the value of down ROTI is set to 0.04 TECU/min. Furthermore, to distinguish from other events which lead not only to TEC depletion but also TEC enhancement (such as MSTIDs), the maximum of detrended TEC should be less than one-third of the depth.    (1) and (2). The blue dotted line and red dotted line mark the thresholds of the up ROTI (0.2 TECU/min) and down ROTI (0.04 TECU/min), respectively. As shown in Figure 3b, the ROTI series during the period of (T 1 , T 2 ) meet the condition of the up and down thresholds (the maximum is 0.49 TECU/min). So, ionospheric irregularity during this period can be marked as an EPB candidate. The partially fitted TEC series, using a third-order polynomial, are also presented in panel (a) for comparison (green line). The red line in panel (c) is the detrended TEC series obtained by subtracting the fitted from the observed TEC series. As shown in Figure 3c, the depths (5.5 TECU) and duration (24 min) of the EPB candidates meet the defined thresholds. So, it can be confirmed that there is an EPB event in the ionosphere during the period of (T 1 , T 2 ).
Remote Sens. 2021, 13, x FOR PEER REVIEW 4 of 9 Figure 3 presents an example of the EPBs detection method in this study using station HKKT and satellite R04 on 28 February 2014. The blue line in panel (a) is the observed TEC series. The black line in panel (b) is the corresponding ROTI series calculated using Equations (1) and (2). The blue dotted line and red dotted line mark the thresholds of the up ROTI (0.2 TECU/min) and down ROTI (0.04 TECU/min), respectively. As shown in Figure 3b, the ROTI series during the period of ( , ) meet the condition of the up and down thresholds (the maximum is 0.49 TECU/min). So, ionospheric irregularity during this period can be marked as an EPB candidate. The partially fitted TEC series, using a third-order polynomial, are also presented in panel (a) for comparison (green line). The red line in panel (c) is the detrended TEC series obtained by subtracting the fitted from the observed TEC series. As shown in Figure 3c, the depths (5.5 TECU) and duration (24 min) of the EPB candidates meet the defined thresholds. So, it can be confirmed that there is an EPB event in the ionosphere during the period of ( , ).   Figure 4a, EPBs generally occur after 20:00 local time and are observed throughout the nighttime. As shown in Figure 4b, the occurrence of EPBs during the spring and autumn equinoxes is significantly greater than their occurrence during the summer and winter solstices. The observed temporal distribution and occurrence of EPBs over the Hong Kong area are also reported in previous studies [12,13], suggesting that the detection method is feasible. As shown in Figure 4a, part of the values of the depletion series is positive. Theoretically, the values of the depletion series should be negative. So, the maximum positive value (MPV) of each depletion series can serve as a reference index to assess the accuracy of the detection result. Figure 5 presents the percent distribution of   Figure 4a, EPBs generally occur after 20:00 local time and are observed throughout the nighttime. As shown in Figure 4b, the occurrence of EPBs during the spring and autumn equinoxes is significantly greater than their occurrence during the summer and winter solstices. The observed temporal distribution and occurrence of EPBs over the Hong Kong area are also reported in previous studies [12,13], suggesting that the detection method is feasible. As shown in Figure 4a, part of the values of the depletion series is positive. Theoretically, the values of the depletion series should be negative. So, the maximum positive value (MPV) of each depletion series can serve as a reference index to assess the accuracy of the detection result. Figure 5 presents the percent distribution of   As indicated in the introduction, an important difference between the EPBs detection method in this study and that of the IBS method is the strategy to derive the detrended TEC series. Compared to the IBS method, the third-order polynomial fitting but not linear fitting (namely first-order polynomial) is employed in this study. Figure 6 presents the detrended TEC series using first-order and third-order partial polynomial fittings for the same example in Figure 3. The green lines and red lines represent the first-order and thirdorder fitting results, respectively. The duration of this EPB event is about 24 min, which can serve as a short-duration case. As shown in Figure 6, the TEC depth difference is insignificant between the different fitting orders. Figure 7 presents a case of a long-duration EPB as observed by station HKKT and satellite G09. As shown in Figure 7, the duration    As indicated in the introduction, an important difference between the EPBs detection method in this study and that of the IBS method is the strategy to derive the detrended TEC series. Compared to the IBS method, the third-order polynomial fitting but not linear fitting (namely first-order polynomial) is employed in this study. Figure 6 presents the detrended TEC series using first-order and third-order partial polynomial fittings for the same example in Figure 3. The green lines and red lines represent the first-order and thirdorder fitting results, respectively. The duration of this EPB event is about 24 min, which can serve as a short-duration case. As shown in Figure 6, the TEC depth difference is insignificant between the different fitting orders. Figure 7 presents a case of a long-duration EPB as observed by station HKKT and satellite G09. As shown in Figure 7, the duration As indicated in the introduction, an important difference between the EPBs detection method in this study and that of the IBS method is the strategy to derive the detrended TEC series. Compared to the IBS method, the third-order polynomial fitting but not linear fitting (namely first-order polynomial) is employed in this study. Figure 6 presents the detrended TEC series using first-order and third-order partial polynomial fittings for the same example in Figure 3. The green lines and red lines represent the first-order and third-order fitting results, respectively. The duration of this EPB event is about 24 min, which can serve as a short-duration case. As shown in Figure 6, the TEC depth difference is insignificant between the different fitting orders. Figure 7 presents a case of a long-duration EPB as observed by station HKKT and satellite G09. As shown in Figure 7, the duration of this EPB event is about 2 h. The fitted TEC series by the third-order polynomial are smoother and more fitting than those by the first-order polynomial (seen Figure 6a). The depth difference is obvious with a magnitude of about 7 TECU between the first-order polynomial (depth 30 TECU) and the third-order polynomial (depth 37 TECU).

Results
Remote Sens. 2021, 13, x FOR PEER REVIEW 6 of 9 of this EPB event is about 2 h. The fitted TEC series by the third-order polynomial are smoother and more fitting than those by the first-order polynomial (seen Figure 6a). The depth difference is obvious with a magnitude of about 7 TECU between the first-order polynomial (depth 30 TECU) and the third-order polynomial (depth 37 TECU).  Detrended TEC (order=3) Detrended TEC (order=1) Figure 6. Same example in Figure 3 with first-order and third-order partial polynomial fitting. (a) Observed total electron content (TEC, blue line) and partially fitted TEC using first-order (green line) and third-order polynomial (red line); (b) detrended TEC by subtracting between the observed and partially fitted TEC (green line for first order and red line for third order).
Remote Sens. 2021, 13, x FOR PEER REVIEW 6 of 9 of this EPB event is about 2 h. The fitted TEC series by the third-order polynomial are smoother and more fitting than those by the first-order polynomial (seen Figure 6a). The depth difference is obvious with a magnitude of about 7 TECU between the first-order polynomial (depth 30 TECU) and the third-order polynomial (depth 37 TECU).

Discussion
In order to adequately extract the TEC depletion by EPBs, the key is to adequately fit the background TEC data during the period of the plasma bubble event. Compared to the IBS method, the ROTI-aided method in this study can more adequately fit the background TEC data during the period of the plasma bubble event. For the IBS method, a partial linear fitting is employed to fit the background TEC during the period of EPB event [17]. In this method, the TEC data during the period of EPB event are not used in the fitting process. Therefore, the effect of EPB on the fitting is removed. However, it is well-known that the variation of the background TEC is not linear. For a short duration, the linear fitting error is small (see Figure 6); the linear fitting error increases with duration (see Figure 7). Therefore, the low-order polynomial fitting is a more suitable fit for the background TEC data. Then, the TEC depletion by the EPB can be adequately extracted (see Figure 7).
It is worth noting that the magnitude of the polynomial order for the ROTI-aided method should not be too high, as in the case of the sixth order. For the EPB with a long duration, the partially fitted TEC may be unreliable with a high-order polynomial. Figure 8 presents the above two cases with a third-order and a sixth-order partial polynomial fitting strategy, respectively. As shown in the figure, for short-duration EPB, both third-order and sixth-order partial polynomial fitting results are suitable. However, for a long-duration EPB, the sixth-order partial polynomial fitting result is not reliable. detrended TEC by subtracting the partially fitted TEC from the observed TEC (green line for first order and red line for third order).

Discussion
In order to adequately extract the TEC depletion by EPBs, the key is to adequately fit the background TEC data during the period of the plasma bubble event. Compared to the IBS method, the ROTI-aided method in this study can more adequately fit the background TEC data during the period of the plasma bubble event. For the IBS method, a partial linear fitting is employed to fit the background TEC during the period of EPB event [17]. In this method, the TEC data during the period of EPB event are not used in the fitting process. Therefore, the effect of EPB on the fitting is removed. However, it is well-known that the variation of the background TEC is not linear. For a short duration, the linear fitting error is small (see Figure 6); the linear fitting error increases with duration (see Figure 7). Therefore, the low-order polynomial fitting is a more suitable fit for the background TEC data. Then, the TEC depletion by the EPB can be adequately extracted (see Figure 7).
It is worth noting that the magnitude of the polynomial order for the ROTI-aided method should not be too high, as in the case of the sixth order. For the EPB with a long duration, the partially fitted TEC may be unreliable with a high-order polynomial. Figure  8 presents the above two cases with a third-order and a sixth-order partial polynomial fitting strategy, respectively. As shown in the figure, for short-duration EPB, both thirdorder and sixth-order partial polynomial fitting results are suitable. However, for a longduration EPB, the sixth-order partial polynomial fitting result is not reliable. As shown in Equation (1), the effect of satellite elevation is calibrated when computing the ROTI value as suggested by Sanz et al. [22]. So, the magnitude of the ROTI values is less than that without the satellite elevation correction [19]. However, it will not affect the extraction of EPB signals because we only used the ROTI values to seek the EPB occurrence time. Certainly, the parameters of the ROTI-aided EPBs detection method should be adjusted when using ROTI values without satellite elevation correction. In addition, the ROTI values were different for various sampling rates [24]. The effect on the detection method was similar to that of the satellite elevation correction. As shown in Equation (1), the effect of satellite elevation is calibrated when computing the ROTI value as suggested by Sanz et al. [22]. So, the magnitude of the ROTI values is less than that without the satellite elevation correction [19]. However, it will not affect the extraction of EPB signals because we only used the ROTI values to seek the EPB occurrence time. Certainly, the parameters of the ROTI-aided EPBs detection method should be adjusted when using ROTI values without satellite elevation correction. In addition, the ROTI values were different for various sampling rates [24]. The effect on the detection method was similar to that of the satellite elevation correction.

Conclusions and Perspectives
In this study, we present a ROTI-aided method to detect the EPBs based on the GNSS TEC data. This method develops the previous IBS method using a low-order polynomial Remote Sens. 2021, 13, 4356 8 of 9 fitting instead of a linear fitting to extract the detrended ionospheric TEC during the period of EPB. The EPBs over the Hong Kong area during 2014 were detected using this technique with a third-order polynomial fitting strategy. The main results and conclusions are listed as follows:

1.
EPBs generally occur after 20:00 local time and are observed throughout the nighttime and the occurrence during the spring and autumn equinoxes is significantly greater than during the summer and winter solstices, which is also reported in previous studies. In addition, most of the TEC depletion error is smaller than 1.5 TECU and the mean error is 0.63 TECU. These findings suggest that the detection method is feasible and highly accurate.

2.
For short-duration EPB, the fitted TEC series are similar between the third-order polynomial and linear fitting; for long-duration EPB, the fitted TEC series via the thirdorder polynomial are smoother and more fitting than those via a linear fitting. This suggests that the third-order polynomial fitting is more suitable to fit the background TEC data, and thus adequately extract the TEC depletion by the EPB.
Overall, this method can acquire high-precision EPBs signals which can help to understand the fine characteristics of EPBs and evaluate their effect on GNSS positioning accurately. Here, only the EPBs over the Hong Kong area are tested. In the future, more cases should be investigated to further check this method.