Method for determining initial value of spacecraft thrust parameter

Abstract

The invention relates to the technical field of satellite orbit control, and particularly discloses a spacecraft thrust parameter initial value determination method comprising the following steps: obtaining satellite observation data; selecting two arc segments with observation arc segments greater than a preset length from the satellite observation data, and performing orbit determination according to the two arc segments and a J2 item perturbation least square method to obtain an orbit determination result; comparing according to the orbit determination result to determine the thrust direction; and obtaining a corresponding parameter pre-identification method according to the thrust direction, and determining an initial value of a thrust parameter according to the parameter pre-identification method.

Description

Technical Field

[0001] The present invention relates to the technical field of satellite orbit control, and particularly relates to a method for determining the initial value of spacecraft thrust parameters. Background Technique

[0002] With the gradual application of electric propulsion technology in various space missions, such as the deployment, orbit maintenance, and deorbiting operations of low-Earth orbit mega constellations like "Starlink", the orbit determination of non-cooperative targets under continuous small-thrust maneuvers has become a complex and challenging task. Mastering its thrust parameters is of utmost importance for achieving space missions such as precise orbit determination, catalog management, and collision warning.

[0003] Different from the precise orbit determination of non-maneuvering satellites, the dynamic model of the orbit determination problem of maneuvering satellites has uncertainties. If unknown maneuvers cannot be correctly compensated, it will seriously affect the orbit determination accuracy. The Kalman filtering algorithm is widely used in the orbit determination problem of maneuvering targets. The commonly used filtering methods can be further divided into three categories: the adaptive Kalman filtering algorithm based on a single model, the multi-model algorithm, and the decision-based adaptive Kalman filtering algorithm. The filtering method can compensate for unknown maneuvers and solve the orbit determination problem of continuous thrust targets, but it cannot provide specific maneuver parameters for better subsequent target tracking, that is, maneuver reconstruction. During the maneuver reconstruction process, the initial value of the thrust parameter used is often a historical experience value or directly set to 0. This method may lead to slow convergence or non-convergence in the precise orbit determination process. For non-cooperative maneuvering spacecraft with unknown orbit and thrust parameters, a set of initial values of orbital elements and thrust parameters relatively close to the true values is required to improve the accuracy and efficiency of precise orbit determination. Summary of the Invention

[0004] Aiming at the above problems, the purpose of the present invention is to provide a method for determining the initial value of spacecraft thrust parameters. For non-cooperative maneuvering spacecraft with unknown orbital elements initial values and thrust information, a set of initial values of orbital elements and thrust parameters relatively close to the true values is provided, which is beneficial to better convergence of precise orbit determination and improves the accuracy and efficiency of precise orbit determination.

[0005] The present invention provides a method for determining the initial value of spacecraft thrust parameters, including: Obtain satellite observation data; Select two arc segments from the satellite observation data whose observation arc lengths are greater than a preset length, and perform orbit determination according to the two arc segments and the least squares method of J2 perturbation to obtain the orbit determination result; Compare according to the orbit determination result to determine the thrust direction; Obtain the corresponding parameter pre-identification method according to the thrust direction, and determine the initial value of the thrust parameter according to the parameter pre-identification method.

[0006] In a possible implementation, the orbit determination according to the least squares method based on the two arc segments and the J2 perturbation includes: Performing orbit determination according to the following formula to obtain the orbit state at the epoch moment and the optimal estimated value : : ; In the formula, is the weight matrix, is the orbit state at the epoch moment , is the actual observation value, is the theoretical observation value.

[0007] In a possible implementation, the orbit determination result is represented by the following formula: ; ; In the formula, is the first orbit determination result, is the second orbit determination result, is the first semi-major axis, is the second semi-major axis, is the first eccentricity, is the second eccentricity, is the first orbital inclination, is the second orbital inclination, is the first longitude of the ascending node, is the second longitude of the ascending node, is the first argument of perigee, is the second argument of perigee, is the first true anomaly, is the second true anomaly.

[0008] In a possible implementation, the comparison based on the orbit determination result to determine the thrust direction includes: Obtaining a preset threshold for the semi-major axis, a preset threshold for the orbital inclination, and a preset threshold for the longitude of the ascending node; When the difference between the average value of the first semi-major axis and the average value of the second semi-major axis is greater than the preset threshold for the semi-major axis, it is determined that there is a tangential thrust; When the difference between the average value of the first orbital inclination and the average value of the second orbital inclination is greater than the preset threshold for the orbital inclination, it is determined that there is a normal thrust for changing the orbital inclination; When the difference between the average value of the first ascending node longitude and the average value of the second ascending node longitude is greater than the preset threshold of the ascending node longitude, it is determined that there is a normal thrust for changing the ascending node longitude.

[0009] In a possible implementation manner, the corresponding parameter pre-identification method according to the thrust direction, and determining the initial value of the thrust parameter according to the parameter pre-identification method includes: When there is a tangential thrust, the thrust parameter is determined according to the following formula : ; In the formula, represents the cumulative amount of the satellite velocity change caused by the tangential thrust, represents the thrust action time, represents the initial semi-major axis of the satellite when the tangential thrust starts to act, represents the final semi-major axis after the tangential thrust acts for a period of time, is the gravitational constant of the central body, is the tangential thrust acceleration.

[0010] In a possible implementation manner, when there is a normal thrust for changing the orbit inclination, the thrust parameter is determined according to the following formula : ; In the formula, is the normal thrust acceleration, is the current value of the orbit inclination, is the initial value of the orbit inclination, is the average angular velocity of the orbit, is the thrust action time.

[0011] In a possible implementation manner, when there is a normal thrust for changing the ascending node longitude, the thrust parameter is determined according to the following formula : ; In the formula, is the change amount of the right ascension of the ascending node caused by the normal thrust, is the initial value of the right ascension of the ascending node, is the change rate of the right ascension of the ascending node under the influence of the J2 term.

[0012] In a possible implementation manner, the cumulative amount of the satellite velocity change caused by the tangential thrust is determined according to the following formula : ; ; In the formula, is the semi-major axis of the orbit.

[0013] In a possible implementation, the rate of change of the right ascension of the ascending node under the influence of the J2 term is determined according to the following formula : ; In the formula, is the first term in the zonal harmonic terms of the Earth's non-spherical gravitational field, is the semi-latus rectum.

[0014] The method for determining the initial values of the spacecraft thrust parameters provided by the present invention can provide a set of initial values of the orbital elements and thrust parameters that are relatively close to the true values in the case of unknown initial values of the orbital elements and thrust information for a non-cooperative maneuvering spacecraft with unknown orbital and thrust parameters. This is beneficial for better convergence of precise orbit determination and improves the accuracy and efficiency of precise orbit determination. BRIEF DESCRIPTION OF THE DRAWINGS

[0015] Figure 1 is a schematic flow chart of the method for determining the initial values provided by the embodiment of the present invention; Figure 2 is a diagram showing the change of the semi-major axis under the action of the tangential thrust provided by the embodiment of the present invention; Figure 3 is a diagram showing the change of the eccentricity under the action of the tangential thrust provided by the embodiment of the present invention; Figure 4 is a diagram showing the change of the orbital inclination under the action of the normal thrust provided by the embodiment of the present invention; Figure 5 is a diagram showing the change of the right ascension of the ascending node under the action of the normal thrust provided by the embodiment of the present invention; Figure 6 is a schematic diagram of two observation arcs provided by the embodiment of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0016] The following further describes in detail the embodiments of the present invention in conjunction with the drawings and examples. The detailed description and drawings of the following examples are used to exemplarily illustrate the principle of the present invention, but cannot be used to limit the scope of the present invention, that is, the present invention is not limited to the described preferred embodiments, and the scope of the present invention is defined by the claims.

[0017] In the description of the present invention, it should be noted that unless otherwise specified, the meaning of "a plurality of" is two or more; the terms "first", "second", etc. are only used for descriptive purposes and cannot be understood as indicating or implying relative importance; for those of ordinary skill in the art, the specific meanings of the above terms in the present invention can be understood according to specific circumstances.

[0018] Figure 1Schematic flowchart of the initial value determination method provided by the embodiments of the present invention, as Figure 1 shown, the present invention provides a method for determining the initial value of spacecraft thrust parameters, including: Step S1, obtaining satellite observation data; Step S2, selecting two arc segments with an observation arc length greater than a preset length from the satellite observation data, and performing orbit determination based on the two arc segments and the least squares method of J2 perturbation to obtain an orbit determination result; Initial orbit determination usually considers the two-body dynamics model, and its result provides an initial value for precise orbit determination. Optimal estimation theory is not used in initial orbit determination. For radar observation data, common initial orbit determination algorithms include the Gibbs method, the Herrick-Gibbs method, etc.; for angle-only measurement data, common initial orbit determination algorithms include the Laplace method, the Gauss method, the Gooding method, etc.

[0019] Initial orbit determination usually uses single-arc measurement data. However, a single observation arc segment is usually short (the arc segment of a low-orbit target is usually only a few minutes to more than ten minutes), and the orbit determination accuracy is limited. In the process of continuous thrust control, the thrust is small and the rate of change of orbital elements is slow. If the orbit determination accuracy is too low, it will lead to a relatively large calculation error of thrust acceleration. To improve the single-arc orbit determination accuracy, the least squares method considering J2 perturbation is used for orbit determination.

[0020] Define the n group of noisy measurement values of distance, azimuth, and elevation angle obtained by ground radar observation equipment as: (1) (2) Let the true value of the measurement sequence be: (3) (4) Assume that the observation error follows a normal distribution with a mean square deviation of , then the statistical characteristics of the observed values are: (5) Define the state quantity of the dynamic system as , where: (6) The system observation model is: (7) Let: (8) The least squares orbit determination problem is described as using the least squares criterion to solve the orbit state at a certain epoch The optimal estimated value that minimizes the weighted sum of squares of residuals between the theoretical and actual observed values.

[0021] Orbit determination is performed according to the following formula to obtain the orbital state at the epoch of the optimal estimated value : (9) where is the weight matrix given considering the accuracies of different observation data. Usually, the reciprocals of the observation data errors are used for weighting, is the orbital state at the epoch , is the actual observed value, is the theoretical observed value.

[0022] The Gaussian perturbation equations for orbital motion considering continuous small thrust are: (10) where are the six orbital elements, , , are the thrust accelerations in the tangential, radial, and normal directions, respectively.

[0023] Currently, the vast majority of on-orbit spacecraft maneuvers are carried out within the orbital plane, manifested as the adjustment of the semi-major axis by tangential acceleration. This method is particularly widely applied during the orbital climb process in the satellite's orbit insertion phase. To improve the orbital climb efficiency, usually only tangential thrust is applied, and the radial thrust is zero. Although applying normal thrust to change the orbital plane consumes a lot of fuel and is rarely used in reality, it does not rule out the situation of adjusting the orbital inclination or the right ascension of the ascending node by applying normal thrust. Therefore, the identification of continuous small thrust parameters includes thrust direction identification and thrust magnitude identification. First, judge the direction of the continuous small thrust, and then adopt different parameter pre-identification methods for thrusts in different directions.

[0024] When , under the action of tangential acceleration, the semi-major axis continuously increases, and the eccentricity changes periodically within a small range, which can be ignored. As shown in Equation (11).

[0025] (11) For a nearly circular orbit with continuous tangential thrust, the thrust acceleration only affects the semi-major axis, and Equation (11) can be simplified to: (12) Simulate and analyze the variations of the semi-major axis and eccentricity under the action of tangential acceleration. The thrust parameters are set with reference to "Starlink" V2.0 mini, the satellite mass is 750 kg, equipped with a Hall electric thruster with a thrust of 170 mN, and the total thrust acceleration generated is 2.26×10 -4 m / s 2 , the initial semi-major axis of the satellite, and the initial eccentricity is 0. Other simulation conditions are shown in Table 1.

[0026] When , the normal acceleration affects the orbital inclination and the right ascension of the ascending node, as shown in Equation (13).

[0027] (13) For a nearly circular orbit, Equation (13) is simplified to: (14) Simulate and analyze the variations of the orbital inclination and the right ascension of the ascending node under the action of normal acceleration. The initial values of the orbital inclination and the right ascension of the ascending node are set to 53° and 353° respectively, and the thrust acceleration is 2.26×10 -4 m / s 2 , and other simulation conditions are the same as those in Table 1.

[0028] If the direction of the normal acceleration remains constant, the orbital inclination and the right ascension of the ascending node will vary periodically, and the integral within one orbital period is 0. To improve the thrust efficiency, if the normal thrust is applied with the aim of changing the orbital inclination, the thrust direction should be reversed at the latitude argument of and to make the orbital inclination change continuously, as shown in Figure 4 , If the normal thrust is applied with the aim of changing the right ascension of the ascending node, the thrust direction should be reversed at the latitude argument of 0 and π to make the right ascension of the ascending node change continuously. To eliminate the influence of the J2 perturbation and more intuitively display the variation of the right ascension of the ascending node, let , represent the average value of the right ascension of the ascending node at the moment t without thrust, and the variation of is as shown in Figure 5 .

[0029] In summary, the tangential thrust parameters can be inversely solved by analyzing the variation of the orbital semi-major axis, and the normal thrust parameters can be inversely solved by analyzing the variations of the orbital inclination and the right ascension of the ascending node.

[0030] The pre-identification results of the continuous thrust parameters are inversely solved from the radar observation data at two intervals with a certain time duration. Let the initial times of the two radar observation arcs be , respectively, and the end times be , , the observation durations are respectively , . The time interval between is denoted as Δ t , as shown in Figure 6 . The orbit determination is respectively performed based on two radar arcs by applying the method in Section II.

[0031] In a possible implementation, the orbit determination result is represented by the following formula: ; ; In the formula, is the first orbit determination result, is the second orbit determination result, is the first semi-major axis, is the second semi-major axis, is the first eccentricity, is the second eccentricity, is the first orbital inclination, is the second orbital inclination, is the first longitude of the ascending node, is the second longitude of the ascending node, is the first argument of perigee, is the second argument of perigee, is the first true anomaly, is the second true anomaly.

[0032] Step S3, compare according to the orbit determination result to determine the thrust direction; In a possible implementation, obtain the preset threshold of the semi-major axis, the preset threshold of the orbital inclination, and the preset threshold of the longitude of the ascending node; When the difference between the average value of the first semi-major axis and the average value of the second semi-major axis is greater than the preset threshold of the semi-major axis, it is determined that there is a tangential thrust; the variation of the semi-major axis and the eccentricity under the action of the tangential thrust is as shown in Figure 2 and Figure 3 .

[0033] When the difference between the average value of the first orbital inclination and the average value of the second orbital inclination is greater than the preset threshold of the orbital inclination, it is determined that there is a normal thrust for changing the orbital inclination; When the difference between the average value of the first longitude of the ascending node and the average value of the second longitude of the ascending node is greater than the preset threshold of the longitude of the ascending node, it is determined that there is a normal thrust for changing the longitude of the ascending node.

[0034] In a possible implementation, convert the orbital elements into mean orbital elements, denoted as: ; ; Compare the two orbit determination results, and focus on the changes in the semi-major axis, eccentricity, orbital inclination, and right ascension of the ascending node. Set the thresholds to respectively. If , it is considered that there is a tangential thrust. If or , then it is considered that there is a normal thrust. Since the right ascension of the ascending node is greatly affected by perturbations, is a variable related to Δ t . Different parameter pre-identification methods are used for different direction thrusts.

[0035] Step S4: Obtain the corresponding parameter pre-identification method according to the thrust direction, and determine the initial value of the thrust parameter according to the parameter pre-identification method.

[0036] In a possible implementation, when there is a tangential thrust, the thrust parameter is determined according to the following formula : ; In the formula, represents the cumulative amount of the satellite velocity change caused by the tangential thrust, represents the thrust action time, represents the initial semi-major axis of the satellite when the tangential thrust starts to act, represents the final semi-major axis after the tangential thrust acts for a period of time, is the gravitational constant of the central body, is the tangential thrust acceleration.

[0037] In a possible implementation, when there is a normal thrust that changes the orbital inclination, the thrust parameter is determined according to the following formula : ; In the formula, is the normal thrust acceleration, is the current value of the orbital inclination, is the initial value of the orbital inclination, is the average orbital angular velocity, , is the thrust action time.

[0038] In a possible implementation, when there is a normal thrust that changes the longitude of the ascending node, the thrust parameter is determined according to the following formula : ; In the formula, is the change in the right ascension of the ascending node caused by the normal thrust, is the initial value of the right ascension of the ascending node, is the rate of change of the right ascension of the ascending node under the influence of the J2 term.

[0039] In a possible implementation, the cumulative amount of the satellite velocity change caused by the tangential thrust is determined according to the following formula : ; ; wherein, is the semi-major axis of the orbit.

[0040] In a possible implementation, the rate of change of the right ascension of the ascending node under the influence of the J2 term is determined according to the following formula : ; wherein, is the first term in the zonal harmonic terms of the Earth's non-spherical gravitational field, is the semi-latus rectum.

[0041] Furthermore, different thrust parameter pre-identification methods are further explained.

[0042] Continuous tangential thrust parameter pre-identification method: The arc segment effect of continuous thrust is reflected by the line integral of Equation (12): (15) (16) wherein, represents the cumulative amount of the satellite velocity change caused by the tangential thrust, and does not represent the direction and magnitude of the actual satellite velocity, represents the initial semi-major axis of the satellite when the tangential thrust starts to act, represents the final semi-major axis after the tangential thrust has acted for a period of time. Let the thrust action time be , then: (17) It can be seen from Equation (17) that the tangential acceleration can be solved from the initial semi-major axis, the final semi-major axis of the arc segment of the tangential thrust action, and the thrust action time.

[0043] Normal thrust parameter pre-identification method for changing the orbit inclination: For a nearly circular orbit, if continuous normal thrust is applied to change the orbit inclination ( at the latitude argument of and in the reverse direction), integrate the first equation of Equation (14) over , : (18) It can be obtained that: (19) If , that is, the semi-major axis a remains unchanged, the orbital inclination will change linearly under the action of continuous normal thrust: (20) If , the semi-major axis a changes under the action of tangential thrust. From Equation (16): (21) Substitute Equation (21) into Equation (19), and it can be obtained that: (22) Integrating both sides gives: (23) is the integration constant. When t = 0, , so: (24) (25) According to Equation (25), the normal acceleration can be solved from the change in orbital inclination, the magnitude of tangential acceleration, and the time interval between two arcs.

[0044] Pre-identification of normal thrust parameters for changing the right ascension of the ascending node: If continuous normal thrust is applied to change the right ascension of the ascending node ( reversing at the argument of latitude of 0 and π), integrate the first equation of Equation (14) over [0, π]: (26) It can be obtained that: (27) Neglect the small-amplitude periodic change in orbital inclination. If , that is, the semi-major axis a remains unchanged, the right ascension of the ascending node will change linearly under the action of continuous normal thrust: (28) If , the semi-major axis a changes under the action of tangential thrust. Following the derivation process in Section 3.2.1, it can be obtained that: (29) Equation (29) describes the change in the right ascension of the ascending node caused by the combined action of the normal acceleration and the tangential acceleration. In fact, the change in the right ascension of the ascending node is greatly affected by perturbations, and the main influencing factor is the perturbation of the J2 term of the Earth's non-spherical gravity. For the long-term perturbation effect, only the influence of the J2 term needs to be considered. Under the influence of the J2 term perturbation, the first-order long-term change rate of the right ascension of the ascending node is: (30) The calculation unit adopts the normalized unit, and the corresponding gravitational coefficient of the central celestial body μ = GM = 1. Let be the corresponding value of the mean elements at time t of the right ascension of the ascending node, then: (31) According to Equation (31), the normal acceleration can be solved from the change in the right ascension of the ascending node, the magnitude of the tangential acceleration, and the time interval between the two arcs.

[0045] For easy understanding, the following is a simulation analysis of the scenario of changing the satellite orbit elements under continuous small thrust.

[0046] Suppose a certain satellite has an initial orbit altitude of 380 km, and continuous small thrust is applied at 12:00 on January 25, 2024. The thrust parameter settings refer to "Starlink" V2.0mini, the specific impulse of the thruster is 2500 s, the electric thrust magnitude is 170 mN, and the thrust acceleration is 2.26×10-4 m / s2. The initial orbit elements of the satellite ignition point are shown in Table 2, and other simulation parameter settings are shown in Table 1.

[0047] Table 1

[0048] Table 2

[0049] Assume that radar devices located at different locations can obtain satellite measurement data, and the distribution position information of the radar devices is shown in Table 3. Assume that the radar range measurement accuracy is 30 m (1σ), the angle measurement accuracy is 0.01° (1σ), and the minimum elevation angle is 5°. The measurement data (azimuth, elevation angle, and range) of the radar for this satellite are obtained through simulation.

[0050] Table 3

[0051] The thrust acceleration result is obtained by inverse solution from the radar observation data at two intervals with a certain time duration. Let the initial moments of the two radar observations be 、 , and the end moments be 、 , and the observation durations be 、 According to the different positional relationships between the satellite and the measurement station, the length of the visible arc segment of a single measurement station varies from 3 to 9 minutes. To ensure the accuracy of the semi-major axis solution, only the arc segments with an observation arc length greater than 6 minutes are selected for orbit determination.

[0052] In the pre-identification of continuous tangential thrust parameters, the time of the two observation arc segments is shown in Table 4.

[0053] Table 4

[0054] Perform 100 simulations on the two observations respectively. Apply different random observation errors to each simulation. Use the single-arc orbit determination method described in Section 3. The semi-major axes of the two arc segments are respectively denoted as 、 。

[0055] According to Equation (17), from and Inverse-solve the tangential acceleration. Solve the tangential acceleration according to the two arc segments with an interval of 23 hours. The average relative error is only 3.2%.

[0056] In the pre-identification of continuous normal thrust parameters for changing the orbit inclination, the time of the two observation arc segments is shown in Table 5.

[0057] Table 5

[0058] Perform 100 simulations on the two observations respectively. Apply different random observation errors to each simulation. Use the single-arc orbit determination method described in Section 3. The orbit inclinations of the two arc segments are respectively denoted as 、 。

[0059] According to Equation (17), from and Inverse-solve the normal acceleration. Solve the normal acceleration according to the two arc segments with an interval of 23 hours. The average relative error is 28.5%.

[0060] In the pre-identification of continuous normal thrust parameters for changing the right ascension of the ascending node, the time of the two observation arc segments is shown in Table 6.

[0061] Table 6

[0062] Perform 100 simulations on the two observations respectively. Apply different random observation errors to each simulation. Use the single-arc orbit determination method described in Section 3. The right ascensions of the ascending nodes of the two arc segments are respectively denoted as 、 。

[0063] According to Equation (17), from and the normal acceleration is solved inversely. The normal acceleration is solved based on two arc segments spaced 23 hours apart, and the average relative error is 30.21%.

[0064] The method for determining the initial values of the spacecraft thrust parameters provided by the present invention, for a non-cooperative maneuvering spacecraft with unknown orbit and thrust parameters, provides a set of initial values of orbit elements and thrust parameters that are relatively close to the true values in the case of unknown initial values of orbit elements and thrust information, which is beneficial to better convergence of precise orbit determination and improves the accuracy and efficiency of precise orbit determination.

[0065] The above is only the specific implementation manner of the present invention, but the protection scope of the present invention is not limited thereto. Any person skilled in the art within the technical scope disclosed by the present invention can easily think of changes or substitutions, which should all be covered within the protection scope of the present invention. Therefore, the protection scope of the present invention shall be subject to the protection scope of the claims.

Claims

1. A method for determining the initial value of the thrust parameter of a spacecraft, characterized in that Including: Obtaining satellite observation data; Selecting two arc segments with observation arcs longer than a preset length from the satellite observation data, and performing orbit determination according to the two arc segments and the least squares method of J2 perturbation to obtain an orbit determination result; Comparing according to the orbit determination result to determine the thrust direction; Obtaining a corresponding parameter pre-identification method according to the thrust direction, and determining an initial value of the thrust parameter according to the parameter pre-identification method.

2. The initial value determination method according to claim 1, wherein The performing orbit determination according to the two arc segments and the least squares method of J2 perturbation includes: Orbit determination is carried out according to the following formula to obtain the orbit state at the epoch of the optimal estimate : ； In the formula, is the weight matrix, is the orbital state at epoch , is the actual observation value, is the theoretical observation value.

3. The initial value determination method according to claim 1, wherein Also including: Representing the orbit determination result according to the following formula: ； ； Wherein, is the first orbit determination result, is the second orbit determination result, is the first semi-major axis, is the second semi-major axis, is the first eccentricity, is the second eccentricity, is the first orbital inclination, is the second orbital inclination, is the first longitude of the ascending node, is the second longitude of the ascending node, is the first argument of periapsis, is the second argument of periapsis, is the first true anomaly, is the second true anomaly.

4. The initial value determination method according to claim 3, wherein The comparing according to the orbit determination result to determine the thrust direction includes: Obtaining a preset threshold for semi-major axis, a preset threshold for orbital inclination, and a preset threshold for longitude of ascending node; When the difference between the average value of the first semi-major axis and the average value of the second semi-major axis is greater than the preset threshold for semi-major axis, it is determined that there is a tangential thrust; When the difference between the average value of the first orbital inclination and the average value of the second orbital inclination is greater than the preset threshold for orbital inclination, it is determined that there is a normal thrust for changing the orbital inclination; When the difference between the average value of the first longitude of ascending node and the average value of the second longitude of ascending node is greater than the preset threshold for longitude of ascending node, it is determined that there is a normal thrust for changing the longitude of ascending node.

5. The initial value determination method according to claim 1, characterized in that The obtaining a corresponding parameter pre-identification method according to the thrust direction and determining an initial value of the thrust parameter according to the parameter pre-identification method includes: When there is a tangential thrust, the thrust parameter is determined according to the following formula :[[]]END]] ； In the formula, represents the cumulative amount of the satellite velocity change caused by the tangential thrust, represents the thrust acting time, represents the initial semi-major axis of the satellite when the tangential thrust starts to act, represents the final semi-major axis after the tangential thrust acts for a period of time, is the gravitational constant of the central celestial body, is the tangential thrust acceleration.

6. The initial value determination method according to claim 5, wherein Also including: When there is a normal thrust for changing the orbital inclination, the thrust parameter is determined according to the following formula : ； In the formula, is the normal thrust acceleration, is the current value of the orbit inclination, is the initial value of the orbit inclination, is the average angular velocity of the orbit, is the thrust application time.

7. The initial value determination method according to claim 6, characterized in that Also including: When there is a normal thrust that changes the longitude of the ascending node, the thrust parameter is determined according to the following formula : ； In the formula, is the change in the right ascension of the ascending node caused by the normal thrust, is the initial value of the right ascension of the ascending node, is the rate of change of the right ascension of the ascending node under the influence of the J2 term.

8. The initial value determination method according to claim 5, characterized in that Also including: Determine the cumulative amount of change in the satellite's velocity caused by the tangential thrust according to the following formula : ； ； In the formula, is the semi-major axis of the orbit.

9. The initial value determination method according to claim 7, wherein Also including: Determine the rate of change of the right ascension of the ascending node under the influence of the J2 term according to the following formula : ； In the formula, is the first term in the zonal harmonics of the Earth's non-spherical gravity, is the semi-latus rectum.

Images