Incorporating Independent Component Analysis and multi-temporal SAR techniques to retrieve rapid postseismic deformation

: This study investigates the ongoing postseismic deformation induced by two moderate mainshocks of Mw 6.1 and Mw 6.0, 2017 Hojedk earthquake in Southern Iran. Available Sentinel-1 TOPS C-band Synthetic Aperture Radar (SAR) images over about one year after the earthquakes are used to analyze the postseismic activities. An adaptive method incorporating Independent Component Analysis (ICA) and multi-temporal Small BAseline Subset (SBAS) Interferometric SAR (InSAR) techniques is proposed and implemented to recover the rapid deformation. This method is applied to the series of interferograms generated in a fully constructed SBAS network to retrieve the postseismic deformation signal. ICA algorithm uses a linear transformation to decompose the input mixed signal to its source components, which are non-Gaussian and mutually independent. This analysis allows extracting the low rate postseismic deformation signal from a mixture of interferometric phase components. The independent sources recovered from the multi-temporal InSAR dataset are then analyzed using a group clustering test aiming to identify and enhance the undescribed deformation signal. Analysis of the processed interferograms indicates a promising performance of the proposed method in determining tectonic deformation. The proposed method works well, mainly when the tectonic signal is dominated by the undesired signals, including atmosphere or orbital/unwrapping noise that counts as temporally uncorrelated components. In contrast to the standard SBAS time series method, the ICA-based time series analysis estimates the cumulative deformation with no prior assumption about elevation dependence of the interferometric phase or temporal nature of the tectonic signal. Application of the method to 433 Sentinel-1 pairs within the dataset reports two distinct deformation patches corresponding to the postseismic deformation. Besides the performance of the ICA-based analysis, the proposed method automatically detects rapid or low rate tectonic processes in unfavorable conditions


INTRODUCTION
The study of the crustal deformation induced by earthquakes provides insight into the rheology of the Earth's crust and active tectonic processes on faults.Understanding the transient response of the lithosphere to the stress changes caused by large or moderate earthquakes could also be useful for the earthquakes prediction and seismic hazard analysis.
Synthetic Aperture Radar (SAR) interferometry is a popular technique to remotely detect centimetric to millimetric deformation over extensive coverage Earth's surface associated with various phenomena including earthquakes, volcanic activity and landslide movements (Bagnardi and Hooper, 2018;Hooper et al., 2004;Huang et al., 2017;Moreno et al., 2018;Wang et al., 2017;Wasowski et al., 2011).One of the main problems facing SAR interferometry is the presence of heterogeneous propagation delays due to atmospheric fluctuation between the two acquisition times.Atmospheric phase delay degrades the quality of SAR interferometry to detect centimetric or millimetric deformation.The temporally uncorrelated components of the atmosphere obscure the weak geophysical signals.It is the issue facing the estimation of deformation in interseismic or postseismic processes (Tarayre and Massonnet, 1996).
There are several methods to mitigate the atmospheric signal, involving stacking a set of short-baseline and long-time span interferograms (Walters et al., 2013), spatiotemporal filtering including applying high-pass and low-pass filtering in time and space, respectively (Ferretti et al., 2001), statistical method to identify and removing the correlated topographic component at single scale (Cavalié et al., 2007), multi scales (Lin et al., 2010), direct estimating of atmospheric delays using independent prediction of water vapor from external sources such as GPS (Onn and Zebker, 2006) or multispectral satellite data, including Moderate Resolution Imaging Spector radiometer (MODIS) and Medium Resolution Imaging Spectrometer (MERIS) (Li et al., 2009).
However, multi-temporal approaches, e.g., stacking or spatiotemporal filtering based methods, are no longer useful where deformation period is much shorter than the measurement interval.It's due to the dominance effects of temporally uncorrelated components, mostly caused by atmosphere error.On the other hand, using the statistical method based on removing the elevation-dependent component from the interferometric phase is not effective where the tectonic signal correlates with topography.Performing direct methods is also limited by the availability of those external sources that predict water vapor content.An adaptive method is proposed here to retrieve the rapid deformation signal.The method includes three steps: first decomposing the multi-temporal signal to its independent sources, second splitting the extracted independent sources into two separate groups, third identify the source containing deformation based on the cluster analysis.The reconstructed interferograms are then generated based on the selected sources resulted from clustering analysis.Finally, a time series map is produced using a modified Small BAseline Subset (SBAS) time series analysis.

Study area and available data
On December 1, 2017, an Mw6.1-earthquake struck a region between Central Iran and Lut block in Southeastern Iran.The mainshock followed by two moderate aftershocks (Mw 5 and Mw 6) occurred on December 12, positioned in close spatial proximity to each other.The area has been hit continuously by the ten tremendous events (Mw > 5) over the eight months.The epicenters of the sequence events were located 25 km Northeast of the town of Hojedk, with a population of 3,000 people in Kerman Province of Iran.Here is this study, available SAR data acquired from Sentinel-1 C-band are collected to examine early near field postseismic deformation ~1 year after the earthquakes.

Independent Component Analysis
Independent Component Analysis (ICA) is an algorithm that aims to decompose a randomly mixed signal to a linear combination of individual components.The decomposition is in the was that the resultant signals are mutually independent, and each of them has a non-Gaussian probability distribution (Comon, 1994;Hyvärinen and Oja, 1997).ICA, as a Blind Source Separation (BSS) algorithm, is widely used to extract blind sources from the mixed-signal (Amato et al., 2008;Barnie and Oppenheimer, 2015;Liu et al., 2015) Suppose there are series of sources, s(t), t=1, …, N where each of them has a fixed probability distribution and the signals were observed indirectly using M individual sensors through M observations, x(t), t=1, …, M. The general objective is to model the mixing process that generates sources through a linear basis transformation as follow: (Hyvärinen and Oja, 1997).
Where x(t) and s(t) are observed and unknown source signal, respectively.The matrix A is an unknown mixing matrix consisting of the coefficients describing the relation between individual source and particular mixed observation.For simplicity let's assume A is a square matrix then independent components (ICs),, are estimated by its inverse,  as: Mixing matrix transforms observed feature vectors to their independent sources such that transformation maximizes the nongaussianity of the resultant ICs.
Before applying ICA, two pre-processing steps, including centering, and whitening, are performed on the observations.The first step, centering, allows simplifying ICA by subtracting the signal's mean value,  = (), from the original signal.The whitening step transforms the mixed signal to uncorrelated variables with variance equal to 1 as satisfies the following equation.

𝐸(𝑥 𝜔 𝑥 𝜔
where   is whitened vector achieved based on Eigen Value Decomposition (EVD) of x by : where V is eigenvectors matrix of (  ) and D is a diagonal matrix consists of eigenvalues.The whitening step minimizes the complexity of ICA by reducing the number of unknown parameters through converting mixing matrix A to an orthogonal matrix   (Hyvärinen and Oja, 2000).

Clustering analysis
Identifying real source (s) related to deformation signal among other resultants ICs is difficult, especially in a case of blind or poorly constrained deformation signal.One can identify the desired source, e.g., tectonic signal in this study, by visual inspection of independent signals, s, or mixing matrix A. Still, the visual inspection requires prior knowledge about the temporal and spatial distribution of deformation.On the other hand, due to the iterative nature of ICA and using a random starting point in estimating each row of the unmixing matrix, , the order of resultant ICs is different in each iteration.It's then difficult to detect the desired component automatically.
Automatic detection of deformation source is of interest in Interferometric SAR (InSAR) studies as hundreds to thousands of pairs are produced to retrieve the tectonic signal.
Automatically solving the source identification can be done by testing the statistical significance of the retrieved ICs sources (Hyvärinen and Ramkumar, 2013).The testing process is done over multiple ICs groups that are produced by randomizing the input data and estimating ICs for each random set (Hyvärinen and Ramkumar, 2013).Testing a group of ICs provides suggestions for identifying blind or poorly constrained signal.
In this study, a modified method inspired by that proposed by Ebmeier (2016), is introduced.The proposed method, (Ebmeier, 2016), was based on using the ICA test that was first developed by (Hyvärinen and Ramkumar, 2013) to examine consistency either in inter-subject or inter-session in neuroscientific research.
The method was based on the following steps.First, it divides a series of interferometric pairs into two groups.The division is done either randomly or systematically, depending on the nature of the source.Second, two sets of IC components are constructed by performing ICA analysis individually on each group.The third cluster of similar components is identified, and finally, in the last step, a set of interferograms is constructed based on the identified cluster corresponds to the deformation signal.
One of the issues associated with the study area is related to the significant contribution of atmospheric noise, mostly the turbulent mixing component.Therefore, the representation of the consistent deformation signal may be affected by the atmospheric noise.For instance, if the turbulent mixing phase delay exists in one image, it may act as a dominant signal in all the single-master interferometric pairs produced by that image.
Therefore, the following procedure is pursued in this study.First, a fully connected SBAS network is constructed.SBAS network is then divided into multiple single-master interferometric subsets.Each of these subsets is further split into two groups.As the rapid postseismic deformation being studied here persists over the whole period of the InSAR observations.Therefore, the splitting step is done such that the interferograms spanning subsequent and separated periods are divided into two different groups.The analysis proceeds with running the ICA separately over each group and identifying clusters with consistent sources.
Two clusters with the highest rate of the similarity are selected for each single-master interferometric subset.Here, one can assume, spatially correlated components related to the atmospheric phase delay, of the master image, might cover all the single-master interferometric pairs and act as the undesired dominant component.The process continues for each singlemaster interferometric subset yields to collect all the clusters representing the authoritative sources, either desired or unwanted components under the study.The resultant ICs components are separated again in two groups, and the final ICA test is done over the groups of ICs to distinguish the real tectonic signal from the undesired dominant signal.
In the last steps, all the InSAR subsets are reconstructed based on the recovered clusters to form modified SBAS interferograms.
Reconstruction is done by calculating the outer product of the relevant row of the unmixing matrix and the identified source (s) obtained from the clustering analysis.The reconstructed SBAS interferograms are used through SBAS analysis to produce the final time series map.

SBAS analysis
The semi-arid condition of the study area increases the probability of producing high-coherence interferograms that allows constructing a fully connected SBAS network.Figure 1 shows the SBAS network that consists of connecting lines representing each interferometric pair.Connections in the SBAS network are differently colored based on the average coherence value of the corresponding interferometric pairs.Figure 1 shows most of the interferograms involved in the time series processing have high coherence with values larger than 0.3, which is related to the semi-arid condition of the study area.
As discussed in session 2.3, SBAS interferograms are first divided into a number of subsets, each of which has a single master date and it is then used as an input for ICA decomposition analysis.After reconstructing SBAS interferograms, SBAS time series analysis is done over the reproduced interferograms.Time series analysis was individually applied for each PS point.It means every point could have its own SBAS network, and it allows increasing the number of observations as high as possible during the time series analysis of each point.It also helps to accurate the time series analysis by using a dense SBAS network as increasing the number of observations could reduce uncorrelated atmospheric effects.

RESULTS AND DISCUSSION
Figure 2 illustrates an example showing the result of the clustering analysis applied to one of the InSAR subsets.Fig. 2a shows components that are separately estimated by performing spatial ICA over two input groups.Based on the visual inspection, IC1, IC7 resulted from group 1 are similar to IC2, IC4 from the second group (Fig. 2a).
Fig. 2b indicates two clusters showing the most similar independent sources retrieved from each group.Fig. 2c represents the original interferograms, including 28 interferometric pairs versus the reconstructed interferograms.A dominant pattern is related to the turbulent mixing noise of the master image, which covers all the original interferograms, the pattern hides the postseismic signal.But when postseismic deformation is accumulated over time, the geophysical signal gets large enough to overcome atmospheric noise.For instance, despite the large atmospheric error, a distinguishable postseismic signal is obvious from the interferogram associated with date 28 (Date number 28 in Fig. 2c).The reconstructed interferograms, illustrated in Fig. 2c, are estimated based on the identified cluster correspond to the tectonic signal.Comparing the original and reconstructed interferograms indicates that applying ICA cluster analysis can monitor the evolution of the cumulative postseismic deformation early after the earthquake.
The resultant SBAS time series maps are represented in Fig. 3.The figure illustrates the comparison of two sets of time series maps.The first subset is derived from original interferograms (Fig. 3a) and the second one estimated from the constructed interferograms (Fig. 3b).
Figure 3. a: Time series maps resulted from the original interferograms correspond to 30 dates that span 1 year after the earthquake.b: Time series maps calculated by reconstructed interferograms.
Fig. 3a slightly shows the accumulation of the postseismic signal even for the early dates after the earthquake, which is not possible to identify based on the time series maps obtained from singlemaster interferograms (shown in Fig. 2c).It's a benefit of using the SBAS technique that could roughly overcome the atmospheric effect.The stacking of a large number of interferograms could reduce the impact of noises; those are temporally uncorrelated.However, some undesired patterns still exist in the cumulated postseismic map (Fig. 4b) that might be misinterpreted as a tectonic signal.As the presence of any undesired signal could affect the seismic modeling analysis, it's better to clean the postseismic deformation map from any undesired effects before performing the postseismic modeling process.Fig. 4d shows the relationship between topography and the cumulative displacement map.It indicates cumulative deformation derived from original interferograms correlates well (with a correlation coefficient of 0.69) with topography that is due to the nature of tectonic signal in case of reverse faulting event.However, topographic dependence is decreased for the resultant cumulative displacement (with a correlation coefficient of 0.002) obtained from the reconstructed interferograms.It proves the ability of ICA based analysis in extracting the valid tectonic signal among the other involved contributions (e.g., elevation-dependent component), considering no prior assumption about elevation dependence of the interferometric phase.
To evaluate the resultant cumulative displacement maps , time series plots are drawn over two points selected in the near-field area of the first (occurred on December 1) and second earthquake (occurred on December 12).P1 and P2 mark the location of the points in Fig. 3b.
The time series plots are depicted in different colors of red and blue corresponding to the time series plot estimated based on original and reconstructed interferograms, respectively (Fig. 5).
The geophysical signal is assumed to be temporally correlated; it means it varies smoothly over time.Comparing the time series plots (Fig. 5a,b) confirm increasing the signal to noise ratio after applying the proposed method as the blue graph have smoother variation over time.However, red lines showing jumps in some dates that indicate the existence of the atmospheric noise in the original interferograms.

CONCLUSION
This study uses ICA analysis to automatically explore the rapid postseismic signal by using multi-temporal InSAR data.An adaptive method is proposed here to retrieve the undescribed deformation source by maximizing the independence of the source component.Analysis of 433 sentinel-1 interferograms using incorporation of ICA and SBAS analysis is able to retrieve two distinct deformation patches corresponding to the postseismic deformation.Placing the coseismic pattern over the postseismic deformation map indicates two patches retrieved as the postseismic deformation surround the coseismic patches associated with the first (December 1) and second (December 12) mainshocks.
The results from SBAS time series analysis show that applying ICA analysis leads to better retrieval of the consistent signal even in the presence of atmospheric noise.ICA based analysis is also able to reduce the elevation-dependent component as topographic dependence of the cumulative displacement is decreased by considering the reconstructed interferograms as input for the time series analysis.The RUPD error of the improved interferograms reconstructed based on cluster analysis decreased by 62 % compared to the original interferograms.

Figure 1 .
Figure 1.SBAS network representing interferometric pair colored based on the average coherence value of each interferogram.Horizontal and

Figure 2 .
Figure 2. Illustration of ICA cluster analysis.a: represents the ICs estimated from the first and second group of the input data.b: Identified clusters of independent components obtained from ICA clustering analysis.c: Comparison of two subsets first includes all the original interferograms produced by a single master and second the subset made of the reconstructed interferograms.Reconstructed interferograms are reproduced based on the optimum cluster that represents best the postseismic signal.

Figure 4
Figure 4. a: Topographic map.b: Cumulative displacement map as a result of time series analysis performed on the original dataset.c: Cumulative displacement map applied to the reconstructed interferograms obtained from ICA based analysis.The black isocontours are superimposed on the cumulative deformation maps that represent coseismic signal resolved from two events of 1 and 12 December 2017.d,e: Scatter plots showing strong and week correlation between topography and cumulative displacements represented in (b) and (c), respectively.The time series maps estimated based on the reconstructed interferograms, represented in Fig.3b, clearly shows the

Figure 5 .
Figure 5.Time series plots illustrating the time-dependent evolution of postseismic deformation on points P1 (a) and P2 (b), which are depicted in Fig. 3b.The last evaluation is done by calculating the RMS of the unmodelled phase delay (RUPD) for each interferogram.It allows defining the effectiveness and reliability of the ICA based analysis proposed here.RUPD obtained for each interferogram by calculating the RMS of the residual, the difference between observed and estimated interferometric phase (in SBAS series model).Fig. 6 shows RUPD estimated for the original and reconstructed interferograms in which red bars correspond to the RUPD error of the time series analysis obtained based on the reconstructed interferogram and blue is related to the original interferogram.The improvement after reducing atmospheric noise can be seen in this figure, Fig. 6. as most of the RUPD errors (showing in blue bars) are decreased by 62 % when time series runs over reconstructed interferograms (red bars in Fig. 6).

Figure 6 .
Figure 6.Representation of RUPD error calculated for each interferogram.Blue corresponds to the time series analysis applied on the original interferograms, red to the time series analysis performed on the reconstructed interferograms.