Abstract
On 3 July 2015, the Mw 6.5 Pishan earthquake occurred at the junction of the southwestern margin of the Tarim Basin and the northwestern margin of the Tibetan Plateau. To understand the seismogenic mechanism and the postseismic deformation behavior, we investigated the characteristics of the postseismic deformation fields in the seismic area, using 9 Sentinel1A TOPS synthetic aperture radar (SAR) images acquired from 18 July 2015 to 22 September 2016 with the Small Baseline Subset Interferometric SAR (SBASInSAR) technique. Postseismic LOS deformation displayed logarithmic behavior, and the temporal evolution of the postseismic deformation is consistent with the aftershock sequence. The main driving mechanism of nearfield postseismic displacement was most likely to be afterslip on the fault and the entire creep process consists of three creeping stages. Afterward, we used the steepest descent method to invert the afterslip evolution process and analyzed the relationship between postseismic afterslip and coseismic slip. The results witness that 447 days after the mainshock (22 September 2016), the afterslip was concentrated within one principal slip center. It was located 5–25 km along the fault strike, 0–10 km along with the fault dip, with a cumulative peak slip of 0.18 m. The 447 days afterslip seismic moment was approximately 2.65 × 10^{17} N m, accounting for approximately 4.1% of the coseismic geodetic moment. The deep afterslip revealed that a creeping process from steadystate “secondary” creeping to accelerating “tertiary” creep in the deep of fault. The future seismic hazard deserves further attention and research.
Introduction
At 09:07 UTC + 8 on 3 July 2015, an M_{W} 6.5 earthquake occurred at Pishan in Xinjiang (Fig. 1). The epicenter was located in the West Kunlun fault belt, at the junction of the southwestern margin of the Tarim Basin and the northwestern margin of the Tibetan Plateau^{1,2}. Since the Kunlun Mountains Mw 8.1 earthquake occurred in 2001, a series of extremely destructive earthquakes have occurred around the Tibetan Plateau, such as the Mw 7.3 earthquake in the Yutian area in 2008 and 2014^{3}, respectively. Both the earthquakes and the Pishan earthquake hit in the junction zone between the west Kunlun fault and the Altyn Tagh fault. As shown in Fig. 1, the northern margin of the Tibetan Plateau in western China reveals the profile of the Tarim Basin^{1}. The western Kunlun Fault reaches the northwest, the Altyn Tagh Fault reach the southeast encompass the Tarim Basin^{4}. The strikeslip Karakax fault was characterized by verging on the southern flank of the Tarim Basin^{4,5}. This large strikeslip fault was considered by some scholars to be the western continuation of the Altyn Tagh fault, to sever through the foundation of the crust^{4,5,6}. The Poskim fault was considered to be a hidden fault in the mountain front^{3}. The Hetian fault belt in the southwest of the Tarim Basin, extending more than 300 km, was taken for one of the primary fault belts^{6}. The Western Kunlun fault is characterized by tectonic deformation associated with reverse faults and folds, which coordinates with the relative movement between the Pamir front thrust belt and the West KunlunTarim Block^{7,8}. Zubovich et al.^{9} demonstrated that the Western Kunlun fault tectonic belt was dominated by compression, shortening, and deformation, and the dextral strikeslip rate was only 1–3 mm/a, forming a typical thinskinned nappe tectonic system.
Based on the solved focal mechanism that a thrust rupture triggered the M_{W} 6.5 Pishan mainshock, for the Global Centroid Moment Tensor (CMT) Catalog, aftershocks of M_{W} 4.5 and M_{W} 4.6 successively struck 2 h after the mainshock. Zhang et al.^{10} used the digital broadband seismic data recorded by Xinjiang network stations to obtain the aftershock sequence, which extended unilaterally along an NWW direction with a spreading length of approximately 50 km. The mainshock fault inclines towards the SW, showing the characteristics of a shovel thrust fault^{10}. Ainscoe et al.^{11,12,13} used InSAR data and seismic waveforms to solve the fault parameters, to inverse the slip distribution of the coseismic deformation. It was revealed that the earthquake was induced by a preexisting slope, reverse‐dipping toward the Tibetan Plateau, the range of ∼9–13 km in deep blind fault^{12,13}. The Pishan earthquake triggered the root faults to produce small swarms, but the previously accumulated tectonic stress was not completely released^{14,15}. The future seismic hazard deserves further attention.
Since the coseismic displacements recorded by InSAR include the influence of M_{W} 4.6 and M_{W} 4.5 aftershocks, the study of the coseismic deformation and fracture mechanism of Pishan M_{W} 6.5 was insufficient to reveal the intrinsic dynamic correlation of mainshock and aftershocks.
Therefore, to further probe into the seismogenic condition and postseismic deformation mechanism, we adopted Sentinel1A data and employed the Small Baseline Subset InSAR (SBASInSAR) technique with the correction of atmospheric error and orbit error to extract the mainshock deformation within 447 days. Then the SDM method was applied to invert the afterslip distribution. The postseismic deformation transmits pivotal feedback in the temporal evolution of the afterslip. A comparison of time evolution between the afterslip and the aftershocks was conducted to understand the mechanism of aftershock triggering. We comprehensively studied the relationship between afterslip and coseismic slip to examine regional tectonic movement.
Results
Postseismic deformation time series
The entire 9 SAR images were gathered from descending Path136 (hereafter referred to as D136). The earliest image acquisition time is 15 days after the earthquake, i.e., 18 July 2015, and the latest image is 22 September 2016. We calculated a time series and referenced it in time to the image acquired 15 days after the mainshock. Figure 2 shows the postseismic deformation time series. Positive values indicate that the deformation approached the satellite in the direction of lineofsight (LOS); negative values denote deformation in LOS motion turned aside away from the satellite. The aftershock dataset (M_{W} \(\ge\) 3.0) after the mainshock within 480 days was assembled. Figure 2 plots the spatial distribution of the aftershock sequence. The aftershock distribution is most dense in the NWW direction (see Fig S1). The bottom of rupture fault portion remained locked that induced the Pishan earthquake^{15}. A strand of aftershocks following the mainshock primarily focused on the Pishan anticline, which mirrored seismic events evoked by enhancing stress forbye the locked fault portion^{10,13,15}.
We selected 9 dates of accumulated LOS deformation within 447 days of the mainshock, as shown in Fig. 2. The accumulated postseismic deformation improved over time. At 20151115 (135 days), InSAR has shown preponderance in measuring the apparent displacement signal after the mainshock, and its LOSdirection maximum positive values around 18 mm and negative values around 14 mm. At the 20160922 (447 days), postseismic displacement features were most prominent, with significant positive values of up to 25 mm and negative values of up to − 20 mm in the LOSdirection. Positive range occurred around the epicenter and negative range occurred in the southern far away from the epicenter. By analyzing the postearthquake deformation field, we found that the deformation mainly occurred in the fault plane.
Three patches (A, B, and C in Fig. 2h) were selected from the deformation timeseries. To obtain a smoothed time series of deformation, we employed a logarithmic function to fit the discrete data of the three patches^{17}:
where \(y\) denote the accumulated LOS deformation (mm), \(a\) show the logarithmic decay for amplitude, and \(t\) denote days of the postearthquake time. Figure 3 depicts the curve fitting lines by a logarithmic function for patches A, B, and C.
A histogram in Fig. 3 shows aftershocks (≥ Mw 3.0) 480 days after the mainshock. As of August 3, 2015, 98 aftershocks(≥ Mw 3.0) occurred; subsequently, the rate of seismic events decreased sharply. Similarly, the timeseries cumulative deformation shows the law of logarithmic function decay. The similarity between the InSAR observations and aftershock activity may both be driven by the same process. There are three main mechanisms for the deformation of the crust postearthquake: afterslip, poroelastic rebound, and viscoelastic relaxation^{17}. Viscoelastic relaxation is the deformation, happened at the middle and lower crust, attribute to the viscous flow of material. Afterslip is a creep slip dislocation on or around the earthquake generating fault caused by coseismic rupture. In the short term, viscoelastic relaxation, generally speaking, does not play a dominant role in nearfield displacement after the earthquake^{17,18,19}. It is considered as poroelastic rebound when the direction of postseismic motion is opposite to the motion of coseismic, or displacement of postearthquake decreased with the extension of time^{20}. After the earthquake, we considered that the nearfield displacement increased with the extension of time, which was contrary to the deformation characteristics caused by the poroelastic rebound. The above features are most likely caused by afterslip^{20}. It is generally believed that deformation is adapted by cracked (fragile) processes over all the upper crust^{21}. Whereas, ductile (plastic) deformation primarily throughout the lower crust, governed by thermallyinduced dislocation creep^{21}. The mainshock was a characteristic reverse thrust belt earthquake, not only induced folding of the upper crust but also leading to aftershocks on the Pishan anticlines^{15}. The postseismic deformation mainly occurred on the fault plane was a manifest of the stress change caused by coseismic rupture. According to the phenomenon of afterslip and aftershocks, the temporal evolution of both are corresponding, Yang et al.^{17} believed that aftershocks were possessed by afterslip. Perfettini et al.^{21} found that postseismic displacement and the speed of seismic events show a coetaneous decline in the region of aftershocks. For the upper crust, it was considered that the afterslip dominated the elastic deformation, giving rise to displacement change observed on the ground^{21}. The temporal evolution of aftershocks and the displacement feature in the satellite LOS direction was corresponding with previous studies^{17,19,21}. Similarly, we conceived that afterslip undertakes a major impact, in the short term, on nearfield displacement after the Pishan earthquake.
Characteristics of postseismic afterslip distribution
From July 18, 2015, to September 22, 2016, We calculated the cumulative afterslip at 8 intervals, as shown in Fig. 4. We determined the slip direction by solving for slip along two vectors (strike vector and dip vector). The main feature of afterslip is thrusting, which is similar to the characteristics of the coseismic slip^{12,14}. According to the analysis of the afterslip (Fig. 4), the afterslip activity indicated that the shallow slip accumulated gradually up from days 87 to 447 of the mainshock in 0–10 km along with the fault dip, and the deep SSE slip in 30–40 km along with the fault dip accumulates more gradually over the entire period.
Theoretical and academic fault models were derived, some authors have penetrated research and performed triaxial “creep tests” to deduce the entire creep process^{21,22,23}. There has been proved that the cycle of primary, secondary, and tertiary creep, repetitively, throughout the overall creep process^{24}. The first period of “primary” creep is undergoing decelerating promptly, followed by stationary state “secondary” creep for a longer time, up to the accelerated process, “tertiary” creep eventually gives rise to a fault rupture by static fatigue; throughout the interseismic cycle, “secondary” semistable creep, in most cases, it is considered as the predominant process^{21}. We should note that the afterslip was in a steady state before 255 days, which corresponded to some steadystate “secondary” creep. The aftershocks occurred continuously on the fault plane during this period. The static fatigue is that the crack or rupture grows slowly with time under the sustaining loading of the stress, which is a subcritical state of fault fracture^{21,22,23,24}. When the fault slip comes through rateand state friction laws on account of imbalance intrinsic, or static fatigue brittle creep spread and extend on the fault zone, it was considered the potential origin for crustal deformationinduced transiently^{21}. The magnitude of afterslip in the shallow fault correspondingly increased from 399 to 447 days, which was consistent with some accelerating “tertiary” creep. The instability of this creeping process may lead to the aftershock on October 5, 2016. Our results also support the comment that aftershocks discontinuously or serially are induced prevailingly by afterslip^{24,25}.
At 447 days after the mainshock (September 22, 2016), the afterslip formed one slip center(Fig. 4). It was located 5–25 km along the fault strike, 0–10 km along with the fault dip, with a peak magnitude of 0.18 m. Along the fault strike at 18–35 km, along with the dip of 30–38 km, with a peak magnitude of 0.09 m. The simulated postearthquake deformation fields and residuals are shown in Fig. 5. The simulated postseismic deformation fields were similar to the postseismic deformation fields measured by InSAR. The simulation residuals of the main deformation zone after the mainshock were usually small and the larger residuals were mainly distributed in the northeast corners (Fig. 5), which might be related to the viscoelastic relaxation effect after the mainshock or to tectonic motion changes.
Fig. S2 shows the time evolution of the cumulative postseismic moment released by afterslip and CENCrecorded aftershocks. The curve fit the InSAR moment using an exponential equation:
where \(k\), \(p\), \(c\) are constants, \(n\left( t \right)\) is the cumulative postseismic moment released by afterslip, and \(t\) is the time after the earthquake (day).
When the shear modulus is 30 GPa, the afterslip moment of 87 days (28 September 2015) after the mainshock was approximately 0.5 \(\times\) 10^{17} N m, the moment magnitude of statistical equivalently, for M_{W} 5.03 (Fig. S1); the afterslip moment of 447 days (September 22, 2016) after the mainshock was approximately 2.65 \(\times\) 10^{17} N m, amounted to the moment magnitude M_{W} 5.5 (Fig. S2). The rate of the seismic moment released by aftershocks from CENC increased instantly until decayed overtime. The afterslip moment release evolves with exponential decay. It was indicated the afterslip moment is greater than the seismic moment here, to a mass of the postseismic deformation is aseismic. Therefore, to explain the surface displacement estimated, for the upper crust, it is momentous due to deep afterslip dominates the elastic deformation^{21}.
Comprehensive analysis of postseismic afterslip and coseismic slip
Some scholars have made achievements in the study of the Pishan coseismic slip distribution^{12,14} and iteratively searched the model parameters and inverted slip distribution of the constructed model. The same Sentinel1A data as Wen et al.^{14} was used for coseismic deformation and coseismic slip inversion (Table S1). Two interferograms covering the region of the coseismic deformation, from the ascending Path56 (for short A56) to descending Path136 (for short D136) were employed, Table S1 showed the data contents. To obtain coseismic deformation interferograms, the twopass InSAR approach, it was utilized for data processing^{20}.
A method of acquiring upward and eastward units of surface deformation was imported, 2.5dimensional (2.5D)analysis technology, to process from ascending and descending LOS data ^{\* MERGEFORMAT}^{14}. \(U\) as the LOS displacements to can be expressed as:
where \(U_{e}\) represent eastward displacement, \(U_{u}\) represent upward displacement. Besides, in the east and vertical directions, separately, \(s_{e}\) and \(s_{u}\) show the direction cosines.
Figure 6 shows the surface deformation map (2.5D) with upward (a) and eastward (b) components. In the eastward direction (Fig. 6b), the main red displacement zone toward moved eastward, with the top value eastward motion up to 6 cm; the blue displacement zone in the northwest corner with the supreme value closed to 3 cm toward westward. In the upward direction, the maximum uplift displacement (red region in Fig. 6a) was 14 cm, and the maximum subsidence (blue region in Fig. 6a) was 3.3 cm. In the Gobi desert, the random and decorrelation noise interference of InSAR measurement was relatively serious, and the coseismic displacement error was large. We adopted the same crustal layering model Crust1.0 and fault model as postseismic afterslip inversion (Tables S2 and S3).
The fault was enlarged to 36 km strike direction, to 40 km downdip direction. The fault plane was subdivided, alongstrike and downdip separately, to 1 \(\times\) 2 km rectangular pieces. Then, the Lcurve method was adopted to solve the smoothing factor. It was eventually determined to be 0.04 (Fig. S3). The fine slip distribution, for the coseismic rupture fault of the Pishan Mw6.5 earthquake, was simulated by the SDM method, as shown in Fig. 7a. The SDM program calculated that the correlation coefficient was 99.3%, to balance the prediction and observation. To score the reliability, Gaussian noise is applied to the observation data to generate 100 perturbed datasets for the fault slip distribution^{11}. We calculated standard deviation of the slip distribution with 100 perturbed datasets, as shown in Fig. 7b. The average errors of slip are ~ 4 cm, and the maximum error appears in the middle of the model up to 6 cm. The slip distribution errors have no essential influence on the overall slip distribution, it was indicated that our consequences are credible in the paper.
At a dip depth of 24–27 km, 22–23 km width, the remarkable slip was embodied a distinct feature, to the top magnitude of 0.9 m. At 0–7 km of the upper layer closed to the value of ~ 0.2 m, there has a nearly uniform slip (Fig. 7a). However, an absence for the slip between 12.5 and 17.5 km (Fig. 7a) may be due to a slip rate deficit. Although no surface rupture was found by the Institute of Geology, claimed by China Earthquake Administration (CEA), Wu et al. \* MERGEFORMAT^{15} found several tensiletype ground fissures at the top of the seismic anticline, and all of the strikes were the same as the NW direction and the anticline. It is conceived that the amount of dislocation near the deep source is bound up with ground fissures and is gradually transformed into the fold of the stratum during propagation and expansion along the fault from stress in the transition. The geodetic moment was approximate 6.36 N\(\times 10^{18}\) m, calculated from the coseismic slip distribution, corresponding to M_{W} 6.5. From M_{W} 6.2 (USGS) to M_{W} 6.5 (CENC) for seismological computations were compared, matched group estimated in Table S2. The geodetic moment of the 447 days after the mainshock, calculated from the afterslip distribution was approximately 2.65 \(\times\) 10^{17} N⋅m, which accounts for approximately 4.1% of the coseismic geodetic moment.
The postseismic afterslip on September 22, 2016 (447 days after the mainshock) and the coseismic slip of the fault were superimposed, as shown in Fig. 7d. The peak magnitude of afterslip at 447 days after the earthquake is only 0.18 m (Fig. 7c), while the peak magnitude of the coseismic slip reaches 0.9 m (Fig. 7a). Most of the afterslip zone overlaps with the coseismic slip zone, but the slip momentum further strengthened based on the coseismic slip. Our results are consistent with a previous study, which reported that afterslip does not change the characteristics of the fault coseismic slip distribution^{26}. Compared with the characteristics of coseismic slip, the postseismic afterslip extended to the shallow part of the northcentral section of the fault, and the slip momentum was increased over time at the dip depth SSE of 30–38 km.
Summary
We presented a spatial–temporal distribution of postseismic deformations from the M_{W} 6.5 Pishan earthquake on 3 July 2015, which was derived from 9 Sentinel1A TOPS SAR images from 18 July 2015 to 22 September 2016 using the SBASInSAR technique. The evolution process of postseismic afterslip distribution was inverted by the SDM method. We comprehensively analyzed the relationship between postseismic slip and coseismic slip. The main conclusions are as follows:

1.
The postseismic deformation of the Pishan Ms6.5 earthquake was manifested as an uplift in the satellite lineofsight (LOS) direction in the western portion of the epicenter, which had a significant positive value of up to 25 mm in the satellite LOS direction and a negative value of up to 20 mm in the satellite LOS direction. As of August 3, 2015, 98 aftershocks of Mw \(\ge\) 3.0 occurred; subsequently, seismic events decreased sharply. The shape of LOS deformation shows the logarithmic function attenuation law in the time evolution.

2.
At 447 days after the mainshock (September 22, 2016), the afterslip fields form one slip center. It was located 5–25 km along the fault strike, 0–10 km along with the dip, with a peak magnitude of 0.18 m. The geodetic moment that was based on the afterslip distribution 447 days after the mainshock was approximate 2.65 N\(\times 10^{17}\) m, which was equivalent to the moment magnitude Mw 5.5, accounting for approximately 4.1% of the coseismic geodetic moment.

3.
The postseismic mechanism was a dominant updip of the seismogenic portion of a fault zone. The temporal evolution of the deformation feature in the satellite LOS direction and aftershocks supported the scenario that afterslip on the fault plane was the most likely mechanism. The temporal evolution characteristics of afterslip and aftershocks supported the standpoint that aftershocks are dominated primarily by afterslip.
Our study helps understand the characteristics of regional fault activity. It will lay a foundation for us that understanding the spatial distribution of fault locking, slip rate deficit and stress accumulation, to better realize mainshockaftershock activity^{27}. These will be the focus of future studies.
Discussion
It has been determined that the Pishan Mw 6.5 earthquake was the result of the Ndirection extrusion of the Tarim block in the Tibetan Plateau, and the seismogenic fault was concealed in the EarlyMiddle Pleistocene^{10,28}. Meanwhile, the results of InSAR observations show that uplift occurred in the western part of the epicenter, and subsidence occurred in the southern part of the epicenter. There are obvious differential deformations on both sides of the fault after the earthquake. The strong earthquake activity near the northwestern margin of the Tibetan Plateau (the West Kunlun orogenic belt) and the Tarim Basin has a certain dynamic relationship with the collision between the Indian Plate and the Eurasian plate, so, close attention should be paid to this region.
The coseismic slip inversion with different InSAR data^{11,12,14} depicted that the focal depth was mainly distributed in 7–15 km. Nonetheless, the postseismic afterslip extended to 0–10 km in the shallow part of the northern section, and the slip trend was increased overtime at the deep SSE of the fault (Fig. 7b). The slip rate deficit in the lock (or partial lock), which was converted into stress and gradually accumulated until an earthquake (or fault creep)^{25}. There was a slight slip between 15 and 17.5 km dip depths of the fault (Fig. 7d), which was related to the slip deficit caused by the tectonic deformation of reverse faultfold. Local tensile stress was generated at the transition part of fold strata to change the momentum of slip^{15}. The real strong seismogenic area was located on the Tekilik fault at the foot of the mountain in the western Kunlun Mountain, where the strain energy was in a locked state^{10,15}. The magnitude of the deep afterslip was increased gradually over time, implying that a creeping process from steadystate “secondary” creeping to accelerating “tertiary” creep in the deep of fault. The region may still trigger earthquakes, and the danger of future earthquakes deserves further attention and research.
Methods
Data
Thirteen descending Sentinel1A Terrain Observation with Progressive Scans (TOPS) SAR images (Cband) acquired from D136, from 18 July 2015 to 22 September 2016 covering the earthquake zone, were collected to estimate the LOS deformation time series. The Sentinel1A Interferometric Wide (IW) imaging mode uses the TOPS method to obtain three subbands with a total width of approximately 250 km and resolutions of 5 m × 20 m along with the range and azimuth directions^{29,30}. Compared to traditional scanning or spotlight imaging modes, TOPS mode can significantly improve ScanSAR and achieve better image quality^{31,32}. Precision ephemeris orbits (https://qc.sentinel1.eo.esa.int) were used for orbital corrections. Elimination of topographic phase effects using 90 m resolution (https://srtm.csi.cgiar.org) Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) provided by NASA.
Generation of multiple differential interferograms
We adopted the Small Baseline Subset InSAR (SBASInSAR) technique to process 9 Sentinel1A TOPS SAR images over the research area. SBASInSAR is a method that effectively mitigates decorrelation phenomena by analyzing distributed scatters (DS) with high coherence based on an appropriate combination of interferograms^{33,}^{\* MERGEFORMAT}^{34}. It combines pairs of images with multimaster image methods, obtains the time series of surface deformation of each set by using the least square estimation method, and obtains the minimum norm least square solution of the surface deformation rate between image sequences by using singular value decomposition.
Assuming that a set of \(N\) + 1 SAR images is obtained by time series \(\left( {t_{0} , \ldots ,t_{N} } \right)\) covering the study area, and coregistered with other SAR images using an arbitrary image as the main image. Finally, a total of \(M\) differential interferograms are obtained. If \(N\) is an odd number, \(M\) can be expressed as^{35}:
We computed a set of multilook (also referred to, as lowpass (LP) filtered) differential interferograms characterized by constraints on the maximum temporal and perpendicular baseline^{36}. We obtained interferograms involving data pairs with Doppler centroid differences within the SAR image Doppler bandwidth. Twentyeight groups of small baseline interferometric pairs were selected, and the parameters of the interferometric pairs are shown in Table S4. The time baseline distance of all interferograms was less than 300 days. All interferograms are less than 110 m in terms of normal baseline (Fig. 8).
Atmospheric error correction
Atmospheric components have a strong correlation at short distances^{37}. To estimate the atmospheric disturbances, we select high phase coherence (HPC) pixels on the image grid to estimate the atmospheric phase screen (APS), which can be interpolated on the grid. Both operations (filtering and resampling) can be performed at the same time using kriging interpolation^{38}.
Estimate the mean value of the estimated atmospheric components in the \(H\) HPC pixels of differential interferograms as^{38}:
where \(APS_{master}^{T}\) is an estimation of the atmospheric phase contribution relative to the master image. \(K\) the differential interferograms. \(\widehat{{E_{atmo}^{^{\prime}} }}\) contains the residues that include atmospheric effects. \(\widehat{{\underline {a} }}\left[ {K \times 1} \right]\) are constant phase values, \(\underline{{\hat{p}_{\xi } }}\)[\(K \times 1\)] and \(\underline{{\hat{p}_{\eta } }}\)[\(K \times 1\)] contain the slope values of the linear phase components, along the azimuth \(\underline{\xi }\)[\(H \times 1\)] and slant range \(\underline{\eta }\)[\(H \times 1\)] direction, due to the atmospheric phase contribution.
The phase of each slave image can be modified by these estimated quantities. The new set of phases of the modified slave images will be represented as^{38}:
where \( \underline {\Psi } \left[ {K \times H} \right]\) contains the phases of the H HPCs as seen by the \(K\) slave images, \(\underline{{\psi_{m} }}\)[\(H \times 1\)] contains the phases of the H HPCs as seen by the master image, compensated for \(APS_{master}^{T}\), and \(E^{\prime \prime } [K \times H]\) is the residue matrix.
Orbit error correction
Although the orbital parameter correction has been performed on the Sentinel1A data using the precision ephemeris orbit data before the data interference processing, there are still residual orbit errors. In this paper, the polynomial optimization method is used to remove the residual orbit errors based on the atmospheric delay error correction^{39}. Estimate the orbital error phase of a single interference image pair as:
where \({ }ij\) represents the SAR image number that constitutes the interference image pair, \(\left( {x,y} \right)\) represents the pixel position of the interferometric image, and \(a_{ij}\), \(b_{ij}\), \(c_{ij}\), \(d_{ij}\) represent the orbit error correction coefficients of the interference image pair \(ij\). Finally, the orbit error correction coefficient is corrected by the leastsquares method to obtain the optimal solution. Fig. S4 shows an example of atmospheric and orbital error correction for interferometry.
Time series retrieval
The deformation rate was estimated using the leastsquares method^{40}. Assuming that the terrain phase has been removed by DEM, the interference phase on the ith point target can be expressed as^{41,42}:
where \(wrap( \cdot )\) represents the phase wrap operator; \(\varphi_{i}^{k}\) represents the interferometric phase observation on the point target; \(\varphi_{i,atmo}^{k} \) indicates atmospheric unevenness interference phase; \(\varphi_{i,orb}^{k}\) is the interference phase caused by the baseline error; \(\varphi_{i,topo  res}^{k}\) is the interference phase caused by the DEM residual; and \(\varphi_{i,noise}^{k} \) indicates a loss of coherent noise.
The phase of the interference \(\varphi_{i,defo}^{k}\) can be expressed as a polynomial model, as shown in the following equation^{42,43}:
where \(u_{i,p}\) correspond to linear rate acceleration and cubic velocity terms in the Norder polynomial deformation model, respectively; \(t_{i}^{k}\) is the image time difference of the kth interferogram. λ is the radar wavelength. \(\varphi_{i,res  def}^{k}\) is the residual phase of the interferogram \(i\) after the Norder polynomial simulation; it will be a small quantity. Theoretically, ground deformation is seen as a continuous function over time and can be uniformly approximated by polynomials^{41}.
There is a linear relationship between \(\varphi_{i,topo  res}^{k}\) and DEM residual \(\Delta h_{i}\), which can be expressed as^{44}:
where λ is the radar wavelength; \(B_{ \bot ,i}^{k}\) is the vertical baseline distance of the kth interferogram, and r is the slant distance between the ground point and sensor; \(\beta_{i}\) is the angle of incidence.
The interference phase model can be further expressed as:
The \(\varphi_{i,noise}^{k}\) can be eliminated by the phase difference operation between two adjacent point targets^{42}. The highorder polynomial simulates the deformation process;\(\varphi_{i,res  def}^{k}\) will be small, which will ensure that Eq. (9) is established and the correct model solution is obtained. In general, a thirdorder polynomial is a robust choice^{42}. We obtain the deformation time series by solving the above interferometric phase model.
Inversion of postseismic afterslip distribution
To reduce the amount of inversion calculation and improve the signaltonoise ratio of the data, we first performed a quadtree downsampling processing on the postseismic deformation fields for the 8 dates, and the number of deformation points is reduced from 21,927 to 967. Each date of the postseismic deformation fields had the same number of observation points after downsampling, which is beneficial for comparison in the afterslip inversion. In the case of quadtree downsampling, the sampling in the nearfield region was relatively dense, and the sampling in the farfield region was relatively sparse; therefore, the highresolution information of the deformation field near the epicenter was obtained to the greatest extent.
The fault model is an important parameter for postseismic afterslip inversion^{26}. In general, the fault parameters of the afterslip inversion can refer to the fault parameters of the coseismic deformation nonlinear inversion^{43}. The parameters of the fault model were set by referring to the existing coseismic fault model parameters (Table S2) and combining the geometric image features of the coseismic deformation fields and the postseismic deformation fields.
Assume that the fault plane passes through the epicenter and its upper boundary reaches the surface. The length of the fault (in the direction of the strike) is extended to 36 km, and the lower boundary is extended downward (in the direction of the dip) to 40 km along the fault plane. According to the focal mechanism solutions given by USGS, GCMT, and CENC, the range of the fault strike change is set between 95° and 125°, the range of the dip change is set to 18° to 35°, and a search calculation is performed every 1°. After many inversion experiments, the combination of strike and dip angles with minimum root mean square error of fitting residual and the highest modulus correlation coefficient was selected as the optimal value^{44}. Finally, the strike and dip angles of the coseismic faults were selected as 115° and 25°, respectively.
The SDM program developed by Professor Wang Rongjiang is used to constrain the postearthquake afterslip. The SDM program uses the steepest descent method for inversion calculation, which is an iterative algorithm for constrained leastsquares optimization^{45}. When the fault geometry parameters are determined, the slip parameters on the fault plane and the InSAR observed deformation values can be transformed into general linear problems. The expression is:
where \(d\) is the LOS InSAR observation; \(G\) is the Green's function calculated using the dislocation theory according to the elastic semiinfinite space model; \(m\) is the momentum of slip on the fault plane; \(\varepsilon\) represents the error of the observed data.
To obtain a more detailed fault slip distribution, the fault plane was divided into discrete rectangular elements of 1 km \(\times\) 2 alongstrike and downdip, and the medium Poisson's ratio parameter was set to 0.25. We adopted the weight ratio using the residual root mean square (RMS) of each dataset to weigh the contributions of two datasets, and the final weight ratio was 1:0.75. To avoid oscillation or excessive smoothing of the inversion results, Laplace smoothing constraints were added:
where \(m\) represents the amount of slip on each subfault plane of the fault, \(G\) is the Green's function, \(y\) represents the LOS InSAR observation, \({\text{a}}\) represents the smoothing factor, weighing the fit of the observed data and the roughness of the inversion result, \(H\) denotes the finite difference approximation expression of the Laplacian operator multiplied by a weighting factor proportional to the relevant slip, and \(\tau \) indicates the shear stress drop linearly related to the slip distribution on the fault plane.
The smoothing constraint was applied to the slip distribution to overcome the instability of the inversion results^{20}. The smoothing factor is generally obtained by plotting the compromise curve between the roughness and the misfit of the test data. The determined optimum smoothing factor was 0.04 using an Lcurve plot (Fig. S3). Meanwhile, the layered earth structure model derived from Crust1.0 A^{46} model was used for fault slip inversion, and the model parameters are shown in Table S3. As a consequence, the fault model was used to invert the afterslip distribution of the remaining 8 dates of the postseismic deformation fields.
References
 1.
Lu, R. et al. Coseismic and blind fault of the 2015 Pishan Mw 6.5 earthquake: Implications for the sedimentarytectonic framework of the western Kunlun Mountains, northern Tibetan Plateau. Tectonics 35, 956–964 (2016).
 2.
Li, S. & Mooney, W. Crustal structure of China from deep sounding profiles. Tectonophysics 288, 105–113 (1998).
 3.
Shan, X. et al. The surface rupture zone and coseismic deformation produced by the Yutian Ms7.3 earthquake of 21 March 2008 Xinjiang. Acta Geologica Sinica 86, 256–265 (2012).
 4.
Li, X. & Wang, K. On tectonic transformation of mountains and basins at northern margin of west Kunlun Mountains. Chin. J. Xinjiang Geol. 20, 19–25 (2002).
 5.
Huang, G.C.D., Roecker, S. W. & Levin, V. Lowercrustal earthquakes in the west Kunlun range. Geophys. Res. Lett. 38, 314–318 (2011).
 6.
Liu, W. Characteristics of Hetian fault belt and its oil potential. Chin. J. Oil Gas Geol. 6, 326–334 (1985).
 7.
Cowgill, E. S. Cenozoic rightslip faulting along the eastern margin of the Pamir salient, northwestern China. Geol. Soc. Am. Bull. 122, 145–161 (2010).
 8.
Edward, R. S. et al. Manfred. Late miocene–pliocene deceleration of dextral slip between Pamir and Tarim: implications for Pamir orogenesis. Earth Planetary Sci. Lett. 304, 369–378 (2011).
 9.
Zubovich, A. V. et al. GPS velocity field for the Tien Shan and surrounding regions. Tectonics. 29, 14–37 (2010).
 10.
GuangWei, Z., HongYan, Z. & ChangQing, S. Mechanism of the 2015 Pishan, Xinjiang, Ms 6.5 mainshock and relocation of its aftershock sequences. Chin. J. Seismol. Geol. 38, 711–720 (2016).
 11.
Ainscoe, E. A. et al. Blind thrusting, surface folding, and the development of geological structure in the Mw 6.3 2015 Pishan (China) earthquake. J. Geophys. Res. Solid Earth. 122, 9359–9382 (2017).
 12.
Li, Y., Luo, Y., Zhang, J. & Jiang, W. The 2015 Mw 64 Pishan earthquake, China: geodetic modelling inferred from Sentinel1A TOPS interferometry. Surv. Rev. 50, 522–530 (2017).
 13.
Zhang, G. et al. Blind thrust rupture of the 2015 Mw 6.4 Pishan earthquake in the northwest Tibetan Plateau by joint inversion of InSAR and seismic data. J. Asian Earth Sci. 132, 118–128 (2016).
 14.
Yangmao, W., Caijun, X., Yang, L. & Guoyan, J. Deformation and source parameters of the 2015 Mw 65 earthquake in Pishan, western China, from Sentinel1A and Alos2 Data. Remote Sens. 8, 134 (2016).
 15.
ChuanYong, W. et al. A folding earthquake occurred in front of the west Kunlun Mountains. Chin. J. Earthq. Geol. 39, 342–355 (2017).
 16.
China Earthquake Network Center. Available online: https://www.ceic.ac.cn (Accessed on 20 September 2019).
 17.
Yang, C. et al. Co and postseismic deformation mechanisms of the Mw 73 Iran earthquake (2017) revealed by Sentinel1 InSAR observations. Remote Sens. 11, 418–435 (2019).
 18.
Peltzer, G., Rosen, P., Rogez, F. & Hudnut, K. Poroelastic rebound along the Landers 1992 earthquake surface rupture. J. Geophys. Res. Solid Earth. 103, 30131–30145 (1998).
 19.
Bie, L., Ryder, I., Nippress, S. E. J. & Burgmann, R. Coseismic and postseismic activity associated with the 2008 Mw 6.3 Damxung earthquake, Tibet, constrained by InSAR. Geophys. J. Int. 196, 788–803 (2014).
 20.
Qu, C. Y. et al. Coseismic and postseismic deformation fields of the 2010 Yushu, Qinghai Ms 7.1 earthquake and their evolution processes. Chin. J. Geophys. 56, 2280–2291 (2013).
 21.
Perfettini, H. & Avouac, J. P. Postseismic relaxation driven by brittle creep: A possible mechanism to reconcile geodetic measurements and the decay rate of aftershocks, application to the ChiChi earthquake, Taiwan. J. Geophys. Res. Solid Earth. 109, B02304 (2004).
 22.
Marone, C. J., Raleigh, C. B. & Scholz, C. H. Frictional behavior and constitutive modeling of simulated fault gouge. J. Geophys. Res. Solid Earth. 95, 7007–7025 (1990).
 23.
Main, I. G. A damage mechanics model for powerlaw creep and earthquake aftershock and foreshock sequences. Geophys. J. Int. 142, 151–161 (2000).
 24.
Fish, A. M. Thermodynamical model of creep at constant stress and constant strain rate. Cold Reg. Sci. Technol. 45, 143–161 (1984).
 25.
Peng, Z. & Zhao, P. Migration of early aftershocks following the 2004 Parkfield earthquake. Nat. Geosci. 2, 877–881 (2009).
 26.
Hong, S., Dong, Y., Meng, G., Zhang, K. & Chen, L. Postseismic deformation extraction and afterslip inversion of the October 2008 Damxung Mw 6.3 earthquake. Chin. J. Geophys. 61, 4827–4837 (2018).
 27.
Jolivet, R. et al. Spatiotemporal evolution of aseismic slip along the Haiyuan Fault, China: implications for fault frictional properties. Earth Planetary Sci. Lett. 377–378, 23–33 (2013).
 28.
Ainscoe, E. A. et al. Blind thrusting, surface folding, and the development of geological structure in the Mw 6.3 2015 Pishan (China) earthquake. J. Geophys. Res. Solid Earth. 122, 9359–9382 (2017).
 29.
LiHua, F., JianPing, Wu. U., WeiLai, W., Ting, Y. & ChangZai, W. Relocation of the 2014 Ms7.3 earthquake sequence in Yutian Xinjiang. Chin. J. Geophys. 58, 802–808 (2015).
 30.
Sowter, A. et al. Mexico city land subsidence in 2014–2015 with Sentinel1 IW TOPS: Results using the intermittent SBAS (ISBAS) technique. Int. J. Appl. Earth Observ. Geoinf. 52, 230–242 (2016).
 31.
Wang, J., Shan, X., Zhang, G., Lv, J. & Shen, X. Inversion of coseismic deformation field and fault slip distribution of InSAR for the 2017 Jiuzhaigou Ms 7.0 earthquake. Chin. J. North China Earthq. Sci. 36, 1–7 (2018).
 32.
Ruiz Rodon, J., Broquetas, A., Gonzalez Arbesu, J. M., Closa, J. & Labriola, M. Signaltonoise ratio equalization for TOPSAR mode using a nonuniform steering rate. IEEE Geosci. Remote Sens. Lett. 9, 199–203 (2012).
 33.
Meta, A. J., Mittermayer, U., Steinbrecher, P., PratsIraola, P. Investigations on the TOPSAR acquisition mode with TerraSARX. Available online: https://www.researchgate.net/publication/4307574_Investigations_on_the_TOPSAR_acquisition_mode_with_TerraSARX (2007).
 34.
Berardino, P., Fornaro, G., Lanari, R. & Sansosti, E. A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms. IEEE Trans. Geosci. Remote Sens. 40, 2375–2383 (2002).
 35.
Lanari, R. et al. A smallbaseline approach for investigating deformations on fullresolution differential SAR interferograms. IEEE Trans. Geosci. Remote Sens. 42, 1377–1386 (2004).
 36.
Zhou, L. et al. Wuhan surface subsidence analysis in 2015–2016 based on Sentinel1A Data by SBASInSAR. Remote Sensing 9, 982–1003 (2017).
 37.
Williams, S., Bock, Y. & Pang, P. Integrated satellite interferometry: Tropospheric noise, GPS estimates and implications for interferometric synthetic aperture radar products. J. Geophys. Res. 103, 27051–27067 (1998).
 38.
Ferretti, A., Prati, C. & Rocca, F. Permanents catterers in SAR interferometry. IEEE Trans. Geosci. Remote Sens. 38, 2202–2212 (2000).
 39.
Cavalié, O., Doin, M. P., Lasserre, C. & Briole, P. Ground motion measurement in the Lake Mead area, nevada, by differential synthetic aperture radar interferometry time series analysis: Probing the lithosphere rheological structure. J. Geophys. Res. 112, B03403 (2007).
 40.
Hooper, A. A multitemporal InSAR method incorporating both persistent scatterer and small baseline approaches. Geophys. Res. Lett. 35, L16302 (2008).
 41.
Ferretti, A., Prati, C. & Rocca, F. Nonlinear subsidence rate estimation using permanent scatterers in differential SAR interferometry. IEEE Trans. Geosci. Remote Sens. 38, 2202–2212 (2000).
 42.
Zhang, Y., Wu, H. & Sun, G. Deformation model of time series interferometric SAR techniques. Acta Geodaetica et Cartographica Sinica 41, 864–869 (2012).
 43.
Bie, L., Ryder, I. & Nippress, S. Seismic and aseismic activity associated with the 2008 Mw 63 Damxung earthquake. Tibet. Geophys. Res. Lett 15, 3288 (2013).
 44.
Wang, K. & Fialko, Y. Observations and modeling of coseismic and postseismic deformation due to the 2015 Mw78 Gorkha (Nepal)earthquake. J. Geophys. Res. Solid Earth. 123, 761–779 (2018).
 45.
Wang, R., Diao, F. & Andreas, H. Sdm  a geodetic inversion code incorporating with layered crust structure and curved fault geometry. Geophys. Res. Lett. 15, 2411 (2013).
 46.
Laske, G., Masters, G., Ma, Z. & Pasyanos, M. Update on CRUST10A1degree global model of Earth’s crust. Geophys. Res. Abstracts 15, 2685 (2015).
Acknowledgements
We acknowledgment the European Space Agency (ESA) for freely making available the Sentinel1A data, NASA for providing the SRTM3 DEM data, Maps were prepared using the Generic Mapping Tools (GMT) and MATLAB software.
Funding
This research was funded jointly by the National Key Research and Development Program of China (Grant Number 2019YFC1509201), the National Natural Science Foundation of China (Grant Number 41374028) and Fundamental Research Funds for the Central Universities (Grant Number 310826175018).
Author information
Affiliations
Contributions
All authors contributed to the manuscript and discussed the results. S.W., Y.Z., conceived and designed the experiments; S.W., Y.W., J.J., performed the experiments and drafted the manuscript; J.J., Z.J., M.H., contributed to the InSAR data analysis.
Corresponding authors
Ethics declarations
Competing interests
The authors declare no competing interests.
Additional information
Publisher's note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary information
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Wang, S., Zhang, Y., Wang, Y. et al. Postseismic deformation mechanism of the July 2015 MW 6.5 Pishan earthquake revealed by Sentinel1A InSAR observation. Sci Rep 10, 18536 (2020). https://doi.org/10.1038/s41598020752780
Received:
Accepted:
Published:
Further reading

Microanchored borehole fiber optics allows strain profiling of the shallow subsurface
Scientific Reports (2021)
Comments
By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.