Atmospheric Density Model Optimization and Spacecraft Orbit Prediction Improvements Based on Q-Sat Orbit Data

Atmospheric drag calculation error greatly reduce the low-earth orbit spacecraft trajectory prediction fidelity. To solve the issue, the"correction - prediction"strategy is usually employed. In the method, one parameter is fixed and other parameters are revised by inverting spacecraft orbit data. However, based on a single spacecraft data, the strategy usually performs poorly as parameters in drag force calculation are coupled with each other, which result in convoluted errors. A gravity field recovery and atmospheric density detection satellite, Q-Sat, developed by xxxxx Lab at xxx University, is launched on August 6th, 2020. The satellite is designed to be spherical for a constant drag coefficient regardless of its attitude. An orbit prediction method for low-earth orbit spacecraft with employment of Q-Sat data is proposed in present paper for decoupling atmospheric density and drag coefficient identification. For the first step, by using a dynamic approach-based inversion, several empirical atmospheric density models are revised based on Q-Sat orbit data. Depending on the performs, one of the revised atmospheric density model would be selected for the next step in which the same inversion is employed for drag coefficient identification for a low-earth orbit operating spacecraft whose orbit needs to be predicted. Finally, orbit prediction is conducted by extrapolation with the dynamic parameters in the previous steps. Tests are carried out with the proposed method by using a GOCE satellite 15-day continuous orbit data. Compared with legacy"correction - prediction"method in which only GOCE data is employed, the accuracy of the 24-hour orbit prediction is improved by about 171m the highest for the proposed method. 14-day averaged 24-hour prediction precision is elevated by approximately 70m.


Introducton
Highly accurate low-earth orbit (LEO) operating satellite orbit prediction method are urgently demanded in nowadays as it has a wide range applications such as providing trajectory prediction for optical monitoring and satellite laser ranging [1], collision possibility analysis as spacecrafts and space debris in LEO increase rapidly [2]. Most of the LEO spacecraft orbit prediction errors result from inaccurate drag calculation especially atmospheric drag. It is common to use "correction-prediction" strategy to elevate prediction accuracy. For the step of "correction" , spacecraft orbital data is employed for inversion to correct parameters in the drag model such as atmosphere density, drag coefficient, and gas-surface interaction parameters. Usually, one parameter is fixed for correction of other parameters. In the step of "prediction", numerical integration is performed to obtain orbit extrapolation and the corrected parameters in previous step are applied [3]. However, coupling of the parameters in the atmospheric drag model results in convoluted uncertainties for the "correction-prediction" strategy [4]. The corrected parameter "absorbs" the uncertainties in other parameters. Generally, for spacecraft orbits at altitude of 400km-500km in the next 24 hours, the accuracy of the forecast using "correction-prediction" scheme can be reached between 30m to 60m. For the spacecraft at altitude of 300km-400km, the prediction accuracy is around 200m to 300m [5,6]. For lower altitudes, (less than 300km) trajectory prediction accuracy drops to 1000m to 1500m [1]. It can be observed that as the rarefied gas density increases thus atmospheric drag ascends, the efficacy of the "correction-forecast" strategy decreased. To unravel the interdependency, investigations were conducted to evaluate spacecraft drag coefficient individually. Nevertheless, factors such as temperature, molecular composition, degree of particles adsorption on the surface, surface material cause a virational drag coefficient [7,8]. It is reported that because of drag coefficient erroneous, atmospheric models there were developed by satellite orbital behavior overestimate or underestimate rarefied gas density tremendously [9]. To have a better estimation for the drag coefficient, diffuse reflection with incomplete accommodation [10], Cercignani-Lampis-Lord [11] and other scattering models are usually employed [12]. Further-more, selection of empirical atmosphere density model (EADM) are vital in drag calculation and therefore the spacecraft orbit forecast [5,13]. At present, EMADs that are commonly employed are, for example, Jacchia models, NRLMSISE-00, JB2008, DTM2000, etc. Each model perform differently and no model could cover all the aspects (solar flux, geomagnetic index, atmospheric winds, distribution of molecular concentrations, etc.) in space environment [14]. For instance, NRLMSISE-00 is less sensitive to the parameter such as solar activity, geomagnetic activity, latitude, and longitude than that of the actual atmosphere [15]. Jacchia models outperform the MSIS models when considering the impact of solar activity for orbital prediction [16]. For altitude less than 500km, JB2008 and Jacchia71 have a good performance in tasks like density representation and orbit prediction [5,17]. Revising a EADM is an important mean to improve orbit prediction accuracy. Inversion methods based on semi-major axis attenuation or non-conservative forces [18][19][20] are often performed by using a priori orbit data for atmosphere density recovery. However, such methods cannot provide atmosphere density forecasts [21]. Based on the dynamic approach-based inversion method, Wang proposed an technique for atmospheric density forecast [22]. The method was later improved by Zhang et al. [23] as thermal parameters are involved in. Recently, efforts has been devoted for diminishing the convolution effect as can be seen work by Ray et al. [24,25] who tried to estimating atmospheric density and drag coefficient simultaneously by spacecraft tracking data. Fourier series was employed for filtering uncertainties in the drag coefficient.
The Q-sat satellite, also known as Gravity and Atmospheric Science Satellite of xxx University, is a micro spherical satellite developed by the xxxxx Lab (xxxxx) at xxx University [26]. The satellite has a net weight of 21.2kg, a diameter of about 510mm. It operates in a sun-synchronous orbit at initial altitude of approximately 500km. Since it designed to be nearly spherical, the Q-sat can be considered to have a constant drag coefficient of 2.2 regardless its orientation based on molecular dynamics numerical simulation and experiments [17,22]. Dual frequency, carrier phase differential global navigation satellite system (GNSS) receiver is equipped and both Global Positioning System (GPS) and Beidou Navigation Satellite System (BDS) signal can be used for positioning. After post-processing, the orbit data could have accuracy of centimeter level. The xxxxx lab fetches the Q-Sat orbit data on a daily basis. Image of the Q-Sat with separation system is provided in Fig. 1. The separation system is intended to mount on the launch vehicle [27,28]. More detail of the satellite can be found in Ref. [26,27].
In current paper, an orbit prediction method which amis to improve prediction accuracy is proposed. In the first step, multiple EADMs would be modified by using Q-sat orbit data and a dynamic approach-based inversion. One of the EADM that performs the best would be selected. For demonstration of the proposed method, the Jacchia-Roberts atmospheric density model are modified in current paper. Next, by using the same inversion method, the drag coefficient for the LEO spacecraft which orbit needs to be forecast is identified. In this way, transmission error caused by using the same satellite orbit data can be reduced and the updated drag coefficient grounds an improved trajectory prediction accuracy. Tests are carried out by employing a Gravity Field and Steady-State Ocean Circulation Explorer (GOCE) satellite 15-day continuous orbit data. Obit prediction for GOCE satellite was conducted by proposed dynamic approach-based orbit prediction procedure. Legacy "correction-prediction" strategy was also employed for comparison. Discussions and conclusion are provided in the end.

Q-Sat Trajectory Based LEO Spacecraft Orbit Prediction Procedure
The proposed LEO spacecraft orbit prediction method can be outlined in 4 steps as shown in the following.
1. Modification of EADMs: A dynamic approachbased inversion method is first performed by using Q-sat orbit data for EADM revisions. Onboard dual-frequency, carrier phase differential GNSS receiver makes orbit determination centimeter level precision possible. Several EADMs would be modified in current step.
2. Selection of modified EADM: a optimized atmospheric density model that performs the best is selected.
3. Identification of drag coefficient: by involving the revised EADM and orbit data of the spacecraft whose trajectory needs to be predicted, the same dynamic approach-based inversion method employed on step 1 is performed again to obtain a drag coefficient for the spacecraft. The spacecraft orbit data for drag coefficient identification and the Q-Sat orbit data for EADM revision should from the same day.
4. Orbit forecast: based on the revised EADM and drag coefficient obtained in previous steps, orbit prediction is achieved by extrapolation.
It should be noted that because of the space environment variation, the obtained EADM and drag coefficient is considered to be valid for a short period of time.
Inversion for dynamic parameters should be carried out again if it is for orbit prediction long time after.

Inversion for EADM Optimization and Drag Coefficient Identification
A dynamic approached-based inversion would be introduced in this section. It can be employed for inverting different parameters individually by using spacecraft orbital data. Acceleration generated by various forces on a LEO spacecraft which orbits around the earth can be expressed as, where r is the satellite position vector,ṙ is satellite velocity vector,r is total acceleration vector for the spacecraft at epoch t, p is the dynamic parameter vector, and GM is the earth gravitational constant. The first term on the right is two-body gravitational acceleration, and the second term is the acceleration produced by the perturbation. The actual spacecraft motion state at epoch t can be represented as, and can be obtained from spacecraft orbit data. Spacecraft motion state obtained by a dynamic approachbased integration according to Eq. 1 can be presented as,X Motion state calculated by the integration is function of the initial motion state and dynamic parameters, Next, difference between satellite actual motion state and motion state obtained by the integration at epoch t n can be linearized as, where ∂X(t n )/∂X o is transtion matrix, it can be written by Φ tn and the matrix represents how the initial motion state error affect the subsequent motion states, ∂X(t n )/∂p is sensitivity matrix and it can be aliased by S tn . The sensitivity matrix shows the impact of dynamic parameter errors on the subsequent motion states [3]. Finally, corrections for initial state and dynamic parameters can be obtained by the least square method [17] since it is an over-determined system, where, EADM can be used to describe earth atmosphere status and how it changes [29], generally it is a function of temperature, atmospheric particle composition, solar activity, geomagnetic activity and other variables. For current investigation, attentions should be paid on parameter sensitivity thus the orbit prediction accuracy. In present paper, in order to present a EADMs modification process, Jacchia-Roberts model is chosen. It should be noted that in actual spacecraft orbit prediction improvement workflow, several EADMs should be modified and compared in terms of orbit prediction performance.
For temperature at height, h, above 125km, where T ∞ is exospheric temperature, T x is the temperature at 125km, R a is the polar radius of the earth and it is defined as R a = 6356.766km, and L is a corrective factor. Detailed descriptions can be obtain in paper by Roberts [30]. For current EADM modification, L is expressed as a polynomial function of T ∞ [31], in which l i is the dynamic model parameters and l = [l 1 , l 2 , · · ·, l 5 ]. Therefore, the sensitivity matrix for the acceleration by atmospheric drag is, where C d is the spacecraft drag coefficient, A is the cross-sectional area, m is the spacecraft mass, v r is spacecraft relative velocity. In Jacchia-Roberts model, atmospheric density can be obtained as [3,30,31], where ρ s is the calculated standard atmospheric density, ∆ρ is the corrected density under the effects of geomagnetism, half-year period, seasonal latitude, etc., and ∆ρ He is the corrected density for helium. Next, partial derivative respect to corrective factor L is performed as, The term of ∂ρ s /∂L can be expanded as, where i = 1, 2, · · ·, 5 and it represent N 2 , Ar, He, O 2 , and O correspondingly, ρ 6 is the density for hydrogen particles, α i is the diffusion coefficient, M i is molar mass of the gas component i, g o is the gravitational acceleration at sea level, and R is the gas constant, and finally for the last term in Eq. 13 which considers the hydrogen particles, where T s is the temperature at 500km according to Eq. 8 and, For the term of ∂ρ He /∂L which is for Helium density correction, Lastly, by substituting the equations above and, in Eq. 10, the sensitivity matrix, S(t), for atmospheric drag acceleration is obtained. The input for the proposed dynamic approach-based inversion is a spacecraft orbital data. In order to optimize a EADM, drag coefficient is fixed to a certain value and corrective parameters for the EADM would be on the output. The over-determined system is solved by least square fitting as shown in Eq. 6. For demonstration of EADM correction in current paper, corrective factors l i , i = 1, 2, · · ·, 5 for Jacchia-Roberts are on the output of the inversion. The practice is especially considered as appropriate for orbit data from spherical satellite, for example the Q-Sat launched by xxxxx at xxx University as its drag coefficient is not relative to its inclination (angle of attack). After EAMDs are revised and the optimal EAMD is selected, the tracking data of spacecraft which trajectory needs to be predicted is employed as the input of the inversion and revised drag coefficient for the spacecraft will be on the output. Lastly, the optimized EADM and drag coefficient are involved for orbit extrapolation.

EADM Optimization
For demonstration for EADM optimization, the Q-Sat orbit data is employed and data for every 24-hour is considered as a data set. A consecutive 5 days data that expands from September 16 to September 20, 2020 is used. After post-processing the orbit data from Q-Sat on-board dual frequency, carrier phase differential GPS/BDS receiver, the obtained orbit data could have centimeter-level accuracy. In present paper, corrective factors, l i , i = 1, 2, · · ·, 5, for Jacchia-Roberts atmospheric model is optimized by the proposed dynamic approach-based inversion method. Five set of corrective factors are obtained as shown in Tab. 1 By applying the corrective factors, atmospheric density along the Q-Sat orbit is calculated and plotted in Fig. 2. At 16:14:20 (hh:mm:ss) on September 20th, the density is the highest at 1.6807 × 10 −13 kg/m 3 during the 5 days. At the time, the Q-Sat height is about 495.6604km, the geomagnetic index is K p ≈ 2.3nT , and the daily average value of solar radiation index is F 10.7 ≈ 71 × 10 −22 W/(m 2 · Hz), which is the largest in 5 days. At 00:17:50 on September 18, 2020, the density is 2.0842 × 10 −14 kg/m 3 which is the smallest in the 5 days. At the time, the Q-Sat height is about 495.6721km, the geomagnetic index is approximately 0nT , which is the smallest in 5 days, and the daily average value of solar radiation index is about 70 × 10 −22 W/(m 2 · Hz), only slightly higher than 69 × 10 −22 W/(m 2 · Hz) on Sept 16th.
In order to examine the proposed inversion for EADM optimization, an arclength of 27 hours orbit was  employed and atmospheric densities calculated by optimized EADM are compared. For the first 27 hours, orbit data from the first day was added with the first 3 hours data on the second day for a total of 27 hours. For the second 27-hour, orbit data from the last 3 hours of the first day was added at the beginning of the second day. Next, the 27-hour data is employed for EADM optimization. The atmospheric densities for the last 3 hours of the first day calculated by EADMs that are optimized by 27-hour data and by 24-hour data were compared. The same comparisons were done for the atmospheric densities for the first 3 hours of the second day. As a result, there are in total of eight 3-hour rarefied gas density intervals for comparison from September 16 to September 20 and maximum relative errors and mean relative errors over the 3-hour are presented in Tab. 2. As one can observe in the table, in the 8 comparison intervals, the maximum relative error is a little less than 18%, the smallest is around 5%. For the mean relative errors, the largest is around 15% at the smallest is around 4.3%. The total average of the mean relative errors is about 8%.

Drag Coefficient Identification
For current step, the input for the dynamic approachbased inversion is the orbit data of the LEO spacecraft which orbit prediction is on demand, the drag coefficient for the spacecraft is on the output. Revised EADM carried out by Q-Sat orbit behavior is employed in the inversion for LEO spacecraft drag coefficient identification. Not to mention that the revised EADM shall only be used for identifying drag coefficient at the same day. The revised EADM and redefined spacecraft drag coefficient are produced in a daily basis in order to include variations and anomalies in an entire day and would not be accurate for other days with different space environment. This practice lead to a potentially better prediction accuracy. Since it is hard to reach orbital information of a LEO spacecraft that is currently orbiting the earth, for demonstration of the proposed orbit prediction improvement method, the GOCE satellite orbit data is employed. The satellite data is open to the public and can be downloaded from the European Space Agency website. The on-board dual frequency GPS receiver provided position data every 10 seconds and the data could have centimeter level accuracy with postprocessing.

Space Environment Parameters
A 15-day continuous orbit for GOCE satellite is used. The arc extents from September 6th to September 21st in 2010 in which the on-board ion propulsion engine was not active for drag-free flight and the satellite was descending. The arc of first day is utilized for drag coefficient identification and the arc for the second day is employed for examining the fidelity of the orbit prediction. Drag coefficient and space environment parameters obtained in the previous day are used in orbit propagation for the current day as the space situation of the current day often show a good match with those of the previous day [1]. Since the Q-Sat was launched lately in 2020 and the orbit data is employed for EADM optimization, to implement drag coefficient identification for GOCE satellite with orbit data in 2010, the space environment for Q-Sat orbit data this is used for EADM optimization should be analogous to that of GOCE [32]. For the Jacchia-Roberts atmospheric model that is modified in previous section, the geomagnetic index, K p , and the solar radiation index, F 10.7 , are the variables for calculating the exospheric temperature, T ∞ . The GOCE arc for drag coefficient identification should have similar values of k p and F 10.7 with the Q-Sat arc which is employed EADM optimization. If it is for prediction of a spacecraft that is currently orbiting, this step is unnecessary and it is only needed for present demonstration. Since usually it is hard to obtain a strict agreement for environmental values for two arcs that have difference of 10 years, approximation should be made. Daily averaged K p and F 10.7 are considered for GOCE arc. Arc for Q-Sat with comparable values is needed to be found. The difference for the daily averaged value for K p and F 10.7 should be no larger than 0.6 and 15 respectively. In Fig. 3 diation 10 year ago is overall stronger than in year of 2020.
In present paper, as data from satellite that is currently orbiting is hard to reach, GOCE satellite is employed for purpose of demonstration, the orbit prediction procedure is complexified. In real orbit prediction task using the proposed method, environmental parameters matching shown in this section is not required.

Orbit Data Length for Drag Coefficient Identification
The dynamic approach-based inversion is employed for drag coefficient identification and the input for the inversion is the orbit data of the spacecraft which in needs of trajectory prediction. Optimized EADMs by corresponding Q-Sat orbit data (Tab. 3) is implemented in the inversion. To save computational expense and produce a rapid prediction, arclength of the last 1.5, 3,6,9,12,15,18,21,24 hours of the first day are tested for finding an optimal trajectory length for drag coefficient identification. In Fig. 5

Spacecraft Orbit Prediction
With the identified drag coefficient and optimized EADM, orbit predictions are conducted by extrapolation with the dynamic parameters employed for drag coefficient identification. The prediction extend for 24 hours and the predicted orbit is compared with the ac-tual orbit. The 3D error, is calculated where i = 1, 2, 3 which represents the along-track, cross-track and radial direction and subscript 'p' represents for predicted orbit and 'a' for actual orbit. Maximum 3D prediction errors in the 24hour interval are found and plotted in Fig. 6. For comparison, maximum 3D error for orbit prediction with original Jacchia-Roberts model are also included. As can be observed in Fig. 6, prediction with optimized EADM has overall decreased maximum 3D errors. It indicated an overall higher predication accuracy. The difference between the two prediction, ∆e 3D , are plotted in Fig. 7. The GOCE satellite 24-hour prediction accuracy is increased by 171m at the highest. The averaged 24-hour prediction accuracy over 14 days is increased approximately 70m. It can be concluded that the accuracy of the orbit prediction for GOCE satellites was improved significantly.

Conclusion
To diminish the convoluted uncertainties in the correction prediction strategy which based on single spacecraft orbital behavior, a Q-Sat trajectory based LEO spacecraft orbit prediction method is proposed in present paper. The Q-Sat satellite, launched by xxxxx lab at xxx University which initially orbit around 500km, is a spherical micro satellite for gravity field recovery and atmosphere density detection. The satellite is considered to have a constant drag coefficient of 2.2 regardless its attitude and its tracking data is employed for EADM modification by dynamic approachbased inversion introduced in present paper. The revised EADM is later employed in drag coefficient identification process for the spacecraft which orbit needs to be predicted. Convoluted error is thus reduced as EDAM modification and drag coefficient recognition are decoupled. As can be seen in many research that, in actual space, spacecraft drag coefficient changes due to variations such as composition of gas, gas-surface interaction, molecules absorption, etc. In current work, drag coefficient of the GOCE satellite is calculated in a daily averaged fashion, which contains the uncertainties occured in an entire day. As a result, the identified drag coefficient of the GOCE satellite is no longer identical from a day to another. So as to revised EADM that varied from day to day. The obtained drag coefficient and EADM with space environment parameters from the previous day are then applied in the trajectory forecast for next 24 hours as the space environ-ment can be considered similar as of the previous day. For demonstration of the proposed method, space environment parameters k p and F 10.7 for arc of GOCE and Q-Sat should be comparable in order to obtain a higher prediction acurracy. In real prediction workflow, this step in not necessary. However, forecast of the space environment parameters for the next day is of also important for improving prediction accuracy. Otherwise, values from previous day will need to be employed. In present work, the GOCE satellite 24-hour prediction accuracy is increased by 171m at the highest compared to result by legacy "correction-prediction" strategy. The 14 days averaged 24-hour prediction accuracy is increased approximately 70m.

Data Availability
The datasets employed during the current study are available from the corresponding author on reasonable request.