1. Introduction

The inverse scattering problem is the determination of physical parameters and features of an unknown object from a limited set of measurements of scattered field data. Ideally, the solution leads to a quantitative description of the scatterer such as its spatial distribution of permittivity, conductivity or other internal constitutional parameters. There are many important imaging and remote sensing applications which require practical inversion algorithms that are applicable when the object in question is strongly scattering and this is well known to be a very difficult problem to solve. While advances have been made primarily using iterative techniques [1, 2], long computational times are necessary and convergence is not guaranteed, especially if only incomplete and noisy data are measured.

Most computationally feasible algorithms assume a weakly scattering condition in order to linearize the inversion procedure [3]. We have developed a method [4, 5] which is simple to implement and which can be applied to any inverse scattering problem including those with targets in cluttered backgrounds and which can provide a good estimate of the image of the scatterer. The success of the method inevitably depends on the quality and extent of the data and it is advantageous if one has data for several illuminating wavelengths. The quality of the reconstruction depends on a number of factors but at the very least can be useful to provide a quick look at the target features and serve as an initial guess for a more rigorous inversion scheme.

 

Fig. 1. Scattering geometry.

Download Full Size | PDF

 

Fig. 2. Plane wave scattering in the Fourier space (k-space) of the object. The direction of the incident field r̂₀ and the direction r̂_S of a particular plane wave component of the scattered field define a point at the Ewald circle. Changing the wavelength of the incident radiation as well as the directions of incident and measured scattered field components provides sufficient information about k-space to compute the image of the object via an inverse Fourier transform.

Download Full Size | PDF

Here, we describe the use of this inverse scattering approach for design. This means, our goal is to find the permittivity distribution of a strongly scattering structure which exhibits a specified scattering behavior. We also want to impose additional constraints, namely the approximate shape or support of the structure to be designed. This inverse scattering problem defines a rather generic design problem, and its solution relates to a number of important applications. In particular, the design of diffractive optical elements in the subwavelength and resonant regime can be addressed. Other applications include the design of artificial materials with low reflectivity for a finite range of wavelengths, which could be used for stealth technology.

In the case of imaging from scattered field data, it is clearly important that the content of the reconstructed image can be relied upon in order to discern features that might have some practical significance. The structure synthesis problem is not so constrained. In principle, any structure which provides the desired scattering characteristics is a solution to the synthesis problem. The only concern might be the ease with which the predicted scattering structure can be fabricated, i.e. in terms of its material properties and structural features.

The inversion and synthesis method we use is based on diffraction tomography for which inversion techniques are formulated as Fourier inversion procedures [3]. Complex scattered field data are collected at a number of scattering angles in the far field for each angle of illumination as shown in Fig. 1. The scattered field data under these conditions, based on the first-order Born approximation are mapped onto the Ewald sphere in k-space [5] and with sufficient k-space coverage can be inverse Fourier transformed to provide information about the scattering object (Fig. 2).

Diffraction tomography algorithms can be executed with ease, but the physical premise upon which they are based is rarely acceptable in practice [6]. Images obtained from strongly scattering objects when diffraction tomography algorithms are applied blindly, lead to grossly distorted and often meaningless images. However, that reconstruction does incorporate information about the true scattering structure, which merely needs to be extracted. Here a nonlinear filtering method is used, which is applied to the result of a diffraction tomography algorithm in order to recover meaningful images of strongly scattering targets. Some preprocessing of the data is necessary for this to succeed which is described in section 3. The advantages of this approach are that it is computationally fast requiring only two fast Fourier transforms. Being Fourier-based, it permits the easy incorporation of well-known methods for dealing with limited noisy sampled data [7, 8], which, however, is not really an issue in the synthesis problem. The key step in the method is the preprocessing of the scattering data to generate a so-called minimum phase function that permits the nonlinear filtering operation to be executed in a stable fashion to recover the image of the target object. In this paper we are able to assume that scattered field data are specified all around the object which further simplifies the implementation of this approach.

We note that our investigation was carried out for two-dimensional scattering geometries assuming a scalar diffraction model. However, in as much as the concept of diffraction tomography can be translated from two to three dimensions, there is no apparent reason why our approach cannot be applied to three-dimensional geometries as well.

2. Nonlinear inverse scattering algorithm

Consider a penetrable scattering object, V(r) in a homogeneous background of permittivity ε ₀. The object is related to the permittivity by V(r) = k ²[ε(r)-1]. The scattered field Ψ_s(r,k r̂₀), due to the interaction of incident wave Ψ₀(r ,k r̂₀) with the target V(r), is given in two dimensions by the integral equation [7]

where r̂₀ is the unit vector denoting the direction of the incident field. The solution to this equation requires the knowledge of total field Ψ(r,k r̂₀) within the object volume D, which is not possible when V(r) is unknown. Using the first Born approximation, one assumes that the total field Ψ(r,k r̂₀) in the integral can be replaced by the known incident field Ψ₀(r ,k r̂₀), which linearizes this inversion problem and permits one to find a solution in principle [3, 7]. The integral in Eq. (1) reduces to a Fourier transform relation between the scattered field data Ψ_s(r,k r̂₀) and V(r). Consequently, each measurement sample of the scattered far field can be related to one sample of the Fourier transform of V(r) by applying the data mapping illustrated in Fig. 2. Physically this requires that the scattering from the object is extremely weak in order for the total field everywhere within the object to be well approximated by the incident field.

When the first born approximation is not valid, the same data inversion step yields information about the secondary source or contrast source function V_B(r,k r̂₀)i.e.

This first born reconstruction is modulated by the field pattern within D, which will be different for every illumination direction r ₀. In diffraction tomography, the scattered field data for all incident field directions is combined in k-space and a Fourier inversion of that data one expects to provide an estimate for V(r). When the first Born approximation is valid, then Ψ(r,k r̂₀)~Ψ₀(r,k r̂₀) and V_B(r)~V(r) at least to within low-pass spatial filtering effects resulting from the available k-space coverage. When the first Born approximation is not valid then V_B(r)~V(r)⟨Ψ(r)⟩ where ⟨Ψ(r)⟩ is a complex and noise-like term with a characteristic range of spatial frequencies dominated by the bandwidth of the source.

The image of the product V〈Ψ〉 can be separated using homomorphic filtering. This product exhibits spatial fluctuations characteristic of the wavelengths being employed as well as the spatial fluctuations of the permittivity. For incremental illumination wavelength changes, the 〈Ψ〉 term will change quite considerably but V(r) need not. The first step in homomorphic filtering is to take the logarithm of the product of V〈Ψ〉 and then employ standard spatial filtering in the associated Fourier domain to remove the field component. However, since

there are inevitably numerous numerical problems associated with forming log(V〈Ψ〉) and taking its Fourier transform. When the magnitude of the term V〈Ψ〉 is close to zero, its logarithm becomes singular, and when the phase of V〈Ψ〉 has a range that exceeds 2π, the resulting wrapped phase introduces spurious spatial frequencies in the cepstrum. Phase unwrapping is exceedingly difficult, especially in two and higher dimensional problems, because of zeros in the field, which are associated with wavefront dislocations [8]. The Fourier transform of log (V〈Ψ〉) generates the cepstrum of the function [9] and spatial filtering in the cepstral domain is known as homomorphic filtering. Phase discontinuities correspond to abrupt and high spatial frequency features, which make successful homomorphic filtering impossible. Modifications to homomorphic filtering, such as differential cepstral filtering [10] have been developed. With our approach, however, we restore the effectiveness of homomorphic filtering by preprocessing the data, and therefore force them to satisfy the minimum phase condition.

3. Minimum-phase based homomorphic filtering

The concept of minimum phase is well understood in one dimensional problems [8]. A one dimensional signal f(x) is minimum phase if and only if its Fourier transform,

has a zero-free upper half plane, i.e. it has no zeros for v > 0 where z = u + iv. The property of minimum phase has a number of interesting consequences. When f(x) is causal, the phase of F(u) is continuous and lies between -π and +π, i.e. the phase is always unwrapped (non-discontinuous) and f(x) has most of its energy close to the origin; phase recovery from |F(u)| also becomes possible using a logarithmic Hilbert transform [8].

The meaning of a zero free half plane in two dimensions is less easily conveyed. Dudgeon and Mersereau [11] describe a two dimensional minimum phase function which is absolutely summable, and whose inverse and complex cepstrum are absolutely summable and have the same (convex) region of support. Taking the logarithm of a bandlimited function produces a bandlimited function and hence a cepstrum of finite support only for minimum phase functions. An important feature of a minimum phase function is that the phase is a continuous function bounded between -π and π and “minimum” in this sense indicates that the phase is already unwrapped. It is possible to enforce the minimum phase condition on a function by applying Rouche’s theorem [12], a multidimensional version of which can be found in ref [13]. The theorem shows that if a bandlimited function ‘H’ has ‘N’ zeros within some contour, and another bandlimited function ‘F’ has ‘M’ zeros in the same contour, and if |H| > |F| on that contour, then H+F will have ‘N’ zeros in that contour i.e. the sum of the two functions will have the number of zeros equal to the number of zeros of the larger magnitude function. Consequently, adding a sufficiently large background or reference wave to a bandlimited function, F, guarantees that the sum is a minimum phase function.

One can therefore preprocess V〈Ψ〉 by introducing a reference point in the data or k- space domain to make it minimum phase prior to taking its logarithm. The first step is to make the data in k-space causal by moving it to one quadrant of k-space and then place a reference point at the origin of k-space. This is equivalent to adding a reference wave to V〈Ψ〉 which needs an amplitude just large enough to ensure that phase of V〈Ψ〉 is continuous and lies within the bounds of -π and + π. This is readily determined by inspection of the phase of the reconstructed image of V〈Ψ〉, i.e. the image resulting from the direct Fourier inversion of the data in k-space. The implementation of the homomorphic filtering algorithm then requires that a low pass filter be applied in the cepstral domain until the wavelike features associated with 〈Ψ〉 are removed from the resulting image. This spatial filtering is successful to the extent to which the field internal to the scattering structure has spatial frequencies that are distinct from those of log(V). Using different source frequencies can resolve this problem.

4. Reconstructions and structure synthesis

In order to ensure that the structure whose scattering properties we wish to define is physically meaningful, we illustrate our method by modifying the scattered field generated by a target provided by AFRL, from what has become known as the Ipswich data set [14]. The data were taken at 10GHz with 36 incident field directions and scattered field measurements every 10 degrees.

 

Fig. 3. The reconstruction of V〈Ψ〉 (left) and V using cepstral filtering.

Download Full Size | PDF

The image in Fig. 3 is of IPS010, a dielectric wedge. On the left the estimate of V〈Ψ〉 is shown computed with the PDFT algorithm [7], which in essence is the Fourier inverse of the data when prior information about the support constraint in incorporated. On the right the cepstral filtered reconstruction for V is shown. This is obtained by adding the reference as described in the previous section. Then the logarithm of the resulting Fourier inverse is computed and Fourier transformed into the cepstral domain of the signal. After applying a linear filter (low-pass) the Fourier inverse of the resulting cepstrum is computed and the final image is obtained after exponentiation. For a detailed description of this procedure we refer to references [4, 5].

For synthesizing structures with a prescribed scattering behavior the k-space measurements are modified before applying the inverse scattering procedure. As an illustrating example we modified the k-space data as shown on the left side of Fig. 4. In two numerical experiments an annulus of k-space data was forced to zero for different directions of the incident field. This corresponds to the requirement that there be no forward or backscattered radiation from this target for the specific incidence field direction. This also means that the scattering amplitudes for other incident field directions exhibit heavily modulated, low scattering amplitudes if the corresponding locations in k-space intersect the zero-annulus. The radius of the annulus was chosen to be relatively large, with the intended primary effect of suppressing forward scattering from the synthesized structures. After modifying k-space cepstral filtering was performed in each case to obtain an estimate of the permittivity distribution which would be similar in shape to the IPS010 object, and which exhibits the desired scattering behavior.

 

Fig. 4. The top figure shows the modification to the estimated structure for IPS010 which, when illuminated in the -y direction (see above) would have very little if any scattering in the forward direction. The lower figure shows a different structure which would exhibit little forward scatter when illuminated in the +x direction.

Download Full Size | PDF

In order to verify that these modified structures would have the anticipated scattering properties, we computed the scattered fields in the forward direction using a finite difference time domain algorithm. A clean and precise nulling of the forward scattered fields is not expected because of the highly nonphysical nature of the constraint imposed on the scattering distribution, i.e. an abrupt transition to zero scattering for certain wavelengths and certain directions, combined with the fact that we can only claim that the reconstruction of IPS010, in Fig. 3, is an estimate. It is a reasonable estimate however, and the use of measured scattered field data as starting point for synthesizing structures allows us to avoid potential errors and inverse crimes associated with numerically determining the scattered field from a strongly scattering structure. Our approach also allows us the possibility of realizing a physical structure that possesses the desired scattering properties by modifying the actual dielectric wedge, IPS010. The results of the finite difference time domain calculations of the scattered fields are illustrated in Figs. 5 and 6.

 

Fig. 5. Scattering amplitude in the forward direction over a range of angles for illumination in the +x direction. The green curve illustrates the scattered field from (our reconstruction of) IPS010, the red and blue curves from the two synthesized structures.

Download Full Size | PDF

 

Fig. 6. Forward scattering amplitudes for both synthesized objects when illuminated from the -y direction. Here both scattered fields show some reduction in amplitude in the forward direction compared to IPS010. In addition, the fields exhibit strong sidelobes.

Download Full Size | PDF

In Fig. 5 a drop in the field amplitude can be observed for the expected structure (red curve), but the scattering amplitude it is neither symmetric with respect to the x- axis, nor does the minimum occur for the +x direction. The other synthesized object (blue curve) has a scattering amplitude that differs little from that of IPS010.

In Fig. 6 this analysis is repeated with the incident wave pointing in -y direction. We observe that the scattered fields are not quite as expected. The synthesized object designed to have the lowest forward scattering (red curve) shows a somewhat larger amplitude lobe structure than the other object (blue curve). This result obviously requires additional interpretation. Apart from a more distinct suppression of forward scattering, the rather simple low pass filter in the cepstral domain is rather expected to generate scattering structures that exhibit a smoother, i.e. more bandlimited scattering patterns, than that specified by the design.

 

Fig. 7. The k-space data removal imposed for each of the two new synthesized objects is illustrated. The axes indicate the k-space origin and the central image shows the distribution of energy in the original IPS010 measured data.

Download Full Size | PDF

In order to try to understand these results more fully, we evaluated the exact location of the imposed annulus of zeros with respect to the k-space origin and the local pixelation or sampling of the annulus with respect to the k-variable. From Fig. 7, we can see that these zero regions have removed the zero-frequency or dc-level of the contrast source function. This results in the energy at high spatial frequencies for the two synthesized objects being more apparent in these images. From a crude inspection of the k-space energy patterns for the synthesized objects, some residual sinusoidal or cosinusoidal modulation of the forward scattering amplitude is expected. Moreover these annuli were not imposed in a symmetric fashion due to sampling grid limitations, which can explain some asymmetry in the resulting scattering amplitudes. The synthesized object exhibiting reduced scattering in the +y direction (left in Fig. 7) appears to generate a more symmetric scattered field and this is borne out from Fig. 6 (red curve) to some extent. The synthesized object corresponding to the k-space plots on the left in Fig. 7, do show more symmetry than for the case on the right, which is consistent with expectations. While these scattering characteristics might not be ideal, they illustrate the correct and anticipated trend in behavior and we would expect the associated structures to, if nothing else, provide a good initial estimate for an iterative and more rigorous design of a structure that more precisely conforms to the desired scattering patterns.

5. Conclusions

In this paper we presented reconstructions of strongly scattering structures, designed to have reduced scattering in the forward direction for certain wavelengths and certain incident field directions. The approach to solve this inverse synthesis problem was to exploit an inverse scattering method based on homomorphic filtering. This form of filtering allows one to recover an estimate of V(r) from the image of V〈Ψ〉 obtained using conventional diffraction tomography. The quality of the reconstruction depends on a number of factors, including the quantity and quality of the data made available. It is obvious that simple low-pass filtering of the spectrum of V〈Ψ〉 as opposed to the cepstrum, will do nothing. This is illustrated in the sequence of images shown in Fig. 8.

 

Fig. 8. Series of images of V〈Ψ〉 following application of a low-pass filter in k-space with increasing radius (compare with Fig. 3); we emphasize that only filtering in the cepstral domain, after enforcing the minimum phase condition, provides an estimates of V(r).

Download Full Size | PDF

In principle, when solving the scattering-structure synthesis problem, one can specify the entire volume of k-space as precisely as one wishes, provided that choice corresponds to a physically meaningful set of scattered field data. In the examples presented here, rather crude annuli of zeros were imposed in k-space, corresponding to zero scattering over a range of scattering angles and wavelengths. Upon inversion of the resulting k-space data set, object structures were computed which can be seen to be subtle modifications of the original object from which the data were originally acquired. The synthesized structures were validated using a finite difference time domain method for computing the corresponding far field patterns. A number of interesting observations were made. Qualitatively the synthesized structures behaved as designed, but a combination of factors, most notably the simple low-pass filtering employed in the cepstral domain and the relatively crude imposition of zero values in k-space, lead to smoothed structures, which in turn resulted to undulating and rather smooth scattered fields.

Future work will investigate the trade-offs here between the k-space grid density and less abrupt k-space modulations, in order to reconcile the balance between the physical allowability of the desired field characteristics (i.e., does the modified k-space represented a physical field) and the need to synthesize a structure which can be made from known materials (i.e., physically realizable permittivities) and with feature scales that permit fabrication.

MAF and MT gratefully acknowledge the support of DARPA grant W911NF-04-1-0319 in support of this research. We are also very grateful to Dr. A. Kanaev for providing the FDTD graphs.