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

AbstractSynthetic aperture radar interferometry (InSAR) measures the projection of three-dimensional (3D) ground displacement in the range direction and in the azimuth direction through image proce...


Introduction
Synthetic aperture radar interferometry (InSAR) has become a widely applied technique for precise and detailed mapping of the Earth's surface deformation. The range displacement is derived from conventional interferometry, while azimuth displacement can be retrieved from multiaperture interferometry (MAI) or offset tracking. Providing the range and azimuth displacements instead of threedimensional (3D) measurements is one of the main shortcomings of InSAR. Retrieving 3D surface displacement maps from InSAR is important for constraining sources of deformation in geophysical phenomena (Lu and Dzurisin 2014). Reconstruction of 3D displacement vectors from range and azimuth measurements is a challenging issue (Hu et al. 2014b) and several approaches have been proposed in this regard. These include (1) pure InSAR-based methods with combining the multipass range and azimuth measurements without an a priori model (Fialko et al. 2001;Wright et al. 2004;Hu et al. 2014a;Grandin et al. 2016) and assessing the alternative acquisition geometries from different look or heading angles (Wright et al. 2004;Ansari et al. 2016); (2) fusion of InSAR and global navigation satellite system (GNSS) data (Samsonov et al. 2008;Hu et al. 2012a); and (3) InSAR-based methods using a priori model and a hypothesis such as geophysical models (Samieie-Esfahany et al. 2009;Motagh et al. 2010;Guglielmino et al. 2011;Shamshiri et al. 2014;Li et al. 2015;Liu et al. 2018).
From a geometric point of view, it is sufficient that the three range or azimuth displacements be available in at least two independent geometries to retrieve the components of the 3D surface displacement field (Rocca 2003). This can be achieved with the images of ascending and descending orbits of a mission and one acquisition from another mission. However, the viewing geometries of SAR satellites are nearly the same; all fly in a near-polar orbit and they are always right-looking. In this situation, the look and incidence angles and line-of-sight (LOS) directions of two SAR missions will be very close together. This limited angular diversity between the SAR acquisitions in multiple missions (e.g., Sentinel-1 and ALOS-2) weakens the strength of the constructed configuration of these multiple missions and the solution is extremely unstable in that a small change in measurement can lead to an enormous change in the estimated model. Moreover, the near-polar orbiting of SAR missions leads to low sensitivity of the north component with respect to two other components (Rocca 2003;Wright et al. 2004). For large displacement fields, range and azimuth measurements can be combined to retrieve the 3D deformation (Neelmeijer et al. 2014). In the absence of azimuth measurement, however, it is possible to stabilize the inversion process by imposing additional constraints, a process that is generally referred to as regularization.
Another effective issue is the precision of the range and azimuth measurements, which are proportional to various factors of radar wavelength and signal coherence in particular (Bamler and Eineder 2005;Ansari et al. 2016;Yang and Peng 2016), and environmental factors (Hanssen 2001;Hu et al. 2012b). The results of the leastsquares adjustment of InSAR measurements are directly influenced by the weights assigned to the observations. The given analytic error bounds in the literature are merely an estimation and prediction of observation precision. To obtain the best linear unbiased estimators (BLUEs) of the decomposed 3D displacements, a proper stochastic model of the observables is required. Therefore, InSAR measurements with different precision (weights) should adequately contribute to the best solution. Determining the exact stochastic model of InSAR observations is challenging due to inaccurate knowledge of error sources (Hu et al. 2012a;Liu et al. 2018). The proper estimation of the (co)variance components (VCs), and hence the weights of observations, and utilization of all available information are of the utmost importance (Grodecki 1998). Hence, the use of the variance component estimations (VCEs) approach could help determine the proper variances and stochastic model of InSAR, in turn improving the retrieval of 3D displacements when heterogeneous data from multiple geometries and sensors need to be combined.
Methods for estimating VCs have been and are being assessed intensively in the statistical and geodetic literature by focusing on relative weighting of heterogeneous data (Kusche 2003a). In these methods, one variance factor is estimated for each category of observations. The Helmert method (Grafarend 2006), minimum norm quadratic unbiased estimator (MINQUE) (Rao 1971), best invariant quadratic unbiased estimation (BIQUE) (Koch 1978), restricted maximum likelihood (REML) (Koch 1986), and least-squares variance component estimation (LS-VCE) (Teunissen 1988;Teunissen and Amiri-Simkooei 2008) are the most famous methods for VCE. These methods are different in the applied estimation principles, as well as the distributional assumptions that are to be made. Previous studies applied the Helmert VCE to resolve the 3D displacements of the InSAR technique (Xu et al. 2010;Hu et al. 2012aHu et al. , 2014aLiu et al. 2018). The contribution of the VCE in the improvement of precision strongly depends on our knowledge about the reality of the InSAR turbulence errors. Every provided equation for the InSAR measurement precision is the primary information. It is usual to estimate the proper weights through a process that is called VCE. However, according to the primary weights, the improvements may be significant or not.
In this study, the LS-VCE and Tikhonov regularization (TR) methods are applied to dominate the relevant limitations of retrieving components of the 3D surface displacement field through displacements of differential InSAR interferograms in zero and positive functional degrees of freedom (DoFs).

Retrieving 3D Displacement Vectors from Range and Azimuth Displacements through TR and LS-VCE
The functional model of linear Gauss-Markov equations is expressed mathematically for range, azimuth, and 3D displacement vectors as follows (Plackett 1950;Grodecki 1997): where x ¼ U e U n U u ½ T are the retrieved eastern, northern, and vertical components of the 3D displacement vector in the local surface coordinate system; y is the n Â 1 vector of observables; C y is the n Â n variance-covariance matrix of observables; A is the n Â 3 design matrix; n ¼ number of observables made at the reflection points (pixel by pixel); and the operators D and E are the dispersion and mathematical expectation operands, respectively. The rows of matrix A for range and azimuth displacements are as follows (Fialko et al. 2001): where a and k ¼ azimuth of the satellite heading (positive clockwise from the north) and incidence angles, respectively. The vector of 3D displacement and its variance-covariance matrix ðCx Þ can be retrieved through weighted least-square (WLS) adjustment as follows (Mikhail and Ackermann 1982): where P is the n Â n weight matrix of observations; and r 2 0 ¼ primary variance factor, usually assumed as r 2 0 ¼ 1. The root square of the diagonal elements of Cx yields the standard deviation of the estimated parameters and is an appropriate measure in parameters' precision. Retrieving the 3D displacement vector with primary standard deviations, as described by Bamler and Eineder (2005), will be called the conventional method (CM) subsequently. The LS-VCE is adopted to deal with the estimation of different variance components to assign appropriate variances for each category of measurements. Since LS-VCE is based on the leastsquares principle, it inherits all the well-known properties of a leastsquares estimator. Moreover, it works with a user-defined weight matrix (Amiri-Simkooei 2007). Variance components are the unknowns of the equations. Therefore, applying the LS-VCE method needs to have redundant observations and positive DoFs for both functional and stochastic models (Teunissen and Amiri-Simkooei 2008).

LS-VCE Method
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; . . . ; p, are some n Â n symmetrical primary known covariance matrices described by Bamler and Eineder (2005) and Ansari et al. (2016); and r 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 Â p normal matrix N and the p Â 1 vector l are presented as and where j ¼ 1; . . . ; p;ê ¼ P t A l are the estimated residuals; and P t A is the orthogonal projector presented as Here, N À1 naturally becomes the variance-covariance matrix of the variance components (i.e., Cr ¼ N À1 ). The covariance matrix C y and the projector P t A 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.,r) 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 r ! 0 ð Þ 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 df 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 df 2 ¼ ½df 1 ðdf 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.

