Multi-thread fusion in vivo ultrasonic thermal strain imaging method and device

Abstract

The present invention discloses a multi-thread fused in vivo ultrasonic thermal strain imaging method and device, which belongs to the field of medical image processing. The method of the present invention is to collect an ultrasonic image sequence, divide the image into multiple threads according to the correlation, and perform image registration and thermal strain imaging. Furthermore, the biological heat conduction equation solved by the time-domain finite difference method is used as a prediction model, and an adaptive Kalman filter is used to filter the thermal strain sequence. Finally, the thermal strains of each thread are fused according to the scalar weighted linear minimum variance information fusion criterion. The present invention overcomes the shortcomings of existing in vivo thermal strain imaging that the imaging results rely on a single motion phase, physiological motion is difficult to accurately compensate, resulting in thermal strain contour distortion and poor signal-to-noise ratio, and conventional spatial filtering sacrifices spatial resolution. Accurate, robust, and high spatiotemporal resolution in vivo thermal strain imaging results can be obtained.

Description

Technical Field

[0001] The present invention relates to the field of medical image processing, and more specifically to a multi-threaded fusion in vivo ultrasonic thermal strain imaging method and device. Background Art

[0002] Ultrasonic thermal strain imaging is a method that analyzes the movement of tissue scattering points in the image during the heating process, and then obtains the thermal displacement and strain of the tissue, and finally achieves quantitative estimation of temperature. However, in in vivo thermal strain imaging, the in vivo motion between images will cause large-amplitude strain artifacts, and even cause thermal strain imaging failure. In order to overcome the influence of in vivo motion, existing thermal strain imaging methods usually use image registration or spatial filtering methods. However, due to the lack of motion phase analysis, the motion suppression effect of image registration is easily affected by the motion state of the starting image, and the spatial filtering method often sacrifices spatial resolution.

[0003] In the prior art, some technical solutions have been proposed for overcoming the thermal strain imaging of living body motion, such as the invention titled: Ultrasonic method for measuring temperature changes of biological tissues based on thermal expansion and gating algorithm (application date: September 25, 2017; application number: 201710876349), which discloses an ultrasonic method for measuring temperature changes of biological tissues based on thermal expansion and gating algorithm, collecting B-ultrasound image sequences during the heating process of the living body, selecting target frames, calculating the time delay image when the ultrasound passes through the tissue and thereby obtaining a temperature change image; training an adaptive filter based on the image outside the heating area, and performing noise suppression on the obtained temperature change image. The disadvantage of this method is that its motion suppression effect is easily affected by the motion phase of the starting image, has poor robustness, and has low output time density. The invention is titled: An adaptive motion compensation method and system for ultrasonic thermal imaging technology (application date: March 24, 2020; application number: 202010213806). The scheme discloses an adaptive motion compensation method and system for ultrasonic thermal imaging technology, which divides the target area of ​​ultrasonic thermal imaging into blocks, selects local reference points for motion compensation, trains the filter coefficients for motion compensation based on the reference points, and performs motion compensation on the temperature distribution. The shortcoming of this method is that it ignores the local differences in living body motion, and the obtained motion compensation filter cannot accurately filter out motion artifacts in the target area.

[0004] In summary, how to obtain thermal strain images with high signal-to-noise ratio, time resolution, and spatial resolution without being affected by the motion state of the starting image, is an urgent problem to be solved in the existing technology. Summary of the invention

[0005] 1. Technical issues to be solved

[0006] In view of the problems in the prior art that the accuracy of thermal strain imaging depends on a single motion phase and physiological motion is difficult to be fully suppressed, resulting in thermal strain contour distortion and poor signal-to-noise ratio, while conventional spatial filtering sacrifices spatial resolution, the present invention provides a multi-threaded fusion in vivo ultrasonic thermal strain imaging method and device, which can realize thermal strain imaging with multi-threaded image registration, and finally fuse stable, accurate, and high-temporal and spatial resolution in vivo thermal strain images through adaptive Kalman filtering with heat conduction equation as the prediction model.

[0007] 2. Technical solution

[0008] The purpose of the present invention is achieved through the following technical solutions.

[0009] The present invention provides a multi-thread fusion in vivo ultrasonic thermal strain imaging method, comprising the following steps:

[0010] Acquisition steps: using a thermal therapy method to heat the target area, acquiring ultrasound images of the live thermal therapy target area, and obtaining a digital ultrasound image sequence;

[0011] Imaging step: The digital ultrasound image sequence obtained in the acquisition step is divided into multiple quasi-periodic threads n (n = 1, 2, ..., N) corresponding to different motion phases through correlation analysis, where N is a positive integer representing the total number of threads; a region of interest is selected in the image according to the target area position, and image registration and thermal strain imaging are performed independently in each thread based on the region of interest to obtain the thermal strain distribution sequence S of each thread n ;

[0012] Calibration steps: S n The heat source center position of the corresponding thread is obtained by fitting the two-dimensional axisymmetric surface function. By aligning the heat source centers of each thread, S n Updated to the thermal strain distribution sequence after position calibration;

[0013] Modeling steps: Solve the biological heat conduction equation by using the finite difference time domain method and construct the prediction model of the Kalman filter;

[0014] Initialization step: set the heat source intensity distribution and initialize the adaptive Kalman filter parameters;

[0015] Filtering step: Take the thermal strain sequence of each thread as the measured value of the Kalman filter, and take the calculation result of the biological heat conduction equation as the predicted value of the Kalman filter to obtain the thermal strain estimate S' after adaptive Kalman filtering n ;

[0016] Fusion step: Based on the system error after thermal strain filtering of each thread, the thermal strain results of each thread are fused according to the scalar weighted linear minimum variance information fusion criterion to output the final thermal strain image sequence S'.

[0017] Furthermore, S n The specific implementation method for fitting the two-dimensional axisymmetric surface function is: take the natural logarithm of the thermal strain of each thread lnS n , polynomial least squares is used to fit it according to the following formula:

[0018] f(x,y)=ax 2 +by 2 +cx+dy+e, (1)

[0019] Get the fitting parameters a, b, c, d, e; then the heat source center position O n The vertical and horizontal coordinates are determined as

[0020]

[0021] Furthermore, the specific implementation of the modeling step is:

[0022] For any thread n, with the heat source as the center, a discrete grid of size A×L is established, where A is the axial size and L is the lateral size, which is not less than the region of interest, the maximum spatial step Δd, the time span 0~T, and the time step Δt. The value of Δd is less than the smaller value of A / 2 and L / 2, and the value of Δt is less than T / 2. The prediction model is

[0023]

[0024] in, and S n ′(d,t m-1 ) are the predicted value of thermal strain at time tm=mΔt and t m-1 = (m-1)Δt Kalman filter thermal strain estimate; d is the spatial coordinate, which can be in the form of standard Cartesian coordinate system, polar coordinate system (two-dimensional), cylindrical coordinate system or spherical coordinate system (three-dimensional); if the thermal strain of each thread is not defined in the coordinate system d, it is first transformed to d by the standard method to obtain S n (d; t); F = A -1 b is the state transfer matrix of the system, A -1 is the inverse matrix of matrix A, B(d,t m )=A -1 Q(d;t m ) / k is the system input at time tm; A and b are coefficient matrices related to tissue specific heat capacity, thermal conductivity, and density, which are obtained by discretizing the biological heat conduction equation through the time-domain finite difference method; k is a constant specified in the thermal strain temperature measurement model δ-ESF, which describes the linear relationship between temperature change and thermal strain; Q(d; t m) is the heat source term, which can be time-varying or time-invariant.

[0025] Furthermore, in the initialization step, the specific implementation method of setting the heat source distribution is:

[0026] For any thread n, based on the following discrete equation

[0027] AS n (d;t m+1 )=bS n (d;t m )+Q(d;t m ) / k, (4)

[0028] Get t m The heat source term at the moment; if the heat source is time-invariant, the heat source terms Q(d; t m ) to get the heat source term of the current thread; average the heat source terms of all threads to get the heat source term used for the prediction model; if the heat source is time-varying, then at each moment t m , the heat source item Q(d; t m ) average to obtain the heat source term used in the prediction model;

[0029] Furthermore, in the initialization step, the specific implementation method of initializing the adaptive Kalman filter parameters is:

[0030] The initial thermal strain S0(d), the covariance W(d) of the estimation error of the prediction model, the covariance R(d) of the measurement error, and the covariance P(d) of the filter estimation error are set according to the initial temperature and thermal strain accumulation of the tissue. Among them, the starting values ​​W0(d), R0(d), and P0(d) should be greater than the maximum noise amplitude.

[0031] Furthermore, the specific implementation of the filtering step is:

[0032] At t0=0, the thermal strain estimate of the Kalman filter is set to S' n (d, t0) = S0(d); for each thread, independently calculate t based on equation (3) m-1 The estimated value of thermal strain at time t is obtained m The predicted value of thermal strain at time t m-1 The covariance P of the system estimation error at the moment n,m-1 , predict t according to the following formula m Covariance of the system estimation error at time

[0033]

[0034] Among them, F Trepresents the transposed matrix of the state transfer matrix F, W n,m-1 is the covariance matrix of the current estimation error of the prediction model; further, And the measurement value estimation error R = R0(d), the Kalman gain K is calculated as follows n :

[0035]

[0036] Furthermore, K n Acting on the predicted value and the measured value S n (d;t m ), the filtered t is calculated as follows m The estimated thermal strain S' at time n (d;t m ):

[0037]

[0038] in, is the new information, representing the deviation between the measured value and the predicted value; further, the covariance matrix P of the system estimation error is updated according to the following formula n,m :

[0039]

[0040] Furthermore, through the latest M iterations of new information v n,m , calculate the new information sequence C as follows m :

[0041]

[0042] Among them, the value of M is an integer greater than 1; further, the covariance matrix W of the current prediction model estimation error is iteratively updated according to the following formula: n,m :

[0043]

[0044] Adaptively correct the weights of model predictions and measurements in the final filter output;

[0045] For any thread n, according to the iterative calculation process of equations (4)-(10), make t m From 0 to T, the final estimated image sequence S' n That is the thermal strain output image sequence obtained by this thread after adaptive Kalman filtering.

[0046] Furthermore, the specific implementation of the fusion step is:

[0047] For any thread n, at any time tm , the system estimation error covariance matrix P updated at the current time n,m , calculate the weight A of the current thread at the current moment according to the scalar weighted linear minimum variance fusion criterion shown in the following formula n,m :

[0048]

[0049] Among them, tr(P n,m ) represents the matrix P n,m Further, the weight A of each thread n,m , the thermal strain distribution of all threads at the current moment after fusion filtering is as follows:

[0050]

[0051] And the covariance of the system estimation error is updated as:

[0052]

[0053] In formula (12), S'(d,t m ) is t m The thermal strain output image after all threads are fused at the moment; by making t m From 0 to T, according to the iteration process of equations (11)-(13), the final output thermal strain image sequence S' can be obtained after adaptive Kalman filtering and thread fusion.

[0054] The present invention also provides a multi-threaded fusion in vivo ultrasonic thermal strain imaging device, comprising: a computer, a B-ultrasound instrument, a linear array probe and a thermal ablation device;

[0055] The B-ultrasound instrument is connected to the linear array probe and is used to collect ultrasound images of living objects;

[0056] Thermal ablation devices are used to heat living subjects, creating displacement and strain distributions in tissue;

[0057] The computer is connected to the B-ultrasound instrument and the thermal ablation device respectively, and is used to control the B-ultrasound instrument and the thermal ablation device, and is configured to implement the in vivo ultrasonic thermal strain imaging method with multi-thread fusion as described above.

[0058] 3. Beneficial effects

[0059] Compared with the prior art, the advantages of the present invention are:

[0060] The present invention adopts an adaptive Kalman filter that combines the biological heat conduction equation with ultrasonic image thermal strain imaging, so that the artifacts caused by physiological movement are corrected by the heat conduction model; multi-threaded thermal strain imaging is independent of each other, providing more complete measurement information, and the failure of any thread does not affect the work of other threads, which greatly enhances the reliability and robustness of the algorithm; thermal strain fusion based on the scalar weighted linear minimum variance criterion improves the temporal resolution, accuracy and signal-to-noise ratio of thermal strain without losing spatial resolution. BRIEF DESCRIPTION OF THE DRAWINGS

[0061] Figure 1 A schematic diagram of the process of a multi-threaded fusion in vivo ultrasound thermal strain imaging method;

[0062] Figure 2 Schematic diagram of a device for heating visceral fat in a living body using a microwave ablation needle and collecting B-ultrasound images of the heating process in Example 2;

[0063] Figure 3 The thermal strain distribution diagrams of the first and eighth threads at different times obtained through image registration and thermal strain calculation in Example 2;

[0064] Figure 4 Schematic diagram of the process of adaptive Kalman filtering in Example 2;

[0065] Figure 5 The first thread and the eighth thread obtained after the adaptive Kalman filter in Example 2 are Figure 3 Thermal strain distribution diagram at the corresponding time;

[0066] Figure 6 This is the thermal strain distribution diagram at t=20 seconds after the multi-thread thermal strain fusion in Example 2 and the thermal strain-time curve diagram of the heat source position.

[0067] Figure numerals: 01, computer; 02, ultrasound imaging system; 03, imaging probe; 04, microwave heater; 05, microwave ablation needle. DETAILED DESCRIPTION

[0068] In order to make the purpose, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are part of the embodiments of the present invention, not all of the embodiments; moreover, the various embodiments are not relatively independent and can be combined with each other as needed to achieve better results. Therefore, the following detailed description of the embodiments of the present invention provided in the drawings is not intended to limit the scope of the invention claimed for protection, but merely represents selected embodiments of the present invention. Based on the embodiments in the present invention, all other embodiments obtained by ordinary technicians in this field without making creative work are within the scope of protection of the present invention.

[0069] In order to further understand the content of the present invention, the present invention is described in detail in conjunction with the accompanying drawings and embodiments.

[0070] Example 1

[0071] Combination Figure 1 As shown, the multi-thread fusion in vivo ultrasonic thermal strain imaging method of the present invention comprises the following steps:

[0072] Acquisition steps: The target area is heated by a thermal therapy method, and ultrasonic images of the live thermal therapy target area are acquired to obtain a digital ultrasonic image sequence.

[0073] Collect digital ultrasound image sequences of target tissues in living bodies. Specifically, the B-ultrasound instrument works in real-time imaging mode, and the probe angle and imaging system configuration are adjusted to make the target area in the imaging field of view. It is worth noting that, according to the depth from the body surface to the target tissue and the acoustic path conditions in the embodiment, imaging probes of different types, different numbers of array elements, and different geometric structures can be selected, including linear arrays and convex array probes with 128 / 256 array elements, and phased array probes with 80 array elements. At the same time, different imaging modes can be selected, such as fundamental wave imaging, harmonic imaging, and synthetic aperture imaging.

[0074] The target area is heated using thermal therapy methods such as focused ultrasound and radiofrequency ablation. It is worth noting that the heating method can be microwave ablation, radiofrequency ablation, laser ablation, ultrasonic ablation, infrared ablation, etc. During the heating process, when it is difficult for the living body to maintain a stable physiological movement state, the living body can be anesthetized by intramuscular injection, intravenous injection, inhalation of anesthetic gas, etc. There is no limit to the breathing method of the living body, and it can be spontaneous breathing, ventilator or artificially supported breathing, etc.

[0075] The images are collected at equal time intervals, and the collection time T is not less than two respiratory cycles, to obtain a digital ultrasound image sequence. It is worth noting that the present invention supports various types of image signal types, including channel signals, ultrasound radio frequency signals after beamforming, analytical signals after orthogonal demodulation or Hilbert transformation, envelope signals obtained by taking the amplitude of analytical signals, and signals after compression and filtering based on the above signals.

