An Optimized Filtering Method of Massive Interferometric SAR Data for Urban Areas by Online Tensor Decomposition

: The ﬁltering of multi-pass synthetic aperture radar interferometry (InSAR) stack data is a necessary preprocessing step utilized to improve the accuracy of the object-based three-dimensional information inversion in urban area. InSAR stack data is composed of multi-temporal homogeneous data, which is regarded as a third-order tensor. The InSAR tensor can be ﬁltered by data fusion, i.e., tensor decomposition, and these ﬁlters keep balance in the noise elimination and the fringe details preservation, especially with abrupt fringe change, e.g., the edge of urban structures. However, tensor decomposition based on batch processing cannot deal with few newly acquired interferograms ﬁltering directly. The ﬁltering of dynamic InSAR tensor is the inevitable challenge when processing InSAR stack data, where dynamic InSAR tensor denotes the size of InSAR tensor increases continuously due to the acquisition of new interferograms. Therefore, based on the online CANDECAMP / PARAFAC (CP) decomposition, we propose an online ﬁlter to fuse data and process the dynamic InSAR tensor, named OLCP-InSAR, which performs well especially for the urban area. In this method, CP rank is utilized to measure the tensor sparsity, which can maintain the structural features of the InSAR tensor. Additionally, CP rank estimation is applied as an important step to improve the robustness of Online CP decomposition - InSAR(OLCP-InSAR). Importing CP rank and outlier’s position as prior information, the ﬁlter fuses the noisy interferograms and decomposes the InSAR tensor to acquire the low rank information, i.e., ﬁltered result. Moreover, this method can not only operate on tensor model, but also e ﬃ ciently ﬁlter the new acquired interferogram as matrix model with the assistance of chosen low rank information. Compared with other tensor-based ﬁlters, e.g., high order robust principal component analysis (HoRPCA) and Kronecker-basis-representation multi-pass SAR interferometry (KBR-InSAR), and the widespread traditional ﬁlters operating on a single interferometric pair, e.g., Goldstein, non-local synthetic aperture radar (NL-SAR), non-local InSAR (NL-InSAR), and InSAR nonlocal block-matching 3-D (InSAR-BM3D), the e ﬀ ectiveness and robustness of OLCP-InSAR are proved in simulated and real InSAR stack data. Especially, OLCP-InSAR can maintain the fringe details at the regular building top with high noise intensity and high outlier ratio.


