# Time delay estimation by a modified orthogonal matching pursuit method for rough pavement

Jingjing Pan, Meng Sun, Yide Wang, Cédric Le Bastard, Vincent Baltazart

► To cite this version:

Jingjing Pan, Meng Sun, Yide Wang, Cédric Le Bastard, Vincent Baltazart. Time delay estimation by a modified orthogonal matching pursuit method for rough pavement. IEEE Transactions on Geoscience and Remote Sensing, 2021, 59 (4), pp.2973-2981. 10.1109/TGRS.2020.3006509. hal-02877473

HAL Id: hal-02877473

https://hal.science/hal-02877473

Submitted on 25 Oct 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Time-Delay Estimation by a Modified Orthogonal Matching Pursuit Method for Rough Pavement

Jingjing Pan<sup>ID</sup>, Member, IEEE, Meng Sun<sup>ID</sup>, Member, IEEE, Yide Wang<sup>ID</sup>, Senior Member, IEEE, Cédric Le Bastard<sup>ID</sup>, and Vincent Baltazart

**Abstract**—Pavement survey is one of the most important applications for ground penetrating radar (GPR) in civil engineering. In the case of centimeter scale of GPR waves, the influence of interface roughness cannot be neglected and should be taken into account in the radar data model. The objective of this article is to estimate the time-delay in the presence of interface roughness by GPR. Using the property of noncircular signals, we propose a modified orthogonal matching pursuit method to estimate the pavement parameters for both overlapped and nonoverlapped echoes. Compared with subspace-based methods in coherent scenarios, the proposed method can estimate the time delays of backscattered echoes without applying the cumbersome interpolation and spatial smoothing procedures, which are more practical in real applications. The performance of the proposed method is tested on both simulated and experimental data. The estimation results show the good performance of the proposed method.

**Index Terms**—Ground penetrating radar (GPR), interface roughness, orthogonal matching pursuit (OMP), time-delay estimation (TDE).

## I. INTRODUCTION

THE research on electromagnetic (EM) scattering by layered interface has a large number of applications, such as pavement survey, underground exploration, and environmental and agricultural monitoring applications [1]–[5]. In pavement survey, it can be used to measure the thickness of layer, and to detect and evaluate the damaged zones (interface debonding of pavements and of seal coats of highway structures) [6]–[8]. EM waves with centimeter-scale wavelengths can be used for measuring the thickness of different layers [9]–[14].

Manuscript received February 13, 2020; revised May 30, 2020; accepted June 19, 2020. Date of publication July 17, 2020; date of current version March 25, 2021. This work was supported in part by the Shanghai Sailing Program under Grant 19YF1419100 and in part by the Key Laboratory of Dynamic Cognitive System of Electromagnetic Spectrum Space (Nanjing University of Aeronautics and Astronautics), Ministry of Industry and Information Technology under Grant KF20202101. The GPR data were collected at UGE/IFSTTAR within the scope of the research project ACIMP (ANR-18-CE22-0020) supported by the French National Research Agency (ANR). (Corresponding author: Meng Sun.)

Jingjing Pan and Meng Sun are with the Department of Electronic Engineering, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China and also with the Key Laboratory of Dynamic Cognitive System of Electromagnetic Spectrum Space (Nanjing University of Aeronautics and Astronautics), Ministry of Industry and Information Technology, Nanjing 210016, China (e-mail: ee_msun@foxmail.com).

Yide Wang is with the Institut d'Électronique et des Technologies du numérique (IETR), UMR CNRS 6164, Polytech Nantes-Site de la Chantrerie, 44306 Nantes, France.

Cédric Le Bastard is with personal research.

Vincent Baltazart is with the UGE/IFSTTAR/COSYS, 44341 Nantes, France.

Color versions of one or more of the figures in this article are available online at https://ieeexplore.ieee.org.

Digital Object Identifier 10.1109/TGRS.2020.3006509

As a comment tool of nondestructive testing method, ground penetrating radar (GPR) is applied in the measurement, which permits rapid data collection for pavement surveys. In pavement survey, the pavement layers are assumed to be horizontally stratified; the vertical structure of the pavement can be deduced from radar profiles by mean of time-delay and amplitude estimations.

As mentioned in [15], within the scope of centimeter scale of EM waves, or for the ultra wide band (UWB) GPR, with upper frequency ranges up to 8–10 GHz, the influence of interface roughness cannot be neglected. The interface roughness leads to a decrease of echoes' amplitude with frequency, which can be characterized by the frequency behavior of backscattered echoes. The interface roughness measured by GPR is very important for road safety and can be applied for pavement skid resistance estimation. Therefore, this frequency behavior should be considered in the radar data model. Moreover, many signal processing methods based on this model have been developed for time-delay estimation in the presence of interface roughness [15]–[18].

In the case of narrow frequency bands, namely less than 2 GHz, the frequency behavior of the backscattered echoes coming from interface roughness can be approximated as an exponential function [15]. High-resolution methods like MUSIC and ESPRIT combined with spatial smoothing-based techniques can easily be adapted for estimating the pavement parameters like layer thickness, permittivity, and interface roughness [17], [18]. However, the assumption of the exponential function is not valid with the widening of the frequency bandwidth [15].

It has been mentioned in [15] that in large frequency bands ($B > 2$ GHz), the frequency behavior of backscattered echoes is more appropriately approximated as a Gaussian function. In this situation, the conventional spatial smoothing technique cannot be applied, due to the nonlinear frequency behavior in the exponent of the exponential function of echoes [15]. In the work of [15] and [16], an interpolated spatial smoothing technique is applied. First, the interpolation is used to transform the frequency behavior of backscattered echoes into a linear exponential frequency behavior; and then the conventional spatial smoothing technique is used to reduce the correlation between echoes. The interpolation can take several possible frequency behaviors into account, which is more applicable for UWB-GPR. Afterward, the modified MUSIC and ESPRIT are adapted for parameter estimation. The interpolation procedure requires to build a database of transformation matrices to record
the transformations between the true frequency behavior and interpolated one for different pavement geometries. Each pavement geometry corresponds to a unique transformation matrix. Additionally, this transformation matrix depends on many factors (geometric parameters), such as the number of layers, interface roughness, thickness, and permittivities, which are very complicated to be acquired. On the other hand, select the good transformation matrix is another difficult task, due to the unknown frequency behavior. Besides, there is no uniform rules for the transformation matrix selection [16], [19], [20]. Due to the complexity in the GPR measurement, the interpolation may be difficult to be applied in practice.

