Contaminant Source Identification from Finite Sensor Data: Perron–Frobenius Operator and Bayesian Inference

Abstract

Sensors in the built environment ensure safety and comfort by tracking contaminants in the occupied space. In the event of contaminant release, it is important to use the limited sensor data to rapidly and accurately identify the release location of the contaminant. Identification of the release location will enable subsequent remediation as well as evacuation decision-making. In previous work, we used an operator theoretic approach—based on the Perron–Frobenius (PF) operator—to estimate the contaminant concentration distribution in the domain given a finite amount of streaming sensor data. In the current work, the approach is extended to identify the most probable contaminant release location. The release location identification is framed as a Bayesian inference problem. The Bayesian inference approach requires considering multiple release location scenarios, which is done efficiently using the discrete PF operator. The discrete PF operator provides a fast, effective and accurate model for contaminant transport modeling. The utility of our PF-based Bayesian inference methodology is illustrated using single-point release scenarios in both two and three-dimensional cases. The method provides a fast, accurate, and efficient framework for real-time identification of contaminant source location.

1. Introduction

The accidental or intentional release of pollutants in the built environment presents a serious risk factor to occupant health and safety. In addition to common pollutants—like CO${}_2$, volatile organic chemical compounds (VOCs), atmospheric particulate matter (pollens), and microbial contaminants—which affect indoor air quality (IAQ) [1], maliciously released pollutants in enclosed public spaces can cause human fatalities. Such events include release of chemical and biological weapons (CBW) or the transmission of infectious diseases (TID) with examples being the sarin gas attack in the Tokyo subway system in 1995 [2,3], spread of influenza in aircraft [4], spread of SARS and COVID virus [5,6,7] and outbreaks of measles and tuberculosis infection in offices and schools [8,9].

In all these release scenarios, an early detection of the location of the release of the hazardous agent is important for appropriate mitigation to ensure occupant safety. The sensor network plays a critical role in providing streaming data that can be used to identify release location and intensity of the contaminant. Additionally, this data can be used for getting accurate estimates of the spatial distribution of pollutant [10], and for predicting where the pollutant will spread next. This information can help determine evacuation pathways, as well as for triaging areas in terms of the extent of contamination, and prioritizing cleanup locations. The problem of estimating the unknown release location using the sensor measurements is characterized as an inverse modeling problem.Such inverse problems for estimating pollutant release locations occur frequently in the fields of atmospheric sciences [11] and hydrology [12]. Based on the available solution methods for the inverse problem in the literature, the methods can be classified into three categories as; direct methods, optimization-based methods, and probability based methods. The direct model evolves the contaminant concentration back in time to compute the source location and intensity at the time of release. This method is particularly suitable in cases where the source location and release time are known and only source flux is required to be computed. A detailed discussion of this method is provided in work by Zhang and Chen [13,14], where they introduce quasi-reversible and pseudo-reversible approaches for back-propagation. While the approach does not require a lot of prior information, it does suffer from numerical instabilities.

In optimization-based methods, a set of possible source locations are considered. The forward airflow and contaminant transport analysis are performed with this set of possible scenarios. The difference between the monitored observation data and simulated data is then used to identify the most plausible scenario. This is usually framed as a minimization (i.e., optimization) of carefully formulated objective functions to identify the source location. A detailed discussion of this approach is provided in Zhang et al. [15], where they especially use a regularization term to quantify a continuously releasing source. In [16], this approach was extended to a unified inverse modeling framework that accounts for sensor alarming time for source identification. To overcome the computational burden of evaluating the forward airflow and contaminant transport analysis multiple times during the optimization, the authors used a response factor method to construct a linear representation which is primarily valid for impulse release scenarios. Following a similar approach, Wei et al. [17] integrated a regularization framework with a Bayesian method to identify the source location. They assumed the availability of concentration measurements at every discrete location in the domain to ensure that the system is invertible for calculating the release rate. Another optimization method—variational continuous assimilation (VCA)—widely used by the meteorology community was utilized by Matsuo et al. [18,19]. The method assumes a steady-state concentration field which only holds when the release is constant. The method utilizes computational fluid dynamics (CFD) computations as the forward model to generate the data, which complicates real-time application due to the high computational cost. Closely related to the ’classical’ optimization-based methods, data-driven approaches based on artificial neural networks (ANN) were also investigated in [20,21]. The approach elegantly integrated machine learning concepts with a multi-zonal modeling approach. The multi-zone model, however, is based on the well-mixed assumption which limits the precise identification of contaminant sources in the building. More recent work tries to relax this assumption by using surrogate models machine learning concepts (Gaussian models, ANN, deep learning) trained using data from CFD simulations [22,23]. Care has to be chosen to carefully select the dataset to account for all possible scenarios which is a non-trivial task.

Finally, the probability-based methods can be further classified into traditional probability methods and adjoint probability methods. The traditional methods run multiple forward simulations sampled from a distribution of all possible source release locations. Bayesian inference (or variants) is used to compute the likelihood of pollutant source using conditional probability estimates. Sohn et al. [24,25] proposed the Bayesian formulation and extended it for various applications including optimal sensor system design [26,27,28]. Most of these results were showcased by utilizing multi-zonal models, which traditionally make a well-mixed air flow assumption. These assumptions limit how precisely one can localize the contaminant locations as the spatial resolution is limited to zones. This is, however, not a serious limitation since increasing the number of zones (i.e., reducing the size of individual zones) can potentially improve the resolution of these methods. The second kind of probabilistic approach, adjoint probability framework [12,29] was first used for groundwater/hydrology research problems. The method was adopted for indoor pollutant source problem by Lui and Zhai [30] for a multizone office building with known source release time. As an extension, Lui and Zhai [31] derived the adjoint equations for CFD and demonstrated its ability for two-dimensional office space and a three-dimensional aircraft cabin. A similar approach was also used by Zhai et al. [32] for the case of a continuous release scenario, where the total release mass of the contaminant is known a priori. In a recent work by Wang et al. [33] the approach is also extended for a dynamic airflow field.

One of the common limitations of the above-discussed methods categorized under probabilistic/optimization/inverse approaches that preclude them from the widespread usage is that they use computationally intensive forward models which rely on solving PDE using numerical modeling. These PDE solves are time-intensive and costly. One way to circumvent this constraint is by using multi-zonal models and CFD-multizone coupled models. However, the well-mixed assumption associated with zonal models precludes fine resolution of the spatial contaminant distribution. Therefore, they are limited to locating the contaminant source up to a zonal resolution of the building. Estimating release location more precisely requires the contaminant transport analysis model to provide more resolved contaminant distribution in the domain, these distributions also serve the inverse modeling framework.

In this context, there have been recent developments in data-driven approaches as well as Perron–Frobenius (PF) operator-based approaches [34,35,36,37] which open up the possibility of a fast, robust and data-driven methodology for performing contaminant transport analysis. This work particularly focuses on a PF approach which has shown utility in designing an optimal sensor network for monitoring indoor air quality [34,35]. The PF approach essentially transforms the problem of contaminant transport into a problem of simple matrix-vector products. The objective of this paper is to use the fast PF operator as a forward model for identifying indoor contaminant source release location. The robustness and accuracy of the PF approach for quickly computing the contaminant evolution for any arbitrary release in the space allow deploying a probabilistic framework, specifically Bayesian inference (Bayes Monte Carlo) technique with sequential updating to solve for the release location. The Bayesian inference approach has been successfully deployed for various applications such as risk assessment and water quality application [38,39]. It was applied to indoor source identification problems by Sohn et al. [24], where the authors utilize a two-stage approach to effectively and easily implement the method for real-time risk assessment. With the increase of edge computing devices (IoT), Bayesian inference can be done efficiently. In addition the approach provides the uncertainty estimates associated with the predictions. The estimates of uncertainty plays an important role in real-time to narrow the search space, which describes the versatility of the approach. The PF operator-based contaminant transport can further simplify this two-stage strategy to a single-stage as the contaminant distribution for any arbitrary release can be simulated instantly, once the PF operator is constructed. Another important property of the PF operator is the linearity. We, however, do not explicitly exploit this property in this paper. Additionally, the PF operator can also provide a systematic approach for sensor layout design as well as for risk assessment, thus providing a unified approach for tackling risk assessment, identification, and mitigation for IAQ, CBW, and TID applications.

