A tokamak plasma density profile inversion method and apparatus

Abstract

The application provides a tokamak plasma density profile inversion method and device, which adopts a CPU-GPU heterogeneous collaborative computing architecture, takes three-frequency band IQ original data measured by a microwave reflectometer as input, and constructs a collaborative working mechanism of a GPU parallel computing pipeline and a CPU serial processing pipeline. The GPU end is responsible for large-scale parallel spectrum analysis and beat frequency extraction tasks, and adopts a sliding window method and a weighted centroid algorithm to realize batch parallel processing. The CPU end is responsible for data screening, iterative inversion and post-processing tasks, adopts a three-strategy consistency voting mechanism to automatically determine an effective frequency range, and recursively solves refractive index and radial position by an iterative phase integral method. Through GPU acceleration, the time consumption of single time point calculation is reduced to less than 4 milliseconds, full-process automatic data inversion is realized, the processing timeliness is high, the processing procedure is simple, and the engineering requirements of tokamak experiment real-time feedback control are met.

Description

Technical Field

[0001] This invention relates to the field of plasma diagnostic technology, and in particular to a method and apparatus for inverting the density profile of tokamak plasma. Background Technology

[0002] Microwave reflectometers, with their advantages of high spatiotemporal resolution and non-invasive measurement, have become an important diagnostic tool for modern fusion plasma physics research and experimental devices. The plasma electron density and its spatiotemporal distribution obtained by microwave reflectometers are core diagnostic parameters for fusion research. They can characterize the plasma state, evaluate the system's confinement performance, and verify physical models, providing valuable reference for assessing plasma stability and fusion output.

[0003] The core principle of a microwave reflector is the cutoff effect of electromagnetic wave propagation in plasma. A series of microwaves of different frequencies are emitted into the plasma; when the emitted microwave frequency equals the plasma frequency, the wave will be reflected back from the corresponding critical density layer. By measuring the phase delay or frequency shift of the reflected wave, the location of the critical density layer can be calculated, thus allowing for the inversion of an electron density profile from the edge to the core.

[0004] Existing traditional density profile inversion methods based on physical models typically employ a serial CPU processing mode, requiring manual intervention in some steps to process data. This significantly increases the time required for data analysis and makes it difficult to meet real-time processing requirements. Furthermore, while data-driven AI inversion methods offer greatly improved computational speed, they lack physical constraints on the model, making reliability difficult to guarantee. Moreover, future upgrades to fusion device detection systems or significant changes in discharge scenarios will require retraining to adapt to the new environments.

[0005] Therefore, a plasma density profile reconstruction model that combines physical reliability, computational efficiency, and automated processing capabilities is needed to meet the engineering requirements for real-time plasma diagnosis and feedback control in fusion devices. Summary of the Invention

[0006] This solution aims to at least address the technical problems existing in the prior art. To this end, the first aspect of this invention proposes a tokamak plasma density profile inversion method, comprising: The IQ target data in three frequency bands—Q-band, V-band, and W-band—acquired by a tokamak microwave reflector, including in-phase component data and quadrature component data. The IQ target data of the three frequency bands are input into the GPU, and the GPU performs sliding window Welch power spectral density estimation on the IQ target data of the three frequency bands, and outputs the power spectral matrix of each frequency band and each time window. The GPU calculates the beat frequency and its beat error for each time window using a weighted centroid algorithm within a preset beat frequency filtering range for each frequency band, based on the power spectrum matrix. The beat frequency and the beat error are transmitted to the CPU, which calculates the time delay data for each time window of each frequency band based on the difference between the beat frequency and the pre-stored vacuum beat frequency. The time delay data includes microwave propagation time delay and time delay error. The CPU concatenates the frequency data of the three frequency bands, the time delay data, the power spectrum matrix, and the beat frequency error in order of frequency magnitude to form a continuous full-band data sequence. The CPU determines the endpoint of the effective frequency range using three filtering strategies based on the power spectrum data, power spectrum gradient data, and beat frequency error data in the full-band data sequence. The minimum value among the endpoints of each strategy is taken as the endpoint of the global effective frequency range. The effective data interval is then extracted based on the endpoint of the global effective frequency range. The CPU uses an iterative phase integration method to calculate the radial position and electron density of each frequency layer, starting from the zero-density cutoff frequency, based on the frequency and time delay data within the effective data interval. The radial position and the electron density are mapped to a preset standard grid to output an electron density profile.

[0007] A second aspect of the present invention provides a tokamak plasma density profile inversion device, comprising: The target data acquisition module is used to acquire IQ target data in three frequency bands—Q-band, V-band, and W-band—collected by the tokamak microwave reflector; the IQ target data includes in-phase component data and quadrature component data. The power spectrum matrix calculation module is used to input the IQ target data of the three frequency bands into the GPU, and the GPU performs sliding window Welch power spectral density estimation on the IQ target data of the three frequency bands, and outputs the power spectrum matrix of each frequency band and each time window. The beat frequency calculation module is used by the GPU to calculate the beat frequency and its beat frequency error for each time window within a preset beat frequency filtering range for each frequency band, based on the power spectrum matrix. The delay calculation module is used to transmit the beat frequency and the beat error to the CPU, and the CPU calculates the delay data of each time window in each frequency band based on the difference between the beat frequency and the pre-stored vacuum beat frequency. The delay data includes microwave propagation delay and delay error. The splicing module is used by the CPU to splice the frequency data of the three frequency bands, the time delay data, the power spectrum matrix and the beat frequency error in order of frequency magnitude to form a continuous full-band data sequence. The effective interval determination module is used by the CPU to determine the endpoint of the effective frequency range by using three filtering strategies based on the power spectrum data, power spectrum gradient data and beat frequency error data in the full-band data sequence, and take the minimum value among the endpoints of each strategy as the global effective frequency range endpoint, and extract the effective data interval based on the global effective frequency range endpoint. The electron density calculation module is used by the CPU to calculate the radial position and electron density of each frequency layer by recursively calculating the frequency and time delay data within the effective data interval, starting from the zero density cutoff frequency point. The profile output module is used to map the radial position and the electron density to a preset standard grid and output an electron density profile.

[0008] A third aspect of the present invention provides an electronic device comprising a processor and a memory, wherein the memory stores at least one instruction, at least one program, a code set, or an instruction set, wherein the at least one instruction, the at least one program, the code set, or the instruction set is loaded and executed by the processor to implement the tokamak plasma density profile inversion method as described in the first aspect.

[0009] A fourth aspect of the present invention provides a computer-readable storage medium storing at least one instruction, at least one program, a code set, or an instruction set, wherein the at least one instruction, the at least one program, the code set, or the instruction set is loaded and executed by a processor to implement the tokamak plasma density profile inversion method as described in the first aspect.

[0010] The embodiments of the present invention have the following beneficial effects: This invention provides a method for inverting the density profile of tokamak plasma. The method employs a CPU-GPU heterogeneous collaborative computing architecture, using raw IQ data from the Q, V, and W bands measured by a microwave reflector as input. A collaborative working mechanism is constructed using the PyTorch framework to integrate the GPU parallel computing pipeline and the CPU serial processing pipeline. The GPU handles large-scale parallel spectrum analysis and beat frequency extraction tasks, employing the sliding window Welch method and a weighted centroid algorithm for batch parallel processing. The CPU handles data filtering, iterative inversion, and post-processing tasks, automatically determining the effective frequency range using a three-strategy consensus voting mechanism based on multi-feature fusion, and recursively solving for the refractive index and radial position layer by layer using an iterative phase integral method. This invention reduces the computation time per time point to less than 4 milliseconds through GPU acceleration, approximately 10 times faster than a pure CPU solution. It achieves fully automated data inversion, offering high processing efficiency, a simplified processing flow, and automated data filtering, meeting the engineering requirements for real-time feedback control in tokamak experiments. Attached Figure Description

[0011] Figure 1 A flowchart illustrating the steps of a tokamak plasma density profile inversion method provided in this embodiment of the invention; Figure 2 This is a schematic diagram illustrating the extraction of valid data ranges according to an embodiment of the present invention. Figure 3 A schematic diagram of an electron density profile provided in an embodiment of the present invention; Figure 4 A system architecture diagram provided for embodiments of the present invention; Figure 5 A comparison chart of batch inversion time between a traditional model and the model proposed in this invention is provided for an embodiment of the present invention. Figure 6 A graph showing the speedup ratio of a conventional model and the model proposed in this invention as a function of the number of profiles, provided for an embodiment of this invention; Figure 7 This is a structural block diagram of a tokamak plasma density profile inversion device provided in an embodiment of the present invention. Detailed Implementation

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

[0013] Hereinafter, the terms "first" and "second" are used for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the number of indicated technical features. Thus, a feature defined as "first" or "second" may explicitly or implicitly include one or more of that feature. In the description of embodiments of this disclosure, unless otherwise stated, "a plurality of" means two or more. Furthermore, the use of "based on" or "according to" implies openness and inclusiveness, because processes, steps, calculations, or other actions "based on" or "according to" one or more of the stated conditions or values ​​may in practice be based on additional conditions or beyond the stated values.

[0014] Figure 1 A flowchart illustrating the steps of a tokamak plasma density profile inversion method provided in this embodiment of the invention.

[0015] like Figure 1 As shown, the method includes the following steps: Step 101: Acquire IQ target data in three frequency bands (Q-band, V-band, and W-band) collected by the tokamak microwave reflector; the IQ target data includes in-phase component data and quadrature component data.

[0016] From the raw signals acquired by the tokamak microwave reflector, IQ data for the Q-band, V-band, and W-band are read. Here, I represents the in-phase component data, indicating the component where the signal has a 0° phase difference with the reference signal (the real part); Q represents the quadrature component data, indicating the component where the signal has a 90° phase difference with the reference signal (the imaginary part).

[0017] The data comes from three frequency bands of the EAST tokamak device: Q band (33-56GHz), V band (48-76GHz), and W band (72-110GHz). Each band contains IQ complex signals for 5 scan cycles, with a sampling frequency of 62.5MHz and an effective sweep time of 40μs.

[0018] Microwave reflectometers measure the reflection positions corresponding to different frequencies by emitting microwaves. The frequency range covered by a single frequency band is limited; only by combining the Q, V, and W bands can the complete density range from the plasma edge to the core be covered.

[0019] As an optional embodiment, step 101 includes: Step 1011: Based on the time difference between the discharge trigger time and the target time of the tokamak microwave reflector, and combined with the sampling frequency and channel data format, calculate the starting byte offset for reading the binary data file.

[0020] This step calculates the starting byte of the binary file to be read from, based on the time difference between the discharge trigger time and the current target time, combined with the sampling frequency and channel data format. This allows for precise location of the target time within the file, avoiding the reading of irrelevant data. The target time refers to the manually selected calculation time.

[0021] Specifically, the time difference = target time Discharge trigger time; Number of sampling points = Time difference × Sampling frequency; Byte offset = Number of sampling points × Number of bytes per point × Number of channels Among them, "channel data format" refers to the number of channels contained in each sampling point (e.g., 8 channels), the data type of each channel (e.g., how many bytes an integer occupies), and the I / Q two-way complex number structure.

[0022] Step 1012: Starting from the stated starting byte offset, read the raw IQ data containing at least 10 scan cycles.

[0023] The subsequent step 1013 requires peak detection to determine the starting position of the scan. Peak detection requires a certain amount of data margin to ensure reliability. If only 5 cycles are read, it may be impossible to completely capture 5 cycles due to the peak position being too early or too late. Reading 10 cycles provides sufficient buffer to ensure that 5 complete cycles can be captured regardless of where the peak appears.

[0024] Step 1013: Perform peak detection on the voltage channel signal and take the peak position as the scan start position; extract the IQ raw data for 5 complete scan cycles from the scan start position as the IQ target data.

[0025] This step performs peak detection on the voltage channel signal to find the location of the maximum voltage value, which is used as the starting position of the scan. Then, IQ data for 5 complete scan cycles is extracted from that position as the final IQ target data.

[0026] The frequency scan of a microwave reflector is periodic, with each scan cycle having a definite start time. The voltage channel signal is typically a marker signal synchronized with the scan (such as a sawtooth wave or pulse), and its peak value corresponds to the start point of the scan cycle. By detecting the voltage peak value, the start point of the scan cycle can be accurately located.

[0027] Because single-period data can be affected by noise and edge effects, while multiple periods can be used for subsequent segment averaging (frame averaging in the Welch method), improving the stability of spectral estimation, choosing 5 periods is a trade-off between accuracy and data volume.

[0028] Step 1014: Reorganize the IQ target data of the three frequency bands, Q-band, V-band and W-band, into two-dimensional complex matrices respectively; the number of rows of the two-dimensional complex matrix is ​​the number of periods, and the number of columns is the number of sampling points in each period.

[0029] Since the original data is a one-dimensional continuous byte stream, it is not convenient to analyze it by period or time window. After reorganizing it into a two-dimensional matrix: the row index corresponds to the period number (period 1 to period 5); the column index corresponds to the sampling point position within the period; each element is a complex number (I + jQ), which contains both amplitude and phase information. In this way, the original byte stream is organized into a structured complex matrix format that is convenient for subsequent GPU batch parallel processing.

[0030] The Q band, V band, and W band each correspond to a two-dimensional complex matrix, which is a complex matrix with 5 rows and 2500 columns.

[0031] Step 102: Input the IQ target data of the three frequency bands into the GPU, and the GPU performs sliding window Welch power spectral density estimation on the IQ target data of the three frequency bands, and outputs the power spectral matrix of each frequency band and each time window.

[0032] The IQ data of the three frequency bands are sent to the GPU, which calculates the power spectral density for each frequency band and each time window, and outputs a three-dimensional power spectral matrix. The specific dimensions include: number of bands, number of time windows, and number of frequency points.

[0033] The difference frequency signal (beat frequency signal) received by the microwave reflector is submerged in noise. By estimating the power spectral density, the time-domain signal can be converted to the frequency domain, and the frequency location where the signal energy is concentrated can be found, i.e., the beat frequency.

[0034] A sliding window refers to dividing a long IQ signal into many short segments with a fixed window length (e.g., 128 sampling points), with some overlap between the windows (e.g., a step size of 64 points, with 50% overlap). Each window corresponds to a time window, reflecting the spectral characteristics within that time period.

[0035] Welch is a power spectral density estimation method. Its advantages include reducing the variance of power spectral estimation by segmented averaging, making the curve smoother; reducing spectral leakage by windowing; and the fact that each segment is independent, making it suitable for GPU batch parallel processing.

[0036] The reason for using GPU processing here is that this step involves three frequency bands, multiple time windows, and multiple cycles, including a large number of independent FFT operations. The parallel cores of the GPU can process hundreds or thousands of windows at the same time, which has a significant speed-up effect.

[0037] As an optional embodiment, performing sliding window Welch power spectral density estimation on the IQ target data of the three frequency bands includes: Step 1021: Concatenate the IQ target data of the three frequency bands into a three-dimensional tensor along the batch dimension.

[0038] The three frequency bands have identical data structures, all being 5-row × 2500-column complex matrices. They only belong to different frequency bands physically, but are isomorphic in terms of data processing logic. After concatenation along the batch dimension, the GPU can perform the exact same operations on all three frequency bands without needing to write separate code for each band.

[0039] The specific dimensions of a three-dimensional tensor are: First dimension: Number of bands (3, i.e., Q / V / W); Second dimension: Number of cycles (5, i.e., 5 scan cycles); The third dimension: number of sampling points (e.g., 2500, which is the number of sampling points per cycle).

[0040] Step 1022: Perform sliding window framing on the frequency axis of the three-dimensional tensor to generate a four-dimensional window tensor. The dimensions of the four-dimensional window tensor are the number of bands, the number of time windows, the number of periods, and the window length.

[0041] The long signal (2500 sampling points) within each period is divided into multiple short segments with a fixed window length. For example, window length: 128 sampling points; sliding step size: 64 sampling points; overlap rate: 50%. In this way, a 2500-point period can be divided into approximately 38 time windows.

[0042] The specific dimensions of a four-dimensional window tensor are: First dimension: Number of bands (3, i.e., Q / V / W); Second dimension: Number of time windows (e.g., 38); Third dimension: Number of cycles (5, i.e., 5 scan cycles); Fourth dimension: Window length (e.g., 128).

[0043] Step 1023: Merge the period number dimension and window length dimension of the four-dimensional window tensor to obtain a two-dimensional flattened tensor.

[0044] Specifically, "merging" involves flattening the "band" and "time window" dimensions in the row direction, with each row representing a specific combination of (band and time window); and flattening the "period" and "window length" dimensions in the column direction, with each column representing a sampling point within that window.

[0045] For example, the original four-dimensional window tensor is: [number of bands, number of time windows, number of periods, window length] = [3, 38, 5, 128]; the merged two-dimensional flattened tensor is: [number of bands × number of time windows, number of periods × window length] = [3 × 38, 5 × 128] = [114, 640].

[0046] In this way, the data is organized into a two-dimensional matrix with "one window per row", providing a standard input format for subsequent batch signal processing operators.

[0047] Step 1024: Perform mean centering, Hanning windowing, fast Fourier transform, and power spectral density calculation sequentially on the data of each window in the two-dimensional flattened tensor to output a three-dimensional power spectrum matrix; the dimensions of the three-dimensional power spectrum matrix are the number of bands, the number of time windows, and the number of frequency points.

[0048] Here, mean-neutralization refers to subtracting the complex mean from each window to eliminate the DC component. This is because the DC component generates a large peak at the zero frequency of the spectrum, interfering with subsequent beat frequency extraction.

[0049] Adding a Hanning window involves multiplying each window by a periodic Hanning window function. The Hanning window gradually decays to zero at both ends, which can reduce spectral leakage caused by signal truncation.

[0050] The Fast Fourier Transform (FFT) refers to performing a 256-point complex FFT on a windowed signal to convert the time-domain signal to the frequency domain and obtain a two-sided spectrum.

[0051] Power spectral density calculation involves taking the square of the modulus of the FFT result and dividing it by a normalization factor to obtain the power spectral density. The normalization factor is the product of the window function energy and the sampling frequency.

[0052] Since each combination (band, time window) retains data for 5 periods, the power spectrum of these 5 periods is averaged to obtain the final power spectrum estimate for that time window, and the three-dimensional power spectrum matrix is ​​output.

[0053] Step 103: The GPU calculates the beat frequency and its beat error for each time window using a weighted centroid algorithm within the preset beat frequency filtering range of each frequency band, based on the power spectrum matrix.

[0054] In this step, the GPU calculates the beat frequency and its error for each time window within a preset beat frequency filtering range for each frequency band, based on the power spectrum matrix obtained in the previous step.

[0055] In a microwave reflector, the transmitted and received signals are mixed to generate a difference frequency signal, the frequency of which is called the beat frequency. The magnitude of the beat frequency is related to the propagation delay of microwaves in the plasma.

[0056] On a power spectrum curve, the beat frequency does not necessarily correspond to a single peak, but rather a peak with a certain width. The weighted centroid algorithm calculates the beat frequency by weighting the frequencies of all valid signal points according to their power spectrum values. The beat frequency error is then calculated based on the frequency standard deviation.

[0057] Since each band and each time window requires an independent weighted centroid calculation, which is a batch parallel task, a GPU is used for processing.

[0058] Specifically, in the above process, the window length of the sliding window framing is 128 sampling points, the sliding step size is 64 sampling points, and the overlap rate of adjacent windows is 50%; the number of points for the fast Fourier transform is 256; the power spectral density estimation adopts the Welch method and uses a periodic Hanning window function; the power spectral density value is obtained by dividing the square of the Fourier transform modulus of the signal within the window by the product of the window function energy and the sampling frequency.

[0059] As an optional embodiment, step 103 includes: Step 1031: Based on the frequency scanning ranges of the Q band, the V band, and the W band, construct the beat frequency filtering ranges for the three frequency bands respectively, and generate corresponding Boolean mask tensors based on the beat frequency filtering ranges to identify the frequency points falling within the beat frequency filtering ranges of each frequency band.

[0060] This step sets the beat frequency filtering range for each of the three frequency bands (Q, V, and W) based on their respective frequency scanning ranges, and generates a Boolean mask tensor to identify which frequency points fall within the filtering range.

[0061] The microwave reflector emits frequencies linearly within each band, and the reflected beat frequency signal should also be within a certain frequency range. Points outside this range are usually noise, interference, or higher harmonics and cannot be considered valid beat frequencies. Therefore, a reasonable beat frequency filtering range needs to be set for each band, retaining only the frequencies within that range.

[0062] A Boolean mask is a tensor of the same length as the frequency axis, where each element is either True or False. True: The frequency point falls within the beat frequency filtering range; False: The frequency point is not within the beat frequency filtering range.

[0063] The reason for using a Boolean mask is that the filtering ranges of the three bands are different, and the power spectrum matrix is ​​three-dimensional (band × time window × frequency point). Using a Boolean mask can be applied to all time windows simultaneously, keeping the tensor dimension unchanged, which is convenient for parallel processing by the GPU.

[0064] Step 1032: Within the beat frequency filtering range of each frequency band, calculate the maximum power spectrum value, and use 1 / e times the maximum power spectrum value as the signal threshold.

[0065] The maximum power spectrum value represents the location of the strongest signal energy within that time window of that band, and is the "peak" of the beat frequency. Using it as a benchmark, a relative threshold can be adaptively set without the need to manually set an absolute threshold.

[0066] The signal threshold is set to 1 / e times the maximum power spectrum value. Frequency points where the power spectrum exceeds this threshold are considered "signal significant" points; points where the power spectrum is below this threshold are considered noise or minor components and are not included in the beat frequency calculation.

[0067] In this way, based on the filtering range, the power threshold is used to further filter out the frequency points that actually contribute to the signal.

[0068] Step 1033: Select the frequency points that simultaneously fall within the beat frequency filtering range and have a power spectral value greater than the signal threshold as valid signal points.

[0069] This step will use frequency points that simultaneously meet both conditions as valid signal points: (1) It falls within the beat frequency filtering range (i.e., the Boolean mask in step 1031 is True). (2) The power spectrum value is greater than the signal threshold (i.e. the threshold in step 1032).

[0070] Step 1034: Calculate the beat frequency based on the frequency-weighted average value of each effective signal point, where the weight is the power spectral value of the corresponding frequency point.

[0071] The formula for calculating the power spectrum value is as follows:

[0072] in, f b S represents the beat frequency, and S represents the set of valid signal points. f i This represents the frequency of the i-th valid signal point. p i This represents the power spectral density value of the i-th valid signal point.

[0073] The reason for using the power spectrum value as the weight is that the larger the power spectrum value, the stronger the signal energy at that frequency point, and the greater its contribution to the beat frequency. Weighted averaging can make the beat frequency more biased towards the location where the signal energy is concentrated, rather than a simple average.

[0074] Step 1035: Calculate the beat frequency error based on the frequency standard deviation of each effective signal point, with the weight being the power spectrum value of the corresponding frequency point.

[0075] The formula for calculating beat frequency error is as follows:

[0076] Where σ represents the beat frequency error, f b S represents the beat frequency, and S represents the set of valid signal points. f i This represents the frequency of the i-th valid signal point. p i This represents the power spectral density value of the i-th valid signal point.

[0077] Steps 1031-1035 first use the beat frequency filtering range to filter out points outside the frequency range, and then use the signal threshold to filter out points with insufficient intensity, obtaining a set of effective signal points. Then, using the power spectrum value as weight, the beat frequency is calculated using a weighted average, and the beat frequency error is calculated using a weighted standard deviation. The beat frequency obtained in this way not only represents the location of the signal energy concentration, but also includes uncertainty estimation, providing a reliable input for subsequent time delay calculation and data filtering.

[0078] Step 104: Transmit the beat frequency and the beat error to the CPU. The CPU calculates the time delay data for each time window of each frequency band based on the difference between the beat frequency and the pre-stored vacuum beat frequency. The time delay data includes microwave propagation time delay and time delay error.

[0079] This step transmits the beat frequency and beat error calculated by the GPU to the CPU, which then calculates the microwave propagation delay and delay error.

[0080] When microwaves propagate in plasma, their speed differs from that in a vacuum because the plasma's refractive index is less than 1. Compared to propagation in a vacuum, microwaves experience an additional delay in plasma; this additional delay is called the propagation time delay.

[0081] Vacuum beat frequency is the beat frequency of microwaves in the absence of plasma, and it serves as a reference. The difference between plasma beat frequency and vacuum beat frequency reflects the influence of plasma on microwave propagation. The time delay error can be calculated based on the combined standard deviation of the plasma beat frequency error and the vacuum beat frequency error.

[0082] This step involves reading pre-stored data, calculating formulas, and converting units. It is a logic-intensive task and is not suitable for GPU parallel computing. It should be processed using a CPU.

[0083] Specifically, the formula for calculating propagation delay is:

[0084] in, For plasma beat frequency, For vacuum beat frequency, The frequency scan rate is given, and the delay unit is given. .

[0085] The time delay error is calculated based on the combined standard deviation of the plasma beat frequency error and the vacuum beat frequency error.

[0086] Step 105: The CPU concatenates the frequency data of the three frequency bands, the time delay data, the power spectrum matrix, and the beat frequency error in order of frequency magnitude to form a continuous full-band data sequence.

[0087] The Q, V, and W bands each cover different frequency ranges (e.g., Q: 33–56 GHz, V: 48–76 GHz, W: 72–110 GHz). A single band can only retrieve a portion of the density profile; the data from the three bands must be pieced together along the frequency axis to form a continuous curve to cover the complete density range from the plasma edge to the core.

[0088] In this step, the CPU concatenates the frequency data, delay data, power spectrum matrix, and beat frequency error of the three frequency bands in ascending order of frequency to form a continuous full-band data sequence.

[0089] During splicing, the four types of data—frequency, time delay, power spectrum, and beat frequency error—are spliced ​​in pairs synchronously to ensure that any frequency point can uniquely correspond to a time delay value, a power spectrum value, and a beat frequency error value.

[0090] As an optional embodiment, step 105 includes: Step 1051: Obtain the start frequency and end frequency of the V-band.

[0091] This step reads the start and end frequencies of the V-band, which serve as the boundary reference for subsequent screening of Q-band and W-band data.

[0092] The reason for choosing the V-band as the boundary reference here is: The V-band is located in the middle, and its start and end frequencies naturally become the dividing points between the Q-band and the W-band.

[0093] Step 1052: Select data points with frequency values ​​lower than the starting frequency from the Q band as target data points in the Q band; the data in the target data points includes frequency data, time delay data, power spectrum matrix and beat frequency error.

[0094] Of all the data points in the Q band, only those with frequency values ​​lower than the starting frequency of the V band are retained as the target data points for the Q band.

[0095] In this way, the Q band only retains the low-frequency portion not yet covered by the V band, that is, the frequency lower than the starting frequency of the V band. Thus, the Q band is responsible for the lowest frequency range, and the V band continues from where the Q band ended.

[0096] Step 1053: Remove the data points from the first and last time windows of the V-band, and take the data points from the middle time window as the target data points of the V-band.

[0097] In all data points of the V-band, remove the data from the first and last time windows, and keep only the data points in the middle time window as the target data points of the V-band.

[0098] The first time window is at the beginning of the scan, when the microwave source may not be stable and the signal quality is poor; the last time window is at the end of the scan, when edge effects or truncation effects may exist; only the middle time window has stable scanning and the best signal quality.

[0099] Step 1054: Select data points with frequency values ​​greater than the end frequency from the W band as target data points for the W band.

[0100] Of all the data points in the W band, only those points with frequency values ​​greater than the end frequency of the V band are retained as the target data points for the W band.

[0101] In this way, the V-band is responsible for the mid-frequency band, and the W-band starts from where the V-band ends, with the frequency axis continuing to extend towards higher frequencies.

[0102] Step 1055: The Q-band target data points, the V-band target data points, and the W-band target data points are spliced ​​together end to end in ascending order of frequency to form a continuous full-band data sequence.

[0103] The target data points in the Q-band, V-band, and W-band are spliced ​​together end-to-end in ascending frequency order to form a continuous full-band data sequence. The spliced ​​frequency coverage extends continuously from the lowest frequency of the Q-band to the highest frequency of the W-band, without overlap or breaks.

[0104] Step 106: The CPU determines the endpoint of the effective frequency range using three filtering strategies based on the power spectrum data, power spectrum gradient data, and beat frequency error data in the full-band data sequence. The minimum value among the endpoints of each strategy is taken as the endpoint of the global effective frequency range. The effective data interval is then extracted based on the endpoint of the global effective frequency range.

[0105] Not all frequency data are reliable. At high frequencies, signals may become invalid due to plasma cutoff, increased noise, error divergence, or other reasons. This step automatically identifies the upper frequency limit of reliable data, discards invalid data, and obtains the valid data range to avoid deviations in the inversion results. The valid data range refers to the frequency range from the zero-density point to the global endpoint.

[0106] The three screening strategies are based on power spectrum threshold, power spectrum gradient stability, and beat frequency error continuity. Each strategy evaluates data reliability from three independent physical dimensions: signal power level, gradient change characteristics, and measurement confidence. Taking the minimum of the three endpoints signifies the most conservative judgment—if any strategy deems a piece of data unreliable, it is truncated at that point. This reduces the probability of misjudgment based on a single criterion and improves robustness.

[0107] Figure 2 This is a schematic diagram illustrating the extraction of valid data ranges according to an embodiment of the present invention.

[0108] like Figure 2As shown in the figure, taking the data from gun number 155037 at 4.6 seconds as an example, the application effect of the effective data automatic screening strategy based on multi-feature fusion is demonstrated. The shaded area in the figure represents the automatically identified effective data interval. The beat frequency power spectrum within this interval has the characteristics of concentrated signal energy, stable gradient change, and controllable confidence error. The results show that this strategy can identify and extract beat frequency signal intervals without manual intervention, improving the robustness and accuracy of real-time density profile inversion.

[0109] As an optional embodiment, the method employs three screening strategies to determine the endpoint of the effective frequency range, including: Step 1061: Obtain the maximum power spectrum value in the power spectrum matrix of each time window, normalize the maximum power spectrum value corresponding to each time window, and obtain the normalized power spectrum sequence.

[0110] The power spectrum matrix is ​​three-dimensional (band × time window × frequency point). For each time window, the power spectrum curve values ​​vary greatly at different frequency points. Therefore, the maximum value represents the strongest point of the signal in that time window and serves as the reference for normalization.

[0111] The absolute values ​​of the power spectrum can vary significantly across different time windows and frequency bands. Directly using absolute values ​​to set thresholds is not universally applicable. Therefore, the maximum power spectrum value is normalized. After normalization, all power spectrum values ​​are mapped to the [0, 1] interval, allowing for a unified threshold setting. This eliminates dimensional differences between different time windows and frequency bands, preparing for a unified threshold determination.

[0112] The specific normalization method is to divide the power spectrum value of each time window by the maximum power spectrum value of that time window to obtain the normalized power spectrum sequence.

[0113] Step 1062: Using the median of the last 20% of the normalized power spectrum sequence as the noise baseline threshold, scan from the low-frequency end to the high-frequency end. When the power spectrum values ​​of three consecutive frequency points are all lower than the noise baseline threshold, take the window position before the first frequency point among the three frequency points as the first endpoint.

[0114] The reason for taking the median of the last 20% of the data in the sequence is that the last 20% of the data is located at the high frequency end, which is usually a noise-dominated region; taking the median instead of the average value can avoid the influence of extreme values; this median represents the "typical level of noise" and serves as a noise baseline threshold.

[0115] The reason for choosing three frequency points is that if only a single frequency point is considered below the threshold and the signal is truncated, misjudgment may occur due to local noise fluctuations. Requiring three consecutive frequency points to be below the threshold indicates that this is not an occasional fluctuation, but rather that the signal has indeed been submerged in noise.

[0116] When the power spectrum values ​​of three consecutive frequency points are all lower than the noise baseline threshold, the window position before the first of the three frequency points is taken as the first endpoint. The first endpoint represents the position where the signal power level drops below the noise baseline, indicating that "the signal becomes unreliable from this frequency point onwards".

[0117] Step 1063: Calculate the first gradient of the power spectrum. Use 0.08 times the minimum gradient value as the low gradient threshold. Within a sliding window with a length of 5 frequency points, select stable windows where the gradient values ​​are all lower than the low gradient threshold and the standard deviation is less than 0.1. Search for the first stable window in reverse from the high frequency end to the low frequency end, and use the position of the first stable window as the second endpoint.

[0118] The first gradient of the power spectrum reflects how fast the power spectrum curve changes. In reliable signal regions, the power spectrum changes relatively smoothly; in noisy regions or regions where the signal disappears, the power spectrum may fluctuate dramatically. Regions with small and stable gradients indicate good signal quality.

[0119] Using 0.08 times the minimum gradient value as the low gradient threshold means that only regions with very small gradients are retained, i.e., "flat regions" where the power spectrum hardly changes.

[0120] Within a sliding window of 5 frequency points, a stable window is selected where the gradient values ​​are all below the low gradient threshold and the standard deviation is less than 0.1. This indicates that the region formed by these 5 points has a low and stable gradient and is a reliable "stable segment".

[0121] The reason for reverse retrieval from high-frequency to low-frequency data is that high-frequency data represents the least reliable region. The first stable window found by searching from high-frequency to low-frequency data marks the boundary between the unreliable and reliable regions. Data beyond this boundary is reliable.

[0122] The second endpoint represents the point where the power spectrum gradient changes from violent fluctuations to a smooth and stable state, signifying that "the data enters a stable and reliable region from this frequency point."

[0123] Step 1064: Based on the screening of the continuity of the beat frequency error, the beat frequency error is divided into a low error zone, a medium error zone, and an abnormally high error zone. Scanning is performed from the low frequency end to the high frequency end. If an abnormally high error zone or three consecutive medium error zones are detected, the window position before the abnormally high error zone or the three medium error zones is taken as the third endpoint.

[0124] The beat frequency error is the uncertainty estimate calculated in step 1035, and it directly reflects the confidence level of the beat frequency measurement. A small error indicates that the measurement is reliable, while a large error indicates that the measurement is unreliable.

[0125] Scanning from low frequency to high frequency, if an abnormally high error region or three consecutive medium error regions are detected, the window position preceding the abnormally high error region or the three medium error regions is taken as the third endpoint. Thus, the third endpoint represents the point where the beat frequency measurement confidence changes from reliable to unreliable, signifying that "the data is no longer reliable from this frequency point onwards."

[0126] Step 1065: Take the minimum value among the first endpoint, the second endpoint, and the third endpoint as the global effective frequency range endpoint.

[0127] The three strategies described above judge data reliability from different perspectives: the first strategy judges whether the signal power is sufficient, the second strategy judges whether the power spectrum gradient is stable, and the third strategy judges whether the beat frequency error is reliable.

[0128] If any of the three strategies deems a certain point of data unreliable, then that point should be truncated. Taking the minimum value among the first, second, and third endpoints as the endpoint of the global effective frequency range means adopting the most conservative judgment; the earliest occurrence of the three endpoints is the global endpoint.

[0129] Step 107: The CPU uses the iterative phase integration method to calculate the radial position and electron density of each frequency layer by layer, starting from the zero density cutoff frequency point, based on the frequency and time delay data within the effective data interval.

[0130] The iterative phase integration method is a density profile inversion algorithm based on microwave reflector measurement data. Starting from the zero-density cutoff frequency, it first integrates the measured propagation delay along the frequency axis to obtain the measured phase; then, it calculates the refractive index by combining the plasma frequency and electron cyclotron frequency, and integrates along the radial path to obtain the calculated phase; finally, it recursively deduces the radial position layer by layer based on the difference between the measured phase and the calculated phase, and finally obtains the electron density distribution corresponding to each frequency layer.

[0131] The physical principle of the iterative phase integration method is as follows: when microwaves propagate in plasma, the phase change is related to the refractive index, which in turn depends on the electron density. By measuring the phase of the microwaves, the electron density and the reflection position can be deduced.

[0132] The phase integral formula is as follows:

[0133] in, For phase, To delay the transmission time, For frequency.

[0134] This step is a sequential calculation that proceeds recursively layer by layer. The result of each layer depends on the output of the previous layer, making it impossible to parallelize. At the same time, it involves complex physical formulas and branch judgments, making it suitable for CPU execution.

[0135] Step 108: Map the radial position and the electron density to a preset standard grid and output the electron density profile.

[0136] The radial positions and electron densities obtained from the inversion are scattered data with uneven frequency distribution. To facilitate subsequent analysis, display, and comparison with other diagnostic data, they need to be interpolated / mapped onto a uniform standard grid to obtain an electron density profile.

[0137] As an optional embodiment, the step of calculating the radial position and electron density corresponding to each frequency layer by recursively calculating the frequency and time delay data within the effective data interval using an iterative phase integration method, starting from the zero-density cutoff frequency point, includes: S1: Using the zero density point as the inversion starting layer, the radial position and zero electron density corresponding to the zero density point are used as the initial conditions of the 0th layer.

[0138] The outermost electron density of a plasma approaches zero. Cutoff reflection occurs when the microwave frequency equals the plasma frequency, and the cutoff frequency is lowest at zero density. Starting from the zero-density point means that the inversion process gradually progresses from the outermost edge of the plasma towards the core.

[0139] The radial position corresponding to the zero-density point refers to the location of the zero-density point within the device, which can be estimated from the plasma current and equilibrium configuration. The electron density at the zero-density point is 0. With these two quantities, the 0th layer is completely determined and can serve as a reliable starting point for iteration.

[0140] S2: For the time delay data within the effective data interval, perform cumulative trapezoidal numerical integration along the frequency axis from small to large to obtain the measured phase at each frequency point.

[0141] The microwave reflectometer measures the propagation time delay at each frequency point, while the inversion requires the phase. Therefore, the time delay must be accumulated along the frequency axis from small to large by a trapezoidal numerical integral. That is, starting from the lowest frequency point, the trapezoidal area between adjacent frequency points is accumulated point by point to obtain the measured phase at each frequency point. This measured phase is the true value given by the measured data, which does not depend on any density distribution assumptions, and serves as the target value to be approximated in subsequent iterations.

[0142] S3: In the current iteration layer, calculate the plasma frequency parameters based on the electron density of the current layer, calculate the electron cyclotron frequency parameters based on the magnetic field strength at the current radial position, and calculate the refractive index corresponding to each frequency point on the current known radial path according to the X-mode dispersion relation.

[0143] The propagation speed of microwaves in plasma depends on the refractive index, which is affected by both electron density and magnetic field. Therefore, in the current iteration layer, it is necessary to calculate the plasma frequency parameters based on the electron density of the current layer, calculate the electron cyclotron frequency parameters based on the magnetic field strength at the current radial position, and then calculate the refractive index corresponding to each frequency point on the current known radial path according to the X-mode dispersion relation to obtain the refractive index distribution at each point on the path, which prepares for the next step of integrating along the radial path to calculate the phase.

[0144] The formula for calculating the refractive index N is as follows:

[0145] in, Where ω is the plasma oscillation frequency, and ω is the angular frequency of the electromagnetic wave. It is the electron cyclotron frequency.

[0146] S4: Perform trapezoidal numerical integration on the refractive index along the currently known radial path to obtain the calculated phase of the current layer.

[0147] Once the refractive index of each point on the known radial path is known, a trapezoidal numerical integral can be performed on the refractive index along the radial path. That is, the refractive index integrals between two adjacent layers on the path are accumulated segment by segment to obtain the calculated phase of the current layer. This calculated phase is based on the theoretical phase change that microwaves should undergo based on the currently guessed density distribution. It is compared with the measured phase obtained in step S2, and the difference between the two will drive the next step of radial position update.

[0148] S5: Based on the difference between the measured phase and the calculated phase, the refractive index at the end of the known radial path, and the frequency parameter corresponding to the end, calculate the radial position increment of the next layer, add the radial position increment to the radial position of the current layer to obtain the radial position of the next layer, update the magnetic field parameters of the next layer based on the radial position of the next layer, and update the electron density of the next layer based on the difference between the measured phase and the calculated phase. The measured phase is compared with the calculated phase. If there is a difference between the two, it means that the current guessed position of the reflective layer does not match the actual position. It is necessary to calculate the radial position increment of the next layer based on the difference, the refractive index at the end of the current path and the corresponding frequency parameters. The radial position increment of the next layer is added to the radial position of the current layer to obtain the radial position of the next layer. Then, the magnetic field parameters at the new position are updated, and the electron density at the new position is updated based on the phase difference, thus completing one layer of recursion and providing new input conditions for the next round of iteration.

[0149] The radial position recursive formula is as follows:

[0150] Where c is the speed of light in a vacuum. For frequency microwaves from spread to The amount of phase change. For frequency microwaves at position The refractive index at that point.

[0151] S6: Repeat S3 to S5 until all frequency points within the valid data range have been traversed.

[0152] After each layer of recursion is completed, the updated radial position, magnetic field parameters, and electron density of the next layer are used as the new conditions for the current layer. The process of refractive index calculation, phase integral calculation, and radial position update is returned to step S3. This process is repeated layer by layer from the zero density point to the plasma core until all frequency points within the effective data range are traversed, thereby obtaining the estimated values ​​of radial position and electron density of each layer covering the entire effective frequency range.

[0153] S7: After the iteration is completed, the arithmetic mean of the radial position and electron density of the two adjacent layers is taken as the final output result of each layer.

[0154] Since the radial position and electron density of each layer are obtained by gradually correcting the phase difference during the iteration process, there may be numerical fluctuations between adjacent layers. In order to smooth these fluctuations and improve the stability of the output results, after the iteration is completed, the arithmetic mean of the radial position and electron density of adjacent layers is taken as the final output result of each layer, so that the density profile is smoother and more continuous, providing higher quality data for subsequent mapping to the standard grid.

[0155] Figure 3 This is a schematic diagram of an electron density profile provided in an embodiment of the present invention.

[0156] like Figure 3 As shown, the radial position-electron density (R,Ne) profile obtained by the method described in this invention has a high degree of agreement with the experimental reference data, with a correlation coefficient of 99%. This result is achieved by significantly improving the processing speed through a CPU-GPU heterogeneous collaborative computing architecture while maintaining the accuracy of the traditional physical model-driven inversion algorithm. This achieves a balance between computational efficiency and inversion accuracy, meeting the engineering application requirements for real-time diagnostics of fusion plasma.

[0157] This invention employs a CPU-GPU heterogeneous collaborative computing architecture. It uses raw IQ data (Q, V, W bands) measured by a microwave reflector in a tokamak discharge experiment as input. A collaborative working mechanism between the GPU's parallel computing pipeline and the CPU's serial processing pipeline is constructed using the PyTorch framework, enabling fully automated inversion from the raw IQ signal to the plasma density profile. In this architecture, the GPU handles large-scale parallel spectrum analysis and beat frequency extraction, while the CPU handles data filtering, iterative inversion, and post-processing. The two systems achieve efficient collaboration through tensor data transmission.

[0158] Compared with existing technologies, this invention reduces the computation time per time point to milliseconds, or even less than 4ms, through GPU-accelerated batch power spectrum analysis, thus meeting the timeliness requirements of real-time feedback control in tokamak experiments. It achieves full automation by automatically determining the effective frequency range using a three-strategy consensus voting mechanism, eliminating the need for manual intervention and improving the efficiency and objectivity of data processing.

[0159] Figure 4 This is a system architecture diagram provided for an embodiment of the present invention.

[0160] like Figure 4 As shown, the system described in this invention adopts a CPU-GPU heterogeneous collaborative computing architecture. After data input, the GPU parallel processing module and the CPU serial processing module perform collaborative calculations, and finally output the results, mapping the inverted radial position R and electron density Ne to a standard 108-point grid.

[0161] The GPU parallel processing module is built on the PyTorch framework and includes a batch spectrum analysis submodule and a beat frequency extraction submodule. This module receives three-band IQ data tensors and outputs the beat frequency sequence and its error estimate for each time window through parallel operations such as sliding window framing, FFT transformation, power spectrum calculation, and weighted centroid analysis.

[0162] The CPU serial processing module includes a latency calculation submodule, a data filtering submodule, and a density inversion submodule. This module receives beat frequency data output by the GPU and performs tasks such as physical quantity conversion, data filtering, and iterative inversion.

[0163] Figure 5 A comparison chart showing the batch inversion time of a traditional model and the model proposed in this invention, provided as an embodiment of the present invention.

[0164] like Figure 5As shown in the example of experimental data from the EAST tokamak device, the inversion results for processing experimental gun number 156666 demonstrate that the traditional physics model-driven CPU serial inversion method takes 32.2 milliseconds per time point, the AI ​​inference model takes approximately 4 milliseconds, while the method described in this invention takes only about 3.6 milliseconds per time point. Even after introducing a Monte Carlo error analysis method with 20 random perturbations, the processing time of the method described in this invention can still be controlled to around 4 milliseconds.

[0165] Figure 6 The graph shows the speedup ratio of a conventional model and the model proposed in this invention as a function of the number of profiles, provided in an embodiment of this invention.

[0166] like Figure 6 As shown, the vertical axis represents the speedup ratio, which is the time taken by the traditional method divided by the time taken by the method of this invention. The horizontal axis represents the number of profiles. When the number of profiles is small, the speedup ratio is only 9.17 times. This is because the data volume is too small, and the parallel cores of the GPU are not fully utilized. However, as the number of profiles increases, the speedup ratio rises rapidly to about 30 times, and remains stable between 28 and 31 times at 5000, 10000, 15000, and 20000 profiles.

[0167] In other words, compared with the traditional CPU serial algorithm, the method described in this invention can achieve a 9.17-fold speedup in single-time-point inversion calculation. The algorithm speedup effect continues to improve with the increase of the number of parallel computing time points, and the final speedup ratio is maintained at around 30 times.

[0168] It is understood that this invention is not limited to the EAST device. For other tokamak devices equipped with microwave reflector systems, the sliding window parameters, filter range, and screening threshold can be adjusted according to their frequency range, sampling rate, number of bands, and geometric parameters. The CPU-GPU heterogeneous architecture can be extended to other plasma diagnostic systems requiring real-time spectrum analysis.

[0169] In summary, the tokamak plasma density profile inversion method provided in this invention employs a CPU-GPU heterogeneous collaborative computing architecture. It uses raw IQ data from the Q, V, and W bands measured by a microwave reflector as input, and utilizes the PyTorch framework to construct a collaborative working mechanism between the GPU parallel computing pipeline and the CPU serial processing pipeline. The GPU handles large-scale parallel spectrum analysis and beat frequency extraction tasks, employing the sliding window Welch method and weighted centroid algorithm for batch parallel processing. The CPU handles data filtering, iterative inversion, and post-processing tasks, automatically determining the effective frequency range using a three-strategy consensus voting mechanism based on multi-feature fusion, and recursively solving for the refractive index and radial position layer by layer using an iterative phase integral method. This invention reduces the computation time per time point to less than 4 milliseconds through GPU acceleration, approximately 10 times faster than a pure CPU solution. It achieves fully automated data inversion, offering high processing efficiency, a simple processing flow, and automated data filtering, meeting the engineering requirements for real-time feedback control in tokamak experiments.

[0170] Figure 7 This is a structural block diagram of a tokamak plasma density profile inversion device provided in an embodiment of the present invention.

[0171] like Figure 7 As shown, the device includes: The target data acquisition module 201 is used to acquire IQ target data in three frequency bands—Q-band, V-band, and W-band—collected by the tokamak microwave reflector; the IQ target data includes in-phase component data and quadrature component data. The power spectrum matrix calculation module 202 is used to input the IQ target data of the three frequency bands into the GPU, and the GPU performs sliding window Welch power spectral density estimation on the IQ target data of the three frequency bands, and outputs the power spectrum matrix of each frequency band and each time window. The beat frequency calculation module 203 is used by the GPU to calculate the beat frequency and its beat frequency error of each time window within the preset beat frequency filtering range of each frequency band based on the power spectrum matrix. The delay calculation module 204 is used to transmit the beat frequency and the beat error to the CPU, and the CPU calculates the delay data of each time window of each frequency band based on the difference between the beat frequency and the pre-stored vacuum beat frequency. The delay data includes microwave propagation delay and delay error. The splicing module 205 is used by the CPU to splice the frequency data of the three frequency bands, the time delay data, the power spectrum matrix and the beat frequency error in order of frequency magnitude to form a continuous full-band data sequence. The effective interval determination module 206 is used by the CPU to determine the endpoint of the effective frequency range by means of three filtering strategies based on the power spectrum data, power spectrum gradient data and beat frequency error data in the full frequency band data sequence, and take the minimum value among the endpoints of each strategy as the global effective frequency range endpoint, and extract the effective data interval based on the global effective frequency range endpoint; The electron density calculation module 207 is used by the CPU to calculate the radial position and electron density of each frequency layer by using the iterative phase integration method to calculate the frequency and time delay data within the effective data interval, starting from the zero density cutoff frequency point and recursively calculating the radial position and electron density of each frequency layer. The profile output module 208 is used to map the radial position and the electron density to a preset standard grid and output an electron density profile.

[0172] Regarding the apparatus in the above embodiments, the specific manner in which each module performs its operation has been described in detail in the embodiments related to the method, and will not be elaborated upon here.

[0173] In the above embodiments, implementation can be achieved, in whole or in part, through software, hardware, firmware, or any combination thereof. When implemented in software, it can be implemented, in whole or in part, as a computer program product. The computer program product includes one or more computer instructions. When the computer program instructions are loaded and executed on a computer, all or part of the processes or functions described in the embodiments of the present invention are generated. The computer can be a general-purpose computer, a special-purpose computer, a computer network, or other programmable device. The computer instructions can be stored in a computer-readable storage medium or transmitted from one computer-readable storage medium to another. For example, the computer instructions can be transmitted from one website, computer, server, or data center to another website, computer, server, or data center via wired (e.g., coaxial cable, fiber optic, digital subscriber line (DSL)) or wireless (e.g., infrared, wireless, microwave, etc.) means. The computer-readable storage medium can be any available medium accessible to a computer or a data storage device such as a server or data center that integrates one or more available media. The available medium can be a magnetic medium (e.g., floppy disk, hard disk, magnetic tape), an optical medium (e.g., DVD), or a semiconductor medium (e.g., solid-state disk (SSD)).

[0174] It should be noted that, in this document, relational terms such as "first" and "second" are used only to distinguish one entity or operation from another, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or apparatus. Without further limitations, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes said element.

[0175] The various embodiments in this specification are described in a related manner. Similar or identical parts between embodiments can be referred to mutually. Each embodiment focuses on describing the differences from other embodiments. In particular, the system embodiments are basically similar to the method embodiments, so the description is relatively simple; relevant parts can be referred to the descriptions of the method embodiments.

[0176] The above description is merely a preferred embodiment of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention are included within the scope of protection of the present invention.

Claims

1. A method for inverting the density profile of a tokamak plasma, characterized in that, include: The IQ target data in three frequency bands—Q-band, V-band, and W-band—acquired by a tokamak microwave reflector, including in-phase component data and quadrature component data. The IQ target data of the three frequency bands are input into the GPU, and the GPU performs sliding window Welch power spectral density estimation on the IQ target data of the three frequency bands, and outputs the power spectral matrix of each frequency band and each time window. The GPU calculates the beat frequency and its beat error for each time window using a weighted centroid algorithm within a preset beat frequency filtering range for each frequency band, based on the power spectrum matrix. The beat frequency and the beat error are transmitted to the CPU, which calculates the time delay data for each time window of each frequency band based on the difference between the beat frequency and the pre-stored vacuum beat frequency. The time delay data includes microwave propagation time delay and time delay error. The CPU concatenates the frequency data of the three frequency bands, the time delay data, the power spectrum matrix, and the beat frequency error in order of frequency magnitude to form a continuous full-band data sequence. The CPU determines the endpoint of the effective frequency range using three filtering strategies based on the power spectrum data, power spectrum gradient data, and beat frequency error data in the full-band data sequence. The minimum value among the endpoints of each strategy is taken as the endpoint of the global effective frequency range. The effective data interval is then extracted based on the endpoint of the global effective frequency range. The CPU uses an iterative phase integration method to calculate the radial position and electron density of each frequency layer, starting from the zero-density cutoff frequency, based on the frequency and time delay data within the effective data interval. The radial position and the electron density are mapped to a preset standard grid to output an electron density profile.

2. The method according to claim 1, characterized in that, The acquisition of IQ target data in the Q-band, V-band, and W-band frequency bands collected by the tokamak microwave reflector includes: Based on the time difference between the discharge trigger time and the target time of the tokamak microwave reflector, combined with the sampling frequency and channel data format, the starting byte offset of the binary data file is calculated. Starting from the aforementioned starting byte offset, read the raw IQ data, which contains at least 10 scan cycles; Peak detection is performed on the voltage channel signal, and the peak position is used as the scan start position; IQ raw data for 5 complete scan cycles is extracted from the scan start position and used as IQ target data; The IQ target data of the three frequency bands, Q-band, V-band and W-band, are respectively reorganized into two-dimensional complex matrices; the number of rows of the two-dimensional complex matrix is ​​the number of periods, and the number of columns is the number of sampling points in each period.

