A method for determining initial values ​​of spacecraft thrust parameters

Abstract

The present invention relates to the field of satellite orbit control technology, and specifically discloses a method for determining initial values ​​of spacecraft thrust parameters, comprising: acquiring 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 based on the two arc segments and a least squares method of J2 perturbation to obtain an orbit determination result; performing a comparison based on the orbit determination results to determine a thrust direction; acquiring a corresponding parameter pre-identification method based on the thrust direction, and determining an initial value of the thrust parameter based on the parameter pre-identification method.

Description

Technical Field

[0001] The present invention relates to the technical field of satellite orbit control, and in particular to a method for determining initial values ​​of spacecraft thrust parameters. Background Art

[0002] As electric propulsion technology is gradually applied to various space missions, such as the deployment, orbit maintenance and deorbit operations of low-orbit giant constellations such as "Starlink", orbit determination of non-cooperative targets under continuous low-thrust maneuvers has become a complex and challenging task. Mastering its thrust parameters is of paramount importance for achieving space missions such as precise orbit determination, catalog management, and collision warning.

[0003] Unlike precise orbit determination for non-maneuvering satellites, orbit determination for maneuvering satellites involves uncertainties in the dynamics model. Failure to properly compensate for unknown maneuvers will severely impact orbit determination accuracy. Kalman filtering algorithms are widely used for orbit determination of maneuvering targets. Common filtering methods can be further categorized into three categories: single-model adaptive Kalman filtering, multi-model algorithms, and decision-based adaptive Kalman filtering. These filtering methods can compensate for unknown maneuvers and address the orbit determination problem for continuous thrust targets, but they cannot provide specific maneuver parameters for subsequent tracking, known as maneuver reconstruction. During maneuver reconstruction, the initial values ​​of the thrust parameters used are often based on historical experience or simply set to zero. This approach can result in slow or even non-convergence in the precise orbit determination process. For non-cooperative maneuvering spacecraft with unknown orbital and thrust parameters, a set of initial values ​​for orbital elements and thrust parameters that are relatively close to the true values ​​is required to improve the accuracy and efficiency of precise orbit determination. Summary of the Invention

[0004] In response to the above problems, the purpose of the present invention is to provide a method for determining the initial values ​​of spacecraft thrust parameters. For non-cooperative maneuvering spacecraft whose orbit and thrust parameters are unknown, a set of initial values ​​of orbital elements and thrust parameters that are relatively close to the true values ​​is provided under the condition of unknown initial values ​​of orbital elements and thrust information, which is conducive 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 initial values ​​of spacecraft thrust parameters, comprising:

[0006] Obtain satellite observation data;

[0007] Selecting two arcs with observation arcs greater than a preset length from the satellite observation data, and performing orbit determination based on the two arcs and the least squares method of the J2 perturbation to obtain an orbit determination result;

[0008] determining the thrust direction based on the orbit determination results by comparison;

[0009] A corresponding parameter pre-identification method is obtained according to the thrust direction, and an initial value of the thrust parameter is determined according to the parameter pre-identification method.

[0010] In a possible implementation, performing orbit determination based on the two arc segments and the J2 perturbation using a least squares method includes:

[0011] The orbit is determined according to the following formula to obtain the epoch time Track status The best estimate of :

[0012] ;

[0013] Where, is the weight matrix, epoch time The orbital state, is the actual observed value, is the theoretical observation value.

[0014] In one possible implementation, the orbit determination result is expressed according to the following formula:

[0015] ;

[0016] ;

[0017] Where, 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 inclination of the second orbit, is the longitude of the first ascending node, is the longitude of the second ascending node, is the first pericenter argument, is the second pericenter argument, is the first true anomaly, is the second true anomaly.

[0018] In a possible implementation, determining the thrust direction by comparing the orbit determination results includes:

[0019] Obtain the preset thresholds for the semi-major axis, orbit inclination, and ascending node longitude;

[0020] 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 value of the semi-major axis, it is determined that a tangential thrust exists;

[0021] 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 value of the orbital inclination, determining that a normal thrust for changing the orbital inclination exists;

[0022] 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 value of the ascending node longitude, it is determined that there is a normal thrust that changes the ascending node longitude.

[0023] In a possible implementation, 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 include:

[0024] When there is tangential thrust, the thrust parameters are determined according to the following formula: :

[0025] ;

[0026] Where, It represents the cumulative change of satellite velocity caused by tangential thrust. Indicates the thrust action time, represents the initial semi-major axis of the satellite when the tangential thrust begins to act, represents the final semi-major axis after a period of tangential thrust, is the gravitational constant of the central celestial body, is the tangential thrust acceleration.

[0027] In one possible implementation, when there is a normal thrust that changes the orbital inclination, the thrust parameter is determined according to the following formula: :

[0028] ;

[0029] Where, is the normal thrust acceleration, is the current value of the orbital inclination, is the initial value of orbital inclination, is the average angular velocity of the orbit, is the thrust action time.

[0030] In one 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: :

[0031] ;

[0032] Where, is the change in 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 J2.

[0033] In one possible implementation, the cumulative change in satellite velocity due to the tangential thrust is determined according to the following formula: :

[0034] ;

[0035] ;

[0036] Where, is the semi-major axis of the orbit.

[0037] In one possible implementation, the rate of change of the right ascension of the ascending node under the influence of the J2 term is determined by the following formula: :

[0038] ;

[0039] Where, is the first harmonic term of the Earth's non-spherical gravitational belt, It is a semi-diameter.

[0040] The method for determining the initial values ​​of spacecraft thrust parameters provided by the present invention provides a set of initial values ​​of orbital elements and thrust parameters that are relatively close to the true values ​​for a non-cooperative maneuvering spacecraft whose orbit and thrust parameters are unknown, under the condition of unknown initial values ​​of orbital elements and thrust information, which is conducive to better convergence of precise orbit determination and improves the accuracy and efficiency of precise orbit determination. BRIEF DESCRIPTION OF THE DRAWINGS

[0041] Figure 1 A schematic flow chart of a method for determining an initial value provided by an embodiment of the present invention;

[0042] Figure 2 A diagram showing the change of the semi-major axis under the action of tangential thrust provided by an embodiment of the present invention;

[0043] Figure 3 A diagram showing changes in eccentricity under tangential thrust provided by an embodiment of the present invention;

[0044] Figure 4 A diagram showing changes in orbital inclination under normal thrust provided by an embodiment of the present invention;

[0045] Figure 5 A diagram showing changes in the right ascension of the ascending node under the action of normal thrust provided by an embodiment of the present invention;

[0046] Figure 6 A schematic diagram of two observation arc segments provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0047] The following detailed description of the embodiments of the present invention is provided in conjunction with the accompanying drawings and examples. The following detailed description of the embodiments and the accompanying drawings are intended to illustrate the principles of the present invention, but are not intended to limit the scope of the present invention. That is, the present invention is not limited to the preferred embodiments described, and the scope of the present invention is defined by the claims.

[0048] In the description of the present invention, it should be noted that, unless otherwise specified, “plurality” means two or more; the terms “first”, “second”, etc. are used for descriptive purposes only and cannot be understood as indicating or implying relative importance; for ordinary technicians in this field, the specific meanings of the above terms in the present invention can be understood according to the specific circumstances.

[0049] Figure 1 A flow chart of a method for determining an initial value provided by an embodiment of the present invention is shown in FIG. Figure 1 As shown, the present invention provides a method for determining the initial value of a spacecraft thrust parameter, comprising:

[0050] Step S1, obtaining satellite observation data;

[0051] Step S2: Select two arcs with observation arcs greater than a preset length from the satellite observation data, and perform orbit determination based on the two arcs and the least squares method of the J2 perturbation to obtain an orbit determination result;

[0052] Initial orbit determination typically considers a two-body dynamics model, and the results provide initial values ​​for precise orbit determination. Initial orbit determination does not use optimal estimation theory. For radar observation data, commonly used initial orbit determination algorithms include the Gibbs method and the Herrick-Gibbs method. For angle measurement data only, commonly used initial orbit determination algorithms include the Laplace method, the Gaussian method, and the Gooding method.

[0053] Initial orbit determination typically uses single-arc segment measurement data. However, a single observation arc is typically short (typically only a few minutes to a dozen minutes for a low-orbit target), limiting orbit determination accuracy. Continuous thrust control (CTC) requires low thrust during orbit changes, resulting in slow orbital element changes. Low orbit determination accuracy can lead to large relative errors in thrust acceleration calculations. To improve single-arc segment orbit determination accuracy, a least-squares method that considers the J2 perturbation is employed.

[0054] Define the distance, azimuth, and elevation angles obtained by ground radar observation equipment n The group noisy measurements are:

[0055] (1)

[0056] (2)

[0057] Assume the true value of the measurement sequence is:

[0058] (3)

[0059] (4)

[0060] Assuming that the observation error follows a normal distribution, the mean square error is , then the statistical characteristics of the observations are:

[0061] (5)

[0062] The state quantity of the dynamic system is defined as ,in:

[0063] (6)

[0064] The system observation model is:

[0065] (7)

[0066] set up:

[0067] (8)

[0068] The least squares orbit determination problem is described as using the least squares criterion to solve the time at a certain epoch Track status The optimal estimate of , which minimizes the sum of squares of the weighted residuals between the theoretical observations and the actual observations.

[0069] The orbit is determined according to the following formula to obtain the epoch time Track status The best estimate of :

[0070] (9)

[0071] Where, In order to consider the weight matrix given by different observation data accuracy, the observation data error is usually used. The inverse weight of epoch time The orbital state, is the actual observed value, is the theoretical observation value.

[0072] The Gaussian perturbation equation of orbital motion considering continuous small thrust is:

[0073] (10)

[0074] in, is the six element number of the orbit, , , are the thrust acceleration in the tangential, radial and normal directions respectively.

[0075] At present, most maneuvers of on-orbit spacecraft are carried out within the orbital plane, manifested as adjustments to the semi-major axis by tangential acceleration. This method is particularly widely used in the orbit climb process during the satellite's orbit insertion phase. To improve the efficiency of orbit climb, thrust is usually applied only in the tangential direction, 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 possibility 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, the direction of the continuous small thrust is determined, and then different parameter pre-identification methods are used for thrust in different directions.

[0076] when When , under the action of tangential acceleration, the semi-major axis continues to increase, and the eccentricity changes periodically in a small range and can be ignored. As shown in Equation (11).

[0077] (11)

[0078] For a near-circular orbit with continuous tangential thrust, the thrust acceleration only affects the semi-major axis, and Equation (11) can be simplified to:

[0079] (12)

[0080] The changes of the semi-major axis and eccentricity under the action of tangential acceleration were simulated and analyzed. The thrust parameter setting was based on the "Starlink" V2.0mini. The satellite mass was 750kg and it was equipped with a Hall electric thruster with a thrust of 170mN. The total thrust acceleration generated was 2.26×10 -4 m / s 2 , the initial semi-major axis of the satellite and the initial eccentricity are 0. Other simulation conditions are shown in Table 1.

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

[0082] (13)

[0083] For a near-circular orbit, Equation (13) is simplified to:

[0084] (14)

[0085] The changes of orbital inclination and right ascension of ascending node under the action of normal acceleration are simulated and analyzed. The initial values ​​of orbital inclination and right ascension of ascending node are set to 53° and 353° respectively, and the thrust acceleration is 2.26×10 -4 m / s 2 ,Other simulation conditions are the same as Table 1.

[0086] If the direction of normal acceleration remains constant, the orbital inclination and right ascension of the ascending node will change periodically, and the integral within one orbital period will be 0. In order to improve the efficiency of thrust, if normal thrust is applied for the purpose of changing the orbital inclination, the thrust direction should be in the direction of the latitude argument. and The inclination of the orbit changes continuously, such as Figure 4 As shown,

[0087] If the normal thrust is applied for the purpose of changing the right ascension of the ascending node, the thrust direction should be reversed at the latitude arguments of 0 and π, so that the right ascension of the ascending node changes continuously. In order to eliminate the influence of the J2 perturbation and more intuitively show the change of the right ascension of the ascending node, it is assumed that , Indicates that there is no thrust. t The average right ascension of the ascending node at the moment, The changes in Figure 5 shown.

[0088] In summary, the tangential thrust parameters can be inversely solved by analyzing the changes in the semi-major axis of the orbit, and the normal thrust parameters can be inversely solved by analyzing the changes in the orbit inclination and the right ascension of the ascending node.

[0089] The continuous thrust parameter pre-identification results are obtained by inverse analysis of two radar observation data with a certain time interval. Assume that the initial moments of the two radar observation arcs are 、 The end times are 、 The observation time is 、 , and The time interval between them is denoted as Δ t ,like Figure 6 The method in Section 2 is used to determine the orbit based on two radar arcs.

[0090] In one possible implementation, the orbit determination result is expressed according to the following formula:

[0091] ;

[0092] ;

[0093] Where, 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 inclination of the second orbit, is the longitude of the first ascending node, is the longitude of the second ascending node, is the first pericenter argument, is the second pericenter argument, is the first true anomaly, is the second true anomaly.

[0094] Step S3, comparing the orbit determination results to determine the thrust direction;

[0095] In a possible implementation, a semi-major axis preset threshold, an orbit inclination preset threshold, and an ascending node longitude preset threshold are obtained;

[0096] 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 judged that there is tangential thrust; the changes of the semi-major axis and eccentricity under the action of tangential thrust are as follows Figure 2 and Figure 3 shown.

[0097] When the difference between the average value of the first orbital inclination and the average value of the second orbital inclination is greater than a preset threshold value of the orbital inclination, it is determined that a normal thrust for changing the orbital inclination exists;

[0098] 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 a preset threshold value of the ascending node longitude, it is determined that there is a normal thrust that changes the ascending node longitude.

[0099] In one possible implementation, the number of orbital elements is converted to the average number of orbital elements, which is expressed as:

[0100] ;

[0101] ;

[0102] Comparing the two orbit determination results, we focus on the changes in the semi-major axis, eccentricity, orbit inclination and right ascension of the ascending node, and set the thresholds to ,like , it is assumed that there is a tangential thrust, if or , then it is considered that there is normal thrust. Since the ascending node right ascension is greatly affected by perturbations, is with Δ t Different parameter pre-identification methods are used for thrust in different directions.

[0103] Step S4: 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.

[0104] In one possible implementation, when there is tangential thrust, the thrust parameter is determined according to the following formula: :

[0105] ;

[0106] Where, It represents the cumulative change of satellite velocity caused by tangential thrust. Indicates the thrust action time, represents the initial semi-major axis of the satellite when the tangential thrust begins to act, represents the final semi-major axis after a period of tangential thrust, is the gravitational constant of the central celestial body, is the tangential thrust acceleration.

[0107] In one possible implementation, when there is a normal thrust that changes the orbital inclination, the thrust parameter is determined according to the following formula: :

[0108] ;

[0109] Where, is the normal thrust acceleration, is the current value of the orbital inclination, is the initial value of orbital inclination, is the average angular velocity of the orbit, , is the thrust action time.

[0110] In one 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: :

[0111] ;

[0112] Where, is the change in 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 J2.

[0113] In one possible implementation, the cumulative change in satellite velocity due to the tangential thrust is determined according to the following formula: :

[0114] ;

[0115] ;

[0116] Where, is the semi-major axis of the orbit.

[0117] In one possible implementation, the rate of change of the right ascension of the ascending node under the influence of the J2 term is determined by the following formula: :

[0118] ;

[0119] Where, is the first harmonic term of the Earth's non-spherical gravitational belt, It is a semi-diameter.

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

[0121] Continuous tangential thrust parameter pre-identification method: The curve integral of formula (12) is used to reflect the arc effect of continuous thrust:

[0122] (15)

[0123] (16)

[0124] Where, It indicates the cumulative change in satellite velocity 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 begins to act, It represents the final semi-major axis after a period of tangential thrust. Assume that the thrust action time is ,but:

[0125] (17)

[0126] From formula (17), it can be seen that the tangential acceleration can be solved by the initial semi-major axis, final semi-major axis and thrust action time of the tangential thrust action arc.

[0127] Normal thrust parameter pre-identification method for changing orbital inclination: For a near-circular orbit, if a continuous normal thrust is applied to change the orbital inclination ( The latitude angle is and In the reverse direction), the first equation of equation (14) is [ , ]integral:

[0128] (18)

[0129] We can get:

[0130] (19)

[0131] like , that is, the semi-major axis a The orbital inclination will change linearly under the action of continuous normal thrust:

[0132] (20)

[0133] like , semi-major axis a Under the action of tangential thrust, it changes according to formula (16):

[0134] (twenty one)

[0135] Substituting formula (21) into formula (19), we can obtain:

[0136] (twenty two)

[0137] Integrating both sides we get:

[0138] (twenty three)

[0139] is the integration constant, when t=0, , so:

[0140] (twenty four)

[0141] (25)

[0142] According to formula (25), the normal acceleration can be solved from the change in orbital inclination, the magnitude of the tangential acceleration and the time interval between the two arc segments.

[0143] Pre-identification of normal thrust parameters for changing the right ascension of the ascending node: If a continuous normal thrust is applied to change the right ascension of the ascending node ( In the opposite direction at latitude arguments of 0 and π), integrate the first equation of (14) over [0, π]:

[0144] (26)

[0145] We can get:

[0146] (27)

[0147] Ignore the small periodic changes in orbital inclination. , that is, the semi-major axis a The right ascension of the ascending node will change linearly under the action of continuous normal thrust:

[0148] (28)

[0149] like , semi-major axis a Under the action of tangential thrust, the derivation process is the same as that in Section 3.2.1, and we can obtain:

[0150] (29)

[0151] Equation (29) describes the change in the right ascension of the ascending node due to the combined effects of normal acceleration and tangential acceleration. In fact, the change in the right ascension of the ascending node is greatly affected by perturbations, of which the main influencing factor is the J2 perturbation 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 perturbation, the first-order long-term change rate of the right ascension of the ascending node is:

[0152] (30)

[0153] The calculation unit uses the normalized unit, and the corresponding central celestial body gravitational coefficient µ=GM=1. is the corresponding value of the average number of roots at the right ascension of the ascending node at time t, then:

[0154] (31)

[0155] According to formula (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 arc segments.

[0156] To facilitate understanding, the following is a simulation analysis of the scenario of changing the number of satellite orbit elements under the action of continuous small thrust.

[0157] Assume that a satellite has an initial orbital altitude of 380 km and a continuous low thrust is applied at 12:00 on January 25, 2024. The thrust parameters are set based on the "Starlink" V2.0mini: thruster specific impulse of 2500 s, electric thrust of 170 mN, and thrust acceleration of 2.26 × 10⁻⁴ m / s². The initial orbital elements at the satellite ignition point are shown in Table 2, and other simulation parameters are set as in Table 1.

[0158] Table 1

[0159]

[0160] Table 2

[0161]

[0162] Assume that radar equipment at different locations can acquire satellite measurement data. The radar equipment distribution information is shown in Table 3. Assume that the radar's range measurement accuracy is 30 m (1σ), the angle measurement accuracy is 0.01° (1σ), and the minimum elevation angle is 5°. Simulation results show the radar's measurement data (azimuth, elevation, and range) for the satellite.

[0163] Table 3

[0164]

[0165] The thrust acceleration result is obtained by inverse analysis of two radar observation data with a certain time interval. Assume that the initial time of the two radar observations is 、 The end times are 、 The observation time is 、 Depending on the positional relationship between the satellite and the observation station, the length of the visible arc at a single observation station can vary from 3 to 9 minutes. To ensure the accuracy of the semi-major axis solution, only observation arcs longer than 6 minutes are used for orbit determination.

[0166] In the continuous tangential thrust parameter pre-identification, the time of two observation arcs is shown in Table 4.

[0167] Table 4

[0168]

[0169] The two observations were simulated 100 times respectively, with different random observation errors imposed in each simulation. The semi-major axes of the two arcs were obtained by using the single arc segment orbit determination method described in Section 3 and recorded as 、 .

[0170] According to formula (17), and The tangential acceleration was inversely solved and solved based on two arc segments with an interval of 23 hours. The average relative error was only 3.2%.

[0171] In the pre-identification of the continuous normal thrust parameters with changing orbital inclination, the times of the two observation arcs are shown in Table 5.

[0172] Table 5

[0173]

[0174] The two observations were simulated 100 times respectively, with different random observation errors imposed in each simulation. The orbit inclinations of the two arcs were obtained using the single arc segment orbit determination method described in Section 3 and recorded as 、 .

[0175] According to formula (17), and The normal acceleration was inversely solved and solved based on two arc segments separated by 23 hours, with an average relative error of 28.5%.

[0176] In the pre-identification of the continuous normal thrust parameters with the change of the right ascension of the ascending node, the time of the two observation arcs is shown in Table 6.

[0177] Table 6

[0178]

[0179] The two observations were simulated 100 times respectively, with different random observation errors imposed in each simulation. The single arc segment orbit determination method described in Section 3 was used to obtain the right ascension of the ascending node of the two arc segments, which were recorded as 、 .

[0180] According to formula (17), and The normal acceleration was inversely solved and solved based on two arc segments separated by 23 hours, with an average relative error of 30.21%.

[0181] The method for determining the initial values ​​of spacecraft thrust parameters provided by the present invention provides a set of initial values ​​of orbital elements and thrust parameters that are relatively close to the true values ​​for a non-cooperative maneuvering spacecraft whose orbit and thrust parameters are unknown, under the condition of unknown initial values ​​of orbital elements and thrust information, which is conducive to better convergence of precise orbit determination and improves the accuracy and efficiency of precise orbit determination.

[0182] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any modifications or substitutions that can be easily conceived by a person skilled in the art within the technical scope disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be based on the scope of protection of the claims.

Claims

1. A method for determining the initial value of a spacecraft thrust parameter, characterized in that: include: Obtain satellite observation data; Selecting two arcs with observation arcs greater than a preset length from the satellite observation data, and performing orbit determination based on the two arcs and the least squares method of the J2 perturbation to obtain an orbit determination result; determining the thrust direction based on the orbit determination results by comparison; Acquiring 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; The orbit determination according to the least square method of the two arc segments and the J2 perturbation includes: The orbit is determined according to the following formula to obtain the epoch time Track status The best estimate of : Where, is the weight matrix, epoch time The orbital state, is the actual observed value, is the theoretical observation value; The method of 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 include: When there is tangential thrust, the thrust parameters are determined according to the following formula: : Where, It represents the cumulative change of satellite velocity caused by tangential thrust. Indicates the thrust action time, represents the initial semi-major axis of the satellite when the tangential thrust begins to act, represents the final semi-major axis after a period of tangential thrust, is the gravitational constant of the central celestial body, is the tangential thrust acceleration; When there is a normal thrust that changes the orbital inclination, the thrust parameter is determined according to the following formula : Where, is the normal thrust acceleration, is the current value of the orbital inclination, is the initial value of orbital inclination, is the average angular velocity of the orbit, is the thrust action time; When there is a normal thrust that changes the longitude of the ascending node, the thrust parameter is determined according to the following formula : Where, is the change in 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 J2.

2. The initial value determination method according to claim 1, characterized in that: Also includes: The orbit determination result is expressed according to the following formula: Where, is the first orbit determination result, is the second orbit determination result, is the first semi-major axis, is the second semi-major axis, The first eccentricity, is the second eccentricity, is the first orbital inclination, is the inclination of the second orbit, is the longitude of the first ascending node, is the longitude of the second ascending node, is the first pericenter argument, is the second pericenter argument, is the first true anomaly, is the second true anomaly.

3. The initial value determination method according to claim 2, characterized in that: The determining the thrust direction by comparing the orbit determination results includes: Obtain the preset thresholds for the semi-major axis, orbit inclination, and ascending node longitude; 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 value of the semi-major axis, it is determined that a tangential thrust exists; 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 value of the orbital inclination, determining that a normal thrust for changing the orbital inclination exists; 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 value of the ascending node longitude, it is determined that there is a normal thrust that changes the ascending node longitude.

4. The method for determining an initial value according to claim 1, wherein: Also includes: The cumulative change in satellite velocity due to the tangential thrust is determined by the following formula: : Where, is the semi-major axis of the orbit.

5. The method for determining an initial value according to claim 1, wherein: Also includes: The rate of change of the right ascension of the ascending node under the influence of the J2 term is determined by the following formula : Where, is the first harmonic term of the Earth's non-spherical gravitational belt, It is a semi-diameter.

Images