Application of singular spectrum analysis to InSAR time-series for constraining the post-seismic deformation due to moderate magnitude earthquakes: the case of 2019 M_w 6 Mirpur earthquake, NW Himalaya

SUMMARY

Detection and separation of the subtle post-seismic deformation signals associated with moderate magnitude earthquakes from interferometric synthetic aperture radar (InSAR) time-series is often challenging. Singular spectrum analysis (SSA) is a statistical non‐parametric technique used to decompose and reconstruct signals from complex time-series data. We show that the SSA analysis effectively distinguished the post-seismic signal associated with the 2019 M_w 6 Mirpur earthquake from periodic and noise components. The SSA-derived post-seismic deformation signal is smoother and fits better to an exponential model with a decay time of 34 d. The post-seismic deformation is confined to the southeast of the rupture area and lasted for ∼90 d following the main shock. Inversion of the post-seismic deformation suggests an afterslip mechanism with a maximum slip of ∼0.07 m on the shallow, updip portions of the Main Himalayan Thrust. The 2019 Mirpur earthquake and afterslip together released less than 12 per cent of the accumulated strain energy since the 1555 Kashmir earthquake and implies continued seismic hazard in the future.

1. INTRODUCTION

Post-seismic deformation is the transient response of the lithosphere to stress changes that occur after an earthquake. Post-seismic deformation can typically last from few weeks to decades (e.g. Gourmelen & Amelung 2005). It can result from various physical processes, including afterslip on the primary fault plane (Marone et al. 1991; Freed 2007), Viscoelastic relaxation of the lower crust and/or upper mantle (Savage & Prescott 1978; Li & Rice 1987; Pollitz et al. 2000), and poro-elastic rebound (Peltzer et al. 1998; Jónsson et al. 2003). Measurements of the spatio-temporal evolution of post-seismic deformation are important to understand the rheological and mechanical behavior of fault zones and thereby for earthquake hazard assessment. Continental earthquakes of large magnitude (M_w 7+) produce prominent post-seismic deformation (> 1 mm yr⁻¹) spanning a few months to several decades and could be measured using space geodetic techniques such as Global Positioning System (GPS) or interferometric synthetic aperture radar (InSAR, e.g. Pollitz et al. 2001; Burgmann et al. 2002; Ryder et al. 2007). However, accurate measurement of post-seismic deformation following a moderate-magnitude earthquake remains challenging due to lower signal-to-noise ratios (SNRs) and shorter duration compared to larger events (e.g. Satyabala et al. 2012; Wang et al. 2020; Furuya & Matsumoto 2022; Marconato et al. 2022). Advanced InSAR time-series analysis methods, such as the small baseline subset techniques (Berardino et al. 2002; Li et al. 2022), enable precise measurement of surface deformation over time. However, separation of the subtle post-seismic deformation from the composite InSAR signal consisting of seasonal/periodic components along with orbital and tropospheric noises is still challenging.

Singular spectrum analysis (SSA), a statistical data adaptive method provides a powerful alternative tool to extract information from noisy time-series without prior knowledge of the dynamics affecting them (Vautard & Ghil 1989). SSA operates as a principal component technique in the time domain, relying on the singular value decomposition of a specific matrix built upon the time-series. This approach decomposes the original series into interpretable components such as slowly varying trends, periodic components and noise. Though SSA-based techniques have been widely used in the fields of oceanography, climate science and geophysics (Ghil et al. 2002; Hassani & Zhigljavsky 2009; Korotchenko et al. 2014; Choliy 2015; Kumar et al. 2015), their application for InSAR signal decomposition is not well explored. Jiang et al. (2018) used SSA to extract seasonal InSAR deformation signals to characterize the aquifer properties. Recently, Sreejith et al. (2024) showed that the deformation signals related to a slow-moving landslide could be recovered from the seasonal, hydrologic and noise components by applying 1-D SSA to the InSAR time-series.

In this study, we employed the SSA technique to analyse post-seismic InSAR time-series data followed by the 2019 September 24, M_w 6, Mirpur earthquake (2019 Mirpur earthquake, hereafter) occurred in NW Himalaya (Figs 1 and 2). We show that SSA analysis effectively extracted the spatio-temporal post-seismic signals related to the 2019 Mirpur earthquake from periodic and noise components. We also show that the observed post-seismic deformation could be explained by an afterslip mechanism and the SSA analysis significantly improves the source model owing to the noise reduction.

2. 2019 MIRPUR EARTHQUAKE AND POST-SEISMIC DEFORMATION

On 2019 September 24, an unusual M_w 6 earthquake struck near the Mirpur city in the Kashmir Himalaya (Fig. 1). It ruptured the shallow, near-horizontal updip portion of the Main Himalayan Thrust (MHT), a detachment boundary between Indian and Eurasian plate, which absorbs roughly 12–14 mm yr⁻¹ convergence in NW Himalaya over 175 km width of the MHT (Kundu et al. 2014; Sreejith et al. 2021; Meyer et al. 2024). The MHT marks the boundary between the Indian Plate to the south and the Eurasian Plate to the north, and it is associated with significant seismic activity and mountain-building processes in the region (e.g. Bilham et al. 2001; Sreejith et al. 2018; Dal Zilio et al. 2021). Himalayan earthquakes generally originate at the lower part of the MHT, ∼100 km north of the MFT, and then move updip along the locked décollement or steeper fault planes (Sreejith et al. 2016; Bilham 2019). In contrast, the 2019 Mirpur Earthquake ruptured the frontal MHT with an extremely low dip angle of ∼3° (Sreejith et al. 2021; Xie et al. 2021; Barkat et al. 2022).

We use InSAR time-series data acquired in ascending and descending tracks from the Sentinel-1 satellite to constrain post-seismic deformation for six months following the 2019 Mirpur earthquake (Figs S1 and S2 and Table S1, Supporting Information). We use InSAR Scientific Computing Environment software (ISCE) to generate interferogram stacks (Rosen et al. 2012). Topographic phases were removed from each interferogram using Shuttle Radar Topography Mission (SRTM) data of 30 m spatial resolution (Farr et al. 2007). Small baseline subset interferometry based time‐series analysis was carried out using the MintPy software package (Berardino et al. 2002; Yunjun et al. 2019). We invert network of interferograms for raw phase time-series, taking the first Single Look Complex (SLC) image acquired after the event as reference. Atmospheric phase screens generated using the Generic Atmospheric Correction Online Service model (GACOS) (Yu et al. 2017) were applied to reduce tropospheric noise for all data sets. Maps showing cumulative line of sight (LOS) deformation of six months following the 2019 Mirpur earthquake derived using Sentinel-1 Ascending (T100) Descending (T107) data are presented in Figs S3(a) and (c) in the Supporting Information. We infer LOS deformation up to 18 and 25 mm (towards the satellite) along ascending and descending orbits towards the southeast of the coseismic rupture (Figs 3c and d). Time-series analysis suggests that the initial rapid deformation following the earthquake became stable by three months following the earthquake. The post-seismic InSAR time-series appears to be noisy, with large variations in displacements in space and time probably related to unmodelled atmospheric artifacts (Fig. 4 and Fig. S3, Supporting Information). However, by examining the time-series at several pixels, we infer that the systematic decaying deformation signals corresponding to the post-seismic period are only observed for ∼ 17 km² patch southeast of the main shock rupture (Fig. 3 and Fig. S3, Supporting Information).

3. SINGULAR SPECTRUM ANALYSIS

SSA is a statistical non-parametric method to probe information from noisy time-series without prior knowledge of the dynamics affecting the time-series (Vautard & Ghil 1989). SSA is essentially a principal component technique in the time domain that relies on the singular valued decomposition of a specific matrix constructed upon the time-series. SSA decomposes the original series into a sum of a small number of interpretable components such as a slowly varying trend, periodic components and noise. The SSA approach consists of two main steps: decomposition and reconstruction (Fig. 2). Each step comprises two separate phases. In the first step, the series is decomposed through embedding and singular value decomposition; in the second stage, the original series is reconstructed via grouping and diagonal averaging to utilize the reconstructed time-series for the trend analysis (Golyndina et al. 2001). These four phases are discussed in detail in Supporting Information Text 1. Sreejith et al. (2024) showed that the SSA analysis is particularly effective in separating non-linear signals from the InSAR 1-D time-series. In this study, we extend the SSA analysis to a multitemporal InSAR data set to separate the spatio-temporal displacement components from the InSAR stack (Fig. 2).

We decomposed the post-seismic InSAR time-series data of the 2019 Mirpur earthquake into seven components using the maximum theoretically allowed window length (L). A larger L value provides a more detailed decomposition into elementary components, offering better separability. The trade-off curve between the original signal and the number of components in the reconstructed signal suggests that the first two principal components are relevant and primarily contribute to the input series or eigenvalues (Fig. S4, Supporting Information). Plots displaying time-series of all seven principal components for both ascending and descending tracks are presented in Fig. S5 (Supporting Information). We interpret the first two components are representing transient post-seismic signal, while the remaining components represent seasonal variations and noise (Figs 3c and d). Interferograms comparing the transient post-seismic component of SSA with the original interferograms are shown in Fig. 4 and Fig. S3 (Supporting Information). We note that the post-seismic motion was completed in three months following the earthquake, with approximately half of it happening within 30 d.

4. POST-SEISMIC MODEL

In order to understand the post-seismic deformation characteristics, we tested the following exponential (eq. 1) and logarithmic functions (eq. 2) to fit with the observed displacement time-series (Freed 2007; Fattahi et al. 2015).

where y is the LOS displacement at given time, D denotes a static offset, a is a magnification factor, ∆t is the elapsed time since the earthquake and τ_l and τ_e denote the decay time coefficient of exponential and logarithmic functions respectively. We infer that the observed post-seismic time-series could be better explained with an exponential model with a lower RMSE (root mean square error, 0.7 mm for ascending and 0.62 mm for descending) compared to the logarithmic model (1.14 mm for ascending and 1.10 mm for descending; Figs 3c and d). The exponential model requires a decay time of 34 d to fit with the observed data. The exponential fit of post-seismic time-series typically indicates a viscoelastic relaxation mechanism at the lower crust or upper mantle (Savage & Prescott 1978). However, recent results suggest that a shallow afterslip could also result exponential decay characteristics in the post-seismic deformation time-series (e.g. Fattahi et al. 2015; Han et al. 2022; Sadeghi Chorsi et al. 2022). While shallow afterslip dominates the near-field post-seismic deformation, in the far-field, the observed deformations are mainly controlled by deeper viscoelastic processes (Marone et al. 1991). Hence, we suggest that the spatially localized deformation with exponential decay characteristics after the 2019 Mirpur earthquake could have been caused by an afterslip mechanism. We have simulated the post-seismic deformation pattern for the 2019 Mirpur earthquake assuming viscoelastic and poro-elastic models (Pollitz 1997; Jónsson et al. 2003). However, both models are insufficient to explain the pattern and magnitude of the post-seismic displacements (Supporting Information Text S2 and S3 and Fig. S6). These results clearly suggested that afterslip (fault creep or aseismic sliding) on the main fault plane is the primary contributor to post-seismic motion after the 2019 Mirpur earthquake.

We invert three months of transient post-seismic motion derived from SSA to model the afterslip on the causative fault of the 2019 Mirpur earthquake assuming linear superposition of rectangular dislocations in an elastic half-space (Okada 1985). We used steepest descent method iterative algorithm (Wang et al. 2013) for the constrained least-squares optimization to solve for the dip-slip and strike-slip components. We used a quadtree variance downsampling approach to reduce computational burden (Jónsson 2002). The afterslip inversion was carried out using the same coseismic fault geometry parameters of Sreejith et al. (2021). The width of the fault was set to 25 km, and the fault segment was divided into 1 × 1 km² patches. The optimum smoothing factor (γ = 0.2) for the inversion was determined from the trade-off curve between the model roughness and misfit (Fig. S7, Supporting Information). The modelled deformation matches well with the observations with an error of 2.7 and 3.1 mm for the ascending and descending tracks, respectively.

The inverted afterslip model suggests that the post-seismic relaxation is confined to the updip portions of the coseismic asperity with a maximum slip of 0.07 m with a rake angle of 100°. The afterslip was dominated by thrust motion, similar to the coseismic slip (Fig. 5). The total moment released by the afterslip (1.67 × 10¹⁷ N·m, M_w 5.45) only accounts for 13 per cent of the coseismic moment (Sreejith et al. 2021).

5. DISCUSSION

5.1 Applications of SSA technique for InSAR analysis

Development of techniques for separating geophysical displacements from noise components is extremely important for large-scale, multitemporal InSAR analysis. Mitigation of atmospheric contributions either using empirical methods or using weather models are challenging for short-wavelength and/or low-magnitude deformation signals (Doin et al. 2009). Another important challenge is to separate geophysical signals of interest from the composite signals from post-seismic processes, fault creep, landslides, or hydrologic loading without having priory constraints on the characteristics of these signals. Independent component analysis (ICA) has been explored for separating geophysical signals from atmospheric noise by using cluster analysis of the spatial patterns of independent components (Ebmeier 2016). Though ICA is effective in detecting spatially or temporally independent components of the mixed signals, the magnitudes of independent components are ambiguous and cannot be ranked with respect to the amount of variance explained by them (Hyvärinen & Oja 2000). Secondly, ICA assumes that the constituent signals are temporally independent and are not sensitive to temporally uncorrelated signals such as oscillatory signals with varying phases and periods. Unlike ICA, the magnitude of extracted signal components in SSA are not ambiguous and does not require scaling or normalization. Hence, SSA provides a unique method for decomposing multitemporal InSAR data to tectonic, non-tectonic, periodic and noise components for quantitative interpretation.

5.2 Impact of SSA in InSAR post-seismic deformation analysis

In this study, the SSA analysis on the InSAR time-series could separate the subtle post-seismic transient deformation from the 2019 Mirpur earthquake from the noise components (Figs 3 and 4). In order to understand the effect of SSA analysis we compared the standard deviations of InSAR time-series (Fig. S8, Supporting Information). We observe that the SSA analysis resulted upto 25 per cent reduction in standard deviations in the InSAR time-series maps, suggesting a significant improvement in the SNR. This is also evident from the reduction in random noise as observed in the time-series at pixels P1 and P2 in Fig. 3. The SSA-derived post-seismic transients fit well with model functions (Fig. 3) with a lower RMSE (0.7–1.1 mm) and higher R² (0.96–0.99) compared to the original InSAR time-series (RMSE = 2.47–2.86 mm and R² = 0.8–0.9).

To further understand the influence on SSA analysis in post-seismic slip models we have compared after slip models inverted using three months of cumulative displacements derived using SSA and with that of the original InSAR time-series (Fig. 5c). We note that the SSA-derived afterslip model is more compact with slip localized to the SE part of the fault plane. The afterslip model showed ∼10 per cent difference in the slip inverted from the deformation compared to the afterslip model without SSA.

5.3 Post-seismic slip distribution along MHT and seismic hazard

It is well known that afterslip triggered by coseismic stress changes are confined to velocity-strengthening regions within the fault zone (Marone et al. 1991). Occurrence of afterslip updip of the coseismic rupture of the 2019 Mirpur earthquake is consistent with the coseismic Coulomb stress change model (Sreejith et al. 2021). However, the shallow portion of the MHT near Mirpur region is believed to be locked indicating velocity weakening characteristics (Meyer et al. 2024). The observed afterslip updip of the coseismic rupture and its concentration towards southeast portion (Fig. 5b) could be a result of heterogeneity in the rheology (Scholz 1998; Mahsas et al. 2008) and/or local variations in the degree of locking along the shallow MHT (Satyabala et al. 2012).

It is interesting to note that the duration of post-seismic slip associated with the 2019 Mirpur (90 d) is significantly less than > 1-yr afterslip reported for most of the small to moderate magnitude (M_w = 5–6) earthquakes (e.g. Furuya & Satyabala 2008; Murray-Moraleda & Simpson 2009; Bell et al. 2012; Fattahi et al. 2015). Furuya & Matsumoto (2022) observed that the long-lasted post-seismic slip for M_w 5–6 events results larger post-seismic-to-coseismic moment ratios (PCMR) of > 30 per cent than the M_w 6 + events (∼30 per cent). Available PCMR values for thrust, strike-slip and normal earthquakes of magnitudes ranging from 5 to 9.2 complied from previous studies are presented in Fig. 5(d). The observed PCMR for the 2019 Mirpur earthquake (13 per cent) is well within the range of earthquake of similar magnitude earthquakes with thrust mechanisms and shorter after slip durations (Fig. 5d). However, strike-slip earthquakes of magnitude 6 or lower have afterslip for significantly longer periods (1 + yr) and hence larger PCMR (75–560 per cent).

The short-lived afterslip following the 2019 Mirpur earthquake implies a rapid relocking of the coseismic rupture area and surrounding regions at the MHT. The 2019 Mirpur earthquake and the afterslip together released only 12 per cent of the accumulated strain on the décollement since the 1555 Kashmir earthquake, indicating the potential for future seismic events in the region. However, whether the accumulated strain energy will be released by a series of earthquakes with magnitude between 6 and 7 causing partial ruptures, or with an M_w 7.5 + event leading to the complete rupture reaching MFT (Sreejith et al. 2018; Dal Zilio et al. 2019) is unclear. Bilham (2019) proposed that the relic strain energy from incomplete ruptures of the MHT could fuel a larger magnitude earthquake.

6. CONCLUSIONS

InSAR time-series analysis using SSA method is a powerful technique to extract subtle post-seismic deformation signals for moderate magnitude earthquake. We show that the SSA analysis successfully distinguished a short-lived (90 d) post-seismic signal following the 2019 M_w 6 Mirpur earthquake from periodic and noise components. We infer that the SSA analysis enhanced the post-seismic deformation pattern in the time-series InSAR images with a significant noise reduction upto 25 per cent. We further show that the SSA-derived post-seismic deformation signal is smoother and fits better to an exponential model with a decay time of 34 ds. Inversion of the post-seismic deformation suggests an afterslip mechanism with a maximum slip of ∼0.07 m on the shallow, updip portions of the MHT. The short-lived afterslip following the 2019 Mirpur earthquake implies a rapid relocking of the coseismic rupture area and surrounding regions at the MHT. The 2019 Mirpur earthquake and the afterslip together released only 12 per cent of the accumulated strain on the décollement since the 1555 Kashmir earthquake, indicating the potential for future seismic events in the region.

The SSA-based InSAR signal decomposition approach presented in the paper could be further explored for separating tectonic and non-tectonic signals from multitemporal InSAR data sets available from Sentinel-1, ALOS-2/4 and the upcoming NASA-ISRO Synthetic Aperture Radar mission.

ACKNOWLEDGEMENTS

This work was carried out as a part of the Geosciences and Applications Program (GAP) of SAC. MJ is thankful to CSIR, New Delhi for awarding the JRF Fellowship (File No.:09/730(0005)/2019-EMR-I). This work is a part of the PhD thesis of MJ. We thank the editor, Prof. Kosuke Heki, Dr Vineet Gahalaut and an anonymous reviewer for their constructive suggestions.

DATA AVAILABILITY

SAR data from Sentinel-1 are available from Alaska Satellite Facility (https://search.asf.alaska.edu/). SRTM data were retrieved from (https://step.esa.int/auxdata/dem/SRTMGL1/). All software used in this work are openly available: MintPy (https://github.com/insarlab/MintPy); ISCE2 (https://github.com/isce-framework/isce2); GMT (https://github.com/GenericMappingTools). We have made the code used for SSA analysis available at https://github.com/mcmjasir/InSAR_SSA

REFERENCES