Next Article in Journal
Robust Direction-of-Arrival Estimation in the Presence of Outliers and Noise Nonuniformity
Next Article in Special Issue
Coastal Reclamation Embankment Deformation: Dynamic Monitoring and Future Trend Prediction Using Multi-Temporal InSAR Technology in Funing Bay, China
Previous Article in Journal
Comparative Analysis of Machine Learning Models for Tropical Cyclone Intensity Estimation
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Accurate Deformation Retrieval of the 2023 Turkey–Syria Earthquakes Using Multi-Track InSAR Data and a Spatio-Temporal Correlation Analysis with the ICA Method

1
MNR Key Laboratory for Geo-Environmental Monitoring of Great Bay Area & Guangdong Key Laboratory of Urban Informatics, Shenzhen University, Shenzhen 518060, China
2
College of Civil and Transportation Engineering, Shenzhen University, Shenzhen 518060, China
3
Department of Land Surveying and Geo-Informatics, The Hong Kong Polytechnic University, Hung Hom, KLN, Hong Kong 999077, China
4
Guangdong Laboratory of Artificial Intelligence and Digital Economy (SZ), Shenzhen 518000, China
5
School of Architecture & Urban Planning, Shenzhen University, Shenzhen 518060, China
*
Author to whom correspondence should be addressed.
Remote Sens. 2024, 16(17), 3139; https://doi.org/10.3390/rs16173139
Submission received: 12 July 2024 / Revised: 17 August 2024 / Accepted: 21 August 2024 / Published: 26 August 2024

Abstract

:
Multi-track synthetic aperture radar interferometry (InSAR) provides a good approach for the monitoring of long-term multi-dimensional earthquake deformation, including pre-, co-, and post-seismic data. However, the removal of atmospheric errors in both single- and multi-track InSAR data presents significant challenges. In this paper, a method of spatio-temporal correlation analysis using independent component analysis (ICA) is proposed, which can extract multi-track deformation components for the accurate retrieval of earthquake deformation time series. Sentinel-1 data covering the double earthquakes in Turkey and Syria in 2023 are used to demonstrate the effectiveness of the proposed method. The results show that co-seismic displacement in the east–west and up–down directions ranged from −114.7 cm to 82.8 cm and from −87.0 cm to 63.9 cm, respectively. Additionally, the deformation rates during the monitoring period ranged from −137.9 cm/year to 123.3 cm/year in the east–west direction and from −51.8 cm/year to 45.7 cm/year in the up–down direction. A comparative validation experiment was conducted using three GPS stations. Compared with the results of the original MSBAS method, the proposed method provides results that are smoother and closer to those of the GPS data, and the average optimization efficiency is 43.08% higher. The experiments demonstrated that the proposed method could provide accurate two-dimensional deformation time series for studying the pre-, co-, and post-earthquake events of the 2023 Turkey–Syria Earthquakes.

1. Introduction

Interferometric synthetic aperture radar (InSAR) has been widely used to monitor ground deformation relative to geological processes over the past three decades [1,2,3]. In earthquake deformation monitoring, InSAR offers significant advantages over conventional methods, such as leveling and global navigation satellite systems (GNSS), due to its high spatial coverage and fieldwork-free mechanism. Differential InSAR (D-InSAR) is one of the most common methods for monitoring earthquake deformation, functioning by mapping interferogram data according to the difference in SAR images before and after earthquakes, thereby inverting co-seismic deformation patterns [4,5]. Unfortunately, D-InSAR is susceptible to atmospheric and decorrelation factors [6,7], which makes it challenging to accurately capture co-seismic deformation. It also has the limitation of only being able to obtain line-of-sight (LOS) deformation due to the observation geometry, which results in an incomplete description of earthquake deformation characteristics [8].
Many methods to mitigate the impact of atmospheric data on D-InSAR processing have been developed in previous studies, including methods based on external data, such as weather data, Global Positioning Systems [9,10,11], and interferogram modeling [12]. However, these methods are affected by the spatial density of the external data and the low-coherence area in the interferogram, making it difficult to remove atmospheric effects [13]. To obtain multi-dimensional co-seismic deformation, pixel offset tracking (POT) [14] and multiple aperture interferometry (MAI) [15] are typically used. The POT method utilizes cross-correlation between master and slave single-look complex (SLC) images to obtain azimuth and range pixel-offset field maps. The MAI method utilizes split-beam processing to generate backward- and forward-looking SLC images, then generates MAI interferograms to extract deformation along the direction of the track.
Multi-dimensional small baseline subset (MSBAS) can generate a two-dimensional deformation time series in the horizontal and vertical directions through processing small-baseline interferometric datasets from ascending and descending paths [16]. The subset not only overcomes the limitation presented by the fact that single-path D-InSAR can only obtain co-seismic deformation in a single direction, but the earthquake deformation time series results are also conducive to the analysis of pre-, co-, and post-seismic mechanisms [17]. As each deformation map in a time series is solved using multiple interferograms, which can be regarded as a kind of stacking, MSBAS can also overcome decorrelation and alleviate atmospheric delays.
Atmospheric delay is one of the dominant error sources in InSAR, including stratified delay, which is correlated with the topography, and turbulent delay, which can be attributed to random processes in the atmosphere [18]. Stratified delays have a more substantial impact in areas with significant topographic variation, and can be removed using elevation-dependent filtering and external data, such as weather data [19,20]. Turbulent delays are associated with random processes in the atmosphere, including local weather conditions, strong convective effects, variations in local land covers, and different ecosystems [16]. They can be mitigated using methods such as stacking, filtering, stochastic models, and interferogram optimization. However, as a kind of stacking technique, MSBAS can only achieve good atmospheric removal effects in linear deformation scenes, thus, it remains unsatisfactory in many situations [21].
Independent component analysis (ICA) is a computational signal processing method used to separate mixed signal matrices into combinations of several non-Gaussian independent sub-component matrices [22,23]. It can separate deformation and various errors in seismic time series with no prior knowledge, and is suitable for removing atmospheric delays and other error components from interferograms. However, the determination of source types is challenging and usually requires manual identification or external data, making it difficult to remove atmospheric errors when using multi-dimensional monitoring methods. In this study, we propose mitigating atmospheric delays in two-dimensional InSAR data using ICA and spatio-temporal correlation analysis. Taking the double earthquakes in Turkey on 6 February 2023 as an example, two-dimensional deformation components are extracted from ascending and descending Sentinel-1 data and their atmospheric effects are removed using the proposed method.

