About this Journal Submit a Manuscript Table of Contents
ISRN Applied Mathematics
Volumeย 2012ย (2012), Article IDย 465320, 20 pages
http://dx.doi.org/10.5402/2012/465320
Research Article

Inverse Dispersion for an Unknown Number of Sources: Model Selection and Uncertainty Analysis

Defence R&D Canada โ€“ Suffield, P.O. Box 4000 Stn Main, Medicine Hat, AB, Canada T1A 8K6

Received 20 June 2012; Accepted 10 July 2012

Academic Editors: Y.ย Dimakopoulos, J. R.ย Fernandez, and T. Y.ย Kam

Copyright ยฉ 2012 Her Majesty the Queen in Right of Canada. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

A simple recursive method is presented for performing the inverse dispersion modeling of an unknown number of (localized) sources, given a finite number of noisy concentration data acquired by an array of detectors. Bayesian probability theory is used to address the problem of selecting the source model which is most plausible in view of the given concentration dataset and all the available prior information. The recursive algorithm involves subtracting a predicted concentration signal arising from a source model consisting of ๐‘ localized sources from the measured concentration data for increasing values of ๐‘ and examining the resulting residual data to determine if the residuals are consistent with the estimated noise level in the concentration data. The method is illustrated by application to a real concentration dataset obtained from an atmospheric dispersion experiment involving the simultaneous release of a tracer from four sources.

1. Introduction

Significant advances in sensing technology for concentration measurements of contaminants (e.g., toxic industrial materials; chemical, biological, and radiological agents) released into the atmosphere, either accidentally or deliberately, have fostered interest in exploiting this information for the reconstruction of the contaminant source distribution responsible for the observed concentration. Indeed, innovative methods for measuring and collecting concentration data in formal observation (or monitoring) networks of ever more compact (nanotechnology-enabled) sensors, including the adoption of new wireless sensor network technologies for acquisition of data in situ and for rapid delivery of this information through wireless transmission [1], have placed the emphasis on data assimilation and inverse dispersion modeling.

This should not be too surprising owing to the fact that inverse dispersion modeling will enable the full realization of the benefits in this new data context, allowing the innovative fusion of sensor data with predictive models for atmospheric dispersion for determination of the unknown source characteristics. In turn, this will lead to a greatly improved situational awareness and result in a significantly enhanced common operational picture required for making more informed decisions for mitigation of the consequences of a (toxic) contaminant release. In this context, inverse dispersion modeling for source reconstruction has important implications for public safety and security.

In the past, two different methodologies have been used to address the inverse dispersion modeling problem, namely, deterministic optimization and stochastic Bayesian approaches. In the optimization method, the parameters ๐œƒ describing the source model are obtained as the solution of a nonlinear least-squares optimization problem: ฬ‚๐œƒ=argmin๐œƒโˆˆโ„๐‘›๎€ฝ๐‘“๐‘‚(๐œƒ)+๐‘“๐‘…๎€พ,(๐œƒ)(1.1) where ๐‘› is the number of parameters (unknowns) used to parameterize the source model, ๐‘“๐‘‚(๐œƒ) is an objective functional that measures the total mismatch (usually taken as the sum-of-squared differences) between the measured concentration data and the synthetic (predicted) concentration data associated with the current source model ๐œƒ, and ๐‘“๐‘…(๐œƒ) is a regularization functional that is used either to impose an additional constraint on the solution or to incorporate a priori information (which is required to produce a mathematically unique solution). Most efforts in application of the optimization approach to inverse dispersion modeling have been focused on the development of numerical methods for solution of the optimization problem or on the prescription of the regularization functional for incorporation of various forms of a priori information. The optimization method for inverse dispersion modeling has been used by Robertson and Langner [2], Bocquet [3], Thomson et al. [4], Allen et al. [5], and Issartel et al. [6] (among others) based on various numerical methods for solution of the optimization problem (e.g., variational data assimilation, genetic algorithm, conjugate gradient algorithm with line search, augmented Lagrangian method) and a number of different forms of regularization (e.g., energy, entropy, smoothness, or flatness).

While the optimization approach for inverse dispersion modeling seeks to provide a single optimal solution for the problem, the stochastic Bayesian approach seeks to generate multiple solutions to the problem with the evaluation of the degree of plausibility for each solution. In contrast to the optimization approach, the stochastic Bayesian approach allows the quantification of the uncertainty in the source reconstruction (arising as such from the use of incomplete and noisy concentration data and imperfect models of atmospheric dispersion for the inverse dispersion modeling). This approach has been developed and used by a number of researchers (e.g., [7โ€“10]).

In all the previous studies cited above, the inverse dispersion modeling involved the problem of identification of source parameters for a single localized source. The determination of the characteristics of multiple localized sources was briefly addressed by Yee [11] and by Sharan et al. [12] using, respectively, a stochastic Bayesian approach and an optimization (least-squares) approach that were similar to those applied previously for the recovery of a single source. In both these cases, this was possible because it was assumed in these investigations that the number of sources was known a priori. The problem of reconstruction of an unknown number of localized sources using a finite number of noisy concentration measurements is a significantly more difficult problem. A solution for this problem was proposed by Yee [13, 14] who approached the problem as a generalized parameter estimation problem in which the number of localized sources, ๐‘, in the (unknown) source distribution was included explicitly in the parameter vector ๐œƒ, in addition to the usual parameters that characterize each localized source (e.g., location, emission rate, source-on time, source-off time).

Solving the reconstruction of an unknown number of localized sources as a generalized parameter estimation problem posed a number of conceptual and technical difficulties. The primary conceptual difficulty resided in the fact that when the number of sources is unknown a priori, the dimension of the parameter space (or, equivalently, the dimension or length of the associated parameter vector ๐œƒ) is an unknown that needs to be estimated using the noisy concentration data. Yee [13, 14] overcame this conceptual difficulty by using Bayesian probability theory to formulate the full joint probability density function (PDF) for the number of sources ๐‘ and the parameters of the ๐‘ localized sources and then demonstrated how a reversible-jump Markov chain Monte Carlo (RJMCMC) algorithm can be designed to draw samples of candidate source models from this joint PDF, allowing for changes in the dimensionality of the source model (or equivalently, changes in the number of localized sources in the source distribution). Furthermore, it was found that the RJMCMC algorithm sampled the parameter space of the unknown source distribution rather inefficiently. To overcome this technical difficulty, Yee [13, 14] showed how the RJMCMC algorithm can be combined either with a form of parallel tempering based on a Metropolis-coupled MCMC algorithm [13] or with a simpler and computationally more efficient simulated annealing scheme [14] to improve significantly the sampling efficiency (or, mixing rate) of the Markov chain in the variable dimension parameter space.

In this paper, we address the inverse dispersion modeling problem for the difficult case of multiple sources when the number of sources is unknown a priori as a model selection problem, rather than as a generalized parameter estimation problem as described by Yee [13, 14]. The model selection problem here is formulated in the Bayesian framework which involves the evaluation of model selection statistics that gauges rigorously the tradeoff between the goodness-of-fit of the source model structure to the concentration data and the complexity of the source model structure. The model selection approach proposed herein for reconstruction of an unknown number of localized sources is conceptually and algorithmically simpler than that based on a generalized parameter estimation approach [13, 14], leading as such to a simpler and more efficient computer implementation. More specifically, the model selection approach is simpler than the generalized parameter estimation approach in that it does not need to deal with the complexity of the variable dimensionality of the parameter space, which requires generally the development of sophisticated and computationally intensive algorithms.

2. Bayesian Model Selection

Let ๐‘š๐‘ denote a source distribution (model) which consists of a known number, ๐‘, of localized sources. The model ๐‘š๐‘ is characterized by a set of parameters ๐œƒ. Bayesian inference for ๐œƒ, when given the concentration data ๐’Ÿ, is based on Bayesโ€™ theorem (rule): ๐‘๎€ท๐œƒโˆฃ๐‘š๐‘๎€ธ=๐‘๎€ท,๐’Ÿ,๐ผ๐œƒโˆฃ๐‘š๐‘๎€ธ๐‘๎€ท,๐ผ๐’Ÿโˆฃ๐‘š๐‘๎€ธ,๐œƒ,๐ผ๐‘๎€ท๐’Ÿโˆฃ๐‘š๐‘๎€ธ,๐ผ,(2.1) where ๐ผ denotes the background (contextual) information that defines the multiple source reconstruction problem (e.g., background meteorology, atmospheric dispersion model used to determine the source-receptor relationship, etc.) and the vertical bar ๎…ขโˆฃโ€ฒโ€ฒ denotes โ€œconditional uponโ€™โ€™ or โ€œgiven.โ€™โ€™ Furthermore, in (2.1), ๐‘(๐œƒโˆฃ๐‘š๐‘,๐’Ÿ,๐ผ)โ‰ก๐‘(๐œƒ) is the posterior probability distribution of the parameters ๐œƒ, ๐‘(๐œƒโˆฃ๐‘š๐‘,๐ผ)โ‰ก๐œ‹(๐œƒ) is the prior probability distribution of the parameters ๐œƒ before the data ๐’Ÿ was made available, and ๐‘(๐’Ÿโˆฃ๐‘š๐‘,๐œƒ,๐ผ)โ‰กโ„’(๐œƒ) is the likelihood function and defines the probabilistic model of how the data were generated. Finally, ๐‘(๐’Ÿโˆฃ๐‘š๐‘,๐ผ) is the evidence and, in the context of parameter estimation, ensures proper normalization of the posterior probability distribution (assuming ๐œƒโˆˆโ„๐‘›): ๐‘๎€ท๐’Ÿโˆฃ๐‘š๐‘๎€ธ=๎€œ,๐ผโ„๐‘›๎€ท๐‘šโ„’(๐œƒ)๐œ‹(๐œƒ)๐‘‘๐œƒโ‰ก๐‘๐‘๎€ธ.(2.2)

However, in our problem, the number of sources (๐‘) is unknown a priori, so the relevant question that needs to be addressed is as follows: given a set of possible source models ๐‘†โ‰ก{๐‘š0,๐‘š1,โ€ฆ,๐‘š๐‘max} (๐‘š๐‘ denotes the source model consisting of ๐‘ localized sources), which source model is the most plausible (probable) given the concentration data ๐’Ÿ? This question is addressed rigorously in the Bayesian framework by using Bayesโ€™ theorem to compute the posterior probability for source model ๐‘š๐‘ (๐‘=0,1,โ€ฆ,๐‘max) to give ๐‘๎€ท๐‘š๐‘๎€ธ=๐‘๎€ท๐‘šโˆฃ๐’Ÿ,๐ผ๐‘๎€ธ๐‘๎€ทโˆฃ๐ผ๐’Ÿโˆฃ๐‘š๐‘๎€ธ,๐ผ.๐‘(๐’Ÿโˆฃ๐ผ)(2.3) Here, ๐‘(๐‘š๐‘โˆฃ๐ผ) is the prior probability of source model ๐‘š๐‘ given only the information ๐ผ, and ๐‘(๐’Ÿโˆฃ๐‘š๐‘๐ผ) is the global likelihood of the concentration data ๐’Ÿ given the source model ๐‘š๐‘ and ๐ผ and represents the goodness-of-fit of the model to the data. Note that the global likelihood of (2.3) is identical to the evidence ๐‘(๐‘š๐‘) defined in (2.2). Finally, the marginal probability of the data given the background information, ๐‘(๐’Ÿโˆฃ๐ผ), is a normalization constant over all source models: ๐‘(๐’Ÿโˆฃ๐ผ)=๐‘max๎“๐‘=0๐‘๎€ท๐‘š๐‘๎€ธ๐‘๎€ทโˆฃ๐ผ๐’Ÿโˆฃ๐‘š๐‘๎€ธ.,๐ผ(2.4)

