3D holographic high-performance simulation method and system based on multi-scale modeling theory

Abstract

The invention discloses a 3D holographic high-performance simulation method and system based on a multi-scale modeling theory, and the method comprises the steps: achieving the self-adaptive subdivision of a grid according to the change of wavelength through prior physical information, so as to capture the details of a microstructure, a complex organ and a plurality of systems in a human brain. Meanwhile, the wave velocity is incorporated into self-adaptive grid generation, and the limitation on watertight geometry in traditional CAD software is overcome. In order to avoid the overwhelming calculation cost in high-frequency wave simulation, an efficient mass lumped method is adopted under a high-performance calculation framework, and two-dimensional and three-dimensional numerical simulation of the human brain is realized in combination with triangular and tetrahedral grids. The numerical method is high in expandability and suitable for rapid and accurate simulation of whole-brain full-waveform. Through multi-scale whole-brain simulation, the method can be used for researching the influence of the skull on the transmission of the high-frequency ultrasound in the brain, and provides a new direction for the research and clinical application of the ultrasound in the field of neuromedicine.

Description

Technical Field

[0001] The present invention belongs to the field of biomedical simulation and relates to a 3D holographic high-performance simulation method and system based on multi-scale modeling theory, which can be used to study the influence of the skull on the propagation of high-frequency ultrasound in the brain, and provides a new direction for the research and clinical application of ultrasound in the field of neuromedicine. Background Art

[0002] High-frequency ultrasound imaging plays a vital role in brain medical research, and understanding the behavior of ultrasound in the brain at different scales, from systems to tissues, is essential to advancing high-resolution imaging. However, the complex internal structure of the brain poses a major challenge to accurate geometric modeling, and the numerical simulation of high-frequency waves incurs a lot of computational costs because of the need to consider numerous wavelengths in space. Therefore, it is necessary to develop more efficient and accurate computational models and methods. Prior art, the Chinese invention patent "A method for anatomical three-dimensional modeling of tomographic imaging of the striatum and hippocampal structure of the human brain" (publication number: CN109544682A) discloses a method for anatomical three-dimensional modeling of the striatum and hippocampal structure of the human brain. This method uses the 3D fast spoiled gradient echo sequence of the MRI system to obtain a horizontal high-resolution 3D structural image of the whole brain. The obtained image is then reconstructed in the coronal and sagittal planes on the workstation, adjusted to the optimal contrast, and observed simultaneously in the horizontal, coronal and sagittal positions. The relevant anatomical structures of the striatum and hippocampus are measured using the stereological method of quantitative analysis. After that, the MRI-measured striatum and hippocampal structure data of the human brain are analyzed using the SPSS13.0 statistical software, and the normal human brain striatum and hippocampal structure data are compared and analyzed one by one with the data in the disease state. Finally, the tomographic image is imported into the Mimics software to perform 3D modeling of the MRI data. However, this method relies on the limitations of commercial software and cannot achieve large-scale rapid generation of 3D models. At the same time, the generated model is limited by the quality of MRI data, and there are problems such as insufficient local accuracy and modeling deviation caused by artifacts. Moreover, the generated CAD model is not suitable for finite element simulation.

[0003] In pursuit of high-resolution ultrasound imaging and modeling of the human brain, on the one hand, it is necessary to generate high-quality mesh models, and on the other hand, it is necessary to combine advanced numerical calculation methods with high-performance computers to implement large-scale parallel computing. In mesh generation, it is simple and feasible to use unified differential mesh modeling (Guasch L, Calderón Agudo O, Tang MX, et al. Full-waveform inversion imaging of the human brain[J]. NPJ digitalmedicine, 2020, 3(1): 28.). However, due to the use of a unified mesh, this method faces accuracy limitations. In contrast, unstructured meshes provide flexibility for various material properties and facilitate simpler parameterization. For large-scale full-waveform simulations of high-frequency waves, finite element methods based on high-performance computing are effective. However, when executing an explicit time integration strategy, the classical finite element method needs to solve a linear system at each time step. The process of solving each time step greatly reduces the computational efficiency.

[0004] Based on this, the present invention proposes a 3D whole-brain holographic high-performance and precise simulation technology based on multi-scale modeling theory, and applies a grid mapping method based on prior physical information to realize the spatial discretization of the complex model inside the human brain. At the same time, the frequency, solution time step, grid size range and other information are integrated into the prior information to constrain the grid to generate a high-quality grid. In the finite element numerical calculation process, the mass lumping technology is implemented to diagonalize the mass matrix to reduce the computational cost. At the same time, a two-level parallel strategy based on MPI communication is used for large-scale parallel computing to achieve accurate modeling of the human brain model and full-wave numerical simulation of high-frequency ultrasound. Because this method is based on an advanced numerical algorithm solution under the high-performance computing framework, it can be used for full-wave modeling problems of any human brain model. Summary of the invention

[0005] In view of the defects in the prior art, the purpose of the present invention is to provide a 3D whole-brain holographic high-performance and precise simulation method and system based on multi-scale modeling theory.

[0006] The technical solution adopted by the present invention is as follows:

[0007] A 3D full-brain holographic high-performance and precise simulation method based on multi-scale modeling theory includes: based on the prior physical information of the human brain velocity field, adaptive meshing is achieved in combination with the minimum wavelength change adjustment; the generated mesh is used as a calculation model, and the high-frequency ultrasound is solved according to the residual form of the acoustic wave equation; the solution is achieved under a two-level integrated parallel high-performance computing framework, using a high-order finite element algorithm based on mass lumped units to achieve full-wave numerical simulation of high-frequency ultrasound.

[0008] In the above technical solution, further, the following steps are specifically included:

[0009] Step 1: According to the needs of simulation, select the prior physical information corresponding to the human brain velocity field, adjust the size of the constraint unit edge combined with the minimum wavelength change, and generate an adaptive grid as a calculation model;

[0010] Step 2: Add a perfectly matched layer to the model region to absorb the reflected waves from the boundary position;

[0011] Step 3: Implement mass lumping technology in the solver, add bubble functions to the mesh elements, and add high-order basis functions to the mesh elements to diagonalize the mass matrix to reduce the computational time of inverting the mass matrix in each iteration, and obtain a mass lumped high-order finite element solver.

[0012] Step 4: Place the generated grid computing model and the established mass lumped high-order finite element solver in a high-performance computing framework, and realize simulation by parallel computing the numerical solution of high-frequency ultrasonic waves at two levels of integration.

[0013] Furthermore, the a priori physical information includes: the center frequency of the excitation source, the number of grid points contained in each wavelength, the minimum time step, and the model velocity parameters.

[0014] Furthermore, the generating of the adaptive grid is specifically as follows:

[0015] The goal of optimizing a mesh is to maximize numerical accuracy while minimizing computational cost by adapting the mesh to local variations in the physical field properties;

[0016] The mesh size is constrained by the CFL condition, namely:

[0017]

[0018] Where h(x) is the diameter of the inscribed sphere or circle associated with the spatial position element, dim is the spatial dimension, which is 2 or 3, is the simulation time step, is the local acoustic longitudinal wave velocity, The Counrant number is a constant; from this, the minimum grid size can be determined;

[0019] At the same time, the length of the edges that make up the model grid unit Related to the center frequency of the source wavelet and the longitudinal velocity of the sound wave:

[0020]

[0021] Where G is the number of grids per wavelength, is the wavelength of the sound wave, For each unit degrees of freedom, is the center frequency; at the same time, the unit edge length Related to the parameter C per wavelength unit, , parameters C and G satisfy:

[0022]

[0023] in is a constant coefficient that is a function of the spatial polynomial k; this allows the grid size to be adjusted based on the minimum wavelength variation.

[0024] Furthermore, the thickness of the perfectly matched layer is set to be greater than a quarter of the minimum wavelength of the simulation.

[0025] Furthermore, the solver adopts an explicit time-marching method to solve the second-order acoustic wave equation, and uses mass lumping to diagonalize the mass matrix to reduce the amount of calculation for inverse solution of the mass matrix in each iteration. The diagonalization process includes: adding additional bubble functions in high-order grid cells, adding high-order basis functions on the grid cells, and using Lagrange basis functions and Newton-Cauchy integration rules to achieve diagonalization of the mass matrix.

[0026] Furthermore, the two-level integrated parallel computing is to realize two-level integrated parallel computing of the excitation source and the computing model in time and space, specifically:

[0027] In integrated parallelism, suppose n×n processes are called to calculate the model. Each process communicates with each other through MPI. The communicator is split into many sub-communicators. Each sub-communicator is referenced as an integrated member. All sub-communicators are divided into two groups, one for integrated sub-communicators and the other for source sub-communicators. For integrated sub-communicators, the computational grid is divided into multiple regions, each integrated member corresponds to one region, and finite element calculations with spatial parallelism are performed between integrated members. For source sub-communicators, each integrated member is equipped with an excitation source, and MPI communication is used between integrated members to ensure the synchronization of time evolution during the calculation execution.

[0028] A computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements a 3D whole-brain holographic high-performance and precise simulation method based on multi-scale modeling theory as described in any of the above items.

