GNSS independent navigation using radio navigation equipment

(1)

navigation equipment

Pontus Törnberg

Space Engineering, master's level

2020

Luleå University of Technology

(2)

(3)

(4)

(5)

List of Figures

2 Coordinate frames ECEF, ECI and NED . . . 2

3 Coordinate system of Body frame . . . 5

4 Definition of bearing . . . 7

5 All four flight paths, flight path (blue), position of the radio station (x), the measurement range (red and dashed) . . . 14

6 Setup for the linear Kalman filter . . . 16

7 Flow chart of the algorithms KF and EKF . . . 18

8 Setup for the Extended Kalman filter . . . 19

9 Estimated (blue) and real (red) flight paths . . . 25

10 Position error in east direction (blue) and the corresponding σ-bounds (red) 26 11 Position error in north direction (blue) and the corresponding σ-bounds (red) . . . 27

(8)

1 Introduction

1.1 Motivation

When flying an aircraft both civilian and military, the navigation system is heavily dependant on Global Navigation Satellite system (GNSS) signals. It can therefore be a target for interference from hostile sources. If the GNSS signals is interfered, the navigation system will not be able to use the information from the GNSS signal to navigate. An alternative could be to use radio navigation stations on the ground. These radio equipment does not provide any information about the aircraft’s position like the GNSS, the information from these stations is relative distance and/or directions from the station to the aircraft.

1.2 Purpose

The purpose of this thesis is to study algorithms to estimate an aircraft’s position with different information from various radio stations and fusing the data with data from the inertial navigation system.

1.3 Limitations

(9)

2 Theory and Background

2.1 Coordinate Frames

This chapter will discuss different coordinate frames that is often used when working with navigation. An inertial system designer has multiple reference coordinate frames to work with in navigation, when computing the attitude, position and velocity. The choice of coordinate frame will largely depend on the requirements of the mission, diffi-culties during implementation, data storage and speed as well as the complexity of the navigation equations. When acceleration, velocity and position is expressed as vectors one can use vector operations to describe relations between system variables. However, these vectors has no meaning what so ever if they are not expressed with respect to a known reference coordinate frame. When one determines the motion of a vehicle it is necessary to relate the solution to a moving earth. One first defines a convenient iner-tial reference frame with respect to earth, then one determines the motion of both the earth and the vehicle with respect to this reference frame. The inertial system designer requires knowledge about initial position, initial velocity and initial orientation of the reference frame to be able to determine the future position, velocity and orientation. Often the desired outputs for an inertial navigation system is altitude h, longitude λ and latitude ϕ. The inertial navigation system must then be mechanized to provide these outputs directly. In general there are six fundamental coordinate frames that are of interest to the inertial system designer to describe motion relative to some reference frame. These frames are: true inertial, earth-centered earth-fixed, earth-centered inertial, navigation, wander-azimuth and body. These coordinate frames are unless specifically defined, right-handed, orthogonal, Cartesian frames. They differ in the relative motion between frames, the relative orientation of the axes and in the location of the origin. [1] This thesis will work in navigation and body frames but it can easily be expanded to work in other frames as well.

(10)

2.1.1 True Inertial Frame - I

This frame is the only frame that in which Newton’s Laws of motion are completely valid. However, Newton’s Laws are also valid in Galilean frames when they do not rotate with respect to one another and are translating in space. These Galilean frames are too impractical to use as a reference frame. Newton assumed it existed a reference frame that had an absolute motion of zero. Newton then considered this inertial frame to be one frame that had an absolute motion of zero. This made Newton’s Laws of motion valid in this reference frame. Controversies has always existed even when Newton was alive. Controversies if a reference frame like this existed or not led to theories of relativity. For Newtonian mechanics this is a special case. However, the true inertial frame is not very practical frame to use and are used as a visualization tool to represent other reference frames. [1]

2.1.2 Earth-Centered Inertial Frame - ECI

This frame has its origin at the center of the earth. It is not rotating relative to inertial space. The frame accelerates with respect to inertial space because the frame moves with the earth as it is moving and rotating around the sun. The inertial frame will appear to be having a rotational rate which will be a combination of earth’s position relative to the sun and that is the Julian day and earth’s rotational rate.

When one starts the navigation mode the x and y plane lies in the equatorial plane. The x axis are usually aligned pointing towards a star and with the z axis lies within earth’s rotational axis. Worth mentioning is that the frame does not rotate with the earth. In more detail, the axes are not inertial per se. The earth exhibits various modes of motion relative to the fixed space which is the axes fixed to the earth.

This frame is very useful because Newton’s laws are almost completely correct, for vehicles that navigates close to earth and in this frame computations of specific force is preformed. [1]

2.1.3 Earth-Centered Earth-Fixed frame - ECEF

This frame can also be referred to as the geocentric frame. It has its origin at the earth’s center, just like the inertial frame. However, unlike the inertial frame, the geocentric frame rotates with the earth and aligns with the inertial frame once every twenty four hour. The x axis is pointing towards the Greenwich Meridian which is also longitude 0° with 0° latitude. The z axis is pointing north towards the pole. The y axis is orthogonal to the x and z axis and is pointing 90° east longitude with 0° latitude.[1]

2.1.4 Navigation Frame

(11)

