Tumor breathing movement trajectory prediction method and system based on spatiotemporal feature separation LSTM (Long Short Term Memory)

Abstract

The invention discloses a tumor respiratory movement trajectory prediction method and system based on spatio-temporal feature separation LSTM, and the method comprises the following steps: collecting a four-dimensional CT respiratory cycle image, extracting a tumor centroid coordinate, and generating a standardized tumor centroid coordinate sequence through sliding window standardization; the method comprises the following steps: inputting bidirectional LSTM to calculate Euclidean distance screening features and weighting to generate a weighted spatial feature vector sequence, extracting time recursion feature weighted fusion through one-dimensional convolution and unidirectional LSTM, predicting a tumor motion trajectory through autoregression, and comparing real-time coordinate correction to generate an adaptive correction trajectory prediction result. A four-dimensional image sequence is extracted and standardized through registration, a two-way network is combined to calculate two-way distance features to recognize breathing phase changes, spatial-temporal features are fused to establish a nonlinear recursive relation to capture tumor motion differences, real-time correction optimization is achieved through deviation statistics and weight self-adjustment, and prediction stability and trajectory fitting precision are improved.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of artificial intelligence, in particular to a tumor respiratory motion trajectory prediction method and system based on spatiotemporal feature separation LSTM. BACKGROUND

[0002] The field of artificial intelligence involves the research and application of using computer systems to simulate human intelligent behavior, covering machine learning, deep learning, neural networks, natural language processing, computer vision, pattern recognition, and data intelligent analysis. Through learning and feature extraction on a large amount of data, the computer has the ability of perception, understanding, reasoning and decision-making, which is widely used in medical image analysis, speech recognition, autonomous driving and intelligent diagnosis. In medical applications, artificial intelligence technology is particularly important in tumor positioning, radiotherapy planning and motion prediction, which comprehensively analyzes the temporal and spatial features to achieve intelligent modeling and prediction of complex physiological motion patterns.

[0003] Among them, the traditional tumor respiratory motion trajectory prediction method refers to the technology of modeling and predicting the tumor motion trajectory in the process of radiotherapy to reduce the tumor position deviation caused by patient respiratory motion. The traditional method is based on respiratory signal or body surface marker displacement data, using linear regression model, polynomial fitting or traditional recurrent neural network to learn time series, so as to obtain the displacement trajectory of tumor with respiratory changes. This method collects four-dimensional CT data or body surface motion data of patients, extracts time series features, and establishes a unified mapping relationship model to realize the prediction of tumor motion trajectory. However, due to the limitations of traditional models in spatiotemporal feature separation and dynamic coupling processing, the prediction accuracy is insufficient under complex respiratory motion.

[0004] The existing technology relies on one-way time series modeling and static feature extraction in tumor respiratory motion trajectory prediction, and the prediction process lacks collaborative analysis of spatial features and temporal dynamics, resulting in distortion of tumor centroid displacement features. Respiratory motion has individual differences and nonlinear fluctuations, and traditional models are insufficient in responding to abnormal respiratory rhythm or non-uniform motion, causing trajectory prediction lag and accumulated deviation. The model relies on fixed mapping relationship in the training stage, lacks adaptive correction based on real-time data feedback, and is difficult to cope with the dynamic changes of respiratory state in the process of radiotherapy, affecting the stability of trajectory prediction and the accuracy of radiotherapy positioning. SUMMARY

[0005] In order to solve the technical problems that the prior art relies on one-way time series modeling and static feature extraction in tumor respiratory motion trajectory prediction, the prediction process lacks collaborative analysis of spatial features and time dynamics, resulting in distortion of tumor centroid displacement features. Respiratory motion has individual differences and nonlinear fluctuations. Traditional models are insufficient in responding to abnormal respiratory rhythm or non-uniform motion, causing trajectory prediction lag and accumulated deviation. The model relies on fixed mapping relationship in the training stage, lacks adaptive correction based on real-time data feedback, and is difficult to cope with the dynamic changes of respiratory state in the process of radiotherapy, affecting the stability of trajectory prediction and the accuracy of radiotherapy positioning. The embodiment of the present application provides a tumor respiratory motion trajectory prediction method based on spatiotemporal feature separation LSTM.

[0006] In order to achieve the above-mentioned purpose, the tumor respiratory motion trajectory prediction method based on spatiotemporal feature separation LSTM is adopted, which comprises the following steps: S1: Collecting chest and abdominal image sequence data of a patient during a respiratory cycle by a four-dimensional computed tomography device, extracting tumor centroid coordinate values based on B-spline non-rigid registration transformation to align phase image frames, and generating standardized tumor centroid coordinate sequence by sliding window standardization processing; S2: Calling the standardized tumor centroid coordinate sequence, inputting the bidirectional long short-term memory network to extract time series features and calculate Euclidean distance feature values, selecting coordinate points exceeding the preset distance threshold to obtain respiratory phase encoding values, and generating weighted spatial feature vector sequence after dimension consistency and coordinate component weighting; S3: Calling the weighted spatial feature vector sequence to perform one-dimensional convolution operation to calculate the feature difference value of adjacent time, and weightedly obtain the time local change feature, input the unidirectional long short-term memory network to extract the recursive feature and weightedly superimpose to generate the time recursive enhancement feature vector sequence; S4: Calling the time recursive enhancement feature vector sequence to perform linear transformation mapping to unify the dimension, then weightedly summing according to the fusion weight proportion, intercepting the current time window feature sequence to input the full connection layer to perform autoregressive prediction, and outputting the tumor respiratory motion trajectory prediction result.

[0007] As a further scheme of the present application, the standardized tumor centroid coordinate sequence includes spatial position standardized value, respiratory cycle phase index and centroid displacement standard deviation, the weighted spatial feature vector sequence includes spatial coordinate weighting coefficient, phase normalization value and bidirectional Euclidean distance feature value, the time recursive enhancement feature vector sequence includes time difference value, local change amplitude coefficient and recursive feature weight, and the tumor respiratory motion trajectory prediction result includes predicted trajectory coordinate sequence, time step prediction error and respiratory motion change trend parameter. The fusion weight proportion is the weight sequence generated by the attention mechanism layer in the time series prediction network.

[0008] As a further scheme of the present application, the specific steps of S1 are: S101: Obtain chest and abdomen image data in a patient's respiratory cycle through a four-dimensional computed tomography device, sort the image frames according to timestamps, calculate the gray value gradient and judge the phase stability, remove the mutation frames, retain the continuous phase frame interval and aggregate the stable frames, and generate a respiratory phase image frame sequence data; S102: Based on the respiratory phase image frame sequence data, perform B-spline non-rigid registration transformation on the adjacent frame voxel gray distribution, calculate the offset vector and adjust it, count the residual error to judge the registration convergence, and obtain a registration aligned image frame data set; S103: According to the registration aligned image frame data set, detect the tumor region pixel barycenter and extract the three-dimensional coordinates, construct a coordinate set, and use a sliding window to normalize the mean and standard deviation of the coordinates, and obtain a standardized tumor centroid coordinate sequence.

[0009] As a further scheme of the present application, the specific steps of S2 are: S201: Obtain the standardized tumor centroid coordinate sequence, input the bidirectional long short-term memory network to extract the time sequence feature, calculate the positive and negative Euclidean distances of adjacent coordinate points in three-dimensional axis components, perform average processing on the bidirectional distance matrix, and obtain a bidirectional Euclidean distance feature matrix; S202: According to the distance component value in the bidirectional Euclidean distance feature matrix, perform difference judgment with a preset distance threshold value, filter out coordinate points with positive difference, extract the time index sequence to calculate the phase displacement parameter, and perform numerical coding according to the phase interval mapping, to generate a respiratory phase coding set; S203: Call the respiratory phase coding set, according to the index matching standardized tumor centroid coordinate components, normalize the coding value and linearly weight the coordinate components, and obtain a weighted spatial feature vector sequence.

[0010] As a further scheme of the present application, the preset distance threshold value is a quantitative parameter determined according to the overall dynamic amplitude feature of the original tumor centroid coordinate sequence in the main motion direction, and in combination with the statistical quantity of the Euclidean distance of adjacent frames of the standardized sequence; The generation of the respiratory phase coding set refers to extracting adjacent time intervals of the filtered positive difference coordinate points according to the order of time indexes, calculating the phase displacement parameter by using the ratio of the time interval to the average time interval of the overall sequence, and forming the respiratory phase coding set by numerically coding the phase displacement parameter.

[0011] As a further scheme of the present application, the specific steps of S3 are: S301: Obtain the weighted spatial feature vector sequence, pair them according to the time sequence order of the feature vectors at adjacent time points, perform one-dimensional convolution calculation on elements of the same dimension and extract the difference value, aggregate the convolution output according to the time step sequence and update the feature index to generate a temporal difference feature sequence. S302: Call the time-series differential feature sequence, extract time slice feature components, construct time weight parameters based on time step difference, perform weighting operation, adjust the position weight ratio based on time step index, rearrange the weighting result and integrate the feature structure to generate a time local change feature sequence; S303: Call the time local change feature sequence, input it into the unidirectional long short-term memory unit according to the time step sequence, calculate the transition difference between hidden states and perform weighted fusion with the output vector to generate a time recursive enhanced feature vector sequence.

[0012] As a further aspect of the present invention, the specific steps of S4 are as follows: S401: Based on the time recursive enhanced feature vector sequence, perform linear mapping operation on the time step feature components of the sequence, perform dimension calibration according to the mapping matrix parameters, perform linear combination of the time step matrix components according to the mapping coefficients and check the consistency to generate a unified dimension feature set; S402: Call the time step feature vector in the unified dimension feature set, perform weighted summation according to the fusion weight ratio, filter continuous frames in combination with the fusion coefficient threshold, extract window frames according to the time index of the weighted feature set, and generate a time window feature sequence. The fusion coefficient threshold is determined by the statistical regression analysis results of the fusion effect of time series feature vectors in the training sample set; S403: Based on the continuous frame data of the time window feature sequence, input the fully connected layer, perform element-wise difference calculation on the numerical vector, establish an autoregressive prediction relationship based on the time step index, rearrange the output layer displacement vector sequence, and generate the tumor respiratory motion trajectory prediction result.

[0013] As a further aspect of the present invention, the method further includes step S5: S5: Call the tumor respiratory motion trajectory prediction results, compare them with the coordinates collected in real time by the radiotherapy equipment sensor, calculate the Euclidean distance deviation value, count the proportion of deviations exceeding the preset deviation threshold, adjust the fusion weight coefficient, re-execute the prediction, and generate adaptive correction trajectory prediction results; The adaptive correction trajectory prediction results include the fusion weight adjustment coefficient, the prediction error correction amount, and the trajectory deviation distribution characteristics.

[0014] As a further aspect of the present invention, the specific steps of S5 are as follows: S501: Call the predicted coordinate data in the tumor respiratory motion trajectory prediction result, match it with the target tissue three-dimensional motion coordinates collected in real time by the radiotherapy equipment sensor, extract the corresponding coordinate component values ​​and perform difference analysis to generate a coordinate deviation distance sequence set; S502: Call the coordinate deviation distance sequence set, compare the deviation distance value of the time node in the sequence with the deviation threshold, identify and mark the time node whose deviation exceeds the threshold, calculate the proportion of nodes exceeding the limit, organize the deviation change characteristics according to the time distribution, and generate the accuracy deviation proportion distribution. The deviation threshold is dynamically adjusted based on the distribution characteristics of the coordinate deviation distance data after time series pairing; S503: Based on the accuracy deviation ratio distribution, adjust the multi-source data fusion weight coefficients in the trajectory prediction results, use the adjusted weights to correct the prediction parameter set, calculate the trajectory point coordinates, and generate adaptive correction trajectory prediction results.

[0015] A tumor respiratory motion trajectory prediction system based on spatiotemporal feature separation LSTM includes: The image time sequence acquisition module acquires chest and abdominal image sequence data during the patient's respiratory cycle using a four-dimensional computed tomography (CT) scanner. After aligning the phase image frames based on B-spline non-rigid registration transformation, the tumor centroid coordinate values ​​are extracted. The tumor centroid coordinate sequence is generated through sliding window standardization and then transmitted to the centroid feature extraction module. The centroid feature extraction module calls the standardized tumor centroid coordinate sequence, inputs it into a bidirectional long short-term memory network to extract temporal features and calculates Euclidean distance feature values, filters coordinate points that exceed a preset distance threshold to obtain respiratory phase encoding values, performs dimensionality consistency and weights the coordinate components to generate a weighted spatial feature vector sequence, and passes it to the spatial feature construction module. The spatial feature construction module calls the weighted spatial feature vector sequence to perform a one-dimensional convolution operation to calculate the feature difference between adjacent time moments, and obtains the local temporal change features by weighting. It then inputs the recursive features extracted by the unidirectional long short-term memory network and performs weighted superposition to generate a temporal recursive enhanced feature vector sequence, which is then passed to the temporal feature recursion module. The time feature recursive module calls the time recursive enhanced feature vector sequence to perform a linear transformation mapping to a unified dimension, then performs a weighted summation according to the fusion weight ratio, extracts the feature sequence of the current time window and inputs it into the fully connected layer to perform autoregressive prediction, outputs the tumor respiratory motion trajectory prediction result, and passes it to the adaptive correction module. The adaptive correction module calls the tumor respiratory motion trajectory prediction results, compares them with the coordinates collected in real time by the radiotherapy equipment sensors, calculates the Euclidean distance deviation value, counts the proportion of deviations exceeding the preset deviation threshold, adjusts the fusion weight coefficient, re-executes the prediction, and generates the adaptive correction trajectory prediction result.

[0016] Compared with the prior art, the beneficial effects of the present invention are as follows: In this invention, by registering and standardizing four-dimensional image sequences, and combining this with a bidirectional network structure to calculate bidirectional distance features to identify respiratory phase changes, the dynamic response capability of spatial features is improved. A nonlinear recursive relationship is established using weighted fusion of spatial and temporal feature sequences, enabling the model to capture subtle differences in tumor movement under respiratory fluctuations and achieve multi-dimensional spatiotemporal coupling prediction. Through continuous feedback bias statistics and weight self-adjustment mechanisms, real-time correction and optimization of trajectory prediction are achieved, significantly improving prediction stability and trajectory fitting accuracy under complex respiratory conditions, and ensuring positioning consistency and control reliability during radiotherapy. Attached Figure Description

[0017] Figure 1 This is a schematic diagram of the steps of the present invention; Figure 2 This is a detailed schematic diagram of S1 of the present invention; Figure 3 This is a detailed schematic diagram of S2 of the present invention; Figure 4 This is a detailed schematic diagram of S3 of the present invention; Figure 5 This is a detailed schematic diagram of S4 of the present invention; Figure 6 This is a detailed schematic diagram of S5 of the present invention; Figure 7 This is a system module diagram of the present invention. Detailed Implementation

[0018] Specific embodiments of the present invention are given below. These specific embodiments are only used to further illustrate the present invention in detail and do not limit the scope of protection of the present invention.

