Ionized layer total electron content prediction method and equipment

Abstract

The invention discloses an ionosphere total electron content prediction method, which comprises the following steps: converting an extended ionosphere TEC data sequence from a time domain to a frequency domain by using standard time-frequency conversion, and obtaining a standard time-frequency conversion spectrum; extracting periodic signals displayed in the standard time-frequency transform spectrum by using a null method, and independently predicting the periodic signals respectively; carrying out decomposition and reconstruction on the residual non-periodic components by adopting an SSA method, and then establishing an ARIMA model for prediction; and finally, fusing the prediction results of the components to obtain a final TEC prediction sequence. According to the method disclosed by the invention, decomposition and differential modeling of the complex TEC sequence are realized, and a plurality of periodic components coupled in the TEC sequence can be effectively separated; for the aperiodic component, noise interference is suppressed and the stability of the signal is improved by means of an SSA method, so that input characteristics with higher signal-to-noise ratio and quality are provided for the ARIMA model, and the prediction capability for the aperiodic component is remarkably enhanced.

Description

Technical Field

[0001] This invention belongs to the field of ionospheric detection and space weather forecasting technology, specifically relating to a method for predicting the total electron content of the ionosphere, and also to a computer device. Background Technology

[0002] Total electron content (TEC) of the ionosphere is a crucial parameter describing ionospheric characteristics and determining the current total number of free electrons in the ionosphere. Changes in TEC cause delays and refractions in satellite signal propagation, thus affecting positioning accuracy. Furthermore, it poses a significant threat to the quality and stability of radio communications. Therefore, accurate monitoring and prediction of dynamic changes in ionospheric TEC are essential for applications in navigation, positioning, and radio communications, especially in fields heavily reliant on accurate navigation and reliable communication, such as aviation, maritime, and land transportation. With the continuous advancement of Global Navigation Satellite Systems (GNSS) and the gradual improvement of systems like GPS, BeiDou, and Galileo, not only have new development opportunities been brought to satellite navigation, timing, and positioning technologies, but related research in the field of ionospheric TEC has also been greatly promoted.

[0003] In the field of ionospheric TEC prediction research, numerous research institutions worldwide have accumulated fruitful results. Internationally, the International GNSS Service (IGS) has established hundreds of GNSS observation stations globally, utilizing station data combined with mathematical models to construct the widely influential Global Ionospheric Model (GIM). Existing technologies include well-known empirical models such as the International Ionospheric Reference (IRI) model, the Bent model, the Klobuchar model, and the NeQuick model. Domestically, research institutions represented by Wuhan University and the Chinese Academy of Sciences, relying on my country's independently developed BeiDou Global Navigation Satellite System, have independently produced and released global or regional ionospheric TEC grid products with internationally advanced levels. These research results not only provide a solid theoretical foundation and rich data resources for ionospheric TEC prediction but also provide support and guarantees for in-depth research and applications in related fields such as space weather and satellite navigation. From the perspective of the evolution of prediction methods, they have gradually developed from traditional models that rely on physical simplification and empirical parameters (such as the Bent model) and static representation models based on mathematical function expansion (such as the spherical harmonic function model) to dynamic prediction models that are data-driven, such as prediction models based on statistical learning and time series analysis.

[0004] Time series models are based on historical data patterns of TEC (Electro-Temporal Dynamics) for modeling and prediction. Common time series models include autoregressive (AR) models, autoregressive moving average (ARIMA) models, and their extensions. These models are widely used due to their clear structure and high computational efficiency. However, they are essentially limited to modeling linear and stationary components and are not accurate enough for predicting complex time series. Ionospheric TEC sequences have a mixed characteristic of multiple components, including periodic and trend components, and are also affected by various factors such as solar activity and geomagnetic disturbances, exhibiting significant non-stationarity. Moreover, there are complex nonlinear coupling relationships between different components. Directly applying traditional time series methods makes it difficult to effectively separate and characterize these signals with different characteristics, thus limiting prediction accuracy. Summary of the Invention

[0005] The purpose of this invention is to address the aforementioned problems in the prior art by providing a method for predicting the total electron content of the ionosphere, and also to provide a computer device.

[0006] The above-mentioned objectives of the present invention are achieved by the following technical means:

[0007] A method for predicting the total electron content of the ionosphere includes the following steps:

[0008] Step 1: Extract ionospheric TEC data from the global ionospheric TEC dataset for selected geographical locations and time periods to construct an ionospheric TEC data sequence, and extend the ionospheric TEC data sequence to obtain an extended ionospheric TEC data sequence.

[0009] Step 2: Perform standard time-frequency transformation on the extended ionospheric TEC data sequence to obtain the standard time-frequency transformation spectrum corresponding to each sub-period signal, and use the unemotional method to extract each sub-period signal from the standard time-frequency transformation spectrum corresponding to each sub-period signal;

[0010] Step 3: For each sub-period signal extracted in Step 2, use the periodic part prediction model based on standard time-frequency transformation to predict it. Then, superimpose the prediction results corresponding to each sub-period signal to obtain the total prediction result of the periodic part. By setting different prediction step sizes, the total prediction result of the periodic part corresponding to each prediction step size can be obtained.

[0011] Step 4: Subtract all the sub-periodic signals extracted in Step 2 from the ionospheric TEC data sequence to obtain the aperiodic residual sequence. Apply singular spectrum analysis to decompose and reconstruct the aperiodic residual sequence to obtain the aperiodic reconstructed sequence.

[0012] Step 5: By setting different prediction step sizes, the autoregressive integral moving average prediction model is used to predict the aperiodic reconstructed sequence, and the total prediction result of the aperiodic part corresponding to each prediction step size is obtained.

[0013] Step 6: Algebraically add the total prediction results of the periodic part and the prediction results of the non-periodic part for each prediction step to obtain the complete ionospheric TEC prediction sequence.

[0014] The standard time-frequency transformation in step 2 as described above is based on the following formula:

[0015] ;

[0016] In the formula, The input time signal Standard time-frequency transform spectrum For real time, For real-time frequency, For the modulus of instantaneous frequency, For window width parameters, The input time signal, For time, The imaginary unit;

[0017] Set window width parameters The scope will be expanded by substituting the extended ionospheric TEC data sequence. The standard time-frequency transform spectrum corresponding to the extended ionospheric TEC data sequence is obtained, and sub-period signals of different periods are determined in the standard time-frequency transform spectrum corresponding to the extended ionospheric TEC data sequence; the window width parameter is set. Then substitute the extended ionospheric TEC data sequence into The total standard time-frequency transform spectrum corresponding to each sub-period signal is obtained. Then, the standard time-frequency transform spectrum corresponding to the original ionospheric TEC data sequence is extracted from the total standard time-frequency transform spectrum corresponding to each sub-period signal, and the standard time-frequency transform spectrum corresponding to each sub-period signal is obtained respectively.

[0018] For real signals, each sub-period signal is extracted from the corresponding standard time-frequency transform spectrum based on the following formula:

[0019] ;

[0020] In the formula, The extracted sub-period signal, for In frequency The standard time-frequency transform spectrum at that location, for In frequency The modulus of the standard time-frequency transform spectrum at that location. The angular frequency corresponding to each sub-period signal. Standard time-frequency transform spectrum The real part, It is the arctangent function. It is a cosine function.

[0021] For complex signals, each sub-period signal is extracted from the corresponding standard time-frequency transform spectrum based on the following formula:

[0022] ;

[0023] In the formula, Standard time-frequency transform spectrum The imaginary part.

[0024] The periodic prediction model described above is based on the following formula:

[0025] ;

[0026] In the formula, For prediction step size The overall prediction results for the periodic portion, To predict the start time, To predict the step size, For trend items, This is an operation for taking the real part of a complex number. This represents the sequence number of the sub-period signal. The number of sub-periodic signals. , , The first The amplitude, frequency, and phase of each sub-cycle signal;

[0027] Take different prediction step sizes Substituting into the formula, we obtain the total prediction results for the periodic portion with different prediction step sizes.

[0028] As described above, the ionospheric TEC data sequence is extended in the following way: using the original ionospheric TEC data sequence as the source, the ionospheric TEC data sequence is copied twice, and then the two copies of the ionospheric TEC data sequence are added to both ends of the original ionospheric TEC data sequence to construct the extended ionospheric TEC data sequence.

[0029] As described above, step 5 specifically includes the following steps:

[0030] Step 5.1: First, differencing the aperiodic reconstructed sequence, then performing an ADF test on the differencing aperiodic reconstructed sequence; if the p-value of the ADF test is less than a set value, the ADF test is passed; if the p-value of the ADF test is greater than or equal to the set value, the differencing aperiodic reconstructed sequence is differencing again, and the ADF test is performed again; until... The aperiodic reconstructed sequence after the second difference passes the ADF test and yields... The aperiodic reconstructed sequence after the second difference is an aperiodic reconstructed stationary sequence. , The degree of the difference;

[0031] Step 5.2: Construct the autoregressive integral moving average prediction model. The autoregressive integral moving average prediction model is as follows:

[0032] ;

[0033] In the formula, To predict the step size is The prediction results for the non-periodic portion, for Aperiodic reconstructed stationary sequence at time 1 The value or predicted value, It is a constant term. It is the first One autoregressive coefficient, It is the autoregressive order. It is the first A moving average coefficient, It is the order of the moving average, the sequence number. Serial Number , yes Error term at time, for Error term at time;

[0034] Determine the optimal Combination, for optimal Autoregressive order in the combination and moving average order and the corresponding constant term Autoregression coefficient and moving average coefficient The final autoregressive integral moving average prediction model was constructed.

[0035] Take different prediction step sizes Substituting the results into the final autoregressive integral moving average prediction model, we obtain the prediction results for the non-periodic portion with different prediction step sizes.

[0036] As mentioned above, autoregressive order and moving average order Determined in the following ways:

[0037] Preset autoregression order and moving average order The range, traversing all Combination, for each The corresponding autoregressive integral moving average prediction models are constructed by combining these models. Then, the Akaike Information Criterion (AIC) or Bayesian Information Criterion (BIC) for each model is calculated. Finally, the model that minimizes either the AIC or BIC value is selected. Combination as optimal combination;

[0038] Based on the optimal Autoregressive order of the combination and moving average order The optimal value is obtained using the maximum likelihood estimation method. The constant term in the corresponding autoregressive integral moving average prediction model Autoregression coefficient and moving average coefficient .

[0039] As described above, step 4 specifically includes the following steps:

[0040] Step 4.1: Subtract all the sub-periodic signals extracted in Step 2 from the ionospheric TEC data sequence to obtain the aperiodic residual sequence. ,in, arrive These are the numbers 1 to 1 in the non-periodic residual sequence. Ionospheric TEC data for each hour. The length of the non-periodic residual sequence;

[0041] Step 4.2: According to the selected window length Transforming the aperiodic residual sequence into a trajectory matrix : ,in, Let be the number of columns in the trajectory matrix. , arrive These are non-periodic residual sequences. medium length is The time delay segment, the window length Set it to the period of the largest amplitude.

[0042] Step 4.3: Analyze the trajectory matrix. Perform singular value decomposition: ,in, For serial number, , For the first There are several elementary matrices, each corresponding to a component pattern in the aperiodic residual sequence. It is a matrix The Large eigenvalues and Trajectory matrices The left singular vector and the right singular vector of the matrix, for The transpose of the matrix;

[0043] Step 4.4, The elementary matrices are selectively divided into Disjoint subsets Then, the elementary matrices in each subset are summed to obtain the contribution rate matrix for each subset: ,in, For the first A subset, number , For the first Contribution rate matrix corresponding to each subset;

[0044] Then the trajectory matrix for The sum of the contribution rate matrices: ;

[0045] Step 4.5, Each contribution rate matrix is ​​obtained by diagonal averaging. One reconstructed sequence;

[0046] The reconstruction order r is determined based on the singular value contribution rate. The reconstruction sequences are sorted from highest to lowest singular value contribution rate. The first r reconstruction sequences are retained. Finally, the retained r reconstruction sequences are added together to obtain the aperiodic reconstruction sequence.

[0047] A computer device includes a memory and a processor, the memory storing a computer program, and the processor executing the computer program to implement the steps of the method for predicting the total electron content of the ionosphere as described above.

[0048] Compared with the prior art, the present invention has the following advantages:

[0049] (1) The present invention realizes the effective separation and differentiated modeling of periodic signals and non-periodic signals, avoiding mutual interference between signals with different frequency characteristics during the modeling process.

[0050] (2) This invention introduces a noise filtering and signal enhancement mechanism. By using singular spectrum analysis, non-periodic signals are decomposed and reconstructed, effectively filtering out noise and improving signal stability. This provides higher quality signal input for the autoregressive integral moving average prediction model and enhances the prediction accuracy of the non-periodic part.

[0051] (3) This invention constructs a collaborative framework of a periodic prediction model based on standard time-frequency transformation + a non-periodic prediction model based on singular spectrum analysis and autoregressive integral moving average, which gives full play to the advantages of each model and significantly improves the overall prediction accuracy and system robustness. Attached Figure Description

[0052] Figure 1 This is a flowchart of the method of the present invention;

[0053] Figure 2 A schematic diagram of the standard time-frequency transform spectrum of a sub-periodic signal with a period of 24 hours;

[0054] Figure 3 A schematic diagram of the standard time-frequency transform spectrum of a sub-periodic signal with a period of 12 hours;

[0055] Figure 4 A schematic diagram of the standard time-frequency transform spectrum of a sub-periodic signal with a period of 8 hours;

[0056] Figure 5 This is a schematic diagram of the sub-period signals extracted using the non-active method;

[0057] Figure 6 A graph showing the contribution ratio of singular values;

[0058] Figure 7 A comparison diagram of aperiodic reconstructed sequences and aperiodic residual sequences;

[0059] Figure 8 A comparison graph showing the prediction results of the ARIMA model alone, the prediction results of the prediction method of this invention, and the actual TEC sequence. Detailed Implementation

[0060] To facilitate understanding and implementation of the present invention by those skilled in the art, the present invention will be further described in detail below with reference to embodiments. The embodiments described herein are for illustration and explanation only and are not intended to limit the present invention.

[0061] Example 1:

[0062] like Figure 1 As shown, a method for predicting the total electron content of the ionosphere is presented. This invention constructs a prediction method based on NTFT (Standard Time-Frequency Transform)-SSA (Singular Spectrum Analysis)-ARIMA (Autoregressive Integral Moving Average Model), which includes the following steps:

[0063] Step 1: Obtain the global ionospheric TEC dataset. Extract ionospheric TEC data from the global ionospheric TEC dataset for selected geographical locations and time periods to construct an ionospheric TEC data sequence. Then, extend the ionospheric TEC data sequence to obtain an extended ionospheric TEC data sequence. The specific steps include:

[0064] Step 1.1, Data Acquisition and Integrity Check: Download the 2021 global ionospheric TEC dataset from the European Orbit Determination Centre (CODE) and check the global ionospheric TEC data daily to ensure that the data is complete and without missing parts.

[0065] The ionospheric TEC dataset provided by the European Orbit Determination Centre (CODE) is a global ionospheric atlas (GIMs) published in IONEX format (a standard data format for exchanging ionospheric data, especially TEC maps of total electron content in the ionosphere). One CODG file is released daily, and each CODG file includes 24 global ionospheric GIMs, one per hour. The spatial coverage of each global ionospheric GIM is: latitude 87.5°N to 87.5°S (interval 2.5°, total 71 points) and longitude 180°W to 180°E (interval 5°, total 73 points).

[0066] Step 1.2: Read the global ionospheric TEC dataset, extract ionospheric TEC data for selected geographical locations and time periods from the global ionospheric TEC dataset to construct an ionospheric TEC data sequence, specifically including the following steps:

[0067] Step 1.2.1: First, extract the ionospheric TEC dataset for the selected geographical location: Convert the latitude and longitude coordinates of the selected geographical location into the row and column indices corresponding to the Global Ionospheric Map (GIM), and then extract the ionospheric TEC dataset for the selected geographical location from the GIM based on the row and column indices. The specific conversion formula is as follows:

[0068] Latitude Index: ;

[0069] Longitude Index: ;

[0070] In the formula, Latitude Longitude It is a floor function, which means taking the largest integer not greater than x, i.e., rounding down towards negative infinity.

[0071] Step 1.2.2: Extract ionospheric TEC data for a selected time period from the ionospheric TEC dataset of the selected geographical location to construct an ionospheric TEC data sequence: This invention uses a prediction mode of "10 days of training and 5 days of prediction" for experiments. Based on this, this embodiment selects 15 consecutive ionospheric TEC data points from the ionospheric TEC dataset of the selected geographical location in 2021, and takes the first 10 days of ionospheric TEC data to construct an ionospheric TEC data sequence, denoted as... The ionospheric TEC data for the following 5 days were used for subsequent comparison and evaluation with the model prediction results, denoted as . .

[0072] To verify the effectiveness of the proposed method, in this embodiment, an ionospheric TEC data sequence at a specific geographical location was selected for the experiment. The latitude and longitude were set to 5°N, 120°W. Based on the index formula described above, the row and column index corresponding to this location in the grid is calculated as follows: Regarding the selection of time periods, the ionospheric TEC data from the period of 218 to 227 days (August 6 to August 15, 2021) constituted the ionospheric TEC data sequence, and the ionospheric TEC data from the period of 228 to 232 days (August 16 to August 20) were compared and evaluated with the model prediction results.

[0073] Step 1.3, sequence extension to suppress edge effects, specifically includes the following processes:

[0074] When processing finite-length signals, time-frequency transformation cannot obtain the complete neighborhood information required for local spectrum analysis due to data interruption at the boundary. It often requires data extension methods such as zero-padding, which leads to distortion of time-frequency representation in the boundary region, i.e., edge effect.

[0075] To avoid the impact of edge effects on signal extraction, it is necessary to process the ionospheric TEC data sequence. The mirror-symmetric extension is performed as follows: using the ionospheric TEC data sequence as the original ionospheric TEC data sequence, the ionospheric TEC data sequence is copied twice, and then the two copies of the ionospheric TEC data sequence are added to both ends of the original ionospheric TEC data sequence to construct the extended ionospheric TEC data sequence. In subsequent steps, during time-frequency transformation, a standard time-frequency transformation is performed on the extended ionospheric TEC data sequence to obtain the standard time-frequency transformation spectrum corresponding to the extended ionospheric TEC data sequence. Then, only the standard time-frequency transformation spectrum corresponding to the original ionospheric TEC data sequence is extracted to extract the signals of each sub-period, thereby effectively mitigating the impact of edge effects without changing the intrinsic characteristics of the original signal.

[0076] Step 2, as follows Figures 2-5The diagram shown illustrates the sub-period signal extraction process. A standard time-frequency transform (CTF) is performed on the extended ionospheric TEC data sequence to obtain the standard time-frequency transform spectrum corresponding to each sub-period signal. Then, an inactive method (IM) based on NTFT is used to extract each sub-period signal from the standard time-frequency transform spectrum. Specifically, the process includes the following steps:

[0077] Step 2.1: The standard time-frequency transformation (NTFT) can convert the extended ionospheric TEC data sequence in the time domain. Mapping to the time-frequency domain generates the NTFT spectrum. Specifically, this embodiment uses the standard Morlet wavelet transform (NMWT) to convert the extended ionospheric TEC data sequence into a time-frequency domain-mapped spectrum, generating the standard time-frequency transform spectrum (NTFT spectrum) corresponding to the extended ionospheric TEC data sequence. The standard time-frequency transform spectrum includes the standard time-frequency transform spectra corresponding to the 24-hour, 12-hour, and 8-hour sub-period signals, respectively. The standard time-frequency transform specifically includes the following processes:

[0078] For time signals ( (For the set of complex numbers), the expression for the standard time-frequency transform is:

[0079] (1);

[0080] In the formula, The input time signal The standard time-frequency transform spectrum (generated after standard time-frequency transform) displays the original input time signal. The sub-periodic signal present in it, Indicates real time. Indicates instantaneous frequency. For the input time signal, this embodiment To extend the ionospheric TEC data sequence, The time-domain kernel function of the standard time-frequency transform. For time, The time-domain kernel function of time-frequency transformation The conjugate operator, To shift in time The conjugate operator of the subsequent time-domain kernel function, It is the set of real numbers.

[0081] Frequency domain kernel function of standard time-frequency transform NTFT The conditions shown in equation (2) must be met:

[0082] (2);

[0083] in, Represents the frequency domain kernel function The model, For frequency, This represents the Fourier transform operator. Representing the time-domain kernel function Fourier transform, time domain After Fourier transform, it is transformed into the frequency domain. .

[0084] Time-domain kernel function of standard time-frequency transform The expression is:

[0085] (3);

[0086] in, As a time-frequency resolution regulator, it can theoretically take any value and expression except 0. For time-frequency resolution regulator The model, The imaginary unit, The window function can take different expressions, such as the time-domain kernel function. Different time-frequency transforms can be obtained by taking different expressions.

[0087] Next, we derive the NTFT transform used in this embodiment. First, we construct the time-domain kernel function. As one possible implementation, the window function and time-frequency resolution adjuster used in this embodiment are as follows:

[0088] (4);

[0089] in, Use a standard Gaussian window. Substituting the window width parameter into equation (3), we obtain the time-domain kernel function constructed in this embodiment as follows:

[0090] (5);

[0091] In the formula, For instantaneous frequency The model.

[0092] Substituting the time-domain kernel function of formula (5) into formula (1), we obtain the standard Morlet wavelet transform NMWT used in this embodiment, specifically:

[0093] (6);

[0094] Among them, window width parameter The selection of the window width parameter will affect the time resolution and frequency resolution of the standard time-frequency transform spectrum; If the value is too small, the frequency resolution will decrease, causing the signals of each period in the standard time-frequency transform spectrum to become mixed, which in turn affects the reliability of the subsequently extracted sub-period signals; conversely, the window width parameter... If the value is too large, the time resolution will be too low, the edge effect range of the standard time-frequency transform spectrum will be larger, and the edges of the subsequently extracted periodic signals will be significantly distorted.

[0095] Set window width parameters The scope will be expanded by substituting the extended ionospheric TEC data sequence. Based on formula (6), the extended ionospheric TEC data sequence is analyzed. The standard Morlet wavelet transform (NMWT) is performed to obtain the standard time-frequency transform spectrum corresponding to the extended ionospheric TEC data sequence. The standard time-frequency transform spectrum of the extended ionospheric TEC data sequence clearly shows sub-periodic signals with periods of 24 hours, 12 hours, and 8 hours. After determining the periodic signals present in the extended ionospheric TEC data sequence, the sub-periodic signals of multiple periods displayed in the standard time-frequency transform spectrum of the extended ionospheric TEC data sequence are processed into periodic segments. Appropriate window width parameters are set in formula (6). The total standard time-frequency transform spectra (TFTs) corresponding to sub-period signals with periods of 24 hours, 12 hours, and 8 hours in the extended ionospheric TEC data sequence were obtained (e.g., Figures 2-4 (as shown), for use in subsequent signal extraction.

[0096] The standard time-frequency transform spectra of the total standard time-frequency transform spectra of the sub-period signals with periods of 24 hours, 12 hours, and 8 hours in the extended ionospheric TEC data sequence were extracted from the original ionospheric TEC data sequence.

[0097] In this embodiment, the extended ionospheric TEC data sequence The window width parameters are set for the 24-hour, 12-hour, and 8-hour sub-period signals, respectively. The window width parameters are set to 80, 90, and 120 respectively. In practical applications, these parameters... Often determined based on experience, those skilled in the art can adjust the window width parameter independently.

[0098] Step 2.2: Extract each sub-period signal from the standard time-frequency transform spectrum corresponding to each sub-period signal using the non-active method (IM) based on NTFT. This specifically includes the following process:

[0099] For the standard time-frequency transform spectra corresponding to each sub-period signal obtained in step 2.1, the inaction method (IM) is used to extract the sub-period signals. The inaction method is based on the standard time-frequency transform. The maximum line of the standard time-frequency transform spectrum of a signal (the time-frequency ridge in the standard time-frequency transform spectrum) is the signal itself, without the need for inverse transform. The following is the principle and process of the inaction method to extract periodic signals from the standard time-frequency transform spectrum.

[0100] A complex harmonic signal When performing standard time-frequency transformation using formula (7), the instantaneous unbiased characteristics of formulas (8) and (9) can be satisfied.

[0101] (7);

[0102] In the formula, The imaginary unit, Represents a complex harmonic signal The amplitude, Amplitude The model, Represents a complex harmonic signal frequency, Represents a complex harmonic signal The initial phase, Indicates the instantaneous phase.

[0103] (8);

[0104] (9);

[0105] In the formula, To obtain the maximum value, For complex harmonic signals Standard time-frequency transform spectrum For complex harmonic signals The modulus of the standard time-frequency transform spectrum, For complex harmonic signals In frequency The standard time-frequency transform spectrum at that location.

[0106] Complex harmonic signals After performing the standard time-frequency transform, the modulus of the standard time-frequency transform spectrum is taken. , The maximum value of is the instantaneous amplitude of the original signal, and its derivation is as follows:

[0107]

[0108]

[0109] (10);

[0110] make ,but Substituting into formula (10), we get:

[0111]

[0112]

[0113] (11);

[0114] Combining the conditions satisfied by the frequency domain kernel function in equation (2), we know that the Fourier transform of the Gaussian window function satisfies equation (12), so when When , there is a maximum value equal to 1.

[0115] (12);

[0116] In the formula, The Fourier transform of the standard Gaussian window function, The modulus of the Fourier transform of the standard Gaussian window function.

[0117] Therefore, the derivation result can be changed to The above derivation and results show that the standard time-frequency transform spectrum of a complex harmonic signal can unbiasedly describe the instantaneous amplitude and instantaneous frequency characteristics of the original signal. The real and imaginary parts determine the instantaneous phase of the signal (arctan(real / imaginary)).

[0118] A real harmonic signal After standard time-frequency transformation, the standard time-frequency transform spectrum should be twice its modulus, i.e. It can represent the instantaneous amplitude and instantaneous frequency of the original signal. The real part of twice the standard time-frequency transform spectrum is the instantaneous phase of the original signal.

[0119] (13);

[0120] According to Euler's formula ,get:

[0121] (14);

[0122] In the formula, This represents the operation of taking the real part of a complex number.

[0123] The above describes the principle of IM signal extraction. In step 2.1, the standard time-frequency transform spectrum displays sub-signals with periods of 24 hours, 12 hours, and 8 hours. The standard time-frequency transform spectra corresponding to each of the three sub-signals are obtained through periodic segmentation processing. For the 24-hour period sub-signal, the passive method is applied, requiring input of its corresponding standard time-frequency transform spectrum, its upper and lower limits, and its period to extract the signal. Then, the passive method is applied sequentially to the other two periodic signals to extract them from the standard time-frequency transform spectrum.

[0124] In summary, the extracted sub-periodic signal from the complex harmonic signal is:

[0125]

[0126] (15);

[0127] in, These are the angular frequencies corresponding to each sub-cycle.

[0128] In this embodiment, the extracted sub-period signal is:

[0129]

[0130] (16);

[0131] in, Take the angular frequencies corresponding to periods of 24 hours, 12 hours, and 8 hours respectively, such as... Figure 5 As shown.

[0132] Step 3: For each sub-period signal extracted in Step 2, a periodic prediction model based on standard time-frequency transform is used for prediction. Then, the prediction results of each sub-period signal are superimposed to obtain the total prediction result of the periodic part. The specific process includes the following steps:

[0133] To derive the periodic partial prediction model based on the standard time-frequency transform, let... Time signal for:

[0134] (17);

[0135] Signal Including trend items The sum of the signals from each sub-period, where Represents the number of sub-period signals. This represents the sequence number of the sub-period signal. , , Representing the first The amplitude, frequency, and phase of each periodic term.

[0136] Since the signal is transformed into the complex domain after standard time-frequency transformation, it will be converted into complex form when extracting the sub-period signal. .

[0137] In actual forecasting, it is necessary to select a forecast start time. Let the forecast start time be... The predicted step size is Therefore, in the periodic prediction model based on the standard time-frequency transform, It is always there:

[0138] (18);

[0139] In the formula, For prediction step size The overall prediction result of the periodic part (i.e., in (Total prediction result for the periodic portion at each time step); taking different prediction step sizes. Substituting into formula (18), we obtain the total prediction results for the periodic portion of different prediction step lengths.

[0140] Equation (18) is the periodic part prediction model based on the standard time-frequency transform. This model tracks the instantaneous amplitude of each component. ,frequency and phase Make predictions.

[0141] In summary, this embodiment uses a periodic prediction model based on standard time-frequency transform to predict each sub-period signal. The prediction start time is the last point of the extracted sub-period signal, and the prediction step size is 1, 2, 3..., 120 hours (the sampling interval of ionospheric TEC data is one point per hour, and a total of 5 days are predicted, so a total of 120 points are predicted, and the prediction step size is...). (From 1 to 120 hours), then the prediction results of each sub-period signal are added together to obtain the total prediction result of the periodic part, denoted as . .

[0142] Step 4, as follows Figure 6 and Figure 7 The schematic diagram of SSA decomposition and reconstruction of the aperiodic residual sequence illustrates the SSA decomposition and reconstruction of the aperiodic residual sequence: from the original ionospheric TEC data sequence. Subtract all the sub-periodic signals extracted in step 2 from the original sequence to obtain the aperiodic residual sequence. Singular spectral analysis (SSA) is applied to non-periodic residual sequences to filter out noise and extract major trends and fluctuation patterns.

[0143] The process of SSA decomposition and reconstruction specifically includes the following steps:

[0144] Step 4.1: Subtract all the sub-periodic signals extracted in Step 2 from the original ionospheric TEC data sequence to obtain the aperiodic residual sequence. ,in, arrive These are the numbers 1 to 1 in the non-periodic residual sequence. Ionospheric TEC data for each hour. is the length of the non-periodic residual sequence.

[0145] Step 4.2: First, according to the selected window length... Transform the one-dimensional aperiodic residual sequence into a trajectory matrix. trajectory matrix for:

[0146] (19);

[0147] in, Trajectory matrix The number of columns, , arrive These are non-periodic residual sequences. medium length is The time delay segment; in practical applications, the window length The choice of has a significant impact on the decomposition effect. A widely used rule of thumb is to set it to the length of the expected dominant period in the data, where the expected dominant period is the period with the largest proportion (amplitude) in the extraction step, which is 24 hours in this embodiment.

[0148] Step 4.3, Next, process the trajectory matrix. Perform singular value decomposition (SVD) and use the following formula to decompose the trajectory matrix. Decomposed into the sum of a series of basic components (primitive matrices):

[0149] (20);

[0150] In the formula, each elementary matrix This corresponds to a component pattern (such as trend, periodicity, noise, etc.) in the non-periodic residual sequence, where, It is a matrix The Large eigenvalues ​​(for matrices) The eigenvalues ​​arranged in descending order (one eigenvalue) Singular values ​​measure the importance (energy) of the corresponding component. Large singular values ​​usually correspond to strong signals, such as trends and main cycles, while small singular values ​​usually correspond to noise. It is a trajectory matrix The left singular vector reflects the non-periodic residual sequence within the window length. Empirical orthogonal functions on, Trajectory matrix The right singular vector reflects the evolution of the components over time. For matrix The transpose of .

[0151] Step 4.4, then elementary matrices Selectively divided into In the grouping, the first component usually corresponds to the strongest trend component, periodic components usually appear in pairs and their corresponding singular values ​​are close, and noise components are divided by observing the inflection point of the singular value decrease in the eigenvalue spectrum. Small singular values ​​after the inflection point are often regarded as noise; each group is a disjoint subset. Each subset corresponds to a specific physical component, such as a trend component, a periodic component, or a noise component, and the elementary matrices belonging to the same subset are summed:

[0152] (twenty one);

[0153] In the formula, For the first A subset, number , For the first Contribution rate matrix corresponding to each subset.

[0154] The trajectory matrix is ​​then represented as The sum of the contribution rate matrices:

[0155] (twenty two).

[0156] Step 4.5, Finally Each contribution rate matrix is ​​converted into a reconstructed time series of the same length as the aperiodic residual sequence: since each contribution rate matrix is ​​approximately a Hankel matrix, the elements on its antidiagonal are averaged by diagonal averaging to obtain a one-dimensional sequence.

[0157] Diagonal averaging specifically refers to: for the contribution rate matrix any element in , belongs to the The anti-diagonal line. Therefore, the line located at the th... The average of all elements on the anti-diagonal lines is used as the reconstructed sequence. At any moment The value is as follows:

[0158] (twenty three);

[0159] In the formula, For reconstructing the sequence In the moment The value, For elements on the opposite diagonal.

[0160] After the Perform the above operation on each group to obtain Reconstruction sequence And satisfy:

[0161] (twenty four);

[0162] in, This is a non-periodic reconstructed sequence that has not yet had noise removed. For the corresponding trend, oscillations of different periods or noise components, the number of singular values ​​to be retained is determined according to the singular value contribution rate. The number of retained singular values ​​is denoted as r, which is also called the reconstruction order. By adjusting the value of the reconstruction order r, signals containing different components (such as trends, periods or noise) can be selectively reconstructed, thereby extracting specific signal features.

[0163] In this embodiment, the above-described SSA method is used to process the non-periodic residual sequence. Decompose and reconstruct: Window length The time interval is set to 24 hours, consistent with the identified main period length in the ionospheric TEC data sequence, to effectively capture its inherent temporal structure. The reconstruction order r is determined based on the singular value contribution rate proportion map of each group. The reconstructed sequences are sorted from highest to lowest singular value contribution rate, and the first r reconstructed sequences are retained. Finally, the retained r reconstructed sequences are added together to obtain the aperiodic reconstructed sequence, such as... Figure 6 As shown, the contribution rates of components after the 15th component approach zero and can be considered as noise. Therefore, the first 15 reconstructed components are selected for signal recovery, ultimately yielding the denoised aperiodic reconstructed sequence, denoted as . ,like Figure 7 The figure shows a comparison between the aperiodic reconstructed sequence and the aperiodic residual sequence.

[0164] Step 5, prediction of the non-periodic portion, specifically includes the following process:

[0165] For the reconstructed aperiodic reconstructed sequence Establish an autoregressive integral moving average prediction model ( The prediction model is used to make predictions. The construction of the autoregressive integral moving average prediction model specifically includes the following processes:

[0166] Step 5.1: First, differencing the aperiodic reconstructed sequence is performed. Then, a stationarity test is conducted on the differencing aperiodic reconstructed sequence, typically using the ADF test. The ADF test determines the stationarity of the sequence by comparing the test statistic with the critical value at a given significance level. If the p-value of the ADF test is less than the set value (usually 0.05), the ADF test is passed, and the differencing aperiodic reconstructed sequence is determined to be stationary. If the p-value of the ADF test is greater than or equal to the set value, the differencing aperiodic reconstructed sequence is differencing again, and the ADF test is performed again; this process continues until... The aperiodic reconstructed sequence after the second difference passes the ADF test and yields... The aperiodic reconstructed sequence after the second difference is an aperiodic reconstructed stationary sequence. , Let be the degree of the difference.

[0167] Step 5.2: Construct the autoregressive integral moving average prediction model. The expression for the autoregressive integral moving average prediction model is:

[0168] (25);

[0169] in, To predict the step size is The prediction results of the non-periodic part (i.e., in (Total prediction result of the non-periodic portion at time) for Aperiodic reconstructed stationary sequence at time 1 The value or predicted value, It is a constant term. It is the first Each autoregressive coefficient reflects the relationship between the current value and the value at a certain point in the past. It is the autoregressive order. It is the first A moving average coefficient describes the relationship between the current value and the error at a past point in time. It is the order of the moving average, the sequence number. Serial Number , yes Error term at time (white noise) for Error term for time.

[0170] In practice, determine The order method can combine grid search and information criteria: a pre-defined autoregressive order is used. and moving average order Within a reasonable range, traverse all Combination, for each The corresponding autoregressive integral moving average prediction models are constructed by combining them. Then, the Akaike Information Criterion (AIC) or Bayesian Information Criterion (BIC) for each autoregressive integral moving average prediction model is calculated. Finally, the model that minimizes the Akaike Information Criterion (AIC) or Bayesian Information Criterion (BIC) value is selected. Combination as optimal combination.

[0171] Based on the optimal Autoregressive order of the combination and moving average order The constant term in the autoregressive integral moving average prediction model can be estimated using parameter estimation methods such as maximum likelihood estimation (MLE). Autoregression coefficient and moving average coefficient Parameter estimation is performed to construct the autoregressive integral moving average prediction model, and different prediction step sizes are selected. Substituting the results into the final autoregressive integral moving average prediction model, we obtain the prediction results for the non-periodic portion with different prediction step sizes.

[0172] The order of the final autoregressive integral moving average prediction model established in this embodiment is: .

[0173] Step 6: Algebraically sum the total prediction results of the periodic portion of each prediction step with the prediction results of the non-periodic portion to obtain the complete ionospheric TEC prediction sequence, denoted as […]. .

[0174] To verify the superiority of the prediction method proposed in this invention, this embodiment compares the prediction results of the prediction method of this invention with the prediction results of a standalone ARIMA model:

[0175] First, an ionospheric TEC data series from day 228 to 232 of the year was constructed to directly establish a standalone ARIMA model for overall prediction. After stationarity testing, model order determination, and parameter estimation, the optimal standalone ARIMA model was determined to be [model name missing]. The prediction results of the standalone ARIMA model are denoted as .

[0176] like Figure 8 As shown, the prediction results of the ARIMA model alone, the prediction results of the method of this invention, and the real TEC sequences are compared. Figure 8In the figure, the complete ionospheric TEC prediction sequence (blue solid line) obtained by the prediction method of the present invention is closer to the real TEC sequence (black solid line) in terms of phase matching and detail characterization than the prediction result of the ARIMA model alone (red solid line). The prediction result of the prediction method of the present invention is obviously better than the prediction result of the ARIMA model alone.

[0177] In one embodiment, a computer device is also provided, including a memory and a processor, wherein the memory stores a computer program, and the processor executes the computer program to implement the steps in the above method embodiments.

[0178] In one embodiment, a computer-readable storage medium is provided having a computer program stored thereon that, when executed by a processor, implements the steps in the above method embodiments.

[0179] In one embodiment, a computer program product is provided, including a computer program that, when executed by a processor, implements the steps in the above method embodiments.

[0180] It should be noted that the embodiments described in this invention are merely illustrative of the spirit of the invention. Those skilled in the art to which this invention pertains can make various modifications or additions to the described embodiments or use similar methods to substitute them, without departing from the spirit of the invention or exceeding the scope defined by the appended claims.

Claims

1. A method for predicting the total electron content of the ionosphere, characterized in that, Includes the following steps: Step 1: Extract ionospheric TEC data from the global ionospheric TEC dataset for selected geographical locations and time periods to construct an ionospheric TEC data sequence, and extend the ionospheric TEC data sequence to obtain an extended ionospheric TEC data sequence. Step 2: Perform standard time-frequency transformation on the extended ionospheric TEC data sequence to obtain the standard time-frequency transformation spectrum corresponding to each sub-period signal, and use the unemotional method to extract each sub-period signal from the standard time-frequency transformation spectrum corresponding to each sub-period signal; Step 3: For each sub-period signal extracted in Step 2, use the periodic part prediction model based on standard time-frequency transformation to predict it. Then, superimpose the prediction results corresponding to each sub-period signal to obtain the total prediction result of the periodic part. By setting different prediction step sizes, the total prediction result of the periodic part corresponding to each prediction step size can be obtained. Step 4: Subtract all the sub-periodic signals extracted in Step 2 from the ionospheric TEC data sequence to obtain the aperiodic residual sequence. Apply singular spectrum analysis to decompose and reconstruct the aperiodic residual sequence to obtain the aperiodic reconstructed sequence. Step 5: By setting different prediction step sizes, the autoregressive integral moving average prediction model is used to predict the aperiodic reconstructed sequence, and the total prediction result of the aperiodic part corresponding to each prediction step size is obtained. Step 6: Algebraically add the total prediction results of the periodic part and the prediction results of the non-periodic part for each prediction step to obtain the complete ionospheric TEC prediction sequence.

2. The method for predicting the total electron content of the ionosphere according to claim 1, characterized in that, The standard time-frequency transformation in step 2 is based on the following formula: ； In the formula, The input time signal Standard time-frequency transform spectrum For real time, For real-time frequency, For the modulus of instantaneous frequency, For window width parameters, The input time signal, For time, The imaginary unit; Set window width parameters The scope will be expanded by substituting the extended ionospheric TEC data sequence. The standard time-frequency transform spectrum corresponding to the extended ionospheric TEC data sequence is obtained, and sub-period signals of different periods are determined in the standard time-frequency transform spectrum corresponding to the extended ionospheric TEC data sequence; the window width parameter is set. Then substitute the extended ionospheric TEC data sequence into The total standard time-frequency transform spectrum corresponding to each sub-period signal is obtained. Then, the standard time-frequency transform spectrum corresponding to the original ionospheric TEC data sequence is extracted from the total standard time-frequency transform spectrum corresponding to each sub-period signal, and the standard time-frequency transform spectrum corresponding to each sub-period signal is obtained respectively.

3. The method for predicting the total electron content of the ionosphere according to claim 2, characterized in that, For real signals, each sub-period signal is extracted from the corresponding standard time-frequency transform spectrum based on the following formula: ； In the formula, The extracted sub-period signal, for In frequency The standard time-frequency transform spectrum at that location, for In frequency The modulus of the standard time-frequency transform spectrum at that location. The angular frequency corresponding to each sub-period signal. Standard time-frequency transform spectrum The real part, It is the arctangent function. It is a cosine function.

4. The method for predicting the total electron content of the ionosphere according to claim 3, characterized in that, For complex signals, each sub-period signal is extracted from the corresponding standard time-frequency transform spectrum based on the following formula: ； In the formula, Standard time-frequency transform spectrum The imaginary part.

5. The method for predicting the total electron content of the ionosphere according to claim 4, characterized in that, The periodic prediction model is derived based on the following formula: ； In the formula, For prediction step size The overall prediction results for the periodic portion, To predict the start time, To predict the step size, For trend items, This is an operation for taking the real part of a complex number. This represents the sequence number of the sub-period signal. The number of sub-periodic signals. , , The first The amplitude, frequency, and phase of each sub-cycle signal; Take different prediction step sizes Substituting into the formula, we obtain the total prediction results for the periodic portion with different prediction step sizes.

6. The method for predicting the total electron content of the ionosphere according to claim 5, characterized in that, The ionospheric TEC data sequence is extended in the following way: using the ionospheric TEC data sequence as the original ionospheric TEC data sequence, the ionospheric TEC data sequence is copied twice, and then the two copies of the ionospheric TEC data sequence are added to both ends of the original ionospheric TEC data sequence to construct an extended ionospheric TEC data sequence.

7. The method for predicting the total electron content of the ionosphere according to claim 6, characterized in that, Step 5 specifically includes the following steps: Step 5.1: First, differencing the aperiodic reconstructed sequence, then performing an ADF test on the differencing aperiodic reconstructed sequence; if the p-value of the ADF test is less than a set value, the ADF test is passed; if the p-value of the ADF test is greater than or equal to the set value, the differencing aperiodic reconstructed sequence is differencing again, and the ADF test is performed again; until... The aperiodic reconstructed sequence after the second difference passes the ADF test and yields... The aperiodic reconstructed sequence after the second difference is an aperiodic reconstructed stationary sequence. , The degree of the difference; Step 5.2: Construct the autoregressive integral moving average prediction model. The autoregressive integral moving average prediction model is as follows: ； In the formula, To predict the step size is The prediction results for the non-periodic portion, for Aperiodic reconstructed stationary sequence at time 1 The value or predicted value, It is a constant term. It is the first One autoregressive coefficient, It is the autoregressive order. It is the first A moving average coefficient, It is the order of the moving average, the sequence number. Serial Number , yes Error term at time, for Error term at time; Determine the optimal Combination, for optimal Autoregressive order in the combination and moving average order and the corresponding constant term Autoregression coefficient and moving average coefficient The final autoregressive integral moving average prediction model was constructed. Take different prediction step sizes Substituting the results into the final autoregressive integral moving average prediction model, we obtain the prediction results for the non-periodic portion with different prediction step sizes.

8. The method for predicting the total electron content of the ionosphere according to claim 7, characterized in that, The autoregressive order and moving average order Determined in the following ways: Preset autoregression order and moving average order The range, traversing all Combination, for each The corresponding autoregressive integral moving average prediction models are constructed by combining these models. Then, the Akaike Information Criterion (AIC) or Bayesian Information Criterion (BIC) for each model is calculated. Finally, the model that minimizes either the AIC or BIC value is selected. Combination as optimal combination; Based on the optimal Autoregressive order of the combination and moving average order The optimal value is obtained using the maximum likelihood estimation method. The constant term in the corresponding autoregressive integral moving average prediction model Autoregression coefficient and moving average coefficient .

9. The method for predicting the total electron content of the ionosphere according to claim 8, characterized in that, Step 4 specifically includes the following steps: Step 4.1: Subtract all the sub-periodic signals extracted in Step 2 from the ionospheric TEC data sequence to obtain the aperiodic residual sequence. ,in, arrive These are the numbers 1 to 1 in the non-periodic residual sequence. Ionospheric TEC data for each hour. The length of the non-periodic residual sequence; Step 4.2: According to the selected window length Transforming the aperiodic residual sequence into a trajectory matrix : ,in, Let be the number of columns in the trajectory matrix. , arrive These are non-periodic residual sequences. medium length is The time delay segment, the window length Set it to the period of the largest amplitude. Step 4.3: Analyze the trajectory matrix. Perform singular value decomposition: ,in, For serial number, , For the first There are several elementary matrices, each corresponding to a component pattern in the aperiodic residual sequence. It is a matrix The Large eigenvalues and Trajectory matrices The left singular vector and the right singular vector of the matrix, for The transpose of the matrix; Step 4.4, The elementary matrices are selectively divided into Disjoint subsets Then, the elementary matrices in each subset are summed to obtain the contribution rate matrix for each subset: ,in, For the first A subset, number , For the first Contribution rate matrix corresponding to each subset; Then the trajectory matrix for The sum of the contribution rate matrices: ; Step 4.5, Each contribution rate matrix is ​​obtained by diagonal averaging. One reconstructed sequence; The reconstruction order r is determined based on the singular value contribution rate. The reconstruction sequences are sorted from highest to lowest singular value contribution rate. The first r reconstruction sequences are retained. Finally, the retained r reconstruction sequences are added together to obtain the aperiodic reconstruction sequence.

10. A computer device comprising a memory and a processor, wherein the memory stores a computer program, characterized in that, When the processor executes the computer program, it implements the steps of the method for predicting the total electron content of the ionosphere as described in any one of claims 1 to 9.

Images