Diffraction tomography systems and methods.

BACKGROUND OF THE INVENTION

0001This invention relates to systems and methods for reconstructing acoustic and/or electromagnetic properties of two or three-dimensional objects using diffraction tomographic procedures. The systems and methods are useful in the medical diagnostics arts such as ultrasonic and X-ray tomography, as well as in the geophysical prospecting arts, such as well-to-well tomography and subsurface seismic prospecting.

0002In conventional parallel beam transmission computed tomography, an object profile O(x,y) representing, for example, an X-ray attenuation coefficient, is reconstructed from the object's projections<maths id="math0001"><img file="EP0109341A2_D0001.tif" /></maths>where (ξ, η) denote the coordinates in a cartesian coordinate system rotated by the angle θ relative to-the (x,y) system as shown in Figure la. The reconstruction of O(x,y) from P_θ (ξ) (which is called the "reconstruction of the object" or the "reconstruction of the object profile" by those skilled in the art) is made possible by the Projection-slice theorem which states that the one-dimensional Fourier transform of P_θ (ξ) is a slice at angle θ through the two dimensional Fourier transform <img file="EP0109341A2_D0002.tif" />(K_x,K_y) of 0(x,y) as shown in Figure lb. Thus, the projections P_θ (ξ) yield, upon Fourier transformation, an ensemble of slices of the two-dimensional Fourier transform <img file="EP0109341A2_D0003.tif" />(K_x,K_y) of O(x,y) O(x,y) can then readily be reconstructed from this ensemble of slices via its Fourier integral representation expressed in circular polar coordinates. In particular, one finds that<maths id="math0002"><img file="EP0109341A2_D0004.tif" /></maths>where<maths id="math0003"><img file="EP0109341A2_D0005.tif" /></maths>is the one-dimensional Fourier transform of P_θ(ξ). In Equation (2), the K integration is limited to the interval -W to W since, in any application, P_θ(K) will be known only over a finite bandwidth which Ls taken to be -W to _W. Because of this restriction, the resulting reconstruction is a low pass filtered version of 0(x,y).

0003In practice, in place of the Fourier integral representation iven in Equation (2), a technique known in the art as the filtered backprojection algorithm is used to reconstruct O(x,y). Defining the filter function<maths id="math0004"><img file="EP0109341A2_D0006.tif" /></maths>which, when convolved with a low pass filter, is known in the art as the standard X-ray tomography deblurring filter, it can be shown (see, for example, A.C. Kak, "Computerized Tomography with X-ray Emission and Ultrasound Sources"; Proc. IEEE 67, ^pp. 1245-1272 (1974)) that Equation (2) can be written in the form<maths id="math0005"><img file="EP0109341A2_D0007.tif" /></maths>where Q_θ (t) , called "filtered projections", are given by<maths id="math0006"><img file="EP0109341A2_D0008.tif" /></maths>The process of reconstructing O(x,y) according to the prior art thus consists of first filtering the projections with the filter function h(t) and then backprojecting the filtered projections.

0004The process of backprojection used in the prior art consists of assigning any given pixel value (x,y) within the image array the value of the filtered projection Q _θ (t) at the point<maths id="math0007"><img file="EP0109341A2_D0009.tif" /></maths>Backprojection at any given view angle e thus results in a partial image consisting of parallel straight lines as defined by Equation (7) with the grey level of each line being assigned the value Q_a(t). Following the process of backprojection at a fixed view angle, the entire process (filtering followed by backprojection) is repeated at different angles and the resulting partial image arrays are added pixel by pixel onto the previously existing partial image (hence, the integral over de in Equation 5). In this manner, the profile of the object is reconstructed. In practice, of course, the filtered projection is only known for a finite number of angles so that a discrete approximation (e.g. a summation) to the backprojection integral Equation (5) is employed.

0005The success of transmission computed tomography clearly depends on the availability of the projections P_θ(ξ). In X-ray tomography, these projections are available in the form of measurements of the log amplitude [1/2 lo^g(intensity)] of the electromagnetic field produced when a plane, monochromatic X-ray beam propagating along the n axis is incident to the object as is illustrated in Figure 2a. At X-ray wavelengths (⁼lA) the incident plane wave suffers very little scattering (refraction or diffraction) in passing through the object. Absorbtion does occur, however, so that the negative of the log amplitude of the field along the line η=1, is simply the integrated value of the X-ray absorbtion coefficient u(x,y) along the travel path (i.e. along the straight line ξ=constant). The object profile O(x,y) in this case is then µ(x,y) and the projections P_θ (ξ) are simply the negative of the log amplitude of the electromagnetic field as measured along the line n=1,.

0006In ultrasonic tomography, the wavelength is considerably greater (Zlmm.) than that of X-rays so that an acoustic wave experiences scattering, in the form of refraction and diffraction, in the process of propagating through an object. The log amplitude of the acoustic field along the line n=l_o resulting from an insonifying plane wave propagating along the n axis is then not simply the projection of the acoustic absorbtion coefficient along the line ξ=constant as is the case in X-ray tomography. This breakdown in the relationship between measured log amplitude and projections of the absorbtion coefficient is one of the major reasons why conventional computed tomography is of limited use in ultrasonic applications.

0007Although it is not possible, in general, to obtain an exact expression for the acoustic field generated by the interaction of an insonifying plane wave with an acoustic (diffracting) object O(x,y), expressions within the so-called first Born and Rytov approximations are available.

0008The Born and Rytov approximations are essentially "weak scattering" approximations and, as such, will be valid for cases where the deviations of the object profile O(x,y) from zero are small. Which approximation, Born or Rytov, is employed depends upon the specific application, although the Rytov approximation is generally thought to be superior in medical ultrasound tomography (See, for example, M. Kaveh, M. Soumekh and R. K. Mueller, "A Comparison of Born and Rytov Approximations in Acoustic Tomography", Acoustical Imaging, Vol. II Plenum Press New York, (1981).) A more complete treatment of the Born and Rytov approximations is found in A. J. Devaney, "Inverse Scattering Theory Within the Rytov Approximation," Optics Letters 6, pp. 374-376 (1981). Within these approximations, it has been shown by the inventor in "A New Approach to Emission and Transmission CT" 1980 IEEE Ultrasonics Symposium Proceedings, pp. 979-983 (1980); "Inverse Scattering Theory Within the Rytov Approximation" Optics Letters 6, pp. 374-376 (1981); and by R.K. Mueller, et al., "Reconstructive Tomography and Application to Ultrasonics", Proc. IEEE 67, pp. 567-587 (1979), that for an incident monochromatic plane wave propagating along the n axis, the profile O(x,y) is related to the processed signals D_θ0(ξ,ω) (hereinafter referred to as "data") taken along the line n=l, via the equation<maths id="math0008"><img file="EP0109341A2_D0010.tif" /></maths>where the variable κ can assume all values in the range[-k,k]. Here θ, - θ ÷ π/2 is the angle that the wave vector of the insonifying plane wave makes with the x axis, ω is the angular frequency, and k=2π/λ is the wavenumber of the insonifying acoustic field. The data D_θ0 (ξ,ω), are shown in the case of the Born approximation, to be simply related to the complex amplitude of the scattered field portion of the acoustic field evaluated along the line η=l₀ and, in the case of the Rytov approximation, to the difference in the complex phase of the total field and incident field evaluated along this same line.

0009For discussion regarding diffraction tomography, it is convenient to introduce two unit vectors;<maths id="math0009"><img file="EP0109341A2_D0011.tif" /></maths><maths id="math0010"><img file="EP0109341A2_D0012.tif" /></maths>where ξ and η are unit vectors along the ξ, _naxes. The unit vector s, is simply the unit propagation vector of the insonifying plane wave. For values of lκl<k, Equation (8) can be expressed in terms of s and s, in the form<maths id="math0011"><img file="EP0109341A2_D0013.tif" /></maths>From Equation (10) it may be concluded that the one-dimensional Fourier transform <img file="EP0109341A2_D0014.tif" />_θ0 (κ,ω) of D(ξ, ω) is equal to the two-dimensional Fourier transform O(K_x,K_y) of the profile O(x,y) evaluated over the locus of points defined by<maths id="math0012"><img file="EP0109341A2_D0015.tif" /></maths>where the unit vector s assumes all values for which s.s.>0. Equation 10 may thus be seen as the diffraction tomography equivalent of the projection slice theorem. The K values satisfying Equation (11) for fixed so with s·s.>0, are seen to lie on a semicircle, centered at -ks. and having a radius equal to k. By changing the direction of propagation of the insonifying plane wave, it is thus possible to sample O(K_x,K_y) over an ensemble of semicircular arcs such as illustrated in Figure 2b. The ultrasonic tomographic reconstruction problem then consists of estimating O(x,y) from specification of O(K_x,K_y) over this ensemble of arcs.

0010The problem of reconstructing an object profile from specification of its Fourier transform over an ensemble of circular arcs such as shown in Figure 2b is an old one, with proposed solutions having first occurred in X-ray crystallography, H. Lipson and W. Cochran, The Determination of Crvstal Structures (Cornell University Press, Ithica, New York, 1966), and later in inverse scattering applications using both X-rays , B.K. Vainshtein, Diffraction of X-ravs by Chain Molecules (John Wiley, New York, 1977), and visible light, E. Wolf, "Three-dimensional Structure Determination of Semi-Transported Objects from Holographic Data", Optics Comm., pp. 153-156 (1969). More recently, these inverse scattering solutions have arisen in ultrasonics, (see e.g. S.J. Norton and M. Linzer, "Ultrasonic Reflectivity Imaging in Three Dimensions: Exact Inverse Scattering Solutions for Phase, Cylindrical and Spherical Apertures," IEEE Trans Biomed. En^g. BME-28, pp. 202-220 (1981)) and also in ultrasonic tomography which has become known as diffraction tomography to distinguish it from conventional tomography where diffraction effects are not taken into account. In all of these applications, the reconstruction of the profile O(x,y) is made possible by the fact that by varying s_o over the unit circle, one can determine O(K_x,K_y) at any point lying within a circle, centered at the origin and having a radius equal to √2k. A low pass filtered version O_LP(x,y) of 0(x,y) can then be obtained by means of the Fourier integral representation<maths id="math0013"><img file="EP0109341A2_D0016.tif" /></maths>where vector notation is used such that r=(x,y), K=(K ,K_y) and d²K stands for the differential surface element in K .

0011In any application, the profile 0(r) will vanish beyond some finite range so that a Fourier series expansion for O_LP(r) can be used in place of the Fourier integral (12). According to the solutions proposed in the art, in the Fourier series expansion, the transform O(K) must then be known for K values lying on a regular, square sampling grid in K space. To obtain the requisite sample values therefore requires that either one interpolate O(K _x,K_y) from the sample values given on arcs over which this quantity is known (interpolation method) or alternatively, that one select propagation directions s, such that at least one arc falls on every equired sample point in the array (direct sampling method). Discussion and examples of the interpolation method can be found in W.H. Carter, "Computerized Reconstruction of Scattering Objects from Hölograms", J. Opt. Soc Am 60., pp. 306-314 (1976), while the direct sampling method procedure is followed in e.g. A. F. Fercher et al., "Image Formation by Inversion of Scattered Field Data: Experiments and Computational Simulation," Applied Optics 18, pp. 2427-2431 (1979).

0012Both the interpolation and direct sampling methods for obtaining sample values of O(K) over a regular square sampling grid have serious drawbacks. Interpolation of 0(K) from sample values given on semicircular arcs suffers from the same limitations as does interpolation from slices in K space. Reconstructions based on two-dimensional Fourier Transforms from interpolated sample values in conventional tomography are well known to be inferior in accuracy to those obtained via filtered backprojection. Direct sampling methods, on the other hand, are undesirable because they gene,rally require an extremely large number of s, values (i.e., many experiments). In addition, both procedures place severe if not insurmountable demands on the experimentalist in that they require extreme accuracy in registering data collected from repeated tomographic procedures.

0013The present invention is a great improvement over the prior art is that a method and system are provided which allow the low pass filtered object profile O_LP(r) to be reconstructed directly from processed signals D_θ2 (ζ,ω) without the necessity of first determining <img file="EP0109341A2_D0017.tif" />(K) over a regular square sampling grid in K space as was required in the prior art. The provided methods and systems employ the co-called filtered backpropagation algorithm (so named by the invention) which is closely analogous to the filtered backprojection algorithm of conventional transmission tomography. The two techniques differ in one important respect: in the backprojection process the filtered projections Q_θ (t) are continued back through the object space along parallel straight line paths; whereas in the process of backpropagation, the filtered phases are continued back in the manner in which they are diffracted. In other wo'rds, in the backprojection process, the points along a wave of energy may each be considered to follow straight parallel paths through the object. In the backpropagation process, the points are each considered to undergo scattering throughout the object. For example, the Rytov and Born approximations permit an analysis of the problem by suggesting that only first scatterin^gs be taken into account, and that secondary scatterin^gs may be ignored as insignificant. Thus, in backpropagation, quantities Q(t), related to the data D_θ0(ζ,ω) through the filtering operation of equation (6), are mapped onto the object space via an integral transform that is the inverse of the transform that governs phase propagation within the Rytov approximation. It is natural then to call this process backpropagation since it corresponds mathematically to the inverse of the usual forward propagation process. The relationship which is employed and embodied in the filtered backpropagation technique has been conveniently denoted by the inventor as the filtered backpropagation algorithm to emphasize its formal similarity to the filtered backprojection algorithm of conventional tomography.

0014Thus, for purposes of this application and the claims herein, the term "backpropagation" shall be defined to mean that operation that is the inverse or approximate inverse of a forward propagation process. The term "filtered backpropagation technique" shall be defined to describe any diffraction tomographic technique for the partial or complete reconstruction of an object where a filtered real or complex amplitude and/or filtered real or complex phase of a wave is backpropagated into the object space; i.e., is propagated back into object space according to the inverse or approximate inverse of the way in which the wave was originally diffracted. The filtered backpropagation technique is usually implemented in the form of a convolution of filters. For purposes of brevity, such an implementation will identically be called the filtered backpropagation technique. A "backpropagation filter" shall be defined to describe that filter of the filtered backpropagation technique which accounts for diffraction in the backpropagation of the phase; e.g. in ultrasound diffraction tomography, the filter of the filtered backpropagation technique which is not the standard x-ray tomographic filter (or a variation thereon). A "filtered backpropagation operation" shall be defined as any procedure which employs the filtered backpropagation technique.

SUMMARY

0015In accordance with the present invention, a system for reconstruction of an object is comprised of means for obtaining signals which are a function of the phase and amplitude of a diffracting propagating wave which has passed through an object, and filter means for converting the obtained signals by means of a filtered backpropagation technique into an array representing the partial reconstruction of the investigated object. More particularly, the signals required to practice the invention may be obtained from a tape of analog or digital signals gained from one or more tomographic procedures, or from other suitable means.

0016A diffraction tomographic system is comprised of a continuous or pulsed source of wave energy, a detecting system for the measurement of the amplitude and/or phase of the diffracted field resulting from the wave energy interacting with a two or three-dimensional obstacle (object), and a filtering system for processing the measured data according to a filtered backpropagation technique, to thereby generate arrays which are used to reconstruct properties of the object. The wave energy provided may be in the form of plane, spherical or cylindrical waves. The detecting system can be comprised of a single detecting element which is scanned over a detecting surface or, alternatively, can be comprised of an array of elements distributed over the detecting surface.

0017The filtering system disclosed may include a subsystem for preprocessing the received waveforms, a subsystem for implementing a filtered back^propagation technique and thereby performing a partial reconstruction of the object from limited data, a subsystem for integrating a number of partial reconstructions so as to obtain as complete a reconstruction as possible from all available data, and a subsystem for displaying the reconstruction.

0018One set of preferred embodiments of diffraction tomography systems and methods disclosed herein applies to ultrasound transmission computed tomography. An ultrasound source directs a plane wave of ultrasonic energy into the body being examined. The body is interposed between the ultrasound source and the detecting system. The detecting system records the pressure field produced over a detecting plane as a function of time. For continuous sine wave sources, each detector element receives a simple sine wave of specific amplitude and phase. For pulsed wave sources, each detector receives a pulsed wave which will, in general, differ in amplitude and shape from the insonifying pulse. The sources, the object and the detector system are configured so as to allow either the object being insonified to rotate about a specified axis, or alternatively, to allow the ultrasound source and detector surface to rotate about such a specified axis denoted as the axis of rotation of the system.

0019For objects whose acoustic properties do not vary appreciably along the axis of rotation of the system ("two-dimensional" objects), the detector system may comprise a linear array of detector elements aligned along a line lying in the plane of rotation of the system (the measurement line). Alternatively, a single detector element can be scanned along the measurement line. For a two-dimensional object, the detected waveforms are used to reconstruct the acoustic velocity and/or attentuation profiles of the object preferably as follows: The analog electrical outputs from the detecting system are converted to digital signals, Fourier transformed in time, normalized, and the normalized signals are processed in such a way as to obtain their unwrapped complex phase. The unwrapped complex phase is input to the subsystem implementing the filtered backpropagation technique. In this subsystem, the unwrapped phase is filtered with a filter function which, in the frequency domain is the product of the standard x-ray deblurring filter with the backpropagation filter. The backpropagation filter is non-stationary, being dependent on the perpendicular distance of the coordinate point at which a reconstruction is performed from the measurement line. For any given angular orientation between the tomographic system and object (view angle), one filtered backpropagation operation (implementation of the filtered backpropagatio technicue) yields a two-dimensional, image array which represents the partial reconstuction of the object as determined at the particular view angle employed. This image array is then interpolated onto a master image array whose orientation is fixed in space and independent of the view angle. By repeating the entire sequence of steps for a number of view angles and adding the partial reconstructions so obtained in each step to the master array, a complete rendition of the (two-dimensional) object profile is obtained and may be displayed on a computer-graphics display system.

0020Another group of ultrasound transmission computed diffraction tomography systems and methods disclosed herein applies to three-dimensional objects whose acoustic properties vary significantly along the axis of rotation of the system. In this configuration, 3 two-dimensional array of detector elements may be placed along a plane lying perpendicular to the plane of rotation of the system (the measurement plane). Alternatively, one or more detector elements can be scanned over the measurement plane. In the preferred embodiment, the analog electrical outputs from the detecting system are converted to digital signals, Fourier transformed in time and normalized, and the normalized signals are then processed in such a way so as to obtain a one-dimensional projection of the unwrapped complex phase of the transformed signals onto the line formed by the intersection of the measurement plane with the plane of rotation of the system.

0021The projected unwrapped complex phase is input to the subsystem implementing the filtered backpropagation technique which is identical to that used for two-dimensional objects described above. The partial image obtained from the filtered backpropagation operation is then interpolated onto a master image array whose orientation is fixed in space and independent of view angle. By repeating the entire sequence of steps for a number of view angles and adding the partial reconstruction so obtained in each step to the master image array, a complete rendition of the planar projection of the object onto the plane of rotation of the system is obtained. By changing the plane of rotation of the system (the planar orientation), any planar projection of the object can be obtained. The complete three-dimensional object can itself be reconstructed from a sufficient number of planar projections so obtained, if desired. Again, the partial reconstruction or the entire three-dimensional reconstruction may be displayed as desired.

0022The present invention also applies to X-ray transmission computed tomography and especially to what is known in the art as "soft" X-rays which have longer wavelengths than standard X-rays. In this embodiment, a quasi-monochromatic X-ray source is arranged opposite an X-ray detecting system. An object is interposed between the source and detecting system and the-system is configured to allow relative rotation between the object and source/detector(_s). The detected X-ray amplitude is employed to reconstruct a planar cross section of the X-ray attenuation profile of the object being examined.

0023As described above, the invention is embodied in numerous ways. Thus, for ultrasound computed tomography, two or three-dimensional image reconstructions can be obtained. For X-ray transmission computed tomography, a two-dimensional image may be obtained. Moreover, continuous or pulsed sources, and different detector systems may be used.

0024Another set of'embodiments of the invention is provided by electromagnetic or sonic well-to-well tomography. In these geophysical embodiments, one or more sources of electromagnetic or sonic waveforms are placed in a borehole. One or more detectors of electromagnetic or sonic waveforms are placed in a second borehole. The waveforms received in the second borehole are then used to reconstruct the velocity and/or attenuation profile of the formation between the two boreholes. In the borehole to borehole (well to well) situation, only a-two-dimensional reconstruction is possible due to the impossibility of rotating either the object (formation) or the sources and detectors.

0025Yet another set of embodiments of the present invention involves subsurface electromagnetic or seismic exploration where electromagnetic or sonic sources are placed on the earth's surface in the immediate vicinity of a borehole and one or more detectors are placed in the borehole. (A second configuration has the detectors on the surface and the sources in the borehole.) The waveforms received by the detectors are used to reconstruct the velocity and/or attenuation profile of the formation existing between the sources and the detectors. Again, due to rotation limitations, only a two-dimensional image may be reconstructed according to the procedures detailed above.

0026When spherical waves, which are generated by point sources, are used in well-to-well tomography and subsurface seismic exploration, the analog electrical signals produced by the detecting system are converted to digital signals and processed in such a way as to obtain the signal that would have been detected if the transmitted wave were a continuous, cylindrical waveform having a linear phase shift along the transmitter axis. This later signal is then processed so as to obtain its unwrapped complex phase.

0027The unwrapped complex phase is input to the subsystem implementing the filtered backpropagation technique. This subsystem is similar to that used for ultrasound transmission computed tomography using plane waves described above, but employs a variation of the ultrasound backpropagation filter. Additionally, the standard X-ray deblurring filter is altered substantially. Following the filtered backpropa^gation operation, the partial image array is transferred to a master image array which may be displayed by suitable means and the entire sequence of steps is repeated for a different choice of linear phase shift. The final two-dimensional image obtained is a rendition of the electromagnetic or acoustic velocity and attenuation profile of the plane section of the formation lying between the two boreholes, or the borehole and the surface, whichever the case may be.

0028It is thus an object of the instant invention to provide methods and systems for reconstructing acoustic and/or electromagnetic properties of two and three-dimensional objects.

0029It is another object of the invention to provide diffraction tomography methods and systems for two and three-dimensional object reconstruction for use in ultrasonic medical diagnosis.

0030It is yet another object of the invention to provide diffraction tomography methods and systems for sonic and/or electromagnetic reconstruction of subsurface earth formations.

BRIEF DESCRIPTION OF THE DRAWINGS

0031The present invention may be better understood and its numerous objects and advantages will become apparent to those skilled in the art by reference to the detailed description of the preferred embodiments and by reference to the representative drawings wherein:

0032<ul id="ul0001" list-style="none"><li>Figure la illustrates the projection obtained in a cartesian coordinate system when conventional parallel beam transmission computed tomography is applied to an object.</li><li>Figure Ib illustrates a slice through the two-dimensional Fourier transform of an object obtained via the projection slice theorem in conventional parallel beam transmission computed tomography.</li><li>Figure 2a illustrates the projection of an object in the prior art X-ray tomography.</li><li>Figure 2b illustrates the diffraction tomography equivalent of the conventional projection slice theorem where the two-dimensional Fourier transform of an object must be evaluated on an ensemble of circular arcs.</li><li>Figure 3 is a partially schemetic and partially block diagram of an ultrasound transmission computed tomographic system of the present invention.</li><li>Figure 4 is a block diagram of the preprocessing performed on detected signals in the two-dimensional object system embodiment.</li><li>Figure 5 is a block diagram of the preprocessing performed on detected signals in an alternative preferred embodiment.</li><li>Figure 6 is a block diagram of the preprocessing performed on detected signals in the three-dimensional object system embodiment.</li><li>Figure 7 is a block diagram of the filtered backpropagation technique.</li><li>Figure 8a-is an illustration of an array obtained by the ultrasound preprocessing of received signals from an ultrasound transmission computed tomographic system.</li><li>Figure 8b is an illustration of the results obtained after the arrays of Figure 8a are filtered by the filtered backpropagation technique of Figure 7.</li><li>Figure 9 is an illustration of the rotation-interpolation operations used to provide a more complete reconstruction.</li><li>Figure 10 illustrates the rotation-interpolation operations in block diagram form.</li><li>Figure 11 is a block diagram illustrating the method inventions.</li><li>Figure 12A is a partially schematic, partially block diagram of a geophysical tomographic system depicting a source-detector and signal preprocessing arrangement.</li><li>Figure 12b is an illustration of an alternate preferred geophysical tomographic system embodiment.</li><li>Figure 12c illustrates the borehole to borehole geophysical tomographic system embodiment.</li><li>Figure 13 illustrates the formation of a cylindrical wave of a determined phase from a series of spherical waves.</li><li>Figure 14 is a block diagram of the preprocessing and filtering performed on the signals in a geophysical tomographic system embodiment.</li><li>Figure 15 illustrates the geometry and the associated coordinate system for the geophysical tomography embodiments.</li><li>Figure 16 illustrates the regions in K space over which the object's Fourier transform is known for geophysical and ultrasonic tomographic procedures.</li></ul>

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

ULTRASOUND DIFFRACTION TOMOGRAPHY

0033As seen in Figure 3, an object 10 interposed in a fluid bath 12 is investigated by the ultrasound transmission computed tomographic system of the present invention. An ultrasound source 14 and a detector plane 16 are located on opposite sides of the object. The source 14 is such that the acoustic wave entering the object 10 is approximately a plane wave. This can be accomplished by appropriately positioning a small disk-shaped single transducer element, or by means of a suitably phased planar array of transducer elements, or by any other suitable means. The ultrasound signal leaving the source can be of either a continuous wave or pulsed type. The center frequency of the source is optimally such that the wave length of the pressure field induced in the fluid bath at this frequency is significantly (i.e., on the order of one magnitude or more) smaller than the dimensions of the object being investigated. The detector plane 16 is preferably aligned parallel to the plane wavefront cf the insonifying plane wave and contains one or more detector elements 18.- In a typical system the transducer elements used for the source 14 and detectors 18 will be piezo-electric crystals.

0034The system is configured so as to allow rotation between the object, and the detector and source elements of the system, about an axis 0 passing through the object and parallel to the detector plane. For purposes of identification, the plane which is perpendicular to the axis of rotation and which bisects the source, is denoted as the plane of rotation of the system. Thus, in Figure 3, the plane of rotation is the plane of the paper. Likewise, the angle 20 formed between the direction of propagation of the insonifying plane wave 22, and an arbitrary reference line 24 fixed relative to the object and passing through the axis of rotation in the plane of rotation, is denoted as the view angle θ.

0035After the insonifying wave propagates through the object, the resulting diffracted field is detected by one or more detector elements 18 on the detector plane 16. The detectors 18 are preferably transducers which output electrical signals e _m(t) (the voltage at the m'th detector over time). These signals are input into computer and graphics systems 25 which implement a filtered backpropagation technique, generate a partial reconstruction of the object, and display the results.

0036One configuration of the systems disclosed herein is represented in Figure 4, where the acoustic properties of an object do not vary appreciably along the axis of rotation of the system. In these two-dimensional objects, the axis of symmetry of the object is the axis of rotation of both the object and the source-detector system. The detecting system of the preferred embodiment in this case consists of a linear array of M elements evenly spaced over the measurement line -- this line being formed by the. intersection of the detector plane 16 with the plane of rotation of the system. The detector elements 18a are separated by a distance approximately equal to the fluid bath wavelength λ_ocorresponding to the center frequency f₀ of the insonifying plane wave.

0037As shown in Figure 4, the output signals e_m(t) from transducers 18a are preprocessed by inputting them, if required or desired, into analog electronic circuitry and a buffer storage system 28a for prefiltering, calibration and temporarily storing the analog electrical signals in preparation for analog to digital conversion. A/D converter 30a converts the stored analog signals into digital format. In one embodiment of the present invention, the digitized electrical signals are then Fourier transformed in time and normalized by suitable means 32a according to the equation<maths id="math0014"><img file="EP0109341A2_D0018.tif" /></maths>where U₀ is the amplitude of the insonifying wave at its center frequency f_o, and k is the wavenumber in the fluid bath at this frequency. The frequency (transform variable) has been set equal to the center frequency f₀ of the insonifying waveform as those skilled in the art will recognize, because the insonifying wave output by source 14 may be a pulse which is comprised of many sine waves at different amplitudes, phases and frequencies. By choosing the center frequency, a constant frequency is employed which will provide the amplitude and phase of the signal at an optimal signal to noise ratio.

0038After transformation and normalization, electrical output signals E_m are further preprocessed by converting them to a complex phase φ_m by suitable means 34a according to the following operations:<maths id="math0015"><img file="EP0109341A2_D0019.tif" /></maths><maths id="math0016"><img file="EP0109341A2_D0020.tif" /></maths>where the principle value of the arc tangent is computed. For each value of m, the quantity φ_m represents the principle value of the complex phase of the pressure field at frequency f_o at the spatial position of the m'th detector element. The usefulness of complex phase converter 34a may be understood by realizing that when an energy wave goes through a diffracting body, the body has two effects: It attenuates the waves and thus changes the amplitude; and it changes the phase of the wave. Thus, the complex phase converter 34a manipulates the incoming data so that the amplitude information is contained in Real and the phase shift is found in Imag φ_m.

0039After the electrical signals are converted into complex phases φ_m, these complex phases are phase unwrapped at 36a and the unwrapped phase array p_m is used to perform a partial reconstruction of the object. Phase unwrapping is desirable because complex phase converter 34a treats input values such as 0, 2π and 4₁ alike, and therefore the output from converter 34a does not necessarily reflect the proper signal. Phase unwraaper 36a, however, as seen in Figure 4 by its additional inputs, utilizes the knowledge that the shift between adjacent points cannot be too large (i.e. >2π), and therefore corrects the phase array. Any standard phase unwrapper may be employed for these purposes.

0040Another manner of preprocessing signals e_m(t) concerns the cases where the pulse duration of the insonifying plane wave is very short. In this situation the electrical signals e_m(t) from the detector array are processed and stored at 28a, A/D converted at 30a and then input into a preprocessing system which may utilize software to estimate the total delay time ζ_m between the pulse received at the detector and the transmitted pulse. The peak amplitude of the received pulse A_m is also estimated. A block diagram of the preprocessing steps in this embodiment is shown in Figure 5.

0041After and A_m are calculated at 29 and 31 respectively, the unwrapped complex phase p_m is approximately determined at 35 via the relationships<maths id="math0017"><img file="EP0109341A2_D0021.tif" /></maths><maths id="math0018"><img file="EP0109341A2_D0022.tif" /></maths>where ω_o= 2πf_o is the center of the insoifying pulse. The above two relationships are only approximate but, nevertheless, provide good approximations in cases where the total pulse duration is very short. In addition, the approximation obtained for Imag µm can be used, if desired, to aid in the phase unwrapping operation 36a required in the embodiment discussed above.

0042Yet another manner of preprocessing signals e_m(t) concerns cases where the Born approximation for the interaction between the insonifying wave and object is employed. For such cases the normalized, Fourier transformed transducer output signals E are processed according to the equations<maths id="math0019"><img file="EP0109341A2_D0023.tif" /></maths>where k is the wavenumber of the insonifying acoustic field and ℓ_o is the perpendicular distance from the center of rotation of the object to the detector array, to yield the array of complex numbers of p_m.

0043A second configuration of the systems disclosed herein applies to objects whose acoustic properties vary appreciably along the axis of rotation of the system (3-D objects). In the preferred embodiment of this configuration, as seen in Figure 6, the detecting system consists of a rectangular NXM array of transducer elements 18 evenly spaced over the detector plane 16. The centers of the detector elements are separated by a distance on the order of the wavelength λ_o corresponding to the center frequency f_o of the insonifying wave. The electrical signals generated by the detector elements are input into the analog electronic circuitry and a buffer storage system 28b which prefilters, calibrates and temporarily stores the analog electrical signals in preparation for the A/D conversion. After A/D conversion at 30b, the digitized electrical signals from the transducer elements are used to calculate the two-dimensional unwrapped phase array p_nm of the pressure field at the detector locations. Again, as in the case of two-dimensional objects, the preprocessing of the digitized signals can be performed in three different preferred manners.

0044In one of the preferred manners, the signals are digitized, Fourier transformed in time and normalised at 32b and the principle valued complex phase φ_nm corresponding to each detector location on the detector plane is calculated at 34b in the same manner as was calculated in the two-dimensional configuration. After the signals are converted into complex phases φ_nm, the two-dimensional array of complex phases are phase unwrapped at 36b to yield the unwrapped phase array P_nm.

0045A second preferred embodiment of the three-dimensional configuration applies to cases where the pulse duration of the insonifying pulse is very short. In this embodiment, the digitized electrical signals e_nm(t) from the detector array are preprocessed so as to yield the delay time ζ_nm and peak amplitude A_nmof the received pulses. These quantities are then processed according to the equations<maths id="math0020"><img file="EP0109341A2_D0024.tif" /></maths><maths id="math0021"><img file="EP0109341A2_D0025.tif" /></maths>to yield approximate values of the unwrapped complex phase array P_nm. nm

0046A third preprocessing method applies when the Born approximation is employed. For this case, the complex array P_nm is calculated according to the equations<maths id="math0022"><img file="EP0109341A2_D0026.tif" /></maths> and the _Pnm are then further processed as described below.

0047Following the calculation of the two-dimensional unwrapped complex phase array P_nm' this array is projected onto the line formed by the intersection of the detector plane with the plane of rotation of the system. The one-dimensional array representing the projected unwrapped complex phases in the three-dimensional system is denoted by p_m and is calculated by summing the unwrapped phase array column fashion according to<maths id="math0023"><img file="EP0109341A2_D0027.tif" /></maths>as seen at 38. This one-dimensional array p_m is used to perform a partial reconstruction of the object. The projection of the two-dimensional array to arrive at the one-dimensional array p_m is required in performing either a reconstruction of a three-dimensional object profile or a reconstruction of a two-dimensional planar slice of an object.

0048Additionally, it should be noted that the one-dimensional array used in partially reconstructing the object is identically denoted pm in both the two-dimensional and three-dimensional configurations. While the notation is identical, as a one-dimensional array is used in the filtered backpropagation technique, the arrays are arrived at differently and contain different information. Those skilled in the art will appreciate that if the three-dimensional embodiment is used on a two-dimensional object to generate array p_m, the result after suitable normalization would be identical to the p_m array generated by the preferred two-dimensional system.

0049After the above-discussed signal preprocessing, the array of complex numbers p_m (m = 1, 2 ..., M) is input into the backpropagation subsystem 50 as shown in Figure 7. Subsystem 50 implements the filtered backpropagation technique. In this embodiment, the quantities p_m (m = 1; ..., M) are Fourier transformed at 52 according to the relationship<maths id="math0024"><img file="EP0109341A2_D0028.tif" /></maths>(where j is a dummy variable) to yield array <img file="EP0109341A2_D0029.tif" />_m(m = 1, ..., M). The M complex numbers <img file="EP0109341A2_D0030.tif" />_m are then input to the filter bank 54 which performs the filtered backpropagation operation (i.e. implements t_i filtered backpropagation technique) and generates the two-dimensional array <img file="EP0109341A2_D0031.tif" />_qm<img file="EP0109341A2_D0032.tif" />m (q = 1, 2,..., L, m = 1, 2,...,M) where L is the total number of depth samples desired for the object reconstruction. Finally, L one-dimensional Fourier transforms are performed to yield the two-dimensional filtered and backpropagated complex phase X_qm according to the relationship<maths id="math0025"><img file="EP0109341A2_D0033.tif" /></maths>where, again, j is used as a dummy variable.

0050In the filter bank, the two-dimensional array Hq_m (q = 1, 2, ... L) (m= 1, 2, ...M) is produced according to the filtered backpropagation technique<maths id="math0026"><img file="EP0109341A2_D0034.tif" /></maths>The quantities h are the standard X-ray tomography smoothing filter, i.e. the Fourier coefficients of a low-pass smoothing filter convolved with an X-ray tomography deblurring filter as discussed above in the Background. The remainder of the filter (eⁱ(Υ_m k)(z_q - l_o)₎ is the preferred backpropagation filter which performs the backpropagation operation. Υ_m are defined as<maths id="math0027"><img file="EP0109341A2_D0035.tif" /></maths>with k being the wavenumber associated with the center frequency of the insonifying plane wave. The quantities z_g (q = 1, 2,...L) are the depth coordinates of lines over which the filtered and backpropagated phase is being evaluated, and l_o is the distance between the axis 0 of the object and the detector plane. Thus, it is seen that the preferred backpropagation filter which is one of the filters of the filtered backpropagation technique, is a non-stationary filter, as it is dependent on object depth z . The backpropagation filter accounts for the diffraction of waves in the object and permits the object profile to be reconstructed directly from the gathered information without having to use interpolation or direct sampling over a square sampling grid.

0051For two-dimensional objects, the X q_m array represents the partial reconstruction of the object profile as determined at the particular view angle employed in obtaining the data. For three-dimensional objects, this array represents the partial reconstruction of the planar projection of the object profile onto the plane of rotation of the system. The q index indicates the depth coordinate while the m index indicates the coordinate position along the line of intersection of the plane of rotation of the system with the detector'plane. Figure 8(a) illustrates the p_m array at the M sample points along the line of intersection for the case of a 45 view angle. For the purposes of the illustration, p_m is plotted as a single real array of sample values. In general, p_m is a complex array so that two plots (one for the real part of the array and one for the imaginary part) are required. In Figure 8(b), a filtered and backpropagated phase array X is shown for selected object depth values. Note, that again, X_gm is plotted as being a single two-dimensional real array whereas, in practice a two-dimensional array would be required.

0052Referring to Figure 8(b), it is seen that the orientation of the ^q and m axes will always be the same, independent of the view angle. However, the partial reconstruction of the object represented by the X_gm array will depend on the view angle. Thus, as seen in Figure 9 it is necessary to rotate the X q_m array about the axis of rotation of the system by the angle -e, and then interpolate the rotated array onto a master image array which is fixed relative to the object in the plane of rotation of the system. For two-dimensional objects, this master image array represents the ultrasound velocity and attenuation profile of the object. For three-dimensional objects, this array represents the projection of the object's three-dimensional velocity and attenuation profile onto the plane of rotation of the system. A block diagram of the rotation-interpolation operation is given in Figure 10 wherein X_gm is rotated by angle -e, at 111, and interpolated onto the master array at 113. In,practice, any interpolation scheme may be employed. The nearest neighbor interpolation technique is well-suited for these purposes.

0053Finally, by-repeating the entire sequence of operations for a new view angle and adding the rotated and interpolated X array pixel by pixel onto the master array, a new and improved partial reconstruction of the object (or its planar projection) will be obtained. By repeating this for a number of view angles, a complete rendition of the object (or its planar projection) will be generated. For three-dimensional objects, the plane of rotation of the system can then be changed and the entire sequence repeated to obtain a projection of the object onto a new plane. In this way, any number of planar projections can be obtained from which a complete three-dimensional rendition of a three-dimensional object can be generated using known 3-D tomographic reconstruction techniques. Of course, the master array and if desired, individual p arrays.may be displayed in both two-dimensional and three-dimensional situations on a screen or as a hard copy as desired.

0054For two-dimensional objects, the master image array obtained from the filtered backpropagation technique is a complex valued two-dimensional array. Its real part is a rendition of the real part of the object profile while its imaginary part is a rendition of the imaginary part of the object profile. The real' part of the object profile carries information primarily about the acoustic velocity profile of the object and may be displayed on a hard copy or video display system by simply displaying the real part of the master image array. The imaginary part of the object profile carries information primarily about the acoustic attenuation profile of the object and may be displayed on a hard copy or video display system by simply displaying the imaginary part of the master image array. Alternatively, the real (velocity) and imaginary (attenuation) profiles may be displayed jointly using a color graphics display system with different colors being used for velocity and attenuation. It should be noted that both the acoustic velocity and acoustic attenuation profiles of tumors are presently important in medical diagnostics of cancer (see, for example, J. F. Greenleaf and R. C. Bah, "Clinical Imaging With Transmissive Ultrasound Computerized Tomography," IEEE Trans. Biomed. Eng., BME--28, 177-185 (1951).)

0055As in the case of two-dimensional objects, the master image array obtained for any given plane of rotation for a three-dimensional object and system is a complex valued two-dimensional array. Its real part is a rendition of the projection of the real part of the three-dimensional object profile onto the plane of the rotation of the system while the imaginary part is the projection of the imaginary part of the three-dimensional object profile onto the plane of rotation of the system. The real and imaginary parts of the two-dimensional master image array may be displayed separately or jointly on hard copy or video graphics display systems as can be done for two-dimensional objects. The real part of the array carries information primarily about the acoustic velocity profile of the three-dimensional object while the imaginary part carries information primarily about the acoustic attenuation profile of the object.

0056It is sometimes desirable to display sections (planar slices) of objects such as biological media rather than planar projections of such media. Sections can be reconstructed from a sufficient number of planar projections by means of the three-dimensional projection slice theorem. In the present invention, the planar projections of the three-dimensional object profiles as generated by the filtered backpropagation technique can be used in conjunction with the three-dimensional projection-slice theorem to yield planar sections of the object profile of the three-dimensional body being examined. The planar sections so obtained can then be displayed on a graphics display system in a way entirely analogous to that employed for two-dimensional object profiles.

0057The planar projections of the object profiles resulting from the filtered backpropagation technique can also be used to obtain an approximate reconstruction of the full three-dimensional object profile. The real and imaginary parts of the three-dimensional profile can then be displayed using known three-dimensional display methods either separately or jointly. Again, as was the case for two-dimensional objects, the real part of the three-dimensional profile carries information primarily about the three-dimensional acoustic velocity profile while the imaginary part carries information primarily about the three-dimensional acoustic attenuation profile.

0058Figure 11 depicts in block diagram form the best mode embodiment of the method for tomographic reconstruction of an object using ultrasound transmission computed tomography. The methods of the invention are integrally connected with the systems described above as will be recognized by those skilled in the art. At step 114, a wave of energy is emitted from a suitable ultrasound source toward the object 115 being examined. The diffracted pressure field produced by the interaction of the wave with the object is detected at step 118 by one or more detectors. As the detectors are typically transducers, the detectors produce an analog signals which are functions of the detected field over time. The produced signals are preprocessed in steps 121, 128, 130, 131, 132, 134, 136, and 138 where they are respectively calibrated, stored, converted into a digital signal, normalized and Fourier transformed, complex phase calculated, phase unwrapped and for three-dimensional objects, projected onto the line formed by the intersection of the detector plane with the plane of rotation of thø system. As described above, where the pulse duration of the insonifying wave is very short, steps 132, 134 and-136 are omitted, and in their place delay times and peak amplitude are measured and processed to determine the complex phases at 145. Also, for cases involving the Born approximation, these complex phases are readily calculated directly from the output of 132. In all three methods, the preprocessed signals are input into a system which implements the filtered backpropagation technique in steps 152, 154, and 156. Step 152 Fourier transforms the preprocessed signal, which is a one-dimensional array, and step 154 performs the filtered backpropagation technique. Step 156 Fourier transforms the two-dimensional filtered and backpropagated complex phase resulting from step 154 and in doing so generates partial reconstruction of the object as determined at the particular view angle employed in obtaining the data.

0059In steps 160, 162, and 164, the two-dimensional array which is the partial reconstruction of the object is rotated, interpolated, and stored on the master array. By repeating the. method described at a'number of view angles and adding the resulting partial reconstructions obtianed at each view angle to the master image array, a complete rendition of a two-dimensional - object profile or the planar projeotion of a three-dimensional object profile onto the plane of rotation of the system is obtained. A complete three-dimensional object profile or a planar slice of the object profile can be reconstructed from a sufficient number of planar projections so obtained, if desired. The obtained arrays, slices, projections of profiles, or profiles may be displayed at 166 using suitable display means.

0060Those skilled in the art will recognize that many permutations and variations may be made in the methods and systems disclosed for the ultrasound computed tomography invention and the disclosed embodiments are not intended to limit the invention in any manner. For example, the detector system may comprise one detector which is scanned over a line for two-dimensional objects or over a plane for three-dimensional objects. Alternatively, the detector system may include a one-dimensional array of detectors for a two-dimensional object or a one-dimensional array that is scanned in a direction perpendicular to the array for a three-dimensional object. A two-dimensional array of detectors may be used for three-dimensional objects. Another possibility would be to use a single large detector which could rrecognize variations in the variations in the diffracted field. Those skilled in the art will recognize other possible detector schemes.

0061While it is preferable that the detectors be aligned parallel to the plane wavefront of the insonifying plane wave, those skilled in the art will recognize that if desired, the detectors can be placed at an angle and the signals detected can be additionally corrected for the deviation. For example, arrangements such as the hereinafter-discussed geophysical tomography embodiments where the sources and detectors are non-parallel and stationary but where the sources are arranged to produce waves of variable phase delay, may be utilized with proper preprocessing in conjunction with the ultrasound diffraction tomography embodiments disclosed herein. Additionally, while it is preferable to set the center frequency of the source so that the wave length of the pressure field induced in the fluid bath is significantly smaller than the dimensions of the object being investigated, wide latitude is permitted and there is no intention to be bound to a range of frequencies. Moreover, while the fluid bath does not comprise part of the invention, those skilled in the art will appreciate that for ultrasonic frequencies, a liquid bath such as water is preferable. This is not to indicate, however, that air or other fluids cannot be used or might not be preferable for other situations. Of course, it may happen that the natural body fluids could constitute the fluid bath.

0062It will also be appreciated that the analog to digital conversion of the data is not necessary to carry out the invention. If analog signals are used, the filtered backpropagation technique is slightly modified as the Fourier transform and the filter array must be changed to account for the difference in data form. Likewise, other aspects of the preprocessing may be changed, or deleted as desired by those skilled in the art. Additionally, it will be apparent that the filtered backpropagation technique may be altered according to various circumstances. Thus, the filter bank may be split so that array p_n is first filtered by an x-ray tomography deblurring filter and then separately filtered according to a backpropagation filter. Moreover, <img file="EP0109341A2_D0036.tif" />_n can be convolved into a bank of backpropagation filters on a point by point basis rather than as an entire array. The deblurring filtering and/or backpropagation filtering may also be performed in the space domain instead of the frequency domain, as those skilled in the art will appreciate. Further the backpropagation filter may be mathematically altered as required or desired to provide equivalents which may account for changes in the system.

0063Those skilled in the art will also appreciate that a plethora of changes may be made to the numerous embodiments and variations disclosed without deviating from the scope of the invention. Additionally, it will be understood that the ultrasound transmission computed tomography systems and methods embodiments of the disclosed invention have broad application. -For example, the methods and systems disclosed may be used to improve medical ultrasound imaging to provide better detail of possible cancerous growth in body tissue, fetal imaging for pre-natal prognosis, etc. In these situations, the object being investigated is a human body. The source of plane waves is a standard ultrasonic source which is chosen to have a dominant frequency with a corresponding wave length in the human body on the order of one millimeter. The source is placed against or at seme distance from the body and pulsed. In the case of cancer detection, the detectors chosen may be standard piezoelectric crystals (transducers) which are arranged in a linear array located directly against or at some distance from the body on the opposite side of the body from the source. A linear array of detectors may be chosen for use, as a two-dimensional reconstruction of the cancerous tissue is often times of interest (i.e. it is desirable to know whether cancerous tissue exists, and if so, where it is located, but not necessarily its exact shape). The typical analog voltage output of the detector-transducers resulting from the wave diffracting through the body is input into a subsystem which preprocesses the signal which is then filtered according to the filtered backpropagation technique, all as described above. If desired, the preprocessing and the filtered backpropagation technique can both be implemented by computer. The resulting array can be stored in computer memory. Simultaneous with the preprocessing and filtering, the source and detectors are rotated.around the investigated body as is well known in the art. The source is again pulsed and the process repeated. The resulting array is interpolated and stored in memory on a master array, if desired. The entire sequence is repeated until a sufficient number of arrays are produced for proper imaging. In some cases, ten arrays will be enough. In other cases, the sequence will be repeated one hundred times. Of course, the more arrays that are produced, the better the resolution and the profile reconstruction. The master array may then be entered into a graphics system which permits visual display of the reconstruction of the object, either on a screen such as a CRT, or as a hard copy print. The visual display would preferably include both the real and imaginary profiles which could be viewed by one skilled in the art to determine whether cancerous tissue is located in the investigated body.

0064For the investigation of a fetus, it is often desirable to obtain a three-dimensional image. The same procedures outlined above for cancer detection would be utilized except that the source and detectors are rotated in two directions. At each different planar view (orientation), the source and detectors are preferably rotated 360° around the object, stopping at different view angles. Of course, as described above, the preprocessing required is slightly different as a projection calculation is required. The visual display would require a more sophisticated system which might use shading and/or might permit substantial user interface to permit two-dimensional slice and three-dimensional representations.

0065In a technology far removed from the medical arts, the ultrasound diffraction tomography systems and methods disclosed can be used with acoustic rather than ultrasonic 'sources to construct three-dimensional thermal maps of ocean regions by recognizing that the travel time of acoustic waves is effected by water temperatures.

0066Of course, there is no intent to be limited by these applications and those skilled in the art will recognize that the ultrasonic embodiments of the invention, and the many variations thereupon, have broad application for both ultrasonics and sonics. In addition, the methods presented above for ultrasound tomography have direct application in X-ray tomography.

0067In X-ray tomography, the object profile O(x,y) represents the X-ray attenuation profile of a planar section of the object. In the system embodiments of the X-ray tomography invention, the sources are well collimated, quasimonochromatic X-ray sources and the detectors are X-ray detectors. Only the imaginary part of the phase of the detected X-ray field (attenuation) is measured so that phase unwrapping is not necessary. The measured phase is processed by a filtered backpropagation technique identical to that employed in ultrasound tomography. The X-ray tomography systems are applicable to X-ray tomography using both soft and hard radiation. Moreover, in the limiting case.where the wavelength becomes zero, the backpropagation filter reduces to unity and the filtered backpropagation technique reduces to the filtered backprojection . technique of conventional X-ray tomography.

GEOPHYSICAL TOMOGRAPHY

0068The geophysical tomography system and method embodiments of the invention are in many ways closely related to the ultrasound tomography embodiments useful in the medical and other arts. The best modes of the geophysical tomography embodiments which may be used to explore aspects of the subsurface terrain between boreholes or between a borehole and the surface utilize sources, detector-transducers, preprocessing systems, backpropagation filters, and master image arrays. The important differences between the geophysical tomography and ultrasound tomography embodiments stem from the physical situations. Clearly, in geophysical tomography, the formation (body) cannot be rotated relative to the sources and detectors. With the detectors and sources on the surface and in a borehole(s), only two-dimensional approximate reconstructions can be accomplished. Another difference between the systems concerns the sources. In geophysical tomography, the sources are sources'of sonic or electromagnetic energy and not ultrasonics due to the depth of ' investigation required, i.e. the distances from borehole to borehole, or borehole to the surface are much greater than the ultrasonics application and hence larger wavelengths are required to obtain large signal to noise ratios at the receivers. Moreover, in geophysical applications, the sources are sources of spherical waves as the formation (body) is larger in relation to the sources. Due to these and other differences which will become apparent, the preprocessing of the detected signals is different from that employed in ultrasonic tomography. Likewise, the filtered backpropagation technique is altered as the standard X-ray filtering function is radically changed due to the different geometries. The backpropagation filter, however, remains essentially the same.

0069More specifically, Figures 12a, 12b and 12c show three of the preferred geophysical tomography source detector arrangements of the present invention. Figures 12a and 12b show two subsurface electromagnetic or seismic exploration arrangements, while'Figure 12c depicts an electromagnetic or sonic well-to-well tomography arrangement. As seen in Figure 12a, the transmitter sources 214a may be arranged on the surface of the formation 210a. The sources 214a are typically electromagnetic or sonic energy sources and may be considered to be sources of spherical waves as the sources are in contact with the formation (body) 210a. The sources are preferably pulsed in either of two manners: sequentially; or pulsed together with a variable phase delay between adjacent transmitters 214a. In both cases, the electromagnetic or sonic energy pulse propagates through the formation 210a and is detected by detector elements 218a which are located in borehole 215a. The signals detected at the j'th detector (j = 1, 2, ...,-N_j) for the spherical wave transmitted from the _v'th transmitter (_v = 1, 2, ..., N ) are denoted for the sequentially-pulsed case as P_j (t;v); and for the phase delayed case as P_j (t;e ) wherein θ is the phase.tilt of the transmitted' cylindrical wave as shown in Figure 13.

0070Where the sources are pulsed together with a variable phase delay, the detected signals are processed, stored, and A-D converted at break 229a in the same manner as discussed in connection with the ultrasonic embodiments. In the embodiment of the present invention where the sources are pulsed with phase delay; the digitized detected signals P_j(t; θ_n) are Fourier transformed in time and normalized to yield Pj(ω_o, θ_n) according to the equation<maths id="math0028"><img file="EP0109341A2_D0037.tif" /></maths>where U_ois the amplitude of the insonifying cylindrical pulse at its center frequency f_o and ω_o = 2 π f_o is the center frequency in radians/second. In this embodiment, the principle value complex phase θ_j (θ_n) of P_j(ω_o ; e ) is calculated at 234 of Figure 14 using the relations<maths id="math0029"><img file="EP0109341A2_D0038.tif" /></maths><maths id="math0030"><img file="EP0109341A2_D0039.tif" /></maths>The principle valued complex phase array φ_j is then phase unwrapped to yield.the unwrapped complex phase array p_j (θ_n) as indicated at 236.

0071In the situation where the pulse duration of the insonifying pulse is very short, the unwrapped complex phase array p_j (θ_n) is approximately determined by means of the following relations (which are completely analogous to the ultrasonics short pulse situation) :<maths id="math0031"><img file="EP0109341A2_D0040.tif" /></maths><maths id="math0032"><img file="EP0109341A2_D0041.tif" /></maths>In these relationships, A_j (θ_n) and ζ_j (θ_n) are, respectively, the amplitude and phase delay of the pulse received at the j'th detector.

0072Yet another embodiment applies when the Born approximation is employed. In this-case, the array p_j (θ_n) is calculated using the equations<maths id="math0033"><img file="EP0109341A2_D0042.tif" /></maths>where χ is the phase of the insonifying cylindrical wave along the measurement line. The unwrapped complex phase array p_j (θ_n) obtained by either of the above three embodiments is then input to a backpropagation system which will be.described below.

0073When the sources are sequentially pulsed, the detected signals are likewise processed, stored and A-D converted. The digital signals are then Fourier transformed in time and normalized at 231a to generate P. (ω_o ;ν) and the transformed signal is slant stacked at stacker 233a as is well known in the geophysical arts. The signals P_j (ω_o ;ν) are preferably stacked so as to simulate the signal that would be received at the j'th detector if a cylindrical wave had been transmitted (i.e., to simulate <img file="EP0109341A2_D0043.tif" />_j(ω_o ; θ )). Thus, the response of a cylindrical wave is artificially generated by stacking the received signal data. The stacker 233a uses a well-known stacking filter F (θ_n, ν) where θ_n is the phase tilt of the transmitted signal along the transmitter axis as shown in Figure 13. Stacker 233a thus generates <img file="EP0109341A2_D0044.tif" />_j (ω_o; θ_n) which is the same signal generated by the phase delayed array.

0074The calculated unwrapped complex phase array p_j (θ_n) j = 1, 2, ..., N_j corresponds to sample values along the borehole-axis. The total number of sample values N_j is determined by the'number of detecting elements in the receiving array. In order to obtain a good reconstruction of the geological formation lying between the transmitter and receiver axis, it is preferable that a great many sample values of the unwrapped complex phase be determined. This may be accomplished in different ways. If an array of detectors is used, the array must be translated along the measurement line to new locations in order to obtain many sample values unless the array is large enough to span the length of the formation of interest. If the array is moved, the entire set of operations delineated above must be repeated for each location. For the case of a single receiver (N_j=1), the unwrapped complex phase is determined only at a single sample point for any given location of the detector. If a total of M sample values are desired, it is then necessary to perform M procedures with the transducer being successively moved from sample point to sample point between procedures.

0075By means of repeated procedures, the unwrapped complex phase is determined at a total of M sample points along the measurement axis. The M sample values P_m (θ_n), (m = 1, 2, ..., M) play the same role in geophysical tomography as does the unwrapped phase array P_m, m = 1, 2, ... M obtained in ultrasound tomography. In particular, the array p (e ) is input into a backpropagation system 250 to generate a reconstruction of the physical properties of the section of geological formation lying between the transmitter and receiver axes. For the case of sonic energy sources and receivers, the reconstruction yields an acoustic profile analogous to that obtained in ultrasound tomography. For the case of electromagnetic tomography, the reconstruction is of the index of refraction profile of the geological formation.

0076As discussed above, the principle differences between subsurface exploration diffraction tomography and ultrasound diffraction tomography stem from the different geometries involved. The difference in geometries requires a change in the filtered backpropagation technique. In particular, for geophysical applications, the mathematical theory of the prior art discussed in the Background with regard to ultrasonics may be conveniently extended to sonic and elctromagnetic applications. Thus, where it is desirable to reconstruct the profile of a region of geological formation, the Born and Rytav approximations may be employed. Within the approximate models, an equation similar to Equation (10) holds true. In particular, processed signals D_θ0 (y, ω_o) (hereafter referred to as data) taken along the receiver axis the line x = 1_o) as shown in Figure 15, are related to the object profile of the formation 0 (x,y) via the equation<maths id="math0034"><img file="EP0109341A2_D0045.tif" /></maths>In this equation, s, is the unit vector normal to a cylindrical wave radiating out from the transmitter axis (see Figure 13) and 8, is the angle that s, makes with the positive x axis. The unit vectar s is given by<maths id="math0035"><img file="EP0109341A2_D0046.tif" /></maths>while s, is then<maths id="math0036"><img file="EP0109341A2_D0047.tif" /></maths>where again x is restricted to lie in the interval -k&lt;x&lt;k.

0077From Equation (28) it follows that the one-dimensional Fourier transform of the data D_θo (y, ω_o) yields the two-dimensional Fourier transform O (K_x, K_y) of the object profile over semicircular arcs. However, these semicircular arcs differ from the ones occurring in ultrasound tomography. In particular, the locus of points defined by Equation (11) with s and s, given by Equations (29a) and (29b) are semicircles centered at -k s, but for which s_x = <img file="EP0109341A2_D0048.tif" /> · s > 0 rather than s · s, > 0 as was the case in ultrasound tomography.

0078It is interesting to note that the angle θ₀ can only range from 0 to π/2 radians'in the geophysical case. Consequently, the data Dc (y, ω) will not allow O (K_x, K_y) to be sampled throughout a region of K space as large as was obtained in the ultrasound tomography case. This region is shown as area 275 in Figure 16, together with the associated region 285 ( |K | &lt; √2k) obtained in the ultrasound tomography case.

0079In the backpropagation system 250 of the subsurface embodiment, after Fourier transformation, the unwrapped phase array P_m (6_D) is filtered using a filtered backpropagation technique<maths id="math0037"><img file="EP0109341A2_D0049.tif" /></maths>,where<maths id="math0038"><img file="EP0109341A2_D0050.tif" /></maths>and<maths id="math0039"><img file="EP0109341A2_D0051.tif" /></maths>where γ_m is defined in Equation (23), and<img file="EP0109341A2_D0052.tif" /> is the unit vector along the x axis

0080In this case the convolutional (stationary) deblurring filter h_m is no longer the standard X-ray deblurring filter, but depends explicitly on the tilt angle θ_m of the insonifying cylindrical wave. Also, the backpropagation filter, <img file="EP0109341A2_D0053.tif" />_pm' is slightly altered from the ultrasonics case in that it depends on<maths id="math0040"><img file="EP0109341A2_D0054.tif" /></maths>where x_q = qox and y_m = moy are the (x, y) coordinates associated with the m, a sample point in the image array. However, for the special case where θ_n = 0, corresponding to the measurement line being perpendicular to the unit wavevector s₀, the filter functions <img file="EP0109341A2_D0055.tif" />_m and <img file="EP0109341A2_D0056.tif" />_qm reduce to those of the ultrasound diffraction tomography.

0081After the signals have undergone the filtering and backpropagation in a filter bank in backpropagation system 250, the signals are inverse Fourier transformed and interpolated onto a master image array. In the geophysical cases, the arrays need not be rotated as part of the interpolation process as the sources and detectors are not rotated relative to the formation. The master array may be displayed as desired

0082By pulsing the sources 214a at a determined tilt and detecting the produced fields with detector 218a along the borehole axis 215a for the entire area of interest, a partial reconstruction profile of the subsurface formation may be created, as discussed above. In order to reconstruct a more complete profile of the subsurface formation, the phase tilt of sources 214a should cover an entire range of 0 to π/2, wherein the signals for each tilt when separately processed provide a partial reconstruction of the subsurface formation. As those skilled in the art will recognize, a phase tilt of 0 causes the cylindrical wave formed by the spherical waves to propagate straight down into the farnatiar. A phase tilt of π/2 causes the wave to propagate toward the borehole in a manner normal to the surface. The range of phase tilt between 0 and π/2 causes waves entering the formation to propagate at different angles. For each different angle, another partial reconstruction is provided, and the interpolation of those reconstructions onto a master array provides a more complete reconstruction of the subsurface formation.

0083As seen in Figure 12b, the location of the sources 214b and the detectors 218b may be reversed from the previously discussed embodiment shown in Figure lla. The preprocessing at break 229b and thereafter, and the filtering remain identical. The source 214b, however, is moved in the borehole according to well-known techniques. The detectors 218b may be placed on the surface of the formation or they may be buried slightly under the formation surface as desired. Of course, the phase tilt also covers a range of 0 to π/2.

0084Figure 12c shows the borehole to borehole embodiment of the invention. Sources 214c are located in borehole or well 215c and detectors 218c are located in borehole or well 215d. As discussed with regard to the two other geophysical embodiments, the sources may be pulsed with a phase tilt or in a sequential fashion. The preprocessing at 229c and thereafter and the filtered backpropagation technique are identical to the other geophysical cases. Of prime importance in the borehole to borehole (well to well) embodiment is that the phase tilts of the propagating waves range from 0 to r, thereby providing a more complete reconstruction of the formation profile. In fact, as those skilled in the art will appreciate, if both'sources and detectors are placed in both boreholes 215c and 215d, a full two-dimensional reconstruction of the formation profile may be obtained.

0085The best mode embodiments of the methods for reconstructing profiles of subsurface formations are integrally connected with the systems and therefore are analogous to the methods for tomographic reconstruction of an object using ultrasound transmission computed tomography. Thus, at least one wave of electromagnetic or sonic energy is directed at a phase tilt into the formation being investigated. The one or more fields produced by the waves are detected as a function of time and signals representative thereof are produced. The produced signals are preprocessed by calibrating, storing, A/D converting, Fourier transforming, normalizing, and, if necessary, stacking. The preprocessed signals are then Fourier transformed and - filtered by a filtered backpropagation technique which includes a deblurring filter and a backpropagation filter which depend on phase tilt angle. After another Fourier transformation, the resulting array is a partial reconstruction of the subsurface formation profile. By repeating the steps at different phase tilts and interpolating the resulting arrays onto a master array, a more complete reconstruction may be obtained.

0086Those skilled in the art will recognize that many permutations and variations may be made in both the methods and systems disclosed for the geophysical (surface and borehole to borehole) embodiments. For example, not only may detectors and sources be placed in one or more boreholes in the borehole to borehole embodiment, but sources and/or detectors may also be placed on the formation surface. It will be appreciated that if sources are placed on the surface and only detectors are placed in the two boreholes, a 0 to range of phase tilts may be covered, and the reconstruction profile of the formation might be quite adequate. In this manner, the problems associated with locating sources in the borehole can be avoided. In any event, the disclosed embodiments are by way of example only and are not intended to limit the invention in any manner.