A short-term ERP prediction method integrating SSA and cascaded LSTM

Abstract

The present invention discloses an ERP short-term prediction method that integrates SSA and cascaded LSTM, and relates to the field of satellite navigation and positioning technology. The method comprises the following steps: preprocessing ERP data; performing signal decomposition using a singular spectrum analysis method, performing spectrum analysis on the obtained signal components, and determining the periodic terms of the components; eliminating signal components with periodic terms less than n; reconstructing the remaining signal components into new ERP time series data after denoising and optimization; constructing a cascaded LSTM time series prediction model to calculate ERP short-term forecast data for the next m days; and adding correction terms to the ERP short-term forecast data to obtain a final ERP forecast sequence. The present invention achieves high-precision reconstruction of the ERP time series by denoising and optimizing the ERP time series signal. A cascaded LSTM prediction model is constructed in which multiple sub-models are interconnected and transmitted step by step, significantly suppressing the propagation of forecast errors and achieving high-precision and high-reliability forecasts of ERP data.

Description

Technical Field

[0001] The present invention relates to the fields of satellite navigation technology and radio navigation, and in particular to an ERP short-term prediction method integrating SSA and cascade LSTM. Background Art

[0002] Earth Rotation Parameters (ERP) are key parameters of the spacetime benchmark and play an indispensable role in deep space exploration, precise satellite orbit determination, geodynamic research, and high-precision timing services. These parameters primarily include polar motion (PM), length of day (LOD), and the difference between Universal Time (UT1) and Coordinated Universal Time (UTC). High-precision ERP data is crucial for orbit determination of deep space probes and high-precision GNSS positioning. Currently, high-precision ERP data can be obtained by combining space geodetic observation techniques such as very long baseline interferometry, global navigation satellite systems, and satellite laser ranging. However, such data exhibits a certain lag, making it difficult to meet the demands of real-time applications such as coordinate system transformation, precise satellite orbit determination, and deep space exploration. Therefore, the development of high-precision ERP prediction technology is urgently needed to provide precise and reliable ERP information for high-precision navigation, timing, and space missions. This is crucial for establishing an independent and controllable ERP service system.

[0003] Machine learning technology, with its excellent nonlinear processing capabilities, has been widely used in Earth rotation parameter prediction research, improving the forecast accuracy of existing technologies. Taking the polar motion X as an example, the MAE of the forecast results for the first, 10, and 30th days decreased from 1mas to 0.3mas, 4mas to 3.01mas, and 11mas to 6.97mas. The Long Short-Term Memory (LSTM) neural network, a classic recurrent neural network, has significant advantages in time series data prediction and has become a key method for ERP parameter prediction. A series of LSTM-based Earth rotation parameter prediction methods have been developed, including the LS-LSTM model, the CNN-LSTM model, the EWT-CNN-LSTM model, and the CNN-LSTM-Attention model.

[0004] However, the aforementioned forecasting method, based on a recursive LSTM neural network framework, ignores the inconsistency between training and application scenarios, leading to the continuous accumulation of forecast errors and severely impacting forecast accuracy. Specifically, based on the recursive LSTM neural network model, the forecast values ​​for the next 1 to m-1 days are used as model input to predict the data for the next m day. This inevitably results in prediction errors in the model input data for the predictions from days 2 to m, which is inconsistent with the high-precision input data used during training. This can easily lead to the accumulation of prediction errors and reduce the accuracy of subsequent forecasts. Furthermore, as a data-driven modeling approach, the accuracy of the training dataset is crucial to the effectiveness of model training. However, existing machine learning-based ERP forecasting methods directly utilize ERP time series data for training, ignoring the impact of signal noise and easily resulting in low model prediction accuracy. Summary of the Invention

[0005] The present invention provides a short-term ERP prediction method that integrates SSA (Singular Spectrum Analysis) and cascaded LSTM. This method can solve the problem of low ERP prediction accuracy caused by signal noise in ERP time series data and the cumulative prediction bias of the recursive LSTM model in the prior art.

