A detailed investigation of low latitude tweek atmospherics observed by the WHU ELF/VLF receiver: I. Automatic detection and analysis method

As a dispersive wave mode produced by lightning strokes, tweek atmospherics provide important hints of lower ionospheric (i.e., D‐region) electron density. Based on data accumulation from the WHU ELF/VLF receiver system, we develop an automatic detection module in terms of the maximum‐entropy‐spectral‐estimation (MESE) method to identify unambiguous instances of low latitude tweeks. We justify the feasibility of our procedure through a detailed analysis of the data observed at the Suizhou Station (31.57°N, 113.32°E) on 17 February 2016. A total of 3961 tweeks were registered by visual inspection; the automatic detection method captured 4342 tweeks, of which 3361 were correct ones, producing a correctness percentage of 77.4% (= 3361/4342) and a false alarm rate of 22.6% (= 981/4342). A Short‐Time Fourier Transformation (STFT) was also applied to trace the power spectral profiles of identified tweeks and to evaluate the tweek propagation distance. It is found that the fitting accuracy of the frequency–time curve and the relative difference of propagation distance between the two methods through the slope and through the intercept can be used to further improve the accuracy of automatic tweek identification. We suggest that our automatic tweek detection and analysis method therefore supplies a valuable means to investigate features of low latitude tweek atmospherics and associated ionospheric parameters comprehensively.


