Real-Time Seismic Waveforms Estimation of the 2019 MW = 6.4 and Mw = 7.1 California Earthquakes With High-Rate Multi-GNSS Observations

Global Navigation Satellite System (GNSS) is recognized as an effective tool to retrieve high-precision seismic displacements, which has been widely used in earthquake early warning systems such as magnitude estimation and fault slip inversion. In this study, we present two multi-GNSS positioning models for rapidly capturing real-time coseismic displacements using the raw pseudorange, carrier phase and Doppler measurements, called constant velocity (CV) and constant acceleration (CA) dynamic precise point positioning (PPP) models, respectively. The proposed methods are validated with the datasets collected during the 2019 July 4 Mw 6.4 and July 6 Mw 7.1 earthquakes in real-time scenarios. Results show that the two models can provide accurate displacement waveforms with the accuracy of few centimeters. Besides, the multi-GNSS integration can significantly improve the model performances compared with GPS-only solutions in both time- and frequency- domains. The dynamic PPP models can also reliably recover the moment magnitudes and permanent seismic displacements, indicating that the methods have the potential to benefit for earthquake early warning and rapid geohazard assessment.


I. INTRODUCTION
The 2019 July 4 M w 6.4 and July 6 M w 7.1 earthquakes occurred in eastern California, southwest of Searles Valley, which are caused by the shallow strike-slip faulting in the crust of the North America plate [1]. Compared to the accelerometers, whose displacements by the double integration are likely to be distorted by the instrumental rotations and tilts, Global Navigation Satellite System (GNSS) is an effective and powerful tool to retrieve high-precision displacements for earthquake early warning [2], [3]. The current studies chiefly focus on Global Positioning System (GPS) seismology [4], [5]. In the GNSS era, the performances of multi-GNSS positioning for seismological applications are expected to be investigated.
The associate editor coordinating the review of this manuscript and approving it for publication was Chong Leong Gan .
Relative positioning is an early GNSS approach for realtime retrieval of seismic displacements, whereas it needs quantities of stations for simultaneous observations and requires assigning baselines or overlapping networks [6]. Precise point positioning (PPP) is flexible and able to provide absolute displacements in a global reference frame with a stand-alone receiver [7]. Approximately 30-60 minutes convergence time is required for PPP because of the estimable float phase ambiguities and it still takes 15-20 minutes to converge when the ambiguity resolution technology is applied [8].
In the recent past, Colosimo et al. [9] proposed a variometric approach to estimate real-time seismic displacements based on time differences of the GNSS observations and reconstruction from the velocity integration, which has been adopted for many seismology studies [10], [11]. Besides, Li et al. [12] present a temporal point positioning (TPP) method to overcome the long convergence problem in PPP by estimating real-time seismic displacements with a standalone receiver, whose method was successfully verified by the Tohoku-Oki earthquake data in Japan. Tu [13] as well presented the GPS PPP velocity estimation method without the long convergence time to retrieve the displacements, which has been adopted by Jin and Su [14] to capture the real-time seismic displacements in 2018 Alaska earthquake.
In this paper, we extend two dynamic PPP approaches to real-time rapid estimate seismic displacements with a standalone multi-GNSS receiver and equip them into the seismic waveforms and earthquakes magnitude estimation applications [15]. The methods of the constant velocity (CV) and constant acceleration (CA) dynamic PPP models are described in Section 2, together with strategies of precise orbit determination (POD) and precise clock estimation (PCE) in real-time positioning system in Section 3. Section 4 validated the displacement precision of proposed methods before and during the earthquake. The real-time GPS-only and multi-GNSS derived displacements are also compared in the frequency domain. The earthquake magnitudes are estimated and evaluated arisen from the proposed methods and the accuracy of permanent displacements of the earthquakes by the two models is analyzed. Finally, conclusions are given in Section 5.

