Kalman Filtering and LBBG Movement Trajectories

(1)

Kalman Filtering and LBBG Movement

Trajectories

Merel Vogel

June 28, 2020

Bachelor thesis Mathematics

Supervisor: dr. ir. Sonja Cox, dr. ir. Emiel van Loon

Korteweg-de Vries Institute for Mathematics Faculty of Sciences

(2)

Abstract

In this thesis, we study the Kalman filter, both in a theoretical and practical setting. This is done by first establishing sufficient background information on covariance matrices and defining the notion of expectation and conditional expectation in the multivariate scenario, as well as introducing multivariate Gaussians. We subsequently present a rigorous mathematical derivation of the Kalman filter by following and clarifying the proof of Kuman & Varaiya (1986) [3]. Lastly, we analyse given GPS-data from June 2010 describing the movement of eight lesser black-backed gulls (Larus fuscus) by applying the Kalman filter to a proposed state-space model that describes the trajectories of the gulls.

Title: Kalman Filtering and LBBG Movement Trajectories Authors: Merel Vogel, merelmvogel@gmail.com, 11054085 Supervisor: dr. ir. Sonja Cox, dr. ir. Emiel van Loon, Second grader: prof. dr. Sindo Nunez Queija,

End date: June 28, 2020

Korteweg-de Vries Institute for Mathematics University of Amsterdam

Science Park 904, 1098 XH Amsterdam http://www.kdvi.uva.nl

(3)

1. Introduction

It was November 1958 when Rudolf Emil Kalman devised the idea of his now infamous Kalman filter, which would later have an enormous impact in the field of control theory [19]. In engineering and mathematics, control theory attempts to analyse and optimise the behaviour of dynamical systems. The Kalman filter is an algorithm that estimates states of a stochastic system based on indirect and uncertain measurements. The al-gorithm aims to reduce uncertainty in the state as well as in the measurements, and is optimising the system in this way. Kalman was inspired by the optimal estimation methods developed by fellow mathematicians Norbert Wiener and Andrei Kolmogorov in the 1940s, whose ideas he translated into a state-space setting [19]. The filter made its first appearance in journals in 1960 [20]. Many notable applications of the filter have been made ever since, but one worth mentioning occurred during the end of the 1960s. This was when humankind flew a spacecraft to the moon and back with the Apollo pro-gramme, for which the filter was used in the estimation of the trajectory of the shuttle [19]. In a sense, the astronauts’ lives depended on it. The significance of the Kalman filter has been large in many other fields of science as well. To name but a few, these include studies such as carbon cycle data assimilation [21], radar tracking [19], inertial navigation [19] and the estimation of traffic states [22]. It will come as no surprise, that we will be looking at yet another new application of the Kalman filter in this study.

In this thesis, rather than rockets the trajectory of another flying entity is studied, the lesser black-backed gull (Larus fuscus, as seen in figure 1.1 below). A method to implement Kalman filters in the study of the gulls’ movement trajectories is described. GPS-data from several gulls (n = 8) living in and around the north of the Netherlands and Texel is used. The data was retrieved by attaching a ‘GPS-backpack’ to the birds, and the measurements were taken during the entire month of June in 2010. The GPS-device that tracked the gulls is able to measure the birds’ position in longitudinal and latitudinal coordinates. However, the measurements do not reveal the speed of the bird. Besides that, it should be noted that even though GPS-measurements are generally very accurate (and ironically even use Kalman filters [19]), some uncertainty about their exact location is still present – especially since the measurements are taken in discrete-time. Also, there is a certain random perturbation that affects the birds’ trajectories, which ecologists cannot quite explain. The role of the Kalman filter will be to account for this uncertainty, and furthermore to predict an unobserved variable (the gulls’ speed).

(5)

Figure 1.1.: A lesser black-backed gull (Vogelbescherming (2020) [26])

The aim of the study is divided into three parts. Firstly, to derive the Kalman filter and understand it in depth. Secondly, we wish to build a model that describes the trajectories of the lesser black-backed gulls. Following our second aim, we also wish to implement the Kalman filter in R and apply it to the established model. The results will be analysed and discussed.

Determining precise trajectories of the lesser black-backed gulls’ travels is important to learn about the gulls’ behavior. For example, the gulls’ travels can illuminate how much of their time is spent on sea versus land, where they find their food and other such aspects of their behavior. In Chapter 4, further motivation is given as to why it is important to learn about their behavioral patterns and movement trajectories.

The outline of this thesis can be summarised as follows. In Chapter 2 we build a mathematical foundation which will pave the way towards the use of Kalman filters for the analysis of bird movement GPS-data. We do this by providing a suitable toolkit of instruments from probability theory, consisting of the notions of expectation and covariance in the multivariate case, as well as a study of multivariate Gaussians and stochastic systems. This information is required for the rigorous derivation of the Kalman filter. The derivation will play the central role in Chapter 3, where we also examine the required setting for applying the filter. The filter appears in practise in Chapter 4: a model is developed that describes the movement of lesser black-backed gulls, to which we apply the Kalman filter. The model is defined by making use of a stochastic acceleration model and incorporating wind data in the model as well. The results will be analysed and ultimately, this thesis will be concluded by determining reasons why in some circumstances the estimations deviate significantly from the measurements.

(6)

2. Mathematical Preliminaries

In order to derive the functioning of the Kalman filter, some mathematical prerequisites are necessary. Therefore, some important and relevant results are presented in this chapter.

Kalman filters estimate states in a so-called stochastic system. To understand what this is, we will start by stating how to extend definitions of expectation and covariance to the multivariate case. Thereafter some background on multivariate Gaussians will be provided. Using the acquired tools, the chapter is then concluded by introducing the aforementioned notion of a stochastic system.

2.1. Covariance matrices and (conditional) expectation

In this subsection, we introduce the definitions of covariance and expectation in the multivariate framework, as well as the notion of conditional expectation. A few lem-mas about the covariance and expectation are stated, which will be relevant for the calculations and derivations leading towards the Kalman filter.

Definition 2.1.1. Let X be a d-dimensional random vector. Then the expectation of X is a vector in Rd denoted by E[X] and is defined by

E[X] = (E[X1], . . . , E[Xd]).

A consequence of this definition is the following.

Lemma 2.1.1. Let X be a d-dimensional random vector, and let A ∈ Rn×d. Then E[AX] = A E[X].

Proof. Note that by definition,

E[AX] = E     d X j=1 aijXj   n i=1  ,

and now by linearity of the expectation and definition 2.1.1, this is equal to

=   d X j=1 aijE[Xj]   n i=1 = A E[X].

(7)

The following corollary of Lemma 2.1.1 is a direct consequence.

Corollary 2.1.1.1. Let X be a d-dimensional random vector, and let A ∈ Rn×d. Then the following holds.

E[AX] = E[AXX>A>] = A E[XX>]A>.

Definition 2.1.2. For any two random vectors X : Ω → Rn and Y : Ω → Rm, the cross-covariance matrix is a matrix in Rm×n defined as:

cov(X, Y ) = E h

