General Geometric Model of GNSS Position Time Series for Crustal Deformation Studies – A Case Study of CORS Stations in Vietnam

In processing of position time series of crustal deformation monitoring stations by continuous GNSS station, it is very important to determine the motion model to accurately determine the displacement velocity and other movements in the time series. This paper proposes (1) the general geometric model for analyzing GNSS position time series, including common phenomena such as linear trend, seasonal term, jumps, and post-seismic deformation; and (2) the approach for directly estimating time decay of postseismic deformations from GNSS position time series, which normally is determined based on seismic models or the physical process seismicity, etc. This model and approach are tested by synthetic position time series, of which the calculation results show that the estimated parameters are equal to the given parameters. In addition they were also used to process the real data which is GNSS position time series of 4 CORS stations in Vietnam, then the estimated velocity of these stations: DANA (n, e, u = -9.5, 31.5, 1.5 mm/year), HCMC (n, e, u = -9.5, 26.2, 1.9 mm/year), NADI (n, e, u = -10.6, 31.5, -13.4 mm/year), and NAVI (n, e, u = -13.9, 32.8, -1.1 mm/year) is similar to previous studies.


Introduction
Monitoring the Earth's crust deformation, land deformation using continuous GNSS (Global Navigation Satellite System) stations is now very popular in the world, greatly contribute to the prediction and early warning of seismic activities (e.g., earthquakes, landslides) in the US [1,2], Japan [3,4], Taiwan [5], etc.
Continuous GNSS observation allows the calculation of the position time series of GNSS stations. The most basic application of GNSS position time series is in geotectonic studies, such as determining the rates of crustal motion, the seismic displacement, the co-seismic deformation, and the postseismic deformation [6]. Additionally, GNSS position time series are also useful in calculating amplitudes and phases of periodic motions caused by seasonal mass loads; for example, atmospheric, hydrological, or oceanic [7], etc.
Then, the study of GNNS position time series processing in order to determine the most accurate information in the time series such as displacement velocity, seasonal motions, postseismic deformation of the Earth's surface at the continuous GNSS station is very important for the applications mentioned above.
In the processing of the position time series of continuous GNSS station by geodetic approach, the synthetic motion model is often used as a geometric model according to the 3-D coordinate components (XYZ) or (n,e,u) [8]. These geometric models need the parameters describing the basic and the most common motions of the coordinate time series, which are the linear trend, seasonal terms, possibly jumps/offset, and post-seismic deformations. The simpler the model of the position time series, the faster the computation, but the accuracy of the velocity estimated is noticeably reduced, especially with complex position time series. Using a geometric model with only one linear trend in position time series processing is also common such as [9] in determining the crustal motion in Vietnam and in the surrounding area by continuous GPS, [10] in determining the horizontal displacement of the Earth's crust in the Northwest region of Vietnam by cycle measurement GPS, etc. However, this model does not fully describe the nonlinear motions and discontinuities in the time series used for these studies leading to the estimated velocity field not being accurate.
The geometric model, including linear trend, seasonal motion, and offsets, has been used by [11] to study the effects of the GPS position time series caused by higher-order ionospheric corrections. Similarly, [12] has also used this model to research the accuracy of the velocity of continuous GPS station, and [13] used it in detecting the jumps automatically in the GPS time series, etc. This model is only suitable for these position time series used in those studies; however, in the case of analyzing complex time series, a synthetic model is needed. This paper proposes an approach to process GNSS position time series using a synthetic geometric model including linear trend, the seasonal term (annual and semiannual), jumps, and post-seismic deformation, which is suitable for most GNSS position time series cases. This model has also been used in programming the GNSS position time series analyzing software [14,15], in the study of seismic deformation by GPS [16]. In this way, the time decay of postseismic deformation is estimated directly on the GNSS position time series, from which it is possible to determine the best fit line of time series. It means that the velocity, parameters of seasonal term, the amplitude of jumps, postseismic movements are estimated the most accurately. Moreover, the tool was written by Python 3 to analyze these data and further results. This tool has been tested with the synthetic model position time series and is used to analyze the real position time series of 4 CORS (Continuously Operating Reference Station) stations in Vietnam.