The objective of this work is to use the real-time data coming from a limited number of sensors for estimating the source location. We consider scenarios where a single stationary contaminant source is released for a finite time. The linearity property of the PF operator can be exploited to construct an algorithm for source localization with multiple release locations without the associated combinatorial explosion in simulations needed for the multiple release scenario. However, in this paper we restrict the study to the case of a single contaminant release source and leave the case of multiple release scenarios for subsequent work. The sensors are placed optimally based on the PF operator-based sensor placement algorithm developed in Fontanini et al. [34]. The implementation of the approach is shown for both two and three-dimensional problems, with the two-dimensional results shown for a typical office space with a manikin, and the approach illustrated for three-dimensional problems using a typical furnished room space. The outline for the paper is as follows. We discuss the method for constructing the transfer PF operator for the contaminant transport in Section 2.2, the approach for source identification based on the Bayesian approach is discussed in Section 2.3. The results for the problem in 2D and 3D are presented in Section 3. Discussions and conclusions are detailed in Section 4 and Section 5.

2. Methodology

2.1. Problem Definition

The objective is to solve the inverse problem of identifying the contaminant source release location in the domain $\mathsf{\Omega }\subset {\mathbf{R}}^d$ using the contaminant monitoring sensor network measurements. The problem can be posed in different ways such as (a) identifying release duration, intensity, and location, (b) identifying release duration and location (c) identifying only release location. Out of the three cases, identifying a location is of major importance, as it helps in the containment of the hazardous compound release. Therefore, the current approach focuses on estimating the contaminant release location for a contaminant release occurring constantly from an unknown release location. We assume that this domain is outfitted with a sensor network that provides observations every time-step ${\tau }_{obs}$. We next briefly review the construction of a discrete form of Perron Frobenius operator which is used as a forward model for evaluating the concentration distribution for an arbitrary release in the domain.

2.2. Construction of Transfer PF Operator for Contaminant Transport

The PF operator is a linear albeit infinite-dimensional operator which captures the non-linearity of a dynamical system [40,41,42]. Once constructed it can be used for evolving the state of the system from a given initial condition to another time step. The operator can be constructed using a set-based [43] or a data-driven approach [44,45]. The discrete form of the PF operator is called the Markov matrix. We next discuss how this Markov matrix is constructed using the contaminant transport equation.

The modeling of the contaminant transport is done by a transient advection-diffusion partial differential equation given as

$$\frac{\partial \mathsf{\Phi }}{\partial t}+\nabla \left(U\mathsf{\Phi }\right)+{\nabla }^2\left(D\mathsf{\Phi }\right)=S_{\mathsf{\Phi }}$$

(1)

The contaminant density, $\mathsf{\Phi }(X,t)$ is the contaminant at spatial location (X) and time (t), is propagated by the air flow field, U in a domain $\mathsf{\Omega }\subset {\mathbf{R}}^d$. The flow field can be generated either experimentally, or computationally using computational fluid dynamics (CFD). Here D is the diffusion constant and $S_{\mathsf{\Phi }}$ is the source term. The velocity flow field, U, can be (un)steady. To model the sensor measurements we define Equation (2) with ${\chi }_{A_k}\left(x\right)$ for $k=1,\dots ,\tilde{p}$ denoting the indicator function for a set $A_k\subset X$ that corresponds to the location of p number of sensors in the domain as;

$$y={\chi }_{A_k}\left(x\right)\mathsf{\Phi },k=1,\dots ,\tilde{p}$$

(2)

The contaminant evolution Equation (1) is numerically solved by typically discretizing the domain spatially and temporally. In the operator setting the numerical scheme used for solving Equation (1) can be seen as a discrete-time equivalent operator given as

$$\begin{array}{c}\hfill {\mathsf{\Phi }}_{t+\delta t}=L\left({\mathsf{\Phi }}_t\right)\end{array}$$

(3)

Therefore, the continuous space-time evolution of the contaminant transport can be replaced by a discrete-time counterpart using linear transfer PF operator. Equation (1), in the absence of source term in the PF operator setting can be described by following

$$\begin{array}{c}\hfill \left[\mathbb{P}\mu \right]\left(A\right)={\int }_{X}p(x,A)d\mu \left(x\right).\end{array}$$

(4)

