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

Charles-Antoine Guérin, Marc Saillard

► To cite this version:

Charles-Antoine Guérin, Marc Saillard. Electromagnetic scattering on fractional Brownian surfaces and estimation of the Hurst exponent. Inverse Problems, 2001, 17 (3), pp.365-386. 10.1088/0266-5611/17/3/301 . hal-00083199

HAL Id: hal-00083199

https://hal.science/hal-00083199v1

Submitted on 10 Jul 2016

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Electromagnetic scattering on fractional Brownian surfaces and estimation of the Hurst exponent

Charles-Antoine Guérin and Marc Saillard

Institut Fresenel, UMR CNRS 6133, Faculté des Sciences de Saint-Jérôme, case 162, F-13397
Marseille Cedex 20, France

E-mail: caguerin@loe.u-3mrs.fr

## Abstract

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 small-perturbation 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.

## 1. 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 band-limited 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. Semi-analytical 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 long-distance correlations within the process. Nethertheless, the very structure of the surface makes
<figure><img src="image_2.png" /><figcaption>Figure 1. Geometry of the problem.</figcaption></figure>

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.

## 2. Electromagnetic scattering on 1D surfaces

### 2.1. 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 $k_i = (\alpha_i, -\beta_i) = (k \sin \theta_i, -k \cos \theta_i)$ is impinging at incidence $\theta_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,

$$F_i(x, y) = e^{i\alpha_i x - i\beta_i y},$$

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\omega t}$ will always be implicit. Above the highest excursion of the surface, the scattered field admits a Rayleigh expansion:

$$F_d(x, y) = \int_{-\infty}^{+\infty} d\alpha_e B(\alpha_i, \alpha_e) e^{i\alpha_e x + i\beta_e y}, \quad (2.1)$$

wherein $\beta_e = \sqrt{k^2 - \alpha_e^2}$ if $\alpha_e^2 \le k^2$ and $\beta_e = i\sqrt{\alpha_e^2 - k^2}$ if $\alpha_e^2 \ge k^2$. The terms $|\alpha|_e \le k$ in the above integral correspond to propagating waves while $|\alpha|_e > k$ correspond to evanescent waves. When $|\alpha|_e \le k$, $B(\alpha_i, \alpha_e)$ is the far-field scattering amplitude in the direction $\theta_e$, where
$\alpha_e = k \sin \theta_e$. Note that a different normalization is sometimes adopted in the literature for the scattering amplitude, namely $\sqrt{\frac{\beta_e}{\beta_i}} B(\alpha_i, \alpha_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:

$$F(M) = F_i(M) + \int_S G(M, M') \frac{dF}{dn'} dM' - \int_S F(M') \frac{dG}{dn'} (M, M') dM' \quad (2.2)$$

where $G$ is the two-dimensional outgoing Green function:

$$G(M, M') = -\frac{i}{4} H_0^{(1)}(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.

## 2.2. 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 $\frac{dE}{dn} = 2\frac{dE_i}{dn}$ for the normal derivative on the surface. Thus,

$$E(x, y) = E_i(x, y) - \frac{1}{2} \int H_0^{(1)}(kMM')(\alpha_i f'(x') + \beta_i) e^{i\alpha_i x' - i\beta_i f(x')} dx',$$

with $MM' = \sqrt{(x - x')^2 + (y - f(x'))^2}$. The plane wave expansion of the Hankel function $H_0^{(1)}$,

$$H_0^{(1)}(kMM') = \frac{1}{\pi} \int_{-\infty}^{+\infty} \frac{d\alpha_e}{\beta_e} e^{i\alpha_e(x-x') + i\beta_e(y-f(x'))}, \quad (y > f(x)),$$

then provides us with an expression of the scattering amplitude by identification with (2.1):

$$B(\alpha_i, \alpha_e) = -\frac{1}{2\pi\beta_e} \int (\alpha_i f'(x') + \beta_i) e^{i(\alpha_i - \alpha_e)x' - i(\beta_i + \beta_e)f(x')} dx.$$

Note that for $f = 0$ we have

$$B(\alpha_i, \alpha_e) = -\frac{\beta_i}{2\pi\beta_e} \int e^{i(\alpha_i - \alpha_e)x} dx = -\delta(\alpha_i - \alpha_e),$$

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

$$B(\alpha_i, \alpha_e) = A_{KA}(\theta_i, \theta_e) \frac{1}{2\pi} \int e^{i(\alpha_i - \alpha_e)x - i(\beta_i + \beta_e)f(x)} dx, \quad (2.3)$$

where the term $A_{KA}(\theta_i, \theta_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:

$$A_{KA}(\theta_i, \theta_e) = -\frac{1 + \cos(\theta_i + \theta_e)}{\cos\theta_e(\cos\theta_i + \cos\theta_e)}.$$
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

$$A_{KA}(\theta_i, \theta_e) = \frac{\cos 2\theta_i + \cos (\theta_i - \theta_e)}{\cos \theta_e (\cos \theta_i + \cos \theta_e)}.$$

## 2.3. 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

$$B(\alpha_i, \alpha_e) = \varepsilon\delta(\alpha_i - \alpha_e) + khA_{SPM}(\theta_i, \theta_e)\frac{1}{2\pi}\hat{s}(\alpha_e - \alpha) \quad (2.4)$$

with angular factors depending on the polarization

$$A_{SPM}(\theta_i, \theta_e) = 2i \cos \theta_i, \quad \varepsilon = -1$$

in the TE case, and

$$A_{SPM}(\theta_i, \theta_e) = -2i \frac{1 - \sin \theta_i \sin \theta_e}{\cos \theta_e}, \quad \varepsilon = +1$$

in the TM case. Here $\hat{s}$ represents the Fourier transform of s, that is

$$\hat{s}(\alpha) = \int e^{-i\alpha x} s(x) dx.$$

## 2.4. 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

$$B(\alpha_i, \alpha_e) = A_{SSA}(\theta_i, \theta_e) \frac{1}{2\pi} \int e^{i(\alpha_i - \alpha_e)x - i(\beta_i + \beta_e)f(x)} dx, \quad (2.5)$$

with

$$A_{SSA}(\theta_i, \theta_e) = - \frac{2 \cos \theta_i}{\cos \theta_e + \cos \theta_i}$$

in the TE case, and

$$A_{SSA}(\theta_i, \theta_e) = 2 \frac{1 - \sin \theta_e \sin \theta_i}{\cos \theta_e (\cos \theta_e + \cos \theta_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.
## 2.5. 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

$$F_i(x, y) = \int \hat{P}_L(\alpha - \alpha_i) e^{i\alpha x - i\beta y} d\alpha,$$

where

$$\hat{P}_L(\alpha) = (2\pi)^{1/4} L^{1/2} e^{-\alpha^2 L^2/4}, \quad P(x) = (2/\pi)^{1/4} L^{-1/2} e^{-x^2/L^2},$$

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:

$$\int |\hat{P}_L(\alpha)|^2 d\alpha = 1.$$

The scattered amplitude $r(\alpha_i, \alpha_e)$ is obtained by summing the contributions of all possible incidences, that is

$$r(\alpha_i, \alpha_e) = \int \hat{P}_L(\alpha - \alpha_i) B(\alpha, \alpha_e) d\alpha.$$

In the KA this leads to

$$r(\alpha_i, \alpha_e) = A_{KA}(\theta_i, \theta_e) \frac{1}{2\pi} \int d\alpha \int \hat{P}_L(\alpha - \alpha_i) e^{i(\alpha - \alpha_e)x - i(\beta + \beta_e)f(x)} dx.$$

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\theta_e - \sin\theta_i| \gg 1$.

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 $\Delta y = \max|f(x) - f(x')|$ of the surface is small compared to the illuminated length.

**Hypothesis 2.2 (Moderate height).** $\Delta y \ll L$.

With these assumptions we derive easily

$$r(\alpha_i, \alpha_e) = A_{KA}(\theta_i, \theta_e) \frac{1}{2\pi} \int e^{i(\alpha_i - \alpha_e)x - i(\beta_i + \beta_e)f(x)} P_L(x) dx \quad (2.6)$$

in the KA. Under the same assumptions we obtain

$$r(\alpha_i, \alpha_e) = \frac{k}{2\pi} A_{SPM}(\theta_i, \theta_e) \hat{P}_L f(\alpha_e - \alpha_i) \quad (2.7)$$

with the SPM. The SSA yields the same formula as (2.6) with angular factor $A_{KA}(\theta_i, \theta_e)$ replaced by $A_{SSA}(\theta_i, \theta_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.
### 3. Surface modelling

#### 3.1. 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

$$C(x_1, x_2) = \frac{\sigma^2}{2} (|x_1|^{2H} + |x_2|^{2H} - |x_1 - x_2|^{2H}), \qquad (3.1)$$

where $\sigma$ 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: $\sigma = T^{1-H}$). The increments of FBM are stationary and follow the normal distribution $f(x) - f(x') \sim N(0, \sigma^2|x - x'|^{2H})$ with characteristic function

$$\mathbb{E}e^{-i\xi(f(x)-f(x'))} = e^{-\xi^2 \frac{\sigma^2}{2}|x-x'|^{2H}}. \qquad (3.2)$$

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. Netheress, the bi-dimensional Fourier transform of its covariance function can be defined in the sense of distributions. It is written

$$\hat{C}(\xi, \eta) = -C_{2H} \frac{\sigma^2}{2} (|\xi|^{-1-2H} \delta(\eta) + |\eta|^{-1-2H} \delta(\xi) - \delta(\xi + \eta)|\xi|^{-1-2H}), \qquad (3.3)$$

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 $\hat{C}(\xi, \eta) = \delta(\xi + \eta)\Gamma(\xi)$ for some positive measure $\Gamma$ called the power spectrum, we then can associate a pseudo-power-spectrum with FBM by identification with (3.3):

$$\Gamma(\xi) = C_{2H} |\xi|^{-1-2H}, \quad \xi \neq 0. \qquad (3.4)$$

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 $\nu$ (that is $|f(x + h) - f(x)| \le O(|h|^\nu)$ for all $\nu < 2H$). 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/\sqrt{2\pi} \exp(-x^2/2)$ and $g_a(x) = a^{-1}g(x/a)$. Then we have pointwise $g_a * f(x) \to f(x)$ as $a \to 0$ while the scattering amplitude $r_a(\alpha_i, \alpha_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(\alpha_i, \alpha_e) \to r(\alpha_i, \alpha_e)$, $a \to 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(\epsilon) \sim \epsilon^{-D_{HB}}$ as
$\epsilon$ goes to zero, where $N(\epsilon)$ is the minimum number of boxes of size $\epsilon$ 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 $\sigma 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).** $\sigma L^{H-1} \ll 1$.

As mentioned in the introduction, there is no correlation length for FBM. Indeed, it can easily be seen from (3.1) that

$$ \mathbb{E}f(x)f(x+l) \sim l^{\max(2H-1,0)}, \quad l \to \infty, $$

which does not even decrease at long lags.

## 3.2. 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

$$ \mathrm{Cov}_X(j, j') = \mathbb{E}X(j)X(j') = \frac{\sigma^2}{2} \left( \left| \frac{j}{N} \right|^{2H} + \left| \frac{j'}{N} \right|^{2H} - \left| \frac{j - j'}{n} \right|^{2H} \right). $$

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 $\mathrm{Cov}_X = MM^t$, $M^t$ being the Hermitian transpose of the $N \times 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, \dots, 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 \simeq 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 $\tilde{Y}_j$ (of size $M \ge 2N$), which amounts to embedding the $N \times N$ covariance matrix $\mathrm{Cov}_Y$ in a larger matrix $\mathrm{Cov}_{\tilde{Y}}$ (of size $M \times 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 $\tilde{Y}$ with covariance matrix $\mathrm{Cov}_{\tilde{Y}}$, whose restriction $Y = (Y_1, \dots, Y_N)$ has covariance matrix $\mathrm{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
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Example of an FBM path with $H = 0.6$ and $\sigma = 0.027$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Details around the central region.</figcaption></figure>

$\sigma = 0.027$ over a length 2048 with 65536 sampling points. Figure 3 shows the details around the central region. The profile has been recentred around zero (a 'true' FBM should starts from zero: $f(0) = 0$).
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Diffraction pattern on the profile shown in figure 2 for $\lambda = 16$, $L = 256$, and for different degrees of smoothing: $a = 1$ (solid curve), $a = 2$ (long-dashed curve), $a = 10$ (medium-dashed curve), $a = 100$ (short-dashed curve).</figcaption></figure>

### 3.3. 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

$$ \langle f(x') - (ax' + b) \rangle_{[x-n\lambda<x'<x+n\lambda]} \ll 1 $$

where $n$ is some small number, $a = (f(x + \lambda/2) - f(x - \lambda/2))/\lambda$ is the local slope around $x$ and the brackets $\langle \rangle$ represent the RMS mean. Now

$$ \begin{aligned} \mathbb{E}[f(x') - (ax' + b)]^2 &\le 2\mathbb{E}[f(x') - f(x)]^2 + 2\mathbb{E}[a(x - x')]^2 \\ &\le 2\sigma^2(n\lambda)^{2H} + 2\mathbb{E}(a^2)(n\lambda)^2 \\ &= 2n^2\sigma^2(\lambda^{2H} + \lambda^{2H}) = 4n^2\sigma^2\lambda^{2H}. \end{aligned} $$

Hence our criterion for the KA is written as follows.

**Hypothesis 3.2 (KA domain).** $\sigma\lambda^H \ll 1$.
The SSA is correct when the height variations over one wavelength are small compared to the wavelength. A criterion of validity is therefore

$$ \langle f(x + \lambda) - f(x) \rangle \ll \lambda. $$

Since $\langle f(x + \lambda) - f(x) \rangle = \sigma \lambda^H$, this amounts to the following assumption.

**Hypothesis 3.3 (SSA domain).** $\sigma \lambda^{H-1} \ll 1$.

Note that the SSA regime together with the large-beam hypothesis $kL \gg 1$ implies hypothesis 3.1 since $\sigma L^{H-1} \ll \sigma k^{1-H} \ll 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 $\sigma L^H$. Therefore the following criterion ensures the validity of the SPM.

**Hypothesis 3.4 (SPM domain).** $\sigma L^H \lambda^{-1} \ll 1$.

Note that the SPM regime hypothesis 3.4 together with the large-beam hypothesis 2.1 also imply hypothesis 3.1 since $\sigma L^{H-1} \ll \sigma kL^H \ll 1$.

## 4. 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 $\theta_i$ by a Gaussian beam of characteristic half-size $L$ and central wavenumber $k$.

**Theorem 4.1.** Under hypothesis 2.1 and 3.3, we have in the KA and SSA

$$ \mathbb{E} r(\alpha_i, \alpha_e) = 0 \qquad (4.1) $$

and

$$ \mathbb{E} |r(\alpha_i, \alpha_e)|^2 = \sigma^2 I_0(\theta_i, \theta_e) k^{1-2H}, \qquad (4.2) $$

where

$$ I_0(\theta_i, \theta_e) = \begin{cases} \frac{C_{2H}}{8\pi^2} \frac{(1 + \cos(\theta_i + \theta_e))^2}{\cos^2 \theta_e} |\sin \theta_i - \sin \theta_e|^{-2H-1} & \text{KA, TE case} \\ \frac{C_{2H}}{8\pi^2} \frac{(\cos 2\theta_i + \cos(\theta_i - \theta_e))^2}{\cos^2 \theta_e} |\sin \theta_i - \sin \theta_e|^{-2H-1} & \text{KA, TM case} \\ \frac{C_{2H}}{2\pi^2} \cos^2 \theta_i |\sin \theta_i - \sin \theta_e|^{-2H-1} & \text{SSA, TE case} \\ \frac{C_{2H}}{2\pi^2} \frac{(1 - \sin \theta_i \sin \theta_e)^2}{\cos^2 \theta_e} |\sin \theta_i - \sin \theta_e|^{-2H-1} & \text{SSA, TM case.} \end{cases} $$

**Proof.** From (2.6) and (3.2) it follows easily, away from the specular direction, that

$$ \mathbb{E} r(\alpha_i, \alpha_e) = A(\theta_i, \theta_e) \frac{1}{2\pi} \int e^{i(\alpha_i - \alpha_e)x} e^{-(\beta_i + \beta_e)^2 \frac{\alpha_e^2}{2} |x|^{2H}} P_L(x) \, dx, $$

$$ \mathbb{E} |r(\alpha_i, \alpha_e)|^2 = |A(\theta_i, \theta_e)|^2 \frac{1}{4\pi^2} \int e^{i(\alpha_i - \alpha_e)(x_1 - x_2)} e^{-(\beta_i + \beta_e)^2 \frac{\alpha_e^2}{2} |x_1 - x_2|^{2H}} \\ \times P_L(x_1) P_L(x_2) \, dx_1 \, dx_2, $$
where $A(\theta_i, \theta_e)$ represents $A_{KA}(\theta_i, \theta_e)$ or $A_{SSA}(\theta_i, \theta_e)$ as appropriate. This gives, after a simple change of variables,

$$ \mathbb{E}r(\alpha_i, \alpha_e) = 2^{-3/4}\pi^{-5/4}L^{1/2}\tau A(\theta_i, \theta_e) \int e^{i\xi x} e^{-\frac{|x|^2 H}{2}} e^{-\tau^2 x^2} dx, $$

$$ \mathbb{E}|r(\alpha_i, \alpha_e)|^2 = \frac{1}{2\sqrt{2}}\pi^{-5/2}L\tau^2|A(\theta_i, \theta_e)|^2 \int e^{i\xi(x_1-x_2)} e^{-\frac{|x_1-x_2|^2 H}{2}} e^{-\tau^2 x_1^2} e^{-\tau^2 x_2^2} dx_1 dx_2, $$

where we have set

$$ \xi = (\alpha_i - \alpha_e)(\beta_i + \beta_e)^{-1/H} \sigma^{-1/H}, \quad \tau = \xi L^{-1} (\alpha_i - \alpha_e)^{-1}. $$

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|\alpha_e - \alpha_i| = a$ (or $\xi = a\tau$) and look at the limit $a \to \infty$. Note that $\xi \sim k^{1-1/H}\sigma^{-1/H} \gg 1$ by the small-slope hypothesis. The conclusion follows from lemma 10.1 given in the appendix. $\square$

Note that the mean diffuse field is zero, as it is the case for stationary surfaces.

**Theorem 4.2.** Under hypotheses 2.1 and 3.1 we have with the SPM

$$ \mathbb{E}r(\alpha_i, \alpha_e) = 0 \qquad (4.3) $$

and

$$ \mathbb{E}|r(\alpha_i, \alpha_e)|^2 = \sigma^2 I_0(\theta_i, \theta_e)k^{1-2H} \qquad (4.4) $$

with

$$ I_0(\theta_i, \theta_e) = \begin{cases} \frac{C_{2H}}{2\pi^2} \cos^2 \theta_i |\sin \theta_i - \sin \theta_e|^{-2H-1} & \text{TE case} \\ \frac{C_{2H}}{2\pi^2} \frac{(1 - \sin \theta_i \sin \theta_e)^2}{\cos^2 \theta_e} |\sin \theta_i - \sin \theta_e|^{-2H-1} & \text{TM case.} \end{cases} $$

**Proof.** Since FBM is a zero-mean process, we obviously have $\mathbb{E}\hat{f}P_L(\alpha) = 0$ for all $\alpha$. This entails (4.3). For the second moment we have, away from the specular direction,

$$ \mathbb{E}|r(\alpha_i, \alpha_e)|^2 = \frac{k^2}{4\pi^2} |A_{\text{SPM}}(\theta_i, \theta_e)|^2 \int d\alpha_1 d\alpha_2 \\ \times \hat{P}_L(\alpha_1 - \alpha_i) \hat{P}_L(\alpha_2 - \alpha_i) \mathbb{E}\hat{f}(\alpha_e - \alpha_1) \bar{\hat{f}}(\alpha_e - \alpha_1). $$

Now

$$ \mathbb{E}\hat{f}(\alpha_e - \alpha_1)\bar{\hat{f}}(\alpha_e - \alpha_1) = \hat{C}(\alpha_e - \alpha_1, \alpha_2 - \alpha_e), $$

where $\hat{C}(\xi, \eta)$ is the bi-dimensional Fourier transform of the covariance function $C(x_1, x_2)$ of FBM. Outside the specular direction, this leads to

$$ \begin{aligned} \mathbb{E}|r(\alpha_i, \alpha_e)|^2 &= C_{2H} \frac{\sigma^2}{8\pi^2} k^2 |A_{\text{SPM}}(\theta_i, \theta_e)|^2 \int d\alpha \hat{P}_L(\alpha - \alpha_i)^2 |\alpha - \alpha_e|^{-1-2H} \\ &= C_{2H} \frac{\sigma^2}{8\pi^2} k^2 |A_{\text{SPM}}(\theta_i, \theta_e)|^2 \frac{L}{\sqrt{2\pi}} \int d\alpha e^{-(\alpha-\alpha_i)^2/2L^2} |\alpha - \alpha_e|^{-1-2H} \end{aligned} $$

which, for $L|\alpha - \alpha_i| \gg 1$, yields (4.4). $\square$
## 5. Building an estimator

According to the KA, SPM as well as SSA, the scattered intensity has a power-law dependence on the frequency $k$,

$$ \mathbb{E}|r(\alpha_i, \alpha_e)|^2 \sim k^\nu, \qquad (5.1) $$

with $\nu = 1 - 2H$. The most natural way of computing the Hurst exponent $H$ is therefore to work at fixed angles of incidence $\theta_i$ and emergence $\theta_e$ and to perform a log-log regression on $\mathcal{I}_1(k) := |r(\alpha_i, \alpha_e)|^2$. Let $k_1, \dots, k_N$ be the available wavenumbers and set

$$ x_j = \log k_j, \quad y_j = \log \mathcal{I}_1(k_j), \quad S_x = \sum x_j, \quad S_{xx} = \sum x_j^2. \qquad (5.2) $$

Then the first estimator is written

$$ \hat{H}_1 = \frac{1 - \hat{\nu}}{2}, \quad \text{with} \quad \hat{\nu} = \sum_{j=1}^{N} w_j y_j \quad \text{and} \quad w_j = \frac{N x_j - S_x}{NS_{xx} - S_x^2}. \qquad (5.3) $$

(Note that the weights satisfy $\sum w_j = 0$ and $\sum 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.

**Hypothesis 5.1 (SPM exactness).** *The SPM is exact.*

In that case, the estimator $\hat{H}_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 $\hat{H}_1$.

First we need to know the statistical nature of the diffracted intensity $\mathcal{I}_1(k)$.

**Proposition 5.2.** *Under hypothesis 5.1 and away from the specular direction, $\mathcal{I}_1(k)$ is an exponential variable.*

**Proof.** With the SPM, it follows from equation (2.7) that $|r(\alpha_i, \alpha_e)|^2$ is the Fourier transform of a Gaussian variable, hence a complex Gaussian variable. Therefore $|r(\alpha_i, \alpha_e)|^2$ is a chi-squared variate with two degrees of freedom, that is simply an exponential variable (i.e. a variable $X$ such that $P(X \ge t) = e^{-\mu t}, t \ge 0$). $\square$

Next we need to estimate the decorrelation of intensities at different intensities. Define

$$ \operatorname{Cor}(\alpha_i, \alpha_e; \alpha'_i, \alpha'_e) = \frac{\mathbb{E}|r(\alpha_i, \alpha_e)|^2 |r(\alpha'_i, \alpha'_e)|^2 - \mathbb{E}|r(\alpha_i, \alpha_e)|^2 \mathbb{E}|r(\alpha'_i, \alpha'_e)|^2}{\mathbb{E}|r(\alpha_i, \alpha_e)|^2 \mathbb{E}|r(\alpha'_i, \alpha'_e)|^2} $$

as the degree of correlation between the scattered intensities. Then the following result holds.

**Proposition 5.3.** *Under hypotheses 5.1 and 2.1, we have*

$$ \operatorname{Cor}(\alpha_i, \alpha_e; \alpha'_i, \alpha'_e) = e^{-\frac{L^2}{2}[(\alpha_e - \alpha_i) - (\alpha'_e - \alpha'_i)]^2} + e^{-\frac{L^2}{2}[(\alpha_e - \alpha_i) + (\alpha'_e - \alpha'_i)]^2}. \qquad (5.4) $$

This result essentially states that $|r(\alpha_i, \alpha_e)|^2$ and $|r(\alpha'_i, \alpha'_e)|^2$ are decorrelated as soon as $|\Delta\alpha \pm \Delta\alpha'| \gg 1/L$, where $\Delta\alpha = \alpha_e - \alpha_i$ is the momentum transfer. In particular, for a fixed geometry $(\theta_i, \theta_e)$, this means that $\mathcal{I}_1(k)$ and $\mathcal{I}_1(k')$ can be considered as uncorrelated whenever $L|k - k'| \gg 1$.
**Proof.** For notational convenience set $X_L(\alpha) = \int e^{-i\alpha x} f(x) P_L(x) dx$ and $\alpha = \alpha_e - \alpha_i$, $\alpha' = \alpha'_e - \alpha'_i$. In view of (2.7) it suffices to estimate the correlation term

$$ \mathbb{E}|X_L(\alpha)X_L(\alpha')|^2 = \int e^{-i\alpha(x-y)}e^{-i\alpha'(x'-y')}\mathbb{E}[f(x)f(y)f(x')f(y')] \\ \times 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:

$$ C(x_1, \dots, x_n) = \mathbb{E}(f(x_1) \dots f(x_n)) = \sum \prod C(u_l, u_m) \quad (5.5) $$

where the sum is taken over all subdivisions of the set $(x_1, \dots, 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

$$ \mathbb{E}|X_L(\alpha)X_L(\alpha')|^2 = C\hat{P}_L(\alpha, -\alpha)\widehat{CP}_L(\alpha', -\alpha') \\ +\widehat{CP}_L(\alpha, \alpha')\widehat{CP}_L(-\alpha, -\alpha') + \widehat{CP}_L(\alpha, -\alpha')\widehat{CP}_L(-\alpha, \alpha'). $$

Now, for all $\alpha_1, \alpha_2$,

$$ \widehat{CP}_L(\alpha_1, \alpha_2) = \int \hat{C}(\alpha_1 - \alpha'_1, \alpha_2 - \alpha'_2) \hat{P}_L(\alpha'_1) \hat{P}_L(\alpha'_2) d\alpha'_1 d\alpha'_2 = T_1 + T_2 + T_3, $$

where by (3.3),

$$ \begin{aligned} T_1 &= -\frac{\sigma^2}{2} C_{2H} \hat{P}_L(\alpha_2) \int |\alpha_1 - \alpha'_1|^{-2H-1} \hat{P}_L(\alpha'_1) d\alpha'_1 \\ T_2 &= -\frac{\sigma^2}{2} C_{2H} \hat{P}_L(\alpha_1) \int |\alpha_2 - \alpha'_2|^{-2H-1} \hat{P}_L(\alpha'_2) d\alpha'_2 \\ T_3 &= +\frac{\sigma^2}{2} C_{2H} \int |\alpha_1 - \alpha'_1|^{-2H-1} \hat{P}_L(\alpha'_1) \hat{P}_L(\alpha_1 + \alpha_2 - \alpha'_1) d\alpha'_1. \end{aligned} $$

In the regime $\alpha_j L \gg 1$, $j = 1, 2$, which is ensured by hypothesis 2.1, we have $L^{-1/2} \hat{P}_L(\alpha_j) \simeq 0$, $L^{1/2} \hat{P}_L(\alpha_j - \alpha) \simeq (2/\pi)^{1/4}\delta(\alpha - \alpha_j)$ and $\hat{P}_L(\alpha) \hat{P}_L(\alpha_1 + \alpha_2 - \alpha) \simeq \sqrt{2}\delta(\alpha) \exp(-\frac{L^2}{4}(\alpha_1 + \alpha_2)^2)$. Hence $T_1 = T_2 = 0$ and it remains

$$ \widehat{CP}_L(\alpha_1, \alpha_2) = \frac{\sigma^2 C_{2H}}{\sqrt{2}} e^{-\frac{L^2}{4}(\alpha_1 + \alpha_2)^2} |\alpha_1|^{-2H-1} = e^{-\frac{L^2}{4}(\alpha_1 + \alpha_2)^2} \widehat{CP}_L(\alpha_1, -\alpha_1). $$

Noting that

$$ \mathbb{E}|X_L(\alpha)|^2 \mathbb{E}|X_L(\alpha')|^2 = \widehat{CP}_L(\alpha, -\alpha)\widehat{CP}_L(\alpha', -\alpha') $$

the conclusion follows.

<figure><img src="image_6.png" /></figure>

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| \ge 1/L$ are perfectly decorrelated, i.e. $\text{Cor}(\alpha_i, \alpha_e; \alpha'_i, \alpha'_e) = 0$.

Under this assumption it is possible to evaluate the performances of the estimator $\hat{H}_1$.

**Theorem 5.5.** For a fixed geometry $(\theta_i, \theta_e)$, let $k^- = k_1 < \dots < k_N = k^+$ be the different available wavenumbers with $|k_j - k_{j'}| > 1/L$, $j \ne j'$. Suppose the wavenumbers are either in arithmetical (i.e. $k_{j+1} - k_j = \text{cst}$) or geometrical progression (i.e. $\log k_{j+1} - \log k_j = \text{cst}$).
Then under hypotheses 2.1, 3.1, 5.1 and 5.4, the estimator $\hat{H}_1$ is unbiased with asymptotic variance of order $1/N$:

$$ \mathbb{E}\hat{H}_1 = H \qquad (5.6) $$

$$ \mathrm{Var}(\hat{H}_1) \approx C/N, \quad N \to \infty \qquad (5.7) $$

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 $\mathcal{I}_1(k) = \sigma^2/2I_0(\theta_1, \theta_e)k^\nu 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 $\mathbb{E}y_j = \nu \log k_j + \text{constant}$ and $\mathrm{Var}(y_j) = \zeta(2, 1)$. Hence,

$$ \mathbb{E}\hat{v} = \sum w_j \mathbb{E}y_j = \nu \sum w_j x_j = \nu, \\ \mathrm{Var}(\hat{v}) = \sum w_j^2 \mathrm{Var}(y_j) = \frac{\zeta(2, 1)}{\ln^2(10)} \sum w_j^2. $$

It follows at once that $\hat{H}_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 $k_j = k_1 + (k_N - k_1)j/N$) or a regular logarithmic sampling (that is $\log k_j = \log(k_1) + (\log k_N - \log k_1)j/N$). An elementary but tedious calculation shows that $\lim_{N \to \infty} N \sum w_j^2 = C$ with

$$ C^{-1} = C_1^{-1} := \frac{1}{k^+ - k^-} \left[ \left( \int_{k^-}^{k^+} \log^2(k) dk \right) - \frac{1}{k^+ - k^-} \left( \int_{k^-}^{k^+} \log(k) dk \right)^2 \right] $$

for linear sampling and $C^{-1} = C_2^{-1} := \frac{1}{12} \log^2(k^+/k^-)$ for logarithmic sampling. Hence,

$$ \mathrm{Var}(\hat{H}_1) = \frac{\mathrm{Var}(\nu)}{4} \approx \frac{\zeta(2, 1)}{4 \ln^2(10)} \frac{C}{N} = 0.14 \frac{C}{N} \qquad (5.8) $$

for large N. Although we have not been able to prove it analytically, a numerical estimation clearly shows that $C_2 \le C_1$ for all $(k^-, k^+)$, so that a further variance reduction can be achieved by choosing a logarithmic sampling. $\square$

The variance of $\hat{H}_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 $\mathrm{Var}(\hat{H}_1)$ is $O(1/L)$.

## 6. An improved estimator

The main drawback of the estimator $\hat{H}_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 < \dots < k_N$ with a regular sampling on a logarithmic scale. For a fixed frequency $k$ and incidence angle $\theta_i$, we can achieve several values of the momentum transfer $\alpha - \alpha_i$ by varying the scattering angle $\theta$ while keeping, however, the same order of magnitude for the latter. Let $\theta$ vary within a cone of angular aperture $(\theta_e^-, \theta_e^+)$ around some mean non-specular direction $\theta_e$. Now define the angular mean:

$$ \mathcal{I}_2(k) = \frac{1}{N_\theta} \sum_{j=1}^{N_\theta} |r(\alpha_i, \alpha_j)|^2, $$

where $\alpha_j = k \sin \theta_j$, $\theta_j = \theta_e^- + j \Delta\theta$ and $\Delta\theta = (\theta_e^+ - \theta_e^-)/N_\theta$. The angular sampling rate $\Delta\theta$ is taken to be the smallest that allows us to keep the successive values of $|r(\alpha_i, \alpha_j)|^2$ uncorrelated at all wavenumbers. This is the case as soon as $k^{-}(\sin(\theta_e + \Delta\theta) - \sin\theta_e) \gg 1/L$, that is roughly as soon as the criterion

$$ \Delta\theta \gg \frac{1}{k^- L \cos \theta_e} \qquad (6.1) $$

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

$$ \mathbb{E} \mathcal{I}_2(k) = \frac{\sigma^2}{N_\theta} \sum_{j=1}^{N_\theta} I(\theta_i, \theta_j) k^\nu, $$

where the expression of the angular factor $I(\theta_i, \theta_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.

**Hypothesis 6.1 (Constancy of the angular factor).** The angular factor $I(\theta_i, \theta)$ remains constant within the prescribed emergence cone $\theta_e^- \le \theta \le \theta_e^+$.

This, together with the assumption of hypothesis 5.4, ensures that $\mathcal{I}_2(k) = \frac{\sigma^2}{2N_\theta} I_0(\theta_i, \theta_e)k^\nu Z$, where $Z$ is a chi-square variate with $2N_\theta$ degrees of freedom. Again, we need only perform a log-log regression on $\mathcal{I}_2(k)$. The corresponding estimator $\hat{H}_2$ is then given by equations (5.2) and (5.3), where $\mathcal{I}_1(k)$ is to be replaced by $\mathcal{I}_2(k)$. The point is that $\mathcal{I}_2(k)$ has a much smaller variance than $\mathcal{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 $\hat{H}_2$ is unbiased and, for the same choice of wavenumbers ($k_j$), its variance is reduced by a factor $N_\theta$ with respect to $\hat{H}_1$.

**Proof.** From appendix we now have $\mathbb{E}y_j = \nu \log k_j + \text{cst}$ and $\text{Var}(y_j) = \zeta(2, N_\theta/2)$. Hence,

$$ \begin{aligned} \mathbb{E}\nu &= \sum w_j \mathbb{E}y_j = \nu \sum w_j x_j = \nu, \\ \text{Var}(\nu) &= \sum w_j^2 \text{Var}(y_j) = \zeta(2, N_\theta) \sum w_j^2. \end{aligned} $$

Again $\hat{H}_2$ is an unbiased estimator of $H$. Moreover, comparison with the previously calculated variance shows that

$$ \frac{\text{Var}(\hat{H}_2)}{\text{Var}(\hat{H}_1)} = \frac{\zeta(2, N_\theta)}{\zeta(2, 1)} = \mathcal{O}\left(\frac{1}{N_\theta}\right). \qquad (6.2) $$

This completes the proof. $\square$
<table><caption>**Table 1.** Performances of the estimators $\hat{H}_1$ and $\hat{H}_2$ over 268 realizations. $\mu$ and $s$ are the theoretical mean and standard deviation; $\bar{\mu}$ and $\bar{s}$ are the corresponding empirical quantities.</caption><thead><tr><th>Estimator</th><th>μ</th><th>μ̄</th><th>s</th><th>s̄</th></tr></thead><tbody><tr><td>\hat{H}_1 (SPM)</td><td>0.7</td><td>0.724</td><td>0.125</td><td>0.143</td></tr><tr><td>\hat{H}_1 (KA)</td><td></td><td>0.681</td><td></td><td>0.089</td></tr><tr><td>\hat{H}_2 (SPM)</td><td>0.7</td><td>0.689</td><td>0.040</td><td>0.041</td></tr><tr><td>\hat{H}_2 (KA)</td><td></td><td>0.681</td><td></td><td>0.034</td></tr></tbody></table>

<table><caption>**Table 2.** Relevant scale parameters of the FBM surface for $H = 0.5$ and $0.7$ and for different values of $\sigma$. The illuminated patch has characteristic half-size $L = 256$ and the incident wavelength ranges from 2 to 16.</caption><thead><tr><th>H</th><th>σ</th><th>σL<sup>H</sup></th><th>σL<sup>H-1</sup></th><th>σλ<sup>H-1</sup> (SSA)</th><th>Plot</th></tr></thead><tbody><tr><td rowspan="3">0.5</td><td>0.62</td><td>10</td><td>0.04</td><td>0.2–0.4</td><td>H05small</td></tr><tr><td>1.25</td><td>20</td><td>0.08</td><td>0.3–0.9</td><td>H05medium</td></tr><tr><td>2.50</td><td>40</td><td>0.16</td><td>0.6–1.8</td><td>H05large</td></tr><tr><td rowspan="3">0.7</td><td>0.31</td><td>15</td><td>0.06</td><td>0.1–0.3</td><td>H07small</td></tr><tr><td>1.24</td><td>60</td><td>0.23</td><td>0.5–1.0</td><td>H07medium</td></tr><tr><td>2.47</td><td>120</td><td>0.47</td><td>1.1–2.0</td><td>H07large</td></tr></tbody></table>

*Confidence intervals.* Once the variance Var($\hat{H}$) = $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]$.

## 7. 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 $\sigma = 0.05$ and let the incident wavelength $\lambda$ range from 0.1 to 16. We have chosen the backscattering direction with $\theta_i = 25$. With these values, hypotheses 2.1, 3.1 and 3.3 are clearly verified with $L(\alpha_e - \alpha_i) \in [85, 13632]$, $\sigma L^{H-1} = 0.01$, $\sigma\lambda^{H-1} \in [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 \theta_e \Delta\theta$ 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.

## 8. 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
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Log-log diagram of the mean backscattered intensity $k^{-2}\mathbb{E}\mathcal{I}_2(k)$ for an incidence $\theta_i = 20^\circ$ (the factor $k^{-2}$ is for better visibility). Theoretical value of the Hurst exponent $H = 0.5$.</figcaption></figure>

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$, $\sigma$ and $\lambda$, and for a fixed incidence $\theta_i = 20^\circ$; in each case averages over 20 realizations of the surface have been performed in order to obtain the mean intensity $\mathbb{E}|r(\alpha_i, \alpha_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 $\sigma$ 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 $\sigma$. Figures 5 and 6 show the mean backscattered intensity $\mathbb{E}\mathcal{I}_2(k) = \mathbb{E}|r(\alpha_i, -\alpha_i)|^2$ in a log-log diagram for the different values of $H$ and $\sigma$, at incidence $\theta_i = 20^\circ$. A further angular average has been performed within a $5^\circ$ 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 $\nu$. The different values of $\mu = (1 - \nu)/2$ are reported in table 3. From these results we may draw several conclusions:
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Same as figure 5 with $H = 0.7$.</figcaption></figure>

<table><caption>**Table 3.** Values of the slope $\mu$ estimated from a log-log regression on the mean backscattered intensity (figures 5 and 6).</caption><thead><tr><th>H</th><th>σ</th><th>μ (exact)</th><th>μ (SSA)</th></tr></thead><tbody><tr><td rowspan="3">0.5</td><td>0.62</td><td>0.495</td><td>0.513</td></tr><tr><td>1.25</td><td>0.553</td><td>0.729</td></tr><tr><td>2.50</td><td>0.890</td><td>1.120</td></tr><tr><td rowspan="3">0.7</td><td>0.31</td><td>0.675</td><td>0.717</td></tr><tr><td>1.24</td><td>0.698</td><td>0.732</td></tr><tr><td>2.47</td><td>1.002</td><td>0.941</td></tr></tbody></table>

* A non-trivial power-law is observed for the mean backscattered intensity far beyond the SSA regime, that is $EI_2(k) \sim k^{-\nu}$.

* The corresponding exponent $\mu$ increases with the height parameter $\sigma$, starting from the value predicted by the SSA: $\mu = H$.

* The exact solution starts departing from the SSA for $\sigma L^{H-1} \ge 0.5$. For $\sigma L^{H-1} \le 0.5$, both the exact solution and the SSA give a good estimate of the Hurst exponent, namely $\mu = H$.

## 9. 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.

## 10. Appendix: some useful formulae

*Homogeneous distributions.* For $\alpha < -1$, the distribution $\kappa_{\alpha}(x) = |x|^{\alpha}$, which is a priori only defined for $x \neq 0$, can be extended over $\mathbb{R}$ (provided $\alpha$ is not a negative integer) by setting [12]

$$ \kappa_{\alpha}(\varphi) = \int_{-\infty}^{+\infty} |x|^{\alpha} \left[ \varphi(x) - \varphi(0) - x\varphi'(0) - \dots - \frac{x^{(n-1)}}{(n-1)!} \varphi^{(n-1)}(0) \right] dx, $$

with $-n-1 < \text{Re} \alpha < -n$. Its Fourier transform is given by

$$ \hat{\kappa}_{\alpha}(\xi) = C_{\alpha} \kappa_{-\alpha-1}(\xi), \qquad (10.1) $$

with $C_{\alpha} = -2\Gamma(\alpha + 1) \sin(\alpha\pi/2)$. The bi-dimensional Fourier transform of $\kappa_{\alpha}(x_1, x_2) = |x_1 - x_2|^{\alpha}$ is

$$ \hat{\kappa}_{\alpha}(\xi, \eta) = -C_{\alpha} \delta(\xi + \eta) \kappa_{-\alpha-1}(\xi). \qquad (10.2) $$

*Oscillating integrals.*

**Lemma 10.1.** Let

$$ A_a(\xi) = \int e^{i\xi x} e^{-\frac{|x|^\nu}{2}} e^{-\frac{\xi^2}{a^2}x^2} dx \\ I_a(\xi) = \int e^{i\xi(x_1-x_2)} e^{-\frac{|x_1-x_2|^\nu}{2}} e^{-\frac{\xi^2}{a^2}x_1^2} e^{-\frac{\xi^2}{a^2}x_2^2} dx_1 dx_2. $$

Then for large $a$ we have

$$ A_a(\xi) \approx \frac{\pi^{1/2}}{2} C_{\nu} \xi^{-\nu-1} + o(\xi^{-\nu-1}) \\ a I_a(\xi) \approx \frac{\pi^{1/2}}{2\sqrt{2}} C_{\nu} \xi^{-\nu-2} + o(\xi^{-\nu-2}). $$

**Proof.** Let us rewrite

$$ A_a(\xi) = \int e^{i\xi x} \left(1 - \frac{|x|^\nu}{2}\right) e^{-\frac{\xi^2}{a^2}x^2} dx + \int e^{i\xi x} \phi(x) e^{-\frac{\xi^2}{a^2}x^2} dx = T_1 + T_2 $$

with an obvious definition for $\phi$. For the first term we obtain, using the Fourier transform of homogenous distributions,

$$ T_1 = a\xi^{-1}e^{-a^2\xi^2} - \frac{C_{\nu}}{2}a\xi^{-1} \int e^{-\frac{\xi^2 a^2}{\xi^2} |\xi'|^{-\nu} - |\xi|^{-\nu-1}} d\xi' $$

which, in the limit $a \to \infty$, yields $T_1 = \pi^{1/2} \frac{C_{\nu}}{2} \xi^{-\nu-1}$. For the second term we clearly have $T_2 = \hat{\phi}(\xi)$ as $a$ goes to infinity. Now note that $\phi$ is a rapidly decreasing function of regularity
$C^{2\nu}$. Its Fourier transform is thus a bounded function that decays at least in $|\xi|^{-[2\nu]}$, $[2\nu]$ being the integer part of $2\nu$. If $[2\nu] > \nu + 1$, this proves the first statement; else if $[2\nu] \le \nu + 1$ we again decompose $\phi(x) = \frac{x^{2\nu}}{2!} + \cdots + (-1)^p \frac{x^{p\nu}}{p!} + \phi_1(x)$ with $[p\nu] > \nu + 1$, where $\phi_1(x)$ is now a rapidly decreasing function of regularity $C^{p\nu}$. The first $p$ terms are homogeneous distributions with a Fourier transform that is $O(\xi^{-\nu-2})$ at infinity while $\hat{\phi}_1$ decays at least in $|\xi|^{-[p\nu]}$. This completes the proof of the first statement. The proof of the second statement is similar (using the bi-dimensional Fourier transform of the distribution $|x_1 - x_2|^\nu)$, so we omit it. $\square$

*Some statistical results.* If $X_1, \dots, X_n$ are i.i.d standard Gaussian variables, then $Z_n = X_1^2 + X_2^2 + \dots + X_n^2$ is a chi-squared variate with $n$ degrees of freedom. Its density $f_n$ is given by $f_n(x) = (2^{n/2}\Gamma(n/2))^{-1}x^{n/2-1}e^{-x/2}$ for $x \ge 0$. It can be shown [24] that

$$ \mathbb{E} \ln Z_n = \psi(n/2) + \ln 2 \\ \mathrm{Var}(\ln Z_n) = \zeta(2, n/2) \approx \frac{2}{n}, \quad n \to \infty, $$

where $\psi(z) = \Gamma'(z)/\Gamma(z)$ is the psi function and $\zeta(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.

## References

[1] Adler R 1981 *The Geometry of Random Fields* (New York: Wiley)

[2] Agnon Y and Stiassnie M 1991 Remote sensing of the roughness of a fractal sea surface *J. Geophys. Res.* **96** 12773–9

[3] Berizzi F and Dalle-Mese E 1999 Fractal analysis of the signal scattered from the sea surface *IEEE Trans. Antennas Propagat.* **47** 324–38

[4] Berry M V 1981 Diffractal echoes *J. Phys. A: Math. Gen.* **14** 3101–10

[5] Chen J, Lo T K L, Leung H and Litva J 1996 The use of fractals for modeling em waves scattering from rough sea surfaces *IEEE Trans. Geosci. Remote Sens.* **34** 966–72

[6] Elfouhaily T, Chapron B, Katsaros K and Vandemark D 1997 A unified directional spectrum for long and short wind-driven waves *J. Geophys. Res.* **102** 15781–96

[7] Flandrin P 1989 On the spectrum of fractional Brownian motion *IEEE Trans. Inform. Theory* **35** 197–9

[8] Franceschetti G, Iodice A, Migliaccio M and Riccio D 1996 An electromagnetic fractal-based model for the study of fading *Radio Sci.* **31** 1749–59

[9] Franceschetti G, Iodice A, Maddaluno S and Riccio D 2000 A fractal-based theoretical framework for retrieval of surface parameters from electromagnetic backscattering data *IEEE Trans. Geosci. Remote Sens.* **38** 641–50

[10] Franceschetti G, Iodice A, Migliaccio M and Riccio D 1999 Fractals and the small perturbation scattering model *Radio Sci.* **34** 1043–54

[11] Franceschetti G, Iodice A, Migliaccio M and Riccio D 1999 Scattering from natural rough surfaces modeled by fractional Brownian motion *IEEE Trans. Antennas Propagat.* **47** 1405–415

[12] Gelfand I M and Shilov G E 1964 *Generalized Functions* vol 1 (New York: Academic) (in Russian) (Engl. Transl. 1977 *Properties and Operations* transl E Saletan))

[13] Guérin C A, Holschneider M and Saillard M 1997 Electromagnetic scattering on multi-scales rough surfaces *Waves Random Media* **7** 331–49

[14] Jaggard D L and Sun X 1990 Scattering from fractally corrugated surfaces *J. Opt. Soc. Am.* **7** 1131–9

[15] Leu J S and Papamarcou A 1995 On estimating the spectral exponent of fractional Brownian motion *IEEE Trans. Inform. Theory* **41** 233–44

[16] Mandelbrot B B and Van Ness J W 1968 Fractional Brownian motions, fractional noises and applications *SIAM Rev.* **10** 422–37

[17] Mattia F 1999 Backstattering properties of multiscale rough surfaces *J. Electromagn. Waves Appl.* **13** 493–527

[18] Maystre D, Mendez O M and Roger A 1983 A new electromagnetic theory for scattering from shallow rough surfaces *Opt. Acta* **30** 1707–23

[19] Rouvier P, Borderies S and Chenerie I 1997 Ultra-wide band electromagnetic scattering of a fractal profile *Radio Sci.* **32** 285
[20] Saillard M and Maystre D 1990 Scattering from metallic and dielectric rough surfaces *J. Opt. Soc. Am.* **7** 331–49

[21] Savaidis P, Frangos S, Jaggard D L and Hizanidis K 1997 Scattering from fractally corrugated surfaces with use of the extended boundary condition method *J. Opt. Soc. Am.* **14** 475–85

[22] Savaidis S, Frangos P, Jaggard D L and Hizanidis K 1995 Scattering from fractally corrugated surfaces: an exact approach *Opt. Lett.* **20** 2357–9

[23] Thorsos E I and Broschat S L 1995 An investigation of the small-slope approximation for scattering from rough surfaces part I *Theory Waves Random Media* **97** 2082–93

[24] Veitch D and Abry P 1998 A wavelet based joint estimator for the parameters of lrd *IEEE Trans. Inform. Theory*. *Multiscale Statistical Signal Analysis and its Application* special issue

[25] Voronovitch A G 1985 Small-slope approximation in wave scattering from rough surfaces *Sov. Phys.–JETP* **62** 65–70

[26] Voronovitch A G 1994 *Wave Scattering From Rough Surfaces* (Springer Series on Wave Phenomena) (Berlin: Springer)

[27] Wood A T A and Chan G 1994 Simulation of stationary gaussian processes in $[0, 1]^d$ *J. Comput. Stat.* **2** 409–32