Microwave thermoacoustic imaging multi-physical field simulation detection and monitoring method

Abstract

The invention relates to a microwave thermoacoustic imaging multi-physical field simulation detection and monitoring method, which is used for establishing a complete thermoacoustic multi-physical field simulation detection and monitoring process. On the basis of COMSOL, electromagnetic wave, frequency domain, biological heat transfer, solid mechanics and pressure acoustics and transient physical field interfaces are sequentially coupled to achieve visualization of the whole microwave thermoacoustic imaging physical process, and thermoacoustic image reconstruction is conducted by calculating thermoacoustic signals detected by each ultrasonic transducer; on the basis, a steady-state solver and parametric scanning are combined to realize visualization of electromagnetic energy absorption of the biological tissue along with time change, and a plurality of transient solvers are adopted to calculate a thermoacoustic signal at each monitoring time point. According to the invention, microwave excitation and energy deposition, thermoelastic expansion and ultrasonic propagation processes related to microwave thermoacoustic imaging can be described, which has important application values for deep research of a microwave thermoacoustic imaging mechanism and guidance of biological tissue imaging detection and monitoring experiments.

Description

Technical Field

[0001] The present invention relates to the field of microwave thermoacoustic imaging, and in particular to a method for microwave thermoacoustic imaging multi-physical field simulation detection and monitoring. Background Technique

[0002] As a new hybrid imaging mode that combines the advantages of microwave imaging and ultrasonic imaging, microwave thermoacoustic imaging technology has the characteristics of high contrast, deep penetration depth, and high resolution, and can reflect the dielectric function information and structural information of biological tissues simultaneously. With the formal entry of microwave thermoacoustic imaging technology into the field of clinical disease detection, ensuring the reliability and safety of the imaging process has become a huge challenge faced by current researchers. On the one hand, the above challenges can be overcome through a large number of repetitive bio-tissue phantom or in-vivo imaging experimental studies, but this undoubtedly consumes a large amount of research costs; on the other hand, numerical simulation has been widely applied to the research in multiple fields. Through numerical simulation, not only can the costs and risks of experimental research be reduced, but also the progress of experiments can be guided to a certain extent.

[0003] The numerical simulation software used in the field of microwave thermoacoustic imaging mainly includes Matlab, CST, ANSYS, and COMSOL Multiphysics, etc. Among them, due to its excellent integration and modularity, simple and fast operation, and significant advantages in multi-physical field coupling research, more and more researchers tend to use COMSOL Multiphysics software for relevant simulation research on microwave thermoacoustic imaging. Among them, Xi Lei et al. selected the "Electromagnetic Waves, Frequency Domain" physical field interface to explore the artifact problem in microwave thermoacoustic imaging and compared the internal SAR distributions of different morphological feature models. T. George et al. performed microwave thermoacoustic imaging of a two-dimensional breast tumor model and coupled the "Electromagnetic Waves, Frequency Domain" and "Bioheat Transfer" physical field interfaces to evaluate the minimum power required to generate thermoacoustic signals by microwave sources with different frequencies. Y. Deng et al. proposed a two-dimensional numerical simulation model for simulating the generation of near-field microwave thermoacoustic signals and explored the electromagnetic interaction between microwaves and imaging targets and the generation and propagation of acoustic signals through the coupling between the "Electromagnetic Waves, Frequency Domain" and "Pressure Acoustics, Transient" physical field interfaces. M. Soltani et al. studied the acoustic signal intensity generated by breast cancer tissue under the action of 2.45 GHz continuous pulsed waves, and then further explored the influence of the size and shape of tumors in breast tissue on microwave thermoacoustic signals; Huang Lin et al. studied the influence of the Young's modulus of breast cancer tissue on microwave thermoacoustic imaging and explored the relationship between the Young's modulus of tissue and the thermal expansion rate. The relevant results of this study contribute to exploring the possibility of microwave thermoacoustic elastography.

[0004] The above-mentioned related simulation studies based on COMSOL Multiphysics have played a certain role in promoting the development of microwave thermoacoustic imaging technology. However, the selected physical field interfaces cannot completely describe the microwave thermoacoustic imaging process. Therefore, constructing a complete multi-physics coupling model that can accurately describe microwave excitation, energy deposition, thermoelastic expansion, and ultrasonic propagation involved in microwave thermoacoustic imaging has important application value for both in-depth research of basic theories and experimental guidance for biological tissue imaging detection and monitoring.

[0005] After retrieval, the application publication number is CN110147633B, which is a simulation method for ultrasonic elastography based on multi-physics coupling. It uses the coupled simulation of the acoustic field and solid mechanics field in the COMSOL Multiphysics multi-physics simulation software. First, the acoustic field simulation is completed using the acoustic module, and the radio frequency echo signal RF1 is extracted. The acoustic radiation force is calculated through the acoustic characteristic parameters in the solution results of the acoustic field simulation and applied to the solid mechanics field to obtain the responses such as the displacement generated when the acoustic radiation force acts on the medium. The information such as the displacement and shear wave velocity in the solid mechanics field is analyzed to complete the forward problem simulation of the ultrasonic elastography method. Then, the deformed geometry is reapplied to the acoustic module for solution to obtain another set of radio frequency echo signals RF2 after the object field deforms. The two sets of echo signals are used for the solution of the inverse problem.

[0006] The principle of this patent is different from that of this patent. The principle of the simulation method for ultrasonic elastography based on multi-physics coupling is that due to effects such as absorption, scattering, and reflection during the propagation of ultrasonic waves, the energy density changes, thereby generating an acoustic radiation force. The acoustic radiation force acts on a medium with elastic properties, resulting in axial compression and tension and transverse shear changes. By calculating the axial displacement or monitoring the propagation of shear waves, the elastic properties of the measured medium are evaluated, and the elastic distribution of the medium is reconstructed. The principle of the multi-physics simulation detection and monitoring method for microwave thermoacoustic imaging in this patent is based on the microwave thermoacoustic effect, that is, the medium generates thermoelastic expansion after absorbing electromagnetic energy and then radiates ultrasonic waves outward; by detecting the ultrasonic waves generated by the medium and processing them through a specific image reconstruction algorithm, a thermoacoustic image containing the dielectric function information and structural information of the medium can be obtained.

[0007] This patent has a different coupling method for the physical field interface from this patent. The ultrasonic elastography simulation method based on multi - physical field coupling uses the COMSOL Multiphysics simulation software to couple the "Pressure Acoustics, Frequency Domain" and "Pressure Acoustics, Transient" physical field interfaces to calculate the acoustic radiation force generated by ultrasonic waves and the radio - frequency echo signal RF1 carrying the initial position information of the medium; it couples the "Solid Mechanics" and "Pressure Acoustics, Transient" physical field interfaces to calculate the local displacement of the medium caused by the acoustic radiation force and the radio - frequency echo signal RF2 carrying the position information of the medium after displacement. While this patent, a multi - physical field simulation detection and monitoring method for microwave thermoacoustic imaging, realizes the visualization of the entire microwave thermoacoustic effect based on the COMSOL Multiphysics simulation software by successively coupling "Electromagnetic Waves, Frequency Domain", "Bioheat Transfer", "Solid Mechanics" and "Pressure Acoustics, Transient". That is, it couples the "Electromagnetic Waves, Frequency Domain" and "Bioheat Transfer" physical field interfaces to calculate the temperature rise generated by the biological tissue absorbing microwave energy; it couples the "Bioheat Transfer" and "Solid Mechanics" physical field interfaces to calculate the thermo - elastic expansion of the biological tissue due to the temperature increase; it couples the "Solid Mechanics" and "Pressure Acoustics, Transient" physical field interfaces to calculate the ultrasonic waves generated by the mutual extrusion of the biological tissue and the propagation of ultrasonic waves in the surrounding medium. Summary of the Invention

[0008] The present invention aims to solve the above - mentioned problems of the prior art. A multi - physical field simulation detection and monitoring method for microwave thermoacoustic imaging is proposed. The technical solution of the present invention is as follows:

[0009] A multi - physical field simulation detection and monitoring method for microwave thermoacoustic imaging, which comprises the following steps:

[0010] S1. Construct a microwave thermoacoustic multi - physical simulation geometric model and set a microwave excitation function;

[0011] S2. According to the thermoacoustic effect, set the "Electromagnetic Waves, Frequency Domain", "Bioheat Transfer", "Solid Mechanics" and "Pressure Acoustics, Transient" physical field interfaces and their corresponding boundary conditions;

[0012] S3. According to the selected physical field interfaces, set the material parameters of the model;

[0013] S4. Divide the grid according to the simulation requirements;

[0014] S5. Add a steady - state solver and a transient solver for subsequent physical field result calculations;

[0015] S6. Obtain the visualization results of microwave excitation and energy deposition, thermo - elastic expansion and ultrasonic propagation, calculate the thermoacoustic signals detected by each ultrasonic transducer and complete the reconstruction of the thermoacoustic image.

[0016] Further, in step S1, constructing a microwave thermoacoustic multi-physics simulation geometric model and setting a microwave excitation function specifically includes:

[0017] The multi-physics simulation geometric model is set according to the prior information of the structure of the biological tissue to be measured;

[0018] The microwave excitation function an1(t) is as follows:

[0019] an1(t) = rect1(t - t 0 )

[0020] where rect1 represents a rectangular wave function with a pulse width of 550 ns, and t 0 represents the rising edge time of the rectangular pulse;

[0021] Further, in step S2, the physical field interfaces of "electromagnetic wave, frequency domain", "bioheat transfer", "solid mechanics", and "pressure acoustics, transient" are sequentially coupled for solving the microwave thermoacoustic multi-physics process; among them,

[0022] The physical field interfaces of "electromagnetic wave, frequency domain" and "bioheat transfer" are used to solve the temperature rise of the biological tissue due to the absorption of electromagnetic energy; the physical field interfaces of "bioheat transfer" and "solid mechanics" are used to solve the thermoelastic expansion of the biological tissue due to the temperature rise; the "solid mechanics" and "pressure acoustics, transient" physical fields are used to solve the thermoacoustic signal generated by the thermoelastic expansion of the biological tissue.

[0023] Further, in step S2, the boundary conditions include electromagnetic boundary conditions, thermodynamic boundary conditions, and acoustic boundary conditions; among them, the electromagnetic boundary conditions include port boundary conditions, perfect electric conductor boundary conditions, and impedance boundary conditions; among them,

[0024] The input power of the port boundary condition is 1000 kW / m, the port type is set to rectangular, the port mode type is set to transverse electric, and the mode number is selected as 1;

[0025] The perfect electric conductor boundary condition is:

[0026]

[0027] The impedance boundary condition is:

[0028]

[0029] where, represents the unit normal vector, represents the electric field at the boundary of the coupling agent and transformer oil, μ represents the magnetic permeability of the coupling agent and transformer oil, ε represents the dielectric constant of the coupling agent and transformer oil, σ represents the conductivity of the coupling agent and transformer oil, ω represents the microwave angular frequency, represents the magnetic field at the boundary of the coupling agent and transformer oil, represents the source electric field;

[0030] The thermodynamic boundary conditions include: heat flux boundary condition and low reflection boundary condition; where,

[0031] The heat flux boundary condition:

[0032]

[0033] The low reflection boundary condition:

[0034]

[0035] where, represents the conductive heat flux of biological tissue, h represents the heat transfer coefficient between biological tissue and the coupling agent and transformer oil, which is set to 5 W / (m 2 K), where, W represents the power unit watt, m represents the distance unit meter, K represents the temperature unit Kelvin; T ext represents the temperature of the coupling agent and transformer oil, which is set to 310.15 K, T represents the temperature of biological tissue, σ s represents the Cauchy stress tensor of biological tissue, represents the mechanical impedance of biological tissue, represents the displacement occurred in biological tissue;