perpendicular to that sphere/ellipsoid. Often the y axis will be pointing towards east, the z axis down and the x axis will be pointing north. Most of the times the z axis will be pointing down but this will depend on the coordinate convention the designer of the system selects. It is worth mentioning that most inertial navigation systems is the local-level type. The platform is constrained with two of its axes in the horizontal plane. Historically, inertial navigation systems have been using local-level mechanization, due to simplifications of maintaining the alignment to the gravity vector of the platform constant and error compensation. The early local-level inertial navigation systems were north-slaved type. This led to both computational and hardware difficulties in the polar regions due to singularities.[1]

2.1.5 Wander-Azimuth Frame

The wander-azimuth frame is a special case of the navigation frame. This frame can also be called the computational frame. It has its origin at the inertial navigation systems location and it also coincident with the navigation frame’s origin. The wander-azimuth frame’s horizontal axes lies in a plane that is tangent to the navigation frame’s horizontal plane. The reference frame is maintained on a local level and can be defined with respect to earth’s frame with three successive Eulerian rotations. It is worth mentioning that the reference frame is measuring in the geodetic horizon plane, thus if the wander angle α is zero, then the reference frame will align with the geographic frame. When the longitude, latitude and wander angle are zero the reference frame will align with the earth-fixed frame. A lot of today’s inertial navigation systems are mechanized in this reference frame. [1]

2.1.6 Body Frame

(12)

Figure 3: Coordinate system of Body frame

2.2 Inertial measurement unit

This chapter will discuss two of the most crucial sensors in inertial navigation. The inertial measurement unit is an electronic device that consists of two different sensors, accelerometers and gyroscopes. In more advanced IMU’s magnetometers can be found in the unit. The IMU measures and outputs a body’s specific force and angular rate.[2] Combining this unit with a computer, it can calculate the position by dead reckoning and can be called an inertial navigation system. All sensors will be influenced by errors, both systematic and random errors. These sensors are no exception.

2.2.1 Accelerometer

Accelerometers measures the specific force, which is defined as the total non gravita-tional acceleration the accelerometer is influenced by. An accelerometer measures in one direction, multiple accelerometers are therefore placed orthogonally to measure the specific force in each axis in a body. From these measurements velocity and position can be estimated by integrating the acceleration. Because of this integration the estimates will drift due to the integration of measurement errors. [3]

2.2.2 Gyroscope

(13)

2.3 Schuler oscillations

One could describe the generation of oscillations in an inertial navigation system in a sim-ple way by studying the behaviour of an inertial navigation system on a stable platform in response to a tilt error. Consider the platform control system to be programmed to keep the stable platform horizontal with one of the axes of the platform. For example use the North,East and Down coordinate system. Then let the x-axis point towards north, the y-axis pointing towards east and the z-axis pointing down in the direction of the local nadir. [4] If there is a tilt error the accelerometer measuring eastward acceleration will then measure an acceleration as

f = −g ∗ sin(ϕ) . (2.1)

Where g is the gravity and ϕ is the tilt error. Then the inertial navigation system will integrate this f to obtain velocity. Which leads to error in velocity that will lead to an error in position. By the derivation in Appendix B one can obtain the equation for the time period for these oscillations as

TIN S= 2 ∗ π

s RE

g , (2.2)

where RE is the radius of the Earth. This period is known as the Schuler period. It is

very well known in the inertial navigation community. 2.4 Radio Navigation Equipment

This chapter will discuss the three different types of radio navigation equipment that is used in this thesis.

2.4.1 Distance measurement equipment (DME)

Distance measurement equipment measures the distance between a ground facility and an aircraft. The distance is determined by measuring the total travel time for the radio signal to go from the aircraft to the ground-system and reply to the aircraft. Then ad-justing for the known reply delay and then convert time to distance. This process starts with the avionics ”interrogating” the ground-system, the ground-system shall then de-code the valid interrogations and reply with a known delay, the avionics on the aircraft receives and decodes the reply and calculates the distance.[5] If one has two or more measurements from different DME stations one can triangulate the position of the air-craft.

2.4.2 Very High Frequency omnidirectional range (VOR)

(14)

identifier in Morse code and information that the aircrafts avionics/receiver can derive to bearing from the beacon and the aircraft.[6] There are different kinds of bearings, in this thesis it is defined as the angle between the radiostation’s north and the aircraft and it is shown in fig.4.

Figure 4: Definition of bearing

The station transmits two signals, a directional signal which rotates with 0.03 seconds per cycle through 360° and a signal called the reference phase with a constant omnidirectional signal. The two signals are only in phase once per rotation, this occurs when direction signal is aligned with the magnetic north. If one has two or more measurements from different VOR stations one can triangulate the position of the aircraft.[7]

2.4.3 TACAN

TACAN is an acronym for tactical air navigation system. This system provides an aircraft information about distance and bearing from the beacon to the aircraft. TACAN is primarily used for military services and can operate from ground-to-air or air-to-air mode. This thesis will focus on ground-to-air mode.

2.5 Bayesian filtering

(15)

If one has a parameter xk that is affected by noise and is non stationary it can be

estimated by dynamic model and observations. The observation model can be describes as the following probability density function

p(zk|xk) , (2.3)

and the probability density function for the dynamic model can be described as