Introduction
The interferometric synthetic aperture radar (InSAR) is an effective remote sensing technique to obtain the multi-baseline interferograms and monitor the surface deformation via repeated observation [1]. Multi-pass InSAR with the characteristics of high resolution and wide detection range has gained a great achievement in monitoring land subsidence, landslides, and small surface Many attempts are devoted to resolving the problem of noise suppression in a single interferometric pair. The boxcar filter [10], as a common spatial domain denoising method, removes noise meanwhile sacrificing interferogram details. The Goldstein filter [11] is established in the frequency domain, and the filter parameters can be adjusted by the coherence within an interferometric pair to improve the performance of reducing noise [12]. Based on the self-similarity property of images, the nonlocal means (NL-means) [13,14] denoises by weighted averaging center pixels in several similar patches of an image. Recently, NL-means has been extended to filter an interferometric pair [15,16]. Nonlocal InSAR (NL-InSAR) [15] combines non-local thought with the appropriate phase-oriented method. Nonlocal SAR (NL-SAR) [16] further improves the robustness by adaptively selecting filter parameters in case of filtering the interferogram. Furthermore, combining with Wiener filtering [17], a nonlocal block-matching 3-D (BM3D) algorithm [18,19] is proposed and performs well in SAR interferometric phase recovery with an appropriate phaseoriented method called InSAR-BM3D [20]. These filters may obtain satisfactory denoising results in most situations, but in the region of rapid phase change, their performance seriously deteriorates [20].
In addition to the above-mentioned filtering methods for a single interferometric pair, tensor decomposition has obtained great improvements as an effective filter to fuse and make full use of the information on multiple interferometric pairs. In fact, the multi-pass InSAR stack data conforms to the three-dimensional tensor mathematical representation [21], where the first and second dimensions represent the spatial distribution of an interferogram, and the third dimension denotes the temporal variation of the interferograms, i.e., the number of interferometric pairs. The Many attempts are devoted to resolving the problem of noise suppression in a single interferometric pair. The boxcar filter [10], as a common spatial domain denoising method, removes noise meanwhile sacrificing interferogram details. The Goldstein filter [11] is established in the frequency domain, and the filter parameters can be adjusted by the coherence within an interferometric pair to improve the performance of reducing noise [12]. Based on the self-similarity property of images, the nonlocal means (NL-means) [13,14] denoises by weighted averaging center pixels in several similar patches of an image. Recently, NL-means has been extended to filter an interferometric pair [15,16]. Nonlocal InSAR (NL-InSAR) [15] combines non-local thought with the appropriate phase-oriented method. Nonlocal SAR (NL-SAR) [16] further improves the robustness by adaptively selecting filter parameters in case of filtering the interferogram. Furthermore, combining with Wiener filtering [17], a nonlocal block-matching 3-D (BM3D) algorithm [18,19] is proposed and performs well in SAR interferometric phase recovery with an appropriate phase-oriented method called InSAR-BM3D [20]. These filters may obtain satisfactory denoising results in most situations, but in the region of rapid phase change, their performance seriously deteriorates [20].
In addition to the above-mentioned filtering methods for a single interferometric pair, tensor decomposition has obtained great improvements as an effective filter to fuse and make full use of the information on multiple interferometric pairs. In fact, the multi-pass InSAR stack data conforms to the three-dimensional tensor mathematical representation [21], where the first and second dimensions represent the spatial distribution of an interferogram, and the third dimension denotes the temporal Remote Sens. 2020, 12, 2582 3 of 26 variation of the interferograms, i.e., the number of interferometric pairs. The information on the interferograms in the InSAR stack data with the inner connection has a low sparsity meanwhile the noise has a high sparsity. Therefore, the authors of [21,22] fuse InSAR stack data and decompose them into low rank tensor and outlier tensor. The tensor decomposition model is transferred to the optimization problem by using Tucker rank [23] to measure the sparsity of the InSAR tensor. Recently, in [24] the authors proposed KBR-InSAR which further decomposes the InSAR tensor into the low-rank, Gaussian noise, and sparse noise tensors. Besides, Kronecker based representation (KBR) [25,26] is imported to definite InSAR tensor rank. KBR is a combination of the zero norm of the core tensor acquired by higher order singular value decomposition (HoSVD) [27] and the relaxation of the Tucker rank [28]. These filters have achieved a good performance in filtering InSAR tensor because of multi-pass data fusion, especially at the region of phase jump.
However, the accuracy can be further improved by a suitable measure of tensor sparsity. Tucker rank and KBR, utilized in HoRPCA and KBR-InSAR, both use kernel norm of mode-i tensor unfolding to measure the sparsity of tensor [23,25]. Tensor unfolding expands tensor into the matrix form and destroys the structure of tensor, which makes it difficult to extract tensor structural information. Different from them, CANDECAMP/PARAFAC (CP) rank [29] is defined as the number of vector outer products, which keeps the tensor structure. Moreover, these tensor-based filters are based on the batch manner and required memorizing the whole tensor. And their filtering result is affected by the tensor depth (the number of interferograms in tensor). In order to ensure the accuracy of filtering, these methods need to fuse the interferograms as much as possible. Therefore, these filters cannot process the dynamic InSAR tensor, where the dynamic InSAR tensor indicates the InSAR tensor depth increases continuously due to the acquisition of new interferograms. Once a new single interferometric pair is obtained, these tensor-based filters cannot process it directly and need to re-run on the updated whole tensor containing the new pair, which is inefficient and inconvenient.
In fact, the low rank information of the accumulated InSAR tensor provides an effective reference for filtering the new acquired interferograms. Fortunately, tensor decomposition has many online forms [30][31][32] and one of them based on CP rank can maintain structural features. Therefore, in order to filter dynamic InSAR tensor, we propose an unsupervised filtering framework to realize the InSAR data fusion by using online CP decomposition. The filter fuses the multi-pass interferograms and learns the low rank information of InSAR tensoronline, which can provide assistance to process new acquired interferograms. To this end, the contributions of this paper are threefold: (1) Based on the properties of the InSAR tensor and the online CP decomposition, an effective filter is proposed, named as OLCP-InSAR, to handle the dynamic InSAR stack data. Compared with other filters, OLCP-InSAR can not only eliminate noise but also keep the fringe details well, especially for the fringe mutation at the edges of buildings in urban areas.
(2) The properties of the CP rank of InSAR tensor are analyzed and an effective method based on MPCA [33] to estimate CP rank is proposed. Compared with other CP rank estimation algorithms [34,35], this method can estimate rank online and has good robustness with high noise intensity. The estimated CP rank, as an important parameter in OLCP-InSAR, helps to improve the robustness of our filtering method.
(3) The robustness and reliability of the framework are demonstrated by simulated data. Furthermore, 10 interferograms acquired by TerraSAR-X are used as experimental data to prove the effectiveness of the proposed framework. The framework can effectively filter the dynamic InSAR tensor and improve the accuracy of object-based interferometric phase estimation, especially at the regular building top with the high noise intensity and high outlier ratio.
The rest of this paper is organized as follows. Section 2 briefly introduces the method applied in the paper and the overall workflow of online dynamic InSAR tensor filtering. Section 3 presents some mathematical tensor notations and preliminaries. Section 4 introduces OLCP-InSAR applied to the InSAR tensor. Section 5 describes the CP rank and the way to estimate the InSAR tensor CP rank. Section 6 elaborates the proposed online filtering framework in detail. Section 7 provides experimental Remote Sens. 2020, 12, 2582 4 of 26 results by using simulated and real data to evaluate the filter performance. The conclusion is given in the final part.

Materials and Method Section
From the point of view of the data model, there is an important problem to estimate the valuable information from a set of uncertain observations (including noise). The appropriate estimation method depends on the appropriate data model. It is inevitable to choose an accurate mathematical model to describe the observed data and a measurement to judge the valuable information.
Consequently, a reasonable multi-baseline InSAR stack data model is established at first, i.e., tensors [22]. The definitions about tensors are introduced in Section 3. The overall workflow of tensor decomposition is shown in Figure 2, the historical InSAR stack data is decomposed into the low rank information, and the new acquired interferograms are filtered with the assistance of low rank information. Finally, the filtered dynamic InSAR phase tensor is acquired by using this online filtering framework.

Materials and Method Section
From the point of view of the data model, there is an important problem to estimate the valuable information from a set of uncertain observations (including noise). The appropriate estimation method depends on the appropriate data model. It is inevitable to choose an accurate mathematical model to describe the observed data and a measurement to judge the valuable information. Consequently, a reasonable multi-baseline InSAR stack data model is established at first, i.e., tensors [22]. The definitions about tensors are introduced in Section 3. The overall workflow of tensor decomposition is shown in Figure 2, the historical InSAR stack data is decomposed into the low rank information, and the new acquired interferograms are filtered with the assistance of low rank information. Finally, the filtered dynamic InSAR phase tensor is acquired by using this online filtering framework.

Notions and Preliminaries
Tensor is the multi-dimensional extension of the vector and matrix, which equals to a multidimensional array. In the following sections, vectors are denoted as boldface lowercase letters, e.g., (1) Some symbols appear in the following sections, and their definitions are shown as follows.
 represents the outer product of two matrices, e.g., 1 2  C A A , and the element of C is calculated as where 12

 C A A , and the definition of the
denotes the vector outer product, e.g., 12 ...

Notions and Preliminaries
Tensor is the multi-dimensional extension of the vector and matrix, which equals to a multi-dimensional array. In the following sections, vectors are denoted as boldface lowercase letters, e.g., a. Matrices are denoted as boldface capital letters, e.g., A, where a i represents the ith row vector of A and a i represents the ith column vector of A. Tensors are denoted as boldface Euler script letters, e.g., A ∈ R I 1 ×I 2 ×···×I N , where N represents the order of tensor, and I n (n = 1, 2, · · · , N) is the size of nth order. The mode-n unfolding of tensor A is an unfolding matrix along the nth mode, i.e., A (n) ∈ R I n ×(I 1 ×···×I n−1 ×I n+1 ×I···×I N ) . For example, assuming that A ∈ R I 1 ×I 2 ×I 3 , and then the mode-1 unfolding of tensor A is A (1) ∈ R I 1 ×I 2 I 3 .
Some symbols appear in the following sections, and their definitions are shown as follows. ⊗ represents the outer product of two matrices, e.g., C = A 1 ⊗ A 2 , and the element of C is calculated as where A 1 , A 2 ∈ R I 1 ×I 2 . denotes element-wise product of two matrices, e.g., C = A 1 A 2 , and the definition of the element of C is • denotes the vector outer product, e.g., A = a 1 • a 2 • . . . • a n , the element of A is calculated as × n denotes n-mode product, which represents the product of tensor and matrix, e.g., C = A × m B, the element of C ∈ R I 1 is calculated as where C ∈ R I 1 ×I 2 ×...×J m ×...×I n , A ∈ R I 1 ×I 2 ×...×I m ×...×I n , B ∈ R J m ×I m .

of 26
The definitions of some norms applied in this paper are introduced as follows. The Frobenius norm refers to the square root of the sum of all tensor elements squares, and is defined as where a i 1 ,i 2 ,··· ,i N is the element in N-order tensor A.
The L2 norm of a vector a ∈ R n refers to the square root of the sum of all vector elements squares, i.e., a 2 , which is Frobenius norm reduce to vectors.
The notation diag() transfers a vector to a diagonal matrix.

Signal Model
InSAR stack data can be regarded as a three-dimensional tensor [21,22], i.e., T ∈ C I 1 ×I 2 ×I 3 , where I 1 and I 2 represent the spatial distribution, I 3 represents the number of interferograms. The noise in the InSAR tensor [24] includes the system thermal noise with Gaussian distribution and the outliers with uniform distribution caused by under sampling. Therefore, the noisy InSAR tensor decomposition model is shown as (6).
where L is the low-rank tensor, i.e., the filtered result. N is the Gaussian noise tensor and E is the outlier tensor. Topography and deformation are the main factors affecting the phase. Therefore, the clean InSAR tensor L is represented as follows: where E ∈ R I 1 ×I 2 is the elevation of the surface, D ∈ R I 1 ×I 2 is the deformation model, t is the temporal baseline, b is the spatial baseline, λ is the wavelength of the radar signal, and d is the distance between the radar and the observed object.
Based on the signal model, we simulated a terrain model E, a linear deformation model D and a distribution of baselines, i.e., t and b, which is shown as Figure 3. The simulated elevation model contains the most common urban topography, including a square building, two irregular structures, and a cone. The temporal baselines are chosen to be close to the uniform distribution, which is almost 0.5 mm/acquisition. The spatial baselines are randomly selected between [−100m, 100m], as shown in Figure 3c. With these inputs, a simulated InSAR tensor L is generated and its angle is shown as Figure 4a. The circular complex standard Gaussian noise with SNR of 5dB is imported to the simulated InSAR stack data, and 30% pixels in each interferogram of InSAR stack data are randomly selected and replaced by −π or π to simulate outliers. The method of noise simulation in the following section is consistent with that of this section. The angle of the noisy simulated InSAR tensor is shown as Figure 4b.
In consideration of the phase that is wrapped between −π and π, the real and the imaginary part of the InSAR tensor needs to be processed, respectively. As the phase outliers are −π or π, the outliers can be understood as the pixels missing the actual value. The outliers are caused by the influence of both real and imaginary part. The complex Gaussian noise N is diffused in the real and imaginary part. Therefore, the tensor decomposition model is shown as (8).
where R, I are denoted as the real and imaginary part of the InSAR tensor. W indicates the position of outliers, where the missing pixels are 0 and other pixels are 1. W can be roughly considered as the Remote Sens. 2020, 12, 2582 6 of 26 position of phase which equals −π or π. R W represents R W. L R and L I are the low rank tensors of R and I, respectively. N R and N I are the Gaussian noise tensors of R and I, respectively.
Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 28   In consideration of the phase that is wrapped between - and  , the real and the imaginary part of the InSAR tensor needs to be processed, respectively. As the phase outliers are - or  , the outliers can be understood as the pixels missing the actual value. The outliers are caused by the influence of both real and imaginary part. The complex Gaussian noise is diffused in the real and imaginary part. Therefore, the tensor decomposition model is shown as (8).

OLCP-InSAR for Tensor Model
The noise elimination in the InSAR stack data is a significant prerequisite for extracting geophysical information. The noisy InSAR tensor decomposition model is introduced in the previous section in detail. The dynamic InSAR tensor consists of accumulated interferograms and new acquired interferograms, and it is difficult to decompose the dynamic data by using the conventional data fusion framework. Therefore, a framework based on online CP (OLCP) tensor decomposition is proposed to fuse and filter the dynamic InSAR tensor, as shown in Figure 5. On the one hand, the tensor model of OLCP-InSAR is applied to the accumulated InSAR stack data, and the specific workflow of OLCP-InSAR tensor model is shown in Figure 6, where the estimation CP rank part will be introduced in Section 6. On the other hand, for a new acquired interferogram, the selected prior low-rank information of the historical data helps to improve the filtering accuracy of the interferogram (matrix model). The selection depends on the spatial and time baselines to the dynamic data. The detail of OLCP-InSAR tensor model is demonstrated in Section 5.
The mathematical representation of online CP decomposition model is shown as Equation (10), and it is different from the standard CP decomposition model shown as Equation (9). The online CP decomposition deal with the slices in tensor to extract the low-rank vectors instead of processing the overall tensor. The excellent property is that the low-rank vectors can be updated by following slices. The standard CP decomposition can be written as where R is the CP rank of . Then the online CP is represented as   Figure 6. Flow chart of Online CANDECAMP/PARAFAC decomposition -multi-pass SAR interferometry (OLCP-InSAR) for tensor model (Algorithm 1).
In order to solve the online tensor decomposition shown as (10), the low-rank problem of online CP decomposition is modeled as follows: z are initialize as follows: where n U , n Σ are acquired by SVD of the mode-n unfolding of , i.e., (11) 2: Initialize The standard CP decomposition can be written as

Algorithm 1 OLCP-InSAR for tensor model
where R is the CP rank of L R . Then the online CP is represented as where L[t] is the t-th slice of L R . The vectors acquired by decomposing In order to solve the online tensor decomposition shown as (10), the low-rank problem of online CP decomposition is modeled as follows: where R[t] is the t-th slice of R and Ω[t] is the t-th slice of W.

Algorithm 1 OLCP-InSAR for tensor model
calculating M i by (24)  4: calculating λ (1) , λ (2) by (25)  The noisy InSAR stack data is decomposed by Equation (11) to get accurate low rank information, i.e., X[τ], Y[τ], and z[τ]. In this algorithm, Recursive Least Squares (RLS) [36] is applied to solve the low-rank problem by updating the low rank information iteratively, and then the low rank information is combined as Equation (10) to calculate the low rank slice, which is regarded as the filtered interferogram. Therefore, OLCP-InSAR for tensor model is summarized as Algorithm 1.
Noticeably, it is necessary to determine the proper CP rank first for OLCP-InSAR. A large value of CP rank is required for preset in initial processing. The lower and upper bounds of tensor rank was studied in [37]. The CP rank of a tensor T ∈ R I 1 ×I 2 ×I 3 is smaller than max(I 1 , I 2 , I 3 ). Therefore, an algorithm to confirm the rationality of the CP rank estimation is proposed for the tensor model of OLCP-InSAR, which is introduced in detail in Section 6, i.e., Algorithm 3. Moreover, OLCP-InSAR for tensor model gradually updates the low rank information slice by slice, it may not obtain the optimized low rank information until the number of slices is sufficient.
the remaining interferograms participate in the decomposition process until the algorithm converges.

OLCP-InSAR for Matrix Model
In fact, the acquisition of the multi-pass InSAR stack data is dynamic and continuous. Therefore, it is necessary to utilize the historical data to assist in processing new observation (interferogram). It is a challenge and very difficult for the conventional filters operating on matrix model, such as Goldstein, NL-InSAR, and InSAR-BM3D. Fortunately, the fusion of historical InSAR data, i.e., online CP tensor decomposition, can extract the steady and accurate low rank information from the accumulated interferograms. Similarly, this low rank information also exists in the new acquired interferogram. In order to handle the new acquired interferogram, OLCP-InSAR for matrix model is proposed which is summarized in Figure 7. The detail of OLCP-InSAR for matrix model is demonstrated as follows. Based on online CP decomposition, the optimization problem is shown as follows: where γ(τ) is the selection of prior low rank information, and the definition of γ(τ) is shown as follows: where b(t) is spatial baseline of the interferogram obtained at t moment. I(·) is the indicator function. λ ∈ [0, 1], and λ is the trade-off parameter to adjust the importance of the time and the spatial baselines.
Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 28 demonstrated as follows. Based on online CP decomposition, the optimization problem is shown as follows: where ( )   is the selection of prior low rank information, and the definition of ( )   is shown as follows: where () t b is spatial baseline of the interferogram obtained at t moment.   I  is the indicator function.
  01   ， , and  is the trade-off parameter to adjust the importance of the time and the spatial baselines.  The difference between the obtained interferograms are caused by the terrain deformation and the variation of the spatial baseline. The linear terrain deformation depends on the change of the time baseline. Therefore, ( )   is designed to select the closest historical interferogram to help filter the new acquired interferogram.
The closest historical interferogram acquired at t moment is selected by ( ) can be calculated by solving the sub-problem as follows: The difference between the obtained interferograms are caused by the terrain deformation and the variation of the spatial baseline. The linear terrain deformation depends on the change of the time baseline. Therefore, γ(τ) is designed to select the closest historical interferogram to help filter the new acquired interferogram.
The closest historical interferogram acquired at t moment is selected by γ(τ). With respect to X[t ] and Y[t ], z[t] can be calculated by solving the sub-problem as follows: where 1 ≤ h ≤ I 1 , 1 ≤ w ≤ I 2 and Ω[t](h, w) 0. According to (15), z[t] can be updated as follows: with respect to z[t] and Y[t ], X[t] can be calculated by solving the sub-problem as follows: To update X[t] row by row, (16) is rewritten as follows: According to (18), x h [t] can be updated as follows: with respect to z[t] and X[t ], Y[t] can be calculated by solving the sub-problem as follows: when updating Y[t] row by row, (19) is rewritten as follows: According to (21), y w [t] can be updated as follows:

CP Rank Estimation for OLCP-InSAR
In OLCP-InSAR, it is required to preset an appropriate value of CP rank, however, it is a NP-hard problem for calculating CP rank directly [38,39] and it influences the performance of the filtering method by initializing with a random CP rank. Therefore, it is a key step to obtain a suitable CP rank. The authors of [34,35] estimate CP rank by initializing a large value of CP rank and gradually decreasing it to be an appropriate one with iterations. However, the methods need the overall tensor as input, which limits its application in online framework. These methods are sensitive to noise, which made them unsuitable for InSAR tensors. If using these methods to estimate the CP rank of the simulated InSAR tensor shown as Figure 4, the initial CP rank cannot obviously decrease by iterations due to the Gaussian noise and outliers. Therefore, we analyze the properties of the InSAR tensor CP rank and take the real part of the InSAR tensor as an example to explain the proposed method to estimate CP rank in detail.

CP Rank of Interferometric Phase Tensor
CP rank is the number of Kronecker bases acquired by CP decomposition, and it fuses the information in tensor and decomposes the tensor into the sum of Kronecker bases. To assess the performance of OLCP-InSAR with different predetermined CP ranks, a series of InSAR tensors with different signal-to-noise ratios (SNRs) and different ratios of outliers are simulated. MSE is calculated between the reference noise-free tensor and the filtered one, and the result curve with different CP ranks are shown in Figure 8. There is no monotonous relationship between MSE and CP ranks, however, each MSE curve has an extreme point when CP rank varies. The whole curve will deteriorate with the increase of deformation rate, and the same phenomenon also appears under different noise conditions. It is remarkable that the extreme points of these MSE curves in Figure 8b have the same x-axis coordinate, that is the MSE under different noise conditions reaches the minimum value with the same CP rank.
Supposing that the predetermined CP rank is not suitable, some noise will not be removed in the filtered interferogram in the case of higher CP rank, or fringes detail will be damaged with lower CP rank. The CP rank in the InSAR tensor means the correlation between the interferograms, which is not affected by the noise intensity. Fortunately, in view of low deformation rate, there is a closer relationship between the interferograms in urban area, which means that the CP decomposition is appropriate to handle the InSAR tensor of urban area.
Supposing that the predetermined CP rank is not suitable, some noise will not be removed in the filtered interferogram in the case of higher CP rank, or fringes detail will be damaged with lower CP rank. The CP rank in the InSAR tensor means the correlation between the interferograms, which is not affected by the noise intensity. Fortunately, in view of low deformation rate, there is a closer relationship between the interferograms in urban area, which means that the CP decomposition is appropriate to handle the InSAR tensor of urban area.

CP Rank Estimation via Multilinear PCA
A filtered interferogram is regarded as the sum of factor matrices shown as (9). If factor matrices contain a lot of redundant information, it implies the current CP rank is higher than the suitable value. The information redundancy can be measured by the correlation between the factor matrices. Principal component analysis (PCA) [40] can effectively analyze the correlation between the vectors, and multilinear PCA (MPCA) [33] extend PCA to high-dimensional information, i.e., tensor. Therefore, MPCA is applied to analyze the redundant information between the factor matrices

CP Rank Estimation via Multilinear PCA
A filtered interferogram is regarded as the sum of factor matrices shown as (9). If factor matrices contain a lot of redundant information, it implies the current CP rank is higher than the suitable value. The information redundancy can be measured by the correlation between the factor matrices. Principal component analysis (PCA) [40] can effectively analyze the correlation between the vectors, and multilinear PCA (MPCA) [33] extend PCA to high-dimensional information, i.e., tensor. Therefore, MPCA is applied to analyze the redundant information between the factor matrices acquired by OLCP-InSAR, and the process of CP rank estimation is shown in Figure 9. MPCA map an original matrix N i to feature space and acquire a new matrix M i as follows: where U (1) and U (2) contain orthogonal vectors. N i is the center processing result of y t i * (x i • z i ) acquired by OLCP-InSAR.
MPCA realize the map shown in (23) by optimizing (24). Alternate least square (ALS) is applied to acquire U (1) and U (2) . arg max The eigen vectors of M i , i.e., λ (1) and λ (2) , can be calculated as follows: j is the j-th value in λ (1) . λ (2) k is the k-th value in λ (2) . The weight matrix W can be calculated by eigen vectors as follows: where W( j, k) is the distance between M j and M k . When the distances between feature matrices are far, the information similarities are poor and the redundant information is less. Therefore, we can judge whether the CP rank is appropriate by comparing the mean value of the weight matrix, i.e., j k W( j, k)/R 2 , and the preset threshold. The algorithm is summarized as Algorithm 2, which is introduced into OLCP-InSAR to help the CP rank estimation. (1) , 1 arg max The eigen vectors of i M , i.e., (1) λ and (2) λ , can be calculated as follows: (1) j  is the j-th value in (1) λ .

