Real-Time Thermospheric Density Estimation Via Two-Line-Element Data Assimilation

Inaccurate estimates of the thermospheric density are a major source of error in low Earth orbit prediction. To improve orbit prediction, real-time density estimation is required. In this work, we develop a reduced-order dynamic model for the thermospheric density by computing the main spatial modes of the atmosphere and deriving a linear model for the dynamics. The model is then used to estimate the density using two-line element (TLE) data by simultaneously estimating the reduced-order modes and the orbits and ballistic coefficients of several objects using an unscented Kalman filter. Accurate density estimation using the TLEs of 17 objects is demonstrated and validated against CHAMP and GRACE accelerometer-derived densities. Finally, the use of the model for density forecasting is shown.


Introduction
Accurate knowledge of the thermospheric density is essential for orbit prediction in low Earth orbit and in particular for conjunction assessments. The most accurate models of the thermosphere are physics-based models, such as the Global Ionosphere-Thermosphere Model (GITM) (Ridley et al., 2006) and the Thermosphere-Ionosphere-Electrodynamics General Circulation Model (TIE-GCM) (Qian et al., 2014). These models solve the continuity, momentum and energy equations for a number of neutral and charged components. Their modeling and prediction performance comes, however, at a high computational cost. The models are very high-dimensional, solving Navier-Stokes equations over a discretized spatial grid involving 10 4 -10 6 state variables and 12-20 inputs and internal parameters. In addition, to fully exploit the forecasting potential of physics-based models the schemes employed for data assimilation need to be improved (Sutton, 2018).

1/24
Empirical models (Jacchia, 1970;Hedin, 1987;Picone et al., 2002;Bowman et al., 2008;Bruinsma, 2015), on the other hand, capture the average behaviour of the atmosphere using low-order, parameterized mathematical formulations based on historical observations. A major advantage of empirical models is that they are fast to evaluate, making them ideal for drag and orbit computations. The accuracy of these empirical models is however limited (He et al., 2018), especially during space weather events. Improved densities can be obtained by calibrating empirical density models using satellite data. The current Air Force standard is the High Accuracy Satellite Drag Model (HASDM) (Storz et al., 2005), which is an empirical model that is calibrated using observations of calibration satellites. These satellite observations are used to determine atmospheric model parameters based on their orbit determination solutions. Due to the lack of access to space surveillance observations, publicly available two-line element (TLE) data have been used in the past to estimate the thermospheric density (Picone et al., 2005;Emmert et al., 2006) and calibrate empirical models (Cefola et al., 2004;Yurasov et al., 2005;Doornbos et al., 2008;Chen et al., 2019). Cefola et al. (2004) and Yurasov et al. (2005) calibrated the GOST and NRLMSISE-00 density models using two scaling parameters based on the fitted ballistic coefficient (BC) values. Doornbos et al. (2008) estimated spherical harmonics coefficients to calibrate NRLMSISE-00 model using TLE-derived density estimates and Sang et al. (2011); Chen et al. (2019) adjusted the 187 coefficients of the DTM87 density model directly during orbit determination of multiple objects. On the other hand, Crowley and Pilinski (2017) used TLE data for data assimilation in the physics-based Dragster model using Ensemble Kalman filtering. Except for the Dragster model, the calibrated empirical models have limited forecasting capability which reduces their effectiveness for orbit prediction.
The benefits of this approach are that: 1) the modes contain most energy of the system and are orthogonal; 2) the modes can be estimated in real-time thanks to the dynamic model, which is not possible with a static model; 3) the dynamic model can be used for forecasting (for static models only the fitted coefficients can be extrapolated in time).
In this work, the reduced-order modelling technique for density estimation is further developed and TLE data is used to estimate the thermospheric density. The availability of TLE data for thousands of objects make them attractive for density estimation; however, the use of TLEs is challenging due to the limited accuracy of the orbital data. The density estimation using TLE data is achieved by simultaneously estimating the orbits and BCs of several objects and the reduced-order density state using an unscented Kalman filter. The main contributions of the paper are: 1. Nonlinear space weather inputs are introduced to improve ROM prediction.
2. Two new ROM models based on the NRLMSISE-00 and JB2008 models are developed to extend the maximum altitude to 800 km.
3. Modified equinoctial elements are employed to express the orbit measurements.
4. Thermospheric densities are estimated using reduced-order modeling and TLE data.
5. Accurate density estimation over extended periods of time using a limited amount of TLE data through the use of the dynamic density models is demonstrated.
6. The estimated densities are validated against CHAMP and GRACE accelerometer-derived density data.
The paper is structured as follows. First the development of a dynamic reduced-order density model described. After that, the estimation of the density via TLE data assimilation using an unscented Kalman is discussed. Then the performance of the ROM density prediction and estimation is assessed using simulated and real TLE cases, and finally conclusions are drawn.

Methodology
The neutral density estimation approach consists of two main components: 1) the development of a dynamic reduced-order model (ROM) for the thermosphere and 2) the calibration of the ROM through assimilation of TLE data.

Reduced-order modeling
The main idea of reduced-order modeling is to reduce the dimensionality of the state space while retaining maximum information. In our case, the full state space consists of the neutral mass density values on a dense uniform grid in latitude, local solar time and altitude. The goal is to develop a model for the density evolution over time. First, to make the problem tractable, the state space dimension is reduced using POD. Second, a linear dynamic model is derived by applying DMDc.

Proper orthogonal decomposition
The concept of order reduction using POD is to project the high-dimensional system and its solution onto a set of low-dimensional basis functions or spatial modes, while where Z 1 = U T r X 1 and Z 2 = U T r X 2 are the reduced-order snapshot matrices, and A r = U T r AU r and B r = U T r B are the reduced-order dynamic and input matrices. The above equation is modified such that: where Ξ and Ψ are the augmented operator and data matrices, respectively: We now estimate the dynamic and input matrices by minimizing Z 2 − ΞΨ . The augmented operator matrix is then solved for by computing the pseudoinverse of Ψ: where the + subscript indicates the pseudoinverse. In MATLAB, this can be easily computed using the backslash operator: Calculating A r and B r without first computing A and B is an improvement with respect to previous work (Mehta et al., 2018) and is also numerically more stable because the number of entries that needs to be determined for matrix A r is much smaller than for A. Finally, the discrete-time matrices are converted to continuous-time matrices to enable continuous-time propagation of the ROM modes for estimation. This can be achieved using the following relation (DeCarlo, 1989): where A c is the dynamic matrix and B c is the input matrix in continuous time, and T is the sample time, i.e. the snapshot resolution.

ROM density model development
Density training data In this work, we have developed three different ROM density models using three different atmospheric models to obtain the snapshot matrices, namely the empirical NRLMSISE-00 (Picone et al., 2002) and Jacchia-Bowman 2008 (JB2008) models (Bowman et al., 2008) and the physics-based TIE-GCM model (Qian et al., 2014). We first defined a spatial grid s in local solar time, geographic latitude and altitude and computed the density on this grid for every hour over 12 years (one solar cycle), resulting in over 105,000 snapshots. These snapshots were then used to compute a dynamic ROM model, as described in the previous section, using a reduced order of r = 10. It should be noted that we computed the variation of the densityx by first taking the log base 10 of the density and then subtracting the mean: x = log 10 x − log 10x , where x andx are the density and mean density on the spatial grid. Details on the spatial partitioning and 12 year periods applied for generating the ROM models can be found in Table 1. Note that the JB2008 ROM model was computed over the years 1999-2010 instead of 1997-2008, because no continuous space weather data was available in the year 1998. An improvement with respect to previous work is the development of ROM models that are valid above 450 km altitude, which is the limiting altitude for TIE-GCM. The new ROM models based on NRLMSISE-00 and JB2008 extend up to 700 and 800 km altitude, respectively, see Table 1.
Space weather inputs The space weather inputs u k used in the dynamical model are taken from the inputs required by the original density models, see second column in   Table 2. In addition to these default inputs, we added the next-hour values for key space weather indices to improve the DMDc prediction, see third column in Table 2. Finally, a new innovation in this work is the addition of nonlinear space weather terms, such as the square of an index, e.g. ap 2 , or the multiplication of two different indices, e.g. ap · F 10.7 , see nonlinear inputs in Table 2. The improvement of the DMDc model due to adding nonlinear terms will be discussed in the results section.

Density estimation
The neutral mass density is estimated through the assimilation of two-line element orbital data in the dynamic ROM model. This is achieved by simultaneously estimating the ROM state and the orbit and BC of objects using an unscented Kalman filter.

