# Inverse scattering with far-field intensity data: random surfaces that belong to a well-defined statistical class

Demetrio Macías, Gustavo Olague, Eugenio Méndez

► To cite this version:

Demetrio Macías, Gustavo Olague, Eugenio Méndez. Inverse scattering with far-field intensity data: random surfaces that belong to a well-defined statistical class. Waves in Random and Complex Media, 2006, 16 (4), pp.545-560. 10.1080/17455030600830385. hal-02372999

HAL Id: hal-02372999

https://utt.hal.science/hal-02372999

Submitted on 2 Mar 2021

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 Attribution 4.0 International License
# Inverse scattering with far-field intensity data: random surfaces that belong to a well-defined statistical class

DEMETRIO MACÍAS†, GUSTAVO OLAGUE‡ and EUGENIO R. MÉNDEZ*‡

†Institut Charles Delaunay, Université de technologie de Troyes, CNRS FRE 2848, Laboratoire de Nanotechnologie et d'Instrumentation Optique, 12 rue Marie Curie, BP-2060 F-10010 Troyes CEDEX, France

‡División de Física Aplicada, CICESE, Km 107 carretera Tijuana-Ensenada, Ensenada, 22860, Mexico

*(Received 26 October 2005; in final form 8 March 2006)*

We consider the inverse scattering problem of retrieving the surface profile function from far-field angle-resolved intensity data. The problem is approached as a nonlinear constrained optimization problem. The surface, assumed one-dimensional and perfectly conducting, is also assumed to be a realization of a Gaussian random process with a Gaussian correlation function with known standard deviation of heights (δ) and correlation length (a). Starting from rigorously calculated far-field angle-resolved scattered data, we search for the optimum profile using evolutionary strategies. Examples that illustrate the proposed scheme are presented. Aspects of the convergence and lack of uniqueness of the solution are discussed.

## 1. Introduction

The scattering of waves from randomly rough surfaces has a long history and is a subject of importance in many fields of physics. The direct scattering problem consists of finding the scattered field or the scattered intensity from knowledge of the surface profile, its optical properties, and the conditions of illumination. Considerable progress towards an understanding of the phenomenon has been made in the last few decades [1–4].

Rough surface inverse scattering problems have also received some attention, and take a variety of forms. Most of the studies have focused on the retrieval of some basic statistical parameter of the surface, such as the standard deviation of heights, or the height correlation length [5–9]. These methods are based on approximate solutions to the direct problem, like the Kirchhoff approximation [1, 2] or small-amplitude perturbation theory [2], and on simple models for the statistical properties of the surface. A novel approach, based on a reverse Monte Carlo simulation, was reported more recently for the estimation of the power spectrum of surfaces in situations in which multiple scattering is significant [10].

A different kind of inverse problem is that of designing a surface with specified scattering properties. This field has developed rapidly in recent years. Procedures to design and fabricate one-dimensional surfaces that produce rather arbitrary angular scattering distributions have been developed [11–14]. Extensions to the case of general two-dimensional angular

<footnote>*Corresponding author. E-mail: emendez@cicese.mx</footnote>
distributions are not simple, but procedures to generate such diffusers have been recently put forward [15, 16].

The problem of recovering surface profiles from complex amplitude far-field data has also been the subject of some studies [17–21]. An important limitation of these methods is that they are based on approximate models for the interaction between the incident light and the surface, and fail when multiple scattering is important. In addition, since optical detectors are not phase sensitive, the need to have amplitude rather than intensity data constitutes a drawback that these methods share.

In the present work, we investigate the possibility of recovering the surface profile function from far-field intensity data. We approach the problem as a problem of constrained optimization. In previous communications [22, 23] we studied the performance of rough surface inversion algorithms based on evolutionary strategies using two kinds of representations of the objective variables. In [22], we employed a spectral representation of the surface and studied the performance of two selection strategies; the elitist and non-elitist strategies. The search space was reduced by considering surfaces that belong to a well-defined statistical class. This constraint did not permit the use of intermediate recombination, as it would eventually lead to departures from the assumed statistics. To circumvent this problem, in [23] the surface was represented in terms of spline curves, permitting the use of intermediate recombination and the treatment of deterministic profiles. It was found that the use of intermediate recombination helped in the convergence of the algorithm.

In this paper, we consider the problem studied in [22] and explore the effects of using dominant recombination, or crossover, on the performance of the algorithm. Dominant recombination does not produce the kind of problems introduced by the intermediate recombination operator. In addition, we illustrate the robustness of the method presenting data corresponding to five different surfaces belonging to the same statistical class.

The methods described here, and in [22, 23], do not rely on approximate expressions for the field-surface interaction. Furthermore, our experience indicates that the convergence to the optimum improves when multiple scattering occurs. Since, in general, the use of intensity data implies that the solution to the problem is not unique, the improvement in the convergence is possibly due to a reduction in the number of solutions in the presence of multiple scattering.

The paper is organized as follows. Section 2 introduces the notation and presents a brief account of the direct scattering problem for one-dimensional surfaces. In Section 3, we formulate the inverse problem as an optimization problem, describing the evolutionary algorithm in some detail and discussing the differences between the intermediate and dominant recombination. Representative results are presented in Section 4 and, finally, in Section 5 we present our conclusions.

## 2. Scattering by one-dimensional rough surfaces

We consider the scattering of light from a one-dimensional, perfectly conducting, randomly rough surface defined by the equation $x_3 = \zeta(x_1)$. The region $\zeta(x_1) > x_3$ is the vacuum, the region $\zeta(x_1) < x_3$ is a perfect conductor, and the plane of incidence is the $x_1x_3$-plane. With reference to figure 1, the surface is illuminated from the vacuum side with a p- or s-polarized plane wave.

The single non-zero component of the electric or magnetic vector of the incident field has the form

$$ \psi_2(x_1, x_3|\omega)_{\text{inc}} = \psi_0 \exp\{i[kx_1 - \alpha_0(k)x_3]\}, \quad (1) $$

where $\alpha_0(k) = \sqrt{(\omega/c)^2 - k^2}$, $\omega$ is the frequency of the field, and $c$ is the speed of light in vacuum. A time dependence of the form $\exp(-i\omega t)$ is assumed, but explicit reference to it is suppressed.
<figure><img src="image_2.png" /><figcaption>Figure 1. Geometry of the scattering problem.</figcaption></figure>

The field scattered by the surface is

$$ \psi_2(x_1, x_3 | \omega)_{sc} = \int_{-\infty}^{\infty} \frac{dq}{2\pi} R_{p,s}(q | k) \exp\{iqx_1 + i\alpha_0(q)x_3\}, \quad (2) $$

where the scattering amplitude $R_{p,s}(q | k)$ can be written in the form [24]

$$ R_{p,s}(q | k) = \frac{i}{2\alpha_0(q)} \int_{-\infty}^{\infty} dx_1 \exp\{-iqx_1 - i\alpha_0(q)\zeta(x_1)\} \chi(k, \zeta'(x_1)) F(x_1 | \omega). \quad (3) $$

The factor $\chi(k, \zeta'(x_1))$ is

$$ \chi(k, \zeta') = \begin{cases} -i[k\zeta'(x_1) - \alpha_0(k)] & \text{in p polarization} \\ 1 & \text{in s polarization,} \end{cases} \quad (4) $$

and $F(x_1 | \omega)$ represents a source function, given by

$$ F(x_1 | \omega) = \begin{cases} \psi_2(x_1, x_3) & \text{in p polarization} \\ [\partial\psi_2(x_1, x_3)/\partial N]_{x_3 = \zeta(x_1)} & \text{in s polarization,} \end{cases} \quad (5) $$

where $\partial/\partial N = N \cdot \nabla = [-\zeta'(x_1)(\partial/\partial x_1) + (\partial/\partial x_3)]$ is the non-unit normal derivative operator. The angles of incidence $\theta_0$ and scattering $\theta_s$ are related to the components of the wavenumbers $k$ and $q$ that are parallel to the mean surface through the expressions

$$ k = \frac{\omega}{c} \sin \theta_0, \quad q = \frac{\omega}{c} \sin \theta_s. \quad (6) $$

The far-field intensity $I_{p,s}(q | k)$ is defined as

$$ I_{p,s}(q | k) = |R_{p,s}(q | k)|^2. \quad (7) $$

In the present work, the goal is to retrieve the unknown surface profile function $\zeta(x_1)$ from the intensity data $I_{p,s}(q | k)$. With the squaring operation of the scattering amplitude $R_{p,s}(q | k)$ the phase information is lost, and this adds to the difficulties of inverting the data. Among other things, the solution to the problem is not unique. That is, there can be two or more surfaces that give rise to the same intensity scattering pattern. Aspects of the lack of uniqueness of the solution will be discussed in Section 4.

### 3. Inverse scattering as a least-squares approximation problem

From the previous discussion, it is clear that establishing an inversion scheme to retrieve the profile of a rough surface from scattered intensity data is not an easy task. In this section we
reformulate the problem of inverse scattering in terms of a nonlinear least-squares bounds-constrained approximation problem. It is assumed that we have access to far-field angle-resolved scattered intensity data corresponding to several angles of incidence. The goal is to retrieve the unknown surface profile function from these data. Some constraints on the kind of surface that we seek are introduced in order to reduce the search space.

## 3.1. Definition of the fitness functional

The inverse scattering problem can be reformulated in terms of the fitness (objective) functional

$$f(z(x_1)) = \min \left\{ \sum_{i=1}^{N_{\text{ang}}} \| I^{(m)}(q | k_i) - I^{(c)}(q | k_i; z(x_1)) \|_2^2 \right\}, \quad (8)$$

where the symbol $\| \cdot \|_2$ represents the Euclidean norm of the intensity vectors as functions of the scattering wavevector components $q$, $N_{\text{ang}}$ is the number of angles of incidence considered, and the $k'_i$s are related to those angles through equation (6). Then, $I^{(m)}(q | k_i)$ represents an angle-resolved far-field scattered intensity pattern of the surface of interest (measured or calculated) and $I^{(c)}(q | k_i; z(x_1))$ is a calculated intensity pattern obtained by solving the direct problem with a trial surface profile $z(x_1)$. The functional $f(z(x_1))$ can be interpreted as an assessment of the closeness between the angular distributions of intensity $I^{(m)}(q | k)$ and $I^{(c)}(q | k; z(x_1))$. The goal would then be to find a surface for which the condition $I^{(c)}(q | k) = I^{(m)}(q | k)$ is satisfied. When this happens, and if the solution to the problem is unique, the best approximation to the original profile has been retrieved.

Note that in our definition of the fitness functional we require that the proposed surface reproduces the 'measured' scattering data for several angles of incidence. The satisfaction of this constraint should reduce the number of possible solutions and, hopefully, produce a unique one. The inverse scattering problem can be viewed, now, as the problem of minimizing $f(\zeta^c(x_1))$.

## 3.2. Representation of the objective variables

To deal with the scattering problem numerically, the surface must be sampled. From the preceding discussion it seems natural to choose, as the parameters of interest, the surface heights evaluated at the sampling points. Changing these numbers independently, however, would lead to surfaces with abrupt height changes, which does not correspond to the physical situation of interest. One way to avoid this problem is to restrict the search space to randomly rough surfaces that belong to a certain class. We are, thus, faced with a problem of constrained optimization.

In our case, we have chosen the target surface as a realization of a stationary, zero-mean one-dimensional Gaussian random process. With this assumption, the random process is completely characterized by its two-point correlation function, which we also assume to be Gaussian:

$$\langle \zeta(x_1) \zeta(x'_1) \rangle = \delta^2 \exp\{(x_1 - x'_1)^2 / a^2\}. \quad (9)$$

Here, the angle brackets denote an average over the ensemble of realizations of the surface profile function, $\delta = \langle \zeta^2(x_1) \rangle^{1/2}$ is the rms height of the surface, and $a$ its correlation length.
Surfaces belonging to this class can be generated numerically with the spectral method described in [24]. Correlated random numbers that represent the surface heights at the sampling points can be obtained through the expression

$$ \zeta_n = \frac{\delta}{\sqrt{L}} \sum_{j = -N/2}^{N/2 - 1} \frac{[M_j + iN_j]}{\sqrt{2}} \left[ \sqrt{\pi} a \exp \left\{ - \left( \frac{aq_j}{2} \right)^2 \right\} \right]^{1/2} \exp\{iq_j \chi_n\}. \quad (10) $$

Here, N represents the total number of points on the surface, L represents its length, $\chi_n = -L/2 + (n - 0.5)\Delta x$ are the sampling points spaced by $\Delta x$ along $x_1$, $q_j = -\pi/\Delta x + 2\pi(j - 0.5)/L$ are the sampling points in Fourier space, and $\zeta_n = \zeta(\chi_n)$. The random sets $\{M_j\}$ and $\{N_j\}$ contain statistically independent random Gaussian variables with zero mean and unit standard deviation. In order to produce a set of real random numbers $\{\zeta_n\}$, it is required that the complex array $\{M_j + iN_j\}$ be Hermitian. The first and second derivatives of the surface profile function, which are required for the direct scattering calculations, can be obtained by differentiation of equation (10).

### 3.3. Evolutionary inversion procedure

At least in principle, any of the optimization techniques reported in the literature could be employed to minimize equation (8). However, as discussed in [22], the form of equation (8) and the constraint imposed by the representation scheme of equation (10) suggest the use of an algorithm belonging to the class of ‘direct search methods’ [25]. The main characteristic of this kind of technique is that, throughout the entire optimization process, one only needs to know the values of the fitness function and not its derivatives [26–28].

Evolutionary algorithms are a relatively recent set of ‘direct search methods’ that have been successful in the solution of ill-posed inverse problems in different scientific disciplines [29,30–32]. Examples of these heuristic population-based techniques are the genetic algorithms [33], the evolution strategies [34], and genetic [35] and evolutionary programming [36]. Notwithstanding their differences, all evolutionary algorithms are inspired on the Darwinian principles of variation and selection [37]. Given the characteristics of the inverse problem studied here, it was considered that the evolutionary strategies were the best-suited evolutionary algorithms for this task.

Several variations of evolutionary strategies have been proposed [38]. All of them, however, follow the canonical structure shown in the flow diagram of figure 2. The starting point of the optimization process is the generation of a random set $P_{\mu}^{(g)}|_{g=0}$ of $\mu$ possible solutions to the problem which, in the present context, are the Gaussian randomly rough one-dimensional surfaces $\{z_n\}$ generated through equation (10). A secondary population $P_{\lambda}^{(g)}$ of $\lambda$ elements is generated through the application of the ‘genetic’ operations of recombination and mutation over the elements of the initial population $P_{\mu}^{(g)}$. This represents the start of the main evolutionary loop.

At this point, it is pertinent to mention that we have previously used $\lambda$ to denote the wavelength of the light, which is the usual notation in optical work. It is believed that due to the different context in which the two quantities are employed, use of the same symbol to denote both should not lead to much confusion.

Recombination exploits the search space through the exchange of information between different elements of the population. Mutation, on the other hand, explores the search space through the introduction of random variations in the newly recombined elements. Depending on the problem studied, it is possible to exclude the recombination operation from the evolutionary loop, as it is indicated by the broken line in figure 2. That is, in the search for
<figure><img src="image_3.png" /><figcaption>Figure 2. Flow diagram of the evolutionary loop.</figcaption></figure>

the optimum, mutation could be the only operator employed. Once the secondary population $P_{\lambda}^{(g)}$ has been generated, one needs to evaluate the quality of its elements. For this, the direct problem must be solved for each one of the newly generated surfaces of the secondary population. With this, a fitness value is associated to each trial surface. This is done through the comparison between the calculated scattering pattern $I^{(c)}(q | k_i; z_n)$ and the ‘measured’ data $I^{(m)}(q | k_i)$, on the basis of equation (8). Only those elements of the secondary population $P_{\lambda}^{(g)}$ leading to promising regions of the search space will be retained, through some selection scheme, as part of the population $P_{\mu}^{(g+1)}$ for the next iteration of the evolutionary loop. The procedure is repeated until a defined termination criterion has been achieved. The respective sizes of the initial and the secondary populations, $P_{\mu}^{(g)}$ and $P_{\lambda}^{(g)}$, remain constant throughout the entire search process.

### 3.3.1 Recombination and mutation.
In the search for the optimum, new trial surfaces are generated from the initial set of surfaces (population). Surfaces can be generated by selectively
<figure><img src="image_4.png" /><figcaption>Figure 3. Schematic representation of the mutation operation.</figcaption></figure>

combining them (intermediate and dominant recombination), introducing random changes in them (mutation), or both.

In the intermediate recombination, new surfaces are generated through linear combinations of the parent population. In our case, the first time that this operator is used the new surfaces are the result of linear combinations of independent Gaussian random processes. The result is also a Gaussian process and one can adjust the parameters of the combination in such a way that the new surfaces belong to the required statistical class. In successive iterations, however, the surfaces become more and more correlated, and the same operation generates surfaces that do not belong to the assumed statistical class. Intermediate recombination is, thus, not compatible with the representation and constraints imposed on the problem.

In the dominant recombination, or crossover, the surfaces exchange elements of the Hermitian arrays (see equation 10). The resulting surfaces are Gaussian and belong to the same statistical class as the parent population. With this scheme, the degree of correlation between the surfaces is unimportant and the operator itself does not introduce any bias.

Mutations are introduced by changing some of the elements of the Hermitian array employed in the generation of a given surface (see figure 3). In this case, provided that the new numbers, $M_j$ and $N_j$, are zero-mean Gaussian-distributed random numbers with unit standard deviation, and the hermiticity of the array is conserved, the new surface will belong to the statistical class specified for the search space.

**3.3.2 Selection.** The selection operator generates, through a deterministic process, the set of surfaces $P_{\mu}^{(g)}$ that will serve as the population for the next iteration of the algorithm. There are two selection procedures employed in evolution strategies. The first one is known as
the 'elitist' or $(\mu + \lambda)$ strategy, whereas the second one is called the 'non-elitist' or $(\mu, \lambda)$ strategy. In the $(\mu, \lambda)$ scheme, the elements to be selected belong, exclusively, to the secondary population $P_{\lambda}^{(g)}$. An important consequence of this is the possibility that the best elements of the new population $P_{\lambda}^{(g+1)}$ are less fit than the best element of the previous population $P_{\lambda}^{(g)}$. This possible deterioration of the fitness values helps the algorithm avoid regions of attraction that could lead to premature convergence to a local minimum [28]. Of course, if the deterioration persists, the algorithm diverges.

In the elitist or $(\mu + \lambda)$ scheme, the new population $P_{\mu}^{(g+1)}$ is drawn from the two sets; that is, from the initial population $P_{\mu}^{(g)}$ and the secondary one, $P_{\lambda}^{(g)}$. Only the best elements are kept. In this case, there is never a deterioration of the fitness value and, for the case of minimization, the elitist scheme guarantees a monotonic decrease of the fitness values. This fact, however, makes the algorithm prone to premature convergence, as it can be trapped in regions of attraction associated with local minima.

# 4. Results

Ultimately, the data that serves as input to the algorithm should be obtained experimentally. However, in order to study and optimize the performance of the algorithms, in these studies we use data obtained through a rigorous numerical solution of the direct scattering problem [24].

For the two strategies explored, $(\mu + \lambda)$ and $(\mu, \lambda)$, each element of the initial population consisted of a realization of a zero-mean stationary Gaussian-correlated Gaussian random process with a $1/e$-value of the correlation function $a = 2\lambda$ and standard deviation of heights $\delta = 0.5\lambda$. Also, we chose the typical values $\mu = 10$ and $\lambda = 100$ [34]. In both cases, the maximum number of iterations was $g = 300$, which also provided the termination criterion. To illustrate the procedure, we now present a detailed discussion of the recovery of a profile.

The surface profile used to generate the original scattering data is shown in figure 4. The surface was sampled at intervals of $\lambda/10$. Since the time of computation required to find the optimum increases with the number of sampling points on the surface, and the direct problem needs to be solved many times, in order to keep the problem to a manageable size we chose a surface with $N = 128$ sampling points.

The data from which the profile is to be recovered were obtained by illuminating the surface in figure 4 from four different directions, defined by the angles of incidence $\theta_0 = -60^\circ, -30^\circ$,

<figure><img src="image_5.png" /><figcaption>Figure 4. Profile used in the generation of the scattering data.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5. Scattered intensity produced by the surface depicted in figure 4 for the case of normal incidence.</figcaption></figure>

0°, and 40°. In figure 5, we show the scattering pattern produced by the surface shown in figure 4 for the case of normal incidence.

For reference, we first present results without the use of recombination. The target profile was searched starting from 30 different, and randomly chosen initial states. As it can be expected when using a heuristic method, not in all of these attempts to recover the profile the algorithms converged to the target surface. We found, however, that in most cases a low value of $f(z_n)$ corresponded to a profile that was close to the original one (some of the exceptions are discussed below). So, the final value of $f(z_n)$ was used as the main criterion to choose the final reconstructed profile.

In figure 6 we present results obtained with the two evolution strategies studied considering the same initial population. To facilitate the visualization of the results, the original profile is shown with a dotted curve. The profile retrieved with the elitist strategy is shown with a solid curve, whereas the profile retrieved with the non-elitist strategy is drawn with a dashed curve. The vertical and horizontal displacements of the profile are understandable, as the far-field intensity is insensitive to such shifts. On the other hand, such displacements are unimportant for practical profilometric applications. Qualitatively at least, it is considered that in this case the two algorithms were able to retrieve the profile.

Since the Euclidean norm of the difference between {$z_n$} and {$\zeta_n$} is sensitive to vertical and lateral displacements of the profiles, it would not constitute an appropriate measure of the accuracy of the solution. For this, let us first define, respectively, the sample mean and

<figure><img src="image_7.png" /><figcaption>Figure 6. Reconstruction of the surface profile using the ($\mu + \lambda$) strategy (solid curve) and the ($\mu, \lambda$) strategy (dashed curve). The original profile is plotted with a dotted line.</figcaption></figure>
variance of the original profile as

$$ \bar{\zeta} = \frac{1}{N} \sum_{n-1}^{N} \zeta_n, \qquad (11) $$

$$ \sigma_{\zeta}^{2} = \frac{1}{N-1} \sum_{n-1}^{N} (\zeta_{n} - \bar{\zeta})^{2}, \qquad (12) $$

and similarly for the trial profile. Defining the new profiles

$$ \Delta z_n = z_n - \bar{z}, \qquad (13) $$

$$ \Delta \zeta_n = \zeta_n - \bar{\zeta}, \qquad (14) $$

we propose the following metric for the resemblance of the profiles

$$ \epsilon = \frac{1}{2\sigma_{\zeta}^2} \min \left\{ \frac{1}{(N - |l|)} \sum_{n} (\Delta z_{n+l} - \Delta \zeta_n)^2 \right\}. \qquad (15) $$

The subtraction of the sample means and the free subindex $l$ in equation (15) reduce the sensitivity of the error parameter $\epsilon$ to vertical and lateral displacements of the recovered profile. For perfectly correlated surfaces, this parameter is zero. On the other hand, for un-correlated random surfaces that belong to the same statistical class, the average $\epsilon$ should be of order one (forgetting about the shift). The recovered profiles shown in figure 6 have error parameters $\epsilon_{(\mu+\lambda)} = 0.006$ and $\epsilon_{(\mu,\lambda)} = 0.004$ for the elitist and non-elitist strategies, respectively.

The curves in figure 7 represent the convergence behaviour associated with the reconstructions of figure 6. The solid curve shows a monotonic decrease of the fitness value, which is typical of the elitist strategy. The dashed curve shows the convergence of the non-elitist strategy towards the stationary point. One can see that the fitness value can increase in some iterations and that the convergence of the algorithm is slower. However, after 300 iterations, the respective fitness values $f(z_n)_{(\mu+\lambda)} = 10.2$ and $f(z_n)_{(\mu,\lambda)} = 13.7$ reached with the elitist and non-elitist strategies are similar.

An interesting result that demonstrates the lack of uniqueness of the solution when intensity data are used is shown in figure 8. The dotted curve represents the original profile, while the reconstructions obtained with the strategies $(\mu + \lambda)$ and $(\mu, \lambda)$ are represented with solid and broken lines, respectively. For this numerical experiment, the initial population was the same for the two strategies, but different from the one used in the examples of figure 6. One can see that, in this case, the recovered profiles do not resemble the original one. They have

<figure><img src="image_8.png" /><figcaption>Figure 7. Convergence behaviour of the $(\mu + \lambda)$ strategy (solid curve) and the $(\mu, \lambda)$ strategy (dashed curve).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 8. Reconstruction of the surface profile using the $(\mu + \lambda)$ (solid line) and $(\mu, \lambda)$ (dashed line) strategies starting from a different initial population than in figure 6. The original profile is plotted with a dotted line.</figcaption></figure>

associated error parameters and fitness values $\epsilon_{(\mu+\lambda)} = 0.205$ and $f(z_n)_{(\mu+\lambda)} = 23.5$ (elitist) and $\epsilon_{(\mu,\lambda)} = 0.291$ and $f(z_n)_{(\mu,\lambda)} = 21.8$ (non-elitist). As explained below, these results illustrate a curious symmetry property of the scattering problem for situations in which the Kirchhoff approximation is valid.

The profiles recovered in figure 8 are presented in figure 9, but reflected with respect to the $x_1$ and $x_3$ axes; that is, we replace $z(x_1)$ by $-z(-x_1)$. Surprisingly, one observes that the resulting profiles resemble the sought one. The error parameters are now $\epsilon_{(\mu+\lambda)} = 0.027$ (elitist) and $\epsilon_{(\mu,\lambda)} = 0.007$ (non-elitist). To better understand this result, let us consider the direct scattering problem in the Kirchhoff approximation. The far-field scattering amplitude $R_{p,s}(q | k)$ can be written in the form [2]

$$R_{p,s}^{(K)}(q | k) = F_{p,s}(\theta_0, \theta_s) \int_{-\infty}^{\infty} dx_1 \exp\{-iv_x x_1 - iv_z \zeta(x_1)\}, \quad (16)$$

where $F_{p,s}(\theta_0, \theta_s)$ is an angular factor given by

$$F_{p,s}(\theta_0, \theta_s) = \pm \frac{1 + \cos \theta_0 \cos \theta_s - \sin \theta_0 \sin \theta_s}{\cos \theta_0 + \cos \theta_s}, \quad (17)$$

$v_x = k - q$, and $v_z = \alpha_0(k) + \alpha_0(q)$, with $k$ and $q$ defined as in equation (6). The ‘+’ sign in equation (17) is for p polarization, and the ‘-’ sign for s polarization.

<figure><img src="image_10.png" /><figcaption>Figure 9. Curves representing the profiles shown in figure 8 reflected with respect to both axes (circles). The original profile is shown with a dotted curve. The reconstructions obtained with the $(\mu+\lambda)$ and $(\mu, \lambda)$ strategies are represented with a solid and a dashed curve, respectively.</figcaption></figure>
If we now consider the mirror profile $\zeta_m(x_1)$, such that $\zeta_m(x_1) = -\zeta(-x_1)$, we can write its scattering amplitude as

$$[R_{p,s}^{(K)}]_m(q | k) = F_{p,s}(\theta_0, \theta_s) \int_{-\infty}^{\infty} dx_1 \exp\{-iv_x x_1 + iv_z \zeta(-x_1)\}. \quad (18)$$

With the change of variable $u = -x_1$, it can be readily shown that the two intensity patterns are equal:

$$I_{p,s}^{(K)}(q | k) = [I_{p,s}^{(K)}(q | k)]_m. \quad (19)$$

Thus, the validity of the Kirchhoff approximation leads to multiple solutions of the inverse scattering problem. It should be mentioned that, in general, the rigorous solution of the direct problem does not have this kind of symmetry. It is thus tempting to think that multiple scattering effects, which invalidate the results obtained with the Kirchhoff approximation, reduce the number of possible solutions of the inverse problem.

To complete this example, we discuss some results obtained when dominant recombination is included within the evolutionary loop. The number of elements to be combined is $\rho_d = \mu$. Furthermore, we consider the same illumination conditions, number of realizations, and initial states as those employed in the studies that led to figure 6.

With recombination, we observed a deterioration in the convergence behaviour of the algorithms. The reconstructed profiles corresponding to the lowest fitness values, $f(z_n)_{(\mu/\rho_d+\lambda)} = 30.3$ and $f(z_n)_{(\mu/\rho_d, \lambda)} = 29.4$, retrieved with the $(\mu/\rho_d+\lambda)$ and $(\mu/\rho_d, \lambda)$ strategies are shown in figure 10 with solid and broken lines, respectively. To facilitate the comparison, the original profile is drawn with a dotted curve. Although at first sight, there seem to be important differences between the original and recovered profiles, lateral and vertical translations show that the main features are preserved in the reconstructions. This is verified by the calculated error parameters, which are $\epsilon_{(\mu/\rho_d+\lambda)} = 0.042$ and $\epsilon_{(\mu/\rho_d, \lambda)} = 0.043$.

We close this section with a summary of more extensive studies of the retrieval of surface profile functions. First, in figure 11, we present four profiles belonging to the same statistical class as the surface of figure 4 ($\delta = 0.5\lambda$ and $a = 2\lambda$), together with the reconstructed ones using the elitist and non-elitist strategies without recombination. If, in all cases, we keep only the reconstruction with the lowest value of the fitness functional, the error parameters associated with the reconstructions are in the range $0.002 \le \epsilon \le 0.014$. So, at least for this statistical class and the chosen geometry, it can be said that the profile can be retrieved with a high degree of confidence.

One should be careful, however, in trying to extrapolate on these results. As an example of this, we present reconstructions of a surface belonging to different statistical class ($\delta = \lambda$ and

<figure><img src="image_11.png" /><figcaption>Figure 10. Reconstruction of the surface profile using the $(\mu/\rho_d + \lambda)$ strategy (solid curve) and, the $(\mu/\rho_d, \lambda)$ strategy (dashed curve). The original profile is plotted with a dotted line.</figcaption></figure>
<figure><img src="image_12.png" /></figure>

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

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

<figure><img src="image_15.png" /><figcaption>Figure 11. Reconstruction of four different surface profiles using the $(\mu + \lambda)$ (solid curve) and $(\mu, \lambda)$ (dashed curve) strategies without recombination. The original profile is plotted with a dotted line. (a) $f(z_n)_{(\mu+\lambda)} = 17.1$, $\epsilon_{(\mu+\lambda)} = 0.003$, and $f(z_n)_{(\mu,\lambda)} = 27.1$, $\epsilon_{(\mu,\lambda)} = 0.008$; (b) $f(z_n)_{(\mu+\lambda)} = 29.3$, $\epsilon_{(\mu+\lambda)} = 0.015$, and $f(z_n)_{(\mu,\lambda)} = 25.4$, $\epsilon_{(\mu,\lambda)} = 0.003$; (c) $f(z_n)_{(\mu+\lambda)} = 10.7$, $\epsilon_{(\mu+\lambda)} = 0.005$, and $f(z_n)_{(\mu,\lambda)} = 5.7$, $\epsilon_{(\mu,\lambda)} = 0.002$; (d) $f(z_n)_{(\mu+\lambda)} = 14.0$, $\epsilon_{(\mu+\lambda)} = 0.014$, and $f(z_n)_{(\mu,\lambda)} = 22.2$, $\epsilon_{(\mu,\lambda)} = 0.048$.</figcaption></figure>

$a = 2\lambda$). The surface is in fact the same one shown in figure 4, but expanded in the vertical direction by a factor of two. The surface and the attempted reconstructions are shown in figure 12. The reconstructed profiles have fitness values, $f(z_n)_{(\mu+\lambda)} = 33.6$ and $f(z_n)_{(\mu,\lambda)} = 26.0$, and error parameters $\epsilon_{(\mu+\lambda)} = 0.051$ and $\epsilon_{(\mu,\lambda)} = 0.339$ for the elitist and non-elitist strategies, respectively.

<figure><img src="image_16.png" /><figcaption>Figure 12. Reconstruction of a realization of a surface profile with $\delta = \lambda$ and $a = 2\lambda$ using the $(\mu + \lambda)$ (solid curve) and $(\mu, \lambda)$ (dashed curve) strategies without recombination. The original profile is plotted with a dotted line.</figcaption></figure>
From the error parameters and the appearance of the curves, it is clear that the reconstruction obtained with the elitist strategy is a much better approximation to the original profile. However, the fitness value obtained with the non-elitist strategy is lower than the one obtained with the elitist strategy. This is probably an indication of the lack of uniqueness of the solution.

## 5. Conclusions

We have presented a study of two evolutionary algorithms to solve inverse scattering problems. Starting from far-field intensity data and using both the ($\mu + \lambda$) and ($\mu, \lambda$) strategies, we have successfully retrieved the surface profile that generated the scattering data. The time of computation is similar, but the elitist strategy ($\mu + \lambda$) uses twice as much memory as the non-elitist strategy ($\mu, \lambda$).

Although it appears that including the recombination operation degrades the performance of the inversion procedures, the partial results obtained are far from being conclusive. Moreover, apart from the works of Beyer [38] on the ($\mu/\rho_d, \lambda$) strategy, we are not aware of other systematic studies on the subject, and in particular with the ($\mu/\rho_d + \lambda$) strategy.

The numerical evidence presented here suggests that the evolutionary strategies are well adapted for inverse scattering applications. Further work is required to understand the influence of the dominant recombination on the dynamics of the inversion scheme. Another important consideration is the effect that noisy data have on the reconstructions. Some results (not presented here for brevity) have been obtained with input data contaminated with noise; they indicate that the inversion schemes are fairly insensitive to the presence of Gaussian additive noise [39]. More systematic studies would be necessary to assess the tolerances of the inversion method to the presence of more general kinds of noise.

Since the solution of the inverse problem is not necessarily unique, it seems appropriate to use a stochastic method to search for the solution. The fitness function has many local minima, and the initialization of the algorithm plays an important role in its convergence. In most of the cases we have studied, the stationary point at the end of the evolutionary loop not only seems to be associated with a low value of the fitness functional, but it also provides the best approximation to the problem.

We have also presented cases in which the problem has more than one solution and cases in which the best solution does not correspond to the surface with the lowest fitness value. In such situations, one would have to perform further tests to decide on the best one. Although the problem has many facets and can be rather complex, the results obtained so far are encouraging.

The success of an inversion scheme based on intensity information opens the possibility of implementing such a procedure experimentally. However, further work is needed, not just regarding the physical aspects of the problem, which have been simplified by our assumptions, but also regarding aspects related to the performance of the evolutionary algorithms.

## Acknowledgments

DM gratefully acknowledges financial support of the Conseil Régional de Champagne-Ardenne. The work of ERM was supported by CONACYT, under grant 47712-F.

## References

[1] Beckmann, P. and Spizzichino, A., 1963, *The Scattering of Electromagnetic Waves from Rough Surfaces* (London: Pergamon Press).

[2] Ogilvy, J. A., 1991, *Theory of Wave Scattering from Random Rough Surfaces* (Bristol: Institute of Physics Publishing).
[3] DeSanto, J. A. and Brown, G. S., 1986, Analytical techniques for multiple scattering from rough surfaces. In: E. Wolf (Ed.) *Progress in Optics*, vol. XXIII, (Amsterdam: North-Holland), pp. 1-62.

[4] Shchegrov, A. V., Maradudin, A. A., and Méndez, E. R., 2004, Multiple scattering of light from randomly rough surfaces. in E. Wolf (Ed.) *Progress in Optics*, vol. 46 (Amsterdam: Elsevier), pp. 117-233.

[5] Léger, D., Mathieu, E., and Perrin, J. C., 1975, Optical surface roughness determination using speckle correlation technique. *Applied Optics*, **14**, 872-877.

[6] Chandley, P. J., 1976, Determination of the autocorrelation function of height on a rough surface from coherent light scattering. *Optical and Quantam Electronics*, **8**, 329-333.

[7] Welford, W. T., 1977, Optical estimation of statistics of surface roughness from light scattering measurements. *Optical and Quantam Electronics*, **9**, 269-287.

[8] Elson, J. M. and Bennett, J. M., 1979, Relation between the angular dependence of scattering and the statistical properties of optical surfaces. *J. Opt. Soc. Am.*, **69**, 31-47.

[9] Stover, J. C., Serati, S. A., and Gillespie, C. H., 1984, Calculation of surface statistics from light scatter. *Opt. Eng.*, **23**, 406-412.

[10] Malyshkin, V., Simeonov, S., McGurn, A. R., and Maradudin, A. A., 1977, Determination of surface profile statistics from electromagnetic scattering data. *Optics Letters*, **22**, 58-60.

[11] Leskova, T. A., Maradudin, A. A., Novikov, I. V., Shchegrov, A. V., and Méndez, E. R., 1998, Design of one-dimensional band-limited uniform diffusers of light. *Applied Physics Letters*, **73**, 1943-1945.

[12] Méndez, E. R., García-Guerrero, E. E., Escamilla, H. M., Maradudin, A. A., Leskova, T. A., and Shchegrov, A. V., 2001, Photofabrication of random achromatic optical diffusers for uniform illumination. *Applied Optics*, **40**, 1098-1108.

[13] Maradudin, A. A., Simonsen, I., Leskova, T. A., and Méndez, E. R., 2001, Design of one-dimensional Lambertian diffusers of light. *Waves Random Media*, **11**, 529-533.

[14] Méndez, E. R., García-Guerrero, E. E., Leskova, T. A., Maradudin, A. A., Muñoz-López, J., and Simonsen, I., 2002, Design of one-dimensional random surfaces with specified scattering properties. *Applied Physics Letters*, **81**, 798-800.

[15] Maradudin, A. A., Leskova, T. A., and Méndez, E. R., 2003, Two-dimensional random surfaces that act as circular diffusers. *Optics Letters*, **28**, 72-74.

[16] Méndez, E. R., Leskova, T. A., Maradudin, A. A., and Muñoz-López, J., 2004, Design of two-dimensional random surfaces with specified scattering properties. *Optics Letters*, **29**, 2917-2919.

[17] Wombel, R. J. and DeSanto, J. A., 1991, Reconstruction of rough-surface profiles with the Kirchhoff approximation. *Journal of the Optical Society of America A*, **8**, 1892-1897.

[18] Wombel, R. J. and DeSanto, J. A., 1991, The reconstruction of shallow rough-surfaces profiles from scattered field data. *Inverse Problems*, **7**, L7-L12.

[19] Quartel, J. C. and Sheppard, C. J. R., 1996, Surface reconstruction using an algorithm based on confocal imaging. *Journal of Modern Optics*, **43**, 469-486.

[20] Quartel, J. C. and Sheppard, C. J. R., 1996, A surface reconstruction algorithm based on confocal interferometric profiling. *Journal of Modern Optics*, **43**, 591-605.

[21] Macías, D., Méndez, E. R., and Ruiz-Cortés, V., 2002, Inverse scattering with a wave-front-matching algorithm. *Journal of the Optical Society of America A*, **19**, 2064-2073.

[22] Macías, D., Olague, G., and Méndez, E. R., 2002, Surface profile reconstruction from scattered intensity data using evolutionary strategies. In: S. Cagnoni et al. (Eds) *EvoWorkshops*, LNCS 2279 (Berlin: Springer), pp. 233-244.

[23] Macías, D., Olague, G., and Méndez, E. R., 2003, Hybrid evolution strategy-downhill simplex algorithm for inverse light scattering problems. In: S. Cagnoni et al. (Eds) *EvoWorkshops*, LNCS 2611 (Berlin: Springer), pp. 339-409.

[24] Maradudin, A. A., Michel, T., McGurn, A. R., and Méndez, E. R., 1990, Enhanced backscattering of light from a random grating. *Annals of Physics (N.Y.)*, **203**, 255-307.

[25] Kold, T. G., Lewis, R. M., and Torczon V., 2003, Optimization by direct search: new perspectives on some classical and modern methods. *SIAM Review*, **45**, 385-482.

[26] Nelder, J. and Mead, R., 1965, A simplex method for function optimization. *Computer Journal*, **7**, 304-308.

[27] Kirkpatrick, S., Gelatt, C. D., and Vecchi, M. P., 1983, Optimization by simulated annealing. *Science*, **220**, 671-680.

[28] Bäck, Th., Hammel, U., and Schwefel, H.P., 1997, Evolutionary computation: comments on the history and current state. *IEEE Transactions on Evolutionary Computation*, **1**, 3-17.

[29] Weile, D. S. and Michielssen, E., 1997, Genetic algorithm optimization applied to electromagnetics: a review. *IEEE Transactions on Antennas and Propagation*, **45**, 343-353.

[30] Hageman, J. A., Wehrens, R., van Sprang, H. A., and Buydens, L. M. C., 2003, A hybrid genetic algorithm - tabu search approach for optimising multilayer optical coatings. *Analytica Chimica Acta*, **490**, 211-222.

[31] Hodgson, R. J., 2000, Genetic algorithm approach to particle identification by light scattering. *Journal of Colloid and Interface Science*, **229**, 399-406.

[32] Djurisic, A. B., Elazar, J. M., and and Rakic, A. D., 1997, Modeling the optical constants of solids using genetic algorithms with parameter space size adjustment. *Optics Communications*, **134**, 407-414.

[33] Holland, J. H., 1992, *Adaption in Natural and Artificial Systems* (Cambridge, Massachusetts: MIT Press/Bradford Books).

[34] Schwefel, H. P., 1995, *Evolution and Optimum Seeking* (New York: John Wiley).
[35] Koza, J. R., 1998, Genetic programming. In: J.G. Williams and A. Kent (Eds) *Encyclopedia of Computer Science and Technology* (New York: Marcel-Dekker), 29–43.

[36] Fogel, L. J., 1962, Autonomous automata. *Ind. Research*, **4**, 14–19.

[37] Salomon, R., 1997, Evolutionary algorithms and gradient search: similarities and differences. *IEEE Transactions on Evolutionary Computation*, **2**, 45–55.

[38] Beyer, H. G., 2001, *The Theory of Evolution Strategies* (Berlin: Springer-Verlag).

[39] Macías, D., 2003, Estudios numéricos de esparcimiento inverso de ondas electromagnéticas. PhD thesis, Centro de Investigación Científica y de Educación Superior de Ensenada, México.