Modelling Scattering of Electromagnetic Waves in Layered Media: An Up-to-Date Perspective

1 Institute for Electromagnetic Sensing of the Environment (IREA), National Research Council (CNR), Napoli, Italy 2Department of Electrical Engineering and Information Technology, University of Naples Federico II, Napoli, Italy 3Department of Electrical, Electronic, Telecommunications Engineering and Naval Architecture (DITEN), University of Genoa, Genova, Italy 4Icam, School of Engineering, Nantes Campus, Carquefou, France


Introduction
The problem of electromagnetic (EM) wave scattering in layered or stratified media has become an extremely important subject with theoretical and practical relevance.
Indeed, most of the real structures of interest, both those occurring naturally and those fabricated artificially, can be reasonably assimilated to layered structures to some degree. As a matter of fact, the solution of Maxwell equations in such structures poses serious difficulties of a mathematical nature, and there is no general and uniform approach. In fact, methods of studying scattering phenomena in layered media are greatly diversified. They depend on the kind and description of structure (objects embedded in layers, inhomogeneities distributed in a continuous manner or in the form of randomly distributed discrete scattering elements, etc.), on the information at our disposal concerning the structure of the medium, and on the kind of information that is sought about the wave process in question. Depending on the application context, the investigation on wave phenomena in layered structure, in some cases, can be successfully conducted by resorting to deterministic description; however, as far as natural scenarios are concerned, a stochastic description can provide a more adequate description of reality.
Within this framework, a clear understanding of the wave scattering processes taking place in the layered structures still poses great theoretical challenges. Accordingly, modelling EM scattering phenomena in stratified media is an active area of research with many practical applications, and there exists a very extensive literature on the subject. The inherent problem formulation can be treated from several viewpoints. As a first distinction, the existing techniques can be grouped into two main classes: direct and inverse approaches. Direct scattering models for layered structures have a practical relevance in a number of contexts, such as radio-wave propagation, optics, radar imaging, and microwave remote sensing. Direct modelling methods can generally be categorized in analytical and numerical methods. Some of them turn out more appropriate, in terms of accuracy and computational cost, in a specific application context. Stratified structures play a paramount role also in the solution of electromagnetic inverse scattering problems, which are involved in imaging of hidden or buried targets in multilayer media. Nondestructive testing and evaluations represent significant example, as well 2 International Journal of Antennas and Propagation as the modelling of layered structures in through-the-wall procedures for security applications and in the geophysical prospecting of buried cultural heritages.
Numerous recent publications, in the last decades, witness an increasing interest in the study of scattering interactions in layered structures, with advancements including both analytical and computational approaches. Nonetheless, the significant progress on this topic, which has also been driven by emerging applications, deserves to be framed and discussed in an organized way.
Therefore, this paper aims at providing a concise and organized exposition of the existing methods of analysis of EM scattering in layered media. Moreover, the emphasis is placed on conceptual advancements and novel methods that have been recently established, thus outlining important new research directions. The different methodologies are classified on the basis of their characteristics, also discussing pertinent basic principles, advantages, and disadvantages.
The paper is organized as follows. Section 2 focused on existing analytical formulations, with a particular emphasis on recently developed functional forms for rough multilayer media scattering. Numerical methods of interest are addressed in Section 3, thus covering both well-established and innovative methods. Section 4 is devoted to EM inverse scattering problems in layered structures, also delineating current progress and trends.

Analytical Formulations
The analytical evaluation of the scattering in layered media has received considerable attention in last decades , due to its crucial role played in the different applications. Indeed, the analytical approach permits an understanding of the functional dependence of the scattering response on the structure parameters and, in some cases, an intelligible explanation of the underlying phenomenon.
Conversely, the regime of validity inherent to the approximate solutions implies a limited application domain. In the following, the emphasis is placed on scattering from layered media with rough interfaces. In particular, we consider the two main classical approaches for scattering from rough dielectric surface that have recently been extended to deal with layered configurations. The first classical approach is provided by the perturbative theory. At the end of 19th century, Lord Rayleigh originally proposed for the first time the small perturbations method (SPM) for the description of wave scattering from a surface separating two media. SPM has been further developed by a number of authors [1][2][3]. The second traditional approach is referred to as the Kirchhoff approximation, in which the tangent plane approximation is used to compute the tangential fields at the interfaces [1,2]. For each approach, the corresponding extension to a layered configuration is extensively discussed in the following. Further approaches are also considered in Section 2.3.

Perturbative Methods.
We refer to the 3D rough multilayered structure depicted in Figure 1. Each layer is assumed to be homogeneous and characterized by deterministic parameters: the dielectric relative permittivity and the thickness Δ = − −1 . In general, the interfacial roughness is assumed to be described by a random function of the space coordinates.
In such a case, an approximated solution for the actual structure, which can be described by small changes with respect to an idealized (unperturbed) structure, is obtained by suitably taking advantage of the exact solution available for the associated unperturbed problem.
Within this framework, two different systematic formulations have recently been introduced to deal with the analysis of a layered structure with an arbitrary number of rough interfaces. Specifically, the results of the boundary perturbation theory (BPT) [4][5][6] lead to polarimetric, formally symmetric, and physically revealing closed-form analytical solutions. In this case, a suitable perturbation pertinent to the structure geometry is concerned. The volumetric-perturbative reciprocal theory (VPRT) for the evaluation of the scattering from a layered structure with an arbitrary number of rough interfaces considers a perturbation pertinent to the dielectric properties of the structure [7][8][9]. It is important to note that both BPT and VPRT lead to formally identical expressions of the first-order scattered field. The BPT/VRPT fist-order solution is hereinafter detailed for a monostatic configuration in terms of normalized radar scattering cross section (NRCS). Assuming a -polarized incident wave impinging on the structure from the upper half space (see Figure 1), the NRCS of the layered medium with corrugated interfaces can be expressed in compact closed-form as follows: whereiñ, +1 pp denotes the polarimetric coefficient pertinent to the rough interface between the layers and + 1, whose expressions arẽ , +1 hh where = √ 2 − ( ⊥ ) 2 , and the following notation has been used: direction through the regions 0 and , respectively. They can be expressed by recursive relations in terms of the ordinary transmission and reflection coefficients, like in [6,10]. Moreover, in (1)-(3) the adopted notation is as follows: the asterisk symbolizes the complex conjugated operator; Re{ } denotes the real part operator; ∈ {V, ℎ} indicates the incident (vertical or horizontal) polarization; k ⊥ represents the 2D projection of incident wavenumber vector on the plane = 0 (superscript referring to the incident field direction); and ⊥ = |k ⊥ |. Finally, ( ) and ( ) denote the (spatial) power spectral density of th corrugated interface and the cross power spectral density (between the interfaces and ), respectively [1,4].
We stress that (1)-(3) represent the backscattering solution; a more general expression valid for a bistatic configuration can be found in [4,6,7]. As far as the backscattering case is concerned, it is should be noted that the BPT/VPRT crosspolarized scattering coefficients (̃0 with ̸ = ) evaluated in the plane of incidence vanish, in full accordance with the classical first-order SPM method for a rough surface between two different media. Moreover, the dual solution concerning the scattering through the rough multilayer has also been obtained in [5,6].
The domain of validity of the perturbative solution is defined as follows: the height deviation of the rough interfaces, about the unperturbed interface, is everywhere small compared to the wavelength of the incoming wave and the gradient of the interface is small in comparison to unity [11]. A comprehensive discussion on the relevant domain of applicability is provided in [12,13]. Concerning the considered BPT/VPRT solutions, several comments are in order.
First, we discuss the inherent functional form. The closedform solution (1)-(3) makes the functional dependence of scattering properties on layered structure (geometric and electromagnetic) parameters explicit. Specifically, the analytical solution permits the evaluation of the backscattering from the layered rough structure, once the 3D layered structure parameters (shape of the roughness spectra, layer thickness, and complex permittivity) and the incident field parameters (frequency, polarization, and direction of incidence) have been specified. Therefore, the considered perturbative solution provides a forward solver whose recursive formulation results in effective computational time. Furthermore, the NRCS of the layered media is sensitive to the correlation between rough profiles of different interfaces. In particular, when completely uncorrelated random interfaces are concerned ( ( ) = 0), scattered intensity (in the first-order approximation) arises from the incoherent superposition of the field contributions scattered from each rough interface [see (1)]. Interface roughness description (and the corresponding spectral representation) is typically provided in terms of classical parameters (height standard deviation and correlation length); however, it can also be expressed in terms of the fractal parameters [14,15]. Indeed, fractal description for the interfacial roughness has also been adopted for evaluating scattering from multilayered structures (see also Section 2.2.3). A fractional Brownian motion (fBm) process based description has been used in conjunction with the BPT scattering solution in [16].
Second, we provide a perspective on the interpretability of the perturbative solutions. The physical meaning of the perturbative solutions for the scattering from and through the 3D layered structure has been investigated in [17,18], where a physically revealing interpretation involving ray-series representation is obtained by rigorously establishing a functional decomposition of the first-order scattering solutions in terms of basic single-scattering local processes. Accordingly, the fundamental interactions in the multilayer contemplated by the mathematical solutions can be revealed, thus gaining a neat picture of the physical meaning of the theoretical construct. Furthermore, the VPRT procedure can also be reformulated in a more physically sound way by avoiding use of Dirac delta function and distribution theory. Such VPRT reformulation also enables an interesting interpretation in 4 International Journal of Antennas and Propagation terms of internal field approximation, which is consistent with the gently rough assumption, thus providing an additional insight into the first-order approximation (which is slightly different from the usual Born approximation) [1,12]. Another interesting interpretation of the scattering solution can be given in terms of multireaction, by exploiting the fundamental concept of Rumsey reaction [7,10].
Third, the theoretical and experimental consistencies of the perturbative solution deserve to be discussed. Regarding the former, it is important to highlight that all the previous existing perturbative scattering models, introduced by different authors to deal with some simplified layered geometry, with one [19][20][21] or two [22] rough interfaces, can be all rigorously regarded as a special cases of the general BPT/ VPRT solutions [4,7,23]. Accordingly, BPT/VPRT results can be regarded as the generalization of the classical firstorder SPM for rough surface to the case of rough multilayer. As a matter of fact, the validity of a mathematical model resides in the concordance of model-based predictions with experimental observations. In this regard, it is worth mentioning that a first assessment of the BTP solution validity in a real scenario has been provided in [24], thus showing that theoretical predictions are in good accordance with the experimental evidence.
Second and higher-order perturbative developments have also been investigated [9,[26][27][28][29]. It is worth highlighting that pertinent analytical developments following a boundary perturbation approach can be particularly cumbersome [27]; conversely, rigorous second-order volumetric-perturbative developments demand for an appropriate mathematical playground (the distribution theory for discontinuous test functions) [9]. In particular, in [27] the predictions of the fourth-order perturbation have been examined for scattering from two rough surfaces in a layered geometry. Accordingly, the interaction effects between the two surfaces can, in some cases, be the dominant contribution to cross-pol returns. Notice that the first-order solution for the scattered field does not predict the cross-polarized intensity in the incidence plane. In [29] it has been demonstrated that there exists a "strong" condition of energy conservation in that the kernel functions multiplying the spectral density of each interface obey energy conservation exactly. Accordingly, the energy is conserved independent of the roughness spectral densities of the rough surfaces. Finally, we observe that the development of perturbative approaches enabling systematical evaluation of both interfacial and volumetric inhomogeneity in layered media is a main challenge.

Kirchhoff Approach/PO.
Here, focus is made on the socalled Kirchhoff-tangent plane approximation, often called Kirchhoff approximation (KA) for short, in which the tangent plane approximation is used to calculate the tangential fields at the interfaces. It is often rather called physical optics (PO) approximation by numericians or in the radar community [30,31]. This approximation is valid for surfaces that can be considered as locally flat. Thus, at moderate angles, the validity domain is such that the mean surface curvature radius is significantly greater than the EM wavelength. It may then be considered as a high-frequency approximation.
As a consequence, at each surface point, the surface can be replaced by its tangent plane, which is a flat surface whose local slope vector is that of the original random rough surface at considered surface point. Then, each ray of the incident wave is reflected/transmitted in the appropriate local specular direction.
However, it must be emphasized that this model becomes less valid for angles becoming lowly grazing (in particular, for low-grazing incidence): it is due to the phenomenon of shadowing of the surface. In such configurations, it is necessary to introduce a corrective parameter called shadowing function (or illumination function) to overcome this issue [32]. Besides, unless the process of scattered field calculation is iterated, the KA does not take the phenomenon of multiple scattering into account. This reduces the validity of the model to surfaces having small to moderate slopes, typically RMS slopes less than about 0.3-0.5 [33].
Starting from the calculation of the scattered field, we are usually interested in calculating a mean scattered intensity. Then, the KA is often not used as such, as it is usually desired to further simplify the mathematical expression of the average scattered intensity, in order to get faster numerical results. It is particularly true here when we are interested in analytical formulations. Usually, two opposite solutions are elected: the first one considers very rough surfaces (the geometric optics (GO) approximation) and the other one slightly rough surfaces (the scalar Kirchhoff approximation (SKA)). Then, the GO calculates only the incoherent scattered intensity, whereas the SKA generally focuses on the calculation of the coherent scattered intensity.

GO Model.
Let us first focus on the GO model. The classical GO model for single reflection onto a rough surface uses the so-called method of stationary phase (MSP) to simplify the expression of the scattered field (or intensity). This method assumes that only the surface points that specularly reflect/transmit the incident wave into the observation direction contribute to the scattering process. This restricts the contribution to surface slopes checking = −Q H / , with Q = k r − k i called the Ewald vector (k i and k r being the incidence and reflection wave vectors, resp.), Q H and being the horizontal and vertical components of Q, respectively. From the MSP, the GO model further simplifies the scattered intensity by retaining only the surface points that are highly correlated. Thus, the NRCS under the GO is only incoherent and can be expressed as follows [34]: with K being the so-called Kirchhoff kernel, which is a polarization-dependent term proportional to the Fresnel reflection coefficient and related to the projection of the incident wave vector onto the normal to the surface. is the shadowing function [32].
International Journal of Antennas and Propagation The extension of the GO to a layer of two random rough interfaces has been first led for 2D problems [35,36]. By doing so, the NRCS associated with a given order of reflection from the layer (see Figure 2) has the great advantage of being proportional to the product of elementary NRCS corresponding to each scattering in reflection or transmission at a given interface. It may be written in the following form for the second-order NRCS (associated with the field 2 in Figure 2) [35,36]: with 0 12 being the NRCS associated with the transmission of the incident wave into the inner layer Ω 2 , 0 23 the NRCS associated with the reflection of this wave onto the lower interface Σ 2 , and 0 21 the NRCS associated with the transmission of this wave back into the incidence medium Ω 1 . Then, it has been shown that this approach can be generalized to any order n of reflection from the layer (associated with the field in Figure 2) as [36] Then, the number of numerical integrations to be computed is equal to − 1. For the more general case of 3D problems, unfortunately it is not possible to write the NRCS in such a compact way of a product between NRCS, owing to the fact that the polarization terms are not scalar anymore but become square matrices of length 2. This is due to the co-and cross-polarization terms that must be considered for having a correct computation of the physical contributions to the scattering process [37,38]. This approach has been validated by a reference numerical method based on the method of moments [37,38].

SKA Model.
Starting from the KA for calculating the scattered field, the scalar Kirchhoff approximation [39], also called zero-order KA [40], is an alternative to the GO model. This approximation is valid for calculating the scattered field around nadir for slightly rough surfaces. It is generally used for evaluating the coherent scattered intensity.
In the context of scattering from random rough layers, each order contribution must be calculated, or at least the first orders. The first order is well known and corresponds to the scattering from a single rough interface. Under this model, it has been shown that taking the interface roughness into account is simple and consists in replacing the Fresnel reflection coefficient 12 by a modified reflection coefficient 2 ) with ℎ1 being the RMS height of the upper interface Σ 1 and 1 the vertical component of the propagation vector inside Ω 1 (see Figure 2). The term 1 ℎ1 is sometimes referred to as the Rayleigh roughness parameter, which accounts for the attenuation of the field amplitude due to the interface roughness. The extension to the calculation of the secondorder contribution leads to an attenuation of [39,40] exp (−2 2 2 ℎ2 with ℎ2 being the RMS height of the interface Σ 2 and 2 the vertical component of the propagation vector inside Ω 2 . Then, by calculating each order contribution, it is possible to obtain a new equivalent reflection coefficient̃e q of the form [40]̃e 2 ) and 2 being the mean thickness of the layer Ω 2 . For perfectly flat interfaces ℎ1 = ℎ2 = 0, this reduces to the classical formula of generalized reflection coefficient R 1|2 [4,10,40]: This approach has been validated by comparison with numerical methods [39,40].

KA Model and Fractal Surfaces.
It is finally worth mentioning that the KA can be adapted to evaluate electromagnetic scattering from fractal surfaces. Fractal geometry was introduced by Mandelbrot in the '70s of the last century to provide a mathematical tool for the complex and irregular shapes of natural objects [14]. Effectiveness of fractals to describe natural surfaces has been demonstrated in a very impressive way by producing surprisingly realistic computergenerated synthetic landscapes [14]. Accordingly, fractal models, and in particular the fractional Brownian motion (fBm) model, have been extensively used in different fields to describe natural surfaces. The reason of this success is the ability of fractal models to properly account for the statistical scale invariance properties (self-affinity) of natural surfaces. We recall that a set is self-similar if it is invariant (possibly in a statistical sense) with respect to a transformation in which all the coordinates are scaled down by the same factor, whereas it is self-affine if it is invariant (possibly in a statistical sense) with respect to a transformation in which coordinates are scaled down by factors not all equal. Self-affinity of fractal sets is the key property that makes them particularly useful in describing natural surfaces. However, while fractal sets maintain their self-affinity at any (arbitrarily small or large) observation scale, natural surfaces have characteristic inner and outer scales. In other words, they exhibit fractal characteristics only on a wide but limited range of scale lengths (range of fractalness). When an electromagnetic wave impinges on a natural surface, the range of scale lengths involved in the scattering mechanism is limited on one side by the finite dimension of the illuminated surface and on the other by the electromagnetic wavelength. If this range of scales is included in the surface range of fractalness, then the fractal description is appropriate (and can be used) to evaluate the scattered field. This is often the case at microwave frequencies [41].
First studies on application of KA to fBm surfaces (sometimes approximated by Weierstrass-Mandelbrot functions [42]) date back to the last decades of last century [15,41,43,44]. To summarize the obtained results, we have to preliminarily recall that an fBm surface is a 2D random process whose increments over a fixed distance are zeromean Gaussian with variance where is the Hurst coefficient, linked to the surface fractal dimension by the relation = 3 − , and is the surface topothesy. According to KA, the mean square value of the modulus of the field scattered by an fBm surface is [41] The integral in (12) can be evaluated via two different series expansions: one is an asymptotic expansion of the scattering integral for near-specular directions or high frequency and the other for far-from-specular direction or low frequency. The behavior of these series and their practical evaluation are analyzed in detail in [45]. It has been also highlighted that the integral in (12) is in the form of the PDF of an alpha-stable process [46][47][48], and a physical interpretation of this result was provided in [48].

Other Analytical Approaches.
For the sake of completeness, in the following we briefly consider further approaches adopting different schematizations and/or approximations, which deserve to be discussed. The small-slope approximation (SSA) was originally proposed by Voronovich [30,49]. SSA exhibits an extended domain of applicability, which is aimed at including the domain of the small-perturbation method and the domain of the Kirchhoff approximation. SSA extended to the fourthorder terms of the perturbative development has been considered in [50] for studying slab with uncorrelated rough 2D interfaces. A first-order SSA method applied to an arbitrary number of 1D interfaces has been presented in [51]. In particular, the scattering amplitudes under the first-order small-slope approximation are deduced from results derived from the first-order SPM.
Generally speaking, modelling techniques for describing electromagnetic scattering in a random layered medium can be categorized in two main classes: continuous and discrete. For the continuous case, exact equations for the mean and correlation of the electric field are referred to as the Dyson and Bethe-Salteper equations, respectively [1]. These equations are difficult to be solved even under ideal conditions, thus simplified assumptions have to be introduced in order to achieve tractable forms. Discrete scattering models, for describing the response of a layer containing a random collection of discrete scatters (e.g., dielectric cylinders or disks), have also been developed (see, e.g., [52][53][54][55][56]). In particular, the scattering coefficients of the layer can be computed using a distorted Born approximation, with the mean wave evaluated according to the Foldy-Lax approximation [1].
However, the consistent treatment of interfacial and volumetric inhomogeneity in layered media is still posing a challenge.

Numerical Methods
In this section, the focus is on the main numerical methods of interest, which exhibit different degrees of approximation and computational complexity . In particular, we first consider the rigorous numerical techniques referred to as the method of moments (MoM), and then efficient approaches are discussed.

Method of Moments (MoM).
As already mentioned, the electromagnetic field scattered by the rough interfaces of a layered medium can be computed from knowledge of the tangential components of the electric and magnetic fields at the interfaces (i.e., equivalent magnetic and electric surface currents, resp.). To compute such tangential field components, integral equations must be solved [82]. For the complete equations to be solved in the single interface case, the reader is referred to [82] for full 3D problems and to [83] for 2D cylindrical problems; for the equations that apply to the case of two or more interfaces, see Section 3.3 and references therein. Here, to simplify the notation, we report the scalar version of those integral equations, which is a Fredholm equation: ∬ (r, r ) (r ) r = (r) (13) in which is the interface (or a set of interfaces), r and r are two points over , (r) is known (it is a component of the incident field), (r, r ) is an element of the dyadic Green function, and (r ) is the unknown (a component of the tangential fields). This integral equation can be solved numerically by using the method of moments (MoM). The unknown function is expanded by using a finite set of basis functions (r ) and the inner product of both sides of (13) with test, or weighting, functions (r) is evaluated: We then get the system of linear algebraic equations = (16) in which = ∬ (r) ∬ (r, r ) (r ) r r, = ∬ (r) (r) r, (17) and are the unknowns. Once the system is solved for , the solution of the original integral equation can be obtained from (14).
If (r) = (r), that is, the sets of basis and test functions coincide, then we get the Galerkin method. A usual choice is to use rectangular pulse basis functions and Dirac pulse test functions, so obtaining the so-called point-matching method. This choice simplifies evaluations of matrix and vector elements, but it requires a high number of unknowns . Different choices of basis and test functions may lead to a smaller number of unknowns and/or to a sparser matrix. In the framework of the problem of scattering from rough interfaces, the point-matching method is the most used [64-70, 72, 74, 83], but a very popular choice is also the Rao-Wilton-Glisson (RWG) basis [84]. More recently, use of Haar [85], B-spline [85], and Coifman [62,86,87] wavelets as basis and test functions in the Galerkin method has also been demonstrated to lead to efficient evaluation of scattered field.

Iterative
Methods. As shown above, by using the MoM, integral equations to be solved to evaluate induced (real or equivalent) surface currents over the (possibly) rough interfaces are converted into a system of linear algebraic equations, whose general form is given by (16). However, direct methods for the solution of linear systems have a high computational cost. In fact, they require a number of computer operations asymptotically increasing as 3 (i.e., ( 3 )), where is the (usually very large) number of unknowns (i.e., for the point-matching method, of surface points). To reduce computational complexity to ( 2 ), iterative methods can be used. Linear iterative methods can be cast in the following general form: where ( ) and ( ) are linear functions of and . If they actually depend on , that is, they are adaptively changed at each iteration, the iterative method is called nonstationary. Different nonstationary methods have been applied to rough interface scattering problem, namely, the conjugate gradient squared (CGS), biconjugate gradient-stable (BICGSTAB), quasiminimum residual (QMR), general minimal residual (GMRES), and conjugate gradient-normal equation (CGNR) [58][59][60][61][62]. Conversely, if and do not depend on , the iterative method is said to be stationary. In addition, if we let then it can be shown that the iterative procedure (19) converges to the exact solution of (16) if and only if all the eigenvalues of lie within the unit circle on the complex plane [63]. Different choices for the "splitting matrix" define the different iterative stationary methods. For instance, = , where is the diagonal part of , defines the Jacobi (or simple iteration) method. In addition, for we get the symmetric successive overrelaxation (SSOR) and the 2 × 2-block-SSOR methods (with unitary relaxation parameter ), respectively. In (20)- (21), − and − are the lower and upper triangular parts of , and ( ) , − ( ) , and − ( ) are the 2 × 2 block-diagonal, block-lower triangular, and block-upper triangular parts of [63]. It turns out that nonstationary methods are in general more robust, meaning that they converge for a wider range of roughness conditions. However, stationary methods have the advantage of having a physical interpretation when applied to the problem of scattering from rough surfaces, so that it is easier to predict interface roughness conditions under which they rapidly converge. In fact, the Jacobi method coincides with the iterative (or "extended") Kirchhoff approach (IKA) [64] for perfectly conducting surfaces, whereas SSOR and block-SSOR methods are equivalent to the forward-backward (FB) (or the differently formulated, but equivalent, method of ordered multiple interactions, MOMI) for perfectly conducting [65,66] and dielectric [67] interfaces, respectively. It is interesting to note that IKA, FB, and MOMI were first devised in [64][65][66] based on physical considerations, and only subsequently it was shown that they are equivalent to already available stationary methods (in [61] for perfectly conducting surfaces and in [68] for dielectric interfaces).
To illustrate the physical meaning of FB, we need to recall that, as shown in [68], (18)- (19) with (21) can be recast as with ,(0) = 0, where = + . It can be shown [67] that solution of (22)-(23) does not involve any matrix inversion, except for the inversion of N 2 × 2 matrices. Analysis of these equations shows that, at the first step ( = 1), (22) allows calculating the surface currents due to forward scattered waves and (23) based on these currents computes the new currents that are due to both forward and backward scattered waves. Accordingly, computed currents after the first step include all orders of multiple scattering which involve no backwardforward change of direction and no more than one forwardbackward change of direction. By iterating the reasoning, it is not difficult to realise that computed currents after the -th step include all orders of multiple scattering which involve no more than − 1 backward-forward changes of direction and no more than forward-backward changes of direction. Since scattered energy is not expected to change direction more than a few times (except that for extremely rough surfaces), we conclude that the iterative procedure should converge very rapidly. This expectation is experimentally confirmed for 1D nonreentrant rough profiles [67][68][69]; for 2D rough surfaces, proper ordering of matrix elements must be performed to obtain convergence of the method [70], even for moderate roughness.
Extension of FB to layered structures with rough interfaces is illustrated in [74] and can be possibly coupled with methods presented in the next subsection. Both stationary and nonstationary iterative techniques can be further accelerated by identifying a strong-interaction region ( elements with | -| smaller than a prescribed threshold) and a weak-interaction region (all other elements of ) and using some form of approximation when summing up terms corresponding to the weak-interaction region in (18). In this way, it is possible to reduce computational complexity up to ( log ) or even ( ). This class of algorithms includes, for instance, the banded-matrix iterative approach/canonical grid method (BMIA/CG) [71], the spectrally accelerated (SA) BCGSTAB [60], and SA-FB [72]. An approach with some analogies with algorithms of this class is the BIE/SDIM (Boundary Integral Equation/Subdomain Decomposition Iterative Method), which has been recently proposed [77]. It is a rigorous method based on MoM, which splits the problem into subdomains. Similar to the PILE method (see the description in Section 3.3), the diagonal blocks of the global impedance matrix correspond to the subdomain impedance matrices, whereas the off-diagonal blocks correspond to coupling matrices accounting for the interactions between two different subdomains. The principle of SDIM is to invert the impedance matrix by blocks, which significantly reduces the complexity by comparison with a direct LU. In addition, to accelerate the matrix-vector products and to reduce the memory requirement, the adaptive cross approximation (ACA) is applied to compress the subdomain coupling matrices.
Finally, convergence properties of iterative methods can be in some cases improved by using proper preconditioning, that is, by left-multiplying both sides of (16) by the inverse of a preconditioner matrix; see, for example, [73].