2. Methods

In this study, we propose a method for removing atmospheric delay errors in multi-track InSAR data and obtaining accurate two-dimensional earthquake deformation time series. Our method incorporates ICA, spatio-temporal correlation analysis, and MSBAS (Figure 1).

2.1. Independent Component Analysis

ICA is a statistical method designed to separate independent components from a mixed signal. It can be used to separate signal components such as surface deformation, atmospheric delay signals, and noise to improve the accuracy of monitoring, thereby facilitating a better understanding of surface changes induced by earthquakes [24], volcanoes [23], landslides [25], and other phenomena. The observation matrix X = x 1 , x 2 , , x n T consists of the n observation vectors x i . The non-Gaussian independent component matrix S = s 1 , s 2 , , s m T consists of the m known independent component vectors s j , which represent m independent components separated with ICA, such as deformation, turbulence, stratified atmosphere, and noise. The mixing matrix A n × m and independent component matrix S comprise the ICA model (Equation (1)), which can be obtained after centering, whitening, and non-Gaussian maximization [26].
X = A S .
For centering, each observation vector x i subtracts its own mean value x i ¯ to become a zero-mean vector. The centered matrix X c is composed of these vectors (Equation (2)):
X c = x 1 x 1 ¯ , x 2 x 2 ¯ , , x n x n ¯ T .
In order to make the components uncorrelated and ensure unit variance, the second step, whitening, is applied, which is usually achieved by performing eigenvalue decomposition on the centered observation. The whitening matrix X c w is obtained using the eigenvector matrix U and the eigenvalue matrix D of X c (Equation (3)).
X c w = U D 1 2 U T X c .
The independent source vectors can be estimated through maximizing the non-Gaussianity. The greater the non-Gaussianity of a signal, the more likely that it is an independent source component [27]. This is because a linear combination of non-Gaussian distributions is more likely to be composed of different independent source signals. Kurtosis and negentropy are common measures of non-Gaussianity [28]. The non-Gaussian maximum objective function is selected to carry out an iterative process. A matrix W , which is the inverse of the mixing matrix A , is calculated after the function converges and the independent component matrix S = W X can be obtained.
ICA can be divided into two categories (Figure 2): spatial ICA (sICA) and temporal ICA (tICA). The objective of sICA is to maximize the spatial independence of the source component. In this mode, each row of S represents a spatial feature of the independent component, while each of column of A represents the temporal characteristics of the corresponding independent component. The objective of tICA is to maximize temporal independence. In this mode, the observation matrix X is a transposed version of sICA’s matrix. Each row of S represents a temporal feature of the independent component, while each of column of A represents the spatial characteristics of the corresponding independent component.

2.2. Spatio-Temporal Correlation Analysis

ICA can be used to separate the monitoring results into independent components and extract deformation sources from the results that have errors. However, the differences in multi-track monitoring results due to different observation geometries and monitoring dates, along with the property of random ordering of independent sources, lead to challenges in correlating deformation-related sources across different tracks. To overcome this issue, a spatio-temporal correlation analysis method is proposed to identify the corresponding sources from multi-track results.
The spatial distributions and temporal features of the deformation components are stable, while those of the atmospheric components are unstable due to the randomness of the atmosphere. Despite their different observation geometries, the deformation components of different paths have similar spatial distribution and temporal characteristics. Therefore, the deformation components can be identified by calculating the spatial and temporal correlations of the spatial pattern and temporal feature vectors, respectively. The spatio-temporal correlation γ i can be calculated using the spatial vectors or temporal vectors of sources, denoted as a and b , respectively, in Equation (4), where i represents the index of the chosen source vectors and j represents the index of the elements in the vectors [29].
γ i = j = 1 N a i j · b i j j = 1 N a i j 2 j = 1 N b i j 2 .
During calculation, the vector with the highest spatial or temporal correlation value is considered a sufficient relevant source. If two sources from different tracks are relevant to each other, they are considered a source pair. Multiple correlation pairs might be found during spatial and temporal analysis. The intersection of the spatial and temporal source pair can be considered the deformation source pair, which is usually a unique source.
In order to test the statistical significance of deformation component extraction, the F-test is employed to determine whether the pair of sources selected is relatively correlated. The F-test [30] value can be expressed as a ratio of the squares of the zero-mean vectors (Equation (5)).
F = j = 1 N a j a j ¯ 2 j = 1 N b j b j ¯ 2 .
The F statistic follows the F distribution with N 1 , N 1 degrees of freedom. The critical value of the F distribution, F N 1 , N 1 , α , is calculated using the given significance level α and the degree of freedom N 1 , N 1 . If F > F N 1 , N 1 , α , the selected source pair is correlated and is considered the deformation source. If F < F N 1 , N 1 , α , the selected source pair is not considered to be correlated. Then, the number of ICA sources should be added, and the ICA and the spatio-temporal correlation analysis should be performed again.

2.3. Multi-Dimensional Small Baseline Subset

The east–west and vertical direction deformation time series V E W and V U D of the common area of the multi-track observed data Φ can be used to generate two-dimensional deformation time series using the MSBAS algorithm [31]. The set of linear equations of MSBAS can be solved in the least-square sense by applying singular value decomposition (SVD). The rank deficiency can be solved with λ -order Tikhonov regularization, where λ = 0,1 , 2 .
A λ L V E W V U D = Φ 0
According to the azimuth angle θ , incident angle ϕ , and time baseline of the interference pair, the east–west and vertical projection relations s E = c o s θ s i n ϕ and s U = c o s ϕ are obtained, constructing the coefficient matrix A. Equation (7), shown below, is the detailed form of Equation (6).
s E a s c s U a s c s E a s c s U a s c 0 0 0 0 0 0 0 0 0 0 s E a s c s U a s c s E a s c s U a s c 0 0 0 0 0 0 0 0 0 0 s E a s c s U a s c s E d s c s U d s c 0 0 0 0 0 0 0 0 0 0 s E d s c s U d s c s E d s c s U d s c 0 0 0 0 0 0 0 0 0 0 s E d s c s U d s c s E d s c s U d s c V 1 E V 1 U V 2 E V 2 U V 3 E V 3 U V 4 E V 4 U V 5 E V 5 U = Φ 1 a s c Φ 2 a s c Φ 3 a s c Φ 1 d s c Φ 2 d s c Φ 3 d s c .
The two-dimensional deformation results for all monitoring dates are calculated and, finally, the two-dimensional deformation time series of the public area with a higher monitoring temporal resolution is obtained based on MSBAS.

3. Study Area and Dataset

3.1. Study Area

On 6 February 2023, a Mw 7.8 earthquake struck Gaziantep Province, Turkey, with a focal depth of 17.5 km. The earthquake’s epicenter was located at 37.226°N, 37.014°E, close to the East Anatolia Fault Zone (EAFZ). After approximately 9 h, a Mw 7.6 aftershock occurred in Kahramanmaras Province, Turkey, with a focal depth of 13.5 km. This second earthquake was situated at 38.011°N, 37.196°E, close to the Cardak Fault (CF) [32]. These earthquakes resulted in over 50,000 fatalities, making them the deadliest earthquake sequence in Turkey in more than a century [33].
The epicenter of the mainshock was located in the EAFZ fault zone at the boundary of the Anatolian and Arabian plates (Figure 3). The Arabian block squeezes the North Anatolian block northwest, and the Anatolian plate slips westward along the NAFZ [34,35]. Despite experiencing eight earthquakes of Mw 7.0 or higher since the 20th century [36], the EAFZ had not witnessed any earthquakes above Mw 7.0 after being compressed by the Arabian plate [37]. This accumulation of significant seismic stress likely contributed to the occurrence of this earthquake [38]. The derived surface deformation fields of the 2023 Turkey–Syria earthquakes have been obtained based on Sentinel-1A ascending and descending data, using POT technology by Dai et al. [39] and based on Sentinel-1 offset tracking and Sentinel-2 optical image correlation carried out by Chen et al. [38]. The multi-scale seismic and space-geodetic observations with multi-fault kinematic inversions and dynamic rupture modeling have been integrated by Jia et al. [40] to unravel the area’s complex rupture history and stress-mediated fault interactions. Kobayashi et al. conducted a Coulomb force analysis, and the results showed that the subsequent 7.6 Mw aftershock was likely triggered by the 7.8 Mw mainshock [41]. In a deformation time series study of the earthquakes in Turkey, the deformation velocity fields of the block–fault structure were constructed, and the main geodynamic processes in the area of the EAFZ were revealed based on stacking-InSAR [42].
The earthquake-affected area has abundant water vapor and is subject to serious atmospheric delay errors. For this research, we conducted a study on the earthquakes in Turkey, employing spatio-temporal correlation analysis and ICA to remove the atmospheric delay and other errors, as well as using MSBAS to obtain two-dimensional deformation time series.

3.2. Dataset

In this study, ascending and descending Sentinel-1 data covering the epicenters of the Mw 7.8 and Mw 7.6 earthquakes were used to obtain two-dimensional deformation time series. The observation period ranged from 4 January 2022 to 16 July 2023. There are 46 scenes each of ascending and descending SLC images, which were obtained during the monitoring period. Correspondingly, 376 scenes of ascending interferograms and 374 scenes of descending interferograms were generated (Table 1).
There are a total of five GPS stations within the Sentinel-1 monitoring area. GPS data from three stations within both ascending and descending data areas were used to verify the MSBAS results after the removal of atmospheric effects (Figure 3).

4. Results and Discussion

The ascending and descending interferograms of co-seismic deformation were generated, and deformation was mainly distributed over the two sides of the fault zones (Table 2, Figure 4a). The region north of the CF and on the northwest side of the EAZF predominantly deformed forward along the LOS direction. Conversely, the area south of the CF and southeast of the EAZF experienced deformation approaching the satellite, along the LOS direction. The deformation range in the descending map spans from −140.0 cm to 109.8 cm, and the maximum positive and negative deformations in the LOS direction are in the southeast EAZF and north of the CF, respectively. The region north of the CF has positive deformation in the LOS direction, and the rest of the areas are primarily dominated by negative deformation (Figure 4b). The deformation range in the ascending map spans from −134.7 cm to 57.6 cm, and the maximum positive and negative spatial distributions are consistent with the descending map. However, compared with the descending data, the spatial distribution of deformation between the CF and EAFZ is slightly different. This may be due to the atmospheric delay and loss of satellite imagery, resulting in a longer temporal baseline.
In this study, MSBAS was used to generate deformation time series of ascending and descending SAR data, following which ICA was performed to remove the atmospheric delay. The linear deformation velocity maps from descending and ascending paths are probably contaminated with atmospheric delay (Figure 5). It is worth noting that the contribution of co-seismic deformation to the accumulated displacement is consistently much larger than all other contributions, which may result in errors in SVD fitting. Therefore, it is necessary to apply a scale factor to the columns of the design matrix related to co-seismic deformation, in order to enhance the stability of the solution [43].
The MSBAS time series cannot accurately represent the deformation characteristics when it is affected by atmospheric delay and other errors. The sICA is used to separate independent sources and extract deformation components, allowing for restoration of the accurate deformation time series. As the order of the ICA results (Figure 6 and Figure 7) does not include correlation information, empirically identifying the corresponding sources from multi-track data is often challenging and leads to inaccurate results. In this study, a spatio-temporal correlation analysis is proposed, in order to match the corresponding independent component source pairs.
The spatial and temporal correlation matrices were generated, which indicate the spatial and temporal correlation values between two different patterns, respectively (Figure 8). The values of the elements in the matrices range from 0 to 1; the higher the value, the higher the correlation. When the matrix element corresponding to a specific row and column is the maximum value, it suggests that the two sources form a spatial or temporal source pair. In the spatial correlation matrix (Figure 8a), three spatial correlation source pairs were matched, including ascending IC source-1 and descending IC source-4, ascending IC source-3 and descending IC source-3, and ascending IC source-4 and descending IC source-1. In the temporal correlation matrix (Figure 8b), two temporal correlation source pairs were matched, including ascending IC source-3 and descending IC source-2, and ascending IC source-4 and descending IC source-1. However, many values in the matrix were close to the maximum value, and are also worthy of attention. It is also worth noting that the monitoring period and arrangement of the ascending and descending orbit time series should be similar and matched. Otherwise, biases can be introduced during the temporal correlation analysis.
As proposed in Section 2.2, the source pair related to both the spatial and temporal analyses can be considered as the deformation source pair. Therefore, ascending IC source-4 and descending IC source-1 comprised the only pair identified as a spatio-temporal correlation source pair.
After the spatial and temporal correlation analysis and the selection of the deformation source pair, we performed an F-test to verify whether the two sources are well correlated. Based on Equation (5), the calculated F-value was 2.4186, which exceeds the critical value at a significance level of 0.95. Consequently, we concluded that ascending IC-4 and descending IC-1 are fully correlated and can be regarded as a deformation source pair. In addition, while the spatial and temporal correlation values of ascending IC source-1 and descending IC source-2 were also high, the F-value did not exceed the critical value.
Compared with the distribution of the fault zone and seismic time series, the corresponding spatial patterns (shown in Figure 6 and Figure 7) of the identified deformation sources were consistent, and the temporal feature vectors were consistent with the earthquake deformation characteristics. Ascending IC source-1 and descending IC source-2 had similar spatial distribution and temporal features; as such, they may be the pre- or post-seismic deformation components.
Stratified atmospheric delay is terrain-dependent. Ascending IC source-5 and descending IC source-4 have spatial distributions related to the terrain elevation and random temporal features. Therefore, these two sources are thought to be attributed to the stratified atmospheric delay. As shown in descending IC source-3, the Erzin region—which is in the middle-west area—has a temporal feature vector with seasonality, which may be caused by Mediterranean water vapor. Thus, descending IC source-3 may be a seasonal turbulent atmospheric component. IC source-5 from the descending data demonstrates that there is a subsidence area in the Antakya area of Turkey. The subsidence is in a farmland area, and the temporal feature vector indicates that the component is not related to the earthquake that occurred in Turkey in 2023. As such, it may be caused by the extraction of groundwater. The remaining IC sources are suspected to be turbulent atmospheric sources.
A two-dimensional deformation time series was obtained after removing atmospheric delay error, and it was used for comparison against the original time series results (Figure 9). Compared with the ascending and descending LOS direction co-seismic deformation (Figure 4), the two-dimensional results with ICA show details not captured in the single LOS results. For example, in the middle-west region of the common area, the subsidence signal is buried in the individual ascending LOS results, but revealed in the multi-track results. Compared with the two-dimensional co-seismic deformation results (Figure 9a,b), the signal caused by the atmospheric delay error is excluded, and the deformation component signal becomes prominent. In the east–west direction (Figure 9a,c), the signal in the middle-west region of the common area and the southern area far away from the fault zone is optimized, revealing the originally subtle subsidence signal and removing the redundant subsidence signal. The displacement in the east–west direction ranges from −114.7 cm to 82.8 cm. In the up–down direction (Figure 9b,d), the deformation of the original result is buried in the signal and challenging to identify; however, after removing the atmospheric error, the deformation signal is dominated. The displacement in the vertical direction ranges from −87.0 cm to 63.9 cm.
In the east–west direction deformation velocity maps (Figure 10a,c), atmospheric signals obscure the deformation signal north of the EAFZ, and it becomes less noticeable. The southern area far away from the fault shows an erroneous subsidence signal due to the influence of a long-term atmospheric trend signal. The deformation rate during the monitoring period ranged from −137.9 cm/year to 123.3 cm/year in the east–west direction. In the up–down direction (Figure 10b,d), the northern area suffered from a stratified atmospheric signal. The removal of atmospheric errors makes the deformation signal more obvious and improves the accuracy of monitoring. The deformation rate during the monitoring period ranged from −51.8 cm/year to 45.7 cm/year in the up–down direction.
The two-dimensional deformation time series in the east–west and up–down directions for both the pre- and post-seismic periods were obtained using MSBAS (Figure 11). In the pre-seismic velocity maps (Figure 11a,b), the deformation in the area north of the EAFZ is mainly upward and eastward, while in the area south of the EAFZ, it is mainly downward and westward. In the post-seismic results, the east–west deformation distribution closely resembles the co-seismic pattern, indicating that similar deformation behavior continued after the earthquake. In the vertical direction (Figure 11d), the uplift signal was distributed between the EAFZ and the CF.
We then validated the effectiveness of the proposed method using the available GPS data. There are three GPS stations in the ascending data area and five stations in the descending area (Figure 3). The three-component GPS data were projected in the LOS direction, and both the GPS and InSAR results were calibrated using the ANTP station—which is relatively far from the epicenter—as a reference (Figure 12). In the ascending data, the results of the proposed method were similar to the original MSBAS results in terms of deformation, but were closer to the GPS results. In the descending data, many errors were removed using the proposed method, causing the results to become smoother and closer to the GPS results than to the original MSBAS results.
In addition, ICA can improve the accuracy of the deformation results from single-orbit data. Based on our calculated standard deviations of MSBAS and ICA relative to the GPS data (Table 3), the average error optimization efficiency increased to 24.14%.
The validation of the removal efficiency evaluation was conducted by comparing the results obtained through the original MSBAS method, the traditional elevation-dependent atmospheric phase trend fitting method, and the proposed method (Figure 13). In both the east–west and vertical directions, the results of the proposed method had less noise and more closely resembled the GPS time profile, outperforming the original MSBAS and traditional fitting methods. The errors in the results generated using the traditional method may be due to the fact that this method only removes the stratified atmosphere related to the terrain and does not remove the influence of the turbulent atmosphere. The standard deviation indicates a difference from the GPS data (Table 4). The average error optimization ratio of the proposed method was 43.08%, which means that the method has a better atmospheric error removal effect.