Two-line element data
The US Air Force Space Command publicly distributes the orbital data of thousands of Earth-orbiting objects in the form of two-line element sets. From this TLE data, the state of an object (position and velocity) at any epoch can be extracted using the SGP4/SDP4 models (Hoots and Roehrich, 1980;Vallado et al., 2006). Hence, the effect of drag can be observed in TLE orbital data if the drag perturbation is strong enough.
A general concern when using TLE data is the accuracy of the orbital data. In the SGP4/SDP4 models only the largest perturbations are included, while higher-order and short-periodic effects, a dynamic atmosphere and orbital maneuvers are not accounted for (Vallado and Cefola, 2012). As a result, there can be large errors in the orbital data (Kelso, 2007). In addition, since 2013, TLEs are fitted to a higher-order orbital solution that includes a future orbit prediction (Hejduk et al., 2013). This generally improve the TLE accuracy at epoch, but may deteriorate the quality if inaccurate future space weather is used for the orbit prediction. To gain understanding about errors in TLE data, we compared the position according to TLE data against GPS data for a Planet Labs satellite at 494 km altitude, see Figure 1. The position error is largest in the along-track direction and varies with a 12-hour period, which is thought to be due to missing tesseral m-daily terms in the SGP4 theory (Herriges, 1988  errors are expected to be limited when computing the state at epochs prior to the TLE epoch (see top plot in Figure 1), which corresponds with the period of tracking data used to generate the TLE. The objects used for density estimation need to be selected carefully. First of all, the objects preferably have a strong drag signal. In particular, the effect of drag should be strong with respect to other non-conservative force effects, else errors in non-conservative force modelling, such as solar radiation pressure, can result in inaccurate density estimates. Second, it is important to have an accurate estimate of the true BC of the object to minimize errors due to inaccuracies in the BC. Therefore, the variation of the BC over time should preferably be very small or else must be modelled accurately (Bowman, 2002). An overview of the objects used in this work is shown in Table 3. The BC values were taken from Bowman et al. (2004), Emmert et al. (2006) and Lu and Hu (2017). Finally, any orbit maneuvers or outliers in the TLE data must be detected and excluded from the data before estimation.

Square-root unscented Kalman filter
To fuse the model and TLE data, we use the unscented Kalman filter (UKF). The UKF was proposed by Julier and Uhlmann (1997) as an extension of the popular Kalman filter (Kalman, 1960) for application to nonlinear systems. Similar to the extended Kalman filter (EKF), the UKF assumes all variables have Gaussian distributions. However, instead of a first-order linearization of the nonlinear dynamics used by the EKF, the UKF uses an unscented transform (UT) to avoid large errors in the true posterior mean and covariance of the variables. Here, the true posterior mean and covariance are computed to the 3rd order by propagating a carefully selected set of sample points, called sigma points, through the true nonlinear dynamics. The UKF is a popular algorithm that is well documented in literature; therefore, we will only present the algorithm and relevant details. More details about the square-root UKF used in this work can be found in Wan and Van Der Merwe (2001).
Let us assume a random variable x ∈ R L with meanx and covariance P x , that is propagated through a nonlinear function f f f such that y = f f f (x). UT uses a set of sigma points to compute the statistics of y. This is achieved by generating a matrix X X X of 2L+1  Bowman et al. (2004), Emmert et al. (2006) and Lu and Hu (2017 Emmert et al. (2006) and Lu and Hu (2017) sigma vectors X i with corresponding weights W i using the following relationships: where λ = α 2 (L + κ) − L is a scaling parameter, α determines the spread of the sigma points aroundx, κ is a secondary scaling parameter, and β is used to incorporate prior knowledge of the distribution of x. Based on the suggested values of the parameters and prior experience, we set the values as α = 1, β = 2, and κ = 3 − L. The above computed sigma vectors are propagated through the nonlinear function: and the mean and covariance for y are approximated using a weighted sample mean and covariance of the posterior sigma points as follows: The UKF extends the UT to recursive estimation through the algorithm given below (Wan and Van Der Merwe, 2001). The state, dynamics, measurements and noise used to estimate the density with the UKF are described in the following.

Algorithm 1 Square-Root Unscented Kalman Filter
Initialize with:x Sigma point calculation and time update: Measurement update: where Q is the process noise covariance and R is the measurement covariance.

9/24
State The state x that is estimated in the UKF consists of the osculating orbital states (expressed in modified equinoctial elements) and the BCs of the objects plus the reduced-order density state z: x = p 1 , f 1 , g 1 , h 1 , k 1 , L 1 , BC 1 , ... , p n , f n , g n , h n , k n , L n , BC n , z T T where n is the number of objects and the modified equinoctial elements (MEE) are defined as (Walker et al., 1985): f = e cos (ω + Ω), g = e sin (ω + Ω), h = tan (i/2) cos Ω, k = tan (i/2) sin Ω, L = Ω + ω + ν.
where a, e, i, Ω, ω and ν are the classical Keplerian orbital elements. The MEE are nonsingular and tend to behave less nonlinear than the Cartesian coordinates (used in previous work) which benefits the Kalman filter estimation. To initialize the state, we use the objects' orbital states according to TLE data and take the BC values from Table 3. The ROM state is initialized using densities from the JB2008 model.

Dynamic model
The dynamic model f f f (x, t) for evolving the state x consists of propagating the orbital states using orbital dynamics and evolving the ROM state using the continuous-time DMDc model: x where [x, y, z] and [v x , v y , v z ] are the inertial position and velocity, respectively, that are used for orbit propagation and BC = C D A m is the ballistic coefficient. The ROM state z is used to compute the atmospheric density by converting z to the full space (see Eq. 5) and interpolating the density grid. After propagation, the Cartesian state is converted back to MEE, which are used in the UKF.
Orbital dynamics The orbital dynamics used in this work considers: • Geopotential acceleration computed using the EGM2008 model, up to degree and order 20 for the harmonics; • Atmospheric drag considering a rotating atmosphere for computing the velocity relative to the atmosphere. The atmospheric density is computed using the ROM density model.
The orbit propagation is carried out in the inertial J2000 reference frame using Cartesian position and velocity while the geopotential and drag accelerations are computed in the Earth-fixed ITRF93 frame. NASA's SPICE toolbox is used for reference frame and time transformations (ITRF93 and J2000 reference frames and leap-seconds kernel). Perturbations due to solar radiation pressure, gravitational attraction by the Sun and Moon, and higher-order Earth harmonics are not included, because their effect on the considered orbits during density estimation was found to be negligible compared to other modeling and measurement errors.
Measurements The measurements used for estimation are the osculating orbital states extracted from TLE data. At one hour intervals the osculating state of each object is computed using the nearest newer TLE by propagating the TLE backward to the measurement epoch using SGP4. These states are then converted to MEE and used as measurements. The 17 objects used in this work for density estimation are shown in Table 3. Note that for calibrating the ROM-TIEGCM model only 11 objects are used, because 6 of the 17 objects have their perigee above the ROM-TIEGCM maximum altitude of 450 km.
Measurement and process noise The measurements noise R was determined empirically. The variance for the measurements of p, f and g was scaled by the orbit's eccentricity e, because the errors were found to increase with increasing eccentricity: [R p , R f , R g , R h , R k , R L ] = [c 1 · 10 −8 , c 2 · 10 −10 , c 2 · 10 −10 , 10 −9 , 10 −9 , 10 −8 ] where c 1 = 1.5 · max(4e, 0.0023) and c 2 = 3 · max(e/0.004, 1). The process noise variance Q for the state and BC was set to: 2·10 −14 , 10 −14 , 10 −14 , 10 −12 , 10 −16 ] (37) The process noise for the ROM state Q z was computed using the 1-hour ROM prediction error on the training data: As a result of this approach, the Kalman filter will give more confidence to the model prediction with respect to measurements when the ROM prediction on the training data is more accurate. This approach for determining Q z was also found to result in good estimates of the uncertainty in the estimated density. Finally, the initial covariance for the state was set to: [P p , P f , P g , P h , P k , P L , P BC , P z1 , P zn ] = [R p , R f , R g , R h , R k , R L , (0.005 · BC) 2 , 20, 5] (39) where P z1 refers to the covariance for the first reduced-order mode z 1 and P zn to the covariance for all other modes. An overview of the reduced-order model density estimation technique is shown below.
Algorithm 2 ROM density estimation Reduced-order modeling 1. Generate density training data X (hourly density on grid) using physics-based or empirical density model 2. Select reduced order r 3. Compute reduced-order model using a SVD of the snapshots X (Eqs. 1-4) 4. Compute DMDc for reduced-order training data (Eqs. 7-13) Density estimation 5. Download TLE data and estimate true BC 6. Select objects with accurate TLEs (check self-consistency) and stable BC (not maneuvering) 7. Generate measurements (orbital states in MEE) every hour from TLEs using SGP4 using nearest newer TLE 8. Estimate ROM modes z using unscented Kalman filter (Algorithm 1) by simultaneously estimating the modes and the state and BC of objects 11/24

Results
In this section, the performance of ROM model forecasting and density estimation using TLE data is assessed.

ROM density prediction
The performance of the dynamic ROM models is tested by comparing density forecasts with training data. Using the three different ROM density models the density was predicted for 5 days during quiet space weather conditions and during a geomagnetic storm in 2002. The resulting density forecast errors (the root mean square (RMS) percentage error on the three-dimensional spatial grid) and space weather conditions are shown in Figure 2. The predictions using the ROM model based on JB2008 is most accurate. This good performance can be explained by the superior space weather proxies used by the ROM-JB2008 model. Moreover, the figure and the average 1-hour prediction error in Table 4 show that the addition of nonlinear space weather terms improves the prediction accuracy for all models. The nonlinear terms especially improve the prediction during a geomagnetic storm. The ROM-JB2008 provides the best predictions with respect to training data; however, this does not necessarily mean it will perform better in estimating true densities.

Simulated TLE test case
To assess whether accurate density estimation using a ROM model and TLE data is feasible, we first tested the technique using simulated TLE data. In this scenario, the 'true' orbits and densities are first computed by propagating several objects using the ROM-TIEGCM density model. The initial orbital parameters for the true orbits are  Table 5. Based on these 'true' orbits, TLE data is simulated using realistic uncertainties. This simulated TLE data is then used to calibrate the ROM-TIEGCM model. Table 6 shows the 1-σ errors used for simulated TLE measurements. These errors were established empirically based on TLE data of near-circular orbits around 400 km altitude in the years 2017 and 2018. These errors convert to RMS position and velocity errors of 0.88 km and 1.0 m/s, respectively. The initial guesses for estimating the orbits, BCs and ROM-state also include realistic errors. Figure 3 shows the errors in the estimated BCs and densities along the orbits during the estimation period. The errors in density and BC remain below 2% after 12 days estimation. In Figure 4 the true and estimated values and errors of the first four ROM modes are shown. All modes converge to their true values, while a small bias in the first mode remains. The convergence of the modes is also correctly displayed by the 3σ error bounds. This result demonstrates that accurate density estimation using TLE data is feasible and that both the BC and density are observable. Nevertheless, in reality, less accurate density estimates can be expected, because new TLE measurements are not available every hour and errors in TLE data are not random nor Gaussian.

Real TLE
In the following, the density is estimated using real TLE data and compared with CHAMP and GRACE accelerometer-derived density data. Figure 5 shows the orbit-averaged estimated density along CHAMP's orbit as well as the density according to CHAMP data and the NRLMSISE-00 and JB2008 density models during August 2002 (the first month for which we had both CHAMP and GRACE data). All three ROM models perform very well. Especially, the ROM-NRLMSISE and ROM-JB2008 models are very well able to estimate density variations due to changes in solar activity. The RMS error in daily-averaged CHAMP density is only about 6-9% for all models, see Table 7. The wiggles in the estimated orbit-averaged density (particularly visible for   ROM-TIEGCM and ROM-NRLMSISE) have a period of 12 hours, which suggests that these are due to errors in the TLEs due to missing m-dailies. Similar 12-hour variations can be found in the estimated BCs (not shown here). Further tuning of the measurement and process noise can reduce the amplitude of these variations due to TLE errors.
Overall, the ROM model based on JB2008 performs best with a RMS error in orbit-averaged density of only 7.3% and 12.1% along CHAMP's and GRACE-A's orbit, respectively, see Table 7. This shows the high accuracy and temporal resolution that can be achieved by the ROM approach using only TLE data. The error in GRACE-A density is possibly higher because less objects around GRACE's 480 km altitude were used than around CHAMP's altitude of 400 km. Only 17 objects were used for calibrating the ROM-JB2008 and ROM-NRLMSISE models and only 11 for the ROM-TIEGCM model. This is significantly lower than the 36 and 48 objects used by Doornbos et al. (2008) and Shi et al. (2015), respectively, but similar to the 16 objects used by Yurasov et al. (2005). It should also be noted that the errors presented in this paper are with respect to accelerometer-derived density data and not with respect to TLE-derived density data as in some other works (Doornbos et al., 2008;Shi et al., 2015).
A close-up of the density along CHAMP's orbit on August 14, 2002 is shown in Figure 6. In this time window, the ROM-estimated densities are very close to the true density (both the mean and variation are estimated well). It is probably not feasible to exactly match the true density using TLE data due to the highly-dynamic character of the atmosphere and low-temporal resolution of TLE data. Besides the good match, one can see differences between the different ROM models in the density variation over one orbit. These differences stem from the underlying base models which result in different spatial modes for each model. The way that the ROM models mimic their base model can also be seen in Figure 7 that shows maps of the modelled and estimated density on August 8, 2002 at 450 km altitude. For example, the simple density distribution in the JB2008 model is also visible in the ROM-JB2008 density, whereas the NRLMSISE-00 and TIEGCM based ROM models show more complex density distributions. On the other hand, independent of the base model, the ROM-estimated densities have a similar magnitude, which indicates successful calibration. Figures 8 and 9 show the uncertainty in the estimated density. The plots show that the error in the estimated density is smaller for lower altitude and inside the diurnal bulge. This indicates that the density estimation is more accurate when the drag signal is stronger. In addition, Figure 9 shows that the uncertainty grows little at altitudes where measurements are available. The 1-σ error varies between 3 and 11%, while Mehta and Linares (2018b) found a 1-σ error of 5% along CHAMP's orbit and up to 25% error at higher altitudes after assimilating CHAMP density data in a TIEGCM-based ROM model. This indicates that the use of data from multiple objects improves the global density estimates.

Full years 2003 and 2007
The neutral mass density was estimated using the ROM-JB2008 model over the entire years 2003 (high solar activity) and 2007 (low solar activity). The difference in the estimated and true densities along CHAMP's and GRACE-A's orbits are shown in Table 8 and Figure 10. Both the ROM estimation and JB2008 model perform very well in 2003, whereas the ROM estimates are much more accurate than the JB2008 and NRLMSISE-00 models in 2007. The accuracies reported here can be further improved    Figure 10. Error in orbit-averaged density with respect to CHAMP data for ROM-JB2008 estimated density, and JB2008 and NRLMSISE-00 models for the year 2007.
by using more accurate orbital data and by improving the spatial and temporal coverage of the measurements. Figure 11 shows the density along CHAMP's orbit estimated using the ROM-JB2008 model during a major geomagnetic storm in 2003. Both the ROM and empirical models provide good density estimates during the storm. However, the empirical models overestimate the density after storm on day 151 (which results in large relative density errors), whereas the ROM estimates are much more accurate. This example shows that the linear ROM model is able to deal with space weather events even though these events are strongly nonlinear.

Including density data
For some periods of time, one may have access to highly-accurate density data, such as CHAMP, GRACE or GOCE accelerometer-derived densities. This data can be included in the data assimilation to improve the global density estimates. Figure 12 shows the estimated orbit-averaged density along GRACE-A's orbit after assimilating CHAMP density data together with TLE data (here CHAMP was at 360 km and GRACE-A at 480 km altitude). This period coincides with the period of reduced estimation accuracy  shown in Figure 10 around day 50. One density measurement is included each hour at the same time as the TLE orbit measurements. The inclusion of the CHAMP densities significantly improves the GRACE density estimates; the error in daily-averaged density reduced from 16.4% to 11.6%. Therefore, global density estimates can be improved by including accurate density data. This is especially useful in case the drag signal is weak and consequently the density is difficult to observe from orbital data such as during periods of low solar activity.

Density forecast
The density along CHAMP's orbit was predicted for 11 days using the ROM-NRLMSISE model, see Figure 13. The initial ROM state used for the prediction was obtained after 30 days calibration in August 2002, see Figure 5. The ROM does very well in predicting the density, even during two geomagnetic storms. This example shows that ROM models are able to accurately forecast the future density if the initial state of the thermosphere and the future space weather are accurately known. In future work, prediction of the future space weather will also be considered.

Conclusions
In this paper we have presented the development of a dynamic ROM model for the thermosphere and estimated its state using TLE data. Three different models based on TIE-GCM, NRLMSISE-00 and JB2008 were generated with an upper altitude up to 800 km. The prediction performance of the models was improved by including nonlinear space weather terms as inputs for the DMDc. In addition, the estimation using an UKF was improved by expressing the measurements in modified equinoctial elements and by 19/24 calculating the process noise for the ROM model based on training data performance. Densities were estimated using TLE data and compared with CHAMP and GRACE accelerometer-derived density data. The results showed that the dynamic model enables accurate density estimation using the TLEs of only a small number of satellites. Improved global density estimates can be obtained by including accurate density measurements. Finally, the ROM was shown to be able to provide accurate density forecasts. Future work will focus on further improving the ROM dynamic models to deal with nonlinearities by improving the choice of space weather inputs and by applying Koopman operator theory. In addition, we can estimate the global thermospheric neutral density using historic TLE data from the 1960's up to present time to generate a density database.