UNIVERSITY OF GLASGOW
Degrees of MEng, BEng, MSc and BSc in Engineering
NAVIGATION SYSTEMS
(ENG5062 / ENG4184)
Thursday 17th December 2020
Release time: 09:30AM (BST) for 4 hours
This exam should take you: 2 hours to complete
However, you have a 4 hours window to download/complete/upload your submission
Attempt ALL THREE Questions.
The numbers in square brackets in the right-hand margin indicate the marks allotted to the
part of the question against which the mark is shown.
These marks are for guidance only.
A FORMULA SHEET IS PROVIDED AT THE END OF PAPER
A calculator may be used. Show intermediate steps in calculations.
Page 1 of 17 Continued overleaf...
Q1.
(a) The following coordinates were recorded using a hand-held GNSS at the beginning
and end of Sauchiehall Street in Glasgow. Calculate the length of the road.
Buchanan Galleries - 55◦51′51′′N 4◦15′11′′W
Charing Cross - 55◦51′58′′N 4◦16′14′′W
(8)
(b) A radar drone surveillance system is used to detect and track small airborne vehicles
operating close to restricted sites such as airports, military bases and nuclear power
plants. The radar antenna is mounted on a two-axis gimbal that allows it to rotate in
both azimuth and elevation. Each drone within range will reflect a small portion of
the radar energy from favourably aligned reflector points on the fuselage.
Using the navigation kinematic notation, show that the velocity of a reflector point p
on the fuselage as seen from the antenna and resolved into antenna axes is given by,
r˙aap = v
a
eb − vaea +CabΩbeblbbp −Ωaea
(
raeb − raea +Cab lpbp
)
where, F b,F a are free to rotate with respect to F e and lbbp is the position of the
fuselage reflector with respect to the body axes. (12)
oi,e xe
ye
ze
reea
reeb
reab = r
e
eb − reea
xa
ya
za
xb yb
zb
p
lbbp
Figure Q1: Coordinate representation of UAV and radar antenna.
[20 Marks]
Page 2 of 17 Continued overleaf...
Solution:
To solve this problem we need to obtain the location of each end in Cartesian coordi-
nates using,
xeeb = (RE(Lb) + hb) cosLb cosλb
yeeb = (RE(Lb) + hb) cosLb sinλb
zeeb =
[
(1− e2)RE(Lb) + hb
)
sinLb
There are a few variables we need to find first. The geodetic latitudes and longitudes
need to be converted to their decimal equivalents from deg/min/sec format. Using the
standard conversion,
Longitude (deg) Latitude (deg)
Buchanan Galleries -4.2531 55.8642
Charing Cross -4.2707 55.8663 (1)
(Note 1
2
point for each correct coordinate).To ensure correct calculation, the meridian
radius of curvature need to be calculated,
RE(L) =
R0√
1− e2 sin2 L
Using the latitudes for each end we get,
• Buchanan Galleries
RE(L) =
6378137√
1− 0.08182 sin2(55.8642)
= 6.392814× 106m
• Charing Cross
RE(L) =
6378137√
1− 0.08182 sin2(55.8663)
= 6.392814× 106m (2)
(Note 1 point for each correct radius of curvature). We are now able to calculate the
Cartesian position for each end of the road, assuming hb = 0m.
Page 3 of 17 Continued overleaf...
• Buchanan Galleries
xeeb = (6392814) cos(55.8642) cos(−4.2531) = 3577488m
yeeb = (6392814) cos(55.8642) sin(−4.2531) = −266047m
zeeb =
[
(1− 0.0818192)6392814] sin(55.8642) = 5255972m
• Charing Cross
xeeb = (6392814) cos(55.8663) cos(−4.2707) = 3577213m
yeeb = (6392814) cos(55.8663) sin(−4.2707) = −267132m
zeeb =
[
(1− 0.0818192)6392814] sin(55.8663) = 5256103m
(3)
(Note 1
2
point for each correct coordinate).
As this is a relatively short distance, there is no need to consider curvature of the earth
and so to calculate distance, use the simple Euclidian distance.
dist(a, b) = ∥b− a∥ (1)
The vector from ’a’ to ’b’ is then,
b− a =

3577213.4416− 3577488.7381
−267132.4375−−266047.9670
5256103.2354− 5255972.6353
 =

−275.2965
−1084.4704
130.6002

And the length of the road is L = 1126.46m. (1)
[8 marks]
(b) The vector we are interested in is the reflector position with respect to the antenna
resolved into antenna axes,
raap = r
a
ep − raea = Caereep −Caereea
= Cae
(
reeb +C
e
bl
b
bp
)−Caereea (3)
Page 4 of 17 Continued overleaf...
Then, if we differentiate this expression,
r˙aap = C˙
a
e
(
reeb +C
e
bl
b
bp
)
+Cae
(
r˙eeb + C˙
e
bl
b
bp +C
e
b l˙
b
bp
)
− C˙aereea −Cae r˙eea
= C˙ae
(
reeb +C
e
bl
b
bp − reea
)
+Cae
(
r˙eeb − r˙eea + C˙eblbbp
)
(3)
The appropriate strapdown equations are,
C˙ae = −ΩaeaCae , C˙eb = CebΩbeb (2)
Substituting,
r˙aap = −ΩaeaCae
(
reeb − reea +Ceblbbp
)
+Cae
(
r˙eeb − r˙eea +CebΩbeblbbp
)
= vaeb − vaea +CabΩbeblbbp −Ωaea
(
raeb − raea +Cab lbbp
)
(4)
... as required.
[12 marks]
Page 5 of 17 Continued overleaf...
Q2.
The navigation equations expanded and resolved into the local geographic frame are,
v˙neb = C
n
b f
b
ib
+ gn
b
(Lb, hb) − (Ωnen + 2Ωeie) v˙neb
Using Trapeziodal integration to solve for velocity & position, integrate this equation over
one time-step to calculate the position at time tk when the conditions shown in Table Q.1
hold. Use the Somigliana model to calculate g, assume a sampling interval of ∆t = 0.1s
and assume the following initial DCM.
Cnb (tk) =

0.7170 −0.6696 0.1938
0.6905 0.7204 −0.0657
−0.0956 0.1810 0.9788

Variable Values Units
v˙neb(tk−1) [0, 0, 0] ms
−2
vneb(tk−1) [240, 0, 0] m/s
Lb(tk−1), λb(tk−1) [55.87858658,−4.676749498] deg
hb 4000 m
ωbib(tk) [0.1, 0.05,−0.2] rad/s
f b
ib
(tk) [0.0, 0.01,−0.05] ms−2
Table Q.2 - Previous INS estimates and current measurements
[20 marks]
Solution:
Step 1: Convert the measured specific force from body to NED axes,
fn
ib
= Cnb (tk)f
b
ib
(tk) =

−0.0164
0.0105
−0.0471
 (1)
Page 6 of 17 Continued overleaf...
Step 2: Calculate g using the Somigliana model.
g = 9.7803253359
(1 + 0.001931853 sin2 L)√
1− e2 sin2 L
=

0.0000
0.0000
9.8158
 (2)
Step 3: Construct the earth spin and transport rate skew-symmetric matrices.The trans-
port rate is given by the following equations,
ωnen =

vneb,E/ (RE(Lb) + hb)
−vneb,N/ (RN(Lb) + hb)
−vneb,E tan(Lb)/ (RE(Lb) + hb)

which means we need to calculate RE(Lb) and RN(Lb). From the datasheet,
RN(L) =
R0 (1 − e2)(
1− e2 sin2 L) 32 = 6379290.3051
RE(L) =
R0√
1− e2 sin2 L
= 6392818.6824 (2)
The transport equation skew-symmetric matrix is then,
Ωnen =

0.000000 −0.000000 −0.000038
0.000000 0.000000 −0.000000
0.000038 0.000000 0.000000
 (2)
The earth rotation rate skew-symmetric matrix is simply,
Ωeie =

0.000000 0.000060 0.000000
−0.000060 0.000000 −0.000041
−0.000000 0.000041 0.000000
 (1)
Step 4: We now have everything necessary to integrate the navigation equation. The
velocity of vehicle frame F b with respect to ECEF frame F e, resolved in the NED
frameF n is,
v˙ =

−0.0164
0.0395
9.7597
 (3)
Page 7 of 17 Continued overleaf...
Step 5: Using the Trapeziodal rule, we can integrate this new acceleration term to
obtain the velocity at time tk,
vneb = v
n
eb(tk−1) + ∆t
v˙neb(tk) + v˙
n
eb(tk−1)
2
=

239.9992
0.0020
0.4880
 (3)
Step 6: Using the velocity terms to calculate the rate of change of curvilinear coordi-
nates,
L˙b =
vneb,N
RN(Lb) + hb
λ˙b =
vneb,E
(RN(Lb) + hb) cosLb
h˙b = −vneb,D
First, get the old and new curvilinear coordinate angular velocities,
L˙b(tk−1) =
vneb,N(tk−1)
RN(Lb) + hb
= 0.000038
λ˙b(tk−1) =
vneb,E(tk−1)
(RN(Lb) + hb) cosLb(tk−1)
= 0.000000
h˙b(tk−1) = −vneb,D(tk−1) = 0.000000
L˙b(tk) =
vneb,N(tk)
RN(Lb) + hb
= 0.000038
λ˙b(tk) =
vneb,E(tk)
(RN(Lb) + hb) cosLb(tk−1)
= 0.000000
h˙b(tk) = −vneb,D(tk) = −0.487983 (3)
Finally, use the trapezoidal rule to get the updated position
p(tk) =

55.8786
−4.6767
3999.9756
 (3)
[20 marks]
Page 8 of 17 Continued overleaf...
Q3.
(a) GNSS systems make use of techniques developed for radio trilateration to determine
receiver position. Consider the simple situation of trilateration in 2D using 3 beacons,
each positioned at the following (x,y) locations,
B1(300, 250), B2(−50, 175), B3(20,−90)
If a receiver is placed at location uT = (30, 30), use an iterated linear approximation
algorithm to calculate the estimated receiver position when ∥δu∥ < 0.1. Assume a
nominal, initial position estimate of u = (0, 0). (10)
(b) A Kalman Filter is to be used to estimate the trajectories of the following three-state
model,
xˆ(tk−1) =

xˆ1
xˆ2
xˆ3
 , Φ =

0 0 1
1 ∆t 0
0 0 1
 , H = [0 1 0] ,
For this model, calculate the a´ posteriori estimated state vector xˆ+k assuming a sam-
pling interval of ∆t = 0.1s and using the following matrices, (10)
xˆ+(tk−1) =

0.40
1.50
−0.50
 , P+(tk−1) =

8.20 3.50 6.60
2.40 1.90 4.70
9.30 2.50 3.50

Q =

10 0 0
0 5 0
0 0 2
, R = 0.01
[20 Marks]
Page 9 of 17 Continued overleaf...
Solution:
(a) Begin by calculating true ranges
B1 =
√
(300− 30)2 + (250− 30)2 = 348.2815
B2 =
√
(−50− 30)2 + (175− 30)2 = 165.6050
B3 =
√
(20− 30)2 + (−90− 30)2 = 120.4159 (1)
Begin iterative algorithm
Iteration 1
Calculating range estimates from current position estimate
R1 =
√
(0.000− 300)2 + (0.000− 250)2 = 390.5125
R2 =
√
(0.000−−50)2 + (0.000− 175)2 = 182.0027 (1)
R3 =
√
(0.000− 20)2 + (0.000−−90)2 = 92.1954
Construct Jacobian
H =

0.00−300.00
390.51
0.00−250.00
390.51
0.00−−50.00
182.00
0.00−175.00
182.00
0.00−20.00
92.20
0.00−−90.00
92.20
 =

−0.768 −0.640
0.275 −0.962
−0.217 0.976
 (1)
Get range error
δR = (348.281− 390.512, 165.605− 182.003, 120.416− 92.195)
= (−42.23099,−16.39780, 28.22050) (1)
Calculate position update
δu = (HTH)−1HT δR
=
0.7127 0.0159
0.0159 2.2873
−1 −0.7682 0.2747 −0.2169
−0.6402 −0.9615 0.9762


−42.2310
−16.3978
28.2205

=
29.9296
30.5493
 (2)
Page 10 of 17 Continued overleaf...
Therefore, the new position estimate is,
u = (29.92963, 30.54925) (1)
and the norm of the position update vector is 42.7673 (1)
Iteration 2
Calculating range estimates from current position estimate
R1 =
√
(29.930− 300)2 + (30.549− 250)2 = 347.9894
R2 =
√
(29.930−−50)2 + (30.549− 175)2 = 165.0902
R3 =
√
(29.930− 20)2 + (30.549−−90)2 = 120.9575
Construct Jacobian
H =

29.93−300.00
347.99
30.55−250.00
347.99
29.93−−50.00
165.09
30.55−175.00
165.09
29.93−20.00
120.96
30.55−−90.00
120.96
 =

−0.776 −0.631
0.484 −0.875
0.082 0.997

Get range error
δR = (348.281− 347.989, 165.605− 165.090, 120.416− 120.958)
= (0.29208, 0.51478,−0.54156)
Calculate position update
δu = (HTH)−1HT δR
=
0.8435 0.1476
0.1476 2.1565
−1 −0.7761 0.4842 0.0821
−0.6306 −0.8750 0.9966


0.2921
0.5148
−0.5416

=
 0.0702
−0.5494

Therefore, the new position estimate is,
u = (29.99980, 29.99989)
Page 11 of 17 Continued overleaf...
and the norm of the position update vector is 0.5538
Iteration 3
Calculating range estimates from current position estimate
R1 =
√
(30.000− 300)2 + (30.000− 250)2 = 348.2817
R2 =
√
(30.000−−50)2 + (30.000− 175)2 = 165.6049
R3 =
√
(30.000− 20)2 + (30.000−−90)2 = 120.4158
Construct Jacobian
H =

30.00−300.00
348.28
30.00−250.00
348.28
30.00−−50.00
165.60
30.00−175.00
165.60
30.00−20.00
120.42
30.00−−90.00
120.42
 =

−0.775 −0.632
0.483 −0.876
0.083 0.997

Get range error
δR = (348.281− 348.282, 165.605− 165.605, 120.416− 120.416)
= (−0.00022, 0.00000, 0.00012)
Calculate position update
δu = (HTH)−1HT δR
=
0.8412 0.1495
0.1495 2.1588
−1 −0.7752 0.4831 0.0830
−0.6317 −0.8756 0.9965


−0.0002
0.0000
0.0001

=
0.0002
0.0001

Therefore, the new position estimate is,
u = (30.00000, 30.00000)
and the norm of the position update vector is 0.0002
Note: Most of the marks are awarded for knowing the process i.e. within the first
iteration (7). Three marks are awarded for the subsequent iterations
Page 12 of 17 Continued overleaf...
(b) In calculating the a´-posteriori estimate, we must first calculate the a´-priori estimate
in the prediction step.
xˆ−(tk) = Φxˆ
+(tk−1))
=

0.00 0.00 1.00
1.00 0.10 0.00
0.00 0.00 1.00


0.40
1.50
−0.50

=

−0.50
0.55
−0.50
 (1)
Next, we need to calculate the a´ priori error covariance matrix P−(tk).
P−(tk) = ΦP+(tk−1)ΦT +Q
=

13.50 9.55 3.50
7.07 13.81 7.07
3.50 9.55 5.50
 (3)
The equation for the Kalman Gain is
K(tk) = P
−(tk)HT
[
HP−(tk)HT +R
]−1
Consider the term inside brackets,
[
HP−(tk)HT +R
]
=
[
0.00 1.00 0.00
]
13.50 9.55 3.50
7.07 13.81 7.07
3.50 9.55 5.50


0.00
1.00
0.00
+ [0.01]
= 13.8190 (2)
The Kalman gain may now be calculated as follows,
K1 = P
−(tk)HT
[
HP−(tk)HT +R
]−1
=

13.50 9.55 3.50
7.07 13.81 7.07
3.50 9.55 5.50


0.00
1.00
0.00
 /13.8190
=

0.6911
0.9993
0.6911
 (2)
Page 13 of 17 Continued overleaf...
The final equation we require to calculate the a´-posteriori state estimate is,
xˆ1 = xˆ
−
1 +K1
(
y
1
−Hxˆ−1
)
Therefore, the a´ posteriori state estimate is,
xˆ+(tk) =

−0.5000
0.5500
−0.5000
+

0.6911
0.9993
0.6911
 [1.75− 0.5500]
=

0.3293
1.7491
0.3293
 (2)
[10 marks]
Page 14 of 17 Continued overleaf...
Navigation Systems Data Sheet
Useful Constants
R0 = 6, 378, 137m f = 1/298.257223563
RP = 6, 356, 752m e = 0.0818191908425
GM = 3986004.4188m3/s2 ωeie = 7.2921150× 10−5 rad/s
Earth Circumference = 40030.2km
Navigation Mathematics
Properties of Coordinate Transformation Matrices:
Cαβ =
(
Cβα
)−1
=
(
Cβα
)T
Cγα = C
γ
βC
β
α
CβαC
α
β = I
Properties of Skew-Symmetric Matrices:
Ωγβα =

0 −rγβα qγβα
rγβα 0 −pγβα
−qγβα pγβα 0

Ωδβα = C
δ
γΩ
γ
βαC
γ
δ
The Strapdown Equation:
C˙αβ = −CαβΩββα = −ΩαβαCαβ = CαβΩβαβ = ΩααβCαβ
Kinematics:
vγβα = C
γ
β r˙
β
βα
aγβα = C
γ
β r¨
β
βα
r¨γβα = −
(
ΩγβγΩ
γ
βγ + Ω˙
γ
βγ
)
rγβα − 2Ωγβγ r˙γβα + aγβα
Page 15 of 17 Continued overleaf...
Earth Surface Model:
RN(L) =
R0 (1 − e2)(
1− e2 sin2 L) 32
RE(L) =
R0√
1− e2 sin2 L
xeeb = (RE(Lb) + hb) cosLb cosλb
yeeb = (RE(Lb) + hb) cosLb sinλb
zeeb =
[
(1− e2)RE(Lb) + hb
)
sinLb
Somigliana Gravity Model:
g
0
(L) = 9.7803253359
(1 + 0.001931853 sin2 L)√
1− e2 sin2 L
ms−2
Inertial Navigation
Transport rate:
ωnen =

vneb,E/ (RE(Lb) + hb)
−vneb,N/ (RN(Lb) + hb)
−vneb,E tan(Lb)/ (RE(Lb) + hb)

Curvilinear coordinate rates:
L˙b =
vneb,N
RN(Lb) + hb
λ˙b =
vneb,E
(RN(Lb) + hb) cosLb
h˙b = −vneb,D
Global Navigation Satellite Systems
Taylor’s Theorem:
f(x, y) = f(x0, y0) +
∂f
∂x
(x− x0) + ∂f
∂y
(y − y0) + 1
2!
∂2f
∂x2
(x− x0)2
+
1
2!
∂2f
∂y2
(y − y0)2 + . . . h.o.t.
Moore-Penrose pseudo inverse:
Hδu = δr
⇒ δu = (HTH)−1HT δr
Page 16 of 17 Continued overleaf...
Integrated Navigation Systems
Gaussian Random variables:
E[x] =
∫ ∞
−∞
xp(x) dx = µ
var(x) = E[(x− µ)2] =
∫ ∞
−∞
(x− µ)2 p(x) dx = σ2x
p(x) =
1√
2πσ
e−
(x−µ)2
2σ2
The Kalman Filter:
xˆ−k = Φk−1xˆk−1
P−k = Φk−1Pk−1Φ
T
k−1 +Qk
Kk = P
−
kH
T
k
[
HkP
−
kH
T
k +Rk
]−1
xˆk = xˆ
−
k +Kk (zk −Hxˆk)
Pk = (I−KkHk)P−k
TimeUpdate
Measurement
Update