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.
(a)
(b)
(c)
(d)
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)
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.
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.
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.
(a)
(b)
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)
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)
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.
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−9 m/s2); 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 m2/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.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
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)