Performance Comparison of an Extended Kalman Filter and an Iterated Extended Kalman Filter for Orbit Determination of Space Debris with Poor Apriori Information and Intermittent Observations by Jimmy Dale Hicks Jr. A thesis submitted to the Graduate Faculty of Auburn University in partial ful llment of the requirements for the Degree of Master of Science Auburn, Alabama December 8, 2012 Keywords: orbit debris, state estimation, kalman ltering Copyright 2012 by Jimmy Dale Hicks Jr. Approved by David Cicci, Chair, Professor of Aerospace Engineering John E. Cochran Jr., Professor and Head of Aerospace Engineering John Hung, Professor of Electrical and Computer Engineering Abstract The growing space debris population threatens active and future missions bound for low earth orbit. The purpose of this study is to determine if an iterated extended Kalman lter (IEKF) can be used to yield better performance than the extended Kalman lter (EKF) in the orbit determination of space debris in low earth orbit with poor apriori data and intermittent observations. A simulation using two-body and J2 e ects was constructed using actual Space Surveillance Network sensor locations to generate experimental observational data, which contains multiple data outages. This data was then used to compare the performance of the EKF and IEKF. The lters are compared using the di erence in the lter estimate and the true state and the di erence in the estimated observation and true observation. An increase in estimate accuracy will allow for better predictions concerning the interaction of space debris with missions in low earth orbit. The IEKF was chosen for this study due to the low number of observations provided and the large apriori covariance matrix. The state update step in the IEKF includes a local iteration that processes each observation until convergence of the state update is reached, which becomes the new best estimate of the state. Following the local convergence, the covariance matrix is updated using this new estimate, which prevents the covariance matrix from growing small as quickly as that of the EKF. ii Acknowledgments I would like to thank Dr. David Cicci for the opportunity to work on this research as well as his knowledge, invaluable support, and time in assisting in the development of this thesis. Our conversations have lead to a deeper understanding of orbit determination, and also the broader eld of state estimation, which I greatly enjoy and have established as the foundation for my future career. I would also like to thank the Aerospace Engineering Department for its nancial assis- tance and Dr. John Cochran for providing a position in his lab where my interest in precise target tracking grew. I would like to thank Dr. John Hung for the knowledge of what graduate school really meant. My conversations with him concerning his class led me to broader understanding of the types of planning and research needed to succeed. iii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 History . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Previous Research . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.3 Problem Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2 Relevant Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.1 Reference Frames . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.1.1 Earth Centered Inertial Coordinate System (ECI) . . . . . . . . . . . 7 2.1.2 Earth-Centered, Earth-Fixed Coordinate System (ECEF) . . . . . . . 8 2.1.3 Perifocal Coordinate System . . . . . . . . . . . . . . . . . . . . . . . 9 2.2 Two-body Equations of Motion . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.3 J2 Perturbation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.4 Orbital Plane Orientation in Space . . . . . . . . . . . . . . . . . . . . . . . 13 2.5 Coordinate Transformations . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.5.1 Obtain ECI position and velocity vectors from the Keplerian orbital elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.5.2 Calculate the classical Keplerian orbital elements from ECI position and velocity vectors . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.5.3 ECI to ECEF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3 State Estimation and Kalman Filtering . . . . . . . . . . . . . . . . . . . . . . . 20 iv 3.1 Linearization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.2 Kalman Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.2.1 Extended Kalman Filter . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.2.2 Iterated Extended Kalman Filter . . . . . . . . . . . . . . . . . . . . 24 4 Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 4.1 Baseline Initial Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 4.2 Test Cases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 4.3 Constants . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 4.4 Dynamical Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 4.5 Measurement Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 4.6 Tracking Station Network . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.7 Propagator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 4.8 Procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 5.1 Measure of Solution Accuracy . . . . . . . . . . . . . . . . . . . . . . . . . . 35 5.2 Baseline Case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 5.3 Inclination Variation Case !i = 75:0 . . . . . . . . . . . . . . . . . . . . . 43 5.4 Argument of Perigee Case !! = 90 . . . . . . . . . . . . . . . . . . . . . . 51 5.5 Additional Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 6 Conclusions and Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 A Derivation of J2 Equations of motion . . . . . . . . . . . . . . . . . . . . . . . . 64 B Additional Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 B.1 Variations on Inclination . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 B.1.1 Inclination !i = 28:5 . . . . . . . . . . . . . . . . . . . . . . . . . . 69 B.2 Variations on Altitude . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 v B.2.1 Altitudes ! rpra = 800km=1000km . . . . . . . . . . . . . . . . . . . . 73 B.2.2 Altitudes ! rpra = 500km=700km . . . . . . . . . . . . . . . . . . . . 77 B.3 Variations on Longitude of the Ascending Node . . . . . . . . . . . . . . . . 81 B.3.1 Longitude of the Ascending Node ! = 0 . . . . . . . . . . . . . . 81 B.3.2 Longitude of the Ascending Node ! = 90 . . . . . . . . . . . . . 85 B.4 Variations on the Argument of Perigee . . . . . . . . . . . . . . . . . . . . . 89 B.4.1 Argument of Perigee !! = 45 . . . . . . . . . . . . . . . . . . . . . 89 B.5 Variations on the True Anomaly . . . . . . . . . . . . . . . . . . . . . . . . . 93 B.5.1 True Anomaly ! = 90 . . . . . . . . . . . . . . . . . . . . . . . . 93 B.5.2 True Anomaly ! = 270 . . . . . . . . . . . . . . . . . . . . . . . . 97 vi List of Figures 1.1 Debris Environment Trends . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.1 ECI coordinate system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.2 ECEF coordinate system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.3 Perifocal coordinate system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.4 Two-body problem in a relative coordinate system . . . . . . . . . . . . . . . . 11 2.5 Orbit in Spherical Coordinates . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.6 Orbital Elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 5.1 Baseline: Residuals vs. Number of Observations Processed . . . . . . . . . . . . 37 5.2 Baseline: Residuals for Batch #6 . . . . . . . . . . . . . . . . . . . . . . . . . . 38 5.3 Baseline: Residuals for Batch #7 . . . . . . . . . . . . . . . . . . . . . . . . . . 38 5.4 Baseline: vs. Number of Observations Processed . . . . . . . . . . . . . . . . 39 5.5 Baseline: NMSE for Batch #6 and Batch #7 . . . . . . . . . . . . . . . . . . . 40 5.6 Baseline: NMSE for Batch #6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 5.7 Baseline: NMSE for Batch #7 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 5.8 i = 75 : Residuals vs. Number of Observations Processed . . . . . . . . . . . . 45 vii 5.9 i = 75 : Residuals for Batch #7 . . . . . . . . . . . . . . . . . . . . . . . . . . 46 5.10 i = 75 : Residuals for Batch #8 . . . . . . . . . . . . . . . . . . . . . . . . . . 46 5.11 i = 75 : vs. Number of Observations Processed . . . . . . . . . . . . . . . . . 47 5.12 i = 75 : NMSE for Batch #7 and Batch #8 . . . . . . . . . . . . . . . . . . . . 48 5.13 i = 75 : NMSE for Batch #7 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 5.14 i = 75 : NMSE for Batch #8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 5.15 ! = 90 : Residuals vs. Number of Observations Processed . . . . . . . . . . . . 52 5.16 ! = 90 : Residuals for Batch #6 . . . . . . . . . . . . . . . . . . . . . . . . . . 53 5.17 ! = 90 : Residuals for Batch #7 . . . . . . . . . . . . . . . . . . . . . . . . . . 53 5.18 ! = 90 : vs. Number of Observations Processed . . . . . . . . . . . . . . . . 54 5.19 ! = 90 : NMSE for Batch #6 and Batch #7 . . . . . . . . . . . . . . . . . . . 55 5.20 ! = 90 : NMSE for Batch #6 . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 5.21 ! = 90 : NMSE for Batch #67 . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 B.1 i = 28:5 : vs. Number of Observations Processed . . . . . . . . . . . . . . . . 69 B.2 i = 28:5 : Residuals vs. Number of Observations Processed . . . . . . . . . . . 70 B.3 i = 28:5 : Residuals Before Largest Outage vs. Time . . . . . . . . . . . . . . . 70 B.4 i = 28:5 : Residuals After Largest Outage vs. Time . . . . . . . . . . . . . . . 71 B.5 i = 28:5 : NMSE Before Largest Outage vs. Time . . . . . . . . . . . . . . . . 71 viii B.6 i = 28:5 : NMSE After Largest Outage vs. Time . . . . . . . . . . . . . . . . . 72 B.7 rpra = 800km=1000km : vs. Number of Observations Processed . . . . . . . . . 73 B.8 rpra = 800km=1000km : Residuals vs. Number of Observations Processed . . . . 74 B.9 rpra = 800km=1000km : Residuals Before Largest Outage vs. Time . . . . . . . . 74 B.10 rpra = 800km=1000km : Residuals After Largest Outage vs. Time . . . . . . . . . 75 B.11 rpra = 800km=1000km : NMSE Before Largest Outage vs. Time . . . . . . . . . . 75 B.12 rpra = 800km=1000km : NMSE After Largest Outage vs. Time . . . . . . . . . . 76 B.13 rpra = 500km=700km : vs. Number of Observations Processed . . . . . . . . . . 77 B.14 rpra = 500km=700km : Residuals vs. Number of Observations Processed . . . . . 78 B.15 rpra = 500km=700km : Residuals Before Largest Outage vs. Time . . . . . . . . . 78 B.16 rpra = 500km=700km : Residuals After Largest Outage vs. Time . . . . . . . . . 79 B.17 rpra = 500km=700km : NMSE Before Largest Outage vs. Time . . . . . . . . . . 79 B.18 rpra = 500km=700km : NMSE After Largest Outage vs. Time . . . . . . . . . . . 80 B.19 = 0 : vs. Number of Observations Processed . . . . . . . . . . . . . . . . . 81 B.20 = 0 : Residuals vs. Number of Observations Processed . . . . . . . . . . . . 82 B.21 = 0 : Residuals Before Largest Outage vs. Time . . . . . . . . . . . . . . . . 82 B.22 = 0 : Residuals After Largest Outage vs. Time . . . . . . . . . . . . . . . . 83 B.23 = 0 : NMSE Before Largest Outage vs. Time . . . . . . . . . . . . . . . . . 83 ix B.24 = 0 : NMSE After Largest Outage vs. Time . . . . . . . . . . . . . . . . . . 84 B.25 = 90 : vs. Number of Observations Processed . . . . . . . . . . . . . . . . 85 B.26 = 90 : Residuals vs. Number of Observations Processed . . . . . . . . . . . . 86 B.27 = 90 : Residuals Before Largest Outage vs. Time . . . . . . . . . . . . . . . 86 B.28 = 90 : Residuals After Largest Outage vs. Time . . . . . . . . . . . . . . . . 87 B.29 = 90 : NMSE Before Largest Outage vs. Time . . . . . . . . . . . . . . . . . 87 B.30 = 90 : NMSE After Largest Outage vs. Time . . . . . . . . . . . . . . . . . 88 B.31 ! = 45 : vs. Number of Observations Processed . . . . . . . . . . . . . . . . 89 B.32 ! = 45 : Residuals vs. Number of Observations Processed . . . . . . . . . . . . 90 B.33 ! = 45 : Residuals Before Largest Outage vs. Time . . . . . . . . . . . . . . . 90 B.34 ! = 45 : Residuals After Largest Outage vs. Time . . . . . . . . . . . . . . . . 91 B.35 ! = 45 : NMSE Before Largest Outage vs. Time . . . . . . . . . . . . . . . . . 91 B.36 ! = 45 : NMSE After Largest Outage vs. Time . . . . . . . . . . . . . . . . . 92 B.37 = 90 : vs. Number of Observations Processed . . . . . . . . . . . . . . . . 93 B.38 = 90 : Residuals vs. Number of Observations Processed . . . . . . . . . . . . 94 B.39 = 90 : Residuals Before Largest Outage vs. Time . . . . . . . . . . . . . . . 94 B.40 = 90 : Residuals After Largest Outage vs. Time . . . . . . . . . . . . . . . . 95 B.41 = 90 : NMSE Before Largest Outage vs. Time . . . . . . . . . . . . . . . . . 95 x B.42 = 90 : NMSE After Largest Outage vs. Time . . . . . . . . . . . . . . . . . . 96 B.43 = 270 : vs. Number of Observations Processed . . . . . . . . . . . . . . . . 97 B.44 = 270 : Residuals vs. Number of Observations Processed . . . . . . . . . . . 98 B.45 = 270 : Residuals Before Largest Outage vs. Time . . . . . . . . . . . . . . . 98 B.46 = 270 : Residuals After Largest Outage vs. Time . . . . . . . . . . . . . . . 99 B.47 = 270 : NMSE Before Largest Outage vs. Time . . . . . . . . . . . . . . . . 99 B.48 = 270 : NMSE After Largest Outage vs. Time . . . . . . . . . . . . . . . . . 100 xi List of Tables 2.1 J2 changes in the orbital elements . . . . . . . . . . . . . . . . . . . . . . . . . . 12 4.1 Baseline Initial Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 4.2 Variations from Baseline Case . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 4.3 Constants . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 4.4 Tracking Station Coordinates (Geodetic C.S.) . . . . . . . . . . . . . . . . . . . 32 4.5 Tracking Station Con guration Information . . . . . . . . . . . . . . . . . . . . 33 5.1 Baseline: Outage Details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 5.2 Baseline: Orbit Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 5.3 Baseline: Output Before/After Largest Outage . . . . . . . . . . . . . . . . . . 42 5.4 Baseline: All Observations (325 observations - 12915 s) . . . . . . . . . . . . . 43 5.5 i = 75:0 : Outage Details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 5.6 i = 75:0 : Orbit Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 5.7 i = 75:0 : Intermediate Output . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 5.8 i = 75:0 : Output for All Observations (326 - 13645 s) Processed . . . . . . . . 50 5.9 ! = 90 : Outage Details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 5.10 ! = 90 : Orbit Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 5.11 ! = 90 : Intermediate Output . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 5.12 ! = 90 : Output for All Observations (318 - 16460 s) Processed . . . . . . . . . 57 5.13 Additional Results Details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 B.1 i = 28:5 : Intermediate Output . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 xii B.2 i = 28:5 : Output for All Observations (307 - 7830 s) Processed . . . . . . . . . 73 B.3 rpra = 800km=1000km: Intermediate Output . . . . . . . . . . . . . . . . . . . . 76 B.4 rpra = 800km=1000km: Output for All Observations (387 - 20195 s) Processed . . 77 B.5 rpra = 500km=700km: Intermediate Output . . . . . . . . . . . . . . . . . . . . . 80 B.6 rpra = 500km=700km: Output for All Observations (352 - 6125 s) Processed . . . 81 B.7 = 0 : Intermediate Output . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 B.8 = 0 : Output for All Observations (311 - 29325 s) Processed . . . . . . . . . 85 B.9 = 90 : Intermediate Output . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 B.10 = 90 : Output for All Observations (323 - 12915 s) Processed . . . . . . . . . 89 B.11 ! = 45 : Intermediate Output . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 B.12 ! = 45 : Output for All Observations (301 - 13425 s) Processed . . . . . . . . . 93 B.13 = 90 : Intermediate Output . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 B.14 = 90 : Output for All Observations (301 - 16150 s) Processed . . . . . . . . . 97 B.15 = 270 : Intermediate Output . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 B.16 = 270 : Output for All Observations (326 - 41630 s) Processed . . . . . . . . 101 xiii Chapter 1 Introduction The term space debris refers to any human created object in space that is not part of an active mission. This broad scope includes propulsive stages, expired satellites, bolts, pieces of experimental hardware, fuel tanks that could contain unspent fuel, fragments from previous space debris collisions, and even a spare astronaut glove, which are a few examples of space debris cluttering low earth orbit (LEO) and geosynchronous orbit (GEO) paths [1]. Since the launch of Sputnik 1 in 1957 by the Soviet Union, many nations have engaged in multiple space launches, some of which have created space debris. While most modern satellites contain a GPS receiver to provide constant and accurate measurements, radar sites and optical telescopes provide the only tracking information avail- able for space debris. This tracking data is used to determine the orbit of each debris object as accurately as possible. Radar sites and optical telescopes provide limited coverage due to environmental and mechanical constraints that can restrict the amount of observations available for a piece of space debris depending on the position, velocity, orientation, and size of each piece of debris. With the use of this intermittent data, traditional orbit deter- mination methods that use the Extended Kalman Filter(EKF) can produce high residuals, which can lead to divergence and failure of the lter. High residuals for a piece of debris can increase the amount of fuel needed for nearby spacecraft to perform collision avoidance maneuvers. This can also lead to poor predictions for close pass scenarios, which may result in a collision. This research evaluates the conceptual application of an Iterated Extended Kalman Filter (IEKF) to the orbit determination of space debris in order to improve the accuracy of the orbit determination process. The results of the IEKF are compared to those of a 1 standard EKF. This chapter will discuss the history of space debris, previous work done concerning the space debris environment, and a more detailed problem description for this speci c study. 1.1 History The primary concern of a growing space debris population is the increased probability of collision with active missions; however, a secondary concern also exists, which involves debris colliding with other debris, thereby generating even more debris. In 1978 Kessler and Cour-Palais discussed the concern of an exponential growth of space debris due to fragments from previous collisions causing future collisions. This exponential growth phenomenon due to objects being added to the LEO environment faster than they can decay from orbit was later de ned as \The Kessler Syndrome." Kessler and Cour-Palais also discuss hypervelocity impacts, which cover the average impact velocity of 10 km=s of objects in LEO. Kessler points out that a collision at this velocity can produce \very high instantaneous pressures" on the order of 1000 GPa or 1:45 x 108 psi, which will usually result in the catastrophic destruction of both objects [2]. NASA also notes that in LEO a 1 cm diameter object with a few grams of mass contains the kinetic energy of a 250 kg mass moving 100 km=h [3]. Seeing a need to characterize the hazards of orbital debris to spacecraft and recommend mitigation techniques to help minimize the growth of the orbit debris environment, NASA established its orbital debris program in 1979. The orbital debris program grew quickly, and even included the construction of hypervelocity gun facilities in order to develop new shielding standards and technologies for the space shuttle and space station. In addition to providing research for shielding technologies, the orbital debris program also began to estab- lish mitigation standards. One of these important standards is the \25-year low-Earth orbit lifetime limitation." Under this standard, NASA stated that \maneuverable spacecraft that are terminating their operational phases at altitudes of less than 2,000 km above Earth shall be maneuvered to reduce their orbital lifetime, commensurate with 25-year low-Earth orbit 2 lifetime limitations [3]." This is very important in controlling the space debris environment as atmospheric drag will eventually decay the orbits of these expired satellites, removing them from the population [4]. The United States Strategic Command (USSTRATCOM) employs the Space Surveil- lance Network (SSN) to track and catalog objects approximately 10 cm in diameter and larger to aid in the orbit determination of those objects. Debris size, sensor coverage, orbit parameters, and spacecraft/debris attitude (cross-sectional area visible by sensor) each a ect the amount of tracking data available for every object. Kessler mentions the complexity of determining the space debris environment due to the di culty to track objects smaller than 10 cm in diameter [2]. Moreover, large gaps in data, can cause certain ltering methods to become unreliable and often biased, which can cause convergence to a trajectory other than the true trajectory or to divergence. NASA also suggested in 1999 that for every cataloged object of at least 10 cm diameter there are ten 1 cm and 10,000 1 mm diameter untrack- able objects. In 2011, there were approximately 22,000 cataloged objects, of which data is available publicly for about 16,000 of those objects [3]. In 2009, the collision between Iridium-33, an active communications satellite, and Cosmos-2251, a Russian satellite that ceased functioning in 1993, produced 1,275 pieces of debris that have been cataloged. Prediction software estimated a close approach between the two satellites; however, that prediction was only one of many conjunctions that were generated to warn of a pass within 5 km. Improvements in modeling, tracking, or lter- ing the data from measurements could reduce the uncertainty involved in forecasting close similar approaches [4]. In contrast to the accidental collision, the intentional destruction of the Chinese gov- ernment?s Fengyun-1C spacecraft, to test its anti-satellite weapon, produced 3,135 pieces of trackable debris. According to a NASA debris news update [5], the number of orbital debris 1 cm in diameter and larger associated with this satellite destruction is at least 150,000. The amount of debris generated by these two collisions has been said to have e ectively negated 3 the work completed to mitigate the growth of the space debris environment for the past twenty years [3]. 1.2 Previous Research Interactions with the space debris environment have been described as being active or passive. Active methods involve techniques that directly interact with the debris, which include the use of tethered satellites, projected lasers, robotics, and other similar techniques. While these methods are expensive, complicated, and largely untested, they could be used to alter the orbit of space debris to prevent collisions of debris with other debris as well as active missions. Passive methods involve techniques that do not physically interact with the debris. These methods include improvements to ltering methods, radar technology, and mitigation techniques similar to NASA?s 25-year low-Earth orbit lifetime requirement. Through these methods, the position and velocity of the debris can be estimated with less error, which can lead to fewer collisions with active missions. In 2008, Liou and Johnson[6] stated that the space debris environment is unstable, and mitigation techniques involving limiting the lifetime of expired satellites to 25 years will not be able to stabilize the population growth. In order to control the growth of the debris environment, they suggest the removal of large debris that has high collision probability. In their sensitivity study, the authors did not explore the technical or economic details of removing high-risk large space debris. If twenty debris objects are removed per year beginning in 2020, they predict the debris environment will stabilize by approximately 2080 [6]. Kessler revisited this phenomenon in 2010, and added that guidelines created in the last ten years must be adhered to while also deploying an active method to retrieve the most likely future debris sources in an attempt to prevent future spacecraft operators from facing limited spacecraft lifetimes due to collisions [7]. In 2012, Sibert and Borgeson[8] wrote a paper focused on the number of detailed con- junction assessments that must be performed and number of collision avoidance maneuvers 4 required to maintain a safe separation between objects. They began by showing that the orbital debris environment has already reached a stage where The Kessler Syndrome has taken over by presenting data that shows the debris environment growing even with zero fu- ture launches. In their paper, they also show that the number of conjunctions and collision avoidance maneuvers can be signi cantly reduced by reducing any or all of the following: state vector position error at time of update, state vector position error growth rate, and the delay between state vector updates. The most interesting part of their paper shows that im- proving the previously mentioned tracking performance signi cantly outweighs actual debris removal due to \very few debris objects turn[ing] out to be reliable repeat o enders[8]." 1.3 Problem Description The growing debris population, shown in Figure 1.1, provided by the NASA Debris Program O ce, coupled with the unpredictable nature of debris has created a need for orbit determination methods capable of obtaining more accurate estimates from fewer observa- tions. Figure 1.1: Debris Environment Trends 5 The purpose of this study is to determine if an Iterated Extended Kalman Filter can be used to yield better performance than the Extended Kalman lter in the orbit determination of space debris with poor apriori data and intermittent observations. A simulation was constructed using actual USN sensor locations to generate experimental observational data. This data was then used to compare the performance of the EKF and IEKF. Changes to the orbital parameters of the debris were analyzed to determine the e ects on lter performance. These initial conditions included variations to the following: altitude, inclination, longitude of the ascending node, argument of perigee, and true anomaly. It was necessary to determine the performance characteristics in order to de ne the comparison of the two lters. The di erence between the estimated state and actual state was chosen as the rst characteristic. Since the actual state is not available in a realistic application, the second characteristic to be compared was the residuals of the estimates. 6 Chapter 2 Relevant Theory The rst step in developing a simulation is to select the speci c coordinate systems and dynamical model. Orbit determination can use multiple reference frames, and many di erent dynamical models. First, the coordinate systems will be introduced. Following the coodinate systems, the equations of motion, orbital elements, and coordinate transformations will be de ned. A number of assumptions and limitations were taken into account and addressed during this research. Gravity was modeled as the gradient of the geopotential function using two-body and J2 e ects only. Atmospheric drag was not modeled, which had no e ect on the research as the measurements were generated using the same assumption. No other perturbation e ects were involved. 2.1 Reference Frames 2.1.1 Earth Centered Inertial Coordinate System (ECI) The origin of the ECI coordinate system, also referred to as the Geocentric equatorial coordinate system, is the center of the earth. The coordinate system is right-handed and has the following properties: the X axis points towards the vernal equinox with unit vector I, the Z axis points in the direction of the north pole with unit vector K, and the fundamental XY-plane coincides with the equatorial plane, de ning the Y axis with unit vector J. The coordinate system is not xed to the earth, and therefore doesn?t rotate with it, but rather the earth rotates relative to it. The ECI coordinate system is represented in Figure 2.1. 7 X Y Z Aligned with Vernal Equinox center of Earth I J K Figure 2.1: ECI coordinate system 2.1.2 Earth-Centered, Earth-Fixed Coordinate System (ECEF) The origin of the ECEF coordinate system, is the center of the spherical earth. The coordinate system is right-handed and has the following properties: the XE axis points towards the Greenwich meridian with unit vector i, the ZE axis points in the direction of the north pole with unit vector k, and the fundamental plane coincides with the equatorial plane, de ning the YE axis with unit vector j. The ECEF coordinate system rotates with the earth, and therefore is not an inertial coordinate system. Tracking station coordinates are typically provided in this system. See Figure 2.2 for a ECEF coordinate system compared with the ECI system. Angle is known as the sidereal angle. 8 X Y Z ? ? X E Y E Aligned with Greenwich Meridian Aligned with Vernal Equinox center of Earth Z E I i j J K,k Figure 2.2: ECEF coordinate system 2.1.3 Perifocal Coordinate System The fundamental plane of the perifocal coordinate system is the orbital plane. The coordinate system is right-handed and has the following properties: the P axis points towards perigee, the Q axis lies in the orbital plane rotated 90 degrees from P, and the W axis is perpendicular to the orbital plane and aligned with the angular momentum vector. The perifocal coordinate system is shown in Figure 2.3. V Q P V center of Earth r r p perigee object Orbital plane Figure 2.3: Perifocal coordinate system 9 2.2 Two-body Equations of Motion Newton?s second law of motion states that the sum of the external forces acting on a particle with mass, m, position, r, and velocity, v, is equal to the time rate of change of the product of the mass and velocity of the particle, which is represented in Equation 2.1. X F = ddt (mv) (2.1) If the mass of the particle is constant with respect to time, Equation 2.1 becomes Equation 2.2 where r represents the acceleration. X F = m r (2.2) Netwon?s law of universal gravitation states that any two particles exert a gravitational force on one another with magnitude, F = Gm1m2r2 (2.3) where m1 and m2 are the masses of the particles, G is the universal gravitational constant, and r is the magnitude of the position vector between them [9]. In the formulation of the two-body problem, the masses of those two bodies are assumed to be point masses, which satis es the assumption of particles used in the formulation of Newton?s laws. While Newton?s second law is valid in an inertial coordinate system, it can be modi ed to be valid in the geocentric equatorial (ECI) coordinate system shown in Figure 2.1. In this coordinate system, the acceleration of mass m2 relative to the Earth due to two-body e ects can be shown in Equation 2.4. r = r3r (2.4) The gravitational parameter, , is de ned such that 10 = G(m1 +m2) (2.5) The motion of m2 relative to m1 in the ECI coordinate system is shown in Figure 2.4. X Y Z O r m 1 m 2 Figure 2.4: Two-body problem in a relative coordinate system 2.3 J2 Perturbation It is well known that the Earth is a non-spherical, inhomogeneous body and not a point mass. To model this non-uniform mass distribution, the equations of motion canl be presented as the gradient of the gravity potential function U in Equation 2.6. r = ~rU (2.6) The dominating member of this gravity potential function is referred to as the oblateness of the Earth, J2. The oblateness of the Earth, which is a bulge of mass symmetric about the equator, is often modeled in orbit determination due to it being approximately 1000 times larger than any of the other perturbation forces in LEO. The J2 perturbation causes short period, long period, and secular changes in the orbital elements as presented in Table 2.1. 11 Table 2.1: J2 changes in the orbital elements Periodic Secular a (semi-major axis) (Longitude of the Ascending Node) e (eccentricity) ! (Argument of periapsis) i (inclination) (True anomaly) X E Y E Z E ? ? m 2 r Figure 2.5: Orbit in Spherical Coordinates In spherical coordinates, shown in Figure 2.5, the potential function of the J2 perturba- tion is given in Equation 2.7 [10], where J2 is the oblateness constant, re is the mean radius of the Earth, is the latitude, and is the longitude. UJ2 = r J2 re r 2 3 2 sin( ) 1 2 (2.7) The acceleration on m2 due to J2 is obtained by taking the gradient of UJ2. The equations of motion of m2 can then be expressed by Equation 2.8 as r = r3r + ~rUJ2 (2.8) Appendix A provides the transformation of ~rUJ2 into the ECI coordinate frame, for which the result is presented in Equation 2.9, where X1, X2, and X3 are the I, J, and K components of the ECI position vector . 12 ~rUJ 2 = rJ2 = 3 2 J2r 2 e r7 2 66 66 4 X1 (5X23 r2) X2 (5X23 r2) X3 (5X23 3r2) 3 77 77 5 (2.9) The updated equations of motion will consist of the two-body e ects combined with the J2 e ects. This is presented in Equation 2.10. a = r3 2 66 66 4 X1 X2 X3 3 77 77 5 + 3 2 J2r 2 e r7 2 66 66 4 X1 (5X23 r2) X2 (5X23 r2) X3 (5X23 3r2) 3 77 77 5 (2.10) 2.4 Orbital Plane Orientation in Space The use of orbital elements is a common way to describe an orbit in three-dimensional space. Instead of solely dealing with position and velocity vectors, one can use a set of the orbital elements that completely describe the orbit. For this problem, the classical Keplerian orbital elements set will be used. The classical orbital elements are as follows (also see Figure 2.6) [11]: 1. Semi-major axis (a) - speci es the size of the orbit 2. Eccentricity (e) - speci es the shape of the orbit 3. Inclination (i) - speci es the angle between the equatorial plane and orbital plane (one of two elements that specify the orientation of the orbit plane) 4. Longitude of the ascending node ( ) - speci es the angle measured counterclockwise in the XY plane of the ECI system from the X axis to the point where the spacecraft crosses the equatorial plane from south-to-north (one of two elements that specify the orientation of the orbit plane) 13 5. Argument of perigee (!) - speci es the angle measured in the orbit plane in the direction of motion, from the ascending node to perigee (speci es the orientation of the orbit in the orbit plane) 6. True anomaly ( ) - measured from periapsis to r in the direction of travel (speci es where on the orbit the spacecraft is located at any instant in time) ? X Y Z ? Line of Nodes V perigee r p i r V object DN AN Figure 2.6: Orbital Elements 2.5 Coordinate Transformations Orbital data can be represented by di erent values in di erent coordinate systems. Orbit determination deals with multiple types of data from various coordinate frames, and it is important to be able to clearly transform data from one system to another. This section discusses the necessary coordinate transformations used in this research. 14 2.5.1 Obtain ECI position and velocity vectors from the Keplerian orbital ele- ments Given the set of Keplerian orbital elements (a, e, i, , !, ), calculate r and v as follows [10]: Step 1 - Calculate the parameter, or semi-latus rectum, p p = a(1 e2) (2.11) Step 2 - Calculate r using the trajectory equation r = p1 +e cos( ) (2.12) Step 3 - Calculate the perifocal position vector rp = 2 66 66 4 r cos( ) r sin( ) 0 3 77 77 5 (2.13) Step 4 - Calculate the perifocal velocity vector vp = 2 66 66 4 q p sin( )q p (e+cos( )) 0 3 77 77 5 (2.14) 15 Step 5 - Calculate the transformation matrix from the perifocal to ECI [T] = 2 66 66 4 cos( )cos(!) sin( )cos(i)sin(!) cos( )sin(!) sin( )cos(i)cos(!) sin( )sin(i) sin( )cos(!) +cos( )cos(i)sin(!) sin( )sin(!) +cos( )cos(i)cos(!) cos( )sin(i) sin(i)sin(!) sin(i)cos(!) cos(i) 3 77 77 5 Step 6 - Calculate the ECI position and velocity vectors r = [T]rp (2.15) v = [T]vp (2.16) 2.5.2 Calculate the classical Keplerian orbital elements from ECI position and velocity vectors Step 1 - Calculate the orbital eccentricity, e r =j~rj v =j~vj e = v2 1 r r 1 (r v)v e =jej (2.17) 16 Step 2 - Calculate the semi-major axis, a h = r v h =jhj p = h 2 a = p1 e2 (2.18) Step 3 - Calculate the inclination, i hK = h K i = cos 1 h K h (2.19) Step 4 - Calculate the longitude of the ascending node, , and determine the correct quadrant n = K h n = 2 66 66 4 nI nj 0 3 77 77 5 I J K 17 n =jnj = cos 1 nI n (2.20) If nJ < 0, then > 180 ! = 360 (2.21) If nJ > 0, then < 180 (2.22) Step 5 - Calculate the argument of periapsis, !, and determine the correct quadrant ! = cos 1 n e ne (2.23) If eK > 0, then !> 180 !! = 360 ! (2.24) If eK < 0, then !< 180 (2.25) Step 6 - Calculate the true anomaly, , and determine the correct quadrant = cos 1 e r er (2.26) If r v < 0, then > 180 ! = 360 (2.27) 18 If r v > 0, then < 180 (2.28) 2.5.3 ECI to ECEF To convert from the ECI coordinate frame to the ECEF coordinate frame, rotate through the sidereal angle, , about the Z axis, as shown in Figure 2.2. This transformation is produced in Equation 2.29. rECEF = 2 66 66 4 cos( ) sin( ) 0 sin( ) cos( ) 0 0 0 1 3 77 77 5 rECI (2.29) 19 Chapter 3 State Estimation and Kalman Filtering The general formulation of the orbit determination problem includes the following well known relationships [12]: _X = F(X;t); X(tk) Xk (3.1) Yk = G(Xk;tk) + k; k = 1;:::;m (3.2) where Equation 3.1 represents an n-dimensional system of di erential equations, _X, as a nonlinear function, F, of the states, X, and time, t, and Equation 3.2 a measurement model, Yk, with each entry as a nonlinear function, G, of the epoch states, X(tk) Xk, the epoch time, tk, and an additive residual at the epoch time, i, for a number of measurements, m. The orbit determination is clearly a nonlinear estimation problem; however, Tapley [12] shows that if a reference trajectory which remains su ciently close to the true trajectory can be determined, the trajectory can be linearized by using a Taylor series expansion about the reference trajectory [12], X . This is an important step that will allow the application of linear theory. 3.1 Linearization The following linearization is provided by Tapley [12]. The state deviation, x(t), the dif- ference between the true trajectory, X and reference trajectory, X , is given in Equation 3.3. The observation residual, y(t), the di erence between the observation and expected obser- vation, is given in Equation 3.4. 20 x(t) = X(t) X (t) (3.3) y(t) = Y(t) Y (t) (3.4) The derivative of the state deviation, _x(t), will be useful in the linearization process, and is given in Equation 3.5. _x(t) = _X(t) _X (t) (3.5) The Taylor series expansion of Equation 3.1 and Equation 3.2 are given in Equation 3.6 and Equation 3.7 respectively where [ ] indicates that the partial derivative is evaluated at the reference solution, and OF and OG are orders of magnitude of the higher order terms in the expansion that will be ignored due to their small magnitudes. _X(t) = F(X;t) = F(X ;t) + @F(t) @X(t) [X(t) X (t)] +OF[X(t) X (t)]2 (3.6) Yk = G(Xk;tk) + k = G(X k;tk) + @G @X k [X(tk) X k(t)]k +OG[Xk(t) X k(t)]2 + k (3.7) Using the conditions _X = F(X ;t) and Y k = G(X k;tk), Equation 3.6 and Equation 3.7 can be written as Equation 3.10 and Equation 3.11 where the de nitions of A(t) and Hk are given below in Equation 3.8 and Equation 3.9 respectively. A(t) = @F(t) @X(t) (3.8) 21 Hk = @G @X k (3.9) _x(t) = A(t)x(t) (3.10) yk = Hkxk + k; k = 1;:::;m (3.11) 3.2 Kalman Filter The sequential estimation algorithm, often called the Kalman Filter, processes measure- ments individually as they are received. From Tapley [12], the Kalman lter equations are given in Equation 3.12, Equation 3.13, and Equation 3.14. These equations assume that the propagated state vector, ^x k , propagated covariance matrix, P k , current observation, Yk, and current observation covariance, Rk are known. Kk = P kHTk HkP kHTk + Rk 1 (3.12) ^x+k = ^x k + Kk Yk Hk^x k (3.13) P+k = [I KkHk]P k (3.14) 3.2.1 Extended Kalman Filter Large deviations from the reference trajectory can cause lter divergence. The Kalman Filter performs an update to the reference trajectory only after all observations have been 22 processed, regardless of processing each measurement individually. In constrast, the Ex- tended Kalman Filter performs an update to the reference trajectory after each measurement has been processed. The following algorithm is provided by Tapley [12]. 1. Integrate from tk 1 to tk using the dynamical model ^x+k 1!^x k (3.15) P+k 1!P k (3.16) 2. Compute the residual and Hk matrix yk = Yk G(X k;tk) (3.17) Hk = @G@x ^x+k 1 (3.18) 3. Compute the Kalman gain Kk = PkHTk HkPkHTk + Rk 1 (3.19) 4. Compute the updated state and covariance ^x+k = ^x k + Kk Yk Hk^x k (3.20) P+k = [I KkHk]P k (3.21) 23 5. Update the reference trajectory to the current state estimate. If additional measure- ments exist go to 1, otherwise stop. ^X = ^x+k (3.22) 3.2.2 Iterated Extended Kalman Filter In order to help minimize linearization errors in highly nonlinear problems, higher order methods are often applied to a problem. In the linearization process used in Equation 3.6 and Equation 3.7, the measurement equation was expanded about the reference solution. With the use of the extended Kalman lter, the reference solution is updated by propagation, ^x+k 1 ! ^x k , and processing of a measurement, ^x k ! ^x+k , which is the step including the linearization about ^x k . Since ^x+k represents the best estimate of the state xk after the measurement is processed, the measurement equation could be re-linearized about this new best estimate of the state to obtain a more accurate solution through as many iterations as are desired. The application of higher order lters carries the weight of additional computational cost. In 1974, Gelb includes the iterated extended Kalman lter in his book with empha- sis on \computation time required to mechanize the lter [13]." Jazwinski shows the IEKF converges more quickly to the exact solution than the EKF; however,he notes that the EKF will converge to the exact solution given enough measurements [14]. The solutions of both lters are subject to changes in the a priori covariance matrix, measurement covariance, and availability of measurements. As the a priori covariance matrix grows, measurement covari- ance decreases, and number of measurements decreases, the performance IEKF is superior to the EKF [14]. This derivation of the iterated extended Kalman lter algorithm is summarized from Simon [15]. IEKF Algorithm 24 1. Initialize the IEKF to the EKF estimate ^x+k;0 = ^x+k (3.23) P+k;0 = P+k (3.24) 2. Expanding the measurement equation around ^x+k;i 1 yields Equation 3.26. Hk;i = @G@x ^x+k;i 1 (3.25) yk;i = Yk G(x+k;i 1;tk) (3.26) 3. Calculate the Kalman gain, Kk;i for the current iteration. Kk;i = P kHTk;i[Hk;iP kHTk;i + Rk] 1 (3.27) 4. Calculate the covariance matrix, P, and the new best estimate of the state, x+k;i+1, for the current iteration. P+k;i = (I Kk;iHk;i)P k (3.28) ^x+k;i = ^x+k + Kk;i[yk;i Hk;i(^x+k ^x+k;i 1)] (3.29) 5. Check for convergence. If convergence criterion not met, initialize the next iteration and goto 2, otherwise stop. 25 ^x+k;i+1 = ^x+k;i (3.30) P+k;i+1 = P+k;i (3.31) 26 Chapter 4 Simulation A simulation was constructed to compare the performance of the EKF and IEKF in the orbit determination scenario involving poor apriori information and intermittent observation data for a piece of space debris. The simulation consisted of multiple functions and programs to construct the orbit determination problem, and will be discussed in this chapter. The initial conditions and statistical data, which vary for most test cases, will also be presented. 4.1 Baseline Initial Conditions The initial conditions chosen for the baseline case included an eccentricity and inclina- tion similar to the International Space Station, and are presented in Table 4.1. This orbit was an important candidate for testing due to the future tra c and amount of debris existing at that altitude. Table 4.1: Baseline Initial Conditions Orbital Elements Value Perigee Altitude (m) 375134.9875 Apogee Altitude (m) 400249.8545 Inclination Angle (degrees) 51.6164 Longitude of Ascending Node (degrees) 270.1324 Argument of Periapsis (degrees) 0.0135 True Anomaly (degrees) 0.0297 The true initial state is given in Equation 4.1. The initial standard deviations, shown in Equation 4.2, were chosen based on the results by Kelso [16] for approximately 5 days 27 of propagation error of a TLE dataset. From these data, the initial covariance matrix in Equation 4.3 was chosen, where r is the covariance of the position and v is the covariance of the velocity. 2 66 66 66 66 66 66 66 4 X Y Z Vx Vy Vz 3 77 77 77 77 77 77 77 5 = 2 66 66 66 66 66 66 66 4 165109:2914 m 6748615:7425 m 188744:0459 m 4771:0926 m=s 284:9491 m=s 6023:9774 m=s 3 77 77 77 77 77 77 77 5 (4.1) r = 2500m v = 10m=s (4.2) P0 = 2 66 66 66 66 66 66 66 4 2rx 0 0 0 0 0 0 2ry 0 0 0 0 0 0 2rz 0 0 0 0 0 0 2vx 0 0 0 0 0 0 2vy 0 0 0 0 0 0 2vz 3 77 77 77 77 77 77 77 5 = 2 66 66 66 66 66 66 66 4 25002 0 0 0 0 0 0 25002 0 0 0 0 0 0 25002 0 0 0 0 0 0 102 0 0 0 0 0 0 102 0 0 0 0 0 0 102 3 77 77 77 77 77 77 77 5 (4.3) 4.2 Test Cases In order to properly compare the performance of the EKF and IEKF, multiple sets of initial conditions and observational data were evaluated. The altitude, inclination, longitude 28 of the ascending node, argument of periapsis, and true anomaly were altered individually for each case. A summary of non-baseline cases can be found in Table 4.2, which are listed by how each case di ers from the baseline case. Table 4.2: Variations from Baseline Case Case Number Variation of Value of Varied Terms 2 Altitude ra = 500km and rp = 700km 3 Altitude ra = 800km and rp = 1000km 4 Inclination i = 28:5 5 Inclination i = 75 6 Longitude of the Ascending Node = 45 7 Longitude of the Ascending Node = 90 8 Argument of Periapsis ! = 45 9 Argument of Periapsis ! = 90 10 True Anomaly = 90 11 True Anomaly = 270 4.3 Constants Many constants are used when performing orbit determination. For completeness, Ta- ble 4.3 lists the values of each constant used in this simulation. Table 4.3: Constants = 3:9860044150 1014 m3s2 J2 = 0:1082635666 10 2 !e = 7:2921156001 10 5 rads re = 6378:13630 km 29 4.4 Dynamical Model The dynamical model represents the equations of motion of the object, and is used to propagate the state, Equation 4.4, forward through time. The dynamical model in this simulation used Equation 2.10 to form the derivatives of the states, which are given in Equation 4.5. This equation includes accelerations due to two-body motion and the J2 oblateness perturbation e ects. X = 2 66 66 66 66 66 66 66 4 rx ry rz vx vy vz 3 77 77 77 77 77 77 77 5 = 2 66 66 66 66 66 66 66 4 X1 X2 X3 X4 X5 X6 3 77 77 77 77 77 77 77 5 (4.4) _X = 2 66 66 66 66 66 66 66 4 vx vy vz ax ay az 3 77 77 77 77 77 77 77 5 = 2 66 66 66 66 66 66 66 4 X4 X5 X6 X1r4+32 J2r2eX1(5X23 r2) r7 X2r4+32 J2r2eX2(5X23 r2) r7 X3r4+32 J2r2eX3(5X23 3r2) r7 3 77 77 77 77 77 77 77 5 (4.5) 4.5 Measurement Model The measurement consists of the equations for each measurement type. For this simu- lation, range-only measurements, YR, were used, and are represented in Equation 4.6, where XS, YS, and ZS are the ECEF coordinates of the tracking station and XE, YE, and ZE are the ECEF coordinates of the object. Equation 4.7 is used to transform the ECEF coordi- nates into the ECI coordinates of the model, which are also the rst three entries in the ECI 30 state vector. Substituting this result into Equation 4.6 yields Equation 4.8, the measurement model in ECI coordinates. YR = (XE XS)2 + (YE YS)2 + (ZE ZS)2 1 2 (4.6) 2 66 66 4 XE YE ZE 3 77 77 5 = 2 66 66 4 cos( ) sin( ) 0 sin( ) cos( ) 0 0 0 1 3 77 77 5 2 66 66 4 X1 X2 X3 3 77 77 5 (4.7) YR = (X1cos( ) +X2sin( ) XS)2 + ( X1sin( ) +X2cos( ) YS)2 + (X3 ZS)2 1 2 (4.8) 4.6 Tracking Station Network The tracking stations used to generate experimental data in the simulation were taken from Vallado [17]. The sites used represent actual tracking stations that are a part of the US Space Surveillance Network, and are represented in Table 4.4. This network consists of mainly northern and western hemisphere stations, which lead to large gaps in the measure- ment data. To simplify the tracking network, the data in Table 4.5 was implemented. A maximum range of 5,000 km was chosen to ensure range would not prevent the data generator from producing measurements as the maximum range generated is on the order of 1,000 km. To mitigate errors associated with sensing near the horizon, the minimum elevation angle required to take a measurement was established at 20 degrees, which Vallado suggests for optical sensors [17]. The sensors were assumed to have no mechanical, geographical, or political constraints that could limit the azimuth. The range noise and bias were assumed to be 5 m and 0 m respectively. 31 Table 4.4: Tracking Station Coordinates (Geodetic C.S.) Site Latitude(degrees) Longitude(degrees) Altitude(meters) Flyingdales, UK 54.37 -0.67 338.9 Ascension, Atlantic -7.91 -14.40 56.1 Clear, AK 64.29 -149.19 213.3 Antigua, West Indies 17.14 -61.79 0.5 Cape Cod, MA 41.75 -70.54 80.3 Beale, CA 39.14 -121.35 115.7 Shemya, AK 52.74 -174.09 89.8 Thule, Greenland 76.57 -68.30 424.7 Cavalier, ND 48.72 -97.90 347.3 Kaena Point, HI 21.57 -158.27 300.2 Kwajalein Atoll, Paci c Ocean 9.39 167.48 62.7 Diego-Garcia, Indian Ocean -7.41 72.45 -61.2 Moron (MOSS), Spain 37.10 -5.37 287.0 Eglin, FL 30.57 -86.21 34.7 Socorro, NM 33.82 -106.66 1510.2 Globus II, Norway 70.37 31.13 1.2 32 Table 4.5: Tracking Station Con guration Information Maximum Range 5000 km Minimum Elevation 20 Maximum Elevation 90 Minimum Azimuth 0 Maximum Azimuth 360 Range Noise 5 m Range Bias 0 m The simulated measurements were generated with an assumed maximum range and minimum elevation angle from the tracking station to maintain a sense of realism, as well as to add measurement outages. The measurement rate was limited to the integration rate of 15 seconds. The measurements were assumed to originate from tracking stations with zero bias, and 5 m noise in range. To simulate a debris object for which tracking data is available, the initial estimate of the state of the debris being tracked was obtained from a simulated two-line element (TLE) data set on the order of two to ve days old. This method generated a large apriori covariance matrix with which to begin the estimation process. 4.7 Propagator The propagation method chosen for this study was a fourth order Runge-Kutta. The Runge-Kutta o ers increased accuracy over the basic Euler method by calculating additional interior points along the path from beginning of the step, t0, to the end of the step, t0 + t. The fourth order Runge-Kutta is widely documented, and can be found below in Equation 4.9 through Equation 4.13 [10]. 33 k1 = f (X(t0);t0) (4.9) k2 = f X(t0) + t2 k1;t0 + t2 (4.10) k3 = f X(t0) + t2 k2;t0 + t2 (4.11) k4 = f X(t0) + tk3;t0 + t2 (4.12) X(to + t) = t6 (k1 + 2k2 + 2k3 +k4) (4.13) 4.8 Procedure For each case, a set of initial conditions was de ned. Those initial conditions were used with the dynamics model, measurement model, ground station data, propagator, and data generator to generate a set of data for that speci c case. The generated data was then used with the dynamics model, measurement mode, propagator, EKF algorithm, IEKF algorithm, and a post processing script to generate MATLAB plots and data of the performance of each lter. This data will be presented in the next chapter. 34 Chapter 5 Results 5.1 Measure of Solution Accuracy The vectors and represent the di erence in the best estimate of the state given by EKF, ^xEKF, and IEKF, ^xIEKF, respectively and the true solution, xTrue. =j^xEKF xTruej =j^xIEKF xTruej The normalised mean squared error (NMSE), given by Equation 5.1, will be used to de ne a performance parameter, , which is given in Equation 5.2. By de nition, if > 0 the IEKF solution is more accurate; however, if < 0 the EKF solution is more accurate. n = " nX i=1 (^xi xi)2 ^xi #1 2 (5.1) = n;EKF n;IEKF 1 (5.2) 5.2 Baseline Case The baseline case consisted of 325 observations spanning 8 batches and 12,915 seconds, of which details can be found in Table 5.1. The largest outage occurs between batch 6 and batch 7, and lasts for 4490 seconds, which is approximately 80% of an orbit. The orbital 35 parameters for this case can be found in Table 5.2, and represent an orbit similar to that of the ISS. Table 5.1: Baseline: Outage Details Batch Number 1 2 3 4 5 6 7 8 First Measurement No. 1 46 92 115 161 211 259 292 Last Measurement No. 45 91 114 160 210 258 291 325 Start Time (s) 320 1195 1580 5590 6960 7290 12015 12750 End Time (s) 540 1420 1690 5815 7205 7525 12175 12915 Outage Before Batch (s) 320 655 160 3900 1145 85 4490 575 Table 5.2: Baseline: Orbit Parameters Orbital Elements Value Perigee Altitude (m) 375134.9875 Apogee Altitude (m) 400249.8545 Inclination Angle (degrees) 51.6164 Longitude of Ascending Node (degrees) 270.1324 Argument of Periapsis (degrees) 0.0135 True Anomaly (degrees) 0.0297 The residuals, the di erence between the observation and the expected observation, for the full set of data are shown in Figure 5.1. The residuals for both the EKF and IEKF remain small for the entire simulation; however, the residuals for the IEKF remain closer to zero than those of the EKF. This can be attributed to the re-linearization of the observation model that is carried out multiple times for each observation in the IEKF. This trend will remain present in the analysis of each case. To allow closer inspection, Figure 5.2 and Figure 5.3 show the residuals for batch 6 and batch 7 respectively. Batch 6 consisted of 47 measurements spanning 235 seconds, while batch 7 consisted of 32 measurements spanning 160 seconds. Due to the very small residual prior to the outage, the residual does not grow signi cantly in either lter between the nal 36 0 100 200 300 400 ?20 ?10 0 10 20 Number of Observations Processed Residual ? meters EKF IEKF Figure 5.1: Baseline: Residuals vs. Number of Observations Processed measurement in batch 6 to the rst measurement in batch 7. The IEKF residuals for batch 7, Figure 5.3, remain clearly smaller than those for the EKF. 37 7250730073507400745075007550 ?20 ?10 0 10 20 Time Elapsed ? seconds Residual ? meters EKF IEKF Figure 5.2: Baseline: Residuals for Batch #6 1.2 1.205 1.21 1.215 1.22x 104 ?20 ?10 0 10 20 Time Elapsed ? seconds Residual ? meters EKF IEKF Figure 5.3: Baseline: Residuals for Batch #7 38 The performance parameter, , shown in Figure 5.4, remains above the reference zero line for the entire simulation. This serves as a quick reference that the IEKF produced more accurate estimates of the true state than the EKF. The NMSE values for batch 6 and 7 are shown together in Figure 5.5. The NMSE values for batch 6 and batch 7 individually are shown in Figure 5.6 and Figure 5.7 respectively. It should be noted that the y-axis scale of these two gures are equal, which allows a direct comparison. The NMSE of the EKF grows approximately 60% during the outage, while the NMSE of the IEKF grows approximately 30%. This growth from an outage is expected; however, the IEKF recovers to a lower value faster than the EKF. 0 100 200 300 400?50 0 50 100 150 200 250 300 Number of Observations Processed ? ? Ref line Figure 5.4: Baseline: vs. Number of Observations Processed 39 0.7 0.8 0.9 1 1.1 1.2 1.3x 1040 0.2 0.4 0.6 0.8 1x 10?5 Time Elapsed ? seconds ? n EKF IEKF Figure 5.5: Baseline: NMSE for Batch #6 and Batch #7 72507300735074007450750075500 0.2 0.4 0.6 0.8 1x 10?5 Time Elapsed ? seconds ? n EKF IEKF Figure 5.6: Baseline: NMSE for Batch #6 40 1.2 1.205 1.21 1.215 1.22x 1040 0.2 0.4 0.6 0.8 1x 10?5 Time Elapsed ? seconds ? n EKF IEKF Figure 5.7: Baseline: NMSE for Batch #7 41 The performance of the lters was evaluated speci cally for batch 6 and batch 7. Ta- ble 5.3 contains data concerning the state, residual, and NMSE for both lters. For batch 6, the and performance values, which represent the di erence in the estimated state and true state, remained small for both lters prior to the outage. The residuals of both lters were also small; however, the residual of the IEKF was approximately an order of magnitude smaller than that of the EKF. Following the outage, and processing of batch 7, both residuals have increased. The NMSE of the IEKF also remained smaller than that of the EKF. This is due again to the IEKF processing an observation more than once, which leads to an arti cally larger set of observations for the IEKF. Table 5.3: Baseline: Output Before/After Largest Outage End of Batch #6 End of Batch #7 State (EKF) (IEKF) (EKF) (IEKF) X(m) 0.7734 0.7548 5.0188 0.9288 Y(m) 0.5084 0.1123 6.4405 0.5038 Z(m) 2.7568 1.1256 3.8424 0.1541 Vx(m/s) 0.0115 0.0005 0.0065 0.0007 Vy(m/s) 0.0006 0.0001 0.0034 0.0003 Vz(m/s) 0.0061 0.0023 0.0043 0.0006 yk(m) 0.7927 0.0794 7.4162 1.0580 NMSE 4.0298E-006 7.1876E-007 7.0072E-006 8.2265E-007 After batch 8 has been processed, all data has been consumed. The performance values can be found in Table 5.4. Given additional observations, the accuracy of the EKF improved, which can be seen by the reduction in size of the residual and NMSE values when compared to Table 5.3. 42 Table 5.4: Baseline: All Observations (325 observations - 12915 s) State (EKF) (IEKF) X(m) 2.1986 0.1477 Y(m) 4.9208 0.4030 Z(m) 5.7930 0.5380 Vx(m/s) 0.0112 0.0007 Vy(m/s) 0.0028 0.0005 Vz(m/s) 0.0018 0.0012 yk(m) 3.2716 1.0989 NMSE 4.8866E-006 5.3018E-007 While the and values for the EKF continue to remain small, the values for the position components are at least an order of magnitude larger than those of the IEKF. Combining this result with the smaller residual and smaller NMSE value, the IEKF produced a more accurate estimate of the state for the baseline case. 5.3 Inclination Variation Case !i = 75:0 The i = 75:0 case consisted of 326 observations spanning 9 batches and 13,410 seconds, of which details can be found in Table 5.5. The largest outage occurs between batch 7 and batch 8, and lasts for 5895 seconds, which is approximately 105% of an orbit. This outage and the initial outage, the time before the beginning of batch 1, are approximately 30% larger than those speci c outages for the baseline case. These outage increases are due to the nite number of ground stations modeled. The orbital parameters for this case can be found in Table 5.6. This case di ers from the baseline case only by the inclination angle. Table 5.5: i = 75:0 : Outage Details Batch Number 1 2 3 4 5 6 7 8 9 First Measurement No. 1 22 71 117 153 183 202 232 279 Last Measurement No. 21 70 116 152 182 201 231 278 326 Start Time (s) 730 905 1330 2755 6330 6630 6970 13010 13410 End Time (s) 830 1145 1555 2930 6475 6720 7115 13240 13645 Outage Before Batch (s) 730 75 185 1200 3400 155 250 5895 170 43 Table 5.6: i = 75:0 : Orbit Parameters Orbital Elements Value Perigee Altitude (m) 375134.9875 Apogee Altitude (m) 400249.8545 Inclination Angle (degrees) 75.121 Longitude of Ascending Node (degrees) 270.1324 Argument of Periapsis (degrees) 0.0135 True Anomaly (degrees) 0.0297 The residuals for the full set of data are shown in Figure 5.8. The IEKF residuals remain small, compared to the EKF residuals, for all observations; however, the EKF residuals grow to near 200 m for some observations. The 30% increase in the initial outage coupled with the size of the rst batch, 20 observations compared to 45 observations for the baseline case, leads the EKF to a less accurate estimate early in the simulation. For each batch of data, the residual of the EKF begins to shrink; however, this process is interrupted by multiple outages. 44 0 100 200 300 400 ?200 ?100 0 100 200 Number of Observations Processed Residual ? meters EKF IEKF Figure 5.8: i = 75 : Residuals vs. Number of Observations Processed Figure 5.9 and Figure 5.10 show the residuals for batch 7 and batch 8 respectively. Batch 7 consisted of 29 measurements spanning 145 seconds, while batch 8 consisted of 46 measurements spanning 230 seconds. For batch 7, not enough observations existed to allow the EKF to attain an estimate as accurate as the IEKF. The residuals for the IEKF remained small for the entire batch. For batch 8, the EKF processed enough observations to appear to converge to an accurate estimate very close to that of the IEKF, which remained close to zero for the entire batch. 45 6950 7000 7050 7100 7150 ?200 ?100 0 100 200 Time Elapsed ? seconds Residual ? meters EKF IEKF Figure 5.9: i = 75 : Residuals for Batch #7 1.3 1.3051.31 1.3151.32 1.325x 104 ?200 ?100 0 100 200 Time Elapsed ? seconds Residual ? meters EKF IEKF Figure 5.10: i = 75 : Residuals for Batch #8 46 The performance parameter, , shown in Figure 5.11, remains above the reference zero line for the entire simulation. This serves as a quick reference that the IEKF produced more accurate estimates of the true state than the EKF, which is the same result as the baseline case. The NMSE values for batch 7 and 8 are shown together in Figure 5.12. The NMSE of the EKF begins and ends greater than the NMSE of the IEKF for both batches. The NMSE values for batch 7 and batch 8 individually are shown in Figure 5.13 and Figure 5.14 respectively. For batch 7, the NMSE steadily decreases as observations are processed. Similar to the residual, this trend would have continued given additional observations. For batch 8, the NMSE of the EKF remains much higher, although relatively constant, than the NMSE of the IEKF. This shows that the IEKF estimate of the state remains much closer to the true state than the EKF estimate. 0 100 200 300 400?20 0 20 40 60 80 100 Time (s) ? ? Ref line Figure 5.11: i = 75 : vs. Number of Observations Processed 47 0.6 0.8 1 1.2 1.4x 1040 0.5 1 1.5 2 2.5x 10?3 Time Elapsed ? seconds ? n EKF IEKF Figure 5.12: i = 75 : NMSE for Batch #7 and Batch #8 6950 7000 7050 7100 71500 0.5 1 1.5 2 2.5 3x 10?3 Time Elapsed ? seconds ? n EKF IEKF Figure 5.13: i = 75 : NMSE for Batch #7 48 1.3 1.305 1.31 1.315 1.32 1.325x 1040 0.5 1 1.5 2 2.5 3x 10?3 Time Elapsed ? seconds ? n EKF IEKF Figure 5.14: i = 75 : NMSE for Batch #8 49 The performance of the lters was evaluated speci cally for batch 7 and batch 8. Table 5.7 contains data concerning the state, residual, and NMSE for both lters. For batch 7 and performance values remained approximately two orders of magnitude smaller for the IEKF than those of the EKF. This can be attributed to the high residual for the EKF at the end of batch 7. For batch 8 the residual for the EKF approaches a small number as seen in Figure 5.10. However, this low residual does not correspond with a low NMSE value. This is a representation of apparent divergence, which is convergence to the wrong value. Table 5.7: i = 75:0 : Intermediate Output End of Batch #7 End of Batch #8 State (EKF) (IEKF) (EKF) (IEKF) X(m) 32.0360 0.3094 317.7501 7.2187 Y(m) 175.8334 2.6102 109.3774 2.8083 Z(m) 207.6050 1.5862 48.3734 1.2793 Vx(m/s) 0.4586 0.0116 0.3879 0.0098 Vy(m/s) 0.2080 0.0023 0.3018 0.0040 Vz(m/s) 0.0540 0.0021 0.3622 0.0063 yk(m) 110.7678 0.8292 1.7156 1.1389 NMSE 9.0870E-004 2.2975E-005 3.9836E-004 9.3518E-006 After batch 9 has been processed, all data has been consumed. The performance values can be found in Table 5.8. Similar to the results shown in Table 5.7, the IEKF produced a more accurate estimate of the true state, as shown by the values being at least one order of magnitude smaller than the values. The NMSE of the IEKF was approximately two orders of magnitude smaller than that of the EKF after processing all observations. Table 5.8: i = 75:0 : Output for All Observations (326 - 13645 s) Processed State (EKF) (IEKF) X(m) 399.6702 9.7243 Y(m) 234.1288 4.8175 Z(m) 92.1286 2.1426 Vx(m/s) 0.1945 0.0055 Vy(m/s) 0.4273 0.0073 Vz(m/s) 0.3399 0.0064 yk(m) 25.9034 0.8774 NMSE 1.2476E-003 2.9886E-005 50 5.4 Argument of Perigee Case !! = 90 The ! = 90 case consisted of 318 observations spanning 8 batches and 16,460 seconds, of which details can be found in Table 5.9. While the largest outage is smaller than the largest outage of the previous case (i = 75:0 ), the initial outage is 4625 seconds. The poor apriori covariance matrix is propagated for 85% of an orbit before the rst observation is processed. The orbital parameters for this case can be found in Table 5.10. This case di ers from the baseline case only by the argument of periapsis. Table 5.9: ! = 90 : Outage Details Batch Number 1 2 3 4 5 6 7 8 First Measurement No. 1 49 97 133 161 208 222 270 Last Measurement No. 48 96 132 160 207 221 269 318 Start Time (s) 4625 4895 5660 7060 10235 10770 15915 16220 End Time (s) 4860 5130 5835 7195 10465 10835 16150 16460 Outage Before Batch (s) 4625 35 530 1225 3040 305 5080 70 Table 5.10: ! = 90 : Orbit Parameters Orbital Elements Value Perigee Altitude (m) 375134.9875 Apogee Altitude (m) 400249.8545 Inclination Angle (degrees) 51.6164 Longitude of Ascending Node (degrees) 270.1324 Argument of Periapsis (degrees) 90.1943 True Anomaly (degrees) 0.0297 The residuals for the full set of data are shown in Figure 5.15. The IEKF residual re- mains small for all observations; however, the EKF residual becomes large very quickly. The combination of very accurate measurements (5 m measurement noise) and the large covari- ance matrix following the starting outage cause the EKF to diverge quickly. The accurate measurements cause the covariance matrix to approach zero before enough measurements are processed to allow for the lter to converge. The small covariance matrix forces the Kalman gain calculation to also be small. If the Kalman gain is small, the state update 51 will also be small, regardless of the size of the residual. Figure 5.16 and Figure 5.17 provide a closer look at the residuals for batch 6 and batch 7 respectively. Figure 5.16 shows the residual of the EKF growing as more measurements are processed, which shows convergence to an incorrect solution. In Figure 5.17 the EKF residuals become smaller with additional observations; however, the scale of the graph shows that the EKF residuals remain large. 0 100 200 300 400?4 ?3 ?2 ?1 0 1x 107 Number of Observations Processed Residual ? meters EKF IEKF Figure 5.15: ! = 90 : Residuals vs. Number of Observations Processed 52 1.076 1.078 1.08 1.082 1.084x 104?2 ?1 0 1 2x 105 Time Elapsed ? seconds Residual ? meters EKF IEKF Figure 5.16: ! = 90 : Residuals for Batch #6 1.59 1.595 1.6 1.605 1.61 1.615x 104?2 ?1 0 1 2x 107 Time Elapsed ? seconds Residual ? meters EKF IEKF Figure 5.17: ! = 90 : Residuals for Batch #7 53 The performance parameter, , shown in Figure 5.18, remains above the reference zero line for the entire simulation. Due to the scale on the graph, it is di cult to see this for approximately the rst 150 observations. The NMSE values for batch 6 and 7 are shown together in Figure 5.19. The NMSE of the EKF begins and ends greater than the NMSE of the IEKF for both batches. The NMSE values for batch 6 and batch 7 individually are shown in Figure 5.20 and Figure 5.21 respectively. For both batches, the IEKF NMSE value remains near zero, while the EKF NMSE is many orders of magnitude larger. Figure 5.21 shows the EKF NMSE becoming smaller similar to the residual from Figure 5.17. These large EKF NMSE values are a result of the innacuracy of the EKF estimate. 0 100 200 300 400?5 0 5 10 15 20x 105 Number of Observations Processed ? ? Ref line Figure 5.18: ! = 90 : vs. Number of Observations Processed 54 1 1.2 1.4 1.6 1.8x 1040 2 4 6 8 10 Time Elapsed ? seconds ? n EKF IEKF Figure 5.19: ! = 90 : NMSE for Batch #6 and Batch #7 1.076 1.078 1.08 1.082 1.084x 1040 0.05 0.1 0.15 0.2 0.25 0.3 0.35 Time Elapsed ? seconds ? n EKF IEKF Figure 5.20: ! = 90 : NMSE for Batch #6 55 1.59 1.595 1.6 1.605 1.61 1.615x 1040 2 4 6 8 10 Time Elapsed ? seconds ? n EKF IEKF Figure 5.21: ! = 90 : NMSE for Batch #67 56 Table 5.11: ! = 90 : Intermediate Output End of Batch #6 End of Batch #7) - After State (EKF) (IEKF) (EKF) (IEKF) X(m) 407280.0317 4.9537 2340837.9835 8.0329 Y(m) 267709.9240 9.2157 566814.6803 0.7537 Z(m) 184296.0433 5.4234 3025870.0294 0.8134 Vx(m/s) 147.8404 0.0028 4731.1437 0.0081 Vy(m/s) 793.0966 0.0001 11469.8375 0.0044 Vz(m/s) 361.9798 0.0093 3859.5161 0.0169 yk(m) 183061.6412 7.8710 2134791.3973 3.3742 NMSE 3.3191E-001 8.2508E-006 3.0272E+000 6.9350E-006 Similar to the previous cases, the lter performance was compared for the batch before and the batch following the largest outage. This information can be found in Table 5.11. The EKF estimate has diverged from the true solution as clearly noted by the values and residual; however, the performance of the IEKF did not su er due to the initial outage or any outages throughout the simulation. This was due to the formulation of the IEKF algorithm. The covariance matrix is not updated until the iterated solution for each observation has converged, which is found in Equation 3.28. Table 5.12: ! = 90 : Output for All Observations (318 - 16460 s) Processed State (EKF) (IEKF) X(m) 3197640.7249 5.3589 Y(m) 2571119.2984 1.2943 Z(m) 3567720.2319 5.9253 Vx(m/s) 7062.2632 0.0104 Vy(m/s) 8570.047 0.0042 Vz(m/s) 4803.1641 0.0147 yk(m) 3662603.4439 0.3185 NMSE 9.3240E+000 1.8473E-005 Following the processing of all additional measurements, the nal performance values can be found in Table 5.12. The IEKF continues to produce an accurate estimate of the true state, from which it di ers by a maximum of approximately 6 m in any position component. The EKF estimate has diverged completely, which can be stated from the di erence from the true state as well as the residual value. 57 5.5 Additional Results Results for additional cases can be found Appendix B. The results of each case belong to one of three categories, which are de ned from the previously discussed cases. Cases behaving similar to the baseline case include results where both lters provided accurate and similar results. Apparent divergence in the EKF, similar to the i = 75:0 case, de nes the second category of results. EKF divergence de nes the third category, and was evident in the ! = 90 case. The categories of the additional test cases can be found in Table 5.13. Table 5.13: Additional Results Details Test Case Performance i = 28:5 Similar performance rp ra = 800km=1000km EKF divergencer p ra = 500km=700km Similar performance = 0 EKF divergence = 90 Similar performance ! = 45 Apparent divergence = 90 Similar performance = 270 Apparent divergence 58 Chapter 6 Conclusions and Future Work A simulation was built to compare the performance of an Extended Kalman Filter and an Iterated Extended Kalman Filter for a piece of debris with poor apriori data, inter- mittent, and innacurate observations. The altitude, inclination, longitude of the ascending node, argument of periapsis, true anomaly, and measurement covariance were varied to gen- erate fteen cases of initial conditions for gaining the performance data for the two ltering methods. With a constant ground network to generate observations, the longitude of the ascending node, argument of periapsis, and true anomaly generated observation data sets with varying measurement outage times, which tested the ability of the lters to recover following an outage. In many situations, the covariance matrix for the EKF became very small due to the low observation covariance and lack of process noise. This lead to smugness, or lter divergence due to the covariance matrix becoming small, in the EKF, an issue that was not antici- pated during the problem formulation. The IEKF overcame this issue due to the process of calculating the covariance matrix from processing each measurement multiple times. The performance of the IEKF generally allowed for a more accurate estimation of the state of the object with fewer measurements. However, this increase in accuracy occurred with an increase in computing cost. On average, the IEKF uses triple the processing needed by the EKF to produce an accuracy increase of at least an order of magnitude in the NMSE. The changes to the observation residual of the measurement without a change to the given observation residual produced unexpected results. Since the IEKF processes a mea- surement multiple times, it was expected that an increase in the measurement noise to increase the error in the estimate of the state produced by the IEKF. However, the IEKF 59 continued to produce a more accurate estimate of the state than the EKF in this situation due to the ability to continue processing a measurement until local convergence. To enhance this study in the future, many changes could be made to the simulation. Including process noise and aerodynamic drag in the dynamical model, adding intermittent measurements based on the attitude of the debris, and incorporating a detailed analysis of the computing power required are a few changes that could be incorporated to extend this study. The addition of process noise would test the ability of the IEKF to lter unknown accelerations. Aerodynamic drag calculations would allow LEO conditions to be better tested. Intermittent measurements based on the attitude of the debris would allow for the modeling of debris that might not produce a radar echo due to an orientation, which is an issue with lower power radar stations or small debris. A detailed analysis of the increase in computing power would allow for the formation of a ratio of increased power needed to the increase in accuracy to be de ned. This de nition would allow for an in depth study of all higher order lters for orbit determination to be compared more easily. 60 Bibliography [1] Orbital Debris: A Technical Assessment, National Academy, 1995. [2] Kessler, D. J. and Cour-Palais, B. G., \Collision Frequency of Arti cial Satellites: The Creation of a Debris Belt," Orbital Debris Quarterly News, Vol. 83, 1978. [3] Portree, D. S. and Joseph P. Loftus, J., \Orbital Debris: A Chronology," NASA paper nasa/tp-1999-208856, 1999. [4] Limiting Future Collision Risk to Spacecraft: An Assessment of NASA?s Meteoroid and Orbital Debris Programs, National Academies, 2011. [5] Liou, J.-C. and Shoots, D., \Fengyun-1C Debris: One Year Later," Orbital Debris Quarterly News, Vol. 12, No. 1, 2008. [6] Liou, J.-C. and Johnson, N. L., \A Sensitivity Study of the E ectiveness of Active Debris Removal in LEO," Acta Astronautica, Vol. 64, 2009, pp. 236 { 243. [7] Kessler, D. J., Johnson, N. L., Liou, J.-C., and Matney, M., \The Kessler Syndrome: Implications to Future Space Operations," AAS paper 10-016, 2010. [8] Sibert, D., Borgeson, D., Peterson, G., Jenkin, A., and Sorge, M., \Operational Impact of Improved Space Tracking on Collision Avoidance in the Future LEO Space Debris Environment," USAF space command report, 2010. [9] Thomson, W. T., Introduction to Space Dynamics, Wiley, 1961. [10] Montenbruck, O. and Gill, E., Satellite Orbits: Models, Methods and Applications, Springer, 2005. [11] Chobotov, V. A., Orbital Mechanics, American Institute of Aeronautics and Astronau- tics, Inc., 2002. [12] Tapley, B. D., Shutz, B. E., and Born, G. H., Statistical Orbit Determination, Elsevier Academic Press, 2004. [13] Gelb, A., Applied Optimal Estimation, M.I.T. Press, 1974. [14] Jazwinski, A. H., Stochastic Processes and Filtering Theory, Academic Press, 1970. [15] Simon, D., Optimal State Estimation: Kalman, H [in nity] and Nonlinear Approaches, Wiley, 2006. 61 [16] Kelso, T., \Validation of SGP4 and IS-GPS-200D Against GPS Precision Ephemerides," AAS paper 07-127, 2007. [17] Vallado, D. A. and McClain, W. D., Fundamentals of Astrodynamics and Applications, Springer, 2007. 62 Appendices 63 Appendix A Derivation of J2 Equations of motion In chapter 2, the acceleration due to the J2 perturbation was discussed. The full deriva- tion of the J2 acceleration in ECI components will be presented here, which will begin with Equation 2.7 and Equation 2.8, and end with Equation 2.9, which are also presented below as Equation A.1, Equation A.2, and Equation A.3 respectively for completion. UJ2 = r J2 re r 2 3 2 sin( ) 1 2 (A.1) r =rUJ2 (A.2) rJ2 = 3 2 J2r 2 e r7 2 66 66 4 X1 (5X23 r2) X2 (5X23 r2) X3 (5X23 3r2) 3 77 77 5 (A.3) Expanding the gradient operator in Equation A.2 yields Equation A.4, which will be evaluated individually by component. rUJ2 = @UJ2@X 1 ^I + @UJ2 @X2 ^J + @UJ2 @X3 ^K (A.4) The de nition of r in the ECI coordinate system and the norm of r, Equation A.5 and Equation A.6 respectfully, will be used to transform a portion of Equation A.1 from spherical 64 coordinate frame to the ECI coordinate frame. A substitution for sin( ) will complete the transformation, and is found from Figure 2.5. This is presented in Equation A.7. r = [X1 X2 X3] 2 66 66 4 I J K 3 77 77 5 (A.5) r = X21 +X22 +X23 12 (A.6) sin( ) = X3r (A.7) Substitute Equation A.7 into Equation A.1 and distribute. UJ2 = 3 2 J2r 2 eX 2 3 r5 + 1 2 J2r 2 e r3 (A.8) Substitute Equation A.6 into Equation A.8 and invert the denominators. This now gives UJ2 only in terms of constants and the components of the ECI position vector. UJ2 = 32 J2r2eX23 X21 +X22 +X23 52 + 12 J2r2e X21 +X22 +X23 32 (A.9) Combine the constants J2 and r2e to form a new constant C and substitute into Equation A.9. C = J2r2e (A.10) UJ2 = 32CX23 X21 +X22 +X23 52 + 12C X21 +X22 +X23 32 (A.11) 65 ^I Component (rUJ2)i = @@X 1 32CX23 X21 +X22 +X23 52 + 12C X21 +X22 +X23 32 (A.12) Perform the chain rule on each portion of the sum. (rUJ2)i = 154 (2X1)CX23 X21 +X22 +X23 72 + 34 (2X1)C X21 +X22 +X23 52 (A.13) Simplify, obtain a common denominator, factor the numerator and substitute for C. (rUJ2)i = 15 2 CX2 3X1 r7 + 32 CX 1 r5 (A.14) (rUJ2)i = 15 2 CX2 3X1 r7 + 32 CX 1r2 r7 (A.15) (rUJ2)i = 3 2 J2r 2 eX1 (5X 2 3 r 2) r7 (A.16) ^J Component (rUJ2)j = @@X 2 32CX23 X21 +X22 +X23 52 + 12C X21 +X22 +X23 32 (A.17) Perform the chain rule on each portion of the sum. 66 (rUJ2)j = 15 4 (2X2)CX23 X21 +X22 +X23 72 + 34 (2X2)C X21 +X22 +X23 52 (A.18) Simplify, obtain a common denominator, factor the numerator and substitute for C. (rUJ2)j = 15 2 CX2 3X2 r7 + 32 CX 2 r5 (A.19) (rUJ2)j = 15 2 CX2 3X2 r7 + 32 CX 2r2 r7 (A.20) (rUJ2)j = 3 2 J2r 2 eX2 (5X 2 3 r 2) r7 (A.21) ^K Component (rUJ2)k = @@X 3 32CX23 X21 +X22 +X23 52 + 12C X21 +X22 +X23 32 (A.22) Perform the chain rule on each portion of the sum (broken into the three terms for clarity). (rUJ2)k = 32 (2X3)C X21 +X22 +X23 52 (A.23) + 15 4 (2X3)CX23 X21 +X22 +X23 72 (A.24) 67 + 34 (2X3)C X21 +X22 +X23 52 (A.25) Sum Equations A.23 - A.25 and simplify, obtain a common denominator, factor the numerator and substitute for C. (rUJ2)k = 15 2 CX2 3X3 r7 9 2 CX 3 r5 (A.26) (rUJ2)k = 15 2 CX2 3X3 r7 9 2 CX 3r2 r7 (A.27) (rUJ2)k = 3 2 J2r 2 eX3 (5X 2 3 3r 2) r7 (A.28) Result Combining Equation A.16, Equation A.21, and Equation A.28 yields the equations of motion for the J2 perturbation in the ECI coordinate frame. aJ2 = 3 2 J2r 2 e r7 2 66 66 4 X1 (5X23 r2) X2 (5X23 r2) X3 (5X23 3r2) 3 77 77 5 (A.29) 68 Appendix B Additional Results B.1 Variations on Inclination B.1.1 Inclination !i = 28:5 0 100 200 300 400?5 0 5 10 15 20 Number of Observations Processed ? ? Ref line Figure B.1: i = 28:5 : vs. Number of Observations Processed 69 0 100 200 300 400 ?200 ?100 0 100 200 Number of Observations Processed Residual ? meters EKF IEKF Figure B.2: i = 28:5 : Residuals vs. Number of Observations Processed 1800 1820 1840 1860 1880 1900 ?200 ?100 0 100 200 Time Elapsed ? seconds Residual ? meters EKF IEKF Figure B.3: i = 28:5 : Residuals Before Largest Outage vs. Time 70 5550 5600 5650 5700 5750 5800 ?200 ?100 0 100 200 Time Elapsed ? seconds Residual ? meters EKF IEKF Figure B.4: i = 28:5 : Residuals After Largest Outage vs. Time 1800 1820 1840 1860 1880 19000 0.5 1 1.5 2 2.5 3x 10?3 Time Elapsed ? seconds ? n EKF IEKF Figure B.5: i = 28:5 : NMSE Before Largest Outage vs. Time 71 5550 5600 5650 5700 5750 58000 0.5 1 1.5 2 2.5 3x 10?3 Time Elapsed ? seconds ? n EKF IEKF Figure B.6: i = 28:5 : NMSE After Largest Outage vs. Time Table B.1: i = 28:5 : Intermediate Output 101 observations (1900 s) - Before 143 observations (5790 s) - After State (EKF) (IEKF) (EKF) (IEKF) X(m) 696.7952 447.1226 169.4131 78.3011 Y(m) 50.5685 15.1300 476.5338 217.3360 Z(m) 94.0815 71.4660 461.5652 215.0349 Vx(m/s) 1.9554 1.2306 0.5576 0.2643 Vy(m/s) 0.2796 0.1241 0.0613 0.0128 Vz(m/s) 1.1402 0.8525 0.0168 0.0021 yk(m) 82.3730 44.7652 55.6191 23.8096 NMSE 7.4027E-004 5.1267E-004 4.4888E-004 2.0883E-004 72 Table B.2: i = 28:5 : Output for All Observations (307 - 7830 s) Processed State (EKF) (IEKF) X(m) 3.1548 2.4580 Y(m) 21.5692 10.5556 Z(m) 87.5892 45.5551 Vx(m/s) 0.0584 0.0299 Vy(m/s) 0.0095 0.0046 Vz(m/s) 0.2404 0.1244 yk(m) 17.1910 7.7568 NMSE 9.5393E-005 4.9467E-005 B.2 Variations on Altitude B.2.1 Altitudes ! rpra = 800km=1000km 0 100 200 300 4000 500 1000 1500 2000 2500 3000 Number of Observations Processed ? ? Ref line Figure B.7: rpra = 800km=1000km : vs. Number of Observations Processed 73 0 100 200 300 400?8 ?6 ?4 ?2 0 2 4 6 8x 105 Number of Observations Processed Residual ? meters EKF IEKF Figure B.8: rpra = 800km=1000km : Residuals vs. Number of Observations Processed 1.321.331.341.351.361.371.38x 104?3 ?2 ?1 0 1 2 3x 105 Time Elapsed ? seconds Residual ? meters EKF IEKF Figure B.9: rpra = 800km=1000km : Residuals Before Largest Outage vs. Time 74 1.861.871.881.89 1.9 1.911.92x 104?3 ?2 ?1 0 1 2 3x 105 Time Elapsed ? seconds Residual ? meters EKF IEKF Figure B.10: rpra = 800km=1000km : Residuals After Largest Outage vs. Time 1.321.331.341.351.361.371.38x 1040 0.2 0.4 0.6 0.8 1 Time Elapsed ? seconds ? n EKF IEKF Figure B.11: rpra = 800km=1000km : NMSE Before Largest Outage vs. Time 75 1.861.871.881.89 1.9 1.911.92x 1040 0.2 0.4 0.6 0.8 1 Time Elapsed ? seconds ? n EKF IEKF Figure B.12: rpra = 800km=1000km : NMSE After Largest Outage vs. Time Table B.3: rpra = 800km=1000km: Intermediate Output 103 observations (3180 s) - Before 193 observations (13735 s) - After State (EKF) (IEKF) (EKF) (IEKF) X(m) 30803.5776 843.7531 10827.2371 68.2194 Y(m) 21335.0050 2114.0440 19576.7362 58.0021 Z(m) 22438.5898 898.8687 3489.0570 328.0003 Vx(m/s) 34.8935 0.1445 110.3779 0.0079 Vy(m/s) 12.7174 3.5697 84.6951 0.1300 Vz(m/s) 17.0665 2.3051 101.5724 0.1981 yk(m) 126.8850 4.5325 22496.07682 0.0034 NMSE 6.0557E-002 3.9742E-003 2.4168E-001 4.4422E-004 76 Table B.4: rpra = 800km=1000km: Output for All Observations (387 - 20195 s) Processed State (EKF) (IEKF) X(m) 669004.4833 483.2484 Y(m) 580538.8523 86.5897 Z(m) 701111.6741 1743.8944 Vx(m/s) 73.9316 0.9436 Vy(m/s) 991.5891 1.1375 Vz(m/s) 1289.7717 1.4730 yk(m) 82558.1852 102.4162 NMSE 1.7505E+000 1.9810E-003 B.2.2 Altitudes ! rpra = 500km=700km 0 100 200 300 400?5 0 5 10 15 Number of Observations Processed ? ? Ref line Figure B.13: rpra = 500km=700km : vs. Number of Observations Processed 77 0 100 200 300 400?50 0 50 Number of Observations Processed Residual ? meters EKF IEKF Figure B.14: rpra = 500km=700km : Residuals vs. Number of Observations Processed 1550160016501700175018001850 ?10 ?5 0 5 10 Time Elapsed ? seconds Residual ? meters EKF IEKF Figure B.15: rpra = 500km=700km : Residuals Before Largest Outage vs. Time 78 2700 2800 2900 3000 3100 ?10 ?5 0 5 10 Time Elapsed ? seconds Residual ? meters EKF IEKF Figure B.16: rpra = 500km=700km : Residuals After Largest Outage vs. Time 15501600165017001750180018500 0.2 0.4 0.6 0.8 1x 10?5 Time Elapsed ? seconds ? n EKF IEKF Figure B.17: rpra = 500km=700km : NMSE Before Largest Outage vs. Time 79 2700 2800 2900 3000 31000 0.2 0.4 0.6 0.8 1x 10?5 Time Elapsed ? seconds ? n EKF IEKF Figure B.18: rpra = 500km=700km : NMSE After Largest Outage vs. Time Table B.5: rpra = 500km=700km: Intermediate Output 230 observations (1845 s) - Before 293 observations (3090 s) - After State (EKF) (IEKF) (EKF) (IEKF) X(m) 3.9550 2.5127 2.9155 1.3144 Y(m) 1.8746 0.9790 5.2372 2.1131 Z(m) 0.3748 0.3505 1.2955 0.4534 Vx(m/s) 0.0039 0.0030 0.0001 0.0001 Vy(m/s) 0.0027 0.0004 0.0069 0.0025 Vz(m/s) 0.0010 0.0005 0.0019 0.0010 yk(m) 4.1478 4.2310 1.7809 1.6388 NMSE 2.2612E-006 1.5855E-006 4.9449E-006 1.9766E-006 80 Table B.6: rpra = 500km=700km: Output for All Observations (352 - 6125 s) Processed State (EKF) (IEKF) X(m) 1.2736 0.9745 Y(m) 0.5712 0.0809 Z(m) 0.0420 0.1487 Vx(m/s) 0.0003 0.0007 Vy(m/s) 0.0008 0.0004 Vz(m/s) 0.0001 0.0001 yk(m) 1.3117 1.3566 NMSE 8.1530E-007 6.2982E-007 B.3 Variations on Longitude of the Ascending Node B.3.1 Longitude of the Ascending Node ! = 0 0 100 200 300 4000 2 4 6 8 10 12x 104 Number of Observations Processed ? ? Ref line Figure B.19: = 0 : vs. Number of Observations Processed 81 0 100 200 300 400?6 ?4 ?2 0 2 4 6x 106 Number of Observations Processed Residual ? meters EKF IEKF Figure B.20: = 0 : Residuals vs. Number of Observations Processed 2500 2550 2600 2650 2700?4000 ?2000 0 2000 4000 Time Elapsed ? seconds Residual ? meters EKF IEKF Figure B.21: = 0 : Residuals Before Largest Outage vs. Time 82 1.671.6751.681.6851.691.6951.7x 104?1 ?0.5 0 0.5 1x 106 Time Elapsed ? seconds Residual ? meters EKF IEKF Figure B.22: = 0 : Residuals After Largest Outage vs. Time 2500 2550 2600 2650 27000 0.1 0.2 0.3 0.4 0.5 Time Elapsed ? seconds ? n EKF IEKF Figure B.23: = 0 : NMSE Before Largest Outage vs. Time 83 1.671.6751.681.6851.691.6951.7x 1040 0.5 1 1.5 2 2.5 Time Elapsed ? seconds ? n EKF IEKF Figure B.24: = 0 : NMSE After Largest Outage vs. Time Table B.7: = 0 : Intermediate Output 33 observations (2695 s) - Before 79 observations (16970 s) - After State (EKF) (IEKF) (EKF) (IEKF) X(m) 17077.7443 93.4052 123554.1400 47.8484 Y(m) 65766.6982 1981.8738 37632.0560 70.8962 Z(m) 8785.4193 426.7605 54013.8444 90.1187 Vx(m/s) 69.1645 2.5717 5830.9754 0.0089 Vy(m/s) 402.4925 15.1743 1819.2047 0.1625 Vz(m/s) 266.2634 11.8059 7545.5126 0.1870 yk(m) 276.0592 4.9075 72437.4711 2.3038 NMSE 4.0674E-001 1.3350E-002 2.2345E+000 7.4732E-005 84 Table B.8: = 0 : Output for All Observations (311 - 29325 s) Processed State (EKF) (IEKF) X(m) 4816071.3633 84.0402 Y(m) 2461167.8588 933.7677 Z(m) 3802048.3414 1271.3336 Vx(m/s) 1871.2821 0.4234 Vy(m/s) 3623.3970 1.4689 Vz(m/s) 5647.5764 2.1060 yk(m) 5243457.7008 507.2702 NMSE 4.0350E+000 1.3556E-003 B.3.2 Longitude of the Ascending Node ! = 90 Longitude of the Ascending Node 90 degrees 0 100 200 300 400?50 0 50 100 150 200 250 300 Number of Observations Processed ? ? Ref line Figure B.25: = 90 : vs. Number of Observations Processed 85 0 100 200 300 400?500 0 500 Number of Observations Processed Residual ? meters EKF IEKF Figure B.26: = 90 : Residuals vs. Number of Observations Processed 7250730073507400745075007550?200 ?100 0 100 200 Time Elapsed ? seconds Residual ? meters EKF IEKF Figure B.27: = 90 : Residuals Before Largest Outage vs. Time 86 1.2 1.205 1.21 1.215 1.22x 104?200 ?100 0 100 200 Time Elapsed ? seconds Residual ? meters EKF IEKF Figure B.28: = 90 : Residuals After Largest Outage vs. Time 72507300735074007450750075500 0.5 1 1.5 2 2.5x 10?5 Time Elapsed ? seconds ? n EKF IEKF Figure B.29: = 90 : NMSE Before Largest Outage vs. Time 87 1.2 1.205 1.21 1.215 1.22x 1040 1 2 x 10?4 Time Elapsed ? seconds ? n EKF IEKF Figure B.30: = 90 : NMSE After Largest Outage vs. Time Table B.9: = 90 : Intermediate Output 256 observations (7525 s) - Before 289 observations (12175 s) - After State (EKF) (IEKF) (EKF) (IEKF) X(m) 17.1138 2.2958 48.9182 3.0406 Y(m) 9.8902 0.4304 216.8610 0.9163 Z(m) 47.6340 1.9219 9.0753 2.7026 Vx(m/s) 0.0438 0.0100 0.1633 0.0056 Vy(m/s) 0.0226 0.0001 0.0819 0.0006 Vz(m/s) 0.0124 0.0075 0.2202 0.0044 yk(m) 6.6001 1.7129 145.7722 3.2522 NMSE 2.0006E-005 3.8020E-006 2.2582E-004 5.2616E-006 88 Table B.10: = 90 : Output for All Observations (323 - 12915 s) Processed State (EKF) (IEKF) X(m) 68.8861 1.4080 Y(m) 141.9399 0.7742 Z(m) 131.4218 0.9827 Vx(m/s) 0.1679 0.0061 Vy(m/s) 0.1082 0.0001 Vz(m/s) 0.1686 0.0052 yk(m) 114.3858 0.5844 NMSE 1.0019E-004 2.9362E-006 B.4 Variations on the Argument of Perigee B.4.1 Argument of Perigee !! = 45 0 100 200 300 400?5 0 5 10 15 20 Number of Observations Processed ? ? Ref line Figure B.31: ! = 45 : vs. Number of Observations Processed 89 0 100 200 300 400?1 ?0.5 0 0.5 1x 104 Number of Observations Processed Residual ? meters EKF IEKF Figure B.32: ! = 45 : Residuals vs. Number of Observations Processed 300 320 340 360 380 400?50 0 50 Time Elapsed ? seconds Residual ? meters EKF IEKF Figure B.33: ! = 45 : Residuals Before Largest Outage vs. Time 90 5100 5150 5200 5250 5300?2000 ?1000 0 1000 2000 Time Elapsed ? seconds Residual ? meters EKF IEKF Figure B.34: ! = 45 : Residuals After Largest Outage vs. Time 300 320 340 360 380 4000 0.005 0.01 0.015 0.02 0.025 Time Elapsed ? seconds ? n EKF IEKF Figure B.35: ! = 45 : NMSE Before Largest Outage vs. Time 91 5100 5150 5200 5250 53000 0.01 0.02 0.03 0.04 0.05 0.06 Time Elapsed ? seconds ? n EKF IEKF Figure B.36: ! = 45 : NMSE After Largest Outage vs. Time Table B.11: ! = 45 : Intermediate Output 19 observations (395 s) - Before 53 observations (5270 s) - After State (EKF) (IEKF) (EKF) (IEKF) X(m) 2724.5758 2579.9428 2986.3228 685.4223 Y(m) 3020.5290 2587.4315 8828.9447 1967.9265 Z(m) 1015.5812 1550.4209 7485.3752 1595.5966 Vx(m/s) 34.0881 31.3766 8.9811 0.9847 Vy(m/s) 3.0843 2.3657 4.0107 1.1805 Vz(m/s) 13.0670 13.8903 5.7072 3.7943 yk(m) 10.3421 5.3753 272.9182 72.8696 NMSE 2.4088E-002 2.2472E-002 4.4730E-003 1.1370E-003 92 Table B.12: ! = 45 : Output for All Observations (301 - 13425 s) Processed State (EKF) (IEKF) X(m) 2098.7609 214.1174 Y(m) 1943.8066 85.5751 Z(m) 1073.0904 120.8817 Vx(m/s) 3.6188 0.2484 Vy(m/s) 3.2936 0.1270 Vz(m/s) 1.7417 0.1291 yk(m) 129.9064 16.1376 NMSE 2.3449E-003 1.9357E-004 B.5 Variations on the True Anomaly B.5.1 True Anomaly ! = 90 True Anomaly 90 degrees 0 100 200 300 400?20 0 20 40 60 80 Number of Observations Processed ? ? Ref line Figure B.37: = 90 : vs. Number of Observations Processed 93 0 100 200 300 400?1000 ?500 0 500 1000 Number of Observations Processed Residual ? meters EKF IEKF Figure B.38: = 90 : Residuals vs. Number of Observations Processed 0 20 40 60 80 100 120?300 ?200 ?100 0 100 200 300 Time Elapsed ? seconds Residual ? meters EKF IEKF Figure B.39: = 90 : Residuals Before Largest Outage vs. Time 94 4600465047004750480048504900?300 ?200 ?100 0 100 200 300 Time Elapsed ? seconds Residual ? meters EKF IEKF Figure B.40: = 90 : Residuals After Largest Outage vs. Time 0 50 100 1500 0.02 0.04 0.06 0.08 0.1 Time Elapsed ? seconds ? n EKF IEKF Figure B.41: = 90 : NMSE Before Largest Outage vs. Time 95 4600 4700 4800 49000 0.02 0.04 0.06 0.08 0.1 Time Elapsed ? seconds ? n EKF IEKF Figure B.42: = 90 : NMSE After Largest Outage vs. Time Table B.13: = 90 : Intermediate Output 22 observations (110 s) - Before 70 observations (4860 s) - After State (EKF) (IEKF) (EKF) (IEKF) X(m) 1681.0471 583.7434 69.6992 77.4679 Y(m) 400.9202 16.4311 86.4232 63.8825 Z(m) 4112.8134 2594.4797 150.5687 120.9722 Vx(m/s) 7.1026 0.4560 0.0970 0.1351 Vy(m/s) 1.3260 1.1095 0.3631 0.0991 Vz(m/s) 12.1601 7.6690 1.2757 0.3698 yk(m) 0.8732 0.6893 31.4049 9.8177 NMSE 1.7805E-002 9.1377E-003 3.1994E-004 1.0862E-004 96 Table B.14: = 90 : Output for All Observations (301 - 16150 s) Processed State (EKF) (IEKF) X(m) 23.1059 7.0732 Y(m) 32.1170 6.0575 Z(m) 53.9966 26.2700 Vx(m/s) 0.0166 0.0002 Vy(m/s) 0.0214 0.0048 Vz(m/s) 0.0642 0.0337 yk(m) 4.7738 10.1642 NMSE 2.7934E-005 1.2870E-005 B.5.2 True Anomaly ! = 270 True Anomaly 270 degrees 0 100 200 300 400?5 0 5 10 15 20 25 Number of Observations Processed ? ? Ref line Figure B.43: = 270 : vs. Number of Observations Processed 97 0 100 200 300 400?500 0 500 Number of Observations Processed Residual ? meters EKF IEKF Figure B.44: = 270 : Residuals vs. Number of Observations Processed 2.3252.332.3352.342.3452.352.355x 104?500 0 500 Time Elapsed ? seconds Residual ? meters EKF IEKF Figure B.45: = 270 : Residuals Before Largest Outage vs. Time 98 3.565 3.57 3.575 3.58 3.585x 104?500 0 500 Time Elapsed ? seconds Residual ? meters EKF IEKF Figure B.46: = 270 : Residuals After Largest Outage vs. Time 2.3252.332.3352.342.3452.352.355x 1040 0.2 0.4 0.6 0.8 1x 10?3 Time Elapsed ? seconds ? n EKF IEKF Figure B.47: = 270 : NMSE Before Largest Outage vs. Time 99 3.565 3.57 3.575 3.58 3.585x 1040 0.2 0.4 0.6 0.8 1x 10?3 Time Elapsed ? seconds ? n EKF IEKF Figure B.48: = 270 : NMSE After Largest Outage vs. Time Table B.15: = 270 : Intermediate Output 43 observations (2,820 s) - Before 61 observations (2,910 s) - After State (EKF) (IEKF) (EKF) (IEKF) X(m) 14039.59483 2908.0591 2122.2537 65.95151 Y(m) 10357.9379 206.7330 60.4991 29.0566 Z(m) 746.5737 6057.0141 6353.7477 210.5541 Vx(m/s) 9.5159 5.9085 7.2863 0.2481 Vy(m/s) 1.3419 0.5698 0.6326 0.2020 Vz(m/s) 5.3723 0.7508 0.8307 0.1025 yk(m) 110.1015 0.0421 11.4919 1.5893 NMSE 3.4027E-002 1.6058E-002 8.7018E-003 3.1182E-004 100 Table B.16: = 270 : Output for All Observations (326 - 41630 s) Processed State (EKF) (IEKF) X(m) 8.8143 0.8248 Y(m) 15.5944 8.2807 Z(m) 92.1830 45.0967 Vx(m/s) 0.0201 0.0364 Vy(m/s) 0.4172 0.2035 Vz(m/s) 0.4210 0.1500 yk(m) 13.8092 4.3864 NMSE 2.1635E-004 8.2804E-005 101