Apparatus and method for imaging objects with wavefields

[0001] This is a continuation of U.S. patent application Ser. No. 10/615,569, filed Jul. 7, 2003, now U.S. Pat. No. 7,570,742 which is a continuation of U.S. patent application Ser. No. 10/024,035, filed on Dec. 17, 2001, now U.S. Pat. No. 6,636,584, which is a continuation of U.S. patent application Ser. No. 09/471,106, filed on Dec. 21, 1999, now U.S. Pat. No. 6,587,540, which is a continuation-in-part of U.S. patent application Ser. No. 08/706,205, filed on Aug. 29, 1996, now abandoned, which is a continuation-in-part of U.S. patent application Ser. No. 08/486,971, filed on Jun. 22, 1995, now abandoned, which is a continuation-in-part of U.S. patent application Ser. No. 07/961,768, filed on Oct. 14, 1992, now U.S. Pat. No. 5,588,032; which are all herein incorporated by reference.

[0002] Material on one (1) compact disc accompanying this patent and labeled COPY 1, a copy of which is also included and is labeled as COPY 2, is incorporated herein by reference. The compact disc has a software appendix thereon in the form of the following three (3) files:

[0003] Appendix A.txt Size: 16Kbytes, Created: Aug. 13, 2001;

[0004] Appendix B.txt, Size: 19Kbytes, Created: Aug. 13, 2001;and

[0005] Appendix C.txt, Size: 1Kbytes, Created: Aug. 13, 2001.

[0006] A portion of the disclosure of this patent document contains material to which a claim of copyright protection is made. The copyright owner(s) has/have no objection to the facsimile reproduction by any one of the patent document or the patent disclosure as it appears in the Patent and Trademark Office patent file or records, but reserves all other rights in the copyrighted word.

[0007] It has long been known that elastic, electromagnetic or acoustic waves in homogeneous and layered environments in the frequency range of a fraction of a cycle per second up to hundreds of millions of cycles per second and higher can be propagated through many solids and liquids. Elastic waves are waves that propagate through solids, and have components of particle motion both parallel (longitudinal, or pressure, wave) and perpendicular (shear wave) to the direction of propagation of the wave energy itself. In contrast to this, acoustic waves are those waves that generate particle motion that is exclusively parallel to the propagation of wave energy. Electromagnetic waves have components of variation of field strength solely in the direction perpendicular to the direction of propagation. All of these types of waves may be used to image the acoustic longitudinal wavespeed and absorption, the electromagnetic wavespeed and absorption, the shear wavespeed, and the density of the material through which the wave energy has traveled.

[0008] It is also known that scattering is produced not only by spatial fluctuations in acoustic impedance, which is the product of mass density times wavespeed, but also by independent fluctuations in electromagnetic permeability, permittivity and conductivity, elastic compressibility, shear modulus, density, and absorption. These lead to variations in phase speed (which is the speed of propagation of fronts of constant phase) and in impedance (for the electromagnetic case, the ratio of the electric to the magnetic field strength). The net property of an object which describes the phenomenon of scattering in a given modality, is called the “scattering potential”.

[0009] Other imaging methods have been applied to one or the other modality, or restricted to acoustic or elastic media, the method described in this patent is applicable to any type of wave motion, whether electromagnetic, elastic (including shear wave effects) or acoustic (scalar approximation valid in liquid and gases).

[0010] Furthermore, the ambient media may have some forms of structure (layering) or microstructure (porosity) relevant to the medical, geophysical, or nondestructive imaging applications envisioned for this technology. In the prior art the presence of this layering or porosity has greatly diminished the effectiveness of the imaging program. The method of this patent minimizes the obscuring effect of such structures in the ambient media. In addition, we have made several changes to the previous U.S. Pat. No. 4,662,222 that significantly extend the applicability and speed of our algorithm. These changes are described, in part, below:

[0011] As discussed in U.S. Pat. No. 4,662,222, by one of the present authors, these principles can be used to image scattering bodies within a given medium. Some definitions will clarify the procedures and methods used:

[0012] The direct or forward scattering problem is concerned with a determination of the scattered energy or fields when the elastic or electromagnetic properties of the scattering potential are known.

[0013] The inverse scattering problem, by contrast, consists in the use of scattered electromagnetic, elastic, or acoustic waves to determine the internal material properties of objects embedded in a known (ambient) medium. An incident wave field is imposed upon the ambient medium and the scatterer. The scattered field is measured at detectors placed a finite distance from the scattering objects. The material parameters of the scatterer are then reconstructed from the information contained in the incident and scattered fields. In other words, as defined herein, acoustic or electromagnetic imaging using inverse scattering techniques is intended to mean electronic or optical reconstruction and display of the size shape, and unique distribution of material elastic or electromagnetic and viscous properties of an object scanned with acoustic, electromagnetic or acoustic energy, i.e., reconstruction of that scattering potential which, for a given incident field and for a given wave equation, would replicate a given measurement of the scattered field for any source location.

[0014] The description of prior art summarized in U.S. Pat. No. 4,662,222 (referred to as “the previous patent”) is also relevant here. It is sufficient to say here, that the use of the conjugate gradient algorithms discussed in the previous patent brought the attainment of inverse scattering into a practical reality. Previous to that patent (U.S. Pat. No. 4,662,222), only linear approximations to the full inverse scattering problem were amenable to solution, even with the aid of modern high speed digital computers. The resolution achievable with the inverse scattering methods described in the previous patent was far superior to any ultrasound imaging in medicine previous to that time.

[0015] This invention is designed to provide improved imaging of bodies using wave field energy such as ultrasound energy, electromagnetic energy, elastic wave energy or Biot wave energy. In particular it is designed to provide improved imaging for medical diagnostic applications, for geophysical imaging applications, for environmental imaging applications, for seismic applications, for sonar applications, for radar applications and for similar applications where wave field energy is used to probe and image the interior of objects.

[0016] In particular, this invention has important applications to the field of medical diagnosis with specific value to improved breast cancer detection and screening. The present method of choice is mammography. This choice is not supported by the outstanding performance of mammography, but rather because it is the most diagnostically effective for its present cost per exam. Another reason for the wide spread use of mammography is the inertia of changing to other modalities. It is known that hand-probe-based, clinical reflection ultrasound can match the performance of mammography in many cases, but only in the hands of specially trained ultrasound radiologists. Today and in the foreseeable future, there are more mammography machines and radiologists trained to use them than there are trained ultrasound radiologists that have access to high quality ultrasound scanners. Sophisticated breast cancer diagnostic centers use both mammography and ultrasound to compensate for weakness in either approach when used alone. When used alone, ultrasound and mammography require a biopsy to remove a sample of breast tissue for pathology lab analysis to determine whether a sampled lesion is cancerous or benign.

[0017] Even when mammography and ultrasound are used together, the specificity for discriminating between cancer and fibrocystic condition or between cancerous tumor and a fibroadenoma is not high enough to eliminate the need for biopsy in 20 to 30 percent of lesions. Given that early diagnosis of breast cancer can ensure survival and that one woman in eight will have breast cancer in her life, it is important to the general population for cancer to be detected as early as possible. Detection of cancer in an early stage for biopsy is thus very important. However, biopsy of benign lesions is traumatic to the patient and expensive. A mammogram costs about $90 but a biopsy is about ten times more expensive. Thus, it is important that a breast cancer diagnostic system have high specificity and sensitivity to eliminate unnecessary biopsies. Increasing the rate of diagnostic true positives to near 100 percent will identify all lesions as a cancer that are cancer without the need to biopsy. However, it is also necessary to increase the rate of true negatives to near 100 percent to eliminate biopsy of benign lesions. Neither mammography or ultrasound or their joint use has provided the combination of sensitivity and specificity to eliminate biopsy or to detect all cancers early enough to insure long term survival after breast cancer surgery for all women that have had breast exams.

[0018] There does not seem to be any obvious improvements in present mammography or hand-probe-based, clinical reflection ultrasound that can significantly improve these statistics. However, there is reason to believe that inverse scattering ultrasound tomography or electric impedance tomography can provide improved diagnostic sensitivity or specificity when used separately or jointly with themselves and with ultrasound and/or mammography. Inverse scattering ultrasound imaging provides several advantages over present clinical reflection ultrasound that uses hand held ultrasound probes. A hand held probe is a transducer assembly that scans an area of the patient's body below the probe. Probes consisting of mechanical scanning transducers and probes comprising electronically scanned arrays of transducer elements are both well developed technologies.

[0019] Inverse scattering has the following advantages over said clinical reflection ultrasound. Inverse scattering images have the following features: (1) two separate images can be made, one of speed of sound and one of acoustic absorption; (2) these images are a quantitative representation of actual acoustic bulk tissue properties; (3) the images are machine independent; (4) the images are operator independent; (5) the images are nearly free of speckle and other artifacts; (6) all orders of scattering tend to be used to enhance the images and do not contribute to artifacts or degrade the images; (7) the images of speed of sound and absorption create a separation of tissues into classes (fat, benign fluid filled cyst, cancer, and fibroadenoma), although the cancer and fibroadenoma values tend to overlap; (8) the speed of sound and acoustic absorption images have excellent spatial resolution of 0.3 to 0.65 mm at or below 5 MHz; and (9) the speed of sound images can be used to correct reflection tomography for phase aberration artifacts and improve ultrasound reflectivity spatial resolution. Inverse scattering easily discriminates fat from other tissues, while mammography and present clinical ultrasound cannot.

[0020] Because of the similar values of speed of sound and acoustic absorption between cancer and fibrocystic condition (including fibroadenoma), it is not known whether inverse scattering will provide the required high level of specificity to eliminate biopsy. Perhaps this performance could be achieved if inverse scattering were combined with reflection ultrasound (including reflection tomography) or with mammography or with both.

[0021] A traditional problem with inverse scattering imaging is the long computing times required to make an image. Diffraction Tomography is an approximation to inverse scattering that uses a first order perturbation expansion of some wave equation. Diffraction Tomography is extremely rapid in image computation, but suffers the fatal flaw for medical imaging of producing images that are blurred and not of acceptable quality for diagnostic purposes and it is not quantitatively accurate. In our last patent application, we addressed the computational speed problem with inverse scattering and showed how to increase the calculation speed by two orders of magnitude over finite difference methods and integral equation methods. This speed up in performance was achieved by use of a parabolic marching method. This improvement in speed was sufficient to allow single 2D slices to be imaged at frequencies of 2 MHz using data collected in the laboratory. However, making a full 3-D image with a 3D algorithm or even from several, stacked 2-D image slices made with a 2D algorithm would have required computing speed not available then and only becoming available now.

[0022] Another imaging modality that has been investigated for detecting breast cancer is EIT (electrical impedance tomography). This modality has been investigated for many years and many of its features are well known. One of its great advantages for breast cancer detection is its high selectivity between cancer and fibrocystic conditions including fibroadenoma. Cancer has high electrical conductivity (low electrical impedance) while fibrocystic conditions and fibroadenoma have low electrical conductivity (high electrical impedance). However, EIT has poor spatial resolution. Also EIT requires the use of inversion algorithms similar (in a general sense) to those of inverse scattering. In addition EIT algorithms have mostly been used to make 2-D images. Work on making 3-D EIT images is less developed because of the increased computer run time of 3-D algorithms.

[0023] Other problems with mammography may be listed, such as the pain associated with compressing the breast between two plates to reduce the thickness of the breast in order to improve image contrast and cancer detection potential. Another problem is the accumulated ionizing radiation dose over years of mammography. It is true that a single mammogram has a very low x-ray dose, especially with modern equipment. However, if mammograms are taken every year from age 40 or earlier, by age 60 to 70 the accumulated dose begins to cause breast cancer. The effect is slight, and the benefits of diagnosis outweigh this risk. Nevertheless, it would be an advantage to eliminate this effect. Mammography is not useful in many younger women (because of the tendency for greater breast density) or older women with dense breasts (such as lactating women). About 15 percent of women have breasts so dense that mammography has no value and another 15 percent of women have breasts dense enough that the value of mammography is questionable. Thus 30 percent of all mammograms are of questionable or zero value because of image artifacts from higher than normal breast density.

[0024] The present invention describes a transmission wave field imaging method, comprising the transmission of an incident wave field into an object, the incident wave field propagating into the object and, at least, partially scattering. The invention also describes the measuring of a wave field transmitted, at least in part, through the object to obtain a measured wave field, the measured wave field based, in part, on the incident wave field and the object. The invention further describes the processing of the measured wave field utilizing a recursive reconstruction algorithm to generate an image data set representing at least one image of the object.

[0025] The imaging method above also describes a recursive reconstruction algorithm that calculates a predicted total wave field utilizing a single independent variable, the single independent variable representing scatter potential. The recursive reconstruction algorithm represents a gradient-based recursive algorithm that solves for an independent scatter potential variable, the gradient-based recursive algorithm treating the total wave field as a dependent variable depending upon the independent scatter potential variable. The processing calculates a predicted total wave field based on the incident wave field and generates the image data set based on a comparison of the predicted total wave field and the measured wave field.

[0026] The imaging method above also describes the transmitting where transmitting includes multiple incident wave fields from associated multiple source positions. The incident wave field includes multiple frequency components, the recursive reconstruction algorithm separately processing the multiple frequency components.

[0027] The imaging method above also describes an incident wave field representing an ultrasound wave field. The method incorporates the displaying of at least one image as at least two images co-displayed, the at least two images including a speed-of-sound image and an attenuation image. The incident wave field represents an ultrasound wave field and the processing including generating a speed-of-sound data set and an attenuation data set, the speed-of-sound data set representing speed of the ultrasound wave field through the object, the attenuation data set representing attenuation of the ultrasound wave field within the object. The incident wave field comprises at least one of acoustic energy, elastic energy, electromagnetic energy and microwave energy.

[0028] The imaging method above also describes the object as a breast, and a transmission source being tilted to direct the incident wave field toward the chest wall of a patient to image features close to the chest wall. The transmitting and receiving are performed by a transmitter and receiver, respectively, the transmitter and/or the receiver are moved with respect to the object.

[0029] The imaging method above also describes repeating the transmitting and receiving at multiple slices within the object to generate a 3D volume of image data sets, displaying 2-D images representative of slices through the object, and displaying a 3-D image of at least a portion of the object. The measuring is performed at multiple positions about the object.

[0030] The present invention also describes transmission wave field imaging method wherein an incident wave field is transmitted into an object, the incident wave field propagates into the object and, at least, partially scatters. The invention also includes the measuring of the total wave field transmitted, at least in part, through the object to obtain a measured wave field, the measured wave field based, in part, on the incident wave field and the object. Furthermore, the method describes processing the measured wave field utilizing a gradient-based reconstruction algorithm to generate an image data set representing at least one image of the object. The gradient-based algorithm solves for an independent scatter potential variable and treats the measured wave field as a dependent variable depending upon the independent scatter potential variable.

[0031] The present invention also describes a transmission wave field imaging method wherein an ultrasound incident wave field is transmitted into an object, the incident wave field is propagated into the object and, at least, partially scatters. The invention further describes the measuring of an ultrasound total wave field transmitted, at least in part through the object wherein the measured wave field is defined, in part, by the incident wave field and the object. Furthermore, the invention obtains speed-of-sound and attenuation values which correspond to the speed of the incident wave field through the tissue. The attenuation values correspond to an amount of attenuation of the incident wave field which is propagated through the tissue. The invention describes the construction of a speed-of-sound image from the speed-of-sound values, an attenuation image from the attenuation values, and the displaying of the speed-of-sound image and the attenuation image.

[0032] The imaging method described above also calculates a scattering potential based on the measured wave field wherein the attenuation image and the speed-of-sound image are representative of a component of the scattering potential. The speed-of-sound and attenuation values represent ranges of gray scale values for pixels comprising the speed-of-sound and attenuation images.

[0033] The imaging method described above utilizes an inverse scattering algorithm to reconstruct the speed-of-sound and attenuation images based on the measured wave field. The processing calculates a predicted total wave field based on the incident wave field and the constructing generates the speed-of-sound and attenuation images based on a comparison of the predicted total wave field and the measured wave field.

[0034] The present invention also describes a method for displaying transmitted ultrasound breast images. Transmitted ultrasound data sets are obtained representative of a measured wave field transmitted, at least in part, through an object, wherein the measured wave field is defined, in part, by an ultrasound incident wave field, a scattering potential and the object through which the wave is transmitted. Speed-of-sound and attenuation values are also obtained. The speed-of-sound values correspond to speed of ultrasound incident wave fields through tissue, the attenuation values correspond to an amount of attenuation of ultrasound incident wave fields propagated through the tissue. Speed-of-sound images are constructed from the speed-of-sound values and attenuation images are constructed from the attenuation values. The method further describes the displaying of those images.

[0035] The present invention also describes a transmission wave field imaging device. The device comprises a transmitter for transmitting an incident wave field into an object wherein the incident wave field is propagating into the object and, at least, partially scatters. The device includes a receiver for measuring the total wave field transmitted, at least in part, through the object to obtain a measured wave field. The measured wave field is based, in part, on the incident wave field and the object. The device also includes a processing module for processing the measured wave field. The module utilizes a recursive reconstruction algorithm to generate an image data set representing at least one image of the object.

[0036] The imaging device above also describes a processor module comprising a plurality of processors operated in parallel and a PCI bus interconnecting the processing module and the transmitter and/or the receiver. The processor module is implemented on either a personal computer, a hand-held device, a multi-processor system and/or a network PC. The transmitter operates in a pulse echo mode to perform reflective tomography and is rotated about the object to obtain measured data at multiple receive positions around the object. The measured data defines the measured wave field and a display for displaying the image based on the image data set uses the speed of sound image and/or the attenuation image.

[0037] The imaging device above also describes the receiver and/or the transmitter as including transducer elements arranged in at least one of a circular array, a linear array, a 2D planar array and a checkerboard 2D array. The transmitter sequentially transmits first and second incident wave fields into the object having different first and second fundamental frequencies, respectively and sequentially transmits a high frequency incident wave field and a low frequency incident wave field.

[0038] The present invention also describes a transmission wave field imaging system. The system comprises a transmitter for transmitting an incident wave field into an object wherein the incident wave field is propagated into the object and, at least, partially scatters. The system includes a receiver measuring a total wave field transmitted, at least in part, through the object to obtain a measured wave field. The measured wave field is defined, in part, by the incident wave field, a scattering potential and the object. The measured wave field includes at least one of speed of sound information and attenuation information characterizing the object. The system also includes a processing module for processing, on-site and in real-time, the measured wave field to generate an image of the object. The image comprises pixel values based on at least one of the speed of sound information and the attenuation information characterizing the object. The system further comprises a display for displaying, on-site and in real-time, the image containing at least one of the speed of sound information and the attenuation information of the object.

[0039] The wave field imaging system above also describes a breast scanner including the transmitter and receiver, and the processing module that generates the image substantially during a breast scan. The processor module calculates a scattering potential from the measured wave field. The scattering potential includes a real component and an imaginary component, the real component defining the speed of sound information and the imaginary component defining the attenuation information.

[0040] The wave field imaging system above also describes a phase detector joined to the receiver wherein the phase detector outputs real and imaginary signal components. The transmitter and receiver repeat the transmitting and receiving operations at multiple slices within the object to collect a 3D volume of data and the display displays 2-D images representative of slices through the object. The display also displays a 3-D image of at least a portion of the object.

[0041] The wave field imaging system above also describes a movable carriage supporting the transmitter and receiver wherein the movable carriage rotating about an axis of a cylindrical chamber is configured to receive the object. The movable carriage supporting the transmitter and receiver is vertically adjusted along an axis of a cylindrical chamber. The chamber is configured to receive the object and a fluid tank is configured to receive at least one of the object, the transmitter and receiver.

[0042] Additional objects and advantages of the invention will be set forth in the description which follows, and in part will be obvious from the description, or may be learned by practice of the invention. The objects and advantages of the invention may be realized and obtained by means of the instruments and combinations particularly pointed out in the appended claims. These and other objects and features of the invention will become more fully apparent from the following description and appended claims, or may be learned by the practice of the invention as further set forth hereinafter.

[0119] Definition 1. An inverse scattering imaging method is any imaging method which totally or partially uses wave equation modeling and which utilizes a nonlinear operator relating the inversion data to the image components.

[0120] Definition 2. The inversion data is defined by that transformation of the wavefield data which is the range of the nonlinear operator. The domain of the nonlinear operator is the image components.

[0121] Definition 3. A parabolic propagation step is an arithmetic operation which computes the wavefield on a set of surfaces given the wavefield on another, disjoint set of surfaces using a parabolic differential (or difference?) equation in the derivation of the arithmetic operation.

[0122] Definition 4. By λ we mean the wavelength in the imbedding medium for the maximum temporal frequency in the wavefield energy.

[0123] Definition 5. By Δ we mean the average spatial separation of the points used to discretize the discretized wavefield in the computer implementation of the method.

[0124] Definition 6. By a short convolution operation we mean a discrete convolution which is faster if done by direct summation rather than by FFT.

[0125] The following explanation of notation is given to facilitate obtaining an understanding of the invention. Unless otherwise noted, the notation will essentially be the same as U.S. Pat. No. 4,662,222. The scattering potential γ changes from point to point within an object or body as well as changing with frequency of the incident field. Thus
γ_ωj≡γ_ω(x_j), x_j∈R³
or R³ is used to signify the scattering potential at pixel j or point j, for the incident field at frequency ω.γ_ω can be considered to be the vector composed of all values of j, i.e. for all pixels:
γ_ω≡(γ_ω1γ_ω2 . . . γ_ωN)^T,
where T indicates “transpose”. For our purposes of exposition and simplicity, we will hence forth consider the case where γ is independent of frequency, although this is not necessary. This convention also is essentially the same as that employed in the U.S. Pat. No. 4,662,222.
Notation
Vector Field Notation

[0126] The following notation will be used throughout the patent to denote a vector field describing elastic waves, or electromagnetic waves:

[0127] $f\left(r\right)\equiv \left(\begin{array}{c}{f}_x\left(x,y,z\right)\\ f_y\left(x,y,z\right)\\ f_z\left(x,y,z\right)\end{array}\right)\equiv {\left(f_x{f}_y{f}_z\right)}^T$
(T-denotes “transpose”) is used to represent the total field.

[0128] $f^i\left(r\right)\equiv \left(\begin{array}{c}{f}_x^i\left(x,y,z\right)\\ f_y^i\left(x,y,z\right)\\ f_z^i\left(x,y,z\right)\end{array}\right)\equiv {\left(f_x^i{f}_y^i{f}_z^i\right)}^T$
(T-denotes “transpose”) is used to denote the incident field. The incident field is the field that would be present if there was no object present to image. In the case of layering in the ambient medium, it is the field that results from the unblemished, piecewise constant layering.

[0129] The scattered field is the difference between the total field and the incident field, it represents that part of the total field that is due to the presence of the inhomogeneity, i.e. the “mine” for the undersea ordnance locater, the hazardous waste canister for the hazardous waste locator, the school of fish for the echo-fish locator/counter, and the malignant tumor for the breast scanner:

[0130] $f^s\left(r\right)\equiv f\left(r\right)-f^i\left(r\right)$ $f^s\left(r\right)\equiv \left(\begin{array}{c}{f}_x^s\left(x,y,z\right)\\ f_y^s\left(x,y,z\right)\\ f_z^s\left(x,y,z\right)\end{array}\right)\equiv {\left(f_x^s{f}_y^s{f}_z^s\right)}^T$
T-denotes “transpose” f_ωφ^inc(r) denotes the scalar incident field coming from direction (source position) Φ at frequency ω. The r could represent either a 3 dimensional, or a 2 dimensional vector of position.
Scalar Field Notation

[0131] A scalar field is indicated by a nonbold f(r), r∈R³

[0132] This first example is designed to give the basic structure of the algorithm and to point out why it is so fast in this particular implementation compared to present state of the art. The background medium is assumed to be homogeneous medium (no layering). This example will highlight the exploitation of convolutional form via the FFT, the use of the Frechet derivative in the Gauss-Newton and FR-RP algorithms, the use of the biconjugate gradient algorithm (BCG) for the forward problems, the independence of the different view, and frequency, forward problems. It will also set up some examples which will elucidate the patent terminology. The field can be represented by a scalar quantity f. The object is assumed to have finite extent The speed of sound within the object is c=c(x). The speed of sound in the background medium is the constant c_o. The mass density is assumed to be constant. The attenuation is modeled as the imaginary part of the wavespeed. These are simplifying assumptions, which we make so as to focus on the important aspects of the imaging algorithm. By no means is the imaging algorithm restricted to the scalar case, either in theory, or, more importantly, in practical implementation.

[0133] We consider the two dimensional case first of all (equivalently the object and incident fields are assumed to be constant in one direction). Given an acoustic scatterer with constant mass density and constant cross-section in z, illuminated by a time harmonic incident field, f_ωφ^inc, with e^iωt dependence and source position (incident angle), φ, the total field satisfies the following integral equation.

[0134] $\begin{array}{cc}{f}_{\omega \varphi }^{\mathrm{inc}}\left(\stackrel{\_}{\rho }\right)=f_{\omega \varphi }\left(\stackrel{\_}{\rho }\right)-\int \int \gamma \left({\stackrel{\_}{\rho }}^{\prime }\right)f_{\omega \varphi }\left({\stackrel{\_}{\rho }}^{\prime }\right)g_{\omega }\left(1\stackrel{\_}{\rho }-\stackrel{\_}{\rho }1\right)d{x}^{\prime }d{y}^{\prime }\mathrm{where}\gamma =\frac{c_0^2}{c^2}-1& \left(1\right)\end{array}$

[0135] where

[0136] ρ=(x, y) and where f≡f_ωφ is the field internal to the object at frequency ω and resulting from the incident field from source position φ,: f_ωφ^inc. The 2-D Helmholtz Green's function is given by:

[0137] $g_{\omega }\left(1\stackrel{\_}{\rho }-{\stackrel{\_}{\rho }}^{\prime }1\right)=\frac{k_{{\omega }^2}}{4i}{H}_0^{\left(2\right)}\left(k_{\omega }l\stackrel{\_}{\rho }-{\stackrel{\_}{\rho }}^{\prime }l\right),k_{\omega }=\frac{\omega }{c_0}$
where H₀⁽²⁾ is the Hankel function of the second kind, and zeroth order.

[0138] Now it is required to compare the field measured at the detectors, with the field as predicted by a given guess, γ⁽ⁿ⁾ for γ. To accomplish this, first define the scattered field at the detectors as the total field minus the incident field.
f_ωφ^sc( ρ)≡f_ωφ( ρ)−f_ωφ^inc( ρ)

[0139] This represents the field due entirely to the scattering potential. Using this relation to rewrite (1), gives
f_ωφ^sc( ρ)=∫∫γ( ρ′)f_ωφ( ρ′)g_ω(l ρ−pl)dx′dy′  (2)

[0140] These two equations, (1) and (2) are the basis of the imaging algorithm. They must be discretized, then exported to the computer in order to solve practical imaging problems. The purpose of the apparatus and method herein described, is to solve the discretized form of these equations without making any linearizing assumptions (such as used in conventional diffraction tomography), and in real time, with presently available computers.

Discretization of the Acoustic Free-Space Lippmann Schwinger Integral Equation and Green's Function—2D Free Space Case

[0141] Let us for a moment drop the frequency and source position subscripts. The scalar field “γf” is given by

[0142] $\gamma f\left(\begin{array}{c}{x}^{\prime }\\ y^{\prime }\end{array}\right)\equiv \gamma \left(\begin{array}{c}{x}^{\prime }\\ y^{\prime }\end{array}\right)f\left(\begin{array}{c}{x}^{\prime }\\ y^{\prime }\end{array}\right)$

[0143] Discretization of the integral equation is achieved by first decomposing this function into a linear combination of certain (displaced) basis functions

[0144] $\begin{array}{cc}\gamma f\left(\begin{array}{c}{x}^{\prime }\\ y^{\prime }\end{array}\right)\approx \sum _{n^{\prime }=1}^{N_x}\sum _{m^{\prime }=1}^{N_y}{a}_{n^{\prime }m}s\left(\begin{array}{c}{x}^{\prime }-n\delta \\ y^{\prime }-m^{\prime }\delta \end{array}\right)& \left(3\right)\end{array}$
where it has been assumed that the scatterer γ lies within the support [0, N_xδ]x[0, N_yδ]−a rectangular subregion.

[0145] The basis functions S can be arbitrary except that we should have cardinality at the grid nodes:

[0146] $S\left(\begin{array}{c}0\\ 0\end{array}\right)=1$
whereas for all other choices of n′,m′,

[0147] $s\left(\begin{array}{c}{n}^{\prime }\delta \\ m^{\prime }\delta \end{array}\right)=0,n^{\prime }{m}^{\prime }\ne 0.$

[0148] An example of such a function is the 2-D “hat” function. Our algorithm uses the “sinc” function has its basic building block—the Whittaker sinc function which is defined as:

[0149] $\sin c\left(x\right)=\frac{\sin \left(\pi x\right)}{\pi x}$

[0150] The two dimensional basis functions are defined by the tensor product:

[0151] $\begin{array}{cc}s\left(\begin{array}{c}{x}^{\prime }-n^{\prime }\delta \\ y^{\prime }-m^{\prime }\delta \end{array}\right)=\sin c\left(\left(x^{\prime }-n^{\prime }\delta \right)/\delta \right)·\sin c\left(\left(y^{\prime }-m^{\prime }\delta \right)/\delta \right)& \left(4\right)\end{array}$
If the equality in equation (3) is presumed to hold at the grid points x=nδ, y=mδ: the coefficients α_nm can be determined to be precisely the value of the field γf:

[0152] $\begin{array}{cc}\gamma f\left(\begin{array}{c}{x}^{\prime }\\ y^{\prime }\end{array}\right)\equiv \sum _{n^{\prime }=1}^{N_x}\begin{array}{c}{N}_y\\ m^{\prime }=1\end{array}\gamma \left(\begin{array}{c}{n}^{\prime }\delta \\ m^{\prime }\delta \end{array}\right)f\left(\begin{array}{c}{n}^{\prime }\delta \\ m^{\prime }\delta \end{array}\right)s\left(\begin{array}{c}{x}^{\prime }-n^{\prime }\delta \\ y^{\prime }-m^{\prime }\delta \end{array}\right)& \left(5\right)\end{array}$
Now this expression for the product γf can be substituted into (1):

[0153] $\begin{array}{cc}{f}^i\left(\begin{array}{c}x\\ y\end{array}0\right)=f\left(\begin{array}{c}x\\ y\end{array}\right)-\frac{k_{{\omega }^2}}{4i}\sum _{n^{\prime }}\sum _{m^{\prime }}{\gamma }_{n^{\prime }{m}^{\prime }{\mathrm{fn}}^{\prime }{m}^{\prime }}\int H_o^{\left(2\right)}\left(k_{\omega }l\stackrel{\_}{\rho }-{\stackrel{\_}{\rho }}^{\prime }l\right)S\left(\begin{array}{c}{x}^{\prime }-n^{\prime }\delta \\ y^{\prime }-m^{\prime }\delta \end{array}\right)d^2{\rho }^{\prime }& \left(6\right)\end{array}$
In particular this equation holds at

[0154] $\left(\begin{array}{c}x\\ y\end{array}\right)=\left(\begin{array}{c}n\delta \\ m\delta \end{array}\right)$
for which we get:

[0155] $\begin{array}{cc}{f}_{\mathrm{nm}}^i=f_{\mathrm{nm}-}\sum _{n^{\prime }-1}^{N_x}\sum _{m^{\prime }-1}^{N_y}{\gamma }_{n^{\prime }{m}^{\prime }{\mathrm{fn}}^{\prime }{m}^{\prime }\mathrm{fn}-n^{\prime },m-m^{\prime }},& \left(7\right)\end{array}$
n=1, . . . , N_x m=1, . . . , N_y where the 2-D discrete Green's function is defined as:

[0156] $\begin{array}{cc}g\left(n-n^{\prime },m-m^{\prime }\right)\equiv \frac{k^2}{4i}\int H_0^2\left(k\left(\begin{array}{c}n\delta -x^{\prime }\\ m\delta -y^{\prime }\end{array}\right)\right)S\left(\begin{array}{c}{x}^{\prime }-n^{\prime }\delta \\ y^{\prime }-m^{\prime }\delta \end{array}\right)d^2{\rho }^{\prime }& \left(8\right)\end{array}$

[0157] Although it is not clear that the dependence of g on m and n is as so stated, in is indeed the case, since we have by the substitution of the transformation
x′→x′+nδ
y′→y′+m′δ
into this last equation shows explicitly that the discretized Green's function “g”, does in fact depend only on the differences n-n′, and m-m″. Thus the discrete equation (7) has inherited the convolutional form of the continuous equation (1). It is this fact which allows the use of the 2-D FFT to compute the summation in (8) in only order N_XN_Y log₂(N_xN_y) arithmetic operations and is therefore critical to the real-time speed of the algorithm.
Expression of the Discretized Lippmann-Schwinger Equation in Compact Notation

[0158] The discretized Lippmann-Schwinger equation can be rewritten using a compact operator notation. We write (7) as:
f_ωφⁱ=(I−G_ω[γ])f_ωφ  (9)

[0159] where G_ω represents 2-D discrete convolution with the discrete, 2-D Green's function for frequency ω and [γ] denotes the operator consisting of pointwise multiplication by the 2-D array γ. I denotes the identity operator.

[0160] In exactly the same manner, and with the same definitions, we can write the scattered field equation (2) as:
f_ωφ^sc=G_ω[γ]f_ωφ  (10)

[0161] for the scattered field inside of the object, and finally
f_ωφ^meas≡T(f_ωφ^sc)=T(G_ω[γ]f_ωφ)  (11)

[0162] for the measured value of the scattered field at the detectors. The T operator is used to indicate the “truncation” of the calculated scattered field, calculated on the entire convolution range (a rectangle of samples [1,N_x]×[1,N_y] covering γ), onto the detectors which lie outside the support of γ. (see FIGS. 1A/B/C.). In other words T is simply the operation of “picking” out of the rectangular range of the convolution, those points which coincide with receiver locations.

Extension of the Method to Remote Receivers with Arbitrary Characteristics

[0163] In the event that the receivers do not lie within the range of the convolution and/or they are not simple point receivers, we need to modify the measurement equations (11). It is well known that, given source distributions within an enclosed boundary, the scattered field everywhere external to the boundary can be computed from the values of the field on the boundary by an application of Green's theorem with a suitably chosen Green's function, i.e.:

[0164] $f^s\left(\stackrel{\_}{\rho }\right)={\oint }_B{f}^s\left({\stackrel{\_}{\rho }}^{\prime }\right)p\left(\stackrel{\_}{\rho },{\stackrel{\_}{\rho }}^{\prime }\right)d{l}^{\prime },\stackrel{\_}{\rho }$
outside boundary B (12), where P is the Green's function the integral is about the boundary, and d1′ is the differential arclength (in the 3-D case, the integral is on an enclosing surface). Equation (12) allows for the construction of a matrix operator which maps the boundary values of the rectangular support of the convolution (2N_x+2N_y−4 values in the discrete case) to values of the scattered field external to the rectangle. Furthermore, this “propagator matrix” can be generalized to incorporate more complex receiver geometries. For example, suppose that the receiver can be modeled as an integration of the scattered field over some support function, i.e.;

[0165] v_n≡signal from receiver
n=∫∫f^s( ρ)S_n( ρds)  (13)
Where S_n is the support of receiver n. Then from (12):

[0166] $\begin{array}{cc}\begin{array}{c}{v}_n={\oint }_B{f}^s\left({\stackrel{\_}{\rho }}^{\prime }\right)\left\{\int \int S_n\left(\stackrel{\_}{\rho }\right)P\left(\stackrel{\_}{\rho },{\stackrel{\_}{\rho }}^{\prime }\right)ds\right\}d{l}^{\prime },\\ \stackrel{\_}{\rho }\mathrm{outside}\mathrm{boundary}B\\ =\underset{B}{\oint }{f}^s\left({\stackrel{\_}{\rho }}^{\prime }\right)P^{\prime }\left(n,{\stackrel{\_}{\rho }}^{\prime }\right)d{l}^{\prime },P^{\prime }\left(n,P^{\prime }\right)\\ =\int \int S_n\left(\stackrel{\_}{\rho }\right)P\left(\stackrel{\_}{\rho },{\stackrel{\_}{\rho }}^{\prime }\right)ds\end{array}& \left(14\right)\end{array}$
Discretizing the integral gives the matrix equation:

[0167] $\begin{array}{cc}{v}_n=\sum _{i=1}^{N_b}{P}_{{\mathrm{nl}}^{\prime }}{\mathrm{fi}}^s,n=1,\dots ,N_d& \left(15\right)\end{array}$
where N_b=2N_x+2N_y−4 is the number of boundary (border) pixels and N_d is the number of receivers. Equation (15) defines the matrix that we shall henceforth refer to as P or the “propagator matrix” for a given distribution of external receivers. The equation:
v_ωφ^meas=P₁₀₇(G_ω[γ]fωφ)  (16)
includes (11) as a special case for which the receivers are point receivers inside the convolution support. Note that P is a function of frequency, but is not a function of source position. v_ωφ^meas is a vector of dimension N_d.

[0168] The added flexibility of this propagator matrix formulation is particularly advantageous when interfacing our algorithms with real laboratory or field data. Often times the precise radiation patterns of the transducers used will not be known a priori. In this event, the transducers must be characterized by measurements. The results of these measurements can be easily incorporated into the construction of the propagator matrix P allowing the empirically determined transducer model to be accurately incorporated into the inversion. Equations (9) and (16) then provide in compact notation the equations which we wish to solve for the unknown scattering potential, γ. First we consider the forward problem, i.e. the determination of the field f_ωφ for a known object function γ and known incident field, f_ωφⁱ. Then we establish how this forward problem is then incorporated into the solution of the inverse problem, i.e., the determination of γ when the incident fields, and the received signals from a set of receivers are known. Note that the internal field to the object is also—along with the object function γ, an unknown in the inverse problem.

[0169] Since (9) is linear in f_ωφ it can be solved by direct inversion of the linear system:
f=(I−G[γ])⁻¹fⁱ  (17)
In discrete form, this is just a matrix equation (after linearly ordering the 2-D supports of the domain and range into vectors). Since the dimension of the range and domain of the linear system is N_xN_y, the dimension of the linear system is the same giving a total of (N_xN_y)² elements. The arithmetic required to solve the system by standard direct means is thus order (N_xN_y)³. In other words, a doubling of the edge dimension of the problem space will increase the CPU time required by a factor 2⁶=64 times! The arithmetic work required will quickly become intolerable as the size of the problem increases. It is precisely this large growth rate which has convinced many researchers, not to pursue inverse scattering approaches based on integral equations. One could, of course, go to an iterative method. A single iteration of an iterative linear system solver, such as bi-conjugate gradients requires order (N_xN_y)² operations (essentially a matrix-vector product). It can be shown that the growth rate in the number of required iterations for sufficient convergence for the BCG algorithm is order N for this equation. Thus, the overall computational complexity is order N⁵—only one order of N has been saved over direct inversion. Since inverse problems in 2-D generally require order N views, the iteration must be done N times and we are back to order N⁶ computation for BCG.

[0170] The key to overcoming this objection is the convolutional form of (7). If this is exploited by use of the FFT algorithm, the computation needed to perform a matrix-vector product is only order N_xN_y log₂ (N_xN_y). This allows the BCG algorithm to be applied to a single view with order N³ log₂(N) operations and to all views with order N⁴ log₂(N) operations. Due to the slow growth rate of log₂(N), this is essentially a reduction of two orders of N over nonconvolutional methods. Also, this convolutional approach avoids the necessity of storing the order N⁴ computer words needed to contain the linear system in full matrix form. Only order N² words are needed due to the shift invariance of the discrete Green's kernel in (7). It is these two savings, more than any other advance that we have made, that allows us to perform inverse scattering in reasonable time for meaningfully sized problems. The use of BCG over the original CG algorithm also represents a major advance since it converges in essentially the square root of the number of iterations needed by CG. This combination of FFT and CG algorithm was originally developed in [Borup, 1989, Ph.D. dissertation, Univ. of Utah, Salt Lake City], herein incorporated by reference.

The Imaging or Inverse Problem

[0171] In order to solve the imaging problem, we need a set of equations relating the unknown scattering potential, γ, and total fields, f_ωφ, with measurements of the scattered field on a set of external detectors. These detector equations are given in (16). Equation (16), and the internal field equations (9) are the equations which are solved simultaneously to determine γ and the f_ωφ. There are N_xN_y unknowns corresponding to the γ values at each of the grid points, and N_x×N_y×Ω×Φ, unknowns corresponding to the unknown fields, Ω is the number of frequencies, and Φ is the number of source positions (angles). We have improved upon the state of the art by considering the internal field equations to define f_ωφ. for a given γ. Thus the total number of unknowns is reduced to N_xN_y.

[0172] The total number of measurement equations is N_d×Ω×Φ where N_d is the number of detectors. In the general case, where, the sources and detectors do not completely surround the object, the problem of determining γ is “ill-posed”, in a precise mathematical sense, therefore, in order to guarantee a solution, the number of equations N_d×Ω×Φ>N_xN_y, is chosen to be larger than the number of pixel values for over determination. Then the system is solved in the least squares sense. More specifically the solution of (9,16) for γ and the set of fields, f_ωΦ, in the least squares sense is obtained by minimizing the real valued, nonlinear functional:

[0173] $\begin{array}{cc}\underset{\gamma }{\min }\sum _{\omega \varphi }{r_{\omega \varphi }}^2\equiv \begin{array}{c}\min \\ \gamma \end{array}\sum _{\omega \varphi }{v_{\omega \varphi }^{\mathrm{meas}}-P_{\omega }{G}_{\omega }\left[f_{\omega \varphi }\right]\gamma }^2& \left(18\right)\end{array}$
subject to the satisfaction of the total field equations, (9), as constraints. The vector r_ωφ of dimension N_d is referred to as the “residual” for frequency, ω., and angle, φ.

[0174] The methods used to solve the nonlinear optimization problem in our inverse scattering algorithms are, thus far, all “gradient methods” in that second derivative information is not used (as it is in the full Newton and quasi-Newton methods). The principal computation involves the calculation of the gradient vector of (18). A straight forward calculation gives the gradient formula: ∇f(x)=j^H(x)r(x) where the superscript H denotes the Hermitian transpose (complex conjugate transpose of the matrix) and J is the Jacobian of the nonlinear system:

[0175] $\begin{array}{cc}J\left(x\right)=\left[\begin{array}{ccc}\frac{\partial }{\partial x_1}{a}_1& \cdots & \frac{\partial }{\partial x_N}{a}_1\\ ⋮& ⋰& ⋮\\ \frac{\partial }{\partial x_1}{a}_M& \cdots & \frac{\partial }{\partial x_N}{a}_M\end{array}\right].& \left(19\right)\end{array}$

[0176] The simplest gradient algorithm is the Gauss-Newton (GN) iteration. The GN iteration for finding a solution to ∇f=0 is given by:
xⁿ=y−a(x)⁽ⁿ⁾)  (20.1)
δx⁽ⁿ⁾=(J_n^HJ_n)⁻¹J_n^Hr⁽ⁿ⁾  (20.2)
x⁽ⁿ⁺¹⁾=x⁽ⁿ⁾+δx⁽ⁿ⁾  (20.3)
where a is the vector of nonlinear equations. This iteration is well defined assuming that the columns of J_n remain linearly independent. Since (20.2) is equivalent to the quadratic minimization problem:

[0177] $\begin{array}{cc}\begin{array}{c}\min \\ \delta x^{\left(n\right)}\end{array}J_n\delta x^{\left(n\right)}-r^{\left(n\right)}M^2& \left(21\right)\end{array}$
it can be solved by the minimum residual conjugate gradient method (MRCG). This approach also ensures that small singular values in J_n will not amplify noise if care is taken not to overconverge the iteration.

[0178] In the previous patent U.S. Pat. No. 4,662,222 the fields f_ωφ were considered to be independent variables. In this present patent these fields are considered to be dependent variables, with dependence upon γ given implicitly by
f_ωφ⁽ⁿ⁾=(I−G_ω[γ⁽ⁿ⁾])⁻¹f_ωφ^inc

[0179] In order to find the Jacobian expression we must then differentiate the residual vector defined in (18) with respect to γ. The details of this calculation are given in [Borup, 1992.]. The result is:

[0180] $\begin{array}{cc}{J}_{\mathrm{\omega \varphi }}=\frac{\delta r_{\mathrm{\omega \varphi }}}{\mathrm{\delta \gamma }}={P_{\omega }\left(I-\left[\gamma \right]G_{\omega }\right)}^{-1}\left[f_{\mathrm{\omega \varphi }}\right]& \left(22\right)\end{array}$
The final Gauss-Newton algorithm for minimizing (18) subject to equality constraints (9) is:

[0181] GN 1. Select an initial guess, γ⁽⁰⁾

[0182] GN 2 Set n=0

[0183] GN 3 Solve the forward problems using (biconjugate gradients) BCG:

[0184] $f_{\mathrm{\omega \varphi }}^{\left(n\right)}={\left(I-G_{\omega }\left[{\gamma }^{\left(n\right)}\right]\right)}^{-1}{f}_{\mathrm{\omega \varphi }}^{\mathrm{inc}}\begin{array}{c}\omega =1,\dots ,\Omega \\ \varphi =1,\dots ,\Phi \end{array}$

[0185] GN 4 Compute the detector residuals:

[0186] $r_{\mathrm{\omega \varphi }}^{\left(n\right)}=v_{\mathrm{\omega \varphi }}^{\mathrm{meas}}-P_{\omega }{G}_{\omega }\left[f_{\mathrm{\omega \varphi }}^{\left(n\right)}\right]{\gamma }^{\left(n\right),}\begin{array}{c}\omega =1,\dots ,\Omega \\ \varphi =1,\dots ,\Phi \end{array}$

[0187] GN5 If

[0188] $\begin{array}{c}\Sigma \\ \mathrm{\omega \varphi }\end{array}{r_{\mathrm{\omega \varphi }}^{\left(n\right)}}^2<\in ,\mathrm{stop}.$

[0189] GN6 Use MRCG to find the least squares solution to the quadratic minimization problem.

[0190] $\begin{array}{c}\min \\ {\mathrm{\delta \gamma }}^{\left(n\right)}\end{array}\underset{\mathrm{\omega \varphi }}{\Sigma }{r_{\mathrm{\omega \varphi }}^{\left(n\right)}+P_{\omega }{G_{\omega }\left(I-\left[{\gamma }^{\left(n\right)}\right]G_{\omega }\right)}^{-1}\left[f_{\mathrm{\omega \varphi }}^{\left(n\right)}\right]{\mathrm{\delta \gamma }}^{\left(n\right)}}^2$
for δγ⁽ⁿ⁾.

[0191] GN7 Update γ. γ⁽ⁿ⁺¹⁾=γ⁽ⁿ⁾+δγ⁽ⁿ⁾.

[0192] GN8 Set n=n+1, go to GN 3

[0193] The crux of the GN iteration is GN 6 where the overdetermined quadratic minimization problem is solved for the scattering potential correction. This correction is approximated by applying a set of M iterations of the MRCG algorithm. The details of GN 6 are:

[0194] GN.6.1 Initialize the zeroth iterate of MRCG: δγ⁰=0

[0195] GN.6.2 Initialize the MRCG residuals equal to the GN outer loop residuals: r_ωφ⁰=r_ωφ⁽ⁿ⁾

[0196] Note that iterates pertaining to the MRCG iteration are indexed without the ( ) in order to distinguish them from outer GN loop iterates.

[0197] GN.6.3 Compute the gradient of the MRCG quadratic functional:

[0198] $g^0=\underset{\mathrm{\omega \varphi }}{\Sigma }\left[f_{\mathrm{\omega \varphi }}^{*\left(n\right)}\right]{\left(I-G_{\omega }^{*}\left[{\gamma }^{*\left(n\right)}\right]\right)}^{-1}{G}_{\omega }*P_{\omega }^H{r}_{\mathrm{\omega \varphi }}^0$

[0199] GN.6.4 Set the initial search direction equal to minus the gradient: p⁰=−g⁰

[0200] GN.6.5 Set m=0

[0201] GN.6.6 Compute

[0202] $t_{\mathrm{\omega \varphi }}^n=p_{\omega }{G_{\omega }\left(I-\left[{\gamma }^{\left(n\right)}\right]G_{\omega }\right)}^{-1}\left[f_{\mathrm{\omega \varphi }}^{\left(n\right)}\right]p^{m}m\begin{array}{c}\omega =1,\dots ,\Omega \\ \varphi =1,\dots ,\Phi \end{array}$

[0203] GN.6.7 Compute the quadratic step length:

[0204] ${\alpha }^m=\underset{\mathrm{\omega \varphi }}{\Sigma }{t_{\mathrm{\omega \varphi }}^m}^2/{g^m}^2$

[0205] GN.6.8 Update the solution of the quadratic minimization: δγ^m+1=δγ^m+α^mp^m

[0206] GN.6.9 Update the MRCG search residuals: r_ωφ^m+1=r_ωφ^m+α^mt_ωφ^m

[0207] GN.6.10 Compute the gradient of the MRCG quadratic functional:

[0208] $g^{m+1}=\sum _{\omega \varphi }\left[f_{\omega \varphi }^{*\left(n\right)}\right]{\left(I-G_{\omega }^{*}\left[{\gamma }^{*\left(n\right)}\right]\right)}^{-1}{G}_{\omega }^{*}{p}_{\omega }^H{\Gamma }_{\mathrm{\omega \varphi }}^{n+1}$

[0209] GN.6.11 Compute

[0210] ${\beta }^m={g^{m+1}}^2/{g^m}^2$

[0211] GN.6.12 Update the MRCG search direction: p^m+1=−g^m+1+β^mp^m

[0212] GN.6.13 If m=M, go to CN.6.15

[0213] GN.6.14 m=m+1. go to GM.6.8

[0214] GN.6.15 Equate the solution of the quadratic minimization problem with the last iterate of MRCG: δγ⁽ⁿ⁾=δγ^M

[0215] GN.6.16 Return to the GN algorithm at GN 7

[0216] A problem with the algorithm above is the presence of the inverse of the transposed total field equation (I−[γ⁽ⁿ⁾]G_ω)⁻¹ in the computation of the Jacobian and its adjoint in 6.3 and 6.10 Since we do not invert (or even store) these large linear systems, the occurrence of these inverse operators must be dealt with. This problem can however be overcome by computing the action of this operator on a given vector, as needed during the MRCG iteration, by a few iterations of BCG. When performed in this way, the shift invariance of the Green's function is exploited in a maximal way. No large matrices need to be inverted or stored. This is because the Jacobian implementation now consists exclusively of shift invariant kernels (diagonal kernels such as pointwise multiplication by γ or f and the shift invariant kernel composed of convolution with the Green's function) Such shift invariant kernels can be implemented efficiently with the FFT as previously described.

[0217] An even greater increase in numerical efficiency can be obtained in cases for which the approximation:
(I−[γ⁽ⁿ⁾]G_ω)⁻¹≈(I+[γ⁽ⁿ⁾]G_ω)  (23)

[0218] (which is similar to the Born approximation of the forward problem) can be used in the Jacobian. This has been found to be the case for many acoustic scattering problems for biological tissue, and for EM problems for which the contrast in dielectric contrast is small. The other two gradient algorithms that are used to solve the inverse scattering equations are the Fletcher-Reeves and Ribiere-Polak algorithms. For the inverse scattering equations, they are given by the iteration:

[0219] RP 1. Select an initial guess, γ⁽⁰⁾

[0220] RP 2. Solve the forward problems using (biconjugate gradients) BCG

[0221] $f_{\omega \varphi }^{\mathrm{inc}}=\left(I-G_{\omega }\left[{\gamma }^{\left(0\right)}\right]\right)f_{\omega {\varphi }^{\prime }}^{\left(0\right)}\begin{array}{c}\omega =1,\dots ,\Omega \\ \varphi =1,\dots ,\Phi \end{array}$

[0222] for the internal fields, f_ωφ⁽⁰⁾

[0223] RP 3 Computer the detector residuals:

[0224] $r_{\omega \varphi }^{\left(0\right)}=P_{\omega }{G}_{\omega }\left[f_{\mathrm{\omega \varphi }}^{\left(0\right)}\right]{\gamma }^{\left(0\right)}-v_{\mathrm{\omega \varphi }}^{\mathrm{meas}},\begin{array}{c}\omega =1,\dots ,\Omega \\ \varphi =1,\dots ,\Phi \end{array}$

[0225] RP 4 Compute the gradient:

[0226] $r_{\omega \varphi }^{\left(0\right)}=P_{\omega }{G}_{\omega }\left[f_{\mathrm{\omega \varphi }}^{\left(0\right)}\right]{\gamma }^{\left(0\right)}-v_{\mathrm{\omega \varphi }}^{\mathrm{meas}},\begin{array}{c}\omega =1,\dots ,\Omega \\ \varphi =1,\dots ,\Phi \end{array}$

[0227] RP 5 computer the search direction: p⁽⁰⁾=−g⁽⁰⁾

[0228] RP 6 n=0

[0229] RP 7 compute the Jacobian operating on the search direction:

[0230] $t_{\omega \varphi }^{\left(0\right)}=P_{\omega }{G_{\omega }\left(I-\left[{\gamma }^{\left(n\right)}\right]G_{\omega }\right)}^{-1}\left[f_{\mathrm{\omega \varphi }}^{\left(0\right)}\right]p^{\left(0\right).}\begin{array}{c}\omega =1,\dots ,\Omega \\ \varphi =1,\dots ,\Phi \end{array}$

[0231] RP 8 compute the step length (quadratic approximation):

[0232] ${\alpha }_n=-\mathrm{Re}\sum _{\omega \varphi }\left(r_{\mathrm{\omega \varphi }}^{\left(n\right)},t_{\mathrm{\omega \varphi }}^{\left(n\right)}\right)/\sum _{\omega \varphi }{t_{\mathrm{\omega \varphi }}^{\left(n\right)}}^2$

[0233] RP 9 update the solution: γ⁽ⁿ⁺¹⁾=γ⁽ⁿ⁾+α_np⁽ⁿ⁾

[0234] RP 10 Solve the forward problems using (biconjugate gradients) BCG:

[0235] $f_{\mathrm{\omega \varphi }}^{\mathrm{inc}}=\left(I-G_{\omega }\left[{\gamma }^{\left(n+1\right)}\right]\right)f_{{\mathrm{\omega \varphi }}^{\prime }}^{\left(n+1\right)}\begin{array}{c}\omega =1,\dots ,\Omega \\ \varphi =1,\dots ,\Phi \end{array}$

[0236] RP 11 Compute the detector residuals:

[0237] $r_{\mathrm{\omega \varphi }}^{\left(n+1\right)}=P_{\omega }{G}_{\omega }\left[f_{\mathrm{\omega \varphi }}^{\left(n+1\right)}\right]{\gamma }^{\left(n+1\right)}-v_{\mathrm{\omega \varphi }}^{\mathrm{meas}},\begin{array}{c}\omega =1,\dots ,\Omega \\ \varphi =1,\dots ,\Phi \end{array}$

[0238] RP 12 If

[0239] $\sum _{\omega \varphi }r_{\mathrm{\omega \varphi }}^{\left(n\right)}<\in ,\mathrm{stop}$

[0240] RP 13 Computer the gradient:

[0241] $g^{\left(n+1\right)}=\sum _{\omega \varphi }\left[f_{\omega \varphi }{*}^{\left(n+1\right)}\right]{\left(i-G_{\omega }^{*}\left[{\gamma }^{*\left(n+1\right)}\right]\right)}^{-1}{G}_{\omega }^{*}{P}_{\omega }^H{r}_{\mathrm{\omega \varphi }}^{\left(n+1\right)}$

[0242] RP 14 Compute:

[0243] ${\beta }_n=\{\begin{array}{c}{g^{\left(n+1\right)}}^2/{g^{\left(n\right)}}^2\\ g^{\left(n+1\right)},g^{\left(n+1\right)}-g^{\left(n\right)}/{g^{\left(n\right)}}^2,\end{array}$

[0244] Fletcher-Reeves

[0245] Ribere-Pollack

[0246] RP 15 Update the search direction: p⁽ⁿ⁺¹⁾=−g⁽ⁿ⁺¹⁾=β_np⁽ⁿ⁾

[0247] RP 16 n=n+1, go to RP 7

[0248] The distinction between FR and RP lies in the calculation of βin RP 14. It is generally believed that the RP calculation is more rapidly convergent. Our experience indicates that this is true and so we generally use the RP iteration rather than FR.

[0249] Comparison of the linear MRCG iteration and the nonlinear RP iteration reveals that they are very similar. In fact, RP is precisely a nonlinear version of MRCG. Note that the only difference between them lies in the fact that the RP residuals, computed in steps RP.10 and RP.11 involve recomputation of the forward problems while the MRCG residuals are simply recurred in 6.9 (additional, trivial, differences exist in the computation of the α and β's). In other words, the RP iteration updates the actual nonlinear system at each iteration, while the MRCG simply iterates on the quadratic functional (GN linearization). The overall GN-MRCG algorithm contains two loops—the outer linearization loop, and the inner MRCG loop, while the RP algorithm contains only one loop. Sine the RP algorithm updates the forward solutions at each step, it tends to converge faster than GN with respect to total iteration count (number of GN outer iterations times the number, M, of inner loop MRCG iterates). The GN method is, however, generally faster since an MRCG step is faster than an RP step due to the need for forward recomputation in the RP step. The overall codes for the GN-MRCG algorithm and the RP algorithm are so similar that a GN-MRCG code can be converted to an RP code with about 10 lines of modification.

[0250] Before leaving this example it is important to note that the examples given in this patent all assume that the different frequencies, ω, and views, φ are all computed in serial fashion. It is important to note however, that another important link in the real time implementation of our algorithm, is the fact that the different frequencies and different views are independent in both the forward problem and Jacobian calculations, and therefore can be computed independently and in parallel. This should be clear by examining, for example, the computations in GN.3 and GN.6.3. The forward problem calculations, GN.3 are completely independent and in the gradient calculation, GN.6.3, the computations are independent in frequency and view up to the point at which these contributions are summed to give the gradient. The GN and RP algorithms could thus be executed on a multinode machine in parallel with node intercommunication required only 2 to 3 times per step in order to collect sums over frequency and view number and to distribute frequency/view independent variables, such as scattering potential iterates, gradient iterates, etc., to the nodes.

[0251] The examples in the previous patent U.S. Pat. No. 4,662,222 were characterized by the consideration of both the internal fields, and the object function γ to be independent variables of equal importance. The example given above, however, has the distinguishing characteristic that the object function γ is considered to be the sole independent variable, and the internal field resulting from the incident field and γ interaction is considered to be a function of this γ. We may write

[0252] $\begin{array}{cc}{f}_{\theta \omega }^{\mathrm{sc}}\equiv f_{\mathrm{\theta \omega }}^{\mathrm{sc}}\left(\gamma \right)\forall \begin{array}{c}\theta =1,\dots ,\Theta \\ \omega =1,\dots ,\Omega \end{array}& \left(24\right)\end{array}$
for the scattered field.

[0253] The difference between this example, and the previous one is that the background medium, in which the object is assumed to be buried, is assumed to be homogeneous in the previous case, whereas it is here assumed that the background is layered. The residual is defined in the same way as the previous example, with d_θω used to represent the N_d-length vector of measured scattered field at the detector positions. That is, the residual vectors for all ω, θ are defined in the following way.

[0254] $\begin{array}{cc}{r}_{\theta \omega }\equiv f_{\theta \omega }^{\mathrm{sc}}\left(\gamma \right)-d_{\theta \omega }=0\forall \begin{array}{c}\theta =1,\dots ,\Theta \\ \omega =1,\dots ,\Omega \end{array}& \left(25\right)\end{array}$

[0255] The functional F(γ), dependent upon γ is defined in the same way as example 1:
F(γ)≡Σ_θω∥r_θω∥₂²  (26)

[0256] This is the functional, dependent upon the object function γ in a highly nonlinear manner, which must be minimized in order to determine the γ values at each gridpoint in the object space.

[0257] It is again necessary to compute the Jacobian:

[0258] $\begin{array}{cc}\left(\begin{array}{c}\partial f_{\theta \omega }^{\mathrm{sc}}\\ \partial \gamma \end{array}\right)=T_o\left({\left(I-G_{\omega }·\left[\gamma \right]\right)}^{-1}\left(G_{\omega }·\left[f_{\theta \omega }\right]\right)\right)\forall \begin{array}{c}\theta =1,\dots ,\Theta \\ \omega =1,\dots ,\Omega \end{array}& \left(27\right)\end{array}$
where θ again refers to the multiple views and the ω to the multiple frequencies available to us. That is, we again assume that the noise level in the system is zero. We apply the Gauss Newton algorithm to the associated least squares problem and get the same result. This leads to the overdetermined system described above (the notation is identical to the previous example, the difference lies in the actual form of the layered Green's function).

[0259] $\begin{array}{cc}\left[{T_o\left({\left[I-G_k·\left[\gamma \right]\right]}^{-1}{G}_k\right)}_{\otimes }{I}_{\theta x\theta }\right]\left[\begin{array}{c}\left[f_{1k}\right]\\ ⋮\\ \left[f_{\theta k}\right]\end{array}\right]\delta {\gamma }^{\left(n\right)}=-\left[\begin{array}{c}\left[r_{1k}^{\left(n\right)}\right]\\ ⋮\\ \left[r_{\theta k}^{\left(n\right)}\right]\end{array}\right]\forall k\text{:}& \left(28\right)\end{array}$
which must be solved for the γ-update δγ The left hand side of the above formula is given by

[0260] $\begin{array}{cc}\left[{T_o\left({\left[I-G_k·\left[\gamma \right]\right]}^{-1}{G}_k\right)}_{\otimes }{I}_{\mathrm{\theta x\theta }}\right]\equiv \left\{\left[\begin{array}{cc}{{\mathrm{TG}}_k\left(1-\left[{\gamma }^{\left(n\right)}\right]G_k\right)}^{-1}& 0..\\ 0& \dots \\ \dots & \dots \\ \dots & \dots \\ 0& ..0{{\mathrm{TG}}_k\left(I-\left[{\gamma }^{\left(n\right)}\right]G_k\right)}^{-1}\end{array}\right]\right\}& \left(29\right)\end{array}$

[0261] Again, we can use the complex analytic version of the Hestenes overdetermined conjugate gradient algorithm, adapted for least squares problems to iteratively solve this system. This is equivalent to finding the minimum norm solution.

[0262] The formula for the Jacobian in the layered medium situation in the presence of multiple frequencies, is,

[0263] $\begin{array}{cc}\left(\begin{array}{c}\partial f_{\theta \omega }^{\mathrm{sc}}\\ \partial \gamma \end{array}\right)=T_o\left({\left(I-G_{\omega }·\left[\gamma \right]\right)}^{-1}\left(G_{\omega }·\left[f_{\theta \omega }\right]\right)\right)\forall \begin{array}{c}\theta =1,\dots ,\Theta \\ \omega =1,\dots ,\Omega \end{array}& \left(30\right)\end{array}$
where G_ω is the Layered Green's function for the frequency ω. Therefore, in effect, to determine the γ-update, δγ, we merely solve the multiple view problem for each particular frequency, that is, we solve the overdetermined system:

[0264] $\begin{array}{cc}\left[{T_o\left({\left[I-G_k·\left[\gamma \right]\right]}^{-1}{G}_k\right)}_{\otimes }{I}_{\mathrm{\theta x\theta }}\right]\left[\begin{array}{c}\left[f_{1k}\right]\\ ⋮\\ \left[f_{\theta k}\right]\end{array}\right]\delta {\gamma }^{\left(n\right)}=\left[\begin{array}{c}\left[r_{1k}^{\left(n\right)}\right]\\ ⋮\\ \left[r_{\theta k}^{\left(n\right)}\right]\end{array}\right]\forall k=1,\dots ,\Omega & \left(31\right)\end{array}$
which in component form is:

[0265] $\begin{array}{cc}\left[T·\left({\left[I-G_k·\left[\gamma \right]\right]}^{-1}{G}_k\right)\right]\left[f_{\mathrm{jk}}\right]\delta {\gamma }^{\left(n\right)=}-\left[r_{\mathrm{jk}}^{\left(n\right)}\right]\forall \begin{array}{c}k+1,\dots ,\Omega \\ j+1\dots ,\Theta \end{array}& \left(32\right)\end{array}$

[0266] See also the computer code in this patent for the complete description of this procedure.

[0267] For multiple view and multiple frequency inversion then, the inversion scheme in layered media reads:

[0268] STEP 1. Choose an initial guess for the scattering potential, γ^(o).

[0269] STEP 2. Set n=0

[0270] STEP 3. solve the forward problems

[0271] $\left(I-G\left[\gamma \right]\right)f_{\theta \omega }=f_{\theta \omega }^{\mathrm{inc}},\forall \begin{array}{c}\theta =1,\dots ,\Theta \\ \omega =1,\dots ,\Omega \end{array}$
using biconjugate gradients, and use the result to compute the φ-forward maps φ_θω(γ⁽ⁿ⁾) and

[0272] $r_{\theta \omega }^{\left(n\right)}\equiv {\varphi }_{\theta \omega }\left({\gamma }^{\left(n\right)}\right)-d_{\theta \omega }=0\forall \begin{array}{c}\theta =1,\dots ,\Theta \\ \omega =1,\dots ,\Omega \end{array}$

[0273] STEP 4. Determine the L₂ norm of the total residual vector

[0274] ${\mathrm{tr}}^{\left(n\right)}\equiv \begin{array}{c}\left[\begin{array}{c}{r}_{11}^{\left(n\right)}\\ ⋮\\ r_{\theta }^{\left(n\right)}\end{array}\right]\\ ⋰\\ ⋰\\ \left[\begin{array}{c}{r}_{1\Omega }^{\left(n\right)}\\ ⋮\\ r_{\theta \Omega }^{\left(n\right)}\end{array}\right]\end{array}$ $\mathrm{STEP}5.\mathrm{If}{{\mathrm{tr}}^{\left(n\right)}}^2\equiv \sum _{\omega =1}^{\Omega }\sum _{\theta =1}^{\Theta }{r_{\theta \omega }^{\left(n\right)}}^2<\in \mathrm{then}\mathrm{stop}.$

[0275] Step 6. Solve the N_dθΩ by N overdetermined system:

[0276] $T·\left[\begin{array}{ccc}{\left({\left[I-G_1·\left[g\right]\right]}^{-1}{G}_1\right)}_{\otimes }{I}_{\theta ⨯\theta }& \cdots & 0\\ ⋮& ⋰& ⋮\\ 0& \cdots & {\left({\left[I-G_{\Omega }·\left[g\right]\right]}^{-1}{G}_{\Omega }\right)}_{\otimes }{I}_{\theta ⨯\theta }\end{array}\right]$
in the least squares sense for δγ⁽ⁿ⁾, using the conjugate gradient algorithm adapted for least squares problems.

[0277] STEP 7. Update 7 via the formula
γ⁽ⁿ⁺¹⁾=γ⁽ⁿ⁾+δγ⁽ⁿ⁾

[0278] STEP 8. Set n=n+1, go to begin

[0279] The actual values of the angles of incidence θ are chosen to diminish as much as possible the ill conditioning of the problem. Equally spaced values 0<θ<2Π are ideal, but experimental constraints may prohibit such values, in which case the multiple frequencies are critical.

[0280] For biological contrasts (γ<0.15 say) the following approximation, which is essentially identical to the assumption in the previous example, is valid:
[I−G_ω[γ]]⁻¹=[I+G_ω[γ]]  (33)

[0281] This is the same Born-like approximation restricted to the Jacobian calculation only, that is assumed in example 1. This approximation has a much larger radius of convergence than the standard born approximation, but substantially reduces the computational complexity of the problem.

How to Implement the Convolution and Correlation in Layered Medium Imaging

[0282] All of the above calculations have been carried out without using the fact that the operation of matrix multiplication with G is in fact the sum of a convolution and a correlation, which is transformed into a pair of convolutions. Symbolically G^L=G_R+G_V, where G_R is the correlation and G_V is a convolution operator. The actual numerical implementation of the above formulas, therefore, is somewhat different than would naively appear. The implementation of the correlation with convolution is a little more challenging than the convolution in the homogeneous case, in that a change of variables must be used to convert the correlation to a convolution before the Fast Fourier Transform (FFT)s can be applied. A “reflection” operator must be applied to γ before the convolution is applied, this operator is denoted by Z, and is defined by:
(zf)(x,y,z,)≡f(x,y,h=z)  (34)

[0283] $\left[\begin{array}{c}\left[\begin{array}{c}\left[f_{11}\right]\\ ⋮\\ \left[f_{\theta 1}\right]\end{array}\right]\\ ⋰\\ ⋰\\ \left[\begin{array}{c}\left[f_{1\Omega }\right]\\ ⋮\\ \left[f_{\Theta \Omega }\right]\end{array}\right]\end{array}\right]{\delta }_g^{\left(n\right)}=-\left[\begin{array}{c}\left[\begin{array}{c}\left[r_{11}^{\left(n\right)}\right]\\ ⋮\\ \left[r_{\theta 1}^{\left(n\right)}\right]\end{array}\right]\\ ⋰\\ ⋰\\ \left[\begin{array}{c}\left[r_{1\Omega }^{\left(n\right)}\right]\\ ⋮\\ \left[r_{\mathrm{\Theta \Omega }}^{\left(n\right)}\right]\end{array}\right]\end{array}\right]$
The use of this operator is incorporated into the action of the Lippmann-Schwinger equation on the field f That is, the equation for the internal fields, which in example 1 read:
(I−Gγ)f=f^inc(I−Gγ)
Now becomes
(I−G_vγ−G-_RZγ)f=f^inc

[0284] This change has non-trivial ramifications for the construction of the algorithms discussed above. The FFT implementation of this is given below. To prevent the notation from getting out of hand we use [] to indicate that a diagonal matrix is constructed out of an n vector, so that [γ⁽ⁿ⁾]f_φ⁽ⁿ⁾ denotes pointwise multiplication of the vectors γ and f. Also F is used to denote the Fourier transform, and * is used to indicate “convolution”. In this notation, the Lippmann-Schwinger equation becomes, in the symbolic form adopted above,
(I−G_vγ−G_RZγ)f=(f_—G_v*([γ]f)−G_R*Z([γ]f))=(f−F⁻¹{(F(G_v)·F([γ]f))−(F(G_R)·F{Z([γ]f)})})  (35)

[0285] Substantial savings is achieved through the observation that the Fourier Transforms of G_V and G_R need only be done once, then stored, and finally, that the reflection operator, Z, will commute with Fourier transformation in the x and y (horizontal) directions, since it only involves reflection in the z (vertical) direction.

[0286] There are changes also, for example in the computation of f. The biconjugate gradient algorithm requires the use of the adjoint in the solution of the above equation. This adjoint differs from the homogeneous ambient medium case with the introduction of the correlation operator. Also, there are substantial changes in the implementation of the overdetermined conjugate gradient algorithm used to solve (in the least squares sense) the overdetermined system of equations. In this section, as before, we incorporate the k_sc² factor into the G_V and G_R terms. Of course, rather than carry out the inversion of a matrix, we use biconjugate gradients and Fast Fourier Transforms (FFT's), which require the Lippmann Schwinger operator, which we denote by LS, and its adjoint, which are defined by the following formulas:
(LS)f≡(I−F⁻¹└( G_v•F)+( G_R•F∘Z)┘)([γ]f)  (36)
Where G_v≡F(G_v) and G_R≡F(G_R) are the Fourier Transforms of G_V and G_R and need only be calculated once, then stored. The unitarity of F is used below to determine the adjoint (LS)^H.
(LS)^H≡((I−F⁻¹[( G_v·F)+( G_R·F∘Z)])[γ])^H=[ γ](I−[(F⁻¹ G_V^H)+(Z^HF⁻¹ G_R^H)](F⁻¹)^H)=( γ)(I−[(F⁻¹ G_V)+(ZF⁻¹ G_R)]F)  (37)
where we have used the unitarity of Fourier Transformation F:F^H=F⁻¹, and the fact that point wise multiplication is used, and the fact that Z^H=Z. For practical problems it is best to obtain a reasonable first guess using low frequency data, then to build on this using the higher frequency information.
The Actual Construction of the Layered Greens Function

[0287] Up to this point, we have assumed the existence of the layered Green's function, but have not shown explicitly how to construct it. This is an essential step. For simplicity we deal with the scalar case, although the inclusion of shear wave phenomena is conceptually no more difficult than this case. The details are given in [Wiskin, 1991]. The general idea is based upon standard ideas in the literature concerning reflection coefficients in layered media. See for example [Muller, G., 1985, Journal of Geophysics, vol. 58] which is herein incorporated by reference. The novelty is the combination of the reflection coefficients into a bona fide Green's function, and the utilization of this in a forward problem, then more importantly, in the inverse problem solution. The procedure involves decomposing the point response in free space into a continuum of plane waves. These plane waves are multiply reflected in the various layers, accounting for all reverberations via the proper plane wave reflection/transmission coefficients. The resulting plane waves are then resumed (via a Weyl-Sommerfeld type integral) into the proper point response, which in essence, is the desired Green's function in the layered medium. The final result is:

[0288] $\begin{array}{cc}{G}_V\left(x-x^{\prime };\Delta z-\Delta z^{\prime }\right)={\int }_{-\infty }^{\infty }{e}^{i\omega u\left(x-x^{\prime }\right)}C\left(u\right)\left[e^{-i\omega b}{\mathrm{sc}}^{\Delta z-\Delta z^{\prime }}+R^{+}{R}^{-}{e}^{i\omega b}{\mathrm{sc}}^{\Delta z-\Delta z^{\prime }}\right]du,\mathrm{and}& \left(38\right)\\ G_R\left(x-x^{\prime };\Delta z-\Delta z\right)={\int }_{-\infty }^{\infty }{e}^{i\omega u\left(x-x^{\prime }\right)}C\left(u\right)\left[R^{+}{e}^{-i\omega b}{\mathrm{sc}}^{\left(\Delta z-\Delta z^{\prime }\right)}+R^{-}{e}^{i\omega b}{\mathrm{sc}}^{\left(\Delta z-\Delta z^{\prime }\right)}\right]du,\mathrm{where}C\left(u\right)={\left(1-R^{+}{R}^{-}\right)}^{-1}\frac{i\omega u}{b_{\mathrm{sc}}},& \left(39\right)\end{array}$

[0289] R⁻ and R³⁰ are the recursively defined reflectivity coefficients described in Muller's paper, u is the horizontal slowness,

[0290] $u\equiv \frac{\sin \theta }{c}$
Recall that Snell's law guarantees that u will remain constant as a given incident plane wave passes through several layers.

[0291] ω is the frequency

[0292] b_sc is the vertical slowness for the particular layer hosting the scattering potential

[0293] This Green's function for imaging inhomogeneities residing within Ambient Layered Media must be quantized by convolving with the sinc (really Jinc) basis functions described below. This is done analytically in [Wiskin, 1991], and the result is given below.

Discretization of Layered Medium Lippmann-Schwinger Equation

[0294] Unfortunately the basis functions that were used in the free space case (example 1) cannot be used to give the discretization in the layered medium case because of the arbitrary distribution of the layers above and below the layer which contains the scattering potential. The sinc functions may continue to be used in the horizontal direction, but the vertical direction requires the use of a basis function with compact support (i.e. should be nonzero only over a finite region). The sampled (i.e. discrete) Green's operator in the layered case is defined by
G^L(y,k,m)=G^L(jδ,kδ,mδ)  (40)
Where the three dimensional basis functions are of the form:
B(x)=B(x,y,z)=Sj,δ)(x)S_k,δ(y)Λ_m(z),  (41)
(I is the set of integers) and where:

[0295] $x=\left(\begin{array}{c}x\\ y\\ z\end{array}\right)\in R^3{x}_{\mathrm{jkm}}\equiv \left(\begin{array}{c}j\delta \\ k\delta \\ m\delta \end{array}\right)j,k,m\in I$
The basis functions in the horizontal (x-y) plane S_(j,δ) are based upon the Whittaker

[0296] $\sin c\left(x\right)=\frac{\sin \left(\pi x\right)}{\pi x}$
The basis functions in the vertical (z) direction, on the other hand, are given by
Λ_m(Z)=Λ(z−mδ),
where Λ(z) is the “tent” function:

[0297] $\Lambda \left(z\right)=\{\begin{array}{c}\frac{1}{\delta }\left(z+\delta \right)z\in \left[-\delta ,0\right]\\ -\frac{1}{\delta }\left(z-\delta \right)z\in \left[0,\delta \right]\end{array}$
The form of the discretized Lippmann Schwinger equation is

[0298] $\begin{array}{cc}\begin{array}{c}\begin{array}{c}{f}^{\mathrm{inc}}\left(\begin{array}{c}n\delta \\ m\delta \\ l\delta \end{array}\right)\equiv f_{\mathrm{nml}}^{\mathrm{inc}}=f\left(\begin{array}{c}n\delta \\ m\delta \\ l\delta \end{array}\right)-\\ k_{\mathrm{sc}}^2\underset{n^{\prime }{m}^{\prime }{l}^{\prime }}{\Sigma }{G}_{\pi /\delta }^V\left(\begin{array}{c}n-n^{\prime }\\ m-m^{\prime }\\ l-l^{\prime }\end{array}\right){\gamma }_{\int }^L\left(\begin{array}{c}{n}^{\prime }\delta \\ m^{\prime }\delta \\ l^{\prime }\delta \end{array}\right)-\end{array}\\ k_{\mathrm{sc}}^2\underset{n^{\prime }{m}^{\prime }{l}^{\prime }}{\Sigma }{\hat{G}}_{\pi /\delta }^R\left(\begin{array}{c}n-n^{\prime }\\ m-m^{\prime }\\ l-l^{\prime }\end{array}\right){\gamma }_{\int }^L\left(\begin{array}{c}{n}^{\prime }\delta \\ m^{\prime }\delta \\ l^{\prime }\delta \end{array}\right)\end{array}& \left(42\right)\end{array}$
with the vectors G_π/δ^R and G_π/δ^V are the result of first convolving the Green's function with the basis functions, and then sampling at the gridpoints. The superscript R refers to the correlation part whereas the superscript V refers to the convolution part of the Green's function. Both are computed via Fast Fourier Transforms, and are very fast. In the free space case the correlation part is zero. That is, we make the definition for L=V or R—(ie. convolution or correlation):

[0299] $\begin{array}{cc}{G}_{\pi /\delta }^L\left(n,m,l\right)\equiv \int \int {\int }_{R^3}S\left(x^{\prime }\right)S\left(y^{\prime }\right)\Lambda \left(z^{\prime }\right)G^L\left(\begin{array}{c}n\delta -x^{\prime }\\ m\delta -y^{\prime }\\ l\delta -z\end{array}\right)d{x}^{\prime }d{y}^{\prime }d{z}^{\prime },& \left(43\right)\end{array}$
substitution into equation (17) yields the form (42). In matix notation, this equation gives:
f(I−k_sc²G_V*γ_L−k_sc²Ĝ_R*Zγ_I)⁻¹f^inc  (44)
where now * represents discrete convolution in 3D. Now that we have constructed the Integral equations in discrete form, we must derive the dosed form expression for the layered Green's operators. The construction of the closed form expression of the sampled and convolved acoustic layered Green's function G_π/δ^L from the layered Green's function given in equations (38 and 39) above is given by performing the integration in (43) analytically. This process is carried out in [Wiskin, 1991]. The final result is, for the convolutional part,