3. The method according to claim 1, characterized in that, The process of performing sliding window Welch power spectral density estimation on the IQ target data of the three frequency bands includes: The IQ target data of the three frequency bands are concatenated into a three-dimensional tensor along the batch dimension; A sliding window framing operation is performed on the frequency axis of the three-dimensional tensor to generate a four-dimensional window tensor. The dimensions of the four-dimensional window tensor are the number of bands, the number of time windows, the number of periods, and the window length, respectively. The period number dimension and window length dimension of the four-dimensional window tensor are combined to obtain a two-dimensional flattened tensor. The mean centering, Hanning windowing, fast Fourier transform, and power spectral density calculation are sequentially performed on the data of each window in the two-dimensional flattened tensor to output a three-dimensional power spectral matrix; the dimensions of the three-dimensional power spectral matrix are the number of bands, the number of time windows, and the number of frequency points.

4. The method according to claim 1, characterized in that, The step of calculating the beat frequency and its error for each time window using a weighted centroid algorithm within a preset beat frequency filtering range for each frequency band, based on the power spectral matrix, includes: Based on the frequency scanning ranges of the Q band, V band, and W band, beat frequency filtering ranges for the three frequency bands are constructed respectively. Based on the beat frequency filtering ranges, corresponding Boolean mask tensors are generated to identify the frequency points that fall within the beat frequency filtering ranges of each frequency band. Within the beat frequency filtering range of each frequency band, the maximum power spectrum value is calculated, and 1 / e times the maximum power spectrum value is used as the signal threshold. The frequency points that simultaneously fall within the beat frequency filtering range and have a power spectral value greater than the signal threshold are considered as valid signal points. The beat frequency is calculated based on the frequency-weighted average of each effective signal point, where the weight is the power spectral density value of the corresponding frequency point. The beat frequency error is calculated based on the frequency standard deviation of each effective signal point, with the weight being the power spectral value of the corresponding frequency point.

5. The method according to claim 1, characterized in that, The step of concatenating the frequency data and time delay data of the three frequency bands in ascending order of frequency to form a continuous full-band data sequence includes: Obtain the start and end frequencies of the V-band; Data points with frequency values ​​lower than the starting frequency are selected from the Q-band as target data points in the Q-band; the data in the target data points include frequency data, time delay data, power spectrum matrix, and beat frequency error; Remove the data points from the first and last time windows of the V-band, and take the data points from the middle time window as the target data points for the V-band; Data points with frequency values ​​greater than the end frequency are selected from the W band and used as target data points for the W band. The Q-band target data points, V-band target data points, and W-band target data points are spliced ​​together end to end in ascending order of frequency to form a continuous full-band data sequence.

6. The method according to claim 1, characterized in that, The three screening strategies employed to determine the endpoint of the effective frequency range include: Obtain the maximum power spectral value in the power spectral matrix of each time window, normalize the maximum power spectral value corresponding to each time window, and obtain the normalized power spectral sequence. Using the median of the last 20% of the normalized power spectrum sequence as the noise baseline threshold, scanning from the low-frequency end to the high-frequency end, when the power spectrum values ​​of three consecutive frequency points are all lower than the noise baseline threshold, the window position before the first frequency point among the three frequency points is taken as the first endpoint. The first gradient of the power spectrum is calculated, and 0.08 times the minimum gradient value is used as the low threshold of the gradient. Within a sliding window with a length of 5 frequency points, stable windows with gradient values ​​all below the low threshold of the gradient and a standard deviation of less than 0.1 are selected. The first stable window is searched in reverse from the high frequency end to the low frequency end, and the position of the first stable window is used as the second endpoint. Based on the screening of the continuity of the beat frequency error, the beat frequency error is divided into a low error zone, a medium error zone, and an abnormally high error zone. The scan proceeds from the low frequency end to the high frequency end. If an abnormally high error zone or three consecutive medium error zones are detected, the window position before the abnormally high error zone or the three medium error zones is taken as the third endpoint. The minimum value among the first endpoint, the second endpoint, and the third endpoint is taken as the global effective frequency range endpoint.

7. The method according to claim 1, characterized in that, The frequency and time delay data within the effective data interval are calculated using an iterative phase integration method, starting from the zero-density cutoff frequency point and recursively layer by layer, to determine the radial position and electron density corresponding to each frequency layer, including: S1: Using the zero density point as the inversion starting layer, the radial position and zero electron density corresponding to the zero density point are used as the initial conditions of the 0th layer; S2: For the time delay data within the effective data interval, perform cumulative trapezoidal numerical integration along the frequency axis from small to large to obtain the measured phase at each frequency point; S3: In the current iteration layer, calculate the plasma frequency parameters based on the electron density of the current layer, calculate the electron cyclotron frequency parameters based on the magnetic field strength at the current radial position, and calculate the refractive index corresponding to each frequency point on the current known radial path according to the X-mode dispersion relation; S4: Perform a trapezoidal numerical integration on the refractive index along the currently known radial path to obtain the calculated phase of the current layer; S5: Based on the difference between the measured phase and the calculated phase, the refractive index at the end of the known radial path, and the frequency parameter corresponding to the end, calculate the radial position increment of the next layer, add the radial position increment to the radial position of the current layer to obtain the radial position of the next layer, update the magnetic field parameters of the next layer based on the radial position of the next layer, and update the electron density of the next layer based on the difference between the measured phase and the calculated phase. S6: Repeat S3 to S5 until all frequency points within the valid data interval have been traversed; S7: After the iteration is completed, the arithmetic mean of the radial position and electron density of the two adjacent layers is taken as the final output result of each layer.

8. A tokamak plasma density profile inversion device, characterized in that, include: The target data acquisition module is used to acquire IQ target data in three frequency bands—Q-band, V-band, and W-band—collected by the tokamak microwave reflector; the IQ target data includes in-phase component data and quadrature component data. The power spectrum matrix calculation module is used to input the IQ target data of the three frequency bands into the GPU, and the GPU performs sliding window Welch power spectral density estimation on the IQ target data of the three frequency bands, and outputs the power spectrum matrix of each frequency band and each time window. The beat frequency calculation module is used by the GPU to calculate the beat frequency and its beat frequency error of each time window within a preset beat frequency filtering range of each frequency band based on the power spectrum matrix. The delay calculation module is used to transmit the beat frequency and the beat error to the CPU, and the CPU calculates the delay data of each time window in each frequency band based on the difference between the beat frequency and the pre-stored vacuum beat frequency. The delay data includes microwave propagation delay and delay error. The splicing module is used by the CPU to splice the frequency data of the three frequency bands, the time delay data, the power spectrum matrix and the beat frequency error in order of frequency magnitude to form a continuous full-band data sequence. The effective interval determination module is used by the CPU to determine the endpoint of the effective frequency range by using three filtering strategies based on the power spectrum data, power spectrum gradient data and beat frequency error data in the full-band data sequence, and take the minimum value among the endpoints of each strategy as the global effective frequency range endpoint, and extract the effective data interval based on the global effective frequency range endpoint. The electron density calculation module is used by the CPU to calculate the radial position and electron density of each frequency layer by recursively calculating the frequency and time delay data within the effective data interval, starting from the zero density cutoff frequency point. The profile output module is used to map the radial position and the electron density to a preset standard grid and output an electron density profile.

9. An electronic device, characterized in that, The electronic device includes a processor and a memory, wherein the memory stores at least one instruction, at least one program, code set, or instruction set, and the at least one instruction, the at least one program, the code set, or the instruction set is loaded and executed by the processor to implement the tokamak plasma density profile inversion method as described in any one of claims 1-7.

10. A computer-readable storage medium, characterized in that, The storage medium stores at least one instruction, at least one program, code set, or instruction set, wherein the at least one instruction, the at least one program, the code set, or the instruction set is loaded and executed by a processor to implement the tokamak plasma density profile inversion method as described in any one of claims 1-7.

Images