Parameter initialization
Online tensor decomposition The weight matrix W can be calculated by eigen vectors as follows:

Experiment Results
In this section, the quantitative and qualitative results are presented to prove the effectiveness of OLCP-InSAR if filtering InSAR tensors or new acquired interferograms. Experiments are performed on simulated and real InSAR stack data.

Simulated Data Experiment
The first experiment is to compare the performance of OLCP-InSAR and other tensor-based filters (i.e., HoRPCA, WHoRPCA, and KBR-InSAR). The parameters of HoRPCA are set as [41], and the parameters of KBR-InSAR and WHoRPCA are set as [24]. A complex InSAR tensor is generated as experimental data with 25 interferograms and its angle form are shown as Figure 4. We impose uncorrelated complex circular Gaussian noise to the simulated InSAR tensor with an SNR of 3, 5, and 7 dB. In consideration outliers π or −π in the real data, we randomly select 30% pixels in each slice of stimulated InSAR tensor with the value of π or −π, which is denoted as 30% outlier in the following section.
MSE is an effective evaluation applied to measure the filtered results, and the evaluation results of MSE between clean references and filtered tensors with different noise conditions are shown in Table 1. The number of residues that remained in these filtered InSAR tensors is also used to compare the results acquired by these filters, and these evaluation results are shown in Table 2. In addition, for visual observation, a noisy interferogram is selected in the simulated InSAR stack data with 5dB Gaussian noise and 30% outlier to compare the filtered results as shown in Figure 10. The phase value at the positions of white short line in Figure 10 is shown as Figure 11.  Figure 10. The filtered result of tensor-based filtering methods. The white short line on the noisy interferogram is the selected profile and the filtered results of these profiles are shown in Figure 6.  As shown in Table 1, OLCP-InSAR works fairly well on filtering the noisy InSAR tensor for MSE, which means that the result acquired by OLCP-InSAR is closest to the ground truth. The same conclusion is obtained by comparing the filtered interferograms shown in Figure 10 and the filtered phase value in the typical positions shown in Figure 12. HoRPCA generates an over-smooth filtered result with the fewest residues shown in Table 2, losing many details in the interferogram. As WHoRPCA introduces the weight tensor into the optimization equation and enhances the expression ability of the equation, the details of the fringes are preserved better compared to HoRPCA. However, since WHoRPCA filters the product of the weight tensor and the noisy InSAR tensor, there are many residues remaining in its filtered results. KBR-InSAR refines the tensor decomposition model and uses KBR to measure the tensor sparsity, which balances the noise removal and the fringe preservation. However, the filtering effect deteriorates as a result of the influence of outliers, which is obvious around the edge of the square building. In consideration of the auxiliary information about outlier positions, OLCP-InSAR is not sensitive to outliers. Since there are often outliers around the edge of buildings, OLCP-InSAR is the most suitable filter to handle the InSAR tensor of urban area.
The second experiment compares the performance of OLCP-InSAR and other conventional filters operating on a single interferometric pair since OLCP-InSAR can also handle a new acquired interferogram. The parameters of the conventional methods are set as [24]. The simulation data is consistent with the data in the first experiment. The prior low-rank information is extracted by the fusion of 20 simulated interferograms. The other five interferograms in the stimulated InSAR tensor are assumed as the new acquired data and filtered in the matrix model. A slice in these five interferograms in the simulated tensor with 5dB Gaussian noise and 30% outliers and the filtered results are shown in Figure 12. The evaluation results of MSE between the clean reference and filtered interferogram in Figure 12 are shown in Table 3. The white short line on the noisy interferogram in Fgiure 12 is the selected profile and the filtered results of these profiles are shown in Figure 13. The number of residues remained in the filtered interferograms is also used to compare the result of these filters as shown in Table 4. As shown in Table 1, OLCP-InSAR works fairly well on filtering the noisy InSAR tensor for MSE, which means that the result acquired by OLCP-InSAR is closest to the ground truth. The same conclusion is obtained by comparing the filtered interferograms shown in Figure 10 and the filtered phase value in the typical positions shown in Figure 12. HoRPCA generates an over-smooth filtered result with the fewest residues shown in Table 2, losing many details in the interferogram. As WHoRPCA introduces the weight tensor into the optimization equation and enhances the expression ability of the equation, the details of the fringes are preserved better compared to HoRPCA. However, since WHoRPCA filters the product of the weight tensor and the noisy InSAR tensor, there are many residues remaining in its filtered results. KBR-InSAR refines the tensor decomposition model and uses KBR to measure the tensor sparsity, which balances the noise removal and the fringe preservation. However, the filtering effect deteriorates as a result of the influence of outliers, which is obvious around the edge of the square building. In consideration of the auxiliary information about outlier positions, OLCP-InSAR is not sensitive to outliers. Since there are often outliers around the edge of buildings, OLCP-InSAR is the most suitable filter to handle the InSAR tensor of urban area.
The second experiment compares the performance of OLCP-InSAR and other conventional filters operating on a single interferometric pair since OLCP-InSAR can also handle a new acquired interferogram. The parameters of the conventional methods are set as [24]. The simulation data is consistent with the data in the first experiment. The prior low-rank information is extracted by the fusion of 20 simulated interferograms. The other five interferograms in the stimulated InSAR tensor are assumed as the new acquired data and filtered in the matrix model. A slice in these five interferograms in the simulated tensor with 5dB Gaussian noise and 30% outliers and the filtered results are shown in Figure 12. The evaluation results of MSE between the clean reference and filtered interferogram in Figure 12 are shown in Table 3. The white short line on the noisy interferogram in Fgiure 12 is the selected profile and the filtered results of these profiles are shown in Figure 13. The number of residues remained in the filtered interferograms is also used to compare the result of these filters as shown in Table 4.
Without considering the prior knowledge of the historical data along temporal directions, the filtered results of these conventional filters are much worse than OLCP-InSAR when visualizing the performance of these method. The boxcar filter shows a clear loss of resolution in the filtered result.
The Goldstein filter and the NL-SAR have some sparse residues remaining which will influence the phase unwrapping. In the absence of phase jumps, InSAR-BM3D and NL-InSAR provide acceptable results, preserving the structure of the original phase. However, NL-InSAR and InSAR-BM3D are both based on the improvement of mean filter. Although these methods filter most residues, the fringe details are destroyed seriously, especially at the region of phase jumps. Therefore, OLCP-InSAR has obvious advantages in dealing with interferograms of urban areas considering that regular buildings are the most common terrain in cities.

Noisy Interferogram
Clean reference Goldstein Boxcar

NL-SAR NL-InSAR
InSAR-BM3D OLCP-InSAR Figure 12. The result of filters operating on an interferogram. The white short line on the noisy interferogram is the selected profile and the filtered results of these profiles are shown in Figure 13.

Real Data Experiment
The experiment on 10 SAR complex images collected from TerraSAR-X covers a large single building near the Beijing Capital International Airport. The filtering for this InSAR stack data is relatively important because of the large flow of people in this kind of building. As shown in Figure   Figure 13. The filtered phase profiles with OLCP-InSAR (blue line) and other reference filtered methods based on tensor decomposition.  Without considering the prior knowledge of the historical data along temporal directions, the filtered results of these conventional filters are much worse than OLCP-InSAR when visualizing the performance of these method. The boxcar filter shows a clear loss of resolution in the filtered result.
The Goldstein filter and the NL-SAR have some sparse residues remaining which will influence the phase unwrapping. In the absence of phase jumps, InSAR-BM3D and NL-InSAR provide acceptable results, preserving the structure of the original phase. However, NL-InSAR and InSAR-BM3D are both based on the improvement of mean filter. Although these methods filter most residues, the fringe details are destroyed seriously, especially at the region of phase jumps. Therefore, OLCP-InSAR has obvious advantages in dealing with interferograms of urban areas considering that regular buildings are the most common terrain in cities.

Real Data Experiment
The experiment on 10 SAR complex images collected from TerraSAR-X covers a large single building near the Beijing Capital International Airport. The filtering for this InSAR stack data is relatively important because of the large flow of people in this kind of building. As shown in Figure 14, the size of our study area is 425 × 449. We selected nine SAR complex images from March 2012 to August 2014 and the manually cropped the target building area as experimental data. It is worth noting that the distributed target can be extracted by the object-oriented method, such as roof and bridge. The pixels in the target are located according to the longitude and latitude information provided by standard SAR image products. The master acquisition (March 2012) was selected for the coregistration with all slave images. A complex InSAR tensor was formed by these coregistrated SLC images. The filtering for this InSAR stack data was relatively challenging on account of the low SNR in the real data and the limited number of interferograms.
Although no clean reference is provided to compute MSE for real data, the number of the residues remaining in the filtered image also provides some valuable information shown in Table 5. Furthermore, an interferogram in the real tensor is selected and its filtered results are shown as Figure 15. The value of pixels, which are near the building edge in Figure 15, is shown as Figure 16. Although HoRPCA has the fewest residues, the fringe details are sacrificed to smooth the result and remove the noise, which is obvious in Figure 16. The same conclusion can be found in experiments on simulated data. The residues remaining in the filtered result acquired by WHoRPCA are obviously reducing the effectiveness and practicability of these methods. The refine tensor decomposition model applied in KBR-InSAR improves the accuracy of filtering with few outliers. Once the number of outliers in the interferograms increases as shown in the real InSAR tensor, the performance of KBR-InSAR is significantly deteriorated. The filtered result acquired by OLCP-InSAR has relatively few residues while most fringe details are preserved.
Furthermore, an interferogram in the real tensor is selected and its filtered results are shown as Figure  15. The value of pixels, which are near the building edge in Figure 15, is shown as Figure 16. Although HoRPCA has the fewest residues, the fringe details are sacrificed to smooth the result and remove the noise, which is obvious in Figure 16. The same conclusion can be found in experiments on simulated data. The residues remaining in the filtered result acquired by WHoRPCA are obviously reducing the effectiveness and practicability of these methods. The refine tensor decomposition model applied in KBR-InSAR improves the accuracy of filtering with few outliers. Once the number of outliers in the interferograms increases as shown in the real InSAR tensor, the performance of KBR-InSAR is significantly deteriorated. The filtered result acquired by OLCP-InSAR has relatively few residues while most fringe details are preserved.      Three interferograms are selected from the real InSAR tensor as the new acquired interferogram shown as Figure 17, and other interferograms are used to extract low rank information for OLCP-InSAR. Their filtered results acquired by OLCP-InSAR and other filters operating on a single interfermetric pair are shown as Figures 18-20. The residues remaining in the filtered results of the chosen interferograms are shown as Table 6. The fifth slice in real data sensor clearly shows some structural features of the building roof, as shown in Figure 17c. Therefore, the filtered results, which should preserve these features, were selected to compare these filters, shown in Figure 21.
Consistent with the above-mentioned analysis, the boxcar filter acquires a low-resolution version of the filtered interferogram with blurred fringe edges. Goldstein filter has unsatisfied performance in the real data as a result of a relatively high signal-to-noise ratio in our real InSAR tensor. NL-SAR filters have better performance than the boxcar filter because they can preserve the details of the fringes. On the down side, this mechanism of non-local mean filter can also lead to insufficient-smooth result at the areas with many outliers. NL-InSAR provides satisfied noise elimination, offering visually appealing filtered images shown in Figures 18 and 19. However, it is hard to guarantee the detail preservation, which is intuitive in Figures 20 and 21. InSAR-BM3D and OLCP-InSAR both tend to produce an appealing result which seems balance the noise suppression and detail retention. However, as shown in Figure 21, the filtered signal acquired by InSAR-BM3D has many jumps which are obviously incorrect. Therefore, the analysis of the experiment about real data provides further support to the effectiveness of OLCP-InSAR. It is worth noting that OLCP-InSAR still maintains a reliable performance with the fusion of few accumulated interferograms. As Three interferograms are selected from the real InSAR tensor as the new acquired interferogram shown as Figure 17, and other interferograms are used to extract low rank information for OLCP-InSAR. Their filtered results acquired by OLCP-InSAR and other filters operating on a single interfermetric pair are shown as Figures 18-20. The residues remaining in the filtered results of the chosen interferograms are shown as Table 6. The fifth slice in real data sensor clearly shows some structural features of the building roof, as shown in Figure 17c. Therefore, the filtered results, which should preserve these features, were selected to compare these filters, shown in Figure 21.
Consistent with the above-mentioned analysis, the boxcar filter acquires a low-resolution version of the filtered interferogram with blurred fringe edges. Goldstein filter has unsatisfied performance in the real data as a result of a relatively high signal-to-noise ratio in our real InSAR tensor. NL-SAR filters have better performance than the boxcar filter because they can preserve the details of the fringes. On the down side, this mechanism of non-local mean filter can also lead to insufficient-smooth result at the areas with many outliers. NL-InSAR provides satisfied noise elimination, offering visually appealing filtered images shown in Figures 18 and 19. However, it is hard to guarantee the detail preservation, which is intuitive in Figures 20 and 21. InSAR-BM3D and OLCP-InSAR both tend to produce an appealing result which seems balance the noise suppression and detail retention. However, as shown in Figure 21, the filtered signal acquired by InSAR-BM3D has many jumps which are obviously incorrect. Therefore, the analysis of the experiment about real data provides further support to the effectiveness of OLCP-InSAR. It is worth noting that OLCP-InSAR still maintains a reliable performance with the fusion of few accumulated interferograms. As the number of accumulated interferograms increases, OLCP-InSAR can obtain more accurate historical low rank information to help filtering new acquired interferograms, which will be far superior to other methods operating on a single interferometric pair. the number of accumulated interferograms increases, OLCP-InSAR can obtain more accurate historical low rank information to help filtering new acquired interferograms, which will be far superior to other methods operating on a single interferometric pair.

Conclusions
The 3D information inversion in urban areas is an important research direction, but the processing of multi-pass InSAR data is difficult because of the requirement of high accurate phase

Conclusions
The 3D information inversion in urban areas is an important research direction, but the processing of multi-pass InSAR data is difficult because of the requirement of high accurate phase estimation. The multi-pass interferometric stack data can be represented as InSAR tensor, and it can be filtered by solving an optimization problem and decomposed into low-rank, Gaussian noise, and outlier tensors. These tensor-based filters outperform the conventional filtering methods due to the data fusion framework. However, with the fast-growing InSAR data, it is difficult to handle the dynamic InSAR tensor for the existed tensor-based filters. Fortunately, online tensor decomposition proposed recently motivates us to fuse the InSAR stack data and estimate the steady structural feature, i.e., the low rank information, under the online decomposition framework. Therefore, a filter based on online CP decomposition, named as OLCP-InSAR, was proposed in this paper. This novel method requires CP rank and the position of outliers as auxiliary information, where CP rank is confirmed according to the correlation of the low rank information obtained by MPCA. OLCP-InSAR has two effective forms including tensor model and matrix model. One is to deal with the accumulated InSAR tensor and the accuracy can be further improved by cycling input of the interferograms in the tensor. The other is to process the new acquired interferograms by fusing the selected low rank information. The experiments were conducted with simulated data and real InSAR tensor generated from TerraSAR-X images, which proves the effectiveness and robustness of OLCP-InSAR when dealing with the InSAR data of urban areas. Comparing with tensor-based filters, OLCP-InSAR is not sensitive to the noise conditions because of the auxiliary information about outlier positions. Compared with other filters operating on a single interferogram (matrix), OLCP-InSAR can maintain the fringe details, especially at the regular building top with the high noise intensity and high outlier ratio. In conclusion, OLCP-InSAR can effectively filter the dynamic InSAR tensor and improve the accuracy of object-based interferometric phase estimation, which is an indispensable step in 3D surface information inversion. In the future work, we will continuously analyze the structural characteristics of InSAR tensor to improve the filtering accuracy and efficiency of online tensor decomposition.