GNSS Imaging of Strain Rate Changes and Vertical Crustal Motions over the Tibetan Plateau

: In this paper, we perform a comprehensive analysis of contemporary three-dimensional crustal deformations over the Tibetan Plateau. Considering that the coverage of continuous GNSS sites in the Tibetan Plateau is sparse, a newly designed method that mainly contains Spatial Structure Function (SSF) construction and Median Spatial Filtering (MSF) is adopted to conduct GNSS imaging of point-wise velocities, which can well reveal the spatial pattern of vertical crustal motions. The result illustrates that the Himalayan belt bordering Nepal appears signiﬁcant uplift at the rates of ~3.5 mm/yr, while the low-altitude regions of Nepal and Bhutan near the Tibetan Plateau are undergoing subsidence. The result suggests that the subduction of the Indian plate is the driving force of the uplift and subsidence in the Himalayan belt and its adjacent regions. Similarly, the thrusting of the Tarim Basin is the main factor of the slight uplift and subsidence in the Tianshan Mountains and Tarim Basin, respectively. In addition, we estimate the strain rate changes over the Tibetan Plateau using high-resolution GNSS horizontal velocities. The result indicates that the Himalayan belt and southeastern Tibetan Plateau have accumulated a large amount of strain rate due to the Indian-Eurasian plate collision and blockage of the South China block, respectively.


