CN113139310B - Method, simulator and system for simulating focused ultrasound temperature field based on FDTD - Google Patents

Abstract

The present invention provides a method, a simulator, and a system for simulating a focused ultrasound temperature field based on FDTD. The method includes the following steps: S1: by obtaining a numerical simulation equation of the focused ultrasound temperature field; wherein, the focused ultrasound temperature The numerical simulation equation of the field is expressed by the Pennes biological heat conduction equation; S2: The numerical simulation equation in S1 is solved by using the time domain finite difference method with the second-order accuracy of space and the first-order accuracy of time to realize the temperature rise effect caused by the simulation of focused ultrasound. The present invention uses the time domain finite difference algorithm with the second-order accuracy of space and the first-order accuracy of time to solve the Pennes biological heat conduction equation, thereby establishing a method for simulating the temperature field of focused ultrasound based on FDTD. The numerical simulation accuracy of the focused ultrasound temperature field does correspond to the time-domain finite difference method with the second-order accuracy in space and the first-order accuracy in time, which ensures the validity and accuracy of the simulation.

Description

Method, simulator and system for simulating focused ultrasound temperature field based on FDTD

technical field

The present invention relates to the technical field of focused ultrasound temperature field simulation, and more specifically, to a method, simulator and system for simulating a focused ultrasound temperature field based on FDTD.

Background technique

High-intensity focused ultrasound (HIFU) tumor ablation technology is a new tumor treatment method that has become popular in recent years. Its biggest feature is that it uses the high-intensity ultrasound emitted by the external ultrasound transducer to pass through the skin non-invasively and form an ultrasound focus on the tumor in the body. The temperature at the focal point rises instantaneously, causing coagulation and necrosis of the surrounding tumor tissue, which eventually forms scars or is absorbed by metabolism. For example, Chinese Patent Publication No.: CN101108268A, publication date: 2008-01-23, discloses a method for forming a multi-mode thermal field of a phased array focused ultrasound in the field of biomedical engineering, including the following steps: (1) adopting a multi-element array spherical surface Annular array treatment head with unequal spacing of openings; (2) using each sound field mode to correspond to the phase and amplitude of a group of driving signals; (3) using a heuristic method to determine the sound intensity of the hyperthermia mode available in the target tissue; ( 4) The sound intensity of the thermal field in the warming mode is adjusted to ensure that the temperature rise rate of the cold spot in the heating area remains at a predetermined value, and finally achieve the purpose of forming any expected uniform warming mode thermal field.

But at present, the existing technologies all have reasons such as the waste heat of the focus or the inaccurate focus position, etc., and the surrounding healthy tissues will be easily killed, causing unnecessary damage. Therefore, the study of the temperature field caused by focused ultrasound through simulation is of great significance for the design of focused ultrasound equipment, the guarantee of the safety and effectiveness of focused ultrasound therapy, and the formulation of appropriate strategies for focused ultrasound therapy.

Contents of the invention

In order to overcome the problem of inaccurate temperature field caused by simulating focused ultrasound in the prior art, the present invention provides a method, simulator and system for simulating and simulating focused ultrasound temperature field based on FDTD, which can effectively and accurately simulate focused ultrasound Temperature Field.

In order to solve the above technical problems, the technical solution of the present invention is as follows: a method for simulating a focused ultrasound temperature field based on FDTD, the method comprising the following steps:

S1: By obtaining the numerical simulation equation of the focused ultrasound temperature field; wherein, the numerical simulation equation of the focused ultrasound temperature field is expressed by the Pennes biological heat conduction equation;

S2: The time-domain finite difference method with second-order accuracy in space and first-order accuracy in time is used to solve the numerical simulation equation in S1 to simulate the temperature rise effect caused by focused ultrasound.

Preferably, the Pennes biological heat conduction equation:

Among them, T represents the temperature of biological tissue; κ _t represents the thermal conductivity; C _t represents the specific heat of biological tissue; C _b represents the specific heat of blood vessels; W _b represents the blood perfusion rate of capillaries in biological tissue; T _a represents the initial temperature; Q _v Indicates the heat absorbed by the tissue per unit volume and unit time; Q _m indicates the heat generation rate of biological metabolism.

表示温度的升高量，忽略Q_m，一般起始温度是常数温度，将Pennes生物热传导方程表示为：Further, take

Indicates the amount of temperature rise, ignoring Q _m , the general initial temperature is a constant temperature, and the Pennes biological heat conduction equation is expressed as:

For the simulation of the focused ultrasound temperature field, _Qv is the external heat source input brought by the focused ultrasound to the biological tissue; for the single-frequency steady-state focused ultrasound field, the ultrasonic heat source _Qv is expressed as follows:

Q _v =αI _av =αρ ₀ c ₀ <V ² >=1/2×αρ ₀ c ₀ V ₀ ²

Among them, α represents the sound absorption coefficient; I _av represents the average sound intensity; c ₀ represents the wave velocity; V represents the particle vibration velocity; <*> represents the periodic average of the variables inside; V ₀ represents the particle vibration velocity amplitude.

Furthermore, the finite-difference time-domain method with second-order spatial precision specifically uses spatial grid points to divide the simulation area, and the distance between adjacent spatial grid points is 0.1λ, where λ represents the wavelength of the ultrasonic wave .

Furthermore, the specific solution of the time-domain finite difference method of the second-order spatial precision is as follows:

For a one-variable function f(x), where x is an independent variable, the second-order precision central difference format of its second-order derivative is:

Among them, the superscript ^{(2) in f (2)} | _i represents the second derivative of the unary function f(x), and the subscript i in f| _i represents the value of the independent variable corresponding to the i-th grid point (x ₀ +iΔx), x ₀ is the initial value of the independent variable, and Δx is the difference between the independent variables between adjacent grid points;

上标n表示时间坐标(t₀+nΔt)，其中t₀表示起始时间，Δt表示数值模拟的时间步长；下标i、j、k表示笛卡尔空间坐标(x₀+iΔx,y₀+jΔy,z₀+kΔz)，其中(x₀,y₀,z₀)表示起始位置，Δx、Δy、Δz分别表示x、y、z方向的空间步长。The temperature rise at different time points and grid point locations in different spaces is defined as

The superscript n represents the time coordinate (t ₀ +nΔt), where t ₀ represents the starting time, Δt represents the time step of the numerical simulation; the subscripts i, j, k represent the Cartesian space coordinates (x ₀ +iΔx,y ₀ +jΔy,z ₀ +kΔz), where (x ₀ ,y ₀ ,z ₀ ) represents the starting position, and Δx, Δy, and Δz represent the spatial steps in the x, y, and z directions, respectively.

Furthermore, the time domain finite difference method with the first-order precision of time is specifically solved as follows:

The explicit time-discrete format of equation (2) is:

的差分格式采用方程(3)，从而有：The second order spatial derivative in equation (4)

The difference format of adopts equation (3), so that:

A kind of simulator, described simulator comprises acquisition module, the time domain finite difference method solution module of space second-order precision, the time domain finite difference method solution module of time first order precision;

The obtaining module is used to obtain the Pennes biological heat conduction equation;

The finite-difference time-domain method solution module with second-order precision in space and the finite-difference time-domain method solution module with first-order time precision solve the Pennes biological heat conduction equation respectively to realize the temperature rise effect caused by simulating focused ultrasound .

Preferably, the solution realized by the finite-difference time-domain method solving module of the second-order space precision is as follows:

For a one-variable function f(x), where x is an independent variable, the second-order precision central difference format of its second-order derivative is:

Among them, the superscript ^{(2) in f (2)} | _i represents the second derivative of the unary function f(x), and the subscript i in f| _i represents the value of the independent variable corresponding to the i-th grid point (x ₀ +iΔx), x ₀ is the initial value of the independent variable, and Δx is the difference between the independent variables between adjacent grid points;

上标n表示时间坐标(t₀+nΔt)，其中t₀表示起始时间，Δt表示数值模拟的时间步长；下标i、j、k表示笛卡尔空间坐标(x₀+iΔx,y₀+jΔy,z₀+kΔz)，其中(x₀,y₀,z₀)表示起始位置，Δx、Δy、Δz分别表示x、y、z方向的空间步长。The temperature rise at different time points and grid point locations in different spaces is defined as