Given any two source models ๐‘š๐‘0 and ๐‘š๐‘1 from the set ๐‘†, the โ€œoddsโ€™โ€™ in favor of model ๐‘š๐‘1 relative to model ๐‘š๐‘0 are given by the ratio of the posterior probabilities of the two models, which on using (2.3) and (2.2), reduces to ๐พ10โ‰ก๐‘๎€ท๐‘š๐‘1๎€ธโˆฃ๐’Ÿ,๐ผ๐‘๎€ท๐‘š๐‘0๎€ธ=๐‘๎€ทโˆฃ๐’Ÿ,๐ผ๐’Ÿโˆฃ๐‘š๐‘1๎€ธ,๐ผ๐‘๎€ท๐’Ÿโˆฃ๐‘š๐‘0๎€ธโ‹…๐‘๎€ท๐‘š,๐ผ๐‘1๎€ธโˆฃ๐ผ๐‘๎€ท๐‘š๐‘0๎€ธ=๐‘๎€ท๐‘šโˆฃ๐ผ๐‘1๎€ธ๐‘๎€ท๐‘š๐‘0๎€ธโ‹…๐‘๎€ท๐‘š๐‘1๎€ธโˆฃ๐ผ๐‘๎€ท๐‘š๐‘0๎€ธ.โˆฃ๐ผ(2.5) Note that the posterior odds ratio ๐พ10 is the product of the evidence ratio ๐‘(๐‘š๐‘1)/๐‘(๐‘š๐‘0) and the prior odds ratio ๐‘(๐‘š๐‘1โˆฃ๐ผ)/๐‘(๐‘š๐‘0โˆฃ๐ผ). If every source model in the set ๐‘† is equally probable (namely, ๐‘(๐‘š๐‘โˆฃ๐ผ)=1/(๐‘max+1) for ๐‘=0,1,โ€ฆ,๐‘max), then the prior odds ratio is unity, and the posterior odds ratio is identically equal to the evidence ratio. When viewed as a model selection problem, the inverse dispersion modeling of an unknown number of sources is conceptually simple: compute the evidence ๐‘(๐‘š๐‘) for each source model ๐‘š๐‘ in ๐‘† given the input concentration data ๐’Ÿ and the background information ๐ผ, select the model with the largest value of the evidence as the most probable model, and recover the parameters ๐œƒ corresponding to this most probable model.

The key to the model selection problem reduces to the computation of the evidence ๐‘(๐‘š๐‘) for an arbitrary value of ๐‘. A perusal of (2.2) shows that the evidence for the source model structure ๐‘š๐‘ (namely, a source model consisting of exactly ๐‘ localized sources) involves a computation of the overlap integral (or inner product) of the likelihood โ„’(๐œƒ) and the prior distribution ๐œ‹(๐œƒ) over the entire parameter space for model ๐‘š๐‘. This inner product can be interpreted as a metric (statistic) for model selection and, as such, is an intrinsic element of the source model structure ๐‘š๐‘ in the sense that it depends on both the set of parameters to be varied and the prior ranges for those parameters but is independent of the most probable (preferred) values for the parameters.

Note that the model selection metric selects preferentially a source model structure with the largest overlap of the likelihood and prior distribution, implying that in order for a more complex source model (with a necessarily higher-dimensional parameter space) to be selected, this requires the prior information to be smeared over a larger hypervolume in the parameter space in order to overlap significantly the likelihood. This mechanism embodies automatically an Occamโ€™s razor in that a simpler source model structure is preferentially selected unless the likelihood of the data for a more complex source model increases significantly more than what would be expected from simply fitting the noise (residual between the measured concentration and the predicted concentration provided by the model).

The solution of the model selection problem is conceptually simple involving, as such, the computation of ๐‘(๐‘š๐‘) (๐‘š๐‘โˆˆ๐‘†). Unfortunately, this computation is a technically difficult problem for two reasons. Firstly, as ๐‘ (or, equivalently, as the number of localized sources) increases, the overlap integral in (2.2) needs to be evaluated over a parameter space of increasing dimensionality ๐‘›. More specifically as ๐‘ increases, the dimensionality of the parameter space ๐‘› increases linearly with ๐‘. Secondly, the integrand of the overlap integral of (2.2) becomes ever more complex with increasing ๐‘, owing to a counting degeneracy in the problem. In particular, it is noted that โ„’(๐œƒ) is invariant under a reordering (relabelling) of the identifiers used to label the individual localized sources in the source model. This degeneracy simply corresponds to different (but physically equivalent) identifications of what is meant by a particular localized source in the source distribution, implying that the number of modes in โ„’(๐œƒ) increases with ๐‘š๐‘ as ๐‘!. The factor ๐‘! here corresponds simply to the number of possible permutations of the labels for the localized sources in the source model ๐‘š๐‘. In summary, the increase in the dimensionality of the parameter space and in the complexity of โ„’(๐œƒ) with increasing ๐‘ makes the problem of computation of ๐‘(๐‘š๐‘) technically difficult.

In view of the difficulty arising from the numerical evaluation of ๐‘(๐‘š๐‘) for increasing values of ๐‘, we consider an alternative (but closely related) methodology for addressing the model selection problem in the context of inverse dispersion modeling of an a priori unknown number of localized sources. To this purpose, the question that needs to be posed is after the removal of a predicted concentration for a source model ๐‘š๐‘ from the measured concentration data, are the residuals that remain consistent with the estimated noise level or is there further evidence for the existence of an additional localized source in the residuals? This question is easily addressed using Bayesโ€™ theorem as follows.

Towards this purpose, the following model is assumed for the measured concentration data ๐‘‘๐‘– (๐‘–=1,2,โ€ฆ,๐‘๐‘š): ๐‘‘๐‘–=๐ถ๐‘–+๐‘’๐‘–(1)=๎‚๐ถ๐‘–+๐‘’๐‘–(1)+๐‘’๐‘–(2)โ‰ก๎‚๐ถ๐‘–+๐‘’๐‘–,(2.6) where ๐ถ๐‘– is the true (albeit unknown) mean concentration, ๐‘’๐‘–(1) is the measurement error, ๎‚๐ถ๐‘– is the predicted mean concentration obtained from an atmospheric dispersion model, ๐‘’๐‘–(2) is the model error incurred by using the (necessarily) imperfect atmospheric model to predict ๐ถ๐‘–, and ๐‘’๐‘– is the composite error that includes both the measurement and model error. Furthermore, it is assumed that the expectation value โŸจ๐‘’๐‘–โŸฉ of the composite error is zero, and the variance โŸจ๐‘’2๐‘–โŸฉ is assumed (for now) to be known and given by ๐œŽ2๐‘–. Hence, it is assumed that we are given the concentration data ๐’Ÿโ‰ก{๐‘‘1,๐‘‘2,โ€ฆ,๐‘‘๐‘๐‘š} and corresponding uncertainties {๐œŽ1,๐œŽ2,โ€ฆ,๐œŽ๐‘๐‘š} where ๐œŽ๐‘– is the standard deviation (square root of the variance) for the ๐‘–-th concentration datum ๐‘‘๐‘–. Normally, exact values for the uncertainties ๐œŽ๐‘– are unknown so that all we would have available would be some estimated values for these uncertainties ๐‘ ๐‘–, rather than the true values ๐œŽ๐‘–.

If we assume a priori that there are ๐‘ localized sources in the domain, we can analyse the concentration data ๐’Ÿ given that ๐‘š๐‘ is the correct model structure and draw samples of source models with exactly ๐‘ localized sources (encoded by ๐œƒ) from the posterior distribution given by (2.1). Assume that ๐‘๐‘  samples of source models with ๐‘ localized sources have been drawn from the posterior distribution; namely, we have available the samples ๐œƒ๐‘˜(๐‘š๐‘)(๐‘˜=1,2,โ€ฆ,๐‘๐‘ ) where ๐œƒ๐‘˜(๐‘š๐‘) is ๐‘˜th source model sample consisting of exactly ๐‘ localized sources. To determine if ๐‘ is the correct number of localized sources, we can compute the residual data โ„ฐโ‰ก{๐œ–1,๐œ–2,โ€ฆ,๐œ–๐‘๐‘š} where ๐œ–๐‘– is the mean residual for the ๐‘–th concentration datum (with realizations of the residuals obtained by subtracting the predicted concentration ๎‚๐ถ๐‘–(๐‘š๐‘) for a sample of the source model ๐‘š๐‘ encoded in ๐œƒ๐‘˜(๐‘š๐‘) from the measured concentration ๐‘‘๐‘–). The mean residual datum ๐œ–๐‘– is estimated from the ensemble of source model samples ๐œƒ๐‘˜(๐‘š๐‘) (๐‘˜=1,2,โ€ฆ,๐‘๐‘ ) as follows: ๎๐œ–๐‘–=1๐‘๐‘ ๐‘๐‘ ๎“๐‘˜=1๎‚€๐‘‘๐‘–โˆ’๎‚๐ถ๐‘–๎€ท๐œƒ๐‘˜๎€ท๐‘š๐‘๎‚.๎€ธ๎€ธ(2.7) Similarly, the uncertainty ๐‘ ๐œ–๐‘– (error) in the residual datum ๐œ–๐‘– is estimated as ๐‘ ๐œ–๐‘–๎๐œ–=[2๐‘–โˆ’(๎๐œ–๐‘–)2]1/2, where ๎๐œ–2๐‘–=1๐‘๐‘ ๐‘๐‘ ๎“๐‘˜=1๎‚€๐‘‘๐‘–โˆ’๎‚๐ถ๐‘–๎€ท๐œƒ๐‘˜๎€ท๐‘š๐‘๎‚๎€ธ๎€ธ2.(2.8)

Now, given the residual data โ„ฐ, we evaluate the evidence for the presence of an additional source by considering two models (hypotheses) for โ„ฐ: (1)โ€‰โ€‰๐ป0: there is no additional source and the residual data ๐œ–๐‘–=๎‚๐ถ๐‘–(๐‘š0)+๐‘’๐‘–=๐‘’๐‘– (where ๎‚๐ถ๐‘–(๐‘š0) is the predicted concentration for a source model with no localized sources and hence has zero value), and (2)โ€‰โ€‰๐ป1: there is an additional source and the residual data ๐œ–๐‘–=๎‚๐ถ๐‘–(๐‘š1)+๐‘’๐‘– (where ๎‚๐ถ๐‘–(๐‘š1) is the predicted concentration for a source model consisting of one localized source). Under hypothesis ๐ป0 (null hypothesis), the residuals are consistent with the estimated noise level ๐‘’๐‘–. To decide which hypothesis is favoured by the data, we calculate the posterior odds ratio for ๐ป1 relative to ๐ป0. From (2.5) and assuming that hypotheses ๐ป0 and ๐ป1 are equally probable so ๐‘(๐ป0โˆฃ๐ผ)=๐‘(๐ป1โˆฃ๐ผ)=1/2, we get ๐พ๐œ–10โ‰ก๐‘๎€ท๐ป1๎€ธโˆฃโ„ฐ,๐ผ๐‘๎€ท๐ป0๎€ธโ‰ก๐‘๎€ท๐‘šโˆฃโ„ฐ,๐ผ1๎€ธโˆฃโ„ฐ,๐ผ๐‘๎€ท๐‘š0๎€ธ=๐‘๎€ทโˆฃโ„ฐ,๐ผโ„ฐโˆฃ๐‘š1๎€ธ,๐ผ๐‘๎€ทโ„ฐโˆฃ๐‘š0๎€ธโ‰ก๐‘,๐ผ๐œ–๎€ท๐‘š1๎€ธ๐‘๐œ–๎€ท๐‘š0๎€ธ.(2.9) If ๐พ๐œ–10 exceeds some preassigned threshold, then there is evidence for an additional source (namely, hypothesis ๐ป1 is favoured with respect to hypothesis ๐ป0). If this is the case, then the number of sources ๐‘ is increased by one, and we repeat the analysis of the concentration data ๐’Ÿ given that ๐‘š๐‘+1 is the correct model structure. This analysis would involve drawing samples from the posterior distribution of (2.1) with ๐‘š๐‘ replaced by ๐‘š๐‘+1, computing the residual data โ„ฐ (obtained now by subtracting the predicted concentration ๎‚๐ถ๐‘–(๐‘š๐‘+1) from the measured concentration ๐‘‘๐‘–), and then determining once more the evidence in favour of ๐ป1 relative to that of ๐ป0. When this test fails (namely, ๐พ๐œ–10 is less than or equal to the preassigned threshold value), the algorithm terminates. At this point, the number of sources ๐‘โˆ— has been determined, and the samples drawn from the posterior distribution of (2.1) for ๐‘š๐‘โˆ— can be used to estimate any posterior statistic of interest for the source parameters ๐œƒ.

To summarize, the algorithm for the inverse dispersion modeling of an unknown number of sources consists of the following steps.(1)Input concentration data ๐’Ÿ={๐‘‘1,๐‘‘2,โ€ฆ,๐‘‘๐‘๐‘š} and associated estimated uncertainties {๐‘ 1,๐‘ 2,โ€ฆ,๐‘ ๐‘๐‘š}.(2)Set ๐‘=1 and specify a threshold ๐พโˆ— for the posterior odds ratio. Draw ๐‘๐‘  samples of source models (encoded by the parameter vectors ๐œƒ๐‘˜(๐‘š๐‘), ๐‘˜=1,2,โ€ฆ,๐‘๐‘ ) from the posterior distribution ๐‘(๐œƒโˆฃ๐‘š๐‘,๐’Ÿ,๐ผ) given by (2.1).(3)Compute the evidences ๐‘(๐‘š๐‘) and ๐‘(๐‘š๐‘โˆ’1), and use these values to determine ๐พ10 in accordance to (2.5). If ๐พ10โ‰ค๐พโˆ—, stop (no sources are present in the domain, so ๐‘โˆ—=0).(4)Use samples ๐œƒ๐‘˜(๐‘š๐‘) (๐‘˜=1,2,โ€ฆ,๐‘๐‘ ) to estimate the residual data โ„ฐ={ฬ‚๐œ–1,ฬ‚๐œ–2,โ€ฆ,ฬ‚๐œ–๐‘๐‘ } and the associated uncertainty in this data {๐‘ ๐œ–1,๐‘ ๐œ–2,โ€ฆ,๐‘ ๐œ–๐‘๐‘š} in accordance to (2.7) and (2.8), respectively.(5)Using the information from the previous step, compute ๐พ๐œ–10 in accordance to (2.9). If ๐พ๐œ–10โ‰ค๐พโˆ—, set ๐‘โˆ—=๐‘, output ๐œƒ๐‘˜(๐‘š๐‘) (๐‘˜=1,2,โ€ฆ,๐‘๐‘ ), and stop.(6)Increase ๐‘โ†’๐‘+1. Draw ๐‘๐‘  samples ๐œƒ๐‘˜(๐‘š๐‘) (๐‘˜=1,2,โ€ฆ,๐‘๐‘ ) of source models from ๐‘(๐œƒโˆฃ๐‘š๐‘,๐’Ÿ,๐ผ) given by (2.1). Continue from Stepโ€‰โ€‰4.

The algorithm summarized previously determines the minimum number of localized sources needed to represent the concentration data ๐’Ÿ down to the estimated noise level, without the requirement to compute the evidence ๐‘(๐‘š๐‘) for ๐‘โ‰ฅ2 (a technically difficult problem computationally). Steps 2 and 3 in the algorithm correspond to the signal detection phase which addresses the specific question: is there evidence in the concentration data ๐’Ÿ for the existence of one or more localized sources in the domain? On the other hand, Steps 4 to 6 in the algorithm address the estimation phase and answer the question given that there is evidence for one or more localized sources in the domain (detection has been confirmed), how many localized sources are present, and what are the values of the source parameters that best characterize each of these localized sources?

3. Source-Receptor Relationship (Dispersion Modeling)

In this paper, we are interested in the reconstruction of an unknown number ๐‘ of localized sources. In view of this, the model for the source density function has the following explicit form (for source model ๐‘š๐‘): ๐’ฎ(๐ฑ,๐‘ก)=๐‘๎“๐‘—=1๐‘„๐‘—๐›ฟ๎€ท๐ฑโˆ’๐ฑ๐‘ ,๐‘—๐‘ˆ๎€ท๎€ธ๎€บ๐‘กโˆ’๐‘‡๐‘,๐‘—๎€ธ๎€ทโˆ’๐‘ˆ๐‘กโˆ’๐‘‡๐‘’,๐‘—,๎€ธ๎€ป(3.1) where ๐›ฟ(โ‹…) and ๐‘ˆ(โ‹…) are the Dirac delta and Heaviside unit step functions, respectively, and, ๐‘„๐‘—, ๐ฑ๐‘ ,๐‘—, ๐‘‡๐‘,๐‘—, and ๐‘‡๐‘’,๐‘— are the emission rate, vector position (location), source-on time, and source-off time, respectively, of the ๐‘—th source (๐‘—=1,2,โ€ฆ,๐‘). For a source model (or molecule) composed of ๐‘ localized sources (or source atoms), it is convenient to define the parameter vector as follows: ๐œƒโ‰ก(๐ฑ๐‘ ,1,๐‘‡๐‘,1,๐‘‡๐‘’,1,๐‘„1,โ€ฆ,๐ฑ๐‘ ,๐‘,๐‘‡๐‘,๐‘,๐‘‡๐‘’,๐‘,๐‘„๐‘)โˆˆโ„6๐‘.

To apply the Bayesian probability theory, it is necessary to relate the hypotheses of interest about the unknown source model ๐’ฎ to the modeled (predicted) concentration ๎‚๐ถ [cf. (2.6)]. The predicted concentration ๎‚๐ถ can be compared directly to the measured concentration ๐‘‘ โ€œseenโ€™โ€™ by a sensor at receptor location ๐ฑ๐‘Ÿ and time ๐‘ก๐‘Ÿ, averaged over the sensor volume and sampling time. The comparison can be effected by averaging the mean concentration ๐ถ(๐ฑ,๐‘ก) over this sensor volume and sampling time to give ๎‚๐ถ๎€ท๐ฑ๐‘Ÿ,๐‘ก๐‘Ÿ๎€ธโ‰ก๎€œ๐‘‡0๎€œ๐‘‘๐‘กโ„›๎€ท๐‘‘๐ฑ๐ถ(๐ฑ,๐‘ก)โ„Ž๐ฑ,๐‘กโˆฃ๐ฑ๐‘Ÿ,๐‘ก๐‘Ÿ๎€ธ๎€ท๐ฑโ‰กโŸจ๐ถโˆฃโ„ŽโŸฉ๐‘Ÿ,๐‘ก๐‘Ÿ๎€ธ,(3.2) where โ„Ž(๐ฑ,๐‘กโˆฃ๐ฑ๐‘Ÿ,๐‘ก๐‘Ÿ) is the spatial-temporal filtering function for the sensor concentration measurement at (๐ฑ๐‘Ÿ,๐‘ก๐‘Ÿ) and โ„›ร—[0,๐‘‡] corresponds to the space-time volume that contains the source distribution and the receptors (sensors). Note that the mean concentration ๎‚๐ถ(๐ฑ๐‘Ÿ,๐‘ก๐‘Ÿ) โ€œseenโ€™โ€™ by a sensor can be expressed succinctly as the inner or scalar product โŸจ๐ถโˆฃโ„ŽโŸฉ of the mean concentration ๐ถ and the sensor response function โ„Ž.

A source-oriented forward-time Lagrangian stochastic (LS) model can be used to predict ๐ถ(๐ฑ,๐‘ก). In this approach, ๐ถ is estimated from the statistical characteristics of โ€œmarkedโ€™โ€™ particle trajectories described by the following stochastic differential equation [15]: ๎€ท๐ถ๐‘‘๐—(๐‘ก)=๐•(๐‘ก)๐‘‘๐‘ก,๐‘‘๐•(๐‘ก)=๐š(๐—(๐‘ก),๐•(๐‘ก),๐‘ก)๐‘‘๐‘ก+0๎€ธ๐œ–(๐—(๐‘ก),๐‘ก)1/2๐‘‘๐–(๐‘ก),(3.3) where ๐—(๐‘ก)โ‰ก(๐‘‹๐‘–(๐‘ก))=(๐‘‹1(๐‘ก),๐‘‹2(๐‘ก),๐‘‹3(๐‘ก)) is the Lagrangian position and ๐•(๐‘ก)โ‰ก(๐‘‰๐‘–(๐‘ก))=(๐‘‰1(๐‘ก),๐‘‰2(๐‘ก),๐‘‰3(๐‘ก)) is the Lagrangian velocity of a โ€œmarkedโ€™โ€™ fluid element at time ๐‘ก (marked by the source as the fluid element passes through it at some earlier time ๐‘ก๎…ž). The state of the fluid particle at any time ๐‘ก after its initial release from the source density function ๐’ฎ is defined by (๐—,๐•). In (3.3), ๐ถ0 is the Kolmogorov โ€œuniversalโ€™โ€™ constant that is associated with the Kolmogorov similarity hypothesis for the form of the second-order Lagrangian velocity structure function in the inertial subrange; ๐œ– is the mean dissipation rate of turbulence kinetic energy; ๐‘‘๐–(๐‘ก)โ‰ก(๐‘‘๐‘Š๐‘–(๐‘ก))=(๐‘‘๐‘Š1(๐‘ก),๐‘‘๐‘Š2(๐‘ก),๐‘‘๐‘Š3(๐‘ก)) are the increments of a vector-valued (three-dimensional) Wiener process; ๐šโ‰ก(๐‘Ž๐‘–)=(๐‘Ž1,๐‘Ž2,๐‘Ž3) is the drift coefficient vector (or, more precisely, the conditional mean acceleration vector).

Unfortunately, for the current application, the computational demands of the source-oriented approach make it impractical as a direct tool for sampling from the posterior distribution because this involves necessarily exploring a large number of source parameter hypotheses. This is highly computer intensive, as the simulation-based Bayesian inference procedure requires a large number of forward calculations of the mean concentration to be performed, each of which involves the numerical solution of (3.2) and (3.3). In this particular case, it may be useful to construct an emulator for the simulation model and use this emulator as an inexpensive surrogate to replace the forward model [16โ€“18]. Fortunately, applying this type of approximation method for the forward model is not required for the current problem. It was shown by Keats et al. [7] and Yee et al. [10] that an exact computationally efficient procedure (appropriate for use in a Bayesian inference scheme) exists in the form of a receptor-oriented strategy for the representation of the source-receptor relationship.

Towards this purpose, the modeled sensor concentration ๎‚๐ถ(๐ฑ๐‘Ÿ,๐‘ก๐‘Ÿ) can be computed alternatively using the following dual representation: ๎‚๐ถ๎€ท๐ฑ๐‘Ÿ,๐‘ก๐‘Ÿ๎€ธโ‰ก๎€œ๐‘‡0๎€œ๐‘‘๐‘กโ„›๎€ท๐‘‘๐ฑ๐ถ(๐ฑ,๐‘ก)โ„Ž๐ฑ,๐‘กโˆฃ๐ฑ๐‘Ÿ,๐‘ก๐‘Ÿ๎€ธ=๎€œ๐‘ก๐‘Ÿ0๐‘‘๐‘ก0๎€œโ„›๐‘‘๐ฑ0๐ถโˆ—๎€ท๐ฑ0,๐‘ก0โˆฃ๐ฑ๐‘Ÿ,๐‘ก๐‘Ÿ๎€ธ๐’ฎ๎€ท๐ฑ0,๐‘ก0๎€ธโ‰กโŸจ๐ถโˆ—๎€ท๐ฑโˆฃ๐’ฎโŸฉ๐‘Ÿ,๐‘ก๐‘Ÿ๎€ธ,(3.4) where ๐ถโˆ—(๐ฑ0,๐‘ก0โˆฃ๐ฑ๐‘Ÿ,๐‘ก๐‘Ÿ) is the adjunct (or dual) concentration at space-time point (๐ฑ0,๐‘ก0) associated with the concentration datum at location ๐ฑ๐‘Ÿ at time ๐‘ก๐‘Ÿ (with ๐‘ก0โ‰ค๐‘ก๐‘Ÿ). A comparison of (3.2) and (3.4) implies the duality relationship โŸจ๐ถโˆฃโ„ŽโŸฉ=โŸจ๐ถโˆ—โˆฃ๐’ฎโŸฉ between ๐ถ and ๐ถโˆ— through the source functions โ„Ž and ๐’ฎ. Moreover, ๐ถโˆ— is uniquely defined in the sense that it is explicitly constructed so that it verifies this duality relationship.

The adjunct field ๐ถโˆ— in the dual representation can be estimated from the statistical characteristics of โ€œmarkedโ€™โ€™ particle trajectories corresponding to a receptor-oriented backward-time LS trajectory model, defined as the solution to the following stochastic differential equation (with ๐‘‘๐‘ก๎…ž>0): ๐‘‘๐—๐‘๎€ท๐‘ก๎…ž๎€ธ=๐•๐‘๎€ท๐‘ก๎…ž๎€ธ๐‘‘๐‘ก๎…ž,(3.5)๐‘‘๐•๐‘๎€ท๐‘ก๎…ž๎€ธ=๐š๐‘๎€ท๐—๐‘๎€ท๐‘ก๎…ž๎€ธ,๐•๐‘๎€ท๐‘ก๎…ž๎€ธ,๐‘ก๎…ž๎€ธ๐‘‘๐‘ก๎…ž+๎€ท๐ถ0๐œ–๎€ท๐—๐‘๎€ท๐‘ก๎…ž๎€ธ,๐‘ก๎…ž๎€ธ๎€ธ1/2๎€ท๐‘ก๐‘‘๐–๎…ž๎€ธ,(3.6) where at any given time ๐‘ก๎…ž, (๐—๐‘,๐•๐‘) is a point in the phase space along the backward trajectory of the โ€œmarkedโ€™โ€™ fluid element (here assumed to be marked or tagged as a fluid particle which at time ๐‘ก๐‘Ÿ passed through the spatial volume of the detector at location ๐ฑ๐‘Ÿ). Here, (๐—๐‘,๐•๐‘) defines the state of a fluid particle at any time ๐‘ก๎…ž<๐‘ก๐‘Ÿ before its โ€œfinalโ€ release from the sensor space-time volume โ„Ž at time ๐‘ก๐‘Ÿ. It can be shown [15, 19] that ๐ถโˆ— obtained from (3.5) for a detector with the filtering function โ„Ž and ๐ถ obtained from (3.2) for a release from the source density ๐’ฎ, are exactly consistent with the duality relationship โŸจ๐ถโˆฃโ„ŽโŸฉ=โŸจ๐ถโˆ—โˆฃ๐’ฎโŸฉ provided: (1)โ€‰โ€‰๐•๐‘ in (3.5) is related to ๐• in (3.3) as ๐•๐‘(๐‘ก)=๐•(๐‘ก), and (2)โ€‰โ€‰๐š๐‘ in (3.5) is related to ๐š in (3.3) as ๐‘Ž๐‘๐‘–(๐ฑ,๐ฎ,๐‘ก)=๐‘Ž๐‘–(๐ฑ,๐ฎ,๐‘ก)โˆ’๐ถ0๐œ•๐œ–(๐ฑ,๐‘ก)๐œ•๐‘ข๐‘–ln๐‘ƒ๐ธ(๐ฎ),(3.7) where ๐‘ƒ๐ธ(๐ฎ) is the background Eulerian velocity PDF. Thomson [15] provides one particular solution for the drift coefficient vector ๐š that is consistent with the well-mixed criterion for the case where ๐‘ƒ๐ธ(๐ฎ) has a Gaussian functional form (namely, for Gaussian turbulence).

Finally, if we substitute (3.1) into (3.4), the predicted concentration ๎‚๎‚๐ถ(๐œƒ)โ‰ก๐ถ(๐ฑ๐‘Ÿ,๐‘ก๐‘Ÿ;๐œƒ) โ€œseenโ€™โ€™ by the sensor at (๐ฑ๐‘Ÿ,๐‘ก๐‘Ÿ) is given explicitly by the following expression: ๎‚๐ถ(๐œƒ)=๐‘๎“๐‘—=1๐‘„๐‘—๎€œmin(๐‘‡,๐‘‡๐‘’,๐‘—)๐‘‡๐‘,๐‘—๐ถโˆ—๎€ท๐ฑ๐‘ ,๐‘—,๐‘ก๐‘ โˆฃ๐ฑ๐‘Ÿ,๐‘ก๐‘Ÿ๎€ธ๐‘‘๐‘ก๐‘ .(3.8)

4. Prior and Likelihood

4.1. Prior

Assuming that each parameter in ๐œƒ (for source model ๐‘š๐‘) is logically independent, the prior ๐œ‹(๐œƒ) can be factored as follows: ๐œ‹(๐œƒ)=๐‘๎‘๐‘—=1๐œ‹๎€ท๐‘„๐‘—๎€ธ๐œ‹๎€ท๐ฑ๐‘ ,๐‘—๎€ธ๐œ‹๎€ท๐‘‡๐‘,๐‘—๎€ธ๐œ‹๎€ท๐‘‡๐‘’,๐‘—๎€ธ.(4.1) In this paper, the component prior distributions are assigned uniform distributions over an appropriate range for each source parameter. Furthermore, it is noted that the parameter ranges specified in the prior define the region in โ„๐‘› (๐‘›โ‰ก6๐‘) over which the evidence integral is carried out [cf. (2.2)]. More specifically, for ๐‘—=1,2,โ€ฆ,๐‘, ๐‘„๐‘—โˆผ๐’ฐ([0,๐‘„max]), where ๐‘„max is an a priori upper bound on the emission rate; ๐ฑ๐‘ ,๐‘—โˆผ๐’ฐ(โ„›) where โ„›โŠ‚โ„3 is the a priori specified spatial domain that is assumed to contain the source ๐’ฎ; and, ๐‘‡๐‘,๐‘—โˆผ๐’ฐ([0,๐‘‡max]) and ๐‘‡๐‘’,๐‘—โˆผ๐’ฐ([๐‘‡๐‘,๐‘—,๐‘‡max]), where ๐‘‡max is the upper bound on the time at which the source was turned on or turned off. Here, the symbol โˆผ denotes โ€œis distributed asโ€™โ€™ and ๐’ฐ([๐‘Ž,๐‘]) is the uniform distribution defined on the closed interval [๐‘Ž,๐‘]. Finally, note in the specification of the prior for ๐‘‡๐‘’,๐‘— that the lower bound of the parameter range is ๐‘‡๐‘,๐‘—, encoding the fact that the time the source is turned off must necessarily occur after it has been turned on.

4.2. Likelihood

