GPR Antipersonnel Mine Detection Based on Tensor Robust Principal Analysis

Abstract

The ground Penetrating Radar (GPR) is a promising remote sensing modality for Antipersonnel Mine (APM) detection. However, detection of the buried APMs are impaired by strong clutter, especially the reflection caused by rough ground surfaces. In this paper, we propose a novel clutter suppression method taking advantage of the low-rank and sparse structure in multidimensional data, based on which an efficient target detection can be accomplished. We firstly created a multidimensional image tensor using sub-band GPR images that are computed from the band-pass filtered GPR signals, such that differences of the target response between sub-bands can be captured. Then, exploiting the low-rank and sparse property of the image tensor, we use the recently proposed Tensor Robust Principal Analysis to remove clutter by decomposing the image tensor into three components: a low-rank component containing clutter, a sparse component capturing target response, and noise. Finally, target detection is accomplished by applying thresholds to the extracted target image. Numerical simulations and experiments with different GPR systems are conducted. The results show that the proposed method effectively improves signal-to-clutter ratio by more than 20 dB and yields satisfactory results with high probability of detection and low false alarm rates.

1. Introduction

Being capable of sensing discontinuity of electromagnetic properties in the ground, Ground Penetrating Radar (GPR) has been used as nondestructive remote sensing modality in many applications, including tree root detection [1], frozen layer monitoring [2], structural assessment [3], infrastructure inspection [4], and landmine detection [5,6,7,8,9]. However, GPR detection of Antipersonnel Mines (APMs) is always challenging because of clutter contamination, such as antenna coupling, ground reflection, and friendly objects (rocks, voids, etc.) [10,11]. As a major source of clutter, the ground reflection is of large magnitude, and usually varies with position due to surface roughness [11]. Since many APMs are of low-metal content and shallow burial, the responses of APMs are often obscured by ground reflections, rendering target visualization and detection very difficult.

Many efforts have been devoted to developing clutter mitigation methods, and these methods can be briefly categorized into three types: model-based methods, filtering, and subspace projection [11,12,13].

In the model-base methods, clutter models are defined, and objects whose responses significantly deviate from the clutter models are detected as potential targets. These methods vary depending on the clutter model chosen to capture the characteristics of the clutter. In [14], a Wide-Sense Stationary (WSS) random process is adopted to describe the clutter whereas in [15], the background clutter is estimated using the least squares technique. A linear prediction model with adaptive coefficients is proposed to characterize the clutter in [5]. However, the model-based methods need “clean” clutter (received signal that contains clutter only) to estimate model parameters which tends to be contaminated by other factors. In addition, the performances of these methods depend strongly on the clutter models and degrade when clutter models deviate from the real situations.

The filtering methods adopt filers in a certain domain to filter out clutter while preserving the target response. An intuitive technique termed time-gating is to mute the received signal within a time window [16]. The problem is how to set a proper time window that covers clutter without compromising target reflections. In [13], the entropy-based criterion for time window selection is proposed under the assumption that clutter has similar behavior for each trace. Another popular method is mean subtraction, which can be viewed as a high-pass filter in the frequency domain, but the calculation of average background clutter trace can be affected by the target responses [16]. A physically-based filtering method is presented in [11] which uses frequency-dependent linear transfer functions to eliminate the strong ground clutter as well as the antenna coupling by subtraction of a simulated Green function. To obtain the linear transfer functions, calibration under controlled environments is required.

The subspace projection method assumes that the ground clutter and the target responses are in different subspaces. Therefore, clutter and target responses can be separated by projecting the recorded GPR data into clutter and target subspaces, respectively. The Principal Component Analysis (PCA) is a representative subspace projection method [17,18]. Usually, the subspace spanned by the first dominant eigencomponent is considered as a clutter subspace to be removed. In practice, the uneven ground surface results in nonstationary clutter whose rank is greater than one and the PCA method cannot remove the clutter completely. Moreover, the PCA method is so sensitive to strong noise that even a single grossly corrupted observation can lead to error results [19]. As an improvement of the traditional PCA, the Robust Principle Component Analysis (RPCA) is proposed for recovering a low-rank subspace from noisy measurements [19]. For GPR data traces, the ground clutter as well as antenna coupling are considered as low-rank, while the target responses are assumed to be contained in the sparse component [20]. Hence, the clutter suppression problem is cast as a low-rank and sparse representation problem and can be solved by the RPCA-based method [12,21]. The recently proposed Tensor Robust Principal Component Analysis (TRPCA) [22,23] provides an alternative perspective of solving the low-rank and sparse decomposition problem. Extended from the RPCA, the TRPCA operates on tensors to seek a low-rank tensor embedded in sparsely corrupted multidimensional data. Avoiding the information loss caused by restructuring high order tensor into matrix, TRPCA makes full use of the abundant information contained in multidimensional data to guarantee flexibility as well as accuracy in video processing [24], seismic data denoising [25], and infrared target detection [26].

Based on TRPCA, this paper presents a GPR APM detection method utilizing the low-rank and sparse properties of multidimensional data for the separation of target response and background clutter. To solve the clutter suppression and target detection problem under tensor algebra framework, we first propose a tensor formation scheme using a sub-band technique. Inspired by the observation that targets give distinct responses at different frequency bands [5,27,28], we exploit sub-band images to characterize the frequency dependent features of APMs. Specifically, the whole frequency band of the recorded GPR signal is divided into several sub-bands to compute individual GPR image for each sub-band and then an image tensor can be formed with these sub-band images inserted as frontal slices (definition in

Table 1. Metamathematical notations.

Notations	Definition
$\mathcal{A},\mathit{A},\mathit{a},a$	tensor, matrix, vector, scalar
$\mathbb{R},\mathbb{C}$	real numbers and complex numbers
${\mathcal{A}}_{i,j,k}$ or $a_{i,j,k}$	the $(i,j,k)$th entry of $\mathcal{A}$
$\mathcal{A}(i,j,:)$	the tube of $\mathcal{A}$
$\mathcal{A}(i,:,:),\mathcal{A}(:,i,:)$	the ith horizontal or lateral slice of $\mathcal{A}$
$\mathcal{A}(:,:,i)$ or ${\mathcal{A}}^{\left(i\right)}$	the ith frontal slice of $\mathcal{A}$
$\overline{\mathcal{A}}=\mathrm{fft}\left(\mathcal{A},\left[\right],3\right)$	Fast Fourier Transformation (FFT) on $\mathcal{A}$ along the third dimension
$\mathcal{A}=\mathrm{ifft}\left(\overline{\mathcal{A}},\left[\right],3\right)$	Inverse FFT on $\mathcal{A}$ along the third dimension
${\parallel \mathcal{A}\parallel }_1$	${\ell }_1$ norm of $\mathcal{A}$, which is computed as $\sum {}_{i,j,k}\left|a_{i,j,k}\right|$
${\parallel \mathcal{A}\parallel }_{\infty }$	infinity norm of $\mathcal{A}$, which is defined as $\max |a_{i,j,k}|$
${\parallel \mathcal{A}\parallel }_F$	Frobenius norm of $\mathcal{A}$, which is computed as $\sqrt{\sum {}_{i,j,k}{\left|a_{i,j,k}\right|}^2}$
${\parallel \mathit{A}\parallel }_{\ast }$	nuclear norm of $\mathit{A}$, which is computed as sum of singular values
${\parallel \mathit{A}\parallel }_2$	spectral norm of $\mathit{A}$, which is defined as the largest singular value

). By this means, the spatial-spectral information is incorporated in the images tensor, which can be simultaneously handled by a tensor based method to facilitate the detection. Then, TRPCA is used to decompose the obtained image tensor into low-rank component and sparse component to separate the target response (sparse) from the clutter (low-rank). Since most unwanted factors are eliminated, the final detection result can be acquired by applying thresholds on the target image. The proposed method is an extension of the RPCA-based clutter suppression method [21] and the contributions are threefold: (1) a tensor formation scheme is proposed to enable the spatial-spectral processing with TRPCA, and as far as we know, the application of TRPCA in GPR APM detection has not yet been reported; (2) TRPCA is adapted to the APM detection scenario and the problem is solved by a tractable optimization—the Alternating Direction Method of Multiplier (ADMM); (3) our proposed method achieves high Probability of Detection (PD) while effectively reducing the False Alarm Rate (FAR). Numerical simulations and laboratory experiments with impulse GPR and Stepped Frequency Continuous Wave (SFCW) GPR were conducted to validate our proposed method. Both of low-metal APMs and nonmetallic APMs were used as targets and were successfully detected in the experiments. It is noteworthy that the goal of detection is to flag areas that may contain anomalies [29]. Therefore, the subsurface (and surface) anomalies, including APMs and friendly objects, are considered as suspected targets, and the clutter refers to antenna coupling and ground reflection in this paper. The rejection of friendly objects is to be done by a classifier, which is topic of our future work.

The paper is organized as follows. Section 2 introduces the notations and briefly reviews the TRPCA method. The mathematical formulation and the proposed TRPCA-based clutter suppression and target detection method are illustrated in Section 3. Section 4 presents the experiment results as well as discussions of the results to demonstrate the effectiveness of our proposed method. Conclusions is drawn in Section 5.

2. Notations and Preliminaries

2.1. Notations

The notations used in this paper are listed in the Table 1. Below are some necessary definitions.

For a tensor $\mathcal{A}\in {\mathbb{R}}^{n_1\times n_2\times n_3}$, we define a block diagonal matrix $\overline{\mathit{A}}\in {\mathbb{R}}^{n_1{n}_3\times n_2{n}_3}$ whose ith block on diagonal is the ith frontal slice ${\overline{\mathit{A}}}^{\left(i\right)}$ of $\overline{\mathcal{A}}$

$$\overline{\mathit{A}}=\left[\begin{array}{c}{\overline{\mathit{A}}}^{\left(1\right)}\\ {\overline{\mathit{A}}}^{\left(2\right)}\\ \ddots \\ {\overline{\mathit{A}}}^{\left(n_3\right)}\end{array}\right].$$

(1)

A block circular matrix l$\mathbf{bcirc}\left(\mathcal{A}\right)$ of size $n_1{n}_3\times n_2{n}_3$ is denoted as

$$\mathbf{bcirc}\left(\mathcal{A}\right)=\left[\begin{array}{cccc}{\mathit{A}}^{\left(1\right)}& {\mathit{A}}^{\left(2\right)}& \cdots & {\mathit{A}}^{\left(n_3\right)}\\ {\mathit{A}}^{\left(2\right)}& {\mathit{A}}^{\left(1\right)}& \cdots & {\mathit{A}}^{\left(n_3\right)}\\ ⋮& ⋮& \ddots & ⋮\\ {\mathit{A}}^{\left(n_3\right)}& {\mathit{A}}^{(n_3-1)}& \cdots & {\mathit{A}}^{\left(1\right)}\end{array}\right].$$

(2)

Two operators, unfold and fold are respectively defined as

$$\mathbf{unfold}\left(\mathcal{A}\right)=\left[\begin{array}{c}{\mathit{A}}^{\left(1\right)}\\ {\mathit{A}}^{\left(2\right)}\\ ⋮\\ {\mathit{A}}^{\left(n_3\right)}\end{array}\right]$$

