A real-time sound field simulation method and system based on physical hierarchical decoupling
By employing a physically layered decoupled acoustic field simulation method, combined with Rayleigh integral and full-wave numerical simulation, the problem of computational time consumption or insufficient accuracy in traditional methods is solved, achieving real-time and efficient acoustic field simulation, which is suitable for low-intensity transcranial focused ultrasound neuromodulation.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- SICHUAN UNIV
- Filing Date
- 2026-05-26
- Publication Date
- 2026-07-31
AI Technical Summary
Traditional sound field simulation methods suffer from excessive computation time or insufficient accuracy in low-intensity transcranial focused ultrasound neuromodulation, making it difficult to meet the needs of real-time treatment navigation and closed-loop control.
By employing a physical layering and decoupling method, the sound field is rapidly calculated by constructing the transfer matrix offline and combining it with Rayleigh integrals. Layered processing is applied to homogeneous media and complex heterogeneous skull regions, achieving fast and high-precision sound field simulation.
While maintaining high precision, it achieves a speedup of hundreds of times in sound field simulation calculation, with a single prediction time of less than 0.06 seconds, meeting the navigation and closed-loop control requirements of real-time transcranial ultrasound neuromodulation.
Smart Images
Figure CN122490838A_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of medical ultrasound imaging technology, and in particular to a real-time sound field simulation method and system based on physical layer decoupling. Background Technology
[0002] Low-intensity transcranial focused ultrasound (LCI) neuromodulation technology offers advantages such as non-invasiveness and high spatial resolution, but its clinical application is limited by phase distortion caused by the skull. Traditional sound field simulation methods (such as k-Wave full-wave simulation) are highly accurate but computationally expensive, making it difficult to meet the needs of real-time treatment navigation and closed-loop control. While existing deep learning methods have accelerated simulation, they suffer from low accuracy and sensitivity to noise.
[0003] Traditional transcranial ultrasound simulation schemes face a dilemma of balancing speed and accuracy: if numerical simulation (such as k-Wave) is used throughout, although it can accurately simulate the complex effects of the skull, the calculation time is too long and cannot meet the real-time requirements; if analytical approximation (such as simplifying the skull into a homogeneous layer) is used throughout, although the calculation is fast, it cannot accurately characterize the scattering, refraction and attenuation caused by the heterogeneity of the skull (such as the spatial variation of sound velocity, density and porosity), resulting in serious insufficient simulation accuracy. Summary of the Invention
[0004] The purpose of this invention is to provide a real-time sound field simulation method and system based on physical layer decoupling. By constructing the transfer matrix from the virtual source layer to the target region offline and combining Rayleigh integral to quickly calculate the propagation from the transducer to the virtual source layer, the computational burden of online sound field prediction is significantly reduced.
[0005] This invention is achieved using the following technical solution: a real-time sound field simulation method based on physical layer decoupling, comprising the following steps: S1: Establishing an acoustic parameter mesh based on CT images, including homogeneous media regions and complex heterogeneous skull regions: S2: Define a virtual source layer in the uniform medium region between the transducer and the skull as a layered decoupling interface between analytical calculation and numerical calculation. S3: Using the Rayleigh-Sommerfeld integral analytical formula, the first transfer matrix from the actual transducer array elements to each virtual source point on the virtual source layer is calculated, thereby quickly obtaining the complex sound pressure distribution of the virtual source layer; S4: For complex and heterogeneous skull regions, a full-wave numerical simulation method is used to pre-calculate the frequency response coefficients from each virtual source point to each observation point in the target region, and to construct and store the second transfer matrix; S5: The excitation vector of the actual transducer is multiplied sequentially with the first transfer matrix and the second transfer matrix to synthesize the complete sound field distribution in the target area in real time.
[0006] Furthermore, step S1 includes the following sub-steps: S11: Acquire skull CT images and calculate bone porosity distribution based on Hounsfield units of the CT images; S12: Based on the bone porosity distribution, the sound velocity, density, and attenuation coefficient distribution of each voxel point in the skull are mapped to form an acoustic parameter grid that includes both homogeneous medium regions and complex heterogeneous skull regions.
[0007] Furthermore, the density and velocity of sound of the voxel point are calculated using the following formulas: density =phi *ρw+ (1 - phi) *ρb; speed =phi * sw + (1 - phi) *sb; Where density is the density of voxel points; phi is the bone porosity; ρw is the water density; ρb is the bone density; sw is the water sound velocity; and sb is the bone sound velocity.
[0008] Furthermore, step S2 includes the following sub-steps: S21: Based on the spatial coordinates of the CT image, obtain the coordinates of the highest point on the outer surface of the skull on one side of the transducer. S22: Based on the coordinates of the highest point, offset a preset distance towards the transducer to form the plane where the virtual source layer is located; S23: Based on the scaling ratio of the acoustic parameter grid in step S1, convert the physical coordinates of each virtual source point on the virtual source layer into a grid index and record the number.
[0009] Furthermore, the first transfer matrix is calculated as follows: ; in, is the element of the first transfer matrix; k is the wave number; ρ0 and c are the density and speed of sound of the homogeneous medium, respectively; ΔS is the area of the transducer element; Rmn is the distance from the nth transducer element to the mth virtual source point; and i is the imaginary unit.
[0010] Furthermore, the full-wave numerical simulation method includes the following sub-steps: The acoustic parameter mesh obtained in step S1 is used as the simulation medium input; Each virtual source point on the virtual source layer is used as an excitation source in turn, and a reference excitation signal is set to perform steady-state sound field simulation. Extract the time-domain sound pressure signal at each observation point within the target area, and obtain the complex amplitude at the target frequency through Fourier transform, which is the frequency response coefficient; The frequency response coefficients between each virtual source point and each observation point are organized into a second transfer matrix H, with the dimension of the second transfer matrix being Nroi × Npres, where Nroi is the number of observation points in the target region and Npres is the number of virtual source points on the virtual source layer.
[0011] Furthermore, the expression for the complete sound field distribution is: ; Among them, E transducer H is the complex excitation vector of the N×1 dimensional transducer array elements; H is the pre-computed transfer matrix of P×M dimensionality; P target Let be the P×1 dimensional complex sound pressure vector of the target region; T is the M×N dimensional Rayleigh integral transfer matrix.
[0012] Furthermore, the homogeneous medium region is an aqueous medium or a homogeneous coupling agent layer between the transducer and the skull, and the propagation of sound waves in this region follows the homogeneous medium wave equation; the complex heterogeneous skull region contains spatial variations in sound velocity, density, and porosity, resulting in sound wave scattering, refraction, and attenuation effects.
[0013] Furthermore, it also includes an adaptive adjustment step for the transducer excitation signal during the online phase: Receive transducer operating parameters or desired treatment focus location in real time; Based on the acoustic field distribution of the target region synthesized in step S5, the complex excitation vector E of each element of the transducer is optimized in reverse. transducer ; Repeat steps S3 to S5 to quickly update the sound field simulation results until the preset focusing accuracy or sound field distribution requirements are met.
[0014] A real-time sound field simulation system based on physical layer decoupling, used to implement the aforementioned real-time sound field simulation method, includes: Personalized modeling module: used to create acoustic parameter meshes containing homogeneous media regions and complex heterogeneous skull regions based on patient CT images; Virtual source layer definition module: used to set up a virtual source layer near the interface between the homogeneous medium and the skull, and to calculate its mesh index; Online analytical calculation module: used to calculate the first transfer matrix from the transducer array elements to the virtual source layer online using the Rayleigh-Sommerfeld integral analytical formula; Offline pre-calculation module: used to pre-calculate and store the second transfer matrix from the virtual source layer to the target region offline using the full-wave numerical simulation method; Real-time sound field synthesis module: used to multiply the transducer excitation vector with the first transfer matrix and the second transfer matrix in sequence, and output the sound field distribution of the target area in real time.
[0015] The beneficial effects of this invention are as follows: This invention employs a segmented processing and collaborative acceleration design: For the homogeneous medium segment, rapid analytical calculation is used: in the water medium (or homogeneous coupling agent) layer between the transducer and the skull, sound wave propagation follows a simple homogeneous medium wave equation, which can be calculated extremely quickly and accurately using the Rayleigh-Sommerfeld integral analytical formula. For the complex heterogeneous structure segment, pre-calculated numerical simulation is used: for the target region containing a complex heterogeneous skull, high-precision numerical simulation (such as k-Wave) is used for a one-time offline calculation, and the results are compressed and stored as an acoustic transfer matrix. Efficient synthesis is performed in the online stage: during online real-time simulation, only the excitation vector representing the transducer state needs to be mapped to the virtual source layer through a rapid analytical step (Rayleigh integral), and then multiplied once with the pre-calculated transfer matrix to instantly synthesize the complete sound field of the target region. This design places all time-consuming numerical calculations in the offline stage, while the online stage only performs analytical calculations and matrix multiplications, thus achieving a speedup of hundreds of times while maintaining the accuracy of numerical simulation.
[0016] Experimental results show that the present invention achieves more than 250 times the speedup while maintaining an accuracy almost identical to k-Wave (NRMSE < 0.6%), with a single prediction time of less than 0.06 seconds, making it suitable for navigation and closed-loop control of real-time transcranial ultrasound neuromodulation. Attached Figure Description
[0017] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on the structures shown in these drawings without creative effort.
[0018] Figure 1 This is a flowchart of the present invention. Detailed Implementation
[0019] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. The components of the embodiments of the present invention described and shown in the accompanying drawings can generally be arranged and designed in various different configurations.
[0020] It should be noted that similar labels and letters in the following figures indicate similar items. Therefore, once an item is defined in one figure, it does not need to be further defined and explained in subsequent figures.
[0021] The following detailed description of some embodiments of the present invention is provided in conjunction with the accompanying drawings. Unless otherwise specified, the following embodiments and features can be combined with each other.
[0022] See Figure 1 A real-time sound field simulation method based on physical layer decoupling includes the following steps: S1: Establishing an acoustic parameter mesh based on CT images, including homogeneous media regions and complex heterogeneous skull regions: S2: Define a virtual source layer in the uniform medium region between the transducer and the skull as a layered decoupling interface between analytical calculation and numerical calculation. S3: Using the Rayleigh-Sommerfeld integral analytical formula, the first transfer matrix from the actual transducer array elements to each virtual source point on the virtual source layer is calculated, thereby quickly obtaining the complex sound pressure distribution of the virtual source layer; S4: For complex and heterogeneous skull regions, a full-wave numerical simulation method is used to pre-calculate the frequency response coefficients from each virtual source point to each observation point in the target region, and to construct and store the second transfer matrix; S5: The excitation vector of the actual transducer is multiplied sequentially with the first transfer matrix and the second transfer matrix to synthesize the complete sound field distribution in the target area in real time.
[0023] In this embodiment, step S1 specifically includes: Based on CT images, Hounsfield units of the skull are obtained to calculate bone porosity and map the distribution of sound velocity, density, and attenuation coefficient. First, the CT image is read, and the specific Hounsfield values are obtained. These values are then cropped and mapped to a grid of a certain size, for example, a 256^3 voxel size scaled to a 110^3 grid. This yields a scaling factor, allowing us to obtain the specific formula for the phi value (bone porosity). It is assumed that HU=0 corresponds to pure water (phi=1), and HU=1500 corresponds to pure bone (phi=0).
[0024] Hu = double(ct resized ); phi = 1.0 - (Hu / 1500.0); Define water density w-ρ (water-ρ), bone density b-ρ (bone-ρ), water sound velocity wss (water-sound-speed), and bone sound velocity bss (bone-sound-speed). Calculate the acoustic parameters using the formulas. `medium` represents the medium mesh, and `Voxel` is the mesh index, used for subsequent input into the kwave calculation library. A fixed attenuation range of 0.13~5.33 dB / (MHz·cm) and a power-law exponent alpha are used.power =1.5.
[0025] medium.density = phi * w-ρ + (1 - phi) * b-ρ; medium.sound speed = phi * wss + (1 - phi) * bss; alpha base = 0.13 + (5.33 - 0.13) * sqrt(phi); medium.alpha coeff = alpha base ; medium.alpha power = 1.5.
[0026] In this embodiment, step S2 specifically includes: A virtual vibration source layer is set between the transducer and the skull (approximately 5 mm above the skull) as the target surface for Rayleigh integral propagation. The medium grid index is calculated using physical coordinates. The x and y coordinates of the sampled CT image are used to obtain the highest point z in the image at these coordinates. The spatial coordinates of the virtual source layer are z+5 mm. With m virtual sources, there are m physical coordinates (x, y, z+5). Then, based on the specific scaling ratio of the CT image mapping grid in the previous step, the grid positions and index numbers of these m virtual sources in the grid are set.
[0027] In this embodiment, step S3 specifically includes: Assuming a uniform water medium between the actual transducer elements and the virtual source layer, the Green's function between all transducer array elements and the virtual source layer grid points is calculated using vectorization. Multiply by the array element area to obtain the complex transfer relationship between the two, and then quickly calculate the complex sound pressure distribution of the virtual source layer through the Rayleigh integration process.
[0028] Among them, the first formula of the Rayleigh-Sommerfeld diffraction integral gives the relationship between the sound pressure p(r) at the observation point r and the vibration velocity vn(rs) on the source surface S: ; In the formula, k is the wave number. Let and c be the density of the medium and the speed of sound, respectively; R = |r - rs| be the distance from the source point to the observation point; and dS be the infinitesimal area of the source surface. Under the assumption of a uniform water medium, this invention uses this formula to calculate the sound field propagation from the actual array element surface to the virtual source layer. The formula for the transfer matrix from each array element to the virtual point source is then obtained as follows: Where ΔS is the area of the array element, Let m be the distance from the array element m to the virtual point source j. This calculation is a purely analytical operation and is extremely fast.
[0029] In this embodiment, step S4 specifically includes: In the offline phase, each transducer unit is excited separately, and its steady-state complex sound pressure response within the target region (ROI) is calculated using a full-wave simulation method (such as k-Wave). The frequency response coefficients between each unit and each ROI point are extracted to construct a complex transfer matrix. The dimension of matrix H is defined as Nroi × Npres, where Nroi is the number of sensors (or grid points) within the ROI, and Npres is the number of preset array elements on the virtual source layer. Using the previously designed medium parameter grid, the sound field distribution in the ROI space is calculated under the activated state of each preset array element on the virtual source layer (under a defined reference signal, such as a pulse signal of 0.25MHz, 25Pa). The complex amplitude at the target frequency (i.e., the aforementioned frequency response coefficients) is extracted using Fourier transform and stored in matrix H.
[0030] In this embodiment, step S5 specifically includes: During the online operation phase, the excitation signal of the actual transducer is projected onto the virtual source layer using the Rayleigh integral formula (Tmn in step 3), and then multiplied with the pre-calculated transfer matrix (H in step 4) to quickly obtain the complex sound pressure distribution at any point in the target area. The entire process can be completed in milliseconds.
[0031] Combining the above steps, the process of online real-time prediction of the sound field distribution in the target area can be represented by the following concise chain matrix operation: ; ; Among them, E transducer Let H be the complex excitation vector of the N×1 dimensional transducer array; H be the pre-calculated P×M dimensional transfer matrix obtained from full-wave simulation; T be the M×N dimensional Rayleigh integral transfer matrix, whose elements are defined by the Green's function formula on the right side of the above equation, and are calculated quickly; P target The target region is a complex sound pressure vector of dimension P×1.
[0032] This invention also provides a real-time sound field simulation system based on physical layer decoupling to implement the aforementioned real-time sound field simulation method. The system includes: a personalized modeling module for establishing an acoustic parameter mesh containing a homogeneous medium region and a complex heterogeneous skull region based on patient CT images; a virtual source layer definition module for setting a virtual source layer near the interface between the homogeneous medium and the skull, and calculating its mesh index; an online analytical calculation module for calculating the first transfer matrix from the transducer elements to the virtual source layer online using the Rayleigh-Sommerfeld integral analytical formula; an offline pre-calculation module for pre-calculating the second transfer matrix from the virtual source layer to the target region offline using a full-wave numerical simulation method, and storing the result; and a real-time sound field synthesis module for multiplying the transducer excitation vector sequentially with the first and second transfer matrices to output the sound field distribution of the target region in real time.
[0033] Based on the above embodiments, the single prediction time of the present invention is about 0.05 seconds, which is more than 250 times faster than the traditional k-Wave method; and the normalized root mean square error is less than 0.6%, and the Pearson correlation coefficient is greater than 0.99; which can meet the millisecond-level response requirements for sound field feedback in therapeutic navigation.
[0034] For the foregoing embodiments, in order to simplify the description, they are all described as a series of actions. However, those skilled in the art should understand that this application is not limited to the described order of actions, because according to this application, some steps can be performed in other orders or simultaneously. Furthermore, those skilled in the art should also understand that the embodiments described in the specification are preferred embodiments, and the actions involved are not necessarily essential to this application.
[0035] The above embodiments describe the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The embodiments and descriptions in the specification are merely illustrative of the principles of the invention. Modifications and variations made by those skilled in the art without departing from the spirit and scope of the invention should be within the protection scope of the appended claims.
Claims
1. A real-time sound field simulation method based on physical layer decoupling, characterized in that, Includes the following steps: S1: Establishing an acoustic parameter mesh based on CT images, including homogeneous media regions and complex heterogeneous skull regions: S2: Define a virtual source layer in the uniform medium region between the transducer and the skull as a layered decoupling interface between analytical calculation and numerical calculation. S3: Using the Rayleigh-Sommerfeld integral analytical formula, the first transfer matrix from the actual transducer array elements to each virtual source point on the virtual source layer is calculated, thereby quickly obtaining the complex sound pressure distribution of the virtual source layer; S4: For complex and heterogeneous skull regions, a full-wave numerical simulation method is used to pre-calculate the frequency response coefficients from each virtual source point to each observation point in the target region, and to construct and store the second transfer matrix; S5: The excitation vector of the actual transducer is multiplied sequentially with the first transfer matrix and the second transfer matrix to synthesize the complete sound field distribution in the target area in real time.
2. The real-time sound field simulation method based on physical layer decoupling as described in claim 1, characterized in that, Step S1 includes the following sub-steps: S11: Acquire skull CT images and calculate bone porosity distribution based on Hounsfield units of the CT images; S12: Based on the bone porosity distribution, the sound velocity, density, and attenuation coefficient distribution of each voxel point in the skull are mapped to form an acoustic parameter grid that includes both homogeneous medium regions and complex heterogeneous skull regions.
3. The real-time sound field simulation method based on physical layer decoupling as described in claim 2, characterized in that, The density and velocity of sound at a voxel point are calculated using the following formulas: density =phi *ρw+ (1 - phi) *ρb; speed =phi * sw + (1 - phi) *sb; Where density is the density of voxel points; phi is the bone porosity; ρw is the water density; ρb is the bone density; sw is the water sound velocity; and sb is the bone sound velocity.
4. The real-time sound field simulation method based on physical layer decoupling as described in claim 1, characterized in that, Step S2 includes the following sub-steps: S21: Based on the spatial coordinates of the CT image, obtain the coordinates of the highest point on the outer surface of the skull on one side of the transducer. S22: Based on the coordinates of the highest point, offset a preset distance towards the transducer to form the plane where the virtual source layer is located; S23: Based on the scaling ratio of the acoustic parameter grid in step S1, convert the physical coordinates of each virtual source point on the virtual source layer into a grid index and record the number.
5. The real-time sound field simulation method based on physical layer decoupling as described in claim 1, characterized in that, The first transfer matrix is calculated as follows: ； in, is the element of the first transfer matrix; k is the wave number; ρ0 and c are the density and speed of sound of the homogeneous medium, respectively; ΔS is the area of the transducer element; Rmn is the distance from the nth transducer element to the mth virtual source point; and i is the imaginary unit.
6. The real-time sound field simulation method based on physical layer decoupling as described in claim 1, characterized in that, The full-wave numerical simulation method includes the following sub-steps: The acoustic parameter mesh obtained in step S1 is used as the simulation medium input; Each virtual source point on the virtual source layer is used as an excitation source in turn, and a reference excitation signal is set to perform steady-state sound field simulation. Extract the time-domain sound pressure signal at each observation point within the target area, and obtain the complex amplitude at the target frequency through Fourier transform, which is the frequency response coefficient; The frequency response coefficients between each virtual source point and each observation point are organized into a second transfer matrix H, with the dimension of the second transfer matrix being Nroi × Npres, where Nroi is the number of observation points in the target region and Npres is the number of virtual source points on the virtual source layer.
7. The real-time sound field simulation method based on physical layer decoupling as described in claim 1, characterized in that, The expression for the complete sound field distribution is: ； Among them, E transducer H is the complex excitation vector of the N×1 dimensional transducer array elements; H is the pre-computed transfer matrix of P×M dimensionality; P target Let be the P×1 dimensional complex sound pressure vector of the target region; T is the M×N dimensional Rayleigh integral transfer matrix.
8. The real-time sound field simulation method based on physical layer decoupling as described in claim 1, characterized in that, The homogeneous medium region is a water medium or a homogeneous coupling agent layer between the transducer and the skull, and the propagation of sound waves in this region follows the homogeneous medium wave equation; the complex heterogeneous skull region contains spatial variations in sound velocity, density, and porosity, resulting in sound wave scattering, refraction, and attenuation effects.
9. A real-time sound field simulation method based on physical layer decoupling as described in any one of claims 1 to 8, characterized in that, It also includes an adaptive adjustment step for the transducer excitation signal during the online phase: Receive transducer operating parameters or desired treatment focus location in real time; Based on the acoustic field distribution of the target region synthesized in step S5, the complex excitation vector E of each element of the transducer is optimized in reverse. transducer ; Repeat steps S3 to S5 to quickly update the sound field simulation results until the preset focusing accuracy or sound field distribution requirements are met.
10. A real-time sound field simulation system based on physical layer decoupling, used to implement the real-time sound field simulation method according to any one of claims 1 to 9, characterized in that, include: Personalized modeling module: used to create acoustic parameter meshes containing homogeneous media regions and complex heterogeneous skull regions based on patient CT images; Virtual source layer definition module: used to set up a virtual source layer near the interface between the homogeneous medium and the skull, and to calculate its mesh index; Online analytical calculation module: used to calculate the first transfer matrix from the transducer array elements to the virtual source layer online using the Rayleigh-Sommerfeld integral analytical formula; Offline pre-calculation module: used to pre-calculate and store the second transfer matrix from the virtual source layer to the target region offline using the full-wave numerical simulation method; Real-time sound field synthesis module: used to multiply the transducer excitation vector with the first transfer matrix and the second transfer matrix in sequence, and output the sound field distribution of the target area in real time.