Using GNSS Phase Observation Residuals and Wavelet Analysis to Detect Earthquakes

With the progress in technology, Global Navigation Satellite System (GNSS) applications expand in the field of global precise positioning with very short latency. Moreover, recent developments in the Galileo and BDS constellations have been shown to improve positioning performance, which is crucial for accurate and fast measurements, including GNSS seismology (Fang et al. 2021).

Initially, GNSS technology was used to determine plate velocities (Teunissen and Montenbruck 2017) and then to determine terrain deformations (DeLoach 1989). Nowadays, GNSS is used for the determination of both short-term and long-term point displacements. Short-term vibrations, such as those occurring during earthquakes, require high-frequency GNSS data sampling, above 1 Hz, because less frequent data may not record all changes in the geophysical process (Smalley 2009). Similar conclusions were reached by Smalley (2009), who analyzed an earthquake with Mw = 7.0 and a 5-Hz frequency. A good compromise between series usefulness and a reasonable amount of data is a sampling rate of 10 Hz (Avallone et al. 2011). The GNSS observations with such a high sampling are called high-rate GNSS with broad applications in GNSS seismology.

High-rate GNSS measurements can complement classical seismic measuring techniques and may be applied to early warning of earthquakes and tsunamis (Shan et al. 2019). Initially, high-rate GNSS earthquake-related measurements were conducted using short GNSS baselines by determining significant point displacements (Hirahara et al. 1995; Melgar et al. 2013) and with instantaneous positioning (Nikolaidis et al. 2001). Recent earthquake studies are mainly conducted with Precise Point Positioning (PPP, Zumberge et al. 1997) because PPP does not require any additional reference stations. Nowadays, there are a plethora of GNSS receivers and sensor networks dedicated to seismological measurements. Among them, one can distinguish the Japanese network called GPS Earth Observation Network (GEONET) (Sagiya 2004) or the Network of the Americas (NOTA) (Dittmann et al. 2020). Other examples are the network of permanently operating reference stations (CORS) in Turkey and Northern Cyprus (Snay and Soler 2008) and the Central Apennine Geodetic Network (CaGeoNet) of the National Institute of Geophysics and Volcanology (INGV) of Italy (Anzidei et al. 2005). These stations measure GNSS signals and other parameters at high frequencies so that earthquake oscillations can be recorded in satellite observations.

The co-seismic ground vibrations are not always evident in the data recorded by the high-rate GNSS receivers mainly due to the poor signal-to-noise ratio. In general, with the increase of the epicentral distance, the amplitudes of the vibration become smaller, and thus, the possibility to detect the event with the high-rate GNSS decreases. However, strong earthquakes may impact even stations located far away from the earthquake epicenter. Such a case took place for the NUTS station located 1000 km away from the earthquake epicenter (MW 9.0) on December 26, 2004, which moved the station by 15 mm (Kouba 2005).

Currently, the PPP technique is the most commonly used technique for earthquake studies (Kudlacik et al. 2019; Tu et al. 2013) due to the limitations of the double-difference method (Li et al. 2018), where the position changes of the reference station may occur and preclude the use of double-difference method. PPP can accurately and efficiently determine the magnitude (Fang et al. 2013; Geng et al. 2016) and epicenters of earthquakes (Fang et al. 2013). In addition, it has been proven that, over short intervals, PPP solutions are comparable to double-difference solutions (Xu et al. 2019) and can detect movements with millimeter accuracy in the horizontal component and centimeter accuracy in the vertical component (Xu et al. 2013). Other methods used for GNSS seismology are as follows: Temporal Point Positioning (TPP) method which uses single differences to accurately determine displacements from a single GNSS receiver (Li et al. 2013) or variometric approach and its variants (Colosimo, Crespi, and Mazzoni 2011; Li et al. 2015) which operate on a different principle – instead of deriving the absolute displacement values, the high-rate GNSS station velocities are determined based on which earthquakes are detected. The high-rate GNSS is mainly applied to strong natural earthquakes; however, some research studies show the potential of the high-rate GNSS for detecting small amplitudes of natural earthquakes (Saunders et al. 2016).

Although the high-rate GNSS allows for deriving high-accuracy results, it suffers from the noise that must be reduced to achieve reliable waveforms. Hence, various methods for denoising are applied to the high-rate GNSS, such as sidereal filtering (Atkins and Ziebart 2016; Wang et al. 2018), Butterworth filtering (Kudlacik et al. 2019; Kudłacik et al. 2021), or wavelets denoising (Kaloop, Cemal O. Yigit, et al. 2020; Li, Xu, and Yi 2017). The Butterworth filter (Butterworth 1930) is one of the most popular anti-aliasing filters (Smith 2021) mainly due to a flat frequency response that does not distort the time series (Kudlacik et al. 2019). The time series used in GNSS seismology is usually short with a length of 20 minutes (Li et al. 2013) or even 5 minutes (Kudlacik et al. 2019), and due to its non-stationary nature, the Fourier transform is not preferred. Therefore, this work employs the wavelet transform.

In the GNSS seismology, the displacements estimated based on the high-rate GNSS are commonly analyzed (Ge et al. 2000; Kudłacik et al. 2021; Nikolaidis et al. 2001; Paziewski et al. 2020; Riquelme et al. 2016). The high-rate GNSS data might also be utilized to assess earthquake mechanisms as complementary to seismic sensors (Hung et al. 2017; Hung and Rau 2013; Xu et al. 2013).

The second field of applications of the high-rate GNSS is structural monitoring of infrastructure that is prone to deformations (Kaloop, Cemal Ozer Yigit, et al. 2020; Shen et al. 2019). For monitoring purposes, GNSS observation residuals can be used as a side result of static positioning. It is proven that the data collected by the high-rate GNSS receiver located on the bridge towers confirm the usefulness of GNSS residuals as an additional source of information about point displacements (Larocca et al. 2006). Importantly, the measurement results were also found to be influenced by the satellite elevation angle relative to the predicted movement of the object.

The GNSS residuals are pseudo-distance corrections, which are obtained as an intermediate stage in the position determination process and quantify the inconsistency between mathematical assumptions and observations (Teunissen and Montenbruck 2017). All errors and background models used in the processing influence the values of the residuals. To better distinguish the influence of individual errors, conversions of the residuals to the frequency domain are performed (Larocca et al. 2006). Therefore, it is sometimes better to analyze the residuals of the observations rather than the coordinates (Suzuki 2019). Static positioning may not contain information about the short-term deformation or vibrations of the points because of the strong constraints of the estimated GNSS station position, which will not change when the shock reaches the GNSS receiver. In contrast, the long-term deformation can be identified with tectonic movements. For the static GNSS positioning, observation residual values contain information about high-frequency station coordinate vibrations due to the earthquakes, which are not visible in the coordinate time series obtained in the static mode. The preliminary residual analysis may contain information which, due to various factors, will be invisible in the coordinate time series, especially in the static positioning mode.

The purpose of this study is to evaluate the usefulness of GNSS observation residuals for the ground vibration detection. The main scope of the research is to develop the methodology for analyzing the time series of GNSS observation residuals in the context of earthquake detection. The proposed methodology provides the process of acquisition, processing, and evaluation of the observation residuals for the detection of point vibrations and the moment of the earthquake. The methodology is tested on a previously well-elaborated earthquake.

In this study, the earthquake from October 26, 2016, in Italy is scrutinized (USGS ID: us1000725y, https://earthquake.usgs.gov/earthquakes/eventpage/us1000725y/executive). The earthquake was previously analyzed by many authors (e.g., Avallone et al. 2011; De Guidi et al. 2017; Kudlacik et al. 2019; Zhong et al. 2018). Table 1 provides a detailed parameter description of the earthquake parameters which were recorded by 4 sets of sensors located in pairs, including a GNSS receiver and a seismograph (Table 2). Figure 1 shows the position of ACCU, AMAT, GUMA, and MTER stations relative to the epicenter. The stations are located at different azimuths from the epicenter, which means that the stations can record the signal in different ways. Additionally, the GNSS receivers are located close to the epicenter, and thus, it is possible to eliminate the influence of different propagation speeds of P and S seismic waves.

Figure 1.

The distance between GNSS receivers and the seismograph does not exceed 1 km, which allows assuming that both techniques recorded the earthquake at the same time. However, the GNSS receiver records time in GPST, while the seismograph records time in Universal Time Coordinated (UTC), and hence, the time difference between those systems in 2016 was equal to 17 s (Quality Positioning Services 2021).

GNSS stations track GPS and GLONASS with an observation interval of 0.1 s (10 Hz). The only exception is the ACCU station, which observes GPS only. All seismological instruments belong to the Italian Strong Motion Network. The study is performed for all stations listed in Table 2. The analyses are conducted for all satellites of the GPS constellation for which there is continuity of observation and the elevation angle is greater than 10° (G05, G07, G13, G20, G28, G30). For simplicity, the methodology and detailed results are shown for satellite G05 from the ACCU station.

In terms of the PPP static solution, the estimated parameters, such as coordinates, stabilize after the initialization phase. The movements of stations cannot be determined in the coordinate time series when the static PPP is applied. Thus, the station coordinate estimates from the PPP kinematic solution are commonly used for the station movement detection instead of static positioning. The proposed methodology assumes usage of the phase GNSS observation residuals that come from the PPP processing conducted in the static mode. Observation residuals in both relative and absolute positioning incorporate all signal deviations (Larocca et al. 2006). The station movement cannot be detected in the coordinates due to position constraining but can map into the observation residuals which absorb station displacements or vibrations. Figure 2 shows the proposed methodology in the form of a diagram. All steps are described in the following section.

Figure 2.

The Canadian Spatial Reference System Precise Point Positioning (CSRS-PPP) service (Banville 2020) is employed to derive PPP solutions. The processing the in static mode is used to obtain GNSS observation residuals employed in the proposed methodology. Additionally, we process observation data in a kinematic mode to show time series in the position domain. Table 3 shows the basic parameters for the computational strategy. Calculations were performed for L1 and L2 frequencies for nine-hour observation files (11:00 to 20:00) at 10 Hz.

The occurrence of the analyzed earthquake is recorded by seismographs, cataloged, and also analytically confirmed. Figure 3 illustrates the time series of coordinates derived from kinematic PPP. The obtained results allow visually confirming of vibration occurrences with coordinate differences for three coordinate components of the ACCU station. The anomaly of receiver position occurs between 19:18:30 and 19:19:00. Additionally, the solid vertical line in Figure 3 indicates the original seismograph range of runtime which coincides with the coordinate changes confirming that the changes in coordinate values are caused by the earthquake. However, due to the fluctuations in GNSS kinematic positions, automatic detection of the moment of the earthquake is a non-trivial task.

Figure 3.

The first step of the proposed methodology includes data preprocessing and filtering. The CSRS-PPP residual output files are trimmed to the time range of the earthquake occurrence. The trimmed data cover the period of 90 s starting 30 s before the time of the earthquake beginning. Additionally, data for satellites with residual discontinuities and for satellites at elevation angles below 10° are also removed. The trimmed GNSS observation residual time series is next used in the earthquake detection procedure.

Figure 4 shows the Fast Fourier Transform for the unfiltered residuals of the GNSS phase observations at the L1W signal for all satellites. The time series consists of the main signal and noise. The frequency bounds on the main signal were empirically determined. The frequency limit equal to 0.1 Hz was used as the cut-off for the Butterworth high-pass filter.

Figure 4.

The filtering result is visible in the top panel of Figure 5, where the unfiltered residual time series (red) and the residual time series after filtering (blue) are illustrated. For the original series, the occurrence of anomalies in the observational residuals is similar to the coordinate difference obtained in the kinematic calculation mode shown in Figure 3.

Figure 5.

The bottom panel of Figure 5 shows the scalogram of the fast wavelet transform for the filtered time series. The wavelet transform divides the signal into individual frequency components, called scales, and is here illustrated in the form of a scalogram. A scalogram is a graphical representation of a matrix containing the powers for a given epoch. It consists of 3 axes: X-axis – time, Y-axis – scale, Z-axis – scale factor value for each epoch and scale (frequency).

Note that the frequency is inversely proportional to scale, so high-frequency values are shown at the bottom of the scalogram. The wavelet transform allows for the detection of low-frequency motions even in very high-frequency regions (Shan et al. 2019). This makes it possible to determine the frequency range of the anomaly being searched for.

The steps described so far do not allow confirmation of the time of earthquake occurrences. The occurrence of an earthquake can be visually detected on the processed scalogram as shown in the bottom panel of Figure 5. However, from a usability point of view, it is important to detect the earthquake automatically. The developed method assumes that the occurrence of an earthquake is related to a high power on the scalogram, which is observed as an outlier or anomaly, to automate the processing. As an anomaly on the scalogram, we consider values exceeding three standard deviations from the arithmetic mean of all absolute coefficients obtained from the continuous wavelet transform for a given satellite, assuming a normal distribution of the residuals. Those points are identified with scalogram outlying values.

The detected points on the scalogram for the analyzed data are shown in the top panel of Figure 6. The occurrence of these outlier points coincides with the occurrence of the anomalies, which confirms that the function correctly determined the moment of the earthquake. Anomaly is detected in the middle part of the scalogram (frequency: 0.03–0.015 Hz; scales from 20 to 40).

Figure 6.

The bottom panel of Figure 6 shows what percentage of the scale factors in an epoch are classified as outliers. At maximum, 40% of the values in individual epochs are outliers. It is worth mentioning that each epoch is divided into 64 scalars out of the 1201 epochs in total. Thus, a total of 76,864 scalar factors were analyzed.

As shown in Figure 6, the detected epoch of the earthquake coincides with the starting epoch from the earthquake catalog. Next, the scale factor values for all satellites are summed, and outliers are again detected for the summed scale factor matrix to strengthen the anomaly detection procedure. We also calculate the percentage of outliers in each epoch or the sum of residuals from all satellites. The results are shown in Figure 7. The methodology detects three clusters of outliers, but the vast majority coincides with the occurrence of the earthquake. The remaining two clusters are small, and the outliers cover a maximum of 25% of the individual scales. Additionally, those two regions are detected for scales higher than the scales of the actual earthquake.

Figure 7.

The methodological steps described so far confirm that the earthquake was observed in residuals of observations from the single satellite G05 (Figure 6) and grouped residuals of observations from all satellites (Figure 7). Unfortunately, it is unclear from (Figure 7) what is the impact of each satellite on recorded anomalies. There is a possibility that the residuals for some satellites did not detect the earthquake at all and therefore have no residual outliers that would affect the results. Thus, a second approach was tested in which we analyze how many satellites provide consistent series of anomalies. A wavelet transform is performed for each satellite separately, and outlier observations are detected.

Figure 8 shows how many satellites have detected an outlier at a given time. It can be seen that all satellites detected the moment of the earthquake at the same time. Only one satellite has outliers at other time intervals. The additional series of outliers for a single satellite are not related to the analyzed earthquake. Furthermore, the scale range for the points detected incorrectly is in the range above 30, while most of the anomalies detected on the scalogram during the actual earthquake period cover scales from 20 to 40.

Figure 8.

The final stage of the study is to determine whether the detected period coincides with the official time of the earthquake according to the seismic catalog. The start and end of the earthquake times are determined by taking the time of the first and last outliers from the main grouping of outliers from the first approach. For the second approach, the start and end of the earthquake times are determined by taking the time of the first and last outliers detected by a minimum of five satellites. The determined times were plotted for the north component of the seismograph located at the same site (Figure 9). Both proposed approaches correctly detected the moment of the earthquake. According to the first method, the moment of earthquake onset occurs 8 seconds after the catalog earthquake onset time, while according to the second method, it is 9 seconds.

Figure 9.

This article describes a new methodology aimed to detect the moment of a short-period earthquake. The calculations are based on phase residuals of the GNSS observations from static PPP solutions, which have not previously been used in GNSS seismology. It was found that the methodology correctly detects the earthquake moment.

The work shows that the GNSS data can reflect the movement of the GNSS receiver located in the area affected by the earthquake. Even though the static coordinate time series cannot reflect tremors, the GNSS residuals contain this information and could be used for the earthquake detection. The considered approaches have shown that wavelet scalograms prepared for both individual and summed residual time series can reveal information about site displacements. It should be noted here that the proposed method is not intended to undermine the quality of the coordinate series obtained by the PPP technique but to prove that the residual series contains information that can be used as a supplement to the classic PPP technique for the earthquake analysis. Further research studies should be directed toward improving the methodology described, e.g., the influence of time series discontinuities would have to be eliminated. This could be done by reducing the weights for the residuals observed at low elevation angles or by automatically selecting the cut-off frequency for Butterworth filtering. Alternatively, one may employ an iterative approach with the elimination of residuals that do not meet the assumed accuracy criteria. It would also be necessary to consider how to weight the observations based on the satellite block as the quality of the onboard equipment, e.g., atomic clocks, could also affect the results.

Further research studies should also be conducted to determine the suitability of the residuals to determine vibrations caused by earthquakes of different origins or magnitudes. Possibly, the developed method will detect short-term displacements of engineering structures efficiently, thus extending the potential of the conducted analyses. The analyses were performed to prove that the residuals record the same phenomena as the coordinate time series. In addition, some anomalies were only recorded by specific satellite observations but, due to their magnitude, were not included in the coordinate series. However, it is hard to determine the values of the possible point displacement caused by the earthquake. Due to terrain conditions, satellite elevation angle, and equipment quality, the residuals for each satellite vary and do not always correctly identify the earthquake.

Studies have shown that residual values are influenced by measurement errors and monitored events. This confirms the usefulness of GNSS residuals in studies for the time determination of earthquake occurrence.