Time-Delay Estimation by a Modified Orthogonal Matching Pursuit Method for Rough Pavement

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.

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 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ 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τ ≤ 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 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 kth backscattered echo, which can be divided into two part: s k ( f ) = s k φ k ( f ); 4) s k is the reflection coefficient of the kth backscattered echo which is independent of frequency f ; 5) φ k ( f ) represents the frequency behavior of the kth backscattered echo; 6) n( f ) is an additive white Gaussian noise with zero mean and variance σ 2 . For N discrete frequencies within bandwidth B, (1) can be written as with the following notational definitions: received signal vector representing the measurements by a step-frequency radar; . . N ( f 1 is the lowest frequency and f is the frequency step); superscript T denotes the transpose operation.
diagonal matrix, whose diagonal elements represent the frequency behavior of the kth backscattered echo; vector with zero mean and covariance matrix σ 2 I, I is the (N × N) identity matrix. In the following, the data whitening procedure is applied. Divided by the radar pulse, the whitened received signal vector r w can be written as where b is the new noise vector after the data whitening procedure.

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: where matrix J is the (N × N) anti-identity matrix (exchange matrix) and the operator * denotes the complex conjugate. Therefore, the new mode vector of this extended radar data model nc a nc (t) can be expressed as 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) 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 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 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 kth iteration, the correlation between the residual signal r k−1 nc (the initial residual signal is set to be r 0 nc = r nc ) and the i th atom of the updated overcomplete dictionary matrix B(τ i ) can be defined as the absolute value of the sum of the elements of B H (τ i )r k−1 nc : where e = [1, . . . , 1] is a (1 × 2N) vector. The atom with the highest correlation (max i μ k (i )) is the best-fitting atom. Besides, the elements of ψ k (k = 1, 2, . . . , K ) corresponding to the amplitudes of the backscattered echoes (varying with frequency) are real. When τ i = t k , the elements of B H (τ i )r k−1 nc representing the estimated amplitude of the kth 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 B H (τ i )r k−1 nc is considered in the comparison. We can rewrite (7) as Obviously, (7) returns the absolute value of the sum of the real part of vector B H (τ i )r k−1 nc 's components. The time of arrival 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 kth 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 kth echo using LS optimization as follows: where operator + is the Moore-Penrose transpose, a + . Then, the residual signal at the (k + 1)th iteration can be written as 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 (ε r2 ) is about 4.5; the second layer is a base layer with relative permittivity ε 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 A and 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 σ h , correlation length L c . In the first simulation, three different rough pavements are studied as follows.   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τ = 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.

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τ = 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τ = 0.84); and the third echo is nonoverlapped with other echoes.     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 witht 1 = 0.998 ns,t 2 = 1.293 ns, and t 3 = 1.730 ns. It can be concluded that the proposed method can be applied to both overlapped and nonoverlapped echoes.

B. Performance Versus Bτ
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τ then changes from 0.3 to 3. The performance of the proposed method versus the product Bτ is evaluated by a Monte-Carlo process of 100 independent runs. The relative root-meansquare error (RRMSE) of the studied parameter is defined as wherex m represents the estimated parameter for the mth run of the method, and x the true value. In the simulation, parameter x can represent either the time-delay τ (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 ( τ ) by the proposed method and conventional OMP. As expected, it can be seen that the RRMSE on τ decreases with increasing of Bτ for both the proposed method and conventional OMP. As shown in Fig. 7, the conventional OMP method has limitation when Bτ ≤ 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  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 arrivalt 1 andt 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 ont 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 A overlying a very thin sand layer B with thickness less than 1 cm and C is the base-layer. The permittivity of the medium 3 is about 5. The length of 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 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 τ of the second and third echoes. In addition, the layer thickness can be calculated by the following equation: where H 2 is the estimated thickness of the sand layer, c is the speed of light in vacuum, and ε r3 is the permittivity of medium 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.  N 0 N 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 r w The number of backscattered echoes K OUTPUT: Estimated time of arrivalt k , k = 1, 2, . . . , K INITIALIZATION: Initial residual signal r 0 w = r w Initial iteration step k = 1 Initial index set 0 = ∅ while not converged (k < K ) do MATCH: μ k (i ) = a H (τ i )r k−1 w , i = 1, 2, . . . , N 0 Find the indext k that solves the optimization problem: t k = arg max i {μ k (i )} UPDATE: Amplitude of the kth echoŝ k = a + (t k )r k−1 w Residual signal r k w = r k−1 w − a(t k )ŝ k Iteration step k = k + 1 Index set k = k−1 ∪t k end while