Efficient Approaches.
Resolving the scattering by a multilayered medium with a MoM-based method may be time consuming for large surfaces and/or high frequencies, in particular for several interfaces. Then, it is of interest to apply efficient approaches that accelerate the computing process. Most approaches have been described in previous subsection.
Based on the MoM, the main feature of PILE is its ability to split the total scattered field into the contribution of each order of reflection from a layer (see 1 , 2 , . . . , in Figure 2). It has the great advantage of being able to be coupled with acceleration methods existing for scattering from single interfaces, like the BMIA/CG or the FB [76,77].
We start from the linear system in (16), where is a vector of length 2( 1 + 2 ) representing the unknown total field on the surfaces Σ 1 and Σ 2 ( 1 and 2 are the number of surface samples for the MoM calculation on Σ 1 and Σ 2 , resp.). Then, can be expressed as with 1 and 2 being vectors of length 1 and 2 , respectively, containing the unknown fields and their normal derivatives on Σ 1 and Σ 2 , respectively. is a source term vector of length 2( 1 + 2 ) such that where 1 and 2 are vectors of length 2 1 and 2 2 , respectively, containing information of the incident field on Σ 1 and Σ 2 , respectively (inside Ω 1 and Ω 2 , resp.; see Figure 2). Here, there is no incident field inside Ω 2 ; hence 2 = 0. This also holds for the last 1 terms of 1 . is a square impedance matrix of size 2( 1 + 2 ), with = ( 1 21 12 2 ) .
International Journal of Antennas and Propagation 9 is the impedance matrix of size 2 ×2 of surface Σ (with = {1, 2}), and 21 of size 2 1 ×2 2 (for the propagation from Σ 2 to Σ 1 ) and 12 of size 2 2 × 2 1 (for the propagation from Σ 1 to Σ 2 ) are coupling matrices between the two surfaces.
To effectively solve the system = , the PILE method has been developed [79]. It is based on an inversion by blocks of the impedance matrix, by making a decomposition of domains from the Taylor series expansion of the inverse of the Schur complement. This leads to where ( and 1 = − −1 2 12 1 . The above series expansion is then valid if ‖ 1 ‖ < 1, with ‖ ⋅ ⋅ ⋅ ‖ the norm operator, which is also here a so-called spectral radius, that is, the modulus of the highest modulus eigenvalue. Physically, the total currents 1 on the upper surface are the sum of the contributions ( ) 1 , which correspond to successive iterations , corresponding to each order of reflection from the layer = + 1 (with ∈ {1, 2, . . . , } in Figure 2). Thus, by construction, PILE computes the scattered field associated with each order of reflection from the upper layer.
By following the same principle, an extension to 3 interfaces and even a generalization to any number of interfaces have been developed [81]. For more details, the reader is referred to Appendix A of [81].

Inverse Scattering Approaches
Layered structures have been widely considered in the solution of electromagnetic inverse scattering problems for imaging purposes. Buried and concealed targets can be inspected by using inverse scattering concepts in several application areas, ranging from nondestructive evaluation and testing to civil engineering, security and military applications, geophysical prospecting, and biomedical diagnostics [88][89][90].
A possible imaging configuration is sketched in Figure 3, where a target is buried in a stratified medium and is illuminated by a source located in the upper half space (noninvasive imaging). One or more antennas or sensors are used to receive the field scattered by the target. It should be considered that several works have been focused on the detection of the dielectric properties of the various layers (stratified media) [10,91,92], whereas for the imaging of concealed targets, the propagating structure is usually assumed to be made of known layers (e.g., in most through-the-wall inspection techniques [93]). The reconstruction process can be performed in time-domain [94] or in the spatial-domain [95]. Timereversal approaches have also been adopted. In most cases, however, time-harmonic fields are used, with reference to both single-and multifrequency processing [96]. From a general point of view, another classification holds: approaches based on radar and beamforming concepts [97,98], such as those adopted in breast cancer imaging, and methods based on inverse scattering [90]. In geophysical applications, both methods can be used to inspect shallow buried targets in order to extend the diagnostic capabilities of the socalled ground penetrating radar [99]. Moreover, the inversion procedure can be based on deterministic [100] or stochastic concepts [101].
Considering approaches based on inverse scattering in harmonic fields, the following electric field integral equation (EFIE) can be considered [90] (it is worth noting that other approaches can be followed, e.g., those based on the so-called contrast source formulation [102]): in this equation, E , E , and E denote the incident, scattered, and total electric field vectors, respectively. The scattering potential is indicated by and contains the information about the dielectric properties of the unknown target. Finally, G str is proper Green's dyadic function for the considered buried or stratified configuration [103].
As it is well known, this equation is nonlinear and severely ill-posed [90]. It can be linearized in the case of presence of weakly scattering targets. One of the most used approximations is the Born one [104], for which the scattered field due to the targets inside the investigation area is expressed in terms of the incident field only. In some case, the use of the second-order Born approximation has been proposed, too [105]. Also the Rytov approximation, which is applied to the complex phase of the field, has been used [106].
It should be mentioned that the Born approximation can be combined with a numerical evaluation of Green's function resulting in a very effective iterative reconstruction procedure, called the distorted-wave Born approximation [107], which is widely used with a great success in several imaging approaches.
Methods based on certain approximations on the considered model can be called qualitative methods [108,109]. They are aimed at retrieving some fundamental properties of a target (e.g., position and shape inside the test area), without providing complete maps of the dielectric properties of the targets. So far, we mentioned approaches for the imaging of penetrable targets. Different approaches can be adopted for conducting objects [110,111].
Gauss-Newton methods seem to be very effective deterministic techniques for inspecting buried or stratified targets. Recently, they have been developed also with reference to Lp Banach spaces [112,113], where the oversmoothing effects usually associated with regularization procedures developed in the common Hilbert space seem to be mitigated. It should be also noticed that sparsity concepts and compressive sensing can be successfully adopted [114,115].
An example of reconstruction results is reported in the following. The inexact-Newton procedure developed in [96] has been used. It refers to the detection of a two-dimensional buried object in a half space configuration (under transverse magnetic illumination conditions). The considered crosssectional investigation test domain inv has height = 1 m and width = 2 m. The region inv has been partitioned into 40 × 20 square cells for solving the forward problem (by means of the method of moments) and 30 × 15 subdomains for the inversion process. The lower half space medium has dielectric properties = 4 and = 0.01 S/m. The upper half space is air (modelled as vacuum). A dielectric circular cylinder is buried inside inv . In particular, it has center at r 1 = (0.35, −0.3) m, diameter 1 = 0.16 m, dielectric permittivity = 8, and electric conductivity = 0.02 S/m. A GPR B-scan has been simulated by considering = 30 equally spaced measurement positions on a line at height ℎ = 0.05 m from the soil level, = 2 m long, and horizontally centered at the origin. The offset between the transmitter and the receiver is TX = 0.3 m. An additive white Gaussian noise with zero-mean value and signal-to-noise ratio SNR = 10 dB has been used for corrupting the scattered electric field data. Figure 4 reports the B-scan of the estimated field scattered by the buried cylinder. In the qualitative step, 21 frequency samples equally spaced in the band from 150 MHz to 350 MHz have been used. In the quantitative inversion, a subset of = 5 equally spaced frequencies in the same band has been considered. In the inexact-Newton algorithm, IN = 20 and LW = 10 maximum Gauss-Newton and Landweber iterations, respectively, have been considered. The distribution of the normalized qualitative object function Λ inside the investigation domain is presented in Figure 5. In Figure 6, the quantitative reconstruction of the relative dielectric permittivity in inv is shown. The scientific literature on electromagnetic inverse problem is now very wide and the reader can refer to the mentioned papers and the references therein. A significant interest has also been focused on the use of canonical multilayer objects, since they allow analytical or quasianalytical solutions that can be used for fast computations when dealing with configurations that approximate real imaging cases [116]. In addition, they can be used in the validation phase of inversion procedures. Just as an example, in [117] a two-step inversion procedure for detecting the dielectric permittivity and the velocity of multilayer axially moving targets has been developed, with particular emphasis to cylinders with elliptical cross sections, for which related Green's function has been deduced in [118].

Conclusion
Layered structure modelling is an important topic in in modem applied electromagnetics. A review of methodologies for the EM scattering in layered media has been presented in this paper. Different techniques in the literature have been examined and discussed, by emphasizing their basic principles, advantages, and disadvantages. Particularly, special attention has been given to direct methods for modelling scattering in layered rough structure that have been recently established. The role of stratified media in the solution of electromagnetic inverse scattering problems has been discussed, too. In particular, the current trend in the development of inverse procedures for target detection in half spaces and in multilayer materials has been delineated. We emphasize that the model selection generally represents a compromise between accuracy and computational complexity, with the best compromise that may significantly change, depending on the application contexts.
Therefore, our review has been aimed at providing a current perspective on the subject, also highlighting innovative research directions. However, we make no claim to cover all the many topics pertinent to layered media. Indeed, with the rapid growth of the field, such a task would be almost impossible in a single paper.
The list of publications included in the paper is representative for our discussion and will hopefully prove useful to any researcher active in the area. Nonetheless, it does not exhaust the very extensive literature on the subject, which indeed covers several disciplines (optics, applied mathematics, etc.).