The 
 
 
 J
 
 
 2
 
 
 Relative Perturbation Analysis of Satellite Formation under the Requirement of Relative Position Maintenance with Millimeter-Level Accuracy

With high-precision DEM (Digital Elevation Model) and GMTI (Ground Moving Target Indicator) as the demand background, the influence of 
 
 
 
 J
 
 
 2
 
 
 
 zonal harmonic term perturbation on the relative motion of the millimeter-level short-range leader-follower satellites in near-circular orbit is studied through the relative perturbation method. An equation of motion that can describe the motion of the leader-follower satellites under the influence of 
 
 
 
 J
 
 
 2
 
 
 
 perturbation in near-circular orbit is derived, and the characteristics of the trajectory of in-plane periodic motion are analyzed. A study shows that under the influence of the relative perturbation of the 
 
 
 
 J
 
 
 2
 
 
 
 term, the in-plane periodic motion of the leader-follower satellites in near-circular orbit is a symmetrical closed “drop-shaped” trajectory with a period of 
 
 2
 π
 /
 n
 
 . By comparing with the results of numerical simulations, the correctness of the conclusions obtained in this paper is verified. According to the research results, it can be known that only using a thruster as the actuator to maintain the relative position can no longer meet the requirements of the long-term mm-level relative position maintenance. In the future, a new technical approach needs to be explored to achieve the long-term relative position maintenance with millimeter-level control accuracy.


