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: 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): 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 ):
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 ( 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 () to give 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:
Given any two source models and from the set , the “odds’’ in favor of model relative to model are given by the ratio of the posterior probabilities of the two models, which on using (2.3) and (2.2), reduces to Note that the posterior odds ratio is the product of the evidence ratio and the prior odds ratio . If every source model in the set is equally probable (namely, for ), 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 (): where is the true (albeit unknown) mean concentration, is the measurement error, is the predicted mean concentration obtained from an atmospheric dispersion model, 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 is assumed (for now) to be known and given by . Hence, it is assumed that we are given the concentration data and corresponding uncertainties 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 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 where is the mean residual for the 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 () as follows: Similarly, the uncertainty (error) in the residual datum is estimated as , where
Now, given the residual data , we evaluate the evidence for the presence of an additional source by considering two models (hypotheses) for : : there is no additional source and the residual data (where is the predicted concentration for a source model with no localized sources and hence has zero value), and : there is an additional source and the residual data (where is the predicted concentration for a source model consisting of one localized source). Under hypothesis (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 relative to . From (2.5) and assuming that hypotheses and are equally probable so , we get If exceeds some preassigned threshold, then there is evidence for an additional source (namely, hypothesis is favoured with respect to hypothesis ). 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 is the correct model structure. This analysis would involve drawing samples from the posterior distribution of (2.1) with replaced by , computing the residual data (obtained now by subtracting the predicted concentration from the measured concentration ), and then determining once more the evidence in favour of relative to that of . When this test fails (namely, 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 and associated estimated uncertainties .(2)Set and specify a threshold for the posterior odds ratio. Draw samples of source models (encoded by the parameter vectors , ) from the posterior distribution given by (2.1).(3)Compute the evidences and , and use these values to determine in accordance to (2.5). If , stop (no sources are present in the domain, so ).(4)Use samples () to estimate the residual data and the associated uncertainty in this data in accordance to (2.7) and (2.8), respectively.(5)Using the information from the previous step, compute in accordance to (2.9). If , set , output (), and stop.(6)Increase . Draw samples () 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 (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 ): 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 (). For a source model (or molecule) composed of localized sources (or source atoms), it is convenient to define the parameter vector as follows: .
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 where is the spatial-temporal filtering function for the sensor concentration measurement at and 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]: where is the Lagrangian position and 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), 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; are the increments of a vector-valued (three-dimensional) Wiener process; 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: where is the adjunct (or dual) concentration at space-time point associated with the concentration datum at location at time (with ). 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 ): 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: in (3.5) is related to in (3.3) as , and in (3.5) is related to in (3.3) as 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:
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: 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 () over which the evidence integral is carried out [cf. (2.2)]. More specifically, for , , where is an a priori upper bound on the emission rate; where is the a priori specified spatial domain that is assumed to contain the source ; and, and , where 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 , 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: where In (4.2), the notation for the likelihood function was adjusted to include the standard deviation for the noise in the conditioning to emphasize the fact that 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 () 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 ): 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 and variance . 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 where are the estimated standard deviations for the noise and . 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 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 ().
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 and . The assignment ensures that encoding our belief that our estimates of the standard deviation of are unbiased. Furthermore, the assignment 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 , 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 required in the evaluation of in (2.9), we are required to specify a functional form for the likelihood of the estimated mean residual data and the associated uncertainty in this data . To this end, the integrated likelihood function of (4.6) is also used for the specification of : Here, 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 for this test as where (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, specification of a method to draw samples from the posterior distribution (Steps 2 and 6) and specification of a method for computation of the evidence (Step 3, cf. (2.2)) or (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 and . Owing to the fact that the determination of and involves only an evaluation of a multidimensional integral in a low-dimensional space ( 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 (). 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 and , 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: where 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 () with and (), the evidence can be approximated from the likelihood-ordered samples as the following weighted sum:
If the prior mass points are sampled in a logarithmic manner as where 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 , , and (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 . Select the sample having the lowest likelihood (which we label as ) in the ensemble and remove (discard) it. Shrink the prior mass to .(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: .(5)If ( 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, where () 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 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 in (2.9) is a summary of the evidence for (additional source is present in the residuals) against (no source is present in the residuals). To interpret how strong (or weak) is the evidence for against provided by the residual data , we use a reference scale suggested by Jeffreys [30]. In this scale, corresponds to inconclusive evidence for against (or the evidence for 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 ). On this same scale, corresponds to very strong evidence for against (associated, as such, with a posterior odds ratio greater than about 150 to 1). In view of this scale for interpreting , 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 (or, equivalently, ).
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, , was estimated to be cm.
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: source 1 is at m; source 2 is at m; source 3 is at m; source 4 is at 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 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 m s−1, m , m, , , and . 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 , 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 . Furthermore, it is assumed that the height of the sources above ground level ( 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 includes the sensor sampling error variance, , in the measurement of the concentration datum and the model error variance, , in the prediction of . The two contributions are combined in quadrature to give the noise error variance . The measurement error standard deviation is estimated as 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, ). In consequence, the estimated noise error standard deviation is given by .
The algorithm described in Section 2 was used for the source reconstruction, applied to the concentration data obtained from 62 detectors in the array () (see Figure 1 which shows the locations of the 62 detectors as filled squares). The hyperparameters defining the prior have been chosen as follows: with g and with m for .
(a) N = 1
(b) N = 2
(c) N = 3
(d) N = 4
After Steps 2 and 3 of the algorithm were completed, we found implying that the concentration data strongly support the source model (one localized source or “signal’’) against the alternative source model (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 , 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) (), Figure 1(b) (), Figure 1(c) (), and Figure 1(d) (). The caption in Figure 1 summarizes the values for obtained from testing the alternative hypothesis (additional source present in the residual data) against the “null’’ hypothesis (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 with , providing the correct number of localized sources in the source distribution ( in this case).
The samples of source models () 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.
(a) Source 1
(b) Source 2
(c) Source 3
(d) Source 4
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 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 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 (), 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 ).