DF

An Estimator for Satellite Orbits
A Statistical Model for a Ground Based Radar System

Master’s thesis in Physics and Astronomy

Enna Kukic & Carl Kylin

Department of Electrical Engineering
CHALMERS UNIVERSITY OF TECHNOLOGY
Göteborg, Sweden 2020





Master’s thesis

An Estimator for Satellite Orbits

A Statistical Model for a Ground Based Radar System

Enna Kukic
&

Carl Kylin

DF

Department of Electrical Engineering
Chalmers University of Technology

Göteborg, Sweden 2020



An Estimator for Satellite Orbits
A Statistical Model for a Ground Based Radar System
Enna Kukic & Carl Kylin

© Enna Kukic, 2020
© Carl Kylin, 2020

Supervisor: Anders Silander, Saab Group
Examiner: Thomas Eriksson, Department of Electrical Engineering

Master’s Thesis
Department of Electrical Engineering
Chalmers University of Technology
SE-412 96 Göteborg
Telephone +46 31 772 1000

Cover: An illustration of the ground radar searching for a satellite. Images from Adobe
Stock ©.

Typeset in LATEX
Printed by Chalmers Reproservice
Göteborg, Sweden 2020

iv



An Estimator for Satellite Orbits
A Statistical Model for a Ground Based Radar System
Enna Kukic & Carl Kylin
Department of Electrical Engineering
Chalmers University of Technology

Abstract
A model for the signal measured by a radar station from a satellite transiting above is
investigated and presented. A simplified model consisting of one satellite in a circular
equatorial orbit of known direction with a radar on the equator of a non-spinning earth is
investigated. The Maximum Likelihood Estimate (MLE) and Cramér-Rao Bound (CRB)
for the model are found. The orbital parameters, meaning the radius and start angle of
the circular satellite orbit, are found with an average error of ±7 m and ±4e-4◦ for an
incoming SNR of 10 dB and ±7 km and ±0.4◦ for an incoming SNR of -20 dB. Since this
variance is only slightly above the CRB, it is concluded that there is room for improvement
but adding complexity to the model should be of greater interest. What should be noted
is that the CRB presented in the report is an approximation of the actual CRB, since the
modulated signal of choice has a difficult time derivative. By choosing a more suitable
modulation it might be easier to study the relationship between the estimate and the
CRB.

Keywords: Radar, Satellite orbit, Signal processing, MLE, CRB, Tracking.

v





Acknowledgements
We would like to begin by thanking our supervisor Anders Silander at SAAB for all the
helpful insights, discussions and guidance he has provided throughout the entire project.
We also wish to thank our examiner Thomas Eriksson for all the valuable input and ideas
that have progressed our project further.

Finally, we want to thank family and friends for all the support through the years, and
SAAB for giving us this wonderful opportunity to learn and develop at a leading radar
surveillance provider.

Enna Kukic & Carl Kylin, Göteborg, June 2020

vii





Contents

List of Figures xi

List of Tables xiii

1 Introduction 1
1.1 Limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Ethical aspects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

2 Theory 3
2.1 Modelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

2.1.1 Describing an orbit . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.1.2 The rotational motion of the Earth . . . . . . . . . . . . . . . . . . 5

2.1.2.1 Transformation matrix, GCRS to ITRS . . . . . . . . . . 6
2.1.3 Rotating to a specific longitude and latitude . . . . . . . . . . . . . 8
2.1.4 Radar setup and signal generation . . . . . . . . . . . . . . . . . . . 8
2.1.5 Antenna geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

2.2 Signal processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2.1 Maximum Likelihood Estimation . . . . . . . . . . . . . . . . . . . 11
2.2.2 Cramér-Rao Bound . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.2.3 Whittaker–Shannon interpolation . . . . . . . . . . . . . . . . . . . 14

3 Methods 15
3.1 Signal generation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

3.1.1 Signal modulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.1.2 Generation of data . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

3.2 Estimation of orbital parameters . . . . . . . . . . . . . . . . . . . . . . . . 18
3.2.1 Parameter choice and integration time . . . . . . . . . . . . . . . . 18
3.2.2 MLE of orbital parameters . . . . . . . . . . . . . . . . . . . . . . . 19
3.2.3 Variance and Cramér-Rao Bound . . . . . . . . . . . . . . . . . . . 20

4 Results 23
4.1 Estimation of orbital parameters . . . . . . . . . . . . . . . . . . . . . . . . 23
4.2 Detection search . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

4.2.1 Near target search . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
4.3 Variance of estimator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

5 Discussion 29
5.1 Signal modulation and numerical CRB . . . . . . . . . . . . . . . . . . . . 29
5.2 Estimation of parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

ix



Contents

5.2.1 Estimation of angle . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
5.2.2 Estimation of orbital parameters . . . . . . . . . . . . . . . . . . . 29

5.3 Simplifications and Assumptions . . . . . . . . . . . . . . . . . . . . . . . . 30
5.4 Complexity of model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
5.5 Antenna geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
5.6 Variance and Cramér–Rao bound . . . . . . . . . . . . . . . . . . . . . . . 32

6 Conclusion 35

Bibliography 37

A Transformation matrices to relate GCRS to ITRS I
A.1 Transformation matrix, GCRS to ITRS . . . . . . . . . . . . . . . . . . . . I

x



List of Figures

2.1 A satellite orbiting the Earth will follow an elliptical path with the Earth
in one of the foci. An ellipse can be characterised by its semi-major axis a
and semi-minor axis b. The closest passage of the satellite is the periapsis
at rp. The satellite location is defined by rS. . . . . . . . . . . . . . . . . . 4

2.2 The orbit of a satellite will most probably not occur in the geocentric x2y2-
plane, visualised by the grey plane. Therefore, the yellow orbital plane must
be described. This can be done by three rotations Ω, i and ω such that
the orbit occurs counter clockwise around the z1-axis with the periapsis
occurring on the x1-axis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2.3 The rotational motion of the Earth. There is a spin around the rotational
axis tilted ≈ 23.5◦ from the ecliptical pole. The rotational axis also under-
goes precession around the ecliptic pole with a periodic perturbation called
nutation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2.4 The Geocentric Celestial Reference System GCRS with a pole at C0 and
origin at Σ0, and the Celestial Intermediate Reference System CIRS with
a Celestial Intermediate Pole CIP at P and Celestial Intermediate Origin
CIO at σ. Point P expressed in GCRS coordinates is: x′ = sin d cosE,
y′ = sin d sinE and z′ = cos d, where angle d represent the tilt of the
rotation axis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

2.5 To the left (a) a cross sectional view of the satellite passing above the radar
station. From the point of view of the radar station the satellite passes by
at an angle α and distance r. The coordinate system used in this figure
is the one derived in earlier sections, where the satellite orbit has been
calculated from the center of the Earth. To the right (b) the view from
above of the satellite passing by. Depending on the orbital parameters, the
satellite need not pass straight over the radar station. Therefore another
angle γ is needed to describe the transit. Finally, in relating (a) and (b)
ỹ = cos γ z + sin γ y. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2.6 An incoming plane wave reaching a row of antennas with spacing d. The in-
coming wave fronts will reach two adjacent antennas with a time difference
∆t = d

c
sinα, where c is the speed of light. . . . . . . . . . . . . . . . . . . 11

3.1 The autocorrelation peak for a AWPN-modulated signal and a sinus-modulated
signal. The parameters for these modulations were chosen according to
Table 3.1. What can be seen is that the AWPN-modulated signal has a
sharper peak and has lower side lobes for the wrong delay than the sinus-
modulated signal. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

xi



List of Figures

4.1 Correlation plot for five sequences with a sequence time of 1 ms each and
a total time of 10 s. This search was done for the entire possible search
area, and is referred to as a detection search. The correlation measurement,
MLE, was performed for four different lengths of Inline-geometry antennas;
1, 10, 100 and 1000. When increasing the number of antennas the intensity
can be seen to gather in α0-direction. The target value is marked with a
green circle. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

4.2 Zoomed in correlation plot for five sequences during a sequence time of 1
ms each and a total time of 10 s. This search was done for an area of
rS ± 15 km and α0 ± 0.7◦, and is referred to as a near target search. The
correlation measurement, MLE, was performed for three different lengths
of Inline-geometry antennas; 1, 100 and 1000. When increasing the number
of antennas the intensity can be seen to gather in α0-direction. The target
value is marked with a green circle, while the estimated value is marked as
the magenta star asterisk. . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

4.3 Variance of the estimates for rS and α0 as a function of the SNR of the in-
coming signal. The variance is plotted together with the CRB for both the
parameters given one antenna, represented by a blue and red line respec-
tively. The dotted black line represent the error for the resolution limit,
meaning the square step length for the grid search. This can be seen to
be much lower then both the CRB and variance. Note that the scales are
logarithmic. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

xii



List of Tables

3.1 Choice of parameters for the estimator, defined and briefly explained. . . . 18

4.1 Values for the root mean square error of the estimate, square root of the
CRB and the step length for the grid search are presented for the highest
and lowest value of SNRin. . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

xiii



List of Tables

xiv



1
Introduction

Today, over 2000 active satellites orbit the earth according to the UCS Satellite Database
[1]. With an increasing militarisation of space the need to monitor these satellites arises.
For a defence industry company like SAAB an exploration of this area is of interest in
order to present viable concepts. In order to monitor the satellites one can observe the
transit using a ground based radar system and estimate their orbital parameters.

Monitoring the orbit of any object presents new challenges but also significant simpli-
fications compared to monitoring for example aircraft. Whereas aircraft can control their
trajectories, objects in orbit follow well defined paths given by their orbital parameters. A
radar monitoring aircraft typically estimates the distance and angle to the aircraft as well
as it’s velocity. For a satellite these parameters change very rapidly during the transit
and are therefore difficult to estimate. On the other hand, the orbital parameters of a
satellite are constant and typically do not change. This project can be seen as an early
pilot study before building a facility which could actually accomplish this. The end goal
of this type of system is to track if a satellite changes it’s orbit in order to better monitor
some geographical area.

The aim of the project is to create an estimator with which one can identify the orbital
parameters of satellites using a radar system of a suitable geometry. The satellites of
interest are located in the lower orbital region called Low Earth Orbit, LEO, with an
altitude of 500 to 2000 km above the surface of the earth. The radar system will operate
in L-band.

1.1 Limitations

In the model describing the satellite’s orbit, the report will only consider satellite motion in
a known direction around the equator. The satellites orbit is modeled as circular without
any perturbations. It will orbit a non-rotating Earth with the radar station positioned
on the equator at zero longitude. The theory for a more general orbit is presented, it was
simply not executed.

There are three main perturbations of the satellite orbit: the Earth’s deviation from
spherical symmetry, air drag, and perturbations due to the Sun and the Moon. The
most significant perturbation for a LEO satellite is the spherical deviation caused by the
gravitational field, but no such effect was taken into account.

Furthermore, refraction and reflection of the signal in the atmosphere will be omitted
from the scope of the project. Also, reflection, refraction, dispersion and scintillation due
to the Ionosphere will be omitted.

1



1. Introduction

1.2 Ethical aspects
When considering the project many ethical and societal questions can arise. The aim of
the project was to create an estimator with which objects in orbit, mainly satellites, can
be monitored. There was no real product developed, solemnly what one can call a pilot
study was performed. The utilization of future products, benefiting from this report, is
tracing satellite motion to identify suspicious behavior. Although future products would
definitely be in the military realm, they are more of a surveillance nature. One could
develop the product further and use it as a target finder to eliminate satellites, but this
is not within the scope of the project and would need several changes and additions.
Furthermore, this type of ground based radar system for tracking satellites already exists,
but the technique is yet to be developed in Sweden. One example is the United States
Space Surveillance Network.

When it comes to environmental impacts only a few have to be taken into consideration.
The execution of the product will not impact the environment itself, but the emitted EM-
waves from a radar system such as this might. Since the ground based radar system would
emit high intensity radiation, more precisely microwave radiation, the area close by would
have to be secluded. Another possible effect on the environment could be that it would
cause issues for birds for example. Because of the signals emitted an electric field will be
induced causing a magnetic field as well as radiation. This could possibly impact birds’
and other flying creatures sense of direction as well as their health.

2



2
Theory

In this chapter, the theory behind the physical model and performed signal processing
will be presented. To start with, the orbit and the rotational motion of the Earth will
be introduced and subsequently described from the perspective of a ground based radar
system. After this, the signal reaching the radar will be described in terms of the orbital
parameters. Finally, the maximum likelihood estimate and Cramér-Rao Bound of the
orbital parameters is found.

2.1 Modelling
Before going into the estimation-part of the problem one needs to fully describe the
physical model of the system both regarding the orbit and the radar setup. In this section
the theory behind the model of choice will be described in detail and the generation of
signal data will be presented.

First off, the satellite will follow an orbit around a rotating Earth, and since the
satellite has it’s own plane of reference the system must be explained in the eyes of the
radar. Moreover, the radar setup and the generated signal must be defined and modelled.

2.1.1 Describing an orbit
A single satellite orbiting the Earth is a classic example of a two-body problem. The
solutions to these problems are widely known and are presented in for example the lecture
notes written by Klioner [2]. In summary, the two bodies will orbit their common center
of mass in the focal point of an ellipse. Because the Earth is several orders of magnitudes
heavier than a satellite the Earth is approximately stationary in that focal point. To
begin describing this orbit, consider an ellipse with semi-major axis a and semi-minor
axis b as in Figure 2.1. The Earth is then located in one of the foci. In the context
of orbits an important quantity is the closest passage, the periapsis rp. The convention
when parmeterising an orbital ellipse is to use the semi-major axis a and the eccentricity
e =

√
1− b2

a2 . Using this convention the coordinates of the satellite can be parameterised
as

rS,1 = a(1− e2)
1 + e cos ν (cos ν x̂1 + sin ν ŷ1), (2.1)

where ν is the angle between the location of the satellite rS and the periapsis rp. This
angular parameter is known as the true anomaly and is used to locate an object moving
along a Keplerian orbit [2].

When measuring any type of signal, the incoming signal will be measured as a function
of time. It is therefore important to establish how the true anomaly evolves in time. This

3



2. Theory

Earth

Satellite

Oa

b

rp

rS

y1

x1

Figure 2.1: A satellite orbiting the Earth will follow an elliptical path with the
Earth in one of the foci. An ellipse can be characterised by its semi-major axis a
and semi-minor axis b. The closest passage of the satellite is the periapsis at rp.
The satellite location is defined by rS .

can be done by solving the Kepler Equation [2]

k[a(1− e2)]−3/2t =
∫ dν

(1 + e cos ν)2 + const., (2.2)

where k =
√
GME, G is the gravitational constant and ME the mass of the Earth. This

integral can not be solved analytically. To simplify the calculations a new anomaly must
be introduced, the eccentric anomaly E. The true anomaly ν can be expressed in terms
of the eccentric anomaly E and the eccentricity e such that [2]

cos ν = cosE − e
1− e cosE ,

sin ν =
√

1− e2 sinE
1− e cosE .

By introducing the eccentric anomaly the Kepler Equation can be simplified to

E − e sinE = M,

where M = n(t − T ) is the mean anomaly, n = ka−3/2 is the mean velocity and T is an
integration constant representing M = 0 [2]. By solving Equation 2.2 numerically the
time dependence of the elliptical orbit can be fully described.

However, not all satellites orbit the earth in the same plane. It is thus important
to describe the orientation of the orbit. This can be done in infinitely many ways, the
traditional being three Euler-like angles illustrated in Figure 2.2 [2]. The orbital plane of
the satellite is shown in yellow, the orientation being such that the satellite goes around
the z1-axis with the periapsis occurring on the positive part of the x1-axis. To reach this
system from some reference coordinate system x2y2z2, a geocentric system, the system is
first rotated Ω around the z2-axis so that the orbit crosses the x2y2-plane at Ω and Ω +π.
After this the system is rotated an angle i about the instantaneous x′-axis in such a way
that the orbit goes in a counter-clockwise motion around the z1-axis. The only thing

4



2. Theory

remaining is to orient the x1-axis such that the periapsis occurs on the positive part of
the x1-axis. This is done by a final rotation ω around the z1-axis. Having decided how
to describe the orbital plane, the motion in that plane can be translated to a reference
system through three rotations as shown by Klioner [2]

rS,2 = Rz(−Ω)Rx(−i)Rz(−ω)rS,1 (2.3)

where Rj(θ) is a rotational matrix around axis j by angle θ. The three angles Ω, i, and
ω are called longitude of the ascending node, inclination, and argument of periapsis [2].

x2
y2

z2

x1

y1
z1

Ω
ω

i

Figure 2.2: The orbit of a satellite will most probably not occur in the geocentric
x2y2-plane, visualised by the grey plane. Therefore, the yellow orbital plane must be
described. This can be done by three rotations Ω, i and ω such that the orbit occurs
counter clockwise around the z1-axis with the periapsis occurring on the x1-axis.

In summary, an orbiting object can be fully described by six parameters. These can be
chosen in many different ways, the convention being the semi-major axis a, the eccentricity
e, the true anamoly ν, the longitude of the ascending node Ω, the inclination i, and the
argument of periapsis ω.

2.1.2 The rotational motion of the Earth
Now that the satellite orbit has been fully described and modelled one needs to address the
fact that the Earth is also rotating, making the Earth’s plane of reference an accelerated
plane of reference.

The Earth rotates both with respect to the Sun and with respect to its own rotational
axis tilted at 23.43928◦ [3]. During this study only the latter will be of interest, and will
henceforth be referred to as the Earth’s rotational motion. The Earth’s rotational motion
consists mainly of three components: precession and nutation, spin and polar motion.
The first two components are illustrated in Figure 2.3.

The motion of precession is equivalent to the motion of a gyroscope. Precession is
the rotation of the rotational axis with respect to the ecliptical pole and nutation is the
periodic perturbation of the precession. The nutation mainly emerges due to solar and
lunar effects [3]. The spin of the Earth is simply the spin of the rotational axis, while the
polar motion is the movement of the rotational axis relative its crust.

All the motions described above imply that an observational system situated on Earth
is in an accelerated frame of reference. When observing celestial or astronomical objects

5



2. Theory

Ecliptical pole
Earth rotation axis

Precession
Nutation

Spin

Figure 2.3: The rotational motion of the Earth. There is a spin around the
rotational axis tilted ≈ 23.5◦ from the ecliptical pole. The rotational axis also
undergoes precession around the ecliptic pole with a periodic perturbation called
nutation.

this can cause incorrect measurements of velocity, direction and position. Therefore, a
reference system which rotates with the rotation of the Earth is of interest and such an
reference system is the International Terrestrial Reference System, ITRS [4]. ITRS is a
geocentric system which takes spin, nutation, precession and polar motion into account.
The celestial object, if rotating around the Earth, moves in a different reference system
called the Geocentric Celestial Reference System, GCRS [5]. The GCRS is still a geo-
centric reference system, but it’s coordinate-axes do not rotate with respect to those of
the Barycentric Celestial Reference System, BCRS, a non rotating system centered at the
barycenter of the Solar System [5]. Thus, the GCRS must be transformed to the ITRS
and that can be done according to the technical notes published by Gérard and Luzum
[6]. The transformation matrices to relate the GCRS to the ITRS can be written as

[ITRS] = (Rprec+nutRspinRpolar)T [GCRS],

where Rprec+nut, Rspin and Rpolar are the transformation matrices connected to the preces-
sion and nutation, spin and polar motion respectively [6].

2.1.2.1 Transformation matrix, GCRS to ITRS

The first rotation goes from GCRS to the Celestial Intermediate Reference System, CIRS,
which is a geocentric coordinate system based upon current position and date taking
precession and nutation into account [5]. The two reference systems are visualised in
Figure 2.4. However, since the period for precession is around 26 000 years and the
operating time of a generic satellite is in the magnitude of ten years this rotation would
only rotate the axis in the order of 10−1◦ − 10−2◦ during the satellites life span making
the effect insignificant [7].

6



2. Theory

x
y

z

X

Y

Z

- GCRS
- CIRS

σ

P

Σ0

C0

E
d

Figure 2.4: The Geocentric Celestial Reference System GCRS with a pole at C0
and origin at Σ0, and the Celestial Intermediate Reference System CIRS with a
Celestial Intermediate Pole CIP at P and Celestial Intermediate Origin CIO at σ.
Point P expressed in GCRS coordinates is: x′ = sin d cosE, y′ = sin d sinE and
z′ = cos d, where angle d represent the tilt of the rotation axis.

The second transformation will bring CIRS to TIRS, the Terrestrial Intermediate Ref-
erence System, taking spin of the rotational axis into account. This effect has a period of
one day, making it significant for an observational system. The corresponding rotational
matrix is

Rspin = Rz(−ERA),

where ERA is the Earth Rotation Angle and spans between the CIO, the Celestial Inter-
mediate Origin, och the TIO, the Terrestrial Intermediate Origin, along the CIRS equator
[6]. The ERA can be calculated using the following equation presented by Gérard and
Luzum [6]

ERA = 2π(UT1 Julian day fraction
+0.7790572732640 + 0.00273781191135448Tu),

where Tu is the difference between the Julian date of UT1 and the Julian date 2451545,
corresponding to the day starting at 12:00 UTC on January 1, 2000. UT1 is a form of
Universal Time taking polar motion into account, which differs at most with 0.9s to the
Coordinated Universal Time known as UTC [8]. The transformation matrix is written in
accordance to the IAU 2006 resolutions, and the transformation follows the CIO based
procedure [6]. A procedure that realizes the CIRS using the Celestial Intermediate pole,
CIP, as its z-axis and the CIO as its x-axis, which can be seen in Figure 2.4.

The final rotation takes the system from TIRS to ITRS by also accounting for polar
motion. Values for the polar coordinates xp and yp are published daily by the International
Earth Rotation and Reference Systems Service, IERS, and are of magnitude 10−2 − 10−1

arc seconds [8]. Assuming that the satellite operates for ten years and the polar motion
constantly moves the CIRS pole away from the GCRS pole (see Figure 2.4), the motion
would only rotate the coordinate system in the order of 10−4◦ during the satellites life
span, making this effect insignificant as well.

The entire transformation matrix now consists of

[ITRS] ≈ Rz(−ERA)T[GCRS] ⇐⇒ rS,3 = Rz(ERA)rS,2 (2.4)

7



2. Theory

To read more about the exact transformation matrices please refer to Appendix A.

2.1.3 Rotating to a specific longitude and latitude
Having reached a coordinate system centered in the core of the Earth with the x-axis
pointing along longitude ϕ = 0 and geodetic latitude λ = 0, what remains is to rotate
to a general system with the x-axis pointing in the direction of some other longitude and
latitude, where the radar is located. The desired transformation to coordinate (ϕ, λ) can
be described by

rS,4 = Ry(Λ)Rz(ϕ)rS,3, (2.5)
where Λ is the geocentric latitude. With geodetic latitude one means the angle between
the normal and the equatorial plane, while geocentric latitude is the angle between a
radial line and the equatorial plane [9]. Using the definitions of geodetic and geocentric
latitude and an ellipsoid model of the Earth the two can be related as

tan Λ = b2
E

a2
E

tan λ,

where aE and bE are the semi-major and minor axis of the ellipsoid. When taking into
account the elliptic shape of the Earth, what also varies is the radius of the Earth. This
can be shown to be

R2
E = a4

E cos2 λ+ b4
E sin2 λ

a2
E cos2 λ+ b2

E sin2 λ
.

Appropriate values for aE and bE can be taken from the WGS84 model defined in [10],
from which

aE = 6 378 137 km
bE ≈ 6 356 752.3142 km.

In words, the final coordinate system is one where the x-axis points towards zenith, y-axis
points east, and z-axis points north, centered in the core of the earth.

2.1.4 Radar setup and signal generation
From the previous sections a description of the satellite orbit has been established. The
final step in specifying the problem is thus to describe what signal the radar station will
measure. Consider a signal in time described by

S(t) = m(t) ei2πfH t, (2.6)

where fH is the carrier frequency andm(t) is the modulated signal. Modulating the signal,
either the amplitude or the phase of the waveform will make it possible to correctly
determine the distance to an object. The carrier frequency fH is the center frequency
around which the radar operates, and is much higher than any frequency of the modulating
signal m(t) [11]. Using a L-band radar the carrier frequency would be ∼1 GHz, whereas
the modulating signal might have a bandwidth in the neighbourhood of 10-100 kHz due
to the property of the carrier frequency explained above.

The signal presented above is emitted from the radar station, travels to the satellite,
bounces of the satellite, and then travels back to the radar station. The signal is thus
delayed by a factor

τ(t) = 2r(t)
c

(2.7)

8



2. Theory

where r is the distance to the satellite and c is the speed of light. This setup is illustrated
in Figure 2.5. In Figure 2.5a the cross section of a satellite transit is shown. From the
view of the radar station the satellite passes by at an angle α and distance r, which can
be found from RS and θ using trigonometric formulae. For example

r(t) = |rS − rE| =
√
r2
S + r2

E − 2rSrE cos θ(t). (2.8)

In Figure 2.5b the same transit is shown, but from above instead. Here the angle γ
gives the final specification for the satellite position. The description of the orbit with
the parameters rS, θ, and γ is similar to that of a standard spherical coordinate system.
Finally, to relate the two figures, ỹ = cos γ z + sin γ y.

Earth center

rS

r

rE

Radar station

θ

α

x

ỹ

(a)

r sinα

Radar station

γz

y

(b)

Figure 2.5: To the left (a) a cross sectional view of the satellite passing above
the radar station. From the point of view of the radar station the satellite passes
by at an angle α and distance r. The coordinate system used in this figure is the
one derived in earlier sections, where the satellite orbit has been calculated from the
center of the Earth. To the right (b) the view from above of the satellite passing by.
Depending on the orbital parameters, the satellite need not pass straight over the
radar station. Therefore another angle γ is needed to describe the transit. Finally,
in relating (a) and (b) ỹ = cos γ z + sin γ y.

As discussed earlier, the signal that returns to the radar station is s(t− τ). What will
be measured is

s(t; τ(t)) = m(t− τ(t))e−i2πfHτ(t). (2.9)

To understand the different factors in the measured signal, they will be explained in de-
tail. As stated earlier the first factor m(t − τ(t)) comes from the modulation and the
second factor e−i2πfHτ(t) comes from the carrier frequency, but as one can notice the pure
t-dependence is gone in the second factor. When measuring a high frequency signal,
that signal can be down-converted into a more easily measurable signal range. A down
conversion can be explained by mixing the signal with a local oscillator [12]. Imagine a
high frequency signal sin(ωHt) which should be measured. This signal is mixed with a si-
nusoidal signal generated from a local oscillator producing a signal sin(ωHt) sin(ωLOt).
Using standard trigonometric formulae it is easy to show that sin(ωHt) sin(ωLOt) =

9



2. Theory

sin((ωH + ωLO)t) + sin((ωH − ωLO)t). By filtering out the high frequency the output
is only the difference in frequency. In this study the high frequency signal is modeled as
being mixed with a local oscillator with the same frequency, down-converting the signal
to base band and making fH the center frequency.

In this project a pure phase modulation was used, meaning that

s(t; τ(t)) = eiφ(t−τ(t))e−i2πfHτ(t). (2.10)

There are still two effects to take into consideration. Firstly, the intensity of the reflected
signal will decrease with an increasing distance proportional to I ∝ 1/τ(t)4. Since the
signal is modelled as having only a phase modulation the modulation will not impact the
intensity. Still the intensity of the incoming signal will be very small and will depend on
the distance to the object. It will therefore not be constant in time. This effect will be
neglected though, because for short integration times the amplitude will approximately
be constant. Secondly, the reflection of the satellite has been assumed to act like a point
source. This means that the satellite is assumed to reflect the signal of one single point,
and radiate that signal like a single point would. This is of course not the case; an
actual target will have a complex geometry which will generate an array of reflections
from different points in space. This is difficult to model without previous knowledge of
the geometry. Imagine that the signal reaches two points on the satellite which are not
the same distance from the radar. The reflection travelling towards the radar station
will consist of a sum s(t; τ) + s(t; τ + dτ) of these two reflections. If the two points are
relatively close to one another the reflection is approximately s(t; τ) + s(t; τ + dτ) ≈
s(t; τ)(1 + e−2πfHdτ ). A simple model for an object which is not too big is thus to add
an unknown phase β to the signal. The intensity will of course also be affected by the
number of points generating reflections. A simplified case where the amplitude is known
and set to one is considered. Since the amplitude is approximately constant for short
integration times this will mostly affect the signal to noise ratio and is thus an acceptable
simplification. Finally the simplified signal model can be written as

s(t; τ) = eiβeiφ(t−τ(t))e−i2πfHτ(t). (2.11)

2.1.5 Antenna geometry
Observing a signal from only one antenna will generally give a good estimate of the delay
of a signal. If one wishes to obtain more information from the incoming signal the next
step is typically to add more antennas. This will add to the amount of data available
for analysis, but also encode new information. Consider the situation shown in Figure
2.6; N antennas in a row with spacing d and an incoming plane wave from some angle
α measured from zenith. The wave fronts of the incoming signal will reach the antennas
with a time difference ∆t = d

c
sinα, where c is the speed of light. In the case of a signal

with a narrow bandwidth compared to its carrier frequency fH , this time difference can
be approximated as a difference in phase

∆φ = 2πfH
d

c
sinα.

This due to the fact that the emitted frequencies can be approximated to one frequency,
the carrier frequency. The situation described here is almost that of the radar station,
with one significant difference. The wave fronts that comes back from the satellite are not

10



2. Theory

planar, but spherical. However, the wave front is approximately planar in the far field,
which is a valid approximation for a distance much greater than the Fraunhofer distance
rf [13]

r � rf = 2(Nd)2

λ
. (2.12)

The spacing d can thus be chosen freely, provided it is not unreasonably small or so big
that the plane wave approximation does not hold. A common choice is half wavelength
spacing d = λ

2 = c
2fH

, which gives a phase difference

∆φ = π sinα. (2.13)

This is a good choice since all angles α = [−π/2, π/2] give a unique phase shift. For
a carrier frequency fH = 1 GHz and a distance to the satellite of r = 500 km, half
wavelength spacing would mean d = 15 cm and impose a limit on the number of antennas
from Equation 2.12, N �

√
2r
λ
≈ 1800 antennas. Neither of these are unreasonable.

d ...

α

Figure 2.6: An incoming plane wave reaching a row of antennas with spacing d.
The incoming wave fronts will reach two adjacent antennas with a time difference
∆t = d

c sinα, where c is the speed of light.

2.2 Signal processing
After describing the model the next step is to review the incoming signal using signal
processing. Below are the methods and estimations that are used in this study.

2.2.1 Maximum Likelihood Estimation
Having described the model in Equation 2.10, the next step is to estimate the orbital
parameters from an incoming signal. One way of doing this is to use a maximum likelihood
estimate. In an actual system there will be noise generated from everything from non-ideal
electronic components to the cosmic microwave radiation. In this study only internally
generated thermal noise will be considered, meaning that the noise can be generated as
uncorrelated normal noise with mean zero. This gives a signal defined by

x = s(θ) + n (2.14)

11



2. Theory

where n ∼ CN (0, σ2
nI) is a complex normal random vector with mean zero and standard

deviation σn, x = [x(t1), ..., x(tN)]T is the signal that is measured by the radar station,
s(θ) = [s(t1; τ(t1;θ)), ..., s(tN ; τ(tN ;θ))]T is the signal model given by the parameters
θ, and N is the length of the signal. In this study θ are the orbital parameters. The
maximum likelihood estimate of θ is given by the max of the likelihood function p(x|θ)

θ̂ML = arg max
θ
{p(x|θ)} ,

where the likelihood is given by [14]

p(θ|x) = |2πσ2
nI|−1/2 exp

{
−1

2[x− s(θ)]H 1
σ2
n

I−1[x− s(θ)]
}
,

where I is an N by N identity matrix and H denotes the hermitian conjugate. This
maximum is of course the same as the maximum of the logarithm of the likelihood, or
log-likelihood, which is given by [14]

L(θ|x) = −N2 ln 2π − 1
2 ln

∣∣∣σ2
nI
∣∣∣− 1

2[x− s(θ)]H 1
σ2
n

[x− s(θ)]. (2.15)

The two first terms are independent of θ, hence

θ̂ML = arg min
θ
{[x− s(θ)]H 1

σ2
n

I−1[x− s(θ)]}

= arg min
θ
{[x− s(θ)]H [x− s(θ)]}

= arg min
θ
{|x|2 + |s(θ)|2 − 2Re{xHs(θ)}}

= arg max
θ
{Re{2xHs(θ)− |s(θ)|2}}.

The model defined in Equation 2.10 is a complex exponential, with only varying phase.
Hence |s(t; τ(t))| = 1 for all t and τ(t). Thus the estimate can be simplified further, giving
us

θ̂ML = arg max
θ
{Re{xHs(θ)}}. (2.16)

If the signal is modeled as having an unknown phase β as in Equation 2.11, this phase
can also be estimated. In this case the MLE can be written as

θ̂orb,ML = arg max
θorb,β
{Re{eiβxHs(θorb)}},

where s(θorb) is the model without the unknown phase added to it and θorb are the orbital
parameters. This is maximised when eiβxHs(θorb) is completely real, which is true when

eiβ = (xHs(θorb))H
|xHs(θorb)| .

Applying this to the maximum likelihood estimate means that

θ̂orb,ML = arg max
θorb

{
Re

{
(xHs(θorb))H
|xHs(θorb)| x

Hs(θorb)
}}

= arg max
θorb
|xHs(θorb)|, (2.17)

in the case of an unknown phase.

12



2. Theory

2.2.2 Cramér-Rao Bound
Provided an estimator satisfying certain conditions, a lower bound on the variance of the
estimate can be found using the so called Cramér-Rao Bound, CRB. This is useful to
visualize how good an estimate is compared to how good it could be. The CRB for a
multi-variate estimator is [14]

CCRB(θ̃) = J−1(θ̃), (2.18)
where θ̃ are the correct parameters and J is the Fisher information matrix, FIM, with
elements

[J(θ̃)]ij = −E
[
∂2L
∂θi∂θj

(θ̃|x)
]

= E

[
∂L
∂θi

(θ̃|x) ∂L
∂θj

(θ̃|x)
]

(2.19)

and the lower bound is given by

var[θ̂i] ≥ [CCRB]ii.

It can be noted that the FIM is symmetric, [J(θ̃)]ij = [J(θ̃)]ji, hence only the upper or
lower half and the diagonal need be computed. Equation 2.15 gives the log-likelihood

L(θ|x) = −N2 ln 2π − 1
2 ln

∣∣∣σ2
nI
∣∣∣− 1

2[x− s(θ)]H 1
σ2
n

[x− s(θ)].

Since not all the terms in the above expression depend on the parameters θ, the FIM is

[J(θ̃)]ij = −E
[

∂2

∂θi∂θj

−1
2σ2

n

(
xHx+ sH(θ)s(θ)− xHs(θ)− sH(θ)x

) ∣∣∣∣∣
θ=θ̃

]
.

Furthermore, sH(θ)s(θ) = ∑N
n=1 s

H(tn;θ)s(tn;θ) = ∑N
n=1 1 = N is constant in all param-

eters. Using this the FIM simplifies further to

[J(θ̃)]ij = − 1
2σ2

n

E

[
∂2

∂θi∂θj

(
xHs(θ) + sH(θ)x

) ∣∣∣∣∣
θ=θ̃

]
.

The expectation value E[ ] is taken for the random vector x from Equation 2.14, for
which E[x] = s(θ̃). This is independent from the parameters θi since x does not depend
on θ, hence the order of the expectation value and the derivatives can be interchanged
resulting in

[J(θ̃)]ij = − 1
2σ2

n

∂2

∂θi∂θj

(
sH(θ̃)s(θ) + sH(θ)s(θ̃)

) ∣∣∣∣∣
θ=θ̃

= − 1
σ2
n

∂2

∂θi∂θj
Re

{
sH(θ̃)s(θ)

} ∣∣∣∣∣
θ=θ̃

.

(2.20)

The parameters that need to be estimated are the orbital parameters and the unknown
phase β. Let θ = [θorb, β], where θorb are only the orbital parameters. Now, s(θ) =
s(θorb)eiβ, meaning that ∂s

∂β
(θ) = is(θ) and ∂2s

∂β2 (θ) = −s(θ). Consider the column of the
FIM given by θj = β;

[J(θ̃)]θiβ = −i 1
2σ2

n

∂

∂θi

(
sH(θ̃)s(θ)− sH(θ)s(θ̃)

) ∣∣∣∣∣
θ=θ̃

= 1
σ2
n

∂

∂θi
Im

{
sH(θ̃)s(θ)

} ∣∣∣∣∣
θ=θ̃

.

(2.21)

13



2. Theory

If additionally θi = β, that element of the FIM will be

[J(θ̃)]ββ = 1
2σ2

n

(
sH(θ̃)s(θ̃) + sH(θ̃)s(θ̃)

)
= N

σ2
n

. (2.22)

Having shown these two cases, consider now the case when θi, θj ∈ θ \ β = θorb. The
CRB is evaluated at the correct parameters θ̃ and β is not one of the parameters being
differentiated, hence β = β̃. In this case s(θ) = s(θorb)eiβ which leads to sH(θ̃)s(θ) =
sH(θ̃orb)s(θorb)ei(β−β̃) = sH(θ̃orb)s(θorb). Using Equation 2.20 this means that

[J(θ̃)]θiθj
= − 1

σ2
n

∂2

∂θi∂θj
Re

{
sH(θ̃orb)s(θorb)

} ∣∣∣∣∣
θorb=θ̃orb

. (2.23)

Using the same information for the case when θi 6= β, θj = β simplifies Equation 2.21 to

[J(θ̃)]ij = 1
σ2
n

∂

∂θi
Im

{
sH(θ̃orb)s(θorb)

} ∣∣∣∣∣
θorb=θ̃orb

. (2.24)

2.2.3 Whittaker–Shannon interpolation
If the sampling of a continuous-time band limited signal is sufficiently dense, then you can
resample the signal via the Whittaker-Shannon interpolation formula to get even more
samples. The meaning of sufficiently dense is that the band limit must be less than half
the sampling rate

B <
fs
2 ,

where fs/2 is the Nyquist sampling rate and fs is the sampling frequency [15].
Given a signal x[n], sampled in discrete time n, the continuous function can be calcu-

lated to
x(t) =

∞∑
n=−∞

x[n] sinct− nT
T

,

using the Whittaker-Shannon interpolation formula, also called sinc interpolation [15]. T
is equal to 1/fs. This method is used in the project to interpolate bandlimited discrete
sampled signals.

14



3
Methods

Having described all the theory, it is time for the estimation of orbital parameters. To
simplify the estimation problem, a circular orbit around the equator with the radar station
on the equator at longitude ϕ = 0 on a non spinning earth is considered. Additionally,
the satellite is limited to spin in a known direction. Because of these simplifications the
rotational matrices from Equations 2.3, 2.4 and 2.5

rS,4 = Ry(Λ)Rz(Φ)Rz(ERA)Rz(−Ω)Rx(i)Rz(−ω)rS,1

become equal to identity and rS,4 = rS,1. The result is no rotation of coordinates. Fur-
thermore, the circular orbit means that the eccentricity e = 0 which simplifies Equation
2.1 to

rS,4 = rS,1 = a(cos νx̂1 + sin νŷ1) = rS(cos νx̂+ sin νŷ).

This is similar to the case in Figure 2.5 with γ = π
2 , and ν = θ. In the case of a circular

orbit with a spherical gravitational field, the time evolution of the anomaly can be found
from Equation 2.2:

ν̇ = θ̇ =
√
GME

r3
S

.

A real radar system will have some limits on what angles it can perceive incoming signals
from, for example caused by surrounding terrain and limits in the lobe from the emitting
and receiving antenna. To this end only satellites that start with an initial incoming
angle from zenith α0 ∈ [−π/4, π/4] rad are considered (see Figure 2.5a). Additionally,
only satellites in a LEO orbit are considered, giving a limit on rS ∈ rE + [500, 2000] km
(see Figure 2.5a). Note that the orbital parameters that are to be estimated are radius
rS and start angle α0. Start angle refers to the angle at which the satellite is positioned
from zenith at the start of the integration time. The radar system is assumed to have
a number of antennas equally spaced with distance d = λ

2 along the equator, a so called
Inline-geometry. Since the angle α is measured in the equatorial plane, this gives extra
information about the current angle of the satellite.

All the calculations, estimations and generation of data were done in MATLAB and
will be described below. Since this study is a pilot study no prior data was collected.

3.1 Signal generation
Before presenting the methods of estimation, the signal generation must be explained.
When generating the signal, no prior data was collected or available. For the signal
generation, first a signal modulation was established and then the generation of signal
data was performed given a set of orbital parameters in MATLAB.

15



3. Methods

3.1.1 Signal modulation
First of all a suitable signal modulation was of interest. To examine this the autocorrela-
tion was studied [14]

R(τ) =
tmax∑
t=tmin

s(t)s(t− τ(t))H.

The main requirement was a sharp autocorrelation peak with a well-defined main lobe and
low side lobes. After studying multiple modulations the most fitting modulations, due
to their autocorrelation, were a sinusoidal frequency modulation and an autocorrelation
optimised white phase noise modulation. The modulated signals can be written as in
Equation 2.10, with the following signal modulations

mAWPN = eiφ(t−τ(t))

msin = e2πiBsin sin[2πfPRF(t−τ(t))]·(t−τ(t)),

where mAWPN is the periodic autocorrelation optimised white phase noise modulation,
henceforth referred to as AWPN-modulation, and msin is the sinusoidal modulation. The
frequencies fPRF and Bsin are the pulse repetition frequency and the band width frequency
respectively. The values for these are presented in Table 3.1 in Section 3.2.1. The phase
φ is a periodic set of data obtained after optimizing white phase noise. The data set
originally had 2000 data points but was resampled using the Whittaker-Shannon interpo-
lation formula described in Section 2.2.3, increasing the resolution by a factor Fres = 100
presented in Table 3.1.

(a) (b)

Figure 3.1: The autocorrelation peak for a AWPN-modulated signal and a sinus-
modulated signal. The parameters for these modulations were chosen according to
Table 3.1. What can be seen is that the AWPN-modulated signal has a sharper
peak and has lower side lobes for the wrong delay than the sinus-modulated signal.

In Figure 3.1 the autocorrelation for the two modulations is presented. As can be
seen in the figure both of the modulations result in a well-defined main lobe at τ = 0,

16



3. Methods

but the AWPN-modulation has lower side lobes for the wrong delay than the sinusoidal
modulation, meaning it is better at rejecting incorrect correlation matches.

A sharp autocorrelation peak will likely lead to a sharp peak when determining the
parameters. This means incorrect parameter choices will be rejected. When considering a
detection problem a very sharp peak can be difficult to find via a grid search. Hence both
the modulations suit near target searches better than detection problems, due to their
sharp peaks. Since near target searches are more relevant for the purpose of tracking
known satellites this is no issue. Moving forward the AWPN-modulation will be used
since it has lower side lobes.

3.1.2 Generation of data
The orbital parameters are uniform randomly generated within the limits α0 ∈ [−π/4, π/4]
rad and rS ∈ rE + [500, 2000] km. With these parameters, a given integration time and
time step defined in detail in the next section, a corresponding incoming signal can be
generated as described in Section 2.1.4. The randomly generated orbital parameters give
a unique delay τ(t) which can be found from Equation 2.7 using the expression for the
distance from the radar to the satellite r(t) in Equation 2.8. For more details see Section
2.1.4. The incoming signal data is generated according to Equation 2.14 using the signal
model described in Equation 2.11.

The unknown phase β is not randomly generated due to the fact that its value does
not need to be estimated. In the MLE, the unknown phase is what leads to the absolute
value of the correlation measurement instead of the real part as described in Section 2.2.1.
Furthermore, the FIM in Section 2.2.2 does not have any explicit dependence on β, see
Equations 2.20-2.24.

When adding noise to the incoming signal x this is done in MATLAB according to

x = s(θ) + n, ni = 10−SNRin/10
√

2
(u+ iv) = σn√

2
(u+ iv), (3.1)

where u, v ∼ N (0, 1) are normally distributed numbers with standard deviation one, and
SNRin is the signal to noise ratio of the incoming signal and noise.

Having generated the signal for one antenna, consisting of N samples in time, the signal
for all the other antennas need to be generated as well. Suppose the radar consists of
Nantenna. The model for the incoming signal can be organized into a N ×Nantenna matrix:

S =


s(t1; τ(t1))e0·∆φ(t1) · · · s(t1; τ(t1))e(Nantenna−1)·∆φ(t1)

... . . . ...
s(tN ; τ(tN))e0·∆φ(tN ) · · · s(tN ; τ(tN))e(Nantenna−1)·∆φ(tN )

 .
Here s(t; τ(t)) is the signal at the first antenna and ∆φ(t) = ∆φ(α(t;θorb)) is the phase
difference from Equation 2.13 which depends on both time and the orbital parameters.
To add noise to this model a matrix of noise N is generated instead of a vector as in
Equation 3.1 and the signal x is a matrix X instead of a vector. This need not confuse;
the data could be organised into a vector, the matrix is used simply because it simplifies
later calculations.

Note that the random generation of parameters is done by performing a grid search in
the relevant area of interest.

17



3. Methods

3.2 Estimation of orbital parameters
In this section the estimation of orbital parameters will be presented. First the choice of
estimator parameters will be motivated followed by the calculations.

3.2.1 Parameter choice and integration time
In Table 3.1 all the used parameters are defined and briefly explained. To understand the
choices they will be motivated below.

Table 3.1: Choice of parameters for the estimator, defined and briefly explained.

Parameter Value Description
dt 0.1 µs Time step for all estimations
Tseq 1 ms Sequence time
Ttot 10 s Total time
Nseq 5 Number of sequences
Ngrid 200 Number of grid points in each direction
Tp 0.2 s Period of AWPN-modulated signal
fH 1 GHz Carrier frequency for both modulations
Bsin 100 kHz Bandwidth of sinus-modulation
fPRF 30 Hz Pulse repetition frequency for sinus-modulation
BAWPN 10 kHz Bandwidth of AWPN-modulation
Fres 100 Resolution factor for Whittaker–Shannon

interpolation
rlim [500,2000] km Maximum radius from the radar
αlim [-45,45]◦ Maximum angle from the radar

Firstly, the time step dt was chosen so that it could dissolve the Doppler shift properly.
Meaning

dt <
1

fH τ̇max
.

And, when calculating the correlation measurements the grid points were chosen to have
both a reasonable resolution and a reasonable computation time. However, when choosing
the period of the AWPN-modulated signal Tp the width of the correlation peak was taken
into account. For a shorter Tp, meaning a faster modulation, the bandwidth will be
increased and the correlation peak will be more defined in radial-direction making it
harder to detect and therefore demanding a better grid resolution. For the period to give
an unambiguous range of 2000 km back and forth, it has to be at least

Tp >
2rmax

c
≈ 0.01 s.

The period was therefore chosen to be greater than 0.01 s, but large enough to not demand
a higher grid resolution.

Secondly, to completely comprehend the orbit one needs to follow the satellite during
a longer time Ttot of a couple of seconds to collect enough data. This size of integration
time causes a high computation time though, and eliminates the possibility of assuming a
constant amplitude and incoming angle during the integration time. To go around these
issues the integration time will be divided into multiple sequences equally distributed.

18



3. Methods

During a shorter sequence time the delay will not change much, it behaves linearly, and
the incoming angle and amplitude will therefore be approximately constant. Instead of
looking at data from the entire integration time, Nseq sequences of length Tseq during a
total time Ttot are selected. See Table 3.1 for exact values.

Finally, the frequencies were chosen to fit the radar specifications and the sampling of
the AWPN-modulation.

3.2.2 MLE of orbital parameters
The approach for the MLE of the orbital parameters was to measure |xHs(θ)| for different
parameters surrounding the correct set of parameters; to perform a grid search. This was
done according to Section 3.1.2. By iteratively changing the parameters for the incoming
signal and matching it with the emitted signal described in Equation 2.6, according to the
maximum likelihood estimate described in Equation 2.17, most of the intensity gathered
near the correct answer resulting in a peak. This peak could then be identified by simply
finding the maximum of the grid search.

The approach was done both for detection of a satellite, where no prior orbital infor-
mation was assumed, and for a near target search where the approximate orbit was known
prior to the search. The grid search consisted of 200× 200 points for the radius and start
angle respectively as stated in Table 3.1. However, for both these searches the integra-
tion time was divided into shorter sequences as explained and motivated in Section 3.2.1.
This will effect the likelihood function and therefore also the MLE. If the data from the
two sequences are assumed to be independent the likelihood functions for the individual
events will multiply [16]. The two sequences are assumed to be generated from the same
orbital parameters but with different unknown phases β. The likelihood function for two
events x1 and x2 therefore result in

p(θ|x1 ∪ x2) = p(θ|x1)p(θ|x2).

Due to this the log-likelihood function will sum up for the sequences. The orbital param-
eters are identical for each sequence and the unknown phase β is not. This combined with
the logic that gave the MLE in Equation 2.17 means that the MLE is given by

θ̂ML = arg max
θorb,β1,...,βN

{Re{eiβ1xH1 s1(θorb) + ...+ eiβNxHNsN(θorb)}} (3.2)

= arg max
θorb
{|xH1 s1(θorb)|+ ...+ |xHNsN(θorb)|}, (3.3)

where xn, sn, and βn, are the data, model, and unknown phase for sequence n. The
measure that is maximised to find the MLE will be referred to as the correlation measure,
since it is essentially a measure of how correlated the model s(θorb) is with the data x for
a parameter choice θorb.

As mentioned in Section 3.1.2 the data is organised into an N ×Nantenna matrix, where
N is the number of samples in time from one antenna element. Consider the correlation
measure |xHn sn(θorb)| for sequence n; this is not well defined when both xn and sn are
matrices. The correct measure would instead be to take the absolute value of the trace of
this matrix. A simplification is made to reduce the time needed to compute the estimate;
the incoming angle α(t) is assumed to be constant in time during a sequence. Given a set
of orbital parameters where the correlation measure is to be evaluated, the incoming angle
is assumed to be equal to what it would be in the middle of the sequence α(t) ≈ α(tN/2).

19



3. Methods

The correlation measure with this approximation can be simplified to∣∣∣ [e0·∆φ(α(tN/2)), ..., e(Nantenna−1)·∆φ(α(tN/2))
]
xHn sn(θorb)

∣∣∣ . (3.4)

Here xn is the previously mentioned N ×Nantenna matrix whereas sn is the model of the
signal that reaches the first antenna, a column vector of length N . The factor[
e0·∆φ(α(tN/2)), ..., e(Nantenna−1)·∆φ(α(tN/2))

]
describes the difference in phase at the different

antenna elements and is a row vector with length Nantenna. The impact of this approxima-
tion has been investigated by comparing the results with and without the approximation.
For the chosen sequence time Tseq the impact of the approximation was not noticeable.
The advantage of this simplification can be explained by considering the number of op-
erations needed to calculate the correlation measure for some orbital parameters θorb.
Without the simplification two matrices would be generated and then the trace would be
calculated for this matrix product. A more computationally effective way is to look at the
data as a vector of length N ·Nantenna and multiply this with the corresponding model of
equal length, since this effectively calculates only the trace. The computational complex-
ity of this operation is O(N2N2

antenna) [17]. Using the simplification, a 1×Nantenna vector is
multiplied with a Nantenna×N matrix and a N ×1 vector. The computational complexity
of this operation is either O(NNantenna)+O(N2

antenna) or O(NNantenna)+O(N2) depending
on the order of operations [17]. If we assume that MATLAB chooses the computationally
most effective way of doing this, i.e. takes into account which of N and Nantenna is largest,
the complexity of the simplified case is ∼ O(NNantenna). The non-simplified case thus has
a computational complexity which is roughly squared compared to that of the simplified
case.

3.2.3 Variance and Cramér-Rao Bound
In Section 2.2.2 the concept of a lower bound for the variance of an estimate, the Cramér-
Rao Bound, was introduced. To evaluate how good an estimator is, it can be compared
to the CRB. If it is close, it is close to the limit of how good any estimator can be. It
is therefore important to calculate the variance of the estimates for different signal to
noise ratios, SNRin = −10 log10 σn [dB], and compare this to the CRB. In this study
the variance was calculated for a given SNRin 20 times in order to get a relatively good
estimate of the variance while not taking too long to perform. The SNRin ranged between
-20 dB to 10 dB. The variance could then be approximated according to

var[θ̂i] ≈
20∑
i=1

(θ̂i − θ̃i)2

20 ,

where θ̂i is the estimate and θ̃i is the correct value. The CRB is calculated as the inverse
of the FIM according to Equation 2.18 and thus the only remaining issue is to calculate
the FIM. According to Equation 2.19 the elements of the FIM can be calculated according
to

[J(θ̃)]ij = −E
[
∂2L
∂θi∂θj

(θ̃|x)
]

= E

[
∂L
∂θi

(θ̃|x) ∂L
∂θj

(θ̃|x)
]
.

The correlation function is measured using short sequences during a longer integration
time as mentioned before. Since the data produced during this time is considered to
be IID except for the unknown phase β, the log-likelihood will split into terms for each

20



3. Methods

sequence which will add up to the total log-likelihood as explained in the previous section.
Considering two sequences x1 and x2, the FIM can be written as

[J(θ̃)]ij = −E
[
∂2L
∂θi∂θj

(θ̃orb, β1|x1) + ∂2L
∂θi∂θj

(θ̃orb, β2|x2)
]
.

The two random phases β1 and β2 are not assumed to be the same. This means that
for every extra sequence considered an extra parameter βn is added to the problem. One
might ask what the FIM element for two different phases is. Let’s consider θi = β1, θj = β2;

[J(θ̃)]β1β2 = −E
[

∂2L
∂β1∂β2

(θ̃orb, β̃1|x1) + ∂2L
∂β1∂β2

(θ̃orb, β̃2|x2)
]

=
{
L(θ̃orb, β̃1|x1) is independent of β2

}
= −E

[
∂2L

∂β1∂β2
(θ̃orb, β̃2|x2)

]
= −E

[
∂2L

∂β2∂β1
(θ̃orb, β̃2|x2)

]
=
{
L(θ̃orb, β̃2|x2) is independent of β1

}
[J(θ̃)]β1β2 = 0.

The FIM entries for θj = βn and θi ∈ θorb can be calculated with Equation 2.24

[J(θ̃)]θiβn = 1
σ2
n

∂

∂θi
Im

{
sHn (θ̃orb)sn(θorb)

} ∣∣∣∣∣
θorb=θ̃orb

,

where sn is the model for the data from sequence n. The data from the the other sequences
are not included since L(θ̃orb, β̃m|xm) is independent of βn for m 6= n. The FIM entries
for θi = θj = βn can similarly be calculated with Equation 2.22

[J(θ̃)]βnβn = N

σ2
n

,

where N = Tseq

dt
+ 1 is the number of samples in one sequence. Finally, the FIM entries

for θi, θj ∈ θorb can be calculated from Equation 2.23

[J(θ̃)]ij = − 1
σ2
n

∑
n

∂2

∂θi∂θj
Re

{
sHn (θ̃orb)sn(θorb)

} ∣∣∣∣∣
θorb=θ̃orb

.

In the case of the two parameters rS and α0 the FIM can be written in matrix form as

J(θ̃) =



JrSrS
JrSα0 JrSβ1 · · · · · · JrSβN

Jα0Rr Jα0α0 Jα0β1 · · · · · · Jα0βN

Jβ1rS
Jβ1α0

N
σ2

n
0 · · · 0

... ... 0 . . . . . . ...

... ... ... . . . . . . 0
JβNrS

JβNα0 0 · · · 0 N
σ2

n


,

where Jij = [J(θ̃)]ij. Since the chosen AWPN-modulated signal can not be described an-
alytically, neither can the derivatives in the equations above. These are instead calculated
numerically using a finite difference. Additionally, because of the random nature of the

21



3. Methods

chosen modulation, calculating derivatives numerically at each sample is not a good idea
since they may be very big for individual points. It is therefore better to calculate these
derivatives after calculating the scalar products. This is done using MATLABs function
gradient which is used to calculate the numerical gradient of the real and imaginary
part of the correlation field sHn (θ̃orb)sn(θorb). The step lengths for derivatives in each
direction were adapted so that good estimates for the second order derivatives could be
made, meaning that they were small enough to generate low error, but not so small as to
hit the limit in resolution generated by the discrete sampling of the AWPN-signal. The
values were chosen to drS = 60 m, dα0 ≈ 4.000 · 10−3◦.

22



4
Results

For all the figures and numbers presented in this section, uniformly generated parameters
in the limits for α0 and rS have been studied. However, for the sake of comparison all
figures in this Chapter will have the same orbital values; rS ≈ 8300 km and α0 ≈ 36◦.
These values give as representative an illustration of the results as any other. The radar
setup has Inline-geometry antennas of various lengths. No noise was added to the received
signal until looking at the variance, to better visualise the results.

The results will consist of two main parts; the estimation of orbital parameters and
the variance of the estimator.

4.1 Estimation of orbital parameters
Two approaches were conducted when estimating the orbital parameters. The first was to
look at a detection problem, meaning there was no prior information assumed regarding
the orbit of the satellite when proceeding with the grid search. The entire possible area of
interest was searched. The second approach was to perform a near target search, meaning
that prior information of the orbit was assumed when doing a grid search. Thus, only a
small selected area around the parameters of interest was searched.

4.2 Detection search
For the detection search a correlation measurement was done for four different lengths
of Inline-geometry antennas; 1, 10, 100 and 1000. The intensity of the correlation can
be visualised in Figure 4.1. What can be seen from the four subfigures is that for a
higher number of antennas the intensity of the correlation measurement seems to gather
in angular direction. The difference especially stands out for the two later cases using
100 and 1000 antennas. In the subfigures the target value, meaning the correct value, is
marked with a green circle. Inside this circle one can see that where the intensity of the
correlation measurement peaks, the area becomes yellow. For the first two cases, using
one and ten antennas, the intensity does gather but not enough to produce a difference to
the width of the yellow area. For 100 antennas the yellow area can be seen to shrink, but
for 1000 antennas it shrinks to an unnoticeable size. The peak however can still be located
close to the target value. Even though there is no visual difference between one and ten
antennas, ten antennas will produce ten times the intensity which of courses increases the
signal to noise ratio.

What also should be mentioned is that the estimated value is as close as it can be,
given that one does not sample the grid search at the correct value. The estimated value
is within half the step length of the grid search from the correct value. Even for one

23



4. Results

antenna this is true. Using 200 × 200 grid points the step length becomes 0.45◦ for α0
and 7.5 km for rS.

(a) (b)

(c) (d)

Figure 4.1: Correlation plot for five sequences with a sequence time of 1 ms each
and a total time of 10 s. This search was done for the entire possible search area,
and is referred to as a detection search. The correlation measurement, MLE, was
performed for four different lengths of Inline-geometry antennas; 1, 10, 100 and
1000. When increasing the number of antennas the intensity can be seen to gather
in α0-direction. The target value is marked with a green circle.

However, in Figure 4.1 the area of the correlation peak in each of the subfigures is
quite small and hard to properly see. To better resolve and visualise the correlation peak
for the different cases one can zoom in around the target value. In the next section the
correlation measure is shown for the same case as in Figure 4.1, only zoomed in. This is
referred to as a near target search.

24



4. Results

4.2.1 Near target search

For the near target search a new correlation measurement using the same number of grid
points was done for three different lengths of Inline-geometry antennas; 1, 100 and 1000,
visualised in Figure 4.2. By using the same number of grid points on an area of rS ± 15
km and α0 ± 0.7◦, the step length becomes 0.007◦ for α0 and 150 m for rS.

(a) (b)

(c)

Figure 4.2: Zoomed in correlation plot for five sequences during a sequence time of
1 ms each and a total time of 10 s. This search was done for an area of rS±15 km and
α0 ± 0.7◦, and is referred to as a near target search. The correlation measurement,
MLE, was performed for three different lengths of Inline-geometry antennas; 1, 100
and 1000. When increasing the number of antennas the intensity can be seen to
gather in α0-direction. The target value is marked with a green circle, while the
estimated value is marked as the magenta star asterisk.

25



4. Results

Now, the difference between using one and 100 antennas appears to not be that sig-
nificant. However, when using 1000 antennas the intensity gathers quite a bit in angular
direction making the correlation peak more well-defined. The number of antennas still
do not change the estimated values for the orbital parameters, these are still limited by
the step length of the grid. In Figure 4.2 the estimated value is marked with a magenta
asterisk, while the target value is marked with a green circle. The yellow areas indicate
high intensity.

Like in the previous section, more antennas will generate more data and will always
increase the SNR. For 1000 antennas there is a significant amount of information about
the angle encoded in the signal. Increasing the number of antennas of course comes at
the cost of increased computational time.

4.3 Variance of estimator
Until now the result of the correlation measurements have only been shown for no noise.
To analyse how well the estimator performs for different noise levels a sweep over different
SNRin was done for the near target search. The variance was then calculated for one an-
tenna with SNRin ranging between -20 to 10 dB. It was then averaged over 20 estimations
for each SNR, as explained in Section 3.2.3. To compare the results of the estimator the
CRB was calculated. Both the variance and the CRB for the two orbital parameters can
be seen in Figure 4.3. A setup with only one antenna was chosen to be able to get a good
estimate of the variance within a reasonable time.

(a) (b)

Figure 4.3: Variance of the estimates for rS and α0 as a function of the SNR of
the incoming signal. The variance is plotted together with the CRB for both the
parameters given one antenna, represented by a blue and red line respectively. The
dotted black line represent the error for the resolution limit, meaning the square
step length for the grid search. This can be seen to be much lower then both the
CRB and variance. Note that the scales are logarithmic.

In Figure 4.3 the variance of the estimates for rS and α0 are shown as a function of the

26



4. Results

SNR of the incoming signal. To avoid being limited by the resolution of the grid search
the resolution was increased with a factor ∝ 1/σn for ranging SNRin, where σn is defined
in Section 3.1.2, Equation 3.1. A dashed line can be seen in the subfigures which gives
the error for the resolution limit. The variance and the CRB are represented by a blue
and red line respectively. Note that both the scales are logarithmic. The variance for
the estimate, the CRB and the error for the resolution limit at maximum and minimum
SNRin are presented in Table 4.1.

Table 4.1: Values for the root mean square error of the estimate, square root of
the CRB and the step length for the grid search are presented for the highest and
lowest value of SNRin.

SNRin [dB] RMSE Square root of CRB Step length
rS α0 rS α0 rS α0

SNRin = -20 7.07 km 0.366◦ 3.95 km 0.213◦ 1.5 km 0.07◦
SNRin = 10 6.76 m 3.82·10−4◦ 3.95 m 0.213 ·10−4◦ 1.50 m 7.00·10−5◦

The CRB is calculated according to the method explained in Section 3.2.3, and is
based on numerical derivatives which are only approximations of the actual derivative.
The CRB shown in Figure 4.3 is thus only an approximation of the actual CRB. Since
the variance of both the estimate for rS in Figure 4.3a and α0 in Figure 4.3b have the
same slope as the CRB they seem to perform well. There is a small offset to the CRB,
perhaps indicating a small bias.

To conclude this section the estimator performs well, as good as ±7 m and ±4 · 10−4◦

for good incoming SNR and ±7 km and ±0.4◦ for poorer incoming SNR.

27



4. Results

28



5
Discussion

In this chapter a discussion will be presented, covering the signal modulation, parameter
properties, variance and Cramér–Rao bound, future applications and improvements. The
chapter will mainly focus on what can be done in future studies.

5.1 Signal modulation and numerical CRB
The signal was modulated with the AWPN-modulation which has a bandwidth of 10 kHz.
This signal modulation seems to perform well for this problem. What could be improved
is the signal bandwidth, which might be made larger in order to incorporate more in-
formation. The one disadvantage of the AWPN-modulation is the fact that it requires
the CRB to be calculated numerically. The numerical derivatives for the CRB are quite
sensitive to the step length in each direction, meaning it is difficult to know what a good
step length would be. For future studies, one might look at using some other sort of
modulation which has a closed analytical form and for which the derivatives can be cal-
culated analytically. This could be used to see if the estimate is close to the CRB, while
an AWPN-like signal could be used to produce lower variance but without the possibility
of comparing to the CRB.

5.2 Estimation of parameters
This section will discuss the estimation of the two orbital parameters; start angle α0 and
radius rS. The discussion will touch upon both how well the estimator performs and what
could be improved.

5.2.1 Estimation of angle
In this very simplified model, one antenna is sufficient to determine the incoming angle.
This is unique to this simplified model and will not be the case for a more detailed
model. Still, this is very different from a regular radar station and highlights the difference
between this problem and a regular radar.

5.2.2 Estimation of orbital parameters
The estimation of the orbital parameters performs well, as well as the grid search allows
it for zero noise and not far off the CRB for more noise cases as can be seen in Figure
4.3. The estimation works both for detection and for near target search, but could of
course be improved. One thing that could be improved is to avoid being dependent on a
rigorous grid search in order to obtain a good estimate, and this could be done by taking

29



5. Discussion

advantage of the shape of a rough surface. For example, one could do either a gradient
search or fit a 2D-Gaussian surface to the peak, but there are some difficulties. The orbital
parameters have a large difference when it comes to their order of magnitude. rS is much
larger than α0 and that introduces many difficulties regarding analysis of the correlation
peak. When attempting to do a gradient search close by the peak or to fit a 2D-Gaussian
surface to the peak, this difference causes problems. A gradient search would be a great
option to avoid rigorous sampling of grid points, but since rS is so much larger than the
angle, the gradient search algorithm avoids steps in rS-direction. Also, a 2D-Gaussian
fit would serve well as a lowpass filter of the surface, damping the surrounding, but it
experiences the same difficulties. The fitted surface is sensitive to the exact choice of
standard deviation σ in rS-direction. However, both these methods are very interesting
to study further but a normalisation of the coordinates would have to be done.

Also, it would be of interest in future studies to investigate how to properly sample the
grid. There might be an optimised way of distributing the grid points, perhaps having
more in areas where the intensity is high. One idea is to randomly distribute the grid
points after an initial correlation has been calculated for the first sequence, and to have
greater probability of grid point in an area of high correlation intensity.

Except for improving the current estimation, one might also consider a different estima-
tion method. One estimate that could be of interest is a Maximum a posteriori estimate,
MAP, which uses a priori knowledge to find the next estimate. This would mean that
while tracking the satellite this estimate would benefit from prior estimates made.

5.3 Simplifications and Assumptions
This report has considered a simplified version of the actual problem. Apart from not
considering the complete set of parameters the simplifications are listed below;

• Uncorrelated noise: In Section 2.2.1 the noise which reaches an antenna is as-
sumed to be completely uncorrelated in time. When only looking at one antenna
this is a good approximation. When looking at several antennas the only noise that
would be completely uncorrelated is the thermal noise generated in the electrical
components. The noise generated from the outside, e.g. from the cosmic microwave
radiation, will be correlated between the different antennas. This needs to be taken
into account into the model.
Another external effect that will generate time correlated noise is the reflections
from the Ionosphere. The Ionosphere consists of three main layers: D, E and F
[18]. It is characterized by it’s ionization from the solar and cosmic radiation, a
factor which effects radio waves. The D layer at 70 to 90 km is only active at day
time, its free electrons recombine with oxygen ions during the night, but even then
the reflections it causes are small. The E layer at 90 to 160 km functions as the
D layer but does not loose all its free electrons during night, and therefore causes
some signal loss during the night as well. Still, the reflections are considered to be
negligible. The F layer though, at 160 to around 500 km, consists of the greatest
concentration of free electrons which almost do not reduce over the night, they only
redistribute. The F layer reflects radio waves with frequencies up to about 35 MHz,
hence it should not pose an issue. Still there might some noise caused from these
reflections.
The current model is still a reasonable first step, but more information needs to be
added to account for all of the phenomena.

30



5. Discussion

• Constant amplitude: In Section 2.1.4 the amplitude A
τ(t)2 of the signal is con-

sidered constant during a short sequence. In fact it is also considered constant
during all of the sequences. Neither of these is correct, although the first is prob-
ably an acceptable simplification. This needs to be investigated further by adding
the amplitude to the set of parameters that need to be estimated. Incorporating
the amplitude into the parameters also means that the FIM matrix which is used
to calculate the CRB will change.

• Independent sequences: In Section 3.2.2 the data from two sequences are as-
sumed to be independent. The signal which is generated from the reflection is
assumed to have some added uncorrelated noise. The viability of this assumption is
discussed above, with the conclusion that the noise should be independent in time
but not between the different antennas. Hence the data in two sequences should
generate independent data.

• Independent unknown phase: In Section 3.2.2 the unknown phase β is assumed
to be completely independent between the different sequences. The unknown phase
β depends on the geometry of the satellite, more specifically on the part facing
the radar station. A satellite will probably face the surface of the earth or at
least not rotate with respect to the surface of the earth. Imagine that the phase
is measured during an entire transit. It will likely change during the transit and
might be used to further identify a particular satellite or if that satellite is facing
a different direction than before. If the phase were to be assumed to be constant
during all of the sequences the maximum likelihood estimate in Equation 3.3 would
go from arg maxθorb

∑
n |xHn sn(θorb)| to arg maxθorb |

∑
n x

H
n sn(θorb)|. Having tried

both ways, the latter seems to displace the maximum away from the true maximum
when the true maximum is not sampled. This is probably because the e−2πfHτ term
of the model rotates the phase very fast as a function of the parameters. The same
offset in parameters θ might produce different offsets in phase for the correlation
measure of two sequences. Hence adding these and then calculating the absolute
value does not seem to generate the correct answer.

• Constant angle: In Section 3.2.2 the incoming angle is assumed to be constant
during a sequence. This is done to reduce the computational complexity as discussed
in Section 3.2.2. As discussed the validity of this approximation was checked and
worked well during the short sequences.

• Non-spinning earth: In the beginning of Chapter 3 the Earth is considered to
be non-spinning. This is obviously something which will need to be added to the
model. However, the rotation is highly deterministic. Adding this effect would likely
only generate more information about the orbit.

• Spherically symmetrical gravitational field: In the model described in the be-
ginning of Chapter 3 the gravitational field of the earth is assumed to be spherically
symmetrical. In reality the earth is a spheroid; it is slightly pressed together at the
poles. The gravitational field is therefore not spherically symmetrical. This will
cause the orbital parameters to drift over time and is the largest effect on the orbit
between 300 km and 30 000 km [2]. For satellites with a lower of higher orbit other
effect may be important to study as well, such as air drag and perturbations due to
the Sun and the Moon.

31



5. Discussion

5.4 Complexity of model
The model that is considered in this report is a simplified case of the actual orbit of a
satellite around the Earth. Except for the simplifications regrading the model covered in
the previous section, the model that is considered is one where most of the parameters
are set to zero. The next natural step would be to introduce an elliptical orbit with some
orientation given by the angle ω in Figure 2.2. After this the two other angles Ω and i
might be introduced and a radar system which is not located on the equator might be
considered. This however, takes the problem from estimating two parameters to having
to estimate four and then six, i.e. increasing the dimension of the problem. It might
be beneficial to further study this simplified case, removing the simplifications from the
previous section first. This since the problem would still only be two dimensional, making
it easier to study for example the correlation measures in Figures 4.1 and 4.2.

5.5 Antenna geometry
The antenna geometry considered is one where the antennas are standing in a line along
the equator. It is obvious from Figures 4.1 and 4.2 that more antennas give more infor-
mation regarding the orbit. In this simplified case this merely adds information about
the incoming angle. When increasing the complexity and considering an elliptical orbit
around the equator, one antenna will still be sufficient to estimate the parameters. When
increasing the complexity one step more and estimating the angles Ω and i, an array
of antennas will be needed to estimate the parameters. In this case an inline array of
antennas is probably not the optimal antenna geometry. One might investigate different
geometries and their effect on for example the CRB for different parameters given a set
number of antennas.

5.6 Variance and Cramér–Rao bound
As can be seen in Figure 4.3 the variance for the estimates of rS and α0 are close to the
CRB. As mentioned previously, the CRB shown can be thought of as an estimate of the
actual CRB. There might be an slightly different offset to the CRB, since the numerical
derivatives might not produce exactly the correct answers, but the effect is negligible.
Using a different modulation the numerical derivatives could be avoided, making the
calculations more precise and computationally effective. What can be seen without any
doubt is that the variance is decreasing with the same slope as the CRB. This indicates
that the estimator performs quite well.

The offset to the CRB, the bias, could come from the number of grid points or some
other unknown factor. When stating that the bias could depend from the number of grid
points, this is obviously not from the resolution limit in Figure 4.3. Instead, imagine
for example that with the current resolution the MLE is some distance away from the
correct answer. There might be a set of parameters between two grid points that better
explain the data and if the correlation was sampled in this point it would produce a
lower variance. The solution to this problem would be to increase resolution. However,
increasing resolution of course also means that the computational time would increase.
By experimenting with different resolutions, it is possible that a good trade-off could be
found.

32



5. Discussion

In Section 4.3 it is explained that the resolution factor of the grid search is increased
with increasing SNR in order not to hit the limit imposed by a finite resolution. Increasing
the resolution also meant decreasing the limits which the correlation is calculated within,
essentially looking in a smaller area. A reasonable question is thus if the decreasing limit
is perhaps ’pushing’ down the variance of the estimate. The variance has been calculated
without increasing the resolution. The only difference from the current result is that the
variance hits the resolution limit around SNRin = −2 dB, hence the decreasing limit is
probably not causing this behavior. The value for the variance was averaged over 20
iterations for each SNR. The number 20 was chosen by iteratively adding more iterations.
The variance seen in Figure 4.3 seems to follow a clear linear trend, hence the decision to
stop at 20 iterations was taken.

33



5. Discussion

34



6
Conclusion

The conclusion of this study is that the chosen method was suitable for the scope of the
project. The estimator performs well for both detection and near target search, and seems
to come close to the CRB. The average error of the parameters is ±7 m and ±4e-4◦ for
good incoming SNR of 10 dB and ±7 km and ±0.4◦ for poorer incoming SNR of -20 dB.
Only a simplified case is considered, hence future studies might first focus on removing
these simplifications. To be able to better study the effect that certain simplifications
might have the signal modulation should be changed to one where the CRB is not based
on numerical derivatives. After these simplifications have been removed, one should focus
on increasing the complexity of the orbital model.

35



6. Conclusion

36



Bibliography

[1] “UCS Satellite Database,” Apr. 1, 2020. Accessed on: June 1, 2020. [Online]
Available: https://www.ucsusa.org/resources/satellite-database.

[2] S. A. Klioner, “Basic Celestial Mechanics,” 2016. Accessed on: Feb. 16, 2020.
[Online] Available: https://arxiv.org/pdf/1609.00915.pdf.

[3] S.-H. Na, “Earth Rotation – Basic Theory and Features,” in Geodetic Sciences -
Observations, Modeling and Applications (Shuanggen Jin, ed.), ch. 9, Rijeka,
Croatia: InTech, 2013.

[4] “The International Terrestrial Reference System (ITRS),” 2013. Accessed on: Mar.
12, 2020. [Online] Available:
https://www.iers.org/IERS/EN/Science/ITRS/ITRS.html.

[5] I. Ridpath, A Dictionary of Astronomy. Oxford Paperback Reference, New York,
NY, USA: OUP Oxford, 2012.

[6] P. Gérard and B. Luzum, “IERS Conventions (2010),” Tech. Note No. 36, IERS
Convention Centre, Frankfurt, Germany, 2010. Accessed on: Feb. 26, 2020. [Online]
Available: http://www.iers.org/TN36/.

[7] D. P. Stern, “Precession,” Oct. 10, 2016. Accessed on: Mars 22, 2020. [Online]
Available: https://www-istp.gsfc.nasa.gov/stargaze/Sprecess.htm.

[8] “EOP of today,” 2013. Accessed on: Mars 19, 2020. [Online] Available:
https://www.iers.org/IERS/EN/DataProducts/tools/eop_of_today/eop_of_
today_tool.html.

[9] “Latitude,” 2020. Accessed on: June 15, 2020. [Online] Available:
https://en.wikipedia.org/wiki/Latitude.

[10] “United States Department of Defense, World Geodetic System 1984: Its Definition
and Relationships with Local Geodetic Systems,” Tech. Rep. TR8350.2, National
Imagery and Mapping Agency, St. Louis, MO, USA, 2000. Accessed on: May 20,
2020. [Online] Available:
https://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin.pdf.

[11] “Carrier wave,” 2020. Accessed on: June 15, 2020. [Online] Available:
https://en.wikipedia.org/wiki/Carrier_wave.

[12] A. M. Niknejad, “Lecture 15: Introduction to mixers.” EECS142, Lecture Notes,
University of California, Berkeley, USA, 2015. Accessed on: June 15, 2020. [Online]
Available: http:
//rfic.eecs.berkeley.edu/~niknejad/ee142_fa05lects/pdf/lect15.pdf.

[13] T. L. Singal, Wireless communications, pp. 79–80. Tata McGraw-Hill Education,
2010.

[14] U. Spagnolini, Statistical signal processing in engineering, pp. 64–65, 117–120,
143–146. John Wiley & Sons Ltd, 1 ed., 2017.

37

https://www.ucsusa.org/resources/satellite-database
https://arxiv.org/pdf/1609.00915.pdf
https://www.iers.org/IERS/EN/Science/ITRS/ITRS.html
http://www.iers.org/TN36/
https://www-istp.gsfc.nasa.gov/stargaze/Sprecess.htm
https://www.iers.org/IERS/EN/DataProducts/tools/eop_of_today/eop_of_today_tool.html
https://www.iers.org/IERS/EN/DataProducts/tools/eop_of_today/eop_of_today_tool.html
https://en.wikipedia.org/wiki/Latitude
https://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin.pdf
https://en.wikipedia.org/wiki/Carrier_wave
http://rfic.eecs.berkeley.edu/~niknejad/ee142_fa05lects/pdf/lect15.pdf
http://rfic.eecs.berkeley.edu/~niknejad/ee142_fa05lects/pdf/lect15.pdf


Bibliography

[15] “Whittaker–Shannon interpolation formula,” 2020. Accessed on: June 1, 2020.
[Online] Available: https://en.wikipedia.org/wiki/Whittaker%E2%80%
93Shannon_interpolation_formula.

[16] “Likelihood function,” 2020. Accessed on: June 20, 2020. [Online] Available: https:
//en.wikipedia.org/wiki/Likelihood_function#Products_of_likelihoods.

[17] “Computational complexity of mathematical operations,” 2020. Accessed on: June
26, 2020. [Online] Available: https://en.wikipedia.org/wiki/Computational_
complexity_of_mathematical_operations#Matrix_algebra.

[18] M. B. McElroy, “Ionosphere and magnetosphere,” Aug. 8, 2012. Accessed on: Apr.
11, 2020. [Online] Available:
https://www.britannica.com/science/ionosphere-and-magnetosphere.

[19] United States Naval Observatory. Nautical Almanac Office, Great Britain. Nautical
Almanac Office, The Astronomical Almanac. St. Louis, MO, USA: [Department of
Defense], Navy Department, Naval Observatory, Nautical Almanac Office, 2011.
Accessed on: May 20, 2020. [Online] Available:
https://books.google.se/books?id=dtKGqkDBCaIC.

38

https://en.wikipedia.org/wiki/Whittaker%E2%80%93Shannon_interpolation_formula
https://en.wikipedia.org/wiki/Whittaker%E2%80%93Shannon_interpolation_formula
https://en.wikipedia.org/wiki/Likelihood_function#Products_of_likelihoods
https://en.wikipedia.org/wiki/Likelihood_function#Products_of_likelihoods
https://en.wikipedia.org/wiki/Computational_complexity_of_mathematical_operations#Matrix_algebra
https://en.wikipedia.org/wiki/Computational_complexity_of_mathematical_operations#Matrix_algebra
https://www.britannica.com/science/ionosphere-and-magnetosphere
https://books.google.se/books?id=dtKGqkDBCaIC


A
Transformation matrices to relate

GCRS to ITRS

The transformation matrices to relate the GCRS to the ITRS can be written as

[ITRS] = (Rprec+nutRspinRpolar)T [GCRS], (A.1)

where Rprec+nut, Rspin and Rpolar are the transformation matrices connected to the
precession and nutation, spin and polar motion respectively [6]. The different
transformation matrices will transform the system into different reference systems, with
the first transition going from GCRS to the Celestial Intermediate Reference System,
CIRS, which is a geocentric coordinate system based upon current position and date
taking precession and nutation into account [5]. The second transformation will bring
CIRS to TIRS, the Terrestrial Intermediate Reference System, taking spin into account.
Finally, the last transition takes the system from TIRS to ITRS by also accounting for
polar motion.
Note that the transformation matrices are written in accordance to the IAU 2006
resolutions and that the transformation will follow the Celestial Intermediate Origin,
CIO, based procedure [6].

A.1 Transformation matrix, GCRS to ITRS
The first transformation needed is a transformation matrix taking the system from the
GCRS to the CIRS. This is done by taking the precession and nutation into account
into account and the transformation matrix can be written as he corresponding
transformation matrix can be written as

Rprec+nut =

 1− aX2 −aXY −X
−aXY 1− aY 2 −Y
X Y 1− a(X2 + Y 2)

Rz(s), (A.2)

where X = sin d cosE and Y = sin d sinE are the coordinates of the celestial
intermediate pole, CIP, the angles E and d are marked in Figure 2.4,
a = 1/(1 + cos d) = 1/(1 +

√
1−X2 − Y 2), rotational matrix Rz(s) and s is the CIO

locator [6][3]. The CIO locator positions the CIO σ on the CIRS equator. It is the angle
corresponding to the difference in right ascension of the equators node due to precession
and nutation [19].
The second transformation needed is a transformation matrix taking the system from
CIRS to TIRS, accounting for the spin of the rotational axis. This effect has a period of
one day, making it significant for an observational system. The corresponding rotational

I



A. Transformation matrices to relate GCRS to ITRS

matrix is
Rspin = Rz(−ERA), (A.3)

where ERA is the Earth Rotation Angle and spans between the CIO och the TIO, the
Terrestial Intermediate Origin along the CIRS equator [6]. ERA can be calculated using
the following equation

ERA =2π(UT1 Julian day fraction
+ 0.7790572732640 + 0.00273781191135448Tu),

(A.4)

where Tu is the difference between the Julian date of UT1 and the Julian date 2451545,
corresponding to the day starting at 12:00 UTC on January 1, 2000. UT1 is a form of
Universal Time taking polar motion into account, which differs at most with 0.9s to the
Coordinated Universal Time known as UTC [8].
Finally, all that is left is to account for the polar motion moving the reference system
from TIRS to ITRS with the corresponding rotation matrix

Rpolar = Rz(−s′)Ry(xp)Rx(yp), (A.5)

where xp and yp are the polar coordinates of the CIP, Ri are rotational matrices and s′
is the TIO locator which positions the Terrestrial Intermediate Origin, TIO, on the
TIRS equator [6]. It is the angle corresponding to the difference in right ascension of the
equators ITRS and TIRS node due to polar motion.

II


	List of Figures
	List of Tables
	Introduction
	Limitations
	Ethical aspects

	Theory
	Modelling
	Describing an orbit
	The rotational motion of the Earth
	Transformation matrix, GCRS to ITRS

	Rotating to a specific longitude and latitude
	Radar setup and signal generation
	Antenna geometry

	Signal processing
	Maximum Likelihood Estimation
	Cramér-Rao Bound
	Whittaker–Shannon interpolation


	Methods
	Signal generation
	Signal modulation
	Generation of data

	Estimation of orbital parameters
	Parameter choice and integration time
	MLE of orbital parameters
	Variance and Cramér-Rao Bound


	Results
	Estimation of orbital parameters
	Detection search
	Near target search

	Variance of estimator

	Discussion
	Signal modulation and numerical CRB
	Estimation of parameters
	Estimation of angle
	Estimation of orbital parameters

	Simplifications and Assumptions
	Complexity of model
	Antenna geometry
	Variance and Cramér–Rao bound

	Conclusion
	Bibliography
	Transformation matrices to relate GCRS to ITRS
	Transformation matrix, GCRS to ITRS