Next Article in Journal
Initial User-Centered Design of a Virtual Reality Heritage System: Applications for Digital Tourism
Next Article in Special Issue
An Improved Index for Urban Population Distribution Mapping Based on Nighttime Lights (DMSP-OLS) Data: An Experiment in Riyadh Province, Saudi Arabia
Previous Article in Journal
An Improved Conversion Relationship between Tropical Cyclone Intensity Index and Maximum Wind Speed for the Advanced Dvorak Technique in the Northwestern Pacific Ocean Using SMAP Data
Previous Article in Special Issue
A Survey of Change Detection Methods Based on Remote Sensing Images for Multi-Source and Multi-Objective Scenarios
Article

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

School of Artificial Intelligence, Beijing University of Posts and Telecommunications, Beijing 100876, China
*
Author to whom correspondence should be addressed.
Remote Sens. 2020, 12(16), 2582; https://0-doi-org.brum.beds.ac.uk/10.3390/rs12162582
Received: 4 July 2020 / Revised: 29 July 2020 / Accepted: 5 August 2020 / Published: 11 August 2020
(This article belongs to the Special Issue Data Fusion for Urban Applications)

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.
Keywords: multi-pass SAR interferometry (InSAR); data fusion; online tensor decomposition; phase filtering multi-pass SAR interferometry (InSAR); data fusion; online tensor decomposition; phase filtering

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., A 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 ) I n × ( I 1 × × I n 1 × I n + 1 × I × I N ) . For example, assuming that A I 1 × I 2 × I 3 , and then the mode-1 unfolding of tensor A is A ( 1 ) 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
C ( i , j ) = k A 1 ( i , k ) A 2 ( k , j )
where A 1 , A 2 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
C ( i , j ) = A 1 ( i , j ) A 2 ( i , j )
denotes the vector outer product, e.g., A = a 1 a 2 a n , the element of A is calculated as
A ( i , j , , k ) = a 1 ( i ) a 2 ( j ) a n ( k )
× n denotes n-mode product, which represents the product of tensor and matrix, e.g., C = A × m B , the element of C I 1 is calculated as
C ( i 1 , i 2 , , j m , .. , i n ) = i m A ( i 1 , i 2 , , i m , .. , i n ) B ( j m , i m )
where C I 1 × I 2 × × J m × × I n , A I 1 × I 2 × × I m × × I n , B J m × 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
A F = ( i 1 i 2 i N | a i 1 , i 2 , , i N | 2 ) 1 / 2
where a i 1 , i 2 , , i N is the element in N -order tensor A .
The L2 norm of a vector a 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 d i a g ( ) 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., T 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).
T = + N +
where is the low-rank tensor, i.e., the filtered result. 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 ( 4 π λ d E b + 4 π λ D t ) }
where E I 1 × I 2 is the elevation of the surface, D 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 [ 100 m , 100 m ] , 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 - π 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).
W = + N W = + N I
where , 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 position of phase which equals - π or π . W represents W . and are the low rank tensors of and , respectively. N and 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
= i = 1 R x i y i z i
where R is the CP rank of . Then the online CP is represented as
L [ t ] = i = 1 R z i × ( x i y i )
where L [ t ] is the t-th slice of . The vectors acquired by decomposing L [ t ] compose X [ τ ] , Y [ τ ] and z [ τ ] , i.e., X [ τ ] = ( x 1 , x 2 , , x R ) , Y [ τ ] = ( y 1 , y 2 , , y R ) , z [ τ ] = ( z 1 , z 2 , , z R ) , where X [ τ ] I 1 × R , Y [ τ ] I 2 × R , z [ τ ] R .
In order to solve the online tensor decomposition shown as (10), the low-rank problem of online CP decomposition is modeled as follows:
min L [ t ] 1 2 Ω [ t ] [ R [ t ] i = 1 R z i × ( x i y i ) ] F 2
where R [ t ] is the t-th slice of and Ω [ t ] is the t-th slice of W .
X [ τ ] , Y [ τ ] , and z [ τ ] are initialize as follows:
X [ 0 ] = U 1 Σ 1 1 2 Y [ 0 ] = U 2 Σ 2 1 2 z [ 0 ] = ( U 3 Σ 3 1 2 ) ( 1 , : )
where U n , Σ n are acquired by SVD of the mode-n unfolding of , i.e., R ( n ) = U n Σ n ( U n ) 1 , where n = 1 , 2 , 3 .
Algorithm 1 OLCP-InSAR for tensor model
Input:  R [ t ] ( t = 1 , 2 , , n ) ,  ξ
Output:  L [ t ] ,  X [ t ] ,  Y [ t ] ,  z [ t ] ( t = 1 , 2 , , n ) ,  R
1: Initialize X [ 0 ] , Y [ 0 ] , z [ 0 ] by (11)
2: Initialize W = 0 , u p d a t e R = t r u e , s t e p = 5 , t h r e s h o l d = 1 e 3
3:  R max ( I 1 , I 2 )
4: while u p d a t a R = t r u e | | X [ t ] X [ t 1 ] F 2 + Y [ t ] Y [ t 1 ] F 2 + z [ t ] z [ t 1 ] F 2 > t h r e s h o l d do
5:   calculating X [ t ] , Y [ t ] , z [ t ] by Online CP decomposition
6:    L [ t ] X [ t ] d i a g ( z [ t ] ) Y [ t ] T
7:   if u p d a t e R then
8:    if Algorithm 2( X [ t ] , Y [ t ] , z [ t ] , ξ ) then
9:       R R s t e p
10:    else
11:       u p d a t e R f a l s e
12:    end if
13:  end if
14: end while
15: return L [ t ] , X [ t ] , Y [ t ] , z [ t ] , R
Algorithm 2 CP Rank estimation via multilinear PCA
Input:  X [ t ] ,  Y [ t ] ,  z [ t ] ,  ξ
1: for i = 1 R do
2:    decentralized y i t × ( x i z i ) to acquire N i
3:    calculating M i by (24)
4:    calculating λ ( 1 ) , λ ( 2 ) by (25)
5: end for
6: calculating W ( j , k ) by (26)
7: ω W ( j , k ) R 2
8: if ω < = ξ 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 [ τ ] , 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 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. i.e., X [ t ] X [ t 1 ] F 2 + Y [ t ] Y [ t 1 ] F 2 + z [ t ] z [ t 1 ] F 2 > 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 min X , Y , z τ = 1 t γ ( τ ) { Ω [ τ ] ( R [ τ ] X [ τ ] diag ( z [ τ ] ) Y [ τ ] T ) F 2 + μ 2 ( X [ τ ] F 2 + Y [ τ ] F 2 + z [ τ ] 2 2 ) }
where γ ( τ ) is the selection of prior low rank information, and the definition of γ ( τ ) is shown as follows:
γ ( τ ) = I min τ ( 1 λ ) | | b ( t ) | | b ( τ ) | | + λ | t τ | ( τ )
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.
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:
min z [ t ] 1 2 Ω [ t ] [ R [ t ] X [ t ] d i a g ( z [ t ] ) Y [ t ] T ] F 2 + μ 2 z [ t ] 2 2 = min z [ t ] 1 2 h , w [ R [ t ] ( h , w ) ( x h [ t ] y w [ t ] ) T z [ t ] ] 2 + μ 2 z [ t ] 2 2
where 1 h I 1 , 1 w I 2 and Ω [ t ] ( h , w ) 0 .
According to (15), z [ t ] can be updated as follows:
z [ t ] = h , w R [ t ] ( h , w ) ( x h [ t ] y w [ t ] ) h , w ( x h [ t ] y w [ t ] ) ( x h [ t ] y w [ t ] ) T + μ I R
with respect to z [ t ] and Y [ t ] , X [ t ] can be calculated by solving the sub-problem as follows:
min X [ t ] 1 2 Ω [ t ] [ R [ t ] X [ t ] d i a g ( z [ t ] ) ( Y [ t ] ) T ] F 2 + μ 2 X [ t ] F 2
To update X [ t ] row by row, (16) is rewritten as follows:
min x h [ t ] 1 2 w ( R [ t ] ( h , w ) ( x h [ t ] ) T d i a g ( z [ t ] ) y w [ t ] ) 2 + μ 2 x h [ t ] 2 2
According to (18), x h [ t ] can be updated as follows:
x h [ t ] = w R [ t ] ( h , w ) ( d i a g ( z [ t ] ) y w [ t ] ) w ( d i a g ( z [ t ] ) y w [ t ] ) ( d i a g ( z [ t ] ) y w [ t ] ) T + μ I R
with respect to z [ t ] and X [ t ] , Y [ t ] can be calculated by solving the sub-problem as follows:
min Y [ t ] 1 2 Ω [ t ] [ R [ t ] X [ t ] d i a g ( z [ t ] ) Y [ t ] T ] F 2 + μ 2 Y [ t ] F 2
when updating Y [ t ] row by row, (19) is rewritten as follows:
min y w [ t ] 1 2 w ( R [ t ] ( h , w ) ( x h [ t ] ) T d i a g ( z [ t ] ) y w [ t ] ) 2 + μ 2 y w [ t ] 2 2
According to (21), y w [ t ] can be updated as follows:
y w [ t ] = h R [ t ] ( h , w ) ( ( x h [ t ] ) T d i a g ( z [ t ] ) ) h ( ( x h [ t ] ) T d i a g ( z [ t ] ) ) ( ( x h [ t ] ) T d i a g ( z [ t ] ) ) T + μ I R
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 ] , Y [ t ] , and z [ t ] . 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 ] ,  Ω [ t ] ,  R ,  X [ t ] ,  Y [ t ] ,  z [ t ]
Output:  L [ t ] ,  X [ t ] ,  Y [ t ] ,  z [ t ]
1: Updating z [ t ] by (16)
2: for h = 1 I 1 do
3:  Updating x h [ t ] by (19)
4: end for
5: for w = 1 I 2 do
6:  Updating y w [ t ] by (22)
7: end for
8: L [ t ] X [ t ] d i a g ( z [ t ] ) Y [ t ] T
9: 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 × 1 U ( 1 ) T × 2 U ( 2 ) T
where U ( 1 ) and U ( 2 ) contain orthogonal vectors. N i is the center processing result of y i t ( 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 U ( 1 ) , U ( 2 ) i = 1 R M i F 2
The eigen vectors of M i , i.e., λ ( 1 ) and λ ( 2 ) , can be calculated as follows:
λ j ( 1 ) = k V ( j , k ) λ k ( 2 ) = j V ( j , k )
where V = i = 1 R ( M i 1 R i = 1 R M i ) 2 . λ j ( 1 ) is the j-th value in λ ( 1 ) . λ k ( 2 ) is the k-th value in λ ( 2 ) .
The weight matrix W can be calculated by eigen vectors as follows:
W ( j , k ) = λ j ( 1 ) · λ k ( 2 )
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.

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

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

Author Contributions

Conceptualization, Y.Y. and R.W.; methodology, Y.Y. and R.W.; software, R.W.; validation, Y.Y. and R.W.; formal analysis, Y.Y. and R.W.; investigation, Y.Y. and R.W.; resources, Y.Y. and W.Z.; data curation, Y.Y.; writing—original draft preparation, R.W.; writing—review and editing, Y.Y. and R.W.; supervision, Y.Y. and W.Z.; project administration, Y.Y.; funding acquisition, Y.Y. and W.Z. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded in part by the MoE-CMCC Artificial Intelligence Project (MCM20190701), the National Key R&D Program of China under Grant 2019YFF0303300 and under Subject II No. 2019YFF0303302.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Perissin, D.; Wang, T. Time-Series InSAR Applications over Urban Areas in China. IEEE J. Sel. Top. Appl. Earth Obs. Remote. Sens. 2011, 4, 92–100. [Google Scholar] [CrossRef]
  2. Thiele, A.; Cadario, E.; Schulz, K.; Thonnessen, U.; Soergel, U. Building Recognition From Multi-Aspect High-Resolution InSAR Data in Urban Areas. IEEE Trans. Geosci. Remote. Sens. 2007, 45, 3583–3593. [Google Scholar] [CrossRef]
  3. Li, T.; Tao, L.; Jingnan, L.; Mingsheng, L. Monitoring city subsidence by D-InSAR in Tianjin area. In Proceedings of the IEEE International Geoscience & Remote Sensing Symposium, Anchorage, AK, USA, 20–24 September 2004. [Google Scholar]
  4. Mansour, M.; Morgenstern, N.R.; Martin, C.D. Expected damage from displacement of slow-moving slides. Landslides 2010, 8, 117–131. [Google Scholar] [CrossRef]
  5. Dos Santos, S.M.; Cabral, J.J.D.S.P.; Filho, I.D.D.S.P. Monitoring of soil subsidence in urban and coastal areas due to groundwater overexploitation using GPS. Nat. Hazards 2012, 64, 421–439. [Google Scholar] [CrossRef]
  6. Wu, T.; Zhang, H.; Wang, C. Ground deformation retrieval of urban and suburb areas based on multi-baseline DInSAR algorithm: A case study in Cangzhou City (China). In Proceedings of the 2007 IEEE International Geoscience and Remote Sensing Symposium, Barcelona, Spain, 23–28 July 2007. [Google Scholar]
  7. Ferretti, A.; Prati, C.; Rocca, F. Permanent scatterers in SAR interferometry. IEEE Trans. Geosci. Remote. Sens. 2001, 39, 8–20. [Google Scholar] [CrossRef]
  8. Pepe, A.; Sansosti, E.; Berardino, P.; Lanari, R. On the generation of ERS/ENVISAT DInSAR time-series via the SBAS technique. IEEE Geosci. Remote. Sens. Lett. 2005, 2, 265–269. [Google Scholar] [CrossRef]
  9. Yang, B.; Xu, H.; Liu, W.; You, Y.; Xie, X. Realistic lower bound on elevation estimation for tomographic SAR. IEEE J. Sel. Top. Appl. Earth Obs. Remote. Sens. 2018, 11, 2429–2439. [Google Scholar] [CrossRef]
  10. Roscoe, A.J.; Blair, S.M. Choice and properties of adaptive and tunable digital boxcar (moving average) filters for power systems and other signal processing applications. IEEE AMPS 2016, 1–6. [Google Scholar] [CrossRef]
  11. Goldstein, R.M.; Werner, C.L. Radar interferogram filtering for geophysical applications. Geophys. Res. Lett. 1998, 25, 4035–4038. [Google Scholar] [CrossRef]
  12. Baran, I.; Stewart, M.P.; Kampes, B.M.; Perski, Z.; Lilly, P. A modification to the Goldstein radar interferogram filter. IEEE Trans. Geosci. Remote Sens. 2003, 41, 2114–2118. [Google Scholar] [CrossRef]
  13. Buades, A.; Coll, B.; Morel, J.M. A review of image denoising algorithms, with a new one. Multiscale Model. Simul. 2005, 4, 490–530. [Google Scholar] [CrossRef]
  14. Buades, A.; Coll, B.; Morel, J.M. Nonlocal image and movie denoising. Int. J. Comput. Vis. 2008, 76, 123–139. [Google Scholar] [CrossRef]
  15. Deledalle, C.A.; Denis, L.; Tupin, F. NL-InSAR: Nonlocal interferogram estimation. IEEE Trans. Geosci. Remote. Sens. 2010, 49, 1441–1452. [Google Scholar] [CrossRef]
  16. Deledalle, C.A.; Denis, L.; Tupin, F.; Reigber, A.; Jager, M. NL-SAR: A unified nonlocal framework for resolution-preserving (Pol)(In)SAR Denoising. IEEE Trans. Geosci. Remote. Sens. 2014, 53, 2021–2038. [Google Scholar] [CrossRef]
  17. Zha, X.; Fu, R.; Dai, Z.; Liu, B. Noise reduction in interferograms using the wavelet packet transform and wiener filtering. IEEE Geosci. Remote. Sens. Lett. 2008, 5, 404–408. [Google Scholar] [CrossRef]
  18. Coupé, P.; Yger, P.; Prima, S.; Hellier, P.; Kervrann, C.; Barillot, C. An optimized blockwise nonlocal means denoising filter for 3-D magnetic resonance images. IEEE Trans. Med. Imaging 2008, 27, 425–441. [Google Scholar] [CrossRef]
  19. Dabov, K.; Foi, A.; Katkovnik, V.; Egiazarian, K. Image denoising by sparse 3-D transform-domain collaborative filtering. IEEE Trans. Image Process. 2007, 16, 2080–2095. [Google Scholar] [CrossRef]
  20. Sica, F.; Cozzolino, D.; Zhu, X.X.; Verdoliva, L.; Poggi, G. InSAR-BM3D: A Nonlocal Filter for SAR Interferometric Phase Restoration. IEEE Trans. Geosci. Remote. Sens. 2018, 56, 3456–3467. [Google Scholar] [CrossRef]
  21. Kang, J.; Wang, Y.; Körner, M.; Zhu, X.X. Robust object-based multipass InSAR deformation reconstruction. IEEE Trans. Geosci. Remote. Sens. 2017, 55, 4239–4251. [Google Scholar] [CrossRef]
  22. Kang, J.; Wang, Y.; Schmitt, M.; Zhu, X.X. Object-based multipass InSAR via robust low-rank tensor decomposition. IEEE Trans. Geosci. Remote. Sens. 2018, 56, 3062–3077. [Google Scholar] [CrossRef]
  23. Tucker, L.R. Some mathematical notes on three-mode factor analysis. Psychometrika 1966, 31, 279–311. [Google Scholar] [CrossRef] [PubMed]
  24. Yanan, Y.; Wang, R.; Zhou, W. A phase filter for multi-pass InSAR stack data by hybrid tensor rank representation. IEEE Access 2019, 7, 135176–135191. [Google Scholar] [CrossRef]
  25. Xie, Q.; Zhao, Q.; Meng, D.; Xu, Z. Kronecker-basis-representation based tensor sparsity and its applications to tensor recovery. IEEE Trans. Pattern Anal. Mach. Intell. 2017, 40, 1888–1902. [Google Scholar] [CrossRef] [PubMed]
  26. Xie, Q.; Zhao, Q.; Meng, D.; Xu, Z.; Gu, S.; Zuo, W.; Zhang, K. Multispectral images denoising by intrinsic tensor sparsity regularization. In Proceedings of the 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR); Institute of Electrical and Electronics Engineers (IEEE): Las Vegas, NV, USA, 2016; pp. 1692–1700. [Google Scholar]
  27. Kolda, T.; Bader, B.W. Tensor decompositions and applications. SIAM Rev. 2009, 51, 455–500. [Google Scholar] [CrossRef]
  28. Liu, J.; Musialski, P.; Wonka, P.; Ye, J. Tensor completion for estimating missing values in visual data. IEEE TPAMI 2013, 35, 208–220. [Google Scholar] [CrossRef]
  29. Hitchcock, F.L. The expression of a tensor or a polyadic as a sum of products. J. Math. Phys. 1927, 6, 164–189. [Google Scholar] [CrossRef]
  30. Ge, R.; Huang, F.; Jin, C.; Yuan, Y. Escaping from saddle points—Online stochastic gradient for tensor decomposition. In Proceedings of the Conference on Learning Theory, Paris, France, 3–6 July 2015; pp. 797–842. [Google Scholar]
  31. Kasai, H. Online low-rank tensor subspace tracking from incomplete data by CP decomposition using recursive least squares. IEEE ICASSP 2016, 2519–2523. [Google Scholar] [CrossRef]
  32. Narayanamurthy, P.; Vaswani, N. Provable dynamic robust PCA or robust subspace tracking. IEEE ISIT 2018, 65, 376–380. [Google Scholar] [CrossRef]
  33. Lu, H.; Plataniotis, K.; Venetsanopoulos, A. MPCA: Multilinear principal component analysis of tensor objects. IEEE Trans. Neural Netw 2008, 19, 18–39. [Google Scholar] [CrossRef]
  34. Zhao, Q.; Zhang, D.; Cichocki, A. Bayesian CP Factorization of incomplete tensors with automatic rank determination. IEEE Trans. Pattern Anal. Mach. Intell. 2015, 37, 1751–1763. [Google Scholar] [CrossRef]
  35. Zhao, Q.; Zhou, G.; Zhang, L.; Cichocki, A.; Amari, S.I. Bayesian robust tensor factorization for incomplete multiway data. IEEE Trans. Neural Netw. Learn. Syst. 2016, 27, 736–748. [Google Scholar] [CrossRef] [PubMed]
  36. Haykin, S.O. Adaptive Filter Theory; Prentice Hall: Upper Saddle River, NJ, USA, 2002. [Google Scholar]
  37. Alexeev, B.; Forbes, M.A.; Tsimerman, J. Tensor rank: Some lower and upper bounds. In Proceedings of the 2011 IEEE 26th Annual Conference on Computational Complexity; Institute of Electrical and Electronics Engineers (IEEE): San Jose, CA, USA, 2011; pp. 283–291. [Google Scholar]
  38. Håstad, J. Tensor rank is NP-complete. J. Algorithms 1990, 11, 644–654. [Google Scholar] [CrossRef]
  39. Hillar, C.J.; Lim, L.H. Most tensor problems Are NP-Hard. J. ACM 2013, 60, 1–39. [Google Scholar] [CrossRef]
  40. Jolliffe, I.T. Principal Component Analysis; Wiley: Hoboken, NJ, USA, 2002. [Google Scholar]
  41. Sobral, A.; Thierry, B.; El-hadi, Z. LRSLibrary: Low-Rank and Sparse Tools for Background Modeling and Subtraction in Videos; CRC Press, Taylor and Francis Group: Boca Raton, FL, USA, 2015. [Google Scholar]
Figure 1. Flow chart of multi-pass SAR interferometry (InSAR) stack data processing.
Figure 1. Flow chart of multi-pass SAR interferometry (InSAR) stack data processing.
Remotesensing 12 02582 g001
Figure 2. The overall workflow of online dynamic InSAR tensor filtering.
Figure 2. The overall workflow of online dynamic InSAR tensor filtering.
Remotesensing 12 02582 g002
Figure 3. Illustration of a simulated example of (a) elevation model; (b) deformation model; (c) the 2-D distribution of temporal baseline and spatial baseline. The temporal baselines satisfy the uniform distribution. The spatial baselines are chosen as random distribution between [−100,100] m.
Figure 3. Illustration of a simulated example of (a) elevation model; (b) deformation model; (c) the 2-D distribution of temporal baseline and spatial baseline. The temporal baselines satisfy the uniform distribution. The spatial baselines are chosen as random distribution between [−100,100] m.
Remotesensing 12 02582 g003
Figure 4. Illustration of the simulated InSAR phase tensor (partly show). (a) The interferograms in stimulated noise-free phase tensor; (b) the interferograms in stimulated noisy phase tensor.
Figure 4. Illustration of the simulated InSAR phase tensor (partly show). (a) The interferograms in stimulated noise-free phase tensor; (b) the interferograms in stimulated noisy phase tensor.
Remotesensing 12 02582 g004
Figure 5. The filtering framework for the dynamic InSAR tensor.
Figure 5. The filtering framework for the dynamic InSAR tensor.
Remotesensing 12 02582 g005
Figure 6. Flow chart of Online CANDECAMP/PARAFAC decomposition - multi-pass SAR interferometry (OLCP-InSAR) for tensor model (Algorithm 1).
Figure 6. Flow chart of Online CANDECAMP/PARAFAC decomposition - multi-pass SAR interferometry (OLCP-InSAR) for tensor model (Algorithm 1).
Remotesensing 12 02582 g006
Figure 7. Flow chart of OLCP-InSAR for matrix model (Algorithm 2).
Figure 7. Flow chart of OLCP-InSAR for matrix model (Algorithm 2).
Remotesensing 12 02582 g007
Figure 8. MSE between the low rank tensor acquired by CP decomposition and the noisy tensor with different CP ranks. (a) InSAR tensors with different deformation rates. (b) InSAR tensors with different noise conditions.
Figure 8. MSE between the low rank tensor acquired by CP decomposition and the noisy tensor with different CP ranks. (a) InSAR tensors with different deformation rates. (b) InSAR tensors with different noise conditions.
Remotesensing 12 02582 g008
Figure 9. Flow chart of CP Rank estimation (Algorithm 2).
Figure 9. Flow chart of CP Rank estimation (Algorithm 2).
Remotesensing 12 02582 g009
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.
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.
Remotesensing 12 02582 g010
Figure 11. The filtered phase profiles with OLCP-InSAR (blue line) and other references. (a) The filtered phase profiles at row 125 and column 25 to 75. (b) The filtered phase profiles at row 275 and column 360 to 410. (c) The filtered phase profiles at row 400 and column 100 to 150. (d) The filtered phase profiles at row 150 and column 375 to 425.
Figure 11. The filtered phase profiles with OLCP-InSAR (blue line) and other references. (a) The filtered phase profiles at row 125 and column 25 to 75. (b) The filtered phase profiles at row 275 and column 360 to 410. (c) The filtered phase profiles at row 400 and column 100 to 150. (d) The filtered phase profiles at row 150 and column 375 to 425.
Remotesensing 12 02582 g011
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.
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.
Remotesensing 12 02582 g012
Figure 13. The filtered phase profiles with OLCP-InSAR (blue line) and other reference filtered methods based on tensor decomposition. (a) The filtered phase profiles at row 125 and column 25 to 75. (b) The filtered phase profiles at row 275 and column 360 to 410. (c) The filtered phase profiles at row 400 and column 100 to 150. (d) The filtered phase profiles at row 150 and column 375 to 425.
Figure 13. The filtered phase profiles with OLCP-InSAR (blue line) and other reference filtered methods based on tensor decomposition. (a) The filtered phase profiles at row 125 and column 25 to 75. (b) The filtered phase profiles at row 275 and column 360 to 410. (c) The filtered phase profiles at row 400 and column 100 to 150. (d) The filtered phase profiles at row 150 and column 375 to 425.
Remotesensing 12 02582 g013
Figure 14. The selected area in TerraSAR-X image and Google Earth. The selected building is a large single building in the city. (a) TerraSAR-X image. (b) Google Earth.
Figure 14. The selected area in TerraSAR-X image and Google Earth. The selected building is a large single building in the city. (a) TerraSAR-X image. (b) Google Earth.
Remotesensing 12 02582 g014
Figure 15. The filtered results acquired by OLCP-InSAR and other filters based on tensor decomposition. The white short line on the noisy interferogram is the selected profile and the filtered results of these profiles are shown in Figure 16.
Figure 15. The filtered results acquired by OLCP-InSAR and other filters based on tensor decomposition. The white short line on the noisy interferogram is the selected profile and the filtered results of these profiles are shown in Figure 16.
Remotesensing 12 02582 g015
Figure 16. The filtered phase profiles with OLCP-InSAR (red line) and other reference filtered methods based on tensor decomposition. (a) The filtered phase profiles at row 130 and column 320 to 370. (b) The filtered phase profiles at row 120 and column 300 to 350. (c) The filtered phase profiles at row 100 and column 160 to 210. (d) The filtered phase profiles at row 160 and column 60 to 110.
Figure 16. The filtered phase profiles with OLCP-InSAR (red line) and other reference filtered methods based on tensor decomposition. (a) The filtered phase profiles at row 130 and column 320 to 370. (b) The filtered phase profiles at row 120 and column 300 to 350. (c) The filtered phase profiles at row 100 and column 160 to 210. (d) The filtered phase profiles at row 160 and column 60 to 110.
Remotesensing 12 02582 g016
Figure 17. The interferograms in the real InSAR tensor. Color bar of the range value is shown for each channel with pointers to indicate the true values. There are small visual differences between slices. For example, the phase pattern of (c) is strongly related to the structure of the roof. (a) The first slice. (b) The second slice. (c) The fourth slice of the real InSAR tensor.
Figure 17. The interferograms in the real InSAR tensor. Color bar of the range value is shown for each channel with pointers to indicate the true values. There are small visual differences between slices. For example, the phase pattern of (c) is strongly related to the structure of the roof. (a) The first slice. (b) The second slice. (c) The fourth slice of the real InSAR tensor.
Remotesensing 12 02582 g017
Figure 18. The filtered result of the first slice in the real InSAR tensor acquired by the filter operating on the matrix model.
Figure 18. The filtered result of the first slice in the real InSAR tensor acquired by the filter operating on the matrix model.
Remotesensing 12 02582 g018
Figure 19. The filtered result of the second slice in the real InSAR tensor acquired by the filter operating on the matrix model.
Figure 19. The filtered result of the second slice in the real InSAR tensor acquired by the filter operating on the matrix model.
Remotesensing 12 02582 g019
Figure 20. The filtered result of the fifth slice in the real InSAR tensor acquired by the filter operating on the matrix model.
Figure 20. The filtered result of the fifth slice in the real InSAR tensor acquired by the filter operating on the matrix model.
Remotesensing 12 02582 g020aRemotesensing 12 02582 g020b
Figure 21. The filtered phase profiles with OLCP-InSAR and other reference filtered methods operating on the matrix model. (a) The filtered phase profiles at row 108 and column 180 to 230. (b) The filtered phase profiles at row 178 and column 45 to 95. (c) The filtered phase profiles at row 140 and column 275 to 325. (d) The filtered phase profiles at row 117 and column 300 to 350. (e) The filtered phase profiles at row 87 and column 165 to 215. (f) The filtered phase profiles at row 226 and column 80 to 130.
Figure 21. The filtered phase profiles with OLCP-InSAR and other reference filtered methods operating on the matrix model. (a) The filtered phase profiles at row 108 and column 180 to 230. (b) The filtered phase profiles at row 178 and column 45 to 95. (c) The filtered phase profiles at row 140 and column 275 to 325. (d) The filtered phase profiles at row 117 and column 300 to 350. (e) The filtered phase profiles at row 87 and column 165 to 215. (f) The filtered phase profiles at row 226 and column 80 to 130.
Remotesensing 12 02582 g021aRemotesensing 12 02582 g021b
Table 1. The mean square error (MSE) of filtered tensors with different noise conditions.
Table 1. The mean square error (MSE) of filtered tensors with different noise conditions.
Noise IntensityOutlier RatioHoRPCAWHoRPCAKBR-InSAROLCP-InSAR
3 dB30%0.250.200.170.04
5 dB30%0.190.140.090.03
7 dB30%0.160.110.070.03
5 dB20%0.150.080.050.03
5 dB40%0.370.290.220.04
Table 2. The number of remaining residues of filtered tensors with different noise conditions.
Table 2. The number of remaining residues of filtered tensors with different noise conditions.
Noise IntensityOutlier RatioHoRPCAWHoRPCAKBR-InSAROLCP-InSAR
3 dB30%69,525617,7501,039,525287,275
5 dB30%65,275539,300262,225226,875
7 dB30%61,575480,450307,625187,550
5 dB20%65,350453,175363,825222,750
5 dB40%67,850637,4001,268,250251,825
Table 3. MSE of a filtered interferogram with different noise conditions.
Table 3. MSE of a filtered interferogram with different noise conditions.
Noise IntensityOutlier RatioGoldsteinBoxcarNL-SARNL-InSARInSAR-BM3DOLCP-InSAR
(Matrix)
3 dB30%0.780.390.380.220.160.013
5 dB30%0.450.280.590.190.140.011
7 dB30%0.530.260.350.150.100.0086
5 dB20%0.270.130.100.0800.0500.010
5 dB40%0.601.00.590.350.280.012
Table 4. The number of residues remaining in a filtered interferogram with different noise conditions.
Table 4. The number of residues remaining in a filtered interferogram with different noise conditions.
Noise IntensityOutlier RatioGoldsteinBoxcarNL-SARNL-InSARInSAR-BM3DOLCP-InSAR
(Matrix)
3 dB30%43,72913,40621,578500220607595
5 dB30%25,230557427,697182856769
7 dB30%35,692905921,135341517595222
5 dB20%34,097961310,063392418376426
5 dB40%11,96137,55420,891462419066542
Table 5. Residues remaining in the filtered real InSAR tensor.
Table 5. Residues remaining in the filtered real InSAR tensor.
HoRPCAWHoRPCAKBR-InSAROLCP-InSAR
1776842,39356,49428,768
Table 6. Residues remaining in the filtered interferograms.
Table 6. Residues remaining in the filtered interferograms.
Slice ID BoxcarGoldsteinNL-InSARNL-SARInSAR-BM3DOLCP-InSAR
1#517810,8151531418837773063
2#31795727627342021791593
5#446775341267541138063371
Back to TopTop