A three-dimensional ultrasonic vector Doppler imaging method based on row-column addressing array

Abstract

The present invention belongs to the technical field of ultrasonic Doppler blood flow imaging, specifically a three-dimensional ultrasonic vector Doppler imaging method based on a row-column addressing array. The present invention combines a high-frame-rate two-dimensional non-focused wave transmission sequence and a row-column addressing array (RCA) of a least square multi-angle Doppler velocity estimator to achieve three-dimensional ultrasonic vector Doppler imaging; specifically, it includes: acquiring radio frequency data from each tilted non-focused wave transmission, and then performing IQ demodulation, beamforming, and clutter filtering; performing high-frame-rate multi-angle vector Doppler analysis and velocity estimation on the results through RCA; and performing dual-wavelength anti-aliasing processing. The present invention can achieve an imaging speed of more than one thousand frames, significantly improve the sensitivity to low-speed flow and the measurement accuracy of high-speed flow, and has been verified by clinical environment experiments; it is suitable for clinical applications, especially in the dynamic monitoring of blood flow in large blood vessels, and can provide important support for hemodynamic research and clinical diagnosis.

Description

Technical Field

[0001] The invention belongs to the technical field of ultrasonic Doppler blood flow imaging, and in particular relates to a three-dimensional ultrasonic vector Doppler imaging method. Background Art

[0002] A large number of cardiovascular and cerebrovascular diseases are associated with hemodynamic abnormalities. Therefore, three-dimensional blood flow visualization technology with high spatiotemporal resolution and a wide field of view is extremely important. Doppler ultrasound is widely used in clinical practice as one of the primary tools for noninvasively estimating blood flow velocity. Ultrasound vector Doppler combines multi-angle, high-frame-rate Doppler with the least squares principle, improving the signal-to-noise ratio while suppressing estimation variance. This technology has been widely used to measure hemodynamics in large vessels such as the carotid artery and to image small vessels with diameters of several hundred microns.

[0003] With the development of technology, three-dimensional imaging technology has gradually been introduced into blood flow measurement. Generally speaking, two-dimensional matrix array imaging requires thousands of electronic channels to control a single element, and the number of channels increases quadratically with the number of elements in the array side length, which shows that large-aperture arrays have high requirements for multi-channel hardware. In addition, high-density element arrangement remains a challenge in the manufacture of two-dimensional probes. The emergence of row-column addressing arrays provides an economical and effective solution for three-dimensional high-frame-rate blood flow imaging. By addressing rows and columns separately, the number of channels can be effectively reduced, and the reduction in resolution and contrast caused by the small number of channels can be compensated by improved composite strategies. At present, there are no reports on methods for estimating in vivo blood flow velocity using row-column addressing arrays. Summary of the Invention

[0004] The object of the present invention is to provide a three-dimensional ultrasonic vector Doppler imaging method based on a row-column addressing array.

[0005] The present invention provides a three-dimensional ultrasonic vector Doppler imaging method based on a row-column addressable array (RCA), which is based on a row-column addressable ultrasonic array (RCA), combined with a high-frame-rate orthogonal two-dimensional non-focused wave transmission sequence and a multi-angle three-dimensional Doppler velocity estimation algorithm to achieve three-dimensional high-frame-rate ultrasonic imaging and three-dimensional blood flow vector Doppler imaging.

[0006] The row-column addressed ultrasound array (RCA) is a two-dimensional probe array, specifically composed of two orthogonal long linear array elements. It transmits alternately in rows and columns; the received signal is a combination of row and column echo signals. By decomposing the two-dimensional area array into two long linear arrays, interleaved transmission of row and column elements is achieved, resulting in composite acquisition of volumetric images.

[0007] The high-frame-rate orthogonal 2D unfocused wave transmission sequence refers to an RCA array-based system that transmits orthogonal multi-angle tilted unfocused waves at a rate of more than 20 frames per second (frame rate > 20 frames / s), enabling 3D ultrasound imaging at over 20 frames per second. Specifically, the row array transmits unfocused waves tilted along the x-direction, while the column array transmits unfocused waves tilted along the y-direction.

[0008] The multi-angle 3D Doppler velocity estimation algorithm is designed to address the correlation between frequency offset and angle in Doppler measurements. By comparing the Doppler frequency offsets of echo signals measured at different angles, it achieves 3D frequency offset extraction and 3D vectorized blood flow velocity imaging of the echo signal. This includes high-frame-rate multi-angle vector Doppler analysis and velocity estimation based on the RCA array, as well as dual-wavelength anti-aliasing processing.

[0009] The three-dimensional high-frame-rate ultrasound imaging refers to the completion of spatial positioning of the measured object in the row and column directions through the delayed superposition beam synthesis step based on the size of ultrasound data received at different times, and the realization of three-dimensional ultrasound imaging with an imaging frame rate greater than 20 frames / s.

[0010] The three-dimensional blood flow vector Doppler imaging refers to imaging the speed and direction of blood flow by using the delayed superposition beam synthesis step and the multi-angle three-dimensional Doppler velocity estimation algorithm based on the size of ultrasonic data received at different times.

[0011] The present invention provides a three-dimensional ultrasonic vector Doppler imaging method based on a row-column addressing array, which specifically comprises the following steps:

[0012] Step S1, acquiring radio frequency data from each tilted non-focused wave transmission, and then performing IQ demodulation, beamforming, and clutter filtering;

[0013] Step S2, performing high frame rate multi-angle vector Doppler analysis and velocity estimation on the results through the RCA array;

[0014] Step S3: dual-wavelength anti-aliasing processing.

[0015] Further:

[0016] In step S1, the clutter filtering can be specifically SVD-based clutter filtering, which is implemented by using a spatiotemporal filter based on singular value decomposition to extract dynamic blood signals and applying them to each steering angle respectively. Specifically, for each beamformed volume, the complex signals of the beamformation corresponding to each steering angle are integrated into a 4D matrix of dimension (Nx, Ny, Nz, M), where M is the set length of the velocity estimate. The matrix is ​​then reconstructed into a two-dimensional Casorati matrix of dimension (Nx×Ny×Nz, M), and the five largest of the 50 singular vectors are removed. Each comprehensive estimation performs a total of 2*K SVD filtering, where K is the number of steering angles.

[0017] In step S2, the high frame rate multi-angle vector Doppler analysis is specifically: Doppler velocity V θi The lagged first-order autocorrelation phase algorithm is used for estimation. The calculation formula is:

[0018]

[0019] Among them, PRF eff is the effective composite frame frequency, f0 is the center frequency, c is the speed of sound, * represents the conjugate operation, {S IQf (1), S IQf (2), ..., S IQf (M)} is the filtered complex signal with a set length of M, R{n} is its autocorrelation value of the nth order lag, usually 0<n≤5, and arg(-) represents the parameter of the complex number. In principle, the derived value in the above formula represents the velocity along the axis of the sound beam, which satisfies the following relationship:

[0020]

[0021] Where |v| is the actual velocity, α is the angle between the flow direction and the normal of the probe surface (also called the angle between the beam and the flow), and θ and are the transmit beam angle and the receive beam angle, respectively.

[0022] In the present invention, only the unfocused transmission wave tilted in one horizontal direction, i.e., tilted in the x or y direction, is considered. Considering the transmission beam angle θ tilted in the x direction, the Doppler velocity V θ It can be obtained by the following formula:

[0023]

[0024] Receiving beam angle It is usually set to 0, so it can be simplified to:

[0025] 2V θ =v z (cosθ+1)+v xsinθ. (5)

[0026] Two separate velocity components v x and v z It can be solved by constructing a linear equation system that is derived from at least two steering angles. To improve the estimation accuracy, more than two transmit beam angles are used to construct an overdetermined equation system to solve v x and v z :

[0027]

[0028] Solving the velocity matrix can be viewed as a least squares fitting problem, where the mean square error of the estimated value is minimized. Rewriting the matrix symbol in the above equation to Av=u, the flow vector v can be solved by pseudo-inverse operation:

[0029]

[0030] The pseudo-inverse term is a constant matrix that can be pre-calculated. y and v z The process is the same as above. That is, by emitting non-focused waves along the corresponding transverse inclination, the two transverse velocity components v are solved separately. x and v y .

[0031] In step S2, the RCA-based velocity estimation is as follows: Since the row array can only transmit non-focused waves tilted along the x direction, the row transmission and column reception (RC) can only estimate v x and v z Similarly, the column transmit and row receive (CR) can only estimate v y and v z In order to reduce v x and v y The time difference between the estimated values, RC and CR, is exchanged after each tilted unfocused wave is transmitted. The sequence consists of two types of blocks: UVD (ultrasound vector Doppler) blocks and B-Mode Ultrasound (B-Mode Ultrasound) blocks.

[0032] For UVD blocks, each UVD block contains 2K transmissions, where K represents the number of tilt angles. The UVD block is executed M times in succession to mitigate the bias in velocity estimation, i.e., M is the set length. Therefore, the volume rate is determined by PRF / (2K). K is usually small (K≤5) to ensure a high volume rate. Using the autocorrelation method, each Doppler velocity v is estimated separately after a tilted unfocused wave transmission. iHowever, there is no dynamic focusing along the emission direction in RC or CR without multi-angle non-focused wave recombination. This indicates that v i Different targets on the emission plane, i.e., the xz plane corresponding to RC and the yz plane corresponding to CR, cannot be distinguished. The transverse velocity is then solved using the least squares method by combining all Doppler velocities of the unfocused waves at different tilts. The RC and CR waves are then composited using an orthogonal unfocused wave composite scheme to obtain the power Doppler, which serves as a mask for region selection.

[0033] For the B-mode block, which contains more steering angles, it is executed after the UVD block. In order to obtain high-quality B-mode images, the orthogonal non-focused wave recombination technology is also used. Since the B-mode block is only executed once in a velocity estimation cycle, its impact on the volume rate is almost negligible. In addition, when the set length M is too short, it tends to complicate clutter filtering and make the velocity estimation more susceptible to noise interference. Therefore, the blocking rate PRF / (2KM) of the velocity estimation is severely limited by the set length M, making it difficult for the UVD method to capture the rapid changes of high-speed blood flow. Therefore, in order to obtain more velocity details that change over time, a sliding window method is used on the UVD block. The set length M remains unchanged, but the set units overlap with each other, and the sliding step size is L.

[0034] In step S3, the dual wavelength anti-aliasing process is specifically as follows: Due to the limited flow velocity, aliasing will inevitably occur when estimating high flow velocities. The absolute value of the maximum velocity that can be measured without aliasing (called the Nyquist velocity (V N ))By PRF eff and the center frequency f0, as shown below:

[0035]

[0036] That is, the measurable Doppler velocity V D Affected by [-V N :V N ] constraint. Alias-free Doppler velocity With V D The relational expression is as follows:

[0037]

[0038] Among them, n N is the Nyquist number (n N ∈Z), only when n N = 0, there is no aliasing phenomenon in Doppler velocity. N It can be expressed by the following formula:

[0039]

[0040] Where floor(-) represents the largest integer less than or equal to. It is worth noting that if the corresponding Nyquist number is known, the Doppler velocity can be de-aliased. According to the dual-wavelength strategy, the de-aliased velocity V D1 and V D2 The corresponding Nyquist velocity V N1 and V N2 The selection must follow the following rules:

[0041]

[0042] Where p is a positive integer and the Nyquist velocity is expanded to pV N1 Combining the above equations, we can get the following equation:

[0043]

[0044] in, n N1 and n N2 are all unknowns, so it is still impossible to solve the equation for a clear speed. However, since n N1 and n N2 are all integers, so formula (11) can be rewritten as:

[0045]

[0046] Among them, nint(-) is the nearest integer. Nyquist number n N1 and n N2 At the extended maximum Nyquist velocity pV N1 or pV N2 As a bound, we can obtain the following additional constraints:

[0047]

[0048] The above constraints ensure that the solution is unique within the defined range. A lookup table can be constructed to solve for the alias-free velocity. Two bandpass filters are performed on the raw data to produce the Doppler velocity V D1 and V D2 The center frequencies f1 and f2 of the two bandpass filters are within the bandwidth of the ultrasound probe. N =cPRF / (4f0), the selection of f1 and f2 should satisfy the following relationship:

[0049]

[0050] Theoretically, the extended Nyquist speed depends only on the value of p. However, when the actual speed exceeds twice the Nyquist speed, especially when the signal-to-noise ratio is low, the method often deviates significantly. Therefore, in general practice, p is usually chosen to be 3 or 4.

[0051] The technical features and functional advantages of the present invention mainly include:

[0052] The present invention utilizes an RCA array and high-frame-rate non-focused wave emission to achieve an imaging speed of more than one thousand frames, significantly improving the sensitivity to low-speed flow and the measurement accuracy of high-speed flow. Through multi-angle vector Doppler technology, angle-independent flow velocity measurement can be achieved, enhancing the visualization capability of complex flow patterns. This method can be applied to 4D flow velocity imaging and in vivo blood flow imaging based on the RCA array, and can also be extended to pulse wave analysis. The present invention provides a cost-effective and real-time alternative for blood flow measurement, which is suitable for clinical applications, especially in the dynamic monitoring of blood flow in large blood vessels (such as the carotid artery). Preliminary clinical environment experiments have verified that this method has the potential to perform flow velocity analysis on large blood vessels such as the human common carotid artery, has significant technical advantages, and can provide important support for hemodynamic research and clinical diagnosis. BRIEF DESCRIPTION OF THE DRAWINGS

[0053] Figure 1 is a sequence diagram. Among them, (a) is the UVD transmission sequence diagram based on RCA. RC and CR are used to estimate v x 、v z and v y 、v z . Each individual unfocused wave transmission is automatically correlated to obtain a Doppler velocity that is not focused in the transmission plane. The least squares method is used to compensate for the unfocusing effect and solve for the lateral velocity. All M UVD blocks are compounded to obtain the power Doppler. After the M UVD blocks, a B-mode block with a steering angle of H is performed to obtain a high-quality B-mode image. (b) shows the sliding window method, where the set length M remains unchanged, but the set units overlap with each other, and the sliding step size is L.

[0054] Figure 2 This figure shows a 3D velocity vector diagram and a comparison of the estimated velocity with the actual velocity. The estimated results were obtained at α = 75°, β = 60°, and a PRF of 40 kHz. The velocity distribution is consistent with the established laminar flow pattern, and the root mean square error of the estimated velocity is less than 10%.

[0055] Figure 3The results of an in vitro flow model with a flow velocity of 24.5 cm / s are shown in the upper panel. The 3D flow vector diagram and 2D velocity profile are shown, while the lower panel depicts the parabolic profiles for all three sections. The blue dashed line is the true velocity, the red line is the estimated mean velocity profile, and the gray shaded area represents an SD value. The estimated mean deviation is shown in the figure. and

[0056] Figure 4 These are estimates of the carotid artery in a healthy volunteer. The top figure shows the estimated mean velocity, maximum velocity, and pulse wave, all acquired simultaneously. The lower inset shows 3D velocity maps acquired at four different time points. DETAILED DESCRIPTION

[0057] The present invention is further described below through embodiments with reference to the accompanying drawings.

[0058] Example

[0059] A customized model was used. A rubber tube with an inner diameter of 4 mm and a wall thickness of 0.5 mm was fixed in a 3D printed box filled with water, and a blood-simulating fluid (KS205D, DELTA) driven by a micro hydraulic pump system (LSP01-2B, LONGER) circulated in the rubber tube. The position of the transducer relative to the rubber tube, especially the two rotation angles: α and β, was adjusted by a customized fixture. The flow rate was set to 24.5 cm / s and controlled by the pumping system. α is between 60°, 70° and 90°, so the maximum Nyquist frequency is 1.91 kHz. The two sides of the container are covered with sponge absorbent material to reduce boundary reflections. In terms of emission, 5 angles evenly distributed in the range of [-7.5°, 7.5°] are alternately emitted by the row and column array, and the PRF is 10 kHz. Therefore, the effective PRF eff = 10kHz / (5*2) = 1kHz. To reduce the interference of noise on speed estimation, the set length is set to 50.

[0060] In this embodiment, the carotid artery of a healthy volunteer is measured. The number of steering angles is set to 5, and the angle span is 15°. The RF data acquisition process lasts for 6 seconds, and the PRF eff The ultrasonic sequence spatial-peak-time-average (ISPTA.3) and mechanical index (MI) meet the recommendations of the U.S. Food and Drug Administration (FDA). At 20V, MI and ISPTA.3 are 0.42 and 127.9mW / cm, respectively. 2The soft tissue thermal index (TIS) is 2.8 and the acquisition lasts for 6 seconds. The set length of velocity estimation is set to 100 and the sliding step is 20. Therefore, the refresh rate of the estimated velocity is 50 Hz. The dynamic blood signal is extracted using a spatiotemporal filter based on SVD. The complex signal of beamforming is reshaped into a dimension of (N x ×N y ×N z ,N t ) and discard the largest 15 of the 100 singular vectors. Pulse waves were collected simultaneously using a gated pulse sensor (GuoJun, HK-2000C).

[0061] Figure 2 The figure shows a 3D velocity vector diagram and a comparison of the estimated velocity results with the actual velocity for the corresponding embodiment. The estimated results were obtained at α = 75°, β = 60°, and a PRF of 40 kHz. The velocity distribution is consistent with the pre-set laminar flow pattern, and the root mean square error of the estimated velocity is less than 10%.

[0062] Figure 3 The following are the results of the in vitro flow model with a flow velocity of 24.5 cm / s in the example. The upper column shows the three-dimensional flow vector diagram and the two-dimensional velocity profile, while the lower figure describes the parabolic profile of all three parts. The blue dashed line is the true velocity, the red line is the estimated mean velocity profile, and the gray shaded area represents an SD value. The estimated mean deviation is shown. and

[0063] Figure 4 The figure below shows the estimated carotid artery velocity for a healthy volunteer in the example. The top figure shows the estimated mean velocity, maximum velocity, and pulse wave, all acquired simultaneously. The lower part of the inset shows 3D velocity maps acquired at four different time points.

Claims

1. A three-dimensional ultrasonic vector Doppler imaging method based on a row-column addressing array, characterized in that: Using a row-column addressed ultrasound array (RCA), combined with a high-frame-rate orthogonal two-dimensional non-focused wave transmission sequence and a multi-angle three-dimensional Doppler velocity estimation algorithm, three-dimensional high-frame-rate ultrasound imaging and three-dimensional blood flow vector Doppler imaging are achieved; The row-column addressable ultrasonic array is a two-dimensional probe array, specifically an ultrasonic array composed of two orthogonal long array elements, which adopts a row-column alternating transmission mode; the received signal is a combination of row and column echo signals; by decomposing the two-dimensional area array into two long linear arrays, the row and column elements are staggered to transmit, and a composite volumetric image is collected; The high frame rate orthogonal two-dimensional unfocused wave transmission sequence is based on an RCA array to transmit orthogonal multi-angle tilted unfocused waves at more than 20 frames per second, thereby achieving three-dimensional high frame rate ultrasound imaging at more than 20 frames per second. The transmission of the orthogonal multi-angle tilted unfocused waves is specifically as follows: the row array transmits unfocused waves tilted along the x-direction, and the column array transmits unfocused waves tilted along the y-direction. The multi-angle 3D Doppler velocity estimation algorithm compares the Doppler frequency offsets of echo signals measured at different angles to extract the 3D frequency offset of the echo signal and perform 3D vectorized blood flow velocity imaging. This includes high-frame-rate multi-angle vector Doppler analysis and velocity estimation based on the RCA array, as well as dual-wavelength anti-aliasing processing. The three-dimensional high-frame-rate ultrasound imaging is to complete the spatial positioning of the measured object in the row and column directions through the delay-superposition beam synthesis step based on the size of the ultrasound data received at different times, and realize three-dimensional high-frame-rate ultrasound imaging; The three-dimensional blood flow vector Doppler imaging is to image the speed and direction of blood flow by using the delay-superposition beam synthesis step and the multi-angle three-dimensional Doppler velocity estimation algorithm based on the size of the ultrasonic data received at different times. The specific steps of the imaging method are: Step S1, acquiring radio frequency data from each tilted non-focused wave transmission, and then performing IQ demodulation, beamforming, and clutter filtering; Step S2, performing high frame rate multi-angle vector Doppler analysis and velocity estimation on the results through the RCA array; Step S3, dual-wavelength anti-aliasing processing; The high frame rate multi-angle vector Doppler analysis described in step S2 is specifically: Doppler velocity V θi The lagged first-order autocorrelation phase algorithm is used for estimation, and the calculation formula is: Among them, PRF eff is the effective composite frame frequency, f0 is the center frequency, c is the speed of sound, {S IQf (1), S IQf (2), ..., S IQf (M)} is a filtered complex signal with a set length of M, * represents the conjugate operation, R{n} is its autocorrelation value of the nth order lag, 0<n≤5, arg(-) represents the parameter of the complex number; the derived value in the above formula represents the velocity along the axis of the sound beam, which satisfies the following relationship: Where |v| is the actual velocity, α is the angle between the flow direction and the normal of the probe surface, also known as the angle between the beam and the flow, and θ and are the transmit beam angle and the receive beam angle respectively; Consider a non-focused wave transmitted in a transverse direction, that is, tilted in the x or y direction; the transmit beam angle θ tilted in the x direction, the Doppler velocity V θ It is obtained by the following formula: Receiving beam angle It is usually set to 0, so Equation (4) is simplified to: 2V θ =v z (cosθ+1)+v x sinθ (5) Two separate velocity components v x and v z By constructing a linear equation system to solve, specifically using more than two transmit beam angles to construct an overdetermined equation system to solve v x and v z : Solving the velocity matrix can be viewed as a least squares fitting problem, where the mean square error of the estimated value is minimized. Rewrite the matrix symbol in the above equation as Av=u, and solve the flow vector v through pseudo-inverse operation: The pseudo-inverse term is a constant matrix that can be pre-calculated; v is solved by emitting an unfocused wave tilted in the y direction. y and v z The process is the same as above, that is, by emitting non-focused waves along the corresponding transverse inclination, the two transverse velocity components v are solved respectively. x and v y ; In step S2, the RCA-based velocity estimation is as follows: since the row array can only transmit non-focused waves tilted along the x direction, the row transmitting and column receiving RC only estimates v x and v z Similarly, the column transmit and row receive CR only estimates v y and v z ; In order to reduce v x and v y The time difference between the estimated values, RC and CR are interchanged after transmitting each tilted unfocused wave; the sequence consists of two types of blocks: UVD blocks and B-mode blocks; For the UVD block, each UVD block contains 2K transmissions, where K represents the number of tilt angles; the UVD block is executed M times in succession to reduce the deviation of velocity estimation, that is, M is the set length; the volume rate is determined by PRF / (2K), and K takes a small value to ensure a high volume rate; the auto-correlation method is used to estimate each Doppler velocity v after a tilted unfocused wave transmission. i However, there is no dynamic focusing along the emission direction in RC or CR without multi-angle non-focused wave recombination; this indicates that v i It is impossible to distinguish different targets on the emission plane, i.e., the xz plane corresponding to RC and the yz plane corresponding to CR. Subsequently, the lateral velocity is solved by combining all Doppler velocities of unfocused waves with different tilts using the least squares method. The RC and CR are composited using an orthogonal unfocused wave composite scheme to obtain the power Doppler, which is used as a mask for region selection. For the B-mode block, it contains more steering angles and is executed after the UVD block. In order to obtain high-quality B-mode images, the orthogonal non-focused wave recombination technology is also used. Since the B-mode block is only executed once in a velocity estimation cycle, its impact on the volume rate is almost negligible. In addition, when the set length M is too short, it will complicate the clutter filtering and make the velocity estimation more susceptible to noise interference. Therefore, the blocking rate PRF / (2KM) of the velocity estimation is severely limited by the set length M, making it difficult for the UVD method to capture the rapid changes of high-speed blood flow. In order to obtain more velocity details that change over time, the sliding window method is used on the UVD block. The set length M remains unchanged, but the set units overlap with each other, and the sliding step size is L. The dual-wavelength anti-aliasing process in step S3 is as follows: due to the limited flow velocity, aliasing will inevitably occur when estimating high flow velocities; the absolute value of the maximum velocity that can be measured without aliasing is called the Nyquist velocity V N , by PRF eff and the center frequency f0, as shown below: The measurable Doppler velocity V D Affected by [-V N :V N ] constraint; alias-free Doppler velocity With V D The relationship is as follows: Among them, n N is the Nyquist number, n N ∈Z, only when n N = 0, there is no aliasing phenomenon in Doppler velocity; n N It is expressed by the following formula: Where floor(-) represents the largest integer less than or equal to ; if the corresponding Nyquist number is known, the Doppler velocity can be de-aliased; according to the dual-wavelength strategy, from the two de-aliased velocity V D1 and V D2 The Nyquist number is obtained; the corresponding Nyquist velocity V N1 and V N2 The selection follows the following rules: Where p is a positive integer and the Nyquist velocity is expanded to pV N1 ; Combining the above formulas, we get the following equation: in, n N1 and n N2 are all unknown; since n N1 and n N2 are all integers, so formula (11) can be rewritten as: Among them, nint(-) is the nearest integer; Nyquist number n N1 and n N2 At the extended maximum Nyquist velocity pV N1 or pV N2 As a bound, we get the following additional constraints: The above constraints ensure that the solution is unique within the defined range; a lookup table is constructed to solve for the alias-free velocity; two bandpass filters are performed on the raw data to generate the Doppler velocity V D1 and V D2 The center frequencies f1 and f2 of the two bandpass filters are within the bandwidth of the ultrasound probe; according to V N =cPRF / (4f0), the selection of f1 and f2 satisfies the relationship: The extended Nyquist speed depends only on the value of p, which can be chosen as 3 or 4.

2. The three-dimensional ultrasonic vector Doppler imaging method according to claim 1, characterized in that: In step S1, the clutter filtering uses a spatiotemporal filter based on singular value decomposition to extract the dynamic blood signal and applies it to each steering angle separately. Specifically, for each beamformed volume, the complex beamformed signal corresponding to each steering angle is integrated into a 4D matrix of dimension (Nx, Ny, Nz, M), where M is the set length of the velocity estimate. The matrix is ​​then reconstructed into a two-dimensional Casorati matrix of dimension (Nx×Ny×Nz, M), and the five largest of the 50 singular vectors are removed. Each comprehensive estimation performs 2*K SVD filtering times, where K is the number of steering angles.

Images