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

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 × 1018 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.


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-mainshockaftershock 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 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).

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 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 Remote Sens. 2021, 13, 4138 5 of 18 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.
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 99.60°100. 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.
Remote Sens. 2021, 13, x FOR PEER REVIEW 6 of 18 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.

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 northsouth (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]. dLOS = dv cos(θ) − dn sin(θ) cos(α-3π/2) − de sin(θ) sin(α-3π/2) 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.
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

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 northsouth (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].
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.
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

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 nonlinear 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].

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]: 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 ∇ 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.

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: where ∆τ and ∆σ n represent the shear stress change and the normal stress change respectively; and µ 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.

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.  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.

Source Model
Remote Sens. 2021, 13, x FOR PEER REVIEW 9 of 18 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.

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. 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 18 Nm (Mw 6.06), which was slightly smaller than the IGP-CEA seismic results. 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 18 Nm (Mw 6.06), which was slightly smaller than the IGP-CEA seismic results.

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

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. 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.  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 2 to 353 cm/s 2 , and the PGA of IX ranges from 354 cm/s 2 to 707 cm/s 2 . It was found that most of the PGA larger than 178 cm/s 2 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 2 . The maximum PGA of 488 cm/s 2 was located northwest of Yangbi County. Since the simulated PGAs of only several points were bigger than 353 cm/s 2 , we propose that the biggest seismic intensity caused by the Yangbi event without accounting for these points was VIII. 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. 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 2 to 353 cm/s 2 , and the PGA of IX ranges from 354 cm/s 2 to 707 cm/s 2 . It was found that most of the PGA larger than 178 cm/s 2 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 2 . The maximum PGA of 488 cm/s 2 was located northwest of Yangbi County. Since the simulated PGAs of only several points were bigger than 353 cm/s 2 , we propose that the biggest seismic intensity caused by the Yangbi event without accounting for these points was VIII. Remote Sens. 2021, 13, x FOR PEER REVIEW 12 of 18

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 nonlinear 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.

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. Remote Sens. 2021, 13, x FOR PEER REVIEW 13 of 18 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 −9 a −1 [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 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 pullapart 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 −9 a −1 [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ıg 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. 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, Rjb). 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. Figure 10. Comparison between the simulated and observed PGAs.

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 18 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

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 18 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 2 , which suggested that the seismic intensity was at least VIII caused by the 2021 Yangbi event.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/rs13204138/s1, Figure S1: Refined interferograms of ASC and DSC with atmospheric and orbital errors corrections, Figure S2: Downsampled InSAR observations, Figure S3: Statistics of the faulting geometric parameters obtained by non-linear inversion, Figure S4: The spatial relationship between the faulting geometry and the aftershocks, Figure S5: The tradeoff curve between the residual and the roughness, Figure S6: Simulated and residual interferograms generated by model-A, Figure S7: Simulated and residual interferograms generated by model-D,   Data Availability Statement: The European Space Agency owns the copyright of Sentinel-1 SAR data, and the Alaska Satellite Facility provides the downloading service through the https://search. asf.alaska.edu/#/, accessed on 1 June 2021). The strong motion records around the Yangbi area were collected from the CENC.