finally for treating me as family throughout these past years. I would also like to thank my supervisor at the Lebanese University, Mr. Imad Mougharbel, for sharing his ideas and experience, for giving me the freedom to explore my own ideas and for

LIST OF FIGURES

INTRODUCTION

Theme and objectives

A wireless sensor network (WSN) consists of a number of spatially distributed sensor nodes which have capabilities of monitoring a physical condition and of cooperating by passing data through the network. Initially motivated by military applications (e.g. surveillance of a battlefield), the properties associated with WSNs, such as their autonomy, low cost, small size, diversity and effectiveness in hardly accessible regions, expanded their use to consumer and industrial applications. However, the advantages of WSNs are limited by a number of critical issues, the most popular of which relies probably in the inherent restrictions on bandwidth and energy consumptions. Typically, ubiquitous sensing might result in overwhelming volumes of raw data. In order to extract the small amount of useful information from the network without communicating large volumes of possibly irrelevant raw data, local processing of the information at the sensor nodes is often useful, this leads to the need of distributed algorithms when performing estimation and/or inference.

Another important difficulty associated with WSNs applications is caused by the fact that data gathered by a sensor node is usually noise-corrupted and requires special caution when performing tasks such as decision making or inference. In fact, many problems have inherent uncertainties, a famous example is that of self-localisation of a mobile robot which receives noisy measurements from its wheels encoders and gets inaccurate views from its range finders. Uncertainties are typically modelled as statistical and thus probabilistic methods offer an attractive framework for addressing the problem. In these approaches, attributes of interest are represented using collections of random variables. Often probability distributions associated to these random variables are complex and in the presence of non-linearities, sequential Monte Carlo methods offer an attractive solution. However, these techniques commonly require high computational complexity due to the large number of samples (or particles) needed to represent complex probability density functions (PDFs). On the other hand, interval estimation is another promising methodology which allows complex PDFs to be estimated using fewer box-particles and is also independent of non-linearities.

In this thesis, we are interested in deriving fusion algorithms for data collected from a WSN while respecting energy and bandwidth restrictions, and taking into account the noisy aspect of the sensed data. Bayesian inference forms the central theme of this manuscript since it acts as the hub to which various analytical and computational procedures attach in order to solve a variety of problems. Mainly, our work can be structured into two parts.

The first part deals with the problem of inference in distributed systems, such as sensor networks, where the probability model is stored in a distributed fashion. In these applications, each node (e.g. sensor) receives information about its local environment and has a local set of attributes (described in terms of random variables) for which it needs to compute a posterior density. Thus each node only stores a relevant portion of the model and has to collaborate with other nodes in the network in order to compute its marginal posterior of interest, given all the observations available to the system [START_REF] Paskin | Exploiting Locality in Probabilistic Inference[END_REF].

An attractive formalism allowing to exploit conditional independence properties of a joint PDF associated with large scale systems is possible through the use of probabilistic graphical models [START_REF] Koller | Probabilistic Graphical Models, Principles and Techniques[END_REF]. Built upon the graphical model, the message passing (or Belief Propagation) approach relies on exchange of messages between neighbour graph nodes. These messages (characterising local/marginal PDFs) are iteratively exchanged between neighbour nodes until convergence is reached. Each node is then able to estimate a specific marginal of the joint posterior density. This approach has many advantages: it is general and can be applied to any network without necessitating knowledge of the overall graphical structure (since each node only knows about its immediate neighbours), and is a distributed algorithm by construction.

Belief Propagation (BP) can be implemented in a parametrised manner (e.g., based on Gaussian PDFs) or in a non-parametric way (e.g., based on sample representation of PDFs). Whilst in the parametric implementation, parameters of the PDF of interest are estimated, e.g., the mean and variance of the Gaussian PDF, in the non-parametric implementation, the PDF is unknown and represented using samples. The Non-parametric BP (NBP) [START_REF] Sudderth | Nonparametric belief propagation[END_REF] is a particle-based implementation of the BP algorithm [START_REF] Sudderth | Nonparametric belief propagation[END_REF] and is especially useful for representing joint and marginal probability distributions that are multi-modal, non-standard Gaussian and in the presence of non-linearities.

The second part of the thesis considers the problem of detecting and localising several sources emitting a chemical compound in the air using exclusively concentration measurements provided by a network of wireless sensors assumed capable of sensing the concentration of the emitted compound at their positions. Various applications are associated with this scenario especially in environmental monitoring and security fields (e.g. in the case of accidental or intentional release of a toxic agent) [START_REF] Weimer | An approach to leak detection using wireless sensor networks at carbon sequestration sites[END_REF][START_REF] Lushi | An inverse Gaussian plume approach for estimating atmospheric pollutant emissions from multiple point sources[END_REF]. Another application of interest in this thesis is that of land mines detection and localisation [START_REF] Jeremić | Landmine detection and localization using chemical sensor array processing[END_REF].

The localisation problem is formulated as Bayesian inference and a Markov Chain Monte Carlo (MCMC) sampler allows to estimate the properties of the posterior PDF of interest. The performance of the probabilistic scheme is also compared to that of an optimisation least squares approach for two scenarios typically associated with land mine detection and localisation [START_REF] Jeremić | Landmine detection and localization using chemical sensor array processing[END_REF]: the first, considers the mines as instantaneous point sources (it assumes that the evaporation of the explosive chemical is induced), the second considers the mines as continuous release point sources that were buried long enough for a steady state concentration profile to be established.

Contributions Part I

The efficiency of the NBP depends on the number of particles, and typically the imprecision on the available information requires a large number of samples (or particles) which leads to high computational complexity. NBP's operations turn out to be complex and computationally demanding. To deal with the high number of particles, interval-based estimation has arisen as a promising methodology and has been widely employed in state estimation problems [START_REF] Gning | Constraints propagation techniques on intervals for a guaranteed localization using redundant data[END_REF][START_REF] Jaulin | Nonlinear bounded-error state estimation of continuous-time systems[END_REF]. The bounded error framework has proven to be robust to non-linearities and able to estimate complex PDF using a small number of box-particles.

We propose a novel non-parametric Belief Propagation approach, which we call box Belief Propagation (box-BP). In the box-PP, a collection of box particles [START_REF] Gning | Integrated Tracking, Classification, and Sensor Management: Theory and Applications, chapter Particle Filtering Combined with Interval Methods for Tracking Applications[END_REF] represents the information exchange between nodes of the network. Adopting an approximation of the Bayesian solution with mixtures of uniform PDFs, having boxes as support [START_REF] Gning | An introduction to box particle filtering[END_REF], we provide the theoretical Bayesian derivation of the box-BP and the corresponding implementation.

Furthermore, the performance of the novel box-BP is validated on localisation applications described by static and dynamic graphical models. The first version of the box-BP [START_REF] Haj Chhadé | Non parametric distributed inference in sensor networks using box particles messages[END_REF] is focused on a sensor calibration problem, i.e., finding the coordinate of all static sensors in the network knowing a prior position of a few of them, for this scenario the efficiency of the box-BP is demonstrated over simulated data. The second application considers the localisation of a single mobile target moving within a network of wireless range sensors, it is more general and proposes a solution to the more difficult, dynamic cases, with dynamic posterior PDFs (e.g., for sequentially estimating the state of a moving object). Additionally, results are reported over simulated and real data.

In both cases, the novel box-BP is compared with the standard particle-based Nonparametric Belief Propagation (NBP). As the NBP, the box-BP is aimed for cases of non-linearities, multi-modal PDFs and, compared with the NBP, has several key advantages. First, it is shown through the Bayesian derivation that it requires simpler and faster computations: this is essentially due to the fact that, in comparison with NBP's representation with samples, it uses box-particles which are interpreted as an approximation of PDFs based on mixtures of uniform distributions. Secondly, it needs much less box particles to estimate the PDFs which is doubly advantageous in terms of memory optimisation and reduction of energy and bandwidth requirements for the messages sent in the network. Thirdly, the box-BP has properties inherent to bounded errors approaches, e.g. it can deal with quantised measurements and bounded noises with unknown statistics, and is able to propagate, in time, supports that are guaranteed to contain the real solutions [START_REF] Gning | Constraints propagation techniques on intervals for a guaranteed localization using redundant data[END_REF].

Hence, the contribution of part I is threefold: first, a novel box-Belief Propagation (box-BP) approach is proposed where the BP is performed based on box representations of the state PDFs; second, the effectiveness of the box-BP approach is demonstrated over simulated data in the case of static graphical models; third, the advantages of the box-BP are proven, using real dataset, for a more general case including a dynamic structure of the graph.

Part II

All existing methods for land mine detection consider a known number of sources in a specific region [START_REF] Badine | Considerations for implementing a wireless network for landmine detection[END_REF][START_REF] Jeremić | Landmine detection and localization using chemical sensor array processing[END_REF]. However this information rarely exists for real applications. One original contribution of this part of the thesis is to present a solution for an unknown number of land mines or vapour-emitting sources. We propose a two-step method where the objective is to estimate the number of sources and then localise them. In the first step, the set of concentration measurements which have been made by the detection system are grouped in a matrix and a PCA scheme is used in order to determine the number of sources. Once this number is known we are then able, in the second step, to locate explosive sources and to estimate the emission rate of each source. Thus, in difference to most of existing methods (which solve the problem of localising a single source), we solve the problem of localising multiple anti-personnel land mines.

Hence, the main novelty in part II consists in considering the case of multiple land mines whose number is unknown a priori.

Thesis outline

As mentioned earlier, the manuscript revolves around two main parts and is organised in six chapters.

Chapter 1 provides an overview of the prior work and the concepts relevant to the thesis. It basically recalls the Monte Carlo sample-based methods as well as an introduction of graphical models and message-passing algorithms.

Chapters 2, 3 and 4 constitute the core of part I and revolve around the three contributions mentioned in the previous section. Chapter 2 presents the main theoretical contribution of this part: the derivation of the box-BP where the BP is performed based on box representations of the state PDFs. The nodes' instantaneous measurement are fused following a graphical model representation and it is shown that the box-BP provides an efficient solution to the Bayesian inference problem. The derivation of the box-BP relies on the interpretation of a box as a uniform PDF defined over its support, and shows how to implement the BP in the bounded error framework. In chapter 3, the performance of the novel box-BP is tested on the problem of calibration in wireless sensor networks. That is the positioning of a number of randomly deployed sensors, according to some reference defined by a set of anchor nodes for which the positions are known a priori. The effectiveness of the new message-passing scheme is demonstrated over simulated data on a static graphical structure. Chapter 4 investigates the problem of continuously localising an object moving within a wireless sensor network. A wide range sensors are to be deployed in a field of interest such that, at any given time, the moving object can be detected by only a small number of time-of-flight sensors. A cooperative and distributed solution is needed to achieve the common localisation task using, exclusively, the available range measurements. Based on local exchange of information, the Belief Propagation (BP) approach provides a powerful solution to sequentially fuse data collected in a distributed fashion from the sensor nodes. Both the NBP and the box-BP schemes are compared for this scenario over real and simulated data sets.

Part II is organised in two chapters. Chapter 5 provides a literature review on the use of sensor networks for source characterisation in several application fields, as well as an overview of advection-diffusion models for the cases of instantaneous and continuous release point sources. It ends with formulating the forward problem aiming at predicting the concentration measurements at the sensors' positions knowing the sources parameters. Chapter 6 then considers the problem of localising an unknown number of land mines using, exclusively, concentration measurements provided by a wireless sensor network. The two-step approach (mentioned earlier) consisting in sequentially determining the number of sources, then localising them is employed to solve the inverse problem, corresponding to the forward model formulated in Chapter 5. Two technical solutions are reported, a probabilistic Bayesian approach based on a Markov Chain Monte Carlo (MCMC) sampling scheme, and an optimisation least squares technique. The effectiveness of these schemes is tested on simulated data for two scenarios, assuming, respectively, the land mines as instantaneous and continuous release sources.

A general conclusion recapitulating the basic concepts and contributions of the thesis, as well as future work and perspectives, is at last given.

Much of the research in this dissertation has been published or submitted in the form of journal and conference papers. Results reported in chapters 2 and 3 are adapted from [START_REF] Haj Chhadé | Non parametric distributed inference in sensor networks using box particles messages[END_REF] while some of the results provided in chapter 6 have been established in [START_REF] Haj Chhadé | Localisation of an unknown number of land mines using a network of vapour detectors[END_REF][START_REF] Haj Chhadé | Multiple land mines localization using a wireless sensor network[END_REF].

BACKGROUND

Introduction

Message-passing algorithms built upon probabilistic graphical models have become important tools to infer in WSNs applications [START_REF] Paskin | Exploiting Locality in Probabilistic Inference[END_REF]. Furthermore, in order to increase the decision accuracy within a stochastic environment, filtering and inference issues should be involved.

In this chapter, we provide a short overview of concepts relevant to the subsequent parts of the thesis, as well as an introduction of the prior work related to the methods and algorithms derived later on. The first part of the thesis considers the problem of inference on arbitrary graphical models where we derive a novel message-passing algorithm defined in the bounded error framework. A graphical model is a statistical tool which can be associated with a distributed system, where the probability model is stored in a distributed fashion, such as sensor networks. A general description of sensor networks and the common issues involved with them (typically communication and energy constraints) is reported in section 1.2. Section 1.3 then defines the popular problem of filtering which will be seen afterwards as a particular case of the inference problem defined on a graphical model considered in this first part of the dissertation. Different non-parametric representations of probability quantities are employed in the probabilistic procedures we consider, or derive, throughout this manuscript, amongst which are the Monte Carlo particle-based representations, that are briefly explained in section 1.4. Finally, basic knowledge of graphical models and message-passing procedures is available through sections 1.5, 1.6 and 1.7 of this introductory chapter.

Sensor networks: a general overview

Sensor networks have become a popular solution in many civilian and military applications [START_REF] Senocak | Stochastic event reconstruction of atmospheric contaminant dispersion using Bayesian inference[END_REF][START_REF] Weimer | An approach to leak detection using wireless sensor networks at carbon sequestration sites[END_REF], they are useful whenever there is a need to monitor or to control physical quantities, such as temperature, brightness or pressure. In typical scenarios, sensors acquire data from their environment and often communicate the information collected to a processing centre.

Due to recent advances in wireless technology, networks of wireless sensors present several advantages over their wired counterparts [START_REF] Makhoul | Réseaux de capteurs : localisation, couverture et fusion de données[END_REF], mainly: they are cheaper because they require much less infrastructure, they are freed from the constraint of wiring which considerably limited the deployment of a large number of nodes, they are able to convey information to all nodes connected to the network and, finally, they are self-configurable and self-organising networks with almost no need for direct human intervention.

With these features, WSNs are quickly and cheaply deployed in environments which are dangerous for people, such as regions of conflicts [START_REF] Ridley | Decentralised ground target tracking with heterogeneous sensing nodes on multiple UAVs[END_REF], or areas which are difficult to access, as in environmental monitoring applications [START_REF] Lushi | An inverse Gaussian plume approach for estimating atmospheric pollutant emissions from multiple point sources[END_REF]; WSNs also make it possible to deploy a network consisting of a large number of sensors to monitor wider areas.

A WSN consists of a large number of units called sensor nodes. A node is generally composed of three main units [START_REF] Bahi | Improving lifetime and coverage through mobile beacon for high density sensor networks[END_REF][START_REF] Ihler | Inference in Sensor Networks: Graphical Models and Particle Methods[END_REF]: an acquisition unit, a processing unit and a communication unit. The acquisition unit obtains numerical measures on environmental parameters depending on the sensing ability of the sensor (e.g. seismic, optical, acoustic or meteorological) and forwards the recorded information to the processing unit. The processing unit is composed of a CPU and a specific operating system allowing the sensor node to perform some local computations. The communication unit sends and receives data via a radio communication medium thus enabling the sensor to exchange information with other nodes of the network. These three units are powered by a self-contained battery as depicted in figure 1.1. Other components may be added to a sensor node, such as a positioning module (GPS) or a mobility component to make it mobile. Despite the advantages of WSNs, a common concern in almost the majority of the work on sensor networks, is that of reducing energy consumption. Indeed, since wireless sensor networks are usually intended to raise information in hostile or inaccessible environments without human intervention, it is hard to find a source of energy other than self-contained batteries; therefore, the lifetime of a sensor node is equal to the lifetime of its battery. Consequently, limitations on the available energy should be considered when performing tasks such as estimation and inference in WSNs. Typically, communication requires many times more the amount of power required for sensing or computation [START_REF] Min | Energy-centric enabling tecumologies for wireless sensor networks[END_REF]. Algorithms for routing and distributed processing [START_REF] Arbula | Distributed algorithm for anchor-free network localization using angle of arrival[END_REF][START_REF] Langendoen | Distributed localization in wireless sensor networks: a quantitative comparison[END_REF][START_REF] Saad | At-angle: a distributed method for localization using angles in sensor networks[END_REF], only based on information gathered by neighbour nodes, i.e. within limited ranges, have arisen as an important key for energy optimisation of a WSN.

Bayesian filtering

A dynamic system is a physical system that evolves over time and can be described by a set of variables grouped into a state vector. State estimation methods aim at studying the behaviour of a dynamic system by calculating, at each time step, an estimate of the state vector, usually through combining information provided by different sensors. Filtering refers to the problem of sequentially estimating the states, or hidden variables, of a system as a set of observations become available [START_REF] Wan | The unscented kalman filter for nonlinear estimation[END_REF] (this problem is illustrated in figure 1.2). Consider a dynamic model described by the following system:

x k = f(x k-1 , v k ), (1.1) y k = g(x k , w k ), (1.2)
where:

• x k ∈ R nx , y k ∈ R ny
are respectively the state and the observation vectors at time step k.

• f is a transition function which allows to predict the state vector x k at time step k in terms of the previous state x k-1 and an i.i.d. process noise sequence, v k ∈ R nx ; equation (1.1) is usually referred to as the prediction model.

• g is a function that defines the relation between the state and the measurement vectors (x k and y k respectively) taking into consideration a measurement noise, w k ∈ R ny , denoting an i.i.d. measurement noise sequence; equation (1.2) is often called the observation model.

• X k = {x k , x k-1 , . . . , x 1 } and Y k = {y k , y k-1 , .
. . , y 1 } are vectors grouping, respectively, the states and the measurements up to time step k.

Within the Bayesian framework, a description of the state given the measurements, up to time step k, is provided by the posterior probability distribution p(X k |Y k ). Furthermore, in Bayesian filtering it is often enough determining, at each time step k, the posterior marginal distribution denoted as p(x k |Y k ) and given by:

p(x k |Y k ) = 1 α k p(x k |Y k-1 )p(y k |x k ), (1.3)
where

α k = p(y k |x k )p(x k |Y k-1 )dx k is a normalisation factor, p(y k |x k ) is the likelihood, p(x k |Y k-1
) is the prior and can be expressed as:

p(x k |Y k-1 ) = p(x k |x k-1 )p(x k-1 |Y k-1 )dx k-1 . (1.4) In equation (1.4), the PDF p(x k |x k-1
) is called the transition density (and is associated with the prediction model) while p(x k-1 |Y k-1 ) is the posterior at time step k. The recursion given by expressions (1.3) and (1.4), corresponding respectively to the time update and the measurement update steps, allows to sequentially compute the posterior marginals p(x k |Y k ), for each k; the iterative process is initialised with a prior p(x 0 ), which can be chosen as an uninformative uniform PDF over some region of the state space.

Monte Carlo methods

Let p(x) denote a probability density function (PDF) defined for x ∈ R n , and x i p ∼ p(x), i = 1, . . . N , an i.i.d. sample of size N . Monte Carlo methods rely on the fact that p(x) can be empirically represented by a set of independent samples drawn from p [START_REF] Johansen | A Tutorial on Particle Filtering and Smoothing: Fifteen years later[END_REF]. A representation of p(x), denoted as p(x), is given by the following expression:

p(x) = 1 N N i=1 δ x i p (x), (1.5)
where δ a (x) denotes the Dirac delta function located at a ∈ R n . The expectation of any function f of x, under the distribution p, can also be approximated as follows:

E p (f(x)) = f(x)p(x)dx ≈ 1 N N i=1 f(x i p ). (1.6)
It can be easily demonstrated that this estimate is unbiased [START_REF] Johansen | A Tutorial on Particle Filtering and Smoothing: Fifteen years later[END_REF].

Since, in some cases, it is not possible to draw samples directly from p(x), such as when this distribution is not given in a closed parametric form, Monte Carlo approximation allows to represent p using a collection of weighted particles, or samples, using a procedure such as importance sampling.

Let q(x) denote a proposal distribution (the next section shows how it is chosen) that can be evaluated and sampled, and let x i q ∼ q(x), i = 1, . . . , N , be an i.i.d. sample for q. Assuming that q(x) is absolutely continuous with respect to p(x), i.e., if q(x 0 ) = 0, for some x 0 ∈ R n , then p(x 0 ) = 0 as well. The relative likelihood of having drawn x i q from p(x) versus q(x) is given by:

ω i ∝ p(x i q ) q(x i q ) . ( 1.7)
The quantities ω i , i = 1, . . . , N , are called weights and are normalised as to obtain i ω i = 1; to compute these weights, only the ratio of p to q needs to be known. The resulting weighted collection of samples {ω i , x i q } represents the distribution p(x) and similarly to equation (1.6) an approximation of the expectation of a function f(x) is given by

f(x)p(x)dx ≈ 1 N N i=1 ω i f(x i q ).
(1.8)

Particle filtering

Different approaches for solving the filtering problem defined in section 1.3 have been proposed in the literature. Amongst the most popular filtering procedures we find the standard Kalman filter (KF) and its variants such as the extended Kalman filter (EKF) [START_REF] Wan | The unscented kalman filter for nonlinear estimation[END_REF][START_REF] Welch | An introduction to the kalman filter[END_REF], these algorithms make an assumption of Gaussian posterior PDFs to approximate the recursive Bayesian estimation described previously. Particle filters make no assumptions on the form of the probability densities of interest, they provide a Monte Carlo based procedure for sequentially estimating the posterior marginal distributions p(x k |y k , y k-1 , . . . , y 1 ) at each time step. The most basic form of particle filtering operates as follows: given a collection of weighted samples {ω j k-1 , x j k-1 } representing the posterior p(x k-1 |Y k-1 ) at time k -1, in the time update step, these particles are propagated through the prediction model by sampling from the conditional distribution p(x k |x k-1 ):

x j k ∼ p(x k |x j k-1 ), (1.9)
the proposal distribution is thus chosen to be the transition probability function. The resulting collection of particles {ω j k-1 , x j k } represents the prior at time k; the next step, known as the measurement update step, is to correct the weights using the likelihood at k according to what follows

ω j k ∝ ω j k-1 p(y k |x k ). (1.10)
The weighted collection of samples {ω j k , x j k } now represents the posterior p(x k |Y k ) at time k. This iterative procedure is initialised by drawing samples x j 1 ∼ p(x 1 ) from some prior distribution p(x 1 ).

A problem often encountered with this procedure is particle depletion: that is one or more particles end up dominating with a relatively great weight as we walk through the iterations. A popular solution for this problem is resampling. It consists in replacing the collection of weighted particles by a set of equally weighted particles, before propagating the particles through the transition probability as shown in figure 1.3. Intuitively, several copies of a single particle would exist, their number is proportional to the initial particle's weight. These copies, when propagated through the transition density, result in different new particles and thus adds variety to the samples representing the prior at time k.

Graphical models

Probabilistic graphical models [START_REF] Koller | Probabilistic Graphical Models, Principles and Techniques[END_REF] are special graphs that combine probability with graph theories in order to visualise conditional independence properties of joint PDFs [START_REF] Paskin | Exploiting Locality in Probabilistic Inference[END_REF]. They are important and useful tools to reason with probability models defined in large spaces [START_REF] Bishop | Pattern Recognition and Machine Learning, chapter 8: Graphical models[END_REF]. Typically, given a set of random variables, conditional independence properties are crucial for understanding the underlying structures of a joint PDF. Graphical models offer a framework to exploit these conditional independence properties and to represent joint probability distributions as a product of local factors where each factor depends only on a subset of variables. This factorisation offers many computational advantages, mainly a more compact representation of the PDF and an accelerated inference [START_REF] Paskin | Exploiting Locality in Probabilistic Inference[END_REF].

Graphical models have important characteristics [START_REF] Bishop | Pattern Recognition and Machine Learning, chapter 8: Graphical models[END_REF] that make them widely employed in fields such as statistics, statistic physics, computer vision and machine learning, where the problems of uncertainty and complexity are often encountered: (a) they represent the structure of a probability model; (b) they offer a visualisation of the conditional independence properties of the model; (c) they translate complex computations, required to perform tasks such as inference, into visible graphical manipulations.

Formally, if X = {x 1 , . . . , x N } denotes a set of random variables, x i , i = 1, . . . , N , p(X) a joint PDF defined over X, a probabilistic graphical model is a tool to represent the underlying structure of the joint PDF p; it is a special graph G = (V, E) where V is a set of nodes or vertices, and E is a set of edges, each node i ∈ V is associated with a random variable x i and each edge (i, j) represents a probabilistic relation between the random variables x i and x j . The graph structure describes the independence relationships between the variables qualitatively. The parameters Θ will be used to quantify the interdependency between nodes i and j.

Different types of graphical models exist: directed graphical models and undirected graphical models. Directed graphical models, also known as Bayesian Networks (BNs), express causal relationships between random variables. Undirected graphical models, also called Markov Random Fields (MRFs), encode constraints and correlations between variables [START_REF]Continuous State Graphical Models for Object Localization, Pose Estimation and Tracking[END_REF]. BNs are popular in statistics, artificial intelligence and machine learning societies while MRFs are widely employed in computer vision and statistical physics [START_REF] Ben-Gal | Encyclopedia of Statistics in Quality and Reliability, chapter Bayesian Networks[END_REF]. A third class, known as factor graphs, generalises both directed and undirected types of graphical models [START_REF] Kschischang | Factor graphs and the sum-product algorithm[END_REF][START_REF] Sudderth | Graphical Models for visual Object Recognition and Tracking[END_REF].

Directed graphical models

Directed graphical models, i.e. Bayesian Networks, are a class of graphical models that uses directed acyclic graphs (DAGs) to express causal relationships between random variables [START_REF] Ben-Gal | Encyclopedia of Statistics in Quality and Reliability, chapter Bayesian Networks[END_REF]. They characterise the way a joint PDF defined over a set of random variables, factors into a product of conditional probabilities imposed by the structure of the graph.

Let p(X) denote a joint probability distribution defined over a set X of N random variables, X = {x 1 , . . . , x N }. Using the chain rule, p can be expressed as:

