Next Article in Journal
An Omnidirectional Morphological Method for Aerial Point Target Detection Based on Infrared Dual-Band Model
Next Article in Special Issue
A Variational Model for Sea Image Enhancement
Previous Article in Journal
A Double-Sampling Extension of the German National Forest Inventory for Design-Based Small Area Estimation on Forest District Levels
Previous Article in Special Issue
System Noise Removal for Gaofen-4 Area-Array Camera
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Millimeter-Wave InSAR Image Reconstruction Approach by Total Variation Regularized Matrix Completion

1
Purple Mountain Observatory, Chinese Academy of Sciences, Nanjing 210034, China
2
Key Lab of Radio Astronomy, Chinese Academy of Sciences, Nanjing 210034, China
*
Author to whom correspondence should be addressed.
Submission received: 17 May 2018 / Revised: 23 June 2018 / Accepted: 30 June 2018 / Published: 3 July 2018
(This article belongs to the Special Issue Data Restoration and Denoising of Remote Sensing Data)

Abstract

:
Millimeter-wave interferometric synthetic aperture radiometer (InSAR) can provide high-resolution observations for many applications by using small antennas to achieve very large synthetic aperture. However, reconstruction of a millimeter-wave InSAR image has been proven to be an ill-posed inverse problem that degrades the performance of InSAR imaging. In this paper, a novel millimeter-wave InSAR image reconstruction approach, referred to as InSAR-TVMC, by total variation (TV) regularized matrix completion (MC) in two-dimensional data space, is proposed. Based on the a priori knowledge that natural millimeter-wave images statistically hold the low-rank property, the proposed approach represents the object images as low-rank matrices and formulates the data acquisition of InSAR in two-dimensional data space directly to undersample visibility function samples. Subsequently, using the undersampled visibility function samples, the optimal solution of the InSAR image reconstruction problem is obtained by simultaneously adopting MC techniques and TV regularization. Experimental results on simulated and real millimeter-wave InSAR image data demonstrate the effectiveness and the significant improvement of the reconstruction performance of the proposed InSAR-TVMC approach over conventional and one-dimensional sparse InSAR image reconstruction approaches.

1. Introduction

Millimeter-wave interferometric synthetic aperture radiometer (InSAR) is a powerful observation system with high-resolution for many applications across the geographical and life sciences, including remote sensing, atmosphere monitoring, weather and climate forecast, anti-terrorist and security check [1,2,3,4]. It outperforms millimeter-wave real aperture radiometry mainly due to the advantages of high-resolution without very large and massive real aperture antennas and a large field of view without mechanical scanning. Sparsely arranging small antennas, millimeter-wave InSAR synthesizes a larger antenna aperture to achieve spatial high-resolution by acquiring cross-correlation value samples, namely the visibility function [5], between the signals received from pairs of these small antennas.
The imaging principle of millimeter-wave InSAR is based on the Fourier transform between the acquired visibility function and modified millimeter-wave brightness temperature distribution of the imaged scene; thus, FFT method is the basic approach to solve the millimeter-wave InSAR image reconstruction problem [1]. However, this reconstruction has been proven to be an ill-posed inverse problem, for which it is hard to obtain a unique and stable solution and which suffers from system noise and error interference. Furthermore, the high-resolution requirement of millimeter-wave InSAR greatly increases system and calibration complexity, data acquisition time, data processing and calculation amount. These difficulties in InSAR image reconstruction degrade performance and limit the application of InSAR in some scenarios.
Over the past few decades, iterative reconstruction methods have emerged as highly effective approaches to solve ill-posed inverse problems in imaging including denoising, deconvolution and interpolation [6,7,8]. These methods produce successful outcomes in millimeter-wave InSAR image reconstruction [9,10,11], but iterative reconstruction methods are still challenging to deploy in practice due to the high computational cost and time-consumption. These methods are limited by the Nyquist sampling theorem and have no ability to reduce data acquisition amount. In recent years, based on the a priori knowledge that millimeter-wave images are one-dimensional sparse after being formulated as vectors, compressive sensing (CS) techniques [12,13,14] have been successfully applied to InSAR image reconstruction, which significantly reduces the required data measurements and systematic complexity and maintains the advantages of high-resolution in the meantime [15,16,17,18]. Millimeter-wave InSAR reconstruction using CS techniques formulates the imaged scene as sparse vectors and uses sparse signal recovery methods to reconstruct the millimeter-wave brightness temperature distribution in a one-dimensional data space. However, the original data space of the InSAR imaged scene and final brightness temperature distribution output of InSAR are both two-dimensional. The performance of CS methods degrades while using one-dimensional signal recovery methods in a large-scale two-dimensional data space, ignoring the high-dimensional a priori knowledge and signal characteristics of millimeter-wave images.
More recently, matrix completion (MC) theory [19,20,21] has attracted much attention along with the rapid development of high-dimensional sparse representation. MC theory is of interest in cases for how to recover an unknown matrix when only a subset of entries of the unknown matrix is observed, and even these observed entries are corrupted by noise. Just as the l 0 norm is the definition of a vector’s sparsity in one-dimensional sparse representation, the rank norm has been proven to be the sparsity of the matrix in two-dimensional sparse representation. Then, for the above problem, if the unknown matrix holds the low-rank property and the observation satisfies some certain conditions, the MC problem can be solved exactly by solving a rank norm optimization problem based on the subset consisting of observed entries or corrupted entries. Owing to the advantages of MC theory, it has been adopted in many radar scenarios, like multiple-input and multiple-output (MIMO) radar [22] and synthetic aperture radar (SAR) [23]. In those radar scenarios, the low-rank condition is guaranteed by applying matched filter methods or transform domain methods on the receiving echo signal, which is certain and coherent along a large transmission distance.
Due to the important a priori knowledge that the natural millimeter-wave images statistically hold the low-rank property and the interferometric correlation operation of InSAR is a linear transformation operation, millimeter-wave InSAR images could be reconstructed based on a part of the visibility function samples via MC techniques. In this paper, a novel lower complexity millimeter-wave InSAR image reconstruction approach by total variation (TV) regularized matrix completion in two-dimensional data space is proposed, referred to as InSAR-TVMC. TV regularization is adopted here as it is a very powerful method in image processing to preserve edge information, explore the spatial piecewise smooth structure and further enhance the suppression ability of noise and error interference [24,25,26,27]. InSAR-TVMC solves the InSAR image reconstruction problem by simultaneously using MC techniques and TV regularization, based only on undersampled visibility function samples. Importantly, InSAR-TVMC achieves the advantages of InSAR using CS techniques directly in the matrix space by further using two-dimensional signal characteristics and a priori knowledge of millimeter-wave InSAR images. Experimental results show significant improvements of reconstruction performance over conventional and CS-based InSAR image reconstruction methods.
The rest of this paper is structured as follows. Section 2 describes the basic millimeter-wave InSAR imaging principle. The proposed InSAR-TVMC imaging approach is given in Section 3. Experimental results are presented in Section 4, and conclusions are given in Section 5.

2. Millimeter-Wave InSAR Imaging Principle

The millimeter-wave InSAR imaging system could obtain the complex cross-correlation value, namely visibility function samples, between small antennas receiving radiation from the imaged scene. The visibility function could be expressed as:
V ( u , v ) = ξ 2 + η 2 1 T ( ξ , η ) e j 2 π ( u ξ + v η ) d ξ d η
where ( ξ , η ) are the direction cosine coordinates, ( u , v ) are spatial frequency coordinates (baselines) and T M ( ξ , η ) is the modified brightness temperature defined as:
T M ( ξ , η ) = D c D l F c ( ξ , η ) F l * ( ξ , η ) 4 π 1 ξ 2 η 2 r u ξ + v η f 0 T B ( ξ , η )
where T B ( ξ , η ) is the brightness temperature; D c and D l are the maximum antenna directivities associated with the target; F c ( ξ , η ) and F l ( ξ , η ) are the normalized antenna voltage patterns of the elements in the array, which should be identical for all the elements in the ideal case; r c , l ( · ) is the fringe washing function accounting for the spatial decorrelation effects of the receiver’s frequency responses, which is nearly equal to one for narrow bandwidth receivers; f 0 is the center frequency.
In the ideal case, the visibility function is the Fourier transform of the modified brightness temperature image. Therefore, based on the antenna array, such as “T” shaped, “Y” shaped, circle shaped or just changing the baseline (relative position between two antennas), the InSAR imaging system could obtain V ( u , v ) and then reconstruct high-resolution target brightness temperature images T M ( ξ , η ) via algorithms like the fast Fourier transformation (FFT) approach or the G matrix inversion approach [28].
Considering errors in correlation observations and the receiving noise of InSAR, according to Equation (1), the basic signal model of InSAR could be expressed as:
V ( u , v ) = F [ T M ( ξ , η ) ] + e
where e represents the multiplicative errors and the receiving noise. The imaging method that directly takes the inverse Fourier transform on V ( u , v ) always needs extra strict calibration, or its performance degrades with the existence of e . To minimize the interference of error and noise in InSAR imaging, an alternative method replaces the Fourier transform relationship in Equation (3) by the G matrix model. However, the visibility function and InSAR images have to be formulated as vectors in the G matrix model. Based on the G matrix model, Equation (3) could be rewritten in matrix form:
V M × 1 = G M × N T N × 1 + E M × 1
where G M × N is the generalized impact function operator. Based on Equation (4), to obtain T N × 1 in the non-aliasing field of view, the lowest sampling frequency is required according to the Nyquist sampling theorem, which means N > 2 M . Therefore, the matrix Equation (4) is an underdetermined equation. Traditionally, an effective method is using the regularization method to reconstruct T N × 1 by using all the data of V M × 1 , considering the interference of error and noise. In this case, the raw two-dimensional InSAR data are formulated as a vector; consequently, the final brightness temperature distribution output is forced to be represented in the one-dimensional data space based on the G matrix model. This limits the use of some two-dimensional constraint conditions of the regularization method, which are suitable for two-dimensional images and could significantly improve image quality.
CS techniques utilize the one-dimensional sparsity of millimeter-wave InSAR images and consider the vector T N × 1 to be the l 0 norm sparse or sparse in a basis domain. Using the sparse representation and recovery methods, InSAR using CS techniques could maintain the advantages of high-resolution and meanwhile reduce the data measurements by only using part of the visibility function samples. A typical signal model of CS-based InSAR could be expressed as:
min T 1 s . t . V = Φ G Ψ T
where V is the part data of V M × 1 and T is the sparse vector of T N × 1 in the sparse basis matrix Ψ . However, CS-based InSAR is still limited by the one-dimensional imaging model. It formulates the two-dimensional InSAR images as vectors and completes the reconstruction in one-dimensional data space. The output of InSAR is resized to the image in the last step, which means CS-based InSAR could not use high-dimensional signal characteristics and a priori knowledge.

3. InSAR-TVMC Imaging Approach

This section provides a general overview of the MC theory and presents the signal model of the proposed InSAR-TVMC, then the specific recovery method is given.

3.1. Principles of MC

Consider the problem of recovering a low-rank matrix M C m × n based on partial entries observed, if the unknown matrix M holds the low-rank property and the observation is sufficiently random; M could be recovered by solving a rank norm optimization:
min r a n k ( X ) s . t . P Ω ( M ) = P Ω ( X )
where P Ω ( · ) is observation projection operator, Ω is the corresponding subset of the entries and X is the unknown variable. However, the observed results are always corrupted by noise in practice; in this case, the recovery of M is completed by solving the following optimization problem:
min r a n k ( X ) s . t . P Ω ( M ) P Ω ( X ) F 2 < ε
where ε is a parameter related to the observation noise variance σ 2 .
As the rank norm of a matrix is non-convex and discontinuous, the low-rank recoveries of Equations (6) and (7) are NP-hard problems, which are hard to solve.
In order to overcome this difficulty, a common method is to use the approximate convex relaxation. For example, in CS theory (one-dimensional vector sparse representation), the l 0 norm is convex and discontinuous, which is approximately relaxed to its most close convex functional l 1 norm to solve the optimization problem. In the two-dimensional matrix sparse representation theory, the nuclear norm (sum of all the singular values) is the nearest convex relaxation of the rank norm [29]. Therefore, the above low-rank recovery of the Equation (7) could be relaxed to the following nuclear norm optimization based on reasonable convex relaxation approximation:
min X * s . t . P Ω ( M ) P Ω ( X ) F 2 < ε
where · * denotes the nuclear norm. It has been proven that if the MC problem meets the condition of incoherence and the number of observed entries satisfies the probabilistic bound, the original low-rank recovery problem could be solved by nuclear norm relaxation with a high probability. Since the observation projection operator P Ω ( · ) could be simply treated as a linear mapping operator, Equation (8) could be rewritten as a nuclear norm optimization problem of the mapping matrix:
min X * s . t . A ( X ) b F 2 < ε
where X C m × n is the decision variable, b C p × q is the mapping result and A ( · ) : C m × n C p × q denotes the linear mapping operator.

3.2. Signal Model of InSAR-TVMC

The millimeter-wave brightness temperature image is a kind of natural image that has the low-rank property as a priori knowledge, and the data acquisition process of InSAR imaging is the mapping process from millimeter-wave brightness temperature image T C m × n to visibility function V C p × q based on the Fourier transform operator F ( · ) ; therefore, according to Equations (3) and (9), the signal model of millimeter-wave InSAR using MC techniques is proposed:
min T * s . t . F ( T ) V F 2 < ε
where millimeter-wave bright temperature image T C m × n is the decision variable, visibility function V C p × q is the mapping result and ε denotes the parameter related to the level of the observation error and the receiving noise of InSAR.
Using completely random selection of visibility function samples to realize random undersampling and based on the undersampled visibility function data, Equation (10) could be rewritten as:
min T * s . t . F ( T ) V F 2 < ε
where V is the undersampled visibility function data used in the recovery.
As the Fourier transform operator F ( · ) is linear, V and T are both two-dimensional data, Equation (11) covers convex formulations and the decision variable T is a matrix rather than a vector, InSAR imaging will be completed in the two-dimensional data space, which demands much lower computational and memory cost.
Similar to the traditional one-dimensional InSAR imaging method, the solution of Equation (11) could be obtained by regularization methods. The signal characteristics and a priori knowledge about the real solutions are always used as constraint conditions in regularization methods. Therefore, considering the nuclear norm as a kind of constraint condition, the regularization solution model of Equation (11) can be expressed as:
min T L ( T ) = F ( T ) V F 2 + λ T *
where the nuclear norm item T * is the constraint condition and λ is the regularization parameter.
In this two-dimensional regularization imaging model, we can also employ another kind of two-dimensional regularization constraint condition. In image processing, total variation (TV) norm minimization of the image is subjected to a constraint on image fidelity to observed data, which could help to preserve edge information, explore the spatial piecewise smooth structure and further enhance the suppression ability of error and noise interference. Here, we use the anisotropic TV norm [24] as an additional constraint to the nuclear norm, which is defined as:
T V ( T ) = i = 1 m 1 j = 1 n 1 | T i , j T i + 1 , j | + | T i , j T i , j + 1 | + i = 1 m 1 | T i , n T i + 1 , n | + j = 1 n 1 | T m , j T m , j + 1 |
Therefore, the TV norm-based InSAR regularization imaging model is expressed as:
min T L ( T ) = F ( T ) V F 2 + λ T T V
In this paper, we simultaneously use nuclear norm and TV norm regularization to establish a new functional constraint condition; the proposed InSAR-TVMC regularization model is defined as follows:
min T J ( T ) = 1 2 F ( T ) V F 2 + λ 1 T * + λ 2 T T V
The nuclear norm regularization considers the global information of the matrix sequence, while TV norm regularization encourages each matrix to be locally consistent. The proposed InSAR-TVMC regularization model (15) combines both types of a priori knowledge by exploiting two-dimensional sparsity and local signal characteristics to achieve more robust performance.

3.3. Recovery Method of InSAR-TVMC

