Phase distortion correction method for transcranial ultrasonic imaging

Abstract

The invention belongs to the technical field of ultrasonic imaging, and particularly relates to a phase distortion correction method for craniocerebral ultrasonic imaging. The method comprises the following steps: acquiring ultrasonic echo data in a skull scene and a B-mode image obtained after beam forming; identifying the skull boundary of the transcranial B-mode image, and constructing a craniocerebral sound velocity model; solving a sound wave travel time matrix under the current sound velocity model on the basis of a ray theory method or an EiKonal equation, and applying the sound wave travel time matrix to beam forming to obtain an image after phase distortion correction; and a transcranial imaging result after beam forming is optimized, including space-time signal filtering processing, Doppler processing, microbubble detection positioning and tracking processing and compensation processing of signal attenuation caused by skull reverberation. According to the method, the axial blood flow velocity estimation precision can be greatly improved, the spatial resolution and saturation of an imaging result are improved, and high-quality skull imaging and hemodynamics and brain function research are facilitated.

Description

Technical Field

[0001] The present invention belongs to the technical field of ultrasonic imaging, and in particular relates to a phase distortion correction method for transcranial brain ultrasonic imaging. Background Art

[0002] Transcranial ultrasound imaging is a technology that uses ultrasound to perform radiation-free, non-invasive cerebral blood flow imaging and brain function evaluation in the presence of the skull. Transcranial ultrasound Doppler technology measures blood flow velocity by detecting the Doppler frequency deviation of the echo received from red blood cells, and transcranial ultrasound localization microscopy technology reconstructs an accurate cerebral blood flow map by identifying, locating and tracking wandering ultrasound microbubble contrast agents.

[0003] In recent years, ultrafast ultrasound imaging technology has developed rapidly, driving the rapid advancement of ultrafast Doppler technology and ultrasound localization microscopy technology. Unlike traditional focused beam scanning, ultrafast ultrasound Doppler technology uses multi-angle composite plane wave technology to achieve a resolution of hundreds of microns without contrast agents, while ultrasound localization microscopy technology further breaks through the diffraction limit of sound waves and achieves blood vessel detection with a diameter of tens of microns. However, the imaging quality of ultrasound Doppler and ultrasound localization microscopy in transcranial conditions is still affected by the skull structure. There is a significant difference in acoustic impedance between the skull and soft tissue, and sound waves undergo non-negligible refraction at the interface, causing phase distortion, resulting in deviations in the position and shape of the final imaging result. Secondly, the skull also causes sound wave attenuation, making it difficult for high-frequency ultrasound signals to penetrate, and low-frequency excitation will lead to a decrease in imaging resolution. Summary of the invention

[0004] The purpose of this application is to provide a phase distortion correction method for transcranial brain ultrasound imaging, so as to significantly improve the spatial resolution and saturation of the imaging results, and to facilitate high-quality transcranial bone imaging and precise hemodynamic research and brain function research.

[0005] The phase distortion correction method for transcranial brain ultrasound imaging provided by the present invention comprises the following specific steps:

[0006] (1) Using a linear array probe or a phased array probe to emit any type of non-focused sound waves, obtain ultrasound echo data in a transcranial scene and a B-mode image obtained after conventional time-delayed stack beamforming;

[0007] (2) Identify the upper and lower boundaries of the skull in transcranial B-mode images using an improved accelerated Dijkstra algorithm or an improved U-Net deep learning network method, and construct a cranial sound velocity model;

[0008] (3) solving the acoustic wave travel time matrix under the current sound velocity model based on the ray theory method or the Eikonal equation method, and using the acoustic wave travel time matrix for beam synthesis to obtain an image after phase distortion correction; and

[0009] (4) Processing and optimizing the transcranial imaging results after phase distortion correction and beam synthesis, including calculating the Doppler frequency shift using the autocorrelation of frame data to obtain blood flow velocity information, as well as spatiotemporal signal filtering, Doppler processing, microbubble detection, positioning and tracking, and compensation for signal attenuation caused by skull reverberation.

[0010] Further:

[0011] In step (1), the conventional beamforming method refers to delay-addition beamforming, specifically, under the assumption that the sound speed distribution is uniform, the delay is calculated by the Euclidean distance between the array element and the imaging area, and corresponding dynamic focusing imaging, or delay-addition imaging, is performed, including transmit focusing and receive focusing.

[0012] In step (1), the echo data is acquired by transmitting any type of non-focused sound waves through a linear array probe or a phased array probe, including single-angle plane waves, multi-angle composite plane waves, undeflected divergent waves, composite divergent waves and synthetic aperture transmission; wherein the synthetic aperture transmission method includes single-element transmission and single-element reception, single-element transmission and multi-element reception, and multi-element transmission and multi-element reception.

[0013] In step (2), the identification of the skull boundary in the transcranial B-mode image specifically includes the following method:

[0014] Extraction and recognition based on traditional image processing methods, extraction and recognition based on improved accelerated Dijkstra algorithm, extraction and recognition based on improved U-Net deep learning network method; among them:

[0015] The traditional image processing method includes median filtering, Gaussian filtering, mean filtering, pixel gradient operator edge processing and threshold processing.

[0016] The improved accelerated Dijkstra algorithm is specifically as follows: a smaller weight is assigned to pixels with higher pixel values, and then a series of points on the edge are artificially selected at similar intervals, and the shortest path along the bright line in the B-Mode image is automatically found between these points by minimizing the sum of the pixel weights along the path, and finally a complete edge is obtained. When the selected point is not included in the processed point set, the loop stops, and only one path backtracking is performed.

[0017] The improved U-Net deep learning network method is specifically as follows: in order to obtain better real-time performance, the traditional VGG16 encoder in the U-Net structure is replaced with a lightweight MobileNetV2, and an inverted residual structure is adopted. First, the dimension is increased by 1×1 convolution, and then features are extracted using depthwise separable convolution. Finally, the dimension is reduced by 1×1 convolution, thereby reducing the amount of network calculation parameters without sacrificing accuracy.

[0018] In step (2), the construction of the cranial brain sound velocity model is specifically as follows:

[0019] In the B-mode imaging results, the upper surface of the skull is identified by using an improved accelerated Dijkstra algorithm or an improved U-Net deep learning network method; the sound velocity distribution is updated, the sound velocity of the area below the identified upper surface in the imaging area is assumed to be the sound velocity of bone tissue, and the travel time obtained by solving the Eikonal equation using the fast marching method is used for beam synthesis to obtain a beam synthesized signal;

[0020] In the B-mode imaging results, the lower surface of the skull is identified by using an improved accelerated Dijkstra algorithm or an improved U-Net deep learning network method; the sound velocity distribution is updated, the sound velocity of the area below the identified lower surface in the imaging area is assumed to be the water sound velocity, and the travel time obtained by solving the Eikonal equation using the fast marching method is used for beam synthesis to obtain a beam synthesized signal;

[0021] In step (3), the acoustic wave travel time matrix is ​​specifically: the imaging area is divided into a number of grids of equal length and width, and the acoustic wave starts from the grid where the geometric center of the sound source corresponding to the array element is located, and the matrix is ​​composed of the time required for the first arrival wave to reach all the grids in the imaging area.

[0022] In step (3), solving the acoustic wave travel time matrix under the current sound velocity model based on the ray theory method or the Eikonal equation method specifically includes:

[0023] The travel time is obtained by solving the Eikonal equation using the fast marching method (Sethian, JA, Fast marching methods. SIAM review, 1999. 41(2): p. 199-235.), and the travel time is obtained by solving the ray path using the ray equation combined with the shooting method.

[0024] In step (3), the phase distortion refers to the fact that due to the significant difference in sound speed, refraction cannot be ignored in the transmission of sound waves, and the delay of sound waves calculated according to the assumption of uniform medium will result in an erroneous first arrival time, which will lead to distortion of the imaging result.

[0025] In step (3), the travel time calculation is performed using the fast marching method, specifically: the travel time field t(x, z) and the velocity field c(x, z) satisfy the eikonal equation:

[0026]

[0027] The initial conditions of the equation are:

[0028] t(x0,z0)=0, (2)

[0029] Among them, x0, z0 represent the coordinates of the sound source. In order to obtain the solution of the travel time field, it is necessary to solve the equation:

[0030]

[0031] Among them, s ij represents the slowness field, which is the velocity field c ij The reciprocal of D -a and D +a They are the forward difference and backward difference operators in the a direction respectively. In order to ensure accuracy, a second-order three-point finite difference operator is selected. The forward difference operator in the x direction is:

[0032]

[0033] Substituting the difference operator, we can get:

[0034]

[0035] in:

[0036]

[0037] Similarly, the forward difference operator in the z direction is:

[0038]

[0039] The backward difference operator in the x direction is:

[0040]

[0041] The backward difference operator in the z direction is:

[0042]

[0043] Therefore, the original eikonal equation can be modified to i,j The quadratic equation of .

[0044] In transcranial imaging, after the travel time field is obtained, beam synthesis can be performed using the travel time field. s =(x s ,z s), the beamformed signal can be written as:

[0045] s DelayAndSum =∑s(τ TX (X s )+τ RX (X s ,j)), (10)

[0046] Among them, τ RX (X s ,j) represents the sound wave from point X s The travel time of the wavefront to the scattering point is:

[0047] τ TX (X s )=min j＝1,2,...,n τ RX (X s ,j), (11)

[0048] Among them, n is the number of array elements.

[0049] In step (4), the Doppler frequency shift is calculated by using the autocorrelation of the frame data, thereby calculating the blood flow velocity information, specifically: by analyzing the echo signal in the imaging area, the blood flow velocity estimation and volume blood flow measurement are realized. The autocorrelation estimation method can convert the frequency shift calculation in the frequency domain into the autocorrelation calculation in the time domain, and the Doppler frequency shift is calculated from the time domain signal. The frame data autocorrelation is:

[0050]

[0051] Among them, s(x,y,k) is the blood flow echo signal of the kth frame, * represents the conjugate operation, N represents the total number of frames processed, and Δt represents the interval time between adjacent frames. Therefore, the blood flow velocity information can be calculated as:

[0052]

[0053] Where |v| represents the target velocity, θ represents the angle between the ultrasound beam and the flow vector, c0 represents the soft tissue sound velocity, PRF represents the pulse repetition frequency, f0 represents the center frequency of the emission, imag(R) represents the imaginary part of R, and real(R) represents the real part of R.

[0054] In step (4), the transcranial imaging results after beam synthesis are processed and optimized, including spatiotemporal signal filtering, Doppler processing, microbubble detection, positioning and tracking processing, and compensation processing for microbubble signal attenuation in the cerebral cortex, specifically:

[0055] Spatiotemporal signal filtering, including SVD filtering, EVD filtering, and Cauchy-RPCA filtering;

[0056] Doppler processing, including power Doppler processing, color Doppler processing, and vector Doppler processing;

[0057] Compensation for microbubble signal attenuation in the cerebral cortex was performed, including a depth-dependent intensity gain factor method.

[0058] In step (4), the SVD method in the spatiotemporal signal filtering is specifically:

[0059] The echo signal received by the transducer is represented as a set of three-dimensional data sets s(x,z,t), where x represents the coordinates of the receiving array element, z represents the imaging depth, t represents the time, and the number of imaging frames accumulates over time. The echo signal received by the ultrasonic transducer consists of three parts: clutter signal c(x,z,t), blood flow signal b(x,z,t) and noise signal n(x,z,t), which can be expressed as follows:

[0060] s(x,z,t)=c(x,z,t)+b(x,z,t)+n(x,z,t), (14)

[0061] The above three signals have different time-frequency domain characteristics. The noise signal can be regarded as zero-mean Gaussian white noise. The tissue signal has a low frequency and the blood flow signal has a high frequency. Compared with the blood flow signal that changes rapidly and frequently over time, the tissue signal has a stronger spatial coherence. The three-dimensional space-time matrix s(x,z,t) is obtained by sampling through ultrafast ultrasound imaging. The data volume is n x ×n z ×n t , where n x is the number of sampling points in the x-axis direction, n z is the number of sampling points in the z-axis direction, n t is the number of sampling frames in the time domain. Next, the three-dimensional data matrix is ​​converted into a two-dimensional spatiotemporal matrix S, and the matrix dimension is (n x ×n z ,n t ).

[0062] The singular value decomposition method is used for analysis. First, the matrix S is decomposed by singular value:

[0063] S=UΔV * , (15)

[0064] Among them, U and V are standard orthogonal matrices with dimensions (n x ×n y ,n x ×n y ) and (nt ×n t ), V * Represents the conjugate transposed matrix of V. Each column of the U and V matrices and They correspond to the spatial and temporal singular vectors of the S matrix, respectively, as follows:

[0065]

[0066] Analyze the components of the echo signal and regard the matrix S as several matrices S i The weighted result, where S i for and The outer product of, the expression of S matrix is:

[0067]

[0068] In ultrasound blood flow imaging, tissue echo signals have greater energy and stronger spatial coherence than blood flow echo signals. Therefore, after eigenvalue decomposition, the components of tissue echo signals correspond to the first few eigenvectors. Therefore, a threshold can be set to delete all eigenvalues ​​with amplitudes greater than the threshold and filter out tissue echo signals. The signal S after filtering is f for:

[0069] S f =UΔ f V * , (19)

[0070] Among them, Δ f is the result of setting the first N largest singular values ​​in Δ to zero.

[0071] In step (4), the depth-dependent intensity gain factor compensation is specifically: due to the imaging artifacts caused by multiple reflections caused by the skull, the microbubble signal is weak, and it is difficult to detect microbubbles close to the skull. This ultimately leads to the inability to reconstruct a large number of cerebral cortical blood vessels. In the present invention, compensation is adopted, that is, the depth-dependent intensity gain control factor is designed before detection and positioning. Specifically, these image batches are divided into several parts in the axial direction (depth direction), and the gain of each part is determined by the depth gain factor. A larger gain factor is designed for a position close to the cortex, and a smaller gain factor is designed for other areas.

[0072] Function and effect of the present invention

[0073] The present invention can be used for phase distortion correction of transcranial ultrasound Doppler and transcranial ultrasound localization microscopy, as well as transcranial hemodynamic research. After skull recognition and phase distortion correction, the accuracy of axial blood flow velocity estimation is greatly improved, and the spatial resolution and saturation of imaging results are significantly improved. At the same time, the cortical blood vessels affected by reverberation can be successfully reconstructed after gain compensation. The present invention is conducive to high-quality transcranial bone imaging, precise hemodynamic research, and brain function research. BRIEF DESCRIPTION OF THE DRAWINGS

[0074] Figure 1 The present invention is a flow chart of a phase distortion correction method applied to transcranial brain ultrasound imaging.

[0075] Figure 2 It is the improved U-Net network structure and segmentation schematic diagram corresponding to the present invention for segmenting the upper and lower surfaces of the skull.

[0076] Figure 3 This is the process of skull recognition and sound velocity distribution update in the phantom skull and phantom blood flow experiment in the embodiment of the present invention.

[0077] Figure 4 It is the visualization result of the acoustic wave travel time matrix described in the present invention.

[0078] Figure 5 These are the ultrafast power Doppler and ultrafast color Doppler experimental results of the phantom skull and phantom blood flow experiments in the embodiments of the present invention.

[0079] Figure 6 This is the calibrated experimental result of transcranial ultrasound localization microscopic imaging of the rat brain in the embodiment of the present invention.

[0080] Figure 7 These are the experimental results of transcranial ultrasound Doppler and transcranial ultrasound localization microscopic imaging of the rat brain before and after correction in the embodiments of the present invention. DETAILED DESCRIPTION

[0081] The present invention is further described below by way of embodiments in conjunction with the accompanying drawings.

[0082] This embodiment uses a phantom blood flow and a phantom skull to conduct experiments, and performs transcranial ultrafast ultrasound Doppler imaging on the phantom blood flow. The experiment was conducted on a homemade blood flow phantom, and two rubber hoses with a radius of 2 mm were prepared, one of which had a hole cut in the middle, and the other hose was inserted into the hole to form a y-shaped bifurcation. This bifurcated structure was placed in a custom box, and the three ends extended from the circular hole of the box. Subsequently, a mixture of agar and hot water was poured into the setting at a weight ratio of 1:24, allowed to cool, and then solidified. Two rubber hoses were drawn out of the agar model to form a vascular phantom. A 6% starch aqueous solution was used as simulated blood, injected from the bifurcation end, and flowed out from the confluence end. In the experiment, the composite emission frame rate of the multi-angle plane wave was set to 400 frames / second, and the sequence consisted of plane wave emissions at 15 angles per frame, with an angle range of -12° to 12° and an interval of 1.6°. The experiment obtained 200 frames of RF data for beamforming. In order to extract blood flow signals, SVD filtering and eigenvector analysis were performed on 200 frames of image data to remove static tissue signals. The eigenvalue range set in the SVD clutter filter is [15,150], in which the static information corresponding to the maximum 15 eigenvalues ​​and the noise information corresponding to the minimum 50 eigenvalues ​​are filtered out. The image data after SVD filtering was demodulated in phase and orthogonally, and the average Doppler frequency shift of each sampling point was calculated. Finally, the velocity of the scattering point was calculated to obtain the power Doppler image and color Doppler image of the blood flow.

[0083] This embodiment uses adult rats for experiments. After anesthesia, the rats are fixed on a rat brain stereotaxic instrument, and transcranial ultrasound localization microscopy imaging of the rat brain is performed. 7-week-old male Sprague-Dawley rats were selected and placed in a small animal stereotaxic device to restrict head movement and maintain body temperature with a thermal blanket. 2.5% isoflurane was continuously induced and maintained anesthesia to ensure the stability of data acquisition. Only the hair on the skin was removed to ensure adequate coupling and prevent acoustic artifacts in the imaging results. The skin and skull of the rats were intact and covered with ultrasound gel. Important parameters (body temperature, respiratory rate) were checked regularly. The acquisition time was 7min, and 840μL of microbubble solution was injected into the jugular catheter at a constant rate of 2μL / s to keep the microbubble concentration in the rat blood volume constant. The acquisition sequence uses a plane wave imaging sequence with a composite frame rate of 1000Hz. In order to match the width of the rat brain, only the central 80 elements of the probe L11-5v were activated. Seven tilted plane waves with angles ranging from -10° to 10° were transmitted at a pulse repetition rate of 7 kHz, close to the physical frame rate limit. The entire acquisition process consisted of 420 buffers in 7 minutes. Each buffer consisted of 300 frames collected in 0.3 seconds.

[0084] Figure 2The figure shows the improved U-Net architecture used in this study to achieve the segmentation of the interface between skull and soft tissue. In order to obtain better real-time performance, the traditional VGG16 encoder in the U-Net structure is replaced with the lightweight MobileNetV2, which reduces the network calculation parameters to about 3% of the original without sacrificing accuracy. The inverted residual structure is used in MobileNetV2. The structure first increases the dimension through 1×1 convolution, then uses depthwise separable convolution to extract features, and finally reduces the dimension through 1×1 convolution. The purpose of this design is to extract more features in high-dimensional space through depthwise separable convolution after expanding the number of channels, thereby increasing the network depth and complexity without significantly increasing the computational cost. Depthwise separable convolution can be regarded as decomposing the standard convolution into depthwise convolution and pointwise convolution. Compared with traditional convolution, the amount of calculation is reduced by several times, which provides support for the lightweight design of the network. At the end of the inverted residual block, the output low-dimensional feature part uses a linear activation function, which helps to avoid the information loss and destruction that may be caused by the use of nonlinear activation function ReLU in the process of high-dimensional to low-dimensional conversion.

[0085] Figure 5 The results show that the accuracy and quality of Doppler blood flow imaging are significantly improved after phase distortion correction. The blood flow centerline extracted based on morphology, black is the result without skull, blue is the result of phase distortion, and red is the result after correction. Taking the result without skull as reference, the average normalized root mean square error of the axial velocity component on the five axial sections after correction is reduced from 17.67% to 8.02%.

[0086] Figure 7 The left image shows that due to the relatively thick skull of adult rats, the resolution of transcranial Doppler imaging is limited and most blood vessels cannot be clearly observed. Figure 7 In the middle and right images, more blood vessels can be identified using the transcranial localization microscopy method. After aberration correction, the location of the microbubbles is more precise and the trajectory is more accurate. Figure 7 In the right image, part of the sagittal sinus (a vein along the midline close to the cortex) is reconstructed after aberration correction, and the intensity of the cerebral cortical blood vessels is higher. Under the same parameters, the corrected transcranial localization microscopy spatial resolution is improved from 38μm to 35μm, indicating that the phase correction method can improve the spatial resolution and image saturation of transcranial localization microscopy.

Claims

1. A phase distortion correction method for transcranial brain ultrasound imaging, characterized in that: The specific steps are: (1) Using a linear array probe or a phased array probe to emit any type of non-focused sound waves, obtain ultrasound echo data in a transcranial scene and a B-mode image obtained after conventional time-delayed stack beamforming; (2) Identify the upper and lower boundaries of the skull in transcranial B-mode images using an improved accelerated Dijkstra algorithm or an improved U-Net deep learning network method, and construct a cranial sound velocity model; (3) Based on the ray theory method or the Eikonal equation method, the acoustic wave travel time matrix under the current sound speed model is solved, and the acoustic wave travel time matrix is ​​used for beam synthesis to obtain an image after phase distortion correction; as well as (4) Processing and optimizing the transcranial imaging results after phase distortion correction and beam synthesis, including calculating the Doppler frequency shift using the autocorrelation of frame data to obtain blood flow velocity information, as well as spatiotemporal signal filtering, Doppler processing, microbubble detection, positioning and tracking, and compensation for signal attenuation caused by skull reverberation.

2. The phase distortion correction method according to claim 1, characterized in that: In step (1): The conventional beamforming method refers to delay-and-addition beamforming. Specifically, under the assumption that the sound velocity distribution is uniform, the delay is calculated by the Euclidean distance between the array element and the imaging area, and corresponding dynamic focusing imaging, or delay-and-addition imaging, is performed, including transmission focusing and reception focusing. The echo data is acquired by transmitting any type of non-focused sound waves through a linear array probe or a phased array probe, including single-angle plane waves, multi-angle composite plane waves, undeflected divergent waves, composite divergent waves and synthetic aperture transmission; wherein the synthetic aperture transmission method includes single-element transmission and single-element reception, single-element transmission and multi-element reception, and multi-element transmission and multi-element reception.

3. The phase distortion correction method according to claim 2, characterized in that: In step (2), the identification of the skull boundary in the transcranial B-mode image specifically includes: Extraction and recognition based on traditional image processing methods, extraction and recognition based on improved accelerated Dijkstra algorithm, extraction and recognition based on improved U-Net deep learning network method; among them: The traditional image processing methods include median filtering, Gaussian filtering, mean filtering, pixel gradient operator edge processing and threshold processing; The improved accelerated Dijkstra algorithm is specifically as follows: a smaller weight is assigned to pixels with higher pixel values, and then a series of points on the edge are artificially selected at similar intervals, and the shortest path along the bright line in the B-Mode image is automatically found between these points by minimizing the sum of the pixel weights along the path, and finally a complete edge is obtained; when the selected point is not included in the processed point set, the loop stops, and only one path backtracking is performed; The improved U-Net deep learning network method is specifically as follows: in order to obtain better real-time performance, the traditional VGG16 encoder in the U-Net structure is replaced with a lightweight MobileNetV2, and an inverted residual structure is adopted. First, the dimension is increased by 1×1 convolution, and then features are extracted using depthwise separable convolution. Finally, the dimension is reduced by 1×1 convolution, thereby reducing the amount of network calculation parameters without sacrificing accuracy.

4. The phase distortion correction method according to claim 3, characterized in that: The step (2) of constructing the cranial brain sound velocity model is specifically as follows: Identify the upper surface of the skull in the B-mode imaging results by using an improved accelerated Dijkstra algorithm or an improved U-Net deep learning network method; Update the sound velocity distribution, assume the sound velocity of the area below the identified upper surface in the imaging area to be the sound velocity of bone tissue, use the fast marching method to solve the Eikonal equation to obtain the travel time for beam synthesis, and obtain the beam synthesized signal; Identify the lower surface of the skull in the B-mode imaging results by using an improved accelerated Dijkstra algorithm or an improved U-Net deep learning network method; The sound velocity distribution is updated, and the sound velocity in the area below the identified lower surface in the imaging area is assumed to be the water sound velocity. The travel time obtained by solving the Eikonal equation using the fast marching method is used for beam synthesis to obtain the beam synthesized signal.

5. The phase distortion correction method according to claim 4, characterized in that: In step (3): The acoustic wave travel time matrix is ​​specifically: the imaging area is divided into a number of grids of equal length and width, and the acoustic wave starts from the grid where the geometric center of the sound source corresponding to the array element is located, and the matrix is ​​composed of the time required for the first arrival wave to reach all the grids in the imaging area; The method based on the ray theory method or the Eikonal equation method is used to solve the acoustic wave travel time matrix under the current sound velocity model, specifically including: The travel time is solved by solving the Eikonal equations by the fast marching method, and the travel time is solved by solving the ray path by the ray equations combined with the shooting method; The phase distortion refers to the fact that due to the significant difference in sound speed, refraction cannot be ignored in sound wave transmission, and the delay of the sound wave calculated according to the uniform medium assumption will result in an erroneous first arrival time, which will lead to distortion of the imaging result.

6. The phase distortion correction method according to claim 5, characterized in that: The travel time calculation using the fast marching method described in step (3) is as follows: the travel time field t(x, z) and the velocity field c(x, z) satisfy the eikonal equation: The initial conditions of the equation are: t(x0,z0)=0, (2) Among them, x0, z0 represent the coordinates of the sound source; in order to obtain the solution of the travel time field, solve the equation: Among them, s ij represents the slowness field, which is the velocity field c ij The reciprocal of D -a and D +a They are the forward difference and backward difference operators in the a direction respectively; the second-order three-point finite difference operator is selected specifically, and the forward difference operator in the x direction is: Substituting the difference operator, we get: in: Similarly, the forward difference operator in the z direction is: The backward difference operator in the x direction is: The backward difference operator in the z direction is: Therefore, the original eikonal equation is modified to i，j The quadratic equation of In transcranial imaging, the travel time field is obtained and then used to perform beam synthesis. s =(x s ,z s ), the beamformed signal is written as: s DelayAndSum =∑s(T TX (X s )+τ RX (X s ,j)), (10) Among them, τ RX (X s ,j) represents the sound wave from point X s The travel time of the wavefront to the scattering point is: T TX (X s )=min j＝1，2，...，n t RX (X s ,j), (11) Among them, n is the number of array elements.

7. The phase distortion correction method according to claim 6, characterized in that: The step (4) uses the autocorrelation of frame data to calculate the Doppler frequency shift, thereby calculating the blood flow velocity information. Specifically, the blood flow velocity estimation and volume blood flow measurement are achieved by analyzing the echo signal in the imaging area; the autocorrelation estimation method converts the frequency shift calculation in the frequency domain into the autocorrelation calculation in the time domain, and the Doppler frequency shift is calculated from the time domain signal; the frame data autocorrelation is: Among them, s(x, y, k) is the blood flow echo signal of the kth frame, * represents the conjugate operation, N represents the total number of frames processed, and Δt represents the interval time between adjacent frames; therefore, the blood flow velocity information is calculated as: Where |v| represents the target velocity, θ represents the angle between the ultrasound beam and the flow vector, c0 represents the soft tissue sound velocity, PRF represents the pulse repetition frequency, f0 represents the center frequency of the emission, imag(R) represents the imaginary part of R, and real(R) represents the real part of R.

8. The phase distortion correction method according to claim 7, characterized in that: The processing and optimization of the transcranial imaging results after beam synthesis described in step (4) include spatiotemporal signal filtering, Doppler processing, microbubble detection, positioning and tracking processing, and compensation processing for microbubble signal attenuation in the cerebral cortex, specifically: Spatiotemporal signal filtering, including SVD filtering, EVD filtering, and Cauchy-RPCA filtering; Doppler processing, including power Doppler processing, color Doppler processing, and vector Doppler processing; Compensation for microbubble signal attenuation in the cerebral cortex was performed, including a depth-dependent intensity gain factor method.

9. The phase distortion correction method according to claim 8, characterized in that: In step (4): The SVD method in the spatiotemporal signal filtering is specifically as follows: the echo signal received by the transducer is represented as a set of three-dimensional data sets s(x, z, t), where x represents the coordinates of the receiving array element, z represents the imaging depth, t represents the time, and the number of imaging frames is accumulated over time; the echo signal received by the ultrasonic transducer consists of three parts: a clutter signal c(x, z, t), a blood flow signal b(x, z, t) and a noise signal n(x, z, t), which can be expressed as the following formula: s(x,z,t)=c(x,z,t)+b(x,z,t)+n(x,z,t), (14) The above three signals have different time-frequency domain characteristics. The noise signal is regarded as zero-mean Gaussian white noise. The tissue signal has a low frequency, the blood flow signal has a high frequency, and has stronger spatial coherence. The three-dimensional space-time matrix s(x,z,t) is obtained by sampling through ultrafast ultrasound imaging. The data volume is n x ×n z ×n t , where n x is the number of sampling points in the x-axis direction, n z is the number of sampling points in the z-axis direction, n t is the number of sampling frames in the time domain; next, the three-dimensional data matrix is ​​converted into a two-dimensional space-time matrix S, and the matrix dimension is (n x ×n z ,n t ); The singular value decomposition method is used for analysis. First, the matrix S is decomposed by singular value: S=UΔV * , (15) Among them, U and V are standard orthogonal matrices with dimensions (n x ×n y ,n x ×n y ) and (n t ×n t ), V * Represents the conjugate transposed matrix of V; each column of the U and V matrices and They correspond to the spatial and temporal singular vectors of the S matrix, respectively, as follows: Analyze the components of the echo signal and regard the matrix S as several matrices S i The weighted result, where S i for and The outer product of , the expression of S matrix is: In ultrasound blood flow imaging, tissue echo signals have greater energy and stronger spatial coherence than blood flow echo signals. Therefore, after eigenvalue decomposition, the components of tissue echo signals correspond to the first few eigenvectors. Therefore, a threshold is set to delete all eigenvalues ​​with amplitudes greater than the threshold and filter out tissue echo signals. The signal S after filtering is f for: With f =-UΔ f In * , (19) Among them, Δ f is the result of setting the first N largest singular values ​​in Δ to zero.

10. The phase distortion correction method according to claim 9, characterized in that: The depth-dependent intensity gain factor compensation described in step (4) is to design a depth-dependent intensity gain control factor before detection and positioning. Specifically, these image batches are divided into several parts in the axial direction, and the gain of each part is determined by the depth gain factor. A larger gain factor is designed for a position close to the cortex, and a smaller gain factor is designed for other areas.

Images