# Method to Locate Contaminant Source and Estimate Emission Strength

Authors: Qu Hongquan; Wang Yongjiao; Pang Liping; Chang Haijuan; Zhang Changnian
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

People greatly concern the issue of air quality in some confined spaces, such as spacecraft, aircraft, and submarine. With the increase of residence time in such confined space, contaminant pollution has become a main factor which endangers life. It is urgent to identify a contaminant source rapidly so that a prompt remedial action can be taken. A procedure of source identification should be able to locate the position and to estimate the emission strength of the contaminant source. In this paper, an identification method was developed to realize these two aims. This method was developed based on a discrete concentration stochastic model. With this model, a sensitivity analysis algorithm was induced to locate the source position, and a Kalman filter was used to further estimate the contaminant emission strength. This method could track and predict the source strength dynamically. Meanwhile, it can predict the distribution of contaminant concentration. Simulation results have shown the virtues of the method.


---

## Body

## 1. Introduction

A pollution incident has the potential harm on people, especially when it happens in an enclosed ventilation space, such as a cabin in a manned spacecraft, submarine, or aircraft [1, 2]. It is necessary to develop related technology to efficiently identify a sudden contaminant source so that control actions can be taken rapidly. Contaminant source identification is a process that realizes source term estimation reversely by using limited information. However, this inverse identification for a source is extremely difficult. In the recent years, researchers from different countries carried on some related studies on the source identification. Most of the methods used retrospective methods. They used the sensor observed data to estimate the unknown contaminant source and then reconstruct the polluting progress. Methods widely used include the analytical approach, the optimization approach, the probabilistic approach, and the backward computational fluid dynamics approach. Skliar, Ramirez, Mamonov and Tsai carried out much research about point source identification in an enclosed environment [3–5]. Liu, Zhai, Zhang, and Chen thoroughly researched to identify virus source in an aircraft cabin [6, 7]. Wagner, Mahar and Datta proposed an optimization method to identifying the source term with combined forward simulation, optimization and regression methods [8–10]. Neupauer, Wilson and Sohn et al. developed backward location and travel time possibilities to realize point source estimation [11, 12]. The above endeavors have been promoting the development of source identification.In some actual instance, the emission process of contaminant sometimes is a continuous and not instantaneous process, and the measured data is influenced by sensor noise. Considering these two factors, a dynamic method for identifying unknown contaminant source will be discussed in this paper. The rest of this paper is organized as follows. Section2 introduces a discretized concentration stochastic model. Section 3 introduces a sensitivity analysis algorithm for preliminary analysis on the strength and location of an unknown source. A dynamic method by combining implicit and explicit Kalman filters is presented in Section 4 to identify the source strength characteristic and predict the distribution of concentration dynamically. Section 5 contains numerical simulation to demonstrate the method.


## 2. Discretized Concentration Stochastic Model

We will set up a three-dimensional discretized concentration transport model by discretizing the convection-diffusion transport model. In this model, we introduce a stochastic noise term to overcome modeling error caused by fluctuation of air flow or some other factors.For an incompressible fluid, the distribution of air flow velocity,u, in enclosed space can be obtained by using SIMPLE algorithm [13, 14]. Assume that a source emits contaminant in a constant temperature condition and the viscous stresses are very small; then a gaseous contaminant transport process by air is governed by the following equation:
(1)∂C∂t+∇(uC)=∇(d·gradC)+S,
where C is the mass pollutant concentration, mg/m3; u is the velocity, m/s; d is the diffusion coefficient, mg/(m·s); S is the contaminant source; t is the forward time, s.Equation (1) can be solved numerically. When the contaminant concentration distribution at time t is known, the sequent at time t+1 can be obtained:
(2)A1Ct+1=A2Ct+S,
where Ct+1 and Ct are the column vectors of the contaminant concentration for each grid node at t+1 and t, respectively. A is a diagonally dominant banded sparse matrix; hence it must have the inverse matrix A-1.However, the calculation load is heavy to solve (2). For example, the equation set of 8000 dimensions will be solved for a mesh of 20 × 20 × 20. This will take enormous memory and time. The explicit methods do not suffer from this because of their poor stability properties. So the solution scheme used in this paper is the classic alternating direction implicit scheme (ADI) [15]. Then the discretized concentration equations of (2) can be expressed as the following equations:
(3)(-Ax2+aP0)C*=(Ax2+Ay+Az+ap0)Ct+b,(-Ay2+aP0)C**=-Ay2Ct+aP0C*,(-Az2+aP0)Ct+1=-Az2Ct+aP0C**,
where Ax, Ay, and Az are tridiagonal matrices. The equations set in (3) can be represented by a single matrix equation in terms of the state transition matrices, as shown in
(4)MCt+1=NCt+B,
where
(5)Ct=(Ct*T,Ct**T,CtT)T,M=(-Ax2+ap000-ap0-Ay2+ap000-ap0-Az2+ap0),N=(00Ax2+Ay+Az+ap000-Ay200-Az2),B=(bT00)T.Considering the system uncertainty caused by some rounding errors caused by discretizing, fluctuation of air velocity, eddy diffusivity, and so on, we add a stochastic disturbance,ω1, to overcome the system uncertainty. Thus, the model turns into a new discretized stochastic model:
(6)MCt+1=NCt+B+Utω1,
where Ut represents the stochastic disturbance transition matrix; ω1 is an uncorrelated Gaussian white sequence with a zero mean and E[ω1·ω1T]=Q; Q is a diagonal matrix that represents the model noise.A highly accurate model would have a low value ofQ; on the other side, a model that does not reflect the physical process too accurately would have a high value of Q.


## 3. Identification Method Based on the Sensitivity Analysis Algorithm

An identification method for sudden contaminant source should be able to realize two functions: position location and emission strength estimation. The latter cannot be realized without accurate position information of source. Here we locate a sudden source by using the sensitivity analysis algorithm [16]. This algorithm can give an initial guess of emission strength for an unknown source. Then extended Kalman filter (EKF) algorithm will be used to further track the characters of dynamic emission strength.


### 3.1. Preliminary Identification of Source Position and Emission Strength