[0019] In embodiments of the present invention, words such as "exemplarily," "for example," etc., are used to indicate that something is an example, illustration, or description. Any embodiment or design described as "exemplary" in the present invention should not be construed as being more preferred or advantageous than other embodiments or designs. Specifically, the use of the word "exemplary" is intended to present the concept in a concrete manner. Furthermore, in embodiments of the present invention, the meaning expressed by "and / or" can be both, or either one.

[0020] In the embodiments of this invention, the terms "image" and "picture" may sometimes be used interchangeably. It should be noted that, without emphasizing the distinction between them, they convey the same meaning. Similarly, the terms "of," "corresponding (relevant)," and "corresponding" may sometimes be used interchangeably. It should be noted that, without emphasizing the distinction between them, they convey the same meaning.

[0021] In this embodiment of the invention, sometimes the subscript such as W1 is written in a non-subscript form such as W1. When the difference is not emphasized, the meaning they express is the same.

[0022] Please see Figure 1 This invention provides a method for predicting tumor respiratory motion trajectory based on spatiotemporal feature separation LSTM, including the following steps: S1: Acquire chest and abdominal image sequence data during the patient's respiratory cycle using a four-dimensional computed tomography (CT) scanner. After aligning the phase image frames based on B-spline non-rigid registration transformation, extract the tumor centroid coordinate values. Then, generate a standardized tumor centroid coordinate sequence through sliding window standardization. S2: Call the standardized tumor centroid coordinate sequence, input the bidirectional long short-term memory network to extract temporal features and calculate the Euclidean distance feature value, filter the coordinate points that exceed the preset distance threshold to obtain the respiratory phase encoding value, perform dimensionality consistency and weight the coordinate components to generate a weighted spatial feature vector sequence; S3: Call the weighted spatial feature vector sequence to perform one-dimensional convolution operation to calculate the feature difference between adjacent time steps, and obtain the local temporal change features by weighting. Input the recursive features into the unidirectional long short-term memory network, perform weighted superposition, and generate a temporal recursive enhanced feature vector sequence. S4: Call the time recursive enhanced feature vector sequence to perform linear transformation mapping to unify the dimension, then sum the weights according to the fusion weight ratio, extract the feature sequence of the current time window and input it into the fully connected layer to perform autoregressive prediction, and output the tumor respiratory motion trajectory prediction result; S5: Call the tumor respiratory motion trajectory prediction results, compare them with the real-time coordinates collected by the radiotherapy equipment sensors, calculate the Euclidean distance deviation value, count the proportion of deviations exceeding the preset deviation threshold, adjust the fusion weight coefficient, re-execute the prediction, and generate adaptive correction trajectory prediction results. The standardized tumor centroid coordinate sequence includes standardized values ​​of spatial location, respiratory cycle phase index, and centroid displacement standard deviation. The weighted spatial feature vector sequence includes spatial coordinate weighting coefficients, phase normalization values, and bidirectional Euclidean distance feature values. The time-recursive enhanced feature vector sequence includes time difference scores, local change amplitude coefficients, and recursive feature weights. The tumor respiratory motion trajectory prediction results include predicted trajectory coordinate sequence, time step prediction error, and respiratory motion change trend parameters. The adaptive correction trajectory prediction results include fusion weight adjustment coefficients, prediction error correction, and trajectory deviation distribution characteristics.

