Next Article in Journal
Performance Improvement of Human Centrifuge Systems through Multi-Objective Configurational Design Optimisation
Previous Article in Journal
Time-of-Flight Measurements in the Jet of a High-Current Vacuum Arc Thruster
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Initial Identification of Thrust and Orbit Elements for Continuous Thrust Spacecraft in Circular Orbit

1
Graduate School, Space Engineering University, Beijing 101400, China
2
Taiyuan Satellite Launch Centre, Xinzhou 034000, China
*
Authors to whom correspondence should be addressed.
Aerospace 2023, 10(12), 1012; https://doi.org/10.3390/aerospace10121012
Submission received: 15 September 2023 / Revised: 29 November 2023 / Accepted: 29 November 2023 / Published: 1 December 2023

Abstract

:
Continuous thrust spacecraft in circular orbits have had a great influence on the identification and cataloging of space targets. Gaussian-type orbital element variational equations are simplified and approximated. Ground-based radar observation datasets are transformed into orbit elements datasets. The initial thrust and orbit elements are obtained by optimally solving the spatial parameter error sum of squares minimization problem with the Levenberg–Marquardt method. The simulation analysis is carried out under the high-precision orbit model, and the solution error of tangential acceleration is around 5 × 10−7 m/s2, and that of normal acceleration is around 3 × 10−6 m/s2; the accuracy of the semi-major axis is 350 m, and the accuracy of inclination is 0.095°. The method is applicable to the preliminary identification of thrust and orbit elements for circular orbit continuous thrust spacecraft and can provide reliable initial values for the subsequent precision orbit determination of such spacecraft.

1. Introduction

The research and application of continuous thrust spacecraft in circular orbit have developed rapidly in recent years with the increasing demand for low-orbit giant constellations, LEO to GEO continuous thrust orbit transfer, deep space exploration, and other space missions [1,2,3]. Up to now, there are nearly 300 high-orbit satellites applying electric propulsion technology, and more than 4200 “Starlink” satellites applying commercialized electric propulsion technology, accounting for nearly half of the world’s total number of in-orbit spacecraft. SpaceX plans to deploy about 12,000 “Starlink” satellites in near-Earth space by 2024, and about 42,000 “Starlink” satellites by 2027 [4]; this type of spacecraft has the advantages of rapid deployment, high precision, and excellent controllability, which enables it to achieve more complex space missions. As a measure of spacecraft propulsion technology, the United States, Russia, Japan, Europe, and China have all mounted continuous thrust payloads on spacecraft and achieved good performance. Continuous thrust spacecraft are also being used in a wider range of applications, which in the future will include areas such as human spaceflight and interplanetary resource acquisition.
The traditional method of initial orbit determination is to solve the parameter estimation problem for conic curves under perturbation conditions, that is, the spacecraft’s orbit in space can be approximated as a conic curve [5]. However, for a continuous thrust spacecraft, its orbit in space cannot be described by a conic curve, the unknown thrust acceleration increases the dimension of the problem, and the traditional initial orbit determination algorithm [6] is no longer applicable. JianRong Chen et al. [7] used the simplex method to solve the orbit determination problem of a continuous small-thrust maneuvering satellite based on the start and end times of electric propulsion control. However, it uses GPS data and cannot solve the problem for non-cooperative targets. Gary M. Goff et al. [8] used extended and traceless augmented Kalman filters, combining multiple models to estimate orbital continuity maneuvers at different thrusts, which were shown to be more accurate in the high-thrust case and less accurate in the low-thrust case. XingYu Zhou [9] proposed an algorithm that merges LSTM (Long Short-Term Memory Neural Network) and Kalman filtering. The fitting polynomial is used to represent its continuous maneuver process, but its higher-order polynomial will bring a higher computational burden. Qinglin Yang [10] utilized high-precision orbital dynamics models to estimate the velocity and thrust magnitude of an orbiting spacecraft in real time through Markov models and volumetric Kalman filtering. Yong Liu [11] applies an unscented Kalman filtering method to calibrate the tangential thrust of a spacecraft in circular orbit using GNSS data. In addition, it is also possible to perform parameter identification of the continuous thrust using spacecraft attitude information [12,13,14,15,16], but it is not possible to determine the orbit elements from it.
For cooperative targets, thrust parameter identification and orbit elements determination can improve the precision and accuracy of spacecraft motion control, help the spacecraft maintain stability in space environments that contain multiple complex perturbation forces, and achieve complex motion control tasks. For non-cooperative targets, mastering their thrust parameters and orbit elements is crucial for tasks such as precision orbit determination, spacecraft cataloging, and collision warning [5]. For most of the current continuous thrust spacecraft, both the orbit climbing section and the orbit holding section are near-circular orbits, and it is of practical significance to carry out the initial identification of thrust and orbit elements of a circular orbit continuous thrust spacecraft, which can provide a reliable initial value for the next precision orbit fixing. This is also the research focus of this paper.
The significance of this study is to achieve orbit determination and thrust parameter identification of a continuous thrust spacecraft in circular orbit relying only on ground-based radar observations without any prior information. Section 2 derives the continuous thrust perturbation equation in a geocentric inertial coordinate system and derives an approximate solution of the perturbation equation based on the set of thrust and orbit elements. Section 3 proposes a space parameter conversion method that enables the interconversion of thrust and orbit elements from radar observations. Section 4 takes advantage of the coupling relationship between the orbit elements to solve the thrust and orbit parameters by the least squares method. Section 5 conducts experimental analyses to verify that the proposed model can achieve the initial identification of thrust and orbit elements of continuous thrust spacecraft in circular orbit. Section 6 provides a summary and outlook.

2. Modelling of Orbit Motion

2.1. Thrust and Orbit Elements

In the RCN (Radial Circumferential Normal) reference system, the transverse and tangential directions of a circular orbit are the same. As shown in Figure 1, the acceleration in the tangential, radial, and normal directions is noted as [ F ¯ T , F ¯ R , F ¯ N ] T . The tangential direction is the direction of the spacecraft velocity. The radial direction is the opposite direction of the “spacecraft to the center of the earth” in the orbit plane. The normal direction is perpendicular to the orbit plane and is in a right-handed system with the tangential and radial directions.
The continuous thrust acting on the spacecraft is not a conservative force, and the corresponding perturbation equations of motion can be established in the form of perturbation acceleration components. Compared with the general Lagrange perturbation equation, the Gaussian perturbation equation is more intuitive in reflecting the influence of the acceleration component on the change rule of each orbit element, which is helpful for the subsequent derivation of formulas. The Gaussian variational equations of the orbit elements [17] are given to describe the process of continuous thrust orbit maneuver in circular orbit.
{ a ˙ = 2 a F ¯ T v ξ ˙ = 2 F ¯ T cos u v F ¯ R sin u v η ˙ = 2 F ¯ T sin u v + F ¯ R cos u v i ˙ = F ¯ N cos u v Ω ˙ = F ¯ N sin u v sin i u ˙ = n + 2 F ¯ R v F ¯ N sin u v tan i
where ξ = e cos ω , η = e sin ω , u = f + ω , a is the semi-major axis, e is the eccentricity, i is the inclination, Ω is the RAAN, ω is the argument of perigee, M is the mean anomaly, f is the true anomaly, n is the mean angular velocity, and v is the spacecraft velocity. For the present space missions applying continuous thrust, the orbits of the mega-constellation satellites such as Starlink, Oneweb, etc., are mostly maintained in near-circular orbits and no thrust is applied in the radial direction, which means that a = r and F ¯ R = 0 . In circular orbits, where just one single angle element (argument of latitude) is sufficient to describe the position of the spacecraft in orbit, it is not necessary to consider both the perigee magnitude and the true perigee angle. Considering the position of the orbit overall in space, the orbital equations for continuous thrust are simplified.
{ r ˙ = 2 F ¯ T n i ˙ = cos u v F ¯ N Ω ˙ = sin u v sin i F ¯ N u ˙ = n F ¯ N sin u v tan i
Observation of the above equation shows the following:
(a)
The semi-major axis variation is related only to the tangential acceleration and not to the orbit plane normal acceleration.
(b)
The inclination variation is related to the normal acceleration and the spacecraft velocity. There is a fixed mathematical relationship between the spacecraft velocity and the half-length axis, so the change in inclination is related to both the normal acceleration in the expression and the normal acceleration (the term affecting the change in the semi-major axis).
Combining the above conclusions, as shown in Figure 2, solving for δ = ( r , i , Ω , u , F ¯ T , F ¯ N ) leads to the determination of the thrust and orbital parameters of continuous thrust spacecraft in circular orbit.

2.2. Approximate Solution of the Perturbation Equations

2.2.1. Continuous Thrust Motion Control Equations

Observing Equation (2), it can be found that if the normal thrust does not change direction relative to the spacecraft, the change in inclination and RAAN in a period will be almost zero due to the integration effect of the cosine term cos u and sine term sin u, so the application of normal thrust will be meaningless. To maximize the change in inclination or RAAN, the normal thrust needs to change direction at a phase-symmetric position (u = ±π/2 or u = 0, π).
For an in-orbit spacecraft (mass: 1000 kg, inclination: 60°, orbit height: 550 km) under the influence of the J2 term perturbation, the RAAN will progress in a circular orbit [18]. The RAAN progression is different for different semi-major axes, so controlling the semi-major axis can make RAAN separation [19]. Also, the effect of normal thrust is compared to RAAN. Let the normal thrust change sign at u = 0, π, so that the integral of sin u in a period is maximum, at which time the effect of the normal thrust is most obvious. Figure 3 represents the efficiency of separating the RAAN using different semi-major axes and applying normal thrust, respectively.
The separation efficiency of J2 term perturbation on RAAN is much larger than the change efficiency of normal thrust on RAAN, while J2 term perturbation has almost no effect on inclination. Therefore, in engineering practice, the RAAN is mostly changed by J2 term perturbation between different orbit heights, and the inclination is changed by normal thrust. This study focuses on the change in spacecraft orbit height and inclination by continuous thrust. In this case, the normal thrust changes direction at ±π/2 to obtain the maximum change in inclination, which is derived as follows. The orbital averaging method [20] is applied during one period.
π / 2 π / 2 F ¯ N cos u d u + π / 2 3 π / 2 F ¯ N cos u d u 2 π = 2 F ¯ N π
π / 2 π / 2 F ¯ N sin u d u + π / 2 3 π / 2 F ¯ N sin u d u 2 π = 0
The average effect of F ¯ N sin u over a period is 0 and the average effect of F ¯ N cos u over a period is 2 F ¯ N / π . Combined with n = v / r , the above equation simplifies to
{ d r d t = 2 r v F ¯ T d i d t = 2 π v F ¯ N d Ω d t = 0 d u d t = v 3 μ

2.2.2. Analytical Model for Semi-Major Axis

The vis-viva equation of circular orbits (r = a) can be simplified as
v 2 = μ ( 2 r 1 a ) = μ r
Differential of the above equation, with
d v d t = d v d r d r d t = 1 2 μ 1 2 r 3 2 d r d t = 1 2 μ 1 2 r 3 2 2 r v F ¯ T = F ¯ T
Then,
v = v 0 + d v d t t = v 0 F ¯ T t
Deform the radius of the orbit:
d r d t = 2 r v F ¯ T = 2 μ v 3 F ¯ T = 2 μ F ¯ T ( v 0 F ¯ T t ) 3
Integrate the above equation:
r = 2 μ F ¯ T ( v 0 F ¯ T t ) 3 d t = μ ( v 0 F ¯ T t ) 2 + C r
where C r is the constant of integration; when t = 0, r = r 0 = μ / ν 0 2 + C r , so
C r = 0
Then,
r = μ ( μ r 0 F ¯ T t ) 2

2.2.3. Analytical Model of Inclination

The inclination is derived:
d i d t = 2 F ¯ N π v = 2 F ¯ N π 1 v 0 F ¯ T t
v0 is the initial velocity of the spacecraft. Integrate the above equation:
i = 2 F ¯ N π 1 v 0 F ¯ T t d t
i = 2 F ¯ N π F ¯ T ln ( v 0 F ¯ T t ) + C i
where C i is the constant of integration; when t = 0, i = i 0 = 2 F ¯ N π F ¯ T l n v 0 + C i , so
C i = i 0 + 2 F ¯ N π F ¯ T ln v 0
Then,
i = i 0 + 2 F ¯ N π F ¯ T ln ( μ μ F ¯ T t r 0 )

2.2.4. Analytical Model of Argument of Latitude

Derive the argument of latitude:
u ˙ = n = v r = v 3 μ
Then,
d u d t = 1 μ ( v 0 F ¯ T t ) 3
Then,
u = 1 μ ( v 0 F ¯ T t ) 3 d t
u = ( v 0 F ¯ T t ) 4 4 μ F ¯ T + C u
where C u is the Constant of integration; when t = 0, u = u 0 = ν 0 4 4 μ F ¯ T + C u , so
C u = u 0 + v 0 4 4 μ F ¯ T
Then,
u = u 0 + μ 2 r 0 2 ( μ r 0 F ¯ T t ) 4 4 μ F ¯ T

2.2.5. Analysis of Perturbation Effects

In practical engineering tasks, the influence of perturbation in near-Earth space ca not be ignored. Among them, the J2 term perturbation is the main influence term; it is of great practical significance to analyze the orbital motion law under its influence. The J2 term has almost no influence on the semi-major axis, inclination, and eccentricity. The effect on RAAN is as follows [21].
Ω ˙ J 2 = 3 J 2 a E 2 cos i 8 v r 3 = 3 J 2 a E 2 cos i 8 μ 3 ( v 0 F ¯ T t ) 7
Then,
Ω J 2 = Ω 0 + 3 J 2 a E 2 cos i 64 F ¯ T μ 3 ( v 0 F ¯ T t ) 8
The influence of the J2 term on the argument of latitude is as follows:
u ˙ J 2 = 3 J 2 a E 2 ( 3 4 sin 2 i ) 8 v r 3 = 3 J 2 a E 2 ( 3 4 sin 2 i ) 8 μ 3 ( v 0 F ¯ T t ) 7
Then,
u J 2 = u 0 + 3 J 2 a E 2 ( 3 4 sin 2 i ) 64 F ¯ T μ 3 ( v 0 F ¯ T t ) 8
Although the atmosphere in space is already extremely thin, spacecraft travel at speeds of several kilometers per second and atmospheric drag has a definite effect on spacecraft. For a spacecraft, the expression for the acceleration of atmospheric drag is
f d r a g = 1 2 C D S d m ρ v v
where CD is the drag coefficient, Sd is the windward surface area, m is the spacecraft mass, S/m is the surface-to-mass ratio, ρ is the atmospheric density, v is the spacecraft velocity, and v is the spacecraft velocity vector. Recent atmospheric density data published by NASA were used [22]. The values of the parameters involved in Equation (28) are taken: ρ = 1.38 × 10−16 g/cm3 at an altitude of 600 km, CD = 2.2, S/m = 0.02, and v = 7.56 × 103 m/s. The atmospheric drag acceleration is calculated to be 1.736 × 10−7 m/s2 in the opposite direction of the spacecraft velocity. From the performance of existing space electric thrusters [23], the acceleration produced by the thruster is generally in the range of 10−3–10−5 m/s2, which will be 2 to 4 orders of magnitude larger than the acceleration produced by atmospheric drag. The acceleration generated by atmospheric drag in the opposite direction of the velocity can be considered when identifying the thrust parameters of a continuously thrusting spacecraft to compensate for the desired tangential thrust acceleration.
Due to the presence of momentum in photons, sunlight striking the spacecraft generates a radiation pressure called solar pressure. The expression for solar pressure is
f S R P = C R S R m p s v s i s
where CR is the solar pressure coefficient, SR is the equivalent area, m is the spacecraft mass, ps is the solar pressure constant at the location of the spacecraft, and vs is the sun exposure coefficient, which is 1 when the spacecraft is illuminated by sunlight, and conversely 0. is characterizes the direction of the solar pressure perturbation force and is the unit vector of the spacecraft to the center of the sun. The literature [24] shows that the acceleration produced by solar pressure on a spacecraft is on the order of 10−7 and hardly affects the thrust acceleration produced by the electric thrusters.
The Earth’s non-spherical perturbation, atmospheric drag, and solar pressure are the perturbations that have the greatest impact on near-Earth spacecraft. In addition to these, there are several other smaller perturbations, such as third-body gravity, Earth deformation perturbations (oceanic, solid, and atmospheric tides), general relativistic effect perturbations, and Earth antiradiation pressure, which can usually be ignored in orbital calculations for general requirements.
In summary, knowing the initial orbital state δ0 = (r0, i0, Ω0, u0), the state of the spacecraft at time t can be expressed by the following equation.
{ r = μ ( μ r 0 F ¯ T t ) 2 i = i 0 + 2 F ¯ N π F ¯ T ln ( μ μ F ¯ T t r 0 ) Ω = Ω 0 + 3 J 2 a E 2 cos i 64 F ¯ T μ 3 ( μ r 0 F ¯ T t ) 8 u = u 0 + μ 2 r 0 2 ( μ r 0 F ¯ T t ) 4 4 F ¯ T μ + 3 J 2 a E 2 ( 3 4 sin 2 i ) ( μ r 0 F ¯ T t ) 8 64 F ¯ T μ 3
The high-precision orbital propagation model in Orekit is selected and used to analyze the created model in comparison with it, as shown in Figure 4.
The newly established orbital model has a small error with the high-precision orbit propagation model, and after one day of continuous thrusting, the error in the semi-major axis of the orbit between the two is 0.0039 m, the error in the inclination is 0.0006°, the error in the RAAN is 0.007°, and the error in the argument of latitude is 0.0031°. It is found by comparison that the established model approximates the error of the model considering the J2 term perturbation and can satisfy the error requirement of the initial orbit determination. In particular, the variation of the inclination is strongly influenced by the perturbation short period term, as shown in Figure 5.
Multi-arc data were selected for orbit determination and thrust parameter identification, and the perturbation of short-period terms had little effect on the results.

3. Transformation of Orbit Elements

3.1. Semi-Major Axis Transformation

From Equation (2), it can be seen that in the continuous maneuver process, the initial semi-major axis r0 and tangential acceleration F ¯ T can affect each parameter. The first equation in Equation (2) implies that the initial semi-major axis and tangential acceleration can be solved by studying the change law of the semi-major axis, which can simplify the process of solving the other parameters in the following. Under the J2000.2 equatorial reference system (geocentric inertial coordinate system), the radar observation data are transformed into spacecraft semi-major axis data. The spatial relationship between the geocenter, the station, and the spacecraft is shown in Figure 6.
r shows the position of the spacecraft, R shows the position of the station, and ρ shows the position of the spacecraft relative to the station.
r = R + ρ
In selecting the reference plane and its principal direction, the effects of exposure to differential age chapter movements and geopotential shifts are considered, as described in reference [17], and will not be repeated here.

3.2. Inclination Transformation

As is shown in Figure 7, a space Cartesian coordinate system in space is established with the Earth as the prime (O) and the equatorial plane (E) in the xOy plane, which can be expressed as z = 0.
In near-Earth space, when a spacecraft moves along its orbit from south to north, the intersection of the orbit with the equatorial plane is called the ascending node, and when it moves from north to south, the intersection of the orbit with the equatorial plane is called the descending node. The intersection line between the orbit plane and the equatorial plane (the line between the ascending and descending nodes) is called the “node line”. From the third equation of Equation (5), the RAAN is not affected by thrust during the continuous long-term maneuvers of the spacecraft. However, from Equation (24), under the influence of the J2 term perturbation, the node line will move, which means that there is an angular difference between the intersection line l0 and li formed by the intersection of the initial orbit plane L0 and the plane Li with the equatorial plane E, respectively, and the angular difference is ΔΩi. Knowing that the initial RAAN is Ω0, the node line l0 can be expressed as
l 0 : { y = tan Ω 0 x z = 0
The RAAN at time ti is Ωi = Ω0 + ΔΩi, and the intersection line li can be expressed as
l i : { y = tan ( Ω 0 + Δ Ω i ) x z = 0
P0 is the position of the continuous thrust spacecraft at the initial moment t0, and Pi denotes the position of the spacecraft at the moment ti; L0 is the initial orbit plane, Li is the orbit plane at the moment ti, and the expressions for P0:(x0, y0, z0) and Pi:(xi, yi, zi) are as follows:
{ P 0 : ( x 0 , y 0 , z 0 ) P i : ( x i , y i , z i )
The intersection line l0 and the straight line OP0 are both on the initial orbit plane L0 and are not parallel when u ≠ 0. The normal vector h0 to the plane L0 is
h 0 = l 0 × O P 0
h i = l i × O P i
e is defined as the normal vector of the equatorial plane E. E is the xOy plane in the space Cartesian coordinate system, so e = (0,0,1). The angle made by the plane L0 and the equatorial plane E is the initial orbit inclination i0.
i 0 = arccos ( h 0 e | h 0 | | e | )
As the above derivation, the inclination ii at time ti is
i i = arccos ( h i e | h i | | e | )

3.3. Argument of Latitude Transformation

The argument of latitude represents the angular distance between the spacecraft and the ascending node, and the angle between the straight line OPi and the node line li is the argument of latitude u.
u = arccos ( O P i l i | O P i | | l i | )

4. Thrust and Orbit Elements Solving Methods

Ground-based radar observation equipment acquired m sets of noise-containing measurements in distance, azimuth, and elevation:
y ˜ i = ( t i , ρ ˜ i , A z ˜ i , E l ˜ i ) , i = 1 , 2 , 3 , m
Define a sequence of observations of the following form:
Y ˜ = [ ρ ˜ 1 , A z ˜ 1 , E l ˜ 1 , ρ ˜ 2 , A z ˜ 2 , E l ˜ 2 , , ρ ˜ m , A z ˜ m , E l ˜ m ]
Let the true value of the observed sequence (without noise, unknown) be
Y ¯ = [ ρ ¯ 1 , A z ¯ 1 , E l ¯ 1 , ρ ¯ 2 , A z ¯ 2 , E l ¯ 2 , , ρ ¯ m , A z ¯ m , E l ¯ m ]
Assuming that the observation errors follow a normal distribution and that the detector line-of-sight errors have a mean squared error of σ2, the statistical properties of the observations are as follows:
{ ρ ˜ i ~ N ( ρ ¯ i , σ 1 2 ) A z ˜ i ~ N ( A z ¯ i , σ 2 2 ) E l ˜ i ~ N ( E l ¯ i , σ 3 2 )
Define the sequence of observations with respect to the semi-major axis, inclination, and argument of latitude according to Equation (41):
{ Y ˜ r = [ r ˜ 1 , r ˜ 2 , , r ˜ m ] Y ˜ i = [ i ˜ 1 , i ˜ 2 , , i ˜ m ] Y ˜ u = [ u ˜ 1 , u ˜ 2 , , u ˜ m ]
The initial semi-major axis and tangential thrust acceleration are denoted by x 1 0 = [ r 0 , F ¯ T ] . The initial inclination, RAAN, and normal thrust acceleration are denoted by x 2 0 = [ i 0 , Ω 0 , F ¯ N ] . The initial argument of latitude is denoted by x 3 0 = u 0 . Based on the analytical model, a set of theoretical orbital sequences can be obtained:
{ Y ( x 1 0 ) = [ r 1 ( x 1 0 ) , r 2 ( x 1 0 ) , , r m ( x 1 0 ) ] Y ( x 2 0 ) = [ i 1 ( x 2 0 ) , i 2 ( x 2 0 ) , , i m ( x 2 0 ) ] Y ( x 3 0 ) = [ u 1 ( x 3 0 ) , u 2 ( x 3 0 ) , , u m ( x 3 0 ) ]
The sum of the squares of the errors of the theoretical and actual observations characterizes the difference between the theoretical and actual orbits:
{ i = 1 m [ r i ( x 1 0 ) r i ˜ ] 2 j = 1 m [ i j ( x 2 0 ) i j ˜ ] 2 j = 1 m [ u j ( x 3 0 ) u j ˜ ] 2
The orbit determination problem can be described as the following optimization problem:
{ min x 1 0 J ( x 1 0 ) = i = 1 m [ r i ( x 1 0 ) r i ˜ ] 2 min x 2 0 J ( x 2 0 ) = j = 1 m [ i j ( x 2 0 ) i j ˜ ] 2 min x 3 0 J ( x 3 0 ) = j = 1 m [ u j ( x 3 0 ) u j ˜ ] 2
In solving the above optimization problem, the initial semi-major axis r0, tangential thrust acceleration F ¯ T , initial inclination i0, initial RAAN Ω0, normal thrust acceleration F ¯ N , and initial argument of latitude u0 can be obtained sequentially.
In summary, the problem of identifying the thrust and orbit elements of a continuous thrust spacecraft is transformed into a least squares problem, which can be solved for δ 0 = ( r 0 , i 0 , Ω 0 , u 0 , F ¯ T , F ¯ N ) by the above process.

5. Simulation Analysis

In order to identify the thrust and orbit elements of a continuous-thrust spacecraft in circular orbit, a full-model simulation in a real scenario is required, so as to verify the accuracy of the parameter-solving method in this paper. The thrust parameters are set in reference to the current Starlink V2.0mini satellite launched by SpaceX, whose mass is 750 kg, and the thrust of the equipped Hall electric thruster is 170 mN, producing a combined thrust acceleration of 2.27 × 10−4 m/s2. At this time, the yaw angle is set to be 30°, the tangential acceleration is 1.966 × 10−4 m/s2, and the normal acceleration is 1.135 × 10−4 m/s2. All simulation conditions are set as shown in Table 1.
Under the effect of various perturbation forces, the continuously thrusting spacecraft conducts continuous maneuvers along its orbit, and, in consideration of the actual radar observation, the radar observable range is limited to 5° after the spacecraft exits the station horizon and 5° before the inbound station horizon. Figure 8a shows the geometric schematic of ground-based radar visibility, with the orange area as the observable area. Figure 8b shows the ground-based radar visibility time window schematic (the grey area is the observable window period); the low-orbit satellite orbit cycle is short and can be observed in more windows, but the observation time of each window is short, which is generally about 10 min. In the simulation, the sampling step is set to be 60 s, and each sampling time is determined according to the length of the observation window, and all the observable windows in the time from 2020.04.02 00:00:00 UTC to 2020.04.04 00:00:00 UTC are sampled for the observation, with a time span of 48 h and a total of 154 groups of observational data, which are recorded as i = 1, 2, …, 154.
Equation (47) is a typical nonlinear least squares problem. In this paper, the Levenberg–Marquardt algorithm (LM algorithm) is used for least squares parameter estimation. This method combines the most rapid descent method and the Gauss–Newton method. At the beginning of the iteration, it is approximately equivalent to the most rapid descent method, which can gradually approach the optimal solution in the case of bad initial values; at the end of the iteration, it is approximately equivalent to the Gauss–Newton method, which can quickly converge to the optimal solution on the basis of better estimates. Applying the LM algorithm to this problem, the optimal values of orbit elements and thrust parameters can be found quickly.
The proposed model does not require high accuracy of the initial values in problem solving. For the semi-major axis, inclination, RAAN inclination, and argument of latitude, the sampled data at the initial moment are taken as the initial value. The thrust acceleration generated by the Hall electric thrusters carried by the existing on-orbit spacecraft does not exceed the order of 10−4 [20], so the tangential acceleration and normal acceleration can be taken as 0, which is the initial value when simulation analysis is performed. The solution results are shown in Table 2.
Using the above method to calculate the thrust and orbit elements, a total of 1000 Monte Carlo simulations were carried out; the relative error of the spacecraft tangential acceleration solution is about 5 × 10−7 m/s2, and the relative error of the normal acceleration solution is about 3 × 10−6 m/s2. The error of the semi-major axis is about 0.35 km, the error of solving the inclination is about 0.095°, and the error of the RAAN and argument of latitude are about 0.74° and 0.83°, respectively. Figure 9 shows the solution error of tangential acceleration, normal acceleration, semi-major axis, inclination, RAAN, and argument of latitude.
From the above results, the estimation accuracy of spacecraft thrust acceleration parameters is high and can be applied to engineering practice. Due to the observation error and system error, the estimated orbital elements have some errors. The estimated orbital elements can be used as reliable initial values for precise orbit determination of continuous thrust spacecraft. Next, the simulation is set to explore the thrust parameter identification ability of the proposed method.
At present, Hall electric thrusters are the thrusters most commonly carried by continuous-thrust spacecraft, and the thrust that Hall electric thrusters can provide to satellites is in the range of 5 mN–1 N. The thrust was set to 1 N, 500 mN, 100 mN, 20 mN, 5 mN, and 1 mN for simulation, the satellite mass was set to 500 kg, and the thrust yaw angle was set to 60° (refer to Table 2 for other parameter settings). The simulation was carried out under different thrust forces, and the results are shown in Table 3.
The comparison between the simulation values and solutions values is shown in Figure 10.
The solution errors for different thrusts are shown in Figure 11.
From the simulation results, the tangential acceleration and normal acceleration have high identification accuracy when the thrust is greater than 5 mN. When the thrust is about 1 mN, the proposed method is able to identify the magnitude of the tangential acceleration, while the normal acceleration will not be identified. This is because the change in the orbit elements caused by the thrust is too small, resulting in the error in the identification of the thrust parameters within a short observation time.

6. Conclusions

An analytical method using geometrical parameter transformations is proposed, which can perform the solution of thrust and orbit elements of a continuously thrusting spacecraft during climb deployment, orbit maintenance, continuous orbital maneuvers, and re-entry segments. It has the following features:
(1)
The adoption of the analytical orbit model can reflect the relationship between the elements more clearly, and, at the same time, simplify the calculation and quickly solve the problem.
(2)
The decoupling of tangential thrust and normal thrust parameters can simplify the solving process, and the solving accuracy of tangential acceleration and normal acceleration can reach the order of 10−6 and 10−5 m/s2, which can be directly applied to engineering practice.
(3)
In the process of solving the RAAN and argument of latitude, only the long-term effect of the spacecraft by the J2 term perturbation is considered, and its short-term effect with the thrust acceleration is not considered, so the solution accuracy is lower, but it can be used as the initial value to be substituted into the next step of the precise orbit determination.
This study aims to achieve the fast pre-identification of thrust parameters and orbit elements, with the purpose of solving the rough position of the spacecraft in near-Earth space from the short-arc data of the multi-segment single-station radar observation, which can provide the initial value for the next step of the precision orbit determination for the continuous-thrust spacecraft. The near-circular orbit is regarded as a circular orbit, and the eccentricity and perigee angle are neglected in the process, which causes an acceptable error in the solution of the orbit elements. At the same time, there are some constraints in the research work, as the propulsion strategy with fixed acceleration in the near-circular orbit can only satisfy most of the continuous-thrust spacecraft operating conditions and cannot cope with the time-varying thrust conditions. The final parameter solution model does not consider the effects of atmospheric drag and solar pressure. The next step of the research focuses on the use of filtering methods to describe the variation of thrust parameters of spacecraft in real time, as well as the precision orbit determination of continuous-thrust spacecraft.

Author Contributions

Conceptualization, Z.L. and S.Z.; methodology, S.Z.; software, X.T.; validation, S.Z. and X.T.; writing—original draft preparation, S.Z.; writing—review and editing, X.T. and Z.L. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

Data are contained within the article.

Acknowledgments

The authors thank Yasheng Zhang for the helpful guide.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Sung, T.; Ahn, J. Optimal Deployment of Satellite Mega-Constellation. Acta Astronaut. 2023, 202, 653–669. [Google Scholar] [CrossRef]
  2. Anderson, J.F.; Cardin, M.-A.; Grogan, P.T. Design and Analysis of Flexible Multi-Layer Staged Deployment for Satellite Mega-Constellations under Demand Uncertainty. Acta Astronaut. 2022, 198, 179–193. [Google Scholar] [CrossRef]
  3. Fan, Z.; Huo, M.; Qi, N.; Xu, Y.; Song, Z. Fast Preliminary Design of Low-Thrust Trajectories for Multi-Asteroid Exploration. Aerosp. Sci. Technol. 2019, 93, 105295. [Google Scholar] [CrossRef]
  4. Ren, Y. The Development Status of Starlink and Its Countermeasures. Mod. Def. Technol. 2022, 50, 11–17. [Google Scholar] [CrossRef]
  5. Yu, S.; Wang, X.; Zhu, T. Maneuver Detection Methods for Space Objects Based on Dynamical Model—ScienceDirect. Adv. Space Res. 2021, 68, 71–84. [Google Scholar] [CrossRef]
  6. Tao, X.; Li, Z.; Gong, Q.; Zhang, Y.; Jiang, P. Uncertainty Analysis of the Short-Arc Initial Orbit Determination. IEEE Access 2020, 8, 38045–38059. [Google Scholar] [CrossRef]
  7. Chen, J.; Li, J.; Wang, X.; Zhu, J.; Wang, D. A Simplex Method for the Orbit Determination of Maneuvering Satellites. Sci. China Phys. Mech. Astron. 2018, 61, 024511. [Google Scholar] [CrossRef]
  8. Goff, G.M.; Black, J.T.; Beck, J.A. Orbit Estimation of a Continuously Thrusting Spacecraft Using Variable Dimension Filters. J. Guid. Control. Dyn. 2015, 38, 2407–2420. [Google Scholar] [CrossRef]
  9. Zhou, X.; Qin, T.; Meng, L. Maneuvering Spacecraft Orbit Determination Using Polynomial Representation. Aerospace 2022, 9, 257. [Google Scholar] [CrossRef]
  10. Yang, Q.; Zhou, W.; Chang, H. Real-Time On-Orbit Estimation Method for Microthruster Thrust Based on High-Precision Orbit Determination. Int. J. Aerosp. Eng. 2021, 2021, 7733495. [Google Scholar] [CrossRef]
  11. Liu, Y.; Song, Z.; Li, L. Tangential Low Thrust on Orbit Calibration Method for Circular Orbit Based on GNSS Precision Orbit Determination. Patent CN103940431A, 23 July 2014. [Google Scholar]
  12. Yu, J.; Lv, J.; Liu, W.; Tie, L.; Yang, L.; Lu, G.; Liang, Y.; Huang, J. On-Orbit High-Precision Calibration Method of Thruster Thrust Based on Angular Displacement Sensor. Patent CN109018433A, 18 December 2018. [Google Scholar]
  13. Zhang, H.; He, L.; Xuemai, G. On-orbit electric thruster calibration based on MME/KF algorithm. Spacecr. Environ. Eng. 2011, 4, 337–343. [Google Scholar]
  14. Wu, S.-F.; Steyn, W.H.; Bordany, R.E. In-Orbit Thruster Calibration Techniques and Experiment Results with UoSAT-12. Control. Eng. Pract. 2004, 12, 87–98. [Google Scholar] [CrossRef]
  15. Li, H.; Sun, Z.; Chen, X. On-Orbit Calibration Based on Second-Order Nonlinear Filtering. In Proceedings of the 2008 International Conference on Machine Learning and Cybernetics, Kunming, China, 12–15 July 2008; IEEE: Piscataway, NJ, USA, 2008; pp. 2970–2975. [Google Scholar]
  16. Wiktor, P.J. On-Orbit Thruster Calibration. J. Guid. Control. Dyn. 1996, 19, 934–940. [Google Scholar] [CrossRef]
  17. Vallado, D.A. Fundamentals of Astrodynamics and Applications, 5th ed.; Space Technology Library; Springer: Berlin/Heidelberg, Germany, 2022. [Google Scholar]
  18. Pirozhenko, A.V.; Maslova, A.I.; Vasilyev, V.V. About the influence of second zonal harmonic on the motion of satellite in almost circular orbits. Spacecr. Dyn. Control. 2019, 25, 3–11. [Google Scholar] [CrossRef]
  19. Bakhtiari, M.; Abbasali, E.; Daneshjoo, K. Minimum Cost Perturbed Multi-impulsive Maneuver Methodology to Accomplish an Optimal Deployment Scheduling for a Satellite Constellation. J. Astronaut. Sci. 2023, 70, 18. [Google Scholar] [CrossRef]
  20. Prussing, J.E. Primer Vector Theory and Applications. In Spacecraft Trajectory Optimization; Conway, B.A., Ed.; Cambridge University Press: Cambridge, UK, 2010; pp. 16–36. ISBN 978-0-521-51850-5. [Google Scholar]
  21. Abbasali, E.; Kosari, A.; Bakhtiari, M. Effects of oblateness of the primaries on natural periodic orbit-attitude behaviour of satellites in three body problem. Adv. Space Res. 2021, 68, 4379–4397. [Google Scholar] [CrossRef]
  22. NRLMSIS Atmosphere Model. Available online: https://kauai.ccmc.gsfc.nasa.gov/instantrun/nrlmsis/ (accessed on 22 August 2023).
  23. Guangqing, X.; Chang, L. Review and Prospect of Electric Propulsion Acceleration Technology. J. Astronaut. 2022, 43, 143–157. [Google Scholar]
  24. Bronshtein, I.N.; Semendyayev, K.A.; Musiol, G.; Mühlig, H. Handbook of Mathematics, 6th ed.; Springer: Berlin/Heidelberg, Germany, 2015. [Google Scholar]
Figure 1. Continuous thrust on spacecraft.
Figure 1. Continuous thrust on spacecraft.
Aerospace 10 01012 g001
Figure 2. Thrust and orbit elements.
Figure 2. Thrust and orbit elements.
Aerospace 10 01012 g002
Figure 3. Separation efficiency of RAAN in two ways.
Figure 3. Separation efficiency of RAAN in two ways.
Aerospace 10 01012 g003
Figure 4. Comparison between the new model and high precision orbit propagation model.
Figure 4. Comparison between the new model and high precision orbit propagation model.
Aerospace 10 01012 g004
Figure 5. Variation of orbit inclination.
Figure 5. Variation of orbit inclination.
Aerospace 10 01012 g005
Figure 6. The spatial relationship between the geocenter, station, and spacecraft.
Figure 6. The spatial relationship between the geocenter, station, and spacecraft.
Aerospace 10 01012 g006
Figure 7. Space motion of continuous thrust spacecraft.
Figure 7. Space motion of continuous thrust spacecraft.
Aerospace 10 01012 g007
Figure 8. Radar visibility analysis. (a) Visible range geometry; (b) visible window.
Figure 8. Radar visibility analysis. (a) Visible range geometry; (b) visible window.
Aerospace 10 01012 g008
Figure 9. Solution error.
Figure 9. Solution error.
Aerospace 10 01012 g009
Figure 10. Comparison of simulation values and solution values.
Figure 10. Comparison of simulation values and solution values.
Aerospace 10 01012 g010
Figure 11. The solution errors for different thrusts.
Figure 11. The solution errors for different thrusts.
Aerospace 10 01012 g011
Table 1. Simulation parameter settings.
Table 1. Simulation parameter settings.
ProjectParameter
Initial epoch time2023.04.02 04:46:39 UTC
Initial orbit [ r , i , Ω , u ] = [ 6933.5345   k m , 65.011 ° , 14.372 ° , 28.728 ° ]
Initial acceleration [ F ˉ T , F ˉ N ] = [ 0.0001966   m / s 2 , 0.0001135   m / s 2 ]
Station coordinates[−2852.90 km, 3399.95 km, 4565.25 km]
Three body gravitySun, moon, and major planets: JPL DE405
Observation error
(Gaussian distribution)
[ σ 1 , σ 2 , σ 3 ] = [ 0.03 km , 0.1 cos ( E l ) ° , 0.1 ° ]
TideSolid tide: IERS Conventions 2003
Non-spherical gravitational fieldGravity model: EGM2008
Degree:21; order:21
RelativityIERS Conventions 2003
Solar pressureShadow model: Dual Cone
Light pressure coefficient: 1.00
Area–mass ratio: 0.02 m2/kg
Atmospheric dragDensity model: Jacchia–Roberts
Drag coefficient: 2.20
Area–mass ratio: 0.02 m2/kg
Table 2. Solution results (UTC: 2023.04.02–2023.04.04).
Table 2. Solution results (UTC: 2023.04.02–2023.04.04).
Orbit ParameterSimulation ValueSolution ResultError
Tangential acceleration (m/s2)0.00019660.00019610.0000005
Normal acceleration (m/s2)0.00011350.00011060.0000029
Semi-major axis (km)6933.5345196933.8856320.351
Inclination (°)65.01164.9160.095
RAAN (°)14.37213.6290.743
Argument of latitude (°)28.72827.8940.834
Table 3. Results of simulation.
Table 3. Results of simulation.
Thrust1 N500 mN100 mN20 mN5 mN1 mN
Acceleration (m/s2)2 × 10−31 × 10−32 × 10−44 × 10−51 × 10−52 × 10−6
Simulation tangential acceleration (m/s2)1 × 10−35 × 10−41 × 10−42 × 10−55 × 10−61 × 10−6
Solution tangential acceleration (m/s2)1.000 × 10−35.000 × 10−40.999 × 10−42.033 × 10−54.645 × 10−60.486 × 10−6
Simulation normal acceleration (m/s2)1.732 × 10−38.660 × 10−41.732 × 10−43.464 × 10−58.660 × 10−61.732 × 10−6
Solution normal acceleration (m/s2)1.738 × 10−38.624 × 10−41.622 × 10−42.682 × 10−55.890 × 10−61.223 × 10−5
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Zhao, S.; Tao, X.; Li, Z. Initial Identification of Thrust and Orbit Elements for Continuous Thrust Spacecraft in Circular Orbit. Aerospace 2023, 10, 1012. https://doi.org/10.3390/aerospace10121012

AMA Style

Zhao S, Tao X, Li Z. Initial Identification of Thrust and Orbit Elements for Continuous Thrust Spacecraft in Circular Orbit. Aerospace. 2023; 10(12):1012. https://doi.org/10.3390/aerospace10121012

Chicago/Turabian Style

Zhao, Shuailong, Xuefeng Tao, and Zhi Li. 2023. "Initial Identification of Thrust and Orbit Elements for Continuous Thrust Spacecraft in Circular Orbit" Aerospace 10, no. 12: 1012. https://doi.org/10.3390/aerospace10121012

APA Style

Zhao, S., Tao, X., & Li, Z. (2023). Initial Identification of Thrust and Orbit Elements for Continuous Thrust Spacecraft in Circular Orbit. Aerospace, 10(12), 1012. https://doi.org/10.3390/aerospace10121012

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop