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

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.


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][2][3][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][6][7][8][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][12][13][14]. Extensions to the case of general two-dimensional angular 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][18][19][20][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.

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 = ζ (x 1 ). The region ζ (x 1 ) > x 3 is the vacuum, the region ζ (x 1 ) < x 3 is a perfect conductor, and the plane of incidence is the x 1 x 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 where α 0 (k) = (ω/c) 2 − k 2 , ω is the frequency of the field, and c is the speed of light in vacuum. A time dependence of the form exp(−iωt) is assumed, but explicit reference to it is suppressed. The field scattered by the surface is where the scattering amplitude R p,s (q | k) can be written in the form [24] The factor χ (k, ζ (x 1 )) is and F(x 1 | ω) represents a source function, given by where ∂/∂ N = N·∇ = [−ζ (x 1 )(∂/∂ x 1 )+(∂/∂ x 3 )] is the non-unit normal derivative operator. The angles of incidence θ 0 and scattering θ s are related to the components of the wavenumbers k and q that are parallel to the mean surface through the expressions The far-field intensity I p,s (q | k) is defined as In the present work, the goal is to retrieve the unknown surface profile function ζ (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.

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 boundsconstrained approximation problem. It is assumed that we have access to far-field angleresolved 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.

Definition of the fitness functional
The inverse scattering problem can be reformulated in terms of the fitness (objective) functional where the symbol · 2 represents the Euclidean norm of the intensity vectors as functions of the scattering wavevector components q, N 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 (ζ c (x 1 )).

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: Here, the angle brackets denote an average over the ensemble of realizations of the surface profile function, δ = ζ 2 (x 1 ) 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 Here, N represents the total number of points on the surface, L represents its length, χ n = −L/2 + (n − 0.5) x are the sampling points spaced by x along x 1 , q j = −π/ x + 2π ( j − 0.5)/L are the sampling points in Fourier space, and ζ n = ζ (χ 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 {ζ n }, it is required that the complex array {M j + i N 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).

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][27][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][31][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 g µ | g = 0 of µ 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 g λ of λ elements is generated through the application of the 'genetic' operations of recombination and mutation over the elements of the initial population P g µ . This represents the start of the main evolutionary loop.
At this point, it is pertinent to mention that we have previously used λ 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  the optimum, mutation could be the only operator employed. Once the secondary population P 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 g λ leading to promising regions of the search space will be retained, through some selection scheme, as part of the population P 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 g µ and P g λ , remain constant throughout the entire search process.

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

Selection.
The selection operator generates, through a deterministic process, the set of surfaces P 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 (µ + λ) strategy, whereas the second one is called the 'non-elitist' or (µ, λ) strategy. In the (µ, λ) scheme, the elements to be selected belong, exclusively, to the secondary population P g λ . An important consequence of this is the possibility that the best elements of the new population P g+1 λ are less fit than the best element of the previous population P 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 (µ + λ) scheme, the new population P g+1 µ is drawn from the two sets; that is, from the initial population P g µ and the secondary one, P 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 guaranties 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.

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, (µ + λ) and (µ, λ), 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λ and standard deviation of heights δ = 0.5λ. Also, we chose the typical values µ = 10 and λ = 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 λ/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 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 {ζ 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 and similarly for the trial profile. Defining the new profiles we propose the following metric for the resemblance of the profiles The subtraction of the sample means and the free subindex l in equation (15) reduce the sensitivity of the error parameter to vertical and lateral displacements of the recovered profile. For perfectly correlated surfaces, this parameter is zero. On the other hand, for uncorrelated random surfaces that belong to the same statistical class, the average should be of order one (forgetting about the shift). The recovered profiles shown in figure 6 have error parameters (µ+λ) = 0.006 and (µ,λ) = 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 ) (µ+λ) = 10.2 and f (z n ) (µ,λ) = 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 (µ + λ) and (µ, λ) 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  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 (µ+λ) = 0.027 (elitist) and (µ,λ) = 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] where F p,s (θ 0 , θ s ) is an angular factor given by F p,s (θ 0 , θ s ) = ± 1 + cos θ 0 cos θ s − sin θ 0 sin θ s cos θ 0 + cos θ s , v x = k − q, and v z = α 0 (k) + α 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. If we now consider the mirror profile ζ m (x 1 ), such that ζ m (x 1 ) = −ζ (−x 1 ), we can write its scattering amplitude as With the change of variable u = −x 1 , it can be readily shown that the two intensity patterns are equal: 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 ρ d = µ. 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 ) (µ/ρ d +λ) = 30.3 and f (z n ) (µ/ρ d ,λ) = 29.4, retrieved with the (µ/ρ d +λ) and (µ/ρ d , λ) 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 (µ/ρ d +λ) = 0.042 and (µ/ρ d ,λ) = 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 (δ = 0.5λ and a = 2λ), 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 ≤ ≤ 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 (δ = λ and 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.

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 (µ + λ) and (µ, λ) strategies, we have successfully retrieved the surface profile that generated the scattering data. The time of computation is similar, but the elitist strategy (µ + λ) uses twice as much memory as the non-elitist strategy (µ, λ).
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 (µ/ρ d , λ) strategy, we are not aware of other systematic studies on the subject, and in particular with the (µ/ρ d + λ) 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.