Improvement of a Space Surveillance Tracking Analysis Tool

Ismail Ugur DUNDAR

Space Engineering, master's level (120 credits) 2018

Luleå University of Technology

Department of Computer Science, Electrical and Space Engineering

IMPROVEMENT OF A SPACE SURVEILLANCE TRACKING ANALYSIS TOOL

SUBMITTED BY

ISMAIL UGUR DUNDAR

Prepared at

Airbus Defence and Space GmbH 88039 Friedrichshafen, Germany

Examiner:

Prof. Dr. Sergio Montenegro, Julius–Maximilians–University Würzburg, Germany Dr. Victoria Barabash, Lulea University of Technology , Kiruna, Sweden

Supervisor:

Dr. Jens Utzmann

MSc. Oscar Rodriguez Fernandez

December 5, 2018

DISCLAIMER

"This project has been funded with support from the European Commission. This publication reflects the views only of the author, and the Commission cannot be held responsible for any use which may be made of the information contained therein."

i

Abstract

Since the beginning of space exploration, the amount of space debris has increased with the development of new space technologies. In fact, when a collision happens, new space debris are generated. Hence, collision risk between space debris and operational satellites rises. The purpose of a surveillance network system consists of the detection of space objects, their classification and their tracking. To avoid collisions, space debris objects’ orbit must be computed with sufficient accuracy.

The goal of this thesis is the improvement of a pre-existing Space Surveillance and Tracking Analysis Tool. The tool is able to simulate different observation scenarios for radar or optical observer, which can be space-based or ground-based. To enhance the orbit determination, an Extended Square Root Information Filter is implemented and incremented with a Smoother. Smoothers have been implemented for the existing filters as well, such as the Extended Kalman Filter and the Unscented Kalman Filter. A bias estimation method was added as part of the OD for all filter types.

Additionally, different outlier detection methods were implemented for the automatic detection of outliers within the measurement data. To find the optimum orbit determination interval, different scenarios were considered in LEO, MEO and GEO orbits. The methods were implemented and different scenarios for validation will be discussed. A wide discussion on the methods implementation and their validation on different scenarios is presented, together with a comparison of the orbit determination results with the other filters.

All the recently implemented features increase the efficiency of the tool to simulate the different

scenarios and enhance the tracking of space debris objects.

Acknowledgements

I would like to thank my supervisor at Airbus Defence and Space, MSc. Oscar Rodriguez Fernandez, for his patience, enthusiasm, motivation, immense knowledge and continuous support. His guidance helped me through all the research and writing of this thesis. My sincere gratitude also goes to Dr.

Jens Utzmann, who provided me an opportunity to join their team as a master thesis student.

I would also like to thank my university supervisors: Prof. Dr. Sergio Montenegro and Dr. Victoria Barabash, for their support and insightful comments.

iii

Declaration of Authorship

I declare that I completed this thesis on my own and that information which has been directly or indirectly taken from other sources has been noted as such. Neither this nor a similar work has been presented to an examination committee.

Name of Student

Friedrichshafen, December 5, 2018

. . . .

Contents

List of Figures . . . . xii

List of Tables . . . xiv

List of Abbreviations . . . . xv

List of Symbols . . . xvi

1 Introduction 1 1.1 Thesis Objectives . . . . 3

2 Theoretical Background 4 2.1 Coordinate Frames . . . . 4

2.1.1 Geocentric Celestial Coordinate System . . . . 4

2.1.2 International Terrestrial Reference Frame . . . . 5

2.2 Object State . . . . 6

2.2.1 Space-based Objects . . . . 6

2.2.2 Ground-based Observers . . . . 7

2.3 Observation Theory . . . . 7

2.3.1 Optical Observer . . . . 8

2.3.2 Radar Observer . . . . 8

2.4 Dynamic Model - Orbital Perturbation . . . . 9

2.4.1 Non-spherical Gravitational Field . . . . 9

2.4.2 Atmospheric Drag . . . . 10

2.4.3 Solar Radiation Pressure . . . . 10

2.4.4 Third-body Perturbations . . . . 10

2.5 Initial Orbit Determination . . . . 11

2.5.1 Gibbs Algorithm . . . . 11

2.5.2 Herrick-Gibbs Algorithm . . . . 11

2.6 Orbit Determination . . . . 11

2.6.1 The Weighted Least Squares Method . . . . 12

2.6.2 The Sequential Batched Least Squares Algorithm . . . . 13

2.6.3 The Extended Kalman Filter . . . . 13

2.6.4 Unscented Kalman Filter . . . . 14

2.7 Orbital Regions . . . . 16

v

3 SPOOK 18

3.1 SPOOK Structure . . . . 18

3.1.1 Simulation Layer . . . . 19

3.1.2 Interface Layer . . . . 20

3.1.3 Analysis Layer . . . . 21

4 Method 25 4.1 Bias Estimation . . . . 25

4.2 Extended Square Root Information Filter . . . . 26

4.2.1 Correction of Non-Positive Definite Matrices due to Numerical Error . . . . . 31

4.3 Smoother . . . . 31

4.3.1 EKF . . . . 32

4.3.2 UKF . . . . 33

4.3.3 ESRIF . . . . 33

4.4 Outlier Detection . . . . 33

4.4.1 Chi-Square Test . . . . 34

4.4.2 Measurement Variance Test . . . . 34

4.4.3 3σ Test . . . . 34

4.4.4 Data Editing . . . . 35

4.5 Orbit Determination Interval . . . . 35

4.5.1 Credibility Test . . . . 36

5 Evaluation and Discussion 38 5.1 Experiment . . . . 38

5.1.1 GEO Orbit - Object Scenario 1 . . . . 38

5.1.2 LEO Orbit - Object Scenario 2 . . . . 39

5.1.3 GEO Orbit - Object Scenario 3 . . . . 40

5.2 Bias estimation . . . . 41

5.3 Extended Square Root Information Filter . . . . 51

5.4 Smoother . . . . 57

5.4.1 EKF Smoother . . . . 58

5.4.2 ESRIF Smoother . . . . 62

5.4.3 UKF Smoother . . . . 66

5.5 Outlier Detection . . . . 70

5.5.1 EKF Outlier . . . . 71

5.5.2 ESRIF Outlier . . . . 74

5.5.3 UKF Outlier . . . . 76

5.6 Orbit Determination Interval . . . . 78

5.6.1 LEO Orbit . . . . 83

5.6.2 MEO Orbit . . . 107

5.6.3 GEO Orbit . . . 110

6 Conclusion 113

7 Future Work 114

A.4 WLS . . . 122

A.5 SBLS . . . 123

B Performance of Filters 125 B.1 Initial Covariance Propagation (ICP) - GEO Object Scenario 1 . . . 125

B.1.1 EKF . . . 125

B.1.2 ESRIF . . . 126

B.1.3 UKF . . . 126

B.2 Initial Covariance Propagation (ICP) - LEO Object Scenario 2 . . . 127

B.2.1 EKF . . . 127

B.2.2 ESRIF . . . 128

B.2.3 UKF . . . 128

B.3 Initial Orbit Determination (IOD) - GEO Object Scenario 1 . . . 129

B.3.1 EKF . . . 129

B.3.2 ESRIF . . . 130

B.3.3 UKF . . . 130

B.4 Initial Orbit Determination (IOD) - LEO Object Scenario 2 . . . 131

B.4.1 EKF . . . 131

B.4.2 ESRIF . . . 132

B.4.3 UKF . . . 132

B.5 Comparison of the Smoother Results with WLS . . . 133

vii

List of Figures

1.1 Distribution of space debris in the Earth orbit . . . . 1

1.2 Simulation of the collision before and during the first three hours . . . . 2

1.3 Number of objects per orbital region included in the TLE catalogue . . . . 2

2.1 Geocentric Celestial Reference Frame . . . . 5

2.2 The International Terrestrial Reference Frame . . . . 5

2.3 Keplerian Classical Orbital Elements . . . . 7

2.4 Angular measurements for optical observer: right ascension α and declination δ . . . 8

2.5 Angular measurement of the radar observer . . . . 9

3.1 SPOOK functional diagram . . . . 19

3.2 SPOOK visualisation tool . . . . 21

3.3 Main structure of SPOOK. . . . . 23

3.4 Structure of Orbit Determination. . . . . 24

4.1 ESRIF Algorithm Structure. . . . 30

4.2 Forward - Backward Filter . . . . 31

4.3 Filter Time Step . . . . 32

5.1 Errors of position and 3σ-sigma uncertainty from covariance of the orbit determination of the object for case 1. . . . . 42

5.2 Residuals of Orbit Determination of the object for case 1. . . . 43

5.3 Errors of position and 3σ-sigma uncertainty from covariance of the orbit determination of the object for case 2. . . . . 44

5.4 Residuals of Orbit Determination of the object for case 2. . . . 45

5.5 Errors of position and 3σ-sigma uncertainty from covariance of the orbit determination of the object for case 3. . . . . 46

5.6 Residuals of Orbit Determination of the object for case 3. . . . 47

5.7 Bias Estimation for case 3. . . . . 48

5.8 Errors of position and 3σ-sigma uncertainty from covariance of the orbit determination of the object for case 4. . . . . 49

5.9 Residuals of Orbit Determination of the object for case 4. . . . 50

5.10 Bias Estimation for case 4. . . . . 51

5.11 Errors of position and uncertainty of the orbit determination - ESRIF. . . . 52

5.12 Residuals of the Orbit Determination - ESRIF. . . . 53

5.13 Errors of position and uncertainty of the orbit determination - EKF. . . . 54

5.14 Residuals of the Orbit Determination - EKF. . . . 55

5.15 Errors of position and uncertainty of the orbit determination - UKF. . . . 56

5.23 ESRIF - Forward Filter. . . . . 62

5.24 ESRIF - Backward Filter. . . . . 63

5.25 ESRIF - Forward Filter and Backward Filter Together. . . . . 64

5.26 ESRIF - Difference of the Forward Filter and Backward Filter. . . . 64

5.27 ESRIF - Forward Filter Residual. . . . . 65

5.28 ESRIF - Backward Filter Residual. . . . . 66

5.29 UKF - Forward Filter of the Orbit Determination Error. . . . . 66

5.30 UKF - Backward Filter of the Orbit Determination Error. . . . . 67

5.31 UKF - Forward and Backward Filter Together. . . . 68

5.32 UKF - Difference of the Forward Filter and Backward Filter. . . . 68

5.33 UKF - Forward Filter of the Residual. . . . . 69

5.34 UKF - Forward Filter of the Residual with large scale. . . . . 69

5.35 UKF - Backward Filter of the Residual. . . . 70

5.36 EKF - Outlier Detection Algorithm not Active. . . . 71

5.37 EKF - Outlier Detection Method 1 Active . . . . 72

5.38 EKF - Outlier Detection Method 2 Active. . . . 72

5.39 EKF - Outlier Detection Method 3 Active. . . . 73

5.40 ESRIF - Outlier Detection Method 1 Active . . . . 74

5.41 ESRIF - Outlier Detection Method 2 Active. . . . . 75

5.42 ESRIF - Outlier Detection Method 3 Active. . . . . 76

5.43 UKF - Outlier Detection not Active. . . . 76

5.44 UKF - Outlier Detection Method 1 Active. . . . . 77

5.45 UKF - Outlier Detection Method 2 Active. . . . . 77

5.46 UKF - Outlier Detection Method 3 Active. . . . . 78

5.47 Orbit Determination Interval Scenario result part 1. . . . 79

5.48 Orbit Determination Interval Scenario result part 2. . . . 79

5.49 Orbit Determination - Object Number 1. . . . 80

5.50 Residual - Object Number 1. . . . . 80

5.51 Orbit Determination - Object Number 4. . . . 81

5.52 Residual - Object Number 4. . . . . 81

5.53 Orbit Determination - Object Number 13. . . . 82

5.54 Residual - Object Number 13. . . . 82

5.55 Sensor Type 1 - Consistency result in Region 1. . . . 85

5.56 Sensor Type 1 - WRMS result in Region 1. . . . 85

5.57 Sensor Type 1 - Consistency result in Region 2. . . . 86

5.58 Sensor Type 1 - WRMS result in Region 2. . . . 86

5.59 Sensor Type 1 - Consistency result in Region 3. . . . 87

5.60 Sensor Type 1 - WRMS result in Region 3. . . . 87

5.61 Sensor Type 1 - Consistency result in Region 4. . . . 88

5.62 Sensor Type 1 - WRMS result in Region 4. . . . 88

ix

5.63 Sensor Type 1 - Consistency result in Region 7. . . . 89

5.64 Sensor Type 1 - WRMS result in Region 7. . . . 89

5.65 Sensor Type 1 - Consistency result in Region 7. . . . 90

5.66 Sensor Type 1 - WRMS result in Region 7. . . . 90

5.67 Sensor Type 2 - Consistency result in Region 1. . . . 91

5.68 Sensor Type 2 - WRMS result in Region 1. . . . 91

5.69 Sensor Type 2 - Consistency result in Region 2. . . . 92

5.70 Sensor Type 2 - WRMS result in Region 2. . . . 92

5.71 Sensor Type 2 - Consistency result in Region 3. . . . 93

5.72 Sensor Type 2 - WRMS result in Region 3. . . . 93

5.73 Sensor Type 2 - Consistency result in Region 4. . . . 94

5.74 Sensor Type 2 - WRMS result in Region 4. . . . 94

5.75 Sensor Type 2 - Consistency result in Region 7. . . . 95

5.76 Sensor Type 2 - WRMS result in Region 7. . . . 95

5.77 Sensor Type 3 - Consistency result in Region 1. . . . 96

5.78 Sensor Type 3 - WRMS result in Region 1. . . . 96

5.79 Sensor Type 3 - Consistency result in Region 2. . . . 97

5.80 Sensor Type 3 - WRMS result in Region 2. . . . 97

5.81 Sensor Type 3 - Consistency result in Region 3. . . . 98

5.82 Sensor Type 3 - WRMS result in Region 3. . . . 98

5.83 Sensor Type 3 - Consistency result in Region 4. . . . 99

5.84 Sensor Type 3 - WRMS result in Region 4. . . . 99

5.85 Sensor Type 3 - Consistency result in Region 7. . . 100

5.86 Sensor Type 3 - WRMS result in Region 7. . . 100

5.87 Sensor Type 5 - Consistency result for Leo Orbit. . . 101

5.88 Sensor Type 5 - WRMS result for LEO Orbit. . . 101

5.89 Sensor Type 7 - Consistency result in Region 1. . . 102

5.90 Sensor Type 7 - WRMS result in Region 1. . . 102

5.91 Sensor Type 7 - Consistency result in Region 2. . . 103

5.92 Sensor Type 7 - WRMS result in Region 2. . . 103

5.93 Sensor Type 7 - Consistency result in Region 3. . . 104

5.94 Sensor Type 7 - WRMS result in Region 3. . . 104

5.95 Sensor Type 7 - Consistency result in Region 4. . . 105

5.96 Sensor Type 7 - WRMS result in Region 4. . . 105

5.97 Sensor Type 7 - Consistency result in Region 7. . . 106

5.98 Sensor Type 7 - WRMS result in Region 7. . . 106

5.99 Sensor Type 1 - Consistency result in Region 7 and 8. . . 107

5.100Sensor Type 1 - WRMS result in Region 7 and 8. . . 107

5.101Sensor Type 1 - Consistency result in Region 10. . . 108

5.102Sensor Type 1 - WRMS result in Region 10. . . 108

5.103Sensor Type 1 - Consistency result in Region 12. . . 109

5.104Sensor Type 1 - WRMS result in Region 12. . . 109

5.105Optical Observer - Consistency result for GEO Orbit. . . . 110

5.106Optical Observer - WRMS result for GEO Orbit. . . 110

5.107Sensor Type 1 - Consistency result in Region 12 and 13. . . 111

A.5 EKF - Result of the OD using Initial Orbit Determination. . . . 119

A.6 Result of the Bias Estimation - Right Ascension (0") and Declination (0") of bias values.119 A.7 EKF - Result of the OD using Initial Orbit Determination. . . . 120

A.8 Result of the Bias Estimation - Right Ascension (1") and Declination (1") of bias values.120 A.9 ESRIF - Result of the OD using Initial Covariance Propagation. . . . 120

A.10 Result of the Bias Estimation - Right Ascension (0") and Declination (0") of bias values.120 A.11 ESRIF - Result of the OD using Initial Covariance Propagation. . . . 121

A.12 Result of the Bias Estimation - Right Ascension (1") and Declination (1") of bias values.121 A.13 UKF - Result of the OD using Initial Covariance Propagation. . . . 121

A.14 Result of the Bias Estimation - Right Ascension (0") and Declination (0") of bias values.121 A.15 UKF - Result of the OD using Initial Covariance Propagation. . . . 122

A.16 Result of the Bias Estimation - Right Ascension (1") and Declination (1") of bias values.122 A.17 WLS - Result of the OD. . . . 122

A.18 Result of the Bias Estimation - Right Ascension (0") and Declination (0") of bias values.122 A.19 WLS - Result of the OD . . . 123

A.20 Result of the Bias Estimation - Right Ascension (1") and Declination (1") of bias values.123 A.21 SBLS - Result of the OD. . . 123

A.22 Result of the Bias Estimation - Right Ascension (0") and Declination (0") of bias values.123 A.23 SBLS - Result of the OD. . . 124

