We calculated the multi-scale strain parameters of the Luding earthquake and its surroundings using GPS velocity field during 2017–2020 (272 stations, Figure 1), and quantitatively analyzed the evolution. The GPS velocity field was processed using the GAMIT/GLOBK software (Herring et al., 2006), and the computation process refers to Wang (2009). As shown by the green arrow in Figure 1, the GPS velocity field on the International Terrestrial Reference Frame (ITRF 2005) shows remarkable southeastern movement trend in the study area. In order to obtain the differential crustal movement characteristics between different blocks and faults, the regional Euler vector method with three parameters is used to simulate the crustal movement trend term, then the residual velocity field is used to solve strain parameters (Niu et al., 2002; Gao et al., 2018). Due to the small internal deformation of the South China block (SCB), which located in the east of the Sichuan-Yunnan block, it seems as a rigid block to fit the Euler rotation vector. After deducting overall rotations, the velocity field relative to the South China block is obtained (as shown by the red arrow in Figure 1). The result shows that the sinistral-slip movement along the eastern boundary of the Sichuan-Yunnan block is significant, while the movement of the LFZ along the southeastern boundary of Bayan Har is blocked and characterized by extrusion and dextral-slip movement.

The continuous GPS data used in the study are based on station time series that was based on the ITRF 2014, and were processed with the GAMIT, GLOBK and QOCA software packages (Herring et al., 2006). GPS single time series are subject to the influences of reference frame, ionosphere, troposphere, etc., and these influences are correlated to a certain extent. Therefore, it is necessary to find some parameters with physical significance that are less affected by common errors to describe crustal deformation information (Dong et al., 2002). The time series results for the baseline between two GPS stations for any given station pair, and the strain parameters calculated by spherical wavelet-based method are used for the analysis of deformation characteristics in our study.

The secondary blocks within the plate boundaries undergo rigid rotation and elastoplastic deformation under the action of the surrounding blocks. Therefore, after deducting the regional global rotation component fitted by the Euler motion model and velocity field, the spatial distribution of the residual velocity component can highlight the local deformation characteristics, thus reflecting the differentiated movements of different tectonic units. The calculation method of the Euler pole of the block is shown in following formula: V E V N = r − sin ⁡ ϕ cos ⁡ λ − sin ⁡ ϕ sin ⁡ λ cos ⁡ ϕ sin ⁡ λ − cos ⁡ λ 0 ω x ω y ω z

Where (V _E, V _N) is the velocity in eastern and northern direction, (λ, ϕ) is the longitude and latitude of the points inside the block, r is the mean curvature radius, (ω _x , ω _y , ω _z ) is the Euler vector of the block motion.

The methods for calculating the strain parameters according to GPS velocity field can be roughly divided into two types: physical and mathematical methods. The physical method is mainly based on dislocation theory or block kinematic model, and the strain parameter inversion is carried out with the velocity field as the constraint (Zhang et al., 2005; Zhao et al., 2013). The mathematical method is to solve the strain parameters directly from the velocity field according to the displacement-velocity differential relation. The least-squares collocation method (Jiang and Liu, 2010; Wu et al., 2011), spherical wavelet-based method (Tape et al., 2009) and Gussian weighting interpolation method (Shen et al., 2015) are commonly used. Different from deformation of plate edges, the seismogenic and dynamic environment in Chinese mainland is more complex under the combined action of the collision, subduction and extrusion around multiple plates, thus the tectonic deformation presents multi-scale spatial characteristics (Xu and Wu, 2016). In addition, related studies have shown that seismic deformation and strain accumulation within plates are usually concentrated in a narrow range of several kilometers to tens of meters near the fault zone (Du et al., 2010). Due to the localized multi-resolution analysis ability in spatial and frequency domain, spherical wavelet-based method can be applied to calculate the strain parameters, which can highlight the irregularly spaced characteristics of GPS stations and obtain the information of crustal movement and deformation at different spatial scales (Su et al., 2016; Jin et al., 2021).

The basic principle of spherical wavelet-based multiscale approach is to make use of the attenuation characteristics difference of wavelet function at different scales (the larger the scale factor q is, the faster the attenuation will be, and the smaller the range of influence will be), then match the smallest resolved process to the local spatial density of observations. The large scale mainly represents the main information of the velocity field, and the small scale describes the details. The above method reflects the process of constant approximation to the velocity field, so as to realize the calculation of strain parameters with multi-resolution (Cheng et al., 2015). The principle of determining the decomposition scale is as follows: Where stations are dense, wavelets with large scale are available, where stations are sparse, wavelets only with small scale are available. If the equilateral triangle is selected as a graphic unit to represent the whole sphere, the condition for establishing basis function is that there are at least three GPS observation points in the area covered by the graphic unit. The average side lengths of triangles at different scales are shown in Table 1 (Wang and Dahlen, 1995). The scale factor corresponding to the smallest graphic unit containing only three GPS points is the maximum decomposition scale q _max. While the minimum decomposition scale q _min is determined by the distribution range of observation points in study area. The specific approach for calculating q _min is as follows: Firstly, the spherical area S of the region surrounded by four points determined by the maximum and minimum latitude and longitude of all the GPS points is calculated, then the radius L ( L = S / π ) is solved, and the available decomposition scale q _min is determined, so that the spatial support of the spherical wavelet corresponding to the scale is less than 2L (Tape et al., 2009; Cheng et al., 2015). The corresponding relationship between spatial support and minimum decomposition scale is detailed in the literature (Tape et al., 2009).

GPS baselines are usually expressed as geodetic lines between two GPS points, which is independent of the reference frame compared to single point timeseries (Dong et al., 2002; Zhan et al., 2011; Liu et al., 2015). Because it can indicate the extension or contraction of the crust between two station sites, the tension, compression, and strike-slip activity of fault zones are often quantitatively studied by analyzing the baseline timeseries orthogonal or oblique to faults (Guo et al., 2009; Wu et al., 2018; Ma et al., 2022). The baseline distance increases when the slope is positive, whereas it decreases when the slope is negative. In this study, continuous observational data collected by CMONOC GPS stations since 2014 is adopted to analyze the various characteristics in the motion rates of nine GPS baselines across the XFZ.

In order to reduce the influence of baseline length on the analysis of deformation characteristics, the annual rate of change per unit length (ARCPL) of each baseline over the years are calculated. Figure 2A shows the baselines SCGZ-SCLT, SCGZ-SCJL, SCLH-SCJL express tensile changes, with ARCPL values of around 1.0 × 10⁻⁸/a, the baselines SCGZ-SCXJ and SCLH-SCXJ manifest compressive changes, with the ARCPL values of −4.0∼−4.2 × 10⁻⁸/a, indicating that the Ganzi-Luhuo segment of the northwestern XZF was dominated by sinistral-slip movement. The baseline SCLT-SCXJ nearly orthogonal to the XZF shows tensile variation, which indicates that this section has certain tensile strain characteristics. The baseline SCJL-SCXJ oblique crossing the southern section of the XFZ and passes through the northernmost aftershock sequence of the Luding earthquake, which shows tensile changes, with the ARCPL value of about 4.6 × 10⁻⁸/a and significant sinistral shear strain. The baseline QHBM-SCXJ inside the Bayan Har block shows compressive change, indicating the existence of extrusion strain in the block. The baseline SCJL-SCTQ expresses nearly orthogonal to the southern section of the XFZ, and passes through the Luding epicenter manifest tensile variation, with the ARCPL value of approximately 2.0 × 10⁻⁸/a.

As is shown in Figure 2B, between the Maduo M _S7.4 earthquake on 22 May 2021 and the Luding M _S6.8 earthquake on 5 September 2022, the baselines SCGZ-SCXJ, SCLH-SCXJ within the Bayan Har block accelerated compression, and the tension rate of the baselines SCGZ-SCLT, SCGZ-SCJL, SCLH-SCJL within the Sichuan-Yunnan block slowed down, indicating that the northern boundary faults of the Bayan Har block were affected by the coseismic rupture of Maduo earthquake, then the block accelerated eastward movement after decoupling. As a result, the strain accumulation in the divergence area bounded by the southern section of XFZ and the southwestern section of the LFZ is accelerated, which may be related to the preparation of the Lushan M _S6.1, Maerkang M _S6.0 and Luding M _S6.8 earthquakes since 2022.

Based on the GPS velocity fields from 1991 to 2016 (Wang and Shen, 2020) and 2017 to 2020, the spherical wavelet-based multiscale approach is applied to analyze the tectonic deformation at different scales and its spatio-temporal evolution characteristics around the Luding earthquake. As is shown in Figure 1, the green arrows represent the velocity field during 2017–2020 on the ITRF 2005, from which we can see that the southeastern movement trend in the whole study area is quite significant. In this study, the South China block is assumed to be rigid and used to fit the Euler rotation vector. The residual velocity field after subtracting the Euler rotation is the motion vector relative to the South China block, which is used as the input data for spherical wavelet-based approach. The above result shows that the sinistral-slip along the eastern boundary of the Sichuan-Yunnan block is remarkable, while the movement of the LFZ along the southeastern boundary of the Bayan Har block is blocked, indicating the characteristics of compression and dextral-slip.

According to the location of the GPS observation stations, the spatial distribution of the maximum decomposition scale in the study area is obtained by using the spatial scanning method. Figure 1 shows that the maximum decomposition scale of the whole study area is above 5, that is, the whole station spacing is less than 220.4 km. Only a few regions with the scale q _max = 9 and the station spacing less than 13.8 km. The regions with the scale q _max = 8 and the station spacing less than 27.5 km is concentrated in the eastern boundary of Sichuan-Yunnan block and the southwestern section of the LFZ. The scale q _max near the epicenter of the Luding earthquake is 7. Therefore, the maximum decomposition scale of the study area is determined to be 6∼8. 2L = 727 km is obtained from the spatial distribution range of the stations, and the minimum scale can be determined to be 3 refers to Table 1.

When the decomposition scale q is 3∼6, the overall deformation characteristics of the study area can be obtained. The positive value of the area strain rate represents tension, and the negative value represents compression. The calculation results in Figure 3A show that compression regions concentrate around the LFZ, and have a high consistency with the spatial distribution of thrust earthquakes in the region (represented by red beachballs). While extension regions concentrate in the interior of the Sichuan-Yunnan block and the eastern boundary of the Qiangtang block, which have a high consistency with the spatial distribution of the normal-fault earthquakes (represented by green beachballs). The maximum shear strain in Figure 3B reveals that the high value areas are distributed in the vicinity of the block boundary, especially the eastern boundary of the Sichuan-Yunnan block and the middle to southwest segment of the LFZ, which have a high consistency with the spatial distribution of strike-slip earthquakes (represented by blue beachballs). The large-scale strain characteristics obtained by the spherical wavelet-based method are in good agreement with the regional tectonic stress field reflected by the focal mechanism solution of historical earthquakes, which illustrates that the method is feasible to solve the crustal strain parameters. Although the results calculated by different methods vary in details, the overall strain distribution features between our and previous research are consistent (Pan and Shen, 2017; Li et al., 2019; Wang and Shen, 2020; Zhang et al., 2022).

The maximum decomposition scales calculated by the GPS velocity fields during 1991–2016 and 2017–2020 reveal that q _max near the Luding epicenter approximate to 7. Therefore, the decomposition scale is set to 3∼7 to calculate the crustal strain parameters of the epicenter and its vicinity, then the spatio-temporal evolution characteristics of regional strain are analyzed and studied. Compared with the results in Figure 3, the strain field with increased decomposition scale has higher spatial resolution and can reflect the local strain accumulation characteristics in a small area. Figures 4A, B show that the epicenter is in the transition zone between compression and extension regions, and the principal strain rate of the fault surrounding the earthquake is manifested as NNE-trending tension and NWW-trending compression. The area strain rate of different sections on the XFZ shows segmentation characteristics, while the LFZ has a high consistency, showing strong compressive strain. The evolution process of the maximum shear strain rate shown in Figures 4C, D indicate that the epicenter is located at the edge of high value region of the maximum shear strain rate during 1991–2016 and 2017–2020, and the shear deformation is stronger in the northern section. The local weakening of shear strain at the block boundary fault reveals that the fault movement near the epicenter is blocked and there is the possibility of fault locking and strain accumulation, which is basically consistent with the understanding that “strike-slip earthquakes mostly occur at the intersection zone between the edge of high maximum shear strain rate and the large strike-slip fault zone” (Wu et al., 2020).

Crustal deformation under the influence of different spatial ranges can be reflected in the strain field of corresponding scale, and the small-scale strain field has the advantage of detecting the local crustal deformation (Xu and Wu, 2016). Accordingly, the spatial distribution of the NE strain rate (nearly perpendicular to the XFZ) and the maximum shear strain rate with a decomposition scale q = 7 is calculated separately. The spatial range corresponding to this scale is about 55 km, which is roughly equivalent to the co-seismic rupture scale of 51 km of the Luring M _S6.8 earthquake obtained from source process inversion (http://www.itpcas.cas.cn/new_kycg/new_kyjz/202209/t20220906_6509485.html). Figure 5A indicates that the strain characteristics perpendicular to the fault direction near the Luding epicenter shows weak extension in recent years, which is consistent with the extension changes revealed in the GPS baseline SCJL-SCTQ across this region. Figure 5B show that the epicenter is still in the edge of the high shear strain region, and the local shear weakening phenomenon still exists. It is concluded that weakening zones at the edge of high shear strain region around large strike-slip faults may suffer the phenomenon of motion blocked and accelerated accumulation of strain energy, while the strain perpendicular to the fault shows the extension change, indicating that the normal stress of the fault trends to decrease. The fault is prone to rupture under the combined action of the above factors, and thus has a high seismic risk.

Luding M _S6.8 earthquake occurred in the southern section of the XFZ at the boundary of Bayan Har block and Sichuan-Yunnan block. The analysis of the tectonic motion evolution characteristics of the two blocks is of great significance for deepening the understanding of the relationship between crustal deformation and earthquakes, and the regional adjustment movement process after major earthquakes. For this reason, GPS velocity fields of the study region during 1999–2007, 2009–2013, 2013–2018 were collected and uniformly converted to the velocity fields with respect to the South China block (As is shown in Figure 6A). In order to quantitatively research the differential characteristics of tectonic motion evolution between blocks, the northern and eastern velocity components of GPS stations within the southeastern Bayan Har block (east of 93°E) and the Sichuan-Yunnan block were statistically analyzed. The purple represents the east-west velocity V _E, the blue represents the north-south velocity V _N, and the positive value of V _E represents the eastward motion, while the negative value represents the westward motion; The positive value of V _N represents northward motion, while the negative value represents southward motion. The “peaks” in each layer of Figures 6B, C denote the fitted distribution density of velocity components at various observation stations within the block, which can characterize the overall motion features of the block along the north-south or east-west directions.

The Sichuan Yunnan block located on the southeast edge of the lateral extrusion deformation of the Tibetan Plateau, and is the main strong earthquake activity area in the Sichuan Yunnan region (Huang et al., 2007). The overall movement of the block is mainly clockwise rotation. As is shown in the three periods across the past 20 years, the southward movement rate of the block is relatively stable. After the Wenchuan earthquake, the eastward movement rate of the Sichuan Yunnan block slightly slowed down, and the overall motion along the east and south directions did little change, indicating that it was less affected by the Wenchuan earthquake than the Bayan Har block. In combination with seismic activity, moderate and strong seismic activity is frequent inside the block, and the accumulated strain is released in the form of earthquake, and the overall tectonic movement evolution is not obvious.

The comparison between the results during 2009–2013 and 1999–2007 shows that after the 2008 Wenchuan M _S 8.0 earthquake, the eastward and southward movements of the Bayan Har block were significantly enhanced due to the decoupling influence of the LFZ. The distribution density of motion components during 2013–2018 reveals that the eastward and southward motion rates have weakened, which may be related to the rapid healing of the LFZ after the earthquake (Zhao et al., 2020). From the perspective of seismicity, after the 2008 Wenchuan earthquake, the moderately strong earthquakes within the Bayan Har block gradually migrated to the southeast with the accelerated eastward and southward motion (see Figure 7A). Figure 7B manifests that before and after Wenchuan earthquake, the cumulative Benioff strain rate in the Bayan Har block increases significantly, while the Sichuan-Yunnan block, which is less affected by the post-seismic effect, has little change. The relationship among post-seismic effect, block motion state and regional seismicity is further proved quantitatively. It shows that the decoupling of boundary fault zone caused by large earthquake will accelerate block motion, thus speeding up the stress accumulation of the existing locked fault, and strengthening regional seismicity.

When talking about the dynamic causes of moderate and strong earthquakes, the kinematic slip evolution, the deep structure, and the cumulative Coulomb stress change of the seismogenic fault are usually considered. Based on the long-term GPS observations, Guo et al. (2023) investigated the rate of interseismic loading and surface creeping along the Moxi segment, and the elastic dislocation model suggests that the slip rate of the Moxi fault is ∼9.7 mm/yr, implying that the seismogenic fault underwent energy accumulation at a fast-loading rate. Besides, the interseismic shallow creeping is too small to make up for the preseismic shallow slip deficit. A detailed 3-D crustal S-wave velocity structure by Song et al. (2023) revealed that the Luding earthquake and its aftershocks are mainly distributed at depth of 1–15 km, and these earthquakes occurred at the boundary of high-value anomalies, indicating that the relatively high-value anomalies reflecting hard materials in the southern XFZ are easy to accumulate energy and generate earthquakes. Based on the elastic dislocation theory and the stratified viscoelastic model, Li et al. (2020) calculated the cumulative Coulomb stress change on the main boundary faults and in inner blocks in Sichuan-Yunnan region, by combining the influence of coseismic dislocations of the M ≥ 7.0 historical strong earthquakes since 1515, M ≥ 8.0 events in and around this region, and the postseismic viscoelastic relaxation effect of the lower crust and upper mantle. The results showed that the Coulomb stress change increases significantly in the southern XFZ, in which the cumulative Coulomb stress change exceeds 0.1 MPa, resulting in the increase of regional seismic hazard.

Furthermore, our research has shown that large-scale crustal movement and changes in the dynamic environment have a promoting effect on the occurrence of strong earthquakes in the context of strong strain accumulation (Jiang et al., 2009). The 2021 Maduo M _S7.4 was the largest earthquake that occurred in the Bayan Har block since 2008. The baseline timeseries across the faults in the southeast part of XFZ expressed that the southeast margin of Bayan Har block accelerated motion under the influence of post-seismic adjustment. Then accelerated the rate of strain accumulation in this area against the background of strong fault locking, and finally caused fault rupture and triggered the Luding earthquake.