[0029] A 3D whole-brain holographic high-performance and accurate simulation system based on multi-scale modeling theory, including:

[0030] one or more processors;

[0031] A memory for storing one or more programs;

[0032] When the one or more programs are executed by the one or more processors, the one or more processors implement the 3D whole-brain holographic high-performance and precise simulation method based on multi-scale modeling theory as described in any of the above items.

[0033] Compared with the prior art, the present invention has the following significant advantages:

[0034] (1) This invention proposes an adaptive mesh method that goes beyond the challenges of traditional CAD: by adjusting the mesh size according to the wavelength change, a mesh generator is used to generate high-quality hierarchical meshes, where the mesh scale can be adjusted according to the changes in the physical field properties of the brain's internal tissues. High-fidelity, multi-scale human brain modeling is achieved. (2) This invention provides an advanced numerical simulation method based on high-performance computers: for large-scale simulation of high-frequency ultrasound, with high-order accuracy and low computational cost, provided by advanced mass lumped elements, ensuring efficient and accurate simulation of wave propagation. (3) High accuracy: The numerical stability of high-order elements is guaranteed by applying additional bubble functions to the mesh elements, and at the same time, high-order basis functions are used to improve the accuracy of the solution. (4) Scalability: The use of mass lumped elements with diagonal mass matrices eliminates the need to solve linear systems, thereby minimizing MPI communication and significantly improving parallel scalability. BRIEF DESCRIPTION OF THE DRAWINGS

[0035] Figure 1 It is a flow chart of the 3D whole-brain holographic high-performance and precise simulation method based on multi-scale modeling theory of the present invention.

[0036] Figure 2 It is a parallel architecture for high performance computing of the present invention: a structure of two-level integrated parallel processors of space and source.

[0037] Figure 3 This is the first application example of the present invention: a double-layer uniform sphere model.

[0038] Figure 4 These are the simulation results of the uniform sphere model: the high-order mass lumping (k=5) method and the pressure time series generated by COMSOL on the receiver.

[0039] Figure 5 is the signal difference between different solvers and the reference solution. Different degrees of high-order mass lumping methods (k=3,4,5) and the reference solution COMSOL L 2 The norm errors are 2.30%, 0.45%, and 0.15%, respectively. The L2 norm error between COMSOL k=5 and the reference solution is 0.18%.

[0040] Figure 6This is the second application example of the present invention: a human brain adaptive grid distribution map based on prior physical information mapping. The grid includes a 1-inch thick PML.

[0041] Figure 7 It is a two-dimensional human brain regional decomposition under a two-level integration framework.

[0042] Figure 8 These are snapshots of the wave field at different time steps with the skull.

[0043] Fig. 9 These are snapshots of the wave field at different time steps without the skull.

[0044] Fig.10 These are 50 time domain waveforms received with and without the skull (a. skull missing b. skull present).

[0045] Fig.11 This is the third application example of the present invention: three-dimensional human brain non-invasive ultrasound experiment setup.

[0046] Fig.12 It is an adaptive grid distribution map of the brain based on prior physical information in three-dimensional space that only displays local tissues or organs.

[0047] Fig.13 It is a three-dimensional human brain regional decomposition under a two-level integration framework.

[0048] Fig.14 These are the wave field snapshots (sagittal plane) at different time steps of a three-dimensional non-invasive ultrasound experiment on the human brain.

[0049] Fig.15 It is a snapshot of the wave field (coronal plane) at different time steps of a three-dimensional non-invasive ultrasound experiment on the human brain.

[0050] Fig.16 It is a snapshot (cross section) of the wave field at different time steps of a three-dimensional non-invasive ultrasound experiment on the human brain. DETAILED DESCRIPTION

[0051] The technical solution of the present invention is further described in detail below in conjunction with the accompanying drawings and specific examples. The following examples will help those skilled in the art to further understand the present invention, but do not limit the present invention in any form. It should be pointed out that, for those of ordinary skill in the art, several changes and improvements can be made without departing from the concept of the present invention. These all belong to the protection scope of the present invention.

[0052] The present invention proposes a 3D full-brain holographic high-performance and accurate simulation technology based on multi-scale modeling theory, which is used to construct a complex multi-scale structure model of the human brain and perform high-frequency ultrasonic full-wave numerical simulation at the same time. The brain is full of multi-scale structures, and a huge grid is required to capture tiny tissues and complex organs. Aiming at the problem of modeling the complex internal structure of the human brain, the present invention uses prior physical information such as frequency, wavelength, time step, speed, etc. to constrain the edges of the grid, and the generated grid units adapt to the changes in the shortest wavelength. This minimum wavelength constrained grid cell can easily capture the tiny structure of the human brain in space. Aiming at the large amount of computing resource consumption caused by the explicit time evolution strategy to invert the mass matrix at each time step, the present invention adopts the mass lumping method to diagonalize the mass matrix, reducing the amount of calculation for inverting the mass matrix in each iteration. Aiming at the challenge of large-scale problem calculation, the present invention proposes a combination of a high-order finite element method based on the mass lumping method and a high-performance computing framework based on MPI two-level integrated parallelism, which effectively utilizes computing resources and improves the computing efficiency of a wide range of problems.

[0053] According to a specific embodiment of the present invention, a 3D whole-brain holographic high-performance and precise simulation method based on multi-scale modeling theory includes: based on the prior physical information of the human brain velocity field, adaptive mesh division is achieved in combination with the minimum wavelength change adjustment, and the generated mesh is used as a calculation model to solve the high-frequency ultrasound according to the residual form of the acoustic wave equation; the solution is achieved in a two-level integrated parallel high-performance computing framework, using a high-order finite element algorithm based on mass lumped units to achieve full-wave numerical simulation of high-frequency ultrasound.

[0054] like Figure 1 , the method specifically comprises the following steps:

[0055] Step 1: According to the needs of the simulation, select the corresponding prior physical information, including the size of the edge of the constraint unit such as the center frequency of the excitation source, the number of grid points contained in each wavelength, the minimum time step, and the model velocity parameters, and adjust the grid size according to the change of the minimum wavelength to generate an adaptive grid; this method can effectively capture the tiny tissues and complex organs inside the human brain, build a multi-scale whole brain model, and reduce the degrees of freedom required for calculation; this step can adaptively obtain a two-dimensional triangular grid or a three-dimensional tetrahedral grid according to the numerical simulation requirements of the two-dimensional or three-dimensional human brain;

[0056] Step 2: Add a perfectly matched layer to the model region to absorb the reflected waves from the boundary position; the thickness of the perfectly matched layer can be selected according to the simulated wavelength and the geometric size of the model, and the thickness of the perfectly matched layer is preferably set to be greater than one quarter of the minimum simulated wavelength;

[0057] Step 3: Implement mass lumping technology in the solver to simplify the redundancy of linear systems caused by time and space discretization in wave problems, including adding additional bubble functions in the grid cells to ensure the stability of the numerical values ​​of high-order cells, and adding high-order basis functions on the grid cells to ensure the accuracy of the solution while reducing the computational time of inverting the mass matrix in each iteration; use Lagrange basis functions and Newton-Cauchy integration rules to achieve diagonalization of the mass matrix;

[0058] Step 4: Place the generated mesh model and the established mass lumped high-order finite element solver in a high-performance computing framework, and calculate the numerical solution of high-frequency ultrasound through two-level integrated parallel computing, that is: due to the need to achieve large-scale parallel computing, the excitation source and the model are integrated and parallelized in two levels in time and space, by splitting the communicator into many sub-communicators, divided into two groups, and each sub-communicator is regarded as an integrated member. Each integrated member in one group uses spatial parallelism in the sub-communicator in the usual way, and each integrated member in the other group of sub-communicators is set with an excitation source to allow communication between integrated members. A full-wave simulation of high-frequency ultrasound propagation inside the human brain is achieved.

[0059] The present invention can be used to reveal the influence of the skull on the propagation of high-frequency ultrasound in the brain, thereby promoting the research and clinical application of ultrasound in the field of neurology. The following is a detailed description of several key points in the scheme of the present invention:

[0060] 1. Acoustic wave equation applying PML strategy

[0061] Ultrasonic waves are high-frequency sound waves that propagate in space following the wave equation of sound waves. Numerical simulation requires simulating space as an infinite region and reducing the contamination of the received signal due to boundary reflection. In this invention, a fully matched layer is added as an absorbing boundary condition outside the physical calculation region. Considering the propagation of the second-order acoustic wave equation in a three-dimensional physical domain, the physical domain Contains the wave absorbing region surrounded by the outer boundary PML layer and the model solution domain , the acoustic wave equation in the free radiation region can be expressed in the rectangular coordinate system as follows:

[0062] (1)

[0063] In equation (1), v is a free parameter related to the wave velocity in space, called the longitudinal wave velocity. The acoustic wave equation is approximately solved using the continuous Galerkin method, and the residual form is obtained by introducing auxiliary variables and Fourier transform derivation.

[0064] (2)

[0065] In equation (2), and Respectively represent the auxiliary vector and auxiliary variable introduced. Its expression is

[0066] (3)

[0067] (4)

[0068] Assume that and Homogeneous Dirichlet boundary conditions are imposed on the ,right No boundary conditions are imposed. The first term of equation (2) introduces an external excitation source f. In the calculation area, for time , At this moment, we can use the residual form to get The coupling form of the modified second-order acoustic PML equations in the domain is:

[0069] (5)

[0070] (6)

[0071] (7)

[0072] (8)

[0073] (9)

[0074] (10)

[0075] In the formula, Representative at the border The external normal vector on is the source vector, and the excitation source is the time-varying Ricker wavelet. is the Lipschitz boundary, It is an unknown field. is the damping function The coefficient term of correlation.

[0076] 2. High-order mass lumping strategy

[0077] The explicit time marching method is used to solve the second-order acoustic wave equation. In order to avoid the large amount of computational resource consumption caused by the inverse of the mass matrix M at each time step in the explicit time marching strategy, mass lumping is used to diagonalize the mass matrix to reduce the computational complexity of the inverse of the mass matrix in each iteration.

[0078] Equation (3-8) in the three-dimensional open computational physics domain Medium damping function , the weak form of the acoustic wave equation under the action of a vector sound source can be expressed as:

[0079] (9)

[0080] set up Tetrahedral tessellation that satisfies the general properties of finite element meshes, h is The minimum inscribed circle diameter of each grid unit, the finite element space U h is defined as the space of continuous functions.

[0081] (10)

[0082] is a Hilbert space. Representation and The problem with the classical finite element method of fitting is to find satisfy:

[0083] (11)

[0084] use Respectively represent The associated degrees of freedom and the corresponding set of Lagrangian bases. For any function ,definition As the coefficient vector. For any unit, Equation (11) can be expressed as h The linear basis of Rewrite as

[0085] (12)

[0086] in, are the mass matrix and stiffness matrix respectively, expressed as

[0087] (13)

[0088] In the formula is the continuous inner product in L2.

[0089] Mass lumping is usually achieved using nodal basis functions and an inexact orthogonal rule for the mass matrix. When the integration points coincide with the nodes of the basis functions, a diagonal matrix is ​​obtained. During the diagonalization process, the mass matrix is ​​diagonalized using Lagrangian basis functions and Newton-Cauchy integration rules by adding additional bubble functions to the higher-order mesh elements and adding higher-order basis functions to the mesh elements.

[0090] Consider a polynomial finite element space of the form

[0091] (14)

[0092] New finite element subspace satisfy In the formula, , represents the space of high-order surfaces and internal bubble functions. When k = 2, the new finite element consists of 7 Lagrange interpolation points. There are three equivalent class points, including three vertices {S 1 ,S 2 ,S 3}, the midpoints of the three sides {M 1 ,M 2 ,M 3} and center of gravity G.

[0093] Constructing a new finite element subspace Is to use a node base , express The collection of all nodes in . Unit node From the reference node Therefore, by mapping the node basis to the physical unit, the node basis function of the element can be obtained. The integral of the mass matrix at the unit cell is approximated by the orthogonal rule as follows:

[0094] (15)

[0095] In the formula For unit The weight of the whole domain. It can be expressed as:

[0096] (16)

[0097] According to the properties of the node basis function, , represents the Kronecker function. Therefore, the corresponding mass matrix M ij It can be expressed as

[0098] (17)

[0099] Because the unit There are three equivalent classes of points in , so the same weights are assigned to the points of the same class. The integral formula can be written as follows

[0100] (18)

[0101] For higher-order finite elements, the quadrature formula can be derived in the same way. Compared with the classical finite element method numerical method, the finite element method using mass lumping shows superior stability, convergence and efficiency.

[0102] 3. Adaptive grid technology based on prior physical information mapping

[0103] Next, we introduce adaptive meshing techniques to construct the finite element mesh required in the numerical calculation process. In this work, the goal of optimizing the mesh is to maximize the numerical accuracy while minimizing the computational cost by adapting the mesh to local changes in the physical field properties. In order to ensure numerical stability, the mesh size must be constrained by the CFL (Courant-Friedrichs-Lewey) condition, because the time step is affected by the minimum unit. Therefore, the mesh generation program needs to ensure that the minimum unit is as large as possible to avoid the simulation time step being too small. For acoustic waves, the CFL condition that limits the mesh size can be described as

[0104] (19)

[0105] where h(x) is the diameter of the inscribed sphere or circle associated with the spatial position element, dim is the spatial dimension of the problem (2 or 3), is the simulation time step, is the local acoustic longitudinal wave velocity. When the Counrant number is defined as a constant (such as ), for a given and , rearrange equation (19) to find the possible minimum grid size.

[0106] Minimizing the computational cost requires the mesh to adapt to local wavelength changes and adjust the mesh size according to the minimum wavelength change; it is assumed that all triangular elements that make up the model are equilateral, which is necessary for accurate simulation using the finite element method. In this case, the edge length of the triangular element is It can be related to the center frequency of the source wavelet and the longitudinal wave velocity of the sound wave.

[0107] (20)

[0108] Where G is the number of grids per wavelength, is the wavelength of the sound wave, For each unit degrees of freedom, is the center frequency. At the same time, the unit edge length It can be related to the per-wavelength unit parameter c. , parameters C and G are related to each other in the following way:

[0109] (twenty one)

[0110] in is a constant coefficient that is a function of the spatial polynomial k. Mass lumped elements have a higher number of nodes per element, so high-order mass lumped elements have higher values ​​for each polynomial degree than classical finite elements.

[0111] 4. Two-level integrated parallel strategy based on high performance computing framework

[0112] In integrated parallelism, assume that n×n processes are called to compute the model. Each process communicates with each other through MPI, such as Figure 2 As shown. The communicator is split into many sub-communicators, each of which is referenced as an integrated member. All sub-communicators are divided into two groups, one for integrated sub-communicators and the other for source sub-communicators. For integrated sub-communicators, in each integrated member, it is allowed to specify existing functions for solving finite element problems that use spatial parallelism across sub-communicators in the usual way. Communication between integrated members is allowed in the other group of source sub-communicators. By dividing the size of the original communicator by the number of processes in each integrated member, each integrated member will have the same spatial parallelism as the number of integrated members. Therefore, the total number of processes initiated by mpiexec must be equal to the product of the number of ensemble members and the number of processes used by each ensemble member. In the actual calculation process, each integrated member in the source sub-communicator is equipped with an excitation source, and MPI communication is used between the integrated members to ensure the synchronization of time evolution during the calculation execution. The integrated members in the integrated sub-communicator divide the computational grid into multiple regions using sub-MPI communication.

[0113] Three numerical examples are given below:

[0114] (1) This case is to verify the accuracy and computational efficiency of the high-order finite element method using mass lumped elements. First, a two-layer uniform sphere model is set up, such as Figure 3 As shown. The wave velocity of the upper medium of the model is 1500m / s, the wave velocity of the lower medium is 3200m / s, and the wave velocity of the ball is 2450m / s. The model area size is 12-inch×12-inch, the excitation source coordinates are (6, 1.8), and the receiver coordinates are (6, 2.5). In order to verify the numerical accuracy, the model calculation uses high-order mass lumped finite elements of different degrees. The grid is encrypted 10 times, and the solution calculated by the commercial simulation software COMSOL based on the discontinuous Galerkin finite element numerical method is used as the reference solution. Compare the time domain waveforms obtained at the receiver between different methods, as shown Figure 4 shown.

[0115] The relationship between the degree of freedom (DoF) and resource consumption when running the program using different finite element solvers is shown in Table 1. When using high-order mass lumped ML, the numerical solution is close to the reference solution. However, the computational time scale is reduced by about 3 times.

[0116] Solver Degrees of Freedom (DoF) Time Consumption Memory (MB) ML k=3 169848 9266 398 ML k=4 287219 14789 443 ML k=5 560818 10425 515 COMSOL k=5 546798 31065 498 Reference 5172360 66324 4852 Table 1

[0117] Figure 5 The signal differences and L between the high-order mass lumping method and the reference solution are given. 2 Norm error. L 2 The error is calculated using the following formula:

[0118] (twenty two)

[0119] In the formula, and are the acoustic pressures of the numerical solution and the reference solution, respectively. The results show that the relative numerical error between the high-order mass lumping method and the reference solution decreases with the increase of the order. 2 The norm error does not exceed 0.25%, indicating that the high-order mass lumping method has good numerical stability when using high-order units. The numerical method of the present invention has a high-order accuracy of 5th order and has strong scalability, which is suitable for fast and accurate simulation of the whole brain and full waveform.

[0120] (2) This case considers the simulation of a human brain model placed in water, and the speed of sound waves in water is set to a constant 1500m / s. The model area size is 12-inch×12-inch, and the sound wave speed in various tissues inside the human brain varies between 1500m / s-3000m / s. A time-varying Ricker wavelet with a center frequency of 5.9MHz is used as the excitation source, and a 1-inch thick unit grid is wrapped outside to act as a PML layer to absorb the reflected waves from the truncation boundary. The mesh unit size is a function of the P wave velocity of the sound wave as a variable. The constant G is set to 5, and a mass lumped unit with k=3 and a time step parameter of 0.0001 s are used to constrain the unit edge length. The computational domain is approximately 3600 wavelengths along one direction relative to the maximum truncation frequency.