A.24 Result of the Bias Estimation - Range(0.5m), Azimuth (1"), Elevation (1") and Range-Rate (0") of bias values. . . 124

B.1 EKF - Result of the OD. . . 125

B.2 EKF - Result of the residual. . . 125

B.3 ESRIF - Result of the OD. . . . 126

B.4 ESRIF - Result of the residual. . . 126

B.5 UKF - Result of the OD. . . . 126

B.6 UKF - Result of the residual. . . 126

B.7 EKF - Result of the OD. . . 127

B.8 EKF - Result of the residual. . . 127

B.9 ESRIF - Result of the OD. . . . 128

B.10 ESRIF - Result of the residual. . . 128

B.11 UKF - Result of the OD. . . . 128

B.12 UKF - Result of the residual. . . 128

B.13 EKF - Result of the OD. . . 129

B.14 EKF - Result of the residual. . . 129

B.15 ESRIF - Result of the OD . . . 130

B.16 ESRIF - Result of the residual. . . 130

B.17 UKF - Result of the OD. . . . 130

B.18 UKF - Result of the residual. . . 130

B.19 EKF - Result of the OD. . . 131

xi

B.20 EKF - Result of the residual. . . 131

B.21 ESRIF - Result of the OD . . . 132

B.22 ESRIF - Result of the residual. . . 132

B.23 UKF - Result of the OD. . . . 132

B.24 UKF - Result of the residual. . . 132

B.25 WLS - Orbit Determination Result. . . . 133

B.26 WLS - Residual of the Measurements. . . 133

List of Tables

2.1 Definition of Classical Keplerian Orbital Elements . . . . 6

2.2 Orbital Region Classification . . . . 16

2.3 Augmented Orbital Region Classification for Orbit Determination Interval . . . . 17

3.1 Observer Types. . . . . 20

4.1 One of the combinations of the metrics parameter . . . . 36

4.2 Parameters of Metrics . . . . 36

5.1 Scenario 1 - Target Object and Observer . . . . 38

5.2 Scenario 1 Parameters . . . . 39

5.3 LEO Orbit - Target Object . . . . 39

5.4 Ground Radar Observer Location . . . . 39

5.5 Scenario 2 Parameters . . . . 40

5.6 Scenario 3 Parameters . . . . 40

5.7 Ground Optical Observer Location . . . . 40

5.8 Standard Deviation and Mean values for case 1. . . . 43

5.9 Standard Deviation and Mean values for case 2. . . . 45

5.10 Standard Deviation and Mean values for case 3. . . . 48

5.11 Standard Deviation and Mean values for case 4. . . . 50

5.12 Standard Deviation and Mean values. . . . . 57

5.13 RMS Error of OD and RMS Error of the Smoother OD. . . . 70

5.14 Outliers’ time amount in the measurement. . . . 71

5.15 Ground Observer Location. . . . 83

5.16 Sensor Type 1 - Combination Parameters. . . . 84

5.17 Sensor Type 1 - Summary of the ODI. . . . 90

5.18 Sensor Type 2 - Combination Parameters. . . . 91

5.19 Sensor Type 2 - Summary of the ODI. . . . 95

5.20 Sensor Type 3 - Combination Parameters. . . . 96

5.21 Sensor Type 3 - Summary of the ODI. . . 100

5.22 Sensor Type 5 - Combination Parameters. . . 101

5.23 Sensor Type 7 - Combination Parameters. . . 102

5.24 Sensor Type 7 - Summary of the ODI. . . 106

5.25 Orbit Determination Interval for Leo Orbit. . . 107

B.1 RMS Error of OD and Mean value of the Residual using ICP for Geo orbit object Scenario. . . 127

B.2 RMS Error of OD and Mean value of the Residual using ICP for Leo orbit object Scenario. . . 129

xiii

B.3 RMS Error of OD and Mean value of the Residual using IOD for Geo orbit object Scenario. . . 131 B.4 RMS Error of OD and Mean value of the Residual using IOD for Leo orbit object

Scenario. . . 133

B.5 RMS Error of OD and RMS Error of the Smoother OD. . . 134

List of Abbreviations

BL Backward Filter

BC Ballistic Coefficient

COPUOS Committee on the Peaceful Uses of Outer Space

ESA European Space Agency EKF Extended Kalman Filter

ESRIF Extended Square Root Information Fil- ter

FoV Field of View

FL Forward Filter

GCRF Geocentric Celestial Coordinate System GEO Geostationary Earth Orbit

GPS Global Positioning System IOD Initial Covariance Propagation IOD Initial Orbit Determination

ICRF International Celestial Reference Frame IERS International Earth Rotation Service ISS International Space Station

ITRS International Terrestrial Reference Frame

LEO Low Earth Orbit

MSIS Mass Spectrometer Incoherent Scatter

MEO Medium Earth Orbit

NEES Normalized Estimation Error Squared NIS Normalized Innovation Squared ODI Orbit Determination Interval OD Orbit Determination

RTN Radial Tangential Normal RTS Rauch Tung Striebel RMS Root Mean Square

SBLS Sequential Batch Least Square SRP Solar Radiation Pressure

SPOOK Space Object Observations and Kalman Filtering

SST Space Surveillance Tracking SRIF Square Root Information Filter TLE Two Line Element

UKF Unscented Kalman Filter

USKF Unscented Schmidt Kalman Filter WLS Weighted Least Square

WRMS Weighted Root Mean Square WGS84 World Geodetic System 1984

xv

List of Symbols

A D Object’s cross-sectional area

" arcsec

α Semi-major Axis

β Azimuth

χ Sigma Points

χ ² Chi-Square Distribution δ Declination

˙

ρ Slant Range-Rate λ Longitude e Eccentricity i Inclination v True Anomaly

w Argument of periapsis µ Gravitational Parameter

Ω Right Ascension of the Ascending Node φ Latitude

Φ State Transition Matrix ρ Slant Range

σ Standard Deviation τ Truncation Error α Right Ascension bj ~ Residual

~ r Position Vector

~ v Velocity Vector

~a Drag Atmospheric Drag Acceleration

~

v t,winds Thermospheric Wind

~

w _⊕ Earth’s Rotational Speed A j Partial Derivative Matrix C D Drag Coefficient

C R Reflectivity Coefficient D Smoother Kalman Gain

el Elevation

F Partial Derivative of the state rates matrix

h Altitude

H j Observation Matrix

K Kalman Gain MAtrix

k Time Step

m obj Mass of Object P Covariance Matrix p SRP Solar Radiation Pressure Q Process Noise Matrix r Spacecraft’s Position R Measurement Noise Matrix R Information Matrix S k Eigen Vector U Gravity Field W Weighted Matrix

y Measurement

Chapter 1

Introduction

Since the beginning of Space activities, man-made space debris have been created as a result of more than 5400 rocket launches and more than 500 on-orbit fragmentation events [1]. Space activities have grown substantially, and space exploration will increase further in the future. Because of lack of the space laws, several organizations have emerged, which are the United Nations Committee on the Peaceful Uses of Outer Space (COPUOS) and the Inter-Agency Space Debris Coordination Committee. The main purpose of these committees is to develop guidelines to decrease the amount of new space debris, upper stage rockets and defunct satellites [2]. Space debris is defined by IASDCC as follows:

"Space debris are all man made objects, including fragments and elements thereof, in Earth orbit or reentering the atmosphere, that are nonfunctional"[3]

Figure 1.1: Distribution of space debris in the Earth orbit [4].

More than 23.000 objects are currently catalogued by the U.S Space Surveillance Network [5].

Nevertheless, the total real object number can only be estimated by a scientific model. According to the European Space Agency (ESA) scientific model, an estimation of space debris in 2013, 29.000 objects are larger than 10 cm in size, 670.000 objects are between 1 cm to 10 cm and 170 million objects are between 1 mm to 1 cm. These sizes of objects are a problem for an operational spacecraft.

1. INTRODUCTION 1

For instance, in case of collision with a 10-cm object, a catastrophic fragmentation of a spacecraft will be unavoidable [6]. Larger than 10 cm in size might cause catastrophic breakup of a spacecraft.

Objects with diameter of 1 cm to 10 cm could create significant damage to a spacecraft and possible mission failure. Objects with diameter of 1 mm to 1 cm could cause component damage, craters, spallations and degradation od spacecraft surfaces [7]. Micrometeoroids have an average velocity and density of 20 km/sec and 0.5 g/cm ³ respectively whereas space objects average velocity and density are 9 to 10 km/sec and 2.8 g/cm ³ depending on the altitude [8]. Micrometeoroids could impact penetration of the spacecraft, RF transients, sudden electrical disturbances and gyro-destabilizations [7].

Two main events in history lead to the creation of a large number of new space debris in Earth orbit. One of them is when China destroyed Fengyun-1C in an anti-satellite weapon test in 2007 with a ballistic missile. The destruction created more than 3000 new objects [9]. The second is when a collision happened between Cosmos 2251 and Iridium 33 in 2009. These satellites are an inactive Russian communication satellite and an active satellite operated by the U.S. This collision created almost 2,000 space debris objects [10].

Figure 1.2: Simulation of the collision before and during the first three hours [10].

Figure 1.2 shows the two satellites at the time of the collision (left), the position of two satellites

10 minutes after the collision (center) and space debris clouds three hours later (right) [10].

laser observation from space or ground. Space objects’ orbital elements can be estimated from the observation measurements. This is known as orbit determination and it will be the major subject of this thesis.

1.1 Thesis Objectives

Space Object Observations and Kalman Filtering (SPOOK) is a software tool that has been being developed by Airbus Defence and Space. This is a highly configurable and versatile tool for simula- tion, detection as well as producing orbit determination of Earth-orbiting objects using either real or synthetic data. This tool has the ability to predict and detect the movement of the Earth-orbiting objects. It has the ability to use different observers (ground-based and space-based observers) and different sensors (radar and optical sensors). In order to perform orbit determination for space objects, several techniques such as Extended Kalman Filter, Weighted Least Square, Sequential Batch Least Square, Unscented Kalman Filter and Unscented Schmidt Kalman Filter can be used.

The objective of this thesis is to enhance this tool, focusing on the Orbit Determination aspect. The thesis works are defined as follows:

• Bias Estimation in Orbit Determination

• Implementation and validation of an ESRIF (Extended Square Root Information Filter)

• Implementation and validation of filter Smoother (EKF, UKF, ESRIF)

• Implementation and validation of Outlier Detection Algorithms

• Selection of optimum data interval for OD in LEO, MEO and GEO orbits

• Comparison of performance between different OD methods

In Chapter 2, the theoretical background is provided which is based on the orbit determina- tion problem. In Chapter 3, SPOOK characteristics and its situation of the beginning of the project are described. The methods to enhance the orbit determination are described in Chapter 4. Subsequently, the simulation parameters are defined, and the methods’ results are evaluated and discussed with different observer and orbit simulation in Chapter 5. Afterwards, Chapter 6 summarizes the achievements accomplished during the master thesis. Finally, Chapter 7 gives an aspect for eventual future enhancements of the SPOOK.

1. INTRODUCTION 3

Chapter 2

Theoretical Background

This section introduces the theoretical background, which is necessary to understand the Orbit Determination process. Some brief descriptions of astrodynamics fundamentals will be given in this chapter.

2.1 Coordinate Frames

There are two main reference system which are used: The Geocentric Celestial Coordinate Reference Frame (GCRF) and the International Terrestrial Reference Frame (ITRF). These frames are introduced briefly in this chapter.

2.1.1 Geocentric Celestial Coordinate System

The Geocentric Celestial Coordinate System is one of the classical inertial coordinate system for the Earth. It is the geocentric counterpart of the International Celestial Reference Frame (ICRF). This coordinate frame has been the accepted celestial reference frame for the International Earth Rotation Service (IERS) since January 1, 1997 [12]. This reference frame can be visualized as the first f 1

axis aligned with vernal equinox and the third f 3 axis aligned with Earth’s north pole (alignment

is exact only for the date 1 January 2000) [13]. Figure 2.2 indicates the GCRF and an arbitrary

position vector to indicate right ascension and declination [13]. GCRF is defined using the J2000

date (Julian date for 1 January 2000 at 12:00). In SPOOK, the GCRF is the main reference frame.

Figure 2.1: Geocentric Celestial Reference Frame [13].

2.1.2 International Terrestrial Reference Frame

A geocentric coordinate frame fixed to the rotating Earth. This coordinate system’s origin is at the center of the Earth [12]. The X-axis directs to the reference meridian, Z-axis directs to the reference pole and Y-axis completes the right-handed system.

More detail about transformation ITRF to GCRF can be found in [12].

Figure 2.2: The International Terrestrial Reference Frame [14].

2. THEORETICAL BACKGROUND 5

2.2 Object State

The goal of the orbit determination is to determine the state of an object at every time step. There are different ways to represent the object state. Space-based (space debris and observers) and ground objects (ground-based observers) state will describe in next section.

2.2.1 Space-based Objects

The object state is described in two ways, which are:

• Object Coordinates

• Classical Orbital Elements

To represent the object state using Cartesian Coordinates frame needs three positions ˆ r and three velocities ˆ v. For example, The International Space Station’s state vector is represented on the 2nd of September 2018 at midday [15]

X = ~ ~r

~ v



=



 

 

 



 

 

 

 r _I r _J r _K v I

v J

v K



 

 

 



 

 

 



=



 

 

 



 

 

 



5247.52821 Km

−2515.71764 Km

−3499.83983 Km 4.848498831 Km.s ⁻¹ 3.822398159 Km.s ⁻¹ 4.530854374 Km.s ⁻¹



 

 

 



 

 

 



(2.1)

The object state can also be described by the 6 classical Keplerian orbital elements which are described in table 2.1

Table 2.1: Definition of Classical Keplerian Orbital Elements [16].

Symbol Name Description

α Semi-major axis Describes the size of the orbit.

e Eccentricity Describes the shape of the orbit.

i Inclination Describes the orientation of the orbit with respect to the Earth’s equator.

Ω Right Ascension of the Ascending Node

Describes where the low point, perigee, of the orbit is with respect to the Earth’s surface.

w Argument of periapsis Describes the location of the ascending and descending orbit loca- tions with respect to the Earth’s equatorial plane.

v True Anomaly Describes where the satellite is within the orbit with respect to

perigee.

Figure 2.3: Keplerian Classical Orbital Elements [12].

The same state vector for ISS in Equation 2.1 is also represented by classical orbital elements as follow:

X = ~



 

 

 



 

 

 

 α e i Ω w v



 

 

 



 

 

 



=



 

 

 



 

 

 



6784.93557Km 0.0010405 51.66595 ^o 2.78131 ^o 104.38744 ^o 214.53843 ^o



 

 

 



 

 

 



(2.2)

2.2.2 Ground-based Observers

The geodetic coordinates latitude φ, longitude λ and altitude h are used for a ground-based observer.

WGS84 datum is used in SPOOK, which is a reference ellipsoid.

2.3 Observation Theory

In order to estimate the space object’s position, orbit determination algorithms use external measurements taken by observers. The observers can be ground-based observers from Earth or from a satellite which is orbiting around the Earth (space-based observers). Different observations methods have already existed for space objects, but the most common are passive optical and radar-based observers. The next sections will describe these two observers’ operational principles.

2. THEORETICAL BACKGROUND 7

2.3.1 Optical Observer

The working principle of an optical observer is based on the light reflection from the observed object. The reflected light can be translated into observer measurements which are right ascension α and declination δ angles. A space-based telescope was first used for space objects in 1996 [17].

Space-based sensors have advantages compared to ground-based ones such as not being affected by atmospheric conditions or geographic restrictions.

Figure 2.4: Angular measurements for optical observer: right ascension α and declination δ [12].

2.3.2 Radar Observer

Radar observers are ground base observers which actively send radio waves that are reflected from

the object and detected by the observers. They are two-way signals. Radar sensor can produce

angles (elevation el and azimuth β) plus distance measurements or slant range. In addition to the

angular angles and slant range ρ, Doppler measurements can measure the alteration of distance over

time or slant range-rate ˙ ρ.

Figure 2.5: Angular measurement of the radar observer [12].

2.4 Dynamic Model - Orbital Perturbation

The dynamic model plays an important role in calculating the object’s orbit. In order to calculate a successful orbit determination, a correct model of the space object must be included. All relevant forces which are acting upon the objects (e.g. space debris) must be taken into account in the orbit determination. This part of the thesis describes the forces which were already implemented in SPOOK before this thesis.

2.4.1 Non-spherical Gravitational Field

If the Earth’s mass concentration is assumed to be in a single point, the expectation will be for a spherical gravitational field. This case is possible only for a simple propagation, but in reality, the Earth’s mass distribution is not homogeneous, and the Earth shape is not spherical. Therefore, a non-spherical gravity field U for the earth is represented in Equation 2.3. More detail about the derivation of U can be found in [12].

U = µ r

"

1 +

∞

X

l=2 l

X

m=0

 R _⊕ r

 ^l

P l,m [sin(φ gC

)] {C l,m cos(mλ sat ) + S l,m sin(mλ sat )}

#

(2.3)

• µ - the gravitational parameter

• r - the spacecraft’s position

2. THEORETICAL BACKGROUND 9

• l, m - the degree and order of the gravity field

• R _⊕ - Earth’s radius

• P l,m - Legendre functions

• φ gC

and λ sat - the geocentric latitude and longitude of the spacecraft

• C l,m and S l,m are the geopotential coefficients that define the Earth’s mass distribution [18].

2.4.2 Atmospheric Drag

In Leo Earth Orbit (LEO), the strongest perturbation is the atmospheric drag. For precise orbit modeling, this force model is one of the big challenges because of some parameters such as the drag coefficient, object orientation (e.g. space debris) are difficult to determine. The MSIS-00 model used in SPOOK for the computation of the atmospheric drag. This model is developed by the US Naval Research Laboratory [17]. According to [12], the drag acceleration can be computed as follow:

~a Drag = − 1 2 C D

A D

m _obj ρv ² _rel ~ v rel

|~v _rel | (2.4)

where C _D is the drag coefficient, which describes the spacecraft surface and the interaction of the atmosphere, A D is the object’s cross-sectional area, m obj is the mass of the object. The atmospheric density is represented by ρ which is related to the altitude. ~ v rel is the object velocity which is between the atmosphere and the object. To calculate the relative velocity, the Earth’s rotational speed ( ~ w _⊕ ) must be considered so that this is called thermospheric winds.

~ v rel = ~ v − ~ w _⊕ × ~r − ~v t,winds (2.5) where ~ v and ~ r are the object’s respective velocity and position. ~ v t,winds is the thermospheric wind velocity which affects in the upper part of the atmosphere.

2.4.3 Solar Radiation Pressure

If an object is exposed to solar radiation, the absorption and reflection of the photons cause an additional acceleration on the object. This phenomenon is called Solar Radiation Pressure. This perturbation force has bigger effects in higher orbits. The acceleration can be calculated as follows:

~a SRP = −p _SRP C R

A R

m _obj

~ r 

r  (2.6)

where A _R is the object’s cross-sectional area, C _R is the object’s reflectivity coefficient, m _obj is the object’s mass, r  is the relative position of Sun to the object position, and p _SRP is the solar radiation pressure value as 4.57 · 10 ⁻⁷ N.m ⁻² .

2.4.4 Third-body Perturbations

As with the Earth’s gravity field effect on the object, other planets similarly affect the object in an

Earth orbit. There are two main perturbations in the Solar System which are Sun and Moon. To

where ~ S is the relative position between the body and the Earth, and µ _body is the gravitational constant of the body.

2.5 Initial Orbit Determination

To compute Orbit Determination (OD) successfully the Initial Orbit must be known. In some cases, initial orbit values can be defined by a user. However, in some cases there is no initial information about the orbit. For that reason, it is important to imply the initial orbit determination algorithm for Orbit Determination using a few measurements from observers. In 1881, Carl Friedrich Gauss used declination and right ascension to determine Dwarf Planet Cere’s orbit. This method is called the Gauss algorithm. This method is much simpler if the data consists of the time of the flight and two position vectors [12]. This algorithm is developed by Gibbs and Herrick- Gibbs. These methods have already been implemented in SPOOK.

2.5.1 Gibbs Algorithm

In 1889, Josiah Gibbs introduced an improvement of the Gauss algorithm. This method uses three position vectors to derive the velocity at the middle time incident. The derivation can be found in [12] and it is implemented by [17] in SPOOK [19].

2.5.2 Herrick-Gibbs Algorithm

The previous approach (Gibbs algorithm) fails when the observations are separated by small angles.

Distinct from the Gibbs algorithm, the Herrick-Gibbs algorithm uses a Taylor-series expansion to obtain an expression for the middle velocity vector [12]. The mathematical expression can be found in [17].

2.6 Orbit Determination

In SPOOK, five different orbit methods were already implemented, and an additional OD method is implemented in this project, giving six OD methods. All filters have a different way to estimate position and velocity, but the greatest difference between the filters are how the filters use the data.

OD methods applied in SPOOK are Weighted Least Squares (WLS), Sequential Batch Least Squares (SBLS), Extended Kalman Filter (EKF), Unscented Kalman Filter (UKF), Extended Square Root Information Filter (ESRIF). The WLS uses all the available measurements at each iteration, while filter-based algorithms use a single measurement at each update step. The SBLS filter is similar to the WLS filter, but this filter uses smaller batches of data at each step. The EKF, UKF, USKF, and ESRIF filters work in real time. For these filters, initial state and initial covariance should be determined before filtering. These algorithms will be described shortly in the next sections.

2. THEORETICAL BACKGROUND 11

2.6.1 The Weighted Least Squares Method

The WLS filter is a batch filter. The initial state vector ~ X ₀ is propagated when all measurements are obtained by an observer:

X ~ j = f  ~ X 0 , t j



(2.8) After calculating the predicted step, the predicted measurements are calculated:

~ y p,j = h( ~ X j ) (2.9)

In order to find the residuals, the real measurements (~ y r,j ) and the predicted measurements are compared so that the residual ~b j can be found:

~b _j = ~ y _r,j − ~y p,j (2.10)

To calculate the covariance matrix, the algorithm needs to find the partial derivative matrix and the weighted matrix. The covariance matrix is calculated as follows:

P o =

n

X A ^T _j W A j

! −1

(2.11) where W is a weighting matrix that provides all measurements, which are considered with their standard deviation, and where n is the total number of measurements from the observer:

W =







1

σ

0 0

0 . . . 0 0 0 _σ ¹

m





 (2.12)

where m is dependent on the observer type. This is for example when the observer is an optical observer which has right ascension and declination so that m is the number of measurements. A j is the partial derivative matrix of the measurements which are considered with the initial state. The observation matrix H j , and the error state transition matrix are shown as follows:

A j = ∂~ y j

∂ ~ X 0

= ∂~ y j

∂ ~ X j

∂ ~ X j

∂ ~ X 0

= H j φ j (2.13)

The error state transition matrix is calculated by solving the differential equation, which is defined as follows:

φ ˙ _j = F j φ (2.14)

where F j is the matrix of partial derivatives.

F = d X ~ ˙ d ~ X =

 d~ v d~ r

d~ v d~ v d~ a d~ r

d~ v d~ v



(2.15) Equation 2.14 is solved using two ordinary Differential Equation Solvers, which are Runge-Kutta- Fehlberg [20] and Shampine-Gordon [21] in SPOOK. Eventually, the update of the state vector can be calculated as follows:

n ! −1 n

τ = t

nm (2.18)

2.6.2 The Sequential Batched Least Squares Algorithm

There is an advantage to using the SBLS algorithm. If we assume that we have 1000 observation measurements, the state vector and covariance have already been computed by the filter, and we receive new measurements from the observer. We should do the process again for 30 of measurements but the filter has already calculated the 1000 measurements. If the filter does the same process again, it will lose the time needed to calculate. So that the SBLS does not need to calculate again it, just needs to use the new measurements for next iteration. The SBLS considers the two input matrices which are ( P n

(A ^T _j W A j ))

old and ( P n

(A ^T _j W b j ))

old , which come from the previous iteration of the SBLS. For the SBLS algorithm, the WLS’s equations 2.11 and 2.16 are arranged as follows:

P _o =

n

X A ^T _j W A _j +

n

X A ^T _j W A _j

!

old

! −1

(2.19)

δ ~ X 0 =

n

X A ^T _j W A j +

n

X A ^T _j W A j

!

old

! −1 n

X A ^T _j W b j +

n

X A ^T _j W b j

!

old

!

(2.20)

2.6.3 The Extended Kalman Filter

Previous filters are based on the processing spread over a certain time interval which might be minutes, hours, days, or a batch of data. These techniques are not sufficient to completely model the perturbation forces. For instance, the atmospheric drag changes due to solar flux and geomagnetic activity. These filters’ characteristic means that the estimate is considered with a particular epoch.

Extended Kalman filter solves this problem. The main disadvantage of this filter is that the filter needs an initial state and covariance value but this issue is solved by an initial orbit determination approach. The EKF has two steps, which are the prediction step and the update step.

Prediction Step:

In the prediction step, the filter propagates the state vector ˆ X k and the covariance ˆ P k to the next time incident, the predicted ¯ X k+1 and ¯ P k+1 .

X ¯ k+1 = ¯ X (t _k+1| t k ) = Z t

t

˙ˆX k dt + ˆ X k Predicted State (2.21)

F = ∂ ˙ˆX t

∂ ˆ X t

Φ(t ˙ _k+1 , t _k ) = F (t)Φ(t _k+1 , t _v k) (2.22) P ¯ k+1 = Φ ˆ P k Φ ^T + Q Predicted Error Covariance (2.23)

2. THEORETICAL BACKGROUND 13

Update Step:

In time update step, the residual is represented by e b k+1 . H k+1 is the observation partial derivative matrix, R is the measurement noise matrix. The measurements update the state vector and the covariance estimation at each measurement. The Observation Partial Derivative Matrix is introduced by partial derivative of the observations with respect to the state vector [17].

H _k+1 = ∂Obs

∂x (2.24)

e b _k+1 = z − H _k+1 X ¯ _k+1 Residual (2.25) K k+1 = ¯ P k+1 H _k+1 ^T H k+1 P ¯ k+1 H _k+1 ^T + R  −1

Kalman Gain (2.26)

R =







σ ² _i 0 0 0 . . . 0 0 0 σ ² _m





 (2.27)

As stated before, R depends on the observer type. For example, if the observer is an optical observer, it has two measurements which are right ascension and declination angles so that R will be a 2 × 2 matrix. The diagonal elements represent the observer measurement’s characteristic in R.

X ˆ k+1 = ¯ X k+1 + K k+1 e b k+1 State Estimate (2.28) P ˆ k+1 = ¯ P k+1 − K k+1 H k+1 P ¯ k+1 Error Covariance Estimate (2.29)

2.6.4 Unscented Kalman Filter

The UKF approaches the state and covariance prediction using the unscented transformation. This transformation utilizes a set of points and propagates these points through actual nonlinear function.

These points are selected based on their covariance, mean, and also higher order moments [22]. The covariance and mean can be recomputed from the propagated points, utilizing a more precise result check against normal function linearization. This technique reduces the computational complexity, compare to Monte Carlo based propagation methods, where several thousand points must be propagated at each step, and enhances the estimate accuracy compared to linearised propagation at the same time.

Filter Model

The UKF begin by augmenting state vector (x ^α ) to the L dimension vector, where L is the sum of dimensions in the main state vector, measurement noise and model noise vector. The covariance matrix is augmented similarly to a L ² matrix [22]. The covariance matrix P ^α and the state estimate vector x ^α are represented in this augmented form as follows:

x _k−1 ^α =



 x ^α _k−1

0 w

0 v



 (2.30)

columns are computed as follows:

χ ^α _0,k−1 = x _k−1 ^α , i = 0 (2.32)

χ ^α _i,k−1 = x ^α _k−1 + (α q LP _k−1 ^α )

i

i = 1...., L (2.33)

χ ^α _i,k−1 = x ^α _k−1 − (α q LP _k−1 ^α )

i−L

, i = L + 1...., 2L (2.34) where symbol i means the ith column of the covariance matrix (square root of the covariance matrix).

This matrix is computed by means of a lower triangular Cholesky decomposition (P = AA ^T ). The α value is in the interval 0 < α ≤ 1 that defines the sigma-point spread. This value is usually selected to be quite low to prevent non-local effects. The sigma points χ ^α _k−1 can be decomposed into the χ ^x _k−1 rows, which includes the state; the χ ^w _k−1 rows, which includes process noise, and the χ ^v _k−1 rows, which includes the measurement noise. Each set of sigma-points is appointed a weight so that these weights are derived by checking against the instant of the sigma-points with the Taylor series [22].

The weighted mean and covariance are represented as follows:

w ₀ ^(m) = 1 − 1

α ² , i = 0 (2.35)

w ^(c) ₀ = 4 − 1

α ² − α ² , i = 0 (2.36)

w ^(m) _i = w ^(c) _i = 1

2α ² L , i = 1...., 2L (2.37)

The filter firstly predicts the following step by propagating the sigma-points through the state vector and measurement model. After that it computes weighted mean and covariance matrices [22] as follows:

χ ^x _k,k−1 = f (χ ^x _k−1 , u _k , χ ^w _k,k−1 ) (2.38)

x ⁻ _k|k−1 =

2L

X

i=0

w _i ^(m) χ ^x _k|k−1 (2.39)

P _k|k−1 ⁻ = P 2L i=0 w ^(c) _i h

χ ^x _k|k−1 − ˆ x _k|k−1 i h

χ ^x _k|k−1 − ˆ x _k|k−1 i T

(2.40)

Z _k|k−1 = h(χ ^x _k|k−1 , χ ^v _k−1 ) (2.41)

z ˆ _k|k−1 =

2L

X

i=0

w _i ^(m) Z _i,k|k−1 ^x (2.42)

2. THEORETICAL BACKGROUND 15

After completing the prediction step, the filter is updated with the new measurements by calculating the state-measurement, cross-correlation and the measurement covariance matrices, which are then utilized to calculate the Kalman gain (K k ):

P _zz =

2L

X

i=0

w ^(c) _i h

Z _i|k−1 ^x − ˆ z _k|k−1 i h

Z _i|k−1 ^x − ˆ z _k|k−1 i ^T

(2.43)

P _xz =

2L

X

i=0

w ^(c) _i h

χ ^x _i|k−1 − x ⁻ _k|k−1 i h

Z _i|k−1 ^x − ˆ z _k|k−1 i ^T

(2.44)

K _k = P _xz P _zz ⁻¹ (2.45)

The updated the state vector and covariance are represented as follows:

ˆ

x _k|k = x ⁻ _k|k−1 + K k (z k − ˆ z _k|k−1 ) (2.46) P _k|k = P _k|k−1 ⁻ − K k P zz K k T

(2.47)

2.7 Orbital Regions

In this section, the orbit regions will be classified for the objects. This classification will be determined according to object’s orbital parameters. These categories of the orbital regions are based on the ESA Space Situational Awareness program [23].

Table 2.2: Orbital Region Classification [19].

Name Perigee [km] Apogee [km]

Min Max Min Max

Leo resident 0 - 2000 0 - 2000

Leo transient 0 - 2000 2000 - ∞

Low MEO resident 2000 - 16000 2000 - 16000

Low MEO transient 2000 - 16000 16000 - ∞

High MEO resident (GNSS) 16000 - 33786 16000 - 33786

High MEO transient 16000 - 33786 33786 - ∞

GEO resident (i ≤ 20 ^o ) 33786 - 37786 33786 - 33786

GEO resident (i > 20 ^o ) 33786 - 37786 33786 - 33786

GEO transient 33786 - 37786 37786 - ∞

HEO 37786 - ∞ 37786 - ∞

Table 2.3 shows that orbital region is augmented in LEO orbit. To better analyse the influence of

the atmospheric drag for the optimum orbit determination interval, it was decided to subdivide the

LEO region into eight different categories. The difference between the tables can be seen comparing

table 2.2 and 2.3

Leo resident(1000-1500 [km]) 1000 - 1500 1000 - 1500

Leo resident(1500-2000 [km]) 1500 - 2000 1500 - 2000

Leo transient (0-500 [km]) 0 - 500 500 - ∞

Leo transient(500-1000 [km]) 500 - 1000 1000 - ∞

Leo transient(1000-1500 [km]) 1000 - 1500 1500 - ∞

Leo transient(1500-2000 [km]) 1500 - 2000 2000 - ∞

Low MEO resident 2000 - 16000 2000 - 16000

Low MEO transient 2000 - 16000 16000 - ∞

High MEO resident (GNSS) 16000 - 33786 16000 - 33786

High MEO transient 16000 - 33786 33786 - ∞

GEO resident (i ≤ 20 ^o ) 33786 - 37786 33786 - 33786

GEO resident (i > 20 ^o ) 33786 - 37786 33786 - 33786

GEO transient 33786 - 37786 37786 - ∞

HEO 37786 - ∞ 37786 - ∞

2. THEORETICAL BACKGROUND 17

Chapter 3

SPOOK

SPOOK (SPace Object Observations and Kalman filtering) is a configurable and versatile tool developed by Airbus Defence and Space. This tool is developed in Fortran and able to simulate Earth orbiting objects. In this chapter, the main abilities developed in SPOOK before this thesis began and how it works will be briefly explained to the readers. More information is possible in [17], [18], [19], [24], [25], [26].

3.1 SPOOK Structure

A functional diagram of SPOOK can be found in figure 3.3. SPOOK’s software structure can be divided into three main categories. In this section, these categories will be described, explaining the capabilities of SPOOK.

• Simulation Layer: This layer contains the Sensor Simulator Module and the Population Generator Module. The simulation category aims to taking measurements of the objects and simulates a population of these objects in Earth orbit.

• Analysis Layer: This module evaluates the efficiency of the surveillance strategy and orbit determination. This module also estimates the state vector of the space objects. To estimate the state of the object, synthetic data or real measurement can be used.

• Interface Layer: This module provides the interfaces with the program, including a visual- ization tool.

Figure 3.3 shows these layers interact each other. SPOOK is aimed to be able to perform the

whole simulation of the SST system.

Figure 3.1: SPOOK functional diagram [11].

3.1.1 Simulation Layer

Population Generator

This layer contains the Population Generator Module. The Population Generator module specifies the properties of the space objects to be tracked by the observer. The population of the object can be defined with two forms; one of them is with TLE files, another form is that a user can define the object’s characteristics such as orbital elements, state vector.

This module is also responsible for classification of the space objects into the different orbit regions (LEO, MEO, GEO with resident and transient regions). Also, it is in charge if generating the reference trajectory for each space object and propagating their state with regard to the time-span [11].

Sensor Simulator

This module generates the synthetic measurements associated to the simulated object. This module enables the use of different observation strategies and sensors. The sensor parameters include their location, the total number of observers, observer types, etc. Before using the real data, most of the tool in the development process uses synthetic data. Therefore, a measurement generation module is necessary for SPOOK. It also contributes to understanding the space environment. More details can be found on the measurement generation process in [18].

3. SPOOK 19

The Sensor simulator can simulate seven different observer types. These measurement types are as follows:

• Ground-based radar: This observer can generate the following sensor elements: Range (ρ, km), Azimuth (β, deg), Elevation (el, deg), Range-rate ( ˙ ρ, km/s). The location of the ground-based radar is fixed using geodetic coordinates. The location is determined by a user.

• Space-based radar: This observer has the same measurement information as the ground-based radars, but its position is defined by a specific user-defined orbit.

• Ground-based optical measurements: This observer is a telescope and provide Right Ascen- sion (α, deg) and Declination (σ, deg) information.

• Space-based optical measurements: This observer provides Right Ascension (α, deg) and Declination (σ, deg) information. Its orbit is defined using the initial state vector or classical Keplerian orbital elements.

The observer types can be seen in table 3.1. Seven observer types are defined and its sensor elements can be seen.

Table 3.1: Observer Types.

Observer Measurements

Type Name α σ β el ρ ρ ˙

Type 1 Optical X X

Type 2 Radar X X X

Type 3 Doppler Radar X X X X

Type 4 Doppler Radar X X X

Type 5 Radar X

Type 6 Doppler Radar X

Type 7 Doppler Radar X X

The observation strategies are divided in three main categories:

• Accessibility Mode: The observers simulate the objects with unlimited field of view. The geometric constraint is described as the line of sight between the target object and an observer’s sensor Field of View (FoV).

• Tracking Mode: The observer tracks the space objects (the line of sight points straight to the target object).

• Surveillance Mode: This mode uses a pointing pattern which is sourced from a specific file.

Different surveillance strategies can be simulated.

3.1.2 Interface Layer

This interface is designed to indicate the result of the orbit determination. Also, this tool is able to

show the sensor location (both for the space-based and the ground-based observer), the reference

orbit of the space object and an estimated orbit of the object. Figure 3.2 shows an example result

of orbit determination.