Moreover, in the past decade, methods based on compressive sensing (CS) have been proposed for parameter estimation that can handle the coherent signals without decorrelation procedure [21], [22]. For example, the matching pursuit (MP) and orthogonal matching pursuit (OMP) [23], [24] are the two most popular classical greedy algorithms for sparse reconstruction, which can be used for time-delay estimation. However, their capability to distinguish two echoes is restricted by the GPR bandwidth [13]. Similar to the classical fast Fourier transform (FFT)-based methods, the conventional MP and OMP methods have a limited resolution power in the scenario of overlapped echoes ($B \Delta \tau \le 1$) [25].

In this context, the aim of this article is to propose a new time-delay estimation method for rough pavement using GPR, especially for overlapped backscattered echoes. Unlike the previous subspace-based methods in the presence of interface roughness [15], [16], the proposed method does not require additional interpolation and spatial smoothing procedures. The objectives of this work are as follows.

1) Improve the time resolution of GPR by exploiting the noncircularity of backscattered echoes using the extended signal model. Within lossless pavements like in [13] and [15], the backscattered echoes depend on the reflection coefficients of the media. According to [26], the reflection coefficients in lossless media are real and therefore the echoes have the noncircular property, like amplitude modulated (AM) or binary phase shift keying (BPSK) modulated signals. By concatenating the GPR measurements and their conjugate components, the size of the observed data is doubled [27], which contributes to the increase of time resolution.

2) Adapt OMP for time-delay estimation in the presence of interface roughness. There are more than one unknown parameters, including time-delay, that indicates interface roughness. The modified OMP separates the time-delay parameter from other unknown parameters in the estimation.

The rest of this article is organized as follows: Section II presents the radar data model and its extension using the property of noncircular signals. In Section III, we present the modified OMP method for time-delay estimation. Sections IV and V provide, respectively, the simulation and experimental results to show the effectiveness of the proposed method. Conclusions and perspectives are drawn in Section VI.

## II. RECEIVED SIGNAL MODEL

Similar to the work in [15] and [16], the first two or three layers of pavement are considered in this article. Therefore, the studied media is regarded as low-loss media [28]. As mentioned in [26], for low-loss media, the dispersivity of the media is negligible. The used radar data model represents the echoes backscattered from the media separated by rough interfaces. The model in [15] is adopted in this article, which can be written as

$$r(f) = \sum_{k=1}^{K} e(f) s_k(f) \exp(-j2\pi ft_k) + n(f) \quad (1)$$

where

1) $K$ is the number of backscattered echoes, which is assumed to be known or can be estimated by some detection criteria [29], [30];

2) $e(f)$ represents the radar pulse in frequency domain;

3) $s_k(f)$ represents the amplitude of the $k$th backscattered echo, which can be divided into two parts: $s_k(f) = s_k \phi_k(f)$;

4) $s_k$ is the reflection coefficient of the $k$th backscattered echo which is independent of frequency $f$;

5) $\phi_k(f)$ represents the frequency behavior of the $k$th backscattered echo;

6) $n(f)$ is an additive white Gaussian noise with zero mean and variance $\sigma^2$.

For $N$ discrete frequencies within bandwidth $B$, (1) can be written as

$$\mathbf{r} = \Lambda \mathbf{A} \mathbf{s} + \mathbf{n} \quad (2)$$

with the following notational definitions:

1) $\mathbf{r} = [r(f_1) \ r(f_2) \ \dots \ r(f_N)]^T$ is the $(N \times 1)$ received signal vector representing the measurements by a step-frequency radar; $f_i = f_1 + (i - 1)\Delta f$ and $i = 1, 2 \dots N$ ($f_1$ is the lowest frequency and $\Delta f$ is the frequency step); superscript $T$ denotes the transpose operation.

2) $\Lambda = \text{diag}\{e(f_1), e(f_2), \dots, e(f_N)\}$ is a $(N \times N)$ diagonal matrix, whose diagonal elements are the radar pulse in frequency domain.

3) $\mathbf{A} = [\mathbf{\Phi}_1 \mathbf{a}(t_1) \ \mathbf{\Phi}_2 \mathbf{a}(t_2) \ \dots \ \mathbf{\Phi}_K \mathbf{a}(t_K)]$ is the mode matrix.

4) $\mathbf{\Phi}_k \mathbf{a}(t_k)$ is the mode vector of the $k$th echo.

5) $\mathbf{a}(t_k) = [e^{-2j\pi f_1 t_k} \ e^{-2j\pi f_2 t_k} \ \dots \ e^{-2j\pi f_N t_k}]^T$.

6) $\mathbf{\Phi}_k = \text{diag}\{\phi_k(f_1), \phi_k(f_2), \dots, \phi_k(f_N)\}$ is a $(N \times N)$ diagonal matrix, whose diagonal elements represent the frequency behavior of the $k$th backscattered echo; $\phi_k(f_i)$ is real.

7) $\mathbf{s} = [s_1 \ s_2 \ \dots \ s_K]^T$ is the $(K \times 1)$ vector of the reflection coefficients of media.

8) $\mathbf{n} = [n(f_1) \ n(f_2) \ \dots \ n(f_N)]^T$ is the $(N \times 1)$ noise vector with zero mean and covariance matrix $\sigma^2 \mathbf{I}$, $\mathbf{I}$ is the $(N \times N)$ identity matrix.

In the following, the data whitening procedure is applied. Divided by the radar pulse, the whitened received signal vector $\mathbf{r}_w$ can be written as

$$\mathbf{r}_w = \Lambda^{-1} \mathbf{r} = \mathbf{A} \mathbf{s} + \Lambda^{-1} \mathbf{n} = \mathbf{A} \mathbf{s} + \mathbf{b} \quad (3)$$
where **b** is the new noise vector after the data whitening procedure.

### III. METHODOLOGY

#### A. Extended Signal Model

In (2) and (3), vector **s** represents the reflection coefficients of media, whose elements are real. Therefore, the backscattered echoes can be considered as noncircular signals. Using the property of noncircular signals, we concatenate the radar measurements and their conjugate components to obtain an extension of the received radar data model as follows:

$$ 
\begin{aligned}
\mathbf{r}_{\text{nc}} &= \begin{pmatrix} \mathbf{J} \mathbf{r}_{\text{w}}^* \\ \mathbf{r}_{\text{w}} \end{pmatrix} = \begin{pmatrix} \mathbf{J} \mathbf{A}^* & \mathbf{0} \\ \mathbf{0} & \mathbf{A} \end{pmatrix} \begin{pmatrix} \mathbf{s}^* \\ \mathbf{s} \end{pmatrix} + \begin{pmatrix} \mathbf{J} \mathbf{b}^* \\ \mathbf{b} \end{pmatrix} \\
&= \begin{pmatrix} \mathbf{J} \mathbf{A}^* \\ \mathbf{A} \end{pmatrix} \mathbf{s} + \begin{pmatrix} \mathbf{J} \mathbf{b}^* \\ \mathbf{b} \end{pmatrix}
\end{aligned} 
\quad (4) $$

