WO2000016124A1

WO2000016124A1 - Detection of vectorial velocity distribution

Info

Publication number: WO2000016124A1
Application number: PCT/NL1999/000567
Authority: WO
Inventors: Arnoldus Petrus Gerardus Hoeks; Léon Armand Franciscus LEDOUX
Original assignee: Universiteit Van Maastricht
Priority date: 1998-09-10
Filing date: 1999-09-10
Publication date: 2000-03-23
Also published as: NL1010061C2

Abstract

A method is described of remotely detecting the distribution of the vectorial velocity in a diffuse medium (13). By means of at least one transducer (21), a pulse-shaped sound signal (22) is emitted towards the medium (13), a reflection sound signal (23) caused by the medium (13) is detected and converted into an electrical signal (24). The electrical signal (24) is sampled to provide successive signal samples RF (τ, z), which are recorded in a memory (31). The reflection singals RF (τ, z) are rendered analytic, and a correlation function R(T, Z) is calculated for the analytic measured signals. The invention provides a reliable formula for estimating the velocity distribution of the medium on the basis of just a few correlation coefficients of said correlation function R(T, Z).

Description

Detection of vectorial velocity distribution

The present invention relates to a method for trifrp remote detection of the vectorial velocity distribution in a medium, wherein a signal is emitted towards the medium to be studied and wherein the said vectorial velocity distribution is calculated on the basis of a signal received back from the medium. Such a method can be used in particular, but not exclusively, to measure the velocity distribution of blood in a blood vessel, which is why the present invention will hereinafter be explained in the context of this application. It should be noted that the term "velocity distribution" is meant to indicate that the medium to be studied need not move as a plug, in which all the elements of the medium have exactly the same velocity, but that velocity variations can exist in the medium as a whole, so that the velocity components of different elements of the medium will exhibit a certain spread or distribution relative to one another.

Methods of the above-described type are already known per se, for example from the international patent WO 95/32667. This publication describes a method for determining the velocity distribution of blood in a blood vessel, a transducer being employed to emit narrow-beam ultrasonic pulses towards the blood vessel . The signals reflected by the blood are received by the transducer and transferred to a signal-processing arrangement. The reflection signals are first used to estimate the axial velocity component of the blood on the basis of the Doppler shift. This is followed by an estimate of the transverse velocity component of the blood. Combining these two velocity components gives the magnitude of the velocity vector.

The method described in the said publication has several drawbacks .

A first drawback of this known method is that in determining the velocity component of the blood in the direction of the beam axis, use is made of a statistical model to calculate the most probable value of that velocity component. This implies that assumptions are made regarding the stochastic parameters of the signal which was received back. This means a restriction of the applicability of the known method. An additional problem in this context is that sound attenuation in the medium to be studied is a function of the frequency of the sound used, so that a change in spectral bandwidth of the generated signal leads to deviations in the axial velocity component calculated on the basis of the statistical model. A second drawback of the method described in the said publication is that the transverse velocity component is calculated on the basis of the rate at which changes in intensity occur in the signal which was received back. This, however, requires the beam width to be known at the location where the reflections occur, which, however, will in general not be the case .

A third drawback is that the application of this known method requires the existence of flow domains having a laminar flow pattern. A fourth drawback of the method described in the said publication is that only an estimate is provided of the magnitude of the velocity vector, but that regarding the direction of the velocity vector no information is derived from the signal received back. With respect to the axial velocity component, the transverse velocity component, after all, still has a freedom of 360°. In the case of the method according to the said publication, it is therefore necessary for the three-dimensional direction of the velocity vector to be known beforehand, for example by the orientation of the flow duct (blood vessel) to be studied being known.

It is an object of the present invention to provide a method in which the abovementioned drawbacks are absent or at least mitigated. In particular, it is an object of the present invention to provide a method which is also applicable to non- laminar flow.

A further object of the present invention is to provide a method which allows, with as few assumptions as possible, the variables to be measured to be derived directly from the signals received back. In particular, it is an object of the present invention to provide a method which allows not only the magnitude but also the direction of the velocity vector to be derived directly from the signals received back.

The present invention is based on an improved insight into the physical processes which lead to the measurable reflection signal. Based on this improved insight, the present invention provides a model describing the reflection signal, and formulae derived therefrom, with the aid of which the transducer parameters and the variables to be measured can be calculated directly from the measured reflection signals received back. The present invention is thus capable of determining the variables to be measured more rapidly, more accurately and more consistently on the basis of the measured reflection signals, and of deriving more data from the observations. More in particular, the present invention is capable of deriving the local beam width from the observations .

According to an important aspect of the present invention, the shape of the ultrasonic beam generated by the transmitter is assumed to be known. This is a reasonable assumption, since that beam is associated with the transmitter in question in a fixed way, and under normal conditions is reproducible, so that the shape of that beam can be measured in a test arrangement . On the basis of the beam shape, the local width of the beam can then be calculated from the signal received back itself .

An important advantage of the model proposed by the present invention is also that it does not make great demands on the measuring conditions and the measuring setup. With known models it is necessary for the axis of the ultrasonic beam to be oriented perpendicular to the flow direction of the medium to be studied, or that at least the angle between the axis of the ultrasonic beam on the one hand and the flow direction of the medium to be studied on the other hand is accurately known. In the case of the model proposed by the present invention this is not necessary, and the said angle can even be calculated from the measured reflection signals. These and other aspects, characteristics and advantages of the present invention will be clarified in more detail by the following discussion of an embodiment with reference to the drawing, in which: figure 1 schematically illustrates a specific measuring setup for measuring flow velocity in a flow duct; figure 2 schematically illustrates the principle underlying the measurement method; figure 3 shows a two-dimensional representation of an example of the envelope of reflection signals measured successively; figure 4 shows a schematic perspective view of a transducer in order to define a coordinate system.

Figure 1 schematically illustrates a measuring setup 10 to measure the flow velocity in a flow duct 11. In this measuring setup 10, the flow duct 11 forms part of a closed circuit 12 in which a medium 13 is recirculated by means of a pump 16.

The measuring setup 10 further comprises a measuring chamber 20 within which the flow duct 11 is situated. Installed in the measuring chamber 20 is a transducer 21 which is designed to generate a beam of ultrasonic pulses, as will be explained in more detail. The distance between the transducer 21 and the flow duct 11 can be adjusted within certain limits, as can the angle between the beam axis and the axis of the flow duct 11.

Figure 2 schematically illustrates the applied principle underlying the measurement. The transducer 21 generates a brief sound pulse 22 in the direction of the medium 13 to be studied which comprises one or more reflection cores 14. If the emitted sound signal 22 reaches a reflection core 14, this will cause a reflected sound signal 23 to be emitted back to the transducer 21. When the reflected sound signal 23 is intercepted by the transducer 21, the latter generates an electrical signal 24 which is representative of the shape of the sound signal 23 received by the transducer 21; this electrical signal 24 is transferred, as known per se, to a signal-processing arrangement 30 with which a memory 31 is associated. The signal -processing arrangement 30 samples the electrical signal 24 and stores the signal samples in the memory 31.

For further processing of the measurement data, the signal samples are rendered analytic with the aid of a Hubert transformation. The reflection signals are then complex numbers. One advantage of working with analytic signals is that the correlation is not influenced by the carrier wave frequency, as will become evident from the following. This eliminates a source of interference. Possible random or secondary maxima in the correlation (see below) will then occur to a much lesser extent, or at least have a lower level. It will be obvious to those skilled in the art that the Hubert transformation to analytic signal samples can be carried out before the signal samples are stored in the memory 31, or when the signal samples are retrieved from the memory 31 for further processing, or in a separate intermediate step. Hereinafter, no distinction will be made between these options.

If only one reflective body or reflection core 14 is present in the medium, the sound signal 23 received back will normally be in the form of a single pulse, the time between emitting the transmission pulse 22 and receiving the reception pulse 23 being, on the basis of the sound velocity in the respective medium 13, a yardstick for the spatial distance between the transducer 21 and the reflective body or reflection core 14. The medium 13 is a diffuse medium, i.e. it contains a multiplicity of randomly distributed particles 14 which may act as reflection cores. In practice this means that the signal 23 received back by the transducer 21 is a combination of many reflection pulses.

Reflections which arrive almost simultaneously at the transducer 21 may interfere with one another. They may reinforce one another (constructive interference) or on the other hand cancel one another (destructive interference) . Figure 3 illustrates an example of the envelope of a set of successive interference signals of this type, obtained in the example shown by measurement in a human carotid artery. In figure 3 , the strength of the electrical signal 24 emitted by the transducer 21, or alternatively the strenth of the reflection signal 23 received by the transducer 21, is plotted as a function of time (horizontal) . Since, as explained in the above, the time corresponds to distance, and in fact the distance is much more interesting to the observer than the time, the distance scale is shown at the horizontal axis. The vertical axis corresponds to the successive received reflection signals.

In figure 3, the strength of the electrical signal 24 emitted by the transducer 21 is represented by blackness (grey scale value) of a dot in the figure. Herein, a strong signal is represented by a dark dot, while a weak signal is represented by a pale dot .

Each horizontal line of measured points in figure 3 represents the envelope of the electrical response of the transducer 21, after a single transmission pulse 22 has been emitted by the transducer 21. During a measuring procedure, the emission of a single transmission pulse 22 by the transducer 21 is repeated a number of times (at a predefined repeat rate) ; the electrical responses of the transducer 21 which are associated with the successive transmission pulses 22 (a number of horizontal lines of experimental points) are shown in figure 3 underneath one another; placed along the left-hand axis of figure 3 are serial numbers of these transmission pulses. In general terms it can be seen in figure 3 that the constructive and destructive interference pattern (response to a transmission pulse) of the reflection signal changes as a function of time. This is a consequence of the changes in the relative position of the individual reflection cores 14. More in particular, the individual reflection cores are shifted owing to the flow of the medium 13, and the relationship relative to one another between the reflection cores will change. These changes take place on a timescale which is larger than the repeat period of the transmission pulses. This relatively slow aspect of the change finds its expression in the relatively large vertical dimensions of the spots in the graph representation of figure 3. More in particular, this vertical length is indicative for the transverse component of the flow velocity. The oblique orientation of these spots implies a horizontal shift of the interference pattern as a function of time, which corresponds to the axial component of the flow of the medium.

In the following, a model proposed by the present invention will be derived for the purpose of interpreting the measured signals. The formulae used herein are collected together at the end of the description.

For definining a coordinates system, reference is made to figure 4. Figure 4 shows a transduceer array 121 comprising a multiplicity of juxtaposed transducers 21. Each transducer 21 is able to emit a sound beam 22, whose beam axis 25 is assumed to be the Z-axis. The X-axis is defined by a line connecting the transducers 21, so that the sound beams of all the transducers 21 are collectively situated in the XZ-plane. The Y-axis is defined perpendicular to the XZ-plane defined by the transducers 21. In the following derivation it will be assumed that use is made of two different transducers 21_: and 21₂ of the transducer array 121, the relative distance in the X- direction between those two transducers being indicated as ΔX and in the Y-direction as ΔY; in the coordinate system described, ΔY is equal to zero. These two transducers 21_x and 21₂ are chosen such that their sound beams partially overlap one another in space. This results in an initial correlation level between the measured signals RF_X provided by the first transducer 21_x and the measured signals RF₂ provided by the second transducer 21₂. This initial correlation level is used in determining the motion components.

Under the assumption that the local distribution of the reflection cores 14 within the medium 13 is constant during passage of a sound wave, the first reflection signal RF_1#0(t) received back by transducer 21_: (t being the time after a sound pulse has been emitted by said transducer) can be described by formula (1) . Therein, g(x,y,z) denotes the (three-dimensional) spatial distribution of the reflection cores 14 in the medium 13, and f_t(x,y,z) is the three-dimensional sensitivity function of the transducer in question. The sensitivity function describes the three- dimensional pressure wave volume which contributes to the reflection signal at time t. It is further assumed that, even though the distribution of the reflection cores 14 in the medium 13 is stochastic, their movement is coherent within the three-dimensional sensitivity function.

Assume that array 121 emits sound pulses with a spacing of ΔT_EP between successive emission pulses. The magnitude of the displacements of the medium in the X- , Y- and Z-direction which occur during that interval will respectively be denoted by S_x, S_γ and S_z. The corresponding velocity components are V_X=S_X/ΔT_EP, V_Y=S_Y/ΔT_EP and V_Z=S_Z/ΔT_EP, respectively. The (T+l)th captured reception signal RF_{2 T}(t) of the other transducer 21₂ can then be described by formula (2) . The term P_nT+p₁₂ in formula (2) is the time difference (unit ΔT_EP) between the transmission pulses associated with the reflection signals RF_lι0(t) and RF_{2 T}(t). The variable p_xl is the time difference (unit ΔT_EP) between two successive reflection signals captured by the same transducer. If each emitted sound pulse results in a captured signal, it is the case that p_n=l . If, however, a number n of reflection signals must be disregarded in between two reflection signals to be captured, it is the case that p₁₁=n+l. The variable p₁₂ is the time difference (unit ΔT_EP) between the mth reflection signal captured by the transducer 21_x and the mth reflection signal captured by transducer 21₂. If both transducers capture their signals simultaneously, it is the case that p₁₂=0.

It should be noted that the derivation given here applies to measurements carried out with the aid of two transducers 21-_L and 21₂. The derivation is equally valid, however, if only one transducer is employed; the formulae then simplify because ΔX=ΔY=0 and p₁₂=0, so that RF_X is equal to RF₂.

The correlation between the initial reflection signal RF_1/0(t) captured by the first transducer 21-, and the depth-shifted (T+l)th reflection signal RF_{2 τ} (t+Z/f_s) captured by the second transducer 21₂, denoted by R(T,Z), is defined by formula (3), wherein T denotes a time difference (unit ΔT_EP) , Z denotes a spatial differential distance, f_s is the sampling frequency of the reflection signal, * means complex conjugate, and <> denotes mathematical expectation.

By substituting formulae (1) and (2) into formula (3), the correlation R(T,Z) can now be written as formula (4) .

Under the assumption that the reflection cores 14 are δ -correlated in space, the expectation of the spatial distribution of the displaced reflection cores can be written as formula (5) , wherein δ is the Dirac delta function, and

<g²> is the mean square of the spatial distribution of the reflection cores. Thus the correlation function of formula (4) can be simplified to give formula (6) .

As the two sensitivity functions f_t* and f_{t+z fs} of formula (6) are located very close together in space, it is assumed that they are of the same shape.

In the following, first an analytical description will be given of the sensitivity function of a transducer. The shape of the sensitivity function in the axial direction can be derived from the spectral distribution of the sound signal, whereas in the XY-plane local characteristics of the sound beam restrict the boundaries of the sensitivity function. These contributions to function are independent of one another and will be discussed independently of one another.

The spectral power density P(f) of an (analytic) reflection signal has an approximately Gaussian shape and in the normalized form can be written as formula (7) , wherein

P(f) is the normalized (P(f_c)=l) power at a particular frequency f, f_c is a central carrier wave frequency, and BW_EQ is the equivalent bandwidth of the spectral distribution. With "equivalent bandwidth" is meant the width of a rectangular distribution having the same area as the curve defined by the actual spectral distribution, the height of the equivalent function being equal to the height of the original function at f_c . For a Gaussian curve, the equivalent bandwidth can be written as formula (8) , wherein BW_dB denotes the bandwidth at a specified level with respect to the maximum in decibels.

Consider a reflection signal originating from a single reflection core which is located at a distance d (=0.5*c*t_d, wherein c is the sound velocity in the medium) in front of the transducer. If the spectral distribution satisfies formula (7) , a reflection signal of this type can be written as formula (9) , where Ψ is an arbitrarily chosen initial phase. The normalized axial sensitivity function f_t(z) can then be written as formula (10), where z_t (=0.5*c*t) is the axial distance corresponding to t seconds after emission of the transmission pulse.

In the directions perpendicular to the beam axis, the beam dimensions contribute to the sensitivity function in the XY-plane. In these directions, the sensitivity functions can be approximated by Gaussian functions in accordance with formulae (11) and (12). Therein, wx_dB(z_t) and wy_dB(z_t) denote the beam width in the x- and y- direction, respectively, at a depth z_t, at a specified level in dB relative to the peak value. Since the total sensitivity function f_t(x,y,z) can be written as f_t(x,y,z) = f_t (x) • f_t (y) • f_t (z) , this function, by multiplying the formulae (10) , (11) and (12) , can now be written as formula (13) .

Formula (13) only applies to plane waves, for which all the pressure waves in the XY-plane are in phase. In the case of focussed transducers, however, a curvature of the emitted wavefront occurs owing to path length differences (ΔL) between the sound waves emitted from the centre of the transducer and from the edges. This curvature can be approximated parabolically by formula (14) , the curvature at position (X,Y,z_t) being defined by the additional path length Δ _X and ΔL_Y in the x- and y-direction, respectively. The sensitivity function f_t' (x,y,z) corrected for this curvature can be written as formula (15) .

This modified sensitivity function f_t' (x,y,z) can be written, according to formula (16) , as a product of independent x- , y- , and z-components .

Thus the correlation function R(T,Z) in formula (6) can be factored into x- , y- , and z-components according to formula (17) . Herein, R_X(T,Z) is the x-component (transverse direction) of the correlation function, R_y(T,Z) is the y-component (height direction) of the correlation function, and R_Z(T,Z) is the z-component (axial direction) of the correlation function.

By combining the formulae (6) and (17) , the x-component R_X(T,Z) of the correlation function can be written as formula (18) . This can be simplified as indicated in formula (19). Herein, C_Bx(z_t) is a negative constant which is a function of the characteristics of the sound beam in transverse direction at depth z_t .

In a similar manner, the y-component R_y(T,Z) of the correlation function can be simplified as indicated in formula (20), and the z-component R_Z(T,Z) of the correlation function can be written as formula (21) . By filling in the axial component of formula (16) , formula (21) can be rewritten as formula (22) .

By writing out formula (17) , and after normalization (R(0,0) = 1), the correlation function can be written as the analytic model of formula (23) .

The correlation of the reflection signals is affected by noise. Assuming that the noise in a reflection signal also has a Gaussian spectral distribution, and given that the noise between successive reflection signals is uncorrelated, the correlation function can be rewritten as formula (24), wherein the argument of R(T,Z) can be written as formula (25), and wherein the magnitude |R(T,Z) | can be written as formula (26) , wherein S is the normalized (S+N=l) signal power and N_M is the normalized noise power. If ΔX = 0, ΔY = 0, p₁₂ = 0 and T = 0, the normalized noise power is equal to N, and in all other cases the normalized noise power is 0.

In the above, a theoretical model has been described which describes a correlation function for the reflection signals received. Hereinafter an explanation will be given how, conversely, data can be derived if the correlation function is known in a situation in a practice situation. First, however, values will have to be found in practice for that correlation coefficient R(T,Z) at various values of T and Z. More in particular, those values must be calculated from the reflection signals received. By way of example, this can be done as follows.

The reflection signals received are sampled; the successive sampling points are denoted by the serial number ζ . The first sampling point of each reflection signal captured is denoted by ζ = 1. The successive sound pulses are distinguished by the serial number T. Thus the fth sampling point of a reflection signal after a rth sound pulse is denoted by RF(τ,$^") . Prior to processing, the reflection signals are rendered analytic via a Hubert transformation. Such a transformation can be written as formula (27), wherein A(f) is the analytic frequency spectrum corresponding to the frequency spectrum S(f) of an actual signal, as will be explained below.

Assume that S(f) is the frequency spectrum which corresponds to the sampled reflection signal RF(τ,f)

(Fourier transformation) . This frequency spectrum can be regarded as a set of frequency components f_A of relative strength S (complex numbers) . These include negative frequency components, which in the course of conversion into analytic form are eliminated by replacing the respective magnitudes S₁ by magnitudes A₁=0. To retain the power of the signals, the respective magnitudes S₁ are doubled for the positive frequency components, i.e. they are replaced by magnitudes A₁=2S₁. The magnitude associated with frequency f=0 is not changed. By performing an inverse Fourier transformation subsequently, an analytic signal RF_A is obtained from RF. Since two transducers are employed, there are also two of those two-dimensional matrices. Each of the two corresponds to one of the two transducers. If only one transducer is employed, the two matrices are equal to one another . To allow stable correlation coefficients to be determined, it is necessary to average. To this end, an identical data window in the two matrices is defined having NZ sampling points in the f-direction and NT sampling points in the T-direction. The data window associated with the signals captured by transducer 21_x is denoted by w₁(τ,π . The data window associated with the signals captured by transducer 21₂ is denoted by w₂(τ,f) . For the first sampling point of the first signal segment in the data windows it is the case that τ=l and ζ=l . The correlation coefficient R(T,Z) of signals from the two data windows, the difference between the serial numbers of the signals in both data windows being T and between the sampling points being Z, is defined according to formula (28) . When calculating the correlation coefficients

R(T,Z) from the sampled analytic reflection signals it is necessary to find a compromise between precision of the correlation coefficients on the one hand, and spatial and temporal resolution on the other hand. A better resolution is obtained by selecting smaller values for NT and/or NZ, whereas better precision of the correlation coefficients is achieved by taking larger values for NT and/or NZ .

In the above, a derivation has been described of a model which describes the two-dimensional correlation function R(T,Z) in the presence of noise (formulae 24 to 26 inclusive) , or if noise may be neglegted (formula 23) . On the basis of this model, the signal parameters and beam shape parameters can be calculated from the measurement data, taking into account the estimated correlation coefficients, as will be described hereinafter. Herein, the estimated value of a parameter x will be indicated by x. In addition to the presented formulae for estimating the parameters, yet other formulae for estimating the same parameters can be derived from the model . These formulae will, however, lead to the same result.

The central frequency f_c can be estimated according to formula (29), wherein Z_x ≠ Z₂. The best results are obtained if Z_x = 1 and Z₂ = 0.

On the basis of the arguments of the correlation coefficients, the axial movement S_z can be estimated by means of formula (30), where ^ ≠T₂. The best results are obtained if T_x = 1 and T₂ = 0. An important aspect of formula (30) is that it is a relatively simple formula requiring relatively little computation time. This formula, however, is subject to a restriction, as is each formula which is based on the argument of a complex number, namely the so-called "aliasing" effect, caused by the fact that the argument of a complex number can only assume values between -180° and +180°. In practice, this means an upper limit for the velocity which can be estimated with the aid of formula (30) . The present invention provides a solution even for this restriction, achieving this by providing the alternative formula (31) for the axial movement S₂, wherein Z ≠ Z₂ ≠ Z₃ ≠ Z_1# and T ≠ 0. This formula (31) is based on the magnitude of the correlation coefficients and therefore is not subject to an upper limit caused by aliasing. A drawback of this formula (31) , however, is that it is highly sensitive to variations in the amplitude of the reflection signals and consequently is somewhat less reliable. The present invention therefore proposes that the two formulae (30) and (31) be combined with one another by first using formula (31) to estimate the velocity in a first approximation, an order of magnitude thus being obtained therefor, and by then using formula (30) for an estimate of the velocity in a second, more accurate approximation, the aliasing problem being eliminated since the order of magnitude is known.

The equivalent bandwidth BW_EQ of the reflection signal can be estimated by means of formula (32) , wherein Z₁ ≠ Z₂ ≠ Z₃ ≠ Z-_L, and to avoid noise effects it is better to chose T ≠ 0.

The normalized signal power S according to formulae (33) and (34), wherein 0≠T.,≠T₂≠0, ΔX=0, ΔY=0, and p₁₂=0, can only be estimated for the signals captured by the separate transducers. The normalized signal powers of the signals captured by the transducers 21_x and 21₂ are denoted by S_τ and S₂, respectively.

The beam constants in the two transverse directions C_Bx(Z_t) and C_By(Z_t) can be estimated by means of the formulae (35) and (36) , where S = /(S^) is the mean power of the signals captured by the two transducers. C_Bx(z_t) can be determined if the measuring setup is arranged in such a way that there are two transducers which are positioned in the XZ-plane, i.e. that ΔY=0, while their relative X-distance ΔX is known. This is therefore the measuring situation illustrated in figure 4. With the aid of formula (37) an estimate can then be calculated for the displacement S_x in the X-direction, from which the X-component V_x of the velocity vector can be calculated. C_By(z_t) can be determined if the measuring setup is arranged in such a way that there are two transducers which are positioned in the YZ-plane, i.e. that ΔX=0, while their relative Y-distance ΔY is known. With the aid of formula (38) an estimate can then be calculated for the displacement S_y in the Y-direction, from which the

Y-component V_y of the velocity vector can be calculated.

By combining V_x and V_y, the magnitude and direction of the transverse velocity vector V_uτ are then known.

The magnitude of S_z can be estimated with the aid of the formulae (30) and (31) , from which the axial velocity component then follows.

The above-described option assumes that at least three transducers are available, of which two are positioned in the XZ-plane and of which two are positioned in the YZ-plane. If only two transducers are available, only one of the beam constants can be determined, however. If the measuring setup is arranged in such a way that those transducers are positioned in the XZ-plane (the measuring situation illustrated in figure 4) , C_Bx can be estimated with the aid of formula (36) , and S_x can be estimated with the aid of formula (37) , wherein S = /(S₁S₂) is the mean power of the signals captured by the two transducers. C_By can then be calculated, however, if a fixed relationship is known between C_Bx and C_By. This is the case, for example, if the sound beams of those transducers are of known symmetry, so that it is possible to write C_Bx = eC_By, where e is a known or measureable form factor. If the measurements are carried out with transducers having a circular symmetrical sound beam, the relationship e = 1 applies. In other cases, a calibration curve should be measured for the missing beam constant, as will be obvious to those skilled in the art. Once C_By is known, S_y can be estimated with the aid of formula (39) . From S_x and S_y together, the horizontal and vertical components of the velocity vector follow, so that magnitude and direction of the transverse component of the velocity vector are known. The magnitude of S_z can be estimated with the aid of the formulae (30) and (31) , from which the axial velocity component then follows.

Obviously the same applies, mutatis mutandis, if the measuring setup is arranged in such a way that those transducers are positioned in the YZ-plane. If ΔX=0 and C_By(Z_t) has already been estimated, the movement in the Y-direction can be estimated by means of formula (38) , where S = /(S₁S₂) is the mean power of the signals captured by the two transducers .

If these measurements are carried out with only one transducer, all the formulae given above are simplified, as it is the case that p₁₂ = 0, ΔX = 0 , and ΔY = 0. Then, however, it is only possible to measure directly the central frequency, the axial velocity component, the bandwidth, the normalized signal power (and therefore the signal/noise ratio) and the variable C_Bx (Z_t) S_x ²+C_By (z_t) S_γ ² (formula 39) . The magnitude of S_z can be estimated with the aid of the formulae (30) and (31) , from which the axial velocity component then follows.

It is then not possible, however, to estimate a beam constant on the basis of the measured signals, and it is therefore not possible to determine S_x and S_y. Only if the measurement is carried out with a circular-symmetric sound beam is it possible, admittedly after a calibration procedure, to estimate the magnitude of the transverse velocity component. Since in that case the beam constant is the same for all transverse directions, it is the case that C_Bx(z_t) = C_By(z_t) = C_B(z_t), in which case the relationship C_Bx(z_t)S_x ²+C_By(z_t)S_Y ² = C_B(z_t) (S_x ² + S_γ ²) =

applies, where S^ is the movement in the XY-plane. The term C_B(z_t) can be determined by means of calibration. It should be noted, however, that using a single beam it is possible to determine the magnitude of the movement perpendicular to the beam axis, but not its direction in the plane perpendicular to the beam axis.

It will be obvious to a person skilled in the art that the scope of the present invention is not limited to the examples discussed hereinabove, but that various alterations and modifications are possible without deviating from the scope of the invention as defined in the appended claims.

Thus it is possible, for example, that instead of the said transducer, a combination of a separate transmitter and a separate detector is employed.

Moreover, the field of application of the present invention is not limited to the example discussed of blood flow in a blood vessel. Rather, the present invention can be used for remote detection of the velocity of any collection of reflection cores which behave as a more or less coherent entity: an example to be considered in this context is that of water droplets in a cloud.

Furthermore, a Hubert transformation can be performed in a manner different from the example described. Moreover, the definition of different data sets can be carried out either before or after the signal samples have been recorded, and either before or after the reflection signals have been rendered analytic.

FORMULES

RFio ) = j j j - y. z) g(x. y. *) x dy dz (1)

RF_2ιT(t) = J J f f,(x. y, z) g(x - (p„T + p₁₂) S_x - ΔX, y - (_PllT + p₁₂) S_γ - ΔY, z - (_PllT + p₁₂) S₂) dx dy dz (2)

-co-oo-oo

R(T, Z) = /RF₁-₀(t)RF_{2 T} t + (3)

R(T,Z) = J X J J J ft'(xN. ^z) f ₊ z (x^,.y',z^,)(g(x,y,z) g(x'-(p_πT + _Pl2) S_x - ΔX,y'-(p„T + p₁₂) S_Y - ΔY,z'-(_Pl1T _{+ Pl2}) s_z))dx dy dz dx' dy' dz¹ (4

(g(x, y, z) g(x'-(p_uT + p₁₂) S, - ΔX, y'-(p_nT + p₁₂) S_Y - ΛY, z'-(_Pl1T + p₁₂) S_z)) =

(5

(g²) δ(x - x'-t (p„T H p₁₂) S_x + ΔX, y - y' . (p„T ^ p₁₂) S_y + ΛY, z - z¹ , (p„T -.- p,₂) S_z)

R(T,Z) = (g²) f f f f_l ^*(x₁y,2)f_t (x + (p_uT + p„)S_x+ΔX,y + (_Pl1T + p,₂)S_γ+ΔY₁z + (p₁₁T + p₁₂)S_z)dxdydz (6)

)π

BW ^v _cε_nα = , - (8)

^" dBkϊioi ^{B dB}

RF(t) = V2 BW_EQ exp[- 2π BW|_Q (t - t ] expjj (2π f_c (t - t_d) + ψ)] (9)

8π BW_Q 4πf_c f,(z) - exp[ C_z (z, - z)²]exp[j(c_p (z, - z) + ψ)] ; C₂ = C_P = (10)

f,(y) = exp[c_Y(_Zl) y²] ; C_γ(z.) = π^ (12

5 wy5_B(z,)

f_t(x,y,z) = exp[c_x(z_t) x² + C_γ(z_t) y² + C_z (z, - z)²] exp[j JC_P (z, - z) + _Ψ)] (13

IΛ O 00 OT σ

CM

II

ι^ OO O

CM CM CM

CM

CM CO en O CO O O

1 ra

E

M

<C 0 O CO CO

oo C33 co CO

Claims

1. Method of remotely detecting the distribution of the vectorial velocities in a diffuse medium (13) , comprising the following steps:

(a) emitting a pulse-shaped sound signal (22) towards the medium (13) by means of at least one suitable transmitter (21) ;

(b) detecting a reflection sound signal (23) caused by the medium (13) by means of at least one suitable detector (21) ; (c) converting the reflection sound signal (23) detected by the detector (21) into an electrical signal (24), the strength of said electrical signal (24) being at all times representative for the instantaneous sound intensity; (d) sampling the electrical signal (24) to provide successive signal samples RF ( ╧ä , ╬╢ ) , wherein ╬╢ denotes a serial number for the successive signal samples, and wherein T is an integer denoting a serial number of the signal captured for processing; (e) recording the signal samples RF { ╧ä , ╬╢) in a memory (31) ;

(f) repeating the steps (a) to (e) inclusive for successive values of T, the signal samples RF(╧ä,f) being grouped into two data sets, denoted as RF₁(╧ä,(^") and RF₂ ( T , ╬╢) , respectively;

(g) rendering analytic the reflection signals RF₁(╧ä,f) and RF₂(╧ä,f) thus obtained by means of e.g. a Hubert transformation to obtain analytic measured signals RF_1A(╧ä, f) and RF_2A(╧ä,f) ! (h) defining mutually identical data windows w₁(╧ä,╧Ç and w₂(╧ä,f) in each data set RF_1A(╧ä,H and RF_2A(╧ä,f) of the analytic reflection signals;

(i) calculating correlation coefficients R(T,Z) between these two data windows v╬╣₁ ( r , ╬╢ ) and w₂(╧ä,f) ; (j) calculating velocity components from the correlation coefficients thus calculated; wherein the axial component of the vectorial velocity in the medium (13) is calculated with the aid of formula (31) on the basis of the magnitudes of the correlation coefficients.

2. Method according to claim 1, wherein a first reflection measurement is carried out with the aid of a first transmitter/detector combination (21_x) to obtain a first set of measured points w₁(╧ä,f), wherein a second reflection measurement is carried out with the aid of a second transmitter/detector combination (21₂) to obtain a second set of measured points vr₂ ( ╧ä , ╬╢) , ^an< wherein the correlation is carried out between, on the one hand, the first measuring points obtained by means of the first transmitter/detector combination (21_x) and, on the other hand, the second measured points obtained by means of the second transmitter/detector combination (21₂) .

3. Method according to claim 2, wherein the two said sets of measured points are obtained by alternately operating the said first and second transmitter/detector combinations (21_1#- 21₂) .

4. Method according to claim 1, wherein a first reflection measurement is performed with the aid of one transmitter/detector combination (21) to obtain a first set of measured points vj₁ { ╧ä , ╬╢) , wherein a second reflection measurement is performed with the aid of the same transmitter/detector combination (21) to obtain a second set of measured points w₂(╧ä,f)_/ and wherein the correlation is carried out between, on the one hand, the first measured points obtained by means of the one transmitter/detector combination (21) and, on the other hand, the second measured points obtained by means of the same transmitter/detector combination (21) .

5. Method according to claim 4, wherein the two said sets of measured points are obtained by carrying out measurements alternately for the first and for the second set by means of the said one transmitter/detector combination (21) .

6. Method according to any one of the preceding claims, wherein first the order of magnitude of the axial component of the vectorial velocity of the medium (13) is calculated with the aid of formula (31) on the basis of the magnitude of the correlation coefficients, and then a more accurate estimate of the axial component of the vectorial velocity of the medium (13) is calculated with the aid of formula (30) on the basis of the arguments of the correlation coefficients, taking into account the order of magnitude calculated with the aid of formula (31) to eliminate aliasing effects.

7. Method according to any one of the preceding claims, wherein use is made of a single sound beam of circular symmetric shape, so that the relationship C_Bx(z_t) = C_By(z_t) = C_B(z_t) applies for the beam constants, wherein the form function C_B(z_t) is determined by means of calibration; and wherein the variable C_B(z_t)S_LAT ² = C_Bx (z_t) S_x ²4-C_By (z_t) S_╬│ ² = C_B(z_t) (S_x ² 4- S_y ²) , wherein S^ is the motion in the XY-plane, is calculated with the aid of formula (39) in order to thus provide a transverse component of the vectorial velocity in the medium (13) .

8. Method according to any one of claims 1-6, wherein use is made of two transducers (21₁ and 21₂) which define an XZ-plane; wherein an X-component of the vectorial velocity in the medium (13) is calculated with the aid of formula (37) ; wherein M is defined by formula (35) ; wherein S = (S₁S₂) ; wherein S₁ and S₂ are defined by the formula (34) ; and wherein K is defined by formula (33) and (32) .

9. Method according to claim 8, wherein the shape of the sound beams (22) emitted by the transmitters (21) has an XY-symmetry, so that the relationship C_Bx = eC_By applies, where e is a form factor which, for example, is equal to 1 if the sound beam (22) emitted by the transmitter (21) is circular symmetric; wherein C_Bx is calculated with the aid of formula (36) ; wherein C_By is calculated from the relationship C_Bx = eC_By; and wherein a Y-component of the vectorial velocity in the medium (13) is calculated with the aid of formula (39) .

10. Method according to claim 8, wherein a relationship for C_By is determined by means of calibration; wherein C_Bx is calculated with the aid of formula (36) ; and wherein a Y-component of the vectorial velocity in the medium (13) is calculated with the aid of formula (39) .

11. Method according to any one of claims 1-6, wherein use is made of three or more transducers, wherein two of said transducers (2^ and 21₂) define an XZ-plane, and wherein two of said transducers define a YZ-plane; wherein an X-component of the vectorial velocity in the medium (13) is calculated with the aid of formula (37) ; wherein M is defined by formula (35) ; wherein S = (S₁S₂) ; wherein S_x and S₂ are defined by formula (34) ; and wherein K is defined by formula (33) and (32) ; and wherein the Y-component of the vectorial velocity in the medium (13) is calculated with the aid of formula (38) .

12. Method according to claim 11, wherein the beam ccoonn;stants C_Bx and C_By are calculated with the aid of formula (36:

13. Method according to any one of the preceding claims, wherein in step (g) the modification of the frequency spectrum is described by formula (27) .

14. Method according to any one of the preceding claims, wherein the calculation of the correlation function R(T,Z) is described by formula (28).

15. Method according to any one of the preceding claims, wherein the central frequency is calculated in accordance with formula (29) .

16. Method according to any one of the preceding claims, wherein the equivalent bandwidth is calculated in accordance with formula (32) .

17. Method according to any one of the preceding claims, wherein the transmitter and detector are combined into a single transducer.