(3)

$$\mathbf{fold}\left(\mathbf{unfold}\right(\mathcal{A}\left)\right)=\mathcal{A}.$$

(4)

Then, we have a new tensor-tensor production termed as t-product.

Definition 1

(t-product [22]). Let $\mathcal{A}\in {\mathbb{R}}^{n_1\times n_2\times n_3}$ and $\mathcal{B}\in {\mathbb{R}}^{n_2\times n_4\times n_3}$. Then the t-product $\mathcal{C}=\mathcal{A}\ast \mathcal{B}$ is a tensor of size $n_1\times n_4\times n_3$ defined as

$$\mathcal{C}=\mathcal{A}\ast \mathcal{B}=\mathbf{fold}\left(\mathbf{bcirc}\left(\mathcal{A}\right)·\mathbf{unfold}\left(\mathcal{B}\right)\right).$$

(5)

Using the t-product, a tensor factorization technique named tensor-SVD (t-SVD) can be defined.

Definition 2

(t-SVD [30]). For $\mathcal{A}\in {\mathbb{R}}^{n_1\times n_2\times n_3}$, the t-SVD of $\mathcal{A}$ is given by

$$\mathcal{A}=\mathcal{U}\ast \mathcal{S}\ast {\mathcal{V}}^T$$

(6)

where $\mathcal{U}\in {\mathbb{R}}^{n_1\times n_1\times n_3}$ and $\mathcal{V}\in {\mathbb{R}}^{n_2\times n_2\times n_3}$ are orthogonal tensors, and $\mathcal{S}$ is a f-diagonal tensor, T denotes the transpose of a tensor. The detail descriptions of orthogonal tensor, f-diagonal tensor, and tensor transpose are given in [22,30].

Based on t-SVD, the tensor tubal rank as well as the Tensor Nuclear Norm (TNN) is defined.

Definition 3

(tensor tubal rank [22]). The tensor tubal rank of $\mathcal{A}\in {\mathbb{R}}^{n_1\times n_2\times n_3}$, denoted as ${\mathrm{rank}}_t\left(\mathcal{A}\right)$ is defined as the number of nonzero singular tubes of $\mathcal{S}$, where $\mathcal{S}$ is from the t-SVD of $\mathcal{A}=\mathcal{U}\ast \mathcal{S}\ast {\mathcal{V}}^{\mathrm{T}}$

$${rank}_t\left(\mathcal{A}\right)=\#\{i:\mathcal{S}(i,i,:)\ne 0\}.$$

(7)

Definition 4

(TNN [23]). The tensor nuclear norm of tensor $\mathcal{A}\in {\mathbb{R}}^{n_1\times n_2\times n_3}$, denoted as ${\parallel \mathcal{A}\parallel }_{\ast }$, is defined as the average of the nuclear norm of all the frontal slices of $\overline{\mathcal{A}}$

$${\parallel \mathcal{A}\parallel }_{\ast }=\frac{1}{n_3}\sum _{i=1}^{n_3}{\parallel {\overline{\mathit{A}}}^{\left(i\right)}\parallel }_{\ast }.$$

(8)

We can further define the tensor spectral norm as

Definition 5

(tensor spectral norm [23]). The tensor spectral norm of tensor $\mathcal{A}\in {\mathbb{R}}^{n_1\times n_2\times n_3}$, denoted as ${\parallel \mathcal{A}\parallel }_2$, is defined as the spectral norm of the block diagonal matrix $\overline{\mathit{A}}$

$${\parallel \mathcal{A}\parallel }_2={\parallel \overline{\mathit{A}}\parallel }_2.$$

(9)

2.2. Principle of the TRPCA

Extended from the matrix RPCA, TRPCA considers a similar problem that is to recover a low-rank component corrupted by sparse errors. Different from the RPCA, which handles with matrices, TRPCA manipulates tensor data to take the advantage of abundant information contained in the multidimensional data structure [23]. More specifically, TRPCA aims to separate the low tubal rank tensor $\mathcal{L}$ and the sparse component $\mathcal{S}$ from the observation $\mathcal{X}$

$$\mathcal{X}=\mathcal{L}+\mathcal{S}.$$

(10)

Based on t-SVD, the above mentioned problem can be addressed by solving the following convex optimization, where the tensor tubal rank is replaced by the TNN [23]

$$\underset{\mathcal{L},\mathcal{S}}{\min }{\parallel \mathcal{L}\parallel }_{\ast }+\lambda {\parallel \mathcal{S}\parallel }_1,\mathrm{s}.\mathrm{t}.\mathcal{X}=\mathcal{L}+\mathcal{S}.$$

(11)

However, the solution adopted in [23] is aiming at the exactly recovery of tensor. TRPCA need to be adapted to the GPR APM detection application and take the ubiquitous noise in practical measurements into consideration.

Note that when the tensor reduces to matrix, i.e., the third dimension is one, the t-product reduces to standard matrix product and the tensor nuclear norm reduces to matrix nuclear norm. Therefore, RPCA is a special case of TRPCA [23]. The consistency between TRPCA and RPCA helps a lot in our GPR APM detection application, which is to be shown in the later sections.

3. Proposed Method

The proposed ground clutter suppression and target detection method includes three procedures: tensor construction, clutter suppression, and target detection. To begin with, the mathematical formulation of GPR APM detection problem under the tensor framework will be introduced.

3.1. Problem Formulation and Tensor Construction

The scenario of GPR detection is depicted in Figure 1. Emitting and receiving signals at each location, the down-looking GPR scans along a line over the ground surface. The recorded data traces are collected as columns to form a two dimensional matrix (B-scan), ldenoted by $\mathit{X}$