Assume that an unknown source releases CO2 at a constant value in a short time, and we define a sensitivity coefficient as follows [16]:
(7)Z=(Z1T,Z2T,…,ZNT)T,
where
(8)Zi=∂C∂Sui,i=1,2,…,N.S
u
i represents the assumed unknown source at point i; the sensitivity coefficient Zi is a 1×N vector and represents how the source at the ith grid point affects the concentration distribution. Here N is the number of grids.In order to calculateZ, we multiply (1) by ∂/∂Sui and then get
(9)∂∂Sui[∂C∂t+∇(uC)=∇[Γ·grad(Cρ)]+Su].SinceSu is independent of the coordinate axes, we can rewrite (9) as
(10)∂∂t∂C∂Sui+∇(u∂C∂Sui)=∇[Γρ·grad(∂C∂Sui)]+∂Su∂Sui.Replacing∂C/∂Sui by Zi, we get the direct well-posed equation as follows:
(11)∂Zi∂t+∇(u·Zi)=∇[(Γρ)gradZi]+Ii,
where Ii represents the ith column of the identity matrix.The boundary and initial conditions will also be divided by∂/∂Sui to obtain the initial and boundary conditions for sensitivity problem.Because the equation structure is unchanged from the original model equation, the same algorithm and numerical techniques can be used to calculate the sensitivity matrix. In order to save the calculation time, it is worth attention that this computation part can be performed offline before the real contaminant source appears. This will save much time for the source identification.The solution to (11), Z(t), is a function of time. Since we only have the concentration information at sensor positions, we are interested in how the source at different position will influence the concentration at sensor position. Therefore, Z(t) can be divided into a usable part, Zsensors, and a nonusable part, Znon-sensors:
(12)Z=[Zsensors,Znon-sensors].Here we definetemit and τw. temit represents the time when the sudden source appears, and τw is the lag time when a source is detected by a sensor. In order to identify the source characters more accurately, we identify the source by using
(13)Zijc=Zsensors∣τw=t-temit,
where Zijc is an N×m matrix; m is the number of sensors.Z
i
j
c represents the concentration change at the jth sensor position when the source changes one unit strength at the ith grid point.After a sudden unknown source appears, the estimated value will not be consistent with the measurement value. Here we lete=(e1,…,ei,…,em) represent the predicted error at time t, and then
(14)Zijc=∂Cj∂Sui=ejΔSui.So source strength can be written asΔSui=ej/Zijc. Redefine the source strength as
(15)STRij=ejZijc,
where STR is an N×m matrix and STRij represents the source strength estimated by the jth sensor when the source is assumed at ith grid point.Then the estimation algorithm is given as follows:(1)
predict the errore by measured data of sensors at t=temit+τw;(2)
input the sensitivity matrixZijc(τw);(3)
calculate the source strength matrixSTR by using (15);(4)
scale eachm-dimension row vector using the mean of the strength values within the row vector;(5)
calculate the standard deviation of the calculated strength values for eachm-dimension row vector;(6)
choose the row with the minimum standard deviation, and this row is the position of source;(7)
calculate the mean of this row vector and obtain the emission strength.So the identification algorithm based on sensitivity matrix is developed, and we can call it the sensitivity analysis algorithm. Actually, there may be only a part of the sensor that predicts errors exceeding the threshold value at timet=temit+τw. In other words, the sudden unknown source only affects a part of concentration values. Some sensors can emit the concentration change and others cannot. Hence, these sensors are classified as activated sensors andinactivated sensors. In this case, the inactivated sensors need to be eliminated from our analysis. Otherwise, it will affect the accuracy of the sensitivity analysis algorithm result. Therefore, the row dimension of STR matrix we used actually is m or less than m.The modified algorithm is shown in Figure1. After the position and the initial guess strength of unknown source are preliminarily identified by using the sensitivity analysis algorithm, the Kalman filter will be used to track the change of strength versus time dynamically [17].Figure 1
Sensitivity analysis algorithm.


### 3.2. Identifiable Zone

Another issue that should be discussed in the analysis of source location is the identifiable zone. The sensitivity analysis algorithm can identify the source position only when two sensors are activated. Therefore, identifiable zone is the positions where two activated sensors both emit the sudden source appearing. Assume that the measurement error of thejth sensor is δej, which will lead to an estimated error of source strength δej/Zijc based on (15). In order to get an identifiable zone of jth sensor with a fiducial probability p, the estimated error of emission strength must be no more than 1-p. The identifiable zone with a fiducial probability p of a single sensor is defined as
(16)Zijc>δej1-p.The zone where no less than two sensors are satisfied (16) is called identifiable zone.The sensitivity matrix contains useful information about the speed at which a perturbation travels through the system. Meanwhile, it can be used to analyze the identifiable zones of the sensitivity analysis algorithm.There are four main factors that affect the identifiable zone, and they are the sensor placement, the sensor measurement noise, the fiducial probability, and the latency time. We analyzed the change of identifiable zone by using different above parameters. The following conclusions can be drawn.(1)
Different sensor placement leads to different identifiable zone: the identification effect in the case where sensors are placed collectively at the downstream of the flow is better than that in the case where that sensors are placed dispersedly.(2)
Different fiducial probability leads to different identifiable zone: the lower the fiducial probability is, the more extensive the identifiable zone is. So the fiducial probability must be fixed according to the practical situation.(3)
Different latency time leads to different identifiable zone: the longer the latency time is, the bigger the identifiable zone is.(4)
Different measurement noise leads to different identifiable zone: a larger measurement noise will lead to a smaller identifiable zone. Therefore, the measurement noise must be considered when we choose the sensors.


## 3.1. Preliminary Identification of Source Position and Emission Strength

Assume that an unknown source releases CO2 at a constant value in a short time, and we define a sensitivity coefficient as follows [16]:
(7)Z=(Z1T,Z2T,…,ZNT)T,
where
(8)Zi=∂C∂Sui,i=1,2,…,N.S
u
i represents the assumed unknown source at point i; the sensitivity coefficient Zi is a 1×N vector and represents how the source at the ith grid point affects the concentration distribution. Here N is the number of grids.In order to calculateZ, we multiply (1) by ∂/∂Sui and then get
(9)∂∂Sui[∂C∂t+∇(uC)=∇[Γ·grad(Cρ)]+Su].SinceSu is independent of the coordinate axes, we can rewrite (9) as
(10)∂∂t∂C∂Sui+∇(u∂C∂Sui)=∇[Γρ·grad(∂C∂Sui)]+∂Su∂Sui.Replacing∂C/∂Sui by Zi, we get the direct well-posed equation as follows:
(11)∂Zi∂t+∇(u·Zi)=∇[(Γρ)gradZi]+Ii,
where Ii represents the ith column of the identity matrix.The boundary and initial conditions will also be divided by∂/∂Sui to obtain the initial and boundary conditions for sensitivity problem.Because the equation structure is unchanged from the original model equation, the same algorithm and numerical techniques can be used to calculate the sensitivity matrix. In order to save the calculation time, it is worth attention that this computation part can be performed offline before the real contaminant source appears. This will save much time for the source identification.The solution to (11), Z(t), is a function of time. Since we only have the concentration information at sensor positions, we are interested in how the source at different position will influence the concentration at sensor position. Therefore, Z(t) can be divided into a usable part, Zsensors, and a nonusable part, Znon-sensors:
(12)Z=[Zsensors,Znon-sensors].Here we definetemit and τw. temit represents the time when the sudden source appears, and τw is the lag time when a source is detected by a sensor. In order to identify the source characters more accurately, we identify the source by using
(13)Zijc=Zsensors∣τw=t-temit,
where Zijc is an N×m matrix; m is the number of sensors.Z
i
j
c represents the concentration change at the jth sensor position when the source changes one unit strength at the ith grid point.After a sudden unknown source appears, the estimated value will not be consistent with the measurement value. Here we lete=(e1,…,ei,…,em) represent the predicted error at time t, and then
(14)Zijc=∂Cj∂Sui=ejΔSui.So source strength can be written asΔSui=ej/Zijc. Redefine the source strength as
(15)STRij=ejZijc,
where STR is an N×m matrix and STRij represents the source strength estimated by the jth sensor when the source is assumed at ith grid point.Then the estimation algorithm is given as follows:(1)
predict the errore by measured data of sensors at t=temit+τw;(2)
input the sensitivity matrixZijc(τw);(3)
calculate the source strength matrixSTR by using (15);(4)
scale eachm-dimension row vector using the mean of the strength values within the row vector;(5)
calculate the standard deviation of the calculated strength values for eachm-dimension row vector;(6)
choose the row with the minimum standard deviation, and this row is the position of source;(7)
calculate the mean of this row vector and obtain the emission strength.So the identification algorithm based on sensitivity matrix is developed, and we can call it the sensitivity analysis algorithm. Actually, there may be only a part of the sensor that predicts errors exceeding the threshold value at timet=temit+τw. In other words, the sudden unknown source only affects a part of concentration values. Some sensors can emit the concentration change and others cannot. Hence, these sensors are classified as activated sensors andinactivated sensors. In this case, the inactivated sensors need to be eliminated from our analysis. Otherwise, it will affect the accuracy of the sensitivity analysis algorithm result. Therefore, the row dimension of STR matrix we used actually is m or less than m.The modified algorithm is shown in Figure1. After the position and the initial guess strength of unknown source are preliminarily identified by using the sensitivity analysis algorithm, the Kalman filter will be used to track the change of strength versus time dynamically [17].Figure 1
Sensitivity analysis algorithm.


## 3.2. Identifiable Zone

Another issue that should be discussed in the analysis of source location is the identifiable zone. The sensitivity analysis algorithm can identify the source position only when two sensors are activated. Therefore, identifiable zone is the positions where two activated sensors both emit the sudden source appearing. Assume that the measurement error of thejth sensor is δej, which will lead to an estimated error of source strength δej/Zijc based on (15). In order to get an identifiable zone of jth sensor with a fiducial probability p, the estimated error of emission strength must be no more than 1-p. The identifiable zone with a fiducial probability p of a single sensor is defined as
(16)Zijc>δej1-p.The zone where no less than two sensors are satisfied (16) is called identifiable zone.The sensitivity matrix contains useful information about the speed at which a perturbation travels through the system. Meanwhile, it can be used to analyze the identifiable zones of the sensitivity analysis algorithm.There are four main factors that affect the identifiable zone, and they are the sensor placement, the sensor measurement noise, the fiducial probability, and the latency time. We analyzed the change of identifiable zone by using different above parameters. The following conclusions can be drawn.(1)
Different sensor placement leads to different identifiable zone: the identification effect in the case where sensors are placed collectively at the downstream of the flow is better than that in the case where that sensors are placed dispersedly.(2)
Different fiducial probability leads to different identifiable zone: the lower the fiducial probability is, the more extensive the identifiable zone is. So the fiducial probability must be fixed according to the practical situation.(3)
Different latency time leads to different identifiable zone: the longer the latency time is, the bigger the identifiable zone is.(4)
Different measurement noise leads to different identifiable zone: a larger measurement noise will lead to a smaller identifiable zone. Therefore, the measurement noise must be considered when we choose the sensors.


## 4. Dynamic Identification of Source Strength

After a source position and its emission strength have been preliminary identified, we need to dynamically track the change of strength versus time by using Kalman filter.Compared with other static algorithms, Kalman filter can track the dynamic change of emission strength very well. So we use an implicit Kalman filter combined with an explicit Kalman filter to identify a dynamic source emission strength [17, 18].The unknown source vector of dimensionN is given by
(17)Su=(0,…,0,sui,0)T,
where sui represents the source strength i is the source position node.Assume that the source term is relatively unchanging in short time, and it satisfies(18)dSudt=0+ω,
where ω is the uncertainty associated with the unknown source. Integrating and discretizing (18), we will obtain
(19)Su,t+1=Su,t+Ktω,
where Kt is an N×N matrix.The state equation including the unknown source is(20)MaCt+1a=NaCta+(B0)+(ω1tωt),
where Ma=[M00W] and Na=[NW0W]. W is an N×N diagonal matrix and has the following structure:
(21)W=[00O⋱01O0⋱0].The matrix is extremely sparse, with all zeros except for source position, where there is a unit term. This unit term picks off the source strength scalar term on expansion to yield (21).The measurement matrix is as follows:(22)zt+1=[H0]Ct+1a+vt+1,
where vt+1 is an uncorrelated Gaussian white noise sequence that represents the measurement noise with covariance represented by R.The matrix,R, is a diagonal matrix that contains information about the uncertainty of the measurement. The lower values for R indicate that the sensors have more accurate for measurement.DefineH1=[H110], where H11 satisfies H11M=H; then
(23)zt+1=[H110]MCt+1a+vt+1.The predicted estimation ofY=A1C is given by
(24)Y^t+1∣t=NCt∣t+Bt+Su,t,Su,t+1∣t=Su,t.The predicted error covariance matrix is given by(25)Pt+1∣t=[NPNT+Ps+QWPsWT+QWPsWT+QWPsWT+Q].The Kalman gain is given by(26)Kt+1=Pt+1∣tH1T(H1Pt+1∣tH1T+R)-1.K has the structure Kt+1=(K1T,K2T)T.The estimation algorithm consists of an implicit and an explicit part. The implicit part is used in the estimation of the state, while a parallel explicit part is used in the estimation of the emission strength. The gain terms multiply the residual errors in order to update the state estimation for the next time step. Consider(27)MCt+1=Y^t+1∣t+K1(zt+1-H11Y^t+1∣t)(implicit),St+1=St+1∣t+K2(zt+1-H11Y^t+1∣t)(explicit).The actual estimation error covariance matrixPt+1∣t+1 is then calculated from the predicted error covariance matrix by using
(28)Pt+1∣t+1=(I-KH1)P(I-KH1)T+KRKT.From (27)-(28), we can realize the identification of emission strength as well as the prediction of contaminant concentration distribution in the confined space.


## 5. Simulations

Our simulations object is an aircraft cabin. By using SIMPLE algorithm, we get the air flow in the cabin as shown in Figure2. The air flow is obviously divided into two zones. The initial CO2 concentration in the cabin is 0.5 g/m3. The CO2 concentration in the inlet stream is also 0.5 g/m3. A sudden source appears at point A, and its emission strength versus time is
(29)S={15,t≤2015-0.2(t-20),20<t≤607,t>60.Figure 2
Air flow, source, and sensors in the cabin.Four sensors with 1% measurement noise are placed atB1, B2, B3, and B4 as shown in Figure 2.We use the presented identification method to locate the source and estimate its emission strength.The sensitivity analysis algorithm reports that pointA has the least standard deviation, 0.2004. The emission strengths estimated by each of the sensors are listed in Table 1.Table 1
Emission strengths estimated by sensors.
Sensor
B1
B2
B3
B4
Strengthsg/(m3·s)
14.7996
14.9639
15.2792
14.9670Inputting the sensitivity analysis algorithm results into the Kalman Filter, we can predict the dynamic change of source emission strength and contaminant concentration in the cabin, as shown in Figures3–6.Figure 3
Filter performance tracking for strength.Figure 4
RMSE of Kalman filter tracking for strength.Figure 5
Prediction of contaminant concentration at sensor 2.Figure 6
Concentration distribution att=120 s.The real line in Figure3 shows the actual source strength, while the dashed line shows the estimated source strength. Figure 4 shows the estimated error of source strength by Kalman filter. From these two figures, we can observe that our algorithm works well to estimate the source strength.The real line in Figure5 shows the actual contaminant concentration versus time measured by at sensor 2, while the dashed line shows the measurement value, and chalk line shows the filter estimated concentration.Figures6, 7, and 8 show the actual concentration distribution, the estimated concentration distribution, and the estimated error in the cabin at time t=120 s, respectively. From Figure 8, we know that the maximum error is less than 0.05 g/m3. Therefore, the presented algorithm can predict the concentration in the enclosed space very well.Figure 7
Estimated concentration distribution att=120 s by Kalman filter.Figure 8
Concentration estimated error att=120 s.


## 6. Conclusion

In order to identify a sudden contaminant source in an enclosed space, a source identification method is developed by using the sensitivity matrix and Kalman state equation. In this method, a sensitivity analysis algorithm is used to locate a source position and estimate its emission strength preliminarily. Implicit and explicit Kalman filters are used to trace and estimate the source emission strength dynamically and predict the concentration distribution together.Some simulations are conducted to investigate the effect of the identification method. Simulation results show that this presented method can identify an unknown source appearing in the enclosed space. The position and the emission strength of a sudden source can be obtained together as well as the concentration distribution prediction.



---

*Source: 163612-2013-12-17.xml*