Source Model and Simulated Strong Ground Motion of the 2021 Yangbi, China Shallow Earthquake Constrained by InSAR Observations

Abstract

On 21 May 2021, an Mw 6.1 earthquake, causing considerable seismic damage, occurred in Yangbi County, Yunnan Province of China. To better understand the surface deformation pattern, source characteristics, seismic effect on nearby faults, and strong ground motion, we processed the ascending and descending SAR images using the interferometric synthetic aperture radar (InSAR) technique to capture the radar line-of-sight (LOS) directional and 2.5-dimensional deformation. The source model was inverted from the LOS deformation observations. We further analyzed the Coulomb failure stress (CFS) transfer and peak ground acceleration (PGA) simulation based on the preferred source model. The results suggest that the 2021 Yangbi earthquake was dextral faulting with the maximum slip of 0.9 m on an unknown blind shallow fault, and the total geodetic moment was 1.4 × 10¹⁸ Nm (Mw 6.06). Comprehensive analysis of the CFS transfer and geological tectonics suggests that the Dian–Xibei pull-apart basin is still suffering high seismic hazards. The PGA result demonstrates that the seismic intensity of this event reached up to VIII. The entire process from InSAR deformation to source modeling and strong ground motion simulation suggests that the InSAR technique will play an important role in the assessment of earthquake disasters in the case of the shortening of the SAR imaging interval.

1. Introduction

On 21 May 2021, an Mw 6.1 strong earthquake occurred in Yangbi County, in southwest China, about 30 km west of the famous tourist city Dali (Figure 1). According to the China Earthquake Networks Center (CENC), the epicenter was located at [99.87°N, 25.67°E], and initiated at 21:48:34 local time (13:48:34 UTC) with a surface wave magnitude of 6.4, leading to strong ground motion. The mainshock was felt by most residents of Yunnan Province and caused at least three deaths, dozens of injuries, and severe damage to many buildings. From 19 May 2021 to the time the mainshock took place, the total number of earthquakes above Mw 4.0 reached six. This implies that the seismic activity in this area increased dramatically before the mainshock. Similarly, the post-seismic activity was also at a high level, where the number of aftershocks seven days after the mainshock was recorded by the Yunnan seismic network, and relocated by a double-difference location method up to about 2000 [1]. This seismic sequence is thought to be the pattern of foreshock–mainshock–aftershock due to the seismic characteristics [2].

Different institutions including the Institute of Geophysics, China Earthquake Administration (IGP-CEA, http://www.cea-igp.ac.cn/kydt/278248.html, accessed on 5 May 2021), the United States Geological Survey (USGS, https://earthquake.usgs.gov/earthquakes/eventpage/us7000e532/executive, accessed on 5 June 2021), and the Harvard Global CMT catalog (gCMT, https://www.globalcmt.org, accessed on 5 June 2021) released the focal mechanism solutions of this event. The results of IGP-CEA showed that this event is a strike-slip earthquake with a centroid depth of 5 km on a high-angle (~82°) fault striking of 138°. Their fault rupture process model, inverted using far-field body-wave observations, demonstrated that the seismic energy was released within the first 8 s on a length (12–15 km) fault, and the moment magnitude was 6.1. The W-phase moment tensor solution reported by USGS was similar to the results of the IGP-CEA. For the gCMT, they published a different focal mechanism, in which the fault strike direction was 315° and the depth was 12 km.

Additionally, geodetic observations, especially the Global Positioning System (GPS) and the Interferometric Synthetic Aperture Radar (InSAR) techniques play an important role in measuring the ground deformation during the seismic cycle including interseismic [3,4,5,6,7], coseismic [8,9,10,11,12], and postseismic [13,14,15] stages. In the obtainment of the coseismic deformation, geodetic measurements can provide more information of the near-field to constrain the seismic source model than the far-field seismic data [16,17]. GPS and InSAR can obtain high precision deformation results and have been widely used in capturing the coseismic ground deformation fields. However, sparse distribution of the near-field GPS networks limits the constrain of the seismic slip model. Zhang et al. [18] inverted the seismic source model by using three-dimensional deformation recorded by 35 GPS sites, and obtained a fault slip distribution with the depths of 3 km to 12 km. Almost all of their GPS sites were located to the north of the epicentral area, moreover, only four of the GPS sites with obvious deformation were in the near field of the mainshock. The fault geometry was simply determined by analyzing the foreshock sequence, which might raise the uncertainty of the source model [18]. Meanwhile, Wang et al. [19] determined another slip source model of this event using the ascending and descending Sentinel-1 InSAR data and four GPS data used in the previous result, and suggested that the fault dipped to the southwest by the trial-and-error method and the slip depth ranged from 2 km to 10 km. They calculated the Coulomb failure stress (CFS) changes on dextral receiver faults to estimate the seismic hazard. Since the two dextral faults of Weixi–Weishan (WX-WS) and Honghe controlling the Dian–Xibei pull-apart basin [20,21,22] are probably normal slip characteristics, the CFS loading analysis was inadequate. There is still a discrepancy in the slip ranging depth between these two previous results. As a reliable source model is always crucial for better understanding regional hazards and can supply important input information for the simulation of strong ground motion, further investigation is still required to clarify which one is reliable for the mainshock. Moreover, there is still not any simulated strong ground motion, which is the basics of the assessment of earthquake disasters, based on the source model constrained by InSAR observations thus far.

In this study, to better understand the characteristics of the source slip model, its implications for seismic hazard and disaster caused by the Yangbi earthquake, we first processed the Sentinel-1 ascending and descending SAR data using the InSAR technique to obtain the ground deformation field caused by the mainshock. Then, we determined the fault geometry and the slip distribution, respectively, using nonlinear and linear inversion methods based on the constraint of the coseismic deformation obtained by the previous step. Accounting for the weak constraint of few GPS observations [18], we did not incorporate them into the inversion. Moreover, referring to the already known geological structures, we analyzed the CFS changes on the nearby faults to evaluate the seismic hazards on them. Additionally, we simulated the peak ground acceleration (PGA) using the determined slip model to understand the possible strong ground motion. The results in this study will help us comprehend the meaning of InSAR coseismic deformation, which can be used for the source model determination, seismic hazard evaluation, and possible disaster assessment.

2. Tectonic Setting

The Yangbi earthquake is linked to the active tectonics setting of the Yunnan and Tibetan Plateau. The Indian plate started to subduct beneath the Eurasian plate about 50 Ma ago and the subducting rate is approximately 40 mm/yr nowadays (Figure 1a) [23,24]. The NE extrusion due to this continent–continent collision led to the changes in the crustal characteristics of the Tibetan Plateau including the uplift, N–S direction shortening, and thickening [25,26]. Meanwhile, crustal materials escape through the southeastern margin, which has been demonstrated by topographic and numerical simulation studies [23,27]. This complicated process formed many major deep faults across the whole Tibetan Plateau. According to the characteristics of the active tectonics of China since the Late Quaternary, Deng et al. [28] divided the whole land into eight primary fault blocks surrounded by the deep faults, and the area in this study is located in the southeast of the biggest Tibetan Plateau block. Furthermore, they segmented each primary block into several secondary blocks with different faults. Chuan–Dian and Dian–Xinan are two secondary blocks belonging to the Tibetan Plateau primary block (Figure 1a), the boundary faults of which are characterized with lateral strike-slip due to the motion pattern of the southeastern margin of the Tibetan Plateau [23].

The WX-WS fault and the Honghe fault separate the Chuan–Dian and Dian–Xinan blocks and constitute the common geological boundary zone of these two secondary blocks [28]. The WX-WS fault is about 280 km long, its strike direction is north–west, and the dip angle ranges from 60° to 80° [29,30]. It consists of two parts, the northern segment, and the southern segment. Chang et al. [30] suggested that the pattern of the WX–WS fault movement is dextral with an average strike-slip rate of 1.8 mm/a–2.4 mm/a and dip-slip rate of 0.30 mm/a–0.35 mm/a since the Late Pleistocene, based on the field geological survey. The right-lateral strike-slip Honghe fault is far longer than the WX-WS fault. Its northern section is parallel with the southern segment of the WX-WS fault and they have similar geological movement. Even though different studies still argue whether the WX-WS fault belongs to the Honghe fault, these two faults jointly control the seismicity of the common boundary zone. Historical earthquake analysis showed that the parallel zone of these two faults is active and its seismic activity is relatively high [29]. The distribution of historical earthquakes of Mw > 4.0 in the past 50 years recorded by USGS also shows that most strong earthquakes occurred in this boundary zone (Figure 1b).

3. Data and Methods

3.1. SAR Data and Processing

The Sentinel-1 mission of the European Space Agency comprises a constellation of A and B polar-orbiting satellites, these two satellites carry the same C-band synthetic aperture radar sensor to image the Earth’s surface all day long. The SAR data imaged on these two different satellites can be used to implement InSAR processing. This type of design enables shorter intervals and relatively weakens the impact of time decorrelation between imageries. The Sentinel-1 constellation imaged the epicentral area of the Yangbi earthquake before and after the mainshock. We chose four images to construct two interferometric pairs with the imaging mode of the Terrain Observations with Progressive Scans (TOPS) along the ascending and descending orbit directions. These two pairs are the latest and have the smallest revisit time to avoid time decorrelation [31]. One of the four images was imaged on the Sentinel-1B satellite and the others on the Sentinel-1A satellite. The detailed information of these four images is listed in

Table 1. Detailed information of SAR data.

Sensor	Acquisition Time (M-D-Y)	Orbit Direction	Path Number	Frame Number	Heading Angle (°)
Sentinel-1A	20 May 2021	Ascending	99	1265	−12.5
Sentinel-1B	26 May 2021
Sentinel-1A	10 May 2021	Descending	135	508	−167.5
Sentinel-1A	22 May 2021

and the coverage is exhibited in Figure 1b.

The interferometric processing is implemented with the conventional two-pass differential InSAR strategy by the open-source software GMTSAR [32,33], which has been widely used in extracting ground deformation information from SAR images [34,35]. Since a bigger multi-look parameter can mitigate noise effects, we selected a multi-look factor 16:4 in range and azimuth to resample the image pixel spacing size of 100 m. The radar topography mission digital elevation model (SRTM DEM) [36] with 90 m resolution was employed to eliminate the phase contribution of the topography. We also filtered the two interferograms with a modified Goldstein filter algorithm [37] and unwrapped the interferometric area where coherence was higher than 0.07 using the statistical-cost network-flow algorithm [38]. Figure 2a,c shows the initial results of the geocoded ascending (ASC) and descending (DSC) orbit interferograms, where each fringe corresponds to 6 cm ground displacement along the LOS.

The atmospheric phase and the orbital errors are two main factors that affect the reliability of the coseismic deformation [39]. In particular, they usually cannot be neglected in small-scale deformation interferograms [39,40,41]. Here, we used the atmospheric error correction data (Figure S1b,f) generated by the generic atmospheric correction online Service for InSAR (GACOS) [42,43] to remove the atmospheric errors and Figure S1c,g are the deformation interferograms after the atmospheric errors were removed. The possible long-wavelength orbital errors were detrended with a linear ramp method [44]. The refined interferograms of ASC and DSC are shown in Figure 2b,d, respectively, and the final coseismic LOS displacements in Figure 3a,b.

The GACOS data show that the atmospheric errors were obvious for ASC and DSC interferograms in this study area, and they appeared to be related to the topography. For the ASC interferogram, the atmospheric error ranged from −0.03 m to 0.01 m; simultaneously, the positive and negative maximum of the atmospheric errors in the DSC interferograms were 0.02 m and −0.01 m. Therefore, the weakening of the atmospheric errors was notable for the final ground displacement, which showed a maximum value of 0.084 m. The ASC and DSC coseismic LOS deformation fields revealed that the ground displacements occurred surrounding Yangbi County, however, their patterns appeared to be different. In the ASC deformation map, there are two deformation areas. One is located in the NWW direction of the Yangbi and uplifts to the satellite along the LOS. The other enclosing Yangbi County shows a decline, and its deformation is not continuous due to InSAR decorrelations. The DSC deformation map showed another different deformation pattern; the main deformation areas looked like two lobes and were distributed SW–NE, but their displacement directions were the opposite. The ground of the northeastern lobe moved toward the satellite with the maximum value of 0.081 m, whereas, the other lobe moved away from the satellite with the maximum value of 0.070 m. Therefore, according to the above, we can see that the ASC and DSC coseismic deformation fields may present different information of the seismogenic source.

3.2. 2.5-D Surface Deformation

A three-dimensional (3-D) deformation field is important to understand the source characteristics qualitatively [45], however, it is difficult to obtain the real 3-D deformation directly using InSAR observations. According to the imaging geometry of SAR data, the deformation field obtained by a single InSAR observation was the integration of the north-south (N–S), east–west (E–W), and vertical displacements along the LOS of the satellite. The projection of the 3-D displacements to the radar LOS displacement is related to the incidence angle and the heading direction of the satellite, and Equation (1) shows the relationship between them [46].

d_LOS = d_v cos(θ) − d_n sin(θ) cos(α-3π/2) − d_e sin(θ) sin(α-3π/2)

(1)

where θ and α correspond with the incidence angle and the heading angle of the satellite, respectively.

Wright et al. [47] proposed a strategy resolving three-dimensional deformation from multiple LOS interferograms whereas Wright et al. [47] suggested that the E–W and vertical deformation could be derived by assuming the north component is negligible. Wang et al. [48] also demonstrated that the precision of the E–W and vertical components was improved by neglecting the north contribution. In this study, we only obtained two LOS deformation fields and could not resolve the 3-D deformation. The incidence angle and heading angle of the ASC deformation were 34.1° and −12.4°, respectively, and 34.2° and −167.5° for the DSC deformation. Therefore, only −0.120 and −0.122 in the ASC and DSC interferograms of the north component contributed to the LOS deformation, which similarly showed that the InSAR observations were the least insensitive to the north deformation. Assuming no north deformation detected in our interferograms, we resolved the E–W and vertical deformation, also called 2.5-D surface deformation [49], based on Equation (2) using the ASC and DSC deformation.

d_LOS = d_v cos(θ) − d_e sin(θ) sin(α-3π/2)

(2)

The 2.5-D surface deformation caused by the mainshock is shown in Figure 4. The E–W deformation field shows that the relative motion was along a NW–SE boundary line. The maximum eastward displacement was 0.119 m and the westward maximum was 0.062 m. There exists a slight uplift or subsidence in the vertical deformation map. This deformation pattern implies that the mainshock is a dextral slip dominant event.

3.3. Source Modeling Method

The source model is the basis for calculating the coseismic CFS changes and simulating the peak ground acceleration. To better retrieve the source model, we applied a traditional two-step strategy including non-linear and linear inversions to determine the geometry and the slip distribution of the seismogenic fault embedded in a homogeneous elastic half-space [50]. As above-mentioned, the centers of surface deformation were different in the ASC and DSC interferograms, so we carried out the inversion three times separately based on only ASC, only DSC, and joint ASC and DSC observations to retrieve a reliable source model.

The advantage of the high spatial resolution of the InSAR technique results in an enormous quantity of surface deformation points. Since the deformation gradient is not notable in ASC and DSC interferograms, therefore, we downsampled the ASC and DSC interferograms with a uniform approach to obtain fewer ground points as constraining observations before the inversion implementation. The total constraining points used were 2870 and 2877 for ASC and DSC, respectively (Figure S2).

In the first step, the geometric parameters of the seismogenic fault include the fault position (horizontal coordinates), depth, strike angle, dip angle, length, width, strike-slip, and dip-slip. These geometric parameters are non-linear with the surface deformation. In the non-linear inversion, we assumed the fault ruptured with uniform slip and no tensile motion happened. Here, we used a source model inversion software called PSOKINV [51], which can implement the non-linear inversion and has been widely used in inverting fault parameters [51]. The integrated approach in PSOKINV is a hybrid of particle swarm optimization and simplex to complete a global search for the undetermined parameters. The equal weight ratio between the ASC and DSC was used in the inversion. To achieve more reliable geometric parameters, we first used geodetic Bayesian inversion software (GBIS) to estimate the ASC and DSC deformation errors, then simulated 200 datasets with the estimated error information [52,53]. We constrained the dextral slip mechanism in the non-linear inversion according to the 2.5-D surface deformation. After 200 times of non-linear inversion, we obtained the geometric parameters constrained by the joint ASC and DSC observations (Figure S3 and Table S1). The inversion results constrained only by ASC or DSC observations were similar to the parameters by joint observations, therefore, we implemented the non-linear inversion only once for each orbital observation.

We fixed the fault geometry and extended the fault length and width to 30 km and 10 km, respectively, to weaken the effects of the border [31]. Since the relocated aftershocks are often used to delineate the geometric characteristics of the seismogenic fault, we compared our inverted geometric parameters with the relocated aftershocks recorded within seven days after the mainshock [1]. Three new constructed fault planes are shown in Figure S4. The results of these three geometric fault planes are well consistent with the aftershock distribution along the strike and dip directions. This demonstrates that the fault geometric parameters are reliable. Fixing the spatial geometry makes the inversion parameters decrease, and only the strike and dip slips are left, which are linear with the surface deformation. To receive the slip distribution, we divided the fault plane into 300 fault patches with a slip resolution of 1 km × 1 km and set up relationships between the fault distributed unit slips and the surface deformation to construct the Green’s matrix G based on the rectangular dislocation model in a homogeneous elastic half-space [50]. By adding a smoothing constrain to avoid potential oscillations, the second step inversion is resolved with the following equation [8]:

$$\left[\begin{array}{c}{d}_{LOS}\\ 0\end{array}\right]=\left[\begin{array}{c}G\\ k·{\nabla }^2\end{array}\right]·s$$

(3)

where d_LOS is the vector of different deformation observations; G is the Green’s matrix; s is the distributed slip vector; k is a second-order finite difference operator; and ${\nabla }^2$ is the smoothing factor. According to the rapid field investigation after the Yangbi event (https://www.eq-igl.ac.cn/kydt/info/2021/33884.html, accessed on 11 June 2021), there were no apparent surface ruptures that occurred in the epicentral area. Therefore, assuming all the boundary patches did not slip during the mainshock, we employed a bounded constrained least-squares method to minimize the residuals. The final smoothing factor of 3.5 (Figure S5) was determined by the tradeoff between the slip distribution roughness and the residual.

3.4. Calculation of Coulomb Failure Stress Changes

Spatial distribution of the CFS changes yielded by the mainshock are usually correlated with the spatial pattern of aftershocks, where the majority of aftershocks are located at the positive CFS increment area [16] and a small number of aftershocks in the negative area. Furthermore, the CFS increment will promote strong earthquakes on nearby active faults, and the decrease in the CFS delays the possible earthquakes (the area called stress shadow) [54,55]. Therefore, the static CFS changes resulting from the mainshock were usually used to qualitatively tackle the possible distribution of sequent aftershocks, and strong earthquake promotion and delay on nearby faults [56]. Even though the triggering or retard threshold of 0.01 MPa (0.1 bar) [57] has not been tested statistically by adequate events, the value of this CFS change is often put into use for the earthquake promotion or delay evaluation [55].

Under the Coulomb failure assumption, the CFS change on the receiver fault caused by a neighboring fault slip event is defined as:

$$\Delta CFS=\Delta \tau +{\mu }^{\prime }\Delta {\sigma }_n$$

(4)

where $\Delta \tau$ and $\Delta {\sigma }_n$ represent the shear stress change and the normal stress change respectively; and ${\mu }^{\prime }$ is the effective friction coefficient of the receiver fault. In this study, we applied Coulomb 3.4, an open-source software, to investigate CFS changes in the surrounding area and nearby faults based on the preferred source model above. We selected 0.4 as the effective coefficient of friction in the calculation of CFS changes.

3.5. Simulation of Strong Ground Motion

Strong ground motion delineates the pattern and severity of shaking caused by strong earthquakes. It is critical information to understand the extent of the areas affected and potentially hardest hit areas. The maximum value of two horizontal components of peak ground acceleration (PGA) is one of the most important parameters of strong ground motion. To better understand the strong ground motion, we calculated the PGA based on the preferred source model. We adopted the ground motion predication equations (GMPEs) for bedrock of PGA in western China, which was integrated into the fourth-generation zoning map of China [58]. Taking into account the local site effect [59], we employed site-amplification effects determined by the relationship between the topographical gradient (from USGS) and Vs30 (average velocity of the shear wave within 30 m below the surface) [60,61], and transformed the ground motion parameters on the bedrock to the ones in the soil. The source model inverted in this study was used as the input source information including the fault rupture range, fault geometric parameters, and magnitude. Since zero-slip was set for the four boundaries of the fault in the inversion of the source model, therefore, the input rupture range did not include the four boundaries. According to the Chinese Seismic Intensity Scale, which provides suggestions about the corresponding relationship between the PGA variational range and the seismic intensity, we determined the coverage of the seismic intensity VIII.

4. Results

4.1. Source Model

Three source models (named model-A, model-D, model-J) by fitting to the three datasets including ASC, DSC, and joint observations are shown in Figure 5. All these models show that the Yangbi earthquake was dominated by dextral motions. They present similar slip patterns and the rupture concentrated at the fault plane center. Model-A, with a maximum slip of 1.00 m, was more similar to model-J than model-D, and the peak slip as 0.88 m for model-D and 0.94 m for model-J. The simulated deformation and the residuals of model-J are shown in Figure 6, and the results of the simulations and residuals for model-A and model-D are shown in Figures S6 and S7. From the simulated results of model-A and model-D, we found that model-A fitted the ASC observations well to some extent, however, the residuals fitting the DSC deformation were still notable and appeared as whole interferometric fringes. In contrast, model-D fit the DSC deformation well, and the residuals for the ASC observations were obvious. Figure 6 shows that the model-J not only fit the ASC deformation well, but also the DSC observations perfectly, and the root mean square was 1.0 cm.

Many studies [16,62,63] on the spatial relationship between the mainshock and aftershocks showed that aftershocks are usually complementary to the mainshock in space. Therefore, we also used the relocated aftershocks to analyze the three source models. Most of the aftershocks occurred below the depth of 5 km, enclosed the main slip area of the three models, and appeared to be complementary to the mainshock slip pattern. Furthermore, the complementarity between model-J slip distribution and the aftershocks was better than the other two models. Therefore, based on the analysis above, we preferred model-J as the source model of the Yangbi earthquake to implement further analysis. According to the optimal source model-J, the calculated geodetic moment tensor was 1.4 × 10¹⁸ Nm (Mw 6.06), which was slightly smaller than the IGP-CEA seismic results.

4.2. CFS Change

Assuming the receiver fault is the same as the source fault, the mechanism of the receiver fault includes a strike angle of 138.8°, a dip angle of 87.2°, and a slip rake of −174°. The calculated horizontal CFS change field at a depth 5 km is shown in Figure 7a, and the CFS change along the AA’ profile is shown in Figure 7b. Figure 7a shows that the CFS of the two elongation zones and the perpendicular zones of the fault increased, and the maximum positive CFS change was 12.5 bar. The profile CFS change along the AA’ line implies two CFS increase zones distributed in the shallow and deep depths, and one CFS decrease zone in the mid-depth. About half of the aftershocks were distributed in the CFS shadow. The possible reason is that the cumulative stress in the CFS shadow area was not totally released during the mainshock and was subsequently released by the aftershocks.

To predict the possible CFS changes caused by this event on the nearby faults, we constructed the location and directional strike parameters of surrounding faults according to the geological results [28] as the receiver faults. We divided all the faults into different segments with the fault inflection points. For expediency, the dip angle of each fault was set to 90°. Considering the crustal thickness of Chuandian [64], we set all the fault widths to 20 km. The CFS changes on each receiver fault related to dextral and normal mechanisms were calculated respectively. Figure 7c shows the predicted CFS changes on the nearby faults when the mechanism of the receiver fault is dextral. The CFS on most faults decreased, except for the middle segment of the WX-WS fault, on which the peak CFS increment reached up to 0.25 bar, and the northern part of the H–H fault. Figure 7d shows the CFS changes on the normal receiver faults, which imply that the CFS was loaded on the mid-part of the WX–WS fault and the position of the positive CFS change was slightly more southern than the position of the dextral mechanism. The maximum of the CFS increase was 0.5 bar, which is apparently bigger than the threshold 0.1 bar.

4.3. Peak Ground Acceleration (PGA)

Figure 8 shows the result of the simulated PGA. The pattern of the simulated PGA looks like a slim flat disc, which may have resulted from the almost vertical fault geometry. The Chinese Seismic Intensity Scale provides suggestions about the corresponding relationship between the PGA variational range and the seismic intensity. The PGA corresponding to the seismic intensity VIII ranges from 178 cm/s² to 353 cm/s², and the PGA of IX ranges from 354 cm/s² to 707 cm/s². It was found that most of the PGA larger than 178 cm/s² were concentrated along the seismogenic fault zone in Yangbi County and did not exceed the boundary of Yangbi County. The coverage of seismic intensity VIII reached up to 717 km². The maximum PGA of 488 cm/s² was located northwest of Yangbi County. Since the simulated PGAs of only several points were bigger than 353 cm/s², we propose that the biggest seismic intensity caused by the Yangbi event without accounting for these points was VIII.

5. Discussion

Here, we carried out a checkerboard test to verify our inversion reliability and the resolution of observations. In the checkerboard test, the geometric parameters and the subfault size of the constructed fault planar were the same as the previous model inversion. We assumed three asperities ruptured on the fault each with 6 × 4 subfaults and specified with a 1 m purely left-lateral slip (Figure 9a). Synthetic deformation for each surface point of the two downsampled InSAR datasets was calculated. During the checkerboard joint inversion for all the synthetic deformation, we applied the same inversion method and smoothing factor as the real data inversion. The result (Figure 9b) shows that the three asperities with a 1 m slip along the 2021 Yangbi fault can be well identified by the synthetic deformation. The sensitivity of the surface observations along the strike direction is better than the one along the dip direction, which demonstrates that the surface deformation observations were less sensitive to the slip changes along the fault dip direction. Even so, the data used in this study were strong enough to reveal the slip distribution of the 2021 Yangbi earthquake on the fixed fault plane. The slip pattern of our preferred source model was similar to the ones of Zhang et al. [18] and Wang et al. [19], however, the difference compared to their studies was mainly the slip ranging-depth. Zhang et al. [18] suggested that the GPS observations used in their inversion were spatially sparse and less sensitive to the resolution of the slip distribution; additionally, they simply determined the fault geometry. Wang et al. [19] determined the fault position and strike angle by integrating the InSAR deformation, relocated aftershocks, and the seismic focal mechanism. Moreover, the dip direction and the dip angle were fixed by the trial-and-error method. They finally combined the four GNSS observations near the epicenter and InSAR data to obtain the source model. In contrast, we inverted the fault geometry by a non-linear inversion only with InSAR observations. A different method used in the determination of the fault geometry is one possible reason for the discrepancy; different data constraint is the other possible reason. Since all the source models can fit the constraint observations well, it is hard to tell which is most likely close to the truth. The spatial complementary between the mainshock slip area and aftershock distribution is a rough method to deduce the reliability of the source model. Figure 5c shows that the slip area of our preferred model is well complementary to the relocated aftershock distribution. This may partly imply that the source model determined in this study could be reasonable.

The two dextral faults of WX–WS and Honghe enclose the Dian–Xibei Basin with other faults. The characteristics and the evolution suggest that the Dian–Xibei basin is a pull-apart basin [20,21,22]; furthermore, this type of pull-apart basin is always related to the high seismic activity [65]. Our preferred source model inverted from the InSAR observations shows that the 2021 Yangbi earthquake ruptured an unknown shallow blind fault running in a S–E direction. Since the seismogenic fault is away from all the known faults, therefore, we give the name of the Yangbi fault to this blind fault, which is parallel to the WX–WS fault and the perpendicular distance to the WX–WS fault was approximately 20 km. The Yangbi fault has the same motion mechanism as the WX–WS and Honghe faults, therefore, we can infer that the Yangbi, as a secondary fault together with WX–WS and Honghe faults, controlled the seismic activity of the pull-apart basin. The maximum shear strain rate derived from GPS horizontal velocity fields demonstrated that the Dali and Lijiang regions are the two areas with high values larger than 100 × 10⁻⁹ a⁻¹ [66]. However, geological studies found that the horizontal slip rate of the Honghe fault is 5 mm/a [67], similarly, the dextral horizontal slip rate of WX–WS fault ranges from 1.8 to 2.4 mm/a smaller than the Honghe fault [30]. Chang et al. [30] suggested that the earthquake recurrence period of the two faults of WX–WS and Honghe is relatively long. Therefore, we deduced that the accumulated strain was partly allocated to the Yangbi fault during the interseismic period and released on 21 May 2021. More work should be carried out on the nucleation of the 2021 Yangbi earthquake to understand strain accumulation and allocation on different faults during the interseismic period. The CFS results of dextral and normal mechanisms on the middle part of the WX–WS fault and the northern section of the Honghe fault appear to have increased since the Yangbi earthquake. Considering the location of these two fault sections located at the margin of the pull-apart basin, we should pay more attention to their seismic hazard in the future.

A near-field strong-motion network with adequate density is ideal to estimate the PGA field caused by a strong earthquake. In the area distributed with dense strong motion stations, the strong motion records are used to fit the ground motion attenuation representing the attenuation characteristics of the earthquake, then the PGA of the entire epicentral area can be obtained by interpolation [68], since this type of processing strategy to obtain the PGA has been applied in areas with dense strong motion stations [69]. Since strong motion observations are sparse or there are no strong motion records in many active earthquake regions, the source characteristics are an important input to achieve the PGA. The source models are always inverted with far-field or near-field seismic records, even jointly with near-field geodetic high-rate observations [70,71]. However, not all the near-field areas are set up with seismic or GPS stations before earthquakes, additionally, the far-field seismic data may not fully reflect the characteristics of the seismic source. We simulated the PGA caused by the Yangbi event with the source model input not integrating the strong motion station records. To validate the reliability of the simulated PGA results, we collected the strong motion data recorded by CENC up to 200 km distances from the Yangbi earthquake, selected the maximum of two horizontal components of the records as the PGA value, and compared the simulated PGA results with the observed ones. We plotted the observed PGA results and the simulated ones of the corresponding positions in Figure 10 as a function of the distance to the surface projection of the fault (Joyner–Boore distance, R_jb). Figure 10 shows that most of the observed PGAs are distributed among the range of two standard deviations (2σ) of the attenuation model, which demonstrates that the attenuation model is reliable in this study. Furthermore, the differences between the observed and simulated PGAs fell into the range of 2σ. All of these suggest that the simulated PGAs are reliable, to a certain extent. Nowadays, the InSAR observation has become vital data to constrain the source model. Cheloni et al. [72] inverted the source model of the 2021 Mw 6.8 Elazığ earthquake and simulated its strong ground motion. We also predicted the PGA caused by the 2021 Yangbi earthquake using the InSAR deformation. Since timeliness is crucially important for earthquake disaster assessments using the PGA field, the temporal resolution limits its application in seismic disaster assessments. With the shortening of SAR imaging interval, we can obtain timely surface deformation and the source model for predicting strong ground motion. It is predicted that in the future, InSAR technique will provide more timely and abundant basic data for earthquake disaster assessment.

6. Conclusions

The geometric parameters and slip distribution of the 2021 Yangbi earthquake inverted by ASC and DSC InSAR observations showed that the Yangbi earthquake occurred on an unknown blind shallow fault, the dextral faulting rupture mainly occurred below the surface within 10 km, and the total geodetic moment was 1.4 × 10¹⁸ Nm (Mw 6.06). Based on this source model, we carried out the analysis of CFS change and the simulation of surface PGA. Combining the geological studies of the Dian–Xibei pull-apart basin and CFS loading on the nearby faults, we infer that the Dian–Xibei pull-apart basin is still suffering a high seismic hazard. According to the PGA simulation, the maximum PGA value reached up to 488 cm/s², which suggested that the seismic intensity was at least VIII caused by the 2021 Yangbi event.