The Crustal Dynamics and Its Geological Explanation of the Three-Dimensional Co-Seismic Deformation Field for the 2021 Maduo MS7.4 Earthquake Based on GNSS and InSAR

Three-dimensional deformation is an important input to explore seismic mechanisms and geodynamics. The GNSS and InSAR technologies are commonly used to obtain the co-seismic three-dimensional deformation field. This paper focused on the effect of calculation accuracy caused by the deformation correlation between the reference point and the points involved in the solution, to build a high-accuracy three-dimensional deformation field for a detailed geological explanation. Based on the variance component estimation (VCE) method, the InSAR LOS, azimuthal deformation, and the GNSS horizontal and vertical deformation were integrated to solve the three-dimensional displacement of the study area in combination with the elasticity theory. The accuracy of the three-dimensional co-seismic deformation field of the 2021 Maduo MS7.4 earthquake obtained by the method proposed in this paper, was compared with that obtained from the only InSAR measurements obtained using a multi-satellite and multi-technology approach. The results showed the difference in root-mean-square errors (RMSE) of the integration and GNSS displacement was 0.98 cm, 5.64 cm, and 1.37 cm in the east–west, north–south and vertical direction respectively, which was better than the RMSE of the method using only InSAR and GNSS displacement, which was 5.2 cm and 12.2 cm in the east–west, north–south, and no vertical direction. With the geological field survey and aftershocks relocation, the results showed good agreement with the strike and the position of the surface rupture. The maximum slip displacement was about 4 m, which was consistent with the result of the empirical statistical formula. It was firstly found that the pre-existing fault controlled the vertical deformation on the south side of the west end of the main surface rupture caused by the Maduo MS7.4 earthquake, which provided the direct evidence for the theoretical hypothesis that large earthquakes could not only produce surface rupture on seismogenic faults, but also trigger pre-existing faults or new faults to produce surface rupture or weak deformation in areas far from seismogenic faults. An adaptive method was proposed in GNSS and InSAR integration, which could take into account the correlation distance and the efficiency of homogeneous point selection. Meanwhile, deformation information of the decoherent region could be recovered without interpolation of the GNSS displacement. This series of findings formed an essential supplement to the field surface rupture survey and provided a novel idea for the combination of the various spatial measurement technologies to improve the seismic deformation monitoring.


Introduction
The three-dimensional deformation field caused by any geological hazard is an important parameter in the study of crustal movement characteristics, which implies certain geodynamics. Accurately acquiring the complete regional three-dimensional crustal deformation field is an important research direction for the dynamics of regional crustal deformation and understanding its geohazard process. The rapid development of geodetic measurement technologies, such as InSAR (Interferometric Synthetic Aperture Radar) and GNSS (Global Navigation Satellite System), provides an opportunity to acquire an accurate regional three-dimensional crustal deformation field. The InSAR has the advantages of all-weather, all-day, high-spatial resolution, while it only provides one-dimensional deformation information, so that the monitoring of true crustal deformation dynamics is restricted. The GNSS can provide high-precision three-dimensional deformation information [1], but its spatial resolution is influenced by the density of GNSS stations so that the obtained deformations are discrete, which also limits the study of regional crustal dynamics.
Integrating the advantages of these two observations to obtain high-accuracy and highspatial-resolution three-dimensional deformation, is a subject that has been explored. Since Guemundsson [2] proposed the concept of combining GPS and INSAR to study crustal deformation, the studies [3][4][5][6][7][8][9][10][11][12][13] have focused on optimizing the Gibbs energy equation for acquiring three-dimensional deformation. They required interpolating GNSS data to the pixel positions of InSAR images [14] and establishing nonlinear equations without considering the correlation between each adjacent point. However, the interpolation itself had errors [15], which would become greater, especially when the abrupt deformation happened, such as earthquakes, volcanoes, and so on. It is not suitable for high accuracy deformation monitoring.
Guglielmino [16] proposed a non-interpolation method for reconstructing threedimensional deformation fields by integrating GNSS and InSAR, namely, the SISTEM (simultaneous and integrated strain tensor estimation from geodetic and satellite deformation measurements), which took into account the spatial correlation between adjacent points, and the linear relationship between strain and displacement based on the elasticity theory. Guglielmino [16][17][18] solved the area of interest three-dimensional displacement field by the weighted least square (WLS) adjustment by using the GNSS discrete three-dimensional displacements around a reference point and the DInSAR line of sight (LOS) direction displacement at the reference point, which did not consider the correlation between the adjacent points and the reference point. This method could not integrate InSAR and GNSS data in the decoherence regions due to the absence of InSAR data. Luo et al. [19,20] increased the InSAR LOS direction displacement around the reference point and developed the ESISTEM (extended simultaneous and integrated strain tensor estimation from geodetic and satellite deformation measurements) method, which used all GNSS horizontal displacement. The above-mentioned methods did not consider the deformation correlation between the GNSS stations and the reference point. Furthermore, the vertical deformation accuracy of GNSS measurement reached the millimeter level, and InSAR azimuthal direction displacement was sensitive to the north-south deformation, which could contribute to a constraint on the accuracy of the three-dimensional deformation.
The high construction cost and strict setting requirements determine the uneven distribution density of GNSS stations. In many cases, only InSAR was used to obtain threedimensional deformation [21][22][23][24][25][26][27][28][29][30] by integrating displacements obtained from different orbits (ascending and descending) and different InSAR techniques, such as the Differential Interferometric Synthetic Aperture Radar (DInSAR) [31], the multiple aperture interferometry (MAI) [32], the pixel offset-tracking (POT) [33], the burst overlap interferometry (BOI) [34]. However, the InSAR measurements were of relative displacement, which should be constrained by GNSS absolute displacement to obtain the absolute three-dimensional deformation in higher accuracy. is the specific tectonic location and active fault of the black box in (a): the green-yellow font represents the abbreviated name of the secondary tectonic block, the light blue font represents the abbreviated name of the main fault, the cyan line represents the primary block boundary, the blue line represents the secondary tectonic block boundary, and the black line is the active fault divided by Deng et al. [44]. The yellow pentagram represents the epicenter location of the 2021 Maduo MS7.4 earthquake according to the China earthquake networks center (CENC), the green triangle represents the GNSS stations near the area of interest, which shows some stations overlap because of close ranges, the purple dotted frame is the area of interest, and the red line is the surface rupture caused by the earthquake, the light blue square and dot are city and town, respectively. So far, the seismogenic structure and deformation characteristics of the Maduo MS7.4 earthquake had been researched by various means [45][46][47][48][49][50][51][52][53], which showed that this is the specific tectonic location and active fault of the black box in (a): the green-yellow font represents the abbreviated name of the secondary tectonic block, the light blue font represents the abbreviated name of the main fault, the cyan line represents the primary block boundary, the blue line represents the secondary tectonic block boundary, and the black line is the active fault divided by Deng et al. [44]. The yellow pentagram represents the epicenter location of the 2021 Maduo MS7.4 earthquake according to the China earthquake networks center (CENC), the green triangle represents the GNSS stations near the area of interest, which shows some stations overlap because of close ranges, the purple dotted frame is the area of interest, and the red line is the surface rupture caused by the earthquake, the light blue square and dot are city and town, respectively. So far, the seismogenic structure and deformation characteristics of the Maduo M S 7.4 earthquake had been researched by various means [45][46][47][48][49][50][51][52][53], which showed that this earthquake was a special and complex left-lateral slip earthquake. Unlike the previous recognition that the main surface rupture was located at the epicenter, no obvious surface rupture had been found near the epicenter according to the field investigation; the obvious surface rupture was far away from the epicenter, and many fissures and extrusion bulges in very different directions were near the surface rupture [45,54], which showed the earthquake was a complex event. As the most intuitive manifestation of an earthquake, the accurate measurement of the three-dimensional deformation contributes to a better understanding of crustal dynamics and seismic triggering mechanisms.
The GNSS stations from different organizations provided the research conditions for the discrete three-dimensional deformation in the Qinghai-Tibet Plateau. However, the regional co-seismic three-dimensional deformation could not be solved simply by interpolating GNSS data. As a complement for the low spatial resolution of GNSS stations, InSAR measurement could be used to obtain high spatial resolution one-dimensional deformation of the study area.

