An adaptive Bayesian inference algorithm to estimate the parameters of a hazardous atmospheric release

Atmospheric Environment 122 (2015) 748e762

Contents lists available at ScienceDirect

Atmospheric Environment journal homepage: www.elsevier.com/locate/atmosenv

An adaptive Bayesian inference algorithm to estimate the parameters of a hazardous atmospheric release Harizo Rajaona a, b, *, François Septier a, Patrick Armand b, Yves Delignon a, Christophe Olry c, Armand Albergel c, Jacques Moussafir c Institut Mines-T
el
ecom/T
el
ecom Lille/LAGIS UMR CNRS 8219, 59650 Villeneuve-d'Ascq, France CEA-DAM, DIF, F-91297 Arpajon, France c ARIA Technologies, 8 rue de la Ferme, 92100 Boulogne-Billancourt, France a

b

h i g h l i g h t s
We
We
We
We
We

tackle the problem of estimating the location and the distribution of an atmospheric polluting source. use Bayesian inference as a main framework, to be able to quantify the uncertainty. derive an analytical way to estimate the release-rate by exploiting Gaussian assumptions on its prior distribution. estimate the location by using an adaptive Monte Carlo algorithm based on the principles of Importance Sampling. test our algorithm on a combination of synthetic and experimental data taken from the FFT07 experiment.

a r t i c l e i n f o

a b s t r a c t

Article history: Received 9 January 2015 Received in revised form 8 October 2015 Accepted 9 October 2015 Available online 23 October 2015

In the eventuality of an accidental or intentional atmospheric release, the reconstruction of the source term using measurements from a set of sensors is an important and challenging inverse problem. A rapid and accurate estimation of the source allows faster and more efficient action for first-response teams, in addition to providing better damage assessment. This paper presents a Bayesian probabilistic approach to estimate the location and the temporal emission profile of a pointwise source. The release rate is evaluated analytically by using a Gaussian assumption on its prior distribution, and is enhanced with a positivity constraint to improve the estimation. The source location is obtained by the means of an advanced iterative Monte-Carlo technique called Adaptive Multiple Importance Sampling (AMIS), which uses a recycling process at each iteration to accelerate its convergence. The proposed methodology is tested using synthetic and real concentration data in the framework of the Fusion Field Trials 2007 (FFT-07) experiment. The quality of the obtained results is comparable to those coming from the Markov Chain Monte Carlo (MCMC) algorithm, a popular Bayesian method used for source estimation. Moreover, the adaptive processing of the AMIS provides a better sampling efficiency by reusing all the generated samples. © 2015 Elsevier Ltd. All rights reserved.

Keywords: Source term estimation Bayesian inference Monte-Carlo techniques Adaptive multiple importance sampling

1. Introduction The threat of chemical, radiological, biological, and nuclear (CRBN) releases in the atmosphere is a key issue. Such incidents may be the consequence of intentional releases using nonconventional methods in order to create panic. The origin of

* Corresponding author. CEA-DAM, DIF, F-91297 Arpajon, France. E-mail address: [email protected] (H. Rajaona). http://dx.doi.org/10.1016/j.atmosenv.2015.10.026 1352-2310/© 2015 Elsevier Ltd. All rights reserved.

these releases can also be accidental, for example given a leak of hazardous material on an industrial site. Either way, the development of tools to reconstruct the source is of utmost importance for the population safety as well as for the efficiency of the firstresponse teams action. Scientifically speaking, the question of source term estimation (STE) is a challenging inverse problem, due to its ill-posed nature. First, there is the issue of the non-uniqueness of the source reconstruction solution: for a given set of concentration measurements obtained by the sensors, there may be several possible

H. Rajaona et al. / Atmospheric Environment 122 (2015) 748e762

source terms that would fit the observations. This problem is aggravated by the fact that the observations could occasionally be non-representative of the reality, because of the errors introduced in the acquisition process such as the sensor noise, or the potential averaging of the measurements done by the receptors. Next comes the fact that the process of estimating the source may easily be jeopardized, due to the physical mechanisms it relies on: the micrometeorological parameters involve complex processes which may provide quite different outputs in certain conditions, should those parameters vary in small or large amplitudes. To solve the problem of STE, several methods have been developed, using different approaches. The most physically intuitive one is the use of inverse transport, which consists in running an atmospheric dispersion model backward in time, from the observation times to the release times (Robertson and Langner, 1998; Pudykiewicz, 1998). The backward model relies on an adjoint transport equation, which is derived by using the principle of timesymmetry employed in atmospheric transport (Hourdin and Talagrand, 2006): switching from the forward model requires changing the sign of each term but the diffusion in the transport equations. Here, the ill-posed aspect of the STE problem is treated by resorting to a regularization process, allowing the construction of a well-posed unique solution by minimizing a cost function. The idea behind this is to penalize the non-desired solutions of the STE problem by applying constraints on them, and to perform a goodness-of-fit measure between the predicted concentrations given by the model and the real observations. There are several ways to optimize the cost function: one approach is to use genetic algorithms, based on evolutionary computation, that runs a selection-and-mutation process to a population of candidates in order to provide the optimal solution (see Haupt et al., 2007; Allen et al., 2007; Cervone et al., 2010; Cervone and Franzese, 2011; Annunzio et al., 2012 for examples of genetic algorithms used in the STE context). Several alternatives are available and have been studied, such as the concept of illumination (Issartel, 2005), or principles relying on information theory, such as the maximum entropy on the mean (Bocquet, 2005). Inverse modelling can also be viewed as a Bayesian problem (Winiarek et al., 2011) where the expression of the cost function is established using Bayesian inference, introducing probabilistic considerations in the problem in order to obtain a point estimate of the source parameters. One Bayesian inference strength is indeed to be able to estimate the uncertainty factor regarding the estimation. Another way of exploiting the Bayesian approach consists in seeking not just for one optimal solution, but obtaining the probability distribution of the source characteristics using the Bayesian paradigm: in this case, the source is defined by a set of parameters, which are the quantities of interest. By the means of stochastic sampling, the posterior probability distribution of these parameters is evaluated in order to fully describe the parameters of the source and the uncertainty on them. The goal of STE is then to look for the most probable parameters for the source in terms of posterior probability. The most frequently used tool for sampling from this posterior distribution in STE problems is the family of Markov Chain Monte Carlo (MCMC) algorithms. In Keats et al., 2007, MCMC is used to determine the source parameters in a complex urban environment. In Senocak et al., 2008, the MCMC procedure is embedded into a reconstruction tool that also assesses the wind field parameters of the dispersion model. In Chow et al., 2008, MCMC is implemented alongside a Computational Fluid Dynamics (CFD) model that aims at improving the accuracy of the physical model. In Delle Monache et al., 2008, a parallel processing scheme using MCMC is built in order to reconstruct the source term of the Algeciras incident. In Yee et al., 2014, MCMC is used to reconstruct the location and the emission rate of a medical isotope production

749

facility using activity concentration measurements obtained from the International Monitoring System radiological network. Overall, MCMC methods have proven to give good results, and their use has extended to other similar purposes, such as water pollution detection (Hazart et al., 2014) or natural gas seeps mapping (Hirst et al., 2013). The method we introduce in this paper offers an alternative way of tackle this Bayesian inference problem, by developing an adaptive scheme based on the principle of Importance Sampling (IS). This method allows to enhance the estimation process by learning from the previous sampling rounds through a recycling process, thus accelerating the convergence and reconstructing faster the source term parameters. In Section 2, we describe the Bayesian framework and the formalism of our study. In Section 3, we develop the foundations of the methodology around the AMIS algorithm for source location purposes. Finally, in Section 4, we implement our solution to the STE problem using a framework derived from the Fusion Field Trials 2007 experiment and compare its behaviour to the MCMC, before concluding. 2. Bayesian modelling In this paper, we are interested in estimating the unknown characteristics q of a source, given the measurements, d, obtained from all the sensors deployed in the field. More specifically, a Bayesian solution will be designed for the inference in order to take into account all the statistical information given by our problem. As a consequence, the quantity of interest is the posterior distribution given by the Bayes's rule as:

pðqjdÞ ¼

pðdjqÞpðqÞ pðdÞ

(1)

As we can see the main ingredients of Bayesian analysis are both the prior distribution, p(q), and the likelihood distribution, p(djq). The prior distribution represents the beliefs about the unknown state before obtaining any observations. On the other hand, the likelihood distribution gives the probability of obtaining the data given a certain set of parameters values. p(d) represents the marginal likelihood and acts as the normalizing constant of the product of the prior and likelihood distributions in order to obtain the posterior distribution. Let us now specify both the likelihood and the prior distributions for the problem considered in this paper. 2.1. The likelihood model In this work, we consider a point-wise and static source represented by q¼(xs,q) where xs ¼ (xs,ys) is the spatial position of the source and q is the release rate vector resulting from the discretization of the plausible emission time interval into Ts time steps, s , in order to take into account the possible time depenfDtn gTn¼1 dence aspect of the source. Let us assume the concentration is measured at a finite number of time points by a set of Nc sensors deployed in the field. The observed concentration di,j acquired by the i-th sensor ci at location xci and time tj, where i ¼ 1,…,Nc and j ¼ 1,…,Tc with Tc the number of time samples for each sensor, is modelled by:

di;j ¼

Ts X

qn C i;j ðxs ; Dtn Þ þ εi;j

(2)

n¼1

The first term in the right-hand side of this expression corresponds to the mean concentration resulting from the superposition s of the Ts releases on the different time steps fDtn gTn¼1 weighted by s the actual emission rates fqn gTn¼1 of the source. Ci,j(xs,Dtn)

750

H. Rajaona et al. / Atmospheric Environment 122 (2015) 748e762

represents therefore the mean concentration observed by the i-th sensor located at xci and at time tj if an unitary release is made during the time step Dtn at location xs. Finally, εi,j corresponds to a random noise vector encompassing the three classical types of error: 1. the model error, related to the underlying uncertainty of the physical model used to compute Ci,j(xs,Dtn), 2. the observation error, inherent to the sensors and the surrounding instrumental environment, 3. the representativeness error, due to the interpolation processes in time and space implied in the data acquisition scheme. It underlines the fact that a dispersion model may be effective at coarser scales, but has more difficulties taking into account subgrid events such as the turbulence effects. The representativeness error is often considered part of the model error, but in fields such as data assimilation it is formally incorporated into the observation error (Koohkan and Bocquet, 2012). As mentioned in Yee, 2008, the choice of a Gaussian noise is justified by bringing forward the argument of the maximum entropy principle (Jaynes, 2003), which stipulates that such an assumption represents a maximally uninformative state of knowledge. Since the combination of possible errors is quite complex, it is difficult to make any other statistical statement about its true sampling distribution than the finite nature of its variance, hence the choice of the Gaussian distribution. Equation (2) can be rewritten in a matrix form in order to take into account all the collection of available measurements:

2

C 1;1 ðxs ; Dt1 Þ 6 « 6 6 C 1;T ðxs ; Dt1 Þ c 6 6 « 6 6 C N ;1 ðxs ; Dt1 Þ c 6 4 « C Nc ;Tc ðxs ; Dt1 Þ

C 1;1 ðxs ; Dt2 Þ « C 1;Tc ðxs ; Dt2 Þ « C Nc ;1 ðxs ; Dt2 Þ « C Nc ;Tc ðxs ; Dt2 Þ



(3) T

where: d ¼ d1;1 … d1;Tc … dNc ;1 … dNc ;Tc is the vector of observed concentration values given by all the network of sensors, and C(xs) takes the following matrix form, also detailed in Fig. 1:

… … …

  3 C 1;1 xs ; DtTs 7 « 7 C 1;Tc xs ; DtTs 7 7 7 « 7 C Nc ;1 xs ; DtTs 7 7 5 «  C Nc ;Tc xs ; DtTs

(4)

In Seibert and Frank, 2004 and Yee, 2009, such operators are called source-receptor matrices. The computation of C is an important part in a STE process: it is the way to quantify the predicted concentration field given by a dispersion model for a given description of the source. As described in Rao, 2007, there are two ways of computing the source-receptor matrix: either using a forward modelling approach, or a backward modelling approach. Winiarek et al., 2011 emphasizes that the forward approach yields a column-by-column computation of the matrix, whereas the backward approach yields a line-by-line computation. In our context, we choose the forward approach: C is hence an operator that takes a source position xs as input and computes the resulting concentration field measured at the Nc sensors over Tc time steps for a number of Ts instantaneous unitary releases. In consequence, C is a matrix of NcTc rows and Ts columns. As in (Yee, 2009), the error term ε is considered to be a spatially and temporally independent zero-mean Gaussian multivariate random variable, i.e. ε  N ð0; s2ε ITc Nc Þ, thus yielding to the following likelihood distribution:

pðdjqÞ ¼

Nc Y Tc  Nc Y Tc   Y   Y p di;j q ¼ N di;j ; C i;j ðxs Þq; s2ε i¼1 j¼1

d ¼ Cðxs Þq þ ε

…

(5)

i¼1 j¼1

with:

 C i;j ðxs Þ ¼ C i;j ðxs ; Dt1 Þ

…

  C i;j xs ; DtTs

(6)

where the notation N ðdi;j ; C i;j ðxs Þq; s2ε Þ represents the value obtained by evaluating the density function of the normal distribution of parameters ðC i;j ðxs Þq; s2ε Þ at point di,j. 2.2. Specification of the prior distributions Specifying the prior consists in representing our belief about the unknown state q using probability distributions before obtaining any observations. In this work, we consider that a source can be anywhere uniformly in some spatial surveillance area, denoted by U3ℝ2 , i.e.:

pðxs Þ ¼ U U ðxs Þ

(7)

Depending on the scenario, more informative distribution could be used in order to take into account that the source is more likely to be on some specific area (nuclear plants, industrial sites, etc.). As in Winiarek et al., 2011, we use a multivariate Gaussian distribution as prior for the emission rate vector:

  pðqÞ ¼ N q; mq ; Sq

Fig. 1. Illustration of the relation between the source's emission rate (in orange, top) and the source-receptor matrix (in blue, bottom) for a given release. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

(8)

Bocquet, 2008 points out that this is a crude approximation as the source emission rate vector should not have negative values, yet this simple assumption is often used in atmospheric dispersion inverse modelling with reasonable performances (Issartel and Baverel, 2003). However, taking into account the gaussianity of the source emission rates could provide valuable information and lead to better performance results. For example, if the release is known a priori to have a low emission rate, then the value of the

H. Rajaona et al. / Atmospheric Environment 122 (2015) 748e762

covariance Sq would be set according to small values. As a consequence, to ensure the coherence of the results while conserving the Gaussian framework, which will be of primary importance in our computation scheme, we propose to enforce a positivity constraint in the inference procedure as detailed in Section 2.3.2. 2.3. Bayesian inference As discussed above, the aim is to provide an approximation of the posterior distribution p(qjd) in order to derive all the possible quantities of interest such as point estimates or confidence intervals. Unfortunately, due to the non-linear dependency of xs in the likelihood through the computation of the source-receptor matrix C(xs), the analytic computation of p(qjd) is intractable. However, the posterior of interest can be expanded as follows:

pðqjdÞ ¼ pðxs ; qjdÞ ¼ pðqjxs ; dÞpðxs jdÞ

(9)

Owing to the Gaussian assumption of both the likelihood in Equation (5) and the prior distribution of q in Equation (8), the rule of conjugate priors states that the conditional posterior of the source emission rate p(qjxs,d) can be evaluated analytically unlike the second term p(xsjd) which remains intractable. By using such decomposition, instead of having to approximate the full posterior distribution p(xs,qjd), only the posterior marginal distribution p(xsjd) needs to be approximated since an analytical expression for p(qjxs,d) can be obtained. Using Bayes' rule in the same form as in Equation (1), the posterior marginal distribution p(xsjd) can be written as follows:

pðdjxs Þpðxs Þ pðxs jdÞ ¼ pðdÞ

given by:

pðdjxs Þ ¼ N



d; Cðxs Þmq ; Cðxs ÞSq Cðxs ÞT þ s2ε ITc Nc



(12)

Because of the complexity of this likelihood function, we develop in the next section an efficient sampling algorithm based on recent advances in Monte-Carlo methodologies which is designed for approximating the posterior marginal distribution p(xsjd). Let us now discuss the analytical derivation of the conditional posterior distribution of the emission rate vector p(qjxs,d). 2.3.1. Conditional posterior of the emission rate As mentioned previously, the conditional posterior p(qjxs,d) can be derived analytically since both the likelihood in Equation (5) and the prior of q are assumed to be Gaussian multivariate distributions. By using the properties inherent to Gaussian distributions and more especially its conditioning and conjugate prior properties (Gelman et al., 2003), this conditional posterior distribution of q is also a Gaussian distribution:

  ~q pðqjxs ; dÞ ¼ N q; m ~q; S

(13)

with parameters given by:

  ~ q ¼ mq þ K d  Cðxs Þmq m ~ q ¼ Sq  KCðxs ÞSq S

(14)

with:

 1 K ¼ Sq Cðxs ÞT Cðxs ÞSq Cðxs ÞT þ s2ε ITc Nc

(15)

(10)

The denominator, often called evidence, or marginal likelihood, represents a measurement of the suitability of the data model to the actual measurements. Although for some specific problems (such as model selection) this evidence is a quantity of interest, this term just represents in our problem a normalizing constant. Moreover, with Monte-Carlo methods, samples from the posterior distribution in Equation (10) can be obtained by using only its formulation up to a normalizing constant, i.e.:

pðxs jdÞfpðdjxs Þpðxs Þ

751

(11)

In Equation (11), since we consider a Gaussian noise, the marginal likelihood p(djxs) is a multivariate normal distribution

2.3.2. Enforcing the positivity constraint on the emission rate Taking into account constraints on the unknown parameters leads generally to an accuracy improvement of the final estimates. In order to incorporate the inequality constraints on the emission rate vector, i.e. c n ¼ 1; …; Ts ; qn  0, we propose to use a technique originally proposed in Simon and Simon, 2010, based on an iterative procedure that will enforce each of these multiple constraints to be satisfied successively in order to make this problem tractable. Each of the constraints is imposed by using a projection method on the constrained space of the original distribution (see Fig. 2 for an example). The algorithm detailed in the seminal paper (Simon and Simon, 2010) is summarized for the problem we consider in this work in Algorithm (1).

Fig. 2. Illustration of the positivity constraint procedure on a bivariate distribution in 3 different cases (a), (b) and (c): plot of the covariance ellipsoid that encapsulates 90% of the probability mass before (red) and after (black) the procedure. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

752

H. Rajaona et al. / Atmospheric Environment 122 (2015) 748e762

After the projection procedure, the conditional posterior distribution is approximated as:

pðqjxs ; dÞzpc ðqjxs ; dÞ ¼ N



c

~ q; m ~ cq ; S q



(16)

Consequently, if such procedure is used, the marginal likelihood p(djxs) is approximated in the same way, and can therefore be expanded by deriving a similar form of Equation (9):

pðdjxs Þzpc ðdjxs Þ ¼

pðdjq; xs ÞpðqÞ pc ðqjxs ; dÞ

(17)

Resorting to a prior distribution supported on semi-finite interval [0;þ∞) such as the log-prior would be an alternative solution, but at the expense of additional computational cost. Indeed, without the analytical solution for the conditional posterior, q must be added to the parameters estimated with the Monte-Carlo sampling method since its conditional posterior distribution given in Equation (13) could not be derived in an analytical form. Let us now describe the proposed Monte-Carlo algorithm in order to obtain the second term on the right-hand side of Equation (9).

3. Estimating the source location with the AMIS algorithm As mentioned previously, we need to approximate the marginal posterior distribution p(xsjd) since no closed-form expression can be derived. Classical Monte Carlo approaches require the ability to sample from the posterior distribution: this task is far from being easy because the distribution of interest does not have a form that allows direct sampling from it. The alternative is then to use more advanced Monte Carlo algorithms. One of the most frequent Monte Carlo algorithms encountered in similar problems is the Markov Chain Monte Carlo (MCMC) method (e.g. Yee et al., 2014). The main idea of MCMC is to create, by using some proposal distribution 4(,), a Markov chain that has the desired distribution as its stationary distribution. There are numerous MCMC algorithms in the literature, the most common ones being the MetropoliseHastings sampler and the Gibbs sampler (Robert and Casella, 2004). Most MCMC algorithms have to get through a transition regime at the beginning of their execution called the burn-in time, during which the generated samples do not qualify as usable results in the convergence process. In this paper, for comparison purposes with our algorithm, we use a simple MetropoliseHastings algorithm, with a transition kernel based on a random walk. Even though there are several ways to build a more efficient sampler, we choose a basic MCMC implementation in the objective of assessing how

H. Rajaona et al. / Atmospheric Environment 122 (2015) 748e762

our own algorithm fits in the Bayesian framework of STE problem, instead of proving its superiority to MCMC methods. The method we develop in this paper does not rely on the evolution of a stochastic Markov process as MCMC does: it leans towards an adaptive way to perform another kind of sampling within the Monte-Carlo framework called the Importance Sampling. Rather than resorting to an acceptationereject procedure as it is done in the MetropoliseHastings algorithm, Importance Sampling uses all the sampled items in a weighting scheme. 3.1. The importance sampling approach

 

p xðiÞ ci2 1; …; Np wðiÞ ¼  ðiÞ  4 x

(18)

By using the Importance Sampling identity (Robert and Casella, 2004), the target distribution is approximated by the following empirical measure:

pðxÞz

Np X

drawn at the previous iterations by using a recycling scheme to improve both the learning of the proposal and the accuracy of the approximation of the target distribution (in the PMC only the particles form the last iteration are used). Using every single simulated particle at every step of this iterative algorithm could therefore lead to a significant improvement as shown in Cornuet et al., 2012. The AMIS algorithm is used to approximate the posterior marginal distribution of the source location. The target distribution is thus defined as:

pðxs Þfpðdjxs Þpðxs Þ

The goal of Importance Sampling (IS) is to generate samples, also called particles, from a distribution p that is impossible to simulate directly, and is referred to as the target distribution, with the help of another distribution 4, called the proposal distribution (Robert and Casella, 2004). Basically, Np particles ðxð1Þ ; …; xðNp Þ Þ are sampled from 4 and then, since the proposal used to draw the samples is not the target distribution, a correction is performed by affecting to each particle an importance weight defined as follows:

753

(20)

where the prior distribution is defined in Equation (7) and the marginal likelihood is given in Equation (12) or in Equation (17) if the positive constraint of the emission rate is taken into account as proposed in Section 2.3.2. The AMIS algorithm requires to choose a parametric proposal distribution that would be flexible enough and easy to sample
et al., 2008, we from: following the idea presented in Cappe consider a mixture of D distributions. Since we want to estimate a 2-dimensions position, the proposal is a mixture of bi-variate Gaussian distributions:

4ðxs ja; m; SÞ ¼

D X

    aðdÞ 4ðdÞ xs mðdÞ ; SðdÞ

(21)

d¼1

where:

~ ðiÞ dxðiÞ ðdxÞ w

(19)

i¼1

~ ðiÞ are the normalized importance weights (i.e. where the w PNp ðjÞ 1 ðiÞ ðiÞ ~ ¼ w ½ j¼1 w  ) and dð,Þ is the Dirac function. IS can be used w to draw samples from both univariate and multivariate distributions. However, the variability of the importance weights may strongly affect the accuracy of the estimate. In the case that these weights vary widely, only a few points with the highest weights effectively contribute to the estimate. This problem causes the subsequent estimations to be inaccurate. This situation may arise when the importance function generates samples that are concentrated on regions with low probability mass regarding the target posterior. In such cases, the IS approximation is not very robust owing to the high discrepancy between the target density and the importance density. Therefore, for importance sampling to work well, the proposal distribution defined by 4(x) should be chosen as close as possible to the target distribution p(x), so that the ratio p(x)/4(x) does not vary too much. 3.2. The AMIS algorithm: an adaptive method Some adaptive importance sampling algorithms, referred to as Population Monte Carlo algorithms (PMC) in the seminal paper of
et al., 2004, have been proposed with the ability to proCappe gressively learn about the target distribution and thus to adaptively modify the proposal distribution based on the gathered information so that the sampling procedure becomes more efficient. The PMC methods allow to improve the proposal distribution by successive iterations using the formerly produced importance samples and weights. More recently, a novel adaptive IS method, named Adaptive Multiple Importance Sampling (AMIS), has been proposed in Cornuet et al., 2012. The main difference between both approaches is the use in the AMIS of all the collection of particles

P ðdÞ ¼ 1 that
a ¼ fað1Þ ; …; aðDÞ g is a weight vector such as D d¼1 a illustrates the influence of each component of the mixture,
ðm; SÞ ¼ fðmð1Þ ; Sð1Þ Þ; …; ðmðDÞ ; SðDÞ Þg is the set of parameters of all the D components of the mixture (means and covariances) that will be adapted through the iterations. Executing the AMIS algorithm requires an initial tuning of the proposal distribution in order to define how the first round of sampling will be conducted. In our case, since a uniform prior is used for the source location, the proposal is set up so that the generated samples cover all the domain of study in a roughly uniform way. Practically, the components of the Gaussian mixture are initially dispatched along the domain of interest in order to broadly cover it for the first sampling round, with no preference among the D components: c d; aðdÞ ¼ 1=D. The flexibility of our model lies in the fact that if there were any information narrowing down a possible area for the source location, we could adjust this initial step to enhance the importance of the suspected area, by tuning the parameters of the mixture accordingly. One example of application would be the use of a preliminary retro-propagation run to fit the initial proposal distribution, even though in the following applications of our study we did not use any retroprogation results. The k-th iteration of the AMIS is made on the following four main steps: ðiÞ

1. Np particles, xs;k , are generated using the current form of the proposal 4k, each of them representing a potential source location; 2. An importance weight wk is computed and associated to each of these novel particles; 3. All of the previous particles are recycled by using a correction mechanism of their previous importance weights, in order to

754

H. Rajaona et al. / Atmospheric Environment 122 (2015) 748e762

obtain w0:k1; which combined with wk, forms the updated importance vector w0:k, 4. The parameters ða; m; SÞ of the proposal distribution are updated by using all the collection of weighted particles i¼1;…;Np ~ nðiÞ ; xðiÞ fw s;n gn¼0;…;k

 pðxs Þ pðxs Þdxs log 4ðxs Þ

Z D KL ðp k 4Þ ¼

ðdÞ

Np k X X

Raw Averaged

0.8 0.7 0.6 0.5 0.4 0.3 0.2

(22)

When the proposal distribution is chosen as a mixture of Gaussian distributions, the updating mechanism described in
et al., 2008 is given at the end of the k-iteration and for Cappe d ¼ 1,…,D, by:

akþ1 ¼

Raw and averaged concentration measured by sensor 76

x 10

0.9

Concentration (kg/m )

In order to perform this adaptation procedure, it is necessary to
et al., 2008, the define a criterion to optimize. According to Cappe most popular choice is the so-called KullbackeLeibler (KL) divergence, where the objective is to get the set of parameters which minimizes the dissimilarity between the two distribution probabilities that are the target and the proposal distributions. The KL divergence, also called relative entropy, is a way to measure this dissimilarity and is defined by:

1

0.1 0 0

1

2

3 t (1/50s)

4

5

6 x 10

Fig. 3. Results of the averaging process over the raw concentration data of FFT-07's trial 7.

~ nðiÞ rðdÞ w n;i

n¼1 i¼1 ðdÞ mkþ1

Pk ¼

Pk ðdÞ Skþ1

PNp

~ ðiÞ rðdÞ xðiÞ w i¼1 n n;i s;n ðdÞ akþ1

n¼1

   ðiÞ ðdÞ ðiÞ ðdÞ T ~ ðiÞ rðdÞ xs;n  mkþ1 xs;n  mkþ1 w i¼1 n n;i ðdÞ

4.4403

akþ1

1 11

ðdÞ

   ðdÞ ðdÞ ðiÞ  ðdÞ ðdÞ ak 4k xs;n mk ; Sk    ¼P ðmÞ ðmÞ ðiÞ  ðmÞ ðmÞ D xs;n mk ; Sk m¼1 ak 4k

4.4402 y (m)

where rn;i is an intermediary value corresponding to the probability of a given particle belonging to the d-th component of the mixture and is given by: ðdÞ rn;i

Structure of the FFT07 sensor network

x 10

4.4404

PNp

n¼1

¼

4.4405

(23)

2 12 21 31

3 13 22 32 41

4.4401

51

(24)

14 23 33 42 52 61

4.44

71

15 24 34 43 53 62 72 81

4.4399

The overall procedure of the AMIS is summarized in Algorithm (2). The collection of weighted particles obtained as output of the AMIS algorithm therefore approximates the posterior marginal distribution of the source location as follows:

pðxs jdÞz

Np K X X n¼1 i¼1

~ nðiÞ d ðiÞ ðdxs Þ w x

Np K X X n¼1 i¼1

3.314

3.315

16 25 35 44 54 63 73 82 92

17 26 36 45 55 64 74 83 93

3.316

19

18

28

27

38

37

47

46

57

56

66

65 75 84 94

76 85 95

3.317 x (m)

10 20 29 39 48 58 67 77 86 96

30 40 49 59 68 78 87 97

50 60 69 79 88 98

70 80 89 99

90 100

Unused sensor Subdomain sensor Source (trial 7) Source (trial 30) 3.318

3.319

3.32 x 10

Fig. 4. The FFT-07 digPIDs layout: the filled dots represent the 25 sensors at the vicinity of the source which form the network we effectively use for our calculations.

(25)

s;n

where K is the total number of iterations. Finally, by plugging this empirical measure into Equation (9), the main posterior of interest, including both the position and the release-rate parameters, is thus approximated by:

pðxs ; qjdÞz

91

4.4398 3.313

6

5

4

9

8

7

   ðiÞ ~ nðiÞ p qxs;n ; d dxðiÞ ðdxs Þ w s;n

(26)

4. Results and discussions In order to validate the algorithms described in Sections 2 and 3, we apply them to the FFT-07 experiment, a fairly used dataset in the field of STE. The following paragraphs explain what FFT-07 is about, and how we use it to test our method.

H. Rajaona et al. / Atmospheric Environment 122 (2015) 748e762

4.1. The FFT-07 experiment The FUSION Field Trials 2007 (FFT-07) experiment was designed by the American Defense Threat Reduction Agency (DTRA) with the aim of providing a comprehensive dataset of meteorological and concentration measurements to research teams so that they may be able to run their STE algorithms using real field data (Platt and Deriggi, 2010; Annunzio et al., 2012; Singh and Rani, 2014). FFT07 was designed for the study of short-range releases: the covered domain is a 500 m  500 m square. Several release configurations were used, each situation being classified as a trial with properties taking into account the time of day or night, the wind speed and direction, the atmospheric stability, and the number of sources. The data acquisition was performed by a network of 100 digital

755

photoionization detectors (digiPIDs), dispatched over a uniform grid on the domain, with a distance of 50 m between each other, at 2 m from the ground. The recorded concentrations were acquired at a high frequency (50 Hz). In our study, given the considerable size of the raw output given by the sensors, we reduce the size of the observation vector by running a mobile average over the raw observations, and then subsampling it to get one value every 10 s. The size of the resulting observation vector is then smaller, and yet the information given by the original signal remains, as shown by Fig. 3. In the FFT-07 terminology there are two kinds of releases: continuous releases last several minutes, whereas instantaneous releases consist in one single puff. For our calculations involving experimental observations, we use the data set of FFT-07's trial 7, which falls under the continuous type. During trial 7, a single source released continuously a tracer gas for a period of 10 min (8:38 AM

756

H. Rajaona et al. / Atmospheric Environment 122 (2015) 748e762

to 8:48 AM). The digiPIDs were measuring concentrations during a period of 19 min (8:36 AM to 8:55). Originally, the FFT-07 was made of 100 digiPID sensors. However in practical cases the domain of interest is not that much instrumented, and the number of sensors is not as numerous as in these experiments. Nonetheless, we want to work with a sufficient and reliable amount of observation data for our validation process. This is why we choose to focus on a reduced number of sensors, by taking into consideration only the nearest neighbours of the digiPID which record the highest concentration. Practically, this leads to consider a subgrid of 25 sensors, as showed in Fig. 4: we refer to the portion of the domain where these sensors are located as the subdomain. The generality of our method enables the use of any available dispersion model: in this paper we choose to work with a Gaussian puff model to generate synthetic data as well as to compute the source-receptor matrix C. The choice of a puff model instead of a plume model is justified by the need of the time-dependency aspect, if we are to estimate time-varying emission rates. It is also appropriate when meteorological parameters such as the wind speed and direction are not constant. The sensors and the source are on the same altitude level, close to the ground (2 m): since the vertical wind has sufficiently small values during the observation time, we simplify the z-dependent exponential term in the classical Gaussian puff formula, which becomes a scalar factor of 2 that illustrates the total reflection of the emitted puffs on the ground. Following Equation 2, each element of the matrix C represents the result of a unitary release at a given time step. The concentration measured at the i-th sensor xci ¼ ðxci ; yci Þ, at the observation time tj, produced by a source located at xs, and emitting at the time step Dtn is hence given by:

C i;tj ðxs ; Dtn Þ ¼

 1 2 2 L exp  þ L x y 2 ð2pÞ3=2 sx sy sz 2Qu Dtp

where:
Qu is the unitary mass released: since the concentration observations are given in kg/m3, we have Qu ¼ 1 kg/m3,
Dtp is the discretization time step of the puff model, namely the elapsed time between the emission of two consecutive puffs. For a good accuracy we choose Dtp ¼ 1 s,
sx, sy and sz are the dispersion coefficients,
Lx and Ly are given by:

0 xci  @xs þ

tj X

1 ux ðtÞ  dw A

t¼Dtn

Lx ¼

sx

0 yci  @ys þ

tj X

1 vx ðtÞ  dw A

Ly ¼

sy

dw is the time interval between two observed values of the wind vectors ux and vx. We consider the case of an isotropic diffusion in x and y so we have sx¼sy. The values of sy and sz are computed using the following power law: ys estimation 0.4

0.35

0.35

0.3

0.3 p(ys|d)

p(xs|d)

0.25 0.2

0.25 0.2

0.15

0.15

0.1

0.1

0.05

0.05

0

331.4

331.6 X (km)

331.8

0 4439.8

332

4440

xs estimation

4440.2 Y (km)

4440.4

ys estimation 0.45

0.3

0.4 0.35

0.25

0.3 p(ys|d)

p(xs|d)

0.2 0.15

0.25 0.2 0.15

0.1

0.1 0.05 0

0.05 331.4

331.6 X (km)

331.8

(28)

t¼Dtn

xs estimation

0.4

(27)

0 4439.8

332 ð1Þ

Fig. 5. Estimation of the posterior distributions for sources located at xs

4440 ð2Þ

(top) and xs

4440.2 Y (km)

4440.4

(bottom) with synthetic concentration data.

H. Rajaona et al. / Atmospheric Environment 122 (2015) 748e762

3

Q estimation

x 10

2.5

2

Q

1.5

1

0.5

0

−0.5

−1

0

100

200

300

400 t (s)

500

600

700

800

Fig. 6. Synthetic concentrations: comparison of the estimation of the release rate q (in kg s1), if accounting for a positivity constraint (blue) or not (green) with the true values (dashed red). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

s ¼ axbs

(29)

(a,b) are empirical coefficients which depend on the atmospheric stability class (Pasquill and Smith, 1983). Their values differ regarding whether sy or sz is evaluated. xs is the distance between the current puff's center of mass and the source. Here we use an open-country Pasquill D-type stability. 4.2. Results for synthetic data using FFT-07 framework In this first step, in order to validate our method, we choose to use a synthetic framework, where the positions of the sensors and the true sources are the same as in FFT-07 trials 07 and 30, but the observations are obtained via simulations. The meteorological data is optimized so that the synthetic concentrations fit the real FFT-07 concentration data, and the wind is considered to be constant in

757

direction and in speed. We choose to work with a number of D ¼ 4 components in the Gaussian mixture of the proposal distribution. The AMIS algorithm is run for K ¼ 10 iterations, with Np ¼ 100 particles sampled at each iteration. Synthetic data are generated with a noise value of s2ε ¼ 105 to provide a realistic observation vector, and the prior parameters are defined by mq ¼ 0 and Sq ¼ s2q ITs , where s2q ¼ 103 . A smaller value of s2q means that the source is a priori close to the sensor, whereas a greater value would mean a farther potential source: the value we choose offers a compromise between both alternatives. The algorithm was tested with two different source positions, both upwind regarding the sensor grid, located respectively at ð1Þ ð2Þ xs ¼ ð331:759; 4439:960Þ and xs ¼ ð331:825; 4439:972Þ, which are inspired from FFT-07's setup for respective trials 7 and 30. These coordinates are measured in kilometres, as for all of the coordinate references in the rest of this paper. We can see in Fig. 5 that the AMIS is able to make a correct estimation of the source position with a good accuracy in both cases. The representation of the estimation using a probability distribution function provides an explicit representation of the uncertainty around the estimation, emphasizing how useful the Bayesian inference framework is when compared to results provided by other methods such as optimization-based ones in which only a point-estimate is provided. It is also possible to compute a point estimate of the location by considering the minimal mean i¼1;…;Np ~ nðiÞ ; xðiÞ square estimator (MMSE). If fw s;n gn¼0;…;K is the final output of the AMIS algorithm, the MMSE estimator b x s is obtained from Equation (25):

b xs ¼

Np K X X

ðiÞ ~ nðiÞ xs;n w

(30)

n¼0 i¼1

It appears that the inferred source is located less than 1 m far from the true source, for both situations. Fig. 6 illustrates the two different possibilities of estimating the Table 1 Statistical evaluation of the estimations given by 40 runs of the AMIS algorithm with and without positivity constraint (P.C.), where dbx ;x is the distance to the actual s s source. P.C. Yes No

dbx ;x ðmÞ s s 0 10.44

sbx ðmÞ s 1

9  10 4  101

sby ðmÞ s 3.2 1.0

Fig. 7. Synthetic case: estimation of q in kg s1 (blue) and its uncertainty interval þ=  2~ sq ; compared to the true value (dashed red), with a simple release (a), a variable release (b) and a continuous release (c).

758

H. Rajaona et al. / Atmospheric Environment 122 (2015) 748e762

Table 2 Statistical evaluation of the estimations given by 40 runs of the AMIS and MCMC algorithm, both with positivity constraint, for a real source located at ð1Þ xs ¼ ð331:759; 4439:966Þ. Algo.

b xs

sbx ðmÞ s

sby ðmÞ s

AMIS MCMC

(331.759,4439.966) (331.742,4439.981)

0.9 130

3.2 95

release-rate vector q, by resorting or not to the positivity constraint (Equation (16) or Equation (13) with xs ¼ b x s ). The output values are compared with the original release rate: the result that yields a positivity constraint is more accurate and much closer to the theoretical value of q, with respect to the amplitude and the time intervals. In practice, the profile of the source emission can take a nontrivial form, especially when the release is not restricted to a single time interval. The reconstruction algorithm was hence tested on a case where there is a slight variation during the emission time when generating the synthetic data, with the use of the positivity constraint. As illustrated in Fig. 7, our method manages to deal efficiently with a non-uniform release rate, as well as for a constant release, assessing the fact that the algorithm is able to perform correctly with non-trivial emission types. Given the fact that q is estimated ~ q , we can also represent an uncertainty alongside its covariance S margin regarding the estimation: in the three cases we can see that the true value of the emission rate remains within the interval of interest. In the case (c) of Fig. 7, the spike at the upper boundary of the estimated profile comes from the fact that since the release is longer than the observation frame, at the time of the spike there are no more information brought by the observations, hence the high value of the uncertainty margin. Since the AMIS algorithm relies on IS, which is stochastic by nature, its validation requires several trial runs in order to assess its robustness, unlike deterministic methods where the estimate is always the same for identical input parameters. Consequently, the AMIS was run multiple times, with different observation data at each run, with and without positivity constraint, to compare the overall performance of the two options. The results shown in Table 1 confirm that the positivity

constraint (P.C.) does indeed bring more accuracy to the estimation process of the source location, in addition to enhancing the quality of the reconstruction of the release rate, by giving an exact estimate ð1Þ of the source position xs . When it comes to estimating q, the gain of precision can be measured by comparing the mean square error b and the true value q: (MSE) between the estimated vector q

b ; qÞ ¼ MSEð q

Ts 1 X b ðDti Þ  qðDti ÞÞ2 ðq Ts i¼1

(31)

The case where the positivity constraint is applied generates a smaller error, with a MSE value of 1.37  1012, compared to the case without positivity constraint, which yields a MSE of 1.26  109. The performance of the AMIS was next compared to another popular Bayesian inference technique for STE, namely the MCMC. To do so, we implemented a MCMC algorithm using the MetropoliseHastings sampler (Metropolis et al., 1953) with a fixed Gaussian kernel transition. We put it up against the AMIS in a robustness test over multiple runs, by changing the starting point of the Markov chain at each iteration, all starting points being sampled randomly across the subdomain and covering the same area as the initial particles of the AMIS which cover uniformly the subdomain. The MCMC algorithm is implemented with the same positivity constraint as the AMIS, and runs for KMCMC ¼ 1000 iterations. While using MCMC algorithms, it is usual to ignore the first batch of results obtained during the so-called burn-in period, where the outputs are not relevant: in our case, we cast out the first 100 iterations. As shown in Table 2, the AMIS performs better, in terms of accuracy and spread of the estimation results. The MCMC is penalized by the random choice of its starting point: the farther it is from the true source, the longer it takes to converge, and the estimation error may sometimes be important, as highlighted by the presence of several outliers revealed by the values of the standard deviations for b x s and b y s . The fact to be noted here is that the AMIS adapts better to different initialization scenarios: sampling particles is in a way more flexible than relying on a single starting point. One can fix this flaw of the MCMC by running several Markov chains with different starting points, but at a cost of increasing the necessary computational resources.

10

80

MCMC AMIS

9

70

8

Effective Sampling Size (ESS)

Distance to the real source (m)

60

7 6 5 4 3

50

40

30

20

2 10

1 0

1

2

3

4

5 6 Iterations (x100 for MCMC)

7

8

9

10

Fig. 8. Comparison of the convergence speed for the MCMC and the AMIS when estimating the distance between the true source and its estimation.

0

1

2

3

4

5

6

7

8

9

10

Iterations

Fig. 9. Evolution of the ESS related to the population of Np ¼ 100 particles of the AMIS in the synthetic case.

H. Rajaona et al. / Atmospheric Environment 122 (2015) 748e762

759

distribution and the initial sampling round are similarly crucial parameters to be tuned: convergence issues can be raised if the space covered by the particles isn't wide enough to encompass the true source unknown location. When it comes to assessing the convergence, the AMIS differs in its nature from the MCMC in the fact that the samples it produces come directly from the target distribution. The MCMC relies on the convergence of a stochastic process, namely a Markov chain, to a stationary distribution, thus the need for rigorous convergence tests. For algorithms depending on IS such as AMIS, an alternative tool is the Effective Sampling Size (ESS), using the importance weights. It is particularly useful to look out for degeneracy problems, when all but a very few weights have a zero value. The ESS is defined by:

PNp i¼1

ESS ¼

PNp i¼1

Fig. 10. Evolution of the Euclidian distance between two consecutive estimations of the position.

We also compare the convergence speed of the AMIS and the MCMC. In Fig. 8, the estimation results are computed at each iteration for the AMIS, and every 100 iterations for the MCMC, using the MMSE estimator. In our case, this allows to scale the quality of the estimation by the total number of computations of the likelihood function, which is at the k-th iteration, kNp for the AMIS and simply k for the MCMC. By choosing Np ¼ 100 particles per iteration for the AMIS, K ¼ 10 iterations, and KMCMC ¼ 1000 iterations for the MCMC, the total amount of computation remains the same between the two algorithms. The results illustrated in Fig. 8 show a better convergence rate of the AMIS, underlining the advantage of an adaptive process which allows to provide a more accurate estimation quickly. The convergence speed of the MCMC algorithm relies on how the kernel transition is chosen, especially on the value of its variance for the Gaussian case: a too small value may cause a much longer convergence of the Markov chain, whereas an excessive value can lead to bad mixing, where the steps taken by the chain are too wide and the estimation process cannot be achieved with sufficient accuracy. With the AMIS, the choice of the proposal

175

Wind direction (degrees)

170

!2 wi w2i

(32)

The ESS measures the efficiency of the samples at representing the distribution of interest and represents therefore an empirical measure on how the proposal distribution is close to the posterior distribution. The maximum value of the ESS is Np (i.e. the number of particles) and is obtained when the proposal used in the algorithm corresponds exactly to the posterior distribution which is the ideal case. As shown in Fig. 9, this ESS increases with the successive iterations. This result empirically shows that the proposal is improving over the successive iterations of the adaptive procedure of the AMIS. When it comes to checking if the algorithm has reached an optimal estimation, one criterion would be to compute a distance between estimates given by the particles and weights at two consecutive iterations, such as in Fig. 10, and to decide that the estimation is good enough when this distance reaches an acceptable threshold, which would be the case here starting from the 8th iteration. 4.3. Results for real concentration data from the FFT-07 experiment Our methodology was also tested for source location estimation using real concentration measurements obtained from the FFT-07 experiment, in order to assess the performance of the AMIS algorithm while confronted to real experimental data. This was accomplished by substituting synthetic concentrations with real ones in the previously described scenario, for the trials where the data sets were available to us. In case of missing data in the raw files, we have to deal with two kinds of issues:
if there is partially missing data, typically when there is a gap of one or several corrupted or missing values between non-zero measurements, we perform a simple interpolation in order to fill these missing values,
in the case of a sensor giving nothing but corrupted or missing data during all the observation time frame, we leave it out of the observation data.

165

160

155

150

145

0

200

400

600 t (s)

800

1000

1200

Fig. 11. Wind direction measured by the wind sensor closest to the source (trial 7).

Using experimental data also allow to see how the estimation algorithm copes with factors such as non-stationary wind, as it is the case with trial 7, where there is a fairly important variation in wind speed and direction, as shown in Fig. 11. Such variations show the need of using an appropriate model such as the puff model, where the wind variations are taken into account in the spatial evolution of the plume. Consequently, for trial 7 we used the data from the two wind sensors which are

760

H. Rajaona et al. / Atmospheric Environment 122 (2015) 748e762

Fig. 12. Estimation of the posterior distribution for the source location (trial 7) using experimental data, compared with the true position (red). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 13. Estimation of the posterior distribution for the source location (trial 30) using experimental data, compared with the true position (red). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 14. Comparison for source position estimation in experimental situation (trial 7) using AMIS (magenta) and MCMC (blue). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

supposedly along the axis of the plume and we averaged the given values over periods of 1 s. The wind variations for trial 30 are less important, so we considered a constant wind for our dispersion model. Similarly to the synthetic case, the representation of the posterior probability distribution function provides the estimate of the position of the source. Figs. 12 and 13 point out a good quality of estimation, almost as close as with the synthetic case. In terms of distance with the true source, the estimate is still within an

acceptable range, with a distance of 9.8 m for trial 7 and 1.5 m for trial 30, emphasizing the fact that the AMIS algorithm is able to perform good source location estimation confronted with real data. The results from trial 7 were obtained with s2ε ¼ 1  1010 and s2q ¼ 5  102 . The results from trial 30 were obtained with s2ε ¼ 5  108 and s2q ¼ 5  103 . For each trial, the choice of s2ε and s2q is made after a trial and error process where we choose the combination which allows the most accurate estimation. As we did before, we also compare the result of the AMIS with

H. Rajaona et al. / Atmospheric Environment 122 (2015) 748e762

761

Fig. 15. Averaged estimation of q (black) using experimental data (trial 7) and its uncertainty interval þ/ 2~ sq (grey), compared to the true value (red), using AMIS (a) and MCMC (b). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

those of the MCMC. For trial 7, using the same parameters mentioned previously for s2ε and s2q , we can see with Fig. 14 that the AMIS performs as well as the MCMC in reconstructing the source's position. The estimation of q for trial 7 is performed on Fig. 15. The values of q are averaged over 1 min windows. The first non-zero values of the estimated emission profile appear earlier and higher than the true values: this is because the estimated position which was obtained by the algorithms (see Fig. 14) is farther upwind compared to the actual position. In order to compensate and produce equivalent concentrations compared to the sensor observations, since there is more distance to cover between the sensors and the estimated source, the estimated release has to start earlier. This phenomenon is slightly more visible for the MCMC case, where the estimated position is a little more upwind than for AMIS. Nonetheless, both methods demonstrates a good estimation of the end time of the release, as well as statistically sound values for the emission rate. Both methods demonstrate a good estimation of the end time of the release, as well as statistically sound values for the emission rate, most of them being in the confidence interval given by the ~ q. covariance matrix S

mislead the algorithm if not properly set up. From an algorithmic perspective, additional work should be considered in narrowing down a convergence criterion from the AMIS, in order to optimize the number of iterations, and thus the total computational load. Regarding the physical aspects of the problem, we realized that a good estimation of the source parameters implies the use of an appropriate dispersion model, which would take into account such phenomena as a fluctuating wind, which impact on the concentration measurements is not to be neglected. That is why we used a puff model capable of dealing with wind variations over time, which proves itself useful in the case of trial 7. However, when it comes to dealing with finer and more complex atmospheric effects, or considering an urban layout, the dispersion model needs to be even more performance. The AMIS algorithm as we use it can easily handle various models, but the computational load shall increase along with the accuracy of the chosen model. Consequently, a potential extension of the work presented in this paper would be to find optimal ways to build the source-receptor matrix using advanced dispersion models, in order to keep a reasonable computation time. References

5. Conclusion In this paper, we presented a Bayesian inference method based on an adaptive IS procedure, the AMIS algorithm, to estimate a single point-wise source of atmospheric release. The posterior distribution of the parameters (xs,q) of the so-called source term given a set of concentration measurements d is reconstructed using a combination of stochastic sampling and analytical results to compute respectively the marginal posterior of the position and the conditional posterior of the emission rate. Our methodology allows us to consider several release types, with potential variations of the release rate in time, or long continuous emissions. Beyond the good results of the estimation yielded by the solution we proposed, the use of an adaptive scheme has proven to increase the effectiveness of the estimation, allowing to reach a comparable accuracy to the existing MCMC approach. Furthermore, the positivity constraint we imposed has shown correct estimations of the release-rate, despite the additional computation load to the overall procedure. The STE algorithm was validated using synthetic and experimental concentration data from the FFT-07 experiment, giving encouraging results. However, the use of the AMIS has also raised some important points that must be carefully taken into account, such as the choice and the initialization of the proposal distribution which can

Allen, C., Young, G., Haupt, S., Apr. 2007. Improving pollutant source characterization by better estimating wind direction with a genetic algorithm. Atmos. Environ. 41 (11), 2283e2289. URL: http://linkinghub.elsevier.com/retrieve/pii/ S1352231006011113. Annunzio, A.J., Young, G.S., Haupt, S.E., Jan. 2012. Utilizing state estimation to determine the source location for a contaminant. Atmos. Environ. 46, 580e589. URL: http://linkinghub.elsevier.com/retrieve/pii/S1352231011004699. Bocquet, M., Jul. 2005. Reconstruction of an atmospheric tracer source using the principle of maximum entropy. I: theory. Q. J. R. Meteorol. Soc. 131 (610), 2191e2208. URL: http://doi.wiley.com/10.1256/qj.04.67. Bocquet, M., 2008. Inverse modelling of atmospheric tracers: non-Gaussian methods and second-order sensitivity analysis. Nonlinear Process. Geophys. 15, 127e143.
, O., Douc, R., Guillin, A., Marin, J.-M., Robert, C.P., Apr. 2008. Adaptive Cappe importance sampling in general mixture classes. Stat. Comput. 18 (4), 447e459. URL: http://arxiv.org/abs/0710.4242.
, O., Marin, J.-M., Robert, C., 2004. Population Monte Carlo. J. Comput. Graph. Cappe Stat. 13 (4). Cervone, G., Franzese, P., Aug. 2011. Non-Darwinian evolution for the source detection of atmospheric releases. Atmos. Environ. 45 (26), 4497e4506. URL: http://linkinghub.elsevier.com/retrieve/pii/S1352231011004432. Cervone, G., Franzese, P., 2010. Machine learning for the source detection of atmospheric emissions. Am. Meteorol. Soc. Chow, F.K., Kosovi
c, B., Chan, S., 2008. Source inversion for contaminant plume dispersion in urban environments using building-resolving simulations. J. Appl. Meteorol. Climatol. 47, 1553e1572. Cornuet, J.-M., Marin, J.-M., Mira, A., Robert, C.P., Jul. 2012. Adaptive multiple importance sampling. Scand. J. Stat. 39 (4), 798e812. URL: http://arxiv.org/abs/ 0907.1254. Delle Monache, L., Lundquist, J.K., Kosovi
c, B., Johannesson, G., Dyer, K.M.,