$$\mathit{X}={\mathit{X}}_1+{\mathit{X}}_2+{\mathit{X}}_3$$

(12)

where ${\mathit{X}}_1$, ${\mathit{X}}_2$, and ${\mathit{X}}_3$ are the clutter, the target response, and noise, respectively. As mentioned above, the clutter that consist of ground reflection and antenna corss-talk is considered as low-rank component in the B-scan, whereas the target response are contained in sparse component [12,21,31]. Therefore, the low-rank and sparse property of GPR image data is utilized to discriminate the clutter ${\mathit{X}}_1$ and the target response ${\mathit{X}}_2$ by RPCA decomposition.

To further improve the performance of GPR APM detection, we consider incorporating more information contained in the spectrum by constructing an image tensor with multiple sub-band images. As reported in [5,27,28], the sub-band processing takes the distinct features of APMs at different frequency bands into account, thereby increasing detection accuracy. To capture the frequency-dependent characteristics of target responses, sub-band images are formed. Specifically, as shown in Figure 2, the recorded GPR data $\mathit{X}$ is first filtered by two band-pass filters, which divide the whole frequency band into the Low-Frequency (LF) band and the High-Frequency (HF) band while permitting overlap. Two sub-band GPR images, i.e., ${\mathit{X}}_{LF}$ and ${\mathit{X}}_{HF}$, are formed by focusing the LF and the HF GPR data, respectively. Then, we stack the sub-band images as frontal slices, that is, $\mathcal{X}(:,:,1)={\mathit{X}}_{LF}$ and $\mathcal{X}(:,:,2)={\mathit{X}}_{HF}$, to create an image tensor resembling a hyperspectral image cube [32], where the spatial characteristics lies in the first two dimensions and the spectral information is contained in the third dimension. According to the low-rank and sparse prior of GPR data [12,21], the obtained GPR image tensor $\mathcal{X}$ can be represented as

$$\mathcal{X}={\mathcal{X}}_1+{\mathcal{X}}_2+{\mathcal{X}}_3$$

(13)

where the ${\mathcal{X}}_1$, ${\mathcal{X}}_2$, and ${\mathcal{X}}_3$ represent the low-rank clutter tensor, the sparse target image tensor, and noise, respectively. Therefore, the clutter suppression and target detection is cast as a TRPCA problem to separate the sparse target from the low-rank background clutter.

3.2. Clutter Suppression via TRPCA

In this section, clutter is removed by solving the aforementioned TRPCA problem Equation (13) via a convex optimization program. Mathematically,

$$\begin{array}{cc}\underset{{\mathcal{X}}_1,{\mathcal{X}}_2,{\mathcal{X}}_3}{\min }\hfill & \parallel {\mathcal{X}}_1{\parallel }_1+{\lambda }_2\parallel {\mathcal{X}}_2{\parallel }_{\ast }+{\lambda }_3{\parallel {\mathcal{X}}_3\parallel }_F^2\hfill \\ & \mathrm{s}.\mathrm{t}.\mathcal{X}={\mathcal{X}}_1+{\mathcal{X}}_2+{\mathcal{X}}_3.\hfill \end{array}$$

(14)

The constrained problem Equation (14) can be addressed by the augmented Lagrangian function [33]

$$\begin{array}{cc}& \mathcal{L}({\mathcal{X}}_1,{\mathcal{X}}_2,{\mathcal{X}}_3,\mathcal{U};\rho ,{\lambda }_2,{\lambda }_3)\hfill \\ =\hfill & \sum _{i=1}^N{\psi }_i\left({\mathcal{X}}_i\right)+\frac{\rho }{2}{\parallel \mathcal{X}-\sum _{i=1}^N{\mathcal{X}}_i+\mathcal{U}\parallel }_F^2\hfill \end{array}$$

(15)

where $i=1,\cdots ,N$, $N=3$, $\rho >0$, and functions ${\psi }_i,i=1,\cdots ,3$ are defined as

$$\begin{array}{cc}{\psi }_1\left(\mathcal{X}\right)\hfill & ={\parallel \mathcal{X}\parallel }_{\ast }\hfill \end{array}$$

(16a)

$$\begin{array}{cc}{\psi }_2\left(\mathcal{X}\right)\hfill & ={\lambda }_2{\parallel \mathcal{X}\parallel }_1\hfill \end{array}$$

(16b)

$$\begin{array}{cc}{\psi }_3\left(\mathcal{X}\right)\hfill & ={\lambda }_3{\parallel \mathcal{X}\parallel }_F^2\hfill \end{array}$$

(16c)

where ${\lambda }_2>0$ and ${\lambda }_3>0$ are tuning parameters.

We use the ADMM [33] to solve the problem Equation (15) by alternately solving the $({\mathcal{X}}_1,{\mathcal{X}}_2,{\mathcal{X}}_3,\mathcal{U})$

$${\mathcal{X}}_i^{k+1}={\mathrm{argmin}}_{{\mathcal{X}}_i}\left({\psi }_i\left({\mathcal{X}}_i\right)+\frac{\rho }{2}{\parallel {\mathcal{X}}_i-{\mathcal{D}}_i^k\parallel }_F^2\right)$$

(17)

$${\mathcal{U}}^{k+1}={\mathcal{U}}^k+{\tilde{\mathcal{X}}}^{k+1}-\frac{1}{N}\mathcal{X}$$

(18)

where $i=1,\cdots ,N$, $N=3$, the average of tensors ${\tilde{\mathcal{X}}}^k=(1/N){\sum }_{i=1}^N{\mathcal{X}}_i^k$, and ${\mathcal{D}}_i^k$ is defined as

$${\mathcal{D}}_i^k={\mathcal{X}}_i^k-{\tilde{\mathcal{X}}}^k+\frac{1}{N}\mathcal{X}-{\mathcal{U}}^k.$$

(19)

To update the ${\mathcal{X}}_1^{k+1}$, we use the tensor Singular Value Thresholding (t-SVT) [34]

$${\mathcal{X}}_1^{k+1}={\mathcal{F}}_{\frac{1}{\rho }}\left({\mathcal{D}}_1^k\right)$$

(20)

where the t-SVT operator ${\mathcal{F}}_{1/\rho }$ is defined as follows: Let $\mathcal{X}=\mathcal{U}\ast \mathcal{S}\ast {\mathcal{V}}^{\mathrm{T}}$ be the t-SVD of $\mathcal{X}\in {\mathbb{R}}^{n_1\times n_2\times n_3}$. For each $1/\rho >0$, the t-SVT on $\mathcal{X}$ is

$${\mathcal{F}}_{\frac{1}{\rho }}\left(\mathcal{X}\right)=\mathcal{U}\ast {\mathcal{S}}_{\frac{1}{\rho }}\ast {\mathcal{V}}^{\mathrm{T}}$$

(21)

where

$${\mathcal{S}}_{\frac{1}{\rho }}=\mathrm{ifft}\left({(\overline{\mathcal{S}}-\frac{1}{\rho })}_{+},\left[\right],3\right).$$

(22)

The subscript ${\left(x\right)}_{+}$ denotes the positive part of x.

The ${\mathcal{X}}_2^{k+1}$ is updated using the elementwise soft-shrinkage operator [33]

$${\mathcal{X}}_2^{k+1}=\mathbf{shrink}\left({\mathcal{D}}_2^k,\frac{{\lambda }_2}{\rho }\right)$$

(23)

where the $\mathbf{shrink}(x,{\lambda }_2/\rho )$ is

