In-vivo imaging method of radiation dose distribution of high-energy x-ray flash radiotherapy under limited view angle

Abstract

The application discloses a kind of limited view under high-energy X-ray FLASH radiotherapy radiation dose distribution in-vivo imaging method, belong to the technical field of three-dimensional radiation dose distribution, comprising: based on the array ultrasonic transducer of limited view acquisition acoustic signal;Based on acoustic signal, using improved weighted least square time reversal method reconstructs two-dimensional dose distribution chart;The reconstructed two-dimensional dose distribution chart is fused with anatomical structure chart, and the accumulated dose distribution chart is obtained.This embodiment is under the constraint of limited view, in the process of radiotherapy, the two-dimensional radiation dose distribution of the imaged tissue site is real-time and high-precision, so as to feedback subsequent dose planning;Radiation dose chart construction mainly uses improved least square time reversal method, through iterative optimization acoustic signal inversion, accurately draw the accumulated radiation dose distribution chart in target tissue.

Description

Technical Field

[0001] This invention belongs to the technical field of three-dimensional radiation dose distribution, specifically relating to an in vivo imaging method for high-energy X-ray FLASH radiotherapy radiation dose distribution under limited viewing angle. Background Technology

[0002] High-energy X-ray flash radiotherapy is an emerging method for treating malignant tumors while protecting surrounding normal tissues. To ensure the accuracy and safety of radiotherapy, it is essential to visualize the three-dimensional radiation dose distribution in the target tissue or organ in real time during treatment and to simultaneously calculate the cumulative dose. This real-time feedback mechanism is crucial for guiding and regulating the high-energy X-ray irradiation dose rate, aiming to deliver a sufficiently lethal dose to the tumor while maximizing the protection of surrounding normal tissues from radiation damage.

[0003] However, existing image-guided flash radiotherapy techniques face multiple challenges. Traditional computed tomography (CT) scans or deployment of radiation array sensors (such as thermal / optical stimulated dosimeters and scintillators) acquire and image transmitted X-rays, but these techniques are typically limited to measurements on the external surface of the human body. This prevents them from performing real-time three-dimensional dose imaging and dose accumulation calculations during radiotherapy, making it difficult to meet practical clinical needs. Although next-generation radiation detectors can be deployed inside the body, they typically cannot simultaneously provide anatomical information of the target tissue or real-time monitoring of the X-ray beam position, thus failing to provide effective real-time dose feedback guidance.

[0004] When using high-energy X-rays for Flash radiotherapy on target tissues or organs, the target tissues or organs absorb the X-ray beam energy, generating heat and forming an acoustic signal. The intensity of this acoustic signal is directly proportional to the energy or radiation dose absorbed by the target tissue or organ. Because the pulse width of medical linear accelerators used in high-energy X-ray radiotherapy is typically a few microseconds, the frequency of the acoustic signal is usually from tens to hundreds of kHz. Furthermore, the noise level in the radiation room is relatively high. To acquire reliable signals, a transducer is needed to collect the X-ray-induced acoustic signal (XA signal). After acoustic signal acquisition, an image reconstruction algorithm is used to construct the radiation dose distribution in the target tissue or organ. Then, a fusion matching algorithm is used to construct a human radiation dose distribution map, enabling real-time monitoring of human radiation dose during clinical radiotherapy and providing feedback to guide Flash radiotherapy.

[0005] Previous studies have validated the feasibility of ionizing radiation acoustic imaging (IRI) in dose distribution reconstruction under conventional radiotherapy conditions. Multi-channel, full-coverage transducer arrays are commonly used to achieve dose image inversion, but this research is often based on idealized full-view array sensors to acquire acoustic signals. Such full-view coverage is virtually impossible to achieve in real-world clinical settings, especially under the specific equipment setup required for high-energy X-ray flash radiotherapy. In actual clinical settings, due to limited space for equipment placement and fixed beam direction, ultrasound transducer arrays are typically deployed only on a limited detection boundary near the target area (e.g., linear or L-shaped arrays), resulting in a "limited field of view" sampling limitation. When photoacoustic tomography (PAT) uses this planar or linear detection geometry, the limited acoustic detection aperture introduces serious image quality problems. These problems include: streaking artifacts in the image, reduced image contrast, structural distortion, and decreased signal-to-noise ratio (SNR).

[0006] Therefore, it can be seen that existing photoacoustic image reconstruction algorithms have the following technical defects: 1. Traditional photoacoustic image reconstruction algorithms (e.g., filtered back projection or time-delay stacking) typically rely on full-view arrays to acquire data and obtain high-precision images. However, in actual clinical settings, due to limitations in equipment layout space and beam direction, transducer arrays can only achieve limited viewing angle sampling. This results in severe stripe artifacts, reduced contrast, and structural distortion in the reconstructed two-dimensional dose images, making it difficult to meet the accuracy requirements for dose verification in clinical settings.

[0007] 2. While existing X-ray acoustic imaging technology can acquire radiation dose distribution in real time by capturing sound wave signals, overcoming the saturation problem of traditional point dosimeters at ultra-high dose rates, traditional analytical algorithms or simple time-reversal methods struggle to effectively suppress severe artifacts when processing data acquired from limited viewing angles. Furthermore, these algorithms are poorly robust to acoustic non-uniformity in human tissues (e.g., changes in sound velocity), and cannot provide reliable real-time dose feedback guidance.

[0008] 3. From a limited perspective, the inversion of acoustic signals into high-precision dose distribution is a challenging inverse problem. Existing technologies struggle to introduce effective constraints at the physical model level to compensate for missing data information. Simply using image post-processing methods to correct the reconstructed image is insufficient to fundamentally recover the basic dose structure lost due to incomplete geometric information. Therefore, a model-driven iterative optimization system needs to be designed to incorporate the physical constraints of the wave equation into the reconstruction process to achieve high-precision and robust dose inversion. Summary of the Invention

[0009] The purpose of this invention is to address the aforementioned shortcomings in the prior art by providing an in vivo imaging method for high-energy X-ray FLASH radiotherapy radiation dose distribution under limited viewing angles, thereby solving the problem that the limited acoustic detection aperture in existing photoacoustic image reconstruction can introduce serious image quality issues.

[0010] To achieve the above objectives, the technical solution adopted by the present invention is as follows: An in vivo imaging method for high-energy X-ray FLASH radiotherapy radiation dose distribution under limited viewing angles, comprising the following steps: S1. Acquisition of acoustic signals using an array ultrasonic transducer with a limited field of view; S2. Based on acoustic signals, a two-dimensional dose distribution map is reconstructed using an improved weighted least squares time-reversal method. S3. The reconstructed two-dimensional dose distribution map is fused with the anatomical structure map to obtain the cumulative dose distribution map.

[0011] Furthermore, S1 specifically includes: Based on the derivation of the wave equation Time-array ultrasonic transducers in The sound pressure signal detected at the location is represented as follows: In the formula, express Time-array ultrasonic transducers in The sound pressure signal detected at the location; Speed ​​of sound; For Grüneisen parameters; The percentage of energy absorbed by the tissue and converted into heat; Tissue density; for Location in the organization at all times Dose deposition.

[0012] Furthermore, S2 specifically includes: The acoustic signal is amplified and filtered; the transmitted waveform is constructed through time-domain inversion operation; then, coherence weighting factor, geometric sensitivity weighting factor and channel reliability weighting factor are jointly adaptively controlled to coherently superimpose the echo signal to obtain the reconstructed sound pressure field; and the reconstructed two-dimensional dose distribution map is obtained based on the reconstructed sound pressure field.

[0013] Furthermore, in S2, a coherence weighting factor, a geometric sensitivity weighting factor, and a channel confidence weighting factor are used for joint adaptive modulation to coherently superimpose the echo signals and obtain the reconstructed sound pressure field, which is expressed as: in: In the formula, For the array ultrasonic transducer in position and time The reconstructed sound pressure field; This is the geometric sensitivity weighting factor; The total number of array elements of the deployed array transducers; For the first in the array transducer Each array element; As a coherence weighting factor; As a weighting factor for channel reliability; This is a time-domain reversal operation; From position To the Channel impulse response of each array element; For the first Channel signals of each array element; For the first The result of time-domain inversion operation on the channel signal of each array element; For the first The duration of the XA signal recorded by each array element.

[0014] Furthermore, in S2, the improved weighted least squares time reversal method iterative optimization objective function is: In the formula, This is the optimized dose distribution estimate; The dose distribution to be optimized; equal This indicates a reduction in the impact of noise channels; Let be the system matrix, representing the forward propagation operator from the dose distribution to the measured sound pressure; This is the measured sound pressure signal; Norm symbol; For regularization parameters; For anisotropic regularization operators; This is the gradient operator.

[0015] Furthermore, the coherence weighting factor Represented as: In the formula, For array element Backpropagation phase Standard deviation; The phase standard deviation threshold; It is a power exponent.

[0016] Furthermore, the geometric sensitivity weighting factor Represented as: In the formula, The probe normal vector; For the first The spatial position of each array element; This is the probe's directionality function.

[0017] Furthermore, the channel trust weighting factor Represented as: In the formula, The slope parameter represents the sigmoid function; Indicates the first The spectral entropy values ​​of each channel; This represents the spectral entropy threshold.

[0018] Furthermore, S2 also includes: improving the processing speed of the least squares time reversal method by using a dual-domain hierarchical k-space pseudospectral method and parallel acceleration of modern graphics processing units (GPUs), specifically as follows: The computational domain is decoupled into a high-fidelity source domain and a fast transmission domain at the physical and algorithmic levels; In the high-fidelity source domain, the complete nonlinear wave equation is solved only in a local small grid surrounding the radiotherapy target area to capture the sharp gradient of the dose boundary, and a spectrum-matched absorption layer is applied at the boundary to record the outgoing wave field without reflection absorption. In the fast transmission domain, the wave equation is solved using Green's function projection, and the signal from the array ultrasonic transducer is calculated by convolving the boundary wave field with the Green's function. In the formula, To be at the detector position place, time The detected sound pressure signal; The region of interest is the computational region of the high-fidelity source domain. For the boundary position place, time Recorded sound pressure of the outgoing wave field; For the time variable of the wave field at the boundary; It is the Green's function; The spatial position of the array ultrasonic transducer; The direction of the outer normal vector of the boundary; The boundary area element; Finally, the high-fidelity source domain and the fast transfer domain are executed independently and in parallel on different computing cores of modern graphics processing units (GPUs).

[0019] Furthermore, while the dual-domain hierarchical k-space pseudospectral method is iteratively executed, a spectral pruning operator is embedded to improve processing speed by forcing frequency domain cutoff and using a coarser grid. The spectral pruning operator is expressed as: In the formula, For the sound pressure field at wavenumber and time Frequency domain representation at; Pruning operator for spectrum; Wave number; This is the cutoff wave number.

[0020] The in vivo imaging method for high-energy X-ray FLASH radiotherapy radiation dose distribution under limited viewing angle provided by this invention has the following beneficial effects: This invention employs X-ray acoustic imaging technology, using a pulsed X-ray radiation beam of a few microseconds to irradiate the target tissue, generating an acoustic signal. This acoustic signal is then acquired by a transducer under limited viewing angle constraints. Since the peak value of the acoustic signal is proportional to the radiation dose absorbed by the target tissue, the cumulative radiation dose in the target tissue can be obtained through the acoustic signal. In order to significantly improve the accuracy of two-dimensional radiation dose imaging when X-ray irradiates target tissue under limited viewing angle, this invention utilizes an improved least squares time reversal method. Through iterative optimization, it effectively suppresses artifacts caused by incomplete data acquisition, accurately reconstructs the cumulative dose distribution map absorbed by the target tissue, and obtains the radiation dose absorption of different parts of the target tissue, thereby providing feedback to guide subsequent radiotherapy dose planning. Attached Figure Description

[0021] Figure 1 This is a flowchart of an in vivo imaging method for high-energy X-ray FLASH radiotherapy radiation dose distribution under limited viewing angles, as described in this embodiment. Detailed Implementation

[0022] The specific embodiments of the present invention are described below to enable those skilled in the art to understand the present invention. However, it should be understood that the present invention is not limited to the scope of the specific embodiments. For those skilled in the art, various changes are obvious as long as they are within the spirit and scope of the present invention as defined and determined by the appended claims. All inventions utilizing the concept of the present invention are protected.

[0023] This embodiment presents a limited-view, in vivo imaging method for high-energy X-ray FLASH radiotherapy radiation dose distribution. It utilizes an array transducer to acquire acoustic signals, and then uses a model-driven W-LSTR algorithm to generate a high-precision radiation dose deposition map. The constructed two-dimensional radiation dose distribution map is then fused with a target tissue anatomical structure map to obtain the cumulative dose distribution at various locations within the target tissue. This feedback guides radiotherapy dose planning and provides reference... Figure 1 Specifically, it includes the following: S1. Acquisition of acoustic signals using an array ultrasonic transducer with a limited field of view; Target tissue is irradiated with a microwave-grade pulsed X-ray beam. The target tissue absorbs the radiation dose and expands due to heat, generating acoustic signals that propagate outwards. These acoustic signals are acquired by an array of ultrasonic transducers deployed within a finite detection boundary (e.g., a linear or L-shaped array). Specifically, when periodic pulsed X-rays irradiate target tissue, the tissue absorbs the energy, causing instantaneous thermal expansion and generating sound wave signals. This physical phenomenon is known as the photoacoustic effect. The formula for calculating the conversion of X-ray energy absorbed by the target tissue into a sound pressure signal is as follows: In the formula, This is the initial sound pressure signal; It is the coefficient of volumetric thermal expansion; Speed ​​of sound; Specific heat capacity; The percentage of energy absorbed by the tissue and converted into heat; The mass attenuation coefficient of X-rays; Tissue density; This represents the X-ray flux.

[0024] X-ray flux and X-ray absorption are generally related to tissue material and vary with tissue density and mass attenuation coefficient. Assuming that... Dose deposition in tissue at time of time Then, from the wave equation, we can derive that in Ultrasonic transducers at all times The sound pressure signal detected at the location is represented as follows: In the formula, express Time-array ultrasonic transducers in The sound pressure signal detected at the location; For Grüneisen parameters; for Location in the organization at all times Dose deposition; Defined as: In the formula, It is the isothermal bulk modulus; Specific heat capacity.

[0025] Since the photoacoustic signal generated by the target tissue is proportional to the absorbed radiation dose, the photoacoustic signal at the sensor location detected by the ultrasonic array transducer can be used to infer the radiation dose absorbed inside the target tissue. This allows for real-time two-dimensional imaging of the radiation dose distribution map and cumulative radiation dose distribution map absorbed by the target tissue at different times during Flash radiotherapy.

[0026] S2. Based on the acoustic signal, a two-dimensional dose distribution map is reconstructed using an improved weighted least squares time-reversal method, which specifically includes the following: First, the acoustic signal is amplified and filtered to improve the signal-to-noise ratio and imaging accuracy. Subsequently, assuming an array transducer containing multiple elements, time-reversal (TR) beamforming technology first utilizes the channel impulse response (CIR) from the target region to each array element to construct the transmitted waveform through a time-domain reversal operation, and then coherently superimposes the echo signals. This time-domain reversal operation can be mathematically represented as: In the formula, This is a time-domain reversal operation; For the first Channel signals of each array element; For the first The result of time-domain inversion operation on the channel signal of each array element; For the first The duration of the XA signal recorded by each array element; Represents the first in the array transducer Each array element; Traditional time-reversal beamforming technology utilizes the channel impulse response from the target region to each array element, constructing the transmitted waveform through time-domain reversal and coherent superposition. Its positional... and time Reconstructed sound pressure field It can be represented as: Where N represents the total number of array elements of the deployed array transducers, and * represents the convolution operation. It is a weighting factor, that is, from the target position To the The channel impulse response (CIR) of each array element. The standard formula implicitly assumes that all signals are in phase superimposed in the time domain. However, under finite viewing angles, severe sidelobes and grating lobes lead to inconsistent phase of artifact signals at different spatial locations, resulting in superposition distortion. To overcome this limitation, this embodiment employs an improved weighted least squares time reversal method (W-LSTR), introducing a spatially variable weighting factor into the standard LSTR framework to achieve joint adaptive control of coherence, geometric sensitivity, and channel reliability, as detailed below: When the signal is correctly focused, the backpropagating signals of each array element should interfere in phase, while artifacts are caused by sidelobe leakage and phase disorder. Array element definition. The back propagation phase is Its standard deviation is Then, a weighted function is constructed: In the formula, As a coherence weighting factor; For array element Backpropagation phase Standard deviation; The phase standard deviation threshold; It is a power exponent used to adjust the response sensitivity of the weighting factor; When the signal is well coherently focused Approaching 0 Approaching 1; when in the region dominated by artifacts at a limited viewpoint, Almost zero, It can dynamically identify and suppress incoherent artifact signals, significantly improving the image contrast-to-noise ratio (CNR).

[0027] Under a limited field of view, voxels at different positions have different radiation solid angles to the array, resulting in uneven receiving sensitivity. Based on this, we define: In the formula, The probe normal vector; For the first The spatial position of each array element; This is the probe's directionality function; This is the geometric sensitivity weighting factor.

[0028] This geometric weight can compensate for the limitations of the "near-strong, far-weak" effect of limited viewing angles, improve the quantitative consistency of dose distribution, and ensure accurate inversion of dose in deep target areas.

[0029] In environments with strong ionizing radiation, individual channels may be susceptible to transient electromagnetic interference. Therefore, the signal purity is determined by calculating the spectral entropy of each channel, and weights are then constructed accordingly. In the formula, As a weighting factor for channel reliability; The slope parameter of the Sigmoid function is used to control the steepness of the weight changes; Indicates the first The spectral entropy value of each channel is used to measure the purity of the signal; This represents the spectral entropy threshold, used to distinguish between normal signals and interfered signals.

[0030] A higher spectral entropy indicates a larger proportion of noise, and the confidence weight automatically decreases, thus adaptively "shielding" the interfered channel and improving robustness and overall imaging stability.

[0031] After introducing the above weights, the improved reconstruction formula is: In the formula, For the array ultrasonic transducer in position and time The reconstructed sound pressure field; To start from the target location To the Channel impulse response of each array element; For the first Channel signals of each array element; For the first The result of time-domain inversion operation on the channel signal of each array element; For the first The duration of the XA signal recorded by each array element.

[0032] Its corresponding iterative optimization objective function is: In the formula, This is the optimized dose distribution estimate; The dose distribution to be optimized; equal This indicates a reduction in the impact of noise channels; For the system matrix; This is the measured sound pressure signal; Norm symbol; For regularization parameters; For anisotropic regularization operators, sharp edges are allowed in the vertical direction, while strong smoothing constraints are applied in the parallel direction; This is the gradient operator.

[0033] By dynamically balancing signal coherence, geometric sensitivity, and channel reliability during the iteration process, the W-LSTR algorithm effectively suppresses stripe artifacts and sidelobe interference caused by limited viewing angles, significantly improves the signal-to-noise ratio (SNR), contrast-to-noise ratio (CNR), and quantitative accuracy of dose inversion of the reconstructed image, and exhibits better stability and robustness.

[0034] In some embodiments, considering that W-LSTR is an iterative algorithm with higher computational complexity than traditional analytical algorithms, this invention performs in-depth optimization based on the traditional k-space method to meet the real-time feedback requirements of FLASH radiotherapy. Although the k-space pseudospectral method calculates the wavefield spatial derivative through Fast Fourier Transform (FFT), which has higher accuracy and smaller dispersion error compared to the Finite Difference Method (FDTD), it has significant efficiency bottlenecks in FLASH radiotherapy scenarios. First, FFT must be run on a regular Cartesian grid, resulting in global grid redundancy; in order to cover the propagation path from the target area to the detector, the simulation domain must include the entire tissue cross-section (e.g., 40cm × 40cm), while the actual irradiation area is only about 5cm × 5cm. This means that more than 98% of the computational area is empty space, resulting in a huge waste of computational power. Second, although k-space allows for a larger time step, high-frequency components are still limited by the Nyquist sampling theorem, making it difficult to accelerate sufficiently.

[0035] Based on this, this embodiment employs a dual-domain hierarchical k-space pseudospectral method and parallel acceleration by modern graphics processing units (GPUs) to improve the processing speed of the least squares time reversal method. Specifically: We propose a dual-domain layered propagation model, which decouples the computational domain into two regions at the physical and algorithmic levels: a high-fidelity source domain and a fast transport domain.

[0036] In the high-fidelity source domain, the complete nonlinear wave equation is solved only in a local small grid surrounding the radiotherapy target area (PTV). The high precision advantage of k-space is used to accurately capture the sharp gradient of the dose boundary, and a spectral matching absorption layer (Spectral PML) is applied at the boundary to absorb the emitted wave field without reflection. In the fast transmission domain, instead of solving the wave equation on a voxel-by-voxel basis, Green's function projection or the ray-acoustic approximation is used to calculate the detector signal by convolving the boundary wave field with the Green's function. In the formula, To be at the detector position place, time The detected sound pressure signal; The region of interest is the computational region of the high-fidelity source domain. For the boundary position place, time Recorded sound pressure of the outgoing wave field; For the time variable of the wave field at the boundary; It is the Green's function; The spatial position of the array ultrasonic transducer; The direction of the outer normal vector of the boundary; The boundary area element; This design transforms the original three-dimensional volume integral into a two-dimensional area integral, achieving a computational speedup of approximately 10-50 times and significantly improving the system's real-time performance.

[0037] Meanwhile, this invention introduces spectral pruning technology to further improve efficiency. In FLASH radiotherapy, the X-ray-induced acoustic energy is mainly concentrated in the low frequency range. To avoid high-frequency redundant calculations required for numerical stability, the algorithm embeds a spectral pruning operator in the k-space iteration: In the formula, For the sound pressure field at wavenumber and time Frequency domain representation at; Pruning operator for spectrum; Wave number (spatial frequency); This is the cutoff wavenumber; high-frequency components exceeding this value are filtered out.

[0038] This embodiment introduces spectral pruning, which, by forcing frequency domain cutoff, allows the use of a coarser grid without introducing aliasing errors. For example, when the grid step size Δx is increased by a factor of 2, the number of voxels is reduced by a factor of 8, and the time step size Δt is increased by a factor of 2, the theoretical total speedup can reach 16 times.

[0039] Finally, leveraging the high parallelism of GPUs, multi-layer parallel computation is achieved within a dual-domain model and spectral pruning framework. The high-fidelity source domain and fast transport domain are executed independently and in parallel on different GPU cores, while spectral pruning and dispersion correction are implemented in the frequency domain using parallel FFT. Experimental results show that this method achieves an overall computational speed improvement of approximately 20 times while maintaining imaging accuracy, providing real-time dose imaging support for FLASH radiotherapy.

[0040] S3. The reconstructed two-dimensional dose distribution map is fused with the anatomical structure map to obtain the cumulative dose distribution map; Specifically, an anatomical map of the target tissue is constructed using computed tomography (CT) technology. The constructed two-dimensional radiation dose distribution map is then fused and matched with the anatomical map of the target tissue to calculate and obtain the cumulative dose distribution of each part of the target tissue in real time.

[0041] This invention employs an improved Weighted Least Squares Time Reversal (W-LSTR) algorithm, treating dose reconstruction as a model-driven iterative optimization process. It innovatively introduces a triple adaptive weighting mechanism: a coherence weighting factor to dynamically identify and suppress incoherent artifact signals; a geometric sensitivity weighting factor to compensate for the "near-strong, far-weak" effect under limited viewing angles; and a channel reliability weighting factor to adaptively shield channels affected by electromagnetic interference. This spatial-coherence joint adaptive control effectively compensates for information loss due to incomplete data at the physical model level, significantly suppressing artifacts. Quantitative evaluation shows that the contrast-to-noise ratio (CNR) of the image reconstructed by the W-LSTR algorithm is up to 4.2 dB higher than that of traditional algorithms, greatly improving the imaging quality and quantitative accuracy of the constructed radiation dose distribution map.

[0042] Secondly, based on the physical model constraint mechanism, clinical robustness and generalization are significantly improved. Existing methods for improving image quality often employ deep learning for post-processing of reconstructed images. When images are severely distorted due to limited viewing angles, such post-processing methods struggle to recover lost physical structural information. The W-LSTR method of this invention directly optimizes the original acoustic signal at the inverse problem solving level by introducing regularization terms and physical constraints (such as dose nonnegativity and anisotropic regularization operators). The anisotropic regularization operator allows sharp dose edges in the direction perpendicular to the detector and applies strong smoothing constraints in the parallel direction, maintaining the clarity of dose distribution boundaries while effectively suppressing noise. Compared to deep learning post-processing that relies solely on image features, this approach obtains more robust physical information, ensuring higher robustness and generalization of the dose reconstruction results when facing complex acoustically inhomogeneous media (such as changes in sound velocity in human tissue).

[0043] Furthermore, this invention achieves an effective balance between high precision and real-time feedback. Although W-LSTR, as an iterative reconstruction algorithm, has relatively high computational complexity, this invention addresses this challenge through systematic technical optimization. First, by focusing the application on two-dimensional imaging, the computational resource requirements for solving the wave equation are reduced. Second, an innovative dual-domain hierarchical propagation model is proposed, decoupling the computational domain into a high-fidelity source domain and a fast transport domain: the complete wave equation is solved in a small local grid surrounding only the radiotherapy target area to accurately capture the dose boundary; in the transport domain, Green's function projection or ray acoustic approximation is used to transform the three-dimensional volume integral into a two-dimensional area integral, achieving a computational speedup of approximately 10-50 times. Furthermore, spectral pruning technology is introduced, using a coarser grid by forcing frequency domain cutoff without introducing aliasing errors, theoretically achieving a speedup of up to 16 times. Finally, combined with GPU parallel computing, multi-layer parallel execution is achieved within the dual-domain model and spectral pruning framework. Experimental results show that this method can improve the overall calculation speed by about 20 times while maintaining imaging accuracy, ensuring that high-precision dose reconstruction can be achieved while achieving a clinically acceptable near real-time imaging speed, meeting the need for real-time feedback guidance on cumulative dose distribution during Flash radiotherapy.

[0044] Although specific embodiments of the invention have been described in detail with reference to the accompanying drawings, this should not be construed as limiting the scope of protection of this patent. Various modifications and variations that can be made by a person skilled in the art without inventive effort within the scope described in the claims still fall within the scope of protection of this patent.

Claims

1. A method for in vivo imaging of high-energy X-ray FLASH radiotherapy radiation dose distribution under limited viewing angles, characterized in that, Includes the following steps: S1. Acquisition of acoustic signals using an array ultrasonic transducer with a limited field of view; S2. Based on acoustic signals, a two-dimensional dose distribution map is reconstructed using an improved weighted least squares time-reversal method. S3. The reconstructed two-dimensional dose distribution map is fused with the anatomical structure map to obtain the cumulative dose distribution map.

2. The in vivo imaging method for high-energy X-ray FLASH radiotherapy radiation dose distribution under limited viewing angle according to claim 1, characterized in that, S1 specifically includes: Based on the derivation of the wave equation Time-array ultrasonic transducers in The sound pressure signal detected at the location is represented as follows: In the formula, express Time-array ultrasonic transducers in The sound pressure signal detected at the location; Speed ​​of sound; For Grüneisen parameters; The percentage of energy absorbed by the tissue and converted into heat; Tissue density; for Location in the organization at all times Dose deposition.

3. The in vivo imaging method for high-energy X-ray FLASH radiotherapy radiation dose distribution under limited viewing angle according to claim 2, characterized in that, S2 specifically includes: The acoustic signal is amplified and filtered; the transmitted waveform is constructed through time-domain inversion operation; then, coherence weighting factor, geometric sensitivity weighting factor and channel reliability weighting factor are jointly adaptively controlled to coherently superimpose the echo signal to obtain the reconstructed sound pressure field; and the reconstructed two-dimensional dose distribution map is obtained based on the reconstructed sound pressure field.

4. The in vivo imaging method for high-energy X-ray FLASH radiotherapy radiation dose distribution under limited viewing angle according to claim 3, characterized in that, In S2, a joint adaptive modulation of the coherence weighting factor, geometric sensitivity weighting factor, and channel confidence weighting factor is used to coherently superimpose the echo signals to obtain the reconstructed sound pressure field, which is expressed as: in: In the formula, For the array ultrasonic transducer in position and time The reconstructed sound pressure field; This is the geometric sensitivity weighting factor; The total number of array elements of the deployed array transducers; For the first in the array transducer Each array element; As a coherence weighting factor; As a weighting factor for channel reliability; This is a time-domain reversal operation; From position To the Channel impulse response of each array element; For the first Channel signals of each array element; For the first The result of time-domain inversion operation on the channel signal of each array element; For the first The duration of the XA signal recorded by each array element.

5. The in vivo imaging method for high-energy X-ray FLASH radiotherapy radiation dose distribution under limited viewing angle according to claim 4, characterized in that, In S2, the improved weighted least squares time reversal method iterative optimization objective function is: In the formula, This is the optimized dose distribution estimate; The dose distribution to be optimized; equal This indicates a reduction in the impact of noise channels; Let be the system matrix, representing the forward propagation operator from the dose distribution to the measured sound pressure; This is the measured sound pressure signal; Norm symbol; For regularization parameters; For anisotropic regularization operators; This is the gradient operator.

6. The in vivo imaging method for high-energy X-ray FLASH radiotherapy radiation dose distribution under limited viewing angle according to claim 4, characterized in that, Coherence weighting factor Represented as: In the formula, For array element Backpropagation phase Standard deviation; The phase standard deviation threshold; It is a power exponent.

7. The in vivo imaging method for high-energy X-ray FLASH radiotherapy radiation dose distribution under limited viewing angle according to claim 4, characterized in that, Geometric sensitivity weighting factor Represented as: In the formula, The probe normal vector; For the first The spatial position of each array element; This is the probe's directionality function.

8. The in vivo imaging method for high-energy X-ray FLASH radiotherapy radiation dose distribution under limited viewing angle according to claim 4, characterized in that, Channel reliability weighting factor Represented as: In the formula, The slope parameter represents the sigmoid function; Indicates the first The spectral entropy values ​​of each channel; This represents the spectral entropy threshold.

9. The in vivo imaging method for high-energy X-ray FLASH radiotherapy radiation dose distribution under limited viewing angle according to claim 3, characterized in that, The S2 further includes: improving the processing speed of the least squares time reversal method by using a dual-domain hierarchical k-space pseudospectral method and parallel acceleration of modern graphics processing units (GPUs), specifically as follows: The computational domain is decoupled into a high-fidelity source domain and a fast transmission domain at the physical and algorithmic levels; In the high-fidelity source domain, the complete nonlinear wave equation is solved only in a local small grid surrounding the radiotherapy target area to capture the sharp gradient of the dose boundary, and a spectrum-matched absorption layer is applied at the boundary to record the outgoing wave field without reflection absorption. In the fast transmission domain, the wave equation is solved using Green's function projection, and the signal from the array ultrasonic transducer is calculated by convolving the boundary wave field with the Green's function. In the formula, To be at the detector position place, time The detected sound pressure signal; The region of interest is the computational region of the high-fidelity source domain. For the boundary position place, time Recorded sound pressure of the outgoing wave field; For the time variable of the wave field at the boundary; It is the Green's function; The spatial position of the array ultrasonic transducer; The direction of the outer normal vector of the boundary; The boundary area element; Finally, the high-fidelity source domain and the fast transfer domain are executed independently and in parallel on different computing cores of modern graphics processing units (GPUs).

10. The in vivo imaging method for high-energy X-ray FLASH radiotherapy radiation dose distribution under limited viewing angle according to claim 9, characterized in that, While the bi-domain hierarchical k-space pseudospectral method is iteratively executed, a spectral pruning operator is embedded to improve processing speed by forcing frequency domain cutoff and using a coarser grid. The spectral pruning operator is expressed as: In the formula, For the sound pressure field at wavenumber and time Frequency domain representation at; Pruning operator for spectrum; Wave number; This is the cutoff wave number.

Images