Full Wave Modeling of Electromagnetic Scattering by an Object Buried between Two Rough Surfaces: Application to GPR

|In this paper, we present an eﬃcient numerical method to calculate the frequency and time responses of the (cid:12)eld scattered by an object buried between two random rough surfaces for a 2-D problem. This method is called Generalized PILE (GPILE) method because it extends the PILE method which considers only two surfaces or an object buried under a surface. The GPILE method solves the Maxwell equations rigourously by using a simple matrix formulation. The obtained results have a straightforward physical interpretation and allow us to investigate the in(cid:13)uence of the object buried between the two rough surfaces. We distinguish the primary echo of the upper surface, the multiple echoes coming from the lower surface, and those arising from the object. The GPILE method is applied to simulate the Ground Penetrating Radar (GPR) signal at nadir. The resulting time response helps the user to detect the presence of the object buried between the two random rough surfaces.


INTRODUCTION
Electromagnetic scattering from an object buried under a rough surface has attracted much interest for several decades. In order to solve this problem, some asymptotic and rigorous models have been developed. The asymptotic (or approximate) models are based on simplifying assumptions. Rigorous models are less restrictive and thus find a large field of validity. Numerical methods for rigorous models can be classified into two families: differential equations methods, where discretization uses Finite Element Method (FEM) [1,2] or Finite Difference Time Domain (FDTD) method [3,4] and the integral boundary methods, where the discretization uses the Method of Moments (MoM) [5,6]. The MoM offers several advantages over FDTD. Indeed, the MoM discretizes the boundaries (here the interfaces), while the FDTD discretizes the space. For rough interfaces, the mesh grid must be refined in the neighbourhood of the interfaces for the FDTD in order to follow the surface profiles, which increases the memory space requirements, unlike for the MoM. The main drawback of MoM is that it can be applied only to interfaces separated by homogenous media and for a single frequency. Thus, the calculation of a time response requires to apply MoM N f times, where N f is the number of frequencies.
In the case of a single rough surface, fast methods for solving the linear system obtained by the MoM have been developed in order to reduce the number of operations as well as the storage memory space. We can cite, for example, the Method of Ordered Multiple Interactions (MOMI) [7]. The major drawback of this method is its slow convergence for the dielectric case. The Forward-Backward (FB) method [8] adapts better to dielectric interfaces [9]. An accelerated version of the FB uses the spectral representation of Green's function (FB-SA) [10]. The Banded-Matrix-Iterative-Approach (BMIA) method [11,12] aims to accelerate the matrix-vector products. An improved version of this method gave rise to the Banded-Matrix-Iterative-Approach/CAnonical Grid (BMIA-CAG) [13].
In the case of scattering by two rough surfaces, some methods have been devoted to obtaining a rigorous solution: the Extended Boundary Condition Method (EBCM) [14], Forward-Backward with Spectral Acceleration (FBSA) method [15], and Steepest Descent Fast Multipode Method (SDFMM) [16]. However, all these methods have some constraints. In contrast, the Propagation Inside-Layer Expansion (PILE) method [17] is rigorous, of a simple mathematical formulation, and has an intuitive physical interpretation. This method has been accelerated by using the fast algorithm BMIA/CAG to give the PILE-BMIA/CAG method [18]. Another acceleration of the PILE method adapts the spectral acceleration (SA), both to local interactions and to the coupling steps and gives the PILE-FBSA method [19]. The PILE method has been adapted to scattering from an object buried under a rough surface [20]. In the same spirit, the PILE method has been extended to a more general case of the scattering from two illuminated scatterers and gives Extended PILE method [21]. A generalization of the PILE method (GPILE: Generalized Propagation Inside-Layer Expansion) has been presented in [22], to calculate the electromagnetic field scattered by three superimposed surfaces.
In this paper, we consider an object buried inside a layer with two rough interfaces for a 2-D problem. The main difference with the problem of three superimposed surfaces treated in [22] is that the object and the lower surface are both illuminated by the wave transmitted by the upper surface, which constitutes an extension of the previous problem. We assume that all media are homogeneous and construct a rigorous model for this problem from an integral boundary method in the frequency domain, where the discretization uses the MoM. To reduce the complexity of the inversion of the impedance matrix obtained by the MoM, we extend the PILE algorithm [17] to the case of three scatterers by forming a composite scatterer made up of the object and lower surface. The purpose of this paper is to present detailed numerical studies of the GPILE method for an object buried between two rough surfaces. The results are presented both in the frequency domain and in the time domain to highlight the variations in the amplitude of the signal over a given period of time. This will then allow us to distinguish the primary echo from the upper surface and the multiple echoes from the object and/or the lower surface. The study is expected to give a better understanding of the propagation phenomena in complex media such as that composed of an object buried between two rough surfaces.
In what follows, we describe the geometry of the problem in Section 2. Section 3 presents the GPILE method based on the PILE method. Then, we establish detailed calculations from this method. Section 4 presents some numerical simulations for a single frequency. Section 5 provides simulations for Ground Penetrating Radar (GPR) signals for a realistic scenario for the problem under study. Fig. 1 two random rough surfaces S 1 and S 3 invariant along the y axis and an object delimited by the surface S 2 buried between the interfaces S 1 and S 3 . The two surface heights are generated by a Gaussian probability density function (PDF) and a Gaussian autocorrelation function (ACF), and we assume that they do not intercept each other or with the object. The three interfaces separate four homogeneous media: the upper medium, Ω 0 , assumed to be vacuum, the intermediate medium, Ω 1 , which constitutes a layer, the buried object which constitutes the medium Ω 2 , and the lower medium, Ω 3 , which we consider to be dielectric. Consider an incident wave ψ inc (r) in the plane (x,ẑ), at an incidence angle θ inc defined relatively to the axis z counterclockwise, as shown in Fig. 1.

