Multipass SAR Interferometry Based on Total Variation Regularized Robust Low Rank Tensor Decomposition

Multipass SAR interferometry (InSAR) techniques based on meter-resolution spaceborne SAR satellites, such as TerraSAR-X or COSMO-SkyMed, provide 3D reconstruction and the measurement of ground displacement over large urban areas. Conventional methods such as persistent scatterer interferometry (PSI) usually requires a fairly large SAR image stack (usually in the order of tens) to achieve reliable estimates of these parameters. Recently, low rank property in multipass InSAR data stack was explored and investigated in our previous work (J. Kang et al., “Object-based multipass InSAR via robust low-rank tensor decomposition,” IEEE Trans. Geosci. Remote Sens., vol. 56, no. 6, 2018). By exploiting this low rank prior, a more accurate estimation of the geophysical parameters can be achieved, which in turn can effectively reduce the number of interferograms required for a reliable estimation. Based on that, this article proposes a novel tensor decomposition method in a complex domain, which jointly exploits low rank and variational prior of the interferometric phase in InSAR data stacks. Specifically, a total variation (TV) regularized robust low rank tensor decomposition method is exploited for recovering outlier-free InSAR stacks. We demonstrate that the filtered InSAR data stacks can greatly improve the accuracy of geophysical parameters estimated from real data. Moreover, this article demonstrates for the first time in the community that tensor-decomposition-based methods can be beneficial for large-scale urban mapping problems using multipass InSAR. Two TerraSAR-X data stacks with large spatial areas demonstrate the promising performance of the proposed method.

geophysical parameters (namely elevation and deformation parameters) for large areas can be accordingly split into two categories: Persistent Scatterer Interferometry (PSI) [2]- [11] and Distributed Scatterer Interferometry (DSI) [12]- [18].Those methods are the backbone of data analysis based on multipass InSAR stacks and widely exploited for 3D urban reconstruction and surface displacement monitoring.
Generally, the key steps of PSI [2]- [11], [19], [20] involve PS candidate identification and parameter estimation.For example, PS pixels can be selected according to amplitude dispersion index, which can be calculated by the ratio between the temporal standard deviation and mean of the amplitudes [2].By exploiting the spatial correlation of phase measurements, Stanford method for persistent scatterers (StaMPS) [21] is applicable for selecting PS in areas undergoing nonsteady deformation without prior knowledge.Likewise, based on spatial correlation analysis, PS pairs are identified via the construction of PS arc in [22].Sublook coherence approached is proposed in [23] for point-like scatterer identification without the requirement of certain number of temporal SAR images.Methods for estimating geophysical parameters such as topography height and linear deformation rates from PS are usually based on maximum likelihood estimator (MLE) [2].In order to describe the precision of the estimated parameters, least squares ambiguity decorrelation (LAMBDA), which is originally developed for the ambiguity resolution of GPS signal, is adapted to parameter estimation for PS signals in [24].When layover phenomenon is taken into account, differential SAR tomography (D-TomoSAR) [25]- [31] was proposed for efficiently reconstructing the real 3D structure of the scene.Such technique mainly contains two steps: identification of pixels with multiple PSs and parameter estimation based on tomographic inversion.
In order to extract geophysical information from non-urban areas with DS, interferometry techniques for parameter estimation from such stochastic signals have been extensively carried out since a decade ago.Normally, statistically homogeneous pixel (SHP) selection for covariance matrix estimation and optimal phase history retrieval from such covariance matrices are the two key steps in DS interferometry.As introduced in [12], SqueeSAR exploits Kolmogorov-Smirnov (KS) test for selecting SHP with the assumption that the statistics of amplitude data can be seen as a proxy for phase stability.Composed of KS, Anderson-Darling (AD), Kullback-Leibler divergence and generalized likelihood-ratio test (GLRT), different amplitude-based methods for selecting SHP are evalu-ated in [32].Estimating optimal phase histories from covariance matrices built by the selected SHP is the second key step in DSI.The construction of covariance matrices can be considered as the generation of multimaster (MM) interferograms.In order to link all the available interferometric phases, optimal phase histories, i.e.SM phases, are then estimated from such covariance matrices.It is also well-known as phase linking or phase triangulation [12], [16], [17], [33].Then, the corresponding geophysical parameters can be reconstructed in a similar processing chain of PS signals.
Although those conventional techniques for geophysical parameter estimation do exploit information from multiple neighbouring pixels, no explicit semantic and geometric information that might be preserved in the images has been utilized.Recently, several multipass InSAR techniques have been developed based on exploiting semantic and geometric information preserved in SAR images for improving geophysical parameter estimation.In [34], Zhu et al. demonstrated that by introducing building footprints from OpenStreetMap (OSM) as prior knowledge of pixels sharing similar heights into frameworks based on joint sparse reconstruction techniques, a highly accurate tomographic reconstruction can be achieved using just six interferograms, instead of the typically-required 20-100.[35] proposed a method for multi-baseline InSAR phase unwrapping based on combining non-local denoising methods and total variation regularized spectral estimation method.In our previous work, a general framework for objectbased InSAR deformation reconstruction based on a tensormodel with a regularization term is proposed.It makes use of external semantic labels of various objects like bridges, roofs and façades, as an input for the support of the total variation regularizer [36], [37].However, it requires explicit and fairly accurate semantic labels for a reliable performance.Therefore, [1] investigated the inherent low rank property of multipass InSAR phase tensors.It allows loose semantic labels, such as a rectangle covering major part of an object, for object-based geophysical parameter reconstruction in urban areas.
As a follow-on work, we seek to develop a novel method for parameter retrieval from multipass InSAR data stacks by jointly considering the variational prior [36] and the low rank property [1] of InSAR stacks.To this end, a total variation (TV) regularized robust low rank tensor decomposition method in complex domain is proposed in this paper, in order to recover outlier-free InSAR data stacks.

