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

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)=\underset{{\xi }^2+{\eta }^2\le 1}{\int \int }T(\xi ,\eta )e^{-j2\pi (u\xi +v\eta )}d\xi d\eta$$

(1)

where $(\xi ,\eta )$ are the direction cosine coordinates, $(u,v)$ are spatial frequency coordinates (baselines) and ${\mathbf{T}}_M(\xi ,\eta )$ is the modified brightness temperature defined as:

$${\mathbf{T}}_M(\xi ,\eta )=\frac{\sqrt{D_c{D}_l}{F}_c(\xi ,\eta )F_l^{*}(\xi ,\eta )}{4\pi \sqrt{1-{\xi }^2-{\eta }^2}}r\left(-\frac{u\xi +v\eta }{f_0}\right){\mathbf{T}}_B(\xi ,\eta )$$

(2)

where ${\mathbf{T}}_B(\xi ,\eta )$ is the brightness temperature; $D_c$ and $D_l$ are the maximum antenna directivities associated with the target; $F_c(\xi ,\eta )$ and $F_l(\xi ,\eta )$ 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 $\mathbf{V}(u,v)$ and then reconstruct high-resolution target brightness temperature images ${\mathbf{T}}_M(\xi ,\eta )$ 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:

$$\mathbf{V}(u,v)=\mathcal{F}\left[{\mathbf{T}}_M(\xi ,\eta )\right]+\mathbf{e}$$

(3)

