1. Introduction

The mathematics of inverse scattering problems is important to optical imaging and microscopy, medical imaging, radar, acoustics, geophysics, and many other disciplines. Inverse scattering is the inference of properties of an inhomogeneous medium from detection of waves that are scattered by the medium. In general, scattering from inhomogeneous media [1] is a complex process where the scattered fields are nonlinear in the medium property. Because of the difficulty of modeling arbitrary inhomogeneous scatterers, approximations are frequently made that simplify the mathematics and modeling of scatterers in special cases. For various practical inverse scattering problems, specific tools have been developed that enable practical solution of specific problems. Over time, these tools have been synthesized into general methodologies that have enabled investigators in many disciplines to mutually improve the state of the art in inverse scattering. This work presents a new series that synthesizes two well-known scattering series into one. Here it is shown how two approximations typically used for linearized inverse scattering, the Born and Rytov approximations, are two extremes of a more generalized family of hybrid approximations, some of which can achieve better accuracy for particular inverse scattering problems than either the Born or Rytov approximations alone. In addition, it is shown how the hybrid approximations can be applied to more complicated multiple scattering problems.

The applicability and accuracy of the Born and Rytov scattering series has been well studied [2, 3, 4, 5, 6, 7, 8, 9]. The Born approximation is frequently employed in geometries where the radiation is collected in the backscattering direction, such as monostatic radar, B-mode ultrasound, and optical coherence tomography. In the first Born approxmation, the scattered field is modeled as a linear function of the scattering contrast (e.g. refractive index or acoustic impedance). Therefore, the reconstructed contrast is usually estimated as linear in the detected field. The Born approximation is accurate when the product of the index contrast and object size is less than one-quarter wavelength. On the other hand, the Rytov approximation is typically used when the field is detected in the forward scattering direction, as is employed in diffraction tomography, many geophysical problems, and X-ray computed tomography [10]. In the first Rytov approximation, the scattered field is modeled as an exponentially dependent function of the scattering contrast. The Rytov approximation is accurate when the square of the phase gradient is much less than the index contrast divided by the wavelength squared. The approximations are asymptotically equivalent for low-contrast, small objects. However, the areas of applicability of the two approximations is usually seen as mostly disjoint, because of the very different assumptions implicit in each about the relation of the field to the scatterer, either linear or exponential.

By using a common limiting form of the exponential function, it can be seen how the two are related. The Born series may be obtained by writing the field as u=u ₀+u _s , where u ₀ is the incident field and us is the scattered field and a perturbative expansion in orders of the contrast is made for u ₀. On the other hand, the Rytov approximation is derived by first expressing the total field u=u ₀ expϕ, where one identifies ϕ as a complex phase of the scattered field. By an identity,

one can rewrite the Rytov model as

Now consider removing the limit from Eq. (2), so that n becomes a parameter that can be varied from one to infinity. If n=1, then u=u ₀+u ₀ ϕ, which is the additive relationship assumed in the Born approximation, if we identify u _s =u ₀ ϕ. This suggests that by using Eq. (2) as a definition of a generalized complex phase where n can take any value from one to infinity, a range of models are available lying between the Born or Rytov approximations by varying the parameter n. It is possible that intermediate values of n may provide better accuracy in situations where neither the Born and Rytov approximations apply. Lu also considered generalized transformations [11, 12] on the field that yield the Born and Rytov approximations, including an intermediate transform that can be varied between Born and Rytov. Here we explore the validity and utility of a specific field transformation that may be especially useful and adaptable to current inverse scattering methods.

More importantly, such a generalized definition of the complex phase may be useful in inverse scattering problems where multiple scattering can not be neglected. It will be shown a complex phase defined with an intermediate value of n can be used to incorporate both the Born and Rytov models in the search for a solution to inverse scattering problems. This way, a solution to the multiple scattering problem can be consistently defined that uses both the Born and Rytov models. This will add a degree of flexibility that may enable inverse scattering problems to be solved where the scatterer has high contrast and large extent. For this reason, results pertaining to the application of the hybrid approximation to multiple scattering methods such as higher order terms of the hybrid series are included in this work.

In section 2, the generalized definition of the complex phase is introduced and from this the integral equation it must satisfy is derived. Based on this relationship, the forward scattering problem is linearized. Section 3 derives the higher order terms of the hybrid series which can be used for more accurate calculations of forward-scattering problems. The Fréchet derivative of the error in the integral equation is derived in section 4, which is useful in methods of nonlinear inverse scattering. It is shown that this derivative contains both the Born and Rytov components that can be used to search for solutions of the multiple scattering problem. In section 5, it is outlined how present inverse scattering methods can be adapted to the hybrid method. Finally, section 6 tests the first order hybrid approximation by numerically reconstructing the contrast function of a circular cylinder with varying sizes and contrast, and comparing the error in the reconstruction to the conventional Born and Rytov methods.

2. Derivation of a Hybrid Approximation

To derive the hybrid approximation discussed above, we study scattering of a scalar wave in an inhomogeneous medium. The field satisfies the spatially inhomogeneous reduced wave equation:

where k(r) is the spatial frequency of the wave in the medium at position r, and k ₀ is the freespace spatial frequency of the scalar field. We posit that the solution to this equation has the form

where u ₀(r) is an unperturbed field that satisfies the homogeneous wave equation ∇² u ₀+$k_0^2$ u ₀=0, and ϕ(r) is a “complex phase” that will be determined to find the scattered field. The function ϕ will be expanded as a perturbation series $\varphi \left(\mathbf{r}\right)={\displaystyle \sum _{m=0}^{\infty }}{\epsilon }^m{\varphi }_m\left(\mathbf{r}\right)$, where ε is the order parameter proportional to the magnitude of the index perturbation of the medium. Note that ϕ ₀=0 because the field u ₀ is defined as satisfying the unperturbed reduced wave equation, so that the lowest order non-zero term of ϕ is first order. The constant n will be a free parameter [1,∞]. This form is chosen so that u(r)=u ₀(r)+u ₀(r)ϕ(r) when n=1, which is the model used for Born scattering with u ₀(r)ϕ(r) being the scattered field. As n→∞, u(r)=u ₀(r)exp(ϕ(r)), which is the form of the total field in the Rytov approximation. By adjusting n, we can produce an approximation that is intermediate between the two approximations.

Substituting Eq. (4) into Eq. (3), we find that

Expanding out the differential operators, this becomes

If we note that ${\left(1+\frac{\varphi }{n}\right)}^n{\nabla }^2{u}_0+{\left(1+\frac{\varphi }{n}\right)}^n{k}_0^2{u}_0=0$ and subtract this relation from Eq. (6), and then divide the resulting equation by ${\left(1+\frac{\varphi }{n}\right)}^{n-1}$, one finds that

Note that we have inserted the order parameter ε to indicate that the index perturbation k ²-$k_0^2$ is of first-order, and defined the contrast κ=k ²-$k_0^2$ for brevity. Now we note that ∇²(u ₀ ϕ)=u ₀∇² ϕ+2∇u ₀·∇ϕ+ϕ∇² u ₀, and that ∇² u ₀=-$k_0^2$ u ₀, which combined yield ∇²(u ₀ ϕ)+$k_0^2$(u ₀ ϕ)=u ₀∇² ϕ+2∇u ₀·∇ϕ. Substituting this into Eq. (7) and rearranging terms, one finds that

Equation (8) may be recast as an integral equation, with the following result obtained

where G(r ^′,r) is the Green’s function of the homogeneous wave equation. Note that no approximations have been made yet. An interesting feature that follows from the derivation is that when n=1, and the quantity u ₀ is added to both sides of this equation, then this equation is the Lippmann-Schwinger equation. To linearize this integral to form the first-order approximation, we explicitly list the quantities summed inside the integral on the right hand side of Eq. (9) separately to find their respective orders of ε:

As explained earlier, the lowest ε-order of ϕ that is non-zero is first order. Equation (10) is the product of a term ${\left(1+\frac{\varphi }{n}\right)}^{-1}$ which is of order zero (ε ⁰), and ∇_ϕ·∇ϕ of order ε², so this term is ε² order. Equation (11) is of order ε¹, because it does not contain a ϕ, only an ε. Finally, Eq. (12) is of order ε ², because it contains the product of ϕ and ε. Therefore to first order only the term of Eq. (11) needs to be retained, so that to a first-order approximation Eq. (9) is

Note that this linearization does not contain the parameter n, so it has the same form for any value of n. This is why the first Born and Rytov approximations have the same form; only the definition of ϕ is different. Therefore n can be chosen as needed as long as the conditions of the approximation are satisfied. To derive these conditions, we require that the omitted terms produce a contribution to the complex phase much smaller than one:

Approximating u ₀(r) as a constant, and assuming that G(r ^′,r)≈|r-r ^′|^-1 is a worst-possible case of constructive interference between the fields scattered inside the volume, then this can be simplified to

We have discarded the third-order and higher terms in ϕ, denoted the average wave number perturbation by Δk (which is of order ε), and the average distance between points in the object (or alternately, the radius of a sphere enclosing the scatterer) as R. This inequality combines the constraints of both the Born and the Rytov approximations. The Rytov constraint, which is weighted by $\frac{n-1}{n}$, restricts the total integrated squared gradient of the complex phase over the scatterer volume. The Born constraint, which is weighted by $\frac{1}{n}$, restricts the total integrated complex phase magnitude over the volume. By adjusting n, one can trade off the amount of error in the reconstruction produced by either the Born or Rytov models, so that accuracy can be maintained when neither of the models separately applies.

3. Higher-Order Terms of the Hybrid Series

To derive the higher order terms of the hybrid series, we expand the complex phase as a power series, $\varphi \left(\mathbf{r}\right)={\displaystyle \sum _{m=0}^{\infty }}{\epsilon }^m{\varphi }_m\left(\mathbf{r}\right)$. The power series for ϕ is inserted into the following equation, which is Eq. (7) after the equation is multiplied by $\left(1+\frac{\varphi }{n}\right)⁄u_0$

After expanding out the products of the power series, and collecting terms, the following equation found for the coefficients of ε^p for p>1:

By rearranging the terms, and noting that (∇²+$k_0^2$)(u ₀ ϕ _m )=2∇u ₀·ϕ _m +u ₀∇² ϕ_m , the following recurrence relation is obtained for (∇²+$k_0^2$)(u ₀ ϕ_p ):

Evaluating the recurrence relation for the first three terms of the series:

These differential equations for u ₀ ϕ_p can be transformed into integral solutions using the Green’s function in a manner identical to the transformation between Eq. (8) and Eq. (9).

In an alternative expansion [6], we can expand the total field $u={\displaystyle \sum _{m=0}^{\infty }}{\epsilon }^m{u}_m$ as a Born series, and the hybrid series as $u=u_0{\left(1+\frac{\varphi }{n}\right)}^n$ where $\varphi ={\displaystyle \sum _{m=0}^{\infty }}{\epsilon }^m{\varphi }_m$ and equate terms of order ε ^p by expanding the right hand side of the equation with the binomial theorem. If this is done, the relationship between the first three Born and the hybrid series terms are:

Because the leading term of u _p is always u ₀ ϕ_p , the Born terms up to order p can be calculated, and then the term ϕ_p can be calculated from ϕ ₁ to ϕ _p-1 and u _p . This may be a more rapid way to calculate ϕ_p . These higher-order forward-scattering terms can be used to compute an inverse-scattering series as given in [13].

4. Evaluating the Fréchet derivative for Inverse Multiple Scattering

In addition to being useful for improving the accuracy of first-order scattering and inverse-scattering calculations, this hybrid approach may be useful to help account for multiple scattering. Methods such as the Distorted Born Iterative Method, the Born Iterative Method [1, 14], and the Modified Gradient Method [15, 16] attempt to solve the Lippmann-Schwinger equation by iteratively estimating the object contrast and the field scattered from the object in order to produce a solution consistent with the measured scattered field data. These methods numerically scatter off of the estimated object multiple times using the first Born approximation so that over many iterations the true scattered field inside the object is estimated. On the other hand, a method of inverting the Rytov series [13] has also been proposed, which may have an advantage for some diffraction tomography problems because of the different validity conditions for the Rytov approximation. Reference [17] compares various inverse multiple scattering methods which is helpful to understand these approaches. The hybrid approximation may have advantages for solving multiple scattering problems because it incorporates both the Born and Rytov series, so that a combination of both can be used each iteration to search for a solution. To explain this, we propose a method of inverse scattering that seeks to minimize the error in the solution of Eq. (9) for ϕ and κ, which is the analogue of the Lippmann-Schwinger equation for the hybrid complex phase. This method is akin to the Modified Gradient Method, but allows an arbitrary n to be used rather than the Born approximation for which n=1. The functional that will be minimized to estimate a solution to Eq. (9) is L(ϕ,κ) and is given by

This functional represents the γ-norm of the total error in Eq. (9) integrated over the entire scatterer volume, where the exponent γ determines the norm to be minimized. The quantity M(r′) is the error in Eq. (9) at point r′ The equation has been further divided into M _R and M _B components which is the Rytov and Born contributions to the error M. The function q() is an optional regularization term constraining the reconstructed κ and potentially its gradient. Minimizing L finds the minimum-norm solution to Eq. (9) for ϕ and κ. To find this minimum, we find the Fréchet derivatives of L with respect to these quantities so that a gradient descent method can be used to minimize L. These derivatives can be computed using the Euler-Lagrange derivative formula $\frac{\delta \left({\mid M\mid }^{\gamma }\right)}{\delta \varphi }=\frac{\partial \left({\mid M\mid }^{\gamma }\right)}{\partial \varphi }-{\nabla }_{\mathbf{r}\text{'}}.\frac{\partial \left({\mid M\mid }^{\gamma }\right)}{\partial (\nabla \varphi )}$:

where

The functional derivative $\frac{\delta \mid M\mid \gamma }{\delta \varphi }$ has been divided into two components: the Rytov component M _R , which is weighted by $\frac{n}{n-1}$, and the Born component M _B , which is weighted by $\frac{1}{n}$. These two components define two search directions that locally point toward solutions for ϕ that conform to the assumptions of either approximation, and the derivative is a weighted sum of these. One should note that the second term of $\frac{\delta M_R}{\delta \varphi }$ is a cross term that is only present when n is not 1 or ∞, so it effectively couples the Born and Rytov approximations together in a way that would not be achieved if the Born and Rytov search directions were used separately. The complex phase, as defined in Eq. (4), may be advantageous over the conventional Born and Rytov definitions because it enables one to find a solution where both the Born and Rytov search directions can be used consistently to find a solution to Eq. (9).

With the gradient of Eq. (25) and the other functional gradient $\frac{\delta \left({\mid M\mid }^{\gamma }+q\right)}{\delta \kappa }$ given by

one can use gradient descent or conjugate gradient methods to find a sequence of ϕ(r) and κ(r) that successively minimize L. With both the Born and Rytov search directions available, there is greater flexibility in finding a solution that can solve the multiple scattering problem without stagnation. Unlike other methods such as the Born Iterative Method or Newton-Kantorovich, because Eq. (9) is nonlinear in ϕ, a successive approximation approach may be more efficient than directly solving for κ(r) and ϕ(r) each iteration, and it can allow the regularization to be introduced in a consistent way.

5. Inverse Scattering using the Hybrid Series

Implementing this method for inverse scattering is similar to implementing the first Rytov approximation in practice. For the Rytov approximation, measurements of the total field u(r) are made, and from these the complex phase can be computed by ϕ(r)=log(u/u ₀). Because the log function is multiple-valued for complex numbers, there is an ambiguity in ϕ of a multiple of 2πi. To account for this, typically the measured value of ϕ is phase unwrapped as a function of position of the scattered field. The unwrapping is done such that the difference between the imaginary parts of two adjacent samples of ϕ differs by less than π. This complex phase is used to solve Eq. (13) for κ(r).

In the hybrid approximation, the complex phase ϕ is given by

There is a multiple-valued inverse associated with taking the n^th root of a complex number. This leads to an unwrapping problem similar to that which exists in the Rytov approximation, because there are n roots of a given complex number if n is an integer. To properly unwrap the n^th root of a complex function f (x), one can write f(x)=|f(x)|exp[iθ(x)], where θ(x)=arg f(x), wrapped from [-π,π]. If θ^′ (x) is the unwrapped θ(x), then the unwrapped f(x)^1/n=|f(x)|^1/n exp[iθ ^′(x)/n]. Alternatively, one can unwrap the logarithm of $\frac{u}{u_0}$ as would be done in the Rytov case, and use the second form of Eq. (27).

With this change in unwrapping, the inverse scattering solution proceeds the same as it would in the Rytov case. The only change is the defintion of the complex phase and how it is unwrapped. Therefore this method should be easy to implement in cases where the Rytov approximation is already used.

As an example of implementing higher order diffraction tomography using the hybrid series, we adapt the method of inversion of nonlinear Volterra operators, which has been applied to the Born series [18], Rytov series [13], and the eikonal equation [19]. Below we outline the steps of estimating the inhomogeneity of an object with a nonlinear Volterra operator modified to use the hybrid series. For this, several operators will need to be defined. One is a repeated forward-scattering operator G _i given by

which computes the i^th Born-series term u ₀ G _i (κ) of the field u ₀ scattered from an inhomogeneity κ. The integration can often be performed using a Fast Fourier Transform if the kernel G(r ^′,r) is shift-invariant. Next, operators are needed that calculate the terms of the hybrid series (based on Eqs. (20–22)):

In general, there will be many incident fields u ₀ corresponding to many illumination source locations for which the scattered phases ϕ will be measured or computed. Braces around a phase {ϕ} denote the set of all phase functions ϕ that are computed or measured for all of the incident fields u ₀.

The linearized inverse of G ₁(κ) is denoted by the operator B({ϕ}) which computes an estimate of κ from all {ϕ} such that B[{G ₁(κ)}]≈κ. This inverse corresponds to the linearized inverse diffraction tomography operator, implemented by the Fourier diffraction theorem [2] and the filtered backpropagation algorithm [20] methods, among others. Regularization of this operator may be needed to ensure that the inverse is stable.

To implement the inverse, the complex phase {ϕ} is calculated from the measured fields {u} using Eq. (27) including phase unwrapping. From this data function, the first two inverse terms of the hybrid series are calculated as follows:

The total object estimate is then given as κ=κ ₁+κ ₂. As a potentially simpler method, one may wish to explore applying the nonlinear inverse methods of [21] to this series.

6. Simulation

To test the utility of the hybrid approximation, a two-dimensional simulation of a diffraction tomography inverse scattering experiment was performed using the first-order hybrid approximation. The error in the hybrid reconstruction is compared to the error of the first-order Born and Rytov approximations. Simulation computer codes [22] were obtained and modified to implement the hybrid reconstruction, so that the details of the simulation are almost identical to that of [7] and [2]. The simulated objects were circular cylinders of varying contrasts and sizes. The circular cylinder was a useful test object because there is an exact solution for plane wave scattering from a cylinder. The simulated diffracted field from the cylinder was computed on a receiver line located at a distance of 100 wavelengths from the center of the cylinder sampled every 1/4 wavelength, so that 400 points are sampled in the detector field. Only one projection needed to be computed because because of the cylindrical symmetry of the object, and the Fourier space was sampled as a 256×256 grid. Fourier reconstruction was implemented using the Fourier Diffraction Theorem [2] including interpolation in the Fourier domain. The only change in how the inverse was implemented from the original code was to employ Eq. (27) to compute the complex phase of the scattered field, with the unwrapped phase included in the imaginary part of the complex logarithm of log $\frac{u}{u_0}$.

The results of the simulations are shown in Figs. 1 and 2. Each of the figure parts show the computed reconstructions of a cylinder, varying the radii and refractive indices. The reconstructed index of each cylinder as a function of radius is plotted, with the blue line corresponding to the first Born approximation, the magenta to the first Rytov approximation, and the black line to the first hybrid approximation. For each reconstruction, the root-mean-square (RMS) percentage error was computed between the computed reconstruction and the ideal cylinder profile, which is indicated in the inset of each graph. The inset of each graph also specifies the radius and refractive index of its respective cylinder, and the exponent n for which the RMS error was minimized for the hybrid approximation.

The tendency appears to be that the Born approximation reproduces the boundary well but not the interior, while the Rytov fills the interior but removes the boundary. By choosing the optimal hybrid exponent, these two tendencies can be balanced and the result it a more uniformly reconstructed cylinder. It appears that the hybrid reconstruction can achieve the best improvement over both the Born and Rytov reconstructions when the optimal exponent is close to two, so that n=2 may be a good starting point to find the optimal exponent. In addition to this figure, animations are available on-line that show the evolution of the computed reconstruction of three of these cylinders from the Born to the Rytov solutions. [Media 1, Media 2, Media 3]

Based on the reconstructions shown in these figures, empirical formulas for the optimal exponent and RMS error for that exponent are fitted:

with N being the refractive index of the cylinder, and R being the radius of the cylinder in wavelengths. The equation for n _opt suggests that perhaps the optimal exponent is proportional to the product of the index contrast and the area of the object. However, with a sufficiently large index difference and object area neither the first Born, Rytov, or the hybrid approximation can be expected to provide accurate results. In compiling the figures, hybrid reconstructions that produced optimal RMS errors greater than 40% were excluded from the table because it was decided that these points produced results too inaccurate to be of interest. Therefore it can not be expected that the empirical formula will be accurate when the RMS error is over 40%.

We have derived a hybrid approximation that incorporates the advantages of both the Born and Rytov approximations. In cases where either of these approximations individually may not apply, this hybrid method helps obtain additional accuracy, but retains the convenient frame-work of linear inverse scattering. Furthermore, we show that the Born and the Rytov method are two endpoints of a family of linearized scattering approximations.

 

Fig. 1. Profile plots of computed reconstructions of cylinders with various refractive index and radii using the Born, Rytov, and hybrid approximations. The blue curve is the Born reconstruction of the refractive index contrast, the magenta curve the Rytov reconstruction, and the black curve the hybrid reconstruction. The text in the upper left corners of each subfigure (a) to (j) specifies the radius and index contrast of each cylinder, and the optimal exponent for which the error was minimized. The text in the upper right corners is the RMS error for the Born, Rytov, and hybrid reconstructions.

Download Full Size | PDF

 

Fig. 2. More profile plots of computed reconstructions of cylinders with various refractive index and radii using the Born, Rytov, and hybrid approximations. The blue curve is the Born reconstruction of the refractive index contrast, the magenta curve the Rytov reconstruction, and the black curve the hybrid reconstruction. The text in the upper left corners of each subfigure (a) to (h) specifies the radius and index contrast of each cylinder, and the optimal exponent for which the error was minimized. The text in the upper right corners is the RMS error for the Born, Rytov, and hybrid reconstructions.

Download Full Size | PDF