Introduction
The satellite formation system [1] refers to a space system composed of multiple physically disconnected satellites in order to complete the same mission. Satellite formation can not only complete complex tasks that cannot be completed by a single satellite but also have more advantages than a single satellite in terms of reliability, flexibility, and economy. Therefore, satellite formation systems have become a research hotspot in the aerospace field and are widely used in many other fields, such as earth observation, astronomical observation, and deep space exploration. This article takes the high-precision DEM (Digital Elevation Model) and GMTI (Ground Moving Target Indicator) tasks as the background. As shown in Figure 1, the leader-follower satellite means that the deputy satellite and the chief satellite operate in the same orbit in a certain sequence, and the orbit elements of leader-follower satellites are only different in argument of latitude. Therefore, the leader-follower satellite is equivalent to a sparse antenna array, using ground echo interferometric imaging principles and formation flying satellite control technology to complete DEM and GMTI tasks.
Whether the satellite formation system can meet the DEM mission measurement accuracy and the minimum detectable speed of the GMTI mission not only depends on the performance of payload, length of baseline, and data processing but also depends on the control accuracy of relative position maintenance of the payload of the leader-follower satellites. In low-precision observation tasks, it only needs to control the relative position accuracy of the interference baseline at the meter or decimeter level. Therefore, in the study of relative dynamics, only the secular and longperiodic effects of various perturbations on the relative position need to be paid attention to. As users have higher and higher requirements for performance and accuracy, the accuracy requirements for the relative position maintenance of formation flying satellite have also been increased from the meter level and decimeter level to the millimeter level, which raises the higher control accuracy of the relative position maintenance of the satellite formation. Therefore, it is necessary not only to offset the effects of the secular and longperiodic contributions due to perturbation, but also to compensate for the short-periodic impact of perturbation. Hence, it is important to conduct a detailed analysis of the shape of the trajectory of formation flying satellites under the impact of perturbation. Among the many perturbations, the influence of the J 2 term perturbation occupies a major position. So this paper's aim is mainly to conduct a detailed analysis of the influence of the J 2 perturbation on the relative position of the leader-follower satellites.
A lot of research efforts had been done in the field of satellite relative motion under the influence of J 2 perturbation based on state transition matrices (STMs) and relative orbital elements (ROEs) [2]. Schweighart and Sedwick [3] derived a linear model by incorporating the first-order term of the J 2 perturbation into the C-W equation with a meter-level position error when compared with a numerical integration. Izzo et al. [4] considered the atmospheric drag and J 2 perturbation in the extension of the C-W model and used the Floquet theory to predict the relative motion under the perturbations with a higher accuracy than that of the C-W equation. Stringer et al. [5] used the quadratic Volterra method developing a nonlinear model for including the J 2 perturbation with submeter-level accuracy.
Unlike the above research, Schaub and Alfriend [6] and Gim and Alfriend [7] considered that the ROEs vary slowly with time and developed the ROE-based equations for relative motion under J 2 perturbation. Schaub and Alfriend [6] used the J 2 -perturbed Hamiltonian deriving a relative state equation defined by Delaunay element differences and explored the impact of J 2 perturbation on relative orbits. Gim and Alfriend [7] developed a state transition matrix (STM) including the effect of secular, long, and shortperiodic contributions due to J 2 to the first order. Gim and Alfriend [8] then developed a model that can avoid the singularity on equatorial orbits by using differential equinoctial elements. D'Amico [9] developed a linear relative motion model which captures the time evolution of quasinonsingular ROE subjects (relative semimajor axis, the relative mean longitude, the relative eccentricity, and inclination vectors) to first-order J 2 perturbations in near-circular orbits. By using the same ROE state, Gaias et al. [10] developed a more complete relative motion model which corrects the original J 2 formulation shortcomings of D'Amico's [9] work and includes the effects of atmospheric drag, with a higher accuracy than that of D'Amico's [9] model. In addition, Mahajan et al. [11] developed a STM for relative motion under the complete zonal perturbations based on the studies of Gim and Alfriend [7]. More recent work by Russell et al. [12] developed an analytical STM for relative motion including J 2 and J 3 perturbations by using the Vinti theory. Burnett et al. [13] derived a new linearized differential equation model including zonal harmonic perturbations and corresponding perturbation solutions for relative motion without averaging the perturbing acceleration or its gradient. The perturbation solution of their model depends only on time, the initial chief orbital elements, and the deputy initial conditions in the chief Hill frame and has a similar error to the GA-STM for zero initial chief orbit eccentricity.
As a nonlinear extension to the ROE-based linear relative motion, Alfriend [14] and Alfriend and Yan [15] derived a nonlinear model for propagating the mean OED in arbitrarily eccentric orbits while including second-order J 2 effects. Based on the Gim-Alfriend's studies, Sengupta et al. [16] developed a second-order state transition tensor (STT) and a nonlinear solution including the effects of J 2 perturbation by combining the Keplerian STT with Gim-Alfriend's [7] linear model. Then, Yang et al. [17] derived a more complete second-order STT for relative motion under the J 2 -perturbed elliptic orbits based on the Sengupta et al.'s [16] work. The two studies above [16,17] did not consider the effect of J 2 perturbation of the second-order transformation from the rectilinear relative state to the osculating ROEs and the second-order propagation of the mean ROEs. Similar to the methods in [16,17], Zhen et al. [18] developed a second-order analytical STT for relative motion under the J 2 -perturbed elliptic orbits by using the geometric method, and the result is more accurate than that of previous linear or nonlinear analytic methods. Moreover, Mahajan et al. [19] considered the J 2 zonal harmonic and sectorial and tesseral harmonics in their work, presenting a STM for perturbed satellite relative motion with the position errors below 1% of the formation size after two days in both the lowly and highly eccentric reference orbits. Recently, Gaias et al. [20] presented an analytical framework for the precise modeling of the relative motion in low earth orbits. Their model includes the first-order expansion of the effects due to any even zonal harmonics and the second-order expansion of the unperturbed and J 2 terms and can achieve high-precision mean/osculating orbital element conversions.
In addition, in order to achieve high-precision relative orbit maintenance, some scholars studied the hovering formation and orbit determination on the basis of the above dynamic research in recent years. Rao et al. [21] obtained a "teardrop" hovering formation by designing a set of relative orbit elements and proposed a new impulsive control strategy to keep the deputy satellite in the hovering pattern for a long time. Sun et al. [22] proposed an innovative orbit  International Journal of Aerospace Engineering determination method which makes use of gravity gradients for low earth-orbiting satellites. Bai et al. [23] proposed a teardrop hovering formation design for the elliptical reference orbit with a J 2 perturbation. In response to the above-mentioned background and research status, this paper studies the shape of relative motion of the leader-follower satellites in near-circular orbit under the effect of J 2 perturbation. On the basis of the predecessors, an equation of motion that can describe the relative motion of the leader-follower satellites in near-circular orbit under the influence of J 2 perturbation is derived, and the characteristics of in-plane periodic motion of the leaderfollower satellites under the influence of J 2 term perturbation are analyzed. The equation of motion derived in this paper has a very simple form and is more desirable and achievable in onboard computers. By comparing with the result of numerical simulations under the same conditions, the relative position error of this equation is less than 20 mm in the radial direction and 60 mm in the tangential direction in a day. The major contributions of this study are summarized as follows: (1) We derive a high-precision equation of motion that can describe the relative motion of the leaderfollower satellites in near-circular orbit under the influence of J 2 perturbation. The equation has a very simple form and is more desirable and achievable in onboard computers, so it can be used as a state equation for high-precision relative navigation (2) We analyze the characteristics of relative motion of leader-follower satellites under the influence of J 2 perturbation, which can provide a theoretical basis for the design of controllers of mm-level relative orbit maintenance The outline of this paper is as follows: the coordinate system and the mathematical description of formation flying are introduced in Section 2. An equation of motion that can describe the relative motion of the leader-follower satellites in near-circular orbit under the influence of J 2 perturbation is derived, and the characteristics of in-plane periodic motion of the leader-follower satellites under the influence of J 2 term perturbation are analyzed in Section 3. Numerical simulations used to validate the results in this paper are shown in Section 4. In Section 5, we summarize the study and propose a new method for future relative orbit maintenance.