p(x 1 , x 2 , ..., x N ) = p(x 1 ) N j=2
p(x j |x j-1 , . . . , x 1 ).

(1.11)

Many graphs, such as static or dynamic sensor networks, are sparse. This implies that, for the associated probability models, the probabilities of some of the variables x j are not dependent of all of their predecessors (i.e. x 1 , x 2 , . . . , x j-1 ) but only depend on some of them. This restricted set of variables is called the set of parents of x j and is denoted by pa(x j ). Using these independence properties, the following equality holds:

p(x j |x j-1 , . . . , x 1 ) = p(x j |pa(x j )),
and hence, equation (1.11) can be put into the following form:

p(x 1 , x 2 , ..., x N ) = j∈{1,...,N } p(x j |pa(x j )).
(1.12)

The directed graph describing this joint probability distribution is defined by two sets: the set of nodes associated with the random variables, and the set of edges representing a direct causal connection between the variables. Indeed, an arrow from node i, associated with random variable x i , to node j associated with x j , encodes the fact that a value taken by x j is dependent of the value taken by x i . The DAG structure provides a qualitative description of the model, the set of quantitative parameters Θ, which specifies the forms of the relationships between the random variables, groups the quantities Conversely, having the DAG describing the structure of a joint PDF, the conditional independence properties, i.e. the factorised form of the model, can be inferred. In fact, each node/variable is independent of its non-descendants1 in the graph given the state of its parents. The joint PDF can be factorised into a product of conditional probability distributions where each node is only dependent of its parents. In the example of figure 1.4, the joint PDF factors into:

θ x i |pa(x i ) = p(x i |pa(x i )).
p(X) = p(x 1 , x 2 , x 3 , x 4 , x 5 ) = p(x 5 |x 2 , x 3 , x 4 )p(x 4 |x 2 )p(x 3 |x 2 )p(x 2 |x 1 )p(x 1 ).
To ensure that this factorisation holds, an important restriction on the topology of the graph must be maintained, that is, no path, from any node of the graph, can lead back to the node itself along the directed edges; this property is guaranteed by the structure of the acyclic graph.

A final remark is that, although the arrows represent a direct causal connection between the variables, information can propagate in any direction. In fact, a causal relation from node i to node j implies that any information on x i may modify the knowledge one has on x j . Conversely any knowledge on x j may reflect an information about x i [START_REF]Continuous State Graphical Models for Object Localization, Pose Estimation and Tracking[END_REF].

Example 1.5.1. A Markov chain is a simple example of a directed graphical model. A set of random variables {x 1 , x 2 , . . . , x N } define a first order Markov chain if

p(x i+1 |x i , x i-1 , . . . , x 1 ) = p(x i+1 |x i ), ∀i ∈ {1, . . . , N -1}.
(1.13)

The corresponding graph is given in figure 1.5.

Undirected graphical models

Markov Random Fields (MRFs) are a class of graphical models that uses an undirected graph to encode conditional independence relationships between random variables [START_REF] Ben-Gal | Encyclopedia of Statistics in Quality and Reliability, chapter Bayesian Networks[END_REF].

Let G = (V, E) be an undirected graph representing a joint probability distribution, p(X) defined over a set of N random variables X = {x i } N i=1 , A, B and D refer to three disjoint subsets of V , X A , X B , X D the sets of random variables associated with the nodes in sets A, B and D respectively. If every path between A and B passes through some node in D, D is said to separate A and B [START_REF] Sudderth | Graphical Models for visual Object Recognition and Tracking[END_REF]. This graphical criterion, known as graph separation [START_REF] Paskin | Exploiting Locality in Probabilistic Inference[END_REF], encodes conditional independence between X A and X B given X D , thus:

p(X A , X B |X D ) = p(X A |X D )p(X B |X D ).
(1.14)

For a node i ∈ V , let Γ i ⊂ V refer to the set of neighbours of i, Γ i hence consists of the nodes of G connected to i by an edge, and designates a set of variables denoted as

X Γ i = {x j |j ∈ Γ i }.
The graph structure suggests that any random variable x i is conditionally independent of all other variables of the model given X Γ i :

p(x i |X V \i ) = p(x i |X Γ i ), (1.15)
where V \i refers to the set of nodes of G except for node i.

To specify the forms of the relationships between random variables, a parametrisation of MRFs is obtained by referring to the Hammersley and Clifford theorem. Recall beforehand that, in graph theory, a clique is a set of fully connected nodes in a graph (see figure 1.6 for an example). Theorem 1.5.1. Hammersley and Clifford Theorem: (adapted from [START_REF] Sudderth | Graphical Models for visual Object Recognition and Tracking[END_REF]) Let G = (V, E) be an undirected graphical model, where each vertex i ∈ V corresponds to a random variable x i , and C denotes a set of cliques of G; then, a probability distribution defined as the product of normalised positive functions (symmetric in their arguments) defined on the cliques is always Markov with respect to the graph.

p(X) ∝ c∈C ψ c (X c ).
(1.16)

Alternatively, any positive joint density function p(X) > 0, ∀X, which is Markov with respect to the graph, implies that there exist positive functions ψ c (symmetric in their arguments) such that the equality in equation (1.16) holds.

According to theorem 1.5.1, the joint distribution can be expressed as a product of factors, each of which is defined on a clique of the graph.

p(X) = 1 Z c∈C ψ c (X c ), (1.17)
where C is a set of cliques of G, X c designates the set of random variables associated with c ∈ C, ψ c (X c ) is a factor associated with c and is referred to as a potential function, and Z is a normalisation constant. Figure 1.7 shows an example of a MRF with N = 5 nodes/variables. One possible parametrisation of the graph is given by

p(X) = 1 Z ψ 235 (x 2 , x 3 , x 5 )ψ 245 (x 2 , x 4 , x 5 )ψ 12 (x 1 , x 2 ).
p(X) = 1 Z ψ 12 (x 1 , x 2 )ψ 23 (x 2 , x 3 )ψ 24 (x 2 , x 4 )ψ 25 (x 2 , x 5 ).
Example 1.5.3. For the example given in figure 1.7, choosing the cliques to be exclusively pairs of nodes connected by an edge, a different parametrisation of the MRF can be obtained:

p(X) = 1 Z ψ 12 (x 1 , x 2 )ψ 23 (x 2 , x 3 )ψ 24 (x 2 , x 4 )ψ 25 (x 2 , x 5 )ψ 35 (x 3 , x 5 )ψ 45 (x 4 , x 5 ).
The parametrisation of a MRF is not unique.

The problem of inference

In a typical scenario, two sets of random variables are considered: the set X and the set Y of respectively hidden and observable variables, and it is usually desired to obtain an information about the hidden quantities given the observations. Consequently, the set of nodes V , associated with the random variables describing the system, is split into two disjoint sets, V x and V y , associated respectively with the set X = {x 1 , x 2 , . . . , x N } of the hidden variables and the set Y = {y 1 , y 2 , . . . , y M } of the observable ones. Nodes in sets V x and V y are referred to as hidden (or latent) and evidence nodes respectively. The problem of inference on a graphical model refers to finding the posterior distribution, denoted as p(X|Y), of some, or all, hidden variables given the observations [START_REF]Continuous State Graphical Models for Object Localization, Pose Estimation and Tracking[END_REF][START_REF] Ying | Feedback message passing for inference in gaussian graphical models[END_REF]. Without loss of generality, let us develop further the expression of the posterior distribution p(X|Y) in the case of an undirected graphical model whose cliques are chosen to be exclusively pairs of nodes connected by edges, e.g. the case of a pairwise MRF. Let us also assume that the set Y of observations is actually a collection of local observations y i of individual node variables x i , i ∈ V , corrupted by some noise [START_REF] Ihler | Inference in Sensor Networks: Graphical Models and Particle Methods[END_REF] and that y i is independent of all variables in the model given x i . Using the Hammersley and Clifford theorem, the probability distribution p(X, Y) can be factorised into a product of factors, or potential functions, which, under the previous assumptions, can be partitioned into two groups: the first group corresponds to edges existing between the hidden variables and their associated local observations and the second group corresponds to edges interconnecting hidden variables. Hence, the posterior distribution can be expressed as a product:

p(X|Y) = p(X, Y) p(Y) ∝ (i,j)∈E ψ ij (x i , x j ) i ψ i (x i , y i ).
(1.18) Figure 1.9 illustrates the problem of inference on a graphical model defined in this section. In this example, the vector of hidden variables is X = (x 1 , x 2 , x 3 , x 4 , x 5 , x 6 , x 7 ) and the vector grouping the observations is Y = (y 1 , y 3 , y 5 ). Typically, nodes associated with Y should have been added to the graph, with edges interconnecting these evidence nodes with their corresponding hidden variables; for simplicity reasons, and since the observations are local to their associated hidden nodes, nodes for which an observation is available are shaded.

Belief Propagation algorithm

The Belief Propagation (BP), also known as the sum-product algorithm, is a statistical estimation tool for solving inference problems defined on graphical models. It exploits the conditional independence relationships encoded by the graph in order to calculate, in an exact or an approximated way, the marginal posterior distribution at a node of the model. In our exposition of the BP, we focus on pairwise MRF, and hence, we reason within the assumptions made in section 1.5.3. However, the methods underlying BP as described in this section are also adapted to message-passing in Bayesian Networks (BNs) and in factor graphs [START_REF] Wainwright | Graphical models, exponential families, and variational inference[END_REF][START_REF] Cheung-Mon-Chan | Réseaux Bayésiens et Filtres Particulaires pour l'égalisation adaptative et le décodage conjoints[END_REF][START_REF] Sudderth | Graphical Models for visual Object Recognition and Tracking[END_REF].

The inference goal considered here is to compute the posterior marginal distribution associated with a node t given the set of observations Y:

p(x t |Y) = X\xt p(X|Y)dX\x t .
(

Based on the factorisation (1.18), this marginalisation task is captured by BP in an iterative fashion: BP [START_REF] Yedidia | Message-passing algorithms for inference and optimization[END_REF] iteratively exchanges information between nodes of a graphical model through a message-passing procedure, the most common form of which is a parallel update scheme where, at each iteration, each node calculates outgoing messages to all of its neighbours simultaneously. A message, denoted by m i ts (x s ), is sent from a node t to one of its neighbours s, at an iteration i of BP, and is given by the following expression: The quantity

m i ts (x s ) ∝ ψ ts (x t , x s )ψ t (x t ) u∈Γt\s m i-
M i ts (x t ) ∝ ψ t (x t )R i ts (x t ) (1.22)
is referred to as the partial belief, it combines, at an iteration i, all the available information about t, from t itself and t's neighbouring nodes (except for the receiving node s). These two quantities will be helpful when we demonstrate the theoretic derivation of the box-BP in the next chapter.

Next to message-passing, the belief at iteration i, at node t, is computed as follows:

q i t (x t ) ∝ ψ t (x t ) u∈Γt m i ut (x t ). (1.23)
Beliefs are PDFs and hence normalised in order to integrate to unity [START_REF] Ihler | Inference in Sensor Networks: Graphical Models and Particle Methods[END_REF]. Figure 1.12 illustrates BP's belief computations for the example in figure 1. [START_REF] Bahi | Improving lifetime and coverage through mobile beacon for high density sensor networks[END_REF]. The belief at a node t approximates the posterior marginal of interest p(x t |Y). In fact, under some conditions the belief converges to the exact posterior marginal distribution [START_REF] Ihler | Inference in Sensor Networks: Graphical Models and Particle Methods[END_REF].

Loopy Belief Propagation

If the graphical model is tree-structured (i.e. is loop free), as in the case of a pairwise MRF, the Belief Propagation is guaranteed to converge after a finite number of iterations (at most equal to the length of the longest path in the graph [START_REF] Ihler | Inference in Sensor Networks: Graphical Models and Particle Methods[END_REF]). In this case, the belief at a node t will be exactly equal to the posterior marginal of interest p(x t |Y). Nevertheless, BP may also be applied to arbitrary graphical models. The same local Belief Propagation equations are iterated overlooking the presence of loops in the graph. This iterative procedure is referred to as loopy BP. In this case, the convergence of the sequence of messages is not guaranteed. Under some conditions however, fixed points will appear and in practice, these fixed points (beliefs) constitute reasonable approximations of the exact posterior marginal distribution associated with the node. Loopy BP is explained in details in [START_REF] Cheung-Mon-Chan | Réseaux Bayésiens et Filtres Particulaires pour l'égalisation adaptative et le décodage conjoints[END_REF].

Non-parametric Belief Propagation

The BP algorithm is a tool for performing exact or approximate marginalisation on arbitrary graphical models. This algorithm relies on two operations [START_REF] Noorshams | Belief propagation for continuous state spaces: Stochastic message-passing with quantitative guarantees[END_REF]: the first is a product of a collection of messages, the second is a convolution operation of this message product with a pairwise potential function as shown in equation (1.20). In many applications of graphical models, e.g. in computer vision [START_REF] Sudderth | Nonparametric belief propagation[END_REF], the hidden variables of interest are described by continuous non-Gaussian distributions; furthermore, the models governing information exchange between nodes may be seriously non-linear which makes BP's operations intractable. Non-parametric Belief Propagation (NBP) [START_REF] Sudderth | Nonparametric belief propagation[END_REF] is an attractive solution in the presence of such complex situations; it is considered, on one hand, as a generalisation of particle filtering applied to arbitrary graphical models and on the other hand as a stochastic approximation of BP [START_REF] Ihler | Inference in Sensor Networks: Graphical Models and Particle Methods[END_REF][START_REF] Sudderth | Nonparametric belief propagation[END_REF][START_REF] Noorshams | Stochastic belief propagation: A lowcomplexity alternative to the sum-product algorithm[END_REF].

The basic idea behind NBP is to use sample-based representations to approximate the operations of BP [START_REF] Sanghavi | Belief propagation and LP relaxation for weighted matching in general graphs[END_REF]. We explain the computations required by NBP by considering each operation of the message update step (see equation (1.20)) stated earlier. Recall that t and s refer to the sending and the receiving nodes respectively, while x t and x s are their associated random variables.

The message product operation

The first operation in the message update procedure is a message product which makes it possible to compute the partial belief at the sending node t according to equation (1.22). In NBP, messages m i ut (x t ) are represented by a collection {ω j ut , x j ut }, j = 1, . . . , N , of weighted samples. A representation of m i ut (x t ) is hence given by

mi ut (x t ) = N j=1 ω j ut δ(x -x j ut ), (1.24)
where δ is the Dirac impulse function. The product u∈Γt\s mi ut (x t ) is virtually guaranteed to be zero everywhere since the event of drawing exactly the same samples from all of the messages is far unlikely; it is thus necessary to have continuous estimates of the messages. Such an estimation can be determined using non-parametric density estimation methods [START_REF] Silverman | Density estimation for statistics and data analysis[END_REF]. In [START_REF] Ihler | Inference in Sensor Networks: Graphical Models and Particle Methods[END_REF], the effect of each sample is smoothed onto a nearby region using a Gaussian Kernel K h (.), so the message is estimated by a mixture of N Gaussian distributions:

m i ut (x t ) = N j=1 ω j ut K h (x -x j ut ). (1.25)
Performing the product of d messages, each represented by a mixture of N Gaussian distributions, and noting that the product of two Gaussian distributions is Gaussian [START_REF]Continuous State Graphical Models for Object Localization, Pose Estimation and Tracking[END_REF], the result is a Gaussian mixture containing N d components. To avoid complex computations, the author in [START_REF] Ihler | Inference in Sensor Networks: Graphical Models and Particle Methods[END_REF] proposes to sample exactly N particles from this product of Gaussian mixtures using a Markov Chain Monte Carlo (MCMC) method, namely the Gibbs sampler. The message product, i.e. the partial belief, is thus represented by a collection {Ω j ts , x j ts } N j=1 of weighted particles. If N is chosen sufficiently large, this collection of samples is guaranteed to accurately characterise the true product [START_REF] Silverman | Density estimation for statistics and data analysis[END_REF].

The convolution operation

The second operation in the message update expression is a convolution through which the available information about the sending node t, represented by the partial belief, is propagated through the pairwise potential function ψ ts (x t , x s ), specifying the relationship between x t and x s , in order to form a message providing some information about the receiving node s' local state.

The partial belief is represented by a collection {Ω j ts , x j ts } N j=1 of weighted particles. In order to approximate the integral in equation (1.20), a direct application of Monte Carlo approximation given in equation (1.6) can be employed. However, the pairwise potential function ψ ts (x t , x s ) can have a marginal influence on x t [START_REF] Sudderth | Graphical Models for visual Object Recognition and Tracking[END_REF], i.e. the integral

ζ ts (x t ) = ψ ts (x t , x s )dx s , (1.26)
is not equal to unity. NBP accounts for this influence by incorporating ζ ts (x t ) into the message product operation, the samples are thus drawn from the product

ζ ts (x t )M i ts (x t ) ∝ ζ ts (x t )ψ t (x t ) u∈Γt\s m i ut (x t ). (1.27)
To approximate the integral in equation (1.20), the resulting particles are then propagated through the pairwise potential function ψ ts (x t , x s ). Figure 1.13 summarises NBP's computations.

Conclusion

This chapter aimed at familiarising the reader with the material discussed in subsequent chapters of the thesis by providing an overview of the related prior work. It evoked general issues involved with sensor networks applications. It also presented a short explanation of Monte Carlo non-parametric representations of PDFs which are widely used in environments presenting complex and unknown uncertainties. The notion of a probabilistic graphical model as well as the problem of inference which is considered in the first part of the thesis were also defined. The general Belief Propagation message-passing algorithm as well as a non-parametric variant of it employing particle-based representations of probability distributions were also emphasised. This algorithm, known as Non-parametric Belief Propagation (NBP), provides a solution to the problem of inference on arbitrary graphical models associated with systems comprising high-dimensional random variables and/or continuous non-Gaussian uncertainties and severe non-linearities. Each iteration of NBP used an efficient sampling procedure to update kernel-based estimates of the true continuous messages. However, on one hand, the storage of a sufficiently large number of weighted particles in order to accurately represent a PDF, and on the other hand, the overhead assigned to computing Kernel density estimates and sampling from products of Gaussian mixtures, increases the complexity associated with NBP's computations. This difficulty was the motivation behind deriving a novel scheme of message-passing in the next chapter. This procedure, which we will call box-BP, is applicable in the presence of complex uncertainties and non-linearities, and uses interval-based approximations of the probability quantities.

BELIEF PROPAGATION IN INTERVAL ANALYSIS FRAMEWORK

Introduction

In the presence of complex statistical uncertainties and non-linearities, non-parametric particle methods offer impressive advantages over their parametric counterparts. We saw in the previous chapter how the NBP approximates the standard Gaussian BP algorithm using Monte Carlo methods and how it can be efficiently employed to infer on arbitrary graphical models in such complex situations. However, this algorithm requires intensive computations and the use of a large number of samples to accurately approximate the operations of the BP. In the context of temporal inference problems, i.e. filtering, the box particle filter [START_REF] Gning | An introduction to box particle filtering[END_REF] achieved a significant efficiency compared with the generic sampling importance resampling particle filter. In fact, the estimation accuracy provided by this last, using thousands of particles, can be achieved by the box-PF using only a few particles [START_REF] Gning | Integrated Tracking, Classification, and Sensor Management: Theory and Applications, chapter Particle Filtering Combined with Interval Methods for Tracking Applications[END_REF]. Inspired by the box particle filter (box-PF) which combines interval analysis with particle filtering to solve temporal inference problems, this chapter presents an algorithm that uses bounded error methods for solving the inference problem in an arbitrary graphical model. By using interval methods, the overhead associated with the computation of nonparametric density estimates and the procedure of sampling from the product of messages in the NBP is avoided. This results in simpler computations in addition to reducing the memory storage required since only a few box-particles can be used to accurately approximate a probability density function. This chapter is organised as follows. In section 2.2 we provide an overview of the fundamental concepts from the interval analysis framework. We then present a brief depiction of the box-PF in section 2.3. Later, we show how to implement the Belief Propagation procedure when the information is represented using boxes. The corresponding theoretical derivation, presented in section 2.4, is based on interpreting a box as a uniform PDF. Finally, the box-PF is seen as a special inference problem defined on a temporal graphical model in section 2.5 and a conclusion of this chapter is given in section 2.6.

Bounded error methods

For systems presenting complex or unknown statistical bounded uncertainties, interval estimation is a promising methodology. Initially introduced to propagate rounding errors in mathematical computations [START_REF] Moore | Interval Analysis[END_REF], applications of the interval analysis framework to state estimation problems have been widely investigated [START_REF] Abdallah | Adapting particle filter on interval data for dynamic state estimation[END_REF][START_REF] Gning | Mixture of uniform probability density functions for non linear state estimation using interval analysis[END_REF][START_REF] Abdallah | Box particle filtering for nonlinear state estimation using interval analysis[END_REF][START_REF] Jaulin | Nonlinear bounded-error state estimation of continuous-time systems[END_REF].

An interval is a closed and connected subset of R defined as

[x] = [x, x] = {x ∈ R | x ≤ x ≤ x}, ( 2.1)
where x and x refer to the minimal and maximal bounds of [x] respectively.

A box is a vector [x] in R n and can be seen as a Cartesian product of n intervals:

[x] = [x 1 ] × . . . × [x n ]. (2.2)
Henceforth, IR refers to the set of intervals in R, IR n the set of boxes in R n , and |[x]| denotes the length of an interval [x] ⊂ R.

Operations on intervals and boxes

Set-theoretic operations, such as intersection and union, are applicable to intervals. The intersection of two intervals is always an interval, whilst their union is not necessarily an interval. The interval union, denoted by , of two intervals [x], [y] ⊂ R is defined as:

[x] [y] = [[x] ∪ [y]] , (2.3)
where ∪ refers to the set-theoretic union operation. The symbol [.] denotes the interval hull operator returning, for any set S in R, the smallest interval enclosing S [START_REF] Jaulin | Applied Interval Analysis[END_REF]. Binary operations {+, -, ×, /} can be extended to intervals [START_REF] Jaulin | Applied Interval Analysis[END_REF]. Let refer to any binary operation and [x] and [y] denote any two intervals in R, the resulting interval [z] = [x] [y] is defined by the following:

[z] = [x] [y] = [{x y | x ∈ [x], y ∈ [y]}] . (2.4)
If the binary operation is continuous, e.g. in the case of the usual arithmetic operations, the set [

{x y ∈ R|x ∈ [x], y ∈ [y]} is an interval. Thus, [x] [y] = {x y ∈ R|x ∈ [x], y ∈ [y]}.
x] + [y] = [x + y, x + y], [x] -[y] = [x -y, x -y], [x] × [y] = [min(xy, xy, xy, xy), max(xy, xy, xy, xy)], [x] [y] = [x] × [1/y, 1/y], 0 / ∈ [y].
(2.5)

Elementary functions, such as exp, ln, cos and sin, can be also easily extended to intervals. Let f : R -→ R be a function, its interval counterpart will be denoted by [f ] and is defined as

[f ]([x]) = [{f (x)|x ∈ [x]}]. (2.6)
If f denotes a continuous and monotonic function,

[f ]([x]) is simply equal to f ([x]).
Example 2.2.2. ln is a continuous and strictly increasing function over R + * .

∀[x] ⊂ R + * , [ln]([x]) = [ln(x), ln(x)].
All operations on intervals can be extended to boxes. b. The interval union of [x] and [y] is defined as:

Example 2.2.3. Let S ⊂ R n , [x] and [y] two boxes in R n .
[x] [y] = ([x 1 ] [y 1 ]) × . . . × ([x n ] [y n ]).
c. Binary operations on real vectors and functions of real vectors may be extended to interval vectors. For instance,

[x] + [y] = ([x 1 ] + [y 1 ]) × . . . × ([x n ] + [y n ]).

Inclusion functions

Let f be a function

f : R n -→ R m , [x] a box in R n . The image f([x]) of [x]
by f is usually not a box (see figure 2.1), and its expression may be quite difficult to obtain. An inclusion function approximates f([x]) [START_REF] Gning | Integrated Tracking, Classification, and Sensor Management: Theory and Applications, chapter Particle Filtering Combined with Interval Methods for Tracking Applications[END_REF]. If [f] denotes an interval function from IR n to IR m then, by definition, [f] is said to be an inclusion function for f if

f ([x]) ⊆ [f ]([x]), ∀[x] ∈ IR n . (2.7)

Inclusion functions can be very pessimistic, an inclusion function

[f] is minimal if, for any x, [f]([x]) is the interval hull of f([x]
). The minimal inclusion function for f is unique and will be denoted by [f] * . Two main purposes of interval analysis [START_REF] Jaulin | Applied Interval Analysis[END_REF] are, first, finding an inclusion function

[f] such that, for most x, [f]([x]) is close to [f ] * ([x]
), and second, finding it with a convenient computational time. Several algorithms have been proposed to reduce the size of the box enclosing f ([x]). The next section presents one such method. such that the components f j , j ∈ {1, . . . , m}, of f represent constraints linking the variables {x i } n i=1 . The vector x belongs to a known prior domain [x] ∈ IR n . A constraint satisfaction problem (CSP) aims at finding the smallest box enclosing the set of all x ∈ [x] which satisfy the constraints f. A CSP is commonly denoted as H and can be formulated as:

Constraints satisfaction and contraction methods

Constraints satisfaction problems (CSPs)

Let x = (x 1 , x 2 , • • • , x n ) T denote a vector of n variables, x i ∈ R, i ∈ {1, . . . , n}, f = (f 1 , f 2 , • • • , f m ) T a multivalued function which satisfies f(x) = 0, (2.8)
H : (f (x) = 0, x ∈ [x]).
(2.9)

The solution set of the CSP H is given by

S = {x ∈ [x] | f(x) = 0} (2.10)
and is not necessarily a box. Solving a CSP refers to finding the smallest box

[x ] ⊂ [x] such that S ⊆ [x ] ⊆ [x] [56]. Figure 2.
2 illustrates an example of a CSP with two variables x 1 , x 2 , and a single constraint f . The initial box is [x

1 x 1 ] × [x 2 x 2 ]
. The solution of the CSP is given by the box

[x 1 x 1 ] × [x 2
x 2 ] and is highlighted in red, it is included in the initial box and it satisfies the constraint f .

Contractors

Contracting H refers to replacing [x] by a smaller domain

[x] such that S ⊆ [x] ⊆ [x].
Several contracting methods are described in [START_REF] Jaulin | Applied Interval Analysis[END_REF] (chapter 4), each of which is suitable for different types of CSPs. We will present a general contraction method known as Constraints Propagation (CP), this technique is simple, efficient and most importantly independent of non-linearities [START_REF] Gning | Integrated Tracking, Classification, and Sensor Management: Theory and Applications, chapter Particle Filtering Combined with Interval Methods for Tracking Applications[END_REF].

The CP method is based on the use of primitive constraints. A constraint is said to be primitive if it involves a single binary operation (such as +, -, ×, /) or a single elementary function (such as sin, cos, ln, exp). Constraint Propagation technique proceeds by contracting H with respect to each primitive constraint until convergence to a minimal domain. This method can be implemented using a forward/backward propagation (FBP) procedure [START_REF] Gning | Constraints propagation techniques on intervals for a guaranteed localization using redundant data[END_REF] which consists of two steps: the forward propagation and the backward propagation. The first step considers the direct forms of the equations. The second uses the inverse of the functions that appear in the equations.

Example 2.2.4. We consider the problem illustrated in figure 2.2 and contract the initial domain [

[x 1 ] × [x 2 ] = [x 1 x 1 ] × [x 2 x 2 ]
x 2 ] ← [x 2 x 2 ] ∩ [f ]([x 1 x 1 ]).
b. Backward propagation: Consider the inverse form of the constraint

x 1 = f -1 (x 2 ), [x 1 ] ← [x 1 x 1 ] ∩ [f -1 ]([x 2 ]).
Figure 2.3 shows the operations of CP for this example.

Example 2.2.5. This example is adapted from [START_REF] Gning | Integrated Tracking, Classification, and Sensor Management: Theory and Applications, chapter Particle Filtering Combined with Interval Methods for Tracking Applications[END_REF].

Consider a three dimensional CSP with a single constraint z = x exp(y). The initial domains are given by

[x] = [1, 7], [y] = [0, 1] and [z] = [0, 3]
. The constraint involves one arithmetic operation (×) and one elementary function (exp), it is not primitive. Introducing an auxiliary variable a, the constraint is decomposed into two primitive constraints:

a = exp(y), z = a • x.
(2.11)

An initial domain of a is merely [a] = [0, +∞]. The forward propagation step is achieved via the following two steps, denoted by F 1 and F 2 :

• F 1 : [a] ← [a] ∩ [exp]([y]) = [0, +∞] ∩ [1, e] = [1, e] • F 2 : [z] ← [z] ∩ [x] • [a] = [0, 3] ∩ [1, 7] • [1, e] = [1, 3].
The backward propagation step is achieved via three steps, B 3 , B 4 and B 5 :

• B 3 : [x] ← [x] ∩ ([z]/[a]) = [1, 7] ∩ [1, 3]/[1, e] = [1, 3] • B 4 : [a] ← [a] ∩ ([z]/[x]) = [1, e] ∩ [1, 3]/[1, 3] = [1, e] • B 5 : [y] ← [y] ∩ [ln]([a]) = [0, 1] ∩ [0, 1] = [0, 1].
The domains of the variables are reduced to , 3]. The FBP alternates between forward and backward propagation; in this example, it can be checked that the previous domains, obtained after one forward-backward propagation cycle, will no longer change after another iteration of FBP.