Introduction
The Tibetan Plateau is the largest inland plateau in China (covering an area of about 2.5 million km 2 ) and the highest plateau in the world (the average elevation reaches 4500 m), which is known as the "roof of the world" and the "third pole" [1][2][3][4]. Due to the Indian-Eurasian plate collision, the Tibetan Plateau is one of the most tectonically active regions in the world, and several large active fault zones are formed inside the Tibetan Plateau [5][6][7]. The tectonic activities of these faults lead to the frequent occurrence of strong earthquakes over the Tibetan Plateau (e.g., 2001 Mw 8.1 Kunlun earthquake, 2010 Mw 7.1 Yushu earthquake, 2013 Ms 7.0 Lushan earthquake, 2015 Mw 7.8 Nepal earthquake, 2020 Mw 6.3 Nima earthquake, and 2021 M 7.4 Maduo earthquake) [8][9][10][11][12][13]. As a continentcontinent collision region that is still undergoing strong tectonic deformations, the Tibetan Plateau has become a natural laboratory for understanding plate collisions and geological evolution [14][15][16]. Thus, characterizing a high-resolution crustal deformation field over the Tibetan Plateau has important implications not only for understanding its tectonic motions and dynamic mechanism but also for seismic hazard assessment in mainland China.
Extensive studies associated with the crustal deformations over the Tibetan Plateau and its surrounding regions were carried out with the aid of space geodetic technologies. Gan et al. [17] utilized the observations of~726 GPS sites to estimate the rigid rotation of the Tibetan Plateau in the Eurasia-fixed reference frame and revealed the deformation features of the interior, southern, northeastern, and eastern parts of the plateau after removing the rigid rotation. Liang et al. [18] established a high-resolution 3D crustal motion field over the Tibetan Plateau using the observations of 750 GPS sites and revealed In order to explore the properties of crustal motions over the Tibetan Plateau, 153 continuous GNSS sites derived from Crustal Movement Observation Network of China (CMONOC) are selected ( Figure 1). The position time series of GNSS sites are provided by the GNSS data product service platform of China Earthquake Administration (ftp: //ftp.cgps.ac.cn/products/position/, accessed on 15 October 2021), and GAMIT/GLOBK 10.4 is utilized to carry out GNSS data processing [24]. The unconstrained daily solutions for satellite orbits and positions of GNSS sites are obtained by GAMIT, and then aligned to the ITRF reference frame using GLOBK. More details about the GNSS data processing can refer to ftp://ftp.cgps.ac.cn/doc/processing_manual.pdf, accessed on 15 October 2021.
Meanwhile, we also adopt the position time series of 52 continuous GNSS sites located in Nepal and Bhutan near the Tibetan Plateau (Figure 1 Nevada Geodetic Laboratory (NGL) (http://geodesy.unr.edu/, accessed on 15 Octob 2021). The time span of all GNSS position time series ends to 2020.7500. The GNSS da processing is conducted by GIPSY-OASIS-II software. NGL analysis center utilizes iono spheric independent combination to remove ionospheric effects, Global Mapping Fun tion (GMF) model to correct tropospheric delays, FES 2004 to correct tidal loading effect and antenna phase center correction parameters refer to the absolute antenna phase cal bration file provided by the International GNSS Service (IGS) center. Figure 1. Spatial distribution of continuous and campaign GNSS sites. Red triangles represent th continuous GNSS sites derived from CMONOC, yellow circles denote the campaign GNSS sit derived from CMONOC, and blue triangles indicate the continuous GNSS sites derived from oth regional networks.
Due to the impact of several strong earthquakes, some sites have recorded significan coseismic and postseismic deformations. Considering this problem, we model the pos tion time series of continuous GNSS sites according to the modeling scheme described Xiang et al. [25], which can simultaneously extract and model the linear, annual, sem annual, coseismic, and postseismic signals [26][27][28]. Figure 2 displays the modeled signa of GNSS position time series for sites XZAR and XZZF. In addition to these continuou GNSS sites, we also collect some horizontal velocities published in recent studies, whic mainly contain the GNSS horizontal velocities published in Kreemer et al. [19] and Wan et al. [23]. The rigid body rotation is utilized to convert these published velocities to coin cide with our solutions, and all GNSS horizontal velocities are relative to the Eurasia plate [29]. Due to the impact of several strong earthquakes, some sites have recorded significant coseismic and postseismic deformations. Considering this problem, we model the position time series of continuous GNSS sites according to the modeling scheme described in Xiang et al. [25], which can simultaneously extract and model the linear, annual, semiannual, coseismic, and postseismic signals [26][27][28]. Figure 2 displays the modeled signals of GNSS position time series for sites XZAR and XZZF. In addition to these continuous GNSS sites, we also collect some horizontal velocities published in recent studies, which mainly contain the GNSS horizontal velocities published in Kreemer et al. [19] and Wang et al. [23]. The rigid body rotation is utilized to convert these published velocities to coincide with our solutions, and all GNSS horizontal velocities are relative to the Eurasian plate [29].

GRACE Datasets
The Gravity Recovery and Climate Experiment (GRACE) and GRACE Follow-On (GFO) are dual-satellite missions monitoring and mapping Earth's changing gravity field. GRACE produced monthly gravity field solutions from April 2002 to June 2017, and GFO continues this record starting in June 2018. Compared with Spherical Harmonic (SH) solutions, mascon solutions can be applied from the regional to global scales. In addition, mascon solutions do not require de-striping and signal restoration, which also have a higher signal-to-noise ratio [30].

GRACE Datasets
The Gravity Recovery and Climate Experiment (GRACE) and GRACE Follow-On (GFO) are dual-satellite missions monitoring and mapping Earth's changing gravity field. GRACE produced monthly gravity field solutions from April 2002 to June 2017, and GFO continues this record starting in June 2018. Compared with Spherical Harmonic (SH) solutions, mascon solutions can be applied from the regional to global scales. In addition, mascon solutions do not require de-striping and signal restoration, which also have a higher signal-to-noise ratio [30].
In order to estimate the trend and annual amplitude induced by surface mass changes over the Tibetan Plateau, the Goddard Space Flight Center (GSFC) mascon solutions provided by GeoForschungs Zentrum (GFZ) Potsdam are adopted in this paper (https://earth.gsfc.nasa.gov/geo/data/grace-mascons, accessed on 15 October 2021). The GSFC mascon solution consists of 41,168 equal-area, 1 arc-degree mascon cells. GGM05C In order to estimate the trend and annual amplitude induced by surface mass changes over the Tibetan Plateau, the Goddard Space Flight Center (GSFC) mascon solutions provided by GeoForschungs Zentrum (GFZ) Potsdam are adopted in this paper (https: //earth.gsfc.nasa.gov/geo/data/grace-mascons, accessed on 15 October 2021). The GSFC mascon solution consists of 41,168 equal-area, 1 arc-degree mascon cells. GGM05C is chosen as the static gravity field model in the solution, the GOT4.7 model of order 90 is adopted as the tidal model, and the RL06 Atmosphere and Ocean De-aliasing Level-1B (AOD1B) product is regarded as the atmospheric and oceanic de-aliasing model [31]. Meanwhile, the Jet Propulsion Laboratory (JPL) TN-13 is utilized to conduct the Geocenter corrections, the TN-14 is utilized to replace the original C20 terms, and the ICE-6G_D is used to correct the Glacial Isostatic Adjustment (GIA) effects. In general, surface mass changes estimated by GRACE/GFO can be expressed as Equivalent Water Height (EWH).

Strain Rates Estimation
In this paper, the vel2strain code developed by the Institute of Seismology, China Earthquake Administration is utilized to estimate the strain rate changes in the Tibetan Plateau. The vel2strain code is written in Fortran, which mainly adopts the method described in Shen et al. [32] and Shen et al. [33]. Shen et al. [32] proposed an algorithm that utilizes an improved least squares method to model strain as a continuous function, which iterates over a two-dimensional space with arbitrarily small increments to ensure the continuity of solution. The horizontal velocity field is expanded to its first derivative at each interpolation position R, which is represented by a rigid block motion (translation and rotation) model and a uniform strain field. The displacement can be expressed by the following formula.
where A is the partial derivative matrix, d is the data vector, ε is the error vector. Among them, m denotes the vector of strain, translation, and rotation, which can be expressed by the following formula in the case of horizontal strain.
where U x and U y are the translation components in the x and y components, τ xx , τ xy , and τ yy are the horizontal strain components, and ω is rotation. The Formula (1) can also be expressed as: where V xi and V yi denote the displacement components of the i-th site located at r i , ∆ xi and ∆ yi denote the components of ∆R i = r i − R. Assuming that C is the covariance matrix of velocities, we can obtain a least-square solution: However, the above solution is the average strain field of the whole region. In order to give more weight to the sites closer to site R, a uniform Gaussian spatial weighting function (G i = L i * Z i ) is adopted. L i and Z i are functions that depend on distance and space coverage. The distance-dependent weighting L i can be defined as: where D is the spatial smoothing parameter. The larger D (i.e., the greater the weight of far-field velocities), the higher the smoothness of strain rate field and the smaller the surface error of strain rate estimation.

GNSS Imaging Scheme
GNSS imaging method is a newly designed method that can generate the image of vertical crustal motions based on point-wise GNSS velocities, which can effectively describe the spatial patterns of vertical crustal motions in continuous space [34][35][36][37][38]. The Spatial Structure Function (SSF) is utilized to describe the correlation between the velocities of GNSS sites, and then the filtering that considers spatial structure is constructed to eliminate outliers and enhance the regional common signals. After that, the velocity field is encrypted according to the image processing technology, and then realize the spatialtemporal features extraction of vertical crustal motions. In this paper, we briefly introduce (1) Spatial structure function SSF groups the GNSS site pairs based on the distances between sites, and builds the mapping relationship between the distances and velocity correlation based on the statistical results of velocity differences, which can reflect the correlation of different GNSS velocities. Firstly, we calculate the absolute value of velocity difference δv ij (i = 1, 2, . . . , n; j = 1, 2, . . . , n) and corresponding distance ∆ ij of each site pair, and sort them according to the distance of site pair. The central angle of large arc formed by site pair is regarded as the distance of site pair, and the absolute value of velocity difference is calculated by: where v i and v j are the velocity of i-th and j-th sites. Secondly, several data pools are divided based on the distance of site pair. In general, 0.25log 10 ∆ is regarded as the separation distance of data pools to ensure that each data pool contains at least 20 site pairs. The MAD of the absolute value of GNSS velocity difference in each data pool is calculated.
Among them, MAD k is the MAD of k-th data pool, δv i is the absolute value of n velocity difference in each data pool, and δv is the median of velocity difference in all data pool.
Finally, we can construct and standardize the SSF according to the following formula: where ss f i is the SSF value of i-th data pool, and MAD i is the MAD of i-th data pool.
(2) Median spatial filtering In this paper, the velocity value of each site is screened by MSF taking into account the spatial structure, which can enhance the commonality of vertical velocities. Firstly, the spherical Delaunay triangulation is utilized to build the spatial relationship between different sites. Then, one GNSS site is selected as the filtering site. The velocity uncertainty and spatial correlation are simultaneously estimated, and the weights of sites connected to the filtering site in the triangulation network are calculated and standardized.
where ss f i is the spatial structure function value obtained by inputting the distance between the current site and estimated site, σ i is the velocity uncertainty of each site, ω i and ω i denote the weights before and after standardization. In order to encrypt the discrete GNSS network, the research region is divided into grid, and each grid point is regarded as a filtering site. According to the velocities and uncertainties of surrounding sites, MSF interpolation is utilized to obtain the velocity of each grid point, and then the GNSS image of vertical crustal motions is generated.
Compared with other methods, GNSS imaging can effectively eliminate the impact of outliers in velocities, which can accurately describe common features of vertical crustal motions.

GNSS Horizontal Velocity Field
By combining with previous studies, a high-resolution GNSS horizontal velocity field over the Tibetan Plateau and its surrounding regions is established (Figure 3), which contains 2032 GNSS velocities. The GNSS velocity field clearly describes the crustal motions over the Tibetan Plateau. The Indian plate is pushing the Eurasian plate in the NE direction, and the velocities decrease sharply along the NE direction. Meanwhile, the motions of the Tibetan Plateau are blocked by the South China block, resulting in the significant clockwise rotation of the GNSS velocity field in the southeastern Tibetan Plateau. In addition, the GNSS velocity field also appears weak clockwise rotation in the Tarim Basin due to the thrusting of the Tarim Basin. The spatial patterns of the GNSS velocity field reveal the obvious lateral extrusion in the Tibetan Plateau, which is consistent with the conclusions of previous studies [18,21,23].
where is the spatial structure function value obtained by inputting the distance between the current site and estimated site, is the velocity uncertainty of each site, and denote the weights before and after standardization. In order to encrypt the discrete GNSS network, the research region is divided into grid, and each grid point is regarded as a filtering site. According to the velocities and uncertainties of surrounding sites, MSF interpolation is utilized to obtain the velocity of each grid point, and then the GNSS image of vertical crustal motions is generated. Compared with other methods, GNSS imaging can effectively eliminate the impact of outliers in velocities, which can accurately describe common features of vertical crustal motions.

GNSS Horizontal Velocity Field
By combining with previous studies, a high-resolution GNSS horizontal velocity field over the Tibetan Plateau and its surrounding regions is established (Figure 3), which contains 2032 GNSS velocities. The GNSS velocity field clearly describes the crustal motions over the Tibetan Plateau. The Indian plate is pushing the Eurasian plate in the NE direction, and the velocities decrease sharply along the NE direction. Meanwhile, the motions of the Tibetan Plateau are blocked by the South China block, resulting in the significant clockwise rotation of the GNSS velocity field in the southeastern Tibetan Plateau. In addition, the GNSS velocity field also appears weak clockwise rotation in the Tarim Basin due to the thrusting of the Tarim Basin. The spatial patterns of the GNSS velocity field reveal the obvious lateral extrusion in the Tibetan Plateau, which is consistent with the conclusions of previous studies [18,21,23]. . GNSS horizontal velocity field with respect to Eurasian plate over the Tibetan Plateau. Red arrows denote our solutions, green arrows denote the velocities derived from Wang et al. [23], and gray arrows denote the velocities derived from Kreemer et al. [19].
Meanwhile, we also display the spatial distribution of the GNSS velocity field in the EW and NS directions, and the statistical results of velocities and uncertainties ( Figure 4). Most GNSS velocity uncertainties are less than 1 mm/yr, indicating that the estimated GNSS velocity field is reliable. The rate of GNSS velocity ranges from −8.4 to 25.6 mm/yr in the EW direction, and it varies from −15.1 to 39.4 mm/yr in the NS direction. For the Figure 3. GNSS horizontal velocity field with respect to Eurasian plate over the Tibetan Plateau. Red arrows denote our solutions, green arrows denote the velocities derived from Wang et al. [23], and gray arrows denote the velocities derived from Kreemer et al. [19].
Meanwhile, we also display the spatial distribution of the GNSS velocity field in the EW and NS directions, and the statistical results of velocities and uncertainties ( Figure 4). Most GNSS velocity uncertainties are less than 1 mm/yr, indicating that the estimated GNSS velocity field is reliable. The rate of GNSS velocity ranges from −8.4 to 25.6 mm/yr in the EW direction, and it varies from −15.1 to 39.4 mm/yr in the NS direction. For the east component of GNSS velocity, the large rates are concentrated in the southeastern Tibetan Plateau, which coincides with the clockwise rotation of the GNSS velocity field in this region. As for the north component of the GNSS velocity, the rates are obviously reduced along the NE direction, suggesting that the Tibetan Plateau and its surroundings have produced significant crustal shortening along the NE direction due to the collision of the Indian and Eurasian plates. east component of GNSS velocity, the large rates are concentrated in the southeastern Tibetan Plateau, which coincides with the clockwise rotation of the GNSS velocity field in this region. As for the north component of the GNSS velocity, the rates are obviously reduced along the NE direction, suggesting that the Tibetan Plateau and its surroundings have produced significant crustal shortening along the NE direction due to the collision of the Indian and Eurasian plates.

Strain Rate Changes
In order to explore the tectonic activities over the Tibetan Plateau, the strain rate changes in this region are characterized based on the high-resolution GNSS velocity field. The region is divided into a grid by 1° of latitude and longitude, and the strain rate of each grid point is estimated by the method described in Section 2.4.
As depicted in Figure 5a, a large amplitude of principal strain rate is accumulated in the Himalaya arc tectonic belt, including Nepal, Bhutan, northwestern and northeastern India, which is dominated by the extrusion strain rate. The Himalaya tectonic belt is the foreland of the Indian-Eurasian plate collision, and rapid plate convergence leads to significant strain rate accumulation. The principal strain rate in the southeastern Tibetan Plateau is also dominated by the extrusion strain rate, mainly because the motion of the Tibetan Plateau is blocked by the South China block. In addition, the principal strain rate in the Tianshan belt is mainly manifested as extrusion strain rate, implying that the thrusting of Tarim Basin has induced the extrusion deformations in the Tianshan belt to a certain

Strain Rate Changes
In order to explore the tectonic activities over the Tibetan Plateau, the strain rate changes in this region are characterized based on the high-resolution GNSS velocity field. The region is divided into a grid by 1 • of latitude and longitude, and the strain rate of each grid point is estimated by the method described in Section 2.4.
As depicted in Figure 5a, a large amplitude of principal strain rate is accumulated in the Himalaya arc tectonic belt, including Nepal, Bhutan, northwestern and northeastern India, which is dominated by the extrusion strain rate. The Himalaya tectonic belt is the foreland of the Indian-Eurasian plate collision, and rapid plate convergence leads to significant strain rate accumulation. The principal strain rate in the southeastern Tibetan Plateau is also dominated by the extrusion strain rate, mainly because the motion of the Tibetan Plateau is blocked by the South China block. In addition, the principal strain rate in the Tianshan belt is mainly manifested as extrusion strain rate, implying that the thrusting of Tarim Basin has induced the extrusion deformations in the Tianshan belt to a certain extent. The amplitude of principal strain rate in the Sichuan Basin and Tarim Basin is relatively small (within 10 nanostrain/yr (i.e., 10 −8 /yr)), indicating that these regions are relatively stable without too much strain rate accumulation.
regions of Myanmar and Bangladesh) (see the background color in Figure 4a). The pea of maximum shear strain rate is located at the Himalayan tectonic belt, and the maximum reaches 128.5 nanostrain/yr. The second is the eastern Tibetan Plateau and some region of Myanmar and Bangladesh, and the maximum shear strain rate values are all above 8 nanostrain/yr. The maximum shear strain rate inside the Sichuan Basin and Tarim Basi is close to 0 nanostain/yr, which also suggests that the interior of the Sichuan Basin an Tarim Basin is relatively stable.  For the maximum shear strain rate, the large value is accumulated in four regions (i.e., Himalaya arc tectonic belt, eastern Tibetan Plateau, Tianshan Mountains, and some regions of Myanmar and Bangladesh) (see the background color in Figure 4a). The peak of maximum shear strain rate is located at the Himalayan tectonic belt, and the maximum reaches 128.5 nanostrain/yr. The second is the eastern Tibetan Plateau and some regions of Myanmar and Bangladesh, and the maximum shear strain rate values are all above 80 nanostrain/yr. The maximum shear strain rate inside the Sichuan Basin and Tarim Basin is close to 0 nanostain/yr, which also suggests that the interior of the Sichuan Basin and Tarim Basin is relatively stable.
The dilatational strain rate over the Tibetan Plateau varies from −112.3 to 48.2 nanostrain/yr with the mean of −10.5 nanostrain/yr (Figure 5b). Similarly, the large-scale negative dilatational strain rate is focused on the Himalayan arc tectonic belt, suggesting that this region has undergone crustal extrusion deformation. Meanwhile, some small-scale negative dilatational strain rate appears in the Tianshan belt and eastern Tibetan Plateau. The dilatational strain rates inside the Tibetan Plateau at the longitude of 87 • to 95 • are positive, implying that the crust inside the Tibetan Plateau has experienced extension deformation. The junction of Tibetan Plateau and South China block (e.g., some regions of Yunnan and Sichuan provinces) also occurs large-scale positive dilatational strain rate, which can be attributed to the extrusion deformation induced by the blockage of South China block.
The above results demonstrate that the Himalayan orogenic belt has accumulated a large amount of strain rate due to the collision of the Indian plate and Eurasian plate, and the tectonic motions of the entire Himalayan arc belt are quite active, which still belongs to the region with seismic hazard risk. In addition, since the motions of the Tibetan Plateau are blocked by the South China block, a large amount of strain rate has also accumulated at the junction, which is also one of the tectonically active regions in mainland China.

GNSS Imaging of Vertical Crustal Motions
Up to now, few studies have released the GNSS vertical velocity field of the Tibetan Plateau, and most of them explore the vertical crustal motion over the Tibetan Plateau based on sparse continuous GNSS sites, which have relatively abundant observations. Liang et al. [18] had released the GNSS vertical velocity field of the Tibetan Plateau on basis of dense campaign GNSS sites. However, only two phases of campaign GNSS observations are utilized to estimate the vertical velocity, which cannot consider the effects of significant seasonal oscillations induced by surface mass loading. Thus, we adopt the GNSS imaging method to generate an image of vertical velocity based on the discrete-continuous GNSS network and then reveal the spatial patterns of vertical crustal motions over the Tibetan Plateau.
Firstly, we construct the Spatial Structure Function (SSF) to contain information representing the vertical signal variance, which can be attributed to the distance (i.e., ∆) between GNSS sites. In general, GNSS sites close to each other will exhibit more similar crustal motion features. For each pair of GNSS sites, the absolute value of vertical velocity difference is calculated, and Figure 6b presents the distribution of differences along log 10 (∆). According to Formulas (8) and (9), the SSF is constructed and normalized, and the value of SSF as a function of distance is depicted in Figure 6c. The value of SSF decreases monotonically as the distance between GNSS sites increases, and it reaches 0 when the distance between sites is around 13.34 km.
Secondly, the Median Spatial Filtering (MSF) is adopted to perform spatial filtering for GNSS vertical velocities, which can reduce the effect of outliers that are obviously different from neighboring velocity rates. MSF utilizes the spherical Delaunay triangulation to connect the evaluation point (i.e., one GNSS site that the spatial filtering is performed) and its nearest neighboring sites, and the weighted median of neighboring velocity rates are computed according to Formulas (10) and (11). Then, the weighted median of the velocities of neighboring sites can be regarded as the velocity rate of the evaluation point. MSF can reduce the importance of sites with large uncertainty and increase the importance of sites near the evaluation point. Figure 7a,b present the original and filtered GNSS vertical velocities in the Tibetan Plateau, and the differences of original and filtered velocities are depicted in Figure 7c. This result indicates that the outliers in the original GNSS velocities are effectively removed, and thus highlighting the common crustal motion signal supported by multiple sites. Secondly, the Median Spatial Filtering (MSF) is adopted to perform spatial filtering for GNSS vertical velocities, which can reduce the effect of outliers that are obviously different from neighboring velocity rates. MSF utilizes the spherical Delaunay triangulation to connect the evaluation point (i.e., one GNSS site that the spatial filtering is performed) and its nearest neighboring sites, and the weighted median of neighboring velocity rates are computed according to Formulas (10) and (11). Then, the weighted median of the velocities of neighboring sites can be regarded as the velocity rate of the evaluation point. MSF can reduce the importance of sites with large uncertainty and increase the importance of sites near the evaluation point. Figure 7a,b present the original and filtered GNSS vertical velocities in the Tibetan Plateau, and the differences of original and filtered velocities are depicted in Figure 7c. This result indicates that the outliers in the original GNSS velocities are effectively removed, and thus highlighting the common crustal motion signal supported by multiple sites.  Finally, we carry out GNSS imaging of point-wise velocities to generate an image of vertical crustal motions over the Tibetan Plateau. This region is divided into a grid by 0.5° of latitude and longitude. Similarly, we establish the spatial structure relationship be- Finally, we carry out GNSS imaging of point-wise velocities to generate an image of vertical crustal motions over the Tibetan Plateau. This region is divided into a grid by 0.5 • of latitude and longitude. Similarly, we establish the spatial structure relationship between each grid point and neighboring sites using spherical Delaunay triangulation. The velocity of each grid point can be replaced by the weighted median of the velocities of neighboring sites. Figure 8a,c display the GNSS imaging of vertical crustal motions using the filtered velocities and velocities without median spatial filtering applied, respectively. Figure 8e presents the GNSS imaging of vertical crustal motion using the filtered velocities, while each grid is replaced by the weight mean of the velocities of neighboring sites. Compared with the image generated by the filtered GPS velocities, the image of vertical crustal motions generated by the unfiltered velocities is more susceptible to the outliers of velocity rate. Meanwhile, if each grid is replaced by the weighted mean of the velocities of neighboring sites, the estimated image of vertical crustal motions tends to be rough and unreliable. As described in Figure 9, the spatial patterns of vertical crustal motions over the Tibetan Plateau are diverse, which may be associated with the complex geological, topographical, hydrological, and climate features of the Tibetan Plateau. The large-scale crustal uplift is focused on the southern Tibetan Plateau, especially for the Himalaya belt bordering Nepal, where the uplift rates are about 3.5 mm/yr. However, some regions of Nepal and Bhutan near the Tibetan Plateau experience remarkable crustal subsidence as altitude decreases, suggesting that the significant crustal subduction appears under the Himalayas. Apart from these regions, there are slight crustal uplifts in the Tianshan belt and As described in Figure 9, the spatial patterns of vertical crustal motions over the Tibetan Plateau are diverse, which may be associated with the complex geological, topographical, hydrological, and climate features of the Tibetan Plateau. The large-scale crustal uplift is focused on the southern Tibetan Plateau, especially for the Himalaya belt border-Remote Sens. 2021, 13, 4937 13 of 20 ing Nepal, where the uplift rates are about 3.5 mm/yr. However, some regions of Nepal and Bhutan near the Tibetan Plateau experience remarkable crustal subsidence as altitude decreases, suggesting that the significant crustal subduction appears under the Himalayas. Apart from these regions, there are slight crustal uplifts in the Tianshan belt and northern Tibetan Plateau, and the Tarim Basin occurs slight subsidence. In addition, southwest China (i.e., some regions of Yunan and Sichuan provinces) weak crustal subsidence with rates of 0 to 1 mm/yr also appear, which may be related to the gain of hydrological mass induced by the slight increase in annual mean precipitation. Moreover, we extract seven profiles throughout the Tibetan Plateau along AA1, BB1, CC1, DD1, EE1, FF1, and GG1 ( Figure 10) in order to better reveal regional vertical crustal motions. The first three profiles along AA1, BB1, and CC1 clearly exhibit the vertical crustal motions across the southern Himalayas. The vertical crustal motions are dominated by subsidence at the low-altitude regions of Nepal and Bhutan outside the Tibetan Plateau, and rapid crustal uplift with the rates of 1-5 mm/yr is found in the Himalaya belt beyond the subduction of the Indian plate. As latitude increases, the vertical crustal motions are followed by subsidence in the Tarim Basin and slight uplift in the Tianshan Mountains. The three profiles along EE1, FF1, and GG1 pass through the Tibetan Plateau and its surrounding regions at the latitude of 37.5°, 32.5°, and 29°. Similarly, slight subsidence occurs in the Tarim Basin and southwest China in comparison with the Tibetan Plateau. Moreover, we extract seven profiles throughout the Tibetan Plateau along AA1, BB1, CC1, DD1, EE1, FF1, and GG1 ( Figure 10) in order to better reveal regional vertical crustal motions. The first three profiles along AA1, BB1, and CC1 clearly exhibit the vertical crustal motions across the southern Himalayas. The vertical crustal motions are dominated by subsidence at the low-altitude regions of Nepal and Bhutan outside the Tibetan Plateau, and rapid crustal uplift with the rates of 1-5 mm/yr is found in the Himalaya belt beyond the subduction of the Indian plate. As latitude increases, the vertical crustal motions are followed by subsidence in the Tarim Basin and slight uplift in the Tianshan Mountains. The three profiles along EE1, FF1, and GG1 pass through the Tibetan Plateau and its surrounding regions at the latitude of 37.5 • , 32.5 • , and 29 • . Similarly, slight subsidence occurs in the Tarim Basin and southwest China in comparison with the Tibetan Plateau. 7 15 of 21 Figure 10. Seven profiles (a-g) of GNSS-imaged vertical velocities and topography changes throughout the Tibetan plateau along AA1, BB1, CC1, DD1, EE1, FF1, and GG1 shown in (h). The gray area denotes the topography changes, and the shaded area of velocity represents the uncertainties in each profile.

Surface Mass Changes over the Tibetan Plateau
In addition to the seasonal oscillations induced by surface mass loading, the longterm gain and loss of surface mass can also result in the crustal uplift and subsidence [39][40][41][42]. Thus, aiming to better explore the properties of vertical crustal motions over the Tibetan Plateau, we also determine the trends and annual amplitudes of surface mass changes based on GRACE mascon solutions. The values of surface mass change trends over the Tibetan Plateau range from -6.42 to 2.25 cm/yr, and the annual amplitudes of surface mass changes vary from 0 to 33.2 cm/yr. As depicted in Figure 11, three regions that exhibit negative surface mass changes are northwest India, the Himalayas, and Tianshan. Among them, the anomalous surface mass loss in northwest India may be associated with the water resource changes induced by groundwater over-exploitation, and the surface mass loss in the Himalayas and Tianshan is mainly attributed to the glacier ice melting.
The surface mass changes over the northern and central Tibetan Plateau, Tarim Basin, Figure 10. Seven profiles (a-g) of GNSS-imaged vertical velocities and topography changes throughout the Tibetan plateau along AA1, BB1, CC1, DD1, EE1, FF1, and GG1 shown in (h). The gray area denotes the topography changes, and the shaded area of velocity represents the uncertainties in each profile.

Surface Mass Changes over the Tibetan Plateau
In addition to the seasonal oscillations induced by surface mass loading, the long-term gain and loss of surface mass can also result in the crustal uplift and subsidence [39][40][41][42]. Thus, aiming to better explore the properties of vertical crustal motions over the Tibetan Plateau, we also determine the trends and annual amplitudes of surface mass changes based on GRACE mascon solutions. The values of surface mass change trends over the Tibetan Plateau range from −6.42 to 2.25 cm/yr, and the annual amplitudes of surface mass changes vary from 0 to 33.2 cm/yr. As depicted in Figure 11, three regions that exhibit negative surface mass changes are northwest India, the Himalayas, and Tianshan. Among them, the anomalous surface mass loss in northwest India may be associated with the water resource changes induced by groundwater over-exploitation, and the surface mass loss in the Himalayas and Tianshan is mainly attributed to the glacier ice melting.

Tectonic Process of Tibetan Plateau
In this paper, we construct a high-resolution GNSS horizontal velocity field by combining our results and the latest studies, which clearly describe the spatial patterns of crustal motions over the Tibetan Plateau and its surroundings. Meanwhile, we also display the spatial distribution of the east and north components of GNSS horizontal velocities, and it can reveal the changes of horizontal crustal motions more intuitively.
It can be found the Indian plate is pushing the Eurasian plate in the NE direction, and the velocities decrease sharply from 40 to 10 mm/yr along the NE direction, revealing the continuous crustal shortening over the Tibetan Plateau due to the pushing of the Indian plate. In the interior of the Tibetan Plateau, GNSS horizontal velocity starts to "escape" laterally along the east and west directions. Due to the blockage of the South China block, the southeastern Tibetan Plateau occurs in an obvious clockwise rotation. The Tarim Basim displays weak clockwise rotation, and the Tianshan belt also appears crustal shortening due to the thrusting of the Tarim Basin. Many studies have made efforts to reveal The surface mass changes over the northern and central Tibetan Plateau, Tarim Basin, and southwest China (e.g., Yunan and Sichuan provinces) occur large-scale positive trends. One potential reason for the positive trends is the gain of hydrological mass induced by the slight increase in annual mean precipitation in these regions. Unlike the surface mass change trends, the large values of the annual amplitudes of surface mass changes are mainly concentrated in low-latitude coastal regions and gradually decrease as the latitude increases. In general, the annual amplitude can reveal the effects of surface mass loading, and this result suggests that the effects of surface mass loading in coastal regions are obviously stronger than in other regions. Combined with the image of vertical crustal motions, it can be found that the surface mass changes can only partly explain the vertical crustal motions of the Tibetan Plateau. Thus, for most regions of the Tibetan Plateau, the tectonic motions, surface mass changes, and other geophysical factors jointly lead to the crustal uplift and subsidence.

Tectonic Process of Tibetan Plateau
In this paper, we construct a high-resolution GNSS horizontal velocity field by combining our results and the latest studies, which clearly describe the spatial patterns of crustal motions over the Tibetan Plateau and its surroundings. Meanwhile, we also display the spatial distribution of the east and north components of GNSS horizontal velocities, and it can reveal the changes of horizontal crustal motions more intuitively.
It can be found the Indian plate is pushing the Eurasian plate in the NE direction, and the velocities decrease sharply from 40 to 10 mm/yr along the NE direction, revealing the continuous crustal shortening over the Tibetan Plateau due to the pushing of the Indian plate. In the interior of the Tibetan Plateau, GNSS horizontal velocity starts to "escape" laterally along the east and west directions. Due to the blockage of the South China block, the southeastern Tibetan Plateau occurs in an obvious clockwise rotation. The Tarim Basim displays weak clockwise rotation, and the Tianshan belt also appears crustal shortening due to the thrusting of the Tarim Basin. Many studies have made efforts to reveal the patterns of crustal motions over the Tibetan Plateau [8,17,18,21], and our results are consistent with previous studies.
Based on the GNSS velocity field, the strain rate changes over the Tibetan Plateau and its surrounding regions are estimated. In total, four regions that have accumulated a large amount of strain rate are the Himalaya arc belt, Tianshan belt, southeastern Tibetan Plateau, and some regions of Myanmar and Bangladesh. The accumulation of strain rate in the Himalaya arc belt and Tianshan belt can be attributed to the Indian-Eurasian collision and the thrusting of Tarim Basin, respectively. Meanwhile, the strain rate accumulation in the southeastern Tibetan Plateau is due to the fact the crustal motions of the Tibetan Plateau are blocked by the South China block. Combined with the spatial patterns of GNSS velocity, it can be found that tectonic convergence is the dominant factor that leads to the crustal motions, crustal shortening, and strain rate accumulation over the Tibetan Plateau. The strain rate accumulation induced by tectonic convergence suggests that the crustal activities in these regions are relatively active, which belong to seismically active belts in the Tibetan Plateau. Meanwhile, some regions (i.e., most regions inside the Tibetan Plateau, Tarim Basin, and Sichuan Basin) have less strain rate accumulation imply that the interior of these regions is relatively stable.
Up to now, many studies have revealed the crustal motions and strain rate changes over the Tibetan Plateau. At first, due to the lack of GNSS observations, the crustal motions in many regions of the Tibetan Plateau cannot be described in detail, especially for northeastern India, eastern Himalaya, and Tianshan Mountains. Meanwhile, some studies focus on the eastern Tibetan Plateau with dense coverage of GNSS sites, which can only partly reveal the crustal deformation features of the Tibetan Plateau [43][44][45]. With the abundance of GNSS observations, the crustal deformations of the Tibetan Plateau can be better explored (e.g., [46]). On basis of previous studies, we collected a high-resolution GNSS velocity field over the Tibetan Plateau and its surroundings. Compared with previous studies, the crustal motions and strain rate changes of the Tibetan Plateau can be well described based on the dense GNSS velocity field, especially for the margin of the Tibetan Plateau bordering India, which is important for understanding the Indian-Eurasian plate collision. Whereas, GNSS observations are still lacking in the southwestern margin and interior of the Tibetan Plateau, especially for the no man's land inside the Tibetan Plateau, limiting further understanding of the geodynamic mechanisms of tectonic motions in the Tibetan Plateau. Thus, the tectonic process over the Tibetan Plateau needs to be further explored with the aid of richer observations.

Spatial Patterns of Vertical Crustal Motions in the Tibetan Plateau
According to the GNSS imaging scheme, we carry out imaging of the point-wise velocities in the Tibetan Plateau and construct an image of vertical crustal motions using the sparse continuous GNSS sites. Compared with the discrete vertical velocity field published in previous studies [18,21,22], the GNSS-derived image can better reveal the spatial patterns of vertical crustal motions. For the Tibetan Plateau and its surroundings, the southern Tibetan Plateau and Tianshan Mountains occur large-scale uplift, and three regions that exhibit slight subsidence are Tarim Basin, southwest China (i.e., some regions of Yunan and Sichuan provinces), and the low-altitude region of Nepal and Bhutan near the Tibetan Plateau. The details of some subsidence and uplift signals over the Tibetan Plateau are discussed below.
The southern Tibetan Plateau is situated at the collision of the Indian and Eurasian plates, and the Himalayan tectonic belt is the foreland of the Indian-Eurasian plate collision. The image of vertical crustal motions illustrates that the Himalayan tectonic belt bordering Nepal appears significant crustal uplift at the rates of~3.5 mm/yr, and the other regions of the southern Tibetan Plateau are uplifted at the rates of~2 mm/yr. However, the surface mass changes in the Himalayan tectonic belt only occur with slight negative trends (Figure 11), which cannot explain the significant crustal uplift. Meanwhile, the low altitude regions of Nepal and Bhutan near the Tibetan Plateau are undergoing obvious crustal subsidence at the rate of~3 mm/yr. These results jointly indicate that the subduction of the Indian plate is the dominant factor of the crustal uplift and subsidence in the Himalayas and its adjacent regions. Thus, the rapid crustal uplift is observed in the Himalayan belt beyond the subduction of the Indian plate.
According to the image of vertical crustal motions, the slight subsidence with the rates of 0-1.5 mm/yr is found in the Tarim Basin, most regions of Yunan province, and some regions of Sichuan province. The distribution of surface mass changes indicates that southwest China occurs large-scale positive trends at the rates of 0-2 mm/yr (Figure 11), which is mainly induced by the gain of hydrological mass induced by the continuous, slight increase in annual mean precipitation in these regions. The distribution of strain rate changes demonstrates that southwest China is relatively stable without too much strain rate accumulation. Therefore, the slight crustal subsidence in southwest China can be attributed to the gain of surface mass induced by the annual mean precipitation. Some studies have explored the vertical crustal motions over the southeastern Tibetan Plateau based on different geodetic observations (e.g., [47]), which is consistent with our conclusion.
In addition, the slight crustal uplift with the rates of 0-1 mm/yr is also found in the Tianshan belt. For the uplift and subsidence of the Tianshan belt and Tarim Basin, it is more likely to be related to tectonic motions. Zhang et al. [48] explore the properties of vertical crustal motions in Tianshan and its surroundings using 31 GNSS sites and suggested that the crustal uplift of Tianshan and subsidence of the Tarim Basin are attributed to the convergence of the Tarim Basin under the pushing of the Indian plate. The Tianshan belt is simultaneously affected by the far-field effects of the Indian-Eurasian collision and the convergence of the Tarim Basin and the Jungger block, which results in frequent tectonic activity in this region [49][50][51]. The thrusting of the Tarim Basin has an important impact on the growth of the Tianshan belt. The image of vertical crustal motions shows that the weak uplift with the rates of~2 mm/yr in the Tianshan belt, and the Tarim Basin is subsiding at the rates of~1 mm/yr, which coincides with the thrusting process. Thus, the thrusting of the Tarim Basin induced by the pushing of the Indian plate is the main factor of the crustal subsidence in the Tarim Basin and the crustal uplift in the Tianshan belt.

Conclusions
In this paper, we construct a high-resolution GNSS horizontal velocity field over the Tibetan Plateau by combining it with previous studies. The high-resolution GNSS horizontal velocity field clearly reveals the crustal motion over the Tibetan Plateau. The Indian plate is pushing the Eurasian plate in a NE direction, and the velocities decrease sharply along the NE direction. Based on the dense GNSS horizontal velocity, we determine the strain rate changes in the Tibetan Plateau and its surrounding regions, including the changes of principal strain rate, maximum shear strain rate, and dilatational strain rate. Two regions that have accumulated a large amount of strain rate are the Indian-Eurasian collision and the junction of Tibetan Plateau and South China block. Moreover, considering that campaign GNSS sites do not have enough observations to determine reliable vertical motion rates, the GNSS imaging method is adopted to carry out imaging of point-wise velocities. According to the image of vertical crustal motions, the spatial patterns of vertical crustal motions over the Tibetan Plateau are explored. Compared with other regions, the Himalayan arc tectonic belt bordering Nepal has significant crustal uplift, which is mainly attributed to the subduction of the Indian plate.