InTechOpen uses cookies to offer you the best online experience. By continuing to use our site, you agree to our Privacy Policy.

Materials Science » "Non-Destructive Testing", book edited by Fausto Pedro Garcia Marquez, Mayorkinos Papaelias and Noor Zaman, ISBN 978-953-51-2502-0, Print ISBN 978-953-51-2501-3, Published: July 27, 2016 under CC BY 3.0 license. © The Author(s).

Chapter 9

Microwave Non‐Destructive Testing of Non‐Dispersive and Dispersive Media Using High‐Resolution Methods

By Cédric Le Bastard, Khaled Chahine, Yide Wang, Vincent Baltazart, Nicolas Pinel, Christophe Bourlier and Xavier Derobert
DOI: 10.5772/62410

Article top


Horizontally stratified medium with K layers probed by Ground Penetrating Radar and the associated radar profile (A‐scan) with three primary echoes.
Figure 1. Horizontally stratified medium with K layers probed by Ground Penetrating Radar and the associated radar profile (A‐scan) with three primary echoes.
Example of B‐scan of GPR.
Figure 2. Example of B‐scan of GPR.
RRMSE variation of the estimated thickness vs. SNR after 100 Monte Carlo simulations with two uncorrelated echoes; BΔτ = 0.5.
Figure 3. RRMSE variation of the estimated thickness vs. SNR after 100 Monte Carlo simulations with two uncorrelated echoes; BΔτ = 0.5.
RRMSE variation of the estimated thickness vs. BΔτ product after 100 Monte Carlo simulations with two uncorrelated echoes; SNR=15dB
Figure 4. RRMSE variation of the estimated thickness vs. BΔτ product after 100 Monte Carlo simulations with two uncorrelated echoes; SNR=15dB
Backscattered raw data from the metallic plate at different heights, h1, h2, h3, and h4.
Figure 5. Backscattered raw data from the metallic plate at different heights, h1, h2, h3, and h4.
MSE of n1 and n2 versus bandwidth.
Figure 6. MSE of n1 and n2 versus bandwidth.
MSE of τ1 and τ2 versus bandwidth.
Figure 7. MSE of τ1 and τ2 versus bandwidth.
RAP‐MPM MSEs of n1 and n2 versus bandwidth.
Figure 8. RAP‐MPM MSEs of n1 and n2 versus bandwidth.
RAP‐MPM MSEs of τ1 and τ2 versus bandwidth.
Figure 9. RAP‐MPM MSEs of τ1 and τ2 versus bandwidth.
MSE of n1 and n2 versus SNR.
Figure 10. MSE of n1 and n2 versus SNR.
MSE of τ1 and τ2 versus SNR.
Figure 11. MSE of τ1 and τ2 versus SNR.
RAP‐MPM MSEs of n1 and n2 versus SNR.
Figure 12. RAP‐MPM MSEs of n1 and n2 versus SNR.
RAP‐MPM MSEs of τ1 and τ2 versus SNR.
Figure 13. RAP‐MPM MSEs of τ1 and τ2 versus SNR.
The RRMSE versus Bτ2.
Figure 14. The RRMSE versus Bτ2.

Microwave Non‐Destructive Testing of Non‐Dispersive and Dispersive Media Using High‐Resolution Methods

Cédric Le Bastard1, Khaled Chahine2, Yide Wang3, Vincent Baltazart4, Nicolas Pinel5, Christophe Bourlier3 and Xavier Derobert4
Show details


This chapter discusses the principle and application of two model‐based algorithms for processing non‐dispersive and dispersive ground penetrating radar (GPR) data over layered medium under monostatic antenna configuration. Both algorithms have been selected for their super‐time resolution capability and reduced computational burden; they allow GPR to measure a layer thickness smaller than the fraction of the dominant wavelength. For non‐dispersive data, the ESPRIT algorithm is generalized to handle different kinds of data models encountered in experiments and in the literature. For dispersive data, the proposed adaptation of the MPM algorithm allows recovering the full‐time resolution and jointly estimating the time delays and quality factors of a layered medium with reduced bias. Both processing techniques are applied to probe‐layered roadways for NDT&E purposes.

Keywords: ground penetrating radar, high‐resolution signal processing methods, non‐dispersive media, dispersive media, civil engineering

1. Introduction and context

Ground penetrating radar (GPR) is a non‐destructive testing tool widely used in the fields of defense, agriculture, and mainly in geological applications and civil engineering [18]. It allows measuring different surveyed media parameters or detecting and localizing buried objects [915 19]. In this chapter, we focus on horizontally stratified media such as roadways, walls, and soils. The inner structure of the media can be deduced from GPR profiles by means of echo detection and amplitude estimation. Some processing approaches to obtain such geometrical information can use the full‐waveform inversion [6, 7] or migration processing [8]. In this chapter, we focus on estimation‐based methods of low computational burden, which allow estimating only one, two, or three parameters.

Echo detection provides time‐delay estimation (TDE) associated with each layer [15]. TDE is usually performed using conventional FFT‐based methods (inverse FFT or cross‐correlation methods), for which the time resolution depends on system GPR frequency bandwidth. To afford enhanced time resolution, e.g., for measuring thin thickness of layered medium, advanced model‐based signal processing techniques are required.

This chapter shows that subspace‐based methods are relevant for GPR applications for improving the performance of classical TDE techniques. These techniques use simplified assumptions, namely monostatic radar configuration at nadir probing layered non‐dispersive media. Under these hypotheses, the data are conventionally modeled as the sum of shifted and attenuated copies of the Tx pulse. When dealing with experiments, the latter simplified data model is restricted to either lossless or weakly lossy materials, e.g., pavement, and to small depth‐layered structures. Subspace‐based techniques then provide unbiased TDE and super‐time resolution capability as opposed to the conventional FFT‐based techniques. Among the subspace‐based techniques, the ESPRIT algorithm has been selected, thanks to its accuracy and its reduced computational burden.

For dispersive media, the above simplified data model does not hold because the wavelets undergo time spread and distortion with depth. Conventional subspace techniques then provide biased TDE and highly degraded time resolution. A few specific advanced processing techniques have been then developed and tested to extend the application of subspace techniques to dispersive media. Given a particular dispersion–absorption model, namely the Jonscher's parameterization of the medium permittivity, the MPM algorithm has been adapted to dispersive media to estimate unbiased TDE and recover the super‐time resolution performance.

Both subspace methods are based on the frequency domain data model. In this case, each echo is mathematically represented by a complex pole from which the time delay can be readily estimated. For non‐dispersive media, the pole is purely imaginary, i.e., e2jπfTk; for dispersive media, the complex pole includes both real and imaginary parts of which the ratio depends on the Q factor, which characterizes the dispersion of the medium.

Therefore, this chapter develops two model‐based high‐resolution methods to estimate the different parameters of both non‐dispersive and dispersive media, such as time delay, thickness, permittivity, and quality factor Q. For non‐dispersive media, the ESPRIT algorithm is extended to the signal model taking into account the radar pulse and any noise characteristics. Furthermore, six special cases of the ESPRIT method are introduced and discussed. The influence of interface roughness and medium heterogeneity inside the medium is also discussed. For dispersive media, the matrix pencil method is extended to two variants of signal model: depth invariant (DI; Q does not vary from one layer to another) and depth variant (DV; Q varies from one layer to another). The performance of the proposed methods in terms of time resolution, accuracy, and root mean square error is tested with simulated data. The obtained results show that the two proposed methods succeed in improving the time resolution in the case of non‐dispersive data and in jointly estimating the time delays and quality factors in the case of dispersive data.

2. Non-dispersive media

2.1. Data model

In this section, we present the forward model used to process radar data with high‐resolution methods. First, numerous physical considerations are discussed. Then, the forward model is described. The high‐resolution methods are applied on a non‐dispersive medium like the roadway [15]. The medium considered here is assumed to be a horizontally stratified medium with K layers, as shown in Figure 1.


Figure 1.

Horizontally stratified medium with K layers probed by Ground Penetrating Radar and the associated radar profile (A‐scan) with three primary echoes.

A monostatic radar antenna is positioned in the far‐field zone of the probed medium (antenna beyond the Fraunhoffer limit) at nadir. A vertical radar profile (function of time axis, called A‐scan) is then obtained, as shown in Figure 1. When a sequence of A scans is acquired along the scanning direction by collecting radar data at different locations, the radargram image (which is also called B‐scan) is obtained, as shown in Figure 2.


Figure 2.

Example of B‐scan of GPR.

Each layer is composed of aggregates, which are embedded within the hosting material. The size of the aggregates (1–2 cm in size) influences the scattering of waves depending on the dominant radar wavelength. When the aggregates are smaller than the radar wavelength, each layer can be considered as a homogeneous medium with smooth interfaces, which is characterized by its permittivity, conductivity, and thickness. However, when the aggregates are bigger or of the same order of size as the radar wavelength, each layer must then be considered as a heterogeneous medium comprising random rough interfaces. In such case, each layer is characterized not only by its effective permittivity, conductivity, and thickness but also by the two random rough interfaces delimiting each layer. These two rough interfaces are defined by their height probability density function and their surface height autocorrelation function. Note that, following Koudogbo et al. [17] and Pinel et al. [18], at nadir and within studied GPR band, surface scattering dominates in comparison with volume scattering. This is an important simplification that allows us to avoid modeling each aggregate individually, which would involve volume scattering. Instead, we can model these aggregates by considering an effective permittivity. Besides, as thoroughly analyzed by Pinel et al. [18], the roughness of the interfaces plays a non‐negligible role in the backscattering process only for large bands, or in other words, for high frequencies: typically, from about 3 GHz. To probe the roadway, the pulse GPRs are more often used, and the typical frequency is between 0.1 and 3 GHz. Nowadays, some GPRs of “stepped frequency radar” kind are also employed. They work at higher frequencies, but the main drawback is the measurement time, which remains important.

Then, as a first approximation, the roadway structure is considered as a stratified medium consisting of K homogeneous layers with small dielectric contrast, for which the interface roughness can be neglected at GPR frequencies (1–2 GHz). The aggregates, which are included within the material, are supposed to be small enough compared to the dominant GPR wavelength. For example, the GPR pulse at 1.5 GHz has a dominant wavelength of 20 cm in air and about 9 cm within the pavement material, whereas the maximum aggregate size is about 2 cm. Thus, each layer can therefore be considered as a homogeneous medium, which is characterized by an effective dielectric constant and thickness.

According to Adous [20], the dielectric constant of pavement layer is generally between 4 and 8. The dispersivity of the medium is usually neglected because the conductivity is low (between 10-3 and 10-2 S/m [21]). Roadway thicknesses typically range from 1 to 40 cm. For the pavements, the thickness varies from 1 to 10 cm. For ultrathin asphalt surfacing (UTAS) pavement [8], the thickness decreases to 1–3 cm. In the latter case, the backscattered echoes from adjacent interfaces overlap with each other, and it becomes very difficult even impossible to correctly estimate the thickness by conventional FFT‐based methods. We also assume that the GPR measurements are performed in far field (i.e., the antennas are located beyond the Fraunhoffer limit), so that the received signal model simplifies to a one‐dimensional model [13], and an analytical expression of the signal can be deduced from the solutions of Maxwell's equations. In standard radar configuration, i.e., monostatic radar in far field at nadir, the received signal is composed of d backscattered echoes, which are reflected from the medium interfaces.

Providing the latter assumptions on the medium and the antenna configuration, the echoes are simply the time‐shifted and attenuated copies of the transmitted pulse e(t). The echo amplitudes ck are given by the Fresnel coefficients; they are reinforced by a large dielectric contrast between the layers. There are two kinds of echoes: the so‐called “primary” echoes and the “multiple” echoes. The first ones (primary echoes) encounter only one reflection from the interface, whereas the multiple echoes result from several reflections. The multiple echoes show negligible amplitude if the layered medium shows small dielectric contrasts, as previously assumed. Thus, the received signal s(t) is the sum of d echoes e(tTk), whose amplitudes are proportional to the Fresnel coefficients of each medium interface. The received signal, also called signal model, can then be written as follows:


where e(t) is the transmitted signal; d is the number of backscattered echoes, which are mostly composed of K primary echoes; ck is the amplitude of the kth echo; Tk is the time delay of the kth echo; and b(t) is an additive noise, which is assumed to represent the measurement uncertainties.

The dispersivity of the medium is supposed negligible. In frequency domain, the echo amplitudes ck are therefore constant over the frequency bandwidth, so the received signal is a linear combination of complex exponentials that are modulated by the radar pulse e(f) plus an additive noise. The received signal in frequency domain can then be written as follows:


where the .˜ symbol represents the corresponding Fourier transform of the time function (Equation 1), and f stands for the frequency. Suppose that N frequency samples are chosen within the frequency bandwidth B of the GPR. Thus, the received signal can be written in matrix form as follows:


  • s˜=[s˜(f1)....s˜(fN)]T is the (N×1) data vector, which may represent either the Fourier transform of the pulse GPR signal or the measurements from a stepped frequency radar. The superscript T indicates the transpose operator;

  • Ω=diag(e˜(f1)...e˜(fN)) is the (N×N) diagonal matrix whose diagonal elements are the amplitudes of the Fourier transform of the radar pulse. The radar pulse is supposed known; it can be measured by the signal backscattered from a metallic plane;

  • A=[a(t1)...a(td)]T is the (N×d) mode matrix;

  • a(tk)=[e2jπf1Tk...e2jπfNTk]T is the (N×1) mode vector or steering vector;

  • c=[c1cd]T is the (d×1) vector of echoes amplitudes;

  • b˜=[b˜(f1)...b˜(fN)]T is the (N×1) zero mean noise vector;

  • fn=f1+(n1)Δf, where f1 is the beginning of the band, and Δf is the frequency difference between two adjacent samples.

The covariance matrix is defined as Γ=s˜s˜H, where the operator denotes the ensemble average, and the superscript H indicates the complex conjugate transpose. In practice, the covariance matrix Γ is estimated from M independent snapshots of the data vector s˜ as follows:


According to signal model (3) and assuming the noise to be independent of the echoes, the covariance matrix Γ of s˜ can be written as follows:


where Γsig and Σ0 are the covariance matrices of the signal and noise, respectively. The matrix Γc is the (d×d) dimensional covariance matrix of the vector of echo amplitudes. The noise covariance has been normalized such that tr(Σ0)=N[22] and tr(.) is the trace operators.