p(xk+1|xk) . (2.4)

Where zk is the measurement, xk is the state at the time k. When using the Bayesian

approach in filtering one estimates the state xk, it is done by computing two distributions

p(xk|z1:k) and p(xk+1|z1:k) with the model expressed in (p(xk+1|xk)). Using Baye’s rule

the posterior distribution can be obtained and expressed as p(xk|z1:k) =

p(zk|xk) ∗ p(xk|z1:k−1)

p(zk|z1:k−1)

. (2.5)

Knowing the probability of the measurement p(zk|xk) and the prior distribution p(xk|z1:k−1)

one can compute the normalization constant as p(zk|z1:k−1) =

Z

<

p(zk|xk) ∗ p(xk|z1:k−1)dxk. (2.6)

Where IR is a set numbers that includes all real numbers from −∞ to ∞. The prediction is obtained by using the total law of probability and can be expressed as

p(xk+1, xk|z1:k) = p(xk+1|xk) ∗ p(xk, z1:k) . (2.7)

If one integrates both sides of Equation (2.7) with respect to xk one gets the following

expression for the time update p(xk+1|z1:k) =

Z

<

p(xk+1|xk) ∗ p(xk|z1:k)dxk. (2.8)

The measurement update is expressed in Equation (2.5), one can use this to reduce the uncertainty in the estimate of the state x at the time k. The prediction step of the filter, also known as the time update is used to predict the state x at the time k + 1 but the uncertainty will increase. [9]

Z.Chen[10] states that an optimal filter is said to be ”optimal” but only in some specific sense. Which means one shall define a criteria that measures optimality. Chen gives potential criterion for measuring optimality such as the minimum mean-squared error. Which can be defined in terms of the trace of the state-error covariance or filtering or prediction error

IE[||xk− ˆxk||2|zk] =

Z

||x_k− ˆxk||2∗ p(xk|zk)dxk . (2.9)

This aims to find the conditional mean for ˆ

xk= IE[xk|zk] =

Z

(16)

2.6 Kalman Filter

This chapter will discuss the Kalman filter for both linear and non linear systems. For the Kalman filter to work properly it requires that all error states are modeled as a zero mean noise process and has known power spectral densities, variances and time correlation parameters. These error quantities has to be estimated and the associated noise in the measurement are all random processes and its correlation structure can be assumed to be known. The filter obtains estimates of the states of the stochastic process and is described as linear or linearized mathematical model. It is done by capitalizing on the known correlation structure of the different processes involved and measurements that is described by linear combinations of the error states. For this to be accomplished both the measurement process and error propagation in time are expressed in vector form. By using a linear matrix, it could be a convenient way to keep track of relative complex relationships between the quantities of interest. In digital signal processing applications this is one of the main reasons that distinguishes a Kalman filter. When assuming Gaussian noise distributions, the filter minimizes the mean square error in the estimates of the state variables. It is reasonable to assume Gaussian noise in navigation applications because the arising noise effects in the measurements are usually a summation of many smaller random errors. According to the central limit theorem of statistics, this tends towards the Gaussian probability distribution. Many navigation systems has used and implemented Kalman filters by implicitly making this assumption. [3]

2.6.1 State-space model

This thesis will use a discrete time model to describe the system. The system is described by these following equations

xk+1= F ∗ xk+ wk (2.11)

zk = Hk∗ xk+ ek . (2.12)

Where xk+1 describes the state of the system and zk describes the measurements or

observations in terms of the system’s state. F is a state transition matrix and it is describing how the previous state transitions into the new state and wk is the process

noise. Hk is the measurement matrix and it is describing the the state vector relates to

the measurement and ek is the measurement noise. Both wk and ek are assumed to be

Gaussian with zero mean.[11]

2.6.2 Linear Kalman filter - KF

As being said in section (2.6), this filter obtains estimates of the state and is described with a linear model. Using the state space model described in section (2.6.1) the following algorithm is the best linear filter.

Measurement update:

(17)

Kk= Pk|k−1∗ HTk∗ S−1k (2.14) k= zk− Hk∗ xk|k−1 (2.15) xk|k= xk|k−1+ Kk∗ k (2.16) Pk|k= Pk|k−1− Pk|k−1∗ HTk∗ S−1k∗ Hk∗ Pk|k−1 (2.17) Time update: x_k+1|k = F ∗ x_k|k (2.18) P_k+1|k = F ∗ P_k|k∗ FT _{+ Q .} _(2.19)

Where Kk is called the Kalman gain, Sk is called the innovation covariance, k is called

the innovation, P is the covariance of the state and as mentioned before x is the state of the system.

2.6.3 Extended Kalman filter - EKF

The algorithm described in section (2.6.2) does not work for non linear systems. There-fore one invented the Extended Kalman filter to estimate non linear systems. The discrete time state space model for the EKF is described by the following equations

xk+1|k = f (xk|k, wk) (2.20)

zk= H(xk) + ek. (2.21)

To linearize the state the Extended Kalman filter uses Taylor expansion to get a lin-earized model of the non linear models f and H and then one can use the algorithm mentioned in section (2.6.2). The Extendended Kalman filter algorithm is described by the following equations

Where H_k0 is the Jacobian matrix of H(xk) and is defined as

(18)

(19)

3 Modeling of sensors and filter setup

3.1 Sensors

This section will explain how the sensors that is being used are modeled. 3.1.1 Accelerometer

This thesis will model the accelerometer as

ak= at+ ba+ ea. (3.1)

Where ak is the measured acceleration, at is the true acceleration, ba is the bias error

and ea is the noise error.

Under the assumption that there is zero errors and by integrating the acceleration one estimates the velocity as the following equation

vk= vk−1+ ak−1∗ ∆t . (3.2)

Integrating one more time one estimates the position as the following equation xk = xk−1+ vk−1∗ ∆t +

ak−1∗ ∆t2

2 . (3.3)

3.1.2 Gyroscope

By assuming zero errors for the gyroscope one ends up with the following sensor model

ωm = ωt, (3.4)

where ωm is measured angular rate and ωt is true angular rate.

By integrating ω one estimates the attitude as the following equation

Ωk= Ωk−1+ ωk−1∗ ∆t , (3.5)

where Ωk is the attitude in euler angles roll φ, pitch θ, yaw ψ.

3.1.3 DME

The DME model will be described as the following equation ˜

rDM E= rDM E+ δrDM E , (3.6)

where ˜rDM E is the distance that is sent to the aircraft. rDM E is the true distance and

is defined as rDM E = q (pNreal− posN) 2_{+ (p} Ereal − posE) 2 _, _(3.7)

where pNreal and pEreal is the real simulated position in north and east coordinates

respectively, posN and posE is the position of the DME station in north and east

coor-dinates.

δrDM E is the measurement errors and it can be expressed as[12]

(20)

3.1.4 VOR

The VOR model will be described as the following equations ˜

βV OR= βV OR+ δβV OR, (3.9)

where ˜βV OR is the bearing angle that is sent to the aircraft. βV OR is the true bearing

and is defined as βV OR=        −π 2 − arctan( p_Nreal−posN

p_Ereal−posE) if pEreal − posE < 0

π

2 − arctan(

p_Nreal−posN

p_Ereal−posE) if pEreal − posE ≥ 0

0 else

,

and as previously mentioned pNreal and pEreal is the real simulated position in north and

east coordinates respectively, posN and posE is the position of the VOR station in north

and east coordinates. δβV OR is the measurement errors and is defined as

δβV OR∼ N (0, σβ) . (3.10)

3.1.5 TACAN

This thesis will model the TACAN measurements as the following equations ˜

βT ACAN = βT ACAN + δβT ACAN , (3.11)

where ˜βT ACAN is the bearing angle that is sent to the aircraft, βT ACAN is the true

bearing which is defined as

βT ACAN =        −π 2 − arctan( p_Nreal−posN

p_Ereal−posE) if pEreal− posE < 0

π

2 − arctan(

p_Nreal−posN

p_Ereal−posE) if pEreal − posE ≥ 0

0 else

,

where pNreal, pEreal is the real simulated position in north and east coordinates

respec-tively, posN and posEis the position of the TACAN station in north and east coordinates.

δβT ACAN is the measurement errors plus the white noise when sending the information

to the aircraft which is defined as

δβT ACAN ∼ N (0, σβ) . (3.12)

The range that the TACAN station measures can be described as ˜

rT ACAN = rT ACAN + δrT ACAN , (3.13)

where ˜rT ACAN is the range measurement that is sent to the aircraft. rT ACAN is the true

distance and can be expressed as rT ACAN =

q

(pNreal− posN)2+ (pEreal− posE)2, (3.14)

where the parameters pNreal, pEreal,posN and posE is described previously. δrT ACAN is

the measurement errors and can be expressed as[12]

(21)

3.2 Path generator

For this thesis a path generator is set up to generate different flight paths. In flight path one the aircraft is traveling in the east direction towards the radio station. In flight path two the aircraft is traveling east towards the radio station and performs a 360° turn once every 60 seconds. In flight path three the aircraft is traveling east but starts to roll immediately and changing its heading towards north. In flight path four the aircraft is traveling east towards the radio station but goes back and forth in the north direction like a sinus wave. The four different flight paths can be seen in fig. 5.

(a) Flight path 1 (b) Flight path 2

(c) Flight path 3 (d) Flight path 4

Figure 5: All four flight paths, flight path (blue), position of the radio station (x), the measurement range (red and dashed)

(22)

maneuver. Which changes the heading of the plane. The heading is the direction where the nose is pointing and in this thesis it is also the direction that the plane is traveling. It can be described with the following equation from [1]

˙

ψ = g ∗ tan(φ)

v , (3.16)

where ˙ψ is the angular rate in heading, g is the gravity constant, φ as being said in section (3.1.2) is the roll angle and v is the speed of the aircraft. One can integrate

˙

ψ to get the heading angle defined as ψ. One can obtain the velocity in the different directions by adding the velocity to the initial velocity vector v. This is done by using a rotational matrix and multiplying it with the measured acceleration, in north, east and down. It is described as the following equation

ˆ

v = v + Rh(ψ) ∗ a ∗ t m/s, (3.17)

where v is the initial velocity vector and is defined as

v =   0 100 0  m/s, (3.18)

and Rh(α) is the rotational matrix defined as

Rh(ψ) =   cos(ψ) sin(ψ) 0 −sin(ψ) cos(ψ) 0 0 0 1   . (3.19)

