A method and device for estimating the distribution of internal heat sources based on ultrasonic imaging

Abstract

The present invention discloses a method and device for estimating the distribution of internal heat sources based on ultrasonic imaging, belonging to the field of medical image processing. The method of the present invention is to collect a digital ultrasonic image sequence in a living body, perform multi-thread partitioning, image registration, thermo-strain imaging, and denoise the thermo-strain image based on spatial correlation. Fit the position of the heat source in the thermo-strain image based on a surface function, and accordingly correct the thermo-strain distribution and compensate the initial strain values of each thread. Subsequently, the estimation of the heat source distribution from the thermo-strain image sequence is realized by the finite-difference time-domain method. For different hyperthermia processes using the same heating scheme, the real-time prediction of the thermo-strain and temperature field distributions in tissues is realized. The present invention does not require the aid of invasive monitoring devices, can effectively estimate the spatial distribution of heat sources in tissues, and predict the thermo-strain and temperature distributions based on this, providing an important means for clinical hyperthermia to evaluate the quantity-effect relationship.

Description

Technical Field

[0001] The present invention relates to the field of medical image processing, and more specifically, to a method and device for estimating the distribution of heat sources in vivo based on ultrasonic imaging. Background Art

[0002] Currently, widely used clinical cancer treatment methods such as surgical treatment, chemotherapy, and radiotherapy have disadvantages such as large toxic side effects, high recurrence rate, and poor specificity. In line with the purpose of reducing toxic side effects and treatment costs, the development of thermal therapy has received much attention in recent years. The thermal dose is an important indicator related to the effectiveness and safety of thermal therapy. Therefore, clinically, dynamically planning and adjusting treatment plans through real-time temperature monitoring is an urgent problem to be solved.

[0003] The bioheat conduction equation (also known as the Pennes equation) is a heat transfer model generally applicable to biological tissues. Given parameters such as the density, specific heat capacity, and thermal conductivity of biological tissues, and the distribution of heat source intensity, the temperature distribution in tissues can be predicted in real time through the bioheat conduction equation. However, in fact, the distribution of effective heat sources in vivo often deviates from theory and simulation, and due to physiological movements, it is difficult to accurately measure the distribution of heat source intensity.

[0004] In the prior art, some technical solutions have been proposed for measuring heat sources in vivo. For example, the invention title is: Analysis Method and Device for the Heat Source Intensity and Temperature Distribution Inside an Organism Based on a Point Heat Source Model (application date: September 29, 2017; application number: 201710908792.6). This solution discloses an analysis method and device for the heat source intensity and temperature distribution inside an organism based on a point heat source model. Aiming at the difficulty that the heat in an organism often cannot be measured in vivo with instrument equipment, based on the assumptions of a steady-state and point heat source accumulation temperature field, the relationship between the point heat source and the spatial temperature field is derived. Subsequently, an infrared CCD detector is used to measure the body surface temperature, and the heat source distribution information is obtained through Lorentz linear fitting of it. Finally, based on the principle of linear accumulation of the position and intensity of the heat source, a three-dimensional temperature field estimation is obtained. This invention provides an objective basis for biomedical analysis. Its disadvantages are that it only considers the steady-state bioheat conduction model and is not applicable to predicting the real-time changing temperature of thermal therapy. Secondly, it ignores the influence of heat convection on the body surface temperature, which may lead to errors in fitting the in vivo heat source information based on the body surface temperature, thereby resulting in errors in the in vivo temperature field estimation.

[0005] In summary, how to effectively measure the heat source in vivo, obtain the position and intensity information of the heat source, and accurately predict the thermal strain and temperature field based on the heat source information is an urgent problem to be solved in the prior art. Summary of the Invention

[0006] 1. Technical problems to be solved

[0007] In view of the problems in the prior art that it is difficult to accurately measure the effective heat source in living tissues, and due to factors such as living physiological movements, the error in heat source estimation is further caused, resulting in the deficiency of temperature field prediction based on the heat source, the present invention provides a method and device for estimating the distribution of in-vivo heat sources based on ultrasonic imaging, which can effectively estimate the heat source information, and based on the estimated heat source, the thermal strain and temperature field changing with time can be predicted in real time.

[0008] 2. Technical solutions

[0009] The object of the present invention is achieved through the following technical solutions.

[0010] The present invention provides a method for estimating the distribution of in-vivo heat sources based on ultrasonic imaging, including the following steps:

[0011] Acquisition step: Using a hyperthermia method to heat the target area, collecting ultrasonic images of the living hyperthermia target area to obtain a digital ultrasonic image sequence I;

[0012] Imaging step: Selecting an area of interest according to the position of the target area, and dividing the image sequence I into multiple threads u corresponding to different motion phases (u = 1, 2,..., U), where U is the total number of threads, through correlation analysis, image registration based on the area of interest, and thermal strain imaging, to obtain a thermal strain image sequence S of the corresponding thread u ; Constructing a denoising function according to the correlation between threads to denoise S u to obtain an updated S u ;

[0013] Fitting step: Performing two-dimensional axisymmetric surface fitting on S u to obtain the heat source center position P of each thread u ;

[0014] Calibration step: By moving the thermal strain image sequence S corresponding to each thread u , making the heat source center position P of each thread u u coordinate-consistent;

[0015] Averaging step: Based on the finite difference method in the time domain and the bioheat conduction equation, establishing a quantitative relationship between the heat source distribution and the thermal strain distribution, estimating the heat source distribution of each thread, and averaging the estimation results of all threads to obtain the average heat source distribution.

[0016] Prediction step: Based on the position of the actual heat source, the average heat source distribution obtained in the averaging step, and the quantitative relationship between the heat source distribution and the thermal strain distribution, performing real-time prediction on the thermal strain and temperature field distribution changing with time.

[0017] Furthermore, the specific implementation of the fitting step is as follows:

[0018] At any time, the thermal strain distribution of thread u is denoted as S u (x, y; t) (u = 1, 2, …, U, where U is the total number of threads), and x, y, and t are the horizontal, vertical, and time coordinates respectively;

[0019] Perform a two-dimensional axisymmetric surface fitting on S u (x, y; T u ) to obtain the heat source center position P of each thread u ;

[0020] A more efficient fitting method is: perform a well-known polynomial least squares fitting on ln[S u (x, y; T u )] (ln is the natural logarithm function with base e) according to the following function f(x, y),

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

[0022] Obtain the fitting parameters a, b, c, d, e; then the vertical and horizontal coordinates of the heat source center position P u are respectively determined as

[0023]

[0024] Furthermore, the specific implementation of the calibration step is as follows:

[0025] Taking any thread as a reference, move the thermal strain images of other threads in the (x, y) plane until the heat source center coordinates of all threads are the same.

[0026] Furthermore, another specific implementation of the calibration step can also be achieved by:

[0027] Performing cross-correlation calculations based on a sliding window on the digital ultrasonic images and thermal strain image sequences corresponding to each thread to obtain a displacement field, and moving S u according to the displacement field, so that the heat source center position P u coordinates of each thread u are the same.

[0028] Furthermore, in the acquisition step, when collecting ultrasonic images of the in-vivo hyperthermia target area, collect and store ultrasonic image sequences of multiple respiratory cycles before starting heating, so that the first frame image of each thread occurs before heating.

[0029] Further, if the ultrasonic image acquisition in the acquisition step starts only after heating, then at the end of the calibration step, it is also necessary to compensate for the strain at the starting moment of each thread. The compensation method is as follows:

[0030] Taking the thread with the earliest start time as the source thread and other threads as target threads respectively, using the first and second frame thermal strain images of the source thread, perform time-domain linear interpolation at the moment of the first frame thermal strain image of the target thread, and the obtained image is used as the compensation image; or perform time-domain linear extrapolation along the negative direction of t for the moment of the first frame thermal strain of the source thread through the first and second frame thermal strain images of the target thread, and the obtained image is taken as negative as the compensation image; add the compensation image to all the thermal strain images of the target thread, then the thermal strain image sequence S u is updated to the compensated result.

[0031] Further, at the end of the calibration step, the thermal strain image sequence S u is transformed from the Cartesian coordinate system (x, y) to the polar coordinate system (r, θ). It should be noted that calculations with a similar principle and method can be carried out in the rectangular coordinate system, but the process will become more cumbersome. In practical applications, considering the symmetry of the heat source, the polar coordinate system is a more widely used coordinate system, which is conducive to the simplification of derivation and calculation.

[0032] Further, the specific implementation method of the averaging step is as follows:

[0033] For any thread, with the heat source position as the center, establish a discrete grid with dimensions A (axial) × L (transverse) not less than the region of interest, a maximum spatial step size Δd, a time span of 0 to T, and a time step size Δ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;

[0034] The thermal strain sequence S u (r, θ; t) is interpolated in the spatial domain and time domain so that the sequence is sampled on the above discrete grid; based on the following discrete equation, the heat source term is obtained:

[0035] AS u (r, θ; t + Δt) = BS u (r, θ; t) + Q u (r, θ; t) / k, (3)

[0036] where A and B are coefficient matrices related to tissue specific heat capacity, thermal conductivity, and density, which can be obtained by solving the bioheat conduction equation using the finite difference method in the time domain; k is a constant specified in the thermal strain temperature measurement model δ-ESF that describes the linear relationship between the temperature change and the thermal strain; the heat source distribution Q u (r, θ) of this thread is the Q predicted at each time step uThe result after averaging (r, θ; t) over t is obtained by averaging the heat source term Q obtained for all threads u over (r, θ) to obtain the average heat source term Q(r, θ).

[0037] Furthermore, the specific implementation of the prediction step is as follows:

[0038] Adopt the same heating scheme as in the acquisition step to heat the same type of living tissue, set the actual heat source center to the origin, and iteratively calculate the thermal strain distribution S(r, θ; t) in the tissue according to the following formula:

[0039] AS(r, θ; t+Δt) = BS(r, θ; t) + Q(r, θ) / k, (4)

[0040] Predict that the temperature distribution in the tissue is T(r, θ; t) = T 0 + kS(r, θ; t), where T 0 is the baseline temperature in the biological tissue before heating.

[0041] If the heating schemes adopted in different hyperthermia processes are the same, but the types of biological tissues treated are different, then the parameters A and B in the prediction step need to be updated, and can be obtained by discretizing the bioheat conduction equation in the spatial and temporal domains using known methods according to the physical constants of the actual biological tissue.

[0042] The present invention also provides a device for estimating the in-vivo heat source distribution based on ultrasonic imaging, including: a computer, a B-ultrasound instrument, a linear array probe, and a thermal ablation device;

[0043] The B-ultrasound instrument is connected to the linear array probe and is used for collecting ultrasonic images of a living object;

[0044] The thermal ablation device is used for heating a living object to form displacement and strain distributions in the tissue;

[0045] The computer is respectively connected to the B-ultrasound instrument and the thermal ablation device, and is used for controlling the B-ultrasound instrument and the thermal ablation device, and is configured to implement the above method for estimating the in-vivo heat source distribution based on ultrasonic imaging.

[0046] 3. Beneficial effects

[0047] Compared with the prior art, the advantages of the present invention are: adopting the ultrasonic image registration algorithm under multiple motion phases, a thermal strain distribution with good motion suppression effect and high signal-to-noise ratio can be obtained; further, by solving the bioheat conduction equation under the heat source model, the position and intensity information of the heat source inside the tissue can be effectively estimated from the thermal strain. Based on this, the real-time prediction of the three-dimensional thermal strain and temperature field changing with time in the tissue can be further realized, which has important application value for clinical treatment and diagnosis based on thermal effects. Description of the drawings

[0048] Figure 1 It is a schematic flow chart of a method and device for estimating the distribution of internal heat sources based on ultrasonic imaging;

[0049] Figure 2 It is a schematic diagram of ultrasonic image acquisition during the heating process of biological tissue in Example 2;

[0050] Figure 3 It is a schematic diagram of the cumulative thermal strain distribution of 8 threads in Example 2;

[0051] Figure 4 It is the distribution of heat source positions of 8 threads obtained by surface fitting and the normalized amplitude of the thermal strain surface S 1 and the schematic diagram of the thermal strain surface obtained by fitting in Example 2;

[0052] Figure 5 It is the heat source distribution estimated by 8 threads in Example 2 and the schematic diagram of the averaged heat source distribution;

[0053] Figure 6 It is the temperature distribution in the xy plane and the temperature contour lines along the x - direction and y - direction at different times over the heat source position when heating for t = 20s;

[0054] Reference numerals: 01, portable B - ultrasound system; 02, computer; 03, linear array probe; 04, microwave ablation needle; 05, microwave heater. Detailed implementation manners

[0055] To make the objectives, 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 with reference to the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are some, but not all, of the embodiments of the present invention; moreover, the various embodiments are not relatively independent and can be combined with each other as needed to achieve better effects. 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 claimed present invention, but merely represents selected embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts fall within the protection scope of the present invention.

[0056] To further understand the content of the present invention, the present invention will be described in detail with reference to the accompanying drawings and embodiments.

[0057] Example 1

[0058] The schematic flow chart of the method of the example is as Figure 1 , A method and device for estimating the distribution of internal heat sources based on ultrasonic imaging according to the present invention includes the following steps:

[0059] Acquisition step: Heat the target area using a hyperthermia method, perform ultrasonic image acquisition on the in-vivo hyperthermia target area, and obtain a digital ultrasonic image sequence I.

[0060] Perform B-mode ultrasound imaging on the in-vivo hyperthermia target area, and continuously acquire a digital ultrasonic image sequence at equal time intervals; the B-mode ultrasound instrument operates in a real-time imaging mode, adjust the angle of the probe and the parameters of the ultrasonic system to make the target area appear in the imaging field of view; it should be noted that in the embodiment, according to the actual situations such as the coupling state between the body surface of the in-vivo and the probe, the depth of the target area, and the acoustic path from the body surface to the target area, imaging probes with different beam control methods, different numbers of array elements, and different geometric structures can be selected, as well as different B-mode ultrasound instrument imaging parameters and imaging modes.

[0061] Heat the target area using hyperthermia methods such as focused ultrasound and radiofrequency ablation; it should be noted that heating the biological tissue target area can be in a non-invasive or invasive mode, and the specific heating methods include microwave ablation, radiofrequency ablation, laser ablation, ultrasound ablation, and infrared ablation; furthermore, in order to control the heating area to fall within the imaging plane, for invasive heating methods, the relative position between the heating area of the heating needle tip and the imaging plane can be controlled by real-time imaging of the heating needle using an ultrasound instrument, and for non-invasive heating methods, it can be achieved through the coaxial design of the heating device and the imaging system.

[0062] Continuously acquire target area images at equal time intervals, and the acquisition duration T is not shorter than two respiratory cycles to obtain a digital ultrasonic image sequence. It should be noted that the acquired digital ultrasonic image sequence can be a channel signal, an ultrasonic radio frequency (RF) signal, a complex-valued analytic signal after quadrature demodulation or Hilbert transform, an envelope signal obtained by taking the amplitude of the analytic signal, a grayscale image signal after further logarithmic compression, and other signals obtained by compressing and filtering the aforementioned signals according to the signal type.

[0063] Imaging step: Select the region of interest according to the target area position, and divide the image sequence I into multiple threads u corresponding to different motion phases (u = 1, 2,..., U), where U is the total number of threads, through correlation analysis, image registration based on the region of interest, and thermo-strain imaging, to obtain the thermo-strain image sequence S of the corresponding thread u ; Construct a denoising function based on the correlation between threads for S u Denoise to obtain the updated S u ;

[0064] Calculate the correlation coefficient between images to determine the respiratory frequency and the characteristic frame of the expiratory state; specifically, take M frames of images containing at least one respiratory cycle in the ultrasonic image sequence as the template frames X in turn, and take all the acquired N frames of images as the detection frames Y, and obtain the cross-correlation coefficient - time curve γ through normalized cross-correlation calculationm (n),

[0065]

[0066] where n = 1, 2, …, N, m = 1, 2, …, M. COV(·) represents covariance, σ X and σ X are the standard deviations of X and Y. The M curves γ m (n) are subjected to a linear amplitude modulation z-transform to obtain the corresponding spectral curves Γ m (f), where f is the frequency variable, and the maximum amplitude of each spectral curve is denoted as A m . The frequency corresponding to the maximum value in A m is the respiratory frequency f res , and the template frame number m = p corresponding to the curve where it is located, then the p-th frame is the exhalation state feature frame.

[0067] According to the cross-correlation coefficient-time curve γ p (n) corresponding to the exhalation feature frame p, the image sequence is divided into different cycles and multiple motion phase periods of exhalation and inhalation states; the cut-off frequency f C = f res of the finite impulse response low-pass filter is set, and the filter order is n tap . The filter is applied to the cross-correlation coefficient-time curve γ p (n) with the exhalation state feature frame as the template frame to obtain the curve γ' p (n) with the high-frequency components filtered out. Along the time direction, the maxima and minima of the curve γ' p (n) are alternately retrieved. Every two adjacent maxima form the boundary of an inhalation cycle, and every two adjacent minima form the boundary of an exhalation cycle. Each exhalation and inhalation cycle is equally divided into U 1 , U 2 phase periods, and a total of U different motion phase periods are obtained.

[0068] Based on the image registration of the region of interest in the same phase period of different cycles, the registered multi-threaded image sequence Q u (u = 1, 2, …, U) is obtained; the images at the middle moment of the first cycle of each phase period are selected, and according to the position of the target area, the center point coordinates are further selected as (i 0 , j 0 ), and the region of interest with a size of A 0 ×L 0 is used as the reference frame Q ref , and it is placed in the subset Q u . In the same phase period of the next cycle, for each frame of the image, a region of interest with a size of A 0 ×L 0, the area with the center point located at (i 0 +δa, j 0 +δl) is the detection frame Q tar . Respectively vary δa and δl within the ranges of -A 0 / 3 to A 0 / 3 and -L 0 / 3 to L 0 / 3. δa and δl are the displacement amounts in the axial and transverse directions respectively, so that the two-dimensional cross-correlation coefficient between Q tar and Q ref reaches the maximum value. The corresponding δa and δl are the displacement amounts Δa and Δl of Q tar . Shift Q tar by Δa and Δl in the axial and transverse directions respectively to obtain the registered image Q' tar . Compare all the maximum cross-correlation coefficients, and take the Q' tar corresponding to the maximum value as the new Q ref , and place it at the end of the subset Q u . Use the new Q ref to perform the same operation on the next cycle, and so on until all cycles are registered, and traverse all phase periods to obtain the multi-threaded image subset Q u .

[0069] The multi-threaded image subset Q u calculates the thermal strain in parallel, and constructs a denoising function based on the inter-thread correlation to obtain the thermal strain image sequence S u after denoising the thermal strain; after each thread image subset Q u obtains the displacement image sequence based on the well-known block matching algorithm, and obtains the thermal strain image sequence S u by taking the axial difference. Taking the spatial point (i, j) of the image S u as the center, take a rectangular window R 1 ×L 1 with a size of A u , and calculate the average cross-correlation coefficient of its corresponding part R v with the strain S v of other threads (v = 1, 2,..., U, v ≠ u) according to the following formula,

[0070]

[0071] where, COV(R u , R v ) represents the covariance between R u and R v , and σ R1 and σ Rv are the standard deviations of R u and R v . The denoising function q u(i,j) is defined as:

[0072]

[0073] where α is the attenuation coefficient, and its value range can be 10 -2 ~10 2 , p th is the correlation coefficient threshold, and its value range can be 0.4 to 0.9. Multiply the denoising function q u (i,j) by S u , then S u is updated to the denoised thermal strain image.

[0074] Fitting step: Perform a two-dimensional axisymmetric surface fitting on S u to obtain the heat source center position P u of each thread; At any moment, denote the thermal strain distribution of thread u as S u (x,y;t) (u = 1,2,…,U, where U is the total number of threads), and x, y, t are the horizontal, vertical, and time coordinates respectively; The last frame thermal strain image of this thread is S u (x,y;T u ), and T u is the corresponding time coordinate;

[0075] Perform a well-known polynomial least squares fitting on ln[S u (x,y;T u )] (ln is the natural logarithm function with base e) according to the following function f(x,y),

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

[0077] to obtain the fitting parameters a, b, c, d, e; Then the vertical and horizontal coordinates of the heat source center position P u are determined respectively as

[0078]

[0079] It should be noted that the fitting function can also be replaced by other functions with two-dimensional axisymmetric distributions, including but not limited to two-dimensional conical functions, two-dimensional cosine functions, two-dimensional Hamming window functions, two-dimensional Blackman window functions; Further, the listed fitting method can also be replaced by directly using S u (x,y;T u ) for the fitting of the two-dimensional axisymmetric surface function, but using ln[S u (x,y;T u)]Performing fitting has higher computational efficiency; fitting the strain centers in the x and y directions respectively based on the maximum points of two-dimensional thermal strain can obtain the two-dimensional strain center more quickly, but its accuracy is lower than that of the two-dimensional surface fitting method.

[0080] Calibration steps: By moving the thermal strain image sequences S corresponding to each thread u , make the heat source center positions P of each thread u u have the same coordinates;

[0081] By moving the images, make the heat source center positions of each thread the same, and realize the correction of the image positions in the thermal strain sequences of each thread. Specifically, taking any thread as a reference, move the thermal strain images of other threads on the (x, y) plane until the heat source center coordinates of all threads are the same; it should be noted that the correction of the spatial positions of the thermal strains of the above-mentioned threads can also be obtained by calculating the displacement field through the cross-correlation of the B-ultrasound images and thermal strain images corresponding to each thread based on a sliding window, and then shifting the thermal strain distribution.

[0082] If in the acquisition step, the acquisition is not carried out before heating starts, then the first frame images of each thread occur after heating. In this case, after the correction of the image positions in the thermal strain sequences of each thread in the calibration step, it is also necessary to compensate for the strains at the starting moments of each thread. Specifically, taking the thread with the earliest starting time as the source thread and other threads as target threads respectively, using the first and second frame thermal strain images of the source thread, perform time-domain linear interpolation at the moment when the first frame thermal strain image of the target thread is located, and the obtained image is used as the compensation image; add the compensation image to all the thermal strain images of the target thread, then the thermal strain image sequence S of the target thread u is updated to the compensated result; it should be noted that it is also possible to perform time-domain linear extrapolation of the moment when the first frame thermal strain of the source thread is located along the negative direction of t through the first and second frame thermal strain images of the target thread, and take the negative of the obtained image as the compensation image. Further, it should be noted that if the ultrasound image sequences of multiple respiratory cycles are acquired and stored before heating starts, such that the first frame images of each thread occur before heating, then the above-mentioned strain compensation process can be omitted.

[0083] Furthermore, transform the thermal strain images of all threads from the Cartesian coordinate system (x, y) to the polar coordinate system (r, θ). It should be noted that in the present invention, the coordinate transformation is performed according to the characteristics such as the spatial symmetry of the heat source model, so subsequent calculations can also be directly based on rectangular coordinates. Because the calculation process will be more cumbersome in the rectangular coordinate system, transforming to polar coordinates is beneficial to the simplification of derivation and calculation.

[0084] Average step: Based on the finite-difference time-domain method and the bio-heat transfer equation, establish a quantitative relationship between the heat source distribution and the thermal strain distribution, estimate the heat source distribution of each thread, and average the estimation results of all threads to obtain the average heat source distribution.

[0085] Based on the finite-difference time-domain method and the well-known bio-heat transfer equation, establish a quantitative relationship between the heat source distribution and the thermal strain distribution; specifically, for any thread, with the heat source position as the center, establish a discrete grid with dimensions A (axial) × L (transverse) not less than the region of interest, a maximum spatial step Δd, a time span of 0 to T, and a time step Δt. The value of Δd is less than the smaller of A / 2 and L / 2, and the value of Δt is less than T / 2.

[0086] Furthermore, interpolate the thermally strained sequence S u (r, θ; t) in the spatial and time domains so that the sequence is sampled on the above discrete grid; it should be noted that the sampling of the thermally strained sequence of each thread on the discrete grid can be achieved through nearest neighbor, linear, cubic spline, and cubic interpolation. Interpolate the thermally strained sequence S u (r, θ; t) in the spatial and time domains so that the sequence is sampled on the above discrete grid; based on the following discrete equation, obtain the heat source term:

[0087] AS u (r, θ; t + Δt) = BS u (r, θ; t) + Q u (r, θ; t) / k,

[0088] where A and B are coefficient matrices related to tissue specific heat capacity, thermal conductivity, and density, which can be obtained by solving the bio-heat transfer equation using the finite-difference time-domain method; it should be noted that the bio-heat transfer equation can be solved by methods such as finite-difference time-domain, frequency-domain difference, and integral method; furthermore, it should be noted that in the finite-difference time-domain method, the first-order differential term has different forms of forward, backward, and central differences, and the specific solution can be explicit or implicit. k is a constant specified in the well-known thermally strained temperature measurement model δ-ESF that describes the linear relationship between the temperature change and the thermal strain.

[0089] Further, estimate the heat source distribution of each thread, and average the estimation results of all threads to obtain the average heat source distribution; specifically, the heat source distribution Q u (r, θ) of this thread is the result of averaging Q u (r, θ; t) predicted at each time step with respect to t. By averaging the heat source terms Q u (r, θ) obtained for all threads, the average heat source term Q(r, θ) is obtained.

[0090] It should be noted that the δ-ESF model is only applicable to the temperature range of 37 to 50 °C. For a larger range of temperature increases, a quantitative relationship between thermal strain and temperature can be established based on a non-linear temperature increase model.

[0091] Prediction step: Based on the position of the actual heat source, the average heat source distribution obtained in the averaging step, and the quantitative relationship between the heat source distribution and the thermal strain distribution, the time-varying thermal strain and temperature field distribution are predicted in real time.

[0092] For different hyperthermia processes using the same heating scheme, based on the position of the actual heat source, the average heat source distribution obtained in the averaging step, and the relationship between the heat source term and the thermal strain distribution, the time-varying thermal strain and temperature field distribution are predicted in real time. Specifically, using the same heating scheme as in the acquisition step, heating the same type of biological tissue; setting the center of the actual heat source as the origin, and according to the following formula, iteratively calculate the thermal strain distribution S(r,θ;t) in the tissue,

[0093] AS(r,θ;t + Δt) = BS(r,θ;t) + Q(r,θ) / k;

[0094] It should be noted that if different hyperthermia processes use the same heating scheme but different types of biological tissues are treated, the parameters A and B in the above formula need to be updated, which can be obtained by discretizing the bioheat conduction equation in the spatial and temporal domains using known methods according to the physical constants of the actual biological tissue; the above prediction process of thermal strain and temperature involves the inverse operation of the coefficient matrix A, and the inverse operation of a large-size matrix is accompanied by a non-negligible computational burden. In particular, when considering that the tissue heat capacity, thermal conductivity, and density vary with temperature, the coefficient matrix A needs to be updated in real time with temperature, further increasing the computational load. Therefore, combining the characteristic that the coefficient matrix A is a tridiagonal matrix, the present invention adopts the chasing method based on LU decomposition to solve, which has higher computational efficiency.

[0095] Furthermore, it should be noted that the thermal strain calculation based on multi-threading in the present invention can effectively suppress the thermal strain artifacts caused by physiological movement, thereby effectively improving the estimation accuracy of the heat source term. At the same time, the multi-threading heat source estimation fusion avoids the dependence on a certain motion phase and improves the accuracy and robustness of the heat source estimation method; in the present invention, a method and device for estimating the in-vivo heat source distribution based on ultrasonic imaging have the advantage of non-invasiveness, overcome the deficiency in the prior art that it is difficult to accurately measure the effective heat source in vivo, and can estimate the position and intensity of the heat source while predicting the thermal strain and temperature field in vivo under the same heating scheme.

[0096] Example 2:

[0097] This embodiment is based on the method of Embodiment 1. A microwave ablation needle is used to heat the visceral fat of a live pig, estimate the heat source distribution in biological tissues, and predict the thermal strain and temperature field in the organism under the same heating scheme. The specific steps of this method are as follows:

[0098] Step 1. As Figure 2 shown, in this embodiment, a portable B-ultrasound system 01 with a sampling rate of 40 MHz is paired with an L12-5 linear array probe 03 with 128 array elements and a center frequency of 10.5 MHz to image the target visceral fat of a live pig. Ultrasonic radio frequency image data is continuously and equidistantly collected. The imaging depth is set to 4 cm, the line density is set to 128, and the imaging mode is set to the fundamental wave mode; the thermal ablation device is used to heat the live object to form displacement and strain distributions in the tissue. In this embodiment, the thermal ablation device uses a microwave ablation needle 04 and a microwave heater 05; the microwave ablation needle 04 is inserted into the target fat along the normal direction of the imaging plane under the guidance of real-time B-ultrasound images; the computer 02 synchronously controls the radio frequency data acquisition of the B-ultrasound imaging device and the power emission of the microwave heater 05. Each experiment collects an image sequence of the heating process with a length of 20 seconds. The ultrasonic images are averaged at 50 frames per second, and each frame of the image corresponds to a tissue area of 40 mm in the longitudinal direction and 38 mm in the transverse direction. The collected ultrasonic radio frequency data is transmitted to the computer for processing.

[0099] Step 2. Perform two-dimensional cross-correlation calculation on the image sequence to obtain the correlation coefficient-time curve. Through time-domain and frequency-domain analysis, the image sequence is divided into 8 threads in 8 motion phases; an interested region with a size of 8.8*10.8 mm is selected, and each thread performs image registration, calculates thermal strain, and denoises in parallel. The processed thermal strain S u As Figure 3 shown, Figure 3 in (a)- Figure 3 in (h) are the cumulative thermal strain distributions S u of the obtained 8 threads. It can be seen that the thermal strain of each thread overcomes the influence of noise caused by physiological motion, and clear and high-signal-to-noise thermal strain distributions with differences from each other are formed in the lower right corner of the shown region.

[0100] Step 3. Perform surface fitting on the denoised thermal strain S u of each thread to obtain the heat source center position P u of each motion phase. As Figure 4 in (a) is the distribution of the heat source center positions of each thread in the interested region, Figure 4 in (b) is the normalized amplitude of the thermal strain S 1 of thread u = 1 and the fitted thermal strain surface along the x direction at different depths y, Figure 4 in (c) is the thermal strain S 1And the normalized amplitude of the thermally induced strain surface obtained by fitting along the y direction at different lateral positions x.

[0101] Step 4: Using the center position P of the thermally induced strain of thread u = 1 1 as a reference point, move the thermally induced strain images of each thread to correct the image positions in the thermally induced strain sequences of each thread; compensate for the strain at the starting moment of each thread, and transform the thermally induced strain distribution into the polar coordinate system.

[0102] Step 5: Establish a discrete grid within the region of interest, set the spatial step size Δd = 0.01 mm, and within the entire heating time of T = 20 s, set the time step size Δt = 0.1 s. Interpolate the thermally induced strain after the coordinate transformation obtained in Step 4 using cubic splines so that the thermally induced strain is sampled on the discrete grid; based on the spherical heat source assumption of the microwave ablation needle and the fact that adipose tissue is homogeneous and isotropic, use the implicit finite-difference time domain to replace the differential and solve for the quantitative relationship between the thermally induced strain and the heat source in the form of Equation (6). Specifically, the coefficient matrix A is:

[0103]

[0104] where a = κ / Δd 2 , b = ρC / Δt, B = b. Substitute the known adipose tissue parameters, such as density ρ = 937 kg / m 3 , specific heat capacity C = 3258 J / kg / K, thermal conductivity κ = 0.21 W / m / K, and the linear coefficient k = 285.71 °C between the temperature change and the thermally induced strain in the δ-ESF model into the quantitative relationship between the thermally induced strain and the heat source term in Equation (6) to estimate the heat source distribution Q u (r) of each thread, and average the heat source terms obtained for all threads to obtain the average heat source term Q(r). As Figure 5 shown, the dotted line and the solid line respectively represent the heat source distributions Q u (r) estimated for 8 threads and the heat source distribution Q(r) after averaging all threads.

[0105] Step 6: Use the same heating scheme to heat in adipose tissue, and select the coordinates of the microwave ablation needle in the ultrasonic image as the heat source position. Under the same discrete grid, based on the heat source distribution Q(r) estimated in Step 5, substitute the adipose tissue-related parameters, and combine the quantitative relationship between the heat source term and the temperature in Equation (7) to estimate the temperature distribution in real time. As Figure 6 shown, Figure 6 (a) in shows the predicted temperature distribution in adipose tissue after cumulative heating for 20 seconds, Figure 6 (b) in shows the temperature distribution along the x direction at the position of the heat source in the tissue at different heating times, Figure 6In (c), it is the temperature distribution along the y direction of the overheat source position in the tissue at different heating times. It can be observed that the temperature at the heat source position gradually increases with heating, and at the same time, the lateral and longitudinal ranges of the temperature rising area spread around due to heat conduction. The present invention does not require the aid of additional equipment, can estimate the heat source distribution in vivo, and predict the temperature field based on the heat source information, providing a temperature basis for clinical treatments, applications such as treating, assisting in treating, or diagnosing through the temperature change of biological tissues, and has high application value.

[0106] The present invention has been described in detail above in connection with specific exemplary embodiments. However, it should be understood that various modifications and variations can be made without departing from the scope of the present invention as defined by the appended claims. The detailed description and the drawings should be considered illustrative only and not restrictive. If there are any such modifications and variations, then they will all fall within the scope of the present invention described herein. In addition, the background art is intended to illustrate the research and development status and significance of the present technology and is not intended to limit the present invention or the application fields of the present application and the present invention.

Claims

1. A method for estimating the distribution of internal heat sources based on ultrasonic imaging, characterized in that, it includes the following steps: Acquisition step: Using a hyperthermia method to heat the target area, performing ultrasonic image acquisition on the in-vivo hyperthermia target area to obtain a digital ultrasonic image sequence I; Imaging step: Select the region of interest according to the target area. Through correlation analysis, image registration based on the region of interest, and thermal strain imaging, divide the image sequence I into multiple threads u corresponding to different motion phases, where u = 1, 2, …, U, and U is the total number of threads, to obtain the thermal strain image sequence S corresponding to the threads u ; Construct a denoising function based on the correlation between threads for S u to perform noise reduction and obtain the updated S u ; Fitting step: For the updated S u perform two-dimensional axisymmetric surface fitting to obtain the heat source center position P of each thread u ; Calibration steps: By moving the thermal strain image sequences S corresponding to each thread u , so that the heat source center positions P of each thread u u have the same coordinates; Averaging step: Based on the finite-difference time-domain method and the bioheat conduction equation, establishing a quantitative relationship between the heat source distribution and the thermal strain distribution, estimating the heat source distribution of each thread, and averaging the estimation results of all threads to obtain the average heat source distribution; Prediction step: Based on the position of the actual heat source, the average heat source distribution obtained in the averaging step, and the quantitative relationship between the heat source distribution and the thermal strain distribution, performing real-time prediction on the thermal strain and temperature field distributions that change with time; The specific implementation manner of the fitting step is: At any moment, denote the thermal strain distribution of thread u as S u (x, y; t), where u = 1, 2, …, U, U is the total number of threads, and x, y, and t are the horizontal, vertical, and time coordinates, respectively; The thermal strain distribution of the last frame image of thread u is S u (x, y; T u ), T u is the corresponding time coordinate; For S u (x, y; T u ) perform two-dimensional axisymmetric surface fitting to obtain the heat source center position P of each thread u ; or Take ln[S u (x, y; T u )] and perform polynomial least squares fitting according to the following function f(x, y). Here, ln is the logarithmic function with base e, f(x,y) = ax 2 + by 2 + cx + dy + e, (1) Obtain the fitting parameters a, b, c, d, e; then the vertical and horizontal coordinates of the heat source center position P u are respectively determined as The specific implementation manner of the averaging step is: For any thread, taking the heat source position as the center, establishing a discrete grid with an axial dimension A×transverse dimension L not less than the region of interest, a maximum spatial step size Δd, a time span of 0 to T, and a time step size Δt, where 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 sequence S of thermal strain after coordinate transformation u (r, θ; t) is interpolated in the spatial and temporal domains so that the sequence is sampled on the above discrete grid; based on the following discrete equation, the heat source term is obtained: AS u (r, θ; t + Δt) = BS u (r, θ; t) + Q u (r, θ; t) / k, (3) Among them, A and B are coefficient matrices related to tissue specific heat capacity, thermal conductivity, and density, which can be obtained by solving the bioheat conduction equation using the finite difference time domain method; k is the constant specified in the thermal strain temperature measurement model δ-ESF that describes the linear relationship between the temperature change and the thermal strain; the heat source distribution Q u (r, θ) of this thread is the Q u (r, θ; t) obtained by predicting at each time step. After averaging Q u (r, θ) with respect to t, by averaging the heat source terms Q (r, θ) obtained from all threads, the average heat source term Q(r, θ) is obtained.

2. A method for estimating the distribution of internal heat sources based on ultrasonic imaging according to claim 1, characterized in that, The specific implementation manner of the calibration step is: Taking any thread as a reference, moving the thermal strain images of other threads on the (x,y) plane until the heat source center coordinates of all threads are the same.

3. A method for estimating the distribution of internal heat sources based on ultrasonic imaging according to claim 1, characterized in that, The specific implementation manner of the calibration step is: The displacement field is obtained by performing cross-correlation calculation based on a sliding window on the digital ultrasonic images and the thermal strain image sequences corresponding to each thread. The updated S is moved according to the displacement field u to make the heat source center position P u of each thread u have the same coordinates.

4. A method for estimating the distribution of internal heat sources based on ultrasonic imaging according to any one of claims 1-3, characterized in that, In the acquisition step, when performing ultrasonic image acquisition on the in-vivo hyperthermia target area, collecting and storing ultrasonic image sequences of multiple respiratory cycles before starting heating, so that the first frame images of each thread occur before heating.

5. A method for estimating the distribution of internal heat sources based on ultrasonic imaging according to any one of claims 1-3, characterized in that, If the ultrasonic image acquisition starts after heating in the acquisition step, then at the end of the calibration step, it is also necessary to compensate for the strain at the starting moment of each thread, and the compensation method is: Taking the thread with the earliest starting time as the source thread, and other threads as target threads respectively, using the first and second frame thermal strain images of the source thread, performing time-domain linear interpolation at the moment when the first frame thermal strain image of the target thread is located, and using the obtained image as the compensation image; or Performing time-domain linear extrapolation on the moment when the first frame thermal strain of the source thread is located along the negative direction of t through the first and second frame thermal strain images of the target thread, and taking the negative of the obtained image as the compensation image; Add the compensated image to all the thermal strain images of the target thread, and then the thermal strain image sequence S of the target thread u is updated to the compensated result.

6. A method for estimating the distribution of internal heat sources based on ultrasonic imaging according to claim 5, characterized in that, At the end of the calibration step, the sequence S of thermal strain images u is transformed from the Cartesian coordinate system (x, y) to the polar coordinate system (r, θ).

7. A method for estimating the distribution of internal heat sources based on ultrasonic imaging according to claim 1, characterized in that, The specific implementation manner of the prediction step is: Using the same heating protocol as in the acquisition step, heat the same type of living tissue, set the actual heat source center as the origin, and iteratively calculate the thermal strain distribution S(r, θ; t) in the tissue according to the following formula: AS(r, θ; t + Δt) = BS(r, θ; t) + Q(r, θ) / k, (4) The predicted temperature distribution in the tissue is T(r, θ; t) = T 0 + kS(r, θ; t), where T 0 is the baseline temperature in the biological tissue before heating.

8. An apparatus for estimating the in-vivo heat source distribution based on ultrasonic imaging, characterized in that, it comprises: a computer, a B-ultrasound machine, a linear array probe, and a thermal ablation device; the B-ultrasound machine is connected to the linear array probe and is used for collecting ultrasonic images of a living object; the thermal ablation device is used for heating the living object to form displacement and strain distributions in the tissue; the computer is respectively connected to the B-ultrasound machine and the thermal ablation device, and is used for controlling the B-ultrasound machine and the thermal ablation device, and is configured to implement a method for estimating the in-vivo heat source distribution based on ultrasonic imaging as described in any one of claims 1-7.

Images