CN115969411B

CN115969411B - Super-resolution ultrasonic passive cavitation imaging method and system based on single cavitation source separation and positioning

Info

Publication number: CN115969411B
Application number: CN202310060777.6A
Authority: CN
Inventors: 路舒宽; 苏瑞波; 万明习; 万春野
Original assignee: Xian Jiaotong University
Current assignee: Xian Jiaotong University
Priority date: 2023-01-16
Filing date: 2023-01-16
Publication date: 2023-07-11
Anticipated expiration: 2043-01-16
Also published as: CN115969411A

Abstract

The invention discloses a super-resolution ultrasonic passive cavitation imaging method and system based on single cavitation source separation and positioning: setting a time gating interval, processing signals in the interval to obtain a time gating energy diagram, carrying out thresholding on the time gating energy diagram, clustering pixels to obtain a global pixel cluster, carrying out thresholding on the local energy diagram, clustering the pixels to obtain a local pixel cluster, searching a main cluster from the local pixel cluster to obtain a single cavitation source energy diagram so as to separate a single cavitation source, removing the single cavitation source energy diagram with a pixel peak value smaller than a threshold value, carrying out Gaussian fitting to position the single cavitation source, detecting whether the positioning of the cavitation source is normal according to the threshold value of the uncertainty and the fitting goodness of the cavitation source coordinate estimation, calculating a positioning distribution diagram of the cavitation source which is positioned normally, and superposing to obtain a super-resolution imaging result. The invention greatly improves the resolution of ultrasonic passive cavitation imaging and realizes the fine analysis of the spatial distribution of cavitation sources in a small area.

Description

Super-resolution ultrasonic passive cavitation imaging method and system based on single cavitation source separation and positioning

Technical Field

The invention belongs to the technical field of ultrasonic detection and ultrasonic imaging, and relates to a super-resolution ultrasonic passive cavitation imaging method and system for obtaining fine resolution of cavitation source spatial distribution in a small area by separating and positioning a single cavitation source.

Background

Acoustic cavitation refers to the physical phenomenon of cavitation bubble formation, expansion and shrinkage and collapse in a medium under the excitation of ultrasonic waves emitted by a cavitation excitation ultrasonic transducer, enhances the interaction between the ultrasonic waves and organisms, and is a core physical foundation for ultrasonic treatment applications such as tumor ablation, thrombolysis, drug release and the like. Cavitation detection helps to monitor the ultrasound therapy process and further clarify the physical mechanism behind it, with the most common means being optical and acoustic detection methods. Optical detection methods such as high-speed photography and sonoluminescence can intuitively observe cavitation change process and have high spatial resolution, but are limited by transparency of media and cannot be used for detecting cavitation in biological tissues. The acoustic detection method can quantitatively analyze the acoustic characteristics of cavitation in biological tissues, and can be divided into active detection and passive detection according to the working mode of a detection transducer, wherein the latter has the advantage of synchronously detecting cavitation in the ultrasonic excitation process. The passive acoustic detection method is further divided into single array element detection and array detection (namely ultrasonic passive cavitation imaging) according to the types of detection transducers, wherein the latter overcomes the defect that the former cannot provide cavitation spatial position and distribution information, and the former has been paid attention to in recent years.

Ultrasonic passive cavitation imaging is generally to perform beam synthesis processing based on relative transit time information on cavitation signals passively received by an array ultrasonic transducer and calculate acoustic energy according to the obtained signals, wherein the spatial resolution performance of the ultrasonic transducer depends on diffraction modes (including factors such as aperture and receiving bandwidth) of the array ultrasonic transducer. The current array ultrasonic transducer for implementing ultrasonic passive cavitation imaging mainly comprises transducers (such as a linear array, a phased array and the like) commonly used for clinical diagnosis and a hemispherical array transducer special for transcranial imaging, wherein the former is applicable to different parts of a human body, and has wider application prospect. However, the limited aperture of the diagnostic transducer limits the spatial resolution (especially the axial resolution) of the cavitation image, while the effects of tissue heterogeneity, source interactions, etc. often accompany high levels of interference artifacts in the cavitation image. Currently researchers have improved image quality by mainly improving and optimizing imaging methods. A typical imaging method adopts a self-adaptive beam forming technology, including beam forming based on signal second-order or higher-order statistical characteristics, beam forming based on covariance matrix feature decomposition and the like, optimizes weight coefficients applied to array elements and performs beam forming accordingly, so that interference artifacts are obviously eliminated, and imaging resolution is enhanced. In addition, researchers also apply amplitude/phase coherent weighting factors to beam forming signals, perform coupling multiplication on amplitudes at two frequencies with equal increments on both sides of the center frequency of the signals, perform symbol squaring, coupling multiplication and superposition on delay signals, perform complementary apodization processing on array signals, and improve imaging methods in such ways that the beam forming signals are weighted by cross-correlation factors, and the researches show that the methods are beneficial to improving imaging resolution.

Although the various imaging methods currently developed exhibit good performance, they do not fundamentally overcome the problem of diffraction limitations of imaging transducers, so that a breakthrough improvement in imaging resolution is difficult to achieve. In acoustic cavitation and related therapeutic ultrasound applications, detection of cavitation in a localized small region of space is critical to the study of acoustic cavitation physics itself and the explanation and clarification related biophysical mechanisms. In view of this, it is needed to provide a super-resolution ultrasonic passive cavitation imaging method and system capable of performing fine resolution on the spatial distribution of cavitation sources in a small area, and this is also a technical problem that has been needed to be overcome in the field of ultrasonic passive cavitation imaging.

Disclosure of Invention

The invention aims to provide a super-resolution ultrasonic passive cavitation imaging method and system based on single cavitation source separation and positioning.

In order to achieve the above purpose, the present invention adopts the following technical scheme:

a super-resolution ultrasonic passive cavitation imaging method based on single cavitation source separation and positioning comprises the following steps:

1) Establishing an imaging coordinate system and planning a pixel grid of a time gating energy chart, calculating double-pass delay according to the distance from each pixel coordinate to a tangent line or a tangent plane at the center coordinate of the surface of the cavitation excitation ultrasonic transducer in the pixel grid and the distance between each pixel coordinate and each array element coordinate of the array ultrasonic transducer, taking the double-pass delay as the initial time of the time gating interval and setting the time length of the time gating interval to obtain the time gating interval of each array element, extracting the time gating signal of each array element from each frame of original radio frequency signal passively received by the array ultrasonic transducer according to the time gating interval of each array element, superposing the time gating signals of each array element along the array element direction to obtain a time gating composite signal, superposing the square of the time gating composite signal along the sampling point direction to obtain the time gating energy chart corresponding to each frame of original radio frequency signal;

2) After thresholding is carried out on each frame of time gating energy diagram, a global pixel set is established according to all non-zero value pixels, and according to the set closed neighborhood radius and the set closed neighborhood pixel number threshold, the pixels in the global pixel set are clustered by repeating the processes of initializing pixel clusters and expanding pixel clusters, so that the global pixel cluster of each frame of time gating energy diagram is obtained;

3) Defining local clustering windows for each global pixel cluster of each frame time gating energy diagram, carrying out frame selection on the time gating energy diagram by utilizing the local clustering windows to obtain a local energy diagram, establishing a local pixel set according to all non-zero value pixels after thresholding the local energy diagram, clustering pixels in the local pixel set by repeating the processes of initializing the pixel clusters and expanding the pixel clusters according to the set closed neighborhood radius and the set closed neighborhood pixel number threshold value to obtain a local pixel cluster, searching a local pixel main cluster from the local pixel cluster, and setting the pixel values of all pixels which are not contained in the local pixel main cluster in the local energy diagram to be zero to obtain a single cavitation source energy diagram, thereby realizing the separation of single cavitation sources;

4) Extracting pixel peaks of all the single cavitation source energy maps, judging the pixel peak value of each single cavitation source energy map according to the set pixel peak value threshold, and establishing a single cavitation source energy map set by the single cavitation source energy maps with the pixel peak value smaller than the pixel peak value threshold;

5) Establishing a Gaussian distribution function in an imaging coordinate system, performing Gaussian fitting on each single cavitation source energy diagram in a single cavitation source energy diagram set to position a single cavitation source, obtaining an estimated value of a parameter in the Gaussian distribution function, wherein the estimated value of a Gaussian distribution peak value position is an estimated value of a cavitation source coordinate, substituting the estimated value of the parameter into the Gaussian distribution function to calculate a Gaussian fitting result, calculating the uncertainty of the estimation of the cavitation source coordinate under a given confidence level, and calculating the fitting goodness according to the single cavitation source energy diagram and the Gaussian fitting result of the single cavitation source energy diagram;

6) Establishing a cavitation source coordinate estimation uncertainty set and a fitting goodness set according to the cavitation source coordinate estimation uncertainty and the fitting goodness of all the cavitation sources, respectively carrying out ascending arrangement on elements in each set, calculating quartiles and quartiles distances of the elements in each set after ascending arrangement, setting a threshold value of the cavitation source coordinate estimation uncertainty and a threshold value of the fitting goodness according to the quartiles and the quartiles distances, and detecting the normally positioned cavitation source by utilizing the set threshold value of the cavitation source coordinate estimation uncertainty and the threshold value of the fitting goodness;

7) Planning a super-resolution imaging pixel grid, establishing a positioning distribution function aiming at each normally positioned cavitation source, calculating a positioning distribution function value of each pixel coordinate in the pixel grid to obtain a positioning distribution diagram of each normally positioned cavitation source, superposing the positioning distribution diagrams of all normally positioned cavitation sources, and carrying out normalization and logarithmic processing to obtain a super-resolution imaging result.

Preferably, in the step 1, the calculation formula of the two-way delay is expressed as:

wherein,,

distance from any one pixel coordinate (x, z) in a pixel grid of the time-gated energy map to a tangent or tangent plane at a center coordinate of the cavitation-excited ultrasonic transducer surface, +.>

For the pixel coordinates (x, z) and the coordinates (x) of the ith array element of the array ultrasonic transducer _i Distance between 0), ∈0)>

i＝1,2,...,N _SE ，N _SE The number of array elements of the array ultrasonic transducer is the sound propagation speed;

when the central axis of the nulling excitation ultrasound transducer is in the imaging plane of the array ultrasound transducer,

exciting the ultrasonic transducer surface center coordinates (x) for the pixel coordinates (x, z) to cavitation _E ,y _E ,z _E ) Distance of tangent line:

when the central axis of the cavitation excitation ultrasonic transducer is intersected with the imaging plane of the array ultrasonic transducer and is vertical to the array element direction of the array ultrasonic transducer,

exciting the ultrasonic transducer surface center coordinates (x) for the pixel coordinates (x, z) to cavitation _E ,y _E ,z _E ) Distance of tangential plane:

wherein alpha is an included angle between the central axis of the cavitation excitation ultrasonic transducer and the central axis of the array ultrasonic transducer, and |·| represents an absolute value.

Preferably, in the step 1, the time length of the time gating interval is the pulse length of the pulse emitted by the cavitation excitation ultrasonic transducer.

Preferably, in the step 1, the time gating signal of each array element is expressed as:

wherein i=1, 2,.. _SE ，k′＝1,2,...,TGS ^TL ，rf _i,k For the array ultrasonic transducer to passively receive the signal value of the ith array element at the kth sampling point in the obtained frame of original radio frequency signal,

for the sampling point number corresponding to the initial time of the time gating interval of the ith array element, TGS ^TL The sampling point number is corresponding to the time length of the time gating interval.

Preferably, in the step 2 or the step 3, the radius of the closed neighborhood is set as the unit of the number of pixels

The number of pixels in the close neighborhood threshold is set to 2.

Preferably, in the step 2 or the step 3, the condition for initializing the pixel cluster is: and selecting an undetermined pixel from the global pixel set or the local pixel set, if the number of pixels in the closed neighborhood of the selected pixel is smaller than the threshold value of the number of pixels in the closed neighborhood, reselecting the pixel until the number of pixels in the closed neighborhood of the selected pixel is larger than or equal to the threshold value of the number of pixels in the closed neighborhood, wherein the selected pixel is a core pixel at the moment, and initializing a pixel cluster.

Preferably, in the step 2 or the step 3, the expanding pixel cluster specifically includes the following steps:

2.1 -aggregating pixels selected from a global or local set of pixels and determined as core pixels into said cluster of pixels;

2.2 Classifying all pixels in the heart-removing closed neighbor of the core pixel into a pixel set to be judged;

2.3 Optionally selecting an undetermined pixel in the pixel set to be judged, classifying the undetermined pixel in the pixel heart-removing closed neighbor of the selected pixel into the pixel set to be judged if the selected pixel is a core pixel, and converting to the step 2.4; if the selected pixel is not the core pixel, directly turning to step 2.4;

2.4 If the pixel selected in the step 2.3 is not clustered, the pixel is clustered into the pixel cluster;

2.5 Repeating the steps 2.3 and 2.4 until all pixels in the pixel set to be judged are traversed, so as to generate a global pixel cluster or a local pixel cluster.

Preferably, in the step 3, the center of the local cluster window is a pixel corresponding to the centroid of the global pixel cluster, and the size of the local cluster window is determined according to the point spread function.

Preferably, the threshold coefficient used in the thresholding in step 3 is smaller than the threshold coefficient used in the thresholding in step 2.

Preferably, in the step 3, the searching the local pixel main cluster from the local pixel clusters specifically includes the following steps: and setting the label pixel as a pixel corresponding to the mass center of the global pixel cluster, and setting a certain local pixel cluster containing the label pixel as a local pixel main cluster.

Preferably, in the step 4, the pixel peak threshold is 0.4 to 0.6 times of the maximum value of the pixel peaks of all the single cavitation source energy maps.

Preferably, in the step 5, the gaussian distribution function allows angular deflection, and the expression is:

wherein A is Gaussian distribution peak value, mu _x Sum mu _z Gaussian distribution peak position, sigma, in x-axis and z-axis directions, respectively _x Sum sigma _z The standard deviation of the Gaussian distribution in the x-axis direction and the z-axis direction are respectively shown, and θ is the Gaussian distribution deflection angle.

Preferably, in the step 5, the gaussian fitting needs to initialize parameters in a gaussian distribution function: initializing a Gaussian distribution peak value to be the maximum pixel value in the single cavitation source energy diagram, initializing the position of the Gaussian distribution peak value to be the coordinate of the maximum pixel value in the single cavitation source energy diagram, and initializing the Gaussian distribution standard deviation according to a point spread function.

Preferably, in the step 5, the gaussian fitting uses a nonlinear least square method based on a Levenberg-Marquardt algorithm.

Preferably, in the step 6, the threshold values of uncertainty and goodness of fit of the cavitation source coordinate estimation are respectively:

wherein,,

and->

Uncertainty delta is estimated for cavitation source coordinates in the x-axis and z-axis directions, respectively _x And delta _z Threshold of->

Threshold for goodness of fit gof, +.>

Uncertainty set delta is estimated for cavitation source coordinates in the x-axis direction and the z-axis direction, respectively _x 、Δ _z The upper quartile of the elements in the set obtained after the ascending arrangement of the elements in (a), the ++>

Uncertainty set delta is estimated for cavitation source coordinates in the x-axis direction and the z-axis direction, respectively _x 、Δ _z The elements in the set obtained after the ascending arrangement of the elements in the set are quartile range, the ++>

And IQR _GOF The lower quartile and the quartile of the elements in the set obtained after the elements in the goo set GOF are arranged in ascending order are respectively.

Preferably, in the step 6, detecting the cavitation source of normal positioning specifically includes the steps of: if cavitation source coordinates in x-axis direction and z-axis direction estimate uncertainty delta _x And delta _z The goodness of fit gof is simultaneously satisfied

And->

And is also provided with

And if not, normal positioning is performed, otherwise abnormal positioning is performed.

Preferably, in the step 7, the number of pixels of the pixel grid of the super-resolution imaging is λ of the number of pixels of the pixel grid of the time-gated energy map _x ×λ _z Multiple of lambda _x And lambda (lambda) _z The values of the multiples are respectively equal to or more than 10 in the x-axis direction and the z-axis direction.

Preferably, in the step 7, the positioning distribution function is a gaussian distribution function allowing angular deflection, and the peak position, standard deviation, peak value and deflection angle of the gaussian distribution function are respectively estimated values of cavitation source coordinates in the x-axis direction and the z-axis direction, estimated uncertainty of cavitation source coordinates in the x-axis direction and the z-axis direction, estimated values of gaussian distribution peak value and estimated values of gaussian distribution deflection angle obtained by performing gaussian fitting on the single cavitation source energy map.

The system comprises a time gating energy map calculation module, a single cavitation source separation module, a single cavitation source energy image element peak value judgment module, a Gaussian fitting cavitation source positioning module, a cavitation source positioning detection module and a cavitation source positioning distribution map calculation and processing module;

the time gating energy map calculation module is used for executing the step 1, and is mainly used for establishing an imaging coordinate system and planning a pixel grid of a time gating energy map, calculating double-pass delay according to the distance from each pixel coordinate to a tangent line or a tangent plane at the center coordinate of the surface of the cavitation excitation ultrasonic transducer in the pixel grid and the distance between each pixel coordinate and each array element coordinate of the array ultrasonic transducer, setting the time length of the time gating interval by taking the double-pass delay as the initial time of the time gating interval, extracting the time gating signals of each array element from each frame of original radio frequency signals passively received by the array ultrasonic transducer according to the time gating interval of each array element, superposing the time gating signals of each array element along the array element direction, and superposing the squares of the superposed time gating synthesized signals along the sampling point direction, thereby obtaining the time gating energy map corresponding to each frame of original radio frequency signals;

The single cavitation source separation module is used for executing the step 2 and the step 3, and specifically comprises a pixel global clustering sub-module and a pixel local clustering sub-module and a local pixel main cluster searching sub-module;

the pixel global clustering submodule is used for executing the step 2, and is mainly used for establishing a global pixel set according to all non-zero value pixels after thresholding the time-gating energy map of each frame obtained by the time-gating energy map calculation module, and clustering pixels in the global pixel set by repeating the processes of initializing the pixel clusters and expanding the pixel clusters according to the set closed neighborhood radius and the set closed neighborhood pixel number threshold value so as to obtain a global pixel cluster of the time-gating energy map of each frame;

the local clustering of pixels and searching for a local pixel main cluster submodule are used for executing the step 3, and are mainly used for defining a local clustering window for each global pixel cluster of each frame time gating energy image obtained by the pixel global clustering submodule, carrying out frame selection on the time gating energy image by utilizing the local clustering window, carrying out thresholding treatment on the local energy image obtained by frame selection, establishing a local pixel set according to all non-zero value pixels, carrying out clustering on pixels in the local pixel set by repeating the processes of initializing the pixel clusters and expanding the pixel clusters, and searching for a local pixel main cluster from the local pixel clusters obtained by the clustering and setting the pixel values of all pixels which are not contained in the local pixel main cluster in the local energy image to be zero so as to obtain a single cavitation source energy image;

The single cavitation source energy image peak value judging module is used for executing the step 4, and is mainly used for extracting pixel peak values of all the single cavitation source energy images obtained by the local pixel clustering and local pixel main cluster searching submodule, judging the pixel peak value of each single cavitation source energy image according to the set pixel peak value threshold value, and establishing a single cavitation source energy image set by the single cavitation source energy images with the pixel peak value smaller than the pixel peak value threshold value;

the Gaussian fitting cavitation source positioning module is used for executing the step 5, and is mainly used for establishing a Gaussian distribution function in an imaging coordinate system, carrying out Gaussian fitting on each single cavitation source energy image in a single cavitation source energy image peak value judgment module to position a single cavitation source, substituting an estimated value of a parameter in the Gaussian distribution function obtained by fitting into the Gaussian distribution function to calculate a Gaussian fitting result, calculating the uncertainty of cavitation source coordinate estimation under a given confidence level, and calculating the fitting goodness according to the single cavitation source energy image and the Gaussian fitting result of the single cavitation source energy image; in the estimated values of the parameters, the estimated value of the peak position of the Gaussian distribution is the estimated value of the cavitation source coordinates;

The cavitation source positioning detection module is used for executing the step 6, and is mainly used for establishing a cavitation source coordinate estimation uncertainty set and a fitting goodness set according to the cavitation source coordinate estimation uncertainty and the fitting goodness of all the cavitation sources obtained by the Gaussian fitting cavitation source positioning module, respectively carrying out ascending arrangement on elements in each set, calculating quartiles and quartiles of the elements in each set after ascending arrangement, setting a threshold value of the cavitation source coordinate estimation uncertainty and a threshold value of the fitting goodness according to the quartiles and the quartiles, and detecting the cavitation source in normal positioning by utilizing the set threshold value of the cavitation source coordinate estimation uncertainty and the threshold value of the fitting goodness;

the cavitation source positioning distribution map calculating and processing module is used for executing the step 7, and is mainly used for establishing a positioning distribution function for each normal positioning cavitation source obtained by the cavitation source positioning detection module, calculating a positioning distribution function value of each pixel coordinate in a planned super-resolution imaging pixel grid, and carrying out normalization and logarithmic processing after overlapping the calculated positioning distribution maps of all the normal positioning cavitation sources so as to obtain a super-resolution imaging result.

The beneficial effects of the invention are as follows:

according to the super-resolution ultrasonic passive cavitation imaging method, the problems that ultrasonic passive cavitation imaging resolution is limited by an array ultrasonic transducer diffraction mode are solved through time gating energy map calculation, single cavitation source separation based on pixel global clustering and pixel local clustering and local pixel main cluster searching, single cavitation source energy image peak value judgment, gaussian fitting cavitation source positioning, cavitation source positioning detection and cavitation source positioning distribution map calculation and processing, and the ultrasonic passive cavitation imaging resolution is greatly improved. The invention can carry out fine analysis on the spatial distribution of cavitation sources in a small area in a micron scale, provides an important technical means for researching cavitation physical phenomenon and cavitation-mediated various therapeutic ultrasonic applications (such as tissue damage, drug delivery and the like), and lays a foundation for realizing accurate diagnosis and treatment of ultrasonic waves.

Furthermore, the method and the device have the advantages that the time gating interval with the initial time of double-pass delay and the time length of the pulse emitted by the cavitation excitation ultrasonic transducer is set, the time gating signal is extracted from the original radio frequency signal, then the time gating signal is processed to obtain the time gating energy diagram, the interference artifact in the axial direction of the array ultrasonic transducer is effectively restrained, and a high-quality original image is provided for the subsequent pixel clustering cavitation source separation and Gaussian fitting cavitation source positioning processing, so that the positioning accuracy of the cavitation source is improved.

Further, the cavitation source contained in the time-gating energy diagram is primarily identified by thresholding the time-gating energy diagram (the threshold coefficient is larger) and globally clustering the non-zero value pixels; then, a local clustering window is defined, a time gating energy diagram is subjected to frame selection, the local energy diagram is subjected to thresholding (the threshold coefficient is smaller), non-zero value pixels are subjected to local clustering, and a local pixel main cluster is searched, so that the area of a single cavitation source is accurately separated; this makes it possible to locate a single cavitation source directly, thereby improving the accuracy of cavitation source location.

Further, the pixel peak value of the single cavitation source energy map is judged by setting the pixel peak value threshold, and the single cavitation source energy map with the pixel peak value larger than or equal to the pixel peak value threshold is removed, so that the problem of inaccurate cavitation source positioning caused by cavitation source adhesion is solved, and the accuracy of super-resolution ultrasonic passive cavitation imaging is improved.

Further, the invention obtains the estimated value of the cavitation source coordinates and the estimated uncertainty of the cavitation source coordinates by carrying out Gaussian fitting on the single cavitation source energy diagram, and can realize the direct positioning of the single cavitation source; the Gaussian fitting is performed by adopting a Gaussian distribution function allowing angle deflection, and the function has universality for energy beams deflected randomly under different conditions and is beneficial to improving the positioning accuracy of cavitation sources.

Further, the uncertainty and the goodness of fit of the cavitation source coordinate estimation are adopted as evaluation indexes of Gaussian fitting cavitation source positioning, the threshold value of the uncertainty and the threshold value of the goodness of fit of the cavitation source coordinate estimation are respectively set according to the quartile and the quartile distance, and the threshold value of the uncertainty and the threshold value of the goodness of fit of the cavitation source coordinate estimation are combined to detect whether the positioning of the cavitation source is normal or not, so that the abnormally positioned cavitation source is removed, and the accuracy of super-resolution ultrasonic passive cavitation imaging is improved.

Furthermore, the method establishes a positioning distribution function aiming at a cavitation source which is positioned normally, calculates a positioning distribution diagram by using a high-density pixel grid, and obtains a super-resolution imaging result by superposition on the basis, thereby carrying out fine analysis on the spatial distribution of the cavitation source on a micron scale.

Drawings

FIG. 1 is a flow chart of time-gated energy graph computation in an embodiment of the present invention.

FIG. 2 is a flow chart of global clustering of pixels in an embodiment of the present invention.

FIG. 3 is a flowchart of local clustering of pixels and finding a local pixel cluster in an embodiment of the present invention.

FIG. 4 is a flow chart of single cavitation source energy image pixel peak judgment and single cavitation source energy image set up in an embodiment of the invention.

FIG. 5 is a flow chart of Gaussian fitting cavitation source localization in an embodiment of the invention.

FIG. 6 is a flowchart of the hollow source location detection according to an embodiment of the present invention.

FIG. 7 is a flowchart illustrating the calculation and processing of a localization distribution map of a localization source according to an embodiment of the present invention.

FIG. 8 shows ultrasonic passive cavitation imaging results under three cavitation source distributions of random distribution (a, d), straight distribution (b, e) and curved distribution (c, f) in the embodiment of the invention; the method comprises the steps of a, b and c, wherein a traditional imaging method based on time delay superposition integration is used, and d, e and f use a super-resolution imaging method based on single cavitation source separation and positioning.

Detailed Description

The invention is described in further detail below with reference to the drawings and examples. The examples are given solely for the purpose of illustration and are not intended to limit the scope of the invention.

The invention provides a super-resolution ultrasonic passive cavitation imaging method based on single cavitation source separation and positioning, which specifically comprises the following steps:

referring to fig. 1, a time gating interval is obtained by calculating a two-way delay and setting a time length of the time gating interval, a time gating signal of each array element is extracted from an original radio frequency signal, the signal is superimposed along the array element direction, and the square of the obtained time gating composite signal is superimposed along the sampling point direction, thereby obtaining a time gating energy diagram. The specific flow is as follows steps (1.1) - (1.11);

(1.1) establishing an imaging coordinate system by taking the center of an array ultrasonic transducer (such as a linear array transducer) as an origin, taking the array element direction of the array ultrasonic transducer as an x-axis and taking the central axis direction of the array ultrasonic transducer as a z-axis; on the basis of an imaging coordinate system, a three-dimensional coordinate system is established by taking a direction perpendicular to an imaging plane of the array ultrasonic transducer as a y-axis;

(1.2) planning a grid of pixels of the time-gated energy map in the imaging coordinate system of step (1.1), the grid being determined by imaging ranges in the x-axis direction and the z-axis direction (e.g., -4 mm and 36-44 mm) and pixel spacings in the x-axis direction and the z-axis direction (e.g., 0.1mm and 0.1 mm);

(1.3) for any one of the pixel coordinates (x, z) in the pixel grid in step (1.2), calculating the distance of the pixel coordinate to the tangent or tangent plane at the center coordinate of the cavitation-excited ultrasonic transducer surface

The coordinates of the surface center of the cavitation excitation ultrasonic transducer in the three-dimensional coordinate system are (x _E ,y _E ,z _E )；

When the central axis of the cavitation excitation ultrasonic transducer is in the imaging plane of the array ultrasonic transducer, the distance from the pixel coordinate (x, z) to the tangent line at the central coordinate of the surface of the cavitation excitation ultrasonic transducer

The calculation formula of (2) is as follows:

when the central axis of the cavitation excitation ultrasonic transducer intersects with the imaging plane of the array ultrasonic transducer and is perpendicular to the array element direction of the array ultrasonic transducer, the distance from the pixel coordinate (x, z) to the tangential plane at the central coordinate of the surface of the cavitation excitation ultrasonic transducer

The calculation formula of (2) is as follows:

wherein x is _E 、y _E And z _E Respectively, an x-axis coordinate of the surface center of the cavitation excitation ultrasonic transducer in the three-dimensional coordinate system in the step (1.1),The y-axis coordinate and the z-axis coordinate, alpha is an included angle (for example, 0-pi/2) between the central axis of the cavitation excitation ultrasonic transducer and the central axis of the array ultrasonic transducer, and |·| represents an absolute value;

(1.4) calculating the pixel coordinates (x, z) in step (1.3) and the i-th array element coordinates (x) of the array ultrasonic transducer _i Distance between 0)

Wherein i=1, 2, …, N _SE ，N _SE For the number of array elements (e.g., 128), (x) _i 0) the coordinates of the ith array element of the array ultrasonic transducer under the imaging coordinate system;

(1.5) the distance obtained according to steps (1.3) and (1.4)

And->

Calculating the double-pass delay tau _i (x,z)：

Wherein c is the sound propagation velocity;

(1.6) setting the initial time of the time gating interval of the ith array element of the array ultrasonic transducer as the double-pass delay tau obtained in the step (1.5) _i (x, z) setting the time length of the time gating interval of each array element of the array ultrasonic transducer as the pulse length tau of the pulse emitted by the cavitation excitation ultrasonic transducer _TL Thereby obtaining the time gating interval of the ith array element of the array ultrasonic transducer;

(1.7) passively connecting the array ultrasonic transducer according to the time gating interval of the ith array element obtained in the step (1.6) Extracting the time gating signal of the ith array element from the received frame of original radio frequency signal

Wherein k' =1, 2, …, TGS ^TL ，rf _i,k For the signal value of the ith array element at the kth sampling point in the original radio frequency signal,

for the number of sampling points corresponding to the initial time of the time gating interval of the ith array element,

TGS ^TL for the sampling point number corresponding to the time length of the time gating interval, TGS ^TL ＝round(τ _TL Fs)，round[·]Representing a rounding operation, fs is the sampling frequency of the original radio frequency signal (e.g., 50 MHz);

(1.8) superposing the time gating signals of all the array elements obtained in the step (1.7) along the array element direction to obtain a time gating composite signal

(1.9) time-gating the resulting time-gated composite signal of step (1.8)

Is superimposed along the direction of the sampling point to obtain a time-gating energy value E at the pixel coordinate (x, z) ^TG (x,z)：

(1.10) repeating the steps (1.3) - (1.9) until the time gating energy values of all the pixel coordinates are calculated, and forming a time gating energy diagram by the time gating energy values of all the pixel coordinates (namely, the pixel values at each pixel coordinate of the time gating energy diagram are corresponding time gating energy values);

(1.11) transmitting N to an air-excited ultrasound transducer _F Passive reception of the resulting N by an array ultrasound transducer under a pulse _F Repeating the calculation process of the frame original radio frequency signal to obtain N _F The frame time gates the energy map.

Referring to fig. 2, thresholding is performed on each frame of time-gated energy map, a global pixel set is established, a close-neighborhood radius and a close-neighborhood pixel number threshold are set, and a global pixel cluster of each frame of time-gated energy map is obtained by repeating the processes of initializing pixel clusters and expanding pixel clusters. The specific flow is as follows steps (2.1) - (2.7);

(2.1) setting a threshold coefficient γ _G (0＜γ _G < 1, for example, 0.6) and according to this threshold coefficient, for the p-th frame (p=1, 2,) obtained in step (1.11), N _F ) Thresholding the time-gated energy map, i.e. using gamma, which is the maximum pixel value of the frame time-gated energy map _G Multiplying the pixel value by a decision value, and setting the pixel value smaller than the decision value in the time gating energy diagram to be zero;

(2.2) establishing a global pixel set S from all non-zero value pixels obtained in step (2.1) _G ，S _G Is the index (u) of these non-zero value pixels in the p-th frame time-gated energy map _i ,v _i ) Wherein u is _i And v _i Index in x-axis direction and z-axis direction, i=1, 2, …, N _PG ，N _PG The number of non-zero value pixels is denoted by an index below;

(2.3) setting the radius epsilon of the closed neighborhood (in pixels, for example,

) And a close neighborhood pixel number threshold β (e.g E.g., 2), wherein the closed-neighborhood refers to the pixel (u) _i ,v _i ) A circular area centered around ε and having a radius, denoted +.>

(2.4) initializing a global pixel cluster:

(2.4.1) the global pixel set S from step (2.2) _G Optionally one undetermined pixel (u) _i ,v _i )；

(2.4.2) if the pixel (u) _i ,v _i ) Is a closed neighborhood of (2)

If the number of the pixels is smaller than beta, the pixels are not core pixels, and the step (2.4.1) is returned; if the pixel (u _i ,v _i ) Is->

The number of the pixels is greater than or equal to beta, the pixel is taken as a core pixel and is marked as +.>

Initializing a pixel cluster C and marking the heart-removed closed neighborhood of the core pixel (i.e. the core pixel is removed from the closed neighborhood) as +.>

(2.5) expanding to obtain a global pixel cluster:

(2.5.1) the core pixel obtained in step (2.4.2)

Into the pixel cluster C;

(2.5.2) the core pixel obtained in step (2.4.2)

Is a heart-removing closure neighborhood->

All pixels in the set S are classified into a pixel set S to be judged _J In (a) and (b);

(2.5.3) at the pixel set S to be judged _J Optionally one undetermined pixel (u) _j ,v _j ) Judging whether the pixel is a core pixel according to the method in the step (2.4.2); if yes, removing the heart from the neighborhood

The pixels which are not judged in the set S are classified into the set S of pixels to be judged _J And jump to step (2.5.4); if not, directly jumping to the step (2.5.4);

(2.5.4) if the pixel (u) in the step (2.5.3) _j ,v _j ) Not clustered, then clustered into the pixel cluster C;

(2.5.5) repeating the steps (2.5.3) and (2.5.4) until the pixel set S to be judged is traversed _J Thereby generating a global pixel cluster;

(2.6) repeating steps (2.4) and (2.5) until the global pixel set S obtained in step (2.2) is traversed _G To generate a p-th frame time-gated energy map

Global pixel clusters;

(2.7) repeating the steps (2.1) to (2.6) to obtain N _F N of frame time gated energy map _ASC A global pixel cluster, wherein

Referring to fig. 3, a local clustering window is defined for each global pixel cluster of each frame time gating energy map, the local energy map selected by the window frame is subjected to thresholding processing, a local pixel set is established, pixels in the local pixel set are clustered to obtain a local pixel cluster, and a local pixel main cluster is searched from the local pixel cluster, so that a single cavitation source energy map is obtained. The specific flow is as follows steps (3.1) - (3.9);

(3.1) the qth global pixel cluster (p=1, 2, …, N) of the energy map for the p-th frame time gating _F ，q＝1,2,…,

) Defining a local cluster window; the center of the local clustering window is a pixel corresponding to the centroid of the global pixel cluster, and the index of the pixel is +. >

Wherein u is _i 、v _i And E is _i Respectively representing indexes and pixel values of pixels in the global pixel cluster in the x-axis direction and the z-axis direction, wherein round (·) represents rounding operation; the size of the local cluster window is determined from the point spread function (e.g., 5 x 5);

(3.2) framing the p-th frame time-gated energy map with the local cluster window obtained in step (3.1) to obtain a local energy map containing pixels corresponding to the centroids of the global pixel clusters and having a determined size (e.g., 5×5);

(3.3) setting a threshold coefficient γ _L (0＜γ _L ＜γ _G For example, 0.3) and thresholding the local energy map obtained in step (3.2) according to the threshold coefficient, i.e. with gamma of the maximum pixel value of the local energy map _L Multiplying the local energy map by a decision value, and setting a pixel value smaller than the decision value in the local energy map to be zero;

(3.4) creating a local pixel set S from all non-zero value pixels obtained in step (3.3) _L ，S _L Is the index of these non-zero value pixels in the time-gated energy map;

(3.5) the local pixel set S obtained in step (3.4) is subjected to the method in steps (2.4) to (2.6) _L Clustering the pixels in the array to obtain a plurality of local pixel clusters;

(3.6) setting the label pixel as the pixel corresponding to the mass center of the global pixel cluster in the step (3.1), and searching a local pixel main cluster from a plurality of local pixel clusters obtained in the step (3.5) according to the label pixel; judging whether the label pixel is in a local pixel cluster or not, if so, the local pixel cluster is a local pixel main cluster;

(3.7) setting the pixel values of all pixels which are not included in the local pixel main cluster obtained in the step (3.6) in the local energy map obtained in the step (3.2) to be zero, so as to obtain a single cavitation source energy map, and further realize separation of single cavitation sources;

(3.8) repeating the steps (3.1) - (3.7) to obtain the p-th frame time-gating energy diagram

A single cavitation source energy map;

(3.9) repeating the steps (3.1) to (3.8) to obtain N _F N of frame time gated energy map _ASC Single cavitation source energy map.

Referring to fig. 4, the pixel peak values of all the single cavitation source energy maps are extracted, the pixel peak value threshold is set, the pixel peak value of each single cavitation source energy map is judged, and a single cavitation source energy map set is built by the single cavitation source energy maps with the pixel peak values smaller than the pixel peak value threshold. The specific flow is as follows steps (4.1) - (4.5);

(4.1) extracting the pixel peaks of all the energy maps of the single cavitation source, namely extracting N obtained in the step (3.9) respectively _ASC Maximum pixel value of each single cavitation source energy diagram in the single cavitation source energy diagrams;

(4.2) setting a pixel peak threshold, for example, the threshold being set to 0.5 times the maximum value of the pixel peaks of all the single cavitation source energy maps obtained in step (4.1);

(4.3) initializing a single cavitation source energy diagram set to be an empty set;

(4.4) N obtained in the step (3.9) _ASC The mth single cavitation source energy map of the single cavitation source energy maps (m=1, 2,., N _ASC ) If the pixel peak value of the single cavitation source energy diagram is smaller than the pixel peak value threshold value set in the step (4.2), adding the single cavitation source energy diagram into a single cavitation source energy diagram set;

(4.5) repeating step (4.4) until for N _ASC The energy diagrams of the single cavitation sources are all judged, so that a single cavitation source energy diagram collection is established and obtained, and the energy diagrams of the single cavitation sources are collectedThe number of elements of the set is noted as

Referring to fig. 5, a gaussian distribution function allowing angular deflection is established in an imaging coordinate system, gaussian fitting is performed on each single cavitation source energy map in a single cavitation source energy map set, and uncertainty of cavitation source coordinate estimation and goodness of fit are calculated. The specific flow is as follows steps (5.1) - (5.7);

(5.1) establishing a gaussian distribution function allowing angular deflection in the imaging coordinate system of step (1.1):

wherein A is Gaussian distribution peak value, mu _x Sum mu _z Gaussian distribution peak position, sigma, in x-axis and z-axis directions, respectively _x Sum sigma _z The standard deviation of the Gaussian distribution in the x-axis direction and the z-axis direction respectively, and θ is the deflection angle of the Gaussian distribution, and the parameters A and μ are the above _x 、μ _z 、σ _x 、σ _z θ are parameters to be estimated;

(5.2) SCE for the nth single cavitation source energy map in the set of single cavitation source energy maps obtained in step (4.5) _n (n＝1,2,...,

) Initializing parameters in the Gaussian distribution function of step (5.1), wherein A is initialized to the maximum pixel value in the single cavitation source energy map, μ _x Sum mu _z Initializing sigma according to a point spread function by initializing to the coordinate of the maximum pixel value in the single cavitation source energy diagram _x Sum sigma _z (e.g., 0.1mm for both), θ is initialized to 0;

(5.3) mapping SCE of the nth single cavitation source energy according to the Gaussian distribution function of the step (5.1) and the initialization value of each parameter obtained in the step (5.2) _n Performing Gaussian fitting, wherein the fitting process adopts a method based on levelNonlinear least squares method of berg-Marquardt algorithm to obtain the estimated value of the parameter of step (5.1)

And

wherein->

And->

The estimated values of the cavitation source coordinates in the x-axis direction and the z-axis direction are respectively, namely, a single cavitation source is positioned through Gaussian fitting;

(5.4) substituting the estimated values of the parameters obtained in the step (5.3) into the Gaussian distribution function obtained in the step (5.1), and calculating to obtain a single cavitation source energy map SCE _n Gaussian fitting result SCEF of (2) _n ；

(5.5) setting a confidence level (e.g., 0.68) based on confidence intervals for cavitation source coordinate estimation in the x-axis and z-axis directions at the confidence level

And->

Calculating cavitation source coordinate estimation uncertainty delta in the x-axis direction and the z-axis direction of the single cavitation source of step (5.3) respectively _x And delta _z ：

(5.6) the single cavitation source energy map SCE according to step (5.2) _n And the Gaussian fitting result SCEF obtained in the step (5.4) _n Calculating a goodness of fit gof:

wherein, SCE _n (x, z) and SCEF _n (x, z) are SCE respectively _n And SCEF _n A pixel value at any one pixel coordinate (x, z),

is SCE (SCE) _n An average of pixel values at all pixel coordinates;

(5.7) repeating steps (5.2) to (5.6) until the energy map of the single cavitation source is gathered

And (3) finishing the processing of the single cavitation source energy map (namely Gaussian fitting, calculating uncertainty of cavitation source coordinate estimation and goodness of fit).

Referring to fig. 6, a cavitation source coordinate estimation uncertainty set and a fitting goodness set are established, thresholds of the cavitation source coordinate estimation uncertainty and the fitting goodness are set according to quartiles and quartiles of elements in the set, and the obtained thresholds are used for detecting the positioning of the cavitation source in the time-gated energy map per frame. The specific flow is as follows steps (6.1) - (6.7);

(6.1) obtained in step (5.7)

Uncertainty delta is estimated from cavitation source coordinates in the x-axis and z-axis directions of each cavitation source _x And delta _z And the goodness of fit gof is taken as an element, and cavitation source coordinate estimation uncertainty sets delta in the x-axis direction and the z-axis direction are respectively established _x And delta _z And a goodness-of-fit set GOF, the number of elements of these sets being +.>

(6.2) for the set delta obtained in step (6.1) _x 、Δ _z The elements in GOF are arranged in ascending order respectively, and the set after ascending order is still recorded as delta for convenience _x 、Δ _z 、GOF；

(6.3) calculating the set delta after the ascending arrangement obtained in the step (6.2) _x 、Δ _z Quartiles of elements in GOF:

wherein M represents the number of elements

Is specifically set delta _x 、Δ _z Or GOF, Q _D And Q _U Respectively represent the set delta _x 、Δ _z Or the lower quartile and the upper quartile, b, of an element in the GOF _D And b _U Respectively represent the set delta _x 、Δ _z Or the lower quartile position and the upper quartile position of the element in the GOF, ++>

floor (-) represents a downward rounding and ceil (-) represents an upward rounding;

let M in the above two formulas be delta _x Calculating to obtain a set delta _x Lower quartile of middle element

And upper quartile

Let M in the above two formulas be delta _z Calculating to obtain a set delta _z Lower quartile of the middle element->

And upper quartile->

Let M in the above two formulas be GOF, calculate and get the lower quartile of the element in the set GOF ++>

And upper quartile->

(6.4) calculating the sets Deltarespectively based on the results obtained in the step (6.3) _x 、Δ _z Quartile range of elements in GOF

IQR _GOF ：

/>

(6.5) setting the cavitation source coordinate estimation uncertainty delta according to the quartile obtained in the step (6.3) and the quartile range obtained in the step (6.4) respectively _x And delta _z Threshold of (2)

And->

And threshold value of goodness of fit gof +.>

(6.6) positioning of the nth cavitation source (n=1, 2.),

) Detecting if the cavitation source coordinates in the x-axis direction and the z-axis direction of the cavitation source estimate uncertainty delta _x And delta _z The goodness of fit gof is satisfied simultaneously +.>

And->

And->

Then the normal positioning is considered, otherwise the abnormal positioning is considered;

(6.7) repeating the step (6.6) until the detection is completed

Positioning of cavitation sources, the number of cavitation sources positioned normally is recorded as +.>

Referring to fig. 7, a positioning distribution function is established for the cavitation sources which are positioned normally, a positioning distribution diagram is calculated, and normalization and logarithmization processing are performed after the positioning distribution diagrams of all the cavitation sources positioned normally are overlapped, so that a super-resolution imaging result is obtained. The specific flow is as follows steps (7.1) - (7.4);

(7.1) planning a grid of pixels of the super-resolution imaging in the imaging coordinate system of step (1.1), the number of pixels of the grid of pixels being the image of the grid of pixels of step (1.2)Lambda of prime number _x ×λ _z Multiple of lambda _x And lambda (lambda) _z Multiples of x-axis and z-axis directions, respectively (e.g., 16 for each);

(7.2) for step (6.7)

The s-th cavitation source of the plurality of normally located cavitation sources (s=1, 2,>

) Establishing a positioning distribution function LD (x, z):

The function is a Gaussian distribution function allowing angular deflection, and the peak position, the peak value and the deflection angle are estimated values of cavitation source coordinates in the x direction and the z direction obtained by carrying out Gaussian fitting on the corresponding single cavitation source energy diagram

And->

Estimated value of Gaussian distribution peak>

And an estimated value of the gaussian distribution deflection angle +.>

The standard deviation is the uncertainty delta of the cavitation source coordinate estimation of the cavitation source in the x-axis direction and the z-axis direction _x And delta _z ；

(7.3) calculating the positioning distribution function value of each pixel coordinate in the pixel grid obtained in the step (7.1) according to the positioning distribution function obtained in the step (7.2), and obtaining the positioning distribution diagram of the s < th > normal positioning cavitation source;

(7.4) repeating the steps (7.2) to (7.3) to obtainTo the point of

A localization profile of a normally localized cavitation source for the resultant +.>

And superposing the positioning distribution diagrams, and then sequentially carrying out normalization and logarithmization treatment to obtain a super-resolution imaging result.

The performance of the super-resolution ultrasonic passive cavitation imaging method is tested by adopting a numerical simulation method, and three simulated cavitation source distributions are adopted: random distribution, linear distribution, and curvilinear distribution. Fig. 8a, 8b and 8c are respectively ultrasonic passive cavitation imaging results (shown in a linear scale) obtained by using a traditional imaging method based on time-lapse superposition integration under random distribution, linear distribution and curve distribution, and fig. 8d, 8e and 8f are respectively ultrasonic passive cavitation imaging results (shown in a logarithmic scale) obtained by using a super-resolution imaging method based on single cavitation source separation and positioning proposed by the present invention under random distribution, linear distribution and curve distribution. 8a, 8b and 8c, a plurality of bright spots with large size appear in the imaging results, the imaging resolution is poor, and the cavitation source spatial distribution information cannot be effectively reflected; in contrast, the spatial distribution of the cavitation sources in the imaging results shown in fig. 8d, 8e and 8f is resolved finely at a very small scale (up to micron level), the imaging resolution is greatly improved, and the cavitation source spatial distribution information is effectively reflected.

The invention has the following advantages:

(1) The problem of poor spatial resolution exists in ultrasonic passive cavitation imaging, and although various improved imaging methods developed at present can improve the image quality to a certain extent, the limitation of the diffraction mode of an imaging transducer on the imaging resolution cannot be fundamentally overcome. According to the method, super-resolution imaging of space distribution of an cavitation source is achieved through time gating energy map calculation, single cavitation source separation based on pixel global clustering and pixel local clustering and local pixel main cluster searching, single cavitation source energy image peak value judgment, gaussian fitting cavitation source positioning, cavitation source positioning detection and cavitation source positioning distribution map calculation and processing, the problem that imaging resolution is limited by a diffraction mode of a transducer is fundamentally solved, and the resolution of ultrasonic passive cavitation imaging is greatly improved.

(2) High-level interference artifacts along the axial direction of an array ultrasonic transducer can appear in an image generated by a traditional ultrasonic passive cavitation imaging method, and a larger positioning error can be caused by using the image to perform Gaussian fitting-based positioning on a cavitation source. The method extracts the signals in the time gating interval with the initial time of double-pass delay and the time length of the pulse emitted by the cavitation excitation ultrasonic transducer, effectively suppresses interference artifacts in the calculated time gating energy diagram, and provides high-quality original images for subsequent pixel clustering cavitation source separation and Gaussian fitting cavitation source positioning treatment, thereby improving the positioning precision of the cavitation sources.

(3) The method comprises the steps of firstly, primarily identifying cavitation sources contained in a time-gating energy map through a high-thresholding time-gating energy map and a pixel global clustering, and then accurately separating a region where a single cavitation source is located through defining a local clustering window frame time-selecting gating energy map, a low-thresholding local energy map, a pixel local clustering and searching a local pixel main cluster, so that the single cavitation source is allowed to be positioned directly, and the positioning precision of the cavitation source is improved; the invention further solves the problem of inaccurate positioning of the cavitation source caused by the adhesion of the cavitation source by judging the pixel peak value of the single cavitation source energy diagram, thereby improving the accuracy of super-resolution ultrasonic passive cavitation imaging.

(4) The Gaussian distribution function allowing angular deflection is adopted when the Gaussian fitting is carried out on the single cavitation source energy diagram, and has universality on arbitrarily deflected energy beams under different situations (for example, cavitation excitation ultrasonic transducers and array ultrasonic transducers are placed in different modes and the like), so that the improvement of the positioning accuracy of the cavitation source is facilitated.

(5) According to the invention, uncertainty and goodness of fit of the cavitation source coordinate estimation are adopted to evaluate Gaussian fit cavitation source positioning, thresholds of the two indexes are set based on a quartile method, abnormal positioning cavitation sources are detected, and abnormal positioning cavitation sources are removed when a cavitation source positioning distribution map is overlapped (namely, only a positioning distribution map of normal positioning cavitation sources is overlapped), so that the accuracy of super-resolution ultrasonic passive cavitation imaging is improved.

(6) The method establishes a positioning distribution function for a cavitation source which is positioned normally, calculates a positioning distribution diagram in a high-density pixel grid, and further generates a super-resolution imaging result through superposition, so that the spatial distribution of the cavitation source is subjected to fine resolution of a micron scale.

(7) The invention provides powerful technical means for analyzing the spatial distribution of the cavitation source generated by the excitation of the cavitation excitation ultrasonic transducer, in particular for finely analyzing the spatial distribution of the cavitation source in a small area in a micron scale, can finely monitor and image the cavitation activity in various high-intensity or low-intensity therapeutic ultrasonic applications such as tissue damage, drug delivery and the like, further deeply clarifies the biophysical mechanism behind the applications, and is beneficial to promoting the development of ultrasonic accurate diagnosis and treatment integrated technology.

Claims

1. A super-resolution ultrasonic passive cavitation imaging method based on single cavitation source separation and positioning is characterized in that: the method comprises the following steps:

3) Defining local clustering windows for each global pixel cluster of each frame time gating energy diagram, carrying out frame selection on the time gating energy diagram by utilizing the local clustering windows to obtain a local energy diagram, establishing a local pixel set according to all non-zero value pixels after thresholding the local energy diagram, clustering pixels in the local pixel set by repeating the processes of initializing the pixel cluster and expanding the pixel cluster according to the set closed neighborhood radius and the set closed neighborhood pixel number threshold to obtain a local pixel cluster, searching a local pixel main cluster from the local pixel cluster, and setting the pixel value of all pixels which are not contained in the local pixel main cluster in the local energy diagram to be zero to obtain a single cavitation source energy diagram;

5) Establishing a Gaussian distribution function in an imaging coordinate system, performing Gaussian fitting on each single cavitation source energy diagram in a single cavitation source energy diagram set to obtain an estimated value of a parameter in the Gaussian distribution function, wherein the estimated value of a Gaussian distribution peak value position is an estimated value of a cavitation source coordinate, substituting the estimated value of the parameter into the Gaussian distribution function to calculate a Gaussian fitting result, calculating the estimated uncertainty of the cavitation source coordinate under a given confidence level, and calculating the fitting goodness according to the single cavitation source energy diagram and the Gaussian fitting result of the single cavitation source energy diagram;

2. The super-resolution ultrasonic passive cavitation imaging method based on single cavitation source separation and positioning as claimed in claim 1, wherein the method is characterized in that: in the step 1, the calculation formula of the double-pass delay is expressed as:

wherein,,

N _SE The number of array elements of the array ultrasonic transducer is the sound propagation speed;

3. The super-resolution ultrasonic passive cavitation imaging method based on single cavitation source separation and positioning as claimed in claim 1, wherein the method is characterized in that: in the step 1, the time length of the time gating interval is the pulse length of the cavitation excitation ultrasonic transducer transmitting pulse.

4. The super-resolution ultrasonic passive cavitation imaging method based on single cavitation source separation and positioning as claimed in claim 1, wherein the method is characterized in that: said step 2Or in step 3, the radius of the closed neighborhood is set as the unit of the number of pixels

The number of pixels in the closed neighborhood threshold is set to 2;

in the step 2 or the step 3, the conditions for initializing the pixel cluster are as follows: selecting an undetermined pixel from the global pixel set or the local pixel set, if the number of pixels in the closed neighborhood of the selected pixel is smaller than the threshold value of the number of pixels in the closed neighborhood, reselecting the pixel until the number of pixels in the closed neighborhood of the selected pixel is larger than or equal to the threshold value of the number of pixels in the closed neighborhood, wherein the selected pixel is a core pixel at the moment, and initializing a pixel cluster;

in the step 2 or the step 3, the expanding pixel cluster specifically includes the following steps:

5. The super-resolution ultrasonic passive cavitation imaging method based on single cavitation source separation and positioning as claimed in claim 1, wherein the method is characterized in that: in the step 3, the center of the local clustering window is a pixel corresponding to the mass center of the global pixel cluster, and the size of the local clustering window is determined according to a point spread function; the threshold coefficient used by the thresholding in the step 3 is smaller than that used by the thresholding in the step 2; in the step 3, searching the local pixel main cluster from the local pixel clusters specifically includes the following steps: and setting the label pixel as a pixel corresponding to the mass center of the global pixel cluster, and setting a certain local pixel cluster containing the label pixel as a local pixel main cluster.

6. The super-resolution ultrasonic passive cavitation imaging method based on single cavitation source separation and positioning as claimed in claim 1, wherein the method is characterized in that: in the step 4, the pixel peak value threshold is 0.4-0.6 times of the maximum value of the pixel peaks of all the single cavitation source energy diagrams.

7. The super-resolution ultrasonic passive cavitation imaging method based on single cavitation source separation and positioning as claimed in claim 1, wherein the method is characterized in that: in the step 5, the gaussian distribution function allows angular deflection, and the expression is:

wherein A is Gaussian distribution peak value, mu _x Sum mu _z Gaussian distribution peak position, sigma, in x-axis and z-axis directions, respectively _x Sum sigma _z The standard deviation of the Gaussian distribution in the x-axis direction and the z-axis direction respectively, and the theta is the Gaussian distribution deflection angle;

the Gaussian fitting adopts a nonlinear least square method based on a Levenberg-Marquardt algorithm, and the parameters in the Gaussian distribution function are initialized as follows: initializing a Gaussian distribution peak value to be the maximum pixel value in the single cavitation source energy diagram, initializing the position of the Gaussian distribution peak value to be the coordinate of the maximum pixel value in the single cavitation source energy diagram, and initializing the Gaussian distribution standard deviation according to a point spread function.

8. The super-resolution ultrasonic passive cavitation imaging method based on single cavitation source separation and positioning as claimed in claim 1, wherein the method is characterized in that: in the step 6, the threshold values of uncertainty and goodness of fit of the cavitation source coordinate estimation are respectively:

Wherein,,

and->

Threshold for goodness of fit gof, +.>

And IQR _GOF Respectively arranging elements in the fitting goodness set GOF in an ascending order to obtain lower quartile and quartile spacing of the elements in the set;

in the step 6, detecting the cavitation source positioned normally specifically includes the following steps: if cavitation source coordinates in x-axis direction and z-axis direction estimate uncertainty delta _x And delta _z The goodness of fit gof is simultaneously satisfied

And->

And->

Then the normal positioning is performed.

9. The super-resolution ultrasonic passive cavitation imaging method based on single cavitation source separation and positioning as claimed in claim 1, wherein the method is characterized in that: in the step 7, the number of pixels of the pixel grid of the super-resolution imaging is lambda of the number of pixels of the pixel grid of the time-gated energy map _x ×λ _z Multiple of lambda _x And lambda (lambda) _z The values of the multiples are respectively equal to or more than 10 in the x-axis direction and the z-axis direction; the positioning distribution function is a Gaussian distribution function allowing angular deflection, and the peak position, standard deviation, peak value and deflection angle of the Gaussian distribution function are respectively estimated values of cavitation source coordinates in the x-axis direction and the z-axis direction, estimated uncertainty of cavitation source coordinates in the x-axis direction and the z-axis direction, estimated values of Gaussian distribution peak value and estimated values of Gaussian distribution deflection angle, which are obtained by carrying out Gaussian fitting on a single cavitation source energy diagram.

10. A super-resolution ultrasonic passive cavitation imaging system based on single cavitation source separation and positioning is characterized in that: the system comprises a time gating energy map calculation module, a single cavitation source separation module, a single cavitation source energy image element peak value judgment module, a Gaussian fitting cavitation source positioning module, a cavitation source positioning detection module and a cavitation source positioning distribution map calculation and processing module;

the time gating energy map calculation module is used for establishing an imaging coordinate system and planning a pixel grid of a time gating energy map, calculating double-pass delay according to the distance from each pixel coordinate in the pixel grid to a tangent line or a tangent plane at the center coordinate of the surface of the cavitation excitation ultrasonic transducer and the distance between each pixel coordinate and each array element coordinate of the array ultrasonic transducer, setting the time length of the time gating interval by taking the double-pass delay as the initial time of the time gating interval, extracting the time gating signals of each array element from each frame of original radio frequency signals passively received by the array ultrasonic transducer according to the time gating interval of each array element, superposing the time gating signals of each array element along the array element direction, and superposing the square of the superposed time gating synthesized signals along the direction of a sampling point, thereby obtaining the time gating energy map corresponding to each frame of original radio frequency signals;

The single cavitation source separation module specifically comprises a pixel global clustering sub-module and a pixel local clustering sub-module and a local pixel main cluster searching sub-module;

the pixel global clustering submodule is used for establishing a global pixel set according to all non-zero value pixels after thresholding the time-gating energy map of each frame obtained by the time-gating energy map calculation module, and clustering pixels in the global pixel set by repeating the processes of initializing the pixel clusters and expanding the pixel clusters according to the set closed neighborhood radius and the set closed neighborhood pixel number threshold value, so as to obtain a global pixel cluster of the time-gating energy map of each frame;

the local clustering and local pixel main cluster searching submodule is used for defining a local clustering window for each global pixel cluster of each frame time gating energy image obtained by the pixel global clustering submodule, carrying out frame selection on the time gating energy image by utilizing the local clustering window, carrying out thresholding treatment on the local energy image obtained by the frame selection, establishing a local pixel set according to all non-zero value pixels, carrying out clustering on pixels in the local pixel set according to the set closed neighborhood radius and the set closed neighborhood pixel number threshold value through repeating the processes of initializing the pixel cluster and expanding the pixel cluster, searching the local pixel main cluster from the local pixel cluster obtained by the clustering, and setting the pixel values of all pixels which are not contained in the local pixel main cluster in the local energy image to be zero, thereby obtaining a single cavitation source energy image;

The single cavitation source energy image peak value judging module is used for extracting pixel peak values of all the single cavitation source energy images obtained by the local pixel clustering and local pixel main cluster searching submodule, judging the pixel peak value of each single cavitation source energy image according to a set pixel peak value threshold value, and establishing a single cavitation source energy image set by the single cavitation source energy images with the pixel peak value smaller than the pixel peak value threshold value;

the Gaussian fitting cavitation source positioning module is used for establishing a Gaussian distribution function in an imaging coordinate system, carrying out Gaussian fitting on each single cavitation source energy image in the single cavitation source energy image element peak value judgment module, substituting an estimated value of a parameter in the Gaussian distribution function obtained by fitting into the Gaussian distribution function to calculate a Gaussian fitting result, calculating the uncertainty of cavitation source coordinate estimation under a given confidence level, and calculating the fitting goodness according to the single cavitation source energy image and the Gaussian fitting result of the single cavitation source energy image; in the estimated values of the parameters, the estimated value of the peak position of the Gaussian distribution is the estimated value of the cavitation source coordinates;

the cavitation source positioning detection module is used for establishing a cavitation source coordinate estimation uncertainty set and a fitting goodness set according to the cavitation source coordinate estimation uncertainty and the fitting goodness of all the cavitation sources obtained by the Gaussian fitting cavitation source positioning module, respectively carrying out ascending arrangement on elements in each set, calculating quartiles and quartiles distances of the elements in each set after ascending arrangement, setting a threshold value of the cavitation source coordinate estimation uncertainty and a threshold value of the fitting goodness according to the quartiles and the quartiles distances, and detecting the cavitation source in normal positioning by utilizing the set threshold value of the cavitation source coordinate estimation uncertainty and the threshold value of the fitting goodness;

The cavitation source positioning distribution diagram calculating and processing module is used for establishing a positioning distribution function for each normally positioned cavitation source obtained by the cavitation source positioning detection module, calculating a positioning distribution function value of each pixel coordinate in a planned super-resolution imaging pixel grid, and carrying out normalization and logarithmic processing after overlapping the calculated positioning distribution diagrams of all normally positioned cavitation sources.