Three-Dimensional Displacement Fields from InSAR through Tikhonov Regularization and Least-Squares Variance Component Estimation

When attempting to integrate the various range and azimuth displacements with different radar wavelengths (e.g., Sentinel-1, ALOS-2) for 3D displacement vectors' extraction, heterogeneous observations with different precisions come across. For estimating the unknown parameters through heterogeneous data, choosing the proper weights is necessary (Koch and Kusche 2002; Kusche 2003b). The proper weight matrix could be achieved by applying the LS-VCE method in which different noise components are estimated properly. In applying the LS-VCE method, the homogenous observations on the basis of their statistical properties are categorized in one set with their own unknown variance components (Teunissen and Amiri-Simkooei 2008). Here, we consider the stochastic part of the linearized Gauss-Markov model of observations as the following equation:

where $C_k,k=1,\dots ,p$, are some $n\times n$ symmetrical primary known covariance matrices described by Bamler and Eineder (2005) and Ansari et al. (2016); and ${\sigma }_k$ = corresponding unknown variance components that should be estimated through LS-VCE for $p$ categories of data. The symmetrical matrix $C_0$ is assumed to be the known part (if it exists) of the stochastic model. The $p$ vector of unknown variance components are calculated as follows:

where the entries of the $p\times p$ normal matrix $N$ and the $p\times 1$ vector $l$ are presented as

and

where $j=1,\dots ,p,\widehat{e}=P_A^{\perp }l$ are the estimated residuals; and $P_A^{\perp }$ is the orthogonal projector presented as

Here, $N^{-1}$ naturally becomes the variance–covariance matrix of the variance components (i.e., $C_{\widehat{\sigma }}=N^{-1}$). The covariance matrix $C_y$ and the projector $P_A^{\perp }$ are both functions of the unknown variance components. Therefore, the process of estimating unknowns should be iterated until the estimated values stop changing by further iterations (Amiri-Simkooei 2013). After estimating the proper variance components (i.e., $\widehat{\sigma }$) and consequently the covariance matrix of observations (i.e., $C_y$) through LS-VCE, these values are applied in WLS [Eq. (4)] for constructing the 3D displacement vectors.

Occurrence of negative variance components is inevitable for various reasons such as an insufficient number of observations in the functional model and an improper structure of the stochastic model (Sjöberg 1984). To ensure that the estimated variances are nonnegative, one can impose the nonnegativity constraints $\left(\sigma \ge 0\right)$ to the linear (co)variance component model (Shaw and Geyer 1997). Teunissen (1988) suggested the reparameterization of the model in a sense that the nonnegativity of the variance components is automatically ensured, although this approach may turn a linear LS-VCE problem into a nonlinear LS-VCE problem. Moghtased-Azar et al. (2014) applied the positive-valued functions (PVFs) concept for unknown variance components in a stochastic model and the reparameterized approach on the REML with no effect on the unbiasedness of the scheme. Amiri-Simkooei (2016) presented the nonnegative least-squares variance component estimation (NNLS-VCE) based on standard theories on NNLS and LS-VCE.

The DoFs of the functional model in Eq. (1) are$d{f}_1=n-u$, where $n$ and $u$ are the numbers of InSAR displacements and numbers of the unknown components of 3D displacements (i.e., $u=3$), respectively. In the stochastic models [Eq. (5)], the DoFs are $d{f}_2=[d{f}_1(d{f}_1+1)]/2-p$, where $p$ is the number of unknown variances (Teunissen and Amiri-Simkooei 2008). Usually, the number of InSAR observations in relation to the three unknowns of 3D displacements for each pixel is not enough to apply VCE. Therefore, observations of some neighborhood cells are taken into account to increase the redundancy of the stochastic model. A series of experiments were carried out by Liu et al. (2018) with window sizes ranging from 3 × 3 to 25 × 25 pixels for accessing the performance of the strain model VCE. They showed that the root-mean-square errors (RMSEs) of all three components decrease with the increase in the window size. However, selecting the optimum size of the window is a trade-off between the accuracies of 3D displacement estimations and the burden of the computation. In this context, the blocks of 3 × 3, 5 × 5, and 7 × 7 cells, which are equivalent to 9, 25, and 49 neighboring cells, respectively, are assessed in the simulated data sets. Consequently, a block of 3 × 3 cells is chosen as the optimum size of the moving window due to the RMSEs comparison, convergence threshold of repetition procedure of LS-VCE, and computational burden. The DoFs of the equations system before applying the observations of the adjacent cells will be called the primary DoFs. A schematic form of these adjacent pixels with their three LOS displacements is shown in Fig. 1.

