Marc Songolo
email: marc.songolo@gmail.com

Nicolas Pinel

Christophe Bourlier

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

published or not. The documents may come

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) [START_REF] Jin | The Finite Element Method in Electromagnetics[END_REF][START_REF] Volakis | Finite Element Method for Electromagnetics[END_REF] or Finite Difference Time Domain (FDTD) method [START_REF] Yee | Numerical solution of initial boundary value problems involving Maxwell equations in isotropic media[END_REF][START_REF] Cole | Computing the Flow of Light: Nonstandard FDTD Methodologies for Photonics Design[END_REF] and the integral boundary methods, where the discretization uses the Method of Moments (MoM) [START_REF] Harrington | Field Computation by Moment Methods[END_REF][START_REF] Bourlier | Method of Moments for 2-D Scattering Problems[END_REF]. 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) [START_REF] Kapp | A new numerical method for rough-surface scattering calculations[END_REF]. The major drawback of this method is its slow convergence for the dielectric case. The Forward-Backward (FB) method [START_REF] Holliday | Forward-backward: A new method for computing low-grazing angle scattering[END_REF] adapts better to dielectric interfaces [START_REF] Iodice | Forward-backward method for scattering from dielectric rough surfaces[END_REF]. An accelerated version of the FB uses the spectral representation of Green's function (FB-SA) [START_REF] Chou | A novel acceleration algorithm for the computation of scattering from rough surfaces with the forward-backward method[END_REF]. The Banded-Matrix-Iterative-Approach (BMIA) method [START_REF] Tsang | A banded matrix iterative approach to Monte Carlo simulations of scattering of waves by large scale random rough surface problems: TM case[END_REF][START_REF] Tsang | A banded matrix iterative approach to Monte Carlo simulations of large scale random rough surface scattering: TE case[END_REF] 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) [START_REF] Tsang | Monte-Carlo simulations of large-scale problems of random rough surface scattering and applications to grazing incidence with the BMIA/canonical grid method[END_REF].

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) [START_REF] Kuo | Scattering from multilayer rough surfaces based on the extended boundary condition method and truncated singular value decomposition[END_REF], Forward-Backward with Spectral Acceleration (FBSA) method [START_REF] Moss | Forward backward method with spectral acceleration for scattering from layered rough surfaces[END_REF], and Steepest Descent Fast Multipode Method (SDFMM) [START_REF] El-Shenawee | Polarimetric scattering from two-layered two-dimensional random rough surfaces with and without buried objects[END_REF]. However, all these methods have some constraints. In contrast, the Propagation Inside-Layer Expansion (PILE) method [START_REF] Déchamps | Fast numerical method for electromagnetic scattering by rough layered interfaces: Propagation-inside-layer expansion method[END_REF] 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 [START_REF] Déchamps | Electromagnetic scattering from a rough layer: Propagationinside-layer expansion method combined to an updated BMIA/CAG approach[END_REF]. 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 [START_REF] Déchamps | Electromagnetic scattering from a rough layer: Propagationinside-layer expansion method combined to the forward-backward novel spectral acceleration[END_REF]. The PILE method has been adapted to scattering from an object buried under a rough surface [START_REF] Bourlier | A fast method to compute scattering by a buried object under a randomly rough surface: PILE combined to FB-SA[END_REF]. 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 [START_REF] Kubické | Scattering by an object above a randomly rough surface from a fast numerical method: Extended PILE method combined to FB-SA[END_REF]. A generalization of the PILE method (GPILE: Generalized Propagation Inside-Layer Expansion) has been presented in [START_REF] Bourlier | Generalization of PILE method to the EM scattering from stratified subsurface with rough interlayers: Application to the detection of debondings within pavement structure[END_REF], 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 [START_REF] Bourlier | Generalization of PILE method to the EM scattering from stratified subsurface with rough interlayers: Application to the detection of debondings within pavement structure[END_REF] 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 [START_REF] Déchamps | Fast numerical method for electromagnetic scattering by rough layered interfaces: Propagation-inside-layer expansion method[END_REF] 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.

GEOMETRY OF THE PROBLEM

Consider in 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.

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 [START_REF] Thorsos | The validity of the Kirchhoff approximation for rough surface scattering using a Gaussian roughness spectrum[END_REF], with an attenuation parameter g which controls the extent of the incident wave.

METHOD OF MOMENTS

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' principle. From the Method of Moments, the boundary integral equations are discretized on each surface of the scatterer, leading to the linear system ZX = b, where the impedance matrix is

x 0 z sca, 1 sca, 2 sca, 3 sca, 4 à à à à inc inc sca à µ µ : 1 r 1 ( )
: 3 r 3 ( ) : 0 r 0 ( ) S 1 S 3 " " " n3 n1 z = H 13 - n2 : 2 r 2 ( ) " S 2
Z =   Z11 Z21 Z31 Z12 Z22 Z32 Z13 Z23 Z33   . ( 1
)
The impedance matrix is of size 2(

∑ 3 i=1 N i ) × 2( ∑ 3 i=1 N i ),
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 [START_REF] Bourlier | Generalization of PILE method to the EM scattering from stratified subsurface with rough interlayers: Application to the detection of debondings within pavement structure[END_REF]. The unknown vectors X i containing the surface currents and their normal derivatives are written as

X i = [ ψ i (r 1 ) . . . ψ i (r N i ) ∂ψ i (r 1 ) ∂n i . . . ∂ψ i (r N i ) ∂n i ] T (2)
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:

b = [b 1 b 2 b 3 ] T . ( 3
)
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 [START_REF] Déchamps | Fast numerical method for electromagnetic scattering by rough layered interfaces: Propagation-inside-layer expansion method[END_REF], an object below a rough surface [START_REF] Bourlier | A fast method to compute scattering by a buried object under a randomly rough surface: PILE combined to FB-SA[END_REF], and an object above a rough surface [START_REF] Kubické | Scattering by an object above a randomly rough surface from a fast numerical method: Extended PILE method combined to FB-SA[END_REF].

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 .

From the expressions of the impedance matrices for the dielectric cases given in Appendix A, we obtain:

                     TE: Z22 = B2 , Z12 = [ Ā12 ρ 10 B12 ] , Z32 = [ Ā32 B32 ] , Z21 = [ 0 B21 ] T , Z23 = [ B23 0] T . TM: Z22 = Ā2 , Z12 = [ Ā12 ρ 10 B12 ] , Z32 = [ Ā32 B32 ] , Z21 = [ 0 Ā21 ] T , Z23 = [ Ā23 0] T , (4)
where ρ 10 = ϵ r1 /ϵ r0 with ϵ r0 and ϵ r1 the relative permitivities of Ω 0 and Ω 1 , respectively. On the other hand, if the lower surface is a perfect conductor, then

                     TE : Z33 = B3 , Z13 = [ Ā13 ρ 10 B13 ] , Z23 = [ Ā23 B23 ] , Z31 = [ 0 B31 ] T , Z32 = [ B32 0] T . TM : Z33 = Ā3 , Z13 = [ Ā13 ρ 10 B13 ] , Z23 = [ Ā23 B23 ] , Z31 = [ 0 Ā31 ] T , Z32 = [ Ā32 0] T .
(

) 5
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 [START_REF] Déchamps | Fast numerical method for electromagnetic scattering by rough layered interfaces: Propagation-inside-layer expansion method[END_REF] 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:

P11 = Z11 , P22 = ( Z22 Z32 Z23 Z33 ) , P12 = ( Z12 Z13 ) T , P21 = ( Z21 Z31 ). (6)
The PILE method in the case of two scatterers consists in inverting the impedance matrix

Z = ( P11 P21 P12 P22 ) (7)
by decomposition of domains from the Taylor series expansion of the inverse of the Schur complement [START_REF] Press | Numerical Recipes[END_REF]:

Z-1 = ( T Ū V W) (8)
where the matrices T, Ū, V, and W are expressed as functions of four block matrices of Z as follows

             T = ( P11 -P21 P-1 22 P12 ) -1 Ū = - ( P11 -P21 P-1 22 P12 ) -1 P21 P-1 22 V = -P-1 22 P12 ( P11 -P21 P-1 22 P12
) -1

W = P-1 22 + P-1 22 P12 ( P11 -P21 P-1 22 P12 ) -1 P21 P-1 22 . ( 9
)
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:

         Y 1 =   P ∑ p=0 Np c   P-1 11 b 1 Y 2 = -P-1 22 P12 Y 1 (10)
where Nc = P-1 11 P21 P-1 22 P12 ( 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 of Nc . Denoting ∥ Nc ∥ sr , this norm is always strictly smaller than 1 as shown in [START_REF] Déchamps | Fast numerical method for electromagnetic scattering by rough layered interfaces: Propagation-inside-layer expansion method[END_REF], 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 matrix Nc to obtain an expression equivalent to X 1 , with the only issue being to invert the block matrix P22 as shown in [START_REF] Bourlier | Method of Moments for 2-D Scattering Problems[END_REF]. By reusing the inverse of the Schur complement algorithm, we obtain

Nc = Z-1 11 Z21 ( T Z12 + Ū Z13 ) + Z-1 11 Z31 ( V Z12 + W Z13 ) (12)
where the matrices T, Ū, V, and W are given by ( 9) as functions of four blocks of the matrix P22 . Recall in addition that, as the matrix Nc 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:

Nc,21 = Z-1 11 Z21 ( T Z12 + Ū Z13 ) (13) and
Nc,31 = Z-1

11 Z31 ( V Z12 + W Z13 ) . ( 14
)
We will see, after developing these two expressions, that the matrices Nc,21 and Nc,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)

Nc,21 ≈ Q ∑ q=0 M(q) c,21 = Z-1 11 Z21   Q ∑ q=0 Mq c,32   Z-1 22 ( Z12 -Z′ 12 ) (15)
and Nc,31

≈ Q ∑ q=0 M(q) c,31 = Z-1 11 Z31 Z-1 33   Z13 -Z23   Q ∑ q=0 Mq c,32   Z-1 22 ( Z12 -Z′ 12 )   (16)
where

Mc,32 = Z-1 22 Z32 Z-1 33 Z23 (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 matrix Z′ 12 = Z32 Z-1 33 Z13 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 [START_REF] El-Shenawee | Polarimetric scattering from two-layered two-dimensional random rough surfaces with and without buried objects[END_REF] obtained by this method is Nc,31

≈ Q ∑ q=0 M(q) c,31 = Z-1 11 Z31   Q ∑ q=0 Mq c,23   Z-1 33 ( Z13 -Z′ 13 ) (18)
where

Mc,23 = Z-1 33 Z23 Z-1 22 Z32 , ( 19
)
and the matrix Z′ 13 = Z23 Z-1 22 Z12 represents the coupling matrix between the upper surface and lower surface via the object.

Expression [START_REF] Déchamps | Electromagnetic scattering from a rough layer: Propagationinside-layer expansion method combined to an updated BMIA/CAG approach[END_REF] 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 matrices Z12 and Z′ 12 which make it possible to quantify these contributions for the object. The coupling matrices Z13 and Z′ 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:

X 1 = P ∑ p=0   Q ∑ q=0 M(q) c,21 + Q ∑ q=0 M(q) c,31   p Z-1 11 b 1 . ( 20
)
Finally,

X 1 = P ∑ p=0   Q ∑ q=0 M(q) c   p Z-1 11 b 1 . ( 21
)
It is important to note that the orders P and Q are arbitrary. In fact, according to the construction of the method (see expressions [START_REF] Chou | A novel acceleration algorithm for the computation of scattering from rough surfaces with the forward-backward method[END_REF]), 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

X 2 = - ( T Z12 + Ū Z13 ) X 1 = - T ( Z12 -Z′ 12 ) X 1 (22)
from where

X 2 = -   Q ∑ q=0 Mq c,32   Z-1 22 ( Z12 -Z′ 12 ) X 1 (23)
and for the lower surface

X 3 = - ( V Z12 + W Z13 ) X 1 = - V ( Z12 -Z′ 12 ) X 1 -Z-1 33 Z13 X 1 (24)
from where

X 3 = Z-1 33 Z23   Q ∑ q=0 Mq c,32   Z-1 22 ( Z12 -Z′ 12 ) X 1 -Z-1 33 Z13 X 1 . ( 25
)
The calculations of the currents on the lower surface by applying the substitution method (see Appendix B) make it possible to obtain

X 3 = -   Q ∑ q=0 Mq c,23   Z-1 33 ( Z13 -Z′ 13 ) X 1 . ( 26
)
The expressions of the currents on the object [START_REF] Thorsos | The validity of the Kirchhoff approximation for rough surface scattering using a Gaussian roughness spectrum[END_REF] and on the lower surface (26) have a physical interpretation in accordance with the EPILE method [START_REF] Kubické | Scattering by an object above a randomly rough surface from a fast numerical method: Extended PILE method combined to FB-SA[END_REF]. 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 [START_REF] Kubické | Scattering by an object above a randomly rough surface from a fast numerical method: Extended PILE method combined to FB-SA[END_REF]. By considering that the object and the lower surface are illuminated by the fields -Z12 X 1 and -Z13 X 1 respectively, we obtain relations [START_REF] Thorsos | The validity of the Kirchhoff approximation for rough surface scattering using a Gaussian roughness spectrum[END_REF] 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 [START_REF] Chou | A novel acceleration algorithm for the computation of scattering from rough surfaces with the forward-backward method[END_REF] and [START_REF] Kubické | Scattering by an object above a randomly rough surface from a fast numerical method: Extended PILE method combined to FB-SA[END_REF] is that in [START_REF] Chou | A novel acceleration algorithm for the computation of scattering from rough surfaces with the forward-backward method[END_REF], Q = ∞, because the impedance matrix P22 (for the lower surface and the object) is inverted directly, while in [START_REF] Kubické | Scattering by an object above a randomly rough surface from a fast numerical method: Extended PILE method combined to FB-SA[END_REF],

Q = 0, 1, 2, etc. In this way Nc = ∑ Q q=0 M(q) 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 in Z13 = 0 (and therefore Z31 = 0). It would follow that Nc,31 = 0, and one would find a GPILE method equivalent to that developed in [START_REF] Bourlier | Generalization of PILE method to the EM scattering from stratified subsurface with rough interlayers: Application to the detection of debondings within pavement structure[END_REF] for the case of three superimposed interfaces.

For the numerical simulations, we implement expression [START_REF] Kubické | Scattering by an object above a randomly rough surface from a fast numerical method: Extended PILE method combined to FB-SA[END_REF] 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 [START_REF] Moss | Forward backward method with spectral acceleration for scattering from layered rough surfaces[END_REF] and [START_REF] Déchamps | Electromagnetic scattering from a rough layer: Propagationinside-layer expansion method combined to an updated BMIA/CAG approach[END_REF]. 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 norm

RRE = norm(X -X LU ) norm(X LU ) , ( 27
)
(X) = N ∑ i=1 |x i | 2 . ( 28
)
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 via Z′ 12 and Z′ 13 as shown by relations ( 15) and ( 18), respectively.

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 [START_REF] Déchamps | Fast numerical method for electromagnetic scattering by rough layered interfaces: Propagation-inside-layer expansion method[END_REF] 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 [START_REF] Bourlier | Full wave PILE method for the electromagnetic scattering from random rough layer[END_REF]. In the time domain, it is defined as In the Fourier domain, it is defined as

s(t) = ( 2π 2 f 2 c t 2 -1 ) exp ( -π 2 f 2 c t 2 ) . ( 29
)
s(f ) = - 2f 2 f 2 c √ π exp ( - f 2 f 2 c ) , f ≥ 0. ( 30
)
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), the signal is Hermitian, i.e., s(-f ) * = s(+f ), where the symbol * stands for the complex conjugate. The number of frequencies is N f = 117, and the band is [0.2-6.0] GHz. To have a better resolution in the frequency domain, we applied the zero-padding technique to increase N f to 2 12 = 4096. Thus, the sampling step in the frequency domain is ∆f = 0.00142 GHz.

Frequency Responses

Figure 7 plots the frequency response of the modulus of the scattered field versus the frequency f . The simulations parameters are the same as in Fig. 4, except the thickness H = 240 mm and the center of the cylinder C = (0, -0.375H). The observation abscissa and height are x obs = 0 and z obs = 470 mm, respectively. 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:

   T 1 = 2z 0 /c T 2,S = T 1 + 2HRe( √ ϵ r,1 ) T 3,S = T 1 + 4HRe( √ ϵ r,1 ) (31)
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 [START_REF] Moss | Forward backward method with spectral acceleration for scattering from layered rough surfaces[END_REF] that the contribution matrix is composed by the direct contribution, represented by Z12 (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 by Z′ 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 [START_REF] Déchamps | Electromagnetic scattering from a rough layer: Propagationinside-layer expansion method combined to an updated BMIA/CAG approach[END_REF] and represented by the contribution matrices Z13 and Z′ 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 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 [START_REF] Cole | Computing the Flow of Light: Nonstandard FDTD Methodologies for Photonics Design[END_REF][START_REF] Harrington | Field Computation by Moment Methods[END_REF][START_REF] Bourlier | Method of Moments for 2-D Scattering Problems[END_REF] 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 [START_REF] Bourlier | Method of Moments for 2-D Scattering Problems[END_REF][START_REF] Kapp | A new numerical method for rough-surface scattering calculations[END_REF] 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. 1 its derivative. The elements of the matrices Ā2 , B2 , C2 and D2 are obtained from those of Ā1 , B1 , C1 and D1 , 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 , B3 , C3 and D3 by replacing the wave numbers by {k 1 , k 1 , k 3 , k 3 }.

D1,mn = j|∆

n | √ 1 + γ 2 n 4      [ 1 + 2j π ln(0.164k 1 √ 1 + γ 2 n ∆ n ) ] , for m = n, H (1) 0 (k 1 ∥r n -r m ∥), for m ̸ = n. where (x n , z n ) ∈ S i , (x m , z m ) ∈ S i , γ = dz/dx, ∥r n -r m ∥ = √ (x n -z n ) 2 + (x m -z m ) 2 , v = n •
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.

Āii ′ ,nm = - jk 1 v i,

APPENDIX B. RESOLUTIONS OF THE SYSTEM ZX = 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):

     Z11 X 1 + Z21 X 2 + Z31 X 3 = b 1 Z12 X 1 + Z22 X 2 + Z32 X 3 = 0 Z13 X 1 + Z23 X 2 + Z33 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

X 2 = ( Ī -Z-1 22 Z32 Z-1 33 Z23 ) -1 Z-1 22 ( Z′ 12 -Z12 ) X 1 (
X 3 = Z-1 33 Z23   Q ∑ q=0 Mq c,32   Z-1 22 ( Z12 -Z′ 12 ) X 1 -Z-1 33 Z13 X 1 . (B5)
The expressions of X 2 and X 3 are introduced into the first equation and rewritten as Then, we deduce the surface currents on the object and the lower interface

X 1 + Z-
X 2 = Z-1 22 Z32   Q ∑ q=0 Mq c,23   Z-1 33 ( Z13 -Z′ 13 ) X 1 -Z-1 22 Z12 X 1 (B14)
and

X 3 = -   Q ∑ q=0 Mq c,23   Z-1 33 ( Z13 -Z′ 13 ) X 1 . (B15)
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.

Figure 1 .

 1 Figure 1. Electromagnetic scattering from an object buried between two 1-D rough surfaces (2-D problem: plane (x, ẑ)).

Figure 2 .

 2 Figure 2. Surface heights versus their abscissa for Gaussian PDF and ACF, with σ h,1 = 1 mm and σ h,3 = 2.5 mm, L ch,1 = 15 mm and L ch,3 = 30 mm. The center of the cylinder is C = (0, -0.5H), and a = 100 mm and b = 10 mm.

Figure 3 .

 3 Figure 3. Modulus of the scattered field in decibel scale (10 log 10 |ψ sca |) versus the observation abscissa x obs ∈ [-L/2; L/2], for an observation height z obs = 470 mm. The frequency is f = 2 GHz and P = 0.

Figure 4 .

 4 Figure 4. Modulus of the scattered field in decibel scale (10 log 10 |ψ sca |) versus the observation abscissa x obs ∈ [-L/2; L/2] and for an observation height z obs = 470 mm. The frequency is f = 2 GHz and the orders P and Q vary. (a) P = 1 and Q = 0. (b) P = Q = 1. (c) P = 2 and Q = 0. (d) P = 2 and Q = 1.

Figure 5 .

 5 Figure 5. Modulus of the scattered field in decibel scale (10 log 10 |ψ sca |) versus the observation abscissa x obs ∈ [-L/2; L/2] and for an observation height z obs = 470 mm. The frequency is f = 2 GHz and the orders P = 2 and Q = 1. (a) Cylinder with a = 1000 mm and b = 10 mm. (b) Cylinder with a = b = 10 mm. (c) The object is a perfect conductor, with a = 100 mm and b = 10 mm. (d) The lower surface is a perfect conductor, with a = 100 mm and b = 10 mm.

Figure 6 .

 6 Figure 6. (a) The Fourier transform of the signal versus the frequency f (in GHz). (b) Signal versus the time t (in ns).

Figure 7 .

 7 Figure 7. Frequency responses of the modulus of the scattered field. x obs = 0, z obs = 470 mm. (a) Input and Output signals. (b) First echo (P = 0). (c) First two echoes (P = 1). (d) First three echoes (P = 2).

Figs. 8 Figure 8 .

 88 Figure 8. Time responses of the modulus of the scattered field. x obs = 0, z obs = 470 mm. (a) Input and Output signals. (b) First echoes (P = 0). (c) Second echoes (P = 1). (d) Third echoes (P = 2).

Figure 9 .

 9 Figure 9. Time response of the modulus of the scattered field: Direct and Indirect contributions. x obs = 0 and z obs = 470 mm.

  ẑ 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 11 Z21

 11 X 2 + Z-1 11 Z31 X 3 = Z-1 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 Q

  n |∆ i,n | 4 (k 1 ∥r i,n -r i ′ ,m ∥) ∥r i,n -r i ′ ,m ∥ [ γ i,n (x i,n -x i ′ ,m ) -(z i,n -z i ′ ,m ) ]

	(1) 1 and H Bii ′ ,nm = j|∆ i,n | √ 4	1 + γ 2 i,n	H	(1)

0 (k 1 ∥r i,n -r i ′ ,m ∥).

Progress In Electromagnetics Research B, Vol. 96, 2022

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:

where the matrices Zii are the self-impedance matrices of the surfaces S i (i = 1, 2, 3), while the matrices Zij 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 matrices Z21 and Z31 are zero because the upper surface is illuminated only by the object and the lower surface in Ω 1 . Likewise, the lower blocks of Z12 and Z32 are zero because the object is illuminated only from the outside. Finally, in the matrices Z13 and Z23 , 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 by Ā1,mn =

, for m = n.


Rough layer scattering filled by elliptical cylinders from the method of moments combined with the characteristic basis function method and the Kirchoff approximation

In this paper, the electromagnetic field scattered by several 2D scatterers of any shape is calculated rigorously from the boundary integral equations discretized by the method of moments with the point matching method and pulse basis functions. In addition, the resulting linear system is efficiently solved from the domain decomposition method named the characteristic basis function method. To accelerate the computation of the primary basis functions, which requires solving sublinear systems, the Kirchoff approximation is applied for metallic and dielectric objects. The efficiency of the method is tested on several applications met in practice: stack of rough interfaces separating homogeneous media, collection of metallic and dielectric elliptical cylinders, collection of coated elliptical cylinders, and a combination of the previous scenarios.

INTRODUCTION

The study of the wave scattering from several scatterers of any shape is a subject of great interest. The applications of such research concern many areas such as remote sensing, radar surveillance, optics, and ocean acoustics.

For a collection of objects of canonical shape, like circular cylinders (of infinite length) and spheres, the scattered field can be computed analytically by introducing special functions, e.g., Bessel's and spherical Bessel's [START_REF] Frezza | Introduction to electromag-698 netic scattering: tutorial[END_REF]2]. For elliptical cylinders [3-5], Mathieu's functions [START_REF] Abramowitz | Handbook of Mathematical 713 Functions[END_REF] are introduced, but they are difficult to program, unlike well-known Bessel's functions. Adding a boundary, like a smooth plate of infinite area (space divided into two media), the previous formulations can be extended [START_REF] Sommerfeld | Lectures on Theoretical Physics[END_REF][START_REF] Youssif | Scattering by two penetrable cylinders 717 at oblique incidence. I. The analytical solution[END_REF][START_REF] Valle | Electromagnetic wave scat-720 tering from conducting cylindrical structures on flat substrates: study tion[END_REF][10][11][12][13][14][15][16][17], and the difficulty of programming increases. For a stack of rough interfaces separating homogeneous media, asymptotic approaches, in which simplifying assumptions are introduced, have been developed. For small roughness, we can cite the small perturbation method [18][START_REF] Ohlidal | Scattering of light from multilayer systems 750 with rough boundaries[END_REF][START_REF] Fuks | Wave diffraction by rough interfaces in an arbitrary plane-layered medium[END_REF][START_REF] Soubret | Backscattering enhancement of an electromagnetic wave scattered by two dimensional[END_REF][START_REF] Tabatabaeenejad | Bistatic scattering from dielectric structures with two rough boundaries using the small perturbation method[END_REF][START_REF] Sanamzadeh | Scattering of electromagnetic waves from 3D multilayer random rough surfaces based on the second-order small perturbation method: energy conservation, reflectivity, and emissivity[END_REF][START_REF] Dusséaux | Level crossing rate and average fade distance of signal scattering from rough layered interfaces[END_REF] and, for high roughness, the geometric optics approximation [START_REF] Pinel | Bistatic scattering from one-dimensional random rough homogeneous layers in the highfrequency limit with shadowing effect[END_REF][START_REF] Pinel | Scattering from very rough layers under the geometric optics approximation: further investigation[END_REF][START_REF] Pinel | A geometrical optics model of three dimensional scattering from a rough surface over a planar surface[END_REF]. For a complex scenario, like inhomogeneous (dielectric objects are present in the layer) layered rough interfaces, it is difficult to derive a closed-form expression of the scattered field.

The well-known method of moments (MoM) [START_REF] Harrington | Field Computation by Moment Method[END_REF][START_REF] Tsang | Scattering of Electromagnetic Waves, Numerical Simulations[END_REF][START_REF] Bourlier | Method of Moments for 2D Scattering Problems[END_REF][START_REF] Sun | Field-only surface integral equations: scattering from a perfect electric conductor[END_REF] is a way of rigorously solving this type of scattering problem by converting the boundary integral equations into a linear system, in which the impedance matrix must be inverted to determine the surface currents. However, the direct solution of the linear calculating the primary basis functions (PBFs). Next, the coupling between them is accounted for via the computation of the characteristic matrix, which involves the coupling matrices between the subproblems. To accelerate the calculation of PBFs, which requires one to solve sublinear systems, the Kirchoff approximation (KA) is applied for metallic and dielectric objects. The efficiency of the method is tested on several applications met in practice: stack of rough interfaces separating homogeneous media, collection of metallic and dielectric elliptical cylinders, collection of coated elliptical cylinders, and a combination of the previous scenarios.

The paper is organized as follows. Section 2 extends the MoM to P scatterers of any shape. Section 3 briefly summarizes the CBFM. Section 4 addresses the KA approximation for the derivation of the PBFs. Section 5 presents numerical results on four scenarios. The last section gives concluding remarks.

METHOD OF MOMENTS FOR SEVERAL SCATTERERS

This section presents the MoM for solving the electromagnetic wave scattering from several scatterers. For one and two scatterers, this approach is thoroughly explained in the textbook [START_REF] Bourlier | Method of Moments for 2D Scattering Problems[END_REF].

To sum, the boundary integral equations are applied on each scatterer and are discretized from the MoM by using the pointmatching method and the pulse basis function. This leads to the linear system ZX = b, in which Z is the impedance matrix of the two scatterers, b is a vector related to the incident field on the scatterers, and X the surface currents on the scatterers, i.e., the unknowns of the problem. In this section, this approach is generalized to several scatterers.

A. Case of a Single Illuminated Scatterer

For a single scatterer as shown in Fig. 1, the impedance matrix is expressed from four submatrices as

Z11 = Ā11 B11 C11 D11 . (1)
The matrix Ā11 is the matrix obtained from the Neumann boundary condition, i.e., when the scatterer is assumed to be perfectly conducting (metallic), and the transverse electric medium 1 . The elements of these four submatrices are given in

b 1 = b 1 0 =    ψ inc (r 1 ) ... ψ inc (r N 1 ) 0 ... 0 N 1 times b T 1 ,r∈S 1    T , (2)
where ψ inc is the incident wave illuminating the scatterer. The 118 symbol T stands for the transpose operator, and N 1 is the num-119 ber of discretization points on the surface S 1 of the object. This 120 means that the size of the matrix is 2N 1 × 2N 1 .

121

The unknown vector X 1 of length 2N 1 is expressed as 122

X 1 = ψ 1 (r 1 ) ... ψ 1 (r N1 ) ∂ψ1(r1) ∂n ... ∂ψ1(r N 1 ) ∂n T r p∈[1;N1] ∈ S 1 , (3)
where 128

P(r, r 1 ) = Āij Bij r j =r 1 ,r i =r , (5)
and r 1 ∈ S 1 (r = r 1 ). If r ∈ 1 (r = r 1 ); then, Eq. ( 4) is In the far field (kr 1, where k is the wavenumber of the 136 medium), the propagation matrix can be simplified as [START_REF] Bourlier | Method of Moments for 2D Scattering Problems[END_REF] 137

P∞ (r, r 1 ) = j 4 2 πkr e -j (π/4+kr ) j v 1 k sca • n1 e -j k sca •r 1 1 + γ 2 1 1 e -j k sca •r 1 1 + γ 2 1 1 , ( 6
)
where k sca = k( x sin θ sca + ẑ cos θ sca ) (see Fig. 1) stands for the The radar cross section (RCS) is written as [START_REF] Bourlier | Method of Moments for 2D Scattering Problems[END_REF] 142

RCS = lim r →∞ 2πr ψsca (r) 2 = 1 4k j v 1 k sca • n1 e -j k sca •r 1 1 + γ 2 1 1 e -j k sca •r 1 1 + γ 2 1 1 X 1 2 , ( 7
)
where RCS is a vector of the same length as ψsca (r).

The normalized radar cross section (NRCS, dimensionless) is 144 written as [START_REF] Tsang | Scattering of Electromagnetic Waves, Numerical Simulations[END_REF][START_REF] Bourlier | Method of Moments for 2D Scattering Problems[END_REF] 145

NRCS = lim r →∞ r 2η 0 | ψsca (r)| 2 p inc , (8)
where η 0 is the wave impedance of medium 0 and p inc the inci-146 dent power. For an incident plane wave of unitary amplitude, 147 p inc = L cos θ inc /(2η 0 ), so that NRCS = RCS cos θ inc /(2π L). For the case of two illuminated scatterers (Fig. 2, with P = 2), 150 the impedance matrix Z is expressed as [START_REF] Bourlier | Method of Moments for 2D Scattering Problems[END_REF] 151

Z = Z11 Z(1) 12 Z(1) 21 Z22 , (9)
where Z11 [Eq. ( 1)] and Z22 are the self-impedance matrices of The size of Z is reduced to

= Ā12 B12 0N 1 ×N 2 0N 1 ×N 2 , Z (1) 21
= Ā21 B21 0N 2 ×N 1 0N 2 ×N 1 , (
(N 1 + N 2 ) × (N 1 + N 2 ).
171

The excitation vector b of length 2(N 1 + N 2 ) is written as

172 b = b 1 b 2 T , (11)
Fig. 2. Scattering from P illuminated scatterers.

where the vector b i is defined from Eq. (2).

The unknown vector X of length 2(N 1 + N 2 ) is expressed as

X = X 1 X 2 T , (12)
where X i is expressed from Eq. (3). The unknown

X on S 1 ∪ S 2 is computed from X = Z-1 b.
The scattered field ψsca (r) in the medium r ∈ 0 (r / ∈ S i ) is then obtained from the Huygens principle expressed as ψsca (r) = -P(r, r 1 ) P(r, r 2 ) X .

(13)

If r ∈ 1 (r / ∈ S 1 ), then the above equation is applied by taking the plus sign on the matrix P(r, r 1 ) (instead of minus).

If r ∈ 2 (r / ∈ S 2 ), then the above equation is applied by taking the plus sign on the matrix P(r, r 2 ).

In this paper, the formulation is generalized to P scatterers.

The impedance matrix is then expressed as

Z =       Z11 Z (
P 1 Z(1) P 2 ... ZPP       , ( 14
)
where Zii is the self-impedance matrix of the scatterer i, whereas Z(1) ij is the coupling matrix between the objects i and j (propagation of the scattered field from j to i) expressed from Eq. (10).

The size of the matrix Z is N × N where N = 2 i=P i=1 N i .

The excitation vector b of length N is written as

b = b 1 b 2 ... b P T , ( 15
)
where the vector b i is defined by Eq. (2). The unknown vector X

of length N is X = X 1 X 2 ... X P T , ( 16
)
where X i is expressed from Eq. (3) and computed by solving the linear system ZX = b.

The scattered field ψsca (r) in the medium r ∈ 0 (r / ∈ S i ) is then obtained from the Huygens principle expressed as

ψsca (r) = -      P(r, r 1 ) P(r, r 2 ) . . . P(r, r P )      X . ( 17
)
For r ∈ i (i = 0 and r / ∈ S i ), the above equation is applied by taking the plus sign on the matrix P(r, r i ) (instead of minus).

C. Case of P Illuminated Scatterers Where Only One Is Illuminated

For the case of two illuminated scatterers where only one is illuminated (Fig. 3 with P = 2), the impedance matrix Z is expressed as [START_REF] Bourlier | Method of Moments for 2D Scattering Problems[END_REF] Fig. 3. P illuminated scatterers where only one is illuminated (the scatterer 1).

Z = Z11 Z(2) 12 Z(2) 21 Z22 , ( 18
)
where the coupling matrices Z(2) 12 and Z(2) 21 are written as 204

Z(2) 12 = 0N 1 ×N 2 0N 1 ×N 2 Ā12 B12 , Z(2) 21 = Ā21 D21 0N 2 ×N 1 0N 2 ×N 1 , ( 19
)
where 0 is a null matrix. The size of the matrix Z is

2(N 1 + N 2 ) × 2(N 1 + N 2 ).
(2N 1 + N 2 ) × (2N 1 + N 2 ).

213

The excitation vector b of length 2(N 1 + N 2 ) is written as

214 b = b 1 0 1×2N 2 T , (20)
where the vector b 1 is defined by Eq. (2). In comparison with scatterer 2 is not illuminated.

217

The unknown vector X of length 2(

N 1 + N 2 ) is expressed as 218 X = X 1 X 2 T , (21)
where X i is expressed from Eq. (3). The unknown

X on S 1 ∪ S 2 219 is computed from X = Z-1 b. 220
The scattered field ψsca (r) in the medium r ∈ 0 (r / ∈ S 1 ) is 221 then obtained from the Huygens principle expressed as

222 ψsca (r) = -P(r, r 1 )X 1 . ( 22
) If r ∈ 1 (r / ∈ S 1 ), then the scattered field ψsca (r) is 223 ψsca (r) = + P(r, r 1 ) -P(r, r 2 ) X . ( 23
) If r ∈ 2 (r / ∈ S 2 ), then the scattered field ψsca (r) is 224 ψsca (r) = P(r, r 2 )X 2 . ( 24
)
In this paper, the formulation is generalized to P scatterers.

The impedance matrix is then expressed as

Z =            Z11 Z(2) 12 0 0 ... 0 Z(2) 21 Z22 Z(2) 23 0 ... 0 0 Z(2) 32 Z33 Z(2)
P ,P -1 ZPP            , (25)
where Zii (or Zi,i ) is the self-impedance matrix of the scatterer i,

whereas Z(2) ij (or Z(2) i, j
) is the coupling matrix between the objects i and j (propagation of the scattered field from j to i) expressed from Eq. ( 19). Unlike the matrix expressed in Eq. ( 18), Eq. ( 25)

shows that only two adjacent scatterers i and min(|i + 1|, P ) interact, which explains why null matrices appear.

The excitation vector b is given by

b = b 1 0 1×2 i=P i=2 N i T . ( 26
)
It differs from Eq. ( 15) because b i = 0 for i ∈ [2; P ] (these scatterers are not illuminated).

The unknown vector X of length N = 2 i=P i=1 N i is expressed from Eq. ( 16) and computed by solving the linear system ZX = b.

The scattered field ψsca (r) in the medium r ∈ 0 (r / ∈ S 1 ) is obtained from the Huygens principle given by Eq. [START_REF] Dusséaux | Level crossing rate and average fade distance of signal scattering from rough layered interfaces[END_REF]. For r ∈

i (i ∈ [1; P -1], r / ∈ (S i ∪ S i+1
)), the scattered field is

ψsca (r) = + P(r, r i ) -P(r, r i+1 ) X i X i+1 . ( 27
)
For r ∈ P (r / ∈ S P ), the scattered field is ψsca (r) = P(r, r P )X P .

D. Combination of Cases

As an example, in this paragraph, the impedance matrix of the scenario presented in Fig. 4 is determined. It is expressed as

Z =         Z11 Z(1) 12 Z(1) 13 0 0 Z(1) 21 Z22 Z(1) 23 0 0 Z(1) 31 Z(1) 32 Z33 Z(2) 34 Z(2) 35 0 0 Z(2) 43 Z44 Z(1) 45 0 0 Z(2) 53 Z(1) 54 Z55         . ( 28
)
The impedance matrix of scatterers 1, 2, and 3 is expressed from Eq. ( 14), with P = 3. This explains the first three rows and columns of Z with the superscript (1) (Case 1). The impedance matrix of scatterers 4 and 5 is expressed from Eq. ( 9), in which the subscripts {1, 2} are replaced by {4, 5}, and the incident medium is 3 . This explains the last two rows and columns of Z with the superscript (1) (Case 1). The impedance matrix of scatterers 3 and 4 union 5 is expressed from Eq. ( 18). This explains the rows (4,5) and columns (4,5) of Z with the superscript (2) (Case 2). The eight null matrices appear because scatterers 1 and 257 2 do not directly interact with objects 4 and 5.

258

The excitation vector b is given by

259 b = b 1 b 2 b 3 0 1×2(N 4 +N 5 ) T , ( 29
)
for which the scatterers 4 and 5 are not illuminated. It is impor-

    X 1 X 2 X 3   . (30)
For r ∈ 3 (r / ∈ (S 3 ∪ S 4 ∪ S 5 )), the scattered field is 267 ψsca (r) =   + P(r, r 1 ) -P(r, r 4 ) -P(r, r 5 )

    X 1 X 4 X 5   . ( 31
)
For r ∈ 4 (r / ∈ S 4 ), ψsca (r) = P(r, r 4 )X 4 . For r ∈ 5 (r / ∈ 268 S 5 ), ψsca (r) = P(r, r 5 )X 5 .

269

In conclusion, the two cases presented in the previous sections

M =      1 1 1 0 0 1 1 1 0 0 1 1 1 1 1 0 0 1 1 1 0 0 1 1 1      . ( 32
)
If M ij = 0, then the corresponding matrix Zij = 0. The 278 matrix M is symmetric, owing to the reciprocity principle, and 279 becomes

M ii = 1 ∀ i ∈ [1; P ].
M =      1 1 1 0 0 1 1 1 0 0 1 1 1 2 2 0 0 2 1 1 0 0 2 1 1      . ( 33
)
In conclusion, if This is computationally expensive for an electrically large multiscale object or a collection of dielectric objects (i.e., N huge).

M ij = 0 then Zij = 0, Zij = Z(m ij ) ij (i = j )
To tackle this issue, the domain decomposition method, named the characteristic basis function method (CBFM), is applied. A summary of this method can be found in [START_REF] Prakash | Characteristic basis function method: a new technique for efficient solution of method of moments matrix equations[END_REF][START_REF] Bourlier | Two domain decomposition methods, SDIM and CBFM, for the scattering from a two-dimensional perfectly-conducting rough surface: comparison and parametric study[END_REF].

It is important to note that, originally, the CBFM was developed for radiation and 3D scattering problems [START_REF] Prakash | Characteristic basis function method: a new technique for efficient solution of method of moments matrix equations[END_REF] and on a single geometry. Next, it was extended to a single and two 3D dielectric scatterers [START_REF] Laviada | Generation of excitationindependent characteristic basis functions for three-dimensional momogeneous dielectric bodies[END_REF][START_REF] Li | Characteristic basis function method for fast analysis of 3-D scattering from objects buried under rough surfaces[END_REF].

For a 2D problem, some minor changes are needed.

A. Case of Metallic Scatterers

First, the scatterers are assumed to be metallic. The CBFM begins by dividing the geometry of the problem to be analyzed into B blocks, where Zi,i is the self-impedance matrix of the block i. Next, a primary basis function (PBF), Y i ,k IPW , is computed for each block by solving the linear system

Zi ,i Y i ,k IPW = B i ,k IPW , (34)
where the subscript prime indicates that the block i is enlarged and k IPW stands for the k IPW th plane wave (ranging from 1 to The choice of N IPW,i must be relevant to avoid the size of the matrix being too big. The redundant information due to the overestimation of N IPW,i is eliminated via the use of a truncated singular value decomposition (SVD). This means that, from a given threshold CBFM,SVD , the values for which the modulii of the normalized eigenvalues are smaller than CBFM,SVD are removed. The size of J i becomes N i × N IPW,SVD,i , with

N IPW,i
N IPW,SVD,i < N IPW,i .
The last stage of CBFM consists in solving a reduced linear

system ZR a R = b R defined as       ZR 1,1 ZR 1,2 ... ZR 1,B ZR 2,1 ZR 2,2 ... ZR 2,B . . . . . . . . . . . . ZR B,1 ZR B,2 ... ZR B,B            a R 1 a R 2 . . . a R B      =      b R 1 b R 2 . . . b R B      , ( 35
)
where the submatrix ZR i, j and the subvector b R i are defined as

ZR i, j = J H i Zi, j J j N IPW,SVD,i × N IPW,SVD, j b R j = J H j b j N IPW,SVD,i × 1 . (36)
Moreover, the symbol H stands for the conjugate transpose operator, and the indexes i and j range from 1 to B. The unknown vector X i of the block i equals

X i = J i a R i .
The problem is then represented by the characteris- to repeat the procedure, and the surface currents are calculated rapidly. The complexity of CBFM is detailed in [START_REF] Bourlier | Two domain decomposition methods, SDIM and CBFM, for the scattering from a two-dimensional perfectly-conducting rough surface: comparison and parametric study[END_REF].

Originally, the CBFM splits up the single geometry into B blocks, since for a 3D problem the number of unknowns can be huge. In this paper, since the number of unknowns N i of a single scatterer is moderate, it is not required to decompose it into blocks, but it is not a limitation. This means that the number of blocks equals the number of scatterers, B = P , and that the overlapping between the scatterers is not necessary.

B. Case of Dielectric Scatterers

As shown by Eq. ( 1), for a dielectric scatterer, the self-impedance matrix of size 2N i × 2N i is expressed from four submatrices of size N i × N i . Thus, the PFBs are calculated by repeating the procedure four times used for a metallic object and illustrated by Eq. [START_REF] Ye | Fast iterative approach to difference electromagnetic scattering from the target above a rough surface[END_REF]. This leads to the matrix

J i = J i,11 J i,12 J i,21 J i,22 , (37)
where the submatrices { J i, pq } ( p = {1, 2}, q = {1, 2}) are obtained from the four submatices of Zi ,i , i.e., { Āi ,i , Bi ,i , Ci ,i , Di ,i }. This representation implies

ψ(r) = [1 + R(θ )] ψ inc (r)I(r) ∂ψ(r) ∂n = [1 -R(θ )] ∂ψ inc (r) ∂n I(r) , ( 38
)
where R is the Fresnel reflection coefficient defined as 394 R = n 2 cos θ-n 1 cos θ t n 2 cos θ+n 1 cos θ t TM polarization

n 1 cos θ-n 2 cos θ t n 1 cos θ+n 2 cos θ t TE polarization . ( 39
)
In addition, For a metallic surface, R = {+1, -1} for the TM and TE polarizations, respectively. This leads to {ψ = 2Iψ inc , ∂ψ/∂n = 0} for TM and to {ψ = 0, ∂ψ/∂n = 2I∂ψ inc /∂n} for TE.

cos θ = -kinc • n, cos θ t = 1 -n 2 1 (1 -cos 2 θ )/n 2 2 ,
Assuming a plane incident wave ψ inc = e j k inc •r , Eq. [START_REF] Déchamps | Electromagnetic scattering from a rough layer: propagation-inside-layer expansion method combined to the forward-backward novel spectral acceleration[END_REF] becomes

ψ(r) = [1 + R(θ )] e j k inc •r I(r) ∂ψ(r) ∂n = [1 -R(θ )] e j k inc •r I(r) j k inc • n . ( 40
)
In Eq. ( 37), the PBFs { J i,11 , J i,12 } are obtained from Eq. ( 40) by taking R = {+1, -1} (metallic case), respectively, whereas { J i,21 , J i,22 } are obtained from Eq. ( 40).

In other words, it is equivalent to write that Zi,i is a diagonal To calculate the PBFs with KA, the resulting complexity is O(N i ) instead of O(N 3 i ) from a conventional LU decomposition. For the numerical results, keep in mind that there is no overlapping between the P = B scatterers, meaning that First, the scattering from a stack of P = 6 random rough interfaces separating homogeneous media is considered. Table 3 lists the simulation parameters; Fig. 5 shows the scenario.

N i = N i .

NUMERICAL RESULTS

In

To attenuate the edge diffraction by the upper surface, Thorsos's [51] Gaussian-tapered incident wave is applied with g = L 1 /6 = 10λ 0 , and the incident angle is θ inc = 0. The parameter g controls the extend of the incident wave.

To calculate the PBFs, the incident waves [vectors {B i ,k IPW } in Eq. ( 34)] are assumed to be plane, and their incidence angles range from 0 to 2π . In addition, they are spaced equally and their number, N IPW , is assumed to be The rough surface number i obeys a Gaussian height distribution with a height autocorrelation function (ACF) assumed to be Gaussian, with a standard deviation σ z,i and a correlation length L c ,i . Its height mean value (or depth) is h i , the surface lengths are equal to 60λ 0 , where λ 0 is the wavelength in vacuum (medium 0 ), and the sampling step per wavelength λ 0 is equal to 20.

N IPW = k 0 D + 1 n IPW , (41)
Scatterer i σ z,i [λ 0 ] L c,i [λ 0 ] h i [λ 0 ] 0 1 1 2 + 0.
RRE = norm (X LU -X CBFM ) norm (X LU ) . ( 42
)
The norm is the norm 2, which is calculated over the surface 2 shows that the size of the reduced matrix is 4.9 461 times smaller than that of the entire problem and the gain in Fig. 6. Modulus of the surface current currents {ψ 1 , ψ 6 } (upper and lower interfaces) in dB scale versus the surface abscissa (scenario depicted in Fig. 5). The polarization is TM. Fig. 7. NRCS in dB scale versus the scattering angle θ sca (scenario depicted in Fig. 5). The polarization is TM.

saving time is of the order of 5. in Fig. 5. Numerical results, not depicted in the paper, demonstrate that, when R c < λ 0 , the precision of CBFM-KA is lower.

In addition, it is important to underline that the multiple reflections are neglected with KA, which implies that the surface RMS slope must be moderate (typically, smaller than 0.3-0.35).

Figure 9 plots the same variations as in Fig. 7, but the polarization is TE. For LU, CBFM,SVD = 10 -4 , whereas for KA, CBFM,SVD = 10 -6 . As we can see, a good agreement is obtained with LU, which implies that the input parameters are well chosen.

This first study showed that, for a stack of rough interfaces, the CBFM combined with either LU or KA gives very good results. To calculate the PBFs, Eq. ( 41) with n IPW = 3 slightly overestimates the number of incident plane waves. The use of a SVD truncation allows us to reduce this number by a factor ranging from 1.5 to 1.6 (see Table 2 or legends of Figs. 789), which implies that the size of the reduced matrix also decreases.

It is also important to underline that the threshold CBFM,SVD must be divided by 100 for KA in comparison with LU (with CBFM,SVD = 10 -4 ) to obtain a similar value of NIPW,SVD .

B. Scenario 2: Collection of Elliptical Cylinders

In this subsection, the scattering from a collection of 21 elliptical dielectric cylinders is considered. Figure 10 degrees, and the center coordinates are equal to {(-6, -2), (-4, -2), (-2, -2), ... , (2, 2), (4, 2), (6, 2)}λ 0 .

The sampling step per wavelength λ 0 is 20, and the incident wave is assumed to be plane ψ inc (r) = e j k inc •r = e j k 0 (x sin θ inc -z cos θ inc ) , with an incidence angle θ inc = 0.

Figure 11 plots the RCS in dBm scale versus the scattering angle θ sca (scenario depicted in Fig. 10). The polarization is TM. In the legend, the labels mean "CBFM( CBFM,SVD -n IPW )-Method: NIPW -NIPW,SVD , RRE" (Table 1), where the relative residual error (RRE) is defined by Eq. ( 42), in which X is the scattered field in far field versus θ sca . "Method" is the method

C. Scenario 3: Collection of Elliptical Coated Cylinders

In this subsection, the scattering from six elliptical coated dielectric cylinders is considered. Figure 13 shows the scenario. A coated cylinder is composed of three nested cylinders separating four homogeneous media. The scatterers numbered 1, 2, and 3 have semimajor axis a = {4, 3, 2}λ 0 , semiminor axis b = {3, 2, 1}λ 0 , rotation angles {0, 30, 60} degrees, and equal centers of coordinates (-4.5, 7)λ 0 .

The permittivities of the four media { i } (i ∈ (1; 4)) are r = {1, 2 + 0.05 j , 2.5 + 0.06 j , 3 + 0.07 j }. The scatterers {7, 8, 9} and {13, 14, 15} are obtained from scatterers {1, 2, 3} by making a vertical translation of -7λ 0 and -14λ 0 , respectively. The remaining scatterers defined for x > 0 are obtained from those defined for x < 0 by symmetry. The sampling step per wavelength λ 0 is 20, and the incident wave is assumed to be plane ψ inc (r) = e j k inc •r = e j k 0 (x sin θ inc -z cos θ inc ) , with an incidence angle θ inc = 0.

Figure 14 plots the RCS in dBm scale versus the scattering angle θ sca (scenario depicted in Fig. 13). The polarization is TM. As we can see, a good agreement is obtained between LU and CBFM-LU and CBFM-KA. Like previously, for KA, the threshold CBFM,SVD = 10 -6 is divided by 100 in comparison with that of LU to select enough PBFs. In Fig. 14, the legend (see also Figure 17 plots the NRCS in dB scale versus the scattering angle θ sca (scenario depicted in Fig. 16). The polarization is TM. Figure 18 plots the same variation as in Fig. 17, but the polarization is TE. From Eq. ( 41), to calculate N IPW,i , n IPW = 1 surprising and allows us to simplify the calculation of PBFs.

A pq ,mn =      - j k 0 v q ,n | q ,n | 4 H (1) 1 (k0 r q ,n -r p,m ) r q ,n -r p,m γ q ,n (x q ,n -x p,m ) -(z q ,n -z p,m ) for m = n + 1 2 - v q ,n | q ,n | 4π γ q ,n 1+γ 2 q ,n for m = n , ( A1
Table 2 shows that the time difference between CBFM-LU and CBFM-KA is small and increases as the number of unknowns grows. As expected, the time is smaller for CBFM-KA. For CBFM-LU, the PBFs are calculated from a LU inversion of complexity O(N 3 i ), whereas for CBFM-KA, the complexity is O(N i ). Thus, for large N i (3D problem), CBFM-KA should be competitive. In addition, with CBFM, the allocation time to calculate the reduced characteristic matrix (second stage common to CBFM-LU and CBFM-KA) is about 80%-90% of the total time. In Table 2, this explains the small differences between t CBFM,LU and t CBFM,KA .

The advantage of the domain decomposition method is that it is highly parallelizable, which further reduces the computing time. The proposed method is then a powerful electromagnetic computation tool to solve any 2D problem, especially when some scatterers are identical. where r q ,n = (x q ,n , z q ,n ) ∈ S q (coordinates of the point on the surface S q ), r p,m = (x p,m , z p,m ) ∈ S p , γ = d z/d x , γ = d γ /d x , q ,n the sampling step on S q , v q ,n = sgn( nq,n • ẑ)

( nq,n is the unitary vector normal to the surface S q at the point r q ,n ), H (1) 0 the zeroth order Hankel function of the first kind and H (1) 1 its derivative.

For the TE polarization, the variable ρ 01 = 1, whereas for the TM polarization, ρ 01 = r ,0 / r ,1 , where r ,i is the relative permittivity of medium i .

It is important to underline that, for a self (interaction of the same surface) impedance submatrix, p = q and a singularity occurs for m = n, whereas for a coupling matrix p = q , there is no singularity because r p = r q . In addition, any matrix Zpq propagates the field from the source points {r q ,n } toward the observation points {r p,m }.

(Fig. 1 .

 1 Fig. 1. Scattering from a single scatterer.

where ∂ψ 1

 1 /∂n = ∇ψ 1 • n1 is the normal derivative, in which 123 n1 is the unitary vector normal to S 1 . The unknown X 1 on S 1 is 124 computed from X 1 = Z-1 11 b 1 . The scattered field ψsca (r) in the 125 medium r ∈ 0 is then obtained from the Huygens principle 126 expressed as 127 ψsca (r) = -P(r, r 1 )X 1 ,

{

  129 applied by taking the plus sign (instead of minus). The matrices 130 Āij , Bij } are expressed in Appendix A. 131 The matrix P(r, r 1 ) propagates the surface currents 132 {ψ 1 , ∂ψ 1 /∂n} from r 1 to r. Its size is N sca × (2N 1 ), where 133 N sca is the number of observation points and ψsca is a vector of 134 length N sca . 135

  138 direction of observation, γ 1 the slope of S 1 at the point r 1 , 1 139 its sampling step, and v 1 = sgn( n1 • ẑ) (sgn stands for the sign 140 function).

  141

148B.

  Case of P Illuminated Scatterers 149

  10) where 0 is a null matrix. The size of the matrix Z is 161 2(N 1 + N 2 ) × 2(N 1 + N 2 ), where N i (i = {1, 2}) is the 162 number of discretization points on the surface S i of the object. 163 The elements of the submatrices { Ā12 , B12 } can be found in 164 Appendix A. 165 If scatterer 2 is perfectly conducting (metallic), then 166 Z22 = P22 , Z(1) 12 = [ P12 0] T and Z(1) 21 = [ Ā21 B21 ], where 167 Pij = { Āij , Bij } for the TM (∂ψ 2 /∂n = 0) and TE (ψ 2 = 0) 168 polarizations, respectively. In addition, if the scatterer 1 is per-169 fectly conducting, then Z(1) 12,21 = P12,21 and Z11,22 = P11,22 .

  170

  205

Fig. 4 .

 4 Fig. 4. Combination of cases.

  260 tant to note that, if a source exists in the medium 3 , then b 4 261 and b 5 differ from 0. By inverting the matrix Z, the surface 262 currents on the surfaces {S i∈[1;P ] } are computed. 263 The scattered field ψsca (r) in the medium r ∈ 0 264 (r / ∈ (S 1 ∪ S 2 ∪ S 3 )) is then obtained from the Huygens 265 principle expressed as 266 ψsca (r) = -  P(r, r 1 ) P(r, r 2 ) P(r, r 3 )

  270 allow us to generalize to any configuration made up of P scatter-271 ers. To construct the impedance matrix, the following method is 272 proposed. 273 For a given scenario, first a boolean interaction matrix M of 274 size P × P is built. If the scatterer j shares a medium j with 275 the scatterer i, then M ij = 1, 0 otherwise. For the case presented 276 in Fig. 4, M is written as 277

  Fig.2, coupling matrix with the superscript (1)]. Equation[START_REF] Johnson | A numerical study of scattering from an object above a rough surface[END_REF]

  otherwise. For the excitation vector b, if the source belongs to the medium i∈[0;P ] , then the scatterer set J = { j } that shares this medium is illuminated. This means that b J = ψ inc (r ∈ S J ), b J = 01× N J otherwise, where J = J / ∈ [1; P ]. The scattered field in medium i∈[0;P ] is computed by using the same methods as those previously addressed. 3. CHARACTERISTIC BASIS FUNCTION METHOD For P dielectric scatterers, the size of the matrix to be inverted is N × N, where N = 2 i=P i=1 N i . The direct solution of the linear system ZX = b through a direct lower upper (LU) decomposition is usually limited by O(N 3 ) and O(N 2 ) complexities in CPU time and memory requirement, respectively.

  tic square matrix of size (B NIPW,SVD ) 2 instead of a square matrix of size N 2 = (B N) 2 , where N = (1/B) B p=1 N i and NIPW,SVD = (1/B) B p=1 N IPW,SVD, p (mean values over the number of blocks B). Then, the reduction factor is ( N/ NIPW,SVD ) 2 . If multiple excitations {b} (for instance, monostatic case) are considered, then the storing of ZR avoids

  and n i = √ r ,i (i = {1, 2}) the 396 refraction index of the medium i . The unitary vector kinc 397 stands for the incident direction, and n is the unitary vector 398 normal to the surface at the point r. In Eq. (38), I(r) denotes a 399 boolean illumination function. If a point on the surface r ∈ S is 400 viewed (that is, cos θ ≥ 0) by the transmitter, then I(r) = 1; 0 401 otherwise.

  matrix (because only local interactions are accounted for, that is, the multiple reflections are neglected) by blocks. The two upper blocks have N i elements equal to 2I(r m ) and 2I(r m )k inc • nm , respectively, whereas the two lower blocks have N i elements equal to I(r m )[1 + R(θ m )] (θ m depends on r m discretized on S) and I(r m )[1 -R(θ m )] j k inc • nm , respectively.

  this section, numerical examples are exhibited to demonstrate the efficiency of CBFM combined with either LU or KA for the calculation of the PBFs.

Fig. 5 .Figure 6

 56 Fig. 5. Scenario 1: A stack of P = 6 rough interfaces separating homogeneous media. The simulation parameters are listed in Table3.

  452 abscissa x i . CBFM-LU means that the PBFs are computed from 453 a LU decomposition. In addition, the legend "LU" means that 454 the strengths are computed from the brute force MoM (LU 455 decomposition of the entire matrix of the problem, that is, the 456 reference solution). 457 As we can see, for n IPW = 3, the results match well with 458 those obtained from LU, and the comparison is better for 459 CBFM,SVD = 10 -4 , as expected, but the number NIPW,SVD is 460 larger. Table

Figure 6 Figure 7 Figure 8

 678 Figure 7 plots the NRCS in dB scale versus the scattering

Fig. 8 .

 8 Fig.8. NRCS in dB scale versus the scattering angle θ sca (scenario depicted in Fig.5). The PBFs are computed from LU and KA, with n IPW = 3.

  shows the scenario, in which the cylinders are identical but with a different rotation angle. Their permittivities are r ,i = 3 + 0.05 j , semimajor axis a = λ 0 , and semiminor axis b = 2λ 0 . They are numbered from the left to the right going from bottom to top (see Fig. 10). Their rotation angles are spaced equally as {-180, -162, -144, ... , 144, 162, 180}

Fig. 9 .

 9 Fig.9. Same variations as in Fig.7, but the polarization is TE.

Fig. 10 . 542 Figure 12

 1054212 Fig. 10. Scenario 2: Collection of P = 21 elliptical dielectric cylinders.

Fig. 11 .

 11 Fig.11. RCS in dBm scale versus the scattering angle θ sca (scenario depicted in Fig.10). The polarization is TM.

Fig. 12 .

 12 Fig. 12. Same variation as in Fig. 11, but the polarization is TE.

Fig. 13 .

 13 Fig. 13. Scenario 3: Six elliptical dielectric coated cylinders (P = 18 scatterers).

Fig. 14 .

 14 Fig.14. RCS in dBm scale versus the scattering angle θ sca (scenario depicted in Fig.13). The polarization is TM.

Fig. 15 .Fig. 16 .

 1516 Fig. 15. Same variation as in Fig. 14, but the polarization is TE.

Fig. 17 .

 17 Fig.17. NRCS in dB scale versus the scattering angle θ sca (scenario depicted in Fig.16). The polarization is TM.

Fig. 18 .

 18 Fig. 18. Same variation as in Fig. 17, but the polarization is TE.

) 620 addressed for efficiently solving the resulting linear system. 621 Considering four different complex scenarios depicted in 622 Figs. 5 ,

 6216225 Figs. 5, 10, 13 and 16, for both the TM and TE polarizations,

  APPENDIX A: ELEMENTS OF THE MATRICESIn Eq. (1), the elements (m, n) (indexes of the row and column, respectively) of the submatrices { Ā11 , B11 , C11 , D11 } are expressed as ( p = q = 1)B pq ,mn = j | q ,n | 1 + γ 2 0 r q ,nr p,m ) for m = n 1 + 2 j π ln 0.164k 0 1 + γ 2 q ,n | q ,n | for m = n ,(A2)C pq ,mn = A pq ,mn k 0 =k 1 for m = n

Table 1 .

 1 Notations Introduced in This Paper Number of unknowns of scatterer i N IPW, p CBFM plane wavenumber of block i N IPW,SVD, p CBFM plane wavenumber of block i after SVD truncation CBFM,SVD CBFM threshold of the SVD truncation N Mean value of N i over p ∈ [1; B] NIPW Mean value of N IPW,i over i ∈ [1; B] NIPW,SVD Mean value of N IPW,SVD,i over i ∈ [1; B] that the numbers of plane waves {N IPW,i, pq } of { J i, pq } sat-369 isfy N IPW,i,11 = N IPW,i,21 and N IPW,i,12 = N IPW,i,22 . This 370 means that the size ofJ i is 2N i × (N IPW,i,11 + N IPW,i,12). This yields that the new size of J i is 2N i × N i,IPW , with 374 N i,IPW < N IPW,i,11 + N IPW,i,12 . Equations (35) and (36) are IPW , requires solving the linear system Zi ,i Y i ,k IPW =

	Name	Definition
	P	Number of scatterers
	B	Number of blocks (= P )
	N	Total number of unknowns
	N i
			371
	Like a metallic scatterer, the redundant information is elimi-	372
	nated via the use of a truncated singular value decomposition	373
	(SVD). 375
	unchanged.	376
	The SVD decomposition can be applied on each submatrices	377
	J i, pq , and the resulting compressed matrix J i is obtained from	378
	Eq. (37). Numerical tests revealed that this procedure has a lower
			381
	paper.		382
	4. CBFM COMBINED WITH THE KIRCHOFF	383
	APPROXIMATION	384
	For a given plane wave B i ,k IPW , the calculation of a PBF,

379

precision than when the SVD is applied on the whole matrix J i . 380 For convenience, Table 1 lists the notations introduced in this 385 Y i ,k 386 B i ,k IPW , leading to a complexity of O(N 3 i ), where N i is the 387 number of unknowns with overlapping. To reduce the complex-388 ity of this operation to O(N i ), the Kirchoff approximation (KA) 389 is applied. 390 For a dielectric scatterer of surface S separating two homo-391 geneous media ( 1 , 2 ), the surface current and its normal 392 derivative on r ∈ S are expressed as [49,50] 393

Table 2 .

 2 Parameters of the Four Scenarios

a Scenario

CBFM,SVD N N R N R,SVD t LU t CBFM-LU t CBFM-KA

	1	10 -4 14,400 4,392 2,925 98	20	16
	2	10 -5	4,074 1096 1036 15	0.6	0.5
	3	10 -4 11,436 3,336 2,168 48	5	4
	4	10 -5 11,056 3,326 2,690 81	13	11

a

See Table

1

for the notations. In addition, computing times of LU (applied on the entire matrix) t LU , CBFM-LU, t CBFM-LU , and CBFM-KA, t CBFM-KA in seconds. The matrix filling time is included.

Table 3 .

 3 Parameters of the First Scenario: A Stack of P = 6 Random Rough Interfaces Separating Homogeneous Media i of Permittivity r,i

	a

Table 1

 1

	), where the	450

Table 2 .

 2

	546
	547	For a collection of elliptical cylinders, this second study
		showed again that the CBFM is efficient and that the value of

548

N IPW expressed by Eq. (

41

) is well suited. In addition, like for

Table 2

 2

) indicates that NIPW moderately changes in comparison with NIPW,SVD . This means that the value NIPW

Table 2

 2 IPW,i given by 632 Eq.[START_REF] Bourlier | Propagation-inside-layerexpansion method combined with physical optics for scattering by coated cylinders, a rough layer, and an object below a rough surface[END_REF], is well suited. For a rough surface, n IPW = 3 and for an 633 elliptical cylinder, n IPW = 1. To reduce the size of the reduced 634 matrix, a SVD truncation can be applied; for elliptical cylinders, 635 however, it is not useful because N IPW,i,SVD is of the same order 636 of N IPW,i . For a collection of objects invariant by translation 637 and rotation, another advantage of CBFM is that the PBFs 638 and the self-impedance submatrices of the scatterers are equal.

	639
	640	Therefore, the first stage of CBFM is applied only on a single
	641	scatterer.
	642	To accelerate the computation of the PBFs, that is, to avoid
		an LU decomposition on the self-impedance submatrices, KA is

627

summarizes the performances of CBFM-LU. For a large num-628 ber of unknowns, N, the efficiency is even better since only 629 the self-impedance submatrices are stored for CBFM and the 630 complexity of LU is O(N 3 ).

631

The calculation of the number of PBFs, N 643 applied. The numerical results showed that CBFM-KA matches 644 well with LU, but in comparison with CBFM-LU, the threshold 645 CBFM,SVD must be divided by 100 to obtain a comparable value 646 of N IPW,i . It is important to point out that CBFM-KA predicts 647 good results, even on geometries for which KA is not valid, like 648 the elliptical cylinders. In addition, in Eq. (

38

), the boolean 649 illumination function I(r) can be omitted. This statement is 650

Vol. 38, No. 9 / September 2021 / Journal of the Optical Society of America A

Eq. (11), the vector b 2 = 0 because, as shown in Fig.1, the

F. Frezza, F. Mangini, and N. Tedeschi, "Introduction to electromag-

netic scattering, part ii: tutorial," J. Opt. Soc. Am. A 37, 1300-1315

(2021).

G. D. Tsogkas, J. A. Roumeliotis, and S. P. Savaidis, "Electromagnetic

scattering by an infinite elliptic dielectric cylinder with small eccen-

tricity using perturbative analysis," IEEE Trans. Antennas Propag. 58,

107-205 (2010).

G. P. Zouros, "Electromagnetic plane wave scattering by arbitrar-

C. Amra, "Light scattering from multilayer optics. I. Tools of investiga-

Soc. Am. 83, 78-92 (1988).

F. T. Ulaby, R. K. Moore, and A. K. Fung, Microwave Remote Sensing:

Active and Passive (1982), Vol. 2.

R. J. Papa and J. F. Lennon, "Conditions for the validity of physical

optics in rough surface scattering," IEEE Trans. Antennas Propag. 36,

647-650 (1988).

Appendix A. 116 The vector b 1 is defined as

Research Article

For the calculation of the scattered field in the far field 690 (u → ∞), the following expansions [START_REF] Abramowitz | Handbook of Mathematical 713 Functions[END_REF] of the Hankel functions 691 can be applied: 692 H (1) 0 (u)

857

surface scattering using a Gaussian roughness spectrum," J. Acoust.


Mohamed Akkari

Gérard Subsol

Nabil Zemiti

Lucas Lavenir

Charlotte Farah

Florence François

Michel Mondain

Guillaume Captier

Philippe Poignet

Frédéric Venail

Feasibility of Cochlea High-frequency Ultrasound and Microcomputed Tomography Registration for Cochlear Computer-assisted Surgery: A Testbed