where matrix **J** is the $(N \times N)$ anti-identity matrix (exchange matrix) and the operator * denotes the complex conjugate. Therefore, the new mode vector of this extended radar data model $\mathbf{\Phi}_{\text{nc}}\mathbf{a}_{\text{nc}}(t)$ can be expressed as

$$ 
\begin{aligned}
\mathbf{\Phi}_{\text{nc}}\mathbf{a}_{\text{nc}}(t) &= \begin{pmatrix} \mathbf{J}\mathbf{\Phi}^* & \mathbf{0} \\ \mathbf{0} & \mathbf{\Phi} \end{pmatrix} \begin{pmatrix} \mathbf{a}^*(t) \\ \mathbf{a}(t) \end{pmatrix} \\
&= [e^{2j\pi f_N t}\phi^*(f_N) \dots e^{2j\pi f_1 t}\phi^*(f_1) \\
& \qquad \qquad \qquad e^{-2j\pi f_1 t}\phi(f_1) \dots e^{-2j\pi f_N t}\phi(f_N)]^T \\
&= [e^{2j\pi f_N t}\phi(f_N) \dots e^{2j\pi f_1 t}\phi(f_1) \\
& \qquad \qquad \qquad e^{-2j\pi f_1 t}\phi(f_1) \dots e^{-2j\pi f_N t}\phi(f_N)]^T \quad (5)
\end{aligned}
$$

with

$$ \mathbf{a}_{\text{nc}}(t) = [e^{2j\pi f_N t} \dots e^{2j\pi f_1 t} e^{-2j\pi f_1 t} \dots e^{-2j\pi f_N t}]^T $$

and

$$ \mathbf{\Phi}_{\text{nc}} = \text{diag}\{\phi(f_N), \dots, \phi(f_1), \phi(f_1), \dots, \phi(f_N)\}. $$

It can be seen from (5) that this new signal model has increased the used frequency band (the frequency band **B** is doubled, which is equal to $2(N-1)\Delta f$) [27] without increasing the physical frequency sampling points. In consequence, the time resolution of the GPR system can be significantly improved using this extended model.

#### B. Modified Orthogonal Matching Pursuit Method

In real measurement, the backscattered echoes are highly correlated, and their frequency behavior can be of various shapes, whose exponent may be a nonlinear function of frequency [15]. Subspace-based methods cannot be applied without interpolated spatial smoothing technique. However, the interpolation procedure is complicated due to the complex condition in the measurement. To solve the abovementioned problems, in this section, a modified OMP method is introduced for time-delay estimation without either interpolation or spatial smoothing technique.

According to the noncircularity of the backscattered echoes, the signal model can be reformulated as

$$ 
\begin{aligned}
\mathbf{r}_{\text{nc}} &= [\mathbf{B}(t_1), \mathbf{B}(t_2), \dots, \mathbf{B}(t_K)] \begin{pmatrix} \psi_1 \\ \psi_2 \\ \vdots \\ \psi_K \end{pmatrix} + \begin{pmatrix} \mathbf{J}\mathbf{b}^* \\ \mathbf{b} \end{pmatrix} \\
&= \sum_{k=1}^{K} \mathbf{B}(t_k)\psi_k + \mathbf{b}_{\text{nc}} \quad (6)
\end{aligned}
$$

with

$$ \mathbf{b}_{\text{nc}} = \begin{pmatrix} \mathbf{J}\mathbf{b}^* \\ \mathbf{b} \end{pmatrix} $$

$$ \mathbf{B}(t_k) = \text{diag}\{e^{2j\pi f_N t_k}, \dots, e^{2j\pi f_1 t_k}, e^{-2j\pi f_1 t_k}, \dots, e^{-2j\pi f_N t_k}\} $$

and

$$ \psi_k = s_k[\phi_k(f_N) \phi_k(f_{N-1}) \dots \phi_k(f_1) \phi_k(f_1) \dots \phi_k(f_N)]^T. $$

A pre-estimation of the time-delays is carried out using the conventional OMP or beamforming method, which is able to roughly estimate the time-delays before the application of the proposed modified OMP method; the range of time domain can then be set according to the pre-estimated results. The pre-estimation method contributes to the reduction of the size of the overcomplete dictionary and the computational complexity of the proposed method [31]. Afterward, the time domain is sampled as $T = [\tau_1 \tau_2 \dots \tau_{N_0}]$, with $N_0 \gg K$ to build the dictionary matrix. Accordingly, the $(2N \times 2NN_0)$-dimensional overcomplete dictionary is shown as $\mathbf{B}_T = [\mathbf{B}(\tau_1), \dots, \mathbf{B}(\tau_i), \dots, \mathbf{B}(\tau_{N_0})]$.

The OMP method is a greedy compressed sensing recovery algorithm. The principle is to select the best-fitting atom of the dictionary matrix in each iteration, and then a least squares (LS) optimization is performed in the subspace spanned by all the previously picked columns. Moreover, the number of iterations is equal to the number of backscattered echoes. A short description of OMP method is shown in the appendix.

Similar to the conventional OMP, at each iteration, the objective of the proposed modified OMP is to select the best-fitting atom from the overcomplete dictionary, which is the most correlated with the residual signal. The position of this atom in the dictionary corresponds to the estimated time delay. At the $k$th iteration, the correlation between the residual signal $\mathbf{r}_{\text{nc}}^{k-1}$ (the initial residual signal is set to be $\mathbf{r}_{\text{nc}}^0 = \mathbf{r}_{\text{nc}}$) and the $i$th atom of the updated overcomplete dictionary matrix $\mathbf{B}(\tau_i)$ can be defined as the absolute value of the sum of the elements of $\mathbf{B}^H(\tau_i)\mathbf{r}_{\text{nc}}^{k-1}$:

$$ \mu_k(i) = |\mathbf{e}\mathbf{B}^H(\tau_i)\mathbf{r}_{\text{nc}}^{k-1}| \quad (7) $$

where $\mathbf{e} = [1, \dots, 1]$ is a $(1 \times 2N)$ vector. The atom with the highest correlation ($\max_i \mu_k(i)$) is the best-fitting atom. Besides, the elements of $\psi_k$ ($k = 1, 2, \dots, K$) corresponding to the amplitudes of the backscattered echoes (varying with frequency) are real. When $\tau_i = t_k$, the elements of $\mathbf{B}^H(\tau_i)\mathbf{r}_{\text{nc}}^{k-1}$ representing the estimated amplitude of the $k$th backscattered echo should be real in theory. Therefore, at each iteration, to find the best-fitting atom with the highest correlation
with the residual signal, only the real part of $\mathbf{B}^H(\tau_i)\mathbf{r}_{nc}^{k-1}$ is considered in the comparison. We can rewrite (7) as