[0006] The present invention provides an ERP short-term prediction method integrating SSA and cascaded LSTM, comprising:

[0007] Step 1: Preprocess the ERP time series data.

[0008] Step 2: Use singular spectrum analysis to decompose the ERP time series data and determine r signal components.

[0009] Step 3: Perform spectrum analysis on the signal components to determine the periodic terms of the signal components.

[0010] Step 4: According to the contribution rate of the signal components, the signal components with periodic terms less than E are eliminated to obtain the remaining signal components.

[0011] Step 5: Reconstruct the remaining signal components into new ERP time series data after denoising and optimization.

[0012] Step 6: Based on the new ERP time series data of the past t days, the constructed cascade LSTM time series forecast model is used to calculate the ERP short-term forecast time series data for the next m days, where the ERP short-term forecast time series data includes PM forecast data, UT1-UTC forecast data, and LOD forecast data.

[0013] Step 7: Add leap seconds and solid tides to the UT1-UTC forecast data, and add solid tides to the LOD forecast data to obtain the final ERP forecast sequence.

[0014] Preferably, in step 1, the ERP time series data includes PM time series data, UT1-UTC time series data and LOD time series data.

[0015] Eliminate leap seconds and the first solid tide in the UT1-UTC time series data , remove the second solid tide in the LOD time series data :

[0016] .

[0017] Where n1 represents the number of first solid tides or second solid tides in different periods, 、 、 、 is the amplitude of the periodic term, , n2 represents the number of solar and lunar nutation parameters, represents the sun-moon nutation parameter, Represents the sun and moon nutation parameters The integer multiplication factor of .

[0018] Preferably, the step 2 comprises:

[0019] Assign a time series window L, perform sliding lag arrangement on the ERP time series data, and obtain the trajectory matrix X:

[0020] .

[0021] in, Represents the one-dimensional time series corresponding to the ERP time series data, and N represents the length of the ERP time series data. Perform singular value decomposition on the trajectory matrix X:

[0022] .

[0023] Where U represents the left singular matrix, V represents the right singular matrix, and Σ represents the diagonal matrix composed of singular values.

[0024] According to the number of non-zero singular values, the rank r of the trajectory matrix X is determined, and the number of signal components is determined to be r.

[0025] Preferably, the step 3 comprises:

[0026] The signal components are processed using fast Fourier transform to obtain the power spectral density estimate corresponding to the signal components. The power spectrum of the signal components is analyzed to identify the frequency f_peak of the energy peak of the signal components and calculate the periodic term T_peak = 1 / f_peak corresponding to the signal components. The periodic term T_main corresponding to the peak with the highest energy in the power spectrum is selected as the periodic term of the corresponding signal component.

[0027] Preferably, the value range of E in step 4 is 30~60.

[0028] Preferably, the step 5 comprises:

[0029] Reconstruct the remaining signal components to obtain a new trajectory matrix , where q is the number of residual signal components;

[0030] Perform anti-diagonal averaging on the new trajectory matrix to obtain new ERP time series data :

[0031] .

[0032] in, represents the reconstructed component obtained by diagonal averaging, Represents the element of the trajectory matrix X at position (i, j), n=1,2,3,…,N, =max{L,K}, K represents the number of state vectors, and P=N represents the number of anti-diagonals.

[0033] Preferably, the step 6 comprises: constructing m sequential sub-models f1, f2, ..., f for the next m days respectively. m ; Input the ERP time series data of the past t days into the sequential sub-model f1, and use the sequential sub-model f1 to predict the ERP short-term forecast data of the next day 1 ; Use the ERP time series data of the past t-1 day and the ERP short-term forecast data of the next 1 day to input the sequential sub-model f2, and use the sequential sub-model f2 to predict the ERP data of the next 2 days Input the ERP time series data of the past t-2 days, the ERP short-term forecast data of the next 1 day, and the ERP short-term forecast data of the next 2 days into the third sequential sub-model f3, and use the third sequential sub-model f3 to predict the ERP data of the next 3 days. ; Input the ERP time series data of the past t-m+1 days and the ERP short-term forecast data of the future 1 to m-1 days into the mth sequential sub-model f m , using the mth sequential sub-model f m Predict and obtain the ERP short-term forecast data for the mth day , obtain the cascaded LSTM time series prediction model:

[0034] .

[0035] Where, Represents the ERP time series data for the first day in the past. Represents the ERP time series data for the past two days, Represents the ERP time series data for the past 3 days, represents the ERP time series data of the past t day, It represents the ERP short-term forecast data for the first day in the future. It represents the ERP short-term forecast data for the next two days. It represents the ERP short-term forecast data for the next 3 days. Represents the ERP short-term forecast data for the next m-1 day.

[0036] Preferably, in step 6, the value range of t is 20-60, and the value range of m is 10-30.

[0037] Preferably, the value range of L is 0~N / 2.

[0038] An embodiment of the present invention provides an ERP short-term forecasting method that integrates SSA and cascaded LSTM. Compared with the existing technology, the present invention has the following beneficial effects: the present invention uses singular spectrum analysis to decompose ERP time series data, eliminates the noise terms and some high-frequency change terms of the signal components to perform denoising optimization, and realizes high-precision reconstruction of ERP time series data; the present invention fully considers the scene feature changes of different forecast days in the future, realizes the error compensation transmission between the output of the preceding sequential sub-model and the input of the subsequent sequential sub-model through the cascade architecture, and constructs a cascaded LSTM time series forecasting model in which multiple sequential sub-models are interconnected and transmitted step by step, effectively suppresses the propagation of forecast errors, realizes high-precision forecasting of ERP time series data, and has extremely high effectiveness and reliability. BRIEF DESCRIPTION OF THE DRAWINGS

[0039] Figure 1 This is a flowchart provided for an embodiment of the present invention.

[0040] Figure 2 A flow chart of the cascaded LSTM time series forecasting model provided in an embodiment of the present invention.

[0041] Figure 3 A comparison chart of the prediction accuracy of the prediction method integrating SSA and cascaded LSTM provided in an embodiment of the present invention and other prediction methods in the polar motion X direction PMX and the polar motion Y direction PMY. DETAILED DESCRIPTION

[0042] To make the above-mentioned objects, features, and advantages of the present invention more readily apparent, specific embodiments of the present invention are described in detail below with reference to the accompanying drawings. The following description sets forth numerous specific details to facilitate a full understanding of the present invention. However, the present invention can be implemented in many other ways than those described herein, and those skilled in the art may make similar modifications without departing from the scope of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed below.

[0043] See also Figure 1~Figure 2 An embodiment of the present invention provides an ERP short-term forecasting method that integrates SSA and cascaded LSTM, including the following steps: Step 1: Collect EOP C04_20 sequence data provided by the International Earth Rotation and Reference System Service, and construct a dataset using 40 years of ERP time series data from January 1, 1984 to January 2, 2024. The dataset includes a training set and a test set. The training set is the ERP time series data from January 1, 1984 to January 2, 2020, and the test set is the ERP time series data from January 3, 2020 to January 2, 2024. The test set constructs 145 forecast periods, each of which forecasts the ERP short-term forecast time series data for the next 10 days.

[0044] Since leap seconds exist in UT1-UTC time series data, it is not conducive to building a forecast model. Therefore, leap seconds need to be removed. That is, when the difference between UT1-UTC for two consecutive days exceeds 0.9 seconds, all UT1-UTC time series data after the current time point will be subtracted by 1 second. In addition, the UT1-UTC time series data and LOD time series data need to be processed to remove solid tides, and the first solid tide in the UT1-UTC time series data should be removed. , remove the second solid tide in the LOD time series data , use formula (1) to calculate the first solid tide and the second solid tide:

[0045] Formula (1).

[0046] Where n1 represents the number of first solid tides or second solid tides in different periods, 、 、 、 is the amplitude of the periodic term, , n2 represents the number of solar and lunar nutation parameters, represents the sun-moon nutation parameter, Represents the sun and moon nutation parameters The integer multiplication factor of .

[0047] Step 2: Use the singular spectrum analysis method to decompose the ERP time series data. First, set the suitable time series window L value to 2190, and perform sliding lag arrangement on the ERP time series data to obtain the trajectory matrix, and obtain the trajectory matrix X shown in formula (2):

[0048] Formula (2).

[0049] in, Represents the one-dimensional time series corresponding to the ERP time series data, and N represents the length of the ERP time series data.

[0050] Secondly, perform singular value decomposition on the trajectory matrix X, as shown in formula (3):

[0051] Formula (3),

[0052] Where U represents the left singular matrix, V represents the right singular matrix, and Σ represents the diagonal matrix composed of singular values.

[0053] Finally, the rank r of the trajectory matrix X, that is, the number of signal components, is determined according to the number of non-zero singular values.

[0054] Step 3: Perform spectral analysis on the signal components. First, use fast Fourier transform to process the signal components and obtain the power spectral density estimation corresponding to each signal component. Secondly, analyze the power spectrum of each signal component one by one, identify the frequency f_peak of its energy peak, and calculate its corresponding periodic term T_peak=1 / f_peak. Finally, for signal components with multiple peaks, select the periodic term T_main corresponding to the peak with the highest energy in the power spectrum as the periodic term of the signal component.

[0055] Step 4: According to the contribution rate of the signal components, the signal components with periodic terms less than 30 are eliminated.

[0056] Step 5: Reconstruct the remaining 154 signal components to obtain a new trajectory matrix , use formula (4) to perform anti-diagonal mean on the matrix to obtain the new ERP time series data :

[0057] Formula (4),

[0058] Where K represents the number of state vectors, which is N-L+1, and P=N represents the number of anti-diagonals.

[0059] Step 6: Build Figure 2The cascaded LSTM time series forecasting model shown in the figure includes: constructing 10 sequential sub-models for the next 10 days respectively, inputting the ERP time series data of the past 10 days into the sequential sub-model f1, and the sequential sub-model f1 predicts the ERP short-term forecast data for the first day in the future ; Construct the second sequential sub-model f2, use the ERP data of the past 9 days and the ERP forecast value of the next day to form the input data of the sequential sub-model f2, and forecast the short-term ERP forecast data of the next day Similarly, the input data of the third sequential sub-model f3 includes the ERP time series data of the past 8 days, the ERP short-term forecast time series data of the first day in the future, the ERP short-term forecast time series data of the second day, and the ERP short-term forecast data of the third day in the future. ; And so on, until the 10th sequential sub-model f m Forecast the ERP short-term forecast data for the 10th day , the cascade LSTM time series prediction model is shown in formula (5):

[0060] Formula (5).

[0061] Step 7: ERP short-term forecast data includes PM forecast data, UT1-UTC forecast data and LOD forecast data. Leap seconds and solid tides are added to the UT1-UTC forecast data in the ERP short-term forecast time series data, and solid tides are added to the LOD forecast data in the ERP short-term forecast time series data to obtain the final ERP forecast sequence.

[0062] like Figure 3 As shown, in order to verify the feasibility of the method of the present invention, taking the PM time series data in the ERP time series data as an example, we constructed Figure 3 The polar motion in the X direction PMX and Figure 3 The prediction method of the polar motion Y direction PMY fusion SSA and cascade LSTM in the middle right figure is compared with the measured data of the existing CNN-LSTM-Attention method, LS+AR+Kalman method, Bulletin A, and the prediction method of the fusion SSA and cascade LSTM provided by the present invention. The mean absolute error (MAE) of the forecast results for each day in the future is calculated to evaluate the forecast accuracy. The results are shown in Tables 1, 2 and Figure 3 shown.

[0063] Table 1. Comparison of the prediction accuracy of the four prediction methods in the PMX direction (unit: mas)

[0064] Forecast days Prediction method integrating SSA and cascaded LSTM CNN-LSTM-Attention method LS+AR+Kalman method Bulletin A 1 0.88 0.57 0.24 0.26 2 1.00 0.84 0.61 0.56 3 1.07 1.19 1.04 0.90 4 1.21 1.54 1.49 1.30 5 1.43 1.88 1.93 1.70 6 1.64 2.22 2.34 2.05 7 1.92 2.54 2.70 2.31 8 2.24 2.83 3.09 2.57 9 2.48 3.12 3.45 2.82 10 2.93 3.40 3.78 3.07 average value 1.68 2.01 2.07 1.75

[0065] Table 2. Comparison of the prediction accuracy of the four prediction methods in the PMY direction (unit: mas)

[0066] Forecast days Prediction method integrating SSA and cascaded LSTM CNN-LSTM-Attention method LS+AR+Kalman method Bulletin A 1 0.57 0.43 0.21 0.21 2 0.59 0.66 0.48 0.40 3 0.63 0.89 0.74 0.61 4 0.70 1.10 1.01 0.83 5 0.78 1.32 1.27 1.06 6 0.89 1.50 1.53 1.25 7 1.01 1.66 1.78 1.43 8 1.21 1.82 2.01 1.60 9 1.41 1.98 2.27 1.75 10 1.58 2.15 2.50 1.89 average value 0.94 1.35 1.38 1.10

[0067] From Table 1, Table 2 and Figure 3 It can be seen that in the X-direction of polar motion, the prediction accuracy of the SSA and cascaded LSTM fusion method is slightly lower in the first two days. However, in the prediction results for days 4-10, the MAE of the SSA and cascaded LSTM fusion method is significantly lower than that of the other two models and Bulletin A, with the corresponding MAE on the 10th day being 2.93mas. Similarly, a similar trend is observed in the Y-direction of polar motion. The prediction results of the LS+AR+Kalman method and the CNN-LSTM-Attention method are basically the same, and the MAE of the forecast results for the 10th day exceeds 2mas. The prediction results of Bulletin A are slightly better than the above two methods. The prediction method of the SSA and cascaded LSTM fusion method shows good prediction results for days 4-10, with the MAE on the 10th day being only 1.58mas. These comparisons effectively demonstrate the high-precision prediction advantage of the SSA and cascaded LSTM fusion prediction method provided by the present invention.

[0068] The above-described embodiments merely illustrate several implementations of the present invention, and while their descriptions are relatively specific and detailed, they should not be construed as limiting the scope of the patent. It should be noted that a person skilled in the art would be able to make numerous variations and improvements without departing from the spirit of the present invention, all of which fall within the scope of protection of the present invention. Therefore, the scope of protection of the patent for this invention shall be determined by the appended claims.

Claims

1. A short-term ERP prediction method integrating SSA and cascaded LSTM, characterized in that: include: Step 1: Preprocess the ERP time series data; Step 2: Decompose the ERP time series data using singular spectrum analysis to determine r signal components; Step 3, performing spectrum analysis on the signal components to determine the periodic terms of the signal components; Step 4: According to the contribution rate of the signal components, the signal components with periodic terms less than E are eliminated to obtain the remaining signal components; Step 5: Reconstruct the remaining signal components into new ERP time series data after denoising and optimization; Step 6: Based on the new ERP time series data of the past t days, the constructed cascade LSTM time series forecast model is used to calculate the ERP short-term forecast time series data for the next m days. The ERP short-term forecast time series data includes PM forecast data, UT1-UTC forecast data, and LOD forecast data. Step 7: Add leap seconds and solid tides to the UT1-UTC forecast data, and add solid tides to the LOD forecast data to obtain the final ERP forecast sequence; In step 1, the ERP time series data includes PM time series data, UT1-UTC time series data and LOD time series data; Eliminate leap seconds and the first solid tide in the UT1-UTC time series data , remove the second solid tide in the LOD time series data : ， Where n1 represents the number of first solid tides or second solid tides in different periods, 、 、 、 is the amplitude of the periodic term, , n2 represents the number of solar and lunar nutation parameters, represents the sun-moon nutation parameter, Represents the sun and moon nutation parameters The integer multiplication factors of ; The step 6 includes: constructing m sequential sub-models f1, f2, ..., f for the next m days respectively. m ; Input the ERP time series data of the past t days into the sequential sub-model f1, and use the sequential sub-model f1 to predict the ERP short-term forecast data of the next day 1 ; Use the ERP time series data of the past t-1 day and the ERP short-term forecast data of the next 1 day to input the sequential sub-model f2, and use the sequential sub-model f2 to predict the ERP data of the next 2 days Input the ERP time series data of the past t-2 days, the ERP short-term forecast data of the next 1 day, and the ERP short-term forecast data of the next 2 days into the third sequential sub-model f3, and use the third sequential sub-model f3 to predict the ERP data of the next 3 days. ; Input the ERP time series data of the past t-m+1 days and the ERP short-term forecast data of the future 1 to m-1 days into the mth sequential sub-model f m , using the mth sequential sub-model f m Predict and obtain the ERP short-term forecast data for the mth day , obtain the cascaded LSTM time series prediction model: , Where, Represents the ERP time series data for the first day in the past. Represents the ERP time series data for the past two days, Represents the ERP time series data for the past 3 days, represents the ERP time series data of the past t day, It represents the ERP short-term forecast data for the first day in the future. It represents the ERP short-term forecast data for the next two days. It represents the ERP short-term forecast data for the next 3 days. Represents the ERP short-term forecast data for the next m-1 day.

2. The ERP short-term prediction method integrating SSA and cascaded LSTM according to claim 1 is characterized in that: The step 2 comprises: Assign a time series window L, perform sliding lag arrangement on the ERP time series data, and obtain the trajectory matrix X: ， in, represents the one-dimensional time series corresponding to the ERP time series data, and N represents the length of the ERP time series data; Perform singular value decomposition on the trajectory matrix X: ， Where U represents the left singular matrix, V represents the right singular matrix, and Σ represents the diagonal matrix composed of singular values; According to the number of non-zero singular values, the rank r of the trajectory matrix X is determined, and the number of signal components is determined to be r.

3. The ERP short-term prediction method integrating SSA and cascaded LSTM according to claim 1 is characterized in that: The step 3 comprises: The signal components are processed using fast Fourier transform to obtain the power spectral density estimate corresponding to the signal components. The power spectrum of the signal components is analyzed to identify the frequency f_peak of the energy peak of the signal components and calculate the periodic term T_peak = 1 / f_peak corresponding to the signal components. The periodic term T_main corresponding to the peak with the highest energy in the power spectrum is selected as the periodic term of the corresponding signal component.

4. The ERP short-term prediction method integrating SSA and cascaded LSTM according to claim 1 is characterized in that: The value range of E in step 4 is 30~60.

5. The ERP short-term prediction method integrating SSA and cascaded LSTM according to claim 1 is characterized in that: The step 5 comprises: Reconstruct the remaining signal components to obtain a new trajectory matrix , where q is the number of residual signal components; Perform anti-diagonal averaging on the new trajectory matrix to obtain new ERP time series data : ， in, represents the reconstructed component obtained by diagonal averaging, Represents the element of the trajectory matrix X at position (i, j), n=1,2,3,…,N, =max{L,K}, K represents the number of state vectors, and P=N represents the number of anti-diagonals.

6. The ERP short-term prediction method integrating SSA and cascaded LSTM according to claim 5 is characterized in that: In step 6, the value range of t is 20-60, and the value range of m is 10-30.

7. The ERP short-term prediction method integrating SSA and cascaded LSTM according to claim 2 is characterized in that: The value range of L is 0~N / 2.

Images