where $\mu \in \mathcal{M}\left(X\right)$ the space of real-valued measure and $p(x,A)$ is the stochastic transition function and describes the transition probability from point x to set $A\subset X$. To construct this finite-dimensional approximation of the PF operator (Markov matrix), the nonlinear flow field U and the diffusion coefficient D are required. The example of this is shown in Figure 1, where the velocity field is shown by vectors and the color contours as the initial concentration.

The PF operator ($\mathbf{P}$) is defined on the discrete representation of the space X. The space, X, is discretized into a finite number of ${\omega }_k$ cells/states for $k=1,\dots ,N$. The time evolution in this finite dimension discrete setting is given by the dynamical linear system model form in Equation (5).

$$\begin{array}{c}\hfill \mathbb{M}({\mu }_0,\widehat{S},\mathbf{P},\mathbf{C})\left\{\begin{array}{cc}& {\mu }_{t_{i+1}}=\mathbf{P}{\mu }_{t_i}+{\widehat{S}}_{t_i,t_{i+1}}i\in \{0,\cdots ,m\}\hfill \\ & y=\mathbf{C}{\mu }_{t_i}\hfill \end{array}\right.\end{array}$$

(5)

where ${\mu }_{t_i}\in {\mathbf{R}}^N$ as the discrete form of the scalar field at given time $t_i$ and $\mathbf{C}\in {\mathbf{R}}^{N\times \tilde{p}}$ is the observation matrix mapping that maps the state vector to the observation space. The source term, ${\widehat{S}}_{t_i,t_{i+1}}\in {\mathbf{R}}^N$, is vector representation which includes volumetric and inlet sources in the domain $\mathsf{\Omega }\subset {\mathbf{R}}^d$ [34,43]. The ${\mu }_{t_i}$ is defined as the cell volumetric average of $\mathsf{\Phi }$ at a given time $t_i$, with $V_{{\omega }_k}$ as the volume of the $k^{th}$ state and is given by

$${\mu }_{t_i}\left({\omega }_k\right)=\frac{1}{V_{{\omega }_k}}\int \int {\int }_{{\omega }_k}\mathsf{\Phi }(x,y,z,t_i)dV,{\omega }_k=1,\cdots ,N$$

(6)

For a transition matrix $\mathbf{P}$, each row k represents the transition probability of state ${\omega }_k$ into other states in the next time step. These transition values are obtained by releasing a normalized concentration of 1.0 from state ${\omega }_k$ and then computing the spread of this initial concentration to the rest of the states in the next time step. The time evolved concentrations are used to populate each row of the Markov matrix as Equation (7).

$$P_{(D_k,D_k+1)}(i,j)={\mu }_{t_i(D_k+1)}\left(j\right)i=1:N$$

(7)

Additional details on the construction of this matrix for large Courant numbers and unsteady flow fields are provided in Fontanini et al. [34], Fontanini et al. [43]. Once the $\mathbf{P}$ matrix is constructed, it can be used for propagating any initial contaminant distribution in time very efficiently simply by matrix-vector multiplication.

2.3. Bayesian Formulation

The Bayesian inference framework poses the source identification problem in a probabilistic setup. The release location is considered as a random variable rather than a constant. As a result, Bayesian inference not only produces the most likely point estimate of the source location but provides a probability distribution of the release location. This feature provides a rigorous framework for quantifying the uncertainties associated with the prior information of the random variable. The method is based on using available information to provide updated estimates of the source release location. This is called the posterior distribution which is computed using the Bayes rule as

$$p\left({\mathbf{X}}_{\mathbf{s}}\right|y)=\frac{p\left(y\right|{\mathbf{X}}_s)p\left({\mathbf{X}}_s\right)}{p\left(y\right)}\propto p\left(y\right|{\mathbf{X}}_s)p\left({\mathbf{X}}_s\right)$$

(8)

The factor $p\left(y\right|{\mathbf{X}}_s)$, is called the likelihood function. The likelihood represents a relative agreement of the observation y, given some known parameter value, or hypothesis, ${\mathbf{X}}_s$. For the source release estimation application, $p\left(y\right|{\mathbf{X}}_s)$ represents the likelihood of observing a set of measurements based on the modeled release scenarios ${\mathbf{X}}_s$. Bayes’ rule provides a means to estimate the inverse probability, $p\left({\mathbf{X}}_{\mathbf{s}}\right|y)$, which is the probablity of hypothesis, ${\mathbf{X}}_s$ given data, y. Another factor in Bayes rule is $p\left({\mathbf{X}}_s\right)$, called as the prior probability. For the current application, $p\left({\mathbf{X}}_s\right)$ provide an assessment of release occurring at particular location, before the release actually occurs. This probability is estimated before actual observations are available. The denominator of Bayes rule is the probability of observing a particular outcome, also called as normalizing term. For a discrete case this is calculated as follows: $p\left(y\right)={\sum }_{{\mathbf{X}}_s}p\left({\mathbf{X}}_s\right)p\left(y\right|{\mathbf{X}}_s)$. In continuous case it is defines as $p\left(y\right)=\int p\left({\mathbf{X}}_s\right)p\left(y\right|{\mathbf{X}}_s)d{\mathbf{X}}_s$. Because $p\left(y\right)$ is a constant with fixed ${\mathbf{X}}_s$ that normalizes the numerator, a proportionality is often used.

Bayesian inference depends on the data to make an improved posterior estimate. With the sensor measurements becoming available at every observation time-step the prior belief can be updated using the current posterior. This is also called Bayesian sequential updating. Consider data arriving sequentially $y_1,\cdots ,y_m,\cdots$, and we seek to update our inference for an unknown ${\mathbf{X}}_s$ online. In the Bayesian framework we have a prior distribution $p\left({\mathbf{X}}_s\right)$ and at time m we have a density for data conditional on ${\mathbf{X}}_s$ as

$$p(y_1,\cdots ,y_m|{\mathbf{X}}_s)=p\left(y_1\right|{\mathbf{X}}_s)p\left(y_2\right|y_1,{\mathbf{X}}_s)\cdots p\left(y_m\right|y_{m-1},{\mathbf{X}}_s)$$

(9)

where we do not assume $y_1,\cdots ,y_m,\cdots$ to be independently conditioned on ${\mathbf{X}}_s$. At time m we can update estimate of ${\mathbf{X}}_s$ as the posterior

$$p_m\left({\mathbf{X}}_s\right)=p\left({\mathbf{X}}_s\right|y_m)\propto p\left(y_m\right|{\mathbf{X}}_s)p\left({\mathbf{X}}_s\right)$$

(10)

With the next observation $y_{m+1}$, we can either start afresh as

$$p_{m+1}\left({\mathbf{X}}_s\right)=p\left({\mathbf{X}}_s\right|y_{m+1})\propto p\left(y_{m+1}\right|{\mathbf{X}}_s)p\left({\mathbf{X}}_s\right)$$

(11)

or since we know our prior belief of ${\mathbf{X}}_s$, before $m+1$ is assimilated into $p_{m+1}\left({\mathbf{X}}_s\right)$, we just use this as our prior distribution for the new piece of information. This is shown in Equation (12a) and is equivalent to Equation (11) as shown below.

$$\begin{array}{cc}\hfill {\widehat{p}}_{m+1}\left({\mathbf{X}}_s\right)& \propto p\left(y_{m+1}\right|y_m,{\mathbf{X}}_s)p_m\left({\mathbf{X}}_s\right)\hfill \end{array}$$

(12a)

$$\begin{array}{cc}\hfill & \propto p\left(y_n\right|y_m,{\mathbf{X}}_s)p\left(y_{m+1}\right|y_m,{\mathbf{X}}_s)p\left({\mathbf{X}}_s\right)\hfill \end{array}$$

(12b)

$$\begin{array}{cc}\hfill & =p\left(y_{m+1}\right|y_m,{\mathbf{X}}_s)p\left({\mathbf{X}}_s\right)\propto p_{m+1}\left({\mathbf{X}}_s\right)\hfill \end{array}$$

(12c)

To compute the posterior depending on the prior distribution and the likelihood distribution the solutions can be calculated analytically [46]. However, in the absence of an analytical solution, numerical modeling approaches are used. We use the "Bayes Monte Carlo updating" approach in this work. In this approach, samples are drawn from the assumed prior distribution for the model input, which is then used to compute the model outputs. A finite set of samples are chosen and are exhaustive enough to compute the denominator accurately. In the discrete representation of the PF operator, each state can act as a release scenario. To avoid an explosion in the number of samples, we choose to sample the release location by uniformly covering the complete domain. For these samples, we generate an epsilon ball radius of r to define release regions with more than one state in them. The release concentration is normalized in the region.

To generate n sample region in the domain $\mathsf{\Omega }\subset {\mathbf{R}}^d$ as release region $\mathcal{R}=\{{\widehat{S}}_1,\cdots ,{\widehat{S}}_n\}$ with epsilon ball of radius $ϵ=r$, the Halton quasi-random sample are generated in ${\mathbf{R}}^d$. Each release scenario is equally probable and therefore we define a uniform distribution for set $\mathcal{R}$. For these samples, the model sensor measurements can be pre-computed or can be computed online using the PF operator contaminant transport model $\mathbb{M}$ Equation (5). This measurements for individual realization can be written as ${\widehat{y}}_m^n=\mathbb{M}({\mu }_0,{\widehat{S}}_n,\mathbf{P},\mathbf{C})$. Considering these n samples the Bayes Equation (8) can be modified as

$$p_m\left({\widehat{S}}_n\right|y_m)=\frac{p_m\left(y_m\right|{\widehat{y}}_m^n)p_m\left({\widehat{S}}_n\right)}{{\sum }_{i=1}^n{p}_m\left(y_m\right|{\widehat{y}}_m^i)p_m\left({\widehat{S}}_i\right)}$$

(13)

where $p_m\left({\widehat{S}}_n\right|y_m)$ is the posterior probability of the $n^{th}$ realization condition on the measurement $y_m$, $p_m\left(y_m\right|{\widehat{y}}_m^n)$ is the likelihood of observing $y_m$ given the $n^{th}$ model measurements, $p_m\left({\widehat{S}}_n\right)$ prior probability of release in the $n^{th}$ region. For each realization Equation (13) can be used as an update equation as discussed in Equation (12) whenever new sensor measurements become available.

In Equation (13), the likelihood function $p_m\left(y_m\right|{\widehat{y}}_m^n)$ quantifies the model to measurement error. For an unbiased measurements with a normally distributed error the likelihood of observing sensor measurements y given the model prediction of the measurements as ${\widehat{y}}_m^n$ can be written as

$$p\left(y_m\right|{\widehat{y}}_m^n)=\frac{1}{\sigma \sqrt{2\pi }}exp\left(-\frac{1}{2}{\left[\frac{y_m-{\widehat{y}}_m^n}{\sigma }\right]}^2\right)$$

(14)

where ${\sigma }^2$ is the error variance of the measurements. The error variance ${\sigma }^2$, incorporates error in the measurement instrument and also the error associated with comparing the model predictions with sensor measurement with a different spatial and temporal average. The Gaussian likelihood function is inappropriate when the error in the data is correlated. In real time operation of the sensor network, the sensors might have a bias of underestimating or overestimating the concentration. Methods for tackling such situations are discussed in the literature [47]. For our illustration of the approach, we assume the measurement errors are uncorrelated and can be described by the Gaussian function. The calculation for the function parameter $\sigma$ can be done in two ways, first, the standard deviation may be known a priori based upon previous statistical analysis of field sampling and/or laboratory measurement error. Second, the standard deviation of the data error can be estimated using maximum likelihood theory [39]. Here, $\sigma$ is assigned by the standard deviation of synthetic measurement at the corresponding sensor. $\sigma$ both influential and difficult to be determined for the estimation. A sensitivity study can be performed for analyzing the influence of $\sigma$ for the problem setup. The real-time application of the full approach is shown as a flow diagram in Figure 2.

3. Results

3.1. Validation of Contaminant Transport from PF Operator

The method used for the construction of the PF operator (Markov matrix) described in Section 2.2 is validated for an IEA-Annex building geometry [48] with a manikin. The geometry, boundary conditions, and the obtained flow field from the CFD simulation are shown in Figure 3a,b. The CFD computations are carried using the open-source CFD tool OpenFOAM [49] We use parallel computing to overcome the burden of numerically solving PDE for constructing the PF operator. The PF operator constructed from the flow field is then used to transport the contaminant $\mathsf{\Phi }$. The comparison of this approach is made against the numerical solution of the PDE transport Equation (1a). The passive scalar is normalized between (0,1) and initialized as one covering half of the domain. The comparison after evolving the system to the final time of 50 s is shown in Figure 4b,c. It can be observed that the contaminant evolution contours computed by simple matrix-vector product (where the matrix is the PF operator acting on the contaminant vector in Equation 5) are indistinguishable from the numerical solution.

Further, a comparison of the concentration profile along the midplane axes of the building is shown in Figure 4d. The dots in the figure show the Markov results, which closely match and overlap the PDE-based predictions. Similar validation has been carried for the 3D case as well. The validations show the accuracy of the Markov approach and effectiveness of the matrix-vector product-based approach for transporting the contaminant concentrations. Throughout the study, we have not used any actual contaminant concentration measurements. The contaminant concentrations are simulated using CFD normalized between 0 to 1.

3.2. Contaminant Source Identification in 2D Office Space

We first illustrate the approach for a 2D office building as shown in Figure 5a for identifying an unknown source release location. The room is of 9 m × 3 m and is a steadily vented room. It has an occupant producing an average flux of 70 W, a computer producing 200 W, an adiabatic desk, and a window with the incoming heat flux of 100 W. With these as the boundary condition we compute the flow field using the bouyantSimpleFoam OpenFOAM CFD solver (validation of the solver shown in Appendix A). The converged mesh used in the study is shown in Figure 5b, and the flow field in Figure 5c. The transfer operator $\mathbf{P}\in {\mathbf{R}}^{9400\times 9400}$ is constructed using the flow field. The constructed operator is used to obtain the optimal locations of eight sensors using the algorithm shown in [34], with no placement constraints. The transfer operator associated time-step is ${\tau }_{mod}=10$ s. The observation by the sensor are made every ${\tau }_{obs}=50$ s, with measurement error variance $\sigma =0.03$. The final time for the Bayesian inference is set to $t_f=500$ s. This is decided based on the room air change rate, as it is important to detect the release location before the air is changed. In the domain we sample 50 quasi-random samples shown in Figure 5d these are used to generate model measurements for placed twelve sensors. These measurements are used in the likelihood in the Bayesian inference. To test the implementation we randomly pick three unsampled release locations to generate real-time measurement data with an added white noise with mean zero and variance 0.01. For these three random releases, we show plot the posteriors contours in the domain shown in Figure 6. It can be seen that the probability of release rapidly converges to the actual release region within one or two time-steps.

3.3. Contaminant Source Identification in 3D Furnished Room

The method is next illustrated for a case with a three-dimensional furnished room. The number of release states is nearly ten times that of the two-dimensional problem. The illustration shows the application of the methodology in a realistic indoor environment and the easy extension of the approach to a relatively complex geometry. Figure 7a shows the three-dimensional furnished room with a bed, TV-set, lamp, lighting, and cabinet represented in the primitive shapes. The room has dimensions 9 m × 3 m × 3 m size with an inlet and outflow to produce a steady vented room. A window is positioned on the right wall with a constant temperature of 302 K. These boundary conditions are used to solve the flow field using the OpenFOAM CFD solver buoyantBoussinesqSimpleFoam. The mesh convergence analysis is carried to ensure the grid independence of the results and 0.6 M cells are used for CFD computation (validation of the used solver and grid independence shown in Appendix A and Appendix B). As described in Section 2.2 the size of the transfer operator is equivalent to the number of states/cells used in the computation of the flow field. We note that the CFD mesh needs to be well resolved to capture the flow features, but the contaminant transport can be reliably simulated on a coarser mesh. We, therefore, interpolate the flow field onto a coarser mesh size of around 70 K cells. The use of CFD mesh would result in a transfer matrix of size $0.6\mathrm{M}\times 0.6\mathrm{M}$. This large matrix will require a lot of memory to store, and subsequent analysis will become infeasible. Therefore, the transfer operator constructed is of size $\mathbb{P}\in {\mathbf{R}}^{70\mathrm{K}\times 70\mathrm{K}}$. The optimal sensor algorithm shown in [34] is used to place twelve sensors in the room.

The transfer operator has the evolution time step ${\tau }_{mod}=10$ s. The sensors samples every ${\tau }_{obs}=150$ s, the final simulation time $t_f=720$ s based on the air change rate of the room. The room is sampled with 50 quasi-random samples in the domain with an epsilon radius of 0.5 m are shown in, Figure 7d, as release locations. For these samples, the model measurement data are computed used for likelihood calculation during real-time Bayesian inference. Three separate random un-sampled release location. We test more than three unknown random release locations. For presentation purposes, we only show three out of them. is simulated to demonstrate the real-time implementation of the approach. The real-time measurements of the sensors are generated with an additional white noise of zero mean and variance of 0.03. Figure 8 shows the posterior map of the release location for the three releases at every sampling time-step. The posteriors are converging to the actual release since we have sufficiently sampled the space to construct the model measurement data. In addition, the result is a posterior distribution of the probable release location. The actual epsilon ball of the release location is not recovered. It can be seen that the posterior closely represent the actual release location in the domain as probable regions.

4. Discussion

The sequential Bayesian inference method in association with the transfer operator approach for simulating the contaminant evolution is developed. The contaminant release locations are identified as a posterior map in the domain. The approach overcomes the real-time implementation limitation posed by traditional Markov Chain Monte Carlo (MCMC) sampling by relying on a Bayes Monte Carlo approach with pre-sampling from the prior distribution of the release locations. The sampling is exhaustive to ensure that the integration in the Bayes formula is accurate. We simulate contaminant releases from these sampled locations and store sensor time-series data. These actual measurements are used for comparing against the real-time streaming data from sensors using the likelihood function. The actual measurements can be generated concurrently as well by simply evolving each release location till the real-time measurement step. The approach is illustrated for both two and three-dimensional building space.

Earlier works in the literature [24,47,50] have reported the use of similar Bayes formulations for source identification but were either based on zonal models or empirical correlations as the forward models. In addition, some reported works rely on steady-state assumptions for estimating contaminant distribution. These assumptions are not well suited for studying the dynamic nature of contaminant evolution in the domain. Therefore, these approaches are usually sub-optimal for real-time implementation and accurately locating release locations. In the present approach, we try to overcome the limitation of the computational complexity of the contaminant transport model via a PF operator approach that results in a simple matrix-vector multiplication for contaminant evolution. The dynamical system setting used here has previously been used for applications including sensor placement under deterministic or stochastic setting [10,35,51,52]. This PF approach provides a unified framework for analysis of building systems during the design/planning phase as well as for risk assessment.

One of the limitations of the present approach is that a good estimate for $\sigma$ for the likelihood function is needed. However, the choice of sigma can be improved using prior model measurements under hypothetical release cases, to reduce the error in the prediction. Further, with the use of good computational methods like adaptive sampling, approximate Bayesian computation using sequential Monte Carlo can be performed.

The accidental or intentional release of hazardous compounds (even in minuscule quantities for short periods) can result in lethal and catastrophic consequences. Therefore, a fast, accurate and reliable source identification method is important for containment and speedy evacuation of the occupants. It is important to resolve the smallest to the largest time scales of the contaminant evolution which can be accomplished easily using the PF operator approach. Furthermore, these estimated contaminant locations can easily assist the rescue team in deciding the efficient evacuation strategy and can also be used to controlling the HVAC unit for the building. In addition, for complex imperfect dynamical systems described in [53], the approach could also be leveraged for designing robust control.

5. Conclusions

We present a method for identifying contaminant source location using the PF operator-based transport model coupled with a sequential Bayesian inference formulation. The approach PF based contaminant transport model provides a fast, accurate, and robust method to act as a forward model for contaminant evolution. The approach reduces the computation cost of numerically solving a PDE to simple matrix-vector multiplication. This makes Bayesian inference computationally easy to implement in real-time. The data required for computing the likelihood can be pre-computed. This means that if we have past measurements and some intuition of the release location, pre-computation is possible or the calculations can be made on the fly for multiple sampled release scenario. In the current work, we have shown the use of a quasi-random space-filling sampling approach to sampling the probable release locations in the room and have conducted the likelihood estimate on the fly . For these locations, the constant contaminant release was simulated.

The temporal measurements of the concentrations at the placed sensors in the room the temporal were stored for computing the likelihood in real-time operation. The implementation results are illustrated both in two-dimensional and three-dimensional building problems. We show the results of source identification both in 2D and a 3D problem for three unknown random constant contaminant release scenarios and compute the temporal updated posterior of release location. The results show that the posterior converges to a region within two-three time steps. In the future, we plan to extend the approach for computing various other parameters such as, release rate and implementing the approach for multiple release identification.