A method and system for maintaining coplanar flying formation configuration under the influence of space perturbations

Abstract

A method and system for maintaining the coplanar flying formation configuration under the influence of space perturbations is disclosed, relating to the field of space satellite formation and configuration maintenance technology. Addressing the problem that traditional satellite formation configuration maintenance methods struggle to maintain the configuration of coplanar flying formation satellites under the influence of J2 term perturbations and atmospheric drag perturbations, this application directly compensates for the perturbation force in the control channel. Furthermore, the non-singular terminal sliding mode can handle nonlinear state-space equations, thus eliminating the need for further linearization of the dynamic equations incorporating the perturbation force, avoiding inherent errors caused by linearization. Therefore, the technical solution of this application can effectively maintain the configuration of coplanar flying formation satellites under the influence of J2 term perturbations and atmospheric drag perturbations.

Description

Technical Field

[0001] This invention relates to the field of space satellite formation and configuration maintenance technology, specifically to a method and system for maintaining the configuration of a coplanar orbiting formation under the influence of space perturbations. Background Technology

[0002] Today, space mission requirements are becoming increasingly complex and diverse. Traditional single-satellite systems struggle to achieve more powerful functions within limited payload space, while constellation systems suffer from drawbacks such as poor flexibility, high cost, and system complexity. To address this issue, it has become more important to replace traditional solutions with flight formations composed of multiple satellites that operate through precise measurements and coordinated operations. Achieving mission requirements and maintaining a stable long-term formation structure requires overcoming numerous challenges. After launch, satellite formations need to enter and individually adjust their orbits to form the formation. Satellite formation flight must consider the relative motion and Doppler effect of satellites in space. Furthermore, Earth-orbiting satellite formations must also consider orbital constraints. Therefore, satellite formations are more complex in terms of control, navigation, and data processing. Satellites operating in orbit for extended periods are also subject to various perturbations, such as non-spherical perturbations, atmospheric drag perturbations, solar radiation pressure perturbations, and gravitational perturbations from the Sun and Moon. Even in well-designed natural formation configurations, the member satellites will experience orbital drift due to the presence of perturbations. To successfully complete satellite Earth observation, communication, and navigation missions and avoid collisions between formation satellites, it is necessary to maintain their orbits, ensuring that the motion parameters at each point in the orbit change as little as possible compared to the previous cycle. Traditional methods for maintaining satellite formation configuration often require solving the analytical solutions to the relative motion equations between satellites and then using these solutions to derive the orbital elements and motion equations of the orbiting satellites. However, solving the analytical solutions to the relative motion equations with perturbations is extremely complex, requiring linearization of the perturbation terms, which inevitably introduces deviations into the final results. Therefore, traditional methods for maintaining satellite formation configuration are insufficient for maintaining the configuration of coplanar orbiting formation satellites affected by J2 term perturbations and atmospheric drag perturbations. Summary of the Invention

[0003] The purpose of this invention is to address the problem that traditional satellite formation configuration maintenance methods are unable to maintain the configuration of coplanar flying formation satellites under the influence of J2 term perturbation and atmospheric drag perturbation, and to propose a method and system for maintaining the coplanar flying formation configuration under the influence of space perturbation.

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

[0005] A method for maintaining coplanar flying formation configuration under the influence of space perturbations includes the following steps:

[0006] Step 1: Establish mathematical models for the J2 perturbation and atmospheric drag perturbation, respectively;

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

[0008] Step 3: Design a natural formation configuration for coplanar orbital flight using the orbital six-root method;

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

[0010] Step 5: Use the position and velocity of the orbiting satellite relative to the host star in the LVLH coordinate system at each moment as the current coordinates, and design a non-singular terminal sliding mode controller using the relative motion dynamics equations with added perturbation.

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

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

[0013]

[0014] Normalized spherical harmonic function The solution is obtained using the following recursive formula:

[0015]

[0016] Among them, the first term on the right side of the perturbation mathematical model of term J2. The last two terms represent the Earth's central gravity. For Earth's non-spherical gravity, R e μ and μ are the Earth's reference radius and the Earth's gravitational constant, respectively. Let C be the geocentric distance, geocentric latitude, and geocentric longitude of the satellite in the International Earth Reference System; n be the order of the selected gravity field model; and m be the degree of the selected gravity field model. When m ≠ n, C nm and S nm Describing the Earth as a checkerboard pattern with similar concavity and convexity is called the harmonic term of the Earth's nonspherical perturbation; when m = n, C nm and S nm The Earth is described as a sector pattern of convex and concave shapes, known as the sector harmonic term of the Earth's nonspherical perturbation; for C nm and S nm If normalization is performed, then in δ is an intermediate variable. This represents the state in the previous iteration step n-1 and the previous iteration step m-1; This represents the state under the previous iteration step n-1 and the current condition m. This represents the state under the previous iteration step m-1 and the current condition n. This represents the state under the current conditions n and m;

[0017] The mathematical model for atmospheric drag perturbation is expressed as follows:

[0018]

[0019]

[0020]

[0021] Among them, C D The drag coefficient, Where ρ is the surface-to-mass ratio and ρ is the atmospheric density. Let v be the satellite's velocity relative to the atmosphere, and v be the velocity of the satellite. a Let X and Y be the velocities of the satellite and atmosphere relative to the Earth's center, respectively, and φ be the geocentric latitude of the satellite's location. e This is the Earth's angular velocity vector.

[0022] Furthermore, the specific steps of step two are as follows:

[0023] The J2 perturbation mathematical model and the atmospheric drag perturbation mathematical model are integrated into the general form of the CW equations as the state-space equations under study. The general form of the CW equations is expressed as:

[0024]

