## 4 Discussion  
The spatiotemporal relative motion of the Sichuan regional reference station is influenced by the external environment, geological factors and other non-rectonic signals. The horizontal motion of the reference station is mainly concentrated in the southeast direction and the motion trend in the vertical direction is mainly uplift. The baseline components show mainly north-south compression and east-west stretching, which are consistent with the direction of crustal deformation in the Sichuan region, suggesting that the collision and wedging of the Indian plate with the Eurasian plate is still the main source of tectonic motion in mainland China [28]. Differences in processing software and strategies, noise models or details of the sites defining the reference system can affect the uncertainty of the velocity field [10]. The most advanced version was not used in the solution strategy of this study, such as the old one (IERS03) used for the solid tide correction model. ITRF2020 is the latest enhanced terrestrial reference frame that applies motion constraints to the complete time series of the four techniques with more reliable origin positions and accuracy over time [29]. We used ITRF2014 from 2015 to 2021 instead of the latest ITRF2020, which may also have a slight influence on the accuracy of the velocity field. Also, considering the complexity of the topography in the Sichuan region, the GNSS time series of the stations in the region inevitably contain a residual signal, which should be modeled deterministically or stochastically for improving the accuracy of the GNSS time series [30]. In addition, the lower time span (seven years) makes it possible for random noise to affect the amount of vertical coordinate change, making it unreliable to describe and compare velocities at each site in detail. Therefore, this study reflects the uplift in Sichuan region through the overall change. Seasonal variations, load model errors and component shares of thermoelastic variability in horizontal and vertical deformations have been explained [31]. If a cleaner time series is desired, noise models need to be considered in addition to filtered CME [15]. In the future, the rationality of using the PCA method can also be investigated by comparing the percentage of components with the corresponding seasonal variations and load model errors [11; 14]. In terms of modeling variance-covariance functions with random influence, the Gauss, Hirvonen and exponential functions mentioned in Section 2.4 are positively definite, and \begin{table}
\begin{tabular}{c c c c c c} \hline \hline  & **Statistics** & **Gauss** & **Hirvonen** & **Polynomial** & **Exponential** \\ \hline \multirow{4}{*}{N} & MAX & 3.58 & 3.14 & 2.03 & 2.17 \\  & MIN & 0.28 & 0.27 & 0.57 & 0.10 \\  & AVE & 1.78 & 1.52 & 1.22 & 1.20 \\  & RMS & 2.08 & 1.76 & 1.34 & 1.41 \\  & MAX & 4.21 & 3.14 & 1.47 & 3.09 \\  & MIN & 0.24 & 0.65 & 0.16 & 0.37 \\  & AVE & 2.12 & 1.72 & 0.66 & 1.33 \\  & RMS & 2.42 & 1.95 & 0.83 & 1.66 \\  & MAX & 0.44 & 0.51 & 0.25 & 0.38 \\  & MIN & 0.01 & 0.08 & 0.04 & 0.04 \\  & AVE & 0.16 & 0.26 & 0.13 & 0.17 \\  & RMS & 0.25 & 0.30 & 0.16 & 0.21 \\  & MAX & 3.98 & 10.85 & 2.12 & 4.03 \\  & MIN & 0.14 & 0.88 & 0.68 & 0.02 \\  & AVE & 1.53 & 5.02 & 1.38 & 2.40 \\  & RMS & 2.07 & 6.10 & 1.49 & 2.76 \\ \hline \hline \end{tabular}
\end{table} 
Table 4: Error statistics of fitting points (unit mm). thus the parameter k in the function model does not affect the model's positive definiteness to some degree. However, polynomial functions do not have positive definiteness, so it is important to determine the acceptable range of variation of the parameter k. Therefore, we focus on the acceptable range of variation of the parameter k in polynomial functions. The function models fitted in the paper are generally predicated on small regional scales, so baseline lengths of up to 300 km are generally chosen.