International Journal of Aerospace Engineering
Volume 2019 (2019), Article ID 1758956, 15 pages
https://doi.org/10.1155/2019/1758956

Research Article

A Spherical Micro Satellite Design and Detection Method for Upper Atmospheric Density Estimation

Zeyang Zhao,¹ Zhaokui Wang,² and Yulin Zhang1,2

¹College of Aeronautics and Astronautics, National University of Defense Technology, Changsha 410073, China
²School of Aerospace, Tsinghua University, Beijing 100084, China

Correspondence should be addressed to Zhaokui Wang; wangzk@tsinghua.edu.cn

Received 3 June 2018; Accepted 22 November 2018; Published 28 January 2019

Guest Editor: Vaios Lappas

Copyright © 2019 Zeyang Zhao et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

A spherical microsatellite, aimed at technology verification of upper atmospheric density detection and precision orbit determination in low Earth orbit (LEO), has been developed by Tsinghua University and is scheduled to launch in the second half of 2019. In order to reduce the influence of other nonconservative forces, the area-to-mass ratio of the satellite was designed to ensure that the aerodynamic drag of the satellite is more than one to two orders of magnitude greater than the solar pressure perturbation. The influence of the components of the antenna, separation support, and solar cell arrays (SCAs) on the spherical structure for the area-to-mass ratio is reduced by reasonable design of the satellite flight attitude. In addition, a novel method has been developed to estimate and correct the parameters for atmospheric density (correlated to the drag) calculation, using precise orbit data (POD) obtained from an onboard dual-frequency global position system (GPS) receiver. The method is used to obtain the partial derivative of the total acceleration of the satellite to the initial state (position, velocity) and Harris-Priester (HP) density model parameters (antapex and apex density). Further, the least squares estimation method is used to solve the overdetermined linear equations and obtain the change of initial state of the satellite and the HP density model parameters relative to change of the satellite state. The validity of this method has been verified through numerical simulations with the parametric setup of the satellite. The estimated density precision was found to be higher than two to three orders of the initial HP density model with continuous iteration.

1. Introduction

Estimation and correction of upper atmospheric semi-empirical density models have attracted substantial interest over recent decades due to their wide applications for the prediction of spacecraft motion, precise orbit determination, debris collision avoidance, and location estimation of re-entry spacecraft. Density decreases exponentially with orbit height but is affected by sunshine and geomagnetic activity, among others, resulting in uncertainties of 10–30% [1] and as much as 100%. Indeed, Vallado and Finkelman [2] suggested that no single atmospheric model is optimal for all applications. Density model uncertainty has adversely affected some low-orbit missions; hence, substantial in-orbit experimental data is still required to improve upper atmosphere empirical models.

Currently, empirical models of upper atmospheric density mainly rely on in-orbit observation data from upper atmospheric density detection satellites such as Atmospheric Neutral Density Experiment (ANDE) [3], Drag and Atmospheric Neutral Density Explorer (DANDE) [4, 5], SpinSat [6], and QB50 missions [7]. A common feature of upper atmospheric density detection satellites is that they include a change in the area-to-mass ratio along the track direction, which can decrease the density error caused by their aerodynamic shapes. In addition to atmospheric detection satellites, near-Earth space science mission data can be used to estimate and correct atmospheric density, including Challenging Minisatellite Payload (CHAMP) [8], Gravity Recovery and Climate Experiment (GRACE) [9], and Gravity Field and Steady-State Ocean Circulation Explorer (GOCE) [10], for which precise orbit data is measured by onboard global position system (GPS) receivers and the nongravitational acceleration vector (sum of the atmospheric drag, solar pressure, and third-body perturbations) measured by accelerators.

Whether using atmospheric density detection satellites or near-Earth space science mission data, the main payloads for atmospheric density detection include mass spectrometers [11–13], accelerometers [14–16], and orbit data (obtained by GPS, two line orbit elements (TLEs) [17, 18], and satellite laser ranging (SLR) [19]), and others such as drag balance instruments which have flown on the San Marco satellite Missions [20] and UNICubeSAT [21, 22]. A mass spectrometer can obtain the composition of rarefied gas particles in the upper atmosphere around the spacecraft; it is typically used in conjunction with an accelerometer (used to identify non-gravity perturbations of satellites in orbit) to determine the density and perturbation of the atmospheric drag. The advantages of these action items are higher resolution and better real-time performance, but they require precise attitude control and are more expensive.

Using satellite orbit data to determine atmospheric density is also important for atmospheric density estimation and correction. Doornbos [23] established the concept of the scale factor. The scale factor is multiplied by the atmospheric density as its error correction term. The scale factor is expanded by the spherical harmonic series in latitude and local solar time, and the coefficients of the spherical harmonic series are solved by TLE data of different satellites; thus, estimation and correction of the large air tightness model are achieved. However, the resolution of TLE data is usually low. Each satellite collects only a few sets of data every day, and the orbital accuracy is of the order of 100–200 m; therefore, it can only reflect the average density under a certain temporal resolution. In response to this problem, precision orbit data (of cm resolution) acquired by SLR or GPS has been used to estimate the atmospheric density in recent years. Because of the high frequency of data recording, it has a very high temporal resolution. For example, Jeon et al. [24] used precise orbit data obtained by the SLR of Starlette satellite to estimate atmospheric density.

In this study, we used the aerodynamic shape design of a spherical satellite and its precision orbit data obtained by onboard GPS receivers of commercial off-the-shelf products to propose a method for estimating and correcting the upper atmospheric density model and verified the proposed method via numerical simulation. The magnitude of the atmospheric drag and the solar pressure at a 400–600 km altitude with different area-to-mass ratios is analyzed in the next section and provides reference for aerodynamic shape design of spherical satellites. The design of a spherical satellite shape to ensure a small error caused by aerodynamic shape and attitude change is outlined in Section 3. In Section 4, the least squares modeling method for estimating and correcting the density model is presented, as well as simulation validation with parametric setup of the spherical satellite. The conclusions of this study are presented in Section 5.

2. Perturbation Analysis for Density Estimation

The motion equation for a satellite in low orbit can be written as

where the acceleration terms on the right-hand side can be grouped into two parts, where represents gravitational acceleration and represents nongravitational acceleration. The gravitational acceleration terms include central gravitational acceleration, nonspherical gravitational perturbation, and others, including gravitational effects of third bodies (sun, moon, and planets) and the effects of general relativity. The nongravitational acceleration terms mainly include aerodynamic acceleration, solar radiation pressure acceleration, and others such as Earth’s albedo, infrared radiation pressure acceleration, and thrusters. The mathematical representation of solar radiation and atmospheric drag perturbation is as follows:

where represents the aerodynamic acceleration, is the drag coefficient, is the projected area of the satellite along the flight direction, is the satellite mass, is the incoming particle velocity (relative to the satellite), and is the mass density. The Harris-Priester model, which is briefly depicted in Supplementary Materials (available here), is adopted in this study. For radiation pressure acceleration , is the radiation pressure force coefficient, is the solar pressure from the sun in a blackbody of one astronomical unit, is the exposure factor, and is the unit vector from the satellite to the sun. It can be seen from Equation (2) that, for a near-earth spacecraft, the variation in the magnitude of solar pressure can be assumed to only be related to the area-to-mass ratio of the satellite, and the variation in the magnitude of the atmospheric drag is related to the orbit height and the area-to-mass ratio of the satellite. To determine the influence of the orbit height to the area-to-mass ratio on the atmospheric drag and the solar pressure, we analyzed the magnitude of the atmospheric drag and solar pressure at 400–600 km orbit with different area-to-mass ratios.

The atmospheric drag and the solar pressure perturbation (calculated using Equation (2)) with an area-to-mass ratio of 0.001 (red dashed line), 0.005 (blue dash-dotted line), and 0.01 (green circle) are calculated using simulation parameters in Table 1 and shown in Figure 1. The simulation results show that the magnitude of the solar pressure perturbation has almost no effect on the 400–600 km orbit height. For the 600 km orbit, the atmospheric drag perturbation is equal to the magnitude of the solar pressure perturbation. For other heights, the magnitude of the atmospheric drag perturbation is over one to two orders of magnitude higher than that of the solar pressure perturbation (Figure 1). Therefore, in orbits below 600 km, the influence of solar pressure on satellite motion can be neglected.

