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

Abstract

The invention discloses an X-ray thermoacoustic imaging reconstruction method based on a Fourier neural operator, and relates to the technical field of X-ray imaging. The method comprises the following steps: quickly constructing a multi-modal X-ray thermoacoustic data set covering different sampling conditions; the frequency domain global feature extraction advantage of the FNO is fully utilized, and the precision of forward physical model simulation and dose prediction is effectively improved; the double-head Fourier structure effectively improves the gradual correction effect of the absorption coefficient and the dose distribution, so that the reconstructed image is clearer, and the quantification is more accurate; lightweight training and recursive parameter updating reduce consumption of computing resources, so that the system has real-time or quasi-real-time imaging capability, a closed-loop system from initial data simulation and network reconstruction to dose distribution verification is completely realized, and the system has high clinical adaptability.

Description

Technical Field

[0001] The present invention relates to the technical field of X-ray imaging, and in particular to an X-ray thermoacoustic imaging reconstruction method based on a Fourier neural operator. Background Art

[0002] X-ray imaging technology has been a core tool in medical diagnostics for more than a century, but it has inherent limitations, such as the need to strictly align the X-ray detector with the radiation source and the potential carcinogenic risks of ionizing radiation.

[0003] X-ray-induced acoustic computed tomography (XACT), an emerging technology, has emerged in recent years. By utilizing short X-ray pulses to excite tissue and generate acoustic signals for imaging, XACT overcomes the limitations of conventional X-ray 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] The filtered back projection method (FBP) is widely used due to its high computational efficiency. However, since the acoustic signal excited by X-rays has 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 image, especially in low-dose scenarios where the artifact ratio can reach more than 30%. The time reversal method is based on the back propagation of the sound field signal through the wave equation. Although it can partially suppress noise, its performance is highly dependent on the completeness of the sensor layout. In actual imaging scenarios such as sparse sampling (such as a 32-channel annular array) or limited viewing angle (<180° coverage), the resolution of the reconstructed image drops by more than 40%, which cannot 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-based optimization) have become a research hotspot. For example, a framework combining alternating optimization of the absorption coefficient and dose distribution can reduce quantitative error to less than 5%. However, such methods require repeated solutions to the thermoacoustic wave equation to match the measured signal. A single forward simulation can take up to several minutes (based on the finite element method), resulting in overall reconstruction times exceeding 1 hour, making it difficult to meet clinical real-time requirements (such as intraoperative navigation requiring ≤5 seconds per frame). To overcome this bottleneck, researchers have attempted to introduce deep learning into XACT. Deep learning methods have been previously applied to XACT reconstruction. For example, end-to-end mapping based on U-Net can rapidly generate images, but requires a large amount of annotated data, is sensitive to noise (SSIM significantly decreases at low doses), and lacks generalization capability. Hybrid architecture methods (such as PixelInterp-CNN) partially alleviate the overfitting problem by combining time-domain signal encoding with initial back-projected images, but still require manual design of regularization terms and cannot adaptively optimize physical model parameters.

[0006] In recent years, neural operators (such as the Fourier neural operator, FNO) have demonstrated revolutionary advantages in solving partial differential equations (PDEs). By parameterizing the integral kernel in Fourier space, FNO can adapt to arbitrary spatial resolution and geometry, avoiding the repeated discretization calculations of traditional numerical methods (such as finite elements). In photoacoustic imaging, FNO has successfully replaced the finite element solver, increasing the speed of light transmission simulation by 30 times while maintaining an error of less than 3%. However, the difference in the physical model of XACT leads to the failure of direct migration. For example, the existing FNO architecture does not model the secondary electron scattering effect of X-rays in the bone area, resulting in a dose distribution estimation error of >15%. In addition, traditional methods process noise suppression, dose correction and absorption coefficient inversion in steps, ignoring their coupling relationship, resulting in cumulative error amplification. Although iterative reconstruction combined with FNO can partially alleviate the computational pressure, its high memory usage (>12GB) still limits its deployment in embedded medical devices. Therefore, a technical solution that deeply integrates physical models and deep learning is needed. By targeted design of network structure, joint optimization framework and lightweight strategy, efficient and high-precision imaging of XACT can be achieved to promote its practical application in clinical practice. Summary of the Invention

[0007] The purpose of the present invention is to provide an X-ray thermoacoustic imaging reconstruction method based on Fourier neural operators, aiming 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) FNO network structure with an iterative optimization framework, efficient and high-precision XACT image reconstruction and dose correction are achieved, while reducing computing resource requirements and promoting its application in clinical real-time imaging.

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

[0009] S1: constructing an X-ray thermoacoustic data simulation and a multimodal data set, wherein the X-ray thermoacoustic data simulation and the multimodal data set include an initial acoustic pressure image and a corresponding time domain signal;

[0010] S2: Construct a dual-head Fourier neural operator model, namely the dual-head FNO model, and input the signal simulated by the dual-head FNO model into the joint network;

[0011] S3: Joint network construction and training: The simulated ultrasonic signal datasets under different sampling angles and the initial sound pressure images are put into the joint optimization network for training;

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

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

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

[0015] First, we determined the tools and environment and used the K-Wave toolbox in the MATLAB environment to simulate X-ray-induced thermoacoustic signals under different sampling conditions.

[0016] Secondly, signal simulation is performed to finely simulate the X-ray induced thermoacoustic signal under different sampling conditions, wherein the different sampling conditions include at least the number of sensors, the shape of the sensor array, and the sampling angle;

[0017] Then, data processing is performed, clinical CT images are input, and tissues are divided into different regions by HU value segmentation, wherein the different regions include at least bones, soft tissues and tumors;

[0018] Then, a preset threshold is performed to assign the corresponding tissue density and Grünitz coefficient to each area to determine the dose-to-acoustic pressure conversion parameters;

[0019] Then, 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 to convert it into the initial sound pressure distribution, and the sound pressure distribution is stored as a grayscale image with a fixed resolution.

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

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

[0022] Furthermore, the dual-head FNO module uses 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;

[0023] The output of the dual-head FNO module can be both a simulated sensor time domain signal and a direct preliminary estimated dose distribution, so as to serve as a basis for data alignment during subsequent iterative updates.

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

[0025] Dimensionality increase means mapping the input image to a high-dimensional channel feature space through a dimension increase layer (fully connected or 1×1 convolution);

[0026] Introducing multiple layers of Fourier layers (e.g., four layers), each of which intercepts 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 the global propagation characteristics;

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

[0029] The projection layer (1×1 convolution) outputs the simulated time-domain signal tensor corresponding to the sensor array channels and time sampling points, replacing the traditional finite element / Monte Carlo forward solution;

[0030] The input of the inversion branch is the current image and its simulation residual, which are mapped and spliced ​​into H×W×2 channel input;

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

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

[0033] The ultrasonic signal in the ultrasonic simulation signal dataset and the X-ray thermoacoustic image reconstructed by FBP are used as input;

[0034] And the label of the forward branch is the real time domain signal generated by the K-Wave toolbox (or other numerical simulation), while the label of the inversion branch is used to calculate the difference between the real pressure distribution and the preliminary reconstructed image (or the real absorption coefficient minus the initial value);

[0035] Using task-weighted loss:

[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 RMSE or SSIM combination), and α is the balance coefficient;

[0038] Iterate network parameters based on the back-propagation Adam (or AdamW) optimization algorithm;

[0039] Combine the early stopping mechanism with the validation set monitoring to select the optimal parameters.

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

[0041] Apply the joint network in S3 to the reconstruction of actual collected X-ray thermoacoustic signals;

[0042] Preprocess the real ultrasound signals collected in the radiotherapy environment and use simple back-projection (such as time reversal algorithm) to obtain low-resolution initial images;

[0043] The preprocessed data is input into the joint network, and the FNO module is used to achieve high-resolution and high signal-to-noise ratio X-ray thermoacoustic image reconstruction to obtain the reconstructed sound pressure distribution map;

[0044] The reconstructed images are superior to traditional methods in detail recovery, noise suppression, and quantitative indicators (such as absorption coefficient and dose accuracy).

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

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

[0047] Then, the RT PLAN dose distribution (DICOM format) output by the clinical radiotherapy system is imported and rigid / affine registration is performed with the reconstructed 2D dose image on the same anatomical plane to ensure one-to-one correspondence between pixels / voxels.

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

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

[0050] Based on the above difference statistics and DVH deviation, the consistency between the reconstructed dose and the clinical RT PLAN was evaluated, and the feasibility of this method in clinical dose monitoring and feedback was verified.

[0051] Compared with the prior art, the present invention has the following beneficial effects:

[0052] The present invention rapidly constructs a multimodal X-ray thermoacoustic dataset covering different sampling conditions; fully utilizes the frequency domain global feature extraction advantage of FNO to effectively improve the accuracy of forward physical model simulation and dose prediction; the dual-stage iterative optimization structure effectively enhances the gradual correction effect of the absorption coefficient and dose distribution, making the reconstructed image clearer and more quantitatively accurate; lightweight training and recursive parameter updates reduce computing resource consumption, enabling the system to have real-time or quasi-real-time imaging capabilities, fully realizing a closed-loop system from initial data simulation, network reconstruction to dose distribution verification, and having high clinical adaptability. BRIEF DESCRIPTION OF THE DRAWINGS

[0053] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the following briefly introduces the drawings required for describing the embodiments. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without creative work.

[0054] Figure 1 This is a schematic diagram of the initial pressure in the data set of the present invention;

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

[0056] Figure 3 This is a flow chart based on FNO of the present invention;

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

[0058] Figure 5 Schematic diagram of human body X-ray thermoacoustic reconstruction based on 128 sensors according to the present invention. DETAILED DESCRIPTION

[0059] The technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, rather than all the embodiments.

[0060] A method for reconstructing X-ray thermoacoustic imaging based on a Fourier neural operator comprises at least the following steps:

[0061] S1: constructing an X-ray thermoacoustic data simulation and a multimodal data set, wherein the X-ray thermoacoustic data simulation and the multimodal data set include an initial acoustic pressure image and a corresponding time domain signal;

[0062] The process of obtaining the initial sound pressure image in S1 is as follows Figure 1As shown in the figure, tissue characteristic parameter extraction and initial sound pressure generation are performed by collecting clinical CT images in DICOM format, with a resolution of 512×512 pixels and a layer thickness of 1 mm. The image areas are divided into three categories according to the HU value using the grayscale value of the CT image: bone area: HU>300, soft tissue area: 50≤HU≤300, and tumor area: HU<50. Corresponding tissue characteristic parameters are assigned to each area: bone: tissue density ρ=1.85g / cm 3 , Grunisson coefficient Γ=0.15, soft tissue: ρ=1.05g / cm 3 , Γ=0.25, tumor: ρ=1.10 g / cm 3 , Γ=0.30.

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

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

[0065] Convert the radiation dose into the initial sound pressure distribution.

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

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

[0068] S2: Construct a dual-head Fourier neural operator model, namely the dual-head FNO model, and input the signal simulated by the dual-head FNO model into the joint network;

[0069] Input X-ray thermoacoustic data simulation and initial acoustic pressure image in multimodal dataset as input signals;

[0070] A multi-layer Fourier neural operator structure is used 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.

[0071] The output of the FNO module can be both a simulated sensor time domain signal and a direct preliminary estimated dose distribution, so as to serve as a basis for data alignment during subsequent iterative updates.

[0072] Specifically, the dual-head FNO joint network construction 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, and the input is first promoted from low-dimensional space to high-dimensional representation through a fully connected layer to facilitate subsequent global frequency domain processing. A four-layer Fourier neural operator is constructed, and each layer is designed as follows: a two-dimensional Fourier transform F is performed on the input to transfer the spatial domain features to the frequency domain; in the frequency domain, only low-frequency modes are retained (for example, 64 low-frequency modes are retained, covering 0.5–5MHz), and a linear transformation R is applied to the low-frequency modes to enhance the global effective information; an inverse Fourier transform F is performed on the processed frequency domain data -1 , restoring it to the spatial domain; at the same time, a parallel branch is set up to directly apply a 1×1 convolution layer (Wt0) to the input to perform a linear transformation; the output of the upper branch is further enhanced with low-frequency features through two 1×1 convolution layers (Wt1 and Wt2) before being fused with the bottom branch to obtain the final output of this layer. After processing through four layers of Fourier layers, the output v4(x) is mapped back to the target dimension through two 1×1 convolution layers to generate an estimated X-ray dose deposition distribution or a simulated output of the sensor's time domain signal, providing the initial mapping result 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 the input of the forward branch, single channel.

[0077] Inversion input: In the kth round of iterative reconstruction:

[0078] Use the trained FNO_forward to estimate the current image x k (x) Perform a forward modeling to obtain the simulated signal

[0079] and the measured signal y meas (c,t) calculate the residual

[0080]

[0081] to r k(c, t) is interpolated or convolved to map it to the same spatial resolution as the image, and the

[0082] With the current estimated image x k (x) Splice in the channel dimension to form a dual-channel input tensor

[0083]

[0084] 2.2 Network Backbone A unified backbone network is used to process the above two inputs: 1. Dimensionality Upgrading Layer

[0085] Apply 1×1 convolution to the input tensor V0 (single channel in forward modeling, dual channel in inversion modeling) to increase the number of channels to C=64, and output

[0086]

[0087] 2. Four layers of Fourier layers, for each layer t=0,1,2,3, v t (x) Send in:

[0088] Frequency domain branch

[0089] Two-dimensional Fourier transform:

[0090] Only the first M=32 low-frequency patterns are retained, and the high frequencies are filtered out.

[0091] Apply the kernel integral operator (Learnable Kernel) R:V′(k)=R·V(k).

[0092] Inverse Fourier transform:

[0093] Two cascaded 1×1 convolutional layers W t1 ,W t2 and GELU activation, resulting in

[0094] The direct branch is connected to the original input v t (x) Directly apply 1×1 convolution W t0 , we get v′ t (x).

[0095] Fusion

[0096]

[0097] 2.3 Output Head

[0098] Forward branch

[0099] Apply 1×1 convolution to v4(x) to map the channels to C sens=128 (sensor channels) × T = 500 (time steps) dimension, output

[0100]

[0101] Inversion branch

[0102] Apply another set of 1×1 convolutions to the same v4(x) to map the channels into a single channel rectified delta

[0103]

[0104] S3: Joint network construction and training: The simulated ultrasonic signal datasets under different sampling angles and the initial sound pressure images are put into the joint optimization network for training;

[0105] S3 specifically includes the following steps:

[0106] 3.1 Training Data Preparation

[0107] Forward mission

[0108] Dataset in Obtained by k-Wave simulation.

[0109] Inversion task

[0110] For each data:

[0111] Initial value Reconstructed by FBP;

[0112] predict

[0113] residual

[0114] Target Correction Increment

[0115] 3.2 Loss Function

[0116] Forward loss

[0117]

[0118] Inversion loss

[0119]

[0120] Combined total loss

[0121] L total =αL fwd +(1-α)L inv ,α=0.5

[0122] The Adam optimizer is used, the initial learning rate is set to 1×10-3, and it decays 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 to reduce memory usage;

[0124] The entire training cycle is about 500 rounds, with an early stopping mechanism (automatically terminated if there is no improvement for 20 consecutive rounds).

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

[0126] Specifically, in an actual radiotherapy environment, an X-ray thermoacoustic imaging device is used to collect ultrasound signals, and the original data is time-domain sampling data. The collected data is preprocessed: data normalization is performed; a low-resolution coarse reconstructed image is generated using a traditional time reversal or back-projection algorithm as the input of the FNO module. The processed time-domain signal is input into the joint network together with the preliminary image. The network first predicts the acoustic signal of the current estimated image through the forward branch and calculates the residual with the actual measured signal; then, in the inversion branch, the residual is spliced ​​with the current image as input, and the image correction increment is directly output through the learned "residual→image correction" operator. The detail retention, PSNR and SSIM indicators of the network-reconstructed image are compared with those of traditional reconstruction methods (such as direct back-projection or CNN-based reconstruction) to verify the advantages of the present invention under low sampling and low-dose conditions.

[0127] S5 includes at least the following steps: After completing the two-dimensional dose image reconstruction, quantitative comparison is performed with the radiotherapy plan (RTPLAN) provided with the case. The specific steps are as follows:

[0128] Import RT PLAN dose distribution and use DICOM-RT standard parsing tools to read the planned dose matrix D of the corresponding level of RTPLAN in the case. plan (x,y), ensuring pixel spacing and D rec consistent.

[0129] Image registration uses rigid or affine registration algorithms to rec With D plan Perform spatial alignment and output the reconstructed image after registration

[0130] Difference map generation

[0131]

[0132] Visualize ΔD 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] max|ΔD(x,y)|

[0138] Standard deviation

[0139]

[0140] Read the ROI contours of PTV, OAR, etc. from the case RT STRUCT and convert them into masks. With D plan The voxel dose values ​​are counted separately to generate two DVH curves.

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

[0142] The above process allows for intuitive and quantitative evaluation of the consistency between the reconstructed dose according to the present invention and the clinical RT plan, providing a reliable basis for clinical dose monitoring and feedback. Quantitative evaluation metrics such as the gamma index are used to compare the planned dose with the reconstructed dose distribution. Common evaluation criteria include 3% / 3mm, 2% / 3mm, 3% / 5mm, and 2% / 5mm. The gamma pass rate is calculated and, in this example, exceeds 97%, validating the present 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 embodied in other specific forms without departing from the spirit or essential characteristics of the invention. Therefore, the embodiments should be considered in all respects as illustrative and non-restrictive, and the scope of the invention is defined by the appended claims, not the foregoing description, and all variations within the meaning and range of equivalents of the claims are intended to be included therein. Any reference sign in a claim should not be construed as limiting the claim to which it relates.

Claims

1. A Fourier neural operator-based X-ray thermoacoustic imaging reconstruction method, characterized by: At least the following steps are included: S1: constructing an X-ray thermoacoustic data simulation and a multimodal data set, wherein the X-ray thermoacoustic data simulation and the multimodal data set include an initial acoustic pressure image and a corresponding time domain signal; S2: Construct a dual-head Fourier neural operator model, namely the dual-head FNO model, and input the signal simulated by the dual-head FNO model into the joint network; S3: Joint network construction and training: The simulated ultrasonic signal datasets under different sampling angles and the initial sound pressure images are put into the joint optimization network for training; S4: Perform high-quality image reconstruction of real radiotherapy signals. X-ray thermoacoustic signals collected in the radiotherapy environment are put into the trained dual-head FNO network to reconstruct high-quality X-ray thermoacoustic images, that is, to obtain a reconstructed sound pressure distribution map. S5: Convert the reconstructed acoustic pressure distribution into a dose distribution according to a corresponding formula, and compare the dose image with the original radiotherapy plan.

2. The X-ray thermoacoustic imaging reconstruction method based on Fourier neural operator according to claim 1, characterized in that: The construction of X-ray thermoacoustic data simulation and multimodal data set comprises at least the following steps: First, we determined the tools and environment and used the K-Wave toolbox in the MATLAB environment to simulate X-ray-induced thermoacoustic signals under different sampling conditions. Secondly, signal simulation is performed to finely simulate the X-ray induced thermoacoustic signal under different sampling conditions, wherein the different sampling conditions include at least the number of sensors, the shape of the sensor array, and the sampling angle; Then, data processing is performed, clinical CT images are input, and tissues are divided into different regions by HU value segmentation, wherein the different regions include at least bones, soft tissues and tumors; Then, a preset threshold is performed to assign the corresponding tissue density and Grünitz coefficient to each area to determine the dose-to-acoustic pressure conversion parameters; Then, 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 to convert it into the initial sound pressure distribution, and the sound pressure distribution is stored as a grayscale image with a fixed resolution. Then, the acoustic wave propagation simulation was carried out. Under the set ultrasonic sensor parameters, K-Wave was used to simulate the acoustic wave propagation. Gaussian white noise was added to generate ultrasonic simulation signal data under different signal-to-noise ratio conditions. Finally, the dataset is constructed. Based on the above steps, an X-ray thermoacoustic data simulation and multimodal dataset are constructed and divided into a training set and a test set in a ratio of 7:

3.

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

4. The X-ray thermoacoustic imaging reconstruction method based on Fourier neural operator according to claim 3, characterized in that: The dual-head FNO model is used to simultaneously implement acoustic wave forward modeling and inversion updates on the same operator skeleton. The overall architecture of the dual-head FNO model adopts the "dimensionality increase - multi-layer Fourier layer - fusion - projection" process, but in the final projection, two independent projection branches are used, namely the forward branch and the inversion branch. Dimensionality upgrading means mapping the input image to a high-dimensional channel feature space through a dimensionality upgrading layer; Introducing multiple layers of Fourier layers, each of which intercepts low-frequency modes in the frequency domain and applies a learnable kernel integral operator; The inverse transform is then fused with the parallel linear mapping branch to capture the global propagation characteristics.

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

6. The X-ray thermoacoustic imaging reconstruction method based on Fourier neural operator according to claim 4, characterized in that: The S3 at least includes the following steps: The ultrasonic signal in the ultrasonic simulation signal dataset and the X-ray thermoacoustic image reconstructed by FBP are used as input; And the label of the forward branch is the real time domain signal generated by the K-Wave toolbox, while the label of the inversion branch is the difference between the real pressure distribution of the social security stock and the preliminary reconstructed image; Using task-weighted loss: L=αLfwd+(1-α)Linv Among them, Lfwd is the time domain signal prediction error, Linv is the image update error, and α is the balance coefficient; Iterate network parameters based on the back-propagation Adam optimization algorithm; Combine the early stopping mechanism with the validation set monitoring to select the optimal parameters.

7. The X-ray thermoacoustic imaging reconstruction method based on Fourier neural operator according to claim 5, characterized in that: The S4 at least includes the following steps: Apply the joint network in S3 to the reconstruction of actual collected X-ray thermoacoustic signals; The ultrasound signals collected in the radiotherapy environment are preprocessed and a low-resolution initial image is obtained using simple back-projection. The preprocessed data are input into the joint network, and the FNO module is used to achieve high-resolution and high signal-to-noise ratio X-ray thermoacoustic image reconstruction to obtain the reconstructed sound pressure distribution map.

8. The X-ray thermoacoustic imaging reconstruction method based on Fourier neural operator according to claim 6, characterized in that: The S5 at least includes the following steps: The reconstructed acoustic pressure distribution map is converted into the corresponding X-ray dose deposition distribution according to a known physical conversion formula; Then, the RT PLAN dose distribution output by the clinical radiotherapy system is imported and reconstructed into a two-dimensional dose image, and rigid / affine registration is performed on the same anatomical plane to ensure one-to-one correspondence between pixels / voxels. The two-dimensional reconstructed acoustic pressure distribution image is converted into a dose distribution image through analytical formulas and tissue characteristic parameters; Calculate the DVH curves of the reconstructed dose and the planned dose within the target volume and major organ contours respectively; Based on the above difference statistics and DVH deviation, the consistency between the reconstructed dose and the clinical RT PLAN was evaluated, and the feasibility of this method in clinical dose monitoring and feedback was verified.

Images