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

Abstract

The filtering 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 filtered by data fusion, i.e., tensor decomposition, and these filters 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 filtering directly. The filtering 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 filter 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 filter fuses the noisy interferograms and decomposes the InSAR tensor to acquire the low rank information, i.e., filtered result. Moreover, this method can not only operate on tensor model, but also efficiently filter the new acquired interferogram as matrix model with the assistance of chosen low rank information. Compared with other tensor-based filters, e.g., high order robust principal component analysis (HoRPCA) and Kronecker-basis-representation multi-pass SAR interferometry (KBR-InSAR), and the widespread traditional filters 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 effectiveness 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.

1. 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 deformation [2,3]. The damage to urban buildings and bridges caused by slow geologic change seriously threaten people’s lives and production safety [4]. Besides, with economic development and population increase, the groundwater funnels appear in urban areas because of the massive and irrational exploitation of groundwater, which impacts the sustainable development [5]. Therefore, it is significant to apply the multi-pass InSAR to monitor the ground deformation in urban areas [6]. The processing of multi-baseline InSAR stack data involves three important steps, including interferometric phase filtering, phase unwrapping, and 3D surface information inversion, which is shown in Figure 1. The phase filtering is applied to eliminate the phase discontinuity and additive noise in interferograms, and thus realizes high-precision phase estimation. The purpose of the phase unwrapping is to recover the true phase from the wrapped phase. With the unwrapping phase as input, the process of 3D surface information inversion converts the phase components to elevation model, deformation model, tomographic model, and so on. There are many effective methods to achieve 3D surface information inversion. For example, some methods utilize the reliable reference points selected by amplitude or coherent threshold to monitor surface deformation, e.g., persistent scatterer InSAR (PSI) [7] and small baseline subset (SBAS) [8]. The localization of these reference points depends on the latitude and longitude information provided by standard SAR image products. If InSAR stack data has a layover effect caused by some tall buildings, tomographic interferometry is introduced to extract the scatter points and invert them to the ortho-coordinate [9]. It can be seen that the accuracy of interferomeric phases affects the inversion of 3D surface information. Therefore, the high-precision phase estimation is a necessary step in the workflow of InSAR stack data processing. Especially in urban areas, the interferometric phase has complicated patterns and many noise factors, and it is important to filter the noisy InSAR stack data for urban interferometric big data.

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 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 results by using simulated and real data to evaluate the filter performance. The conclusion is given in the final part.

2. 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.

3. 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., $\mathcal{A}\in {ℝ}^{I_1\times I_2\times \cdot \cdot \cdot \times I_N}$, where $N$ represents the order of tensor, and $I_n$ ($n=1,2,\cdots ,N$) is the size of nth order. The mode-n unfolding of tensor $\mathcal{A}$ is an unfolding matrix along the nth mode, i.e., $A_{(n)}\in {ℝ}^{I_n\times \left(I_1\times \cdots \times I_{n-1}\times I_{n+1}\times I\cdots \times I_N\right)}$. For example, assuming that $\mathcal{A}\in {ℝ}^{I_1\times I_2\times I_3}$, and then the mode-1 unfolding of tensor $\mathcal{A}$ is $A_{(1)}\in {ℝ}^{I_1\times I_2{I}_3}$.

Some symbols appear in the following sections, and their definitions are shown as follows.

$\otimes$ represents the outer product of two matrices, e.g., $C=A_1\otimes A_2$, and the element of $C$ is calculated as

$$C(i,j)={\displaystyle \sum _k{A}_1(i,k)}{A}_2(k,j)$$

(1)

where $A_1,A_2\in {ℝ}^{I_1\times I_2}$.

$\odot$ denotes element-wise product of two matrices, e.g., $C=A_1\odot A_2$, and the definition of the element of $C$ is

$$C(i,j)=A_1(i,j)A_2(i,j)$$

(2)

$\circ$ denotes the vector outer product, e.g., $\mathcal{A}=a_1\circ a_2\circ \dots \circ a_n$, the element of $\mathcal{A}$ is calculated as

$$\mathcal{A}(i,j,\dots ,k)=a_1(i)a_2(j)\dots a_n(k)$$

(3)

${\times }_n$ denotes n-mode product, which represents the product of tensor and matrix, e.g., $\mathcal{C}=\mathcal{A}{\times }_{m}B$, the element of $\mathcal{C}\in {ℝ}^{I_1}$ is calculated as

$$\mathcal{C}(i_1,i_2,\dots ,j_m,\mathrm{..},i_n)={\displaystyle \sum _{i_m}^{}\mathcal{A}}(i_1,i_2,\dots ,i_m,\mathrm{..},i_n)B(j_m,i_m)$$

(4)

where $\mathcal{C}\in {ℝ}^{I_1\times I_2\times \dots \times J_m\times \dots \times I_n}$, $\mathcal{A}\in {ℝ}^{I_1\times I_2\times \dots \times I_m\times \dots \times I_n}$, $B\in {ℝ}^{J_m\times I_m}$.

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

$${\Vert \mathcal{A}\Vert }_F={\left({\displaystyle \sum _{i_1}{\displaystyle \sum _{i_2}\cdots {\displaystyle \sum _{i_N}|a_{i_1,i_2,\cdots ,i_N}^{}{|}^2}}}\right)}^{1/2}$$

(5)

where $a_{i_1,i_2,\cdots ,i_N}$ is the element in $N$-order tensor $\mathcal{A}$.

The L2 norm of a vector $a\in {ℝ}^n$ refers to the square root of the sum of all vector elements squares, i.e., ${\Vert a\Vert }_2$, which is Frobenius norm reduce to vectors.

The notation $diag()$ transfers a vector to a diagonal matrix.

4. Filtering Method via Online CP Decomposition for Dynamic InSAR Tensor

4.1. Signal Model

InSAR stack data can be regarded as a three-dimensional tensor [21,22], i.e., $\mathcal{T}\in {ℂ}^{I_1\times I_2\times 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).

$$\mathcal{T}=ℒ+\mathcal{N}+ℰ$$

(6)

where $ℒ$ is the low-rank tensor, i.e., the filtered result. $\mathcal{N}$ is the Gaussian noise tensor and $ℰ$ is the outlier tensor. Topography and deformation are the main factors affecting the phase. Therefore, the clean InSAR tensor $ℒ$ is represented as follows:

$$ℒ=\exp \{-j(\frac{4\pi }{\lambda d}E\otimes b+\frac{4\pi }{\lambda }D\otimes t)\}$$

(7)

where $E\in {ℝ}^{I_1\times I_2}$ is the elevation of the surface, $D\in {ℝ}^{I_1\times I_2}$ is the deformation model, $t$ is the temporal baseline, $b$ is the spatial baseline, $\lambda$ 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 $ℒ$ 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 $-\pi$ or $\pi$ 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 $-\pi$ and $\pi$, the real and the imaginary part of the InSAR tensor needs to be processed, respectively. As the phase outliers are $-\pi$ or $\pi$, 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 $\mathcal{N}$ is diffused in the real and imaginary part. Therefore, the tensor decomposition model is shown as (8).

$$\begin{array}{l}{ℛ}_{\mathcal{W}}={ℒ}_{ℛ}+{\mathcal{N}}_{ℛ}\\ {ℐ}_{\mathcal{W}}={ℒ}_{ℐ}+{\mathcal{N}}_I\end{array}$$

(8)

where $ℛ$, $ℐ$ are denoted as the real and imaginary part of the InSAR tensor. $\mathcal{W}$ indicates the position of outliers, where the missing pixels are 0 and other pixels are 1. $\mathcal{W}$ can be roughly considered as the position of phase which equals $-\pi$ or $\pi$. ${ℛ}_{\mathcal{W}}$ represents $ℛ\odot \mathcal{W}$. ${ℒ}_{ℛ}$ and ${ℒ}_{ℐ}$ are the low rank tensors of $ℛ$ and $ℐ$, respectively. ${\mathcal{N}}_{ℛ}$ and ${\mathcal{N}}_I$ are the Gaussian noise tensors of $ℛ$ and $ℐ$, respectively.

4.2. 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

$${ℒ}_{ℛ}={\displaystyle \sum _{i=1}^R{x}_i\circ y_i\circ z_i}$$

(9)

where $R$ is the CP rank of ${ℒ}_{ℛ}$. Then the online CP is represented as

$$\mathbf{L}[t]={\displaystyle \sum _{i=1}^R{z}_i\times (x_i\circ y_i)}$$

(10)