The Mathematical Description of Formation Flying
In order to derive the equation of the motion of leaderfollower satellites under the influence of the J 2 relative perturbation, it is first necessary to establish a coordinate system and a mathematical model to describe the relative motion [24] of the satellite formation. As shown in Figure 2, denote the chief satellite as S and the deputy satellite as C. The orbital coordinate system O S − xyz of the chief satellite S is defined as follows: the origin O S is located at the center of the chief satellite and moves with it, the x-axis takes the radial direction of the reference satellite and points from the center of the earth to S, the z-axis points to the positive normal direction of the chief satellite S orbital plane, and the y-axis is determined by the right-hand rule.
Assuming that the chief satellite S is moving in a nearcircular orbit, the following C-W equation [25] can be obtained to describe the relative motion while ignoring high-order small quantities (since the distance between the chief satellite and the deputy satellite ρ is much smaller than the geocentric distance of the chief satellite r, we can ignore the ðρ/rÞ 2 and higher power terms): where ðx, y, zÞ is the relative position of the deputy satellite C in the orbital coordinate system of the chief satellite S, n = ffiffiffiffiffiffiffiffi ffi μ/a 3 p is the average angular velocity of the chief satellite S, μ is the geocentric gravitational constant, and a is the semimajor axis of the chief satellite. Δf x , Δf y , and Δf z are the differences between the perturbation acting on the chief satellite and the deputy satellite, that is, the relative perturbation.
Taking the initial moment as the moment when the chief satellite is located at the RAAN (Right Ascension of Ascending Node), in a near-circular orbit, the argument of latitude of the chief satellite is Then, there is In the context of this paper, the effect of the perturbation  (1); then, the independent variable t of the C-W equation can be changed to u, and we have In order to facilitate the analysis of the motion deviation of the deputy satellite in the relative motion orbit under the influence of relative perturbation, let ðΔx,Δy,ΔzÞ replace ðx, y, zÞ in the above formula; then, equation (4) can be changed to where Δx = x − x p , Δy = y − y p , Δz = z − z p , ðx, y, zÞ is the position of the deputy satellite in the relative motion orbit only under the influence of two-body gravity, and ðx p , y p , z p Þ is the actual position of the deputy satellite in the relative motion orbit under the influence of relative perturbation and two-body gravity. The solution ðΔx,Δy,ΔzÞ of equation (5) can describe the relative motion only under the influence of relative perturbation.

J 2 Relative Perturbation Analysis
In this section, the relative perturbation method is used to derive the equation of relative motion of the leader-follower satellites under the influence of J 2 perturbation. The gravitational potential function of the zonal harmonic term of the earth [26] is where R e is the radius of the earth, ðr, φÞ is the geocentric distance and geocentric latitude of the spacecraft in the groundfixed coordinate system, and J 2 , J 3 , and J 4 are the coefficients of zonal harmonic term perturbation. Then, the perturbation potential function of J 2 (the central gravitational term is omitted) is The components of the perturbation force in spherical coordinates are From the relationship between satellite spherical coordinates and orbital coordinates, it can be seen that the radial perturbation forces under satellite spherical coordinates and orbital coordinates are the same, and the tangential and normal perturbation forces f y and f z have the following conversion relationship where β is the angle between the orbital plane and the meridian plane of the point where the satellite is located. From the spherical triangle relationship, we have From equations (6), (7), (8), (9), and (10), we have the perturbation force of the J 2 zonal harmonic term in orbital coordinates: where k = ð3/2ÞJ 2 ðμR 2 e /r 4 Þ and i and u are the inclination and In equation (11), by taking partial differentiation of r, i, and u, we can obtain the following equation: In the above formula, Δr is the difference of the geocentric distance between the deputy satellite and the chief satellite, Δi is the difference of the inclination between the deputy satellite and the chief satellite, and Δu is the difference of the argument of latitude between the deputy satellite and the chief satellite.
As we all know, the orbit elements of leader-follower satellites are only different in argument of latitude. Therefore, in a near-circular orbit, there are no differences in the geocentric distance and inclination between the chief and deputy satellite, that is, Δr = Δi = 0. Then, equation (12) can be changed to Δf x = 3k sin 2 i sin 2u•Δu, Assuming that the deputy satellite has no relative position and velocity deviation between the influence of relative perturbation and two-body gravity at the initial moment, the initial value is Substituting equations (14) and (15) into equation (5), we can get the following analytical solution: In formulas (16), (17), and (18), ðΔx,ΔyÞ is the in-plane motion trajectory of the leader-follower satellites in nearcircular orbit only under the influence of J 2 perturbation, and Δz reflects the out-of-plane motion only under the influence of J 2 perturbation. It can be seen from the above formula that the relative position deviation in the three directions under the influence of J 2 perturbation is positively correlated with the difference of argument of latitude Δu between the deputy satellite and the chief satellite. Therefore, the greater the distance between the leader-follower satellites, Table 1: Special point parameters in the process of the in-plane motion only under the influence of J 2 relative perturbation. Because the relative motion of short-range leaderfollower satellites in the z direction under the influence of J 2 perturbation is less than 1 mm, we mainly discuss inplane motion. For in-plane motion, let Then, formulas (16) and (17) can be simplified as The parametric equation represented by equation (20) describes the in-plane motion of leader-follower satellites in near-circular orbit only under the influence of J 2 perturbation. In equation (20), by finding the first-order and second-order derivatives for t, we can obtain the velocity and acceleration of the in-plane motion only under the influence of J 2 perturbation as follows: By analyzing equations (20), (21), and (22), it can be obtained that the in-plane motion of leader-follower satellites in near-circular orbit only under the influence of J 2 perturbation has the following law: (a) The curve described by the parametric equation (20) is a symmetrically closed "drop-shaped" line as shown in Figure 3, and Δx = 0 is the symmetry axis of the curve. Let t = 0 be the starting time of the inplane motion, which corresponds to Oð0, 0Þ in the figure. The deputy satellite moves periodically in the clockwise direction from point O, and the period is 2π/n. The detailed parameters at the special points in the motion are shown in Table 1 -5    (20), the parameter "A" determines the magnitude of the in-plane motion only under the influence of J 2 perturbation. When the distance between the leader-follower satellites is certain, the magnitude of "A" is also related to the orbital altitude and the inclination of the chief satellite. As the orbital altitude increases, "A" continues to decrease, and the in-plane relative position deviation of the leaderfollower satellites also decreases. When 0°< i < 90°, "A" increases with the increase in inclination i, and the in-plane relative position deviation of leaderfollower satellites also increases. When i = 90°, the in-plane relative position deviation reaches its maximum value. When 90°< i < 180°, "A" decreases as the inclination i increases, and the in-plane relative position deviation of the leader-follower satellites also decreases (c) Figures 4 and 5 show the velocity and acceleration of the deputy satellite's in-plane motion in near-circular orbit only under the influence of J 2 relative perturbation at an orbital altitude of 500 km. It can be seen from the figure that in order to control this in-plane motion within the mm level, the required acceleration is of the order of 10 −7 m/s 2 , and the thruster is required to provide a continuous thrust of about 0.1 mN. In order to eliminate this movement caused by the relative perturbation of J 2 , more frequent thruster control must be relied on.
The long-term frequent jets of the thruster will increase fuel consumption and greatly reduce the service life of the satellites. Therefore, at this stage, simply using the thruster as the actuator cannot achieve mm-level long-term relative position maintaining control

Simulation
According to the background, the simulation parameters select the commonly used orbital parameters for earth observation using synthetic aperture radar satellites. The two satellites are, respectively, selected for simulation in the form of forward and backward flying formation with a distance of 100 m and 200 m. Select a low-earth circular orbit with an orbit altitude of 500 km, and the remaining initial orbit parameters are shown in Tables 2 and 3. The numerical simulator takes the initial analytic equations of motion as described in equation (23) and integrates the absolute motion of each satellite in time. The motion of each satellite is then differenced and converted to a coordinate system that was mentioned in Section 2.
In order to eliminate the error of integration, we take the initial analytic equations of motion as described in equation (24) as a numerical simulator and integrate the absolute motion of each satellite in time once more and perform the same processing as above. By subtracting the above two results, the relative motion result of the satellites only under the J 2 perturbation can be obtained, which can be compared with the analytical solution result in this paper. By comparing the conclusions obtained in this paper with the results of the numerical simulation, the correctness of the results obtained in this paper is verified.
The simulation results shown in Figures 6(b) and 7(b) show that under the influence of the relative perturbation of J 2 , the in-plane periodic motion of the leader-follower satellites is a symmetrical closed "drop-shaped" trajectory, which is coincident with the results of numerical simulations (the results shown in Figures 6(a) and 7(a)), and the relative position error of this equation is less than 20 mm in the radial direction and 60 mm in the tangential direction in a day, as shown in Figure 8.

Conclusions
In this paper, the relative perturbation method is used to study the influence of the relative perturbation of J 2 zonal harmonic term perturbation on the relative motion of leader-follower satellites under the requirement of millimeter-level accuracy in near-circular orbit. Studies have shown that under the influence of the J 2 relative perturbation, the in-plane periodic motion of the leader-follower satellites in near-circular orbit is a symmetrical closed "dropshaped" trajectory with a period of 2π/n. In order to eliminate this motion caused by the J 2 relative perturbation, more frequent thruster control must be relied on. The long-term frequent jets of the thruster will increase fuel consumption and greatly reduce the service life of the satellite. Therefore, only using a thruster as the actuator to maintain the relative position can no longer meet the requirements of the longterm mm-level relative position maintenance. In the future, we can follow the principle of flywheel angular momentum exchange and use the principle of linear momentum exchange for relative position control. By installing a twodimensional sliding mechanism on the two satellites and controlling the sliding mechanism to move, the linear momentum exchange caused by the continuous motion of the sliding mechanism can change the relative position of the formation satellites. Since the relative motion under the influence of J 2 perturbation is a continuous closed periodic trajectory and the motion amplitude is limited, this method can well offset this "drop-shaped" motion. Moreover, this method does not need to carry additional fuel and can provide continuous control power only by relying on electric energy. It has advantages over the use of thrusters in terms of accuracy and fuel saving.

Data Availability
The data used to support the findings of this study are available from the author upon request.

Conflicts of Interest
The authors declare that they have no conflicts of interest.