## 3 Results
### 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 non-correlation 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  
Figure 5: Time series analysis of precipitation at SCPZ station. (**a**) The blue histogram is the daily average of precipitation and the red line is the extracted seasonal term. (**b**) Morlet wavelet transform, where red means peak and blue means trough. Figure 6: Annual average results of precipitation superimposed on 24 stations. (**a**) Annual and semi-annual seasonal terms extracted by Morlet wavelet transform. (**b**) Monthly average value of precipitation and percentage of annual precipitation. 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:  
\[r=\frac{\sum\limits_{i=1}^{n}{(x(i)-\overline{x})(y(i)-\overline{y})}}{\sqrt{ \sum\limits_{i=1}^{n}{(x(i)-\overline{x})^{2}\sqrt{\sum\limits_{i=1}^{n}{(y(i) -\overline{y})^{2}}}}}} \tag{5}\]  
where \(\overline{x}\), \(\overline{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_{m}^{2}\). For all the combinations, the correlation coefficients is \(r_{j}\), \(j=1,2,\ldots\ldots,M\). The cross-correlation coefficient is defined as:  
\[R=\frac{\sum\limits_{j=1}^{n}{r_{j}}}{M} \tag{6}\]  
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. Figure 7: GNSS and GRACE time series of SCPZ stations. Figure 8: Time series averaging results of GNSS **(a)** and GRACE (**b**) seasonal terms for 24 stations. (The colored lines are the average seasonal terms of each station, and the black dashed lines are the average seasonal terms of the 24 stations stacked).