Automated Detection and Characterization of Wave Structures Obtained From GNSS Measurements

A simple technique designed to automatically identify and characterize wave structures from total electron content (TEC) data obtained from Global Navigation Satellite System (GNSS) satellites and multiple receiver stations is presented. We used 11 years of GNSS data (one complete solar cycle) to detect and characterize the traveling ionospheric disturbances (TIDs) (or wave structures) over high latitude (65°N ± 5°, 147°W ± 5°—Poker Flat, AK) and middle latitude (40°N ± 5°, 117°W ± 5°—Mt. Moses, NV) regions. The algorithm is capable of automatically detecting basic wave parameters such as wave period, horizontal phase velocity, amplitude, wave propagation direction, and wavelength. The designed algorithm can be applicable in the following areas: (a) global GNSS‐TEC (b) climatology of waves (c) input into innovative machine learning algorithms, such as ABCGAN (e.g., Valentic, 2023, https://doi.org/10.5281/zenodo.7747377 designed to characterize background ionospheric plasma and predict wave‐like high/low‐frequency perturbations) (d) anomaly detection for real‐time scenarios. Furthermore, we apply the developed wave detect and characterization technique to investigate three rockets launched on 26 and 28 January 2015. We examine the distribution of different scales of TIDs, and how they varied from high (Poker Flat) to middle (Mt. Moses) latitudes. Lastly, we show that the wave structures at the high‐latitude regions of Poker Flat are substantially affected by auroral processes and those from the middle‐latitude regions of Mt. Moses are impacted by AGWs coupling from below.

AGWs-TIDs can be generated in auroral regions from Joule heating caused by enhanced geomagnetic storm condition.Azeem et al. (2017) and Vadas et al. (2012) show how severe meteorological events such as thunderstorms and tornadoes can generate AGWs-TIDs which can appear as concentric rings at higher altitudes if the intervening neutral winds are small.The signatures of concentric AGWs-TIDs have also been observed in the ionosphere by Azeem et al. (2017Azeem et al. ( , 2015) ) and Nishioka et al. (2013).Several authors have shown that deep convection in the troposphere is one of the primary drivers of concentric AGWs-TIDs that can propagate upward into the mesopause and ionospheric regions (Alexander, 1996;Holton & Alexander, 1999;Jonah et al., 2016;Lane et al., 2001;Vadas & Fritts, 2009).Recent theoretical studies have also shown that some AGWs from deep convection can also propagate efficiently into the thermosphere (Kherani et al., 2009;Vadas, 2007;Vadas & Crowley, 2010;Vadas et al., 2014).Anthropogenic events or man-made activities; nuclear detonations as well as earthquakes-tsunamis and other transient events (such as rocket lunching, earthquakes, and tornado events) are potential sources of smaller scale AGWs-TIDs in Earth's atmosphere (Artru et al., 2005;Hines, 1967;Jonah et al., 2021;Komjathy et al., 2016;Liu et al., 2006;Lognonńe et al., 2006;Rolland et al., 2010).Furthermore, the investigation of AGWs-TIDs using dense Global Navigation Satellite System (GNSS) receiver network, which can be used measure the total electron content (TEC) along satellite-receiver paths with high temporal resolution, has been proven to be a practicable means of observing and studying different kinds of wave structures and their associated irregularities in the upper atmosphere (Coster et al., 2017;Jonah et al., 2018Jonah et al., , 2021;;Komjathy et al., 2016;Zhang et al., 2017).Giving the densely distributed GNSS-receivers and easy access of the measurement globally, GNSS observations have been widely used to study the spatial variations and other characteristics of TIDs and other wave structures.However, most studies determine wave properties by manually examining a series of keogram/hodogram figures or by inspecting time varying density/TEC maps.This approach is not only prone to human errors but also liable to information overload and mental bias and is often subject to ambiguity or decision fatigue.As a result, we have seen evidence of discrepancies in observational results and consequent interpretations of TIDs or wave structure analysis in the literature.For example, the climatology study by Kotake et al. (2006) over Southern California showed daytime MSTIDs with the velocity, period, and wavelength of 80-180 m/s, 20-35 min, and 100-250 km, respectively.Whereas the study by Ding et al. (2011) over central China reported daytime MSTIDs with velocity and period of 100-400 m/s, 20-60 min.On the other hand, the study by Figueredo et al. (2018) over South America reported daytime MSTIDs with velocity, period, and wavelength of 323 ± 81 m/s 24 ± 5 min, 452 ± 107 km, respectively.While the case study of daytime MSTIDs over similar locations in South America by Jonah et al. (2016) shows that MSTIDs are characterized by 255-389 km wavelength, 122-189 m/s, and 20-55 min, periodicity.These discrepancies in the wave parameters do not only result from different sources or locations of the waves but could also be due to the approaches used in identifying and classifying the wave parameters by the many authors.Furthermore, Belehaki et al. (2020) worked on Warning and Mitigation Technologies for Traveling Ionospheric Disturbances , where they exactly identified and tracked TIDs with different scales.Recent paper by Borries et al. (2023) also worked on a new TID index (ATID), which is based on correlation analysis with upper atmospheric drivers and approach for LSTID detection.These studies dealt with identification and prediction of wave but not classifications of their parameters.The present study focuses on the determination and classification of wave parameters which is important for wave source identification and characterization.In this study, we describe the simple approach of automatically detecting wave parameters, including the direction of propagation, wavelength, velocity, and period of the wave.We also present some case studies of small-scale TIDs from transient events (rocket launching events).Lastly, we present the climatology and statistics analysis of MSTIDs and LSTIDs from 11-year data computation.The climatology and statistics studies as well as different analyses of various driving sources will be detailed in future work.

Method of Processing GNSS-TEC From Multiple RINEX Observations
The technique for processing TEC from multiple dual-frequency Global Navigation Satellite System (GNSS) receivers has been developed by numerous authors (e.g., Komjathy, 2005;Mannucci et al., 1998;Rideout & Coster, 2006).The MIT Automated Processing of GPS (MAPGPS) developed by Rideout and Coster (2006) is employed for the current study because of the efficient method of calculating the receiver biases.The method of calculating the individual receiver biases is essential for accurate TEC estimation because one of the largest sources of error in estimating TEC from GPS data is the determination of the unknown receiver biases.Thus, the bias determination procedure for the MAPGPS has recently been updated by using newly developed weighted 10.1029/2023EA003183 3 of 17 linear least squares of independent differences as described in Vierimen et al. (2016).TEC results from MAPGPS algorithm have been validated/used by many researchers in numerous publications (e.g., Coster et al., 2017;Jonah et al., 2018Jonah et al., , 2020;;Zhang et al., 2019, etc.).There are several steps involved in the automated processing of MAPGPS TEC: 1. Download and read all versions of RINEX, and other data formats.
The RINEX files consist of the observables such as the L1 and L2 frequencies as well as the satellite position which are used to estimate the TEC (Mannucci et al., 1998).
2. The ionospheric delay (∆ρ) on GPS signals can be described in terms of TEC, inversely proportional to a square of frequency.
3. Calculate STEC by integrating electron density along the signal path from each receiver to all satellites using a combination of processed f1 and f2 frequencies on L1 and L2 pseudorange and phase data.Screen for and correct loss of lock in the carrier-phase observables (Blewitt, 1990), and use the carrier-phase to smooth the pseudorange values: where f 1 and f 2 are the GPS L1 and L2 frequencies, Φ is the carrier phase measurement.
4. Compute the receiver and satellite biases, b, c, respectively.The step-by-step method of calculating instrument bias is explained in Vierimen et al. (2016).Biases represent the additional delay between measured L1 and L2 GPS signals at the receiver due to both satellite and receiver hardware.The TEC with instrument bias considered is called the absolute TEC and is given as: where b is the receiver bias, c is the satellite bias, and  ℑ is the measurement noise.The measurement is scaled to TEC units, that is, 10 16 m −2 .(Rideout & Coster, 2006;Vierimen et al., 2016) 5. Compute the mapping function, which is the multiplicative factor used to convert line of sight TEC to vertical TEC (VTEC).The TEC obtained here is referred to as VTEC.The mapping function used is given as: where and h = 335km and RE is the radius of the earth given as 6,378.1 km.
Lastly, a 30° cutoff elevation for ground-satellite ray paths is used to eliminate data close to the horizon.For the further description of the above TEC estimation processing, see Rideout and Coster (2006) and Vierimen et al. (2016).

Differential-TEC Estimation Method
After obtaining the STEC and VTEC in the processes explained in Section 2.1, we developed an algorithm to compute the differential-TEC (∆TEC) and an automated procedure to detect and identify seven different wave structure parameters.In the analysis discussed below, we first described the estimated method of obtaining the ∆TEC, followed by the wave parameter detection procedure.
In this study, we focused on GNSS TEC data around Poker Flat (within +5°) to extract different wave structures.More than 40 GNSS receiver stations found within the radar perimeter (65°N ± 5° latitude, 147°W ± 5° longitudes) are used in our analysis (see Figure 1).Each receiver station is used to collect measurements from all satellites in view and are passed through our designed wave processing procedure discussed below.First, determine the kind of filtering software best suited for our analysis.There are three main types of filtering applications that are commonly used by researchers.(a) The polynomial filtering method (Jonah et al., 2016); (b) the low-pass Savitzky-Golay filter (e.g., Savitzky & Golay, 1964;Zhang et al., 2019) and band-pass filter approach (e.g., Jonah et al., 2020).We disregard the polynomial filter as it is well known that it becomes unstable with high degrees.
After examining the two other most used method that is, the bandpass and the Savitzky and Golay (1964) low pass filter (refer to the supporting material), we found the bandpass method most suitable for this application.Detail analysis of this experiment is discussed in the Supporting Information S1.

Procedures for Wave Identifications
There are three main stages: (stage 1) the computation of TEC: This can be obtained from Madrigal website or processed using alternative methods as discussed above.Stage 2 involves the computation of wave structure 10.1029/2023EA003183 5 of 17 according to the following steps: First, it is important to note that 3 days (centered around the day of interest) are needed for the wave parameter computation.This helps to limit/eliminate edge effect problems which usual occur in filtering analysis.We reject signal-to-noise ratio that are not up to 1.1 or 0.04 TECU.(c) For medium scale waves, we reject periods less than 5 min, which is the Brunt Vaisala peak period for wave propagations (Snively & Pasko, 2003).
Figure 3 illustrates examples of how edge effects and artifacts can occur if gaps and discontinuities in the data are not handled carefully following the above procedure.Figure 3a shows the data gap observed from signal between receiver "ac65" and prn 30 and the when filter analysis is applied without considering the data gaps.This artificially create a large amplitude perturbation that can be mistaken for proper wave structure particularly during irregular events such as geomagnetic storm.Figure 3b demonstrates the result when proper measures are taken, before applying the analysis.Furthermore, Figure 3c shows that the bandpass filter works well both at when smaller window data are used to avoid data gap and when large window (without data gap) is used for satellite with PRN 31 and receiver station "ac53."The blue and cyan curves illustrate the application of whole data and 15 min window data set.Lastly, in Figure 3d, PRN 8 and receiver station "spbg" are used; we show an example of edge effect problem in the wave analysis.The blue and green curves show the results with and without considering the edge problem.The blue curve at the onset of wave could easily be mistaken for higher amplitude of the wave.This is just filtering error as a result of edge problems.

Estimation of Different Wave Parameters
The third stage (stage 3) describes how our algorithm compute wave parameters (periodicity, horizontal wavelength, wave propagation direction, and wave velocity) as observed one satellite-receiver-pair from GNSS-TEC data.Detail discussion is given below.At least two receiver stations and all available satellite-receiver paired are need.Two receiver stations are particularly needed to calculate the wave velocity.All the prn (satellite) are averaged out to compensate for the satellite motion's effect on the estimated wave velocity.More detailed explanations about this procedure will be discussed later in this section.
1. Period: The wave period is obtained by taking the spectra analysis of estimated wave structure as a function of time using Morlet wavelet transform (Figure 4c).Over a given time window, the wavelet transformation is collapsed to 1D by finding the maximum wave power for that period (Figure 4d).The dominate period and range of dominate periods are found from the peak and full-width half-max of Figure 4d, respectively.2. Wavelength (λ): To calculate the λ, we first calculate the Euclidean distance (EUD) using the latitude (lat) and longitude (lon) information from the satellite PRNs.EUD is calculated from the Earth-Centered, Earth-Fixed (ECEF) coordinates of the PRN location relative to the ECEF coordinates reference point using the Pythagorean theorem.
Note: reference point is the epicenter coordinates of the particular event of interest.
The wavelength is finally determined by taking the spectra analysis of wave structure as a function of estimated EUD using Morlet wavelet transform.

Velocity (vel):
To calculate the velocity, two adjacent stations are required, and two different approaches were used.a. From Figure 4b, the algorithm finds two dominant consecutive peaks (point A and point B) from the two wave structures obtained from each satellite-receiver-pair, measures the distance between points A and B (red curve) and points C and D (blue curve) which are the two dominant peaks and take their average.The average time difference associated to these points are also noted.Lastly, the change in displacement (ΔX), and change in time (ΔT), are used to compute the wavelength (λ), period (τ) and consequently, the velocity (λ/τ) of the waves.This process is implemented for all specific receiver stations and all PRNs in view, then average-out to minimize the error from satellite motion.b.The relationship between frequency ( f ), and wavelength (λ) of wave property (i.e., f λ), where f = 1/τ is applied to estimate the wave speed.A continuous wavelet transform was performed and applied to analyze the wave structure as a function of time and EuD.This technique allowed us to identify the dominant periods and dominant wavelengths in a particular wave structure (as shown in Figure 4d).Then, the λ/τ relation is applied to estimate the speed of the observed waves.The percentage difference between the "first approach" and the "second approach" for measuring wave speed is calculated to be ∼0.56%,which indicate good consistency in both approaches.4. Wave Propagation direction (∅): To determine the orientation of wave propagation (with bearing angle from north), we applied the Morlet wavelet transfer to the wave oscillation in the latitudinal and longitudinal directions.This allows us to estimate wavelength in both in the latitudinal and longitudinal directions (Haralambous & Paul, 2023;Hocking, 2001).Assuming the meridional wavelength is related to the wave observed in the latitudinal wavelength, which is given by  2∕   and the horizontal wavelength is related to the wave observed in the longitudinal wavelength given by  2∕   .Thus, the direction of wave propagation direction can be defined as where k y and k x are waves in the longitudinal and latitudinal directions respectively.
The four-quadrant inverse tangent is used to resolve the phase ambiguity of a tan  (
Furthermore, given that wavelet function output cannot be negative values, we obtained the negative sign of the waves (i.e., westward or southward directions of the waves) by identifying the dominant wavelength that falls on  10.1029/2023EA003183 9 of 17 the southern and western side of the reference station from the wave oscillation and apply the sign to the λ x (latitude) or λ y (longitude).

Specific Case Studies
The frequency of rocket lunches has increased over the years.Recently, SpaceX successfully launched its massive Falcon Heavy rocket for the first time on 6 February 2018.This signals a new era of space exploration and facilitate the study of space and atmospheric science.Many researchers have shown that rocket launch events are potential sources of ionospheric disturbances around the lunch sites (Afraimovich et al., 2002;Bowling et al., 2013;Lin et al., 2017).The blast impact generated as a consequence of these launch events can cause strong sharp changes in pressure, temperature, density and electrical conductivity (Drobzheva et al., 2003).According to Mendillo (1981), the dominant cause of atmospheric perturbations due to rocket exhaust is mainly associated with the variety of chemical reactions that can occur between molecular materials and the neutral and ionized components of the atmosphere.Lin, Chen, Matsumura, Lin, & Kakinami (2017) first reported the rocket induced shock waves and concentric TIDs (CTIDs) using Global Positioning System TEC over California-Pacific region.They suggested that the CTIDs are the manifestation of concentric GWs that were originated from the mesopause region.Here, we use the ground-based GNSS network over Poker Flat to observe the ionospheric responses to three different rocket launch.The GNSS's wide coverage with dense TEC observations provides an opportunity to monitor the onset and evolution of rocket-induced atmospheric waves.In this study, we investigate the consequent effect of the rocket launch by using GNSS-TEC data set to examine the atmospheric changes resulting from these transient events during and after the launches.We analyzed three different rocket launch events from the Poker Flat Rocket Range in Alaska: the NASA 46.009 UE, NASA 46.010 UE both launched on 26 January with apogee of 160 km, and the NASA 49.002 UE launched on 28 January with an apogee of 590 km.Detailed information about these events are indicated in Table 1.We observed wave propagation in the direction of the lunch, relatively small amplitude, short-period ionospheric and velocity in the range of sound waves.The wave characteristics of rocket event that reached 160 km is similar with that of the rocket event that reached 590 km as shown in Figures 6 and 7.
Powerful rocket launching has been demonstrated to cause disturbances in the upper atmosphere (Afraimovich et al., 2002;Mendillo, 1981).Mendillo (1981) observed ionospheric hole phenomenon, which was considered to result from the following three coupled processes: (a) diffusion of an exhaust cloud of highly reactive neutral molecules through a tenuous, multi-constituent atmosphere, (b) chemical reactions between the expanding cloud and the various species (ions and neutral) in the atmosphere, and (c) solar and/or dynamical replenishment processes.While Afraimovich et al. (2002) recorded shock acoustic wave (SAW) with azimuth angle of the wave vector that varies from 30° to 60°, and the SAW phase velocity of 900-1,200 m/s coordinated with the sound velocity at heights of the ionospheric F region peak.Review paper by Karlov et al. (1980), shows that the oscillation period of the ionospheric response, obtained from the Apollo rocket launching mission, varied from 6 to 90 min, and the propagation velocity was in the range from 600 to 1,670 m/s. Figure 5 indicates that there was a minor geomagnetic storm on 26 January 2015, around the time of the NASA 46.009UE and NASA 46.009UE rocket launches.To eliminate the influence of the geomagnetic condition during the launch period, we use the bandpass filter with a period of 5-15 min (Chou et al., 2018).This allows us to extract only wave perturbations generated by the launch activities.Previous studies have shown that the periodicity of wave structures induced by rocket launch, or other anthropogenic or transient events ranges between 6 and 8 min, and the waves structures generated from geomagnetical activities are usually above 60 min periodicity (Hunsucker, 1982;Jonah et al., 2018Jonah et al., , 2020;;Richmond, 1978;Zhang et al., 2019).Therefore, by using the above filter band width we have ultimately eliminated the effect of geomagnetic storm on our data.started appearing after 16 and 14 min of the launch the first and second launch, respectively.Figure 6f shows that the wave traveled in the southeastward direction with average velocity of ∼500 m/s.The wave properties such as velocity, amplitude, and propagation direction on 26 January events are very analogous with the 28 January event except for the time delay.The difference time delay could be as a result of different lunch apogees.The rocket lunches on 26 January events only reached 160 km altitude, while that of 28 January lunch event reach 590 km altitude.

Eleven Years Wave Distribution Parameters Over High and Low Latitude Regions
The processing of GNSS-TEC data to extract and analyze multiple wave parameters are carried out for two locations as shown in Figure 1.One at high latitude (around the auroral region, 65° ± 5° latitudes) and the other at middle (around 40° ± 5° latitudes).About 50 GNSS receivers were identified around each region, and their RINEX data were downloaded and processed.
Figures 8a and 8b show distribution of hourly detected waves over high latitude (Poker Flat region) and low latitude (Mt.Moses region) for MSTID and LSTID bands.We found similar wave distribution of MSTIDs (blue  Ding et al., 2007;Tsugawa et al., 2004).Fuller-Rowell et al. (1996) used modeling analysis to explain that the energy input in the auroral region can heat the thermosphere, drive equatorward wind surges, which can greatly contribute to the seeding of equatorward-traveling LSTIDs.Similarly, Horvath and Lovell (2010) showed that energy input from the auroral region can heat the thermosphere and propel an equatorward wind, providing a driving force for LSTIDs.Ding et al. (2007) used GPS-TEC data to detect LSTIDs linked to the westward auroral electrojet as detected through decreases in the H and X components of the magnetic field.Jonah et al. (2018) revealed that the observed growth in the LSTID was interconnected with this intermittent energy input from the auroral source into the ionospheric system in a clear indication of magnetosphere-ionosphere coupling.Thus, LSTIDs are commonly generated during geomagnetic storms from the period auroral energy input (Jonah et al., 2020;Tsugawa et al., 2003) and large subauroral polarization stream-induced ionospheric flows (Zakharenkova et al., 2016;Zhang et al., 2019).Therefore, the high occurrence of LSTIDs over Poker Flat (high latitudes) region is mostly as a result of auroral activity prominent in the region.On the other hand, Mt.
Moses, a middle latitude location, mountain waves are very prominent because of the land mountains and valleys topography associated to the region.This mountain wave waves are known to generate GWs that can consistently seed MSTIDs (Heale et al., 2020;McLandress et al., 2012).
Figures 9a-9d and 10a-10d show the summaries of distributions of different parameters extracted around high latitude Poker Flat, AK, and middle latitude Mt.Moses, NV.The period, wavelength, and velocity have a more Gaussian distribution shape for MSTIDs than for LSTIDs.This is because MSTIDs can be generated very often by regular sources such as weather convection from lower atmosphere (Azeem et al., 2017(Azeem et al., , 2015;;Jonah et al., 2018Jonah et al., , 2016) ) and other localized sources (e.g., mountain waves).However, the LSTIDs wave parameters (such as velocity and wavelength) for Poker Flat are right shewed, while the LSTIDs wavelength distribution for Mt.Moses is characterized with double Gaussian distribution compared to the wavelength distribution at Poker Flat which formed only single bell shape.This indicates that the LSTIDs at Mt. Moses are affected primarily by strong winds created by mountain waves associated with the region.Furthermore, the wave direction of propagation (shown in the polar plot in Figure 9) of both MSTIDs and LSTIDs over Poker Flat have predominantly southward component which is consist with the effect of prominent LSTIDs auroral source over the region.Moreover, the direction of wave propagation at Mt. Moses (shown in the polar plot in Figure 10) consists of two clear components: the southward and northward.This discrepancy can also be associated to the different source of AGW at Mt. Moses, middle latitude region.The AGW source at Mt. Moses is predominantly the mountain waves and strong winds compared to the AGW source at Poker Flat which result from persistent auroral energy associated to the region.

Limitations of the Automated Wave Parameter Computation Approach
1.The approach used to obtain the wave parameters is based on hourly data.This is done by iterating on each hour with 2 hr before and 2 hr after, which essentially sum-up to a 4-5-hr data set which is equivalent to one satellite arc.This means that the signal processing is applied to each satellite pass for a single satellite to receiver pair.The limitation here is that hours at the edges are completed by day before (for 0 hr) and day after (for the 23rd hour).Sometimes these hours are not available for the either the day before or after and therefore the data set is discarded.This contributed to why the distribution at the edges in Figures 9 and 10 is low compared to other hours.2. Since only 4-5-hr data set is considered in our approach for the signal processing or filter analysis, the wave period above 2 hr may not be captured by our method.That is, a different of computing the wave parameter based on all 24 hr for each may show a slightly different result for the LSTIDs.This approach will be tested in future studies and in the part II of the research which is in progress.3. Bias in the GNSS receiver distribution over Poker Flat could be a contributing factor to the observed propagation direction wave structure.This is an important question which may be difficult to answer in this current   study.However, probably by using other instruments like PFISR and magnetometers over the region, more wave sources can be studied in future endeavors to further illuminating our understanding.This is out of scope of the current paper and has been left to future studies.

Conclusion
In this study, we designed an automated algorithm that is capable of automatically detecting anomalous wave structures and computing basic wave parameters such as wave period, horizontal phase velocity, amplitude, wave propagation direction, and wavelength.We analyzed three rocket lunch events of 2015 as showed in Table 1 and compute different wave parameters to validate our approach.Our analyses demonstrate that the rocket lunch events trigger ionospheric perturbation with wave structures moving with 420-600 m/s in the southeastward direction.We also specify that lunch apogee could have significant impact on the time of wave arrival at the ionospheric altitudes (300 km).The wave detection algorithm was also used to study the morphological and statistical distribution of LSTIDs and MSTIDs from 11 years of GNSS-TEC data over Poker Flat and Mt.Moses.This analysis reveals that MSTIDs over Poker Flat and Mt.Moses have similar occurrence morphology; however, the morphology of LSTIDs over Poker Flat is very different from that of Mt.Moses, with Poker Flat exhibiting higher occurrence rate of LSTIDs compared to Mt. Moses.This is expected because Poker Flat is a location over high latitude where there are intermittent auroral energy inputs which are principal source of LSTIDs (see Jonah et al., 2018;Zhang et al., 2019).The directions of propagation observed in this study are consistent with literatures.The automatic wave detection algorithm described in this paper can be applied to many areas, including

Figure 1 .
Figure 1.A and B represent examples of the Global Navigation Satellite System receiver station distribution within ±5° (for 28 January 2015) over Poker Flat and Mount Moses, respectively.AA and BB represent the zoomed images and station names of A and B. The colored horizontal lines represent lines of constant magnetic latitude over North America.
Figure2illustrates the flowchart of the approach and measures taken to ensure cleaned input data before the zero-phase Butterworth bandpass filtering analysis is applied.These procedures are explained as follows: (a) 5 hr of data centered around the period of interest which make approximately one satellite arc period are collected and arranged in a time sequence.(b) All multiple NaN values (i.e., large data gaps) as a result of data arrangement in time are removed and data set with one or two NaN values are interpolated using nearest neighbor interpolation method.(c) Remove data where NaNs cover over half the 5-hr window.(d) All forms of discontinuity for example, data gaps due to cycle slips, multipaths, or power outages are detected and removed from the analysis.(e) Next, the data is detrended to remove any trend in the data set that can make the filter bias toward normal day-to-day variation of the data.(f) The edge effects are avoided by repeating 10% of the data set to both ends of the original data before applying the filter.This is known as padding technique.(g) The bandpass filter analysis is then carried out and the 10% padding data are removed from the boundaries before we continue with the rest of the process.After the filter analysis is done, the following measures are taken to avoid misinterpreting noise as signal or wave structures.(a) We reject wave amplitudes that are less than 10% of the TEC amplitude.(b)

Figure 2 .
Figure 2. Flowchat illustrating different processing stages of the wave detection algorithm.

Figure 3 .
Figure 3. Wave structure analysis on 19 June 2015.(a) ∆TEC obtained when filter analysis is applied without considering the data gaps (i.e., data gap included in the filter input), (b) ∆TEC obtained when data gap is removed before applying the filter analyze.Blue and red curves represent the ∆TEC and total electron content (TEC), respectively.(c) Blue color shows ∆TEC.The cyan color is derived from bandpass filter with a 15-min window interval.(d) Blue color shows the wave analysis without edge effect treatments. Green color shows wave analysis with edge effect treatments.

Figure 4 .
Figure 4. (a) Wave structures from multiple receiver stations on 2 January 2017.Left panel shows diff-total electron content (TEC) as a function of time to obtain the period.Right panel shows similar wave structure but plotted as a function of distance to obtain the wavelength.(b) The red and blue curves represent wave structures from single PRN, but two different receiver stations used for the computation of phase velocity.(c) Wavelet analysis showing example of dominant periodicity.(d) The collapsed 1D wavelet transformation, used to identify the dominate periods from (c).The vertical red lines give the range of dominate periods for this event.

Figures
Figures 6a-6d represent hourly spatial-temporal variation of differential TEC (DTEC) variations for two rocket launch events on 26 January 2015.A strong disturbance can be clearly observed around the launch hour (panel a) at 9:13 UT.It is still possible to observe similar perturbations around same location (at panel b) which could be from the second rocket launch at 9:46 UT. Figure 6e represents the hodogram with latitude as a function of time.The red and blue vertical dash line indicate the time of first and second rocket launches.The wave perturbations

Figures
Figures 7a-7d also represent hourly temporal-spatial variation of DTEC (∆TEC) variations but for NASA 41.002 UE rocket launch event on 28 January 2015.A strong disturbance can be clearly observed around the launch hour (panel b). Figure 8e represents the hodogram with latitude as a function of time.The dash vertical red line indicates the specific time of the rocket launch.It is possible to see that after about 6-8 min of the launch the wave structures started appearing close to the to the location of the launch.Panel (f) shows that the wave traveled in the southeastward direction with average velocity of ∼500 m/s.

Figure 5 .
Figure 5. Top, middle, and bottom panels represent the Bz, SymH, and the AE indies for days 26-28 of January 2015.Day 26 show a minor geomagnetic storm activity while day 28 is relatively calm.The dash-lines indicate the rocket lunch onset.The 5-15 min bandpass window was effectively used to minimized geomagnetic storm source from our results.

Figure 6 .
Figure 6.(a-b) Represents hourly temporal-spatial variation of differential TEC (∆TEC) variations.The star marker represents the location of the launch.Panels (e and f) represent the hodogram (time/distance diagram) and the polar plot that indicates the direction of wave propagation.

Figure 7 .
Figure 7. Same as Figure 6 but for day 28 January 2015.

Figure 8 .
Figure 8. (a and b) Distribution of hours detected waves as a function of UT hour of the day for Poker Flat, AK, and Mt.Moses, NV, respectively.

Figure 9 .
Figure 9. (a-d) Represent wave period, horizontal wavelength, horizontal phase velocity, and wave propagation direction, respectively for high latitude (Poker Flat) region.The blue and brown colors represent medium-scale traveling ionospheric disturbance (MSTID) and large scale traveling ionospheric disturbance (LSTID), respectively.
(a) apply to other GNSS TEC data from receivers anywhere around the world (b) easily apply in climatology study of wave analyses.(c) Serves as input into innovative machine learning algorithms ABCGAN (e.g., Valentic, 2023) designed to characterize background plasma parameters of the ionosphere and predict wave-like high/low-frequency perturbations during set of conditional drivers.The actual radar data and wave parameters from GNSS data and drivers used to train the ABCGAN has been published as ABCdata (e.g., Reyes et al., 2023).(d) Serves as anomaly detection for real-time scenarios.This material is based on work supported by the Defense Advanced Research Projects Agency (DARPA) under Contract No. HR001121C0026 and the "Atmosphere as a Sensor (Atmosense)" program.Special thanks to MIT Haystack and Bill Rideout for providing the access to TEC processing resources and Madrigal database which facilitate the TEC data processing in this study.

Table 1
Information About Rocket Launches Over Poker Flat on 26 and 28 January 2015