GNSS Data Estimated the Three-Dimensional Deformation
Using the GAMIT/GLOBK software, the time series of high accuracy GNSS threedimensional deformation could be obtained before and after the earthquake. This paper applied the continuous GNSS observation stations from different sources, including the Crustal Movement Observation Network of China, the provincial earthquake agency, enterprises, and more than 200 IGS stations uniformly distributed worldwide [55].
Firstly, all the stations used were loosely constrained, and their daily solutions by the GAMIT software were obtained, where the baseline solution used the RELAX mode, the satellite elevation cutoff angle was 10 • , the atmospheric zenith wet delay was solved every 2 h, and mapped to the satellite by the GMF mapping function. All files, including the receiver antenna file, the satellite antenna file, the receiver and satellite antenna phase center models, the solid tide model, the solar ephemeris table, the lunar ephemeris table, and the ocean tide model, were updated. The daily solutions were combined with the calculated global distribution of 200 IGS stations. then, we chose some global, even, stable, and selfconsistent IGS stations as constraints, and performed a seven-parameter transformation to tie the above solution to the International Terrestrial Reference Frame 2014 (ITRF14) by the GLOBK software. We could obtain the accurate three-dimensional coordinates time series under the ITRF 2014.
To avoid the after-slip effect and increase the GNSS measurement reliability, the coseismic displacement of GNSS stations before and after the earthquake was chosen based on the principle that the ratio of displacement to error was greater than 1. Some data after the earthquake was removed, so we used selected data, including 3 days before the earthquake and 1~3 days after the earthquake, to calculate the co-seismic deformation. The station time series were averaged, respectively, and the differential before and after the earthquake, the errors in each direction could be obtained by the error propagation law.
In Figure 2, we can see that more displacements were near the epicenter. In Figure 2a, the horizontal co-seismic displacement was symmetrically distributed in four quadrants. Moreover, some stations' horizontal co-seismic displacement directions were parallel to the seismogenic fault, which indicated it was a left-lateral strike-slip earthquake. In Figure 2b, we could see an uplift to the west and subsidence to the east on the northern side of the seismogenic fault, and subsidence and small displacement could be found on the southern side of the fault, which indicated it was uneven deformation in the vertical direction.  Figure 2. The co-seismic displacement obtained by GNSS data. The yellow lines are the seco tectonic block boundary, (a) represents the co-seismic displacement in horizontal direction, th and the red arrows represent the different displacement ratio, and (b) represents the co-seism placement in the vertical direction. The blue represents the uplift displacement, and the red sents the subsidence displacement: the red circle represents the error ellipse, the purple star sents the epicenter of the Maduo MS7.4 earthquake, the red curve near the epicenter represen

InSAR Data Estimated the One-Dimensional Deformation
InSAR was the space-to-earth observation combining synthetic aperture radar ( technology with interferometric technology. The SAR satellites are in an almost polar around the north and south poles of the Earth, thus it was very insensitive to the n south (NS) deformation of InSAR LOS direction, while the azimuthal direction c make up for the shortcoming of the InSAR LOS direction. The DInSAR technique c acquire the LOS deformation, and the MAI and POT techniques could acquire the and azimuthal deformation, in which the azimuthal accuracy of the MAI was signific better than that of the POT in areas with higher coherence, while the POT could decoherence and obtain the fault trace, which all belong to the InSAR measurement The SAR satellite provided two different geometry SAR images by ascending descending orbits, respectively. The research shows [56,57] that it was of benefit t proving the accuracy of the three-dimensional deformation field by both ascending descending at the same time. That was, the measured deformation should include asc ing LOS, descending LOS, ascending azimuth, and descending azimuth.
The European Space Agency (ESA) launched two C-band SAR satellites, the Sen 1A satellite, and Sentinel-1B satellite, in April 2014 and April 2016, respectively. The r period was 12 days for a single satellite, while the revisit period could be shortened days for two satellites. The Sentinel-1 observation data could be downloaded freely ESA official website (https://scihub.copernicus.eu/dhus/#/home, accessed on 12 July by registered users, and the corresponding precise orbit file can be downloaded f about 21 days after measurement on the other website (https://scihub.cop cus.eu/gnss/#/home, accessed on 12 July 2021). We downloaded the Sentinel-1A/B ferometric wide (IW) swath mode single look complex (SLC) data before and afte 2021 Maduo MS7.4 earthquake, as well as the precise orbit files at the corresponding Figure 2. The co-seismic displacement obtained by GNSS data. The yellow lines are the secondary tectonic block boundary, (a) represents the co-seismic displacement in horizontal direction, the blue and the red arrows represent the different displacement ratio, and (b) represents the co-seismic displacement in the vertical direction. The blue represents the uplift displacement, and the red represents the subsidence displacement: the red circle represents the error ellipse, the purple star represents the epicenter of the Maduo M S 7.4 earthquake, the red curve near the epicenter represents.

InSAR Data Estimated the One-Dimensional Deformation
InSAR was the space-to-earth observation combining synthetic aperture radar (SAR) technology with interferometric technology. The SAR satellites are in an almost polar orbit around the north and south poles of the Earth, thus it was very insensitive to the north-south (NS) deformation of InSAR LOS direction, while the azimuthal direction could make up for the shortcoming of the InSAR LOS direction. The DInSAR technique could acquire the LOS deformation, and the MAI and POT techniques could acquire the LOS and azimuthal deformation, in which the azimuthal accuracy of the MAI was significantly better than that of the POT in areas with higher coherence, while the POT could resist decoherence and obtain the fault trace, which all belong to the InSAR measurement.
The SAR satellite provided two different geometry SAR images by ascending and descending orbits, respectively. The research shows [56,57] that it was of benefit to improving the accuracy of the three-dimensional deformation field by both ascending and descending at the same time. That was, the measured deformation should include ascending LOS, descending LOS, ascending azimuth, and descending azimuth.
The European Space Agency (ESA) launched two C-band SAR satellites, the Sentinel-1A satellite, and Sentinel-1B satellite, in April 2014 and April 2016, respectively. The revisit period was 12 days for a single satellite, while the revisit period could be shortened to 6 days for two satellites. The Sentinel-1 observation data could be downloaded freely from ESA official website (https://scihub.copernicus.eu/dhus/#/home, accessed on 12 July 2021) by registered users, and the corresponding precise orbit file can be downloaded freely about 21 days after measurement on the other website (https://scihub.copernicus.eu/gnss/#/home, accessed on 12 July 2021). We downloaded the Sentinel-1A/B interferometric wide (IW) swath mode single look complex (SLC) data before and after the 2021 Maduo M S 7.4 earthquake, as well as the precise orbit files at the corresponding time. The specific information is shown in Table 1. In this paper, the InSAR measurements, including the DInSAR, MAI and POT techniques, were used to obtain interferograms for the ascending and descending Sentinel-1A/B satellites data, respectively. In the processing, the STRM DEM [58] with 30 m resolution was used to eliminate terrain error effect, the multi-looking was carried out 2:10 by azimuth resolution and range resolution to reduce the speckle noise of the SAR images and improve signal-to-noise ratio, the phase unwrapping was carried out by Goldstein filtering and the minimum cost flow solved the problem of the phase ambiguity, the precision orbit files could be used to re-flatten and eliminate orbital errors by selecting the ground control points and the polynomial refinement method, then the geocoding was carried out to obtain the LOS and azimuthal co-seismic deformation in the WGS84 coordinate system of the 2021 Maduo M S 7.4 earthquake. The results were expressed as DInSAR_LOS_des, DInSAR_LOS_as, MAI_LOS_des, MAI_LOS_as, MAI_AZI_des, MAI_AZI_as, POT_LOS_des, POT_ LOS_as, POT_AZI_des, and POT_AZI_as. For example, DInSAR_LOS_des represented the descending LOS deformation obtained by the DInSAR technique, and MAI_AZI_as represented the ascending azimuthal deformation obtained by the MAI technique. The specific co-seismic deformation results from different techniques are shown in Figure 3. 26 May 2021 (Sentinel-1B) IW 5. 6 6 In this paper, the InSAR measurements, including the DInSAR, MAI and POT techniques, were used to obtain interferograms for the ascending and descending Sentinel-1A/B satellites data, respectively. In the processing, the STRM DEM [58] with 30 m resolution was used to eliminate terrain error effect, the multi-looking was carried out 2:10 by azimuth resolution and range resolution to reduce the speckle noise of the SAR images and improve signal-to-noise ratio, the phase unwrapping was carried out by Goldstein filtering and the minimum cost flow solved the problem of the phase ambiguity, the precision orbit files could be used to re-flatten and eliminate orbital errors by selecting the ground control points and the polynomial refinement method, then the geocoding was carried out to obtain the LOS and azimuthal co-seismic deformation in the WGS84 coordinate system of the 2021 Maduo MS7.4 earthquake. The results were expressed as DIn-SAR_LOS_des, DInSAR_LOS_as, MAI_LOS_des, MAI_LOS_as, MAI_AZI_des, MAI_AZI_as, POT_LOS_des, POT_ LOS_as, POT_AZI_des, and POT_AZI_as. For example, DInSAR_LOS_des represented the descending LOS deformation obtained by the DIn-SAR technique, and MAI_AZI_as represented the ascending azimuthal deformation obtained by the MAI technique. The specific co-seismic deformation results from different techniques are shown in Figure 3.  The azimuthal accuracy of the POT technique was 1/30~1/10 of the pixel azimuthal resolution used. However, the azimuthal resolution was slightly low in the Sentinel-1 satellite, so that a lot of noise occurred in the deformation result from the POT. In this paper, the POT only was used to delineate the fault traces and eliminated non-homogeneous points by its resisting decoherence advantage. For the InSAR azimuthal deformation, we only used the MAI azimuthal direction deformation, which was more accurate than the POT.

Multi-Source Integrating the Three-Dimensional Deformation
According to the respective advantages of the GNSS and InSAR, we built the highaccuracy and high-spatial-resolution regional three-dimensional deformation fields, and provided the foundation for detailed analysis of the dynamic process in the area.
Quantitative description in seismology is based on solid continuum mechanics and the infinitesimal and homogeneous unit, often known as elasticity mechanics or elasticity theory. The displacement between a point r and the reference point 0 r was 0 d r r Δ = − . Based on the elasticity theory, strain was the spatial gradient of displacement, and was a relative value. The relation between the displacement and the strain could be expressed as When a geodynamic process occurs, such as earthquakes or volcanoes, and deforms a portion of the earth's surface, the surface is considered to be a uniform strain [16]. An arbitrary point P locates the position ( ) comprises unknown components, and the surrounding N points' The azimuthal accuracy of the POT technique was 1/30~1/10 of the pixel azimuthal resolution used. However, the azimuthal resolution was slightly low in the Sentinel-1 satellite, so that a lot of noise occurred in the deformation result from the POT. In this paper, the POT only was used to delineate the fault traces and eliminated non-homogeneous points by its resisting decoherence advantage. For the InSAR azimuthal deformation, we only used the MAI azimuthal direction deformation, which was more accurate than the POT.

Multi-Source Integrating the Three-Dimensional Deformation
According to the respective advantages of the GNSS and InSAR, we built the highaccuracy and high-spatial-resolution regional three-dimensional deformation fields, and provided the foundation for detailed analysis of the dynamic process in the area.
Quantitative description in seismology is based on solid continuum mechanics and the infinitesimal and homogeneous unit, often known as elasticity mechanics or elasticity theory. The displacement between a point r and the reference point r 0 was ∆d = r − r 0 . Based on the elasticity theory, strain was the spatial gradient of displacement, and was a relative value. The relation between the displacement and the strain could be expressed as ε = r−r 0 r . When a geodynamic process occurs, such as earthquakes or volcanoes, and deforms a portion of the earth's surface, the surface is considered to be a uniform strain [16]. An arbitrary point P locates the position x 0 e , x 0 n , x 0 u , and its three-dimensional displacement d 0 e , d 0 n , d 0 u comprises unknown components, and the surrounding N points' position x i e , x i n , x i u and the three-dimensional displacement d i e , d i n , d i u are known (i = 1, 2, . . . , N).
For the infinitesimal deformation, we could perform the Taylor expansion for each of the three components without taking into account the second orders and the above terms. Thus, the equation for the linear relationship between strain and displacement is as follows.
In which the ∆x i = x i − x 0 represents the vector distance between the point i and point P in the corresponding east (E), north (N), and vertical (U) directions, which are needed to convert the coordinates of all points to the same geodetic coordinate system. The H represented the displacement gradient tensor, which could be expressed by the symmetric strain tensor ε and anti-symmetric rigid body rotation tensor ω [16].
By bringing the Equation (2) into Equation (1), it could be expressed as In which, the unknown matrix was expressed as follows.
The coefficient matrix was as follows.
For the GNSS measurement, the GNSS stations discrete co-seismic three-dimensional deformation can be directly obtained by displacement difference before and after an earthquake. It can be expressed as follows.
L i e , L i n , L i u were the GNSS stations deformations in east-west, north-south and vertical directions, respectively.
For the InSAR measurement, only one-dimensional deformation could be obtained in the azimuthal or LOS direction, and needed to be decomposed into three-dimensional direction by a certain formula. The relationship between three-dimensional and onedimensional displacement of the same point was as follows.
In which azides represent the InSAR LOS and azimuthal displacement of the ascending and descending by the different measurement techniques, respectively. α i as , α i des , θ i as , θ i des represented the pixel i heading angles and elevation angles of the ascending and descending orbits, respectively. Equation (8) was a unified expression of the orbit and direction, which could be obtained by different techniques, for example, the L i losas was the ascending LOS displacement which may be obtained by the DInSAR or the MAI. Then, the above corresponding relationship of the five deformation types (DInSAR_LOS_des, DInSAR_LOS_as, MAI_LOS_des, MAI_LOS_as, and MAI_AZI_as) involved in the solution can be brought in Equation (8).
According to Equations (6)- (8), the relationship between the observations and the calculated values was as follows.
The error equations could be listed as follows.
Establishing the observation equation involved two critical steps: the selection of the surrounding observations that were related to the reference point, and the integration of the data by an appropriate weighting allocation. We needed to determine the deformation correlation distance of each reference point and adjust the relevant window size around the reference point, in order to establish the observation equations with the correlation conditions. On this basis, determining the reasonable weights for various observations was essential to obtaining high accuracy deformation results in the process of integrating different measurement technologies.

Adaptive Selection of the Relevant Points
The research showed that [25] the number of observations involved in the solution had a great influence on the calculation accuracy and efficiency. The larger the window size, the more pixels involved in the solution, the less the difference the root-mean-square error (RMSE) between the calculated result and GNSS displacement in the three components, and the three-dimensional deformation was more accurate. However, the deformation correlation between the observation and the reference point became weakened or even irrelevant with the increase in the window size, which would not satisfy the requirement of the same deformation model. Since it had an impact on the deformation correlation and the computational efficiency, the window size should have a certain constraint threshold. To this end, we established a variogram [59,60] by the InSAR pixels to determine the correlation distance in the area of interest. The geostatistical approach used in references [59,60] demonstrated that it was the exponential relationship between the covariance and distance. The equation was as follows. where C(h) was the covariance function [56] between pixels with a distance h in the study area, C 0 was the variance of the study area, and s was the unknown correlation distance and proportional to h. According to the relationship between the pixel spatial resolution and the distance, the h and s could be directly expressed as the window size at the same time.
When more than 200 pixels were involved in the solution, the RMSE tended to be steady [25]. We linearized formula (11) and established the observation equations by the 15 × 15 pixels as the initial window size. The correlation distance threshold could be determined by the least squares fitting with the window sizes being gradually increased, for each technological observation. To guarantee that each technological observations could fulfill the correlation distance, the final window size threshold was the minimum of the correlation distance thresholds of the various measurement technologies. The window size was automatically expanded until that more than 200 pixels were involved in the solution.
The heterogeneous points near the fault must be considered in the process of adaptively adjusting the window sizes, due to the deformation characteristics being different on both sides of the earthquake fault. To search heterogeneous points around the reference point and remove them, we obtained the specific earthquake fault trace and the boundary of the heterogeneous points, due to the POT being insensitive to decoherence. It was different from the empirical approach that all GNSS stations were involved in the solution; we chose the GNSS stations within the correlation distance in order to ensure the deformation correlation between the used observations and the reference point. It meant that the number of the GNSS stations and the InSAR pixels involved in the solution were variable pixel by pixel.
After determining the observations involved in the solution by adaptively adjusting the window size, we could establish the error equation. We considered weighting reasonably for the different observations.

Integrating Data Based on the VCE Weighting
The observations were classified according to their measurement accuracy. The GNSS observations were weighted as two types of observations, due to the different measurement accuracy in the horizontal and vertical directions. InSAR observations were weighted according to the directions and the orbit from different measurement technologies. Each type of observations was considered to be independent.
Firstly, the prior variances reciprocal was as the initial weights. For the GNSS measurement, the prior variances were the error square in the corresponding directions. For InSAR measurement, the prior variances were calculated by using the pixels within the adaptive window size base on the ergodicity assumption; the formula is as follows.
where σ 2 0j represented the prior variance of the type j observations, j represented the type j observations, d i j represented the pixel i deformation in the type j observations, k j represented the total number of the type j observations within the window size, d j represented the average of the k j pixels.
After estimating the unit weight variance of each type of observation by the VCE method, the corresponding new weight was given by iterating the initial weight, that was as follows.
where P j was the initial weight value of type j of observation, P j was the corresponding new weight that would participate in the new adjustment calculation, c represented a certain value generally, such as σ 2 j . The new unit weight variance could be obtained by the VCE method, iterating until the unit weight variances of all types of observations were approximately equal as follows.
In the weighting iteration process, the difference between the different types of unit weight variance should be less than a convergence threshold ε 2 , which depended on the observation accuracy. Comprehensively considering the observation accuracy of the GNSS and InSAR measurement, for example, the GNSS measurement accuracy was millimeter or submillimeter, the DInSAR measurement accuracy was centimeter or millimeter, and the MAI measurement accuracy was decimeter or centimeter. We chose the unit weight variance convergence threshold of 0.0001 m 2 . That is, the end-of-iteration condition was satisfied σ 2 max − σ 2 min < 0.0001. σ 2 max denoted the maximum value of the unit weight variances for the various types of observations, and σ 2 min denoted the minimum value.

Results and Analysis
In this paper, the three-dimensional surface deformation field in the hundreds of kilometers along the seismic fault near the epicenter was obtained by effectively integrating the GNSS and InSAR measurements. Figure 4 showed the regional co-seismic threedimensional deformation field.
The strike of the earthquake fault was NWW and locally near EW. Since the angle between the strike of the earthquake fault and the east-west direction was small, the eastwest deformation was dominant in the earthquake fault deformation. We considered that the east-west component result could indicate the earthquake fault movement mode. As shown in Figure 4, when the displacement was decomposed into three components, the displacement was large in the east-west direction with a wide range (See Figure 4a,d), and obvious left-lateral shear stress happened on both sides of the fault (the displacement to the west in the north of the earthquake fault, the displacement to the east in the south of the earthquake fault), which indicated that it was a typical left-lateral strike-slip earthquake. The deformations at both ends of the main surface rupture far away from the epicenter were larger than that near the epicenter. In the vertical, the deformation on both sides of the rupture was smaller and anti-symmetrically alternating between subsidence and uplift. These deformation characteristics were consistent with the results of the field geological survey [44,45,53] and were in agreement with the focal mechanism solution [61].
The field scientific investigation of the 2021 Maduo M S 7.4 earthquake indicated [44] that no obvious main surface rupture had occurred near the epicenter (shown in Figure 4d-f), while there was secondary crack development in a wide range (mainly seismic fractures). Away from the epicenter, obvious main surface ruptures had been found along the eastern and western sides of the earthquake fault. Continuing eastward along the fault, no main surface rupture was found but some small secondary cracks (mainly seismic fractures), until the main surface rupture appeared at the northeast and southwest bifurcation phenomenon. However, at the western end of the main surface rupture, the rupture strike also shifted suddenly to a southern deflection.
We overlaid the main surface rupture of the geological mapping with the threedimensional deformation field. In Figure 4a-c, the black line was the fitted rupture from the filed investigation [44] and the InSAR deformation field in this study. In Figure 4d-f, the black line is the surface rupture from the fine interpretation from the unmanned aerial vehicle (UAV). The east-west directional deformation in Figure 4a was in good agreement with the rupture strike. The apparent displacement to the west on the north of the earthquake fault and to the east on the south, also indicated that the earthquake was a NWW left-lateral strike slip (the local was near the EW strike slip). The displacement was mainly manifested in the east-west direction; we have calculated the displacement was 1.6 m to the east and 2.3 m to the west. Given a certain point as a reference, the maximum displacement was 1.6 m plus 2.3 m, about 4 m, which was consistent with the calculated maximum displacement of 4 m according to the statistical relationship between the maximum slip D max and the earthquake magnitude M W : M W = 6.81 + 0.781 · log(D max ) [62]. The strike of the earthquake fault was NWW and locally near EW. Since the ang between the strike of the earthquake fault and the east-west direction was small, the eas (c,f) the displacement of the vertical direction. In Figure 4a-c, the black line is the fitted rupture from the filed investigation [44] and the InSAR deformation field in this study. In Figure 4d-f, the black line is the surface rupture from the fine interpretation from the unmanned aerial vehicle (UAV), the red star represents the epicenter of the Maduo M S 7.4 earthquake.
The deformation range in the south of the earthquake fault were significantly larger than that in the north of the earthquake fault, which may indicate that the primary power source causing this earthquake was on the south of the earthquake fault, which was consistent with the knowledge that the Qinghai-Tibet Plateau has been driven northward by the Indian plate. Figure 4b showed a clear and intermittent north-south (NS) tensional movement along the main surface rupture, which was consistent with the field investigation [44,53] and the USGS providing the left-lateral strike-slip with a few normal faults. The strike suddenly deflected along the main surface rupture to the west, where the displacement of the north-south direction reached the maximum; we explained that it was likely blocked by some geological body (obstacle body) so that the internal accumulated energy had to be released by the southward change in the rupture strike, and the material was suddenly squeezed southward on the deflecting area. This deflection phenomenon was similar to that of the 2001 Kunlun Mountain Pass West M S 8.1 earthquake [63]. According to the field investigation, the continuity and large surface rupture were found here, including a series of the alternate tensile fissure and squeeze drums, as well as en echelon arrays of shear rupture [48,53]. The above energy release means could also be reasonably explained in the vertical direction. Figure 4c showed that the deformation in the vertical component was alternate subsidence (the maximum reached 0.3 m) and uplift (the maximum reached 0.68 m), which was consistent with the results of the vertical displacement being about 1 m by the unmanned aerial vehicle (UAV) measurement [45]. Where the main surface rupture started from the epicenter to the east, obvious subsidence ( Figure 4c) and the small horizontal displacement (Figure 4a) formed some left-stepping tensional cracks; those were consistent with the phenomenon of the small pull-apart basin found in the field [48].
The surface ruptures caused by the Maduo M S 7.4 earthquake were from the epicenter to the east and to the west and had obvious broom-shaped bifurcations at both ends of the rupture. This phenomenon was a typical manifestation of energy release and conformed with the principles of seismology. However, it was a relatively rare deformation feature within the Bayan Har block.
An important finding was that the vertical component had special deformation characteristics at the west end of the main surface rupture, called the near east-west main surface rupture. The boundary between the uplift and subsidence was located on the north side of the near east-west main surface rupture (See Figure 5), while the main surface rupture passed through the uplift areas. This meant that this main surface rupture segment dominated the pure strike-slip movement, while the vertical deformation might have been caused by a pre-existing fault (see the red line in Figure 5). Moreover, the result was consistent with the position of the presumed rupture in the field from Pan et al. [53]. A few deformation traces, but no large-scale surface rupture was obtained from the surface survey [53,64]. Comparing the position with the actual geological map and the geological field survey results from Xiong et al. [65][66][67], a strike-slip fault extending to the east bank of Eling Lake was coincident with the boundary between the subsidence and uplift (See Figure 6). The aftershock sequence relocation also showed that a dense belt of aftershocks occurred along the northwest on the north side of the near east-west main surface rupture [68]. The result provided the evidence for the theoretical hypothesis proposed by Peterson [69] and Livio et al. [70] that the large earthquake could not only cause rupture on the seismogenic fault, but also might trigger pre-existing faults or new faults to produce a surface rupture or weak deformation in areas far away from seismogenic faults.
Combining the deformation result with the field survey, the above analysis results proved the reliability of the calculated three-dimensional deformation, and illustrated the importance of the three-dimensional deformation field for the fine interpretation of the surface rupture and the fault kinematics characteristics. east bank of Eling Lake was coincident with the boundary between the subsidence and uplift (See Figure 6). The aftershock sequence relocation also showed that a dense belt of aftershocks occurred along the northwest on the north side of the near east-west main surface rupture [68]. The result provided the evidence for the theoretical hypothesis proposed by Peterson [69] and Livio et al. [70] that the large earthquake could not only cause rupture on the seismogenic fault, but also might trigger pre-existing faults or new faults to produce a surface rupture or weak deformation in areas far away from seismogenic faults. Combining the deformation result with the field survey, the above analysis results proved the reliability of the calculated three-dimensional deformation, and illustrated the importance of the three-dimensional deformation field for the fine interpretation of the surface rupture and the fault kinematics characteristics.

The Integrating Mutil-Source Data
While the InSAR could not directly provide three-dimensional deformation, the GNSS three-dimensional deformation played an important role in solving the regional three-dimensional deformation field. Compared to the GNSS measurement discreteness, the InSAR measurement was important in spatial scale. The integration of GNSS and In-SAR complement each other to maximize their advantages.
The differences in root-mean-square errors (RMSE) between calculated three-dimensional deformation and the corresponding GNSS stations were 0.98 cm, 5.64 cm, and 1.37 cm in the east-west, north-south, and vertical directions. The accuracy was better than that of only InSAR measurement, whose RMSEs were 5.2 cm and 12.2 cm in the east-west and north-south directions, but none in the vertical direction [30]. The accuracy of the integrating GNSS and InSAR increased by 81.2% and 53.8% in the east-west direction and

The Integrating Mutil-Source Data
While the InSAR could not directly provide three-dimensional deformation, the GNSS three-dimensional deformation played an important role in solving the regional threedimensional deformation field. Compared to the GNSS measurement discreteness, the InSAR measurement was important in spatial scale. The integration of GNSS and InSAR complement each other to maximize their advantages. The differences in root-mean-square errors (RMSE) between calculated three-dimensional deformation and the corresponding GNSS stations were 0.98 cm, 5.64 cm, and 1.37 cm in the east-west, north-south, and vertical directions. The accuracy was better than that of only InSAR measurement, whose RMSEs were 5.2 cm and 12.2 cm in the east-west and north-south directions, but none in the vertical direction [30]. The accuracy of the integrating GNSS and InSAR increased by 81.2% and 53.8% in the east-west direction and north-south direction, respectively, according to the more constrained condition, just as increasing the number of GNSS stations, considering deformation correlation, removing heterogeneous points around the reference points and the variance component estimation weighting. Even if the pixel resolution and the number of SAR satellites in this paper were not as good as those of only InSAR measurement using a multi-satellite and multi-technology approach [30], this result reflected the advantage of integrating the rich high-precision GNSS measurement and the high-spatial resolution InSAR measurement. Since the nadir angle of the SAR satellite was between 20 • and 40 • , the InSAR LOS direction deformation was more sensitive in the vertical direction than in the horizontal direction, while the InSAR azimuthal direction deformation had no component in the vertical direction and dominated the north-south deformation. Since the GNSS measurement within the correlation distance provided sufficient precision support, the accuracy of the north-south deformation was considerably improved. The accuracy of the north-south direction was still improved, even if the azimuthal direction displacement used only one type of InSAR measurement (MAI_AZI_as) (see Figure 7). The displacement of the GNSS stations was a little larger than that of the integrated result, because GNSS stations used 1~3 days of data before and after the earthquake; however, InSAR observations were from the 21 May, before the earthquake, and the 26 May, after the earthquake in 2021, so the InSAR deformation decreased slightly. correlation distance provided sufficient precision support, the accuracy of the north-south deformation was considerably improved. The accuracy of the north-south direction was still improved, even if the azimuthal direction displacement used only one type of InSAR measurement (MAI_AZI_as) (see Figure 7). The displacement of the GNSS stations was a little larger than that of the integrated result, because GNSS stations used 1~3 days of data before and after the earthquake; however, InSAR observations were from the 21 May, before the earthquake, and the 26 May, after the earthquake in 2021, so the InSAR deformation decreased slightly. For different types of observations, the calculation precision of the variance component estimation (VCE) weighing was better than the weighted least squares (WLS), which used the observation prior variance as the weight with some uncertainty assumptions, which made it difficult to accurately reflect the observation weight. We took the reciprocal of the observation prior variances within the adaptive window size as the initial weights, then carried out several iterations based on the VCE to obtain the reasonable weights of Figure 7. Comparison of GNSS corresponding three-dimensional deformation and the integrated three-dimensional deformation by GNSS and InSAR (the x-axis represents the names of the GNSS station for verification, the y-axis represents the displacement of the deformation, the red line represents the integrated displacement by the GNSS and InSAR in three directions, and the blue line represents the displacement of the corresponding GNSS stations in three directions).
For different types of observations, the calculation precision of the variance component estimation (VCE) weighing was better than the weighted least squares (WLS), which used the observation prior variance as the weight with some uncertainty assumptions, which made it difficult to accurately reflect the observation weight. We took the reciprocal of the observation prior variances within the adaptive window size as the initial weights, then carried out several iterations based on the VCE to obtain the reasonable weights of the observations. The analysis result verified that it was feasible to obtain the high accuracy three-dimensional deformation field by the VCE weighting.

The Geodynamic Characteristic
The Bayan Har block located in the central and eastern regions of the Qinghai-Tibet Plateau is one of the strongest seismically active regions in China nowadays [71]. Its southern and northern boundaries were controlled by two large strike-slip fault zones, the Ganzi-Yushu-Xianshuihe Fault and the East Kunlun Fault, and its eastern boundary consisted of the middle-southern segment of the Longmenshan fault, the Minjiang thrust fault, and the Huya thrust fault. This has been the area of many destructive earthquakes of magnitude 7.0 and above in history [65]. Although a large number of destructive earthquakes had occurred along the boundary of the Bayan Har block, some earthquakes have also occurred within the block, such as the 1947 M7 3 4 earthquake and the 2021 Maduo M S 7.4 earthquake.
The three components interpretation for the co-seismic deformation was helpful to realize the kinematic processing of the seismogenesis structure, and provided the reference for understanding the geological structure and the geodynamics in this region. The three-dimensional deformations showed that the seismogenesis fault had a typical leftlateral strike-slip movement, and the two ends of the main surface rupture had obvious displacement and strike deflections. It was speculated that the fault was blocked by some hard material (obstacle body) during the strike-slip process, which led to the main surface rupture strike southward deflection at the west end and formed two branches at the east end. This was more conducive to energy release in some way. The energy dispersion led to surface distributed deformation [73].
The east-west deformation dominating the NWW strike-slip displacement, had a good indicator of the regional geodynamic environment. The active source of this earthquake might lie in the south of the earthquake fault, which was consistent with the fact that the Qinghai-Tibet Plateau is continuously pushed by the southern Indian plate to the Eurasian plate [74]. The northern boundary and eastern boundary of the Qinghai-Tibet Plateau being blocked by the hard Ordos block and the Sichuan Basin, thus resulted in its southeast escape movement [75]. Several field investigations also indicated that there was a widely distributed deformation range in the south of the earthquake fault [72,[76][77][78]. Under the pushing action of the Indian plate, in the south of the Kunlun Maintain Pass-Maduo-Gande fault and in the north of the Qinghai-Tibet Plateau, the material moved faster to SEE than that in the north of the fault, so that the long period of stress accumulation led to this typical strike-slip earthquake.
The vertical deformation in both sides of the main surface rupture (See Figure 4c,f) being in opposite directions, showed that there were some little pivotal movement characteristics in this earthquake.
Intuitively, the southward deflection at the west end and the bifurcation rupture at the east end of the main surface rupture showed that the 2021 Maduo M S 7.4 earthquake caused at least three different strike ruptures. As field investigation from Li et al. [45] and Pan et al. [54] showed, the surface rupture caused by this earthquake had obvious space segmentation. XIAO et al. [79] used the Monte Carlo method and the steepest descent method (SDM) to study the slip distribution of this earthquake and divided the fault into three segments, which also verified the viewpoint of this paper. The result was consistent with the focal mechanism and rupture process research from Deng et al. [80], who believed that the seismogenic fault plane caused by the strong earthquake was not a single plane Sensors 2023, 23, 3793 18 of 25 structure. The above results indicated that the underground structure of the earthquake was complex [81].
The continuous surface rupture was the most obvious in the south of the Eling Lake. The rupture's geometric distribution was consistent with the bifurcation phenomenon of the aftershock sequence relocation [68] (see Figure 8), which might be blocked by some hard geological body (obstacle body) in the rupture process from the epicenter to both ends. The pre-existing fault could explain the vertical deformation on the north side of the near east-west main surface rupture. In the geological field survey, Pan et al. [54] and Han et al. [82] found that there was no obvious strike-slip dislocation on the north side of the west end of the rupture, but there were several parallel cracks extending to the northwest and disappearing into the Eling Lake. In the research on the late Quaternary activity of the related faults in the north of the Bayan Har block, we also found that there were several NW-NWW strike-slip active faults near main surface rupture of this earthquake [67] (see Figure 6). ends. The pre-existing fault could explain the vertical deformation on the north side of the near east-west main surface rupture. In the geological field survey, Pan et al. [54] and Han et al. [82] found that there was no obvious strike-slip dislocation on the north side of the west end of the rupture, but there were several parallel cracks extending to the northwest and disappearing into the Eling Lake. In the research on the late Quaternary activity of the related faults in the north of the Bayan Har block, we also found that there were several NW-NWW strike-slip active faults near main surface rupture of this earthquake [67] (see Figure 6). The surface rupture and deformation characteristics of the Maduo earthquake were significantly different from those of several strike-slip earthquakes that had occurred in the Qinghai-Tibet Plateau in the past, and this indicates the complexity of the earthquake. For example, the main surface rupture of the 2001 Kunlun Mountain Pass West MS8.1 earthquake [63], located on the northern boundary of the Bayan Har block, was somewhat similar to that of the Maduo earthquake; it also showed strike deflection at the west end of the rupture. However, the rupture process of the 2010 Yushu MS7.1 earthquake [83], which occurred on the southern boundary of the Bayan Har block, was a typical unilateral rupture, and the deformation characteristic was significantly different from that of the Maduo earthquake in the vertical direction, which was uplift and subsidence, but with no obvious alternation on both sides along the surface rupture of the Yushu earthquake.
In the surface rupture, we could see some left-stepping tensional cracks and pullapart basins; the mechanism of these phenomena was caused by the Riedel shear. The major strike-slip faults, which may be thousands of kilometers long and tens of kilometers wide, are in domains of simple shear, and the displacement can be monitored in hundreds of kilometers along a fault. Within the domain of simple shear, the most recently active The surface rupture and deformation characteristics of the Maduo earthquake were significantly different from those of several strike-slip earthquakes that had occurred in the Qinghai-Tibet Plateau in the past, and this indicates the complexity of the earthquake. For example, the main surface rupture of the 2001 Kunlun Mountain Pass West M S 8.1 earthquake [63], located on the northern boundary of the Bayan Har block, was somewhat similar to that of the Maduo earthquake; it also showed strike deflection at the west end of the rupture. However, the rupture process of the 2010 Yushu M S 7.1 earthquake [83], which occurred on the southern boundary of the Bayan Har block, was a typical unilateral rupture, and the deformation characteristic was significantly different from that of the Maduo earthquake in the vertical direction, which was uplift and subsidence, but with no obvious alternation on both sides along the surface rupture of the Yushu earthquake.
In the surface rupture, we could see some left-stepping tensional cracks and pull-apart basins; the mechanism of these phenomena was caused by the Riedel shear. The major strike-slip faults, which may be thousands of kilometers long and tens of kilometers wide, are in domains of simple shear, and the displacement can be monitored in hundreds of kilometers along a fault. Within the domain of simple shear, the most recently active strand might be a zone of active faulting that was only a few meters wide. The simple shear had a donoclinic symmetry of strain because it was rotational, and a greater variety of structure forms in simple shear than in pure shear [84]. The en echelon arrays composed of earthquake faults, extension fractures, and pressure ridges are very common in a strikeslip fault zone. Left-stepping and right-stepping in a left-lateral strike-slip fault zone is generally derived from extensional tectonics and compressional tectonics, respectively (see Figure 9a,b). Geometrically, an idealized Riedel shear zone was composed of six principal elements (Figure 9c), which were X shear and Y shear faults, Riedel (R) and R conjugate shears, T tensional fractures, and P shears, which were all oriented at specific angles with respect to the general trend of the shear zone [85]. The co-seismic Riedel shear structures revealed a left-lateral strike-slip sense for the seismogenic fault of the 2021 Maduo M S 7.4 earthquake [86]. The focal mechanism data for the Maduo earthquake showed two nodal faults trending nearly EW and NS, which was in accordance with the field investigations [75][76][77][78][79][80][81][82][83][84][85][86][87]. The Riedel shear model primarily controlled the earthquake surface ruptures of a strike-slip earthquake and dominated the formation and evolution of strike-slip faults from a wider perspective. of structure forms in simple shear than in pure shear [84]. The en echelon arrays composed of earthquake faults, extension fractures, and pressure ridges are very common in a strikeslip fault zone. Left-stepping and right-stepping in a left-lateral strike-slip fault zone is generally derived from extensional tectonics and compressional tectonics, respectively (see Figure 9a,b). Geometrically, an idealized Riedel shear zone was composed of six principal elements (Figure 9c), which were X shear and Y shear faults, Riedel (R) and R′ conjugate shears, T tensional fractures, and P shears, which were all oriented at specific angles with respect to the general trend of the shear zone [85]. The co-seismic Riedel shear structures revealed a left-lateral strike-slip sense for the seismogenic fault of the 2021 Maduo MS7.4 earthquake [86]. The focal mechanism data for the Maduo earthquake showed two nodal faults trending nearly EW and NS, which was in accordance with the field investigations [75][76][77][78][79][80][81][82][83][84][85][86][87]. The Riedel shear model primarily controlled the earthquake surface ruptures of a strike-slip earthquake and dominated the formation and evolution of strike-slip faults from a wider perspective. The calculated three-dimensional deformation field was in agreement with several geological field surveys [44,45,53,[66][67][68]72,[75][76][77][80][81][82][88][89][90][91][92][93][94][95][96][97]. There was no obvious surface rupture near the epicenter, while there were some large-scale forked-broom-shaped co-seismic surface ruptures on both ends of the rupture. This verified the research result [84] that the surface deformation and the rupture characteristics with the multiple bifurcation, tail folding, and spalling at the ends of the rupture, were very different from those of the main strike-slip section of the strike-slip fault.
Combining the surface rupture field investigation with the regional three-dimensional deformation field by integrating the GNSS and InSAR, we considered that the 2021 Maduo MS7.4 earthquake was a typical left-lateral strike-slip earthquake; the revealed kinematic characteristics of the seismogenic fault showed that there was a relatively obvious dip-slip component at the ends of the rupture.
In the inversion of the three-dimensional co-seismic deformation field, we also made the important discovery that there was the obvious left-lateral strike-slip deformation but no vertical deformation at the west end of the main surface rupture, while the pre-existing The calculated three-dimensional deformation field was in agreement with several geological field surveys [44,45,53,[66][67][68]72,[75][76][77][80][81][82][88][89][90][91][92][93][94][95][96][97]. There was no obvious surface rupture near the epicenter, while there were some large-scale forked-broom-shaped coseismic surface ruptures on both ends of the rupture. This verified the research result [84] that the surface deformation and the rupture characteristics with the multiple bifurcation, tail folding, and spalling at the ends of the rupture, were very different from those of the main strike-slip section of the strike-slip fault.
Combining the surface rupture field investigation with the regional three-dimensional deformation field by integrating the GNSS and InSAR, we considered that the 2021 Maduo M S 7.4 earthquake was a typical left-lateral strike-slip earthquake; the revealed kinematic characteristics of the seismogenic fault showed that there was a relatively obvious dip-slip component at the ends of the rupture.
In the inversion of the three-dimensional co-seismic deformation field, we also made the important discovery that there was the obvious left-lateral strike-slip deformation but no vertical deformation at the west end of the main surface rupture, while the pre-existing fault dominated the vertical deformation. This discovery was an important supplement to the co-seismic deformation characteristics obtained from the field investigation. It provided an important reference for further exploring this earthquake's deformation mechanism and regional geodynamics.

Conclusions
The deformation correlation between the reference point and the points involved in the solution was fully considered, and it was taken as the decorrelation distance range for selecting the points used. An adaptive method was proposed to obtain a high accuracy regional three-dimensional deformation field. Taking the 2021 Maduo MS7.4 earthquake as an example, based on the elasticity theory, we established the accurate regional co-seismic three-dimensional deformation field by effectively and rationally integrating GNSS and InSAR, and gave the crustal dynamics and geological explanation of the three-dimensional co-seismic deformation field. The conclusions could be drawn as follows.
(1) The accuracy obtained by integrating GNSS and InSAR increased by 81.2% and 53.8% in the east-west direction and north-south direction, respectively. Additional constraints were helpful to obtain a high accuracy regional three-dimensional deformation field, such as increasing the GNSS stations, considering deformation correlation, removing heterogeneous points around the reference points, and the variance component estimation weighting. (2) The spatial position and the left-lateral strike-slip motion of the earthquake were consistent with those of the Jiangcuo section of the NWW strike of the Kunlun Mountain Pass-Jiangcuo fault in the south of the East Kunlun Fault. The mechanism of the left-lateral left-stepping was caused by the Riedel shear, which dominated the formation and evolution of strike-slip faults. (3) The maximum slip of the earthquake was about 4 m, which was difficult to investigate in the geological field survey due to the lack of reliable markers, and the slip was consistent with the empirical statistical results. (4) The ruptures started from the epicenter to both ends and formed bifurcations, which were rarely found in the same type of previous earthquakes that occurred in the Qinghai-Tibet Plateau. (5) The pre-existing fault controlled the vertical deformation at the north side of the near east-west surface rupture, which only dominated the east-west strike-slip. The three-dimensional deformation was consistent with the surface rupture trace. This provided direct evidence for the theoretical hypothesis that large earthquakes can not only produce a surface rupture on seismogenic faults, but also trigger pre-existing faults or new faults to produce surface ruptures or weak deformation in areas far from seismogenic faults.
We recognized that the new seismic surface rupture did not necessarily completely control the distribution of seismic deformation. For example, in the western end of the surface rupture of the Maduo earthquake, the strike-slip deformation was distributed along the new near east-west surface rupture, while the vertical deformation was controlled by the pre-existing structure in the northwest direction.
The deformation distribution range in the south of the fault was wider than that in the north, which might indicate that the physical properties of the geological bodies were different on both sides of the fault. The northward extrusion of the Indian plate led to the formation of large strike-slip faults in the Qinghai-Tibet Plateau, and these faults constituted the boundary of the active block. The epicenter was located on the south of the East Kunlun fault, which was the northern boundary of the Bayan Har block. The co-seismic deformation field indicated that the thrusting effect from the southern Indian plate was transmitted to the epicenter. The block in the south of the earthquake fault moved faster in the SEE direction than the block in the north. In future, we would obtain the displacement distribution on both sides of the surface rupture through a series of deformation profiles, and compare it with the field investigation to further discuss the changing trend of deformation and its mechanism. In addition, the relationship between the earthquake fault and the East Kunlun fault, and the relationship between the new rupture and the pre-existing fault are topics to be further studied.
This series of results formed an important supplement to the field survey and provided a reference for the risk assessment of geological hazards [98] in the area of interest. It provided a novel idea for the integration of multi-source measurement technologies to improve seismic deformation monitoring. However, it should be noted that it will be a modeling challenge for the real earth.

Data Availability Statement:
The other data used to support the findings of this study are included within the article.