Relative attitude stability analysis of double satellite formation for gravity field exploration in space debris environment

Spacecraft operating in low orbit are at risk of being hit by space debris. In the debris environment, the impact of debris is likely to cause the double satellite formation to exit science mode or even lead to the divergence of the control system, thus affecting the scientific exploration mission. In this paper, the attitude stability of the double satellite formation for gravity field in the near circular and polar orbit in the space debris environment is studied. Firstly, based on Lyapunov control and LQR, two sets of control models of stochastic collision for two satellites aligned with each other were proposed, and the actuators were modelled and assigned. Secondly, models of collision probability and momentum are developed. The distribution law of space debris is obtained according to the international common debris software. Meanwhile, probability density function of two independent collisions is gained. Finally, through Monte Carlo simulation and statistics, the changes of relative attitude and thrust torque are simulated when the satellite obtains the angular momentum for a short period of time due to being impacted by space debris. During the 10-year mission period, the number of times that the space debris impact makes the satellite attitude out of the science mode and the number of times that the control system diverges are obtained, which provides a reference for the normal manner of the double satellite formation for gravity field exploration.

exploration launched by NASA and GFZ.The satellites in the formation are about 220 km apart in near circular and polar orbits.They communicate with each other through microwave, and attitude pointing accuracy is better than 3 mrad.The actuator is composed of magnetorquers and cold air thrusters 10,11 .In 2018, they jointly launched GRACE-FO with the distance between two satellites is 50 km.The change of double satellite distance was measured by laser interferometry, and the attitude pointing accuracy was higher, reaching 0.24 mrad 12 .The satellites for gravity field exploration from single satellite 13,14 to double satellite formation [10][11][12] and from microwave to laser interferometry requires more and more attitude stability, supporting inversion of higher precision on the Earth's gravity field.Afterwards, ESA and NASA jointly proposed Next Generation Gravity Missions-NGGM with two formations, one with an orbital inclination of 90° and the other of 63°1 5 and China proposed TianQin-2 test satellite.
Since the two satellites of the double satellite formation are identical, the laser beams emitted from the two satellites must be aligned with each other to ensure that the two laser beams can interfere.In other words, the target of the attitude control of the formation is that the outgoing beam of one satellite is aligned with the receiving end of other satellite, and the outgoing beam of the other satellite is aligned with the receiving end of this satellite.If the attitude of one satellite is misaligned, the optical power received on the interferometric signal of other satellite is relatively low, causing inefficient interference efficiency.While the information of angle jitter will be coupled into the distance measurement of the formation, causing measurement bias 16 .Therefore, the maintaining the stability of the relative attitude of the double satellite formation is a key step in ensuring laser interferometric ranging.However, GRACE's satellite development company, Airbus Defense and Space, has only briefly described the control of its attitude 11 , not to mention research on its stability maintenance in the face of complex space environments.Reference 17 developed two control algorithms based on Lyapunov control and LQR on the condition of attitude control accuracy of GRACE and GRACE-FO.Simulation results show that the controller designed by Lyapunov control algorithm has better comprehensive control effect.However, no relevant research has been found on the impact of space debris on the relative attitude stability of the double satellite formation.
Since the control of two satellites in the double satellite formation for gravity field exploration is almost the same, this paper takes one of the double satellite formation as the research object and establishes a satellite relative attitude dynamics model.The disturbance torque in this model includes the gravity gradient torque and torque caused by the difference of tensor of inertia.The control torques of magnetorquers are described, and the thrust model of cold air thruster is established.The serial link control is adopted for the two actuators.Two control algorithms based on Lyapunov and LQR have been cited in the space debris environment.The number of space debris making satellite attitudes exit science mode that caused by impact, as well as making the control system diverged that obtained under certain control accuracy, and their probability of normal manner has been analysed.This article first focuses on and analyzes the normal manner of low orbit complex high-precision formation detector systems in debris environments.With the increasing number of space debris, this issue will become more remarkable.Establishing a set of analysis methods and means for the normal manner ability of detector systems has important reference value for the scientific measurement and in orbit operation management of detectors.This paper uses the control algorithm that meets the task requirements to obtain the critical collision of space debris that causes the control system to diverge.A rich supply of data that meets the law of debris distribution is selected, which conducts a Monte Carlo simulation on them.Creatively combining the control algorithm with a Monte Carlo simulation.The control system diverges due to be impacted by space debris in the current task period is simulated through a large amount of data.
The second part is the description of satellite angle motion.The third part is the design of the formation attitude controller and control assignment of the actuator.The fourth part is the probability and momentum modeling of space debris impact, and the fifth part is the simulation and discussion.

Description of the angular motion
The double satellite formation for the gravity field exploration requires that the two satellites align with each other, and one satellite always points to the other accurately.The double satellite formation is composed of two satellites with almost identical motion and control.Therefore, this paper takes one of the satellites as the research object and controls its body frame to coincide with the reference frame within the error range.To simplify the study, we assume here that the reference frame is known and not affected by space debris impacts.The definition of coordinate frame is as follows.The reference frame of the double satellite formation is shown in Fig. 1a, and the simplified shape of the satellite is shown in Fig. 1b i 1 , i 2 , i 3 and e 1 , e 2 , e 3 are the unit vectors of the x, y, z axis under the orbital and reference frames, respectively, and the unit vector is defined as follows The rotation matrix and its derivatives from the inertial frame to the reference system are The reference angular velocity expressed without collision with space in the body frame is The skew symmetric matrix of cross product is

Lyapunov attitude control
The attitude dynamics equation of the satellite moving around the center of mass is where,J = diag(J 1 , J 2 , J 3 ),ω abs is angular velocity of satellite relative inertial frame, M ctrl ∈ R n×1 is the control torque applied to the satellite, M ext ∈ R n×1 is the external disturbance torque, and M impact ∈ R n×1 is the instan- taneous torque obtained after the satellite is impacted.Q is the rotation matrix from the inertial frame to body frame expressed in terms of Euler angle.The relative angular velocity and angular acceleration from body frame to the reference frame are ω = ω abs − Dω ref and ω = ωabs − D ωref + ω × Dω ref respectively, where D is the rotation matrix from the reference frame to the body frame.
The derivative of this function with respect to time is where, S = (D 23 − D 32 , D 31 − D 13 , D 12 − D 21 ),D ij is the element corresponding to the rotation matrix D.

If the expression satisfies the following equation
The control torque is (5) www.nature.com/scientificreports/

Dynamic equation in vicinity of equilibrium
It is difficult to accurately calculate all the external interference torque M ext for the satellite, so the external torque in this paper includes gravity gradient torque, and torque caused by the change of tensor of inertia, which are respectively The rotation matrix D rotates in the order of 2-3-1.α 1 , α 2 and α 3 represent roll, pitch and yaw respectively.
Linearized in vicinity of equilibrium, the expression omits that the second order minima To sum up, the relative angular motion equation of the satellite omitting high-order small quantities is where

Control torques generated by magnetorquers
The magnetorquers are installed along the body frame of the satellite.Three magnetic torques in mutually perpendicular directions are generated by the action of the external magnetic field.However, the control accuracy of the magnetorquers can only reach the order of degrees, which is far from meeting the requirements of attitude control accuracy.It also needs to be combined with other actuators to meet the requirements of the mission 19,20 .
The control torque generated by the magnetorquers is where m is the magnetic dipole vector and B is the magnetic induction intensity of the Earth.According to the control torque equation, the control torque is always orthogonal to the direction of magnetic induction intensity.When the satellite is near the equator, the magnetic torque can only control the pitch and yaw angle.When the satellite is near the poles, the magnetic torque can only control the roll and pitch angle.This means that at any time, there is always a direction where the magnetorquers cannot generate control torque.The dipole moment is generated by the electromagnetic coil in the magnetorquers, which is composed of electromagnetic materials with high permeability.Here, we assume that the input current does not exceed ±110 mA and the maximum magnetic dipole moment does not exceed ±27.5 A m 221 .
The magnetic induction intensity of the Earth is B , which is expressed in the orbital frame 22,23 where µ B is the geomagnetic constant, i is the orbital inclination, u is the latitude argument and its expression is u = ω 0 t + u 0 , ω 0 is the angular velocity of the satellite.The control torque generated by the magnetorquer is where, e B = B √ (B,B) .

Control algorithm of LQR
For optimal control of linear quadratic regulator (LQR), the cost function and system state equation are as follows ( 6) www.nature.com/scientificreports/where, x is the state vector, T f is the final time, P is a positive definite symmetric constant matrix, Q and R are posi- tive definite symmetric time-varying matrices, and u is the control torque.A is the dynamics matrix, B ctrl is the control matrix, w is the model noise, B d is the noise coefficient, and B is the magnetic induction intensity.
The optimal control is The magnetic dipole vector matrix is m = u .The differential Riccati equation of P(t) is The boundary condition of the equation is P T f = P f .

Control distribution of actuator
For a system composed of two or more actuators, how to reasonably assign virtual expected instructions to each actuator to satisfy the stability requirements of spacecraft attitude is the problem of actuator control allocation.The common method is to incorporate the control allocation of the actuator into the design of the control law.
Considering the applicability of the project, this paper adopts a simple and practical string link allocation rule.This distribution rule assumes that different actuators provide control torque according to different priorities.
After the actuators with higher priorities reach saturation, the remaining execution instructions are completed by the next level of actuators.As shown in Fig. 2, the control command of virtual expectation can be converted into the control command of each actuator only after control allocation.The control allocation expression form is as follows , where, v(t) is the desired virtual control instruction, u(t) is the input instruction of the actuator, and B eff is the control efficiency matrix.The control allocation problem can be transformed into The actuator uses a magnetic torquer for control distribution first, that is, We know that it is dif- ficult to meet the control requirements only with a magnetic torquer, that is In order to control a certain attitude of the satellite, two heterogeneous actuators need to cooperate with each other to complete the control task, that is, the input torque u(t) of the actuator is not unique.Since the magnetic dipole vector is limited to 27.5 A m 2 and the maximum thrust of the thruster is 10 mN , the actual per- formance index of the actuator is limited by physical conditions, and the input variable u(t) satisfies the inequality u min ≤ u(t) ≤ u max .Different actuators have different corresponding rates of instructions.This characteristic is represented by the rate u(t) of output instruction u(t) , i.e. ρ min ≤ u(t) ≤ ρ min .Therefore, the restrictions of actuators are

Control torque generated by thruster
It is difficult to achieve the attitude control accuracy by using only the torque generated by the magnetorquers.Therefore, micro thrusters with a maximum thrust of 10 mN need to be installed on the satellite, in order to achieve the attitude control accuracy of the satellite.Because it will produce a large torque in a short time, the attitude and angular velocity of the satellite will change almost simultaneously when the micro thruster is working.When the relative attitude and relative angular velocity reach the allowable boundary, the actuator starts to work.It should be noted that the switch of thruster does not consider the delay phenomenon.The control torque M ctrl2 generated by the micro thruster is where In order to solve the approximate values of k α , k ω , using Floquet theory, the approximate initial condition of system (7) is �(0) = I 6×6 , and ρ k is the characteristic root of characteristic equation det (�(T) − ρ k I 6×6 ) = 0 .To make the linear system (α, ω) T = 0 asymptotically stable in a large range near the equilibrium position, the inequality kRe(ln ρ k ) < 0 for any ρ k is constant.When max [Re(ln ρ k )] is the smallest, the system returns to the equilibrium position faster.Therefore, the most suitable control parameters for system stability are obtained through max k [Re(ln ρ k )] → min.

Probability and momentum of collision
In this paper, the average orbital height of the formation is h = 450km , eccentricity is e = 0.001 , orbit inclina- tion is i = 89.5 • , and average effective cross-sectional area is S ≈ 3m 2 .Space debris was assessed using the international common software ORDEM2000.From the software ORDEM2000, the flux of space debris with a size greater than or equal to 0.1 mm over a 10-year mission period is = 1427.56393 m 2 10 year , the flux of debris with a size greater than or equal to 1 mm is = 3.0216 3 m 2 10 year , and the flux of debris impacted with a size greater than or equal to 10 mm is = 3.0734 × 10 −4 3 m 2 10 year。Therefore, the probability of the formation being impacted by space debris of centimetre size is relatively small, and that of millimetre or smaller (15) www.nature.com/scientificreports/size is relatively high, the impact of space debris of 0.1 mm size and above on the double satellite formation for gravity field exploration is investigated and analysed.Two independent collisions obey an exponential distribution, then the probability density function of time is 9 The number of space debris with size greater than or equal to 0.1 mm colliding the formation per unit time is = 4.5314 × 10 −6 s −1 .Therefore the probability density function of two independent collision times can be obtained from Eq. ( 20), as shown in Fig. 3.
The space debris impact the formation is a random event, obeying Poisson distribution, then the probability of double satellite formation being impacted N times is If the probability density function f X (x) is a continuous function, then the probability of random variable X in interval a and b is If the size of space debris (or its speed relative to the satellite) is within the interval [a, b] , the probability of collision with the formation in this interval satisfies Eq. ( 22), and the angle function of impact is where, f θ,φ x, y represents the probability density function of altitude and azimuth angle, the range of altitude angle is [−90 • , 90 • ] , and the range of azimuth angle is [−180 Before and after the double satellite formation is impacted by space debris, the momentum of the system composed of the debris and the formation is conserved.Because the actual motion of the satellite after being hit by debris is more complex-both translational and rotational motion.For the convenience of the study, only the translational motion of the satellite is studied when the impact passes through the satellite center of mass.Only the rotation of the satellite is studied when it does not pass through the center of mass.
When the extension of the debris velocity passes the satellite center of mass, the momentum lost by the debris is equal to the increased momentum of the satellite.However, when the extension of the debris velocity does not pass through the satellite center of mass, the momentum lost by the debris is converted into the increased angular momentum of the satellite.
Assuming that at time τ , the i-th space debris hits the satellite formation for a sustained impact time of ε , and that the change in linear momentum of the satellite during that time is P i and the change in angular momentum is H i , the increase in force or moment of the satellite upon impact is where δ(t) is the unit step function.The relationship between angular momentum and momentum is H = r × P. www.nature.com/scientificreports/When the satellite is hit by debris at high speed, the direction of debris momentum does not pass through the satellite center of mass, the moment M impact is obtained immediately after the satellite is hit, and the attitude and angular velocity of the satellite will change accordingly.When the direction of debris momentum passes through the satellite center of mass, the force F impact will be obtained immediately after the satellite is hit, and the position and velocity of the impacted satellite will change accordingly.
This paper mainly studies the attitude stability analysis of formation after being impacted by debris.The angular velocity of the satellite after impact is where, ω A is the angular velocity of the satellite after impact, ω B is the angular velocity of the satellite before impact, and dω is the increased angular velocity of the satellite after impact.
Since the angular momentum transferred to the satellite is instantaneous when the space debris collides with the satellite.When the momentum of the debris after the collision with the satellite is all converted into the increased angular momentum of the satellite, where, r is the vector from the satellite center of mass to the impact point, p is the momentum of the debris, and J is the rotational inertia of the satellite.
From Eqs. ( 25) and ( 26), the angular velocity of the spacecraft after impact is Assuming that the satellite formation is impacted by the i-th space debris whose the mass is m i , the velocity is v i , the height angle is θ i , and an azimuth angle is φ i , and the linear momentum of the space debris is p i in the orbital frame can be expressed as Space debris has various shapes and types.For the convenience of research, this paper assumes that the debris is a sphere with a density of ρ , and the mass of the debris is where, l i is the size of the debris, subject to software ORDEM2000.
The impact of space debris on satellite is an inelastic collision process, in which momentum is conserved but energy is lost.Due to the complexity of the high-speed impact of space debris on satellites, this paper assumes that the impact of space debris on satellites is an inelastic collision process, in which momentum is conserved but energy is greatly lost.
The change of momentum after the satellite is impacted by debris is In order to accurately calculate the angular momentum of the satellite after being impacted, it is necessary to know the vector r .However, it is complicated to accurately calculate its specific expression.In this paper, it is assumed that the impact position is subject to uniform distribution in the impact plane, then the height angle and azimuth angle of the i-th space debris impact in plane A j are θ ji and φ ji respectively.
The probability of the j ∈ {1, • • • , 6} plane of the satellite being impacted is where, A j is the area of the j-th plane, and ψ j is the angle between the momentum direction of the debris and the normal in the j-th plane before impact.

The expression of vector r in body frame is
The two coordinates in x ji , y ji , z ji are uniformly distributed, and the other coordinate can be easily deter- mined according to impact surface A j .
In the process of high-speed collision, the momentum reduced by space debris is transferred to the satellite.At the same time, the sputter generated by debris collision will be ejected against the direction of the formation operation, taking away a small part of momentum.The change in the momentum of the system occurs as follows (25) www.nature.com/scientificreports/where, ξ is the momentum enhancement factor, indicating the impact of the sputter on the transferred momen- tum after the satellite is impacted by space debris.

Motion control simulation without debris impact in scientific mode
In order to verify the attitude controller developed in this paper, the parameters in Tables 1 and 2  Figure 4 shows that the control algorithm designed by Lyapunov has higher control accuracy, and the control accuracy is within ± 1.5°, when the actuator is only the magnetorquer and all initial values are 1°, while the  control algorithm designed by LQR is accurate within ± 6°.Both control algorithms show that the pitch axes have the highest control accuracy, as that there is always a magnetic torque control effect on this axis.Figure 5a shows that the magnetic dipole vector transitions frequently between saturated states ± 27.5 A m 2 , while Fig. 5b shows that the magnetic dipole vector can initially reaches saturation, but later it is below ± 10 A m 2 , indicating that the Lyapunov-based control algorithm is more effective.The control results of the control algorithms based on Lyapunov and LQR when the thruster provides thrust torque are shown in Figs. 6, 7 and 8. Figure 6 shows the relative attitude control of the satellite achieved by the two control algorithms, and three axes can meet the attitude control accuracy.Figure 7 shows the relative angular velocity control of the satellite implemented by the two control algorithms, three axes also satisfy the angular velocity control requirement.Figure 8 shows the thrust torque required by the two control algorithms.Obviously, the torque required based on Lyapunov control algorithm is small, and the switching firing frequency is also low.Over a 24-h period, the total firing frequency of Lyapunov are 26 less than that of LQR, which is similar to the results of literature 17 .

Monte Carlo simulation and critical momentum statistics
There are three changes to the control system after the double satellite formation being impacted by space debris.Firstly, the accuracy of the control does not vary significantly because the control system designed is somewhat robust.Secondly, the control accuracy exceeds the maximum allowed value, causing the system to exit the science mode.Thirdly, it causes the control system to diverge and the task to fail.Scientific mode means that the spacecraft is in the normal operation stage.Exit from scientific mode means that the spacecraft can not work normally, but may return to normal operation after relevant control.
The logical relationship of the change in the control system of the double satellite formation after a high-speed impact of space debris is shown in Fig. 9.
The critical momentum referred to is the momentum that the impact of space debris happens to cause the satellite to exit the science mode.As soon as one degree of freedom exits the science mode during the collision,  www.nature.com/scientificreports/we consider that the momentum gained by this impact on the satellite exceeds the critical momentum.The relationship between the Monte Carlo simulation and the control system is shown in Fig. 10.The density of space debris is taken as ρ = 2.8 g cm 3 , assuming a collision duration of ε = 0.1 s , which occurs at τ = 18, 000 s .The area of each surface of the satellite is shown in Table 3.
The thrustless attitude control of the satellite is shown in Fig. 11 when the satellite is hit by debris with an increased angular momentum of L 1 = −1.8016× 10 −4 kg m 2 s , L 2 = 1.6135 × 10 −4 kg m 2 s and L 3 = 0.8253 × 10 −4 kg m 2 s respectively.Figure 11a shows that the Lyapunov-based control algorithm has a large attitude change within a short period of time after impact, but quick plateaus.While Fig. 11b shows that the LQR-based control algorithm takes a long time to plateau after impact, which further illustrates the control law of high control accuracy and low robustness.
The change in attitude and thrust moment of the satellite when it acquires angular momentum of L 1 = −1.7921× 10 −2 kg m 2 s , L 2 = 1.1543 × 10 −2 kg m 2 s and L 3 = 0.7861 × 10 −2 kg m 2 s after impact are shown in Figs. 12 and 13 respectively.Figure 12 shows that when space debris with the same momentum hits the satellite formation, the attitude will exceed the requirements under the different control.Figure 13 shows that the control algorithm based on Lyapunov control can return to normal thrust conditions within 25 min, while the control algorithm based on LQR takes 48 min for the thrusters to return to normal thrust conditions under the impact of space debris of same momentum.
The process of Monte Carlo simulation and statistics is shown in Fig. 10.In this paper, 100,000 data satisfying the size distribution, velocity distribution and angle distribution of space debris software ORDEM2000 are selected.These data and the double satellite formation at the average orbital altitude randomly collide, and the probability density of two collision times conforms to the rule of Fig. 3.The results are as shown in Fig. 14 when       www.nature.com/scientificreports/Lyapunov control is adopted, the number of impacts greater than or equal to the critical angular momentum is 31, of which 5 times cause the control system to diverge.When LQR control is adopted, the number of impacts greater than or equal to the critical angular momentum is 38, of which 7 times cause the control system to diverge.
Figure 15 shows that the probability of the double satellite formation exiting the science mode and divergence due to at least one impact is as shown in Fig. 15a. Figure 15b shows that the Lyapunov-based control algorithm has higher stability after being hit.In the 10-year mission period, the formation for gravity field exploration has been impacted by 1428 space debris with size greater than or equal to 0.1 mm.For the two control algorithms, debris impact can cause the system to exit the science mode and cause the control system to diverge.There is a slight difference in the probability of normal manner, but the difference is small when using different control algorithms.

Conclusion
In this paper, the stability of relative attitude of the double satellite formation for gravity field in space debris environment is studied.We established the dynamical equation of random collision of relative attitude, and adopted two control algorithms, Lyapunov control and LQR.Under the corresponding conditions, the space debris distribution function is established with the international space debris software, and the probability density function of the time of two independent collisions is obtained according to the corresponding conditions.The impact of debris on the attitude control system was simulated by Monte Carlo simulation with 100,000 data satisfying the model.The results show that during the 10-year mission period, using the control algorithm designed by Lyapunov control, 31 impacts caused the satellite to exit science mode, 5 impacts caused its control to diverge.Based on the LQR, 38 impacts caused the satellite to exit the science mode, and 7 impacts caused the control system to diverge.This shows that the probability of satellite being knocked over vary in a small range due to different control algorithms, the two algorithms that meet the attitude control accuracy will both exit the science mode and unstability in the space debris environment.The control system still has the risk of interrupting the scientific detection mode.It is necessary to consider the satellite operation and maintenance technology and further study the countermeasures.

•
O I − X I Y I Z I inertial frame-IF.The origin is located in the Earth center of mass, O I Z I is the rotation axis of the Earth and O I X I points to the vernal equinox of J2000 epoch.• o O − x O y O z O orbit frame-OF.The origin is located in the satellite center of mass, o O z O points to the Earth center, and o O x O is located in the orbit plane perpendicular to o O z O and pointing to the direction of motion.• o B − x B y B z B body frame-BF.Its axes are the principal axes of inertia for the satellite.oB x B is the sight line of the laser, and o B z B is perpendicular to the bottom of the satellite.• o R − x R y R z R reference frame-RF.The origin is the midpoint of the line from following satellite to the main satellite, where o R x R points from the follow satellite to the main satellite, and o R z R is perpendicular to o R x R in the orbit plane.

Figure 3 .
Figure 3. Probability density functions of time between independent collisions.

Figure 6 .
Figure 6.Control of satellite attitude with thruster action.

Figure 7 .
Figure 7.Control of angular velocity with thruster action.

Figure 10 .
Figure 10.Relationship between attitude control and Monte Carlo simulation.

Figure 11 .
Figure 11.Attitude without thruster control after space debris impact.

Figure 12 .
Figure 12.Satellite attitude after space debris impact.

Figure 14 . 8 Figure 15 .
Figure 14.The number of times the impact that causes the formation to exit the scientific mode and diverge under different control algorithms.

Table 2 .
Attitude control accuracy and maximum thrust torque.

Table 3 .
Area of each surface of the satellite.