The gravity is subtracted by using an additional rotational matrix defined as

Rr(φ) =   1 0 0 0 cos(φ) sin(φ) 0 −sin(φ) cos(φ)   . (3.20)

Multiplying Rr(φ) with the gravity vector g which is defined as

g =   0 0 9.81  m/s2. (3.21)

One can obtain the vector ˆg as the following

(23)

With the ˆg vector one can easily subtract the gravity from the measurements from the accelerometer and the remaining values from the accelerometer is acceleration, noise and the effects of Schuler oscillations. Which leaves the acceleration a as

ˆ a = ak+ w − g Re ∗ sin( r g Re ∗ t) , (3.23)

where w is noise from the accelerometer, g is the gravitational constant and Re is the radius of the Earth. Then by integrating Equation (3.2) one obtains the position of the aircraft and is determined by Equation (3.3) but in continuous time and it is defined as the following equation

x = x0+ ˆv ∗ t + ˆa ∗

t2

2 , (3.24)

where x0 is the starting position of the aircraft.

3.3 Converted TACAN measurements

As it can be seen in section (3.1.5) the TACAN measurements are non linear which are in region of the Extended Kalman filter. However, because one knows the position of the TACAN station one can derive a pseudo position of the aircraft.

Figure 6: Setup for the linear Kalman filter

One can therefore use the linear Kalman filter when dealing with TACAN measurements. Which gives the setup as seen in fig. 6. Assume the TACAN station has a position

pT ACAN =pNT ACAN

pET ACAN

. (3.25)

Because one knows that there is errors in the measurements from the TACAN station one can define the pseudo position as

˜ pT ACAN = pNtrue+ δpNT ACAN pEtrue+ δpET ACAN . (3.26)

Knowing the bearing and the distance to the aircraft one can compute the pseudo posi-tion as

˜

pT ACAN =

pNT ACAN + ˜rT ACAN ∗ cos( ˜βT ACAN)

pET ACAN + ˜rT ACAN ∗ sin( ˜βT ACAN)

(24)

The position obtained from the inertial navigation system can be defined as ˜ pIN S = pNIN S pEIN S =pNtrue+ δpNIN S pEtrue+ δpEIN S . (3.28)

Let the state vector x be defined as

x =     δpNIN S δpEIN S δvNIN S δvEIN S     . (3.29)

The measurement zk used is obtained by subtracting the pseudo position from the

esti-mated position from the inertial navigation system which can be written as

zk= ˜pIN S− ˜pT ACAN , (3.30)

which in fact is

zk= ptrue+ δpins− (ptrue+ δpT ACAN) . (3.31)

As it can seen in Equation (3.31) the true position ptrue cancel each other and leaves the

following expression for zk

zk= δpIN S− δpT ACAN . (3.32)

Which has been said previously is the error from the inertial navigation system minus the measurement error from the TACAN station.

Recall the state space model mentioned in section (2.6.1)

xk+1= F ∗ xk+ G ∗ wk (3.33)

zk= H ∗ xk+ ek. (3.34)

The state transition matrix F can be derived and expressed as

F =     1 0 T 0 0 1 0 T 0 0 1 0 0 0 0 1     . (3.35)

As well as the state transition matrix F , the measurement matrix H can be derived expressed as

H =1 0 0 0

0 1 0 0

, (3.36)

and wk is the noise measured from the accelerometer and it can be assumed to be

Gaussian noise with zero mean

(25)

The process noise covariance matrix can be expressed as Qk= Cov(wk) = G ∗ σ2 w 0 0 σ2_w ∗ GT , (3.38) where G is defined as G =     T2 2 0 0 T₂2 T 0 0 T     . (3.39)

ek is the noise measured and it is assumed to be uncorrelated and is as well modeled as

Gaussian noise with zero mean which gives the following measurement covariance matrix

R =cos( ˜β) −sin( ˜β) sin( ˜β) cos( ˜β) ∗σ 2 N 0 0 σ_E2 ∗cos( ˜β) −sin( ˜β) sin( ˜β) cos( ˜β) T , (3.40)

where σ_N2 and σ_E2 are defined as the variance in north and east position respectively and ˜

β is the measured bearing. However, this is only an approximation because the pseudo position are correlated and biased. A proof of this can be seen in [13]. It is worth mentioning that the paper also shows a way to account for this. However, this thesis will not account for these errors in the covariance matrix for the measurement.

Then one can implement the algorithm mentioned in section (2.6.2).

Figure 7: Flow chart of the algorithms KF and EKF

As it can be seen in fig. 7 if there is no measurement available the filter should only make a time updates until next measurement is available.

3.4 DME measurements

(26)

used when one need to use the Extended Kalman filter. This setup is being used for all measurements except for the converted TACAN measurement scenario.

Figure 8: Setup for the Extended Kalman filter Let the state vector be defined as

x =         pNest pEest vNest vEest aNest aEest         . (3.41)

Because the non linearity in the DME measurements the measurement vector zk is

defined as mentioned in section (2.6.3)

zk= H(xk) + ek , (3.42)

where ek is the noise measured, it is assumed to be uncorrelated and white Gaussian

