PUMA Applied to Time Delay Estimation for Processing GPR Data over Debonded Pavement Structures

: In this paper, principal-singular-vector utilization for modal analysis (PUMA) was adapted to perform time delay estimation on ground-penetrating radar (GPR) data by taking into account the shape of the transmitted GPR signal. The super-resolution capability of PUMA was used to separate overlapping backscattered echoes from a layered pavement structure with some embedded debondings. The well-known root-MUSIC algorithm was selected as a benchmark for performance assessment. The simulation results showed that the proposed PUMA performs very well, especially in the case where the sources are totally coherent, and it requires much less computational time than the root-MUSIC algorithm.

This paper focuses on a pavement survey using air-coupled GPR. Several defects that appear on a road surface such as cracks may be caused by debondings (delaminations) between the layers of the road pavement. Therefore, an early detection and characterization of debondings before any visual deterioration on the road surface can help maintenance services take appropriate measures at the right time, increasing the longevity of the road pavement and reducing the maintenance cost. In this paper, GPR aims at characterizing thin debonding layers (<1 cm thick) embedded between the first two top layers of the pavement structure [30][31][32]. The thickness of the debonding layers can be estimated with the help of the permittivity of the material found within the debonding interfaces and the time delay difference of the echoes reflected from the debonding interfaces [2,32,33]. The thin debonding thickness (<1 cm) is much smaller than the resolution capacity of the classical GPR algorithms, resulting in overlapping echoes for which the time delay difference ∆T is lower than the inverse of the GPR system bandwidth B. There are two ways by which the resolution of a GPR can be improved [3]: (I) by increasing the signal bandwidth of the GPR or (I I) by applying advanced signal processing methods to the recorded GPR data. The latter solution was preferred in this paper given that operators usually have little control over the bandwidth of conventional GPRs.
To this aim, two subspace-based algorithms were considered for high-resolution time delay estimation (TDE) in this paper, namely the conventional root-MUSIC algorithm [34] and the more recent PUMA algorithm [35]. Root-MUSIC was successfully used for the TDE of overlapping and nonoverlapping echoes using synthetic and field GPR data in [3], and it has been shown to be sensitive to the correlated multipath background. Specific preprocessing, namely sub-band averaging techniques, are then required to recover the super-time-resolution performance at the expense of additional computing time.
By contrast, PUMA is a more recent algorithm. It is an iterative subspace method based on linear prediction (LP) and weighted least squares (WLS) that was first introduced for spectral analysis in [36,37] and then for direction of arrival (DOA) estimation in [35,38,39]. PUMA has been shown to have a performance close to that of the maximum likelihood (ML) method, but with a lower computational cost.
In this paper, the proposed adaptation of PUMA to TDE is appealing for GPR application mainly because PUMA has been shown to be more efficient in the presence of coherent sources compared to other subspace algorithms [35,36]. In contrast to unitary-ESPRIT [40] and unitary root-MUSIC [41], PUMA can handle more than two coherent sources. Nevertheless, the proposed adaptation of PUMA to GPR data processing requires taking into account the shape of the GPR pulse in the mathematical formalism. PUMA is compared to the widespread root-MUSIC algorithm, for which new sub-band averaging techniques are introduced to handle correlated echoes.
The rest of the paper is organized as follows. The GPR data model taking the shape of the transmitted radar pulse into account is presented in Section 2. Section 3 briefly reviews the principle of root-MUSIC for TDE and details the adaptation of PUMA to process the data model in Section 2 for TDE. In Sections 4 and 5, the performance of both algorithms is compared in terms of the robustness to noise, time resolution, and average computational time. Conclusions and perspectives are given in Section 6.
Notation: Throughout this paper, vectors and matrices are denoted by lowercase and uppercase boldface letters, respectively. (·) * , [·] T , (·) H , (·) −1 , and (·) † represent the complex conjugate, transpose, conjugate-transpose, matrix inverse, and pseudo inverse operations, respectively. E[·] is the statistical mean. ⊗ is the Kronecker product. vec(·) is the vectorization operator. tr(·) is the trace operator.â represents the estimate of a. diag(·) denotes a diagonal matrix. I m is the (m, m) identity matrix, and J m is the (m, m) exchange matrix with ones on its antidiagonal and zeros elsewhere.
[a] m refers to the m th element of any vector a.

Data Model
Let us consider a road pavement probed by a GPR at vertical incidence under the following assumptions:

•
The road pavement is a stratified medium consisting of K smooth and homogeneous layers with low permittivity and negligible electrical conductivity; • The impinging EM waves on the road pavement are assumed to be plane waves emitted from an antenna in the far field; • Each layer contributes only one echo to the data model (single scattering assumption).
The backscattered signal recorded by the GPR under the above assumptions is the sum of the shifted and attenuated copies of the emitted radar pulse. The amplitude of each echo depends on the permittivity of the layers through Fresnel reflection coefficients. The received signal can be expressed as follows in the frequency domain [3,33]: where s k is the amplitude of the k th echo, which is considered to be independent of the frequency because the pavement layers have smooth interfaces, andẽ( f i ) andñ( f i ) are the Fourier transforms of the emitted and noise signals, respectively. Equation (1) can be written in matrix form as: with: .,ẽ( f N ) ): (N, N) diagonal matrix whose elements are the Fourier transforms of the emitted pulse; f i = f 1 + (i − 1) ∆ f : the frequency samples where f 1 is the beginning of the bandwidth, ∆ f is the frequency difference between two adjacent samples, and i = 1, 2, · · · , N.
The data covariance matrix Γ can be written as: where Γ s = E[ss H ] and σ 2 I N = E[nn H ] are the covariance matrices of the source vector and the noise vector, with dimensions (K, K) and (N, N), respectively.
Compared to the conventional model in the literature, the magnitude variations of the radar pulse over the bandwidth are taken into account in the data model and the covariance matrix in (1) and (3), respectively. For the GPR configuration used in this application, the backscattered echoes consist of multipath from the same transmitter (Tx). The echoes in the source vector s are then expected to be highly correlated with each other, and the associated covariance matrix Γ s may be rank deficient, i.e., 1 ≤ rank(Γ s ) ≤ K. Different correlation scenarios are thus considered in Section 4 for testing the TDE algorithms presented in the next section.

Time Delay Estimation Algorithms
This section presents two root-finding algorithms that are used in Section 4 to process GPR data for TDE. Firstly, the well-known root-MUSIC algorithm is recalled, and then, the proposed adaptation of PUMA to TDE is presented. Both algorithms take advantage of the eigendecomposition of the covariance matrix in different ways in order to estimate the time delays of the echoes.

Eigendecomposition
Assuming that Γ s is of full rank, without loss of generality (in the case where the sources are totally coherent, the sub-band averaging techniques can be used for recovering the rank of Γ s ), the covariance matrix of the received signal can be expressed as follows after eigendecomposition: where the columns of E s = [u 1 , · · · , u K ], which are the eigenvectors associated with the first K largest eigenvalues in D s = diag(λ 1 , · · · , λ K ), span the signal subspace and the columns of E n = [u K+1 , · · · , u N ], which are the eigenvectors associated with the N − K smallest eigenvalues in D n = diag(λ K+1 , · · · , λ N ), span the noise subspace. The eigenvalues are arranged in decreasing order: Using (3) and (4), From (5), where T is a nonsingular matrix of dimension (K, K). Equation (6) means that span(ΛA) = span(E s ) ⇒ span(A) = span(Λ −1 E s ).

Root-MUSIC
MUSIC is a well-known subspace method, which was adapted to TDE for processing GPR data in [3]. It is based on the orthogonality between the signal subspace, spanned by the columns of ΛA, and the noise subspace generated by the columns of E n = [u K+1 , · · · , u N ]. For the data model in (3), MUSIC searches for the K time delays that minimize the following function: The polynomial version of MUSIC, namely root-MUSIC, achieves better performance with a lower computational burden. Root-MUSIC searches for the complex roots z of the following 2N − 2-degree polynomial: The time delays are estimated from the phases of the K rootsẑ k of F RM (z) closest to and inside the unit circle as follows:

PUMA for TDE
In contrast to root-MUSIC, PUMA is an iterative subspace method [35,36]. According to the literature, PUMA may be directly applied to raw GPR data without preprocessing because it is more robust in the presence of correlated sources. In addition, it achieves a similar performance to that of conventional subspace methods in the presence of uncorrelated sources and has a lower average computational cost.
PUMA relies on linear prediction theory to estimate the time delays in the mode matrix A. A has a Vandermonde structure; thus, the following orthogonal relation holds [35] where: is an (N − K, N) Toeplitz matrix. The coefficients of B define a polynomial with roots lying on the unit circle as follows [35]: where c 0 = 0. Let c 0 = 1 and c = [c 1 , c 2 , · · · , c K ] T ; the time delays are estimated from the phases of the K roots of (12) with the use of (9). Equation (10) also implies that each m th element (m ≥ K + 1) in the k th column of A can be written as a linear combination of the previous K elements in that column as: In practice, we do not have access to A. However, following (6), the columns of matrix Λ −1 E s span the same subspace generated by the columns of the mode matrix A in the absence of noise. Let: thus, (13) can be rewritten as: Apart from the modification in (14), the rest of the algorithm is similar to that in [35]. For the completeness of presentation, the different steps of the algorithm are briefly recalled.
Equation (15) can be written in matrix form as: with: We can then write: and: where: Vector c can be found by solving: The least squares estimate of c can be obtained as: However, in the presence of noise, (21) gives a biased estimate of c, which eventually leads to inaccurate time delay estimates.
PUMA proposes to use weighted least squares to find c. It solves [35,36]: , and it is defined as: where T is a diagonal matrix defined as: (22) is given by [35,36,38]: It can be seen that (25) cannot be used directly because W depends on c; thus, an initial value of c needs to be found. The following steps summarize PUMA for TDE: 1. Initialize c with the least squares solution in (21): c = c LS ; 2. Calculate the weighting matrix W with c using (11) and (23) It was found that 3 iterations of Steps 2 and 3 were enough to satisfy the condition in Step 4 [36]. The time delays were estimated from (9) after finding the K roots of (12) with c obtained from the steps described above.

Simulated Data
In this section, the parameters used to simulate the GPR data over a simplified debonded pavement structure are presented. They serve to compute the data covariance matrices that are processed by the two selected TDE algorithms. In order to test the robustness of the selected algorithms within a multipath background, three different correlation scenarios were introduced, and two modified sub-band averaging techniques are proposed in Appendices A and B to mitigate the influence of coherent sources.

Parameters for the Pavement Survey
As shown in Figure 1, the defected pavement was modeled as a layered structure in which the debonding (Layer 2) was a thin layer with relative permittivity r2 = 10 embedded between Layer 1 and Layer 3, with relative permittivities r1 = 4 and r3 = 7, respectively. The pavement structure was probed with an air-coupled step-frequency GPR, of which the pulse was modeled as a Ricker pulse with processed bandwidth B = 3 GHz ([0.7 − 3.7] GHz) and f c = 2.2 GHz as the central frequency. N = 150 frequency samples were considered in the data vector. The noise was assumed to be white Gaussian. The thickness of Layer 1 was taken as d 1 = 7 cm, and that of Layer 3 was considered infinite. In the RMSE versus SNR simulations, the debonding region was set to d 2 = 0.5 cm; thus, the times of arrival associated with the backscattered echoes from the top surfaces of Layer 1, the debonding region, and Layer 3 were [τ 1 , τ 2 , τ 3 ] = [2.0, 2.9, 3.0] ns, respectively. The amplitudes of the backscattered echoes from Layer 1, the debonding, and Layer 3 were found to be s 1 = −0.333, s 2 = −0.200, and s 3 = 0.075, respectively, from the Fresnel equations. Figure 2 displays the shape of the Ricker pulse and the noiseless simulated signals over the debonded (d 2 = 0.5 cm) and the healthy areas.  To test the time resolution capability of the selected TDE algorithms in Section 5, the debonding thickness d 2 was varied from 0.05 cm to 3.97 cm, corresponding to a variation of B∆T 2 from 0.03 to 2.51, with ∆T 2 = τ 3 − τ 2 . Conventional FFT-based methods can easily resolve echoes when B∆T ≥ 1, but they fail to detect overlapping echoes, i.e., B∆T < 1.

Correlation Scenarios
The two selected algorithms were performed on the simulated data covariance matrix in (3). To test the robustness of the algorithms, the following three simulation scenarios were considered by varying the correlation coefficients in the source covariance matrix Γ s in (3), which is defined for three echoes as follows: where s k with k = {1, 2, 3} are the amplitudes of the echoes and ρ ij are the correlation coefficients between the i th and j th echoes (with i = {1, 2, 3}, j = {1, 2, 3}, i = j).

Scenario 1: Uncorrelated Echoes
This scenario investigated the best performance of the selected algorithms by assuming that the 3 backscattered echoes were fully uncorrelated. This correlation scenario was believed to match to the performance achieved by the Cramer-Rao bound. The full rank data covariance matrix corresponding to this case was simulated by setting ρ 12 = ρ 13 = ρ 23 = 0 in (26).

Scenario 2: Full Correlation between Overlapping Echoes
In this scenario, the last two echoes were correlated with each other while uncorrelated with the first one; thus, ρ 12 = ρ 13 = 0 and ρ 23 = 1 in (26). This scenario is similar to the first simulation carried out by Qian et al. in [35] within the scope of DOA application. It was expected to be closer to what can be observed in practice for the GPR application under scope, because the second and third echoes strongly overlapped.

