Analysis of Crustal Movement and Deformation in Mainland China Based on CMONOC Baseline Time Series

: In this paper, we propose a method for the analysis of tectonic movement and crustal deformation by using GNSS baseline length change rates or baseline linear strain rates. The method is applied to daily coordinate solutions of continuous GNSS stations of the Crustal Movement Observation Network of China (CMONOC). The results show that: (a) The baseline linear strain rates are uneven in space, which is prominent in the Tianshan, Sichuan-Yunnan, Qinghai-Tibet Plateau, and Yanjing areas, with a maximum value of 1 × 10 − 7 a − 1 , and about two orders smaller in the South China block, the Northeast block, and the inner area of the Tarim basin, where the average baseline linear strain rates are 1.471 × 10 − 9 a − 1 , 2.242 × 10 − 9 a − 1 , and 3.056 × 10 − 9 a − 1 , respectively; (b) Active crustal deformation and strong earthquakes in the Xinjiang area are mainly located in the north and south sides of the Tianshan block; the compression deformations both inside the Tarim block and in the southern Tianshan fault zone are all increasing from east to west, and the Tarim block is not a completely “rigid block”, with the shrinkage rate in the west part at about 1~2 mm/a; (c) The principal directions of crustal deformation in the Xinjiang, Tibet, and Sichuan-Yunnan regions are generally in the north—south compression and east—west extension, indicating that the collision and wedging between the Indian and Eurasian plates are still the main source of tectonic movements in mainland China.


Introduction
China is located in the southeast of the Eurasian plate. Due to the continuous collision and pushing between the Indian plate and the Eurasian plate in the southwest, a strong deformation area with a width of thousands of kilometers is formed in mainland China, and it also affects the hinterland of Central Asia, Southeast Asia, and Eastern China [1][2][3][4]. The unique tectonic background has made China one of the countries with the most severe seismic disasters in the world, and has also made the formation and evolution of tectonic deformation a hot spot in earth science research.
One of the effective means of studying the mechanism of plate movement is to measure the displacement of plate boundaries [5]. By the end of the 20th century, the main sources of data for understanding tectonic deformation were field fault surveys, satellite images, and seismic moment tensor inversions [6][7][8]. In recent years, with the rapid development of modern space geodesy, mainly represented by GPS, it has become possible to detect crustal deformation with high precision and high spatiotemporal resolution [9][10][11][12][13]. The United States has deployed a GPS network with hundreds of continuous stations in Southern California to monitor crustal deformation at the boundary between the Pacific and North American plates [14]. Japan has set up a GPS monitoring network with an average interval of about 30 km and more than 1200 stations covering the whole territory, mainly for monitoring crustal deformation [15]. China established a crustal movement monitoring network in 2000, which included 27 continuous GPS stations and more than 1000 campaign stations. On this basis, the Crustal Movement Observation Network of China (CMONOC) was built in 2012. At present, CMONOC has 267 continuous GPS stations and more than 2000 campaign stations to realize comprehensive monitoring of multiple sub-tectonic block movements, and also the atmosphere of near-earth surface [16]. Currently, CMONOC is mainly used for earthquake mechanism research, geodetic frame refinement, disaster warning, and so on [11,[17][18][19].
The complexity of continental tectonic deformation is well known. Whether it is a local geological map recording the evolution of tectonic movement or a large-scale seismic activity map, it is difficult to quantify crustal deformation accurately [20]. Based on a large amount of geological survey data, Ding et al. [6] estimated the average movement rate of the main tectonic blocks in mainland China earlier and gave the preliminary boundary of each tectonic block. Armoijo et al. [21] and Avouac et al. [22] used field survey data and Spot satellite images to deeply study the geological and tectonic movements in Central Asia and deduced the contraction and sliding rates of the Tianshan Mountain and Kunlun Mountain faults since the Holocene. Liu et al. [23] calculated the rate of horizontal crustal movement on the Tibetan Plateau, which is in better agreement with the results solved by geological models.
With the improvement of GPS observation precision and the increasing density of GPS stations, the large-scale strain field has become an important indicator for describing crustal deformation. Researchers have carried out a series of studies in crustal strain, stress, and geodynamics based on GPS observations, and have achieved many meaningful results [5,[24][25][26][27][28][29]. The strain rate field in mainland China was obtained by Shen et al. [13] and Jiang et al. [30] based on GPS velocity field data using regression least squares [31] and the least squares collocation method [32], respectively. Most of these studies are based on the assumption of uniform infinitesimal strain, using discrete velocities or modeled velocity fields to calculate the crustal strain rates. In the boundary of the block or the area where the GPS station is relatively sparse, the reliability of the obtained strain rate is low. On the other hand, any deformation in a solid medium is essentially related to the variation of distance; therefore, it is more straightforward to study crustal deformation by using the variation of baseline length. In addition, the analysis of the spatiotemporal variation of the baseline length is not limited by the uniform infinitesimal strain assumption, and the baseline can also be freely chosen according to our specific research purpose.
In this paper, by using the daily coordinate solutions from the continuous GPS stations of CMONOC, the time series of baseline length between any two chosen stations can be formulated, after which the baseline length change rate or the baseline linear strain rates can be calculated. Based on the spatial distribution of baseline linear strain rates, the strain level of major tectonic blocks in mainland China is analyzed, and also the relative movements of some active tectonic blocks in Xinjiang are given so as to show the applicability and flexibility of the baseline length change analysis method in exploring crustal deformation and tectonic movements.

GPS Data Set
There are 246 continuous GPS stations in CMONOC, 27 of which are the first phase stations. The starting epoch of these first phase stations was around 1998-240 (the 240th day in 1998), and the ending epoch we used was 2019-273, with about 21 years of observations. The other 219 GPS stations were built in the second phase, evenly covering mainland China, especially the areas where tectonic movements are active, most of which started in 2010-2011. The last observation epoch we used was the same as that of the first phase-i.e., 2019-273-with the overall observation time period at 9-10 years. Figure 1 shows the distribution of all these continuous GPS stations and also the division of the first-order tectonic blocks. Generally, the velocity of station caused by crustal movement is about several millimeters per year [29], which requires the accuracy of daily coordinate solutions of the GPS station to reach millimeter level. Therefore, it is necessary to ensure not only highquality observations from GPS receivers, but also the processing of the GPS data following stringent procedures [33]. At present, there are three kinds of high-precision GPS data processing software available: GAMIT/GLOBK [34], GIPSY [35], and Bernese [36]. For this paper, the Bernese software (version 5.2) was used to obtain high-precision daily coordinate solutions of GPS stations. Double-phase difference observations of chosen baselines were solved first for precise baseline solutions. Then, network adjustment with CDR (Control Datum Remove) and MC (Minimum Constrain) were used to obtain daily coordinate solutions from all GPS stations. Details of our GPS data processing is presented in [37]. After that, the daily baseline length between any two stations in the GPS network could be calculated.

Ellipsoidal Baseline Length Calculation
Considering the relatively low accuracy in the vertical components of daily coordinate solutions, the ellipsoidal distance between two stations are used in the following baseline length change analysis. Assuming that the three-dimensional coordinates of stations A and B are A(X 1 , Y 1 , Z 1 ), B(X 2 , Y 2 , Z 2 ), respectively, then the ellipsoidal baseline length S between the two stations would be as follows [38]: In order to facilitate the error estimation of the baseline length by using the law of error propagation, we approximate as follows: Partial derivative on both sides of the Equation (2): Based on the law of error propagation, Among them, σ S is the standard deviation of the baseline length S and σ X , σ Y , and σ Z are the standard deviations of coordinate components X, Y, and Z, respectively, which are obtained from the daily GPS network adjustment.

Baseline Linear Strain Rate
Assuming that S 1 is the baseline length before the change, and S 2 is the baseline length after the change, then the linear strain ε is defined as: The linear strain rate . ε is: .
where ∆t is the time interval of the baseline length change.
In practice, the linear strain rate is obtained by estimating the trend of the baseline length time series. Frequently used methods for estimating the trend of time series include least squares, Theil-Sen, and MIDAS (Median Interannual Difference Adjusted for Skewness). However, the least squares method is not robust; when outliers appear in the time series, the estimated trends are often biased [39,40]. The MIDAS method is a non-parametric and robust trend estimation method based on the improvement of the Theil-Sen method, which is insensitive to sudden jumps, outliers, and seasonal signals with a period of less than one year [40]. Here we give a brief introduction to Theil-Sen and MIDAS estimation methods; see Blewitt [40] for details. Suppose the linear model of the baseline length time series is: where S is the daily length of the baseline,m is the change rate of the baseline length to be estimated (trend term), t is the time, andd is the offset to be estimated. The Theil-Sen method [41,42] can be expressed as: where t j , t i are the time corresponding to the j and i epochs, respectively, while S j and S i are the baseline lengths corresponding to the j and i epochs. Evidently, this method selects the changes of all possible baseline pairs in the time series to calculate the corresponding rate of change and takes the median as the final trend item estimate. It can be seen from Equation (8) that the time step of the Theil-Sen method takes all possible values, making it sensitive to periodic signals. In response to this problem, the improved Theil-Sen method sets the time step interval at one year, which is: But in practical application, it is found that when the time series contains sudden jumps, the trend item estimated by the improved Theil-Sen method is biased. Therefore, based on the improved Theil-Sen method, MIDAS uses the following steps to attain a more robust estimation:

1.
Combining Equations (8) and (10), the trend term estimation and its standard deviation can be obtained as follows: Remove the value greater than 2 times of the standard deviation from m(k) to obtain m(p), that is: 3.
Then take the median of m(p) to get the final trend item estimatem, that is: In the above equation, σ is the standard deviation, which is given under the assumption that most of the difference, m(k) −m, obeys the Gaussian distribution. The scale coefficient of 1.4826 is an empirical value. Finally, use m(p) and Equation (13) to obtain the standard deviation of the trend term estimation.   Figure 2 shows that although the baseline length varies slightly from year to year, even pronounced periodic signals and small sudden jumps exist, and the trend of the baseline length change can still be found throughout the long period of observations. Comparing the trend terms obtained by different estimation methods, we can see that when the periodic signal in the baseline length series is significantly regular and the sudden jumps are less, then the trend terms obtained by the three methods do not differ much, as shown in Figure 2a,c. When there are short-term surges in the time series (as in Figure 2b), the result of Theil-Sen method and MIDAS method are almost unaffected and both can estimate the trend term reliably; on the contrary, the least squares method is significantly affected and the estimated trend term is biased. When longer-term sudden jumps occur in the baseline time series, such as the sudden jump in 2015 in YNYM-YNJD (Figure 2d), it leads to some bias in the estimates of both least squares and Theil-Sen, with least squares especially suffering the most, while the trend term estimated by the MIDAS method looks even and realistic. Therefore, in the following analysis, only the MIDAS method is used for the estimation of the baseline length change rate.

Baseline Selection and Pre-Analysis of Daily Coordinate Time Series
In order to estimate the baseline length change rate of CMONOC, the Delaunay triangle is used to connect GPS stations; a total of 681 baselines are obtained (see Figure 3). Since the observation start and end time for each baseline may be different, and also some stations may lack one or more days of observations as even gross error may occasionally occur in daily coordinate solutions, these factors will adversely affect the trend item estimation of the baseline length change. Therefore, we perform the following pre-analysis on the daily coordinate time series: 1.
Aiming at the problem of missing observation data at each station, Lagrange interpolation is used to fill in the missing data.

2.
In view of the difference of the start and end observation times between the stations, only coordinate solutions with a common observation period is chosen.

3.
For the outliers in the daily coordinate solutions, if they are directly filtered (such as sliding filtering), the correlation between the daily solution coordinates will be enhanced, which will be harmful to the subsequent analyses. Therefore, we compare the original coordinate time series with that of the sliding window filtering (window size is seven days); If the absolute value of the difference between original coordinate and the filtered one is greater than two times of the standard deviation, the original coordinate is considered to be an outlier, which is then replaced by the corresponding filter value.

CMONOC Baseline Linear Strain Rate Distribution
According to Equation (6), ∆S = S · ∆t · . ε, ∆S is the cumulative baseline length change in the ∆t period when the linear strain rate of the baseline with length S is . ε. From Equation (4), it can be seen that due to the precision of the baseline measurement, if the cumulative length change of a baseline during the observation time is less than the baseline measurement error σ S , it is considered not suitable for the calculation of the linear strain rate and should be eliminated from the following data processing. Considering that the daily coordinate solution of a GPS network has a precision of 3~5 mm in the horizontal direction and 6~8 mm in the vertical direction [43]. Here, we assume that the precision of each horizontal component of the daily coordinate solution is 3 mm, 4 mm, 5 mm, and the corresponding precision of the vertical component is 6 mm, 8 mm, and 10 mm, respectively. The 681 baselines in Section 3.1 are grouped into six categories according to their length, and the average baseline standard deviation in each group is calculated using Equation (4). The results are shown in Table 1 (unit: mm). As show in Table 1, under the same standard deviation of coordinate component, there is no significant difference in standard deviation of baseline within different length groups. However, in the same baseline length group, the standard deviation of baseline positively correlated with the standard deviation of coordinate components. That is, the precision of baseline length is positively correlated with the precision of the coordinate component solved by the daily GPS network adjustment, but has no significant relation with the length of baseline.
Combined with the average precision of the coordinate solution, the average standard deviation of the baseline length is estimated to be about 7.0 mm in our data set; we therefore set 7.0 mm to be the threshold of cumulative length change of baseline. That means that only the baseline length change that exceeds the threshold will be used for further analysis. From Equation (4), the time required for the cumulative baseline length change to exceed the threshold value depends on both the baseline length and the baseline linear strain rate. Figure 4 shows the time required for the cumulative baseline length change ∆S to reach the threshold of 7.0 mm, and the baseline linear strain rate is supposed to be in the range [1, 10] (10 −8 /a). The longer the baseline and the faster the linear strain rate, the less time is required for the accumulation of observable baseline length change.
Existing studies have shown that the principal strain rate in mainland China is mainly at the order of 10 −8 /a [1,13,30]. In this case, we take the smaller value of 1 × 10 −8 /a (the filtering conditions are more stringent), and calculate the minimum observation time period required for the 681 baselines in Section 3.1. If the time period of the baseline time series is less than the corresponding minimum time, the baseline will be eliminated from further analysis, and a total of 653 baselines will be left. Based on the MIDAS method, linear strain rates are estimated and drawn in Figure 5.  From Figure 5, we find that the baseline linear strain rates in South China, Northeast China, and the Tarim block are relatively small, and that the average baseline linear strain rates are 1.471 × 10 −9 /a, 2.242 × 10 −9 /a, 3.056 × 10 −9 /a, respectively. This is consistent with the research results based on the GPS velocity field [1,13,30], indicating that the tectonic movement of these blocks is relatively stable, and that the crustal deformation is subtler. The large baseline linear strain rates are mainly located in the Qinghai-Tibet Plateau, Tianshan Mountains, Sichuan-Yunnan, and Yanjing regions, where the recent tectonic movement is active [7,25,28]. The baseline linear strain rate in the Qinghai-Tibet Plateau region reaches 2 × 10 −8 /a, and manifested as north-south shortening and eastwest extension, which indicates that part of the north-south shortening caused by the Indian plate pushing northward against the Eurasia plate is absorbed by the east-west extension within the Plateau. Tianshan Mountains and Junggar Basin in Xinjiang exhibit compression movements on the whole, as the baseline linear strain rates increase from east to west, and reach up to 4.761 × 10 −8 /a in the far west. Li et al. [26] obtained the principal strain rate of 4.387 × 10 −8 /a based on block partitioning and strain analysis in this area, which agrees well with the results of this study.
The Sichuan-Yunnan region is located in the southeastern corner of the Tibetan Plateau, i.e., the Sichuan-Yunnan rhombus block, which is one of the regions with the strongest seismic activity in mainland China. The results of this paper show that the baseline linear strain rate in the Sichuan-Yunnan region is large, up to 10 −7 /a, and generally shows compression in the northwest-southeast direction and extension in the southwestnortheast direction. Active tectonic studies show that the area has a trend of movement in the southeast direction [7,13,28], which is basically consistent with the results of this paper. However, due to the more intricate distribution of faults in the area, the crustal deformation is somewhat heterogeneous, which is one of the main reasons for the anomalous changes in some of the baselines in the area. Figure 6a shows that the magnitude of the baseline linear strain rate is mainly distributed in the interval [−4, 4] (10 −8 /a), with a maximum value of 10 −7 /a, which looks similar to a normal distribution, indicating that the linear strain rate of most baselines is small and induced by the movement of rigid blocks [20]. It can be seen from Figure 6b that the standard deviations of the baseline linear strain rates are mainly concentrated in the order of 10 −11~1 0 −10 , and a few larger ones reach the order of 10 −9 . Hence, the reliability of the baseline linear strain rate estimation is high.

Analysis of Baseline Length Change Rates and Crustal Deformation in Tianshan Area
The Tianshan area is located in the Cenozoic regenerative orogenic belt within the continent, under the compression of the tectonic blocks. The far-field effect of the Indian plate pushing and the strong compression in the Pamir Plateau are the basic factors of the deformation in this region [3,[44][45][46]. As the boundary zone of tectonic blocks, the north and south flanks of Tianshan Mountains are transition zones for the Tarim block and the Junggar block where the main regional deformation occurs [2,8,47,48]. Since 2010, a total of nine earthquakes of magnitude six or above have occurred in Xinjiang, indicating that the present crustal movement in this region remains high. In order to further study the crustal deformation in Tianshan area, we calculated some chosen baseline length change rates; the lengths of these baselines are similar in general, and the results are shown in Figure 6. Figure 7 shows that the baseline length change rates and the distribution of earthquakes with magnitude greater than six from 1967 to 2020 (Earthquake catalog data source: www.globalcmt.org (accessed on 23 Novenber 2020)). In order to study the tectonic movement of this area, baselines were intentionally chosen across or along the block boundaries or fault zones. We found that the chosen baselines were mainly shortened, and only a slight tensile deformation was found in the northeast of the Tianshan block (see upper right of Figure 7), which is consistent with the results obtained based on the strain field and GPS velocity [13,49]. In addition, the baseline length change rates inside the Tianshan block and Tarim block were relatively small, but the baseline length change rates of those baselines across the Tianshan block and across the Junggar block were relatively large. The baseline length change rates of those baselines across the South Tianshan fault increased from east to west. Earthquakes are mainly distributed at the boundary between the northern part of the West Kunlun block and the Tianshan block, as well as in the inner part of the Tianshan block, and also in the northern Tianshan fault zone. In order to show the shrinkage of the Junggar block in more detail, the baseline length time series of those baselines across the Junggar block are drawn in Figure 8.
In Figure 8, the time span of each baseline is from 2010 to 2019. Although there are a few sudden jumps in the time series, they all have a significant negative linear trend. Combined with Figure 7, it can be noted that the length change rates of the baselines crossing the Junggar block show different amplitudes of extrusion. In Figure 8, the baseline length change rates from east to west are −1.23 mm/a, −3.22 mm/a, −1.33 mm/a, and −3.29 mm/a; the average shrinkage rate of the Junggar block is about 3 mm/a. The Tianshan block adjacent to the Junggar block stretches 2500 km from east to west and is one of the largest inland orogenic belts in the world [8]. Figure 9 shows baseline length change time series of four chosen baselines across the Tianshan block (XJSS-XJML, XJKE-XJSH, XJKC-XJDS, WUSH-XJWQ).     The length change of the baseline XJSS-XJML, which is located in the Tarim block and the eastern part of the South Tianshan Fault, has obvious periodicity. Even though the length change rate is small, the weak contraction can still be detected through the long time series. The baseline WUSH-XJWQ has the largest length change rate, reaching −10.50 mm/a. Combined with Figure 7, it can be observed that the length change rates of the four chosen baselines gradually increase from east to west, indicating that the crustal convergence amplitude of the Tianshan block increases from east to west. The South Tianshan fault between the Tarim block and the Tianshan block has been an intensive earthquake zone. Three baselines-XJKE-XJML, XJKC-XJBY, and XJBC-XJZS-across the South Tianshan Fault also exhibit compressive shortening with magnitude gradually increasing from east to west. Figure 10 shows the time series of baseline length change of six baselines from east to west inside the Tarim block. The direction of baselines is approximately north-south (see Figure 7), and the time span is from 2010 to 2019. It can be seen that these six baselines all show compressive shorting. Except for the baseline XJBC-XJYC, the other baselines all have recognizable periodic length changes. According to Figure 7, the baselines XJJJ-GSAX and XJSS-GSDH are located in the easternmost part of the Tarim block, which are with very small length change rates, and almost no change at all in the past 10 years. The length change rates of the baselines XJBC-XJHT and XJBC-XJYC in the west of the Tarim block are relatively large, which shortened by about 15-20 mm in the past 10 years; this means that the west of the Tarim Block is still shrinking at a rate of 1~2 mm/a, and the compression deformation in the west is stronger than that in the east.

Discussion
The average compression rate of the Junggar block in the Quaternary is 2-6 mm/a [46], and the shrinkage rate since the Holocene is about 2-5 mm/a [50]. Based on the time series of baseline length changes in the Junggar Basin during the past 10 years, this paper estimates that the shrinkage rate of the Junggar block is about 3 mm/a, indicating that the shortening deformation of the Junggar block has been relatively stable since Quaternary. On the other hand, judging from the baseline length change rates of those baselines across the Junggar block, the spatiotemporal tectonic deformation of the Junggar block is relatively uniform and stable. The contraction rate of the baseline WUSH-XJWQ across the central Tianshan Mountains (79 •~8 2 • E) is 10.5 mm/a, and the baselines XJKE-XJSH and XJKC-XJDS across the East Tianshan Mountains (82 •~8 7 • E) have a contraction rate of 3~4 mm/a, which is consistent with the convergence rate of the Tianshan Mountains based on GPS velocity field and seismic moment tensors [2,47,48,51]. The shrinkage rate of the Tianshan Mountains as a whole, and the convergence rate between the Tianshan Mountains and the Tarim block all present the characteristics of uneven distribution from east to west, i.e., gradually increasing from east to west.
The Tianshan Mountain is located in the hinterland of central Asia. Its north and south sides are clamped by the Tarim block and Junggar block, respectively. The two sides, including the interior of the Tianshan block, are the main areas where structural deformation occurs [7,28]. The distribution of earthquakes in Figure 7 shows that strong earthquakes mainly occur in the fault zones on the north and south sides of the Tianshan mountain. Especially at the intersection of the western part of the Tarim block, the northern part of the West Kunlun block, and the southern part of the Tianshan Mountains. There have been one earthquake of a magnitude greater than seven and six earthquakes of a magnitude greater than six since 1967. The focal mechanism solutions show that the earthquakes at the north and south sides of the Tianshan Mountains are mainly thrust, and the earthquakes in the Tianshan Mountains are mainly strike-slip, which is consistent with the others research results [4,51]. The nodal plane is approximately parallel to the extension direction of the Tianshan block, and the principal compressive stress axis is almost perpendicular to the Tianshan block, corresponding to the north-south shortening and the east-west extension. Figure 9 shows that the convergence rate of the Tianshan block gradually increases from east to west and reaches 10.5 mm/a in the western region. Based on geological data, Ding et al. [6] deduced that the convergence rate in western Tianshan is about 10 mm/a, which is consistent with the results of this paper.
The Tarim block is not a completely "rigid body", and the area from its middle inner up to the west still shrinks at a rate of about 1-2 mm per year. However, the strong stress that results in the north-south contraction of the Tarim "rigid body" is clearly not caused by the compression of the Pamir Plateau in the west. For a rigid block on a spherical surface, its motion can be described by the rotation around the Euler pole; that is, the rate of motion inside the block increases regularly with the increase of the radius of rotation [4]. Combined with the difference of baseline length change between the east and the west part of the southern Tianshan fault zone, it is inferred that the Qinghai-Tibet Plateau compacts the Tarim block northward and causes it to rotate clockwise, which is one of the leading factors to the deformation of the Tianshan Mountain, and is also the main cause of the occurrence of strong earthquakes and tectonic deformation in this area.
We also obtained the velocities of continuous GPS stations through the analysis of coordinate time series, and displayed them in Figure 11. The image shows the overall trend of the plate movement.  Figure 11 shows that mainland China is characterized by significant regional movements, with the most intense tectonic movements mainly concentrated in the intersection of the Indian and Eurasian plates, i.e., the southern border of the Qinghai-Tibet Plateau. Both Xinjiang and the Qinghai-Tibet Plateau regions move northward. The velocity of the eastern part of the Qinghai-Tibet Plateau and the Sichuan-Yunnan region rotates 90 • clockwise from southeast to southwest. In addition, southwest China has a clockwise rotation centered on the Qinghai-Tibet Plateau, while northwest China has a clockwise rotation trend centered on the Tarim block, which is consistent with others research results [8,24,45]. However, although the GPS velocity field can reflect the overall block motion, it cannot reflect the strain state intuitively. Here, we plot the distribution of the baseline length change rate in mainland China, as shown in Figure 12. In Figure 12, the baseline length change rates are larger in the Tianshan mountain, the Qinghai-Tibet Plateau, and the Sichuan-Yunnan regions, and smaller in the South China and Northeast China regions. In the southern margin of the Qinghai-Tibet Plateau, the baselines mainly exhibit east-west elongation, with the extension rate being smaller near 80 • E and larger at 83 • E and 78 • E, at about 5 mm/a, and consistent with the literature [1]. Figure 12 shows that the rate of northward motion at this margin is about 20-30 mm/a, and the length change rates of the north-south oriented baselines across the Tibetan Plateau are about 5-15 mm/a, which means that a large part of the collision between the Indian plate and the Qinghai-Tibet Plateau is absorbed by the east-west elongation and inter-block interaction within the Tibetan Plateau. The remains of the collision are transmitted to the Tianshan region through the Tarim block. In addition, comparing Figures 5 and 12, it can be seen that both baseline linear strain rates and baseline length change rates have relatively larger values at the block boundaries, and on the contrary, they are smaller inside the blocks. Therefore, in theory, the classification of blocks can be performed by clustering analysis of baseline strain rates or baseline length change rates, when the density of GPS observation stations is high enough.

Conclusions
Based on the GPS observations of CMONOC from 1998 to 2019, the daily coordinate solutions of continuous GPS stations were obtained by using high-precision GPS data processing methods and the Bernese software; baseline linear strain rates of the chosen baselines of CMONOC were subsequently estimated by using the MIDAS method. According to the baseline linear strain rates and/or the baseline length change rates, when the spatial and temporal evolution of crustal deformation in mainland China were studied and compared with the use of the conventional strain analysis method, they generally agreed with each other. However, by using our proposed baseline length change analysis method, with the intentional choice of baselines, the relative motion between tectonic blocks and the spatiotemporal distribution of crustal strain can be analyzed visually and effectively.
Author Contributions: J.W. and X.S. conceived and designed the methodology and experiments; W.W. processed the CMONOC data; X.S. wrote the main manuscript text; G.M. and Y.R. analyzed tectonic movements and earthquakes. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement:
The study did not involve humans or animals.

Informed Consent Statement:
The study did not involve humans or animals.

Data Availability Statement:
The study did not report any data.