noise. The measurement is obtained with measurements from the DME station and the accelerometers and the DME measurements are as mentioned in section (2.4.1) the distance to the aircraft. This will provide the following measurement noise covariance matrix R =           σ2_˜_r 0 0 0 0 0 0 0 σ2_p_N 0 0 0 0 0 0 0 σ_p2_E 0 0 0 0 0 0 0 σ2 vN 0 0 0 0 0 0 0 σ2_v_E 0 0 0 0 0 0 0 σ_a2_N 0 0 0 0 0 0 0 σ_a2 E           , (3.43)

(27)

H(xk) is the measurement function vector and it is defined as H(xk) =           ˜ rDM E pNIN S pEIN S vNIN S vEIN S aNIN S aEIN S           , (3.44)

where ˜rDM E and it is defined in section (3.1.3). pNIN S and pEIN S is defined in the

previous section (3.3). vNIN S and vEIN S can be described in the same way as the inertial

navigation system’s estimated position such as vIN S = vNIN S vEIN S =vNest+ δvNIN S vEest+ δvEIN S , (3.45)

where δvNIN S and δvEIN S is the velocity error in north and east direction respectively.

aNIN S and aEIN S can be described in the same way as the inertial navigation system’s

estimated position and velocity as aIN S= aNIN S aEIN S =aNest+ δaNIN S aEest+ δaEIN S , (3.46)

where δaNIN S and δaEIN S is the acceleration error in north and east direction

respec-tively.

However, this thesis uses a linear motion model so the state x is still described as mentioned in section (2.6.1). Therefore there is no use to linearize F. Knowing that the motion model is linear one can therefore define the state transition matrix as

F =          1 0 T 0 T₂2 0 0 1 0 T 0 T₂2 0 0 1 0 T 0 0 0 0 1 0 T 0 0 0 0 1 0 0 0 0 0 0 1          . (3.47)

(28)

As being said in section (2.6.3) one needs to derive the Jacobian matrix for H(xk)

because the measurements are non linear. Computing this Jacobian matrix with the expressions for the measurements and knowing the position of the DME stations one obtains the following Jacobian matrix

H0 =            ∆pN √ (∆pN)2+(∆pE)2 ∆pE √ (∆pN)2+(∆pE)2 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1            , (3.50)

where ∆pN and ∆pE is the difference between the estimated true position and the DME

station’s position, defined as

∆pN = pNest− pNDM E (3.51)

∆pE = pEest− pEDM E . (3.52)

Then implement the algorithm mentioned in section (2.6.3) and as mentioned previously and seen in fig. 7 if there is no measurement available the filter should only make a time updates until next measurement is available.

3.5 VOR measurements

As mentioned in section (3.1.4) the VOR measurements are non linear as well, therefore the use of the Extended Kalman filter is needed. As well for these measurements the state vector x and the process noise matrix is defined as in section (3.4). Being men-tioned earlier the measurements from the VOR station are non linear and therefore the measurement model zk is defined as the previous section (3.4). The noise measured ek

can as well be assumed to be uncorrelated and Gaussian with zero mean. The mea-surement is obtained with meamea-surements from the VOR station and the accelerometers and the VOR measurements are described in section (3.1.4) the bearing angle to the air craft. This will provide the following measurement noise covariance matrix

R =            σ2_˜ β 0 0 0 0 0 0 0 σ_p2_N 0 0 0 0 0 0 0 σ2_p_E 0 0 0 0 0 0 0 σ2_v_N 0 0 0 0 0 0 0 σ_v2_E 0 0 0 0 0 0 0 σ2_a N 0 0 0 0 0 0 0 σ2_a_E            , (3.53)

(29)

H(xk) is the measurement function vector and can be expressed as H(xk) =           ˜ βV OR pNIN S pEIN S vNIN S vEIN S aNIN S aEIN S           , (3.54)

where ˜βV ORis the measured bearing angle and it is defined in section (3.1.4). pNIN S and

pEIN S is defined in the section (3.3), vNIN S, vEIN S, aNIN S and aEIN S is defined in the

previous section. As mentioned before the motion model is linear and therefore the state transition matrix F and the process noise matrix Q can be expressed as the previous section (3.4). Knowing that the measurements are non linear for the VOR measurements as well one needs to compute the Jacobian matrix for the measurement vector H(xk).

It is defined by the following matrix

H0 =            −∆pE (∆pN)2+(∆pE)2 ∆pN (∆pN)2+(∆pE)2 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1            , (3.55)

where ∆pN and ∆pE is defined as the difference between the VOR station and the

aircraft. It can be expressed as

∆pN = pNest− pNV OR (3.56)

∆pE = pEest− pEV OR . (3.57)

Then again implement the algorithm mentioned in section (2.6.3) and being explained previously and seen in fig. 7 if there is no measurement obtained the filter should only make a time updates until next measurement is obtained.

3.6 TACAN measurements

As being said previously the TACAN measurements are non linear therefore one need to use the Extended Kalman filter. Let the state vector x and the process noise matrix be defined as in section (3.4). Knowing that the TACAN measurements are non linear one can describe the measurement vector zk as mentioned in section (2.6.3) as

(30)

ek is as mentioned in previous section the noise measured and it is assumed to be