Scenario 3: Full Correlation between Echoes
In this case, the three backscattered echoes were considered to be fully correlated; thus, ρ 12 = ρ 13 = ρ 23 = 1, reducing the rank of the covariance matrix of sources Γ s to 1. Compared to Scenarios 1 and 2, this scenario was expected to be the closest to the GPR experiment, where the backscattered echoes originated from the same transmitter.

Computed Data Covariance Matrix
The signal parameters in Section 4.1 were used to compute the data covariance matrix Γ in (3) that were processed by the two selected TDE algorithms. The noise covariance matrix is estimated from I independent snapshots as follows: 27) where N = [n 1 , · · · , n I ] is the noise data matrix.
To mitigate the influence of coherent sources in Scenarios 2 and 3, TDE was also performed on the modified FB (forward-backward) data covariance matrix (Γ FB Λ ) and the modified MSSP (modified spatial smoothing processing) data covariance matrix (Γ MSSP Λ ), which are detailed in Appendices A and B, respectively. Both matrices rely on the wellknown FB and MSSP sub-band averaging techniques [41,42]. The two proposed techniques avoid some limitations of the conventional FB and MSSP techniques when they are applied to the GPR data model in (2).
FB applies the averaging technique only once over the whole GPR bandwidth in both the forward and backward directions. For MSSP, the averaging technique was performed over overlapping frequency sub-bands of size L. According to [3], MSSP performs best when the effective bandwidth ratio L/N is between 80 % and 90 %. In this paper, L/N = 80 % was chosen.

Results
In this section, the performance of the proposed PUMA algorithm is compared to root-MUSIC for the three correlation scenarios introduced in the previous section.

Evaluation Criteria
The selected algorithms were compared in terms of the root-mean-squared error (RMSE) versus the signal-to-noise ratio (SNR), the RMSE versus B∆T 2 , and average computational time. The RMSE was focused on the two overlapping echoes to expect more sensitive results, and it is computed as follows: whereτ k,i is the k th time delay estimate obtained from the i th test and τ k is the true k th time delay. I = 100 is the number of Monte Carlo tests. For this application, the SNR is defined as the ratio between the sum of the squared magnitude of the two overlapping echoes to that of the noise power: Figures 3 and 4 show the RMSE versus the SNR and the RMSE versus B∆T 2 plots, respectively. It can be seen from Figure 3 that root-MUSIC is more robust to noise than PUMA. At SNRs ≥ 10 dB, the performance of PUMA approaches that of root-MUSIC. From Figure 4, it can be concluded that the selected algorithms can successfully resolve the two overlapping echoes from the debonding with root-MUSIC having better resolution capabilities than PUMA for SNR = 25 dB. The RMSE plots of root-MUSIC are used as the benchmark in the next scenarios.

Scenario 2: Full Correlation between the Overlapping Echoes
The two selected algorithms were successively applied on the data covariance matrix in (3) and on the FB data covariance matrix Γ FB Λ defined in (A6). During the processing of the former and the latter, (8) and (14) were computed using the pulse matrix Λ in (2) and Λ FB derived from (A4), respectively. The RMSE versus the SNR and the RMSE versus B∆T 2 plots are shown in Figures 5 and 6, respectively.
As expected, root-MUSIC failed to work on the covariance data matrix in (3), because it had a reduced rank equal to two. By contrast, PUMA succeeded in TDE for SNRs ≥ 20 dB, as shown in Figure 5, and likely achieved the best performance at a very high SNR.
The FB-root-MUSIC and FB-PUMA approach had the best performance, with FBroot-MUSIC performing better than FB-PUMA. The former and the latter outperformed root-MUSIC and PUMA because forward-backward averaging restored the rank of Γ s . It is worth noticing that under this scenario, the performance of PUMA in [35] corresponds to that of FB-PUMA in our paper.
The time resolution of both algorithms is displayed in Figure 6. Figure 6 shows that only PUMA can separate the overlapping echoes, i.e., for B∆T 2 < 1, when the algorithms are performed on the data covariance matrix. The time resolution of PUMA and root-MUSIC improved once they were applied to Γ FB Λ . At B∆T 2 = 0.3 (i.e., 0.5 cm), FB-PUMA and FB-root-MUSIC showed similar performance, and they both outperformed PUMA and root-MUSIC. However, the ultimate time resolution stayed lower than the one achieved with fully uncorrelated sources.

Scenario 3: Full Correlation between Echoes
As for Scenario 2, PUMA and root-MUSIC were successively applied on the data covariance matrix in (3), the FB data covariance matrix Γ FB Λ in (A6), and finally, the MSSP data covariance matrix Γ MSSP Λ in (A19). For MSSP, (8) and (14) were computed using Λ MSSP derived from (A17).  Figure 7, FB-root-MUSIC failed to work because the FB averaging technique can handle at most two correlated sources. FB-PUMA worked for SNRs > 6 dB, thus improving the performance of PUMA, which is limited to SNRs > 22 dB. MSSP-PUMA and MSSP-root-MUSIC had a similar performance, and they both outperformed FB-PUMA, FB-root-MUSIC, and PUMA because Γ MSSP Λ was of full rank. Surprisingly, MSSP-PUMA and MSSP-root-MUSIC performed better than the reference algorithm at SNRs < 10 dB, and MSSP was believed to improve the robustness to noise with its sub-band average.   Figure 8 shows that MSSP-PUMA and MSSP-root-MUSC outperformed FB-PUMA, PUMA, and FB-root-MUSIC when B∆T 2 ≤ 0.3 at an SNR of 25 dB. The "bumps" observed were due to the oscillations of the effective correlation coefficients obtained after the MSSP sub-band averaging technique. The magnitude of the "bumps" was smaller for PUMA owing to a better robustness against correlation.
Finally, in Figure 9, the average computational time of PUMA and root-MUSIC for SNR = 25 dB is plotted against the size of the data vector N. N was varied from 50 to 300 samples, and the computational time at each number of samples was averaged over 100 runs. The computational time increased with respect to N for both algorithms, with PUMA being about 10-times faster than root-MUSIC. All the simulations were carried out on a PC with Intel Core i5, 4 GB RAM, and 2.5 GHz.

Conclusions
In this paper, PUMA was successfully adapted to the processing of GPR data for time delay estimation by directly taking into account the radar pulse form. Two important modified sub-band averaging techniques were introduced. The simulation results illustrated the use of PUMA for the detection of thin debondings within pavement structures from air-coupled GPR data. The performance of the proposed method was compared with the classical root-MUSIC method. In terms of robustness against noise and time resolution capability, PUMA achieved a similar performance to that of the reference in the ideal scenario, i.e., fully uncorrelated echoes, and lower performance compared to the reference in more realistic scenarios with correlated echoes. Besides, PUMA was found to be about 10-times faster than root-MUSIC. In any case, the achieved time resolution clearly proves that PUMA can meet the needs of the application under scope by detecting subcentimeter-thick debondings within a pavement structure. However, the simulation results suggested that the PUMA algorithm can be further enhanced in the case of fully correlated echoes, and this can be done by extending the PUMA algorithm to directly exploit the properties of the rank-deficient source covariance matrix. In the near future, it is planned to test PUMA on the existing field database in [43].

Conflicts of Interest:
The authors declare no conflict of interest.
are matrices with dimensions (K, K) and (N, N), respectively. Γ FB s is the FB covariance matrix of the source vector [41], and Σ FB is the noise metric obtained after FB averaging. Let: where Λ FB is the modified pulse matrix with dimension (N, N). Thus, the noise metric, Σ FB , can be rewritten as: As opposed to the general assumption in the literature, the FB noise covariance matrix in (A2) corresponds to a colored noise, which limits the performance of the TDE algorithms.
To recover the ideal noise metric and then achieve optimal TDE performance, the following modified FB data covariance matrix is proposed: It is easily seen that the final formalism of Γ FB Λ in (A6) is similar to the initial formalism of the data covariance matrix in (3).

(K,K)
, MSSP is performed on the M overlapping sub-bands as follows: are matrices with dimensions (K, K) and (L, L), respectively. Γ MSSP s is the MSSP data covariance matrix of the source vector [42], and Σ MSSP is the noise metric obtained after applying MSSP. Let: where Λ MSSP is the modified pulse matrix with dimension (L, L). Therefore, the noise metric, Σ MSSP , can be written as: Hence, (A14) can be rewritten as follows: Given that the noise covariance matrix of Γ MSSP w is colored, in order to recover the ideal noise metric and achieve the optimal performance of the TDE algorithms, the following modified MSSP data covariance matrix, which whitens the noise and takes the shape of the radar pulse into account, is proposed: It can be observed that the formalism of the modified MSSP covariance matrix, Γ MSSP Λ , is similar to that of the data covariance matrix, Γ, in (3).