In general, retrieving the 3D displacement vectors from range displacements is an ill-posed inverse problem, where any small change in measurement can lead to a great change in the estimated displacements, especially when the functional primary DoF is zero. Usually, it is possible to stabilize the inversion process truncating the small singular values or imposing additional constraints generally referred to as regularization (Tikhonov 1963; Tikhonov et al. 1977; Schaffrin 1980; Xu 1992; Aster et al. 2013). Here, the TR solution is applied due to its efficiency and low DoFs of the problem in increasing the stability of normal equations. The regularization of the inverse problem introduces a trade-off between the norms of the regularized solutions and the observation residuals. This solution is based on minimizing the following objective function (Hansen 1990):

where $\alpha$ = positive regularization parameter; and $\parallel .{\parallel }_{L_2}$ is the Euclidean $L_2$ norm. The $\alpha$ parameter can be estimated through the L-curve criterion (Hansen 1992, 2007), generalized cross validation (Wahba 1976), and VCE (Arsenin and Krianev 1992). The regularized solution of Eq. (10) and its variance–covariance matrix $(x_{\text{reg}}\text{ and }{C}_{x_{\text{reg}}})$ would yield

where $P$ and $I$ are the weight and identity matrices, respectively. Adding the regularization parameter biases the results of the unknown parameters (Hoerl and Kennard 1970; Xu 1992; Xu et al. 2006; Shen et al. 2012) and the unbiased retrieved 3D displacements are obtained as follows:

provided that the inverse of the estimated covariance matrix of observations [i.e., Eq. (5)] through the LS-VCE method is applied as a weight matrix $(P)$ in the regularized solution. To be simple, the term of regularized least-squares variance component estimation (RLS-VCE) is applied subsequently to present the use of weights of the LS-VCE method in retrieving the 3D displacement vectors through Eq. (12).

In this sample, three 500 × 500 pixel images were synthetized along the range direction (LOS) of ascending and descending orbits of Sentinel-1 and descending orbit of ALOS-2 based on the local incidence and azimuth angles of the missions for the Illapel region. The heading angle values varied from 188.7°, 194.8°, and 343.8° to 190.9°, 195.8°, and 344.7° for descending ALOS-2 and descending and ascending Sentinel-1, respectively. The values of incidence angle for these three missions varied from 38.2°, 31.7°, and 37.8° to 49.3°, 43.6°, and 45.7°, respectively. In this context each pixel was observed with three independent geometries where the primary functional DoF is zero and a block of $3\times 3$ cells including nine pixels around the centered pixel was implemented in the stochastic model for estimating the variance components of observations through LS-VCE. The LOS displacements were applied to estimate the 3D displacements through CM Eq. (4) and RLS-VCE method Eq. (12).

In this sample, the number of unknown parameters was five—three for displacement components and two for variance components of Sentinel-1 and ALOS-2—while the number of observations for LS-VCE was 27, that is, three independent observables for nine cells. The results of CM and RLS-VCE in the three directions of eastern, northern, and vertical directions are illustrated in Figs. 2(b and c), respectively. The results indicate that the RLS-VCE method retrieves the eastern, northern, and vertical components better than CM. The overall root-mean-square errors (overall RMSEs for all 3D components simultaneously) of conventional and unbiased RLS-VCE methods with respect to known 3D values were 11.15% and 2.98 cm, respectively. This indicates approximately 73% improvement in the accuracy of estimating 3D displacement vectors when the RLS-VCE method is implemented instead of CM. The RMSEs of the eastern, northern, and vertical components were 0.49, 18.58, and 5.27 cm for CM and 0.24, 4.87, and 1.67 cm for RLS-VCE, respectively; that is, improvements of about 51%, 74%, and 68% for eastern, northern, and vertical directions when using RLS-VCE instead of CM; the results are summarized in Table 1. The expected fact that the northern displacement component is extraordinarily more sensitive to the errors of the LOS measurements of InSAR than the other two components is clearly visible in the northern component in Fig. 2. Compared to the eastern and vertical components, the northern component of the displacement field is not well resolved when only LOS observations are used.

To evaluate the efficacy of the LS-VCE method, five 500 × 500 pixel images of range and azimuth displacements in both ascending and descending orbits of Sentinel-1 and descending range of ALOS-2 were synthetized according to Eq. (13). Therefore, the primary DoFs of the functional model were 2. The incidence and azimuth angles of Sentinel-1 and ALOS-2 were synthetized approximately similar to the values of the Illapel region. These synthetized displacements were applied in estimating the 3D displacements through the CM and LS-VCE methods.

