Medical image virtual simulation method and system based on digital twinning

Abstract

The invention relates to the field of medical image processing, and discloses a medical image virtual simulation method and system based on digital twinning. The method comprises the following steps: acquiring original physiological data and carrying out sliding window time sequence sampling to obtain a preliminary physiological movement sequence; extracting position deviation characteristics from the preliminary physiological movement sequence to obtain a dynamic interference distribution matrix; abnormal fluctuation in the dynamic interference distribution matrix is detected, nonlinear supplementary interference samples are generated and fused, and an enhanced interference data field is obtained; performing superposition operation on standard nuclear magnetic imaging data by using the enhanced interference data field to obtain a preliminary simulation image sequence; extracting time dynamic features of the preliminary simulation image sequence and a real clinical artifact image, and comparing the time dynamic features to generate a feature difference chart; and according to the feature difference image, calculating an adaptive fusion weight and performing residual feature completion to obtain a final medical image simulation result. According to the method, high-fidelity reproduction and closed-loop optimization of the dynamic artifacts of the medical image are realized.

Description

Technical Field

[0001] This invention relates to the field of medical image processing, and in particular to a method and system for virtual simulation of medical images based on digital twins. Background Technology

[0002] Currently, with the rapid development of medical imaging technology, magnetic resonance imaging (MRI), with its excellent soft tissue contrast and non-ionizing radiation characteristics, has become a key tool for clinical disease diagnosis and pathological research. Therefore, computer simulation-based virtual generation technology for medical images has become a current research hotspot.

[0003] In a current technique, a linear simulation method based on geometric transformations is typically used to generate motion artifact images. This type of method assumes that human organs are rigid bodies and simulates motion effects by performing simple translations and rotations on a clear image in the image domain, or by adding phase to the data in the frequency domain (K-space). For example, a sine function is used to simulate the periodicity of respiratory movements and directly applied to the K-space filling trajectory, thereby perturbing the final image generation process. While this method is computationally simple, its physical model is overly idealized, simplifying complex physiological movements to a single-dimensional linear displacement, ignoring the non-rigid deformations of human tissues during movement (such as compression and stretching) and the nonlinear decay characteristics of magnetic resonance signals over time.

[0004] Existing technologies suffer from problems such as insufficient realism in generated simulation images and limited texture features. Summary of the Invention

[0005] This invention provides a method and system for virtual simulation of medical images based on digital twins, in order to solve the problems of insufficient realism and monotonous texture features in the generated simulation images in the prior art.

[0006] In a first aspect, to address the aforementioned technical problems, the present invention provides a method for virtual simulation of medical images based on digital twins, comprising: Raw physiological data is acquired, and sliding window temporal sampling is performed on the raw physiological data to obtain a preliminary physiological movement sequence; The preliminary physiological motion sequence is input into a preset convolutional neural network for multidimensional feature mapping, and the positional deviation features in the time dimension are extracted and amplified to obtain a dynamic interference distribution matrix. The pixel offset amplitude in the dynamic interference distribution matrix is ​​detected. If the pixel offset amplitude exceeds a preset amplitude threshold, a nonlinear supplementary interference sample is generated and fused into the dynamic interference distribution matrix to obtain an enhanced interference data field. Standard MRI data is acquired, the enhanced interference data field is mapped and superimposed onto the pixel matrix of the standard MRI data to obtain a preliminary simulated image sequence containing multiple time frames. Acquire real clinical artifact images, extract deep texture features from the preliminary simulated image sequence and the real clinical artifact images, compare them, and generate a feature difference map reflecting spatial distribution differences. The adaptive fusion weights are calculated based on the feature difference map. The feature difference map is then weighted and mapped back to the preliminary simulation image sequence based on the adaptive fusion weights to complete the residual features, thus obtaining the final medical image simulation result.

[0007] Secondly, the present invention provides a medical image virtual simulation system based on digital twins, comprising: The physiological sequence acquisition module is used to acquire raw physiological data and perform sliding window temporal sampling on the raw physiological data to obtain a preliminary physiological movement sequence. The interference matrix generation module is used to input the preliminary physiological motion sequence into a preset convolutional neural network for multidimensional feature mapping, extract and amplify the positional deviation features in the time dimension, and obtain a dynamic interference distribution matrix. The interference data enhancement module is used to detect the pixel offset amplitude in the dynamic interference distribution matrix. If the pixel offset amplitude exceeds a preset amplitude threshold, a nonlinear supplementary interference sample is generated and fused into the dynamic interference distribution matrix to obtain an enhanced interference data field. The physical simulation overlay module is used to acquire standard nuclear magnetic resonance imaging data, map the enhanced interference data field, and overlay it onto the pixel matrix of the standard nuclear magnetic resonance imaging data to obtain a preliminary simulation image sequence containing multiple time frames. The differential feature analysis module is used to acquire real clinical artifact images, extract deep texture features of the preliminary simulation image sequence and the real clinical artifact images, compare them, and generate a feature difference map reflecting spatial distribution differences. The final simulation refinement module is used to calculate adaptive fusion weights based on the feature difference map, and then weight and map the feature difference map back to the preliminary simulation image sequence based on the adaptive fusion weights to complete the residual features and obtain the final medical image simulation result.

[0008] Compared with the prior art, the present invention has the following beneficial effects: (1) This invention acquires raw physiological data and performs sliding window temporal sampling, and combines convolutional neural networks to extract and amplify positional deviation features to construct a digital twin motion model that conforms to the real physiological mechanism of the human body. This method can accurately capture the nonlinear continuous change law of physiological signals such as breathing and heartbeat in the time dimension, and map one-dimensional physiological data into a multi-dimensional dynamic interference distribution matrix. It effectively solves the problem of motion feature distortion caused by traditional methods based only on simple geometric transformations or linear rigid body assumptions, and significantly improves the physical consistency and dynamic coherence of the source of medical image simulation. (2) This invention detects abnormal fluctuation regions in the dynamic interference distribution matrix and introduces a generative adversarial network to generate nonlinear supplementary interference samples for pixel-level fusion, thereby achieving deep enhancement of complex artifact textures. This hybrid mechanism of "physical rule guidance + deep learning filling" can automatically supplement high-frequency texture details that conform to the laws of real fluid dynamics for areas with strong interference, effectively overcoming the limitation that simple physical modeling is difficult to simulate extreme or complex motion artifacts (such as blood flow pulsation artifacts), eliminating the artificial sense of synthesis in the simulation results, and thus achieving high-quality medical artifact image generation; (3) This invention extracts and compares the deep temporal dynamic feature vectors of the preliminary simulated image sequence with those of real clinical artifact images, and uses the calculated feature difference map and adaptive fusion weights to complete the residual features of the simulation results, thus constructing a closed-loop feedback mechanism from preliminary simulation to fine calibration. This mechanism can quantify and automatically fill the domain gap between the simulated images and real clinical data in terms of spatial distribution and texture features, ensuring that the final output medical image simulation results are not only visually highly realistic, but also consistent with the clinical gold standard in terms of statistical properties, providing high-confidence data support for the training of medical AI diagnostic models and equipment calibration. Attached Figure Description

[0009] Figure 1 This is a schematic diagram of the process of the medical image virtual simulation method based on digital twin provided in the first embodiment of the present invention.

[0010] Figure 2 This is a schematic diagram of the structure of a medical imaging virtual simulation system based on digital twins provided in the second embodiment of the present invention. Detailed Implementation

[0011] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0012] Reference Figure 1 The first embodiment of the present invention provides a method and system for virtual simulation of medical images based on digital twins, including the following steps: S11, acquire raw physiological data, perform sliding window temporal sampling on the raw physiological data to obtain a preliminary physiological movement sequence; S12, the preliminary physiological motion sequence is input into a preset convolutional neural network for multidimensional feature mapping, and the positional deviation features in the time dimension are extracted and amplified to obtain a dynamic interference distribution matrix; S13, detect the pixel offset amplitude in the dynamic interference distribution matrix. If the pixel offset amplitude exceeds a preset amplitude threshold, generate a nonlinear supplementary interference sample and fuse the nonlinear supplementary interference sample into the dynamic interference distribution matrix to obtain an enhanced interference data field. S14, acquire standard MRI imaging data, map the enhanced interference data field, and superimpose it onto the pixel matrix of the standard MRI imaging data to obtain a preliminary simulation image sequence containing multiple time frames; S15, acquire real clinical artifact images, extract deep texture features of the preliminary simulation image sequence and the real clinical artifact images, compare them, and generate a feature difference map reflecting spatial distribution differences. S16, Calculate the adaptive fusion weights based on the feature difference map, and then weight and map the feature difference map back to the preliminary simulation image sequence based on the adaptive fusion weights to complete the residual features and obtain the final medical image simulation result.

[0013] In step S11, raw physiological data is acquired, and sliding window temporal sampling is performed on the raw physiological data to obtain a preliminary physiological movement sequence, including: Obtain raw physiological data, and extract amplitude numerical points by discretizing the raw physiological data according to a preset sampling time interval; Discrete physiological data segments are obtained by sliding and truncating the amplitude value points using a data segmentation window. Local feature vectors are extracted from the discrete physiological data segments, and multiple local feature vectors are aligned and connected according to the overlapping area of ​​the data segmentation window to construct a motion trend map; The motion trend mapping is interpolated and fitted to obtain the preliminary physiological motion sequence containing information on breathing and heartbeat.

[0014] It should be noted that the acquisition of raw physiological data is achieved using multimodal physiological signal acquisition equipment (such as photoplethysmography (PPG) and piezoelectric breathing straps). The reading of amplitude values ​​at preset sampling time intervals is achieved using linear interpolation resampling technology. Since the acquisition frequency of the raw data may be uneven or out of sync with the system clock, this operation first constructs a uniform time axis based on the preset sampling time intervals. For each standard moment, the two adjacent acquisition points in the raw data stream are found. The reciprocal of the distance weight is calculated based on the time distance between the standard moment and the two acquisition points. Subsequently, the amplitude values ​​of the two acquisition points are weighted and summed to obtain the amplitude value at that standard moment. This process ensures a uniform distribution of data along the time dimension.

[0015] It should be noted that the sliding truncation using data segmentation windows is implemented using matrix slicing operations. The system sets a fixed window length and a sliding step size smaller than that window length, thereby creating overlapping regions between adjacent windows. The system moves the index on the time axis according to the sliding step size, and each time truncates a one-dimensional array of a specified window length as the discrete physiological data segment.

[0016] It should be noted that the extraction of local feature vectors and the construction of motion trend maps from the discrete physiological data segments are achieved using the overlap-add method. For each segment, the waveform vector after baseline drift removal is extracted as the local feature vector. During connection, a linearly gradual weighting strategy is used for the overlapping region of two adjacent segments. Specifically, for each sampling point in the overlapping region, the preceding segment is assigned a weight that linearly decreases from 1 to 0, while the following segment is assigned a weight that linearly increases from 0 to 1. The values ​​of the two segments at the same position are multiplied by their corresponding weights and then summed to obtain the final value at that position. This operation eliminates the discontinuities at the segmentation points, smoothly connecting all segments to form a continuous time series, i.e., the motion trend map.

[0017] It should be noted that the interpolation fitting process for the motion trend mapping is implemented using the cubic spline interpolation algorithm. This algorithm constructs a cubic polynomial within each data interval and strictly requires the continuity of the first and second derivatives at the connection nodes. This operation can eliminate the slight high-frequency noise that may be generated due to superposition and fusion, generating a smooth preliminary physiological motion sequence that retains local extremum features.

[0018] It is worth noting that the determination of the preset sampling time interval is based on a comprehensive analysis of the Nyquist sampling theorem and the spectral characteristics of physiological signals. The highest effective cutoff frequency of the signal is determined by performing spectral analysis on historical heartbeat (typically below 3Hz) and respiration (typically below 0.5Hz) signals. To ensure accurate reconstruction of phase information, the reciprocal of 10 times this highest effective cutoff frequency is selected as the preset value (e.g., 10ms) to ensure the precise location of the R-wave peak is captured. The determination of the window length is based on statistical analysis of historical human resting cardiac cycles. The average RR interval in the historical data is calculated, and the number of sampling points corresponding to 1.5 times this average value is selected as the window length to ensure that each window contains at least one complete cardiac cycle.

[0019] For example, a segment of raw respiratory waveform data is read from the database. The preset sampling interval is 0.1 seconds. After resampling, the first 10 amplitude values ​​are [0.20, 0.25, 0.35, 0.50, 0.65, 0.75, 0.80, 0.75, 0.65, 0.50]. The window length is set to 5 points, and the sliding step size is 3 points. The first discrete segment is [0.20, 0.25, 0.35, 0.50, 0.65]. The second discrete segment is extracted starting from the 3rd point (index 3, i.e., 0.50). Normally, in programming, the index starts from 0. Assuming subsequent data continues, the extracted segment is [0.50, 0.65, 0.75, 0.80, 0.75]. The overlap region length is 5-3=2 points. The overlapping values ​​of the first segment are [0.50, 0.65] (corresponding to the original indices 3 and 4). The overlapping values ​​of the second segment are [0.50, 0.65] (corresponding to the original indices 3 and 4; for simplicity, this example assumes the values ​​overlap, but there may be slight deviations in reality). During alignment and connection, for the first point in the overlapping region, the weight allocation is 0.67 for the preceding segment and 0.33 for the following segment (assuming linear stepping), calculating 0.50 × 0.67 + 0.50 × 0.33 = 0.50. For the second point in the overlapping region, the weight allocation is 0.33 for the preceding segment and 0.67 for the following segment, calculating 0.65 × 0.33 + 0.65 × 0.67 = 0.65. The connected sequence maintains numerical continuity. Finally, after cubic spline interpolation smoothing, a continuous curve describing respiratory fluctuations is obtained, which is the preliminary physiological movement sequence.

[0020] In step S12, the preliminary physiological motion sequence is input into a preset convolutional neural network for multidimensional feature mapping, extracting and amplifying the positional deviation features in the time dimension to obtain a dynamic interference distribution matrix, including: The preliminary physiological movement sequence is mapped into a multidimensional temporal tensor; The multidimensional temporal tensor is slid-scanned using the multi-layer convolutional kernels of the preset convolutional neural network to extract the positional deviation values ​​and obtain the original feature map. A nonlinear mapping calculation is performed on the original feature map to amplify the positional deviation value and obtain dynamically changing detailed features. The signal fluctuation intensity is calculated based on the dynamic change details, and the signal fluctuation intensity is mapped to pixel grayscale values ​​to obtain the dynamic interference distribution matrix.

[0021] It should be noted that mapping the preliminary physiological motion sequence into a multidimensional temporal tensor is achieved using dimension expansion and data stacking techniques. This operation first acquires different physiological signal channels (such as respiratory and cardiac data) from the sequence, using the one-dimensional time series of each channel as an independent base feature layer. Subsequently, a sliding time window technique is used to truncate fixed-length local sequences, and these local sequences are stacked along the new channel dimension, thereby constructing a three-dimensional tensor containing a time dimension, a signal channel dimension, and a feature depth dimension. This tensor structure can simultaneously preserve the temporal dependence of physiological signals and the correlation between multi-source signals.

[0022] It should be noted that the preset convolutional neural network can adopt an encoder structure, for example, it may sequentially include an input layer, two one-dimensional convolutional layers with a kernel size of 3, each convolutional layer followed by a ReLU activation function and a batch normalization layer, a group of dilated convolutional layers with dilation rates of 1, 2 and 4, and an output convolutional layer; this structure captures local temporal features through shallow convolutions and expands the receptive field through dilated convolutions to associate long-range periodic information, and finally the output layer maps the learned features into an intensity estimate of the positional deviation.

[0023] It is worth noting that the construction and parameter determination of the multi-layer convolutional kernels of the pre-defined convolutional neural network are based on a hierarchical receptive field matching strategy. Specifically, the network includes shallow feature extraction kernels and deep temporal correlation kernels. The size parameter of the shallow feature extraction kernel is determined based on statistical analysis of the duration and sampling frequency of the smallest effective waveform in physiological signals (such as the ECG R wave), i.e., the number of sampling points that can completely cover the smallest effective waveform is selected as the size parameter. The deep temporal correlation kernel adopts a dilated convolution structure, and its dilation rate parameter increases exponentially with the number of network layers (e.g., according to the sequence 1, 2, 4, 8) to ensure that the equivalent receptive field duration of the top layer of the network can cover at least one complete respiratory cycle. This construction method enables the convolutional kernel to accurately capture transient positional jumps and effectively extract long-term periodic fluctuation patterns.

[0024] It should be noted that the extraction of positional deviation values ​​using multi-layer convolutional kernels is implemented using one-dimensional convolution operations. This operation involves stepping a set of pre-trained filter kernels along the time axis of the multi-dimensional temporal tensor. At each position of the slide, the weights of the filter kernels and the corresponding local tensor data are calculated, and a bias term is added. This process is equivalent to using a specific set of bandpass filters to separate fluctuation components of different frequency bands from the original signal. The response values ​​of these fluctuation components in the feature space constitute the original feature map representing the probability of positional shift.

[0025] It is worth noting that the construction and training process of the pre-defined convolutional neural network is based on a "physical-data dual-driven" strategy. This network employs a fully convolutional network (FCN) architecture, primarily consisting of multiple layers of dilated 1D convolutions to expand the receptive field and capture long-term physiological dependencies. Its training dataset is generated using a high-precision MRI physical simulation platform (such as the Bloch Simulator). Specifically, the system inputs a large number of randomly generated simulated physiological signals (including respiratory and heartbeat waveforms of different frequencies) and uses a physical model to calculate the corresponding ground truth (Ground Truth) based on the true K-space phase error and image domain pixel offset. Using mean squared error (MSE) as the loss function, the network weights are iteratively updated through backpropagation until the network can accurately fit the nonlinear mapping relationship from one-dimensional time-series signals to two-dimensional positional deviation distributions. The reason this network can extract positional deviation features is that, after training with massive amounts of data, its convolutional kernels automatically learn the coupling relationship between the phase and amplitude changes of physiological signals and the MRI encoding gradient.

[0026] It is worth noting that the construction and training of the pre-defined generative adversarial network is based on a local texture transfer strategy. The training dataset consists of two parts: one part is real high-frequency artifact texture patches extracted from high-resolution clinical medical images, and the other part is smoothed simulated data at the corresponding locations. The network uses Wasserstein distance as the loss function. Through adversarial training, the generator can learn the mapping relationship from low-frequency distributed feature vectors to high-frequency realistic textures, ensuring that the texture information stored in the network weights conforms to the statistical laws of biofluid dynamics.

[0027] It should be noted that the nonlinear mapping calculation and amplification of deviations on the original feature map are implemented using an exponentially weighted activation function. The system traverses every feature value in the original feature map. For background noise features with small values, a linear function is used to preserve their original amplitude or suppress them; while for significant feature values ​​exceeding a certain benchmark, an exponential function is used to perform a nonlinear transformation. This differentiated mapping strategy can significantly expand the dynamic range between the effective signal and background noise, allowing minute physiological fluctuations to be amplified into significant structured features in the feature space, thereby obtaining the aforementioned dynamic change details.

[0028] It is worth noting that the determination of the specific benchmark is based on statistical analysis of the response distribution of historical physiological data in the feature space. Specifically, the system pre-collects a large number of standardized physiological motion sequence samples and inputs them into the network to obtain the corresponding feature map set. Subsequently, histogram statistics are performed on the response values ​​of all pixels in the set to construct a probability density curve. Typically, the response of background noise follows a low-amplitude normal distribution, while effective physiological signals exhibit a high-amplitude long-tailed distribution. The system determines the specific benchmark by calculating the intersection of these two distributions or selecting the value corresponding to the cumulative probability density distribution function reaching a preset proportion (e.g., 85%). This process ensures that only statistically significant physiological fluctuation features are activated and amplified, while random noise is effectively shielded.