[0025] Where ρ=r c -r s =[xyz] T Let be the position coordinate vector of the accompanying spacecraft relative to the host spacecraft. Let Δf be the acceleration coordinate vector of the accompanying spacecraft relative to the host spacecraft. x , Δf y , Δf z Let ε be the three-axis component of the acceleration of the accompanying spacecraft relative to the host spacecraft in the relative coordinate system, caused by all perturbations and control forces other than the gravitational force at the Earth's center, and let ε be the average orbital speed.

[0026] To incorporate perturbation forces into the general CW equations, we write the J2 term perturbation mathematical model and the atmospheric drag perturbation mathematical model in the form of accelerations and substitute them into Δf. x , Δfy , Δf z Of the three terms, considering the mathematical models of the perturbation force J2 and atmospheric drag, the acceleration generated by the non-spherical perturbation J2 term 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 satellite's position vector, J2 is the second spherical harmonic term of the Earth's gravitational potential energy, x, y, z represent the spacecraft coordinates in the LVLH coordinate system, and R e Let θ be the Earth's average equatorial radius, i be the orbital inclination, θ be the ascending angular distance, and x, y, z be the unit direction vectors of the three axes of the corresponding coordinate system.

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

[0031]

[0032] Where, r ref Let i be the constant radius of the circular reference orbit. ref The initial orbital inclination of the reference spacecraft, t being a given time, and parameters c, k, l, and q being used to correct for nodal drift caused by the J2 effect, are used for the circular reference orbit period.

[0033]

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

[0035]

[0036] Among them, s, γ0, Φ0, and ΔΩ0 are intermediate variables. Let be the first derivative of Δz0, where Δz0 is the difference between the z-axis coordinates of the primary spacecraft and the accompanying spacecraft in the LVLH system in the initial state. sThe initial orbital inclination of the main spacecraft is given by i. Since the main spacecraft is initially in a circular reference orbit, i is... s Set to i ref i c The initial orbital inclination of the reference spacecraft;

[0037] Based on the improved equations of relative motion, substituting the acceleration caused by atmospheric drag, it can be expressed as:

[0038]

[0039] Among them, f Drag,x f Drag,y f Drag,z Represent the additional accelerations caused by aerodynamic drag in the radial, in-orbit, and orbital normal directions, respectively; solve the simultaneous equations. Using a circular reference orbit, the velocity expression of the spacecraft relative to the rotating atmosphere is as follows:

[0040]

[0041] in, Let x, y, and z be the first derivatives, respectively. Substitute into the equation The following state-space equations with perturbation force were then obtained:

[0042]

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

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

[0045] p = ae A , s=aΔi, α=σ-ψ, l=aΔλ, can be obtained from the following expressions:

[0046]

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

[0048]

[0049] Where arctan is the arcsine function, and ψ and Δλ are intermediate variables;

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

[0051] The expression for the orbital inclination mapping of a star is as follows:

[0052]

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

[0054] The perigee argument mapping expression for an orbiting star is as follows:

[0055] The expression for the orbital near-point angle mapping is:

[0056]

[0057] Where 'a' is the semi-major axis of the orbit, representing the size of the orbit, and 'e' is the semi-major axis of the orbit. A e represents the eccentricity at the relative ascending node of two satellites. cir and e ref The orbital eccentricities of the surrounding star and the primary star, respectively, ω cir and ω ref M represents the perigee argument of the orbiting star and the primary star, respectively. cir and M ref The mean anomalies of the orbiting star and the primary star, respectively, i cir Ω represents the orbital inclination around the star. cir and Ω ref Δλ and ΔΩ are the right ascensions of the ascending nodes of the orbiting star and the primary star, respectively. ψ and Δλ are intermediate parameters. ΔΩ is the difference in right ascensions of the ascending nodes of the two satellite orbits. Δi is the angle between the two satellite orbits. arccos is the inverse cosine function.

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

[0059] Write the state-space equation as The form is expressed as:

[0060]

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

[0062]

[0063] We can obtain:

[0064]

[0065] Tracking is performed based on the satellite's relative position and relative velocity errors; therefore, x = [x1 x2] T From e = [e1e2] T If we substitute, and e1 is the positional error and e2 is the velocity error, then the sliding surface is represented as follows: The non-singular terminal sliding mode controller u is designed as follows:

[0066]

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

[0068] A coplanar flying formation maintenance system under the influence of spatial perturbation, the system comprising a mathematical model construction module, a state space equation construction module, a formation configuration construction module, a two-body motion module, and a sliding mode control module;

[0069] The mathematical model construction module is used to establish mathematical models of J2 perturbation force and atmospheric drag perturbation force, respectively, for J2 perturbation and atmospheric drag perturbation.

[0070] The state-space equation construction module is used to substitute the J2 term perturbation mathematical model and the atmospheric drag perturbation mathematical model into the CW equation to obtain the state-space equation with perturbation added.

[0071] The formation configuration construction module is used to design coplanar flying natural formation configurations using the orbital six-root method;

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

[0073] The sliding mode control module is used to take the position and velocity of the orbiting satellite relative to the host star in the LVLH coordinate system at each moment as the current coordinates, and to design a non-singular terminal sliding mode controller using the relative motion dynamics equations with added perturbation. Finally, the non-singular terminal sliding mode controller is used to maintain the coplanar flying formation configuration.

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

[0075]

[0076] Normalized spherical harmonic function The solution is obtained using the following recursive formula:

[0077]

[0078] Among them, the first term on the right side of the equation The last two terms represent the Earth's central gravity. and For Earth's non-spherical gravity, R e μ and μ are the Earth's reference radius and the Earth's gravitational constant, respectively. Let C be the geocentric distance, geocentric latitude, and geocentric longitude of the satellite in the International Earth Reference System; n be the order of the selected gravity field model; and m be the degree of the selected gravity field model. When m ≠ n, C nm and S nm Describing the Earth as a checkerboard pattern with similar concavity and convexity is called the harmonic term of the Earth's nonspherical perturbation; when m = n, C nm and S nm The Earth is described as a sector pattern of convex and concave shapes, known as the sector harmonic term of the Earth's nonspherical perturbation; for C nm and S nm If normalization is performed, then in δ is an intermediate variable. This represents the state in the previous iteration step n-1 and the previous iteration step m-1; This represents the state under the previous iteration step n-1 and the current condition m. This represents the state under the previous iteration step m-1 and the current condition n. This represents the state under the current conditions n and m;

