Master thesis, 30 hp
Master of Science programme in Engineering Physics, 300 hp
Spring term 2021
Estimation of satellite orbits using ground based
radar concept
Jonas Gabrielsson
A big thank you to my supervisor Carl Kylin for being enthusiastically helpful and providing invaluable guidance during this project. Thank you to all of ”Teknisk Fysik” at Ume˚a University for being such a great programme. Thank you to my family and friends not least the ones I have gotten to know during my studies. Finally a big thank you to SAAB for giving me this opportunity, it has been an amazing learning experience.
Master thesis
”Estimation of satellite orbits using ground based radar concept”
Jonas Gabrielsson
Abstract
Today an abundance of objects are circulating in earth captured orbit. Monitoring these objects is of national security interest. One way to map any object in orbit is with their Keplerian elements. A method for estimating the Keplerian elements of a satellite orbit simulating a ground based radar station has been investigated. A frequency modulated continuous wave radar (FMCW) with a central transmitter antenna and a grid of receivers was modeled in MATLAB. The maximum likelihood estimator (MLE) was obtained to estimate the parameters from the received signal. The method takes advantage of the relations between the Cartesian position and velocity and the Keplerian elements to confine the search space. For a signal to noise ratio (SNR) of 10dB, the satellite was followed during a time period of 0.1s where the positions were found within average error of range ±1.4m, azimuth ±2.0·10−6rad and elevation
±8.4 · 10−7 rad. Using a linear approximation of the velocity the Keplerian elements were found within average error of i : ±0.0050 rad, Ω : ±0.0050 rad, ω : ±0.0058 rad, a : ±2.60 · 105m, e: ±0.0021 and ν : ±0.24 rad. A method to obtain more accurate estimates is proposed.
3
2 Ethical considerations 5
3 Theory 5
3.1 Orbital Dynamics. . . 5
3.2 Radar geometry . . . 7
3.3 Signal . . . 7
3.4 Direction of arrival . . . 8
3.5 Signal matching. . . 8
3.6 Maximum likelihood estimation . . . 9
3.7 Orbital parameters . . . 10
4 Method 10 4.1 Assumptions . . . 10
4.2 Constraints . . . 11
4.3 Signal construction . . . 11
Generating orbital parameters•Dividing the information 4.4 Strategy . . . 11
4.5 Estimation . . . 12
5 Results 12 5.1 Position estimation . . . 13
5.2 Orbital parameter estimation . . . 15
6 Discussion 16 6.1 Estimation . . . 16
6.2 Signal modulation . . . 16
6.3 Simplifications . . . 16
7 Conclusions 17
References 18
1. Background
The atmosphere of the Earth is becoming more and more populated by satellites. Keeping track of these satellites is of importance from a national security perspective. For an aerospace and defence company like SAAB, there is understandably an interest in de- veloping a tracking system for satellite traffic. The idea is to use a ground based radar station to track the satellites. One way to map orbiting objects is via their so called ”Keplerian elements” also known as ”orbital parameters”. The aim of this thesis has been to estimate these orbital parameters from sim- ulated radar signals. This study is based on work done in a previous master thesis project by Kylin, Kukic [11] where the satellites where restricted to a specific orbit. Now the task is to extend this model allowing any orbit to be found.
2. Ethical considerations
Since the study aims to develop a method which would be of use within the military realm, the need to contemplate ethical aspects arises. The end prod- uct considered in this study would be a radar which would allow surveillance of satellites passing over.
The main goal is thus to keep track of potential threats to national intelligence, hence a defensive purpose. The proposed method could not be directly used for e.g shooting down or harming satellites, that would require a different type of study.
3. Theory
3.1 Orbital Dynamics
For an object in elliptical orbit around a much heav- ier object, one of the foci can be assumed to be cen- tered in the heavier body as in fig. 1. Parametrising such an orbit can be done using Keplerian orbital parameters which describes the orbit’s shape, size, orientation and the position of the object. Assuming a spherically symmetric gravitational field, the posi- tion vector and its time derivatives can be described in a so called perifocal frame as [1]
r = h2
µ (1 + e cos ν )(Pcosν + Qsinν)
Figure 1. A satellite in elliptical orbit around Earth. The closest passage of the satellite is the periapsis at rrrp. The satellite location is defined by rrrwhich is at an angle, true anomaly (ν), from rrrp.
Where h is the specific angular momentum and µ is the standard gravitational parameter. a and b are the semi-major and minor axis respectively, which describe the size of the ellipse and e =
q 1 −b2
a2
is the eccentricity which essentially describes how
”non-circular” the orbit is, where e = 0 represents a perfect circle. ν is the true anomaly which is an angle representing the objects position on the orbit and is the only time dependent parameter. The directional vectors P and Q are orthonormal vectors that lie in the orbital plane. Using the fact that the specific angular momentum is given by h = pµ a(1 − e2) we may obtain
r = a(1 − e2)
1 + ecosν(Pcosν + Qsinν) (1) taking the time derivative can be shown to be (a more complete derivation is available at [1] chap.
2.10)
˙r =
r µ
a(1 − e2)(−Psinν + Q(e + cosν)). (2) furthermore from the gravitational force and New- tons II law,
¨r =−µ
|r|3r. (3)
and its time derivative ...r =3µr · ˙r
|r|5 r − µ ˙r
|r|3. (4)
The normal to the orbital plane is given by the vec- tor W.
An Earth Centered Inertial frame (ECI) is a coor- dinate system with origin at the Earth’s center of mass. Consider a J2000 ECI [2] with basis vectors I, J, K where K is parallel to the Earth’s rotational axis, I to the mean equinox and J 90° East about the celestial equator. I, J constitute the equatorial plane as in fig. 2. The transformation between these two coordinate systems is given by the following rotations.
Figure 2. The figure shows an ECI frame relative to an orbital plane frame
P = (cosΩ cosω − sinΩ cosi sinω)I + (sinΩ cosω + cosΩ cosi sinω)J + sini sinωK
Q = (−cosΩ sinω − sinΩ cosi cosω)I+
(−sinΩ sinω + cosΩ cosi cosω)J+
sini cosωK
W = sinΩ sini I +
−sini cosΩ J + cosi K
The orbital plane is oriented relative to the reference I − J plane. The angles i, Ω, ω are called inclination,
longitude of the ascending node and argument of periapsis, these are visualized in fig. 3. Together with a, e and ν, these are what make up the orbital parameters.
Longitude of ascending node
Argument of periapsis True anomaly
Inclination Ascending node
Reference direction Celestial body
Plane of r
eference
Orbit
Ω ω
ν
i
☊
و
Figure 3. The figure shows angles that can
describe the orientation of an orbital plane relative to a reference ECI frame.
I am throughout this report going to use
ψ = [i Ω ω a e ν ] as a vector containing all orbital parameters.
3.2 Radar geometry
Since this is study aims to describe the orbit from a radar signal, we must model the position from the radar’s point of view. Let q be the vector going from the Earth’s core to the radar station and ρρρ be the vector going from the radar station to the satellite as seen in fig. 4. The Cartesian position vector can thus be described as r = q + ρρρ.
Figure 4. The figure shows vectors that describe a satellite positions relative both the Earth’s center and a radar station
Since the coordinate system is not co-rotating with the Earth, the vector q will be time dependant and is described as
q = q[cos(φ )sin(θ )I + sin(φ )sin(θ )J + cos(θ )K]
where θ and φ are angles as in standard spherical coordinates giving the radar stations position. Since the Earth is not simply rotating about a single axis both angles will have a time dependence [4].
The vector ρρρ can be described using an alternative coordinate system as seen in fig.5with basis vec- tors ˆθ , ˆφ , ˆq all being the spherical coordinates as described above with radial direction being in q direction. Hence we get
ρ
ρρ = ρ[cosαsinδ ˆφ + sinα sinδ ˆθ + cosδ ˆq] (5) where ρ is the distance to the satellite or ”range”, α and δ are angles known as ”azimuth” and ”eleva- tion”. The purpose of this new coordinate system is so that we may easily link these angles to the position vector and thus the orbital parameters as ρρρ = r(ψ) − q. From eq. 5we can obtain the angles
α = tan−1(ρθ
ρφ) (6)
Figure 5. The figure shows a way to describe a satellite positions from a radars point of view
δ = tan−1( ρq q
ρ2φ+ ρ2
θ
) (7)
3.3 Signal
The signal I have chosen to study is a frequency modulated continuous wave (FMCW) signal with triangular modulation, where the phase is the in- tegral of the instantaneous frequency. The signal time t is modelled as, for some total signal time T, t ∈ [−T2 ,T2]. Using a passband signal, it’s com- plex baseband representation [5] of the transmitted signal will be
st(t,tf) = Atei2π( fct±(∆ f t2f−∆ f tfTm)),
where At, ∆ f , fc, Tm,tf is the amplitude, frequency sweep, carrier wave frequency, modulation period time and fast time respectively. Fast time is the time between each modulation pulse as tf ∈ [0, Tm].
The plus term is when 0 < tf < T2m and the negative term when T2m < tf < Tm. The signal time t will start from−T2 and go linearly toT2. Fast time tf will
start at some time (given by the initial phase of the modulation) tf0between 0 and Tmand will linearly increase until a signal time of t =−T2 + Tm−tf0has passed, then tf will be set to 0. This is repeated for each ∆t = Tm.
The time delay τ for the signal to bounce off a satellite and come back to the receiver is given by the range ρ as
τ (t) = 2ρ(t) c .
The high velocity of a satellite demands that the delay has to be a function of time. The received signal will be
sr(t) = Arei2π( fc(t−τ)+ fDt± (∆ f t2f−∆ f tf(Tm+τ)), Where Ar is the amplitude of the received signal.
The Doppler frequency shift fD of the signal is given by the range rate
fD=2 fc˙ρ(t)
c .
Nyquist sampling theorem tells us that the sample rate of a signal needs to be twice the highest fre- quency component to not lose information [6]. To avoid having to sample the signal at twice the carrier wave frequency fcone can use a down converting process. Using the expression for the Delay and the Doppler shift, mixing the transmitted and received signals and applying a low pass filter gives [7]
s(t, ρ(t), ˙ρ(t)) = ArAtexp{
i4π
c [ fcρ + fc˙ρt + ˙ρt2±∆ f
Tmρtf]}. (8) 3.4 Direction of arrival
Consider several antennas in a two dimensional N× M grid with positions relative some corner an- tenna given by rantenna= nd ˆφ + md ˆθ where n = 0, 1, 2, ..., N − 1, m = 0, 1, 2, ..., M − 1 and d is the spacing between each adjacent antenna.
The returning signal wave will travel in the direction of −ρρρ. Assuming a plane wave, it will hit the grid
d ...
α
Figure 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.
of antennas with angles α and δ as in eq.5. A side view of this is seen in fig. 6. The time difference between the first antenna hit and the others is given by
∆t = nd
ccos(α)sin(δ ) + md
csin(α)sin(δ ).
This difference will show up as a phase difference in the signal received by each antenna as
∆Φnm= 2π fc
d
c(ncos(α)sin(δ )+msin(α)sin(δ )).
one can choose d = λc/2 which is half the carrier wavelength and which reduces the equation to
∆Φnm= π(ncos(α)sin(δ ) + msin(α)sin(δ )). (9) 3.5 Signal matching
The measured signal will include noise, from, e.g., thermal effects. Each antenna is assumed to have uncorrelated noise w which is a complex normal random vector with mean zero and standard devia- tion σw. Recall that the phase difference between each antenna is given by ∆Φnmas in eq. 9, since the satellite is moving the angles α and δ are functions of time. The received data is modeled as
xnm= s(t, ρ(t), ˙ρ(t))ei∆Φnm(α(t),δ (t))+ wnm(t)
Attempting to match this signal to extract informa- tion about ρ, ˙ρ, α and δ is hard since we would have to guess both the initial values and the time develop- ment of each parameter. This is also of little interest since the end goal is to obtain the orbital parame- ters. Since the motion of the satellite is completely deterministic one could instead guess the position of the satellite at the middle of the segment as in fig.
7.
Figure 7. The figure illustrates how a satellite orbit segment can be described from the middle position if the orbital parameters are known
Recall that ψ = [i Ω ω a e ν], one can use the Taylor expansion of the immediate position at the center of the signal and its derivatives of eq. 1- eq. 4,
dnρρρ(ψ) dtn = dn
dtn(r(ψ) − q) (10)
ρ0= |(r(ψ) − q)|
˙ρ = ρρρ · ˙ρρρ ρ0
¨ρ = (ρρρ · ˙ρρρ)2
ρ30 +| ˙ρρρ|2+ ρρρ · ¨ρρρ ρ0
...ρ =3(ρρρ · ˙ρρρ)3
ρ50 − 3(ρρρ · ˙ρρρ)(| ˙ρρρ|2+ ρρρ · ¨ρρρ)
ρ30 +
3( ˙ρρρ · ¨ρρρ) + (ρρρ ·...
ρρρ ) ρ0
ρ(t) ≈ ρ0+ ˙ρt + ¨ρt2 2 +
...ρ t3
6 (11)
˙ρ(t) ≈ ˙ρ + ¨ρt + ...ρ t2
2 . (12)
The range is expanded to fourth order so that the satellite motion can be accurately described for up to 10 seconds [9]. The same procedure is done for the angles α and δ using eq.6and7. This allows us to model the signal information as dependent directly on one set of orbital parameters as
xnm= s(t, ρ(ψ), ˙ρ(ψ))ei∆Φnm(ψ)+ wnm(t) (13) 3.6 Maximum likelihood estimation
To estimate parameters from data one can use the maximum likelihood estimation (MLE) method.
The Maximum Likelihood Estimation (MLE) is defined as
ψˆML= arg max
ψ
p(ψ|x)
here x is the received signal for one antenna as in eq. 13and p is the likelihood function and is given by [12]
p(ψ|x) = |2πσn2I|−1/2
× exp
−1
2[x − s(ψ)]† 1
σn2I−1[x − s(ψ)]
.
where s(ψ) = s(t, ρ(ψ), ˙ρ(ψ)) is the signal model as in eq.13. From previous work [11] it is showed that the MLE function can be simplified to
ψˆorb,ML= arg max
ψ
|
N,M
∑
n,m=0
x†nms(ψ)ei∆Φnm(ψ)| (14)
which means that the parameters in ψ that maxi- mizes eq.14is the estimation of the parameters.
3.7 Orbital parameters
The orbital parameters can be obtained if the Carte- sian position and velocity vector is known at any point [10]. Assume that for any vector x, the norm is x = ||x||.
The specific angular momentum vector is h =
q
µ a(1 − e2)W = r × ˙r (15) and the eccentricity vector
e = ˙r × h µ −r
r, (16)
the vector pointing towards the ascending node is n = (0, 0, 1) × h.
The true anomaly can then be obtained from
ν =
(arccos(e·rer), for r · ˙r ≥ 0 2π − arccos(e·rer), otherwise (17) the inclination
i= arccoshz
h, (18)
the right ascenscion of the ascending node
Ω =
(arccos(nnx), for ny≥ 0
2π − arccos(nnx), otherwise (19) the argument of periapsis
ω =
(arccos(n·ene), for ez≥ 0
2π − arccos(n·ene), otherwise (20) and finally the semi-major axis
a= 1
2 r− ˙r2
µ
. (21)
This equation can be rearranged as r˙2= µ(2
r−1
a) (22)
from the specific angular momentum in eq.15we may form [3]
r · ˙r = ±p
r2˙r2− h2
which is r · ˙r = ±
r r2µ (2
r−1
a) − µa(1 − e2) (23) Since W is perpendicular to both r and ˙r we get r · W = 0 and ˙r · h = 0 which can be used to form
˙r ·
rKsin(Ω)
−rKcos(Ω) rJcos(Ω) − rIsin(Ω))
= 0. (24) Eqs. 22, 23 and 24, defines a sphere with three planes intersecting as seen in fig. 8creating four possible velocities. Which gives four possible val- ues of ˙r obtained from a set of r, Ω, a and e.
Figure 8. The figure shows a system of equations which solution is given by a sphere with three planes intersecting giving four solutions
A more thorough derivation of these equations and clarifications regarding how they correspond to the sphere and the planes is available in Appendix: Clar- ification of equations.
4. Method
4.1 Assumptions
I have chosen to use the ECI coordinate system which does not co-rotate with the earth, this means that the radar station will have significant motion from the earths rotation. From eq. 10we see that this motion can easily be implemented as the time derivative of q. I have however chosen to simplify the problem by assuming a non-rotating Earth and thus setting q as constant in time. I have also as- sumed constant amplitude of the signal and that the reflection off the satellite will act like a point.
4.2 Constraints
The orbits considered will all be in lower earth orbit (LEO) thus the range ρ will be constrained within 250 km< ρ < 2000 km. I am also assuming that the range rate ˙ρ will be within −3kms < ˙ρ < 3kms , this is simply from observing what the values usually where during the study. To ensure that there is no ambiguity in our range measurement the pulse rep- etition frequency has to be fpr f < c
2·106m [13]. The lower earth orbits will likely not be highly eccentric hence I have considered eccentricities e ≤ 0.2. For the satellite to be within radar vision the angles are constrained as π2 > δ > 0 and 2π > α > 0 . 4.3 Signal construction
A signal s(t, ρ(t), ˙ρ(t)) can be constructed as eq. 8, the dependence of ρ and ˙ρ has an implicit depen- dence on the position vector, as seen in eq. 10, which has dependence of the orbital parameters r(ψ). Thus we can model a signal to be dependent of the orbital parameters directly.
4.3.1 Generating orbital parameters
Creating the recieved signal xnm as in eq. 13 can be done by randomly generating orbital parameters.
But these parameters have to correspond to a pass over the radar station. One can create a position r by using randomly generated values of ρ, α and δ within the limits described above. This ensures that the satellite is within the radars vision. Taking a random eccentricity e also within the limit, gives a limit as to what the semi major axis a can be as 1+er < a < 1−er from which a can be randomly drawn. Now simply generating a value of Ω ∈ [0, 2π] one can obtain four valid velocities using eqs.
22,23and24. From these one can choose randomly and the rest of the parameters are obtained from eqs.
17,18and20.
4.3.2 Dividing the information
The signal as described in eq. 13will be an N×M×L matrix, where L is the amount of time samples, con- taining the signal for each antenna element. To pro- cess this data as a whole would be computationally heavy, so an approximation is done. The informa- tion about the phase is divided into an N×M matrix where the angles are assumed to be constant during
the integration time so we get
Xnm=s(tL/2, ρ(tL/2), ˙ρ(tL/2)) (25) ei∆Φnm(α(tL/2),δ (tL/2))+ wnm(tL/2).
Any single antenna will receive the signal as x = s(t, ρ(t), ˙ρ(t)) + w(t). (26) 4.4 Strategy
All of the orbital elements can be directly estimated from the received signal as described above, there are however several issues. The search space is massive and a sweep of e.g. 100 different values of each parameter would mean processing 1012sig- nals which is not even a high enough resolution.
I investigated of some of the parameters could be estimated while keeping others constant but to no avail. Since the parameters seem to all be heavily cross correlated all need to be estimated together. A better approach is to obtain the Cartesian position from the range, azimuth and declination via eq.5.
If one assumes linear velocity between small time steps, one can construct both r and ˙r1≈r2−r1
∆t from these estimates as seen in fig.9.
Figure 9.The figure illustrates how short segments of a satellite trajectory can be used to obtain information about each middle position
This allows for an estimate of the orbital parameters as described in section3.7. However these estimates are not going to be great since they depend heavily
on ˙r which is obtained via the linear approximation.
Better estimates can be obtained by estimating the velocity directly from a longer segment. If one assumes the position is well known, an estimate of the velocity as a function of Ω, a and e is possible using22 - 24 yielding four possible velocities as described in fig. 8. The correct velocity is chosen as the one that gives the highest value in the MLE function. The signal can now be constructed as s = s(t, r, ˙r(Ω, a, e)) and we have a method which only requires estimation of three parameters. Obtaining the rest of the parameters can be done either from section3.7or from direct estimation using eq. 1- 4.
4.5 Estimation
The radial distance ρ and velocity ˙ρ can be esti- mated directly from any single antenna as described in eq. 26. The time period T during which these are estimated is assumed to be so short that the range velocity time dependence is completely lin- ear, using: ρ(t) ≈ ρ0+ t ˙ρ. One could obtain better estimates by also estimating higher order terms as well and using longer integration times, this is just a matter of computational cost.
Then a grid containing possible values of ρ and ˙ρ is constructed, the values which maximizes eq.14 is then used in the next iterations where the search will be centered around this value and the interval to be searched is smaller. This is done iteratively until it is expected that the error from noise is greater than error from resolution size. The angles α and δ are estimated from eq. 25where there is one point in time and the integration is done solely over the number of antennas N and M. The Cartesian posi- tion can then be obtained from ρρρ(ρ, α, δ ) in eq. 5 and then r = ρρρ + q.
This process is repeated for several segments, all separated by a small time difference ∆t, this allows for an estimate of the Cartesian velocity vector
˙ri≈ ri+1−ri
∆t . The values of ρ, α and δ are collected for every segment. The true value of these change almost linearly over short times thus fitting the es- timated values to a line allows for further error re- duction.
From eq. 19,21and16we may obtain a first esti-
mate of Ω, a and e using the previously obtained r and ˙r. This is used for setting up a grid around these values to find the ˙r that maximizes eq.14.
5. Results
During the study all orbital parameters are created from a randomly generated valid orbit on a position within the radar field of view. For the sake of proper comparison the same set of parameters will be used, see table1.
Table 1. Generated orbital parameters
Parameter Value
i 1.055
Ω 0.613
ω 0.059
a 9.660 ·106m
e 0.200
ν0 0.649
The signal parameter values that were used during the simulations is seen in table2.
Table 2. Signal parameters
Parameter Description Value
fc Carrier frequency 200 MHz
T Segment duration 10 ms
∆ f Frequency sweep 1 kHz
fpr f Pulse repition frequency 300 Hz
dt Time step 2 ·10−6 s
N Antenna grid width 200
M Antenna grid length 200
Variance is defined in this section as the mean square error σ2= ∑
n
(θ0−θn)2
N . Signal to noise ratio SNR is defined (in decibels) as SNR = −10log10(σw) where σw is the standard deviation of the noise, since the amplitude of the signal is assumed to be a constant 1. The satellite was followed for 10 seg- ments, each 10 ms long, with 10 ms time separation between each segment. This process was repeated for several SNR’s.
5.1 Position estimation
The angles α and δ where estimated separately from the range and range rate. In fig.10we can see a contour plot of the MLE field when estimating ρ and ˙ρ at a SNR of -10dB, where the red circle and the green star are the correct and estimated values respectively. The amount of grid points is 500 × 500 giving an initial grid resolution of 3500m and 12m/s. The max is found at some point (ρMLE,
˙ρMLE). Fig. 11shows the next iteration where the field is recalculated but this time around ρ ∈ ρMLE± 45km and ˙ρ ∈ ˙ρMLE± 150m/s giving a resolution of 90m and 0.6 m/s. This procedure is repeated 10 times resulting in a final resolution of 1.4 · 10−5m and 4.4 · 10−9m/s, since this is smaller than the carrier wave length this guarantees that no error is from resolution size. The same procedure is done for α and δ as seen in fig. 12 and 13 but with 100 × 100 grid points, resulting in a final resolution of 4.6 · 10−10and 3.7 · 10−10. This whole procedure is also repeated for each of the 10 segments each being separated by 10ms.
Figure 10. A contour-plot of a MLE field when estimating ρ and ˙ρ when SNR = -10dB
Figure 11. A zoomed in contour-plot of a MLE field when estimating ρ and ˙ρ when SNR = -10dB
Figure 12. A contour-plot of a MLE field when estimating α and δ when SNR = -10dB
Figure 13. A zoomed in contour-plot of a MLE field when estimating α and δ when SNR = -10dB All obtained values were gathered for each segment as ρi, αiand δi, they were then fitted to first order using MATLAB’s polyfit function. Figures14and 15shows the variance with both the first estimates and the fitted estimates. As expected the estimates are worse for smaller SNRs. The polyfit is able to reduce the error at certain SNRs. At SNR > - 20dB we have a variance of ρ of σ2< 108m2which means that we can expect an error of less than 104m at a range of 1.6 · 106m. At SNR > -20dB the angles have a variance of σ2< 10−5which means that we can expect an error of less than 10−2.5.
Figure 14. The variance of ρ both with and without being fitted to a first order polynomial, at different SNRs
Figure 15. The variance of α and δ both with and without being fitted to a first order polynomial, at different SNRs
5.2 Orbital parameter estimation
Using the fitted values I could create rifor each seg- ment and then ˙ri=ri+110ms−ri. The orbital parameters where then calculated as described in3.7. Since all parameters except ν are constant in time I took the average e.g. Ω = ∑iΩi/10 and ν = ν1. Figures16, 17and18shows the variance of all the angles, semi major axis and eccentricity. It seems from these results that a, e, ν and ω have very poor estimates for SNR < 10dB. The estimates for i and Ω seems quite good for SNR > -10dB.
Figure 16. The logarithm variance of the angles i, Ω, ω and ν at different SNRs
Figure 17. The logarithm variance semi major axis aat different SNRs
Figure 18. The logarithm variance eccentricity e at different SNRs
In section4.4I describe how the estimates could be improved by modelling the signal as s(t, r, ˙r(Ω, a, e)).
Taking an estimated position r and extending the signal time around this position from 10 ms to 100ms should introduce more information about the velocity in the signal. I set Ω as the value from the first estimation and calculated the MLE field as a function of a and e. Fig. 19shows the resulting field. Since no peak could be discerned I could not obtain a reasonable estimate of the parameters and thus no further estimation were done.
Figure 19. A contour-plot of a MLE field when estimating a and e when SNR = -10dB
6. Discussion
6.1 Estimation
The estimates of the position vector seems quite good for SNR > −20 dB, and it seems to translate into very good estimates for the two angles i and Ω. The rest of the parameters however have quite unreliable estimates, especially a. This would not be an issue if the complete method had worked out.
The ambiguity in the MLE field when estimating a and e makes it impossible to discern a peak. I attempted to do a gradient ascent search on the field but there are too many oscillations so the search always ends up stuck in a saddle point. It would however be of interest to further investigate other methods of estimation. Unfortunately I arrived at this particular solution at a rather late stage of the project and did not have much time to analyze the problems. It would be interesting to study longer integration times e.g T = 10s, the issue for me was that this would mean both that the constant angle approximation is invalid and also longer integration times will demand more samples in time. This is too computationally heavy for my PC. Studying the effect of other signal modulations and other radar setups, such as two radars separated by a large distance in a bi-static configuration, would also be of interest.
In the pre-study I considered a machine learning approach to the problem which might have worked quite well for the simulations but if the model is trained on simulated data and then doesn’t work on real data, we haven’t really learned anything.
6.2 Signal modulation
I’ve studied three different signal frequency modu- lations, the saw tooth (linear chirp), the sinusoidal and finally the triangular. The problem with the linear chirp is that the Doppler frequency cannot be separated from the signal and just adds error to the range measurement [14]. Since the Doppler shift can be quite significant for a satellite I discarded this modulation. The sinusoidal modulation seemed promising but I had issues finding a clear peak in the MLE field. Thus I chose to work with the triangular modulation.
6.3 Simplifications
I have made several simplifications during this study, these have been made due to both lack of time and resources such as computing power. As well as to make the model work as is. The simplifications are listed below with a description and discussion of the implications of each.
• Non rotating Earth
I assumed the Earth to be non rotating, in a real application the motion of the Earth would obviously have to be added. It should not be a problem however since the dynamics are well known.
• Point source
The reflection off the satellite is assumed to come from one point, it was suggested in the previous work done by Kylin, Kukic, that the geometry of the satellite be modeled as some random phase term that could be estimated separately.
• Constant amplitude
The signal is assumed to have constant ampli- tude, this would not be the case in a real life situation since the intensity would decrease with distance. This would have to be added to the model.
• Uncorrelated noise
The noise between all antennas are assumed to be from completely uncorrelated thermal effects. In reality there would be effects such as background cosmic radiation and atmo- spheric disturbances to consider.
• Constant angle
The constant angle approximation seems to be fine for short sequences.
• Round Earth
The effects from gravity is assumed to be completely spherically symmetric, when in fact the Earth is slightly oblate.
• Single satellite
I have not considered the effect several ob-
jects would have on a single signal. This would have to be investigated.
• Homogeneous beamforming
I have assumed that the transmitted signal is a homogeneous spherically symmetric sig- nal that can occupy the full 90◦from zenith space above the radar. In reality a more nar- row beamforming would most likely be used, from which new problems may arise with finding and tracking the satellite.
7. Conclusions
I laid out a foundation for how a received radar signal returning from a satellite depends on the Ke- pler elements, however extracting that information is complicated. Finding the Cartesian position is far easier and with knowledge of both position and velocity at one point in time, all elements can be obtained. A method that works well for finding the position of the satellite is used with a linear approxi- mation of the velocity. This leads to good estimates of the inclination and longitude of the ascending node i and Ω, although not great estimates for the rest of the parameters. I proposed a method that directly estimates the velocity without the linear approximation to obtain better estimates of these, which I believe shows promise. There are some simplification done that would need to be adjusted to obtain a more realistic model.
References
[1] H. Curtis, ”Orbital Mechanics for Engineering Students” pp. 76-78 Butterworth-Heinemann, first ed. (2005)
[2] Wikipedia.org, ”Earth-centered inertial”,
[Online] available: https://en.wikipedia.org/wiki/Earth-centered inertial [accessed]: 05/20 - 2021
[3] B. Hennessy, M. Rutten, S. Tingay, and R. Young,
”Orbit Determination Before Detect: Orbital Parameter Matched Filtering for Uncued Detection”
2020 IEEE International Radar Conference (RADAR), (Mars 2020)
[4] Wikipedia.org, ”Earths rotation”,
[Online] available: https://en.wikipedia.org/wiki/Earth%27s rotation [accessed]: 05/20 - 2021
[5] Mathuranathan, GaussianWaves, ”Complex Baseband Equivalent Models”,
[Online]: https://www.gaussianwaves.com/2017/10/complex-baseband-equivalent-models/
[accessed]: 05/24 - 2021
[6] Smyth, ”Nyquist Sampling Theorem”,
[Online]: http://musicweb.ucsd.edu/∼trsmyth/digitalAudio171/Nyquist Sampling Theorem.html [accessed]: 05/24 - 2021
[7] Parrish, ”An Overview of FMCW Systems in MATLAB”, (May 2010)
[8] Wikipedia.org, ”Continuous-wave radar”,
[Online] available: https://en.wikipedia.org/wiki/Continuous-wave radar [accessed]: 05/24 - 2021
[9] B. Hennessy, S. Tingay, et al. ”Improved Techniques for the Surveillance of the Near Earth Space Environment with the Murchison Widefield Array”,
2019 IEEE Radar Conference (RadarConf)(April 2019)
[10] Schwarz, ”Cartesian State Vectors − > Keplerian Orbit Elements”
[Online] available:https://downloads.rene-schwarz.com/download/M002-Cartesian State Vectors to Keplerian Orbit Elements.pdf
[accessed]: 05/06 - 2021
[11] Kylin, Kukic
”An Estimator for Satellite Orbits, A Statistical Model for a Ground Based Radar System”
Department of Electrical Engineering, CHALMERS UNIVERSITY OF TECHNOLOGY, G¨oteborg, Sweden 2020
[12] U. Spagnolini, ”Statistical signal processing in engineering”, pp. 64–65, 117–120,143–146.
John Wiley & Sons Ltd, 1 ed., 2017.
[13] Wolff, Radartutorial.eu ”Maximum Unambiguous Range”
[Online] available: https://www.radartutorial.eu/01.basics/Maximum%20Unambiguous%20Range.
en.html
[accessed]: 05/06 - 2021
[14] Wolff, Radartutorial.eu ”Frequency-Modulated Continuous-Wave Radar (FMCW Radar)”
[Online] available: https://www.radartutorial.eu/02.basics/Frequency%20Modulated%20Continuous%
20Wave%20Radar.en.html [accessed]: 05/06 - 2021
r˙2= µ(2 r−1
a)
setting µ(2r−1a) = R2we have r˙2= ˙r2I + ˙rJ2+ ˙rK2 = R2 which gives a sphere with radius R.
Two of the planes intersecting the sphere comes from eq. 23 r · ˙r = ±
r r2µ (2
r −1
a) − µa(1 − e2) To show how this is a plane, it is the same as
rI˙rI+ rJ˙rJ+ rK˙rK= ± r
r2µ (2 r −1
a) − µa(1 − e2) putting rI = A, rJ= B, rK = C and
q
r2µ (2r −1a) − µa(1 − e2) = D we get
A˙rI+ B ˙rJ+C ˙rK = ±D. (27)
Which is the equation for two planes one for +D and one for −D.
Furthermore to show how to obtain eq.24, starting from that the normal is orthogonal to the plane r · W = ˙r · W = 0
where
W = sin Ω sin i I − sin i cos Ω J + cos i K then
r · W = rIsin Ω sin i − rJsin i cos Ω + rKcos i since r · W = r ·sin iW = 0
r · W = rIsin Ω − rJcos Ω + rKcot i = 0 solving for cot i
cot i = −rIsin Ω + rJcos Ω rK
Plugging in the expression for cot i
˙r · W = ˙rIsin Ω − ˙rJcos Ω + ˙rK−rIsin Ω + rJcos Ω
rK = 0
multiplying by rK gives
˙r · W = ˙rIrKsin Ω − ˙rJrKcos Ω + ˙rK(−rIsin Ω + rJcos Ω) = 0. (28) This is also the equation for a plane. Which is the same as eq. 24
˙r ·
rKsin(Ω)
−rKcos(Ω) rJcos(Ω) − rIsin(Ω)
= 0
So does the planes intersect the sphere? The signed distance from the planes in eq.27to the center of the sphere is
ρ = ± q
r2µ (2r−1a) − µa(1 − e2) r
and for an intersection we have to have ρ2< R2 r2µ (2r −1a) − µa(1 − e2)
r2 < µ(2 r −1
a).
Multiplying by µ−1and r2gives 2r −r2
a − a(1 − e2) < 2r −r2 a Which is
0 < a(1 − e2)
which is true for all e < 1 and since this study doesn’t consider hyperbolic orbits this is fine and we can assume that there is an intersection for all orbits considered. The signed distance from eq. 28to the center of the sphere is 0, thus we clearly have an intersection.
