X-ray thermoacoustic imaging reconstruction method based on fourier neural operator

Abstract

This invention discloses an X-ray thermoacoustic imaging reconstruction method based on Fourier neural operators, belonging to the field of X-ray imaging technology. This invention rapidly constructs a multimodal X-ray thermoacoustic dataset covering different sampling conditions; fully utilizes the frequency domain global feature extraction advantages of Fourier operators to effectively improve the accuracy of forward physics model simulation and dose prediction; the dual-head Fourier structure effectively enhances the stepwise correction effect of absorption coefficient and dose distribution, making the reconstructed image clearer and more quantitatively accurate; lightweight training and recursive parameter updates reduce computational resource consumption, enabling the system to achieve real-time or near-real-time imaging capabilities, fully realizing a closed-loop system from initial data simulation, network reconstruction to dose distribution verification, and exhibiting high clinical adaptability.

Description

Technical Field

[0001] This invention relates to the field of X-ray imaging technology, specifically to an X-ray thermoacoustic imaging reconstruction method based on Fourier neural operators. Background Technology

[0002] X-ray imaging technology has been a core tool for medical diagnosis for more than a century, but it has inherent limitations, such as the need for strict alignment of the X-ray detector with the radiation source and the potential carcinogenic risk of ionizing radiation.

[0003] In recent years, X-ray induced acoustic computed tomography (XACT), as an emerging technology, has broken through the limitations of traditional X-ray imaging by using short-pulse X-rays to excite tissues to generate acoustic signals for imaging. Based on the linear relationship between X-ray dose and X-ray induced acoustic (XA) signals, XACT offers new possibilities for real-time monitoring of radiotherapy and quantitative imaging of deep tissues. However, existing XACT reconstruction methods still face significant challenges.

[0004] Filtered back projection (FBP) is widely used due to its high computational efficiency. However, because the acoustic signals excited by X-rays have complex positive and negative intensity time-series characteristics (such as multi-frequency oscillations and scattering interference), its radial projection summation mechanism cannot accurately quantify the dose distribution, resulting in significant noise in the reconstructed images, especially in low-dose scenarios where artifacts can account for more than 30%. Time-reversal methods, based on the wave equation to propagate the acoustic field signal backward, can partially suppress noise, but their performance is highly dependent on the completeness of the sensor layout. In practical imaging scenarios such as sparse sampling (e.g., a 32-channel ring array) or limited viewing angles (<180° coverage), the resolution of the reconstructed images decreases by more than 40%, failing to meet the needs of fine diagnosis such as tumor boundary delineation.

[0005] To improve imaging quality in complex scenarios, iterative reconstruction methods (such as variational regularization and model basis optimization) have gradually become a research hotspot. For example, the combined optimization framework of absorption coefficient and dose distribution can reduce the quantitative error to within 5%. However, such methods require repeated solving of the thermoacoustic wave equation to match the measured signal, and a single forward simulation takes several minutes (based on the finite element method), resulting in an overall reconstruction time of more than 1 hour, which is difficult to meet the real-time requirements of clinical practice (such as intraoperative navigation requiring ≤5 seconds / frame). To overcome this bottleneck, researchers have attempted to introduce deep learning into XACT. Deep learning methods have been applied to XACT reconstruction. For example, end-to-end mapping based on U-Net can quickly generate images, but it requires a large amount of labeled data and is sensitive to noise (SSIM decreases significantly at low doses), and its generalization ability is insufficient. Hybrid architecture methods (such as PixelInterp-CNN) partially alleviate the overfitting problem by combining temporal signal encoding with the initial backprojected image, but still require manual design of regularization terms and cannot adaptively optimize the physical model parameters.

[0006] In recent years, neural operators (such as Fourier neural operators, FNO) have demonstrated revolutionary advantages in solving partial differential equations (PDEs). FNOs, through Fourier spatial parameterization of the integral kernel, can adapt to arbitrary spatial resolution and geometry, avoiding the repetitive discretization calculations of traditional numerical methods (such as the finite element method). In photoacoustic imaging, FNOs have successfully replaced finite element solvers, increasing the speed of light transmission simulation by 30 times while maintaining an error of <3%. However, differences in the physical model of XACT lead to direct transfer failure. For example, existing FNO architectures do not model the secondary electron scattering effect of X-rays in the bone region, resulting in a dose distribution estimation bias >15%. Furthermore, traditional methods process noise suppression, dose correction, and absorption coefficient inversion step-by-step, neglecting their coupling relationship, leading to amplified cumulative errors. Although iterative reconstruction combined with FNOs can partially alleviate computational pressure, its high memory footprint (>12 GB) still limits its deployment in embedded medical devices. Therefore, a technical solution that deeply integrates physical models and deep learning is needed to achieve efficient and high-precision XACT imaging through targeted network structure design, joint optimization framework and lightweight strategy, so as to promote its practical application in clinical practice. Summary of the Invention

[0007] The purpose of this invention is to provide an X-ray thermoacoustic imaging reconstruction method based on Fourier neural operators, which aims to solve the problems of low iterative reconstruction efficiency, insufficient dose correction accuracy, and high noise sensitivity in existing X-ray thermoacoustic imaging (XACT) technology. By integrating a dual-head (dual-branch) dual-head FNO model structure with an iterative optimization framework, efficient and high-precision XACT image reconstruction and dose correction are achieved, while reducing the computational resource requirements and promoting its application in clinical real-time imaging.

[0008] To achieve the above objectives, the present invention provides the following technical solution: an X-ray thermoacoustic imaging reconstruction method based on Fourier neural operators, comprising at least the following steps:

[0009] S1: Construct an X-ray thermoacoustic data simulation and multimodal dataset, which includes an initial sound pressure image and the corresponding time-domain signal;

[0010] S2: Construct a two-headed Fourier neural operator model, also known as a two-headed FNO model, and input the signal simulated by the two-headed FNO model into the joint optimization network;

[0011] S3: Joint optimization network construction and training: The simulated ultrasound signal datasets from different sampling perspectives and the initial sound pressure images are fed into the joint optimization network for training;

[0012] S4: Perform high-quality image reconstruction of real radiotherapy signals. Put the X-ray thermoacoustic signals collected in the radiotherapy environment into the trained dual-head FNO model to reconstruct high-quality X-ray thermoacoustic images, that is, obtain the reconstructed sound pressure distribution map.

[0013] S5: Convert the reconstructed sound pressure distribution into a dose distribution according to the corresponding formula, and compare the dose distribution with the original radiotherapy plan.

[0014] Furthermore, the construction of the X-ray thermoacoustic data simulation and multimodal dataset includes at least the following steps:

[0015] First, determine the tools and environment. In the MATLAB environment, use the K-Wave toolbox to simulate X-ray induced thermoacoustic signals under different sampling conditions.

[0016] Secondly, signal simulation is performed to conduct detailed simulation of X-ray induced thermoacoustic signals under different sampling conditions. The different sampling conditions include at least the number of sensors, the shape of the sensor array, and the sampling angle.

[0017] The data is then processed by inputting clinical CT images and dividing the tissue into different regions using HU value segmentation. These different regions include at least bone, soft tissue, and tumor.

[0018] Then, a preset threshold is set, and corresponding tissue density and Grunison coefficient are assigned to each region to determine the dose-sound pressure conversion parameters.

[0019] Subsequently, the initial sound pressure distribution is generated. According to the physical conversion formula, the radiation dose distribution in the radiotherapy plan is combined with the tissue characteristic parameters and converted into the initial sound pressure distribution. The initial sound pressure distribution is then stored as a grayscale image with a fixed resolution.

[0020] Subsequently, sound wave propagation simulation was performed. Under the set ultrasonic sensor parameters, K-Wave was used to simulate sound wave propagation, and Gaussian white noise was added to generate ultrasonic simulation signal data under different signal-to-noise ratio conditions.

[0021] Finally, the dataset was constructed. Based on the above steps, an X-ray thermoacoustic data simulation and multimodal dataset was built and divided into a training set and a test set in a 7:3 ratio.

[0022] Furthermore, the dual-head FNO model, by employing a dual-head multilayer Fourier neural operator structure, extracts global information from the input signal in the frequency domain, simulating the signal propagation effect caused by X-ray dose deposition in different tissues;

[0023] The output of the dual-head FNO model can be used as both a simulated sensor time-domain signal and a preliminary estimate of the dose distribution, which can then be used as a basis for data alignment during subsequent iterative updates.

[0024] Furthermore, the dual-head FNO model is used to simultaneously realize acoustic wave forward and inverse updates on the same operator skeleton. The overall architecture of the dual-head FNO model adopts a "dimensionality increase – multi-layer Fourier layer – fusion – projection" process, but in the final projection, two independent projection branches are used, namely the forward projection branch and the inverse projection branch.

[0025] Dimensional upscaling refers to mapping the input image to a higher-dimensional channel feature space through a dimensional upscaling layer (fully connected or 1×1 convolution);

[0026] Introduce multiple Fourier layers (e.g., four layers), each layer truncates low-frequency modes in the frequency domain and applies a learnable kernel integral operator;

[0027] The inverse transform is then fused with the parallel linear mapping branch to capture global propagation characteristics;

[0028] Furthermore, the input to the forward branch is the current initial sound pressure image;

[0029] The simulated time-domain signal tensor, which corresponds to the number of sensor array channels and time sampling points, is output through a projection layer (1×1 convolution) to replace the traditional finite element / Monte Carlo forward modeling solution.

[0030] The input to the inversion branch is the current image and its simulation residual, which are then mapped and stitched together to form an H×W×2 channel input.

[0031] The function of the inversion branch is to directly learn the "residual → image correction" operator, replacing the traditional CNN.

[0032] Furthermore, S3 includes at least the following steps:

[0033] The ultrasonic signals and FBP-reconstructed X-ray thermoacoustic images in the ultrasonic simulation signal dataset are used as inputs.

[0034] Furthermore, the labels for the forward branch are real time-domain signals generated by the K-Wave toolbox (or other numerical simulations), while the labels for the inverse branch are used to calculate the difference between the real pressure distribution and the initial reconstructed image (or the real absorption coefficient minus the initial value).

[0035] Task-weighted loss is used:

[0036] L=αLfwd+(1−α)Linv

[0037] Where Lfwd is the time-domain signal prediction error (such as RMSE), Linv is the image update error (such as a combination of RMSE or SSIM), and α is the balance coefficient;

[0038] The Adam (or AdamW) optimization algorithm based on backpropagation iterates over network parameters;

[0039] By combining the early stopping mechanism with validation set monitoring, the optimal parameters are selected.

[0040] Furthermore, S4 includes at least the following steps:

[0041] The joint optimization network in S3 was applied to the reconstruction of the actual acquired X-ray thermoacoustic signals;

[0042] The ultrasound signals acquired in the actual radiotherapy environment are preprocessed, and a low-resolution initial image is obtained by using simple back projection (such as time reversal algorithm).

[0043] The preprocessed data is input into a joint optimization network, and a dual-head FNO model is used to reconstruct high-resolution and high signal-to-noise ratio X-ray thermoacoustic images to obtain the reconstructed sound pressure distribution map.

[0044] The reconstructed images outperform traditional methods in detail restoration, noise suppression, and quantitative metrics such as absorption coefficient and dose accuracy.

[0045] Furthermore, S5 includes at least the following steps:

[0046] The reconstructed sound pressure distribution map is converted into the corresponding X-ray dose deposition distribution map according to the known physical conversion formula;

[0047] Then, the RT PLAN dose distribution (DICOM format) output from the clinical radiotherapy system is imported and rigidly / affinely registered with the reconstructed two-dimensional dose distribution on the same anatomical plane to ensure a one-to-one correspondence between pixels and voxels.

[0048] The two-dimensional reconstructed sound pressure distribution image is converted into a dose distribution image using analytical formulas and tissue characteristic parameters;

[0049] Calculate the DVH curves of reconstructed dose versus planned dose within the target area and the contours of major organs (such as PTV, bladder, rectum, etc.);

[0050] Based on statistical differences and DVH deviation, the consistency between the reconstructed dose and the clinical RT PLAN was evaluated, and the usability of this method in clinical dose monitoring and feedback was verified.

[0051] Compared with the prior art, the beneficial effects of the present invention are:

[0052] This invention rapidly constructs a multimodal X-ray thermoacoustic dataset covering different sampling conditions; fully utilizes the frequency domain global feature extraction advantages of FNO to effectively improve the accuracy of forward physics model simulation and dose prediction; the two-stage iterative optimization structure effectively enhances the gradual correction effect of absorption coefficient and dose distribution, making the reconstructed images clearer and the quantification more accurate; lightweight training and recursive parameter updates reduce computational resource consumption, enabling the system to have real-time or near-real-time imaging capabilities and fully realize a closed-loop system from initial data simulation, network reconstruction to dose distribution verification, with high clinical adaptability. Attached Figure Description

[0053] To more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings used in the description of the embodiments 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 these drawings without creative effort.

[0054] Figure 1 This is a schematic diagram of the initial pressure in the dataset of this invention;

[0055] Figure 2 This is a schematic diagram of the arrangement of the 128 sensors of the present invention;

[0056] Figure 3 This is a flowchart of the FNO-based process architecture of the present invention;

[0057] Figure 4 This is a diagram of the internal framework of the FNO module of the present invention;

[0058] Figure 5 This is a schematic diagram of the human body X-ray thermoacoustic reconstruction based on 128 sensors according to the present invention. Detailed Implementation

[0059] The technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments.

[0060] An X-ray thermoacoustic imaging reconstruction method based on Fourier neural operators includes at least the following steps:

[0061] S1: Construct an X-ray thermoacoustic data simulation and multimodal dataset, which includes an initial sound pressure image and the corresponding time-domain signal;

[0062] The process of obtaining the initial sound pressure image in S1 is as follows: Figure 1As shown; Tissue characteristic parameter extraction and initial sound pressure generation were performed by acquiring clinical CT images in DICOM format with a resolution of 512×512 pixels and a slice thickness of 1 mm. Using the grayscale values ​​of the CT images, the image regions were divided into three categories according to HU values: bone region: HU > 300, soft tissue region: 50 ≤ HU ≤ 300, and tumor region: HU < 50. Corresponding tissue characteristic parameters were assigned to each region: bone: tissue density ρ = 1.85 g / cm³, Glunnison coefficient Γ = 0.15, soft tissue: ρ = 1.05 g / cm³, Γ = 0.25, tumor: ρ = 1.10 g / cm³, Γ = 0.30.

[0063] Based on the radiation dose distribution D(r) and local X-ray absorption coefficient μ provided in the radiotherapy plan (r) (selected or estimated based on organization type), using the following conversion formula:

[0064] p0(r)=ρΓηthD(r)

[0065] The radiation dose is converted into an initial sound pressure distribution.

[0066] The obtained sound pressure distribution was resampled into a 128×128 grayscale image and used as label data for subsequent use.

[0067] Specifically, the K-Wave toolbox (MATLAB version) is used to simulate sound wave propagation. The sensor array is set to a circular array, with the number of sensors selectable as 32, 64, or 128; the array radius is set to 40 mm; the ultrasonic sensor parameters are: center frequency 1 MHz, bandwidth 80%, sampling rate 8 MHz; the medium sound velocity is set to 1540 m / s. Gaussian white noise is added during the simulation to simulate signal-to-noise ratios of 15 dB, 20 dB, 25 dB, and 30 dB. Each set of data includes the initial sound pressure level image and the corresponding ultrasonic time-domain signal (approximately 512 sampling points). It is recommended to generate approximately 10,000 sets of data in total, divided into training and test sets in a 7:3 ratio.

[0068] S2: Construct a two-headed Fourier neural operator model, also known as a two-headed FNO model, and input the signal simulated by the two-headed FNO model into the joint optimization network;

[0069] The initial sound pressure image from the X-ray thermoacoustic data simulation and the multimodal dataset is used as the input signal;

[0070] A multi-layer Fourier neural operator structure is used to extract global information from the input signal in the frequency domain to simulate the signal propagation effect caused by X-ray dose deposition in different tissues;

[0071] The output of the FNO module can be either a simulated sensor time-domain signal or a preliminary estimate of the dose distribution, which can then be used as a basis for data alignment during subsequent iterative updates.

[0072] Specifically, the construction of the dual-head FNO joint optimization network in step S2 is as follows: Figure 3 and Figure 4 As shown, the initial sound pressure image a(x) (128×128) generated in S1 is input. First, a fully connected layer elevates the input from a low-dimensional space to a high-dimensional representation to facilitate subsequent global frequency domain processing. A four-layer Fourier neural operator is constructed, with each layer designed as follows: a two-dimensional Fourier transform F is performed on the input to transfer spatial domain features to the frequency domain; in the frequency domain, only low-frequency modes are retained (e.g., 64 low-frequency modes are retained, covering 0.5–5 MHz), and a linear transform R is applied to the low-frequency modes to enhance the globally effective information; an inverse Fourier transform is performed on the processed frequency domain data. This restores the input to the spatial domain; simultaneously, a parallel branch is set up to directly apply a 1×1 convolutional layer to the input. The upper branch output is further enhanced with low-frequency features through two 1×1 convolutional layers (Wt1 and Wt2) and then fused with the lower branch to obtain the final output of this layer. After processing through four Fourier layers, the output v4(x) is mapped back to the target dimension through two 1×1 convolutional layers to generate a simulated output of the estimated X-ray dose deposition distribution or sensor time-domain signal, providing initial mapping results for subsequent iterative optimization. S2 can be divided into the following steps:

[0073] 2.1 Input Preprocessing

[0074] 1. Forward input: Initial sound pressure distribution image

[0075]

[0076] As input to the forward branch, it is a single channel.

[0077] Inversion input: In the k-th round of iterative reconstruction:

[0078] Utilize the trained FNO_forward on the current estimated image Perform a forward modeling to obtain the analog signal. .

[0079] With measured signal Calculate residuals

[0080]

[0081] right Perform pixel interpolation or a small convolution to map it to the same spatial resolution as the image, and obtain... .

[0082] Compared with the current estimated image By splicing the channels together, a dual-channel input tensor is formed.

[0083]

[0084] 2.2 The network backbone adopts a unified backbone network to process the above two types of inputs: 1. Upgrading layer

[0085] For input tensor (Single channel during forward modeling, dual channel during inversion) Applications Convolution increases the number of channels to Output

[0086]

[0087] 2. Four Fourier layers for each layer ,Will Send in:

[0088] Frequency domain branch

[0089] Two-dimensional Fourier transform: .

[0090] Keep only the previous one A low-frequency mode filters out high frequencies.

[0091] Applying the learnable kernel operator : .

[0092] Inverse Fourier Transform: .

[0093] Two series Convolutional layer And GELU activation, to obtain .

[0094] Direct branch to the original input Direct application convolution ,get .

[0095] Fusion

[0096]

[0097] 2.3 Output Head

[0098] Forward branch

[0099] right application Convolution maps channels to (Sensor Channel) (Time step) dimension, output

[0100]

[0101] Inversion branch

[0102] For the same Apply another set Convolution maps channels to single-channel correction increments.

[0103]

[0104] S3: Joint optimization network construction and training: The simulated ultrasound signal datasets from different sampling perspectives and the initial sound pressure images are fed into the joint optimization network for training;

[0105] S3 specifically includes the following steps:

[0106] 3.1 Training Data Preparation

[0107] Actual mission

[0108] Dataset ,in Obtained through k-Wave simulation.

[0109] Inversion task

[0110] For each data entry:

[0111] initial value Reconstructed by FBP;

[0112] predict ;

[0113] residual ;

[0114] Target correction increment .

[0115] 3.2 Loss Function

[0116] Orthogonal loss

[0117]

[0118] Inversion loss

[0119]

[0120] Total losses

[0121]

[0122] The Adam optimizer was used, with an initial learning rate of 1×10⁻³, which decayed by 0.5 every 50 rounds.

[0123] The batch size is set to 12, and a recursive training strategy is used, that is, the network parameters are fixed in some iterations and only the input features are updated, thereby reducing memory usage;

[0124] The entire training cycle consists of approximately 500 rounds, with an early stop mechanism (automatic termination if there is no improvement after 20 consecutive rounds).

[0125] S4: Perform high-quality image reconstruction of real radiotherapy signals. Put the X-ray thermoacoustic signals collected in the radiotherapy environment into the trained dual-head FNO model to reconstruct high-quality X-ray thermoacoustic images, that is, obtain the reconstructed sound pressure distribution map.

[0126] Specifically, in a real radiotherapy environment, ultrasound signals are acquired using X-ray thermoacoustic imaging equipment, with the raw data being time-domain sampled data. The acquired data undergoes preprocessing: data normalization is performed; a low-resolution coarse reconstructed image is generated using traditional time reversal or backprojection algorithms, serving as input to the FNO module. The processed time-domain signal and this preliminary image are then input into a joint optimization network. The network first predicts the acoustic signal of the currently estimated image through a forward branch and calculates the residual with the actual measured signal; then, in the inversion branch, the residual is concatenated with the current image as input, and the learned "residual → image correction" operator directly outputs the image correction increment. The detail preservation, PSNR, and SSIM metrics of the network-reconstructed image are compared with those of traditional reconstruction methods (e.g., direct backprojection or CNN-based reconstruction), verifying the advantages of this invention under low sampling and low-dose conditions.

[0127] S5 includes at least the following steps: After completing the two-dimensional dose distribution reconstruction, a quantitative comparison is performed with the patient's own radiotherapy plan (RTPLAN). The specific steps are as follows:

[0128] Import the RT PLAN dose distribution and use the DICOM-RT standard analysis tool to read the planned dose matrix corresponding to the RT PLAN in the case. Ensure pixel spacing and Consistent.

[0129] Image registration uses rigid or affine registration algorithms for... and Perform spatial alignment and output the registered reconstructed image. .

[0130] Difference map generation

[0131]

[0132] Will Visualized as a pseudo-color difference map for qualitative observation.

[0133] Error Statistics

[0134] Calculate the mean absolute error (MAE) of the entire image.

[0135]

[0136] Maximum error

[0137]

[0138] Standard deviation

[0139]

[0140] The contours of ROIs such as PTV and OAR are read from the RT STRUCT of the case and converted into masks. For each ROI, in and Voxel dose values ​​were statistically analyzed to generate two DVH curves.

[0141] The above error statistics and DVH differences are summarized to generate a quantitative assessment report. The output includes: mean error (MAE), maximum error, standard deviation; difference plot and DVH comparison plot.

[0142] The above process allows for a direct and quantitative assessment of the consistency between the reconstructed dose of this invention and the clinical RT PLAN, providing a reliable basis for clinical dose monitoring and feedback. Quantitative evaluation indicators such as the gamma index are used to compare the planned dose with the reconstructed dose distribution. Commonly used evaluation criteria include 3% / 3 mm, 2% / 3 mm, 3% / 5 mm, and 2% / 5 mm. Gamma pass rate is calculated; in the example, the gamma pass rate can reach over 97%, thus validating the invention.

[0143] It will be apparent to those skilled in the art that the present invention is not limited to the details of the exemplary embodiments described above, and that the invention can be implemented in other specific forms without departing from its spirit or essential characteristics. Therefore, the embodiments should be considered in all respects as exemplary and non-limiting, and the scope of the invention is defined by the appended claims rather than the foregoing description. Thus, all variations falling within the meaning and scope of equivalents of the claims are intended to be included within the present invention. No reference numerals in the claims should be construed as limiting the scope of the claims.

Claims

1. A method for X-ray thermoacoustic imaging reconstruction based on Fourier neural operators, characterized in that: At least the following steps are included: S1: Construct an X-ray thermoacoustic data simulation and multimodal dataset, which includes an initial sound pressure image and the corresponding time-domain signal; S2: Construct a two-headed Fourier neural operator model, also known as a two-headed FNO model, and input the signal simulated by the two-headed FNO model into the joint optimization network; S3: Joint optimization network construction and training: The simulated ultrasound signal datasets from different sampling perspectives and the initial sound pressure images are fed into the joint optimization network for training; S4: Perform high-quality image reconstruction of real radiotherapy signals. Put the X-ray thermoacoustic signals collected in the radiotherapy environment into the trained dual-head FNO model to reconstruct high-quality X-ray thermoacoustic images, that is, obtain the reconstructed sound pressure distribution map. S5: Convert the reconstructed sound pressure distribution into a dose distribution according to the corresponding formula, and compare the dose distribution with the original radiotherapy plan.

2. The X-ray thermoacoustic imaging reconstruction method based on Fourier neural operators according to claim 1, characterized in that: The construction of the X-ray thermoacoustic data simulation and multimodal dataset includes at least the following steps: First, determine the tools and environment. In the MATLAB environment, use the K-Wave toolbox to simulate X-ray induced thermoacoustic signals under different sampling conditions. Secondly, signal simulation is performed to conduct detailed simulation of X-ray induced thermoacoustic signals under different sampling conditions. The different sampling conditions include at least the number of sensors, the shape of the sensor array, and the sampling angle. The data is then processed by inputting clinical CT images and dividing the tissue into different regions using HU value segmentation. These different regions include at least bone, soft tissue, and tumor. Then, preset thresholds are applied, and corresponding tissue densities and Grunison coefficients are assigned to each region to determine the dose-sound pressure conversion parameters. Subsequently, the initial sound pressure distribution is generated. According to the physical conversion formula, the radiation dose distribution in the radiotherapy plan is combined with the tissue characteristic parameters and converted into the initial sound pressure distribution. The initial sound pressure distribution is then stored as a grayscale image with a fixed resolution. Subsequently, sound wave propagation simulation was performed. Under the set ultrasonic sensor parameters, K-Wave was used to simulate sound wave propagation, and Gaussian white noise was added to generate ultrasonic simulation signal data under different signal-to-noise ratio conditions. Finally, the dataset was constructed. Based on the above steps, an X-ray thermoacoustic data simulation and multimodal dataset was built and divided into a training set and a test set in a 7:3 ratio.

3. The X-ray thermoacoustic imaging reconstruction method based on Fourier neural operators according to claim 2, characterized in that: The dual-head FNO model employs a dual-head multi-layer Fourier neural operator structure to extract global information from the input signal in the frequency domain, simulating the signal propagation effect caused by X-ray dose deposition in different tissues. The output of the dual-head FNO model can be used as both a simulated sensor time-domain signal and a preliminary estimate of the dose distribution, which can then be used as a basis for data alignment during subsequent iterative updates.

4. The X-ray thermoacoustic imaging reconstruction method based on Fourier neural operators according to claim 3, characterized in that: The dual-head FNO model is used to simultaneously realize acoustic forward and inverse updates on the same operator skeleton. The overall architecture of the dual-head FNO model adopts the process of "dimensionality increase - multi-layer Fourier layer - fusion - projection", but in the final projection, two independent projection branches are used, namely the forward branch and the inverse branch. Dimensional upscaling refers to mapping the input image to a higher-dimensional channel feature space through a dimensionality upscaling layer; Multi-layer Fourier layers are introduced, with each layer extracting low-frequency modes in the frequency domain and applying a learnable kernel integral operator; The inverse transform is then fused with the parallel linear mapping branch to capture global propagation characteristics.

5. The X-ray thermoacoustic imaging reconstruction method based on Fourier neural operators according to claim 4, characterized in that: The input to the forward modeling branch is the current initial sound pressure level image; The simulated time-domain signal tensor, which corresponds to the number of sensor array channels and time sampling points, is output through the projection layer to replace the traditional finite element / Monte Carlo forward modeling solution. The input to the inversion branch is the current image and its simulation residual, which are then mapped and stitched together to form an H×W×2 channel input. The function of the inversion branch is to directly learn the "residual → image correction" operator, replacing the traditional CNN.

6. The X-ray thermoacoustic imaging reconstruction method based on Fourier neural operators according to claim 4, characterized in that: The S3 includes at least the following steps: The ultrasonic signals and FBP-reconstructed X-ray thermoacoustic images in the ultrasonic simulation signal dataset are used as inputs. Furthermore, the labels of the forward branch are real time-domain signals generated by the K-Wave toolbox, while the labels of the inverse branch are the difference between the real pressure distribution of social security stocks and the preliminary reconstructed image. Task-weighted loss is used: L=αLfwd+(1−α)Linv Where Lfwd is the time-domain signal prediction error, Linv is the image update error, and α is the balance coefficient; The Adam optimization algorithm based on backpropagation iterates the network parameters; By combining the early stopping mechanism with validation set monitoring, the optimal parameters are selected.

7. The X-ray thermoacoustic imaging reconstruction method based on Fourier neural operators according to claim 5, characterized in that: The S4 includes at least the following steps: The joint optimization network in S3 was applied to the reconstruction of the actual acquired X-ray thermoacoustic signals; The ultrasound signals acquired in the actual radiotherapy environment were preprocessed, and a low-resolution initial image was obtained by simple back projection. The preprocessed data is input into a joint optimization network, and a dual-head FNO model is used to reconstruct high-resolution and high signal-to-noise ratio X-ray thermoacoustic images to obtain the reconstructed sound pressure distribution map.

8. The X-ray thermoacoustic imaging reconstruction method based on Fourier neural operators according to claim 6, characterized in that: The S5 includes at least the following steps: The reconstructed sound pressure distribution map is converted into the corresponding X-ray dose deposition distribution map according to the known physical conversion formula; Then, the RT PLAN dose distribution output from the clinical radiotherapy system is imported and rigidly / affinely registered with the reconstructed two-dimensional dose distribution on the same anatomical plane to ensure a one-to-one correspondence between pixels and voxels. The two-dimensional reconstructed sound pressure distribution image is converted into a dose distribution image using analytical formulas and tissue characteristic parameters; Calculate the DVH curves of reconstructed dose versus planned dose within the target area and organ contour, respectively; Based on statistical differences and DVH deviation, the consistency between the reconstructed dose and the clinical RT PLAN was evaluated, and the usability of this method in clinical dose monitoring and feedback was verified.

Images