B. Contributions
The contributions of this paper are summarized as follows: • Based on the prior knowledge of low rank and smoothness of multipass InSAR data stacks, a novel tensor decomposition method in complex domain is proposed, i.e. a total variation (TV) regularized robust low rank tensor decomposition, for recovering outlier-free InSAR data stacks.• The proposed method not only takes advantages of both variational prior [36] and the low rank property [1] of InSAR stacks, but also it can avoid the requirement of explicit semantic labels for object-based geophysical parameter reconstruction.

C. Structure of this paper
The rest of this paper is organized as follows.Section II introduces the notations utilized in this paper and recaps our previous work.In Section III, the proposed total variation regularized robust low rank tensor decomposition in complex domain is introduced, together with its optimization procedure.Simulated experiments are conducted in Section IV.Case studies of large-scale real data in Berlin and Las Vegas are performed in Section V. Section VI draws the conclusion of this paper.

A. Notations and Tensor Model of Multipass InSAR Data Stacks
A tensor can be considered as a multi-dimensional array.The order of a tensor is the number of its modes or dimensions.A tensor of order N in the complex domain can be denoted as X ∈ C I 1 ×I 2 •••×I N and its entries as x i 1 ,i 2 , ••• ,i N .Specifically, vector x is a tensor of order one, and matrix X can be represented as a tensor of order two.Fibers are the higherorder analogy of matrix rows and columns, which are defined by fixing every index but one.Slices of a tensor are obtained by fixing all but two indices.Matricization, also known as unfolding, is the process of reordering the elements of a tensor into a matrix.Specifically, the mode-n unfolding of tensor X is defined by X (n) that is obtained by arranging the mode-n fibers as the columns of the matrix.The utilized tensor notations are summarized in Table I.
As proposed in our previous work [1], [36], an InSAR data stack can be represented by a 3-mode tensor: G ∈ C I 1 ×I 2 ×I 3 , where I 1 and I 2 represent the spatial dimensions in range and azimuth, and I 3 denotes the number of SAR interferograms.The InSAR data tensor can be modeled by where G is the modeled InSAR data tensor, A denotes the modeled amplitude tensor, b ∈ R I 3 is the vector of the spatial baselines, τ ∈ R I 3 is a warped time variable [28], e.g.τ = t for a linear motion, and τ = sin(2π(t − t 0 )) for a seasonal motion model with temporal baseline t and time offset t 0 .S ∈ R I 1 ×I 2 and P ∈ R I 1 ×I 2 are the unknown elevation and deformation maps to be estimated, respectively, λ is the wavelength of the radar signals and r denotes the range between radar and the observed area.

