Electromagnetic scattering on fractional Brownian surfaces and estimation of the Hurst exponent

Fractional Brownian motion is known to be a realistic model for many natural rough surfaces. It is defined by means of a single parameter, the Hurst exponent, which determines the fractal characteristics of the surface. We propose a method to estimate the Hurst exponent of a fractional Brownian profile from the electromagnetic scattering data. The method is developed in the framework of three usual approximations, with different domains of validity: the Kirchhoff approximation, the small-slope approximation of Voronovitch and the smallperturbation method. A universal power-law dependence upon the incident wavenumber is shown to hold for the scattered far-field intensity, irrespective of the considered approximation and the polarization, with a common scaling exponent trivially related to the Hurst exponent. This leads naturally to an estimator of the latter based on a log–log regression of the far-field intensity at fixed scattering angle. We discuss the performance of this estimator and propose an improved version by allowing the scattering angle to vary. The theoretical performance of these estimators is then checked by numerical simulations. Finally, we present a rigorous numerical computation of the scattered intensity in the resonance domain, where none of the aforementioned approximations applies. The numerical results show the persistence of a power-law behaviour, but with a different and still non-trivial exponent.


Introduction
Fractal models are known to be of particular relevance in the description of natural rough surfaces in spite of their simplicity.In view of remote-sensing applications, it is therefore important to understand how fractal characteristics of such surfaces impact on wave scattering.Fractal models have the ability of describing complicated and chaotic media by means of some statistical parameters, essentially fractal dimensions.This leads to a new class of inverse problems: instead of reconstructing the exact shape of the surface from the scattering data, it is sought to recover its fractal dimensions only.
The most widely used model of fractal surfaces in the literature is based on bandlimited Weierstrass functions and is deterministic.It has the advantage of providing an analytical expression involving a few parameters together with explicit formulae for the fractal dimensions.Several approaches have been developed to compute the far-field amplitude scattered by this type of profile: the Kirchhoff approximation (KA) in [2,3,5,8,14], the extended boundary condition method of Waterman in [21,22] and the rigorous numerical method in [19].Some qualitative results relating the scattering amplitude to fractal dimensions of the surface have been exhibited but fail to give a precise and general way of computing the latter.Recently, the present authors proposed a method to compute at least one fractal dimension (the so-called correlation dimension) of a deterministic rough surface by the sole knowledge of the far-field intensity [13].The technique relies mainly on the small-perturbation method (SPM) and its simple connection to Fourier analysis.The method was tested numerically on Weierstrass profiles.
The Weierstrass model is convenient for the determination of the scattered field but is far from being realistic.A more appropriate model is given by the class of 1/f processes.These random processes are self-similar in the sense that their spectra exhibit a power law, at least in some wide range of frequencies.The prototype of such processes is fractional Brownian motion (FBM).It is described by a single parameter, called the Hurst exponent, which enters the power law.In contradistinction to Weierstrass functions, FBM turns out to be surprisingly accurate in the representation of real sea or soil surfaces.As a matter of fact, many surfaces obtained by either direct measurement or theoretical considerations exhibit a 1/f spectrum over a wide range of scales, with some non-trivial characteristic exponent (see e.g.[6] for sea spectra).Even though the FBM model remains an idealization, it is relevant in so far as the interrogating wave probes the scales at which a fractal regime is observed.The inverse scattering problem in that case is therefore particularly simple in that it reduces to finding a unique parameter.Berry and Blackwell in their pioneering work [4] were the first to establish a (simple) relation between the scattering data from such fractal surfaces and the Hurst exponent.Using the KA and working under normal incidence, they derived a non-trivial power law, governed by the Hurst exponent, for the echo power of a quasi-monochromatic pulse.Recently, several works have been devoted to the study of electromagnetic scattering from FBM or 1/f surfaces.Semianalytical formulae in form of a series expansion have been obtained for the backscattering coefficient using the KA [11] and the integral equation method [17], making possible its numerical estimation and comparison with experimental data.A dependence upon the value of the Hurst exponent at various incidences and wavenumbers has been evidenced but remains qualitative only, the formulae being too involved for a general retrieval method.More precise information can be obtained in the framework of the SPM [9,10]; in that case the backscattering coefficient follows a power law governed by the Hurst exponent and a numerical estimation of the latter is straightforward.A priori confidence intervals on the estimated value are, however, not provided.
In this paper, we set up a quantitative method for the estimation of the Hurst exponent from a single FBM profile.The method will not be bound to the backscattering direction or to the SPM and the estimation error will be controlled.
For this first step we will restrict ourselves to one-dimensional surfaces (that is with an invariance direction), the principle of the method remaining essentially the same in higher dimension.Such surfaces possess several features that render their study more arduous than 'classical' rough surfaces: they are non-stationary, so their roughness cannot be controlled by a single height parameter; they have no characteristic scale, so it is difficult to single out a resonance domain; they have no meaningful correlation length, due to the persistence of longdistance correlations within the process.Nethertheless, the very structure of the surface makes its identification possible in the scattering data.We will consider three usual approximations: KA, the SPM and the small-slope approximation (SSA).In the framework of either of these approximations and either polarization case (TE or TM), the scattered intensity will be shown to follow a universal power law in frequency governed by the Hurst exponent.On the basis of this result we will build estimators for the latter and discuss their performances.The theoretical properties of these estimators will be checked numerically.Finally, we will test numerically the robustness of the method as one leaves the domain of validity of the different approximations; new phenomena as one enters the 'resonance domain' (so-called) will be unveiled.

Scattering geometry
A z-invariant surface y = f (x) separates the vacuum from a perfectly conducting medium located at y < f (x).A linearly polarized time-harmonic plane wave with wavevector ) is impinging at incidence θ i on the top of the surface (see figure 1).The total field F in the upper medium is written F = F i + F d , where F i is the incident field, and F d is the field scattered by the surface.Here F stands for the electric field E or magnetic field H according to whether the polarization is TE or TM.The boundary conditions on the surface have to be adapted accordingly.The harmonic time dependence e −iωt will always be implicit.Above the highest excursion of the surface, the scattered field admits a Rayleigh expansion: wherein The terms |α| e k in the above integral correspond to propagating waves while |α| e > k correspond to evanescent waves.When |α| e k, B(α i , α e ) is the far-field scattering amplitude in the direction θ e , where α e = k sin θ e .Note that a different normalization is sometimes adopted in the literature for the scattering amplitude, namely β e β i B(α i , α e ), which corresponds to sending plane waves with unit Poynting flux through a unit horizontal surface.
The field above the surface is given by the Kirchhoff-Helmholtz formula: where G is the two-dimensional outgoing Green function: H (1)  0 (kMM ).This equation cannot be solved analytically except for very peculiar cases and for theoretical predictions on the field one has to resort to approximations.To obtain rise of convergence problems in the following integrations, we will assume the height profile to be of finite extent, that is f (x) = 0 for x large enough.As we will see, this is no restriction since only a finite part of the surface will be illuminated.

Kirchhoff approximation
KA is the most popular approximation in scattering from rough surfaces.It relies on the physical optics approximation, which amounts to identifying the surface with its local tangent plane.Although its accuracy is difficult to estimate, KA is known to be reliable as the incident wavenumber becomes much larger than the maximum surface wavenumber, provided the local incidence angle remains low.
In the TE case, the (electric) field satisfies a Dirichlet boundary conditions on the surface.In the KA this yields the condition dE dn = 2 dE i dn for the normal derivative on the surface.Thus, The plane wave expansion of the Hankel function H (1)  0 , H (1) β e e iα e (x−x )+iβ e (y−f (x )) , (y > f (x)), then provides us with an expression of the scattering amplitude by identification with (2.1): Note that for f = 0 we have which coincides, as expected, with the scattering amplitude of a plane mirror.The calculation can be carried a little bit further.A straightforward calculus using integration by parts yields where the term A KA (θ i , θ e ) is an angular factor which does not depend on the surface or on the wavenumber k but only on the geometry of the problem: In the TM case, the (magnetic) field satisfies a Neumann boundary condition on the surface which is easily shown to imply H = 2H i .Calculations similar to the previous case lead to the same expression (2.3), with a different angular factor cos θ e (cos θ i + cos θ e ) .

Small-perturbation method
The SPM consists in expanding the scattering coefficient in a perturbation series with respect to a height parameter, namely the root mean square (RMS) height of the profile.It is exact to the first order in height as the size of the asperities goes to zero (in contrast to the KA, see [18] and discussion therein).Setting f (x) = hs(x), where s has unit RMS, we have at first order in h with angular factors depending on the polarization in the TE case, and in the TM case.Here ŝ represents the Fourier transform of s, that is

Small-slope approximation
The SSA was introduced more recently by Voronovitch [25] to bridge the gap between the KA and the SPM.Its consists in expanding the scattered field in a perturbation series with respect to a slope parameter (essentially the Fourier transform of the profile), assuming the latter is sufficiently small.We refer to [23] for a detailed computation in the TE case for one-dimensional surfaces and to [25,26] for the general case.The method is consistent with the SPM in that the small-slope series reduces to the small-(height-) perturbation series as the quantity kh goes to zero.Moreover, the first two orders coincide with the KA in the high-frequency limit.To first order in slope we have with 2 cos θ i cos θ e + cos θ i in the TE case, and 1 − sin θ e sin θ i cos θ e (cos θ e + cos θ i ) in the TM case.Note that the SSA is consistent with the SPM as the maximum amplitude of the profile goes to zero.Indeed, the linearization of the exponential term in (2.5) brings us back to (2.4).This is not the case for the KA, as can be seen using the same procedure.Note also that the expression of the SSA is similar to the KA, apart from an angular factor that does not depend on the profile.

Gaussian beams
In numerical as well as physical applications, the illuminated area is of finite extent; the incident beam is thus no longer a plane wave but is assumed to have a Gaussian distribution around some mean incidence.Therefore the incident field takes the form where and L characterizes the half-size of the illuminated region.The normalization of the incident beam is chosen so that the total intensity sums to unity: The scattered amplitude r(α i , α e ) is obtained by summing the contributions of all possible incidences, that is In the KA this leads to Here we will make two natural assumptions.First, we assume that the characteristic half-size L of the illuminating beam satisfies the following.

Hypothesis 2.1 (Large beam). kL| sin θ
This means at the same time that the size of the illuminated patch is much larger than the incident wavelength and that the scattering angle is taken away from the specular region.Second, we suppose that the maximum height variation y = max |f (x) − f (x )| of the surface is small compared to the illuminated length.

Hypothesis 2.2 (Moderate height). y L.
With these assumptions we derive easily in the KA.Under the same assumptions we obtain with the SPM.The SSA yields the same formula as (2.6) with angular factor A KA (θ i , θ e ) replaced by A SSA (θ i , θ e ).Note that the integrands appearing in both formulae are windowed by Gaussian functions and thus the effective range of integration is of finite extent, thereby justifying our assumption that the profile f (x) be set to zero at long distances.

Fractional Brownian motion (FBM)
For a given 0 < H < 1 called the Hurst exponent, FBM [16] can be defined as the unique H -self-similar Gaussian process.This means f (at) = a H f (t) for all t > 0 and a > 0, this equality holding for all finite-dimensional distributions.Equivalently, this is the unique zero-mean Gaussian process with covariance where σ is some height parameter (in the physics literature, this parameter is sometimes expressed in terms of the so-called 'topothesy' T , which is the distance over which chords joining points on the surface have RMS slope equal to unity: σ = T 1−H ).The increments of FBM are stationary and follow the normal distribution The case H = 1/2 corresponds to usual Brownian motion.Note that FBM is non-stationary and as such does not have a properly defined power spectrum.Nethertheless, the bi-dimensional Fourier transform of its covariance function can be defined in the sense of distributions.It is written where C 2H is some positive constant (see the appendix).Recalling that the bi-dimensional Fourier transform of a stationary covariance is of the form Ĉ(ξ, η) = δ(ξ + η) (ξ ) for some positive measure called the power spectrum, we then can associate a pseudo-power-spectrum with FBM by identification with (3.3): A less ad hoc definition of the spectrum of FBM can be given by means of the wavelet transform [7].
The sample paths of FBM are continuous but not differentiable; precisely they are Hölderian of order ν (that is |f Hence the normal derivative on the surface does not exist and the electromagnetic boundary problem is not defined.We expect, however, the diffracted field to be insensitive to the details of the profile that are small compared to the incident wavelength.A way to give a meaning to the scattering problem is therefore to build a sequence of smoothed versions of the profile at finer and finer resolution and to consider the limit of the corresponding diffraction diagrams, provided this limit exists.Precisely, let g(x) = 1/ √ 2π exp (−x 2 /2) and g a (x) = a −1 g(x/a).Then we have pointwise g a * f (x) → f (x) as a → 0 while the scattering amplitude r a (α i , α e ) of the smoothed profile g a * f is well defined.We are not aware of any mathematical result stating the convergence of the scattering amplitude and we will assume this actually holds, that is r a (α i , α e ) → r(α i , α e ), a → 0. This hypothesis permits us to carry over the computations as if the surface were an exact mathematical FBM.This has the advantage of considerably simplifying the statistical estimations.Note, however, that a convergence test can be performed numerically.An illustration is given in figure 4, where the diffraction diagrams of smoothed FBMs profiles become indistinguishable as the resolution 1/a increases.
FBM belongs to the class of fractal processes, in the sense that it has a non-integer Hausdorff-Besicovitch dimension D HB .Recall that the latter is defined by N( ) ∼ −D HB as goes to zero, where N( ) is the minimum number of boxes of size that are necessary to cover the graph of the process over a finite interval.The dimension D HB is trivially related to the Hurst exponent by D HB = 2 − H (see [1]).Hence the parameter H can also be seen as a fractal dimension.
Being non-stationary, the height variations of FBM cannot be controlled by a single parameter.Note, however, that the maximum excursion over a distance L is of order σ L H . Hence, the relevant height parameter is now the combination of two parameters: a 'classical' RMS term and a length parameter.For further reference we also note at once that FBM paths over a finite distance satisfy the growth requirement hypothesis 2.2 provided.

Hypothesis 3.1 (Moderate height). σ L H
As mentioned in the introduction, there is no correlation length for FBM.Indeed, it can easily be seen from (3.1) that which does not even decrease at long lags.

Simulation techniques
In view of a numerical validation of the forthcoming results, a reliable simulation procedure is needed.Consider the problem of generating an exact sample of FBM on the interval, say, [0, 1], that is a sequence (X(j ) = f (j/N) j = 0, N − 1), with prescribed covariance matrix They are several simulation techniques to generate such samples.The simplest method in principle is Cholevsky decomposition.This amounts to factorizing the covariance matrix as Cov X = MM t , M t being the Hermitian transpose of the N × N matrix M. Once the corresponding matrix M has been found, one only needs to generate N Gaussian i.i.d variables W = (W 1 , . . ., W N ) and to set X = MW , which by construction has the desired covariance matrix.This method is exact but suffers from dramatic numerical limitations since the storage requirement is in O(N 2 ) (in practice this limits the sample size to N 1200).Many approximate methods have been proposed to overcome this limitation.Without being exhaustive, let us mention the spectral methods, based on the inversion of the (known) Fourier or wavelet coefficients of FBM via an efficient algorithm (fast Fourier transform or fast pyramidal algorithm); bisection methods (random midpoint displacement and conditionalized random midpoint displacement), where points are constructed iteratively by interpolating existing points according to a 'top-down' algorithm; aggregation methods such as superposition of ON/OFF sources.There exists, however, a simple, exact and at the same time efficient method to generate arbitrary stationary Gaussian processes with prescribed covariance [27], a method that can be used to generate the increments Y j = X j +1 − X j of FBM (which are stationary).It consists in extending the original stationary process Y j to a cyclostationary process Ỹj (of size M 2N ), which amounts to embedding the N × N covariance matrix Cov Y in a larger matrix Cov Ỹ (of size M × M) that is chosen to be symmetric and circulant.Now the diagonalization (and consequently the Cholevsky decomposition) of such matrices can be performed in an efficient manner by means of the fast Fourier transform.This way one can generate a vector Ỹ with covariance matrix Cov Ỹ , whose restriction Y = (Y 1 , . . ., Y N ) has covariance matrix Cov Y .This method requires little storage space and can comfortably (and quickly!) deal with samples up to N = 2 20 .Figure 2 shows an example of simulated FBM with H = 0.6 and

Domain of validity of the approximations
For 'classical' random rough surfaces the domains of validity of the aforementioned approximations are well characterized in terms of relative magnitude of the wavelength, the RMS height, the correlation length, the mean square slope, the radius of curvature or the profile etc.For fractal and non-stationary surfaces such as FBM there are no such criteria and we wish to provide some heuristic scaling arguments to establish the different regimes of validity.
Using the KA for non-differentiable surfaces is a priori impossible since no tangent plane can be defined.However, as mentioned earlier, only the details that are comparable to or greater than the incident wavelength will contribute to the scattering, so the effective scattering profile is a smoothed version of the mathematical FBM at the resolution of the wavelength.Note also that the integral in (2.3) is well defined as soon as the profile is continuous, so Kirchhoff's formula can always be used a posteriori.Intuitively, the KA is valid if the effective scattering surface does not deviate much from a mean plane over some wavelengths.This can be rephrased mathematically as where n is some small number, a = (f (x + λ/2) − f (x − λ/2))/λ is the local slope around x and the brackets represent the RMS mean.Now Hence our criterion for the KA is written as follows.
The SSA is correct when the height variations over one wavelength are small compared to the wavelength.A criterion of validity is therefore Since f (x + λ) − f (x) = σ λ H , this amounts to the following assumption.

Hypothesis 3.3 (SSA domain). σ λ H
Note that the SSA regime together with the large-beam hypothesis kL 1 implies hypothesis 3.1 since σ L H −1 σ k 1−H 1.The SPM holds when the RMS height is small compared to the wavelength.For an illuminated patch of typical size L, the maximum RMS height with respect to its mean plane is σ L H . Therefore the following criterion ensures the validity of the SPM.

Hypothesis 3.4 (SPM domain). σ L
Note that the SPM regime hypothesis 3.4 together with the large-beam hypothesis 2.1 also imply hypothesis 3.1 since σ L H −1 σ kL H 1.

Statistical properties of the diffracted field
We will now study the statistical properties of the diffracted field in the framework of the KA, the SSA and the SPM.Since the KA and SSA have the same expression, apart from an angular factor, they will be treated simultaneously.From now on, the surface under consideration will be a perfect FBM surface with given Hurst exponent, illuminated at incidence θ i by a Gaussian beam of characteristic half-size L and central wavenumber k.where Proof.From (2.6) and (3.2) it follows easily, away from the specular direction, that where A(θ i , θ e ) represents A KA (θ i , θ e ) or A SSA (θ i , θ e ) as appropriate.This gives, after a simple change of variables, where we have set We do not know of analytical expressions for these integrals, so we rather look at their asymptotic behaviour.In view of the large-beam hypothesis, we set L|α e − α i | = a (or ξ = aτ ) and look at the limit a → ∞.Note that ξ ∼ k 1−1/H σ −1/H 1 by the small-slope hypothesis.The conclusion follows from lemma 10.1 given in the appendix.
Note that the mean diffuse field is zero, as it is the case for stationary surfaces.and Proof.Since FBM is a zero-mean process, we obviously have E f P L (α) = 0 for all α.This entails (4.3).For the second moment we have, away from the specular direction, where Ĉ(ξ, η) is the bi-dimensional Fourier transform of the covariance function C(x 1 , x 2 ) of FBM.Outside the specular direction, this leads to which, for L|α − α i | 1, yields (4.4).

Building an estimator
According to the KA, SPM as well as SSA, the scattered intensity has a power-law dependence on the frequency k, with ν = 1 − 2H .The most natural way of computing the Hurst exponent H is therefore to work at fixed angles of incidence θ i and emergence θ e and to perform a log-log regression on I 1 (k) := |r(α i , α e )| 2 .Let k 1 , . . ., k N be the available wavenumbers and set Then the first estimator is written w j y j and w j = Nx j − S x NS xx − S 2 x .
(5.3) (Note that the weights satisfy w j = 0 and x j w j = 1.)We shall now study the properties of this estimator in the framework of SPM.We have not been able to do the same for the KA and only numerical simulations will be presented for the latter.We will suppose the following.
In that case, the estimator Ĥ1 is essentially a spectral estimator (that is, Fourier based).We do not know of rigorous results on the estimation of the Hurst parameter of FBM from the continuous Fourier transform of a sample path.Precise informations have been given in [15] for the spectral estimator based on the discrete Fourier transform of the increments of FBM, showing that the latter is asymptotically unbiased with a variance of order log 2 (N)/N, where N is the number of sampling points on the profile.As we will see, some simplifying but reasonable assumptions make it possible to derive comparable performances for the estimator Ĥ1 .First we need to know the statistical nature of the diffracted intensity I 1 (k).
Proposition 5.2.Under hypothesis 5.1 and away from the specular direction, I 1 (k) is an exponential variable.
Proof.With the SPM, it follows from equation (2.7) that |r(α i , α e )| 2 is the Fourier transform of a Gaussian variable, hence a complex Gaussian variable.Therefore |r(α i , α e )| 2 is a chisquared variate with two degrees of freedom, that is simply an exponential variable (i.e. a variable X such that P (X t) = e −μt , t 0).
Next we need to estimate the decorrelation of intensities at different intensities.Define as the degree of correlation between the scattered intensities.Then the following result holds.Proof.For notational convenience set X L (α) = e −iαx f (x)P L (x) dx and α = α e − α i , α = α e − α i .In view of (2.7) it suffices to estimate the correlation term

Proposition 5.3. Under hypotheses 5.1 and 2.1, we have
×P L (x, y)P L (x , y ) dx dy dx dy where P L (x, y) = P L (x)P L (y).Now recall the well known formula for the higher-order correlation functions of a Gaussian process: where the sum is taken over all subdivisions of the set (x 1 , . . ., x n ) into pairs (u l , u m ) and the product is to be taken over all pairs of the corresponding subdivision.For n = 4 this yields Now, for all α 1 , α 2 , where by (3.3), In the regime α j L 1, j = 1, 2, which is ensured by hypothesis 2.1, we have ). Hence T 1 = T 2 = 0 and it remains Noting that the conclusion follows.
To make the evaluation of statistical performances possible, we will adopt the following crucial assumption.

Hypothesis 5.4 (Perfect decorrelation). The intensities corresponding to two wavenumbers
k, k such that |k − k| 1/L are perfectly decorrelated, i.e.Cor(α i , α e ; α i , α e ) = 0.Under this assumption it is possible to evaluate the performances of the estimator Ĥ1 .Theorem 5.5.For a fixed geometry the different available wavenumbers with |k j − k j | > 1/L, j = j .Suppose the wavenumbers are either in arithmetical (i.e.k j +1 − k j = cst) or geometrical progression (i.e.log k j +1 − log k j = cst).
Then under hypotheses 2.1, 3.1, 5.1 and 5.4, the estimator Ĥ1 is unbiased with asymptotic variance of order 1/N : for some constant C depending only on k − , k + and the type of progression.This constant is smaller if the wavenumbers are in geometrical progression.
Proof.Under hypotheses 2.1, 3.1 and 5.1, the power-law regime (5.1) is ensured by theorem 4.2.This, together with proposition 5.2, makes it possible to rewrite I 1 (k) = σ 2 /2I 0 (θ i , θ e )k ν Z, where Z is an exponential variable with mean 2 (equivalently, a standard chi-square variate with two degrees of freedom).It follows from the results recalled in appendix that Ey j = ν log k j + constant and Var(y j ) = ζ(2, 1).Hence, It follows at once that Ĥ1 is an unbiased estimator of H . Let k − and k + be the lowest and highest available wavenumbers.The point is now how to choose the intermediate wavenumbers so as to minimize the variance.There are essentially two choices: a regular linear sampling (that is ).An elementary but tedious calculation shows that lim for large N.Although we have not been able to prove it analytically, a numerical estimation clearly shows that C 2 C 1 for all (k − , k + ), so that a further variance reduction can be achieved by choosing a logarithmic sampling.
The variance of Ĥ1 is thus no better than O(1/N) in the most favourable case where the SPM is exact and the scattered intensities at different wavenumbers are perfectly decorrelated.For a fixed maximum frequency k + , note that the maximum number N of uncorrelated wavenumbers is of order O(L) as the size of the illuminated patch is increased and thus Var( Ĥ1 ) is O(1/L).

An improved estimator
The main drawback of the estimator Ĥ1 is that it requires the sampling of a large set of wavenumbers.Now, from a practical point of view, it is sometimes difficult to obtain a large set of wavenumbers, a fortiori a large set of uncorrelated wavenumbers.For a laser beam experiment, for instance, only a sparse set of wavenumbers is available.Numerically, working at varying frequency is computationally very demanding.Sampling over different emergence angles, however, is experimentally viable and at the same time of low computational cost in numerical simulations.Therefore, it seems preferable to design an estimator on a varying-angle rather than varying-frequency scheme.Suppose we have only a few wavenumbers available, say k 1 < • • • < k N with a regular sampling on a logarithmic scale.For a fixed frequency k and incidence angle θ i , we can achieve several values of the momentum transfer α − α i by varying the scattering angle θ while keeping, however, the same order of magnitude for the latter.Let θ vary within a cone of angular aperture (θ − e , θ + e ) around some mean non-specular direction θ e .Now define the angular mean: where α j = k sin θ j , θ j = θ − e + j θ and θ = (θ + e − θ − e )/N θ .The angular sampling rate θ is taken to be the smallest that allows us to keep the successive values of |r(α i , α j )| 2 uncorrelated at all wavenumbers.This is the case as soon as k − (sin (θ e + θ) − sin θ e ) 1/L, that is roughly as soon as the criterion is satisfied.Note that the number of available emergence angles is then of order O(L) as the size of the illuminated patch is increased.In the KA, SSA or SPM we have where the expression of the angular factor I (θ i , θ j ) is given by (4.4) or (4.2).Here we will need an additional simplification, namely that the angular factor does not vary appreciably within the prescribed angular cone.2N θ I 0 (θ i , θ e )k ν Z, where Z is a chi-square variate with 2N θ degrees of freedom.Again, we need only perform a log-log regression on I 2 (k).The corresponding estimator Ĥ2 is then given by equations (5.2) and (5.3),where I 1 (k) is to be replaced by I 2 (k).The point is that I 2 (k) has a much smaller variance than I 1 (k) since it is the result of summing identical uncorrelated variables.Theorem 6.2.Under hypotheses 2.1, 3.1, 5.1, 5.4 and 6.1, the estimator Ĥ2 is unbiased and, for the same choice of wavenumbers (k j ), its variance is reduced by a factor N θ with respect to Ĥ1 .
Proof.From appendix we now have Ey j = ν log k j + cst and Var(y j ) = ζ(2, N θ /2).Hence, Eν = w j Ey j = ν w j x j = ν, Again Ĥ2 is an unbiased estimator of H .Moreover, comparison with the previously calculated variance shows that Var( Ĥ2 ) This completes the proof.Confidence intervals.Once the variance Var( Ĥ ) = s 2 of the estimators has been estimated, the 95% confidence interval for the estimated value of H is [H − 1.96s, H + 1.96s].

Numerical performances of the estimator
As an illustration we have computed the scattering intensity on a simulated FBM with H = 0.7 according to KA (2.6) and SPM (2.7).We have set the width of the incident beam to L = 256 and the height parameter to σ = 0.05 and let the incident wavelength λ range from 0.1 to 16.We have chosen the backscattering direction with θ i = 25.With these values, hypotheses 2.1, 3.1 and 3.3 are clearly verified with L(α e −α i ) ∈ [85, 13 632], σ L H −1 = 0.01, σ λ H −1 ∈ [0.02, 0.1].Note that the unit here is arbitrary (all that matters is the ratio between the different quantities).The mean intensity I 1 (k) has been computed for 20 logarithmically sampled wavenumbers.The angular-averaged intensity I 2 (k) has been computed for 20 wavenumbers via an angular mean of the scattered intensity over ten emergence angles in a 5 • cone about the backscattering direction.With these values kL cos θ e θ ranges from 0.8 to 127, and therefore the decorrelation condition (6.1) is also satisfied.To check out the calculations on the theoretical mean and variance of the estimators given in (5.8) and (6.2), we have computed the corresponding empirical quantities on 268 samples.The result is summarized in table 1.
The discrepancy between the theoretical and empirical values for the SPM is less than 3%, thereby making very realistic the simplifying assumptions that have been adopted.

Exact results in the resonance domain
A challenging question is whether the power-law behaviour that is observed for the scattered intensity under the different single-scattering approximations (KA, SSA and SPM) continues to hold in the resonance domain.The term 'resonance domain' is somewhat inappropriate since no typical scale can be singled out on a fractal surface and it will be simply understood as a domain where the multi-scattering phenomena become important.To derive general results in this regime seems a difficult task since no analytical formula can be relied on to guide the study.We wish, however, to present some numerical results supporting the conjecture that a power-law behaviour persists far beyond the domain of validity of the aforementioned approximations.
A rigorous numerical solution of the scattering problem can be achieved by solving the contour integral equation (2.2) satisfied by the field.For a description of the numerical algorithm and its accuracy we refer to [20].The diffraction diagrams of FBM surfaces in the TE case have been computed for various values of H , σ and λ, and for a fixed incidence θ i = 20 • ; in each case averages over 20 realizations of the surface have been performed in order to obtain the mean intensity E|r(α i , α e )| 2 .The characteristic half-size of the illuminated patch has been set to L = 256 and the wavelength ranges from 2 to 16, so that the large-beam condition (2.1) is clearly satisfied.Two typical values of H have been chosen, namely H = 0.5 and 0.7.The values of the height parameter σ have been chosen so as to leave progressively the domain of validity of the different approximations (according to our heuristic criteria (3.2)-(3.4)).A numerical comparison has been made with the SSA in order to identify the different regimes (single or multi-scattering).
The summary provided by table 2 gives a complete recap of the values or ranges of values of the various parameters.The corresponding diffraction diagrams have been labelled 'small', 'medium' or 'large' according to the increasing values of the height parameter σ .Figures 5  and 6 show the mean backscattered intensity EI 2 (k) = E|r(α i , −α i )| 2 in a log-log diagram for the different values of H and σ , at incidence θ i = 20 • .A further angular average has been performed within a 5 • cone about the backscattering direction to smooth out the diagrams.The curves with suffix '.ex' correspond to the exact field while the '.ss' curves correspond to the SSA.For each plot we have performed a least-squares linear regression to estimate the power exponent ν.The different values of μ = (1 − ν)/2 are reported in table 3. From these results we may draw several conclusions:  • A non-trivial power-law is observed for the mean backscattered intensity far beyond the SSA regime, that is EI 2 (k) ∼ k −ν .• The corresponding exponent μ increases with the height parameter σ , starting from the value predicted by the SSA: μ = H . • The exact solution starts departing from the SSA for σ L H −1 0.5.For σ L H −1 0.5, both the exact solution and the SSA give a good estimate of the Hurst exponent, namely μ = H .

Conclusion
In this paper we have studied the scattering from a perfectly conducting FBM-shaped profile.
We have derived heuristically the different domains of validity of the usual approximations.
The scattered intensity has been computed under each of these approximations, resulting in a remarkable generic feature: a power law with respect to the wavenumber with the same universal exponent 2H − 1, albeit with different multiplicative constant.This has suggested a method to estimate the Hurst exponent from the scattering diagrams; the performances of the corresponding estimator have been investigated both theoretically and numerically.Finally, we have computed the exact diffraction diagrams in the resonance domain and evidenced the persistence of a power law far beyond the domain of validity of the approximations.The numerical computation of the scattering diagram on a surface with a large number of wavelengths imposes severe time limitations and we restricted the study to two characteristic values of H and three increasing values of the height parameter σ .A reasonable conjecture is that the observed phenomena are generic with respect to the value of H .However, the behaviour of the power-law exponent μ as one continues to increase σ remains totally unknown and it is an open question whether it eventually reaches some non-trivial limit value depending on H .This problem is left for further research.

Figure 1 .
Figure 1.Geometry of the problem.

Figure 5 .
Figure 5. Log-log diagram of the mean backscattered intensity k −2 EI 2 (k) for an incidence θ i = 20 • (the factor k −2 is for better visibility).Theoretical value of the Hurst exponent H = 0.5.

Table 1 .
Performances of the estimators Ĥ1 and Ĥ2 over 268 realizations.μ and s are the theoretical mean and standard deviation; μ and s are the corresponding empirical quantities.

Table 2 .
Relevant scale parameters of the FBM surface for H = 0.5 and 0.7 and for different values of σ .The illuminated patch has characteristic half-size L = 256 and the incident wavelength ranges from 2 to 16.

Table 3 .
Values of the slope μ estimated from a log-log regression on the mean backscattered intensity (figures 5 and 6).
= (z)/ (z) is the psi function and ζ(z, n) a generalized Rieman zeta (or Hurwitz zeta) function.Note that for n = 2, Z 2 is simply an exponential variable with mean 2.