[0076] Imaging step: The digital ultrasound image sequence obtained in the acquisition step is divided into multiple quasi-periodic threads n (n = 1, 2, ..., N) corresponding to different motion phases through correlation analysis, where N is a positive integer representing the total number of threads; a region of interest is selected in the image according to the target area position, and image registration and thermal strain imaging are performed independently in each thread based on the region of interest to obtain the thermal strain distribution sequence S of each thread n .

[0077] Calculate the correlation coefficient between images to determine the respiratory frequency and exhalation state feature frames; specifically, take the M frames of the ultrasound image sequence containing at least one respiratory cycle as the template frame X, take all the acquired N frames of the image as the detection frame Y, and calculate the normalized cross-correlation to obtain the cross-correlation coefficient-time curve γ m (n),

[0078]

[0079] Where n=1,2,…,N,m=1,2,…,M. COV(·) represents covariance, σ X and σ X is the standard deviation of X and Y. m (n) Through linear amplitude modulation z transformation, the corresponding spectrum curve Γ is obtained m (f), f is the frequency variable, and the maximum amplitude of each spectrum curve is recorded as A m . A m The frequency corresponding to the maximum value is the respiratory frequency f res , the template frame number m=p corresponding to the curve, then the pth frame is the exhalation state feature frame.

[0080] According to the cross-correlation coefficient-time curve γ corresponding to the exhalation feature frame p p (n), divide the image sequence into different cycles of exhalation and inhalation states and multiple motion phase periods; set the cutoff frequency f of the finite impulse response low-pass filter C =f res , the filter order is n tap The filter is applied to the cross-correlation coefficient-time curve γ of the exhalation state feature frame as the template frame p (n), and the high-frequency component is filtered out to obtain the curve γ' p(n). Alternately search for the curve γ' along the time direction p The maximum and minimum of (n), every two adjacent maximum values ​​constitute the boundary of an inhalation cycle, and every two adjacent minimum values ​​constitute the boundary of an exhalation cycle. Each exhalation and inhalation cycle is divided into N1 and N2 phase periods at equal time intervals, and a total of N different motion phase periods are obtained.

[0081] In the same phase period of different cycles, the image registration based on the region of interest is obtained to obtain the registered multi-threaded image sequence Q n (n=1,2,…,N); select the image at the middle moment of the first cycle of each phase period, and further select the region of interest with the center point coordinates (i0,j0) and the size A0×L0 as the reference frame Q according to the target area position ref , put it into the subset Q n In the same phase of the next cycle, for each frame of the image, the area with a size of A0×L0 and a center point at (i0+δa,j0+δl) is taken as the detection frame Q tar . Change δa and δl in the range of -A0 / 3~A0 / 3 and -L0 / 3~L0 / 3 respectively, where δa and δl are the axial and lateral displacements respectively, so that Q tar With Q ref The two-dimensional correlation coefficient of the two-dimensional correlation coefficient reaches the maximum value, and the corresponding δa and δl are Q tar The displacements Δa and Δl. tar After axial and lateral shifts of Δa and Δl respectively, the registered image Q' is obtained. tar . Compare all the maximum mutual correlation coefficients and take the Q' corresponding to the maximum value tar As the new Q ref , put into subset Q n At the end of the Q ref The same operation is performed for the next cycle, and the cycle is repeated until all cycles are registered, and all phase periods are traversed to obtain the multi-threaded image subset Q n .

[0082] Multithreaded Image Subsetting Q n The thermal strain is calculated in parallel, and a denoising function is constructed based on the correlation between threads to denoise the thermal strain and obtain the thermal strain image sequence S n ; Each thread image subset Q n After obtaining the displacement image sequence based on the well-known block matching algorithm, the thermal strain image sequence S is obtained by axial difference. n . Take image S n The spatial point (i, j) is taken as the center and a rectangular window R with a size of A1×L1 is taken. n , and compare it with the strain S of other threads vThe corresponding part R of (v=1,2,…,N,v≠n) v The average cross-correlation coefficient was calculated as follows:

[0083]

[0084] Among them, COV(R n ,R v ) represents R n With R v The covariance of R1 and σ Rv For R n and R v The standard deviation of . Denoising function q n (i,j) is defined as:

[0085]

[0086] Among them, α is the attenuation coefficient, and the value range is 10 -2 ~10 2 , p th is the correlation coefficient threshold, which can range from 0.4 to 0.9. n (i,j) and S n Multiply, then S n Updated to the denoised thermal strain image.

[0087] Calibration steps: S n The heat source center position of the corresponding thread is obtained by fitting the two-dimensional axisymmetric surface function. By aligning the heat source centers of each thread, S n Updated to the position-calibrated thermal strain distribution sequence.

[0088] Based on the two-dimensional surface fitting method, the heat source center coordinates of each thread are determined. n is the center coordinate of the heat source of thread n, and based on this, the spatial position of the thermal strain image of each thread is corrected to obtain the thermal strain distribution sequence after all threads are spatially aligned, S n is the spatially aligned thermal strain distribution sequence of thread n, with the heat source position O n For invasive heating, the coordinates of the heating needle can be extracted from the image, and for non-invasive heating, the spatial position of the preset focus can be obtained; S n Fit the two-dimensional axisymmetric surface function to obtain the heat source center position O of the corresponding thread n , by aligning the heat source centers of each thread, S n Update to the thermal strain sequence after position calibration. Specifically, take the natural logarithm of the thermal strain of each thread lnS n , the well-known polynomial least squares method is used to fit it according to the following formula:

[0089] f(x,y)=ax 2 +by 2 +cx+dy+e, (4)

[0090] Get the fitting parameters a, b, c, d, e; then the heat source center position O n The vertical and horizontal coordinates are determined as

[0091]

[0092] It is worth noting that the surface fitting function includes but is not limited to a two-dimensional cone function, a two-dimensional raised cosine function, a two-dimensional Hamming window function, and a two-dimensional Blackman window function; another feasible spatial correction method uses the ultrasound image corresponding to a certain thread as a template frame, and the ultrasound images of other threads as detection frames, and adopts a two-dimensional cross-correlation based on a sliding window algorithm. Based on the calculated displacement field, the thermal strain distribution of each thread is shifted point by point in space.

[0093] Modeling steps: Solve the biological heat conduction equation through the time-domain finite difference method and build the prediction model of the Kalman filter.

[0094] For any thread n, with the heat source as the center, a discrete grid with size A (axial) × L (lateral) not less than the area of ​​interest, maximum spatial step Δd, time span 0~T, and time step Δt is established. The value of Δd is less than the smaller value of A / 2 and L / 2, and the value of Δt is less than T / 2. The prediction model is

[0095]

[0096] in, and S n ′(d,t m-1 ) are t m = predicted value of thermal strain at time mΔt and t m-1 = (m-1)Δt Kalman filter thermal strain estimate; d is the spatial coordinate, which can be in the form of standard Cartesian coordinate system, polar coordinate system (two-dimensional), cylindrical coordinate system or spherical coordinate system (three-dimensional); if the thermal strain of each thread is not defined in the coordinate system d, it is first transformed to d by the standard method to obtain S n (d; t); F = A -1 b is the state transfer matrix of the system, A -1 is the inverse matrix of matrix A, B(d,t m )=A -1 Q(d;t m ) / k is t mA and b are coefficient matrices related to tissue specific heat capacity, thermal conductivity, and density, which can be obtained by discretizing the biological heat conduction equation using the time-domain finite difference method by a known method; k is a constant specified in the known thermal strain temperature measurement model δ-ESF, which describes the linear relationship between temperature change and thermal strain; Q(d; t m ) is the heat source term, which can be time-varying or time-invariant.

[0097] It is worth noting that the present invention adopts the time-domain finite difference method to solve the biological heat conduction equation, and can also be solved by frequency domain difference, integration method, finite element method and other methods. According to the form of the solved state transfer equation, a linear Kalman filter or an extended Kalman filter can be selected; δ-ESF provides a linear relationship between temperature change and thermal strain. For a larger range of temperature rise, a nonlinear relationship between temperature change and thermal strain can be used.

[0098] Initialization step: Set the heat source intensity distribution and initialize the adaptive Kalman filter parameters.

[0099] For any thread n, based on the following discrete equation

[0100] AS n (d;t m+1 )=bS n (d;t m )+Q(d;t m ) / k, (7)

[0101] Get t m The heat source term at the moment; if the heat source is time-invariant, the heat source terms Q(d; t m ) to obtain the heat source term of the current thread; average the heat source terms of all threads to obtain the heat source term used in formula (6); if the heat source is time-varying, then at each time t m , the heat source item Q(d; t m ) are averaged to obtain the heat source term used in equation (6).

[0102] The initial thermal strain S0(d), the covariance W(d) of the estimation error of the prediction model, the covariance R(d) of the measurement error, and the covariance P(d) of the filter estimation error are set according to the initial temperature and thermal strain accumulation of the tissue. Among them, the starting values ​​W0(d), R0(d), and P0(d) should be greater than the maximum noise amplitude to include the possible noise amplitude.

[0103] Filtering step: Take the thermal strain sequence of each thread as the measured value of the Kalman filter, and take the calculation result of the biological heat conduction equation as the predicted value of the Kalman filter to obtain the thermal strain estimate S' after adaptive Kalman filtering n ;

[0104] At t0=0, the thermal strain estimate of the Kalman filter is set to S' n (d, t0) = S0(d); for each thread, independently calculate t based on equation (6) m-1 The estimated value of thermal strain at time t is obtained m The predicted value of thermal strain at time t m-1 The covariance P of the system estimation error at the moment n,m-1 , predict t according to the following formula m Covariance of the system estimation error at time

[0105]

[0106] Among them, F T represents the transposed matrix of the state transfer matrix F, W n,m-1 is the covariance matrix of the current estimation error of the prediction model;

[0107] Further, by And the measurement value estimation error R = R0(d), the Kalman gain K is calculated as follows n :

[0108]

[0109] Furthermore, K n Acting on the predicted value and the measured value S n (d;t m ), the filtered t is calculated as follows m The estimated thermal strain S' at time n (d;t m ):

[0110]

[0111] in, is the new information, representing the deviation between the measured value and the predicted value;

[0112] Furthermore, the covariance matrix P of the system estimation error is updated as follows: n,m :

[0113]

[0114] Furthermore, through the latest M iterations of new information v n,m , calculate the new information sequence C as follows m :

[0115]

[0116] Among them, the value of M is an integer greater than 1; further, the covariance matrix W of the current prediction model estimation error is iteratively updated according to the following formula: n,m :

[0117]

[0118] That is, the weights of the model prediction values ​​and the measured values ​​in the final filter output can be adaptively corrected;

[0119] For any thread n, according to the iterative calculation process of equations (7)-(13), make t m From 0 to T, the final estimated image sequence S' n That is the thermal strain output image sequence obtained by this thread after adaptive Kalman filtering.

[0120] It is worth noting that in the process of adaptive parameters, the present invention keeps the covariance matrix R of the measurement error constant and adaptively updates the prediction model error W n,m In fact, we can also keep W constant and adaptively update R n,m .

[0121] It is further worth noting that, when the measurement error covariance R of thermal strain imaging is effectively measured and estimated, a non-adaptive Kalman filter, that is, a Kalman filter in which the covariance matrices W and R are constant during the process, can also be used; in the present invention, a Kalman filter with adaptive parameters based on the properties of new information is used, and an adaptive Kalman filter designed by a method such as multi-mode filtering weighted summation can also be used;

[0122] Fusion step: Based on the system error after thermal strain filtering of each thread, the thermal strain results of each thread are fused according to the scalar weighted linear minimum variance information fusion criterion to output the final thermal strain image sequence S'.

[0123] For any thread n, at any time t m , the system estimation error covariance matrix P updated at the current time n,m , calculate the weight A of the current thread at the current moment according to the scalar weighted linear minimum variance fusion criterion shown in the following formula n,m :

[0124]

[0125] Among them, tr(P n,m ) represents the matrix P n,m Further, the weight A of each thread n,m , the thermal strain distribution of all threads at the current moment after fusion filtering is as follows:

[0126]

[0127] And the covariance of the system estimation error is updated as:

[0128]

[0129] In formula (15), S'(d,t m ) is t m The thermal strain output image after all threads are fused at the moment; by making t m From 0 to T, according to the iteration process of equations (14)-(16), the final output thermal strain image sequence S' can be obtained after adaptive Kalman filtering and thread fusion.

[0130] It should be noted that the weight factor of scalar weighted fusion is obtained by estimating the error P n,m In fact, the fusion of thermal strains of each thread in the present invention is not limited to scalar weighted fusion. Considering the correlation of thermal strain estimation errors of different threads and the affordable complexity of calculation, matrix weighted fusion and component scalar weighted fusion can also be used.

[0131] It is further worth noting that, compared with the known motion compensation, denoising, spatial filtering and other solutions for ultrasonic thermal strain imaging, this embodiment uses the biological heat conduction equation as the prediction model of the adaptive Kalman filter, and iteratively corrects the artifacts and noise in the thermal strain imaging along the time direction. Furthermore, through the linear minimum variance fusion of multi-source information, the thermal strains of multiple threads are combined, and when the strain error of a certain thread is large, an accurate thermal strain image with a good signal-to-noise ratio can still be obtained.

[0132] Example 2

[0133] This embodiment adopts the method in Embodiment 1, uses microwave ablation to heat the visceral fat of living pigs, and calculates the thermal strain during the heating process. The specific steps of this method are as follows:

[0134] Step 1: Figure 2 As shown, during the process of biological tissue heating and ultrasonic image acquisition, the computer 01 controls the image acquisition of the B-ultrasound imager 02 and the power emission of the microwave heater 04; the sampling rate of the B-ultrasound imager is 40MHz, and it is equipped with a linear array probe 03 with 128 array elements and a center frequency of 10.5MHz. The imaging depth is set to 4cm, and focused beam imaging with column-by-column scanning is adopted; the microwave ablation needle 05 is inserted into the fat tissue along the normal direction of the B-ultrasound imaging plane, and its power emission area is located in the ultrasonic imaging plane; 20 seconds of ultrasonic images are continuously acquired at an average frame rate of 50Hz, and each frame of the image corresponds to a spatial range of 40mm in axial depth and 38mm in lateral width.

[0135] Step 2: Calculate the correlation between images through two-dimensional cross-correlation, divide the acquired image sequence into N = 8 threads, and perform image registration and thermal strain imaging on each thread independently to obtain the thermal strain distribution S of each thread. n . Figure 3 a-3c is the cumulative thermal strain distribution of thread n = 1 at t = 5.1, 10.3, and 18.3 seconds. It can be seen that the thermal strain gradually increases with time and the profile is in line with the expected one in the lower left corner of the area of ​​interest of 10 mm in the axial direction and 12.9 mm in the lateral direction; Figure 3 d-3f is the cumulative thermal strain distribution of thread n = 8 at t = 4.5, 9.9, and 19.9 seconds. Figure 3 f The effective thermal strain in the lower left corner is almost swamped by the artifacts above the region.

[0136] Step 3: Use a two-dimensional axisymmetric surface function to calculate the accumulated thermal strain S of each thread n Fitting heat source center O n Finally, align the heat source centers of each thread and place S n Updated to the thermal strain series after position calibration;

[0137] Step 4: Based on the assumption that the microwave ablation needle is a spherical heat source, S n Transform to polar coordinate system S n (r, t); discrete grid points are set along the radial direction in the region of interest, where the spatial step Δr = 0.01 mm, and discrete grid points are set during the heating time, where the time step Δt = 0.1 s; in the prediction model, the state transfer matrix F = A -1 b, system input B(r,t m )=A -1 Q(r,t m ) / k, the coefficient matrix A is:

[0138]

[0139] Where a = κ / Δd 2 ,b=ρC / Δt. Substitute the known fat parameters, such as density ρ=937kg / m 3 , specific heat C = 3258 J / kg / K, thermal conductivity κ = 0.21 W / m / K, and the linear coefficient of temperature change and thermal strain in the δ-ESF model k = 285.71 °C, then the prediction model equation of the Kalman filter was established.

[0140] Step 5. During the entire heating process, the power of the microwave ablation needle is constant. The heat source distribution Q(r) in the fat with the same heating scheme is estimated by thermal strain imaging. The covariance of the initial estimation error of the prediction model is W0=5E, the covariance of the initial measurement error is R0=5E, the covariance of the initial estimation error of the filter is P0=0.1E, the spatial thermal strain S0(r)=0 at the initial moment, and E is a unit matrix of the same size as the coefficient matrix A.

[0141] Step 6: Take the thermal strain sequence of each thread as the measured value of the Kalman filter, and take the calculation result of the biological heat conduction equation as the predicted value of the Kalman filter to obtain the thermal strain estimate S' after adaptive Kalman filtering. n , the process of adaptive Kalman filtering is as follows Figure 4 ,in for Sn,m is the abbreviation of S n (d,t m ). Thermal strain S' after adaptive Kalman filtering n like Figure 5 , Figure 5 a-5c is the filtered thermal strain distribution of thread n=1 at t=5.1, 10.3, and 18.3 seconds. Figure 5 d-5f is the filtered thermal strain distribution of thread n=8 at t=4.5, 9.9, and 19.9 seconds. The noise above and below the right of the shown area is effectively suppressed, which significantly improves the signal-to-noise ratio of the thermal strain image and presents a thermal strain distribution that is consistent with the characteristics of the microwave ablation needle heat source.

[0142] Step 7: According to the scalar weighted linear minimum variance information fusion criterion, the thermal strains of each thread are fused to obtain the final thermal strain estimate S'. Figure 6 a is the thermal strain distribution S' at t = 20s. It can be seen that the final thermal strain is not affected by the large amount of noise in the thermal strain of thread n = 8. The noise in the shown area is basically filtered out, and the thermal strain has a high signal-to-noise ratio and spatial resolution. Figure 6 b is the thermal strain-time curve at the heat source position. It can be seen that with heating, the thermal strain curve densely outputs thermal strain values ​​at intervals of Δt = 0.1s, and its amplitude gradually increases, showing a similar trend to the experimental measurement and numerical simulation results of the biological heat transfer model.

[0143] The present invention has been described in detail above in conjunction with specific exemplary embodiments. However, it should be understood that various modifications and variations may be made without departing from the scope of the present invention as defined by the appended claims. The detailed description and the accompanying drawings should be considered only as illustrative and not restrictive, and any such modifications and variations, if any, will fall within the scope of the present invention described herein. In addition, the background art is intended to illustrate the current status and significance of the present technology and is not intended to limit the present invention or the application field of the present application and the present invention.

Claims

1. A multi-threaded in vivo ultrasonic thermal strain imaging method, characterized in that: The following steps are involved: Acquisition steps: heating the target area, acquiring ultrasound images of the in vivo hyperthermia target area, and obtaining a digital ultrasound image sequence; Imaging step: The digital ultrasound image sequence obtained in the acquisition step is divided into multiple quasi-periodic threads corresponding to different motion phases through correlation analysis n ,in n =1,2,…, N , N A positive integer represents the total number of threads; According to the target area position, the region of interest is selected in the image, and image registration and thermal strain imaging are performed independently in each thread based on the region of interest to obtain the thermal strain distribution sequence of each thread S n ; Calibration steps: S n The heat source center position of the corresponding thread is obtained by fitting the two-dimensional axisymmetric surface function. By aligning the heat source centers of each thread, S n Updated to the thermal strain distribution sequence after position calibration; Modeling steps: Solve the biological heat conduction equation by using the finite difference time domain method and construct the prediction model of the Kalman filter; Initialization steps: set the heat source intensity distribution; initialize the adaptive Kalman filter parameters; Filtering step: Take the thermal strain sequence of each thread as the measured value of the Kalman filter, and take the calculation result of the biological heat conduction equation as the predicted value of the Kalman filter to obtain the thermal strain estimation after adaptive Kalman filtering. S' n ; Fusion step: Based on the system error after thermal strain filtering of each thread, the thermal strain results of each thread are fused according to the scalar weighted linear minimum variance information fusion criterion to output the final thermal strain image sequence. S' ; The specific implementation of the modeling step is: For any thread n , with the heat source as the center, establish the size A × L , A The axial dimension ，L The horizontal dimension ， Not less than the region of interest, the maximum spatial step length Δ d 、Time span 0~ T , time step Δ t The discrete grid, Δ d The value is less than A / 2 and L / 2, the smaller value, Δ t The value is less than T / 2; the prediction model is (3) in, and They are t m = m Δ t The predicted thermal strain value and The Kalman filter thermal strain estimate at time t; d The coordinates in the space are in the form of standard Cartesian coordinate system, two-dimensional polar coordinate system, three-dimensional cylindrical coordinate system or three-dimensional spherical coordinate system; when the thermal strain of each thread is not defined in the coordinate system d , then first transform it to d ,get S n ( d ; t ); F = A -1 b is the state transition matrix of the system, A -1 For the matrix A The inverse matrix of for t m System input at time A and b is the coefficient matrix related to tissue specific heat capacity, thermal conductivity, and density, which is obtained by discretizing the biological heat conduction equation using the finite-difference time-domain method; k It is a constant specified in the thermal strain temperature measurement model δ-ESF, which describes the linear relationship between temperature change and thermal strain; is the heat source term; The specific implementation of the filtering step is: exist t At time 0=0, the thermal strain estimate of the Kalman filter is set to S' n ( d , t 0)= S 0( d ); for each thread, independently based on (3) by t m−1 The estimated value of thermal strain at time t m The predicted value of thermal strain at time t m−1 Covariance of the system estimation error at time P n,m−1 , according to the following prediction t m Covariance of the system estimation error at time : ，(5) in, F T Represents the state transfer matrix F The transposed matrix of W n,m−1 is the covariance matrix of the current estimated error of the prediction model Depend on and the measurement value estimation error R = R 0( d ), the Kalman gain is calculated as follows K n : (6) Will K n Acting on the predicted value and measured values S n ( d ; t m ), and the filtered t m Estimated thermal strain at time S' n ( d ; t m ): ，(7) in, is the new information, representing the deviation between the measured value and the predicted value; further, the covariance matrix of the system estimation error is updated as follows: : ；(8) By recent M The new information of the iteration v n,m , calculate the new information sequence according to the following formula C m : ，(9) in, M The value of is an integer greater than 1; further, the covariance matrix of the current prediction model estimation error is iteratively updated according to the following formula: W n,m : ，(10) Adaptively correct the weights of model predictions and measurements in the final filter output; For any thread n , according to the iterative calculation process of formula (4)-(10), t m From 0 to T , the final estimated image sequence S' n That is the thermal strain output image sequence obtained by this thread after adaptive Kalman filtering.

2. The multi-thread fusion in vivo ultrasonic thermal strain imaging method according to claim 1, characterized in that: The target area is heated in the acquisition step by one or more of microwave ablation, radiofrequency ablation, laser ablation, ultrasonic ablation and infrared ablation.

3. The multi-thread fusion in vivo ultrasonic thermal strain imaging method according to claim 1, characterized in that: The S n The specific implementation method for fitting the two-dimensional axisymmetric surface function is: take the natural logarithm of the thermal strain of each thread ln S n , polynomial least squares is used to fit it according to the following formula: (1) Get the fitting parameters a , b , c , d , e ; then the heat source center position O n The vertical and horizontal coordinates are determined as (2)。 4. The multi-thread fusion in vivo ultrasonic thermal strain imaging method according to claim 1, characterized in that: In the initialization step, the specific implementation method of setting the heat source distribution is: For any thread n , based on the following discrete equation (4) get t m The heat source term at the moment; if the heat source is time-invariant, the heat source terms of all time serial numbers Average to get the heat source item of the current thread; average the heat source items of all threads to get the heat source item used for the prediction model; If the heat source is time-varying, then at each moment t m , the heat source items of each thread The average is used to obtain the heat source term used in the prediction model.

5. The multi-thread fusion in vivo ultrasonic thermal strain imaging method according to claim 1 or 2, characterized in that: In the initialization step, the specific implementation method of initializing the adaptive Kalman filter parameters is: Set initial thermal strain based on initial temperature and thermal strain accumulation of tissue S 0( d ), the covariance of the prediction model estimation error W ( d ), covariance of measurement error R ( d ) and the covariance of the filter estimation error P ( d ), where the starting value W 0( d ), R 0( d ), P 0( d ) should be greater than the maximum noise amplitude.

6. The multi-thread fusion in vivo ultrasonic thermal strain imaging method according to claim 1, characterized in that: The specific implementation method of the fusion step is: For any thread n , at any time t m , the system estimation error covariance matrix updated at the current time P n,m , calculate the weight of the current thread at the current moment according to the scalar weighted linear minimum variance fusion criterion shown in the following formula A n,m : (11) Among them, tr( P n,m ) represents the matrix P n,m trace; further, the weight of each thread A n,m , the thermal strain distribution of all threads at the current moment after fusion filtering is as follows: ，(12) And the covariance of the system estimation error is updated as: ；(13) In formula (12), S' ( d , t m ) is t m The thermal strain output image after all threads are fused at the moment; by using t m From 0 to T , according to the process iteration of equations (11)-(13), the final output thermal strain image sequence is obtained after adaptive Kalman filtering and thread fusion. S' .

7. A multi-threaded in vivo ultrasonic thermal strain imaging device, characterized in that: include: Computer, B-ultrasound machine, linear array probe and thermal ablation equipment; The B-ultrasound instrument is connected to the linear array probe and is used to collect ultrasound images of living objects; Thermal ablation devices are used to heat living subjects, creating displacement and strain distributions in tissue; The computer is connected to the B-ultrasound instrument and the thermal ablation device respectively, and is used to control the B-ultrasound instrument and the thermal ablation device, and is configured to implement the multi-threaded fusion in vivo ultrasonic thermal strain imaging method as described in any one of claims 1-6.

Images