762

H. Rajaona et al. / Atmospheric Environment 122 (2015) 748e762

Aines, R.D., Chow, F.K., Belles, R.D., Hanley, W.G., Larsen, S.C., Loosmore, G. a., Nitao, J.J., Sugiyama, G. a., Vogt, P.J., Oct. 2008. Bayesian inference and markov chain Monte Carlo sampling to reconstruct a contaminant source on a Continental scale. J. Appl. Meteorol. Climatol. 47 (10), 2600e2613. URL: http:// journals.ametsoc.org/doi/abs/10.1175/2008JAMC1766.1. Gelman, A., Carlin, J., Stern, H., Rubin, D., 2003. Bayesian Data Analysis, second ed. Chapman and Hall/CRC. Haupt, S.E., Young, G.S., Allen, C.T., 2007. A genetic algorithm method to assimilate sensor data for a toxic contaminant release. J. Comput. 2 (6), 85e93. Hazart, A., Giovannelli, J.-F., Dubost, S., Chatellier, L., Mar. 2014. Inverse transport problem of estimating point-like source using a Bayesian parametric method with MCMC. Signal Process. 96, 346e361. URL: http://linkinghub.elsevier.com/ retrieve/pii/S016516841300323X.
lez del Cueto, F., Randell, D., Kosut, O., Aug. 2013. Hirst, B., Jonathan, P., Gonza Locating and quantifying gas emission sources using remotely obtained concentration data. Atmos. Environ. 74, 141e158. URL: http://linkinghub.elsevier. com/retrieve/pii/S1352231013002203. Hourdin, F., Talagrand, O., Jan. 2006. Eulerian backtracking of atmospheric tracers. I: adjoint derivation and parametrization of subgrid-scale transport. Q. J. R. Meteorol. Soc. 132 (615), 567e583. URL: http://doi.wiley.com/10.1256/qj.03.198. A. Issartel, J., 2005. Physics emergence of a tracer source from air concentration measurements, a new strategy for linear assimilation. Atmos. Chem. Phys. 5, 249e273. Issartel, J., Baverel, J., 2003. Inverse transport for the verification of the comprehensive nuclear test ban treaty. Atmos. Chem. Phys. 3, 475e486. Jaynes, E., 2003. Probability Theory: the Logic of Science. Cambridge University Press. Keats, A., Yee, E., Lien, F.-S., Jan. 2007. Bayesian inference for source determination with applications to a complex urban environment. Atmos. Environ. 41 (3), 465e479. URL: http://linkinghub.elsevier.com/retrieve/pii/ S1352231006008703. Koohkan, M.R., Bocquet, M., 2012. Accounting for representativeness errors in the inversion of atmospheric constituent emissions : application to the retrieval of regional carbon monoxide fluxes. Tellus 64, 1e17. URL: http://dx.doi.org/10. 3402/tellusb.v64i0.190471. Metropolis, N., Rosenbluth, A.W., Rosenbluth, M.N., Teller, A.H., Teller, E., 1953. Equation of state calculations by fast computing machines. J. Chem. Phys. 21 (6), 1087. URL: http://scitation.aip.org/content/aip/journal/jcp/21/6/10.1063/1. 1699114. Pasquill, F., Smith, F., 1983. Atmospheric Diffusion, third ed. Wiley. Platt, N., Deriggi, D., 2010. Comparative investigation of source term estimation

algorithms using Fusion Field Trial 2007-Linear regression analysis. In: 13th Conference On Harmonisation Within Dispersion Modelling for Regulatory Purposes. No. June, pp. 901e905. Pudykiewicz, J. a, Sep. 1998. Application of adjoint tracer transport equations for evaluating source parameters. Atmos. Environ. 32 (17), 3039e3050. URL: http:// linkinghub.elsevier.com/retrieve/pii/S1352231097004809. Rao, K.S., Oct. 2007. Source estimation methods for atmospheric dispersion. Atmos. Environ. 41 (33), 6964e6973. URL: http://linkinghub.elsevier.com/retrieve/pii/ S1352231007004086. Robert, C.P., Casella, G., 2004. Monte Carlo Statistical Methods. Springer. Robertson, L., Langner, J., Dec. 1998. Source function estimate by means of variational data assimilation applied to the ETEX-I tracer experiment. Atmos. Environ. 32 (24), 4219e4225. URL: http://linkinghub.elsevier.com/retrieve/pii/ S1352231098001769. Seibert, P., Frank, a, Jan. 2004. Source-receptor matrix calculation with a Lagrangian particle dispersion model in backward mode. Atmos. Chem. Phys. 4 (1), 51e63. URL: http://www.atmos-chem-phys.net/4/51/2004/. Senocak, I., Hengartner, N.W., Short, M.B., Daniel, W.B., Oct. 2008. Stochastic event reconstruction of atmospheric contaminant dispersion using Bayesian inference. Atmos. Environ. 42 (33), 7718e7727. URL: http://linkinghub.elsevier.com/ retrieve/pii/S1352231008005189. Simon, D., Simon, D.L., Feb. 2010. Constrained Kalman filtering via density function truncation for turbofan engine health estimation. Int. J. Syst. Sci. 41 (2), 159e171. Singh, S.K., Rani, R., Aug. 2014. A least-squares inversion technique for identification of a point release: application to Fusion Field Trials 2007. Atmos. Environ. 92, 104e117. URL: http://linkinghub.elsevier.com/retrieve/pii/S1352231014002842. Winiarek, V., Vira, J., Bocquet, M., Sofiev, M., Saunier, O., Jun. 2011. Towards the operational estimation of a radiological plume using data assimilation after a radiological accidental atmospheric release. Atmos. Environ. 45 (17), 2944e2955. URL: http://linkinghub.elsevier.com/retrieve/pii/ S1352231010010514. Yee, E., Mar. 2008. Theory for reconstruction of an unknown number of contaminant sources using probabilistic inference. Bound. Layer Meteorol. 127 (3), 359e394. URL: http://link.springer.com/10.1007/s10546-008-9270-5. Yee, E., 2009. Validation of a Sensor-driven Modeling Paradigm for Multiple Source Reconstruction with FFT-07 Data. Tech. Rep. May. Defence Research and Development Canada, Suffield. Yee, E., Hoffman, I., Ungar, K., 2014. Bayesian inference for source reconstruction: a real-world application. Int. Sch. Res. Notices 2014 (1), 1e12. URL: http://www. hindawi.com/journals/isrn/2014/507634/.