[x] = [1, 3], [y] = [0, 1] and [z] = [1

The box Particle Filter

Sequential Bayesian filtering (see section 1.3) consists in estimating the states of a system as a set of observations become available [START_REF] Wan | The unscented kalman filter for nonlinear estimation[END_REF]. Unlike the standard Bayesian solutions such as the Kalman Filter (KF) [START_REF] Welch | An introduction to the kalman filter[END_REF] and its variants [START_REF] Julier | A new approach for filtering nonlinear systems[END_REF], sequential Monte Carlo [START_REF] Johansen | A Tutorial on Particle Filtering and Smoothing: Fifteen years later[END_REF] techniques, also known as Particle Filtering (PF), make no assumptions on the form of the probability densities of interest. We have seen in section 1.4.1 that PFs [6] use sample-based representations in order to construct Monte Carlo approximations of the required integrals (e.g. in equation (1.4)) in the case of seriously non-linear and complex posterior distributions. However, the efficiency of particle filters (PFs) depends on the number of samples and on the proposal functions used for importance sampling. In fact, the imprecision on the available information requires potentially a large number of particles and thus a high computational complexity. The aim of the box-PF is to generalise particle filtering into the bounded error context; instead of propagating a large number of weighted point particles representing the posterior distribution, the Box-PF propagates a small number of weighted box particles through bounded error models and within the interval analysis framework.

The box-PF operates in four steps, for which the Bayesian justification was established by interpreting each box as a uniform PDF in [START_REF] Gning | Integrated Tracking, Classification, and Sensor Management: Theory and Applications, chapter Particle Filtering Combined with Interval Methods for Tracking Applications[END_REF]. These steps are: box particle initialisation, time update, measurement update and resampling.

k-1 ]} N =1
representing the state at time step k -1, and assuming that the system noise is enclosed in [v k ], the boxes at step k -1 are propagated using interval analysis tools through the transition function f (refer to equation (1.1)) [x

( ) k ] = [f]([x ( ) k-1 ], [v k ]) f or = 1, . . . , N,
where [f] is an inclusion function for f (refer to section 2.2.2 to recall the definition of an inclusion function). c) Measurement update step: Similarly to particle filtering, the weights of the predicted box particles are updated using the new measurement at time step k. For the box-PF, likelihood factors are calculated using innovation quantities [START_REF] Gning | An introduction to box particle filtering[END_REF]: the innovation for the -th box particle reflects the proximity between the measured box and the predicted box measurement. Hence, the innovation can be represented using the intersection between these two boxes. For each box particle, the predicted box measurement has the following expression (see equation (1.2)):

[y

( ) k ] = [g]([x ( ) k ], [w k ]), for = 1, . . . , N,
where [g] is an inclusion function for g and [w k ] is a bounded measurement noise. The innovation is given by: [r

( ) k ] = [y ( ) k ] ∩ [y k ].
In the bounded error context, the likelihood is calculated based on the following idea: if the predicted box measurement does not intersect with the corresponding measured box, this box particle has a likelihood factor equal to zero. In contrast, if the predicted box measurement is included in the corresponding measured box, this box particle has a likelihood close to one [START_REF] Gning | An introduction to box particle filtering[END_REF].

Furthermore, a contraction step (see section 2.2.3) is performed in order to eliminate the inconsistent part of the box particles with respect to the measured boxes, and to preserve an appropriate size of the boxes. The box likelihood is thus given by:

L ( ) k = nx j=1 L ( ),j k ,
where n x represents the dimension of the state and the likelihood factor according to a dimension j of x is given by

L ( ),j k = |[x ( ) k+1 (j)]| |[x ( ) k+1 (j)]| . The term [x ( ) k+1 (j)]
represents the new -th box particle after the contraction step. d) Resampling step: Similarly to the PF algorithm, a resampling step is also added in order to introduce variety into the box particles. Different resampling algorithms exist; in [START_REF] Gning | Integrated Tracking, Classification, and Sensor Management: Theory and Applications, chapter Particle Filtering Combined with Interval Methods for Tracking Applications[END_REF] for instance, the multinomial resampling combined with a new subdivision step is employed, in this scheme, after the resampling step, each box is divided by its corresponding number of realisations in order to obtain smaller boxes around regions with high likelihoods.

The computations of the box-PF are recapitulated in figure 2.4.

Theoretic derivation of the box-Belief Propagation

This section presents the main theoretical contribution of part I of the thesis: the derivation of the box-Belief Propagation (box-BP), which is a message-passing algorithm used to infer on graphical models when probability quantities are represented using boxes. The box-BP is indeed an extension of the BP algorithm into the bounded error framework.

Inspired by the theoretic derivation of the box-PF presented in [START_REF] Gning | Integrated Tracking, Classification, and Sensor Management: Theory and Applications, chapter Particle Filtering Combined with Interval Methods for Tracking Applications[END_REF][START_REF] Gning | An introduction to box particle filtering[END_REF], we show the theoretic derivation of BP in interval analysis framework by interpreting a box as a uniform PDF. Note, beforehand, that an advantage offered by this interpretation is that, while strictly speaking, a set of samples constitute a "representation" of a PDF, a set of boxes constitute a direct approximation of the PDF.

Going back to equation (1.20) representing BP's message update (for a message sent from a node t to its neighbour s, at iteration i of BP), recall that this iteration can be divided into two operations: the message product, or the computation of the partial belief given in equation (1.22), and the convolution operation. We demonstrate the extension to intervals of each operation in the following two subsections. A summary of the operations of the novel box-BP is illustrated in figure 2.5.

The message product operation

Let Γ t denote the set of neighbours of the sending node t, u l ∈ Γ t \s a node connected to t by an edge, except for the receiving node s, and d = card(Γ t \s); then the message product can be expressed as:

R i-1 ts (x t ) ∝ u∈Γt\s m i-1 ut (x t ) = d l=1 m i-1 u l t (x t ). (2.12)
In interval analysis framework, the messages are represented using N weighted boxes, thus, a message m u l t (x t ) received by node t from its neighbour u l , l = 1, . . . , d, is represented by a cloud of boxes {ω p l u l t , [x p l u l t ]} N p l =1 . Note that, for simplicity reasons, in this representation we skip the index i referring to the message passing iteration number on the particles and their weights. We are also assuming that all messages are represented using the same number of boxes, N . This reasoning may however be easily extended to a more general case with varying number of box particles per node/message.

Since each box is interpreted as a uniform probability distribution, we have:

m i-1 u l t (x t ) = N p l =1 ω p l u l t U [x p l u l t ] (x t ), for l = 1, . . . , d, (2.13)
where U [x] denotes the uniform PDF over the box [x]. Let P denote the vector of indexes (p 0 , p 1 , . . . , p d ) and I = {1, 2, . . . , N }. Replacing expression (2.13) in equation (2.12), we obtain