TR Method
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 a ¼ positive regularization parameter; and k : k L2 is the Euclidean L 2 norm. The a parameter can be estimated through the L-curve criterion (Hansen 1992(Hansen , 2007, generalized cross validation (Wahba 1976), and VCE (Arsenin and Krianev 1992). The regularized solution of Eq. (10) and its variance-covariance matrix ðx reg and C x 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).

Results and Discussion
The performance of our proposed methods in this research was evaluated through two synthetic samples and one real data set. In order to make a more realistic simulation, the parameters of the samples were synthesized for Sentinel-1 and ALOS-2 missions, resembling the real data set of the 2015 Illapel (Chile) earthquake. In order to process the synthetic samples, the following complex functions were assumed as a known 3D displacement field in local surface coordinate system, the ðe; n; uÞ: uðx; yÞ ¼ xe Àðx 2 þy 2 Þ À2:5 x 2:5; À2:5 y 2:5 The vector of ðe; n; uÞ is the theoretical value of 3D displacement for each pixel at ðx; yÞ coordinates. The samples with 500Â 500 pixels of eastern, northern, and vertical directions of synthesized displacement components of Eq. (13) are illustrated in Fig. 2(a).
The known values of 3D displacements together with the synthesized incidence and azimuth angles were applied in simulating the range and azimuth displacements. For simplification and practical purposes, the exact values of the synthesized displacements were contaminated by spatially correlated noises, which were obtained by bivariate Gaussian random numbers with zero mean and standard deviations of 2 and 3 cm for variables, respectively. The covariance of the variables was assumed to be 0.5 cm 2 . The simulated errors were projected to azimuth and range directions on the entirety of the images and added to the exact synthesized displacements. The values of displacements in Eq. (13) are in meters; consequently, the values of the color bars in Figs. 2 and 3 are in meters. For these two figures, the horizontal and vertical axes represent the eastern and northern directions, respectively.

Case I: Evaluation of the RLS-VCE Methods through a Simulation Data Set
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 Â 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 fivethree 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.

Case II: Evaluation of LS-VCE Method through a Simulation Data Set
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.

Illapel 2015 Earthquake Data Set
In this section, the results of our methods for the 2015 M 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 descend-  ing 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 forwardlooking 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 longwavelength 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 df 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 ob- servations. 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 Â 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.
First Scenario: Evaluation of RLS-VCE Methods through Illapel 2015 Earthquake Data Set When Primary DoFs 5 0 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 illposed 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 ðaÞ was estimated through the L-curve criterion for each evendetermined 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, e ¼ 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 . Comparing the values of displacement in  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 conven- (Reprinted from Grandin et al. 2016, with permission.) tional 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°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, illus-  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.

Conclusion
The 3D displacement vectors in InSAR measurements can be retrieved through multiple InSAR measurements acquired from at least three displacements in two independent imaging geometries in a theoretical manner. However, this is a mathematically ill-posed inverse problem because SAR mission geometry renders poor sensitivity of range direction interferometry to the northern component of displacement compared to east-west and vertical. Therefore, the retrieval process of 3D displacements becomes more sensitive to observation errors.
Combining different data sets to resolve this issue requires proper treatment of the weight of observations, which otherwise will have a negative effect on both precision and accuracy of the estimated 3D displacement field. The TR and LS-VCE methods were introduced here to overcome the instability of the geometric configuration of the multiple observations and to integrate heterogeneous data by assigning the proper variance components, respectively. Assigning proper weight of observations through LS-VCE and RLS-VCE instead of CM using primary weight improves the precision and accuracy of the retrieved 3D displacement vectors.
Experimental results, obtained through synthetic and real data sets, reveal the capabilities of the TR and LS-VCE methods in improving the precision and accuracy of retrieving 3D displacement vectors by up to 80% by combining range and azimuth displace-ments of InSAR instead of only range displacements with primary weights through CM.