where $\mathbf{L}[t]$ is the t-th slice of ${ℒ}_{ℛ}$. The vectors acquired by decomposing $\mathbf{L}[t]$ compose $X[\tau ]$, $Y[\tau ]$ and $z[\tau ]$, i.e., $X[\tau ]=(x_1,x_2,\dots ,x_R)$, $Y[\tau ]=(y_1,y_2,\dots ,y_R)$, $z[\tau ]=(z_1,z_2,\dots ,z_R)$, where $X[\tau ]\in {ℝ}^{I_1\times R}$, $Y[\tau ]\in {ℝ}^{I_2\times R}$, $z[\tau ]\in {ℝ}^R$.

In order to solve the online tensor decomposition shown as (10), the low-rank problem of online CP decomposition is modeled as follows:

$$\begin{array}{cc}\underset{\mathbf{L}[t]}{\min }& \frac{1}{2}\end{array}{\Vert \Omega [t]\odot \left[\mathbf{R}[t]-{\displaystyle \sum _{i=1}^R{z}_i\times (x_i\circ y_i)}\right]\Vert }_F^2$$

(11)

where $\mathbf{R}[t]$ is the t-th slice of $ℛ$ and $\Omega [t]$ is the t-th slice of $\mathcal{W}$.

$X[\tau ]$, $Y[\tau ]$, and $z[\tau ]$ are initialize as follows:

$$X[0]=\begin{array}{cc}{U}_1{\Sigma }_1{}^{\frac{1}{2}}& Y[0]=\end{array}{U}_2{\Sigma }_2{}^{\frac{1}{2}}\begin{array}{cc}& z[0]=(U_3{\Sigma }_3^{\frac{1}{2}})(1,:)\end{array}$$

(12)

where $U_n$, ${\Sigma }_n$ are acquired by SVD of the mode-n unfolding of $ℛ$, i.e., $R_{(n)}=U_n{\Sigma }_n{(U_n)}^{-1}$, where $n=1,2,3$.

Algorithm 1 OLCP-InSAR for tensor model

Input:  $R[t](t=1,2,\dots ,n)$,  $\xi$

Output:  $\mathbf{L}[t]$,  $X[t]$,  $Y[t]$,  $z[t](t=1,2,\dots ,n)$,  $R$

1: Initialize $X[0]$, $Y[0]$, $z[0]$ by (11)

2: Initialize $W=0$, $updateR=true$, $step=5$, $threshold=1e-3$

3: $R\leftarrow \max (I_1,I_2)$

4: while $updataR=true||{\Vert X[t]-X[t-1]\Vert }_F^2+{\Vert Y[t]-Y[t-1]\Vert }_F^2+{\Vert z[t]-z[t-1]\Vert }_F^2>threshold$ do

5:   calculating $X[t]$, $Y[t]$, $z[t]$ by Online CP decomposition

6:   $\mathbf{L}[t]\leftarrow X[t]diag(z[t])Y{[t]}^T$

7:   if $updateR$ then

8:    if Algorithm 2($X[t]$, $Y[t]$, $z[t]$, $\xi$) then

9:      $R\leftarrow R-step$

10:    else

11:      $updateR\leftarrow false$

12:    end if

13:  end if

14: end while

15: return $\mathbf{L}[t]$, $X[t]$, $Y[t]$, $z[t]$, $R$

Algorithm 2 CP Rank estimation via multilinear PCA

Input:  $X[t]$,  $Y[t]$,  $z[t]$,  $\xi$

1: for $i=1\to R$ do

2:    decentralized $y_i^t\times (x_i\circ z_i)$ to acquire $N_i$

3:    calculating $M_i$ by (24)

4:    calculating ${\lambda }^{(1)}$, ${\lambda }^{(2)}$ by (25)

5: end for

6: calculating $W(j,k)$ by (26)

7: $\omega \leftarrow \frac{{\displaystyle \sum {\displaystyle \sum W(j,k)}}}{R^2}$

8: if $\omega <=\xi$ then

9:    return true

10:else

11:    return false

12:end if

The noisy InSAR stack data is decomposed by Equation (11) to get accurate low rank information, i.e., $X[\tau ]$, $Y[\tau ]$, and $z[\tau ]$. 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 $\mathcal{T}\in {ℝ}^{I_1\times I_2\times 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. i.e., ${\Vert X[t]-X[t-1]\Vert }_F^2+{\Vert Y[t]-Y[t-1]\Vert }_F^2+{\Vert z[t]-z[t-1]\Vert }_F^2>\mathrm{threshold}$, the remaining interferograms participate in the decomposition process until the algorithm converges.

4.3. 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:

$$\arg \underset{X,Y,z}{\min }{\displaystyle \sum _{\tau =1}^t\gamma (\tau )\cdot \left\{{\Vert \Omega [\tau ]\odot (R[\tau ]-\mathbf{X}[\tau ]\mathrm{diag}(\mathbf{z}[\tau ])\mathbf{Y}{[\tau ]}^{\mathrm{T}})\Vert }_F^2+\frac{\mu }{2}({\Vert \mathbf{X}[\tau ]\Vert }_F^2+{\Vert \mathbf{Y}[\tau ]\Vert }_F^2+{\Vert \mathbf{z}[\tau ]\Vert }_2^2)\right\}}$$

(13)

where $\gamma (\tau )$ is the selection of prior low rank information, and the definition of $\gamma (\tau )$ is shown as follows:

$$\gamma (\tau )=I_{\underset{\tau }{\min }(1-\lambda )\left|\left|b(t)\right|-\left|b(\tau )\right|\right|+\lambda \left|t-\tau \right|}(\tau )$$

(14)

where $b(t)$ is spatial baseline of the interferogram obtained at $t$ moment. $I(\cdot )$ is the indicator function. $\lambda \in \left[0,1\right]$, and $\lambda$ 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, $\gamma (\tau )$ is designed to select the closest historical interferogram to help filter the new acquired interferogram.

The closest historical interferogram acquired at $t^{\prime }$ moment is selected by $\gamma (\tau )$. With respect to $X[t^{\prime }]$ and $Y[t^{\prime }]$, $z[t]$ can be calculated by solving the sub-problem as follows:

$$\begin{array}{l}\begin{array}{cc}\underset{z[t]}{\min }& \frac{1}{2}\end{array}{\Vert \Omega [t]\odot \left[\mathbf{R}[t]-\mathbf{X}[t^{\prime }]diag(\mathbf{z}[t])\mathbf{Y}{[t^{\prime }]}^{\mathrm{T}}\right]\Vert }_F^2+\frac{\mu }{2}{\Vert \mathbf{z}[t]\Vert }_2^2\\ =\begin{array}{cc}\underset{\mathbf{z}[t]}{\min }& \frac{1}{2}\end{array}{\displaystyle \sum _{h,w}^{}{[\mathbf{R}[t](h,w)-{({\mathbf{x}}_{}^h[t^{\prime }]\odot {\mathbf{y}}_{}^w[t^{\prime }])}^{\mathrm{T}}\mathbf{z}[t]]}_{}^2}+\frac{\mu }{2}{\Vert \mathbf{z}[t]\Vert }_2^2\end{array}$$

(15)

where $1\le h\le I_1$, $1\le w\le I_2$ and $\Omega [t](h,w)\ne 0$.

According to (15), $z[t]$ can be updated as follows:

$$\mathbf{z}[t]=\frac{{\displaystyle \sum _{h,w}\mathbf{R}[t](h,w)({\mathbf{x}}_{}^h[t^{\prime }]\odot {\mathbf{y}}_{}^w[t^{\prime }])}}{{\displaystyle \sum _{h,w}^{}({\mathbf{x}}_{}^h[t^{\prime }]\odot {\mathbf{y}}_{}^w[t^{\prime }]){({\mathbf{x}}_{}^h[t^{\prime }]\odot {\mathbf{y}}_{}^w[t^{\prime }])}^T+\mu {\mathbf{I}}_R}}$$

(16)

with respect to $z[t]$ and $Y[t^{\prime }]$, $X[t]$ can be calculated by solving the sub-problem as follows:

$$\begin{array}{cc}\underset{X[t]}{\min }& \frac{1}{2}\end{array}{\Vert \Omega [t]\odot \left[\mathbf{R}[t]-\mathbf{X}[t]diag(\mathbf{z}[t]){(\mathbf{Y}[t^{\prime }])}^{\mathrm{T}}\right]\Vert }_F^2+\frac{\mu }{2}{\Vert \mathbf{X}[t]\Vert }_F^2$$

(17)

To update $X[t]$ row by row, (16) is rewritten as follows:

$$\begin{array}{cc}\underset{x_{}^h[t]}{\min }& \frac{1}{2}\end{array}{\displaystyle \sum _w^{}(}\mathbf{R}[t](h,w)-{({\mathbf{x}}_{}^h[t])}^{T}diag(\mathbf{z}[t]){\mathbf{y}}_{}^w[t^{\prime }]{)}^2+\frac{\mu }{2}{\Vert {\mathbf{x}}_{}^h[t]\Vert }_2^2$$

(18)

According to (18), ${\mathbf{x}}_{}^h[t]$ can be updated as follows:

$${\mathbf{x}}^h[t]=\frac{{\displaystyle \sum _w\mathbf{R}[t](h,w)(diag(\mathbf{z}[t]){\mathbf{y}}^w[t^{\prime }])}}{{\displaystyle \sum _w^{}(diag(\mathbf{z}[t]){\mathbf{y}}^w[t^{\prime }]){(diag(\mathbf{z}[t]){\mathbf{y}}^w[t^{\prime }])}^T+\mu {\mathbf{I}}_R}}$$

(19)

with respect to $z[t]$ and $X[t^{\prime }]$, $Y[t]$ can be calculated by solving the sub-problem as follows:

$$\begin{array}{cc}\underset{Y[t]}{\min }& \frac{1}{2}\end{array}{\Vert \Omega [t]\odot \left[\mathbf{R}[t]-\mathbf{X}[t^{\prime }]diag(\mathbf{z}[t])\mathbf{Y}{[t]}^{\mathrm{T}}\right]\Vert }_F^2+\frac{\mu }{2}{\Vert \mathbf{Y}[t]\Vert }_F^2$$

(20)

when updating $Y[t]$ row by row, (19) is rewritten as follows:

$$\begin{array}{cc}\underset{y^w[t]}{\min }& \frac{1}{2}\end{array}{\displaystyle \sum _w^{}(}\mathbf{R}[t](h,w)-{({\mathbf{x}}^h[t^{\prime }])}^{T}diag(\mathbf{z}[t]){\mathbf{y}}^w[t]{)}^2+\frac{\mu }{2}{\Vert {\mathbf{y}}^w[t]\Vert }_2^2$$

(21)

According to (21), ${\mathbf{y}}^w[t]$ can be updated as follows:

$$y^w[t]=\frac{{\displaystyle \sum _h\mathbf{R}[t](h,w)({({\mathbf{x}}^h[t^{\prime }])}^{\mathrm{T}}diag(\mathbf{z}[t]))}}{{\displaystyle \sum _h^{}({({\mathbf{x}}^h[t^{\prime }])}^{\mathrm{T}}diag(\mathbf{z}[t])){({({\mathbf{x}}^h[t^{\prime }])}^{\mathrm{T}}diag(\mathbf{z}[t]))}^{\mathrm{T}}+\mu {\mathbf{I}}_R}}$$

(22)

In conclusion, OLCP-InSAR for matrix model is shown in Algorithm 3.

If the tth new acquired interferogram is processed, $X[t]$, $Y[t]$, and $z[t]$ are acquired by updating $X[t^{\prime }]$, $Y[t^{\prime }]$, and $z[t^{\prime }]$. Instead of dealing with the whole tensor, OLCP-InSAR filter a new acquired interferogram by the low-rank information of a selected historical interferogram without reprocessing all the historical interferograms.

Algorithm 3 OLCP-InSAR for matrix model

Input:  $R[t]$,  $\Omega [t]$,  $R$,  $X[t^{\prime }]$,  $Y[t^{\prime }]$,  $z[t^{\prime }]$

Output:  $\mathbf{L}[t]$,  $X[t]$,  $Y[t]$,  $z[t]$

1: Updating $z[t]$ by (16)

2: for $h=1\to I_1$ do

3:  Updating $x^h[t]$ by (19)

4: end for

5: for $w=1\to I_2$ do

6:  Updating $y^w[t]$ by (22)

7: end for

8: $\mathbf{L}[t]\leftarrow X[t]diag(z[t])Y{[t]}^T$

9: $\mathbf{L}[t]$, $X[t]$, $Y[t]$, $z[t]$

5. 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.

5.1. 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.

5.2. 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:

$$M_i=N_i{\times }_1{U}^{{(1)}^T}{\times }_2{U}^{{(2)}^T}$$

(23)

where $U^{(1)}$ and $U^{(2)}$ contain orthogonal vectors. $N_i$ is the center processing result of $y_i^t\ast (x_i\circ 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 \underset{U^{(1)},U^{(2)}}{\max }{\displaystyle \sum _{i=1}^R{\Vert M_i\Vert }_F^2}$$

(24)

The eigen vectors of $M_i$, i.e., ${\lambda }^{(1)}$ and ${\lambda }^{(2)}$, can be calculated as follows:

$${\lambda }_j^{(1)}={\displaystyle \sum _k^{}V(j,k)}\begin{array}{cc}& \end{array}{\lambda }_k^{(2)}={\displaystyle \sum _j^{}V(j,k)}$$

(25)

where $V={\displaystyle \sum _{i=1}^R(M_i-\frac{1}{R}{\displaystyle \sum _{i=1}^R{M}_i{)}^2}}$. ${\lambda }_j^{(1)}$ is the j-th value in ${\lambda }^{(1)}$. ${\lambda }_k^{(2)}$ is the k-th value in ${\lambda }^{(2)}$.

The weight matrix $W$ can be calculated by eigen vectors as follows:

$$W(j,k)=\sqrt{{\lambda }_j^{(1)}·{\lambda }_k^{(2)}}$$

(26)

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., ${\displaystyle {\sum }_j{\displaystyle {\sum }_{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.

6. 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.

6.1. 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 $\pi$ or $-\pi$ in the real data, we randomly select 30% pixels in each slice of stimulated InSAR tensor with the value of $\pi$ or $-\pi$, 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 mean square error (MSE) of filtered tensors with different noise conditions.

Noise Intensity	Outlier Ratio	HoRPCA	WHoRPCA	KBR-InSAR	OLCP-InSAR
3 dB	30%	0.25	0.20	0.17	0.04
5 dB	30%	0.19	0.14	0.09	0.03
7 dB	30%	0.16	0.11	0.07	0.03
5 dB	20%	0.15	0.08	0.05	0.03
5 dB	40%	0.37	0.29	0.22	0.04

. 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. The number of remaining residues of filtered tensors with different noise conditions.

Noise Intensity	Outlier Ratio	HoRPCA	WHoRPCA	KBR-InSAR	OLCP-InSAR
3 dB	30%	69,525	617,750	1,039,525	287,275
5 dB	30%	65,275	539,300	262,225	226,875
7 dB	30%	61,575	480,450	307,625	187,550
5 dB	20%	65,350	453,175	363,825	222,750
5 dB	40%	67,850	637,400	1,268,250	251,825

. 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.

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. MSE of a filtered interferogram with different noise conditions.

Noise Intensity	Outlier Ratio	Goldstein	Boxcar	NL-SAR	NL-InSAR	InSAR-BM3D	OLCP-InSAR (Matrix)
3 dB	30%	0.78	0.39	0.38	0.22	0.16	0.013
5 dB	30%	0.45	0.28	0.59	0.19	0.14	0.011
7 dB	30%	0.53	0.26	0.35	0.15	0.10	0.0086
5 dB	20%	0.27	0.13	0.10	0.080	0.050	0.010
5 dB	40%	0.60	1.0	0.59	0.35	0.28	0.012

. 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. The number of residues remaining in a filtered interferogram with different noise conditions.

Noise Intensity	Outlier Ratio	Goldstein	Boxcar	NL-SAR	NL-InSAR	InSAR-BM3D	OLCP-InSAR (Matrix)
3 dB	30%	43,729	13,406	21,578	5002	2060	7595
5 dB	30%	25,230	5574	27,697	1828	5	6769
7 dB	30%	35,692	9059	21,135	3415	1759	5222
5 dB	20%	34,097	9613	10,063	3924	1837	6426
5 dB	40%	11,961	37,554	20,891	4624	1906	6542

.

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.

6.2. 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. Residues remaining in the filtered real InSAR tensor.

HoRPCA	WHoRPCA	KBR-InSAR	OLCP-InSAR
17768	42,393	56,494	28,768

. 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 Figure 18, Figure 19 and Figure 20. The residues remaining in the filtered results of the chosen interferograms are shown as

Table 6. Residues remaining in the filtered interferograms.

Slice ID	Boxcar	Goldstein	NL-InSAR	NL-SAR	InSAR-BM3D	OLCP-InSAR
1#	5178	10,815	1531	4188	3777	3063
2#	3179	5727	627	3420	2179	1593
5#	4467	7534	1267	5411	3806	3371

. 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 Figure 18 and Figure 19. However, it is hard to guarantee the detail preservation, which is intuitive in Figure 20 and Figure 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.

7. 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.