CN114972567A - Medical ultrasonic CT multi-parameter image reconstruction method based on wave equation - Google Patents
Medical ultrasonic CT multi-parameter image reconstruction method based on wave equation Download PDFInfo
- Publication number
- CN114972567A CN114972567A CN202210604492.XA CN202210604492A CN114972567A CN 114972567 A CN114972567 A CN 114972567A CN 202210604492 A CN202210604492 A CN 202210604492A CN 114972567 A CN114972567 A CN 114972567A
- Authority
- CN
- China
- Prior art keywords
- parameter
- deconvolution
- sound pressure
- wave equation
- pressure signal
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T11/00—2D [Two Dimensional] image generation
- G06T11/003—Reconstruction from projections, e.g. tomography
- G06T11/005—Specific pre-processing for tomographic reconstruction, e.g. calibration, source positioning, rebinning, scatter correction, retrospective gating
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T3/00—Geometric image transformation in the plane of the image
- G06T3/40—Scaling the whole image or part thereof
- G06T3/4053—Super resolution, i.e. output image resolution higher than sensor resolution
Abstract
The invention provides a wave equation-based medical ultrasonic CT multi-parameter image reconstruction method, which comprises the following steps of: transmitting and receiving ultrasound by using ultrasound acquisition equipment to acquire an observation sound pressure signal of a receiving array element in a real three-dimensional space; simulating the process by using a multi-parameter viscoacoustic wave equation to obtain a synthetic sound pressure signal corresponding to a receiving array element in a real three-dimensional space; comparing the difference between the observed sound pressure signal and the synthesized sound pressure signal by using a deconvolution target function and a two-norm target function, and solving a first-order gradient; solving a new parameter matrix by using the first-order gradient; and inputting the new parameter matrix into the multi-parameter viscoacoustic wave equation to carry out iterative operation, and obtaining the next new parameter matrix to obtain a multi-parameter reconstructed image meeting the conditions. The method adopts an FWI algorithm, introduces a multi-parameter viscoacoustic wave equation and a deconvolution target function, and utilizes full waveform information to accurately reconstruct images of various human tissue parameters to obtain a high-resolution medical ultrasonic image superior to a traditional ray method.
Description
Technical Field
The application relates to the technical field of medical ultrasonic imaging, in particular to a medical ultrasonic CT multi-parameter image reconstruction method based on a wave equation.
Background
The conventional ultrasound CT image reconstruction is generally classified into echo intensity image reconstruction based on echo and sound velocity and attenuation coefficient reconstruction based on transmitted wave. The echo intensity image is of higher resolution but can only assist the physician in diagnosing when there is significant lesion in the tissue, a lesion with differences in acoustic properties has formed. The transmitted wave method can reconstruct functional parameters such as sound velocity, attenuation coefficient and the like, has important significance for monitoring the early functional parameter change of the pathological tissue, but has lower resolution, and is difficult to distinguish the characteristic difference between the background tissue and the pathological tissue under extreme conditions. Therefore, the reconstruction of high-resolution functional parameters in combination with multiple types of acoustic signals, including but not limited to echo and transmitted waves, is of great significance for early imaging and diagnosis of lesions.
The ultrasonic CT image reconstruction of the multi-type signals has two methods based on ray theory and fluctuation theory. The multi-mode reconstruction based on the ray theory generally adopts a method for fusing the two images, which is limited by the difference between the resolutions of the two images, and the fusion process needs proper processing, otherwise, wrong focus information may be introduced. The reconstruction based on the wave equation can directly utilize all waveform data to reconstruct multiple parameter images such as sound velocity, acoustic impedance, attenuation coefficient and the like, has higher resolution and better effect of assisting in distinguishing focuses. The imaging method of Full Waveform Inversion (FWI) combines the calculated amount of the leading edge computing device within an acceptable range, and properly designs the processing flow to ensure the robustness of the method, so that the method is a stable, practical and high-resolution multi-feature image reconstruction method at present.
Due to the fact that a complicated multi-parameter wave equation (such as sound velocity, acoustic impedance, attenuation and the like) needs to be introduced into the FWI-based multi-parameter image reconstruction method, the sound velocity single-parameter reconstruction is mainly focused on the FWI-based ultrasonic CT image reconstruction method at home and abroad at present. Moreover, FWl is a local optimization method, and has the problem of local minima, which affects the accuracy of the results. Therefore, a precise FWI ultrasonic CT image reconstruction method based on a multi-parameter wave equation is needed to obtain a high-resolution multi-parameter image to assist lesion identification and diagnosis.
Disclosure of Invention
The application aims to solve the problem of obtaining high-resolution multi-parameter ultrasound CT accurate images based on FWI.
In order to achieve the above object, the present application provides a wave equation-based medical ultrasound CT multi-parameter image reconstruction method, comprising the following steps:
transmitting and receiving ultrasound by using ultrasound CT data acquisition equipment to acquire real three-dimensional space data; the data is an observed sound pressure signal of the position of the receiving array element in the real three-dimensional space;
simulating the process of transmitting and receiving ultrasound by the acquisition equipment by using a multi-parameter viscoacoustic dynamic equation to obtain a synthetic sound pressure signal corresponding to the position of the receiving array element in the real three-dimensional space; the input variable of the multi-parameter viscous sound wave equation is a parameter matrix consisting of an original sound velocity graph, an original acoustic impedance graph and an original attenuation quality factor graph;
solving a difference between the observed sound pressure signal and the synthesized sound pressure signal by using a deconvolution target function and a two-norm target function;
calculating the deconvolution target function and the two-norm target function by a adjoint state method to obtain a first-order gradient of the parameter matrix; calculating to obtain a new parameter matrix based on the first-order gradient and the parameter matrix; converting the new parameter matrix into a new sound velocity diagram, a new acoustic impedance diagram and a new attenuation quality factor diagram;
inputting the new parameter matrix into the multi-parameter viscoacoustic wave equation, and combining the deconvolution target function and the first-order gradient to perform iterative operation to obtain a next new parameter matrix; the iterative operation stop condition is that the synthesized sound pressure signal and the observed sound pressure signal are consistent.
Preferably, the ultrasonic CT data acquisition equipment is provided with a multilayer annular array, and the multilayer array simultaneously acquires data after single excitation to acquire three-dimensional volume data with a longitudinal aperture;
the data acquisition equipment is an observation array transducer containing a plurality of array elements; each array element of the plurality of array elements can transmit and receive ultrasonic signals;
optionally, the observation array transducer is selected from an annular probe, a through linear array probe, a through planar array probe and a through concave array probe.
Preferably, the observed sound pressure signal of the position of each of the plurality of array elements is d obs (x r ，t；x s ) Wherein x is s Is the coordinate position, x, of the transmitting array element r T represents the reception time for the coordinate position of the receiving array element.
Preferably, the multi-parameter viscoacoustic wave equation is a forward function, and the three-dimensional equation formula is as follows:
wherein the proton velocity describing the acoustic field is v i Stress of σ ij The synthesized sound pressure p ═ Σ k σ kk A/3; the parameters describing the physical properties of the medium are rho, sound speed V p Acoustic impedance I p ＝ρV p Described, the attenuation coefficient is determined by the quality factor Q through the hysteresis coefficient
Introducing that the hysteresis elastic coefficient adopts L at a reference frequency
Memory variable of
If attenuation is not considered, the last two equations and
the coefficients of the correlations are truncated and degenerate back to the conventional isotropic acoustic wave equation, whose two-dimensional form can be formed by truncating the z-related term, the variable v i ，σ ij ，ρ，V p ，
Are both images defined in the computed imaging space x and time t;
given a transmit array element (position x) s ) The information of (2) is input into a parameter matrix of a sound velocity diagram, an acoustic impedance diagram and an acoustic attenuation diagram, sound pressure information p (x, t) on the whole calculation space x and time t can be solved discretely through a numerical simulation method (such as finite difference, finite element and the like), and then a receiving array element x is extracted r Synthesized sound pressure signal d cal (x r ，t；x s )＝R(x，x r )p(x，t；x s ) The maximum recording time T should be such that the transmitted wave signal through the target body tissue resulting from the farthest apart transmit and receive array element combination can be recorded.
Preferably, the deconvolution target function is an inversion function, so that the result of the FWI algorithm is closer to a global minimum, a time domain wiener filter between the measured data and the simulated data is solved, the target function is designed by weighting the filter, the transit time difference between the two data based on the wave equation can be effectively highlighted, the influence of inaccurate forward amplitude on inversion is weakened, and the deconvolution target function can be expressed as:
wherein A (tau) is a weighting function, w represents a wiener filter,
in inversion, the selection of the weight function A (tau) is crucial, and the weight function can determine the robustness of the algorithm, | tau |/tau |, and max 、τ 2 and gaussian functions may be chosen. The multi-scale inversion can be realized by setting an upper limit of the tau value in the weight function number, because the upper limit represents the maximum transit time difference between the observed data and the simulation data which can be considered by the objective function. The upper limit ofShould combine frequency division multiscale inversion selection, between 3-5 times of reciprocal of dominant frequency of corresponding frequency band data, because frequency is gradually increased, the effect of gradually decreasing upper limit of tau is also realized, and then the traditional two-norm objective function C is utilized l2 The iteration increases the resolution of the image,
the conventional two-norm objective function still needs to be improved in the aspects of stability, accuracy, high efficiency and the like: under the conditions that low-frequency information is lacked in data and an initial model is far away from a real model, an optimization algorithm is difficult to effectively converge to a global minimum value to obtain accurate ultrasonic CT image reconstruction.
Preferably, said deconvolution objective function is said first order gradient with respect to said parameter matrix
It can be calculated by the adjoint state method,
wherein the adjoint wave field operator mu is an equation
The solution of (1). Changing the deconvolution objective function does not affect the calculation of the forward wavefield w. However, the adjoint wavefield μ is transmitted by reverse propagation
And obtaining the partial derivative, wherein the partial derivative changes with the change of the deconvolution objective function. In the case of a two-norm objective function, the partial derivative is:
under deconvolution conditions, the partial derivative is:
will be provided with
The adjoint wave field operator mu can be calculated by wave equation back propagation.
Preferably, the iteration can be regarded as a process of solving an optimization problem, and the initial parameter matrix m is input in the 0 th iteration 0 I.e. an estimation of the acoustic parameter model based on a priori knowledge;
alternatively, in the case of detection of soft tissue such as breast, it may be assumed that the parameter matrix m is initiated 0 Equivalent to a homogeneous aqueous medium.
Preferably, the iterative updating generally adopts a gradient descent method, and the values of the parameter matrix are iteratively updated along a descent direction d,
m k+1 ＝m k +α k d k ，
wherein alpha is k Typically obtained using a line search method, and d k There are different calculation methods in different optimization strategies (e.g., steepest descent method, conjugate gradient method, L-BFGS method, Newton method), and the steepest descent method generally uses the negative direction of the objective function with respect to the gradient of the model as the descent direction, i.e., the steepest descent method
The conjugate gradient method, the L-BFGS method and other methods need to improve the current gradient by utilizing gradient information in the pre-iteration to obtain the descending direction, and the Newton method needs to estimate a Hatheri matrix, namely a second derivative of an objective function about the gradient, to restrict the current gradient to obtain the descending direction.
Preferably, the iterative operation stop condition may be a threshold value at which the deconvolution objective function value decreases or a number of iterations.
Preferably, the codes for solving the multi-parameter viscoacoustic wave equation and the deconvolution objective function can be verified by a CPU and then transplanted to a GPU to improve efficiency so as to meet the requirement of a time-sensitive ultrasonic CT medical imaging clinical diagnosis task.
The application provides a medical ultrasound CT high-resolution image reconstruction method based on a wave equation, introduces a viscoacoustic wave equation (including sound velocity, acoustic impedance and attenuation parameters) and a deconvolution target function aiming at the problem of accurate multi-parameter image reconstruction, provides a method for accurately reconstructing various human tissue parameter images by utilizing full waveform information from the angle of an algorithm, ensures high-resolution image quality which is obviously superior to that of a traditional ray method, provides an auxiliary diagnosis basis of iconography and functional parameters for focus diagnosis, and comprises but not limited to a scene which is difficult to diagnose in an early stage by a traditional standard molybdenum target method, such as screening micro tumors in dense mammary glands.
Drawings
In order to more simply explain the technical solution of the embodiment of the present invention, the drawings needed to be used in the description of the embodiment will be briefly introduced below, and it is obvious that the drawings in the following description are only some embodiments of the present invention, and it is obvious for those skilled in the art to obtain other drawings based on these drawings without creative efforts.
FIG. 1 is a schematic flow chart of a wave equation-based medical ultrasound CT multi-parameter image reconstruction method in an embodiment of the present application;
FIG. 2 is a flow chart of an FWI-based ultrasound CT image reconstruction process according to an embodiment of the present application;
fig. 3 is a human tissue phantom and FWI-based ultrasound CT multi-parameter image in an embodiment of the present application.
Detailed Description
To make the objects, technical solutions and advantages of the embodiments of the present invention clearer and more complete, the technical solutions in the embodiments of the present invention will be described below with reference to the drawings in the embodiments of the present invention, it is obvious that the described embodiments are a part of the embodiments of the present invention, but not all embodiments, and all other embodiments obtained by those skilled in the art without any inventive work belong to the scope of the present invention
In fig. 1, the steps of the medical ultrasound CT multi-parameter image reconstruction method based on the wave equation provided by the present invention are as follows:
step S110: transmitting and receiving ultrasound by using ultrasound CT data acquisition equipment to acquire real three-dimensional space data; the data is an observation sound pressure signal of the position of the receiving array element in a real three-dimensional space;
step S120: simulating the process of transmitting and receiving ultrasound by using a multi-parameter viscoacoustic wave kinetic equation to acquire a synthetic sound pressure signal corresponding to the position of a receiving array element in a real three-dimensional space; the input variable of the multi-parameter viscous sound wave equation is a parameter matrix composed of an original sound velocity graph, an original acoustic impedance graph and an original attenuation quality factor graph;
step S130: solving a difference between the observed sound pressure signal and the synthesized sound pressure signal by using a deconvolution target function and a two-norm target function;
step S140: calculating the deconvolution target function and the two-norm target function through a companion state method to obtain a first-order gradient of the parameter matrix; calculating to obtain a new parameter matrix based on the first-order gradient and the parameter matrix; converting the new parameter matrix into a new sound velocity diagram, a new acoustic impedance diagram and a new attenuation quality factor diagram;
step S150: inputting the new parameter matrix into a multi-parameter viscoacoustic wave equation, and combining a deconvolution target function and a first-order gradient to perform iterative operation to obtain a next new parameter matrix; the iterative operation stop condition is that the synthesized sound pressure signal and the observed sound pressure signal are consistent.
In one embodiment, the ultrasonic CT data acquisition equipment is provided with a multilayer annular array, and the multilayer array acquires data simultaneously after single excitation to acquire three-dimensional volume data with a longitudinal aperture;
the data acquisition equipment is an observation array transducer containing a plurality of array elements; each array element of the plurality of array elements can transmit and receive ultrasonic signals;
optionally, the observation array transducer is selected from an annular probe, a through linear array probe, a through planar array probe and a through concave array probe.
In one embodiment, the observed acoustic pressure signal at the location of each of the plurality of array elements is d obs (x r ，t；x s ) Wherein x is s Is the coordinate position, x, of the transmitting array element r T represents the reception time for the coordinate position of the receiving array element.
In one embodiment, the multi-parameter visco-acoustic wave equation is a forward function, and the three-dimensional equation is as follows:
wherein the proton velocity describing the acoustic field is v i Stress is σ ij The synthesized sound pressure p ═ Σ k σ kk A/3; the parameters describing the physical properties of the medium are rho, sound velocity V p Acoustic impedance I p ＝ρV p Described, the attenuation coefficient is determined by the quality factor Q through the hysteresis coefficient
Introducing that the hysteresis elastic coefficient adopts L at a reference frequency
Memory variable of
If attenuation is not considered, the last two equations and
the coefficients of the correlations are truncated and degenerate back to the conventional isotropic acoustic wave equation, whose two-dimensional form can be formed by truncating the z-related term, the variable v i ，σ ij ，ρ，V p ，
Are both defined in the computed imaging space x and time tAn image;
given a transmit array element (position x) s ) The information of (2) is input into a parameter matrix of a sound velocity diagram, an acoustic impedance diagram and an acoustic attenuation diagram, sound pressure information p (x, t) on the whole calculation space x and time t can be solved discretely through a numerical simulation method (such as finite difference, finite element and the like), and then a receiving array element x is extracted r Synthesized sound pressure signal d cal (x r ，t；x s )＝R(x，x r )p(x，t；x s ) The maximum recording time T should be such that the transmitted wave signal through the target body tissue resulting from the farthest apart transmit and receive array element combination can be recorded.
In one embodiment, the deconvolution target function is an inversion function, so that the result of the FWI algorithm is closer to a global minimum, a time domain wiener filter between the measured data and the simulated data is solved, the target function is designed by weighting the filter, the transit time difference between the two data based on the wave equation can be effectively highlighted, the influence of inaccurate forward amplitude on inversion is weakened, and the deconvolution target function can be expressed as:
wherein A (tau) is a weighting function, w represents a wiener filter,
in inversion, the selection of the weight function A (tau) is crucial, and the weight function can determine the robustness of the algorithm, | tau |/tau |, and max 、τ 2 and gaussian functions may be chosen. The multi-scale inversion can be realized by setting an upper limit of the tau value in the weight function number, because the upper limit represents the maximum transit time difference between the observed data and the simulation data which can be considered by the objective function. The upper limit is selected by combining frequency division multi-scale inversion, and the effect of gradually reducing the upper limit of tau is realized because the frequency is gradually increased between 3 and 5 times of the reciprocal of the dominant frequency of the data of the corresponding frequency bandIf desired, the conventional two-norm objective function C is reused l2 The iteration increases the resolution of the image,
the traditional two-norm objective function still needs to be improved in the aspects of stability, accuracy, high efficiency and the like: under the conditions that low-frequency information is lacked in data and an initial model is far away from a real model, an optimization algorithm is difficult to effectively converge to a global minimum value to obtain accurate ultrasonic CT image reconstruction.
In one embodiment, the deconvolution objective function has a first order gradient with respect to the parameter matrix
It can be calculated by the adjoint state method,
wherein the adjoint wave field operator mu is an equation
The solution of (1). Changing the deconvolution objective function does not affect the calculation of the forward wavefield w. However, the adjoint wavefield μ is transmitted by reverse propagation
And obtaining the partial derivative, wherein the partial derivative changes with the change of the deconvolution objective function. In the case of a two-norm objective function, the partial derivative is:
under deconvolution conditions, the partial derivative is:
will be provided with
The adjoint wave field operator mu can be calculated by wave equation back propagation.
In one embodiment, the iteration may be viewed as a process of solving an optimization problem, with an initial parameter matrix m being input at iteration 0 0 I.e. estimation of acoustic parameter models based on a priori knowledge;
alternatively, in the case of detecting soft tissue such as breast, the initial parameter matrix m may be assumed 0 Equivalent to a homogeneous aqueous medium.
In one embodiment, the iterative update generally employs a gradient descent method, iteratively updating the values of the parameter matrix along a descent direction d,
m k+1 ＝m k +α k d k ，
wherein alpha is k Typically obtained using a line search method, and d k There are different calculation methods in different optimization strategies (e.g., steepest descent method, conjugate gradient method, L-BFGS method, Newton method), and the steepest descent method generally uses the negative direction of the objective function with respect to the gradient of the model as the descent direction, i.e., the steepest descent method
The conjugate gradient method, the L-BFGS method and other methods need to improve the current gradient by utilizing gradient information in the pre-iteration to obtain the descending direction, and the Newton method needs to estimate a Hatheri matrix, namely a second derivative of an objective function about the gradient, to restrict the current gradient to obtain the descending direction.
In one embodiment, the iteration stop condition may be a threshold value for decreasing the deconvolution objective function value or a number of iterations.
In one embodiment, the code for solving the multi-parameter visco-acoustic wave equation, deconvolution objective function can be validated by the CPU and then transplanted to the GPU to improve efficiency to meet the clinical diagnosis task requirements of time-sensitive ultrasound CT medical imaging.
In fig. 2, a single circular array is taken as an example to detect a mammary gland phantom by full matrix acquisition, and fig. 3(a) is a top view of a human tissue phantom in the present example, and the embodiment of the present application is specifically described under the condition of assuming water immersion observation. An observation system is selected as a single annular array, acoustic parameters of the mammary gland phantom are detected in a water coupling mode, the phantom contains 1.5-2.0mm of micro defects, the center frequency of a system emission pulse is 1.0MHz, and the specific flow is as follows:
step 1: after the ring array is used for surrounding the imitation, N is obtained by utilizing the acquisition mode of a full matrix s ×N e ×N t Radio frequency observation data d of size obs And (4) carrying out pretreatment. Wherein the data acquisition equipment contains N e Array transducer of individual array elements, N s For the number of shots, N t Is the number of time sampling points, and N s ≤N e 。
Step 2: two-dimensional acoustic wave equations based on sound velocity density parameters are selected, partial differential equation solutions of second-order time difference and fourth-order space staggered grid difference are adopted, the convolution perfect matching layer is combined to eliminate and calculate boundary echoes, the time and space grid sizes are determined through the library lambert, and observation data are resampled to guarantee the same time sampling rate. Then, the model is initialized by using the medium attribute of water, and the sound velocity in the initial model
Is uniform 1490m/s, density rho 0 Is 1000g/cm 3 I.e. by
Finally, estimating a wavelet time function s (t) at the excitation array element based on the initial model, then injecting the wavelet function at the excitation array element to generate a forward transmission wave field w, and then extracting simulation data d at a receiving position cal . Here, the compressible forward wavefield, or its boundary, may be recorded and the boundary may be used to reconstruct the forward wavefield in reverse.
And step 3: selecting an objective function, calculating partial derivatives of the objective function with respect to the simulation data
This partial derivative is also referred to as the residual. Calculating the adjoint wave field mu by the residual error through a backward propagation wave equation method, and on the corresponding time step, combining the backward transmission wave field mu with the recorded or reconstructed forward transmission wave field
Cross-correlation with zero lag is made and acted upon by a coefficient matrix. After the back propagation is finished, the gradient of the target function about the model parameter can be obtained
And 4, step 4: and (3) calculating the descending direction by using the gradient obtained in the step (3) through an optimization algorithm (such as a steepest descent method, a conjugate gradient method, LBFGS (length-weighted variance), a Newton method and the like), acting with a pre-condition operator matrix, and obtaining the updating step length by using a line search method. At this point, the existing parametric model, m, can be updated with the processed descent direction at k iterations k+1 ＝m k +α k d k . The loop is then iterated back to step 3 until a predetermined termination condition is reached.
And 5: the steps 2-4 are internal circulation, and in actual operation, observation data with different center frequencies can be formed by filtering the observation data. Then, the frequency is gradually increased from the low frequency, the processing of step 2-4 can be performed once in each frequency band, and the output model of the previous frequency band is adopted as the initial model in the next frequency band. Different frequency bands may also employ different objective functions.
Step 6: and stopping iteration when the termination condition is reached.
And 7: and outputting a multi-parameter image, wherein fig. 3(b) is an ultrasonic CT acoustic impedance image based on FWI in the present embodiment, and fig. 3(c) is an ultrasonic CT sound velocity image based on FWI in the present embodiment.
The invention provides an ultrasonic CT imaging based on a wave equation, in particular to an ultrasonic CT multi-parameter image reconstruction system and method based on FWI. The numerical solution of the multi-parameter wave equation can be developed on a CPU and deployed on a GPU for actual detection tasks. And constructing an innovative multi-scale deconvolution target function and proposing a specific use strategy of the corresponding target function in actual reconstruction. The robustness of model reconstruction is improved by emphasizing the idea of comparing the dynamic characteristics of simulation data and observation data.
It should be noted that the methods provided herein are not inherently related to any particular computer, virtual machine, or other apparatus. Various general purpose devices may be used with the teachings herein. The required structure for constructing such a device will be apparent from the description above. In addition, the present invention is not directed to any particular programming language. It is to be appreciated that a variety of programming languages may be used to implement the teachings of the present invention as described herein, and that the specific languages, calls for system function blocks, are provided for disclosure as preferred embodiments of the present invention.
In the description provided herein, numerous specific details are set forth. It is understood, however, that embodiments of the invention may be practiced without these specific details. In some instances, well-known methods, structures and techniques have not been shown in detail in order not to obscure an understanding of this description.
It will be apparent to those skilled in the art that various changes and modifications may be made in the present invention without departing from the spirit and scope of the invention. Thus, it is intended that the present invention cover the modifications and variations of this invention provided they come within the scope of the appended claims and their equivalents.
Claims (10)
1. A medical ultrasound CT multi-parameter image reconstruction method based on a wave equation is characterized by comprising the following steps:
transmitting and receiving ultrasound by using ultrasound CT data acquisition equipment to acquire real three-dimensional space data; the data is an observed sound pressure signal of the position of the receiving array element in the real three-dimensional space;
simulating the process of transmitting and receiving ultrasound by the acquisition equipment by using a multi-parameter viscous sound wave kinetic equation to obtain a synthetic sound pressure signal corresponding to the position of the receiving array element in the real three-dimensional space; the input variable of the multi-parameter viscous sound wave equation is a parameter matrix composed of an original sound velocity graph, an original acoustic impedance graph and an original attenuation quality factor graph;
solving a difference between the observed sound pressure signal and the synthesized sound pressure signal by using a deconvolution target function and a two-norm target function;
calculating the deconvolution target function and the two-norm target function by a adjoint state method to obtain a first-order gradient of the parameter matrix; calculating to obtain a new parameter matrix based on the first-order gradient and the parameter matrix; converting the new parameter matrix into a new sound velocity map, a new acoustic impedance map and a new attenuation quality factor map;
inputting the new parameter matrix into the multi-parameter viscoacoustic wave equation, and combining the deconvolution target function and the first-order gradient to perform iterative operation to obtain a next new parameter matrix; the iterative operation stop condition is that the synthesized sound pressure signal and the observed sound pressure signal are consistent.
2. The multi-parameter image reconstruction method according to claim 1, wherein the ultrasound CT data acquisition equipment is provided with a multi-layer annular array, and the multi-layer array acquires data simultaneously after a single excitation to acquire three-dimensional volume data with a longitudinal aperture;
the data acquisition equipment is an observation array transducer containing a plurality of array elements; each array element of the plurality of array elements can transmit and receive ultrasonic signals;
optionally, the observation array transducer is selected from an annular probe, a through linear array probe, a through planar array probe and a through concave array probe.
3. Method for multi-parameter image reconstruction according to claim 2, characterized in thatWherein the observed sound pressure signal at the position of each of the plurality of array elements is d obs (x r ，t；x s ) Wherein x is s To the coordinate position of the transmitting array element, x r T represents the reception time for the coordinate position of the receiving array element.
4. The multi-parameter image reconstruction method according to claim 1, wherein the multi-parameter visco-acoustic wave equation is a forward function, and the formula of the three-dimensional equation is as follows:
wherein the proton velocity describing the acoustic field is v i Stress of σ ij The synthesized sound pressure p ═ Σ k σ kk A/3; the parameters describing the physical properties of the medium are rho, sound velocity V p Acoustic impedance I p ＝ρV p Described, the attenuation coefficient is determined by the quality factor Q through the hysteresis coefficient
Introducing that the hysteresis elastic coefficient adopts L at a reference frequency
Memory variable of
If attenuation is not considered, the last two equations and
the coefficients of the correlations are truncated and degenerate back to the conventional isotropic acoustic wave equation, whose two-dimensional form can be formed by truncating the z-related term, the variable v i ，σ ij ，ρ，V p ，
Are both images defined in the computed imaging space x and time t;
given a transmit array element (position x) s ) The information of (2) is input into a parameter matrix of a sound velocity diagram, an acoustic impedance diagram and an acoustic attenuation diagram, sound pressure information p (x, t) on the whole calculation space x and time t can be solved discretely through a numerical simulation method (such as finite difference, finite element and the like), and then a receiving array element x is extracted r Synthesized sound pressure signal d cal (x r ，t；x s )＝R(x，x r )p(x，t；x s ) The maximum recording time T should be such that the transmitted wave signal through the target body tissue resulting from the farthest apart transmit and receive array element combination can be recorded.
5. The multi-parameter image reconstruction method according to claim 1, wherein the deconvolution target function is an inversion function, so that the result of the FWI algorithm is closer to a global minimum, a time domain wiener filter between the measured data and the simulated data is solved, the filter is weighted to design the target function, so that the transit time difference between the two data based on the wave equation can be effectively highlighted, the influence of inaccurate forward amplitude on inversion is weakened, and the deconvolution target function can be expressed as:
wherein A (tau) is a weighting function, w represents a wiener filter,
in inversion, the selection of the weight function A (tau) is crucial, and the weight function can determine the robustness of the algorithm, | tau |/tau |, and max 、τ 2 and gaussian functions may be chosen. The characteristics of multi-scale inversion can be realized by setting the upper limit of the tau value in the weight function number, because the upper limit represents the target functionThe maximum time difference of flight for observed and simulated data can be considered. The upper limit is selected by combining frequency division multi-scale inversion, the effect of gradually reducing the upper limit of tau is realized between 3 and 5 times of the reciprocal of the dominant frequency of corresponding frequency band data due to the gradual increase of the frequency, and then the traditional two-norm objective function C is utilized l2 The iteration increases the resolution of the image,
the conventional two-norm objective function still needs to be improved in the aspects of stability, accuracy, high efficiency and the like: under the conditions that low-frequency information is lacked in data and an initial model is far away from a real model, an optimization algorithm is difficult to effectively converge to a global minimum value to obtain accurate ultrasonic CT image reconstruction.
6. Method for multi-parametric image reconstruction according to claim 1, characterized in that the deconvolution objective function is the first order gradient with respect to the parameter matrix
It can be calculated by the adjoint state method,
wherein the adjoint wave field operator mu is an equation
The solution of (1). Changing the deconvolution objective function does not affect the calculation of the forward wavefield w. However, the adjoint wavefield μ is transmitted by reverse propagation
And obtaining the partial derivative, wherein the partial derivative changes with the change of the deconvolution objective function. In the case of a two-norm objective function, the partial derivative is:
under deconvolution conditions, the partial derivative is:
will be provided with
The adjoint wave field operator mu can be calculated by wave equation back propagation.
7. The method of claim 1, wherein the iteration is regarded as a process of solving an optimization problem, and the initial parameter matrix m is input at the 0 th iteration 0 I.e. estimation of acoustic parameter models based on a priori knowledge;
alternatively, in the case of detection of soft tissue such as breast, it may be assumed that the parameter matrix m is initiated 0 Equivalent to a homogeneous aqueous medium.
8. The multi-parameter image reconstruction method of claim 1, wherein said iterative updating iteratively updates values of said parameter matrix along a descent direction d, typically using a gradient descent method,
m k+1 ＝m k +α k d k ，
wherein alpha is k Typically obtained using a line search method, and d k There are different calculation methods in different optimization strategies (e.g., steepest descent method, conjugate gradient method, L-BFGS method, Newton method), and the steepest descent method generally uses the negative direction of the objective function with respect to the gradient of the model as the descent direction, i.e., the steepest descent method
The conjugate gradient method, the L-BFGS method and other methods need to improve the current gradient by utilizing gradient information in the pre-iteration to obtain the descending direction, and the Newton method needs to estimate a Hatheri matrix, namely a second derivative of an objective function about the gradient, to restrict the current gradient to obtain the descending direction.
9. The multi-parameter image reconstruction method according to claim 1, wherein the iterative operation stop condition is selected from a threshold value of a decrease in the deconvolution objective function value or a number of iterations.
10. The multi-parameter image reconstruction method according to claim 1, wherein the code for solving the multi-parameter visco-acoustic wave equation and the deconvolution objective function can be verified by a CPU and then transplanted to a GPU to improve efficiency for meeting the clinical diagnosis task requirement of time-sensitive ultrasonic CT medical imaging.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210604492.XA CN114972567B (en) | 2022-05-30 | 2022-05-30 | Medical ultrasonic CT multi-parameter image reconstruction method based on wave equation |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210604492.XA CN114972567B (en) | 2022-05-30 | 2022-05-30 | Medical ultrasonic CT multi-parameter image reconstruction method based on wave equation |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114972567A true CN114972567A (en) | 2022-08-30 |
| CN114972567B CN114972567B (en) | 2023-02-03 |
Family
ID=82957516
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210604492.XA Active CN114972567B (en) | 2022-05-30 | 2022-05-30 | Medical ultrasonic CT multi-parameter image reconstruction method based on wave equation |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114972567B (en) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116327250A (en) * | 2023-02-13 | 2023-06-27 | 中国科学院地质与地球物理研究所 | Mammary gland ultrasonic three-dimensional imaging method based on full waveform inversion technology |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106037795A (en) * | 2016-05-17 | 2016-10-26 | 南京大学 | Method for reconstructing sound velocity image by virtue of optimized full wave inversion method |
| CN109584323A (en) * | 2018-10-18 | 2019-04-05 | 天津大学 | The information constrained coeliac disease electrical impedance images method for reconstructing of ultrasonic reflection |
| CN113689545A (en) * | 2021-08-02 | 2021-11-23 | 华东师范大学 | 2D-to-3D end-to-end ultrasonic or CT medical image cross-modal reconstruction method |
-
2022
- 2022-05-30 CN CN202210604492.XA patent/CN114972567B/en active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106037795A (en) * | 2016-05-17 | 2016-10-26 | 南京大学 | Method for reconstructing sound velocity image by virtue of optimized full wave inversion method |
| CN109584323A (en) * | 2018-10-18 | 2019-04-05 | 天津大学 | The information constrained coeliac disease electrical impedance images method for reconstructing of ultrasonic reflection |
| CN113689545A (en) * | 2021-08-02 | 2021-11-23 | 华东师范大学 | 2D-to-3D end-to-end ultrasonic or CT medical image cross-modal reconstruction method |
Non-Patent Citations (5)
| Title |
|---|
| GANG YAO等: "A review on reflection-waveform inversion", 《PETROLEUM SCIENCE》 * |
| 余锦华等: "医学超声成像的模拟研究", 《声学技术》 * |
| 刘玉柱等: "基于Born敏感核函数的VTI介质多参数全波形反演", 《地球物理学报》 * |
| 张宇: "从成像到反演:叠前深度偏移的理论、实践与发展", 《石油物探》 * |
| 李庆洋等: "基于Student′s t分布的不依赖子波最小二乘逆时偏移", 《地球物理学报》 * |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116327250A (en) * | 2023-02-13 | 2023-06-27 | 中国科学院地质与地球物理研究所 | Mammary gland ultrasonic three-dimensional imaging method based on full waveform inversion technology |
| CN116327250B (en) * | 2023-02-13 | 2023-08-25 | 中国科学院地质与地球物理研究所 | Mammary gland ultrasonic three-dimensional imaging method based on full waveform inversion technology |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114972567B (en) | 2023-02-03 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110381845B (en) | Ultrasound imaging system with neural network for deriving imaging data and tissue information | |
| US20190336107A1 (en) | Ultrasound imaging system with a neural network for image formation and tissue characterization | |
| US10743837B2 (en) | Ultrasound waveform tomography method and system | |
| US20060184021A1 (en) | Method of improving the quality of a three-dimensional ultrasound doppler image | |
| KR20010014491A (en) | Semi-automated segmentation method for 3-dimensional ultrasound | |
| AU2017375899B2 (en) | Method of, and apparatus for, non-invasive medical imaging using waveform inversion | |
| JP2022503946A (en) | Image reconstruction method based on learning nonlinear mapping | |
| CN110772281B (en) | Ultrasonic CT imaging system based on improved ray tracing method | |
| CN106618635A (en) | Shear wave elastic imaging method and device | |
| JP2005514997A (en) | Observation system having means for processing an ultrasound image sequence for performing a quantitative calculation of flow in a body organ | |
| CN105120761A (en) | Ultrasound vector flow imaging (VFI) with curve tracing | |
| CN114972567B (en) | Medical ultrasonic CT multi-parameter image reconstruction method based on wave equation | |
| Kretzek et al. | GPU-based 3D SAFT reconstruction including attenuation correction | |
| JP2023521817A (en) | System and method for non-invasive pressure measurement | |
| KR20180031559A (en) | Apparatus and method for displaying ultrasound image | |
| Benjamin et al. | Reflection ultrasound tomography using localized freehand scans | |
| Chumchob et al. | An improved variational model and its numerical solutions for speckle noise removal from real ultrasound images | |
| Long et al. | Deep Learning Ultrasound Computed Tomography under Sparse Sampling | |
| KR102388950B1 (en) | Method for measuring shear wave speed and ultrasonic device for performing the same | |
| US20210177380A1 (en) | Method and apparatus for quantitative ultrasound imaging using single-ultrasound probe | |
| US20230228873A1 (en) | Systems and methods for generating color doppler images from short and undersampled ensembles | |
| Seliverstova et al. | Cardiac shear wave speed estimation in 3D: an in silico and in vivo study | |
| KR20240015562A (en) | Method and system of linking ultrasound image data associated with a medium with other image data associated with the medium | |
| Mendizabal Ruiz | LUMEN SEGMENTATION IN INTRAVASCULAR ULTRASOUND DATA | |
| Nguyen | Information theoretic design of breast sonography |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |