Simulation of subsurface imaging for remote sensing and buried object detection from airborne platforms

ABSTRACT This paper addresses the problem of efficient simulations of the signal backscattered from buried structures. The proposed modelling of subsurface soil structure is based on a geologically – and hydrologically-motivated approach by using of real geological database. Six different Soil Typological Units have been used to model the inhomogeneous and geological soil structures as well as for moisture profile simulation. We use the semi-empirical pedotransfer function to calculate the Complex Relative Effective Permittivity in frequency range of 0.3–1.3 GHz. The subsurface penetration capabilities have been analysed with the regard to attenuation along the propagation path in a layered structure. We use the backward-propagation matrix and ray-tracing method to compute the backscattered signals from five targets located at five different depths in order to analyse the defocusing effects due to EM wave propagation. The computation was performed for three different radar systems with various bandwidths and carrier frequencies.


Introduction
Subsurface remote sensing, soil parameter estimation and buried object detection are a constantly growing field of interest in the remote sensing domain. Nowadays, the ground penetrating radar (GPR) (Jol, 2008) is the leading radar technology for non-destructive microwave imaging of the subsurface structure. Several imaging methods are applied for buried object imaging and localization, in particular, the synthetic aperture radar (SAR) technique (Gilmore, Jeffrey, & LoVetri, 2006). Both GPR and SAR imaging processes are formed using antenna movements. The essential attribute of the airborne ground penetrating radar (AGPR) system lies in combining the advantages and key features of these two mapping techniques for imaging large areas and adequately penetrating the soil with the scope of providing high resolution imaging capabilities. Some early attempts of AGPR have already been performed with the shuttle imaging radar (SIR) imaging instruments in the 80s in order to reveal the subsurface valleys and sand-buried landscape as well as for snow mapping. The French agency ONERA has developed the RAMSES, which has been used for L-band subsurface imaging in south-central Egypt (Paillou et al., 2003). In these AGPR pioneering experiments, the estimated penetration depth varies from 1.5m down to several meters which is promising for subsurface sensing.
However, these AGPR experiments have been carried out over arid areas, the structures of soils/snow are strongly similar to a homogeneous structure (dry sand, snow or ice) with almost constant EM properties. A similar homogeneous approach of the soil geological structure was used for numerical simulations (Fu, Liu, & Liu, 2012;Liu & Feng, 2011;Martinez-Lorenzo, Rappaport, & Quivira, 2011;Zhai, Zhang, & Zhang, 2006). Unfortunately, such homogeneous soil structures do not correspond to the majority of lands, which are a mixture of different soil types, arranged in horizontallayered structures with varying texture and content inside a geological soil horizon. Due to the inhomogeneous, dispersive and lossy nature of the soil, the backscattered signal derived from multi-reflections of the electromagnetic waveform is disturbed by dispersion and lossy effects. This leads to energy dissipation and received pulse distortion that strongly degrad the range (depth) resolution. In order to fully understand the defocalization mechanisms and the possible refocussing algorithms (e.g. time reversal), numerical simulations provide a better insight. Many methods have been applied for numerical computations and simulations of the signal backscattered from buried structures (Fu et al., 2012;Le Caillec, Redadaa, Sintes, Solaiman, & Benslama, 2011;Liu & Feng, 2011;Martinez-Lorenzo et al., 2011;Zhai et al., 2006). Some experiments were performed in laboratory-scale model (Elsherbini & Sarabandi, 2010;Fu, Liu, Liu, & Lei, 2014).
The two main scopes of this paper follows: • We propose a realistic model of geologic, hydrologic and electromagnetic parameters for simulating the EM propagation into the soil. Based on six Soil Topological Units (STUs) chosen from the CONTACT J-M Le Caillec jm.lecaillec@imt-atlantique.fr Institut Mines Télécom -Télécom Bretagne, CNRS UMR 6285 LabSTICC, Plouzané, France SPADE-2 European project from the European Soil Data Centre (namely BE 320,202,FI 358,005,GM 490,149,HU 360,084,IT 390,925,PO 3,511,382). • For such a soil structure, the moisture profile and electromagnetic model are calculated. Then, we use the ray-tracing method (Spencer & Murty, 1962) to compute the signal backscattered from point-like target located under the surface. The backward-propagation matrix is calculated for highly inhomogeneous structures consisting of layers with random geological texture and moisture content, the final aim of this paper being to study the defocusing effects of the SAR pulse.
This paper is organized as follows. The second section describes the model of backscattered signal for the AGPR systems and the third models the geologic and hydrologic parametres. The fourth section discusses the computed results for the six structures and compares the backscattered signals for the three AGPR systems for the most inhomogeneous structure, the fifth section section concluding the paper.

Airborne ground penetrating radar
From the geometrical point of view, the AGPR imaging configuration is based on the Side Looking Airborne Radar (SLAR) configuration ( Figure 1), with parameters of conventional airborne SARs (Table 1). While airborne SAR systems are devoted to imaging the surface with either man-made or natural objects, the AGPR intends to image the optically invisible soil structure located under the surface. We assume a linearly polarized (TE or TM) plane wave (radar pulse) as Figure 1. Reflection and transmission of the obliquely incident EM wave in stratified isotropic medium made up from N z homogeneous, planar and thin layers.
where A 0 is the constant amplitude, τ the time duration and t the time, f 0 the initial frequency, f ¼ f 0 þ α f t the instant frequency and ω ¼ 2πf its corresponding pulsation, ϕ 0 an initial phase, α f ¼ AE B f =τ the LFM slope, B f the bandwidth and f c ¼ f 0 þ α f τ=2, the central (carrier) frequency. Table 1 presents the radar pulse parameters of three SAR systems simulated in Section 4. For all radars, we have assumed τ ¼ 1:00μs and a radar signal peak power P 0 ¼ 1kW, (P 0 ¼ 8kW, for RAMSES Dubois-Fernandez et al., 2002). Moreover, we set the platform altitude H p ¼ 4000m and the incident angle θ i ¼ 40°( at the swath centre) for all the considered radars. The theoretical slant range resolution in free space In a rich scattering medium, as the subsurface, the dominant reflection echo comes from the boundary between the radar half-space and target half-space. The response from the soil structure contains direct amplitude-scaled, time-delayed and distorted echo from point scatterers, superposition of indirect multipath reflections from the point scatterers as well as of an infinite number of multipath reflections from the soil inhomogeneity. Hence, the target response s r t ð Þ is expressed as (Kedzierawski, Le Caillec, & Czarnecki, 2014) where σ c is surface reflection coefficient proportional to the surface clutter area A c , i.e., σ c ¼ A c σ 0 , τ d is the two-way time delay between radar and surface, h i t ð Þ is impulse response of ith propagation path inside the soil, t i is propagation time of ith path and * denotes time convolution. Empirical model of σ 0 can be found in Oh, Sarabandi, and Ulaby (1992). n t ð Þis the additive antenna noise with zero mean, the power of which, denoted P N , depends on the bandwidth.

Subsurface modelling
From the geological point of view, the soil is made up of the mineral particles, organic matter, water and air (Behari, 2005; Food and Agriculture Organization of the United Nations, 1985). Soil composition is a key point to derive the dielectric properties for realistic simulations of the backscattered signal. Geological databases, such as the real geological database SPADE2, are a valuable source of information about the soil textural composition and inhomogeneity (Kedzierawski, Le Caillec, Czarnecki, & Pasternak, 2012). In many works, for instance (Kedzierawski et al., 2012;Kuo & Moghaddam, 2007;Lambot, Slob, van den Bosch, Stockbroeckx, & Vanclooste, 2004;Le Caillec et al., 2011), the soil has been modelled as homogeneous and planar layers with different heights, inner textural composition and thus EM properties.
The problem of defining an accurate model for the EM wave propagation lies in the proper choice of the number and the height of the soil layers. The moisture profile is also dependent on the type of soil, and its shape results from the long-time internal diffusion processes. Therefore, the natural phenomenon of soil inhomogeneity has to be taken into account in order to derive the spatial distribution of the dielectric constants.

Geological layer
Files from the SPADE-2 database contain the content of sand s d clay c l , silt s t and the bulk density at non uniform depths (named horizons). In order to perform the EM propagation simulations, we first define the z depth vector (of random step), for which the soil texture is approximated (see below). It starts at z 0 ¼ 0:0m and is iteratively sampled as follows: where n is the layer number. Total number of layers is N Z ¼ 1509 and maximum considered depth of the soil is Z max ¼ 15:0317m. The individual soil layer l n is modelled as a geologically uniform, thin and planar layer located between two following interfaces. The last (absorbing) soil layer l N Z starts at the z max; (see Figure 1). In order to convey the inhomogeneity within a geological horizon, we assume that h n layer height has a Gaussian distribution of mean μ h ¼ 0:01m and standard deviation σ h % 0:0017m. This value of μ h is tradeoff between realistic simulations and the computational load. In fact, as seen in Section 3.2, the Richards' equation used for the soil moisture simulation, is a highly nonlinear differential equation and leads to strong soil moisture profile variations in the first few centimetres of the subsurface (as shown in Section 4.2). The value of σ h allows some randomness of the thickness but avoid negative values of h n for the considered number of layers. An additional refinement for the EM wave propagation could be to take into account the roughness between two layers (instead of planar interfaces), but except for the surface, such roughness models, which depend on the composition of the two layers, are not available.
Since c l z ð Þ þ s t z ð Þ þ s d z ð Þ ¼ 100%; , the c l z ð Þ clay and s d z ð Þ sand content are considered as two independent variables while s t z ð Þ silt content is calculated as fulfilment to 100% at z n .
Clay c l z ð Þ and sand s d z ð Þ contents are given by Gaussian distributions, where μ c (μ s ) is a mean value of clay (sand) content in the horizon containing the n th layer (values extracted from the database file) and σ c (σ s ) is the standard deviation of clay c l z ð Þ (sand Additionally, the STU contains information about the ρ b bulk density average of the dry soil in each horizon. However, due to the approximation of soil texture at z n , ρ b z ð Þ and p z ð Þ, the porosity, are calculated based on clayc l z ð Þ and sand s d z ð Þ content as follows (Behari, 2005): where ρ s is the density of the rock in the soil. These equations are used to simulate the water content profile (next section) and EM propagation into the soil (Section 3.3).

Moisture content and profile
Besides the mineral particles, the soil is mostly made up of water and air. Analytical methods can be applied to estimate the moisture profile Θ v in one-dimension (Varado, Brauda, Rossd, & Haverkampa, 2006). In fact, Θ v is the water flow in porous media, usually described by the Richards' equation (Richards, 1931): Þis a source function accounting for water exchange process at the soil-air interface such as precipitation, evaporation and run-off as discussed again below. K t; r ð Þ is the hydraulic conductivity and h t; r ð Þ is soil water pressure head. Furthermore, both K t; r ð Þ and h t; r ð Þ are highly non-linear functions of the soil moisture, and are expressed as Þ is soil water tension and K s , Θ s Ψ s and b are the saturated hydraulic conductivity, saturated water content, saturated soil water tension and curve-filtering (shape) parameter. All these parameters depend on the soil types (textural class) (see Cosby, Hornberger, Clapp, & Ginn, 1984, Table 3). The vertical water flow (thus r ¼ z in the previous equations) is calculated at the z n points, for which the soil texture has been approximated. Three additional features have to be taken into account in the simulation, i.e., rainfall, evaporation coefficient and initial moisture profile. The influence of both rainfall and evaporation is included in the source function F t; Θ v t; 0 ð Þ ð Þ , which determines the soil non-homogeneity in the first layers. Depending from the soil textural classes, the maximum amount of water theoretically possible to store differ significantly, and is limited by Θ s the saturated water content or p z ð Þ porosity (see Equation (4)).

Dielectric properties
Besides the difference in the number of factors taken into account, models for deriving the dielectric properties are mostly based on the pedotransfer approach allowing the calculation of EM properties as a function of soil textural composition c l z ð Þ, s d z ð Þ, ρ b z ð Þ the bulk density, p z ð Þ the porosity, Θ v z ð Þ the moisture profile, f the frequency and T c the temperature (in this section the z dependence is omitted for the sake of notation clarity). Several models or mixing equations of EM properties have been proposed for soil-water mixture, for instance (Behari, 2005;Hallikainen, Ulaby, Dobson, El-Rayes, & Lil-Kun, 1985;Peplinski, Ulaby, & Dobson, 1995;Wang, 1980). These models give both real and imaginary parts of the complex relative permittivity as a function of frequency and water content. for dielectric constant modelling, since it allows the calculation of these constants in a wide GPR frequency range, i.e., the whole UHF band and the lower part of L band assumed for AGPR (Table 1).
The common features of the aforementioned models are, first, the permeability is assumed the same as in free space, i.e., the relative permeability equals 1 and second the soil is characterized as a lossy material, the Complex Relative Effective Permittivity (CREP) is expressed as ε r ¼ ε 0 r À jε 00 r . ε 0 r is the relative dielectric constant of an equivalent lossless soil and ε 00 r is the relative dielectric loss factor related to energy absorption by the moist soil (Behari, 2005). In the Peplinski's model, they are given by where ρ s ¼ 2:66 g cm À Á 3 is the specific density of the solid soil particles (assumed constant for all soil types), ζ P ¼ 0:65. β 0 ¼ 1:2748 À 0:519s d À 0:152c l and β 00 ¼ 1:33797 À 0:603s d À 0:166c l are empirically determined soil-type dependent constants, and finally ε s is the dielectric constant of soil particles (rocks), is given by the empirical model ε s ¼ 1:01 þ 0:44ρ b ð Þ 2 À 0:062 The water mixture complex permittivity (common to the other models) is given by where ε w0 is the low-frequency limit of permittivity for water, ε w1 is the high-frequency limit of ε 0 fw , τ w is relaxation time of free water, σ w is the water conductivity. According to Wang (1980), ε w1 varies with water profile and porosity as follows: where ε rock % ε s and ε a ¼ 1 are the dielectric constants of rocks and air. The temperature dependence of soilwater mixture is given by the temperature relationship of ε w0 and τ w as follows (Behari, 2005): þ 1:07510 À5 T 3 c 2πτ w ¼ 1:1109:10 À10 À 3:82410 À12 T c þ 6:93810 À14 T 2 c þ 5:09610 À16 T 3 c and at room temperature T c ¼ 20 C, ε w0 % 80.1 and τ w % 9:276410 À12 s. The Peplinski's model, σ w takes into account the porosity as σ w ¼ σ eff :p=Θ v with σ eff ¼ 0:0467 þ 0:2204ρ b À 0:4111s d þ 0:6614c l leading to the static conductivity σ dc ¼ σ eff :p:Θ β 00 ζ P À1 v . This model clearly shows the dependence of ε r on the soil texture, moisture profile, frequency and temperature as well as the dependence of σ dc on soil texture and moisture profile only. Since, both soil texture and moisture profile can be expressed as a function of depth, i.e.,c l z ð Þ, s d z ð Þ, ρ b z ð Þ and Θ v z ð Þ (Equations 4 and 5), these quantities depend also on depth. Gathering all the features, the CREP can be written as and its space-time discretized version is used for numerical simulations EM field in multi-layered medium In layer n, of height h n , the local relative permittivity and permability are denoted ε n ¼ ε 0 ε r;n μ n . We denote τ n the EM travelling time within this layer. Since the soil layer l n is modelled as thin, planar and homogeneous ( Figure 1) the electric field in it can be written as (Kong, 1986) E n ¼ A n e Àjk z n Áz þ C n e jk z n Áz À Á e jk x n Áx where A n is magnitude of the field crossing downward and C n is magnitude of fields crossing upward, k x n and k z n , are the wave number components along x and z axis. Assuming a TE plane wave, which is obliquely incident at soil surface, and the field above the surface is E a ¼ E 0 e Àjk z a Áz þR a;0 E 0 e Àjk z a Áz À Á e jk x a Áx whereR a;0 represents the generalized reflection coefficient. In turn, the total field in the last layer l N z (see Figure 1) becomes Nz :x C N Z ¼ 0 since this last layer is assumed to be absorbing (Kong, 1986). With, the boundary conditions of continuity of the tangential electric and magnetic fields at each interface, the magnitudes of upwards and downwards wave components between layers l n and l nþ1 , i.e., A n , C n , A nþ1 , C nþ1 are coupled as follows: A n e Àjk z n :z n C n e jk z n :z n The EM field within each layer can be thus expressed through the propagation matrices and boundary conditions (much deeper than the location of the object as seen in Section 4) as follows: whereR andT are two unknowns. B n;nþ1 ð Þ is derived from the dispersion relation linking (k x n ,k z n ). In fact, using the tangential field continuity (phase matching conditions) at the interface of two layers (Kong, 1986), we have: where k i is the magnitude of wave vector of incident wave in the air (see Figure 1). Then, the k z n can be deduced from the equation of dispersion We observe that k z n (unlike k x n ) is complex, its real part gives the direction of propagation, in particular the local incident angle θ i;n and its imaginary part the EM losses. This computational formalism can be easily applied for TM polarization by replacing μ n by ε n in backward-propagation matrix (Equation (6)). After identifyingR andT by solving Equation (6), A n and C n can be calculated and the EM wave attenuation/ time of travel can be derived in each layer.
In a GPR approach, for a lossy medium, a usual quantity, denoted δ 1=e , is the depth of penetration (Behari, 2005;Richards, 2009;Zhang & Li, 2008) or skin depth that is the depth over which the wave amplitude decays to 1=e % 0:369 (−8.69 dB) of its value just under the interface. For homogeneous and loss half-space, δ 1=e increases with wavelength (the lower frequency, the higher penetration depth). This useful concept of δ 1=e , shows the effect of soil EM losses (determined by the soil composition and moisture profile). Shortest wavelengths and driest soils provide the greatest depths of propagation but, in practice, an exact value of δ 1=e ; is difficult to compute or estimate.
From a SAR point of view, the attenuation of the penetrating wave (using A n ) at depth z can be calculated, conveying the cumulative losses (penetration losses at the surface and EM losses) and is denoted L α z ð Þ. While δ 1=e is defined with respect to immediate sub-surface value, the transmission losses at airsoil interface are taken into account in L α z ð Þ and thus the difference of dielectric constants between air-soil as well as the incident angle interfere in this quantity.

Soil parameter calculation
The soil texture approximation at z n for the six chosen STU profiles is shown in Figure 2. For these structures, the soil type frequently changes, in the first layers down to 0.55 m, between sandy loam and sandy clay loam, and later from 1.2 m between loamy sand and sand, and from 1.4 m between clay and clay loam, respectively. Comparing the six approximated STU-based soil structures, HU 360,084 contains six different soil types. Although, FI 3,580,005 has also strong variations of structure, but it contains only three different soil types. On the other hand, IT 395025 has regular structure and has also three different soil types, the same as PO 3511382. BE 320202 and GM 490143 jointly contain the same two types of soil, i.e., silt loam in the upper part and loam in the lower part of the structure.

Soil moisture simulation
The soil moisture simulation period was T r ¼ 4 days, with simulation step ΔT r ¼ 1s, and it was assumed as a rainless period occurring after very long rainfall. The evaporation coefficient was simply chosen as proportional to water content at the surface. Results of Θ v z ð Þ soil moisture profile simulation for six chosen STU profiles are shown in Figure 3.
As it can be seen in this figure, for all STU profiles the initial moisture profiles are smaller than 50% and convey the inhomogeneous nature of the soil structure. Due to the surface evaporation and infiltration into the soil (even below 15 m) the total water content decreases with time according to the hydrological properties as, K s ; Θ s , Ψ s , and b (Cosby et al., 1984). The numerically calculated Θ v z ð Þ, the soil moisture profiles are plotted with solid lines. Only, HU 360084 exhibits curves in the initial layers while FI 3580005 and PO 3511382 are slightly curved and other profiles are fairly smoothed. BE 320202 and GM 490143 have almost identical moisture profiles. Except HU 360084, other profiles tend to close water content profiles. Infiltration mostly depends on the hydraulic properties, especially on K s , that describes the speed with which the water can penetrate through porous soil.

Dielectric constant calculation
The CREPs for the six STU profiles based on Peplinski et al. model, are shown in Figure 4 for ε 0 r z; f ð Þ and in Figure 5 for ε 00 r z; f ð Þ. A focus on the variations these two constants (with σ dc z ð Þ) with depth at 800 MHz is given in Figure 6 while these constants are depicted in the first layer over the 0. Comparing all the plots of these figures, the difference between free space and first layer are evident, in the case of ε 0 r z; f ð Þ (Figure 4). The real and imaginary parts of the CREP in the first layer are marked at each plots (see red dots) for 0.8 GHz. The frequency bandwidth of the three systems of Table 1 is also marked by vertical lines. The value of the CREP in the first layer for 0.3 GHz, 0.8 GHz and 1.3 GHz are presented in Table 2 for the STU profiles. As seen in this table, ε 00 r has slightly higher variations than ε 0 r in 0.3-1.3 GHz range. As stated before, the shape of both parts of ε r reproduce the shape of Θ v z ð Þ. Further, ε 0 r is more impacted by the textural composition than ε 00 r , since the soil texture only affects ε 00 r through β 00 (see Section 3.3), while ε 0 r is impacted through β 0 and ρ b z ð Þ bulk density.
When Θ v z ð Þ the moisture profile ( Figure 3) is a monotonically increasing function of depth, then ε r shows piecewise nature of length associated with the height of consecutive soil horizons (Figures 5 and 6).
Although, a high reflection coefficient is desired for conventional surface imaging systems, while for subsurface imaging a high value results in a less amount of energy penetrating into the soil. Thus, the slight change of ε r between air and surface (low EM reflection coefficient), as well as between near-surface layers, results in a high EM transmission coefficient of energy. σ dc z ð Þ 2πfε 0 strongly varies in the considered frequency range, compared to ε 00 r z; f ð Þ. The static conductivity σ dc (Figure 6(c)) is the main factor of EM losses and is associated with sand s d z ð Þ and clay c l z ð Þ content as seen in Section 3.1. For instance, FI 3580005 has the dominant s d z ð Þ content in all horizons, thus, it has the lowest σ dc . On the other hand, IT 390925 has clay for dominant content, hence its σ dc is the highest.
Attenuation and penetration capabilities L α z ð Þ is depicted in Figure 8 for the STU-based soil structures only for 40°and both polarizations. The estimated penetration depths, δ 1=e , are presented in Table 3 and are marked (in thin vertical lines) for various carrier frequencies in Figure 8. In the simulations, θ i was equal to 20°, 40°and 60°, respectively, for the three airborne SAR systems (Table 1). For 20°and 60°, L α z ð Þ has the same shape as for 40°, but is shifted due to different transmission loss at air-soil interface, (small decrease for 60°and relatively high increase for 20°and L α z ð Þ is not depicted for these values). δ 1=e does not exceed few centimetres below soil surface limit marked by dashed horizontal line. The attenuation is nearly linear through soil horizon and it changes significantly at interface between the soil horizons. The deflection of curves occurred at the boundary between soil horizons, in particular, when the soil composition strongly changes between two different types. Since ε 00 r þ σ dc z ð Þ 2πfε 0 decreases with frequency, and then the transmission losses are smaller for Radar 4 than for Radar 3 and Radar 2, the attenuation in the lowest layers is significantly higher for Radar 4. Hence, although L α z 0À ð Þ is lower, just below the surface for Radar 4 than for Radar 2 and 3, at the top of second soil layer L α z 1À ð Þ, is higher for Radar4 than for other SARs. An analysis of the total attenuation along the propagation path in the soil has to be equated to the wave attenuation caused by vegetation or tree canopy (stalks or trunks, and leafs, see Figures 11-19, Figures 11-53 Ulaby, 2013), with the difference that soil surface is a uniform interface between air and soil structure, whereas vegetation or tree crown are not uniform and cohesive surfaces.
The estimated values of δ 1=e differs for various carrier frequencies, are similar for Radar2 and Radar3 and lower for Radar4. For the sake of comparison we also calculate Δδ 1=e the difference of skin depth between radar 2 and Radar 4. In addition, throughout the considered frequency bandwidth, Δδ 1=e penetration depth deviation takes negative values as previously stated. Those negative values of Δδ 1=e show that higher frequency are attenuated strongly. As results from Table 3, the values of δ 1=e and Δδ 1=e depend on soil composition in consecutive layers. However, Δδ 1=e for BE 320202 and HU 360084 as well for GM 490143 and IT 390925 takes similar values, and is smaller than 10 cm or 7 cm, respectively. Hence, small deviations of the penetration depths can possibly occur in the highly lossy soil, even if the penetration depth is much shallower than in the "'easily penetrating'" soils, where Δδ 1=e may be an order of magnitude greater (compare maximum and minimum of δ 1=e and Δδ 1=e for FI 3580005 and IT 390925).

Backscattered signal
In this section, we consider the EM response of a buried point scatterer. Assuming an isotropic point scatterer at depth d n , (see Figure 1) where 0 < d n ( Z max , the upward backscattering component of the electric field can be obtained by changing the wavenumber sign (Le Caillec et al., 2011) when reflecting on the target at d n .
The ray-tracing method allows one to easily understand the distortion, i.e. irregular attenuation and dispersion effect, leading to envelope fluctuations of the (direct) backscattered signal. However, the ray tracing does not take into account the multipath propagation, but it is a first basis to validate more complex methods such as finite difference methods However, the advantage of this method is its ability to easily calculate both EM amplitude and propagation time τ n of the radar chirp through layer l n . Hence, based on τ n the two-way delay time, denoted t d , inside the soil is derived by summing the τ n .
Thus, if we discard the multipath backscattering, then Equation (1) can be rewritten only with time- delayed surface reflection and direct backscattered, amplitude-scaled, time-delay and distorted echo from point scatterers buried in the soil at a specified depth, as follows The two first terms in Equation (7) can be considered separately without lost of accuracy for the backscattered signal envelope analysis, since in the next section we focus on the target imaging.

Defocussing effects
The focussing, in the range domain, is performed by convoluting the received pulse (Equation (7)) with the emitted pulse (Equation (1)) leading to a cardinal sine under the assumption that the received pulse is the weighted/delayed emitted pulse. The defocusing problems, i.e., irregular attenuation, pulse distortion, envelope modulation and energy misplacement, are presented in regard to one individual probing over an identical soil structure (i.e. HU 360084). The AGPR simulations, in order to image a point scatterer at depth d n , have been carried out taking into account the entire propagation path (radarairsubsurface structureairradar) resulting from mission geometry, radar parameters and calculated EM properties of the soil structure ( Figure 1). The point scatterers lie at different depths in order to take into account the pulse distortion with depth. For each target the surface refraction point X TSR (see Figure 1, the TSR abbreviates target surface refraction) is different. The wave travel time in the soil (t d in Equation (7)) strongly depends on its propagation velocity inside soil (i.e., thus on ε r;n see Section 3.3) and propagation angle. Strongly inhomogeneous media induce curved propagation paths and rapid phase changes, which are highly dependent on the incident angle θ i and on the difference of EM constants between successive thin homogeneous layers. Figure 9 shows the normalized envelopes of the backscattered signals. Pulse distortion depends on the centre (carrier) frequency f c , bandwidth B f , and polarization as expected. The irregular modulation of the envelope occurs along the backscattered signal duration, and it increases when the signal bandwidth increases (as seen in Equation (1) the emitted pulse envelop is constant). Obviously, the larger the signal bandwidth, the stronger the irregularities of the pulse envelope. However, although the bandwidth of the Radar 3 is slightly narrower than that of Radar 4, the envelope of the Radar 3 backscattered signal is much more deformed. As already observed in Le Caillec et al. (2011), the lowest frequencies, although having higher penetration capabilities, have not necessary the highest effective bandwidth (that is a frequency range over which the pulse envelop is almost constant) and the highest focussing capabilities.
As previously stated, attenuation increases with the carrier frequency. Table 4 presents the average level of backscattered power in relation to its average level at soil surface (P v for vertical and P H for horizontal). While attenuations in free space are similar due to slightly different slant range to X TSR , the two-way attenuation in soil irregularly increases with the target depth. T D the effective time duration of backscattered pulse is slightly different than the initial pulse (see Table 4 and Figure 9). In fact, the backscattered pulse is longer than the emitted radar pulse duration, and elongation increases non-linearly for the deepest targets. However, similarly to the envelope modulation, the relative elongation of the pulse is longer for the widest bandwidths centralized around the lowest carrier frequencies.
The corresponding normalized outputs of matched filtering are shown in Figure 10. Depending on the distortion level of the backscattered signal, the main lobe of the matched filtering output is wider, due to the mismatch between backscattered pulse and replica of transmitted signal. Additionally, the maximum of output is shifted proportionally to the changes of the effective time duration of backscattered pulse. Since for the shallowest scatterer, the output of the matched filtering, except a slight shift, is almost unaffected. The outputs for deepest scatterers are enriched with irregular side lobes and thus ghost targets appear. As seen in this simulation, the target response is highly dependent on the soil composition. Then the (re)-focussing of buried target requires an adaptive waveform approach such as time reversal for instance.

Conclusions
This paper addresses the EM wave propagation in a highly inhomogeneous, dispersive and lossy medium, as the subsurface. We detail the EM propagation within this (layered) subsurface. For this, a geologicallyand hydrologically-motivated approach is a more realistic, when the soil cannot  be assumed to be an arid or semi-arid area. Numerical calculations show that moisture profile is strongly dependent on the soil type as well as its distribution along the depth. The dispersion changes the effective duration time of the pulse, energy dissipation and pulse distortion strongly affecting the output of the matched filtering and significantly degrading the spatial (range/depth) Figure 8. Total one-way total attenuation L α z ð Þ for various carrier frequency and two polarizations for the STU-based soil structures, θ i ¼ 40 . resolution. We calculate the envelope modulation of the backscattered signal that leads to ghost targets, which introduce ambiguity in the evaluation of the scatterer depth (location). The correct assessment of the true target location (both in depth dimension as well in ground plane) is a major problem, in particular, for subsurface imaging at oblique incidence, where the relative depth misplacement can be relatively high compared to the true target location.