[0036] The acoustic boundary conditions include: acoustic-structure boundary condition and plane wave radiation boundary condition; where,

[0037] The acoustic-structure boundary condition:

[0038]

[0039] The plane wave radiation boundary condition:

[0040]

[0041] where, p represents the sound pressure generated by the thermal expansion of biological tissue, ρ represents the density of the coupling agent and transformer oil, which is set to 860 kg / m 3 , v s represents the sound speed in the coupling agent and transformer oil, which is set to 1420 m / s.

[0042] Furthermore, in step S3, according to the selection of the physical field interface, the material parameters of the model are set; specifically including:

[0043] The physical field interface of "electromagnetic wave, frequency domain" needs to define the permittivity, permeability and conductivity of the geometric model;

[0044] The physical field interface of "bioheat transfer" needs to define the constant pressure heat capacity, thermal conductivity and density of biological tissue;

[0045] The physical field interface of "solid mechanics" needs to define the thermal expansion coefficient, bulk modulus and shear modulus of biological tissue;

[0046] The "pressure acoustics, transient" needs to define the sound speed and density of the coupling agent transformer oil.

[0047] Furthermore, in step S4, the maximum mesh size h max needs to satisfy the following conditions:

[0048]

[0049] where f max represents the maximum ultrasonic frequency to be solved by the physical field of "pressure acoustics, transient". Considering that the microwave pulse width is set to 550 ns in the simulation, the high-frequency component of the frequency of the thermoacoustic signal excited does not exceed about 2 MHz. Therefore, f max is set to 2 MHz.

[0050] Furthermore, step S5 specifically includes the following steps:

[0051] S501. Multi-physics simulation solution settings;

[0052] S502. Multi-physics monitoring solution settings;

[0053] Among them, in step S501, a steady-state solver is used to calculate the physical field interface of "electromagnetic wave, frequency domain", and the center frequency is set to 3 GHz; three transient solvers are used to solve the physical field interfaces of "bioheat transfer", "solid mechanics" and "pressure acoustics, transient"; among them, for the physical field interface of "bioheat transfer", the "backward difference formula" method is selected in the solver configuration, the time step is set to 0.1 μs, and the solution time is set to 100 μs; for the physical field interface of "solid mechanics", the "generalized α" method is selected in the solver configuration, the time step is set to 0.1 μs, and the solution time is set to 100 μs; for the physical field interface of "pressure acoustics, transient", the "generalized α" method is selected in the solver configuration, and the time step is set to The solution time is set to 100 μs;

[0054] In step S502, the difference in the solver used for physical field solution and the specific configuration of the solver lies in that for the physical field interface of "electromagnetic wave, frequency domain", in addition to adding a steady-state solver, it is also necessary to additionally define a variable as the monitoring time parameter and add parametric scanning to calculate the electromagnetic energy absorption distribution corresponding to any time; for the physical field interfaces of "bioheat transfer", "solid mechanics", and "pressure acoustics, transient", every 3 transient solvers are used to calculate the thermoacoustic signals generated at a single monitoring moment. Then, the number of transient solvers required to monitor the thermoacoustic signals generated at multiple moments is:

[0055] n = 3n t

[0056] where n represents the number of transient solvers, and n t represents the number of monitoring moments.

[0057] Furthermore, the "backward difference formula" method is a numerical calculation method used to approximate the derivative of a function, specifically expressed as:

[0058]

[0059] where ▽ represents the difference operator. For the solution of the physical field interface of "bioheat transfer", f(t i ) represents the temperature of biological tissue at the t i -th solution moment, m represents the maximum BDF order, which does not exceed 4 at most, is the combination number, representing the temperature coefficient of biological tissue at the t i-j -th solution moment;

[0060] The "generalized α" method is a numerical calculation method for solving second-order differential equations, specifically expressed as:

[0061]

[0062] where α f , α m , γ and β represent correction coefficients used to control the accuracy and stability of equation solution, which are derived from the numerical dissipation coefficient ρ ∞ ; and are set to 0.75 in the solutions of "solid mechanics" and "pressure acoustics, transient" physical fields.

[0063] Furthermore, in step S6, the visualization results of microwave excitation, energy deposition, thermoelastic expansion, and ultrasonic propagation are obtained as follows:

[0064] After simulation calculation, add a 2D plotting group to the results and select the surface to analyze the specific absorption rate distribution, temperature distribution, displacement distribution, and acoustic field distribution of the simulation; add a 1D plotting group and select the point result graph to analyze the changes in temperature, displacement, and pressure over time at a certain point within the biological tissue.

[0065] The calculation of the photoacoustic signals detected by each ultrasonic transducer and the completion of photoacoustic image reconstruction are specifically as follows:

[0066] After simulation calculation, add a 2D cut point to the dataset, input the coordinate positions of each ultrasonic transducer, and change the dataset to "Pressure Acoustics, Transient"; on this basis, add a point calculation to the derived values, change the dataset to the 2D cut point, and change the expression to actd.p_t, then the photoacoustic signals detected by each ultrasonic transducer can be calculated; photoacoustic image reconstruction is performed using the delay superposition algorithm based on the time-domain information and spatial position relationship of the photoacoustic signals detected by each ultrasonic transducer.

[0067] Furthermore, in the microwave photoacoustic imaging, the visualization of the complete processes of microwave excitation, energy deposition, thermoelastic expansion, and ultrasonic propagation is achieved, and photoacoustic signals are obtained; furthermore, for the photoacoustic image artifacts caused by the presence of bones, lungs, and other organs with highly mismatched sound speeds, the time of the photoacoustic signals of the organs with highly mismatched sound speeds detected by the ultrasonic transducer can be calculated through the simulation results, so as to weight the subsequent signal intensity to achieve the effect of removing artifacts.

[0068] The advantages and beneficial effects of the present invention are as follows:

[0069] A microwave photoacoustic imaging multi-physical field simulation detection method provided by the present invention, step S2 is based on the COMSOL coupling of multi-physical field interfaces of "Electromagnetic Waves, Frequency Domain", "Bioheat Transfer", "Solid Mechanics", and "Pressure Acoustics, Transient", and cooperates with step S501 to achieve the visualization of the complete physical process in microwave photoacoustic imaging, and detect and obtain photoacoustic signals for photoacoustic image reconstruction, which has important application value for both the in-depth study of subsequent photoacoustic basic theory and the guidance of detection schemes for biological tissue imaging experiments.

[0070] A microwave photoacoustic imaging multi-physical field simulation monitoring method provided by the present invention, step S502 realizes the visualization of the change results of the electromagnetic energy absorption density over time within any monitoring time by defining a monitoring time variable and combining a steady-state solver and parametric scanning, and subsequent photoacoustic signals generated at any monitoring moment can be calculated by adding multiple transient solvers; this method provides certain simulation guidance for photoacoustic imaging technology in monitoring research schemes for some dielectric and kinetic characteristics such as epilepsy and acupuncture.

[0071] The ingenuity of the present invention lies in: (1) Compared with traditional microwave thermoacoustic imaging simulations, the present invention couples the "electromagnetic waves, frequency domain", "bioheat transfer", "solid mechanics", and "pressure acoustics, transient" physical field interfaces in COMSOL, more completely and accurately describing the visualization of the microwave excitation, energy deposition, thermoelastic expansion, and ultrasonic propagation processes involved in microwave thermoacoustic imaging; (2) The present invention uses only one simulation software to implement the entire physical process of microwave thermoacoustic imaging, avoiding compatibility problems between different software, improving the matching between various steps in the simulation, reducing the simulation error, and providing more realistic simulation guidance for microwave thermoacoustic imaging experiments; (3) The present invention couples the "solid mechanics" physical field interface to calculate the mechanical properties of biological tissues, and combined with corresponding signal processing methods, can remove thermoacoustic image artifacts caused by bones, lungs, and other organs with highly mismatched sound speeds from the perspective of thermoacoustic signals; (4) The microwave thermoacoustic imaging simulation monitoring method proposed by the present invention can accurately capture the dynamic changes in the dielectric properties of biological tissues from the perspective of thermoacoustic signals, providing efficient and reliable simulation guidance for future research on dielectric property dynamics of this technology. BRIEF DESCRIPTION OF THE DRAWINGS

[0072] Figure 1 is a flowchart of a microwave thermoacoustic multi-physical field simulation detection and monitoring method provided by a preferred embodiment of the present invention;

[0073] Figure 2 is a schematic diagram of a simulation geometric model and physical field boundary conditions provided by an embodiment of the present invention;

[0074] Figure 3 is a flowchart of a method for suppressing bone artifacts based on the results of thermoacoustic multi-physical field simulation provided by another embodiment of the present invention;

[0075] Figure 4 is a model and thermoacoustic signal artifact suppression diagram provided by another embodiment of the present invention; (a) a diagram of a complex biological tissue model containing bones; (b) the original thermoacoustic signal detected by the 20th ultrasonic transducer; (c) the thermoacoustic signal after suppressing the trailing of the bone signal in (b); (d) the original thermoacoustic signal detected by the 105th ultrasonic transducer; (e) the thermoacoustic signal after suppressing the trailing of the bone signal in (d);

[0076] Figure 5 is a thermoacoustic imaging result diagram provided by another embodiment of the present invention; Figure 5 (a) The thermoacoustic image reconstructed from the original thermoacoustic signal, Figure 5 (b) The thermoacoustic image after suppressing bone artifacts. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0077] Next, the technical solutions in the embodiments of the present invention will be clearly and detailedly described in conjunction with the accompanying drawings in the embodiments of the present invention. The described embodiments are only a part of the embodiments of the present invention.

[0078] The technical solution of the present invention to solve the above technical problems is as follows:

[0079] Please refer to the attached Figure 1 , the present invention provides a microwave thermoacoustic imaging multi-physical field simulation detection and monitoring method. By sequentially coupling "electromagnetic wave, frequency domain", "bioheat transfer", "solid mechanics" and "pressure acoustics, transient" in COMSOL, the visualization of the complete physical process of microwave thermoacoustic imaging is realized, and the thermoacoustic signal of biological tissue obtained through multi-physical field simulation is acquired. Furthermore, the thermoacoustic image is reconstructed through an image algorithm; on this basis, by combining a steady-state solver and parametric scanning, the visualization of the change of electromagnetic energy absorption in biological tissue over time is realized, and multiple transient solvers are used to calculate the thermoacoustic signals at each monitoring time point. The specific steps are as follows:

[0080] S1. Construct a microwave thermoacoustic multi-physical simulation geometric model and set a microwave excitation function;

[0081] According to the experimental configuration of actual microwave thermoacoustic imaging, without affecting the results, in order to simplify the calculation amount and improve the calculation efficiency, a two-dimensional model is selected to be established based on the imaging layer for calculation. Please refer to the attached Figure 2 , the simulation geometric model consists of a rectangular waveguide, a coupling agent transformer oil, biological tissue, and an annular ultrasonic transducer array. Among them, the size of the rectangular waveguide is set to 100mm×72.14mm, and the size of the coupling agent transformer oil is set to 160mm×200mm; the ultrasonic transducers are evenly distributed on a circle with a radius of 40mm and the center of the circle located at the geometric center of the coupling agent transformer oil. According to the parameters of the thermoacoustic imaging system, the microwave excitation function an1(t) is as follows:

[0082] an1(t) = rect1(t - t 0 ) (1)

[0083] where rect1 represents a rectangular wave function with a pulse width of 550ns, and t 0 represents the rising edge time of the rectangular pulse;

[0084] S2. According to the thermoacoustic effect, set the physical field interfaces of "electromagnetic wave, frequency domain", "bioheat transfer", "solid mechanics" and "pressure acoustics, transient" and their corresponding boundary conditions.

[0085] The thermoacoustic effect involves the following four physical processes:

[0086] (1) Microwave irradiates biological tissue, and the biological tissue absorbs electromagnetic energy to generate a temperature rise;

[0087] (2) The biological tissue undergoes thermoelastic expansion due to the temperature rise;

[0088] (3) The biological tissues are squeezed against each other to generate ultrasonic waves;

[0089] (4) The ultrasonic waves propagate into the surrounding medium.

[0090] Therefore, the physical field interfaces selected in the present invention include the "electromagnetic wave, frequency domain", "bioheat transfer", "solid mechanics", and "pressure acoustics, transient" physical field interfaces. Please refer to the appendix Figure 2 , after the selection of the physical field interfaces is completed, corresponding boundary conditions need to be defined for each physical field interface according to the actual experimental situation.

[0091] Specifically as follows:

[0092] The "electromagnetic wave, frequency domain" physical field interface acts on the entire geometric model; among them, the right opening of the rectangular waveguide is used to radiate pulsed microwaves to the biological tissue, and the port boundary conditions need to be satisfied. The upper and lower boundaries of the waveguide are set as perfect electric conductors; considering that the heating process of the biological tissue by microwaves in microwave thermoacoustic imaging is different from the heating process of the biological tissue by an induction cooker, therefore, except for the boundary where the coupling agent transformer oil is in direct contact with the waveguide, the remaining boundaries are set as impedance boundaries; specifically as follows:

[0093] The input power of the port boundary condition is 1000 kW / m, the port type is set as rectangular, the port mode type is set as transverse electric, and the mode number is selected as 1;

[0094] The perfect electric conductor boundary condition:

[0095]

[0096] The impedance boundary condition:

[0097]

[0098] Among them, represents the unit normal vector, represents the electric field at the boundary of the coupling agent transformer oil, μ represents the magnetic permeability of the coupling agent transformer oil, ε represents the dielectric constant of the coupling agent transformer oil, σ represents the conductivity of the coupling agent transformer oil, ω represents the microwave angular frequency, represents the magnetic field at the boundary of the coupling agent transformer oil, represents the source electric field.

[0099] The interface of the physical fields of "bio-heat transfer" and "solid mechanics" acts on biological tissues. Considering that there is a certain heat exchange between biological tissues and the external coupling agent, transformer oil, and the pressure wave generated by the thermal expansion of biological tissues can smoothly propagate from the tissue boundary, the outermost boundary of biological tissues needs to satisfy both the "heat flux" boundary and the "low reflection" boundary simultaneously. Specifically as follows:

[0100] The heat flux boundary condition:

[0101]

[0102] The low reflection boundary condition:

[0103]

[0104] Among them, represents the conductive heat flux of biological tissues, h represents the heat transfer coefficient between biological tissues and the coupling agent, transformer oil, which is set to 5 W / (m 2 K), where W represents the power unit watt, m represents the distance unit meter, and K represents the temperature unit Kelvin; T ext represents the temperature of the coupling agent, transformer oil, which is set to 310.15 K, T represents the temperature of biological tissues, and σ s represents the Cauchy stress tensor of biological tissues, represents the mechanical impedance of biological tissues, represents the displacement occurred in biological tissues.

[0105] The interface of the physical field of "pressure acoustics, transient" acts on the coupling agent, transformer oil. Considering that biological tissues generate ultrasonic signals through vibration, therefore, the outermost boundary of biological tissues is set as the "acoustic-structure" boundary; in order to minimize the reflection of ultrasonic waves by the boundary as much as possible, the boundary of the coupling agent, transformer oil, is set as the "plane wave radiation" boundary. Specifically as follows:

[0106] The acoustic-structure boundary condition:

[0107]

[0108] The plane wave radiation boundary condition:

[0109]

[0110] Among them, p represents the sound pressure generated by the thermal expansion of biological tissues, ρ represents the density of the coupling agent, transformer oil, which is set to 860 kg / m 3 , v s represents the sound speed in the coupling agent, transformer oil, which is set to 1420 m / s.

[0111] S3. Set the model material parameters according to the selected physical field interface;

[0112] The "Electromagnetic wave, frequency domain" physical field interface needs to define the permittivity, permeability, and conductivity of the geometric model;

[0113] The "Bioheat transfer" physical field interface needs to define the constant pressure heat capacity, thermal conductivity, and density of biological tissues;

[0114] The "Solid mechanics" physical field interface needs to define the thermal expansion coefficient, bulk modulus, and shear modulus of biological tissues;

[0115] The "Pressure acoustics, transient" needs to define the sound speed and density of the coupling agent transformer oil.

[0116] S4. Reasonably divide the mesh according to the simulation requirements.

[0117] Considering that the microwave pulse width in the simulation is set to 550 ns and the high-frequency component of the generated thermoacoustic signal does not exceed 2 MHz, the maximum frequency f to be solved by the "Pressure acoustics, transient" physical field interface is set max = 2 MHz; in addition, to ensure the convergence of the calculation results, the maximum mesh size h max needs to meet the following conditions:

[0118]

[0119] S5. Add a steady-state solver and a transient solver for subsequent physical field result calculations; specifically, it includes the following steps:

[0120] S501. Multi-physics simulation detection solution settings;

[0121] S502. Multi-physics simulation monitoring solution settings;

[0122] Specifically, in step S501, the "Electromagnetic wave, frequency domain" physical field interface is calculated using a steady-state solver, and the center frequency is set to 3 GHz; three transient solvers are used to solve the "Bioheat transfer", "Solid mechanics", and "Pressure acoustics, transient" physical field interfaces; among them, for the "Bioheat transfer" physical field interface, the "Backward difference formula" method is selected in the solver configuration, the time step is set to 0.1 μs, and the solution time is set to 100 μs; for the "Solid mechanics" physical field interface, the "Generalized α" method is selected in the solver configuration, the time step is set to 0.1 μs, and the solution time is set to 100 μs; for the "Pressure acoustics, transient" physical field interface, the "Generalized α" method is selected in the solver configuration, the time step is set to The solution time is set to 100 μs;

[0123] Specifically, in step S502, the solver used for solving the physical field and the specific configuration of the solver are generally the same as those in step S501. The difference is that for the "electromagnetic wave, frequency domain" physical field interface, in addition to adding a steady-state solver, an additional variable needs to be defined as the monitoring time parameter, and parametric scanning needs to be added to calculate the electromagnetic energy absorption distribution corresponding to any time. For the "bioheat transfer", "solid mechanics", and "pressure acoustics, transient" physical field interfaces, every 3 transient solvers are used to calculate the thermoacoustic signals generated at a single monitoring moment. Then, the number of transient solvers required to monitor the thermoacoustic signals generated at multiple moments is:

[0124] n = 3n t (9)

[0125] where n represents the number of transient solvers added, and n t represents the number of monitoring moments;

[0126] More specifically, the "backward difference formula" method is a numerical calculation method used to approximate the derivative of a function, specifically expressed as:

[0127]

[0128] where ▽ represents the difference operator. For the solution of the "bioheat transfer" physical field interface, f(t i ) represents the temperature of the biological tissue at the t i th solution moment, m represents the maximum BDF order, which does not exceed 4 at most, is the combination number, representing the temperature coefficient of the biological tissue at the t i-j th solution moment;

[0129] More specifically, the "generalized α" method is a numerical method for solving second-order differential equations, specifically expressed as:

[0130]

[0131]

[0132] where α f , α m , γ, and β represent correction coefficients used to control the accuracy and stability of the equation solution, which are derived from the numerical dissipation coefficient ρ ∞ . They are set to 0.75 in the solutions of the "solid mechanics" and "pressure acoustics, transient" physical fields.

[0133] S6. Obtain the visualization results of microwave excitation, energy deposition, thermoelastic expansion, and ultrasonic propagation, calculate the thermoacoustic signals detected by each ultrasonic transducer, and complete the reconstruction of the thermoacoustic image; where

[0134] Obtain the visualization results of microwave excitation, energy deposition, thermoelastic expansion, and ultrasonic propagation, specifically as follows:

[0135] After simulation calculation, add a 2D plot group in the results and select the surface to analyze the specific absorption rate distribution, temperature distribution, displacement distribution, and acoustic field distribution of the simulation; add a 1D plot group and select the point result graph to analyze the changes in temperature, displacement, and pressure over time at a certain point within the biological tissue;

[0136] Calculate the photoacoustic signals detected by each ultrasonic transducer and complete photoacoustic image reconstruction, specifically as follows:

[0137] After simulation calculation, add a 2D cut point in the dataset, input the coordinate positions of each ultrasonic transducer, and change the dataset to "Pressure Acoustics, Transient"; on this basis, add a point calculation in the derived values, change the dataset to the 2D cut point, and change the expression to actd.p_t, then the photoacoustic signals detected by each ultrasonic transducer can be calculated; perform photoacoustic image reconstruction using the delay superposition algorithm based on the time-domain information and spatial position relationship of the photoacoustic signals detected by each ultrasonic transducer;

[0138] A microwave photoacoustic imaging multi-physics simulation detection and monitoring method proposed by the present invention includes coupling the "Solid Mechanics" physics field interface to calculate the mechanical properties of biological tissues and having the ability to analyze the characteristics of photoacoustic signals generated by hard tissues such as bones and stones; therefore, in another embodiment of the present invention, a method for suppressing bone artifacts in combination with the above simulation method is provided, and the specific steps are as follows:

[0139] Construct a complex biological tissue model containing bones and perform photoacoustic multi-physics field simulation;

[0140] Analyze the relative positions of each ultrasonic transducer set in the simulation and the characteristics of photoacoustic signals to determine the weight function w i (t) applied to each ultrasonic transducer;

[0141] Update the photoacoustic signals according to the weight function w i (t) and reconstruct the image with bone artifact suppression using the updated photoacoustic signals.

[0142] Specifically, please refer to Appendix Figure 4 , estimate the time t i when the ultrasonic transducer detects the bone signal according to the relative position between the ultrasonic transducer and the bone and in combination with the characteristics of the photoacoustic signal; specifically as follows:

[0143]

[0144] Among them, r i is the distance from the ultrasonic transducer to the bone, and v sThe sound velocity of the thermoacoustic signal in the coupling agent transformer oil is set to 1420 m / s.

[0145] Specifically, please refer to the appendix Figure 4 , the presence of bones causes a long tailing of the thermoacoustic signal, which appears as outwardly divergent circular artifacts in the thermoacoustic image. Therefore, the role of the weight function w i (t) is to suppress the signal tailing caused by bones from the perspective of the thermoacoustic signal, so as to achieve the purpose of eliminating circular artifacts in the image; specifically as follows:

[0146]

[0147] Among them, a is a constant, which affects the suppression effect of w i (t) and is affected by the relative position of the bone and the ultrasonic transducer; the closer the relative position of the ultrasonic transducer and the bone, the smaller the value of a, and the faster w i (t) suppresses the initial thermoacoustic signal. In this embodiment, the value range of a is (20, 40);

[0148] Specifically, use the updated thermoacoustic data to reconstruct the image with bone artifact suppression; specifically as follows:

[0149]

[0150] Among them, p' sum is the superimposed sound pressure obtained from the updated thermoacoustic data, and p i (t) is the intensity of the original thermoacoustic signal detected by the i-th ultrasonic transducer;

[0151] More specifically, please refer to Figure 5 , which is the reconstructed thermoacoustic image corresponding to the complex biological tissue model with bones provided in this embodiment. Among them, Figure 4 (a) is the thermoacoustic image reconstructed from the original thermoacoustic signal, Figure 5 (b) is the thermoacoustic image after suppressing bone artifacts; the results show that after suppressing the thermoacoustic signal tailing caused by bones, the artifacts brought by bones can be significantly reduced. Figure 5 (b) is the thermoacoustic image after suppressing bone artifacts; the results show that after suppressing the thermoacoustic signal tailing caused by bones, the artifacts brought by bones can be significantly reduced.

[0152] In the description of the present invention, it should be understood that the terms "including" and "comprising" indicate the presence of the described features, wholes, steps, operations, elements and / or groups, but do not exclude the presence or addition of one or more other features, wholes, steps, operations, elements and / or their collections.

[0153] The systems, devices, modules or units illustrated in the above embodiments can be specifically implemented by computer chips or entities, or by products with certain functions.

[0154] It should also be noted that the term "comprising", "including" or any other variant thereof is intended to cover non-exclusive inclusion, such that a process, method, commodity or device comprising a series of elements not only includes those elements but also includes other elements not expressly listed, or further includes elements inherent to such process, method, commodity or device. Without further limitation, an element defined by the phrase "comprising an..." does not exclude the presence of additional identical elements in the process, method, commodity or device comprising said element.

[0155] The above embodiments should be understood as being only for the purpose of illustrating the present invention and not for limiting the scope of protection of the present invention. After reading the content described in the present invention, those skilled in the art can make various changes or modifications to the present invention, and these equivalent changes and modifications also fall within the scope defined by the claims of the present invention.

Claims

1. A microwave thermoacoustic imaging multi-physics field simulation detection and monitoring method, characterized in that: The following steps are involved: S1. Construct a microwave thermoacoustic multi-physics simulation geometry model and set a microwave excitation function; S2. According to the thermoacoustic effect, set up the "Electromagnetic Waves, Frequency Domain", "Bioheat Transfer", "Solid Mechanics" and "Pressure Acoustics, Transient" physical field interfaces and their corresponding boundary conditions; S3. Set the material parameters of the model according to the selected physical field interface; S4. Divide the grid according to simulation requirements; S5. Add steady-state solver and transient solver for subsequent physical field result calculation; S6. Obtain visualization results of microwave excitation and energy deposition, thermoelastic expansion, and ultrasonic propagation, calculate the thermoacoustic signal detected by each ultrasonic transducer, and complete thermoacoustic image reconstruction.

2. A microwave thermoacoustic imaging multi-physics field simulation detection and monitoring method according to claim 1, characterized in that: The S1, constructing a microwave thermoacoustic multi-physics simulation geometric model and setting a microwave excitation function, specifically includes: The multi-physics simulation geometric model is set according to structural prior information provided by anatomy or imaging of the biological tissue being measured, including size and shape and relative positional relationships between tissues; The microwave excitation function an1(t) is as follows: an1(t)=rect1(t-t0) Among them, rect1 represents a rectangular wave function with a pulse width of 550ns, and t0 represents the rising edge time of the rectangular pulse.

3. A microwave thermoacoustic imaging multi-physics field simulation detection and monitoring method according to claim 1, characterized in that: In step S2, the "Electromagnetic Waves, Frequency Domain", "Bioheat Transfer", "Solid Mechanics" and "Pressure Acoustics, Transient" physical field interfaces are coupled in sequence to solve the microwave thermoacoustic multi-physics process; wherein, The "Electromagnetic Waves, Frequency Domain" and "Bioheat Transfer" physics field interfaces are used to solve the temperature rise caused by the absorption of electromagnetic energy in biological tissues; the "Bioheat Transfer" and "Solid Mechanics" physics field interfaces are used to solve the thermoelastic expansion of biological tissues caused by temperature rise; the "Solid Mechanics" and "Pressure Acoustics, Transient" physics fields are used to solve the thermoacoustic signals generated by the thermoelastic expansion of biological tissues.

4. According to claim 1, a microwave thermoacoustic imaging multi-physics field simulation detection and monitoring method is characterized in that: In step S2, the boundary conditions include electromagnetic boundary conditions, thermodynamic boundary conditions and acoustic boundary conditions; wherein the electromagnetic boundary conditions include port boundary conditions, ideal conductor boundary conditions and impedance boundary conditions; wherein, The input power of the port boundary condition is 1000kW / m, the port type is set to rectangular, the port mode type is set to transverse electric, and the modulus is selected as 1; The ideal conductor boundary condition is: The impedance boundary condition is: in, represents the unit normal vector, represents the electric field at the boundary of the coupling agent transformer oil, μ represents the magnetic permeability of the coupling agent transformer oil, ε represents the dielectric constant of the coupling agent transformer oil, σ represents the conductivity of the coupling agent transformer oil, ω represents the microwave angular frequency, represents the magnetic field at the boundary of the couplant transformer oil, represents the source electric field; The thermodynamic boundary conditions include: heat flux boundary conditions and low reflection boundary conditions; wherein, The heat flux boundary condition: The low-reflection boundary condition: in, represents the conductive heat flux of biological tissue, h represents the heat transfer coefficient between biological tissue and coupling agent transformer oil, which is set to 5 W / (m 2 K), where W represents the power unit watt, m represents the distance unit meter, and K represents the temperature unit Kelvin; T ext represents the temperature of the coupling agent transformer oil, which is set to 310.15K, T represents the temperature of the biological tissue, σ s represents the Cauchy stress tensor of biological tissue, represents the mechanical impedance of biological tissue, It indicates the displacement of biological tissue; The acoustic boundary conditions include: acoustic-structure boundary conditions and plane wave radiation boundary conditions; wherein, The acoustic-structure boundary condition: The plane wave radiation boundary condition: Where p represents the acoustic pressure generated by the thermal expansion of biological tissue, and ρ represents the density of the coupling agent transformer oil, which is set to 860 kg / m 3 , v s Represents the speed of sound in the coupling agent transformer oil and is set to 1420m / s.

5. The microwave thermoacoustic imaging multi-physics field simulation detection and monitoring method according to claim 1 is characterized in that: In the step S3, the material parameters of the model are set according to the selection of the physical field interface; specifically, the steps include: The Electromagnetic Waves, Frequency Domain physics interface requires the definition of the permittivity, magnetic permeability, and electrical conductivity of the geometry. The Bioheat Transfer physics interface requires the definition of the constant-pressure heat capacity, thermal conductivity, and density of biological tissues. The Solid Mechanics physics interface requires the definition of the thermal expansion coefficient, bulk modulus, and shear modulus of biological tissue; The "Pressure Acoustics, Transient" section requires the definition of the sound velocity and density of the couplant transformer oil.

6. The microwave thermoacoustic imaging multi-physics field simulation detection and monitoring method according to claim 1 is characterized in that: In step S4, the maximum grid size h max The following conditions need to be met: Among them, f max represents the maximum ultrasonic frequency that needs to be solved for the "Pressure Acoustics, Transient" physics field. Considering that the microwave pulse width in the simulation is set to 550 ns, the high-frequency component of the excited thermoacoustic signal does not exceed about 2 MHz. Therefore, f max is set to 2MHz.

7. The microwave thermoacoustic imaging multi-physics field simulation detection and monitoring method according to claim 1 is characterized in that: The step S5 specifically comprises the following steps: S501, multi-physics simulation solution setting; S502, multi-physics monitoring solution setting; In step S501, a steady-state solver is used to calculate the "Electromagnetic Waves, Frequency Domain" physical field interface, and the center frequency is set to 3 GHz; three transient solvers are used to solve the "Bioheat Transfer", "Solid Mechanics" and "Pressure Acoustics, Transient" physical field interfaces; for the "Bioheat Transfer" physical field interface, the "Backward Difference Formula" method is selected in the solver configuration, the time step is set to 0.1 μs, and the solution time is set to 100 μs; for the "Solid Mechanics" physical field interface, the "Generalized α" method is selected in the solver configuration, the time step is set to 0.1 μs, and the solution time is set to 100 μs; for the "Pressure Acoustics, Transient" physical field interface, the "Generalized α" method is selected in the solver configuration, and the time step is set to The solution time is set to 100 μs; In step S502, the solvers used for solving the corresponding physical fields and the specific configurations of the solvers are different in that, for the "Electromagnetic Waves, Frequency Domain" physical field interface, in addition to adding a steady-state solver, it is also necessary to define an additional variable as a monitoring time parameter, and add a parametric sweep to calculate the electromagnetic energy absorption distribution corresponding to any time; for the "Bioheat Transfer", "Solid Mechanics" and "Pressure Acoustics, Transient" physical field interfaces, every three transient solvers are used to calculate the thermoacoustic signal generated at a single monitoring moment, then the number of transient solvers required to monitor the thermoacoustic signals generated at multiple moments is: <h2 style=";text-align:left;direction:ltr">n=3n<h2 style=";text-align:left;direction:ltr"> t Where n is the number of transient solvers, t Indicates the number of monitoring moments.

8. A microwave thermoacoustic imaging multi-physics field simulation detection and monitoring method according to claim 7, characterized in that: The "backward difference formula" method is a numerical calculation method used to approximate the derivative of a function, which is specifically expressed as: in, represents the difference operator. For the solution of the Bioheat Transfer physics interface, f(t i ) represents the tth i The temperature of the biological tissue at the solution moment, m represents the maximum BDF order, which does not exceed 4. is the number of combinations, indicating the tth i-j The temperature coefficient of biological tissue at the solution time; The "generalized α" method is a numerical calculation method for solving second-order differential equations, which is specifically expressed as: Among them, α f , α m , γ and β are correction coefficients used to control the accuracy and stability of the solution of the equation, which is determined by the numerical dissipation coefficient ρ ∞ Exported; set to 0.75 in solvers for the Solid Mechanics and Pressure Acoustics, Transient physics.

9. The microwave thermoacoustic imaging multi-physics field simulation detection and monitoring method according to claim 1, characterized in that: In step S6, the visualization results of microwave excitation and energy deposition, thermoelastic expansion and ultrasonic propagation are obtained as follows: After the simulation calculation, add a 2D plot group in the results and select the surface to analyze the simulated specific absorption rate distribution, temperature distribution, displacement distribution, and acoustic field distribution; add a 1D plot group and select the point result graph to analyze the changes in temperature, displacement, and pressure over time at a certain point in the biological tissue; The calculation of the thermoacoustic signal detected by each ultrasonic transducer and the completion of thermoacoustic image reconstruction are specifically as follows: After the simulation calculation, add a two-dimensional cutoff point to the data set, enter the coordinate position of each ultrasonic transducer, and change the data set to "Pressure Acoustics, Transient"; on the current basis, add point calculations in the derived values, change the data set to a two-dimensional cutoff point, and change the expression to actd.p_t to calculate the thermoacoustic signal detected by each ultrasonic transducer; reconstruct the thermoacoustic image using the delayed superposition algorithm based on the time domain information and spatial position relationship of the thermoacoustic signal detected by each ultrasonic transducer.

10. A microwave thermoacoustic imaging multi-physics field simulation detection and monitoring method according to claim 9, characterized in that: The complete process of microwave excitation and energy deposition, thermal elastic expansion and ultrasonic propagation in the microwave thermoacoustic imaging is visualized, and thermoacoustic signals are obtained; furthermore, for thermoacoustic image artifacts caused by the presence of bones, lungs and other organs with highly mismatched sound velocities, the time when the ultrasonic transducer detects the thermoacoustic signal of the organ with highly mismatched sound velocities can be calculated through simulation results, thereby weighting the subsequent signal intensity to achieve the effect of removing artifacts.

Images