Consider in
Since surfaces have a finite length, edge diffraction occurs because the incident field does not vanish at the edges of the surfaces. To reduce this phenomenon in the study of scattering from rough surfaces, a tapered incident wave must be used instead of a plane incident wave. Here, we will apply the Thorsos wave, as described in [23], with an attenuation parameter g which controls the extent of the incident wave.

Impedance Matrix
The integral equation method consists in evaluating the currents ψ i and their normal derivatives ∂ψ i /∂n i on each surface S i (i = {1, 2, 3}) and in deducing the scattered fields in each medium by using Huygens' θ sca ψ sca, 2 ψ sca, 3 ψ sca, 4 n1n 2n 3 Figure 1. Electromagnetic scattering from an object buried between two 1-D rough surfaces (2-D problem: plane (x,ẑ)).
principle. From the Method of Moments, the boundary integral equations are discretized on each surface of the scatterer, leading to the linear systemZX = b, where the impedance matrix is The impedance matrix is of size 2( where N i is the number of samples on S i . All the impedance matrices are also block matrices and are detailed in Appendix A. We will see later that these matrices are simplified in the perfect conductor cases. It is valuable to notice that all the coupling matrices are non-zero, since all the interfaces have direct interaction between them, unlike the case of three interfaces treated in [22]. The unknown vectors X i containing the surface currents and their normal derivatives are written as where r p∈{1;N i } ∈ S i , i = {1, 3} for the two surfaces, and r p∈{1;N 2 } ∈ S 2 for the object. The term b is the source and contains information on the incident field: In the absence of one of the scatterers, the impedance matrix reduces to the case of the other two scatterers. This remark establishes a theoretical and rapid validation of the MoM compared to existing work on electromagnetic scattering, in the case of two rough surfaces [17], an object below a rough surface [20], and an object above a rough surface [21].

Perfect Conductor Cases
We can consider the cases where the object and/or the lower surface are perfect conductors. In these cases, the impedance matrices are simplified. If the object is a perfect conductor, then ψ 2 vanishes on the interface of the object in Transverse Electric polarization (TE, Dirichlet boundary conditions), and the only unknown is ∂ψ 2 /∂n 2 . On the other hand, in Transverse Magnetic polarization (TM, Neumann boundary conditions), ∂ψ 2 /∂n 2 vanishes on the interface of the object, and the only unknown is ψ 2 .
On the other hand, if the lower surface is a perfect conductor, then ] , ] T , ] , From these two cases, we can find the case where both the object and the lower surface are perfect conductors.

