An ultrasonic super-resolution imaging method based on adaptive radial basis function interpolation

Abstract

This invention provides an ultrasound super-resolution imaging method based on adaptive radial basis function interpolation, belonging to the field of medical ultrasound imaging and signal processing technology. The invention includes the following steps: Step 1: Acquire ultrasound radio frequency data, obtain low frame rate (IQ) data for clinical ultrasound contrast imaging through orthogonal demodulation, and process it; Step 2: Based on the adaptive radial basis function, interpolate the low frame rate IQ data of clinical ultrasound contrast imaging processed in Step 1 to generate target high frame rate IQ data; Step 3: Perform clutter filtering, microbubble localization and tracking on the interpolated high frame rate IQ data from Step 2 to obtain microbubble trajectories, reconstruct and output a super-resolution vascular structure map. This invention enables conventional low frame rate clinical ultrasound equipment to achieve high-quality micron-level ultrasound super-resolution imaging, possessing extremely high clinical equipment versatility.

Description

Technical Field

[0001] This invention relates to an ultrasound super-resolution imaging method based on adaptive radial basis function interpolation, belonging to the field of medical ultrasound imaging and signal processing technology. Background Technology

[0002] Traditional ultrasound imaging is limited by the fundamental acoustic diffraction limit, making it difficult to resolve the microvascular networks deep within the human body at the micrometer scale. To overcome this physical bottleneck, Ultrasound Localization Microscope (ULM) imaging technology has emerged. ULM injects ultrasound-contrast microbubbles into the bloodstream, using these microbubbles as independent strong point reflectors. It achieves precise centering and trajectory tracking of sparse microbubbles within a continuous high-temporal-resolution image sequence, thus successfully breaking the diffraction limit and enabling micrometer-scale microvascular imaging. However, microbubbles in blood flow are highly dynamic, and the accuracy of ULM is extremely dependent on high frame rates (typically greater than 1000Hz) to avoid spatial aliasing and ensure the temporal continuity of microbubble localization and tracking. In practical applications, however, due to limitations in system hardware transmission bandwidth or traditional line-by-line scanning imaging modes, conventional clinical ultrasound equipment typically operates at low frame rates below 200Hz. At such low sampling rates, microbubbles can undergo excessive displacement between adjacent frames, leading to serious problems such as low vascular structure reconstruction, microbubble tracking failure, and distortion of microhemodynamic estimation.

[0003] Existing techniques often employ linear or spline interpolation to improve frame rates, but these methods struggle to capture the nonlinear motion characteristics of microbubbles in complex blood flow environments. Recent studies have proposed using radial basis functions (RBF) for two-dimensional spatiotemporal interpolation to recover missing frames, demonstrating RBF's advantages in fitting microbubble trajectories. However, traditional RBF interpolation methods typically use globally fixed shape parameters. When processing non-uniformly distributed microbubble signals, dense microbubble regions are prone to over-smoothing, while sparse regions are susceptible to artifacts. Furthermore, in low-frame-rate clinical data, fixed interpolation matrices often face matrix singularity and inversion instability due to excessively large condition numbers. Therefore, a high-precision interpolation scheme capable of adaptively adjusting parameters is urgently needed to balance clinical equipment compatibility and imaging accuracy. Summary of the Invention

[0004] The purpose of this invention is to solve the problems of low frame rate ULM imaging accuracy and difficulty in balancing nonlinear motion fitting and numerical stability in the existing technology, and to provide an ultrasound super-resolution imaging method based on adaptive radial basis function interpolation.

[0005] The objective of this invention is achieved through the following technical solution:

[0006] An ultrasound super-resolution imaging method based on adaptive radial basis function interpolation includes the following steps:

[0007] Step 1: Acquire ultrasound radiofrequency data, obtain low frame rate IQ data for clinical ultrasound contrast imaging through orthogonal demodulation, and process the data.

[0008] Step 2: Based on the adaptive radial basis function, interpolate the low frame rate IQ data of clinical ultrasound contrast imaging after Step 1 to generate the target high frame rate IQ data;

[0009] Step 3: Perform clutter filtering, microbubble localization and tracking on the high frame rate IQ data after interpolation in Step 2 to obtain microbubble trajectories, reconstruct and output super-resolution vascular structure maps.

[0010] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0011] 1. This invention reduces the dimensionality of complete data into multiple independent... By processing and recombining the planar data, the spatial and temporal information of ultrasound is fully utilized to improve the interpolation accuracy of microbubble motion, while avoiding the high computational cost brought about by global three-dimensional spatiotemporal matrix operations.

[0012] 2. This invention introduces adaptive shape parameters. When there are dense and high-speed moving microbubbles in the local target area, the corresponding quadratic function is a small-scale sharpening kernel to preserve high-frequency motion details; when the microbubbles in the local area are extremely sparse and the motion is gentle, the corresponding quadratic function is a large-scale smoothing kernel to ensure trajectory continuity and improve trajectory fitting accuracy.

[0013] 3. This invention solves the singularity problem that traditional RBF interpolation matrices are prone to when processing low frame rate and high noise clinical data by introducing a dynamic regularization term that adapts to the time interval, thereby enhancing the stability of numerical calculation.

[0014] 4. This invention, by designing a multi-scale kernel function fusion algorithm, breaks through the limitation of a single kernel function, and can simultaneously take into account the macroscopic trend and microscopic pulsation of microbubble motion, thereby improving the ability to capture multi-scale features.

[0015] 5. Experimental results show that the method of the present invention can upscale low frame rate data of about 100Hz to more than 800Hz with high fidelity, enabling conventional low frame rate clinical ultrasound equipment to achieve high-quality micron-level ultrasound super-resolution imaging, and has extremely high clinical equipment versatility. Attached Figure Description

[0016] Figure 1 This is a flowchart of an ultrasound super-resolution imaging method based on adaptive radial basis function interpolation according to the present invention.

[0017] Figure 2 for A schematic diagram of local feature extraction using a sliding window for two-dimensional spatiotemporal plane data.

[0018] Figure 3 The image shows a comparison of the super-resolution vessel reconstruction results before and after interpolation using the method of this invention; wherein:

[0019] Figure 3 (a) is a reconstruction image of the original 800Hz high frame rate data using an ultrasound positioning microscope;

[0020] Figure 3 (b) is a reconstruction image from a 100Hz low frame rate ultrasound positioning microscope;

[0021] Figure 3 (c) is a reconstruction image of ultrasonic positioning microscope based on 100Hz downsampled data and interpolated back to 800Hz using the method of this invention. Detailed Implementation

[0022] The present invention will be further described in detail below with reference to the accompanying drawings: This embodiment is implemented under the premise of the technical solution of the present invention, and detailed implementation methods are given, but the protection scope of the present invention is not limited to the following embodiments.

[0023] like Figure 1 As shown in the figure, the ultrasound super-resolution imaging method based on adaptive radial basis function interpolation involved in this embodiment includes the following steps:

[0024] Step 1: Processing of low frame rate IQ data from clinical ultrasound contrast imaging:

[0025] By emitting ultrasonic waves and receiving echo signals through the probe of an ultrasound imaging system, continuous ultrasound radio frequency data of the target area after injection of ultrasound contrast microbubbles is acquired. The received radio frequency data is then orthogonally demodulated to obtain a low-frame-rate IQ data matrix. .

[0026] in, This represents the horizontal coordinate in space, i.e., the orientation of the probe array elements. Represents the spatial axial coordinates. This represents time-series coordinates. Unlike traditional envelope or ultrasound image data, this IQ data preserves complete phase information of the echo signal. Due to the hardware transmission bandwidth limitations of conventional low-frame-rate clinical equipment, the sampling frame rate of this data is typically far lower than the ideal high-frequency acquisition standard for microbubble localization tracking.

[0027] Considering that the interpolation of the three-dimensional spatiotemporal motion of microbubbles in complex blood flow requires extremely high computational resources, in order to avoid the high computational cost and memory overhead caused by global three-dimensional spatiotemporal matrix operations, this invention addresses the IQ data matrix. Along the horizontal axis of space (i.e. Perform line-by-line slicing on the axis. Specifically, fix the horizontal coordinate. Extracting data from all depth directions and time dimensions at that location, thus decomposing the complete 3D data block laterally into... A completely independent Two-dimensional planar data ,in This represents the total number of horizontal scan lines, and This dimensionality reduction decomposition step transforms the complex 3D interpolation problem into multiple independent 2D interpolation problems by processing only a single 2D plane slice for each independent algorithmic step, significantly improving the system's computational efficiency.

[0028] Based on this, for each independent element after dimensionality reduction decomposition Two-dimensional planar data Using data sampling points as the basic unit, a spatiotemporal vector containing spatial depth and temporal features is constructed. Specifically, the spatial depth position of the data points within the slice plane. With the time frame Extract and combine to construct a two-dimensional spatiotemporal coordinate vector. The essence of this operation is to unify spatial distance parameters and temporal span parameters, which originally had completely different physical dimensions, into the same mathematical metric space, thereby forming the basic input parameters for subsequent dynamic calculation of adaptive multiple quadratic function (MQ) shape parameters and interpolation node distances. This unified spatiotemporal vectorization representation ensures that the system can simultaneously evaluate the spatial displacement changes of microbubbles and their evolution characteristics over time.

[0029] Step 2: Interpolate low frame rate IQ data from clinical ultrasound contrast imaging based on adaptive radial basis functions:

[0030] Step 2.1: Extract local spatiotemporal signal features and dynamically calculate the adaptive shape parameters of the corresponding multiple quadratic functions (MQ). :

[0031] For each independent element obtained from the dimensionality reduction decomposition in step 1 Two-dimensional planar data ,like Figure 2 As shown, the present invention sets a size of Local spatiotemporal sliding window (in Window length representing spatial depth. (Representing the window length in the time dimension). Using this sliding window... Perform a non-overlapping sliding traversal on the plane, targeting the center point within the window. Local spatial energy density and temporal envelope gradient are extracted as local spatiotemporal features, respectively.

[0032] First, calculate the local space energy density. This characteristic parameter is used to quantify the spatial density of microbubbles and the intensity of scattered echoes within the local observation area. Its calculation formula is as follows:

[0033]

[0034] In the formula, For the sliding window Neighborhood coordinates within; The signal envelope amplitude of the IQ data at this neighborhood point; A weighting function for a predefined two-dimensional space. Weighting function A two-dimensional Gaussian kernel function is used to assign higher energy weights to pixels closer to the central target point, thereby highlighting the signal contribution of the local core region and reducing the interference of edge noise.

[0035] The time envelope gradient is then calculated. This characteristic parameter is used to quantify the rate of change of microbubble echo signals over time, thereby objectively reflecting the intensity of microbubble blood flow in the local area. Its calculation formula is as follows:

[0036]

[0037] Through the aforementioned sliding window mechanism and mathematical feature modeling, this invention can accurately and in real-time capture the differences in physical properties of microbubbles at different anatomical locations, such as the center of large blood vessels with dense microbubbles and high flow rates, or the ends of capillaries with sparse microbubbles and slow flow rates. When the local spatial energy density... and temporal envelope gradient When both values ​​are high, it indicates that there is a high-density and high-speed nonlinear movement of microbubble groups in the current window region; conversely, it indicates that the microbubbles in the region are sparsely distributed or in a slow laminar flow state.

[0038] In obtaining the local spatial energy density and temporal envelope gradient Subsequently, this invention constructs a nonlinear negative exponential mapping model based on local features, used to dynamically calculate the adaptive shape parameters of the corresponding multiple quadratic functions (MQ). The specific calculation formula is as follows:

[0039]

[0040] In the formula, These are the preset upper and lower limits of the system physical parameters for the MQ shape parameters, and their specific values ​​are calibrated based on the center frequency and imaging depth of the ultrasound probe. Each is currently being processed The global maximum spatial energy density and maximum temporal envelope gradient within two-dimensional planar data are used as normalization factors to uniformly map input features to a dimensionless model. Within the range; This is an empirical adjustment coefficient, with values ​​ranging from [value range to value]. This is used to control the sensitivity weights of spatial and temporal features to the final shape parameters, respectively.

[0041] When there are dense and high-speed moving microbubbles in a local target area (i.e., local spatial energy density) and temporal envelope gradient When all values ​​tend to be high, the negative exponent term in the formula rapidly decays and approaches 0, resulting in the calculated shape parameters. Approaching the lower realm In this case, the corresponding quadratic function manifests as a small-scale sharpening kernel, primarily relying on effective data points in the nearest neighborhood for interpolation, thereby accurately preserving high-frequency local motion details. Conversely, when microbubbles in a local region are extremely sparse and their motion is gentle (i.e., local spatial energy density...), the... and temporal envelope gradient When all values ​​tend to be low, the negative exponent term approaches 1, and the shape parameter... Automatically approaching the upper limit In this case, the corresponding quadratic function manifests as a large-scale smoothing kernel, which expands the scope of spatial features to cross data blind spots, thus ensuring the global continuity of sparse microbubble trajectories. This dynamic mechanism can accurately capture the subtle and complex changes in blood flow within microvessels, and effectively overcomes the problems of oversmoothing or trajectory discontinuity that are easily caused by traditional fixed-parameter interpolation, significantly improving the spatiotemporal fitting accuracy of microbubble motion trajectories.

[0042] Among them, the small-scale sharpening kernel (corresponding to shape parameters approaching the lower bound): In radial basis functions, when the shape parameter is extremely small, the function image appears very "sharp" (strong local effect). This corresponds to the main vascular region with dense microbubbles and fast flow velocity. In this case, interpolation only strongly depends on the sampling points at very close distances, thus avoiding smoothing out the surrounding irrelevant microbubble trajectories, with the aim of accurately preserving motion details such as high-frequency turbulence or pulsation. The large-scale smoothing kernel (corresponding to shape parameters approaching the upper bound): When the shape parameter is large, the function image appears relatively "smooth and broad" (strong global effect). This corresponds to the capillary region with sparse microbubbles and slow flow velocity. Due to the small number of effective sampling points, it is necessary to expand the range of action of the basis function to cross the spatiotemporal blind zone, with the aim of filling the gaps, ensuring the global continuity of microbubble motion trajectories, and preventing discontinuity.

[0043] Step 2.2: Introduce a regularization operator that dynamically adjusts with the sampling time interval. Based on the adaptive shape parameters in Step 2.1, construct the MQ interpolation matrix with dynamic regularization term and solve the equation:

[0044] In practical clinical applications, ultrasound equipment is limited by low frame rates below 200Hz, causing large-scale displacement of microbubbles between adjacent frames. Traditional radial basis interpolation typically uses globally fixed empirical regularization parameters. When dealing with highly dynamic nonlinear motion of microbubbles, fixed parameters face severe challenges: parameters that are too small cannot suppress interpolation oscillations and matrix singularities over long time intervals; parameters that are too large will lead to excessive smoothing of high-frequency motion details within large blood vessels. To address this issue, this invention overcomes the limitations of traditional static global regularization by introducing a spatiotemporal dynamic regularization operator that is deeply coupled with both the hardware sampling time interval and the local hemodynamic characteristics of microbubbles (spatial energy density and temporal envelope gradient).

[0045] First, a basic multiple quadratic function (MQ) interpolation matrix is ​​constructed. Based on the adaptive shape parameters extracted in the previous steps, for the... The basic interpolation matrix of a local observation set consisting of discrete sampling points. elements The calculation formula is as follows:

[0046]

[0047] In the formula, and These represent the two-dimensional spatiotemporal coordinates of any two valid echo signal data points within a local spatiotemporal sliding window; The adaptive shape parameters corresponding to the center coordinates of the window, calculated in the previous step, are deeply coupled by the basic matrix, which incorporates spatial distance, temporal span, and local blood flow physical characteristics.

[0048] Based on this, a dynamic regularization coefficient based on spatiotemporal uncertainty and physical priors is constructed. The increased uncertainty in the actual motion trajectory of microbubbles necessitates enhanced regularization penalties to ensure smoothness; while the local spatial energy density... and temporal envelope gradient At higher levels (characterized by a region of dense microbubble density and intense blood flow in the core of large blood vessels), the regularization penalty needs to be reduced to preserve high-frequency hemodynamic details. Combining these opposing physical constraints, this invention designs a dynamic regularization coefficient. The calculation formula is as follows:

[0049]

[0050] In the formula, It is the basic environmental noise regularization parameter, which is calibrated according to the thermal noise level of the ultrasonic system substrate and is usually set to 1e-6. This is the local sampling time interval involved in the current interpolation calculation; This is the preset reference high sampling time interval threshold for this method; This is the time interval penalty control coefficient; These are the global maximum spatial energy density and the maximum temporal envelope gradient, respectively. These are the regularization sensitivity coefficients for spatial and temporal characteristics, respectively.

[0051] Based on the above formula for calculating the dynamic regularization coefficient, when processing capillary regions with sparse microbubbles and large frame intervals, the feature term decays, and the time interval penalty term dominates. The feature term significantly increases, thereby strongly suppressing noise and interpolation oscillations; when processing the main vascular region with dense microbubbles and high flow velocity, the feature term partially offsets the time interval penalty. Adaptive reduction to prevent over-smoothing.

[0052] Finally, the dynamic regularization coefficients and penalty operators are introduced into the basic interpolation matrix to construct the final MQ interpolation system equations with dynamic regularization terms:

[0053]

[0054] In the formula, For the matrix An identity matrix of the same dimension is used as a standard regularization operator; The input is a column vector of low frame rate IQ observation data; Let be the high-precision interpolation coefficient vector to be solved. By introducing a dynamic regularization term that adaptively adjusts with time intervals, this invention solves the singularity problem that traditional RBF interpolation matrices are prone to when processing low frame rate, high-noise clinical data at a mathematical level. This mechanism can find the optimal solution between protecting realistic high-frequency motion characteristics and suppressing numerical calculation errors, thus enhancing the stability of numerical calculations.

[0055] Step 2.3: Solve for the interpolation coefficients and use a multi-scale kernel function fusion algorithm to interpolate the IQ data to generate the target high frame rate IQ data.

[0056] Based on the multi-quadratic function interpolation system equations with dynamic regularization terms constructed in step 2.2, the column vector of low frame rate IQ observation data within the local observation window is... Let these be known quantities. Due to the introduction of the dynamic regularization operator, the coefficient matrix of this system already possesses good positive definiteness and well-state characteristics. In practical engineering calculations, robust numerical algebraic methods such as Singular Value Decomposition (SVD) or Cholesky Decomposition are preferred for solving this system of linear equations, thereby obtaining the interpolation coefficient vector within the current local spatiotemporal window. coefficient vector It contains the spatiotemporal distribution weights of the microbubble scattering signal in this local region.

[0057] Building upon this foundation, a multi-scale kernel function fusion algorithm is constructed. Existing interpolation methods based on a single kernel function struggle to simultaneously capture both the macroscopic laminar flow trend and the microscopic local turbulent fluctuations of microbubble motion. In real microvascular networks, blood flow velocity and direction exhibit significant spatial heterogeneity. To overcome the limitations of a single kernel function, this invention designs a multi-scale kernel function fusion algorithm incorporating different scale scaling factors, constructing a multi-scale spatial fusion basis function. The calculation formula is as follows:

[0058]

[0059] In the formula, The number of scale layers set for the system is taken as follows: These respectively characterize the macroscopic main blood flow, mesoscopic branch laminar flow, and microscopic terminal turbulent flow characteristics of microvessels; The spatiotemporal Euclidean distance between the target interpolation point and the known sampling point; Scaling factors at different scale levels (0.5, 1.0, and 2.0 respectively) are used to adapt the shape parameters. Perform non-linear broadening or sharpening in different dimensions; For the normalized fusion weights corresponding to the kernel functions at each scale, based on The size must be set proportionally and must meet the constraints. By using the weighted fusion of this multi-scale kernel function, the algorithm can accurately capture the minute high-frequency pulsation features of microbubbles while preserving the continuity of their macroscopic motion trajectories, thus improving the ability to capture multi-scale features.

[0060] Finally, based on the set target high frame rate, target time grid points are inserted between the existing time sampling sequences. Using the interpolation coefficient vector obtained from the solution... With the above multi-scale fusion kernel function Interpolation calculations are performed on the IQ data at the new time grid points. Spatiotemporal reconstruction is performed at this location:

[0061]

[0062] In the formula, This represents the total number of valid sampling points within the local window. The spatiotemporal coordinates of the original known sampling points; This refers to the reconstructed local target high frame rate IQ data. For all independent... Repeat steps 2.1 to 2.3 for the two-dimensional planar local window data, recombining them into a complete form. Two-dimensional planar data. Finalized according to the original horizontal spatial coordinates. For complete multiple copies By reconstructing two-dimensional planar data, the complete data can be restored into a three-dimensional high frame rate data volume without loss. This step achieves spatiotemporal oversampling at the underlying complex signal level, completely solving the technical pain points of microbubble spatial aliasing and trajectory disconnection at low frame rates.

[0063] Step 3: Perform clutter filtering, microbubble localization and tracking on the interpolated high frame rate IQ data to obtain microbubble trajectories, reconstruct and output a high-quality super-resolution vascular structure map:

[0064] After acquiring the target high frame rate IQ data, this invention achieves super-resolution imaging through B-mode data conversion, clutter filtering, localization tracking, and trajectory reconstruction. First, the reconstructed high frame rate IQ data undergoes envelope extraction and logarithmic compression to convert it into a high temporal resolution B-mode image sequence. Since blood echoes contain strong tissue clutter, this invention employs a dual clutter filtering technique based on singular value decomposition (SVD) and temporal filtering. By constructing a spatiotemporal matrix from the B-mode image sequence and performing SVD, an eigenvalue truncation threshold is set to filter out low-order principal components corresponding to slowly moving tissues. Subsequently, a Butterworth temporal filter is introduced in the temporal dimension for secondary purification, filtering out residual micro-moving tissue artifacts and extracting a pure microbubble signal sequence.

[0065] For clean microbubble image sequences, this invention employs a radially symmetric (RS) algorithm for center extraction and sub-pixel-level localization of individual microbubbles. This algorithm leverages the approximate Gaussian point spread function of microbubbles, calculating the radially symmetric center of the image grayscale gradient to achieve high-precision localization exceeding the physical diffraction limit, obtaining the spatial coordinate set of microbubbles within each frame. Subsequently, this invention utilizes the Hungarian algorithm for globally optimal data association and coordinate matching between adjacent frames. To address the flickering or cross-occlusion of microbubbles in complex blood flow, a Kalman filter is embedded in the tracking framework, using its state equation for trajectory prediction and error updating, compensating for missed disconnected microbubbles, and ensuring the spatiotemporal continuity of the microbubble motion trajectory. The specific steps are as follows:

[0066] For clean microbubble image sequences, this invention employs a radially symmetric (RS) algorithm to extract the center of individual microbubbles and locate them at the sub-pixel level, obtaining the spatial coordinate set of microbubbles within each frame. After obtaining the discrete microbubble coordinates, this invention constructs a multi-target tracking framework based on a "prediction-association-update" closed loop to concatenate microbubble trajectories between adjacent frames. This tracking framework specifically includes the following three collaborative steps, where Kalman filtering is embedded as the core state estimator in the prediction and update stages:

[0067] One is prior state prediction: in processing the first When viewing a frame of an image, the prediction equation using Kalman filtering is based on the first frame. Given the established trajectory and velocity of the microbubble in the frame, calculate and output the current position of the microbubble. The prior predicted coordinates and predicted covariance matrix of the frame. This step predefines the possible range of microbubbles.

[0068] Secondly, there is the globally optimal data association: the "prior predicted coordinates" output by the Kalman filter are compared with the radial symmetric algorithm at the... The actual observed coordinates extracted from the frame are used to calculate the spatial distance and construct a cost matrix. Then, the Hungarian algorithm is used to perform bipartite graph matching on the cost matrix to obtain the globally optimal correspondence between the predicted trajectory and the actual observed points.

[0069] Thirdly, there is the posterior state and error update: Based on the matching results of the Hungarian algorithm, if a trajectory successfully matches a real observation point, the real observation coordinates are input into the Kalman filter to calculate the Kalman gain, and the prior predicted coordinates are weighted and corrected to complete the posterior update of the microbubble position for that frame, and then the loop for the next frame begins; if a trajectory fails to match a real observation point due to microbubble flickering or occlusion, the "prior predicted coordinates" of the Kalman filter are directly used as the compensation coordinates for that frame to fill in the trajectory, thereby effectively compensating for missed or disconnected microbubbles and ensuring the spatiotemporal continuity and smoothness of the overall microbubble motion trajectory.

[0070] After completing microbubble localization and trajectory tracking for all time frames, all extracted effective microbubble motion trajectories are mapped onto a preset ultra-high resolution spatial grid. This invention effectively compensates for the information loss in the temporal dimension of low-frame-rate clinical data through the aforementioned adaptive radial basis function interpolation process, eliminating microbubble spatial aliasing and trajectory discontinuity caused by insufficient sampling. Finally, by accumulating and superimposing these highly continuous microbubble trajectory points and performing spatial density rendering, a high-contrast micron-level super-resolution vascular structure map is reconstructed and output, achieving high-fidelity representation of the fine vascular morphology of the target region on low-frame-rate data.

[0071] Step 4, Experimental Verification and Application:

[0072] To verify the effectiveness and superiority of the super-resolution imaging method described in this invention, experimental verification was conducted using mouse brain microvascular ultrasound contrast imaging data. The experiment first acquired high-frame-rate data with a sampling rate of 800Hz as a real reference standard. To simulate the low-frame-rate acquisition environment of conventional clinical equipment, this high-frame-rate data was downsampled 8 times in the time dimension to construct a low-frame-rate test dataset with an equivalent frame rate of 100Hz. Specifically, RMSE (Root Mean Square Error) and SSIM (Structural Similarity Score) evaluation metrics were used to quantitatively assess image quality. In the verification of the ultrasound super-resolution imaging method, the motion of in vivo microbubbles is highly transient. If two independent physical acquisitions are performed in the in vivo experiment (e.g., first acquiring a segment of 800Hz data, then switching equipment parameters to acquire a segment of native 100Hz data), due to real-time changes in hemodynamics, the spatial distribution, motion trajectory, and instantaneous flow velocity of the microbubbles captured in these two acquisitions will have fundamental physical differences. This physiological bias caused by 'non-homogeneous data' completely masks the errors in the interpolation algorithm itself, rendering pixel-by-pixel quantitative comparison (RMSE / SSIM) meaningless. Therefore, based on the standard quantitative verification paradigm in the field of super-resolution image processing, a downsampling factor is applied to high frame rate data to simulate low frame rate data acquisition. By rigorously and uniformly extracting the same batch of 800Hz high frame rate raw data in the time dimension (i.e., downsampling by 8 times to construct an equivalent 100Hz test set), the random interference of microbubble distribution caused by different acquisitions is eliminated, preserving the original high frame rate data as an absolute reference standard. Thus, the ultra-high frequency restoration capability of the method of this invention for clinical low frame rate data is objectively and purely demonstrated.

[0073] First, the ultrasound super-resolution imaging method based on adaptive radial basis function interpolation proposed in this invention was applied to the 100Hz dataset, and then compared with low frame rate test data for microvascular reconstruction. Figure 3As shown, the low frame rate data reconstruction results exhibit significant microvascular trajectory breaks and blur artifacts, with RMSE and SSIM values ​​of 24.02 and 0.6607, respectively. Conversely, the high frame rate reconstruction results obtained based on the method of this invention successfully and faithfully recover the missing vascular information, reducing the RMSE to 11.08 and improving the SSIM to 0.8724.

[0074] Experimental results show that, in terms of data, the method of this invention can interpolate low frame rate data of around 100Hz to an equivalent high frame rate of over 800Hz (the interpolation upsampling factor can be freely adjusted as needed, generally not exceeding 10 times to balance effect and computational efficiency); simultaneously, in terms of image quality, the reconstructed super-resolution vascular structure map highly approximates the real reference image generated from the original 800Hz data in terms of morphological smoothness, global trajectory continuity, and preservation of minute branch details. This experiment fully demonstrates that the method described in this invention can effectively overcome the technical bottleneck caused by low sampling rates, enabling conventional 100Hz-level low frame rate clinical ultrasound equipment to achieve micron-level hemodynamic estimation and fine structural imaging comparable to high-frequency equipment, significantly improving the clinical applicability and translational value of this technology.

[0075] The above description is merely a preferred embodiment of the present invention. These specific embodiments are different implementations based on the overall concept of the present invention, and the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. An ultrasound super-resolution imaging method based on adaptive radial basis function interpolation, characterized in that, Includes the following steps: Step 1: Acquire ultrasound radiofrequency data, obtain low frame rate IQ data for clinical ultrasound contrast imaging through orthogonal demodulation, and process the data. Step 2: Based on the adaptive radial basis function, interpolate the low frame rate IQ data of clinical ultrasound contrast imaging after Step 1 to generate the target high frame rate IQ data; Step 3: Perform clutter filtering, microbubble localization and tracking on the high frame rate IQ data after interpolation in Step 2 to obtain microbubble trajectories, reconstruct and output super-resolution vascular structure maps.

2. The ultrasound super-resolution imaging method based on adaptive radial basis function interpolation according to claim 1, characterized in that, The specific steps of step one are as follows: Step 11: The ultrasound imaging system emits ultrasound waves and receives echo signals to acquire continuous ultrasound radio frequency data of the target area after the injection of ultrasound contrast microbubbles. The received ultrasound radio frequency data is then orthogonally demodulated to obtain a low frame rate IQ data matrix. ,in, This represents the horizontal coordinate in space, i.e., the orientation of the probe array elements. Represents the spatial axial coordinates. Represents time series coordinates; Steps 1 and 2: The IQ data matrix obtained in Step 1. Along the horizontal axis of space Perform line-by-line slicing on the axis, fixing the horizontal coordinate. Extracting data from all depth directions and time dimensions at that location, thus decomposing the complete 3D data block laterally into... A completely independent Two-dimensional planar data ,in This represents the total number of horizontal scan lines, and .

3. The ultrasound super-resolution imaging method based on adaptive radial basis function interpolation according to claim 1, characterized in that, The specific steps of step two are as follows: Step 21: Extract local spatiotemporal signal features and dynamically calculate the adaptive shape parameters of the corresponding quadratic functions. ; Step 22: Introduce a regularization operator that dynamically adjusts with the sampling time interval. Based on the adaptive shape parameters in Step 21, construct a multi-quadratic function interpolation matrix with dynamic regularization term to solve the equation. Step 23: Solve for the interpolation coefficients in Step 22, and use the multi-scale kernel function fusion algorithm to interpolate the IQ data to generate the target high frame rate IQ data.

4. The ultrasound super-resolution imaging method based on adaptive radial basis function interpolation according to claim 3, characterized in that, The specific steps of step two-one are as follows: For each independent element obtained from the dimensionality reduction decomposition Two-dimensional planar data Set a size as Local spatiotemporal sliding window ,in Window length representing spatial depth. The window length represents the time dimension; this sliding window is used to... Perform a non-overlapping sliding traversal on the plane, targeting the center point within the window. Local spatial energy density and temporal envelope gradient are extracted as local spatiotemporal features, respectively. First, calculate the local space energy density. This characteristic parameter is used to quantify the spatial density of microbubbles and the intensity of scattered echoes within the local observation area. Its calculation formula is as follows: In the formula, For the sliding window Neighborhood coordinates within; The signal envelope amplitude of the IQ data at this neighborhood point; The weighting function is a preset two-dimensional space weighting function. A two-dimensional Gaussian kernel function is used; The time envelope gradient is then calculated. This characteristic parameter is used to quantify the rate of change of microbubble echo signals over time, thereby objectively reflecting the intensity of microbubble blood flow in the local area. Its calculation formula is as follows: When the local space energy density and temporal envelope gradient When all values ​​are high, it indicates that there is a high-density and high-speed nonlinear movement of microbubble groups in the current window region; conversely, it indicates that the microbubbles in the region are sparsely distributed or in a slow laminar flow state. In obtaining the local spatial energy density and temporal envelope gradient Subsequently, a nonlinear negative exponential mapping model based on local features is constructed to dynamically calculate the adaptive shape parameters of the corresponding multiple quadratic functions. The specific calculation formula is as follows: In the formula, These are the upper and lower limits of the system physical parameters for the preset multi-quadratic function shape parameters, and their specific values ​​are calibrated based on the center frequency and imaging depth of the ultrasonic probe; Each is currently being processed The global maximum spatial energy density and maximum temporal envelope gradient within two-dimensional planar data are used as normalization factors to uniformly map input features to a dimensionless model. Within the range; These are empirical adjustment coefficients used to control the sensitivity weights of spatial and temporal features to the final shape parameters, respectively. When dense, high-speed microbubbles exist in the local target region, meaning both the local spatial energy density and temporal envelope gradient tend to be high, the negative exponential term in the formula rapidly decays and approaches 0, and the calculated shape parameter approaches the lower bound. At this point, the corresponding quadratic function manifests as a small-scale sharpening kernel, relying on effective data points in the nearest neighborhood for interpolation, thereby accurately preserving high-frequency local motion details. Conversely, when microbubbles in local regions are extremely sparse and their motion is gentle, i.e., when both are low, the negative exponent term approaches 1, and the shape parameter automatically approaches the upper bound. At this point, the corresponding quadratic function manifests as a large-scale smoothing kernel, which expands the scope of spatial features to cross data blind spots, thereby ensuring the global continuity of sparse microbubble trajectories.

5. The ultrasound super-resolution imaging method based on adaptive radial basis function interpolation according to claim 4, characterized in that, The specific steps of step two are as follows: First, construct the basic multiquadratic function interpolation matrix. Based on the adaptive shape parameters extracted in step two, for the... The basic interpolation matrix of a local observation set consisting of discrete sampling points. elements The calculation formula is as follows: In the formula, and These represent the two-dimensional spatiotemporal coordinates of any two valid echo signal data points within a local spatiotemporal sliding window; The adaptive shape parameters corresponding to the center coordinates of the window, calculated in the previous step; Based on this, a dynamic regularization coefficient based on spatiotemporal uncertainty and physical priors is constructed. Dynamic regularization coefficient The calculation formula is as follows: In the formula, These are the basic environmental noise regularization parameters. This is the local sampling time interval involved in the current interpolation calculation; This is the preset reference high sampling time interval threshold for this method; This is the time interval penalty control coefficient; These are the global maximum spatial energy density and the maximum temporal envelope gradient, respectively. These are the regularization sensitivity coefficients for spatial and temporal characteristics, respectively; Based on the above formula for calculating the dynamic regularization coefficient, when processing capillary regions with sparse microbubbles and large frame intervals, the feature term decays, and the time interval penalty term dominates. The feature term significantly increases, thereby strongly suppressing noise and interpolation oscillations; when processing the main vascular region with dense microbubbles and high flow velocity, the feature term partially offsets the time interval penalty. Adaptive reduction to prevent over-smoothing; Finally, the dynamic regularization coefficients are introduced into the basic interpolation matrix to construct the final equations for the multi-quadratic function interpolation system with dynamic regularization terms: In the formula, To be interpolated with the basic interpolation matrix An identity matrix of the same dimension is used as a standard regularization operator; The input is a column vector of low frame rate IQ observation data; The vector of high-precision interpolation coefficients to be solved.

6. The ultrasound super-resolution imaging method based on adaptive radial basis function interpolation according to claim 5, characterized in that, The specific steps of steps two and three are as follows: Based on the multi-quadratic function interpolation system equations with dynamic regularization constructed in step 22, the low frame rate IQ observation data column vector within the local observation window will be used. Given these quantities, we solve the system of linear equations using numerical algebra methods to obtain the interpolation coefficient vector within the current local spatiotemporal window. coefficient vector It contains the spatiotemporal distribution weights of the microbubble scattering signal within this local region; Based on this, a multi-scale kernel function fusion algorithm incorporating different scaling factors is designed, and a multi-scale spatial fusion basis function is constructed. The calculation formula is as follows: In the formula, The number of scale layers set for the system is taken as follows: These respectively characterize the macroscopic main blood flow, mesoscopic branch laminar flow, and microscopic terminal turbulent flow characteristics of microvessels; The spatiotemporal Euclidean distance between the target interpolation point and the known sampling point; Scaling factors at different scale levels are used to adapt shape parameters. Perform non-linear broadening or sharpening in different dimensions; For the normalized fusion weights corresponding to the kernel functions at each scale, based on The size must be set proportionally and must meet the constraints. ; Finally, based on the set target high frame rate, target time grid points are inserted between the existing time sampling sequences. Using the interpolation coefficient vector obtained by solving With the above multi-scale fusion kernel function Interpolation calculations are performed on the IQ data at the new time grid points. Spatiotemporal reconstruction is performed: High frame rate IQ data of local targets generated by the reconstruction. The calculation formula is as follows: In the formula, This represents the total number of valid sampling points within the local window. For the original known sampling points, the spatiotemporal coordinates are given; for all independent Repeat steps two through three on the two-dimensional local window data, recombine them into a complete array. Two-dimensional planar data, ultimately based on the original horizontal spatial coordinates For complete multiple copies By reconstructing two-dimensional planar data, the complete data can be restored into a three-dimensional high frame rate data volume without loss.

7. The ultrasound super-resolution imaging method based on adaptive radial basis function interpolation according to claim 1, characterized in that, The specific steps of step three are as follows: Step 31: For the interpolated high frame rate IQ data, perform envelope extraction and logarithmic compression to convert it into a high temporal resolution B-mode image sequence; adopt a dual clutter filtering technique based on singular value decomposition and temporal filtering, construct a spatiotemporal matrix for the B-mode image sequence and perform singular value decomposition, set an eigenvalue truncation threshold, and filter out the low-order principal components corresponding to slow-moving tissues. Subsequently, a Butterworth time-domain filter was introduced in the time dimension for secondary purification to filter out residual micro-motion artifacts and extract a pure microbubble signal sequence. Step 32: For the pure microbubble signal sequence extracted in Step 31, the radial symmetric algorithm is used to extract the center of each microbubble and locate it at the sub-pixel level, so as to obtain the spatial coordinate set of microbubbles in each frame; Subsequently, the Hungarian algorithm was used to perform globally optimal data association and coordinate matching between adjacent frames. To address the flickering or cross-occlusion of microbubbles in complex blood flow, a Kalman filter was embedded in the tracking framework. Its state equation was used for trajectory prediction and error update to compensate for missed disconnected microbubbles and ensure the spatiotemporal continuity of the microbubble motion trajectory. Step 33: After completing the microbubble localization and trajectory tracking for all time frames, all the extracted effective microbubble motion trajectories are mapped onto a preset ultra-high resolution spatial grid. By accumulating and superimposing these highly continuous microbubble trajectory points and performing spatial density rendering, a high-contrast micron-level super-resolution vascular structure map is reconstructed and output.

Images