Acoustic velocity non-uniform spiral wave spectrum passive cavitation imaging method, device, equipment and medium

Abstract

The invention provides a sound velocity non-uniform spiral wave spectrum passive cavitation imaging method, device and equipment and a medium, and the method comprises the steps: a convex array ultrasonic transducer passively receives an ultrasonic time domain signal p (r, phi, t) in the sound velocity non-uniform medium, array elements of the convex array ultrasonic transducer are distributed on an arc with a radius R, r is the extrapolation radius of a current layer, t is the extrapolation radius of the current layer, and t is the extrapolation radius of the current layer; phi is a circumferential angle in the coordinate system, t is time, and R is a distribution surface where an array element of the convex array ultrasonic transducer is located; determining sound velocity distribution of a propagation medium based on the ultrasonic time domain signal p (r, phi, t), and aligning the sound velocity distribution with a probe of the convex array ultrasonic transducer to generate an aligned sound velocity diagram; the ultrasonic time domain signal p (r, phi, t) is subjected to signal conversion in the time t and the circumferential direction phi, a time frequency domain wave spectrum # imgabs0 # and a spatial frequency domain wave spectrum Pn (r, omega) are obtained, omega represents frequency, and n is spatial frequency corresponding to Fourier transform with phi as an independent variable; and on the basis of the aligned sound velocity graph, selecting sound field intensity in a specific frequency band omega, and performing superposition to obtain a passive cavitation image I (r, phi) under a coordinate system.

Description

Technical Field

[0001] The present invention belongs to the technical field of ultrasonic imaging, and in particular relates to a method, device, equipment and medium for passive cavitation imaging of sound velocity non-uniform spiral wave spectrum. Background Art

[0002] Focused ultrasound therapy technology based on the acoustic cavitation effect (ultrasound cavitation therapy) can improve the permeability of biological barriers, deliver therapeutic genes or drugs to lesions, produce the desired biological effects to achieve the purpose of treatment, and has important application prospects in tumor diseases, central nervous system diseases, etc. However, acoustic cavitation activity is affected by the dynamic changes of factors such as the acoustic pressure and cavitation threshold in the treatment area, and is difficult to predict and control. Therefore, there is a risk of off-target treatment in ultrasonic cavitation therapy, and efficient methods are needed to monitor microbubble cavitation in real time and provide precise positioning to ensure the safety and reliability of the ultrasonic cavitation therapy process. Summary of the invention

[0003] In view of the above-mentioned shortcomings of the prior art, the object of the present invention is to provide a method, device, equipment and medium for passive cavitation imaging of sound velocity non-uniform spiral wave spectrum to solve the above-mentioned problems.

[0004] The present invention provides a method for passive cavitation imaging of sound velocity non-uniform spiral wave spectrum, comprising:

[0005] A convex array ultrasonic transducer passively receives an ultrasonic time domain signal p(r,φ,t) in a medium with non-uniform sound velocity, wherein the array elements of the convex array ultrasonic transducer are distributed on an arc of radius R, wherein r is the extrapolated radius of the current layer, φ is the circumferential angle in the coordinate system, t is time, and R is the distribution surface of the array elements of the convex array ultrasonic transducer;

[0006] Based on the ultrasonic time domain signal p(r,φ,t), determine the sound velocity distribution of the propagation medium, and align it with the probe of the convex array ultrasonic transducer to generate an aligned sound velocity map;

[0007] The ultrasonic time domain signal p(r,φ,t) is converted at time t and circumferential φ to obtain a time-frequency domain spectrum and spatial frequency domain spectrum P n (r,ω), where ω represents the frequency and n is the spatial frequency corresponding to the Fourier transform with φ as the independent variable;

[0008] Based on the aligned sound velocity diagram, the sound field intensity within a specific frequency band Ω is selected and superimposed to obtain a passive cavitation image I(r,φ) in a coordinate system.

[0009] In one embodiment of the present invention, the ultrasonic time domain signal is converted in time and circumference, including:

[0010] Perform Fourier transform on the ultrasonic time domain signal p(r,φ,t) with respect to time t to convert the ultrasonic time domain signal p(r,φ,t) into a time-frequency domain spectrum

[0011] Time-frequency spectrum Perform Fourier transform on the circumferential direction φ to convert the time-frequency spectrum Converted into spatial frequency domain spectrum P n (r,ω).

[0012] In one embodiment of the present invention, based on the aligned sound velocity map, the sound field intensity within a specific frequency band Ω is selected and superimposed to obtain a passive cavitation image I(r,φ) in a coordinate system, including:

[0013] Determining the slowness u in the medium based on the aligned sound velocity map;

[0014] If r = r m The corresponding spatial frequency domain spectrum P n (r m ,ω), probe distribution surface r m+1 =r m +Δr, [r m ,r m+1 The medium slowness u in the interval is decomposed into the average slowness u0 of the interval and the deviation Δu between the slowness in the interval and the average slowness;

[0015] Based on the interval average slowness u0 and the deviation Δu between the slowness in the interval and the average slowness, the sound field intensity in a specific frequency band Ω is selected and superimposed to obtain a passive cavitation image I(r,φ) in the coordinate system.

[0016] In one embodiment of the present invention, based on the aligned sound velocity map, the sound field intensity within a specific frequency band Ω is selected and superimposed to obtain a passive cavitation image I(r,φ) in a coordinate system, including:

[0017] Calculation of the spiral spectrum by first-order direct integration P n (r m+1 ,ω) to perform inverse Fourier transform and obtain the frequency domain sound field

[0018] or

[0019] Based on the first-order direct integration method, r m The distributed surface wave field propagates to r m+1 Distribution surface, determine r m+1 Spiral spectrum of distribution surface P n (r m+1 ,ω), for the spiral spectrum P n (r m+1,ω) to perform inverse Fourier transform and obtain the frequency domain sound field

[0020] or

[0021] Determination of the frequency domain acoustic field by step-by-step spiral spectroscopy

[0022] In one embodiment of the present invention, the spiral spectrum is calculated by a first-order direct integration method. include:

[0023]

[0024] in,

[0025] k0=ωu0;

[0026]

[0027] is the second kind nth order Hankel equation;

[0028] is the nth-order Hankel equation of the first kind.

[0029] In one embodiment of the present invention, based on the first-order direct integration method, r m The distributed surface wave field propagates to r m+1 Distribution surface, determine r m+1 Spiral spectrum of distribution surface P n (r m+1 ,ω), including:

[0030]

[0031] in,

[0032]

[0033] In one embodiment of the present invention, the frequency domain sound field is determined by step-by-step spiral spectroscopy. include:

[0034]

[0035] The present invention provides a non-uniform spiral wave spectrum passive cavitation imaging device, comprising:

[0036] A convex array ultrasonic transducer, used for passively receiving an ultrasonic time domain signal p(r, φ, t) in a medium with non-uniform sound velocity, wherein the array elements of the convex array ultrasonic transducer are distributed on an arc of radius R, wherein r is the extrapolated radius of the current layer, φ is the circumferential angle in the coordinate system, t is time, and R is the distribution surface of the array elements of the convex array ultrasonic transducer;

[0037] An aligned sound velocity map generation module is used to determine the sound velocity distribution of the propagation medium based on the ultrasonic time domain signal p(r, φ, t), and align it with the probe of the convex array ultrasonic transducer to generate an aligned sound velocity map;

[0038] A processing module is used to convert the ultrasonic time domain signal p(r,φ,t) in time t and circumferential direction φ to obtain a time-frequency domain spectrum. and spatial frequency domain spectrum P n (r,ω), where ω represents the frequency;

[0039] The imaging module is used to obtain a passive cavitation image I(r,φ) in a coordinate system by selecting the sound field intensity within a specific frequency band Ω based on the aligned sound velocity map.

[0040] The electronic device provided by the present invention comprises:

[0041] one or more processors;

[0042] The storage device is used to store one or more programs. When the one or more programs are executed by the one or more processors, the electronic device implements the passive cavitation imaging method of sound velocity non-uniform spiral wave spectrum.

[0043] The computer-readable storage medium provided by the present invention stores a computer program thereon, and when the computer program is executed by a processor of a computer, the computer is enabled to execute the passive cavitation imaging method of sound velocity non-uniform spiral wave spectrum.

[0044] Beneficial effects of the present invention: Compared with the existing spiral spectrum method, the sound velocity heterogeneous spiral spectrum passive cavitation imaging method of the present invention can accurately locate the sound source in the sound velocity heterogeneous medium and more effectively suppress distortion errors and artifacts. Compared with the time delay accumulation integration method, it has the advantages of low algorithm complexity and low calculation time consumption.

[0045] It should be understood that the foregoing general description and the following detailed description are exemplary and explanatory only and are not restrictive of the present application. BRIEF DESCRIPTION OF THE DRAWINGS

[0046] The drawings herein are incorporated into the specification and constitute a part of the specification, showing embodiments consistent with the present application, and together with the specification, are used to explain the principles of the present application. Obviously, the drawings described below are only some embodiments of the present application, and for those of ordinary skill in the art, other drawings can be obtained based on these drawings without creative work. In the drawings:

[0047] Figure 1 It is a reference coordinate system shown in an exemplary embodiment of the present application.

[0048] Figure 2 It is a flow chart of a method for passive cavitation imaging of sound velocity non-uniform spiral wave spectrum shown in an exemplary embodiment of the present application.

[0049] Figure 3 It is a schematic diagram of a coordinate system for imaging of a convex array ultrasonic transducer in an inhomogeneous medium, as shown in an exemplary embodiment of the present application.

[0050] Figure 4 It is an experimental schematic diagram showing that a sound source is distributed in the shape of a blood vessel according to an exemplary embodiment of the present application.

[0051] Figure 5 It is an experimental schematic diagram showing that the sound sources are distributed as grid points according to an exemplary embodiment of the present application.

[0052] Figure 6 It is a schematic diagram showing an exemplary embodiment of the present application of using a translation stage to move a focusing probe to collect data at different positions.

[0053] Figure 7 It is a comparative display diagram of an HHWS image result of using sound velocity information to correct distortion and an HWS result using average sound velocity, shown as an exemplary embodiment of the present application.

[0054] Figure 8 1 is a PAM image of HWS and HHWS of a multi-cavitation source experiment shown in an exemplary embodiment of the present application.

[0055] Fig. 9 1 is a comparative diagram showing an exemplary embodiment of the present application showing a HHWS image result using an exact Hankel function and a HWS result using an average sound velocity without correcting distortion.

[0056] Fig.10 1 is a comparison diagram of HHWS image results using an exact Hankel function and an approximate function, shown in an exemplary embodiment of the present application. DETAILED DESCRIPTION

[0057] The following will describe the embodiments of the present invention with reference to the accompanying drawings and preferred embodiments. Those skilled in the art can easily understand other advantages and effects of the present invention from the contents disclosed in this specification. The present invention can also be implemented or applied through other different specific embodiments, and the details in this specification can also be modified or changed in various ways based on different viewpoints and applications without departing from the spirit of the present invention. It should be understood that the preferred embodiments are only for illustrating the present invention, not for limiting the scope of protection of the present invention.

[0058] It should be noted that the illustrations provided in the following embodiments are only schematic illustrations of the basic concept of the present invention, and thus the drawings only show components related to the present invention rather than being drawn according to the number, shape and size of components in actual implementation. In actual implementation, the type, quantity and proportion of each component may be changed arbitrarily, and the component layout may also be more complicated.

[0059] It is worth mentioning that passive acoustic mapping (PAM) reconstructs images by using ultrasonic transducers to passively receive acoustic cavitation signals, and visually monitors the cavitation activity generated during ultrasound treatment. Abdominal organs are an important application scenario for ultrasonic cavitation treatment, and abdominal ultrasound imaging mostly uses convex arrays to obtain a wider imaging field of view. However, the types of tissues in the abdominal cavity are complex, and there are heterogeneous components of sound velocity such as skin, fat, muscle and other organs in the abdominal wall. The distortion of sound waves propagating in the abdominal cavity will affect the quality of image reconstruction. At present, the method suitable for convex array sound velocity inhomogeneous media is mainly based on the passive cavitation imaging algorithm based on "delay accumulation and integration". In the sound velocity inhomogeneous medium, the delay accumulation and integration method calculates the delay by integrating the slowness along the sound propagation path, adjusts the waveform according to the delay, and then coherently accumulates the adjusted signal and calculates the signal energy as the image intensity information. This method has high algorithm complexity and time-consuming calculation. Frequency-wavenumber domain methods such as the angular spectrum method have lower computational complexity. In the rectangular coordinate system, the sound velocity heterogeneity angular spectrum method is applied to linear array transcranial passive cavitation imaging to overcome the distortion problem caused by sound velocity heterogeneity. However, the sound velocity heterogeneity angular spectrum method cannot be directly applied to convex arrays. In recent years, the helical wave spectrum passive cavitation imaging method (HWS) has been proposed as an extension of the angular spectrum method in the cylindrical coordinate system. It is suitable for passive cavitation imaging based on convex arrays in uniform media, but not for sound velocity inhomogeneous media.

[0060] In the following description, numerous details are discussed to provide a more thorough explanation of the embodiments of the present invention. However, it is obvious to those skilled in the art that the embodiments of the present invention can be implemented without these specific details. In other embodiments, well-known structures and devices are shown in the form of block diagrams rather than in detail to avoid making the embodiments of the present invention difficult to understand.

[0061] For example, the distribution surface of the array elements of the convex array ultrasonic transducer can be a cylindrical surface or a circumferential surface. For a convex array probe with one-dimensional array elements, y is always 0, the distribution surface is the circumferential surface, and the passive cavitation imaging method of the non-uniform spiral wave spectrum of sound velocity can be regarded as in the polar coordinate system. Next, establish (refer to the attached Figure 1 As shown). For a convex array probe (i.e., a probe of a convex array ultrasonic transducer) in which array elements are arranged two-dimensionally on a cylindrical surface, the y value is variable, and the distribution surface is a cylindrical surface, which can be regarded as being established in a cylindrical coordinate system.

[0062] In the embodiment of the present application, a one-dimensional arrangement of array elements is used as an example to illustrate the algorithm flow of the proposed convex array sound velocity heterogeneous passive cavitation imaging method as follows: Figure 2 As shown, when the sound field distribution of the peripheral surface r is known, the sound field of the peripheral surface r+Δr is calculated, where ω is the angular frequency corresponding to the Fourier transform of time t.

[0063] Step S100: data collection.

[0064] Specifically, refer to the attached Figure 3 As shown in the figure, in order to monitor the acoustic cavitation activity, a convex array ultrasonic transducer is used to passively receive the ultrasonic time domain signal p(r,φ,t) in a medium with non-uniform sound velocity. The array elements of the convex array ultrasonic transducer are distributed on a circular arc with a radius of R. Among them, r is the extrapolated radius of the current layer when acquiring images at different imaging depths, φ is the circumferential angle in the polar coordinate system, and t is time.

[0065] Step S200: generating an aligned sound velocity map.

[0066] Specifically, a B-ultrasound image is generated according to the ultrasonic time domain signal p(r,φ,t). The sound velocity distribution of the propagation medium is obtained according to the B-ultrasound image or sound velocity reconstruction, and the position of the ultrasonic probe (referring to the probe of the convex array ultrasonic transducer, referred to as the ultrasonic probe) is aligned to obtain an aligned sound velocity map, and the sound velocity of each sampling point is obtained by sampling according to the aligned sound velocity map, so as to determine the slowness of the sound wave in the medium.

[0067] Step S300, calculating the sound field distribution.

[0068] Step S310: Perform Fourier transform on the ultrasonic time domain signal p(r,φ,t) received by the ultrasonic probe with respect to time t. To convert the ultrasonic time domain signal p(r,φ,t) into the time-frequency domain spectrum Denoted as:

[0069]

[0070] Step S320: time-frequency spectrum Perform Fourier transform on the circumferential direction φ to convert the time-frequency spectrum Converted into spatial frequency domain spectrum P n (r,ω), the process is expressed as:

[0071]

[0072] Among them, n is the spatial frequency corresponding to the Fourier transform with φ as the independent variable, that is, the circumferential wave number.

[0073] For example, assuming that the circumference r=r m The corresponding spatial frequency domain spectrum is P n (r m ,ω), the convex array can only receive the sound waves emitted inward by the cavitation source located outside, so we only consider the inward waves. By aligning the sound velocity diagram, we can know that the slowness in the medium u(l) = 1 / c(l), where c(l) is the sound velocity in the medium with spatial coordinates l = (r,φ). Let r m+1 =r m +Δr, for [r m ,r m+1 ] interval, the medium slowness in the interval is decomposed into two parts, u0 and Δu, denoted as u(r,φ)=u0+Δu(r,φ), where u0 is the average slowness of the interval (constant), which is independent of the position coordinates φ and r, and Δu is the deviation of the slowness in the interval from the average slowness.

[0074] Exemplarily, the sound velocity heterogeneous spiral spectrum method suitable for convex arrays (i.e., the sound velocity non-uniform spiral spectrum passive cavitation imaging method) can be divided into three methods: first-order, second-order direct integration method and step-by-step spiral spectrum method, according to the different sound field propagation calculation processes. Among them, the step-by-step spiral spectrum method has the simplest algorithm and the highest robustness, but it has a slight error compared to the direct integration method. Therefore, it is recommended to use the step-by-step spiral spectrum method for preliminary verification of the imaging effect, and then replace it with the direct integration method when the imaging effect needs to be improved. The second-order direct integration method is used when the axial interval Δr<λ / 2, and the first-order direct integration method is used when Δr<λ / 2, where λ is the wavelength of the sound wave.

[0075] Method 1: The first-order direct integration method uses the first-order forward Euler method for calculation, as follows:

[0076]

[0077] in,

[0078] k0=ωu0, is the nth-order Hankel equation of the second kind, is the nth-order Hankel equation of the first kind.

[0079] Method 2: If the calculation result of the first-order direct integration method will be used in subsequent calculations, in order to facilitate distinction, the calculation result of the first-order direct integration method is recorded as It can be used in the second-order direct integration method for calculation.

[0080] Second-order direct integration method Based on the first-order direct integration method, the second-order prediction-correction method is used to convert r m The surface wave field propagates to r m+1 The trapezoidal numerical integration method is used to calculate the surrounding surface and obtain r m+1 The spiral spectrum P of the circumference n (r m+1 ,ω). The specific process is as follows:

[0081]

[0082] in,

[0083] is the inverse Fourier transform.

[0084] For the spiral spectrum P calculated by the first-order direct integration method or the second-order direct integration method n (r m+1 ,ω), perform inverse Fourier transform to obtain the frequency domain sound field It is expressed as:

[0085]

[0086] Method 3: Determine the frequency domain sound field based on the step-by-step spiral spectrum method Its calculation is expressed as:

[0087]

[0088] For the above three methods, the Hankel equation in the first-order direct integration method, the second-order direct integration method or the step-by-step spiral spectrum method can be simplified by using asymptotic analytical solutions (i.e., using asymptotic formulas), including:

[0089]

[0090] in,

[0091] for The conjugation of.

[0092] The calculation involving the Hankel equation simplifies to:

[0093]

[0094] in,

[0095]

[0096] For the non-evanescent wave part, the asymptotic formula can be accurately approximated, reducing the computational cost.

[0097] Starting from the ultrasonic probe surface r = R, the above method is used to iteratively calculate the spectrum distribution of the entire imaging space, and the sound field intensity in a specific frequency band Ω is selected and superimposed to obtain the passive cavitation image I(r,φ) in the polar coordinate system, as shown below:

[0098]

[0099] Step S400: coordinate system conversion

[0100] Specifically, the polar coordinate sampling points are converted into images in a rectangular coordinate system through an interpolation algorithm for easy observation and evaluation. The intersection of the arc where the probe of the convex array ultrasonic transducer is located and the symmetry axis is taken as the origin, the axial direction is the x-axis, and the transverse direction is the z-direction.

[0101] This embodiment starts from the sound velocity weak heterogeneity wave equation in the polar coordinate system, derives the sound velocity heterogeneous spiral wave spectrum propagation method in the frequency-wavenumber domain, and applies it to the passive cavitation imaging of convex arrays. The sound velocity heterogeneous spiral wave spectrum method (Heterogeneous Helical Wave Spectrum Method, HHWS) can be divided into the first-order direct integration method, the second-order direct integration method and the step-by-step spiral wave spectrum method according to the different sound field propagation calculation processes. Compared with the spiral wave spectrum method, the sound velocity heterogeneous spiral wave spectrum passive cavitation imaging method can accurately locate the sound source in the sound velocity non-uniform medium and more effectively suppress distortion errors and artifacts. Compared with the time-delay cumulative integration method, it has the advantages of low algorithm complexity and low calculation time. The method proposed in this embodiment can correct the phase distortion caused by sound velocity heterogeneity, and provide key technical support for real-time and accurate monitoring and guidance of abdominal ultrasonic cavitation treatment. The sound velocity heterogeneous spiral wave spectrum propagation method derived in this embodiment is a general solution to the sound velocity weak heterogeneity wave equation. As a fast solution method for the circumferential wave field, it can also be used to reconstruct B-mode ultrasound images or as a forward model for solving the circumferential sound velocity field.

[0102] Experiment 1: Simulation experiment

[0103] Step S100: Data collection

[0104] k-Wave was used to simulate a 128-channel convex array ultrasonic transducer array (aperture radius R = 49.6 mm, scanning angle 74.6°) to receive point source acoustic signals in a medium with non-uniform sound velocity. The signal sampling frequency was 13.89 MHz. The acoustic characteristics of the medium were set to simulate human soft tissue (1540 m / s, 1050 kg / m 3 , 0.54dB / cm / MHz), fat (1430m / s, 928kg / m 3 , 0.60dB / cm / MHz) and human abdominal wall muscles (1580m / s, 1041kg / m 3 , 0.57dB / cm / MHz). In soft tissue, 200 point sources were randomly placed in a vascular shape area about 50mm deep from the convex array (Figure (4(a)). In addition, 9×13 point sources were placed in a rectangular grid with lateral and axial position ranges of -30 to 30mm and 40 to 80mm, respectively (Figure (5(a)). These point sources sequentially emitted 100-cycle sine signals of 2MHz.

[0105] Step S200: Acquire an aligned sound velocity map

[0106] The sound velocity diagram shows the sound velocity distribution of the simulation settings.

[0107] Step S300: Calculation of sound field distribution

[0108] Taking the second-order direct integration method as an example, according to step S310, the received microbubble cavitation signal is Fourier transformed in the time t direction. Then according to step S320, Fourier transform is performed in the circumferential direction to obtain the spatial frequency domain spectrum P n (r,ω). From r = R to R + 100 mm depth, the spatial frequency domain spectrum P is derived by iterating layer by layer with an interval of Δr = 0.2 mm using the first-order direct integration method and the second-order direct integration method. n (r,ω). Inverse Fourier transform to time-frequency domain spectrum After that, the sound field intensity of the single frequency component corresponding to the selected sound source is superimposed to obtain the passive cavitation image in polar coordinates, and the sampling grid is 384×501 (φ×r).

[0109] Step S400: Coordinate system conversion

[0110] After conversion to the rectangular coordinate system, the simulated sound source setting position is used as the positioning reference (l r ). Compare the HHWS image results using sound velocity information to correct the distortion with the HWS results using the average sound velocity. For the HWS method, the positioning errors of the 2MHz, 3MHz, and 4MHz signal sources (positioned at the peak position of the sound intensity, positioning error: |ll r|) were 4.8±2.1mm, 4.2±2.0mm and 3.4±1.3mm respectively. After using the HHWS method, the errors were reduced to 0.5±0.3mm, 0.7±0.3mm and 0.8±0.3mm respectively. After correction, the positioning signal source and the shape of the blood vessel can be well matched. For details, please refer to the attached Figure 4 (b)-(c).

[0111] Attached Figure 4 (a) is a schematic diagram of a simulation experiment in which sound sources are distributed in the shape of blood vessels. (b) compares the localized (orange dots) and true (grey dots) sound source locations. The color of the localized points is normalized according to the maximum intensity at that frequency (2MHz, 3MHz, and 4MHz). (c) PAM image of the site indicated by the arrow (the red cross marks the true sound source location).

[0112] For the experiment where the sound source is distributed in a grid, the average localization error of the HHWS method is reduced from 5.6±3.0mm to 0.7±0.5mm compared with the HWS method. The average energy spread area (ESA, the area in the image where the intensity is greater than half of the maximum value) is reduced from 20.9±18.0mm to 1.3±0.5mm. 2 Down to 8.8±4.0mm 2 , refer to the attached Figure 5 (b)-(c), attached Figure 5 (a) Schematic diagram of the simulation experiment in which the sound source is distributed as a grid point. (b) Positioning error of PAM images reconstructed using the HWS method and the HHWS method, and (c) ESA data.

[0113] Experiment 2: Single cavitation source experiment

[0114] Step S100: Data collection

[0115] A convex array ultrasonic transducer with 128 channels, a curvature radius of 49.57 mm, and an imaging field of view of 74.6° was used to receive the acoustic cavitation signals generated by the stimulation of ultrasonic microbubbles for reconstructing passive cavitation images. The center frequency of the transducer was 3.57 MHz, and the sampling frequency was 13.89 MHz. An elliptical cylindrical phantom with a speed close to that of muscle sound (1580 m / s) was placed under the convex array ultrasonic transducer in water (1493 m / s). A 3D-printed bracket was used to fix the relative position of the focused ultrasonic probe with a center frequency of 1 MHz and the rectangular wallless phantom, which was connected to the translation stage and placed under the elliptical phantom. The focused ultrasonic probe excited microbubbles on the channel at a repetition frequency of 30 Hz to produce cavitation activity. The focused sound pressure amplitude was set to 800 kPa, and the focused ultrasonic probe was moved using a translation stage to collect data at 9 different positions centered on [0 mm, 60 mm] ( Figure 6 ).

[0116] Step S200: Acquire an aligned sound velocity map

[0117] The elliptical phantom is made by 3D printing mold, the semi-axis length is fixed, and the sound velocity map is reconstructed by identifying the position of the elliptical phantom through B-mode image, such as Figure 4 shown.

[0118] Step S300: Calculation of sound field distribution

[0119] Taking the second-order direct integration method as an example, according to step S310, the received microbubble cavitation signal is Fourier transformed in the time direction. According to step S320, Fourier transform is performed in the circumferential direction to obtain the spatial frequency domain spectrum P n (r,ω). At depths from r = R to R + 100 mm, the spectrum P is derived by iterating layer by layer with intervals of Δr = 0.2 mm using the first-order direct integration method and the second-order direct integration method. n (r,ω). Inverse Fourier transform to time-frequency domain spectrum Finally, the acoustic field intensity of each frequency in the selected frequency band Ω is superimposed to obtain the passive cavitation image in polar coordinates, and the sampling grid is 384×501 (φ×r). All frequency components (1 to 6 MHz) covering the working bandwidth of the ultrasound probe are selected to reconstruct the panoramic image.

[0120] Step 4: Coordinate system conversion

[0121] After conversion to the rectangular coordinate system, the reconstruction result of the HWS method in a uniform medium (water) with the elliptical high-speed sound body removed is used as the positioning reference (l r ). Select the data with the sound source position at (0.8mm, 48.8mm) and (-8.8mm, 58.4mm), and compare the HHWS image results using sound velocity information to correct the distortion with the HWS results using the average sound velocity, as shown below: Figure 7 Attached Figure 7 Example of a single cavitation source experiment (a) PAM image of the sound source located at (0.8mm, 48.8mm). (b) PAM image of the sound source located at (-8.8mm, 58.4mm). The cross mark is the reference position of the sound source.

[0122] The localization errors and ESA data of the 9 sound source positions are listed in Table (1). In the absence of correction (HWS), the average localization error (localization at the sound intensity peak position l, localization error: |ll r |), the energy spread area (ESA, the area in the image where the intensity is greater than half of the maximum value) ranges from 6.4-12.8mm and 32.3-95.5mm, respectively 2 After correction (HHWS), these values ​​range from 0.8-1.9mm and 10.3-22.8mm respectively. 2After correction, both the image quality and positioning accuracy are improved.

[0123] Table 1 Positioning error of single cavitation source experiment and ESA data

[0124]

[0125]

[0126] Experiment 3: Multiple cavitation source experiment

[0127] Step S100: Data collection

[0128] A convex array ultrasonic transducer (128 channels, a radius of curvature of 49.57 mm, and an imaging field of view of 74.6°) was used to passively receive microbubble cavitation signals (i.e., ultrasonic time domain signals) and reconstruct passive cavitation images. The center frequency of the transducer (short for convex array ultrasonic transducer) was 3.57 MHz, and the sampling frequency was 13.89 MHz. Two elliptical phantoms with a sound velocity of 1580 m / s were placed under the transducer in water (23°, 1493 m / s). A focused ultrasonic probe with a center frequency of 1 MHz was fixed to a dual-channel rectangular wallless phantom under the elliptical phantom, and microbubbles were excited to produce cavitation activity on the channel at a repetition frequency of 30 Hz and a sound pressure of 400 kPa.

[0129] Step S200: Acquire an aligned sound velocity map

[0130] Specifically, the sound velocity map is reconstructed by identifying the position of the elliptical phantom through B-mode images.

[0131] Step S300: Calculation of sound field distribution

[0132] Taking the second-order direct integration method as an example, according to step S310, the received microbubble cavitation signal is Fourier transformed in the time direction. According to step S320, Fourier transform is performed in the circumferential direction to obtain the spatial frequency domain spectrum P n (r,ω). At a depth of 0 to 100 mm, the formula in the first-order direct integration method and the formula in the second-order direct integration method are iterated layer by layer at intervals of Δr = 0.2 mm to derive the spectrum P n (r,ω). Inverse Fourier transform to time-frequency domain spectrum Finally, the acoustic field intensity of each frequency in the selected frequency band Ω is superimposed to obtain the passive cavitation image in polar coordinates, and the sampling grid is 384×501 (φ×r). All frequency components (1 to 6 MHz) covering the ultrasound probe frequency band are selected to reconstruct the panoramic image.

[0133] Step 4: Coordinate system conversion

[0134] After conversion to the rectangular coordinate system, the reconstruction result of the HWS method in a uniform medium (water) with the elliptical high-speed sound body removed is used as the positioning reference (l r ), the HHWS image results using sound velocity information to correct distortion are compared with the HWS results using the average sound velocity, as shown in Figure 8 Compared with the corrected case, the positioning error (positioning at the peak position of the sound intensity l, positioning error: |ll r |) decreased from 4.3±2.3mm (8.6±9.1mm) to 1.6±3.4mm (1.9±2.6mm) at the left (right) site, and ESA decreased from 47.0±30.1mm at the left (right) site 2 Reduced to 14.7±8.5mm 2 (from 79.8±48.8mm 2 Reduced to 20.0±16.0mm 2 ).

[0135] Experiment 4: Dual-modal experiment

[0136] Step S100: Data collection

[0137] A convex array ultrasonic transducer (128 channels, radius of curvature of 49.57 mm, imaging field of view of 74.6°) was used to passively receive microbubble cavitation signals and reconstruct passive cavitation images. The center frequency of the transducer was 3.57 MHz and the sampling frequency was 13.89 MHz. A fan-shaped wallless phantom was embedded in a high-speed elliptical phantom and placed under the transducer. The focused ultrasonic probe excited microbubbles on the channel to produce cavitation activity at a repetition frequency of 30 Hz and a focused sound pressure amplitude of 800 kPa.

[0138] Step S200: Acquire an aligned sound velocity map

[0139] The sound velocity map is reconstructed by identifying the position of the elliptical phantom through B-mode images.

[0140] Step S300: Calculation of sound field distribution

[0141] Taking the second-order direct integration method as an example, according to step S310, the received microbubble cavitation signal is Fourier transformed in the time direction. According to step S320, Fourier transform is performed in the circumferential direction to obtain the spatial frequency domain spectrum P n (r, ω). In the imaging range from r = R to r = R + 100 mm, iterate layer by layer with an interval of Δr = 0.2 mm. The spatial frequency domain spectrum P is derived by using the formula in the second-order direct integration method and the formula in the second-order direct integration method after replacing the asymptotic formula. n (r,ω). Inverse Fourier transform to time-frequency domain spectrum Afterwards, the acoustic field intensity of each frequency (1-6 MHz) in the selected frequency band Ω is superimposed to obtain the passive cavitation image in polar coordinates, and the sampling grid is 384×501 (φ×r).

[0142] Step 4: Coordinate system conversion

[0143] After conversion to the rectangular coordinate system, the HHWS image results using the exact Hankel function (Matlab "besselh" function) are compared with the HWS results using the average sound velocity without distortion correction, as shown in the figure. Fig. 9 As shown in Figure 2, the passive cavitation image is correctly focused on the channel when compared with the uncorrected and corrected cases. The HHWS image results using the exact Hankel function (the formula in the second-order direct integration method) and the approximate function (the formula in the second-order direct integration method after replacing the asymptotic formula) are compared, as shown in Figure 2. Fig.10 The difference between the PAM image result of the approximate function and the exact Hankel function result is 10 -3 The HHWS method can achieve real-time imaging. For a 1(739) frequency window, the calculation time on the CPU using the accurate Hankel function is 102.1±3.1ms (74.4±1.2s) and on the GPU is 18.9±0.5ms (66.9±0.7ms). After using the approximate Hankel function, the calculation time on the CPU and GPU in the 739 frequency window is reduced to 27.8±0.4s and 46.7±0.8ms, respectively.

[0144] Reference Fig. 9 As shown, attached Fig. 9 The PAM image and B-mode image are superimposed to show the sound source at (a) and the sound source position moved 10 mm to the right (b) using the accurate Hankel function reconstruction. The orange arrow in the figure indicates the flow direction of the microbubble (MB) suspension.

[0145] Reference Fig.10 As shown, attached Fig.10 PAM images reconstructed using (a) the exact Hankel function and (b) the approximate function. (c) The difference between the PAM images reconstructed using the exact Hankel function and the approximate function. All images were evaluated by normalizing them using the maximum value in the image formed using the exact Hankel function.

[0146] In view of the above-mentioned sound velocity non-uniform spiral wave spectrum passive cavitation imaging method, the present application also provides a sound velocity non-uniform spiral wave spectrum passive cavitation imaging device, comprising:

[0147] A convex array ultrasonic transducer, used for passively receiving an ultrasonic time domain signal p(r, φ, t) in a medium with non-uniform sound velocity, wherein the array elements of the convex array ultrasonic transducer are distributed on an arc of radius R, wherein r is the extrapolated radius of the current layer, φ is the circumferential angle in the polar coordinate system, and t is time;

[0148] An aligned sound velocity map generation module is used to determine the sound velocity distribution of the propagation medium based on the ultrasonic time domain signal p(r, φ, t), and align it with the ultrasonic probe of the convex array ultrasonic transducer to generate an aligned sound velocity map;

[0149] A processing module is used to convert the ultrasonic time domain signal p(r,φ,t) in time t and circumferential direction φ to obtain a time-frequency domain spectrum. and spatial frequency domain spectrum P n (r,ω), where ω represents the frequency;

[0150] An imaging module, used to derive a spatial frequency domain spectrum P of an imaging space based on the aligned sound velocity map n (r, ω), and the sound field intensity in a specific frequency band Ω is selected and superimposed to obtain the passive cavitation image I(r, φ) in the polar coordinate system.

[0151] The present application also provides a device, comprising:

[0152] one or more processors and memory,

[0153] The memory stores a computer program, and when the one or more processors execute the computer program, the device executes the passive cavitation imaging method of sound velocity non-uniform spiral wave spectrum.

[0154] Another aspect of the present application also provides a computer-readable storage medium on which a computer program is stored. When the computer program is executed by a processor of a computer, the computer executes the passive cavitation imaging method of non-uniform spiral wave spectrum of sound velocity as described above. The computer-readable storage medium may be included in the electronic device described in the above embodiment, or may exist independently without being assembled into the electronic device.

[0155] Another aspect of the present application also provides a computer program product or a computer program, which includes a computer instruction stored in a computer-readable storage medium. A processor of a computer device reads the computer instruction from the computer-readable storage medium, and the processor executes the computer instruction, so that the computer device executes the passive cavitation imaging method of the sound velocity non-uniform spiral wave spectrum provided in the above-mentioned embodiments.

[0156] The above embodiments are merely illustrative of the principles and effects of the present invention, and are not intended to limit the present invention. Anyone familiar with the technology may modify or change the above embodiments without violating the spirit and scope of the present invention. Therefore, all equivalent modifications or changes made by a person of ordinary skill in the art without departing from the spirit and technical ideas disclosed by the present invention shall still be covered by the claims of the present invention.

Claims

1. A method for passive cavitation imaging of sound velocity non-uniform spiral wave spectrum, characterized in that: include: A convex array ultrasonic transducer passively receives an ultrasonic time domain signal p(r,φ,t) in a medium with non-uniform sound velocity, wherein the array elements of the convex array ultrasonic transducer are distributed on an arc of radius R, wherein r is the extrapolated radius of the current layer, φ is the circumferential angle in the coordinate system, t is time, and R is the distribution surface of the array elements of the convex array ultrasonic transducer; Based on the ultrasonic time domain signal p(r,φ,t), determine the sound velocity distribution of the propagation medium, and align it with the probe of the convex array ultrasonic transducer to generate an aligned sound velocity map; The ultrasonic time domain signal p(r,φ,t) is converted at time t and circumferential φ to obtain a time-frequency domain spectrum and spatial frequency domain spectrum P n r,ω), where ω represents the frequency and n is the spatial frequency corresponding to the Fourier transform with φ as the independent variable; Based on the aligned sound velocity diagram, the sound field intensity within a specific frequency band Ω is selected and superimposed to obtain a passive cavitation image I(r,φ) in a coordinate system.

2. The passive cavitation imaging method of sound velocity non-uniform spiral wave spectrum according to claim 1, characterized in that: The ultrasonic time domain signal is converted in time and circumference, comprising: Perform Fourier transform on the ultrasonic time domain signal p(r,φ,t) with respect to time t to convert the ultrasonic time domain signal p(r,φ,t) into a time-frequency domain spectrum Time-frequency spectrum Perform Fourier transform on the circumferential direction φ to convert the time-frequency spectrum Converted into spatial frequency domain spectrum P n (r,ω).

3. The passive cavitation imaging method of sound velocity non-uniform spiral wave spectrum according to claim 1, characterized in that: Based on the aligned sound velocity map, the sound field intensity in a specific frequency band Ω is selected and superimposed to obtain a passive cavitation image I(r,φ) in the coordinate system, including: Determining the slowness u in the medium based on the aligned sound velocity map; If r = r m The corresponding spatial frequency domain spectrum is P n (r m ,ω), probe distribution surface r m+1 =r m +Δr, [r m ,r m+1 The medium slowness u in the interval is decomposed into the average slowness u0 of the interval and the deviation Δu between the slowness in the interval and the average slowness; Based on the interval average slowness u0 and the deviation Δu between the slowness in the interval and the average slowness, the sound field intensity in a specific frequency band Ω is selected and superimposed to obtain a passive cavitation image I(r,φ) in the coordinate system.

4. The method for passive cavitation imaging of sound velocity non-uniform spiral wave spectrum according to claim 3, characterized in that: Based on the aligned sound velocity map, the sound field intensity in a specific frequency band Ω is selected and superimposed to obtain a passive cavitation image I(r,φ) in the coordinate system, including: Calculation of spiral wave spectra by first-order direct integration method P n (r m+1 ,ω) to perform inverse Fourier transform and obtain the frequency domain sound field or Based on the first-order direct integration method, r m The distributed surface wave field propagates to r m+1 Distribution surface, determine r m+1 Spiral spectrum of distribution surface P n (r m+1 ,ω); for the spiral spectrum P n (r m+1 ,ω) to perform inverse Fourier transform and obtain the frequency domain sound field or Determination of the frequency domain acoustic field by step-by-step spiral spectroscopy 5. The method for passive cavitation imaging of sound velocity non-uniform spiral wave spectrum according to claim 4, characterized in that: Calculation of spiral wave spectra by first-order direct integration method include: in, k0=ωu0; is the second kind nth order Hankel equation; is the nth-order Hankel equation of the first kind.

6. The method for passive cavitation imaging of sound velocity non-uniform spiral wave spectrum according to claim 4, characterized in that: Based on the first-order direct integration method, r m The distributed surface wave field propagates to r m+1 Distribution surface, determine r m+1 Spiral spectrum of distribution surface P n (r m+1 ,ω), include: in, 7. The passive cavitation imaging method of sound velocity non-uniform spiral wave spectrum according to claim 1, characterized in that: Determination of the frequency domain acoustic field by step-by-step spiral spectroscopy include:

8. A passive cavitation imaging device for non-uniform spiral wave spectrum of sound velocity, characterized in that: include: A convex array ultrasonic transducer, used for passively receiving an ultrasonic time domain signal p(r, φ, t) in a medium with non-uniform sound velocity, wherein the array elements of the convex array ultrasonic transducer are distributed on an arc of radius R, wherein r is the extrapolated radius of the current layer, φ is the circumferential angle in the coordinate system, t is time, and R is the distribution surface of the array elements of the convex array ultrasonic transducer; An aligned sound velocity map generation module is used to determine the sound velocity distribution of the propagation medium based on the ultrasonic time domain signal p(r, φ, t), and align it with the probe of the convex array ultrasonic transducer to generate an aligned sound velocity map; A processing module is used to convert the ultrasonic time domain signal p(r,φ,t) in time t and circumferential direction φ to obtain a time-frequency domain spectrum. and spatial frequency domain spectrum P n r,ω), where ω represents the frequency; The imaging module is used to obtain a passive cavitation image I(r,φ) in a coordinate system by selecting the sound field intensity within a specific frequency band Ω based on the aligned sound velocity map.

9. A device, characterized in that: include: one or more processors and memory, A computer program is stored in the memory, and when the one or more processors execute the computer program, the device is caused to perform the method according to any one of claims 1 to 7.

10. A computer-readable storage medium, characterized in that: A computer program is stored thereon, which, when executed by one or more processors, causes the device to perform the method according to any one of claims 1 to 7.

Images