[0029] It should be noted that the calculation of signal fluctuation intensity and its mapping to pixel grayscale values ​​is achieved using Min-Max Normalization and grayscale quantization techniques. This operation first statistically analyzes the absolute value distribution of dynamic change details across the entire time period, calculating the L2 norm of the feature vector at each time point as the signal fluctuation intensity. Subsequently, this fluctuation intensity value is linearly mapped to a standard grayscale range (e.g., 0 to 255), generating a two-dimensional matrix. In this matrix, high grayscale value regions correspond to periods of intense physiological movement, while low grayscale value regions correspond to resting periods; this matrix is ​​the dynamic interference distribution matrix.

[0030] It is worth noting that the base parameter of the exponential function in the nonlinear mapping calculation is determined based on statistical analysis of artifact intensity in historical MRI images. The system extracts the difference signal between the historical clear image and the artifact image and calculates their signal-to-noise ratio (SNR). A magnification factor that can increase the SNR to twice the original value is selected, and the exponent base is determined by logarithmic back-calculation (e.g., set to 1.2) to ensure that the magnified features can be effectively recognized by the subsequent network without introducing oversaturation distortion.

[0031] For example, when performing nonlinear mapping calculations on the original feature map, an activation function with saturation characteristics, such as the sigmoid function or a variant of the hyperbolic tangent function, can be used. This allows smaller feature values ​​to be processed linearly or weakly nonlinearly, while larger feature values ​​are amplified and asymptotically approach a preset upper limit value. This amplifies the effective signal while preventing the output value from growing indefinitely. For instance, when the input feature value is 0.5, after such function mapping, the output value may be amplified to around 1.24; if the input value is 1.0, the output value may be limited to a certain saturation value between 3.0 and 5.0, rather than increasing indefinitely.

[0032] In step S13, the pixel offset amplitude in the dynamic interference distribution matrix is ​​detected. If the pixel offset amplitude exceeds a preset amplitude threshold, a nonlinear supplementary interference sample is generated and fused into the dynamic interference distribution matrix to obtain an enhanced interference data field, including: Scan the dynamic interference distribution matrix, extract the pixel grayscale value of each pixel in the dynamic interference distribution matrix as the pixel offset amplitude, and filter out abnormal fluctuation areas where the pixel offset amplitude exceeds the preset amplitude threshold. The distribution feature vector is extracted from the abnormal fluctuation region, and the distribution feature vector is input into the generator of the preset generative adversarial network. The nonlinear supplementary interference samples are output using the generator of the preset generative adversarial network. The nonlinear supplementary interference samples are embedded into the corresponding regions of the dynamic interference distribution matrix and then superimposed and fused at the pixel level. The enhanced interference data field is constructed based on the fused global data distribution.

[0033] It should be noted that scanning the dynamic interference distribution matrix and filtering abnormal fluctuation regions is achieved using a full image traversal and threshold comparison method. The system reads the grayscale value of each pixel in the matrix row by row and column by column, comparing it with a preset amplitude threshold. When a consecutive set of pixels is found to have values ​​higher than the threshold, a connected component labeling algorithm (such as the Two-Pass algorithm) is used to aggregate these pixels into an independent region of interest (ROI), and its coordinate range in the matrix is ​​recorded, which is the abnormal fluctuation region.

[0034] It should be noted that the extraction of the distribution feature vector and the generation of nonlinear supplementary interference samples using a generative adversarial network (GAN) are achieved through the inference mechanism of a Conditional Generative Adversarial Network (GAN). First, global average pooling or flattening is performed on the abnormal fluctuation region to extract statistical features reflecting the intensity and texture roughness of the region, forming a distribution feature vector. This vector is then used as a conditional input to a pre-trained generator network. The generator of the GAN can employ a structure of one fully connected layer followed by several deconvolutional layers. This generator network, composed of multiple stacked deconvolutional layers, can map the input low-dimensional feature vector back to a high-dimensional image space, outputting a matrix segment with the same size as the abnormal fluctuation region. When extracting the distribution feature vector, in addition to global average pooling, gray-level statistical features within the region, such as mean, variance, and gradient histogram, can be preferably combined to form a more representative feature description. This matrix segment not only retains the original intensity features but also adds random texture details that conform to real physical laws through the network's nonlinear fitting capability; this is the nonlinear supplementary interference sample.

[0035] It should be noted that the embedding and pixel-level overlay fusion of the nonlinear supplementary interference samples is achieved using a weighted blending technique. To avoid boundary artifacts caused by direct replacement, the system constructs a weight mask with the same size as the sample, where the center weight is 1 and gradually decays to 0 towards the edges. During the fusion process, for each pixel within the abnormal fluctuation region, the gray value in the original matrix is ​​multiplied by the retained weight (1 minus the mask weight), and the corresponding gray value in the supplementary interference sample is multiplied by the injected weight (i.e., the mask weight). Finally, the two results are added together to obtain the updated pixel value. The updated overall matrix, obtained by traversing all regions, constitutes the enhanced interference data field.

[0036] It is worth noting that the construction of the weight mask is based on the spatial decay law of a two-dimensional Gaussian distribution. Specifically, the system first determines the geometric center of the sample region as a reference point and sets the weight value of this reference point to the maximum value (i.e., 1). Subsequently, for each pixel position within the mask, the Euclidean distance between it and the reference point is calculated. The system determines the weight value based on this Euclidean distance, and its change follows a negative exponential decay characteristic, that is, as the distance between the pixel and the reference point increases, the corresponding weight value decreases non-linearly and smoothly. The rate of decrease is controlled by a preset diffusion parameter, thereby ensuring that the weight distribution can naturally converge to a small value close to zero when approaching the sample edge, forming a smooth transition shape similar to a bell curve.

[0037] It is worth noting that the determination of the preset diffusion parameters is based on statistical optimization analysis of the gradient continuity of historical fusion boundaries. Specifically, the system constructs a historical test fusion dataset, traverses a set of discrete diffusion coefficient candidate values ​​to conduct simulated fusion experiments, and calculates the gradient difference between sample edge pixels and background image pixels. The parameter value that minimizes the global mean of this gradient difference, i.e., achieves the smoothest visual transition (without obvious edge abruptness), is selected as the preset value (statistical results show that this value usually converges to one-quarter of the sample side length), thereby ensuring that the weight distribution can naturally converge to a small value close to zero when near the sample edge, forming a smooth transition shape similar to a bell curve.

[0038] It is worth noting that the preset amplitude threshold was determined based on statistical analysis of the intensity of typical artifacts in historical MRI scans. The system collected a large number of clinical images containing motion artifacts and extracted gray-level histograms of the artifact regions. The cumulative distribution function (CDF) of this histogram was calculated, and the gray value at which the cumulative probability reached 95% was selected as the preset amplitude threshold. This setting ensures that depth enhancement is performed only on those strongly interfering regions with significant intensity and a major impact on image quality, while ignoring weak random noise in the background, thus balancing simulation performance and computational resources.

[0039] For example, the input dynamic interference distribution matrix is ​​256×256 in size. The preset amplitude threshold is 200 (8-bit grayscale space). The system scan finds that within the coordinate range [100:110, 100:110], the pixel grayscale values ​​are all between 210 and 220, which is determined to be an abnormal fluctuation area. After extracting the feature vector of this area, it is input into the generator. The generator outputs a 10×10 nonlinear supplementary interference sample, whose internal values ​​present a texture that simulates the pulsation of real blood flow, with the values ​​distributed between 205 and 225. During fusion, a point at the center of the sample is selected, with an original value of 215, a sample value of 220, and a mask weight of 0.8. The fusion value is calculated as 215×(1-0.8)+220×0.8=43+176=219. The value of this point is updated to 219, and the surrounding values ​​are also subjected to similar smoothing processing. The final constructed enhanced interference data field adds realistic texture details while retaining the original high intensity.

[0040] In step S14, standard MRI data is acquired, the enhanced interference data field is mapped and superimposed onto the pixel matrix of the standard MRI data to obtain a preliminary simulated image sequence containing multiple time frames, including: A global mapping index is established based on the enhanced interference data field, and the interference intensity values ​​are extracted using the global mapping index to obtain a normalized interference coefficient matrix. The normalized interference coefficient matrix is ​​dimensionally aligned with the pixel matrix of the standard MRI data, and a nonlinear superposition operation is performed to obtain the intensity of the disturbed mixed signal. The intensity of the disturbed mixed signal is collected and assembled into a composite pixel data stream. The composite pixel data stream is then remapped in grayscale space to obtain the preliminary simulated image sequence.

[0041] It should be noted that acquiring standard MRI data involves loading artifact-free, high signal-to-noise ratio T1 or T2 weighted image data from a pre-set medical imaging standard database. The global mapping index based on the enhanced interference data field is established using coordinate gridding technology. The system constructs a two-dimensional coordinate system with the same size as the enhanced interference data field, traversing each coordinate point. The system reads the pixel grayscale value at that location as the interference intensity value. Then, it processes these values ​​using a Min-Max Normalization algorithm. Specifically, the system searches for the maximum and minimum interference intensity values ​​across the entire field, calculates the difference between the current value and the minimum value, and divides this difference by the range (i.e., the difference between the maximum and minimum values), thereby mapping all interference intensity values ​​to the closed interval [0,1], thus forming the normalized interference coefficient matrix.

[0042] The enhanced interference data field should be processed into a time-series interference field set according to the time length of the target simulation, for example, covering at least two complete respiratory cycles; accordingly, this sequence of interference fields is superimposed on the same standard MRI data in sequence to output a temporally continuous preliminary simulation image sequence.

[0043] It is worth noting that the pre-built medical imaging standard database is constructed in advance during the system initialization phase. First, raw head MRI scan data are extensively collected from publicly available research datasets (such as ADNI or OASIS) and hospital image archiving and communication systems (PACS). Second, an automated image quality assessment algorithm, combined with manual expert sampling, is used to rigorously screen out clean samples with a signal-to-noise ratio higher than the preset standard and confirmed to be free of visible motion artifacts. Finally, the screened samples are uniformly processed with skull dissection, field correction, and grayscale histogram matching, and then resampled to a uniform spatial resolution, thereby forming the medical imaging standard database containing standardized, high signal-to-noise ratio T1 or T2 weighted image data.

[0044] It should be noted that dimensional alignment of the normalized interference coefficient matrix with standard MRI data is achieved using a bicubic interpolation algorithm. This is because the resolution of standard MRI data (e.g., ...) is limited. This may be related to the resolution of the generated interference data field (e.g., To address the inconsistency, this operation scales the dimension of the normalized interference coefficient matrix to perfectly match the standard NMR data through interpolation. The nonlinear superposition operation is implemented using a signal modulation coupling model. This model simulates the nonlinear attenuation effect of motion interference on the magnetic resonance signal intensity. For each pixel in the image, the system first reads the original signal intensity of the standard NMR data, and then reads the corresponding normalized interference coefficient. The operation logic is as follows: the normalized interference coefficient is squared to enhance the nonlinear characteristics; then, 1 is subtracted from the squared value to obtain the attenuation factor; finally, the original signal intensity is multiplied by the attenuation factor to obtain the disturbed mixed signal intensity modulated by the motion interference.

[0045] It should be noted that the aggregation of disturbed mixed signal intensities into a composite pixel data stream is achieved using row-major scanning serialization technology. The system reads the calculated floating-point values ​​of the mixed signal intensities sequentially according to the row order of the image and stores them in a contiguous memory buffer. Gray-scale space remapping of the composite pixel data stream is achieved using dynamic range compression technology. Since the values ​​after nonlinear superposition may be floating-point numbers and the dynamic range is not fixed, this operation first counts the distribution range of values ​​in the data stream, then remaps these floating-point values ​​back to the bit depth space of standard medical images (e.g., 0-65535 in the 16-bit DICOM standard or 0-255 in the 8-bit display standard) using a linear stretching formula, and performs a rounding operation to finally generate the visualized preliminary simulation image sequence.

[0046] It is worth noting that the selection of the source database for the standard MRI data was based on a disease matching analysis of the target application scenario. The system statistically analyzed the number of samples and image quality scores of the diseases to be simulated (such as Alzheimer's disease) in various public datasets, selecting the database with the largest sample size and the highest average signal-to-noise ratio as the preset data source. The use of squares as the basis for the nonlinear operator in the nonlinear superposition operation is based on a simplified simulation of the attenuation law of the transverse magnetization vector in the Bloch equation. Through polynomial fitting analysis of the signal attenuation curves of historical real motion artifact images, it was found that the quadratic function can approximate the real signal loss trend with minimal computational cost.

[0047] For example, a pixel in standard MRI data The original grayscale value is 180 (8-bit). The corresponding enhanced interference data field has an interference intensity value of 200 at this location, a maximum value of 250 across the entire field, and a minimum value of 0. Calculate the normalized interference coefficient. Next, a nonlinear superposition operation is performed, first calculating the square of the coefficients. ; Calculate the attenuation factor ; Calculate the intensity of the mixed signal The pixel, initially bright at 180, is reduced to a darker 64.8, simulating signal loss due to motion. After grayscale remapping (assuming the mapping range remains unchanged), the output grayscale value for this pixel is 65. The entire image thus exhibits dark spots or blurred areas consistent with the distribution of interference.

[0048] In step S15, a real clinical artifact image is acquired, and the deep texture features of the preliminary simulated image sequence and the real clinical artifact image are extracted and compared to generate a feature difference map reflecting spatial distribution differences, including: Three-dimensional convolution operations are performed on the preliminary simulated image sequence and the real clinical artifact image respectively to extract deep temporal dynamic feature vectors; Calculate the feature distance vector between the deep temporal dynamic feature vector of the preliminary simulated image sequence and the reference feature centroid of the real clinical artifact image; The residual energy value of the feature distribution is calculated based on the feature distance vector, and the feature difference map is obtained by projecting the residual energy value onto the image space.

[0049] It should be noted that the real clinical artifact images are loaded from a pre-built clinical artifact gold standard database. This database is built during the system initialization phase. The construction process involves collecting a large number of clinical MRI sequences labeled as typical motion artifacts based on historical radiology data, covering different body parts (such as the head and abdomen) and different motion types (such as swallowing and breathing), and then performing desensitization and standardization processing. It is particularly important to note that, to accommodate the time dimension requirements of three-dimensional convolution operations, the preliminary simulated image sequence mentioned here refers to a dynamic image sequence containing a complete motion cycle, generated by repeatedly executing steps S12 to S14 at different time points based on the preliminary physiological motion sequence generated in step S11.

[0050] It should be noted that the extraction of deep temporal dynamic feature vectors through 3D convolution is achieved using a spatiotemporal feature extraction network (such as a variant of C3D or I3D). This operation treats the dynamic image sequence as a dimensional... The system utilizes a four-dimensional tensor (channel, time, height, width). It employs a convolutional kernel with a three-dimensional receptive field that slides simultaneously across the spatiotemporal dimension to capture the spatial texture distribution of pixels and its temporal variations. After multiple convolutional and downsampling processes, a dense feature map that preserves spatial structural information is output. Each spatial location on this feature map corresponds to a high-dimensional vector, which is the deep temporal dynamic feature vector, encoding the dynamic texture pattern of that local region in spatiotemporal space.

[0051] It should be noted that the calculation of the feature distance vector is implemented using the Euclidean distance algorithm. Prior to this, the system needs to load a pre-defined baseline feature centroid. This centroid is determined as follows: In the offline phase, the aforementioned 3D convolutional network is used to extract features from all real clinical artifact images in the database, and the extracted feature vectors are then clustered using K-means in a high-dimensional feature space. The geometric center of the main cluster containing the largest number of samples (i.e., the most typical artifact pattern) is selected as the baseline feature centroid. During online calculation, the system subtracts each local feature vector of the initial simulated image sequence element-wise from this baseline feature centroid to obtain the feature distance vector reflecting the deviation between the simulation details and the real standard.

[0052] It should be noted that the calculation of the residual energy value and the projection to obtain the feature difference map are achieved using vector norm calculation and up-sampling techniques. For the feature distance vector at each spatial location, the system calculates its L2 norm (i.e., vector length) as the residual energy value at that location. This value quantifies the difference in realism between the simulated artifacts and the actual artifacts. Because the size of the feature map (e.g., after convolutional downsampling) is significantly increased... Smaller than the original image size (e.g.) The system uses bilinear interpolation to magnify the matrix composed of residual energy values ​​to the original image resolution. This magnified matrix visually shows which regions in the image have artifact textures that are not realistic, which is the feature difference map.

[0053] It is worth noting that the determination of the temporal depth parameter of the three-dimensional convolutional kernel is based on autocorrelation analysis of the physiological motion cycle. The system calculates the autocorrelation function of physiological signals (such as respiratory waves) and determines the time span at which the signal correlation decays to half. The number of frames corresponding to this span (e.g., 16 frames) is selected as the temporal depth of the convolutional kernel to ensure that the network can capture complete local motion patterns rather than instantaneous noise.

[0054] For example, the initial simulated image sequence has a size of (16 frames). After processing by a 3D convolutional network, the output size is... The feature map has a feature vector of length 128 at each location. For each location on the feature map... eigenvectors Read the corresponding reference feature centroid Calculate the distance vector. Assuming the sum of squares of the elements of this vector is 4.0, then the residual energy value is... This location The corresponding residual energy value is 2.0. After 8x upsampling projection, the approximate region in the original image is... It was assigned a difference value close to 2.0. In the final generated feature difference map, this area is highlighted, indicating that there is a significant difference between the simulated artifacts and the real situation, and further correction is needed.

[0055] In step S16, adaptive fusion weights are calculated based on the feature difference map. The feature difference map is then weighted and mapped back to the preliminary simulation image sequence based on these adaptive fusion weights to complete residual features, resulting in the final medical image simulation result, including: A nonlinear mapping matrix is ​​constructed based on the residual data in the feature difference map; Calculate the deviation intensity of the feature difference map in different regions, and assign adaptive fusion weights based on the deviation intensity; The nonlinear mapping matrix is ​​weighted and modulated using the adaptive fusion weights to obtain a weighted and modulated nonlinear mapping matrix. The weighted and modulated nonlinear mapping matrix is ​​reprojected onto the initial simulation image sequence and superimposed to obtain the final medical image simulation result.

[0056] It should be noted that the nonlinear mapping matrix constructed based on the residual data in the feature difference map is achieved using a nonlinear activation transformation technique. The system traverses each pixel in the feature difference map and reads its residual energy value (i.e., the L2 norm calculated in S15). Since the residual energy value represents the abstract distance in the feature space, its numerical range may be large and its distribution uneven. Therefore, it needs to be mapped to a pixel-level adjustment amplitude using a nonlinear function. Specifically, the hyperbolic tangent function (Tanh) is used to construct the mapping relationship, and the formula is as follows: ,in This represents the residual energy value. The maximum adjustment range coefficient. This is the sensitivity coefficient. This operation compresses the unbounded energy value into a finite pixel adjustment range, generating a matrix. This is the nonlinear mapping matrix.

[0057] It should be noted that the calculation of bias intensity and the allocation of adaptive fusion weights are implemented using Local Variance Analysis and a Sigmoid Gain Function. The system sets a sliding window (e.g., The algorithm slides across the feature difference map and calculates the local variance of the residual data within the window, defining this variance value as the deviation intensity of the current center pixel. Then, a sigmoid function is used to map the deviation intensity to weight values ​​between 0 and 1. The mapping logic is as follows: for flat regions with small local variance (i.e., regions with good simulation results), weights close to 0 are assigned to avoid over-modification; for complex texture regions with large local variance (i.e., regions with significant simulation deviations), weights close to 1 are assigned to enhance the completion effect. The matrix formed by these weights is the adaptive fusion weight.

[0058] It is worth noting that the specific parameters of the S-shaped gain function are adaptively set based on the statistical distribution characteristics of historical simulation deviation data. Specifically, the function includes two core control parameters: an inflection point threshold parameter and a gain slope parameter. The inflection point threshold parameter is determined based on the median of the deviation intensity in the historical feature difference graph. This parameter determines the center position of the weight mapping curve; that is, when the deviation intensity reaches this value, the output weight is exactly 0.5. The gain slope parameter is calculated based on the standard deviation of the deviation intensity distribution, aiming to control the sensitivity of weight changes. The system adjusts this parameter so that within a preset neighborhood near the inflection point threshold (e.g., within one time deviation of the standard deviation), the weight value can smoothly and quickly transition from low to high values, thereby ensuring a rapid response to significant deviations while maintaining smoothness.

[0059] It should be noted that the weighted modulation and reprojection of the nonlinear mapping matrix are achieved using the Hadamard product and image arithmetic operations. The system multiplies the nonlinear mapping matrix element-wise with the adaptive fusion weight matrix to obtain the weighted modulated nonlinear mapping matrix. This matrix actually represents the high-frequency texture details or local contrast correction values ​​that need to be added to the original image. Finally, this matrix is ​​directly superimposed with the pixel matrix of the preliminary simulation image sequence generated in step S14, and a clipping operation is performed on the result to ensure that the final pixel values ​​fall within the standard grayscale range (e.g., 0-255), thereby obtaining the final medical image simulation result with realistic texture and rich details.

[0060] It is worth noting that the maximum adjustment range coefficient... The determination of the threshold was based on statistical analysis of the Just Noticeable Difference (JND) of human vision in medical images. The system used psychophysical experiments to determine the minimum perceptual threshold for grayscale changes in MRI images for radiologists, and selected five times this threshold as the threshold. (For example, a value of 20) is used to ensure that the completed details are visible without producing false artifacts. The sensitivity coefficient... The determination of is based on the probability density distribution of the residual energy values ​​in the feature difference map. The median of the residual energy values ​​is calculated such that... This ensures that the dynamic range of the mapping function in the linear region can cover most common deviations.

[0061] For example, a pixel in the feature difference map The residual energy value is 1.5. Setting parameters... Construct a nonlinear mapping matrix. Calculate this point. The local variance (bias strength) of the neighborhood is 0.8. Assume the weight mapping function is... The weights are calculated. Next, weighted modulation will be performed. The initial simulated image sequence has a grayscale value of 100 at this point. After superposition... The final output pixel value is 112. This process shows that, due to the large differences in features in this area, the system automatically superimposed approximately 12 grayscale levels of detail information, enhancing the realism of the simulation.

[0062] In summary, this invention, by acquiring multimodal physiological data and performing sliding window temporal sampling, combined with the multidimensional feature mapping technology of convolutional neural networks, constructs a dynamic motion model that conforms to human physiological mechanisms, effectively overcoming the technical challenge of poor physical consistency in traditional geometric transformation simulation methods. By combining the nonlinear sample generation and physical superposition mechanism of generative adversarial networks, high-frequency texture details are automatically supplemented for abnormal fluctuation regions, accurately reproducing complex physiological artifacts such as blood flow pulsation and eliminating the artificial synthetic feel of the simulation images. Furthermore, by utilizing the comparison feedback and adaptive fusion weight refinement of real clinical artifact images, the bottleneck of domain gap between simulation data and real data is overcome, achieving pixel-level realistic reproduction of medical image artifacts. This method can provide large-scale, high-confidence standard data support for the training of medical imaging AI models and the calibration of imaging equipment, significantly improving the robustness and clinical application value of medical image analysis technology.

[0063] Reference Figure 2 The second embodiment of the present invention provides a medical image virtual simulation system based on digital twins, comprising: The physiological sequence acquisition module is used to acquire raw physiological data and perform sliding window temporal sampling on the raw physiological data to obtain a preliminary physiological movement sequence. The interference matrix generation module is used to input the preliminary physiological motion sequence into a preset convolutional neural network for multidimensional feature mapping, extract and amplify the positional deviation features in the time dimension, and obtain a dynamic interference distribution matrix. The interference data enhancement module is used to detect the pixel offset amplitude in the dynamic interference distribution matrix. If the pixel offset amplitude exceeds a preset amplitude threshold, a nonlinear supplementary interference sample is generated and fused into the dynamic interference distribution matrix to obtain an enhanced interference data field. The physical simulation overlay module is used to acquire standard nuclear magnetic resonance imaging data, map the enhanced interference data field, and overlay it onto the pixel matrix of the standard nuclear magnetic resonance imaging data to obtain a preliminary simulation image sequence containing multiple time frames. The differential feature analysis module is used to acquire real clinical artifact images, extract deep texture features of the preliminary simulation image sequence and the real clinical artifact images, compare them, and generate a feature difference map reflecting spatial distribution differences. The final simulation refinement module is used to calculate adaptive fusion weights based on the feature difference map, and then weight and map the feature difference map back to the preliminary simulation image sequence based on the adaptive fusion weights to complete the residual features and obtain the final medical image simulation result.

[0064] It should be noted that the digital twin-based medical image virtual simulation system provided in this embodiment of the invention is used to execute all the process steps of the digital twin-based medical image virtual simulation method in the above embodiment. The working principles and beneficial effects of the two are one-to-one, so they will not be described again.

[0065] It should be noted that the system embodiments described above are merely illustrative. The units described as separate components may or may not be physically separate, and the components shown as units may or may not be physical units; that is, they may be located in one place or distributed across multiple network units. Some or all of the modules can be selected to achieve the purpose of this embodiment according to actual needs. Furthermore, in the accompanying drawings of the system embodiments provided by this invention, the connection relationships between modules indicate that they have communication connections, which can be specifically implemented as one or more communication buses or signal lines. Those skilled in the art can understand and implement this without any creative effort.

[0066] The specific embodiments described above further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above descriptions are merely specific embodiments of the present invention and are not intended to limit the scope of protection of the present invention. In particular, it should be noted that any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention for those skilled in the art.

Claims

1. A method for virtual simulation of medical images based on digital twins, characterized in that, include: Raw physiological data is acquired, and sliding window temporal sampling is performed on the raw physiological data to obtain a preliminary physiological movement sequence; The preliminary physiological motion sequence is input into a preset convolutional neural network for multidimensional feature mapping, and the positional deviation features in the time dimension are extracted and amplified to obtain a dynamic interference distribution matrix. The pixel offset amplitude in the dynamic interference distribution matrix is ​​detected. If the pixel offset amplitude exceeds a preset amplitude threshold, a nonlinear supplementary interference sample is generated and fused into the dynamic interference distribution matrix to obtain an enhanced interference data field. Standard MRI data is acquired, the enhanced interference data field is mapped and superimposed onto the pixel matrix of the standard MRI data to obtain a preliminary simulated image sequence containing multiple time frames. Acquire real clinical artifact images, extract deep texture features from the preliminary simulated image sequence and the real clinical artifact images, compare them, and generate a feature difference map reflecting spatial distribution differences. The adaptive fusion weights are calculated based on the feature difference map. The feature difference map is then weighted and mapped back to the preliminary simulation image sequence based on the adaptive fusion weights to complete the residual features, thus obtaining the final medical image simulation result.

2. The medical image virtual simulation method based on digital twins according to claim 1, characterized in that, The process of acquiring raw physiological data and performing sliding window temporal sampling on the raw physiological data to obtain a preliminary physiological movement sequence includes: Obtain raw physiological data, and extract amplitude numerical points by discretizing the raw physiological data according to a preset sampling time interval; Discrete physiological data segments are obtained by sliding and truncating the amplitude value points using a data segmentation window. Local feature vectors are extracted from the discrete physiological data segments, and multiple local feature vectors are aligned and connected according to the overlapping area of ​​the data segmentation window to construct a motion trend map; The motion trend mapping is interpolated and fitted to obtain the preliminary physiological motion sequence containing information on breathing and heartbeat.

3. The medical image virtual simulation method based on digital twins according to claim 1, characterized in that, The step of inputting the preliminary physiological motion sequence into a preset convolutional neural network for multidimensional feature mapping, extracting and amplifying the positional deviation features in the time dimension, and obtaining a dynamic interference distribution matrix includes: The preliminary physiological movement sequence is mapped into a multidimensional temporal tensor; The multidimensional temporal tensor is slid-scanned using the multi-layer convolutional kernels of the preset convolutional neural network to extract the positional deviation values ​​and obtain the original feature map. A nonlinear mapping calculation is performed on the original feature map to amplify the positional deviation value and obtain dynamically changing detailed features. The signal fluctuation intensity is calculated based on the dynamic change details, and the signal fluctuation intensity is mapped to pixel grayscale values ​​to obtain the dynamic interference distribution matrix.

4. The medical image virtual simulation method based on digital twins according to claim 1, characterized in that, The process involves detecting the pixel offset amplitude in the dynamic interference distribution matrix. If the pixel offset amplitude exceeds a preset amplitude threshold, a nonlinear supplementary interference sample is generated and fused into the dynamic interference distribution matrix to obtain an enhanced interference data field. This includes: Scan the dynamic interference distribution matrix, extract the pixel grayscale value of each pixel in the dynamic interference distribution matrix as the pixel offset amplitude, and filter out abnormal fluctuation areas where the pixel offset amplitude exceeds the preset amplitude threshold. The distribution feature vector is extracted from the abnormal fluctuation region, and the distribution feature vector is input into the generator of the preset generative adversarial network. The nonlinear supplementary interference samples are output using the generator of the preset generative adversarial network. The nonlinear supplementary interference samples are embedded into the corresponding regions of the dynamic interference distribution matrix and then superimposed and fused at the pixel level. The enhanced interference data field is constructed based on the fused global data distribution.

5. The medical image virtual simulation method based on digital twin according to claim 1, characterized in that, The process involves acquiring standard MRI imaging data, mapping the enhanced interference data field, and superimposing it onto the pixel matrix of the standard MRI imaging data to obtain a preliminary simulated image sequence containing multiple time frames, including: A global mapping index is established based on the enhanced interference data field, and the interference intensity values ​​are extracted using the global mapping index to obtain a normalized interference coefficient matrix. The normalized interference coefficient matrix is ​​dimensionally aligned with the pixel matrix of the standard MRI data, and a nonlinear superposition operation is performed to obtain the intensity of the disturbed mixed signal. The intensity of the disturbed mixed signal is collected and assembled into a composite pixel data stream. The composite pixel data stream is then remapped in grayscale space to obtain the preliminary simulated image sequence.

6. The medical image virtual simulation method based on digital twin according to claim 1, characterized in that, The process of acquiring real clinical artifact images, extracting deep texture features from the preliminary simulated image sequence and the real clinical artifact images, comparing them, and generating a feature difference map reflecting spatial distribution differences includes: Three-dimensional convolution operations are performed on the preliminary simulated image sequence and the real clinical artifact image respectively to extract deep temporal dynamic feature vectors; Calculate the feature distance vector between the deep temporal dynamic feature vector of the preliminary simulated image sequence and the reference feature centroid of the real clinical artifact image; The residual energy value of the feature distribution is calculated based on the feature distance vector, and the feature difference map is obtained by projecting the residual energy value onto the image space.

7. The medical image virtual simulation method based on digital twin according to claim 1, characterized in that, The step of calculating adaptive fusion weights based on the feature difference map, and then weighting and mapping the feature difference map back to the preliminary simulation image sequence based on the adaptive fusion weights to complete residual features, thereby obtaining the final medical image simulation result, includes: A nonlinear mapping matrix is ​​constructed based on the residual data in the feature difference map; Calculate the deviation intensity of the feature difference map in different regions, and assign adaptive fusion weights based on the deviation intensity; The nonlinear mapping matrix is ​​weighted and modulated using the adaptive fusion weights to obtain a weighted and modulated nonlinear mapping matrix. The weighted and modulated nonlinear mapping matrix is ​​reprojected onto the initial simulation image sequence and superimposed to obtain the final medical image simulation result.

8. A medical imaging virtual simulation system based on digital twins, characterized in that, include: The physiological sequence acquisition module is used to acquire raw physiological data and perform sliding window temporal sampling on the raw physiological data to obtain a preliminary physiological movement sequence. The interference matrix generation module is used to input the preliminary physiological motion sequence into a preset convolutional neural network for multidimensional feature mapping, extract and amplify the positional deviation features in the time dimension, and obtain a dynamic interference distribution matrix. The interference data enhancement module is used to detect the pixel offset amplitude in the dynamic interference distribution matrix. If the pixel offset amplitude exceeds a preset amplitude threshold, a nonlinear supplementary interference sample is generated and fused into the dynamic interference distribution matrix to obtain an enhanced interference data field. The physical simulation overlay module is used to acquire standard nuclear magnetic resonance imaging data, map the enhanced interference data field, and overlay it onto the pixel matrix of the standard nuclear magnetic resonance imaging data to obtain a preliminary simulation image sequence containing multiple time frames. The differential feature analysis module is used to acquire real clinical artifact images, extract deep texture features of the preliminary simulation image sequence and the real clinical artifact images, compare them, and generate a feature difference map reflecting spatial distribution differences. The final simulation refinement module is used to calculate adaptive fusion weights based on the feature difference map, and then weight and map the feature difference map back to the preliminary simulation image sequence based on the adaptive fusion weights to complete the residual features and obtain the final medical image simulation result.

Images