B. Multipass InSAR with TV Regularizer
By integrating smoothness prior knowledge of deformation map into the parameter retrieval, [36] introduces a joint reconstruction model of object-based deformation parameters mode-n unfolding of tensor X X, Y inner product of tensor X and Y, i.e. the sum of product of their entries X F Frobenius norm of tensor X, i.e.X F = X, X vec(X) vectorization of X X 1 L 1 norm of tensor X, i.e.X 1 = vec(X) 1 X * matrix nuclear norm: the sum of its singular values, i.e.X * := i σ i ⊗ outer product element-wise product by exploiting TV regularization.Correspondingly, the objectbased model can be summarized as: where G is the observed InSAR data stack, W denotes a weighting tensor, η is the penalty parameter for balancing the two terms in (2) and f (S, P) denotes the penalty term which represents the spatial prior of S and P. Specifically, smoothness prior, such as TV norm, can be considered for urban area reconstruction.

C. Low Rank Tensor Decomposition in Multipass InSAR
Moreover, seeking to magnify the power of object-based method for multipass InSAR, we investigate the low rank property inherent in InSAR data stacks [1], according to the following knowledge: • It can be generally assumed that the elevation and deformation maps, S and P, follow certain regular structure or homogeneous pattern, because of the regular man-made structures in urban areas.
• The observed SAR images of urban object areas are usually highly correlated along the temporal dimension.By exploiting the low rank property, object-based InSAR data stacks can be robustly recovered based on robust low rank tensor decomposition: where X and Ê are the recovered outlier-free InSAR data tensor and the estimated outlier tensor, respectively.Based on this model, [1] demonstrates that reliable parameter estimation can be maintained, given loose semantic labels of objects.However, smoothness structures of multipass InSAR data stacks are not exploited in the model 3.As introduced in [36], [38], [39], geophysical paramter estimation can be enhanced by considering smoothness structures of elevation and deformation maps.

III. COMBINING TOTAL VARIATION REGULARIZED ROBUST LOW RANK TENSOR DECOMPOSITION A. Total Variation Regularized Robust Low Rank Tensor Decomposition
To this end, we develop a novel tensor decomposition method in complex domain, which jointly optimizes low rank and TV terms for recovering outlier-free InSAR data stacks.Given the observed InSAR data tensor G, it can be decomposed into two parts: a low rank tensor X and a sparse outlier tensor E. To maintain smoothness structure of InSAR stacks, the decomposition can be regularized by a TV term.Correspondingly, the proposed total variation regularized robust low rank tensor decomposition method is described by: where X 3DTV is the 3D TV term for the three-mode tensor, X * denotes the tensor nuclear norm, E 1 is the tensor L 1 norm of sparse outliers and α, β and γ are the associated parameters for controlling the balance of the three terms.X * can be calculated by the sum of the N nuclear norms of the mode-n unfoldings of X, i.e.X * = n X (n) * .The 3D TV term can be defined as:

B. Optimization by Alternating Direction Method of Multipliers (ADMM)
In order to solve the optimization problem with a TV term, we first introduce auxiliary variables Z and F , and rewrite (4) as: where ) is the first-order difference operator with respect to the i n dimension of InSAR data stack.
The optimization problem ( 6) can be solved by the framework of Alternating Direction Method of Multipliers (ADMM) [40]- [42].The corresponding constraint optimization problem can be converted into an augmented Lagrangian function, yielding where T 1 , T 2 , T 3 are the introduced dual variables and µ is the penalty parameter.ADMM takes advantage of splitting one difficult optimization problem into several subproblems, where each of them has a closed-form solution.Accordingly, the minimization of L(X, E, F , Z, T 1 , T 2 , T 3 ) with respect to each variable can be solved by optimizing the following subproblems: 1) X subproblem: By fixing the other variables, the subproblem of L with respect to X is: It can be solved by the Singular Value Thresholding (SVT) operator [43], [44] on the mode-n(n = 1, 2, 3) unfolding of the tensor , where SVT operator is defined as S µ (A) := Udiag(max(σ i − µ, 0))V with U, V and σ i obtained from Singular Value Decomposition (SVD) of the matrix A.
2) Z subproblem: By fixing the other variables, the subproblem of L with respect to Z has the following form: Then, by calculating the gradient of L with respect to Z and setting it as zero, we have: where D * (•) is the adjoint operator of D(•).According to the block-circulant structure of the matrix D * D, this inverse problem can be efficiently solved by exploiting 3D Fast Fourier Transform (FFT) and its inverse transform [45], [46].
3) F subproblem: By fixing the other variables, the subproblem of L with respect to F can be written as: This L 1 -norm-induced subproblem can be efficiently solved by applying the soft-thresholding operator defined as R γ (A) := sign(A) max(|A| − γ, 0), where denotes the elementwise product (Hadamard product) of two tensors, and |A| = sign(A) A.
4) E subproblem: By fixing the other variables, the subproblem of L with respect to E is: Likewise, this subproblem can also be solved by softthresholding operator.
5) Multiplier Updating: All the dual variables can be updated by: The detailed ADMM pseudocode for solving (6) is summarized in Algorithm 1.
Using a predefined convergence condition, the solution ( X and Ê) can be obtained, i.e. the outlier-free InSAR data tensor and the sparse outlier tensor, respectively.To this end, by applying conventional multipass InSAR techniques, e.g.PSI [2], on X, we can robustly retrieve the geophysical parameters.

IV. SIMULATIONS A. Simulation Results
We simulated a multipass InSAR data stack of 200 × 250 pixels by 29 images with the true elevation and linear deformation rate shown in Fig. 1.The simulation is comparable to real scenario of urban areas.The flat background of the elevation map and different blocks on it represent the ground and buildings with different heights, respectively.Also, as shown by the simulated deformation map, gradually increasing displacement is often observed in real data.The linear deformation rates range from −15mm/year to 15mm/year and elevations are from −100m to 100m.The spatial baseline and the temporal baseline were chosen to be comparable to those of TerraSAR-X.Uncorrelated complex circular Gaussian noise was added to the simulated stack with an signal-to-noise ratio (SNR) of 0dB i.e. following the PS model.To simulate sparse outliers in the stacks, 20% of pixels randomly selected from the stack were replaced with uniformly distributed phases.
As illustrated in Fig. 1, we compared the geophysical parameters estimated by PSI, Robust Multipass InSAR technique via Object-based low rank tensor decomposition (RoMIO) [1] and the proposed method.The parameters of the proposed method are set as α = 0.1, β = 2 and γ = 0.48, respectively.The parameter selection is discussed in the following subsection.Furthermore, as shown in Fig. 2. in order to test the capability of the proposed method for handling small stacks, we calculated standard deviation (SD) of the residuals between the estimated parameters and the ground truth with respect to decreasing number of interferograms down to 9. Besides, the performance of the proposed method against different values of SNR and percentages of outliers were tested and plotted in Fig. 3 and 4.