2.2. High‐resolution method ESPRIT applied to TDE

When enhanced time resolution is required for separating overlapping echoes (e.g., probing thin layered media), conventional FFT‐based methods are not efficient. Advanced signal processing methods, such as subspace‐based methods, have emerged in the 1980s to elegantly solve this problem. These methods are called subspace‐based methods because they exploit the properties of the eigenstructure of the covariance matrix Γ [22]. Among the existing methods, ESPRIT has been selected because it affords direct parameter estimation with a lower computational load.

The ESPRIT algorithm was firstly proposed for array signal processing by Roy et al. [23]. Several extensions of the ESPRIT algorithm have been proposed for TDE in the literature. N‐ESPRIT by Saarnisaari [24] can handle the colored noise but assumes an ideal radar pulse, whereas P‐ESPRIT by Swindlehurst [25] takes into account the radar pulse but supposes a temporally white Gaussian noise. Compared to the latter extensions, the M‐ESPRIT algorithm, which is presented in this section, is matched to the general data model in Equation (3) by taking both the radar pulse shape and any noise into account [16].

The principle of the ESPRIT method is to exploit some properties of rotation invariance [23]. The first step of the algorithm consists in dividing the data vector s˜ into two overlapping subvectors. Each subvector has N1 elements and overlaps with each other by N2 elements. Thus, the ((N1)×d) dimensional mode matrices of each sub‐band, A1 and A2, are related to each other by the (d×d) diagonal matrix Ψ, whose elements depend on the time‐delay parameter to be estimated as follows:

such that:




The matrix Ψ defined in Equation (8) cannot be estimated from data. On the basis of the generalized singular value decomposition (GSVD) of the covariance matrix Γ, the diagonal elements of Ψ can be retrieved from the similar matrix Φ defined thereafter in Equation (17c). From the GSVD of the covariance matrix Γ, we have the following equation:




where the generalized eigenvectors vk associated with the signal subspace are arranged in the matrix Vsig as columns; λi is the ith generalized eigenvalue of covariance matrix Γsig; and σ2 is the smallest generalized eigenvalue of the noise subspace. ESPRIT exploits the linear relation in Equation (6) within the GSVD of Γ. Thus, the GSVD of Equation (9) can be written in each data sub‐band as follows:


where Γsig,1 and Γsig,2 are the two (N1×N) dimensional matrices defined as follows:


where Ω1=diag(e˜(f1)...e˜(fN1)) and Ω2=diag(e˜(f2)...e˜(fN)) are the two ((N1)×(N1)) diagonal matrices. Σ0,1 and Σ0,2 are the two overlapping band matrices of size ((N×1)×N), and defined as the N1 upper lines and the N1 lower lines of the noise covariance matrix as follows:


Equations (6), (12), and (14) are used to obtain the following expression:


where T=ΓcAHΩHVsigΥ1 a (K×K) invertible matrix. Then, Equations. (11), (13), and (16) are used, and after some mathematical manipulations, we obtain the following equation:



Equation (17c) means that the matrices Φ and Ψ are similar, and thus they have the same eigenvalues, which are expressed as e2jπΔfTk. As opposed to Ψ,Φ can be retrieved from data. Thus, the d time delays of the echoes can be retrieved from the argument of the d eigenvalues of the matrixΦ In practice, the covariance matrix is estimated, and only the matrices Σ^0,1, Σ^0,2, and V^sig are estimated. The matrices Ω1 and Ω2 are obtained by the measurements of the radar pulse over a metallic plane. From Equation (17a), there exists two ways to calculate the matrix Φ. The first one uses the least square (LS) solution as follows:


The second way is the total least square (TLS) solution, which is known to improve the robustness of the solution especially at low signal‐to‐noise ratio (SNR) [26]. The TLS solution is computationally more demanding and requires performing the GSVD of the couple of matrices [Σ^0,2V^sig,ΘΣ^0,1V^sig].

2.3. M‐ESPRIT applied to specific data models

M‐ESPRIT is a general algorithm, which can handle different data configurations encountered in the literature and in GPR experiments. This section reviews the simplifications in Equation (17a), which corresponds to the six specific data models in Table 1. The assumptions on the data model imply the modifications of both the noise metric and the signal subspace. The corresponding equations to be solved for matrix Φ are also derived in Table 1, column 6. The LS and the TLS solutions for matrix Φ can then be readily established.

Case Assumptions Data model (s˜)
and Covariance matrix (Γ)
Equation (17c)
to solve for K
Dirac pulse and  
white noise
Radar pulse and
any noise
Radar pulse
and white
Dirac pulse
and any noise
Variant of
case 4
Whitened data
and any
With Σw=Ω1Σ0ΩH
Variant of
case 5
to deal with
specific pre‐
In case of

Table 1.

Expressions of 6 data models, the corresponding covariance matrices and signal subspaces, and the associated solutions to Equation (17a).

Case 1 refers to the conventional ESPRIT algorithm [23] with the ideal data model. In comparison with the data model in Equation (3), both the radar pulse matrix Ω and the covariance matrix of noise Σ0 reduce to the (N×N) identity matrix. In this case, Vsig,1 and Vsig,2 are the two overlapping sub‐band matrices of size ((N×1)×N), and they are defined as the N1 upper lines and the N1 lower lines of the matrix Vsig. Case 2 in Table 1 describes the general data model, which is supported by M‐ESPRIT presented in Section 2.2.

In comparison, Cases 3 and 4 describe intermediate solutions for data modeling. For Case 3, the radar pulse is taken into account in the data model by matrix Ω. The data model in Case 4 is composed of time‐shifted Dirac pulses with any noise (Σ0).

The two last data models, Cases 5 and 6, have a practical use when further preprocessing is performed. The whitening of the data by the radar pulse in Case 5 makes it possible to recover the situation of ideal Dirac pulses like for conventional ESPRIT, but in this case, the noise is modified. In this case, the Case 4 can be used.

In Table 1, the sixth data model provides a way to deal with the situation of totally correlated backscattered echoes. When the echoes are totally correlated, the full rank property of the covariance matrix of the sources Γc does not hold anymore, and the ESPRIT algorithm fails to work. This situation is likely to be encountered within a multipath environment because the source vector c in Equation (3) represents the different paths from the same original source. To resolve this problem, the covariance matrix may be estimated by using an averaging preprocessing technique [27]. This technique first whitens the data by the radar pulse and then organizes the frequency bandwidth (N samples) into M sub‐bands of L data each, overlapping with each another by L  1 samples. Consequently, the covariance matrix to be processed by the subspace algorithm is reduced to the size L × L. L will be referred to below as the effective frequency bandwidth; obviously, it must remain larger than the number of echoes d. In Table 1, matrix A1 is defined as the L first lines of the mode matrix in Equation (3). The matrices ΓcM and ΣM are the covariance matrices of sources and noise, respectively, after using a preprocessing method.

Furthermore, to decrease the computational burden of the ESPRIT method, linear versions have also been introduced by Zha [28] and Le Bastard et al. [29]. These methods use linear operations on the data and do not require the costly eigendecomposition of the covariance matrix. Further details are available in Refs. [2831].

2.4. Simulations, experiments, and discussion

2.4.1. Simulation

This section illustrates with simulated data the performance of the M‐ESPRIT method, in comparison with N‐ESPRIT, P‐ESPRIT, and the conventional ESPRIT. The performance of the algorithms has been evaluated with respect to both the time resolution capability and the relative root mean square error (RRMSE) on the estimated thickness. RRMSE on the estimated thickness is defined as follows:


where H is the true value of a parameter, and H^i is the corresponding estimated value in the ith simulation. I is the number of Monte Carlo trials. The simulations have been carried out for the thickness measurement of a two‐layered stratified medium. The time resolution describes the processing capability to separate two echoes within a limited band. It is thereafter defined as the BΔτ product, where B is the frequency bandwidth of the radar device, and Δτ is the smallest time shift between two echoes that the processing technique is able to distinguish. In this section, different data sets have been generated by varying the thickness from 1 to 41.2 mm, such that the corresponding BΔτ product spans the interval [0.027; 1], as shown in Figure 4.


Figure 3.

RRMSE variation of the estimated thickness vs. SNR after 100 Monte Carlo simulations with two uncorrelated echoes; BΔτ = 0.5.

The performance is established by a Monte Carlo process, which consists of I = 100 independent runs of algorithms. For each run, the covariance matrix was estimated from 100 independent snapshots. To analyze the influence of the noise, the simulations were carried out with a noise whose covariance matrix is non‐diagonal. The radar pulse e(t) was estimated in practice by the signal backscattered from a metallic plane. For the simulations, it is modeled as a Ricker pulse defined as the second derivative of a Gaussian pulse [31]. The parameters of the pulse have been chosen according to the characteristics of conventional GPR: central frequency of 1.5 GHz with frequency bandwidth of 2 GHz. The data set is composed of N = 41 equi‐spaced frequency samples within a 2 GHz bandwidth. According to Adous [20], the dielectric constants of the first and second layers were set to 4 and 7, respectively, and the conductivity of the two layers was neglected. The SNR was fixed with respect to the second echo and defined as the ratio between the power of the second echo and the noise variance.

Figure 3 shows the RRMSE variations of the estimated thickness versus SNR in the context of overlapping uncorrelated echoes. It can be seen that the RRMSE continuously decreases with increasing SNR for all algorithms (except for N‐ESPRIT). M‐ESPRIT allows obtaining the best performance for all SNRs. In contrast, the algorithms N‐ESPRIT and conventional ESPRIT provide a large and constant RRMSE meaning unrealistic biased TDE; both algorithms are not adapted to the application because they do not take into account the radar pulse or/and the noise metric. In this simulation, P‐ESPRIT performs better than the two latter algorithms but is expected to reach a constant RRMSE at higher SNR, meaning an asymptotically biased TDE. To further analyze the performance, N‐ESPRIT was used after data whitening by the radar pulse, referred to in Figures 3 and 4 as N‐ESPRIT with whitening. This modification of the original algorithm proposed by Saarnisaari [24] strongly improves the performance, which comes closer to that of M‐ESPRIT. The latter result means that it is more important to take into account the radar pulse in the mode vector than in the noise metric.


Figure 4.

RRMSE variation of the estimated thickness vs. BΔτ product after 100 Monte Carlo simulations with two uncorrelated echoes; SNR=15dB

Figure 4 shows the RRMSE variations of the estimated thickness versus BΔτ product at medium SNR for overlapping uncorrelated echoes. It can be seen that the RRMSE increases continuously with decreasing BΔτ product. Like in Figure 3, M‐ESPRIT shows the best performance followed by N‐ESPRIT with whitening. In contrast, the other algorithms (conventional ESPRIT, N‐ESPRIT, and P‐ESPRIT) show poorer performance.

2.4.2. Experiments

The performance of the algorithms was also tested in laboratory on smooth materials. A wideband exponential tapered slot antenna (ETSA) [32] is combined with a broadband vector network analyzer (VNA) to provide the experimental ultra‐wide band (UWB) stepped‐frequency radar (SFR). Experiment 1

First, the algorithms presented in the previous sections were tested on a virtual two‐layered medium filled in with air. This virtual medium was created with the backscattered data recorded by the SFR from the metallic plates at two different heights. The data from the virtual two‐layered medium is then considered as the sum of the two backscattered data as follows: s˜(f)=s˜hi(f)+s˜hj(f), where s˜hi(f) and s˜hj(f) are the backscattered data from the metallic plate at the heights hiand hj, respectively, as shown in Figure 5.

First, the null correlation case is considered to provide the best performance of the algorithms. The full correlation case will be presented in the following subsection. Thus, the null correlation matrix Γnull was formed as the sum of the covariance matrices of each individual echo as follows: Γnull=s˜his˜hiH+s˜hjs˜hjH. The following four heights are considered to form three two‐layered virtual media: h1, h2= h1 +2.25 cm, h3 = h1+4.5 cm, h4 = h1+ 6.75 cm. They span three BΔτ products: BΔτ=0.3,BΔτ=0.6,andBΔτ=0.9, respectively. The frequency band (2–4) GHz is spanned with N= 201 samples. The measurements are carried out at high SNR (about 40 dB). The transient signal in the antennas is filtered with a time window, as shown in Figure 5 (grey curves).


Figure 5.

Backscattered raw data from the metallic plate at different heights, h1, h2, h3, and h4.

BΔτ = 0.3BΔτ = 0.6BΔτ = 0.9
N‐ESPRIT with whitening3.240.110.06

Table 2.

Relative bias (in %) on the thickness for 3 BΔτ products with conventional ESPRIT, M‐ESPRIT, and N‐ESPRIT with whitening.

Table 2 shows that the relative bias decreases with increasing BΔτ product. The conventional ESPRIT algorithm provides the largest error on the thickness estimation, followed by N‐ESPRIT and M‐ESPRIT. At high SNR, M‐ESPRIT performs better than N‐ESPRIT at small BΔτ product; both algorithms give similar results for BΔτ beyond 0.6. Experiment 2

Second, SFR comprising two ETSA antennas (whose transmitter and receiver are very close to each other, about 20 cm) is tested on a smooth PVC slab. The PVC slab has a dielectric constant of εr=2.97+0.015j and a thickness of 4 cm. The height of the antenna is about 70 cm (far‐field condition). The frequency band is (1-2.6) GHz with N = 81 samples. Thus, BΔτ product is 0.728. In this case, the echoes are totally correlated, thus the Case 6 of the Table 1 is applied. We used the MSSP preprocessing [8, 20] with a sub‐band size of about 90% of the whole frequency bandwidth (L=73). The bias on the estimated thickness is 3.11%, which highlights the good performance of the M‐ESPRIT algorithm.

3. Dispersive media

3.1. Constant‐Q data model for dispersive media

In Ref. [33], Bano modeled the frequency dependence of permittivity for radar wave propagation in constant‐Q media using Kjartansson's formulation for seismic wave propagation in viscoelastic media [34]. This formulation is based on a complex power law of frequency for the effective dielectric permittivity and results in a frequency‐independent quality factor Q of the form: 1Q=εeεe=tan[π2(1n)], where εe and εe are the real and imaginary parts of the effective dielectric permittivity, respectively, and the index n (0<n<1) is related to Q as follows: n=2πarctan(Q). The case n=1 corresponds to lossless propagation (Q=, i.e., no attenuation) presented in Section 2.

A backscattered complex signal, represented as a linear combination of d echoes, results from probing a horizontally stratified medium, under normal incidence of plane waves. An echo originates from the interface between two horizontally superposed layers. Each layer is assumed to be homogeneous and having a thickness e, a constant quality factor Q (or, equivalently, a dispersion index n), and a dielectric constant ε0. The stratified structure is said to be DI (Depth Invariant) if Q does not vary from one layer to another; otherwise, it is called a DV (Depth-variant) structure [35]. Considering Kjartansson's absorption–dispersion model and substituting for each layer of a DV structure the expression of the corresponding wavenumber k in the equation of a plane wave exp{-2jke} the following backscattered signal is obtained:


where cm=amexp{jθm} is the frequency‐independent complex amplitude of the mth primary echo, τl=2el/Vrl is the traveling time of the wave in the lth layer, Vrl is the reference phase velocity at an arbitrary reference angular frequency, and b(f) is a complex white Gaussian noise with zero mean and variance 2σ2. After sampling, the frequency variable f is replaced byfk=kfs, where fs is the sampling period or frequency shift. The backscattered signal becomes:

s(k)=m=1dcmΠl=1mzlk(1+nl)/2+b(k)    k=1,2,...,N

where the pole of the lth layer is given by the following equation:


The product in the signal model represents the cumulative effect of the (m1) traversed layers of the stratified medium on the mth echo.

The assumption of a DI structure allows rewriting Equation (21) in the following compact form:


where the pole of the mth primary echo is:




is the equivalent time delay of the mth echo. Note that, compared with Equation (21), the subscript of n disappears in Equation (23). From Equations (21) and (23), it can be seen that the DI case is a special case of the more general DV one. Under matrix form, the signal model of both structures is expressed by the following equation:

with the following notational definitions:

s=[s(1) s(2)  s(N)]T: N×1 data vectorA=[a1 a2  ad]:N×d matrix of mode vectorsam=[l=1mzl1(1+nl)/2 l=1mzl2(1+nl)/2  l=1mzlN(1+nl)/2]T:N×1 DV mode vectoram=[zm1(1+n)/2 zm2(1+n)/2  zmN(1+n)/2]T:N×1 DI mode vectorc=[c1 c2 cd]T:d×1 vector of complex amplitudesb=[b(1) b(2)  b(N)]T:N×1 noise vector.

Note that because of the DV structure, the product of the poles cannot be put in the form of zm as shown in Equation (24) and, consequently, mode vectors of different echoes have different forms.

3.2. Nonlinearity compensation through predistortion linearization

Using the first and second extension theorems in [36, Section 4] yields, the following geometric progression of data samples for a single echo is obtained [36]:

{s˜(k)}k=1N={cz,cz2(1+n)/2, cz2×2(1+n)/21, , cz(N1)×2(1+n)/2(N2)}

Evidently, any Hankel data matrix constructed from the above data sequence admits Vandermonde decomposition. Moreover, we note that the geometric progression arises from probing the medium with non‐uniform frequencies given by the following equation:

v(k)=[(k1)2(1+n)/2(k2)]2/(1+n)fs, k=1,2, , N.

Hence, v(k) can be considered as the transfer function or the nonlinearity of a predistorter. For a non‐dispersive medium (n=1), we point out that v(k) becomes equal to kfs. Probing a dispersive medium with v in lieu of fresults in a linear frequency dependency of the phase.

3.3. Depth‐invariant method

In this section, we develop interpolation‐based MPM for the problem of parameter estimation of a DI structure characterized by the data model in Equation (23). For the DI structure, the parameter estimation problem can now be stated as follows: Given the data sequence {s(k)}k=1N estimate the time delays {Tmeq}m=1d and the index n (or, equivalently, the quality factor Q) of the stratified medium.

The required data sequence is obtained by spline interpolation. This is achieved by first determining the spline interpolant from the uniform data sequence provided by the measurement device and then evaluating the spline interpolant at non‐uniform frequencies {v(k)}k=1N obtained for a given value of n. Assuming a priori knowledge of n and neglecting the interpolation errors, the discrete backscattered signal becomes as follows:


Then s˜=A˜c+b with the following definitions:


It is readily verified that the transformed mode vector a˜m admits a Vandermonde structure (its samples follow a geometric progression with as the common ratio). As a consequence, the matrices S˜1 and S˜2, obtained from the Hankel matrix of s˜ by deleting its last and first columns, respectively, admit the following decomposition:

S˜1=Z˜1CZ˜2   S˜2=Z˜1CZ˜0Z˜2


Z˜1=[z1z2zdz12(1+n)/2z22(1+n)/2zd2(1+n)/2z1(NL1)×2(1+n)/2(NL2)z2(NL1)×2(1+n)/2(NL2)zd(NL1)×2(1+n)/2(NL2)]Z˜2=[1z12(1+n)/21z1(2(1+n)/21)(L1)1z22(1+n)/21z2(2(1+n)/21)(L1)1zd2(1+n)/21zd(2(1+n)/21)(L1)]Z˜0=diag{z12(1+n)/21, z22(1+n)/21, zd2(1+n)/21}C=diag{c1, c2, cd}.

As in conventional MPM, each value of is a rank reducing number of the pencil S˜2λS˜1. The estimates of the index n and the time delays {Tmeq}m=1d are then:


These estimates allow determining the complex amplitudes using a LS fit whose solution is given by the following equation:


In practice, however, the value of n is unknown and needs to be estimated. To solve this problem, we propose a recursive approach summarized in the following steps [37]:

  1. Step 1. Use the data sequence s˜ provided by the spline interpolant at frequencies selected according to Equation (28) to construct a Hankel matrix S˜. Start with n=1 as initial value, which corresponds to the initial uniform data sequence provided by real measurements.

  2. Step 2. Run MPM and determine n and {Tmeq}m=1d by using Equations (31) and (32).

  3. Step 3. Compare the estimated value of n with the value used in Step 1. If they are different, i.e., Δn0, substitute the new estimate of n in Equation (28) to obtain a new set of non‐uniform frequencies at which the spline interpolant is to be evaluated. Repeat Steps 1–3 until convergence, i.e., Δn=0. This termination condition might never be satisfied for noisy data, and so it will be replaced by a predefined practical threshold, e.g., 106.

  4. Step 4. Use the estimates of n and {Tmeq}m=1d obtained in the last iteration to estimate the complex amplitudes using Equation (33).

3.4. Depth‐variant method

This section is concerned with the problem of parameter estimation of a DV structure characterized by the data model in Equation (21). The proposed algorithm can be seen as an extension or a generalization of the interpolation‐based MPM described in the above section to a DV data model. It first unifies the form of the mode vectors through an equivalent‐medium approximation. Indeed, the forward Q filtering effect experienced by each of the dechoes of a DV Q structure in a horizontally stratified earth can be approximated by the effect of a one‐layered medium characterized by its equivalent parameters  (Tmeq, nmeq)Tmeq is given by Equation (25), and  nmeq is given by the following equation [31]:


The equivalent‐medium approximation then allows writing:





The above approximation becomes exact for a DI structure, i.e., z^m=zm. Note also that upon setting {nl=1}l=1d, the above model reduces to the undamped exponential model. Second, it uses a projection scheme to strip the mode vectors alternately. The DV parameter estimation problem can be stated as follows: Given the data sequence {s(k)}k=1N estimate the time delays {τl}l=1d and the dispersion indices {nl}l=1d of the stratified medium.

The proposed algorithm, called recursive and alternately projected MPM (RAP‐MPM), is developed for the case of two echoes [38]. The extension for d>2 is straightforward; however, it requires more effort and notations and is not considered in this work. Suppose, for now, that n1eq and n2eq are known (n1eq is practically equal to n1) and consider the application of the nonlinear transformation in Equation (28) to the uniform frequencies. Setting n=n1eq in ν(k)  and probing the medium with {νk}1N give rise to the following discrete backscattered signal:


Nevertheless, from a practical point of view, as it is impossible to probe the medium with non‐uniform frequencies, the data sequence in Equation (38) is obtained, as for interpolation‐based MPM, by cubic spline interpolation. Neglecting the interpolation errors, the data model in matrix form becomes:

with the following notational definitions:


It is readily verified from the equation of a˜1 that the transformed mode vector of the first echo admits a Vandermonde structure, and so its data matrix S˜1 becomes of rank one. Consequently, the matrices S˜11 and S˜12 admit the following decomposition:




The corresponding matrix pencil can then be written as follows:


As in conventional MPM, the value of λ=z^12(1+n1eq)/21 reduces the rank of the pencil in Equation (41) to zero and ought to be a GEV of the total matrix pair containing both echoes. However, due to the DV structure of the dispersive medium, the Vandermonde structure can be restored for one mode vector at a time. In such a case, estimating λ as the GEV of the total matrix pair is valid provided that the contribution of other echoes is filtered out after each restoration. Therefore, λ will be first estimated coarsely as a GEV of the total matrix pair. The estimation will be refined in later steps of the algorithm. The estimates of n1eq and T1eq can then be deduced from λ as follows:

n1eq=14πarctan[(log z^12(1+n1eq)/21)(logz^12(1+n1eq)/21)]
T1eq=[(log z^12(1+n1eq)/21)2πfr(2(1+n1eq)/21)(fsfr)(1+n1eq)/2]

Note that we know a rough estimate of the pole, we can filter out the contribution of the associated mode vector to the Hankel data matrix S˜=S˜1+S˜2 through an orthogonal projection as in the work of Chen et al. [39]. However, contrary to Chen et al. [39] that used a QR decomposition of the known mode vector to determine the orthogonal subspace, we carry out a singular value decomposition of the Hankel data matrix constructed from a˜1 as S˜1=U1Σ1V1H. Using MATLAB notation, the orthogonal complement is then given by U1orth=U1(:  ,rank(Σ1)+1:end) The definition of U1orth reflects the fact that the rank of S˜1 or, equivalently, Σ1 would be overestimated, unless a˜1 admits a Vandermonde structure guaranteed by setting n=n1eq in Equation (28). More importantly, such a definition of U1orth enables, once n1eq and T1eq are estimated, to filter out the contribution of a˜1 even if the rank is overestimated or, equivalently, nn1eq. To summarize, the Vandermonde structure of the mode vector is essential only to parameter estimation as seen in Equations (42) and (43) but not to filtering as shown by the definition of the orthogonal complement. Left multiplying S˜ by the conjugate transpose of U1orth gives S˜pr:


As pointed out in Ref. [40], the multiplication from the left destroys the shift‐invariance property in the column space but maintains it in the row space. This means that S˜pr can still be used to estimate the parameters of the second echo. When the processed data sequence is noisy, the orthogonal projection is only an approximation due to the perturbation introduced by noise in the signal subspace, and thus only the principal left singular vectors are to be considered when determining the orthogonal complement. The number of these vectors must be chosen in a way that ensures the best orthogonal projection on the one hand and maximum filtering of the first mode vector on the other hand. Obviously, these requirements are contradictory; orthogonal projection improves with decreasing number of singular vectors because the largest singular values are less affected by noise, while the filtering effect improves with increasing number of vectors. The optimum number of principal singular vectors as to minimizing estimation errors is hence found heuristically to be two, i.e., U1orth=U1(:  ,2+1:end). Similar to Equation (38), setting n=n2eq in Equation (28) restores the Vandermonde structure, only this time for the second mode vector. The estimates of n2eq and T2eq are then determined from the GEV of the matrix pair [S˜2pr,S˜1pr] using Equations (42) and (43). The following steps of this alternately projected approach are to filter out the contribution of the second mode vector and re‐estimate the parameters of the first and so on until practical convergence is achieved. The estimated values of {nmeq}m=1d and {Tm}m=1d then allow determining the complex amplitudes by means of a LS fit having the following solution:


Note that the development of the algorithm has till now been based on the a priori knowledge of n1eq and  n2eq. In practice, however, the values of the indices are unknown and need to be estimated. To solve this problem, we propose a recursive scheme of the previously described approach. The RAP-MPM (Recursive and Alternately Projected Matrix Pencil Method) can now be summarized in the following steps for the case of d=2 [38]:

  • Step 1. Form a Hankel data matrix S˜ from a data sequence corresponding to frequencies given in Equation (28). If the model order M<d (d is known), left multiply S˜ by the conjugate transpose of  Umorth. Start with M=d  and n=1 as initial values, which correspond to the exact uniform data sequence provided by real measurements.

  • Step 2. Run MPM and determine {nmeq}m=1d and {Tmeq}m=1d by using Equations (42) and (43).

  • Step 3. Compare the estimated value of nmeq with the value used in step 1. If they are different, i.e., Δn0 substitute the new estimate of nmeq in Equation (28) to obtain a new set of non‐uniform frequencies at which the spline interpolant is to be evaluated. Repeat steps 1, 2 and 3 until convergence, i.e., Δn=0. This termination condition might never be satisfied for noisy data, and so it will be replaced by a predefined practical threshold, e.g., 106

  • Step 4. Retain the estimates of nmeq and Tmeq of the last iteration and use them each time Umorth is to be determined from the mode vector a˜m for a given value of n.

  • Step 5. Set the model order M to d1 and repeat steps 1 to 4 until the relative change of the estimates between two consecutive alternated projections is smaller than a predefined value. Two to three projections were found to be sufficient.

  • Step 6. Use the estimates of {nmeq}m=1d and {Tm}m=1d to determine the complex amplitudes using Equation (45).


Figure 6.

MSE of n1 and n2 versus bandwidth.


Figure 7.

MSE of τ1 and τ2 versus bandwidth.

3.5. Performance evaluation

This section first examines the performance of both algorithms for two simulated backscattered echoes of a DI structure under scenarios of varying bandwidth, SNR, andBτ2. Then, it examines the performance of RAP‐MPM in the case of a DV structure. We consider a reference frequency fr=1.5GHz and a frequency shiftfs=10MHz. The frequency‐independent complex amplitudes are  c1=0.5exp(j0.7) and c2=1c12. Unless otherwise defined, SNR represents the average SNR over the data interval. The convergence threshold for both algorithms is fixed at 106, and the number of alternate projections for RAP‐MPM is set to three. For all simulations, the estimates of n and τ are computed for 500 independent trials. The mean squared errors (MSEs) are compared with the Cramer‐Rao lower bounds (CRLBs).