Introduction
The majority of the energy of natural lightning discharge radiates as electromagnetic waves in the extremely-low-frequency (ELF: 0.3-3 kHz) and very-low-frequency range (VLF: 3-30 kHz) (Kumar et al., 2008;Volland et al., 1987), most of which propagate with multiple reflections through the Earth-ionosphere waveguide (EI-WG). During waveguide propagation, dispersive effects occur near the cutoff frequencies for a portion of VLF sferics, which are called 'tweek atmospherics'. As tweeks reflect from the lower ionosphere, ground-based detectors can monitor these electromagnetic emissions. Information obtained regarding the frequency spectrum characteristics and cutoff frequencies of tweeks can subsequently be adopted to evaluate the wave reflection height in the ionospheric D-region and the corresponding electron density (Kumar et al., 1994;Reeve and Rycroft, 1972). In this way, distinct from a number of conventional high-frequency (HF) devices that aim to detect the ionospheric E-and F-region actively, ground-based ELF/VLF receivers provide an economically useful tool to monitor passively the lower ionospheric D-region.
Networks of ground-based ELF/VLF receivers have been deployed over many of the world's regions, including the Atmospheric Weather Electromagnetic System for Observation, Modeling, and Education (AWESOME) (Cohen et al., 2010;Carpenter et al., 2012), the World Wide Lightning Location Network (WWLLN) (Dowden et al., 2002), and the Antarctic-Arctic Radiation-belt (Dynamic) Deposition-VLF Atmospheric Research Konsortium (AARD-DVARK) (Clilverd et al., 2009), to observe natural whistler wave emissions and man-made VLF transmitter signals. Based upon all these ELF/VLF receiver observation campaigns, much progress has been achieved in the past decade, in particular a number of useful applications in the context of geophysical studies of the ionosphere and magnetosphere, such as radio atmospherics, ionospheric remote sensing, and ELF/VLF generation and detection (e.g., Ohya et al., 2008;Kumar et al., 2009;Yusop et al., 2013;Maurya et al., 2012;Singh et al., 2011).
Since July 2015 a new digital low-frequency receiver system has been developed at Wuhan University (WHU) (Chen YP et al., 2016, 2017Wang YP et al., 2019;Yi J et al., 2019) for sensitive reception of low-latitude broadband ELF/VLF radio waves between 1-50 kHz originating from either natural or artificial sources. Figure 1 shows the structure diagram of WHU-ELF/VLF system and Figure 2 exhibits the antenna, the pre-amplifier box, and the chip box of the WHU-ELF/VLF receiver system. As shown in Figure 1, this system distinguishes itself by using the field programmable gate array (FPGA) in the digital receiver and adopting the Universal Serial Bus 2.0 (USB 2.0) in the transmission module. The type of FPGA used in WHU-ELF/VLF receiver system is FPGA (EP3C55F484). The receiver system uses a 16 bit A/D converter chip with 250 kSPS (kilo samples per second) throughput. A standard GPS receiver card (M12M + TIMING ONCORE) is implemented in the system to provide 1-pulse-per-second (1 PPS) timing signals. The phase locked loop modulation circuit insures minimal deviation of the system clock, i.e., < 0.005 Hz, which is shorter than the data sampling rate. Significantly, the timing accuracy of GPS timing signals is ~100 ns, which is much less than the 4 ms timing accuracy of the WHU-ELF/VLF receiver sampling chip. And according to Dowden et al. (2002) and Lay et al. (2004), the timing accuracy of the WWLLN hardware is ~20 ms, corresponding to a sampling frequency of 48 kHz. Therefore, the WHU ELF/VLF receiver system, when a VLF receiver network becomes available, has the potential for locating lightning activities as well as does the WWLLN. Some key parameters of the WHU-ELF/VLF receiver system are given in Table 1. Readers are referred to Chen YP et al. (2016, 2017 for more details.   Using the first experimental results from the WHU ELF/VLF receiver system at the Half Hill Station of Wuhan University (30.54°N, 114.37°E), Chen YP et al. (2016) reported the clear signature of low latitude tweek atmospherics over central China at multiple modes, up to n = 6, and demonstrated a strong correlation of both the reflection height and mean cut-off frequency with mode number. Since January 2016 the Suizhou station of the WHU ELF/VLF receiver system has become operational at a rural village in Suizhou City, China (31.57°N, 113.32°E), whose environment of electromagnetic interference, which is encountered by power-line harmonic radiation, is ideally quiet for the purpose of clear measurements. The receiver operates routinely to measure ELF/VLF emissions in two alternative modes: daytime (08-18 LT) mode (recording the first 10 seconds in each minute) and nighttime (18-08 LT of next day) mode (recording 0-10 s, 15-25 s, 30-40 s, and 45-55 s in each minute). Both modes produce data files that record for 10-second intervals. It is worthwhile to point out that other ELF/VLF receiver systems, such as those operated by Nagoya University (Ohya et al., 2011(Ohya et al., , 2012(Ohya et al., , 2015, record much less data per hour than this system, and at sampling rates much lower than is used in our system. Therefore, the accumulation of high quality data by the WHU ELF/VLF receiver provides a unique opportunity to investigate in detail the spectral characteristics of low-latitude tweek atmospherics and their occurrence patterns over central China, which can be further adopted to explore the features of associated lower ionospheric parameters. Due to the substantially frequent occurrences of tweek atmospherics resulting from lightning strokes, it becomes both labor and time prohibitive to establish a robust tweek database by manual identification. Instead, a reasonable automatic detection method is required to go through all available datasets and identify tweeks automatically and reliably, which is the focus of the present study. The outline of this paper is as follows. In Section 2 we develop a new detection procedure, based upon the maximum-entropy-spectral-estimate (MESE) method, to automatically capture tweek atmospherics. This procedure is applied to the data observed at the Suizhou Station on 17 February 2016, and the automatic detection results are validated with results from careful visual observations. Section 3 presents the construction of the analysis module for identified tweek atmospherics by applying a Short-Time Fourier Transformation (STFT) to the traces of tweek power spectra and by evaluating the reflection height, lower ionospheric electron density, and propagation distance from each source lightning stroke. Section 4 discusses our ELF/VLF system and developed analysis method; Section 5 presents concluding remarks.

Atmospherics
To identify tweek atmospherics events, we focus on one of the two antenna channels (i.e., the North-South channel) to extract wave data obtained by the WHU ELF/VLF receiver. For each 10second dataset, we acquire 2.5×10 6 data points with 250 kSPS. For each dataset, the maximum-entropy-spectral-estimation (MESE) method is adopted to process each group of 512 consecutive data points allocated starting from the first point. The MESE method, as a method of spectral density estimation, chooses the spectrum that corresponds to the most random or the most unpredictable time series whose autocorrelation function agrees with the known values (Burg, 1975). By improving the spectral quality based on the principle of maximum entropy, the MESE method is capable of deriving, with enhanced accuracy, a better estimate of the spectrum than can the Fourier transformation method. In our automatic detection procedure, we choose the function of the Burg's method provided by the MATLAB platform to estimate the autoregressive power spectral density. Figure 3 presents a representative 10-s time-frequency spectrogram with two zoom-in plots recorded at 22:50:00-22:50:10 LT on 17 February 2016 at the Suizhou station, processed and obtained by the MESE method. This method makes it easy to set the noise threshold level for data processing. Visualized in Figure 3a, we can see many tweek activities associated with sferics, which can be seen in the two bottom panels in Figure 3, featured by the 'impulsive' structures that cover almost the entire frequency range.
Due to the high occurrence rate of tweeks, it is impossible to identify manually (i.e., by vision) every tweek in a long-term run. In order to investigate efficiently the statistical properties of tweek atmospherics over central China, we develop a new automatic detection technique, on the basis of the MESE method, to determine tweek events captured by the WHU ELF/VLF receiver.
The major procedures of the automatic detection method are described as follows: (1) Application of the MESE method to the data. This procedure separates every 10-second data file into 4882 equal data segments. Every segment has 512 data points, and is processed by the MESE method to obtain the energy spectrum. The combination of all energy spectra constitutes a matrix [E ij ] for each 10second data file, which can then be adopted to plot the time-frequency spectrogram as shown by Figure 3. The subscripts i and j denote the time segment T i and frequency segment F j respectively.
(2) Registration of useful signals. First, we set a threshold value, thd, dependent on the local magnetic environment at the receiver station. Note that the pre-set number thd is not time dynamic but a constant that is adopted to deal with the inference from background noise. In order to reach a good balance between maintenance of most useful data points and inhibition of background noise interference, after many test runs we chose for purposes of our analysis the value of thd to be 20 dBFS above the averaged background noise intensity, over the frequency range of 1-10 kHz. It is worth noting that the averaged background noise intensity is measured on a weak interference day. Thus, the averaged background noise intensity and thd are both constant in our method. If E ij > thd, we set it as a useful point. If not, we set no signal. Since both tweeks and sferics have similar spectral profiles that cover almost the entire frequency range of 1-10 kHz within a very short time interval, in this step we register both wave activities at the time stamp T i as long as the ratio of identified useful signal points to total points is greater than a pre-set value β p .
(3) Tweek identification. While similar to sferics, tweeks are dis- persive waves with a dispersive 'tail'. We introduce the parameter t d to denote the shortest time duration over which a tweek event can be identified in the dispersive tail of the observed wave frequency spectrum. Given the time stamp T i , if there are useful signals between 1-3 kHz during the time interval between T i and T i +t d , then we identify that interval as including a tweek event. In order to avoid repeated registration of the same tweek, we introduce the parameter Δt to denote the smallest time interval above which two identified tweeks can be regarded as isolated events.
We then further specify that no tweek exists during the time interval between T i and T i +Δt.
As described above, the prescribed parameters β p , t d and Δt can affect more or less the performance of the automatic detection method. To explain these parameters clearly, Figure 4 graphically represents the meaning of parameters β p , t d and Δt with black arrows and dash lines. According to Figure 4, t d indicates the dispersion time duration of a tweek, Δt exhibits the time interval between two closed tweeks and β p displays the frequency range that a tweek covers in 1-10 kHz. Notably, β p is been normalized by the total frequency range (10-1 = 9 kHz).
In order to check the sensitivity to these parameters and then to adopt an appropriate set of parameters for automatic detection, we concentrate our detailed analysis on a 10-min dataset recor-  Table 2 shows the performance results of our automatic detection method corresponding to 9 groups of parameters β p , t d , Δt for the considered time period. Our visual check identifies 144 tweeks in total. In Table 2, the number of tweeks that are found by both automatic detection and visual check is marked as 'correctness'. If tweeks are recognized by visual check and not recognized via automatic detection, they will count as 'missing tweeks'. 'False alarm' tweeks are signals captured by automatic detection while visual check results do not consider them as tweeks. The total number of tweeks refers to tweeks observed visually during the entire time period under consideration. We can see that, while changes in β introduce only small variations in the automatic detection performance, the latter is more sensitive to the values of and . As and decrease, our detection method can capture a larger number of correct tweeks but also a larger number of false alarm tweeks. Therefore, we have to select an appropriate set of the three parameters to balance between high correctness and low false alert rate. After careful consideration, we choose for subsequent analyses the values β p = 2/3, t d = 10.24 ms, and Δt = 6.144 ms, to give the best overall performance. To evaluate carefully the feasibility of our developed automatic detection method, we collect the actual tweeks that cannot be identified by our method for further exploration. Figure 5 displays zoom-in plots of the tweeks occurring in the 10-second period corresponding to Figure 3. In each panel, the red arrows denote the tweeks captured by our automatic detection method and the blue arrows denote the tweeks missed by our method. For the sake of robustness of our developed tweek database for quantitative analysis, we demand that a good tweek event should be strong enough in both the 'hiatus' and dispersive tail portions to stand out from the background noise level. We also require that the dispersive tail of a useful tweek should have a duration longer than 10 ms. It is seen that the missing tweek events are caused either by wave strength hiatus at frequencies of 1-10 kHz and weak wave intensity (e.g., Figure 5d) or by short time intervals between two neighboring events (e.g., Figures 5f and 5h). However, it is worthwhile to note that weak tweeks or ones with an obscure trace are difficult to measure accurately and very likely to provide poor quality data with misleading information that cannot be included. In addition, when one tweek has a very short interval (< 10 ms) to another, the second tweek will affect the res-ult of curve fitting of the first tweek (the curve fitting process will be mentioned in Section 3). And this will lead to a bad correlation coefficient for the first tweek, which means the first tweek can be excluded from the database and has little influence in our analysis result. Under these considerations, we judge that application of our new automatic detection method to the WHU ELF/VLF receiver data during the entire month of February 2016 can provide a robust database of tweek atmospherics for a detailed analysis of its occurrence features at low latitudes over central China.
In order to further evaluate the applicability of our automatic detection method, we apply it to the 1-day period of 17 February 2016, which provides 4-hour data sets for investigation, and we compare the automatic results with those manually analyzed.
Based on Table 3, which gives detailed information of the analysis results of our automatic method on February 17, we evaluate the reliability and validity of the automatic detection procedure. During these four hours, the quantity of actually occurring (visually observed) tweeks in each hour varies from 500 to 1350, compared to which the number of auto-detected tweeks varies between 640 and 1800. This produces a correctness percentage of   well above 74% and a missing percentage that is less than 17%. In total, 4342 signals are captured by the automatic detection method with 981 false alert events and 600 missing tweeks, producing a 77.4% overall percentage of correctness and an overall missing rate of 15.1%.
In general, about a quarter of the automatically-detected tweek events are falsely alerted by the automatic detection method and nearly one-sixth of the confirmed events are missed. However, as

Earth and Planetary Physics
doi: 10.26464/epp2020018 125 discussed above, the missing tweeks with low intensity or short duration time or unclear dispersive features are insufficient to provide reliable information for subsequent analysis. As a consequence, we conclude that our automatic detection method can offer an effective and efficient tool to reasonably identify welldefined tweek events for more detailed investigation of the properties of low-latitude tweek atmospherics.

Analysis Module for Tweek Atmospherics
To analyze the tweeks captured by our automatic detection method, our analysis module has to obtain accurately the dispersive features and the cut-off frequency. Due to the limits of the MESE and Fourier transformation methods, neither can guarantee the required precision in both the time and the frequency domain. Thus, the Short-time Fourier Transform (STFT) method, described by Chen YP et al. (2016), is applied to every 0.1-second time segment in which a tweek occurs. Although the maximum entropy method and the wavelet method can be more suitable for rapid signals like tweeks, both methods would be much more computationally expensive, based on the MATLAB platform, than the STFT method adopted in this study, which can produce results good enough for our automatic detection method. Since we use only the first mode dispersive features, the upper and lower frequency limits of the spectral range are 3.5 kHz and 1.5 kHz respectively.
Given the spectral peak dataset defined as (t i , f i ), the tracing procedure first finds the time stamp t 1 with the strongest spectral intensity and spectral peak (t 1 , f 1 ) at t 1 . Then our analysis module extracts the spectral peaks at each time step. This extraction process does not stop until at least one of the following judgment conditions is satisfied: (1) the occurrence of two consecutive frequency drops that meet the criterion ; (2) the occurrence of one frequency increase that meets the criterion ; (3) the spectral peak (t i , f i ) is 60 dB lower than the maximum spectral peak. Obtained datasets (t i , f i ) are subsequently fitted by a linear curve via the least-mean-square method, and the mean value of the five lowest frequencies is adopted as the cut-off frequency, which can be used for further evaluations of associated ionospheric parameters.
Based on the dispersion characteristics of tweeks in the Earth-Ionosphere Waveguide (EIWG) (Yamashita, 1978;Outsu, 1960), we follow the dispersion relation for tweek atmospherics as where v gn is the group velocity of electromagnetic waves propagating in the waveguide, f is the wave frequency, f cn is the cut-off frequency of nth mode of the tweek, and c is the velocity of light in free space. On account of the planar waveguide model and ray tracing method (Ramachandran et al., 2007;Dowden et al., 2002), the above dispersion relation can be transformed into a time-frequency relation as where d is the tweek propagation distance, t 0 is the actual occurrence time of the tweek, and t is the time that the maximum spectral intensity of the signal is detected by our ELF/VLF receiver. To better illustrate the fitting curve, we use A to represent the independent variable as follow, (3) Subsequently, Equation (2) yields Apparently, Equation (4) introduces a linear fitting curve to reveal the relation between time and transformed frequency (i.e., A). In turn, we can judge the data reliability from the fitting performance.
On the basis of the relation between t, t 0 , f cn and the propagation velocity, there are two different ways to calculate the tweek propagation distance. One is calculated through the slope as where k is the slope of the fitting curve; the other is acquired by the intercept as It is worthwhile to point out that the group velocity of an electromagnetic wave in the free space waveguide is much closer to the light velocity if the wave frequency is much higher than the lower cut-off frequency. Therefore, the arrival time of tweeks, evaluated using Equations (2)-(4), should be more accurate when their high frequency parts are adopted. However, because the frequency range of the pre-filter and pre-amplifier of the WHU-ELF/VLF receiver is confined between 1-50 kHz, the signals above 50 kHz have been unavailable. Additionally, as seen in Figure 3, the receiver captures a number of strong signals of VLF transmitters that occupy the frequencies above 10 kHz, and relatively much weaker emissions at tens of kHz that are hardly usable. Consequently, we use the wave data below 10 kHz to estimate tweek arrival time. Figure 6 depicts an example of tweeks automatically identified from the broadband ELF/VLF signal data on 9 February 2016. The observed wave spectra are then fitted for the first and second modes of tweek atmospherics. The dispersive structure of the tweek power spectrum is obvious in Figure 6a, exhibiting a strong time and frequency dependence. For the two modes of tweeks, the last several white circles along the dispersive track approach the cut-off frequencies at 1.7529 kHz and 3.8391 kHz, respectively. Only the last 20 points are used to fit the dispersion relationship because points that fit the slope part of tweeks (marked by the black arrow) might be discontinous. For the fitting curves in Figure 6b, a larger slope corresponds to wave propagation that has reached a greater distance from the source. Specifically, the blue line with circles presents the linear relation between time and A (defined by Euqation (3)) for the first mode tweek, producing the tweek occurrence time of 1.6063 sec and the propagation distance of 4266.9 km. The red line with circles for the second It is worth pointing out that theoretically there should be little difference in estimating propagation distance from the dispersion of n = 1 and n = 2 tweek mode via Equation (5) and (6) if both modes follow the same propagation channel. In reality, however, due to the fact that the signals of n = 2 generally suffer from more attenuation than those of n = 1 (e.g., Figure 6 indicates that the dispersion trace of n = 2 is much shorter than the trace of n = 1), the differences between estimated distances from the dispersion of the n = 1 and n = 2 can be non-negligible. Therefore, in this study, instead of using the traces of both modes, only the traces of the n = 1 mode of tweeks are used to estimate the propagation distance.
By fitting the datasets (t i , f i ) in terms of Equations (3) and (4), we can obtain the related parameters of tweeks atmospherics. However, as mentioned in Section 2, there is a possibility that some events identified by the automatic detection method can be false, thereby leading to improper curve fittings and unreliable parameter values. Theoretically, for a good tweek event with clear trace information, the tweek propagation distances (i.e., d 1 and d 2 ) evaluated by Equations (5) and (6) should be identical or close. Notably, both propagation distances, i.e., d 1 and d 2 , are evaluated from the first mode tail instead of different mode tails (e.g., higher order harmonics) of one tweek. Therefore, calculations of the correlation coefficient for the fitting and relative difference between d 1 and d 2 provide a useful means to further evaluate the performance of our automatic detection method and to improve the reliability and robustness of obtained tweek database. Figure 7 illustrates the histograms of fitting curve correlation coefficient (CC) and distance relative difference (DRD) for the tweek datasets established in Section 2 using our automatic detection method. The tweek events that fail to provide enough information for propagation distance calculations are excluded, resulting in 3944 events in total for analysis. It is clear that in the entire database a large number (41.23%) of events exhibit distance relative differences > 1, suggesting that they have large uncertainty associated with them. In contrast, 51.70% of the events have distance relative differences ≤ 0.3, which strongly justifies their reliability, together with the fact that 54.34% of the events have correlation coefficients ≥ 0.9. For the tweek database with DRD ≤ 1, we can see that the distance relative difference ≤ 0.3 and the correlation coefficients ≥ 0.9 for the majority of the events (i.e., > 88%), indicates that these events are robust tweeks that can   (3)) for the first mode and second mode tweeks.

Earth and Planetary Physics
doi: 10.26464/epp2020018 127 provide valuable information for detailed analyses. Consequently, we reach a reasonable criterion, that is, that requiring the distance relative difference to be ≤ 0.3 and the correlation coefficients to be ≥ 0.9 between the two different methods given by Equations (5) and (6) should guarantee the reliability of an automatically identified tweek event.

Discussion
While our study has adopted the model described by Equations (1)-(6) to estimate the tweek propagation distances and then obtain estimates of source lightning locations, it is important to compare our model results with direct observations of lightning activity, e.g., data from the WWLLN. Our companion paper (Yi J et al., 2020) has performed such an analysis for the period of this study, February 2016. It concludes that our ELF/VLF receiver at the Suizhou station captures a majority of tweek events with propagation distances < 4000 km, but does not detect as many lightning events at greater distances from Suizhou (i.e., > 4000 km) as are found in WWLLN data. The reasons for the difference may can be many-fold and will be further analyzed in our following investigation. But it should be pointed out that the WWLLN is a networked world wide lightning detection system while in this study we only use one station for tweeks detection.
The WHU ELF/VLF receiver implements a triangle antenna instead of a circle for the sake of its measurement stability and easiness of deployment. The two orthogonal magnetic loop antennas therefore record the wave signals from both the N-S and the E-W directions; thus the available datasets make it possible to estimate the direction of arrival (DoA) of tweeks. Although the present study uses only the E-W channel data to develop the automatic detec-tion and analysis method, goals of our future work include using the two-channel data to evaluate the DoA profile and to obtain improved descriptions of tweek properties, associated ionospheric parameters and causative lightning location.
We also note that the present study utilizes only nighttime (i.e., 22-02 LT) data for analysis. The major reason is that the number of daytime tweeks is much lower than that of nighttime tweeks, as shown in Figure 1 of the companion paper (Yi J et al., 2020). Such a day-night asymmetry of the occurrence pattern of tweeks can be well explained by the fact that the diurnal variation of the lower ionosphere is predominantly controlled by the abundance of sunlight.
In addition, while Figure 6 illustrates an example of tweeks that have both the first and second mode tails for analysis, this study adopts only the data of first mode tails of tweeks to evaluate the cut-off frequency and propagation distance, since not all the observed tweeks have the second mode profile. We leave the analysis of higher order harmonics of tweek emissions as a subject of our future investigation.

Concluding Remarks
Based on data from the WHU ELF/VLF receiver system, we have developed an automatic detection module that utilizes the maximum-entropy-spectral-estimation (MESE) method to identify reliable low latitude tweek events. The feasibility of our procedure has been well justified through a detailed analysis of the data observed at the Suizhou Station on 17 February 2016. Corresponding to 3961 tweeks in total registered by careful vision check, the automatic detection method captures 4342 tweeks, including . Histograms of fitting curve correlation coefficient (CC) and distance relative difference (DRD) for the tweek datasets established in Section 2 using our automatic detection method.

128
Earth and Planetary Physics doi: 10.26464/epp2020018 centage of 77.4% and a missing rate of 15.1%. The short-time Fourier Transformation (STFT) has been further applied to trace the power spectral profile of identified tweeks and to construct a tweek analysis module.
To exclude less reliable events, we have implemented two different ways to calculate the tweek propagation distance (i.e., through the slope, Equation (5); through the intercept, Equation (6)), for which the correlation coefficients of the fit to the frequency-time profile and the relative difference of propagation distance between the two methods are adopted to further improve the reliability of the analysis. By defining reasonable criteria, that is, that the distance relative difference be ≤ 0.3 and the correlation coefficients be ≥ 0.9, our established automatic detection and analysis method supplies a valuable means to investigate comprehensively the occurrence characteristics of low latitude tweek atmospherics. Its underlying sources and associated ionospheric parameters are evaluated in detail in a companion paper of Yi J et al. (2020).