Table 1: Simulation conditions.

Simulation time (s)
Initial orbital element

Orbital height (km)

Eccentricity
Inclination (°)
Right ascension of ascending node (°)
Argument of perigee (°)
True anomaly (°)

1500
412.18
0
96.8213
340.251
0
0

512.18
0
96.8213
340.251
0
0

612.18
0
96.8213
340.251
0
0

(a)

(b)

(c)

(d)

Figure 1: Atmospheric drag at heights of (a) 400 km, (b) 500 km, (c) 600 km, and (d) solar pressure for different area-to-mass ratios.

3. Spacecraft Design for Density Detection

In this section, we outline the design of a spherical microsatellite, including the platform and primary payload, considering the influence of separation supports, antennas, and the installed SCAs on the aerodynamic shape of the satellite to ensure that the area-to-mass ratio experiences minimal changes.

3.1. Mission Overview

The satellite (mass: 23 kg, diameter: 626 mm, height: 510 mm) has a payload of a dual-frequency commercial off-the-shelf (COTS) GPS receiver (Figure 2(a)) to verify the precision orbit determination technology (Figure 2(b)) and upper atmospheric density estimation (Figure 2(c)) without an accelerometer and achieve long-wave gravity field recovery (Figure 2(d)) and as a pre-research platform for the future inner-formation flying system (IFS) [25–27]. A flywheel and three magnetorquers were adopted for three-axis attitude stability for the earth. Both uniform S-band (USB) and industrial, scientific and medical (ISM) bands are adopted for telemetry, telecontrol, and data transmissions. The satellite will be launched in the second half of 2019 as the secondary payload in a 500 km sun synchronous orbit.

(a)

(b)

(c)

(d)

Figure 2: (a) Dual-frequency GPS, (b) precision orbit determination technology based on GNSS, (c) upper atmospheric density detection, and (d) long-wave gravity field detection.

3.2. Spacecraft Platform and Payload Design

The satellite [28] is a spherical microsatellite, and its main components are shown in Table 2 and its exploded-view drawing shown in Figure 3.

Table 2: Main components of the satellite.

Subsystem
Main components
Key specifications

Structure
Hemispherical frame (HF), payload cabin (PC), central connection ring (CCR), releasing device (RD)
2A12-H112 aluminum alloy

Payload
Dual-frequency GPS (D-GPS), single-frequency GPS (S-GPS), micro camera (MC)
Real-time positioning (velocity) precision: less than 20 m (0.2 m/s)

Attitude control
Inertial measurement unit (IMU), integrated full field sun sensor (IFFSS), magnetometer, flywheel, magnetorquer
Attitude determination precision: 0.5°

Attitude control precision: 2°

Power
Solar cell array (SCA), power controller unit (PCU), battery (BAT)
SCA: 931 cells of triple-junction GaAs (29.5%)

BAT capacity: 8.8 Ah

Bus voltage: 6.9–8.2 V

Secondary voltage: 12 V, 3.3 V, 5 V

Power consumption: 11.45 W (steady state)

Integrated electronic
Onboard computer (OBC), attitude control computer (ACC), integrated electronic circuit board (IECB)
Processing chip: LPC2294

Telemetry, telecontrol, and data transition
Unified S-band (USB), industrial scientific medical (ISM) band, power amplifier (PA), combiner
USB: Telemetry: 2000 bps

Telecontrol: 4096 bps

Data transition: 0.5 Mbps

ISM: 115.2 kbps

Figure 3: Exploded view of the satellite.

The satellite uses a central load-bearing device (CCR and PC) with two hemispherical frames to assemble the configuration. A releasing device enables separation of satellite and rocket via a locking and ejecting mechanism on the CCR. The heavier components and the large power consumption devices are directly installed in the ninth palace configuration, which improves the mechanical performance and heat dissipation environment. The light components on the satellite are connected through the IECB (Figure 4), which is responsible for the management of satellite power supply control, information collection of component and payload, and data processing. This design eliminates the need for complex cabling on the satellite and improves reliability. COTS technologies and industrial-grade devices are used in the assembly of the satellite.

Figure 4: D-GPS.

The primary payload of the satellite is the 80 mm long, 52 mm wide, 20 mm high, dual-frequency COTS GPS receiver (Figure 4). The microstrip antenna (Figure 5) is placed on a pentagon plate and is covered with a glass steel to resist space radiation. The position precision for navigation is less than 10 m, and the precision orbit could be less than 10 cm based on the carrier and pseudo range data. An adaptive design of the space environment for GPS has already been completed and verified via a thermal vacuum test (Figure 6) to ensure the capability of the GPS.

Figure 5: Antenna of D-GPS.

(a)

(b)

Figure 6: Thermal vacuum test for D-GPS.

3.3. Aerodynamics Configuration Design

Although the satellite is spherical, its characteristics will be diluted owing to separation of the support, antenna, and SCAs. In this section, the influence of the satellite antenna and separation support on the surface quality ratio of the satellite is analyzed. The flying attitude design and SCA installation are considered for reducing the error of the area-to-mass ratio caused by the aerodynamic configuration of the satellite. Part of the spherical satellite design has been granted patent authorization [29]. Figure 7(a) shows that the components affecting the spherical configuration of the satellite include the antenna and the separation support, which influences the area-to-mass ratio. In order to reduce the deviation of the satellite area-to-mass ratio, a method of adjusting the flight attitude of the satellite is adopted so that the separation support part of the antenna contributes to the windward area. The satellite axis points to the ground while the flight is stable, points toward the flight direction and and , and follows the right-handed rule along the negative direction of the orbital plane (Figure 7(b)). This design can reduce the influence of the antenna and the separation support on the spherical structure of the satellite.

(a)

(b)

Figure 7: (a) Components at satellite surface and (b) flying attitude of the satellite.

The SCAs are another major factor affecting the spherical structure of the satellite. In the initial design of the satellite, the installation of the soleplate was designed as a plane (Figure 8(a)), considering the installation of the SCAs. Subsequently, the battery was installed directly on the spherical surface and the SCA installation floor was redesigned (Figure 8(b)). The concrete method was first used to attach the polyimide film to the insulation at the ball floor (Figure 8(c)); then, the battery sheet was placed on the spherical surface (Figure 8(d)). In this manner, the spherical structure of the satellite is guaranteed.

(a)

(b)

(c)

(d)

Figure 8: (a) Prototype of satellite, (b) new surface for the satellite, (c) the polyimide film, and (d) a sample of the SCAs.

4. Dynamic Method of Atmospheric Density Model Estimation and Correction

When a low-orbit satellite is in orbit, the inaccuracy of the empirical atmospheric density model leads to deviation between the real orbit and the reference orbit integrated with the dynamic model. The current empirical model is mainly based on the hypothesis of hydrostatic balance, and the parameters of model are fitted by track data, mass spectrometer data, and acceleration data. In this section, a method for modifying the parameters of the empirical density model based on the precision orbit data is proposed, and the parameter correction model of the Harris-Priester (HP) density model is derived using this method. Finally, the validity of the corrected model is verified by a simulation test.

4.1. Estimation and Correction Model

Montenbruck and Gill [30] proposed the parameter estimation-based orbit determination. In the J2000 coordinate system, the satellite state vector including the position and velocity is shown as

The satellite state vector is also relative to the initial state, dynamics parameters, and time and is defined as

According to Figure 9, the real orbit can be assumed a precise orbit with cm-scale errors obtained by processed GPS precision data that reflect the actual perturbation model. The reference orbit is integrated with the dynamic model.

Figure 9: Comparison of the reference orbit and real orbit.

At time , the difference between and can be written as

reflects the bias of the reference model and the real model, so the equation can be written as

where is the state transition matrix and is the parameter sensitivity matrix. The satellite dynamic equation can be written as

where is the position of the satellite, is the velocity of the satellite, and is the vector of the parameter. We derive to as follows:

We define , , and .

Further, Equation (8) can be written as

We define , and Equation (9) can then be written as

We then create variational equations:

The derivation equation (Equation 12) at time is as follows:

where the initial condition is , and and can be calculated by solving the variational equations (Equation 12). After obtaining the state transition matrix and parameter sensitivity matrix, can then be obtained using observation Equation (13).

In order to achieve the effective solutions, the number of observation equations should be much greater than the number of unknowns to be solved. Then, the solution to Equation (13) needs to be used in the least squares method to solve the overdetermined linear equations in Equation (14):

4.2. Parameter Correction Model for the HP Density Model

The parameters of density are obtained from the density model and influence the accuracy of the model. The precise orbit data can be applied to derive density data by estimating the aerodynamic scale factor parameters [31], which assume that errors in the measurements, measurement modeling, and force models are negligible, and the effect of the true density on the drag force will be close to that of the scale factor times the modeled density over the time span. The second method of parameter estimation uses empirical acceleration parameters, which are considered as the whole non-gravity accelerometers [32]. The third method employs accelerometer calibration parameters and determines the scale and bias parameters for each axis using precise orbit data [33]. All of these methods do not analyze the parameters in the density model itself, as with the tabulated minimum and maximum density values and in the HP density model, which are used to obtain the antapex density and the apex density by exponential interpolation in the HP density model. Estimation of the and method is proposed here.

For a given orbit altitude, , the tabulated minimum and maximum density values, , , , and , can be used to compute the coefficient matrix, , , and of the variational equations, which can be written as

where can be further expanded to

Then,

where can be further expanded to

and can be further expanded by

Then,

Thus,

Equations (15), (16), (17), (18), and (19) describe the solution to the coefficient matrices , , and of Equation (12). Equation (12) can be further solved to obtain the state transition matrix and the parameter sensitivity matrix , and finally the correction value of atmospheric density parameters could be obtained by the least squares method.

4.3. Simulation Validation for Correction of the HP Model

The real orbit should be represented by observation data produced by precision orbit data. However, in order to validate the method of parameter estimation and correction, we assume the following: the real orbit is an orbit added into the error-based reference orbit; the nonspherical and atmospheric drag perturbations are considered both in the reference orbit and in the real orbit (additionally, three-body gravitational perturbations and solar radiation perturbation are considered as a constant disturbance force of 10⁻⁹ m/s²); the initial HP density model in the real orbit is accurate; the same magnitude bias is added to the HP density model for the reference orbit as to the initial HP density model for the real orbit; a fix drag coefficient [34] of 2.2 is generally used; and an area-to-mass ratio of the satellite of 0.008 m²/kg is adapted for the simulations (the design parameters of the spherical satellite).

Because calculation of the dynamic inversion method is very large and time-consuming, the parallel computation method is considered. The computation progress is shown in Figure 10. First, an orbit only considering the center gravity is achieved and employed as the initial value of the reference orbit and real orbit. Second, the reference orbit and the real orbit are generated by the dynamic models. Third, we obtain the residual of the real orbit and reference orbit as well as the variational equations generated by the reference orbit for computation of the norm equation, whose purpose is to obtain the correction of the density parameters. The iteration continues until the bias density model is close to the initial HP density model. The simulation conditions are listed in Table 3. The 5 cm and 10 cm orbit errors are considered here as they are general precision levels for orbit determination [35], and the numbers of iterations are 6 and 8, respectively. The correction results are given in Tables 4 and 5 and Figures 11 and 12.

Figure 10: Computation progress of the dynamic inversion.

Table 3: Simulation parameters for density estimation.

Orbit height (km)
Error of orbit determination (cm)
Sample time/total time (s)
Order of the gravity field
Arc length (s) (parallel program)
Total simulation times (s)
Number of iterations

512
5
1
30
400
9600
6

8

10
1
30
800
14400
6

8

Table 4: Comparison of corrected parameters and initial parameters for an arc length of 400 s and a height of 500 km.

(10⁻¹² kg/m³)
(10⁻¹² kg/m³)
(10⁻¹² kg/m³)
(10⁻¹² kg/m³)

Bias HP model
0.8916
1.5042
0.7819
1.105

No. of iterations
Corrected HP
Initial HP [29]
Corrected HP
Initial HP [29]
Corrected HP
Initial HP [29]
Corrected HP
Initial HP [29]

6
0.7531
0.3916
1.9648
2.042
0.2861
0.2819
1.9127
1.605

8
0.3927
2.0369
0.2823
1.6081

Table 5: Comparison of corrected parameters and initial parameters for an arctime of 800 s and a height of 500 km.

(10⁻¹² kg/m³)
(10⁻¹² kg/m³)
(10⁻¹² kg/m³)
(10⁻¹² kg/m³)

Bias HP model
0.8916
1.5042
0.7819
1.105

No. of iterations
Corrected HP
Initial HP [29]
Corrected HP
Initial HP [29]
Corrected HP
Initial HP [29]
Corrected HP
Initial HP [29]

6
0.5073
0.3916
2.2606
2.042
0.3384
0.2819
1.6452
1.605

8
0.3916
2.0403
0.2820
1.6063

(a)

(b)

(c)

(d)

Figure 11: Initial density model (red dotted line) and the corrected density model (green dashed line) indicating the corrected density minus the initial density with (a and b) 6 and (c and d) 8 iterations for a 5 cm orbit error.

(a)

(b)

(c)

(d)

Figure 12: Initial density model (red dotted line) and the corrected density model (green dashed line) indicating the corrected density minus the initial density with (a and b) 6 and (c and d) 8 iterations for a 10 cm orbit error.

The parallel calculation is used here in order to improve the effectiveness of calculation. The total simulation times were divided into some equal time (arc length) which is assigned to different CPU processors for solving the orbital integral and differential equations. Arc length is an important computation parameter, because atmospheric correction is a nonlinear process, and currently we can only use the linear approximation method. Arc length has an influence on the linearization error and thus the corrected density model. Here, different arc lengths were considered to reflect the performance of the corrected density model. First, the 5 cm orbit error was considered and the initial bias value for the HP model and the final correction of the density parameter for six and eight iterations are shown in Table 4. From Figures 11(a) and 11(b), we can see that increasing the iteration enhances the correction effect of the density model. Figure 11(c) indicates that the corrected density model and the initial density model are the same magnitude, and Figure 11(d) indicates that the initial density model is smaller by two orders of magnitude after eight iterations. Thus, the performance of the corrected parameter is improved by increasing the amount of iterations. For the 10 cm orbit error of orbit determination, the arc length becomes 800 s, and the correction of parameters with six iterations is clearly better than that for the 5 cm orbit error (Figure 12(c)). Figure 12(d) indicates that the initial density model is smaller by three orders of magnitude after eight iterations. Thus, a longer orbit integrated time enhances the performance of correction. The initial bias value for the HP model and the final correction of the density parameter for six and eight iterations are shown in Table 5.

5. Conclusions

In this paper, the design of a spherical microsatellite, which is aimed at upper atmospheric density and precision orbit determination, was presented. The satellite will be launched in 2019 as a secondary payload. The rational aerodynamic shape and flight attitude design of the satellite reduces any interference of the area-to-mass ratio change in atmospheric density estimation precision. Thus, using a spherical microsatellite to detect the upper atmosphere density is superior to the traditional satellite. The precision orbit data obtained by GPS can improve the temporal resolution of the density model estimation. The estimation and correction method achieves density detection, without accelerometer data, and establishes the relationship between the orbit change and density model parameters and numerical simulation results, illustrating the effectiveness of the proposed method. Note that the orbit determination error, arc length, and number of iterations affect the performance of parameter correction in the HP density model. From the perspective of orbit simulation, the magnitude of the orbit error directly affects the number of iterations; for the same number of iterations, a longer integral arc length produces a better model correction result.

Data Availability

The Harris-Priester model data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This satellite mission and research was supported by the National Natural Science Foundation of China (No. 11002076). The author is grateful to the many colleagues and associates who helped with this research.

Supplementary Materials

The expression of the Harris-Priester (HP) model which is mentioned in Dynamic Method of Atmospheric Density Model Estimation and Correction. (Supplementary Materials)

References

F. A. Marcos, B. R. Bowman, and R. E. Sheehan, “Accuracy of Earth’s thermospheric neutral density models,” in AIAA/AAS Astrodynamics Specialist Conference, Keystone, CO, USA, August 2006.

D. A. Vallado and D. Finkleman, “A critical assessment of satellite drag and atmospheric density modeling,” Acta Astronautica, vol. 95, pp. 141–165, 2014.

A. C. Nicholas, G. C. Gilbreath, S. E. Thonnard et al., “The Atmospheric Neutral Density Experiment (ANDE) and Modulating Retroreflector in Space (MODRAS): combined flight experiments for the space test program,” in Proceedings Volume 4884, Optics in Atmospheric Propagation and Adaptive Systems V, Crete, Greece, March 2003.

M. D. Pilinski and S. E. Palo, “An innovative method for measuring drag on small satellites,” in 23th Annual AAIA/USU Conference on Small Satellite, Logan, Utah, 2009.

M. D. Pilinski, “Analysis of a novel approach for determining atmospheric density from satellite drag,” Tech. Rep., MSc Thesis. University of Colorado, 2008.

A. Nicholas, T. Finne, and J. Galysh, “SpinSat mission overview,” in 27th Annual AAIA/USU Conference on Small Satellite, Logan, Utah, September 2013.

J. Thoemel, F. Singarayar, T. Scholz et al., “Status of the QB50 cubesat constellation mission,” in 65th International Astronautical Congress, Toronto, Canada, January 2014.

P. Moore, J. F. Turner, and Z. Qiang, “CHAMP orbit determination and gravity field recovery,” Advances in Space Research, vol. 31, no. 8, pp. 1897–1903, 2003.

Z. Wei, H. Hou-Tse, Z. Min, and Y. Mei-Juan, “Accurate and rapid error estimation on global gravitational field from current GRACE and future GRACE follow-on missions,” Chinese Physics B, vol. 18, no. 8, pp. 3597–3604, 2009.

D. Muzi and A. Allasio, “GOCE: the first core earth explorer of ESA’s earth observation programme,” Acta Astronautica, vol. 54, no. 3, pp. 167–175, 2004.

A. O. Nier, W. E. Potter, D. R. Hickman, and K. Mauersberger, “The open-source neutral-mass spectrometer on Atmosphere Explorer-C,-D, and -E,” RadioScience, vol. 8, no. 4, pp. 271–276, 1973.

J. H. Hoffman, W. B. Hanson, C. R. Lippincott, and E. E. Ferguson, “The magnetic ion-mass spectrometer on Atmosphere Explorer,” Radio Science, vol. 8, no. 4, pp. 315–322, 1973.

H. B. Niemann, “An atomic oxygen beam system for the investigation of mass spectrometer response in the upper atmosphere,” The Review of Scientific Instruments, vol. 43, no. 8, pp. 1151–1161, 1972.

J. T. Emmert, “Thermospheric mass density: a review,” Advances in Space Research, vol. 56, no. 5, pp. 773–824, 2015.

K. S. W. Champion and F. A. Marcos, “The triaxial-accelerometer system on atmosphere explorer,” Radio Science, vol. 8, no. 4, pp. 297–303, 1973.

E. Doornbos, Thermospheric Density and Wind Determination from Satellite Dynamics (Ph.D. Dissertation), University of Delft, 2011.

D. A. Vallado and P. Crawford, “SGP4 orbit determination,” in AIAA/AAS Astrodynamics Specialist Conference and Exhibit, Honolulu, Hawaii, August 2008.

C. Shi, W. Li, M. Li, Q. Zhao, and J. Sang, “Calibrating the scale of the NRLMSISE00 model during solar maximum using the two line elements dataset,” Advances in Space Research, vol. 56, no. 1, pp. 1–9, 2015.

S. Krzysztof, “Impact of the atmospheric drag on Starlette, Stella, Ajisai, and Lares Orbits,” Artificial Satellites, vol. 50, no. 1, pp. 1–18, 2015.

C. Arduini, L. Broglio, and U. Ponzi, “Drag balance measurements in the San Marco D/L Mission,” Advances in Space Research, vol. 13, no. 1, pp. 185–200, 1993.

F. Santoni, F. Piergentili, and F. Graziani, “Broglio drag balance for neutral thermosphere density measurement on UNICubeSAT,” Advances in Space Research, vol. 45, no. 5, pp. 651–660, 2010.

L. Felicetti, F. Piergentili, and F. Santoni, “Thermosphere density and wind measurements in the equatorial region using a constellation of drag balance nanospacecraft,” Advances in Space Research, vol. 54, no. 3, pp. 546–553, 2014.

E. Doornbos, H. Klinkrad, and P. Visser, “Use of two-line element data for thermosphere neutral density model calibration,” Advances in Space Research, vol. 41, no. 7, pp. 1115–1122, 2008.

H. S. Jeon, S. Cho, Y. S. Kwak et al., “Mass density of the upper atmosphere derived from Starlette’s precise orbit determination with satellite laser ranging,” Astrophysics and Space Science, vol. 332, no. 2, pp. 341–351, 2011.

Z. Wang and Y. Zhang, “Acquirement of pure gravity orbit using precise formation flying technology,” Acta Astronautica, vol. 82, no. 1, pp. 124–128, 2013.

H. Liu, Z. Wang, and Y. Zhang, “Modeling and analysis of radiometer effect on the inner-satellite,” Chinese Journal of Space Science, vol. 30, pp. 243–249, 2010, (in Chinese).

H. Liu, Z. Wang, and Y. Zhang, “Analysis of radiometer effect on proof mass in purely gravitational orbit,” Applied Mathematics and Mechanics, vol. 33, no. 5, pp. 583–592, 2012.

W. Zhaokui, Z. Zeyang, and Z. Yulin, “Mission overview and latest progress of the Tsinghua scientific satellite,” in 11th IAA Symposium on Small Satellites for Earth Observation, Berlin, Germany, April 2017.

Y. Zhang, W. Zhaokui, and Z. Zeyang, “The basic satellite platform of spherical configuration,” Tech. Rep., China Patent. No. CN106143951A, 2016.

O. Montenbruck and E. Gill, Satellite Orbits: Models, Methods, and Applications, Springer, Berlin, 2000.

E. Doornbos, H. Klinkrad, R. Scharroo, and P. Visser, “Thermosphere density model calibration in the orbit determination and prediction of ERS-2 and Envisat,” in Proceeding of Envisat symposium, H. Lacoste, Ed., Montreux, Switzerland, April 2007.

J. van den IJssel and P. Visser, “Determination of non-gravitational accelerations from GPS satellite-to-satellite tracking of CHAMP,” Advances in Space Research, vol. 36, no. 3, pp. 418–423, 2005.

T. van Helleputte, E. Doornbos, and P. Visser, “CHAMP and GRACE accelerometer calibration by GPS-based orbit determination,” Advances in Space Research, vol. 43, no. 12, pp. 1890–1896, 2009.

P. M. Mehta, A. Walker, C. A. McLaughlin, and J. Koller, “Comparing physical drag coefficients computed using different gas–surface interaction models,” Journal of Spacecraft and Rockets, vol. 51, no. 3, pp. 873–883, 2014.

D. Gu, Y. Liu, B. Yi, J. Cao, and X. Li, “In-flight performance analysis of MEMS GPS receiver and its application to precise orbit determination of APOD-A satellite,” Advances in Space Research, vol. 60, no. 12, pp. 2723–2732, 2017.