Analyzing the Seasonal Deformation of the Sichuan–Yunnan Region Using GNSS, GRACE, and Precipitation Data

: The Global Navigation Satellite System (GNSS) time series non-constructive deformation shows signiﬁcant seasonal variations, and the study of its periodic term components and possible physical mechanisms has important theoretical signiﬁcance and application value for the accurate use of the data and more in-depth analysis. In this paper, wavelet transform (WT) is used to extract the seasonal terms of 24 GNSS continuous station time series, Gravity Recovery and Climate Experiment (GRACE) displacement time series and precipitation data in the Sichuan–Yunnan area, and the three data sets are compared and analyzed in terms of amplitude, phase and cross-correlation coefﬁcient (CC), and the results show that the seasonal deformation in the area is strongly related to the precipitation variation. GNSS and GRACE have good consistency in the vertical component, and the seasonal variation is mainly related to the hydrological load; the difference in the horizontal component is obvious, and the amplitude of the seasonal term of GNSS is larger than that of GRACE, indicating that the resolution of GRACE in the horizontal component is lower than that of the vertical component, and the overall estimation accuracy is lower than that of GNSS. There are signiﬁcant seasonal terms of annual and semi-annual cycles for the three GNSS components, and the vertical component mainly shows the seasonal deformation of the annual cycle with the strongest seasonality; the horizontal component mainly shows the deformation of the semi-annual cycle, and the seasonality of the N component is stronger than that of the E component, and they are negatively correlated with a coefﬁcient of − 0.90.


Introduction
Changes in surface mass (e.g., atmospheric pressure, non-tidal oceans, and hydrological loading) can change the gravitational field and three-dimensional displacement of the elastic crust [1][2][3]. Many studies have found that seasonal variations in hydrological loading, which is majorly affected by abundant precipitation, are the main cause of seasonal crustal deformation [4][5][6].
GNSS, which accurately responds to tectonic and non-tectonic movements, is a highly precise multi-dimensional method that can be applied on large scales to observe crustal deformation. It is not only subject to large-scale regional tectonic movements, but is also particularly sensitive to changes in local observation environments; thus, it contains more complex motion information. Related research results show that the Chinese GNSS time series are mainly dominated by linear motion in the horizontal direction with periodic variations [7,8], while the vertical direction shows complex periodic oscillations and nonlinear motion. Moreover, the periodic fluctuations differ among different stations. Gravity Recovery and Climate Experiment (GRACE) satellite (include GRACE and GRACE-FO) observations comprehensively respond to terrestrial gravity field anomalies. Global timevarying gravity data provided by GRACE have been widely used to study terrestrial water changes [9][10][11], large watersheds and rivers [12], and changes in the mass of polar ice sheets [13,14]. According to the elastic load theory [1], the time-varying gravity data of GRACE can not only represent gravity field changes caused by changes in surface mass but can also estimate the corresponding three-dimensional crustal load deformation. However, affected by the resolution, it provides information on large-scale crustal deformation characteristics. Precipitation, as the main influencing factor of hydrological loading change, provides important information for the study of surface quality change. The theoretical displacements obtained by GRACE inversion combined with precipitation information can be used to better analyze seasonal deformation characteristics of GNSS time series from the perspective of seasonal loading studies.
Sichuan and Yunnan provinces are located on the southeastern margin of the Tibetan Plateau in China, a mountainous plateau terrain, high in the north and low in the south, and most of the region is located in the middle altitude area, mostly hills and mountainous geography. As the Eurasian Plate and Indian Plate collide, these provinces experience some of the strongest tectonic and seismic activities in mainland China. In addition, the region is located in southwestern China and experiences a subtropical monsoon climate with a small annual variation in temperature, a high daily temperature difference, distinct wet and dry seasons, and high precipitation during the rainy season, with a significant annual period ( Figure 1). Existing research on the seasonality of this region has mainly focused on vertical features, and annual cycles [15][16][17][18]. Our research found that the time series of GNSS continuous stations not only have significant annual periodic changes in the vertical direction but also in the horizontal direction; additionally, evident semi-annual cycle signals exist. However, studies on the seasonality and horizontal direction of semi-annual cycles are few. To bridge this research gap, this study aimed to comprehensively analyze the GNSS three-dimensional displacement coordinate time series, GRACE inversion modelling results of crustal displacement, and main periodic components of precipitation observation data; moreover, the possible interactions between these factors are discussed.

GNSS Data and Pre-Processing
The data integrity rate of the 24 continuous observation stations of the Crustal Movement Network of China (CMONOC), GNSS, selected in this study was >95%. The observation period of most stations was January 2011-June 2021. GAMIT/GLOBK software was used for calculations [19]. First, the single-day relaxation solutions of the International GNSS Service (IGS) station and regional GNSS stations were obtained using GAMIT. Antenna phase center was corrected using the absolute antenna phase center model. The atmospheric refraction model used Vienna Mapping Functions 1. Global Earth tide and ocean tide corrections followed the IERS 2003 Conventions. The IGS station was used as the fixed-frame station, and GLOBK was used to adjust and estimate the time series and errors of point coordinates of regional stations under the International Terrestrial Reference Framework 2014. Then, the steps of the time series caused by antenna replacement and co-seismic displacement of strong earthquakes were repaired, and outliers were eliminated. Finally, the surface loading model [20][21][22] is used for atmospheric and non-tidal ocean load model corrections, but not for hydrographic load corrections.

GRACE Data and Processing
In this paper, we use the Green's function convolution method to calculate the influence of regional grid point loading on the stations. Assuming that the mass load on the Earth's surface expressed in terms of the equivalent water height h w (θ, λ) is known, then the surface fluid mass load can be expanded into the summation form of the spherical associated Legendre function [23]: where θ, λ denote the residual latitude and longitude east of the load point; ∆C h nm , ∆S h nm are the spherical harmonic coefficients of the global surface load equivalent water height expansion; P nm is the 4π normalized associated Legendre function; n and m denote the order and number, respectively; N max is the maximum order of truncation, which is related to the spatial resolution. Assume that ρ e is the average density of the earth, ρ w is the density of water, R is the average radius of the earth, (θ, λ) and (θ , λ ) denote the residual latitude and longitude east of the point to be sought and the load point, respectively, and P n is the Legendre function; h n and l n are the load Leff number. According to the elastic load theory, the equation of the displacement of the surface in 3 directions caused by the mass load in the local Cartesian coordinates coordinate system: We used the load Love number provided by Guo et al. [24] and the RL05 equivalent water height data product with a resolution of 0.5 • × 0.5 • × 10 days provided by the Centre National d'Etudes Spatiales/Groupe de Recherche de Géodésie Spatiale (CNES/GRGS) [25] to calculate the crustal displacement of gravity load that could facilitate comparisons with GNSS data. Atmospheric and non-tidal ocean loading models were corrected while processing these data. To extract and analyze data periodic terms using the WT method and ensure the continuity of the signal frequency domain, we used the cubic-polynomial interpolation method to provide the missing data. Seasonal analysis was not performed for long periods of missing data between GRACE and GRACE-FO during 2017-2018.

Precipitation Data
To analyze the relationship between precipitation and crustal deformation observation data, we obtained daily mean surface precipitation for county-level administrative areas in Sichuan and Yunnan provinces for 2010-2021 from the China Meteorological Science Data Center. The selected 24 GNSS stations were located in these administrative areas.

WT Method
GNSS, GRACE, and precipitation time series data include different periodic components, such as observational errors, noise, seasonal loading deformations, and trend signals. WT can decompose a time domain signal into several stacked wavelet functions, which are obtained by translation and scaling of basis functions. Compared with the traditional Fourier transform method, WT can conduct multi-resolution analysis, and simultaneously characterize the local signal characteristics in the time and frequency domains. WT, a time-frequency analysis method, allows changes in both time and frequency windows. We used Morlet basis functions [26] and Fourier transform to analyze the periodic components of the time series; subsequently, we used Daubechies (db) basis functions [27,28] to extract the periodic terms, which were used for seasonal intensity analysis of different data at each station. The periodic term extracted by the seventh order db8 basis function corresponded to a daily sampling rate of 0.35-0.70 years, including semi-annual periodic components, whereas the seasonal term extracted from the eighth order corresponded to 0.70-1.40 years, including annual periodic components. The sum of the periodic components represented both annual and semi-annual cycles.

GNSS Data
Considering the SCPZ station as an example, we used the Morlet wavelet transform and Fourier transform to examine the main periodic components of the three components of GNSS time series NS, EW and UD in the time-frequency domain and the frequency domain, respectively ( Figure 2). We found 0.5 years and 1.0 year periods for the three components of the SCPZ time series. Other stations in the area showed that the semi-annual cycle of the horizontal component was stronger than that of the annual cycle, and the N component was stronger than the E component, indicating that the semi-annual cycle seasonal term should not be overlooked. In addition, other significant periodic components, such as 1.5 years, 2.0 years, and 3.5 years, were observed. Figure 3 shows the main frequency components of the 24 stations in three directions extracted using Fourier transform. As the effective observation time of most stations is approximately 10 years, based on the Nyquist sampling theorem, the longest period that can be accurately determined is five years; therefore, we only focused on periodic components of less than five years. The N and E periodic components were scattered and complex, with the main periods being 0.5 years and 1.0 year; additionally, a 1.1 years periodic term was observed. However, every station did not have an annual (e.g.,

GRACE Data
Based on the equivalent water height data, we obtained the GRACE displacement time series for different stations using Green's function. Obvious annual and semi-annual seasonal terms were observed in all three components (N, E, and U); however, the semiannual cycle of the E direction showed less seasonality, and the vertical component showed significantly more seasonality than the horizontal component. In addition, some stations had seasonal terms significantly longer than 1.0 year, such as 2.0 years, 3.0-4.0 years, and 6.0-8.0 years; however, these terms differed considerably between stations. Few studies have been analyzed this aspect, and thus, further in-depth analysis is required. Figure 4 shows the time series of the SCPZ station and its WT and Fourier spectra. To ensure consistency, this study selected GRACE data from 2010 to 2021 for analysis, which was similar to the GNSS observation period.

Precipitation Data
The region is rich in precipitation, characterized by less in the northwest and more in the surrounding areas, with a strong seasonality. The annual and semi-annual seasonal terms were significant, and repeatability was strong. The end of July experienced the highest precipitation each year. Figure 5 shows the daily mean precipitation at the SCPZ station and its WT. As precipitation differed considerably throughout the year, and the data values were consistently positive, the amplitude of the extracted periodic term was much smaller than the actual value (the red curve in Figure 5a). In this study, we mainly focused on the temporal changes in the phase. Figure 6 shows the annual means for the 24 stations after stacking the daily mean precipitation data, and the results show that the rainy season generally commenced in May in the Sichuan-Yunnan region and the daily mean precipitation peaked on 23 July (day of year, DOY: 204), with a standard deviation of ±12 days. Precipitation in July accounted for 20.5% of the total annual precipitation, and precipitation from May to October accounted for 85.6%. The maximum (YNRL) and minimum precipitation (YNYM) values were 1250 mm and 794 mm, respectively, with a mean value of 987.5 mm.

Seasonality Analysis
After extracting the seasonal terms of each station by WT, we analyzed their seasonal intensities from the perspectives of amplitude, phase, and CC. Amplitude and phase were obtained by stacking and averaging the annual periodic terms. Amplitude indicated the periodic displacement magnitude of a station and it was expressed as the root mean square of the displacement time series of each component. Phase was the DOY corresponding to the highest displacement value. The CC indicated the seasonal repeatability of each station, and it was expressed as a mean value of the annual periodic terms of each station.
Stacking can obtain more stable and accurate amplitude and phase of the periodic term of the station, which mainly takes advantage of the correlation of the signal and the noncorrelation of the noise. Some abnormal noise information (such as non-periodic changes in the load around the station, volcanic eruption, flooding, etc.) that cannot be corrected or verified by mathematical models can also be weakened by stacking. Theoretically, the amplitude of the interference information will be weakened by overlapping times, and the phase shift will be smaller with more stacking times, and the signal-to-noise ratio will improve at the rate of the square root of the stacking times.
Suppose the length of a station time series is m years, and the signals of any two years are x(i) and y(i) respectively, then the correlation coefficients are as follows: where x, y are the average values of x(i) and y(i), respectively, n is the signal length. The problem is a combination of any two years in m. All possible combinations are M = C 2 m . For all the combinations, the correlation coefficients is r j , j = 1, 2, . . . . . . , M. The cross-correlation coefficient is defined as: Figure 7 shows the GNSS time series without the trend term and the GNSS and GRACE periodic term curves were extracted by WT using the SCPZ station as an example. Table 1 shows the seasonal parameters of each station in the study area. Figure 8 shows the mean seasonal curve of the three components seasonal time series of GNSS and GRACE and the mean curve of the 24 stations both stacked by year. Both results showed seasonality. The amplitude and CC of the horizontal component was smaller than that of the vertical component. The seasonality of the N component was stronger than that of the E component, and GRACE showed stronger seasonality in the three components. These trends may be related to the linear inversion based on equivalent water height.   On comparing the N component of GNSS and GRACE, the GNSS results revealed obvious semi-annual periodic changes, consistent seasonality among different stations, and 93 and 269 as the DOYs of the two peaks of the mean seasonal curves, respectively, with an interval of approximately 0.50 years. The GNSS phases, appearing at the end of March or April, occurred before the GRACE phases, which were concentrated in May.
Regarding the E component, the GNSS results of different stations also showed significant semi-annual periodic seasonal terms, but due to the large differences between the annual periods of different stations, the stacked amplitudes negated each other to some extent, and the amplitude of the mean seasonal term curve reduced. Furthermore, DOYs of GNSS and GRACE also changed to 8 and 202, respectively, with an interval of approximately 0.53 years. The GNSS phase was unstable and occurred later than the GRACE phase, which occurred from February to April. The phase and amplitude of XIAG differed from those of other stations, probably because it was located southwest of and close to (<300 m) the Erhai Lake. This implied that the E component displacement was more likely to be affected by the water storage load of the lake, which had minimum and maximum water levels in May and October, respectively, each year [29]. Regarding the U component, the annual and semi-annual seasonal terms of the GNSS results differed significantly across the different stations, and the morphological consistency was strong. On stacking the GNSS and GRACE results, seasonal changes were observed in the annual cycle. The GNSS phase which generally appeared at the end of May and the beginning of June (except for SCXC, where it appeared in January), was later than that of GRACE, which appeared at the end of April and the beginning of May. Both phases were relatively stable.
The mean CC of the N, E, and U components for GNSS were 0.63, 0.41, and 0.87, respectively, whereas the values were all >0.93 for GRACE, indicating that GNSS and GRACE showed strong seasonality regarding the N and U components. Regarding the E component, the mechanisms affecting GNSS observations may be more complex; therefore, the annual periodic term repeatability was poor, and phase differences between stations were considerable.
Furthermore, precipitation varied greatly between the different stations in the study area, but the precipitation time is concentrated in May to October. The maximum and minimum CC of the stations were 0.97 and 0.89, respectively, with a mean value of 0.94. Thus, strong seasonality was observed.

Discussion
This study used WT to extract the seasonal terms of GNSS and GRACE displacement time series and analyzed the phase, amplitude, and CC of 24 stations in the Sichuan-Yunnan region. Precipitation data was also considered for analysis. The three data types showed notable seasonal variability and a certain correlation, but with clear systematic differences.
First, the GRACE results intuitively responded to changes in the surface loading within the study area, especially the vertical component, which was negatively correlated with changes in the surface mass. Referring to the U component in Figure 8b, precipitation in April was found to gradually increase with the beginning of the rainy season, and crustal displacement gradually changed from a positive maximum value to a negative maximum value due to load application. As precipitation reduced towards the end of the rainy season, surface uplift occurred, and similar changes occurred in the horizontal component. The difference between the two phases may be because precipitation records an instantaneous change, while GRACE records long-term cumulative changes in surface load. Consequently, the maximum displacement of the GRACE-stacked mean seasonal term (DOY: 254) was 50 d (1.7 months) later than the peak precipitation.
Second, the GNSS amplitude was greater than that of GRACE in both horizontal and vertical components, indicating that stations were not only affected by regional load as monitored by GRACE, but also by local factors to some extent. Thus, GRACE could not effectively resolve this problem due to the local variations in the estimation accuracy. Moreover, the periods and phases of GRACE and GNSS differed largely in the horizontal component, but marginally in the vertical component, indicating that the horizontal resolution of GRACE was lower than its vertical resolution [10]. In the vertical component, the stacked mean seasonal correlation coefficients of GRACE and GNSS were as high as 0.92, but the GNSS phase occurred 29 days later than that the GRACE phase, and their troughs were 20 days apart. This may be because the instantaneous hydrological loading deformation caused by precipitation was recorded simultaneously by both GRACE and GNSS, but changes in pore pressure in the subsurface medium caused by water diffusion were detected at different times [30]. If the pore pressure increases or decreases, pore volume may increase or decrease, accompanied by an increase or a decrease in the horizontal displacement of the surface, which can be recorded only by GNSS, consequently, leading to changes in the amplitude and phase. Pore pressure changes cause changes in the properties of the subsurface medium and cause differences in the seismic wave velocity changes from the GRACE phase. Liu et al. [30] indicated that the seasonal variation of seismic wave velocity occurs approximately 20 days later than the high phase of equivalent water table, which also supports this view. The horizontal component was strongly correlated between GNSS and GRACE, indicating that the actual GNSS seasonal deformation is more complicated in the horizontal component and cannot be explained well by a single load accumulation model.
Third, on comparing the results of different GNSS stations, high or low differences in the coordinate time series of stations were observed, indicating that GNSS is not only affected by large-scale non-structural factors but also considerably by local observation environments. For example, for non-single load areas around the observation stations, multiple loads usually increase the amplitude of the vertical component, but they may offset the amplitude and deviate the phase in the horizontal direction [10]; therefore, consistency between stations was better in the vertical direction, whereas larger differences were observed in the horizontal direction.
Fourth, in the GNSS horizontal direction, the N and E components exhibited a strong negative correlation (r = −0.90), and seasonal deformation was possibly controlled by changes in both hydrological loading and pore pressure. Precipitation was lower in the N component from January to March; therefore, the influence of hydrological loading gradually disappeared, and stations moved northward simultaneously as surface uplift occurred. In April, when precipitation increased, hydrological loading increased, and stations moved southward until precipitation peaked in July. Deformation also reached its extreme negative point. Subsequently, precipitation decreased, and instantaneous precipitation loading deformation weakened. Stations resumed their northward movement until October. Seasonal deformation from April to October might largely be due to an increase or decrease in an instantaneous load from precipitation. Precipitation accumulation in October peaked within a short period, and under the combined effect of other hydrological loadings, stations moved southward until December. Later, precipitation decreased sharply, and surface water was continuously lost due to surface runoff, evaporation, and other factors. The surface load gradually disappeared, and stations started moving northward again until March the following year. From October to March, the movement of the sites may be controlled by the dual effects of hydrological loading and pore pressure. The action mechanism of deformation features and influencing factors for the E component was opposite to that for the N component.
Fifth, during GNSS and GRACE data processing, in addition to semi-annual and annual terms, we found longer terms, such as >1 year, which may also be related to precipitation. Many studies [6,[31][32][33][34] have stated that precipitation in different areas of Sichuan and Yunnan has a periodic component from 1.0 year to 8.0 years or longer. This study used the mean precipitation of administrative areas where stations were located; however, this was not reflected well in the seasonal analysis of precipitation.

Conclusions
In this study, GNSS, GRACE, and precipitation data were combined to analyze the seasonal deformation of the Sichuan-Yunnan region. The following main conclusions were drawn: 1.
The three data types for Sichuan and Yunnan have significant annual seasonality. The GNSS time series not only has a notable annual cycle but also an evident semiannual cycle, and thus, the impact of the semi-annual cycle on the horizontal component should not be overlooked. Moreover, there are other shorter and longer periodic signals.

2.
A comparative analysis of GNSS and GRACE revealed that the vertical components are largely consistent, but the horizontal components differed considerably. GRACE provides comprehensive observation results for large areas, and has better vertical resolution than horizontal resolution. Regarding the actual observations, GNSS was affected by non-structural factors of large areas, and was easily influenced by the local observation environment.

3.
The seasonality of the vertical component of the GNSS three components was strong. The horizontal N component was stronger than the E component, and both were strongly negatively correlated. Seasonal deformation was mainly affected by instantaneous loading deformation and changes in the pore pressure due to precipitation.

4.
Regarding seasonal deformation, this study mainly discussed the influence of precipitation. Nevertheless, the GNSS observation stations were also affected by other factors, such as atmospheric pressure load, non-tidal ocean load, thermal expansion effect of bedrock at observation points, and data processing errors, which are also important factors that cause differences between GNSS and GRACE in terms of amplitude and phase. 5.
The correction of seasonal deformation of GNSS time series has important theoretical significance and application value for obtaining accurate information of surface movement and accurate use of observation data.

6.
Through the integrated application of GNSS and GRACE, it can also effectively monitor the regional land water storage changes or precipitation changes, and make up for the lack of weakness of professional technical means in some research areas.