$$\mu_k(i) = |\text{real}\{\mathbf{e}\mathbf{B}^H(\tau_i)\mathbf{r}_{nc}^{k-1}\}|.$$

Obviously, (7) returns the absolute value of the sum of the real part of vector $\mathbf{B}^H(\tau_i)\mathbf{r}_{nc}^{k-1}$'s components. The time of arrival $\hat{t}_k$ can then be retrieved by searching the best-fitting atom.

The second step is to calculate a new residual signal for the $(k+1)$ iteration. To start, we calculate the amplitude of the $k$th echo. As mentioned in [15], the echoes' amplitudes change with frequency due to interface roughness, which are not constant. To estimate the following residual signal, we calculate the mean of amplitude of $k$th echo using LS optimization as follows:

$$\hat{s}_k = \mathbf{a}_{nc}^{+}(\hat{t}_k)\mathbf{r}_{nc}^{k-1} \quad (8)$$

where operator + is the Moore-Penrose transpose, $\mathbf{a}_{nc}^{+}(\hat{t}_k) = (\mathbf{a}_{nc}^H(\hat{t}_k)\mathbf{a}_{nc}(\hat{t}_k))^{-1}\mathbf{a}_{nc}^H(\hat{t}_k)$. Then, the residual signal at the $(k+1)$th iteration can be written as

$$\mathbf{r}_{nc}^{k} = \mathbf{r}_{nc}^{k-1} - \mathbf{a}_{nc}(\hat{t}_{k})\hat{s}_{k}. \quad (9)$$

A summary of the proposed modified OMP method is shown in Table I. Using the proposed modified OMP method, the time-delays of backscattered echoes are estimated. Afterward, we apply maximum-likelihood estimation (MLE) method for interface roughness estimation with the estimated time-delays. This procedure has been well discussed, the details can be found in [15].

## IV. SIMULATION EXPERIMENT

In this section, the simulation results of the proposed modified OMP are presented. In the simulation, we only study the results of time-delay estimation. As previously mentioned, the interface roughness estimation by MLE has been well discussed in [15], which is not the focus of this article. The simulated data are provided by the method of moments, which represents the GPR backscattered echoes at nadir. It can be seen from Fig. 1 that the backscattered echoes are reflected from a rough pavement made up of two uncorrelated random rough interfaces separating homogeneous media. The studied pavement has two layers: the first layer is an ultra-thin asphalt surfacing (UTAS), whose relative permittivity ($\varepsilon_{r2}$) is about 4.5; the second layer is a base layer with relative permittivity $\varepsilon_{r3} = 7$. The thickness of the UTAS is approximately 15 mm, therefore, the times of arrival of the backscattered echoes $s_1$ and $s_2$ (see Fig. 1) from UTAS are 1 and 1.21 ns, respectively. The interfaces $\Sigma_A$ and $\Sigma_B$ are considered to follow a Gaussian height probability density function and an exponential height autocorrelation function. The interface roughness can then be defined by two parameters: the root-mean square height $\sigma_h$, correlation length $L_c$. In the first simulation, three different rough pavements are studied as follows.

1) Case I: $\sigma_{hA} = 1.0$ mm, $L_{cA} = 6.4$ mm, $\sigma_{hB} = 2.0$ mm, $L_{cB} = 15$ mm, lossless media.

2) Case II: $\sigma_{hA} = 1.5$ mm, $L_{cA} = 6.4$ mm, $\sigma_{hB} = 2.5$ mm, $L_{cB} = 15$ mm, lossless media.

<table><caption>TABLE I
PROCEDURE OF MODIFIED OMP METHOD</caption><tr><td colspan="2"><b>INPUT:</b></td></tr><tr><td colspan="2">(2N × 2NN<sub>0</sub>)-dimensional overcomplete dictionary B<sub>T</sub></td></tr><tr><td colspan="2">(2N × 1)-dimensional received signal vector r<sub>nc</sub></td></tr><tr><td colspan="2">The number of backscattered echoes K</td></tr><tr><td colspan="2"><b>OUTPUT:</b></td></tr><tr><td colspan="2">Estimated time of arrival t&#x0302;<sub>k</sub>, k = 1, 2, ..., K</td></tr><tr><td colspan="2"><b>INITIALIZATION:</b></td></tr><tr><td>Initial residual signal r<sub>nc</sub><sup>0</sup> = r<sub>nc</sub></td><td></td></tr><tr><td>Initial iteration step k = 1</td><td></td></tr><tr><td>Initial index set Λ<sup>0</sup> = ∅</td><td></td></tr><tr><td>Initial matrix of chosen atoms B<sub>0</sub> = ∅</td><td></td></tr><tr><td colspan="2"><b>while</b> not converged (k &lt; K) do</td></tr><tr><td colspan="2"><b>MATCH:</b></td></tr><tr><td>μ<sub>k</sub>(i) = |real{eB<sup>H</sup>(τ<sub>i</sub>)r<sub>nc</sub><sup>k-1</sup>}|, i = 1, 2, ..., N<sub>0</sub></td><td></td></tr><tr><td colspan="2">Find the index t&#x0302;<sub>k</sub> that solves the optimization problem:</td></tr><tr><td>t&#x0302;<sub>k</sub> = arg max<sub>i</sub> {μ<sub>k</sub>(i)}</td><td></td></tr><tr><td colspan="2"><b>UPDATE:</b></td></tr><tr><td>Amplitude of the kth echo s&#x0302;<sub>k</sub> = a<sub>nc</sub><sup>+</sup>(t&#x0302;<sub>k</sub>)r<sub>nc</sub><sup>k-1</sup></td><td></td></tr><tr><td>Residual signal r<sub>nc</sub><sup>k</sup> = r<sub>nc</sub><sup>k-1</sup> - a<sub>nc</sub>(t&#x0302;<sub>k</sub>)s&#x0302;<sub>k</sub></td><td></td></tr><tr><td>Iteration step k = k + 1</td><td></td></tr><tr><td>Index set Λ<sup>k</sup> = Λ<sup>k-1</sup> ∪ t&#x0302;<sub>k</sub></td><td></td></tr><tr><td>Matrix of chosen atoms B<sub>k</sub> = [B<sub>k-1</sub>, B(t&#x0302;<sub>k</sub>)]</td><td></td></tr><tr><td colspan="2"><b>end while</b></td></tr></table>

<figure><img src="image_2.png" /><figcaption>Fig. 1. Rough pavement configuration.</figcaption></figure>

3) Case III: $\sigma_{hA} = 1.0$ mm, $L_{cA} = 6.4$ mm, $\sigma_{hB} = 2.0$ mm, $L_{cB} = 15$ mm, low-loss media (conductivity $\rho_A = 5 \times 10^{-3}$ S/m, conductivity $\rho_B = 10^{-2}$ S/m).

In the simulations, the proposed modified OMP method is tested in the frequency band 0.5-3.5 GHz for overlapped echoes, with 0.1 GHz frequency step (31 frequency samples). Therefore, the backscattered echoes are overlapped with $B \Delta \tau = 0.63$. The number of independent snapshots is 100. The signal-to-noise ratio (SNR) is defined as the ratio between the powers of the first echo and the noise variance.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Case I: Time-delay estimation (TDE) with SNR= 0 dB, the red dashed lines represent the true times of arrival, nonoverlapped echoes.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Case I: TDE with SNR = 0 dB, the red-dashed lines represent the true times of arrival, overlapped echoes.</figcaption></figure>

<table><caption>TABLE II
ESTIMATED TIMES OF ARRIVAL BY MODIFIED OMP, $\hat{t}_k$ REPRESENTS THE ESTIMATED TIME OF ARRIVAL OF THE kTH ECHO</caption><thead><tr><th rowspan="2">Time</th><th colspan="2">Modified OMP</th><th colspan="2">OMP</th></tr><tr><th>t̂<sub>1</sub> (ns)</th><th>t̂<sub>2</sub> (ns)</th><th>t̂<sub>1</sub> (ns)</th><th>t̂<sub>2</sub> (ns)</th></tr></thead><tbody><tr><td>case a non-overlapped</td><td>1.000</td><td>1.210</td><td>0.998</td><td>1.206</td></tr><tr><td>case a overlapped</td><td>1.000</td><td>1.208</td><td>0.972</td><td>1.312</td></tr><tr><td>case b overlapped</td><td>1.000</td><td>1.210</td><td>0.975</td><td>1.310</td></tr><tr><td>case c overlapped</td><td>0.998</td><td>1.214</td><td>0.970</td><td>1.322</td></tr></tbody></table>

## A. Performance Under Different Rough Pavements

In the first simulation, SNR is fixed at 0 dB for the different rough pavements. Moreover, the results of the conventional OMP method are taken as a comparison. To begin with, a nonoverlapped situation (the frequency bandwidth 0.5–6.5 GHz, $B\Delta\tau = 1.26$) is tested with 31 frequency samples for Case I. Fig. 2 and Table II show the estimated times of arrival of the nonoverlapped backscattered echoes $s_1$ and $s_2$. It can be seen that for nonoverlapped echoes, time-delay is well estimated by both the proposed modified OMP and conventional OMP. However, the conventional OMP method suffers performance degradation with overlapped echoes. In the following, Figs. 3–5 and Table II provide the estimated times of arrival of the overlapped backscattered echoes $s_1$ and $s_2$ (the frequency bandwidth 0.5–3.5 GHz) by the proposed modified OMP and conventional OMP for different rough pavements. As shown in Figs. 3–5, the proposed method can detect the true time-delays of the backscattered echoes for the three different cases (two peaks corresponding to the estimated times of arrival are clearly detected by the proposed method). The times of arrival of $s_1$ and $s_2$ are well estimated by the proposed method with a small error, see Table II. However, the conventional OMP fails in time-delay estimation, due to the limited resolution power. From the results of the first simulation, we can conclude that the proposed modified OMP can be applied on either lossless or low-loss media.

In addition, a more generalized situation (Case IV) is considered: three backscattered echoes with pavement made up of three rough interfaces (three layers). The simulation parameters are shown in Table III. The roughness parameters are obtained from the received signal model in (1) with frequency behavior being Gaussian function [15]. Clearly, the first two echoes are overlapped with each other ($B\Delta\tau = 0.84$); and the third echo is nonoverlapped with other echoes.

<figure><img src="image_5.png" /><figcaption>Fig. 4. Case II: TDE with SNR = 0 dB, the red-dashed lines represent the true times of arrival, overlapped echoes.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Case III: TDE with SNR = 0 dB, the red-dashed lines represent the true times of arrival, overlapped echoes.</figcaption></figure>

<table><caption>TABLE III
SIMULATION PARAMETERS APPLIED IN CASE IV</caption><thead><tr><th>Parameter</th><th>1</th><th>2</th><th>3</th></tr></thead><tbody><tr><td>Permittivity of layer</td><td>4.5</td><td>7.0</td><td>10.0</td></tr><tr><td>Thickness of layer (cm)</td><td>2.0</td><td>2.5</td><td>∞</td></tr><tr><td>Time of arrival (ns)</td><td>1.00</td><td>1.28</td><td>1.72</td></tr><tr><td>Roughness parameter (GHz<sup>-2</sup>)</td><td>0.0016</td><td>0.0169</td><td>0.0235</td></tr></tbody></table>

Fig. 6 shows the simulation results of the time-delay estimation by both the proposed method and conventional OMP for three echoes. The conventional OMP fails to detect the true time-delays of echoes. Nevertheless, the proposed method remains robust in this situation. Obviously, the three peaks corresponding to the times of arrival of the three echoes are well estimated with $\hat{t}_1 = 0.998$ ns, $\hat{t}_2 = 1.293$ ns, and $\hat{t}_3 = 1.730$ ns. It can be concluded that the proposed method can be applied to both overlapped and nonoverlapped echoes.
<figure><img src="image_7.png" /><figcaption>Fig. 6. Case IV: TDE with SNR= 0 dB, the red-dashed lines represent the true times of arrival.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. RRMSE on the estimated time-delay $\hat{\Delta\tau}$ versus the product $B\Delta\tau$, with 100 Monte-Carlo simulations.</figcaption></figure>

## B. Performance Versus $B\Delta\tau$

To determine the limitation of the proposed modified OMP in terms of the resolution power, in the second simulation, different data sets are generated by varying the layer thickness $H$ from 0.71 to 7.1 cm, the corresponding product $B\Delta\tau$ then changes from 0.3 to 3. The performance of the proposed method versus the product $B\Delta\tau$ is evaluated by a Monte-Carlo process of 100 independent runs. The relative root-mean-square error (RRMSE) of the studied parameter is defined as

$$RRMSE(x) = \frac{\sqrt{\frac{1}{M} \sum_{m=1}^{M} (\hat{x}_m - x)^2}}{x} \quad (10)$$

where $\hat{x}_m$ represents the estimated parameter for the $m$th run of the method, and $x$ the true value. In the simulation, parameter $x$ can represent either the time-delay $\Delta\tau$ ($t_2-t_1$) or the time of arrival $t_k$ ($k=1,2$). Only Case I is considered. The proposed method is also compared with the conventional OMP.