where $\mathbf{e}$ represents the multiplicative errors and the receiving noise. The imaging method that directly takes the inverse Fourier transform on $\mathbf{V}(u,v)$ always needs extra strict calibration, or its performance degrades with the existence of $\mathbf{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:

$${\mathbf{V}}_{M\times 1}={\mathbf{G}}_{M\times N}{\mathbf{T}}_{N\times 1}+{\mathbf{E}}_{M\times 1}$$

(4)

where ${\mathbf{G}}_{M\times N}$ is the generalized impact function operator. Based on Equation (4), to obtain ${\mathbf{T}}_{N\times 1}$ in the non-aliasing field of view, the lowest sampling frequency is required according to the Nyquist sampling theorem, which means $N>2M$. Therefore, the matrix Equation (4) is an underdetermined equation. Traditionally, an effective method is using the regularization method to reconstruct ${\mathbf{T}}_{N\times 1}$ by using all the data of ${\mathbf{V}}_{M\times 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 ${\mathbf{T}}_{N\times 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{∥{\mathbf{T}}^{\prime }∥}_1\mathrm{s}.\mathrm{t}.\stackrel{\wedge }{\mathbf{V}}=\stackrel{\wedge }{\mathsf{\Phi }}{\mathbf{G}\mathsf{\Psi }\mathbf{T}}^{\prime }$$

(5)

where $\stackrel{\wedge }{\mathbf{V}}$ is the part data of ${\mathbf{V}}_{M\times 1}$ and ${\mathbf{T}}^{\prime }$ is the sparse vector of ${\mathbf{T}}_{N\times 1}$ in the sparse basis matrix $\mathsf{\Psi }$. 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 $\mathbf{M}\in {\mathbb{C}}^{m\times n}$ based on partial entries observed, if the unknown matrix $\mathbf{M}$ holds the low-rank property and the observation is sufficiently random; $\mathbf{M}$ could be recovered by solving a rank norm optimization:

$$minrank\left(\mathbf{X}\right)\mathrm{s}.\mathrm{t}.P_{\mathsf{\Omega }}\left(\mathbf{M}\right)=P_{\mathsf{\Omega }}\left(\mathbf{X}\right)$$

(6)

where $P_{\mathsf{\Omega }}(·)$ is observation projection operator, $\mathsf{\Omega }$ is the corresponding subset of the entries and $\mathbf{X}$ is the unknown variable. However, the observed results are always corrupted by noise in practice; in this case, the recovery of $\mathbf{M}$ is completed by solving the following optimization problem:

$$minrank\left(\mathbf{X}\right)\mathrm{s}.\mathrm{t}.{∥P_{\mathsf{\Omega }}\left(\mathbf{M}\right)-P_{\mathsf{\Omega }}\left(\mathbf{X}\right)∥}_F^2<\epsilon$$

(7)

where $\epsilon$ is a parameter related to the observation noise variance ${\sigma }^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{∥\mathbf{X}∥}_{*}\mathrm{s}.\mathrm{t}.{∥P_{\mathsf{\Omega }}\left(\mathbf{M}\right)-P_{\mathsf{\Omega }}\left(\mathbf{X}\right)∥}_F^2<\epsilon$$

(8)

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_{\mathsf{\Omega }}(·)$ 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{∥\mathbf{X}∥}_{*}\mathrm{s}.\mathrm{t}.{∥\mathcal{A}\left(\mathbf{X}\right)-\mathbf{b}∥}_F^2<\epsilon$$

(9)

where $\mathbf{X}\in {\mathbb{C}}^{m\times n}$ is the decision variable, $\mathbf{b}\in {\mathbb{C}}^{p\times q}$ is the mapping result and $\mathcal{A}(·):{\mathbb{C}}^{m\times n}\to {\mathbb{C}}^{{}^{p\times 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 $\mathbf{T}\in {\mathbb{C}}^{m\times n}$ to visibility function $\mathbf{V}\in {\mathbb{C}}^{p\times q}$ based on the Fourier transform operator $\mathcal{F}(·)$; therefore, according to Equations (3) and (9), the signal model of millimeter-wave InSAR using MC techniques is proposed:

$$min{∥\mathbf{T}∥}_{*}\mathrm{s}.\mathrm{t}.{∥\mathcal{F}\left(\mathbf{T}\right)-\mathbf{V}∥}_F^2<\epsilon$$

(10)

where millimeter-wave bright temperature image $\mathbf{T}\in {\mathbb{C}}^{m\times n}$ is the decision variable, visibility function $\mathbf{V}\in {\mathbb{C}}^{p\times q}$ is the mapping result and $\epsilon$ 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{∥\mathbf{T}∥}_{*}\mathrm{s}.\mathrm{t}.{∥\mathcal{F}\left(\mathbf{T}\right)-\stackrel{\wedge }{\mathbf{V}}∥}_F^2<\epsilon$$

(11)

where $\stackrel{\wedge }{\mathbf{V}}$ is the undersampled visibility function data used in the recovery.

As the Fourier transform operator $\mathcal{F}(·)$ is linear, $\stackrel{\wedge }{\mathbf{V}}$ and $\mathbf{T}$ are both two-dimensional data, Equation (11) covers convex formulations and the decision variable $\mathbf{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:

$$\underset{\mathbf{T}}{min}L\left(\mathbf{T}\right)={∥\mathcal{F}\left(\mathbf{T}\right)-\stackrel{\wedge }{\mathbf{V}}∥}_F^2+\lambda {∥\mathbf{T}∥}_{*}$$

(12)

where the nuclear norm item ${∥\mathbf{T}∥}_{*}$ is the constraint condition and $\lambda$ 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:

$$\begin{array}{cc}\hfill TV\left(\mathbf{T}\right)=& {\displaystyle \sum _{i=1}^{m-1}}{\displaystyle \sum _{j=1}^{n-1}}\left\{|T_{i,j}-T_{i+1,j}|+|T_{i,j}-T_{i,j+1}|\right\}\hfill \\ & +{\displaystyle \sum _{i=1}^{m-1}}|T_{i,n}-T_{i+1,n}|+{\displaystyle \sum _{j=1}^{n-1}}|T_{m,j}-T_{m,j+1}|\hfill \end{array}$$

(13)

Therefore, the TV norm-based InSAR regularization imaging model is expressed as:

$$\underset{\mathbf{T}}{min}L\left(\mathbf{T}\right)={∥\mathcal{F}\left(\mathbf{T}\right)-\stackrel{\wedge }{\mathbf{V}}∥}_F^2+\lambda {∥\mathbf{T}∥}_{TV}$$

(14)

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:

$$\underset{\mathbf{T}}{min}J\left(\mathbf{T}\right)=\frac{1}{2}{∥\mathcal{F}\left(\mathbf{T}\right)-\stackrel{\wedge }{\mathbf{V}}∥}_F^2+{\lambda }_1{∥\mathbf{T}∥}_{*}+{\lambda }_2{∥\mathbf{T}∥}_{TV}$$

(15)

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 $\mathbf{D}\mathbf{1}$ and $\mathbf{D}\mathbf{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 ${\mathbf{T}}_{N\times N}$ are distributed in the standard two-dimensional rectangular grid, which can be rewritten as the following two-dimensional matrix form:

$${\mathbf{V}}_{C\times L}=\mathbf{D}{\mathbf{1}}_{C\times N}{\mathbf{T}}_{N\times N}\mathbf{D}{\mathbf{2}}_{N\times L}$$

(16)

where the visibility function ${\mathbf{V}}_{C\times L}$ is also distributed in the two-dimensional rectangular grid to form the matrix directly, its element $V_{c\times 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 $\mathbf{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 $\mathcal{F}(·)$ as a two-dimensional linear transformation directly in matrix form; therefore, Equation (15) could be rewritten as:

$$\underset{\mathbf{T}}{min}J\left(\mathbf{T}\right)=\frac{1}{2}{∥\mathbf{D}\mathbf{1}·\mathbf{T}·\mathbf{D}\mathbf{2}-\stackrel{\wedge }{\mathbf{V}}∥}_F^2+{\lambda }_1{∥\mathbf{T}∥}_{*}+{\lambda }_2{∥\mathbf{T}∥}_{TV}$$

(17)

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:

$$\underset{\mathbf{T}}{min}\underset{\mathbf{Y}}{max}\frac{1}{2}{∥\mathbf{D}\mathbf{1}·\mathbf{T}·\mathbf{D}\mathbf{2}-\stackrel{\wedge }{\mathbf{V}}∥}_F^2+{\lambda }_1{∥\mathbf{T}∥}_{*}+{\lambda }_2\Re (<\nabla \mathbf{T},\mathbf{Y}>)-I_{B\infty }\left(\mathbf{Y}\right)$$

(18)

where $\mathbf{Y}$ is the dual variable and $\Re (·)$ is the real part operator, while $I_{B\infty }\left(\mathbf{Y}\right)$ is the indicator function of the $l_{\infty }$ unit norm ball defined as:

$$I_{B\infty }\left(\mathbf{Y}\right)=\left\{\begin{array}{cc}\hfill 0& ∥\mathbf{Y}∥\le 1\hfill \\ \hfill +\infty & otherwise\hfill \end{array}\right.$$

(19)

Then, the min-max problem (18) can be solved by a splitting scheme as:

$$\begin{array}{cc}\hfill {\mathbf{T}}^{n+1}=& argmin\frac{1}{2}{∥\mathbf{T}-{\mathbf{T}}^n∥}_F^2+\frac{t_1}{2}{∥\mathbf{D}\mathbf{1}·\mathbf{T}·\mathbf{D}\mathbf{2}-\stackrel{\wedge }{\mathbf{V}}∥}_F^2\hfill \\ & +t_1{\lambda }_1{∥\mathbf{T}∥}_{*}+t_1{\lambda }_2\Re (<\nabla \mathbf{T},{\mathbf{Y}}^n>)\hfill \end{array}$$

(20)

$${\mathbf{Y}}^{n+1}=argmin\frac{1}{2}{∥\mathbf{Y}-{\mathbf{Y}}^n∥}_F^2-t_2{\lambda }_2\Re (<\nabla (2{\mathbf{T}}^{n+1}-{\mathbf{T}}^n),\mathbf{Y}>)+I_{B\infty }\left(\mathbf{Y}\right)$$

(21)

where ${\mathbf{T}}^n$ and ${\mathbf{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\left(\mathbf{T}\right)$ is defined as:

$$f\left(\mathbf{T}\right)=\frac{1}{2}{∥\mathbf{D}\mathbf{1}·\mathbf{T}·\mathbf{D}\mathbf{2}-\stackrel{\wedge }{\mathbf{V}}∥}_F^2$$

(22)

Obviously, $f\left(\mathbf{T}\right)$ is convex and continuously differentiable with a smooth Lipschitz gradient, then we have:

$$\nabla f\left(\mathbf{T}\right)=\mathbf{D}{\mathbf{1}}^T(\mathbf{D}\mathbf{1}·\mathbf{T}·\mathbf{D}\mathbf{2}-\stackrel{\wedge }{\mathbf{V}})\mathbf{D}{\mathbf{2}}^T$$

(23)

Based on Equation (23), the Lipschitz constant is given by $L={\lambda }_{max}(\mathbf{D}{\mathbf{1}}^T·\mathbf{D}\mathbf{1}·\mathbf{D}\mathbf{2}·\mathbf{D}{\mathbf{2}}^T)$. By extending this gradient updating mechanism to Equation (20), we obtain:

$$\begin{array}{cc}\hfill {\mathbf{T}}^{n+1}=& argmin\frac{1}{2}{∥\mathbf{T}-{\mathbf{T}}^n∥}_F^2+t_1{\lambda }_1{∥\mathbf{T}∥}_{*}+t_1{\lambda }_2\Re (<\nabla \mathbf{T},{\mathbf{Y}}^n>)\hfill \\ & \frac{t_1}{2}\nabla f\left({\mathbf{T}}^n\right)+\frac{{\mathrm{t}}_{1}L}{2}{∥\mathbf{T}-{\mathbf{T}}^n∥}_F^2+t_1\Re (<\nabla f\left({\mathbf{T}}^n\right),\mathbf{T}-{\mathbf{T}}^n>)\hfill \end{array}$$

(24)

By neglecting constant terms, Equation (24) can be written as:

$$\begin{array}{cc}\hfill {\mathbf{T}}^{n+1}=& argmin\frac{1}{2}{∥\mathbf{T}-({\mathbf{T}}^n-\frac{t_1}{1+t_{1}L}\nabla f\left({\mathbf{T}}^n\right)∥}_F^2+\frac{t_1{\lambda }_2}{1+t_{1}L}\Re (<\nabla \mathbf{T},{\mathbf{Y}}^n>)\hfill \\ & +\frac{t_1{\lambda }_1}{1+t_{1}L}{∥\mathbf{T}∥}_{*}\hfill \end{array}$$

(25)

With the existence of $\nabla \mathbf{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 ${\mathcal{L}}^T$ that maps matrix pairs $\left(\mathbf{p},\mathbf{q}\right)$ to an m-by-n matrix as:

$${({\mathcal{L}}^T(\mathbf{p},\mathbf{q}))}_{\mathbf{i},\mathbf{j}}=p_{i,j}-p_{i-1,j}+q_{i,j}-q_{i,j-1}\mathrm{for}1\le i\le m,1\le j\le n$$

(26)

where $\mathbf{p}\in {\mathbb{C}}^{m-1\times n}$, $\mathbf{q}\in {\mathbb{C}}^{m\times n-1}$ and $p_{0,j}=p_{m,j}=q_{i,0}=q_{i,n}$.

The adjoint of the operator ${\mathcal{L}}^T$, denoted by $\mathcal{L}$, is a linear map from an m-by-n matrix $\mathbf{T}$ to a matrix pair $\left(\mathbf{p},\mathbf{q}\right)$ and is defined by:

$$\mathcal{L}\left(\mathbf{T}\right)=(p,q)$$

(27)

where:

$$\left\{\begin{array}{c}{p}_{i,j}=T_{i,j}-T_{i+1,j}\mathrm{for}1\le i\le m-1,1\le j\le n\hfill \\ q_{i,j}=T_{i,j}-T_{i,j+1}\mathrm{for}1\le i\le m,1\le j\le n-1\hfill \end{array}\right.$$

(28)

Based on the definition of the adjoint of operator $\mathcal{L}$ and the difference operator, the linear operator ${\mathcal{L}}^T$ satisfies:

$$<\nabla \mathbf{T},\mathbf{Y}>=<\mathbf{T},{\mathcal{L}}^T\left(\mathbf{Y}\right)>$$

(29)

Therefore, the dual variable $\mathbf{Y}$ is constructed by the matrix pair $\left(\mathbf{p},\mathbf{q}\right)$, where:

$${({\mathcal{L}}^T(\mathbf{Y}))}_{i,j}={({\mathcal{L}}^T(\mathbf{p},\mathbf{q}))}_{\mathbf{i},\mathbf{j}}=p_{i,j}-p_{i-1,j}+q_{i,j}-q_{i,j-1}$$

(30)

Based on Equations (29) and (30), we could simplify Equation (25) to the nuclear norm minimization problem:

$${\mathbf{T}}^{n+1}=argmin\frac{1}{2}{∥\mathbf{T}-{\mathbf{T}}_c^n∥}_F^2+\frac{t_1{\lambda }_1}{1+t_{1}L}{∥\mathbf{T}∥}_{*}$$

(31)

where:

$${\mathbf{T}}_c^n={\mathbf{T}}^n-\frac{t_1}{1+t_{1}L}\nabla f\left({\mathbf{T}}^n\right)-\frac{t_1}{1+t_{1}L}{\mathcal{L}}^T\left({\mathbf{Y}}^n\right)$$

(32)

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_{\tau }\left(\mathbf{M}\right)=arg\underset{\mathbf{X}}{min}\{\frac{1}{2}{∥\mathbf{X}-\mathbf{M}∥}_F^2+\tau {∥\mathbf{X}∥}_{*}\}$$

(33)

where $S_{\tau }(·)=diag(max{({\sigma }_i-\tau ,0)}_{1\le i\le 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 $\tau =\frac{t_1{\lambda }_1}{1+t_{1}L}\ge 0$, we have:

$${\mathbf{T}}^{n+1}=S_{\tau }\left({\mathbf{T}}_c^n\right)$$

(34)

Then, we consider the $\mathbf{Y}$ subproblem in Equation (21):

$${\mathbf{Y}}^{n+1}=argmin\frac{1}{2}{∥\mathbf{Y}-{\mathbf{Y}}^n∥}_F^2-t_2{\lambda }_2\Re (<\nabla (2{\mathbf{T}}^{n+1}-{\mathbf{T}}^n),\mathbf{Y}>)+I_{B\infty }\left(\mathbf{Y}\right)$$

(35)

After simplification, it becomes:

$${\mathbf{Y}}^{n+1}=argmin\frac{1}{2}{∥\mathbf{Y}-{\mathbf{Y}}_c^n∥}_F^2+I_{B\infty }\left(\mathbf{Y}\right)$$

(36)

where:

$${\mathbf{Y}}_c^n={\mathbf{Y}}^n+t_2{\lambda }_2\nabla (2{\mathbf{T}}^{n+1}-{\mathbf{T}}^n)$$

(37)

The solution of Equation (36) can be obtained by the Euclidean projection of ${\mathbf{Y}}_c^n$ onto a $l_{\infty }$ unit norm ball, which can be evaluated by:

$${\mathbf{Y}}^{n+1}=\mathrm{sgn}\left({\mathbf{Y}}_c^n\right)·min({\mathbf{Y}}_c^n,1)$$

(38)

where $\mathrm{sgn}(·)$ is the sign operator.

Now, the $\mathbf{T}$ and $\mathbf{Y}$ subproblems have been solved. The proposed InSAR-TVMC approach is summarized in

Table 1. InSAR-total variation regularized matrix completion (TVMC) imaging approach.

Input: D1, D2, V, ${\mathit{\lambda }}_{\mathbf{1}}$, ${\mathit{\lambda }}_{\mathbf{2}}$, N

Output: brightness temperature $\mathbf{T}$

Algorithm:

Step 1: initialize ${\mathbf{Y}}^{\mathbf{0}}$, ${\mathbf{T}}^{\mathbf{0}}$, and undersample $\mathbf{V}$ to obtain $\stackrel{\wedge }{\mathbf{V}}$;

Step 2: update ${\mathbf{T}}^n$ by using Equations (31) and (34);

Step 3: update ${\mathbf{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 ${\mathbf{T}}^N$

.

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. Simulation InSAR system parameters.

Simulation Parameters	Value
center frequency	50.3 GHz
bandwidth	200 MHz
image pixel size	128 × 128
image gray	(0, 255)
brightness temperature	2.73~350 K
antenna array (“T” shaped)	64 + 63
visibility function samples size	128 × 128
G matrix size	16,384 × 1
D matrix size	128 × 128

. 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 $\lambda =0.085$) and the InSAR-TVMC approach (the regularization parameters ${\lambda }_1=0.9$ and ${\lambda }_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 ${\sigma }^2$ representing all the noise and error in real interferometric measurements. Using a corrupted visibility function samples with two different levels of noise intensity (${\sigma }^2=0.05$ and ${\sigma }^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 (${\sigma }^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 ${\sigma }^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 (${\sigma }^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. Computation time comparison of CS and InSAR-TVMC.

Undersampling Rate	40%	50%	60%	70%	80%	90%
$t_{CS}\left(s\right)$	65.87	81.86	119.34	141.58	166.53	180.45
$t_{TVMC}\left(s\right)$	33.91	46.76	65.31	80.33	94.22	103.74

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 (${\sigma }^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.