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

Abstract

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.

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 [1–8]. It allows measuring different surveyed media parameters or detecting and localizing buried objects [9–15 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., $e^{-2j\pi f{T}_k}$; 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\left(t\right)$. The echo amplitudes c_k 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\left(t\right)$ is the sum of d echoes $e(t-T_k)$, 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:

$s(t)={\displaystyle {\sum }_{k=1}^d{c}_{k}e(t-T_k)+b(t)}$E1

where $e\left(t\right)$ is the transmitted signal; $d$ is the number of backscattered echoes, which are mostly composed of K primary echoes; c_k is the amplitude of the k^th echo; $T_k$ is the time delay of the k^th echo; and $b\left(t\right)$ 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 c_k 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\left(f\right)$ plus an additive noise. The received signal in frequency domain can then be written as follows:

$\tilde{s}(f)=\tilde{e}(f){\displaystyle {\sum }_{k=1}^d{c}_k{e}^{-2j\pi f{T}_k}+\tilde{b}(f)}$E2

where the $\widetilde{.}$ 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:

$\tilde{s}=\Omega Ac+\tilde{b}$E3

where:

• $\tilde{s}\text{=[}\tilde{s}\text{(}{f}_{\text{1}})\mathrm{....}\tilde{s}(f_N){]}^T$ is the $(N\times 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;

• $\Omega =\text{diag}(\tilde{e}(f_1)\mathrm{...}\tilde{e}(f_N))$ is the ($N\times 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(t_1)\mathrm{...}a(t_d)]}^T$ is the $(N\times d)$ mode matrix;

• $a(t_k)={[e^{-2j\pi f_1{T}_k}\mathrm{...}{e}^{-2j\pi f_N{T}_k}]}^T$ is the $(N\times 1)$ mode vector or steering vector;

• $c={\left[\begin{array}{ccc}{c}_1& \dots & c_d\end{array}\right]}^T$ is the $(d\times 1)$ vector of echoes amplitudes;

• $\tilde{b}={[\tilde{b}(f_1)\mathrm{...}\tilde{b}(f_N)]}^T$ is the $(N\times \text{1)}$ zero mean noise vector;

• $f_n=f_1+(n-1)\Delta f,$ where $f_1$ is the beginning of the band, and $\Delta f$ is the frequency difference between two adjacent samples.

The covariance matrix is defined as $\Gamma =⟨\tilde{s}{\tilde{s}}^H⟩$, where the $⟨\cdot ⟩$ operator denotes the ensemble average, and the superscript H indicates the complex conjugate transpose. In practice, the covariance matrix $\Gamma$ is estimated from M independent snapshots of the data vector $\tilde{s}$$$as follows:

$\widehat{\Gamma }=\frac{1}{M}{\displaystyle {\sum }_{m=1}^M{\tilde{s}}_m}{\tilde{s}}_m^H$E4

According to signal model (3) and assuming the noise to be independent of the echoes, the covariance matrix $\Gamma$ of $\tilde{\text{s}}$ can be written as follows:

$\Gamma ={\Gamma }_{sig}+{\sigma }^2{\sum }_0=\Omega A{\Gamma }_c{A}^H{\Omega }^H+{\sigma }^2{\sum }_0$E5

where ${\Gamma }_{sig}$ and ${\displaystyle {\Sigma }_0}$ are the covariance matrices of the signal and noise, respectively. The matrix ${\Gamma }_c$ is the $(d\times d)$ dimensional covariance matrix of the vector of echo amplitudes. The noise covariance has been normalized such that $tr({\Sigma }_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 $\Gamma$ [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 $\tilde{\text{s}}$ into two overlapping subvectors. Each subvector has $N-1$ elements and overlaps with each other by $N-2$ elements. Thus, the $((N-1)\times d)$ dimensional mode matrices of each sub‐band, $A_1$ and $A_2$, are related to each other by the $(d\times d)$ diagonal matrix $\Psi \text{,}$ whose elements depend on the time‐delay parameter to be estimated as follows:

$A_2=A_1\Psi$E6

such that:

$A=\left(\begin{array}{l}{A}_1\\ -\end{array}\right)=\left(\begin{array}{l}-\\ A_2\end{array}\right)$E7

and

$\Psi \text{=diag}(e^{-2j\pi \Delta f{T}_1}\mathrm{...}{e}^{-2j\pi \Delta f{T}_d})$E8

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

${\Gamma }_{sig}{V}_{sig}={\Sigma }_0{V}_{sig}\Upsilon$E9

with

$\Upsilon =diag({\lambda }_1-{\sigma }^2,\mathrm{...},{\lambda }_d-{\sigma }^2)$E10

where the generalized eigenvectors $v_k$ associated with the signal subspace are arranged in the matrix $V_{sig}$ as columns; ${\lambda }_i$ is the i^th generalized eigenvalue of covariance matrix ${\Gamma }_{sig}$; and ${\sigma }^2$ is the smallest generalized eigenvalue of the noise subspace. ESPRIT exploits the linear relation in Equation (6) within the GSVD of $\Gamma$. Thus, the GSVD of Equation (9) can be written in each data sub‐band as follows:

${\Gamma }_{sig,1}{V}_{sig}={\Sigma }_{0,1}{V}_{sig}\Upsilon$E11

${\Gamma }_{sig,2}{V}_{sig}={\Sigma }_{0,2}{V}_{sig}\Upsilon$E12

where ${\Gamma }_{sig,1}$ and ${\Gamma }_{sig,2}$ are the two $(N-1\times N)$ dimensional matrices defined as follows:

${\Gamma }_{sig,1}={\Omega }_1{A}_1{\Gamma }_s{A}^H{\Omega }^H$E13

${\Gamma }_{sig,2}={\Omega }_2{A}_2{\Gamma }_s{A}^H{\Omega }^H$E14

where ${\Omega }_1=\text{diag}(\tilde{e}(f_1)\mathrm{...}\tilde{e}(f_{N-1}))$ and ${\Omega }_2=\text{diag}(\tilde{e}(f_2)\mathrm{...}\tilde{e}(f_N))$ are the two $((N-1)\times (N-1))$ diagonal matrices. ${\Sigma }_{0,1}$ and ${\Sigma }_{0,2}$ are the two overlapping band matrices of size $((N\times 1)\times N),$ and defined as the $N-1$ upper lines and the $N-1$ lower lines of the noise covariance matrix as follows:

${\Sigma }_0=\left(\begin{array}{c}{\Sigma }_{0,1}\\ -\end{array}\right)=\left(\begin{array}{c}-\\ {\Sigma }_{0,2}\end{array}\right)$E15

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

${\Sigma }_{0,2}{V}_{sig}={\Omega }_2{A}_1\Psi T$E16

where $\mathbf{T}={\mathbf{\Gamma }}_c{\mathbf{A}}^H{\mathbf{\Omega }}^H{\mathbf{V}}_{sig}{\mathbf{{\rm Y}}}^{-1}$ a (K×K) invertible matrix. Then, Equations. (11), (13), and (16) are used, and after some mathematical manipulations, we obtain the following equation:

${\Sigma }_{0,2}{V}_{sig}=\Theta {\Sigma }_{0,1}{V}_{sig}\Phi$E17a

with

$\Theta ={\Omega }_2{\Omega }_1{}^{‐1}$E17b

$\Phi =T^{-1}\Psi T$E17c

Equation (17c) means that the matrices $\Phi$ and $\Psi$ are similar, and thus they have the same eigenvalues, which are expressed as $e^{-2j\pi \Delta f{T}_k}$. As opposed to $\Psi ,\Phi$ 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 ${\widehat{\Sigma }}_{0,1},{\widehat{\Sigma }}_{\text{0,2}}{\text{, and }\widehat{\text{V}}}_{sig}$ are estimated. The matrices ${\Omega }_1$ and ${\Omega }_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:

$\widehat{\Phi }={\left[{\left({\widehat{\Sigma }}_{0,1}{\widehat{V}}_{sig}\right)}^H{\widehat{\Sigma }}_{0,1}{\widehat{V}}_{sig}\right]}^{-1}{\left({\widehat{\Sigma }}_{0,1}{\widehat{V}}_{sig}\right)}^H{\Theta }^{‐1}{\widehat{\Sigma }}_{0,2}{\widehat{V}}_{sig}$E18

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 $\left[{\widehat{\Sigma }}_{0,2}{\widehat{V}}_{sig},\Theta {\widehat{\Sigma }}_{0,1}{\widehat{V}}_{sig}\right]$.

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 $\Phi$ are also derived in Table 1, column 6. The LS and the TLS solutions for matrix $\Phi$ can then be readily established.

Case	Assumptions	Data model ($\tilde{s}$); and Covariance matrix (Γ)	Signal; subspace	Noise; metric	Equation (17c); to solve for K
1; Conventional; ESPRIT	Dirac pulse and; white noise	$\tilde{s}=Ac+\tilde{b}$; $\Gamma ={A\Gamma }_c{A}^H+{\sigma }^{2}I$	$A$	$I$	$V_{sig,2}=V_{sig,1}\Phi$; $V_{sig}=\left(\begin{array}{c}{V}_{sig,1}\\ -\end{array}\right)=\left(\begin{array}{c}-\\ V_{sig,2}\end{array}\right)$
2; M‐ESPRIT	Radar pulse and; any noise	$\tilde{s}=\Omega Ac+\tilde{b}$; $\Gamma ={\Omega A\Gamma }_c{A}^H{\Omega }^H+{\sigma }^2{\Sigma }_0$	$\Omega A$	${\Sigma }_0$	${\Sigma }_{0,2}{V}_{sig}=\Theta {\Sigma }_{0,1}{V}_{sig}\Phi$; ${\Sigma }_0=\left(\begin{array}{c}{\Sigma }_{0,1}\\ -\end{array}\right)=\left(\begin{array}{c}-\\ {\Sigma }_{0,2}\end{array}\right)$; $\Theta =\text{diag}\left(\begin{array}{ccc}\frac{\tilde{e}\left(f_2\right)}{\tilde{e}\left(f_1\right)}& \cdots & \frac{\tilde{e}\left(f_N\right)}{\tilde{e}\left(f_{N-1}\right)}\end{array}\right)$
3; P‐ESPRIT	Radar pulse; and white; noise	$\tilde{s}=\Omega Ac+\tilde{b}$; $\Gamma ={\Omega A\Gamma }_c{A}^H{\Omega }^H+{\sigma }^{2}I$	$\Omega A$	$I$	$V_{sig,2}=\Theta V_{sig,1}\Phi$
4; N‐ESPRIT	Dirac pulse; and any noise	$\tilde{s}=\Omega Ac+\tilde{b}$; $\Gamma ={A\Gamma }_c{A}^H+{\sigma }^2{\Sigma }_0$	$A$	${\Sigma }_0$	${\Sigma }_{0,2}{V}_{sig}={\Sigma }_{0,1}{V}_{sig}\Phi$
5; Variant of; case 4	Whitened data; and any; noise	$\tilde{s}=Ac+{\Omega }^{-1}\tilde{b}$; $\Gamma ={A\Gamma }_c{A}^H+{\sigma }^2{\Sigma }^w$; With${\Sigma }^w={\Omega }^{-1}{\Sigma }_0{\Omega }^{-H}$	$A$	${\Sigma }^w$	${\Sigma }_{0,2}^w{V}_{sig}={\Sigma }_{0,1}^w{V}_{sig}\Phi$; ${\Sigma }^w=\left(\begin{array}{c}{\Sigma }_{0,1}^w\\ -\end{array}\right)=\left(\begin{array}{c}-\\ {\Sigma }_{0,2}^w\end{array}\right)$
6; Variant of; case 5; to deal with; specific pre‐; processing	In case of; sub‐band; averaging; techniques; (correlated; echoes)	${\tilde{s}}_m^w=A_1{\Psi }^{m-1}c+{\Omega }_m^{-1}{\tilde{b}}_m$; ${\Gamma }_{\text{M}}=A_1{\Gamma }_c^M{A}_1{}^H+{\sigma }^2{\Sigma }^M$	$A_1$	${\Sigma }^M$	${\Sigma }_{0,2}^M{V}_{sig}={\Sigma }_{0,1}^M{V}_{sig}\Phi$; ${\Sigma }^M=\left(\begin{array}{c}{\Sigma }_{0,1}^M\\ -\end{array}\right)=\left(\begin{array}{c}-\\ {\Sigma }_{0,2}^M\end{array}\right)$

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 $\Omega$ and the covariance matrix of noise ${\Sigma }_0$ reduce to the $(N\times N)$ identity matrix. In this case, $V_{sig,1}$ and $V_{sig,2}$ are the two overlapping sub‐band matrices of size $\left((N\times 1)\times N\right),$ and they are defined as the $N-1$ upper lines and the $N-1$ lower lines of the matrix $V_{sig}.$ 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 $\Omega .$ The data model in Case 4 is composed of time‐shifted Dirac pulses with any noise $\left({\Sigma }_0\right).$

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 ${\Gamma }_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\times 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 $A_{\text{1}}$ is defined as the $L$ first lines of the mode matrix in Equation (3). The matrices ${\Gamma }_c^M$ and ${\Sigma }^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. [28–31].

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:

$RRMSE(\%)=100\times \frac{\sqrt{\frac{1}{I}{\displaystyle {\sum }_{i=1}^I({\widehat{H}}_i}-H{)}^2}}{H}$E19

where H is the true value of a parameter, and ${\widehat{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\Delta \tau$ product, where $B$ is the frequency bandwidth of the radar device, and $\Delta \tau$ 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\Delta \tau$ 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; $\text{BΔτ}$ = 0.5.

The performance is established by a Monte Carlo process, which consists of $I=\text{ 1}00$ 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\left(t\right)$ 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. $\text{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\Delta \tau$ product at medium SNR for overlapping uncorrelated echoes. It can be seen that the RRMSE increases continuously with decreasing $B\Delta \tau$ 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).

2.4.2.1. 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: $\tilde{s}(f)={\tilde{s}}_{hi}(f)+{\tilde{s}}_{hj}(f),$ where ${\tilde{s}}_{hi}(f)$ and ${\tilde{s}}_{hj}(f)$ are the backscattered data from the metallic plate at the heights $h_i$and $h_j$, 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: ${\Gamma }_{null}=⟨{\tilde{s}}_{hi}{\tilde{s}}_{hi}{}^H⟩+⟨{\tilde{s}}_{hj}{\tilde{s}}_{hj}{}^H⟩.$ The following four heights are considered to form three two‐layered virtual media: $h_1, h_2= h_{\text{1}}+2.\text{25 cm}, h_{\text{3}}= h_1+4.\text{5 cm}, h_{\text{4}}= h_1+\text{ 6}.\text{75 cm}.$ They span three $B\Delta \tau$ products: $B\Delta \tau =0.3,B\Delta \tau =0.6,\text{and}B\Delta \tau =0.9,$ respectively. The frequency band (2–4) GHz is spanned with $N=\text{ 2}0\text{1}$ 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, h₁, h₂, h₃, and h_4.

	BΔτ = 0.3	BΔτ = 0.6	BΔτ = 0.9
ESPRIT	76.5	16.41	2.46
M‐ESPRIT	1.92	0.34	0.04
N‐ESPRIT with whitening	3.24	0.11	0.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\Delta \tau$ 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\Delta \tau$ product; both algorithms give similar results for $B\Delta \tau$ beyond 0.6.

2.4.2.2. 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 ${\text{ε}}_{\text{r}}\text{=}\text{2}\text{.97+0}\text{.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=\text{ 81}$ samples. Thus, $B\Delta \tau$ 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 $\left(\text{L}=\text{73}\right).$ 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: $\frac{1}{Q}=\frac{{{\text{ε}}^{″}}_{\text{e}}}{{{\text{ε}}^{\prime }}_e}=\tan \left[\frac{\pi }{2}(1-n)\right]$, where ${{\text{ε}}^{\prime }}_e$ and ${{\text{ε}}^{″}}_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=\frac{2}{\pi }\arctan (Q).$ The case $n=1$ corresponds to lossless propagation ($Q=\infty$, 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 ${\text{ε}}^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:

$s(f)={\displaystyle {\sum }_{m=1}^d{c}_m}{\Pi }_{i=1}^m\exp \left\{-2\pi f_r{\tau }_l\left(j+\tan \left[\frac{\pi }{4}(1-n_l)\right]\right){(\frac{f}{f_r})}^{(1+n_l)/2}\right\}+b(f)$E20

where $c_m=a_m\exp \left\{j{\theta }_m\right\}$ is the frequency‐independent complex amplitude of the mth primary echo, ${\tau }_l=2e_l/V_{rl}$ is the traveling time of the wave in the l^th layer, $V_{rl}$ 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{\sigma }^2$. After sampling, the frequency variable $f$ is replaced by$f_k=k{f}_s,$ where $f_s$ is the sampling period or frequency shift. The backscattered signal becomes:

$s(k)={\displaystyle {\sum }_{m=1}^d{c}_m{\Pi }_{l=1}^m}{z}_l^{k^{(1+n_l)/2}}+b(k)k=1,2,\mathrm{...},N$E21

where the pole of the l^th layer is given by the following equation:

$z_l=\exp \left\{-2\pi f_r{\tau }_l\left(j+\tan \left[\frac{\pi }{4}(1-n_l)\right]\right){\left(\frac{f_s}{f_r}\right)}^{(1+n_l)/2}\right\}.$E22

The product in the signal model represents the cumulative effect of the $(m-1)$ traversed layers of the stratified medium on the m^th echo.

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

$s\left(k\right)={\displaystyle {\sum }_{m=1}^d{c}_m{\stackrel{⌣}{z}}_m^{k^{\left(1+n\right)/2}}}+b\left(k\right)$E23

where the pole of the m^th primary echo is:

${\stackrel{⌣}{z}}_m={\displaystyle {\prod }_{l=1}^m{z}_l}=\exp \left\{-2\pi f_r{T}_m^{eq}\left(j+\tan \left[\frac{\pi }{4}\left(1-n\right)\right]\right){\left(\frac{f_s}{f_r}\right)}^{\left(1+n\right)/2}\right\}$E24

with

$T_m^{eq}={\displaystyle {\sum }_{l=1}^m{\tau }_l}$E25

is the equivalent time delay of the m^th 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:

$s=Ac+b$E26

with the following notational definitions:

$\begin{array}{ll}s\hfill & ={\left[s\left(1\right)s\left(2\right)\dots s\left(N\right)\right]}^T:N\times 1\text{ data vector}\hfill \\ A\hfill & =\left[a_1{a}_{\text{2}}\dots a_d\right]:N\times d\text{ matrix of mode vectors}\hfill \\ a_m\hfill & ={\left[{\displaystyle \prod _{l=1}^m{z}_l^{1^{\left(1+n_l\right)/2}}}{\displaystyle \prod _{l=1}^m{z}_l^{2^{\left(1+n_l\right)/2}}}\cdots {\displaystyle \prod _{l=1}^m{z}_l^{N^{\left(1+n_l\right)/2}}}\right]}^T:N\times 1\text{ DV mode vector}\hfill \\ a_m\hfill & ={\left[{\stackrel{⌣}{z}}_m^{1^{\left(1+n\right)/2}}{\stackrel{⌣}{z}}_m^{2^{\left(1+n\right)/2}}\dots {\stackrel{⌣}{z}}_m^{N^{\left(1+n\right)/2}}\right]}^T:N\times 1\text{ DI mode vector}\hfill \\ c\hfill & ={\left[c_1{c}_{\text{2}}\dots c_d\right]}^T:d\times 1\text{ vector of complex amplitudes}\hfill \\ b\hfill & ={\left[b\left(1\right)b\left(2\right)\dots b\left(N\right)\right]}^T:N\times 1\text{ noise vector}\text{.}\hfill \end{array}$E101

Note that because of the DV structure, the product of the poles cannot be put in the form of ${\stackrel{⌣}{z}}_m$ 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]:

${\left\{\tilde{s}\left(k\right)\right\}}_{k=1}^N=\left\{c\stackrel{⌣}{z},c{\stackrel{⌣}{z}}^{2\left(1+n\right)/2},c{\stackrel{⌣}{z}}^{2\times 2^{\left(1+n\right)/2}-1},\dots ,c{\stackrel{⌣}{z}}^{\left(N-1\right)\times 2^{\left(1+n\right)/2}-\left(N-2\right)}\right\}$E27

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\left(k\right)={\left[\left(k-1\right)2^{\left(1+n\right)/2}-\left(k-2\right)\right]}^{2/\left(1+n\right)}{f}_s,k=1,2,\dots ,\text{ N}\text{.}$E28

Hence, $v\left(k\right)$ can be considered as the transfer function or the nonlinearity of a predistorter. For a non‐dispersive medium $\left(n=1\right)$, we point out that $v\left(k\right)$ becomes equal to $k{f}_s$. Probing a dispersive medium with $v$ in lieu of $f$results 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 ${\left\{s(k)\right\}}_{k=1}^N$ estimate the time delays ${\left\{T_m^{eq}\right\}}_{m=1}^d$ 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 ${\left\{v\left(k\right)\right\}}_{k=1}^N$ 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:

$\tilde{s}\left(k\right)={\displaystyle {\sum }_{m=1}^d{c}_m{\stackrel{⌣}{z}}_m^{\left(k+1\right)2^{\left(1+n\right)/2}-\left(k-2\right)}}+\tilde{b}\left(k\right).$E29

Then $\tilde{s}=\tilde{A}c+b$ with the following definitions:

$\begin{array}{ll}\tilde{s}\hfill & ={\left[\begin{array}{llll}\tilde{s}(1)\hfill & \tilde{s}(2)\hfill & \dots \hfill & \tilde{s}(N)\hfill \end{array}\right]}^T\hfill \\ \tilde{A}\hfill & =\left[\begin{array}{llll}{\tilde{a}}_1\hfill & {\tilde{a}}_2\hfill & \dots \hfill & {\tilde{a}}_d\hfill \end{array}\right]\hfill \\ {\tilde{a}}_m\hfill & ={\left[\begin{array}{lllll}{\stackrel{⌣}{z}}_m\hfill & {\stackrel{⌣}{z}}_m^{2^{(1+n)/2}}\hfill & {\stackrel{⌣}{z}}_m^{2\times 2^{(1+n)/2}-1}\hfill & \dots \hfill & {\stackrel{⌣}{z}}_m^{(N-1)\times 2^{(1+n)/2}-(N-2)}\hfill \end{array}\right]}^T\hfill \\ \tilde{b}\hfill & ={\left[\begin{array}{llll}\tilde{b}(1)\hfill & \tilde{b}(2)\hfill & \dots \hfill & \tilde{b}(N)\hfill \end{array}\right]}^T.\hfill \end{array}$E102

It is readily verified that the transformed mode vector ${\tilde{a}}_m$ admits a Vandermonde structure (its samples follow a geometric progression with as the common ratio). As a consequence, the matrices ${\tilde{S}}_1$ and ${\tilde{S}}_2$, obtained from the Hankel matrix of $\tilde{s}$ by deleting its last and first columns, respectively, admit the following decomposition:

${\tilde{S}}_1={\tilde{Z}}_{1}C{\tilde{Z}}_2{\tilde{S}}_2={\tilde{Z}}_{1}C{\tilde{Z}}_0{\tilde{Z}}_2$E30

where

$\begin{array}{c}{\tilde{Z}}_1=\left[\begin{array}{llll}{\stackrel{⌣}{z}}_1\hfill & {\stackrel{⌣}{z}}_2\hfill & \cdots \hfill & {\stackrel{⌣}{z}}_d\hfill \\ {\stackrel{⌣}{z}}_{{}_1}^{2^{\left(1+n\right)/2}}\hfill & {\stackrel{⌣}{z}}_2^{2^{\left(1+n\right)/2}}\hfill & \cdots \hfill & {\stackrel{⌣}{z}}_d^{2^{\left(1+n\right)/2}}\hfill \\ ⋮\hfill & ⋮\hfill & \ddots \hfill & ⋮\hfill \\ {\stackrel{⌣}{z}}_{{}_1}^{\left(N-L-1\right)\times 2^{\left(1+n\right)/2}-\left(N-L-2\right)}\hfill & {\stackrel{⌣}{z}}_{{}_2}^{\left(N-L-1\right)\times 2^{\left(1+n\right)/2}-\left(N-L-2\right)}\hfill & \dots \hfill & {\stackrel{⌣}{z}}_d^{\left(N-L-1\right)\times 2^{\left(1+n\right)/2}-\left(N-L-2\right)}\hfill \end{array}\right]\\ {\tilde{Z}}_2=\left[\begin{array}{llll}1\hfill & {\stackrel{⌣}{z}}_1^{2^{\left(1+n\right)/2}-1}\hfill & \cdots \hfill & {\stackrel{⌣}{z}}_1^{\left(2^{\left(1+n\right)/2}-1\right)\left(L-1\right)}\hfill \\ 1\hfill & {\stackrel{⌣}{z}}_2^{2^{\left(1+n\right)/2}-1}\hfill & \cdots \hfill & {\stackrel{⌣}{z}}_2^{\left(2^{\left(1+n\right)/2}-1\right)\left(L-1\right)}\hfill \\ ⋮\hfill & ⋮\hfill & \ddots \hfill & ⋮\hfill \\ 1\hfill & {\stackrel{⌣}{z}}_d^{2^{\left(1+n\right)/2}-1}\hfill & \cdots \hfill & {\stackrel{⌣}{z}}_d^{\left(2^{\left(1+n\right)/2}-1\right)\left(L-1\right)}\hfill \end{array}\right]\\ {\tilde{Z}}_0=\text{diag}\left\{{\stackrel{⌣}{z}}_1^{2^{\left(1+n\right)/2}-1},{\stackrel{⌣}{z}}_2^{2^{\left(1+n\right)/2}-1},\dots {\text{, }\stackrel{⌣}{\text{z}}}_d^{2^{\left(1+n\right)/2}-1}\right\}\\ C=\text{diag}\left\{c_1,c_2,\dots \text{, }{c}_d\right\}.\end{array}$E103

As in conventional MPM, each value of is a rank reducing number of the pencil ${\tilde{S}}_2-\lambda {\tilde{S}}_1$. The estimates of the index n and the time delays ${\left\{T_m^{eq}\right\}}_{m=1}^d$ are then:

$n=1-\frac{4}{\pi }\arctan \left[\frac{\Re \left(\log {\stackrel{⌣}{z}}_m^{2^{\frac{1+n}{2}}-1}\right)}{\Im \left(\log {\stackrel{⌣}{z}}_m^{2^{\frac{1+n}{2}}-1}\right)}\right]$E31

$T_m^{eq}=-\left[\frac{\Im \left(\log {\stackrel{⌣}{z}}_m^{2^{\left(1+n\right)/2}-1}\right)}{2\pi f_r\left(2^{\left(1+n\right)/2}-1\right){\left(\frac{f_s}{f_r}\right)}^{\left(1+n\right)/2}}\right]$E32

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

$c={\left({\tilde{A}}^H\tilde{A}\right)}^{-1}{\tilde{A}}^H\tilde{s}.$E33

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 $\tilde{s}$ provided by the spline interpolant at frequencies selected according to Equation (28) to construct a Hankel matrix $\tilde{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 ${\left\{T_m^{eq}\right\}}_{m=1}^d$ 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., $\Delta n\ne 0$, 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., $\Delta 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., ${10}^{-6}$.

4. Step 4. Use the estimates of $n$ and ${\left\{T_m^{eq}\right\}}_{m=1}^d$ 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 $d$echoes 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 $\left(T_m^{eq}, n_m^{eq}\right)T_m^{eq}$ is given by Equation (25), and $n_m^{eq}$ is given by the following equation [31]:

$n_m^{eq}=1-\frac{4}{\pi }\text{arctan}\frac{{{\displaystyle \sum }}_{l=1}^m{\tau }_l\text{tan}\left[\frac{\pi }{4}\left(1-n_l\right)\right]}{T_m^{eq}}.$E34

The equivalent‐medium approximation then allows writing:

$A\approx \left[\begin{array}{llll}{\widehat{a}}_1\hfill & {\widehat{a}}_2\hfill & \dots \hfill & {\widehat{a}}_d\hfill \end{array}\right]$E35

where

${\widehat{a}}_{\text{m}}={\left[\begin{array}{llll}{\widehat{\text{z}}}_{\text{m}}^{1^{(1+{\text{n}}_{\text{m}}^{\text{eq}})/2}}\hfill & {\widehat{\text{z}}}_{\text{m}}^{2^{(1+{\text{n}}_{\text{m}}^{\text{eq}})/2}}\hfill & \dots \hfill & {\widehat{\text{z}}}_{\text{m}}^{{\text{N}}^{(1+{\text{n}}_{\text{m}}^{\text{eq}})/2}}\hfill \end{array}\right]}^{\text{T}}$E36

${\widehat{\text{z}}}_m=\text{exp}\left\{-2\pi f_r{T}_m^{eq}\left(j+\text{tan}\left[\frac{\pi }{4}\left(1-n_m^{eq}\right)\right]\right){\left(\frac{f_s}{f_r}\right)}^{(1+n_m^{eq})/2}\right\}.$E37

The above approximation becomes exact for a DI structure, i.e., ${\widehat{\text{z}}}_{\text{m}}={\stackrel{⌣}{z}}_{\text{m}}$. Note also that upon setting${\left\{n_l=1\right\}}_{l=1}^d$, 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 ${\left\{s(k)\right\}}_{k=1}^N$ estimate the time delays ${\left\{{\tau }_l\right\}}_{l=1}^d$ and the dispersion indices ${\left\{n_l\right\}}_{l=1}^d$ 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 $n_1^{eq}$ and $n_2^{eq}$ are known ($n_1^{eq}$ is practically equal to $n_1$) and consider the application of the nonlinear transformation in Equation (28) to the uniform frequencies. Setting $n=n_1^{eq}$ in $\nu \left(k\right)$ and probing the medium with ${\left\{{\nu }_k\right\}}_1^N$ give rise to the following discrete backscattered signal:

$\tilde{s}(k)=c_1{\widehat{z}}_1^{{\left[\frac{{\nu }_k}{f_s}\right]}^{(1+n_1^{eq})/2}}+c_2{\widehat{z}}_2^{{\left[\frac{{\nu }_k}{f_s}\right]}^{(1+n_2^{eq})/2}}+\tilde{b}(k)$E38

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:

$\tilde{s}=\tilde{A}c+\tilde{b}$E39

with the following notational definitions:

$\begin{array}{c}\tilde{s}={\left[\begin{array}{llll}\tilde{s}(1)\hfill & \tilde{s}(2)\hfill & \dots \hfill & \tilde{s}(\text{N})\hfill \end{array}\right]}^{\text{T}}\\ \tilde{A}=\left[\begin{array}{cc}{\tilde{a}}_1& {\tilde{a}}_2\end{array}\right]\\ {\tilde{a}}_1={\left[\begin{array}{ccc}{\widehat{\text{z}}}_1{\widehat{\text{z}}}_1^2{}^{\left(1+n_1^{eq}\right)/2}& \dots & {\widehat{\text{z}}}_1^{\left(N-1\right)\times 2^{\left(1+n_1^{eq}\right)/2}-\left(N-2\right)}\end{array}\right]}^T\\ {\tilde{a}}_2={\left[\begin{array}{ccc}{\widehat{\text{z}}}_2{\widehat{\text{z}}}_2^{2\left(1+n_2^{eq}\right)/2}& \dots & {\widehat{\text{z}}}_2^{{\left[\left(N-1\right)\times 2^{\left(1+n_1^{eq}\right)/2}-\left(N-2\right)\right]}^{\left(1+n_2^{eq}\right)/\left(1+n_1^{eq}\right)}}\end{array}\right]}^T\\ \tilde{b}={\left[\begin{array}{ccc}\tilde{b}\left(1\right)\tilde{b}\left(2\right)& \dots & \tilde{b}\left(N\right)\end{array}\right]}^T\end{array}$E105

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

${\tilde{S}}_{11}={\text{c}}_1{\tilde{z}}_{11}{\tilde{z}}_{21}\text{and}{\tilde{S}}_{12}={\text{c}}_1{\tilde{z}}_{11}{\tilde{z}}_1^{2^{\left(1+n_1^{eq}\right)/2}-1}{\tilde{z}}_{21}$E40

where

$\begin{array}{ll}{\tilde{z}}_{11}\hfill & ={\left[\begin{array}{llll}{\widehat{z}}_1\hfill & {\widehat{z}}_1^{2^{\left(1+n_1^{eq}\right)/2}}\hfill & \dots \hfill & {\widehat{z}}_1^{\left(N-L-1\right)\times 2^{\left(1+n_1^{eq}\right)/2}-\left(N-L-2\right)}\hfill \end{array}\right]}^T\hfill \\ {\tilde{z}}_{21}\hfill & =\left[\begin{array}{llll}1\hfill & {\widehat{z}}_1^{2^{\left(1+n_1^{eq}\right)/2}-1}\hfill & \dots \hfill & {\widehat{z}}_1^{\left(2^{\left(1+n_1^{eq}\right)/2}-1\right)\left(L-1\right)}\hfill \end{array}\right]\hfill \end{array}$E106

The corresponding matrix pencil can then be written as follows:

${\tilde{s}}_{12}-\text{λ}{\tilde{s}}_{11}=c_1{\tilde{z}}_{11}\left[{\widehat{z}}_1^{2^{\left(1+n_1^{eq}\right)/2}-1}-\text{λ}\right]{\tilde{z}}_{21}$E41

As in conventional MPM, the value of $\text{λ}\text{=}{\widehat{z}}_1^{2^{\left(1+n_1^{eq}\right)/2}-1}$ 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 $\text{λ}$ as the GEV of the total matrix pair is valid provided that the contribution of other echoes is filtered out after each restoration. Therefore, $\text{λ}$ 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 $n_1^{eq}$ and $T_1^{eq}$ can then be deduced from $\text{λ}$ as follows:

$n_1^{eq}=1-\frac{4}{\pi }\text{arctan}\left[\frac{\Re (\text{log }{\widehat{z}}_1^{2^{\left(1+n_1^{eq}\right)/2}-1})}{\Im (\text{log}{\widehat{z}}_1^{2^{\left(1+n_1^{eq}\right)/2}-1})}\right]$E42

$T_1^{eq}=-\left[\frac{\Im (\text{log} {\widehat{z}}_1^{2\left(1+n_1^{eq}\right)/2_{-1}})}{2\pi f_r(2^{\left(1+n_1^{eq}\right)/2}-1){\left(\frac{f_s}{f_r}\right)}^{\left(1+n_1^{eq}\right)/2}}\right]$E43

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 $\tilde{S}={\tilde{S}}_1+{\tilde{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 ${\tilde{a}}_1$ as ${\tilde{S}}_1=U_1{\Sigma }_1{V}_1^H$. Using MATLAB notation, the orthogonal complement is then given by $U_1^{\text{ort}h}=U_1(:,\text{rank}({\Sigma }_1)+1:\text{end})$ The definition of $U_1^{\text{ort}h}$ reflects the fact that the rank of ${\tilde{S}}_1$ or, equivalently, ${\Sigma }_1$ would be overestimated, unless ${\tilde{a}}_1$ admits a Vandermonde structure guaranteed by setting $n=n_1^{eq}$ in Equation (28). More importantly, such a definition of $U_1^{\text{ort}h}$ enables, once $n_1^{eq}$ and $T_1^{eq}$ are estimated, to filter out the contribution of ${\tilde{a}}_1$ even if the rank is overestimated or, equivalently, $n\ne n_1^{eq}$. 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 $\tilde{S}$ by the conjugate transpose of $U_1^{\text{ort}h}$ gives ${\tilde{S}}^{\text{pr}}$:

${\tilde{S}}^{\text{pr}}={\left(U_1^{\text{ort}h}\right)}^{\text{H}}\tilde{S}={\left(U_1^{\text{ort}h}\right)}^{\text{H}}\left({\tilde{S}}_1+{\tilde{S}}_2\right)={\left(U_1^{\text{ort}h}\right)}^{\text{H}}{\tilde{S}}_2$E44

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 ${\tilde{S}}^{\text{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., $U_1^{\text{ort}h}=U_1(:,2+1:\text{end})$. Similar to Equation (38), setting $n=n_2^{eq}$ in Equation (28) restores the Vandermonde structure, only this time for the second mode vector. The estimates of $n_2^{eq}$ and $T_2^{eq}$ are then determined from the GEV of the matrix pair $\left[{\tilde{S}}_2^{\text{pr}},{\tilde{S}}_1^{\text{pr}}\right]$ 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 ${\left\{n_m^{eq}\right\}}_{m=1}^d$ and ${\left\{T_m\right\}}_{m=1}^d$ then allow determining the complex amplitudes by means of a LS fit having the following solution:

$c={({\tilde{A}}^{\text{H}}\tilde{A})}^{-1}{\tilde{A}}^{\text{H}}\tilde{s}.$E45

Note that the development of the algorithm has till now been based on the a priori knowledge of $n_1^{eq}$ and $n_2^{eq}$. 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]:

1. Step 1. Form a Hankel data matrix $\tilde{S}$ from a data sequence corresponding to frequencies given in Equation (28). If the model order $M<d$ (d is known), left multiply $\tilde{S}$ by the conjugate transpose of $U_{\text{m}}^{\text{ort}h}$. Start with $M=d$ and $n=1$ as initial values, which correspond to the exact uniform data sequence provided by real measurements.

2. Step 2. Run MPM and determine ${\left\{n_m^{eq}\right\}}_{m=1}^d$ and ${\left\{T_m^{eq}\right\}}_{m=1}^d$ by using Equations (42) and (43).

3. Step 3. Compare the estimated value of $n_m^{eq}$ with the value used in step 1. If they are different, i.e., $\Delta n\ne 0$ substitute the new estimate of $n_m^{eq}$ 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., $\Delta 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., ${10}^{-6}$

4. Step 4. Retain the estimates of $n_m^{eq}$ and $T_m^{eq}$ of the last iteration and use them each time $U_{\text{m}}^{\text{ort}h}$ is to be determined from the mode vector ${\tilde{a}}_m$ for a given value of n.

5. Step 5. Set the model order $M$ to $d-1$ 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.

6. Step 6. Use the estimates of ${\left\{n_m^{eq}\right\}}_{m=1}^d$ and ${\left\{T_m\right\}}_{m=1}^d$ to determine the complex amplitudes using Equation (45).

Figure 6.

MSE of n₁ and n₂ versus bandwidth.

Figure 7.

MSE of τ₁ and τ₂ 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, and$B{\tau }_{\text{2}}$. Then, it examines the performance of RAP‐MPM in the case of a DV structure. We consider a reference frequency $f_r=1.5\text{GHz}$ and a frequency shift$f_s=10\text{MHz}$. The frequency‐independent complex amplitudes are $c_1=0.5\text{exp}(j0.7)$ and$c_2=1-c_1^2$. Unless otherwise defined, SNR represents the average SNR over the data interval. The convergence threshold for both algorithms is fixed at${10}^{-6}$, 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 about$f_{\text{r}}$. Figures 6 and 7 show the MSEs of the parameters for two ($n_1=n_2=0.83,{\tau }_{\text{1}}=\text{ 3}.\text{2 ns},{\tau }_{\text{2}}=\text{ 1}.\text{2 ns}$ and $\text{SNR}\text{=}\text{30}\text{dB}$) 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 ($n_{\text{1}}=0.\text{87},n_{\text{2}}=0.\text{83},{\tau }_{\text{1}}=\text{ 3}.\text{2 ns},{\tau }_{\text{2}}=\text{ 1}.\text{2 ns}$ and $\text{SNR}\text{=}\text{30}\text{dB}$), 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 n₁ and n₂ versus bandwidth.

Figure 9.

RAP‐MPM MSEs of τ₁ and τ₂ 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 ($n_1=n_2=0.83,{\tau }_{\text{1}}=\text{ 3}.\text{2 ns},{\tau }_{\text{2}}=\text{ 1}.\text{2 ns}$, and $\text{BW}\text{=}\text{2}\text{GHz}$) 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 ($n_1=0.87, n_2=0.83,{\tau }_1=3.2\text{ns,}{\tau }_2=1.2\text{ns}$, and $\text{BW}\text{=}\text{2}\text{GHz}$), the MSEs of the indices and the time delays approach their CRLBs starting, respectively, from $\text{SNR}\text{=}\text{5}\text{and 10}\text{dB}$ as shown in Figures 12 and 13.

Figure 10.

MSE of n₁ and n₂ versus SNR.

Figure 11.

MSE of τ₁ and τ₂ versus SNR.

Figure 12.

RAP‐MPM MSEs of n₁ and n₂ versus SNR.

Figure 13.

RAP‐MPM MSEs of ${\tau }_1$ and ${\tau }_2$ versus SNR.

3.6. RRMSE versus $B{\tau }_2$

This simulation examines the resolution capability of the algorithms by varying ${\tau }_{2 }$ at $\text{BW }=\text{ 2 GHz}$ so that $0.2\le B{\tau }_2\le 2$. Figure 14 shows the RRMSE of the estimates of ${\tau }_{2 }$ as a function of $B{\tau }_2$ along with the CRLB and a 5% error threshold for $n_1=n_2=0.83,{\tau }_{\text{1}}=\text{ 3}.\text{2 ns}$, and $\text{SNR}\text{=}\text{30}\text{dB}$. 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{\tau }_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{\tau }_2<1$, which indicates their high‐resolution capability.

Figure 14.

The RRMSE versus Bτ₂.

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τ₂. 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.

Acknowledgments

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