A high-precision space target orbit prediction method based on numerical fitting

Abstract

The application discloses a high-precision space target orbit prediction method based on numerical fitting, and relates to the technical field of space target orbit prediction.The frequency domain decomposition of the application is combined with a harmonic compensation model to completely eliminate seasonal system drift caused by the gravitation of the sun and the moon; high-frequency components are stripped of attitude noise through light pressure parameterization and smooth difference to avoid the transmission of perturbation coupling error; a harmonic order adaptive mechanism captures the main disturbance term with the least amount of calculation; least square batch estimation realizes one-time solution of coefficients; template library preloading technology effectively compresses the time consumption of new target initialization; in addition, the application introduces a geomagnetic modulation factor to dynamically suppress compensation distortion during magnetic storms, and a sliding window mechanism guarantees real-time updating of coefficients, and stable prediction is still maintained under extreme events; photometric curve inversion is adopted instead of attitude measurement to make non-cooperative target prediction possible; tangential component directional compensation ensures the stability of the orbit plane and reduces satellite position fuel consumption.

Description

Technical Field

[0001] This invention relates to the field of space target orbit prediction technology, and in particular to a high-precision space target orbit prediction method based on numerical fitting. Background Technology

[0002] In key fields such as communications and meteorological monitoring, geostationary orbit satellites need to maintain millimeter-level positional stability. Current high-precision forecasting schemes have integrated numerical integration and perturbation compensation models, such as empirical parameters of the Earth's gravity field and solar radiation pressure in JGM-3, but they still show significant defects in scenarios involving gravitational perturbations by the Sun and Moon. Especially during the spring / autumn equinox, when the Sun, Earth, and satellite form a near-straight configuration, the gravitational forces of the Sun and Moon generate a continuous torque on the orbital tangential. Traditional spherical harmonic function models accumulate tangential velocity errors because they do not separate the perturbation coupling effect.

[0003] In recent years, mainstream solutions have included enhancing the numerical integration order, such as using the RKF8(9) integrator. Although this improves accuracy, it also increases the computational load, making it difficult to meet the real-time early warning requirements of constellations. Alternatively, data-driven compensation such as LSTM can be introduced, but this relies on massive amounts of historical ephemeris data and is slow to respond to newly built satellites and sudden orbital disturbances. Neither of these solutions has solved the problem of frequency domain aliasing caused by perturbation. The gravitational forces of the Sun and Moon generate harmonic interference with the non-spherical gravity of the Earth in the 0.1-0.3Hz frequency band, and traditional time-domain algorithms cannot effectively decouple them. A test conducted by a European space agency in 2023 showed that in a dual-satellite co-position control scenario, such interference resulted in a collision probability misjudgment rate as high as 1 / 1000. Therefore, a high-precision space target orbit prediction method based on numerical fitting is urgently needed to solve this problem. Summary of the Invention

[0004] In view of the aforementioned existing problems, the present invention is proposed.

[0005] This invention provides a high-precision space target orbit prediction method based on numerical fitting to solve the problem of tangential drift caused by the coplanar gravitational interference of the Sun and Moon in geostationary orbits. Existing models have difficulty balancing accuracy and real-time performance due to perturbation coupling and computational load.

[0006] To solve the above-mentioned technical problems, the present invention provides the following technical solution:

[0007] This invention provides a high-precision space target orbit prediction method based on numerical fitting, which includes:

[0008] Step S1: Obtain the initial orbital parameters and real-time space environment data of the space target;

[0009] Step S2: Perform frequency domain decomposition on the perturbation acceleration to separate the original acceleration observation signal into low-frequency components, mid-frequency components, and high-frequency components;

[0010] Step S3: Perform cubic spline fitting on the intermediate frequency component to generate a perturbation compensation curve;

[0011] Step S4: Combine the low-frequency component, the perturbation compensation curve, and the high-frequency component to reconstruct the acceleration;

[0012] Step S5: Perform orbital numerical integration based on the reconstructed acceleration, and output the orbital state vector of the space target at the preset time point.

[0013] As a preferred embodiment of the high-precision space target orbit prediction method based on numerical fitting described in this invention, the frequency domain decomposition in step S2 specifically includes:

[0014] The wavelet packet transform algorithm is used to perform multi-scale decomposition of the acceleration observation signal with predefined basis functions;

[0015] The low-frequency component corresponds to the frequency band dominated by the Earth's non-spherical gravity, the mid-frequency component corresponds to the frequency band dominated by the Sun and Moon's gravity, and the high-frequency component corresponds to the frequency band dominated by solar radiation pressure.

[0016] As a preferred embodiment of the high-precision space target orbit prediction method based on numerical fitting described in this invention, wherein: in step S3, when performing cubic spline fitting on the mid-frequency components:

[0017] The distribution of fitted nodes is dynamically adjusted based on the instantaneous rate of change of the right ascension of the ascending node of the orbit.

[0018] Increase node density during periods when the rate of change of right ascension at the ascending node exceeds a preset threshold.

[0019] As a preferred embodiment of the high-precision space target orbit prediction method based on numerical fitting described in this invention, it further includes a periodic compensation term generation step:

[0020] Construct a compensation acceleration term based on orbital period characteristics;

[0021] The compensation term is superimposed on the tangential component of the reconstructed acceleration;

[0022] The construction of the compensation acceleration term includes:

[0023] Obtain the periodic characteristic parameters of the target orbit;

[0024] Construct a combined expression that includes sine and cosine functions;

[0025] The amplitude coefficients in the expression were determined by fitting historical ephemeris data.

[0026] As a preferred embodiment of the high-precision space target orbit prediction method based on numerical fitting described in this invention, the periodic compensation term generation step includes:

[0027] The formula for calculating the fundamental frequency angular velocity is:

[0028]

[0029] Where, ω p The fundamental frequency angular velocity of the target orbit, expressed in rad / s. -1 T p The target orbital period, in seconds;

[0030] Determine the harmonic order adaptively:

[0031]

[0032] Where, N h This represents the highest order of harmonics, and its value is a positive integer. The up rounding operator is represented, κ1 represents the eccentricity weighting coefficient (dimensionless), e represents the orbital eccentricity (dimensionless), κ2 represents the altitude weighting coefficient (dimensionless), and H represents the average orbital altitude (in km).

[0033] Construct a harmonic compensation model, expressed as:

[0034]

[0035] Among them, a c (t) represents the value of the periodic compensation acceleration at time t, in milliseconds (ms). -2 k is the harmonic order index, which is dimensionless, and A k The amplitude coefficient of the k-th harmonic sinusoidal wave is given in milliseconds (ms). -2 sin() is the sine function, B k The amplitude coefficient of the k-th harmonic cosine is given in milliseconds (ms). -2 cos() is the cosine function, t represents time in seconds;

[0036] Batch estimation of amplitude coefficients is performed, expressed as:

[0037]

[0038] Where C represents the magnitude coefficient vector, according to Arrangement, unit: milliseconds -2 Φ represents the residual observation matrix, with dimensions M×2N. h Dimensionless This represents the matrix transpose operator, where Δa represents the column vector of historical acceleration residuals, with dimension M×1 and unit ms. -2 M is the number of residual samples, which is dimensionless;

[0039] In the formula, the elements of the observation matrix are defined as:

[0040] Φ i,2k-1 =sin(kω) p t i ),Φ i,2k =cos(kω) p t i ),

[0041] Where, Φ i,j Let t be the element in the i-th row and j-th column of matrix Φ, dimensionless. i Let be the time label of the i-th residual sample, in seconds;

[0042] Modulated geomagnetic disturbance:

[0043] γ=1-αΔK p ,A′ k =γA k ,B′ k =γB k ,

[0044] Where γ is the geomagnetic modulation factor, dimensionless, α is the geomagnetic sensitivity coefficient, dimensionless, and ΔK p The latest instantaneous change in the geomagnetic index is dimensionless, A′. k ,B′ k This represents the harmonic amplitude coefficient after geomagnetic modulation, in milliseconds (ms). -2 ;

[0045] The formula for superimposing tangential component compensation is as follows:

[0046]

[0047] in, This represents the tangential acceleration vector after compensation. a represents the original tangential acceleration vector. c (t) represents the periodic compensation acceleration scalar obtained in step S3, in milliseconds (ms). -2 .

[0048] As a preferred embodiment of the high-precision space target orbit prediction method based on numerical fitting described in this invention, the coefficient update of the compensation acceleration term includes:

[0049] When the change in the real-time spatial environment index exceeds the threshold value...

[0050] The coefficients were refitted using the track residual data within the sliding time window;

[0051] The width setting of the sliding time window includes:

[0052] The window width is dynamically adjusted based on the magnitude of sudden changes in the geomagnetic index.

[0053] The magnitude of the mutation is negatively correlated with the window width.

[0054] As a preferred embodiment of the high-precision space target orbit prediction method based on numerical fitting described in this invention, wherein: in step S5, the initialization of the orbit numerical integration includes:

[0055] Match a pre-stored template library based on the target orbit inclination and eccentricity;

[0056] The perturbation fitting parameters of the loaded matching template are used as initial values ​​for numerical integration.

[0057] As a preferred embodiment of the high-precision space target orbit prediction method based on numerical fitting described in this invention, the processing of the high-frequency components includes:

[0058] Establish a parameterized model of the equivalent area of ​​solar radiation pressure;

[0059] Based on the inversion model parameters of historical luminosity curves;

[0060] The establishment of the parameterized model includes:

[0061] Simplify the space target into an ellipsoidal geometry;

[0062] Establish the trigonometric relationship between the illuminated area and the solar incidence angle.

[0063] As a preferred embodiment of the high-precision space target orbit prediction method based on numerical fitting described in this invention, the high-frequency component processing steps include:

[0064] The formula for calculating the angle of incidence of sunlight is:

[0065]

[0066] Where θ(t) represents the angle of incidence of sunlight, in rad. This represents the unit vector along the principal inertial axis of the target body; it is dimensionless. This represents a unit vector representing the direction of sunlight incidence, which varies with time t and is dimensionless.

[0067] The equivalent illuminated area is parameterized and expressed as:

[0068] A eff (θ)=A0[1-λ1sin 2 θ+λ2sin 4 θ],

[0069] Among them, A eff (θ) represents the equivalent area exposed to solar radiation, in meters (m²). 2A0 represents the normal incident cross-sectional area, in meters. 2 λ1 represents the quadratic modulation coefficient, with a value range of 0-1, sin() represents the sine function, and λ2 is the quartic modulation coefficient, with a value range of 0-1.

[0070] The formula for calculating solar radiation pressure acceleration is:

[0071]

[0072] Among them, a sp (t) represents the solar radiation pressure acceleration vector, in milliseconds (ms). -2 C R The reflection and absorption coefficient is dimensionless, and P0 represents the solar radiation pressure at 1 AU, a constant of 4.56 × 10⁻⁶. -6 N m -2 m t The target mass in space is expressed in kg.

[0073] Parameter inversion based on photometric curves is expressed as follows:

[0074]

[0075] in, Let Ψ be the column vector of the three parameters to be estimated, and Ψ represent the photometric observation matrix. The element in the i-th row is:

[0076] Ψ i1 =1,Ψ i2 =-sin 2 θ(t i ),Ψ i3 =sin 4 θ(t i (), dimensionless, F represents the historical photometric conversion flux column vector, unit W / m -2 , t i This represents the time of the i-th observation epoch, in seconds.

[0077] Extracting high-frequency acceleration:

[0078]

[0079] Among them, a hf (t) represents the high-frequency component acceleration vector, in milliseconds (ms). -2 , Indicates that for a sp (t) The result after applying cubic spline smoothing with one orbital period, in milliseconds. -2 .

[0080] As a preferred embodiment of the high-precision space target orbit prediction method based on numerical fitting described in this invention, when outputting the orbit state vector:

[0081] The position and velocity vectors are in a geocentric inertial coordinate system.

[0082] The timestamps are in UTC format and include millisecond-level timestamps.

[0083] The beneficial effects of this invention are as follows: the frequency domain decomposition combined with the harmonic compensation model completely eliminates the seasonal system drift caused by the gravitational pull of the sun and moon; the high-frequency components are stripped of attitude noise through optical pressure parameterization and smoothing differential, avoiding the propagation of perturbation coupling errors; the harmonic order adaptive mechanism captures the main perturbation term with minimal computation; least squares batch estimation achieves one-time solution of coefficients; the template library preloading technology effectively compresses the initialization time of new targets; in addition, this invention introduces a geomagnetic modulation factor to dynamically suppress compensation distortion during geomagnetic storms, and the sliding window mechanism ensures real-time updates of coefficients, maintaining stable forecasts even under extreme events; photometric curve inversion is used to replace attitude measurement, making non-cooperative target forecasting possible; tangential component directional compensation ensures orbital plane stability and reduces satellite position maintenance fuel consumption. Attached Figure Description

[0084] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the following description of the embodiments will be briefly introduced. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0085] Figure 1 This is a flowchart illustrating the high-precision space target orbit prediction method based on numerical fitting in Example 1. Detailed Implementation

[0086] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings.

[0087] Many specific details are set forth in the following description in order to provide a full understanding of the invention. However, the invention may also be practiced in other ways different from those described herein, and those skilled in the art can make similar extensions without departing from the spirit of the invention. Therefore, the invention is not limited to the specific embodiments disclosed below.

[0088] Secondly, the term "one embodiment" or "embodiment" as used herein refers to a specific feature, structure, or characteristic that may be included in at least one implementation of the present invention. The phrase "in one embodiment" appearing in different places in this specification does not necessarily refer to the same embodiment, nor is it a single or selective embodiment that is mutually exclusive with other embodiments.

[0089] Example 1, referring to Figure 1 This embodiment provides a high-precision space target orbit prediction method based on numerical fitting, including:

[0090] Step S1: Obtain the initial orbital parameters and real-time space environment data of the space target;

[0091] Step S2: Perform frequency domain decomposition on the perturbation acceleration to separate the original acceleration observation signal into low-frequency components, mid-frequency components, and high-frequency components;

[0092] In step S2, frequency domain decomposition specifically includes:

[0093] The wavelet packet transform algorithm is used to perform multi-scale decomposition of the acceleration observation signal with predefined basis functions;

[0094] Among them, the low-frequency component corresponds to the frequency band dominated by the Earth's non-spherical gravity, the mid-frequency component corresponds to the frequency band dominated by the Sun and Moon's gravity, and the high-frequency component corresponds to the frequency band dominated by solar radiation pressure.

[0095] In step S3, when performing cubic spline fitting on the mid-frequency components:

[0096] The distribution of fitted nodes is dynamically adjusted based on the instantaneous rate of change of the right ascension of the ascending node of the orbit.

[0097] Increase node density during periods when the rate of change of right ascension at the ascending node exceeds a preset threshold;

[0098] Step S3: Perform cubic spline fitting on the mid-frequency components to generate the perturbation compensation curve;

[0099] Step S4: Combine the low-frequency component, the perturbation compensation curve, and the high-frequency component to reconstruct the acceleration;

[0100] Step S5: Perform orbital numerical integration based on the reconstructed acceleration, and output the orbital state vector of the space target at the preset time point;

[0101] In step S5, the initialization of the orbital numerical integration includes:

[0102] Match a pre-stored template library based on the target orbit inclination and eccentricity;

[0103] The perturbation fitting parameters of the matching template are loaded as initial values ​​for numerical integration;

[0104] Processing of high-frequency components includes:

[0105] Establish a parameterized model of the equivalent area of ​​solar radiation pressure;

[0106] Based on the inversion model parameters of historical luminosity curves;

[0107] The establishment of a parametric model includes:

[0108] Simplify the space target into an ellipsoidal geometry;

[0109] Establish the trigonometric relationship between the illuminated area and the solar incidence angle;

[0110] The processing steps for high-frequency components include:

[0111] The formula for calculating the angle of incidence of sunlight is:

[0112]

[0113] Where θ(t) represents the angle of incidence of sunlight, in rad. This represents the unit vector along the principal inertial axis of the target body; it is dimensionless. This represents a unit vector representing the direction of sunlight incidence, which varies with time t and is dimensionless.

[0114] The equivalent illuminated area is parameterized and expressed as:

[0115] A eff (θ)=A0[1-λ1sin 2 θ+λ2sin 4 θ],

[0116] Among them, A eff (θ) represents the equivalent area exposed to solar radiation, in meters (m²). 2 A0 represents the normal incident cross-sectional area, in meters. 2 λ1 represents the quadratic modulation coefficient, with a value range of 0-1, sin() represents the sine function, and λ2 is the quartic modulation coefficient, with a value range of 0-1.

[0117] The formula for calculating solar radiation pressure acceleration is:

[0118]

[0119] Among them, a sp (t) represents the solar radiation pressure acceleration vector, in milliseconds (ms). -2 C R The reflection and absorption coefficient is 1.0–1.8, dimensionless, and P0 represents the solar radiation pressure at 1 AU, a constant of 4.56 × 10⁻⁶. -6 N m -2 m t The target mass in space is expressed in kg.

[0120] Parameter inversion based on photometric curves is expressed as follows:

[0121]

[0122] in, Let Ψ be the column vector of the three parameters to be estimated, and Ψ represent the photometric observation matrix. The element in the i-th row is:

[0123] Ψ i1 =1,Ψ i2 =-sin 2 θ(t i ),Ψ i3 =sin 4 θ(t i (), dimensionless, F represents the historical photometric conversion flux column vector, unit W / m -2 , t i This represents the time of the i-th observation epoch, in seconds.

[0124] Extracting high-frequency acceleration:

[0125]

[0126] Among them, a hf (t) represents the high-frequency component acceleration vector, in milliseconds (ms). -2 , Indicates that for a sp (t) The result after applying cubic spline smoothing with one orbital period, in milliseconds. -2 ;

[0127] Specifically, by defining the incident angle and performing a fourth-order sine expansion on the illuminated area, the rapid area oscillations caused by the geometric coupling of the target attitude can be explicitly captured. These oscillations directly drive the high-frequency fluctuations of solar radiation pressure. The three parameters are obtained by least-squares inversion of the photometric curve, making the model match the actual reflection characteristics without the need for attitude priors. Then, the instantaneous radiation pressure is differiated from its smoothed value to obtain a pure high-frequency component, thereby extracting the disturbances caused by rapid attitude jitter or structural elastic vibrations. This component can be injected at high frequency in subsequent numerical integration without changing the low-frequency energy distribution, effectively suppressing the periodic sub-periodic oscillations of the orbital state vector, and improving the robustness of short-term forecasts and long-term convergence.

[0128] This method outputs the orbital state vector:

[0129] The position and velocity vectors are in a geocentric inertial coordinate system.

[0130] The timestamps are in UTC format and include millisecond-level timestamps;

[0131] The method also includes a periodic compensation term generation step:

[0132] Construct a compensation acceleration term based on orbital period characteristics;

[0133] The compensation term is superimposed on the tangential component of the reconstructed acceleration;

[0134] The construction of the compensation acceleration term includes:

[0135] Obtain the periodic characteristic parameters of the target orbit;

[0136] Construct a combined expression that includes sine and cosine functions;

[0137] The amplitude coefficients in the expression were determined by fitting historical ephemeris data;

[0138] The steps for generating periodic compensation terms include:

[0139] The formula for calculating the fundamental frequency angular velocity is:

[0140]

[0141] Where, ω p The fundamental frequency angular velocity of the target orbit, expressed in rad / s. -1 T p The target orbital period, in seconds, is obtained through TLE data parsing.

[0142] Determine the harmonic order adaptively:

[0143]

[0144] Where, N h This represents the highest order of harmonics, and its value is a positive integer. The up-rounding operator is indicated; κ1 represents the eccentricity weighting coefficient, with an empirical range of 8-12 and is dimensionless; e represents the orbital eccentricity, dimensionless; κ2 represents the altitude weighting coefficient, with an empirical range of 1-3 and is dimensionless; and H represents the average orbital altitude, in km.

[0145] Construct a harmonic compensation model, expressed as:

[0146]

[0147] Among them, a c (t) represents the value of the periodic compensation acceleration at time t, in milliseconds (ms). -2 k is the harmonic order index, which is dimensionless, and A k The amplitude coefficient of the k-th harmonic sinusoidal wave is given in milliseconds (ms). -2 sin() is the sine function, B k The amplitude coefficient of the k-th harmonic cosine is given in milliseconds (ms). -2 cos() is the cosine function, t represents time in seconds;

[0148] Batch estimation of amplitude coefficients is performed, expressed as:

[0149]

[0150] Where C represents the magnitude coefficient vector, according to Arrangement, unit: milliseconds -2 Φ represents the residual observation matrix, with dimensions M×2N. h Dimensionless This represents the matrix transpose operator, where Δa represents the column vector of historical acceleration residuals, with dimension M×1 and unit ms. -2 M is the number of residual samples, which is dimensionless;

[0151] In the formula, the elements of the observation matrix are defined as:

[0152] Φ i,2k-1 =sin(kω) p t i ),Φ i,2k =cos(kω) p t i ),

[0153] Where, Φ i,j Let t be the element in the i-th row and j-th column of matrix Φ, dimensionless. i Let be the time label of the i-th residual sample, in seconds;

[0154] Modulated geomagnetic disturbance:

[0155] γ=1-αΔK p A k =γA′ k ,B′ k =γB k ,

[0156] Where γ is the geomagnetic modulation factor, dimensionless, α is the geomagnetic sensitivity coefficient, calibrated in the range of 0-0.05, dimensionless, and ΔK p The latest instantaneous change in the geomagnetic index is dimensionless, A′. k ,B′ k This represents the harmonic amplitude coefficient after geomagnetic modulation, in milliseconds (ms). -2 ;

[0157] The formula for superimposing tangential component compensation is as follows:

[0158]

[0159] in, This represents the tangential acceleration vector after compensation. a represents the original tangential acceleration vector. c (t) represents the periodic compensation acceleration scalar obtained in step S3, in milliseconds (ms). -2 ;

[0160] Specifically, the fundamental frequency angular velocity is derived from the orbital period, and the harmonic order is adaptively determined based on eccentricity and altitude, thus conforming to the dynamic characteristics of different orbital types in terms of algorithm structure. The harmonic compensation model transforms the long-term residual into a set of sine and cosine basis functions. The least squares batch estimation completes the solution of all amplitudes through a single matrix inversion, and the computational complexity increases linearly with the sample size, meeting the requirements of real-time processing. The geomagnetic disturbance modulation stage introduces an external space environment index, which can dynamically compress the harmonic amplitude during extreme geomagnetic storms and suppress the risk of miscompensation. Finally, the compensation amount is limited to tangential direction superposition, which effectively improves orbital energy error without destroying orbital plane stability.

[0161] The coefficient update for the compensation acceleration term includes:

[0162] When the change in the real-time spatial environment index exceeds the threshold value...

[0163] The coefficients were refitted using the track residual data within the sliding time window;

[0164] The width settings for the sliding time window include:

[0165] The window width is dynamically adjusted based on the magnitude of sudden changes in the geomagnetic index.

[0166] The magnitude of the mutation is negatively correlated with the window width.

[0167] It should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention, and all such modifications or substitutions should be covered within the scope of the claims of the present invention.

Claims

1. A high-precision space target orbit prediction method based on numerical fitting, characterized in that, include, Step S1: Obtain the initial orbital parameters and real-time space environment data of the space target; Step S2: Perform frequency domain decomposition on the perturbation acceleration to separate the original acceleration observation signal into low-frequency components, mid-frequency components, and high-frequency components; Step S3: Perform cubic spline fitting on the intermediate frequency component to generate a perturbation compensation curve; Step S4: Combine the low-frequency component, the perturbation compensation curve, and the high-frequency component to reconstruct the acceleration; Step S5: Perform orbital numerical integration based on the reconstructed acceleration, and output the orbital state vector of the space target at the preset time point; It also includes the step of generating periodic compensation terms: Construct a compensation acceleration term based on orbital period characteristics; The compensation term is superimposed on the tangential component of the reconstructed acceleration; The construction of the compensation acceleration term includes: Obtain the periodic characteristic parameters of the target orbit; Construct a combined expression that includes sine and cosine functions; The amplitude coefficients in the expression were determined by fitting historical ephemeris data; The periodic compensation term generation step includes: The formula for calculating the fundamental frequency angular velocity is: ， in, Represents the fundamental frequency angular velocity of the target orbit, in units of , The target orbital period, in seconds; Determine the harmonic order adaptively: ， in, This represents the highest order of harmonics, and its value is a positive integer. This represents the floor function. This represents the eccentricity weighting coefficient, which is dimensionless. Represents orbital eccentricity, dimensionless. This represents a highly weighted coefficient, which is dimensionless. This indicates the average altitude of the track, in kilometers. Construct a harmonic compensation model, expressed as: ， in, For periodic compensation acceleration at time Value, unit , This is a harmonic order index, dimensionless. For the first First harmonic sinusoidal amplitude coefficient, unit , It is a sine function. For the first First harmonic cosine amplitude coefficient, unit , It is a cosine function. Indicates time, in seconds (s); Batch estimation of amplitude coefficients is performed, expressed as: ， in, Represents the amplitude coefficient vector, according to Arrangement, Unit , Represents the residual observation matrix, dimension Dimensionless This represents the matrix transpose operator. Represents the column vector of historical acceleration residuals, with dimension [missing information]. ,unit , The number of residual samples is dimensionless. In the formula, the elements of the observation matrix are defined as: ， in, For matrix The Line number Column elements, dimensionless. For the first Time labels for each residual sample, in seconds; Modulated geomagnetic disturbance: ， in, The geomagnetic modulation factor is dimensionless. The geomagnetic sensitivity coefficient is dimensionless. This is the instantaneous change in the latest geomagnetic index, dimensionless. Represents the harmonic amplitude coefficient after geomagnetic modulation, in units of ; The formula for superimposing tangential component compensation is as follows: ， in, This represents the tangential acceleration vector after compensation. Represents the original tangential acceleration vector. The periodic compensation acceleration scalar obtained in step S3, in units .

2. The high-precision space target orbit prediction method based on numerical fitting as described in claim 1, characterized in that, The frequency domain decomposition in step S2 specifically includes: The wavelet packet transform algorithm is used to perform multi-scale decomposition of the acceleration observation signal with predefined basis functions; The low-frequency component corresponds to the frequency band dominated by the Earth's non-spherical gravity, the mid-frequency component corresponds to the frequency band dominated by the Sun and Moon's gravity, and the high-frequency component corresponds to the frequency band dominated by solar radiation pressure.

3. The high-precision space target orbit prediction method based on numerical fitting as described in claim 2, characterized in that, In step S3, when performing cubic spline fitting on the mid-frequency components: The distribution of fitted nodes is dynamically adjusted based on the instantaneous rate of change of the right ascension of the ascending node of the orbit. Increase node density during periods when the rate of change of right ascension at the ascending node exceeds a preset threshold.

4. The high-precision space target orbit prediction method based on numerical fitting as described in claim 1, characterized in that, The coefficient update of the compensation acceleration term includes: When the change in the real-time spatial environment index exceeds the threshold value... The coefficients were refitted using the track residual data within the sliding time window; The width setting of the sliding time window includes: The window width is dynamically adjusted based on the magnitude of sudden changes in the geomagnetic index. The magnitude of the mutation is negatively correlated with the window width.

5. The high-precision space target orbit prediction method based on numerical fitting as described in claim 1, characterized in that, In step S5, the initialization of the orbital numerical integration includes: Match a pre-stored template library based on the target orbit inclination and eccentricity; The perturbation fitting parameters of the loaded matching template are used as initial values ​​for numerical integration.

6. The high-precision space target orbit prediction method based on numerical fitting as described in claim 1, characterized in that, The processing of the high-frequency components includes: Establish a parameterized model of the equivalent area of ​​solar radiation pressure; Based on the inversion model parameters of historical luminosity curves; The establishment of the parameterized model includes: Simplify the space target into an ellipsoidal geometry; Establish the trigonometric relationship between the illuminated area and the solar incidence angle.

7. The high-precision space target orbit prediction method based on numerical fitting as described in claim 6, characterized in that, The processing steps for the high-frequency components include: The formula for calculating the angle of incidence of sunlight is: ， in, Indicates the angle of incidence of sunlight, in rad. This represents the unit vector along the principal inertial axis of the target body; it is dimensionless. Represents the unit vector of the direction of sunlight incidence, as a function of time. Change, dimensionless The equivalent illuminated area is parameterized and expressed as: ， in, Represents the equivalent area exposed to solar radiation, in units of , Represents the normal incident cross-section, in units of , This represents the quadratic modulation coefficient, with a value range of 0-1. Represents the sine function. This is the fourth-order modulation coefficient, with a value range of 0-1; The formula for calculating solar radiation pressure acceleration is: ， in, Represents the solar radiation pressure acceleration vector, in units of , The reflection and absorption coefficient is dimensionless. Represents the solar radiation pressure at 1 AU, a constant. , The target mass in space is expressed in kg. Parameter inversion based on photometric curves is expressed as follows: ， in, Let be the column vector of the three parameters to be estimated. Represents the photometric observation matrix, the first... Line element: Dimensionless Represents the historical photometric conversion flux column vector, in units of , Indicates the first Each observation epoch, in seconds; Extracting high-frequency acceleration: ， in, For high-frequency component acceleration vectors, units , Indicates to The result after applying a cubic spline smoothing with a one-orbit period, in units .

8. A high-precision space target orbit prediction method based on numerical fitting as described in any one of claims 1 to 7, characterized in that, When outputting the orbital state vector: The position and velocity vectors are in a geocentric inertial coordinate system. The timestamps are in UTC format and include millisecond-level timestamps.

Images