5. Conclusions

In this study, a spatio-temporal correlation analysis using the ICA method was proposed to identify and extract deformation components from multi-track InSAR data. The 2023 earthquake sequences in Turkey were adopted to test the proposed method. Two-dimensional deformation time series of the 2023 earthquake sequences in Turkey were generated using Sentinel-1 data from ascending and descending tracks, which is beneficial for analyzing pre-, co-, and post-seismic mechanisms. Specifically, in the pre-seismic velocity maps, the deformation north of the EAFZ is mainly upward and eastward, while south of the EAFZ, it is predominantly downward and westward. In the co-seismic deformation maps, the northern region of the CF and the northwestern part of the EAFZ show primarily upward and eastward deformation, whereas the southern areas of CF and EAFZ exhibit downward and westward deformation. In the post-seismic results, the deformation distribution closely resembles the co-seismic pattern in the east-west direction, while in the vertical direction, the uplift signal is concentrated between the EAFZ and the CF. The two-dimensional deformation results obtained with the original, traditional elevation-dependent atmospheric phase trend fitting method and proposed method were compared with GPS data from three stations. The average error optimization ratio of the proposed method was 43.08%, which is 38.49% higher than the traditional surface fitting method and means a better atmospheric error removal effect.
Our study helps in obtaining accurate two-dimensional earthquake deformation time series and providing reliable information for the study of pre-, co-, and post-earthquake deformation mechanisms. Furthermore, this method can be applied to capture accurate multi-dimensional deformation in other scenarios, such as volcanoes, landslides, and glaciers.

Author Contributions

S.W. and Y.L. conceived and designed the experiments; Y.L. performed the experiments, conducted data analysis and validation, and drafted the manuscript; B.Z. contributed to review, editing, and supervision; S.W., S.X. and C.W. provided supervision. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded the National Natural Science Foundation of China grants number 42304014 and 42374018, the Guangdong Basic and Applied Basic Research Foundation grants number 2022A1515110861 and 2024A1515011679, the University Grants Council of the Hong Kong Polytechnic University grant number P0045896, and The Shenzhen Research Institute of the Hong Kong Polytechnic University grant number C2023A022.

Data Availability Statement

The original data presented in the study are openly available in the European Space Agency (ESA) at https://scihub.copernicus.eu accessed on 20 August 2024.

Acknowledgments