B. Parameter Selection
There are totally four parameters introduced in the proposed method, i.e. α, β, γ, µ, where α, β, γ control the balances of the three optimization terms and µ comes with the Lagrange multiplier terms.µ can be initially set as 10 −2 and updated in each iteration by µ := min(ηµ, µ max ), where η = 1.1.As introduced in [45], [47], [48], γ can be set as 100/ √ I 1 I 2 .In our experiences, α is selected in a range from 0 to 0.2 and β can be chosen between 0 to 10.As shown in Fig. 5, based on the simulation of Fig. 4, we performed the estimation accuracy  [1] and the proposed method.Uncorrelated complex circular Gaussian noise was added to the simulated InSAR stack with an SNR of 0dB, i.e. according to PS model.To simulate sparse outliers in the stacks, 20% of pixels randomly selected from the stack were replaced with uniformly distributed phases.It can be seen that most points cannot be correctly estimated by PSI.Especially for the estimates of the ground deformation, the increasing trend from top left to the bottom right corner is not clearly visible in the PSI result.As a comparison, both the patterns of elevation and deformation maps from RoMIO and the proposed method are more clearly displayed than PSI.However, without TV regularization, the reconstruction of some "building blocks" is more blurred in RoMIO than the proposed method, e.g. the area indicated by the red circle.As the number of interferograms utilized for the reconstruction decreases, the performances of all the methods decline, but our method can still maintain the best enhancement of the estimation accuracy.
of the parameters with respect to different values of α and β.
It can be seen that optimal α and β for this simulation are around 0.11 and 1, respectively.

C. Performance Analysis
According to the visualization results shown in Fig. 1, under SNR=0dB and 20% outliers, most points cannot be correctly estimated by PSI.Especially for the background of deformation map, the increasing trend from top left to the bottom right corner is not clearly visible in the PSI result.As a comparison, both the patterns of elevation and deformation maps from RoMIO and the proposed method are more clearly displayed than PSI.However, without TV As SNR grows, the efficiency improvement of RoMIO is less prominent than the proposed method.One plausible reason may be owing to the mitigation effect of Gaussian noise by TV regularization in the proposed method.
regularization, the reconstruction of some "building blocks" is more blurred in RoMIO than the proposed method, e.g. the area indicated by the red circle in Fig. 1, since piece-wise smoothness cannot be maintained by RoMIO.Besides, as displayed by the deformation results, non-piece-wise smoothness information can also be preserved in the proposed method.
As shown in Fig. 2, under this simulation, the improvement of the estimation accuracy by both RoMIO and the proposed method can achieve ten times better than PSI.Besides, as the number of interferograms utilized for the reconstruction decreases, the performances of all the methods decline, but our method can still maintain the best estimation accuracy.Based on Fig. 3 and 4, we can see that the proposed method can mitigate the influences from both complex Gaussian noises and outliers in the InSAR stack and accomplishes more accurate reconstruction than the other two methods.As SNR grows, the efficiency improvement of RoMIO is less prominent than the proposed method.One plausible reason may be owing to the mitigation effect of Gaussian noise by TV regularization in the proposed method.It can be also observed that the performance of PSI is more severely impacted by outliers than complex Gaussian noise.The reason lies in the fact that Periodogram exploited in PSI are only the Maximum Likelihood (ML) estimator under complex Gaussian noise.It is not robust to outliers.

1) Las Vegas:
The first study area is in Las Vegas, as demonstrated in Fig. 6.The InSAR stack contains 29 TerraSAR-X interferograms in total, with the spatial dimension of 1950 × 1950 pixels.In order to test the performance of the proposed method under a low number of interferograms, a substack with 11 interferograms were selected from the full stack.The interferograms were selected so that their spatial and temporal baselines are close to uniform distribution, which is illustrated in Fig. 7. Since this spatial area is relatively large, RoMIO and the proposed method were conducted in a slidingwindow manner, with a patch size of 100 × 100 pixels.The parameters of our method were set as α = 0.12, β = 5 and  γ = 1.The estimated elevations and linear deformation rates by PSI, RoMIO and the proposed method are displayed in Fig. 8 and 9, respectively.
2) Berlin: Another study area is in Berlin, as shown in Fig. 11.The InSAR stack totally contains in total 41 TerraSAR-X interferograms, with the spatial dimension of 3000 × 2500 pixels.A substack with 15 interferograms were selected from the full stack and the associated baselines were plotted in Fig. 12.Likewise, the patch size used in the slidingwindow processing is chosen as 200 × 200 pixels.For this area, the parameters of our method were set as α = 0.12, Proposed PSI Elevations RoMIO Fig. 8.Estimated elevation maps by PSI, RoMIO and the proposed method with 11 interferograms of one area in Las Vegas.Consistent with the simulations, the tensor-decomposition-based methods, i.e.RoMIO and the proposed method, can achieve more robust performances than PSI, since many noisy points are observed in the result of PSI.For a detailed comparison, profiles of building façade (indicated by the white arrows) are plotted in Fig. 16.

PSI Proposed
Linear deformation rates RoMIO Fig. 9.Estimated linear deformation rates by PSI and the proposed method with 11 interferograms of one area in Las Vegas.Obviously, tensor-decompositionbased methods, RoMIO and the proposed one, can better maintain the smoothness of the reconstructed deformation maps.The reconstruction results of the convention center (white rectangular) are displayed in Fig. 10.
β = 3 and γ = 0.5.The estimated elevations and amplitudes of seasonal motions by PSI, RoMIO and the proposed method are displayed in Fig. 13 and 15, respectively.

B. Performance Analysis
1) Las Vegas: As shown in Fig. 8 and 9, consistent with the simulations, the tensor-decomposition-based methods, i.e.RoMIO and the proposed method, can achieve more robust performances than PSI.In particular, both of them can maintain reliable reconstruction results with a substack of 11 interferograms.Illustrated by the deformation estimates of Las Vegas Convention Center (Fig. 10), many incorrectly estimated pixels of the central area on the roof exist in the PSI result.Compared to RoMIO, the proposed method can better estimate the flat roof areas.As marked by the red dashed circle, the group of noisy points is mitigated in the result of the proposed method.Moreover, the geometric structure of building façade can also be well preserved by the proposed method.As illustrated in Fig. 16, the elevation profiles are extracted from the results in Fig. 8 (indicated by the white arrows).It is obvious that more noisy points exist in the result of PSI than RoMIO and the proposed method.It also gives us a hint that more accurate 3D models of urban areas can be obtained by the point cloud generated from our method.Besides, the histograms of temporal coherences are displayed in Fig. 18 (Left) based on the three reconstructed results.We can see that the fitness between the filtered InSAR data stack by our method and the model does apparently increase and there are more highly coherent points in the proposed method than RoMIO.Moreover, to further assess the reconstruction  quality of the proposed method, the parameters estimated by the proposed method on the full InSAR stack were regarded as the reference, in order to compare the results of the three methods applying on a smaller InSAR stack with 11 interferorgams.The performance is demonstrated in Table II.
It can be seen that the proposed method can achieve more reliable estimates of geophysical parameters than both RoMIO and PSI.
2) Berlin: From the study area shown in Fig. 11, we can see it is mainly composed by building blocks and high-rise buildings.As demonstrated in Fig. 13 and one zoom-in area in Fig. 14, more outliers appear in the 3D reconstruction by PSI than RoMIO and the proposed method.Compared to RoMIO, the proposed method can better reconstruct road areas, since smoothness structure is able to be preserved by TV regularization.As an example shown in Fig. 13 (Middle), one road profile indicated by the red curve is extracted from the results of RoMIO and the proposed method, respectively, and displayed in Fig. 17.Obviously, piecewise smooth property can be better maintained in the proposed method than RoMIO.Moreover, Fig. 15 shows that the proposed method can produce the smoothest map of deformations than RoMIO and PSI, which indicates that incorrectly estimates can be mitigated by the proposed method.Consistent with the previous experiment, the filtered InSAR stack by the proposed method can best fit the model among the three comparing methods, which is displayed by the histograms of temporal coherences in Fig. 18 (Right).Besides, the numerical analysis is done in the same manner as the above experiment.As illustrated in Table III, the estimates from the proposed method are much closer than the other two methods given the estimates from the full stack.Fig. 16.The extracted elevation profiles from the results shown in Fig. 8 (indicated by white arrows).Besides flat areas, the geometric structure of building façade can also be well preserved by the proposed method.It also gives us a hint that more accurate 3D models of urban areas can be obtained by the point cloud generated from our method.The object-based approach in [36] contains two separate stages for the geophysical parameter estimation: tensor robust principle component analysis and the TV regularized parameter estimation.Differently with the previous approach, the proposed method integrates the two prior knowledge, i.e. the variation and low rank, into a single-stage processing.To compare the efficiencies of the two methods, we choose the same real dataset used in [36], i.e. one bridge in the central area in Berlin.As illustrated in Fig. 19, it can be seen that the two methods can achieve comparable performance.Such result in turn supports our motivation of this paper: the separated optimization steps of low rank tensor decomposition and TV regularization in [36] can be merged into a single optimization.Afterwards, the estimation of height and deformation can be done pixel by pixel, which can avoid the requirement of explicit semantic masks required in [36].This is an advantage for code parallelization in large areas processing.