The superscript n represents the time coordinate (t ₀ +nΔt), where t ₀ represents the starting time, Δt represents the time step of the numerical simulation; the subscripts i, j, k represent the Cartesian space coordinates (x ₀ +iΔx,y ₀ +jΔy,z ₀ +kΔz), where (x ₀ ,y ₀ ,z ₀ ) represents the starting position, and Δx, Δy, and Δz represent the spatial steps in the x, y, and z directions, respectively.

Further, the solution implemented by the finite-difference time-domain method solution module with the first-order accuracy of time is as follows:

The explicit time-discrete format of equation (2) is:

的差分格式采用方程(3)，从而有：The second order spatial derivative in equation (4)

The difference format of adopts equation (3), so that:

A computer system, comprising a memory, a processor, and a computer program stored in the memory and operable on the processor. When the processor executes the computer program, the computer program described in any one of claims 1 to 6 is realized. steps of the method described above.

Compared with the prior art, the beneficial effects of the technical solution of the present invention are:

The invention solves the Pennes biological heat conduction equation by using the time-domain finite difference algorithm with the second-order accuracy of space and the first-order accuracy of time, thereby establishing a method for simulating the temperature field of focused ultrasound based on FDTD. The numerical simulation accuracy of the focused ultrasonic temperature field corresponds to the second-order accuracy of space and the first-order accuracy of time domain finite difference method, which improves the effectiveness and accuracy of simulation.

Description of drawings

Fig. 1 is a flow chart of the steps of the method described in this embodiment.

Fig. 2 is the variation of the 2-norm and ∞-norm of the error of the normalized temperature rise with the space step size in this embodiment.

Fig. 3 shows the variation of the 2-norm and ∞-norm of the error of the normalized temperature rise along with the time step in this embodiment.

Fig. 4 is the change of the temperature rise of the center of the region with time in this embodiment.

Detailed ways

The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, and are only for illustrative purposes and cannot understood as a limitation on this patent. Based on the embodiments of the present invention, all other embodiments obtained by persons of ordinary skill in the art without creative efforts fall within the protection scope of the present invention.

The technical solutions of the present invention will be further described below in conjunction with the accompanying drawings and embodiments.

Example 1

As shown in Figure 1, a method for simulating a focused ultrasound temperature field based on FDTD, the method includes the following steps:

S1: By obtaining the numerical simulation equation of the focused ultrasound temperature field; wherein, the numerical simulation equation of the focused ultrasound temperature field is expressed by the Pennes biological heat conduction equation;

S2: The time-domain finite difference method with second-order accuracy in space and first-order accuracy in time is used to solve the numerical simulation equation in S1 to simulate the temperature rise effect caused by focused ultrasound.

In a specific embodiment, the Pennes biological heat conduction equation:

Among them, T represents the temperature of biological tissue; κ _t represents the thermal conductivity; C _t represents the specific heat of biological tissue; C _b represents the specific heat of blood vessels; W _b represents the blood perfusion rate of capillaries in biological tissue; T _a represents the initial temperature; Q _v Indicates the heat absorbed by the tissue per unit volume and unit time; Q _m indicates the heat generation rate of biological metabolism, and Q _m can generally be ignored during calculation.

表示温度的升高量，忽略Q_m，一般起始温度是常数温度，将Pennes生物热传导方程表示为：Pick

Indicates the amount of temperature rise, ignoring Q _m , the general initial temperature is a constant temperature, and the Pennes biological heat conduction equation is expressed as:

For the simulation of the focused ultrasound temperature field, _Qv is the external heat source input brought by the focused ultrasound to the biological tissue; for the single-frequency steady-state focused ultrasound field, the ultrasonic heat source _Qv is expressed as follows:

Q _v =αI _av =αρ ₀ c ₀ <V ² >=1/2×αρ ₀ c ₀ V ₀ ²

Among them, α represents the sound absorption coefficient; I _av represents the average sound intensity; c ₀ represents the wave velocity; V represents the particle vibration velocity; <*> represents the periodic average of the variables inside; V ₀ represents the particle vibration velocity amplitude.

In a specific embodiment, the finite-difference time-domain method with second-order spatial precision specifically uses spatial grid points to divide the simulation area, and the distance between adjacent spatial grid points is 0.1λ, where λ Indicates the wavelength of ultrasound.

Furthermore, the specific solution of the time-domain finite difference method of the second-order spatial precision is as follows:

For a one-variable function f(x), the second-order precision central difference format of its second-order derivative is:

Among them, the superscript ^{(2) in f (2)} | _i represents the second derivative of the unary function f(x), and the subscript i in f| _i represents the value of the independent variable corresponding to the i-th grid point (x ₀ +iΔx), x ₀ is the initial value of the independent variable, and Δx is the difference between the independent variables between adjacent grid points;

上标n表示时间坐标(t₀+nΔt)，其中t₀表示起始时间，Δt表示数值模拟的时间步长；下标i、j、k表示笛卡尔空间坐标(x₀+iΔx,y₀+jΔy,z₀+kΔz)，其中(x₀,y₀,z₀)表示起始位置，Δx、Δy、Δz分别表示x、y、z方向的空间步长。The temperature rise at different time points and grid point locations in different spaces is defined as

The superscript n represents the time coordinate (t ₀ +nΔt), where t ₀ represents the starting time, Δt represents the time step of the numerical simulation; the subscripts i, j, k represent the Cartesian space coordinates (x ₀ +iΔx,y ₀ +jΔy,z ₀ +kΔz), where (x ₀ ,y ₀ ,z ₀ ) represents the starting position, and Δx, Δy, and Δz represent the spatial steps in the x, y, and z directions, respectively.

In a specific embodiment, the time-domain finite-difference method of the first-order precision of the time is specifically solved as follows:

The explicit time-discrete format of equation (2) is:

的差分格式采用方程(3)，从而有：The second order spatial derivative in equation (4)

The difference format of adopts equation (3), so that:

In order to confirm the reliable operation of the method described in this embodiment, carry out following verification:

1. The initial value problem of temperature rise in a limited area

则方程(2)具有如下形式：Ignoring the heat source term Q _v on the right side of the biological heat conduction equation (2) and the terms related to blood vessels

Then equation (2) has the following form:

For a cuboid finite area with side lengths a, b, and c, it is assumed that the surface temperature rise of this area remains constant:

The initial temperature rise distribution is:

Since then, the change of temperature rise with time and space is: ^[9]

The numerical simulation settings are as follows. The side lengths of cuboids are equal, a=b=c=10cm, density ρ ₀ =1000kg/m ³ , thermal conductivity κ _t =0.6W/(m·K), specific heat C _t =4180J/(kg·K),

2. Use the 2-norm and ∞-norm to characterize the numerical calculation error, as follows:

表示温升数值解，

表示方程(9)的温升解析解。in,

represents the numerical solution of temperature rise,

represents the temperature-rise analytical solution of equation (9).

Set the time step to satisfy Δt=0.0001s and keep it constant, and analyze the influence of the space step on the numerical calculation results of temperature rise. The 2-norm and ∞-norm of the temperature rise error change with the space step as shown in Figure 1, and the numerical calculation error of the temperature rise increases with the increase of the space step. The fitting curve of the error represented by the norm along with the spatial step size distribution is also drawn in Figure 1. The result of the fitting curve shows that the numerical calculation result of the temperature rise is indeed the second-order accuracy in space.

Similarly, set the spatial step size to Δx=0.001m and keep it constant, and analyze the influence of the time step size on the numerical calculation results of the temperature rise. The 2-norm and ∞-norm of the error of the normalized temperature rise change with the time step as shown in Figure 2. The numerical calculation error of the temperature rise also increases with the increase of the time step. Moreover, the result of fitting the curve also confirms that the numerical calculation result of temperature rise is indeed the first-order accuracy in time.

Set the space step and time step as Δx=0.005m and Δt=0.1s respectively, and the total calculation time is 2500s. Figure 3 shows the curve of the temperature rise of the area center changing with time calculated by numerical simulation and the curve of the temperature rise of the area center changing with time calculated theoretically by using equation (9). It can be clearly seen from the figure The fit of the two curves is very high, and the numerical simulation results are in good agreement with the analytical results.

Example 2

A kind of simulator, described simulator comprises acquisition module, the time domain finite difference method solution module of space second-order precision, the time domain finite difference method solution module of time first order precision;

The obtaining module is used to obtain the Pennes biological heat conduction equation;

The Pennes biological heat conduction equation:

Among them, T represents the temperature of biological tissue; κ _t represents the thermal conductivity; C _t represents the specific heat of biological tissue; C _b represents the specific heat of blood vessels; W _b represents the blood perfusion rate of capillaries in biological tissue; T _a represents the initial temperature; Q _v Indicates the heat absorbed by the tissue per unit volume and unit time; Q _m indicates the heat generation rate of biological metabolism, and Q _m can generally be ignored during calculation.

表示温度的升高量，忽略Q_m，一般起始温度是常数温度，将Pennes生物热传导方程表示为：Pick

Indicates the amount of temperature rise, ignoring Q _m , the general initial temperature is a constant temperature, and the Pennes biological heat conduction equation is expressed as:

For the simulation of the focused ultrasound temperature field, _Qv is the external heat source input brought by the focused ultrasound to the biological tissue; for the single-frequency steady-state focused ultrasound field, the ultrasonic heat source _Qv is expressed as follows:

Q _v =αI _av =αρ ₀ c ₀ <V ² >=1/2×αρ ₀ c ₀ V ₀ ²

Among them, α represents the sound absorption coefficient; I _av represents the average sound intensity; c ₀ represents the wave velocity; V represents the particle vibration velocity; <*> represents the periodic average of the variables inside; V ₀ represents the particle vibration velocity amplitude.

The time-domain finite-difference method solution module with second-order accuracy in space and the finite-difference time-domain method solution module with first-order accuracy in time are respectively coupled to solve the Pennes biological heat conduction equation to realize the temperature rise effect caused by simulating focused ultrasound.

Wherein, the solution realized by the time-domain finite difference method solution module of the second-order space precision is as follows:

For a one-variable function f(x), where x is an independent variable, the second-order precision central difference format of its second-order derivative is:

Among them, the superscript ^{(2) in f (2)} | _i represents the second derivative of the unary function f(x), and the subscript i in f| _i represents the value of the independent variable corresponding to the i-th grid point (x ₀ +iΔx), x ₀ is the initial value of the independent variable, and Δx is the difference between the independent variables between adjacent grid points;

上标n表示时间坐标(t₀+nΔt)，其中t₀表示起始时间，Δt表示数值模拟的时间步长；下标i、j、k表示笛卡尔空间坐标(x₀+iΔx,y₀+jΔy,z₀+kΔz)，其中(x₀,y₀,z₀)表示起始位置，Δx、Δy、Δz分别表示x、y、z方向的空间步长。The temperature rise at different time points and grid point locations in different spaces is defined as

The superscript n represents the time coordinate (t ₀ +nΔt), where t ₀ represents the starting time, Δt represents the time step of the numerical simulation; the subscripts i, j, k represent the Cartesian space coordinates (x ₀ +iΔx,y ₀ +jΔy,z ₀ +kΔz), where (x ₀ ,y ₀ ,z ₀ ) represents the starting position, and Δx, Δy, and Δz represent the spatial steps in the x, y, and z directions, respectively.

The solution realized by the finite-difference time-domain method solving module of the first-order accuracy of time is as follows:

The explicit time-discrete format of equation (2) is:

的差分格式采用方程(3)，从而有：The second order spatial derivative in equation (4)

The difference format of adopts equation (3), so that:

Example 3

A computer system includes a memory, a processor, and a computer program stored in the memory and operable on the processor. When the processor executes the computer program, the steps of the method are as follows:

S1: By obtaining the numerical simulation equation of the focused ultrasound temperature field; wherein, the numerical simulation equation of the focused ultrasound temperature field is expressed by the Pennes biological heat conduction equation;

S2: The time-domain finite difference method with second-order accuracy in space and first-order accuracy in time is used to solve the numerical simulation equation in S1 to simulate the temperature rise effect caused by focused ultrasound.

Apparently, the above-mentioned embodiments of the present invention are only examples for clearly illustrating the present invention, rather than limiting the implementation of the present invention. For those of ordinary skill in the art, other changes or changes in different forms can be made on the basis of the above description. It is not necessary and impossible to exhaustively list all the implementation manners here. All modifications, equivalent replacements and improvements made within the spirit and principles of the present invention shall be included within the protection scope of the claims of the present invention.

Claims (6)

1. A method for simulating a focused ultrasound temperature field based on FDTD is characterized in that: the method comprises the following steps:

s1: obtaining a numerical simulation equation of a focused ultrasound temperature field; wherein, the numerical simulation equation of the focused ultrasound temperature field is expressed by a Pennes biological heat conduction equation;

s2: solving a numerical simulation equation in the S1 by adopting a time domain finite difference method of spatial second-order precision and temporal first-order precision to realize the temperature rise effect caused by simulated focused ultrasound;

the Pennes biological heat conduction equation:

wherein T represents a biological tissue temperature; kappa _t Represents a thermal conductivity coefficient; c _t Represents the specific heat of the biological tissue; c _b Represents the vascular specific heat; w _b Representing the blood perfusion rate of capillaries within the biological tissue; t is _a Represents the initial temperature; q _v Represents the amount of heat absorbed by the tissue per unit volume and unit time; q _m Representing the rate of heat generation of biological metabolism;

get

Indicating the rise in temperature, ignoring Q _m Generally, the starting temperature is a constant temperature, and the Pennes biological heat transfer equation is expressed as:

for focused ultrasound temperature field simulation, a single-frequency steady-state focused ultrasound field, Q, is used _v Is represented as follows:

Q _v ＝αI _av ＝αρ ₀ c ₀ <V ² >＝1/2×αρ ₀ c ₀ V ₀ ²

wherein α represents a sound absorption coefficient; I.C. A _av Represents the average sound intensity; c. C ₀ Represents the wave velocity; v represents the particle vibration velocity;<*>the periodic average of the variables in the device is shown; v ₀ Representing the particle vibration velocity magnitude.

2. The FDTD-based simulation focused ultrasound temperature field method of claim 1, wherein: the time domain finite difference method of the spatial second-order precision specifically adopts spatial grid points to divide an analog simulation area, and the distance between the adjacent spatial grid points is 0.1 lambda, wherein lambda represents the wavelength of ultrasonic waves.

3. The FDTD-based analog simulation focused ultrasound temperature field method of claim 2, wherein: the time domain finite difference method of the spatial second-order precision concretely solves the following steps:

for a unary function f (x), where x is an argument, the second-order center-of-precision difference format for its second derivative is:

wherein f is ⁽²⁾ | _i The superscript (2) in (1) represents the second derivative of the unary function f (x), f $ _i The index i in (a) indicates the value (x) of the argument corresponding to the ith grid point ₀ +iΔx)，x ₀ Is the initial value of the argument, Δ x is the difference of the arguments between adjacent grid points;

the temperature rises at different time and different space grid point positions are respectively defined as

Superscript n denotes the time coordinate (t) ₀ + n Δ t), where t ₀ Denotes the starting time, Δ t denotes the time step of the numerical simulation; the indices i, j, k denote the cartesian space coordinates (x) ₀ +iΔx,y ₀ +jΔy,z ₀ + k Δ z) where (x) ₀ ,y ₀ ,z ₀ ) Denotes the starting position, and Δ x, Δ y, Δ z denote the spatial step in the x, y, z directions, respectively.

4. The FDTD-based simulation focused ultrasound temperature field method according to claim 3, wherein: the time domain finite difference method of the time first-order precision concretely solves the following steps:

the time discrete format is:

second spatial derivative in equation (4)

The difference format of (3) uses equation (3) such that there is:

5. a simulator based on the method of any one of claims 1 to 4, characterized in that: the simulator comprises an acquisition module, a time domain finite difference method solving module with spatial second-order precision and a time domain finite difference method solving module with temporal first-order precision;

the acquisition module is used for acquiring a Pennes biological heat conduction equation;

the module for solving the Pennes biological heat conduction equation is respectively solved by the module for solving the time domain finite difference method with the spatial second-order precision and the module for solving the time domain finite difference method with the temporal first-order precision, so that the temperature rise effect caused by the simulated focused ultrasound is realized.

6. A computer system comprising a memory, a processor, and a computer program stored on the memory and executable on the processor, wherein: the processor, when executing the computer program, performs the steps of the method according to any of claims 1 to 4.

Images