GPILE
To reduce the complexity of the inversion of the impedance matrix obtained by the MoM, we extend the PILE algorithm [17] to the case of three scatterers by considering on the one hand the upper surface as scatterer 1, and on the other hand the object and the lower surface as a single scatterer, which we call composite scatterer 2. Let us rewrite impedance matrices for this structure as: ) , The PILE method in the case of two scatterers consists in inverting the impedance matrix by decomposition of domains from the Taylor series expansion of the inverse of the Schur complement [24]:Z where the matricesT,Ū,V, andW are expressed as functions of four block matrices ofZ as follows The PILE method then makes it possible to obtain the surface currents Y 1 and Y 2 on scatterer 1 and composite scatterer 2 (the object and the lower surface) respectively: (11) is the coupling matrix of the composite scatterer 2 towards the scatterer 1 and thus brings the contributions of the composite scatterer 2 to the scatterer 1.
The convergence of this method is based on the value of the spectral norm ofN c . Denoting ∥N c ∥ sr , this norm is always strictly smaller than 1 as shown in [17], which ensures the convergence of the algorithm.
We now want to detail the calculations of the currents on each surface in the spirit of the original PILE method. To do this, we denote the surface currents on the upper surface, the object and the lower surface by X 1 , X 2 , and X 3 , respectively. Note that X 1 = Y 1 , because scatterer 1 in (10) corresponds to the upper surface. It is then sufficient to detail the calculations of the characteristic matrixN c to obtain an expression equivalent to X 1 , with the only issue being to invert the block matrixP 22 as shown in (6). By reusing the inverse of the Schur complement algorithm, we obtain where the matricesT,Ū,V, andW are given by (9) as functions of four blocks of the matrixP 22 .
Recall in addition that, as the matrixN c is the sum of the contributions of the object and of the lower surface on the upper surface, we can explicitly highlight the two contributions by noting: We will see, after developing these two expressions, that the matricesN c,21 andN c,31 are characteristic matrices for the principal interactions (those coming from the upper surface towards the object or the lower surface) and bring back the contributions of the object and the lower surface towards the upper surface, respectively. By injecting relations (9) into the characteristic matrices (13) and (14), we obtain their approximations by the Taylor series expansion (with truncation to the order Q) (17) is the characteristic matrix which brings back the contributions from the lower surface to the object, taking into account a secondary interaction at the order q. The matrixZ ′ 12 =Z 32Z −1 33Z 13 represents the coupling matrix between the upper surface and the object via the lower surface.
We present in Appendix B the calculations of surface currents by the substitution method. An equivalent expression to the characteristic matrix (16) obtained by this method is and the matrixZ ′ 13 =Z 23Z −1 22Z 12 represents the coupling matrix between the upper surface and lower surface via the object.
Expression (18) brings out the symmetry of the problem with respect to (15) instead of (16). It will be used for the rest of the paper. Moreover, the contribution matrix of the object and the lower surface each consists of two contributions, with one of them being direct and the other being indirect. It is the coupling matricesZ 12 andZ ′ 12 which make it possible to quantify these contributions for the object. The coupling matricesZ 13 andZ ′ 13 play the same role for the lower surface. By regrouping in one expression the two contributions to calculate the total currents on the upper surface, we obtain: Finally, It is important to note that the orders P and Q are arbitrary. In fact, according to the construction of the method (see expressions (10)), the main order p represents the round trips of the wave between the upper surface and the object or the lower surface. As for the secondary order q, it represents the round trips of the wave between the object and the lower surface, and vice versa.
The calculations of the currents on the other interfaces are done by the same types of algebraic manipulations, and we obtain for the object from where and for the lower surface from where The calculations of the currents on the lower surface by applying the substitution method (see Appendix B) make it possible to obtain The expressions of the currents on the object (23) and on the lower surface (26) have a physical interpretation in accordance with the EPILE method [21]. Remark 1. We can calculate the currents on the upper surface according to the expression (21) and calculate them on the object and the lower surface by using the EPILE algorithm [21]. By considering that the object and the lower surface are illuminated by the fields −Z 12 X 1 and −Z 13 X 1 respectively, we obtain relations (23) and (26). The sign (−) in the fields comes from the fact that the two scatterers are in the transmitted medium with respect to the upper surface.
Remark 2. The fundamental difference between (10) and (21) is that in (10), Q = ∞, because the impedance matrixP 22 (for the lower surface and the object) is inverted directly, while in (21), c , and the expression (21) allows us to study numerically the contributions of the object and of the lower surface for each scattering.
Remark 3. It is valuable to note that if instead of an object we had a surface, there would not be direct interaction between the upper surface and lower surface, which would result inZ 13 =0 (and thereforeZ 31 =0). It would follow thatN c,31 =0, and one would find a GPILE method equivalent to that developed in [22] for the case of three superimposed interfaces.
For the numerical simulations, we implement expression (21) to calculate the currents on the upper surface and expressions (23) and (26) for the calculation of the currents on the object and on the lower surface, respectively. The characteristic matrices will be approximated by expressions (15) and (18). To avoid calculating products of matrices when approximating the characteristic matrices, we will start the process by calculating the currents and their normal derivatives on the upper surface for P = 0, which amounts to calculating a matrix-vector product. Having the current vector on the upper surface at this order, we can then calculate it for the order P = 1 by matrix-vector products and so on. This procedure is well suited for calculating currents for any value of Q. In this way, the currents on the object and on the lower surface will also be calculated for any values of P and Q by matrix-vector products.

NUMERICAL RESULTS FOR A SINGLE FREQUENCY
The GPILE method was derived from the PILE method by coupling the upper surface, and the object plus the lower surface as a single composite scatterer. Thus, by decoupling the object and the lower surface, we can numerically study the impact of the orders p and q on the convergence of the algorithm, as well as the contributions of the object and of the lower surface on the scattered field in the incidence medium. In addition, the orders p and q make it possible to analyse the influence of the primary echo of the upper surface and the multiple echoes in the intermediate medium.

Numerical Parameters
Consider an elliptical cylinder of semi-major axis a = 100 mm and semi-minor axis b = 10 mm buried between two rough surfaces of length 3000 mm, which are generated by a Gaussian probability density function (PDF) and a Gaussian autocorrelation function (ACF). Moreover, the profiles of two generated surfaces do not have to intersect each other. So, the distance between either the upper or lower surface and the object has been chosen to avoid this possibility. The same caution is used for the two interfaces. Numerically, during the program execution, a check is made before the electromagnetic computation. The three interfaces separate four homogeneous media (Ω 0 , Ω 1 , Ω 2 , Ω 3 ) as described in Section 2. The standard deviations of the heights of the two random rough surfaces are σ h,1 = 1 mm and σ h,3 = 2.5 mm, respectively. Their correlation lengths are L ch,1 = 15 mm and L ch,3 = 30 mm. The mean thickness between the two surfaces is H = 40 mm. In this paper, we consider only one realisation of the two surfaces. The cylinder is centered at the point of coordinates C = (0, −0.5H). Fig. 2 plots the heights (in mm) versus their abscissa (in mm). The medium relative permittivities are ϵ r,0 = 1, ϵ r,1 = 2, ϵ r,2 = 7, and ϵ r,3 = 4. As we will focus on GPR applications, we consider a Thorsos wave at normal incidence (θ inc = 0 • ) of frequency f = 2 GHz whose attenuation parameter is g = L/6. The number of points n i (i = {1, 2, 3}) per wavelength is 10, and the wavelength in vacuum is λ 0 = 150 mm.

Convergence of GPILE: Impact of Orders P and Q
To test the convergence of the GPILE algorithm, we plot the scattered field and calculate the residual relative error (RRE). In fact, the MoM with the classical lower-upper (LU) inversion algorithm contains an infinite number of interactions among the three scatterers (p = q = ∞), whereas for GPILE we only consider a limited number of interactions. We have defined two types of interactions, those called principal (between the upper surface and lower surface or the object) and those called secondary (between the object and the lower surface). In this way, the orders p and q of GPILE are the principal order and the secondary order, respectively. The residual relative error is calculated by where the norm of a vector of components x i and of length N is given by The vector X can represent the surface currents ψ i or their normal derivatives ∂ψ i /∂n i calculated by GPILE on all surfaces. In the context of this study, X will rather represent the scattered field. The LU label means that the vector is calculated by the LU inversion algorithm. We present the following figures where we compare both the GPILE and LU methods as well as the PILE methods for the cases of two superimposed surfaces or an object buried under a surface. This allows us to study the performance of GPILE and also to show the contributions of the two buried scatterers. We present the moduli of the scattered fields in decibels (10 log 10 |ψ sca |) with respect to the observation abscissa (x obs ∈ [−L/2; L/2]). In the legend, S+O+S represents the case of an object buried between two surfaces, while S+O represents an object buried under a surface, and S+S the case of two superimposed surfaces. Then, we have the orders P and Q (or only P for PILE). The number in parenthesis is the RRE with respect to the reference LU method.
We first look at the impact of the principal order p. Thus, in Fig. 3, the PILE and GPILE methods are presented for P = 0. We observe that the scattered fields calculated by the three methods are exactly the same, because at this order, the three methods calculate only the first echo, with account of only the upper surface. In Figs. 4(a)-(b), P = 1, differences among the three methods appear because for the PILE methods, the upper surface only interacts with the lower surface or the object, but for GPILE, the upper surface interacts with both the lower surface and the object. This makes it possible to highlight the presence of the object between the two surfaces. In Figs. 4(c)-(d), with P = 2, we observe a perfect agreement between the LU method and GPILE method when Q = 1, which shows that P = 2 is sufficient to take account of all the significant contributions. Note also that P = 2 means that the wave makes two round trips between the upper surface and lower surface or the object.
To evaluate the impact of the secondary order q, several numerical simulations, which we do not present here, show that Q = 0 or 1 is sufficient for the convergence of the algorithm when P ≥ 1 and that the higher orders of q do not have a significant impact on the precision of GPILE. Besides, Q = 0 means that the wave does not perform a round trip between the object and the lower surface, and vice versa. However, even for Q = 0, there are interactions between these two scatterers viaZ ′ 12 andZ ′ 13 as shown by relations (15) and (18)

Variations of Some Parameters
We show in Fig. 5 the impact of some parameters of the object on the scattered field. Fig. 5(a) plots the case of a large elliptical cylinder with semi-major axis a = 1000 mm while Fig. 5(b) plots the case of a small elliptical cylinder with semi-major axis a = 10 mm. The semi-minor axis is b = 10 mm for the two cylinders. We can see on the levels of curves as well as on the RRE that the contribution of the object varies proportionally to its size. We also observe that the scattered field of an object buried between two rough surfaces tends towards that of two rough surfaces when the object becomes smaller. Likewise, regarding the relative permittivities, we show in Figs. 5(c)-(d) that the contributions of the object or the lower surface are much more important when these two scatterers are perfect conductors. Remark also that, in all these cases, the GPILE method converges for P = 2 and Q = 1, despite changing parameters. This remark also holds for the simulations presented in Section 5, in which the frequency band is [0.2-6.0] GHz. In the following, we consider the GPILE method for P = 2 and Q = 1 as a reference for the study of electromagnetic scattering by an object buried between two random rough surfaces.
In addition, it has already been established in [17] that the scattered field in the case of two superimposed rough interfaces depends on the geometry and physical parameters, such as the permittivities of the media, the statistical parameters of the interfaces, and the thickness between these interfaces. Turning now to the study of an object buried between two rough surfaces, several simulations that we do not present here show that in addition to the surface parameters, the scattered field also depends on the parameters of the object such as its size, shape, position, and permittivity. In this way, the sensitivity of the scattered field to the variation of these parameters could allow the characterization of the object in an inverse problem.

FREQUENCY AND TIME RESPONSE RESULTS
In this section, PILE and GPILE methods are used to provide simulations of GPR data for an object buried between two random rough surfaces. To compute the GPR responses, both methods are performed on N f equally spaced frequencies covering the GPR band.
A Ricker pulse is considered as an input signal [25]. In the time domain, it is defined as  In the Fourier domain, it is defined as The absolute value of the spectrum is maximum for f = f c , and the value of the maximum is s 0 = 2/(ef c √ π). For our simulations, we take f c = 2 GHz. Fig. 6 plots the Fourier transform versus the frequency f and the signal versus the time t. As shown in Fig. 6(a) Fig. 7(a) shows the input signal and scattered field. In Fig. 7(b), the GPILE and PILE methods give the same results because for P = 0, the three methods are calculated for the first echo only, i.e., the incident wave interacts only with the upper surface. In Figs. 7(c)-(d), the difference among the three methods appears because for the PILE methods, the upper surface interacts only with the lower surface or the object, but for the GPILE method, the upper surface interacts with both the lower surface and the object. On the other hand, the unremarkable difference between these two figures shows that the third echo is weak. We will see later that this is the case.

Echo Times
We plot in Fig. 8 the time responses of the scattered field. The simulations parameters are the same as in Fig. 7. For smooth interfaces, the first three echoes of the lower surface arise at the times T 1 , T 2,S , and T 3,S defined as: where H is the thickness of the layer. We can obtain the second and third echo times for the object by replacing H by H 12 in (31), where H 12 is the depth of the object. Those echoes are named T 2,O and T 3,O . Note that the depth of the object can be calculated in three different ways, first relative to the center of the object, second relative to the upper vertex, and finally relative to the lower vertex. Each way of calculating this depth gives a different result. In this paper, we will focus on calculating the depth of the object from the center. Remember that these echoes are calculated by the PILE methods. Regarding the GPILE method, the second and third echoes result from multiple scattering between the object and the lower surface in the medium Ω 1 . Thus, for P = 1 and Q = 0, the second echo is made up of two contributions: one from the object and the other from the lower surface. Each of these contributions has two components: one direct and the other indirect. For the object for example, we have seen in relation (15) that the contribution matrix is composed by the direct contribution, represented byZ 12 (taking account for the wave coming from the scattering on the object of the wave transmitted by the upper surface) and the indirect contribution represented byZ ′ 12 (taking account for the wave coming from the scattering on the object of the wave coming from the scattering on the lower surface of the wave transmitted by the upper surface). For the lower surface, the direct and indirect contributions are highlighted in relation (18) and represented by the contribution matricesZ 13 andZ ′ 13 , respectively. Moreover, the direct contribution echoes of the object and of the lower surface arise at different times, because for the first one the wave goes back and forth to the object while for the second one the wave goes back and forth to the lower surface. In fact, the indirect contribution echoes of the two scatterers arise at the same time (see Ind. Object and Ind. Lower Surf. in Fig. 9). For the object, the wave goes first to the lower surface, then passes through the object before returning to the upper surface. For the lower surface, the wave passes through the same path but in the opposite direction. Theoretically, the indirect contribution echoes arrive at the same time as the direct contribution echo from the lower surface. However, in practice a time lag will be created for the indirect contribution echoes because of the scattering from the object. For Q = 1, it suffices to add the time of a round trip between the object and the lower surface, given by T OS = T 2,S − T 2,O .

