EFFECTS OF DIFFERENT DRAG FREE SYSTEM ACCELERATION NOISE LEVELS FOR FUTURE SATELLITE GEODESY MISSIONS By SEONG HYEON HONG A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2017
201 7 Seong Hyeon Hong
To my m other in h eaven and my lovely wife
4 ACKNOWLEDGMENTS First of all, I owe many thanks to Dr. John W. Conklin for being my advisor. I know this research we worked on for past few years was very adventurous and require d a lot of patience. In fact I could not have completed this work without h is teaching, guidance and patience. I also thank all my other committee members; Dr. Norman Fitz Coy, Dr. Prabir Barooah and Dr. Guido Mueller for all of their time and considerations. I thank Stephen Apple from UF Precision Torsion Pendulum team for providi ng noise data and help plotting them. In addition, I thank all other colleagues from Precisions Space Systems Laboratory for all their support and friendship. I thank David Wiese from JPL for his advises through emails and teleconferences and his comments during the conferences. I would also like to acknowledge Thomas Gruber from Technical University of Munich and Elisa Fagiolini from GFZ for providing additional information about their gravity field models. This research was funded by Florida Space Instit ute and I am very grateful of their support. Numerous thanks to my parents for all their love and support. I am really proud to be your son. Special thanks to my lovely wife for all her care beside s me. We did this together. I thank all my friends all ove r the world for their encouragements and prayers. Lastly and most importantly, I thank GOD for all his guidance, support and mercy he has shown me throughout my life. I came this far, all because of your allowance and the grace is all yours.
5 TABLE OF CON TENTS page ACKNOWLEDGMENTS ................................ ................................ ................................ .. 4 LIST OF TABLES ................................ ................................ ................................ ............ 7 LIST OF FIGURES ................................ ................................ ................................ .......... 8 LIST OF ABBREV IATIONS ................................ ................................ ........................... 10 ABSTRACT ................................ ................................ ................................ ................... 12 C HAPTER 1 INTRODUCTION ................................ ................................ ................................ .... 14 History of Satellite Geodesy ................................ ................................ .................... 14 Gravity Measurements before Satellite Geodesy ................................ ............. 14 Satellite Geodesy ................................ ................................ ............................. 15 Early Satellite Geodesy Missions ................................ ................................ ..... 15 The CHAMP Mission ................................ ................................ ........................ 16 The GRACE Mission ................................ ................................ ........................ 17 The GOCE Mission ................................ ................................ ........................... 19 The GRACE Follow On Mission ................................ ................................ ....... 19 Applications and Resolution Requirements of Gravity Measurements ................... 20 Applications of Gravity Measurements ................................ ............................. 20 Resolution Requirements ................................ ................................ ................. 21 Goals for Future Geodesy Missions ................................ ................................ ........ 24 Motivation ................................ ................................ ................................ ............... 25 Dissertation Contr ibution and Outline ................................ ................................ ..... 27 2 BASIC THEORY IN SATELLITE GEODESY AND PREVIOUS STUDIES .............. 29 Representation of Gravity (or gravitational) Potential ................................ ............. 29 Analysis of Gravity Field ................................ ................................ ......................... 31 Sa tellite Formations ................................ ................................ ................................ 34 Singe Pair LL SST ................................ ................................ ............................ 35 Double Pair LL SST ................................ ................................ .......................... 36 Optimal Configurations ................................ ................................ ..................... 37 Drag f ree Systems ................................ ................................ ................................ .. 38 Major Source of Errors ................................ ................................ ............................ 40 3 MISSION ARCHITECTURE ................................ ................................ .................... 42 Orbit Selection ................................ ................................ ................................ ........ 42 Ground Track Patterns ................................ ................................ ............................ 45
6 Technology and Associated Noise ................................ ................................ .......... 48 Position and Range Noise Models ................................ ................................ .... 48 Acceleration Noise Models ................................ ................................ ............... 49 Generating Colored Noise ................................ ................................ ...................... 51 Summ ary of the Simulation Design ................................ ................................ ......... 52 4 ESTIMATION AND SIMULATION PROCEDURE ................................ ................... 55 Gravity Model Selection ................................ ................................ .......................... 55 Crea ting Measurement Data ................................ ................................ ................... 58 Estimation Method ................................ ................................ ................................ .. 59 Simulation Procedure ................................ ................................ .............................. 66 5 RESULT AND ANALYSIS ................................ ................................ ....................... 68 Repres entation of Solutions ................................ ................................ .................... 68 Geoid Degree Difference Error ................................ ................................ ............... 69 Equivalent Water Height ................................ ................................ ......................... 73 Disc ussion ................................ ................................ ................................ .............. 82 6 IMPROVEMENTS IN ESTIMATION ................................ ................................ ....... 88 State and Measurement Uncertainty ................................ ................................ ...... 88 ARMA Filter ................................ ................................ ................................ ............ 89 Accounting for A cceleration Noise in the Estimation Routine ................................ 93 Applying both the ARMA Approximation and the ANC ................................ ........... 94 Result ................................ ................................ ................................ ...................... 95 Applying the ARMA Approximation and the ANC for Higher Acceleration Noise Level Mission ................................ ................................ ................................ ..... 101 Discussion ................................ ................................ ................................ ............ 105 7 TOOLBOX VALIDATION USING THE UF PRECISION TORSION PENDULUM 108 UF Precision Torsion Pendulum ................................ ................................ ........... 108 Accumulated Acceleration Noise ................................ ................................ .......... 110 Result ................................ ................................ ................................ .................... 113 Discussion ................................ ................................ ................................ ............ 117 8 CONCLUSIONS AND RECOMMENDATIONS ................................ ..................... 119 Conclusions ................................ ................................ ................................ .......... 119 Recommendations for Future Research ................................ ............................... 120 LIST OF REFERENCES ................................ ................................ ............................. 122 BIOGRAPHICAL SKETCH ................................ ................................ .......................... 127
7 LIST OF TABLES Table page 1 1 Applications of g ravity m easurements. ................................ ............................... 20 2 1 Technology and associated residual acceleration noise. ................................ .... 41 3 1 Summary of orbit parameters used for both the polar pair and inclined pair of satellites. ................................ ................................ ................................ ............ 47 3 2 Measurement technologies used in this analysis. ................................ ............... 48 3 3 LRI and Acceleration noise models used in this analysis. ................................ .. 49 3 4 Summary of simulation design with respect to error sources and launched missions. ................................ ................................ ................................ ............ 54 4 1 Stat ic and dynamic gravity models used in this analysis. ................................ ... 56 5 1 Spatial RMS errors for all cases in units of cm EW H ................................ ......... 80 5 2 Differences between current and Wiese simulations. ................................ ......... 85 6 1 Spatial RMS error for Case A2 in units of cm of EW H with various MNC matrices. ................................ ................................ ................................ ........... 100 6 2 Spatial RMS error for Case A1 in units of cm of EW H with various MNC matrices. ................................ ................................ ................................ ........... 105 6 3 Spatial RMS error for Cases A1 and A2 in units of cm of EW H ....................... 107 7 1 Spatial RMS error for pendulum cases and Case A2 in units of cm of EW H ... 117
8 LIST OF FIGURES Figure page 1 1 Time and spatial resolution requirements. ................................ .......................... 22 1 2 V ariou s static gravity solutions in geoid accuracy. ................................ .............. 23 2 1 Gravitational potential and spherical coordinate system ................................ .... 29 2 2 Reference ellipsoid, geoid and topography. ................................ ....................... 31 2 3 High low and low low satellite to satellite tracking. ................................ ............. 34 2 4 Single pair LL SST configurations. ................................ ................................ ..... 36 2 5 Double pair LL SST configurations. ................................ ................................ .... 37 2 6 Drag free system. ................................ ................................ ............................... 39 3 1 Separation distance between the lead ing and trailing satellites as a function of time. ................................ ................................ ................................ ................ 44 3 2 Ground tracks of polar pair and inclined pairs of satellites. ................................ 46 3 3 The noise model and numerically generated noise. ................................ ........... 53 4 1 Block diagram of the simulation procedure. ................................ ........................ 67 5 1 Errors and Gaussian true models. ................................ ................................ ................................ ........ 70 5 2 Errors and Gaussian smoothed (GS) errors expressed as GDD true models. ................................ ................................ ................................ ........ 71 5 3 and the figure include s hydrology and ice mass variations. ................................ 74 5 4 EWH error versus latitude plots for Case s A1 and A2. ................................ ....... 75 5 5 and the figure includes hydrology and ice mass variations. ................................ 78 5 6 EWH error versus latitude plots for Case s B 1 and B 2. ................................ ....... 79 6 1 ASD of the LRI generated noise, the ARMA filtered noise, and the noise model. ................................ ................................ ................................ ................. 93
9 6 2 ASD of the GOCE level drag free acceleration generated noise, the ARMA filtered noise, and the noise model. ................................ ................................ .... 95 6 3 Errors of the ARMA and the ANC applied single pair solutions for Case A2 expressed as GDD. ................................ ................................ ............................ 97 6 4 Estimated single pair solutions for Case A2 in units of cm of EWH with various MNC matrices. ................................ ................................ ....................... 99 6 5 GS applied single pair solutions for Case A2 in units of cm of EWH with various MNC matrices. ................................ ................................ ..................... 100 6 6 Errors of the ARMA and the ANC applied single pair solutions for Case A1 expressed as GDD. ................................ ................................ .......................... 102 6 7 Estimated single pair solutions for Case A1 in units of cm of EWH with various MNC matrices. ................................ ................................ ..................... 104 7 1 CAD model and actual apparatus of a torsion pendulum taken from Physics Space Laboratory at University of Florida in 2015. ................................ ........... 109 7 2 Inertial members of a torsion pendulum taken from Physics Space Laboratory at University of Florida in 2017. ................................ ...................... 110 7 3 ASD of raw acceleration noise data. ................................ ................................ 111 7 4 ASD of extracted acceleration noise data. ................................ ........................ 112 7 5 ASD of actual pendulum noise, model based simulated pendulum noise, and GOCE level simulated acceleration noise. ................................ ....................... 113 7 6 Errors of the actual and the simulated UF pendulum noise applied single pair solutions expressed as GDD. ................................ ................................ ........... 114 7 7 Estimated single pair solutions in units of cm of EWH applying the pendulum noise. ................................ ................................ ................................ ................ 115
10 LIST OF ABBREVIATIONS ACS A uto covariance sequence ANC A cceleration noise covariance AOD A tmosphere and O cean D e aliasing AR Auto regressive ARMA A uto regressive moving average ASD Amplitude spectral density BLUE B est L inear U nbiased E stimator CHAMP Challenging Mini Satellite Payload CSR Center for Space Research DFAC Drag free Attitude Control DOF Degree of freedom DORIS Doppler Orbitography and Radiopositioning Integrated by Satellite ECMWF European Center for Medium range Weather Forecast EGG Electrostatic Gravity Gradiometer ESA European Space Agency EW H E quivalent w ater h eight FFT fast Fourier transform FTG French Tidal Group GDD Geoid degree difference GFZ Geo Forschungs Zentrum GNSS Gl obal Navigation Satellite System GOCE Gravity field and steady state Ocean Circulation Explorer Mission GPS Global Positioning System GRACE Gravity Recovery and Climate Experiment
11 GRACE FO GRACE Follow On GRS Gravitational reference sensor GS Gaussian Smoothing HI Hydrology and Ice HL SST High low satellite to satellite tracking IFFT Inverse fast Fourier transform JPL Jet Propulsion Laboratory KBR K Band microwave ranging LAGEOS Laser Geodynamic Satellites LISA Laser Interferometer Space Antenna LL SST Low low satellite to satellite tracking LRI Laser r anging i nterferometer MA Moving average MNC Measurement noise covariance OMCT Ocean Model for Circulation and Tides PCCG P reconditioned conjugate gradient PPHA Pacanowski, Ponte, Hirose and Ali PSD Power spectral density RAAN R ight ascension of the ascending node RK4 Runga Kutta 4 th order RMS Root mean square SLR Satellite Laser Ranging UF University of Florida
12 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy EFFECTS OF DIFFERENT DRAG FREE SYSTEM ACCELERATION NOISE LEVELS FOR FUTURE SATELLITE GEODESY MISSIONS By Seong Hyeon Hong May 2017 Chair: John W. Conklin Major: Mechanical Engineering This research evaluates the impact of residual acceleration noise on the varying gravity field for future low low satellite to satellite tracki ng missions. The Gravity Recovery And Climate Experiment (GRACE) has provided monthly average gravity field solutions in spherical harmonic coefficients for more than a decade. It provides information about land and ocean mass variations with a spatial re solution of 3 5 0 km and with an accuracy within 2 cm throughout the entire Earth. GRACE Follow on is expected to be launched in 2017 to advance the work of GRACE and to test a new laser ranging interferometer, which measures the range between the two satell ites with higher precision than the K Band ranging system used in GRACE. Moreover, there have been simulation studies that show, an additional pair of satellites in an inclined orbit increases the sampling frequency and thus reduces temporal aliasing error s. Given the fact that future missions will likely continue to use the low low satellite to satellite tracking formation with laser ranging interferometry, it is expected that the residual acceleration noise caused by non gravitational forces and aliasing errors will
13 become limiting factors for the time varying gravity field solution. We evaluate three different levels of residual acceleration noise based on demonstrated drag free systems to find an optimal drag free performance target for upcoming geodesy missions. We analyze both a single collinear polar pair and the optimal double collinear pair of drag free satellites and assume the use of a laser ranging interferometer, performing estimation with two different sets of de aliasing models to represent di fferent levels of aliasing errors. A partitioned best linear unbiased estimator is used to compute the solutions in terms of spherical harmonics. The estimation routine was improved with precise modeling of the measurement noise covariance matrix and the e ntire simulation toolbox was tested by using the actual residual acceleration noise data from the UF precision torsion pendulum. It was found that the suitable residual acceleration noise is around 2 10 12 ms 2 Hz 1/2 Decreasing the acceleration noise below this level did not result in more accurate gravity field solutions for the chosen mission architecture.
14 CHAPTER 1 INTRODUCTION History of Satellite Geodesy Gravity Measurements before Satellite Geodesy Geodesy is a field of study that examines the measurements, which determine the physical quantities of the Earth. Gravity field (or potential) is one of the key measurements in geodesy, which quantify the mass distribution and hence the approximate shape o f the Earth. Before the satellite geodesy era, the gravity field was measured locally by mechanical devices such as pendulums, torsion balances, static springs and etc. [ Nerem et al. 1995]. Each of these devices works in different ways, and as a consequen ce they have different equations of motion. Classical pendulums for instance can measure the gravity by direct measurements of the period of oscillation with known mass and length of the pendulum. From the equation of motion, the natural frequency of the p endulum can be calculated and this can be converted to the local gravitational strength. The method is similar for torsion balances and static springs. If the characteristics of mechanical devices are known and the angular or translational movements are me asured, signals produced by these devices can be converted to local gravity measurements. These measurements have been improved by many years of work to increase the accuracy. However, each measurement can only provide information from the specific point w here the measurement was taken. Consequently with local measurements, it becomes almost impossible to measure the entire gravity field of the Earth since it requires numerous devices and efforts. Another complication is that each device ha s to be calibrate d absolutely with high accuracy so that it can be compared with results from different measurement devices. Furthermore, some places
15 are prohibited by environmental conditions or political boundaries. New technology was desired to measure the gravity field of the entire Earth. Satellite Geodesy Satellite geodesy is a branch of geodesy that uses satellites orbiting around the E arth, mainly determining the gravit y potential with the measurements taken from the satellites. These measurements include position, velocity and acceleration of the satellites, which can be used to estimate the gravity field of the Earth using the equation of motion of the satellites. Since the Earth does not have a uniform mass, it exerts different gravitational forces on the satelli tes at different positions and effects their flight paths. Therefore, by measuring and observing the path of satellites, one can recover the gravitational potential of the Earth. Satellite geodesy became the best way to measure the global gravity field of the Earth since satellites do not have any boundaries and can observe the entire Earth within short period of time. Early Satellite Geodesy Missions Satellite geodesy era began with the first satellite, Sputnik in 1957. It was limited by a discontinuous tr acking system using cameras, radio Doppler and radio interferometry. Satellite Laser Ranging (SLR) was introduced in mid 1960s and it was demonstrated in Laser Geodynamic Satellites (LAGEOS) launched in 1976 [ Nerem et al. 1995]. The accuracy of SLR with t he LAGEOS mission was about 2 cm, which was a big improvement from camera tracking system which had an accuracy of about 10 m [ Nerem et al. 1995]. However, there were still limitations in satellite geodesy caused by the discontinuous tracking of SLR due t o high cost of ground stations and weather uncertainty. In 1990s, continuous tracking by the French Doppler Orbitography and Radiopositioning Integrated by Satellite (DORIS) system was developed [ Nerem et al.
16 1995]. DORIS measures the continuous range rat e between the satellites and the ground stations by emitting radio signals from 50 ~ 60 ground stations, and satellites receiving these signals and measuring the Doppler shifts. The accuracy of range rate measurements were about 0.3 to 0.5 mm/s and these m easurements were used to estimate the orbits of satellites with accuracy of 2.5 cm [ Seeber 2003]. Around that time, the U.S. developed another continuous tracking system called the Global Positioning System (GPS), which was first demonstrated by TOPEX/Pos eidon in 1992 [ Nerem et al. 1995]. GPS tracking of a low Earth orbiting satellite for the purpose of Earth geodesy measurements is known as high low satellite to satellite tracking (HL SST) [ Kim 2000]. GPS can achieve the orbit determination accuracy of as low as 1 cm [ Kang et al ., 2006]. This class of measurements provide continuous position and velocity of satellites with respect to the Earth. GPS along with DORIS initiated an escalation of accurate gravity field measurements. Nevertheless, these techno logies were not sufficient to measure high spatial and temporal frequency gravity field signals with high accuracy. The CHAMP Mission Challenging Mini Satellite Payload ( CHAMP ) (2000 ~ 2010) mission was launched in a near polar orbit and altitude of 450 km This was a pre GRACE technology and science mission, launched 2 years before the actual GRACE mission which is explained in the following paragraph CHAMP consisted of a GPS receiver for continuous tracking and an electrostatic accelerometer to measure in orbit non gravitational forces [ Reigber et al. 2003]. Different types of non gravitational forces are explained later in this dissertation. CHAMP was the first satellite geodesy mission that enabled 550 km spatial resolution gravitational field solutio ns called EIGEN 2. These
17 CHAMP only solutions were at least an order of magnitude improvement compared to previous satellite only solutions [ Reigber et al. 2003]. Also, before the CHAMP mission, all the gravity solutions acquired by satellite missions wer e mostly used to build static models, which means that the solutions represent average gravity potential of the Earth which do not change with time. On the contrary, the gravity potential of the Earth varies with the time due to the mass variations of the Earth. Dynamic gravity solutions represent the variations of the gravity potential with respect to the static potential. Static and dynamic solutions are explained more in the following section of this chapter Combination of the data provided by CHAMP and SLR from other five geodetic satellites, Lageos 1 and 2, Stella, Starlette, and Ajisai produced the annual and semiannual dynamic gravity field solutions with high spatial resolution [ Moore et al. 2005]. The G RACE Mission Gravity Recovery and Climate Exp eriment (GRACE) (2002 ~ present) started a new era of gravity field measurement. GRACE is a mission committed to measuring the E arth once per month This mission consists of t wo identical satellites flying at the initial altitude of around 500 km (currently around 400 km) in nearly circular polar orbits one trailing the other by about 220 km (170 km ~ 270 km). Each satellite is an identical design from the CHAMP mission but GRA CE has the updated electrostatic accelerometer s and an additional K Band microwave ranging (KBR) system between the two. This measures the range rate between the satellites [ Tapley et al. 2004 ]. This type of tracking is known as low low satellite to satellite tracking (LL SST) which is explained more in Chapter 2 The KBR system not only provides additional
18 measurements from another satellite but also keep s track of range between the pair which plays a n important role in increasing the time and spatial resolution and accuracy of gravity field models [ Kim 2000]. The KBR system provides range measurements with an accuracy of less than 10 m which is more accurate by about three orders of magnitude than the range measurement acquired by G PS [ Tapley et al. 2004 ]. GRACE solutions are produced monthly with spatial resolution on the surface of the Earth of ~ 35 0 km. With advancements in estimation techniques the temporal resolution has decreased to as low as one day although this requires pr ior information for a hydrological model [ Kurtenbach et al. 2009 ; Sakumura et al. 2016 ]. Moreover, the accuracy of the monthly solutions were improved at least by 20% over the ocean through numerical modeling of the non conservative acceleration and atti tude errors [ Chen et al. 2016]. On the other hand, a new calibration method for the GRACE accelerometer data that considers temperature effects on the measurements, showed an improvement in estimating the C20 spherical harmonic coefficients [ Klinger and M ayer Grr 2016]. Apart from these previous research efforts, the most important lesson from GRACE was that accelerometer noise was conjectured to be the dominating uncertainty for measurement frequencies between ~0.2 mHz and 2~3 mHz, and at frequencies ab ove 2~3 mHz, the noise of the microwave ranging system was speculated to be the limiting factor for the gravity field solution [ Sheard et al ., 2012]. Therefore, a Laser Ranging Interferometer (LRI) was developed, which is expected to reduce the satellite to satellite ranging noise by at least factor of 5 compared to the noise of the GRACE microwave ranging system [ Sheard et al ., 2012]. The LRI is going
19 to be tested by GRACE Follow On (GRACE FO), which is scheduled to be launched in 2017 to continue the work of GRACE [ Flechtner et al ., 2014]. The G OCE Mission After GRACE, t he Gravity field and steady state Ocean Circulation Explorer (GOCE) (2009 ~ 2013) was launched [ Floberghazen et al ., 2011]. The m ission objective was to recover static gravity field and reach geoid accuracy of 1~2 cm with spherical harmonics of degree 200 which is equivalent to spatial resolution of 100 km. In order to reach a high sp atial resolution solution, GOCE was flying at the low altitude of 279 km and inclination of 96.7 in near circular orbit [ Floberghazen et al ., 2011]. G OCE was equipped with an Electrostatic Gravity Gradiometer (EGG) and Drag free Attitude Control (DFAC) sy stem [ Canuto and Massotti 2009]. These two comprised a drag free system (system that provides satellites to fly in a purely gravitational orbit ), which will be explained in Chapter 2 The EGG reduced the acceleration noise down to 2 10 12 ms 2 Hz 1/2 f or frequency of 5 mHz ~ 0.1 Hz [ Touboul et al 2012]. However, the actual residual acceleration noise of t his drag free system was higher due to the actuator performance Overall, the performance of GOCE drag free system was better than that of the electrost atic accelerometer used in GRACE. The G RACE Follow On Mission The GRACE Follow On (GRACE FO) mission, implemented by a US German partnership will continue the work of GRACE with minimum gap in the data stream for 5 more years and it is s cheduled to be la unched in 2017 [ Fle ch t n er et al ., 2014]. GRACE FO will be launched in to a circular orbit with 490 km altitude and 89 inclination angle similar to GRACE. There are two main objectives for GRACE FO. One is to continue
20 the measurements and observations of th e current GRACE mission and another is to demonstrate the performance of a Laser Ranging Interferometer (LRI) which is expected to improve low low satellite to satellite tracking (LL SST) measurement, and increasing the accuracy of gravity field estimatio n with the same spatial and time resolution as GRACE [ Sheard et al ., 2012]. The instru ment of GRACE FO is based on that of GRACE and the major change is having LRI along with KBR (ranging instrument used in GRACE) to guarantee the performance of GRACE. The number of star c ameras to measure the attitude of GRACE has been increased from two to three and the configuration of theses cameras ha s been modified to improve attitude determination. SuperSTAR accelerometers used in GRACE are plann ed to be used in GRACE FO to measure the effects of non gravitational forces. Applications and Resolution Requirements of Gravity Measurements Applications of Gravity Measurements Table 1 1. Applications of gravity measurements [ Sneeuw et al. 2005 ]. Scienc e Area Theme Solid Earth glacial isostatic adjustment, co /post seismic deformation, slow/silent earthquakes, plate tectonics, mantle convection, volcanos, core motion, seismic normal modes Hydrology snow, precipitation, ground water, dams, soil moisture, run off, evapo transpiration Ocean mean flow, coastal currents along shelf edges, interaction mean and eddy flow, ocean fronts position, bathymetry, basin scale mass change, deep water formation, bottom currents Sea Level global sea level change monitoring Ice ice mass balance, bottom topography, ice compaction Geodesy precise heights for engineering, GNSS levelling, coastal height reference, sea level monitoring, inertial navigation, atmosphere, planets
21 There has been a great developme nt and improvement in the gravity measurement s of the Earth for the past few decades due to the various satellite missions as introduced in previous section. This would not be accomplished if the gravity measurements were not useful. In fact, these measure ments are widely used in scientific areas as shown in Table 1 1 [ Sneeuw et al. 2005]. The core motivation of these listed scientific areas arises from the idea of observing and studying the planet that we are living in such that we can predict and prepare for the future. Resolution Requirements For each scientific themes listed in Table 1 1, there is a resolution requirement for the gravity solution which is quantified by spatial and time (temporal) resolutions. There are limitations on both spatial and time resolutions with current technology and each space mission aims to achieve specific resolutions for different purposes. Resolution requirement for each science theme is shown in Figure 1 1. This figure also shows the targeting resolutions for GRACE an d GOCE missions. As explained in a previous section, GOCE focused on increasing the spatial resolution by lowering the altitude with a drag free system, whereas GRACE was focused more on increasing the time resolution using low low satelli te to satellite t racking within a reasonable band of spatial resolution. The actual mission measurements showed that GRACE achieved a sp atial resolution of about 3 5 0 km with roughly 2 cm geoid accuracy for monthly solutions which was better than what was expected [ Save et al. 2016]. This wa s due to the advances in the estimation technique and continuously updated a priori information. Also as explained earlier, time resolution can be reduced to daily solutions. However, when daily solutions were computed the spatial resol ution was only 500 km [ Kurtenbach
22 et al. 2009]. In sum, there is a tradeoff between the time resolution and the spatial resolution. Figure 1 1 T ime and spatial resolution requirements [ Sneeuw et al ., 2005 ] It is important to note what limits resoluti on The current limit of spatial resolution is about 100 km and the time resolution is semi diurnal for single pair LL SST missions.
23 To achieve higher spatial resolution, satellites must fly at low altitudes of less than 250 km. This is achievable by apply ing drag free system s but it will require a lot of fuel to maintain the orbit and therefore shorten the mission lifetime dramatically. In terms of time resolution, it is difficult to achieve less than semi diurnal resolution, because to get a full coverag e of the Earth with single satellite, requires more than half a day, even when the satellite is in a low altitude polar orbit. Multiple satellites can be used to solve this problem but it will require a larger amount of cost. Figure 1 2 V arious static gravity solutions in geoid accuracy Nevertheless, it is essential to understand that the accuracy of solution is different from the resolutions. Gravity solutions need to achieve certain accuracy with given time and spatial resolutions. Geoid accuracy is often calculated to evaluate the solution. Geoid accuracy simply means the gravity solution expressed in terms of a distance The definition and details of the geoid calculation will be explained in C hapter 2 Figure 1 2 shows the accuracy of the static g ravity solutions, in terms of geoid
24 accuracy versus spatial resolution, derived from the measurements of different missions or combinations of missions. The CHAMP based solution is called EIGEN 2 which has the accuracy below 10 cm [ Reigber et al ., 2003]. This solution is only based on the measurements taken by CHAMP up to the spatial resolution of 500 km. Solution s beyond this resolution depends on other exis ting models. GGM03 is the GRACE based solution which used four years of accumulated monthly soluti ons to recover the static gravity field with spatial resolution close to 100 km [ Tapley et al. 2007]. EIGEN 5C is a static solution based on combining the data from LAGEOS and GRACE with spatial resolution of less than 100 km [ Foerste et al. 2008]. EGM_G OC is a GOCE based solution but also used the data from LAGEOS and GRACE [ Bruinsma et al. 2013]. This also has the spatial resolution of less than 100 km. Figure 1 2 shows that the GOCE drag free system which compensate s for the effects of non gravitatio nal forces, provides more accurate gravity solution. Furthermore, this figure shows that the static solution can be recovered by using a combination of data from different missions. This implies that there are already a lot of data collected to recover the static gravity field. Therefore, future gravity missions will focus on precisely recovering dynamic gravity solutions. Goals for Future Geodesy Missions Starting from the GRACE mission, the low low satellite to satellite tracking (LL SST) configuration ha s become a standard tracking system for satellite geodesy missions aiming for spatial resolutions down to 100 km and time resolution as short as one day Therefore, the GRACE FO mission will be launched with same LL SST system. This GRACE FO mission will s oon represent the state of the art for satellite geodesy missions. It is speculated that the LRI will reduce the high frequency
25 measurement noise in GRACE FO, but acceleration noise and aliasing errors will continue to dominate the uncertainty of the solut ion. It is therefore likely that drag free systems or more accurate electrostatic accelerometers will be employed with laser ranging interferometry in future missions beyond GRACE FO. Loomis (2012) already showed, through simulation, that LRI does not impr ove the solution with the same accelerometers used in GRACE, but modest improvements can be made by applying the LRI with a drag free system and lowering the altitude. In contrast to this work, a recent simulation done by Flechtner (2016) showed improvemen ts of 23% on a global scale using LRI and GRACE like accelerometers, although different acceleration noise levels were not considered. Indeed, many other simulation studies for future missions, including those of Wiese (2012) and Elsaka (2014), assume the acceleration noise of the GOCE drag free system and the accuracy of LRI based on the same speculations regarding the limiting error sources. Moreover, the benefits of multiple pairs of satellite have also been evaluated These analysis determine design considerations for optimizing gravity recovery mission s for two pairs of satellites in in line configurations, also known as a Bender formation and this showed a positive result [ Wiese et al ., 2012]. The main purpose of using multiple pa irs of satellites is to increase the sampling frequency to detect changes in the gravity field with higher accuracy. Using multiple pai rs of satellites might not be feasible due to its high cost. However, development of small satellites is actively pursued in various groups since small satellites a re less expe nsive, although relatively low accuracy of measurements is still an issue. Motivation Given the necessity of gravity solutions in the broad science community, it is required to continue the work of sat ellite geodesy. For every gravity field recovery, there
26 are spatial and time resolutions to satisfy and a certain level of accuracy is required for the given resolutions. Conseq uently, there has been continuous technology development to improve the gravity solutions. In advance of the development of GPS, the estimation of gravity potential was limited by the discontinuous tracking system. When the HL SST system was established, the gravity solution was limited by the non gravitational forces. The e lectrosta tic accelerometer was developed to measure th e acceleration caused by the non gravitational forces starting from CHAMP. After CHAMP, GRACE initiated the LL SST which measured the range between the two satellites. This new type of data played a huge role in increasing the resolutions and accuracies of gravity solutions. However, the microwave ranging system from GRACE i s speculated to be the limiting factor for high frequency gravity solutions. Therefore, the Laser Ranging Interferometer (LRI) was introduced which is expected to operate with noise smaller by at least one order of magnitude compare d to the noise of the microwave ranging system. The LRI is soon going to be tested by GRACE FO, and if it operates as expected, temporal aliasing error and residual acceleration noise are likely to be the limiting factor s for the gravity field recovery. Drag free system s decrease the residual acceleration noise and therefore increase the accuracy of gravity solution s Consequently, d rag free systems will be required in future geodesy missions and the necessity for drag free systems was already proved through the GOCE mission. Also, t emporal aliasing error s, caused by the under sampling of the dynamic signals, needs to be dealt with. Numerical filter s can be used to r emove this temporal aliasing error but these filters act as a low pass filter s and remove the true high frequency gravity signals along with the error. Thus, multiple pairs of satellites or updat ed de aliasing models may
27 be used to decrease the error caus ed by the under sampling. Because these new technologies and methods are being considered for missions beyond GRACE FO, there have been several simulation studies including Wiese (2012) and Elsaka (2013, 2014) to evaluate them In these simulations, the ef fects of applying an additional pair of satellites and different configurations were explored. However, there is no study done yet to find the optimal acceleration noise level for the drag free system although the GOCE level of acceleration noise was used in many simulations. Therefore, it is necessary to explore the effects of the acceleration noise level of drag free system s for different magnitudes of aliasing error. Dissertation Contribution and Outline Based on the motivation above, there are three m ain goals for this research The f irst is to develop a toolbox to simulate gravity field recovery for satellite geodesy missions. This includes creating satellite data and an estimation routine. We describe this new optimal estimation routine, incorporatin g several novel features that was developed from the ground up for this study. Second goal is to utilize the actual device noise from UF precision torsion pendulum to validate the performance of the toolbox and observe how these technologies can change the performance of LL SST mission s The l ast goal is to find the optimal residual acceleration noise level for future LL SST missions. Here, we present the results of a new analysis of drag free LL SST missions that consider s a range of residual acceleration noise levels, where all other mission parameters are held constant, corresponding to the latest measurement technology and optimal orbit configurations. We examine drag free acceleration noise performance equivalent to that of the GRACE electrostatic accel erometers, the GOCE drag free system, and the LISA Pathfinder drag
28 gravity potential, expanded in spherical harmonics up to degree and order 60, and evaluate the geopotential estimation error as a function of acceleration noise. Both single polar pair (in line) and two pairs (in line Bender) optimal configurations are analyzed. Moreover, the case s of having realistic and idealistic de aliasing models w ere also explored to observe the case s with different aliasi ng error levels Furthermore, estimation routine was improved by considering the correlation of the noise and acceleration noise effects on the measurement data. There are eight chapters in this dissertation, seven excluding this introduction chapter. Chap ter 2 deals with basic theory and previous studies o f satellite geodesy and drag free systems. Chapter 3 presents the mission architecture which is equivalent to designing a future geodesy mission. Chapter 4 explains the estimation routine and the simulation steps that are developed. Chapter 5 presents the result s and analysis Chapter 6 revisits the estimation routine and shows steps to improve the solution. The UF precision torsion pendulum is introduced in Chapter 7 to validate the toolbox using the actual devise noise data. Chapter 8 provides final remarks and recommendations for future research.
29 CHAPTER 2 BASIC THEORY IN SATELLITE GEODESY AND PREVIOUS STUDIES Gravity field solutions are often represented in spherical harmonic s This chapter explains the relation between these coefficients and the equation of motion of drag free satellites. Also, there are descriptions of how the coefficients are evaluated in different ways. Major concepts that are necessary in the field of satellite geodesy, such as satellite configurations, drag free systems, and major error sources are also explained in this chapter. Representation of Gravity (or gravitational) Potential The gravitati onal potential is a measure of how much work is applied by gravitational force at given distance per unit mass. F igure 2 1 represents the gravitational potential at point with respect to the unknown mass [ Heiskanen and Moritz 1967] The equation for this gravity potential at point is defined by Equation ( 2 1 ) [ Heiskanen and Moritz 1967]. Figure 2 1 Gravit ational potential and spherical coordinate system [ Heiskanen and Moritz 1967 ]
30 (2 1) Here, is the gravitational potential at point is gravitational constant and is density of the unknown mass element Therefore, if this unknown mass is the mass of the Earth and is the average density of the Earth, then is the gravit ational potential at point caused by the gravitational force of the Earth. It is easier to analyze the gravit ational potential in a spherical coordinate system as shown in Figure 2 1 and the converted equation is shown below [ Heiskanen a nd Moritz 1967]. (2 2) (2 3) Now, is expressed in spherical coordinates instead of Cartesian coord i n ates where is the distance between the point and center of mass of the Earth is colatitude and is longitude. Unknown constants, and are the mass and radius of the E arth respectively Operator, is the associated Legendre polynomial of degree n and order m. The associated Legendre polynomial, is defined as shown in Equation ( 2 3 ) [ Heiskanen and Moritz 1967] Most importantly, and are the geopotential spherical harmonic coefficients which define the m agnitude of gravity potential. These coefficients are parameters that need to be estimated and are often referred to as the gravity field solution. The g ravity field represents the gravitational force applied per unit mass which is similar to gravit ational
31 potential but it represent the force itself instead of the work done by that force. The d egree and order of these coefficients define the number of parameters used to define the gravit ational potential of the Earth. Hence, the larger the degree and order of coefficients, the higher the spatial resolution. The a cceleration of point P due to the gravitational force of the Earth can be calculated by taking the gradient of gravit ational potential as shown in Equation ( 2 4 ) [ Heiskane n and Moritz 1967]. This equation is equivalent to the equation of motion for pure drag free satellites. (2 4) Analysis of Gravity Field Figure 2 2 Reference ellipsoid, geoid and topography [ Barthelmes 2013 ] The d efinition of the geoid is essential in understanding gravity field solutions which is represented in Figure 2 2 [ Barthelmes 2013 ]. The g eoid is an equipotential surface that meets the und isturbed surface of the ocean. An e quipotential surface is not the surface where the magnitude of g ravitational potential is constant. I nstead it means
32 that all the gravity vectors are orthogonal to the surface. The g eoid height is often calculated to eva luate the accuracy of gravity solutions. The g eoid height is the height between the geoid and the reference ellipsoid where the reference ellipsoid is the fixed ellipsoid which best approximates the geoid. The g eoid height, can be calculated by t aking the difference between the gravit ational potential of geoid and the reference ellipsoid, and then dividing the difference by T he equation for the geoid height [ Barthelmes 2013 ] is, (2 5) The coefficients and are of the geoid and and are coefficients of reference ellipsoid. Equation ( 2 5 ) shows that the geoid height can be computed for all range s of colatitude and longitude once the spherical harmonic coefficients are computed. For analyzing dynamic gravity field solutions, the geoid height variation must be calculated which is similar to the geoid hei ght calculation. Since the dynamic gravity field solution is defined by the magnitude of change from the static model, and represent dynamic solutions as defined in Equation ( 2 6 ) [ Barthelmes 2013 ]. (2 6)
33 Here, is the variation of the geoid height, which represents the geoid height change due to the dynamic solution. Geoid height and geoid height variations are functions with respect to the colatitude and longitude. In addition to viewing the result in terms of colatitude and longitude, it is also required to view the result in terms of the degree of the coefficients. The g eoid deg ree difference (GDD) can be calculate d by taking the root sum square of the dynamic coefficients and multiplying by the radius of the Earth. The g eoid degree difference ( ) provides the result in terms of length (distance) and the equation is provided below [ Kim 2000 ] (2 7) C hange in surface mass density and equivalent water height ( EW H ) are also often used to represent the solutions of the gravity field. The c hange in surface mass density ( ) shows how much the mass has changed over a specified surface area. Equation (2 8) provides the calculation of surface mass density [ Wahr et al. 1998 ] (2 8) Here, is the average density of the earth and is a load Love number of degree n which is a defined constant. F rom here, EW H ( ) can be simply calculated by dividing the change in surface mass density by density of water ( ) as in Equation 2 9 [ Wahr et al. 1998 ]
34 (2 9) EW H shows how much the mass has changed in units of water height variations. This computed EW H is usually presented by a color map of the Earth which directly helps us to visualize the gravity field solution. The same EW H can also be converted to a numerical result by taking the root mean square of the EW H over the entire or specific area of the map. However, this is not equally weighted RMS. At high latitudes, surface area of a grid is smaller compared to surface area of a grid at low latitudes. Therefore, when spatial RMS is computed, the appropriate weights must be given according to the surfa ce area of that grid. Satellite Formations (a) (b) Figure 2 3 (a) H igh low satellite to satellite tracking and (b) L ow low satell ite to satellite tracking There are different types of satellite formations considered to be used in satellite geodesy. These formations are all based on the continuous tracking system involving GPS satellites. GPS is the most accurate tracking system that exists and provides absolute position measurements. As seen in Figure 2 3, a single satellite using GPS
35 tracking is called high low satellite to satellite tracking (HL SST). Low low satellite to satellite tracking (LL SST) requires additional satellite as seen from Figure 2 3. LL SST continuously measures the range b etween the two satellites This relative position measurement is more precise than the GPS absolute position measurement accuracy. Moreover, the range measurement also provides essential information related to the gravity field of the Earth since the rang e is directly affected by the mass variations. Singe Pair LL SST There are numerous types of LL SST system configurations. One of the important fact about one pair LL SST system s is that for any configuration, the satellites must fly in a polar orbit in or der to get full coverage of the Earth. The three most common configurations are in line, pendulum and cartwheel as shown in Figure 2 4 [ Elsaka et al. 2013]. The in line configuration is a satellite pair follow ing one after another, in the same orbit such that the range between the two satellites varies only along the flight direction of the satellites. The GRACE mission uses this in line configuration. A big disadvantage of th e in line configuration is the lack of East West sensitivity, which causes error stripes when viewed on the surface of the Earth Therefore, to overcome this issue, the pendulum configuration was introduced. As seen from the figure, both the pendulum and the cartwheel configurations use two different orbits for each satellite. In the c ase of the pendulum configuration, the orbit for the following satellite is shifted by the right ascension of the ascending node (RAAN) angle from the orbit of the leading satellite. This pendulum configuration measures the range such that the range measur ements are sensitive to both East West and North South directions. The c artwheel configuration has two satellites flying in different orbits but
36 with the same inclination and RAAN angles, such that their flight direction has less effect on the range of the satellites. I n line (GRACE) Pendulum ( GRACE FO ) Cartwheel Figure 2 4 Single pair LL SST configurations [ Elsaka et al. 2013] Double Pair LL SST As was introduced earlier, it is likely that future geodesy missions will use two pairs of satellites in an LL SST configuration This will increase the sampling frequency such that the aliasing error will be reduced. Also, by using two pairs of satellites, the coverage of the Earth can be complete d more thoroughly or in less time eve ntually opening the possibility of improving the time and spatial resolution s Similar to the single pair case, there are numerous double pair LL SST system configurations and the most common configurations are shown in Figure 2 5 [ Wiese et al. 2012; Els aka et al. 2013; Elsaka et al. 2014]. Three of these configurations are based on the in line configuration. Wiese (2012) applied a Monte Carlo analysis to find the optimal orbit for the second pair of satellites which covered all of the different double in line configurations. The i n line M configuration is when the two pairs are in the same orbit
37 but separated from each other by different mean anomaly. Similarly, the in line configuration has two pairs in different orbits that are separated by different RAAN angles. The i n line Bender configuration has an inclined orbit for the second pair along with the polar pair. Lastly, the advanced pendulum configuration is similar to a single pair pendul um configuration but uses pairs of single pair pendulums across from each other as shown from the figure. In line M In line In line Bender Adv. Pendulum Figure 2 5 Double pair LL SST configurations [ Wiese et al. 2012; Elsaka et al. 2013; Elsaka et al. 2014] Optimal Configurations A s imulation study was done by Elsaka (2013) to find which of these configurations bring s the most accurate solution. It was found that the pendulum configuration is the optimal configuration for a single pair L L S ST. However, a Monte Carlo analysis must be done for this configuration to find the optimal orbits, similar to how the in line Bender configuration was optimized by Wiese (2012). For the double pair LL SST, the in line Bender configuration was found to be the configuration that results in the least error in the solution.
38 Although the pendulum configuration was more suitable than the in line configuration for a single pair LL SST, this diss ertation research used the in line configuration. First of all, finding the optimal configuration is out of the scope of this research. Second, th is work also uses the in line Bender configuration for the two pair LL SST. Therefore, the in l ine inclined pa ir can be added to the in line polar pair greatly simplifying the computation Lastly, in line configuration is the only configuration tested configuration by the GRACE mission and this research desires to show the e ffects on the solutions for mi ssions si milar to GRACE. Drag free System s Drag free systems [ DeBra et al., 2011] offer the most promising technique for reducing acceleration noise below the level of accelerometers by eliminating the noisy electrostatic suspension force s applied to the test mass es. A micro propulsion system is used to keep the host satellite centered with respect to its test mass so that the test mass and spacecraft follow a near perfect gravitational orbit. An example of drag free system in shown in Figure 2 6 [ Nguyen 2016]. Th e t est mass is shielded by a housing from non gravitational forces and thruster s or a suspension system keeps the test mass centered. A drift mode accelerometer is a new hybrid technology that reduces the acceleration noise by cycling the electrostatic sus pension system and uses a laser gravitational acceleration without the need for external propulsion or drag free control [ Conklin, 2015]. Drag free technology was used by the Gravity Field and Steady State Ocea n Circulation Explorer (GOCE) (2009 ~ free system directly cancels the perturbing effects of atmospheric drag, it was able to fly at a low altitude of about 250 km. The
39 measured acceleration noise of the drag free GOCE mission was two orders of Touboul et al. 2012]. Advanced drag free systems with accuracies even better than that of GOCE have been developed for future gravitational wave observatories [ Dolesi et al., 2003]. One such system is being tested in 2016 by the LISA Pathfinder mission [ Armano et al., 2015]. The first results of the LISA Pathfinder mission s howed that the measured acceleration noise was more accurate than the acceleration noise requirement [ Armano et al. 2016]. Drag free small satellite platforms are also being examined, which may reduce the cost of this technology and make multiple satellit e pairs more feasible [ Nguyen et al. 2015]. In the future, drag free low low satellite to satellite missions could improve time varying gravity field estimates beyond those provided by GRACE and GRACE FO. Consequently, it is important to explore the effec ts of different drag free systems on the gravity field solutions to determine the acceleration noise level that is optimal with regard to both performance and cost. Figure 2 6 Drag free system [ Nguyen 2016 ]
40 Major Source of Errors There are three major sources of error in gravity field recovery : (1) device measurement noise, (2) temporal aliasing error and (3) non gravitational forces. For LL SST system s device measurement noise comes from measuring the position and the range of the satellite pai r. The GPS tracking of CHAMP, GRACE and GOCE had 1 cmHz 1/2 orbit determination accuracy. The K Band Microwave Ranging (KBR) system in GRACE had about 1 mHz 1/2 range accuracy. The Laser Ranging Interferometer (LRI) is expected to perform with accuracy of about 80 nmHz 1/2 All three numerical values are true for the frequencies between 1 mHz ~ 10 mHz. Temporal aliasing error is caused by under sampling of an uncertain dynamic model. Since the mass of the Earth is always changing there is no way to avoid t his error. There are four possible ways to decrease this error. The f irst way is to change the configurations of satellites from collinear shape (in line) to other shapes in a way that it increases East West sensitivity such as the cartwheel or the pendulum configurations as it was explained in the previous section. A s econd way is to develop more accurate nominal dynamic models (or a priori information), also known as de aliasing models. A t hird way is to increase the sampling frequency by adding mo re pairs of satellites. There were several analyses including Wiese (2012) and Elsaka (2013) which showed a reduction in the aliasing error by adding another pair of satellites. The l ast way is to use numerical filters to smooth the solution Wahr (2002). In this dissertation to decrease the temporal aliasing error, a second pair of satellites will be added and numerical filters will be used. De aliasing models will also be used and both case s of using realistic and idealistic de aliasing models are explored.
41 D eveloping more accurate de aliasing models and searching for optimal configurations of satellites are out of the scope of this research. Non gravitational forces also known as a cceleration noise affects t he equation of motion of the satellites directly. Three major non gravitational forces are solar radiation pressure, atmospheri c drag and thrust forces. As was explained in earlier section, acceleration noise can be compensated for, by using electrostatic accelerometers or drag free systems. Even after acceleration error is compensated, a residual noise of the drag free system persists This residual noise is caused by the position error of proof mass and noise from the actuators among other sources Table 2 1 shows the noise level s for each existing drag free technology. Table 2 1. Technology and associated residual acceleration noise [ Flury et al ., 2008; Touboul et al. 2012; Antonucci et al. 2012] Technology Residual Acceleration Noise ms 2 Hz 1/2 Electrostatic Accelerometer ~10 9 ( CHAMP, 2000) ~10 10 ( GRACE, 2002 ) Gradiometer drag free control ~ 2 10 12 ( GOCE, 2009 ) True drag free control ~ 3 10 14 ( LISA Pathfinder, 2015)
42 CHAPTER 3 MISSION ARCHITECTURE Before we begin the estimation, there are m any parameters to decide. M ainly we must choose the orbits of the satellites. They must be chosen to match those considered in future missions. Therefore, previous missions and simulation studies are surveyed to choose such parameters. Also, technology for the mission must be chosen based on the latest existing ones with provided noise magnitudes. The mission design has been motivated throughout th is chapter and the method of creating simulation noise will be introduced. Orbit Selection varying gravity field, it is most likely that future Earth geodesy missions will continue to use the low low satellite to satellite tracking approach. We therefore choose to use a measurement scheme Since the microwave ranging system will soon be replaced by the LRI of GRACE FO, we will also assume a satellite to satellite ranging system with a measurement noise equal to that expected for the LRI. Moreover, to focus on future missions, gravity field solutions of two pairs of satellites are considered together with the single pair solution to compare the effects of acceleration noise with differ ent magnitudes of temporal aliasing error. Orbit selection for two pairs of satellites must be chosen for the simulation. There have been simulation studies performed to determine the optimal orbits for two pairs of satellites in gravity field recovery [ Wi ese et al. 2012]. It was concluded that the inclination angle is the most significant parameter that affects the gravity field solution.
43 According to this study, if the first pair of satellites is flying in a polar orbit, then the optimal orbit for the se cond pair is inclined at an angle of 72 degrees. This is based on the fact that both pairs are flying in near circular orbits (eccentricity of 0.001). The optimal altitude (semi major axis) is more difficult to choose since there is a tradeoff between mis sion lifetime and the spatial resolution. Lower altitudes exhibit higher atmospheric drag and therefore require greater fuel consumption to compensate for the drag. At the same time, lower altitudes enable measurements to be made closer ace, making them more sensitive to higher order spherical harmonics of the gravitational potential. Future Earth geodesy missions are likely to focus on spatial resolutions up to 100 ~200 km, because the needs of different science fields are concentrated in this region [ Sneeuw et al 2005 ; Pail et al ., 2015 ]. An altitude of 350 km appears to be a good compr om ise between maximizing the lifetime of the mission and maximizing the resulting spatial resolution. At an altitude of 350 km, the lifetime of the satel lites can be over 10 years with the same mass propellant fraction as GRACE [ Wiese et al. 2012] Temporal resolution of GRACE based gravity field models has been improved up to a single day using a Kalman filter approach, although this requires a priori i nformation derived from geophysical models [ Kurtenbach et al. 2009]. Institutions responsible for GRACE solutions; Geo Forschungs Zentrum (GFZ) at Postdam Germany [ Dahle 2013], Center for Space Research (CSR) at Austin Texas [ Bettadpur 2012], and the Je t Propulsion Laboratory (JPL) at Pasadena California [ Watkins and Yuan 2012], compute monthly mass variations but with greater spatial resolution than the daily solutions. Considering the current state of gravity field solutions, it is plausible to think of time
44 resolutions shorter than one month since we are considering additional pair of satellites, and the science needs of most fields is focused on solutions with time resolutions of two weeks or longer [ Sneeuw et al 2005]. Hence, the time resolution (o r repeat period) of the gravity field solution, which is equivalent to the accumulated data length used to recover the gravity model, is chosen to be two weeks (half a month) for this simulation. Other orbital parameters such as Right Ascension of the Asc ending Node (RAAN) and argument of perigee do not have much effect on the solution [ Wiese et al. 2012] and therefore they are chosen to be similar to those of the actual GRACE mission as seen from the Table 3 1. (a) (b) Figure 3 1 Separation distance between the leading and trailing satellites as a function of time for ( a ) the polar pair and ( b ) for the inclined pair The separation distance between the two satellites is another parameter that must be chosen. The two GRACE satellites are separ ated by a distance of 220 km on average. For this simulation work, we assume that the separation distance is periodically corrected to keep a separation distance of 220 km equivalent to GRACE. The separation distance of both pairs of satellites are shown i n Figure 3 1. The graphs show the separation distance between polar and inclined pairs of simulation satellites
45 for the research, and the range varies from 210 km to 230 km. In average, the range variation is 16.7 m at frequency of 100 mHz and 729.1 m at f requency of 0.1 mHz. Ground Track Patterns It is important to make sure that the ground track patterns of satellites have the desired repeat period. This type of ground track is known as a repeating ground track [ Wiese et al ., 2012]. This way the data will not have a bias towards any specific region of the Earth. Repeating ground tracks are defined by over the surface of the Earth inclination, eccentricity, and semi major axis, as well as with the targeting repeat period in sidereal days, which is slightly shorter than the repeat period in calendar days. Vtipil and Newman (2012) discussed several methods to modify the ground tracks using the above parameters, and the method chosen for this analysis among the discussed methods is called the Flower Constellati on. This is a simple method that considers only term) in a two body (Earth and satellite) approximation. This method computes the semi major axis, when the desired inclination, eccentricity and repeat period are pro vided. The semi major axis can be computed by solving the following equation [ Vtipil and Newman 2012] (3 1)
46 Here, is the semi major axis, is the inclination angle, is the Earth radius, is the second order zonal effects, is the n umber of sidereal days the Earth completes during the period of repetition, is the number of revolutions along the orbit in one period of repetition, and is the and are intermediate variables. (a) (b) Figure 3 2 Ground tracks of (a) polar pair and (b) inclined pairs of satellites. As seen from the Equation (3 1) this is a nonlinear problem. Therefore an initial value for the semi major axis is required. Appropriate values for inclination, eccen tricity and repeat period are discussed in the previous section. The rest of the parameters are well known constants. Since this method is an approximation, the resulting semi major axis might not provide the exact repeating ground track. The most accurate repeating
47 ground track was selected empirically by varying the semi major axis, calculated using Equation ( 3 1), in increments of 2 km and actually computing the positions of satellite with full spherical harmonic coefficients. The ground track patterns f ound are shown in Figure 3 2, and from this figure, it can be seen that they have repeat periods close to 14 days. The final set of orbit parameters used in this analysis is summarized in Table 3 1, including the altitudes calculated by ground track adjust ment. It is worth noting that a Modified Flower Constellation method that considers the additional term, also exists [ Vtipil and Newman 2012]. However, this modified method is more complex in terms of computation and produces a semi major axis that is only different by few meters compared to the result from the standard Flower Constellation method. Consequently, the modified approach was not used in this analysis. Table 3 1. Summary of orbit parameters used for both the polar pair and inclined pair of satellites Parameters Values Semi major Axis of Pair 1 (average altitude) 6734.58 km (356.45 km) Semi major Axis of Pair 2 (average altitude) 6732.69 km (354.56 km) Eccentricity 0.001 Inclination of Pair 1 90 Inclination of Pair 2 72 Argument of Perigee 133 Right Ascension of the Ascending Node (RAAN) 155.8 Separation Distance 220 km Repeat Period (Temporal Resolution) 14 days
48 Technology and Associated Noise The analysis presented here assumes that state of the art measurement technologies will be used in future geodesy missions. Since this analysis is based on a simulation, appropriate noise models for each measurement technology are required. Table 3 2 lists the three measurements required for low low satellite to satellite tracking missio ns and the technologies that are assumed for each measurement in the simulation. Attitude measurements and associated noise are neglected for this research. Table 3 2 Measurement technologies used in this analysis Measurement Technology Satellite position Global Positioning System (GPS) Satellite to satellite range Laser Ranging In t erferometer (LRI) Spacecraft non gravitational acceleration Drag free System Position and Range Noise Models GPS is the most accurate continuous position tracking system available, and it is the most widely used technology to track positions of satellites. The orbit determination accuracy for CHAMP, GRACE and GOCE was computed to be close to 1 cm [ Kang et al. 2006]. Accordingly, we modeled GPS position measurement noise simply as white (Gaussian) noise with a standard deviation of 1 cm Hz 1/2 This noise is added to the simulated spacecraft positions, to all three directions; along track, radial and cross track. The microwave ranging system used in GRACE has an accuracy of ~5 m Hz 1/2 for frequencies between 1 and 10 mHz [ Sheard et al. 2012]. The measurement noise of the LRI is expected to be ~80 nm Hz 1/2 for frequencies between 2 and 10 mHz [ Heinzel
49 et al. 20 12]. The official noise requirement which is used as the LRI noise in this simulation is shown in Table 3 3. As shown in Table 3 3, the model was assumed to be true for larger band of frequencies and also we are assuming that the LRI performance reaches it s requirement. The method used to generate colore d noise will be discussed in th e following section Acceleration Noise Models Table 3 3 LRI and Acceleration n oise models used in this analysis Technology Noise Accuracy ( ) Laser Ranging Interferometer GRACE electrostatic accelerometer (along track and radial directions) (cross track direction) GOCE drag free system ( all directions ) LISA Pathfinder drag free system ( all directions ) T wo major non gravitational forces acting on low Earth orbiting satellit es are solar radiation pressure and atmospheric drag. These forces must be measured by accelerometers or directly compensated for with drag free systems. When accelerometers are used, acc eleration noise is caused by the uncertainties in measuring these non gravitational forces acting on the spacecraft. However, when a drag free system is employed, the acceleration noise is caused by residual forces acting on the test mass. This residual no ise is caused by residual coupling to the relative
50 motion of the spacecraft and test mass (stiffness), by noise from the test mass actuators, and by temperature, pressure, and electric and magnetic field fluctuations. To examine the effects of different dr ag free residual acceleration noise levels, we apply three spectra of acceleration noise, each representing a distinct, existing technology. The technologies evaluated are, the electrostatic accelerometer used in GRACE, the drag free system used in GOCE, a nd the LISA Pathfinder drag free system. The noise levels used in the simulation for the electrostatic accelerometer [ Flury et al. 2008] and the two different drag free systems [ Touboul et al. 2012; Antonucci et al. 2012 ] are shown in the T able 3 3 The drag free system with an acceleration noise equivalent to that of the GRACE electrostatic accelerometers, which has an accuracy of 10 9 ms 2 Hz 1/2 in cross track direction and 10 10 ms 2 Hz 1/2 in along track and radial directions, was applied in the simul ation. The GOCE drag free system was applied only on the along track direction. The acceleration measurement from GOCE had the accuracy of 2 10 12 ms 2 Hz 1/2 with a frequency bandwidth from 5 to 100 mHz [ Touboul et al. 2012]. However, the residual acceleration noise of the drag free system was larger. Nevertheless for this simulation, the GOCE drag free system was applied with the accuracy of acceleration measurement on all directions. For LISA Pathfinder drag free sys tem, a requirement for the noise magnitude wa s 3 10 14 ms 2 Hz 1/2 over a frequency range of 1 to 100 mHz [ Antonucci et al. 2012]. Although the actual residual acceleration noise was found to be smaller, this noise level requirement was used in the simul ation. Since the noise spectra is not provided for other frequencies, it was assumed for both drag free systems that the noise magnitude slightly increases for frequencies less than 1 mHz. For each simulation, we considered a frequency band that
51 ranges fro m 0.1 to 100 mHz. The lower bound limited by the size of the data segments and the higher bound is set by the sampling rate. All the acceleration noise models used in this simulation are shown in Table 3 3 Generating Colored Noise Most measurement noise s pectra can be modeled as where GPS, LRI, and accelerometer noises can be modeled accurately with linear combinations of such curves. Kasdin (1995) and Stoyanov (2011) have provided useful methods to create colored noise of the form Th is method is described by Equations (3 2) through (3 3) [ Kasdin 1995; Stoyanov et al. 2011]. (3 2) for (3 3) H ere, is the impulse response function, is the white noise random field, is the resulting random field, are continuous independent variable s and is the frequency. Equation (3 3) is the discrete form of Equation (3 2), which is the discrete convolution. The discrete convo lution can be easily computed by taking the fast Fourier transform (FFT), multiplying the impulse response function and the white noise random field in the frequency domain, and finally transferring the product back to the time domain using the inverse fas t Fourier transform (IFFT). The simulated noises using this approach and the associated models for the LRI and drag free systems are shown in
52 Figure 3 3. We tried to make the noise that accurately represents the model using the method explained in this sec tion. However, there are some discrepancies especially at low frequencies for the LRI noise, where the slope is The generated LRI noise is about one order of magnitude smaller than the provided model at low frequencies. Summary of the Simulation Desi gn T here has been development s in various technologies and methods that were mentioned above Since the main objective of this research is to explore the effects of different acceleration noise levels and hence find the optimal acceleration noise level for future LL SST geodesy missions, it is required to set other tec hnologies in the simulation as the most a ppropriate ones. Furthermore, the optimal orbit and configuration must be chosen to reduce the error as much as possible. GPS is the most accurate continuous position tracking system available and this was used in recent missions such as GRACE and GOCE. For range measurements, the microwave ranging system was used in GRACE. However, the laser rangin g interferometer (LRI) known to be more accurate and soon will be tested in GRACE FO. In order to reduce the temporal aliasing error, additional in line pair will be simulated in an in line Bender configuration and be compared with a single in line pair. A lso, both realistic and idealistic de aliasing models will be used to explore the effects of acceleration noise for both cases. A d rag free system was used i n GOCE but it had only 1 degree of freedom (DOF) This simulation will consider 3 DOF drag free sys tem s with different levels of acceleration noise. Table 3 4 summarizes the simulation design with respect to the error sources and launched satellite geodesy missions.
53 (a) ( b ) ( c ) (d) (e) Figure 3 3 T he noise model and numerically generated noise for (a) LRI, (b) GRACE along track and radial acceleration noise, (c) GRACE cross track acceleration noise, (d) GOCE acceleration noise, and (e) LISA Pathfinder acceleration noise
54 Table 3 4 Summary of simulation design with respect to error sources and launched missions Error Source Launched Missions Simulation Device m easurement GPS K Band Ranging GPS LRI Temporal a liasing In line configuration Realistic de aliasing model In line & in line Bender configurations Realistic & idealistic de aliasing models Acceleration Electrostatic accelerometer Drag free system (1 DOF) D rag free system (3 DOF)
55 CHAPTER 4 ESTIMATION AND SIMULATION PROCEDURE A toolbox to simulate satellite tracking for geodesy missions was created that can compute measurement data of the satellites, and also estimate the gravity field of the Earth using these measurements. We followed classical approaches that exist in several text books and journal articles that fit into our research environment. This kind of toolbox already exists in other institutions where they have been worki ng on satellite geodesy for decades However, to understand the full view of this field, and also to bring unique research in the future, this toolbox was created. G ravity models used in this simulation are explained. The methods of creating data and the estimation process are presented throughout the chapter. Gravity Model Selection Since this analysis is based on simulation rather than actual mission data, we mission. Nominal data is also required, which is the initial approximation of the actual data based on a priori i nformation, which also serve as de aliasing models. As mentioned earlier, both static and dynamic gravity field solutions exist. Typically, five independent dynamic models comprise the major sources of mass variations of the Earth. They are the ocean tide, ocean, atmosphere, hydrology and ice models. Ocean tide is the ocean variation mainly caused by the Moon and the Sun. The ocean and atmosphere models represent mass variations of the ocean and atmosphere respectively. Hydrology and ice models refer to the water variations over the land that include soil moisture, groundwater, rivers, reservoirs, snow, etc.
56 Fortunately, there are several different sources for each of the five dynamic gravity field models, some of which are generated by independent research groups. This allows us to use one such source as the true model and another for the nominal model. The models chosen to create true and nominal data are shown in Table 4 1 Table 4 1 Static and dynamic gravity models used in this analysis True A True B Nominal Static EIGEN 5C EIGEN 5C EIGEN 5C Dynamic Ocean Tide EOT11a FES2004 FES2004 Ocean AOD1B 05 (OMCT) AOD1B 01 (PPHA) AOD1B 01 (PPHA) Atmosphere AOD1B 05 (ECMWF) AOD1B 01 (ECMWF) AOD1B 01 (ECMWF) Hydrology Updated ESA Updated ESA None Ice Updated ESA Updated ESA None The static gravity field model, EIGEN 5C is produced by Geo Forschungs Zentrum at Potsdam (GFZ) in 2008 [ Foerste et al. 2008], and EIGEN 5C is derived from a combination of GRACE and LAGEOS mission data with 0.5 0.5 deg gravimetry and altimetry surface data [ Foerste et al. 2008]. The ocean tide model EOT11a is derived from an empirical harmonic analysis based on multi mission satellite altimetry data and published by Technische Universitt Mnchen (TUM ) [ Savcenk o and Bosch 2012; Rieser et al. 2012], and FES2004 is derived by a modeling approach based on hydrodynamic and assimilation models and published by French Tidal Group (FTG ) [ Lyard et al. 2006]. The atmosphere and ocean de aliasing product AOD1B is used for the GRACE data analysis. This model consists of the gravity variations over six hour
57 periods in atmosphere and ocean with respect to the mean gravitational field. The AOD1B Release 1 product (AOD1B 01) used the ocean model called Pacanowski, Ponte, Hir ose and Ali (PPHA) to obtain the ocean variations and European Center for Medium range Weather Forecast (ECMWF) to obtain the atmosphere variations [ Flechtner 2007]. The AOD1B Release 5 product (AOD1B 05) used the ocean model called Ocean Model for Circul ation and Tides (OMCT) to obtain the ocean variations and ECMWF data to derive the atmosphere variations [ Dobslaw et al ., 2013; Flechtner et al. 2014]. Both AOD1B 01 and AOD1B 05 used ECMWF to obtain the atmosphere variations, but ECMWF was updated betwee n the two releases and some unknown changes were made, resulting in different atmosphere models. The AOD1B 05 is the product of five generations of model improvements from the AOD1B 01. The AOD1B product is provided in spherical harmonic coefficients up to degree and order 100. Hydrology and Ice models are provided by the European Space Agency (ESA) [ Gruber et al. 2011; Dobslaw et al. 2015]. They created their first version of simulated time variable gravity field with a period of 6 hours that includes at mosphere, ocean, hydrology, ice and solid earth components. This data is provided in spherical harmonic coefficients up to degree and order 180. In 2014, ESA generated an updated version of their time variable gravity field models. The hydrology and ice mo dels used were extracted from this updated ESA model. We also use two sets of true models, denoted A and B. The True A set represents the realistic case where we assume that de aliasing models have the accuracy of existing products. The True B set of models were used to generate measurements, assuming that the de aliasing models are more accurate such that the
58 aliasing error is only caused by the dyna mics of the hydrology and ice signals. The accuracy of these idealistic de aliasing models is difficult to achieve in practice, but such models have continued to improve over the past several years. It is therefore enlightening to examine solutions where t here are less aliasing errors caused by the inaccuracy of these de aliasing products. Hydrology and ice models are not used as a part of the de aliasing product, because this research evaluates the accuracy in estimating the true hydrology and ice signals and we want to remove any bias that comes from pre knowledge of the solution. Creating Measurement Data There are six orbital elements; semi major axis, eccentricity, mean anomaly, RAAN, inclination and argument of perigee These orbital elements are chose n according to the mission architecture from Chapter 3 and they can be used to compute the initial position and velocity vectors. Therefore starting from this initial position and velocity, satellite measurements can be obtained. First, acceleration cause d by the gravitational field of the Earth can be computed for given position and velocity vectors by applying the following equations [ Montenbruck and Gill, 2005 ] , for ( m = 0 ) for ( m = 0 ) (4 1)
59 Once the acceleration is computed, this acceleration is integrated using Runga Kutta 4 th order (RK4) method to compute the position and velocity vectors of the next epoch T his procedure is repeated until all the measurements are computed. Estimation Method T he most standard estimation technique to apply is the best linear unbiased estimator (BLUE) which uses equations equivalent to the weighted least squares estimation. However, considering the size of the data and the estimating parameters, it is not possible to perform the estimation with a single block of data and Jacobian (transformatio n) matrix. This would require a large amount of RAM space and will require a long computation time. Accordingly, the partitioned BLUE is a viable candidate. Using the partitioned BLUE, we can avoid the lengthy computation of large matrices by simply dividi ng the data and Jacobian matrices into smaller blocks, and
60 parallel computing the estimation process using these smaller blocks. This reduces, not just the computation time for the estimation, but also the computation time needed to create Jacobian matri ce s For this simulation, each block was chosen to span 3 hours of data. In other words, the arc length was chosen to be 3 hours. Different institutions depending on their estimation method use various arc lengths. JPL, CSR and GFZ institutions use daily arc lengths. We decided to use a shorter arc length because, we estimate fewer parameters compared to the institutions mentioned above. We have reduced the arc length up to the point where the data size for each arc still exceeds the number of parameters and still allows the satellites to complete more than one orbit so that this will not affect the low degree coefficients. There are two types of measurement data as discussed from previous section, inter satellite range and position. The sampling time is chose n to be 10 seconds, which provides a valid frequency for recovering the solution given the desired spatial resolution. During each time step, seven measurements are made, denoted as vector ( 4 2) Here, is the inter satellite range, and ( ) and ( ) are the positions of the following and the leading satellites respectively. The interferometer only measures the relative range, so GPS data must be used to provide absolute positions of the satellites. As mentioned above, the measurement data is divided into 3
61 hour arcs (batches). This leads to 112 batches of data for the entire 14 day repeat period. The measurement equation is shown in Equation (4 3). ( 4 3 ) Here, is the measurement data vector is the Jacobian (state to measurement transformation matrix) is the state vector (estimating parameters) at the initial epoch The variable represents the initial state of the satell ites and represents the spherical harmonic coefficients. The state transition matrix is which converts the static state vector at epoch to the current time Lastly, is the measurement noise vector These equations show that even though the gravity field is time varying, we are estimating the best average approximation over a given 14 day period. The variational equations are used to create the state transition matrix. This method is highly computa tional but it is known to be straight forward and accurate. The general form of the variational equations are as follows [ Montenbruck and Gill 2005] ( 4 4 ) ( 4 5 )
62 Here, is the identity matrix, , and are the position, velocity and total acceleration of the satellite respectively, and is the number of spherical harmonic coefficients. Subscript refers to the portion of the state transition matrix associated with position and velocity. We let and since the problem is nonlinear, we solve it using the Gauss Newton method. The estimation equation then becomes ( 4 6 ) or more compactly, (4 7) Here, represents the measurements created by a priori information and for 14 days of analysis. The measurement noise model is provided by the GPS and LRI instrument performance given in previous sections. The measurement noise covariance matrix used in the estimation process was computed based on these noise models. On the other hand, the non gravit ational acceleration of the spacecraft is not measured directly by drag free satellites. Instead, disturbance accelerations are compensated by thrusters that keep the position of the spacecraft centered with respect to a proof mass. Therefore, the accelera tion noise is not a measurement noise for drag free systems This is in contrast to GRACE or GOCE where the spacecraft acceleration is measured by the on board accelerometers where acceleration noise was measurement noise. In addition to measurement noise, the estimation process also uses state noise to account for the uncertainty of the state equation. For GRACE and GOCE data analysis,
63 empirical or geophysical models were used to estimate the state noise [ Kurtenbach et al. 2009]. In this analysis, an appr opriate uncorrelated state noise covariance matrix , was computed by first averaging the nominal set of dynamic models from Table 4 1 for a time span of two weeks and then scaling this value empirically, such that the estimation process was stable. A ll parameters were given the same weights to avoid any bias. The term which is the initial guess of the parameters, also known as a priori information, is set equal to zero for all arcs. Now Equation ( 4 7 ) can be modified as follows, ( 4 8 ) Here, is a priori information of the parameter vector with the state noise Creating the measurement covariance matrix, that accurately represent the measurement variation is challenging due to the noise correlation. Two approaches were co nsidered to compute these covariance matrices. The first was to choose a white noise covariance that bounds the appropriate colored noise. The second was to use a filter to convert a white noise covariance into the appropriate colored noise covariance. Add ing a filter will produce more accurate estimates and the resulting posteriori statistical error will be more representative of the true estimation error [ Petovello et al. 2008]. On the other hand, a filter will require extra computations. Since this will increase the computation time and cost, a bounding white noise covariance was used. A filter is considered in later chapter of this dissertation. Equation s (4 9) and (4 10) are the BLUE solutions with a priori information. They are equivalent but Equatio n (4 10) is used to perform the estimation.
64 ( 4 9 ) ( 4 10 ) From this point, the best linear unbiased estimation can be performed simply by using each set of matrices to perform the estimation and then combining the result using the posteriori covariance matrices and estimates to compute the final solution. In this combining process, it was assumed that the each estimates are correlated only for one day period. This is to mimic the daily solutions computed by the other groups as mentioned above. Ideally, this result is equivalent to performing a single best linear unbiased estimation using the whole concatenated matrices. There were two different way s applied to compute the BLUE. F i rst is the well known QR factorization method. Second is the preconditioned conjugate gradient (PCCG) method introduced by Klees (2003). Equation (4 11) presents QR factorization, where is the orthonormal matrix and i s the upper triangular matrix. ( 4 11 ) When QR factorization is applied to BLUE, the estimation equation reduces to as follows. ( 4 12 )
65 T h e advantage of QR factorization is that the complexity of computation reduces dramatically and since is an upper triangular matrix, taking an inverse is avoided. Another big advantage is that the size of matrix reduces from number of rows to columns. In most cases, there are more rows than co lumns, meaning that there are more observations than the parameters need to be estimated. Reducing the matrix size is essential when dealing with a large amount of data. The PCCG method is an iterative process which is presented below with series of equations [ Klees et al. 2003]. F or until ( 4 13 )
66 Here, is an approximation of the matrix This iterative process is truncated while and is chosen empirically. A b ig advantage of PCCG over QR factorization is that it uses less RAM space since all the terms being iterated are vectors, not matrices. Large data computation such as gravity field estimation, requires a lot of RAM space so PCCG becomes a very efficient method. The t wo different approaches produced equivalent results with similar computation time. However, we chose to use, QR factorization since it is easier to implement and modify Also, since the computation was not limited by the R AM space there was no need for the use of the PCCG method This does not mean that the PCCG method is less efficient than QR factorization. The number of iterations we chose empirically took same time as the QR factorization method, but there is a possibility that PCCG can be more time efficient with the same accuracy if it were tuned properly Simulation Procedure The entire simulation procedure is shown in Figure 4 1. There are four parts in this simulation work The f irst part create s the measurement data using true set of g ravity models There are seven parameters to generate, the inter satellite range and the positions of two satellites. These measurement data can be computed simply by integrating the equation of motion of pure drag free satellites and adding GPS, LRI and a cceleration noise as explained thoroughly in previous sections This part must be computed in a sequence, but it does not require a lot of memory. The s econd part of the simulation is creating the transformation matrices and reference measurements using no minal gravity models. The reference data is created in a same manner as the measurement data. However, computing transformation matrices require s a lot of
67 memory and takes long time, requiring as many CPUs as available Fortunately, this part can be comput ed in parallel and therefore, the University of Florida HiPerGator supercomputer was used to implement this part. Third part of the simulation is the estimation. The BLUE is applied to find the gravity solution. Daily solutions were computed first, and the n these daily solutions were used to estimate the final biweekly solution. There are two parts to the estimated solution. One is the orbital elements, which is equivalent to the initial positions and velocities of satellites. Another is the spherical harmo nic coefficients, whic h represents the size of the gravit ational potential. In total, 120,960 observation sets were used to estimate 3,729 parameters. Finally, the last part is to evaluate the solution using techniques brought up in Chapter 2. Figure 4 1 Block diagram of the simulation procedure
68 CHAPTER 5 RESULT AND ANALYSIS Representation of Solution s The solutions were estimated in terms of spherical harmonics up to degree and order 60 using the steps provided in the previous section. The most common way to show such results is by plotting the geoid degree difference (GDD) error and the equivalent water height (EWH). The GDD error is the empirical error between the true and estimated coefficients. It is the root sum square of the differ ence with respect to each degree of coefficients, such that it shows the error of the result in terms of geoid height for every degree of the spherical harmonic coefficients [ Wiese 2011]. Equivalent water height is a physical quantity that shows the mass variations over the entire globe in terms of water height difference [ Wahr et al. 1998]. Spatial root mean square (RMS) error of the EWH solution can also be calculated by subtracting true EWH from the estimated EWH, and then taking the root mean square o f the entire set with appropriate weights based on the surface area of each grid. For all solutions, a numerical filter radius of 300 km or less (depending on the solutions), to remove the aliasing error in the spherical harmonic coefficient estimates. Both non filtered and filtered results are presented throughout the section. Gravity field solution s were computed for the low low satellite to satellite tracking system for three different cases. Case 1, 2, and 3 refers to GRACE, GOCE and LISA Pathfinder residual acceleration noise levels respectively, which is shown in Figure 3 3 and Table 3 3 As described in Table 4 1 there are two different sets of true models, A and B, to app ly different levels of aliasing errors. From this point we will present results
69 using two different sets of true models for all three cases mentioned above and denote them as (A1, A2, A3) and (B1, B2, B3), where A and B refer to the different sets of true models and 1, 2, and 3 refer to different acceleration noise levels. Geoid Degree Difference Error Figure 5 1 shows the errors in geoid degree difference for Cases A1, A2 and A3. The true signal in Figure 5 1 is the true hydrology and ice mass variations converted to geoid height. Here, only hydrology and ice signals were drawn as the true signal since the goal of simulation is to recover hydrology and ice signals among all dynamic signals. As seen from this f igure, all the estimated errors are below the true signal for the low degree coefficients. This is not a new finding, but it shows that the solutions from this simulation are reliable. We can see the increas e in error towards the high degree coefficients, which is caused by various error sources including aliasing error, measurement noise and residual acceleration noise. By applying Gaussian smoothing (GS), this error can be reduced, resulting in more accurate solutions throughout the entire range of coeffi cients as shown by the dashed lines in the figure. A 300 km radius was applied to all single pair solutions. For two pair solutions, same 300 km radius was applied to Case A1, but both 200 and 300 km radiuses were applied to Cases A2 and A3 to find the rad ius that minimizes the error. GS between these two radiuses produced similar results, however 300 km radius had smaller error and so this radius was chosen. Applying the same GS radius for all three cases, allows us to compare the effects of GS. However, i t is true that this radius does not provide best solutions for all cases.
70 Figure 5 1 Errors and Gaussian smoothed (GS) errors expressed as GDD for one pair (top ) and two pairs ( bottom ) using set A true models.
71 Figure 5 2 Errors and Gaussia n smoothed (GS) errors expressed as GDD for one pair (top ) and two pairs ( bottom ) using set true models.
72 From Figure 5 1 it can be seen that the high degree errors are greatly reduced when GS is applied. On the other hand, by applying GS, errors were increased for some low degrees. This is expected since GS removes some of the true signals along with the errors. Also, the solid curves for the single pair solutions on the left have greater magnitudes than the solid curves for the two pair solutions on t he right. This indicates that adding a second pair of satellites reduces errors. In order to quantitatively determine the magnitude of the error reduction, the spatial RMS error is computed later in this chapter Most importantly, Figure 5 1 shows that the error for Case A1 (GRACE acceleration noise), whether it is filtered or not, is always higher than the errors associated with the other two cases. This shows that by reducing the residual acceleration noise level from the GRACE to the GOCE level, the solution is improved. Moreover, Case s A2 and A3 (GOCE and LISA Pathfinder acceleration noise) solutions overlap almost exactly, indicating that there is not much improvement in the solution when the acceleration noise is reduced to below the GOC E level. Figure 5 2 shows the estimation errors in terms of GDD for Cases B1, B2 and B3, as well as the true signal. The true signal in Figures 5 1 and 5 2 are the same, since the same hydrology and ice models are being estimated. Figure 5 2 shows the resu lts of applying the idealistic de aliasing models. By reducing the aliasing error with the idealistic de aliasing models, the contributions of the aliasing error at different acceleration noise levels were revealed. Compared to Figure 5 1 Case B1 does not show an improvement, whereas the other two cases show a large reduction of the error. This means that at the GRACE acceleration noise level, reducing the aliasing error does not improve the solution, indicating that the acceleration noise is dominating. O n the
73 other hand, at the GOCE acceleration noise level, applying idealistic de aliasing models removes a significant amount of the error. This reveals that at GOCE acceleration noise level, the aliasing error is one of the dominating error sources. When ac celeration noise is further reduced to the LISA Pathfinder level, the solution did not show significant improvement. This means that the Case B3 solution is limited by either the measurement noise from LRI and GPS tracking or the aliasing error caused by t he hydrology and ice signals. Gaussian Smoothing with a 300 km radius was applied for Case B1. For Cases B2 and B3, a 100 km radius was applied since the estimated results seem to have very small errors based on the fact that the estimated signal is only s lightly above the true signal beyond degree 50. In other words, GS was almost unnecessary, since the error curves were below the true signal for most of the degree range. Figure 5 2 also shows that when the in line Bender configuration was applied, the sol utions were improved for all three cases. This further reduction of the aliasing error does not change the general conclusion that was drawn for a single pair of satellites. Even when an additional pair of satellites and idealistic de aliasing models were used to reduce the aliasing error, acceleration noise is still not a dominating source of error at GOCE acceleration noise level or below. Equivalent Water Height The geoid degree difference error (or geoid height error) shows the global result in a single one dimensional curve, and therefore it is difficult to draw a complete conclusion from it. To obtain a deeper insight into the estimation errors, we visualize th e solution by drawing global equivalent water height (EWH) maps. The EWH graphs for the solutions discussed above are shown in Figures 5 3 and 5 5 Table 5 1 shows the
74 respective spatial root mean square error of the EWH, which we use to quantify the magni tude of the error reduction. Figure 5 3 pair solution (row 2), double pair solution (row 3), single pair solution with GS (row 4) and double pair solution with GS (row 5). From row 2 to row 5 Case A1 (left), Case A2 (center) and Case A3 (right). Units are in cm of EW H and the figure includes hydrology and ice mass variations.
75 Figure 5 3 Continued. Figure 5 4 EWH error versus latitude plots; Case A1 on left and Case A2 on right. Figure 5 3 Cases A1, A2 and A3. These color maps are drawn with a range of 15 cm. Row 1 is the true hydrology and ice mass variation in terms of EWH. The goal of the estimation process is to recover these signals with minimum error. Rows 2 through 5 in Figure 5 3 have three columns each. Left, center, and right each refers to Cas es A1, A2 and A3 respectively. Row 2 is the estimated solutions for a single pair of satellites and row 3 shows the solutions for two pairs. Row 4 shows the GS applied result for a single pair of satellites and row 5 shows the GS applied result for two pai rs. Both single and double pair solutions for Case A1 show discrepancies compared to the solutions for the other
76 two cases. The solutions for Cases A2 and A3 are almost identical even when a second pair of satellites is used. This result remains true even after applying GS. When compared with the true data (Row 1), we see that true signals are recovered for all cases. However, for Case A1, errors mask most of the true signal. The recovered signals become more pronounced and much of the aliasing error stripe s are removed when an additional pair of satellites is used. However, this addition increases the errors at very high latitudes (both North and South) for all cases, which is shown by the red and blue regions in Figure 5 3 To show this affect more clearly in Figure 5 4 we plot the average EWH errors versus latitude for Cases A1 and A2. Case A3 is omitted because this result is almost identical to Case A2. From the figure, for both cases, when latitudes are greater than 70, the error increased when the seco nd pair of satellites were added. This is because the inclined pair does not cover the high latitudes. Nonetheless, for most other latitudes, the error decreases when applying the Bender formation, because these latitudes are covered by both pairs of satel lites. In addition, Figure 5 4 shows that the effects of GS are greater for Case A1 compared to other two cases even when the same radius is used. Overall, with these EWH results, we are able to observe the amount of error reduction when the acceleration n oise is reduced from the GRACE (Case A1) to the GOCE (Case A2) level. Moreover, the solution with LISA Pathfinder (Case A3) noise, did not show further reduction of the error anywhere in the map. There is one additional comment on GS. Rows 4 and 5 of Figur e 5 3 show that GS eliminates most of the errors although it also diminishes the actual signal. For example, in the center of Greenland, there is an increase in the EWH (red color) over a relatively small area. This signal was recovered in the unfiltered s olutions, but it was
77 removed when GS is applied, because of averaging with the nearby surface. This shows that there is a tradeoff between reducing the aliasing errors and reducing the true signals when applying Gaussian Smoothing. This tradeoff motivated the use of the in line Bender formation and idealistic de aliasing models in this study. They allow the influence of acceleration noise on the solution to be clearly observed without needing to use GS, which reduces the true signal. Similar to Figure 5 3 Figure 5 5 models, which are Cases B1, B2 and B3. The arrangement of Figure 5 5 is the same as that of Figure 5 3 models appear to be more accurate (less discrepancies compared to the true EWH), which is expected since there are less errors caused by the imperfection of de aliasing models. However, Case B1 and Case A1 results appear to be very similar, indicating that the solution did not im prove by applying idealistic de aliasing models at GRACE acceleration noise level. On the other hand, the solutions for Cases B2 and B3 improved compared to Cases A2 and A3 where realistic de aliasing models were applied. This implies that for acceleration noise levels smaller than that of GOCE, the acceleration noise is not a dominating source of error, but the aliasing error is at least applying the in line Bender formation improved all the solutions. However, this improvement is seemingly negligible for Cases B2 and B3, because the single pair solutions exhibit very little errors. Spatial RMS EWH errors can be computed to see if this is actually true. Similar to Figure 5 4 Figure 5 6 shows the average EWH errors for Cases B1 and B2 versus latitude. Case B3 is omitted for the same reason as that for
78 Case A3. Applying the Bender formation again increased the error at high latitudes (70 or above) but decreased the error at mid low latitudes. Also according to the same figure, GS generally improved the solution for Case B1 but not for Case B2, because some of the true signals were filtered. However, this may not be a good comparison since different radii were used when applying G S. Figure 5 5 pair solution (row 2), double pair solution (row 3), single pair solution with GS (row 4) and double pair solution with GS (row 5). From row 2 to row 5 Case B1 (left), Case B2 (center) and Case B3 (right). Units are in cm of EW H and the figure includes hydrology and ice mass variations.
79 Figure 5 5 Continued. Figure 5 6 EWH error versus latitude plots; Case B 1 on left and Case B 2 on right. Figures 5 3 and 5 5 clearly help to visualize the solutions. However numerical values are needed to quantitatively assess these solutions. Table 5 1 shows the spatial root mean square (RMS) error of the EWH solution, shown graphically in Figures 5 3 and 5 5 The hydrolog y and ice (HI) error represents the error over the land and total error is the error over the entire map. According to the table, the total error is smaller
80 than HI error for most cases, which was not visible from the EWH maps. This is because the nominal models for the atmosphere and ocean bottom pressure were used to compensate for the variations over the ocean, but no nominal models were used to compensate for the HI variations. However, this may not always be true because the error of ocean and atmosphe re models may still be greater than the HI variations. Even though the goal of this simulation is to estimate HI signals, it is suitable to look at the total errors because estimation errors can occur over both the land and the oceans. Table 5 1 Spatial RMS errors for all cases in units of cm of EW H HI Error (True: 12.44) Total Error (True: 7.85) Max. Error Raw GS Raw GS Raw GS Case A1 1 Pair 11.72 4.50 11.85 3.30 83.15 55.96 2 Pairs 11.90 4.32 11.58 3.19 97.66 57.96 Case A2 1 Pair 5.47 3.59 4.45 2.11 82.93 52.70 2 Pairs 4.62 3.39 3.86 1.97 46.59 40.88 Case A3 1 Pair 5.45 3.58 4.44 2.10 82.46 52.75 2 Pairs 4.62 3.39 3.84 1.97 44.84 40.84 Case B1 1 Pair 10.32 4.02 10.88 2.92 67.33 52.32 2 Pairs 10.73 3.97 10.75 2.93 91.87 55.70 Case B2 1 Pair 1.97 1.99 1.29 1.21 31.55 49.85 2 Pairs 1.87 1.91 1.18 1.13 28.26 31.56 Case B3 1 Pair 1.97 1.99 1.29 1.21 31.71 49.86 2 Pairs 1.87 1.91 1.18 1.13 28.21 31.54 The numerical values in Table 5 1 show that the errors for Case A1 are larger than for Cases A2 and A3 in all categories. Focusing on the total errors, Case A1 error is more than 7 cm (59%) larger than Case A2. Moreover, Case A2 error is only 0.1 mm (0.2%) larger than Case A3 error. From this we can state that reducing the acceleration
81 noise level from the GRACE level to the GOCE level improves the solution and furthermore, reducing the acceleration noise level from the GOCE level to the LISA Pathfinder level improves the solution by only a negligible amount. The Bender formation improved the solutions for all three acceleration noise cases. However, Case A1 did not improve as much as Cases A2 and A3. To be more specific, Case A1 improved by ~3 mm and the other two cases improved by ~6 mm. This implies that the Bender formation has a larger impact when the solution is less dominated by the acceleration noise. The spatial RMS errors for Cases B1, B2 and B3 show that by applying idealistic de aliasing models, the errors are reduced. In additio n, by applying Bender configuration, this aliasing error was further reduced for all three cases by ~1 mm. Reduction of these errors is expected since improved de aliasing models were used. Again, looking at the numerical results and only focusing on the t otal errors, the Case B1 error is only 1 cm (8%) smaller than the Case A1 error. For other two cases, the errors were reduced more than 3 cm (67%) by applying idealistic de aliasing models and an inclined pair of satellites. This enhances the conclusion th at was made previously, which is that at the GRACE acceleration noise level, acceleration noise is a dominating source of error and when acceleration noise is below the level of GOCE, it is no longer a dominating source of error. Hence, w e have determined that the aliasing error becomes one of the dominating sources once the acceleration noise is reduced to the GOCE level. In addition, even when the idealistic de aliasing models are applied and an inclined pair of satellites is added, the acceleration noise still appears to be non dominating, because the Case B3 solution showed no improvement compared to the
82 Case B2 solution. It is also interesting to compare raw and GS applied solutions. Errors for Cases A2 and A3 with GS are less than 1 c m larger than the unfiltered solutions for Cases B2 and B3. This indicates the importance of GS, even though the effects of GS is limited by the magnitude of the original errors. Discussion We examined the results of a new low low satellite to satellite tr acking mission simulation using a range of residual acceleration noise levels and single and double pairs of satellites (in line and in line Bender formations) All other measurements noise and orbit parameters are kept constant so that the results reflect only the effects of acceleration noise. These residual acceleration noise levels were applied to find the suitable drag free system accuracy for future gravity missions. To determine if acceleration noise levels below that of GOCE would be of value, the o ptimal in line Bender formation by Wiese (2012) and idealistic de aliasing models were used at the same time to reduce aliasing errors. By observing all the results above, including geoid degree difference (GDD) error curves, equivalent water height ( EWH ) graphs and associated spatial root mean square (RMS) error values, we could draw several conclusions. First, it is clear that residual acceleration noise for drag free systems in future geodesy missions should be lower than that of the electrostatic accele rometers used in GRACE, which is close to 1 10 10 ms 2 Hz 1/2 The results above show that with this acceleration noise level, there are larger errors compared to other solutions. Moreover, when the additional pair of satellites and the idealistic de al iasing models were used,
83 the solution only improved by 8%. This indicates that acceleration noise is the dominating source of error when it is equivalent to that of the GRACE accelerometers. Using the GOCE and LISA Pathfinder drag free residual acceleratio n noise levels improved the solutions, indicating that, as expected, such drag free systems should be employed for future geodesy missions to increase the accuracy of the solution. However, the GDD error curves and the EWH plots for these two cases were ne arly identical. Even though the spatial RMS error values indicated that the error for the LISA P athfinder case was actually smaller than the error for the GOCE case, the difference between the two cases is less than 1%. This means that, it is unnecessary t o further reduce the acceleration noise below that of the GOCE level. Furthermore, when the acceleration noise is equivalent to that of the GOCE drag free system, reducing the aliasing error by applying idealistic de aliasing models and the Bender configuration improves the solution by 67%. This shows that at this acceleration noise level, aliasing is one of the dominating source of error, whereas acceleration noise is no longer a limiting one. In addition, at this low level of aliasing error, further reduction of the acceleration noise down to the level of the LISA Pathfinder drag free system did not improve the solution. This indicates that the acceleration noise remains non limiting even when the aliasing error i s reduced. Therefore, when aliasing error is reduced and the acceleration noise is below the level of the GOCE drag free system, the solution is limited by the measurement noise and the variations in the hydrology and ice signals. In sum, these results ind icate that residual acceleration noise lower than the level of GOCE (2 10 12 ms 2 Hz 1/2 ) is not yet necessary with the simulated technology. In
84 other words, residual acceleration noise is no longer the limiting factor once it is near the level of the GOCE drag free system for low low satellite to satellite tracking with LRI, even when the aliasing errors are reduced using the optimal in line Bende r formation and idealistic de aliasing models. However, the EWH spatial RMS errors indicate that when the overall error drops below 0.1 mm, there might be a need for a lower residual acceleration noise. T his conclusion may be related to the conclusion of L oomis (2012). Loomis showed from his simulation that reducing the acceleration noise level from GRACE to the drag free system, the solution was improved although applying LRI instead of the microwave ranging interferometer (MRI) did not show any improvemen ts. However, one cannot make direct comparisons with this work, since the acceleration noise magnitude for the drag terms of acceleration noise. However, it is clear that the same accele ration noise models were used for the case of the GRACE electrostatic accelerometers. But the conclusion that reducing the residual acceleration noise improves the solution is consistent with this conclusion, even though the altitude and range between the two Another conclusion that can be drawn from this analysis is that adding a second pair of satellites only modestly improve d the estimation results Moreover, when idealistic de ali asing models were applied, second pair only improved the result by 1 mm for all cases. Also, it was shown that adding an inclined pair of satellites, increased error at latitudes above 70 degrees. Overall, this might be in contrast to the work of Wiese (20 12), which showed significant improvements to the single polar pair solutions
85 when an inclined pair of satellites was added. However, there are significant differences between the mission architecture analyzed here and that of Wiese (2012), which is listed in Table 5 2 Table 5 2 D ifference s between current and Wiese simulations Current Wi e se Batch Size 3 hours 1 day Estimated Coefficients degree and order 60 degree and order 100 Altitude ~ 350 km 290~300 km Separation Distance 220 km 100 km Repeat Period 14 days 13 days L RI noise 80 nm Hz 1/2 5 nm Hz 1/2 GOCE Drag free noise 2 10 12 ms 2 Hz 1/2 10 10 12 ms 2 Hz 1/2 In general, our gravity field estimates are better than those of Wiese (2012) for analysis the second pair of satellites reduced the error by ~50%, whereas, for this a nalysis, the second pair reduced the error by only ~16% (~1 cm EWH). From Table 5 2 it is clear that the biggest differences between the two simulations are the size of the data batch, number of estimating parameters, and the noise levels. Increasing from degree and order 60 to 100 is equivalent to estimating more than double the number of parameters. Estimating more parameters would cause greater uncertainty in the estimates of each parameter, especially since there are more aliasing errors at higher degr ee coefficients. Moreover, t he LRI noise and altitude would tend to reduce the
86 ater error caused by the higher acceleration noise. There are several explanations for the discrepancy of the results for the double pair case. One of them is that since the LRI noise used in this analysis was more than factor of 10 times higher than that used in Wiese (2012), our solution may have been limited by the LRI noise, such that the second pair could not further improve the solution. Another possible explanation is the difference in the batch size. Reducing the batch size from 1 day to 3 hours may have a large effect on the solution, since having a smaller batch size would make the estimation process more sensitive to the high frequency signals. This may reduce the aliasing error which appears mostly in high degree coefficients. However, the repea t periods used in two simulations are roughly the same, so if the batch size is influencing the result, it would be due to the nonlinear nature of the model or numerical issues. In addition, smaller arc length allow less temporal aliasing errors since the orbits are readjusted every 3 hours. These explanations may obfuscate the benefit of the second pair of satellites. This explanation is consistent with the phenomenon of having lower single pair solution error and higher double pair solution error, compare d with Wiese (2012). Flechtner (2016) applied acceleration noise equivalent to that of the GRACE accelerometers to simulate the expectations of GRACE FO mission. Comparing his result with Case A1 single pair solution from Table 5 1 we find that there is a significant difference between the two results. There are many variables that define the two simulations, but clearly there are several factors that can explain the large discrepancy between the two results. First, altitude for the GRACE FO dedicated simu lation was
87 more than 100 km higher than the altitude used for this research. The repeat period was kept constant only for this research, which is another advantage of drag free systems compared with accelerometer in addition to having lower acceleration no ise. Also, similar to the analysis of Wiese (2012), there were more parameters estimated for the GRACE FO simulation compared to this work. The plots of EWH error verses spherical essively as the degree of the coefficients increases. Therefore, it is likely that if the solution was computed only up to degree and order of 60, as was done in this research, the total average error would decrease compared to the Case A1 solution. Moreov er, having different de aliasing models and forward models would also cause discrepancies between the two analyses. Finally, we conclude that for low low satellite to satellite tracking with LRI, the suitable residual acceleration noise for the drag free s ystem is found to be around 2 10 12 ms 2 Hz 1/2 which is equivalent to the drag free system used in the GOCE mission. Further simulation analysis is recommended to determine the exact upper and lower bounds of the residual acceleration noise that can ma intain the solution errors within say a few percent, as well as whether or not a non drag free electrostatic accelerometer with this level of performance would provide similar results.
88 CHAPTER 6 IMPROVEMENTS IN ESTIMATION State and Measurement Uncertainty It was brought up in Chapter 4, that the state uncertainty was bounded by the number found empirically which is applied to all parameters equivalently. This is because most of the state uncertainty is caused by the unknown uncertainty of de aliasing models and the actual Hydrology and Ice signal. Therefore, it is impossible to design the state noise covariance matrix that accurately represents the state uncertainty unless there exists a model that precisely repr esent s this true signal. This is why the state noise covariance is designed such that the estimation depends mostly on the measurement data, rather than a priori information, which is actually set equal to zero for all arcs. Also, it seems plausible to rep eat the entire estimation process using the state noise covariance matrix computed by the previously estimated result. However, since the estimated solution is not perfect, the errors will restrict the solution from recovering the true signals. Therefore, it is better to allow the estimator to depend mostly on the measurements. On the other hand, knowledge of measurement uncertainty is known This allows the computation of the precise measurement noise covariance (MNC) matrix. However, in Chapter 4, there w ere two assumptions made to simplify the process. These assumption s are (a) considering measurement noise as white and (b) assuming that acceleration noise has no effect on the measurement data. In this chapter, these assumptions are removed to make a more precise MNC matrix and so improve the estimation At first, auto regressive moving average (ARMA) filter i s designed and used to consider the correlation of the noise which
89 was previously ignored. Moreover, the known acceleration noise spectrum i s used t o modify the measurement covariance matrix to account for the unmeasured acceleration noise that actually caused additional uncertainty to the measurements. Methods and results of this improved estimation routine are presented in this chapter. ARMA Filter As seen from Table 3 3 and Figure 3 3, the LRI and acceleration noise s are colored, although up to this point all the noise covariance were considered as white in the estimation routine This is because of the complexity that arises when dealing with correlation of the noise. This noise correlation can be applied in the estimation using an auto regressive moving average (ARMA) filter and Toeplitz systems with the same noise power spectral density (PSD) that were used to generate the measurements [ Klees et al ., 2003]. In other words, if the colored noise and the PSD of the noise model are provided, an ARMA filter can be designed that can filter white noise into colored noise. The ARMA filter is given in Equations (6 1), where is the filtered output and is the white noise with variance [ Stoica and Moses 1997] ( 6 1 ) Equation ( 6 1) can be rewritten as Equation ( 6 2) where and are parameters of the filter [ Stoica and Moses 1997] Once we find these parameters, then we use these values to compute the correlated measurement covariance matrix. There are two special cases. T he filter is called auto
90 regressive (AR) when and when the filter is called moving average (MA) ( 6 2 ) Before estimating the filter parameters, the auto covariance sequence ( ACS ) of the signal must be calculated. This can be done by applying the inverse Fourier transform to the sampled power spectral density (PSD) of the correlated noise. The equation for the ACS computation is provided below [ Klees et al ., 2003]. ( 6 3 ) Here is the auto covariance of lag is the sampled PSD noise function with frequency is the length of the noise samples in the time domain used to estimate the PSD, is the sampling time, and i s the maximum lag which cannot exceed Now, an assumption is made that an AR filter with very large , is equivalent to an ARMA filter as shown in Equation (6 4) [ Stoica and Moses 1997]. With this idea, parameters will be estimated, and then these estimates will be used to estimate the ARMA paramete rs and ( 6 4 ) This AR filter with large is also known as a long AR filter and its parameter s can be computed by applying the Yule Walker method using the ACS computed
91 previously. The relation between the ACS and AR filter parameters is shown in Equation (6 5) [ Stoica and Moses 1997]. ( 6 5 ) Using Equation (6 5), and are computed using a least squares estimation. T he ARMA filter parameters can be computed by the following steps Equation (6 2) can be rearranged as follows [ Stoica and Moses 1997]. ( 6 6 ) If we let then Equation (6 6) can be written as [ Stoica and Moses 1997], ( 6 7 ) In Equation ( 6 7 ) is the provided colored noise vector, and is the approximated white noise vector which was computed by the l ong AR filter from Equation (6 4). Then, by applying least squares estimation of the above equation, can be estimated, which is denoted as Finally, white noise variance is estimated as shown in Equation (6 8) [ Stoica and Moses 1997].
92 ( 6 8 ) Now the ARMA filter parameters are computed by applying Toe plitz systems, to create a well defined measurement covariance Two Toeplitz matrices are defined as follows [ Klees et al ., 2003], ( 6 9 ) These defined Toeplitz matrices can represent the filter that takes white noise ( ) as an input and colored noise ( ) as an output, which is shown below [ Klees et al ., 2003]. ( 6 10 ) Lastly, the MNC matrix is computed as presented below. ( 6 11 ) It is important to make sure that is a positive definite symmetric matrix for estimation to work without numerical error Before applying the above covariance matrix in the estimation process, the designed ARMA filter was tested. The generated noise shown in Figure 3 3, was used to design the ARMA filter through the pro cess explained in this section. Once the ARMA filter is designed, the white noise is created using the variance found during the process. This white noise is passed through the ARMA filter to create the colored noise. This ARMA filtered noise is compared w ith the
93 simulated noise from the previous chapter. Figure 6 1 shows the ASD of generated and ARMA filtered LRI noise with the noise model. The figure shows that the ARMA filtered noise ASD resembles the generated noise. This confirms that the ARMA filter w as designed well enough to be applied in the estimation routine. Figure 6 1 ASD of the LRI generated noise, the ARMA filtered noise, and the noise model. Accounting for Acceleration Noise in the Estimation Routine During the simulation process, several different acceleration noise levels were used each representing a specific case. However, since acceleration is not part of the measurement vector, the information about the acceleration noise was not included when c reating the covariance of the measurement uncertainty, even though the acceleration noise actually affect s the measurements. Therefore, it is worth exploring the effects of building an MNC matrix that accounts for the acceleration noise.
94 The f irst task i s to create the transformation matrix, that relates the acceleration and the measurement data. It is difficult to create a transformation matrix that performs double integration such that the acceleration can be converted to positions and range measurem ents. Instead, the transformation matrix is computed in two steps. First using Equation (4 1), the transformation matrix that relates acceleration and the spherical harmonic coefficients is computed. Once this matrix is created, it is multiplied by the Ja cobian matrix which transfers the state parameters to the measurements. From this point, a cceleration noise covariance (ANC) matrix , is computed as shown in Equation ( 6 1 2 ) and the resulting new MNC matrix , is simply computed by the summation presented in Equation ( 6 1 3 ) = ( ( 6 12 ) ( 6 13 ) Applying both the ARMA Approximation and the ANC In the previous two sections, the ARMA approximation and the ANC w ere explained. Each of these methods are independent approach es to improving the model of measurement covariance matrix. The ARMA approximation consider s the correlation of the LRI noise and the ANC account s for the acceleration noise on the me asurements. These two methods can be applied together Moreover, all generated acceleration noise are colored noise as shown in Figure 3 3. Therefore, the same ARMA approximation can be applied to create the ANC matrix, and this ANC can be added to the ARMA approximated measurement noise as shown in Equation (6 14).
95 ( 6 1 4 ) To verify that the ARMA approximation of the acceleration colored noise is correct the ASD of filtered noise was drawn with the generated noise in Figure 6 2. This figure shows that both the filtered and generated noise s have the same amplitude over the entire spectrum. This means that the ARMA filter was designed well and convert s white noise to the colored noise that is similar to the act ual noise. Figure 6 2 ASD of the GOCE level drag free acceleration generated noise, the ARMA filtered noise, and the noise model. Result Case A2 was recomputed three times in the same manner as explained in Chapter 4 to see how the solution was affected by different modified MNC matrices. At first, the ARMA approximation method was applied. Then, the ANC was added to MNC without the ARMA approximation. Lastly, both the ARMA approximation and the ANC were used to estimate the solution. The n ew
96 sol utions are analyzed and compared with the previous solution from Chapter 5. I t is important to make sure that MNC matrices are always positive definite. Fortunately, all the MNC matrices created were positive definite, and so the estimation was performed w ithout running into any numerical issue s. The GDD error curves for the different solutions are drawn in Figure 6 3. From the figure, it can be seen that geopotential coefficients beyond degree 10, show comparably larger differences than those below degree 10 It is clear that, the ARMA approximated solution has largest error and the ANC applied solution has the smallest error for coefficients greater than degree 10. When both the ARMA and the ANC were applied, the error level was between the ARMA only and the ANC only levels but almost always higher than the previous solution. The result f or coefficients below degree 10 need to be looked more closely to observe which solution is mo re accurate. Figure 6 3 b shows the zoom ed in view of Figure 6 3 a for degre es below 10. It can be seen from this zoom ed in view, that the ARMA approximation improves the solution for degrees less than 10, although it makes the solution worse for coefficients greater than this degree. GS applied GDD errors are shown in Figure 6 3 c. A radius of 300 km was applied to all GS solutions. Since the radius was 300 km, GS does not affect solutions below degree 20 and therefore the conclusion for low degree terms remain s constant. Even for higher degrees, the conclusion that the ARMA approximation increases error, is still applicable to filtered solutions. When GS is applied, the ANC error is almost identical to the initial error.
97 ( a ) ( b ) Figure 6 3 Errors of the ARMA and the ANC applied single pair solutio ns for Case A2 expressed as GDD; (a) raw errors, (b) raw errors for low degrees, and (c) GS applied errors.
98 ( c ) Figure 6 3 Continued. In addition to the GDD error curves, EW H plots are d rawn for each estimation method and they are shown in Figures 6 4 and 6 5. The corresponding spatial RMS values are shown in Table 6 1. Figure 6 4 shows different EW H plots that represent raw estimated solutions and Figure 6 5 include EW H plots that are filtered by GS with a 300 km radi us. For both figures, the top left panel is the initial solution, the top right panel is the ARMA approximated solution, the bottom left panel is the ANC applied solution and the bottom right panel is the both ARMA approximated and ANC applied solution. From the figures, it can be clearly observed that when ARMA approximation is applied with or without ANC, more error stripes appear. Other than this, it is difficult to tell the difference between th e plots. When GS is applied, all four solutions look alike. From Table 6 1, it is clear that the ARMA approximation
99 increases the error by ~1 cm compared to the initial solution. On the other hand, the ANC improves the initial solution by 1~2 mm. When both methods are used at the same time, the error increased but not as much as the ARMA case. When GS is applied to the solutions, the magnitude of the changes are decreased but the general trends remain constant. Figure 6 4 Estimated single pair solutions for Case A2 in units of cm of EW H with various MNC matrices. Initial (top left), ARMA approximation (top right), ANC (bottom left) and ARMA approximation & ANC (bottom right).
100 Figure 6 5 GS applied single pair solutions for Case A2 in units of cm of EW H with various MNC matrices. Initial (top left), ARMA approximation (top right), ANC (bottom left) and ARMA approximation & ANC (bottom right). Table 6 1 Spatial RMS error for Case A2 in uni ts of cm of EW H with various MNC matrices. HI Error (True: 12.44) Total Error (True: 7.85) Max. Error Raw GS Raw GS Raw GS Initial 5.47 3.59 4.45 2.11 82.93 52.70 ARMA 6.49 3.72 5.34 2.26 115.81 52.39 ANC 5.31 3.55 4.27 2.08 82.39 52.36 ARMA & ANC 5.63 3.58 4.61 2.13 115.2 51.72
101 Applying the ARMA Approximation and the ANC for Higher A cceleration Noise Level Mission In the previous section, the result of applying the ANC showed small improvement s for the case of using the GOCE level residual acceleration noise (Case A2). We have seen from Chapter 5 that at this noise level, the residual acceleration noise probably has a small effect compared to other error sources since decreasing the noise beyond this level did not change the solution. Howev er, applying the ANC still improved the solution, which indicates that for higher acceleration noise levels the ANC will be more effective. Therefore, in this section the estimation i s repeated for Case A1 with the ARMA approximation and the ANC. It can be expected that the ARMA approximation will not have large impact on the solution since the overall error is large for Case A1. However, a big improvement is expected by applying the ANC. The r esults are shown in a similar manner as was done previously starting with the GDD errors in Figure 6 6. It is very clear from the figure that the ANC improved the general solution. In fact, beyond degree 13, the ANC based errors are always smaller than the initial solution. However, for low degrees, it is difficult to tell which solution is the best one. Figure 6 6 b shows a close r view near low degree coefficients and the ARMA and initial solution s are better than the ANC solution. This is consistent with the result for Case A2. Nevertheless, unlike Case A2, the ARMA approximation does not improve the solution for low degree coefficients. The GS applied GDD errors are shown in Figure 6 6 c. The ANC solutions are still distinctively better than the other two solutions for degrees beyond 10. For degrees below 10, the ARMA approximation seem s to be slightly better than the initial solution.
102 ( a ) ( b ) Figure 6 6 Errors of the ARMA and the ANC applied single pair solutions for Case A1 expressed as GDD ; (a) raw errors, (b) raw errors for low degrees, and (c) GS applied errors.
103 ( c ) Figure 6 6 Continued. The EW H color maps for the different methods of Case A1 are shown in Figure 6 7. The f irst row shows the initial solutions, the second row is the ARMA approximated solutions and the third row shows the ANC applied solutions. The s olutions on the left column are the raw and the right column are the GS applied solutions From this figure, it is somewhat difficult to notice the difference between the initial solutions and the ARMA approximated solutions. Nonetheless, it is readily seen that the ANC applied solutions have less error stripes for both the filtered and the non filtered solutions. T able 6 2 shows the corresponding spatial RMS errors. These numerical values show that the ARMA approximation increased the error about 1 cm, whereas the ANC decreased the error by almost 4 cm. Even when GS is applied, the ANC showed 1 cm improvement.
104 Figure 6 7 Estimated single pair solutions for Case A1 in units of cm of EW H with various MNC matrices. Initial (row 1), ARMA approximated (row 2) and ANC (row 3), non filtered (left), and filtered (right) solutions.
105 Table 6 2 Spatial RMS error for Case A1 in units of cm of EW H with various MNC matrices. HI Error (True: 12.44) Total Error (True: 7.85) Max. Error Raw GS Raw GS Raw GS Initial 11.72 4.50 11.85 3.30 83.15 55.96 ARMA 12.91 4.71 12.85 3.51 106.58 55.97 ANC 7.78 3.71 8.11 2.52 46.48 54.67 Discussion In this chapter, the estimation routine was modified by removing some of the simplifying assumptions that were made regarding the measurement noise in the previous chapter. These assumptions were (a) having no correlation on the measurement noise and (b) no effect s of acceleration noise accounted for the measurement noise. By applying the ARMA approximation, the measurement noise was considered to be correlated and by adding the ANC to the MNC, the eff ect of acceleration nose on the measurement data was included. This was first applied to Case A2, where the residual acceleration noise is equivalent to the GOCE drag free system. The result showed that the ANC generally improved the solution by 1~2 mm (3% ) and the ARMA increased the overall error by ~1 cm (18%) However, when the errors were observed closely, it was seen that the ARMA improved the solution for low degree coefficients (less than 10) and ANC improve the solution for high degree coefficients (greater than 10). Th ese modified estimation methods were also tried for Case A1, GRACE level residual acceleration noise. As expected from the result s of Case A2, the ANC improved the initial solution dramatically by 3 ~4 cm (30%) and ARMA approximation ag ain increased the error by ~1 cm (8.5%) Even for Case A1, the ANC improved the
106 solution for high degree coefficients (greater than 10) but the ARMA showed improvement in low degree coefficients (less than 10) only for the GS applied solution. However, fo r both cases, it was clear that the ARMA approximated solutions were better than the ANC solutions at low degrees and vice versa for high degrees. When the ARMA approximation was applied, the knowledge on the correlation of the noise informed the estimato r that the noise actually ha s larger errors for lower spatial frequencies. This additional information about the low frequency noise allowed the estimator to estimate the low degree coefficients more accurately although the improvement was small. The resu lts for Case A1 showed a larger improvement compared to Case A2 when the ANC was applied This is because based on the result from Chapter 5, it is clear that Case A1 was dominated by the acceleration noise. In order to make the best use of the ARMA approximation and the ANC, the solutions from the two methods were combined. The GDD error curves showed that the ARMA approximation improved the solution for degrees below 10, whereas the ANC decreased the error for degrees above 10 for both cases. T herefore, the results from the ARMA approximation and the ANC were combined to see how much the solutions are improved. For up to degree 10, the ARMA approximation solution and for higher degrees, the ANC solution were taken. The EW H for th ese combined sol utions were computed and the corresponding spatial RMS is shown in Table 6 3. The GDD error curves and EW H color maps are neglected here According to the table, the combined solution has the minimum error for Case A2, even smaller than the ANC only
107 soluti on but the difference is small enough to be negligible. For Case A1, the combined solution did not show any improvements over the ANC solution. Hence, t he ARMA approximation is not necessary for gravity estimation, but acceleration noise should be conside red in the MNC, especially when the residual acceleration noise is large. Table 6 3 Spatial RMS error for Cases A1 and A2 in units of cm of EW H HI Error (True: 12.44) Total Error (True: 7.85) Max. Error Raw GS Raw GS Raw GS Case A1 Initial 11.72 4.50 11.85 3.30 83.15 55.96 A NC 7.78 3.71 8.11 2.52 46.48 54.67 Combined 7.78 3.71 8.11 2.53 4 6 .73 54.82 Case A2 Initial 5.47 3.59 4.45 2.11 82.93 52.70 A NC 5.31 3.55 4.27 2.08 82.39 52.36 Combined 5.30 3.55 4.26 2.07 82.42 52.34
108 CHAPTER 7 TOOLBOX VALIDATION USING THE UF PRECISION TORSION PENDULUM Although instrument noise can be created with given models and be used in simulations as demonstrated in previous chapters, this simulated noise is only a representation of the actual noise. Measu rement noise from actual instruments are often far richer and more complex than the noise models used in previous chapters. Therefore, it is desirable to test the toolbox and run the simulation with the actual noise from the existing instruments. The r esid ual acceleration noise of a drag free system can be obtained from UF precision torsion pendulum which simulates free fall of a drag free test mass in the laboratory environment This chapter introduces the UF precision torsion pendulum and compare s the gr avity solutions estimated using noise data taken from the pendulum and the simulated noise UF Precision Torsion Pendulum Cooperative work between the P hysics and A erospace E ngineering department s at the University of Florida ( UF ) is currently developing a facility that is designed to test new technologies for the upcoming Laser Interferometer Space Antenna (LISA) mission, a dedicated space mission aimed to detect gravitational wave s in space [ Ciani et al. 2017]. Figure 7 1 shows the CAD model and the ac tual apparatus of the torsion pendulum. As seen from the figure, the system is placed inside the vacuum chamber and the entire system is enclosed in a thermal room to mitigate disturbances due to air currents and thermal fluctuations
109 CAD model Actua l Apparatus Figure 7 1 CAD model and actual apparatus of a torsion pendulum taken from Physics Space Laboratory at University of Florida in 2015 [ C ia ni et al. 2017 ]. UF precision torsion pendulum measure s the residual acceleration noise of a proof mass in one degree of freedom by measuring the angle variations along the vertical axis. Figur e 7 2 shows the inertial member of the pendulum, which are held by a tungsten fiber. As seen from this figure, two displacements and are measured using both ca pacitive sensor and laser interferometer, which are then used to compute the angle variation, This angle measurement is first converted to a torque measurement and then to an acceleration measurement
110 Figure 7 2 Inertial members of a torsion pendul um taken from Physics Space Laboratory at University of Florida in 2017 [ C ia ni et al. 2017]. Accuracy of the acceleration measurement is aimed to be comparable to the LISA requirement which is [ Chilton et al. 2015] ( 7 1 ) A g ravitational reference sensor (GRS) is a key technology of a drag free system that measures the displacement of a test mass relative to its housing. Therefore, the noise measurement from UF precision torsion pendulum can be used as the actual acceleration noise of a drag free system. Moreover, LISA is focusing on frequencies between 0.1 mHz and 100 mHz, the same frequency band that satellite g eodesy missions are interested in. Accumulated Acceleration Noise The acceleration noise measured by the UF precision torsion pendulum in November and December 2016 were accumulated and rearranged to be used in this research. Figure 7 3 shows the amplitud e spectral density (ASD) curves of the accumulated noise. The noise data is in 3 hour segments with a sampling
111 frequency of 8.6 Hz. As seen from the figure, there is a peak caused by the actual swing mode of the pendulum at frequency near 0.4 Hz. This kind of error caused by the pendulum natural dynamics does not exist in an actual drag free system. Reg ardless of this peak error, the errors above around 3 mHz are greatly dominated by the measurement system of the pendulum and below this frequency, the noise is dominated by the acceleration noise. Therefore, in order to extract only the acceleration noise, the low pass filter must be applied at frequency of 3 mHz and then appropriate acceleration noise must be added for frequency above 3 mHz. Acceleration noi se models for drag free systems shown in Figure 3 3, tend to decrease and level off at higher frequencies. Accordingly, white noise was added only for the frequencies higher than 3 mHz after applying the low pass filter Therefore, this added white noise d oes not alter the noise at low frequencies. Figure 7 3 ASD of raw acceleration noise data
112 In addition, the sampling frequency of the measurements in this research is 0.1 Hz, which is smaller than the sampling frequency of the pendulum. The accumulated acceleration data therefore needs to be reorganized such that the noise data is provided at sl ower rate. This can be done simply by taking the average of the filtered data for time interval of 10 sec. This new extracted acceleration noise from the accumulated data at 0.1 Hz are shown in Figure 7 4. Figure 7 4 ASD of extracted acceleration noise data. To be able to compare the actual noise with the model based simulated noise, a simulated noise was created using the method presented in Chapter 3. Figure 7 5 shows ASD curves of actual acceleration noise, simulated acceleration noise and the model used to generate the noise. At the bottom of this figure, the model for the pendulum noise is presented. Also, the GOCE level drag free system simulated noise is shown in the same figure. The actual pendulum noise is steeper at low frequency, but other tha n this, the model based simulated noise seems to be a good approximation of the actual pendulum noise.
113 The pendulum noise generally has a similar amplitude as the GOCE level drag free system, but the GOCE noise is slightly lower than the pendulum noise for all frequenc ies At high frequencies, the pendulum noise has the magnitude of 8 10 12 ms 2 Hz 1/2 whereas the GOCE noise has the magnitude of 2 10 12 ms 2 Hz 1/2 Figure 7 5 ASD of actual pendulum noise, model based simulated pendulum noise, and GOCE level simulated acceleration noise. The model of the simulated pendulum noise is shown below the graph. Result In total there are 26 sets of pendulum acceleration noise acquired. These noise data sets were used randomly and evenly to create new measurement data. Also, equal number of simulated acceleration noise sets were created and used to compute measurement data for co mparison. The solutions were computed using the same method as discussed in Chapter 4.
114 Figure 7 6 Errors of the actual and the simulated UF pendulum noise applied single pair solutions expressed as GDD. Raw errors (top) and GS applied errors (bottom).
115 The GDD error curves are shown in Figure 7 6. The raw estimates are shown on top and the filtered solutions are shown on bottom with 300 km radius GS. It is clear from this figure that both cases with actual and generated pendulum noise are somewhat similar. As one would expect, the actual pendulum noise produced a slightly larger error compared to the simulated solutions even when the filter was applied. This is because the actual noise acquired from the pendulum is less consistent compared to the simulated noise In other words, there are more fluctuations in the magnitude of the ASD of the real noise compared to the simulated noise However, the main purpose of applying the actual noise is to validate the toolbox. The GDD error plots clearly show that the estimator produces solutions even when the ASD of the noise is more complex than the simulated noise The solution using the simulated nois e of GOCE is also drawn for comparison. I t is shown that when the acceleration noise is slightly higher than the GOCE level, it worsens the solution. Figure 7 7 Estimated single pair solutions in units of cm of EW H applying the pendulum noise. Actual pendulum noise (row 1), generated pendulum noise (row 2) and generated GOCE noise (row 3), non filtered (left), and filtered (right) solutions.
116 Figure 7 7 Continued. Figure 7 7 shows the solutions in EW H N on filtered solutions are on left and filtered solutions are on right. The f irst row shows the solutions using the actual pendulum noise and right below them are solutions using the simulated pendulum noise. It is quite clear that more errors are shown in the first ro w than the second row. This is consistent with what we have seen from the GDD error plots in Figure 7 6. The last row shows the result using the GOCE simulated acceleration noise, and the discrepanc ies between the GOCE and the pendulum solutions are hard to observe for both raw and GS applied solutions. Table 7 1 show s the numerical values of the spatial RMS EW H for all three solutions plotted in Figure 7 7. The table values tell us that o n average, the actual
117 pendulum noise case has about 2 cm (26%) larger error compared to the generated pendulum noise. This might seem relatively large, but the actual pendulum noise ha s a larger variations in the noise data. For this reason, a larger error is expected when actual noise is used compared to simulated n oise. If the generated noise was created more carefully constructing a more detailed acceleration noise model then the two solutions might look more alike. However, the goal was to check the validity of the toolbox used in this simulation by the actual no ise and this result proves that the toolbox operates well even when acceleration noise taken from actual hardware is used The p endulum simulated noise i s then compared with the GOCE level simulated noise. This comparison shows that the pendulum simulated noise case produce s a solution with an error of 4~5 mm (10%) higher than the GOCE level simulated noise (Case A2). Table 7 1 Spatial RMS error for pendulum cases and Case A2 in units of cm of EW H HI Error (True: 12.44) Total Error (True: 7.85) Max. Error Raw GS Raw GS Raw GS Pendulum: Actual 7.65 3.82 6.36 2.43 80.6 51.87 Pendulum: Simulated 5.64 3.61 4.90 2.18 83.16 52.5 0 GOCE: Simulated 5.47 3.59 4.45 2.11 82.93 52. 70 Discussion In this chapter the UF precision pendulum was introduced and the residual acceleration noise measured by the pendulum was used to validate the toolbox that was developed and used through out this research. The a ccumulated noise from the pendulum had to be modified to be a pplied in this simulation. Since the raw noise data is the combination of the acceleration and the
118 measurement noise, the acceleration noise from the raw noise data had to be ext racted. This was done by applying a low pass filter and adding appropriate hig h frequency noise. Then simulated acceleration noise was generated based on this actual noise to compare th e solutions. The result showed that the solution was successfully estimated with the actual pendulum noise. This means that the toolbox is stable and produces valid results even w hen the actual noise is used Also, the solution was comparable to the solution from the simulated pendulum noise. Although the error magnitude decreased with the simulated noise, if the noise was generated more carefully the two solutions might become more alike. This shows that the generated noise can be a good representation of the actual noise. Finally, by comparing the Case A2 (GOCE acceleration noise) solution with the simulated pendulum noise, we showed that GOCE level acceleration noise is very close to the optimal value. The acceleration noise difference between Case A2 and A3 are about 2 orders of magnitude. However, pendulum acceleration noise was only 4 times larg er than the GOCE acceleration noise and the solution error increased by 4~5 mm (10%) whereas the difference between Case A2 an d A3 solutions is only about 0.1 mm (0.2%) This result confirms that in order to reduce the effect of varying residual accelerat ion noise below 1 mm EWH the noise magnitude must be reduced below 8 10 12 ms 2 Hz 1/2 close to 2 10 12 ms 2 Hz 1/2
119 CHAPTER 8 CONCLUSION S AND RECOMMENDATIONS Conclusions This research accomplished main objectives, which are to build a toolbox for satellite geodesy simulation and to explore the effects of different drag free systems for future satellite geodesy missions. This toolbox was built using MATLAB from the ground up. There are two parts to this tool box. The f irst part creates the measurem ent data of the satellites and the second part uses this measurement to estimate the gravity field solutions. There are six different gravity models; static, ocean tide, ocean bottom pressure, atmosphere, hydrology and ice used in the simulation. The m easu rement noise and the residual acceleration noise are generated using the noise models provided from the various technology teams, and are utilized in the process of creating measurement data. The estimation part was difficult to achieve due to the large am ount of data and parameters. Various methods were applied and the most efficient one in terms of computation time and a convenien ce was the partitioned best linear unbiased estimator (weighted least squares estimator). With this unique toolbox, different levels of residual acceleration noise were applied, each corresponding to an existing technology. The simulation was designed that can repres ent future geodesy mission beyond GRACE FO. In line and in line Bender configurations were both tried in the simulation to evaluate the different levels of aliasi ng errors. Moreover, both realistic and idealistic de aliasing models were applied to discover the relation between the acceleration noise and the aliasing error which are the main concerns for the next generation geodesy missions. Through repeated simulation, it was concluded that the acceleration noise level
120 of the drag free system for future LL SST mission is close to 2 10 12 ms 2 Hz 1/2 and any further reduction of this acceleration noise is not necessary even when aliasing error is reduced dramatically using idealistic de aliasing models and the in line Bender configuration. In addition, t he actual residual acceleration noise from the UF precision torsion pendulum was used to validate the toolbox. The noise from UF precision torsion pendulum not only validated the toolbox, but it also explored an additional level of acceleration noise which was 4 times larger than the GOCE drag free system noise. The solution using the UF pendulum acceleration noise showed that the acceleration noise must be reduced down to GOCE level, in order to reduce the effect of acceleration noise below 1 mm of EW H M oreover, the estimation routine was revisited by removing some of the assumptions regarding the noise. The ARMA approximation and the ANC methods were introduced to modify the measurement noise covariance matrix The ARMA approximation was applied to consi der the correlation of the noise and the ANC was applied to add the effects of acceleration noise on the measurement. It was shown that the ANC improved the EW H results by ~ 3 0% at GRACE level and ~ 3% at GOCE level acceleration noise. However, the ARMA approximation increased the EW H error for both cases. Recommendations for Future Research There are several recommendations for the future research. Updating the toolbox is highly recommended. First, the toolbox can be converted to different programm ing language to increase the speed of computation because MATLAB is inefficient at processing large volumes of data This way, more parameters can
121 be estimated increasing the spatial res olution of the solution and saving computation time. Also, this toolb ox uses spherical harmonic coefficients to represent the gravity field I nstead of spherical harmonic coefficients, mascon parameters can be used [ Save et al. 2016] This new class of parameters are more accurate because each of parameter represent s a sin gle independent grid of the map, whereas one spherical harmonic coefficient effect s the entire globe. Optimizing the arc length should also be explored. Currently different institutions use different arc length However, there is no dedicated study done to find the optimal arc length. It would be valuable to do a research on Cube Sat based geodesy missions Accurate mathematical model s and data from the actual hardware can be applied to this toolbox and the efficiency of Cube Sat s on geodesy missions could be explored. Multiple Cube Sat formations can be considered, and the optimal orbits for multiple CubeSats could be found and compared to the one and two pair GRACE like missions.
122 LIST OF REFERENCES Antonucci, F., Armano, M., Audley, H., Auger, G., Benedetti, M., Binetruy, P., ... & Caleno, M. (2012). The LISA pathfinder mission. Classical and Quantum Gravity 29 (12), 124014. Armano, M., Audley, H., Auger, G., Baird, J., Binetruy, P., Born, M., ... & Cavalleri, A. (2015). The LISA pathfinder mission. In Journal of Physics: Conference Series (Vol. 610, No. 1, p. 012005). IOP Publishing. Armano, M., Audley, H., Auger, G., Baird, J. T., Bassan, M., Binetruy, P., ... & Carbone, L. (2016). Sub femto g free fall for space based gravitational wave observatories: LISA pathfinder results. Physical Review Letters 116 (23), 231101 Barthelmes, F. (2009). Definition of functionals of the geopotential and their calculation from spherical harmonic models. http://publi cations. iass potsdam. de/pubman/item/escidoc 104132 (3), 0902 2. Bettadpur, S. (2012). UTCSR Level 2 processing standards document for Level 2 product release 5. Bruinsma, S. L., Frste, C., Abrikosov, O., Marty, J. C., Rio, M. H., Mulet, S., & Bonvalot, S. (2013). The new ESA satellite only gravity field model via the direct approach. Geophysical Research Letters 40 (14), 3607 3612. Canuto, E., & Massotti, L. (2009). All propulsion design of the drag free and attitude control of the European satellite GOC E. Acta Astronautica 64 (2), 325 344. Chen, Q., Shen, Y., Chen, W., Zhang, X., & Hsu, H. (2016). An improved GRACE monthly gravity field solution by modeling the non conservative acceleration and attitude observation errors. Journal of geodesy 90 (6), 503 523. Chilton, A., Shelley, R., Olatunde, T., Ciani, G., Conklin, J. W., & Mueller, G. (2015). The UF Torsion Pendulum, a LISA Technology Testbed: Sensing System and Initial Results. In Journal of Physics: Conference Series (Vol. 610, No. 1, p. 012038). IOP Publishing. Ciani, G., Aitken, M., Apple, S., Chilton, A., Olatunde, T., Mueller, G., & Conklin, J. W. (2017). A New Torsion Pendulum for Gravitational Reference Sensor Technology Development. arXiv preprint arXiv:1701.08911 Conklin, J. W. (2015). Drift mode accelerometry for spaceborne gravity measurements. Journal of Geodesy 89 (11), 1053 1070. Dahle, C. (2013). GFZ Level 2 processing standards document for Level 2 product release 5
123 DeBra, D. B., & Conklin, J. W. (2011). Measurement of drag and its cancellation. Classical and Quantum Gravity 28 (9), 094015. Dobslaw, H., Flechtner, F., Bergmann Wolf, I., Dahle, C., Dill, R., Esselborn, S., ... & Thomas, M. (2013). Simulating high frequency atmosphere ocean mass variability for dealiasing of satellite gravity observations: AOD1B RL05. Journal of Geophysical Research: Oceans 118 (7), 3704 3711. Dobslaw, H., Bergmann Wolf, I., Dill, R., Forootan, E., Klemann, V., Kusche, J., & Sasgen, I. (2015). The upda ted ESA Earth System Model for future gravity mission simulation studies. Journal of Geodesy 89 (5), 505 513. Dolesi, R., Bortoluzzi, D., Bosetti, P., Carbone, L., Cavalleri, A., Cristofolini, I., ... & Hoyle, C. D. (2003). Gravitational sensor for LISA an d its technology demonstration mission. Classical and quantum gravity 20 (10), S99. Elsaka, B., Raimondo, J. C., Brieden, P., Reubelt, T., Kusche, J., Flechtn er, F., ... & Mller, J. (2013 ). Comparing seven candidate mission configurations for temporal gra vity field retrieval through full scale numerical simulation. Journal of Geodesy 88 (1), 31 43. Elsaka, B. (2014). Sub monthly gravity field recovery from simulated multi GRACE mission type. Acta Geophysica 62 (1), 241 258. Flechtner, F. (2007). AOD1B Prod uct Description Document for Product Releases 01 to 04 (Rev. 3.1, April 13, 2007). GRACE project document 327 750. Flechtner, F., Dobslaw, H., & Fagiolini, E. (2014). AOD1B Product Description Document for Product Release 05 (Rev. 4.2, May 20, 2014). Tech. Rep. GFZ German Research Centre for Geosciences, Potsdam Flechtner, F., Morton, P., Watkins, M., & Webb, F. (2014). Status of the GRACE follow on mission. In Gravity, geoid and height systems (pp. 117 121). Springer International Publishing. Flecht ner, F., Neumayer, K. H., Dahle, C., Dobslaw, H., Fagiolini, E., Raimondo, J. C., & Gntner, A. (2016). What can be expected from the GRACE FO Laser Ranging Interferometer for Earth Science applications?. In Remote Sensing and Water Resources (pp. 263 280) Springer International Publishing. Floberghagen, R., Fehringer, M., Lamarre, D., Muzi, D., Frommknecht, B., Steiger, C., ... & Da Costa, A. (2011). Mission design, operation and exploitation of the gravity field and steady state ocean circulation explore r mission. Journal of Geodesy 85 (11), 749 758. Flury, J., Bettadpur, S., & Tapley, B. D. (2008). Precise accelerometry onboard the GRACE gravity field satellite mission. Advances in Space Research 42 (8), 1414 1423.
124 Foerste, C., Flechtner, F., Schmidt, R., Stubenvoll, R., Rothacher, M., Kusche, J., ... & Bruinsma, S. (2008, April). EIGEN GL05C A new global combined high resolution GRACE based gravity field model of the GFZ GRGS cooperation. General Assembly European G eosciences Union (Vienna, Austria 2008), Geophysical Research Abstracts 10 (EGU2008 A 06944). Gruber, T., Bamber, J. L., Bierkens, M. F. P., Dobslaw, H., Murbck, M., Thomas, M., ... & Visser, P. N. A. M. (2011). Simulation of the time variable gravity field by means of coupled geophysical models. Earth System Science Data 3 (1), 19 35. Heinzel, G., Sheard, B., Brause, N., Danzmann, K., Dehne, M., Gerberding, O., ... & Klipstein, W. (20 12). Laser ranging interferometer for GRACE follow on. In Proceeding of International Conference on Space Optics (ICSO)(Ajaccio, Corsica, France, 2012) Heiskanen, W. A., & Moritz, H. (1967). Physical geodesy. Bulletin Godsique (1946 1975) 86 (1), 491 49 2. Kang, Z., Tapley, B., Bettadpur, S., Ries, J., Nagel, P., & Pastor, R. (2006). Precise orbit determination for the GRACE mission using only GPS data. Journal of Geodesy 80 (6), 322 331. Kasdin, N. J. (1995). Discrete simulation of colored noise and stoc hastic processes and 1/f/sup/spl alpha//power law noise generation. Proceedings of the IEEE 83 (5), 802 827. Kim, J. (2000). Simulation study of a low low satellite to satellite tracking mission (Doctoral dissertation, The University of Texas at Austin). K lees, R., Ditmar, P., & Broersen, P. (2003). How to handle colored observation noise in large least squares problems. Journal of Geodesy 76 (11 12), 629 640. Klinger, B., & Mayer Grr, T. (2016). The role of accelerometer data calibration within GRACE gravity field recovery: Results from ITSG Grace2016. Advances in Space Research 58 (9), 1597 1609. Kurtenbach, E., Mayer G rr, T., & Eicker, A. (2009). Deriving daily snapshots of the Earth's gravity field from GRACE L1B data using Kalman filtering. Geophysical Research Letters 36 (17). Loomis, B. D., Nerem, R. S., & Luthcke, S. B. (2012). Simulation study of a follow on gravity mission to GRACE. Journal of Geodesy 86 (5), 319 3 35. Lyard, F., Lefevre, F., Letellier, T., & Francis, O. (2006). Modelling the global ocean tides: modern insights from FES2004. Ocean Dynamics 56 (5 6), 394 415. M ontenbruck, O., & Gill, E. (2005 ). Satellite orbits Models, Methods and Applications (3 rd Edition) New York: Springer
125 Moore, P., Zhang, Q., & Alothman, A. (2005). Annual and semiannual variations of the Earth's gravitational field from satellite laser ranging and CHAMP. Journal of Geophysical Research: Solid Earth 110 (B6). Nerem, R. S., Jeke li, C., & Kaula, W. M. (1995). Gravity field determination and characteristics: retrospective and prospective. Journal of Geophysical Research: Solid Earth 100 (B8), 15053 15074. Nguyen, A. (2016). Small satellite technologies for precision drag free contr ol and drag force recovery (Doctoral dissertation, The University of Florida). Nguyen, A. & Conklin, J. W. (2015). Three Axis Drag Free Control and Drag Force Recovery of a Single Thruster Small Satellite. Journal of Spacecraft and Rockets 52 (6), 1640 16 50. Pail, R., Bingham, R., Braitenberg, C., Dobslaw, H., Eicker, A., Gntner, A., ... & Wouters, B. (2015). Science and user needs for observing global mass transport to understand global change and to benefit society. Surveys in geophysics 36 (6), 743 772. time correlated errors in a Kalman filter applicable to GNSS. Journal of Geodesy 83 (1), 51 56. Reigber, C., Schwintzer, P., Neumayer, K. H., Barthelmes, F., K nig, R., Frste, C., ... & Bruinsma, S. (2003). The CHAMP only Earth gravity field model EIGEN 2. Advances in Space Research 31 (8), 1883 1888. Rieser, D., Mayer Grr, T., Savcenko, R., Bosch, W., Wnsch, J., Dahle, C., & Flechtner, F. (2012). The ocean t ide model EOT11a in spherical harmonics representation. URL: http:// portal. tugraz. at/ portal/ page/ portal/ Files/ i5210/ files/ projekte/ COTAGA/ TN_ EOT11a. pdf Sakumura, C., Bettadpur, S., Save, H., & McCullough, C. (2016). High freque ncy terrestrial water storage signal capture via a regularized sliding window mascon product from GRACE. Journal of Geophysical Research: Solid Earth 121 (5), 4014 4030. Savcenko, R., & Bosch, W. (2012). EOT11a empirical ocean tide model from multi mission satellite altimetry. Save, H., Bettadpur, S., & Tapley, B. D. (2016). High resolution CSR GRACE RL05 mascons. Journal of Geophysical Research: Solid Earth 121 (10), 7547 7569. Seeber, G. (2003). Satellite geodesy: foundations, methods, and applications Walter de Gruyter.
126 Sheard, B. S., Heinzel, G., Danzmann, K., Shaddock, D. A., Klipstein, W. M., & Folkner, W. M. (2012). Intersatellite laser ranging instrument for the GRACE follow on mission. Journal of Geodesy 86 (12), 1083 1095. Sneeuw, N., Flury, J., & Rummel, R. (2005). Science requirements on future missions and simulated mission scenarios. In Future Satellite Gravimetry and Earth Dynamics (pp. 113 142). Springer New York. Stoica, P., & Moses, R. L. (1997). Introduction to spectral analysis (Vol. 1, pp. 3 4). Upper Saddle River: Prentice hall. effect on solutions of differential equations. International journal for uncertainty quantification 1 (3). Swenson, S., & Wahr, J. (2002). Methods for inferring regional surface mass anomalies from Gravity Recovery and Climate Experiment (GRACE) measurements of time variable gravity. Journal of Geophysical Research: Solid Earth 107 (B9). Tapley, B. D., Bettadpur, S., Watkins, M., & Reigber, C. (2004). The gravity recovery and climate experiment: Mission overview and early results. Geophysical Research Letters 31 (9). Tapley, B. D. Ries, J., Bettadpur, S., Chambers, D., Cheng, M., Condi, F., & Poole, S. (2007, December). The G GM03 mean earth gravity model from GRACE. In AGU Fall Meeting Abstracts (Vol. 1, p. 03). Touboul, P., Foulon, B., Christophe, B., Marque, J.P., (2012). CHAMP, GRACE, GOCE instruments and beyond. Geodesy for planet earth In: Kenyon, S., et al. (Eds.), IAGS, vol. 136, pp. 215 221. Vtipil, S., & Newman, B. (2012). Determining an Earth Observation Repeat Ground Track Orbit for an Optimization Methodology. Journal of Spacecraft and Rockets 49 (1), 157 164. Wahr, J., Molenaar, M., & Bryan, F. (1998). Time va riability of the Earth's gravity field: Hydrological and oceanic effects and their possible detection using GRACE. Journal of Geophysical Research: Solid Earth 103 (B12), 30205 30229. Watkins, M. M., & Yuan, D. N. (2012). JPL Level 2 processing standards d ocument for Level 2 product release 5. Wiese, D. N. (2011). Optimizing two pairs of GRACE like satellites for recovering temporal gravity variations (Doctoral dissertation, The University of Colorado). Wiese, D. N., Nerem, R. S., & Lemoine, F. G. (2012). Design considerations for a dedicated gravity recovery satellite mission consisting of two pairs of satellites. Journal of Geodesy 86 (2), 81 98.
127 BIOGRAPHICAL SKETCH Seong Hyeon Hong received h degree in mechanical and electronic control e ngineering as a dual degree in Feb. 2012. He did his entire undergraduate at the Handong Global University (HGU) starting from Mar. 2005. Between his undergraduate years, he served in the Republic of Korea milit ary for two years from 2007 to 2009. At HGU, he was a member of robotics group and participated in the national robot competition with unmanned vehicle s He worked on building a two wheeled self balancing vehicle with a linear quadratic controller for his thesis. Seong Hyeon began his graduate study from Aug. 2012 at University of Florida. He was awarded with Master of Science (M.S.) degree in Dec. 2013 in mechanical e ngineering. He continued his graduate study at University of Florida in a Ph.D. program. During his time at University of Florida, he worked as a research assistant in Precision Space Systems Laboratory from Jan. 2013 to Apr. 2017. His research field was satellite geodesy, specifically work ing on the effects of drag free system acceleration no ise levels for future low low satellite to satellite tracking missions. For this research, he built a new toolbox that can be used to simulate satellite ge odesy missions and estimate gravity field solutions. He also served as a teaching assistant for Mecha nical and Aerospace Engineering D epartment from Jan. 2013 to Apr. 2017. He served as a mentor for a student with a learning disability with the same major from Aug. 2014 to Apr. 2015. In the summer 2016, he did internship at Korea Helicopter Company in Sou th Korea. He worked on building commercial quadrotors, which spray agricultural chemicals.
128 Seong Hyeon will be awarded with Doctor of Philosophy (Ph.D.) degree as a m echanical e ngineering major in May 2017 at University of Florida. He will continue to wor k on advanced research with experience he gained at University of Florida.