[0023] Please see Figure 2 The specific steps of S1 are as follows: S101: Acquire chest and abdominal image data of the patient during the respiratory cycle using a four-dimensional computed tomography (CT) scanner. Sort the image frames by timestamp, calculate the gray value gradient and determine the phase stability, remove abrupt frames, retain continuous phase frame intervals and aggregate stable frames to generate respiratory phase image frame sequence data. Based on chest and abdominal images of patients during respiratory cycles acquired using a four-dimensional computed tomography (CT) scanner, the specific execution process first involves extracting the acquisition timestamp of each image frame from the DICOM file header. For example, the timestamp is obtained from the tag (0008, 0032). This process is then applied to an image sequence containing 50 frames (5 respiratory cycles, 10 frames per cycle). The time sequence is precisely ordered, and then adjacent image frames f are sorted. i and f i+1 The grayscale gradient is calculated by comparing the differences in grayscale values ​​between two frames, voxel by voxel. For example, for two 512×512×80 voxel images, the grayscale values ​​of corresponding voxel positions are taken. and Calculate the sum of the absolute differences in gray levels of voxels, i.e. If the SAD calculation value between frames f3 and f4 is 280000, and the SAD calculation value between frames f4 and f5 is 890000, then the phase stability is determined. This determination requires setting a phase stability gradient threshold T. grad The threshold is set with reference to the average SAD value μ of adjacent frames throughout the sequence. SAD and standard deviation σ SAD Assuming the mean of the 49 SAD values ​​calculated from 50 frames of image data is 300,000 and the standard deviation is 150,000, then the following can be set: Inter-frame changes with grayscale gradients below 210,000 are considered stable. Gradients between 210,000 and 600,000 are considered normal respiratory motion changes, while gradients above 600,000 are considered phase abrupt changes, such as those caused by coughing or body swaying. Therefore, the SAD value between frames f3 and f4 (280,000) is less than 600,000, indicating a stable transition. However, the SAD value between frames f4 and f5 (890,000) is greater than 600,000, so frame f5 is identified as a phase abrupt change and removed. After removing the phase abrupt change frame, the remaining image frames form multiple consecutive phase frame intervals. For example, frames in the original sequence... At this point, stable image frames in the same respiratory phase across different respiratory cycles are aggregated. For example, the first frame representing the end of expiration (0% phase) in the first respiratory cycle, the 11th frame in the second cycle, and the 21st frame in the third cycle are aggregated into the same phase group, "end-expiration group." Finally, the stable frames are arranged according to 10 respiratory phases ( The data is then reorganized to generate respiratory phase image frame sequence data.

[0024] S102: Based on respiratory phase image frame sequence data, B-spline non-rigid registration transformation is performed on the voxel grayscale distribution of adjacent frames, the offset vector is calculated and adjusted, and the statistical residual is used to determine the registration convergence, thus obtaining the registered and aligned image frame dataset. Select two adjacent frames within a phase group, such as the end-expiratory phase frame (0% phase), as a fixed image I. fixed And 10% of the frames are used as floating images I moving In fixed image I fixed A control point grid is defined with a grid spacing of 20×20×15 voxels. This is achieved by moving the floating image I. moving The control point positions are determined, and the displacement of voxels between grids is calculated using B-spline basis function interpolation, thereby deforming the floating image. Then, a non-rigid B-spline registration transformation is performed on the voxel grayscale distribution of adjacent frames. During registration, the offset vector is calculated and adjusted based on the similarity measure between the two images. Here, Normalized Mutual Information (NMI) is used as the metric. The mutual information value is obtained by calculating the joint probability distribution p(i,j) of grayscale values ​​in the overlapping area of ​​the two images, as well as the edge probability distributions p(i) and p(j). For example, the calculated NMI value is 1.215 in the k-th iteration. After adjusting the control point positions, the calculated NMI value in the (k+1)-th iteration is 1.218. This increment of 0.003 indicates that the adjustment enhances the statistical correlation of grayscale values ​​between the two images. Next, the residual is statistically analyzed to determine whether the registration has converged. Here, the residual is defined as the change in the NMI value between two consecutive iterations. And set a convergence threshold T. convThis threshold is set with reference to the target registration accuracy. For example, for a registration accuracy requirement of 1mm, T can be set as follows: conv Setting it to 0.0005, in the example above, This value is greater than the threshold T conv Therefore, the iteration continues. When the change in NMI value for three consecutive iterations is less than 0.0005, for example, the NMI is 1.241 for the nth iteration, 1.2414 for the (n+1)th iteration, 1.2417 for the (n+2)th iteration, and 1.2419 for the (n+3)th iteration, with changes of 0.0004, 0.0003, and 0.0002 respectively, all less than 0.0005, the registration process is considered to have converged, the iteration stops, and the transformed version of the floating image at this time is saved. This process is repeated for image frames of adjacent time phases to obtain the registered and aligned image frame dataset.

[0025] S103: Based on the registered and aligned image frame dataset, the centroid of the tumor region pixels is detected and the three-dimensional coordinates are extracted. A coordinate set is constructed, and the coordinates are normalized by mean and standard deviation using a sliding window to obtain a standardized tumor centroid coordinate sequence. Based on the generated registration and alignment image frame dataset, an initial 3D region of interest (ROI) containing the tumor is first manually drawn or determined using semi-automatic tools on the image frame at a reference time phase (e.g., 0% end-tidal). Then, the centroid of the tumor pixels within this ROI is detected and its 3D coordinates are extracted. Specifically, a grayscale threshold range is first set to segment the tumor. Referring to the characteristic of lung cancer solid tumors in CT images where the Henlein unit (HU) value is in the range of [-20, 60], voxels within this ROI with HU values ​​in this range are selected. and its corresponding grayscale value I(V) i The tumor voxel was identified, and then the gray-weighted average coordinates of the voxel were calculated as the centroid. For example, if there are three tumor voxels within a ROI, with coordinates and gray values ​​of (101, 255, 83) where the gray value is 30, (102, 255, 83) where the gray value is 45, and (101, 256, 83) where the gray value is 35, then... Axial centroid coordinates C is calculated in this way. y and C z This process is applied to each frame of the image sequence to construct a coordinate set consisting of N three-dimensional coordinate points. Then, a sliding window of size W=5 is used to normalize the coordinates in the coordinate set by mean and standard deviation. When the window slides to the i-th position of the sequence (covering 5 coordinate points from i to i+4), the mean and standard deviation of these 5 coordinate points on the X, Y, and Z axes are calculated. For example, the X-axis coordinate within the window is... (Unit: mm), with a mean value of μ x It is 22.16 mm, with a standard deviation σ. x The value is 0.48mm. Then, the normalized coordinate points are calculated (here, the last point in the window, i.e., the (i+4)th point). Coordinates are Perform the same operation on the Y-axis and Z-axis coordinates, and slide the window one position backward. Repeat this process until the entire coordinate sequence has been traversed. Specific data can be found in Table 1.

[0026] Table 1: Examples of Tumor Centroid Coordinates and Normalized Representations

[0027] As shown in Table 1, the table lists the original three-dimensional coordinates of the tumor at some timestamps and their coordinate values ​​after sliding window normalization. The first four timestamps could not be calculated because the window was not full. Normalization results were generated starting from the fifth timestamp, and finally a standardized tumor centroid coordinate sequence that can be used for subsequent model training was obtained.

[0028] Please see Figure 3 The specific steps of S2 are as follows: S201: Obtain the standardized tumor centroid coordinate sequence, input it into a bidirectional long short-term memory network to extract time series features, calculate the forward and reverse Euclidean distances between adjacent coordinate points in the three-dimensional axis components, perform averaging on the bidirectional distance matrix, and obtain the bidirectional Euclidean distance feature matrix. Obtain a standardized sequence of tumor centroid coordinates, for example, a sequence containing 100 time points, from t1 to t2. 100 sequence ,in This sequence is simultaneously fed into a forward processing unit and a backward processing unit, and the forward unit is processed according to... Processed sequentially, backward units are processed according to... The sequential processing extracts time-series features that integrate past and future information at each time point. Then, it calculates the coordinates of any two points in the sequence. and The Euclidean distance between the three-dimensional axial components is calculated by extracting the multi-axis coordinate components of the two points. and This is achieved by calculating the square root of the sum of the squares of their differences, for example, calculating the coordinates of time point t5. Coordinates of time point t6 The positive Euclidean distance d(5,6) between them is calculated as follows: Simultaneously, the reverse Euclidean distance d(6,5) is calculated, which is equal to d(5,6), both being 0.5408. Subsequently, the forward distance d(i,j) and the reverse distance d(j,i) calculated for point pair (i,j) are averaged. Specifically, d(i,j) and d(j,i) are added together and then divided by 2. Since d(i,j) = d(j,i), this average value is the Euclidean distance value calculated in one direction itself. The distance values ​​d(i,j) are combined into a 100×100 matrix to obtain the bidirectional Euclidean distance feature matrix.

[0029] S202: Based on the distance component values ​​in the bidirectional Euclidean distance feature matrix, perform difference judgment with the preset distance threshold, filter out coordinate points with positive difference, extract the time index sequence to calculate the phase displacement parameters, and perform numerical encoding based on the phase interval mapping to generate a respiratory phase code set; Based on the neighbor distance component values ​​d(i,i+1) in the bidirectional Euclidean distance feature matrix, and a preset distance threshold T, dist The difference judgment is performed, and the preset distance threshold is determined based on the dynamic amplitude statistics of the original, unstandardized tumor centroid coordinate sequence. Specifically, the determination process involves first calculating the maximum and minimum displacements of the entire coordinate sequence on the axis with the largest motion amplitude (the head-to-feet direction, i.e., the Z-axis). For example, if the Z-axis coordinate range is [-25.2, -15.8] mm, then the maximum dynamic amplitude is... mm, then calculate the average Euclidean distance μ between adjacent frames in the normalized space. d Suppose that the value of d(i,i+1) is extracted from the bidirectional Euclidean distance feature matrix, and its mean is calculated to be μ. d =0.45, standard deviation is σ d =0.18, then the threshold T dist Set as This setting aims to filter out inflection points where the movement speed is above average. Subsequently, it filters out coordinate points with positive differences, i.e., for each adjacent point distance d(i,i+1), it calculates d(i,i+1)-T. dist If the result is greater than 0, then time point i+1 is marked as a key point. For example, the result calculated above... The difference between it and the threshold is Therefore, time point t6 is selected. Assume the time index sequence of consecutive key points selected using this method is... Next, the time interval between adjacent key points is extracted, i.e. Then, the average time interval between the selected key points is calculated. For example, if the average interval is 5.5 time units, the phase displacement parameter is calculated using the ratio of the time interval to this average time interval. The phase displacement parameter for the first interval is... The second one is The third one is Finally, numerical encoding is performed through phase interval mapping. The mapping rule here is as follows: phase displacement parameter values ​​less than 0.85 are encoded as 1 (indicating acceleration), values ​​between 0.85 and 1.15 are encoded as 2 (indicating constant speed), and values ​​greater than 1.15 are encoded as 3 (indicating deceleration). According to this rule, the three phase displacement parameters 0.91, 0.73, and 1.27 are encoded as 2, 1, and 3 respectively, generating a respiratory phase encoding set. .

[0030] S203: Call the respiratory phase encoding set, match the standardized tumor centroid coordinate components according to the index, normalize the encoded values ​​and linearly weight them with the coordinate components to obtain the weighted spatial feature vector sequence; Calling the respiratory phase encoding set Based on the time index corresponding to each code, for example, code 2 corresponds to time index 11, code 1 corresponds to time index 15, and code 3 corresponds to time index 22, the corresponding coordinate components in the obtained standardized tumor centroid coordinate sequence are matched, as shown in Table 2. The extracted code values ​​are then normalized using the min-max normalization method, with the code value range set to [1,3]. The normalization calculation is (code value - 1) / (3 - 1). Therefore, the normalized value of code 2 is... The normalized value of code 1 is The normalized value of code 3 is Then, the normalized encoded values ​​are linearly weighted with the corresponding coordinate components. This weighting process requires preset weight coefficients; here, the coordinate component weights are set. and phase encoding weights The sum of these two weight values ​​is 1. This is based on the premise that the spatial location of the tumor is the primary information, while phase encoding serves as secondary, supplementary temporal information. Taking the data at time index 11 as an example, its standardized coordinates are... The corresponding phase encoding is 2, which becomes 0.5 after normalization. Therefore, the weighted spatial feature vector... The calculation is as follows: The X component is Similarly, the Y component is The Z component is Therefore, the weighted vector is Repeat this process for the selected key points to obtain a weighted spatial feature vector sequence.

[0031] Table 2: Example Table of Phase Encoding and Coordinate Weighting

[0032] As shown in Table 2, this table displays the original phase code matched based on the time index, its normalized value, the corresponding standardized coordinate components, and the weighted spatial feature components obtained by the final linear weighted calculation, demonstrating the fusion process of spatial information and temporal phase information.

[0033] Please see Figure 4 The specific steps of S3 are as follows: S301: Obtain the weighted spatial feature vector sequence, pair them according to the time sequence order of the feature vectors at adjacent time points, perform one-dimensional convolution calculation on elements of the same dimension and extract the difference value, aggregate the convolution output according to the time step sequence and update the feature index to generate the temporal difference feature sequence. Obtain a generated sequence containing 20 three-dimensional vectors. First, pairing is performed based on the time sequence order of feature vectors at adjacent times. Specifically, the t-th vector in the sequence is... With the (t+1)th vector Form a pair For example, a vector with time index 11 With the vector of time index 12 After pairing, a one-dimensional convolution is performed on elements of the same dimension, and the difference value is extracted. This convolution is equivalent to applying a one-dimensional convolution kernel of size 2 with a kernel of [-1, 1] and sliding it across the paired vectors. For the X dimension, the calculation involves performing a dot product operation between the kernel [-1, 1] and the X dimension component [-0.725, -0.750]. This result is the difference between the two moments in the X dimension. For the Y-dimensional components Perform the same operation to obtain Performing the same operation on the Z-dimensional components [-0.564, -0.590] yields the following result. By performing this calculation on the dimension, the first difference vector is obtained. Then advance the time step forward to... Perform the same convolution difference calculation to obtain This process is applied to adjacent pairings in the entire weighted space feature vector sequence, and the calculated difference vector is obtained according to the original time step sequence. Aggregate the sequences and assign new feature indices from 1 to 19 to generate temporal difference feature sequences.

[0034] S302: Call the time-series differential feature sequence, extract time-slice feature components, construct time weight parameters based on the time step difference, perform weighted calculations, and adjust the position weight ratio according to the time step index, using the following formula:

[0035] Calculate the location weight ratio parameter, rearrange the weighted results, and integrate the feature structure to generate a temporally local variation feature sequence; where β t ΔX represents the position weight ratio parameter at time step t. i Represents the temporal difference feature component at time i. W represents the absolute value of the temporal difference characteristic components. i γ represents the time weight parameter corresponding to time i, λ represents the normalization adjustment constant, and γ represents the time weight parameter corresponding to time i. t The baseline weight offset represents time step t, n represents the number of feature components in the time slice, and i represents the time slice index. Call The positional weighting parameter β at time step t is calculated. t First, a feature component and its corresponding time weight parameter are extracted within a time slice of length n ending at time t. Here, the time slice length n=3. Then, when calculating β5, the feature component and its corresponding time weight parameter are extracted. Three vectors are used, and the dimension component with the most significant change in motion is selected as ΔX. i Assuming the Z-axis is the primary direction of motion, then the Z-axis component is extracted, i.e. And extract the corresponding time weight parameter W. i The weighting parameter is set based on the time step difference; features closer to the current time step have higher weights. A linearly decaying weight is used here. For a time slice of length 3, the weights are set as follows: Their sum is 1. Then, after performing a weighted operation, the position weight ratio is adjusted according to the time step index, and the position weight ratio parameter β is calculated using the formula. t This formula is achieved by dividing the sum of the absolute values ​​of the weighted time-series difference features by a denominator smoothed by the L2 norm and adjustment constant λ, and then adding a baseline weight offset γ. t To dynamically adjust the weights at each time step, where β t ΔX is the position weight ratio parameter for the final calculated time step t. i It refers to the temporal difference feature component at the i-th original moment relative to the current time step t within the time slice. It is the absolute value of that component, W i It is the time weight assigned to the i-th original time step, where λ is a small positive number, let's call it λ. γ is used to maintain numerical stability when the denominator is close to zero. t This is the baseline weight offset, which is set with reference to the average fluctuation of the entire sequence difference values. If the overall fluctuation of the sequence is small, a lower baseline value is set. Here, based on the statistics of the original data, γ is set. t=0.1, n is the number of feature components in the time slice, here n=3, and i is the index in the time slice, from 1 to n.

[0036] The advantage of the formula lies in the fact that it uses time weights W i With difference feature ΔX i The absolute values ​​of the features are multiplied and normalized with the L2 norm of the features to dynamically calculate the weight ratio β at each time step. t This allows the weights to consider not only the drasticness of the change, but also its temporal relationship.

[0037] Table 3: Example Table for Calculating Location Weight Ratio Parameters

[0038] As shown in Table 3, it is assumed that the extracted Z-axis difference components and their corresponding weights are as described above when calculating β5.

[0039] Based on the data in Table 3, β5 is calculated: First, the numerator term is calculated: ; Then calculate the sum of squares in the denominator: ; Next, calculate the complete denominator: ; Finally, calculate β5: ; The result shows that the position weighting parameter for time step 5 is 0.5919. Applying this calculation process to each time step yields a series of β values. t Value, using each β t The value is multiplied by the original time-series difference vector at the corresponding time point. ,For example and the newly generated vector Arranged chronologically, the feature structures are integrated to generate a sequence of features showing local temporal changes.

[0040] S303: Call the time local change feature sequence, input it into the unidirectional long short-term memory unit according to the time step sequence, calculate the transition difference between hidden states and combine it with the output vector for weighted fusion to generate a time recursive enhanced feature vector sequence; Call a by The resulting sequence of temporal local variation features is input into the unidirectional long short-term memory unit step by step according to the time step sequence. Specifically, during execution, for time step t, the feature vector... Compared with the hidden state vector h of the previous time step t-1 and cell state vector C t-1Together, they are sent to the computation unit. Inside the unit, the forget gate first determines which information to discard from the cell state. This calculation will... and h t-1 After concatenation, the result is multiplied by the forgetting gate weight matrix, a bias term is added, and then passed through a sigmoid activation function to obtain a forgetting factor between 0 and 1. This factor is related to the cell state C at the previous time step. t-1 Element-wise multiplication is performed, followed by an input gate to determine which new information is stored in the cell state. This process consists of two parts: one part uses a sigmoid activation function to determine which values ​​will be updated, and the other part uses a tanh activation function to create a new candidate value vector. The result is then multiplied by the forgotten C... t-1 Adding them together yields a new cell state C. t Then, the output gate is calculated to determine the hidden state h of the output. t ,Will and h t-1 The new cell state C is then activated using a sigmoid activation function. t The new hidden state h is obtained by multiplying the two states using a tanh activation function. t and output vector o t Then, the transition differences between hidden states are calculated, i.e. For example, if the previous hidden state The currently calculated hidden state Then the transfer difference Then, this difference vector is compared with the output vector o at the current time. t Perform weighted fusion, assuming the current output vector And set the fusion weight coefficient, which is determined based on experimental tuning, and the differential weight w diff Set to 0.4, output weight w out If we set it to 0.6, then the weighted fused vector E t Calculated as E calculated at time step t By combining them in sequence, a time-recursive enhanced feature vector sequence is generated.

[0041] Please see Figure 5 The specific steps of S4 are as follows: S401: Based on time recursion, enhance the feature vector sequence, perform linear mapping operation on the feature components of the sequence time step, perform dimension calibration according to the mapping matrix parameters, perform linear combination of the time step matrix components according to the mapping coefficients and check the consistency, and generate a unified dimension feature set. Based on the generated sequence containing 19 three-dimensional vectors Where vector E1 is First, a linear mapping operation is performed on the feature components at each time step in the sequence. This operation is implemented using a pre-defined mapping matrix M. The parameters of this matrix are learned during the model training phase through the backpropagation algorithm. Its function is to transform the input three-dimensional feature vector into a higher-dimensional five-dimensional feature vector. In this example, the mapping matrix M is a 5×3 matrix with the following values: ; Dimensional calibration is performed based on the parameters of this mapping matrix, specifically by adjusting the input 3×1 vector E. t Multiplying it by a 5×3 matrix M yields a 5×1 output vector F. t For example, for vectors The mapping is performed, and the calculation of its first component is as follows: By analogy, five components are calculated to obtain... Next, the time step matrix components are linearly combined according to the mapping coefficients and consistency is checked. This process applies the above linear mapping operation to each vector E in the sequence. t Using the same mapping matrix M for calculation, the resulting five-dimensional vector These are combined to generate a unified dimensional feature set.

[0042] S402: Call the time step feature vector in the unified dimension feature set, perform weighted summation according to the fusion weight ratio, filter continuous frames by combining the fusion coefficient threshold, extract window frames by time index of the weighted feature set, and generate time window feature sequence. The fusion coefficient threshold is determined by the statistical regression analysis results of the fusion effect of time series feature vectors in the training sample set; Call the time-step feature vectors from the unified-dimensional feature set, for example First, based on a preset fusion weight ratio sequence A weighted summation is performed; this weight sequence, calculated by the attention mechanism layer in the network, reflects the importance of features at each time step. The weight sequence shown here is... Next, a fusion coefficient threshold T is combined. fusion To filter out consecutive frames with higher weights, the fusion coefficient threshold T is used. fusion The fusion weights were determined through statistical regression analysis on the relationship between the fusion weights of time-series feature vectors and the final prediction error in a training sample set containing 500 different patients. The analysis results showed that when the fusion weights a t When T is greater than 0.12, the corresponding eigenvectors significantly contribute to the prediction results, therefore T is... fusion The value is set to 0.12, and the filtering process involves comparing the weight 'a' one by one.t With T fusion ,For example If so, then discard feature F1. Then feature F2 is retained. If feature F3 is retained, as shown in Table 4, a fixed-size window frame is extracted from the retained weighted feature set after filtering, based on its original time index. Here, the window size is set to 5. From the filtered feature sequence... In the middle, the first window is captured as The second window then slides backward to capture... The captured windows are combined to generate a time window feature sequence.

[0043] Table 4: Example Table of Feature Filtering and Window Construction

[0044] As shown in Table 4, the table lists the fusion weights for some time steps and their selection decisions based on a threshold of 0.12. The decision directly determines which feature vectors will be included in subsequent time windows.

[0045] S403: Based on the continuous frame data of the time window feature sequence, input the fully connected layer, perform element-wise difference calculation on the numerical vector, establish an autoregressive prediction relationship based on the time step index, rearrange the output layer displacement vector sequence, and generate the tumor respiratory motion trajectory prediction result. Based on the time window feature sequence of consecutive frame data, such as the first window , where each F t Each vector is a five-dimensional vector. First, the vector within this window is flattened into a single input vector with dimensions 5×5=25. This flattened vector is then input to the input node of a fully connected layer. This fully connected layer contains 25 input nodes and 3 output nodes. Each output node corresponds to the predicted displacement along the X, Y, and Z axes at a future time step. Inside the fully connected layer, the 25-dimensional input vector is multiplied by a 25×3 weight matrix and a 3×1 bias vector is added. Then, an activation function is applied to obtain the output. Next, element-wise interpolation is performed on the output numerical vector. This process is reflected in the network's output as the displacement vector. Instead of absolute coordinates, an autoregressive prediction relationship is established based on the time step index, that is, using the raw information within the time window. To predict the displacement at the next time point immediately following, for example, the normalized displacement vector output by the fully connected layer is... Finally, the output layer displacement vector sequence is rearranged. This rearrangement refers to converting the predicted normalized displacement vectors into displacements in the real-world coordinate system through an inverse normalization operation. This operation requires calling the mean and standard deviation of the corresponding time window, assuming the last known point P... last If the standard deviation of the Z-axis in the window is 1.5 mm, then the predicted true displacement of the Z-axis is 0.55 × 1.5 = 0.825 mm. The calculated true displacement of the three-dimensional axis is added to the true coordinates of the last known point to obtain the predicted coordinates of the next time point. This process is repeated to connect the prediction results of each time window to generate the prediction result of the tumor respiratory motion trajectory.

[0046] Please see Figure 6 The specific steps of S5 are as follows: S501: Call the predicted coordinate data in the tumor respiratory motion trajectory prediction results, match it with the target tissue three-dimensional motion coordinates collected in real time by the radiotherapy equipment sensor, extract the corresponding coordinate component values ​​and perform difference analysis to generate a coordinate deviation distance sequence set; Access the predicted coordinate data from the tumor respiratory motion trajectory prediction results, such as a predicted coordinate sequence containing 100 time points. Simultaneously, the time series of the target tissue's three-dimensional motion coordinates, acquired in real-time by the radiotherapy equipment's sensors, is paired with the coordinates acquired in real-time. This pairing operation will predict the coordinates of the same time index t. with actual coordinates For example, at time t=50, the predicted coordinates are... mm, while the true coordinates obtained through real-time monitoring using electromagnetic beacons implanted near the tumor are Next, the corresponding coordinate component values ​​of the two paired coordinates are extracted and a difference analysis is performed. This difference analysis is achieved by calculating the three-dimensional Euclidean distance between the two points. For the coordinate pair at time point t=50, the difference between the components is first calculated, i.e. mm, mm, Then, calculate the square root of the sum of squares of the differences to obtain the coordinate deviation distance at that time point. mm, applying this calculation process to paired coordinates at 100 time points yields a sequence containing 100 offset distance values, for example ,in mm, mm, ..., mm, ..., mm, this sequence is the coordinate deviation distance sequence set.

[0047] S502: Call the coordinate deviation distance sequence set, compare the deviation distance value of the time node in the sequence with the deviation threshold, identify and mark the time node whose deviation exceeds the threshold, calculate the proportion of nodes exceeding the limit, organize the deviation change characteristics according to the time distribution, and generate the accuracy deviation proportion distribution. The deviation threshold is dynamically adjusted based on the distribution characteristics of the coordinate deviation distance data after time series pairing; Call It contains 100 deviation distance values. First, the deviation distance value at each time point in the sequence is compared with a dynamically adjusted deviation threshold T. prec The dynamic adjustment of this deviation threshold is based on the calculation of the mean μ of the deviation distance data within the most recent respiratory cycle (assuming 30 time points). dev and standard deviation σ dev And set the threshold to For example, for the current time point, the mean of the first 30 deviation distance values ​​is μ. dev =1.2mm, standard deviation is σ dev =0.5mm, then the current deviation threshold mm, then, identify and mark the time points where the deviation exceeds the threshold. Specifically, this involves setting the distance value D for each deviation in the sequence. t Compared to 2.2mm, if D t If the value is greater than 2.2 mm, then time node t is marked as an out-of-limit node. For example, when D... 55 When D = 2.31mm, node 55 is marked because 2.31 > 2.2, while when D 54 When the value is 2.15mm, node 54 is not marked because 2.15 < 2.2. After judging all 100 nodes, the proportion of nodes exceeding the limit is calculated. Assuming that a total of 8 nodes exceeding the limit are marked, the proportion of nodes exceeding the limit is 8 / 100 = 8%. Then, the deviation change characteristics are sorted out according to the time distribution, as shown in Table 5. The respiratory cycle is divided into three phase intervals: peak inspiratory phase (phase 0-0.25), expiratory phase (phase 0.25-0.75), and trough expiratory phase (phase 0.75-1.0). By analyzing the time position of the 8 nodes exceeding the limit in the original respiratory signal, it was found that 6 nodes fell in the expiratory phase, accounting for 75% of the total nodes exceeding the limit. This information is integrated to generate the accuracy deviation proportion distribution.

[0048] Table 5: Phase Distribution of Deviation Nodes

[0049] As shown in Table 5, this table statistically analyzes the distribution of time points when deviations exceed the limit in different respiratory phase intervals, revealing that the deviations are mainly concentrated in the expiratory phase.

[0050] S503: Based on the distribution of accuracy deviation ratio, adjust the weight coefficients of multi-source data fusion in the trajectory prediction results, use the adjusted weights to correct the prediction parameter set, calculate the trajectory point coordinates, and generate adaptive correction trajectory prediction results. Based on the distribution of accuracy deviation proportions, namely, the proportion of out-of-limit nodes is 8% and 75% of the deviations are concentrated in the expiratory phase, the weighting coefficients of multi-source data fusion in the trajectory prediction results are adjusted. The adjustment target here is the coordinate component weight w defined above. coord and phase encoding weight w phase Since the deviation mainly occurs in specific respiratory phases (expiration), this indicates that the temporal dynamic information introduced by phase encoding fails to accurately reflect the actual motion during this phase. Therefore, it is necessary to reduce its influence. The adjustment rule is set as follows: when the total exceedance ratio is below 10% but above 5%, and the concentration of deviation in a specific phase region exceeds 70%, the encoding weight w of that phase is adjusted. phase Reduce its current value by 20%, and correspondingly increase the coordinate component weight w. coord To keep the sum at 1, the original weights are w. coord =0.7, w phase =0.3, the adjustment process is to calculate the change in weight. The adjusted new weights are , The prediction parameter set is corrected using the adjusted weights, and the new weights are used. and Recalculate the weighted spatial feature vector sequence; for example, for the vector at time index 11, the newly calculated weighted X-coordinate is... After recalculating the feature vector with new weights, the resulting corrected feature sequence is used as input to perform temporal difference feature extraction, generation of temporal local change feature sequence, generation of temporal recursive enhanced feature vector sequence, generation of unified dimension feature set, generation of time window feature sequence, and final autoregressive prediction calculation. Through this entire set of subsequent processes, the trajectory point coordinates are recalculated to generate adaptive correction trajectory prediction results.

[0051] Please see Figure 7 A tumor respiratory motion trajectory prediction system based on spatiotemporal feature separation LSTM includes: The image time sequence acquisition module acquires chest and abdominal image sequence data during the patient's respiratory cycle using a four-dimensional computed tomography (CT) scanner. After aligning the phase image frames based on B-spline non-rigid registration transformation, the tumor centroid coordinate values ​​are extracted. The tumor centroid coordinate sequence is generated through sliding window standardization and then transmitted to the centroid feature extraction module. The centroid feature extraction module calls a standardized tumor centroid coordinate sequence, inputs a bidirectional long short-term memory network to calculate bidirectional Euclidean distance feature values, filters coordinate points that exceed a preset distance threshold to obtain corresponding respiratory phase encoding values, performs normalization transformation and weights the coordinate components to generate a weighted spatial feature vector sequence, and passes it to the spatial feature construction module. The spatial feature construction module calls the weighted spatial feature vector sequence to perform one-dimensional convolution operation to calculate the feature difference between adjacent time steps, and obtains the local temporal change features by weighting. It inputs the recursive features extracted by the unidirectional long short-term memory network, performs weighted superposition, generates a temporal recursive enhanced feature vector sequence, and passes it to the temporal feature recursion module. The time feature recursion module calls the time recursion enhanced feature vector sequence to perform a linear transformation mapping to a unified dimension, and then performs a weighted summation according to the fusion weight ratio. It extracts the feature sequence of the current time window and inputs it into the fully connected layer to perform autoregressive prediction, outputting the tumor respiratory motion trajectory prediction result, which is then passed to the adaptive correction module. The adaptive correction module calls the tumor respiratory motion trajectory prediction results, compares them with the real-time coordinates collected by the radiotherapy equipment sensors, calculates the Euclidean distance deviation value, counts the proportion of deviations exceeding the preset deviation threshold, adjusts the fusion weight coefficients, re-executes the prediction, and generates the adaptive correction trajectory prediction result.

[0052] Any aspects not covered in this invention are applicable to existing technologies.

Claims

1. A tumor respiratory motion trajectory prediction method based on spatiotemporal feature separation LSTM, characterized in that, Includes the following steps: S1: Acquire chest and abdominal image sequence data during the patient's respiratory cycle using a four-dimensional computed tomography (CT) scanner. After aligning the phase image frames based on B-spline non-rigid registration transformation, extract the tumor centroid coordinate values. Then, generate a standardized tumor centroid coordinate sequence through sliding window standardization. S2: Call the standardized tumor centroid coordinate sequence, input the bidirectional long short-term memory network to extract temporal features and calculate Euclidean distance feature values, filter coordinate points that exceed the preset distance threshold to obtain respiratory phase encoding values, perform dimensionality consistency and weight the coordinate components to generate a weighted spatial feature vector sequence; S3: Call the weighted spatial feature vector sequence to perform one-dimensional convolution operation to calculate the feature difference between adjacent time steps, and obtain the local temporal change features by weighting. Input the recursive features into the unidirectional long short-term memory network to extract recursive features and perform weighted superposition to generate a temporal recursive enhanced feature vector sequence. S4: Call the time recursive enhanced feature vector sequence to perform a linear transformation mapping to a unified dimension, then sum the weights according to the fusion weight ratio, extract the feature sequence of the current time window, input it into the fully connected layer to perform autoregressive prediction, and output the tumor respiratory motion trajectory prediction result.

2. The tumor respiratory motion trajectory prediction method based on spatiotemporal feature separation LSTM according to claim 1, characterized in that, The standardized tumor centroid coordinate sequence includes spatial location standardized values, respiratory cycle phase index, and centroid displacement standard deviation; the weighted spatial feature vector sequence includes spatial coordinate weighting coefficients, phase normalization values, and bidirectional Euclidean distance feature values; the time recursive enhanced feature vector sequence includes time difference scores, local change amplitude coefficients, and recursive feature weights; and the tumor respiratory motion trajectory prediction result includes predicted trajectory coordinate sequence, time step prediction error, and respiratory motion change trend parameters. The fusion weight ratio is a weight sequence generated by the attention mechanism layer in the temporal prediction network.

3. The tumor respiratory motion trajectory prediction method based on spatiotemporal feature separation LSTM according to claim 1, characterized in that, The specific steps of S1 are as follows: S101: Acquire chest and abdominal image data of the patient during the respiratory cycle using a four-dimensional computed tomography (CT) scanner. Sort the image frames by timestamp, calculate the gray value gradient and determine the phase stability, remove abrupt frames, retain continuous phase frame intervals and aggregate stable frames to generate respiratory phase image frame sequence data. S102: Based on the respiratory phase image frame sequence data, perform B-spline non-rigid registration transformation on the voxel grayscale distribution of adjacent frames, calculate and adjust the offset vector, statistically determine the registration convergence by the residual, and obtain the registered and aligned image frame dataset. S103: Based on the registered and aligned image frame dataset, detect the centroid of the tumor region pixels and extract the three-dimensional coordinates, construct a coordinate set, and use a sliding window to normalize the coordinates by mean and standard deviation to obtain a standardized tumor centroid coordinate sequence.

4. The tumor respiratory motion trajectory prediction method based on spatiotemporal feature separation LSTM according to claim 3, characterized in that, The specific steps of S2 are as follows: S201: Obtain the standardized tumor centroid coordinate sequence, input it into a bidirectional long short-term memory network to extract time series features, calculate the positive and negative Euclidean distances between adjacent coordinate points in the three-dimensional axis components, perform averaging on the bidirectional distance matrix, and obtain the bidirectional Euclidean distance feature matrix. S202: Based on the distance component values ​​in the bidirectional Euclidean distance feature matrix, perform difference judgment with a preset distance threshold, filter out coordinate points with positive differences, extract time index sequences to calculate phase displacement parameters, and perform numerical encoding based on phase interval mapping to generate a respiratory phase encoding set; S203: Call the respiratory phase encoding set, match the standardized tumor centroid coordinate components according to the index, normalize the encoded values ​​and linearly weight them with the coordinate components to obtain a weighted spatial feature vector sequence.

5. The tumor respiratory motion trajectory prediction method based on spatiotemporal feature separation LSTM according to claim 4, characterized in that, The preset distance threshold is a quantitative parameter determined based on the overall dynamic amplitude characteristics of the original tumor centroid coordinate sequence in the main motion direction, combined with the statistical value of the Euclidean distance between adjacent frames of the standardized sequence. The generation of the respiratory phase coding set refers to extracting adjacent time intervals from the selected positive difference coordinate points according to the order of time index, calculating the phase displacement parameter by the ratio of the time interval to the average time interval of the overall sequence, and forming the respiratory phase coding set by numerically encoding the phase displacement parameter.

6. The tumor respiratory motion trajectory prediction method based on spatiotemporal feature separation LSTM according to claim 4, characterized in that, The specific steps for S3 are as follows: S301: Obtain the weighted spatial feature vector sequence, pair them according to the time sequence order of the feature vectors at adjacent time points, perform one-dimensional convolution calculation on elements of the same dimension and extract the difference value, aggregate the convolution output according to the time step sequence and update the feature index to generate a temporal difference feature sequence. S302: Call the time-series differential feature sequence, extract time slice feature components, construct time weight parameters based on time step difference, perform weighting operation, adjust the position weight ratio based on time step index, rearrange the weighting result and integrate the feature structure to generate a time local change feature sequence; S303: Call the time local change feature sequence, input it into the unidirectional long short-term memory unit according to the time step sequence, calculate the transition difference between hidden states and perform weighted fusion with the output vector to generate a time recursive enhanced feature vector sequence.

7. The tumor respiratory motion trajectory prediction method based on spatiotemporal feature separation LSTM according to claim 6, characterized in that, The specific steps of S4 are as follows: S401: Based on the time recursive enhanced feature vector sequence, perform linear mapping operation on the time step feature components of the sequence, perform dimension calibration according to the mapping matrix parameters, perform linear combination of the time step matrix components according to the mapping coefficients and check the consistency to generate a unified dimension feature set; S402: Call the time step feature vector in the unified dimension feature set, perform weighted summation according to the fusion weight ratio, filter continuous frames in combination with the fusion coefficient threshold, extract window frames according to the time index of the weighted feature set, and generate a time window feature sequence. S403: Based on the continuous frame data of the time window feature sequence, input the fully connected layer, perform element-wise difference calculation on the numerical vector, establish an autoregressive prediction relationship based on the time step index, rearrange the output layer displacement vector sequence, and generate the tumor respiratory motion trajectory prediction result.

8. The tumor respiratory motion trajectory prediction method based on spatiotemporal feature separation LSTM according to claim 1, characterized in that, The method also includes step S5: S5: Call the tumor respiratory motion trajectory prediction results, compare them with the coordinates collected in real time by the radiotherapy equipment sensor, calculate the Euclidean distance deviation value, count the proportion of deviations exceeding the preset deviation threshold, adjust the fusion weight coefficient, re-execute the prediction, and generate adaptive correction trajectory prediction results; The adaptive correction trajectory prediction results include the fusion weight adjustment coefficient, the prediction error correction amount, and the trajectory deviation distribution characteristics.

9. The tumor respiratory motion trajectory prediction method based on spatiotemporal feature separation LSTM according to claim 8, characterized in that, The specific steps of S5 are as follows: S501: Call the predicted coordinate data in the tumor respiratory motion trajectory prediction result, match it with the target tissue three-dimensional motion coordinates collected in real time by the radiotherapy equipment sensor, extract the corresponding coordinate component values ​​and perform difference analysis to generate a coordinate deviation distance sequence set; S502: Call the coordinate deviation distance sequence set, compare the deviation distance value of the time node in the sequence with the deviation threshold, identify and mark the time node whose deviation exceeds the threshold, calculate the proportion of nodes exceeding the limit, organize the deviation change characteristics according to the time distribution, and generate the accuracy deviation proportion distribution. S503: Based on the accuracy deviation ratio distribution, adjust the multi-source data fusion weight coefficients in the trajectory prediction results, use the adjusted weights to correct the prediction parameter set, calculate the trajectory point coordinates, and generate adaptive correction trajectory prediction results.

10. A tumor respiratory motion trajectory prediction system based on spatiotemporal feature separation LSTM, characterized in that, The system is used to implement the tumor respiratory motion trajectory prediction method based on spatiotemporal feature separation LSTM as described in any one of claims 1-9, and the system comprises: The image time sequence acquisition module acquires chest and abdominal image sequence data during the patient's respiratory cycle using a four-dimensional computed tomography (CT) scanner. After aligning the phase image frames based on B-spline non-rigid registration transformation, the tumor centroid coordinate values ​​are extracted. The tumor centroid coordinate sequence is generated through sliding window standardization and then transmitted to the centroid feature extraction module. The centroid feature extraction module calls the standardized tumor centroid coordinate sequence, inputs it into a bidirectional long short-term memory network to extract temporal features and calculates Euclidean distance feature values, filters coordinate points that exceed a preset distance threshold to obtain respiratory phase encoding values, performs dimensionality consistency and weights the coordinate components to generate a weighted spatial feature vector sequence, and passes it to the spatial feature construction module. The spatial feature construction module calls the weighted spatial feature vector sequence to perform a one-dimensional convolution operation to calculate the feature difference between adjacent time moments, and obtains the local temporal change features by weighting. It then inputs the recursive features extracted by the unidirectional long short-term memory network and performs weighted superposition to generate a temporal recursive enhanced feature vector sequence, which is then passed to the temporal feature recursion module. The time feature recursive module calls the time recursive enhanced feature vector sequence to perform a linear transformation mapping to a unified dimension, then performs a weighted summation according to the fusion weight ratio, extracts the feature sequence of the current time window and inputs it into the fully connected layer to perform autoregressive prediction, outputs the tumor respiratory motion trajectory prediction result, and passes it to the adaptive correction module. The adaptive correction module calls the tumor respiratory motion trajectory prediction results, compares them with the coordinates collected in real time by the radiotherapy equipment sensors, calculates the Euclidean distance deviation value, counts the proportion of deviations exceeding the preset deviation threshold, adjusts the fusion weight coefficient, re-executes the prediction, and generates the adaptive correction trajectory prediction result.

Images