(X − E[X])((Y − E[Y ])> i

. Typically, cov(X, Y ) is also denoted by ΣXY.

Definition 2.1.3. Let X ∈ Rd be a random vector. Then the covariance matrix of X is a matrix in Rd×d defined as:

cov(X) := cov(X, X). Typically, cov(X) is also denoted by ΣX.

In this way, the covariance matrix describes the covariance between the entries of a random vector, whereas the cross-variance matrix gives the covariance between entries from two separate input random vectors. The latter can be viewed as the multidimen-sional version of the covariance.

Lemma 2.1.2. For any two random vectors X : Ω → Rn and Y : Ω → Rm and any two matrices A, B ∈ Rm×n it holds that

cov(X, Y ) = cov(Y, X)>, (2.1) cov(AX, BY ) = A cov(X, Y )B>. (2.2) Proof. The first statement follows almost immediately; observe that

[cov(X, Y )]i,j = E[(Xi− E[Xi])(Yj − E[Yj])] = [cov(Y, X)]j,i.

The second statement is a consequence of definition 2.1.3 in combination with Lemma 2.1.1. Also, recall the following auxiliary result from linear algebra: for any two matrices A, B it holds that (AB)> = B>A>. With this information, and also by using Lemma 2.1.1 again, we can now prove (2.2) as follows:

cov(AX, BY ) = E h

(AX − E[AX])((BY − E[BY ])> i

= Eh(AX − A E[X])((BY − B E[Y ])>i = EhA(X − E[X])((B(Y − E[Y ]))>i = E

h

A(X − E[X])((Y − E[Y ])>B> i

= A Eh(X − E[X])((Y − E[Y ])>iB> = A cov(X, Y )B>.

(8)

Lemma 2.1.3. For any two random vectors X : Ω → Rn and Y : Ω → Rn that are independent, it holds that

cov(X + Y ) = cov(X) + cov(Y ).

Proof. By definition 2.1.3 and the second equation in Lemma 2.1.2 it follows that: [cov(X + Y )]i,j = [cov(X)]i,j+ [cov(Y )]i,j+ cov(Xi, Yi) + cov(Xj, Yi). (2.3)

Since we know that X and Y are independent, which means that E[XY ] = E[X] E[Y ], a straightforward calculation yields the following:

cov(Xi, Yj) = E[(Xi− E[Xi])(Y − E[Yi])]

= E[XiYj − E[Yi]Xi− E[Xi]Yi+ E[Xi] E[Yi]]

= E[Xi] E[Yi] − E[Yi] E[Xi] − E[Xi] E[Yi] + E[Xi] E[Yi]

= 0.

Besides independence of X and Y , we have used the linearity of the expectation. By putting this result back into (2.3), we may conclude that cov(X + Y ) = cov(X) + cov(Y ).

We will now introduce conditional expectation in a measure theoretic fashion. That is, rather than conditioning on the value of a random variable, we condition on a given σ-algebra. In definition (b) which follows we relate it to conditioning on a random variable again.

Definition 2.1.4. (a) Let (Ω, F , P) be a probability space and let G ⊂ F be a σ-algebra. For all X ∈ L1(Ω, G, P), let E[X | G] denote the unique element of L1(Ω, F , P) such that for all G ∈ G it holds that

Z G E[X | G] d P = Z G X d P .

Here the notation X ∈ L1(Ω, G, P) means that X is P-integrable, i.e., the integral of |X| is finite with respect to the measure P.

(b) Now, E[X | Y ] is the (almost surely) unique random variable such that for all y ∈ R it holds that

E[X | σ(Y )] = E[X1{Y ≤y}] = E[E[X | Y ]1{Y ≤y}],

and E[X | Y ] = ϕ(Y ) for some function ϕ : R → R.

To see why E[X | G] is in fact unique, one may use Lemma 8.3 in [16, Chapter 8]. Existence of the conditional expectation is guaranteed by Lemma 8.6 in [16, Chapter 8]. Lemma 2.1.4. Let X, Y ∈ L1(Ω, F , P), and let G be a sub-σ-algebra of F. The following properties hold.

(9)

1. If X ≥ 0 almost surely, then E[X | G] ≥ 0;

2. If X ≥ Y almost surely, then E[X | G] ≥ E[Y | G]; 3. E[E[X | G]] = E[X];

4. If a, b ∈ R, then a E[X | G] + b E[Y | G] = E[aX + bY | G].

Proof. The proof will not be given here; instead [16, Chapter 8] can be consulted. It is important to realize that the conditional expectation is in fact a random variable.

2.2. Random vectors and Gaussians

In this subsection we will define a random vector and some of its properties. A random vector is essentially the multidimensional version of a random variable. The notion of Gaussian random vectors is also introduced.

Definition 2.2.1. Let (Ω, F , P) be a probability space, where Ω is the (non-empty) sample space, F is a σ-algebra on Ω, and P is a probability measure on (Ω, F). A random vector or multivariate random variable is a vector X = [X1, . . . , Xd]>: Ω → Rd, where

each Xi (i ∈ {1, . . . , d}) is a random variable itself on the probability space (Ω, F , P).

We now introduce the notion of Gaussians in the multivariate case due to Lapidoth, A. (2009) [14, Chapter 23]. In order to do so, we will first define a standard Gaussian. Definition 2.2.2. An n-dimensional random vector Z = [Z1, . . . , Zn]>is said to be

stan-dard Gaussian if its components Z1, . . . , Znare independent and each Zi, i ∈ {1, . . . , n}

is Gaussian with mean zero and unit variance.

Definition 2.2.3. An n-dimensional random vector X = [X1, . . . , Xk]> is called

multi-variate Gaussian (or normal ) if there exists some m ∈ N, some matrix A ∈ Rn×m, some µ ∈ Rn and an m-dimensional standard Gaussian Z such that

X = AZ + µ.

Lemma 2.2.1. Let X be an n-dimensional multivariate Gaussian and let m ∈ N, A ∈ Rn×mand µ ∈ Rnbe such that X = AZ + µ for some standard m-dimensional Gaussian. Then E[X] = µ and cov(X) = AA>.

Proof. By linearity of the expectation and Lemma 2.1.1, it holds that E[X] = E[AZ + µ] = A E[Z] + E[µ]. Since E[Z] = 0 because it is a standard Gaussian, it follows that E[X] = E[µ]. Furthermore, cov(X) = cov(AZ + µ) = A cov(Z)A> = AA> by Lemma 2.1.2 and the fact that Z is a standard Gaussian.

If the random vector X is multivariate Gaussian, with mean vector µ and covariance matrix Σ, this is commonly denoted by X ∼ N (µ, Σ). We now proceed by proving a lemma to characterise multivariate Gaussians in more depth. In order to derive this lemma, we first need to bring to mind two properties of one-dimensional random variables that are normally distributed.

(10)

Lemma 2.2.2. Let X and Y be random variables on the probability space (Ω, F , P). 1. If X ∼ N (µX, σX2) and Y ∼ N (µY, σY2), and X and Y are independent, then

X + Y ∼ N (µX+ µY, σ_X2 + σ2_Y).

2. If X ∼ N (µ, σ2) and c, d ∈ R are some constants, then cX + d ∼ N (cµ + d, c2σ2). Proof. The proof of (1) is rather large and irrelevant for the rest of this thesis, and is therefore omitted. However, the article by Eisenberg, B. & Sullivan, R. (2008) [17] gives an overview of possible approaches for the proof. The proof of (2) can be formulated as follows. Since X is Gaussian, we may write X = σZ + µ by definition 2.2.4 and Lemma 2.2.1, where σ, µ ∈ R. Now set Y = cX + d, from which it then follows that Y = c(σZ + µ) + d = (cσ)Z + (cµ + d). This means that we may conclude that Y is again Gaussian and Y ∼ N (cµ + d, c2σ2).

Lemma 2.2.3. A random vector X : Ω → Rd is multivariate Gaussian if and only if for every y ∈ Rd the variable X>y is Gaussian.

Proof. We only prove the forward implication. Since X is multivariate Gaussian, we may write X = AZ + µ, so X>y = Z>A>y + µ>y. By Lemma 2.2.2(2), it suffices to verify that ZTv is Gaussian for all v ∈ Rm, but this follows immediately from Lemma 2.2.2(1). For a proof of the reverse implication see e.g. Jacod & Protter (2004) [18, Theorem 16.2].

Lemma 2.2.4. Let X : Ω → Rd be multivariate Gaussian and let A ∈ Rn×d. Then AX ∈ Rn is multivariate Gaussian.

Proof. Let m ∈ N, B ∈ Rn×m, µ ∈ Rn be such that X = BZ + µ for some standard m-dimensional Gaussian Z. Then AX = ABZ + Aµ, so by definition X is multivariate Gaussian.

Lemma 2.2.5. Suppose the random variable X and Y are jointly Gaussian. Then X and Y are independent if and only if E[XY ] = E[X] E[Y ].

Proof. See Subsection 6.8 in Grimmett & Welsh (2014) [12].

Another concept we need to introduce before we can proceed with Chapter 3, is the notion of a conditional density. Let us recall what this means by stating the following definition from Grimmett & Welsh (2014) [12, Subsection 6.6].

Definition 2.2.4. Let X and Y be jointly continuous random variables, with joint density function pX,Y. The conditional density function of Y given X = x is denoted by

fY |X(· | x) and is defined by

p_{Y |X}(y | x) = pX,Y(x, y) pX(x)

, for y ∈ R and for all values of x satisfying pX(x) > 0.

In the case that X and Y of Definition 2.2.4 are both multivariate Gaussians, we can say something about the mean and covariance of the conditional distribution. This is reflected in Lemma 3.2.1

(11)

2.3. Stochastic system

We define a state space model (SSM) of a stochastic system.

Definition 2.3.1. A discrete time stochastic system is a dynamical system of the fol-lowing form:

xk+1= fk(xk, uk, wk), (2.4)

yk= hk(xk, vk). (2.5)

Equation (2.4) describes the state equation. It tells us that we can find the next state xk+1 by putting together information we have about the previous state xk, some input

uk and noise wk. To be more specific, xk∈ Rn is the state of the system at time k ∈ N,

uk ∈ Rm is the input at time k, and wk ∈ Rg is the (unknown) disturbance or noise at

time k, modelled by a random variable. The second equation (2.5) is the measurement equation. It is assumed that the output yk ∈ Rp can be physically measured. The

unknown measurement error is modelled by the random variable vk∈ Rh. The functions

f0, f1, . . . and h1, h2, . . . are assumed to be known. Lastly, x0 is the initial state of the

system. When yk6= xk, we are dealing with partial observation and it is unreasonable to

assume an exact state of x0. Therefore, in the case of partial observation, x0 is assumed

to be a random variable. Otherwise, x0 can be set equal to the measurement.

It is perhaps appropriate to foreshadow what will happen in Chapter 4, where a stochastic system appears in practise. Therefore, we give a rough sketch of the meaning of these equations in the context of the given data that will be studied. A particular stochastic system will be studied to which the Kalman filter will be applied. The Kalman filter is used to generate the conditional density of the state xk, given the measurements

yk.

In the analysis of gull movement that will follow, xkis the state-vector of the modelled

bird. This state consists of the location and speed of the bird. Also, uk will model the

wind, having an obvious impact on a bird’s trajectory. Furthermore, yk will consist

of the measurements performed by the GPS-device. GPS tracking is subject to some measurement error, and it will play the role of the random variable vk. The variable wk

(12)

3. The Kalman Filter

Now that we have established sufficient background information in Chapter 2, the Kalman Filter can be reviewed in detail. The filter will play the leading role in this chapter. To get started, some intuition will be shaped with an example based on the launch of a rocket. Subsequently a rigorous mathematical derivation will follow. The mathematical background presented in Chapter 2 will be frequently addressed.

3.1. Intuition

As mentioned in the introduction, the Kalman filter estimates states of a stochastic system that includes indirect and uncertain measurements. It is the process of finding the best estimate from noisy data, hence the term ‘filter’: the noise is filtered out, so to say. A common example of applying Kalman filters can be found in the navigation systems of cars. For example, GPS signal-loss occurs when one drives through a tunnel, but you still want to know where the car is located. This results in a need for This can be accomplished with the usage of Kalman filters.

We now describe a more explicit example. In the introduction it was noted that one of the first projects that made use of the Kalman filter was the Apollo programme. To stick to the same theme, we paraphrase an example due to MATLAB [23]. The example shows some intuition behind the usage of the Kalman filter and motivation for what it can be used.

Example 3.1. Imagine you want to launch a spacecraft into space. This reveals many difficulties, one of which is that you will need to monitor the temperature of your engines. They may not become too hot or too cold. Monitoring of the temperature is also necessary in order to adjust fuel usage. In this example the fuel will describe the input equation, as seen in Chapter 2. Since the engines are very hot and a sensor would melt if you install it directly in an engine, the true temperature cannot be measured. Therefore, you decide to install a temperature sensor on the outside of the engine. The measurement of the external temperature will define the output equation. The true temperature will be the state you wish to estimate, as it can not be measured directly. In order to do so, a model is constructed. Since we can observe the true external temperature, we would like the estimated external temperature retrieved from the model to be the equal to the true external temperature. If the estimated external temperature converges to the true external temperature, the internal temperature of the engine we wish to know about will accordingly also converge to the truth. That is, given that all noise is Gaussian as in definition 2.2.3, the Kalman filter minimises the mean square error of the estimated parameters [3, Chapter 7, p. 97].

(13)

Hopefully, some intuition and motivation for using the filter have been shaped by these two examples. In the next section, we investigate the filter from a more precise mathematical point of view and present the technicalities of the filter.

3.2. Derivation of the Kalman filter

In this section, we will derive the Kalman filter in a rather rigorous manner. We follow the proof constructed by Kuman & Varaiya (1986) [3, Chapter 7]. Firstly, we will have a look at the stochastic system that we are going to analyse. Consider the following linear discrete-time stochastic system given by the following state and measurement equations: Xk+1 = AkXk+ BkUk+ GkWk; (3.1)

Yk = CkXk+ HkVk. (3.2)

In these equations, Xk, Uk, Yk, Wk, Vkare random vectors of respective dimensions n, m, p, g

and h. The matrices Ak, Bk, Gk, Ckand Hk are of appropriate dimension, which may or

may not depend on time. Both (Wk)∞_k=0 and (Vk)∞_k=0 are an i.i.d.-sequence, for which

it holds that W0 ∼ N (0, Q) for some given Q ∈ Rg×g and V0 ∼ N (0, R) for some

given R ∈ Rk×k. The random variables X0, (Wk)∞k=0 and (Vk)∞k=0 are assumed to be

independent, and normally distributed:

X0 ∼ N (0, Σ0), Wk∼ N (0, Q), Vk∼ N (0, R). (3.3)

After stating definition 2.3.1 it has already explained what the different random vectors describe, but let us shortly repeat that here: Xkis the state of the system at time k ∈ N,

Wk and Vk describe random perturbations of the system, Yk models the measurements

and lastly Uk is the input of the system. In this case, the input vector can be anything

that you do not wish to estimate and is not part of the state of the object that is being studied, but you still want to incorporate it in the model because it has an impact on the system dynamics.

In this case, the conditional distribution of the state Xkgiven X0 and (Wl)kl=0, (Vl)kl=1

is provided by the Kalman filter. The rest of this subsection will concern the derivation of the filter.

Throughout the next part we will assume that Uk= 0 unless stated otherwise.

Equa-tions (3.1) and (3.2) then simplify to:

Xk+1= AkXk+ Gkwk; (3.4)

Yk= CkXk+ Hkvk, . (3.5)

The random variables Xk, Xk+1 and Ykare jointly Gaussian. This means the conditional

densities are also Gaussian:

p_k|k(Xk| Yk) ∼ N (Xk|k, Σk|k),

(14)

Here we used the following notation:

Xk|k = E[Xk | Yk]; (3.6)

Σk|k = E[(Xk− Xk|k)(Xk− Xk|k)>) | Yk]. (3.7)

That is to say, the conditional expectation of Xkconditioned on Ykis Xk|k(see definition

2.1.4), and the conditional covariance is Σ_k|k. Analogously, we define the conditional mean and covariance in the next timestep as follows:

Xk+1|k = E[Xk+1| Yk]; (3.8)

Σ_k+1|k = E[(Xk+1− Xk+1|k)(Xk+1− Xk+1|k)>) | Yk]. (3.9)

Using the introduced assumptions and definitions, we can now state a lemma. The following lemma is required in the derivation of the Kalman filter.

Lemma 3.2.1. Let X and Y be jointly Gaussian random variables, with X ∼ N (X, ΣX),

Y ∼ N (Y , ΣY), with X ∈ Rd, ΣX ∈ Rd×d, Y ∈ Rn, ΣY ∈ Rn×n and use the notation

ΣXY = cov(X, Y ) as in definition 2.1.3. Let ˆX = E[X | Y ], ˜X = X − ˆX. Then

ˆ

X = X + ΣXYΣ−1_Y (Y − Y ).

Moreover ˜X is independent of Y , and thus of ˆX, and ˜X ∼ N (0, Σ_X˜), with

Σ_X˜ = ΣX − ΣXYΣ−1_Y (Y − Y ).

Proof. Define ˆv = X + ΣXYΣ−1_Y (Y − Y ), and ˜v = X − ˆv. By repeatedly using the

linearity of the expectation and Lemma 2.1.1, we then obtain: E[˜v] = E[X − ˆv]

= E[X − X − ΣXYΣ−1Y (Y − Y )]

= E[X] − E[X] − ΣXYΣ−1_Y E[Y − Y ]

= E[X] − E[X] − ΣXYΣ−1_Y (E[Y ] − E[Y ]) = 0.

Likewise, again by using linearty and Lemma 2.1.1, observe that: E[˜v(Y − Y )>] = E[(X − ˆv)(Y − Y )>]

= E[(X − X − ΣXYΣ−1_Y (Y − Y ))(Y − Y )>]

= E[((X − X)(Y − Y )>− Σ_XYΣ−1_Y (Y − Y ))(Y − Y )>] = E[(X − X)(Y − Y )>] − ΣXYΣ−1Y E[(Y − Y ))(Y − Y )

>_].

= ΣXY − ΣXYΣ−1_Y ΣY

= 0.

This shows that ˜v and Y are independent, since they are jointly Gaussian (see Lemma 2.2.5). Hence we see that

ˆ

(15)

This proves the first assertion. Next, we examine Σ_X˜. As a first observation, note

that since ˆX and ˜X are independent, Lemma 2.1.3 can be applied and it follows that cov(X) = cov( ˜X) + cov( ˆX). Now recall two common results from linear algebra: for any two matrices A, B it holds that (AB)> = B>A>, and (A + B)> = A>+ B>. In addition, as seen in Lemma 2.1.2, covariance matrices are symmetric: Σ>_X = ΣX for

any random vector X. It should also be observed that cov(X, Y ) = (cov(Y, X))> and thus ΣXY = Σ>Y X. Lastly, note that E[ ˜X] = 0. Taking all of this knowledge into

consideration, as a final step we can derive the following: Σ_X˜ = E[ ˜X ˜X>]

= E[(X − X + ΣXYΣ−1_Y (Y − Y ))(X − X + ΣXYΣ−1_Y (Y − Y ))>]

= E[(X − X)(X − X)>− (X − X)(ΣXYΣ−1_Y (Y − Y )>)

− (ΣXYΣ−1_Y (Y − Y ))(X − X)>+ (ΣXYΣ−1_Y (Y − Y ))(ΣXYΣ−1_Y (Y − Y ))>]

= ΣX − E[(X − X)(Y − Y )>Σ−1Y ΣY X] − E[(ΣXYΣ−1Y (Y − Y ))(X − X) >_] + E[(ΣXYΣ−1Y (Y − Y ))((Y − Y ) >_Σ−1 Y ΣY X)] = ΣX − ΣXYΣ−1_Y ΣY X− ΣXYΣ−1_Y ΣY X+ ΣXYΣ−1_Y ΣYΣ−1_Y ΣY X = ΣX − ΣXYΣ−1_Y ΣY X,

which is exactly what we were looking for.

Corollary 3.2.1.1. Let X, Y and Z be jointly Gaussian random variables with X ∼ N (X, Σ_X), Y ∼ N (Y , ΣY) and Z ∼ N (Z, ΣZ), with X ∈ Rd, ΣX ∈ Rd×d, Y ∈ Rn, ΣY ∈

Rn×n, Z ∈ Rk and ΣZ ∈ Rk×k. In addition, let Y and Z be independent. Define

ˆ

X = E[X | Y, Z] and ˜X = X − ˆX. Then ˆ

X = X + ΣXYΣ−1_Y (Y − Y ) + ΣXZΣ−1_Z (Z − Z),

Σ_X˜ = ΣX− ΣXYΣ−1_Y Σ>XY − ΣXZΣ−1_Z Σ>XZ.

Proof. Apply lemma 3.2.1 to random vectors X and [Y, Z]>, and note that

ΣY Z =

ΣY 0

0 ΣZ

, because Y and Z are independent.

Using Lemma 3.2.1 and its corollary, we can now derive the Kalman filter. It consists of equations (3.10)–(3.14), as stated in the following theorem.

Theorem 3.2.2 (Kalman filter). Let Xk, Ak, Ck, Gk, Q and R be as described in the

(16)

Here we have used the following notation:

Lk= Σk|k−1Ck>(CkΣk|k−1Ck>+ HkRHk>) −1

, (3.15)

L0= Σ0C0>(C0Σ0C0>+ H0RH0>)−1. (3.16)

Proof. The proof is divided into three steps. Step 1

For the purpose of convenience, define ˜ X_k|k= Xk− Xk|k, (3.17) and ˜ X_k+1|k = ( Xk+1− Xk+1|k for k ≥ 0, E[X0] = 0 for k = 0. (3.18)

From (3.4) and Lemma 2.1.1, we get

E[Xk+1 | Yk] = E[AkXk+ GkWk| Yk]

= AkE[Xk| Yk] + GkE[Wk | Yk],

and since Wk and Yk are independent, the right part becomes 0, and we may also write

this as

X_k+1|k = AkXk|k. (3.19)

Combining (3.17) and (3.19) gives us ˜

Xk+1|k = AkXk+ GkWk− AkXk|k= AkX˜k|k+ GkWk.

Now we can take the covariance on both sides:

Σ_k+1|k = cov( ˜X_k+1|k) = cov(AkX˜k|k+ GkWk) (3.20)

= Akcov( ˜Xk|k)A>k + Gkcov(Wk)G>k (3.21)

(17)

where (3.21) follows from Corollary 2.1.1.1 and because ˜X_k|k and wk are independent,

meaning we can apply Lemma 2.1.3. Also, (3.22) follows from the definition made at the beginning of this section in (3.3) and (3.7).

Equation (3.22) is the result from this step that we will need to complete the proof later on; we now proceed with step 2.

Step 2 Define

Yk|k−1= E[Yk| Yk−1].

Inserting the definition of Yk in (3.5), we get a new expression for Yk|k−1:

Yk|k−1= E[CkXk+ HkVk| Yk]

Again, using Lemma 2.1.1, the fact that Vkand Yk−1 are independent, and that E[Vk] =

˜

Yk|k−1= CkXk+ HkVk− CkXk|k−1

= Ck(Xk− Xk|k−1) + HkVk

= CkX˜k|k−1+ HkVk. (3.24)

The covariance of ˜Yk|k−1, which we shall denote by Σ k|k−1

˜

Y , is equal to:

Σk|k−1_˜

Y = cov(CkX˜k|k−1+ HkVk).

Because ˜Xk|k−1 and Vk are independent, this covariance can also be written as

Σ_Yk|k−1= CkΣk|k−1Ck>+ HkRHk>. (3.25)

The expression retrieved for ˜Y_k|k−1in equation (3.24) also leads to the observation that E[XkY˜k|k−1> ] = E[Xk(CkX˜k|k−1+ HkVk)>], (3.26)

= E[XkX˜k|k−1> C >

k] + E[XkVk>Hk>]. (3.27)

From the definition for ˜Xk+1|k made in Step 1, we get that Xk = Xk|k−1 + ˜Xk|k−1.

(18)

E[XkX˜_k|k−1> ] = Σk|k−1. Moreover, Xk and Vk are independent by definition. This means

that the last term in (3.27) vanishes, and as a last result in this step we obtain

E[XkY˜k|k−1> ] = Σk|k−1Ck>. (3.28)

Step 3

In the previous step we have seen that Yk and (Yk−1, ˜Yk|k−1) are functions of each other.

This means that we may condition on them in stead of on Yk, without changing anything:

Xk|k= E[Xk | Yk] = E[Xk | Yk−1, ˜Yk|k−1].

In addition, by using lemma 3.2.1, we acquire independence of Yk−1 and ˜Yk|k−1. Hence

in accordance with corollary 3.2.1.1 we see that Xk|k = Xk+ ΣXkYk−1Σ −1 Yk−1(Y k−1_{− Y}k−1_{) + Σ} YkY˜k|k−1Σ −1 ˜ Yk|k−1 ( ˜Yk|k−1− ˜Yk|k−1) = E[Xk| Yk−1] + E[XkY˜_k|k−1> ](ΣY_k|k−1)−1Y˜k|k−1.

Similarly, we see that

Σ_k|k = Σ_k|k−1− E[XkY˜k|k−1> ](ΣYk|k−1) −1

E[XkY˜k|k−1> ] >_.

Finally, we implement the values found in (15) and (16), and by using (7) and (8) we get

Xk|k = Xk|k−1+ Σk|k−1Ck>(CkΣk|k−1Ck>+ HkRHk>)−1Y˜k|k−1, (3.29)

Σ_k|k = Σ_k|k−1− Σ_k|k−1(CkΣk|k−1Ck>+ HkRHk>) −1_C

kΣk|k−1. (3.30)

This concludes Step 3.

The proof is now almost an immediate consequence of Steps 1–3. Let us first show how (3.10) follows. By (3.19) and the definition of ˜Yk+1|k, we see that:

˜

Y_k+1|k = Yk+1− Yk+1|k

= Yk+1− Ck+1Xk+1|k

= Yk+1Ck+1AkXk|k.

Substituting in the result from (3.29) and (3.19) again, this gives us (3.10). Then, (3.11) is derived by combining (3.18), (3.29) and (3.16). Subsequently, using (3.15) in (3.30) gives us the result in (3.12). Furthermore, the result in (3.13) is exactly what we derived as a last result in Step 1. Lastly, (3.23) in combination with (3.30) leads us to (3.14). Corollary 3.2.2.1. Assume the setting of Theorem 3.2.2. Then Xk|k and Σk|k, k ∈ N

can be obtained given (Ak)k∈N, (Ck)k∈N, (Lk)k∈N0, (Yk)k∈N0 and R by recursively solving equations (3.10)–(3.14).

(19)

4. Analysis of LBBG Movement Data

In this chapter, a GPS-dataset retrieved from tracking a population of 8 lesser black backed gulls (LBBG) is studied and analysed. The data shows local movements during breeding season. We propose a model for the movement by combining a stochastic acceleration model and also incorporating the wind. We then apply the Kalman filter.

To start off the chapter, some contextual information about the study of bird move-ment is given. Next, the method of the analysis that was performed is motivated and illustrated. The analysis was done in R, using the recipe described for the Kalman filter as in Chapter 3. The chapter is then concluded by presenting the found results.

4.1. Bird movement and trajectories

Animal movement is broadly studied in today’s scientific world. The understanding of migratory behaviour can produce insight into a population’s behavioral dynamics. Even in the oldest of times, back in hunter-gatherer days, knowledge of how and when and where animal populations move around was an area of interest. This is because people had to capture animals as a source of nutrition, and hence needed to know where to find them. Nowadays, however, scientists’ intrinsic motivation for studying animal trajectories lies elsewhere: in an attempt to care for nature and our co-inhabitants of the earth, understanding animal migration is important to take into account when making any kind of decision regarding the management and conservation of the environment [2]. Furthermore, the understanding of animal movement provides a powerful tool to explain behavioral patterns and the ecology in question [15].

To illustrate, the lives of many birds are lost when they get caught up in the strong current caused by wind turbines, as seen in [10], and many other studies. If researchers understand the choices of birds more thoroughly, so they can perhaps predict when they decide to migrate, municipalities could opt to turn off the turbines when they expect the birds to cross them. In such a way a considerable amount of animals could be saved. Another more simple example is the protection of breeding areas once they have been revealed through observational studies.

However, modelling animal movement is a challenging task. State-space models (com-monly denoted in abbreviated form as SSMs) have been largely accepted as a trustworthy framework to study and analyse animal tracking data [9, 15]. Such state-space models are linear stochastic systems, as stated in Chapter 2, definition 2.3.1. The measurement noise as well as the system noise are assumed to have a Gaussian distribution. Therefore, SSMs allow the usage of Kalman filters, as described in Chapter 3.

Kalman filters have been applied to determine animal trajectories widely, particularly in the study of marine animals, such as sea lions and tuna (as seen in [1], [7] and [11]). It

(20)

is an important tool in the field. Few applications have been made to bird trajectories, but one example can be found in [6]. In their study, veeries are tracked in western Canada making use of geolocators, which show relatively large measurement errors. The trajectories were analysed with the help of the Kalman filter. In this study, due to to the higher accuracy of GPS-location, there is a smaller measurement error. This was taken into account in the construction of the model. Further elaboration on the model and data can be found in the next sections.

4.2. Preparing the data

GPS telemetry-data retrieved from lesser black backed gulls tracked with a GPS-device in the time interval 2010-06-01 to 2010-06-30 was used. The tracking took place during breeding season and consequently the birds were not migrating. The available database used contains tables of various observation variables, of which in particular the ‘tracking data’ was analysed. This in turn contains 66146 observations of the following cate-gories: ‘device info serial’, ‘date and time’, ‘latitude’, ‘longitude’, ‘altitude’, ‘tempera-ture’, ‘satellites used’, ‘positiondop’, ‘2d speed’ and ‘direction’. The column ‘device info serial’ states which gull the measurement concerns; the gulls were named 298, 311, 317, 320, 327, 329, 344 and 355. It was in particular the latitudinal, longitudinal and date and time variables that were considered in this analysis.

The analysis was done in R 4.0.0 (R Core Team, 2020) [24]. The data was first arranged in date and time order, according to the ‘date and time’ observations, with the function arrange from the lubridate-package [25]. Using lubridate again, the duration between measurements was determined (in seconds) and added as a variable in the dataframe.

The birds were assumed to be in their colony if the latitudinal coordinate of the measurement was between 53.0070 and 53.0115, and the longitudinal coordinate was between 4.7150 and 4.7200. The data does not only contain measurements showing that the birds were within their colony, but the data also shows that the birds are leaving the colony. This data was of our particular interest, since we are looking to investigate their movement. The measurements spent outside the colony were subdivided into ‘trips’, in order to make it more insightful by chopping the measurements into manageable, more interpretable parts. A ‘trip’ contains data from the moment the bird leaves the colony until it enters the colony again. For each of the eight birds, five of its trips were used in the analysis.

Since wind is also incorporated in the model, this was also added to the dataframes. Data describing wind speed in m/s during the observation period was retrieved from ESA databases by using the R-package RNCEP, as seen in [13]. Data taken from the pressure level at 1000 hPa was used, since this is where the birds are usually located. The wind speeds were then added to the different dataframes of the trips.

To further make the data ready for analysis, the measured longitude-latitude coordi-nates were converted into coordicoordi-nates in the Dutch coordinate system, in such a manner that every variable in the model is now displayed in matching units (metres). The Dutch coordinate system is designed such that the latitudinal coordinate, now in metres, ranges

(21)

from 300 to 600, and the longitudinal coordinate ranges from 0 to 289. The center point (0,0) is Paris; hence, everything is measured in relative distance to Paris.

4.3. The model

State space models rely on the definition of a stochastic system as seen in Definition 2.3.1. We recite the equations here:

Xk+1= fk(Xk, Uk, Wk); (4.1)

Yk= hk(Xk, Vk). (4.2)

This system can also be considered in the following form:

Xk+1 = AkXk+ BkUk+ GkWk; (4.3)

Yk = CkXk+ HkVk. (4.4)

In either case, the first equation describes the evolution of the animal’s position and state xt, which usually consists of the animal’s position and/or speed. The second equation

shows how xtis partially observed (usually not all states can be observed or measured).

To construct the model for LBBG movement accordingly, the matrices and random vectors in these equations need to be defined. The model was built attempting to take into account as much of the contextual information of the birds as possible. The idea on how to model the location of the birds was primarily based on the random accelerations model as described in Chapter 18, subsection 3 of [4]. The model constructed for this study describes the state vector Xt∈ R6 as follows:

X_t>=X1t X2t X˙1t X˙2t X¨1t X¨2t . (4.5)

Here, X1t is the easting location of the bird, X2t is the northing location of the bird,

˙

X1t is the easting air speed in m/s of the bird and ˙X2t is the northing air speed of the

bird. The last two entries of the state vector, ¨X1t and ¨X2t represent the ground speed

of the gull in easting and northing direction, respectively. The ground speed is equal to the sum of the air speed and wind speed. Furthermore, it is assumed that the system is perturbed by random Gaussian noise Wt∼ N (0, Q). In a standard random accelerations

model, the information given so far is all that is required. Yet in our model the impact of the wind on the bird trajectory was also implemented. This was done by adding another vector to the equation. This vector Ut∈ R2 contains the wind speed in northing

direction in the first entry, and the second entry is the wind speed in easting direction. It is denoted as follows:

(22)

Thus altogether, the given state equation is described as follows:          X1t X2t ˙ X1t ˙ X2t ¨ X1t ¨ X2t          =         1 0 ∆ 0 0 0 0 1 0 ∆ 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 1 0 0                  X1,t−1 X2,t−1 ˙ X1,t−1 ˙ X2,t−1 ¨ X1,t−1 ¨ X2,t−1          +         ∆ 0 0 ∆ 0 0 0 0 ∆ 0 0 ∆         U1t U2t +         W1t W2t W3t W4t W5t W6t         , (4.7)

where ∆ is the period of time between consecutive GPS-measurements in seconds. Note that ˙X1,t−1and ˙X2,t−1 represent the air speed of the bird at the previous timestep, and

the GPS measurements are the ground speed of the bird. The model dynamics thus say that the new location Xj,t is the previous location Xj,t−1 plus ∆ times the previous

velocity ˙Xj,t−1, plus the input from the wind Uj,t times ∆ and random noise Wj,t. The

new velocity ˙Xj,t is the previous velocity ˙Xj,t−1 plus random noise Wj,t. To summarise,

the model is a random acceleration model, which also accounts for another variable, namely the wind.

As the next part of the model, the observation equation was defined. The location and ground speed of the bird can be directly observed through the retrieved GPS-data. It is self-evident that the GPS-measurements give the location, and since they are also marked with a timestamp, the duration between measurements can be computed with which we can calculate the ground speed. However, direct observation of their air speed is not possible. Hence the air speed plays the role of the unknown variable and is to be estimated in this analysis. Also, the observations are subject to some measurement error due to the (relatively low) inaccuracy of the GPS-data, which we implement in the equation by adding a vector Vt ∼ N (0, R). The inaccuracy of GPS-data can roughly

attain a maximum of 50m [2]. Everything included, the measurement equation is thus given by:     Y1t Y2t ¨ Y1t ¨ Y2t     =     1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1              X1,t−1 X2,t−1 ˙ X1,t−1 ˙ X2,t−1 ¨ X1,t−1 ¨ X2,t−1          +     V1t V2t V3t V4t     . (4.8)

It was assumed that the initial state is normally distributed and centered, i.e. X0 ∼

N (0, Σ0), to complete the model. This assumption is required to apply the Kalman

filter.

4.4. Kalman filter in this case

As the GPS-measurements come in discrete times (the GPS-tracker sends a signal every few minutes), a certain uncertainty arises in what the true position of the birds may be. As mentioned before, the measurements can have a measurement error of up to 50m as

(23)

well, which contributes even more to the uncertainty. This is where the Kalman filter comes into play. The Kalman filter was implemented to correct for noise and uncertainty in the data, as well as modelling an unobserved variable, namely the air speed of the bird. In order to apply the Kalman filter, it was decided to write own code which matches the theory given in Chapter 3. This was done after reviewing the available software to implement the Kalman filter in [5]. The used code can be found in appendix A.

4.5. Results

In this section the results from the study are presented. To start, five of the fourty trajectories are shown alongside their filtered trajectory, to illustrate the behavior of the filter. We also list the filtered and unfiltered ground speeds, and the estimates of the air speed.

4.5.1. Filtered and unfiltered trajectories

Below, we have listed a selection of five filtered and unfiltered trips to compare. The unfiltered and filtered version are shown in separate plots because they are very similar and the graphs would overlap when shown in the same plot. It can be observed that the trips start and end in the top right, where the colony is located. Further details on the filtered and unfiltered positions are provided in appendix A, where the standard deviations are listed.

(24)

Figure 4.2.: Filtered and unfiltered trajectory for trip 5 of bird 298

(25)

Figure 4.4.: Filtered and unfiltered trajectory for trip 2 of bird 320

(26)

4.5.2. Ground speed

The ground speed of the birds was computed of both the filtered and unfiltered trajec-tories. It was computed entry-wise by the Pythagorean theorem as follows:

ground speed = s x1,t+1− x1,t ∆ 2 + x2,t+1− x2,t ∆ 2 . As before in the previous sections, we have used the following notation:

x1,t: longitudinal position at time t

x2,t: latitudinal position at time t

∆ : duration between measurements at time t + 1 and t

The mean of all of the filtered trips’ ground speed was found to be 2.40 m/s; the mean of the unfiltered trips’ ground speed was 2.42 m/s. Their Pearson correlation was found to be 0.62, as computed with the cor function in R. Below in table 4.1 the results are listed per bird.

bird mean ground speed (unfiltered) mean ground speed (filtered)

298 3.17 3.15 311 3.13 3.10 317 1.96 1.90 320 3.92 4.00 327 2.27 2.57 329 1.48 1.43 344 2.61 2.52 355 1.18 1.13

Table 4.1.: Mean ground speed per bird in m/s, from filtered and unfiltered trajectories

4.5.3. Air speed

The mean of all of the filtered trips’ air speed was estimated to be 4.33 m/s. Below in table 4.2 the results are listed per bird. The mean is given as the mean of all the different estimated speeds of all five trips per birds.

(27)

bird mean air speed (filtered) 298 4.42 311 5.80 317 3.75 320 5.50 327 3.41 329 4.24 344 4.07 355 3.44

(28)

5. Conclusion and Discussion

5.1. Conclusion

We have provided a framework which allows understanding the exact functioning of the Kalman filter. Besides this, we have built a model for lesser black-backed gull trajectories such that the Kalman filter was applicable. With the established model, the Kalman filter was implemented in R in order to analyse the gulls’ trajectories. It can be concluded that the model provided seems to model the trajectories as expected, and the Kalman filter proves to be a suitable tool in the analysis of the gull’s trajectories.

5.1.1. Conclusion on results from Chapter 4

It can observed that the filtered versions very much resemble the unfiltered versions of the trajectories, as expected – in like manner, both when the trajectories are relatively straight (such as in Figure 4.1 and Figure 4.5) and as well as when they are more wiggly (such as in Figure 4.2 and Figure 4.3). It can also be seen that the filtered trajectories are sometimes shorter.

As a first observation, it is noted that in 24 out of 40 (60%) of cases, the filtered location has higher standard deviation, both for northing and easting locations. The highest differences in standard deviation occur in easting location of bird 344’s third and fifth trip, as seen in Table B.7, which shows a difference of 643 and 487m respectively.

Secondly, the mean ground speed of the filtered version was lower for six out of eight birds (75%). This is also reflected in the correlation. This could be due to the fact that the filtered trajectories are shorter (although this has to be verified in future conducted research), and accordingly the birds would have had the same amount of time to travel shorter distance, resulting in a lower ground speed.

5.2. Discussion

We will not be discussing about Chapter 1, 2 and 3 in much detail, since they are part of a literary study and there is not much to be said about that in terms of discussion. However, there are a few things that are important to mention regarding the analysis done in Chapter 4. First of all, this experiment was conducted with fixed matrices Q and R. These are the covariance matrices of the process and measurement error vectors. These matrices have a notable influence on the filtering, in the following way. Implementing a large-valued covariance matrix Q means it is assumed that the state equation is subject to large error, and in essence, choosing a large Q is thus like telling the filter to correct more (and vice versa: if you pick small values for Q, the filter

(29)

will correct less). Analogously, picking large values for the covariance matrix R of the measurement noise expresses the belief that the measurements are not very accurate, and results in the filter correcting less with the measurement update step. Although it is known that GPS-measurements are accurate up to around 50m and hence do not show large errors, it has not quite been established how much random perturbation the state of the gulls is subject to. Therefore, a fixed choice of R probably makes sense, but it would be interesting to experiment with the values of Q in future studies. The fixed choice of Q in this particular experiment has undeniably had an influence on the results of the filtering.

Another decision that was made to set up the experiment is the choice of the covariance matrix Σ0 of the initial state of the gull (its position and velocity). This one is of

less importance, since it does not play a role in the algorithm itself. Is is merely an assumption used during the first filter step and the filter adjusts well to it.

A note on data pre-processing is perhaps also in place. GPS-measurements are gen-erally less accurate when the speed of the tracking object is below 2 m/s, and this was not taken into account while selecting trips from the data of the birds. This is a feasible source of error in the obtained trajectories, and could be incorporated in future studies to reduce error even more.

Furthermore, other than the few conclusions that we have drawn from the output of applying the Kalman filter here, much more can possibly be said or investigated about the results. We have now only compared ground speed. However, other variables such as air speed can also be analysed. In addition, one could look at the length of the different filtered and unfiltered trajectories. It could make sense if the filtered trajectories are less long, since the filter performs a kind of ‘regression’ (it is perhaps unnecessary to mention, but the filter is not really doing regression – it is making a time-step prediction). The different lengths of filtered and unfiltered trajectories can be seen in the exemplary figures 4.1–4.5 and could be investigated further. Differences in length could be quantified in an exact manner, rather than concluding they are shorter by only looking at the trajectory graph.

Another question that could be addressed is the following: how do the gulls adjust their speed? It could be looked into where big changes in velocity appear and if this can be linked to certain circumstances.

Moreover, the 40 selected trips of the birds could be subdivided into an outgoing and returning part, and see if there is anything to be said about that. For example, it could be possible that when the gull is flying in direction home, it would prefer a shorter route over a longer one, whereas if the gull is still outgoing it is perhaps on a search for food and is not on the way to a particular destination. But these are mere speculations and could be researched further in other studies.

(30)

Bibliography

[1] Sibert, J. R., Musyl, M. K., & Brill, R. W. (2003). Horizontal movements of bigeye tuna (Thunnus obesus) near Hawaii determined by Kalman filter analysis of archival tagging data. Fisheries Oceanography, 12(3), 141–151. https://doi.org/10.1046/j. 1365-2419.2003.00228.x

[2] Hooten, M. B., Johnson, D. S., McClintock, B. T. & Morales, J. M. (2017). An-imal movement: Statistical models for telemetry data. https://doi.org/10.1201/ 9781315117744

[3] Kumar, P.R. & Varaiya, P. (1987). Stochastic Systems: Estimation, Identifica-tion and Adaptive Control. Prentice-Hall, Inc. (Vol. 24). https://doi.org/10.1177/ 002072098702400424

[4] Borovcnik, M., Bentz, H.J. & Kapadia, R. (1991). Machine Learning: A Probabilistic Perspective. Chance Encounters: Probability in Education. https://doi.org/10. 1007/978-94-011-3532-02

[5] Tussel, F. (2011). Kalman Filtering in R. Journal of Statistical Software, 39(2). [6] Hobson, K. A. & Kardynal, K. J. (2015). Western Veeries use an eastern

shortest-distance pathway: New insights to migration routes and phenology using light-level geolocators. The Auk: Ornithological Advances, 132(3), 540–550. https://doi.org/ 10.1642/auk-14-260.1

[7] Royer, F. & Lutcavage, M. (2008). Filtering and interpreting location errors in satel-lite telemetry of marine animals. Journal of Experimental Marine Biology and Ecology, 359(1), 1–10. https://doi.org/10.1016/j.jembe.2008.01.026

[8] Fleming, C. H., Sheldon, D., Gurarie, E., Fagan, W. F., LaPoint, S. & Calabrese, J. M. (2017). K´alm´an filters for continuous-time movement models. Ecological Infor-matics, 40(April), 8–21. https://doi.org/10.1016/j.ecoinf.2017.04.008

[9] Breed, G. A., Costa, D. P., Jonsen, I. D., Robinson, P. W. & Mills-Flemming, J. (2012). State-space methods for more completely capturing behavioral dynamics from animal tracks. Ecological Modelling, 235–236, 49–58. https://doi.org/10.1016/j. ecolmodel.2012.03.021

[10] Barrios, L. & Rodr´ıguez, A. (2004). Behavioural and environmental correlates of soaring-bird mortality at on-shore wind turbines. Journal of Applied Ecology, 41(1), 72–81. https://doi.org/10.1111/j.1365-2664.2004.00876.x

(31)

[11] Heerah, K., Andrews-Goff, V., Williams, G., Sultan, E., Hindell, M., Patterson, T. & Charrassin, J. B. (2013). Ecology of Weddell seals during winter: Influence of environmental parameters on their foraging behaviour. Deep-Sea Research Part II: Topical Studies in Oceanography, 88–89, 23–33. https://doi.org/10.1016/j.dsr2. 2012.08.025

[12] Grimmett, G. & Welsh, D. (2014). Probability: An Introduction. Oxford University Press, Second Edition. ISBN 978-0-19-870997-8.

[13] Kemp, M. U., van Loon, E.E., Shamoun-Baranes, J. & Bouten, W. (2012). RNCEP: Global weather and climate data at your fingertips. Methods in Ecology and Evolution, 3(1), 65–70. https://doi.org/10.1111/j.2041-210X.2011.00138.x

[14] Lapidoth, A. (2009). A Foundation in Digital Communication. Cambridge Univer-sity Press. ISBN 978-0-521-19395-5.

[15] Patterson, T. A., Thomas, L., Wilcox, C., Ovaskainen, O. & Matthiopoulos, J. (2008). State-space models of individual animal movement. Trends in Ecology and Evolution, 23(2), 87–94. https://doi.org/10.1016/j.tree.2007.10.009

[16] Spreij, P.J.C. & Cox, S.G. (2019). Measure Theoretic Probability, lecture notes from course Measure Theoretic Probability at the University of Amsterdam.

[17] Eisenberg, B. & Sullivan, R. (2008). Why Is the Sum of Independent Normal Random Variables Normal? Mathematics Magazine. Retrieved 12-05-2020. https: //www.maa.org/sites/default/files/Eisenberg12-0825740.pdf

[18] Jacod, J. & Protter, P. (2004). Probability Essentials. Springer, Second Edition. ISBN 978-3-540-43871-7.

[19] Grewal, M. S. & Andrews, A. P. (2010). Applications of Kalman Filtering in Aerospace 1960 to the Present [Historical Perspectives]. IEEE Control Systems Mag-azine, vol. 30, no. 3, pp. 69-78. https://doi.org/10.1109/MCS.2010.936465. [20] Kalman, R. E. (1960). A new approach to linear filtering and prediction problems.

ASME Trans., J. Basic Eng., ser. D., vol. 82, pp. 35–45.

[21] Kang, J.S. et al (2011). Variable localization in an ensemble Kalman filter: Ap-plication to the carbon cycle data assimilation. J. Geophys. Res., 116, D09110. https://doi.org/10.1029/2010JD014673.

[22] Tampere, C. M. J. & Immers, L. H. (2007). An Extended Kalman Filter Application for Traffic State Estimation Using CTM with Implicit Mode Switching and Dynamic Parameters. IEEE Intelligent Transportation Systems Conference, Seattle, WA, pp. 209-216. https://doi.org/10.1109/ITSC.2007.4357755.

[23] MATLAB (2017). Understanding Kalman Filters, Part 1: Why Use Kalman Fil-ters? Retrieved 06-04-2020. https://www.youtube.com/watch?v=mwn8xhgNpFY

(32)

[24] R Core Team (2020). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. http://www.R-project. org/

[25] Grolemund, G. & Wickham, H. (2011). Dates and Times Made Easy with lubridate. Journal of Statistical Software, 40(3), 1–25. http://www.jstatsoft.org/v40/i03/. [26] Vogelbescherming (n.d.). Kleine mantelmeeuw. Retrieved 06-06-2020.

https://www.vogelbescherming.nl/ontdek-vogels/kennis-over-vogels/ vogelgids/vogel/kleine-mantelmeeuw.

(33)

Populaire samenvatting

Het bestuderen van trajecten van dieren is een belangrijk onderzoeksgebied binnen de ecologie. Migreer- en bewegingsdata kunnen ons veel leren over het gedrag van de onderzochte dieren. In deze scriptie wordt onder andere een manier om het gedrag van de kleine mantelmeeuw te onderzoeken beschreven aan de hand van bewegingstrajecten en het toepassen van de Kalman filter.

We hebben GPS-data gebruikt van acht kleine mantelmeeuwen, die allemaal gedurende de maand juni 2010 een soort ‘GPS-rugzakje’ hebben gedragen. Het GPS-apparaat gaf ongeveer elke 10 minuten een seintje door waar de vogel is geweest, in longitudinale en longitudinale co¨ordinaten. Dit seintje kunnen we zien als een vector [x1, x2]>.

Gebruikmakend van de gegeven GPS-data, kunnen we het formuleren van een model voor de beweging van de kleine mantelmeeuwtrajecten als volgt aanpakken. Het eerste stukje informatie dat we hebben is de GPS-meting. We beschouwen eerst de meting van de longitudinale positie, laten we die gemeten positie x1 noemen. We willen nu met ons

model een voorspelling doen over wat de volgende positie van de vogel zal zijn. Het is dan vrij logisch om te bedenken dat de volgende positie wordt gegeven door een bepaalde verplaatsing bij de vorige positie op te tellen. Die verplaatsing wordt bepaald door de snelheid van de vogel (in m/s) te vermenigvuldigen met een bepaalde tijdseenheid (in s), om op die manier weer iets in meters te verkrijgen. Die bepaalde tijdseenheid in dit geval is de tijd die het heeft geduurd voor het GPS-apparaat de volgende meting doorgaf. Dat is als het ware de hoeveelheid tijd die de vogel heeft gehad om zich te verplaatsen. Laten we die hoeveelheid tijd aanduiden met ∆, en de snelheid met v1. We

kunnen dan zeggen dat de volgende (longitudinale) positie van de vogel, die we voor nu aanduiden met x0₁, wordt gegeven door de volgende vergelijking:

x0₁= x1+ ∆v1.

Het op deze manier beschrijven van de dynamica heet een stochastisch acceleratie model. Dezelfde redenering kunnen we maken voor de latitudinale positie van de vogel.

Omdat de GPS-metingen onderhevig zijn aan een bepaalde fout, is het wenselijk om hiervoor te corrigeren in het model dat we gaan gebruiken. Hiervoor is een handige methode ontwikkeld die de Kalman filter heet. Het algoritme is niet zo ´e´en, twee drie uit te leggen; in deze scriptie is dan ook een heel hoofdstuk gewijd aan afleiden van de Kalman filter. We kunnen echter wel het idee van de filter schetsen door middel van een analogie met een koffiefilter. Om koffie te zetten, leg je gemalen koffie over de filter en giet je er heet water overheen. De gemalen koffie en heet water kunnen we zien als de ‘input’. Vervolgens filtert de koffiefilter dit en is het eindresultaat een kop koffie. Dit kunnen we zien als de ‘output’. Als residu houden we wat koffiedrap over. De Kalman filter werkt aan de hand van eenzelfde soort mechanisme. De input die we nu gebruiken

(34)

zijn de GPS-metingen (die dus bepaalde fouten bevatten), en bepaalde aannames die we doen over de positie en snelheid van de vogel. Als we dit nu door de Kalman filter laten lopen, krijgen we als output (onze kop koffie) de beoogde posities en snelheden, en als residu (de koffiedrap) hebben we de foutbevattende metingen. Dergelijke systemen met in- en output variabelen worden bestudeerd in het vakgebied dat control theory heet.

De Kalman filter hebben we toen gebruikt om de bewegingsdynamica van de vo-gel te beschrijven. Kalman filters toepassen op diertrajecten is al vaker gedaan, maar nog niet veel in studies die vogels betreffen – het ging dan vooral om dieren als tonijn en zeeleeuwen. Met behulp van de filter is het ook mogelijk om een ongeobserveerde variabele te voorspellen. De ongeobserveerde variabele in dit experiment was de lucht-snelheid van de vogel.

(35)

A. R-code

Below, trip329 1 denotes that we are using data from the first trip of bird 329. #load data

trip329_1 <- readRDS("trip329_1")

#initialize required vectors

#tot is the total amount of measurements to determine length of for-loop tot <- length(trip329_1$X)

#delta is the duration between measurements (in seconds) delta <- trip329_1$duration

l_x <- trip329_1$X[1] #location in easting, m. inital position is the first observation l_y <- trip329_1$Y[1] #location in northing, m

dxt <- numeric(1) #speed in easting, m/s dyt <- numeric(1) #speed in northing, m/s

#intitialize the output values so we can save them later Zx <- numeric(1) Zy <- numeric(1) Px <- numeric(1) Py <- numeric(1) Pdx <- numeric(1) Pdy <- numeric(1) Rdx <- numeric(1) Rdy <- numeric(1) Kx <- numeric(1) Ky <- numeric(1) Kdx <- numeric(1) Kdy <- numeric(1) out <- data.frame(l_x,l_y,dxt,dyt,Zx,Zy,Px,Py,Pdx,Pdy,Rdx,Rdy,Kx,Ky,Kdx,Kdy)

#append initial values

x=c(l_x[1], l_y[1], 0, 0) #state vector

u=c(trip329_1$uwind[1], trip329_1$vwind[1]) #input vector, the wind in m/s

#initialize required matrices

(36)

nrow = 4, ncol = 4)

B = matrix(c(delta[1], 0, 0, 0, 0, delta[1], 0, 0), nrow = 4, ncol = 2) H = matrix(c(1,0,0,1,0,0,0,0), nrow=2, ncol=4)

P = matrix(c(rnorm(1), 0, 0, 0, 0, rnorm(1), 0, 0, 0, 0, rnorm(1), 0, 0, 0, 0, rnorm(1)), nrow=4, ncol=4) G = t(c(0.5*delta[1]^2, 0.5*delta[1]^2, delta[1], delta[1])) Q = G%*%t(G)*0.3

R = matrix(c(rnorm(1), 0, 0, rnorm(1)),nrow=2,ncol=2) #measurement noise I=diag(4)

#the for loop which describes the Kalman Filter for (i in 1:tot){

x=A%*%x+B%*%u P=A%*%P%*%t(A)

#define the Kalman gain S=H%*%P%*%t(H)+R

K=P%*%t(H)%*%pseudoinverse(S)

#update the state estimate

Z=c(trip329_1$X[i],trip329_1$Y[i]) y=Z-(H%*%x)

x=x+K%*%y

#update the error covariance matrix P=(I-K%*%H)%*%P

#update matrices with new delta and new random values from normal distrib

A = matrix(c(1, 0, 0, 0, 0, 1, 0, 0, delta[i+1], 0, 1, 0, 0, delta[i+1], 0, 1), nrow = 4, ncol = 4)

B = matrix(c(delta[i+1], 0, 0, 0, 0, delta[i+1], 0, 0), nrow = 4, ncol = 2) H = matrix(c(1,0,0,1,0,0,0,0), nrow=2, ncol=4)

P = matrix(c(rnorm(1), 0, 0, 0, 0, rnorm(1), 0, 0, 0, 0, rnorm(1),

0, 0, 0, 0, rnorm(1)), nrow=4, ncol=4)

G = t(c(0.5*delta[i+1]^2, 0.5*delta[i+1]^2, delta[i+1], delta[i+1])) Q = G%*%t(G)*0.3

R = matrix(c(rnorm(1), 0, 0, rnorm(1)),nrow=2,ncol=2) u = c(trip329_1$uwind[i+1], trip329_1$vwind[i+1])

#save the acquired states in each loop so we can plot them

new <- c(x[1],x[2],x[3],x[4],Z[1],Z[2],P[1,1],P[2,2],P[3,3],P[4,4], R[1,1],R[2,2],K[1,1],K[2,1],K[3,1],K[4,1])

out <- rbind(out,new) }

(37)

B. Standard deviation of unfiltered and

filtered positions

In the following 8 tables, standard deviations of filtered and unfiltered locations are shown, and the difference between them. The columm heading ‘trip’ denotes the trip of the bird, StDev x1t denotes the standard deviation of unfiltered easting location, StDev

ˆ

x1tdenotes the standard deviation of filtered easting location, ‘difference’ denotes StDev

x1t minus StDev ˆx1t, StDev x2t denotes the standard deviation of unfiltered northing

location, StDev ˆx2t denotes the standard deviation of the filtered northing location, and

finally ‘difference’ denotes StDev x2t minus StDev ˆx2t.

trip StDev x1t StDev ˆx1t Difference StDev x2t StDev ˆx2t Difference

1 5523.61 5447.70 75.91 5897.19 5767.60 129.59

2 9542.43 9650.28 -107.85 8217.70 8407.34 -189.65

3 7197.73 7267.84 -70.10 13176.14 13336.0 -159.895

4 10274.43 10151.25 123.18 11544.57 11535.46 9.11

5 9828.77 9971.35 -142.58 10607.97 10659.82 -51.86

Table B.1.: Standard deviation of filtered and unfiltered location for different trips of bird 298

1 9429.95 9564.98 -135.03 5764.40 5874.77 -110.37

2 5392.91 5312.29 80.62 12791.71 12905.89 -114.18

3 6627.39 6613.73 13.65 16806.83 16660.34 146.48

4 3035.42 3118.07 -82.66 11771.54 12157.94 -386.40

5 5900.29 6167.35 -267.06 7683.07 7779.02 -95.94

(38)

trip StDev x1t StDev ˆx1t Difference StDev x2t StDev ˆx2t Difference 1 13357.60 13358.89 -1.29 6090.90 6086.69 4.22 2 4907.47 5230.72 -323.25 699.11 738.01 -38.90 3 12939.43 12921.52 17.91 7055.53 7101.69 -46.16 4 5840.74 5844.42 -3.68 14661.12 14696.94 -35.81 5 6545.15 6937.20 -392.06 3887.67 3850.94 36.73

1 11636.07 11806.54 -170.46 8075.05 8173.02 -97.98

2 8042.44 8346.87 -304.43 5664.41 5880.30 -215.90

3 7468.15 7483.00 -14.85 11495.57 11680.55 -184.98

4 3202.64 3126.64 76.00 4165.71 4075.41 90.30

5 5835.98 6027.17 -191.19 3641.47 3739.95 -98.50

1 7016.27 6931.69 84.59 4903.95 4604.14 299.80

2 2624.16 2622.31 1.85 728.83 733.39 -4.56

3 1883.60 1941.09 -57.50 1199.73 1157.80 41.92

4 3063.86 2995.61 68.25 2343.33 2458.94 -115.60

5 1285.15 1268.36 16.79 3767.48 3819.39 -51.91

1 4288.18 4349.88 -61.69 3013.63 3010.68 2.94

2 5457.02 5586.56 -129.54 3728.52 3914.87 -186.35

3 3352.98 3330.31 22.67 1780.18 1788.83 -8.66

4 5714.39 5716.97 -2.58 3952.22 3907.19 45.03

5 8357.50 8333.37 24.13 9017.75 8975.30 42.45

(39)

trip StDev x1t StDev ˆx1t Difference StDev x2t StDev ˆx2t Difference 1 12245.98 12222.50 23.49 2090.93 2078.74 12.18 2 6791.67 6707.38 84.29 7403.21 7373.30 29.91 3 4687.08 4951.79 -264.71 13996.74 14639.69 -642.95 4 6738.26 6989.74 -251.48 8047.92 8235.71 -187.79 5 4975.85 5114.73 -138.88 12619.78 13106.41 -486.62 Table B.7.: Standard deviation of filtered and unfiltered location for different trips of

bird 344

1 3499.19 3534.33 -35.15 3183.67 3129.26 54.41

2 4407.85 4443.19 -35.34 3026.46 3043.40 -16.94

3 3416.17 3382.47 33.69 2900.41 2798.24 102.16

4 3863.90 3880.47 -16.57 3165.06 3221.05 -55.99

5 1204.58 1340.56 -135.98 1444.84 1584.15 -139.30

Kalman Filtering and LBBG Movement Trajectories