II. REAL-TIME POSITIONING SYSTEM AND METHODS
The Multi-GNSS Experiment (MGEX), set up by the International GNSS service (IGS), can track and analyze the available GNSS signals, including the GPS, Beidou, GLONASS and Galileo systems, as well as Quasi-Zenith Satellite System (QZSS) and Indian Regional Navigation Satellite System (IRNSS). Over the past years, quantities of multi-GNSS stations were established and integrated into the MGEX reference stations [16]. For a complete real-time earthquake monitoring and early warning system, three aspects are necessary, including POD, PCE and precise positioning, respectively. Under the system of IGS MGEX, real-time data from the global network of reference stations are available and precise orbits and clocks generating can be guaranteed in real-time scenarios [17].
For the user-end, when the real-time orbit and clock products and the high-rate GPS, GLONASS and Galileo observations are available, the seismic waveforms can be retrieved using the precise positioning technology. CV motion model is considered as a simple and ideal model, where the position noise is a random variable that changes with time and the state vector follows the first-order Gaussian-Markov process. When the receiver epoch interval t is less than 1 s, the CV motion model assumption can be generally recognized, whose linear state equation can be written as [18]: where x andẋ denote the three-dimension position and velocity vectors, receptively; I m denotes the m × m identity matrix; w denotes the process noise matrix; k denotes the corresponding epoch; q v denotes the dynamic noise of velocity; ⊗ denotes the Kronecker product. Further than that, CA motion model assumes that a steady acceleration exists between two adjacent epochs by the second-order Gaussian-Markov process, in which the variations of displacement, velocity, and acceleration are random. In general, the time interval of CA motion model is 1-3 s and the corresponding linear state equation can be expresses as [19]: whereẍ denotes the acceleration vector; q a denotes the dynamic noise of the acceleration. The GNSS pseudorange, carrier phase and Doppler observables of n satellites at a single epoch for CV and CA dynamic PPP can be expressed as an unified form, that is [20]: where P, and D denote the observed minus computed pseudorange, carrier phase and Doppler observables vectors; dx and dẋ denote three-dimension receiver position and velocity increments vectors; B denotes the corresponding design matrix; Z w denotes the zenith wet delay (ZWD) vector; W denotes the tropospheric wet mapping function vector; dt r and dṫ r denote the receiver clock offset and drift vectors; e n is the n-row vector with all 1's values; DCB denotes the receiver differential code bias (DCB) vector between pseudoranges on the first and second frequency; τ denotes the slant ionosphere parameters vector; τ 0 denotes the vector of ionospheric prior observations and can be derived from external ionospheric products, such as global ionospheric maps (GIMs); a denotes the vector of the ambiguities; a 1 and a 2 are the vector of the first and second frequency carrier phase ambiguities; f 1 and f 2 denote the frequency values; = diag(λ 1 , λ 2 ) denotes the waveforms matrix; ε denotes the corresponding observations noises; σ denotes the corresponding variance factors; Q 0 denotes the cofactor matrix; The most significant contributors to the Doppler observation are the receiver velocity and clock drift. Generally speaking, the satellite clock drift is very small and atmospheric and hardware delays do not contribute to the Doppler for they vary slowly with time so that they are neglected in the adjustment.
In two dynamic PPP models, the estimated displacements are discarded for they are biased by the dynamic noises of the models and the operational displacements are obtained by the integration and can be expressed, respectively, as:  [ where x denotes the integrated displacements and k 0 denotes the beginning epoch. Thus far, the real-time dynamic PPP solutions can be guaranteed with the estimated orbit and clock products. The data flow of the real-time positioning system for the serverand user-end, including the POD, PCE and PPP solutions, is shown in Fig. 1.

III. DATA AND PROCESSING STRATEGIES
The University Navstar Consortium (UNAVCO) provides access to GNSS datasets from thousands of globally distributed stations. The multi-GNSS observations from MGEX network stations are used to estimate the real-time orbit and clock products, which are shown in Fig. 2, together with the data of high-rate UNAVCO GNSS stations within 300 km from the M w 6.4 and M w 7.1 earthquakes' epicentres. The selected UNAVCO GNSS stations are used to perform realtime PPP solutions. The coordinates of MGEX stations are tightly constrained to the reference values in POD and then the satellite orbit and station coordinates are tightly constrained in PCE processing. The CA and CV dynamic PPP solutions are applied for 1 Hz and 5 Hz data, respectively, according to the assumptions in Section 2. Tab. 1 summarizes the GNSS data processing strategies for POD, PCE and PPP. Meanwhile, the dynamic models in multi-GNSS POD can refer to Li et al. [21].

IV. RESULTS AND ANALYSIS
The section starts with the precision assessment of displacements for CA and CV dynamic PPP models. Then, the coseismic waveforms from two PPP approaches during the  earthquakes are compared. The power spectral densities (PSDs) on the frequency bands are also analyzed. Next, the extractive peak ground displacements (PGDs) are used to estimate the earthquake Richter magnitudes. Finally, the accuracy of the dynamic PPP derived permanent displacements is analyzed.

A. PRECISION ASSESSMENT IN REAL-TIME SCENARIOS
In the place near Ridgecrest, California, an M w 6.4 earthquake occurred at 17:33:49 Coordinated Universal Time (UTC) on July 4, 2019 and then an M w 7.1 earthquake occurred again at 03:19:53 UTC, July 6. The simulated realtime CA and CV dynamic PPP solutions are carried out to retrieve the displacements during the events. Taking the multi-GNSS station P616 during the M w 6.4 earthquake as a case, Fig. 3 shows the estimated displacements derived from the conventional kinematic, CA and CV dynamic PPP models in the 3 h interval from 16:00:00 to 19:00:00 UTC on July 4, 2019. For conventional kinematic PPP, it needs dozes of minutes to achieve few centimeters accuracy. As a contrast, after subtracting the systematic bias arisen from initial integrated errors, the two models can both determine the precise displacements for high-rate data without suffering the long convergence time. This happens because the corresponding motion changing parameters are dependent with the state equations and irrelevant for unconverged ambiguities [14]. As expected, the precision improvement of multi-GNSS PPP is obvious, when compared with the GPS-only solutions. To explain the reason, the visible GPS/GLONASS/Galileo satellites at station P616 during 16:00:00 to 19:00:00 UTC are shown in Fig. 4. It can be seen that more than 30 GNSS satellites can be observed, compared to about 15 GPS visible satellites. The average observed GPS, GLONASS and Galileo satellites number are 9.5, 6.9 and 4.6, respectively. Increasing satellite number can contribute to the better satellite geometry condition and result in the more accurate displacements estimates [27].
To assess the precision of displacements by the two dynamic PPP models, the GPS-only and GPS/GLONASS/ Galileo solutions using datasets collected on day of year (DOY) 185 and 187 during the period of the earthquakes are compared. In this study, one-hour positioning solutions prior to each earthquake are used to analyze the displacement precision. Altogether 802 and 577 sets of two-day results are utilized to get statistical values on the displacement precision for 1 Hz and 5 Hz datasets, respectively. Fig. 5 depicts the boxplot of the displacement repeatability in the north, east and up components for GPS-only and GPS/GLONASS/Galileo CA and CV dynamic PPP solutions, where the boxplot denotes the 25th to 75th percentiles with the median percentile and the whiskers denote the 5th and 95th percentiles. The median displacement repeatability values are also shown in Fig. 5. It is apparent that the displacement precision performances of two models agree well with each other, for instance, of which the median values of GPS-only dynamic PPP displacement repeatability are (1.78, 1.65, 3.52) cm and (1.59, 1.68, 3.38) cm for CA and CV solutions, respectively. The derived vertical components have the noisiest level due to the correlation between up component and neutral atmospheric delay and the satellite geometry. In general, a displacement precision of better than 2 cm and 4 cm in horizontal and vertical parts can be achieved with a confidence 95% level by the two dynamic PPP solutions and the multi-GNSS integration can exhibit the higher precision displacements.

B. COMPARISION OF COSEISMIC WAVEFORMS DURING THE EARTHQUAKES
To better validate the dynamic PPP capability of capturing the seismic waveforms during the earthquakes, 1 Hz and 5 Hz high-rate data of stations P595, CCCC, P580, RAMT and P593 are selected as the examples to analyze the horizontal and vertical seismic waveforms by dynamic PPP in real-time scenarios. Fig. 6 shows the acceleration, velocity and displacement waveforms with 300 s intervals around the M w 6.4 and M w 7.1 earthquake events with multi-GNSS observations using CA and CV dynamic PPP solutions for 1 Hz and 5 Hz data, respectively. The selected five stations for both events are arranged in terms of their hypocentral distances. The P and S waves expected arrival times are also marked with red and blue colors, respectively, in Fig. 6, which are calculated using the empirical values of compressional and shear crust velocities (Vp = 5.5 km/s and Vs = 3.2 km/s) [28]. As shown, P waves arrived first with a higher speed, followed by S waves. The CA dynamic PPP displacement waveforms agree with that of the CV solution in horizontal components with similar amplitudes and aligned phases, whereas the vertical waveforms are likely to disagree with some stations. The phenomena are reasonable for the 1 Hz and 5 Hz data are applied by two models. Besides, the analysis of the arrival times of calculated seismic wave reveals that the horizontal peak GNSS derived displacement waveforms are caused by the S wave rather than the P wave.
To verify the reliability and precision of the displacements by the dynamic PPP approaches during the earthquake, the solutions calculated by postprocessed multi-GNSS kinematic PPP are used as the ground truth, which has been adopted from Geng et al. [11]. For the selected five stations, the standard deviation (STD) of differences of the north, east and up components are (0.22-0.35, 0.44-0.58, 0.58-0.82) cm and (0.08-0.27, 0.16-0.27, 0.21-0.69) cm for dynamic CA and CV PPP solutions, respectively. The dynamic CA and CV PPP approaches are capable of capturing displacement waveforms with a reliable and convincible precision of better than 1 cm during the earthquake seismic period. The narrow differences of the methods are caused by the corresponding dynamic noises set in the state equations of the PPP models.

C. COMPARISION OF THE DISPLACEMENTS IN FREQUENCY DOMAIN
To investigate the dynamic PPP derived 1 Hz and 5 Hz displacements in the frequency domain, the PSDs of the displacements by GPS-only and multi-GNSS solutions at stations P616 and TONO from 02:00:00 to 05:00:00 UTC with the 3 h length in the north, east and up components are shown in Fig. 7. The PSDs of two station displacement by different methods exhibit nearly the same general shape. The PSDs decrease with increasing slope at the low frequency range. Then, the PSDs start to decline with an approximately linear trend. Exceeding the above transition frequency states, another swinging variation status takes place again. Comparable noise levels of the displacements exist in the north and east components. The noise levels of the vertical component are higher than the horizontal components, which are affected by the precision of displacements. Besides, the multi-GNSS integration can reduce the noise level during the entire frequency bands in the three components, especially obvious during the frequencies from approximately 10 −3 to 10 −1 Hz. In contrast, the decline of PSDs on the higher frequencies (larger than 10 −1 Hz) is slight and PSDs have the slight fluctuations. The PSDs analysis also indicates that the multi-GNSS integration have the better displacement precision with the lower noise level. The analysis of the displacement in the frequency domain shows that the proposed methods are reliable.

D. MAGNITUDES ESTIMATES OF THE EATHQUAKES
The magnitude is an important source parameter to measure the earthquake energy in the seismology research.
The earliest magnitude scales defined by Charles Richter referred to the Richter magnitude [29]. GNSS can provide the displacements with a high sign-to-noise ratio (SNR) level, whereas the broadband seismometers are likely to clip in the near field and the strong-motion accelerometers could clip in case of ground accelerations of 2-3 g depending on the instrumentation in theory. Based on the quantities of earthquakes recorded by the strong-motion accelerometers, Gutenberg [30] given an empirical formula to get the Richter magnitude value, which can be expressed as: where M denotes the estimated Richter magnitude, A denotes the peak ground displacement (PGD) during the earthquake in unit of micrometre, whose extracting method can refer to Fang et al. [31]. denotes the hypocentral distance in unit of degree.
In term of the Equation (7), the parameters relationship can also be organized as: For the earthquakes, the set of the PGD values from the corresponding GNSS stations is retrieved in real-time scenarios. The distance of the area affected by the earthquakes is related to the earthquake magnitude. In view of the status and reflection of GNSS stations, the stations with hypocentral distances VOLUME 8, 2020  less than 100 and 300 km, respectively, are selected for magnitude determination of the M w 6.4 and M w 7.1 earthquakes. Fig. 8 shows the PGD values of the selected stations plotted as the function of hypocentral distances in a log-log scale for the M w 6.4 and M w 7.1 earthquakes with the dynamic CA and CV solutions, respectively. The distribution of the calculated magnitudes with two methods is also plotted in Fig. 8 Fig. 9 presents in map view the permanent displacements corresponding to the M w 6.4 and M w 7.1 earthquake events on the horizontal and vertical components released by USGS using the Kalman estimation. The offsets are estimated from the accelerometer and GNSS daily solutions between both sides of the events from Kalman filter time series analyzes [1]. The black scale arrow represents 100 mm length. For July 4 M w 6.4 earthquake, maximum horizontal and vertical displacements are approximately 10 cm and 2.5 cm, towards to the east-north and descent directions, respectively. For July 6 M w 7.1 earthquake, the permanent displacements are much larger and affect farther regions than the previously event. For instance, the GNSS station P580, approximately 40 km away from the epicenter, moved -2.4, 19.9 and 0.5 cm for the north, east and up components in the earthquake, respectively. The maximum displacement of the GNSS station is more than 0.5 m at station P595. The both events are east-west extensional and north-south contractional, indicating a right lateral strike slip event on a northwest striking plane.
To evaluate the reliability of permanent displacement estimation by the dynamic PPP solutions, the offsets are computed by the subtracting the average of 30 min length displacements of 20 min before and after the events. The unavailable offsets for the stations are assumed static in that they are not affected by the earthquake, most of them are concentrated in the M w 6.4 event. Displacements of them are compared and the distribution of differences is plotted on Fig. 10. In terms of the RMS, agreements of the north, east and up components are (1.54, 2.72, 2.89) cm and (1.72, 2.77, 2.94) cm for the CA and CV dynamic PPP solutions, respectively. All plots of Fig. 10 are likely to exhibit normal distribution. It is vividly demonstrated that the permanent displacement derived by the two methods has a reliable accuracy, which certainly can be employed for deriving offsets for real-time early warning and geohazard monitoring with an accuracy of few centimeters.

V. CONCLUSION
This study presents two dynamic PPP approaches in realtime scenarios for the stand-alone multi-GNSS receivers. The high-rate GNSS observations recorded during 2019 July 4 Mw 6.4 and July 6 M w 7.1 Earthquakes near Ridgecrest, California are used to analyze the capability of GPS and the benefit of multi-GNSS for rapid determining precise displacements. We demonstrated that the GPS-only PPP can achieve a displacement precision of (1.78, 1.65, 3.52) cm and (1.59, 1.68, 3.38) cm in the north, east and up components for CA and CV PPP solutions, respectively, and multi-GNSSderived solutions show higher precision, which can contribute to the better geometry. The analysis of the displacements in the frequency domain also indicates that the multi-GNSS integration can improve the displacement precision. During the earthquakes, the seismic displacement waveforms derived from GNSS dynamic PPP have a precision of better than 1 cm. The Richter magnitudes derived from extracted PGDs from the proposed methods agree well with the reference magnitudes with the precisions of 0.3. Moreover, the derived permanent displacements agree well with the theoretical displacement with an accuracy of few centimeters. The above results indicate that the high-accuracy seismic waveforms can be captured in real-time scenarios using the GNSS observations by our proposed methods, which has a wide prospect in GNSS seismology applications such as fault slip inversion and magnitude estimation.