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