Correction of phasefront aberrations and pulse reverberations in medical ultrasound imaging

[0002] This application claims priority from Provisional application Ser. No. 60/179,196, filed Jan. 31, 2000.

[0003] 1. Field of the Invention

[0004] The present invention is directed to methods for estimating corrections for the image degradation produced in medical ultrasound images by phasefront aberrations and reverberations. The method hence has applications to all situations were ultrasound imaging is used in medicine, and also other similar situations of ultrasound imaging.

[0005] 2. Description of the Related Art

[0006] With ultrasound imaging of objects through complex structures of tissue, the following effects will degrade the images

[0007] i) Variations of the acoustic velocity within the complex tissue structures produce aberrations of the acoustic wavefront, destroying the focusing of the beam mainlobe and increasing the beam sidelobes.*

[0008] ii) Interfaces between materials with large differences in acoustic properties can give so strong reflections of the ultrasound pulse that multiple reflections get large amplitudes. Such multiple reflections are termed pulse reverberations, and add a tail to the propagating ultrasound pulse, which shows as noise in the ultrasound image.*

[0009] The reduced focusing of the beam main lobe reduces the spatial resolution in the ultrasound imaging system. The increase in beam side lobes and the pulse reverberations, introduce additive noise in the image, which is termed acoustic noise as it is produced by the transmitted ultrasound pulse itself. Increasing the transmitted pulse power will hence not improve the power ratio of the signal to the noise of this type, contrary to what is found with electronic receiver noise.

[0010] The materials with largest differences in acoustic properties are muscles, fat, connective tissue, cartilage, bone, air, and the ultrasound transducer itself. Mixtures of fat, muscles, connective tissue and cartilage in the body wall can therefore produce very large phase front aberrations and reverberations. Especially, one will get strong reverberations from the transducer reflections of the returning signals from interfaces of such tissues in the body wall. The mixtures of such tissues in the body wall are therefore the major cause of the degradation found with non-invasive ultrasound imaging in many patients. Reducing the effect of the reverberations and phase front aberrations in the body wall is hence much needed in many applications of medical ultrasound imaging.

[0011] With a two-dimensional transducer array, the effect of the phasefront aberrations can in many situations be reduced by adding corrective delays and gain factors to the signals for the individual array elements, in the following referred to as element signals. This has been presented in many papers. In more complex situations of tissue mixtures, the phasefront aberrations and pulse reverberations can produce modifications of the pulse form. It is less known that such pulse modifications can be corrected by a filter for each of the element signals. Such correction filters gives the most general correction method, and delay/amplitude corrections can be considered as a special case or an approximation of correction filters.

[0012] 2^ndharmonic ultrasound imaging reduces the body-wall reverberations and the phasefront aberrations. The basis for this method is that the wave propagation velocity increases with the pressure amplitude, so that the positive pressure swing gets a higher propagation velocity than the negative pressure swing. This produces a non-linear propagation distortion of the forward propagating pulse that increases with the transmitted pulse amplitude. The distortion first increases with the propagation depth, producing higher harmonic frequency bands in the pulse oscillation. However, because of the power absorption of ultrasound in tissue increases with frequency, the distortion reaches a maximum and finally reduces for large depths.

[0013] Transmitting a pulse with center frequency around f₀≈1.5 MHz, one gets adequate power in the 2^ndharmonic band around 2f₀≈3 MHz for imaging of the heart and other organs in the 3-15 cm range. The power in the higher than the 2^ndharmonic component is so low that it is not useful for imaging.

[0014] The 2^ndharmonic imaging has two advantages above first harmonic imaging around the same receive frequency:

[0015] i) The 2^ndharmonic amplitude is very low as the outbound pulse passes the body wall, so that the 2^ndharmonic components in the body-wall reverberations are low.*

[0016] ii) As the transmitted frequency f₀is low, the first harmonic transmitted beam is less affected by phase aberrations in the body wall. The 2^ndharmonic beam is also generated over a certain volume, which makes the 2^ndharmonic beam less sensitive to phase aberrations, and the field past the focus is also more collimated than for a 1^stharmonic beam at 2f₀.* These two effects hence produce less body-wall reverberations and phase aberrations in the 2^ndharmonic pulse compared to the 1^stharmonic pulse at the same frequency.

[0017] However, when the back-scattered signal passes the body wall on its return, the 2^ndharmonic components are subject to the same amount of phasefront aberrations and pulse reverberations as the 1^stharmonic pulse. Corrections for phasefront aberrations and pulse reverberations in the receive beamformer, hence improves the receiver resolution and acoustic noise, in the same way as for the 1^stharmonic image. Also, there will be some residual phasefront aberrations and body-wall reverberations in the transmitted 2^ndharmonic pulse, which can further be reduced by corrections in the transmit beam former.

[0018] Although the principle of correction for phasefront aberrations and pulse reverberations is well defined, it is in the imaging situation generally difficult to determine the correction filters or the simplified delay and amplitude corrections. The present invention devices two solutions to this problem:

[0019] i) The backscattered signal from point scatterers that are artificially introduced into the body serves as beacons to determine the correction filters. Such point scatterers are either connected to intervention tools that are introduced into the body, for example a biopsy needle, or ultrasound contrast agent bubbles in such a dilute concentration that the signals from individual bubbles can be discriminated from each other. The point scatterers must be so spaced apart that the signals from different point scatterers are clearly differentiable from each other and be so strong that they are differentiable from the tissue signal. To maintain adequate distance between the bubbles, the invention devices to use high transmitted pulse amplitudes to destroy the bubbles in selected image regions, so that for an adequate interval of time after this bubble destruction, new inflow bubbles have adequate distance to each other. Consecutive transmission of high amplitude pulses into the region can be used for repeated destruction of the bubbles so that a continuously changing set of point scatterer bubbles in the image region is obtained. To discriminate the bubble signal from tissue signal one can typically use backscattered frequency components in a band around the the 3^rdor 4^thharmonic component of the transmit frequency band, or sub-harmonic components. Transmission of coded sequences with pulse compression of the received signal can also be used to improve the signal to noise ratio in the received signal from the contrast agent bubbles, and hence the detection of the signal. Spaced apart contrast agent bubbles can also be used as point scatterers on the intervention tool.*

[0020] ii) With the other method one uses stochastic analysis of the back-scattered signal from distributed scatterers with short correlation length compared to the wave length. One general problem for such analysis is that the backscattered signal is corroborated with acoustic noise from pulse reverberations and phase front aberrations, and the invention devices two methods for reduction of such acoustic noise before the determination of the correction filters:

[0021] a) The 2^ndharmonic component of the backscattered element signals, which has reduced acoustic noise, is used for the analysis. However, this requires one filter per element signal, and the invention therefore also devices a simplified method of using the 2^ndharmonic component of the backscattered signal, where filtering of the individual element signals is avoided. In this method, the element signals are compared with a reference signal obtained from the element signals, the reference signal being modified so that the 1st harmonic band in the reference signal is highly attenuated, for example through filtering or pulse inversion techniques that are commonly known.*

[0022] b) The body wall pulse reverberations are fairly stationary in time. By using the backscattered signal from moving or time varying scatterers obtained with multiple transmit pulses with the same focus and beam direction, the temporally stationary acoustic noise is suppressed by highpass filtering each range in the backscattered signal along the pulse number coordinate, so that mainly the signals from the moving scatterers passes the filter for further processing. Typical moving scatterers can be the myocardium or an arterial wall, or scatterers found in blood or other body fluids. To enhance the scattering from such fluids, the invention also devices the use of ultrasound contrast agent to be injected into the body fluid. Time varying scatterers can be ultrasound contrast agent where so high transmit pulses are used that one get destruction of at least some of the contrast agent bubbles between the pulses.*

[0023] The stochastic analysis commonly contains an averaging of signal parameters, where averaging of signal parameters from different depths or possibly also different beam directions typically can be used. Such methods often provides a limited number of samples to average, which gives variance noise in the estimates.

[0024] To improve the estimation robustness and reduce the variance in the estimates, the invention devices a method that uses the backscattered signal from moving or time varying scatterers acquired with multiple transmit pulses with the same focus and beam direction. Signal parameters obtained for each transmit pulse are then averaged for many transmit pulses, possibly in combination with averaging over depth and beam direction, to reduce the variance in the estimates.

[0025] In addition to these basic principles, the invention devices several detailed methods for estimation of the correction filters.

[0026] Strong scatterers off the beam axis, can introduce interference in the correction estimates, and the invention devices methods to reduce the effect of such scatterers, using spatial lowpass filtering of the received signal across the transducer surface, or highpass filtering of estimated correction phases or delays across the transducer surface. Such highpass filtering can conveniently be done by expanding the correction delays in a generalized Fourier series, for example using Legendre polynomials as basis functions, and leaving out the lowest coefficients that relates to offset direction and possible offset focusing of the scatterer.

[0027] The correlation length of the phase aberrations and pulse reverberations along the transducer array surface has a lower bound. One can therefore also truncate the generalized Fourier series at the upper end, reducing the total number of coefficients in the series. The information carrying coefficients in the generalized Fourier series is hence a reduced parameter set that represents the correction filters, and is conveniently estimated in a parameter estimation scheme.

[0028] Often one also find that the correlation length of the phase aberrations and pulse reverberations along the transducer array surface is larger than the array element dimensions, as for example with phased arrays for sector steering of the beam. For this situation the invention devices combination of the element signals from neighboring elements before estimation of the correction filters. This combination reduces the total number of signal channels used in the estimation, hence simplifying the processing and increasing the signal to noise ratio in the resulting channel.

[0029] Other objects and features of the present invention will become apparent from the following detailed description considered in conjunction with the accompanying drawings. It is to be understood, however, that the drawings are designed solely for purposes of illustration and not as a definition of the limits of the invention, for which reference should be made to the appended claims. It should be further understood that the drawings are not necessarily drawn to scale and that, unless otherwise indicated, they are merely intended to conceptually illustrate the structures and procedures described herein.

[0030] In the drawings:

[0031] FIG. 1 shows a typical measurement situation where the ultrasound array in contact with the skin transmits ultrasound pulses through the body wall, which is a mixture of fat, muscles, and connective tissue. The outbound pulse is modified by phase aberrations and pulse reverberations as it passes the body wall. A point scatterer, scatters a spherical wave which is modified by reverberations and phase front aberrations as it returns through the body wall.

[0032] FIG. 2 shows an example of a practical focal plane beam profile that is obtained for a phase aberrated pulse together with a corrected beam profile and the ideal beamprofile obtained without phase aberrations.

[0033] FIG. 3 In a) shows the use of artificially introduced point scatterers on an intervention tool for estimation of correction filters. In b) is shown an example block diagram of a unit for estimation of correction filters both from point scatterers and distributed scatterers, and how the correction filters are used for corrections of phase aberrations and pulse reverberations both in the receiver and the transmit beams.

[0034] FIG. 4 shows in a) a block diagram of a FIR filter for corrections of element signal no. n in the receiver beam former, while in b) is shown a unit for correction of transmitted signal on element no. n in the transmit beam former.

[0035] FIG. 5 shows the use of artificially introduced point scatterers in the form of contrast agent bubbles that are introduced at low concentration in the blood or other body fluids.

[0036] FIG. 6 shows the coordinate system used for calculating the expressions of the element signals.

[0037] FIG. 7 shows in a) the prism type refraction effect of the ultrasound beam produced by a wedge shaped fat layer in the body wall, and in b) the added effect of miss-estimating the location of the center of gravity of the scatterer distribution by a strong scatterer in the beam side lobe.

[0038] FIG. 8 shows a method of hierarchical integration of the neighbor difference in the phase or delay correlations, to reduce the effect of drift in the integration process produced by estimation errors.

[0039] FIG. 1 shows a typical measurement situation where an ultrasound transducer array (101) driven by a transmit beam former (102), transmits a pulsed and focused ultrasound beam through the body wall (103) towards an object (104) to be imaged. The body wall is a heterogeneous mixture of fat, muscles, and connective tissue with differences in acoustic velocity and characteristic impedance. Multiple reflections within the body wall and between structures in the body wall and the transducer, produces a reverberation tail (105) to the transmitted pulse (106) as it passes the wall. Similarly, the variations in the acoustic velocity produce aberrations of the phasefront (107) as the pulse passes the wall.

[0040] At reflection of the pulse from a point scatterer (108) within the object, a wave with spherical wave front (109) is reflected, while the temporal variation of the pulse with reverberation tail is preserved (110). In passing the body wall on its path back to the transducer array, aberrations of the pulse phasefront (111) occur, and additional reverberations of the pulse is found (112).

[0041] In order to focus the receiver beam onto the point scatterer, a standard receiver beam former will delay the individual array signals a time interval calculated from the assumption that the wavefront that arrives from the point scatterer has a spherical shape. Amplitude apodization of the array element signals is also used to reduce the sidelobes in the focal plane. When the real phase front (111) is not ideally spherical due to aberrations, the standard beam former produces a less than ideal focus of the receiver beam. This is illustrated by the focal beam profile 201 of FIG. 2. Using corrective delays on the element signals, one is able to produce an improved focal beam profile as 202 in the Figure. For comparison, the Figure also shows the focal beam profile 203 that is obtained for wave propagation in water.

[0042] In the linear approximation of wave propagation in the tissue, all aspects of the propagation is contained in the Green's function. This function gives the wave-field from a unit point source located at all source points r_fin space. In the temporal frequency domain, we describe the Green's function as g(r,r_f;ω), and it gives the field at the field-point r from a unit point source located at the source-point r_fand radiating continuous, time harmonic waves with temporal angular frequency ω. In the following analysis we relate the Green's function for heterogeneous tissue mixtures to the free-space Green's function g_h(r−r_f;ω) for a homogeneous material as

[0043] g(,_f;ω)=s(,_f;ω)g_h(−_f;ω)  (1)

[0044] where $\begin{array}{cc}{g}_h\left(\underset{\_}{r}-{\underset{\_}{r}}_f;\omega \right)=\frac{e^{-\mathrm{ik}\left|\underset{\_}{r}-{\underset{\_}{r}}_f\right|}}{4\pi \left|\underset{\_}{r}-{\underset{\_}{r}}_r\right|}k=\omega /c& \left(2\right)\end{array}$

[0045] c is the wave propagation velocity in the homogeneous tissue. s(r,r_f; ω) is called the frequency dependent phase-amplitude screen for a scatterer at location r_f. We hence see that s(r,r_f;ω) represents both the phase aberration and reverberation distortion of the body wall on the wave propagation.

[0046] When there is no power absorption in the tissue, it has been shown by M. Fink [1], that time-reversal of the received element signals from a point source in a two-dimensional array, produces the optimal correction of the element transmit signals, for focusing the transmit beam onto the point source. This is equivalent to setting up the transmit beam former in the standard way assuming homogeneous tissue with constant wave velocity, and then filtering the transmit pulses at the element location r on the transducer surface, with the time reversal filter

[0047] H_tr(,_f;ω)=s*(,_f;ω)  (3)

[0048] where * denotes complex conjugation. By reciprocity of the transmit and the receive beams, we also get optimal focusing of the receive beam onto the point target at r_f, by filtering the element receive signals at location r with H_tr(r,r_f;ω).

[0049] The reason one must use a two-dimensional array for corrections of the body wall disturbances is that the body wall produces a two-dimensional disturbance pattern over the transducer array. The correlation length of the disturbances has however a lower limit, so that the width of the elements used in the correction can be larger than typical element widths in many electronic arrays, for example the phased array. For efficient processing, one can in such situations combine the element signals from neighboring, narrow elements to give fewer element signals from larger elements, that are still smaller than the correlation length of the disturbances over the array.

[0050] With uniform acoustic power absorption in the tissue, the gain and phase variations of s* will still correct for the spatial variations in propagation velocity in the tissue and hence reduce the effect of phase aberrations and pulse reverberations. The absorption will, however, attenuate the pulse and also produce some pulse stretching since the power absorption increases with frequency. These effects can be compensated for by a depth variable gain and filter. Strongly heterogeneous absorption will produce spatial amplitude variations in the Green's function that requires special compensation not accounted for by H_trin Eq.(3).

[0051] From the definition in Eq.(1) we see that one can determine s(r,r_f;ω) from the signal from a point scatterer located at r_f, with a received signal that is clearly differentiable from that of the surrounding tissue. However, it is rare to find such point scatterers in the tissue, and it is hence a challenge to determine s(r,r_f;ω).

[0052] In a first aspect, the invention devices the introduction of artificial point scatterers with signal that is clearly differentiable in amplitude or frequency content from the tissue signal, so that g(r,r_fk;ω)=s(r,r_fk;ω)g_h(r−r_fk;ω) can be determined, where r_fkare the positions of the point scatterers that are used for corrected focusing of the transmit and receive beams onto all r_fk's according to Eq.(3). One then can calculate s(r,r_fk;ω)=g(r,r_fk;ω)=g_h(r−r_fk;ω). The distance between the point scatterers then must be so large that the signals from the point scatterers also are differentiable from each other.

[0053] According to the invention, the point scatterers can be introduced on an intervention tool, like a biopsy needle, for improved imaging of the tissue around the needle, especially the biopsy target around the needle tip and also to avoid critical areas around the needle tip when it is introduced, like blood vessels. An example of such a situation is shown in FIGS. 3aand 3b, where 301 shows the intervention tool imbedded in the tissue 302, with two point scatterers 303 and 304 attached to the tool. The point scatterers could also be made as a ring or point indentation in the tool. A two-dimensional transducer array 305 at the surface of the body wall 306 is driven by a transmit beam former 307 shown in FIG. 3b, that is connected to the array via a transmit/receiver switch 308 to emit focused, pulsed ultrasound beams that are direction scanned into the tissue. In this particular depiction the array is shown to transmit a beam 309 that is directed against the point scatterer 303, which backscatters a spherical wave 310. When this spherical wave passes the body wall, it is modified due to phase front aberrations and pulse reverberations as exemplified by 311, before it is picked up by the array 305, that shortly after the pulse transmission is connected via the switch 308 to receiver amplifiers 312 and a correction estimation unit 313, in parrallel to a receive correction unit 314, followed by a receiver beamformer 316 that feeds the output to an ultrasound image processing and display unit 318.

[0054] The correction estimation unit 313 estimates s(r,r_fl;ω) where r_flin this particular example is the location of point scatterer 303. The correction filter transfer functions according to Eq.(3) are transferred via 315 to the receive correction unit 314, where the receive element signals are correction filtered and possibly combined to fewer element signals as described below, before they are transmitted to the receiver beam former 316. This beam former can be of a standard type and has as its output standard ultrasound receive RF-lines 317, possibly for a multitude of neighboring parallel receive directions. The received RF-lines are fed to an ultrasound image processing and display unit 318, that generates and displays full ultrasound tissue and velocity images, as well as ultrasound Doppler spectra, according to known principles.

[0055] The received ultrasound RF-lines 317 could also conveniently be transferred to the correction estimation unit 313 to be used in some estimation algorithms like the reference signal correlation method described in Eqs.(68-74, 90-94) below. Such reference signals could also conveniently be generated within the correction estimation unit itself.

[0056] As the first point scatterer 303 is located near the tip of the tool 301, one will get a corrected image around the tool tip as it is introduced into the tissue. When the beam direction scans, the correction estimation unit 313 will also detect the other point scatterers and produce corrected images in a region around these point scatterers also.

[0057] Corrections of the receive element signals could typically be done in a filter unit shown for one single element signal no n as 401 in FIG. 4a. The frequency response of this filter is for focusing onto r_fgiven as H_n(r_f;ω)=H_tr(r_n,r_f;ω) according to Eq.(3), and the impulse response h_n(r_f;τ) is the inverse Fourier transform of H_n(r_f;ω). This unit takes as input an unmodified signal x_n(t) from element no n, shown as 402 in the Figure, and produce the corrected signal y_n(t) (403) for further processing in the receive beam former. The simplest form of the correction filter is a pure delay function and possible amplitude modification according to the discussion below, in particular in relation to Eqs.(78-94). The receive correction unit could also conveniently combine several channels with minimal nominal delays, for example focusing without direction steering in the elevation direction of the array (y-direction in FIG. 6). The element signals in the elevation direction (y-direction) can then be added in the correction unit for each azimuth position (x-direction of FIG. 6) to reduce the required number of channels in the receiver beamformer 316.

[0058] The output 315 of the correction estimation unit could also conveniently be transferred to the transmit beamformer 307 to modify the transmit pulses to obtain a phase aberration and reverberation corrected beam. The simplest correction is then a delay of the transmit pulse and possibly an amplitude modification. Complex correction according to the full filter in Eq.(3) can for example be obtained for each channel by the unit in FIG. 4b, where the digital representation of the base pulse is stored in a unit 405, and the correction filter impulse response as the inverse Fourier transform of Eq.(3) is set as the coefficients 406 of a transversal filter unit 407. The output of this filter unit is fed to a D/A converter 408 and and the transducer drive amplifier 409 that drives the transducer element no n, 410, possibly via a cable.

[0059] According to the invention, the point scatterers can also be introduced as a dilute concentration of contrast agent bubbles, so that the bubbles in the tissue are so spaced apart that their signals can be discriminated from each other, as illustrated in FIG. 5. This Figure shows as in FIG. 3 the same 2D transducer array 305, body wall 306, tissue 302, and transmit beam 309, which are given the same numeration, where the intervention tool is removed and substituted with a set of contrast agent bubbles 501-505 for the sake of example. To obtain adequately spaced apart contrast bubbles, one can first destroy all contrast bubbles in selected image regions by transmitting high amplitude pressure pulses. The low concentration of contrast bubbles in the inflow arteries will then introduce new bubbles at random locations in the tissue so that in an adequate time interval after the destruction, the new bubbles have adequate distance to each other. With moderate transmit amplitudes one can obtain destruction of a fraction of the contrast agent bubbles for each pulse, which together with the continuous inflow through the tissue blood perfusion produces a time variable set of point scatterers.

[0060] Each contrast bubble reflects when they are hit by an incident ultrasound pulse, a set of non-linearly distorted pulses that propagates as spherical waves, where one in the example illustrates the wave 506 reflected from contrast bubble 505. When propagating through the body wall, the wave is distorted by phase aberrations and pulse reverberations, indicated as 507 in the Figure, similar to 311 in FIG. 3.

[0061] To discriminate contrast bubble signal from the surrounding tissue signal, one can use the sub, 2^nd, 3^rdor 4^thharmonic component of the backscattered signal, that are weak for the tissue signal, but still strong for the contrast bubble signal. The correction filter estimator unit 313 then contains extra filter units that select the selected harmonic components of the back scattered signal in each element signal.

[0062] To avoid unwanted destruction of the contrast agent bubbles, the incident pulse amplitude in the tissue must be kept limited. Maintaining this limit in the incident pulse, one can improve the signal to noise ratio in the received element signals by transmission of coded pulse sequences with subsequent pulse compression in the receiver, and hence reduce the noise in the estimated s(r,r_f;ω).

[0063] Such contrast bubbles are also good candidates for point scatterers to be attached at localized sites on intervention tools.

[0064] In the second aspect of the invention, the present patent devices a method for estimating s(r,r_f;ω) from the back-scattered signal from a distribution of scatterers that can be modeled by a stochastic ensemble with short correlation length so that δ-correlation of the scatterers is an adequate approximation. The acoustic noise in the signal is heavily attenuated before the estimation of s by using 2^ndharmonic processing or signal from moving or time varying scatterers. The signal from such scatterers can also be used to reduce variance and noise in the estimates. Ultrasound contrast agent can also be used in these materials to enhance the scattering from the materials. Destruction of the contrast agent by the incident pulses can also be used to generate time varying scatterers in a region. Specially designed algorithms can be used to estimate parametric representations of s(r,r_f;ω).

[0065] A presupposition is that the scatterers used to estimate s(r,r_f;ω) can be modeled as a δ-correlated ensemble of scatterers. This means that each outcome of the ensemble represents a spatial distribution of the scatterers ν(r) so that averaging over the ensemble we get the following δ-correlation of the scatterer distribution

[0066] <ν(₁)ν(₂)>=σ_ν²(₁)δ(₂−₁)  (4)

[0067] < > denotes ensemble averaging over a mathematically defined outcome ensemble to be further discussed below, and σ_ν²(r₁) is the variance parameter of the distribution at location r₁. We should note that the δ-correlation is an approximation when the correlation length of the scatterer distribution is adequately short to do the approximations from Eq.(14) to Eqs.(15, 35, 70), etc.

[0068] The outcome ensemble is a selected set of scattering distribution functions, where the particular scatterer distribution found in an experiment is considered an outcome or realization of the ensemble. As one in the experimental situation can not observe all outcomes of the ensemble, it is with practical experiments generally a problem to carry through an ensemble averaging as defined above. For spatially stationary scatterer distributions, one can to a limited degree substitute ensemble averaging with spatial averaging over the regions of scatterer positions r₁where s(r,r₁;ω) is practically independent of r₁.

[0069] The present invention also devices an additional method for estimating ensemble averaging of signal parameters, where one uses the signal from moving scatterers observed with consecutive transmit pulses with the same focus and amplitude, separated with so large intervals in time that the scatterers are exchanged in the actual region for each transmit pulse. Hence, the signal from each transmit pulse can be considered as an outcome of the ensemble, and averaging of signal parameters for a particular region can be done by averaging over the signal from the consecutive transmit pulses, which is further described below.

[0070] To obtain time varying scatterers in a tissue region that otherwise has limited movement, and where adequately large blood vessels where the blood can be used as moving scatterers, are missing, one can insert ultrasound contrast agent into the blood vessels, which subsequently flows into the capillaries of the tissue region. With a high amplitude of the incident imaging pulses, some or all of the contrast agent bubbles within the beam can be destroyed. This destruction combined with the new inflow of contrast agent bubbles to the tissue region can be used to produce variation of the contrast agent bubbles in the imaging region between pulses. Temporal averaging of signal parameters of the backscattered signal from the contrast agent for a set of transmitted pulses, can then be used as an estimate of ensemble averages.

[0071] Using the expression for the Green's function in Eq.(1), we get the following expression for the spatial frequency response of the transmitted beam at r₁, when the beam is focused at r_f $\begin{array}{cc}{H}_t\left({\underset{\_}{r}}_1{\underset{\_}{r}}_f;\omega \right)={\int }_{S_t}{d}^2{r}_0\frac{e^{-\mathrm{ik}\left|{\underset{\_}{r}}_1-{\underset{\_}{r}}_0\right|-i\omega {\tau }_{\mathrm{sf}}\left({\underset{\_}{r}}_0,{\underset{\_}{r}}_f\right)}}{2\pi \left|\underset{\_}{r}-{\underset{\_}{r}}_0\right|}s\left({\underset{\_}{r}}_0,{\underset{\_}{r}}_1;\omega \right)\omega {\tau }_{\mathrm{sf}}\left({\underset{\_}{r}}_0,{\underset{\_}{r}}_f\right)={\mathrm{kr}}_f-k\left|{\underset{\_}{r}}_0-{\underset{\_}{r}}_f\right|& \left(5\right)\end{array}$

[0072] where S_tis the transducer surface. τ_sf(r₀,r_f) are the steering and focusing delays for an element r₀at to focus the beam at r_fin a homogeneous medium with wave propagation velocity of c=ω/k.

[0073] The coordinates are illustrated in FIG. 6, where the transducer surface S_tis shown as 601 with a source coordinate at the transducer r₀shown as 602. The vector coordinate r_ffor the focus is shown as 603 with the focal surface S_fshown as the spherical surface 604 with radius r_fcentered at the transducer center. The field-point determined by the field vector r₁is shown as 605. For the calculations below it is convenient to decompose r₁into a component r_1∥ along r_f, shown as 606, and a component r_1⊥ orthogonal to r_fshown as 607, so that r₁=r_1∥+r_1⊥. The focused beam is indicated as 608, where a selected region V_faround the focus is shown as 609.

[0074] Around the focus r_f, one can use the paraxial approximation of the exponent in Eq.(5) that gives $\begin{array}{cc}-\left|{\underset{\_}{r}}_1-{\underset{\_}{r}}_0\right|-r_f+\left|{\underset{\_}{r}}_0-{\underset{\_}{r}}_f\right|\approx -r_1+{\underset{\_}{e}}_1·{\underset{\_}{r}}_0-\frac{r_0^2}{2r_1}-r_f+r_f-{\underset{\_}{e}}_f·{\underset{\_}{r}}_0+\frac{r_0^2}{2r_f}=-r_1+\left({\underset{\_}{e}}_1-{\underset{\_}{e}}_f\right)·{\underset{\_}{r}}_0+r_0^2\left(1/2r_f-1/2r_1\right)& \left(6\right)\end{array}$

[0075] where e₁=r₁/r₁and e_f=r_f/r_fare the unit vectors in the direction of r₁and r_f, respectively. This expression allows us to approximate the transmitted field in Eq.(5) as $\begin{array}{cc}{H}_t\left({\underset{\_}{r}}_1,{\underset{\_}{r}}_f;\omega \right)=\frac{e^{-{\mathrm{ikr}}_1}}{2\pi r_1}{\int }_{S_t}{d}^2{r}_0{e}^{\mathrm{ik}\left({\underset{\_}{e}}_1-{\underset{\_}{e}}_f\right){\underset{\_}{r}}_0+{\mathrm{ikr}}_0^2\left(1/2r_f-1/2r_1\right)}s\left({\underset{\_}{r}}_0,{\underset{\_}{r}}_1;\omega \right)& \left(7\right)\end{array}$

[0076] In the focal surface S_fshown as 604 in FIG. 6, we have r₁=r_f, which gives $\begin{array}{cc}{H}_t\left({\underset{\_}{r}}_1,{\underset{\_}{r}}_f;\omega \right)=\frac{e^{-{\mathrm{ikr}}_1}}{2\pi r_1}{H}_{\mathrm{et}}\left({\underset{\_}{r}}_1,{\underset{\_}{r}}_f;\omega \right)& \left(8\right)\\ H_{\mathrm{et}}\left({\underset{\_}{r}}_1,{\underset{\_}{r}}_f;\omega \right)={\int }_{S_t}{d}^2{r}_0{e}^{\mathrm{ik}\left({\underset{\_}{e}}_1-{\underset{\_}{e}}_f\right){\underset{\_}{r}}_0}s\left({\underset{\_}{r}}_0,{\underset{\_}{r}}_1;\omega \right)& \end{array}$

[0077] The focal beam profile H_etis hence the Fourier transform of the phase-amplitude screen. Close to the focal plane we see from the above that we can approximate H_tas

[0078] H_t(₁,_f;ω)≈e^−ikrF(r₁−r_f)H_et(₁,_f;ω)  (9)

[0079] which will be used in the following analysis.

[0080] The temporal Fourier transform of the received signal at location r on the receiving transducer array, shown as 610 in FIG. 6, from the region V_fshown as 609 around the focus, is $\begin{array}{cc}\begin{array}{c}y\left(\underset{\_}{r};\omega \right)=P\left(\omega \right){\int }_{V_f}{d}^3{r}_{1}2g\left(\underset{\_}{r},{\underset{\_}{r}}_1;\omega \right)\sigma \left({\underset{\_}{r}}_1;\omega \right)\\ =P\left(\omega \right){\int }_{V_f}{d}^3{r}_{1}s\left(\underset{\_}{r},{\underset{\_}{r}}_1;\omega \right)2g_h\left(\underset{\_}{r}-{\underset{\_}{r}}_1;\omega \right)\sigma \left({\underset{\_}{r}}_1;\omega \right)\end{array}& \left(10\right)\end{array}$

[0081] where P(ω) is the temporal Fourier transform of the received pulse from a point scatterer, and σ(r₁;ω) is a source density given as the product of the transmit beam and the scattering density as

[0082] σ(₁;ω)=−k²H_t(₁,_f;ω)ν(₁)  (11)

[0083] We note from Eqs.(4, 11) that σ is δ-correlated in space, i.e.

[0084] <σ*(₁;ω)σ(₂;ω)>=σ²_σ(₁;ω)δ(₂−₁) σ²_σ(₁;ω)=k⁴|H_t(₁,_f;ω)|²σ²_ν  (12)

[0085] For steering and focusing the receive beam onto r_f, we correct y as

[0086] y_f(;ω)=y(;ω)2πr_fe^−iωτ() ωτ_sf(,_f)=kr_f−k|_f|  (13)

[0087] where τ_sf(r,r_f) is the steering and focusing delays that focuses the receive beam onto r_ffor the homomgeneous medium, as in Eq.(5).

[0088] We shall now analyze a situation where the phase-amplitude screen varies so slowly with the scatterer position r₁εV_fthat we can approximate s(r,r₁;ω)≈s(r,r_f;ω), which allows us to take s outside the integral in Eq.(10). Eq.(13) can then be rewritten as $\begin{array}{cc}a)\begin{array}{c}{y}_f\left(\underset{\_}{r};\omega \right)=s\left(\underset{\_}{r},{\underset{\_}{r}}_f;\omega \right)f\left(\underset{\_}{r},{\underset{\_}{r}}_f;\omega \right)\\ f\left(\underset{\_}{r},{\underset{\_}{r}}_f;\omega \right)=P\left(\omega \right){\int }_{V_f}{d}^3{r}_1{g}_f\left(\underset{\_}{r}-{\underset{\_}{r}}_1;\omega \right)\sigma \left({\underset{\_}{r}}_1;\omega \right)\end{array}& \left(14\right)\\ b)\begin{array}{c}{g}_f\left(\underset{\_}{r},{\underset{\_}{r}}_1;\omega \right)=4\pi r_f{e}^{-i\omega {\tau }_{\mathrm{sf}}(r,r_f)}{g}_h\left(\underset{\_}{r}-{\underset{\_}{r}}_1;\omega \right)\\ =\frac{r_f{e}^{-{\mathrm{ikr}}_f+\mathrm{ik}\left|\underset{\_}{r}-{\underset{\_}{r}}_f\right|-\mathrm{ik}\left|\underset{\_}{r}-{\underset{\_}{r}}_1\right|}}{\left|\underset{\_}{r}-{\underset{\_}{r}}_1\right|}\end{array}& \end{array}$

[0089] As σ is a δ-correlated process with zero mean, y_fand f will by the Central Limit Theorem be zero mean Gaussian processes. Hence, all information is contained in the 2^ndorder moments, i.e. the auto-correlation function in space for y_f(r,ω) which is obtained from Eqs.(12, 14) as $\begin{array}{cc}\begin{array}{c}{R}_{\mathrm{yf}}\left(\underset{\_}{r},\underset{\_}{\xi };\omega \right)=⟨y_f^{*}\left(\underset{\_}{r};\omega \right)y_f\left(\underset{\_}{r}+\underset{\_}{\xi };\omega \right)⟩\\ =s^{*}\left(\underset{\_}{r},{\underset{\_}{r}}_f;\omega \right)s\left(\underset{\_}{r}+\underset{\_}{\xi },{\underset{\_}{r}}_f;\omega \right)R_f\left(\underset{\_}{r},\underset{\_}{\xi };\omega \right)\\ R_f\left(\underset{\_}{r},\underset{\_}{\xi };\omega \right)={\left|P\left(\omega \right)\right|}^2{\int }_{V_f}{d}^3{r}_1{g}_f^{*}\left(\underset{\_}{r},{\underset{\_}{r}}_1;\omega \right)g_f(\underset{\_}{r}+\\ \underset{\_}{\xi },{\underset{\_}{r}}_1;\omega ){\sigma }_{\sigma }^2\left({\underset{\_}{r}}_1;\omega \right)\end{array}& \left(15\right)\end{array}$

[0090] where the averaging is done over an infinite number of different outcomes of the scatterer distribution, i.e. ensemble averaging. From Eq.(14) we note that we can use the following approximation $\begin{array}{cc}\begin{array}{c}a)R_{\mathrm{yf}}\left(\underset{\_}{r},\underset{\_}{\xi };\omega \right)=⟨y_f^{*}\left(\underset{\_}{r};\omega \right)y_f\left(\underset{\_}{r}+\underset{\_}{\xi };\omega \right)⟩=e^{-\mathrm{ik}\left|\underset{\_}{r}-{\underset{\_}{r}}_f\right|+\mathrm{ik}\left|\underset{\_}{r}+\underset{\_}{\xi }-{\underset{\_}{r}}_f\right|}\\ e^{\mathrm{ik}\left|\underset{\_}{r}-{\underset{\_}{r}}_1\right|-\mathrm{ik}\left|\underset{\_}{r}+\underset{\_}{\xi }-{\underset{\_}{r}}_1\right|}\\ {\approx }^{\mathrm{ik}\left({\underset{\_}{e}}_1-{\underset{\_}{e}}_f\right)\underset{\_}{\xi }}\end{array}& \left(16\right)\\ \begin{array}{c}b)R_f\left(\underset{\_}{r},\underset{\_}{\xi };\omega \right)\approx {\left|P\left(\omega \right)\right|}^2{\int }_{V_f}{d}^3{r}_1{e}^{\mathrm{ik}\left({\underset{\_}{e}}_1-{\underset{\_}{e}}_f\right)\underset{\_}{\xi }}{\sigma }_{\delta }^2\left({\underset{\_}{r}}_1;\omega \right)\\ =k^4{\sigma }_{\upsilon }^2{\left|P\left(\omega \right)\right|}^2{\int }_{V_f}{d}^3{r}_1{e}^{\mathrm{ik}\left({\underset{\_}{e}}_1-{\underset{\_}{e}}_f\right)\underset{\_}{\xi }}{\left|H_t\left({\underset{\_}{r}}_1,{\underset{\_}{r}}_f;\omega \right)\right|}^2\end{array}& \end{array}$

[0091] Inserting the approximation of H_tclose to the focal plane as given in Eq.(8), we get $\begin{array}{cc}{R}_f\left(\underset{\_}{r},\underset{\_}{\xi };\omega \right)=k^4{\sigma }_{\upsilon }^2{\left|P\left(\omega \right)\right|}^2{\left\{{\int }_{V_f}{d}^3{r}_1{e}^{-\mathrm{ikJ}\left({\underset{\_}{e}}_1-{\underset{\_}{e}}_f\right)\xi }{\left|H_{\mathrm{et}}\left({\underset{\_}{r}}_1,{\underset{\_}{r}}_f;\omega \right)\right|}^2{\left|F\left({\underset{\_}{r}}_1-{\underset{\_}{r}}_f\right)\right|}^2\right\}}^{*}& \left(17\right)\end{array}$

[0092] We note that the two-dimensional, lateral part of the integral is the inverse Fourier transform of |H_et|², while the one-dimensional, longitudinal part of the integral gives a constant. As by Eq.(8) H_etis the Fourier transform of s(r₀,r_f;ω), the above integral can be expressed as $\begin{array}{cc}{R}_f\left(\underset{\_}{r},\underset{\_}{\xi };\omega \right)=k^4{\sigma }_{\upsilon }^2{\left|P\left(\omega \right)\right|}^2{A}^2{\int }_{S_t}{d}^2{r}_{a}s\left({\underset{\_}{r}}_a-\underset{\_}{\xi },{\underset{\_}{r}}_f;\omega \right)s^{*}\left({\underset{\_}{r}}_a,{\underset{\_}{r}}_f;\omega \right)& \left(18\right)\end{array}$

[0093] where A is a constant. Especially we note that for zero lag we have $\begin{array}{cc}\begin{array}{c}{R}_f\left(\underset{\_}{r},\underset{\_}{0};\omega \right)={\left|P\left(\omega \right)\right|}^2{\int }_{V_f}{d}^3{r}_1{\left|g_f\left(\underset{\_}{r},{\underset{\_}{r}}_1;\omega \right)\right|}^2{\sigma }_{\sigma }^2\left({\underset{\_}{r}}_1;\omega \right)\\ ={\left|P\left(\omega \right)\right|}^2{\int }_{V_f}{d}^3{r_1\left(\frac{r_f}{\left|\underset{\_}{r}-{\underset{\_}{r}}_1\right|}\right)}^2{\sigma }_{\sigma }^2\left({\underset{\_}{r}}_1;\omega \right)\\ \approx {\left|P\left(\omega \right)\right|}^2{\int }_{V_f}{d}^3{r}_1{\sigma }_{\sigma }^2\left({\underset{\_}{r}}_1;\omega \right)\\ =k^4{\sigma }_{\upsilon }^2{\left|P\left(\omega \right)\right|}^2{\int }_{V_f}{d}^3{r}_1{\left|H_t\left({\underset{\_}{r}}_1,{\underset{\_}{r}}_f;\omega \right)\right|}^2\end{array}& \left(19\right)\end{array}$

[0094] and utilizing Eq.(18) we get $\begin{array}{cc}{R}_f\left(\underset{\_}{r},\underset{\_}{0};\omega \right)=k^4{\sigma }_{\upsilon }^2{\left|P\left(\omega \right)\right|}^2{A}^2{\int }_{S_t}{d}^2{r}_a{\left|s\left({\underset{\_}{r}}_a,{\underset{\_}{r}}_f;\omega \right)\right|}^2& \left(20\right)\end{array}$

[0095] To find the amplitude-phase screen we express it by its amplitude a_sand phase θ_sas

[0096] s(,_f;ω)=a_s(,_f;ω)e^iθ(;ω)  (21)

[0097] The amplitude of the screen is obtained from Eq.(15) with zero space displacement ξ $\begin{array}{cc}\left|s\left(\underset{\_}{r},{\underset{\_}{r}}_f;\omega \right)\right|=\sqrt{\frac{R_{\mathrm{yf}}\left(\underset{\_}{r},\underset{\_}{0};\omega \right)}{R_f\left(\underset{\_}{r},\underset{\_}{0};\omega \right)}}& \left(22\right)\end{array}$

[0098] We note that R_yf(r,;ω) is the average power of the received element signal. From Eqs.(15, 20) we have that $\begin{array}{cc}{\int }_{S_1}{d}^2{\mathrm{rR}}_{\mathrm{yf}}\left(\underset{\_}{r},\underset{\_}{0};\omega \right)=k^4{\sigma }_{\upsilon }^2{\left|P\left(\omega \right)\right|}^2{A}^2{\left\{{\int }_{S_t}{d}^2{r}_a{\left|s\left({\underset{\_}{r}}_a,{\underset{\_}{r}}_f;\omega \right)\right|}^2\right\}}^2& \left(23\right)\end{array}$

[0099] Solving for the integral of |s|²and inserting into Eq.(20) we get $\begin{array}{cc}{R}_f\left(\underset{\_}{r},\underset{\_}{0};\omega \right)=k^2{\sigma }_{\upsilon }\left|P\left(\omega \right)\right|A\sqrt{{\int }_{S_t}{d}^2{\mathrm{rR}}_{\mathrm{yf}}\left(\underset{\_}{r},\underset{\_}{0};\omega \right)}& \left(24\right)\end{array}$

[0100] which inserted into Eq.(22) gives $\begin{array}{cc}\left|s\left(\underset{\_}{r},{\underset{\_}{r}}_f,\underset{\_}{0};\omega \right)\right|=\sqrt{\frac{R_{\mathrm{yf}}\left(\underset{\_}{r},\underset{\_}{0};\omega \right)}{k^2{\sigma }_{\upsilon }\left|P\left(\omega \right)\right|A\sqrt{{\int }_{S_1}{d}^2{\mathrm{rR}}_{\mathrm{yf}}\left(\underset{\_}{r},\underset{\_}{0};\omega \right)}}}& \left(25\right)\end{array}$

[0101] We note that for Gaussian variables we have $\begin{array}{cc}{R}_{\mathrm{yf}}\left(\underset{\_}{r},\underset{\_}{0};\omega \right)=\frac{\pi }{2}{⟨\left|y_f\left(\underset{\_}{r};\omega \right)\right|⟩}^2& \left(26\right)\end{array}$

[0102] This gives the following modified estimation of the screen amplitude in Eq.(25) $\begin{array}{cc}\left|s\left(\underset{\_}{r},{\underset{\_}{r}}_f;\omega \right)\right|=⟨\left|y_f\left(\underset{\_}{r};\omega \right)\right|⟩\sqrt[4]{\frac{\pi /2}{k^4{\sigma }_{\upsilon }^2{\left|P\left(\omega \right)\right|}^2{A}^2{\int }_{S_t}{d}^{2}r{⟨\left|y_f\left(\underset{\_}{r};\omega \right)\right|⟩}^2}}& \left(27\right)\end{array}$

[0103] To find the phase of the screen, we note that $\begin{array}{cc}\nabla s=\left\{\frac{\nabla a_s}{a_s}+i\nabla {\theta }_s\right\}S\Rightarrow \nabla {\theta }_s=\mathrm{Im}\left\{\frac{\nabla s}{s}\right\}=\mathrm{Im}\left\{\frac{s^{*}\nabla s}{{\left|s\right|}^2}\right\}& \left(28\right)\end{array}$

[0104] We utilize the spatial gradient of the auto-correlation function as

[0105] ∇_ξR_yf(,s*(,_f;ω)∇_ξs(+_f;ω)R_f(,

[0106] +s*(,_f;ω)s(+_f;ω)∇_ξR_f(,  (29)

[0107] Noting that ∇_ξs=∇_rs, we get for zero displacement ξ

[0108] ∇_ξR_yf(,;ω)=s*(,_f;ω)∇_rs(,_f;ω)R_f(,0;ω)+

[0109] |s(,_f; ω)|²∇_ξR_f(,0;ω)  (30)

[0110] which allows us to calculate $\begin{array}{cc}\begin{array}{c}{\nabla }_r{\theta }_s\left(\underset{\_}{r},{\underset{\_}{r}}_f;\omega \right)=\mathrm{Im}\left\{\frac{s^{*}\left(\underset{\_}{r},{\underset{\_}{r}}_f;\omega \right){\nabla }_{r}s\left(\underset{\_}{r},{\underset{\_}{r}}_f;\omega \right)}{{\left|s\left(\underset{\_}{r},{\underset{\_}{r}}_f;\omega \right)\right|}^2}\right\}\\ =\mathrm{Im}\left\{\frac{{\nabla }_{\xi }{R}_{\mathrm{yf}}\left(\underset{\_}{r},\underset{\_}{0};\omega \right)}{R_{\mathrm{yf}}\left(\underset{\_}{r},\underset{\_}{0};\omega \right)}\right\}-\mathrm{Im}\left\{\frac{{\nabla }_{\xi }{R}_f\left(\underset{\_}{r},\underset{\_}{0};\omega \right)}{R_f\left(\underset{\_}{r},\underset{\_}{0};\omega \right)}\right\}\end{array}& \left(31\right)\end{array}$

[0111] We express the auto-correlation function by its amplitude and phase as

[0112] R_yf(,)=A_R(,)e^iθ(a) R_f(,)=A_Rf(,)e^iθ(b)  (32)

[0113] In analogy with Eq.(28) we get $\begin{array}{cc}a){\nabla }_{\xi }{\theta }_R\left(\underset{\_}{r},\underset{\_}{\xi };\omega \right)=\mathrm{Im}\left\{\frac{{\nabla }_{\xi }{R}_{\mathrm{yf}}\left(\underset{\_}{r},\underset{\_}{\xi };\omega \right)}{R_{\mathrm{yf}}\left(\underset{\_}{r},\underset{\_}{\xi };\omega \right)}\right\}& \left(33\right)\\ b){\nabla }_{\xi }{\theta }_{\mathrm{Rf}}\left(\underset{\_}{r},\underset{\_}{\xi };\omega \right)=\mathrm{Im}\left\{\frac{{\nabla }_{\xi }{R}_f\left(\underset{\_}{r},\underset{\_}{\xi };\omega \right)}{R_f\left(\underset{\_}{r},\underset{\_}{\xi };\omega \right)}\right\}& \end{array}$

[0114] Comparing with Eq. (31) we hence get

[0115] ∇_rθ_s(,_f;ω)=∇_ξθ_R(,0;ω)−∇_ξθ_Rf(,0;ω)  (34)

[0116] Integrating this expression with respect to r hence gives the phase of the screen. We note that R_f(r,0; ω) from Eq.(20) is real, so that the imaginary part of the fraction for ∇_ξθ_Rf(r,0;ω) is given by the imaginary part of ∇_ξR_f(r,0;ω). From Eq.(15) we note that $\begin{array}{cc}{\nabla }_{\xi }{R}_f\left(\underset{\_}{r},\underset{\_}{0};\omega \right)={\left|P\left(\omega \right)\right|}^2{\int }_{V_f}{d}^3{r}_1{g}_f^{*}\left(\underset{\_}{r},{\underset{\_}{r}}_1;\omega \right){\nabla }_r{g}_f\left(\underset{\_}{r},{\underset{\_}{r}}_1;\omega \right){\sigma }_{\sigma }^2\left({\underset{\_}{r}}_1;\omega \right)& \left(35\right)\end{array}$

[0117] From Eqs.(15, 33) we hence get $\begin{array}{cc}\begin{array}{c}{\nabla }_{\xi }{\theta }_{\mathrm{Rf}}\left(\underset{\_}{r},\underset{\_}{0};\omega \right)=\mathrm{Im}\left\{\frac{{\nabla }_{\xi }{R}_f(\underset{\_}{r},\underset{\_}{0};\omega }{R_f\left(\underset{\_}{r},\underset{\_}{0};\omega \right)}\right\}\\ =\frac{{\int }_{V_f}{d}^3{r}_1\mathrm{Im}\left\{g_f^{*}\left(\underset{\_}{r},{\underset{\_}{r}}_1;\omega \right){\nabla }_r{g}_f\left(\underset{\_}{r},{\underset{\_}{r}}_1;\omega \right)\right\}{\sigma }_{\sigma }^2\left({\underset{\_}{r}}_1;\omega \right)}{{\int }_{V_f}{d}^3{r}_1\left|g_f\left(\underset{\_}{r},{\underset{\_}{r}}_1;\omega \right)\right|{\sigma }_{\sigma }^2\left({\underset{\_}{r}}_1;\omega \right)}\\ =\frac{{\int }_{V_f}{d}^3{r}_1{\nabla }_r{\theta }_{\mathrm{gf}}\left(\underset{\_}{r},{\underset{\_}{r}}_1;\omega \right){\left|g_f\left(\underset{\_}{r},{\underset{\_}{r}}_1;\omega \right)\right|}^2{\sigma }_{\sigma }^2\left({\underset{\_}{r}}_1;\omega \right)}{{\int }_{V_f}{d}^3{r}_1{\left|g_f\left(\underset{\_}{r},{\underset{\_}{r}}_1;\omega \right)\right|}^2{\sigma }_{\sigma }^2\left({\underset{\_}{r}}_1;\omega \right)}\\ \approx \frac{{\int }_{V_f}{d}^3{r}_1{\nabla }_r{\theta }_{\mathrm{gf}}\left(\underset{\_}{r},{\underset{\_}{r}}_1;\omega \right){\sigma }_{\sigma }^2\left({\underset{\_}{r}}_1;\omega \right)}{{\int }_{V_f}{d}^3{r}_1{\sigma }_{\sigma }^2\left({\underset{\_}{r}}_1;\omega \right)}\end{array}& \left(36\right)\end{array}$

[0118] where θ_gfis the phase of g_fand we have used the relation in Eq.(28) for g_f. From Eq.(14) we have approximated |g_f|²=(r_f/|r₁−r|)²≈1 for r₁part of V_f. From Eq.(14b) we note that

[0119] ∇_rθ_gf(,₁;ω)≈k(_r1−r−_rf−r  (37)

[0120] where e_r1−r−(r₁−r)/|r₁−r| and e_rf−r=(r_f−r)/|r_f−r| are the unit vectors pointing from r to r₁and r to r_f, respectively. Inserting this expression into Eq.(36) gives $\begin{array}{cc}\begin{array}{c}{\nabla }_{\xi }{\theta }_{\mathrm{Rf}}\left(\underset{\_}{r},\underset{\_}{0};\omega \right)\approx k\frac{{\int }_{V_f}{d}^3{r}_1\left({\underset{\_}{e}}_{\mathrm{r1}-r}-{\underset{\_}{e}}_{\mathrm{rf}-r}\right){\sigma }_{\sigma }^2\left({\underset{\_}{r}}_1;\omega \right)}{{\int }_{V_f}{d}^3{r}_1{\sigma }_{\sigma }^2\left({\underset{\_}{r}}_1;\omega \right)}\\ =k\frac{{\int }_{S_f}{d}^2{r}_{1\perp }\left({\underset{\_}{e}}_{\mathrm{r1}-r}-{\underset{\_}{e}}_{\mathrm{rf}-r}\right){\left|H_{\mathrm{et}}\left({\underset{\_}{r}}_1,{\underset{\_}{r}}_f;\omega \right)\right|}^2}{{\int }_{S_t}{d}^2{r}_{1\perp }{\left|H_{\mathrm{et}}\left({\underset{\_}{r}}_1,{\underset{\_}{r}}_f;\omega \right)\right|}^2}\end{array}& \left(38\right)\end{array}$

[0121] where we in the last expression have substituted Eq.(12) with constant σ_ν²and carried though the integration of r₁along r_1∥ defined as 606 in FIG. 6.

[0122] The above expression can be written in the short form as $\begin{array}{cc}{\nabla }_{\xi }{\theta }_{\mathrm{Rf}}\left(\underset{\_}{r},\underset{\_}{0};\omega \right)\approx k\left({\underset{\_}{e}}_{\mathrm{rc}-r}-{\underset{\_}{e}}_{\mathrm{rf}-r}\right)\begin{array}{c}{\underset{\_}{e}}_{\mathrm{rc}-r}=\frac{{\int }_{V_f}{d}^3{r}_1{\underset{\_}{e}}_{\mathrm{r1}-r}{\sigma }_{\sigma }^2\left({\underset{\_}{r}}_1;\omega \right)}{{\int }_{V_f}{d}^3{r}_1{\sigma }_{\sigma }^2\left({\underset{\_}{r}}_1;\omega \right)}\\ =\frac{{\int }_{S_f}{d}^2{r}_{1\perp }{\underset{\_}{e}}_{\mathrm{r1}-r}{\left|H_{\mathrm{et}}\left({\underset{\_}{r}}_1,{\underset{\_}{r}}_f;\omega \right)\right|}^2}{{\int }_{S_f}{d}^2{r}_{1\perp }{\left|H_{\mathrm{et}}\left({\underset{\_}{r}}_1,{\underset{\_}{r}}_f;\omega \right)\right|}^2}\end{array}& \left(39\right)\end{array}$

[0123] where e_rc−r=(r_c−r)/|r_c−r| is the unit vector pointing from r to the center of gravity r_cof σ²(r₁;ω) over V_f. In the last expressions of Eqs.(38, 39) the integration is done laterally over the focal surface S_f. Inserting Eq.(39) into Eq.(34) gives

[0124] ∇_rθ_s(,_f;ω)=∇_ξθ_R(,0;ω)−k(_rc−r−_rf−r  (40)

[0125] This equation can be integrated as $\begin{array}{cc}\begin{array}{c}{\theta }_s\left(\underset{\_}{r},{\underset{\_}{r}}_f;\omega \right)={\int }_0^{x_1}d{\rho }_1{\nabla }_{\mathrm{\xi 1}}{\theta }_R\left(\underset{\_}{\rho },\underset{\_}{0};\omega \right)-\\ k\left(\left|{\underset{\_}{r}}_c-\underset{\_}{r}\right|-\left|{\underset{\_}{r}}_f-\underset{\_}{r}\right|\right)+\mathrm{const}\\ ={\int }_0^{x_2}d{\rho }_2{\nabla }_{\mathrm{\xi 2}}{\theta }_R\left(\underset{\_}{\rho },\underset{\_}{0};\omega \right)-\\ k\left(\left|{\underset{\_}{r}}_c-\underset{\_}{r}\right|-\left|{\underset{\_}{r}}_f-\underset{\_}{r}\right|\right)+\mathrm{const}\end{array}& \left(41\right)\end{array}$

[0126] We note that the term −k(|r_c−r|−|r_f−r|) produces a net refraction and possible change in range of the focus, as illustrated in FIG. 7awhere 701 illustrates a transducer array at the skin 702. The focus delays of the transmit beam is set for focusing at r_f(703) using the assumption of a homogeneous tissue. A wedge shaped fat layer 704 in the body wall gives a prism like refraction of the focus due to its varying thickness and also modifies the focal distance due to the curvature of the layer surface. These effects move the center of gravity of the transmit beam to r_cat 705. The directional unit vectors e_rf−rand e_rc−rfrom a receiver point r at 706 on the transducer surface to r_fand r_care now found as 707 and 708.

[0127] With a homogeneous scatterer distribution the center of gravity of σ²(r₁;ω) over V_fis now found at r_c. With a heterogeneous scatterer distribution, for example a scatterer distribution that contains a strong scatterer shown as 709 in FIG. 7bat a distance from the predicted beam axis 711, the center of gravity of σ²(r₁;ω) over V_fis further moved to r_cat 710.

[0128] Hence, deviation of r_cfrom r_fcan be produced both by a prism like effect of the body wall and by strong scatterers off the beam axis. The last component in θ_swill therefore produce a direction refraction and modify the focusing of the uncompensated beam, but will not increase the beam width and sidelobes. Such a component hence produces a distortion of the uncompensated image. To further analyze this term, we make the following approximation $\begin{array}{cc}\begin{array}{c}k\left({\underset{\_}{e}}_{\mathrm{rc}-r}-{\underset{\_}{e}}_{\mathrm{rf}-r}\right)\approx k\left({\underset{\_}{e}}_{\mathrm{rc}}-{\underset{\_}{e}}_{\mathrm{rf}}\right)\\ =\frac{{\int }_{S_f}{d}^2{r}_{1\perp }\left({\underset{\_}{e}}_{\mathrm{r1}}-{\underset{\_}{e}}_{\mathrm{rf}}\right){\left|H_{\mathrm{et}}\left({\underset{\_}{r}}_1,{\underset{\_}{r}}_f;\omega \right)\right|}^2}{{\int }_{S_f}{d}^2{r}_{1\perp }{\left|H_{\mathrm{et}}\left({\underset{\_}{r}}_1,{\underset{\_}{r}}_f;\omega \right)\right|}^2}\end{array}& \left(42\right)\end{array}$

[0129] where e_rcand e_rfare the unit vectors pointing from the center of the transducer array onto r_cand r_f, respectively. Inserting Eq.(8) for H_et, we get $\begin{array}{cc}=\frac{{\int }_{S_f}{d}^2{r_{1\perp }\left({\underset{\_}{e}}_{\mathrm{r1}}-\underset{\_}{e}\right)}_{\mathrm{rf}}{\int }_{S_t\times S_t}{d}^2{r}_a{d}^2{r}_b{e}^{\mathrm{ik}\left({\underset{\_}{e}}_{\mathrm{r1}}-{\underset{\_}{e}}_{\mathrm{rf}}\right)\left({\underset{\_}{r}}_b-{\underset{\_}{r}}_a\right)}{s}^{*}\left({\underset{\_}{r}}_a,{\underset{\_}{r}}_f;\omega \right)s\left({\underset{\_}{r}}_b,{\underset{\_}{r}}_f;\omega \right)}{{\int }_{S_f}{d}^2{r}_{1\perp }{\int }_{S_t\times S_t}{d}^2{r}_a{d}^2{r}_b{e}^{\mathrm{ik}\left({\underset{\_}{e}}_{\mathrm{r1}}-{\underset{\_}{e}}_{\mathrm{rf}}\right)\left({\underset{\_}{r}}_b-{\underset{\_}{r}}_a\right)}{s}^{*}\left({\underset{\_}{r}}_a,{\underset{\_}{r}}_f;\omega \right)s\left({\underset{\_}{r}}_b,{\underset{\_}{r}}_f;\omega \right)}& \left(43\right)\end{array}$

[0130] We notice that the integration over r₁gives i4π²∇δ(r_a−r_b) in the numerator and 4π²δ(r_a−r_b) in the denominator, which gives $\begin{array}{cc}=\frac{i{\int }_{S_f}{d}^2{r}_a{s}^{*}\left({\underset{\_}{r}}_a,{\underset{\_}{r}}_f;\omega \right){\nabla }_{\mathrm{ra}}s\left({\underset{\_}{r}}_a,{\underset{\_}{r}}_f;\omega \right)}{{\int }_{S_f}{d}^2{r}_a{\left|s\left({\underset{\_}{r}}_a,{\underset{\_}{r}}_f;\omega \right)\right|}^2}& \left(44\right)\end{array}$

[0131] Inserting ∇s from Eq.(28) gives $\begin{array}{cc}\begin{array}{c}=-\frac{{\int }_{S_t}{d}^2{r}_a{\nabla }_{\mathrm{ra}}{\theta }_s\left({\underset{\_}{r}}_a,{\underset{\_}{r}}_f;\omega \right){\left|s\left({\underset{\_}{r}}_a,{\underset{\_}{r}}_f;\omega \right)\right|}^2}{{\int }_{S_t}{d}^2{r}_a{\left|s\left({\underset{\_}{r}}_a,{\underset{\_}{r}}_f;\omega \right)\right|}^2}+\\ i\frac{{\int }_{S_t}{d}^2{r}_a{\nabla }_{\mathrm{ra}}{\left|s\left({\underset{\_}{r}}_a,{\underset{\_}{r}}_f;\omega \right)\right|}^2}{2{\int }_{S_t}{d}^2{r}_a{\left|a\left({\underset{\_}{r}}_a,{\underset{\_}{r}}_f;\omega \right)\right|}^2}\end{array}& \left(45\right)\end{array}$

[0132] The numerator of the last term can be expanded by Gauss theorem as $\begin{array}{cc}{\int }_{S_t}{d}^2{r}_a{\nabla }_{\mathrm{ra}}{\left|s\left({\underset{\_}{r}}_a,{\underset{\_}{r}}_f;\omega \right)\right|}^2={\oint }_{{\Gamma }_t}d{r}_a\underset{\_}{n}{\left|s\left({\underset{\_}{r}}_a,{\underset{\_}{r}}_f;\omega \right)\right|}^2=0& \left(46\right)\end{array}$

[0133] where the last integral is taken along the boundary of S_t, and we have assumed that the integration over S_tis enlarged an ε-area outside the transducer where the excitation and hence s is zero. Going back to the starting equation, Eq.(42), we now see that $\begin{array}{cc}k\left({\underset{\_}{e}}_{\mathrm{rc}-r}-{\underset{\_}{e}}_{\mathrm{rf}-r}\right)\approx -\frac{{\int }_{S_t}{d}^2{r}_a{\nabla }_{\mathrm{ra}}{\theta }_s\left({\underset{\_}{r}}_a,{\underset{\_}{r}}_f;\omega \right)s{\left|\left({\underset{\_}{r}}_a,{\underset{\_}{r}}_f;\omega \right)\right|}^2}{{\int }_{S_t}{d}^2{r}_a{\left|s\left({\underset{\_}{r}}_a,{\underset{\_}{r}}_f;\omega \right)\right|}^2}& \left(47\right)\end{array}$

[0134] Inserting Eq.(40) we see that $\begin{array}{cc}\frac{{\int }_{S_t}{d}^2{r}_a{\nabla }_{\xi }{\theta }_R\left(\underset{\_}{r},\underset{\_}{0};\omega \right){\left|s\left({\underset{\_}{r}}_a,{\underset{\_}{r}}_f;\omega \right)\right|}^2}{{\int }_{S_t}{d}^2{r}_a{\left|s\left({\underset{\_}{r}}_a,{\underset{\_}{r}}_f;\omega \right)\right|}^2}=0& \left(48\right)\end{array}$

[0135] We do not in general know r_c, and the last term in Eqs.(40, 41) is hence unknown. The measurement gives ∇_ξθ_R, and to compensate for beam widening and increase in side lobes by the screen, we use an estimate of θ_sbased on the first term in Eqs.(40, 41). As the measurement is subject to noise, we subtract the DC component over r of ∇_ξθ_R, and use the following estimate for ∇_rθ_s

[0136] ∇_r{circumflex over (θ)}_s(,_f;ω)=∇_ξ{circumflex over (θ)}_R(,;ω)  (49)

[0137] where the hat for θ_Rindicates that any DC-bias has been removed as described. The expression can be integrated to $\begin{array}{cc}\begin{array}{c}{\hat{\theta }}_s\left(\underset{\_}{r},{\underset{\_}{r}}_f;\omega \right)={\int }_0^{x_1}d{\rho }_1{\nabla }_{\xi 1}{\hat{\theta }}_R\left(\underset{\_}{\rho },\underset{\_}{\theta };\omega \right)\\ ={\int }_0^{x_2}d{\rho }_2{\nabla }_{\mathrm{\xi 2}}{\hat{\theta }}_R\left(\underset{\_}{\rho },\underset{\_}{0};\omega \right)\end{array}& \left(50\right)\end{array}$