In the absence of a detailed knowledge of the noise distribution ๐‘’๐‘– in (2.6), other than the fact that it has a finite variance ๐œŽ2๐‘–, the application of the principle of maximum entropy [20] informs us that a Gaussian distribution is the most conservative choice for the direct probability (likelihood) of the concentration data ๐’Ÿ. In consequence, the likelihood function โ„’(๐œƒ) has the following form: 1โ„’(๐œƒ)โ‰กโ„’(๐œƒโˆฃ๐’Ÿ,๐œ)=โˆ๐‘๐‘š๐‘–=1โˆš2๐œ‹๐œŽ๐‘–๎‚€โˆ’1exp2๐œ’2๎‚,(๐œƒ)(4.2) where ๐œ’2(๐œƒ)โ‰ก๐‘๐‘š๎“๐‘–=1๎ƒฉ๐‘‘๐‘–โˆ’๎‚๐ถ๐‘–(๐œƒ)๐œŽ๐‘–๎ƒช2.(4.3) In (4.2), the notation for the likelihood function was adjusted to include the standard deviation for the noise ๐œโ‰ก(๐œŽ1,๐œŽ2,โ€ฆ,๐œŽ๐‘๐‘š) in the conditioning to emphasize the fact that {๐œŽ๐‘–}๐‘๐‘š๐‘–=1 are assumed to be known quantities.

Unfortunately, as mentioned previously, the noise term ๐‘’๐‘– in (2.6) is rather complicated, arising as such from a superposition of both model and measurement errors. In consequence, reliable estimates for ๐œŽ๐‘– (๐‘–=1,2,โ€ฆ,๐‘๐‘š) are difficult to obtain in practical applications. In view of this, let us denote by ๐‘ ๐‘– the quoted estimate of the standard deviation for the noise term ๐‘’๐‘– for which the true (but unknown) standard deviation is ๐œŽ๐‘–. Now, let us characterize the uncertainty in the specification of the standard deviation of ๐‘’๐‘– with an inverse gamma distribution of the following form (or, equivalently, prior distribution for ๐œŽ๐‘–): ๐œ‹๎€ท๐œŽ๐‘–โˆฃ๐‘ ๐‘–๎€ธ๐›ผ,๐›ผ,๐›ฝ=2๐›ฝ๎‚ต๐‘ ฮ“(๐›ฝ)๐‘–๐œŽ๐‘–๎‚ถ2๐›ฝ๎ƒฉ๐‘ expโˆ’๐›ผ2๐‘–๐œŽ2๐‘–๎ƒช1๐œŽ๐‘–,๐‘–=1,2,โ€ฆ,๐‘๐‘š,(4.4) where ฮ“(๐ฑ) denotes the gamma function and ๐›ผ and ๐›ฝ are scale and shape parameters, respectively, that define the inverse gamma distribution. The inverse gamma distribution for ๐œŽ๐‘– has mean โŸจ๐œŽ๐‘–โŸฉ=๐‘ ๐‘–๐›ผ1/2ฮ“(๐›ฝโˆ’1/2)/ฮ“(๐›ฝ) and variance Var[๐œŽ๐‘–]=๐‘ 2๐‘–๐›ผ[1/(๐›ฝโˆ’1)โˆ’ฮ“2(๐›ฝโˆ’1/2)/ฮ“2(๐›ฝ)]. Again, the parameters ๐›ผ and ๐›ฝ have been added to the PDF of the noise uncertainty in (4.4) to indicate that the values for these parameters are assumed to be known.

In view of (4.4), the true but unknown standard deviations ๐œŽ๐‘– of ๐‘’๐‘– that appear in (4.2) can be treated as nuisance parameters and eliminated by considering the integrated likelihood ๎€œ=๎€œโ„’(๐œƒโˆฃ๐’Ÿ,๐ฌ,๐›ผ,๐›ฝ)=โ„’(๐œƒโˆฃ๐’Ÿ,๐œ)๐œ‹(๐œโˆฃ๐ฌ,๐›ผ,๐›ฝ)๐‘‘๐œโ„’(๐œƒโˆฃ๐’Ÿ,๐œ)๐‘๐‘š๎‘๐‘–=12๐›ผ๐›ฝ๎‚ต๐‘ ฮ“(๐›ฝ)๐‘–๐œŽ๐‘–๎‚ถ2๐›ฝ๎ƒฉ๐‘ expโˆ’๐›ผ2๐‘–๐œŽ2๐‘–๎ƒช1๐œŽ๐‘–๐‘‘๐œ,(4.5) where ๐ฌโ‰ก(๐‘ 1,๐‘ 2,โ€ฆ,๐‘ ๐‘๐‘š) are the estimated standard deviations for the noise (๐‘’1,๐‘’2,โ€ฆ,๐‘’๐‘๐‘š) and ๐‘‘๐œโ‰ก๐‘‘๐œŽ1๐‘‘๐œŽ2โ‹ฏ๐‘‘๐œŽ๐‘๐‘š. Now, substituting the form for โ„’(๐œƒโˆฃ๐’Ÿ,๐œ) from (4.2) and (4.3) into (4.5) and performing the integration with respect to ๐œ, we obtain the integrated likelihood function given by โ„’(๐œƒโˆฃ๐’Ÿ,๐ฌ,๐›ผ,๐›ฝ)=๐‘๐‘š๎‘๐‘–=1๐›ผ๐›ฝฮ“(๐›ฝ+1/2)โˆš2๐œ‹๐‘ ๐‘–1ฮ“(๐›ฝ)๎‚ธ๎‚€๐‘‘๐›ผ+๐‘–โˆ’๎‚๐ถ๐‘–๎‚(๐œƒ)2/(2๐‘ 2๐‘–)๎‚น๐›ฝ+1/2.(4.6) The integrated likelihood of (4.6) can be interpreted simply as an average over all conditional likelihoods given ๐œ, weighted by their prior probabilities. In so doing, the integrated likelihood incorporates the uncertainties regarding the standard deviations for ๐‘’๐‘– (๐‘–=1,2,โ€ฆ,๐‘๐‘š).

The integrated likelihood function given by (4.6) depends explicitly on the hyperparameters ๐›ผ and ๐›ฝ for which values need to be assigned. In this paper, the values of ๐›ผ and ๐›ฝ are assigned as ๐›ผ=๐œ‹โˆ’1 and ๐›ฝ=1. The assignment ๐›ผ=๐œ‹โˆ’1 ensures that โŸจ๐œŽ๐‘–โŸฉ=๐‘ ๐‘– encoding our belief that our estimates ๐‘ ๐‘– of the standard deviation of ๐‘’๐‘– are unbiased. Furthermore, the assignment ๐›ฝ=1 results in a very heavy-tailed distribution for ๐œ‹(๐œŽ๐‘–โˆฃ๐‘ ๐‘–,๐›ผ,๐›ฝ) which allows significant deviations of the noise uncertainty from the quoted value of ๐‘ ๐‘– (provided by the user). Indeed, with the choice ๐›ฝ=1, the variance associated with ๐‘(๐œŽ๐‘–โˆฃ๐‘ ๐‘–,๐›ผ,๐›ฝ) in (4.4) becomes infinite. The heavy tail of the distribution is chosen to account for possibly significant under-estimations of the actual uncertainty (namely, the quoted uncertainty ๐‘ ๐‘–โ‰ช๐œŽ๐‘–). This could arise from inconsistencies in the model concentration predictions owing to structural model error or to โ€œoutliersโ€™โ€™ in the measured concentration data owing to either measurement error or perhaps distortion of the measured concentration data due to some unrecognized spurious (background) source.

In order to compute ๐‘๐œ–(๐‘š1) required in the evaluation of ๐พ๐œ–10 in (2.9), we are required to specify a functional form for the likelihood โ„’(๐œƒ(๐‘š1)โˆฃโ„ฐ,๐ฌ๐œ–) of the estimated mean residual data โ„ฐโ‰ก{ฬ‚๐œ–1,ฬ‚๐œ–2,โ€ฆ,ฬ‚๐œ–๐‘๐‘ } and the associated uncertainty in this data ๐ฌ๐œ–โ‰ก{๐‘ ๐œ–1,๐‘ ๐œ–2,โ€ฆ,๐‘ ๐œ–๐‘๐‘š}. To this end, the integrated likelihood function of (4.6) is also used for the specification of โ„’(๐œƒ(๐‘š1)โˆฃโ„ฐ,๐ฌ๐œ–): โ„’๎€ท๐œƒ๎€ท๐‘š1๎€ธโˆฃโ„ฐ,๐ฌ๐œ–๎€ธ๎€ท๐œƒ๎€ท๐‘šโ‰กโ„’1๎€ธโˆฃโ„ฐ,๐ฌ๐œ–๎€ธ=,๐›ผ,๐›ฝ๐‘๐‘š๎‘๐‘–=1๐›ผ๐›ฝฮ“(๐›ฝ+1/2)โˆš2๐œ‹๐‘ ๐œ–๐‘–ฮ“1(๐›ฝ)๎‚ธ๎‚€๐›ผ+ฬ‚๐œ–๐‘–โˆ’๎‚๐ถ๐‘–๎€ท๐œƒ๎€ท๐‘š1๎‚๎€ธ๎€ธ2/(2๐‘ ๐‘–๐œ–2)๎‚น๐›ฝ+1/2.(4.7) Here, ๐œƒ(๐‘š1) are the parameters of a source model consisting of a single localized source, used in the test to determine if the mean residual data โ„ฐ contains evidence for the presence of an additional source. Towards this objective, the likelihood function of (4.7) is used to calculate the evidence ๐‘๐œ–(๐‘š1) for this test as ๐‘๐œ–๎€ท๐‘š1๎€ธ=๎€œโ„๐‘›โ„’๎€ท๐œƒ๎€ท๐‘š1๎€ธโˆฃโ„ฐ,๐ฌ๐œ–๎€ธ๐œ‹๎€ท๐œƒ๎€ท๐‘š,๐›ผ,๐›ฝ1๎€ท๐‘š๎€ธ๎€ธ๐‘‘๐œƒ1๎€ธ,(4.8) where ๐‘›=6 (corresponding to the dimension of the parameter space for one source).

5. Computational Methods

There are two major computational problems that need to be addressed in order to fully specify the algorithm summarized in Section 2, namely, (1) specification of a method to draw samples from the posterior distribution ๐‘(๐œƒโˆฃ๐‘š๐‘,๐’Ÿ,๐ผ) (Steps 2 and 6) and (2) specification of a method for computation of the evidence ๐‘(๐‘š1) (Stepโ€‰โ€‰3, cf. (2.2)) or ๐‘๐œ–(๐‘š1) (Stepโ€‰โ€‰5, cf. (4.8)).

The most popular method that can be used for sampling from ๐‘(๐œƒโˆฃ๐‘š๐‘,๐’Ÿ,๐ผ) (see (2.1)) is a stochastic algorithm referred to as Markov chain Monte Carlo (MCMC). There has been considerable effort expended to improve the efficiency and the chain mixing efficacy of MCMC sampling procedures in order to allow rapid chain convergence to and efficient sampling from the posterior distribution of ๐œƒ. Improvements to MCMC sampling methods include: introduction of an adaptive Metropolis algorithm by Haario et al. [21], formulation of the differential evolution Monte Carlo algorithm by Ter Braak [22], development of a differential evolution adaptive Metropolis algorithm (DREAM) by Vrugt et al. [23], and design of multiple-try Metropolis algorithms by Liu et al. [24].

The are a number of options that are available for the computation of the evidences ๐‘(๐‘š1) and ๐‘๐œ–(๐‘š1). Owing to the fact that the determination of ๐‘(๐‘š1) and ๐‘๐œ–(๐‘š1) involves only an evaluation of a multidimensional integral in a low-dimensional space (โ„6 in this case), it is in principle possible to apply a brute force numerical integration [25] to address this problem. An alternative would be to apply a methodology that is applicable to the evaluation of the evidence in the general case, involving the evaluation of an overlap integral in a high-dimensional parameter space โ„๐‘› (๐‘›โ‰ซ1). These methodologies include importance sampling estimators [26] and thermodynamic integration [27].

Although there are various alternatives (leading to many different combinations) that could be used potentially to draw samples from ๐‘(๐œƒโˆฃ๐‘š๐‘,๐’Ÿ,๐ผ) and to evaluate ๐‘(๐‘š1) and ๐‘๐œ–(๐‘š1), we have used instead a single methodology to do both. The methodology that is used in this paper to address both these problems is nested sampling developed by Skilling [28] for the efficient evaluation of the evidence ๐‘ for the general case. The nested sampling algorithm transforms the multidimensional evidence integral of (2.2) into a simple one-dimensional representation: ๎€ท๐‘š๐‘โ‰ก๐‘๐‘๎€ธ=๎€œ10โ„’(๐œ’)๐‘‘๐œ’,(5.1) where ๐œ’๎€ทโ„’โˆ—๎€ธ=๎€œโ„’(๐œƒ)>โ„’โˆ—๐œ‹(๐œƒ)๐‘‘๐œƒ(5.2) is the prior mass in the parameter (hypothesis) space enclosed with likelihood greater than โ„’โˆ— and โ„’(๐œ’) is the inverse which labels the likelihood contour that encloses a prior mass ๐œ’. If we evaluate the likelihood โ„’(๐œ’) at a sequence of ๐‘š points ๐œ’๐‘– (๐‘–=1,2,โ€ฆ,๐‘š) with 0<๐œ’๐‘š<๐œ’๐‘šโˆ’1<โ‹ฏ<๐œ’2<๐œ’1<1 and โ„’(๐œ’๐‘–)>โ„’(๐œ’๐‘—) (๐‘–>๐‘—), the evidence ๐‘ can be approximated from the likelihood-ordered samples as the following weighted sum: ๎๐‘=๐‘š๎“๐‘–=1โ„’๎€ท๐œ’๐‘–๎€ธ๐›ฟ๐œ’๐‘–,๐›ฟ๐œ’๐‘–=๐œ’๐‘–โˆ’1โˆ’๐œ’๐‘–.(5.3)

If the prior mass points ๐œ’๐‘– are sampled in a logarithmic manner as ๐œ’๐‘–=โˆ๐‘–๐‘—=1๐‘ก๐‘— where ๐‘ก๐‘—โˆˆ(0,1) is a shrinkage ratio, then the nested algorithm consists of the following steps. The reader is referred to Skilling [28] for further details of the algorithm.(1)Set ๐‘–=0, ๐œ’0=1, ๐‘0=0 and ๐‘“=0.5 (preset fraction used in the stopping criterion). Randomly draw ๐‘€ samples ๐œƒ from the prior ๐œ‹(๐œƒ) to give an ensemble of samples. Evaluate the likelihood โ„’(๐œƒ) for each of the samples in the ensemble.(2)Increase ๐‘–โ†’๐‘–+1. Select the sample having the lowest likelihood (which we label as โ„’๐‘–) in the ensemble and remove (discard) it. Shrink the prior mass to ๐œ’๐‘–=๐œ’๐‘–โˆ’1๐‘’โˆ’1/๐‘€.(3)Draw a new sample ๐œƒ from the prior ๐œ‹(๐œƒ) subject to the hard likelihood constraint โ„’(๐œƒ)>โ„’๐‘–, and add this sample to the ensemble of samples.(4)Increment the evidence: ๐‘๐‘–=๐‘๐‘–โˆ’1+โ„’๐‘–(๐œ’๐‘–โˆ’1โˆ’๐œ’๐‘–).(5)If โ„’max๐œ’๐‘–<๐‘“๐‘๐‘– (โ„’max is the largest value of the likelihood in the ensemble of samples), add in contributions to the evidence ๐‘๐‘– from the ensemble of samples (remaining ๐‘€ samples that have not been discarded), and stop; otherwise, continue from Stepโ€‰โ€‰2.

In Stepโ€‰โ€‰5 of the algorithm, if the stopping criterion is satisfied, the estimate for the evidence ๎๐‘ is completed by adding the contribution of the remaining ๐‘€ samples in the ensemble to ๐‘๐‘–; namely, ๎๐‘โ†’๐‘๐‘–+๐‘€โˆ’1(โ„’(๐œƒ1)+โ„’(๐œƒ2)+โ‹ฏ+โ„’(๐œƒ๐‘€))๐œ’๐‘– where โ„’(๐œƒ๐‘˜) (๐‘˜=1,2,โ€ฆ,๐‘€) are the likelihood values evaluated at the remaining ๐‘€ samples. It is noted that the algorithm summarized above for estimation of ๐‘ (overlap integral of โ„’(๐œƒ) and ๐œ‹(๐œƒ)) automatically provides weighted samples ๐œƒ๐‘˜(๐‘š๐‘) drawn from the posterior distribution ๐‘(๐œƒโˆฃ๐‘š๐‘,๐’Ÿ,๐ผ)โˆโ„’(๐œƒ)๐œ‹(๐œƒ) (cf. (2.1) and (2.2)). More specifically, the ๐‘–th discarded sample in Stepโ€‰โ€‰2 of the algorithm can be interpreted as a sample drawn from the posterior distribution of ๐œƒ with weight given by ๐‘๐‘–=โ„’๐‘–(๐œ’๐‘–โˆ’1โˆ’๐œ’๐‘–๎๐‘)/ where ๎๐‘ is the estimate of the evidence obtained in Stepโ€‰โ€‰5 on termination of the algorithm.

The key part of the algorithm is Stepโ€‰โ€‰3 involving drawing a sample from the prior ๐œ‹(๐œƒ) within a prescribed hard likelihood constraint โ„’(๐œƒ)>โ„’โˆ—. To this purpose, Feroz et al. [29] developed a very efficient algorithm (which the authors refer to as MULTINEST) for sampling from a prior within a hard likelihood constraint involving a very sophisticated procedure for decomposition of the support of the likelihood above a given bound โ„’โˆ— into a set of overlapping ellipsoids and then sampling from the resulting ellipsoids. The algorithm is appropriate for sampling from posterior distributions with multiple modes and with pronounced curving degeneracies in a high-dimensional parameter space.

Finally, the posterior odds ratio ๐พ๐œ–10 in (2.9) is a summary of the evidence for ๐ป1 (additional source is present in the residuals) against ๐ป0 (no source is present in the residuals). To interpret how strong (or weak) is the evidence for ๐ป1 against ๐ป0 provided by the residual data โ„ฐ, we use a reference scale suggested by Jeffreys [30]. In this scale, log(๐พ๐œ–10)<1 corresponds to inconclusive evidence for ๐ป1 against ๐ป0 (or the evidence for ๐ป1 is not worth more than a bare mention and corresponds to a posterior odds ratio of less than about 3 to 1 in favor of ๐ป1). On this same scale, log(๐พ๐œ–10)>5 corresponds to very strong evidence for ๐ป1 against ๐ป0 (associated, as such, with a posterior odds ratio greater than about 150 to 1). In view of this scale for interpreting ๐พ๐œ–10, we choose the threshold ๐พโˆ— for the posterior odds ratio (required to be specified in Stepโ€‰โ€‰2 of the algorithm summarized in Section 2) to be log(๐พโˆ—)=1 (or, equivalently, ๐พโˆ—โ‰ˆ2.7).

6. Example: An Application of the Methodology

In this section, we apply the proposed methodology for inverse dispersion modeling of an unknown number of sources to a real dispersion data set obtained from a specific experiment conducted under the FUsing Sensor Information from Observing Networks (FUSION) Field Trial 2007 (FFT-07). The scientific objective of this field campaign was to acquire a comprehensive meteorological and dispersion dataset that can be used to validate methodologies developed for source reconstruction. Details of the instrumentation deployed and the experiments conducted in FFT-07 are given in Storwold [31], so only a brief summary of FFT-07 will be presented here. In particular, only the relevant details of the particular experiment that are required for the interpretation of the results in this paper are emphasized.

The experiments in FFT-07 were carried out in September 2007 at Tower Grid on US Army Dugway Proving Ground, Utah about 2โ€‰km west of Camel Back Ridge on the Great Salt Lake Desert. The easterly through southerly drainage flows that predominate during the early morning hours at this site originate on the higher terrain to the southeast and are channeled by Camel Back Ridge. Generally, the terrain was flat, uniform, and homogeneous consisting primarily of short grass interspersed with low shrubs with a height between about 0.25 and 0.75โ€‰m. The momentum roughness length, ๐‘ง0, was estimated to be ๐‘ง0=1.3ยฑ0.2cm.

In all the experiments, propylene (C3H6) was used as the tracer gas. The concentration detectors used were fast-response digital photoionization (dPID) detectors. In the experiments, a plume was formed in the atmospheric surface layer by releasing propylene from one or more (up to a maximum of four) purpose-designed gas dissemination systems. The network (or array) of concentration detectors consisted of up to 100โ€‰dPIDs, arranged in a staggered configuration of 10 rows of 10 detectors, with the rows of detectors separated by 50โ€‰m and the detectors along each row spaced 50โ€‰m apart. The concentration detectors along the ten sampling lines in the array were placed at a height, ๐‘ง๐‘‘, of 2.0โ€‰m.

In the experiment used to test the inverse dispersion modeling methodology proposed herein, tracer gas was released continuously over a period of 10โ€‰min from four source locations at a height, ๐‘ง๐‘ , of 2.0โ€‰m: (1) source 1 is at (๐‘ฅ๐‘ ,๐‘ฆ๐‘ )=(33.0,171.0)โ€‰m; (2) source 2 is at (๐‘ฅ๐‘ ,๐‘ฆ๐‘ )=(33.8,240.7)โ€‰m; (3) source 3 is at (๐‘ฅ๐‘ ,๐‘ฆ๐‘ )=(30.0,312.9)โ€‰m; (4) source 4 is at (๐‘ฅ๐‘ ,๐‘ฆ๐‘ )=(26.0,384.4)โ€‰m. The coordinates reported here for the source locations are referenced with respect to a local Cartesian coordinate system. Unfortunately, in this experiment, only the mass flow controller for source 3 functioned properly. The mass flow controllers for sources 1, 2, and 4 failed to properly regulate the flow, so the emission rates from these sources were unknown.

A three-dimensional sonic anemometer, placed at the 2โ€‰m level of a lattice tower located upwind of the array of concentration detectors, was used to characterize the background micrometeorological state of the atmospheric surface layer. For this experiment, the horizontal mean wind speed ๐‘†2 at the 2โ€‰m level, the friction velocity (๐‘ขโˆ—), the atmospheric stability (Obukhov length ๐ฟ), and the standard deviations of the velocity fluctuations in the alongwind (๐œŽ๐‘ข), crosswind (๐œŽ๐‘ฃ), and vertical (๐œŽ๐‘ค) directions were ๐‘†2=3.61โ€‰mโ€‰sโˆ’1, ๐‘ขโˆ—=0.282m sโˆ’1, ๐ฟ=โˆ’27.3m, ๐œŽ๐‘ข/๐‘ขโˆ—=2.33, ๐œŽ๐‘ฃ/๐‘ขโˆ—=1.86, and ๐œŽ๐‘ค/๐‘ขโˆ—=1.10. The mean wind direction in the experiment was normally incident on the detector array (namely, the mean wind direction corresponded to a wind from the +๐‘ฅ-direction and, hence, is perpendicular to the sampling lines of detectors along the ๐‘ฆ-axis).

The wind velocity and turbulence statistics were used in conjunction with Monin-Obukhov similarity theory relationships [32] as input to the atmospheric dispersion model (namely, the backward-time LS particle trajectory model described briefly in Section 3) for provision of the predicted concentration ๎‚๐ถ (cf. (3.4)). More specifically, the backward-time LS model of (3.5) was applied with the Kolmogorov constant ๐ถ0=4.8, a value which was recommended by Wilson et al. [33] from a calibration of the model against concentration data obtained from the Project Prairie Grass atmospheric dispersion experiments.

The example used here to illustrate the inverse dispersion methodology involve continuously emitting sources, so the relevant source parameters for source model ๐‘š๐‘ are ๐œƒโ‰ก๐œƒ(๐‘š๐‘)=(๐ฑ๐‘ ,1,๐‘„1,โ€ฆ,๐ฑ๐‘ ,๐‘,๐‘„๐‘). Furthermore, it is assumed that the height of the sources above ground level (๐‘ง๐‘ =2.0โ€‰m) is known a priori, so the only unknown location parameters are the (๐‘ฅ๐‘ ,๐‘ฆ๐‘ ) coordinates of the sources. Owing to the horizontal homogeneity of the mean flow and turbulence statistics in the current example, the adjunct concentration ๐ถโˆ—(๐ฑ๐‘ โˆฃ๐ฑ๐‘‘) in (3.4) can be precalculated for one detector location ๐ฑ๐‘‘ (for the known source height ๐‘ง๐‘ ), with the adjunct concentration ๐ถโˆ— (considered as a function of ๐ฑ๐‘ ) at all other detector locations obtained simply by a linear translation of ๐ถโˆ—(๐ฑ๐‘ โˆฃ๐ฑ๐‘‘) in the horizontal (๐‘ฅ,๐‘ฆ)-plane.

In order to calculate the likelihood function given by (4.6), we need to provide values for the estimated standard deviations ๐‘ ๐‘– for the noise ๐‘’๐‘–. The noise error variance ๐œŽ2๐‘– includes the sensor sampling error variance, ๐œŽ2๐‘‘,๐‘–, in the measurement of the concentration datum ๐‘‘๐‘– and the model error variance, ๐œŽ2๐‘š,๐‘–, in the prediction of ๎‚๐ถ๐‘–. The two contributions are combined in quadrature to give the noise error variance ๐œŽ2๐‘–=๐œŽ2๐‘‘,๐‘–+๐œŽ2๐‘š,๐‘–. The measurement error standard deviation is estimated as ๎๐œŽ๐‘‘,๐‘–โ‰ก๐‘ ๐‘‘,๐‘–=max(0.05,0.02๐‘‘๐‘–)โ€‰ppm (parts per million by volume) where the lower limit of 0.05โ€‰ppm represents the precision in the concentration measurements using the dPIDs. The model error standard deviation is estimated to be 20% of the predicted value of the concentration ๎‚๐ถ๐‘–(๐œƒ) (namely, ๎๐œŽ๐‘š,๐‘–โ‰ก๐‘ ๐‘š,๐‘–๎‚๐ถ=0.20๐‘–(๐œƒ)). In consequence, the estimated noise error standard deviation is given by ๎๐œŽ๐‘–โ‰ก๐‘ ๐‘–=(๐‘ 2๐‘‘,๐‘–+๐‘ 2๐‘š,๐‘–)1/2.

The algorithm described in Section 2 was used for the source reconstruction, applied to the concentration data ๐’Ÿโ‰ก{๐‘‘1,๐‘‘2,โ€ฆ,๐‘‘๐‘๐‘š} obtained from 62 detectors in the array (๐‘๐‘š=62) (see Figure 1 which shows the locations of the 62 detectors as filled squares). The hyperparameters defining the prior ๐œ‹(๐œƒ) have been chosen as follows: ๐‘„๐‘—โˆผ๐’ฐ([0,๐‘„max]) with ๐‘„max=100g sโˆ’1 and (๐‘ฅ๐‘ ,๐‘—,๐‘ฆ๐‘ ,๐‘—)โˆผ๐’ฐ(โ„›) with โ„›=[0,100]mร—[0,500]โ€‰m for ๐‘—=1,2,โ€ฆ,๐‘.

fig1
Figure 1: Four iterations of the algorithm of Section 2 showing the samples of source model ๐‘š๐‘ drawn from the posterior distribution ๐‘(๐œƒโˆฃ๐‘š๐‘,๐’Ÿ,๐ผ) and the computation of ๐พ๐œ–10 for (a) ๐‘=1, log(๐พ๐œ–10)=62.4ยฑ0.2; (b) ๐‘=2, log(๐พ๐œ–10)=57.3ยฑ0.3; (c) ๐‘=3, log(๐พ๐œ–10)=39.1ยฑ0.3; (d) ๐‘=4, log(๐พ๐œ–10)=0.92ยฑ0.33. A solid square denotes a detector.

After Steps 2 and 3 of the algorithm were completed, we found log(๐พ10)=log(๐‘(๐‘š1)/๐‘(๐‘š0))=64.0ยฑ0.2 implying that the concentration data strongly support the source model ๐‘š1 (one localized source or โ€œsignalโ€™โ€™) against the alternative source model ๐‘š0 (no localized source or โ€œno signalโ€™โ€™). In other words, at this point in the algorithm, a source for the concentration has been detected in the data ๐’Ÿ. Owing to the fact that log(๐พ10)>log(๐พโˆ—)=1, the algorithm continues to the iterative loop involving Steps 4 to 6. This iterative loop executes four times before it is terminated. The results of the execution of each of these four iterations are exhibited in Figure 1(a) (๐‘=1), Figure 1(b) (๐‘=2), Figure 1(c) (๐‘=3), and Figure 1(d) (๐‘=4). The caption in Figure 1 summarizes the values for log(๐พ๐œ–10)=log(๐‘๐œ–(๐‘š1)/๐‘๐œ–(๐‘š0)) obtained from testing the alternative hypothesis ๐ป1 (additional source present in the residual data) against the โ€œnullโ€™โ€™ hypothesis ๐ป0 (no additional source present in the residual data). Each panel in Figure 1 also displays a density plot of samples of the source model ๐‘š๐‘ drawn from the posterior distribution ๐‘(๐œƒโˆฃ๐‘š๐‘,๐’Ÿ,๐ผ) and projected onto the (๐‘ฅ,๐‘ฆ) horizontal plane (namely, each point in the plot corresponds to the (๐‘ฅ,๐‘ฆ) location of a sample of the source model ๐‘š๐‘ drawn from the posterior distribution). Note that the algorithm terminates with ๐‘โˆ—=๐‘=4 with log(๐พ๐œ–10)=0.92โ‰คlog(๐พโˆ—)=1, providing the correct number of localized sources in the source distribution (๐‘โˆ—=4 in this case).

The samples of source models ๐‘š๐‘ (๐‘=4) drawn from the posterior distribution ๐‘(๐œƒโˆฃ๐‘š๐‘,๐’Ÿ,๐ผ) in the last iteration of the algorithm (see Figure 1(d)) can be used to determine the characteristics (location, emission rate) of the four localized sources. Towards this objective, Figure 2 displays the univariate (diagonal) and bivariate (off-diagonal) marginal posterior distributions for the parameters of each source. For each univariate marginal posterior distribution of a parameter for a given source, the solid vertical line indicates the true value of the parameter (for this source), and the dashed vertical line corresponds to the best estimate of the parameter (for this source) obtained as the posterior mean. For the bivariate marginal posterior distribution of parameters of a given source, the solid square represents the position of the true source parameter values, and the solid circle indicates the best estimate of the true source parameter values obtained as the posterior means. The posterior mean, posterior standard deviation, and the lower and upper bounds for the 95% highest posterior distribution (HPD) interval of the source parameters for each identified source are summarized in Table 1. Comparing these estimated values of the source parameters, it can be seen that the proposed algorithm has adequately recovered the true parameters for each source (if these are known) to within the stated uncertainties.