3.5.1. MSE versus bandwidth

This simulation is concerned with the ability of the developed algorithms to handle wideband signals. The bandwidth is varied symmetrically aboutfr. Figures 6 and 7 show the MSEs of the parameters for two (n1=n2=0.83,τ1 = 3.2 ns, τ2 = 1.2 ns and SNR=30dB) backscattered echoes along with the corresponding CRLBs as a function of bandwidth. It can be seen that the MSEs of both algorithms are comparable and mirror the CRLB as the bandwidth varies. This is in contrast with the behavior of conventional algorithms where their MSEs deviate from the CRLB beyond a certain bandwidth. Hence, the obtained results prove that the developed algorithms do not suffer from the bandwidth limitation of conventional estimation methods. A similar conclusion can be drawn for the DV structure (n1=0.87,n2=0.83, τ1= 3.2 ns, τ2 = 1.2 ns and SNR=30dB), as Figures 8 and 9 show that the estimates of RAP‐MPM have MSEs that decrease as bandwidth increases.


Figure 8.

RAP‐MPM MSEs of n1 and n2 versus bandwidth.


Figure 9.

RAP‐MPM MSEs of τ1 and τ2 versus bandwidth.

3.5.2. MSE versus SNR

This simulation investigates the robustness of the algorithms with respect to noise. Figures 10 and 11 depict the MSEs of the parameters for two backscattered echoes (n1=n2=0.83,τ1 = 3.2 ns, τ2 = 1.2 ns, and BW=2GHz) along with the corresponding CRLBs as a function of SNR. It can be seen that both algorithms are unbiased and robust to noise, as their MSEs approach the CRLBs even for low SNRs.

For the DV structure ( n1=0.87, n2=0.83,τ1=3.2ns,τ2=1.2ns, and BW=2GHz), the MSEs of the indices and the time delays approach their CRLBs starting, respectively, from SNR=5and 10dB as shown in Figures 12 and 13.


Figure 10.

MSE of n1 and n2 versus SNR.


Figure 11.

MSE of τ1 and τ2 versus SNR.


Figure 12.

RAP‐MPM MSEs of n1 and n2 versus SNR.


Figure 13.

RAP‐MPM MSEs of τ1 and τ2 versus SNR.

3.6. RRMSE versus Bτ2

This simulation examines the resolution capability of the algorithms by varying τ2  at BW = 2 GHz so that 0.2Bτ22. Figure 14 shows the RRMSE of the estimates of τ2  as a function of Bτ2  along with the CRLB and a 5% error threshold for n1=n2=0.83,τ1 = 3.2 ns, and SNR=30dB. We note here that defining SNR with respect to both echoes leads to a drastic decrease in the SNR of the second echo as Bτ2  increases. SNR is hence fixed with respect to the second echo. The obtained results illustrate that the algorithms afford less than 5% error for  Bτ2<1, which indicates their high‐resolution capability.


Figure 14.

The RRMSE versus Bτ2.

4. Conclusion

This chapter deals with model‐based high‐resolution signal processing algorithms applied to non‐dispersive and dispersive GPR data. For non‐dispersive data, a modified version of the conventional ESPRIT algorithm, called M‐ESPRIT, is presented for estimating the time delays. M‐ESPRIT allows taking the radar pulse and any noise signature into account and improves the time resolution of the GPR for thin pavement thickness estimation with conventional GPR of 2 GHz bandwidth. Contrary to the other versions of ESPRIT, M‐ESPRIT enables processing raw data directly without any preprocessing. Simulation and experimental results showed that M‐ESPRIT provides better performance. For dispersive data, a wideband extension of the matrix pencil method to the joint parameter estimation of dispersive media obeying the constant‐Q model is presented. The first and second extension theorems were used to derive a data sequence that enables the singular extension of the Hankel data matrix and, as a consequence, restores its effective rank. The desired data sequence, however, was found to correspond to a set of non‐uniform frequencies which are impossible to acquire from real measurements. To solve this problem, two MPM‐based algorithms were then developed. IB‐MPM and RAP‐MPM are based on interpolating the data sequence to reach the desired Vandermonde structure of the mode matrix. To handle the DV structure, RAP‐MPM makes use of the equivalent medium approximation to unify the form of the mode vectors and then uses a projection scheme to strip them alternately. The algorithms were subjected to simulation scenarios of varying bandwidth, SNR, and Bτ2. The results confirmed the wideband attribute of the algorithms, their unbiasedness over the considered SNR range, robustness to noise, rapid convergence, and high‐resolution capability.


This work may contribute to the European Cooperation in Science and Technology (COST) Action TU1208 Civil Engineering Applications of Ground Penetrating Radar ( [Accessed: 2016‐01‐26]).


1 - In: Proceedings of the 15th IEEE International Conference on Ground Penetrating Radar; 30 June–4 July 2014; Brussels. New York: IEEE; 2014.
2 - Daniels D.J. Ground‐Penetrating Radar. 2nd ed. London, UK: Institute of Electrical Engineers; 2004. ISBN: 086341360.
3 - Jol H. Ground Penetrating Radar Theory and Applications. Amsterdam (The Netherlands) and Oxford (UK): Elsevier; 2008. ISBN: 978‐0‐444‐53348‐7.
4 - Persico R. Introduction to Ground Penetrating Radar: Inverse Scattering and Data Processing. New Jersey (USA): Wiley‐IEEE Press; 2014. 392 p. ISBN: 978‐1‐118‐30500‐3.
5 - Benedetto A., Pajewski L. Civil Engineering Applications of Ground Penetrating Radar. New Jersey (USA): Springer; 2015. 371 p. ISBN: 978‐3‐319‐04812‐3.
6 - Van Der Kruk J., Streich R., Green A. Properties of surface waveguides derived from separate and joint inversion of dispersive TE and TM GPR data. Geophysics. 2006; 71: K19–K29. DOI: 10.1190/1.2168011.
7 - Lambot S., Slob E.C., Minet J., Jadoon K.Z., Vanclooster M., Vereecken H. Full‐waveform modeling and inversion of ground penetrating radar data for non‐invasive characterization of the soil hydrogeophysical properties. In: Viscara Rossel R.A., McBratney A.B., and Minasny B., editors. Proximal Soil Sensing, Developments in Soil Science Series. New Jersey (USA): Springer; 2010. p. 299–313.
8 - Radzevicius S. Practical 3‐D migration and visualization for accurate imaging of complex geometries with GPR. Journal of Environmental & Engineering Geophysics. 2008; 13(2): 99–112. DOI: 10.2113/JEEG13.2.99.
9 - AL‐Qadi I.L., Lahouar S. Measuring layer thicknesses with GPR theory to practice. Construction and Building Materials. 2005; 19(10): 763–772. DOI: 10.1016/j.conbuildmat.2005.06.005.
10 - Benedetto A., Benedetto F., De Blasiis M.R., Giunta G. Reliability of signal processing technique for pavement damages detection and classification using ground penetrating radar. IEEE Sensors Journal. 2005; 5(3): 471–480. DOI: 10.1109/JSEN.2005.846176.
11 - Lahouar S., Al‐Qadi I.L. Automatic detection of multiple pavement layers from GPR data. NDT&E International. 2008; 41(2): 69–81. DOI: 10.1016/j.ndteint.2007.09.001.
12 - Spagnolini U., Rampa V. Multitarget detection/tracking for monostatic ground penetrating radar: Application to pavement profiling. IEEE Transactions on Geoscience and Remote Sensing. 1999; 37(1): 383–394. DOI: 10.1109/36.739074.
13 - Spagnolini U. Permittivity measurements of multilayered media with monostatic pulse radar. IEEE Transactions on Geoscience and Remote Sensing. 1997; 35(2): 454–463. DOI: 10.1109/36.563284.
14 - Saarenketo T., Scullion T. Road evaluation with ground penetrating radar. Journal of Applied Geophysics. 2000; 43(2–4): 119–138. DOI: 10.1016/S0926‐9851(99)00052‐X.
15 - Le Bastard C., Baltazart V., Wang Y., Saillard J. Thin pavement thickness estimation using GPR with high‐resolution and super resolution methods. IEEE Transactions on Geoscience and Remote Sensing. 2007; 45(8): 2511–2519. DOI: 10.1109/TGRS.2007.900982.
16 - Le Bastard C., Baltazart V., Wang Y. Modified ESPRIT (M‐ESPRIT) algorithm for time delay estimation in both any noise and any radar pulse context by a GPR radar. Signal Processing. 2010; 90: 173–179. DOI: 10.1016/j.sigpro.2009.06.007.
17 - Koudogbo F., Mametsa H., Combes P. Surface and volume scattering from natural and manmade rough surfaces in the process of setting up data base coefficients. In: IEEE, editor. International Geoscience and Remote Sensing Symposium; July 2003; Toulouse, France; 2003. p. 4211–4213. DOI: 10.1109/IGARSS.2003.1295466.
18 - Pinel N., Le Bastard C., Baltazart V., Bourlier C., Wang Y. Influence of layer roughness for road survey by ground penetrating radar at nadir: Theoretical study. IET Radar, Sonar and Navigation. 2011; 5(6): 650–656. DOI: 10.1049/iet‐rsn.2010.0197.
19 - Bourlier C., Le Bastard C., Baltazart V. Generalization of PILE method to the EM scattering from stratified subsurface with rough interlayers – Application to the detection of debondings within pavement structure. IEEE Transactions on Geoscience and Remote Sensing. 2015; 53(7): 4104–4115. DOI: 10.1109/TGRS.2015.2390677.
20 - Adous M. Caractérisation electromagnetique des materiaux traites du genie civil dans la bande de frequence 50MHz‐13GHz. Thesis. University of Nantes; 2006.
21 - Fauchard C. Utilisation de radars tres hautes frequences: application a l'auscultation non destructive des chaussées. Thesis. University of Nantes; 2001.
22 - Schmidt R.O. A signal subspace approach to multiple emitter location and spectral estimation. Thesis. Stanford University; 1981.
23 - Roy R., Paulraj A., Kaikath T. ESPRIT—A subspace rotation approach to estimation of parameters of cisoids in noise. IEEE Transactions on Acoustics, Speech and Signal Processing. 1986; 34(5): 1340–1342. DOI: 10.1109/TASSP.1986.1164935.
24 - Saarnisaari H. TLS‐ESPRIT in a time delay estimation. In: Proceedings of the IEEE 47th Vehicular Technology Conference; Phoenix, USA; May 1997. p. 1619–1623. DOI: 10.1109/VETEC.1997.605832.
25 - Swindlehurst A.L. Time delay and spatial signature estimation using known asynchronous signals. IEEE Transactions on Signal Processing. 1998; 46(2): 449–462. DOI: 10.1109/78.655429.
26 - Roy R., Kailath T. Total least squares ESPRIT. In: Proceedings of the Asilomar Conference on Circuits, Systems and Computers; Pacific Grove, CA, USA; 1987. p. 297–301.
27 - Sun M., Le Bastard C., Wang Y., Pinel N. Time‐delay estimation using ESPRIT with extended improved spatial smoothing techniques for radar signals. IEEE Geoscience and Remote Sensing Letters. 2016; 13(1): 73–77. DOI: 10.1109/LGRS.2015.2497378.
28 - Zha H. Fast algorithms for direction‐of‐arrival finding using large ESPRIT arrays. Signal Processing. 1996; 48: 111–121. DOI: 10.1016/0165‐1684(95)00128‐X.
29 - Le Bastard C., Wang Y., Baltazart V. Algorithmes haute résolution linéaires pour l'estimation des temps de retard dans un contexte de bruit quelconque. Traitement du Signal. 2008; 25(4): 253–263. Available from: [2016‐01‐19]
30 - Marcos S., Marsal A., Benidir M. The propagator method for source bearing estimation. Signal Processing. 1995; 42: 121–138. DOI: 10.1016/0165‐1684(94)00122‐G.
31 - Schneider J.B. Plane waves in FDTD simulations and a nearly perfect. IEEE Transaction Antenna and Propagation. 2004; 52(12): 3280–3287. DOI: 10.1109/TAP.2004.836403.
32 - Guillanton E., Dauvignac J.‐Y., Pichot C., Cashman J. A new design tapered slot antenna for ultra‐wide band application. Microwave and Optical Technology Letters. 1998; 19(4): 286–289. DOI: 10.1002/(SICI)1098‐2760(199811)19:4<286::AID‐MOP12>3.0.CO;2‐0.
33 - Bano M. Constant dielectric losses of ground‐penetrating radar waves. Geophysical Journal International. 1996; 124(1): 279–288. DOI: 10.1111/j.1365‐246X.1996.tb06370.x.
34 - Kjartansson E. Constant Q‐wave propagation and attenuation. Journal of Geophysical Research. 1979; 84(B9): 4737–4748. DOI: 10.1029/JB084iB09p04737.
35 - Wang Y. Inverse Q‐filter for seismic resolution enhancement. Geophysics 2006; 71(3): 51–60. DOI: 10.1190/1.2192912.
36 - Chahine K., Baltazart V., Wang Y. Subspace leakage suppression for joint parameter estimation of quality factors and time delays in dispersive media. Circuits, Systems, and Signal Processing. 2015; 1: 15. DOI: 10.1007/s00034‐015‐0180‐8.
37 - Chahine K., Baltazart V., Wang Y. Interpolation‐based matrix pencil method for parameter estimation of dispersive media in civil engineering. Signal Processing. 2010; 90(10): 2567–2580. DOI: 10.1016/j.sigpro.2010.03.003.
38 - Chahine K., Baltazart V., Wang Y. Parameter estimation of damped power‐law phase signals via a recursive and alternately projected matrix pencil method. IEEE Transactions on Antennas and Propagation. 2011; 59(4): 1207–1216. DOI: 10.1109/TAP.2011.2109359.
39 - Chen H., Van Huffel S., Vandewalle J. Improved methods for exponential parameter estimation in the presence of known poles and noise. IEEE Transactions on Signal Processing. 1997; 45(5): 1390–1393. DOI: 10.1109/78.575717.
40 - Chen H., Van Huffel S., Van Den Boom A., Van Den Bosch P. Subspace‐based parameter estimation of exponentially damped sinusoids using prior knowledge of frequency and phase. Signal Processing. 1997; 59: 129–136. DOI: 10.1016/S0165‐1684(97)00085‐6.