uncorrelated and Gaussian with zero mean. However this measurement is obtained from the accelerometers and the TACAN station. The TACAN measurements are as mentioned in section (3.1.5) the distance and the bearing angle to the aircraft. This provides the following measurement noise covariance matrix

R =              σ_r2_˜ 0 0 0 0 0 0 0 0 σ2_˜ β 0 0 0 0 0 0 0 0 σ2 pN 0 0 0 0 0 0 0 0 σ2_p_E 0 0 0 0 0 0 0 0 σ2_v_N 0 0 0 0 0 0 0 0 σ_v2_E 0 0 0 0 0 0 0 0 σ2_a_N 0 0 0 0 0 0 0 0 σ2_a_E              , (3.59)

in which σ2_˜_r is the variance in distance, σ2_˜

β is the variance in the bearing angle. σ 2 pN,

σ_p2_E, σ_v2_N, σ_v2_E, σ_a2_N and σ2_a_E is defined in section (3.4).

H(xk) is the measurement function vector and can be expressed as

H(xk) =             ˜ rT ACAN ˜ βT ACAN pNIN S pEIN S vNIN S vEIN S aNIN S aEIN S             , (3.60)

where ˜rT ACAN and ˜βT ACAN is defined in section (3.1.5). pNIN S, pEIN S, vNIN S, vEIN S,

aNIN S and aEIN S is defined in the previous sections (3.3) and (3.4). As explained in

section (3.4) the motion model is linear therefore one can express the state transition matrix F and the process noise matrix Q is defined as in section (3.4). Mentioned as well in section (2.6.3) one needs to derive the Jacobian matrix for H(xk) + ek because

(31)

following Jacobian matrix H0 =               ∆pN √ (∆pN)2+(∆pE)2 ∆pE √ (∆pN)2+(∆pE)2 0 0 0 0 −∆pE (∆pN)2+(∆pE)2 ∆pN (∆pN)2+(∆pE)2 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1               , (3.61)

where ∆pN and ∆pE is the difference between the estimated true position and the

TACAN station’s position and is defined as

∆pN = pNest− pNT ACAN (3.62)

∆pE = pEest− pET ACAN . (3.63)

(32)

4 Results

4.1 Converted TACAN measurements

In fig. 9 it can be observed that the linear Kalman filter estimates the position fairly well. There are no estimates that deviates very much from the true position.

Figure 9: Estimated (blue) and real (red) flight paths

(33)

(34)

Figure 11: Position error in north direction (blue) and the corresponding σ-bounds (red) 4.2 DME measurements

(35)

(36)

(37)

Figure 14: Position error in north direction (blue) and the corresponding σ-bounds (red) 4.3 VOR measurements

(38)

(39)

(40)

Figure 17: Position error in north direction (blue) and the corresponding σ-bounds (red) 4.4 TACAN measurements

(41)

(42)

(43)

(44)

5 Discussion

The converted TACAN measurement algorithm estimates is good in both north and east direction. However, the filter could use a better tuning because as one can see in fig. 10a and fig. 10c some estimates is not within the σ-bounds. This is done by changing some parameters of the R and Q matrices. Worth mentioning, this method will only work with TACAN measurements.

As mentioned in section (4.2) the filter for DME measurements estimates the easting position very well and it estimates the northing position very poorly, except for Flight path 3. It appears that the filter estimates the position in the direction of the DME station very well when the aircraft is traveling towards the station. However, the filter makes poor estimates in the orthogonal direction of the traveling direction. As it can be seen in fig. 14 there is poor estimates for the position the north direction, except for Flight path 3. Which indicates that the algorithm for DME measurements can be sufficient enough to estimate a position but only for one direction. The estimates follows the σ-bounds very well and as it can be seen in fig. 13 and fig. 14 the estimates is well within these bounds.

In section (4.3) it is mentioned that the filter estimates the easting position very poorly and it estimates the northing position fairly well, except for Flight path 3 in which it is the other way around. It appears that the filter estimates the position in the orthogonal direction of the traveling direction very well when the aircraft is traveling towards the VOR station. Just as in the DME case the estimates of the filter seems to follow the σ-bounds very well and is within these bounds. However, as it can be seen in fig. 16d the σ-bounds are quite large.

Looking into both the DME and VOR algorithm. It can be seen that the VOR al-gorithm preforms better in one direction and the DME alal-gorithm in the other direction. This suggest that if one combine these measurements, it could result in better estimates in both directions. The estimates could improve as well by adding more VOR or DME measurements from different stations at different locations.

(45)

6 Conclusion

6.1 Results

The converted TACAN measurement and TACAN measurement algorithm can be used for navigation purposes by its own measurements. However, the DME and VOR mea-surement algorithms need to be combined or one need to change the flight path depending on which station the aircraft is approaching, if there is only one radio station. If there would be multiple stations in different directions this would improve the performance for these filters. All of the filter could use some better tuning to get the optimal filter, but it is not necessary.

6.2 Future work

One could improve the flight simulator to be more realistic. It is now somewhat in three dimensions but in the measurements the altitude dimension is remove due to the assumption that the aircraft is on a constant altitude. Therefore one could expand it to three dimensions and work in different coordinate systems. For example by switch-ing between the ECEF, NED and body frames and make the proper calculations in the proper coordinate system. However, one will need airdata to estimate the altitude, which will increase the computational power. The inertial measurement unit can be improved by adding noise to the gyroscopes and therefore estimate the aircrafts attitude with an Extended Kalman filter for example. One then need to add more states to the state space model and it will get more computationally heavy. With the improved inertial measurement unit and the improved flight simulator one can implement a more correct error model for the Schuler feedback.

(46)

[1] G. M.Siouris, Aerospace Avionics Systems, A Modern Synthesis. Academic Press, 1993, isbn: 0-12-646890-7.

[2] E. Foxlin, “Pedestrian tracking with shoe-mounted inertial sensors”, IEEE Com-puter Graphics and Applications, vol. 25, no. 6, pp. 38–46, 2005.

[3] W. F. Myron Kayton, Avionics navigation systems. John Wiley Sons, 1997, isbn: 0-417-54795-6.

[4] P. G. Savage, “Schuler oscillations”, Jun. 2014. [Online]. Available: http : / / strapdownassociates.com/Schuler%20Oscillations.pdf.

[5] R. Lilley and R. Erikson, “Dme/dme for alternate position, navigation, and timing (apnt)”, APNT White Paper, 2012.

[6] “Vhf omnidirectional radio range (vor)”, [Online]. Available: https : / / www . skybrary.aero/index.php/VHF_Omnidirectional_Radio_Range_(VOR), (ac-cessed: 10.05.2020).

[7] “Vhf omnidirectional range - vor”, [Online]. Available: https://mediawiki.ivao. aero/index.php?title=VHF_Omnidirectional_Range_- _VOR_(Instrument), (accessed: 10.05.2020).

[8] J. G.Minkler, Theory and Application of Kalman Filtering. Magellan Book Com-pany, 1993, isbn: 0-9621618-2-9.

[9] C. Rapp, Unmanned aerial vehicle positioning using a phased array radio and gnss independent sensors, 2019. [Online]. Available: http://urn.kb.se/resolve?urn= urn:nbn:se:liu:diva-157414.

[10] Z. Chen, “Bayesian filtering: From kalman filters to particle filters, and beyond”, Statistics, vol. 182, Jan. 2003. doi: 10.1080/02331880309257.

[11] J. B. M. Brian D. O. Anderson, Optimal Filtering. Prentice Hall, 1979, isbn: 0-13-638122-7.

[12] S. Lo, Y. Chen, P. Enge, B. G. K. Peterson, R. Erikson, and R. Lilley, “Distance measuring equipment accuracy performance today and for future alternative posi-tion navigaposi-tion and timing (apnt)”, Proceedings of the 26th Internaposi-tional Techni-cal Meeting of the Satellite Division of The Institute of Navigation (ION GNSS+ 2013), Nashville, TN, September 2013, pp. 711-721, 2013.

(47)

A.1 Radio station’s position

For flight path 1,2,4 the station’s position is set to posstation= 0 888960 m For flight path 3 the station’s position is set to

posstation=

888960 −5000

m

A.2 Filter Parameters • x0 is the initial state.

• P0 is the initial error covariance matrix.

• Q is the process noise matrix. • R is the measurement noise matrix.

• f is the update frequency of the filter and is set to 10 Hz A.2.1 Converted TACAN measurement

For all flight paths

(48)

For flight path 1,2,4 x0 =         5000 0 0 100 0 0        

For flight path 3

x0 =         5000 0 100 0 0 0        

(49)

For flight path 1,2,4 x0 =         5000 0 0 100 0 0        

For flight path 3

x0 =         5000 0 100 0 0 0        

(50)

For flight path 1,2,4 x0 =         5000 0 0 100 0 0        

For flight path 3

x0 =         5000 0 100 0 0 0        

(51)

Derivation of the Schuler period

Let the acceleration f measured by tilt error be defined as the following equation f = −g ∗ sin(ϕ)

. Using arc length formulas one can obtain the linear acceleration a along an arc. s is the arc length

s = l ∗ ϕ . Differentiating s one obtains the velocity v

v = ds dt = l ∗

dϕ dt . Differentiating v one arrives at the linear acceleration a

a = dv dt = l ∗ d2ϕ d2_t . Let a be equal to f a = f , l ∗d 2_ϕ d2_t = −g ∗ sin(ϕ)

. Adding f to both sides of the equation one arrives at the differential equation that describes the motion of a simple pendulum

l ∗d

2_ϕ

d2_t + g ∗ sin(ϕ) = 0

. Using small angle approximation

sin(ϕ) ≈ ϕ , one obtains the equation for a harmonic oscillator

d2ϕ d2_t +

g

l ∗ ϕ = 0

. Using initial conditions as ϕ(0) = ϕ0 and dϕ_dt(0) = 0 will give the following solution

ϕ(t) = ϕ0∗ cos(

(52)

ωIN S=

r g l

. The formula for angular velocity is defined as the following equation ωIN S =

2π TIN S

. Let the two equations above be equal r g

l = 2π TIN S

, and solve for TIN S

TIN S = 2π ∗

s Re

GNSS independent navigation using radio navigation equipment

navigation equipment

Pontus Törnberg

Space Engineering, master's level

2020

Contents

List of Figures

1

Introduction

2

Theory and Background

3

Modeling of sensors and filter setup

4

Results

5

Discussion

6

Conclusion