3.1. Time-Series InSAR Processing Method for Coseismic Deformation Field Extraction
In this study, we used a single master image interferometric configuration, specifically conventional PS-InSAR, and applied a spatio-temporal algorithm for phase unwrap**. Because the extracted deformation signals were abrupt, the main improvement of the method suggested in this study was the segmented time-series process for estimating various errors. The phase of the coseismic deformation contained in the post-seismic time-series interferogram can affect error estimates, resulting in incorrect estimates that affect the extraction of the final coseismic deformation signal. Therefore, this method devides unwrapped interferograms into pre-seismic and post-seismic sequences. To estimate the spatially correlated look angle (SCLA) error and the atmospheric and orbit error (AOE) due to the master image, we only used the pre-seismic time-series interferogram for spatio-temporal filtering. This prevented the phase of the coseismic deformation in the post-seismic interferograms from affecting the error estimation. Finally, AOE due to the slave images were separated from the interferogram by performing temporal filtering on the pre-seismic and post-seismic interferograms, respectively. In a general analysis, filtering was performed uniformly for all time-series interferograms and weak coseismic deformation signals may be erroneously filtered out as noise, which must be accounted for in any processing flow. The complete process flow is illustrated in
Figure 2.
An SAR image acquired just before the earthquake was used as the master image, and co-registered with other pre-seismic and post-seismic images that were used as the slave images. The master and slave images were co-registered individually to perform differential interferometry. In this case, the coseismic deformation signal is not present in the pre-seismic interferogram but rather in the post-seismic interferogram. In StaMPS, the threshold is first set using the amplitude deviation method combined with phase stability estimates to select PS candidates, and then the wrapped phase is corrected to filter out spatially uncorrelated look angle errors. Finally, the phase of the PS candidates were unwrapped to recover the absolute deformation phase using the 3D spatio-temporal method. The phases in the unwrapped interferograms were divided into five terms:
where
indicates the unwrapped phase,
indicates the surface deformation phase,
indicates the error due to atmospheric effect,
is the residual orbital error,
is the residual phase due to the SCLA error, and
is the term for noise due to other factors, for example, co-registration and phase unwrap**.
Because the pre-seismic and post-seismic time-series interferograms needed to be processed separately, two different expressions were used for the unwrapped phases before and after the earthquake. The spatially uncorrelated look angle errors were screened out before the interferogram was unwrapped, so that the remaining residual errors in the unwrapped interferometric phase of the digital elevation model (DEM) only included the SCLA errors. The other spatially correlated errors were divided into AOEs due to the master and slave images, where the AOEs due to the master image are present in each interferogram. The deformation phase and other error terms in the pre-seismic and post-seismic interferograms can be expressed by the following equations:
where
denotes the master image,
represents the unwrapped phase of the pre-seismic interferogram,
represents the unwrapped phase of the post-seismic interferogram, and the superscripts m and s denote the contributions of the master and slave images to each spatially correlated error, respectively.
denotes the surface deformation phase that does not include coseismic deformation signals and is present in all interferograms, and
represents the coseismic deformation phases included in the post-seismic interferogram.
When using a DEM for terrain phase estimation, DEM errors can lead to other terrain-related residual phases in the interferometric phase. DEM errors tend to be partially spatially correlated and mapped to look angle errors; therefore, the SCLA error accounts for most DEM errors [
22]. An SCLA error is present in each pair of interferograms. Because the terrain undergoes substantial deformation after an earthquake, the post-seismic time-series interferograms containing the coseismic deformation signals were not involved in the estimation of the SCLA error in this study. The following equation represents the SCLA error due to the DEM:
where
denotes the radar wavelength;
denotes the baseline distance between the slave image and master image sensor position;
denotes the look angle in the master image geometry;
denotes the angle between the baseline and horizontal vectors; and
denotes the vertical component of the baseline.
Atmospheric conditions are unlikely to be identical when satellites make repeated observations of the ground; therefore, the path of the electromagnetic waves changes during propagation and produces phase delays during radar imaging, resulting in substantial errors in the measurement of surface deformation [
23,
24]. Both the positioning accuracy of the satellite orbit and atmospheric effects during operation can lead to biases when a satellite revisits the same position. Even with orbital parameters, it is not possible to completely remove orbital error streaks, and residual errors propagate with subsequent processing [
25]. Especially in coseismic deformation measurements, orbital errors propagate into the interferometric phase and affect the final estimation of the deformation. The orbital and atmospheric errors are similar in their spatial and temporal properties; therefore, we adopted spatiotemporal filtering to separate the coseismic deformation signal from the AOEs. The AOE due to the master image in the interferometric phase was estimated first. Considering the effect of post-earthquake coseismic deformation on the estimation, only pre-seismic interferograms were selected to estimate this error. Low frequencies characterize the nonlinear deformation signal in both the time and space domains, whereas high frequencies characterize the AOE in the time domain [
12]. Therefore, we performed low-pass filtering in the time dimension of the unwrapped phase for the pre-seismic interferogram to separate the surface deformation signal from the AOE due to the master image. The AOE due to the master image was removed from the interferogram, and an estimate of the master image AOE was obtained:
where
is the temporal low-pass filter operator.
To estimate the AOE due to slave images, separate spatio-temporal filtering was applyied to pre-seismic and post-seismic time-series interferograms. Given the high-frequency characteristics of the AOE in the time dimension, high-pass filtering was performed for the remaining phases in the pre-seismic and post-seismic time-series interferograms. While the AOE has low-frequency characteristics in the spatial dimension, the high-pass filter allows spatial low-pass filtering for the signal output, and the spatial filter outputs the final estimate of the AOE due to the slave image:
where
is the temporal high-pass filtering operator and
is the spatial low-pass filtering operator.
The phase of the interferogram formed by the SAR images on the nearest date before and after the earthquake date contains the most accurate coseismic deformation signal; therefore, the final coseismic deformation estimate is only based on the interferometric phase in this interferogram. The DEM-induced SCLA error, AOE due to the master image, and estimated AOE due to the slave image were first removed from the interferometric phase, and other elevation-related errors were then removed from the remaining phases. Finally, the extracted coseismic deformation phases were spatially filtered to remove other spatial noises and obtain a final coseismic deformation phase estimate with high precision. The coseismic deformation estimate was given by the following equation:
where
indicates the first post-seismic image to the date of the earthquake. As other surface deformation
is negligible,
represents the final estimate of coseismic deformation.
3.2. Inversion of Fault Geometry Parameters and Slip Distribution
Taking the coseismic deformation fields extracted by the above time-series InSAR processing method as constraints, we performed an inversion of the 2021 Yangbi M
W 6.1 earthquake to obtain the fault geometry parameters and fault slip distribution. Based on this, we studied the fault activity and rupture characteristics of this earthquake in detail. In this study, we used the Geodetic Bayesian Inversion Software (GBIS) developed by Bagnardi and Hooper in 2018 to perform nonlinear homogeneous and slip inversion of the best fault geometry parameters for this earthquake [
26]. GBIS uses a Bayesian probabilistic inversion algorithm for InSAR observations derived from the original datasets. In the Bayesian framework, an optimal set of source model parameters is extracted from the posterior probability density function (PDF) by quickly estimating the best model parameters and associated uncertainties through efficient sampling the posterior PDF. In the process of seismic inversion using this method, the Markov chain Monte Carlo (MCMC) method was used to perform sampling [
27], which was controlled in combination with the Metropolis Hastings algorithm [
28,
29]. Automatic step-size selection can control the random walk step of each model parameter to ensure sampling efficiency. After multiple iterations, the set of model parameters with the maximum a posterior probability after pre-assigned iterations was selected.
We inverted the fault based on the elastic half-space dislocation model proposed by Okada in 1985, which provides the relationship between the surface deformation and fault slip at depth [
30]. This elastic dislocation model provides nine source model parameters, including the length, width, depth, X and Y coordinates of the lower edge midpoint, dip and strike of the fault, and dip slip and strike slip. First, the seismic region was modeled as a uniform rectangular fault plane and the geometric parameters of the fault were accurately estimated using the GBIS algorithm. To obtain the dislocation distribution across a fault, it is necessary to invert the slip distribution across the fault plane. Because the slip distribution is nonuniform, the inversion of the fault slip distribution is equivalent to solving the optimal value of a linear problem. Based on the parameters of the fault geometric model, the fault plane was discretized into several uniform sub-fault planes, and the slip distribution of each sub-fault plane was then inverted to further obtain the fine slip distribution of the entire fault. The slip parameter model coefficient matrix is called Green’s function matrix, and we chose the constrained least-squares algorithm to estimate slip distribution [
31]. As Green’s function matrix is severely rank deficient, the Laplace smoothing matrix was added to the inversion for regularization, and the L-curve method was used to select the best regularization parameters [
32].