Fig. 7 plots the RRMSE on the estimated time-delay ($\hat{\Delta\tau}$) by the proposed method and conventional OMP. As expected, it can be seen that the RRMSE on $\hat{\Delta\tau}$ decreases with increasing of $B\Delta\tau$ for both the proposed method and conventional OMP. As shown in Fig. 7, the conventional OMP method has limitation when $B\Delta\tau \le 1$, which is restricted to identify overlapped echoes. The proposed modified OMP has no such limitation, which can solve the cases for both overlapped and nonoverlapped echoes. Moreover, the proposed method has better accuracy and outperforms the conventional OMP.

## C. Performance Versus SNR

In the third simulation, the performance of the proposed method versus SNR is studied with 500 independent runs

<figure><img src="image_9.png" /><figcaption>Fig. 8. RRMSE on the estimated time of arrival $\hat{t}_k$ ($k=1,2$) versus SNR, after 500 Monte-Carlo simulations.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. Experimental device.</figcaption></figure>

of the method for the smooth interfaces. Moreover, the Cramér-Rao bound (CRB) results are also provided [32]. The RRMSEs on the estimated times of arrival $\hat{t}_1$ and $\hat{t}_2$ are calculated. In this simulation, SNR varies from -10 to 5 dB. Fig. 8 shows the RRMSE on the estimated times of arrival by the proposed method. Similar to the second simulation, the RRMSE on $\hat{t}_k$ ($k=1,2$) is continuously decreasing when SNR increases. Moreover, the RRMSE also depends on the amplitude of the backscattered echo. Normally, under same condition, the echo with large amplitude (the first echo) has smaller RRMSE than that of the weak echo (the second echo).

## V. REAL DATA EXPERIMENT

In this section, the proposed modified OMP is tested on the experimental data [33]. An UWB step-frequency radar is used, which is composed of vector network analyzer (VNA) and a bistatic antenna device whose transmitter (Tx) and receiver (Rx) are close to each other and fixed at the same distance along the scanning direction, see Fig. 9. Tx and Rx are of the ETSA A5 antennas. The antennas are about 70.0 cm above the tested pavement satisfying the far-field condition. The radar frequency bandwidth ranges from 0.8 to 10.8 GHz, with 0.025 GHz frequency step (401 frequency samples).

In the experiment, we study a rough asphalt pavement, see Fig. 10, which is made of three rough interfaces separating media, as shown in Fig. 9. The studied pavement structure is made of an asphalt layers $\Sigma_A$ overlying a very thin sand layer $\Sigma_B$ with thickness less than 1 cm and $\Sigma_C$ is the base-layer. The permittivity of the medium $\Omega_3$ is about 5. The length of
<figure><img src="image_11.png" /><figcaption>Fig. 10. Studied pavement.</figcaption></figure>

the tested surface is about 50 cm with a sample step 1 cm (50 sample points). For each sample point, a single snapshot is carried out.

## A. Data Set

The antenna is moved slightly between various sample points to generate independent spatial measurements. Fig. 11(a) and (b) displays the raw experimental data for both A-scan and B-scan, respectively. The multiple waves as the wave coming from the test bed and the air wave between the Tx and Rx devices can be canceled by a time filter. The two peaks in Fig. 11(a) and (b) corresponds to the first three echoes backscattered from interfaces: the first peak represents the echo backscattered from the top surface (first echo) and the second peak corresponds to the backscattered echoes from the debonding layer with a very small embedded sand thickness (second and third echoes). It can be seen from both A-scan and B-scan that the first backscattered echoes is clearly visible; however, the echoes from the second and third interfaces of the sand layer are overlapping and cannot be distinguished in the time data. In the following, we focus on the second and third backscattered echoes, namely, the small thickness of sand layer. The first echo will be removed after it is detected by the proposed method.

## B. Preprocessing of the Data

Two preprocessing techniques are applied combined with the modified OMP method: time filtering and data whitening.

1. Time Filtering. The time filtering is used to remove the echoes outside the GPR working time window or interested region. For example, by applying time filtering, the multiple echoes and air wave can then be eliminated; moreover, the time filtering may also be used to remove the impact of the residual of detected echoes in time-delay estimation (the echo backscattered from surface layer), see Fig. 11(c).

2. Data Whitening. To apply the proposed algorithms, a whitening procedure by the pulse is necessary. The

<figure><img src="image_12.png" /><figcaption>Fig. 11. Backscattered echoes from studied pavement. (a) Raw GPR data (A-scan). (b) Raw GPR data (B-scan). (c) Time filtered data.</figcaption></figure>

radar pulse is measured as the backscattered echo from a metallic plane.

After the two preprocessing techniques, the proposed modified OMP is then applied for the time-delay estimation.

## C. Time-Delay and Thickness Estimation

In the experiment, the proposed method is tested in frequency bandwidth $B = [0.8, 5.8]$ GHz (201 samples) at the sixth sample point. The conventional OMP is taken as comparison. By applying the proposed method, the times of arrival of the second and third backscattered echoes are estimated. Table IV presents the estimated time-delay $\hat{\Delta}\tau$ of the second and third echoes. In addition, the layer thickness can be calculated by the following equation:

$$ \hat{H}_2 = \frac{c \hat{\Delta}\tau}{2\sqrt{\epsilon_{r3}}} $$

where $\hat{H}_2$ is the estimated thickness of the sand layer, c is the speed of light in vacuum, and $\epsilon_{r3}$ is the permittivity of medium $\Omega_3$. The estimated layer thickness is also shown in Table IV. Compared with the conventional OMP, the proposed method gives a relatively good performance on thickness estimation as the estimated thickness is close to the real value.
<table><caption>TABLE IV
ESTIMATED TIME-DELAY AND THICKNESS BY THE
MODIFIED OMP AND OMP</caption><thead><tr><th>Method<br>Parameter</th><th>Modified OMP</th><th>OMP</th></tr></thead><tbody><tr><td>Δτ̂</td><td>0.130 ns</td><td>0.180 ns</td></tr><tr><td>Ĥ<sub>2</sub></td><td>0.87 cm</td><td>1.21 cm</td></tr></tbody></table>

Moreover, the computational burdens of the conventional OMP and proposed method are $O(KN_0N)$ [34] and $O(4KN_0N^2)$, respectively. The computational load of the proposed method is a bit heavier than that of the conventional OMP method.

## VI. CONCLUSION

In this article, we propose a novel method for estimating the time-delays within rough pavements using GPR. The proposed method begins with a double extension of the received signal model of the pavement by exploiting the noncircularity of backscattered echoes. This extension greatly improves the GPR time resolution without increasing the physical frequency bandwidth. Then, based on the extended model, the OMP method is adapted for TDE in the presence of interface roughness. In theory, the proposed method can directly deal with coherent backscattered echoes and has improved time resolution. Unlike the prior subspace-based researches which take the interface roughness into account, the proposed method does not require any interpolation or spatial smoothing procedures. The performance of the proposed method is shown not only with numerical data but also with a field experiment. This article provides a new approach for the enhancement of the time resolution of GPR in the survey of stratified media. Besides, the modified OMP contributes to the theoretical development of 2-D or multidimensional nonlinear parameter estimation using compressed sensing techniques. Future work would focus on the reduction of computational complexity for the joint estimation of multidimensional parameters.

## APPENDIX

In this appendix, a brief introduction of the conventional OMP method is presented.

To start with, let us form a dictionary matrix such that the entire time domain can be sampled as $T = [τ_1 τ_2 ... τ_{N_0}]$, with $N_0 ≫ K$, and the $(N × N_0)$-dimensional overcomplete dictionary is shown as $A_T = [a(τ_1) a(τ_2) ... a(τ_{N_0})]$. Then, the description of the OMP method is given as follows:

**INPUT:**

$(N × N_0)$-dimensional overcomplete dictionary $A_T$

$(N × 1)$-dimensional received signal vector $\mathbf{r}_w$

The number of backscattered echoes $K$

**OUTPUT:**

Estimated time of arrival $\hat{t}_k$, $k = 1, 2, ..., K$

**INITIALIZATION:**

Initial residual signal $\mathbf{r}_w^0 = \mathbf{r}_w$

Initial iteration step $k = 1$

Initial index set $\Lambda^0 = \emptyset$

while not converged ($k < K$) do

**MATCH:**

$\mu_k(i) = \mathbf{a}^H(\tau_i)\mathbf{r}_w^{k-1}$, $i = 1, 2, ..., N_0$

Find the index $\hat{t}_k$ that solves the optimization problem:

$\hat{t}_k = \arg\max_i \{\mu_k(i)\}$

**UPDATE:**

Amplitude of the $k$th echo $\hat{s}_k = \mathbf{a}^+(\hat{t}_k)\mathbf{r}_w^{k-1}$

Residual signal $\mathbf{r}_w^k = \mathbf{r}_w^{k-1} - \mathbf{a}(\hat{t}_k)\hat{s}_k$

Iteration step $k = k + 1$

Index set $\Lambda^k = \Lambda^{k-1} \cup \hat{t}_k$

**end while**

## ACKNOWLEDGMENT

The authors would like to thank to Dr. J. M. Simonin and Dr. X. Dérobert for their contributions in the measurement experiment, the database is available to the GPR community at https://doi.org/10.25578/N0RXSQ.

## REFERENCES

[1] M. Sun, J. Pan, C. Le Bastard, Y. Wang, and J. Li, "Advanced signal processing methods for ground-penetrating radar: Applications to civil engineering," *IEEE Signal Process. Mag.*, vol. 36, no. 4, pp. 74–84, Jul. 2019.

[2] M. Ambrosanio, M. T. Bevacqua, T. Isernia, and V. Pascazio, "The tomographic approach to ground-penetrating radar for underground exploration and monitoring: A more user-friendly and unconventional method for subsurface investigation," *IEEE Signal Process. Mag.*, vol. 36, no. 4, pp. 62–73, Jul. 2019.

[3] M. T. Falconi and F. S. Marzano, "Weather radar data processing and atmospheric applications: An overview of tools for monitoring clouds and detecting wind shear," *IEEE Signal Process. Mag.*, vol. 36, no. 4, pp. 85–97, Jul. 2019.

[4] X. Liu et al., "Ground penetrating radar (GPR) detects fine roots of agricultural crops in the field," *Plant Soil*, vol. 423, nos. 1–2, pp. 517–531, Feb. 2018.

[5] A. Benedetto and L. Pajewski, *Civil Engineering Applications of Ground Penetrating Radar*. Springer, 2015.

[6] D.-H. Chen, F. Hong, W. Zhou, and P. Ying, "Estimating the hotmix asphalt air voids from ground penetrating radar," *NDT E Int.*, vol. 68, pp. 120–127, Dec. 2014.

[7] D. Kurrant and E. Fear, "Technique to decompose near-field reflection data generated from an object consisting of thin dielectric layers," *IEEE Trans. Antennas Propag.*, vol. 60, no. 8, pp. 3684–3692, Aug. 2012.

[8] J. S. Lee, C. Nguyen, and T. Scullion, "A novel, compact, low-cost, impulse ground-penetrating radar for nondestructive evaluation of pavements," *IEEE Trans. Instrum. Meas.*, vol. 53, no. 6, pp. 1502–1509, Dec. 2004.

[9] S. Wang, S. Zhao, and I. L. Al-Qadi, "Real-time density and thickness estimation of thin asphalt pavement overlay during compaction using ground penetrating radar data," *Surveys Geophys.*, vol. 41, pp. 431–445, 2020.

[10] J. Pan, C. Le Bastard, Y. Wang, and M. Sun, "Time-delay estimation using ground-penetrating radar with a support vector regression-based linear prediction method," *IEEE Trans. Geosci. Remote Sens.*, vol. 56, no. 5, pp. 2833–2840, May 2018.

[11] S. Wang, S. Zhao, and I. L. Al-Qadi, "Continuous real-time monitoring of flexible pavement layer density and thickness using ground penetrating radar," *NDT E Int.*, vol. 100, pp. 48–54, Dec. 2018.

[12] M. Sun, C. Le Bastard, Y. Wang, and N. Pinel, "Time-delay estimation using ESPRIT with extended improved spatial smoothing techniques for radar signals," *IEEE Geosci. Remote Sens. Lett.*, vol. 13, no. 1, pp. 73–77, Jan. 2016.

[13] C. Le Bastard, V. Baltazart, Y. Wang, and J. Saillard, "Thin-pavement thickness estimation using GPR with high-resolution and superresolution methods," *IEEE Trans. Geosci. Remote Sens.*, vol. 45, no. 8, pp. 2511–2519, Aug. 2007.

[14] I. L. AL-Qadi and S. Lahouar, "Measuring layer thicknesses with GPR-theory to practice," *Construct. Building Mater.*, vol. 19, no. 10, pp. 763–772, Dec. 2005.

[15] M. Sun et al., "Estimation of time delay and interface roughness by GPR using modified MUSIC," *Signal Process.*, vol. 132, pp. 272–283, Mar. 2017.

[16] M. Sun et al., "Time delay and interface roughness estimation using modified ESPRIT with interpolated spatial smoothing technique," *IEEE Trans. Geosci. Remote Sens.*, vol. 56, no. 3, pp. 1475–1484, Mar. 2018.
[17] M. Sun, N. Pinel, C. Le Bastard, V. Baltazart, A. Ihamouten, and Y. Wang, "Time delay and interface roughness estimations by GPR for pavement survey," *Near Surf. Geophys.*, vol. 13, no. 3, pp. 279-287, Jun. 2015.

[18] M. Sun, Y. Wang, N. Pinel, J. Li, and C. Le Bastard, "Road surface layers geometric parameters estimation by ground penetrating radar using estimation of signal parameters via rotational invariance techniques method," *IET Radar, Sonar Navigat.*, vol. 10, no. 3, pp. 603-609, Mar. 2016.

[19] B. Friedlander and A. J. Weiss, "Direction finding using spatial smoothing with interpolated arrays," *IEEE Trans. Aerosp. Electron. Syst.*, vol. 28, no. 2, pp. 574-587, Apr. 1992.

[20] A. J. Weiss and B. Friedlander, "Performance analysis of spatial smooching with interpolated arrays," *IEEE Trans. Signal Process.*, vol. 41, no. 5, pp. 1881-1892, May 1993.

[21] M. Guo, Y. D. Zhang, and T. Chen, "DOA estimation using compressed sparse array," *IEEE Trans. Signal Process.*, vol. 66, no. 15, pp. 4133-4146, Aug. 2018.

[22] J. Li, C. Le Bastard, Y. Wang, G. Wei, B. Ma, and M. Sun, "Enhanced GPR signal for layered media time-delay estimation in low-SNR scenario," *IEEE Geosci. Remote Sens. Lett.*, vol. 13, no. 3, pp. 299-303, Mar. 2016.

[23] J. A. Tropp and A. C. Gilbert, "Signal recovery from random measurements via orthogonal matching pursuit," *IEEE Trans. Inf. Theory*, vol. 53, no. 12, pp. 4655-4666, Dec. 2007.

[24] M. A. Davenport and M. B. Wakin, "Analysis of orthogonal matching pursuit using the restricted isometry property," *IEEE Trans. Inf. Theory*, vol. 56, no. 9, pp. 4395-4401, Sep. 2010.

[25] P. Stoica, *Introduction to Spectral Analysis*. Upper Saddle River, NJ, USA: Prentice-Hall, 1997.

[26] D. J. Daniels, *Ground Penetrating Radar*, 1st ed. Edison, NJ, USA: IET, 2004.

[27] P. Chargé, Y. Wang, and J. Saillard, "A non-circular sources direction finding method using polynomial rooting," *Signal Process.*, vol. 81, no. 8, pp. 1765-1770, Aug. 2001.

[28] C. Fauchard, "Utilisation de radars tres hautes frequences: Application a l'auscultation non destructive des chaussees," Ph.D. dissertation, Dept. Lab. Central Ponts Chaussées (LCPC), Univ. Nantes, Nantes, France, 2001.

[29] H. Akaike, "A new look at the statistical model identification," in *Selected Papers of Hirotugu Akaike*. Cham, Switzerland: Springer, 1974, pp. 215-222.

[30] M. Wax and T. Kailath, "Detection of signals by information theoretic criteria," *IEEE Trans. Acoust., Speech, Signal Process.*, vol. 33, no. 2, pp. 387-392, Apr. 1985.

[31] A. C. Gurbuz, O. Teke, and O. Arikan, "Sparse ground-penetrating radar imaging method for off-the-grid target problem," *J. Electron. Imag.*, vol. 22, no. 2, Jan. 2013, Art. no. 021007.

[32] P. Stoica, E. G. Larsson, and A. B. Gershman, "The stochastic CRB for array processing: A textbook derivation," *IEEE Signal Process. Lett.*, vol. 8, no. 5, pp. 148-150, May 2001.

[33] X. Dérobert et al., "Radar database collected over artificial debonding pavement structures during APT at the IFSTTAR's fatigue rollover," in *Proc. Accelerated Pavement Test. Transp. Infrastruct. Innov.*, in Lecture Note for Civil Engineering. Springer, 2020.

[34] A. M. Bruckstein, D. L. Donoho, and M. Elad, "From sparse solutions of systems of equations to sparse modeling of signals and images," *SIAM Rev.*, vol. 51, no. 1, pp. 34-81, Feb. 2009.

<figure><img src="image_13.png" /></figure>

**Meng Sun** (Member, IEEE) received the B.S. and M.S. degrees from Northwest University, Xi'an, China, and South China University of Technology, Guangzhou, China, in 2010 and 2013 respectively. He received, the Ph.D. degree from IETR Laboratory at Polytech Nantes, Université de Nantes, Nantes, France, in 2016.

He is with Electronic Engineering Department, Nanjing University of Aeronautics and Astronautics, Nanjing, China. His main research focus is signal processing techniques on parameter estimation of civil engineering materials.

<figure><img src="image_14.png" /></figure>

**Yide Wang** (Senior Member, IEEE) received the B.S. degree in electrical engineering from the Beijing University of Post and Telecommunication, Beijing, China, in 1985, and the M.S. and the Ph.D. degrees in signal processing and telecommunications from the University of Rennes, France, in 1986 and 1989, respectively.

He is a Professor with the Ecole Polytechnique de l'Université de Nantes. His research interests include array signal processing, spectral analysis, and mobile wireless communication systems.

<figure><img src="image_15.png" /></figure>

**Cédric Le Bastard** received the B.S. and M.S. degrees in electronic engineering from the University of Rennes, Rennes, France, in 2001 and 2003, respectively, and the Ph.D. degree and French Research Habilitation (D.Sc.) in signal processing from the Université de Nantes, Nantes, France, in 2007 and 2017, respectively.

<figure><img src="image_16.png" /></figure>

**Jingjing Pan** (Member, IEEE) received the B.S. and M.S. degrees in cartography and geographic information science from East China Normal University, Shanghai, China, in 2012, and Beijing Normal University, Beijing, China, in 2015, and the Ph.D. degree from Polytech Nantes, Université de Nantes, Nantes, France, in 2019.

She is a Post-Doctoral Fellow with the Electronic Engineering Department, Nanjing University of Aeronautics and Astronautics, Nanjing, China. Her research interests include array signal processing, ground penetrating radar, and nondestructive testing.

<figure><img src="image_17.png" /></figure>

**Vincent Baltazart** received the Ph.D. degree in signal processing from the University of Rennes, Rennes, France, in 1994.

From 1992 to 1993, he worked on ionospheric modeling and propagation at IPS Radio and Space Services, Sydney, Australia. From 1994 to 1996, he worked on microwave remote-sensing techniques at the Université Catholique de Louvain, Belgium. In 1996, he joined the Laboratoire Central des Ponts et Chaussées, France, as a Researcher in the field of optical sensors. He is involved in nondestructive testing and evaluation techniques for civil engineering applications with IFSTTAR, Nantes, France.