Signal Analysis
In Fig. 8, PILE methods show the impact of the object and of the lower surface while GPILE method takes also account the multiple interactions between the two buried scatterers. As can be seen in Figs. 8(a) and 8(b), the first echo is correlated to the input signal. Fig. 8(c) shows that the second echo of the object has two peaks between [3.5-4.5] ns. The first peak comes from the signal emanating directly from the antenna while the second peak comes from the closed shape of the object and the roughness of the upper surface. The second echo of the lower surface has a peak between [5-5.5] ns. This peak is followed by a few low intensity peaks linked to the roughness of the two surfaces. The times calculated for the second echoes from the two buried scatterers are consistent with the observed times. We calculated the depth of the object in relation to its center, and we obtained T 2,O = 3.98 ns, which corresponds to the point of zero amplitude (inflexion point). However, if we calculated the depth from the upper vertex, we would have T 2,O = 3.88 ns, which corresponds to the maximum amplitude peak. On the other hand, if we calculated the depth from the lower vertex, we would have T 2,O = 4.07 ns, which corresponds to the minimum amplitude peak.
The second echo for GPILE 1/1 is composed of the sum of the direct and indirect contribution echoes of the two buried scatterers, and one interaction between these scatterers. Observe that GPILE and PILE(S+O) coincide around T 2,O as shown in Fig. 8(c), because the direct contribution echo for the object for GPILE is exactly PILE(S+O). However, for the lower surface, there is a time lag between PILE(S+S) and GPILE around T 2,S due to the scattering from the object. There is also a difference of amplitude between PILE(S+S) and GPILE around T 2,S , because the direct contribution echo of the lower surface arrives at the same time as the indirect contribution echoes of the two buried scatterers, and as Q = 1, one must add the contribution of one round trip of the wave between the object and the lower surface for the direct contribution of the object. This last contribution also justifies the difference in the amplitude between GPILE 1/0 and GPILE 1/1 around T 2,S (see Fig. 9). The contribution of the lower surface for Q = 1 arrives at around T 2,S + T OS = 6.81 ns and can be observed in Fig. 8(c). This contribution justifies the difference in amplitude between GPILE 1/0 and GPILE 1/1 between [6.5-7.5] ns as shown in Fig. 9. We now want to analyze the different contributions that make up the third echo of GPILE. To do this, remember that the second echo is made up of four contributions, namely, a direct contribution and another indirect for each of the buried scatterers. Each of the four contributions of the second echo will create four contributions for the third echo, making a total of 16 contributions for the third echo, for a fixed value of Q. We are first interested in the cases of third echoes which do not contain secondary interactions, i.e., for which Q = 0. We have seen that for the second echo, the direct contribution of the object arrives before its indirect contribution, which arrives at the same time as the two contributions of the lower surface. Thus, we can group the 16 contribution echoes that make up the third echo of GPILE into three classes. In the first class, we have the echo which involves only the direct contribution of the object, for both the second echo and third echo. This echo arises at the time T SO + T 2,O = 4.64 ns where T SO is the time of a round trip of the wave between the upper surface and the object, which is why PILE(S+O) coincides with GPILE between [4][5][6] ns as shown in Fig. 8(d). In the second class, we have all the other echoes that involve the direct contribution of the object, whatever the contribution of the lower surface. These echoes are six, arise at the time T SO + T 2,S = 6.15 ns, and can be observed in the interval [6][7] ns. In the third class, we have all the echoes which do not involve the direct contribution of the object, whatever the contribution of the lower surface. These echoes are nine and arrive at the time T SS + T 2,S = T 3,S = 7.65 ns where T SS is the time of a round trip of the wave between the upper surface and the lower surface, and can be observed in the interval [7-8.5] ns.
Note that for Q = 1, it suffices to add the time of a round trip of the wave between the object and the lower surface, and recalculate all the contribution echo times. However, we cannot detail these calculations here as long as these contributions are very weak as we have already underlined.

CONCLUSION
In this paper, the GPILE method has been developed to rigorously calculate the frequency and time responses of the scattered field for an object buried between two random rough surfaces. This method has a straightforward physical interpretation and allows us to distinguish the primary echo of the upper surface, the multiple echoes coming from the lower surface, and those arising from the object. The GPILE method is applied to simulate the Ground Penetrating Radar signal at nadir in both frequency and time domains. The simulated signal was used to detect the presence of the object buried between the two random rough surfaces.
The next step of this research program is to adapt the physical and geometric properties of the object to that of a crack buried under the road surface, in order to apply the results to the detection of non-emerging cracks under the base of a road.

ACKNOWLEDGMENT
The authors acknowledge the support of the French Agence Nationale de la Recherche (ANR) through the project ACIMP (Improvement of knowledge of complex media composed of cracks and multiparameter inversion) under reference ANR-18-CE22-0020.

APPENDIX A. SUB-MATRICES OF THE IMPEDANCE MATRIX
In the case of electromagnetic wave scattering by an object buried between two rough surfaces, the sub-matrices of the impedance matrix are: ) ,

00
) , where the matricesZ ii are the self-impedance matrices of the surfaces S i (i = 1, 2, 3), while the matrices Z ij are the coupling matrices between interfaces S i and S j (i ̸ = j).
Since an interface cannot be illuminated from both above and below (from the outside and the inside for an object) by the same scatterer, this mathematically results in the fact that the upper blocks of matricesZ 21 andZ 31 are zero because the upper surface is illuminated only by the object and the lower surface in Ω 1 . Likewise, the lower blocks ofZ 12 andZ 32 are zero because the object is illuminated only from the outside. Finally, in the matricesZ 13 andZ 23 , the lower blocks are zero because the lower surface is illuminated only in Ω 1 .
The elements of the sub-matrices of the impedance matrix are given bȳ , for m = n.
, for m = n.
, for m = n, , v =n ·ẑ is the direction of the normal to the surface, |∆ n | the discretization step size, H (1) 0 a zero-order Hankel function of the first kind, and H (1) 1 its derivative. The elements of the matricesĀ 2 ,B 2 ,C 2 andD 2 are obtained from those ofĀ 1 ,B 1 ,C 1 andD 1 , in which the wave numbers {k 0 , k 0 , k 1 , k 1 } are replaced by {k 1 , k 1 , k 2 , k 2 } respectively. In the same way, we obtain the elements ofĀ 3 ,B 3 ,C 3 andD 3 by replacing the wave numbers by {k 1 , k 1 , k 3 , k 3 }.
Moreover, the elements of the coupling matrices among the three interfaces are given ∀i, i ′ = 1, 2, 3, by:Ā The indexes i, n and i ′ , m indicate that the points considered are on two different surfaces.

APPENDIX B. RESOLUTIONS OF THE SYSTEMZX = B BY THE SUBSTITUTION METHOD
Consider the linear system obtained by the MoM for the electromagnetic scattering by an object buried between two surfaces (or two objects buried under a surface):     Z 11 X 1 +Z 21 X 2 +Z 31 X 3 = b 1 Z 12 X 1 +Z 22 X 2 +Z 32 X 3 = 0 Z 13 X 1 +Z 23 X 2 +Z 33 X 3 = 0.

(B1)
We have two possibilities to solve the system.
(i) The expression of X 3 is introduced in the second equation, which makes it possible to write X 2 as a function of X 1 so that Then, the expressions of X 2 are introduced into the third equation, which makes it possible to write X 3 as a function of X 1 so that The expressions of X 2 and X 3 are introduced into the first equation and rewritten as and we get where the characteristic matrices are given by Finally, Note that the obtained results are exactly the same as those of the GPILE method presented in this paper. (ii) We first take X 2 in the second equation, which we inject into the third equation. This allows us to express X 3 as a function of X 1 . The expression of X 3 is then introduced in the second equation and gives explicitly X 2 as a function of X 1 . The expressions of X 2 and X 3 introduced in the first equation give the relation (B10) where the characteristic matrices become Then, we deduce the surface currents on the object and the lower interface and The obtained results are also equivalent to the GPILE method presented in this paper.
To have expressions which show the symmetry of the problem, we choose to take the expressions (B3) and (B15) for the surface currents on the object and on the lower surface (therefore (B8) and (B12) for the characteristic matrices). Thus, we obtain expressions which have a physical interpretation conforming to the PILE, EPILE, and GPILE methods.