Coplanar fly-around formation configuration maintaining method and system under influence of space perturbation

Abstract

The invention discloses a coplanar flying-around formation configuration maintenance method and system under the influence of space perturbation, relates to the technical field of space satellite formation and configuration maintenance, and aims to solve the problem that a conventional satellite formation configuration maintenance method is difficult to maintain the configuration of coplanar flying-around formation satellites under the influence of J2 item perturbation and atmospheric drag perturbation. According to the method, the perturbation force is directly compensated on the control channel, and the nonsingular terminal sliding mode can process the nonlinear state-space equation, so that linearization processing does not need to be performed on the kinetic equation added with the perturbation force again, and inherent errors caused by linearization are avoided. Therefore, according to the technical scheme, the configuration of the coplanar flying-around formation satellite under the influence of J2 item perturbation and atmospheric drag perturbation can be effectively maintained.

Description

Technical Field

[0001] The present invention relates to the technical field of space satellite formation and configuration maintenance, and in particular to a method and system for maintaining a coplanar flyby formation configuration under the influence of space perturbations. Background Art

[0002] Nowadays, space mission requirements are becoming increasingly complex and diversified. Traditional single-star systems are difficult to achieve more powerful functions within a limited payload space. Constellation systems also have disadvantages such as poor flexibility, high cost, and complex systems. In order to address this problem, it has become more important to replace traditional solutions with flight formations composed of multiple satellites through precise measurement and collaborative operation. To achieve mission requirements and maintain a long-term stable formation structure, many challenges need to be overcome. After the launch, the satellite formation needs to enter and adjust the orbit one by one to form a formation structure. Satellite formation flight needs to consider the relative motion and Doppler effect of satellites in space. In addition, the formation of satellites orbiting the earth must also consider the orbital constraints of satellites. Therefore, satellite formations are more complex in terms of control, navigation, and data processing. Satellites operating in orbit for a long time in space will also be affected by a variety of perturbations: such as non-spherical perturbations, atmospheric drag perturbations, solar light pressure perturbations, and solar and lunar gravitational perturbations. Even if it is a well-designed natural formation configuration, its member satellites will drift in the satellite orbit due to the existence of perturbations. In order to successfully complete satellite Earth observation, communication, navigation and other flight missions and avoid collision risks between formation member satellites, it is necessary to maintain the orbits of formation member satellites so that the motion parameters corresponding to each point on the orbit change as little as possible compared with the parameters of the previous cycle. Traditional satellite formation configuration maintenance methods often require the analytical solution of the relative motion equation between satellites, and the orbital elements and motion equations of the orbiting satellites are calculated based on the analytical solution. However, the analytical solution process of the relative motion equation with perturbation force is very complicated, and the perturbation term needs to be linearized, which will inevitably lead to deviations in the final result. Therefore, traditional satellite formation configuration maintenance methods are difficult to maintain the configuration of coplanar orbiting formation satellites under the influence of J2 perturbation and atmospheric drag perturbation. Summary of the invention

[0003] The purpose of the present invention is to propose a method and system for maintaining the configuration of a coplanar flyby formation under the influence of space perturbations in order to address the problem that traditional satellite formation configuration maintenance methods are difficult to maintain the configuration of coplanar flyby formation satellites under the influence of J2 perturbations and atmospheric drag perturbations.

[0004] The technical solution adopted by the present invention to solve the above technical problems is:

[0005] A method for maintaining a coplanar fly-around formation configuration under the influence of space perturbations comprises the following steps:

[0006] Step 1: For the J2 perturbation and the atmospheric drag perturbation, establish the J2 perturbation force mathematical model and the atmospheric drag perturbation force mathematical model respectively;

[0007] Step 2: Substitute the mathematical model of the J2 perturbation force and the mathematical model of the atmospheric drag perturbation force into the CW equation to obtain the state space equation with the perturbation force added;

[0008] Step 3: Use the orbit six-element method to design a natural formation configuration for coplanar flyby;

[0009] Step 4: Based on the natural formation configuration of coplanar flyby, and using the fourth-order Runge-Kutta method, the position and velocity of the orbiting satellite relative to the primary satellite in the LVLH coordinate system at each moment under the condition of two-body motion are obtained;

[0010] Step 5: Take the position and velocity of the orbiting satellite relative to the host satellite in the LVLH coordinate system at each moment as the current coordinates, and use the relative motion dynamics equation with perturbation force to design a non-singular terminal sliding mode controller;

[0011] Step 6: Use a non-singular terminal sliding mode controller to maintain the coplanar flyby formation configuration.

[0012] Furthermore, the J2 perturbation force mathematical model is expressed as:

[0013]

[0014] Normalized spherical harmonics The specific solution is obtained by the following recursive formula:

[0015]

[0016] Among them, the first term on the right side of the mathematical model of J2 perturbation dynamics is represents the gravity at the center of the earth, and the last two terms is the non-spherical gravity of the Earth, R e and μ are the reference radius of the Earth and the Earth's gravitational constant, respectively. is the geocentric distance, geocentric latitude and geocentric longitude of the satellite in the international earth reference system, n is the order of the selected gravity field model, m is the degree of the selected gravity field model, when m≠n, C nm and S nm The earth is described as a chessboard pattern with similar concavities and convexities, which is called the field harmonic term of the earth's non-spherical perturbation; when m = n, C nm and S nm The earth will be described as a convex and concave fan-shaped pattern, which is called the fan harmonic of the earth's non-spherical perturbation; nm and S nm Normalize it, then in δ is an intermediate variable, Indicates the state of the previous iteration step n-1 and the previous iteration step m-1; represents the state of the previous iteration step n-1 and the current condition m, Represents the state of the previous iteration step m-1 and the current condition n, Indicates the state under the current conditions n and m;

[0017] The mathematical model of atmospheric drag force is expressed as:

[0018]

[0019]

[0020]

[0021] Among them, C D is the drag coefficient, is the surface-to-mass ratio, ρ is the atmospheric density, is the speed of the satellite relative to the atmosphere, v and v a are the speeds of the satellite and atmosphere relative to the center of the earth, respectively; X and Y are the speeds of the satellite and atmosphere relative to the center of the earth, respectively; φ is the geocentric latitude of the satellite's location, ω e is the angular velocity vector of the Earth.

[0022] Furthermore, the specific steps of step 2 are:

[0023] The mathematical models of the J2 perturbation force and the atmospheric drag perturbation force are integrated into the general form of the CW equation as the state space equation of the study. The general form of the CW equation is expressed as:

[0024]

[0025] Where ρ = r c -r s =[xyz] T is the position coordinate vector of the companion spacecraft relative to the main spacecraft, is the acceleration coordinate vector of the companion spacecraft relative to the main spacecraft, Δf x , Δf y , Δf z are the three-axis components of the acceleration caused by all the perturbations and control forces on the companion spacecraft relative to the main spacecraft in the relative coordinate system except the gravity of the center of the earth, and ε is the average orbital velocity;