VI. CONCLUSION
This paper proposed a novel tensor decomposition method in complex domain based on the prior knowledge of the low rank property and smoothness structure in multipass InSAR data stacks.Based on the proposed method, geophysical parameter estimation can be improved in real data cases, compared with conventional methods, such as PSI, and also recently proposed method -RoMIO.Demonstrated by the case study, compared with PSI, the proposed method can improve the parameter estimation by a factor of more than seven for Berlin, and ten for Las Vegas.Furthermore, this work is the first to demonstrate tensor-decomposition-based multipass InSAR techniques can be beneficial for large-scale urban mapping problems using InSAR, including 3D urban reconstruction and surface displacement monitoring.Since the proposed method is based on the assumption that the signals are similar both in spatial and time domain, the reconstruction for irregular signals e.g. a breakpoint or a sudden jump in deformation signals, may not be satisfied to the proposed optimization model.The results based on the proposed method favor the smoothness reconstruction and such sparse signals may be "inpainted" according to the neighboring signals in spatial and time domain.
As a future work, we will combine the proposed method with more advanced multipass InSAR method such as D-TomoSAR, in order to produce more accurate 3D reconstruction in urban areas.The improved 3D reconstruction can be a great input to the urban 3D model reconstruction [49].Moreover, since atmosphere phase screens (APS) in multi-temporal InSAR stack partly fulfill the assumption of the proposed model (i.e.APS is only spatially correlated but not temporally), it would be interesting to systematically investigate the performance of atmosphere signal removal based on such tensor-decomposition based method.

Fig. 1 .
Fig.1.The simulated ground truth maps of linear deformation rate and elevation, along with the estimated results by PSI, RoMIO[1] and the proposed method.Uncorrelated complex circular Gaussian noise was added to the simulated InSAR stack with an SNR of 0dB, i.e. according to PS model.To simulate sparse outliers in the stacks, 20% of pixels randomly selected from the stack were replaced with uniformly distributed phases.It can be seen that most points cannot be correctly estimated by PSI.Especially for the estimates of the ground deformation, the increasing trend from top left to the bottom right corner is not clearly visible in the PSI result.As a comparison, both the patterns of elevation and deformation maps from RoMIO and the proposed method are more clearly displayed than PSI.However, without TV regularization, the reconstruction of some "building blocks" is more blurred in RoMIO than the proposed method, e.g. the area indicated by the red circle.

Fig. 2 .
Fig.2.Plot of the estimation accuracy with respect to different numbers of interferograms.As the number of interferograms utilized for the reconstruction decreases, the performances of all the methods decline, but our method can still maintain the best enhancement of the estimation accuracy.

Fig. 3 .
Fig. 3. Plot of the estimation accuracy with respect to different values of SNR.As SNR grows, the efficiency improvement of RoMIO is less prominent than the proposed method.One plausible reason may be owing to the mitigation effect of Gaussian noise by TV regularization in the proposed method.

Fig. 4 .Fig. 5 .
Fig. 4. Plot of the estimation accuracy with respect to different percentages of outliers.It can be seen that both RoMIO and the proposed method can robustly estimate geophysical parameters.

Fig. 6 .
Fig. 6.The study area of Las Vegas shown by the mean amplitude (log scale) of a TerraSAR-X InSAR stack.

Fig. 7 .
Fig. 7.The 2D distribution of spatial and temporal baselines of the selected 11 interferograms for reconstruction.The master baseline is shown in red.

Fig. 10 .
Fig. 10.The cropped zoom-in areas of the results in Fig. 9 by the dashed white rectangular.Compared to RoMIO, the proposed method can better estimate the flat roof areas, since the group of noisy points (indicated by red dashed circle) is eliminated in the result of the proposed method.

Fig. 11 .
Fig. 11.The study area of Berlin shown by the mean amplitude (log scale) of a TerraSAR-X InSAR stack.

Fig. 12 .
Fig. 12.The 2D distribution of spatial and temporal baselines of the selected 15 interferograms for reconstruction.The master baseline is shown in red.

Fig. 13 .Fig. 14 .
Fig.13.Estimated elevation maps by PSI, RoMIO and the proposed method with 15 interferograms of one area in Berlin.Besides the reconstruction of flat areas as Las Vegas, the proposed method can also achieve the robust retrieval of this complex area composed by building blocks and high-rise buildings.For a better comparison of the three methods, one zoom-in area and one road profile are displayed in Fig.14 and 17, respectively.

Fig. 17
Fig.17.The extracted elevation profiles from the results shown in Fig.13(indicated by red curve).Obviously, the proposed TV regularized tensor decomposition method can better preserve piecewise smoothness for the 3D reconstruction of roads than RoMIO.

Fig. 18 .
Fig.17.The extracted elevation profiles from the results shown in Fig.13(indicated by red curve).Obviously, the proposed TV regularized tensor decomposition method can better preserve piecewise smoothness for the 3D reconstruction of roads than RoMIO.

Fig. 19 .
Fig. 19.(Top) Amplitudes of seasonal motion estimation on one bridge in Berlin based on the method in [36].(Down) The result based on the proposed method in this paper.

TABLE II QUANTITATIVE
STUDY FOR THE RESULTS OF LAS VEGAS DATA.THE PARAMETERS ESTIMATED BY THE PROPOSED METHOD ON THE FULL INSAR STACK WERE REGARDED AS THE REFERENCE, IN ORDER TO COMPARE THE RESULTS OF THE THREE METHODS APPLYING ON A SMALLER INSAR STACK WITH 11 INTERFERORGAMS.
15MIOFig.15.Estimated amplitudes of seasonal motions by PSI, RoMIO and the proposed method with 15 interferograms of one area in Berlin.Smoothness structure can be well maintained in the reconstructed deformation map by the proposed method.

TABLE III QUANTITATIVE
STUDY FOR THE RESULTS OF BERLIN DATA.THE PARAMETERS ESTIMATED BY THE PROPOSED METHOD ON THE FULL INSAR STACK WERE REGARDED AS THE REFERENCE, IN ORDER TO COMPARE THE RESULTS OF THE THREE METHODS APPLYING ON A SMALLER INSAR STACK WITH 15 INTERFERORGAMS.