CN107684437B - Ultrasonic temperature imaging method combining echo time shifting technology and 2D heat conduction equation - Google Patents

Abstract

An ultrasonic temperature imaging method combining echo time-shifting technique and 2D heat conduction equation, including a tissue heating stage and a tissue cooling stage, in the tissue cooling stage: obtain the original radio frequency signal at the time of cauterization stop; use the echo time-shift-based ultrasonic temperature The imaging technology obtains the tissue temperature distribution map at time tk; a sliding window is designed to estimate the tissue temperature distribution at tk _{+ 1} .

Description

Ultrasonic temperature imaging method combining echo time-shifting technique and 2D heat conduction equation

technical field

The invention relates to the technical fields of medical image processing, ultrasonic temperature imaging and radio frequency ablation, in particular to monitoring tissue radio frequency ablation in combination with echo time-shift technology and 2-D heat conduction equation in radio frequency ablation surgery.

Background technique

Primary liver cancer is the fifth most common cancer worldwide, and the mortality rate ranks third among all cancers [1]. Hepatocellular carcinoma (HCC) accounts for 85% to 90% of primary liver cancer patients. Surgical resection and liver transplantation are currently recognized as the two main treatment modalities for liver cancer. However, only a small number of patients who are physically able to receive these two modalities of treatment [2]. Therefore, the development of other more applicable liver cancer treatment methods is imminent.

Radiofrequency ablation (RFA) is a minimally invasive and invasive tumor treatment method, and it is currently one of the preferred tumor cautery modes in most institutions [3-4]. When patients with early-stage hepatocellular carcinoma are not suitable for surgical resection and liver transplantation, RFA is recognized as one of the best alternatives [5-6]. During radiofrequency ablation, a radiofrequency electrode (RF electrode) is inserted into the cancerous part of the patient. The radiofrequency electrode releases a strong current pulse, and the charged ions in the tissue are affected by the current and oscillate and produce Heat causes tissue temperature to rise. When the temperature in the tissue reaches about 50 °C, it will cause protein denaturation, resulting in tissue coagulation-type necrosis, so as to achieve the purpose of tumor treatment [7-8]. One of the keys to radiofrequency ablation is how to use image-guided radiofrequency electrodes to correctly insert the patient's cancerous site. Ultrasound imaging technology has become one of the most popular image guidance methods due to its flexibility, low price, and real-time feedback [9-11].

In order to better control the cautery dose and cause minimal harm to the patient, doctors must monitor the temperature distribution of the cautery area in real time during the radiofrequency ablation process. The temperature change inside the tissue will lead to the change of the scattering medium characteristics of the tissue, thus changing some acoustic characteristics in the ultrasonic transmission process, such as ultrasonic sound speed, attenuation coefficient, inverse scattering coefficient. Using these changes in acoustic properties, ultrasound can estimate temperature changes inside the tissue. The ultrasonic temperature imaging technology based on echo time-shift is an ultrasonic temperature imaging technology that uses the ultrasonic sound velocity and the thermal expansion effect of the tissue caused by the temperature change to estimate the temperature change of the cauterized tissue. The echo time-shifting technique has received extensive attention due to its high accuracy in the temperature range of 20 to 43 °C. In the past ten years, many studies have successfully applied echo time-shift-based ultrasonic temperature imaging technology to monitor radiofrequency ablation procedures.

However, there are still many defects in the echo time-shifting technology. Previous studies have shown that when the tissue temperature exceeds 50 °C, the ultrasonic sound speed will not change significantly with the increase of temperature. For adipose tissue, the ultrasonic sound velocity will decrease with the increase of tissue temperature. The reason for this phenomenon is that the coagulation-type necrosis of the tissue brings irreversible changes to the acoustic properties of the tissue, resulting in echo time shift in The temperature change inside the tissue cannot be effectively estimated after tissue necrosis. In the cauterization and cooling stage when the tissue has been necrotic, the echo time-shift technique has a big defect.

references:

1.Parkin D M.Global cancer statistics in the year 2000.[J].LancetOncology,2001, 2(9):533-543.

2. El-Serag H B, Rudolph K L. Hepatocellular carcinoma: epidemiology and molecular carcinogenesis[J]. Gastroenterology, 2007, 132(7): 2557-2576.

3.Lencioni R,Crocetti L.Radiofrequency ablation of liver cancer[J].Techniques in Vascular&Interventional Radiology,2007,10(1):38-46.

4. Bosch F X, Ribes J, Díaz M, et al. Primary liver cancer: Worldwide incidence and trends[J]. Gastroenterology, 2004, 127(5Suppl 1):S5-S16.

5. Sokka S, King R, Mcdannold N, et al. Gas-bubble enhanced heating inrabbit thigh in vivo[C]. Ultrasonics Symposium, 2002. Proceedings. IEEE, 2002:vol.2 1415-1418.

6. Varghese T, Techavipoo U, Zagzebski J A, et al. Impact of gas bubbles generated during interstitial ablation on elastographic depiction of in vitro thermal lesions[J]. Journal of Ultrasound in Medicine Official Journal of the American Institute of Ultrasound in Medicine, 2004, 23 (4): 545-546.

7. Goldberg S N. Radiofrequency Tumor Ablation: Principles and Techniques[J]. European Journal of Ultrasound Official Journal of the European Federation of Societies for Ultrasound in Medicine & Biology, 2001, 13(2): 129-147.

8. Jiao L R, Hansen P D, Havlik R, et al.Clinical short-term results of radiofrequency ablation in primary and secondary liver tumors[J].American Journal of Surgery,1999,177(4): 303-306.

9.Machi J, Oishi A J, Mossing A J, et al.Hand-assisted laparoscopicultrasound-guided radiofrequency thermal ablation of liver tumors:a technicalreport[J].Surgical Laparoscopy Endoscopy&Percutaneous Techniques,2002,12(3):160-164.

10. Chiou S Y, Liu J B, Needleman L. Current status of sonographicallyguided radiofrequency ablation techniques[J]. Journal of Ultrasound in Medicine Official Journal of the American Institute of Ultrasound in Medicine, 2007, 26(4): 487-499.

11.Solbiati L,Ierace T,Tonolini M,et al.Guidance and monitoring of radiofrequency liver tumor ablation with contrast-enhanced ultrasound[J].European Journal of Radiology,2004,51 Suppl(Suppl):S19-S23.

SUMMARY OF THE INVENTION

In order to better monitor the temperature distribution of the tissue burning area in the entire radiofrequency ablation process, the present invention provides an ultrasonic temperature imaging method. In the present invention, the echo time-shift technique is used to estimate the temperature distribution of the tissue in the tissue heating stage, and in the cooling stage when the cautery is stopped, the heat conduction equation in physics is used to estimate the temperature distribution of the tissue in the cooling process on the basis of the existing ultrasonic temperature image. Tissue temperature changes. The technical solution is as follows:

An ultrasonic temperature imaging method combining echo time-shifting technique and 2D heat conduction equation, including tissue heating stage and tissue cooling stage:

During the tissue warming phase:

A. Before the cauterization, acquire the ultrasonic radio frequency signal r(x, z, t ₀ ), and the ultrasonic sound speed at this time is c ₀ =1540m/s;

B. Start burning;

累加的回波信号时间偏移图可以通过下式获得C. At the time point t=t _i , acquire the ultrasonic radio frequency signal r(x, z, t _i ), and perform one-dimensional cross-correlation analysis on the radio frequency signals of the two time points t _i-1 and t _i to obtain the radio frequency signal at Time offset of two time points

The accumulated echo signal time offset map can be obtained by the following formula

D. In order to reduce the ripple effect in the temperature estimation, use an 8th-order low-pass Butterworth filter (Butterworth) to perform axial and lateral filtering on the obtained echo signal time shift map;

E. According to the following formula:

in

为时间偏移对轴向的微分，α(z)为是热膨胀系数，β(z)代表超声波声速变化与组织温度变化之间的线性关系；然后，将D步骤的计算结果对时间偏移图进行轴向微分，再乘上系数k和

即可得到组织温度分布图，大小为m×n；where T(z) is the temperature, c ₀ (z) and T ₀ are the initial sound velocity and initial temperature,

is the differential of the time migration to the axial direction, α(z) is the thermal expansion coefficient, β(z) represents the linear relationship between the ultrasonic sound velocity change and the tissue temperature change; then, the calculation results of the D step are plotted against the time migration diagram Differentiate axially, multiply by the coefficient k and

The tissue temperature distribution map can be obtained, and the size is m×n;

During the tissue cooling phase:

A. First obtain the original radio frequency signal at the cauterization stop time, denoted as r(x, z, t _k );

B. Obtain the tissue temperature distribution map at time t _k using the ultrasonic temperature imaging technology based on echo time shift;

C. According to the following formula:

代表位置(x_i,y_j)在时刻t_k的组织温度，l＝α²(t_k+1-t_k)，

为热扩散系数，k为热传导系数，c为介质比热容，ρ为介质密度，取α²＝0.586W/m·k；对于位置为(x_i,y_j)的像素点在时刻t_k+1的温度

可以通过该点以及该点周围的点在时刻t_k的温度u^k估计出：设计一个滑动窗

使用该滑动窗扫描整个图像矩阵，估计时刻t_k+1的温度图像中每一个像素点的温度值，另外，由于图像边缘的像素点使用滑动窗的时候会超出图像矩阵的索引范围，使用室内温度值T₀对图像矩阵进行扩充，使图像的大小变为(m+2)×(n+2)；in

represents the tissue temperature of the position (x _i , y _j ) at time t _k , l=α ² (t _k+1 −t _k ),

is the thermal diffusivity, k is the thermal conductivity, c is the specific heat capacity of the medium, ρ is the density of the medium, take α ² =0.586W/m·k; for the pixel at the position of (x _i , y _j ) at time t _k+1 temperature

It can be estimated by the temperature u ^k of this point and the points around this point at time t _k : Design a sliding window

Use the sliding window to scan the entire image matrix, and estimate the temperature value of each pixel in the temperature image at time t _k+1 . In addition, since the pixels at the edge of the image will exceed the index range of the image matrix when using the sliding window, use indoor The temperature value T ₀ expands the image matrix so that the size of the image becomes (m+2)×(n+2);

D. After step C, the tissue temperature distribution at time t _k+1 can be obtained; then step B is repeated for the image to obtain the tissue temperature distribution image at the next time.

Description of drawings

Fig. 1 is the flow chart of the ultrasonic temperature imaging algorithm estimated using the echo time-shift technique in the heating stage.

Figure 2 is a flow chart of estimating the temperature distribution in the cooling process using the two-dimensional heat conduction equation.

Figure 3 is: the experimental equipment architecture diagram.

Detailed ways

The present invention will be described in several aspects below.

Ultrasound Thermography

The effect of temperature change on ultrasonic echoes was proposed and widely discussed as early as the 1880s. The main effects include the intensity of the echo signal and the time offset of the echo signal. In this study, we mainly discuss the principles and assumptions of this technique. In 1996, Roberto Maass-Moreno believed that the time shift of ultrasonic echo caused by temperature change mainly comes from two physical changes: one is the change of ultrasonic sound speed caused by temperature change; the other is the thermal expansion of tissue caused by temperature change. According to the analysis, the relationship between the temperature change δT and the time delay δt can be described as:

where z is the axial depth along the ultrasonic wave propagation direction, α(z) is the thermal expansion coefficient, β(z) represents the linear relationship between the change in the ultrasonic sound velocity and the change in tissue temperature, and c ₀ (z) is dependent on the Ultrasonic sound velocity at tissue depth z. To simplify the analysis, the following assumptions are made: (i) Tissue properties do not change as tissue temperature increases. (ii) The thermal expansion coefficient α is the same inside the tissue. The above two assumptions are satisfied in the temperature range in this study, so the nonlinear effect caused by thermal expansion can be ignored, and Equation (1) is simplified as:

in

为时间偏移对轴向的微分，k为两种效应相互作用的结果。从公式(2)和(3)可以看出，通过获取超声波回波信号时间偏移并对轴向微分，再乘以系数k，就可以得到组织的温度变化，而信号的时间偏移可以通过互相关分析法获得。where T(z) is the temperature, c ₀ (z) and T ₀ are the initial sound velocity and initial temperature,

is the differential of the time offset with respect to the axial direction, and k is the result of the interaction of the two effects. It can be seen from formulas (2) and (3) that the temperature change of the tissue can be obtained by obtaining the time offset of the ultrasonic echo signal and differentiating it in the axial direction, and then multiplying it by the coefficient k, and the time offset of the signal can be obtained by obtained by cross-correlation analysis.

The main steps for estimating tissue temperature using echo time-shifting techniques are as follows:

1. First, before cauterization, acquire the ultrasonic radio frequency signal r(x,z,t ₀ ). The ultrasonic sound velocity at this time is taken as c ₀ =1540 m/s.

2. Start cauterizing.

累加的回波信号时间偏移图可以通过下式获得3. At time point t=t _i , acquire the ultrasonic radio frequency signal r(x,z,t _i ). Perform one-dimensional cross-correlation analysis on the RF signals at two time points t _i-1 and t _i to obtain the time offset of the RF signal at the two time points

The accumulated echo signal time offset map can be obtained by the following formula

4. In order to reduce the ripple effect in the temperature estimation, we use an 8th-order low-pass Butterworth filter (Butterworth) to perform axial and lateral filtering on the obtained echo signal time-shift map.

即可得到组织温度分布图。5. Finally, according to the theoretical model proposed before, the time migration map is axially differentiated, and then multiplied by the coefficient k and

The tissue temperature distribution map can be obtained.

There are still many inherent defects in the echo time-shift method. During the cooling phase at the end of cauterization, the acoustic properties of the tissue are also changed due to coagulation-type necrosis of the tissue. In this case, it will be very difficult to continue to use the cross-correlation algorithm to analyze the time shift of the RF signal before and after the temperature change. More seriously, the cross-correlation algorithm is very sensitive to tissue movement and air bubbles in the tissue, and the large number of air bubbles generated by cauterization will cause large errors in the calculated time offset values. Therefore, it is no longer a good choice to continue to use the echo time-shift method to estimate the ultrasonic temperature image in the cooling stage.

2-D heat conduction equation to estimate temperature distribution

At present, a large number of studies have discussed and solved the two-dimensional heat conduction equation, and here we only discuss the application of the two-dimensional heat conduction equation in estimating ultrasonic temperature images. When the cauterization is stopped, the temperature distribution map of the tissue at this time is estimated by using the echo time-shift technique, assuming that the temperature is the maximum temperature during the cauterization process. A rectangular coordinate system is established with the left vertex of the image as the origin, the direction along the ultrasonic probe is the vertical axis z, and the direction perpendicular to the ultrasonic probe is the horizontal axis x. u(x,z,t) represents the temperature of the point (x,z) at time t. According to Fourier's experimental law in heat transfer, the total heat that flows into this surface from time t ₁ to t ₂ is equal to the heat that should be absorbed in the same time interval. Since the RFA is over, there is no heat source in the tissue, so we can get the homogeneous equation for heat conduction:

为热扩散系数，k为热传导系数，c为介质比热容，ρ为介质密度。in,

is the thermal diffusivity, k is the thermal conductivity, c is the specific heat capacity of the medium, and ρ is the density of the medium.

The basic idea of solving partial differential equations by the finite difference method is to convert the continuous distribution function in the solution area into discrete function values in the discrete grid through discretization, and use finite difference equations to approximate the continuous differential equations, so as to obtain Approximate solutions to partial differential equations. In the RFA cooling stage, assuming that the room temperature is constant, the temperature distribution map at the end of heating is obtained through the echo time-shift technique. Each pixel in the temperature distribution map represents the tissue temperature value at that location, so the number of discrete equations is limited. . As the cooling process progresses, the tissue temperature will infinitely approach room temperature. The initial conditions and termination conditions are known, so the difference equation must have a solution. The finite difference format of the heat conduction equation is:

代表位置(x_i,y_j)在时刻t_k的组织温度。x_i+1-x_i＝z_j+1-z_j＝1，则公式(6)可以化简为：in

represents the tissue temperature at time t _k at the location ( _xi , y _j ). x _i+1 -x _i =z _j+1 -z _j =1, then formula (6) can be simplified as:

where m=α ² (t _k+1 −t _k ). Since the human body is mostly composed of water, α ² =0.586W/m·k.

The present invention uses the two-dimensional heat conduction equation to estimate the RFA cooling process, which can be divided into the following steps:

A. First obtain the original radio frequency signal at the cauterization stop time, denoted as r(x, z, t _k ).

B. The tissue temperature distribution map at time t _k is obtained using the echo time-shift-based ultrasonic temperature imaging technique, with a size of m×n.

可以通过该点以及该点周围的点在时刻t_k的温度u^k估计出。因此我们设计了一个滑动窗

使用该滑动窗扫描整个图像矩阵，估计时刻t_k+1的温度图像中每一个像素点的温度值。另外，由于图像边缘的像素点使用滑动窗的时候会超出图像矩阵的索引范围，因此我们使用室内温度值T₀对图像矩阵进行扩充，使图像的大小变为(m+2)×(n+2)。C. According to formula (7), consider a point with coordinates (x _i , y _j ), the temperature of this point at time t _k+1

It can be estimated from the temperature u ^k of this point and the points surrounding it at time t _k . So we designed a sliding window

Use the sliding window to scan the entire image matrix, and estimate the temperature value of each pixel in the temperature image at time tk ₊₁ . In addition, since the pixels at the edge of the image will exceed the index range of the image matrix when using the sliding window, we use the indoor temperature value T ₀ to expand the image matrix, so that the size of the image becomes (m+2)×(n+ 2).

D. After step C, the tissue temperature distribution at time t _k+1 can be obtained. Then repeat step B for that image.

Before radiofrequency ablation, the porcine liver sample was cut to an appropriate size, placed in a plastic box filled with normal saline, and then the ablation electrode was inserted into the isolated liver through a small hole. Clay material is used to avoid leakage of brine solution. The ultrasound probe is placed in the liver and immersed in saline solution, and the distance between the transducer and the sample depends on the focal length of the transducer, which is adjustable. Thus, the sample can be located in the focal region of the ultrasound scan. The ultrasound system is then turned on, and the location of the electrodes can be found. In the process of radiofrequency ablation, the transducer of the ultrasound system continuously captures the backscattered signal scattered by the pig liver tissue, and after capturing it by the transducer, the software system that comes with the ultrasound automatically saves it as binary data. Then we can get ultrasound images by processing these binary data.

The RFA system operates in the default automatic mode, starting at 50W/min and then automatically increasing by 10W/min, because the output of the high impedance and RF systems is the same for different electrode lengths. During the heating phase (12 min), raw RF data were obtained separately from the tissue, including 256 scanning ultrasound inverse scatter signal lines. The sampling frequency and pulse length were set to 30 MHz and 0.7 mm. After the heating phase, the radiofrequency ablation system automatically stops working. An additional experiment was then performed for each electrode length (0.5, 1, and 1.5 cm) and 5 pig liver samples (n=15).

The obtained ultrasonic data is read into the Matlab program, and then processed according to the algorithm described above, and then the final ultrasonic image is obtained, which is then compared with the actual burning situation during the test, and then the feasibility of this algorithm is verified. .

Claims (1)

1. An ultrasonic temperature imaging method combining an echo time shift technology and a 2D heat conduction equation is used for imaging in vitro pork liver tissues and comprises a tissue heating stage and a tissue cooling stage:

in the tissue warming stage:

A. prior to cauterization, an ultrasonic radio frequency signal r (x, z, t) is acquired₀) At this time, the ultrasonic sound velocity is c₀＝1540m/s；

B. Starting to burn;

C. at time t ═ t_iAcquiring an ultrasonic radio frequency signal r (x, z, t)_i) For two time points t_i-1And t_iThe radio frequency signals are subjected to one-dimensional cross-correlation analysis to obtain the time offset of the radio frequency signals at two time points

The time-shift map of the accumulated echo signal is obtained by

D. In order to reduce ripple effect appearing in temperature estimation, axial and lateral filtering is carried out on the obtained echo signal time offset graph by using an 8-order low-pass Butterworth filter;

E. according to the following formula:

wherein

Wherein T (z) is temperature, c₀(z) is the ultrasonic sound velocity dependent on the tissue depth z,

α (z) is thermal expansion coefficient for differentiating the time deviation to the axial direction, β (z) represents the linear relation between the ultrasonic sound velocity change and the tissue temperature change, the calculation result of the step D is subjected to axial differentiation on the time deviation graph, and then the coefficients k and k are multiplied

Obtaining a tissue temperature distribution map with the size of m × n;

in the tissue cooling stage:

A. firstly, the original radio frequency signal at the moment of cauterization stop is obtained and recorded as r (x, z, t)_h)；

B. Obtaining a time t using an echo time-shift based ultrasound thermography technique_hTissue temperature profile of;

C. according to the following formula:

wherein

Representative position (x)_i,y_j) At time t_hTissue temperature of (l) ═ α²(t_h+1-t_h)，

Q is a thermal diffusivity, q is a thermal conductivity, c is a specific heat capacity of a medium, ρ is a density of the medium, and q is 0.586, unit: W/m.K, for the position (x)_i,y_j) At time t_h+1Temperature of

Can pass through the point and points around the point at the time t_hTemperature u of^hEstimating that: design a sliding window

Using the sliding window to scan the entire image matrix, the time t is estimated_h+1In addition, because the pixel points at the edge of the image exceed the index range of the image matrix when using the sliding window, the indoor temperature value T is used₀Expanding the image matrix so that the size of the image becomes (m +2) × (n + 2);

D. the moment t can be obtained through the step C of the tissue cooling stage_h+1Tissue temperature distribution of (a); and then repeating the step B of the tissue cooling stage on the image to obtain the tissue temperature distribution image at the next moment.

Images