[0079] The mathematical model for atmospheric drag perturbation is expressed as follows:

[0080]

[0081]

[0082]

[0083] Among them, C D The drag coefficient, Where ρ is the surface-to-mass ratio and ρ is the atmospheric density. Let v be the satellite's velocity relative to the atmosphere, and v be the velocity of the satellite. a Let X and Y be the velocities of the satellite and atmosphere relative to the Earth's center, respectively, and φ be the geocentric latitude of the satellite's location. e This is the Earth's angular velocity vector.

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

[0085] The J2 perturbation mathematical model and the atmospheric drag perturbation mathematical model are integrated into the general form of the CW equations as the state-space equations under study. The general form of the CW equations is expressed as:

[0086]

[0087] Where ρ=r c -r s =[xyz] T Let be the position coordinate vector of the accompanying spacecraft relative to the host spacecraft. Let Δf be the acceleration coordinate vector of the accompanying spacecraft relative to the host spacecraft. x , Δf y , Δf z Let ε be the three-axis component of the acceleration of the accompanying spacecraft relative to the host spacecraft in the relative coordinate system, caused by all perturbations and control forces other than the gravitational force at the Earth's center, and let ε be the average orbital speed.

[0088] To incorporate perturbation forces into the general CW equations, we write the J2 term perturbation mathematical model and the atmospheric drag perturbation mathematical model in the form of accelerations and substitute them into Δf. x , Δf y , Δf z Of the three terms, considering the mathematical models of the perturbation force J2 and atmospheric drag, the acceleration generated by the non-spherical perturbation J2 term 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 satellite's position vector, J2 is the second spherical harmonic term of the Earth's gravitational potential energy, x, y, z represent the spacecraft coordinates in the LVLH coordinate system, and R e Let θ be the Earth's average equatorial radius, i be the orbital inclination, θ be the ascending angular distance, and x, y, z be the unit direction vectors of the three axes of the corresponding coordinate system.

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

[0093]

[0094] Where, r ref Let i be the constant radius of the circular reference orbit. ref The initial orbital inclination of the reference spacecraft, t being a given time, and parameters c, k, l, and q being used to correct for nodal drift caused by the J2 effect, are used for the circular reference orbit period.

[0095]

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

[0097]

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

[0099] Based on the improved equations of relative motion, substituting the acceleration caused by atmospheric drag, it can be expressed as:

[0100]

[0101] Among them, f Drag,x f Drag,y f Drag,z Represent the additional accelerations caused by aerodynamic drag in the radial, in-orbit, and orbital normal directions, respectively; solve the simultaneous equations. Using a circular reference orbit, the velocity expression of the spacecraft relative to the rotating atmosphere is as follows:

[0102]

[0103] in, Let x, y, and z be the first derivatives, respectively. Substitute into the equation The following state-space equations with perturbation force were then obtained:

[0104]

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

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

[0107] p = ae A , s=aΔi, α=σ-ψ, l=aΔλ, can be obtained from the following expressions:

[0108]

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

[0110]

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

[0112]

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

[0114] The perigee argument mapping expression for an orbiting star is as follows: The expression for the orbital near-point angle mapping is:

[0115]

[0116] Where 'a' is the semi-major axis of the orbit, representing the size of the orbit, and 'e' is the semi-major axis of the orbit. A e represents the eccentricity at the relative ascending node of two satellites. cir and e ref The orbital eccentricities of the surrounding star and the primary star, respectively, ω cir and ω ref M represents the perigee argument of the orbiting star and the primary star, respectively. cir and M ref The mean anomalies of the orbiting star and the primary star, respectively, i cir Ω represents the orbital inclination around the star. cir and Ω ref Δλ and ΔΩ are the right ascensions of the ascending nodes of the orbiting star and the primary star, respectively. ψ and Δλ are intermediate parameters. ΔΩ is the difference in right ascensions of the ascending nodes of the two satellite orbits. Δi is the angle between the two satellite orbits. arccos is the inverse cosine function.

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

[0118] Write the state-space equation as The form is expressed as:

[0119]

[0120] Where u x ,u y ,u z The input is the control input for the three axes, and 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 obtain:

[0124]

[0125] Tracking is performed based on the satellite's relative position and relative velocity errors; therefore, x = [x1 x2] T From e = [e1e2] T If we substitute, and e1 is the positional error and e2 is the velocity error, then the sliding surface is represented as follows: The non-singular terminal sliding mode controller u is designed as follows:

[0126]

[0127] Where f(x) is the unperturbed term in the state-space equation, g(x) represents the perturbed disturbance, b(x) is the gain matrix, β, and e = [e1 e2]. T l gBoth η and β are adjustment parameters in sliding mode control, and β is taken as 1e -3 , l g ≥g(x), and l g Take 2e -10 η takes 1e -11 .

[0128] The beneficial effects of this invention are:

[0129] This application directly compensates for the perturbation force in the control channel, and the non-singular terminal sliding mode can handle nonlinear state-space equations. Therefore, it is not necessary to linearize the dynamic equations incorporating the perturbation force again, avoiding the inherent errors caused by linearization. Thus, the technical solution of this application can effectively maintain the configuration of coplanar flying formation satellites affected by J2 term perturbations and atmospheric drag perturbations. Attached Figure Description

[0130] Figure 1 This is a flowchart of the application;

[0131] Figure 2 A schematic diagram showing the comparison of the absolute error of the x-axis in three-axis STK simulation data;

[0132] Figure 3 A schematic diagram showing the comparison of the absolute error of the y-axis in three-axis STK simulation data;

[0133] Figure 4 A schematic diagram showing the comparison of the absolute error of the z-axis in three-axis STK simulation data;

[0134] Figure 5 A schematic diagram of a coplanar flying configuration;

[0135] Figure 6 The diagram shows the corresponding results of the PD control method and the method without control under the same simulation time and step size when comparing the expected value of the three-axis tracking of the central satellite and the orbiting satellite, as well as the error between the control method and the control method.

[0136] Figure 7 This diagram illustrates the corresponding results of the proposed method and the method without control under simulation conditions of the same simulation duration and step size when comparing the expected values ​​of the three-axis tracking of the central satellite and the orbiting satellite, as well as the errors with and without control. Detailed Implementation

[0137] It should be noted that, where there is no conflict, the various embodiments disclosed in this application can be combined with each other.

[0138] Specific Implementation Method 1: A method for maintaining coplanar flying formation configuration under the influence of spatial perturbations as described in this implementation method;

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

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

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

[0142] Step 4: Numerical analysis using the fourth-order Runge-Kutta method to solve for the desired coordinates of the orbiting satellite's position and velocity relative to the host star in the LVLH coordinate system at each moment under the condition of two-body motion (i.e., unaffected by any perturbation). The simulation duration is t, and the step size is t. step ;

[0143] Step 5: Design a non-singular terminal sliding mode controller, using 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, expressed as... Sliding surface is represented as Where e1 is the position error, e2 is the velocity error, β > 0, and p and q (p > q) are positive odd numbers. The controller is designed as follows:

[0144]

[0145] Step Six: Numerical analysis using a fourth-order Runge-Kutta controller is performed to solve the control effect of the non-singular terminal sliding mode controller on satellite configuration maintenance under J2 term perturbation and atmospheric drag perturbation. The simulation duration is t, and the step size is t. step .

[0146] This application can be used for satellite formation systems subjected to any space perturbation, and can simplify the state equations and optimize the final control results in the design of satellite configuration maintenance methods. The method specifically includes the following steps:

[0147] Step 1: Determine the type of perturbation to be considered. Commonly considered types are J2 perturbation and atmospheric drag perturbation, and establish their mathematical models. For J2 perturbation, which originates from the Earth's non-spherical shape, the normalized gravitational potential function can be written in the following form:

[0148]

[0149] In this equation, the first term on the right-hand side represents the Earth's central gravitational force, and the last two terms represent the Earth's non-spherical gravitational force, R. e μ and μ are the Earth's reference radius and the Earth's gravitational constant, respectively. This represents the satellite's geocentric distance, geocentric latitude, and geocentric longitude in the International Earth Reference System (ITRS). and The normalized gravitational coefficients are measured using gravity on the ground or in the air. n is the order of the selected gravitational field model, and m is the degree of the selected gravitational field model. The normalized spherical harmonic function. The solution can be obtained using the following recursive formula:

[0150]

[0151] The initial values ​​for the recursion are:

[0152]

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

[0154] For Earth, only term J2 (n=2, m=0) has a normalization coefficient of 10. -3 The magnitude is [value], and the coefficients of other terms are generally 10. -7 ~10 -6 On the order of magnitude. The farther the spacecraft is from Earth, the more... The smaller the value of , the smaller the effect of non-spherical perturbation, especially for higher-order terms. When n is relatively large, its effect is even smaller; therefore, only the J2 term perturbation is considered. For atmospheric drag perturbation, it is generally written in the following form:

[0155]

[0156] Among them, C D The drag coefficient is set to 2.2 in this implementation scheme. The surface-to-mass ratio, in relation to atmospheric drag perturbations, is defined as the ratio of the satellite's cross-sectional area perpendicular to its direction of motion to its mass. It is related to the satellite's shape, area, and mass; in this implementation scheme, it is taken as 0.0056m. 2 / kg. ρ is the atmospheric density, taken as 4.79651e-13kg / m³. 3 . Let be the satellite's velocity relative to the atmosphere, which, due to the Earth's rotation, can be expressed as... Where v and va v represents the velocity of the satellite and the atmosphere relative to the Earth's center, respectively. a The expression is as follows:

[0157]

[0158] Where X and Y are the velocities of the satellite and the atmosphere relative to the Earth's center, respectively, 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, integrate them into the general form of the CW equations as the state-space equations under study. The general form of the CW equations is as follows:

[0160]

[0161] The position coordinate vector of the accompanying spacecraft relative to the host spacecraft is [xyz]. T , Δf x , Δf y , Δf z Let ε be the three-axis components of the acceleration of the accompanying spacecraft relative to the host spacecraft in the relative coordinate system, resulting from all perturbations and control forces other than the Earth's gravity. Let ε be the average orbital speed. When the accompanying spacecraft is in a circular orbit, the values ​​of ε are as follows:

[0162]

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

[0164] Adding a perturbation force to the general CW equations involves rewriting the perturbation force model established in step two as an expression of acceleration and substituting it into Δf. x , Δf y , Δf z Of the three items.

[0165] Consider the mathematical models 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 satellite's position vector, and J2 is the second spherical harmonic term of the Earth's gravitational potential energy, with a value of 1.082638e. -3 ; i is the orbital inclination angle, with a value of 0.001; θ is the ascending angular distance, with a value of 0; x, y, and z are the unit direction vectors of the three axes of the corresponding coordinate system.

[0168] Substituting equation (9) into equation (7) and linearizing it, we can finally obtain the improved relative motion dynamics equation as follows:

[0169]

[0170] Where x, y, and z represent the spacecraft coordinates in the LVLH coordinate system; R e This represents the Earth's average equatorial radius, taken as 6378.1363 km; r ref The constant radius of the circular reference orbit is taken as 6900 km in this implementation scheme; ref The initial orbital inclination of the reference spacecraft is set to 0.001; t is a given time, and the simulation duration is set to t = 10000s.

[0171] Based on equation (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 Let v represent the additional accelerations due to aerodynamic drag in the radial, in-orbit, and orbital normal directions, respectively. Combining equations (5) and (6), we get v. a The expression is as follows:

[0174]

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

[0176]

[0177] in The term can be removed since it is a circular reference orbit; orbital speed. When nc is converted to nc, it represents the constant orbital speed of the circular reference orbit.

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

[0179] Use N ref N cir M represents the ascending node of the primary star and the orbiting star, respectively; ref M cirLet A and B represent the mean anomalies of the primary star and the orbiting star, respectively; let A represent the relative ascending node of the two satellites; let ΔΩ represent the difference in right ascension of the ascending nodes of the two satellite orbits, ΔΩ = Ω. cir -Ω ref ; Indicates from N ref The geocentric angle from A; k represents the distance from N. cir The geocentric angle to A; Δi represents the angle between the two satellite orbits; u ref u cir These represent the geocentric angles from the relative ascending intersection of the two orbits to the current position.

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

[0181]

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

[0183]

[0184] in:

[0185]

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

[0187] The expression for the orbital inclination mapping of a star is as follows:

[0188]

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

[0190] The perigee argument mapping expression for an orbiting star is as follows:

[0191] The expression for the mapping of the perihelion angle around the star is as follows:

[0192] In the coplanar flying 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 is taken as the altitude of the geosynchronous orbit, i.e., a = 6900 km. The six elements of the final designed coplanar flying formation configuration are shown in Table 1.

[0193] Step 4: Numerical analysis using the fourth-order Runge-Kutta method. Under the condition of two-body motion (i.e., unaffected by any perturbation), the position and velocity of the orbiting satellite relative to the host star in the LVLH coordinate system at each moment are taken as the desired coordinates, expressed as follows: The simulation duration is t = 10000s, and the step size is t. step =0.1s;

[0194] Step 5: Design a non-singular terminal sliding mode controller, using 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, expressed as: Sliding surface is represented as Where e1 is the position error, expressed as e1 = [x d -xy d -yz d -z] T e2 represents the speed error, expressed as β > 0, p and q (p > q) are positive odd numbers;

[0195] Regarding the dynamic equations of the relative motion of the satellite addressed in this application, the CW equations with perturbation established in step two can be written as follows: Format:

[0196]

[0197] Where u x ,u y ,u z The input is the control input for the three axes, and t is the simulation time, which is set to 10000s.

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

[0199]

[0200] We can obtain:

[0201]

[0202] Since this application requires tracking the relative position error and relative velocity error of the satellite, therefore x = [x1 x2] T From e = [e1 e2] T The substitution is performed, and e1 is the position error, and e2 is the velocity error. Then the sliding surface is represented as... The non-singular terminal sliding mode controller is designed as follows:

[0203]

[0204] Where β takes the value of 1e- 3 , lg ≥g(x) and l g Take 2e -10 η takes 1e -11 .

[0205] Step Six: Numerical analysis using a fourth-order Runge-Kutta controller is performed to solve the control effect of the non-singular terminal sliding mode controller on satellite configuration maintenance under J2 term perturbation and atmospheric drag perturbation. The simulation duration is t = 10000 s, and the step size is t. step =0.1s.

[0206] Experiment: The proposed scheme was tested in a simulation environment according to Specific Implementation Method 1. The proposed scheme can be divided into three parts: dynamic modeling, formation configuration design, and configuration maintenance. To verify the accuracy of the dynamic model, using STK data as a benchmark, the absolute accuracy of the CW equation simulation data and the relative motion equation simulation data under J2 perturbation and atmospheric drag perturbation was compared under the same initial conditions. The simulation duration was 20000s, with a step size of 0.1s. From... Figure 5 It can be seen that the CW equations with perturbation are more accurate than those without, and are more advantageous in describing the relative motion between satellites. Table 2 shows the numerical results for comparing absolute errors. To verify the correctness of the formation design results, the following are presented: Figure 6 This is a schematic diagram of a coplanar flying configuration, including a three-dimensional schematic diagram, an xy-plane schematic diagram, and a yz-plane schematic diagram; to verify the feasibility and superiority of the proposed scheme, a comparison is made. Figure 6 and Figure 7 As shown in Table 3, the proposed solution successfully maintains the satellite formation configuration. Using only PD control or only natural fly-around formation design results in reduced accuracy and stability. This application, by successfully maintaining the satellite formation configuration with added perturbation forces, effectively improves the stability and accuracy of the configuration, demonstrating significant practical value.

[0207] Table 1. Orbital parameters of six elements for coplanar flying configuration satellites

[0208]

[0209] Table 2 shows the comparison results of the absolute error with STK data.

[0210]

[0211]

[0212] Table 3 Comparison of absolute errors between the method and PD results in this application

[0213]

[0214] This application integrates the perturbation equations into the satellite's relative motion dynamics equations. It is applicable not only to motion equations influenced by J2 term perturbations and atmospheric drag perturbations, but also to cases influenced by a single perturbation force. It is also applicable to other types of perturbations, such as solar radiation pressure perturbations and lunar gravitational perturbations. As long as a mathematical model of the relative perturbation forces can be established, the model can be incorporated into the space equations, providing a general method for maintaining satellite configuration under perturbation environments.

[0215] This application addresses the problem of maintaining the configuration of space satellite formations by proposing a method that combines natural formation design with a controller for further maintenance. First, a natural formation configuration without perturbations is designed, and the position and velocity of the orbiting star relative to the host star at each moment are calculated as the desired position and velocity values ​​for the current moment. Then, a sliding mode controller is designed using the current position and velocity differences for control. This solves the problem of needing Taylor expansion, discarding higher-order terms to reconstruct the mathematical model, and further obtaining analytical solutions when designing satellite formations using CW equations with perturbations, thus simplifying the solution process.

[0216] This application addresses the nonlinear and time-varying characteristics of the dynamic equations with added perturbation by employing the features of nonsingular terminal sliding mode. A nonlinear function is introduced into the sliding hyperplane, ensuring the convergence speed during equation solving using fourth-order Runge-Kutta and avoiding singularities during integration. This ensures the optimization of the control process while maintaining the satellite formation configuration.

[0217] It should be noted that the specific embodiments are merely explanations and illustrations of the technical solution of the present invention and should not be used to limit the scope of protection. Any modifications made in accordance with the claims and specification of the present invention that are only partial should still fall within the protection scope of the present invention.

Claims

1. A method for maintaining a coplanar flying formation configuration under the influence of spatial perturbations, characterized in that... Includes the following steps: Step 1: Establish mathematical models for the J2 perturbation and atmospheric drag perturbation, respectively; Step 2: Substitute the J2 term perturbation mathematical model and the atmospheric drag perturbation mathematical model into the CW equation to obtain the state-space equation with perturbation added; Step 3: Design a natural formation configuration for coplanar orbital flight using the orbital six-root method; Step 4: Based on the natural formation configuration of coplanar orbiting, and using the fourth-order Runge-Kutta method, obtain the position and velocity of the orbiting satellite relative to the host star in the LVLH coordinate system at each moment under the two-body motion condition; Step 5: Use the position and velocity of the orbiting satellite relative to the host star in the LVLH coordinate system at each moment as the current coordinates, and design a non-singular terminal sliding mode controller using the relative motion dynamics equations with added perturbation. Step 6: Maintain coplanar flying formation configuration using a non-singular terminal sliding mode controller; The mathematical model of the J2 term perturbation force is expressed as follows: Normalized spherical harmonic function The solution is obtained using the following recursive formula: Among them, the first term on the right side of the perturbation mathematical model of term J2. The last two terms represent the Earth's central gravity. and For Earth's non-spherical gravity, R e μ and μ are the Earth's reference radius and the Earth's gravitational constant, respectively. Let C be the geocentric distance, geocentric latitude, and geocentric longitude of the satellite in the International Earth Reference System; n be the order of the selected gravity field model; and m be the degree of the selected gravity field model. When m ≠ n, C nm and S nm Describing the Earth as a checkerboard pattern with similar concavity and convexity is called the harmonic term of the Earth's nonspherical perturbation; when m = n, C nm and S nm The Earth is described as a sector pattern of convex and concave shapes, known as the sector harmonic term of the Earth's nonspherical perturbation; for C nm and S nm If normalization is performed, then in δ is an intermediate variable. This represents the state in the previous iteration step n-1 and the previous iteration step m-1; This represents the state under the previous iteration step n-1 and the current condition m. This represents the state under the previous iteration step m-1 and the current condition n. This represents the state under the current conditions n and m; The mathematical model for atmospheric drag perturbation is expressed as follows: Among them, C D The drag coefficient, Where ρ is the surface-to-mass ratio and ρ is the atmospheric density. v is the satellite's velocity relative to the atmosphere. s and v a Let X and Y be the velocities of the satellite and atmosphere relative to the Earth's center, respectively, and φ be the geocentric latitude of the satellite's location. e This is the Earth's angular velocity vector.

2. The method for maintaining coplanar flying formation configuration under the influence of space perturbation according to claim 1, characterized in that... The specific steps of step two are as follows: The J2 perturbation mathematical model and the atmospheric drag perturbation mathematical model are integrated into the general form of the CW equations as the state-space equations under study. The general form of the CW equations is expressed as: Where ρ=r c -r s =[xyz] T Let be the position coordinate vector of the accompanying spacecraft relative to the host spacecraft. Let Δf be the acceleration coordinate vector of the accompanying spacecraft relative to the host spacecraft. x , Δf y , Δf z Let ε be the three-axis component of the acceleration of the accompanying spacecraft relative to the host spacecraft in the relative coordinate system, caused by all perturbations and control forces other than the gravitational force at the Earth's center, and let ε be the average orbital speed. To incorporate perturbation forces into the general CW equations, we write the J2 term perturbation mathematical model and the atmospheric drag perturbation mathematical model in the form of accelerations and substitute them into Δf. x , Δf y , Δf z Of the three terms, considering the mathematical models of the perturbation force J2 and atmospheric drag, the acceleration generated by the non-spherical perturbation J2 term 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 satellite's position vector, J2 is the second spherical harmonic term of the Earth's gravitational potential energy, x, y, z represent the spacecraft coordinates in the LVLH coordinate system, and R e Let θ be the Earth's average equatorial radius, i be the orbital inclination, θ be the ascending angular distance, and x, y, z be the unit direction vectors of the three axes of the corresponding coordinate system. Substituting J2(r) into the general form of the CW equations and linearizing it, we obtain the improved equations of relative motion dynamics, expressed as: Where, r ref Let i be the constant radius of the circular reference orbit. ref The initial orbital inclination of the reference spacecraft, t being a given time, and parameters c, k, l, and q being used to correct for nodal drift caused by the J2 effect, are used for the circular reference orbit period. Φ0=cos -1 [Like this c Like this s +sini c sins s cosΔΩ0] Among them, s, γ0, Φ0, and ΔΩ0 are intermediate variables. Let be the first derivative of Δz0, where Δz0 is the difference between the z-axis coordinates of the primary spacecraft and the accompanying spacecraft in the LVLH frame in the initial state, and i is the initial orbital inclination of the primary spacecraft. Since the primary spacecraft is initially in a circular reference orbit, i is... s Set to i ref i c The initial orbital inclination of the reference spacecraft; Based on the improved equations of relative motion, substituting the acceleration caused by atmospheric drag, it can be expressed as: Among them, f Drag,x f Drag,y f Drag,z Represent the additional accelerations caused by aerodynamic drag in the radial, in-orbit, and orbital normal directions, respectively; solve the simultaneous equations. Using a circular reference orbit, the velocity expression of the spacecraft relative to the rotating atmosphere is as follows: in, Let x, y, and z be the first derivatives, respectively. Substitute into the equation The following state-space equations with perturbation force were then obtained: in, express Calculations are performed along the x-axis. express Calculations are performed in the y-axis direction. express The calculation is performed in the z-axis direction.

3. The method for maintaining coplanar flying formation configuration under the influence of space perturbations according to claim 2, characterized in that... The formation parameters of the natural formation configuration of the coplanar flight include: the minor radius p of the flight ellipse, the motion amplitude S in the direction perpendicular to the orbital plane, the initial phase difference α, the initial phase σ of the flight ellipse, and the distance l along the track from the center of the flight ellipse. p = ae A , s=aΔi, α=σ-ψ, l=aΔλ, can be obtained from the following expressions: The expression for the eccentricity mapping of an orbiting star is as follows: Where arctan is the arcsine function, and ψ and Δλ are intermediate variables; For a circular reference orbit, e ref =0, then we have The expression for the orbital inclination mapping of a star is as follows: The expression for the right ascension mapping of the ascending node of an orbiting star is Ω. cir =Ω ref +ΔΩ; The perigee argument mapping expression for an orbiting star is as follows: The expression for the orbital near-point angle mapping is: Where 'a' is the semi-major axis of the orbit, representing the size of the orbit, and 'e' is the semi-major axis of the orbit. A e represents the eccentricity at the relative ascending node of two satellites. cir and e ref The orbital eccentricities of the surrounding star and the primary star, respectively, ω cir and ω ref M represents the perigee argument of the orbiting star and the primary star, respectively. cir and M ref The mean anomalies of the orbiting star and the primary star, respectively, i cir Ω represents the orbital inclination around the star. cir and Ω ref Δλ and ΔΩ are the right ascensions of the ascending nodes of the orbiting star and the primary star, respectively. ψ and Δλ are intermediate parameters. ΔΩ is the difference in right ascensions of the ascending nodes of the two satellite orbits. Δi is the angle between the two satellite orbits. arccos is the inverse cosine function.

4. The method for maintaining coplanar flying formation configuration under the influence of space perturbations according to claim 3, characterized in that... The specific steps of step five are as follows: Write the state-space equation as The form is expressed as: Where u x ,u y ,u z Given the control input for the three axes and t as the simulation time, the above equation can be written as a second-order nonlinear dynamic system in the following form: We can obtain: Tracking is performed based on the satellite's relative position and relative velocity errors; therefore, x = [x1 x2] T From e = [e1 e2] T If we substitute, and e1 is the positional error and e2 is the velocity error, then the sliding surface is represented as follows: The non-singular terminal sliding mode controller u is designed as follows: Where f(x) is the unperturbed term in the state-space equation, g(x) represents the perturbed disturbance, b(x) is the gain matrix, β, and e = [e1 e2]. T l g Both η and β are adjustment parameters in sliding mode control, and β is taken as 1e -3 , l g ≥g(x), and l g Take 2e -10 η takes 1e -11 .

5. A coplanar flying formation maintenance system under the influence of spatial perturbations, characterized in that... The system includes a mathematical model construction module, a state-space equation construction module, a formation configuration construction module, a two-body motion module, and a sliding mode control module; The mathematical model construction module is used to establish mathematical models of J2 perturbation force and atmospheric drag perturbation force, respectively, for J2 perturbation and atmospheric drag perturbation. The state-space equation construction module is used to substitute the J2 term perturbation mathematical model and the atmospheric drag perturbation mathematical model into the CW equation to obtain the state-space equation with perturbation added. The formation configuration construction module is used to design coplanar flying natural formation configurations using the orbital six-root method; The two-body motion module is used for natural formation configuration based on coplanar orbiting, and the fourth-order Runge-Kutta method is used to obtain the position and velocity of the orbiting satellite relative to the host star in the LVLH coordinate system at each moment under the two-body motion condition. The sliding mode control module is used to take the position and velocity of the orbiting satellite relative to the host star in the LVLH coordinate system at each moment as the current coordinates, and to design a non-singular terminal sliding mode controller using the relative motion dynamics equation with added perturbation. Finally, the non-singular terminal sliding mode controller is used to maintain the coplanar flying formation configuration. The mathematical model of the J2 term perturbation force is expressed as follows: Normalized spherical harmonic function The solution is obtained using the following recursive formula: Among them, the first term on the right side of the equation The last two terms represent the Earth's central gravity. and For Earth's non-spherical gravity, R e μ and μ are the Earth's reference radius and the Earth's gravitational constant, respectively. Let C be the geocentric distance, geocentric latitude, and geocentric longitude of the satellite in the International Earth Reference System; n be the order of the selected gravity field model; and m be the degree of the selected gravity field model. When m ≠ n, C nm and S nm Describing the Earth as a checkerboard pattern with similar concavity and convexity is called the harmonic term of the Earth's nonspherical perturbation; when m = n, C nm and S nm The Earth is described as a sector pattern of convex and concave shapes, known as the sector harmonic term of the Earth's nonspherical perturbation; for C nm and S nm If normalization is performed, then in δ is an intermediate variable. This represents the state in the previous iteration step n-1 and the previous iteration step m-1; This represents the state under the previous iteration step n-1 and the current condition m. This represents the state under the previous iteration step m-1 and the current condition n. This represents the state under the current conditions n and m; The mathematical model for atmospheric drag perturbation is expressed as follows: Among them, C D The drag coefficient, Where ρ is the surface-to-mass ratio and ρ is the atmospheric density. v is the satellite's velocity relative to the atmosphere. s and v a Let X and Y be the velocities of the satellite and atmosphere relative to the Earth's center, respectively, and φ be the geocentric latitude of the satellite's location. e This is the Earth's angular velocity vector.

6. A coplanar flying formation maintenance system under the influence of space perturbations according to claim 5, characterized in that... The state-space equation construction module specifically performs the following steps: The J2 perturbation mathematical model and the atmospheric drag perturbation mathematical model are integrated into the general form of the CW equations as the state-space equations under study. The general form of the CW equations is expressed as: Where ρ=r c -r s =[xyz] T Let be the position coordinate vector of the accompanying spacecraft relative to the host spacecraft. Let Δf be the acceleration coordinate vector of the accompanying spacecraft relative to the host spacecraft. x , Δf y , Δf z Let ε be the three-axis component of the acceleration of the accompanying spacecraft relative to the host spacecraft in the relative coordinate system, caused by all perturbations and control forces other than the gravitational force at the Earth's center, and let ε be the average orbital speed. To incorporate perturbation forces into the general CW equations, we write the J2 term perturbation mathematical model and the atmospheric drag perturbation mathematical model in the form of accelerations and substitute them into Δf. x , Δf y , Δf z Of the three terms, considering the mathematical models of the perturbation force J2 and atmospheric drag, the acceleration generated by the non-spherical perturbation J2 term 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 satellite's position vector, J2 is the second spherical harmonic term of the Earth's gravitational potential energy, x, y, z represent the spacecraft coordinates in the LVLH coordinate system, and R e Let θ be the Earth's average equatorial radius, i be the orbital inclination, θ be the ascending angular distance, and x, y, z be the unit direction vectors of the three axes of the corresponding coordinate system. Substituting J2(r) into the general form of the CW equations and linearizing it, we obtain the improved equations of relative motion dynamics, expressed as: Where, r ref Let i be the constant radius of the circular reference orbit. ref The initial orbital inclination of the reference spacecraft, t being a given time, and parameters c, k, l, and q being used to correct for nodal drift caused by the J2 effect, are used for the circular reference orbit period. Φ0=cos -1 [Like this c Like this s +sini c sins s cosΔΩ0] Among them, s, γ0, Φ0, and ΔΩ0 are intermediate variables. Let be the first derivative of Δz0, where Δz0 is the difference between the z-axis coordinates of the primary spacecraft and the accompanying spacecraft in the LVLH system in the initial state. s The initial orbital inclination of the main spacecraft is given by i. Since the main spacecraft is initially in a circular reference orbit, i is... s Set to i ref i c The initial orbital inclination of the reference spacecraft; Based on the improved equations of relative motion, substituting the acceleration caused by atmospheric drag, it can be expressed as: Among them, f Drag,x f Drag,y f Drag,z Represent the additional accelerations caused by aerodynamic drag in the radial, in-orbit, and orbital normal directions, respectively; solve the simultaneous equations. Using a circular reference orbit, the velocity expression of the spacecraft relative to the rotating atmosphere is as follows: in, Let x, y, and z be the first derivatives, respectively. Substitute into the equation The following state-space equations with perturbation force were then obtained: in, express Calculations are performed along the x-axis. express Calculations are performed in the y-axis direction. express The calculation is performed in the z-axis direction.

7. A coplanar flying formation maintenance system under the influence of space perturbations according to claim 6, characterized in that... The formation parameters of the natural formation configuration of the coplanar flight include: the minor radius p of the flight ellipse, the motion amplitude S in the direction perpendicular to the orbital plane, the initial phase difference α, the initial phase σ of the flight ellipse, and the distance l along the track from the center of the flight ellipse. p = ae A , s=aΔi, α=σ-ψ, l=aΔλ, can be obtained from the following expressions: The expression for the eccentricity mapping of an orbiting star is as follows: Where arctan is the arcsine function, and ψ and Δλ are intermediate variables; for a circular reference orbit, e ref =0, then we have The expression for the orbital inclination mapping of a star is as follows: The expression for the right ascension mapping of the ascending node of an orbiting star is Ω. cir =Ω ref +ΔΩ; The perigee argument mapping expression for an orbiting star is as follows: The expression for the orbital near-point angle mapping is: Where 'a' is the semi-major axis of the orbit, representing the size of the orbit, and 'e' is the semi-major axis of the orbit. A e represents the eccentricity at the relative ascending node of two satellites. cir and e ref The orbital eccentricities of the surrounding star and the primary star, respectively, ω cir and ω ref M represents the perigee argument of the orbiting star and the primary star, respectively. cir and M ref The mean anomalies of the orbiting star and the primary star, respectively, i cir Ω represents the orbital inclination around the star. cir and Ω ref Δλ and ΔΩ are the right ascensions of the ascending nodes of the orbiting star and the primary star, respectively. ψ and Δλ are intermediate parameters. ΔΩ is the difference in right ascensions of the ascending nodes of the two satellite orbits. Δi is the angle between the two satellite orbits. arccos is the inverse cosine function.

8. A coplanar flying formation maintenance system under the influence of space perturbations according to claim 7, characterized in that... The sliding mode control module specifically performs the following steps: Write the state-space equation as The form is expressed as: Where u x ,u y ,u z The input is the control input for the three axes, and t is the simulation time. The above equation can be written as a second-order nonlinear dynamic system in the following form: We can obtain: Tracking is performed based on the satellite's relative position and relative velocity errors; therefore, x = [x1 x2] T From e = [e1 e2] T If we substitute, and e1 is the positional error and e2 is the velocity error, then the sliding surface is represented as follows: The non-singular terminal sliding mode controller u is designed as follows: Where f(x) is the unperturbed term in the state-space equation, g(x) represents the perturbed disturbance, b(x) is the gain matrix, β, and e = [e1 e2]. T l g Both η and β are adjustment parameters in sliding mode control, and β is taken as 1e -3 , l g ≥g(x), and l g Take 2e -10 η takes 1e -11 .

Images