Scattering on rough surfaces with alpha-stable non-Gaussian height distributions

Abstract We study the electromagnetic scattering problem on a random rough surface when the height distribution of the profile belongs to the family of alpha-stable laws. This allows us to model peaks of very large amplitude that are not accounted for by the classical Gaussian scheme. For such probability distributions with infinite variance the usual roughness parameters such as the RMS height, the correlation length or the correlation function are irrelevant. We show, however, that these notions can be extended to the alpha-stable case and introduce a set of adapted roughness parameters that coincide with the classical quantities in the Gaussian case. Then we study the scattering problem on a stationary alpha-stable surface and compute the scattering coefficient under the first-order Kirchhoff and small-slope approximations. An analytical formula is given in the high-frequency limit, which generalizes the well known geometrical optics approximation. Some numerical results are given and discussed.


Introduction
In the theory of scattering by random rough surfaces, a systematic yet crucial assumption is the Gaussianity of the height distribution. This hypothesis is adopted for technical reasons but is often far from being realistic. Gaussian distributions miss essentially two important features of natural surfaces. First, they are unskewed, that is unable to explain an eventual asymmetry between crests and troughs. Second, the tail of these distributions is exponentially decreasing, thereby not allowing the occurrence of 'big' events; Gaussian distributions are therefore also irrelevant for rough surfaces that can exhibit peaks of very large amplitude with respect to the mean level. In this paper, we shall focus on the second aspect only and discard for this first attempt the modellization of skewness.
Surfaces with a non-negligible probability of very large peaks are well described by heavy-tailed height distributions, that is distributions whose tail is polynomially (instead of exponentially) decreasing. The archetype of such distributions is the stable distribution. This family of distributions shares the essential property of being stable under linear combination of independent variables with the Gaussian distribution, and also a generalized central limit theorem for normalized sums of independent identically distributed variables. This is why they are found in the descriptions of many physical and economical systems. To obtain at once a better feeling for how stable distributions can improve the description of rough surfaces with contrasted peaks, consider the sample surfaces plotted in figure 1. Both surfaces have the same first absolute moment (i.e. |h| ) and the same correlation length (these notions will be explained in the following), but have different height distributions: one is Gaussian while the other is alpha stable with α = 1.5. As can be seen, the alpha-stable surface is globally flatter but has much larger peaks than the Gaussian one.
To the best of our knowledge, very few works have been concerned with non-Gaussian surface scattering. Chen and Fung have proposed a numerical simulation method for (some) non-Gaussian profiles with given correlation function [1] and performed a numerical study of the one-dimensional scattering problem for a particular height distribution (the so-called modified exponential distribution) [2]. Tatarskii [3] has developed a specific non-Gaussian surface model that allows us to compute all the multi-dimensional characteristic functions, thereby making scattering computations possible [4]. Jakeman [5] has considered the case of a non-Gaussian random phase screen in optical scattering and computed the scattered intensity for some particular distributions. We take another step in the direction of non-Gaussian surfaces by studying the family of alpha-stable surfaces. In contrast to the previous models, they are not second order (that is they have infinite variance) and thus require a complete redefinition of the roughness parameters. We shall show that, in spite of their apparent complexity, alpha-stable surfaces yield to simple expressions for the scattering computations.
To avoid some confusion about the scope of this work, let us make it clear that we do not seek to improve ocean models, as this has been the main motivation for developing non-Gaussian models in scattering from rough surfaces (e.g. [3,4]). Our aim is to study the impact of a strong non-Gaussianity of the height distribution on the scattering amplitude, and this for an arbitrary correlation structure.
We have chosen to work in the one-dimensional case, which means that the surface is assumed to be invariant in one direction and hence only described by its profile. This restriction has been adopted to make the presentation clearer in introducing new concepts and to avoid some technical difficulties. The extension to the two-dimensional case is, however, straightforward.

The scattering problem
Although the one-dimensional scattering problem is well known, we briefly recall the main results for the convenience of the reader and to fix the notations.

The scattering amplitude
A z-invariant infinite surface y = h(x) separates the vacuum (upper medium) from a perfectly conducting medium. A linearly polarized time-harmonic plane wave with wavevector K 0 = (k 0 , −q 0 ) and wavenumber K = k 2 0 + q 2 0 = 2π/λ is impinging at incidence θ 0 on the top of the surface. The geometry of the problem is depicted in figure 2. The total field F in the upper medium is written F = F i + F s , where F i is the incident field: and F s is the scattered field. Here F stands for the electric field E or magnetic field H according to whether the polarization is s or p. It satisfies the Helmholtz equation in the upper domain, and vanishes in the lower domain, F = 0, y < f (x), with Dirichlet (Neumann) boundary conditions on the interface in the s-(p-) polarization case. The harmonic time dependence e −iωt will always be implicit. The normalization factor q −1/2 0 is set to obtain a constant energy flux (which is proportional to q 0 |F i | 2 , see [6] for details) through a unit horizontal surface. Above the highest excursion of the surface, the scattered field admits a Rayleigh expansion: (2.1) Here and everywhere, k and q represent the horizontal and vertical components of the wavevector K, where the square root defining q is taken to be positive imaginary when k > K (corresponding to evanescent waves). The complex coefficient S(k, k 0 ) is called the scattering amplitude. There is in general no closed analytical formula for the latter, which has to be evaluated numerically. However, under some physical assumptions, the scattering amplitude admits simple analytical approximations. In this paper we shall rely on two widely used approximations, namely the first-order Kirchhoff approximation (KA) and the small-slope approximation (SSA) introduced recently by Voronovich [7,8]. The former is known to be valid when the surface is smooth at the resolution of the electromagnetic wavelength, roughly speaking, while the latter becomes accurate when the mean surface slope goes to zero. In the KA as well as SSA, the scattering amplitude can be written where the optical factor B(k, k 0 ) depends only on the geometry, the polarization and the considered approximation, and the structure factor (k, k 0 ; h) contains the dependence in the profile. Precisely, Note that B(k 0 , k 0 ) = 1 and that |S(k, k 0 )| does not depend on the polarization under the KA.

Finite beams
Although the scattering amplitude is defined in the basis of plane waves, only finite beams can be sent in reality or dealt with numerically. The incident beam then assumes the form where the spectral window functionŴ L = L 1/2Ŵ (kL) is centred about zero. We find it convenient to define this window function as the Fourier transform of a family of functions The function W L is a window function in the space variable and delimits the illuminated region. The parameter L characterizes the typical size of the illuminated patch. Usually the window function is chosen to be Gaussian, The normalization is again chosen so as to obtain a unit incident energy flux E in through any infinite horizontal plane.
Above the surface, the scattered field is again given by a Rayleigh expansion (2.1), defining some other scattering amplitude, say S L (k, k 0 ). Note that the scattering amplitude for an incident beam is the linear superposition of the plane wave scattering amplitude weighted by the spectral window function: This slightly complicates the expression of the structure factor (2.3) since it involves an additional integration in the spectral variable. This, however, can be greatly simplified when the illuminated region becomes large with respect to the electromagnetic wavelength, that is kL 1. In this case we have where the optical factor (2.2) is unchanged and the structure factor is a tapered integral in the space variable, The total scattered energy flux E out through an infinite horizontal plane is given by and the conservation of energy reads E in = E out = 1. The quantity q|S L (k, k 0 )| 2 can be seen as the amount of energy scattered by unit angle. In the statistical case where the surface is a random process and therefore the scattering amplitude a random variable, the ensembleaveraged quantity is called the scattering coefficient and can be measured experimentally. The use of finite beams is not mandatory in this theoretical study. However, all the formulas pertaining to the case of incident plane waves can be directly recovered in the limit KL → ∞, and thus this framework is more general. Also it prevents some mathematical improperness (such as dealing with non-convergent integrals) and makes the numerical checks possible.

The Gaussian picture
Stationary centred Gaussian processes h(x) are completely described by their second-order properties, that is either their correlation function or their so-called power spectrum (or spectral density) Here and everywhere the brackets represent the ensemble average andf the Fourier transform of a function f . When the power spectrum is an absolutely integrable function (as in most cases), the process admits two equivalent representations in the spatial and spectral domains. It can be written as the moving average (MA) of a Brownian motion B, where F is some square integrable function (henceforth referred to as a 'filter') satisfying and dB is white noise, i.e. formally dB(x) dB(y) = δ(x − y). It can also be written as a Fourier integral, where Z is a stationary Gaussian process with independent increments, formally Roughness parameters for Gaussian surfaces. The roughness of Gaussian surfaces is usually described by two parameters. The vertical scale is given by the RMS height σ : while the lateral scale is determined by the correlation length l, namely the shortest distance for which The behaviour of the correlation function at zero determines alone the differentiability or regularity properties of the process. Mean-square differentiable surfaces are those for which the correlation function is twice differentiable at zero or, what amounts to the same, for which the spectrum has a finite second moment: The above expression then defines the root-mean square slope s = h 2 1/2 . In terms of filter, this is also equivalent to the requirement that the function F have a square-integrable derivative F : Such surfaces are termed mono-scale in the literature, due to the high-frequency cut-off of their spectrum. When this condition is not fulfilled, the surface is said to be multi-scale or fractal, because of the slow (i.e. power-like) decrease of the spectrum. Given a correlation function C 2 (r), it is difficult in general to find an analytical expression for the corresponding filter F by deconvolution of (3.6). This is, however, possible in the particular but important case of Gaussian correlation functions: which admits the simple filter (3.10)

The stable case
We refer to the by now classical book [9] for an exhaustive presentation of stable variables and processes. We shall only briefly recall here the notations and properties that are necessary for the following. There are many equivalent definitions of stable variables. The most constructive definition, although the less intuitive, is the following. A random variable X is said to be stable with stability index α, mean µ, skewness β and scale parameter γ (written X ∼ S α (γ, β, µ)) if its characteristic function satisfies The term 'stable' originates from the stability property of such variables under linear combinations. Precisely, if X 1 and X 2 are independent with X i ∼ S α (γ i , β i , µ i ), then X 1 + X 2 ∼ S α (γ, β, µ) with γ α = γ α 1 + γ α 2 , βγ α = β 1 γ α 1 + β 2 γ α 2 and µ = µ 1 + µ 2 . The stability index satisfies 0 < α 2 and the skewness parameter −1 β +1. The case α = 2 corresponds to the Gaussian case. In this case the skewness parameter is irrelevant. In this paper, we shall consider for the sake of simplicity only symmetric stable processes, that is we shall assume β = µ = 0. Following the standard convention, we shall write X ∼ S α S to say that a random variable X follows a symmetric stable distribution. The characteristic function of such a random variable reduces to exp(iuX) = exp(−|u| α γ α ). (3.11) Hence an S α S variable has a probability distribution function Note that an S 2 S variable with scale parameter γ corresponds to a Gaussian variable N(0, 2γ 2 ).
There are no closed analytical formulae for the probability distribution functions of stable variables, except in the particular cases α = 2 (Gaussian distribution), α = 1 (Cauchy distribution) and α = 1/2 (Levy distribution). However, the following asymptotic behaviour can be shown for 0 < α < 2: with c α = 1 π (α) sin(απ/2). Hence stable distributions are heavy tailed with a power-law decay related to their stability index. An essential difference between non-Gaussian stable and Gaussian variables is the occurrence of infinite moments for the former, as can be seen from the previous formula. If X ∼ S α S, then In particular stable variables have infinite variance for α < 2, and infinite absolute mean for α < 1. In this paper we shall assume 1 < α 2 for technical reasons. In this case the first absolute moments are finite and given by (see [9, p 18]) where γ is the scale parameter. A stable process is a process whose finite multi-dimensional distributions are stable (we refer to [9] for the definition of stable vectors since this would bring us too far off the subject). Unlike Gaussian processes, a stable process cannot be determined by its two-point distribution function only. However, one can take advantage of the stability property to construct stable processes by linear summation of stable noise. For this, recall the definition of the S α S Levy motion, which is the stable analogue of Brownian motion. This is a process M(x) with independent increments such that M(x) − M(y) ∼ S α (|x − y| 1/α , 0, 0) and M(0) = 0 almost surely. Then any process of the form where F is some measurable function with |F | α < ∞, is a stationary S α S process, with characteristic function hence with scale parameter The stochastic integral (3.15) is well defined when F is piece-wise constant. In the general case, it is defined by approximating F by such functions and taking the limit in probability. The 1/ √ 2 factor in the MA representation has been set so as to coincide with the Gaussian case (3.15) for α = 2, since in that case 2 −1/2 dM(x) = dB(x) corresponds to standard Brownian motion.

Roughness parameters for stable surfaces. The usual RMS height is no longer relevant for stable non-Gaussian surfaces since the latter have infinite variance. However, a natural extension of the RMS height is the (renormalized) scale parameter
2σ, 0, 0), so σ is actually a good indicator of the vertical scale of the process. Likewise, the correlation function is not defined for α < 2 but can be replaced by its stable analogue for 1 < α 2, the so-called covariation function, which in the case of an MA process (3.15) can be defined as Here a p is the signed power of a, that is a p = sign(a)|a| p . The covariation has in general weaker properties than the covariance; in particular a vanishing covariation does not necessarily imply independence of the corresponding variables. However, in the case of an MA representation (3.15) with non-negative filter, note that C α (r) = 0 implies F (x)F α−1 (x + r) = 0 a.e., which in turn implies the independence of h(x) and h(x + r) (see e.g. [9, p 129]). Note also that which is a generalization of (3.7), and as can be easily seen using the Hölder inequality. Also we have C α (r) → 0 as r → ∞, as follows from lemma A.1(ii). Hence, the covariation function is a good indicator of the dependence of points with respect to their mutual distance. A natural extension of the correlation length is then the shortest distance l for which As in the Gaussian case, the differentiability properties of stable processes depend on the regularity of the filter. For 1 < α 2, the sample paths are almost surely differentiable as soon as F has a derivative F in L α (R) (see the appendix). The differentiated process h is then again an S α S process with scale parameter F α . By analogy with the Gaussian case, we shall identify the quantity F α with the slope parameter. Again, it is in general difficult to find closed forms for the filter F corresponding to a given covariation function. However, for Gaussian covariation functions C α (r) = σ α exp(−r 2 /l 2 ), the filter is easily found to be . (3.20) Remark. An alternative way of defining the vertical scale in the case 1 < α 2 is to use the so-called average roughness |h| (also known in surface metrology as arithmetic average, centre line average or arithmetical mean deviation), that is finite. Note that it is proportional to the scale parameter by (3.14). However, we do not like this solution, since it does not coincide with the usual roughness parameter h 2 in the Gaussian case. Also, it is restricted to the case 1 < α 2.
Numerical generation. Sample paths of stable surfaces can be obtained by approaching the filter by piece-wise continuous functions in (3.15). If the filter has compact support, say [−L/2, +L/2], we may approximate where the M j are i.i.d. S α ((L/N ) 1/α , 0, 0) variables. Note that this is a discrete convolution, which can be performed efficiently by fast Fourier transform. As to the S α S random variables, they can be generated numerically through the use of uniform variables thanks to the transformation recalled in the appendix (proposition A.3). This method has been used to obtain the plots in figure 1.

The scattering coefficient for stable surfaces
We shall now assume that the surface profile h(x) is an S α S process with 1 < α 2 and compute the corresponding scattering coefficient. It will be convenient to consider a family of filter functions derived from a normalized template: where the filter F has unit scale parameter and correlation length, while the corresponding parameters for F σ,l are σ and l. Note that the different slope parameters are related to one another by the formula The process will thus be given by the MA representation (3.15) with filter F σ,l . Under the KA or SSA, the first two moments of the structure factor are given by Now the one-and two-point characteristic functions of the process can be computed explicitly using property (3.16), yielding Note that in the Gaussian case α = 2, this last quantity is none other than the structure function: By lemma A.1 we have σ α cross (x) → 2σ α as |x| → ∞, which suggests introducing the function . Using the convolution theorem and recalling that the two-dimensional Fourier transform of a diagonal function f (x − y) is 2πf (ξ)δ(ξ + η), the second moment can be reduced to a simple integral, namely In the plane wave limit KL → ∞ we haveŴ L (k) → 0 and |Ŵ L (k)| 2 → δ(k). This yields the following expression for the scattering coefficient: This formula allows for numerical estimation but in general the analytical computations cannot be pushed further, even in the Gaussian case α = 2. However, a good analytical approximation can be derived under some additional physical assumption. For |x| l we have σ α cross (x) 2σ α while for |x| l we may approximate Since σ α l −α |x| α σ α ( σ α ) whenever |x| l ( |x| l), we have the global approximation This entails the following approximate formula for the scattering coefficient: where p α is the S α S probability distribution function (3.12) and is the slope parameter. For α = 2 we recover the well known formula (e.g. [6, p 127]) which is often referred to as the geometrical optics (GO) approximation of the scattering coefficient, as it corresponds to the high-frequency limit. The error committed on approximating ϕ(x; q + q 0 ) is of the order exp{−[(q + q 0 )/2] α σ α }. Thus, the above approximation becomes reliable for a large Rayleigh parameter q 0 σ , as in the Gaussian case (see e.g. [10] for a complete review in the Gaussian case).
The last results require some comments. First note that in the GO approximation (4.25) the scattering coefficient depends only on the slope parameter s, as in the Gaussian case. Furthermore the scattering diagram (i.e. the intensity I as a function of the scattering angle θ , for fixed incidence θ 0 ) has the following remarkable property: its shape is given by the height distribution function p α (providing the change of representation θ → (k − k 0 )/(q + q 0 )). This can be interesting for an estimation procedure in the inverse problem (that is to estimate α from the scattering coefficient, or to determine whether the surface is alpha stable or not) but will not be discussed here. Figures 3 and 4 present some numerical simulation. The scattering coefficient has been plotted in decibels (i.e. 10 log I (k, k 0 )), according to the usual convention. The diagram in figure 3 shows how the scattering pattern is affected by the stability index α for a given slope parameter s and a fixed correlation length l in the case of a Gaussian filter. Note that keeping a constant slope and correlation length for varying αs also imposes varying the scale parameter σ in view of (4.26). In the stable case, the curve is narrower about the specular direction but has a slower decay than in the Gaussian case at large angles.
For Gaussian surfaces the degree of roughness is generally quantified in terms of RMS height or average roughness. In the non-Gaussian case it is, however, dangerous to use the moments as an indication of roughness, as far as the scattering properties are concerned. This is illustrated in figure 4, where surfaces with the same Gaussian filter, the same correlation length l = 3 and same average roughness |h| = 0.5/ √ 2 have been considered for varying αs (hence also varying σ s according to (3.14)). Two of these sample surfaces have been plotted in figure 1. The scattering pattern concentrates more and more about the specular direction as α becomes smaller. This means that the surface appears less rough for smaller αs, even though the average roughness is kept constant. This is easily understandable: for smaller values of α, a given level of average roughness is reached thanks to the contribution of a small number of large peaks in a flat environment. Note also that the scattered intensity decreases more slowly at grazing incidence in the case of stable surfaces. Again, this is quite natural: large and sparse peaks can be seen as isolated antennas that radiate cylindrical waves in all directions of space.

Conclusion
In this paper we have studied the scattering problem on a one-dimensional rough surface modelled by a stationary stable process. We have shown that the usual roughness parameters that are used for Gaussian surfaces can be extended to the larger family of stable surfaces and are a particular case thereof. Finally, we have derived analytical approximations for the scattering coefficient in the framework of first-order approximations (KA and SSA). This is only the beginning of a fruitful field of investigation: important questions remain such as finding an estimation procedure for the stability index in the inverse problem. Also, this work has to be extended to the two-dimensional case and the model to be confronted with real data.
Lemma A.1. Let F be a function in L α (R) with 1 < α 2. Then Proof. Denote K R = [−R, +R]. By the Hölder inequality with conjugate powers α and α/(α − 1) we have Since |F | α is integrable, one integral involved in each product goes to zero while the other goes to a constant, proving (i). The proof of (ii) is similar using the Cauchy-Schwarz inequality.