[0026] Add the perturbation force to the general CW equation, that is, write the mathematical model of the J2 perturbation force and the mathematical model of the atmospheric drag perturbation force into the form of acceleration and bring it into Δf x , Δfy , Δf z Among the three items, considering the mathematical model of the J2 perturbation force and the mathematical model of the atmospheric drag perturbation force, the acceleration generated by the non-spherical perturbation J2 item is as follows:

[0027] J2(r)=-(3 / 2)(J2μRe 2 / r 4 )[(1-3sin 2 isin 2 θ)x

[0028] +(2sin 2 isinθcosθ)y+(2sinicosisinθ)z]

[0029] Where r is the position vector of the satellite, J2 is the second spherical harmonic term of the Earth's gravitational field potential, x, y, z represent the spacecraft coordinates in the LVLH coordinate system, and R e is the mean equatorial radius of the Earth, i is the orbital inclination, θ is the ascending node angular distance, and x, y, z are the unit direction vectors of the three axes of the corresponding coordinate system;

[0030] Substituting J2(r) into the general form of the CW equation and linearizing it, we get the improved relative motion dynamics equation, which is expressed as:

[0031]

[0032] Among them, r ref is the constant radius of the circular reference orbit, i ref is the initial orbital inclination of the reference spacecraft, t is a given time, and the parameters c, k, l, q are used to correct the circular reference orbit period caused by the nodal drift caused by the J2 effect:

[0033]

[0034] Φ0=cos -1 [cosi c cosi s +sini c sini s cosΔΩ0]

[0035]

[0036] Among them, s, γ0, Φ0, ΔΩ0 are intermediate variables, is the first-order derivative of Δz0, Δz0 is the difference in the z-axis coordinates between the main spacecraft and the companion spacecraft in the LVLH system in the initial state, i sis the initial orbital inclination of the main spacecraft. Since the main spacecraft is initially located in a circular reference orbit, i s Set to i ref ;i c is the initial orbital inclination of the reference spacecraft;

[0037] Based on the improved relative motion dynamics equation, the acceleration caused by atmospheric drag is substituted and expressed as:

[0038]

[0039] Among them, f Drag,x 、f Drag,y 、f Drag,z They represent the additional acceleration due to aerodynamic drag in the radial direction, in the orbit, and in the normal direction of the orbit respectively; the simultaneous equations Using a circular reference orbit, the velocity of the spacecraft relative to the rotating atmosphere is given by:

[0040]

[0041] in, are the first-order derivatives of x, y, and z respectively, Substituting into the equation Then we get the following state space equation with the perturbation force added:

[0042]

[0043] in, express Calculate in the x-axis direction, express Calculate in the y-axis direction, express The calculation is performed in the z-axis direction.

[0044] Furthermore, the formation configuration parameters of the coplanar flyby natural formation configuration include: the short radius p of the flyby ellipse, the motion amplitude S in the direction perpendicular to the orbital plane, the initial phase difference α, the initial phase σ of the flyby ellipse, and the distance l along the track from the center of the flyby ellipse;

[0045] p=ae A , s = aΔi, α = σ-ψ, l = aΔλ, obtained by the following expression:

[0046]

[0047] The mapping expression of the orbiting star eccentricity is as follows:

[0048]

[0049] Among them, arctan is the inverse sine function, ψ and Δλ are intermediate variables;

[0050] For a circular reference orbit, e ref =0, then

[0051] The mapping expression of the orbital inclination of the orbit around the star is as follows:

[0052]

[0053] The mapping expression of the right ascension of the ascending node of the orbiting star is Ω cir =Ω ref +ΔΩ;

[0054] The mapping expression of the perigee angle of the orbiting star is:

[0055] The mapping expression of the orbiting star mean anomaly is:

[0056]

[0057] Among them, a is the semi-major axis of the orbit, indicating the size of the orbit, and e A represents the eccentricity of the two satellites relative to the ascending node, e cir and e ref are the orbital eccentricities of the orbiting star and the primary star, ω cir and ω ref are the perigee arguments of the orbiting star and the primary star, M cir and M ref are the mean anomalies of the orbiting star and the primary star, i cir is the orbital inclination of the star, Ω cir and Ω ref are the right ascension of the ascending node of the orbiting star and the primary star respectively, ψ and Δλ are intermediate parameters, ΔΩ is the difference in the right ascension of the ascending node of the two satellite orbits, Δi is the angle between the two satellite orbits, and arccos is the inverse cosine function.

[0058] Furthermore, the specific steps of step five are:

[0059] The state space equation can be written as The form is expressed as:

[0060]

[0061] where u x ,u y ,u z is the control input of the three axes, t is the simulation time, and the above equation is written as a second-order nonlinear dynamic system in the following form:

[0062]

[0063] We can get:

[0064]

[0065] The relative position error and relative velocity error of the satellite are tracked, so x = [x1 x2] T From e = [e1e2] T Substitute, and e1 is the position error, e2 is the speed error, then the sliding surface is expressed as The non-singular terminal sliding mode controller u is designed as follows:

[0066]

[0067] Where f(x) is the non-perturbation term in the state space equation, g(x) represents the perturbation force, b(x) is the gain matrix, β, e = [e1 e2] T , l g and η are both adjustment parameters in sliding mode control, and β is 1e -3 , l g ≥g(x), and l g Take 2e -10 , η is 1e -11 .

[0068] A coplanar fly-around formation configuration maintenance system under the influence of space perturbations, the system comprising a mathematical model building module, a state space equation building module, a formation configuration building module, a two-body motion module and a sliding mode control module;

[0069] The mathematical model building module is used to establish a J2 perturbation force mathematical model and an atmospheric drag perturbation force mathematical model respectively for the J2 perturbation and the atmospheric drag perturbation;

[0070] The state space equation construction module is used to bring the mathematical model of the J2 term perturbation force and the mathematical model of the atmospheric drag perturbation force into the CW equation to obtain the state space equation with the perturbation force added;

[0071] The formation configuration building module is used to design a natural formation configuration for coplanar flyby using the orbit six-element method;

[0072] The two-body motion module is used to obtain the position and velocity of the orbiting satellite relative to the primary satellite in the LVLH coordinate system at each moment under the two-body motion condition based on the natural formation configuration of coplanar flyby and the fourth-order Runge-Kutta method;

[0073] The sliding mode control module is used to take the position and velocity of the orbiting satellite relative to the main satellite in the LVLH coordinate system at each moment as the current coordinates, and use the relative motion dynamics equation with perturbation force to design a non-singular terminal sliding mode controller, and finally use the non-singular terminal sliding mode controller to complete the maintenance of the coplanar flyby formation configuration.

[0074] Furthermore, the J2 perturbation force mathematical model is expressed as:

[0075]

[0076] Normalized spherical harmonics The specific solution is obtained by the following recursive formula:

[0077]

[0078] The first term on the right side of the equation is represents the gravity at the center of the earth, and the last two terms and is the non-spherical gravity of the Earth, R e and μ are the reference radius of the Earth and the Earth's gravitational constant, respectively. is the geocentric distance, geocentric latitude and geocentric longitude of the satellite in the international earth reference system, n is the order of the selected gravity field model, m is the degree of the selected gravity field model, when m≠n, C nm and S nm The earth is described as a chessboard pattern with similar concavities and convexities, which is called the field harmonic term of the earth's non-spherical perturbation; when m = n, C nm and S nm The earth will be described as a convex and concave fan-shaped pattern, which is called the fan harmonic of the earth's non-spherical perturbation; nm and S nm Normalize it, then in δ is an intermediate variable, Indicates the state of the previous iteration step n-1 and the previous iteration step m-1; represents the state of the previous iteration step n-1 and the current condition m, Represents the state of the previous iteration step m-1 and the current condition n, Indicates the state under the current conditions n and m;

[0079] The mathematical model of atmospheric drag force is expressed as:

[0080]

[0081]

[0082]

[0083] Among them, C D is the drag coefficient, is the surface-to-mass ratio, ρ is the atmospheric density, is the speed of the satellite relative to the atmosphere, v and v a are the speeds of the satellite and atmosphere relative to the center of the earth, respectively; X and Y are the speeds of the satellite and atmosphere relative to the center of the earth, respectively; φ is the geocentric latitude of the satellite's location, ω e is the angular velocity vector of the Earth.

[0084] Furthermore, the state space equation construction module specifically performs the following steps:

[0085] The mathematical models of the J2 perturbation force and the atmospheric drag perturbation force are integrated into the general form of the CW equation as the state space equation of the study. The general form of the CW equation is expressed as:

[0086]

[0087] Where ρ = r c -r s =[xyz] T is the position coordinate vector of the companion spacecraft relative to the main spacecraft, is the acceleration coordinate vector of the companion spacecraft relative to the main spacecraft, Δf x , Δf y , Δf z are the three-axis components of the acceleration caused by all the perturbations and control forces on the companion spacecraft relative to the main spacecraft in the relative coordinate system except the gravity of the center of the earth, and ε is the average orbital velocity;

[0088] Add the perturbation force to the general CW equation, that is, write the mathematical model of the J2 perturbation force and the mathematical model of the atmospheric drag perturbation force into the form of acceleration and bring it into Δf x , Δf y , Δf z Among the three items, considering the mathematical model of the J2 perturbation force and the mathematical model of the atmospheric drag perturbation force, the acceleration generated by the non-spherical perturbation J2 item is as follows:

[0089] J2(r)=-(3 / 2)(J2μR e 2 / r 4 )[(1-3sin 2 isin 2 θ)x

[0090] +(2sin 2 isinθcosθ)y+(2sinicosisinθ)z]

[0091] Where r is the position vector of the satellite, J2 is the second spherical harmonic term of the Earth's gravitational field potential, x, y, z represent the spacecraft coordinates in the LVLH coordinate system, and R e is the mean equatorial radius of the Earth, i is the orbital inclination, θ is the ascending node angular distance, and x, y, z are the unit direction vectors of the three axes of the corresponding coordinate system;

[0092] Substituting J2(r) into the general form of the CW equation and linearizing it, we get the improved relative motion dynamics equation, which is expressed as:

[0093]

[0094] Among them, r ref is the constant radius of the circular reference orbit, i ref is the initial orbital inclination of the reference spacecraft, t is a given time, and the parameters c, k, l, q are used to correct the circular reference orbit period caused by the nodal drift caused by the J2 effect:

[0095]

[0096] Φ0=cos -1 [cosi c cosi s +sini c sini s cosΔΩ0]

[0097]

[0098] Among them, s, γ0, Φ0, ΔΩ0 are intermediate variables, is the first-order derivative of Δz0, Δz0 is the difference in the z-axis coordinates between the main spacecraft and the companion spacecraft in the LVLH system in the initial state, i s is the initial orbital inclination of the main spacecraft. Since the main spacecraft is initially located in a circular reference orbit, i s Set to i ref ;i c is the initial orbital inclination of the reference spacecraft;

[0099] Based on the improved relative motion dynamics equation, the acceleration caused by atmospheric drag is substituted and expressed as:

[0100]

[0101] Among them, f Drag,x 、f Drag,y 、f Drag,z They represent the additional acceleration due to aerodynamic drag in the radial direction, in the orbit, and in the normal direction of the orbit respectively; the simultaneous equations Using a circular reference orbit, the velocity of the spacecraft relative to the rotating atmosphere is given by:

[0102]

[0103] in, are the first-order derivatives of x, y, and z respectively, Substituting into the equation Then we get the following state space equation with the perturbation force added:

[0104]

[0105] in, express Calculate in the x-axis direction, express Calculate in the y-axis direction, express The calculation is performed in the z-axis direction.

[0106] Furthermore, the formation configuration parameters of the coplanar flyby natural formation configuration include: the short radius p of the flyby ellipse, the motion amplitude S in the direction perpendicular to the orbital plane, the initial phase difference α, the initial phase σ of the flyby ellipse, and the distance l along the track from the center of the flyby ellipse;

[0107] p=ae A , s = aΔi, α = σ-ψ, l = aΔλ, obtained by the following expression:

[0108]

[0109] The mapping expression of the orbiting star eccentricity is as follows:

[0110]

[0111] where arctan is the inverse sine function, ψ and Δλ are intermediate variables; for a circular reference orbit, e ref =0, then The mapping expression of the orbital inclination of the orbit around the star is as follows:

[0112]

[0113] The mapping expression of the right ascension of the ascending node of the orbiting star is Ω cir =Ω ref +ΔΩ;

[0114] The mapping expression of the perigee angle of the orbiting star is: The mapping expression of the orbiting star mean anomaly is:

[0115]

[0116] Among them, a is the semi-major axis of the orbit, indicating the size of the orbit, and e A represents the eccentricity of the two satellites relative to the ascending node, e cir and e ref are the orbital eccentricities of the orbiting star and the primary star, ω cir and ω ref are the perigee arguments of the orbiting star and the primary star, M cir and M ref are the mean anomalies of the orbiting star and the primary star, i cir is the orbital inclination of the star, Ω cir and Ω ref are the right ascension of the ascending node of the orbiting star and the primary star respectively, ψ and Δλ are intermediate parameters, ΔΩ is the difference in the right ascension of the ascending node of the two satellite orbits, Δi is the angle between the two satellite orbits, and arccos is the inverse cosine function.

[0117] Furthermore, the sliding mode control module specifically performs the following steps:

[0118] The state space equation can be written as The form is expressed as:

[0119]

[0120] where u x ,u y ,u z is the control input of the three axes, t is the simulation time,

[0121] The above equation can be written as a second-order nonlinear dynamic system in the following form:

[0122]

[0123] We can get:

[0124]

[0125] The relative position error and relative velocity error of the satellite are tracked, so x = [x1 x2] T From e = [e1e2] T Substitute, and e1 is the position error, e2 is the speed error, then the sliding surface is expressed as The non-singular terminal sliding mode controller u is designed as follows:

[0126]

[0127] Where f(x) is the non-perturbation term in the state space equation, g(x) represents the perturbation force, b(x) is the gain matrix, β, e = [e1 e2] T , l gand η are both adjustment parameters in sliding mode control, and β is 1e -3 , l g ≥g(x), and l g Take 2e -10 , η is 1e -11 .

[0128] The beneficial effects of the present invention are:

[0129] The present application directly compensates the perturbation force on the control channel, and the non-singular terminal sliding mode can handle nonlinear state space equations, so there is no need to linearize the dynamic equations with the perturbation force added again, thus avoiding the inherent error caused by linearization. Therefore, the technical solution of the present application can effectively maintain the configuration of the coplanar flyby formation satellite under the influence of the J2 perturbation and the atmospheric drag perturbation. BRIEF DESCRIPTION OF THE DRAWINGS

[0130] Figure 1 A flowchart of this application;

[0131] Figure 2 This is a schematic diagram of the x-axis absolute error comparison of the three-axis STK simulation data;

[0132] Figure 3 This is a schematic diagram of the y-axis absolute error comparison of the three-axis STK simulation data;

[0133] Figure 4 This is a schematic diagram of the comparison of the absolute error of the z-axis of the three-axis STK simulation data;

[0134] Figure 5 Schematic diagram of the coplanar flyby configuration;

[0135] Figure 6 When comparing the expected value results of the three-axis tracking of the central satellite and the orbiting satellite and the errors between the cases of not applying control and applying control, it is a schematic diagram showing the corresponding results of the PD control method and the simulation conditions of not applying control under the same simulation time and step length;

[0136] Figure 7 When comparing the expected value results of three-axis tracking of the central satellite and the orbiting satellite and the errors between not applying control and applying control, it is a schematic diagram showing the corresponding results of the control method of the present application and the simulation conditions of not applying control under the same simulation time and step size. DETAILED DESCRIPTION

[0137] It should be particularly noted that, in the absence of conflict, the various embodiments disclosed in this application can be combined with each other.

[0138] Specific implementation method 1: A method for maintaining a coplanar fly-around formation configuration under the influence of space perturbations described in this implementation method;

[0139] Step 1: Determine the types of perturbations to be considered: J2 perturbation and atmospheric drag perturbation, and establish mathematical models of the two perturbations;

[0140] Step 2: Based on the mathematical model of the perturbation force, derive the relative motion dynamics equation with the perturbation force added;

[0141] Step 3. Use the six orbital elements method (the six orbital elements are semi-major axis a, eccentricity e, orbital inclination i, ascending node right ascension Ω, perigee argument ω and mean anomaly M) to design a natural formation configuration for coplanar flyby. The formation configuration parameters include: short radius p of the flyby ellipse, motion amplitude S in the direction perpendicular to the orbital plane, initial phase difference α, initial phase θ of the flyby ellipse, and distance l from the center of the flyby ellipse along the track.

[0142] Step 4: Use the fourth-order Runge-Kutta method to perform numerical analysis to solve the position and velocity of the orbiting satellite relative to the primary star in the LVLH coordinate system at each moment under the condition of two-body motion (i.e., not affected by any perturbation force) as the expected coordinates, expressed as The simulation duration is t and the step length is t step ;

[0143] Step 5: Design a non-singular terminal sliding mode controller, taking the position and velocity of the orbiting satellite relative to the primary satellite in the LVLH coordinate system at each moment as the current coordinates, expressed as The sliding surface is expressed as Where e1 is the position error, e2 is the speed error, β>0, p and q (p>q) are positive odd numbers. The controller is designed as follows:

[0144]

[0145] Step 6: Use the fourth-order Runge-Kutta algorithm to perform numerical analysis to solve the control effect of the non-singular terminal sliding mode controller on the satellite configuration maintenance under the condition of J2 perturbation and atmospheric drag perturbation. The simulation time is t and the step size is t. step .

[0146] The present application can be used for a satellite formation system subjected to any space perturbation force, and can simplify the state equation and optimize the final control result in the design of a satellite configuration maintenance method. The method specifically includes the following steps:

[0147] Step 1: Determine the type of perturbation force to be considered. The most commonly considered are J2 perturbation and atmospheric drag perturbation, and establish their mathematical models. For J2 perturbation, which comes from the non-spherical perturbation of the earth, the normalized gravitational potential function can be written as follows:

[0148]

[0149] The first term on the right side of the equation represents the gravity of the center of the earth, and the last two terms represent the non-spherical gravity of the earth. e and μ are the reference radius of the Earth and the Earth's gravitational constant, respectively. is the satellite's geocentric distance, geocentric latitude and geocentric longitude in the International Terrestrial Reference System (ITRS). and is the normalized gravity coefficient, measured by ground or air gravity. n is the order of the selected gravity field model, and m is the degree of the selected gravity field model. Normalized spherical harmonic function It can be solved by the following recursive formula:

[0150]

[0151] The recursive initial value is:

[0152]

[0153] When n=1, m=1 The value of When n=2, m=1 The value of When n=2, m=2 These three values ​​are the initial values ​​of the recursive formula.

[0154] For the Earth, only the J2 term (n = 2, m = 0) has a normalization factor of 10 -3 The coefficients of other terms are generally 10 -7 ~10 -6 The farther the spacecraft is from the Earth, The smaller the value of , the smaller the impact of non-spherical perturbation, especially for high-order terms. When n is large, its impact is smaller, so only J2 perturbation is considered. For atmospheric drag perturbation, it is generally written as follows:

[0155]

[0156] Among them, C D is the drag coefficient, which is 2.2 in this embodiment, The satellite surface mass ratio is defined as the ratio of the cross-sectional area of ​​the satellite perpendicular to the direction of satellite motion to the satellite mass. It is related to the shape, area and mass of the satellite. In this implementation, it is taken as 0.0056m 2 / kg. ρ is the atmospheric density, which is 4.79651e-13kg / m 3 . is the speed of the satellite relative to the atmosphere. Due to the influence of the earth's rotation, it can be expressed as where v and va are the velocities of the satellite and atmosphere relative to the center of the Earth, v a The expression is as follows:

[0157]

[0158] Where X and Y are the velocities of the satellite and atmosphere relative to the center of the Earth, and φ is the geocentric latitude of the satellite’s location.