In the traditional way, all the data of InSAR are processed in a one-dimensional vector data space based on the G matrix model. To solve the proposed InSAR-TVMC regularization problem directly in the two-dimensional matrix data space, a new two-dimensional imaging model of InSAR [30] is adopted in this paper.
Inspired by the generalized impact function operator used in the G matrix model, we use generalized impact function matrices D 1 and D 2 in a rectangular field of view formed by a rectangular antenna array configuration, such as “T” shaped, “U” shaped and “L” shaped. The pixel rules of the millimeter-wave bright image T N × N are distributed in the standard two-dimensional rectangular grid, which can be rewritten as the following two-dimensional matrix form:
V C × L = D 1 C × N T N × N D 2 N × L
where the visibility function V C × L is also distributed in the two-dimensional rectangular grid to form the matrix directly, its element V c × l is the visibility function sample value between the antennas (antenna#c and antenna#l) in different locations. It should be noted that for visibility function V distributed in hexagon and circular grids formed by “Y” shaped and circular shaped antenna array configurations, the new two-dimensional imaging model cannot be used for them directly.
The new two-dimensional imaging model (16) represents the Fourier transform operator F ( · ) as a two-dimensional linear transformation directly in matrix form; therefore, Equation (15) could be rewritten as:
min T J ( T ) = 1 2 D 1 · T · D 2 V F 2 + λ 1 T * + λ 2 T T V
Even though TV and nuclear norm minimization are both convex functions, the joint TV and nuclear norm minimization (17) is still very hard to solve as TV and the nuclear norm are non-separated and non-smooth. Inspired by the method [31] of solving TV regularization via its dual form, we solve a primal-dual form of the InSAR-TVMC regularization problem (17) instead of directly solving the original problem. Using the primal-dual form of the total variation norm by the Legendre–Fenchel transformation, the original problem (17) could be converted as:
min T max Y 1 2 D 1 · T · D 2 V F 2 + λ 1 T * + λ 2 ( < T , Y > ) I B ( Y )
where Y is the dual variable and ( · ) is the real part operator, while I B ( Y ) is the indicator function of the l unit norm ball defined as:
I B ( Y ) = 0 Y 1 + o t h e r w i s e
Then, the min-max problem (18) can be solved by a splitting scheme as:
T n + 1 = arg min 1 2 T T n F 2 + t 1 2 D 1 · T · D 2 V F 2 + t 1 λ 1 T * + t 1 λ 2 ( < T , Y n > )
Y n + 1 = arg min 1 2 Y Y n F 2 t 2 λ 2 ( < ( 2 T n + 1 T n ) , Y > ) + I B ( Y )
where T n and Y n are the primal and dual variables in the n-th iteration, respectively; t 1 and t 2 denote the corresponding iteration step sizes. To simplify Equation (20), a convex optimization using the gradient iterative solver is adopted by the InSAR-TVMC approach to approximate the least squares term [32]. Here, an auxiliary function f ( T ) is defined as:
f ( T ) = 1 2 D 1 · T · D 2 V F 2
Obviously, f ( T ) is convex and continuously differentiable with a smooth Lipschitz gradient, then we have:
f ( T ) = D 1 T ( D 1 · T · D 2 V ) D 2 T
Based on Equation (23), the Lipschitz constant is given by L = λ max ( D 1 T · D 1 · D 2 · D 2 T ) . By extending this gradient updating mechanism to Equation (20), we obtain:
T n + 1 = arg min 1 2 T T n F 2 + t 1 λ 1 T * + t 1 λ 2 ( < T , Y n > ) t 1 2 f ( T n ) + t 1 L 2 T T n F 2 + t 1 ( < f ( T n ) , T T n > )
By neglecting constant terms, Equation (24) can be written as:
T n + 1 = arg min 1 2 T ( T n t 1 1 + t 1 L f ( T n ) F 2 + t 1 λ 2 1 + t 1 L ( < T , Y n > ) + t 1 λ 1 1 + t 1 L T *
With the existence of T , the closed-form solution of Equation (25) is still unclear. To continue simplifying the problem to get a simple nuclear norm regularization form, an adjoint operator of the difference operator is introduced here. To describe the processing, we need to define the following items:
Define linear operator L T that maps matrix pairs ( p , q ) to an m-by-n matrix as:
( L T ( p , q ) ) i , j = p i , j p i 1 , j + q i , j q i , j 1 for 1 i m , 1 j n
where p C m 1 × n , q C m × n 1 and p 0 , j = p m , j = q i , 0 = q i , n .
The adjoint of the operator L T , denoted by L , is a linear map from an m-by-n matrix T to a matrix pair ( p , q ) and is defined by:
L ( T ) = ( p , q )
where:
p i , j = T i , j T i + 1 , j for 1 i m 1 , 1 j n q i , j = T i , j T i , j + 1 for 1 i m , 1 j n 1
Based on the definition of the adjoint of operator L and the difference operator, the linear operator L T satisfies:
< T , Y > = < T , L T ( Y ) >
Therefore, the dual variable Y is constructed by the matrix pair ( p , q ) , where:
( L T ( Y ) ) i , j = ( L T ( p , q ) ) i , j = p i , j p i 1 , j + q i , j q i , j 1
Based on Equations (29) and (30), we could simplify Equation (25) to the nuclear norm minimization problem:
T n + 1 = arg min 1 2 T T c n F 2 + t 1 λ 1 1 + t 1 L T *
where:
T c n = T n t 1 1 + t 1 L f ( T n ) t 1 1 + t 1 L L T ( Y n )
As the nuclear norm minimization is a convex optimization problem and has the closed-form solution, we could solve it by using the singular value thresholding (SVT) algorithm [33], which is simply expressed as:
S τ ( M ) = arg min X { 1 2 X M F 2 + τ X * }
where S τ ( · ) = d i a g ( max ( σ i τ , 0 ) 1 i r ) is the singular value thresholding operator. The SVT algorithm is a simple first-order iteration method, and the computational cost is low. Furthermore, all the data processing of the SVT algorithm is completed in matrix form, and the required storage space is minimal during every iteration. The specific SVT iteration algorithm is expressed as follows: for fixed parameter τ = t 1 λ 1 1 + t 1 L 0 , we have:
T n + 1 = S τ ( T c n )
Then, we consider the Y subproblem in Equation (21):
Y n + 1 = arg min 1 2 Y Y n F 2 t 2 λ 2 ( < ( 2 T n + 1 T n ) , Y > ) + I B ( Y )
After simplification, it becomes:
Y n + 1 = arg min 1 2 Y Y c n F 2 + I B ( Y )
where:
Y c n = Y n + t 2 λ 2 ( 2 T n + 1 T n )
The solution of Equation (36) can be obtained by the Euclidean projection of Y c n onto a l unit norm ball, which can be evaluated by:
Y n + 1 = sgn ( Y c n ) · min ( Y c n , 1 )
where sgn ( · ) is the sign operator.
Now, the T and Y subproblems have been solved. The proposed InSAR-TVMC approach is summarized in Table 1.

4. Experiments and Results

In this section, numerical simulation experiments are carried out on simulated and real millimeter-wave InSAR brightness temperature images to demonstrate the effectiveness and performance of the proposed InSAR-TVMC imaging approach, compared to the conventional FFT approach and CS approach. The InSAR imaging system uses a “T”-shaped antenna configuration with 127 antenna elements to form a rectangular visibility function distribution, shown in Figure 1. The specific InSAR system parameters are listed in Table 2. The fringe washing function in Equation (1) is set to one because the InSAR imaging system satisfies the narrow bandwidth system condition. The millimeter-wave InSAR brightness temperature images are simulated for each baseline of InSAR according to Equation (1) to form the measured visibility function samples, which are used to reconstruct InSAR images via the FFT approach, the CS approach (the regularization parameter λ = 0.085 ) and the InSAR-TVMC approach (the regularization parameters λ 1 = 0.9 and λ 2 = 0.01 ).
In the first experiments, the performances of all approaches are evaluated using a simulated millimeter-wave brightness temperature image of the Earth viewed from geostationary orbit, shown in Figure 2a, the brightness temperature of which is 250~350 K and the cold sky brightness temperature of which is 2.73 K. In the ideal case where the aforementioned systematic receiving noise and radiometric error are ignored, reconstructed images of Figure 2a by the three approaches are shown in Figure 2b–d, respectively. It should be pointed out that in this experiment, the FFT approach uses the entire visibility function samples to complete the reconstruction, while the CS approach and InSAR-TVMC approach only use 70% of the visibility function samples to realize undersampling based on one-dimensional and two-dimensional sparsity, respectively.
Comparing the reconstruction result in Figure 2, the FFT approach, the CS approach and the InSAR-TVMC approach all could obtain InSAR reconstruction results well ignoring system receiving noise and interferometric measurement error. However, InSAR-TVMC approach produces smoother results and weaker oscillation rings on the edge of the Earth and sky, showing the better suppression ability of the Gibbs effect.
In real interferometric measurements of the InSAR imaging system, the system receiving noise and interferometric measurement error are inevitably caused by hardware and position uncertainty. To test the sensitivity of the InSAR-TVMC approach and the compared approaches, the visibility function samples are corrupted by adding a white Gaussian noise with zero mean and variance σ 2 representing all the noise and error in real interferometric measurements. Using a corrupted visibility function samples with two different levels of noise intensity ( σ 2 = 0.05 and σ 2 = 0.1 ), reconstructed images of Figure 2a via FFT using entire samples, CS and InSAR-TVMC using 70% samples are shown in Figure 3.
It is clear that the reconstructed image in Figure 3a,b via FFT was seriously destroyed by noise, and it is almost impossible to distinguish any detailed information of the Earth. Even on the edge of the Earth and sky, oscillation rings are expanded. That is because FFT itself does not contain any denoising processing and suppression ability on the Gibbs effect. The quality of the reconstructed image in Figure 3c via CS is much better than FFT under a low level of noise intensity ( σ 2 = 0.05 ), but CS also suffers from the noise and error interference, so the performance of the reconstructed images in Figure 3d degrades. Compared to the reconstructed images via FFT and CS, InSAR-TVMC still produces visually better results, which retain much more detailed information of the Earth and exhibit better suppression of oscillation rings caused by the Gibbs effect, demonstrating the effectiveness of the TV norm and nuclear norm constraints.
To quantitatively analyze the sensitivity of the above three methods, the peak signal to noise ratio (PSNR) performance with different noise levels is shown in Figure 4.
Figure 4 exhibits the robustness of InSAR-TVMC to noise and error interference compared to FFT and CS objectively, demonstrating that InSAR-TVMC could better improve the accuracy of InSAR image reconstruction.
To further analyze the robustness of InSAR-TVMC to the undersampling ratio (usage rate of visibility function samples), comparison experiments of CS and InSAR-TVMC were carried out with different usage rates of visibility function samples corrupted by white Gaussian noise with variance σ 2 = 0.05 . Recovery results of a 40%, 50% and 60% undersampling ratio of CS and InSAR-TVMC are shown in Figure 5.
Figure 5a shows that much detailed information is missing, and the outline information of the Earth is distorted in the reconstructed image by CS using 40% of samples, demonstrating that the reconstructed images via CS is seriously distorted with a low undersampling ratio. From Figure 5d–f, we can find that InSAR-TVMC produces a much more robust result with different undersampling ratios.
Moreover, quantitative evaluation of the robustness to different undersampling ratios using the root-mean square error (RMSE) performance of InSAR-TVMC and CS is presented in Figure 6. From Figure 6, we can see that InSAR-TVMC has a lower RMSE at any undersampling ratio, especially having significant advantages over CS at low undersampling rates, demonstrating the robustness of InSAR-TVMC objectively.
After the comparison experiments on simulated millimeter-wave brightness temperature data of the Earth, we explored the InSAR-MC approach on real millimeter-wave InSAR data, which were acquired by a Geostationary Interferometric Microwave Sounder (GIMS) demonstrator in the near field by the National Space Science Center of China [34]. The test radiometric image data are the 50.3-GHz brightness temperature data shown in Figure 7a, and the light image data are in Figure 7b. Reconstructed images via CS and InSAR-TVMC with 70% of visibility function samples corrupted by white Gaussian noise ( σ 2 = 0.05 ) are presented in Figure 7c,d, respectively.
Figure 7 shows that InSAR-TVMC produces a much better reconstructed result containing clear details, while the result of CS lost much detailed information compared to the original radiometric image data, demonstrating the effectiveness of the InSAR-TVMC approach.
What is more, as the InSAR-TVMC approach completes the reconstruction in the two-dimensional data space, so it demands much lower computational and memory cost, while the dimensionality of the vector converted from the matrix data by the CS approach could be quite large; especially, the scale of the original matrix data themselves is very large, as well. Table 3 reports the MATLAB runtime (MATLAB2016a with a 3.2-GHz Intel i7-8700 processor, 32 GB of memory) of the InSAR-TVMC and CS approach with different undersampling ratios. The result demonstrates that the two-dimensional sparse imaging approach, InSAR-TVMC, is more suitable for InSAR image reconstruction compared to one-dimensional sparse imaging approaches.

5. Conclusions

In this paper, we proposed a novel millimeter-wave InSAR image reconstruction approach by total variation regularized matrix completion for high-resolution imaging with undersampled data. Based on the a priori knowledge that millimeter-wave InSAR images hold the low-rank property, the proposed InSAR-TVMC approach represented the object images as low-rank matrices and formulated the data acquisition of InSAR in the two-dimensional data space directly to undersample visibility function samples. Then, the optimal solution of the InSAR image reconstruction problem was obtained using the undersampled visibility samples by simultaneously solving total variation and nuclear norm regularization via convex optimization. Using corrupted visibility function samples with high levels of noise intensity ( σ 2 = 0.1 ), the PSNR of the InSAR-TVMC reconstructed result using 70% of samples is 25.8 dB compared to 22.3 dB for the CS result using 70% of samples and 17.1 dB for the FFT result using all of the samples. The reconstruction time of InSAR-TVMC could be shortened to 80.33 s compared to 141.58 s for CS, both with a 70% undersampling ratio. Experimental results demonstrate the effectiveness and the significant improvement of the reconstruction performance of the proposed InSAR-TVMC approach over conventional and one-dimensional sparse InSAR image reconstruction approaches.

Author Contributions

Y.Z. and S.S. conceived of and designed the experiments. Y.Z. performed the experiments. All authors analyzed the data and contributed to writing the paper.

Funding

This work was supported in part by the National Key Basic Research and Development Program 2017YFA0304003 and 2018YFA0404701, in part by the Chinese Academy of Sciences (CAS) under Grant Nos. QYZDJ-SSW-SLH043 and GJJSTD20180003, in part by the National Postdoctoral Program for Innovative Talents of China under Grant BX201700309 and in part by the China Postdoctoral Science Foundation under Grant 2017LH049.

Acknowledgments

The authors would like to thank the anonymous reviewer and editors for their helpful comments and suggestions

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Ruf, C.S.; Swift, C.T.; Tanner, A.B.; Le Vine, D.M. Interferometric synthetic aperture microwave radiometry for the remote sensing of the Earth. IEEE Trans. Geosci. Remote Sens. 1988, 26, 597–611. [Google Scholar] [CrossRef]
  2. Kerr, Y.H.; Waldteufel, P.; Wigneron, J.P.; Martinuzzi, J.; Font, J.; Berger, M. Soil moisture retrieval from space: The Soil Moisture and Ocean Salinity (SMOS) mission. IEEE Trans. Geosci. Remote Sens. 2002, 39, 1729–1735. [Google Scholar] [CrossRef]
  3. Reising, S.C.; Gaier, T.C.; Kummerow, C.D.; Padmanabhan, S.; Lim, B.H.; Brown, S.T.; Heneghan, C.; Chandra, C.V.; Olson, J.; Berg, W. Temporal Experiment for Storms and Tropical Systems Technology Demonstration (TEMPEST-D): Reducing risk for 6U-Class nanosatellite constellations. In Proceedings of the IEEE International Geoscience and Remote Sensing Symposium (IGARSS), Beijing, China, 10–15 July 2016; pp. 5559–5560. [Google Scholar] [CrossRef]
  4. Chen, H.M.; Lee, S.R.; Rao, M.; Slamani, M.A.; Varshney, P.K. Imaging for concealed weapon detection: A tutorial overview of development in imaging sensors and processing. IEEE Signal Proc. Mag. 2005, 22, 52–61. [Google Scholar] [CrossRef]
  5. Corbella, I.; Duffo, N.; Vall-llossera, M.; Camps, A.; Torres, F. The visibility function in interferometric aperture synthesis radiometry. IEEE Trans. Geosci. Remote Sens. 2004, 42, 1677–1682. [Google Scholar] [CrossRef] [Green Version]
  6. Chang, S.G.; Yu, B.; Vetterli, M. Adaptive wavelet thresholding for image denoising and compression. IEEE Trans. Image Process. 2000, 9, 1532–1546. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  7. Mignotte, M. A segmentation-based regularization term for image deconvolution. IEEE Trans. Image Process. 2006, 15, 1973–1984. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  8. Ramani, S.; Thevenaz, P.; Unser, M. Regularized interpolation for noisy images. IEEE Trans. Med. Imag. 2010, 29, 543–558. [Google Scholar] [CrossRef] [PubMed]
  9. Picard, B.; Anterrieu, E. Comparison of regularized inversion methods in synthetic aperture imaging radiometry. IEEE Trans. Geosci. Remote Sens. 2005, 43, 218–224. [Google Scholar] [CrossRef]
  10. Wu, L.; Hu, F.; He, F.; Li, J.; Peng, X.; Zhu, D. Statistical regularization in synthetic aperture imaging radiometry. In Proceedings of the IEEE International Geoscience and Remote Sensing Symposium (IGARSS), Milan, Italy, 26–31 July 2015; pp. 4726–4729. [Google Scholar] [CrossRef]
  11. Wu, L.; Hu, F.; He, F.; Li, J. Bayesian inference for inversion in synthetic aperture imaging radiometry. IEEE Geosci. Remote Sens. Lett. 2016, 13, 1049–1053. [Google Scholar] [CrossRef]
  12. Candes, E.J.; Romberg, J.; Tao, T. Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information. IEEE Trans. Inform. Theory 2006, 52, 489–509. [Google Scholar] [CrossRef]
  13. Candes, E.J.; Tao, T. Near-optimal signal recovery from random projections: Universal encoding strategies. IEEE Trans. Inf. Theory 2006, 52, 5406–5425. [Google Scholar] [CrossRef]
  14. Candes, E.J.; Wakin, M.B. An introduction to compressive sampling. IEEE Signal Process. Mag. 2008, 25, 21–30. [Google Scholar] [CrossRef]
  15. Li, S.Y.; Zhou, X.; Ren, B.; Sun, H.J.; Lv, X. A compressive sensing approach for synthetic aperture imaging radiometers. Prog. Electromagn. Res. 2013, 135, 583–599. [Google Scholar] [CrossRef]
  16. Chen, J.F.; Li, Y.H. The CS-Based imaging algorithm for near-Field synthetic aperture imaging radiometer. Prog. Electromagn. Res. 2014, 97, 911–914. [Google Scholar] [CrossRef]
  17. Zhang, Y.L.; Li, Y.H.; Zhu, S.J.; Li, Y. A robust reweighted L1-minimization imaging algorithm for passive millimeter-wave sair in near field. Sensors 2015, 15, 24945–24960. [Google Scholar] [CrossRef] [PubMed]
  18. Kpre, E.L.; Decroze, C. Passive Coding Technique Applied to Synthetic Aperture Interferometric Radiometer. IEEE Geosci. Remote Sens. Lett. 2017, 14, 1193–1197. [Google Scholar] [CrossRef]
  19. Candes, E.J.; Recht, B. Exact matrix completion via convex optimization. Found. Comput. Math. 2008, 9, 717–772. [Google Scholar] [CrossRef]
  20. Candes, E.J.; Plan, Y. Matrix completion with noise. Proc. IEEE 2010, 98, 925–936. [Google Scholar] [CrossRef]
  21. Candes, E.J.; Tao, T. The power of convex relaxation: Near-optimal matrix completion. IEEE Trans. Inf. Theory 2010, 56, 2053–2080. [Google Scholar] [CrossRef]
  22. Sun, S.Q.; Bajwa, W.U.; Petropulu, A.P. MIMO-MC radar: A MIMO radar approach based on matrix completion. IEEE Trans. Aerosp. Electron. Syst. 2015, 51, 1839–1852. [Google Scholar] [CrossRef] [Green Version]
  23. Yang, D.; Liao, G.S.; Zhu, S.Q.; Yang, X.; Zhang, X. SAR imaging with undersampled data via matrix completion. IEEE Geosci. Remote Sens. Lett. 2014, 14, 1539–1543. [Google Scholar] [CrossRef]
  24. Chambolle, A. An algorithm for total variation minimization and applications. J. Math. Imaging Vis. 2004, 20, 89–97. [Google Scholar] [CrossRef]
  25. Sidky, E.Y.; Pan, X. Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization. Phys. Med. Biol. 2008, 53, 4777–4807. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  26. Wang, Y.; Peng, J.; Zhao, Q.; Leung, Y.; Zhao, X.L.; Meng, D. Hyperspectral image restoration via total variation regularized low-rank tensor decomposition. IEEE J-STARS. 2017, 99, 1–19. [Google Scholar] [CrossRef]
  27. Yao, J.W.; Xu, Z.; Huang, X.L.; Huang, J. An efficient algorithm for dynamic MRI using low-rank and total variation regularizations. Med. Image Anal. 2018, 44, 14–27. [Google Scholar] [CrossRef] [PubMed]
  28. Anterrieu, E. A resolving matrix approach for synthetic aperture imaging radiometers. IEEE Trans. Geosci. Remote Sens. 2004, 42, 1649–1656. [Google Scholar] [CrossRef]
  29. Fazel, M.; Pong, T.K.; Sun, D.; Tseng, P. Hankel matrix rank minimization with applications to system identification and realization. SIAM J. Matrix Anal. Appl. 2013, 34, 946–977. [Google Scholar] [CrossRef]
  30. Chen, J.F.; Li, Y.H.; Zhang, Y.L. An accurate imaging algorithm for millimeter-wave synthetic aperture imaging radiometer in near-field. Prog. Electromagn. Res. 2013, 141, 517–535. [Google Scholar] [CrossRef]
  31. Chambolle, A.; Pock, T. On the ergodic convergence rates of a first-order primal-dual algorithm. Math. Program. 2016, 159, 253–287. [Google Scholar] [CrossRef]
  32. Beck, A.; Teboulle, M. A fast iterative shrinkage-thresholding algorithm for linear inverse problems. SIAM J. Imaging Sci. 2009, 2, 183–202. [Google Scholar] [CrossRef]
  33. Cai, J.F.; Candes, E.J.; Shen, Z. A Singular Value Thresholding Algorithm for Matrix Completion. SIAM J. Optimiz. 2008, 20, 1956–1982. [Google Scholar] [CrossRef]
  34. Zhang, C.; Liu, H.; Wu, J.; Zhang, S.; Yan, J.; Niu, L.; Sun, W.; Li, H. Imaging analysis and first results of the geostationary interferometric microwave sounder demonstrator. IEEE Trans. Geosci. Remote Sens. 2015, 53, 207–218. [Google Scholar] [CrossRef]
Figure 1. “T”-shaped antenna array and rectangular visibility function distribution.
Figure 1. “T”-shaped antenna array and rectangular visibility function distribution.
Remotesensing 10 01053 g001
Figure 2. Result of the reconstruction without noise interference. (a) Original simulated millimeter-wave InSAR image. (b) Reconstructed by FFT. (c) Reconstructed by compressive sensing (CS). (d) Reconstructed by InSAR-TVMC.
Figure 2. Result of the reconstruction without noise interference. (a) Original simulated millimeter-wave InSAR image. (b) Reconstructed by FFT. (c) Reconstructed by compressive sensing (CS). (d) Reconstructed by InSAR-TVMC.
Remotesensing 10 01053 g002
Figure 3. Result of the reconstruction. (a) Reconstructed by FFT, σ 2 = 0.05 . (b) Reconstructed by FFT, σ 2 = 0.1 . (c) Reconstructed by CS, σ 2 = 0.05 . (d) Reconstructed by CS, σ 2 = 0.1 . (e) Reconstructed by InSAR-TVMC, σ 2 = 0.05 . (f) Reconstructed by InSAR-TVMC, σ 2 = 0.1 .
Figure 3. Result of the reconstruction. (a) Reconstructed by FFT, σ 2 = 0.05 . (b) Reconstructed by FFT, σ 2 = 0.1 . (c) Reconstructed by CS, σ 2 = 0.05 . (d) Reconstructed by CS, σ 2 = 0.1 . (e) Reconstructed by InSAR-TVMC, σ 2 = 0.05 . (f) Reconstructed by InSAR-TVMC, σ 2 = 0.1 .
Remotesensing 10 01053 g003
Figure 4. PSNR performance with different noise levels of FFT, CS and InSAR-TVMC.
Figure 4. PSNR performance with different noise levels of FFT, CS and InSAR-TVMC.
Remotesensing 10 01053 g004
Figure 5. Result of the reconstruction at different undersampling ratios. (a) Reconstructed by CS at a 40% ratio. (b) Reconstructed by CS at a 50% ratio. (c) Reconstructed by CS at a 60% ratio. (d) Reconstructed by InSAR-TVMC at a 40% ratio. (e) Reconstructed by InSAR-TVMC at a 50% ratio. (f) Reconstructed by InSAR-TVMC at a 60% ratio.
Figure 5. Result of the reconstruction at different undersampling ratios. (a) Reconstructed by CS at a 40% ratio. (b) Reconstructed by CS at a 50% ratio. (c) Reconstructed by CS at a 60% ratio. (d) Reconstructed by InSAR-TVMC at a 40% ratio. (e) Reconstructed by InSAR-TVMC at a 50% ratio. (f) Reconstructed by InSAR-TVMC at a 60% ratio.
Remotesensing 10 01053 g005
Figure 6. RMSE performance with different undersampling ratios of CS and InSAR-TVMC.
Figure 6. RMSE performance with different undersampling ratios of CS and InSAR-TVMC.
Remotesensing 10 01053 g006
Figure 7. Result of the reconstruction. (a) Original real millimeter-wave InSAR image. (b) Optical image. (c) Reconstructed by CS. (d) Reconstructed by InSAR-TVMC.
Figure 7. Result of the reconstruction. (a) Original real millimeter-wave InSAR image. (b) Optical image. (c) Reconstructed by CS. (d) Reconstructed by InSAR-TVMC.
Remotesensing 10 01053 g007
Table 1. InSAR-total variation regularized matrix completion (TVMC) imaging approach.
Table 1. InSAR-total variation regularized matrix completion (TVMC) imaging approach.
Input: D1, D2, V, λ 1 , λ 2 , N
Output: brightness temperature T
Algorithm:
Step 1: initialize Y 0 , T 0 , and undersample V to obtain V ;
Step 2: update T n by using Equations (31) and (34);
Step 3: update Y n by using Equations (37) and (38);
Step 4: n + 1 , if n < N , repeat from Step 2; otherwise, go to Step 5;
Step 5: output T N
Table 2. Simulation InSAR system parameters.
Table 2. Simulation InSAR system parameters.
Simulation ParametersValue
center frequency50.3 GHz
bandwidth200 MHz
image pixel size128 × 128
image gray(0, 255)
brightness temperature2.73~350 K
antenna array (“T” shaped)64 + 63
visibility function samples size128 × 128
G matrix size16,384 × 1
D matrix size128 × 128
Table 3. Computation time comparison of CS and InSAR-TVMC.
Table 3. Computation time comparison of CS and InSAR-TVMC.
Undersampling Rate40%50%60%70%80%90%
t C S ( s ) 65.8781.86119.34141.58166.53180.45
t T V M C ( s ) 33.9146.7665.3180.3394.22103.74

Share and Cite

MDPI and ACS Style

Zhang, Y.; Miao, W.; Lin, Z.; Gao, H.; Shi, S. Millimeter-Wave InSAR Image Reconstruction Approach by Total Variation Regularized Matrix Completion. Remote Sens. 2018, 10, 1053. https://0-doi-org.brum.beds.ac.uk/10.3390/rs10071053

AMA Style

Zhang Y, Miao W, Lin Z, Gao H, Shi S. Millimeter-Wave InSAR Image Reconstruction Approach by Total Variation Regularized Matrix Completion. Remote Sensing. 2018; 10(7):1053. https://0-doi-org.brum.beds.ac.uk/10.3390/rs10071053

Chicago/Turabian Style

Zhang, Yilong, Wei Miao, Zhenhui Lin, Hao Gao, and Shengcai Shi. 2018. "Millimeter-Wave InSAR Image Reconstruction Approach by Total Variation Regularized Matrix Completion" Remote Sensing 10, no. 7: 1053. https://0-doi-org.brum.beds.ac.uk/10.3390/rs10071053

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop