A high-precision phase-keeping method for low-orbit constellation systems considering long-period terms

Abstract

The present invention discloses a high-precision phase holding method for a low-orbit constellation system considering a long-period term, including: determining the relationship between the change of the satellite reference flat phase angle and time; sampling the position and speed of the satellite once at a fixed time interval, converting it into an actual flat phase angle, and calculating the long-period term of the flat phase angle at this time; calculating the flat phase angle deviation, using the long-period term to make corrections, and obtaining the flat phase angle deviation coefficient; predicting the flat phase angle deviation at the next moment, determining the moment when the flat phase angle deviation exceeds a threshold; determining the initial solution of the change in the semi-major axis of the phase-keeping orbit maneuver; predicting the flat phase angle deviation of the next phase-keeping cycle, using the long-period term to make corrections, and determining the corrected solution of the change in the semi-major axis of the phase-keeping orbit maneuver; applying a speed increment when the flat phase angle deviation reaches a second threshold, and repeating the above steps until the phase holding is completed. The high-precision phase holding of the low-orbit constellation system is achieved in a perturbation environment including the non-spherical shape of the earth, atmospheric drag, the gravity of the three bodies of the sun and the moon, and the solar light pressure.

Description

Technical Field

[0001] The present invention relates to a high-precision phase keeping method for a low-orbit constellation system, and more specifically, to a high-precision phase keeping method for a low-orbit constellation system taking into account long-period terms. Background Art

[0002] Use low-orbit constellations to provide satellite Internet services covering the entire world. As a huge space system, a stable configuration is the basis for the normal functioning of the constellation. Due to differences in the orbital deviations and perturbations of each satellite, the constellation configuration tends to diverge, thus affecting the coverage characteristics and mission realization of the constellation. In addition, the constellation system will face the threat of collisions between satellites and between satellites and space debris. In order for the low-orbit constellation to meet the coverage characteristics requirements and space orbit safety requirements, the satellite phase of the low-orbit constellation needs to be kept within a certain small range (for example, within ±0.1° of the reference phase).

[0003] In the research of high-precision phase keeping of low-orbit constellations, the existing technology, based on the assumption that the first-order rate of change of the semi-major axis is mainly caused by atmospheric drag and is approximately constant, establishes a model in which the deviation of the flat phase angle is a quadratic function over time, calculates the change in the phase-keeping orbital altitude, and achieves the phase keeping of the low-orbit constellation system within ±0.1°. However, this method does not consider the influence of the long-period term of the flat phase angle caused by perturbations such as the gravity of the sun, the moon and the three bodies. Therefore, although this method can keep the phase within the range of ±0.1°, it does not fully utilize the phase keeping range, and the phase keeping orbit maneuvers are relatively frequent.

[0004] Therefore, it is necessary to provide a high-precision phase-keeping method for a low-orbit constellation system that takes into account long-period terms, and use the long-period term of the flat phase angle caused by the perturbations of the earth's non-sphericity, atmospheric drag, the gravitational force of the sun and the moon, and the solar pressure as a correction term for the flat orbit value, separate it from the flat phase angle, and then fit the flat phase angle deviation coefficient to accurately predict the flat phase angle deviation of the current cycle and the next cycle, and achieve precise control of the second threshold of the flat phase angle deviation by iteratively calculating the moment when the flat phase angle deviation reaches the second threshold; and achieve precise control of the first threshold of the flat phase angle deviation in the next phase-keeping cycle by iteratively calculating the size of the change in the semi-major axis. The method proposed in the present invention can achieve high-precision phase keeping for a low-orbit constellation system in a complex perturbation environment including the earth's non-sphericity, atmospheric drag, the gravitational force of the sun and the moon, and the solar pressure. Summary of the invention

[0005] The object of the present invention is to provide a high-precision phase keeping method for a low-orbit constellation system taking into account long-period terms.

[0006] In order to achieve the above object, the present invention adopts the following technical solutions:

[0007] The present invention provides a high-precision phase keeping method for a low-orbit constellation system considering a long-period term, the method comprising the steps of:

[0008] S1. Determine the relationship between the satellite reference flat phase angle and time;

[0009] S2, sampling the position and velocity of the satellite at fixed intervals, converting them into actual average phase angles, and calculating the average phase angle long-period term at the corresponding moment;

[0010] S3, subtracting the actual flat phase angle from the reference flat phase angle to calculate the flat phase angle deviation, and using the flat phase angle long period term for correction to obtain the flat phase angle deviation coefficient;

[0011] S4, predicting the flat phase angle deviation at the next sampling moment, and if it exceeds the second threshold, calculating the moment when the flat phase angle deviation exceeds the second threshold;

[0012] S5, determining an initial solution for the change in the semi-major axis of the phase-keeping orbit maneuver according to the moment when the second threshold value is exceeded and the flat phase angle deviation coefficient;

[0013] S6. Predict the flat phase angle deviation of the next phase holding period according to the flat phase angle deviation coefficient of the current phase holding period, and use the long period term for correction to determine the correction solution of the change of the semi-major axis of the phase holding orbit maneuver;

[0014] S7, applying a velocity increment when the flat phase angle deviation value reaches the second threshold, taking the orbital maneuver end time as the initial time of the next control cycle, and repeating S2 to S7 until the phase holding of the entire phase holding mission cycle of the satellite is completed.

[0015] Preferably, the first threshold is the value of the left endpoint of the flat phase angle deviation interval, and the second threshold is the value of the right endpoint of the flat phase angle deviation interval.

[0016] Preferably, the step S1 comprises:

[0017] Determine the rate of change of the reference flat phase angle over time based on the gravitational field at the center of the Earth, the non-spherical perturbation of the Earth, and the gravitational perturbation of the Sun and the Moon;

[0018] The reference flat phase angle of the satellite at a certain moment is determined according to the reference flat phase angle of the satellite at the initial moment of the satellite phase keeping mission cycle and the rate of change of the reference flat phase angle over time.

[0019] Preferably, the step S2 comprises:

[0020] Substituting the average orbital elements into the formula of the long-period term caused by the non-sphericity of the earth, atmospheric drag, the gravity of the sun and the moon, and the solar light pressure perturbation, the flat phase angle long-period term is obtained.

[0021] Preferably, step S3 comprises:

[0022] The flat phase angle long-period term is used for correction to obtain a flat phase angle deviation correction value, and the flat phase angle deviation coefficient is calculated by fitting according to the variation trend of the flat phase angle deviation correction value over time.

[0023] Preferably, the flat phase angle deviation coefficient is a quadratic polynomial coefficient of the flat phase angle deviation correction value changing with time, including the flat phase angle deviation at the initial moment, the first-order term of the flat phase angle drift caused by the semi-major axis deviation at the initial moment, and the second-order term of the flat phase angle drift caused by the linear change of the semi-major axis with time.

[0024] Preferably, the step S4 further comprises:

[0025] If the flat phase angle deviation forecast value does not exceed the second threshold, steps S2 to S3 are repeated to update the flat phase angle deviation coefficient until the flat phase angle deviation at this moment exceeds the second threshold;

[0026] If the flat phase angle deviation forecast value exceeds the second threshold, a free drift duration less than one sampling period is introduced to calculate the time when the flat phase angle deviation exceeds the second threshold.

[0027] Preferably, the step S5 comprises:

[0028] According to the flat phase angle deviation coefficient, calculating the deviation between the actual semi-major axis and the reference semi-major axis at the initial moment and the first-order linear change rate of the orbital semi-major axis of the satellite from the initial moment to the moment when the flat phase angle deviation exceeds the second threshold, thereby determining the deviation between the actual semi-major axis and the reference semi-major axis at the moment when the flat phase angle deviation exceeds the second threshold;

[0029] According to the flat phase angle deviation interval and the second-order term of the flat phase angle drift, the duration between two orbital altitude increases of the next phase is estimated;

[0030] An initial solution for the change in the semi-major axis of the orbital maneuver is determined based on the semi-major axis drift between two orbital altitude increases and the deviation between the actual semi-major axis and the reference semi-major axis when the flat phase angle deviation exceeds a second threshold.

[0031] Preferably, step S6 comprises:

[0032] Predict the sun, moon and satellite orbital elements of the next phase holding period, calculate the corresponding long-period term of the phase angle deviation, obtain the mean phase angle deviation of the next phase holding period, and determine a corrected solution for the change in the semi-major axis of the phase holding orbit maneuver that makes the minimum value of the mean phase angle deviation of the next phase holding period just reach the first threshold.

[0033] Preferably, the orbital elements of the sun, moon and satellite are different in different phase holding periods, and the long-period terms of the mean phase angle caused by the non-sphericity of the earth, atmospheric drag, gravity of the sun and moon, and solar light pressure perturbations are also different.

[0034] The beneficial effects of the present invention are as follows:

[0035] 1. The high-precision phase keeping method for a low-orbit constellation system considering long-period terms of the present invention separates the long-period term of the flat phase angle from the flat phase angle, and then fits the flat phase angle deviation coefficient, accurately predicts the flat phase angle deviation of the current cycle and the next cycle, thereby accurately controlling the flat phase angle deviation to reach the upper and lower boundaries (thresholds), thereby achieving high-precision phase keeping for a low-orbit constellation system in a complex perturbation environment considering the non-spherical shape of the earth, atmospheric drag, gravity of the sun and the moon, and solar pressure.

[0036] 2. The high-precision phase holding method of a low-orbit constellation system of the present invention can achieve high-precision phase holding when the satellite position and velocity determination accuracy is not high and the positioning data sampling rate is low, thereby reducing the requirements of high-precision phase holding on the satellite orbit determination accuracy and orbit determination data sampling and storage, and has great engineering application value. BRIEF DESCRIPTION OF THE DRAWINGS

[0037] The specific implementation modes of the present invention are further described in detail below in conjunction with the accompanying drawings.

[0038] Figure 1 A flow chart of a high-precision phase keeping method for a low-orbit constellation system taking into account long-period terms proposed by the present invention is shown.

[0039] Figure 2 The time history of the short-period term of the flat phase angle deviation caused by various perturbation factors is shown.

[0040] Figure 3 The time history of the long-period term of the flat phase angle deviation caused by various perturbation factors is shown.

[0041] Figure 4 It shows the time history of the instantaneous track value and the flat track value of the phase angle deviation during the phase holding task cycle. DETAILED DESCRIPTION

[0042] In order to more clearly illustrate the present invention, the present invention is further described below in conjunction with preferred embodiments and accompanying drawings. Similar components in the accompanying drawings are represented by the same reference numerals. It should be understood by those skilled in the art that the content specifically described below is illustrative rather than restrictive, and should not be used to limit the scope of protection of the present invention.

[0043] The first embodiment of the present invention describes in detail a high-precision phase keeping method for a low-orbit constellation system taking into account long-period terms proposed by the present invention.

[0044] The present invention proposes to use the long-period term of the mean phase angle caused by the perturbations of the non-spherical shape of the earth, atmospheric drag, gravity of the sun and the moon, and solar pressure as a correction term for the mean orbit value, separate it from the mean phase angle, and then fit the mean phase angle deviation coefficient, accurately predict the mean phase angle deviation of the current cycle and the next cycle, and achieve precise control of the second threshold of the mean phase angle deviation by iteratively calculating the moment when the mean phase angle deviation reaches the second threshold; and achieve precise control of the first threshold of the mean phase angle deviation in the next phase holding cycle by iteratively calculating the size of the change in the semi-major axis. The method proposed in the present invention can achieve high-precision phase holding of a low-orbit constellation system in a complex perturbation environment including the non-spherical shape of the earth, atmospheric drag, gravity of the sun and the moon, and solar pressure.

[0045] By adopting the absolute phase holding method, each satellite in the constellation has the same phase holding strategy, but only has different reference phases. By using the control strategy of the present invention, the actual flat phase angle of each satellite tracks its reference flat phase angle, and the phase holding of the entire constellation can be achieved. Since the phase holding strategy of each satellite is the same, the following description will not be expanded for the constellation, but will be directly expanded for the satellite.

[0046] like Figure 1 As shown, the method mainly includes the following 7 steps:

[0047] S1. Determine the relationship between the satellite reference flat phase angle and time;

[0048] Any velocity and position state of a satellite corresponds to a set of orbital elements, which is called the instantaneous orbital elements. This set of orbital elements is constantly changing, and using it as the orbital elements of a satellite over a period of time will produce large errors. The analytical method using the averaging idea decomposes the changes in the instantaneous orbital elements into three categories with different properties: long-term terms, long-period terms, and short-period terms. Generally speaking, the average orbital elements (abbreviated as average elements, average orbit values, etc.) refer to a set of orbital elements corresponding to the average of the short-period terms. The average orbital elements can reflect the motion laws of the satellite at a higher level. Usually in satellite orbit design, orbital maneuvers and other tasks, the average orbital elements are used as design and control parameters.

[0049] The horizontal phase angle λ in the present invention is defined as the sum of the horizontal orbit value of the perigee argument ω and the horizontal orbit value of the horizontal anomaly M, that is, λ = ω + M. The change of the reference horizontal phase angle over time takes into account the gravitational field of the center of the earth, the non-spherical perturbation of the earth and the gravitational perturbation of the sun and the moon. The rate of change of the reference horizontal phase angle λ over time is recorded as n ref , the rate of change of the phase angle caused by the gravitational field at the center of the earth is recorded as n r , the rate of change of the phase angle caused by the non-spherical perturbation of the earth is recorded as n z , the rate of change of the phase angle caused by the sun's gravity is recorded as n s, the rate of change of the phase angle caused by the moon's gravity is recorded as n m ,but:

[0050] n ref =n r +n z +n s +n m (1)

[0051] In a possible implementation, assuming that at the start of the mission, the satellite's phase angle has been adjusted to within the holding interval, the start time of the entire mission is taken as the initial time of the first phase holding cycle. The reference phase angle of the satellite at the initial time t0 is denoted as λ ref (t0), then the reference flat phase angle of the satellite at time t is:

[0052] λ ref (t) = λ ref (t0)+n ref ·(t-t0) (2)

[0053] S2, sampling the position and velocity of the satellite at fixed intervals, converting them into actual average phase angles, and calculating the average phase angle long-period term at the corresponding moment;

[0054] Starting from the initial moment, the satellite position and velocity are sampled at fixed intervals within the same phase holding period, converted into actual mean phase angles, and the long-period term of the mean phase angle at the corresponding moment is calculated.

[0055] In a possible implementation, the position and velocity sampling interval is recorded as Δt, and the position vector r and velocity vector v values ​​of L+1 sampling points are obtained from the initial time, where L is the initial sampling interval number of a phase holding period. In the L·Δt time period, it is not necessary to judge whether the flat phase angle deviation exceeds the second threshold. The actual flat phase angle λ of the satellite at time t=t0+l·Δt (l∈[0,1,2,…,L]) is calculated by the method of converting position and velocity to orbit elements and the method of converting instantaneous orbit elements to average orbit elements. real (t).

[0056] In one possible implementation, the present invention substitutes the average orbital elements into the formula of the long-period term caused by the non-sphericity of the earth, atmospheric drag, the gravity of the sun and the moon, and the solar pressure perturbation, and obtains the sum λ1(t) of the long-period terms of the flat phase angle caused by each perturbation factor at the corresponding moment.

[0057] S3, subtracting the actual flat phase angle from the reference flat phase angle to calculate the flat phase angle deviation, and using the flat phase angle long period term for correction to obtain the flat phase angle deviation coefficient;

[0058] The actual flat phase angle of the satellite at the corresponding time is subtracted from its corresponding reference flat phase angle to obtain the flat phase angle deviation at time t = t0 + l·Δt:

[0059] Δλ(t)=λ real (t)-λ r ef(t) (3)

[0060] In a possible implementation, the present invention proposes to use the flat phase angle long period term for correction to obtain the flat phase angle deviation correction value Then fit and calculate the flat phase angle deviation coefficient. Calculate the flat phase angle long period term at the corresponding moment above to obtain the flat phase angle deviation correction value:

[0061]

[0062] λ l (t0) is the sum of the long-period terms of the mean phase angle caused by various perturbation factors at the initial time t0. After removing the influence of the long-period terms caused by the non-spherical shape of the earth, atmospheric drag, the gravitational force of the sun and the moon, and the solar light pressure perturbation, under the action of atmospheric drag, the semi-major axis changes approximately in a first-order linear manner over time, and the mean phase angle deviation correction value changes in a second-order nonlinear manner over time. The changing trend over time t is approximately:

[0063]

[0064] Among them: Δλ0 is the flat phase angle deviation at the initial time t0, Δλ1 is the first-order term of the flat phase angle drift caused by the semi-major axis deviation at the initial time t0, Δλ2 is the second-order term of the flat phase angle drift caused by the linear change of the semi-major axis over time. These three terms are the quadratic polynomial coefficients of the flat phase angle deviation changing with time.

[0065] According to the time t=t0+l·Δt(l∈[0,1,2,…,L]), the value of the satellite’s flat phase angle deviation Δλ(t), and expression (5), the least squares method is used to fit the flat phase angle deviation coefficients Δλ0, Δλ1, and Δλ2.

[0066] S4, predicting the flat phase angle deviation at the next sampling moment, and if it exceeds the second threshold, calculating the moment when the flat phase angle deviation exceeds the second threshold;

[0067] After determining the values ​​of the flat phase angle deviation coefficients Δλ0, Δλ1 and Δλ2, the flat phase angle deviation correction value at the forecast time t = t0 + (L + H) · Δt (H = 1, 2, 3, ...) is Simultaneously forecast the orbital elements of the satellite and the sun and moon at time t = t0 + (L + H) · Δt, obtain the forecast value of the long-period term of the mean phase angle, and thus obtain the forecast value Δλ of the mean phase angle deviation at time t = t0 + (L + H) · Δt p (t).

[0068]

[0069] In a possible implementation, if the flat phase angle deviation prediction value at time t=t0+(L+H)·Δt does not exceed the second threshold Δλ max , then in the time period [t0+(L+H-1)·Δt t0+(L+H)·Δt], no phase-holding orbit control is applied, and the actual flat phase angle and the flat phase angle long-period term at the moment t0+(L+H)·Δt are obtained by the method of step S2, and the quadratic polynomial coefficient (deviation coefficient) of the flat phase angle deviation varying with time is updated by the method of step S3 after considering the flat phase angle deviation correction value at the moment t0+(L+H)·Δt.

[0070] If the flat phase angle deviation prediction value Δλ at time t = t0 + (L + H) · Δt p (t) exceeds the second threshold Δλ max , then let:

[0071] Δλ p (t) = Δλ max (7)

[0072] Combined with equation (6), the phase angle deviation reaches the second threshold Δλ max The time t f

[0073] t f =t0+(L+H-1)·Δt+Δt n (8)

[0074] Where Δt n <Δt.

[0075] In a possible implementation, the present invention introduces a free drift time Δt that is less than one sampling period. n , which makes the second threshold control of the flat phase angle deviation more accurate. l (t)-λ l (t0)] involves the prediction of satellite orbit elements, and an iterative method is used to determine the time when the phase angle deviation exceeds the second threshold t f .

[0076] S5, determining an initial solution for the change in the semi-major axis of the phase-keeping orbit maneuver according to the moment when the second threshold value is exceeded and the flat phase angle deviation coefficient;

[0077] Determine the satellite from time t0 to time t f The coefficients Δλ0, Δλ1 and Δλ2 of the flat phase angle deviation correction value change with time. The deviation Δa(t0) between the actual semi-major axis and the reference semi-major axis at the initial time t0 can be calculated by the following formula.

[0078]

[0079] In the formula, μ is the gravitational constant of the earth, a c The first-order linear rate of change of the satellite's orbital semi-major axis over time during this period can be calculated by the following formula:

[0080]

[0081] So we can determine the time t f The deviation Δa(t f ).

[0082]

[0083] According to the phase angle deviation interval [Δλ min Δλ max ], and the second-order term of the flat phase angle drift Δλ2, the duration of the next phase between two orbital height increases can be estimated as:

[0084]

[0085] The semi-major axis drift between two orbital altitude increases is:

[0086]

[0087] Therefore, the initial solution for the change in the semi-major axis required for the phase-keeping orbit maneuver is:

[0088]

[0089] S6. Predict the flat phase angle deviation of the next phase holding period according to the flat phase angle deviation coefficient of the current phase holding period, and use the long period term for correction to determine the correction solution of the change of the semi-major axis of the phase holding orbit maneuver;

[0090] The phase angle deviation of the next phase holding period is predicted based on the phase angle deviation coefficient of the current phase holding period, and the long period term is used for correction. A correction solution for the change in the semi-major axis of the phase holding orbit maneuver is determined so that the minimum value of the phase angle deviation of the next phase holding period just reaches the first threshold.

[0091] In a possible implementation, the present invention takes into account the differences in the sun, moon and satellite orbital elements in different phase-holding periods, and the long-period term of the flat phase angle caused by the non-spherical shape of the earth, atmospheric drag, the gravity of the sun and moon, and the solar light pressure perturbation is also different. Therefore, the semi-major axis change amount of the phase-holding orbit maneuver obtained in step S5 is only an initial solution. To accurately control the flat phase angle deviation of the next phase-holding period to reach the first threshold, the semi-major axis change amount obtained above should also be corrected.

[0092] In a possible implementation, assuming that the current phase holding period is the kth phase holding period, the flat phase angle deviation coefficients Δλ0, Δλ1, and Δλ2 defined in the formula are respectively denoted as Δλ 0，k , Δλ 1，k , Δλ 2，k , then the semi-major axis deviation Δa at the end of the kth phase holding period is calculated according to the formula k (t f ), calculate the initial solution δa′ of the semi-major axis change of the kth phase holding period according to the formula k , let the semi-major axis correction be Δa m , then the modified solution of the change in the semi-major axis of the kth phase holding period is δa k ,satisfy:

[0093] δa k =δa′ k +Δa m (15)

[0094] Then the semi-major axis deviation value Δa at the initial moment of the k+1th phase holding period is k+1 (t0),

[0095] Δa k+1 (t0) = Δa k (t f )+δa k (16)

[0096] Predict the phase angle deviation coefficient Δλ of the k+1th phase holding period 0，k+1 , Δλ 1，k+1 , Δλ 2，k+1 for:

[0097]

[0098] Predict the sun, moon and satellite orbit elements of the k+1th phase holding period and calculate the corresponding long-period term of phase angle deviation [λ l (t)-λ l (t0)], and obtain the phase angle deviation Δλ(k+1, t) of the k+1th phase holding period

[0099] Δλ(k+1,t)=Δλ 0，k+1 +Δλ 1，k+1 (t-t0)+Δλ 2，k+1 (t-t0) 2 +[λ l (t)-λ l (t0)] (18)

[0100] Determine the corrected solution δa of the semi-major axis change of the phase-keeping orbit maneuver by iterative method k , so that the minimum phase angle deviation min[Δλ(k+1, t)] of the next phase holding period just reaches the first threshold Δλ min .

[0101] S7, applying a velocity increment when the flat phase angle deviation value reaches the second threshold, taking the orbital maneuver end time as the initial time of the next control cycle, and repeating steps S2 to S7 until the phase holding of the entire phase holding mission cycle of the satellite is completed.

[0102] In one possible implementation, for a near-circular orbit, the corresponding semi-major axis change δa is k The required lateral velocity increment is:

[0103]

[0104] The phase-keeping orbit maneuvering thruster working time is:

[0105]

[0106] Where m is the satellite mass, F p The phase angle deviation obtained in step S4 reaches the second threshold time t f =t0+(L+H-1)·Δt+Δt n Apply tangential thrust to obtain the lateral velocity increment Δv T The total duration of the phase holding period is defined as (L+H-1)·Δt+Δt n +Δt m The orbital maneuver end time t = t0 + (L + H-1) · Δt + Δt n +Δt m It serves as the initial moment of the next control cycle, and repeats S2 to S7 until the phase holding of the entire phase holding mission cycle of the satellite is completed.

[0107] Another embodiment of the present invention is combined with the attached Figure 2-4 The present invention is further described.

[0108] like Figure 2As shown in the figure, the time history of the short-period term of the average phase angle deviation caused by each perturbation factor is given. The amplitude of the instantaneous orbit value of the average phase angle deviation caused by the main term of the earth belt harmonic is about 5.5×10 -2 (°), the amplitude of the instantaneous orbit value of the phase angle deviation caused by the main term of the earth field harmonics is about 5.3×10 -3 (°), the amplitude of the instantaneous orbit value of the phase angle deviation caused by atmospheric drag is about 9.0×10 -8 (°), the amplitude of the instantaneous orbit value of the phase angle deviation caused by the gravity of the sun is about 5.8×10 -6 (°), the amplitude of the instantaneous orbit value of the phase angle deviation caused by the lunar gravity is about 1.5×10 -5 (°), the amplitude of the instantaneous orbit value of the phase angle deviation caused by solar pressure is about 3.6×10 -6 (°).

[0109] Assuming that the task requires the instantaneous orbit value of the phase angle deviation to be kept within the range of ±0.1°, the short-period terms caused by the above perturbations and the high-order short-period terms not considered in the instantaneous-to-flat conversion are combined, and a certain control margin is retained to keep the flat-track value of the flat phase angle deviation within the range of [-0.03° 0.03°], that is, Δλ min =-0.03°, Δλ max =0.03°.

[0110] In this embodiment, the non-spherical perturbation of the earth, the atmospheric drag perturbation, the gravitational perturbation of the sun and the moon, and the solar pressure perturbation are considered. The non-spherical perturbation of the earth adopts the WGS84-EGM96 model, the Degree and Order are set to 21, and the atmospheric drag perturbation adopts the US 1976 Standard Atmospheric Density Model. The satellite reference orbit altitude is taken as 800km, and the initial orbit elements are shown in Table 1. The satellite's surface mass ratio is taken as 0.0175, the drag coefficient is taken as 2.2, the phase holding mission period is taken as 90 days, the orbit determination data sampling interval is taken as 2 hours, the three-axis position determination error is taken as 5m (1σ), and the three-axis velocity determination error is taken as 0.02m / s (1σ). The scene epoch time is 1 Jan 2022 00:00:00:000 UTCG.

[0111] Table 1 Initial average orbital elements

[0112] Number of orbital elements value a 7178.137km e 0.001 i 87° Ω 20° ω 90° f 10° λ 99.98°

[0113] Determine the relationship between the rate of change of the satellite's reference flat phase angle over time and the change of the reference flat phase angle over time.

[0114] When the satellite nominal semi-major axis, eccentricity, orbital inclination and epoch time are given, the rate of change of the reference phase angle with time can be directly obtained by analyzing the low-order terms of the Earth's center gravitational field and the Earth's non-spherical perturbation, as well as the long-term terms of the Sun-Moon three-body gravitational perturbation:

[0115] n ref =5.941×10 -2 (° / s)

[0116] Assume that at the start of the mission, the satellite's phase angle deviation is exactly at the upper boundary Δλ max The start time of the entire mission is taken as the initial time of the first phase holding cycle. Then the reference flat phase angle λ of the satellite at the initial time t0 is ref (t0) should be 99.95°. The reference phase angle of the satellite at time t is:

[0117] λ ref (t) = 99.95° + 5.941 × 10 -2 t(°)

[0118] Starting from the initial moment, the satellite position and velocity are sampled at fixed intervals within the same phase holding period, converted into actual mean phase angles, and the long-period term of the mean phase angle at the corresponding moment is calculated.

[0119] Take L = 30, and start from the initial time, sample the value of the satellite's position vector r and velocity vector v every Δt = 2h, and calculate the actual flat phase angle λ of the satellite at time t = l·7200sec (l∈[0,1,2,…,30]) by the method of converting position and velocity to orbit elements and the method of converting instantaneous orbit elements to average orbit elements. real Substituting the average orbital elements into the long-period term formula caused by the non-spherical shape of the earth, atmospheric drag, the gravitational force of the sun and the moon, and the solar light pressure perturbation, the sum of the long-period terms of the average phase angle at the corresponding moment λ is obtained l (1·7200sec).

[0120] like Figure 3 As shown in the figure, the time history of the long-period term caused by various perturbation factors is given. The amplitude of the instantaneous orbit value of the flat phase angle deviation caused by the main term of the earth belt harmonic is about 3.0×10 -5 (°), the amplitude of the instantaneous orbit value of the phase angle deviation caused by atmospheric drag is about 2.0×10 -16 (°), the amplitude of the instantaneous orbit value of the phase angle deviation caused by the gravity of the sun is about 6.0×10 -3 (°), the amplitude of the instantaneous orbit value of the phase angle deviation caused by the lunar gravity is about 3.0×10 -3(°), the amplitude of the instantaneous orbit value of the phase angle deviation caused by solar pressure is about 5.0×10 -6 (°).

[0121] It can be seen that the amplitude of the instantaneous orbit value of the flat phase angle deviation caused by the combined effect of various perturbations is about 1.0×10 -2 (°), which is mainly composed of the long-period terms caused by the gravity of the sun and the moon. Therefore, the long-period term caused by the gravity of the sun and the moon plays a dominant role in the long-period terms caused by various perturbation factors. Also, since the long-period term caused by the gravity of the sun and the moon is already of the same order of magnitude as the accuracy of the mean phase angle deviation, it is very necessary to use the long-period term to correct the mean phase angle deviation in the perturbation model considering the influence of the sun and the moon.

[0122] The flat phase angle deviation is calculated and corrected with the flat phase angle long period term to obtain the flat phase angle deviation coefficient.

[0123] The actual satellite phase angle at the corresponding time is subtracted from its corresponding reference phase angle to obtain the phase angle deviation Δλ(l·7200sec) at time t=l·7200sec, and the long-period term [λ l (l·7200sec)-λ l (0·7200sec)] to make corrections and obtain the correction value of the flat phase angle deviation Table 2 gives the flat phase angle deviation, long period correction term and flat phase angle deviation correction term at different simulation times, where Δλ(l·Δt), λ l (t)-λ l (t0), The unit is degree (°).

[0124] According to the variation law of the flat phase angle deviation correction value with time t, the flat phase angle deviation coefficient can be calculated as:

[0125]

[0126] The flat phase angle deviation at the next sampling moment is predicted, and the moment when the flat phase angle deviation exceeds the second threshold is determined.

[0127] According to the factor parameters of the flat phase angle deviation changing with time, it can be predicted that the flat phase angle at t = 31.7200s is -0.0011°, which does not reach the boundary value of 0.03°. No phase-holding orbit control is applied in the time period [(31)·7200 (32)·7200]. The actual flat phase angle deviation and the flat phase angle long-period term at time (31)·7200s are continuously obtained, and the flat phase angle deviation coefficient is updated.

[0128] The above steps are used to sequentially predict the flat phase angle deviation at time (32)·7200s and (33)·7200s until the flat phase angle deviation at the corresponding time exceeds the drift boundary (second threshold) Δλ max , determine the satellite uncontrolled free drift time starting from time t0 = 0 as t end =(30+159)·7200s+4479s, that is, H=160, Δt n =4479s, and the factor parameters of the flat phase angle deviation during this period of time.

[0129]

[0130] Table 2 Flat phase angle deviation, long period correction term and flat phase angle deviation correction term at different simulation times

[0131]

[0132] Determine the initial solution for the change in the semi-major axis of the phase-keeping orbit maneuver.

[0133] Calculate the deviation Δa(t0) between the actual semi-major axis and the reference semi-major axis at the initial time t0.

[0134] Δa(0)≈14.886m

[0135] Calculate the first-order linear rate of change of the satellite's orbital semi-major axis over time during this period

[0136]

[0137] So we can determine the time t end The deviation Δa(t f ):

[0138] Δa(t f )≈-14.333m

[0139] According to the flat phase angle deviation interval [-0.03° 0.03°] and the second-order term of the flat phase angle drift Δλ2, the duration between two orbital altitude increases in the next phase can be estimated as:

[0140] ΔT=15.548day

[0141] The semi-major axis drift between two orbital altitude increases is:

[0142] Δa d =-28.748m

[0143] Therefore, the change in the semi-major axis required for the phase-keeping orbit maneuver is:

[0144] δa′=28.707m

[0145] The phase angle deviation of the next phase holding period is predicted based on the phase angle deviation coefficient of the current phase holding period, and the long period term is used for correction. A correction solution for the change in the semi-major axis of the phase holding orbit maneuver is determined so that the minimum value of the phase angle deviation of the next phase holding period just reaches the first threshold.

[0146] Assume that the current phase holding cycle is the kth phase holding cycle, and the phase angle deviation coefficient is Δλ 0，k , Δλ 1，k , Δλ 2，k , determine Δλ 2，k =2.321×10 -15 / s 2 , determine the semi-major axis deviation value Δa at the end of the kth phase holding period k (t f )=-14.333m, and the initial solution of the semi-major axis change of the kth phase holding period δa′ k =28.707m. The semi-major axis correction solution is calculated by iterative method. First, assume that the semi-major axis correction value Δa m =0, then the correction solution of the semi-major axis change is:

[0147] δa k =28.707m

[0148] Then the semi-major axis deviation value at the initial moment of the k+1th phase holding period is:

[0149] Δa k+1 (t0) = 14.374m

[0150] Predict the phase angle deviation coefficient Δλ of the k+1th phase holding period 0，k+1 , Δλ 1，k+1 , Δλ 2，k+1 for:

[0151]

[0152] Predict the sun, moon and satellite orbit elements of the k+1th phase holding period and calculate the corresponding long-period term of phase angle deviation [λ l (t)-λ l (t0)], and obtain the phase angle deviation Δλ(k+1, t) of the k+1th phase holding period, where the minimum value is:

[0153] min[Δλ(k+1,t)]=-0.291°>-0.03°

[0154] It shows that the applied semi-major axis change is slightly small, and the semi-major axis change should be increased. After several iterations, the semi-major axis correction Δa is determined.m =0.106m, then the correction solution of the semi-major axis change is:

[0155] δa k =28.813m

[0156] The semi-major axis deviation value at the initial moment of the k+1th phase holding period:

[0157] Δa k+1 (t0) = 14.480m

[0158] The updated phase angle deviation coefficient Δλ of the k+1th phase holding period 0，k+1 , Δλ 1，k+1 , Δλ 2，k+1 for:

[0159]

[0160] It is predicted that the minimum value of the phase angle deviation in the k+1th phase holding period is exactly -0.03°.

[0161] When the calculated flat phase angle deviation reaches the upper boundary, a velocity increment is applied. The orbit maneuver end time is used as the initial time of the next control cycle, and the above steps are repeated until the phase holding of the entire mission cycle of the satellite is completed.

[0162] For a near-circular orbit, the corresponding change in the semi-major axis is δa k The required lateral velocity increment is:

[0163] Δv T =1.496×10 -2 m / s

[0164] Assuming the satellite mass is 200kg and the thrust value is 80mN, the working time of the phase-keeping orbit maneuvering thruster is:

[0165] Δt m =37.390s

[0166] At the moment t when the obtained flat phase angle deviation reaches the upper boundary f = (30 + 159) · 7200s + 4479s Apply tangential thrust to obtain the lateral velocity increment Δv T The total duration of the phase hold period is defined as t end =(30+159)·7200s+4479s+37.390s. The orbital maneuver end time t = t end As the initial moment of the next control cycle, the above steps are repeated until the phase holding of the entire task cycle is completed.

[0167] Figure 4The simulation results of 6 phase holding cycles are given. It can be seen that the instantaneous orbit value of the phase angle deviation is kept within ±0.1°, and the flat orbit value of the phase angle deviation is basically controlled within ±0.03°. The method proposed in the present invention realizes high-precision phase holding of the low-orbit constellation system under complex perturbation environment.

[0168] Obviously, the above embodiments of the present invention are merely examples for clearly illustrating the present invention, and are not limitations on the implementation methods of the present invention. For ordinary technicians in the relevant field, other different forms of changes or modifications can be made based on the above description. It is impossible to list all the implementation methods here. All obvious changes or modifications derived from the technical solution of the present invention are still within the protection scope of the present invention.

Claims

1. A high-precision phase keeping method for a low-orbit constellation system considering long-period terms, characterized in that: The method includes: S1. Determine the relationship between the satellite reference flat phase angle and time; S2, sampling the satellite position and velocity once every fixed time, converting it into the actual flat phase angle, and calculating the flat phase angle long-period term at the corresponding moment; The step S2 comprises: Substituting the average orbital elements into the long-period term formula caused by the non-sphericity of the earth, atmospheric drag, the gravitational force of the sun and the moon, and the solar light pressure perturbation, the flat phase angle long-period term is obtained; S3, subtracting the actual flat phase angle from the reference flat phase angle at the corresponding moment, calculating the flat phase angle deviation, and using the flat phase angle long period term for correction to obtain the flat phase angle deviation coefficient; The step S3 comprises: The flat phase angle long-period term is used for correction to obtain a flat phase angle deviation correction value, and a flat phase angle deviation coefficient is calculated by fitting according to a trend of the flat phase angle deviation correction value over time; S4, predicting the flat phase angle deviation at the next sampling moment, and if it exceeds the second threshold, calculating the moment when the flat phase angle deviation exceeds the second threshold; S5, determining an initial solution for the change in the semi-major axis of the phase-keeping orbit maneuver according to the moment when the second threshold value is exceeded and the flat phase angle deviation coefficient; S6. Predict the flat phase angle deviation of the next phase holding period according to the flat phase angle deviation coefficient of the current phase holding period, and use the long period term for correction to determine the correction solution of the change of the semi-major axis of the phase holding orbit maneuver; The step S6 comprises: Predicting the sun, moon and satellite orbital elements of the next phase-holding period, calculating the corresponding long-period term of the phase angle deviation, obtaining the mean phase angle deviation of the next phase-holding period, and determining a correction solution for the change in the semi-major axis of the phase-holding orbit maneuver that makes the minimum mean phase angle deviation of the next phase-holding period just reach a first threshold; S7, applying a velocity increment when the flat phase angle deviation value reaches the second threshold, taking the orbital maneuver end time as the initial time of the next control cycle, and repeating steps S2 to S7 until the phase holding of the entire phase holding mission cycle of the satellite is completed.

2. The method according to claim 1, characterized in that The second threshold is the value of the right endpoint of the flat phase angle deviation interval, and the value of the left endpoint of the flat phase angle deviation interval is the first threshold.

3. The method according to claim 1, characterized in that: The step S1 comprises: Determine the rate of change of the reference flat phase angle over time based on the gravitational field at the center of the Earth, the non-spherical perturbation of the Earth, and the gravitational perturbation of the Sun and the Moon; The reference flat phase angle of the satellite at a certain moment is determined according to the reference flat phase angle of the satellite at the initial moment of the satellite phase keeping mission cycle and the rate of change of the reference flat phase angle over time.

4. The method according to claim 1, characterized in that The flat phase angle deviation coefficient is a quadratic polynomial coefficient of the flat phase angle deviation correction value changing with time, including the flat phase angle deviation at the initial moment, the first-order term of the flat phase angle drift caused by the semi-major axis deviation at the initial moment, and the second-order term of the flat phase angle drift caused by the linear change of the semi-major axis with time.

5. The method according to claim 1, characterized in that The step S4 further comprises: If the flat phase angle deviation forecast value does not exceed the second threshold, steps S2 to S3 are repeated to update the flat phase angle deviation coefficient until the flat phase angle deviation at this moment exceeds the second threshold; If the flat phase angle deviation forecast value exceeds the second threshold, a free drift duration less than one sampling period is introduced to calculate the time when the flat phase angle deviation exceeds the second threshold.

6. The method according to claim 1, characterized in that The step S5 comprises: According to the flat phase angle deviation coefficient, calculating the deviation between the actual semi-major axis and the reference semi-major axis at the initial moment and the first-order linear change rate of the orbital semi-major axis of the satellite from the initial moment to the moment when the flat phase angle deviation exceeds the second threshold, thereby determining the deviation between the actual semi-major axis and the reference semi-major axis at the moment when the flat phase angle deviation exceeds the second threshold; According to the flat phase angle deviation interval and the second-order term of the flat phase angle drift, the duration between two orbital altitude increases of the next phase is estimated; An initial solution for the change in the semi-major axis of the orbital maneuver is determined based on the semi-major axis drift between two orbital altitude increases and the deviation between the actual semi-major axis and the reference semi-major axis when the flat phase angle deviation exceeds a second threshold.

7. The method according to any one of claims 1 to 6, characterized in that: The orbital elements of the sun, moon and satellite are different in different phase holding periods, and the long-period terms of the mean phase angle caused by the non-spherical shape of the earth, atmospheric drag, the gravitational force of the sun and moon, and the solar light pressure perturbation are also different.

Images