[0159] Step 2: Based on the mathematical models of J2 perturbation and atmospheric drag perturbation, they are integrated into the general form of the CW equation. As the state space equation under study, the general form of the CW equation is as follows:

[0160]

[0161] Among them, the position coordinate vector of the companion spacecraft relative to the main spacecraft is [xyz] T , Δf x , Δf y , Δf z are the three-axis components of the acceleration caused by all the perturbations and control forces on the companion spacecraft relative to the main spacecraft in the relative coordinate system except the gravity of the center of the earth, and ε is the average orbital velocity. When the companion spacecraft is in a circular orbit, the value of ε is as follows:

[0162]

[0163] in, is the first-order derivative of ε.

[0164] Adding the perturbation force to the general CW equation means writing the perturbation force model established in step 2 in the form of acceleration and substituting it into Δf x , Δf y , Δf z Among the three items.

[0165] Consider the mathematical model of J2 perturbation and atmospheric drag perturbation. The acceleration generated by the non-spherical perturbation J2 term is as follows:

[0166]

[0167] Where r is the position vector of the satellite, J2 is the second spherical harmonic term of the earth's gravitational field potential energy, and its value is 1.082638e -3 ; i is the orbital inclination, which is 0.001; θ is the ascending node angular distance, which is 0; x, y, z are the unit direction vectors of the three axes of the corresponding coordinate system.

[0168] Substituting equation (9) into equation (7) for linearization, the improved relative motion dynamics equation can be obtained as follows:

[0169]

[0170] Where x, y, and z represent the spacecraft coordinates in the LVLH coordinate system; R e represents the average equatorial radius of the earth, which is 6378.1363 km; r ref is the constant radius of the circular reference orbit, which is 6900 km in this embodiment; i ref is the initial orbital inclination of the reference spacecraft, which is taken as 0.001; t is the given time, and the simulation duration is taken as t = 10000s.

[0171] On the basis of formula (10), substituting the acceleration caused by atmospheric drag, its form is as follows:

[0172]

[0173] where f Drag,x 、f Drag,y 、f Drag,z are the additional accelerations due to aerodynamic drag in the radial direction, in the orbit, and in the normal direction of the orbit, respectively. Combining equations (5) and (6), we can get v a The expression is as follows:

[0174]

[0175] For this scheme, using a circular reference orbit, the velocity of the spacecraft relative to the rotating atmosphere is simplified to:

[0176]

[0177] in The term can be removed because it is a circular reference orbit, and the orbital velocity becomes nc, which is the constant orbital velocity of the circular reference orbit.

[0178] Step 3. Based on the formation configuration parameters: the short radius of the flyby ellipse p, the motion amplitude S in the direction perpendicular to the orbital plane, the initial phase difference α, the initial phase θ of the flyby ellipse and the distance l along the track of the center of the flyby ellipse, a coplanar flyby formation configuration is designed. The mapping equation between the primary star (denoted by the subscript ref) and the orbiting star (denoted by the subscript cir) can be written in the form of (Δe, Δi, ΔΩ, Δω, ΔM) = f(p, s, α, θ, l); set the following parameters:

[0179] Use N ref 、N cir Respectively represent the ascending nodes of the main star and the surrounding star; M ref , M cirRespectively represent the mean anomaly angle of the primary star and the orbiting star; A represents the relative ascending node of the two satellites; ΔΩ represents the difference in right ascension of the ascending node of the two satellite orbits, ΔΩ=Ω cir -Ω ref ; Indicates that from N ref to the geocentric angle of A; k represents the angle from N cir The geocentric angle to A; Δi represents the angle between the two satellite orbits; u ref 、u cir They respectively represent the geocentric angle from the relative ascending node of the two orbits to the current position.

[0180] The calculation formulas for orbital configuration parameters are p = ae A , s = aΔi, α = θ-ψ, l = aΔλ, which can be obtained by the following expression:

[0181]

[0182] The mapping expression of the orbiting star eccentricity is as follows:

[0183]

[0184] in:

[0185]

[0186] For the circular reference orbit in this experiment, e ref =0, then

[0187] The mapping expression of the orbital inclination of the orbit around the star is as follows:

[0188]

[0189] The mapping expression of the right ascension of the ascending node of the orbiting star is Ω cir =Ω ref +ΔΩ.

[0190] The mapping expression of the perigee angle of the orbiting star is:

[0191] The mapping expression of the orbiting satellite mean anomaly angle is:

[0192] In the coplanar flyby configuration, the z-axis component of the orbiting star in the reference coordinate system of the primary star is 0; the semi-major axis of the orbit takes the altitude of the geosynchronous orbit, i.e., a = 6900 km. The six orbit numbers of the coplanar flyby formation configuration finally designed are shown in Table 1.

[0193] Step 4: Use the fourth-order Runge-Kutta method to perform numerical analysis and solve the position and velocity of the orbiting satellite relative to the primary star in the LVLH coordinate system at each moment under the condition of two-body motion (i.e., not affected by any perturbation force) as the expected coordinates, expressed as The simulation duration is t=10000s, and the step length is t step =0.1s;

[0194] Step 5: Design a non-singular terminal sliding mode controller, taking the position and velocity of the orbiting satellite relative to the primary satellite in the LVLH coordinate system at each moment as the current coordinates, expressed as The sliding surface is expressed as Where e1 is the position error, expressed as e1 = [x d -xy d -yz d -z] T , e2 is the speed error, expressed as β>0, p and q (p>q) are positive odd numbers;

[0195] For the dynamic equation of relative motion of satellites targeted by this application, the CW equation with perturbation force established in step 2 can be written as In the form of:

[0196]

[0197] where u x ,u y ,u z is the control input of the three axes, and t is the simulation time, which is 10000s.

[0198] Formula (18) can be written as a second-order nonlinear dynamic system in the following form:

[0199]

[0200] We can get:

[0201]

[0202] Since this application needs to track the relative position error and relative velocity error of the satellite, x = [x1 x2] T From e=[e1 e2] T Substitute, and e1 is the position error, e2 is the speed error. Then the sliding surface is expressed as The non-singular terminal sliding mode controller is designed as follows:

[0203]

[0204] Among them, β is 1e- 3 , lg ≥g(x) and l g Take 2e -10 , η is 1e -11 .

[0205] Step 6: Use the fourth-order Runge-Kutta numerical analysis to solve the control effect of the non-singular terminal sliding mode controller on the satellite configuration maintenance under the condition of J2 perturbation and atmospheric drag perturbation. The simulation time is t = 10000s and the step size is t step =0.1s.

[0206] Experiment: The solution of this application is tested in a simulation environment according to the first specific implementation method. The solution of this application can be divided into three parts: dynamic modeling, formation configuration design, and configuration maintenance. In order to verify the accuracy of the dynamic model, the absolute accuracy of the CW equation simulation data and the relative motion equation simulation data under J2 perturbation & atmospheric drag perturbation are compared under the same initial conditions, taking STK data as the benchmark. The simulation time is 20,000 seconds and the step size is 0.1 seconds. From Figure 5 It can be seen that the CW equation with perturbation force is more accurate than the CW equation without perturbation force, and has more advantages in describing the relative motion between satellites. Table 2 shows the numerical results of absolute error comparison. In order to verify the correctness of the formation design results, the following are given: Figure 6 is a schematic diagram of the coplanar fly-around configuration, including a three-dimensional schematic diagram, an xy plane schematic diagram, and a yz plane schematic diagram; in order to verify the feasibility and superiority of the application plan, the comparison Figure 6 and Figure 7 As can be seen from Table 3, the solution of the present application successfully maintains the satellite formation configuration. If only PD control is used or only natural fly-around formation configuration design is performed, there will be a decrease in accuracy and stability. The present application effectively improves the stability and accuracy of the configuration while successfully maintaining the satellite formation configuration with the addition of perturbation force, and has a very high practical value.

[0207] Table 1 Six parameters of satellite orbit in coplanar flyby configuration

[0208]

[0209] Table 2. Absolute error comparison results with STK data

[0210]

[0211]

[0212] Table 3. Comparison of absolute errors between the proposed method and PD results

[0213]

[0214] This application integrates the perturbation force equation into the satellite relative motion dynamics equation, which is not only applicable to the motion equation under the influence of J2 term perturbation and atmospheric drag perturbation, but also to the situation under the influence of a single perturbation force. It is also applicable to other types of perturbations such as solar light pressure perturbation, solar and lunar gravitational perturbation, etc. As long as a relative mathematical model of perturbation force can be established, the model can be written into the space equation, which is a general method for maintaining satellite configuration under perturbation environment.

[0215] This application aims at the problem of maintaining the configuration of space satellite formations, and proposes a method for natural formation design and further maintenance by a controller. First, a natural formation configuration is designed without perturbation force, and the position and velocity of the orbiting satellite relative to the main satellite at each moment are calculated as the expected position and velocity values ​​at the current moment. Subsequently, the current position difference and velocity difference are used to design a sliding mode controller for control. This solves the problem of needing to perform Taylor expansion in the process of designing satellite formations using CW equations with perturbation force, discarding high-order small terms to reconstruct the mathematical model and further obtain an analytical solution, thereby simplifying the solution process.

[0216] In view of the nonlinear and time-varying characteristics of the dynamic equations with the addition of perturbation forces, the present application adopts the characteristics of non-singular terminal sliding mode and introduces nonlinear functions in the sliding hyperplane, thereby ensuring the convergence speed in the process of solving the equations using the fourth-order Runge-Kutta method and avoiding singular phenomena generated in the integration process, thereby ensuring that the control process is optimized while maintaining the satellite formation configuration.

[0217] It should be noted that the specific implementation is only an explanation and description of the technical solution of the present invention, and cannot be used to limit the scope of protection of the rights. Any partial changes made according to the claims and description of the present invention should still fall within the scope of protection of the present invention.

Claims

1. A method for maintaining a coplanar fly-around formation configuration under the influence of spatial perturbations, characterized in that The following steps are involved: Step 1: For the J2 perturbation and the atmospheric drag perturbation, establish the J2 perturbation force mathematical model and the atmospheric drag perturbation force mathematical model respectively; Step 2: Substitute the mathematical model of the J2 perturbation force and the mathematical model of the atmospheric drag perturbation force into the CW equation to obtain the state space equation with the perturbation force added; Step 3: Use the orbit six-element method to design a natural formation configuration for coplanar flyby; Step 4: Based on the natural formation configuration of coplanar flyby, and using the fourth-order Runge-Kutta method, the position and velocity of the orbiting satellite relative to the primary star in the LVLH coordinate system at each moment under the condition of two-body motion are obtained; Step 5: Take the position and velocity of the orbiting satellite relative to the host satellite in the LVLH coordinate system at each moment as the current coordinates, and use the relative motion dynamics equation with perturbation force to design a non-singular terminal sliding mode controller; Step 6: Use a non-singular terminal sliding mode controller to maintain the coplanar flyby formation configuration.

2. The method for maintaining a coplanar fly-around formation configuration under the influence of space perturbations according to claim 1, characterized in that The mathematical model of the J2 perturbation force is expressed as: Normalized spherical harmonics The specific solution is obtained by the following recursive formula: Among them, the first term on the right side of the mathematical model of J2 perturbation dynamics is represents the gravity at the center of the earth, and the last two terms and is the non-spherical gravity of the Earth, R e and μ are the reference radius of the Earth and the Earth's gravitational constant, respectively. is the geocentric distance, geocentric latitude and geocentric longitude of the satellite in the international earth reference system, n is the order of the selected gravity field model, m is the degree of the selected gravity field model, when m≠n, C nm and S nm The earth is described as a chessboard pattern with similar concavities and convexities, which is called the field harmonic term of the earth's non-spherical perturbation; when m = n, C nm and S nm The earth will be described as a convex and concave fan-shaped pattern, which is called the fan harmonic of the earth's non-spherical perturbation; nm and S nm Normalize it, then in δ is an intermediate variable, Indicates the state of the previous iteration step n-1 and the previous iteration step m-1; represents the state of the previous iteration step n-1 and the current condition m, Represents the state of the previous iteration step m-1 and the current condition n, Indicates the state under the current conditions n and m; The mathematical model of atmospheric drag force is expressed as: Among them, C D is the drag coefficient, is the surface mass ratio, ρ is the atmospheric density, v srel is the speed of the satellite relative to the atmosphere, v and v a are the speeds of the satellite and atmosphere relative to the center of the earth, respectively; X and Y are the speeds of the satellite and atmosphere relative to the center of the earth, respectively; φ is the geocentric latitude of the satellite's location, ω e is the angular velocity vector of the Earth.

3. The method for maintaining a coplanar fly-around formation configuration under the influence of space perturbations according to claim 2, characterized in that The specific steps of step 2 are: The mathematical models of the J2 perturbation force and the atmospheric drag perturbation force are integrated into the general form of the CW equation as the state space equation of the study. The general form of the CW equation is expressed as: Where ρ = r c -r s =[xyz] T is the position coordinate vector of the companion spacecraft relative to the main spacecraft, is the acceleration coordinate vector of the companion spacecraft relative to the main spacecraft, Δf x , Δf y , Δf z are the three-axis components of the acceleration caused by all the perturbations and control forces on the companion spacecraft relative to the main spacecraft in the relative coordinate system except the gravity of the center of the earth, and ε is the average orbital velocity; Add the perturbation force to the general CW equation, that is, write the mathematical model of the J2 perturbation force and the mathematical model of the atmospheric drag perturbation force into the form of acceleration and bring it into Δf x , Δf y , Δf z Among the three items, considering the mathematical model of the J2 perturbation force and the mathematical model of the atmospheric drag perturbation force, the acceleration generated by the non-spherical perturbation J2 item is as follows: J2(r)=-(3 / 2)(J2μR e 2 / r 4 )[(1-3sin 2 isin 2 θ)x+(2sin 2 isinθcosθ)y+(2sinicosisinθ)z] Where r is the position vector of the satellite, J2 is the second spherical harmonic term of the Earth's gravitational field potential, x, y, z represent the spacecraft coordinates in the LVLH coordinate system, and R e is the mean equatorial radius of the Earth, i is the orbital inclination, θ is the ascending node angular distance, and x, y, z are the unit direction vectors of the three axes of the corresponding coordinate system; Substituting J2(r) into the general form of the CW equation and linearizing it, the improved relative motion dynamics equation is obtained, which is expressed as: Among them, r ref is the constant radius of the circular reference orbit, i ref is the initial orbital inclination of the reference spacecraft, t is a given time, and the parameters c, k, l, q are used to correct the circular reference orbit period caused by the nodal drift caused by the J2 effect: Φ0=cos -1 [Like this c Like this s +sini c sini s cosΔΩ0] Among them, s, γ0, Φ0, ΔΩ0 are intermediate variables, is the first-order derivative of Δz0, Δz0 is the difference in the z-axis coordinates between the main spacecraft and the companion spacecraft in the LVLH system in the initial state, i s is the initial orbital inclination of the main spacecraft. Since the main spacecraft is initially located in a circular reference orbit, i s Set to i ref ;i c is the initial orbital inclination of the reference spacecraft; Based on the improved relative motion dynamics equation, the acceleration caused by atmospheric drag is substituted and expressed as: Among them, f Drag,x 、f Drag,y 、f Drag,z They represent the additional acceleration due to aerodynamic drag in the radial direction, in the orbit, and in the normal direction of the orbit respectively; the simultaneous equations Using a circular reference orbit, the velocity of the spacecraft relative to the rotating atmosphere is given by: in, are the first-order derivatives of x, y, and z respectively, Substituting into the equation Then we get the following state space equation with the perturbation force added: in, express Calculate in the x-axis direction, express Calculate in the y-axis direction, express The calculation is performed in the z-axis direction.

4. The method for maintaining a coplanar fly-around formation configuration under the influence of space perturbations according to claim 3 is characterized in that The formation configuration parameters of the coplanar flyby natural formation configuration include: the short radius p of the flyby ellipse, the motion amplitude S in the direction perpendicular to the orbital plane, the initial phase difference α, the initial phase σ of the flyby ellipse, and the distance l along the track from the center of the flyby ellipse; p=ae A , s = aΔi, α = σ-ψ, l = aΔλ, obtained by the following expression: The mapping expression of the orbiting star eccentricity is as follows: Among them, arctan is the inverse sine function, ψ and Δλ are intermediate variables; For a circular reference orbit, e ref =0, then The mapping expression of the orbital inclination of the orbit around the star is as follows: The mapping expression of the right ascension of the ascending node of the orbiting star is Ω cir =Ω ref +ΔΩ; The mapping expression of the perigee angle of the orbiting star is: The mapping expression of the orbiting star mean anomaly is: Among them, a is the semi-major axis of the orbit, indicating the size of the orbit, and e A represents the eccentricity of the two satellites relative to the ascending node, e cir and e ref are the orbital eccentricities of the orbiting star and the primary star, ω cir and ω ref are the perigee arguments of the orbiting star and the primary star, M cir and M ref are the mean anomalies of the orbiting star and the primary star, i cir is the orbital inclination of the star, Ω cir and Ω ref are the right ascension of the ascending node of the orbiting star and the primary star respectively, ψ and Δλ are intermediate parameters, ΔΩ is the difference in the right ascension of the ascending node of the two satellite orbits, Δi is the angle between the two satellite orbits, and arccos is the inverse cosine function.

5. The method for maintaining a coplanar fly-around formation configuration under the influence of space perturbations according to claim 4, characterized in that The specific steps of step five are: The state space equation can be written as The form is expressed as: where u x ,u y ,u z is the control input of the three axes, t is the simulation time, and the above equation is written as a second-order nonlinear dynamic system in the following form: We can get: The relative position error and relative velocity error of the satellite are tracked, so x = [x1 x2] T From e=[e1 e2] T Substitute, and e1 is the position error, e2 is the speed error, then the sliding surface is expressed as The non-singular terminal sliding mode controller u is designed as follows: Where f(x) is the non-perturbation term in the state space equation, g(x) represents the perturbation force, b(x) is the gain matrix, β, e = [e1 e2] T , l g and η are both adjustment parameters in sliding mode control, and β is 1e -3 , l g ≥g(x), and l g Take 2e -10 , η is 1e -11 .

6. A coplanar fly-around formation configuration maintenance system under the influence of space perturbations, characterized in that The system includes a mathematical model building module, a state space equation building module, a formation configuration building module, a two-body motion module and a sliding mode control module; The mathematical model building module is used to establish a J2 perturbation force mathematical model and an atmospheric drag perturbation force mathematical model respectively for the J2 perturbation and the atmospheric drag perturbation; The state space equation construction module is used to bring the mathematical model of the J2 term perturbation force and the mathematical model of the atmospheric drag perturbation force into the CW equation to obtain the state space equation with the perturbation force added; The formation configuration building module is used to design a natural formation configuration for coplanar flyby using the orbit six-element method; The two-body motion module is used to obtain the position and velocity of the orbiting satellite relative to the primary satellite in the LVLH coordinate system at each moment under the two-body motion condition based on the natural formation configuration of coplanar flyby and the fourth-order Runge-Kutta method; The sliding mode control module is used to take the position and velocity of the orbiting satellite relative to the main satellite in the LVLH coordinate system at each moment as the current coordinates, and use the relative motion dynamics equation with perturbation force to design a non-singular terminal sliding mode controller, and finally use the non-singular terminal sliding mode controller to complete the maintenance of the coplanar flyby formation configuration.

7. The system for maintaining a coplanar fly-around formation configuration under the influence of space perturbations according to claim 6, characterized in that The mathematical model of the J2 perturbation force is expressed as: Normalized spherical harmonics The specific solution is obtained by the following recursive formula: The first term on the right side of the equation is represents the gravity at the center of the earth, and the last two terms and is the non-spherical gravity of the Earth, R e and μ are the reference radius of the Earth and the Earth's gravitational constant, respectively. is the geocentric distance, geocentric latitude and geocentric longitude of the satellite in the international earth reference system, n is the order of the selected gravity field model, m is the degree of the selected gravity field model, when m≠n, C nm and S nm The earth is described as a chessboard pattern with similar concavities and convexities, which is called the field harmonic term of the earth's non-spherical perturbation; when m = n, C nm and S nm The earth will be described as a convex and concave fan-shaped pattern, which is called the fan harmonic of the earth's non-spherical perturbation; nm and S nm Normalize it, then in δ is an intermediate variable, Indicates the state of the previous iteration step n-1 and the previous iteration step m-1; represents the state of the previous iteration step n-1 and the current condition m, Represents the state of the previous iteration step m-1 and the current condition n, Indicates the state under the current conditions n and m; The mathematical model of atmospheric drag force is expressed as: Among them, C D is the drag coefficient, is the surface-to-mass ratio, ρ is the atmospheric density, is the speed of the satellite relative to the atmosphere, v and v a are the speeds of the satellite and atmosphere relative to the center of the earth, respectively; X and Y are the speeds of the satellite and atmosphere relative to the center of the earth, respectively; φ is the geocentric latitude of the satellite's location, ω e is the angular velocity vector of the Earth.

8. The system for maintaining a coplanar fly-around formation configuration under the influence of space perturbations according to claim 7, characterized in that The state space equation construction module specifically performs the following steps: The mathematical models of the J2 perturbation force and the atmospheric drag perturbation force are integrated into the general form of the CW equation as the state space equation of the study. The general form of the CW equation is expressed as: Where ρ = r c -r s =[xyz] T is the position coordinate vector of the companion spacecraft relative to the main spacecraft, is the acceleration coordinate vector of the companion spacecraft relative to the main spacecraft, Δf x , Δf y , Δf z are the three-axis components of the acceleration caused by all the perturbations and control forces on the companion spacecraft relative to the main spacecraft in the relative coordinate system except the gravity of the center of the earth, and ε is the average orbital velocity; Add the perturbation force to the general CW equation, that is, write the mathematical model of the J2 perturbation force and the mathematical model of the atmospheric drag perturbation force into the form of acceleration and bring it into Δf x , Δf y , Δf z Among the three items, considering the mathematical model of the J2 perturbation force and the mathematical model of the atmospheric drag perturbation force, the acceleration generated by the non-spherical perturbation J2 item is as follows: J2(r)=-(3 / 2)(J2μR e 2 / r 4 )[(1-3sin 2 isin 2 θ)x+(2sin 2 isinθcosθ)y+(2sinicosisinθ)z] Where r is the position vector of the satellite, J2 is the second spherical harmonic term of the Earth's gravitational field potential, x, y, z represent the spacecraft coordinates in the LVLH coordinate system, and R e is the mean equatorial radius of the Earth, i is the orbital inclination, θ is the ascending node angular distance, and x, y, z are the unit direction vectors of the three axes of the corresponding coordinate system; Substituting J2(r) into the general form of the CW equation and linearizing it, the improved relative motion dynamics equation is obtained, which is expressed as: Among them, r ref is the constant radius of the circular reference orbit, i ref is the initial orbital inclination of the reference spacecraft, t is a given time, and the parameters c, k, l, q are used to correct the circular reference orbit period caused by the nodal drift caused by the J2 effect: Φ0=cos -1 [Like this c Like this s +sini c sini s cosΔΩ0] Among them, s, γ0, Φ0, ΔΩ0 are intermediate variables, is the first-order derivative of Δz0, Δz0 is the difference in the z-axis coordinates between the main spacecraft and the companion spacecraft in the LVLH system in the initial state, i s is the initial orbital inclination of the main spacecraft. Since the main spacecraft is initially located in a circular reference orbit, i s Set to i ref ;i c is the initial orbital inclination of the reference spacecraft; Based on the improved relative motion dynamics equation, the acceleration caused by atmospheric drag is substituted and expressed as: Among them, f Drag,x 、f Drag,y 、f Drag,z They represent the additional acceleration due to aerodynamic drag in the radial direction, in the orbit, and in the normal direction of the orbit respectively; the simultaneous equations Using a circular reference orbit, the velocity of the spacecraft relative to the rotating atmosphere is given by: in, are the first-order derivatives of x, y, and z respectively, Substituting into the equation Then we get the following state space equation with the perturbation force added: in, express Calculate in the x-axis direction, express Calculate in the y-axis direction, express The calculation is performed in the z-axis direction.

9. A coplanar fly-around formation configuration maintenance system under the influence of space perturbations according to claim 8, characterized in that The formation configuration parameters of the coplanar flyby natural formation configuration include: the short radius p of the flyby ellipse, the motion amplitude S in the direction perpendicular to the orbital plane, the initial phase difference α, the initial phase σ of the flyby ellipse, and the distance l along the track from the center of the flyby ellipse; p=ae A , s = aΔi, α = σ-ψ, l = aΔλ, obtained by the following expression: The mapping expression of the orbiting star eccentricity is as follows: Among them, arctan is the inverse sine function, ψ and Δλ are intermediate variables; For a circular reference orbit, e ref =0, then The mapping expression of the orbital inclination of the orbit around the star is as follows: The mapping expression of the right ascension of the ascending node of the orbiting star is Ω cir =Ω ref +ΔΩ; The mapping expression of the perigee angle of the orbiting star is: The mapping expression of the orbiting star mean anomaly is: Among them, a is the semi-major axis of the orbit, indicating the size of the orbit, and e A represents the eccentricity of the two satellites relative to the ascending node, e cir and e ref are the orbital eccentricities of the orbiting star and the primary star, ω cir and ω ref are the perigee arguments of the orbiting star and the primary star, M cir and M ref are the mean anomalies of the orbiting star and the primary star, i cir is the orbital inclination of the star, Ω cir and Ω ref are the right ascension of the ascending node of the orbiting star and the primary star respectively, ψ and Δλ are intermediate parameters, ΔΩ is the difference in the right ascension of the ascending node of the two satellite orbits, Δi is the angle between the two satellite orbits, and arccos is the inverse cosine function.

10. The system for maintaining a coplanar fly-around formation configuration under the influence of space perturbations according to claim 9, characterized in that The sliding mode control module specifically performs the following steps: The state space equation can be written as The form is expressed as: where u x ,u y ,u z is the control input of the three axes, t is the simulation time, The above equation can be written as a second-order nonlinear dynamic system in the following form: We can get: The relative position error and relative velocity error of the satellite are tracked, so x = [x1 x2] T From e=[e1 e2] T Substitute, and e1 is the position error, e2 is the speed error, then the sliding surface is expressed as The non-singular terminal sliding mode controller u is designed as follows: Where f(x) is the non-perturbation term in the state space equation, g(x) represents the perturbation force, b(x) is the gain matrix, β, e = [e1 e2] T , l g and η are both adjustment parameters in sliding mode control, and β is 1e -3 , l g ≥g(x), and l g Take 2e -10 , η is 1e -11 .

Images