The results of CM and LS-VCE in three directions are shown in Figs. 3(b and c), respectively. The overall RMSE of the CM and LS-VCE methods with respect to known 3D values were 8.50 and 5.20 cm, respectively. Therefore, the use of the LS-VCE method that estimates the best value of observational weights improves the accuracy of retrieving 3D displacement vectors by approximately 39% compared to CM. The componential RMSEs of eastern, northern, and vertical components were 0.4, 13.1, and 6.7 cm for CM and 0.3, 8.5, and 2.9 cm for RLS-VCE, respectively; that is, improvements of about 25%, 35%, and 57% for eastern, northern, and vertical directions when using LS-VCE instead of CM; the results are summarized in Table 1. The fact that adding the azimuth displacements drastically stabilizes the solution and leads to a better estimate of the northern component of displacement is clearly highlighted in Fig. 3 (in comparison with Fig. 2), as well as in the condition number (as a measure of the well-posedness) of the normal matrix, which decreased from 5,000 to 5. Consequently, the problem will be solved without any regularization process.

In this section, the results of our methods for the 2015 $M_{\text{W}}$ 8.3 Illapel (Chile) earthquake are examined and discussed. This event occurred on September 16, 2015 (22:54:33 GMT), 46 km offshore from Illapel, Chile, at the depth of 25.0 km. The InSAR data for this earthquake were obtained from a stack of three interferogram pairs in the following sequence: (1) three adjacent descending frames of terrain observation by progressive scan (TOPS) Sentinel-1, (2) two adjacent ascending frames of TOPS Sentinel-1, and (3) a pair of descending frames of phased array type L-band synthetic aperture radar (PALSAR-2) ALOS-2 (Table 2).

The repeat-pass interferometry technique was applied to process the SAR data. The small time span of the acquiring pairs, low levels of precipitation, sparse vegetation, and magnitude of the deformation signal make a significant part of the interferometric phase coherent in the coseismic interferograms, except for the mountainous area to the east.

The two adjacent ascending frames and three adjacent descending frames of Sentinel-1 were processed using GAMMA software. The frames were first concatenated to produce long frames in ascending and descending orbits. Following this, wrapped and unwrapped interferograms of ascending and descending frames of Sentinel-1 were constructed using a 30-m resolution Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) for topography removal and geocoding. Moreover, the azimuth components of displacements of ascending and descending directions were deduced through burst overlap interferometry through double difference between backward- and forward-looking interferograms within the region of burst overlap (Grandin et al. 2016). The PALSAR-2 ALOS-2 interferometric processing was performed using ENVI SARscape software.

Different types of errors exist in the interferograms that make undesired phases. Spatiotemporal variations of the atmosphere are the main cause of the turbulent phase delay, which may be on the order of several centimeters and often overwhelm the deformation signal of interest (Hanssen 2001; Jolivet et al. 2014). The 3D displacement vectors of GNSS stations allow detrending (unification of the datum) of the InSAR-derived results from different pairs in various environmental and geometric noises. Eliminating the long-wavelength trend from the interferograms by fitting them to the GNSS data yields a corrected result (Hanssen 2001; Normand and Heggy 2015). The best-fitted surface among all the possible surfaces (linear or nonlinear) for detrending the InSAR results could be obtained through the analysis of least-squares residuals (Samsonov et al. 2008). In some cases, with or without GNSS data, this detrending could be made through data in the far fields of the interferograms, away from the deforming area (Wright et al. 2004; Lindsey et al. 2015).

In this article, the GNSS displacement vectors of the Integrated Plate Boundary Observatory Chile (IPOC) network have been applied. A total of 10 permanent GNSS displacement vectors surrounding the Illapel earthquake epicenter were applied in this study, as illustrated in Fig. 4.

The values of the GNSS displacement vectors were derived from subtracting the recorded coordinates of stations on August 24, 2015, and September 17, 2015, based on the acquisition date of the Sentinel-1 and ALOS-2 missions covering the September 16, 2015, Illapel earthquake. According to the time interval of data acquisition, the applied displacement vectors contain co- and maybe postseismic slips, which is not a matter of concern in this study. Fig. 5 illustrates the coseismic interferograms associated with the 2015 Illapel earthquake. Each fringe in the TOPS Sentinel-1 corresponds to 2.8 cm of LOS displacement [Figs. 5(a and b)], while each one corresponds to 11.8 cm of LOS displacement in PALSAR-2 ALOS-2 [Fig. 5(c)].

As shown in Figs. 5(a–c), the deformation field comprises dense, concentric semicircular fringes that are convex toward the east; closer to the coast, fringes become more dense, which implies increasing deformation gradients. The along-track displacements of descending and ascending orbits of Sentinel-1, which were produced through burst overlap interferometry, are shown in Fig. 6.

After producing the unwrapped interferograms, the range and azimuth displacements were used to decompose 3D displacement vectors in two cases: first, when only the range displacements were applied, and second, when the images of azimuth displacements were added to the models. In the first case, the primary DoFs of the functional model with three LOS measurements and three components of displacements were zero. Consequently, the primary DoFs of the stochastic model were −2 due to $d{f}_1=0$ and the existence of two unknown variance components (one for ALOS-2 and another for Sentinel-1). In the second case, which includes three range and two azimuth measurements, and three unknown variance components (one for ALOS-2, one for range of Sentinel-1, and another for azimuth of Sentinel-1 missions), the primary DoFs of the functional and stochastic models were 2 and 0, respectively.

To check the validity of the results in each state, we compared the result of our methods with displacement values of GNSS observations. The values of range and azimuth displacements of GNSS stations were extracted and applied to decompose 3D displacement vectors for GNSS stations.

The LS-VCE is not applicable in the two aforementioned states in a traditional manner for each pixel due to lack of redundancy. To increase the redundancy of the stochastic model, the displacements in a block of $3\times 3$ pixels surrounding the pixel of the GNSS station were considered as the vectors of observations, as discussed in a previous section (Fig. 1). The DoFs of the stochastic model in this way increase manyfold in a significant manner, and consequently the estimated values of variance components converge into a precise solution in the iterated process of LS-VCE.

In this test, only range displacements of two ascending and descending orbits of Sentinel-1 and descending orbit of ALOS-2 were applied for producing 3D components of displacements. Therefore, the primary DoFs of the functional model were zero (three LOS displacements for 3D unknown displacements and two variance components). The average condition number of the covariance matrix of the parameters was about 75,000. The large condition numbers indicate that the problem of retrieving 3D displacements is an ill-posed problem. Therefore, the TR method was applied to stabilize the estimation process. In this context, these three displacements were classified into two groups on the basis of their statistical properties. The corresponding ascending and descending LOSs of Sentinel-1 were considered in the first group, while the descending ALOS-2 was in the second group. Moreover, to estimate the best values of weights for observations, the RLS-VCE method was implemented in nine neighborhood pixels. For this purpose, the regularization parameter $(\alpha )$ was estimated through the L-curve criterion for each even-determined equation system of pixels. The condition number of the covariance matrix of the TR method decreased to 50 from the 75,000 of CM on average. The primary precision of the range observations of ALOS-2 and Sentinel-1 were calculated with Bamler and Eineder's (2005) equation and the initial values of variance components were assumed to be 1. The threshold of iterations, $ϵ={10}^{-8}$, was used to estimate the variance factors. The individual variances were estimated through Eqs. (5)–(9); experimental iteration was up to 8 times. The primary standard deviations of the Sentinel-1 and ALOS-2 range displacements were approximately 0.58 and 1.61 cm, respectively. The average standard deviations, after applying the LS-VCE method, changed to 0.37 and 1.35 cm, respectively. This is what we would expect because estimating two separate variance components for the displacements of two missions will affect the contribution of the observations on each axis to the final least-squares solution. The next stage is regularizing the solution through TR. The L-curve method was applied to compute the regularization parameter for each cell solution. The result of the regularization parameter was 0.04 on average for the 10 cells of GNSS stations. The decomposed eastern, northern, and vertical components of 3D displacement vectors are shown in Figs. 7(a–c). Comparing the values of displacement in Figs. 7(a–c) indicates that the dominant component of displacement is westward with a maximum of 210 cm close to the coast. The vertical displacement map shows uplift of about 30 cm along the coast in an approximately circular region near the station CNBA, surrounded by dominant subsidence of about 20 cm. The RMSEs between the InSAR and GNSS were 0.5, 13.5, and 3.7 cm for the eastern, northern, and vertical directions, respectively, in CM, while they were 0.12, 7.8, and 1.3 cm in the RLS-VCE method. Comparing these values indicates a 75%, 40%, and 65% improvement in estimating the eastern, northern, and vertical components of the deformation field, respectively. In particular, for the northern component a decrease in error estimation from 13.5 to 7.8 cm was observed. The overall RMSE of conventional, biased, and unbiased TR methods with primary weights of measurements were 8.11, 5.76, and 5.26 cm, respectively. Compared to CM, the overall RMSEs in the biased and unbiased TR method show approximately 29% and 35% improvement, respectively. After estimating the proper weights of measurements through the LS-VCE method and applying the TR method (i.e., applying the RLS-VCE method), the RMSEs decreased to 4.92 and 4.54 cm for biased and unbiased states, respectively, equivalent to 39% and 44% improvement in retrieving the 3D displacement vector (Table 3). The standard deviation of the eastern, northern, and vertical components of GNSS stations retrieved through conventional and RLS-VCE methods were calculated through the stochastic section of Eqs. (4) and (11), respectively. These componential precision results for each station and their related average standard deviations are given in Table 4. The average of the componential standard deviation of the eastern, northern, and vertical components were 0.7, 39.5, and 9 cm for CM, while these values were 0.1, 4.0, and 0.9 cm in the RLS-VCE method. The results indicate 85%, 90%, and 90% improvement in the precision of eastern, northern, and vertical directions, respectively; consequently, the standard deviation of the northern component decreased from 39.5 to 4.0 cm. The results are listed in Table 3. Besides the standard deviations, comparing the eigenvalues of the variance–covariance matrices of the retrieved 3D displacements confirms the 75% improvement in the precision of the LS-VCE method with respect to CM.

In the second case, in which the azimuth displacements were included in the model, the condition numbers of covariance matrices dropped to about 5, indicating that adding the azimuth displacements increases the geometric strength of the configuration of the problem and changes it to a well-posed situation. Therefore, there is no need to regularize the problem and LS-VCE would work.

The retrieved 3D displacement fields of the Illapel 2015 earthquake from displacements and the LS-VCE method are shown in Fig. 8. In the case of five displacements for nine adjacent pixels, the overall RMSE of CM with primary weights and the LS-VCE method will be 1.46 and 1.34 cm, respectively. Deriving the overall RMSEs shows approximately 8% improvement in the accuracy achieved by the LS-VCE method in comparison with CM with primary weights. The standard deviation of the componential displacements through the LS-VCE method, in this scenario, are given in Table 4. The GNSS displacement vectors, retrieved 3D displacement vectors of GNSS stations with only range displacements and retrieved 3D displacement vectors with both range and azimuth displacements are shown in Fig. 9. Comparison of the arrows with range displacements and those with both azimuth and range displacements clearly shows the positive effect of adding the azimuth displacements, where those with both azimuth and range displacements are closer to the GNSS vectors than those with only range displacement. In this case, the componential RMSE of northern displacements decreased drastically to 2.2 cm from 13.5 cm (for range displacements and primary weights). This is approximately equivalent to 80% improvement in the accuracy of estimating the northern component of displacement. Achieving this accuracy for the northern component is of major interest for all disciplines of geoscience dealing with 3D surface deformation analysis.

The primary standard deviations of range (ALOS-2) and range and azimuth (Sentinel-1) measurements according to Bamler and Eineder (2005) are equal to 0.97, 0.16, and 4.5 cm, respectively. After estimating the proper weights of measurements through the LS-VCE method, these values become 2.7, 0.05, and 0.4 cm on average, respectively. These values so far are the best estimations for the precision of ALOS-2 and Sentinel-1 in range and azimuth modes. The RMSE of the RLS-VCE method decreased from 4.54 cm in the first scenario to 1.34 cm for the LS-VCE method in the second scenario; almost a 70% increase in the accuracy of estimation when the azimuth displacements are incorporated in retrieving 3D displacement vectors without any regularization.

The smaller RMSE of the inversion with primary DoFs = 2 compared to the case with primary DoFs = 0 highlights the more precise decomposition of the 3D displacements, in particular for the northern component. We observed a clear northern-southern convergent boundary around ${31.2}^{°}\text{S}$ latitude, while in the first case it is not possible to distinguish this boundary. This detected convergent boundary line is shown in Fig. 8(b), which is clearly distinguishable by considering the northern displacement of GNSS vectors, illustrated in Fig. 10. Comparing Figs. 9 and 10 with Figs. 8(a–c) clearly shows that the retrieved 3D displacement field of InSAR is consistent with GNSS displacement vectors. Furthermore, displacements of the deprecated InSAR (D-InSAR) technique provide a great volume of information about ground motion in relation to GNSS.