https://orcid.org/0000-0003-2988-4850
SUMMARY
Near-field large-amplitude seismograms are essential for the rapid inversion of earthquake source parameters using waveform inversion methods such as the Cut And Paste (CAP) method for disaster assessment. High-rate Global Navigation Satellite System (GNSS) relative positioning (RP) provides precise, rapid, and real-time measurement of near-field large-amplitude displacements. However, RP records motion with respect to a reference station, and the reference station's movements become part of the relative displacement waveforms. Therefore, seismic source parameter estimates may be inaccurate if the reference station's motion is not taken into consideration, and doing so affects some basic assumptions made in the CAP method. To overcome this problem, we develop an expanded differential CAP inversion approach specifically for high-rate GNSS RP (CAP-RP) that accounts for the motion of the reference station. Two methods are proposed to implement CAP-RP: an expanded differential CAP (D-CAP) and an iterative post-processing CAP (P-CAP). We assess the performance of CAP-RP with different data sets, using the July 2019 Mw 6.4 earthquake in California as a case study. Both CAP-RP techniques produce accurate source parameters in synthetic data inversion tests, indicating the feasibility of the strategy. However, P-CAP is more time-efficient than D-CAP, making it the better option. Generally, results from high-rate GNSS RP, broad-band seismographs, and their inverted combinations exhibit consistency in observational data inversion testing. Our results also demonstrate that more accurate source parameters can be obtained by combining sensitive far-field broad-band seismograph data with large amplitude near-field GNSS RP waveforms.
1. INTRODUCTION
High-rate Global Navigation Satellite System (GNSS) data can record large amplitude seismic waves, playing a crucial role in seismological research (Freymueller et al. 1994; Bock et al. 2004; Blewitt et al. 2009; Guo et al. 2019; Xu et al. 2019). With improvements in theory and in the precision of satellite orbit and clocks products, the accuracy of high-rate GNSS solutions has commensurately improved significantly (Kouba & Héroux 2001; Li et al. 2020). This enhancement has facilitated the successful application of high-rate GNSS in seismology, gradually evolving into the field of GNSS seismology (Hirahara et al. 1994; Bilich et al. 2008; Larson 2009). With rigorous algorithmic processing, high-rate GNSS can obtain precise displacement, velocity and acceleration solutions (Laurichesse & Privat 2015; Wang et al. 2020; Xu et al. 2021), as well as ionospheric total electron content (Heki 2011, 2024). High-rate GNSS is not easy to clip and the displacement estimates are not affected by tilt or rotation, making it suitable for capturing near field strong ground motion in the frequency range of 0 Hz to the Nyquist frequency (Larson et al. 2003; Bilich et al. 2008). These strong ground motion observations can effectively record the complex rupture process of earthquakes, which plays a key role in the study of source parameter inversion (Jiang et al. 2019; Shan et al. 2019; Xu et al. 2021; Chen et al. 2022).
In high-rate GNSS data processing, precise displacements can be obtained using two basic methods: Relative Positioning (RP) and Precise Point Positioning (PPP), and both depend on the availability of precise satellite orbits. Today, near-real time and even predicted satellite orbits are readily available from a variety of sources. RP typically uses differences between the observations of two stations to determine their relative motions. PPP (Zumberge et al. 1997) uses precise satellite clock estimates from a prior solution to enable single station positioning, and is widely used in source parameter inversion (Larson et al. 2003; Sheng et al. 2020; Chen et al. 2022). When combined with Ambiguity Resolution (AR) techniques (Collins et al. 2008; Ge et al. 2008; Laurichesse et al. 2009; Bertiger et al. 2010), PPP-AR can provide post-processed positioning services with mm∼cm accuracy (Geng et al. 2019; Zhang et al. 2023), once precise orbit and clock products are available; however, real-time PPP-AR is limited to a precision of a few to several centimeters because of the lack of precise real-time clock estimates.
High-rate PPP is often used to capture large amplitude ground motions in the near field, and the PPP waveforms can be directly applied to source parameter inversion without any conversion. However, unlike data from seismographs, there can be a considerable delay in obtaining accurate high-rate PPP waveforms. Rapid orbit and clock products are available only one to two days after data collection, and final orbit and clock products up to one to two weeks later. This is because high-rate PPP waveforms require support from products such as precise satellite orbits and clocks, while PPP-AR additionally necessitates AR products (Laurichesse & Privat 2015; Zhang et al. 2023). These products are produced and released by analysis centres, which collect observations from globally distributed GNSS stations and process them to obtain precise satellite products (Kouba & Héroux 2001; Schaer et al. 2021). This delay hinders high-rate PPP from meeting the application requirements for the rapid acquisition of focal mechanisms.
The RP method can rapidly and even in real time obtain relative displacements at the mm∼cm precision. Relative displacements are the difference between the absolute displacements of the rover station and a reference station. This method eliminates various errors such as satellite orbit errors, satellite clock errors, and atmospheric errors through differential processing (Dong & Bock 1989; Bisnath 2020). However, the choice of the reference station in high-rate GNSS RP is critical, as this will directly affect the relative displacement results. On the one hand, the baseline length between the rover and reference station affects the accuracy of the RP results. When the baseline is excessively long, the correlation of errors (such as atmospheric delay) decreases. As a result, differential processing becomes less effective at eliminating these errors, which can significantly reduce the accuracy of RP results (Blewitt 1989; Dong & Bock 1989). On the other hand, the motion of the reference station can impact the measurement of absolute displacement at the rover station. If the reference station is situated close to the epicentre, its significant displacement can interfere with the waveforms at the rover station (Zheng et al. 2012; Yin & Wdowinski 2014; Guo et al. 2015). The accuracy of inverted source parameters will diminish when applying the disturbed relative displacements directly. Therefore, research on inversion methods that adapt to the characteristics of high-rate GNSS RP displacement time-series is imperative.
The Cut And Paste (CAP) method is widely used in studies of earthquake source parameters (Zhao & Helmberger 1994; Zhu & Helmberger 1996). This approach cuts seismograms into five segments, fitting them individually while correcting for time-shifts between the different segments caused by velocity model errors. This enables CAP to achieve more stable and reliable source parameters (Chen et al. 2015; Guo et al. 2015; He and Ni 2017; Sheng et al. 2020). High-rate GNSS RP waveforms have long been applied to CAP source parameter inversion (Bock et al. 2011; Zheng et al. 2012; Zhang et al. 2021). To reduce or eliminate interference resulting from the movement of the reference station, a station that is far from the epicentre can be chosen as the reference station. If the reference station experiences minimal motion, then the motion of the reference station can be disregarded (Bock et al. 2011; Zhang et al. 2021). However, the long baselines that result from such a choice can negatively impact the accuracy of RP results. When the reference station moves significantly, the differing arrival times of seismic waves at the rover and reference stations may allow for the waveforms to be distinguished, but the motion of the reference station always will appear somewhere in the RP displacement time-series (Zheng et al. 2012). To address this concern, previous research has concentrated on data processing techniques aimed at minimizing or eliminating the effects of reference station movement. For instance, Bock et al. (2011) uses a dense GNSS network solution to obtain displacements relative to stations situated outside the deformation area. Yin & Wdowinski (2014) proposes a stacking spatial filtering method that eliminates the influence of reference station motion. However, the effectiveness of these methods diminishes when the GNSS network is sparse (Blewitt 1989; Yin & Wdowinski 2014). Moreover, employing spatial filtering may lead to errors due to waveform differences at various rover stations, and manually detecting reference station movement can be challenging. Therefore, from the perspective of inversion method, we propose an extended differential CAP source parameter inversion method for high-rate GNSS RP (CAP-RP).
This paper begins by presenting the theoretical framework for earthquake source parameter inversion specifically for RP waveforms, followed by an implementation of the CAP-RP method with two different approaches. Next, inversion tests are performed using both synthetic seismic waves and actual observational seismic waves to assess the method's feasibility and reliability. A comparative analysis of the different methods is then conducted using both synthetic and observational data sets. Finally, inversion tests are conducted by combining high-rate GNSS RP waveforms with broad-band seismograph data, highlighting the advantages of integrating different observations.
2. THE CAP METHOD AND THE CAP-RP APPROACH
The CAP method was initially proposed by Zhao & Helmberger (1994), to determine the magnitude, fault angles (strike/dip/rake) and focal depth for a point source model of earthquake. Later, this method was enhanced by Zhu & Helmberger (1996) and Zhu & Ben-Zion (2013), and it has been widely adopted in studies on earthquake source parameter inversion (e.g. Bock et al. 2011; Zheng et al. 2012; Guo et al. 2015). CAP cuts local three-component seismograms (|$r,t,z$|) into five segments (the Cut step), which are matched with corresponding synthetic waveform segments. Each segment is allowed to shift in time to account for the inaccuracy in the arrival times of synthetic seismograms due to errors in the crustal velocity models (the Paste step). The synthetic waveforms are calculated using the average 1-D velocity model through the frequency–wavenumber (FK) integration method (Zhu & Rivera 2002) with given locations of epicentre and stations (or distances and azimuths from the earthquake source). The magnitude, fault angles and focal depth are determined with grid search algorithm.
As shown in Fig. 1, PPP waveforms (Fig. 1b) can be directly integrated into CAP inversion, analogous to seismograph data. However, RP waveforms (Fig. 1c) include the reference station motion, which complicates the use of RP displacements in CAP inversion, and compromises the reliability of the inversion results except in the case where the reference station is far enough away that the motion of reference station can be ignored in any time window of interest. From the perspective of the inversion method, synthetic waveforms should ideally be structured in the form of RP displacements (Fig. 1c) to avoid introducing artificial waveforms. Thus, differencing the synthetic waveforms at the rover and reference stations is necessary, which requires introducing a relative time-shift, as these time-shifts are independently determined at each station in CAP. In addition, the rover station and reference station will in general be located at different azimuths with respect to the source, so the radiation pattern effect on the synthetic waveforms will be different at these two stations. The following analysis details the differential process and determination of relative time-shifts in CAP-RP.
(a) Locations of the earthquake and GNSS stations. (b) and (c) are observations and synthetics in PPP and RP. S1-S2 represents the RP waveform with S2 as the reference station. While precise PPP requires a delay of as much as 1–2 weeks, RP offers rapid, near real-time waveform (S1–S2) with comparable accuracy. PPP waveforms can be directly matched with synthetics, whereas RP waveform should be matched with differential synthetics.
The key to achieve the CAP-RP is to match the synthetic waveform with the observed waveform. As shown in eq. (1), the high-rate GNSS RP displacement is the differential displacement time-series between the rover station and the reference station. Correspondingly, as shown in eq. (2), the differential synthetic waveform should be obtained. In this way, we can match the RP displacements with differential synthetics. Then, we analyse how to match the differenced synthetic waveforms in CAP method, so as to achieve the purpose of CAP inversion with relative motion records.
Where |$\mathrm{ obs}$| and |$\mathrm{ syn}$| represent observed and synthetic waveforms respectively. The superscript |$RP$| represents the relative waveform, m and b represent the rover and reference stations. The subscripts |$r,t,z$| represent radial, tangential and vertical components, respectively. And subscripts |$n,e,z$| represent north, east and zenith components, respectively. |${T_{{\theta _m}}}$| is the rotation matrix from |$n,e,z$| to |$r,t,z$| (https://service.iris.edu/irisws/rotation/), and |${\theta _m}$| is the azimuth of the rover station.
In addition to the time-shift due to the different distances and azimuths from the source to the reference and rover stations, there also is a relative time-shift implied in eq. (2) because of potential errors in the velocity model. If the seismic velocity model was perfectly known, then the application of eqs (1) and (2) would be simple. However, the velocity model is not perfectly known, and because the rover and reference stations will in general have different distances and azimuths from the source, the observations and synthetics will be time-shifted relative to each other and will have different azimuthal dependence (the radiation pattern effect). This means the relative positioning data (and relative synthetics) cannot simply be used in the standard CAP inversion except when the reference station is very far away.
CAP synthesizes three-component (|$r,t,z$|) seismic waveforms based on the 9 component Green's functions, fault angles and station's azimuth. The azimuths of each station are known and the Green's functions are generated by FK before inversion. The fault angles are varied while travelling through each grid search point (this grid is defined by a combination of magnitude, fault angles and time-shifts). Therefore, the synthesis method can be described by eq. (3). The functions |${g_i}$| represent each component of the Green's functions, |${g_{0,3,6}}$|, |${g_{1,4,7}}$| and |${g_{5,8}}$| represent the vertical, radial and tangential Green's functions, respectively, from three fundamental fault types: strike-slip, dip-slip and 45° dip-slip (Helmberger 1983). The |${a_i}$| are the synthesis coefficients, which are determined by the fault angles and station's azimuth (as depicted in eq. 4), reflecting the source geometry and radiation pattern, respectively. Based on the three-component seismic waveform synthesis method described in eq. (3), the differential synthetic waveform in eq. (2) can be expressed as eq. (5). In the final form of eq. (5), we re-factor the differential Green's functions expressions in terms of the ratios of the synthesis coefficients |${a_i}$| between the reference and rover stations, to account for the difference in radiation pattern. The advantage of this re-factoring is that, if the ratio is known, the form of eq. (5) becomes the same as that of eq. (3), except that the Green's functions are now for differential motion.
Where |${a_i}$| are the synthesis coefficients, which are the intermediate products, |${g_i}$| represent the components of the Green's functions, |$\alpha $| represents the station azimuth from the source minus the fault strike, |$\lambda $| is the fault dip, and |$\delta $| is the fault rake. A key point here is that the |${a_i}$| all depend on |$\alpha $| (except for |${a_3}$|), and thus the reference and rover stations have different |${a_i}$| values. |$T_{{\theta _b}}^{\rm{^{\prime}}}$| is rotation matrix from |$r,t,z$| to |$n,e,z$|, and |$g^{\prime}$| is Green's function transformed with |${T_{{\theta _m}}}T_{{\theta _b}}^{\rm{^{\prime}}}$|.
In this paper, we propose two approaches to realize CAP-RP. One is an extended differential CAP named D-CAP; the other is an iterative post-processing CAP named P-CAP. To realize the form of eq. (2) with the unknown relative time-shifts, these shifts must be estimated as well. Generally, the relative time-shifts at different rover stations are not the same. However, independently searching for all relative time-shifts at each station during the inversion process can be very time-consuming, as the search time grows exponentially with the number of stations. Therefore, we assume that the time-shifts between different rover stations and the reference station are the same. This approach reduces the added search dimension to a manageable range within the given relative time-shift limits.
While both the D-CAP and P-CAP methods work for CAP-RP and ultimately give similar results, they differ in how the differential synthetics are computed. In the D-CAP method, a grid search for relative time-shifts (between stations) is incorporated within the original CAP method's grid search. The extended search grid includes ranges for magnitude, fault angles, time-shifts and relative time-shifts, with the relative time-shifts being searched at each possible grid point. Conversely, in the P-CAP method, the differencing of Green's functions is decoupled from the grid search process, reducing its dimensionality. Instead, the differencing operation is performed during the preparation of the Green's functions library, pre-computed with various relative time-shifts, where these shifts are explored for each focal depth.
As shown in Fig. 2, D-CAP conducts the differential operation in the grid search process according to eq. (5). The relative shifts are searched together with magnitude, fault angles, time-shifts. All parameters are known at each grid search point, so the differential form of synthetic waveforms can be obtained directly during the grid search, but at the cost of a heavy computational burden. To enhance computational efficiency, the original CAP separates the computationally intensive cross-correlation calculations and the grid search. However, the D-CAP method has to integrate the Green's functions differencing and cross-correlation operations within the grid search process to account for the relative time-shifts. As a result, the computational cost of D-CAP is much higher than that of CAP, and the computational cost greatly increases with a decrease of the parameter search interval or an increase in the data sampling rate.
Extending differential CAP with source code for high-rate GNSS relative positioning solutions. This approach puts the differencing of Green's functions and the calculation of cross-correlation inside the grid search process.
In order to solve the problem of high computational load and long processing times in D-CAP, we develop an iterative post-processing method we name P-CAP. As shown in Fig. 3, this approach separates the differencing of Green's functions from the grid search operation to reduce the computational load. In P-CAP, the relative time-shifts are searched together with focal depths, while the magnitude, fault angles and time-shifts continue to be searched by the core code. The differencing process is completed during the preparation of the Green's functions library with different relative time shifts. Based on eq. (5), the differenced Green's functions library can be computed using eq. (6).
Developing an iterative post-processing CAP procedure for high-rate GNSS relative positioning solutions. This approach put the differencing of Green's functions (with different relative time-shifts) in the Green's functions library preparation. The coefficients |${a_i}$| and Green's functions library are updated iteratively.
In eq. (6), the ratios of coefficients |${a_i}$| depend on both the fault angles and the station's azimuth, so these ratios are unknown before the inversion. Fortunately, the ratio of |${a_3}$| between reference station and rover station is always 1 because |${a_3}$| is independent of the station's azimuth, but the other ratios will differ from 1 unless the rover and reference stations are located at the same azimuth with respect to the source (eq. 4). According to this characteristic, we initialize all the ratios of |${a_i}$| to 1. The |${a_i}$| ratios can be initialized to different values if a preliminary solution is available, but the inversion converges quickly even if no such information is available. The |${a_i}$| ratios are updated at each iteration, so that the differenced Green's functions library becomes more accurate with each iteration. The convergence conditions for the iterative method will be explored in later tests, and we find that the solution converges rapidly.
3. SOURCE PARAMETERS INVERSION
In July 2019, a Mw 6.4 earthquake occurred in California, USA, the first major event in the 2019 Ridgecrest earthquake sequence. There are abundant broad-band seismograph observations and high-rate GNSS observations in this area, which provide good data support for the source parameter inversion. Previous studies on this event (Liu et al. 2019; Jia et al. 2020) offer a reliable reference for analysing the results. The velocity model in this area is available from Shearer et al. (2005). In this study, 24 high signal-to-noise ratio broad-band seismograph recordings and 23 good-quality 1 Hz GNSS observations are selected. The earthquake location and station distribution are shown in Fig. 4. The broad-band seismograph data were downloaded from IRIS (https://ds.iris.edu/wilber3/find_event), and the high-rate GNSS data were downloaded from UNAVCO (EarthScope) (https://gage-data.earthscope.org/archive/gnss/highrate/1-Hz). We processed the GNSS data using both RP and PPP methods. To obtain stable and reliable high-rate GNSS waveforms, the RP results were processed using the well-known RTKLIB software (https://www.rtklib.com/), while the PPP results were processed using PPP-ARISEN (Zhang et al. 2023), which supports multiple AR products. These results are shown in Fig. 5. In this study, the station with the shortest average baseline length to other stations is selected as the reference station (the green triangle in Fig. 4). The RP solutions are the displacements of the rover stations relative to the reference station, which are equivalent to the displacements difference between other stations and the reference station (red curves in the right panels of Fig. 5) in the PPP solutions. We note that, in Fig. 5, the effect of the reference station motion on the waveforms is significant (as highlighted in the green box). In the PPP displacements, the reference station has a prominent upward arrival in the North component (Fig. 5, top right) at about 32 s; it appears just beneath the grey scale bar in the figure. In the RP displacements, every station has a negative arrival at this same time, and in some cases the negative arrival is the first arrival. This negative arrival results from the fact that the RP solution is explicitly the differential motion of two stations. Although this may be visually confusing to anyone not familiar with relative positions, this study accounts for this effect and accurately applies the RP solutions to the rapid inversion of source parameters.
Earthquake location and station distribution. The red pentagram indicates the epicentre (source information comes from the Global Centroid-Moment-Tensor (GCMT) project); the blue triangles represent GNSS rover stations; the green triangle denotes the GNSS reference station; the red triangles indicate broad-band seismograph stations.
High-rate GNSS results. RP represents the relative positioning results, while PPP denotes the precise point positioning results. In PPP results, the red curves are the displacement time-series of the reference station in RP. The RP results are equal to the PPP results with the red curves subtracted. Consequently, as highlighted by the green box in the N component of RP, the RP results are affected by the motion of the reference station.
3.1 Synthetic inversion tests
Based on the above method and observations, we first performed inversion tests with synthetic data. The synthetic data were generated by the FK method. We generated data for nine virtual stations well distributed in azimuth from 0° to 360°, with epicentral distances randomly distributed between 50 and 150 km. The synthetic waveforms were generated based on a magnitude 6.4 earthquake with a strike/dip/rake of 227°/86°/3°, as reported by GCMT, and a focal depth of 12 km. To obtain synthetic waveforms consistent with high-rate GNSS RP solutions, the sampling rate was set to 1 Hz. One of the virtual stations was selected as the reference station, and the rover stations’ waveforms were obtained by differencing their waveforms with those of the reference station. Based on these synthetic relative waveforms, inversion tests were performed for the above two methods with and without RP noise. Note that in this synthetic test, because we use the same velocity model for the synthetic data and the inversion, the relative time-shifts are all exactly zero because there is no velocity model error.
In the noise-free inversion test, we set the filter frequency range to 0.01–0.45 Hz and the search interval for the fault angles to 3°. The inversion results of the D-CAP and P-CAP methods are listed in Table 1. Both methods yielded source parameters close to the synthetic ones, and the waveforms are fit well. This result preliminarily validates the feasibility of our proposed theories and methods.
Results of synthetic inversion tests.
| Method | Variance | Mw | Depth(km) | Strike(°) | Dip(°) | Rake(°) |
|---|---|---|---|---|---|---|
| Syn | – | 6.4 | 12 | 227 | 86 | 3 |
| D-CAP | 95.3 | 6.38 | 12 | 227 | 86 | 3 |
| P-CAP | 95.5 | 6.37 | 12 | 227 | 86 | 5 |
| Method | Variance | Mw | Depth(km) | Strike(°) | Dip(°) | Rake(°) |
|---|---|---|---|---|---|---|
| Syn | – | 6.4 | 12 | 227 | 86 | 3 |
| D-CAP | 95.3 | 6.38 | 12 | 227 | 86 | 3 |
| P-CAP | 95.5 | 6.37 | 12 | 227 | 86 | 5 |
Results of synthetic inversion tests.
| Method | Variance | Mw | Depth(km) | Strike(°) | Dip(°) | Rake(°) |
|---|---|---|---|---|---|---|
| Syn | – | 6.4 | 12 | 227 | 86 | 3 |
| D-CAP | 95.3 | 6.38 | 12 | 227 | 86 | 3 |
| P-CAP | 95.5 | 6.37 | 12 | 227 | 86 | 5 |
| Method | Variance | Mw | Depth(km) | Strike(°) | Dip(°) | Rake(°) |
|---|---|---|---|---|---|---|
| Syn | – | 6.4 | 12 | 227 | 86 | 3 |
| D-CAP | 95.3 | 6.38 | 12 | 227 | 86 | 3 |
| P-CAP | 95.5 | 6.37 | 12 | 227 | 86 | 5 |
Next, we discuss the convergence conditions for the P-CAP Method. In this inversion, we perform five iterations, and the results of each iteration are shown in Table 2. The behaviours of each |${a_i}$| are illustrated in Fig. 6, depicting the |${a_i}$| ratios between the reference station and rover stations. As shown, stable source parameters are achieved after the second iteration, with variance reduction reaching 95, after which the ratios and variance reduction are stable. Most of the ratios align closely with the theoretical values, as indicated by the proximity of the estimated ratios to the dashed lines in Fig. 6. However, the ratio of |${a_5}$| in Fig. 6(b) shows a large difference from the theoretical value (the theoretical value is 15.62, and the corresponding dashed line is outside of the plot box). This discrepancy does not impact the inversion results, likely due the near-zero value of |${a_5}$| for this particular synthetic waveform. P-CAP converges after two iterations, suggesting that setting the number of iterations to two may be an appropriate convergence criterion. However, further testing and analysis are required. Fig. 7 shows the search results for the relative time-shift and the depth in P-CAP. The optimal inversion results correspond to a time-shift of 0 s and a depth of 12 km, which are consistent with the synthetic parameters. Notably, there is a small difference (within one step) in the slip angle for P-CAP, which may be related to the iterative process and will be further analysed in the discussion section. In this test, D-CAP took approximately 142.74 hr, while P-CAP took about 0.80 hr.
Behaviour of |${a_i}$| ratios at six virtual rover stations. Since all eight virtual rover stations exhibit similar characteristics, only the first (azimuth order) six are plotted. The circles connected by solid lines represent the ratios at each iteration and the dashed lines are the theoretical values. The values stabilize after 1–2 iterations. Not all Green's function components are equally important to reconstructing the final synthetic waveform; their relative importance depends on the specific earthquake source.
Time-shift and source depth grid search. The optimal result is indicated by the pink dot. The left panel shows the results for the entire grid, while the right panel details the variations around the optimal result.
Iteration process of P-CAP.
| Iteration | Variance | Mw | Depth(km) | Strike(°) | Dip(°) | Rake(°) |
|---|---|---|---|---|---|---|
| 1 | 55 | 6.22 | 12 | 234 | 66 | 30 |
| 2 | 95.4 | 6.37 | 12 | 227 | 86 | 5 |
| 3 | 95.5 | 6.37 | 12 | 227 | 86 | 5 |
| 4 | 95.5 | 6.37 | 12 | 227 | 87 | 4 |
| 5 | 95.5 | 6.37 | 12 | 227 | 86 | 5 |
| Syn | – | 6.4 | 12 | 227 | 86 | 3 |
| Iteration | Variance | Mw | Depth(km) | Strike(°) | Dip(°) | Rake(°) |
|---|---|---|---|---|---|---|
| 1 | 55 | 6.22 | 12 | 234 | 66 | 30 |
| 2 | 95.4 | 6.37 | 12 | 227 | 86 | 5 |
| 3 | 95.5 | 6.37 | 12 | 227 | 86 | 5 |
| 4 | 95.5 | 6.37 | 12 | 227 | 87 | 4 |
| 5 | 95.5 | 6.37 | 12 | 227 | 86 | 5 |
| Syn | – | 6.4 | 12 | 227 | 86 | 3 |
Iteration process of P-CAP.
| Iteration | Variance | Mw | Depth(km) | Strike(°) | Dip(°) | Rake(°) |
|---|---|---|---|---|---|---|
| 1 | 55 | 6.22 | 12 | 234 | 66 | 30 |
| 2 | 95.4 | 6.37 | 12 | 227 | 86 | 5 |
| 3 | 95.5 | 6.37 | 12 | 227 | 86 | 5 |
| 4 | 95.5 | 6.37 | 12 | 227 | 87 | 4 |
| 5 | 95.5 | 6.37 | 12 | 227 | 86 | 5 |
| Syn | – | 6.4 | 12 | 227 | 86 | 3 |
| Iteration | Variance | Mw | Depth(km) | Strike(°) | Dip(°) | Rake(°) |
|---|---|---|---|---|---|---|
| 1 | 55 | 6.22 | 12 | 234 | 66 | 30 |
| 2 | 95.4 | 6.37 | 12 | 227 | 86 | 5 |
| 3 | 95.5 | 6.37 | 12 | 227 | 86 | 5 |
| 4 | 95.5 | 6.37 | 12 | 227 | 87 | 4 |
| 5 | 95.5 | 6.37 | 12 | 227 | 86 | 5 |
| Syn | – | 6.4 | 12 | 227 | 86 | 3 |
The synthetic inversion test was repeated with the addition of realistic noise to the RP waveforms. The RP noise comes from the displacement time-series of the rover stations prior to the earthquake (Fig. 3), when there was no real movement. With eight virtual stations’ synthetic relative displacements, noise time-series from the eight rover stations were randomly selected and added to these waveforms to obtain synthetic relative waveforms. Similarly, we set the filter frequency range from 0.01 to 0.45 Hz to include more waveform information, and set the fault angle search interval to 3°. We use the P-CAP method to obtain the source parameters quickly. While the high frequency noise is significant in high-rate GNSS, we use the Radial (R) and Tangential (T) components with high signal-to-noise ratios for the inversion test. After three iterations, we obtain the source parameters and waveform fitting results shown in Fig. 8. These source parameters are consistent with the synthetic parameters, with slight differences (within two steps). This result further validates the feasibility of the CAP-RP method. The slight differences may be related to the iterative convergence process of P-CAP or/and the noise in high-rate GNSS. The RP noise is complex (black curve in Fig. 8), which significantly reduced the variance reduction compared to the noise-free test. However, when the station amplitudes are large, the waveform fitting remains reliable.
Waveform fitting of P-CAP. The three numbers following ‘FM’ represent the strike, dip and rake angles of the fault plane, respectively. The left side sequentially indicates the epicentral distance, station name and azimuth. The black curves are the observed waveforms, and the red curves are the synthetic waveforms.
3.2 Observational inversion tests
The source parameters of this earthquake were inverted using data from broad-band seismographs and high-rate GNSS observations. First, we perform the inversion using only the 24 high signal-to-noise ratio broad-band seismographs. The filter frequency bands of body wave and surface wave are set to 0.02–0.1 and 0.02–0.08 Hz, respectively. For comparison, we use D-CAP and P-CAP to invert for the source model using the 22 high-rate GNSS RP displacements only. Compared to the broad-band seismograph observations, the noise in high-rate GNSS RP results is significantly higher. Therefore, we also only used the R and T components of the surface waves in this inversion, with the filter frequency band set to 0.02–0.08 Hz. Finally, we conducted a joint inversion combining broad-band seismograph observations with high-rate GNSS RP waveforms. Four broad-band seismograph observations and ten high-rate GNSS RP solutions were used in this inversion to simulate sparse station coverage. The filter frequency bands were set to 0.02–0.1 Hz for body waves and 0.02–0.08 Hz for surface waves. The search interval for the fault plane parameters was set to 5°. Note that the joint inversion of broad-band seismograph observations and high-rate GNSS RP results was performed using only the P-CAP method due to its computational efficiency.
The inversion results from different methods are summarized in Table 3 and compared with the source parameters published by GCMT and USGS (United States Geological Survey). The variance reductions for all the inversions are reached 74, and they give consistent magnitude and similar fault plane parameters. The source parameters estimated by D-CAP and P-CAP are slightly different (within two step). This may be related to the iterative process of P-CAP and will be further analysed in the discussion section. In this inversion test, D-CAP took approximately 70.77 hr, while P-CAP took about 0.17 hr. Compared to the source parameters released by GCMT and USGS, the magnitudes are consistent. However, there are differences in the source depth and fault plane angles. These discrepancies could be attributed to differences in the inversion methods and data used (Kagan 2003). The fault plane parameters are different between the CAP-RP methods and CAP with seismographs. The CAP-RP results are more consistent with those from GCMT, while the CAP with seismograph results align more closely with those from USGS. This earthquake occurred on an orthogonal fault system, and this discrepancy may be related to the complexity of this fault system. Jia et al. (2020) and Liu et al. (2019) analysed this earthquake with subevent and finite fault methods respectively, and their major fault segments are shown in Table 3. Compared to other methods, the P-CAP-Joint method produces the results most consistent results with theirs. The waveform fitting of the P-CAP-Joint method is shown in Fig. 9. In the joint inversion, the broad-band seismograph observations provide far field constraints, while high-rate GNSS RP solutions provide near field constraints. This complementarity between the two data sources contributed to more reliable inversion results.
Waveform fitting of the P-CAP-Joint inversion. The first four rows are from broad-band seismographs, and the following ten rows are high-rate GNSS RP solutions. The three numbers following ‘FM’ represent the strike, dip and rake angles of the fault plane, respectively. The left side sequentially indicates the epicentral distance, station name and azimuth. The black curves are the observed waveforms, and the red curves are the synthetic waveforms.
We further discuss the convergence conditions during the iterative process for the P-CAP and P-CAP-Joint inversions. In the synthetic inversion tests, stable inversion results were obtained after two iterations. In the observational inversion tests, we set the number of iterations to 3, and the resulting variance reduction and its variation curves are shown in Fig. 10. The variance growth is the ratio of the variance reduction increase in the current iteration to the variance reduction of the previous result. Fig. 10 indicates that P-CAP achieves stable source parameters by the second iteration, which supports convergence conditions in the synthetic inversion tests. Variance reduction of P-CAP-Joint further demonstrates that broad-band seismograph and high-rate GNSS joint inversion does not affect the convergence process of P-CAP.
Variance reduction and its variation in P-CAP and P-CAP-Joint. The variance growth is the ratio of the variance reduction increase in the current iteration to the variance reduction of the previous result.
4. DISCUSSION
The CAP-RP matches observational waveforms with theoretical waveforms computed by differencing Green's functions, so as to realize the CAP inversion for high-rate GNSS RP solutions. Here, we discuss the impact of using undifferenced Green's functions for high-rate GNSS RP solutions in CAP inversion. The test was conducted with the Green's functions of each rover stations and the RP solutions, with input parameters consistent with those in the observational inversion tests. The undifferenced results are as follows: the magnitude is Mw 6.46, the focal depth is 9 km, and the strike/dip/rake angles are 50°/68°/19°, respectively. Compared to P-CAP-Joint, the magnitude is slightly larger, the focal depth shows one step difference, and the strike, dip and rake angles show a maximum difference of four steps. On the one hand, the inversion results suggest that the undifferenced approach can still yield reasonably close source parameters. On the other hand, the difference of results shows the validity of differential Green's functions when using high-rate GNSS RP for source parameter inversion. It is likely that the rapid convergence of the CAP-RP method results from the fact that the CAP method overall is relatively insensitive to errors in the Green's functions.
In both the synthetic and observational inversion tests, the iterative process of P-CAP indicates that stable source parameters are obtained after the second iteration although there is some incremental improvement beyond that. To ensure the reliability of the results, we recommend setting the number of iterations for P-CAP to 3. An additional iteration can further validate the reliability of the inversion. It is noteworthy that some minor differences are observed in the P-CAP inversion tests. We analyse the effect of the iterative convergence process and high-rate GNSS noise in these differences. Based on the synthetic inversion tests, we substitute the synthetic parameters into eq. (8) to generate the differenced Green's functions library and conduct both noise-free and noise-added inversion tests with this library. In the noise-free inversion test, the obtained source parameters are the same as the synthetic parameters. The result indicates that the iterative convergence process may introduce errors with one or two interval steps. In the noise-added inversion test, the fault strike angle still has an error within one step. This result suggests that high-rate GNSS noise has a small impact on the source parameter inversion.
In GNSS RP data processing, the baseline length between the rover and reference stations can affect the solution accuracy. A longer baseline may degrade the accuracy of solutions, which directly impacts the source parameter inversion results. To quantify this effect, we calculated the RMS of the North-East-Zenith observations (Fig. 5) 50 s before the earthquake and examined their relationship with baseline length. We fit a line to the RMS values and baseline lengths to show the trends. As shown in Fig. 11, all components exhibit positive correlations, where RMS increases with baseline lengths. Compared to the horizontal components, the Z component shows a more rapid increase, reaching centimeter level error. To ensure the precise RP solutions, the station with the smallest average distance to all other stations is selected as the reference station in this study. The initialization of synthetic coefficients in P-CAP is crucial. Based on the characteristic that |${a_3}$| is independent of the station azimuth, the synthetic coefficients are initialized to 1 for all components. The behaviours of the |${a_i}$| ratios and results from both synthetic and observational inversion tests have validated the feasibility of this initialization.
Impact of baseline length on the accuracy of RP results. FN/E/Z represents the fitted lines for N/E/Z components.
When differencing the Green's functions, there is a time-shift between the synthetic waveforms of the rover and reference stations (relative time-shift). In this study, a simplified approach is applied to handle the relative time-shift. The CAP method allows for a time-shift between synthetic and observed waveforms to account for potential errors in the velocity model. Typically, the time-shift at each station is searched independently, and the time-shifts are not the same across all stations. If the CAP-RP method considers every possible combination of time-shifts at each station, the number of Green's functions to be searched would increase significantly, which will lead to a very time-consuming inversion process. We assume that the relative time-shifts are the same across all rover stations. The advantage of this approach is that it only adds a search dimension within the time-shift range. However, the disadvantage is that the differences between different rover stations are ignored. How to better balance the difference of time-shift and the time consuming of inversion still needs further consideration.
5. CONCLUSION
We have developed a modified CAP source parameter inversion method for high-rate GNSS relative positioning records (CAP-RP). This method can use relative positioning waveforms obtained rapidly or in real time to a CAP earthquake source inversion, which meets the application requirements of high-rate GNSS in the rapid acquisition of source parameters. CAP-RP takes into account the interference in the seismic waveforms caused by the reference station's movement. The observational waveforms and synthetic waveforms are matched by differencing Green's functions. Through derivation, we establish the method for calculating the Green's functions relative to the reference station (as shown in eqs 7 and 8). The CAP-RP method is implemented using two approaches: the extended differential CAP (D-CAP) and the iterative post-processing CAP (P-CAP). In synthetic inversion tests, both methods demonstrate consistent performance. With the Green's functions library constrained by synthetic parameters from eq. (8), P-CAP is able to obtain the same source parameters as synthetic parameters. Compared to the D-CAP, P-CAP is more time-efficient, and thus is recommended. In observational inversion tests, CAP-RP is able to produce source parameters that are very consistent with those obtained from broad-band seismographs. Combining the large amplitude near field GNSS RP results with the highly sensitive far field broad-band seismograph data can yield more accurate and reliable earthquake source parameters.
We note that this paper conducts the source parameter inversions of high-rate GNSS RP records with a 1-D velocity model. However, the simplification of a 1-D velocity model may ignore the complexities of velocity variations, which could introduce potential errors, particularly in regions with complex terrain and tectonic activity. A 3-D velocity model can more accurately describe the intricate velocity variations. Future research should consider incorporating a 3-D velocity model to improve inversion accuracy.
ACKNOWLEDGMENTS
This work was supported by the National Natural Science Foundation of China (42488301, 42030311, 42374077), and the State Key Laboratory of Geodesy and Earth's Dynamics (S22L620104). Thanks to GCMT (https://www.globalcmt.org/CMTsearch.html) and the USGS (https://www.usgs.gov/programs/earthquake-hazards) for providing the earthquake source parameters.
DATA AVAILABILITY
The high-rate GNSS observations, precise satellite products and broad-band seismograph observations used in this study are from UNAVCO (https://www.unavco. org/data/gps-gnss/gps-gnss.html), CODE (http://ftp.aiub.unibe.ch/) and IRIS (https://ds.iris.edu/wilber3/find_event), respectively. The D-CAP and P-CAP source code is available on GitHub (https://github.com/zhangcf20/CAP-RP).
