A method for automatic calculation of control quantity under customized retaining ring constraints

Abstract

The present disclosure is about a method for automatically calculating the control amount under the constraints of a customized holding loop. The method includes: determining the western boundary value of the spacecraft orbit holding loop, and calculating the initial value of the control amount of the semi-major axis of the spacecraft orbit; using the precise ephemeris data within the preset time as pseudo-observation data, improving the spacecraft orbit, and obtaining the corrected average atmospheric damping coefficient; establishing a trajectory loop reference sequence according to the spacecraft regression characteristics; predicting the trajectory of the spacecraft according to the corrected average atmospheric damping coefficient and the state motion equation of the spacecraft in the inertial reference system, and judging whether to enter the interval coarse screening calculation stage; if entering the interval coarse screening calculation stage, determining the coarse screening interval of the control amount of the semi-major axis of the spacecraft orbit; based on the coarse screening interval of the control amount, iteratively optimizing the optimal control amount of the semi-major axis of the spacecraft orbit. The present disclosure can automatically give a control strategy that satisfies the western boundary of the trajectory network to reach the specified range.

Description

Technical Field

[0001] The present disclosure relates to the field of aerospace measurement and control technology, and in particular to a method for automatically calculating a control quantity under a customized retaining ring constraint. Background Art

[0002] Low-orbit Earth observation spacecraft usually adopts a sun-synchronous regression orbit. Before launching into orbit, the periodic coverage days, number of coverage circles, sub-satellite target points, and ascending and descending orbit calculation methods are selected according to the payload performance. Due to the needs of the observation mission, such spacecraft have strict requirements on the sub-satellite point trajectory. In actual operation on orbit, due to the influence of perturbations such as atmospheric drag, the non-spherical shape of the earth, and solar light pressure, the trajectory of the spacecraft will deviate from the design value, and regular trajectory network maintenance control is required to keep the trajectory drift deviation within the required range.

[0003] In related technologies, the sub-satellite point trajectory maintenance method based on the semi-major axis offset of the orbit is often used to make the deviation of the sub-satellite point trajectory from the nominal trajectory change in a parabolic shape, so as to achieve the purpose of reducing the number of maintenance times. The selection of the offset is usually completed by the calculation personnel through repeated trial calculations. Usually, the orbit period variability is estimated based on the average root number, and the semi-major axis offset is calculated using an analytical expression. This method is difficult to separate the average root number calculation error, and there are differences with the actual position and speed of the spacecraft. Especially in the case of a narrow trajectory maintenance loop, it is difficult to determine the offset. In addition, the selection of the trajectory network control quantity is usually based on manual trial calculation. The semi-major axis offset is calculated by an analytical expression, combined with historical control records, and repeated trial calculations are performed until the requirements are met. This method requires high experience from the control personnel, because the analytical method roughly considers the influence of atmospheric drag perturbations, and does not comprehensively consider solar gravitational perturbations (orbital inclination change rate), space environment changes, etc., which deviates greatly from the actual trajectory situation, requires multiple trial calculations, and takes a long time. With the increasing number of spacecraft in orbit, the current mode that mainly relies on manual trial calculations has been severely challenged. It is necessary to study the automatic calculation method of trajectory network maintenance strategy to adapt to the management requirements of large numbers of spacecraft in the future.

[0004] Therefore, it is necessary to provide a new technical solution to improve one or more problems existing in the above solutions.

[0005] It should be noted that the information disclosed in the above background technology section is only used to enhance the understanding of the background of the present disclosure, and therefore may include information that does not constitute the prior art known to ordinary technicians in the field. Summary of the invention

[0006] The purpose of the present disclosure is to provide a customized method for automatically calculating the control amount under the constraint of a retaining ring, thereby overcoming one or more problems caused by the limitations and defects of the related art at least to a certain extent.

[0007] According to a first aspect of an embodiment of the present disclosure, a method for automatically calculating a control amount under a customized retaining ring constraint is provided, the method comprising:

[0008] Determine the western boundary value of the spacecraft orbit holding loop and calculate the initial value of the control variable of the semi-major axis of the spacecraft orbit;

[0009] The precise ephemeris data within the preset time is used as pseudo observation data to improve the spacecraft orbit and obtain the corrected average atmospheric damping coefficient;

[0010] Establishing a trajectory loop reference sequence according to the spacecraft regression characteristics; wherein the regression characteristics include regression days and regression circles, and the trajectory loop reference sequence includes: time and longitude;

[0011] Predicting the trajectory of the spacecraft according to the corrected atmospheric damping coefficient and the state motion equation of the spacecraft in the inertial reference system, and determining whether to enter the interval coarse screening calculation stage;

[0012] If the interval rough screening calculation phase is entered, a control amount rough screening interval of the semi-major axis of the spacecraft orbit is determined;

[0013] Based on the coarse screening interval of the control quantity, an iterative search is performed for the optimal control quantity of the semi-major axis of the spacecraft orbit.

[0014] In the embodiment of the present disclosure, the initial value of the control amount of the semi-major axis of the spacecraft orbit is:

[0015]

[0016] Among them, Δa is the initial value of the control variable of the semi-major axis of the spacecraft orbit, λ f is the western boundary value of the spacecraft orbit holding ring, a is the semi-major axis of the spacecraft orbit, is the average attenuation value of the semi-major axis of the spacecraft orbit, R e is the radius of the earth and π is the circumference of a circle.

[0017] In an embodiment of the present disclosure, the state motion equation of the spacecraft in the inertial reference frame is:

[0018]

[0019] Where X is the n-dimensional state vector to be improved, and the state vector is composed of the spacecraft position, velocity, atmospheric damping coefficient and dynamic model parameters, t is the motion time of the spacecraft, F is the n-dimensional nonlinear function of the state vector X, t 0 is the motion time of the spacecraft in the initial state, X 0 is the state vector at t 0 The value of the moment;

[0020] The actual observed quantity is:

[0021]

[0022] Among them, t i is a moment in the spacecraft's motion, X i t i The state vector of the spacecraft at time , Y i t i The actual observed value at the time, G(X i ,t i ) is the n-dimensional nonlinear function of the observation vector, From the initial state X 0 Starting from the state equation of motion, t is calculated i The calculated observation quantity at time, ε i is random noise.

[0023] In the embodiment of the present disclosure, the time and longitude in the trajectory loop reference reference sequence are respectively:

[0024]

[0025]

[0026] Where D is the return day, t 1 is the time for the spacecraft to pass the nominal target point in the first lap, λ 0 is the trajectory target point, λ 1 is the longitude of the spacecraft's first orbit over the nominal target point, J 2 constant, cosi is the orbital inclination, ΔL is the angular distance between two consecutive orbits in the regression period, is the angular distance between two adjacent tracks, and T is the intersection period of the track motion.

[0027] In an embodiment of the present disclosure, the step of predicting the trajectory of the spacecraft according to the atmospheric damping coefficient and the trajectory loop reference standard sequence and determining whether to enter the interval coarse screening calculation stage includes:

[0028] The position of the spacecraft is predicted according to the revised atmospheric damping coefficient, and the position is converted into the sub-satellite point, and the trajectory longitude deviation λ of the sub-satellite point at the start time is calculated according to the trajectory loop reference sequence. 1 If the trajectory longitude deviation is greater than the target trajectory ring longitude, the interval coarse screening calculation stage is entered.

[0029] In an embodiment of the present disclosure, if the interval rough screening calculation stage is entered, the step of determining the control amount interval rough screening of the semi-major axis of the spacecraft orbit includes:

[0030] Data fitting is performed with time as the independent variable and trajectory longitude deviation as the variable. If the trajectory longitude deviation satisfies the fitting function, the minimum value of the fitting function is calculated.

[0031] In an embodiment of the present disclosure, if the interval rough screening calculation stage is entered, the step of determining the control amount interval rough screening of the semi-major axis of the spacecraft orbit also includes:

[0032] If Δλ 1 -kλ f >0, then Δa 1 =2Δa 0 , iterate the prediction trajectory loop until (Δθ i+1 -kλ f )(Δλ i -kλ f )<0, the control quantity rough screening interval calculation ends, and the control quantity rough screening interval is [Δa i ,Δa i+1 ], where Δλ 1 is the minimum value of the fitting function, Δa 1 is the first iteration semi-major axis control variable, Δλ i is the minimum value of the fitting function calculated for the i-th time, Δλ i+1 is the minimum value of the fitting function calculated for the i+1th time, Δa i is the semi-major axis control value of the i-th iteration, Δa i+1 is the semi-major axis control value for the i+1th iteration, kλ f is the target trajectory ring longitude, and k is the customization coefficient.

[0033] In the embodiment of the present disclosure, if the interval rough screening calculation stage is entered, the step of determining the control amount interval rough screening of the semi-major axis of the spacecraft orbit also includes: if Δλ 1 -kλ f <0, then Δa 1 =0.5Δa 0 , iterate the prediction trajectory loop until (Δλ i+1 -kλ f )(Δλ i -kλ f )<0, the control quantity rough screening interval calculation ends, and the control quantity rough screening interval is [Δa i ,Δa i+1 ].

[0034] In an embodiment of the present disclosure, the step of iteratively optimizing the optimal control amount of the semi-major axis of the spacecraft orbit based on the control amount calculation interval includes:

[0035] The control amount coarse screening interval [Δa i ,Δa i+1] as the initial threshold, and [Δa i ,(Δa i+1 -Δa i ) / 2] is used as the control quantity to predict the trajectory and calculate the extreme point of the trajectory longitude deviation Δλ 1 ', adjust the control amount coarse screening interval according to the extreme point of the trajectory longitude deviation.

[0036] In an embodiment of the present disclosure, the step of iteratively optimizing the optimal control amount of the semi-major axis of the spacecraft orbit based on the control amount calculation interval also includes:

[0037] If Δλ 1 ′-kλ f >0, the control threshold is reduced to [(Δa i+1 -Δa i ) / 2,Δa i+1 ];

[0038] If Δλ 1 ′-kλ f <0, the control threshold is reduced to [Δa i ,(Δa i+1 -Δa i ) / 2], and given the error threshold ε, iterative calculation search until abs(Δλ i ′-kλ f )<ε, stop the calculation and get the optimal control amount Δa opm , where abs means absolute value, Δλ i ′ is the function extreme value calculated at the i-th iteration after entering the automatic optimization.

[0039] The technical solution provided by the embodiments of the present disclosure may have the following beneficial effects:

[0040] In one embodiment of the present disclosure, through the above method, the historical precise ephemeris of the spacecraft is used as pseudo-observation data to obtain actual observation quantities and calculated observation quantities, thereby obtaining observation residuals, and the optimal initial state vector of the spacecraft is calculated based on the observation residuals, and the corrected average atmospheric damping parameters are inverted to obtain the corrected average atmospheric damping parameters, and a trajectory loop reference benchmark sequence is established according to the regression characteristics of the spacecraft, and the corrected average atmospheric damping parameters and state motion equations are used in trajectory network forecasting to improve the forecast accuracy; and it is determined whether to enter the interval coarse screening calculation stage, and the control quantity coarse screening interval is obtained, and the optimal control quantity is iteratively optimized, and a control strategy that satisfies the western boundary of the trajectory network to reach the specified range is automatically given.

[0041] It is to be understood that the foregoing general description and the following detailed description are exemplary and explanatory only and are not restrictive of the present disclosure. BRIEF DESCRIPTION OF THE DRAWINGS

[0042] The accompanying drawings herein are incorporated into the specification and constitute a part of the specification, illustrate embodiments consistent with the present disclosure, and together with the specification are used to explain the principles of the present disclosure. Obviously, the accompanying drawings described below are only some embodiments of the present disclosure, and for ordinary technicians in this field, other accompanying drawings can be obtained based on these accompanying drawings without creative work.

[0043] Figure 1 A flowchart schematically illustrates the steps of a method for automatically calculating a control amount under a customized retaining ring constraint in an exemplary embodiment of the present disclosure;

[0044] Figure 2 Schematically illustrating the automatic iterative convergence process of the control quantities of three satellites in an exemplary embodiment of the present disclosure;

[0045] Figure 3 The automatic iterative convergence process of trajectory deviations of three satellites in an exemplary embodiment of the present disclosure is schematically illustrated. DETAILED DESCRIPTION

[0046] Example embodiments will now be described more fully with reference to the accompanying drawings. However, example embodiments can be implemented in a variety of forms and should not be construed as limited to the examples set forth herein; rather, these embodiments are provided so that the disclosure will be more comprehensive and complete and to fully convey the concepts of the example embodiments to those skilled in the art. The described features, structures, or characteristics may be combined in any suitable manner in one or more embodiments.

[0047] This example implementation provides a customized method for automatically calculating the control amount under the constraint of the retaining ring. Figure 1 As shown in , the method may include:

[0048] Step S101: Determine the western boundary value of the spacecraft orbit holding loop and calculate the initial value of the control amount of the semi-major axis of the spacecraft orbit;

[0049] Step S102: using the precise ephemeris data within a preset time as pseudo observation data, improving the spacecraft orbit, and obtaining a corrected average atmospheric damping coefficient;

[0050] Step S103: establishing a trajectory loop reference sequence according to the spacecraft regression characteristics; wherein the regression characteristics include regression days and regression circles, and the trajectory loop reference sequence includes: reference time and longitude;

[0051] Step S104: predicting the trajectory of the spacecraft according to the corrected average atmospheric damping coefficient and the state motion equation of the spacecraft in the inertial reference system, and determining whether to enter the interval coarse screening calculation stage;

[0052] Step S105: if the interval rough screening calculation phase is entered, a control amount rough screening interval of the semi-major axis of the spacecraft orbit is determined;

[0053] Step S106: Based on the coarse screening interval of the control amount, iteratively search for the optimal control amount of the semi-major axis of the spacecraft orbit.

[0054] Through the above method, the spacecraft's historical precise ephemeris is used as pseudo-observation data to obtain actual observations and calculated observations, thereby obtaining observation residuals. The optimal initial state vector of the spacecraft is calculated based on the observation residuals, and the corrected average atmospheric damping parameters are inverted to obtain the corrected average atmospheric damping parameters. The trajectory loop reference sequence is established according to the spacecraft regression characteristics, and the corrected average atmospheric damping parameters and state motion equations are used in trajectory network forecasting to improve the forecast accuracy. It is also determined whether to enter the interval coarse screening calculation stage, and the control quantity coarse screening interval is obtained. The optimal control quantity is iteratively optimized, and a control strategy that satisfies the western boundary of the trajectory network to reach the specified range is automatically given.

[0055] Next, we will refer to Figure 1 Each step of the above method in this example implementation is described in more detail.

[0056] In step S101, the western boundary value of the spacecraft orbit holding loop is determined, and the initial value of the control amount of the semi-major axis of the spacecraft orbit is calculated. Specifically, in order to control the drift of the ground track, the adjustment amount of the control amount of the semi-major axis of the spacecraft orbit and the adjustment time interval are related to the allowable drift range. Set the western boundary value of the spacecraft orbit holding loop to control the trajectory within the range of the western boundary values ​​of the east and west holding loops of the nominal trajectory, then the total allowable drift is twice the western boundary value of the holding loop. In the analytical model based on the square root number, the initial value of the control amount of the semi-major axis of the spacecraft orbit is:

[0057]

[0058] Among them, Δa is the initial value of the control variable of the semi-major axis of the spacecraft orbit, λ f is the western boundary value of the spacecraft orbit holding ring, a is the semi-major axis of the spacecraft orbit, is the average attenuation value of the semi-major axis of the spacecraft orbit, R e is the radius of the earth, and π is the pi. Specifically, the initial value of the control amount of the semi-major axis of the spacecraft orbit is also the adjustment amount of the semi-major axis of the spacecraft orbit each time.

[0059] In step S102, the precise ephemeris data within the preset time is used as pseudo observation data, and the orbit of the spacecraft is improved to obtain the corrected average atmospheric damping coefficient. Specifically, the precise ephemeris data within the preset time is used as pseudo observation data to obtain the actual observation amount and the calculated observation amount of the spacecraft, and the observation residual is obtained based on the actual observation amount and the calculated observation amount, and the optimal initial state vector of the spacecraft is calculated based on the observation residual, and the corrected average atmospheric damping parameters are obtained by inversion; wherein, the optimal state vector of the spacecraft is composed of the spacecraft position, velocity, corrected average atmospheric damping coefficient and dynamic model parameters. Among them, the precise ephemeris data of the spacecraft within a period of time is used as pseudo observation data, and the state vector of the spacecraft is determined by the batch least squares method. The components of the state vector include the position, velocity, atmospheric damping coefficient, mechanical model and measurement model parameters of the spacecraft. The dynamic equation of the spacecraft is integrated from a reference epoch to each observation time to generate the predicted observation amount. Due to the difference between the predicted observation amount and the actual observation amount, the observation residual is generated. In the least squares sense, the components of the spacecraft's state vector are adjusted to minimize the observation residuals.

[0060] In one embodiment, the state motion equation of the spacecraft in the inertial reference frame is:

[0061]

[0062] Where X is the n-dimensional state vector to be improved, which is composed of the spacecraft position, velocity, atmospheric damping coefficient and some dynamic model parameters. F is the n-dimensional nonlinear function of the state vector X. 0 is the state vector at t 0 The value of the moment.

[0063] The observed quantity is a nonlinear function of the state quantity, which can be expressed as follows:

[0064]

[0065] Where Y i t i The actual observed value at time (denoted as O), From the initial state X 0 Starting from the perturbation model, t is calculated i The calculated observation quantity (denoted as C) at the time of ε is random noise. In fact, due to the initial state X 0 Inaccurate, the calculated observations do not match the actual observations, and there will be observation residuals (OC) between the actual observations and the calculated observations. Through the pseudo-observation data of the spacecraft in a certain arc segment (precise ephemeris of the spacecraft), the weighted sum of squares of the observation residuals is minimized as an indicator to obtain the optimal initial state quantity, and the corrected atmospheric damping parameters are inverted and used as the input parameters of the trajectory network forecast.

[0066] In step S103, a trajectory loop reference sequence is established according to the spacecraft regression characteristics; wherein the regression characteristics include regression days and regression circles, and the trajectory loop reference sequence includes: time and longitude.

[0067] In one embodiment, the time and longitude in the trajectory loop reference sequence are:

[0068]

[0069]

[0070] Where D is the return day, t 1 ′ is the time when the spacecraft passes the nominal target point for the first time, t i ′ is the time it takes for the spacecraft to pass the nominal target point in the i-th circle, λ 0 is the trajectory target point, λ 1 is the longitude of the spacecraft's first orbit over the nominal target point, J 2 is a constant, cosi is the orbital inclination, ΔL is the angular distance between two consecutive orbits in the regression period, is the angular distance between two adjacent tracks, and T is the intersection period of the track motion. It should be noted that when the track target point λ 0 After the regression characteristics are determined, the reference longitude of each orbit of the sun-synchronous regression orbit spacecraft can be determined. Assuming that the drift to the west is negative and the drift to the east is positive relative to the nominal trajectory, λ 1 is the longitude of the spacecraft's first orbit over the nominal target point, t 1 ′ is the time when the spacecraft passes the nominal target point in the first circle, and the nominal time and longitude of the target point after the i-th circle are obtained. The longitude of the sub-satellite point of each circle of the spacecraft must fall between the two nominal reference longitudes. It is difficult to determine the actual trajectory deviation of the spacecraft only from the orbit calculation data of the current circle. In engineering practice, the "track jumping" phenomenon is easy to occur, and the wrong trajectory deviation result is calculated. Take the establishment of the nominal trajectory loop reference sequence (t i ′,λ i ), the nominal trajectory reference longitude sequence is predicted during the calculation process, the sub-satellite point trajectory automatically matches the nominal reference longitude through time, the deviation of the trajectory loop from the reference datum is determined, and track hopping is avoided.

[0071] It should also be noted that the angular distance ΔL between two consecutive tracks and the angular distance between two adjacent tracks in the trajectory loop reference sequence within the regression period are It is obtained by the following formulas (5), (6) and (7).

[0072] In one embodiment, the sub-satellite point trajectory of the sun-synchronous regression orbit repeats periodically. Since the lateral movement of the sub-satellite point trajectory is a combination of the earth's rotation, the orbital pitch line precession and the orbital movement of the spacecraft, within one orbital period, the angular distance ΔL between two consecutive orbits of the sub-satellite point within the regression period is:

[0073]

[0074] Where T is the intersection period of the orbital motion, ω e is the angular velocity of the Earth's rotation, is the orbital precession rate.

[0075] Assuming that on D day, the spacecraft runs N orbits, the westward displacement of the sub-satellite point satisfies:

[0076]

[0077] Assume that the angular distance between two adjacent tracks is for

[0078]

[0079] From formulas (6) and (7), we can see that the angular distance ΔL between two consecutive tracks in the regression period is equal to the angular distance between two adjacent tracks. The relationship is,

[0080]

[0081] In step S104, the trajectory of the spacecraft is predicted based on the corrected average atmospheric damping coefficient and the state motion equation of the spacecraft in the inertial reference frame, and it is determined whether to enter the interval coarse screening calculation stage.

[0082] In one embodiment, the position of the spacecraft is predicted according to the revised atmospheric damping coefficient, and the position is converted into the sub-satellite point, and the trajectory longitude deviation of the sub-satellite point at the start time is calculated according to the trajectory loop reference sequence. If the trajectory longitude deviation is greater than the target trajectory loop longitude, the interval coarse screening calculation stage is entered. Specifically, based on the spacecraft keeping the loop western boundary value λ f , and the boundary of the spacecraft keeping ring is [-λ f ,λ f ], the target trajectory ring longitude is kλ f , k is the customization coefficient of the holding loop, and the trajectory longitude deviation λ of the sub-satellite point at the start-up time is calculated according to the trajectory loop reference sequence 1 , if λ 1 >kλ f , it is judged that the semi-major axis of the spacecraft has entered the interval coarse screening calculation stage. Otherwise, it is prompted that the trajectory loop has exceeded the expected western boundary and the calculation is terminated.

[0083] In step S105, if the interval coarse screening calculation stage is entered, the control amount coarse screening interval of the semi-major axis of the spacecraft orbit is determined.

[0084] In one embodiment, data fitting is performed with time as the independent variable and the trajectory longitude deviation as the variable. If the trajectory longitude deviation satisfies the fitting function, the minimum value of the fitting function is calculated. Specifically, if the trajectory deviation satisfies the linear fitting f(x)=mx+n, the minimum value f(x) is calculated. min =Δλ 1 ; If the trajectory deviation satisfies the quadratic fitting f(x) = mx 2 +nx+l, calculate the extreme point f(x) min =Δλ 1 , where m, n, and l are constants.

[0085] In one embodiment, if Δλ 1 -kλ f >0, then Δa 1 =2Δa 0 , iterate the prediction trajectory loop until (Δλ i+1 -kλ f )(Δλ i -kλ f )<0, the control quantity rough screening interval calculation ends, and the control quantity rough screening interval is [Δa i ,Δa i+1 ], where Δλ 1 is the minimum value of the fitting function, Δa 1 is the first iteration semi-major axis control variable, Δλ i is the minimum value of the fitting function calculated for the i-th time, Δλ i+1 is the minimum value of the fitting function calculated for the i+1th time, Δa i is the semi-major axis control value of the i-th iteration, Δa i+1 is the semi-major axis control value for the i+1th iteration, kλ f is the target trajectory ring longitude, and k is the ring-keeping customization coefficient. Specifically, if the trajectory longitude deviation satisfies the minimum value of the fitting function Δλ 1 Greater than the target trajectory ring longitude kλ f , then take 2 times the initial value of the control amount as the control amount, and perform the iterative prediction trajectory loop until (Δλ i+1 -kλ f )(Δλ i -kλ f )<0, the control amount rough screening interval calculation is completed, and the control amount rough screening interval is determined to be [Δa i ,Δa i+1 ]

[0086] In one embodiment, if Δλ 1 -kλf <0, then Δa 1 =0.5Δa 0 , iterate the prediction trajectory loop until (Δλ i+1 -kλ f )(Δλ i -kλ f )<0, the control quantity rough screening interval calculation ends, and the control quantity rough screening interval is [Δa i ,Δa i+1 Specifically, if the trajectory longitude deviation satisfies the minimum value of the fitting function Δλ 1 Less than the target trajectory ring longitude kλ f , then take half of the initial value of the control quantity as the control quantity, and perform the iterative prediction trajectory loop until (Δλ i+1 -kλ f )(Δλ i -kλ f )<0, the control amount rough screening interval calculation is completed, and the control amount rough screening interval is determined to be [Δa i ,Δa i+1 ].

[0087] In step S106, based on the coarse screening interval of the control amount, an iterative search is performed for the optimal control amount of the semi-major axis of the spacecraft orbit.

[0088] In one embodiment, the control amount coarse screening interval [Δa i ,Δa i+1 ] as the initial threshold, and [Δa i ,(Δa i+1 -Δa i ) / 2] is used as the control quantity to predict the trajectory and calculate the extreme point of the trajectory longitude deviation Δλ 1 ′, adjust the control quantity coarse screening interval according to the extreme point of the trajectory longitude deviation. Specifically, based on the control quantity coarse screening interval, when iteratively optimizing the optimal control quantity of the semi-major axis of the spacecraft orbit, by judging the positive and negative signs of the function values ​​of the middle values ​​of the control quantity coarse screening interval, the finite interval is gradually reduced in half until the interval is reduced to within the allowable error range, and then the midpoint of the interval is taken as the approximate value. The control quantity coarse screening interval [Δa i ,Δa i+1 ] as the initial threshold, and [Δa i ,(Δa i+1 -Δa i ) / 2] is used as the control quantity to predict the trajectory and calculate the extreme point of the trajectory longitude deviation Δλ 1 ′, and adjust the control amount coarse screening interval according to the size of the extreme point of the trajectory longitude deviation and the target trajectory longitude ring.

[0089] In one embodiment, if Δλ 1 ′-kλf >0, the control threshold is reduced to [(Δa i+1 -Δa i ) / 2,Δa i+1 ];

[0090] If Δλ 1 ′-kλ f <0, the control threshold is reduced to [Δa i ,(Δa i+1 -Δa i ) / 2], and given the error threshold ε, iterative calculation search until abs(Δλ i ′-kλ f )<ε, stop the calculation and get the optimal control amount Δa opm , where abs means absolute value, Δλ i ′ is the function extreme value calculated at the i-th iteration after entering the automatic optimization. Specifically, if the trajectory longitude deviation extreme point Δλ 1 ′ is greater than the target trajectory longitude ring kλ f , then the control threshold is reduced to [(Δa i+1 -Δa i ) / 2,Δa i+1 ]. If the trajectory longitude deviation extreme point Δλ 1 ′ is less than the target trajectory longitude ring kλ f , then the control threshold is reduced to [Δa i ,(Δa i+1 -Δa i ) / 2], and based on the given error threshold ε, perform iterative calculation and search until abs(Δλ i ′-kλ f )<ε, stop the calculation and get the optimal control amount Δa opm .

[0091] Through the above method, the spacecraft's historical precise ephemeris is used as pseudo-observation data to obtain actual observations and calculated observations, thereby obtaining observation residuals. The optimal initial state vector of the spacecraft is calculated based on the observation residuals, and the corrected average atmospheric damping parameters are inverted to obtain the corrected average atmospheric damping parameters. The trajectory loop reference sequence is established according to the spacecraft regression characteristics, and the corrected average atmospheric damping parameters and state motion equations are used in trajectory network forecasting to improve the forecast accuracy. It is also determined whether to enter the interval coarse screening calculation stage, and the control quantity coarse screening interval is obtained. The optimal control quantity is iteratively optimized, and a control strategy that satisfies the western boundary of the trajectory network to reach the specified range is automatically given.

[0092] The present disclosure is further described below through the following examples.

[0093] Assuming three satellites with different trajectory ring constraints, the western boundary trajectory deviation is set to maintain 40% of the ring constraint, as shown in Table 1. The control optimization algorithm given in this paper is used to calculate the control quantity and the implementation process is given.

[0094] Table 1 Orbital altitudes of the three satellites and the constraints of the retaining ring

[0095]

[0096] Step 1: Calculation of control objectives and initial values

[0097] (1) According to the three satellites’ ring constraints and customized coefficients in Table 1, the automatically given satellite western boundary value kλ can be calculated. f They are -1.12 km, -4 km, and -8 km respectively, with a convergence error of 4%.

[0098] (2) Calculate the average attenuation value of the semi-major axis of the satellite orbit Using Formula 1, according to the semi-major axes of the three satellite orbits and the constraints of the retaining ring, the initial values ​​of the iterative calculation of the control quantities of the three satellites are obtained, which are 72.76 meters, 142.86 meters, and 502.11 meters respectively.

[0099] Step 2: Determine the average atmospheric damping coefficient

[0100] (1) Download satellite precise orbits and generate GPS pseudo-observation data;

[0101] (2) Using the first point ephemeris data of the precise orbit as the initial orbit, the orbit is improved and the average atmospheric damping coefficient is solved. The average atmospheric damping coefficients of the three satellites are 2.04, 1.89, and 2.37 respectively.

[0102] Step 3: Create a reference sequence for the trajectory loop

[0103] (1) Based on the return days and return circles of the three satellites, the sub-satellite point interval between two adjacent orbits of the satellite is calculated using Formula 6 The sub-satellite point interval between two adjacent orbits of three satellites They are 0.01038, 0.007 and 0.01038 respectively.

[0104] (2) Based on the Kepler orbit of the satellite, the reference longitude sequence of the nominal trajectory of the three satellites is established using formulas 3 and 4 (t i ,λ i ).

[0105] Step 4: Rough screening and calculation of control quantity interval

[0106] (1) Interval rough screening judgment. Use the precise orbit prediction model to predict the trajectory of the spacecraft and calculate the longitude of the sub-satellite trajectory deviation at the start-up time of the three satellites as λ1 , respectively 2484 meters, 9125 meters, 19362 meters, if λ 1 >kλ f , enter the interval coarse screening calculation stage.

[0107] (2) Interval rough screening calculation. Take the initial value Δa of the three satellite control quantities 0 The trajectory loop is predicted within a specified time period, and the data is fitted with time as the independent variable and the trajectory longitude deviation as the variable. The extreme point Δλ 1 They are 2172 meters, 6801 meters, and 19022 meters respectively. 1 -kλ f > 0, then continue with Δa 1 =2Δa 0 Predict the trajectory loop within the specified time period, the extreme point Δλ 2 They are -6518 meters, -33462 meters, and -25657 meters respectively. i+1 -kλ f )(λ i -kλ f )<0, the rough screening intervals of the three satellite control quantities can be determined as [Δa 0 ,Δa 1 ].

[0108] Step 5: Iterative optimization calculation of control quantity

[0109] Control volume rough screening interval [Δa i ,Δa i+1 ] is the initial threshold, and (Δa i+1 -Δa i ) / 2 is the control quantity for trajectory prediction and calculation of the trajectory deviation extreme point Δλ 1 ′, if Δλ 1 ′-kλ f >0, the control threshold is reduced to [(Δa i+1 -Δa i ) / 2,Δa i+1 ], if Δλ 1 ′-kλ f <0, the control threshold is reduced to [Δa i ,(Δa i+1 -Δa i ) / 2], given the error threshold ε, iterative calculation search until abs(λ i ′-kλ f )<ε, stop the calculation and get the optimal control amount Δa opm The specific iterative convergence process of the three satellites is shown in Table 2-Table 4, and the automatic iterative convergence process of the control quantities of the three satellites is shown in Figure 2The automatic iterative convergence process of the trajectory deviations of the three satellites is shown in Figure 3 .

[0110] Table 2. SAT-1 iterative convergence process

[0111]

[0112] Table 3. SAT-2 iterative convergence process

[0113] Iterations schedule Control volume (m) Deviation extreme value (m) 1. Interval rough screening judgment 0 9125 2. Interval rough screening calculation 142.86 6801 3. Interval rough screening calculation 285.72 -33462 4. Control optimization module 214.29 -6112 5. Control optimization module 178.58 1843 6. Control optimization module 196.44 -1802 7. Control optimization module 205.36 -3842 8. Control optimization module 209.82 -4974 9. Control optimization module 207.59 -4974 10. Termination threshold 206.48 -4131

[0114] Table 4. SAT-3 iterative convergence process

[0115]

[0116]

[0117] Those skilled in the art will readily appreciate other embodiments of the present disclosure after considering the specification and practicing the invention disclosed herein. This application is intended to cover any modification, use or adaptation of the present disclosure, which follows the general principles of the present disclosure and includes common knowledge or customary techniques in the art that are not disclosed in the present disclosure. The specification and examples are intended to be exemplary only, and the true scope and spirit of the present disclosure are indicated by the appended claims.

Claims

1. A customized method for automatically calculating the control quantity under the constraint of the retaining ring, characterized in that: The method includes: Determine the western boundary value of the spacecraft orbit holding loop and calculate the initial value of the control variable of the semi-major axis of the spacecraft orbit; The precise ephemeris data within the preset time is used as pseudo observation data to improve the spacecraft orbit and obtain the corrected average atmospheric damping coefficient; Establishing a trajectory loop reference sequence according to the spacecraft regression characteristics; wherein the regression characteristics include regression days and regression circles, and the trajectory loop reference sequence includes: time and longitude; Predicting the trajectory of the spacecraft according to the corrected average atmospheric damping coefficient and the state motion equation of the spacecraft in the inertial reference system, and determining whether to enter the interval coarse screening calculation stage; If the interval rough screening calculation phase is entered, a control amount rough screening interval of the semi-major axis of the spacecraft orbit is determined; Based on the coarse screening interval of the control quantity, an iterative search is performed for the optimal control quantity of the semi-major axis of the spacecraft orbit.

2. According to the automatic calculation of the control quantity under the customized retaining ring constraint of claim 1, it is characterized in that: The initial value of the control quantity of the semi-major axis of the spacecraft orbit is: Among them, Δa0 is the initial value of the control variable of the semi-major axis of the spacecraft orbit, λ f is the western boundary value of the spacecraft orbit holding ring, a is the semi-major axis of the spacecraft orbit, is the average attenuation value of the semi-major axis of the spacecraft orbit, R e is the radius of the earth and π is the circumference of a circle.

3. The method for automatically calculating the control amount under the customized retaining ring constraint according to claim 2 is characterized in that: The state motion equation of the spacecraft in the inertial reference frame is: Where X is the n-dimensional state vector to be improved, and the state vector is composed of the spacecraft position, velocity, average atmospheric damping coefficient and dynamic model parameters, t is the motion time of the spacecraft, F is the n-dimensional nonlinear function of the state vector X, t0 is the motion time of the spacecraft initial state X0 is the value of the state vector at time t0; The actual observed quantity is: Among them, t i is a moment in the spacecraft's motion, X i t i The state vector of the spacecraft at time , Y i t i The actual observed value at the time, G(X i ,t i ) is the n-dimensional nonlinear function of the observation vector, is t calculated under the state motion equation starting from the initial state X0 i The calculated observation quantity at time, ε i is random noise.

4. The method for automatically calculating the control amount under the customized retaining ring constraint according to claim 3 is characterized in that: The time and longitude in the trajectory loop reference sequence are: Where D is the return day, t1′ is the time when the spacecraft passes the nominal target point for the first time, and t i ′ is the time when the spacecraft passes the nominal target point in the i-th circle, λ0 is the trajectory target point, λ1 is the longitude of the spacecraft passing the nominal target point in the first circle, J2 is the constant, cosi is the orbit inclination, ΔL is the angular distance between two consecutive tracks in the regression period, is the angular distance between two adjacent tracks, and T is the intersection period of the track motion.

5. The method for automatically calculating the control amount under the customized retaining ring constraint according to claim 4 is characterized in that: The step of predicting the trajectory of the spacecraft according to the corrected average atmospheric damping coefficient and the state motion equation of the spacecraft in the inertial reference system and determining whether to enter the interval coarse screening calculation stage includes: The position of the spacecraft is predicted according to the revised average atmospheric damping coefficient, and the position is converted into the sub-satellite point. The trajectory longitude deviation of the sub-satellite point at the start-up time is calculated according to the trajectory ring reference benchmark sequence. If the trajectory longitude deviation is greater than the target trajectory ring longitude, the interval coarse screening calculation stage is entered.

6. The method for automatically calculating the control amount under the customized retaining ring constraint according to claim 5 is characterized in that: If the interval rough screening calculation stage is entered, the step of determining the control amount interval rough screening of the semi-major axis of the spacecraft orbit includes: Data fitting is performed with time as the independent variable and trajectory longitude deviation as the variable. If the trajectory longitude deviation satisfies the fitting function, the minimum value of the fitting function is calculated.

7. The method for automatically calculating the control amount under the customized retaining ring constraint according to claim 6 is characterized in that: If the interval rough screening calculation stage is entered, the step of determining the control amount interval rough screening of the semi-major axis of the spacecraft orbit also includes: If Δλ1-kλ f >0, then Δa1=2Δa0, iterate the prediction trajectory loop until (Δλ i+1 -kλ f )(Δλ i -kλ f )<0, the control quantity rough screening interval calculation ends, and the control quantity rough screening interval is [Δa i ,Δa i+1 ], where Δλ1 is the minimum value of the fitting function, Δa1 is the semi-major axis control value of the first iteration, and Δλ i is the minimum value of the fitting function calculated for the i-th time, Δλ i+1 is the minimum value of the fitting function calculated for the i+1th time, Δa i is the semi-major axis control value of the i-th iteration, Δa i+1 is the semi-major axis control value for the i+1th iteration, kλ f is the target trajectory ring longitude, and k is the ring-keeping customization coefficient.

8. The method for automatically calculating the control amount under the customized retaining ring constraint according to claim 7 is characterized in that: If the interval rough screening calculation stage is entered, the step of determining the control amount interval rough screening of the semi-major axis of the spacecraft orbit also includes: If Δλ1-kλ f <0, then Δa1=0.5Δa0, iterate the prediction trajectory loop until (Δλ i+1 -kλ f )(Δλ i -kλ f )<0, the control quantity rough screening interval calculation ends, and the control quantity rough screening interval is [Δa i ,Δa i+1 ].

9. The method for automatically calculating the control amount under the customized retaining ring constraint according to claim 8 is characterized in that: The step of iteratively optimizing the optimal control amount of the semi-major axis of the spacecraft orbit based on the control amount calculation interval includes: The control amount coarse screening interval [Δa i ,Δa i+1 ] as the initial threshold, and [Δa i ,(Δa i+1 -Δa i ) / 2] is used to predict the trajectory of the control amount, calculate the extreme point of the trajectory longitude deviation Δλ1′, and adjust the coarse screening interval of the control amount according to the extreme point of the trajectory longitude deviation.

10. The method for automatically calculating the control amount under the customized retaining ring constraint according to claim 9 is characterized in that: The step of iteratively optimizing the optimal control amount of the semi-major axis of the spacecraft orbit based on the control amount calculation interval also includes: If Δλ1′-kλ f >0, the control threshold is reduced to [(Δa i+1 -Δa i ) / 2,Δa i+1 ]; If Δλ1′-kλ f <0, then reduce the control threshold to [Δa i ,(Δa i+1 -Δa i ) / 2], and given the error threshold ε, iterative calculation search until abs(Δλ i ′-kλ f )<ε, stop the calculation and get the optimal control amount Δa opm , where abs means absolute value, Δλ i ′ is the function extreme value calculated at the i-th iteration after entering the automatic optimization.

Images