R i-1 ts (x t ) ∝ l ( N p l =1 ω p l u l t U [x p l u l t ] (x t )) ∝ P ∈I d ω p 1 u 1 t . . . ω p d u d t U [x p 1 u 1 t ] . . . U [x p d u d t ] (x t ). ( 2

The product U

[x p 1 u 1 t ] . . . U [x p d u d t ] (x t
) is then different to zero if x t belongs to the intersection of the supports of its terms. This product may hence be modelled using a uniform PDF given by

U [x p 1 u 1 t ] . . .U [x p d u d t ] (x t ) = U [x p 1 u 1 t ]∩...∩[x p d u d t ] × |[x p 1 u 1 t ] ∩ . . . ∩ [x p d u d t ]| |[x p 1 u 1 t ]| . . . |[x p d u d t ]| . (2.15) Then R i-1 ts (x t ) ∝ P ∈I d ω p 1 u 1 t . . . ω p d u d t U [x p 1 u 1 t ]∩...∩[x p d u d t ] × |[x p 1 u 1 t ] ∩ . . . ∩ [x p d u d t ]| |[x p 1 u 1 t ]| . . . |[x p d u d t ]| . (2.16)
The number of possible assignments of P is N d . Thus, the sum above contains at most N d terms. In practice, this number is much less since some assignments result in empty intersections while others result in coincident boxes.

Let us denote by Q the set of assignments of P , for which [

x p 1 u 1 t ] ∩ . . . ∩ [x p d u d t ] = ∅.
We will also adopt the following notations for an assignment

P k ∈ Q, where k = 1, . . . , card(Q): ω p 1 u 1 t . . . ω p d u d t × |[x p 1 u 1 t ] ∩ . . . ∩ [x p d u d t ]| |[x p 1 u 1 t ]| . . . |[x p d u d t ]| = ω k ts , (2.17)
[x p 1 u 1 t ] ∩ . . . ∩ [x p d u d t ] = [x k ts ]. (2.18)
The message product given in equation (1.21) is now represented by the cloud of boxes {ω k ts , x k ts }. Based on equation (2.16), we can define an algorithm to perform the message product operation.

Example 2.4.1. Algorithm 1 describes the method for combining two messages according to equation (2.16); note that both operations, (×) and (∩) are associative operations. In this algorithm, each message is represented using a certain number of weighted boxes, e.g.,

{ω i 1 , [x i 1 ]} N 1 i=1 , {ω j 2 , [x j 2 ]} N 2
j=1 denote the collection of weighted boxes representing message 1 and message 2 respectively. The cloud {ω k res , [x k res ]} Nres k=1 denotes the resulting set of boxes and the corresponding weights according to equations (2.17) and (2.18), iscoincident() is a function that takes as arguments two boxes and returns 1 if the boxes are coincident according a threshold and 0 if not.

Computation of the partial belief

The partial belief is now given by the following expression:

M i-1 ts (x t ) ∝ ψ t (x t ) card(Q) k=1 ω k ts U [x k ts ] (x t ). (2.19)
where, for simplicity and without loss of generality, ψ t (x t ) replaces ψ t (x t , y t ) under the assumption of an observation y t of the same nature as the state x t as stated in section 1.7. Within the bounded error framework, we assume that y t is known within a bounded box [y t ] and the potential ψ t (x t ) is approximated equal to the uniform PDF U [y t ] (x t ).

Algorithm 1 Message product

Input: {ω i 1 , [x i 1 ]} N 1 i=1 , {ω j 2 , [x j 2 ]} N 2 j=1 Output: {ω k res , [x k res ]} Nres k=1 k ← 0 for i = 1, . . . , N 1 do for j = 1, . . . , N 2 do if [x i 1 ] ∩ [x j 2 ] = ∅ then k = k + 1 [x k res ] = [x i 1 ] ∩ [x j 2 ] ω k res = ω i 1 × ω j 2 × |[x i 1 ]∩[x j 2 ]| |[x i 1 ]|.|[x j 2 ]| end if end for end for loop k 1 , k 2 ∈ {1, . . . , k} if iscoincident([x k 1 res , [x k 2 res ]
) then keep only one occurrence sum up the corresponding weights end if end loop Weights normalisation:

ω k res ←-ω k res /sum(w res ) N res = length(ω res )
In the general case, the potential function ψ t (x t , y t ) can be represented with the uniform PDF U [y t ] (g(x t )), where g is a function linking the local observation y t at node t to the local state x t (observation model), y t = g(x t ). Using contraction techniques [START_REF] Jaulin | Applied Interval Analysis[END_REF], the resulting weighted boxes

{ω k ts , [x k ts ]} card(Q) k=1
can be contracted and re-weighted. Under the assumption of an observation of the same nature as the state, the contraction reduces into an intersection between the resulting boxes and the observation. Algorithm 2 summarises the procedure of contracting the boxes

{ω k ts , [x k ts ]} card(Q) k=1
using the measurement [y t ], the resulting weighted boxes are denoted by {ω k ts,c , [x k ts,c ]} Nc k=1 where subscript c stands for contracted. In this procedure [g] denotes an inclusion function for g.

As a final remark, note that the belief at node t as given in equation (1.23) can be computed using the same procedure derived in this section.

Algorithm 2 Contraction using the local observation

Input: {w k ts , [x k ts ]} card(Q) k=1 , [y t ] Output: {w k ts,c , [x k ts,c ]} Nc k=1 for k = 1, . . . card(Q) do Predicted measurement: [y k ] = [g]([x k ts ] Innovation: [r k ] = [y k ] ∩ [y t ] Box particle contraction: if [r k ] = ∅ then contract [x k ts ] using CP algorithm the resulting box is denoted [x k ts,c ] else [x k ts,c ] = ∅ end if Re-weighting: w k ts,c ←-w k ts × |[x k ts,c ]| |[x k ts ]|
Weights normalisation: w k ts,c ←-w k ts,c /sum(w ts,c ). end for N c = length(w ts,c = 0)

The convolution operation

The next step in the message update iteration is the convolution of the partial belief with the pairwise potential function:

m i ts (x s ) = ψ ts (x t , x s )M i-1 ts (x t )dx t = ψ ts (x t , x s ) Nc k=1 ω k ts,c U [x k ts,c ] (x t ) = Nc k=1 ω k ts,c [x k ts,c ] ψ ts (x t , x s )U [x k ts,c ] dx t = Nc k=1 ω k ts,c 1 |[x k ts,c ]| [x k ts,c ] ψ ts (x t , x s )dx t .
(2.20)

To evaluate the integral in equation (2.20), consider the general case where the potential ψ ts (x t , x s ) can be represented as a transition function by which we can pass from x t to x s , i.e. x s = f (x t , v; e) where v denotes a noise and e refers to a vector of constants (e.g. some known parameters of the model), and let [f ] be an inclusion function for f . We assume that the noise v is bounded in the box [v].

Then, by definition of an inclusion function, we have

∀x t ∈ [x k ts,c ], x s ∈ [f]([x k ts,c ], [v]; [e]) f or k ∈ {1, . . . , N c }, i.e. ψ ts (x t , x s )U [x k ts,c ] (x t ) = 0, ∀x s / ∈ [f]([x k ts,c ], [v]; [e]).
(2.21) Equation (2.21) shows that for any transition function f, using interval analysis techniques, the support for the PDF terms

[x k ts,c ] ψ ts (x t , x s )U [x k ts,c ] dx t can be approximated by [f]([x k ts,c ], [v]; [e]
) and thus:

[x k ts,c ] ψ ts (x t , x s )U [x k ts,c ] dx t ≈ U [f]([x k ts,c ],[v],[e]) . (2.22)
Note that, it is shown in [START_REF] Gning | Integrated Tracking, Classification, and Sensor Management: Theory and Applications, chapter Particle Filtering Combined with Interval Methods for Tracking Applications[END_REF] that the approximation (2.22) can be done more precisely at a computation cost using a mixture of uniform boxes, e.g., a cloud of box particles.

For simplicity and without loss of generality only one box is used in this chapter.

Based on equations (2.20) and (2.22), we can perform the message update procedure, for a message sent from node t to its neighbour s, according to Algorithm 3: once the incoming messages to the sending node t are combined as depicted in Algorithm 1 and contracted using Algorithm 2, the collection of weighted boxes {ω k ts,c , [x k ts,c ]} Nc k=1 , combining all information about x t , are propagated through the model x s = f(x t , v, e) to obtain an information about x s . The resulting set {ω ts,r , [x ts,r ]} N =1 , where the subscript r stands for "resulting", represent the message sent from t to s. Note that for the first iteration of message-passing between neighbouring nodes, messages are initialised using the local marginal PDF since, obviously, no previous messages exist. If no observation is available for the sending node t, no message can be sent, at the first iteration, from t to any of its neighbours. If, however, an observation is available for t, we use it to generate uniformly N boxes mutually disjoint and equally weighted with

ψ k-1,k (x k-1 , x k ) = p(x k |x k-1 ), (2.23) ψ k (x k , y k ) = p(y k |x k ).
(2.24)

In the first equation, we recognise the transition probability density while in the second the likelihood at time instant k.

In the first iteration, a message m 01 (x 1 ) will be sent from node x 0 to node x 1 , it is initialised using a prior domain [x 0 ] for x 0 : [x 0 ] is split into N , equally weighted and mutually disjoint, boxes {[x 0 ]} N =1 . This corresponds to the box-particle initialisation step. To form the message m 01 (x 1 ), the resulting boxes are propagated through the transition function describing the pairwise potential function ψ 01 (x 0 , x 1 ) = p(x 1 |x 0 ) according to Algorithm 3. Thus the message m 01 (x 1 ) corresponds to the prior for x 1 . This corresponds to the time update step.

The belief at node x 1 is computed by the product of the only message received by x 1 from x 0 , that is the prior at time instant 1, with ψ 1 (x 1 , y 1 ) = p(y 1 |x 1 ) which is the likelihood, the result is the posterior at time instant 1. In this step, innovation quantities are calculated, boxes are contracted and re-weighted according to Algorithm 2. This corresponds to the measurement update step. In order to introduce variety into the box particles, a resampling step can be added. Note that the partial belief at time instant 1 is equal to the belief at this node.

In general, the interpretation of the messages and (partial) beliefs at time step k is given by

M k-1,k (x k-1 ) = p(x k-1 |y k-1 , . . . , y 1 ), the posterior at k -1, m k-1,k (x k ) = p(x k |y k-1 , . . . , y 1 )
, the prior at k.

In terms of messages and beliefs, the box particle filter can be described as follows: Given a collection of weighted samples

{ω k-1 , x k-1 } N =1 representing the partial belief M k-1,k (x k-1 ), which is the posterior at time step k-1, the message m k-1,k (x k ) is obtained by propagating {x k } N
=1 through the transition density, given by the process model and representing the pairwise potential function. The resulting cloud of boxes represents the prior at time step k. Next, the belief at node x k is computed by multiplying the only message received by

x k with the local likelihood ψ k (x k , y k ) = p(y k |x k ).
The result is the posterior at k.

Consequently, the box-particle filter is seen as a special case of the box-Belief Propagation, applicable to temporal inference problems. Equivalently, the box-BP generalises the box-PF techniques to arbitrary structured graphical models.

Conclusion

This chapter presented the main theoretical contribution of part I of this dissertation. It introduced a novel scheme of message-passing using interval representation of probability quantities. The purpose was to derive a BP applicable for graphical models associated with general statistical probability densities while avoiding the overhead of storing thousands of particles and the use of Kernel density estimates and the Gibbs sampler at each iteration of the NBP introduced in section 1.7. The efficiency of this novel box-BP will be tested in the following two chapters, on two different applications: the problem of self localisation in static WSNs and that of continuously localising a mobile target moving within a network of wireless range sensors.

SELF LOCALISATION IN WIRELESS SENSOR NETWORKS

Introduction

Many industrial, scientific and domestic applications employ sensor networks whenever there is a need to monitor, and possibly control, physical quantities. Sensors acquire information from their environment (e.g. temperature or pressure, depending on their sensing ability) and usually communicate the data collected to a processing centre. However, a vast majority of applications in sensor networks deploys a large number of sensors randomly, usually due to the hostility of the area to be monitored, or its immensity. The sensors are thus randomly dropped, from a flying agent for instance, and fall into the region of interest with practically no information about their positions. The localisation of the sensors is often necessary to make the data collected informative.

In this chapter, we consider the problem of self localisation in wireless sensor networks. Since all sensors cannot be equipped with a positioning unit, e.g. GPS module, due to cost and energy constraints, each sensor is alternatively equipped with a transmitterreceiver module and communicates with neighbouring sensors. Hence, each sensor is within the vicinity of only a few nodes of the network, and its position can be estimated by exchanging information with several nearby sensor nodes, with which a communication can be established. An attractive formalism to represent this system while visualising the local independence conditional relationships and the underlying structure of the associated joint PDF is possible through the use of graphical models. A literature review on self localisation in sensor networks is provided in section 3.2 then the focus is oriented towards MDS (multidimensional scaling) based localisation algorithms in section 3.3. Since the problem of auto-localisation in wireless sensor networks is pre-eminently a distributed inference problem, we formulate it as inference on a graphical model [START_REF] Ihler | Nonparametric belief propagation for self-localization of sensor networks[END_REF] and use the BP algorithm to perform the task of inference in section 3.4, the effectiveness of the novel box-BP approach is then tested on simulated data. Section 3.5 at last concludes this chapter.

Literature review

The problem of auto-localisation in WSNs have been widely investigated and a great deal of methods were proposed in the literature, most of which rely on the premise that a sensor node is able to communicate with its nearby counterparts in the network. Several classifications of localisation methods have been proposed. We cite for instance range-based versus range-free approaches. Range-based methods assume that sensors have capabilities to estimate distances with their neighbours, using technologies such as RSSI (Received Signal Strength Indicator), TOA (Time of Arrival) and TDoA (Time Difference of Arrival) [START_REF] Ihler | Nonparametric belief propagation for self-localization of sensor networks[END_REF], or angles, using AoA (angle of Arrival) [START_REF] Saad | At-angle: a distributed method for localization using angles in sensor networks[END_REF]. These techniques vary in their accuracy and complexity and can be divided into radio frequency (RF) ranging and acoustic ranging. The RF ranging relies on the premise that, by measuring the received signal strength, a receiver can determine its distance to a transmitter. The second class of rang-ing schemes measures the time difference of arrival of emitted and received signals [START_REF] Gustafsson | Localization in sensor networks based on log range observations[END_REF]. In contrast to range-based algorithms, range-free [START_REF] He | Range-free localization schemes for large scale sensor networks[END_REF] schemes do not need point-to-point distance or angle estimation for positioning. Approaches in [START_REF] Essoloh | Localisation distribuée dans les réseaux de capteurs sans fil par résolution d'un problème quadratique[END_REF][START_REF] Shang | Localization from connectivity in sensor networks[END_REF] for instance do not involve any range or angle measurement but rather rely on measures of connectivity.

Another classification subdivides the approaches for self localisation in wireless sensor networks to anchor-based and anchor-free. In both schemes, nodes must themselves determine their respective positions through cooperation techniques. In fact, the goal of a self localisation algorithm is to calculate the coordinates of each sensor based on proximity measures. The computed coordinates can be global, which requires the position of a number of anchors to be known a priori, these approaches are known as anchorbased [START_REF] Bahi | A mobile beacon based approach for sensor network localization[END_REF]. Other methods create a relative map without the use of anchors and are called anchor-free [START_REF] Arbula | Distributed algorithm for anchor-free network localization using angle of arrival[END_REF]. In range-based techniques, the placement of the anchors can often have a significant impact on the solution: it was found that the location accuracy improves if the anchors form a convex polygon around the network [START_REF] Langendoen | Distributed localization in wireless sensor networks: a quantitative comparison[END_REF]. Additional anchors placed at the centre of the network are nonetheless useful. Note that, other algorithms, related to auto-localisation in WSNs, were proposed to allow each node to locate its neighbours qualitatively; work in [START_REF] Heurtefeux | Localisation collaborative pour réseaux de capteurs[END_REF] offers one such example of a localised algorithm whose purpose is, for each node, to classify its neighbours into one of three categories: very close, near and far, rather than estimating distances separating them.

Each sensor has limited resources (e.g. bandwidth, battery energy, memory capacity, emission power), it can detect and communicate with other nodes in the network only within some maximum span. Early approaches to solving the problem of localisation in wireless sensor networks were proposed in a centralised environment. However, this strategy does not comply the principal energy constraint and is not appropriate for largescale networks. Indeed, local data processing has low energy cost, (see [START_REF] Cook | Soc issues for rf smart dust[END_REF] for details about energetic cost of RF communication in function of distance). In addition, the reliability of the centralised approach is low because a failure in the main processing unit causes the entire system to fail. Therefore, distributed approaches [START_REF] Ihler | Nonparametric belief propagation for self-localization of sensor networks[END_REF][START_REF] Essoloh | Localisation distribuée dans les réseaux de capteurs sans fil par résolution d'un problème quadratique[END_REF] were proposed. A third approach is to assist the nodes of the network by a mobile anchor enabling them to locate themselves. This approach [START_REF] Bahi | A mobile beacon based approach for sensor network localization[END_REF] has many advantages in terms of energy and location accuracy but its implementation can be expensive and is not always feasible as some areas can be hostile to moving robots.

We can also distinguish optimisation methods [START_REF] Essoloh | Localisation distribuée dans les réseaux de capteurs sans fil par résolution d'un problème quadratique[END_REF][START_REF] Shang | Localization from connectivity in sensor networks[END_REF] and probabilistic schemes [START_REF] Rong | Angle of arrival localization for wireless sensor networks[END_REF]. While optimisation approaches provide a single optimal solution that minimises or maximises some criterion, the solution in probabilistic methods takes the form of a probability distribution therefore allowing to quantify the uncertainty on the estimated positions. Work in [START_REF] Rong | Angle of arrival localization for wireless sensor networks[END_REF] presents a probabilistic method for locating sensors based on AoA technology, in this article the error on the measured angles is modelled using a Gaussian PDF. However, for some environments, uncertainties which are usually modelled as probabilistic can be complex; non-parametric methods offer an attractive framework to address these kind of problems [START_REF] Honeine | Auto-localisation dans les réseaux de capteurs sans fil par régréssion de matrices de gram[END_REF]. In [START_REF] Ihler | Nonparametric belief propagation for self-localization of sensor networks[END_REF] for instance, the auto-localisation problem is formulated as inference on a graphical model. The Non-parametric Belief Propagation (NBP) algorithm described in section 1.7 is used to combine the information obtained from a global positioning system, with measures of relative distances between neighbouring sensors. Recall that the NBP algorithm is itself a variant of BP, where a set of particles is used to represent probability quantities in a non-parametric way.

Furthermore, in some problems, uncertainties or noise characteristics are unknown or complex; instead, only minimum and maximum values of the noise are available, e.g. quantised measurements. The interval analysis framework offers promising methodologies for reasoning in the presence of unknown or complex statistical but bounded noises [START_REF] Gning | Integrated Tracking, Classification, and Sensor Management: Theory and Applications, chapter Particle Filtering Combined with Interval Methods for Tracking Applications[END_REF].

In this chapter, we propose a Bayesian probabilistic scheme for self localisation in WSNs based on a variant of BP algorithm where information is represented using a collection of boxes (intervals in the case of real variables). The use of this approach involves simplicity in modelling the information and memory optimisation. In fact, we will demonstrate that the use of interval representations in our approach offers many advantages, basically: a. a reduction of the memory space required to store a PDF (thousands of particles are needed to efficiently represent a PDF using Monte Carlo methods while only a few box-particles are required to approximate a probability distribution).

b. energy saving since less information is exchanged between communicating sensors, this also implies a reduction of the required bandwidth.

c. using interval techniques results in simpler and faster computations and thus more time saving.

d. a set of boxes constitute a direct approximation of the PDF while a set of samples constitute a representation of a PDF. This fact allows to save the energy and time needed to estimate the PDF using Kernel Density Estimation (KDE) techniques in the case of the standard sample-based NBP. This also introduces less approximations since the kernel is artificially added to the particles.

Multidimensional scaling for self localisation in WSNs

Multidimensional scaling (MDS) is a data analysis technique that transforms proximity information into a geometric embedding [START_REF] Michael | Springer Handbooks Comp.Statistics, Handbook of Data Visualization, chapter Multidimensional Scaling[END_REF]. Formally, given n points in a p-dimensional space E and a measure of dissimilarity reflecting, for a pair (i, j), i, j ∈ {1, . . . , n}, the remoteness between points i and j, the purpose of the classical Multidimensional scaling is to re-place these n points in a new subspace (which is often Euclidean) of a size q ≤ p in such a manner that the inter-point distances in this subspace approximate at best the dissimilarities. MDS is at the heart of several approaches for self localisation in WSNs, the most basic form of which is the MDS-MAP(C) algorithm introduced in [START_REF] Shang | Localization from connectivity in sensor networks[END_REF] in a centralised context. In a typical scenario for auto-localisation of wireless sensors, each sensor has a limited communication range allowing it to only compute the distances separating it from its neighbours. In [START_REF] Shang | Localization from connectivity in sensor networks[END_REF], MDS is used to solve the localisation problem after a preliminary treatment which makes it possible to determine a dissimilarity matrix. The dissimilarity matrix is in fact a distance matrix constructed by calculating the shortest path between each pair of nodes in the network; for instance, the algorithm of Dijkstra can be employed to determine optimal paths between pairs of nodes along the arcs of the network within a time complexity O(N 3 ). MDS-MAP(C) operates in three steps: a. Compute the dissimilarity, i.e. distance, matrix given the network connectivity information.

b. Use classical MDS to derive relative node coordinates that fit those distances. The network is reconstructed by only considering the q eigenvectors having the highest energy (q = 2 or 3 depending on whether the localisation procedure operates in 2D or 3D). Note that MDS offers the best approximation by a lower rank matrix in terms of least squares [START_REF] Essoloh | Localisation distribuée dans les réseaux de capteurs sans fil par résolution d'un problème quadratique[END_REF].

c. Reset the absolute positions of the nodes using the known coordinates of a number m of anchors within O(m) + O(n) time complexity.

MDS-MAP(C) offers a good performance on networks with relatively uniform node density. On irregular networks, its performance degrades since the shortest path does not correspond well to the Euclidean distance. To tackle this difficulty, the same reference [START_REF] Shang | Localization from connectivity in sensor networks[END_REF] proposes a variant of the MDS-MAP(C) employing a distributed implementation of it. The resulting algorithm is referred to as MDS-MAP(P) and operates by building, for each node, a local map of the small sub-network within the vicinity of it, and then merging all local maps together to form a global one.

Other approaches employing MDS for auto-localisation in WSNs can be found in [START_REF] Costa | Distributed weighted-multidimensional scaling for node localization in sensor networks[END_REF][START_REF] Deleskog | Robust nls sensor localization using mds initialization[END_REF][START_REF] Vivekanandan | Ordinal mds-based localization for wireless sensor networks[END_REF]. In [START_REF] Costa | Distributed weighted-multidimensional scaling for node localization in sensor networks[END_REF] for instance, a scalable and distributed weighted-multidimensional scaling (dwMDS) algorithm is proposed to adaptively emphasise the most accurate range measurements and thus account for communication constraints within the network. The sensors localisation task is formulated in [START_REF] Deleskog | Robust nls sensor localization using mds initialization[END_REF] as a non-linear least squares (NLS) problem and is robustly initialised using MDS. A different type of MDS (called ordinal MDS) is used in [START_REF] Vivekanandan | Ordinal mds-based localization for wireless sensor networks[END_REF] for self localisation in WSNs, the ordinal MDS differs from its classical counterpart by that it only requires, for each pair of nodes, a monotonicity constraint between the actual Euclidean distance and the shortest path distance.

Self localisation in WSNs as inference on a graphical model

In this section, the problem of self localisation in a WSN is formulated as inference on a graphical model whose nodes are associated with randomly deployed sensors. Only few nodes of the network, called anchors, have a prior information about their positions; these sensors are placed manually at well known locations or equipped with a GPS module and randomly dropped, along with other sensors, if manual placement is impractical. Each sensor is able to communicate with nearby nodes and to estimate the distance separating them. The purpose is to form a global map of the nodes by combining the information available at the anchors, e.g. obtained from a global positioning system, with measures of relative distances, by exchanging local messages between neighbouring sensors using the box-BP scheme introduced in chapter 2. After convergence of this iterative procedure, each sensor is left with an estimation of its 2D position. This formulation (using a graphical model) is compliant with the distributed aspect of the problem we are considering since it allows to exploit the local independence conditional relationships between the sensor nodes. The use of a message-passing BP scheme is also advantageous since BP is a distributed algorithm which can be employed to infer on a graphical model even when the whole structure of the graph is not known a priori. In the following, the performance of the box-BP will be also compared to that of the standard particle-based NBP and the MDS-MAP(C) algorithm reported in the previous section.

Problem formulation

A number n of sensors is randomly deployed in a planar region. Each sensor computes noisy measurements of its distances from its neighbouring counterparts, the estimated distance is supposed symmetric for both communicating nodes. Two sensors are able to communicate if the distance separating them is less than some maximum range denoted as R, in this case they are said to be neighbours. The position of a sensor t is denoted as x t , t = 1, ..., n, and only a small number of anchor nodes have significant a priori information, p t (x t ), about their positions.

Let d ts denote the noisy measurement of the distance between sensors s and t, then:

d ts = ||x t -x s || + v ts , v ts ∼ p v (x t , x s ), ( 3
P (o ts = 1) = 1 ||xt-xs||≤R . ( 3.2)
The goal is, for each sensor node t of the network, to compute its marginal distribution of interest p(x t |X\x t , {d ts } (t,s):ots=1 ), where X = {x 1 , . . . , x n } is a vector grouping the positions of all the sensors in the network.

Looking back at equation (3.1), for the joint probability distribution p(X, {d ts } (t,s):ots=1 ), the quantities d ts only depend on the random variables x t and x s , and thus, using the chain rule, this joint PDF can be factorised as follows:

p(x 1 , . . . , x n , {d ts } (t,s):ots=1 ) = (t,s):ots=1 p(d ts |x t , x s ) t p t (x t ). ( 3.3)
The undirected graph describing this joint probability distribution is defined by the following two sets: the set V of nodes associated with each sensor and the set E of edges defined as E = {(t, s), o ts = 1}. Each edge thus models an established neighbourhood between two sensors. This factorisation indeed justifies the intuitive representation of this system by a graph wherein each edge represents an actual physical link established between neighbouring sensors/nodes allowing them to communicate and to mutually detect each other.

Implemented solutions

Based on equation (3.3), the set of potential functions associated with the sensor network, i.e. the corresponding graphical model, can be divided into two subsets: the set of pairwise potential functions ψ ts (x t , x s ) defined over two neighbouring/communicating nodes and associated with the set E of edges of the graph, and the set ψ t (x t ) of local likelihoods at each sensor, where:

ψ ts (x t , x s ) = p(d ts |x t , x s ) = p v (d ts -||x t -x s ||), (t, s) ∈ E, ψ t (x t ) = p t (x t ), t ∈ V, ( 3.4)
then, equation (3.3) can be rewritten as:

p(x 1 , . . . , x n , {d ts } (t,s):ots=1 ) = t ψ t (x t ) t,s ψ ts (x t , x s ). (3.5)
Both algorithms, the novel box-BP with interval-based representations of probability quantities and the standard particle-based NBP, are tested for the scenario described here above. For the novel box-BP algorithm, message passing between two neighbouring nodes, namely t and s, is conducted as follows: for the first iteration, messages are initialised using the local observations, only anchor nodes are thus able to send informative messages to their neighbours, since they are the only nodes with significant a priori information, p t (x t ), about their positions. The prior information being represented by a box, N boxes [x j t ] 0 , j = 1, . . . , N, mutually disjoint and weighted equally with of N weighted box-particles, obtained using Algorithm 1 (see section 2.4.1). These boxes are contracted and re-weighted using the local evidence at t as described in Algorithm 2 in the same section. Note that in this application, we are in the special case of an observation y t of the same nature as the state x t , the contraction reduces into an intersection between [y t ] and [x t ] i . The forwarded message from t to s, separated by noisy distance d ts , is then computed by propagating these boxes {[x t ] (i) } through the following model, as shown in the approximation (2.22):

(i) t , [x t ] (i) }
x (i) ts = x (i) t + d ts cos(θ) sin(θ) , ( 3.6)
where θ ∈ [0, 2π]. Algorithm 4 recapitulates BP's message update operation for each sensor t ∈ V . Figure 3.1 illustrates the result of propagating a box through the model (3.6). Let

[x t ] = [x t x t ] × [y t y t ] and [d ts ] = [d d]. Then: [d ts ] × cos[θ] = [d d] × [-1 1] = [-d d], [x t ] + [d ts ] × cos[θ] = [x t x t ] + [-d d] = [x t -d x t + d].
Similarly,

[y t ] + [d ts ] × sin[θ] = [y t -d y t + d].
The implementation of NBP for self localisation in sensor networks is described next:

for the first iteration, N particles {ω 0 t , x 0 t } are sampled from p t (x t ); at an iteration i of NBP, the message product at node t is represented by a collection {ω i t , x i t } of weighted particles. The weights are corrected using the local evidence at t. The forwarded message from t to s, separated by noisy distance d ts , is then computed by propagating the samples {x i t } through the model (3.6) where θ ∼ U ([0, 2π]). In our implementation of the NBP, we used the NBP toolbox available from the work of Leonid Sigal [START_REF]Continuous State Graphical Models for Object Localization, Pose Estimation and Tracking[END_REF].

Simulation results

In this section, we present some simulation results. Both box-BP and NBP algorithms are tested on simulated scenarios and compared to the MDS-MAP(C) described in section 3.3. In our simulations, 100 sensors are randomly deployed in a planar region L × L. We study the performance of the auto localisation algorithms with two different configurations of the anchors: a grid-like and a random placement, we also study the effect of the number of the anchors on the resulting estimation accuracy.

Grid-like placement of the anchors

Basically, two parameters can have impact on the solution provided by a self localisation algorithm: the range within which the sensors can communicate and the number of anchor nodes. Varying the communication range R leads to a variation of the number of neighbourhoods established. Thus, by increasing R, more links between nearby nodes are established and a more complex and dense graph is obtained. In [START_REF] Ihler | Nonparametric belief propagation for self-localization of sensor networks[END_REF], only three anchor nodes are used for calibration (since three is the minimal required number of anchors to be able to determine the global coordinates of the sensors in a two dimensional space), in order to obtain accurate results for NBP, the authors chose to increase the number of neighbourhoods/links established between the sensors. This resulted in a very dense graphical model and thus complex computations and more energy and time consumption. Drawing on scenarios presented in [START_REF] Essoloh | Localisation distribuée dans les réseaux de capteurs sans fil par résolution d'un problème quadratique[END_REF][START_REF] Langendoen | Distributed localization in wireless sensor networks: a quantitative comparison[END_REF], we choose to place nine anchors in a grid-like position, eight of them are located on the contour of the region, these can be placed manually, and one anchor is at the centre of the network. In fact, positioning algorithms perform better when anchors surround the nodes with unknown positions [START_REF] Langendoen | Distributed localization in wireless sensor networks: a quantitative comparison[END_REF]; intuitively, nodes at the edges of the graph are less likely to be connected making their localisation more difficult. The range of communication is set to R = L/4, the noise standard deviation is set to 0.005L. Figure 3.2 shows the distribution of the sensors where anchor nodes are marked by circles, while figure 3.3 shows the corresponding graphical model for R = L/4. Table 3.1 summarises the results obtained for both NBP and box-BP algorithms. The simulations were carried out using Matlab on Intel Core i7-3520M processor (2.90GHz-4MB Cache, Dual-core) for L = 100m. The error refers to the root mean squared error (RMSE), it shows the mean distance between the true and the estimated position of a sensor and is given as a percentage of L. For the box-BP algorithm, the error is calculated using the distance between the centre of the boxes and the true position of the sensors. Note that all three algorithms achieve comparable accuracy (see table 3.1). However, in order to represent the posterior PDF, NBP stores 200 weighted particles, which corresponds, in our 2D application to a total of 600 floating points values. The box-BP algorithm uses only 9 box particles, that is equivalently 45 floating points values. This reduced storage capability is a great enhancement in terms of energy saving and bandwidth needed for the information exchanged in the network. Furthermore, the box-PF is about 10 times faster than NBP (see table 3.1). Regrading the centralised MDS-MAP(C) algorithm, the computational time is naturally less than that of the iterative BP procedures, however the centralised aspect of the MDS-MAP(C) does not comply with the network energy constraint. It also requires the knowledge of the whole network structure (estimating the distance between each pair of nodes) without offering an improved accuracy. Figures 3.4

Random placement of the anchors

We now test both message-passing algorithms on different anchor layouts and compare their performance to the MDS-MAP(C) procedure. Rather than placing the anchors in a grid-like position so as to enclose the randomly spread sensors, the anchors are placed randomly as shown in figure 3.7. We also study the performance of the localisation algorithms while varying the number of anchors within the range 6 to 9. The results are grouped in table 3.2 and table 3.3 for the box-BP and the NBP respectively. Table 3.2 shows that the novel box-BP algorithm presents comparable results in term of the root mean squared error for different number of anchors. However, this error rate is achieved within less time, i.e. less iterations, if the number of anchors is increased. The No only condition is for the anchors to be uniformly deployed across the region. Figure 3.8 shows the results obtained using the box-BP algorithm in the case of a layout with 6 anchor nodes. By examining table 3.3, the NBP turns out to be more sensitive to the anchors layout. Figure 3.9 shows the results obtained using the NBP algorithm for a 6-anchor layout. It can be seen that the sensors on the edge of the region could not be well positioned, the error on the position of some of these reached around L/3 (which justifies the large value of the root mean squared error noticed in table 3.3), this is observed for all layouts (defined by different numbers of anchor nodes).

Finally, by examining table 3.4 the MDS-MAP(C) centralised scheme offers an accuracy that is comparable to that provided by the distributed box-BP algorithm, again with practically no difference between different anchor layouts. To sum up, the novel box-BP algorithm offers an accuracy that is comparable to that of the centralised MDS-MAP(C) without requiring the whole structure of the network and using local information exchange No of anchors in a distributed scheme, which is compliant with the energy constraint in a WSN. The box-BP also offers a better performance than that of the standard NBP while requiring less memory storage, bandwidth and computational time as was discussed previously.

Conclusion

In this chapter, we tested the box-BP algorithm, that uses interval representations of probability quantities, on the problem of self localisation in wireless sensor networks. The simulation results showed that for a grid-like placement of the anchors the estimation accuracy provided by NBP algorithm is achieved by the novel box-BP using much less particles and within less computational time. For a random placement of the anchors, the NBP algorithm failed to accurately locate sensors on the edge of the network. The box-BP provided in this case more accurate results. The advantages offered by the new algorithm are reducing the required storage memory, bandwidth and energy needed to exchange information between nodes of a network. Another advantage of the box-BP over the standard NBP is decreasing the complexity of the computations and the mean computational time. Both BP schemes were also compared to the standard centralised MDS-MAP(C) localisation algorithm. The accuracy offered by the novel box-BP is comparable to that of the centralised approach without the necessity to knowing the whole structure of the network and within a distributed context which complies the energy constraint in a WSN. Nevertheless, in this chapter, the performance of the box-BP was tested on a static wireless sensor network and over simulated data. In the next chapter, we consider a different application involving a dynamic network and test the novel message-passing algorithm on the application of continuously localising a mobile object over simulated and real datasets.

LOCALISATION OF A MOBILE OBJECT IN A WIRELESS SENSOR NETWORK

Introduction

Wireless sensor networks (WSNs) built up with cheap and simple sensor nodes have become a popular solution in a range of areas involving monitoring and decision making [START_REF] Li | Detection, classification, and tracking of targets[END_REF].

From environmental monitoring [START_REF] Lushi | An inverse Gaussian plume approach for estimating atmospheric pollutant emissions from multiple point sources[END_REF] to surveillance [START_REF] Weimer | An approach to leak detection using wireless sensor networks at carbon sequestration sites[END_REF], WSNs applications also include target detection, localisation and tracking [START_REF] Li | Detection, classification, and tracking of targets[END_REF][START_REF] Schiff | Nonparametric belief propagation for distributed tracking of robot networks with noisy inter-distance measurements[END_REF][START_REF] Mao | Multiple objects device-free passive tracking using wireless sensor networks[END_REF][START_REF] Mourad | Anchor-based localization via interval analysis for mobile ad-hoc sensor networks[END_REF]. Many applications need to locate people or objects, e.g., location aware end-users, intrusion detection and emergency services for finding persons in need, especially in rescue operations based on mobile phones [START_REF] Sayed | Network-based wireless location: challenges faced in developing techniques for accurate wireless location information[END_REF]. Localisation algorithms can be classified into different categories as was discussed in the previous chapter. In this chapter, we consider a range-based algorithm for mobility tracking of a single target using information collected from a network of range sensors. At each time step, the target can be detected only by a small number of time-of-flight (or range) sensors, and estimates of the distances separating them from the target are calculated. Thus, at each time instant, the position of the target can be evaluated knowing only the positions of the sensors within the vicinity of it. Since usually even the sensors are not perfectly localised, i.e. their positions are known up to some error, the sensors' coordinates are described using random variables, and it is desired to improve the accuracy on these positions while performing the common task of mobility tracking. In order to exploit the local independence conditional relationships, we will formulate the localisation problem as inference on a dynamic graphical model and use a BP scheme to continuously localise the mobile object. Using the BP, we will also be able to reduce the error on the sensors' positions.

In the following, we provide an overview of the prior related work in section 4.2 and then formulate our problem as inference on a graphical model in section 4.3. The box-BP will be employed to perform the mobility tracking task and its performance will be compared to that of the standard NBP. The effectiveness of the novel interval-based message-passing procedure will be demonstrated on simulated and real datasets. Similarly to what was observed in the previous chapter (over the static network), the box-BP will offer an impressive reduction of the memory, bandwidth and computational time required to achieve the localisation goal. Conclusions about this chapter are drawn at last in section 4.4.

Literature review

Location awareness is fundamental in many applications of WSNs. Localisation and object tracking in sensor networks have become important research problems and have been widely studied. For instance the purpose of the work of Ermis et al [START_REF] Ermis | Detection and localization in sensor networks using distributed fdr[END_REF][START_REF] Ermis | Robust distributed detection with limited range sensors[END_REF][START_REF] Ermis | Adaptive statistical sampling methods for decentralized estimation and detection of localized phenomena[END_REF] is to detect localised distributed events or sources in a WSN while minimising inter-network communications cost. The problem of detecting and localising a source emitting a signal with unknown power is considered in [START_REF] Ermis | Detection and localization in sensor networks using distributed fdr[END_REF], whilst [START_REF] Ermis | Robust distributed detection with limited range sensors[END_REF] develops a multi-target detection al-gorithm and [START_REF] Li | Detection, classification, and tracking of targets[END_REF][START_REF] Mao | Multiple objects device-free passive tracking using wireless sensor networks[END_REF][START_REF] Sheng | Distributed particle filter with gmm approximation for multiple targets localization and tracking in wireless sensor network[END_REF] present multiple target tracking results. In [START_REF] Li | Detection, classification, and tracking of targets[END_REF], multiple objects moving in a sensor field are tracked by fulfilling several tasks in a collaborative signal processing framework. Decision making requires target detection, classification, localisation and future position prediction. The localisation process employs commonly data provided by a number of sensors such as acoustic, seismic or thermal sensing modalities and is based on the time-varying space-time signature an object generates when moving in a geographical region. In [START_REF] Mao | Multiple objects device-free passive tracking using wireless sensor networks[END_REF] light sensors (along with light sources) are used to track multiple moving objects. A distributed particle filter (DPF) is proposed in [START_REF] Sheng | Distributed particle filter with gmm approximation for multiple targets localization and tracking in wireless sensor network[END_REF] to track multiple targets in a WSN. In this approach, the sensors are first clustered into disjoint and uncorrelated cliques. Next, unlike a centralised particle filter (CPF), a DPF updates the weights of the particles, resulting from the prediction step of a particle filter [START_REF] Johansen | A Tutorial on Particle Filtering and Smoothing: Fifteen years later[END_REF], sequentially or in parallel, using only local measurements (i.e. observations within a clique of sensors).

In [START_REF] Baggio | Monte carlo localization for mobile sensor networks[END_REF][START_REF] Dil | Range-based localization in mobile sensor networks[END_REF] sequential Monte Carlo (SMC) localisation is presented that improves the approach introduced in [START_REF] Hu | Localization for mobile sensor networks[END_REF]. The SMC algorithm described in [START_REF] Hu | Localization for mobile sensor networks[END_REF] considers a scenario with mobile nodes and anchors, in an environment where no range measurements can be obtained and under the only assumption that all nodes and anchors have a known maximum speed and the same radio range. In [START_REF] Dil | Range-based localization in mobile sensor networks[END_REF], range measurements are included into the recursive SMC computations. Other SMC algorithms are developed in [START_REF] Mihaylova | Localization of mobile nodes in wireless networks with correlated in time measurement noise[END_REF][START_REF] Mihaylova | Sequential Monte Carlo methods for localization in wireless networks[END_REF]. In contrast to [6], the algorithm introduced in [8] remains range-free and works on constricting the area from which the SMC procedure draws samples. In [START_REF] Mourad | Anchor-based localization via interval analysis for mobile ad-hoc sensor networks[END_REF], an algorithm is proposed in the bounded error framework for self-localisation of a mobile node in an ad-hoc sensor network, while in [START_REF] Schiff | Nonparametric belief propagation for distributed tracking of robot networks with noisy inter-distance measurements[END_REF] the purpose is to solve the problem of tracking multiple robots by formulating it as inference on a graphical model and by using a variant of Belief Propagation (BP) algorithm, called Non-parametric Belief Propagation (NBP) [START_REF] Sudderth | Nonparametric belief propagation[END_REF]. While providing more accurate location estimations, the NBP requires more computations. This observation was the motivation behind using the box-BP approach we introduced in chapter 2 to localise an object moving within a WSN. The objectives are to achieve the accuracy provided by the NBP within less computational time. A similar approach is given in [START_REF] Zhexu | An efficient message passing algorithm for multi-target tracking[END_REF] where a message-passing algorithm is used to solve the problem of tracking multiple targets. This algorithm propagates information from future data about past hypotheses using particle representations of the messages.

The scenario we consider is similar to what was proposed in [START_REF] Sayed | Network-based wireless location: challenges faced in developing techniques for accurate wireless location information[END_REF][START_REF] Thomas | Revisiting trilateration for robot localization[END_REF]. The purpose is to locate an object using exclusively its distance measurements to other known points or anchors with no assumptions or limitations on the object movement. We assume that a set of range sensors are deployed in the field of interest and are able to detect and estimate the distance separating them from an object that falls into their sensing range, i.e., their maximal detectable distance. At each time step, the moving object will be in the range of a limited set of sensors in the network. The sensors provide noisy measurements of the distance separating them from the target. Since the dynamics of the moving object is not known, the localisation at each step is based exclusively on exchanging information with the sensors that are able to detect the target. We formulate the problem as inference on a graphical model and employ a variant of Belief Propagation (BP) algorithm, the box-BP which uses box representations for probability quantities, to localise the target over time.

Localisation of a mobile target as inference on a graphical model 4.3.1 Simulated Scenario

In this section, we consider the problem of localisation of a mobile target, moving within a region of interest, in the presence of randomly distributed sensors and data arriving sequentially in time. A number n of range sensors are deployed in the considered region and are assumed to be able to detect and estimate the distance separating them from an object that falls into their sensing range. This range denotes the maximal detectable distance and will be referred to as R. At each time step (whose value will be denoted as τ ), the mobile object will be in the range of a limited set of sensors. Thus, estimates of the positions of the moving object will computed each τ sec using only available measurements of the distances separating it from the limited number of sensor nodes that were able to detect it.

The position of the target at time kτ , k ∈ N, is denoted as x t,k while that of the j th sensor is referred to as x s j , j = 1, . . . , n. The sensors are static, their positions are fixed in time and are assumed to be known a priori (up to some error). The inference purpose in this application is to determine the posterior distributions p(x t,k |{x s j } n j=1 ), k = 1, . . . , T where T refers to the number of time steps.

Implemented Solutions

Let c j,k denote a boolean variable which indicates whether or not the target t is in the range of sensor j at time kτ . Thus,

c j,k =    1 if ||x t,k -x s j || ≤ R, 0 otherwise.
In the following, Γ k designates the set of sensors that are able to detect the target at

d t ,k = ||x t,k -x s || + v t ,k , v t ,k ∼ p v (x t,k , x s ), (4.1)
where v t ,k denotes a statistical measurement noise of a probability distribution p v (x t,k , x s ), and || . || denotes the Euclidean distance between the target t at kτ and s ∈ Γ k . Given the position x s of sensor s , the target will be at position

x t,k = x s + ||x t,k -x s || × cos(θ t ,k ) sin(θ t ,k ) ,
where θ t ,k is the angle between x s and x t,k . Thus, the joint probability distribution p(x t,k , x s 1 . . . , x sn ) can be factorised as follows:

p(x t,k , {x s j } n j=1 ) = p(x t,k |{x s j } n j=1 ) j p(x s j ) = p(x t,k |{x s } s ∈Γ k ) j p(x s j ), (4.2)
where p(x s j ) refers to the prior probability distribution of sensor j's position.

According to (4.2), the directed graph describing this joint probability distribution is a graph whose nodes represent the target at time kτ and the set Γ k of its neighbouring sensors, and whose edges are directed from nodes {s ∈ Γ k } to the node associated with the target.

Figure 4.2 shows the directed graphs representing the joint PDF in (4.2) at two consecutive time steps, k -1 and k. Note that if the dynamics of the target is known, thus a function linking x t,k-1 to x t,k is available, information would also propagate between time steps which translates graphically into an edge directed from node t k-1 towards t k . Another point is that although the arcs of the graph are directed, information circulates in either direction and the operations of BP as described in the previous sections can be used.

In our simulations, we fixed the time period to τ = 1 sec; hence, each second, new links are established between the target and the sensors that, at the considered time instant, have the moving object within their sensing range, while links that were made at the previous time step would be broken when the moving object steps out of their range: the resulting graphical model is dynamic as can be seen in Figure 4.2. In order to solve the inference problem for the simulated scenario described here (which consists in estimating the posteriors p(x t,k |{x s j } n j=1 )), two algorithms, the novel NBP with box representations -referred to as box-BP-and NBP with sample representation -referred to simply as NBP -are compared. During a time period τ , three iterations of messagepassing are run, messages are exchanged along the edges of the graph and the target node computes an estimation of its current position. Recall that the BP is a tool to perform marginalisation of joint PDFs represented by graphical structures, and that it allows to estimate the sensors' positions wile performing the mobility tracking task. However, in our simulations, the purpose is to continuously localise the mobile object and thus the marginal PDFs of interest are given by p(x t,k |{x s j } n j=1

). An example of message-passing computations is given in the following.

Algorithm 5 BP's message update for mobility tracking at time instant kτ

Γ k = {s , c ,k = 1}
Message passing from the sensors to the mobile target for s ∈ Γ k do use the available observation at sensor s to generate uniformly N boxes mutually disjoint and weighted equally with 1

N

propagate the boxes through the model x

(i) s ,t k = x (i) s + d tl,k [cos(θ) sin(θ)]
T end for combine the messages sequentially using the series of intersections in Algorithm 1 normalise the weights the resulting boxes referred to as {[x t,k ] (i) , ω t,k } N i=1 represent the belief at the target at time kτ Message passing from the target to the sensors for s ∈ Γ k do propagate the resulting boxes through the model x (i)

t k ,s = x (i) t,k + d tl,k [cos(θ) sin(θ)]
T combine the messages sequentially using the series of intersections in Algorithm 1 contract the resulting boxes using the observation at sensor s normalise the weights the resulting boxes {[x s ] (i) , ω s } N i=1 now represent the position of sensor s end for

In this example, message-passing between two neighbouring nodes (which will be denoted as 1 and 2 for the simplicity of the notation) is explained. For the novel box-BP algorithm, the prior information at each sensor is represented by a box, N boxes mutually disjoint and weighted equally with 1 N are uniformly generated from the available information about sensors positions. For NBP, N particles are sampled from the prior p(x s ). At an iteration i, in the bounded error context, the message product at node 1 is represented by a collection {ω (which is merely a disk equation):

x (i) 12 = x (i) 1 + d 12 [cos(θ) sin(θ)] T , ( 4.3)
where θ ∈ [0, 2π]. Algorithm 5 describes the message passing procedure to perform the mobility tracking task.

For NBP algorithm, the message product at the sending node 1 being represented by a collection {ω i 1 , x i 1 } of weighted particles, the forwarded message from 1 to 2, separated by noisy distance d 12 , is computed by propagating the samples {x i 1 } through the model (4.3) where θ ∼ U ([0, 2π]).

Furthermore, both probabilistic algorithms are also compared to the optimisation approach introduced by Sayed et al. in [START_REF] Sayed | Network-based wireless location: challenges faced in developing techniques for accurate wireless location information[END_REF] and to which we will be referring as Sayed05. A brief review of the Sayed05 algorithm is given in the next subsection.

The Sayed05 Localisation Algorithm

Sayed05 is an algorithm for 2D landmark-based localisation that was proposed in [START_REF] Sayed | Network-based wireless location: challenges faced in developing techniques for accurate wireless location information[END_REF] in response to a mandate issued by the Federal Communications Commission (FCC) that aimed to solve a public safety problem: the need to determine the geographic coordinates of a mobile subscriber in a cellular network knowing that an important proportion of emergency calls originates from cell phones. This algorithm is based on the knowledge of distances separating the target, i.e. the non-localised object, from three known points or stations also called landmarks. This problem is commonly known as trilateration and is illustrated in Figure 4.4. The Sayed05 algorithm works as follows. Let r t = [x t , y t ] T denote the position of the target and r i = [x i , y i ] T , i = 1, 2, 3 the coordinates of three landmarks. The squared distance between r t and r i is given by Subtracting the first equation (i =1) from the i-th equation (i=2, 3) in (4.4), we obtain

d 2 i = (x i -x t ) 2 + (y i -y t ) 2 , i = 1, 2, 3. (4.4)
d 2 i -d 2 1 = x2 i -2x i xt + ỹ2 i -2ỹ i ỹt , i = 2, 3, (4.5)
where r = [x, ỹ] T = rr 1 . The constraints in (4.5) can be expressed in matrix notation as x2 ỹ2 x3 ỹ3

xt ỹt = 1 2 x2 2 + ỹ2 2 + d 2 1 -d 2 2 x2 3 + ỹ2 3 + d 2 1 -d 2 3 , ( 4.6)
or, in a compact form, Hr t = b. (4.7)

The algorithm then uses a least squares approach to determine rt = bH -1 .

Simulation Results

The simulations are run under Matlab on Intel Core i7-3520M processor (2.90 GHzâĂŞ4 MB Cache, Dual-core). The sensing range is set to R = 15 m. We run two series of simulations: the first series considers a white Gaussian noise, of standard deviation σ, on the measured distances; the second considers a bounded measurement noise modelled as a uniform PDF defined over [-a a], (σ, a ∈ R + * ). The performance of both box-BP and NBP are compared to that of the Sayed05 localisation algorithm that was briefly recapitulated in the previous subsection. We also study the effect of the measurement noise by varying σ and a defined earlier. Tables 4.1 and 4.2 summarize the results obtained for all three algorithms in the case of a Gaussian measurement noise and a bounded noise respectively. The error refers to the root mean squared error, that is the mean distance between the true and the estimated positions. The time refers to the mean computational time per time step. These tables show that the novel box-BP algorithm provides a better accuracy (or a similar accuracy in the case of a Gaussian measurement noise for large σ) than the NBP, the optimisation Sayed05 algorithm while offering similar performance to the probabilistic NBP and box-BP algorithms for low measurement noise, provides less accurate estimates as the noise magnitude increases. An interpretation of these results is based on the fact that the box-BP provides a direct approximation of the posterior distribution using a mixture of uniform PDFs, and not a representation of it using particles as in the case of NBP. To calculate the messages and beliefs, there is no need to introduce kernels, which is not necessarily optimal for the NBP. Also, compared with the NBP, the box-BP and the Sayed05 algorithms are simpler and introduce less approximations. Furthermore, the novel box-BP uses only 9 box particles to calculate the probability quantities, that is equivalently 45 floating points values, versus 200 weighted particles for NBP, which corresponds, in our 2D application to a total of 600 floating points values. This reduced storage requirement is a great enhancement in terms of energy saving and bandwidth needed for the information exchanged in the network. The box-PF is about 30 times faster than the NBP. While the mean computational time for the Sayed05 optimisation algorithm is about 10 times less than that of the box-BP, the performance of this optimisation technique degrades in noisy environments since it does not take into account a measurement error. These results can be seen by examining Figures 4.5

σ = 0.5 m σ = 1 m σ = 1.5 m Algorithm N o of

Results Obtained on Real Data

In this section, the box-BP algorithm is tested on real data available from the UTIAS multi-robot cooperative localisation and mapping dataset [START_REF] Leung | The utias multi-robot cooperative localization and mapping dataset[END_REF]. This data set is produced using 5 robots moving in an indoor workspace with an area of 15 × 8 m 2 and where 15 landmarks were placed as illustrated in Figure 4.11. The landmarks consisted of cylindrical tubes of 30 cm in height. Figure 4.8 shows the cylindrical landmarks and the robots.

The robots are identical in construction, equipped with a monocular camera serving as a primary sensing module and interfacing with a netbook. The camera is mounted in a way that the robot's body frame and the measurement coordinate frame are coincident in two dimensions, as shown in Figure 4.9.

During the experiments the robots were driving to randomly generated waypoints in the workspace. When a robot or landmark falls in a robot's field of view, a range and bearing measurements are taken. Throughout the data collection process, the groundtruth data for the robots and the landmarks were made available through the use of a 10-camera Vicon motion capture. Figure 4.10 illustrates a scene of the experiments conducted while producing the UTIAS dataset.

The original dataset is transformed as in [START_REF] Choi | Triangulation toolbox: Open-source algorithms and benchmarks for landmark-based localization[END_REF] in order to adapt to the scenario we consider in this paper: the localisation of a single moving target using exclusively range measurements from landmarks. We thus extract from the original dataset, data associated with poses, of any of the five robots, for which range measurements to landmarks are available. A table of real groundtruth data for a single target is obtained, along with a table indicating for each pose/step identifiers of the landmarks within the vicinity of the non localised object and a real noisy measurement of the distance separating them. The performance of the box-BP is again compared to that of the NBP and Sayed05 algorithms.

The results are grouped in Table 4.3.

The cumulative distribution function of the error is given in Figure 4.12.

The novel box-BP algorithm outperforms both NBP and Sayed05 algorithms. The performance of the Sayed05 algorithm turns out to degrade in noisy environments; though this optimisation algorithm requires less computations than the remaining two probabilistic algorithms. The advantages of storage, memory and bandwidth reduction offered by the box-BP over the NBP remain outstanding as with the simulated scenarios.

Conclusion

In this chapter, we apply the novel box-BP to perform the localisation of a moving target, at each time step, using only range measurements from landmarks. The effectiveness of the new box-BP algorithm is shown on simulated and on real data. Its performance is also compared to that of a probabilistic and an optimisation algorithms, the NBP and the Sayed05 algorithms, respectively. The simulation results shows that the estimation accuracy provided by the novel box-BP algorithm is better (or similar) to that of the NBP, and is achieved using much less particles and within less computational time. The performance of the novel algorithm also remains outstanding even in noisy environments unlike the fast optimisation Sayed05 algorithm. The accuracy of the box-BP turns out to be better on the real data set. In summary, the advantages offered by the new algorithm are a reduction of the required storage memory, bandwidth and energy needed to exchange information between the nodes as compared to the NBP and a better accuracy in noisy scenarios when compared to the Sayed05 algorithm. The box-BP in fact introduces less approximations than the NBP (which uses kernels to approximate the PDFs) and thus offer a better accuracy, it is also applicable in the presence of quantised measurements, this would require less bandwidth for information exchange since float type measurements are often replaced by integers.

LITERATURE REVIEW ON SOURCE CHARACTERISATION

AND ADVECTION-DIFFUSION MODELS

Introduction

Recently, wireless sensor networks (WSNs) are being successfully used in source identification applications, which are approaches for automatic detection and characterisation of sources of dispersion (i.e. sources releasing a chemical or biological agents into their environment) based on data collected from a WSN. In this thesis, we are interested in the application of land mine detection and localisation given the humanitarian aspect of this topic, and we refer, throughout the following chapters, to the task of determining the land mines/ sources unknown parameters, including their positions and emission rates, as the problem of source characterisation. An anti-personnel mine is a device designed to injure or to kill whomever comes into contact with it through direct pressure or a trip-wire [START_REF] Habib | Mine detection and sensing technologies-new development potentials in the context of humanitarian demining[END_REF]. Land mine detection, localisation and clearance is of great importance due to the danger that buried land mines still represent to people all over the world. It is pointed out in [START_REF] Robledo | A survey of land mine detection technology[END_REF] that about 50-100 million anti-personnel mines exist in more than 80 countries and that more than 20,000 people are killed or injured every year due to the explosion of buried land mines. Though the Ottawa treaty prohibited in 1997 the use of this weapon, some countries have not signed the treaty, and nearly two million land mines are laid per year. The dangers are particularly acute for those responsible for localising and decommissioning land mines. To limit the number of victims, land mine detection and clearance actions have taken place since the end of the Second World War [START_REF] Cardona | Landmine detection technologies to face the demining problem in antioquia[END_REF][START_REF] Ivashov | A review of the remote sensing laboratory's techniques for humanitarian demining[END_REF].

Historically, the most common method for land mine detection is metal detection. Although this technique has proved to be effective with primary land mines, it fails with many modern land mines, which are fabricated from sophisticated non-metallic materials, such as plastic and wood [START_REF] Habib | Mine detection and sensing technologies-new development potentials in the context of humanitarian demining[END_REF], making them invisible to the metal detector. Therefore, many other methods have been developed. These include the use of trained dogs and several physical detection techniques based on ground penetrating radar (GPR), X-ray, infrared (IR) imaging [START_REF] Lopez | Antipersonnel mine detection on infrared images[END_REF][START_REF] Paik | Image processing-based mine detection techniques: A review[END_REF], neutron activation (TNA) and nuclear quadrupole resonance (NQR) [START_REF] Bruschini | A survey of current sensor technology research for the detection of landmines[END_REF][START_REF] Robledo | A survey of land mine detection technology[END_REF]. However, a common problem with all these techniques is that the probability of false positives is high [START_REF] García-González | Sensors: From biosensors to the electronic nose[END_REF]. Other approaches employ unmanned vehicles for landmine detection [START_REF] Gooneratne | A review of sensing technologies for landmine detection: Unmanned vehicle based approach[END_REF]. This technique requires sophisticated and rather expensive equipment and control.

In this second part of the thesis, we address the problem of localising an unknown number of land mines using data collected from a network of wireless sensors capable of detecting the concentration of the explosive chemicals in the air. The motivation behind the proposed framework is that buried land mines can be considered as vapour-emitting sources, based on the observation that some explosive chemicals, such as trinitrotoluene (TNT) or dinitrotoluene (DNT), leak out from buried land mines into the surrounding environment and are transported through the air by mechanisms such as advection and diffusion; the idea of using a WSN for land mine localisation then seems intuitively feasible. the idea is to deploy a network of wireless sensors, capable of sensing the concentration of the explosive chemical at their positions, in the contaminated region. Sensors for these types of chemical explosive materials exist [START_REF] Yu | Ionic liquid thin layer EQCM explosives sensor[END_REF], thus, having a vector of concentration measurements provided by a WSN, by expressing the concentration of the explosive as a function of the land mines' locations and by solving the inverse problem, we will show in the following chapters that the proposed framework is able to detect, locate and find the emission rates of several land mines.

One primary advantage of the proposed technique consists in the fact that existing land mine detection and localisation techniques (e.g. those reported here above) involve either human (or animals) intervention or the use of bulky and expensive equipments. This translates into threatening risks on the operators lives and/or high operational costs. Using a network of low-cost small-sized wireless sensors allows an automatic and less expensive detection of buried land mines. Furthermore, it was suggested in [START_REF] Dias | Investigating The Viability Of Mems Vapor Sensors For Detecting Land Mines[END_REF] that using this technique reduces false positives since the sensors respond to the key element of a landmine, the explosive chemical.

In this chapter, we first provide a literature review on the use of sensor networks for source characterisation in section 5.2. Then, an overview of advection-diffusion models in the cases of instantaneous and continuous release point sources is available in section 5.3. The forward model, aiming at predicting the concentration measurements at the sensors' positions given the sources parameters, is formulated in section 5.4. Section 5.5 then concludes the chapter.

Related works

The advances in sensing technologies [START_REF] García-González | Sensors: From biosensors to the electronic nose[END_REF] increase the use of sensor networks in a vast range of applications [START_REF] Badine | Considerations for implementing a wireless network for landmine detection[END_REF][START_REF] Senocak | Stochastic event reconstruction of atmospheric contaminant dispersion using Bayesian inference[END_REF][START_REF] Lu | Space encoding based compressive multiple human tracking with distributed binary pyroelectric infrared sensor networks[END_REF]. Recently, wireless sensor networks (WSNs) have become popular in source identification applications. In fact, risk management applications in the fields of environment [START_REF] Borysiewicz | Stochastic algorithm for estimation of the model's unknown parameters via Bayesian inference[END_REF][START_REF] Keats | Bayesian inference for source determination in the atmospheric environment[END_REF][START_REF] Lushi | An inverse Gaussian plume approach for estimating atmospheric pollutant emissions from multiple point sources[END_REF] and security [START_REF] Senocak | Stochastic event reconstruction of atmospheric contaminant dispersion using Bayesian inference[END_REF][START_REF] Weimer | An approach to leak detection using wireless sensor networks at carbon sequestration sites[END_REF] rely on data collected from a WSN in order to characterise a source of dispersion, e.g., in the case of an accidental or intentional release of a chemical or biological substance in the air. In [START_REF] Weimer | An approach to leak detection using wireless sensor networks at carbon sequestration sites[END_REF], an algorithm is derived to detect CO2 leaks at several potential locations at a carbon sequestration site. The aim in [START_REF] Lushi | An inverse Gaussian plume approach for estimating atmospheric pollutant emissions from multiple point sources[END_REF] is to study the emissions of a number of contaminant sources, located at well-known positions, at a large lead-zinc smelter. In [START_REF] Borysiewicz | Stochastic algorithm for estimation of the model's unknown parameters via Bayesian inference[END_REF][START_REF] Keats | Bayesian inference for source determination in the atmospheric environment[END_REF][START_REF] Senocak | Stochastic event reconstruction of atmospheric contaminant dispersion using Bayesian inference[END_REF], probabilistic Bayesian approaches are used to determine the unknown position and possibly other model parameters, such as the emission rate and the diffusion coefficient, of a single dispersion source using data collected from a WSN. In [START_REF] Ram | Localization and intensity tracking of diffusing point sources using sensor networks[END_REF], a recursive algorithm based on a state space representation of the system is developed to estimate a single diffusion source position and to track its intensity in time using concentration measurements provided by a sensor network. In [START_REF] Ristic | Cramer-rao bound for source estimation using a network of binary sensors[END_REF] theoretical results are derived characterising the accuracy of the location estimate of a single gas emitting source using a network of binary sensors. The measurements are quantised and a single bit of information is generated depending on whether the sensed value is lower or higher than some threshold. In [START_REF] Ristic | Bayesian likelihood-free localisation of a biochemical source using multiple dispersion models[END_REF] a computation method is proposed to overcome the difficulties associated with the choice of an adequate dispersion model and the calculation of the likelihood function in a Bayesian framework in order to solve the problem of localisation of a source of toxic release.

For landmine detection, the idea is to spread wireless sensors over the contaminated field, e.g. using a flying robot, these sensors measure the concentration of the explosive material at their positions and communicate with a processing unit. The measured concentrations depend on the land mines and sensors positions and on the perturbation caused by multiple environmental effects, e.g. the wind. The complexity of the problem is seen through the following reasons: a. the sensors positions are unknown since they are randomly dropped. This problem was however dealt with in chapter 3, b. the mines are buried and diffusion constants are different in the ground and in the air, c. diffusion constants are not known accurately, d. environmental disturbance is hardly modelled, e. an error exists between the real and the measured concentrations.

The idea of using a sensor network for land mine localisation is addressed in [START_REF] Dias | Investigating The Viability Of Mems Vapor Sensors For Detecting Land Mines[END_REF][START_REF] Jeremić | Landmine detection and localization using chemical sensor array processing[END_REF], but both consider the case of a single land mine. In [START_REF] Dias | Investigating The Viability Of Mems Vapor Sensors For Detecting Land Mines[END_REF], the problem of localisation of a single land mine is considered using an analytical solution of the inverse problem, not taking into account a model or measurement noise. In [START_REF] Jeremić | Landmine detection and localization using chemical sensor array processing[END_REF], a maximum likelihood estimation algorithm is derived in order to locate a single land mine and find its emission rate. The performance of the estimator is evaluated by computing the Cramer-Rao bound.

The case of parameter estimation for multiple sources is briefly addressed in [START_REF] Yee | Bayesian probabilistic approach for inverse source determination from limited and noisy chemical or biological sensor concentration measurements[END_REF] and [START_REF] Sharan | Least square data assimilation for identification of the point source emissions[END_REF] in a Bayesian probabilistic and an optimisation least squares frameworks and for a known number of sources. However, in the context of land mine detection, this information rarely exists for real applications. A more difficult case when the number of sources is unknown is addressed in [START_REF] Yee | Theory for reconstruction of an unknown number of contaminant sources using probabilistic inference[END_REF][START_REF] Yee | Inverse dispersion for an unknown number of sources: Model selection and uncertainty analysis[END_REF][START_REF] Yee | Probability theory as logic: data assimilation for multiple source reconstruction[END_REF]. In [START_REF] Yee | Theory for reconstruction of an unknown number of contaminant sources using probabilistic inference[END_REF][START_REF] Yee | Probability theory as logic: data assimilation for multiple source reconstruction[END_REF], the problem is formulated as a generalised parameter estimation problem, where the number of sources is included in the vector of unknown parameters. This approach implements a reversible-jump MCMC algorithm and requires intensive computations since the dimensionality of the unknown parameters' vector is variable. In [START_REF] Yee | Inverse dispersion for an unknown number of sources: Model selection and uncertainty analysis[END_REF], Yee formulates the problem of characterising an unknown number of sources as a model selection problem. While this approach is less complex than the previous one, it is also computationally demanding.

The application of interest in this chapter is the localisation of multiple anti-personnel land mines using a WSN, where we consider the case of an unknown number of sources. In our method, the objective is to estimate the number of sources and then characterise them. Briefly speaking, the set of concentration measurements which have been made by the detection system are grouped in a matrix and a PCA scheme is used in order to determine the number of sources. The use of this strategy makes the problem of localisation less complex and more efficient. The localisation problem is addressed after in a probabilistic Bayesian framework and a Markov chain Monte Carlo (MCMC) algorithm, namely, that of slice sampling, is used in order to sample from the posterior density of interest. This probabilistic approach is tested and compared to an optimisation technique, i.e., the popular least squares approach [START_REF] Tarantola | Inverse Problem Theory and Methods for Model Parameter Estimation[END_REF]. The advantages and limitations of both techniques are discussed in detail. There are two main advantages of using a probabilistic approach. First, the solution provided takes on the form of a probability distribution, so the uncertainty on the estimated position can be quantified [START_REF] Keats | Bayesian inference for source determination in the atmospheric environment[END_REF], rather than approximated, e.g., by computing the Cramer-Rao bound, as in [START_REF] Jeremić | Landmine detection and localization using chemical sensor array processing[END_REF]. Another important advantage of the proposed Bayesian technique is that it overcomes the convergence problems (to local minima) that the least squares approach could face.

Advection-diffusion models

Processes that move a compound from the location at which it is generated through the air (or water) are called mass transport processes and can be divided into two categories: diffusion and advection. Diffusion works to eliminate discontinuities in concentration profiles by transporting a compound from a region where its concentration is high to a region of a lower concentration through the action of random motions. A simple example is that of a cigarette smoke which diffuses into the air and spreads throughout the room or a drop of ink which, once dropped in a glass of water, diffuses into the whole glass. Advection on the other hand refers to the transport of a compound due to the mean fluid velocity or flow. For instance, a wind blowing to the north carries a compound present in the atmosphere towards the north by advection.

In the following, we start by considering the case of a single emitting source located at (x s , y s , z s ) T in an unbounded domain (x, y, z). We deal first with the case of an instantaneous point release source: we derive the expression of the concentration C(r, t) at position r = (x, y, z) T at time t with an assumption that the emitted compound is exclusively transported by diffusion (i.e. the advection time scale is slow compared to the diffusive time scale); then using the previous result we generalise to the case of a transport by advection and diffusion. Second, we derive a concentration profile by solving the advection-diffusion equation in the case of a continuous release source emitting a substance at a constant rate.

Instantaneous release source

Transport by diffusion

Consider that an instantaneous release of a total mass m [Kg] occurs at time t 0 at point r s = (x s , y s , z s ) T . The emitted compound molecules are diffused in the atmosphere, the variation in time and space of the compound concentration C(r, t) [Kg/m 3 ] is governed by the equation of mass conservation [START_REF] Kathirgamanathan | Source term estimation of pollution from an instantaneous point source[END_REF]:

∂C ∂t = -∇q, (5.1)
where q = -K.∇C [Kg/s/m 2 ] (

is the mass flux per unit area and K is a tensor grouping the diffusivities (m 2 /s) in the x, y and z directions and is assumed to be diagonal, i.e. K = diag(K x , K y , K z ). Under the assumption

K x = K y = K z = K, equation (5.1) becomes ∂C ∂t = K ∂ 2 C ∂x 2 + ∂ 2 C ∂y 2 + ∂ 2 C ∂z 2 , ( 5.3)
and is called the diffusion equation.

Let us begin with a one-dimensional system in x, e.g. the total amount of substance m is diffusing in a cylinder of an infinite length and a cross-sectional area A yz [m 2 ]. The diffusion equation reduces to

∂C ∂t = K ∂ 2 C ∂x 2 . ( 5.4)
The boundary conditions for an unbounded domain are

C(x, t), ∂C ∂x -→ 0 if x -→ ±∞,
while the initial conditions are

C(r, t 0 ) = C(x, t 0 ) = m A yz δ(x -x s ).
The solution of equation ( 5.4) has the form [START_REF] Crank | The mathematics of diffusion[END_REF]:

C(x, t) = B (t -t 0 ) (1/2) exp - (x -x s ) 2 4K(t -t 0 ) , ( 5.5)
where B is a constant to be determined using the following mass conservation condition:

+∞ -∞ C(x, t)dx = m.
Setting λ 2 = x 2 /4K(t -t 0 ) and integrating with respect to λ we can find:

C(x, t) = m/A yz (4πK(t -t 0 )) (1/2) exp - (x -x s ) 2 4K(t -t 0 ) . ( 5.6)
In the case of a three-dimensional system, the initial and the boundary conditions for an unbounded domain become, respectively:

C(r, t 0 ) = C(x, y, z, t 0 ) = mδ(x -x s , y -y s , z -z s ), C(r, t), ∂C ∂x , ∂C ∂y , ∂C ∂z -→ 0 if x, y, z -→ ±∞.
The solution to equation (5.3) can be found using the variable separation technique, thus, we try to determine a solution having the form

C(x, y, z, t) = C 1 (x, t)C 2 (y, t)C 3 (z, t) which satisfies the condition of mass conservation +∞ -∞ +∞ -∞ +∞ -∞ C(x, y, z, t)dxdydz = m.
This expression can be easily deduced from expression (5.6) and is given by [START_REF] Crank | The mathematics of diffusion[END_REF][START_REF] Kathirgamanathan | Source term estimation of pollution from an instantaneous point source[END_REF]:

C(x, y, z, t) = m (4πK(t -t 0 )) (3/2) exp - (x -x s ) 2 + (y -y s ) 2 + (z -z s ) 2 4K(t -t 0 ) . (5.7)
Note that the concentration is in fact expressed in terms of the distance ||rr s || between location r and the source. Figure 5.1 illustrates the variation in space and in time of the concentration of the emitted chemical in function of the distance ||rr s ||, separating the point where the concentration is estimated from the source, according to equation (5.7). In this example, the parameters are fixed to t 0 = 0, K = 25 m 2 /s and m = 1 mg.

Transport by advection-diffusion

In this section, we consider the combined transport of the emitted substance by advection and diffusion for a single instantaneous point release. Considering a mean wind velocity v = (V x , V y , V z ), the mass flux per unit area is now given by: q = Cv -K.∇C.

(5.8)

Under the same assumptions made in the previous section, the advection-diffusion equation is deduced from condition (5.1) and is written as:

∂C ∂t + V x ∂C ∂x + V y ∂C ∂y + V z ∂C ∂z = K ∂ 2 C ∂x 2 + ∂ 2 C ∂y 2 + ∂ 2 C ∂z 2 .
(5.9)

Again, we begin with a one-dimensional system as described in the previous section, a total mass release m occurs at time t 0 at location x s . The advection-diffusion equation reduces to

∂C ∂t + V x ∂C ∂x = K ∂ 2 C ∂x 2 , ( 5.10)
with the same initial and boundary conditions for an unbounded domain.

To solve equation (5.10), let us change the frame of reference from the stationary coordinate system (x, t) to the moving frame (θ, t) where θ = x -x s -V x (t -t 0 ) which moves with a velocity equal to that of the mean wind flow V x (see figure 5.2). Using the (a) A one-dimensional system with an instantaneous release at t 0 at position x s .

(b) The new frame of reference indicated in blue is moving with a velocity V x . chain rule of differentiation, we have

∂C(θ, t) ∂t = ∂C ∂θ ∂θ ∂t + ∂C ∂t ∂t ∂t = -V ∂C ∂θ + ∂C ∂t , ∂C(θ, t) ∂x = ∂C ∂θ ∂θ ∂x + ∂C ∂t ∂t ∂x = ∂C ∂θ , ∂ 2 C(θ, t) ∂x 2 = ∂ 2 C ∂θ 2 .
(5.11)

Replacing expressions (5.11) in the differential equation (5.10), we have

∂C ∂t = K ∂ 2 C ∂θ 2 , ( 5.12)
with the initial conditions C(θ, t 0 ) = m Ayz δ(θ -θ s ), θ s = 0, and the same boundary conditions in θ. The solution is merely given by equation (5.6):

C(θ, t) = m/A yz (4πK(t -t 0 )) (1/2) exp - θ 2 4K(t -t 0 ) , ( 5.13)
and thus, in function of x,

C(x, t) = m/A yz (4πK(t -t 0 )) (1/2) exp - (x -x s -V x (t -t 0 )) 2 4K(t -t 0 ) . ( 5.14)
In the case of a three-dimensional system, the solution to equation (5.9) can be deduced from the previous expression as shown previously (using variable separation technique) and is given by [START_REF] Kathirgamanathan | Source term estimation of pollution from an instantaneous point source[END_REF]:

C(x, y, z, t) = m (4πK(t -t 0 )) (3/2) exp - (x -x s -V x t) 2 + (y -y s -V y t) 2 + (z -z s -V z t) 2
4K(t -t 0 ) .

(5.15) Figure 5.3 shows the variation of the concentration profile in a one-dimensional system caused by an instantaneous point source located at x s = 0 and emitting a total mass m = 1mg at t 0 = 0, the wind flow is directed along the x-axis and its mean velocity is set to V x = 15m/s.

Diffusion in a semi-infinite medium: method of reflection and superposition

Consider the one-dimensional system in x described here above. In the case of an unbounded domain we have yet considered, half of the emitted material moves in the direction of x > x s and the other half along x < x s .

Assuming an instantaneous point source located at x s of a semi-infinite cylinder extending over x > x s , e.g. imagine having an impermeable boundary at x = x s , all of the emitted substance diffuses along the direction x > x s . We can consider [START_REF] Crank | The mathematics of diffusion[END_REF] that the solution for x < x s is reflected on the boundary x = x s and superposed on the original distribution for x > x s . Note that the solution in the case of an infinite medium given in equation (5.6) is symmetrical about x s , thus the concentration in a semi-infinite medium is given by:

C(x, t) = 2 × m/A yz (4πK(t -t 0 )) (1/2) exp - (x -x s ) 2 4K(t -t 0 ) . ( 5.16)
The principle of reflection and superposition also allows us to deduce the solution to the advection-diffusion equation given in (5.9) for the case of a semi-infinite medium as illustrated through the following example.

Example 5.3.1. Consider a Cartesian coordinate system (x, y, z) where the x-axis is chosen to be oriented in the direction of the mean wind velocity, the y-axis represents the horizontal cross-wind direction and the z-axis is oriented in the upward vertical direction. The mean wind velocity will be denoted as v = (V, 0, 0) [m/s] where V is constant. Suppose that an instantaneous point source is located at (x s , y s , z s ) on an impermeable surface z = z s = 0. The solution in the case of an infinite medium is given in expression (5.15) and can be written for this example as

C(x, y, z, t) = m (4πK(t -t 0 )) (3/2) exp - (x -x s -V t) 2 + (y -y s ) 2 + z 2 4K(t -t 0 ) . ( 5.17)
Using the principle of reflection and superposition, the solution for z < 0 is reflected on the impermeable boundary z = z s = 0 and superposed on the concentration for z > 0, since the original solution is symmetric about z s = 0, the resultant concentration in a semi-infinite medium is given by

C(x, y, z, t) = 2 × m (4πK(t -t 0 )) (3/2) exp - (x -x s -V t) 2 + (y -y s ) 2 + z 2 4K(t -t 0 ) . (5.18)
Figure 5.4 shows the concentration profile at level z = 0 in the case of an instantaneous point source placed at (0, 0) on the impermeable surface z = 0. The wind velocity is chosen to be V = 20m/s and the other parameters K = 25m 2 /s, t 0 = 0 and m = 1mg.

Continuous release source

Let us first consider the case of a single source placed at position r s = (x s , y s , z s ) T on the impermeable plane boundary z = 0. The source emits a chemical with a constant rate Q (g/s) starting at time t 0 . The mean wind velocity is modelled as a constant vector directed along the x-axis and of a constant magnitude V (m/s). The differential equation [START_REF] Badine | Considerations for implementing a wireless network for landmine detection[END_REF][START_REF] Matthes | Optimal weighting of networked electronic noses for the source localization[END_REF] governing the variation of the concentration C(r, t) of the emitted compound at time t ≥ t 0 and at position r = (x, y, z) T in the semi-infinite medium z ≥ 0 is given by:

∂C ∂t -K ∂ 2 C ∂x 2 + ∂ 2 C ∂y 2 + ∂ 2 C ∂z 2 + V ∂C ∂x = 2Q
.u(t -t 0 ).δ(x -x s ).δ(y -y s ).δ(z -z s ) (5.19) where K (m 2 /s) is the isotropic air diffusion coefficient; δ(.) denotes the Dirac delta function and u(t -t 0 ) refers to the step function vanishing for t < t 0 and equal to unity for t ≥ t 0 . With an initial condition C(r, t) = 0, for t < t 0 , the solution [START_REF] Carslaw | Conduction of heat in solids[END_REF] for t ≥ t 0 at position r is given by:

C(Q, r s ,r, t, t 0 ) = Q π 3 2 Kd exp V (x -x s ) 2K × ∞ d 2 √ K(t-t 0 ) exp -u 2 - V 2 d 2 16K 2 u 2 du (5.20)
where d = d(r s , r) = ||r s -r|| is the Euclidean distance between positions r = (x, y, z) T and r s = (x s , y s , z s ) T . Under the assumptions of constant V and Q, and after a sufficiently long time (t → ∞), a stationary concentration profile [START_REF] Badine | Considerations for implementing a wireless network for landmine detection[END_REF][START_REF] Matthes | Optimal weighting of networked electronic noses for the source localization[END_REF], given by the following concentration, is established:

C ∞ (Q, r s , r) = Q 2πKd exp   - V. d -(x -x s ) 2K  
(5.21)

Case of multiple point sources

Let us consider now the case of N point sources where the i -th source, denoted as S i , i = 1, . . . , N , is at position r s i = (x s i , y s i , z s i ) T . Let us further assume that all sources are on the same impermeable plane boundary z = 0 and that the wind velocity is directed along the x-axis.

Using the principle of superposition, and noting that the general advection-diffusion equation and its boundary conditions are linear, the resultant concentration [START_REF] Lushi | An inverse Gaussian plume approach for estimating atmospheric pollutant emissions from multiple point sources[END_REF] for N instantaneous point sources each emitting a mass m i , i = 1, . . . , N at time t 0 evaluated at position r = (x, y, z) T is deduced from equation (5.18) by adding together (superposing) individual solutions and is hence given by:

C(x, y, z, t) = N i=1 m i 4(πK(t -t 0 )) (3/2) exp - (x -x s i -V (t -t 0 )) 2 + (y -y s i ) 2 + z 2
4K(t -t 0 ) .

(5.22) Similarly, for the case of N continuous point sources each emitting the compound at a constant rate Q i , i = 1, . . . , N , the resultant stationary concentration evaluated at position r = (x, y, z) T can be written as:

C r {Q i , r s i } N i=1 , r = N i=1 C ∞ (Q i , r s i , r) (5.23)
where C ∞ (Q i , r s i , r) is given in equation (5.21).

The forward problem for land mine localisation

Recall that the purpose of part II of the thesis is to develop algorithms to locate vapouremitting land mines exclusively using concentration measurements provided by a network of vapour detectors/sensors. The forward model is used to compute an estimated concentration of the explosive chemical at a certain location given a vector of parameters consisting of land mine locations, emission rates and environmental conditions, such as the diffusivity of the air and the wind velocity. It describes the transport of the explosive chemical emitted by the land mines due to the advection and diffusion processes. Note that numerical solutions for modelling the transport of TNT emanating from land mines [START_REF] Irrazábal | Modeling of TNT transport from landmines: Numerical approach[END_REF] were proposed. In our approach, we model a land mine as a point source placed on an impermeable planar surface z = 0 and diffusing an explosive chemical, such as the TNT, and an analytical model for the transport of the explosive vapours is used. We keep the notations and assumptions of section 5.3.3. Two possible scenarios are associated with the land mine localisation problem [START_REF] Jeremić | Landmine detection and localization using chemical sensor array processing[END_REF]:

a. the first scenario considers that after a long period of time, buried land mines no longer emit important amounts of explosive vapours and consequently the corresponding concentrations are not detectable. To overcome this problem, it was proposed to induce the evaporation of the explosive chemical (i.e. using microwaves). In this case, the land mines are considered as instantaneous point sources. This scenario is appropriate to model long-buried land mines.

b. the second scenario assumes that the land mines have been buried since a sufficiently long time for a stationary profile to be established. Here, the land mines are considered as continuous point sources and this scenario is more adapted to recently buried land mines.

Both scenarios will be considered in what follows. Furthermore, in the remaining sections, the sources/mines and the sensors are considered to be in the same plane z = 0, and we omit thus the third coordinate z in all position vectors.

Consider M sensors which are placed at known positions r j = (x j , y j ) T , j = 1, . . . , M . For scenario (a), let us set t 0 = 0, the instant when the evaporation of the explosive chemical is induced and consider that the measurements are taken at time t 1 . Referring to equation (5.22), the estimated concentration at position r j can be written as:

C j {m i , r s i } N i=1 , r j = N i=1 m i × h ji (r s i , r j ) (5.24)
where

h ji (r s i , r j ) = 1 4(πKt 1 ) (3/2) exp - 1 4Kt 1 (x j -x s i -V t 1 ) 2 + (y j -y s i ) 2 .
(5.25)

For scenario (b), we refer to equations (5.21) and (5.23), the estimated concentration at position r j of sensor j can be expressed as:

C j {Q i , r s i } N i=1 , r j = N i=1 Q i × g ji (r s i , r j ) (5.26)
where

g ji (r s i , r j ) = 1 2πKd i j exp   - V. d i j -(x j -x s i ) 2K   (5.27)
and d i j = ||r s ir j ||.

Formulation of the forward problem

Let P denote the vector of parameters of the sources, i.e., sources' positions and emission rates. For N land mines,

P = [{x s i , y s i , m i } N i=1 ] T , (scenario (a)) P = [{x s i , y s i , Q i } N i=1 ] T , (scenario (b))
The concentration measurements provided by the sensor network are grouped in an array denoted by where e j denotes the quantity characterising the difference between measured and mod- elled/estimated concentrations. Equivalently, in vector notation:

Y m = [C m 1 , . . . , C m M ] T . Let Y t = [C t 1 , . . . , C t M ] T denote
Y m = Y t + E m Y e = Y t + E e (5 .31) and:
Y m = Y e + E (5.32) where E refers to an error. Based on (5.31) and (5.24) for scenario (a), (5.26) for scenario (b), we can define an operator F, such that: Y m = F(P, A) (5.33) where A is a constant vector grouping information about the environment, sensors characteristics and model applicability, e.g., wind velocity, air diffusion coefficient and noise variances.

Conclusion

This chapter formulates the forward problem for land mines localisation. This model describes the transport of the explosive chemical emitted by the source into the surrounding area due to advection and diffusion. It allows to predict the concentration at the sensors' positions if the sources parameters (coordinates and emissions) are known. A literature review on the solution of differential advection-diffusion equations was provided and two scenarios, applicable to long-buried and recent-buried land mines, were considered.

THE INVERSE PROBLEM FOR SOURCE CHARACTERISATION

Introduction

In this chapter, we consider the problem of localising an unknown number of land mines exclusively using concentration measurements provided by a wireless sensor network. We adopt a two-step approach consisting in sequentially determining the number of sources, and then characterising them (i.e. computing their coordinates and emissions) by solving the inverse problem associated with the forward model formulated in the previous chapter. We start in section 6.2 with a general overview of inverse problems and their common solution techniques. The first step consisting in determining the unknown number of sources is dealt with in section 6.3 where we use the principal component analysis (PCA), a well known dimension reduction technique, on a matrix of concentration measurements.

Once the number of sources is determined, the inverse problem describing the source characterisation problem is formulated in section 6.4. Two solution techniques are reported, a probabilistic Bayesian approach based on a Markov Chain Monte Carlo (MCMC) sampling scheme, and an optimisation least squares technique. The effectiveness of these schemes is tested on simulated data for two scenarios, considering respectively instantaneous and continuous release sources, in section 6.5. Section 6.6 at last recapitulates the conclusions drawn about this chapter.

6.2 Generalities on inverse problems

General formulation

The inverse problem can be defined as the process of inferring causes, conditioned on knowledge of the effects, as opposed to the forward or direct problem, allowing to determine the effects knowing the causes. Typically, the forward and corresponding inverse problems are linked via a possibly non-linear operator H which represents the relation between a system's input and output. An example of an inverse problem, which will be of main concern in this chapter, is the one of parameter identification [START_REF] Tarantola | Inverse Problem Theory and Methods for Model Parameter Estimation[END_REF]. The corresponding forward problem consists of determining the output of a system knowing the system's parameters. If P denotes a vector of parameters which is linked by an operator H to vector Z grouping the system's outputs (possibly noise-corrupted), then the forward model can be written as Z = H(P), and the inverse problem refers to solving P = H -1 (Z).

A practical difficulty in the study of an inverse problem is that it is often ill-posed [START_REF] Keats | Bayesian inference for source determination in the atmospheric environment[END_REF], meaning that an inverse transformation of the direct model may not exist, may not be unique, and might be unstable.

Common solution techniques

The solution of an inverse problem can be approached in several distinct ways. Amongst the most general and popular techniques are least squares and regularisation approaches.

The key idea for both schemes is to reformulate the inverse problem as an optimisation problem, usually consisting of minimising a functional error between actual measurements and predicted ones obtained by resolving the direct problem. However, for regularisation, the objective function is better-conditioned through the use of regularisation parameters which allows a stable approximation of the inverse operator.

Probabilistic approaches can also be used to address the problem of parameter identification. A primary advantage of probabilistic methods is that the solution takes on the form of a probability distribution rather than a point solution, optimal in terms of a given criterion [START_REF] Keats | Bayesian inference for source determination in the atmospheric environment[END_REF]. In this chapter, we compare the optimisation least squares and the probabilistic Bayesian approaches while citing the advantages and limitations of each scheme.

Determining the unknown number of sources

The purpose of this chapter is to determine unknown land mines/sources' parameters including their positions and emissions. We also consider the case where the number of buried land mines is unknown a priori. In this section, we use the principal component analysis (PCA) in order to determine the number of sources based on a set of measurements. PCA [START_REF] Abdi | Principal component analysis[END_REF] is a popular statistical method that has been widely applied in the analysis of multidimensional data sets, which are usually represented by tables of observations of many possibly inter-correlated variables. Since the information provided by these variables is often redundant, PCA attempts to replace the original set of variables by a smaller number of new variables, called principal components, without losing too much information. This technique considers that the new variables are linear combinations of the original ones and that they are linearly uncorrelated.

Mathematically, PCA transforms the data to a new coordinate system, such that the original set of observations is expressed in terms of the principal components. A technique for performing PCA is to compute the covariance matrix calculated from the available measurements and then to determine its eigenvalues. The corresponding normalised eigenvectors, ordered according to decreasing eigenvalues, define the new coordinate system. A smaller dimensional coordinate system, which is supposed to conserve most of the information, can be obtained by only retaining the eigenvectors associated with the largest eigenvalues.

In our application, M sensors are randomly deployed in a region of interest where the unknown number of vapour-emitting sources is denoted by N . Each sensor provides measurements of the concentration of the explosive material emitted by the sources and transported to the sensor's position due to advection and diffusion processes. Supposing that concentration measurements are taken on each sensor at T time instants, the measurements are grouped in an M × T matrix C = (C jt ), j = 1, . . . , M, t = 1, . . . , T , where the j -th row groups the concentration measurements recorded by the j -th sensor at different time instants and the t -th column groups measurements recorded by all M sensors at time instant t. We explain in the following how the matrix C is obtained for both scenarios (a) and (b).

Scenario (a)

In this case, the emissions of the sources are considered to be induced and thus the land mines are considered as instantaneous sources. We assume that the evaporation of the chemical is induced T times and that, each time, the concentration measurements are taken on the sensors after a constant period of time t 1 . In the following, let m it denote the emission [kg] of source i, i = 1, . . . , N, induced at time t, t = 1, . . . , T ; -M = (m it ) an N × T matrix grouping the emissions of the sources, where each row i is associated with source i and each column t groups the emissions of the sources at time t.

-H = (h ji ) an M × N matrix grouping the factors h ji (r s i , r j ) associated with instantaneous source i and sensor j, as given in equation (5.25) If C t denotes the t-th column of matrix C, i.e., the vector of concentration measurements provided by the sensors at time t and M t denotes the t -th column of matrix M, i.e., the vector of sources emission rates at time t, referring to Equation (5.24), we have in matrix notation:

C t = H × M t , t = 1, . . . , T. (6.1)
Taking into account an additive noise, as in Equation (5.30), we obtain:

C = H × M + , ( 6.2)
where = (e jt ), 1 ≤ j ≤ M, 1 ≤ t ≤ T is an additive noise matrix.

Scenario (b)

For the case of continuous release point sources, since the land mines emissions are likely to change slowly over time, we model these emissions using piecewise constant functions. Referring to figure 5.6, the stationary concentration profile is established within only a few minutes; thus, the recorded concentration measurements can be considered as stationary concentrations established for the constant emission rates between time instants t -1 and t similarly to what was suggested in [START_REF] Borysiewicz | Stochastic algorithm for estimation of the model's unknown parameters via Bayesian inference[END_REF]. Figure 6.1 recapitulates the assumptions we make for scenario (b). It illustrates the variation of the emissions of three sources and indicates, using red boxes, the regions where stationary concentration profiles are established. The measurements grouped in matrix C are, for instance, taken at time instants t 1 , t 2 , etc. Note that these measurements can be chosen by simply examining the signals provided by the sensors and detecting stationary points, i.e. the steady-state concentration values. -Q = (Q it ) is an N × T matrix grouping the emission rates of the sources, where each row i is associated with source i and each column t groups the emissions of the sources at time t, t = 1, . . . , T ;

-G = (g ji ) refers to an an M × N matrix grouping the factors g ji (r s i , r j ) associated with continuous source i and sensor j, as in equation (5.27).

Similarly, if C t (respectively Q t ) denotes the vector of measurements taken at the sensors (the sources' emission rates) at time t, looking back at equation (5.26), we can write

C t = G × Q t , t = 1, . . . , T (6.3)
and again

C = G × Q + , ( 6.4)
where is an additive noise matrix. Equations (6.2) and (6.4) shows that the concentrations measured on the sensors are linear combinations of the sources emissions/emission rates. In our application, the only available information is the matrix of concentrations C. Performing the PCA technique on matrix C and retaining only the eigenvalues, which are larger than some threshold λ th , allows to transform the data to a new coordinate system whose dimension is d ≤ M . An ideal choice of the threshold would recover precisely the number of sources (d = N ). The condition d ≤ M implies that the number of sensors cannot be less than the number of mines.

In our simulations, the number of concentration records is fixed to T = 10, and at each time instant, the emissions (emission rates) of the sources are drawn uniformly from the interval [5 10]µg (µg/s); and a white Gaussian noise of standard deviation σ is added to the estimated concentrations. The sources and the sensors are randomly placed in an L × L planar region.

In order to develop the relationship when N ≤ M , we studied what is the minimal number of sensors that can detect the number N of the true hidden sources using the PCA technique for the simulation conditions described above. In our study, we fixed first the number N of sources and initiated M = N . A hundred source-sensor configurations, obtained by randomly deploying N sources and M sensors in an L × L region, are generated. The PCA technique is tested on the different configurations, and the number of sensors is iteratively increased until the true number of sources is obtained for at least η% of the configurations. For the simulation conditions recapitulated in tables 6.1 and 6.3, an adequate threshold λ th = 10 -7 is empirically determined. The results for scenario (b) are given below in table 6.4. Note that for T = 10 and for the simulation conditions we fixed, the maximal detectable number is N = 8 sources. This technique is tested on the examples given in figures 5.8 and 5.9. First, a matrix of concentration measurements is obtained as described above and a principal component analysis is conducted.

T L [m] K [m 2 /s] V [cm/s] t
The largest eigenvalues for scenario (a) illustrated in figure 5.8 are λ 1 = 0.1428 × 10 -4 , λ 2 = 0.0348 × 10 -4 , λ 3 = 0.0299 × 10 -4 , λ 4 = 0.0135 × 10 -4 ; their number is equal to four, which is the number of the sources. The sum of the remaining eigenvalues is 2.1818×10 -8 .

For scenario (b), we obtain λ 1 = 0.1116×10 -3 , λ 2 = 0.0751×10 -3 , λ 3 = 0.0345×10 -3 , λ 4 = 0.0159 × 10 -3 ; the sum of the remaining eigenvalues is 4.1672 × 10 -8 . Thus, we were able at this first step to successfully determine the number of sources in the considered region. The next step is to localise the sources given concentration measurements.

Formulation of the inverse problem for source characterisation

Referring to Equation (5.33) where A is a constant vector, we can define an operator H, such that: Y m = H(P) (6.5) and introduce the inverse source characterisation problem as finding:

P H -1 (Y m ) (6.6)
The exact solution for Equation (6.6) is usually not tractable. This is primarily due to the existence of model and measurement noises. In other words, the forward model is not perfectly known. Furthermore, due to the commutative property of the addition in Equations (5.24) and (5.26), the solution is not unique; for instance, in the case of two land mines, i.e., N = 2, if

P 1 = [x s 1 , ỹs 1 , Q1 , xs 2 , ỹs 2 , Q2 ] T is a solution to the inverse problem, then P 2 = [x s 2 , ỹs 2 , Q2 , xs 1 , ỹs 1 , Q1 ]
T is also a solution. In general, any parameter vector obtained by only permuting the labels of the sources in a solution vector is also a solution.

The estimation problem is thus over a set of land mines.

Bayesian Inference for Solving the Inverse Problem

In this section, the inverse problem is solved within a probabilistic Bayesian framework. Based on Equation (5.33), the problem consists in determining P having the vector of concentration measurements Y m and some prior knowledge gathered in the constant vector A. In a Bayesian framework, this refers to finding the posterior distribution p(P|Y m , A).

According to Bayes theorem:

p(P|Y m , A) = p(Y m |P, A)p(P|A) p(Y m |A) (6.7)
where p(Y m |P, A) is the measurement likelihood, p(P|A) is the prior distribution and p(Y m |A) is the evidence. The evidence measures the suitability of the model (depending on the number of sources) to the available data [START_REF] Keats | Bayesian inference for source determination in the atmospheric environment[END_REF]. The evidence values are calculated and compared for different models in order to determine the most probable number of sources. The higher the evidence, the better the model can predict the data. Since in our method, the number of sources is determined a priori as described in section 6.3, the evidence is considered as a normalisation factor. Let us consider a bounded domain denoted as Ω for the land mine locations, i.e., the sources lie within a bounded region [x min x max ]×[y min y max ], and the emissions/emission rates are also bounded within lower and upper bounds, ), e e j ∼ N (0, σ 2 e,j ), j = 1, . . . , M then, it can be shown that:

m i ∈ [m min m max ], Q i ∈ [Q min Q max ] i = 1, . . . N .
p(Y m |P, A) ∝ exp   - 1 2 M j=1 (C j -C e j (P)) 2 σ 2 m,j + σ 2 e,j   (6.9)
Sampling directly from this distribution is difficult, and approximate numerical techniques must be used. A widely used approach for estimating the properties of the posterior distribution given in (6.8) is to perform Markov chain Monte Carlo (MCMC) sampling [START_REF] Andrieu | An introduction to MCMC for machine learning[END_REF]. In MCMC algorithms, samples are drawn from the target distribution in the form of a Markov chain where each sample depends on the previous one in the chain. The earliest MCMC algorithm is the random walk Metropolis (RWM) [START_REF] Chib | Understanding the Metropolis-Hastings algorithm[END_REF]. Its basic principle is to sample a candidate value from a proposal distribution depending on the current position of the chain. The candidate is then accepted or rejected according to the Metropolis acceptance probability as will be defined using the following example. Consider sampling from a PDF π(.). If x i-1 denotes the current state of the Markov chain, a trial state z is sampled according to z = x i-1 + u, where u ∼ N (0, Σ), for instance, and Σ denotes a covariance matrix. The candidate z is accepted or rejected according to the Metropolis acceptance probability a given by:

a(x i-1 , z) = min π(z) π(x i-1 ) , 1 if π(x i-1 ) > 0 1 if π(x i-1 ) = 0.
If the candidate is accepted, the chain moves to x i = z; otherwise, the chain remains at x i = x i-1 . The procedure only requires the choice of a proposal function f (.).

While the early RWM algorithm requires the proposal distribution f to be symmetric, i.e., f (x i-1 , z) = f (z, x i-1 ), the Metropolis-Hastings (MH) algorithm [START_REF] Chib | Understanding the Metropolis-Hastings algorithm[END_REF] generalises the approach to non-symmetric proposals. Obviously, the choice of the proposal distribution is crucial to the algorithm convergence. Thus, several procedures were proposed in order to improve the algorithm's convergence. These include, for instance, the adaptive Metropolis algorithm [START_REF] Haario | An adaptive Metropolis algorithm[END_REF], the differential evolution Markov chain Monte Carlo (DE-MC) [START_REF] Ter | A Markov Chain Monte Carlo version of the genetic algorithm differential evolution: easy Bayesian computing for real parameter spaces[END_REF] and the differential evolution adaptive Metropolis (DREAM) algorithm [START_REF] Vrugt | Accelerating Markov Chain Monte Carlo simulation by differential evolution with self-adaptive randomized subspace sampling[END_REF]. Another class of MCMC sampling techniques is the slice sampling technique [START_REF] Neal | Slice sampling[END_REF] and will be used in this paper in order to draw samples from the posterior distribution given in Equation (6.8).

The slice sampling algorithm relies on the observation that sampling from a probability distribution, e.g., π(.) in the case of a univariate distribution, can be done by drawing samples uniformly from the region under the plot of π(.) [START_REF] Neal | Slice sampling[END_REF]. It has an advantage over other MCMC methods, such as the Gibbs sampler and the RWM, in that the magnitude of the changes made to move from one element to the next in the chain is chosen adaptively. . Several methods can be used at this step. The approach adopted here (and illustrated in Figure 6.2b) is called "stepping-out". It requires fixing, a priori, an interval width W and operates as follows: first, set an interval of width W randomly around x 0 . Then iteratively expand this interval in steps of size W and stop when both interval ends become outside the slice {x, f (x) > y}. (c) Draw a value x 1 from the part of the slice that is within the interval determined in (b). The technique used here is referred to as "shrinkage" (see Figure 6.2c) because it picks points uniformly from the determined interval, shrinks this last using points that are outside the slice, and stops whenever finding a point inside it.

Slice sampling can also be used to sample from multivariate distributions. This can be done by updating each variable in turn. It is useful though to note that slice sampling methods, which update all variables of a multivariate distribution simultaneously, do exist [START_REF] Neal | Slice sampling[END_REF].

The Least Squares Technique for Source Characterisation

In this section, we formulate the source characterisation problem as an optimisation problem and propose to solve it using the least squares (LS) approach. LS is a popular method for solving the inverse problems [START_REF] Tarantola | Inverse Problem Theory and Methods for Model Parameter Estimation[END_REF]; it seeks an optimal point solution usually by minimising a quadratic error or cost function between actual measurements and synthetic ones estimated using the forward model.

Referring to Equation (5.32), we propose to solve the inverse source characterisation problem by minimising over P, the vector of unknown parameters, the functional:

J = (Y m -Y e (P)) T (Y m -Y e (P)) (6.10)
The solution is thus given by: P opt = argmin P J (6.11)

Despite the wide applicability, the ease of use and ease of understanding associated with the least squares technique, this method presents a well-known problem: it is sensitive to the convexity of the cost function [START_REF] Matthes | Optimal weighting of networked electronic noses for the source localization[END_REF] and can converge to local optima, thus diverging from the true solution. 5.9. Again, the values of J are calculated as a function of (x s 3 , y s 3 ) after fixing the remaining unknown parameters to their true values. Note that the functional J is not convex and presents local maxima at the sensors' positions. This will cause the convergence of the least squares search algorithm to a local minimum, as we show in the next section, where we test the least squares technique on a simulated scenario and compare its performance to the probabilistic Bayesian approach introduced earlier.

Case studies 6.5.1 Case of an instantaneous release

In this section, we consider the problem of localising N = 4 land mines by randomly deploying M = 30 sensors in a 50 × 50 m 2 planar region according to the scenario shown in figure 5.8. Model and measurement noises are considered to be white Gaussian with an identical standard deviation equal to 0.001µg/m 3 . We fixed the wind velocity to V = 5cm/s and the air diffusion coefficient to K = 25m 2 /s, as in [START_REF] Jeremić | Landmine detection and localization using chemical sensor array processing[END_REF]. Table 6.5 shows the true values of the unknown parameters to be determined; these are the sources positions and emissions. Tab. 6.5: True parameters for scenario (a).

Probabilistic Bayesian Approach

The probabilistic Bayesian approach as described in Section 6.4.1 is tested first. The slice sampling scheme was used in order to draw N p = 4, 000 particles/samples from the posterior distribution defined in Equation 6.7. Figures 6.5 and 6.6 show, respectively, the variation of the log-likelihood of the samples and the evolution of the Markov chain through the iterations. The dimension of the parameter vector is 12. Note on the graphs of figure 6.6 that there is a transition phase (where the samples likelihood is low) before the chain converges to the posterior distribution of interest. This phase is referred to as the burn-in.

In theory, the effect of the initial values tends to zero if the Markov chain is run for an infinite amount of time. In practice, however, an infinite number of samples cannot be drawn, so it is generally assumed that only after a certain number of iterations, the chain reaches the target distribution. Thus, in order to minimise the effect of initial values on the posterior inference, an initial portion, N burn , of a Markov chain samples is discarded and the remaining samples are used to estimate the properties of the posterior distribution. The number N burn of the iterations that will be discarded is called the burn-in number.

Figure 6.7 shows the normalised histograms of the samples corresponding to the different parameters. The empirical distributions of the parameters are also estimated (using kernel density estimation KDE) and shown in red on the same graphs. Note that the empirical PDFs are centred near the true values of the parameters. Figure 6.8 illustrates the true and the estimated positions, which were determined using MCMC slice sampling algorithm in order to solve the Bayesian inference problem. The positions are estimated by computing the mean value. The simulations were carried out using MATLAB on an Intel Core i7-3520M processor (2.90 GHz, 4-MB Cache, Dualcore). The computational time for this approach is 244.60 s. Tab. 6.6: Estimated parameters using the slice sampling for scenario (a). Fig. 6.8: Estimated positions using Bayesian inference and MCMC sampling in the case of instantaneous point sources.

Convergence Diagnostic

Two common critical issues when using an MCMC sampler in order to estimate the properties of a PDF are, first, how to decide when to stop sampling and use the available samples in order to estimate the characteristics of the posterior distribution of interest, and second, how to determine the number of iterations that correspond to the burn-in and should be discarded [START_REF] Cowles | Markov Chain Monte Carlo convergence diagnostics: A comparative review[END_REF]. While it is difficult to predict the number of iterations N iter after which it is safe to stop sampling and the number N burn of initial samples to be discarded, diagnostic tools can be applied to the output of the MCMC samplers in order to address the convergence problem.

In order to decide if the resulting samples accurately estimate the posterior distribution of interest, we apply in this section a convergence diagnostic to the chain outputted by our sampler. We use the StatLib implementation of the Raftery and Lewis diagnostic (1992) [START_REF] Cowles | Markov Chain Monte Carlo convergence diagnostics: A comparative review[END_REF]. This test requires as inputs a posterior quantile of interest q, an acceptable tolerance r for q and a probability s of being within this tolerance. It outputs, amongst other parameters, the number of iterations N iter and burn-ins N burn necessary to satisfy the specified conditions. The diagnostic was run on the resultant Markov chain for q = 0.5, r = 0.01 and s = 0.95, which means we want to measure the 0.5 quantile with an accuracy of 0.01. The output was a total number of iterations N iter = 3672 to be run, of which the first N burn = 14 samples correspond to the burn-in and should be discarded. Thus, we are 95% sure that the true quantile is within ±0.01 from the corresponding estimated value.

Least Squares Technique

Next, the probabilistic Bayesian approach is compared to the generic least squares optimisation approach.

Looking back at Figure 6.3, the cost function to be minimised is convex and has one global minimum. Recall that this figure illustrates the variation of the functional J given by equation (6.10) as a function of Source 3 coordinates after fixing the remaining parameters to their exact values. Figure 6.9 illustrates the optimal land mine positions, resulting in minimising the functional J given in equation (6.10). The least squares technique provides in this example an accurate estimation of the unknown parameters. Both the least squares technique and the Bayesian probabilistic approach manage to provide accurate estimates of the sources' parameters. The least squares approach however offers a better performance and requires less computational time. This was predicted since the cost function to be minimised was convex with one global minimum. Nevertheless, we will see in the next section that the performance of the least squares degrades when the cost function is no longer convex.

Case of a continuous release

In this section, we consider the problem of localising N = 4 land mines, considered as continuous release point sources, by randomly deploying M = 40 sensors in a 20 × 20 m 2 planar region according to the scenario shown in Figure 5.9. Model and measurement noises are considered to be white Gaussian with an identical standard deviation equal to 0.001µg/m 3 . The wind velocity was set to V = 5cm/s and the air diffusion coefficient to K = 25m 2 /s, as in the previous scenario.

In this section, the aim is to localise the continuous release sources given concentration measurements. Table 6.9 shows the true values of the unknown parameters to be determined; these are the sources positions and emission rates. Tab. 6.9: True parameters.

Probabilistic Bayesian Approach

The probabilistic Bayesian approach explained in Section 6.4.1 is tested now. The slice sampling scheme is used in order to draw N p = 4, 000 particles/samples from the posterior distribution defined in the same section. Figures 6.10 and 6.11 show, respectively, the variation of the log-likelihood of the samples and the evolution of the Markov chain through the iterations. In this example, the dimension of the parameter vector is 12. Note again on the graphs of figures 6.10 and 6.11 the burn-in phase before the chain converges to the posterior distribution of interest, the initial samples likelihood is low. Figure 6.12 shows the normalised histograms of the samples corresponding to the different parameters. The empirical PDFs of the parameters are also estimated and shown on the same graphs, they are centred near the true values of the parameters. Figure 6.13 illustrates the true and the estimated positions, which were determined using the slice sampling scheme. The computational time for this approach is 20.35 s. Table 6.10 shows the estimated parameters for the Bayesian probabilistic technique.

Convergence Diagnostic

In order to decide if the resulting samples accurately estimate the posterior distribution of interest, we apply a convergence diagnostic to the chain outputted by the slice sampler. Again, we use the Raftery and Lewis diagnostic.

The test was run on the resultant Markov chain for q = 0.5, r = 0.01 and s = 0.95, which means we want to measure the 0.5 quantile with an accuracy of 0.01. The output was a total number of iterations N iter = 2, 655 to be run, of which the first N burn = 10 samples correspond to the burn-in and should be discarded. Thus, we are 95% sure that the true quantile is within ±0.01 from the corresponding estimated value.

Least Squares Technique

Next, the probabilistic Bayesian approach is compared to the generic least squares optimisation approach.

Looking back at Figure 6.4, the cost function to be minimised is not convex and has multiple local minima. This figure illustrates the variation of the functional J given by Equation (6.10) as a function of Source 3 coordinates after fixing the remaining parameters to their exact values. A least squares search algorithm might fall into some local minimum and, thus, diverge from the true global minimum located near Source 3's true position. Figure 6.14 shows the solution provided by the least squares technique when the vector of parameters is randomly initialised. The algorithm stopped at a local minimum, so the solution provided diverges from the true parameters.

The choice of the initial point for the search algorithm is crucial. If this start point is situated in the restrained convex region around the global minimum (see Figure 6.4), the least squares approach is likely to converge to the true solution. In order to overcome this problem, we choose the start position parameters to be the positions of the three sensors indicating the greatest concentration measurements. These sensors are likely to be the closest to the sources. Figure 6.15 illustrates the optimal land mine positions, resulting in minimising the functional J given in Equation (6.10). The start position parameters for the search algorithm are chosen to be the positions of the four sensors, indicating the maximal concentration measurements. The least squares technique provides in this case an accurate estimation of the unknown parameters. Both the least squares technique and the Bayesian probabilistic approach offer a similar performance provided an adequate choice of the start point for the optimisation search algorithm. The least squares approach also requires less computational time. However, it is important to note that, as the sensors are randomly deployed, even the ones with the greatest concentration measurements might not always be close enough to the sources, so as to find the global minimum. This optimisation technique is very sensitive to the choice of the initial point. On the other side, the probabilistic Bayesian approach together with an efficient sampling algorithm turns out to be more robust and less sensitive to the choice of the initial sample of the Markov chain.

Conclusion

While previous work on land mine localisation using sensor networks solves the problem of locating a single source, this chapter considers the problem of locating several land mines. It also deals with the more difficult scenario of an unknown number of sources to be characterised. First, the PCA technique is used in order to determine the number of land mines. Second, the inverse problem consisting of locating and estimating the emission rates of the land mines is solved in a probabilistic Bayesian framework. In our simulations, we compare the results obtained using this approach with those provided by the least squares optimisation technique. Both methods localise successfully the sources and provide an accurate estimate of the emission rates of multiple land mines. The optimisation least squares offers a good performance when the functional to be minimised is convex as in the case of the instantaneous release point sources. For the scenario considering continuously emitting sources, the cost function presents several local minima. The main advantage of the probabilistic technique is that, using an efficient sampling scheme, it turned out to be less sensitive to the choice of the initial point of the chain, in contrast with the optimisation technique for which the choice of the start point of the search algorithm is crucial. The probabilistic method also makes it possible to quantify uncertainty on the estimated positions, since a PDF of the unknown parameters is obtained, rather than a single optimal point solution.

In this thesis, we derived fusion algorithms for data collected from a WSN. A variety of problems were solved mainly in a Bayesian inference framework. Our work was structured in two parts.

In the first part, we dealt with the problem of inference on graphical models (describing distributed systems). We provided an overview of different classes of graphical models and introduced a message-passing algorithm which provides a distributed scheme for solving the inference problem in a Bayesian framework. This algorithm known as Belief Propagation (BP), also the sum of product algorithm, is a tool for performing exact or approximate marginalisation of joint PDFs represented by graphical structures.

However, for some applications, the problem is not always described by linear models and the uncertainty on the available information is statistically complex. These conditions make the computations of BP in its closed parametric form intractable. A first non-parametric approach known as the Non-parametric Belief Propagation (NBP) was introduced in [START_REF] Sudderth | Nonparametric belief propagation[END_REF] and utilised sample-based representations for probability quantities (as in particle filtering) in order to approximate BP's operations using Monte Carlo techniques. A novel BP scheme, inspired by the box-PF, was developed in this part of the dissertation. It is based on interval representations of the PDFs. This algorithm, which we called box Belief Propagation (box-BP) showed several advantages over the standard Monte Carlo based NBP when tested on both static and dynamic distributed systems. Mainly, it used a smaller number of box-particles to approximate the PDFs which leaded to simpler and faster computations, and to a more compact representation of the information to be exchanged wich in turn implies a reduction of the memory storage, bandwidth an energy consumption necessary to perform inference.

However, for all BP approaches stated earlier, we adopted a parallel update scheme in which, at each iteration, each node sends information to all of its neighbours simultaneously. As future work we can find an optimal message-passing scheme, or an optimal path for information circulation, that results in a reasonable approximation of the posterior PDF of interest with minimal communication requirements.

Perspectives to part I of the thesis also include the application of the box-BP on the localisation of multiple mobile targets, with communications between the mobile targets themselves and between the targets and stationary landmarks. In this case, the system is represented using a dynamic graphical model whose nodes are associated with the targets and the landmarks. At each time step, links are established between each target and the nodes that are within the vicinity of it and thus information can exchanged between a mobile target and its other counterparts in addition to the landmarks. This scenario extends the case we treated in chapter 4 to the problem of mobility tracking of several mobile targets using range measurements. Similarly to what was proposed in chapter 4, messages are exchanged exclusively between nodes representing the local states at one time instant k (even though the graphical model associated with this system is dynamic).

Another perspective is to consider the dynamics of the robots known, thus at each time step a prediction of the future pose can be made and information will circulate between time steps. For this scenario, the nodes of the graphical model describing the system include, in addition to the landmarks and the pose of the mobile targets at time instant k, the positions of the targets at time instants k -1 and/or k + 1; i.e. the dynamic graphical model consists of nodes describing the system's states at consecutive time steps. Thus, range information is exchanged between nodes representing the local states at one time instant k, and message passing is also possible between nodes representing the local state of one mobile target at consecutive time steps since a model for the evolution of the targets is available.

The second part of the thesis considered the problem of characterising (i.e. determining the unknown parameters including the positions and emissions of) an unknown number of sources emitting a chemical/biological agent in their environment. The approach is based on deploying a number of sensors capable of estimating (up to an error level) the concentration of the emitted compound transported to its position due to advection and diffusion. Given a matrix of concentration measurements, we managed to determine the number of sources to be localised using a popular dimension reduction technique. Then the localisation problem was formulated as Bayesian inference and a MCMC sampler was used in order to estimate the properties of the PDF of interest. In this part, the application of interest is that of localisation of anti-personnel land mines. Since two scenarios can be typically associated with this application [START_REF] Jeremić | Landmine detection and localization using chemical sensor array processing[END_REF], the effectiveness of the proposed two-step technique was demonstrated on two distinct scenarios: the first considering instantaneous release sources, the second continuously emitting point sources. However, in our context, sources and sensors were assumed to be in the same plan.

For future work, a model considering a three-dimensional position for the sources and sensors can be employed since the land mines are buried and a third coordinate indicating the depth of the sources is necessary for a more realistic scenario. In this case, the emitted compound is transported by advection-diffusion via a two-layered medium (ground and air). The forward model that predicts the concentration of the emitted chemical at the sensors positions should thus take into consideration the existence of two different transport media with distinct diffusion coefficients.

The validation of the proposed scheme on real data sets is also important especially that using analytical models to describe the transport of the emitted chemical always necessitates fixing simplifying hypothesis (such as homogeneous and/or isotropic media, horizontal separating surface between the layers, constant wind velocity) which is not necessary compliant with the real scenario. Having a real dataset, a neural network could, for instance, be employed in order to predict concentration estimates at any position.

Furthermore, a recent article [START_REF] Sun | Mobile target scenario recognition via low-cost pyroelectric sensing system: Toward a context-enhanced accurate identification[END_REF] presents a hierarchical model to find the groundtruth source bases using Non-negative Matrix Factorisation (NMF) [START_REF] Sun | Non-informative hierarchical bayesian inference for non-negative matrix factorization[END_REF]. Since in our approach, PCA was only used to determine the number of sources, the use of source separation techniques (amongst which are PCA and NMF) for source localisation forms an interesting subject of research. Looking back at equations (6.2) and (6.4), if we are able to reconstruct matrices H (G) and M (Q), having the concentration measurement at a sensor's position (due to all of the sources) it would be possible to determine the contribution of each source to the resultant measured concentration. Thus, the complex problem of localising multiple sources is reduced to several simpler tasks consisting of localising a single source.

1. 1

 1 Main components of a sensor node. . . . . . . . . . . . . . . . . . . . . . . 1.2 The general problem of filtering. . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Particle filtering. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Example of a Bayesian network: the DAG encodes the structure of the distribution p(X) = p(x 1 , x 2 , x 3 , x 4 , x 5 ). . . . . . . . . . . . . . . . . . . . . 1.5 Example of a Bayesian network: a first order Markov chain. . . . . . . . . 1.6 (a) this configuration does not define a clique because no edge interconnects nodes 2 and 3. (b) an example of a clique. . . . . . . . . . . . . . . . . . . 1.7 Example of a Markov Random Field. . . . . . . . . . . . . . . . . . . . . . 1.8 Example of a pairwise MRF. . . . . . . . . . . . . . . . . . . . . . . . . . . 1.9 The inference problem refers to determining the posterior distribution p(X|Y). Shaded nodes are associated with variables for which an observation is available. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.10 A sample graph. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.11 Message update. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.12 Belief computation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.13 Non-parametric Belief Propagation. . . . . . . . . . . . . . . . . . . . . . . 2.1 The resulting image of a box [x] by a function f, an inclusion function [f ] for f and the minimal inclusion function [f ] * . . . . . . . . . . . . . . . . . . 2.2 An example of a CSP with two variables and a single constraint. . . . . . . 2.3 An example of CP contraction technique. . . . . . . . . . . . . . . . . . . . 2.4 The different steps of the box-PF (adapted from [41]). . . . . . . . . . . . . 2.5 The box-BP algorithm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 The propagation model: the small box [x t ] represents the position of sensor t, the red box an inclusion function for the model (3.6) which is merely a disk equation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 A scenario with 100 sensors/nodes and 9 grid-like placed anchors. . . . . . 3.3 The corresponding graph for R = L/4. The sensors form the nodes and an edge indicates an established communication between two sensor nodes. . . 3.4 Results for the NBP algorithm in the case of a grid-like placement of the anchors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Results for the box-BP algorithm in the case of a grid-like placement of the anchors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Results for the MDS-MAP(C) algorithm in the case of a grid-like placement of the anchors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7 A scenario with 100 sensors and 9 randomly spread anchors. The anchors are marked as circles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8 Results for the box-BP algorithm in the case of 6 randomly spread anchors. 3.9 Results for the NBP algorithm in the case of 6 randomly spread anchors. .

3. 10

 10 Results for the MDS-MAP(C) algorithm in the case of 6 randomly spread anchors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Simulated scenario. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 The directed acyclic graphs representing the joint probability distribution at two consecutive time steps k -1 and k. . . . . . . . . . . . . . . . . . . 4.3 The result of propagating a box through the disk equation given by (4.3). An inclusion function for this propagation model is also given. . . . . . . . 4.4 Trilateration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Cumulative distribution function of the error: case of a Gaussian measurement noise. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 Cumulative distribution function of the error: case of a bounded measurement noise. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.7 Real and estimated trajectories for box-BP, NBP and Sayed05 (case of a bounded measurement noise with a = 2m). . . . . . . . . . . . . . . . . . . 4.8 The cylindrical landmarks and robots used for the generation of the UTIAS dataset. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.9 The workspace for the real data set. . . . . . . . . . . . . . . . . . . . . . . 4.10 The workspace for the real dataset. . . . . . . . . . . . . . . . . . . . . . . 4.11 Positions of the landmarks in the real workspace. . . . . . . . . . . . . . . 4.12 Cumulative distribution function of the error for the real data set. . . . . . 5.1 Concentration profile: case of an instantaneous source and assuming the transport of the emitted material occurs exclusively by diffusion. . . . . . . 5.2 The stationary and the new moving coordinate systems. . . . . . . . . . . 5.3 Variation of the concentration in space and time in the case of a 1D instantaneous source and assuming the transport of the emitted material occurs by advection-diffusion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 Concentration profile in the case of an instantaneous source placed in a semi-infinite medium and assuming the transport of the emitted material occurs by advection-diffusion. . . . . . . . . . . . . . . . . . . . . . . . . . 5.5 An illustrative example with a single continuous point source. . . . . . . . 5.6 Concentration profiles are given in solid line at positions 1, 2 and 3 for the example given in figure 5.5. Dotted lines show the stationary concentrations. 5.7 The stationary concentration profile for the example in figure 5.5. . . . . . 5.8 Illustration of the forward model for scenario (a). . . . . . . . . . . . . . . 5.9 Illustration of the forward model for scenario (b). . . . . . . . . . . . . . . 6.1 The emission rates are piecewise constant, and a series of stationary concentration profiles are established. . . . . . . . . . . . . . . . . . . . . . . . 6.2 Slice sampling (adapted from [78]). . . . . . . . . . . . . . . . . . . . . . . 6.3 Functional to be minimised in the logarithmic scale, case of instantaneous release sources. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.4 Functional to be minimised in the logarithmic scale, case of continuous point sources. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.5 Variation of the log-likelihood in the Markov chain for scenario (a). . . . . 6.6 Evolution of the Markov chain for scenario (a). . . . . . . . . . . . . . . . 6.7 Samples empirical distributions for scenario (a). . . . . . . . . . . . . . . . 6.8 Estimated positions using Bayesian inference and MCMC sampling in the case of instantaneous point sources. . . . . . . . . . . . . . . . . . . . . . .3.3 Multidimensional scaling for self localisation in WSNs . . . . . . . . . . . . 3.4 Self localisation in WSNs as inference on a graphical model . . . . . . . . . 3.4.1 Problem formulation . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.2 Implemented solutions . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.3 Simulation results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4. Localisation of a mobile object in a wireless sensor network . . . . . . . . . . . . 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Literature review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Localisation of a mobile target as inference on a graphical model . . . . . . 4.3.1 Simulated Scenario . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.2 Results Obtained on Real Data . . . . . . . . . . . . . . . . . . . . 4.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Part II Source characterisation problem 5. Literature review on source characterisation and advection-diffusion models . . . 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Related works . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Advection-diffusion models . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3.1 Instantaneous release source . . . . . . . . . . . . . . . . . . . . . . 5.3.2 Continuous release source . . . . . . . . . . . . . . . . . . . . . . . 5.3.3 Case of multiple point sources . . . . . . . . . . . . . . . . . . . . . 5.4 The forward problem for land mine localisation . . . . . . . . . . . . . . . 5.4.1 Formulation of the forward problem . . . . . . . . . . . . . . . . . . 5.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6. The inverse problem for source characterisation . . . . . . . . . . . . . . . . . . 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Generalities on inverse problems . . . . . . . . . . . . . . . . . . . . . . . . 6.2.1 General formulation . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2.2 Common solution techniques . . . . . . . . . . . . . . . . . . . . . . 6.3 Determining the unknown number of sources . . . . . . . . . . . . . . . . . 6.4 Formulation of the inverse problem for source characterisation . . . . . . .

Fig. 1 . 1 :

 11 Fig. 1.1: Main components of a sensor node.

Fig. 1 . 2 :

 12 Fig. 1.2: The general problem of filtering.

Fig. 1 . 3 :

 13 Fig. 1.3: Particle filtering.

Fig. 1 . 4 :

 14 Fig. 1.4: Example of a Bayesian network: the DAG encodes the structure of the distribution p(X) = p(x 1 , x 2 , x 3 , x 4 , x 5 ).

Fig. 1 . 5 :

 15 Fig. 1.5: Example of a Bayesian network: a first order Markov chain.

Figure 1 .

 1 [START_REF] Andrieu | An introduction to MCMC for machine learning[END_REF] illustrates an example of a BN associated with a joint PDF defined over N = 5 random variables.

Fig. 1 . 6 :

 16 Fig. 1.6: (a) this configuration does not define a clique because no edge interconnects nodes 2 and 3. (b) an example of a clique.

Fig. 1 . 7 :Example 1 . 5 . 2 .

 17152 Fig. 1.7: Example of a Markov Random Field.

Fig. 1 . 8 :

 18 Fig. 1.8: Example of a pairwise MRF.

Fig. 1 . 9 :

 19 Fig. 1.9: The inference problem refers to determining the posterior distribution p(X|Y). Shaded nodes are associated with variables for which an observation is available.

Fig. 1 .

 1 Fig. 1.10: A sample graph.

Fig. 1 . 11 :

 111 Fig. 1.11: Message update.

Fig. 1 . 12 :

 112 Fig. 1.12: Belief computation.

Fig. 1 . 13 :

 113 Fig. 1.13: Non-parametric Belief Propagation.

Example 2 . 2 . 1 .

 221 The extension to intervals of the usual arithmetic operations {+, -, ×, /} is given by the following expressions where [x] = [x, x] and [y] = [y, y] are two intervals in R.

Fig. 2 . 1 :

 21 Fig. 2.1: The resulting image of a box [x] by a function f, an inclusion function [f ] for f and the minimal inclusion function [f ] * .

Fig. 2 . 2 :

 22 Fig. 2.2: An example of a CSP with two variables and a single constraint.

  using the CP technique. a. Forward propagation: Consider the direct form of the constraint x 2 = f (x 1 ), (a) Initial domain. (b) Forward propagation. (c) Backward propagation.

Fig. 2 . 3 :

 23 Fig. 2.3: An example of CP contraction technique.

Fig. 2 . 4 :

 24 Fig. 2.4: The different steps of the box-PF (adapted from [41]).

a)

  Box particle initialisation: A prior bounded state space region is split into N , equally weighted and mutually disjoint, boxes {[x ( ) 0 ]} N =1 . This initialisation using boxes allows to explore a large prior uncertainty region using only a few box particles. b) Time update step: Knowing the cloud of box particles {[x ( )

  .14) A uniform PDF U [x] is constant over its support (and equal to 1/|[x]|) and equal to zero everywhere else. Note that, by abuse of notation, |[x]|, where x ∈ IR nx and n x is the state space dimension, denotes the product nx =1 |[x( )]|, e.g. |[x]| denotes the area of the rectangle defined by x in 2D, the volume of the rectangular cuboid in 3D, etc.

Fig. 2 . 5 :

 25 Fig. 2.5: The box-BP algorithm.

2. 5

 5 The box-BP: a generalisation of the box-PF Same as the NBP generalises the PF to general inference problems on arbitrary graphical models, wed can demonstrate that the box-BP is also a generalisation of the box-PF.Consider the special case of a temporal Markov chain. At each time step k an observation y k is available (refer to figure 1.2). A convenient parametrisation of the graph is as follows:

Fig. 3 . 1 :

 31 Fig. 3.1: The propagation model: the small box [x t ] represents the position of sensor t, the red box an inclusion function for the model (3.6) which is merely a disk equation.

Fig. 3 . 2 :

 32 Fig. 3.2: A scenario with 100 sensors/nodes and 9 grid-like placed anchors.

Fig. 3 . 3 :

 33 Fig. 3.3: The corresponding graph for R = L/4. The sensors form the nodes and an edge indicates an established communication between two sensor nodes.

Tab. 3 . 1 :

 31 Simulation results for a grid-like placement of the anchors.

Fig. 3 . 4 :

 34 Fig. 3.4: Results for the NBP algorithm in the case of a grid-like placement of the anchors.

  , 3.5 and 3.6 illustrate the results obtained for the NBP, box-BP and MDS-MAP(C) algorithms respectively. The dots denote the true sensors positions whereas the lines indicate the error on these positions (the distance between the true and the estimated positions).

Fig. 3 . 5 :

 35 Fig. 3.5: Results for the box-BP algorithm in the case of a grid-like placement of the anchors.

Fig. 3 . 6 :

 36 Fig. 3.6: Results for the MDS-MAP(C) algorithm in the case of a grid-like placement of the anchors.

Fig. 3 . 7 : 87 Tab. 3 . 3 :

 378733 Fig. 3.7: A scenario with 100 sensors and 9 randomly spread anchors. The anchors are marked as circles.

30 Tab. 3 . 4 :

 3034 1.70 1.71 1.71 time (sec) 0.30 0.29 0.30 0.Results for the MDS-MAP(C) algorithm in the case of randomly spread anchors.

Fig. 3 . 8 :

 38 Fig. 3.8: Results for the box-BP algorithm in the case of 6 randomly spread anchors.

Fig. 3 . 9 :

 39 Fig. 3.9: Results for the NBP algorithm in the case of 6 randomly spread anchors.

Fig. 3 . 10 :

 310 Fig. 3.10: Results for the MDS-MAP(C) algorithm in the case of 6 randomly spread anchors.

Fig. 4 . 1 :

 41 Fig. 4.1: Simulated scenario.

Figure 4 .

 4 1 shows a simulated scenario with n = 200 sensors deployed in a planar L × L region with L = 100 m. The trajectory of the mobile target is also shown.

Fig. 4 . 2 :

 42 Fig. 4.2: The directed acyclic graphs representing the joint probability distribution at two consecutive time steps k -1 and k.

(i) 1 ,

 1 [x 1 ] (i) } of N weighted box-particles obtained using Algorithm 1. The forwarded message from 1 to 2, separated by noisy distance d 12 , is denoted as m (i) 12 and is represented by the cloud of boxes {ω (i) 12 , x (i) 12 } obtained, as shown in (2.22), by propagating the collection {[x 1 ] (i) } through the following model

Fig. 4 . 3 :

 43 Fig. 4.3: The result of propagating a box through the disk equation given by (4.3). An inclusion function for this propagation model is also given.

Figure 4 .

 4 Figure 4.3 illustrates the result of propagating a box through the model (4.3). In this example, the position of the sending node is [x 1 ], the receiving node 2 is at a distance d 12 = [d d] from the sending node 1.Algorithm 5 describes the message passing procedure to perform the mobility tracking task.For NBP algorithm, the message product at the sending node 1 being represented by a collection {ω i 1 , x i 1 } of weighted particles, the forwarded message from 1 to 2, separated by noisy distance d 12 , is computed by propagating the samples {x i 1 } through the model (4.3) where θ ∼ U ([0, 2π]).Furthermore, both probabilistic algorithms are also compared to the optimisation approach introduced by Sayed et al. in[START_REF] Sayed | Network-based wireless location: challenges faced in developing techniques for accurate wireless location information[END_REF] and to which we will be referring as Sayed05. A brief review of the Sayed05 algorithm is given in the next subsection.

Fig. 4 . 4 :

 44 Fig. 4.4: Trilateration.

  (a) σ = 0.5 m. (b) σ = 1 m. (c) σ = 1.5 m.

Fig. 4 . 5 :

 45 Fig. 4.5: Cumulative distribution function of the error: case of a Gaussian measurement noise.

  and 4.6 illustrating the empirical cumulative distribution function of the error for all three algorithms and for both Gaussian and bounded measurement error environments respectively.

Figure 4 .

 4 7 illustrates the real trajectory and the estimated trajectories obtained for box-BP, NBP and Sayed05 algorithms, respectively. The presented results correspond to a bounded error environment with a = 2 m.

  (a) a = 0.5 m. (b) a = 2 m. (c) a = 4 m.

Fig. 4 . 6 :

 46 Fig. 4.6: Cumulative distribution function of the error: case of a bounded measurement noise.

Fig. 4 . 7 :

 47 Fig. 4.7: Real and estimated trajectories for box-BP, NBP and Sayed05 (case of a bounded measurement noise with a = 2m).

Fig. 4 . 8 :

 48 Fig. 4.8: The cylindrical landmarks and robots used for the generation of the UTIAS dataset.

Fig. 4 . 9 :

 49 Fig. 4.9: The workspace for the real data set.

Fig. 4 . 10 :Tab. 4 . 3 :

 41043 Fig. 4.10: The workspace for the real dataset.

Fig. 4 . 11 :

 411 Fig. 4.11: Positions of the landmarks in the real workspace.

Fig. 4 . 12 :

 412 Fig. 4.12: Cumulative distribution function of the error for the real data set.

Fig. 5 . 1 :

 51 Fig. 5.1: Concentration profile: case of an instantaneous source and assuming the transport of the emitted material occurs exclusively by diffusion.

Fig. 5 . 2 :

 52 Fig. 5.2: The stationary and the new moving coordinate systems.

Fig. 5 . 3 :

 53 Fig. 5.3: Variation of the concentration in space and time in the case of a 1D instantaneous source and assuming the transport of the emitted material occurs by advection-diffusion.

Fig. 5 . 4 :

 54 Fig. 5.4: Concentration profile in the case of an instantaneous source placed in a semi-infinite medium and assuming the transport of the emitted material occurs by advectiondiffusion.

Fig. 5 . 5 :

 55 Fig. 5.5: An illustrative example with a single continuous point source.

Fig. 5 . 6 :

 56 Fig. 5.6: Concentration profiles are given in solid line at positions 1, 2 and 3 for the example given in figure 5.5. Dotted lines show the stationary concentrations.

Figures 5 .

 5 Figures 5.5 and 5.6 show respectively an illustrative example with one continuous point source placed at position (x s = 5, y s = 5) in the plane z = 0, and the variation of the concentration at three different positions in a 20 m × 20 m planar region, i.e. at level z = 0, as given in equation (5.20). The parameters are set to Q = 5mg/s, K = 25m 2 /s and V = 5cm/s. The stationary concentrations determined using expression (5.21) are also shown in dotted lines. The graph shows that the concentration change within 5 minutes reaches 97.50, 95.96 and 91.40 percent of the stationary concentration that would be established at the positions 1, 2 and 3, respectively. Figure 5.7 illustrates the stationary concentration profile for this example.

Fig. 5 . 7 :

 57 Fig. 5.7: The stationary concentration profile for the example in figure 5.5.

  the vector of true concentrations at positions r j , andE m = [e m 1 , . . . , e m M ]T a vector of measurements error, thus:C m j = C t j + e m j , j = 1, . . . , M (5.28) Furthermore, if Y e = [C e 1 , . . . , C e M ]T refers to the vector of estimated concentrations at r j , i.e., the concentrations obtained by resolving the forward problem according to equations(5.24) or(5.26), and E e = [e e 1 , . . . , e e M ] T a vector of model error, then:C e j = C t j + e e j , j = 1, . . . , M (5.29)Based on equations (5.28) and (5.29), we can write:C m j = C e j + e j , j = 1, . . . , M(5.30)

( a )

 a Concentration profile at t 1 = 2s for a scenario with N = 4 instantaneous point sources and M = 30 sensors. (b) The corresponding perspective view.

Fig. 5 . 8 :

 58 Fig. 5.8: Illustration of the forward model for scenario (a).

Figure 5 .

 5 8 shows the concentration profile at t 1 = 2s of the emitted chemical determined using equation (5.22) for a scenario with N = 4 instantaneous point sources and M = 30 sensors, which are randomly deployed in a 50 × 50 m 2 planar region. Emissions are randomly fixed to m 1 = 9.78 µg, m 2 = 8.43 µg, m 3 = 7.58 µg, and m 4 = 9.94 µg.

Figure 5 .

 5 9 on the other hand shows the stationary concentration profile of the diffused explosive chemical determined using equation (5.23) for a scenario with N = 4 continuous point sources and M = 40 sensors, which are randomly deployed in a 20 × 20 m 2 planar region. Emission rates are randomly fixed to Q 1 = 8.01 µg/s, Q 2 = 5.38 µg/s, Q 3 = 5.87 µg/s and Q 4 = 9.39 µg/s. For both examples, the environmental parameters are set to V = 5cm/s and K = 25m 2 /s.

( a )

 a Stationary concentration profile for a scenario with N = 4 continuous point sources and M = 40 sensors. (b) The corresponding perspective view.

Fig. 5 . 9 :

 59 Fig. 5.9: Illustration of the forward model for scenario (b).

Fig. 6 . 1 :

 61 Fig. 6.1: The emission rates are piecewise constant, and a series of stationary concentration profiles are established.

Tab. 6 . 4 :

 64 The minimal number of sensors required in order to detect the true number of continuous release sources for the simulation conditions fixed in table 6.3.

Fig. 6 . 2 :

 62 Fig. 6.2: Slice sampling (adapted from [78]).

Figure 6 .

 6 Figure 6.2 illustrates the operations of the slice sampling algorithm in the case of a univariate target distribution π(.). The procedure requires only the knowledge of a function f (.) that is proportional to π(.). It operates iteratively in three steps: (a) Starting from the current position of the chain, denoted as x 0 , and such that f (x 0 ) > 0, draw a value y uniformly from the interval [0, f (x 0 )]. The horizontal slice defined by y consists of the values of x for which f (x) > y (see Figure 6.2a). (b) Find an interval around x 0 comprising the majority, or the totality, of the slice defined in (a). Several methods can be used at this step. The approach adopted here (and

Fig. 6 . 3 :

 63 Fig. 6.3: Functional to be minimised in the logarithmic scale, case of instantaneous release sources.

Figure 6 .

 6 Figure 6.3 illustrates the variation of the criterion J in terms of Source 2 coordinates for the scenario shown in figure 5.8. The values of J are calculated as a function of (x s 2 , y s 2 ) after fixing the remaining unknown parameters, i.e., [{x s i , y s i , Q i } i=1,3,4 , Q 2 ] T , to their true values. Note that the functional J is convex and presents one global maximum near Source 2's true position. The least squares search algorithm is expected to provide accurate position estimates in this case.

Fig. 6 . 4 :

 64 Fig. 6.4: Functional to be minimised in the logarithmic scale, case of continuous point sources.

Figure 6 .

 6 Figure 6.4 illustrates the variation of the criterion J in terms of Source 3 coordinates for the scenario shown in Figure5.9. Again, the values of J are calculated as a function of (x s 3 , y s 3 ) after fixing the remaining unknown parameters to their true values. Note that the functional J is not convex and presents local maxima at the sensors' positions. This will cause the convergence of the least squares search algorithm to a local minimum, as we show in the next section, where we test the least squares technique on a simulated scenario and compare its performance to the probabilistic Bayesian approach introduced earlier.

Fig. 6 . 5 :

 65 Fig. 6.5: Variation of the log-likelihood in the Markov chain for scenario (a).

Fig. 6 . 6 :

 66 Fig. 6.6: Evolution of the Markov chain for scenario (a).

Fig. 6 . 7 :

 67 Fig. 6.7: Samples empirical distributions for scenario (a).

Fig. 6 . 9 :

 69 Fig. 6.9: Optimal solution provided by the least squares approach (scenario (a)).

  i

Fig. 6 . 10 :

 610 Fig. 6.10: Variation of the log-likelihood in the Markov chain for scenario (b).

Fig. 6 . 11 :

 611 Fig. 6.11: Evolution of the Markov chain for scenario (b).

Fig. 6 . 12 :

 612 Fig. 6.12: Samples empirical distributions for scenario (b).

Fig. 6 . 13 :

 613 Fig. 6.13: Estimated positions of the continuous release sources using Bayesian inference and MCMC sampling.

Fig. 6 . 14 :

 614 Fig. 6.14: Solution provided by the least squares approach for randomly chosen start parameters.

Fig. 6 . 15 :

 615 Fig. 6.15: Optimal solution provided by the least squares approach in the case of continuous release sources.

  .1) where v ts refers to a noise, p v (x t , x s ) denotes a noise probability distribution and || . || is the Euclidean distance between t and s. Note that both sensors t and s are assumed to have the same estimation d ts of the distance separating them. The initialisation of the sensor network for the self localisation procedure is recapitulated hereafter, mainly, each sensor node t has to: a. obtain local information p t (x t ), if available. b. broadcast its ID and listen for other sensors' broadcasts. c. estimate the distance d ts separating it from any sensor s amongst the set of received broadcasts/IDs. d. communicate with neighbouring nodes to symmetrise the estimated distances. Consider a binary random variable o ts indicating whether the distance d ts is observed or not, i.e. whether a communication between sensors t and s was established or not, thus, according to the assumption above:

  1 N are uniformly generated from the available observations at the anchors. Algorithm 4 BP's message update for each node t in the sensor network Γ

	end if
	end for
	Message product
	if M sgCounter == 0 then
	if t is an anchor node then
	use the available measurement at t to generate uniformly N boxes mutually disjoint
	and weighted equally with 1 N
	else
	no message can be forwarded from t to s
	end if
	else
	for i = 1, . . . , M sgCounter do
	combine the messages sequentially using the series of intersections in Algorithm 1
	contract the resulting boxes by intersecting them with the available observation at
	t (if any)
	end for
	end if
	propagate the boxes through the model (3.6)
	normalise the weights

t ← {u ∈ V, o tu = 1} s ∈ Γ t M sgCounter ← 0 for u ∈ Γ t and u = s do if a message is incoming from u to t then M sgCounter ← M sgCounter + 1

At an iteration i, the message product at node t is represented by a collection {ω

  particles error [m] time [ms] error [m] time [ms] error [m] time [ms] Simulation results in the case of a Gaussian measurement noise. Algorithm N o of particles error [m] time [ms] error [m] time [ms] error [m] time [ms]

	box-BP	9	0.3284	2	0.8287	2.8	1.6878	3.3
	NBP	200	0.8317	59.1	1.0704	57.1	1.5069	58.4
	Sayed05	N/A	0.1923	0.18	0.8567	0.17	1.9835	0.16
	Tab. 4.1: a = 0.5 m	a = 2 m	a = 4 m
	box-BP	9	0.3808	1.8	0.6964	2.3	1.3885	3
	NBP	200	0.9325	59.2	2.2367	57.7	1.5032	57.3
	Sayed05	N/A	0.2339	0.18	0.9488	0.16	2.0321	0.24

Tab. 4.2: Simulation results in the case of a bounded measurement noise.

  The minimal number of sensors required in order to detect the true number of instantaneous sources for the simulation conditions considered in table 6.1.

	1 [s] η% 2 80 Tab. 6.1: Simulation parameters for scenario (a). 10 50 25 5 The results for scenario (a) are given in table 6.2. N 2 3 4 5 6 7 M 6 11 22 45 100 206 10 20 25 5 90 Tab. 6.2: T L [m] K [m 2 /s] V [cm/s] η% Tab. 6.3: Simulation parameters for scenario (b).

  Choosing a non informative distribution for the prior, i.e., a uniform PDF, we can write:

	e m j	∼ N (0, σ 2 m,j

p(P|Y m , A) ∝ 1 P∈Ω p(Y m |P, A) (6.8)

where 1 P∈Ω denotes the indicator function taking on a value 1 if P ∈ Ω and 0 if not. Additionally, if the measurement and model noises are assumed to be white and Gaussian, i.e.:

Table 6 .

 6 6 shows the estimated parameters using the slice sampling technique.

	i	(x s i , y s i )	m i [µg]
	1 (0.53, 3.65)	10.08
	2 (44.90, 30.82)	8.33
	3 (4.95, 26.30)	8.20
	4 (44.49, 13.29) 10.22

Table 6 .

 6 7 groups the optimal results obtained using the least squares technique. The computational time is 1.99 s.Tab. 6.7: Estimated parameters using the Least Squares search algorithm for scenario (a).

	i	(x s i , y s i )	m i [µg]
	1 (0.95, 5.94)	9.59
	2 (44.78, 30.50)	8.29
	3 (5.04, 27.03)	7.45
	4 (44.35, 13.14) 10.04

Table 6 .

 6 8 groups the mean squared errors on sources positions and emission rates for both probabilistic Bayesian and optimisation least squares approaches. Mean squared errors on positions and emission rates in the case of instantaneous release sources.

	Approach	Error on Position Emission Rate
	Slice sampling	1.207	0.140
	Least Squares	0.088	0.020
	Tab. 6.8:

  (x s i , y s i )

		Q i [µg/s]
	1 (6.99, 2.18)	8.01
	2 (1.73, 12.27)	5.38
	3 (16.52, 6.45)	5.87
	4 (14.11, 13.01)	9.39

Table 6 .

 6 11 groups the optimal results obtained using the least squares technique. The computational time is 6.34 s.Tab. 6.11: Estimated parameters using the Least Squares search algorithm, case of continuous sources.

	i	(x s i , y s i )	Q i [µg/s]
	1 (6.99, 2.22)	8.53
	2 (1.78, 12.49)	5.06
	3 (16.54, 6.54)	5.77
	4 (14.26, 13.01)	8.98

Table 6 .

 6 [START_REF] Bishop | Pattern Recognition and Machine Learning, chapter 8: Graphical models[END_REF] groups the mean squared errors on sources positions and emission rates for both approaches.Tab. 6.12: Mean squared errors on positions and emission rates for scenario (b).

	Approach	Error on Position Emission Rate
	Slice sampling	0.040	0.223
	Least Squares	0.021	0.137

Résumé de la thèse en françaisL'objectif de la thèse est de développer des algorithmes de fusion de données recueillies à l'aide d'un réseau de capteurs sans fil afin de localiser plusieurs sources émettant un agent chimique ou biologique dans l'air. Ces capteurs détectent la concentration de la substance émise, transportée par advection et diffusion, au niveau de leurs positions et de communiquer cette information à un centre de traitement. L'information recueillie de façon collaborative est d'abord utilisée pour localiser les capteurs déployés au hasard et ensuite pour localiser les sources. Les applications comprennent, entre autres, la surveillance environnementale et la surveillance de sites sensibles ainsi que des applications de sécurité dans le cas d'une libération accidentelle ou intentionnelle d'un agent toxique. Toutefois, l'application considérée dans la thèse est celle de la détection et la localisation de mines terrestres. Dans cette approche, les mines sont considérées comme des sources émettrices de produits chimiques explosifs. La thèse comprend une contribution théorique où nous étendons l'algorithme de propagation de la croyance, un algorithme de fusion de données bien connu et largement utilisé pour l'estimation collaborative d'état dans les réseaux de capteurs, au cadre des méthodes à erreurs bornées. Le nouvel algorithme est testé sur le problème de l'auto-localisation dans les réseaux de capteurs statiques ainsi que l'application de suivi d'un objet mobile en utilisant un réseau de capteurs de distance. Autres contributions comprennent l'utilisation d'une approche probabiliste Bayésienne avec des techniques d'analyse de données pour localiser un nombre inconnu de sources émettrices de vapeur.

In graph theory, the set of descendants of a node comprises the set of nodes that can be reached from this node on a direct path[START_REF] Ben-Gal | Encyclopedia of Statistics in Quality and Reliability, chapter Bayesian Networks[END_REF].

ACKNOWLEDGEMENTS

This PhD dissertation represents not only my efforts, but also the time, energy, support, and encouragement of my mentors, colleagues, and friends. I am grateful for this opportunity to acknowledge all those who made this thesis possible.

First and foremost I thank my supervisor at the UTC, Mr. Fahed Abdallah, a person who is exceptional in all ways, scientifically, professionally and especially on the human side, for sharing his vision and experience, for putting me on track in various fields, and Part II SOURCE CHARACTERISATION PROBLEM

, Mr. Amadou Gning, for sharing his knowledge and talent, for making complex material accessible and for all the help and goodwill he showed whenever I met a difficulty. I also acknowledge the supports from the Lebanese CNRS and from the U.K. Engineering and Physical Sciences Research Council (EPSRC) through the efforts and collaboration of Mrs. Lyudmila Mihaylova from the University of Sheffield.

Looking back, I owe a tremendous debt of gratitude to the director and my professor at the faculty of engineering, Mr.

LIST OF TABLES