Method to Locate Contaminant Source and Estimate Emission Strength

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.

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. Section 2 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, 𝐮, 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: 𝜕𝐶()𝜕𝑡+∇𝐮𝐶=∇𝑑⋅grad𝐶+𝑆, (1) where 𝐶 is the mass pollutant concentration, mg/m³; 𝐮 is the velocity, m/s; 𝑑 is the diffusion coefficient, mg/(m·s); 𝑆 is the contaminant source; 𝑡 is the forward time, 𝑠.

Equation (1) can be solved numerically. When the contaminant concentration distribution at time 𝑡 is known, the sequent at time 𝑡+1 can be obtained: 𝐴₁𝐶^𝑡+1=𝐴₂𝐶^𝑡+𝑆, (2) where 𝐶^𝑡+1 and 𝐶^𝑡 are the column vectors of the contaminant concentration for each grid node at 𝑡+1 and 𝑡, respectively. 𝐀 is a diagonally dominant banded sparse matrix; hence it must have the inverse matrix 𝐀⁻¹.

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: −𝐴_𝑥2+𝑎⁰_𝑃𝐶^∗=𝐴_𝑥2+𝐴_𝑦+𝐴_𝑧+𝑎⁰_𝑝𝐶^𝑡−𝐴+𝑏,_𝑦2+𝑎⁰_𝑃𝐶^∗∗𝐴=−^𝑦2𝐶^𝑡+𝑎⁰_𝑃𝐶^∗,−𝐴_𝑧2+𝑎⁰_𝑃𝐶^𝑡+1𝐴=−^𝑧2𝐶^𝑡+𝑎⁰_𝑃𝐶^∗∗, (3) where 𝐴_𝑥, 𝐴_𝑦, and 𝐴_𝑧 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 𝑀𝐂^𝑡+1=𝑁𝐂^𝑡+𝐵, (4) where 𝐂^𝑡=𝐶^∗_𝑡^𝑇,𝐶_𝑡^𝑇∗∗,𝐶^𝑇_𝑡^𝑇,−𝐴𝑀=^𝑥2+𝑎⁰_𝑝00−𝑎⁰_𝑝−𝐴_𝑦2+𝑎⁰_𝑝00−𝑎⁰_𝑝−𝐴_𝑧2+𝑎⁰_𝑝,𝐴𝑁=00^𝑥2+𝐴_𝑦+𝐴_𝑧+𝑎⁰_𝑝𝐴00−^𝑦2𝐴00−^𝑧2,𝐵=𝑏^𝑇00^𝑇. (5)

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, 𝜔₁, to overcome the system uncertainty. Thus, the model turns into a new discretized stochastic model: 𝑀𝐶^𝑡+1=𝑁𝐶^𝑡+𝐵+𝑈_𝑡𝜔₁, (6) where 𝑈_𝑡 represents the stochastic disturbance transition matrix; 𝜔₁ is an uncorrelated Gaussian white sequence with a zero mean and 𝐸[𝜔₁⋅𝜔^𝑇₁]=𝑄; 𝑄 is a diagonal matrix that represents the model noise.

A highly accurate model would have a low value of 𝑄; on the other side, a model that does not reflect the physical process too accurately would have a high value of 𝑄.

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 CO₂ at a constant value in a short time, and we define a sensitivity coefficient as follows [16]: 𝑍=𝑍^𝑇₁,𝑍^𝑇₂,…,𝑍^𝑇_𝑁^𝑇, (7) where 𝑍_𝑖=𝜕𝐶𝜕𝑆_𝑢𝑖,𝑖=1,2,…,𝑁. (8)

𝑆_𝑢𝑖 represents the assumed unknown source at point 𝑖; the sensitivity coefficient 𝑍_𝑖 is a 1×𝑁 vector and represents how the source at the 𝑖th grid point affects the concentration distribution. Here 𝑁 is the number of grids.

In order to calculate 𝑍, we multiply (1) by 𝜕/𝜕𝑆_𝑢𝑖 and then get 𝜕𝜕𝑆_𝑢𝑖𝜕𝐶()𝐶𝜕𝑡+∇𝐮𝐶=∇Γ⋅grad𝜌+𝑆_𝑢. (9)

Since 𝑆_𝑢 is independent of the coordinate axes, we can rewrite (9) as 𝜕𝜕𝑡𝜕𝐶𝜕𝑆_𝑢𝑖+∇𝐮𝜕𝐶𝜕𝑆_𝑢𝑖Γ=∇𝜌⋅grad𝜕𝐶𝜕𝑆_𝑢𝑖+𝜕𝑆_𝑢𝜕𝑆_𝑢𝑖. (10)

Replacing 𝜕𝐶/𝜕𝑆_𝑢𝑖 by 𝑍_𝑖, we get the direct well-posed equation as follows: 𝜕𝑍_𝑖𝜕𝑡+∇𝐮⋅𝑍_𝑖Γ=∇𝜌grad𝑍_𝑖+𝐼_𝑖, (11) where 𝐼_𝑖 represents the 𝑖th column of the identity matrix.

The boundary and initial conditions will also be divided by 𝜕/𝜕𝑆_𝑢𝑖 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), 𝑍(𝑡), 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, 𝑍(𝑡) can be divided into a usable part, 𝑍_sensors, and a nonusable part, 𝑍_non-sensors: 𝑍=𝑍_sensors,𝑍_non-sensors. (12)

Here we define 𝑡_emit and 𝜏_𝑤. 𝑡_emit represents the time when the sudden source appears, and 𝜏_𝑤 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 𝑍^𝑐_𝑖𝑗=𝑍_{sensors∣𝜏𝑤=𝑡−𝑡emit}, (13) where 𝑍^𝑐_𝑖𝑗 is an 𝑁×𝑚 matrix; 𝑚 is the number of sensors.

𝑍^𝑐_𝑖𝑗 represents the concentration change at the 𝑗th sensor position when the source changes one unit strength at the 𝑖th grid point.

After a sudden unknown source appears, the estimated value will not be consistent with the measurement value. Here we let 𝑒=(𝑒₁,…,𝑒_𝑖,…,𝑒_𝑚) represent the predicted error at time 𝑡, and then 𝑍^𝑐_𝑖𝑗=𝜕𝐶_𝑗𝜕𝑆_𝑢𝑖=𝑒_𝑗Δ𝑆_𝑢𝑖. (14)

So source strength can be written as Δ𝑆_𝑢𝑖=𝑒_𝑗/𝑍^𝑐_𝑖𝑗. Redefine the source strength as STR_𝑖𝑗=𝑒_𝑗𝑍^𝑐_𝑖𝑗, (15) where STR is an 𝑁×𝑚 matrix and STR_𝑖𝑗 represents the source strength estimated by the 𝑗th sensor when the source is assumed at 𝑖th grid point.

Then the estimation algorithm is given as follows:(1)predict the error 𝑒 by measured data of sensors at 𝑡=𝑡_emit+𝜏_𝑤;(2)input the sensitivity matrix 𝑍^𝑐_𝑖𝑗(𝜏_𝑤);(3)calculate the source strength matrix STR by using (15);(4)scale each 𝑚-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 each 𝑚-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 time 𝑡=𝑡_emit+𝜏_𝑤. 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 𝑚 or less than 𝑚.

The modified algorithm is shown in Figure 1. 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 the 𝑗th sensor is 𝛿𝑒_𝑗, which will lead to an estimated error of source strength 𝛿𝑒_𝑗/𝑍^𝑐_𝑖𝑗 based on (15). In order to get an identifiable zone of 𝑗th sensor with a fiducial probability 𝑝, the estimated error of emission strength must be no more than 1−𝑝. The identifiable zone with a fiducial probability 𝑝 of a single sensor is defined as 𝑍^𝑐_𝑖𝑗>𝛿𝑒_𝑗.1−𝑝 (16)

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 dimension 𝑁 is given by 𝑆_𝑢=0,…,0,𝑠_𝑢𝑖,0^𝑇, (17) where 𝑠_𝑢𝑖 represents the source strength 𝑖 is the source position node.

Assume that the source term is relatively unchanging in short time, and it satisfies 𝑑𝑆_𝑢𝑑𝑡=0+𝜔, (18) where 𝜔 is the uncertainty associated with the unknown source. Integrating and discretizing (18), we will obtain 𝑆_𝑢,𝑡+1=𝑆_𝑢,𝑡+𝐾_𝑡𝜔, (19) where 𝐾_𝑡 is an 𝑁×𝑁 matrix.

The state equation including the unknown source is 𝑀^𝑎𝐂^𝑎_𝑡+1=𝑁^𝑎𝐂^𝑎_𝑡+𝐵0+𝜔_1𝑡𝜔_𝑡, (20) where 𝑀^𝑎=^𝑀0_0𝑊 and 𝑁^𝑎=^𝑁𝑊_0𝑊. 𝑊 is an 𝑁×𝑁 diagonal matrix and has the following structure: 0⋱01⋱0𝑊=0𝑂𝑂0. (21)

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: 𝑧_𝑡+1=[𝐻0]𝐂^𝑎_𝑡+1+𝑣_𝑡+1, (22) where 𝑣_𝑡+1 is an uncorrelated Gaussian white noise sequence that represents the measurement noise with covariance represented by 𝑅.

The matrix, 𝑅, is a diagonal matrix that contains information about the uncertainty of the measurement. The lower values for 𝑅 indicate that the sensors have more accurate for measurement.

Define 𝐻₁=[𝐻₁₁0], where 𝐻₁₁ satisfies 𝐻₁₁𝑀=𝐻; then 𝑧_𝑡+1=𝐻₁₁0𝑀𝐶^𝑎_𝑡+1+𝑣_𝑡+1. (23)

The predicted estimation of 𝑌=𝐴₁𝐶 is given by 𝑌_{𝑡+1∣𝑡}=𝑁𝐂_𝑡∣𝑡+𝐵_𝑡+𝑆_𝑢,𝑡,𝑆_{𝑢,𝑡+1∣𝑡}=𝑆_𝑢,𝑡. (24)

The predicted error covariance matrix is given by 𝑃_{𝑡+1∣𝑡}=𝑁𝑃𝑁^𝑇+𝑃_𝑠+𝑄𝑊𝑃_𝑠𝑊^𝑇+𝑄𝑊𝑃_𝑠𝑊^𝑇+𝑄𝑊𝑃_𝑠𝑊^𝑇+𝑄. (25)

The Kalman gain is given by 𝐾_𝑡+1=𝑃_{𝑡+1∣𝑡}𝐻^𝑇₁𝐻₁𝑃_{𝑡+1∣𝑡}𝐻^𝑇₁+𝑅⁻¹. (26)

𝐾 has the structure 𝐾_𝑡+1=(𝐾₁^𝑇,𝐾₂^𝑇)^𝑇.

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 𝑀𝐶_𝑡+1=𝑌_{𝑡+1∣𝑡}+𝐾₁𝑧_𝑡+1−𝐻₁₁𝑌_{𝑡+1∣𝑡},𝑆implicit_𝑡+1=𝑆_{𝑡+1∣𝑡}+𝐾₂𝑧_𝑡+1−𝐻₁₁𝑌_{𝑡+1∣𝑡}explicit. (27)

The actual estimation error covariance matrix 𝑃_{𝑡+1∣𝑡+1} is then calculated from the predicted error covariance matrix by using 𝑃_{𝑡+1∣𝑡+1}=𝐼−𝐾𝐻₁𝑃𝐼−𝐾𝐻₁^𝑇+𝐾𝑅𝐾^𝑇. (28)

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 Figure 2. The air flow is obviously divided into two zones. The initial CO₂ concentration in the cabin is 0.5 g/m³. The CO₂ concentration in the inlet stream is also 0.5 g/m³. A sudden source appears at point A, and its emission strength versus time is 𝑆=()15,𝑡≤2015−0.2𝑡−20,20<𝑡≤607,𝑡>60. (29)

Figure 2 

Air flow, source, and sensors in the cabin.

Four sensors with 1% measurement noise are placed at B₁, B₂, B₃, and B₄ 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 point A 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.

Inputting 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 Figures 3–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 at 𝑡=120 s.

The real line in Figure 3 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 Figure 5 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.

Figures 6, 7, and 8 show the actual concentration distribution, the estimated concentration distribution, and the estimated error in the cabin at time 𝑡=120 s, respectively. From Figure 8, we know that the maximum error is less than 0.05 g/m³. Therefore, the presented algorithm can predict the concentration in the enclosed space very well.

Figure 7 

Estimated concentration distribution at 𝑡=120 s by Kalman filter.

Figure 8 

Concentration estimated error at 𝑡=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.