General geometric model of GNSS time series
GNSS position time series represents the temporal evolution of a point in the 3-D space in the geocentric coordinate system (X, Y, Z) or local geodetic coordinate system (north -n, east -e, up -u). In other words, the GNSS position time series represents the time variation of the coordinates of the GNSS station.
Continuous GNSS stations are built on the land surface, so it is also affected by land surface changes such as displacement, periodic movement due to changes in groundwater levels, geo-tidal, etc., by geotectonic activities such as earthquakes. These influences are evident in the position time series of the GNSS station. The simplest geometric model for a position time series in terms of coordinate components n, e, u is a linear trend [8]: Where i t is the time of the series positions in a decimal year, a is the initial position at time reference, b is the linear velocity, and v is the measurement errors. The more common model is a linear trend with the addition of seasonal term (includes annual and semiannual periodic motion) and jumps in the equation (1) [12,17,18] Where c, d, and e, f are harmonic components of annual motion and semi-annual motion, respectively.
A jump (step or offset) in the position time series is a sudden change in the mean coordinates [18], which can occur one or many times in a GNSS position time series at different times with various magnitudes. The cause of the jumps is probably GNSS receiver changes or receiver's firmware update [19], or geophysical phenomena such as land subsidence, earthquakes [17].
The seasonal term is sinusoidal motions with periods of 12 months (annual motion) and six months (semi-annual motion). This is due to cyclical changes over time, such as the effect of the continental water storage loading [20], or the effect of seasonal change of groundwater level [3], etc. The general geometric model proposed for the GNSS position time series is equation (2) with the addition of the post-seismic motion [8, 14, Where j h is rate change after an earthquake at the time hj T ( h n is the number of rate changes). j k is the magnitude of postseismic relaxation at time kj T ( k n is the number of postseismic).  is decay time (relaxation time) of postseismic deformation in the decimal year.
Thus, the parameters of the general geometric model of GNSS position time series are a, b, c, d, e, f, In equation (3), the postseismic deformations that occur after the jumps of earthquakes are described in the square brackets ([]), which includes the rate changes (change in both direction and value) and postseismic relaxation (described as an exponential decay function with magnitude and decay time).
The GNSS position time series record the seismic, co-seismic, inter-seismic, and postseismic relaxation [21]. Postseismic relaxation can last for years or more. In addition, after an earthquake, the rate of displacement of the GNSS station can change in both direction and value. The earthquake jumps and postseismic relaxations depend on the epicenter, earthquake's magnitude, and the distance between the epicenter to GNSS station [8].

Determining the parameters of GNSS position time series
Assuming that the times gj T of jumps, hj T of rate changes, and kj T of postseismic relaxation are known. Model (3) is linear with respect to the parameters: v is a vector of residual position time series, and y is a vector of position time series, 1 2 [ ( ) ( ) ...
( )] T n y y t y t y t = (8) The positions in the series are independent observations, then the observation covariance matrix is defined by individual variances i  of the position i: The least square for the best linear unbiased estimates of the parameters is: With parameter covariance: Where  is weighted root mean square error: From equation (5), the modelled postions (predicted position) are: And residual positions are calculated as: The amplitude, the phase of the annual and semi-annual motions are calculated from the parameters c, d, e, f [22]: The flowchart describes the process of calculating and analyzing the GNSS position time series in Fig. 1. In the flowchart (Fig. 1), among the input coefficients for calculation, the decay time (relaxation time) of postseismic deformation τ is known. This coefficient is proposed by estimating directly from the time series and presented in the following section 2.3.
The program has been developed for analyzing GNSS position time series in Python, named pygps_ts (GPS time series by Python), which uses the method described above. Pygps_ts helps to solve the most common problems of time series quickly and accurately.

Determining the time decay of postseismic deformation
The postseismic relaxation time (decay time) is usually determined based on seismic models [23], or the physical process seismicity, or by Omori's law [24,25]. In this paper, we propose a method to calculate the decay time directly from the position time series based on the equation of the general geometric model (3). τ is iteratively calculated in 2 steps: First step, choosing the initial value τ 0 (here, we choose τ 0 = 0.0027 yr, equivalent to 01 day [15]). Substituting τ0 and (time known of seismic jumps) into equation (3) The system of equations (18) in matrix form: A y = Where, Unknown vector: Vector of free coefficient: And design matrix A as equation (24).
Then, the unknown vector  is esimated using least squares estimation: Second step, calculating the decay time τ in the next iteration i: The iterative process stops when the maximum of ∂τj is less than the threshold value, and therefore, convergence is achieved. In practice, we choose the threshold value at 0.003 yr [15].
The process of calculating the decay time of postseismic deformation is summarized in the flowchart      In pygps_ts script, we have programmed the τ determination module using this algorithm which will be used if the input data has the times of seismic jumps.

Synthetic data model
It is examined that the performance of this algorithm by a synthetic model produced as a north component of the daily position time series of the GNSS station named SYNT (it mean synthetic) for the period 2015.0-2021.0.
The synthetic model coordinate time series of SYNT is introduced: seasonal motion, coseismic jump at 2018.0, and postseismic relaxation, whose coefficients/parameters are in Tab. 1. Then, this model time series is perturbed by the Normal distribution law σ = 1.0 mm to create the time series of random coordinates which are illustrated as blue points in Fig. 3 on the left.
Using pygps_ts to process SYNT's coordinate time series, the estimated parameters are obtained and shown in Tab. 1. In addition, the predicted and residual coordinate time series are shown in Fig. 3. Note that the number of iterative estimations of decay time τ is 4, and the processing duration is small, only seconds. Tab. 1 shows that the parameters estimated by using pygps_ts are almost equal to the given parameters of the model time series. The difference between them is mainly in the parameters of seismic motion but is considered insignificant. The test result and also our processing results in practice have proved that the proposed geometric model, method, and pygps_ts program for processing GNSS position time series are appropriate and ensure accuracy.

Some CORS stations in Vietnam
GNSS CORS stations are not only used as a control network for surveying services [26] but also widely used in tectonic plate movement monitoring [27,28], in landslide monitoring [29] (landslide due to mining is regularly monitoring in Quang Ninh [30,31]), etc. This study processes the GNSS position time series of 4 CORS stations in Vietnam: NAVI (Ha Noi), NADI (Nam Dinh), DANA (Da Nang), and HCMC (Ho Chi Minh), which are actual data to verify the algorithm and pygps_ts program.
These CORS stations measure continuously from January 2016 to December 2018, of which 3 stations (DANA, HCMC, and NADI) are designed and operated by Tuong Anh Science and Technology Equipment Joint Stock Company, and one station (NAVI) belongs to The Japan Aerospace Exploration Agency (JAXA). The data of these CORS stations has been [32] processed by the PPP method using PPPC (Precise Point Positioning by using C language software) software to obtain the position time series.
The pygps_ts is used to process the position time series of these 4 CORS stations. The estimated parameters of all these time series only have the linear trend and the seasonal motion according to the 3-D coordinate components n, e, u, but no jumps, no postseismic deformation (illustrated in Fig. 4). The estimated parameters and their accuracy are shown in Tab. 2.
From the estimated parameters describing the seasonal motion, pygps_ts calculates the phase and amplitude of annual and semi-annual motions are shown in Tab. 3. The position time series, the fit line, and the residual position time series of 4 CORS stations are shown in Fig. 4.

Discussions and conclusions
The general geometric model mentioned in this study includes common motion phenomena such as linear trend, seasonal motion (including annual and semi-annual displacement), jumps, and postseismic deformations (including rate change and postseismic relaxation) is the complete model for GNSS position time series. However, it is rare that the transient deformation is recorded in the GNSS position time series in areas with complex geotectonic activities such as American Basin and Range [33], Guerrero Mexico [34], etc. Transient deformation is defined as a nonperiodic, nonsecular accumulation of strain in the crust [35], so it is difficult to model geometry like other motions. Note that this phenomenon has not been recorded in GNSS station Vietnam; therefore, it is not mentioned in this paper.
To calculate the magnitude of jumps in the position time series, it is necessary to provide the time of the jumps. The study of automatically determining the jump has also been researched [13,15] but needs to be improved and examined more closely.
The proposed method of determining the decay time directly from the coordinate time series allows a more accurate description of postseismic deformation. This method's algorithm is programmed in Python, allowing fast and accurate calculations with the data model. However, the GNSS position time series of CORS stations in Vietnam do not record postseismic deformation to verify the method. Therefore, it needs to be tested with real complex data to improve the method as well as the pygps_ts program.
The general geometric model of the GNSS position time series describes most of the common motions, which is true for most of the position time series in the study of land surface deformation, crustal deformation due to seismic, mining, etc. The calculating results of 4 CORS stations in Vietnam are highly accurate, and the processing time is fast because the data of these CORS stations are quite "beautiful". In the upcoming time, the study will be conducted with more CORS station data in Vietnam.
Vietnam is in the process of building and developing continuous GNSS networks for geotectonic study, surveying services, etc. the study of methods of processing GNSS position time series is essential to promote the role of GNSS networks in science and practical applications.

Acknowledgments
This research has been supported by the Ministry of Education and Training of Vietnam (under grant B2020-XDA-05). The CORS data is provided by Tuong Anh Science and Technology Equipment Joint Stock Company and Japan Aerospace Exploration Agency. We acknowledge Ngoc Lau NGUYEN for providing us the GNSS position time series of 4 CORS stations.
The paper was presented during the 6th VIET -POL International Conference on Scientific-Research Cooperation between Vietnam and Poland, 10-14.11.2021, HUMG, Hanoi, Vietnam.

Data Availability Statement
All data, code that support the findings of this study are available from the corresponding author upon reasonable request: -Code of pygps_ts in Python 3.
-Position time series of SYNT and 4 CORS station.