We thank the three anonymous reviewers for their constructive suggestions. We acknowledge the European Space Agency (ESA) for freely making available the Sentinel-1A data. Maps were prepared using the QGIS (3.10.4 A Coruña) and MATLAB (R2021a) software.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Rosen, P.A.; Hensley, S.; Joughin, I.R.; Li, F.K.; Madsen, S.N.; Rodriguez, E.; Goldstein, R.M. Synthetic Aperture Radar Interferometry. Proc. IEEE 2000, 88, 333–382. [Google Scholar] [CrossRef]
  2. Wang, C.; Chang, L.; Wang, X.-S.; Zhang, B.; Stein, A. Interferometric Synthetic Aperture Radar Statistical Inference in Deformation Measurement and Geophysical Inversion: A Review. In IEEE Geoscience and Remote Sensing Magazine; IEEE: Piscataway, NJ, USA, 2024. [Google Scholar]
  3. Wang, C.; Wei, M.; Qin, X.; Li, T.; Chen, S.; Zhu, C.; Liu, P.; Chang, L. Three-Dimensional Lookup Table for More Precise SAR Scatterers Positioning in Urban Scenarios. ISPRS J. Photogramm. Remote Sens. 2024, 209, 133–149. [Google Scholar] [CrossRef]
  4. Massonnet, D.; Rossi, M.; Carmona, C.; Adragna, F.; Peltzer, G.; Feigl, K.; Rabaute, T. The Displacement Field of the Landers Earthquake Mapped by Radar Interferometry. Nature 1993, 364, 138–142. [Google Scholar] [CrossRef]
  5. Zhang, B.; Ding, X.; Amelung, F.; Wang, C.; Xu, W.; Zhu, W.; Shimada, M.; Zhang, Q.; Baba, T. Impact of Ionosphere on InSAR Observation and Coseismic Slip Inversion: Improved Slip Model for the 2010 Maule, Chile, Earthquake. Remote Sens. Environ. 2021, 267, 112733. [Google Scholar] [CrossRef]
  6. Agram, P.; Simons, M. A Noise Model for InSAR Time Series. J. Geophys. Res. Solid Earth 2015, 120, 2752–2771. [Google Scholar] [CrossRef]
  7. Zebker, H.A.; Rosen, P.A.; Hensley, S. Atmospheric Effects in Interferometric Synthetic Aperture Radar Surface Deformation and Topographic Maps. J. Geophys. Res. Solid Earth 1997, 102, 7547–7563. [Google Scholar] [CrossRef]
  8. Fuhrmann, T.; Garthwaite, M.C. Resolving Three-Dimensional Surface Motion with InSAR: Constraints from Multi-Geometry Data Fusion. Remote Sens. 2019, 11, 241. [Google Scholar] [CrossRef]
  9. Li, Z.; Muller, J.-P.; Cross, P. Tropospheric Correction Techniques in Repeat-Pass SAR Interferometry; ESA ESRIN: Frascati, Italy, 2003; pp. 1–5. [Google Scholar]
  10. Hopfield, H.S. Tropospheric Effect on Electromagnetically Measured Range: Prediction from Surface Weather Data. Radio Sci. 1971, 6, 357–367. [Google Scholar] [CrossRef]
  11. Bock, Y.; Williams, S. Integrated Satellite Interferometry in Southern California. Eos Trans. Am. Geophys. Union 1997, 78, 293–300. [Google Scholar] [CrossRef]
  12. Chaabane, F.; Avallone, A.; Tupin, F.; Briole, P.; Maître, H. A Multitemporal Method for Correction of Tropospheric Effects in Differential SAR Interferometry: Application to the Gulf of Corinth Earthquake. IEEE Trans. Geosci. Remote Sens. 2007, 45, 1605–1615. [Google Scholar] [CrossRef]
  13. Ding, X.; Li, Z.; Zhu, J.; Feng, G.; Long, J. Atmospheric Effects on InSAR Measurements and Their Mitigation. Sensors 2008, 8, 5426–5448. [Google Scholar] [CrossRef]
  14. Fielding, E.J.; Lundgren, P.R.; Taymaz, T.; Yolsal-Çevikbilen, S.; Owen, S.E. Fault-slip Source Models for the 2011 M 7.1 Van Earthquake in Turkey from SAR Interferometry, Pixel Offset Tracking, GPS, and Seismic Waveform Analysis. Seismol. Res. Lett. 2013, 84, 579–593. [Google Scholar] [CrossRef]
  15. Wang, X.; Liu, G.; Yu, B.; Dai, K.; Zhang, R.; Chen, Q.; Li, Z. 3D Coseismic Deformations and Source Parameters of the 2010 Yushu Earthquake (China) Inferred from DInSAR and Multiple-Aperture InSAR Measurements. Remote Sens. Environ. 2014, 152, 174–189. [Google Scholar] [CrossRef]
  16. Samsonov, S.V.; d’Oreye, N. Multidimensional Small Baseline Subset (MSBAS) for Two-Dimensional Deformation Analysis: Case Study Mexico City. Can. J. Remote Sens. 2017, 43, 318–329. [Google Scholar] [CrossRef]
  17. Yang, C.; Han, B.; Zhao, C.; Du, J.; Zhang, D.; Zhu, S. Co-and Post-Seismic Deformation Mechanisms of the MW 7.3 Iran Earthquake (2017) Revealed by Sentinel-1 InSAR Observations. Remote Sens. 2019, 11, 418. [Google Scholar] [CrossRef]
  18. Kirui, P.K.; Reinosch, E.; Isya, N.; Riedel, B.; Gerke, M. Mitigation of Atmospheric Artefacts in Multi Temporal InSAR: A Review. Remote Sens. Geoinf. Sci. 2021, 89, 251–272. [Google Scholar] [CrossRef]
  19. Xiao, R.; Yu, C.; Li, Z.; He, X. Statistical Assessment Metrics for InSAR Atmospheric Correction: Applications to Generic Atmospheric Correction Online Service for InSAR (GACOS) in Eastern China. Int. J. Appl. Earth Obs. Geoinf. 2021, 96, 102289. [Google Scholar] [CrossRef]
  20. Jolivet, R.; Agram, P.S.; Lin, N.Y.; Simons, M.; Doin, M.; Peltzer, G.; Li, Z. Improving InSAR Geodesy Using Global Atmospheric Models. J. Geophys. Res. Solid Earth 2014, 119, 2324–2341. [Google Scholar] [CrossRef]
  21. Li, Z.; Duan, M.; Cao, Y.; Mu, M.; He, X.; Wei, J. Mitigation of Time-Series InSAR Turbulent Atmospheric Phase Noise: A Review. Geod. Geodyn. 2022, 13, 93–103. [Google Scholar] [CrossRef]
  22. Ebmeier, S.K. Application of Independent Component Analysis to Multitemporal InSAR Data with Volcanic Case Studies: ICA Analysis of InSAR Data. J. Geophys. Res. Solid Earth 2016, 121, 8970–8986. [Google Scholar] [CrossRef]
  23. Maubant, L.; Pathier, E.; Daout, S.; Radiguet, M.; Doin, M.-P.; Kazachkina, E.; Kostoglodov, V.; Cotte, N.; Walpersdorf, A. Independent Component Analysis and Parametric Approach for Source Separation in InSAR Time Series at Regional Scale: Application to the 2017–2018 Slow Slip Event in Guerrero (Mexico). J. Geophys. Res. Solid Earth 2020, 125, e2019JB018187. [Google Scholar] [CrossRef]
  24. Peng, M.; Motagh, M.; Lu, Z.; Xia, Z.; Guo, Z.; Zhao, C.; Liu, Q. Characterization and Prediction of InSAR-Derived Ground Motion with ICA-Assisted LSTM Model. Remote Sens. Environ. 2024, 301, 113923. [Google Scholar] [CrossRef]
  25. Cohen-Waeber, J.; Bürgmann, R.; Chaussard, E.; Giannico, C.; Ferretti, A. Spatiotemporal Patterns of Precipitation-modulated Landslide Deformation from Independent Component Analysis of InSAR Time Series. Geophys. Res. Lett. 2018, 45, 1878–1887. [Google Scholar] [CrossRef]
  26. Hyvärinen, A.; Oja, E. Independent Component Analysis: Algorithms and Applications. Neural Netw. 2000, 13, 411–430. [Google Scholar] [CrossRef]
  27. Comon, P. Independent Component Analysis, A New Concept? Signal Process. 1994, 36, 287–314. [Google Scholar] [CrossRef]
  28. Kumar, M.; Jayanthi, V. Blind Source Separation Using Kurtosis, Negentropy and Maximum Likelihood Functions. Int. J. Speech Technol. 2020, 23, 13–21. [Google Scholar] [CrossRef]
  29. Liang, H.; Zhang, L.; Lu, Z.; Li, X. Nonparametric Estimation of DEM Error in Multitemporal InSAR. IEEE Trans. Geosci. Remote Sens. 2019, 57, 10004–10014. [Google Scholar] [CrossRef]
  30. Kanji, G.K. 100 Statistical Tests; SAGE Publications Ltd.: Thousand Oaks, CA, USA, 2006; pp. 1–256. [Google Scholar]
  31. Samsonov, S.; d’Oreye, N. Multidimensional Time-Series Analysis of Ground Deformation from Multiple InSAR Data Sets Applied to Virunga Volcanic Province. Geophys. J. Int. 2012, 191, 1095–1108. [Google Scholar]
  32. He, L.; Feng, G.; Xu, W.; Wang, Y.; Xiong, Z.; Gao, H.; Liu, X. Coseismic Kinematics of the 2023 Kahramanmaras, Turkey Earthquake Sequence from InSAR and Optical Data. Geophys. Res. Lett. 2023, 50, e2023GL104693. [Google Scholar] [CrossRef]
  33. Zhao, J.-J.; Chen, Q.; Yang, Y.-H.; Xu, Q. Coseismic Faulting Model and Post-Seismic Surface Motion of the 2023 Turkey–Syria Earthquake Doublet Revealed by InSAR and GPS Measurements. Remote Sens. 2023, 15, 3327. [Google Scholar] [CrossRef]
  34. Emre, Ö.; Kondo, H.; Özalp, S.; Elmacı, H. Fault Geometry, Segmentation and Slip Distribution Associated with the 1939 Erzincan Earthquake Rupture along the North Anatolian Fault, Turkey; Geological Society: London, UK, 2021. [Google Scholar]
  35. Reilinger, R.; McClusky, S.; Vernant, P.; Lawrence, S.; Ergintav, S.; Cakmak, R.; Ozener, H.; Kadirov, F.; Guliev, I.; Stepanyan, R. GPS Constraints on Continental Deformation in the Africa-Arabia-Eurasia Continental Collision Zone and Implications for the Dynamics of Plate Interactions. J. Geophys. Res. Solid Earth 2006, 111. [Google Scholar] [CrossRef]
  36. Marza, V.I. On the Death Toll of the 1999 Izmit (Turkey) Major Earthquake; ESC General Assembly Papers; European Seismological Commission: Potsdam, Germany, 2004. [Google Scholar]
  37. Hubert-Ferrari, A.; Lamair, L.; Hage, S.; Schmidt, S.; Çağatay, M.N.; Avşar, U. A 3800 Yr Paleoseismic Record (Lake Hazar Sediments, Eastern Turkey): Implications for the East Anatolian Fault Seismic Cycle. Earth Planet. Sci. Lett. 2020, 538, 116152. [Google Scholar] [CrossRef]
  38. Chen, J.; Zhou, Y. Coseismic Slip Distribution of the 2023 Earthquake Doublet in Turkey and Syria from Joint Inversion of Sentinel-1 and Sentinel-2 Data: An Iterative Modeling Method for Mapping Large Earthquake Deformation. Geophys. J. Int. 2024, 237, 636–648. [Google Scholar] [CrossRef]
  39. Dai, X.; Liu, X.; Liu, R.; Song, M.; Zhu, G.; Chang, X.; Guo, J. Coseismic Slip Distribution and Coulomb Stress Change of the 2023 MW 7.8 Pazarcik and MW 7.5 Elbistan Earthquakes in Turkey. Remote Sens. 2024, 16, 240. [Google Scholar] [CrossRef]
  40. Jia, Z.; Jin, Z.; Marchandon, M.; Ulrich, T.; Gabriel, A.-A.; Fan, W.; Shearer, P.; Zou, X.; Rekoske, J.; Bulut, F. The Complex Dynamics of the 2023 Kahramanmaraş, Turkey, M w 7.8-7.7 Earthquake Doublet. Science 2023, 381, 985–990. [Google Scholar] [CrossRef]
  41. Li, S.; Wang, X.; Tao, T.; Zhu, Y.; Qu, X.; Li, Z.; Huang, J.; Song, S. Source Model of the 2023 Turkey Earthquake Sequence Imaged by Sentinel-1 and GPS Measurements: Implications for Heterogeneous Fault Behavior along the East Anatolian Fault Zone. Remote Sens. 2023, 15, 2618. [Google Scholar] [CrossRef]
  42. Bondur, V.; Chimitdorzhiev, T.; Dmitriev, A. Assessment of Anomalous Geodynamics before the 2023 Mw 7.8 Earthquake in Turkey by Stacking-InSAR Method. Izv. Atmos. Ocean. Phys. 2023, 59, 1001–1008. [Google Scholar] [CrossRef]
  43. Zhang, L.; Hu, J.; Ding, X.; Wen, Y.; Liang, H. Estimation of Coseismic Deformation with Multitemporal Radar Interferometry. IEEE Geosci. Remote Sens. Lett. 2021, 19, 1–5. [Google Scholar] [CrossRef]
Figure 1. Flow chart of the proposed atmospheric delay removal algorithm.
Figure 1. Flow chart of the proposed atmospheric delay removal algorithm.
Remotesensing 16 03139 g001
Figure 2. Schematic diagram of the spatial and temporal ICA on InSAR time series.
Figure 2. Schematic diagram of the spatial and temporal ICA on InSAR time series.
Remotesensing 16 03139 g002
Figure 3. Locations of 2023 Turkey–Syria earthquake sequences. SAR dataset coverage is outlined in green rectangles. Locations of epicenters and GPS stations are marked with circles and triangles, respectively.
Figure 3. Locations of 2023 Turkey–Syria earthquake sequences. SAR dataset coverage is outlined in green rectangles. Locations of epicenters and GPS stations are marked with circles and triangles, respectively.
Remotesensing 16 03139 g003
Figure 4. Co-seismic deformation from (a) descending and (b) ascending SAR data. The epicenters of the Mw 7.8 and Mw 7.6 earthquakes are marked with red diamonds.
Figure 4. Co-seismic deformation from (a) descending and (b) ascending SAR data. The epicenters of the Mw 7.8 and Mw 7.6 earthquakes are marked with red diamonds.
Remotesensing 16 03139 g004
Figure 5. Deformation velocity from (a) descending and (b) ascending SAR data.
Figure 5. Deformation velocity from (a) descending and (b) ascending SAR data.
Remotesensing 16 03139 g005
Figure 6. Independent component analysis results of ascending data. The left column contains the spatial patterns that illustrate the spatial distributions of signals, and the right column shows the corresponding temporal feature vectors, representing the temporal contribution of each source.
Figure 6. Independent component analysis results of ascending data. The left column contains the spatial patterns that illustrate the spatial distributions of signals, and the right column shows the corresponding temporal feature vectors, representing the temporal contribution of each source.
Remotesensing 16 03139 g006
Figure 7. The ICA results of descending data. The left column contains the spatial patterns that illustrate the spatial distributions of signals, and the right column shows the corresponding temporal feature vectors, representing the temporal contribution of each source.
Figure 7. The ICA results of descending data. The left column contains the spatial patterns that illustrate the spatial distributions of signals, and the right column shows the corresponding temporal feature vectors, representing the temporal contribution of each source.
Remotesensing 16 03139 g007
Figure 8. Spatio-temporal analysis matrix: (a) spatial correlation matrix; and (b) temporal correlation matrix.
Figure 8. Spatio-temporal analysis matrix: (a) spatial correlation matrix; and (b) temporal correlation matrix.
Remotesensing 16 03139 g008
Figure 9. Two-dimensional co-seismic deformation maps of (a) original west–east, (b) original up–down, (c) ICA west–east, and (d) ICA up–down directions.
Figure 9. Two-dimensional co-seismic deformation maps of (a) original west–east, (b) original up–down, (c) ICA west–east, and (d) ICA up–down directions.
Remotesensing 16 03139 g009
Figure 10. Two-dimensional deformation velocity: (a) original west–east; (b) original up–down; (c) ICA west–east; and (d) ICA up–down.
Figure 10. Two-dimensional deformation velocity: (a) original west–east; (b) original up–down; (c) ICA west–east; and (d) ICA up–down.
Remotesensing 16 03139 g010
Figure 11. Two-dimensional deformation velocity during pre-seismic period in (a) east–west and (b) up–down directions, and during post-seismic period in (c) east–west and (d) up–down directions.
Figure 11. Two-dimensional deformation velocity during pre-seismic period in (a) east–west and (b) up–down directions, and during post-seismic period in (c) east–west and (d) up–down directions.
Remotesensing 16 03139 g011
Figure 12. Time series deformation in LOS direction from InSAR and GPS datasets.
Figure 12. Time series deformation in LOS direction from InSAR and GPS datasets.
Remotesensing 16 03139 g012
Figure 13. Time series deformation in horizontal and vertical directions from InSAR, traditional elevation-dependent atmospheric phase trend fitting method, proposed method, and GPS datasets.
Figure 13. Time series deformation in horizontal and vertical directions from InSAR, traditional elevation-dependent atmospheric phase trend fitting method, proposed method, and GPS datasets.
Remotesensing 16 03139 g013
Table 1. The parameters of the SAR dataset used.
Table 1. The parameters of the SAR dataset used.
PathPass DirectionFrameDateNo. of SLCsNo. of InterferogramsIncident AngleAzimuth Angle
116Ascending11920220104–202307104637633.82°−13.23°
21Descending465,47120220110–202307164637422.83°−166.72°
Table 2. Information of co-seismic image pairs.
Table 2. Information of co-seismic image pairs.
PathPass DirectionMasterSlavePerpendicular
Baseline (m)
Temporal
Baseline (Days)
116Ascending2023020420230229113.5424
21Descending2023012920230210111.1912
Table 3. Standard deviation of the LOS displacement between GPS and InSAR.
Table 3. Standard deviation of the LOS displacement between GPS and InSAR.
PathStationOriginal MSBAS STD (cm)Proposed Method STD (cm)Ratio
AscADY10.00220.00207.10%
AscANTPRefRefRef
AscKAHR0.01770.011236.86%
DscADY10.01370.010424.46%
DscANTPRefRefRef
DscFEEK0.01270.007739.20%
DscKAHR0.02740.022816.50%
DscONIY0.01750.013920.70%
Table 4. Standard deviation of the two-dimensional displacement between GPS and InSAR.
Table 4. Standard deviation of the two-dimensional displacement between GPS and InSAR.
StationDirectionOriginal MSBAS STD (cm)Traditional Method STD (cm)RatioProposed Method STD (cm)Ratio
ADY1East–West0.00350.003110.94%0.001266.83%
ANTPEast–WestRefRefRefRefRef
KAHREast–West0.02240.02245.20%0.017223.29%
ADY1Up–Down0.00200.00194.95%0.001620.94%
ANTPUp–DownRefRefRefRefRef
KAHRUp–Down0.00360.00351.29%0.001361.25%
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Liu, Y.; Wu, S.; Zhang, B.; Xiong, S.; Wang, C. Accurate Deformation Retrieval of the 2023 Turkey–Syria Earthquakes Using Multi-Track InSAR Data and a Spatio-Temporal Correlation Analysis with the ICA Method. Remote Sens. 2024, 16, 3139. https://doi.org/10.3390/rs16173139

AMA Style

Liu Y, Wu S, Zhang B, Xiong S, Wang C. Accurate Deformation Retrieval of the 2023 Turkey–Syria Earthquakes Using Multi-Track InSAR Data and a Spatio-Temporal Correlation Analysis with the ICA Method. Remote Sensing. 2024; 16(17):3139. https://doi.org/10.3390/rs16173139

Chicago/Turabian Style

Liu, Yuhao, Songbo Wu, Bochen Zhang, Siting Xiong, and Chisheng Wang. 2024. "Accurate Deformation Retrieval of the 2023 Turkey–Syria Earthquakes Using Multi-Track InSAR Data and a Spatio-Temporal Correlation Analysis with the ICA Method" Remote Sensing 16, no. 17: 3139. https://doi.org/10.3390/rs16173139

APA Style

Liu, Y., Wu, S., Zhang, B., Xiong, S., & Wang, C. (2024). Accurate Deformation Retrieval of the 2023 Turkey–Syria Earthquakes Using Multi-Track InSAR Data and a Spatio-Temporal Correlation Analysis with the ICA Method. Remote Sensing, 16(17), 3139. https://doi.org/10.3390/rs16173139

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop