Estimation of vector velocity

Abstract

A method and apparatus are provided for estimating the velocity vector of a remotely sensed object or group of objects using either ultrasound or electromagnetic radiation. The movement of the object is determined by emitting and receiving a pulsed field with spatial oscillations in both the axial and transverse directions. Using a number of pulsed emissions and the transverse spatial oscillations the received signal is influenced by transverse motion and a new autocorrelation estimator is used for determining the velocity vector.

Description

FIELD OF THE INVENTION

The invention relates to a method and an apparatus for estimating the velocity vector of a remotely sensed object or group of objects using either sound, in particular ultrasound, or electro-magnetic radiation, eg radar. The movement of the object is determined by emitting and receiving a pulsed field with spatial oscillations in both the axial direction of the transducer and in one or two directions transverse to the axial direction. By using a number of pulse emissions, the inter pulse movement can be estimated and the velocity found from the estimated movement and the time between pulses. The invention is based on the principle of using transverse spatial oscillations for making the received signal influenced by transverse motion, and then using a new autocorrelation estimator for determining the velocity vector. The estimator uses fourth order moments rather than second order moments. The estimation process can also be optimized by using lags different from 1 (one), compensating for the lateral and axial frequencies, and by averaging RF samples.

BACKGROUND OF THE INVENTION

Medical ultrasound is extensively used for studying flow dynamics in the human body by using color flow mapping. The technique displays a color image of the flow superimposed on the normal anatomic B-mode image. The velocity component along the ultrasound beam direction is measured, and a flow transverse to the beam is not displayed. This is shown in FIG. 1, where the flow in the carotid artery and the jugular vein is displayed. The image is acquired with a convex array, and the angles between flow direction and the ultrasound beam changes over the image. Notice the chance of estimated flow direction around the dashed line in both vessels due to the change of angle between the flow and the ultrasound beam. This is one of the main limitations of current ultrasound flow systems, since most vessels are parallel to the skin surface, and therefore it is a problem to get a sufficiently small angle between the flow and the beam. Also the flow is often not parallel to the vessel surface, and it is therefore difficult, if not impossible, to estimate the correct angle and compensate for it [1].

Several authors have attempted to remedy this artifact. Fox [2] suggested using two beams to find the transverse component. The system works well for large transducers and investigations close to the transducer, but the variance of the transverse component increases for situations with large depths and smaller transducers as used in cardiac scanning through the ribs. Trahey and co-workers [3] have suggested using speckle tracking in which a small search region in one image is correlated or compared to a subsequent image. This approach has problems in terms of frame rate, since images are compared, and the resolution of the velocity estimates can be low. Newhouse et al [4] developed a method in which the total bandwidth of the received signal is affected by the transverse velocity. It is, however, often difficult to find this bandwidth due to the inherent noise in the signal.

In this invention a new and improved estimator is presented for the approach described previously in [5] and [6], which makes it possible to estimate the flow vector using a transversally modulated probing field.

SUMMARY OF THE INVENTION

A new estimator for finding the velocity transverse to the ultrasound beam has been developed. The estimator takes into account the influence from the axial velocity both through a subtraction of phase shifts and through tracking as described below under the headline Using Axial Velocity Compensation. Tracking is done between consecutive lines to minimize the axial range of tracking, and thereby minimize the effect of velocity dispersion. The estimator also partly compensates for the difference in the modulation period of the axial and transverse modulation by incorporating a lag different from 1 as used in the traditional autocorrelation approach. The effect of noise is taken into account through averaging RF samples over the pulse length. Estimating the actual center frequency before averaging also takes into account the effect from attenuation. The new estimator is unbiased since the actual mean modulation frequency of the transverse field is estimated before being applied in the estimator. This can be done since the estimator finds the mean velocity in the transverse direction.

Measurement of Transverse Velocities

Conventional velocity estimation systems measure or estimate axial velocity only. By axial velocity is understood the component of the velocity vector in the direction of propagation of ultrasound energy from the ultrasound transducer. In the conventional systems the measurement is performed by emitting a sinusoidal ultrasound pulsed field in one direction a number of times. The returned signal is then sampled at the depth of interest do. The sampled signal for a nearly monochromatic pulse is given by [7]$\begin{array}{cc}r\left(i\right)=\cos \left(2\pi \frac{2v_z}{c}{f}_0{\mathrm{iT}}_{\mathrm{prf}}+\varphi \right)& \left(1\right)\end{array}$where c is the speed of sound, v_z the blood velocity component in the ultrasound direction, f₀ the emitted center frequency, i the pulse number, T_prf the time between pulse emissions, and φ is an arbitrary phase factor that depends on the depth. The frequency of the returned signal $\begin{array}{cc}{f}_p=\frac{2v_z}{c}{f}_0& \left(2\right)\end{array}$is, thus, proportional to the blood velocity, and can be estimated as either a mean frequency in the spectrum of the received signal or a phase shift.

Velocity transverse to the ultrasound beam cannot be estimated from the sampled signal, and a signal influenced by the transverse velocity must be used. The underlying mechanism making it possible to perform axial velocity estimation is the oscillations in the emitted signal. Introducing a transverse spatial oscillation in the ultrasound field causes the transverse velocity to influence the received signal as described in [5] and [6]. The received signal can then be written as: $\begin{array}{cc}{r}_t\left(i\right)=\cos \left(2\pi f_p{\mathrm{iT}}_{\mathrm{prf}}\right)\cos \left(2\pi \frac{v_x}{d_x}{\mathrm{iT}}_{\mathrm{prf}}\right)& \left(3\right)\end{array}$where v_x is the transverse velocity and d_x is the lateral modulation period. The frequency due to the transverse motion can then be written $\begin{array}{cc}{f}_x=\frac{v_x}{d_x}& \left(4\right)\end{array}$Such an approach has been suggested in [6], [5], [8].

The velocities can be both positive and negative and a signal with a one-sided spectrum should be employed to probe the region of interest. This can be found by performing a Hilbert transform on the signals, and for the axial velocity estimation the sampled signal is then $\begin{array}{cc}{r}_q\left(i\right)=\exp \left(\mathrm{j2}\pi \frac{2v_z}{c}{f}_0{\mathrm{iT}}_{\mathrm{prf}}+\varphi \right)& \left(5\right)\end{array}$

A spatial Hilbert transform must be employed, when finding the transverse velocity, and this can be approximated by having two parallel probing beams displayed a distance d_x/4 to yield the spatial quadrature field.

Derivation of the Estimator

The received and sampled in phase spatial quadrature field can be written asr_sq(i)=cos(2πf_piT_prf)exp(j2πf_xiT_prf)  (6) assuming that both the temporal and spatial fields are monochromatic and of unit amplitude. The received field is, thus, influenced by both the axial and the transverse velocity. The influence from the axial velocity on the transverse estimate has previously been compensated for by using tracked data, but any error in tracking due to a poor axial velocity estimate can influence the transverse velocity estimation.Basic Estimator

The axial velocity compensation by tracking can be avoided by employing the more advanced estimator developed in this section.

The temporal Hilbert transform of (6) yields the temporal quadrature spatial quadrature field signal:r_sqh(i)=sin(2πf_piT_prf)exp(j2πf_xiT_prf)  (7)

Rewriting (6) and (7) using Euler's equations gives $\begin{array}{cc}{r}_{\mathrm{sq}}=\frac{1}{2}\left(\exp \left(\mathrm{j2}\pi i{T}_{\mathrm{prf}}\left(f_x-f_p\right)\right)+\exp \left(\mathrm{j2}\pi i{T}_{\mathrm{prf}}\left(f_x-f_p\right)\right)\right)r_{\mathrm{sqh}}=\frac{1}{2}\left(\exp \left(\mathrm{j2}\pi i{T}_{\mathrm{prf}}\left(f_x+f_p\right)\right)-\exp \left(\mathrm{j2}\pi i{T}_{\mathrm{prf}}\left(f_x-f_p\right)\right)\right)& \left(8\right)\end{array}$

Two new signals are then formed from $\begin{array}{cc}{r}_1=r_{\mathrm{sq}}\left(i\right)+{\mathrm{jr}}_{\mathrm{sqh}}\left(i\right)=\frac{1}{2}\left(\exp \left(\mathrm{j2}\pi {\mathrm{iT}}_{\mathrm{prf}}\left(f_x+f_p\right)\right)+\exp \left(\mathrm{j2}\pi {\mathrm{iT}}_{\mathrm{prf}}\left(f_x-f_p\right)\right)\right)+j\frac{1}{2j}\left(\exp \left(\mathrm{j2}\pi {\mathrm{iT}}_{\mathrm{prf}}\left(f_x-f_p\right)\right)-\exp \left(\mathrm{j2}\pi {\mathrm{iT}}_{\mathrm{prf}}\left(f_x-f_p\right)\right)\right)=\exp \left(\mathrm{j2}\pi {\mathrm{iT}}_{\mathrm{prf}}\left(f_x+f_p\right)\right)& \left(9\right)\\ r_2=r_{\mathrm{sq}}\left(i\right)-{\mathrm{jr}}_{\mathrm{sqh}}\left(i\right)=\frac{1}{2}\left(\exp \left(\mathrm{j2}\pi {\mathrm{iT}}_{\mathrm{prf}}\left(f_x-f_p\right)\right)+\exp \left(\mathrm{j2}\pi {\mathrm{iT}}_{\mathrm{prf}}\left(f_x-f_p\right)\right)\right)-j\frac{1}{2j}\left(\exp \left(\mathrm{j2}\pi {\mathrm{iT}}_{\mathrm{prf}}\left(f_x+f_p\right)\right)-\exp \left(\mathrm{j2}\pi {\mathrm{iT}}_{\mathrm{prf}}\left(f_x-f_p\right)\right)\right)=\exp \left(\mathrm{j2}\pi {\mathrm{iT}}_{\mathrm{prf}}\left(f_x-f_p\right)\right)& \left(10\right)\end{array}$

Finding the changes in phase as a function of sample number for the two signals give: $\begin{array}{cc}d{\Theta }_1=\frac{{\mathrm{\Delta \Theta }}_1}{\Delta i}=2\pi T_{\mathrm{prf}}\left(f_x+f_p\right)d{\Theta }_2=\frac{{\mathrm{\Delta \Theta }}_2}{\Delta i}=2\pi T_{\mathrm{prf}}\left(f_x-f_p\right)& \left(11\right)\end{array}$

Adding the two phase changes gives $\begin{array}{cc}d{\Theta }_1+d{\Theta }_2=2\pi 2T_{\mathrm{prf}}{f}_x=4\pi T_{\mathrm{prf}}\frac{v_x}{d_x}& \left(12\right)\end{array}$and subtracting them gives $\begin{array}{cc}d{\Theta }_1-d{\Theta }_2=2\pi 2T_{\mathrm{prf}}{f}_p=4\pi T_{\mathrm{prf}}\frac{2v_x}{c}{f}_x& \left(13\right)\end{array}$

The transverse velocity can thus be found directly from: $\begin{array}{cc}{v}_x=\frac{\left(d{\Theta }_1+d{\Theta }_2\right)d_x}{2\pi 2T_{\mathrm{prf}}}& \left(14\right)\end{array}$and the axial velocity from $\begin{array}{cc}{v}_x=\frac{\left(d{\Theta }_1+d{\Theta }_2\right)c}{2\pi 4T_{\mathrm{prf}}{f}_0}& \left(15\right)\end{array}$

The combination of signals, thus, makes it possible to automatically compensate for the axial and transverse movements, respectively. This is especially important for the transverse estimation due to the rapid variation in phase for an axial movement compared to a transverse movement. An alternative to (15) is to find the axial velocity through a traditional estimator by forming a received beam with conventional beam forming without the transverse modulation. This can potentially yield a higher precision at the expense of an extra receive beam former.

An alternative to combining the signals in (10) would be to find the autocorrelation functions of the two signals and then perform the subtraction or addition on the autocorrelations. The velocity components are then found from the phase shifts in the combined autocorrelations.

The determination of the phase changes for the complex signal can, eg, be-done using the standard autocorrelation estimator [9], [7]. Having the complex signalr(i)=x(i)+jy(i)  (16) the phase change is determined by $\begin{array}{cc}d\Theta \exists =\arctan \left(\frac{\sum _{i=0}^{N-1}y\left(i\right)x\left(i-1\right)-y\left(i-1\right)x\left(i\right)}{\sum _{i=0}^{N-1}x\left(i\right)x\left(i-1\right)+y\left(i\right)y\left(i-1\right)}\right)& \left(17\right)\end{array}$

Using the estimated complex autocorrelation of the signal $\begin{array}{cc}R\exists \left(m\right)=\frac{1}{N-m}\sum _{i=0}^{N-m}{r}^{*}\left(i\right)r\left(i+m\right)& \left(18\right)\end{array}$this can also be stated as $\begin{array}{cc}d\Theta \exists =\arctan \left(\frac{𝒥\left\{R\left(1\right)\right\}}{ℛ\left\{R\left(1\right)\right\}}\right)& \left(19\right)\end{array}$where ℑ{R(1)} denotes the imaginary part of the complex autocorrelation, {R(1)} the real part, and 1)1 is the lag in the autocorrelation. This is equivalent to finding the mean frequency in the power density spectrum given by [7] $\begin{array}{cc}\stackrel{\_}{f}=\frac{{\int }_{-\mathrm{fprf}}^{{\mathrm{fprf}}^{/2}}\mathrm{fP}\left(f\right)df}{{\int }_{-\mathrm{fprf}}^{{\mathrm{fprf}}^{/2}}\mathrm{fP}\left(f\right)df}& \left(20\right)\end{array}$where R(m)⇄P(f) from which the axial velocity is determined by $\begin{array}{cc}{v}_z=\frac{c}{2}·\frac{\stackrel{\_}{f}}{f_0}& \left(21\right)\end{array}$

This estimator, thus, finds the mean velocity. The estimator is also unbiased for white noise added to the input signal r (i) [7].

Finding the phase change by (17) entails finding the arctan of the argument, and the transverse velocity estimation through (14), thus, depends on two arctangents. This creates problems when the phase aliases and a better calculation can be found from using the relation $\begin{array}{cc}\tan \left(A+B\right)=\frac{\tan \left(A\right)+\tan \left(B\right)}{1-\tan \left(A\right)\tan \left(B\right)}\mathrm{then}& \left(22\right)\\ \begin{array}{c}\tan \left(d{\Theta }_1+d{\Theta }_2\right)=\tan (\arctan \left(\frac{𝒥\left\{R_1\left(1\right)\right\}}{ℛ\left\{R_1\left(1\right)\right\}}\right)+\\ \arctan \left(\frac{𝒥\left\{R_2\left(1\right)\right\}}{ℛ\left\{R_2\left(1\right)\right\}}\right))\\ =\frac{\frac{𝒥\left\{R_1\left(1\right)\right\}}{ℛ\left\{R_1\left(1\right)\right\}}\frac{𝒥\left\{R_2\left(1\right)\right\}}{ℛ\left\{R_2\left(1\right)\right\}}}{1-\frac{𝒥\left\{R_1\left(1\right)\right\}}{ℛ\left\{R_1\left(1\right)\right\}}\frac{𝒥\left\{R_2\left(1\right)\right\}}{ℛ\left\{R_2\left(1\right)\right\}}}\\ =\frac{𝒥\left\{R_1\left(1\right)\right\}ℛ\left\{R_2\left(1\right)\right\}+𝒥\left\{R_2\left(1\right)\right\}ℛ\left\{R_1\left(1\right)\right\}}{ℛ\left\{R_1\left(1\right)\right\}ℛ\left\{R_2\left(1\right)\right\}-𝒥\left\{R_1\left(1\right)\right\}𝒥\left\{R_2\left(1\right)\right\}}\end{array}& \left(23\right)\end{array}$where R₁(1) is the complex lag one autocorrelation value for r₁(i) and R₂(1) is the complex lag one autocorrelation value for r₂(i) A similar expression can be derived for the axial velocity, and the estimators are: $\begin{array}{cc}{v}_x=\frac{d_x}{2\pi 2T_{\mathrm{prf}}}\arctan \left(\frac{𝒥\left\{R_1\left(1\right)\right\}ℛ\left\{R_2\left(1\right)\right\}+𝒥\left\{R_2\left(1\right)\right\}ℛ\left\{R_1\left(1\right)\right\}}{ℛ\left\{R_1\left(1\right)\right\}ℛ\left\{R_2\left(1\right)\right\}-𝒥\left\{R_1\left(1\right)\right\}𝒥\left\{R_2\left(1\right)\right\}}\right)\mathrm{and}& \left(24\right)\\ v_z=\frac{c}{2\pi 4T_{\mathrm{prf}}{f}_0}\arctan \left(\frac{\begin{array}{c}𝒥\left\{R_1\left(1\right)\right\}ℛ\left\{R_2\left(1\right)\right\}-\\ 𝒥\left\{R_2\left(1\right)\right\}ℛ\left\{R_1\left(1\right)\right\}\end{array}}{\begin{array}{c}ℛ\left\{R_1\left(1\right)\right\}ℛ\left\{R_2\left(1\right)\right\}+\\ 𝒥\left\{R_1\left(1\right)\right\}𝒥\left\{R_2\left(1\right)\right\}\end{array}}\right)& \left(25\right)\end{array}$These two new estimators compensate for the phase shift in the transverse direction to the velocity component estimated by using fourth order moments compared to the second order estimator used in [6].Compensation for Different Wavelengths

The lateral modulation period will in general be larger than the wavelength of the probing ultrasound pulse. For a given velocity the change in phase for the transverse signal will, thus, be smaller than the change in phase for the axial signal. Optimizing the pulse repetition time for both measurements simultaneously is therefore not possible, and a larger change in phase for the transverse motion must be artificially introduced. This can be attained by using a lag different from one in the autocorrelation estimator as $\begin{array}{cc}d\hat{\Theta }=\frac{1}{k}\arctan \left(\frac{𝒥\left\{R\left(k\right)\right\}}{ℛ\left\{R\left(k\right)\right\}}\right)\hat{R}\left(k\right)=\frac{1}{N-k}\sum _{i=0}^{N-k}{r}^{*}\left(i\right)r\left(i+k\right)& \left(26\right)\end{array}$which can be directly used in (24) and (25).

Using the condition that the phase shift should be the same for v_x=v_z, the lag can be roughly determined by: $\begin{array}{cc}k\approx \frac{d_x}{\lambda }\lambda =\frac{c}{f_0}& \left(27\right)\end{array}$

Often this equation will give a large value for k, and the calculation of the autocorrelation in (26) will include too few values for a low variance estimate. A compromise can be attained by reducing k to obtain both a larger phase shift and a sufficient number of data values for the calculation of {circumflex over (R)}(k)

Optimization in the Case of Noise and Attenuation

The scattering of ultrasound form blood is weak, and a poor signal-to-noise ratio is often found after stationary echo canceling. A prime concern is therefore to make the estimator robust in a noisy environment. This can be attained by averaging the autocorrelation estimate over the length of the interrogating pulse corresponding to performing a matched filtration. The length of the pulse in terms of RF samples is given by $\begin{array}{cc}{N}_p=\frac{M}{f_0}{f}_s& \left(28\right)\end{array}$where M is the number of periods in the pulse (typically 4 to 8) and f_s is the sampling frequency. The autocorrelation estimate is then calculated by $\begin{array}{cc}R\exists \left(k\right)=\frac{1}{\left(N-k\right)N_p}\sum _{i=0}^{N-k}\sum _{n=-N_p/2}^{N_p/2-1}{r}^{*}\left(i,n+N_0\right)r\left(i+k,n+N_0\right)& \left(29\right)\end{array}$

Here r(i,ii) denotes the received signal for the i'th line and its ii'th RF sample. No is the sample number for the position of the velocity estimation, and the averaging is done symmetrically around this position to reduce the effects of velocity dispersion.

The center frequency, f₀ of the ultrasound pulse will change as a function of depth due to the frequency dependent attenuation of the tissue. The number of samples N_p should therefore also changer. The actual mean center frequency can be determined by $\begin{array}{cc}{\stackrel{\_}{f}}_0\left(z\right)=\frac{{\int }_{-f_s/2}^{f_s/2}{\mathrm{fP}}_{\mathrm{rf}}\left(f,z\right)df}{{\int }_{-f_s/2}^{f_s/2}{P}_{\mathrm{rf}}\left(f,z\right)}& \left(30\right)\end{array}$where f_s is the sampling frequency, z is the depth of interest, and P_st(f,z) is the spectrum of the received RF signal around the depth z. {overscore (f)}_st(z) can then be used in (28) to fit the filtration to the measurement situation.Using Axial Velocity Compensation

The estimation process can also be optimized by partially compensating for the axial velocity, when doing the transverse velocity estimation. This can easily be done since the developed estimator essentially subtracts out the phase shift from the axial motion during the transverse estimation process. Whether this is the actual velocity or a smaller phase shift is of no importance. Also the autocorrelation estimator finds the phase shift from one line to the next and therefore a fixed phase shift or equivalently delay can be used. The autocorrelation estimator is then given by. $\begin{array}{cc}R\exists \left(k\right)=\frac{1}{\left(N-k\right)N_p}\sum _{i=0}^{N-k}\sum _{n=-N_p/2}^{N_p/2-1}{r}^{*}\left(i,n+N_0-n_s/2\right)r\left(i+k,n+N_0+n_s/2\right)& \left(31\right)\end{array}$where n_s is the axial movement compensation delay given by $\begin{array}{cc}{n}_s=\mathrm{round}\left(k\frac{2v_z}{c}{T}_{\mathrm{prf}}{f}_s\right)& \left(32\right)\end{array}$rounded off to the nearest number of samples. This ensures that the estimator has the smallest phase shift from the axial movement to compensate for.Estimating the Lateral Modulation Period

The transverse velocity estimate is directly proportional to the lateral modulation period, and a wrong modulation period gives rise to a biased estimate.

The lateral modulation does not have so sharply defined a band pass spectrum as found for the axial pulse. This can be seen from FIGS. 2 and 3. The autocorrelation method, however, estimates the mean frequency, and an estimate of the mean modulation period can, thus, be used to ensure unbiased estimates. This is obtained by first simulating or measuring both the in-phase and quadrature field from a point scatterer moving in front of the transducer at the depth of interest. The spatio-temporal Fourier transform H(f_time, f_space) of the complex probing field is then found as shown on FIG. 3.

The real part of the transformed signal is the in-phase field and the imaginary part is the quadrature component. The spectrum is approximately one sided due to the Hilbert transform relation between the real and imaginary part of the signal. The mean spatial frequency can be then be found from. $\begin{array}{cc}{\stackrel{\_}{f}}_{\mathrm{space}}=\frac{{\int }_{-f_s/2}^{+f_s/2}{\int }_{-f_{\mathrm{sx}}/2}^{+f_{\mathrm{sx}}/2}{f}_{\mathrm{space}}{H\left(f_{\mathrm{time}},f_{\mathrm{space}}\right)}^{2}d{f}_{\mathrm{time}}d{f}_{\mathrm{space}}}{{\int }_{-f_s/2}^{+f_s/2}{\int }_{-f_{\mathrm{sx}}/2}^{+f_{\mathrm{sx}}/2}{H\left(f_{\mathrm{time}},f_{\mathrm{space}}\right)}^{2}d{f}_{\mathrm{time}}d{f}_{\mathrm{space}}}& \left(33\right)\end{array}$assuming that the scattering signal from blood has a white spectrum and is homogeneous over the interrogated region. Here f_sx is the lateral spatial sampling frequency. The mean lateral modulation period is then $\begin{array}{cc}{\stackrel{\_}{d}}_x=\frac{1}{{\stackrel{\_}{f}}_{\mathrm{space}}}& \left(34\right)\end{array}$

FIG. 3 and (33) also indicate how to optimize the measurement field. Ideally a narrow band, single-sided spectrum should be used for having a well-defined measurement situation and thereby a precise velocity estimate. A measure of the spectral spread is obtained by. $\begin{array}{cc}{\sigma }^2\left[i\right]=\frac{{\int }_{-f_s/2}^{+f_s/2}{\int }_{-f_{\mathrm{sx}}/2}^{+f_{\mathrm{sx}}/2}{\left(f_{\mathrm{space}}-{\stackrel{\_}{f}}_{\mathrm{space}}\right)}^2{H\left(f_{\mathrm{time}},f_{\mathrm{space}}\right)}^{2}d{f}_{\mathrm{time}}d{f}_{\mathrm{space}}}{{\int }_{-f_s/2}^{+f_s/2}{\int }_{-f_{\mathrm{sx}}/2}^{+f_{\mathrm{sx}}/2}{H\left(f_{\mathrm{time}},f_{\mathrm{space}}\right)}^2}& \left(35\right)\end{array}$

This can form the basis for the optimization of the spatial quadrature field and the lowest possible value for σ²_jspace will give the best result.

FUNCTIONALITY OF THE INVENTION

The first example in FIG. 4 shows the performance of the new estimator for a purely monochromatic field, where the field is given by equation (6). The lines from origo to ‘+’ indicates the true velocities. The lines from origo to ‘x’ indicate the respective mean values of all the estimates. The ellipses each represent one standard deviation of the estimates in both the transverse and axial direction. The simulation parameters are given in Table 1. Gaussian white noise has been added to the signals. It is seen that the estimator correctly finds the velocity without error. No axial velocity compensation is used in the estimator for this plot, but a similar performance is obtained with axial velocity compensation.

In the second example a pulsed field generated by the Field 11 simulation program described in [10] and [11] is used. This field has been convolved with a two-dimensional random, Gaussian signal for generating a signal with speckle characteristics. The received radio frequency (RF) signals are then generated by selecting the appropriate data from this two-dimensional image according to the transverse and axial velocities. The signals are then shifted an integer number of samples between pulse emissions. The simulation parameters are given in Table 2.

The results of the simulation in the second example are shown in FIG. 5 for a signal-to-noise ratio of 20 dB, and for 100 different estimates. The axial velocity was found through the normal autocorrelation estimation with RF averaging, and the new estimator was used for the transverse velocity. Axial velocity compensation as described was used in the new estimator.

Increasing the noise to yield as signal-to-noise ratio of 0 dB gives the results shown in FIG. 6. Satisfactory results are still seen at this low signal to noise ratio.

To avoid moving an integer number of samples in the venerated speckle pattern, a full Field 11 simulation has been made using roughly 35,000 point-scatterers. The same parameters as before are used. Making plug flow estimation for different angles gives the estimates shown in FIG. 7. No noise has been added to these data, but similar results are obtained for signal-to-noise ratios above 10 dB.

A full simulation with a parabolic flow profile using 36,000 point-scatterers was performed for a vessel with a radius of 5 mm. The peak velocity was 0.5 m/s, and the vessel was perpendicular to the ultrasound beam. A conventional color flow mapping system would in this situation show a velocity of 0 m/s and, thus, show that no velocity is present at this position in the image. The estimates from the new method are shown in FIG. 8. Gaussian noise has been added to the simulation result to obtain a signal-to-noise ratio of 20 dB otherwise the same parameters as before are used. The true velocity profile is shown as the dashed line. Comparing with this the standard deviation of the result is 0.050 m/s, and the mean deviation of the whole profile is −0.0067 m/s. The standard deviation relative to the maximum: velocity is 10.1%. FIG. 9 shows when the simulation has been repeated 50 times. The mean value±one standard deviation is shown.

REFERENCES

• [1] D. J. Phillips, K. W. Beach, and J. Primozich D. E. Strandness. Should results of ultrasound Doppler studies be reported in units of frequency or velocity? Ultrasound Med. Biol., 15:205-212, 1989.
• [2] M. D. Fox. Multiple crossed-beam ultrasound Doppler velocimetry. IEEE Trans. Son. Ultrason., SU-25:281-286, 1978.
• [3] G. E. Trahey, J. W. Allison, and O. T. von Ramm. Angle independent ultrasonic detection of blood flow. IEEE Trans. Biomed. Eng., BME-34:965-967, 1987.
• [4] V. L. Newhouse, D. Censor, T. Vontz, J. A. Cisneros, and B. B. Goldberg. Ultrasound Doppler probing of flows transverse with respect to beam axis. IEEE Trans. Biomed. Eng., BNE-34:779-788, 1987.
• [5] J. A. Jensen and P. Munk. A new method for estimation of velocity vectors. IEEE Trans. Ultrason., Ferroelec., Freq. Contr., 45:837-851, 1998.
• [6] Apparatus and method for determining movements and velocities of moving objects, International patent application, publication number WO 98/00719, published 8 Jan. 1998.
• [7] J. A. Jensen. Estimation of Blood Velocities Using Ultrasound: A Signal Processiiug Approach. Cambridge University Press, New York, 1996.
• [8] P. Munk and J. A. Jensen. Performance of a vector velocity estimator. In Proc. IEEE Ultrason. Symp., pages 1489-1493, 1998.
• [9] C. Kasai, K. Namekawa, A. Koyano, and R. Omoto. Real-time two-dimensional blood flow imaging using an autocorrelation technique. IEEE Trans. Son. Ultrason., 32:458-463, 1985.
• [10] J. A. Jensen and N. B. Svendsen. Calculation of pressure fields from arbitrarily shaped, apodized, and excited ultrasound transducers. IEEE Trans. Ultrason. Ferroelec., Freq. Contr., 39:262-267, 1992.
• [11] J. A. Jensen. Field: A program for simulating ultrasound systems. Med. Biol. Eng. Comp., 10th Nordic-Baltic Conference oil Biomedical Imaging, Vol. 4, Supplement 1, Part 1:351-353, 1996b.

TABLE 1
Simulation parameters for monochromatic simulation.
Parameter          Value
ƒs                 50 MHz
ƒ0                 3 MHz
M                  4
dx                 2 mm
Tprƒ               100 μs
N                  10
snr                {square root over (10)}
Lag compensation k 4

TABLE 2
Simulation parameters for simulations
using the Field II simulation program.
Parameter                 Value
ƒs                        100 MHz
ƒ0                        6 MHz
M                         4
{overscore (d)}x          2.26 mm
Tprƒ                      100 μs
N                         20
snr                       10
Lag compensation k        4
Transducer elements       128
Transducer element width  λ/2
Transducer element height 5 mm
Transducer element kerf   λ/10
Transmit focus            70 mm
Transmit apodization      van Hann window
Receiver focus            fixed at 38.5 mm
Receiver apodization      double sinc
Center of vessel          38.5 mm

Claims

1. An apparatus for measuring the velocity of a moving object or a collection of moving objects, the apparatus comprising: a generator for generating pulses of excitation signals, an emitting transducer for transforming said excitation signals into a beam of wave energy and for emitting said wave energy in a predetermined direction of propagation, a receiving transducer for receiving signals from said moving object or objects generated by interaction with said wave energy emitted from said emitting transducer, wherein said emitting transducer and said receiving transducer having respective sensitivities which in combination give a resulting sensitivity oscillating spatially in a direction transverse to said direction of propagation, wherein the velocity component transverse to the beam is estimated using vx=dx2⁢π⁢ ⁢k2Tprf⁢arctan⁡(𝒥⁢{R1⁡(k)}⁢ℛ⁢{R2⁡(k)}+𝒥⁢{R2⁡(k)}⁢ℛ⁢{R1⁡(k)}ℛ⁢{R1⁡(k)}⁢ℛ⁢{R2⁡(k)}-𝒥⁢{R1⁡(k)}⁢𝒥⁢{R2⁡(k)})(36) as an estimator, and the axial velocity component along the beam is estimated using vz=c2⁢π⁢ ⁢k4Tprf⁢f0⁢arctan⁡(𝒥⁢{R1⁡(k)}⁢ℛ⁢{R2⁡(k)}-𝒥⁢{R2⁡(k)}⁢ℛ⁢{R1⁡(k)}ℛ⁢{R1⁡(k)}⁢ℛ⁢{R2⁡(k)}+𝒥⁢{R1⁡(k)}⁢𝒥⁢{R2⁡(k)})(37) as an estimator, where ℑ{R(k)} denotes the imaginary part of the complex autocorrelation and ℑ{R(k)} the real part, and where R1(k) is the complex lag k autocorrelation value for the signal r1(i) given by r1(i)=rsq(i)+jrsqh(i)  (38)  where rsq(i) is the received and sampled spatial quadrature field and rsqh(i) is the spatial Hilbert transform of rsq(i) and R2(1) is the complex lag one autocorrelation value for the signal r2(i) given by r2(i)=rsq(i)−jrsqh(i),  (39)  where i denotes the pulse-echo line number, c is the speed of the wave energy, f0 is the center frequency of the emitted wave energy, dx is a period of the lateral oscillation of the sensitivity, and Tprf is the time between pulse emissions.

2. An apparatus according to claim 1 where a lag greater than one is used in the estimators.

3. An apparatus according to claim 1 where the received and sampled signal is summed as R^⁡(k)=1(N-k)⁢Np⁢∑i=0N-k⁢∑n=-Np/2Np/2-1⁢r*⁡(i,n+N0)⁢r⁡(i+k,n+N0)(40)

4. An apparatus according to claim 1 where the center frequency is estimated using f_0⁡(z)=∫-fs/2fs/2⁢fPrf⁡(f,z)⁢ⅆf∫-fs/2fs/2⁢Prf⁡(f,z)⁢ⅆf(41)

5. An apparatus according to claim 1 where the axial velocity component is compensated for by using R^⁡(k)=1(N-k)⁢Np⁢∑i=0N-k⁢∑n=-Np/2Np/2-1⁢r*⁡(i,n+N0-ns/2)⁢r⁡(i+k,n+N0+ns/2)(42)

6. An apparatus according to claim 1 where the transverse, modulation frequency used in the velocity estimator is found according to f_space=∫-fs/2+fs/2⁢∫-fsx/2+fsx/2⁢fspace⁢H⁡(ftime,fspace)2⁢ⅆftime⁢ⅆfspace∫-fs/2+fs/2⁢∫-fsx/2+fsx/2⁢H⁡(ftime,fspace)2⁢ⅆftime⁢ⅆfspace(44)

7. A method of estimating the velocity of a moving object or a collection of moving objects, the method comprising emitting an excitation signal of pulses of wave energy as a beam in a predetermined direction of propagation, whereby at least part of the wave energy will interact with the moving object or collection of moving objects, receiving, from said moving object or objects, reflected signals resulting from interaction of emitted wave energy with the moving object or collection of moving objects, wherein the emission of the excitation signal and the reception of reflected signals have respective sensitivities which in combination give a resulting sensitivity oscillating spatially in a direction transverse to the direction of propagation, estimating the velocity component transverse to the predetermined direction, using vx=dx2⁢π⁢ ⁢k2Tprf⁢arctan⁡(𝒥⁢{R1⁡(k)}⁢ℛ⁢{R2⁡(k)}+𝒥⁢{R2⁡(k)}⁢ℛ⁢{R1⁡(k)}ℛ⁢{R1⁡(k)}⁢ℛ⁢{R2⁡(k)}-𝒥⁢{R1⁡(k)}⁢𝒥⁢{R2⁡(k)})(36) as an estimator, and estimating the axial velocity component along the predetermined direction using vz=c2⁢π⁢ ⁢k4Tprf⁢f0⁢arctan⁡(𝒥⁢{R1⁡(k)}⁢ℛ⁢{R2⁡(k)}-𝒥⁢{R2⁡(k)}⁢ℛ⁢{R1⁡(k)}ℛ⁢{R1⁡(k)}⁢ℛ⁢{R2⁡(k)}+𝒥⁢{R1⁡(k)}⁢𝒥⁢{R2⁡(k)})(37) as an estimator, where ℑ{R(k)} denotes the imaginary part of the complex autocorrelation and {R(k)} the real part, and where R1(k) is the complex lag k autocorrelation value for the signal r1(i) given by r1(i)=rsq(l)+jrsqh(i)  (38)  where rsq(i) is the received and sampled spatial quadrature field and rsqh(i) is the spatial Hilbert transform of rsq(i) and R2(1) is the complex lag one autocorrelation value for the signal r2(i) given by r2(i)=rsq(i)−jrsqh(i),  (39)  where i denotes the pulse-echo line number, c is the speed of the wave energy, f0 is the center frequency of the emitted wave energy, dx is a period of the lateral oscillation of the sensitivity, and Tprf is the time between pulse emissions.

8. A method according to claim 7 where the received and sampled signal is summed as R^⁡(k)=1(N-k)⁢Np⁢∑i=0N-k⁢∑n=-Np/2Np/2-1⁢r*⁡(i,n+N0)⁢r⁡(i+k,n+N0)(40)

9. A method according to claim 7 where the center frequency is estimated using f_0⁡(z)=∫-fs/2fs/2⁢fPrf⁡(f,z)⁢ⅆf∫-fs/2fs/2⁢Prf⁡(f,z)⁢ⅆf(41)

10. A method according to claim 7 where the axial velocity component is compensated for by using R^⁡(k)=1(N-k)⁢Np⁢∑i=0N-k⁢∑n=-Np/2Np/2-1⁢r*⁡(i,n+N0-ns/2)⁢r⁡(i+k,n+N0+ns/2)(42)

11. A method according to claim 7 where the transverse modulation frequency used in the velocity estimator is found according to f_space=∫-fs/2+fs/2⁢∫-fsx/2+fsx/2⁢fspace⁢H⁡(ftime,fspace)2⁢ⅆftime⁢ⅆfspace∫-fs/2+fs/2⁢∫-fsx/2+fsx/2⁢H⁡(ftime,fspace)2⁢ⅆftime⁢ⅆfspace(44)