Analysis of Seismic Deformation from Global Three-Decade GNSS Displacements: Implications for a Three-Dimensional Earth GNSS Velocity Field

Abstract

With the rapid development of Global Navigation Satellite System (GNSS) technology, the long-term accumulated GNSS observations of global reference stations have provided valuable data for geodesy and geodynamics studies since the 1990s. Acquiring the precise velocity of GNSS stations is very important for the study of global plate movement, crustal deformation, etc. However, the seismic activities nearby some GNSS observation stations may seriously change the station’s motion trajectory. Therefore, our research was motivated to propose a method allowing for station seismic deformation, and apply it to construct an updated global GNSS velocity field. The main contributions of this work included the following. Firstly, we improved the GNSS data processing procedures and seismic data selection strategies to obtain GNSS coordinate time series with mm-level precision (3–5 and 6–8 mm in the horizontal and vertical, respectively) and information of each site impacted by seismic events, which provides necessary input data for further analysis. Secondly, an Integrated Time Series Method (ITSM) concerning the effect of seismic deformation was proposed to model the station’s nonlinear motion accurately. Distinguished with existing studies, all parameters including seismic relaxation time can be simultaneously estimated by ITSM, which improves the accuracy and reliability of GNSS station velocity significantly. Thirdly, to optimize the ITSM-based model, the influences of seismic relaxation time (a. 0.1 × true, b. 10 × true, c. true), parameterization mode (a. Offset + Velocity, b. Offset + Velocity + PSD, c. Offset + Velocity + PSD + Period), and the Post-Seismic Deformation (PSD) model (a. None, b. Exp, c. Log, d. Exp + Log) on results of GNSS time series analyzing were discussed. The results showed that the fitting accuracy of GNSS displacements was better than 5 mm and 10 mm in the horizontal and vertical, respectively. Finally, the global GNSS station velocity field (referred to as GGV2020 hereafter) was refined by ITSM using global GNSS observations and seismic data during 1990–2020. This not only helps interpret plate tectonic motion, establish and maintain a Dynamic Terrestrial Reference Frame (DTRF) but also contributes to better investigating geodynamic processes. GGV2020 results showed that the accuracy of global velocity was better than 1 mm/a, and the averages of Root Mean Square Error (RMSE) were 0.19 mm/a, 0.19 mm/a, and 0.33 mm/a in the north, east, and up direction, respectively. Besides, the RMSE obeys normal distribution. Compared with ITRF2014, there was a difference of about 1–2 mm/a between them due to differences in terms of observation span, processing model, and geodetic technology. Moreover, GGV2020 is expected to enrich and update the existing velocity field products to describe the characteristics of regional crustal movement in more detail, especially in Antarctica.

1. Introduction

In 1994, Strange et al. first proposed the concept of Continuous Operation of Reference Stations (CORS) [1]. Some countries, regions, and relevant organizations in the world have successively established high-precision GNSS CORS networks. Among them, IGS CORS is the most widely distributed and the largest continuous station network [2,3]. With the enrichment of GNSS observations and the improvement of its processing strategies, the long-term accumulated GNSS observations of global reference stations have provided basic data for geodesy and geodynamics studies [4,5,6]. In these related studies, the theory and application of GNSS coordinate time series analysis have become a research hotspot in the field of geodesy and geophysics. GNSS station coordinate time series usually contain both linear and nonlinear variation signals. The former, generally expressed as velocity, describes the long-term tectonic movement trend in a regional area, while the latter, particularly the periodical variation, mainly reflects the geophysical effects [7,8,9]. Additionally, there are offsets caused by equipment replacement and other reasons [10,11] or co-seismic and post-seismic deformation impacted by the earthquakes [12,13,14] in the time series. Using the long-term GNSS coordinate time series to study the various mechanisms of nonlinear motion, we can accurately separate the linear and nonlinear motion of the station. It not only helps explain the plate tectonic movement reasonably, establish and maintain the Dynamic Earth Reference Frame (DTRF), but also better study the geodynamic processes such as post-ice rebound, sea-level change, and inversion of ice-snow mass change.

In the past 30 years, many violent earthquakes have occurred in the world. The GNSS sites within the range of hundreds or thousands of meters around the epicenter have been affected, and the seismic deformation will last for months to years [15,16]. Therefore, it is necessary to extract and model the seismic deformation in GNSS coordinate time series. Large earthquakes cause not only significant Co-Seismic Deformation (CSD) in the surrounding area but also time-varying Post-Seismic Deformation (PSD) through the transmission of surrounding geological media [17]. PSD mainly includes short-term deformation caused by the rebound of pore fluid [18], fault slip [19,20], and long-term deformation caused by viscous relaxation [21]. In the PSD model of long-scale GNSS coordinate time series, the short-term deformation of several hours to several days can be ignored. The long-term deformation of several months to several years is often represented as logarithmic or exponential relaxation terms [22]. The extraction of seismic deformation in GNSS time series, on the one hand, plays a vital role in studying the physical mechanism of the seismic fault and crustal rheology. On the other hand, more accurate station velocity can be estimated by eliminating the influence of seismic deformation, which is the basis for studying plate motion, stress-strain accumulation, earthquake preparation, and prediction. Besides, modelling seismic deformation is also crucial to improve the precision of global DTRF realization [23,24,25]. Therefore, it is indispensable to select the sites affected by the earthquakes and establish a perfect time series model including seismic deformation.

There are three elements in the PSD model: earthquake occurrence time, amplitude coefficient, and relaxation time. For the time of earthquake occurrence, some scholars have proposed automatic detection methods of offset and transient deformation [10,26,27,28]. Nevertheless, the automatic detection of offset and transient deformation without a specific epoch is challenging, which also has a high missing rate and misjudged rate. Comparatively, the results on manual detecting offsets in GNSS time series could be better [24]. The amplitude coefficient can be estimated as an unknown parameter in the GNSS time series model. For relaxation time, it usually needs to be estimated separately. Some studies directly set it to one year [29], and assumed that all stations in a nearby grid area obey the same decay [30], or estimated it through the maximum likelihood method [31]. Some scholars have also proposed the method of trial and error to detect the CSD and PSD of the GNSS station [32]. Li et al. adopted a nonlinear least square model to estimate post-seismic logarithmic relaxation time in the GNSS time series [33]. In the above studies, the relaxation time was estimated individually and then substituted back to the model, where the correlations between the relaxation time and other parameters were disregarded and the implementation is more complicated.

For the above reasons, we propose an Integrated Time Series Method (ITSM) including seismic deformation to realize a comprehensive solution for GNSS time series analysis, which improves the accuracy of parameter estimation and the fitting result of time series effectively. Section 2 briefly presents the global data sets during 1990–2020 and processing strategies used in this paper to obtain high-precision GNSS coordinate time series and information of sites impacted by seismic events. Then, an integrated time series analysis method adding seismic deformation is introduced to accurately describe the station’s nonlinear motion trajectory, which is an ordered combination of the Basic Time Series Model (BTSM), Iteratively Reweighted Least Squares (IRLS), and Constrained Nonlinear Optimization (CNO). In Section 3, we optimize the model in detail, analyze the fitting accuracy of the time series model with different combination modes, and verify the reliability of the ITSM using case results. Then, we construct and evaluate GGV2020 using global GNSS observations and seismic datasets during the past three decades, and make a detailed comparative analysis of ITRF2014 and GGV2020 to validate the accuracy and reliability of the latter. In Section 4, we further discuss the impact of velocity discontinuity on ITSM results. In Section 5, we formulate some valuable conclusions.

2. Data and Methods

This paper aimed to study the influence of major earthquakes on long-term global GNSS coordinate time series. For this purpose, we collected the GNSS raw observations of global CORS (data source: https://garner.ucsd.edu/archive/garner/rinex/, the accessed date on 1 January 2021) and the complete records of the global seismic event (Mw ≥ 5.0) (data source: https://earthquake.usgs.gov/, the accessed date on 1 January 2021) from 1990 to 2020. Then, we processed the collected GNSS observations and seismic records to obtain high-precision GNSS coordinate time series as a clean and complete input for the subsequent analysis. In addition, we proposed an integrated time series analysis method adding seismic deformation to accurately model the station’s nonlinear motion trajectory in this section, and the ITSM is an ordered combination of BTSM, IRLS, and CNO.

2.1. GNSS Observations and Seismic Records

The GNSS station observation data recorded for more than 1000 global CORS sites (Figure 1) in this paper are mainly from the Scripps Orbit and Permanent Array Center (SOPAC).

Further, the longer the observation span of GNSS stations is, the more reliable results from time series analysis and modeling are. Because GNSS time series usually contain annual and semi-annual signals, it is required that the time series should span more than one year [34]. Considered the effects of data defect and other factors (offset, earthquake, noise, etc.) on results of the GNSS time series analysis, only those stations with observations of more than two years were selected for the following experiment. Table 1 gives the statistics results of global GNSS observations’ duration. Among them, the length of the shortest time series is less than one year, and the longest is more than 28 years. The sites with observations of less than 2 years accounted for about 20%, and those with observations of 2–10 years and more than 10 years were both about 40%. The average length of time series of all GNSS stations was more than 10 years.

For seismic activities, we collected the records (the location of the epicenter, magnitude, occurrence time, etc.) of global seismic events (Mw ≥ 5.0) in the past 30 years. The distribution of all recorded seismic activities is displayed in Figure 2. The distribution of these seismic events presents the main seismic belts (the Pacific, Eurasian, and ridge seismic belts) in the world [35,36]. Seismic zones are mainly at the junction of global plates, and the distribution of earthquakes is concentrated in the seismic zone and scattered outside the seismic zone.

There occurred 51,689 significant earthquakes between 1990 and 2020; thereof, there were 51,222 with Mw 5–8, 437 with Mw 7–8, and 30 with Mw 8–10. It is concluded that the global crustal activities have indeed been intense in the last 30 years, and especially the GNSS sites around the significant earthquakes will be affected, resulting in seismic deformation. If the structural deformation is not considered, the accuracy and reliability of the time series model will be seriously reduced.

2.2. Processing Strategies for GNSS and Seismic Data

At present, with the improvement of GNSS observation technology and the data processing model, the accuracy of the GNSS daily coordinate solution could respectively achieve 3–5 mm and 6–8 mm in the horizontal and vertical direction [3,37], which made it possible for good modeling of seismic deformation in long-term GNSS time series. In the field of GNSS data processing, GAMIT/GLOBK [38], GIPSY/OASIS [39], and Bernese [40] have been widely used to obtain high-precision daily solutions of GNSS station coordinates. Compared with GIPSY/OASIS packages and Bernese software, the GAMIT/GLOBK software is open source, and has the following advantages: high precision calculation, fast calculation efficiency, timely updates, and high automation. Additionally, the relative accuracy of the long baseline solution of the latter can reach about 10⁻⁹ [38]. Figure 3 shows the flow of high-precision GNSS data processing using GAMIT/GLOBK in our works.

The more detailed GNSS data processing steps are as follows:

• GNSS data preparation and preprocessing: We need to download navigation ephemeris (BRDC) and precise ephemeris (SP3) corresponding to observation files (RINEX) in advance, and the tables file used for GAMIT/GLOBK, and preprocess RINEX (Receiver Independent Exchange Format) files by TEQC [41] or GFZRNX tools [42];

• Stations subnet division: To improve data processing efficiency, we divide the global sites into multiple subnets. Any subnet should have 3–6 common sites with its immediate neighboring subnet;

• Baseline processing: Daily GNSS baseline solutions for all subnets are carried out by GAMIT/GLOBK software with a distributed processing approach [38];

• Subnet fusion: The common sites between subnets are used to merge all subnets into a global network by GLRED [38] software;

• Comprehensive adjustment: We will finally obtain the high-precision GNSS coordinate time series of global sites using GLRED software.

As for the seismic records, although we obtained all global significant earthquake records in the past 30 years, it is impossible to perform a model and analysis of all seismic events. One reason is the great number of earthquakes, and the other is that it is unnecessary to model for every earthquake because not every earthquake will cause the GNSS station displacements. Therefore, we selected the earthquakes and sites impacted by seismic deformation, and the others were ignored. Figure 4 shows the selection flow of the earthquakes and sites impacted by seismic deformation.

More steps are as follows:

• Input the seismic records, including the location of the epicenter, magnitude, epoch, etc.;

• Search for all GNSS sites within a specific range of the epicenter (difference in longitude ≤ 3 × Mv, the difference in latitude ≤ 2 × Mv);

• Extract the two subsequences called sub1 and sub2 at the epoch of a seismic event. Then, find the median difference between sub1 and sub2. If the difference is greater than the default threshold, mark the earthquake and site;

• Repeat step 3 until the sites of step 2 are traversed;

• Repeat steps 1–4 until all earthquakes are traversed. We finally obtain all earthquakes and sites impacted by seismic deformation.

2.3. An Integrated Time Series Method for GNSS Displacements

A GNSS coordinate time series usually contains the secular trend (linear velocity), seasonal variation (annual and semi-annual signals), offsets caused by non-seismic factors (equipment replacement, antenna height measurement error, phase center modeling error, or other human and software errors) or seismic factors (co-seismic deformation), post-seismic deformation, and other unmodeled errors. Without consideration of the unmodeled error term, the BTSM of GNSS stations can be composed of four sub-models:

(1)$\mathrm{X}\left(t\right)=X_{velocity}+X_{period}+X_{offset}+X_{psd}$

Because of the weak correlation between the coordinate components of different epochs [43], an individual component time series (∆N, ∆E, or ∆U) at discrete epochs $t_i$ can be modeled as [44]:

(2)$$\begin{gathered} y\left(t_i\right)=a+b{t}_i+csin\left(2\pi t_i\right)+dcos\left(2\pi t_i\right)+esin\left(4\pi t_i\right)+fcos\left(4\pi t_i\right) \\ +{\displaystyle \sum }_{j=1}^{n_g}{g}_{j}H\left(t_i-T_{g_j}\right)+X_{psd}+{\epsilon }_i \end{gathered}$$

where $a$ is the value at the initial epoch $t_0$, $t_i$ denotes the time elapsed from $t_0$, the linear velocity $b$ represents the secular tectonic motion, and the coefficients c, d, e, and f denote annual and semi-annual variations in GNSS coordinate time series. Annual and semi-annual terms can be estimated accurately when enough GNSS Observations have been collected. Amplitude and phase of seasonal signals are expressed according to the sine function. The magnitudes g of $n_g$ offsets are due to co-seismic deformation or non-coseismic changes at epochs $T_{g_j}$. H denotes the discrete Heaviside function. $X_{psd}$ is the function of PSD, and ${\mathsf{\epsilon }}_i$ is the observation noise.

Where H is expressed as:

(3)$H=\{\begin{array}{c}0, t_i-T_{k_j}<0\\ 1, t_i-T_{k_j}\ge 0\end{array}$

$X_{psd}$ can usually be expressed as an exponential or logarithmic function. The logarithmic model is another parameterization associated with afterslip on the fault surface; the exponential model is associated with motion below the crust [20]. The exponential model is expressed as [45]:

(4)${\displaystyle \sum }_{j=1}^{n_k}{k}_j\left(1-e^{\left[-\left(\frac{t_i-T_{k_j}}{{\tau }_j}\right)\right]}\right)H\left(t_i-T_{k_j}\right)$

Coefficients $k_j$ are for $n_k$ PSD starting at epochs $T_{k_j}$ and decay exponentially with a relaxation time ${\tau }_j$. The logarithmic model is expressed as [21]:

(5)${\displaystyle \sum }_{j=1}^{n_k}{k}_{j}log\left(1+\frac{t_i-T_{k_j}}{{\tau }_j}\right)H\left(t_i-T_{k_j}\right)$

2.3.1. The Estimation of Seismic Relaxation Time Factor

The existence of the seismic relaxation time factor in the PSD model causes the nonlinear problem of the objective function. It should be noted that there is a difference between linear programming and nonlinear programming. If the optimal solution exists, the optimal solution of the former can only be achieved at the boundary of the feasible region, while the optimal solution of the latter can be performed at any point of the feasible region.

It can be seen from Figure 5 that the nonlinear optimal solution is sensitive to the initial value and interval of the parameter to be estimated. Therefore, giving appropriate parameter constraints, we can use nonlinear programming to obtain the parameters of interest.

In this paper, we adopted the CNO model taking the seismic relaxation time as the nonlinear objective function parameter to obtain the optimal global solution of the objective function using an iterative algorithm. Considered the velocity, period, offset, and PSD model, the CNO model of GNSS time series is expressed as:

(6)$$\begin{gathered} \underset{x}{min}\Vert f{\left(x\right)\Vert }_2^2=\underset{x}{min}\left(f_{velocity}{\left(x\right)}^2+f_{period}{\left(x\right)}^2+f_{offset}{\left(x\right)}^2+f_{psd}{\left(x\right)}^2\right) \\ s. t. \{\begin{array}{c}{x}_0=Initial value\\ \begin{array}{c}x \ge Lower bounds\\ x \le Upper bounds\end{array}\end{array} \end{gathered}$$

The Levenberg-Marquardt and trust-region-reflective algorithms are based on nonlinear programming. We used the trust-region-reflective method, a subspace trust-region method, based on the interior-reflective Newton method [46,47]. Each iteration involves the approximate solution of an extensive linear system using preconditioned conjugate gradients. The Levenberg-Marquardt method was described in references [48,49,50]. It should be noted that the trust-region-reflective algorithm cannot solve uncertain systems, and it requires that the number of equations (the row dimension of function) be at least as great as the number of unknown parameters. In uncertain cases, we can use the Levenberg-Marquardt algorithm.

2.3.2. The Estimation of Other Time Series Parameters

The occurrence time of each event at GNSS station can be available from the earthquake catalog, observation log, automatic detection algorithm, or visual inspection. Because the automatic detection algorithm is unreliable and the visual inspection method is complicated, we determined the offset epoch and the PSD epoch in the GNSS time series mainly according to the seismic catalog and the station log. The seismic relaxation time is typically estimated separately by other methods (see Section 2.3.1) so that estimation of the remaining GNSS time series parameters can be expressed as a linear problem:

(7)$\{\begin{array}{c}\mathit{y}=\mathit{B}\widehat{\mathit{X}}+\mathit{\epsilon }\\ {\mathit{E}}_{\epsilon }=0\\ {\mathit{D}}_{\epsilon }={\sigma }_0^{2}c\end{array}$

where $\mathit{B}$ is the design matrix and $\widehat{\mathit{X}}$ is the parameter vector to be estimated, $\widehat{\mathit{X}}={\left[a b c d e f \mathit{g} \mathit{k}\right]}^T$, ${\mathit{E}}_{\epsilon }$ is the statistical expectation, ${\mathit{D}}_{\epsilon }$ is the statistical dispersion, ${\mathit{C}}_{\epsilon }$ is the covariance matrix of observation errors, $\mathit{P}={\mathit{C}}_{\epsilon }{}^{-1}$ is the weight matrix, and ${\mathsf{\sigma }}_0^2$ is an a priori variance factor.

The error equation is:

(8)$\mathit{V}=\mathit{B}\widehat{\mathit{X}}-\mathit{L}$

According to the least square fitting estimation, the first estimation is obtained under the estimation criterion of Equation (9).

(9)${\mathit{V}}^{\mathit{T}}\mathit{P}\mathit{V}=\min$

(10)$\widehat{\mathit{X}}={\left({\mathit{B}}^{\mathit{T}}\mathit{P}\mathit{B}\right)}^{-1}\left({\mathit{B}}^{\mathit{T}}\mathit{P}\mathit{L}\right)$

The weight of each observation is re-determined according to the residuals $\mathit{V}$ and the weight factor Equation (11), and the iterative calculation ends until the difference between the two last estimated values is less than the tolerance.

(11)$\rho \left(v\right)=\{\begin{array}{c}1, v\ge C\\ \frac{1}{\left|v\right|}, v<C\end{array} C=3{\sigma }_0$

The above three models (BTSM, IRLS, and CNO) are introduced respectively, while the three are used together in the practical application process. ITSM describes the combination and optimization of the three models to obtain the optimal solution of the objective function parameter estimation. Figure 6 shows the analysis flow of ITSM.

More steps are as follows:

• Set the initial value of seismic relaxation time (tau), with the default value of tau being one year;

• Perform the first IRLS solver and obtain the initial values of other parameters;

• The estimated parameters in step 2 are used as the initial values of the parameters of the CNO model. Set appropriate parameter constraints;

• Perform the CNO solver and obtain the tau’s optimal solution;

• Perform the second IRLS solver and obtain the time series parameter’s optimal solution;

• Perform the BTSM solver for accuracy evaluation and parameter analysis.

3. Results and Analysis

In this section, we mainly carried out the comparison and evaluation of time series analysis results and verify the accuracy and reliability of the ITSM. Finally, we constructed and evaluated GGV2020 using global GNSS observations and seismic datasets of three decades.

3.1. Analysis of GNSS Time Series Preprocessing

The purpose of time series preprocessing is to obtain a clean and reliable time series. Generally, two steps are required to eliminate large outliers and mark the values with low precision. Besides, to simplify the GNSS time series analysis steps, it is necessary to mark the offset epoch and seismic epoch in advance. The seismic events and sites impacted by seismic deformation were selected using the strategy in Section 2.2. Figure 7 illustrates in red the location of 208 earthquake epicenters that caused significant seismic deformation, green denotes the 365 impacted stations located at these global GNSS sites. The results provided the basic materials for the construction of seismic deformation models in the following research. Additionally, the offsets in GNSS displacements, caused by equipment replacement, were automatically detected with the station log.

The default criteria in Table 2 were applied to detect and reject outliers in the GNSS time series. However, we can also specify individual criteria for each station if needed:

• Weak observation criteria based on the formal errors. If the formal sigma of one site is more significant than the criteria at one epoch, the solution of this site at this epoch will be ignored;

• Outliers criteria based on the post-fit residuals. If the residuals of one site are bigger than the criteria at one epoch, the solution of this site at this epoch will be ignored;

• Bad observation criteria based on the outlier threshold. If the values have outliers, the initial adjustment will be biased. They will be removed before the adjustment.

Due to a large number of global GNSS sites, we only took the ANTC site as a case study. The result of preprocessing of ANTC was shown in Figure 8. The bad values, weak values, offsets, and seismic events are respectively marked in the north, east, and up direction. The offsets and seismic events are synchronously marked, but not for bad and weak values.

According to the station logs, there are 10,832 offsets of more than 1000 sites involved in this study, which were caused by equipment replacement. Therefore, we must model the offsets in the GNSS time series. Otherwise, we may get terrible results; we paid particular attention to the velocity results. In addition, it should be noted that some equipment replacements may only cause minor offsets or even no offsets. For these minor offsets caused by equipment replacement, if the offset value was less than 3 mm, we judged it as invalid, and it was not treated as an offset in the subsequent analysis.

3.2. Feasibility Analysis of Relaxation Time Selection

It has been mentioned in Section 2.3 that the appropriate settings of initial value constraint conditions of the parameters have positive significance for nonlinear programming. If the initial value of tau is set to zero in this experiment, the parameter solution cannot be obtained by nonlinear programming. Bevis et al. mentioned that the tau value is insensitive to the time series model. We further analyzed how the selection of tau will affect the time series model under different model conditions to provide a good foundation and explanation for the selection of initial value in nonlinear programming.

Figure 9 shows the influence of tau on PSD with different combination models. In this experiment, the PSD was considered for all combinations, and the first model only included the offset. The second included offset and velocity, and the third included offset velocity and period. The results showed that: RMS (tau = true × 1.0) < RMS (tau = true × 10) < RMS (tau = true × 0.1) of the three models. With the improvement of the model, the fitting accuracy was higher no matter what the value of tau was: RMS (PSD + Offset + Velocity + Period) < RMS (PSD + Offset + Velocity) < RMS (PSD + Offset). It can be found that the RMS (PSD + Offset, tau = true × 0.1) was the largest, and the accuracy was the worst among the nine combinations, which was more than 10 mm. The RMS difference of the first model using three taus was more than 5 mm, that of the second one was no more than 1 mm, and that of the third was the smallest and no more than 0.5 mm.

From the above analysis, it can be concluded that the value of seismic relaxation time has no significant influence on the overall trajectory and variation characteristics of the single GNSS coordinate time series. After the relaxation time, offsets, velocity, periodic signals, and other coefficients in the time series model were estimated, the station displacements can then be accurately represented, and the fitting accuracy can reach 10 mm.

In addition, to illustrate the validity of setting the tau’s initial value to 1.0 in the CNO model, we analyzed all the GNSS time series impacted by seismic deformation, using ITSM with the tau’s initial value of 1.0. Table 3 provides the fitting RMS of the GNSS sites impacted by seismic deformation. The mean RMS of these stations was about 2–3 mm and 6–7 mm in the horizontal and vertical direction, respectively. Setting the tau’s initial value to 1.0 in the CNO model is appropriate.

3.3. Analysis of Case Result on Different Model

In the ITSM model, various factors are considered comprehensively. However, to simplify the calculation and express the parameter results more concisely, some parameter items with less influence may be ignored. Hence, we analyzed the impact of each parameter on the time series fitting accuracy, mainly for PSD and periodic terms. Figure 10 shows the results of time series fitting models.

In this experiment, the first model included the offset and velocity, the second included offset, velocity, and PSD, and the third included offset, velocity, PSD, and periodic signals. The results of fit1 (red curve) deviated from the actual trajectory the most. The differences between the fitting velocities and the actual values were 3, 6, and 4 mm/a in three directions and the RMS of time series residuals were more than 19, 87, and 29 mm, respectively. The second model’s results (blue curve) were consistent with the actual trajectory. The results of fit2 and fit3 showed that the RMS of residuals in the horizontal direction was the same, while the RMS difference in the vertical direction was about 1 mm, which shows that the periodic variation in the horizontal direction of the site is not apparent, while that in the vertical is apparent. The velocity solutions of fit2 and fit3 were the same, indicating that the accuracy of the velocity solution is hardly affected by the periodic term. If the parameter of interest is the velocity, the periodic term may be ignored. It can be seen from Figure 10 that the site’s trajectory in the U direction presented a noticeable periodic movement.

We also analyzed other stations, and most of them followed the above analysis results:

• The accuracy of the velocity solution was hardly affected by the periodic term because the length of the time series we used was more than two years. However, if the duration is less than one year, the accuracy of velocity (especially vertical velocity solution) may decrease;

• The periodic variation in the horizontal direction was not obvious, while the periodic variation in the vertical direction was noticeable. This is because the periodic change caused by the change of surface load which is mainly reflected in the vertical direction.

The next question is how to determine whether the sites existed PSD and the specific expression form of PSD (i.e., exponential or logarithmic form). In case of a time series with a unique earthquake causing PSD, four different models were tried: None (PSD1), Exp (PSD2), Log (PSD3), and Log + Exp (PSD4), each combined with the offsets, velocity, and periodic terms. Among all the models, we finally selected the model with the best fitting accuracy.

As shown in Figure 11, we used four PSD models for analysis, and we found that RMS (PSD1) > RMS (PSD2) > RMS (PSD3) > RMS (PSD4). The fitting results of PSD4 (green curve) had the best accuracy, whose RMS of the three-dimensional residuals were about 4.7, 12.4, and 7.2 mm, respectively. The results of PSD1 (red curve) deviated from the actual trajectory the most. The difference between PSD3 and PSD4 was minimal, and the difference of velocity solution was less than 1 mm/a, and the difference of residual RMS was less than 1 mm, which indicates that the PSD model is mainly explained by logarithmic form.

We also analyzed other stations. The results showed that some sites had no obvious PSD after the earthquake. Among the sites impacted by PSD, most of the PSD models are explained by logarithmic changes, and some are explained by exponential changes. Additionally, a few ones based on the logarithmic and exponential combination models achieve the best fitting accuracy.

So far, we have solved a series of problems encountered in the GNSS time series and began to use the ITSM to analyze the global GNSS time series. Figure 12 shows the fitting results of ANTC’s time series using the ITSM.

It can be seen from Figure 12 that modeling the offset, velocity, PSD, and seasonal variation (annual and semi-annual) in GNSS time series, and the ITSM model can almost perfectly describe the trajectory of the site in each direction. Table 4 shows the specific parameter result of ANTC’s time series using the ITSM.

To further verify the accuracy and reliability of the ITSM, in addition to the ANTC site, the time series fitting results of some sites distributed worldwide with different change trajectories are listed in Figure 13. It can be seen from Figure 13 that the fitting results of these sites were good, which verified the correctness and reliability of the ITSM. The slop accuracy of these sites was better than 1 mm/a, which provides data support and a theoretical basis for analyzing the global GNSS velocity field.

To further study whether the different PSD combination patterns have certain spatial distribution characteristics, Figure 14 shows the PSD combination patterns of the GNSS sites that impacted seismic deformation. It illustrates the sites of PSD expressed with Exp in green, and that of the PSD expressed with Log in blue, and that of the PSD expressed with Exp and Log in red. Unfortunately, it was not found whether the PSD combination patterns had obvious distribution characteristics, and the exponential and logarithmic model’s distribution was random. Among all these sites’ impacted seismic deformation, the PSD models of most sites were composed of Exp and Log. The specific combination mode of the PSD is explained by the characteristics of each earthquake and its surrounding geological environment. If we want to further explore the hidden features, it is necessary to analyze individual earthquakes in future research.

Then, we performed the ITSM model on global GNSS coordinate time series of more than 1000 sites. To evaluate the accuracy and reliability of the ITSM model, Figure 15 shows the fitting RMS of all GNSS time series using ITSM. Among all GNSS time series fitting results, the maximum RMS was less than 10, 10, and 15 mm in the north, east, and up direction, respectively. The mean RMS was 2.21, 2.34 and 6.33 mm in the north, east and up direction, respectively. On the whole, the fitting accuracy in the horizontal direction was better than that in the vertical. Most sites’ fitting accuracy was better than 5 and 10 mm in the horizontal and vertical direction, respectively.

In summary, the ITSM proposed in this study is suitable for time series analysis of global GNSS sites, and we can acquire the coordinate fitting accuracy with a horizontal direction better than 5 mm and a vertical direction better than 10 mm, and velocity fitting accuracy with a three-dimensional direction better than 1 mm/a.

3.4. Comparison and Analysis of the New Global GNSS Velocity Field

In Section 3.3, we performed the ITSM model on all global GNSS coordinate time series, and then we extracted the global high-precision GNSS velocity field products from 1900 to 2020 (GGV2020). To evaluate the accuracy of GGV2020, Figure 16 showed the RMS of global GNSS sites’ velocity field in the north, east, and up direction. The medians of GGV2020 RMS were 0.19, 0.19, and 0.33 mm/a in each direction. The vertical velocity accuracy of GGV2020 was lower than the horizontal because the accuracy of the former coordinate solution was lower than that of the latter. We can demonstrate GGV2020’s RMS distribution characteristics using the half-normal distribution shown with the gray curve in Figure 16.

ITRF2014 is the latest version of the international terrestrial reference frame, which was released in 2014. However, as the years went by, the global crustal movement is also changing, so the accuracy of ITRF products has declined, but it can still keep the dynamic cm-level. To verify the reliability and accuracy of GGV2020, we compared it with ITRF2014. Figure 17a,b represent the global GNSS three-dimensional velocities of ITRF2014, Figure 17c,d the GGV2020’s global GNSS velocities, and Figure 17e,f directly show the comparison of both velocity fields. In those figures, the red arrow represents the GGV2020’s velocities and the blue one represents the ITRF14’s velocities. It can be seen from Figure 17 that the two products differ in the number of GNSS sites, and the two products show a similar crustal movement trend. As GGV2020 includes more GNSS sites than ITRF2014, it is expected that a more robust global plate motion model can be derived, compared with the ITRF2014 one, involving more sites on stable areas of the tectonic plates.

GGV2020 not only enriches and updates the existing velocity field products but also is expected to have more reliable accuracy through the analysis of seismic deformation from global three-decade GNSS displacements. GGV2020 horizontal velocities illustrated by Figure 17 show the clear regional movement trend. It represents a clockwise rotation trend in North America and South America, while in Europe, Asia, and Africa there exist the opposite characteristics, that is, a counter-clockwise rotation trend. Due to its geographical location, there have unique characteristics of crustal movement in Oceania. Compared with ITRF2014, GGV2020 reveals the characteristics of clockwise rotation in Antarctica through more GNSS sites. In addition, GGV2020 vertical velocities illustrated by Figure 17 also represent the clear regional patterns, especially in Greenland, North America, and Fennoscandia, reflecting the uplift caused by glacial isostatic adjustment and recent ice melting. In the Antarctic region, there is no obvious feature of decline or uplift.

Figure 18 shows the difference in the same sites’ velocity fields between ITRF2014 and GGV2020. The velocity differences of the two products in three directions were less than 2 mm/a, and the RMS of differences were 1.73, 1.18, and 1.83 mm/a, respectively. The difference between the velocity fields obeys normal distribution. The above analysis results showed that the two products are similar, but there is still a difference of about 1–2 mm/a on the whole. The specific reasons are as follows:

• The observation’s duration of the GNSS site is different. ITRF14 uses GNSS observation data util 2014, while GGV2020 uses all data from 1990 to 2020, so the latter has a larger amount of data, and the dataset is updated;

• The solution models of time series analysis are different. Compared with ITRF2014, GGV2020 is obtained by using the ITSM model, which considers as many factors as possible in GNSS coordinate time series, especially considering the nearby major seismic activities;

• The geodetic techniques are different. In addition to GNSS technology, ITRF2014 also integrates VLBI, SLR, and DORIS technologies. Therefore, the results of the two will inevitably have a minor difference.

Through the above comparative analysis of both global GNSS velocity fields, we can see that the two global GNSS velocity fields difference in the three-dimensional direction was minor, about 1–2 mm/a, and they represent the similar characteristics of regional crustal movement. The GGV2020 is expected to enrich and update the existing velocity field products to describe the characteristics of regional crustal movement in more detail. Compared with ITRF2014, GGV2020 is based on more GNSS sites, updated observations, and the ITSM model. In addition, GGV2020 represents the characteristics of crustal movement in Antarctica.

4. Discussion

In the BTSM model, we regard the station’s velocity as a constant. However, it can be seen from Figure 12 that the fitting accuracy of the ANTC site was low up to 12 mm in the east direction, but the normal one should be 3–5 mm.

Figure 19 shows the residuals of ANTC using the ITSM, and we can find that the slope of the residuals was discontinuous, especially in the east direction. In the same series, as time changed, the site’s velocity also changed slightly, which is called the velocity discontinuity. To solve this problem, we added a sub-model to optimize the model further:

(12)$$\begin{gathered} y\left(t_i\right)=a+b{t}_i+csin\left(2\pi t_i\right)+dcos\left(2\pi t_i\right)+esin\left(4\pi t_i\right)+fcos\left(4\pi t_i\right) \\ +{\displaystyle \sum }_{j=1}^{n_g}{g}_{j}H\left(t_i-T_{g_j}\right)+{\displaystyle \sum }_{j=1}^{n_h}{h}_{j}H\left(t_i-T_{h_j}\right)t_i+X_{psd}+{\mathsf{\epsilon }}_i \end{gathered}$$

The possible $n_h$ changes in velocity are denoted by the new velocity values $h_j$ at epochs of $T_{h_j}$. We re-performed the new ITSM on the ANTC time series to acquire its residual results, as shown in Figure 20. After adopting the new model, the fitting accuracy of the time series in the three directions was improved, which are 3.15, 2.90, and 6.44 mm, respectively, with an increase of about 14%, 35%, and 7% compared with the previous model.

We also analyzed other stations. Most of the stations had no apparent discontinuity in speed. Some of the stations that had very little difference are also regarded as the same. A small number of stations had a noticeable discontinuity in velocity, and the velocity difference was less than 1 mm/a.

5. Conclusions

This research mainly contributed to an integrated time series analysis model adding station seismic deformation for global GNSS observations and seismic activities over three decades period, and we applied the ITSM to construct a new global GNSS velocity field.

The ITSM model directly estimates seismic relaxation time by using the CNO model to overcome the problem of separating relaxation time from other parameters, which can significantly improve the accuracy and reliability of velocity parameters. Experimental results showed that the ITSM proposed in this study is suitable for time series analysis of global GNSS sites, and we can achieve a coordinate fitting accuracy better than 5 and 10 mm in the horizontal and vertical, and velocity fitting accuracy better than 1 mm/a in each direction. In the process of model optimization, we also formulated some practical conclusions. (1) The initial value of seismic relaxation time has no significant influence on the overall trajectory and variation characteristics of the single GNSS coordinate time series; (2) The accuracy of horizontal velocity solution is hardly affected by the periodic term; (3) The exponential and logarithmic combination model mainly explains the PSD of most GNSS stations, while the others are explained by logarithmic or exponential models; additionally, the combination pattern of the PSD model presents the characteristics of random distribution in geographical space.

The GGV2020 production of this research showed that the accuracy of global velocity was better than 1 mm/a, and the average of RMSE were 0.19, 0.19, and 0.33 mm/a in the north, east, and up direction, respectively, and the RMSE obeyed half-normal distribution. Compared with ITRF2014, there was a difference of about 1–2 mm/a between them due to the difference in data duration, processing model, and geodetic technology. Although both products represent similar characteristics of regional crustal movement, GGV2020 is expected to enrich and update the existing velocity field products to describe the characteristics of regional crustal movement in more detail, especially in Antarctica.

The ITSM model based on constrained nonlinear programming to estimate seismic relaxation time is the core theoretical achievement of this research, and the global GNSS site velocity field obtained from this model is a crucial product. Since this research used the global GNSS observations and seismic data from 1990 to 2020 as the original data input, with a large amount of observation data and an updated dataset, the research results of this paper are anticipated to provide a reference for GNSS time series analysis, and the GGV2020 is expected to enrich and update the existing velocity field products.

Author Contributions

Conceptualization, Y.R. and L.L.; methodology, Y.R.; software, Y.R.; validation, L.L. and J.W.; formal analysis, Y.R.; investigation, Y.R.; resources, Y.R. and L.L.; data curation, Y.R.; writing—original draft preparation, Y.R.; writing—review and editing, L.L., J.W.; visualization, Y.R.; supervision, L.L. and J.W.; project administration, L.L. and J.W.; funding acquisition, L.L. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Youth Science Fund Project of National Natural Science Foundation of China, grant number 11903065.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

The authors gratefully acknowledge the SOPAC (Scripps Orbit and Permanent Array Center) for providing GNSS raw observations and USGS (United States Geological Survey) for Seismic records. The authors express thanks to MIT (Massachusetts Institute of Technology) and SOPAC for providing high precision GNSS data post-processing software (GAMIT/GLOBK). And some of the figures were plotted by the GMT (Generic Mapping Tools) graphics package [51].

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Distribution of global GNSS CORS network.

Figure 2. Distribution of global seismic events (Mw ≥ 5.0).

Figure 3. High-precision GNSS data processing flow based on the GAMIT/GLOBK software.

Figure 4. Selection flow of earthquakes and sites impacted by seismic deformation.

Figure 5. Linear and nonlinear function graphs.

Figure 6. Analysis flow of the integrated time series model.

Figure 7. Distribution of earthquake epicenters and sites impacted by seismic deformation.

Figure 8. Preprocessing results of the raw coordinate time series.

Figure 9. Influence of tau selection (a. true × 0.1, b. true × 10, c. true × 1.0) on the post-seismic deformation model with different combinations (a. Offset + Velocity, b. Offset + Velocity + PSD, c. Offset + Velocity + PSD + Period).

Figure 10. Analysis results of time series fitting models with different combinations (a. Offset + Velocity, b. Offset + Velocity + PSD, c. Offset + Velocity + PSD + Period).

Figure 11. Analysis results of the PSD model with different forms (None, Exp, Log, Exp + Log).

Figure 12. Fitting results of ANTC’s time series using ITSM.

Figure 13. Fitting results of some sites with different trajectory characteristics.

Figure 14. The sites’ impacted seismic deformation with different PSD combination patterns.

Figure 15. The fitting RMSE of all GNSS time series using ITSM.

Figure 16. RMSE of global GNSS sites’ velocity fields in the north, east, and up directions.

Figure 17. Comparison of the ITRF2014 and GGV2020 global GNSS velocity fields. (a) ITRF2014 horizontal velocity field; (b) ITRF2014 vertical velocity field; (c) GGV2020 horizontal velocity field; (d) GGV2020 vertical velocity field; (e) ITRF14 and GGV2020 horizontal merging velocity field; (f) ITRF14 and GGV2020 vertical merging velocity field.

Figure 18. Difference of velocity field between ITRF2014 and GGV2020.

Figure 19. Fitting Residuals of ANTC based on the ITSM.

Figure 20. Residuals of ANTC using ITSM adding velocity discontinuity.

Table 1

Global GNSS observations’ duration.

Time Span (a)	Number	Percentage (%)	Mean (a)
0–2	215	19.03	10.48
2–10	467	41.33
10–30	448	39.65

Table 2

Default criteria of the ITSM.

Criteria	North (mm)	East (mm)	Up (mm)
weak data	20	20	40
outlier	20	20	40
very bad data	1000	1000	3000

Table 3

The fitting RMSE (Root Mean Square Error) of the GNSS sites impacted by seismic deformation.

RMSE	North (mm)	East (mm)	Up (mm)
min	1.19	1.62	4.21
max	7.27	7.03	13.64
mean	2.83	2.90	6.72

Table 4

Parameter fitting result of ANTC’s time series using ITSM.

	Velocity	Earthquake	CSD	PSD	Annual	Semi-Annual
	Slope (mm/a)Sigma (mm/a)	Time (a)Tau (a)	Offset (mm)Sigma (mm)	Exp and Log Coefficient (mm)Sigma (mm)	Amplitude (mm)Sigma (mm)Phase (rad)	Amplitude (mm)Sigma (mm)Phase (rad)
N	9.98	2010.1575	196.02	−6.94, 32.72	0.41	0.34
	0.70	0.2601	0.64	1.12, 0.39	0.14	0.14
					−0.0489	0.3514
E	13.07	2010.1575	−880.54	65.75, −165.70	1.13	0.50
	0.09	0.2601	1.01	1.65, 0.50	0.10	0.10
					0.2633	0.1044
U	3.53	2010.1575	−28.1	−17.67, 49.19	5.81	0.98
	0.13	0.2601	1.22	2.13, 0.75	0.15	0.15
					0.0042	0.6022