$$\mathbf{shrink}\left(x,\frac{{\lambda }_2}{\rho }\right)=\left\{\begin{array}{ccc}x\hfill & -\frac{{\lambda }_2}{\rho }\hfill & x>\frac{{\lambda }_2}{\rho }\hfill \\ 0\hfill & \hfill \left|x\right|\le \frac{{\lambda }_2}{\rho }\\ x\hfill & +\frac{{\lambda }_2}{\rho }\hfill & x<-\frac{{\lambda }_2}{\rho }.\hfill \end{array}\right.$$

(24)

For ${\mathcal{X}}_3^{k+1}$, we have the closed form solution as [33]

$${\mathcal{X}}_3^{k+1}=\frac{1}{1+{\lambda }_3/\rho }{\mathcal{D}}_3^k.$$

(25)

Finally, the solving procedures are described in Algorithm 1.

Algorithm 1 The ADMM for TRPCA

Input:    Image tensor $\mathcal{X}$ Output:    Target image tensor ${\mathcal{X}}_2$ Initialize:    ${\mathcal{X}}_1^0={\mathcal{X}}_2^0={\mathcal{X}}_3^0=0,{\mathcal{U}}^0=0,\rho =1,ϵ=1e^{-6},$ $m=200,{\lambda }_2={0.15\ast \parallel \mathcal{X}\parallel }_2,{\lambda }_3=0.15\ast {\parallel \mathcal{X}\parallel }_{\infty }$ 1: while not converged do 2:    ${\mathcal{X}}_1^{k+1}\leftarrow$${\mathcal{F}}_{1/\rho }\left({\mathcal{D}}_1^k\right)$ 3:    ${\mathcal{X}}_2^{k+1}\leftarrow$$\mathrm{shrink}\left({\mathcal{D}}_2^k,{\lambda }_2/\rho \right)$ 4:    ${\mathcal{X}}_3^{k+1}\leftarrow {\mathcal{D}}_3^k/(1+{\lambda }_3/\rho )$ 5:    ${\mathcal{U}}^{k+1}\leftarrow {\mathcal{U}}^k$$+{\tilde{\mathcal{X}}}^{k+1}-\mathcal{X}/N$ 6:     $c_1\leftarrow {\parallel {\mathcal{X}}_2^{k+1}-{\mathcal{X}}_2^k\parallel }_{\infty }$     $c_2\leftarrow {\parallel {\mathcal{X}}_3^{k+1}-{\mathcal{X}}_3^k\parallel }_{\infty }$     $c_3\leftarrow {\parallel {\mathcal{U}}^{k+1}-{\mathcal{U}}^k\parallel }_{\infty }$ 7:    if $\max (c_1,c_2,c_3)<ϵ\mathrm{or}k>m$ then 8:      Break 9:    else 10:      $k\leftarrow k+1$ 11:    end if 12: end while

3.3. Target Detection

After the TRPCA decomposition, a clear target image tensor with most of the clutter removed is obtained, yet, there is residual clutter that may raise false alarms. To facilitate the detection, we sum up the frontal slices of the target image tensor to obtain a target image, that is, ${\mathit{X}}_2={\sum }_{i=1}^2{\mathcal{X}}_2(:,:,i)$.

With a clear target image, detection can be accomplished as follows. First, an energy detector is applied to the target image. Top ${\alpha }_1$ percent of pixels are retained. Then, remove isolated detections whose areas are smaller than ${\alpha }_2$ (pixels). At last, cluster the detections and recalculate the center of mass. Those closely spaced detections with distances between each other that are smaller than ${\alpha }_3$ (centimeters) are clustered.

As shown in Figure 2, the procedures of our method can be summarized as follows:

• Tensor construction: Equally divide the bandwidth into two sub-bands and form a sub-band GPR image on each sub-band data using the Layered Range Migration (LRM) method [35]. Create the three-dimensional image tensor by inserting the sub-band images as frontal slices.
• TRPCA: Decompose the formed image tensor to separate the low-rank component ${\mathcal{X}}_1$ and the sparse component ${\mathcal{X}}_2$. The target response is contained in the sparse component ${\mathcal{X}}_2$.
• Target detection: Calculate the target image ${\mathit{X}}_2$ by summing up the frontal slices of ${\mathcal{X}}_2$. Remove clutter residue with thresholds and finally locate the target. The center coordinates of the remaining nonzero regions in the target image are output as the detection results.

4. Results and Discussion

To coroborate the proposed clutter suppression and target detection method, both numerical simulations and experiments using different GPR systems were carried out. The simulation data were synthesized using the gprMax [36]. In the experiments, two datasets were collected and tested: one was acquired with an impulse GPR developed by our research group; the other was collected by a SFCW GPR and published by the Georgia Technology Institute [37,38].

4.1. Metrics and Baselines

The Signal to Clutter Ratio (SCR) is adopted to assess the effectiveness of clutter suppression processing. The SCR is calculated as

$$\mathrm{SCR}=\frac{N_c}{N_t}\sum _{p\in R_t}{\left|I\left(p\right)\right|}^2/\sum _{p\in R_c}{\left|I\left(p\right)\right|}^2$$

(26)

where $I\left(p\right)$ is the p-th pixel in image. $N_t$ and $N_c$ denotes the number of pixels in the target region $R_t$ and the number of pixels in the clutter region $R_c$, respectively. The target region is indicated by a box (or boxes in a scenario with multiple targets) enclosing the location of target. The size of the box is manually set to cover the target image. The entire image excluding the target region is defined as the clutter region.

The metrics of evaluating detection performance are PD $P_d$ and FAR $F_a$, which are defined as

$$P_d=\frac{\mathrm{number} \mathrm{of} \mathrm{correct} \mathrm{detections}}{\mathrm{number} \mathrm{of} \mathrm{actual} \mathrm{targets}}$$

(27)

$$F_a=\frac{\mathrm{number} \mathrm{of} \mathrm{false} \mathrm{alarms}}{\mathrm{total} \mathrm{length} \mathrm{of} \mathrm{the} \mathrm{surveyed} \mathrm{scenario}}$$

(28)

Note that $F_a$ is calculated on a meter basis. For example, a 1% FAR means that, on average, there is one false alarm per 100 m.

For comparison purposes, two representative clutter suppression and target detection methods, the PCA-CFAR method and the RPCA-based method, are adopted as baseline methods. The PCA-CFAR method utilizes the PCA [17] to remove clutter from a GPR image by discarding the largest principal component, and then detects target by Cell Averaging Constant FAR (CA-CFAR) technique [39]. The RPCA-based detection method [40] first extracts a target image from clutter using the Go-Decomposition (GoDec) [41] and then applies thresholds to reveal the target location. Same thresholds ${\alpha }_1$, ${\alpha }_2$, and ${\alpha }_3$ are used here.

4.2. Numerical Simulations

Shown as Figure 1, point source and receiver spaced 2 cm apart were 15 cm distant from the ground and scanned with 1 cm step length. A Ricker waveform with 2 GHz center frequency and 1.6 GHz bandwidth were adopted. The time window is set to 3.11 ns and the sampling interval is 2.59 ps. The soil is considered homogeneous, and the relative permittivity and conductivity of soil are $ϵ=5$ and $\sigma =0.01$ S/m, respectively. As shown in Figure 1, the APM was modeled as a hollow cylinder containing another cylinder inside. The outer cylinder was of radius $r_1=3.5$ cm and height $h_1=4.5$ cm, while the inner cylinder was of $r_2=2.5$ cm radius and $h_2=3$ cm height. The relative permittivity of the outer and the inner cylinders were ${ϵ}_1=2.9$ and ${ϵ}_2=3.6$, respectively.

One thousand B-scans were synthesized with the parameters listed above. Each B-scan is of the same size with 60 columns and 1198 rows, that is, it covers a length of 60 cm. The rough ground surfaces in these simulations were randomly generated using the fractal correlated noise model [42]. The fractal dimension was set to 1.5 and variation in height was 6 cm. In addition, Gaussian noise was added to the data and thus the Signal to Noise Ratio (SNR) was reduced to 0 dB.

Figure 3a displays one of the simulated noise-free B-scans where an APM is flush buried. Therefore, a part of ground reflections coincides with the target responses in time domain and results in waveform distortion. After noise is added, the target responses in Figure 3b are overwhelmed and hardly any hyperbola signature of the target are extracted. Using the LRM to focus the GPR data yields a full band imaging result (Figure 3c), where the target size and location are indicated by a dash-dotted box. Due to obscuration of the dominant ground clutter, the target image cannot be observed. Similar to [27], we equally divide the frequency bandwidth into two sub-bands while allowing 50% overlap: the LF band of 0.1–3.8 GHz and the HF band of 1.9–5.6 GHz. The optimal selection of the overlapping range will be discussed at the end of this section. As shown Figure 3d,e, two sub-band images, i.e., the LF image and the HF image are computed using these sub-band GPR data. Note that the uneven surface has similar shapes in these images, whereas the obscuration caused by the ground clutter renders target detection difficult. Then, we create an image tensor with these sub-band images and decompose the image tensor with TRPCA to obtain the target image tensor. In order to display the decomposition results, frontal slices of the obtained target image tensor are summed up to form a two-dimensional image, as shown in Figure 3f. The target image in Figure 3f is clear enough to recognize the APM. As listed in

Table 2. SCR of different clutter suppression methods.

	Imaing	The Proposed Method	PCA-CFAR	RPCA-Based Method
Simulated data 1	5.6	29.7	12.4	20.3
Simulated data 2	6.7	31.5	11.4	21.0
Impulse GPR data 1	−4.0	22.3	0.1	–
Impulse GPR data 2	−0.5	30.0	12.6	19.8
SFCW GPR data	−10.6	26.8	0.4	17.3

, the SCR of Figure 3f is improved by 24 dB compared with the full band image shown in Figure 3c. Since the major part of the clutter is removed, we can use thresholds to remove clutter residues and get the final detection result. As shown in Figure 3g, the target is indicated by a point.

The reasons why TRPCA performs well in the environment with heavy clutter are due to two aspects: (1) the low-rank and sparse properties of the GPR images, and (2) the frequency dependency of target response. Since the model and theory of TRPCA are consistent with RPCA [23], the low-rank and sparse structure of GPR image can be utilized by the TRPCA for the separation of target image and clutter in each sub-band image (frontal slice), which is similar to the RPCA-based methods [12,21]. The proposed method uses imaging to enhance the sparsity and intensity of a target response, which generally leads to a better decomposition [19]. Moreover, our method takes it one step further than the RPCA-based method by making use of the spectral information of target response lies in the third dimension. Though the differences caused by frequency dependency of target response are hard to discern by naked eyes in the sub-band images (Figure 3d,e), they can be captured and revealed by the TRPCA, which is intended for extraction of the variations in frontal slices [22].

A different situation is shown in Figure 4a where an APM is buried deeper (7 cm beneath the ground) due to the raised surface. The ground clutter and target responses are separated in time, but two pits in the ground at $x=14$ cm and $x=42$ cm causes strong reflections with hyperbolic curve signatures. After imaging, the target images can be found in the centers of sub-band images (Figure 4d,e). However, the clutter, especially images of the pits are of large magnitude that can raise false alarms. To remove the ground clutter and avoid false alarms, our proposed method is applied to improve SCR of the target image in Figure 4f by 24.8 dB, compared to the full band image Figure 4c. The target location is marked by a point in Figure 4g after the threshold-based detection procedure.

Processing all of the one thousand simulated B-scans with our proposed method, the detection performance in the low SNR (0 dB) environment is of 99.5% PD and 10% FAR, as shown in

Table 3. Comparison of detection performance (PD/FAR).

	The Proposed Method	PCA-CFAR	RPCA-Based Method
Simulated dataset	99.5%/10.0%	87.0%/67.0%	93.6%/27.7%
Impulse GPR dataset	100%/12.3%	89.4%/54.5%	93.3%/10.8%
SFCW GPR dataset	100%/7%	84.6%/19.2%	100%/8.4%

.

The PCA-CFAR method and the recently proposed RPCA-based method were tested on the same synthesized dataset. The full band images (Figure 3c and Figure 4c) were used as the inputs for these two methods. Figure 3h shows the detection result of the PCA-CFAR method in the shallow-burying scenario, where the target and false alarms are indicated with point and triangles, respectively. Due to the uneven ground surface, the rank of clutter is bigger than one and thus cannot be fully removed because only the first principal component is removed. In addition, the traditional PCA is so sensitive to strong noise that the processed image is contaminated by remaining noise. Therefore, together with the target marked by point, there is residual clutter and noise which degrades SCR and causes false alarms. This is similar to the case with deeply buried targets; there is clutter and noise residue in Figure 4h due to the deficiencies of PCA-CFAR. The improvements of SCR of the PCA-CFAR method in these two scenarios are listed in Table 2, which are 6.8 dB and 4.7 dB, respectively. The detection performance on the whole dataset is low PD (87.0%) with a high FAR of 67.0%. It shows that the PCA-CFAR method cannot handle the clutter caused by an uneven ground surface very well and its performance degrades severely under low SNR.

Thanks to the low-rank and sparse signal model, the RPCA-based method outperforms the PCA-CFAR method in the simulation by yielding higher SCRs (20.3 dB in Figure 3i, 21 dB in Figure 4i) and better detection performance (93.6% PD and 27.7% FAR). However, the variation of surface height results in nonstationary ground clutter whose low-rank property is degraded and thus, the RPCA-based method performs inadequately. This is the reason why both bulge and cavity in the ground cause false alarms in detection results Figure 3i and Figure 4i. Thus, in the numerical simulation, the detection performance of RPCA-based method is inferior to that of the proposed TRPCA-based method.

The time consumptions of these three methods are listed in

Table 4. Time consumption of different methods.

	The Proposed Method	PCA-CFAR	RPCA-Based Method
Time consumption (s)	0.24	3.8	0.14

. The solving programs of these methods are running on a same computer with 2.4 GHz CPU and 8 GB memory. It takes 0.24 s for the proposed method to complete the detection on a single B-scan, whereas the PCA-CFAR method takes 3.7 s. By using a accelerated solution program [41], time consumption of the RPCA-based method is 0.14 s. Note that the aforementioned optimization program ADMM can be further improved by some techniques [32,43,44]. However, these discussions are out of the scope of this paper.

4.3. Laboratory Experiments

Laboratory experiments using different GPR systems were conducted to demonstrate the effectiveness of our proposed method on real data.

In the first experiment, an impulse GPR of 1.6 GHz center frequency and 1.2 GHz bandwidth was used to scan a sand pit. The horn antennas were 15 cm above the surface. To deal successfully with the strong ground reflections, a dynamic range of over 69 dB is expected [45]. The space sampling step length was 5 mm and the time sampling interval was 9.77 ps. As Figure 5 shows, the pit was filled with dry sand, whose relative permittivity was 4. The surface of sand was flat. Seventy-two B-scans containing 90 APMs were collected. Each B-scan covered a length of 0.9 m, where only one APM or at most two APMs were buried. Both of metal APMs and nonmetallic APMs were used as targets in the experiment. We take two B-scans out from the dataset as examples to show the feasibility of our method.

As shown in Figure 6a, one of the B-scans contains a plastic APM buried 5 cm under the ground. Indicated by a dash-dotted box, the APM is of 4 cm height and 3.5 cm radius. We can hardly notice the subtle target responses around 2.0 ns because the direct waves and ground reflections (overlapped between 0.2 ns and 1.8 ns) are dominant in the B-scan (Figure 6a). The imaging result of the B-scan is displayed in Figure 6b where an ambiguous dot appears at 15 cm location. Though focused and enhanced by the imaging processing, the images of target in sub-band images Figure 6c,d are still obscured by the clutter. Processing the sub-band images with our method yields much more clearer target image (Figure 6e) that are brighter and larger than the remaining clutter. Comparing Figure 6e and Figure 6b, the SCR is improved by 26.3 dB. After applying thresholds, clutter residues are eliminated and the target is detected, as shown in Figure 6f. The results of competitive methods are shown in Figure 6g,h. Since the target is small and contains no metal, the responses are too weak to be detected by the PCA-CFAR method and the RPCA-based method. In addition, clutter residues in Figure 6g cause multiple false alarm in the detection results.

In the other B-scan, there are two shallowly buried (3 cm) APMs containing both metal and nonmetal. These APMs are identical and their heights and radiuses are 5.4 cm and 6 cm, respectively. The size and location of the APMs are indicated by dash-dotted boxes. The reflections from metal are strong enough to be observed in Figure 7a. However, the target image in Figure 7b is overlapped with strong ground clutter and the SCR is insufficient to confirm the target. Similar are the sub-band images shown in Figure 7c,d, where the target images are of small magnitude compared with the ground reflections. Note that more details of the target can be found in the HF image Figure 7d than in LF image Figure 7c, which demonstrates the frequency dependency of target responses in different sub-bands. Applying our method to this GPR data, we obtain the prominent target image with a much higher SCR (30 dB) in Figure 7e. Then, the targets’ locations are successfully detected and marked by points in Figure 7f. Figure 7g,h display the results of the PCA-CFAR method and the RPCA-based method, and the SCR of these images are 12.6 dB and 19.8 dB, respectively.

As listed in Table 3, the performance of our method on this dataset outperforms the PCA-CFAR with more accurate detection (100% PD) and less false alarms (12.3% FAR), whereas the PCA-CFAR method yields a higher FAR (54.5%) to achieve 89.4% PD. The performance of RPCA-based method is of 93.3% PD and 10.8% FAR.

The second experiment was conducted on the dataset published by the Georgia Technology Institute. Seven types of mines were buried in a 1.8 m × 1.8 m region. The test site was different from the actual minefield because the targets were close to each other, but this dataset can still be used to demonstrate the applicability of our method to the SFCW GPR data. The experimental setup is detailed in [37,38]. Fifty B-scans acquired along the x axis were selected to test our method. As shown in the burial map Figure 8a, the coverage area of these B-scans is enclosed by a rectangle, and a dash line across the center indicates location of the example B-scan. The locations and sizes of these four target are indicated by four dash-dotted boxes. As we can see in the imaging result (Figure 8b), the horizontal stripe-like clutter is overlapped with the images of targets. Two deep buried targets can be found due to their big sizes and strong reflections while the other two small targets buried shallowly near the center are hidden in the clutter. In the LF image Figure 8c, two targets of the bigger sizes can be noticed, while the smallest simulant mine is prominent in the HF image Figure 8d. Removing clutter using our method, the result is shown in Figure 8e, where the SCR is increased to 26.8 dB. Finally, all of the four targets are correctly detected and indicated with points in Figure 8f after the thresholds are applied. After processing the data with the PCA-CFAR method and the RPCA-based method, SCRs are improved by 11 dB and 27.9 dB, respectively. Also, the targets are detected by these two methods, as shown in Figure 8g,h. However, there is a false alarm indicated by a triangle in Figure 8h, which may be caused by the residual ground clutter.

After processing this dataset, our proposed method detects all of the targets shown in Figure 8a with an FAR of 7.0%, while the PCA-CFAR miss two small targets and the FAR is 19.2%. The RPCA-based method achieves 100% PD and 8.4% FAR. The performances are summarized in Table 3.

4.4. Discussion

• Reconstruction of the target image. As a part of our proposed method, GPR data is migrated using the LRM algorithm. The resolution of the signal is not enough for an accurate shape reconstruction of our mine-like targets but improves the SCR of the image, making the targets look like points with high energy.
• Overlap of sub-band. The idea of improving performance by allowing some overlap comes form [27]. According to our investigation on the synthesized dataset, the best performance is obtained with an overlap close to 50–60% (see Figure 9).

5. Conclusions

In this paper, we propose a novel TRPCA-based method for clutter suppression and target detection. A tensor formation scheme capturing the target sub-band features is firstly presented. Then, taking advantage of the low-rank and sparse property of the multidimensional tensor, our method formulates the clutter removal and target detection as a TRPCA problem and then solves the problem by decomposing the obtained tensor with a tractable method—ADMM. The feasibility and effectiveness of our proposed method for APM detection applications are verified by both numerical simulations and laboratory experiments with different GPR systems. Our proposed method is distinguished from traditional ones by the advantage of multidimensional data structure, which enables exploiting of the information contained in different sub-bands. Since the sub-band technique is a rough description of the frequency-dependent features, our future works will focus on designing a precise characterization of target features to guide the tensor construction and applying classification to discard false positives. Considering the conditions of laboratory experiments are different from those of the actual minefields, the analysis should be extended to field experiments under realistic clutter conditions.