[0121] Figure 6The adaptive mesh of the two-dimensional human brain generated by mapping the prior physical information is presented. This minimum wavelength constrained mesh unit can easily capture the tiny structure of the human brain in space, but it has a large cross-scale property. For example, after a local magnification of the cerebellum region by 256 times, the wave velocity changes between the cerebellar cortex, medulla, and brainstem lead to scale changes within the mapping mesh space. In addition, the adaptive mesh conforms to the complex geometric structure of a variety of brain tissues, reducing the degrees of freedom of the computational space compared to the uniform structured mesh. This method meets the accuracy requirements of specific solutions and adapts to local changes in physical properties.

[0122] Figure 7 The regional decomposition diagram of the two-dimensional human brain is given under the two-level integrated parallel high-performance computing framework. The model uses 80 cores for parallel computing, including 10 excitation sources, and each excitation source divides the space into 8 regions for distributed memory parallel computing during the time evolution process. This method greatly improves the computing efficiency.

[0123] Figure 8 , Fig. 9 The ultrasound propagation process between tissues in the human brain at different time steps calculated using a high-order mass lumped element combined with a two-level integrated parallel framework is shown. Figure 8 ) and without skull ( Fig. 9 ) in the case of transcranial ultrasound. Fig.10 The time domain waveforms received by the receiving array with and without the skull are given.

[0124] (3) This case further considers a non-invasive ultrasound experiment of the human brain in three-dimensional space. The simulation places a multimodal human brain model generated based on prior physical information mapping in a water tank. The model area size is 12-inch×12-inch×8.75-inch, the excitation source is placed on the right side of the water tank, and a 400×400 hydroacoustic array for receiving ultrasound signals is placed on the left side, as shown in Fig.11 The model area is discretized using mass lumped elements with a degree of 3 (k=3), and the mesh model generated using adaptive mesh technology based on prior physical information is shown in Fig.12 As shown, the number of mesh elements in the discrete region is 16,450,340, resulting in 114,731,163 degrees of freedom.

[0125] The software runs on a server cluster consisting of six computing nodes under a two-level integrated parallel framework. The server cluster includes four computing nodes with 56-core Intel Gold 6348 CPUs, one computing node with 112 cores, and one management node with 24 cores. The actual execution process uses 300 cores and performs parallel calculations across multiple nodes. The regional decomposition of the model during operation is as follows Fig.13 shown.

[0126] Fig.14 , Fig.15 , Fig.16 The wavefield snapshots of high-frequency ultrasound in the human brain at different time steps are shown. By providing wavefield snapshots in the sagittal, coronal and transverse planes, the propagation of ultrasound in three-dimensional human brain tissue at different time points is shown. These wavefield snapshots provide valuable insights into the propagation laws and behavioral characteristics of ultrasound in the complex brain environment.

[0127] It will be appreciated by those skilled in the art that embodiments of the present invention may be provided as methods, systems, or computer program products. Therefore, the present invention may take the form of a complete hardware embodiment, a complete software embodiment, or an embodiment combining software and hardware. Furthermore, the present invention may take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0128] The present invention is described with reference to flowcharts and / or block diagrams of methods, devices (systems), and computer program products according to embodiments of the present invention. It should be understood that each process and / or block in the flowchart and / or block diagram, as well as the combination of processes and / or blocks in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing device to produce a machine, so that the instructions executed by the processor of the computer or other programmable data processing device generate instructions for implementing the processes in the flowchart and / or block diagram. Figure 1 A process or multiple processes and / or boxes Figure 1 A device that provides the functions specified in a block or multiple blocks.

[0129] These computer program instructions may also be stored in a computer-readable memory capable of directing a computer or other programmable data processing device to operate in a specific manner, so that the instructions stored in the computer-readable memory produce an article of manufacture comprising an instruction device, which implements the process Figure 1 A process or multiple processes and / or boxes Figure 1 A function specified in one or more boxes.

[0130] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operating steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing instructions for implementing the process in the computer or other programmable device. Figure 1 A process or multiple processes and / or boxes Figure 1 A step that specifies a function in one or more boxes.

[0131] The above describes the specific embodiments of the present invention. It should be understood that the present invention is not limited to the above specific embodiments, and those skilled in the art may make various changes or modifications within the scope of the claims, which will not affect the essence of the present invention.

Claims

1. A 3D holographic high-performance simulation method based on multi-scale modeling theory, characterized in that: include: Based on the prior physical information of the human brain velocity field, adaptive mesh generation is achieved in combination with the minimum wavelength change adjustment. The generated mesh is used as a calculation model, and the high-frequency ultrasonic wave is solved according to the residual form of the acoustic wave equation. The solution is achieved in a two-level integrated parallel high-performance computing framework using a high-order finite element algorithm based on mass lumped units to achieve full-wave numerical simulation of high-frequency ultrasonic waves.

2. The 3D holographic high-performance simulation method based on multi-scale modeling theory according to claim 1 is characterized in that: The specific steps include: Step 1: According to the needs of simulation, select the prior physical information corresponding to the human brain velocity field, adjust the size of the constraint unit edge combined with the minimum wavelength change, and generate an adaptive grid as a calculation model; Step 2: Add a perfectly matched layer to the model region to absorb the reflected waves from the boundary position; Step 3: Implement mass lumping technology in the solver, add bubble functions to the mesh elements, and add high-order basis functions to the mesh elements to diagonalize the mass matrix to reduce the computational time of inverting the mass matrix in each iteration, and obtain a mass lumped high-order finite element solver. Step 4: Place the generated grid computing model and the established mass lumped high-order finite element solver in a high-performance computing framework, and realize simulation by parallel computing the numerical solution of high-frequency ultrasonic waves at two levels of integration.

3. The 3D holographic high-performance simulation method based on multi-scale modeling theory according to claim 2 is characterized in that: The prior physical information includes: the central frequency of the excitation source, the number of grid points contained in each wavelength, the minimum time step, and the model velocity parameters.

4. The 3D holographic high-performance simulation method based on multi-scale modeling theory according to claim 2 is characterized in that: The generating of the adaptive grid is specifically as follows: The goal of optimizing a mesh is to maximize numerical accuracy while minimizing computational cost by adapting the mesh to local variations in the physical properties; The mesh size is constrained by the CFL condition, namely: ， Where h(x) is the diameter of the inscribed sphere or circle associated with the spatial position element, dim is the spatial dimension, which is 2 or 3, is the simulation time step, is the local acoustic longitudinal wave velocity, The Counrant number is a constant; from this, the minimum grid size can be determined; At the same time, the length of the edges that make up the model grid unit Related to the center frequency of the source wavelet and the longitudinal velocity of the sound wave: ， Where G is the number of grids per wavelength, is the wavelength of the sound wave, For each unit degrees of freedom, is the center frequency; at the same time, the unit edge length Related to the parameter C per wavelength unit, , parameters C and G satisfy: ， in is a constant coefficient that is a function of the spatial polynomial k; this allows the grid size to be adjusted based on the minimum wavelength variation.

5. The 3D holographic high-performance simulation method based on multi-scale modeling theory according to claim 2 is characterized in that: The thickness of the perfectly matched layer is set to be greater than a quarter of the minimum wavelength of the simulation.

6. The 3D holographic high-performance simulation method based on multi-scale modeling theory according to claim 2 is characterized in that: The solver adopts the explicit time marching method to solve the second-order acoustic wave equation, and adopts mass lumping to diagonalize the mass matrix to reduce the amount of calculation for inverse solution of the mass matrix in each iteration. The diagonalization process includes: adding additional bubble functions in high-order grid cells, adding high-order basis functions on the grid cells, and using Lagrange basis functions and Newton-Cauchy integration rules to realize the diagonalization of the mass matrix.

7. The 3D holographic high-performance simulation method based on multi-scale modeling theory according to claim 2 is characterized in that: The two-level integrated parallel computing is to realize two-level integrated parallel computing of the excitation source and the computing model in time and space, specifically: In integrated parallelism, suppose n×n processes are called to calculate the model. Each process communicates with each other through MPI. The communicator is split into many sub-communicators. Each sub-communicator is referenced as an integrated member. All sub-communicators are divided into two groups, one for integrated sub-communicators and the other for source sub-communicators. For integrated sub-communicators, the computational grid is divided into multiple regions, each integrated member corresponds to one region, and finite element calculations with spatial parallelism are performed between integrated members. For source sub-communicators, each integrated member is equipped with an excitation source, and MPI communication is used between integrated members to ensure the synchronization of time evolution during the calculation execution.

8. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the program is executed by a processor, it implements the 3D holographic high-performance simulation method based on multi-scale modeling theory as described in any one of claims 1-7.

9. A 3D holographic high-performance simulation system based on multi-scale modeling theory, characterized in that: include: one or more processors; A memory for storing one or more programs; When the one or more programs are executed by the one or more processors, the one or more processors implement the 3D holographic high-performance simulation method based on multi-scale modeling theory as described in any one of claims 1-7.

Images