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

Authors: Eugene Yee
Journal: ISRN Applied Mathematics
 (2012)

---

## 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 ofN localized sources from the measured concentration data for increasing values of N 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.


---

## Body

## 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:
(1.1)θ^=argminθ∈ℝn{fO(θ)+fR(θ)},
where n is the number of parameters (unknowns) used to parameterize the source model, fO(θ) 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 fR(θ) 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 asingle 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, N, 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 unknowna 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 N and the parameters of the N 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 unknowna 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

LetmN denote a source distribution (model) which consists of a known number, N, of localized sources. The model mN is characterized by a set of parameters θ. Bayesian inference for θ, when given the concentration data 𝒟, is based on Bayes’ theorem (rule):
(2.1)p(θ∣mN,𝒟,I)=p(θ∣mN,I)p(𝒟∣mN,θ,I)p(𝒟∣mN,I),
where I 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), p(θ∣mN,𝒟,I)≡p(θ) is the posterior probability distribution of the parameters θ, p(θ∣mN,I)≡π(θ) is the prior probability distribution of the parameters θ before the data 𝒟 was made available, and p(𝒟∣mN,θ,I)≡ℒ(θ) is the likelihood function and defines the probabilistic model of how the data were generated. Finally, p(𝒟∣mN,I) is the evidence and, in the context of parameter estimation, ensures proper normalization of the posterior probability distribution (assuming θ∈ℝn):
(2.2)p(𝒟∣mN,I)=∫ℝn‍ℒ(θ)π(θ)dθ≡Z(mN).However, in our problem, the number of sources (N) is unknown a priori, so the relevant question that needs to be addressed is as follows: given a set of possible source models S≡{m0,m1,…,mNmax} (mN denotes the source model consisting of N 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 mN (N=0,1,…,Nmax) to give
(2.3)p(mN∣𝒟,I)=p(mN∣I)p(𝒟∣mN,I)p(𝒟∣I).
Here, p(mN∣I) is the prior probability of source model mN given only the information I, and p(𝒟∣mNI) is the global likelihood of the concentration data 𝒟 given the source model mN and I 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 Z(mN) defined in (2.2). Finally, the marginal probability of the data given the background information, p(𝒟∣I), is a normalization constant over all source models:
(2.4)p(𝒟∣I)=∑N=0Nmax‍p(mN∣I)p(𝒟∣mN,I).Given any two source modelsmN0 and mN1 from the set S, the “odds’’ in favor of model mN1 relative to model mN0 are given by the ratio of the posterior probabilities of the two models, which on using (2.3) and (2.2), reduces to
(2.5)K10≡p(mN1∣𝒟,I)p(mN0∣𝒟,I)=p(𝒟∣mN1,I)p(𝒟∣mN0,I)·p(mN1∣I)p(mN0∣I)=Z(mN1)Z(mN0)·p(mN1∣I)p(mN0∣I).
Note that the posterior odds ratio K10 is the product of the evidence ratio Z(mN1)/Z(mN0) and the prior odds ratio p(mN1∣I)/p(mN0∣I). If every source model in the set S is equally probable (namely, p(mN∣I)=1/(Nmax+1) for N=0,1,…,Nmax), 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 Z(mN) for each source model mN in S given the input concentration data 𝒟 and the background information I, 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 evidenceZ(mN) for an arbitrary value of N. A perusal of (2.2) shows that the evidence for the source model structure mN (namely, a source model consisting of exactly N 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 mN. 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 mN 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 ofZ(mN) (mN∈S). Unfortunately, this computation is a technically difficult problem for two reasons. Firstly, as N (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 n. More specifically as N increases, the dimensionality of the parameter space n increases linearly with N. Secondly, the integrand of the overlap integral of (2.2) becomes ever more complex with increasing N, 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 mN as N!. The factor N! here corresponds simply to the number of possible permutations of the labels for the localized sources in the source model mN. In summary, the increase in the dimensionality of the parameter space and in the complexity of ℒ(θ) with increasing N makes the problem of computation of Z(mN) technically difficult.In view of the difficulty arising from the numerical evaluation ofZ(mN) for increasing values of N, 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 mN 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 datadi (i=1,2,…,Nm):
(2.6)di=C¯i+ei(1)=C~i+ei(1)+ei(2)≡C~i+ei,
where C¯i is the true (albeit unknown) mean concentration, ei(1) is the measurement error, C~i is the predicted mean concentration obtained from an atmospheric dispersion model, ei(2) is the model error incurred by using the (necessarily) imperfect atmospheric model to predict C¯i, and ei is the composite error that includes both the measurement and model error. Furthermore, it is assumed that the expectation value 〈ei〉 of the composite error is zero, and the variance 〈ei2〉 is assumed (for now) to be known and given by σi2. Hence, it is assumed that we are given the concentration data 𝒟≡{d1,d2,…,dNm} and corresponding uncertainties {σ1,σ2,…,σNm} where σi is the standard deviation (square root of the variance) for the i-th concentration datum di. Normally, exact values for the uncertainties σi are unknown so that all we would have available would be some estimated values for these uncertainties si, rather than the true values σi.If we assumea priori that there are N localized sources in the domain, we can analyse the concentration data 𝒟 given that mN is the correct model structure and draw samples of source models with exactly N localized sources (encoded by θ) from the posterior distribution given by (2.1). Assume that Ns samples of source models with N localized sources have been drawn from the posterior distribution; namely, we have available the samples θk(mN)(k=1,2,…,Ns) where θk(mN) is kth source model sample consisting of exactly N localized sources. To determine if N is the correct number of localized sources, we can compute the residual data ℰ≡{ϵ1,ϵ2,…,ϵNm} where ϵi is the mean residual for the ith concentration datum (with realizations of the residuals obtained by subtracting the predicted concentration C~i(mN) for a sample of the source model mN encoded in θk(mN) from the measured concentration di). The mean residual datum ϵi is estimated from the ensemble of source model samples θk(mN) (k=1,2,…,Ns) as follows:
(2.7)ϵi^=1Ns∑k=1Ns‍(di-C~i(θk(mN))).
Similarly, the uncertainty siϵ (error) in the residual datum ϵi is estimated as siϵ=[ϵi2^-(ϵi^)2]1/2, where
(2.8)ϵi2^=1Ns∑k=1Ns‍(di-C~i(θk(mN)))2.Now, given the residual dataℰ, we evaluate the evidence for the presence of an additional source by considering two models (hypotheses) for ℰ: (1)  H0: there is no additional source and the residual data ϵi=C~i(m0)+ei=ei (where C~i(m0) is the predicted concentration for a source model with no localized sources and hence has zero value), and (2)  H1: there is an additional source and the residual data ϵi=C~i(m1)+ei (where C~i(m1) is the predicted concentration for a source model consisting of one localized source). Under hypothesis H0 (null hypothesis), the residuals are consistent with the estimated noise level ei. To decide which hypothesis is favoured by the data, we calculate the posterior odds ratio for H1 relative to H0. From (2.5) and assuming that hypotheses H0 and H1 are equally probable so p(H0∣I)=p(H1∣I)=1/2, we get
(2.9)K10ϵ≡p(H1∣ℰ,I)p(H0∣ℰ,I)≡p(m1∣ℰ,I)p(m0∣ℰ,I)=p(ℰ∣m1,I)p(ℰ∣m0,I)≡Zϵ(m1)Zϵ(m0).
If K10ϵ exceeds some preassigned threshold, then there is evidence for an additional source (namely, hypothesis H1 is favoured with respect to hypothesis H0). If this is the case, then the number of sources N is increased by one, and we repeat the analysis of the concentration data 𝒟 given that mN+1 is the correct model structure. This analysis would involve drawing samples from the posterior distribution of (2.1) with mN replaced by mN+1, computing the residual data ℰ (obtained now by subtracting the predicted concentration C~i(mN+1) from the measured concentration di), and then determining once more the evidence in favour of H1 relative to that of H0. When this test fails (namely, K10ϵ is less than or equal to the preassigned threshold value), the algorithm terminates. At this point, the number of sources N* has been determined, and the samples drawn from the posterior distribution of (2.1) for mN* 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𝒟={d1,d2,…,dNm} and associated estimated uncertainties {s1,s2,…,sNm}.(2)
SetN=1 and specify a threshold K* for the posterior odds ratio. Draw Ns samples of source models (encoded by the parameter vectors θk(mN), k=1,2,…,Ns) from the posterior distribution p(θ∣mN,𝒟,I) given by (2.1).(3)
Compute the evidencesZ(mN) and Z(mN-1), and use these values to determine K10 in accordance to (2.5). If K10≤K*, stop (no sources are present in the domain, so N*=0).(4)
Use samplesθk(mN) (k=1,2,…,Ns) to estimate the residual data ℰ={ϵ^1,ϵ^2,…,ϵ^Ns} and the associated uncertainty in this data {s1ϵ,s2ϵ,…,sNmϵ} in accordance to (2.7) and (2.8), respectively.(5)
Using the information from the previous step, computeK10ϵ in accordance to (2.9). If K10ϵ≤K*, set N*=N, output θk(mN) (k=1,2,…,Ns), and stop.(6)
IncreaseN→N+1. Draw Ns samples θk(mN) (k=1,2,…,Ns) of source models from p(θ∣mN,𝒟,I) 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 Z(mN) for N≥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 numberN of localized sources. In view of this, the model for the source density function has the following explicit form (for source model mN):
(3.1)𝒮(x,t)=∑j=1N‍Qjδ(x-xs,j)[U(t-Tb,j)-U(t-Te,j)],
where δ(·) and U(·) are the Dirac delta and Heaviside unit step functions, respectively, and, Qj, xs,j, Tb,j, and Te,j are the emission rate, vector position (location), source-on time, and source-off time, respectively, of the jth source (j=1,2,…,N). For a source model (or molecule) composed of N localized sources (or source atoms), it is convenient to define the parameter vector as follows: θ≡(xs,1,Tb,1,Te,1,Q1,…,xs,N,Tb,N,Te,N,QN)∈ℝ6N.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 C~ [cf. (2.6)]. The predicted concentration C~ can be compared directly to the measured concentration d “seen’’ by a sensor at receptor location xr and time tr, averaged over the sensor volume and sampling time. The comparison can be effected by averaging the mean concentration C(x,t) over this sensor volume and sampling time to give
(3.2)C~(xr,tr)≡∫0T‍dt∫ℛ‍dxC(x,t)h(x,t∣xr,tr)≡〈C∣h〉(xr,tr),
where h(x,t∣xr,tr) is the spatial-temporal filtering function for the sensor concentration measurement at (xr,tr) and ℛ×[0,T] corresponds to the space-time volume that contains the source distribution and the receptors (sensors). Note that the mean concentration C~(xr,tr) “seen’’ by a sensor can be expressed succinctly as the inner or scalar product 〈C∣h〉 of the mean concentration C and the sensor response function h.A source-oriented forward-time Lagrangian stochastic (LS) model can be used to predictC(x,t). In this approach, C is estimated from the statistical characteristics of “marked’’ particle trajectories described by the following stochastic differential equation [15]:
(3.3)dX(t)=V(t)dt,dV(t)=a(X(t),V(t),t)dt+(C0ϵ(X(t),t))1/2dW(t),
where X(t)≡(Xi(t))=(X1(t),X2(t),X3(t)) is the Lagrangian position and V(t)≡(Vi(t))=(V1(t),V2(t),V3(t)) is the Lagrangian velocity of a “marked’’ fluid element at time t (marked by the source as the fluid element passes through it at some earlier time t′). The state of the fluid particle at any time t after its initial release from the source density function 𝒮 is defined by (X,V). In (3.3), C0 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; dW(t)≡(dWi(t))=(dW1(t),dW2(t),dW3(t)) are the increments of a vector-valued (three-dimensional) Wiener process; a≡(ai)=(a1,a2,a3) 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 concentrationC~(xr,tr) can be computed alternatively using the following dual representation:
(3.4)C~(xr,tr)≡∫0T‍dt∫ℛ‍dxC(x,t)h(x,t∣xr,tr)=∫0trdt0∫ℛdx0C*(x0,t0∣xr,tr)𝒮(x0,t0)≡〈C*∣𝒮〉(xr,tr),
where C*(x0,t0∣xr,tr) is the adjunct (or dual) concentration at space-time point (x0,t0) associated with the concentration datum at location xr at time tr (with t0≤tr). A comparison of (3.2) and (3.4) implies the duality relationship 〈C∣h〉=〈C*∣𝒮〉 between C and C* through the source functions h and 𝒮. Moreover, C* is uniquely defined in the sense that it is explicitly constructed so that it verifies this duality relationship.The adjunct fieldC* 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 dt′>0):
(3.5)dXb(t′)=Vb(t′)dt′,(3.6)dVb(t′)=ab(Xb(t′),Vb(t′),t′)dt′+(C0ϵ(Xb(t′),t′))1/2dW(t′),
where at any given time t′, (Xb,Vb) 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 tr passed through the spatial volume of the detector at location xr). Here, (Xb,Vb) defines the state of a fluid particle at any time t′<tr before its “final” release from the sensor space-time volume h at time tr. It can be shown [15, 19] that C* obtained from (3.5) for a detector with the filtering function h and C obtained from (3.2) for a release from the source density 𝒮, are exactly consistent with the duality relationship 〈C∣h〉=〈C*∣𝒮〉 provided: (1)  Vb in (3.5) is related to V in (3.3) as Vb(t)=V(t), and (2)  ab in (3.5) is related to a in (3.3) as
(3.7)aib(x,u,t)=ai(x,u,t)-C0ϵ(x,t)∂∂uilnPE(u),
where PE(u) is the background Eulerian velocity PDF. Thomson [15] provides one particular solution for the drift coefficient vector a that is consistent with the well-mixed criterion for the case where PE(u) has a Gaussian functional form (namely, for Gaussian turbulence).Finally, if we substitute (3.1) into (3.4), the predicted concentration C~(θ)≡C~(xr,tr;θ) “seen’’ by the sensor at (xr,tr) is given explicitly by the following expression:
(3.8)C~(θ)=∑j=1N‍Qj∫Tb,jmin(T,Te,j)‍C*(xs,j,ts∣xr,tr)dts.


## 4.  Prior and Likelihood

### 4.1. Prior

Assuming that each parameter inθ (for source model mN) is logically independent, the prior π(θ) can be factored as follows:
(4.1)π(θ)=∏j=1Nπ(Qj)π(xs,j)π(Tb,j)π(Te,j).
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 ℝn (n≡6N) over which the evidence integral is carried out [cf. (2.2)]. More specifically, for j=1,2,…,N, Qj~𝒰([0,Qmax]), where Qmax is an a priori upper bound on the emission rate; xs,j~𝒰(ℛ) where ℛ⊂ℝ3 is the a priori specified spatial domain that is assumed to contain the source 𝒮; and, Tb,j~𝒰([0,Tmax]) and Te,j~𝒰([Tb,j,Tmax]), where Tmax 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 𝒰([a,b]) is the uniform distribution defined on the closed interval [a,b]. Finally, note in the specification of the prior for Te,j that the lower bound of the parameter range is Tb,j, 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 distributionei in (2.6), other than the fact that it has a finite variance σi2, 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:
(4.2)ℒ(θ)≡ℒ(θ∣𝒟,ς)=1∏i=1Nm‍2πσiexp(-12χ2(θ)),
where
(4.3)χ2(θ)≡∑i=1Nm‍(di-C~i(θ)σi)2.
In (4.2), the notation for the likelihood function was adjusted to include the standard deviation for the noise ς≡(σ1,σ2,…,σNm) in the conditioning to emphasize the fact that {σi}i=1Nm are assumed to be known quantities.Unfortunately, as mentioned previously, the noise termei in (2.6) is rather complicated, arising as such from a superposition of both model and measurement errors. In consequence, reliable estimates for σi (i=1,2,…,Nm) are difficult to obtain in practical applications. In view of this, let us denote by si the quoted estimate of the standard deviation for the noise term ei for which the true (but unknown) standard deviation is σi. Now, let us characterize the uncertainty in the specification of the standard deviation of ei with an inverse gamma distribution of the following form (or, equivalently, prior distribution for σi):
(4.4)π(σi∣si,α,β)=2αβΓ(β)(siσi)2βexp(-αsi2σi2)1σi,i=1,2,…,Nm,
where Γ(x) denotes the gamma function and α and β are scale and shape parameters, respectively, that define the inverse gamma distribution. The inverse gamma distribution for σi has mean 〈σi〉=siα1/2Γ(β-1/2)/Γ(β) and variance Var[σi]=si2α[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 σi of ei that appear in (4.2) can be treated as nuisance parameters and eliminated by considering the integrated likelihood
(4.5)ℒ(θ∣𝒟,s,α,β)=∫‍ℒ(θ∣𝒟,ς)π(ς∣s,α,β)dς=∫ℒ(θ∣𝒟,ς)∏i=1Nm2αβΓ(β)(siσi)2βexp(-αsi2σi2)1σidς,
where s≡(s1,s2,…,sNm) are the estimated standard deviations for the noise (e1,e2,…,eNm) and dς≡dσ1dσ2⋯dσNm. 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
(4.6)ℒ(θ∣𝒟,s,α,β)=∏i=1NmαβΓ(β+1/2)2πsiΓ(β)1[α+(di-C~i(θ))2/(2si2)]β+1/2.
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 ei (i=1,2,…,Nm).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 〈σi〉=si encoding our belief that our estimates si of the standard deviation of ei are unbiased. Furthermore, the assignment β=1 results in a very heavy-tailed distribution for π(σi∣si,α,β) which allows significant deviations of the noise uncertainty from the quoted value of si (provided by the user). Indeed, with the choice β=1, the variance associated with p(σi∣si,α,β) 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 si≪σi). 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 computeZϵ(m1) required in the evaluation of K10ϵ in (2.9), we are required to specify a functional form for the likelihood ℒ(θ(m1)∣ℰ,sϵ) of the estimated mean residual data ℰ≡{ϵ^1,ϵ^2,…,ϵ^Ns} and the associated uncertainty in this data sϵ≡{s1ϵ,s2ϵ,…,sNmϵ}. To this end, the integrated likelihood function of (4.6) is also used for the specification of ℒ(θ(m1)∣ℰ,sϵ):
(4.7)ℒ(θ(m1)∣ℰ,sϵ)≡ℒ(θ(m1)∣ℰ,sϵ,α,β)=∏i=1NmαβΓ(β+1/2)2πsiϵΓ(β)1[α+(ϵ^i-C~i(θ(m1)))2/(2siϵ2)]β+1/2.
Here, θ(m1) 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 Zϵ(m1) for this test as
(4.8)Zϵ(m1)=∫ℝn‍ℒ(θ(m1)∣ℰ,sϵ,α,β)π(θ(m1))dθ(m1),
where n=6 (corresponding to the dimension of the parameter space for one source).


## 4.1. Prior

Assuming that each parameter inθ (for source model mN) is logically independent, the prior π(θ) can be factored as follows:
(4.1)π(θ)=∏j=1Nπ(Qj)π(xs,j)π(Tb,j)π(Te,j).
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 ℝn (n≡6N) over which the evidence integral is carried out [cf. (2.2)]. More specifically, for j=1,2,…,N, Qj~𝒰([0,Qmax]), where Qmax is an a priori upper bound on the emission rate; xs,j~𝒰(ℛ) where ℛ⊂ℝ3 is the a priori specified spatial domain that is assumed to contain the source 𝒮; and, Tb,j~𝒰([0,Tmax]) and Te,j~𝒰([Tb,j,Tmax]), where Tmax 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 𝒰([a,b]) is the uniform distribution defined on the closed interval [a,b]. Finally, note in the specification of the prior for Te,j that the lower bound of the parameter range is Tb,j, 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 distributionei in (2.6), other than the fact that it has a finite variance σi2, 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:
(4.2)ℒ(θ)≡ℒ(θ∣𝒟,ς)=1∏i=1Nm‍2πσiexp(-12χ2(θ)),
where
(4.3)χ2(θ)≡∑i=1Nm‍(di-C~i(θ)σi)2.
In (4.2), the notation for the likelihood function was adjusted to include the standard deviation for the noise ς≡(σ1,σ2,…,σNm) in the conditioning to emphasize the fact that {σi}i=1Nm are assumed to be known quantities.Unfortunately, as mentioned previously, the noise termei in (2.6) is rather complicated, arising as such from a superposition of both model and measurement errors. In consequence, reliable estimates for σi (i=1,2,…,Nm) are difficult to obtain in practical applications. In view of this, let us denote by si the quoted estimate of the standard deviation for the noise term ei for which the true (but unknown) standard deviation is σi. Now, let us characterize the uncertainty in the specification of the standard deviation of ei with an inverse gamma distribution of the following form (or, equivalently, prior distribution for σi):
(4.4)π(σi∣si,α,β)=2αβΓ(β)(siσi)2βexp(-αsi2σi2)1σi,i=1,2,…,Nm,
where Γ(x) denotes the gamma function and α and β are scale and shape parameters, respectively, that define the inverse gamma distribution. The inverse gamma distribution for σi has mean 〈σi〉=siα1/2Γ(β-1/2)/Γ(β) and variance Var[σi]=si2α[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 σi of ei that appear in (4.2) can be treated as nuisance parameters and eliminated by considering the integrated likelihood
(4.5)ℒ(θ∣𝒟,s,α,β)=∫‍ℒ(θ∣𝒟,ς)π(ς∣s,α,β)dς=∫ℒ(θ∣𝒟,ς)∏i=1Nm2αβΓ(β)(siσi)2βexp(-αsi2σi2)1σidς,
where s≡(s1,s2,…,sNm) are the estimated standard deviations for the noise (e1,e2,…,eNm) and dς≡dσ1dσ2⋯dσNm. 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
(4.6)ℒ(θ∣𝒟,s,α,β)=∏i=1NmαβΓ(β+1/2)2πsiΓ(β)1[α+(di-C~i(θ))2/(2si2)]β+1/2.
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 ei (i=1,2,…,Nm).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 〈σi〉=si encoding our belief that our estimates si of the standard deviation of ei are unbiased. Furthermore, the assignment β=1 results in a very heavy-tailed distribution for π(σi∣si,α,β) which allows significant deviations of the noise uncertainty from the quoted value of si (provided by the user). Indeed, with the choice β=1, the variance associated with p(σi∣si,α,β) 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 si≪σi). 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 computeZϵ(m1) required in the evaluation of K10ϵ in (2.9), we are required to specify a functional form for the likelihood ℒ(θ(m1)∣ℰ,sϵ) of the estimated mean residual data ℰ≡{ϵ^1,ϵ^2,…,ϵ^Ns} and the associated uncertainty in this data sϵ≡{s1ϵ,s2ϵ,…,sNmϵ}. To this end, the integrated likelihood function of (4.6) is also used for the specification of ℒ(θ(m1)∣ℰ,sϵ):
(4.7)ℒ(θ(m1)∣ℰ,sϵ)≡ℒ(θ(m1)∣ℰ,sϵ,α,β)=∏i=1NmαβΓ(β+1/2)2πsiϵΓ(β)1[α+(ϵ^i-C~i(θ(m1)))2/(2siϵ2)]β+1/2.
Here, θ(m1) 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 Zϵ(m1) for this test as
(4.8)Zϵ(m1)=∫ℝn‍ℒ(θ(m1)∣ℰ,sϵ,α,β)π(θ(m1))dθ(m1),
where n=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 Section2, namely, (1) specification of a method to draw samples from the posterior distribution p(θ∣mN,𝒟,I) (Steps 2 and 6) and (2) specification of a method for computation of the evidence Z(m1) (Step  3, cf. (2.2)) or Zϵ(m1) (Step  5, cf. (4.8)).The most popular method that can be used for sampling fromp(θ∣mN,𝒟,I) (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 evidencesZ(m1) and Zϵ(m1). Owing to the fact that the determination of Z(m1) and Zϵ(m1) 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 ℝn (n≫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 fromp(θ∣mN,𝒟,I) and to evaluate Z(m1) and Zϵ(m1), 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 Z for the general case. The nested sampling algorithm transforms the multidimensional evidence integral of (2.2) into a simple one-dimensional representation:
(5.1)Z≡Z(mN)=∫01‍ℒ(χ)dχ,
where
(5.2)χ(ℒ*)=∫ℒ(θ)>ℒ*‍π(θ)dθ
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 m points χi (i=1,2,…,m) with 0<χm<χm-1<⋯<χ2<χ1<1 and ℒ(χi)>ℒ(χj) (i>j), the evidence Z can be approximated from the likelihood-ordered samples as the following weighted sum:
(5.3)Z^=∑i=1m‍ℒ(χi)δχi,δχi=χi-1-χi.If the prior mass pointsχi are sampled in a logarithmic manner as χi=∏j=1i‍tj where tj∈(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)
Seti=0, χ0=1, Z0=0 and f=0.5 (preset fraction used in the stopping criterion). Randomly draw M samples θ from the prior π(θ) to give an ensemble of samples. Evaluate the likelihood ℒ(θ) for each of the samples in the ensemble.(2)
Increasei→i+1. Select the sample having the lowest likelihood (which we label as ℒi) in the ensemble and remove (discard) it. Shrink the prior mass to χi=χi-1e-1/M.(3)
Draw a new sampleθ from the prior π(θ) subject to the hard likelihood constraint ℒ(θ)>ℒi, and add this sample to the ensemble of samples.(4)
Increment the evidence:Zi=Zi-1+ℒi(χi-1-χi).(5)
Ifℒmaxχi<fZi (ℒmax is the largest value of the likelihood in the ensemble of samples), add in contributions to the evidence Zi from the ensemble of samples (remaining M 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 evidenceZ^ is completed by adding the contribution of the remaining M samples in the ensemble to Zi; namely, Z^→Zi+M-1(ℒ(θ1)+ℒ(θ2)+⋯+ℒ(θM))χi where ℒ(θk) (k=1,2,…,M) are the likelihood values evaluated at the remaining M samples. It is noted that the algorithm summarized above for estimation of Z (overlap integral of ℒ(θ) and π(θ)) automatically provides weighted samples θk(mN) drawn from the posterior distribution p(θ∣mN,𝒟,I)∝ℒ(θ)π(θ) (cf. (2.1) and (2.2)). More specifically, the ith discarded sample in Step  2 of the algorithm can be interpreted as a sample drawn from the posterior distribution of θ with weight given by pi=ℒi(χi-1-χi)/Z^ where Z^ 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 ratioK10ϵ in (2.9) is a summary of the evidence for H1 (additional source is present in the residuals) against H0 (no source is present in the residuals). To interpret how strong (or weak) is the evidence for H1 against H0 provided by the residual data ℰ, we use a reference scale suggested by Jeffreys [30]. In this scale, log(K10ϵ)<1 corresponds to inconclusive evidence for H1 against H0 (or the evidence for H1 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 H1). On this same scale, log(K10ϵ)>5 corresponds to very strong evidence for H1 against H0 (associated, as such, with a posterior odds ratio greater than about 150 to 1). In view of this scale for interpreting K10ϵ, we choose the threshold K* for the posterior odds ratio (required to be specified in Step  2 of the algorithm summarized in Section 2) to be log(K*)=1 (or, equivalently, K*≈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,z0, was estimated to be z0=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, zd, 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,zs, of 2.0 m: (1) source 1 is at (xs,ys)=(33.0,171.0) m; (2) source 2 is at (xs,ys)=(33.8,240.7) m; (3) source 3 is at (xs,ys)=(30.0,312.9) m; (4) source 4 is at (xs,ys)=(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 speedS2 at the 2 m level, the friction velocity (u*), the atmospheric stability (Obukhov length L), and the standard deviations of the velocity fluctuations in the alongwind (σu), crosswind (σv), and vertical (σw) directions were S2=3.61 m s−1, u*=0.282m s-1, L=-27.3m, σu/u*=2.33, σv/u*=1.86, and σw/u*=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 +x-direction and, hence, is perpendicular to the sampling lines of detectors along the y-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 C~ (cf. (3.4)). More specifically, the backward-time LS model of (3.5) was applied with the Kolmogorov constant C0=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 modelmN are θ≡θ(mN)=(xs,1,Q1,…,xs,N,QN). Furthermore, it is assumed that the height of the sources above ground level (zs=2.0 m) is known a priori, so the only unknown location parameters are the (xs,ys) coordinates of the sources. Owing to the horizontal homogeneity of the mean flow and turbulence statistics in the current example, the adjunct concentration C*(xs∣xd) in (3.4) can be precalculated for one detector location xd (for the known source height zs), with the adjunct concentration C* (considered as a function of xs) at all other detector locations obtained simply by a linear translation of C*(xs∣xd) in the horizontal (x,y)-plane.In order to calculate the likelihood function given by (4.6), we need to provide values for the estimated standard deviations si for the noise ei. The noise error variance σi2 includes the sensor sampling error variance, σd,i2, in the measurement of the concentration datum di and the model error variance, σm,i2, in the prediction of C~i. The two contributions are combined in quadrature to give the noise error variance σi2=σd,i2+σm,i2. The measurement error standard deviation is estimated as σ^d,i≡sd,i=max(0.05,0.02di) 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 C~i(θ) (namely, σ^m,i≡sm,i=0.20C~i(θ)). In consequence, the estimated noise error standard deviation is given by σ^i≡si=(sd,i2+sm,i2)1/2.The algorithm described in Section2 was used for the source reconstruction, applied to the concentration data 𝒟≡{d1,d2,…,dNm} obtained from 62 detectors in the array (Nm=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: Qj~𝒰([0,Qmax]) with Qmax=100g s-1 and (xs,j,ys,j)~𝒰(ℛ) with ℛ=[0,100]m×[0,500] m for j=1,2,…,N.Four iterations of the algorithm of Section2 showing the samples of source model mN drawn from the posterior distribution p(θ∣mN,𝒟,I) and the computation of K10ϵ for (a) N=1, log(K10ϵ)=62.4±0.2; (b) N=2, log(K10ϵ)=57.3±0.3; (c) N=3, log(K10ϵ)=39.1±0.3; (d) N=4, log(K10ϵ)=0.92±0.33. A solid square denotes a detector.
(a)
N = 1
(b)
N = 2
(c)
N = 3
(d)
N = 4After Steps 2 and 3 of the algorithm were completed, we foundlog(K10)=log(Z(m1)/Z(m0))=64.0±0.2 implying that the concentration data strongly support the source model m1 (one localized source or “signal’’) against the alternative source model m0 (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(K10)>log(K*)=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) (N=1), Figure 1(b) (N=2), Figure 1(c) (N=3), and Figure 1(d) (N=4). The caption in Figure 1 summarizes the values for log(K10ϵ)=log(Zϵ(m1)/Zϵ(m0)) obtained from testing the alternative hypothesis H1 (additional source present in the residual data) against the “null’’ hypothesis H0 (no additional source present in the residual data). Each panel in Figure 1 also displays a density plot of samples of the source model mN drawn from the posterior distribution p(θ∣mN,𝒟,I) and projected onto the (x,y) horizontal plane (namely, each point in the plot corresponds to the (x,y) location of a sample of the source model mN drawn from the posterior distribution). Note that the algorithm terminates with N*=N=4 with log(K10ϵ)=0.92≤log(K*)=1, providing the correct number of localized sources in the source distribution (N*=4 in this case).The samples of source modelsmN (N=4) drawn from the posterior distribution p(θ∣mN,𝒟,I) 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.Table 1
The posterior mean, posterior standard deviation, and lower and upper bounds of the 95% HPD interval of the parametersxs,j (m), ys,j (m), and Qj≡qs,j (g s-1) for j=1, 2, 3, and 4 calculated from samples of source model mN with N=4 drawn from the posterior distribution p(θ|mN,𝒟,I) [cf. Figure 1(d)].
Parameter
Mean
Standard deviation
95% HPD
Actual
x
s (m)
34.9
4.3
(27.6, 44.4)
33.0
j
=
1
y
s (m)
170.8
0.6
(169.5, 171.8)
171.0
q
s (g s−1)
8.9
0.4
(8.0, 9.3)
—
x
s (m)
33.7
1.1
(31.6, 35.9)
33.8
j
=
2
y
s (m)
240.7
0.3
(240.2, 241.4)
240.7
q
s (g s−1)
7.5
0.4
(6.7, 8.1)
—
x
s (m)
25.4
3.9
(17.1, 31.2)
30.0
j
=
3
y
s (m)
313.2
0.6
(311.9, 314.6)
312.9
q
s (g s−1)
4.2
0.2
(3.7, 4.6)
3.8
x
s (m)
25.1
3.0
(
18.9
,
30.3
)
26.0
j
=
4
y
s (m)
384.5
0.5
(383.5, 385.5)
384.4
q
s (g s−1)
6.3
0.3
(5.8, 6.8)
—Univariate (diagonal) and bivariate (off-diagonal) marginal posterior distributions of the source parameters, namely, alongwind locationxs, crosswind location ys, and emission rate qs that characterize sources 1 (a), 2 (b), 3 (c), and 4 (d) obtained from samples of source model mN (N=4) drawn from the posterior distribution p(θ∣mN,𝒟,I) ([cf. Figure 1(d)).
(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 evidenceZ(mN) for N>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 mN 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 Zϵ(mN) for N=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 p(θ∣mN,𝒟,I) 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 (N*=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 mN of different dimensions involving different numbers of sources N).



---

*Source: 465320-2012-08-29.xml*