tab1
Table 1: The posterior mean, posterior standard deviation, and lower and upper bounds of the 95% HPD interval of the parameters ๐‘ฅ๐‘ ,๐‘— (m), ๐‘ฆ๐‘ ,๐‘— (m), and ๐‘„๐‘—โ‰ก๐‘ž๐‘ ,๐‘— (gsโˆ’1) for ๐‘—=1, 2, 3, and 4 calculated from samples of source model ๐‘š๐‘ with ๐‘=4 drawn from the posterior distribution ๐‘(๐œƒ|๐‘š๐‘,๐’Ÿ,๐ผ) [cf. Figure 1(d)].
fig2
Figure 2: Univariate (diagonal) and bivariate (off-diagonal) marginal posterior distributions of the source parameters, namely, alongwind location ๐‘ฅ๐‘ , crosswind location ๐‘ฆ๐‘ , and emission rate ๐‘ž๐‘  that characterize sources 1 (a), 2 (b), 3 (c), and 4 (d) obtained from samples of source model ๐‘š๐‘ (๐‘=4) drawn from the posterior distribution ๐‘(๐œƒโˆฃ๐‘š๐‘,๐’Ÿ,๐ผ) ([cf. Figure 1(d)).

7. Conclusions

In this paper, we have proposed a Bayesian inference approach for addressing the inverse dispersion of an unknown number of sources using model comparison (selection). The necessary integrations to compute the evidence ๐‘(๐‘š๐‘) for ๐‘>1 can be very computationally demanding (as well as technically difficult), and, as a consequence, we developed an efficient and robust algorithm for model comparison that recursively removes the influence of a source model ๐‘š๐‘ from the measured concentration data ๐’Ÿ and tests the resulting residual data โ„ฐ to determine if the residual data are consistent with the estimated noise level. This test requires nothing more than the computation of the evidence ๐‘๐œ–(๐‘š๐‘) for ๐‘=1 which is usually computationally simple. The procedure finds the minimum number of localized sources necessary to represent the concentration signal in the data ๐’Ÿ down to the estimated noise level. Furthermore, the uncertainty in the estimated noise level (which includes contributions from both measurement and model errors) is treated by using an integrated likelihood, which was implemented using a specific prior (inverse gamma distribution) to represent the uncertainty in the estimated noise level. Nested sampling is used for the evidence computation, as well as for sampling from the posterior distribution ๐‘(๐œƒโˆฃ๐‘š๐‘,๐’Ÿ,๐ผ) in the proposed algorithm for inverse dispersion modeling.

The new algorithm has been applied successfully to a real concentration dataset obtained from an atmospheric dispersion experiment conducted in the FFT-07 field campaign consisting of the simultaneous release of a tracer from four sources. It is shown that the proposed algorithm performed well for this example: the number of sources was determined correctly (๐‘โˆ—=4), and for each of the identified sources the corresponding parameters (e.g., location, emission rate) were estimated reliably along with the determination of the uncertainty in the parameter (in the form of either a standard deviation or a credible interval that encloses a prespecified posterior probability mass). The methodology proposed herein offers a simpler alternative for the inverse dispersion modeling of an unknown number of sources that was addressed previously [13, 14] as a generalized parameter estimation problem, the latter of which involves necessarily the complexities in the design of an appropriate reversible-jump MCMC sampling algorithm to allow transdimensional jumps in the generalized parameter space (namely, jumps between the parameter space of source models ๐‘š๐‘ of different dimensions involving different numbers of sources ๐‘).

References

  1. L. Porta, T. H. Illangasekare, P. Loden, Q. Han, and A. P. Jayasumana, โ€œContinuous plume monitoring using wireless sensors: proof of concept in intermediate scale tank,โ€ Journal of Environmental Engineering, vol. 135, no. 9, pp. 831โ€“838, 2009. View at Publisher ยท View at Google Scholar ยท View at Scopus
  2. L. Robertson and J. Langner, โ€œSource function estimate by means of variational data assimilation applied to the ETEX-I tracer experiment,โ€ Atmospheric Environment, vol. 32, no. 24, pp. 4219โ€“4225, 1998. View at Scopus
  3. M. Bocquet, โ€œReconstruction of an atmospheric tracer source using the principle of maximum entropy. I: theory,โ€ Quarterly Journal of the Royal Meteorological Society, vol. 131, no. 610 B, pp. 2191โ€“2208, 2005. View at Publisher ยท View at Google Scholar ยท View at Scopus
  4. L. C. Thomson, B. Hirst, G. Gibson et al., โ€œAn improved algorithm for locating a gas source using inverse methods,โ€ Atmospheric Environment, vol. 41, no. 6, pp. 1128โ€“1134, 2007. View at Publisher ยท View at Google Scholar ยท View at Scopus
  5. C. T. Allen, G. S. Young, and S. E. Haupt, โ€œImproving pollutant source characterization by better estimating wind direction with a genetic algorithm,โ€ Atmospheric Environment, vol. 41, no. 11, pp. 2283โ€“2289, 2007. View at Publisher ยท View at Google Scholar ยท View at Scopus
  6. J. P. Issartel, M. Sharan, and S. K. Singh, โ€œIdentification of a point source by use of optimal weighted least squares,โ€ Pure and Applied Geophysics, vol. 169, pp. 467โ€“482, 2012.
  7. A. Keats, E. Yee, and F. S. Lien, โ€œBayesian inference for source determination with applications to a complex urban environment,โ€ Atmospheric Environment, vol. 41, no. 3, pp. 465โ€“479, 2007. View at Publisher ยท View at Google Scholar ยท View at Scopus
  8. A. Keats, E. Yee, and F. S. Lien, โ€œEfficiently characterizing the origin and decay rate of a nonconservative scalar using probability theory,โ€ Ecological Modelling, vol. 205, no. 3-4, pp. 437โ€“452, 2007. View at Publisher ยท View at Google Scholar ยท View at Scopus
  9. F. K. Chow, B. Kosović, and S. Chan, โ€œSource inversion for contaminant plume dispersion in urban environments using building-resolving simulations,โ€ Journal of Applied Meteorology and Climatology, vol. 47, no. 6, pp. 1533โ€“1572, 2008. View at Publisher ยท View at Google Scholar ยท View at Scopus
  10. E. Yee, F. S. Lien, A. Keats, and R. D'Amours, โ€œBayesian inversion of concentration data: source reconstruction in the adjoint representation of atmospheric diffusion,โ€ Journal of Wind Engineering and Industrial Aerodynamics, vol. 96, no. 10-11, pp. 1805โ€“1816, 2008. View at Publisher ยท View at Google Scholar ยท View at Scopus
  11. E. Yee, โ€œBayesian probabilistic approach for inverse source determination from limited and noisy chemical or biological sensor concentration measurements,โ€ in Chemical and Biological Sensing VIII, A. W. Fountain III, Ed., vol. 6554, 65540W, of Proceedings of the SPIE, p. 12, April 2007. View at Publisher ยท View at Google Scholar ยท View at Scopus
  12. M. Sharan, S. K. Singh, and J. P. Issartel, โ€œLeast square data assimilation for identification of the point source emissions,โ€ Pure and Applied Geophysics, vol. 169, pp. 483โ€“497, 2012.
  13. E. Yee, โ€œTheory for reconstruction of an unknown number of contaminant sources using probabilistic inference,โ€ Boundary-Layer Meteorology, vol. 127, no. 3, pp. 359โ€“394, 2008. View at Publisher ยท View at Google Scholar ยท View at Scopus
  14. E. Yee, โ€œProbability theory as logic: data assimilation for multiple source reconstruction,โ€ Pure and Applied Geophysics, vol. 169, pp. 499โ€“517, 2012.
  15. D. J. Thomson, โ€œCriteria for the selection of stochastic models for particle trajectories in turbulent flows,โ€ Journal of Fluid Mechanics, vol. 180, pp. 529โ€“556, 1987. View at Scopus
  16. M. C. Kennedy and A. O'Hagan, โ€œBayesian calibration of computer models,โ€ Journal of the Royal Statistical Society B, vol. 63, no. 3, pp. 425โ€“464, 2001. View at Publisher ยท View at Google Scholar
  17. J. Sacks, W. J. Welch, T. J. Mitchell, and H. P. Wynn, โ€œDesign and analysis of computer experiments,โ€ Statistical Science, vol. 4, no. 4, pp. 409โ€“435, 1989.
  18. D. Higdon, J. Gattiker, B. Williams, and M. Rightley, โ€œComputer model calibration using high-dimensional output,โ€ Journal of the American Statistical Association, vol. 103, no. 482, pp. 570โ€“583, 2008. View at Publisher ยท View at Google Scholar
  19. T. K. Flesch, J. D. Wilson, and E. Yee, โ€œBackward-time Lagrangian stochastic dispersion models and their application to estimate gaseous emissions,โ€ Journal of Applied Meteorology, vol. 34, no. 6, pp. 1320โ€“1332, 1995. View at Scopus
  20. E. T. Jaynes, Probability Theory: The Logic of Science, Cambridge University Press, Cambridge, Mass, USA, 2003. View at Publisher ยท View at Google Scholar
  21. H. Haario, E. Saksman, and J. Tamminen, โ€œAn adaptive Metropolis algorithm,โ€ Bernoulli, vol. 7, no. 2, pp. 223โ€“242, 2001. View at Publisher ยท View at Google Scholar
  22. C. J. F. Ter Braak, โ€œA Markov chain Monte Carlo version of the genetic algorithm differential evolution: easy Bayesian computing for real parameter spaces,โ€ Statistics and Computing, vol. 16, no. 3, pp. 239โ€“249, 2006. View at Publisher ยท View at Google Scholar
  23. J. A. Vrugt, C. J. F. Ter Braak, C. G. H. Diks, B. A. Robinson, J. M. Hyman, and D. Higdon, โ€œAccelerating Markov chain Monte Carlo simulation by differential evolution with self-adaptive randomized subspace sampling,โ€ International Journal of Nonlinear Sciences and Numerical Simulation, vol. 10, no. 3, pp. 273โ€“290, 2009. View at Scopus
  24. J. S. Liu, F. Liang, and W. H. Wong, โ€œThe multiple-try method and local optimization in Metropolis sampling,โ€ Journal of the American Statistical Association, vol. 95, no. 449, pp. 121โ€“134, 2000. View at Publisher ยท View at Google Scholar
  25. J. Berntsen, T. O. Espelid, and A. Genz, โ€œAn adaptive algorithm for the approximate calculation of multiple integrals,โ€ ACM Transactions on Mathematical Software, vol. 17, no. 4, pp. 437โ€“451, 1991. View at Publisher ยท View at Google Scholar
  26. M. A. Newton and A. E Rafferty, โ€œApproximate Bayesian inference by the weighted likelihood bootstrap,โ€ Journal of the Royal Statistical Society B, vol. 3, pp. 3โ€“48, 1994.
  27. A. Gelman and X.-L. Meng, โ€œSimulating normalizing constants: from importance sampling to bridge sampling to path sampling,โ€ Statistical Science, vol. 13, no. 2, pp. 163โ€“185, 1998. View at Publisher ยท View at Google Scholar
  28. J. Skilling, โ€œNested sampling for general Bayesian computation,โ€ Bayesian Analysis, vol. 1, no. 4, pp. 833โ€“859, 2006.
  29. F. Feroz, M. P. Hobson, and M. Bridges, โ€œMultiNest: an efficient and robust Bayesian inference tool for cosmology and particle physics,โ€ Monthly Notices of the Royal Astronomical Society, vol. 398, no. 4, pp. 1601โ€“1614, 2009. View at Publisher ยท View at Google Scholar ยท View at Scopus
  30. H. Jeffreys, Theory of Probability, Oxford University Press, Oxford, UK, 3rd edition, 1961.
  31. D. P. Storwold, Detailed Test Plan for Fusing Sensor Information from Observing Networks (FUSION) Field Trial 2007 (FFT-07), WDTC-TP-07-078, West Desert Test Center, US Army Dugway Proving Ground, 2007.
  32. R. B. Stull, An Introduction to Boundary-Layer Meteorology, Kluwer Academic Publishers, 1988.
  33. J. D. Wilson, T. K. Flesch, and L. A. Harper, โ€œMicro-meteorological methods for estimating surface exchange with a disturbed windflow,โ€ Agricultural and Forest Meteorology, vol. 107, no. 3, pp. 207โ€“225, 2001. View at Publisher ยท View at Google Scholar ยท View at Scopus