Method to Locate Contaminant Source and Estimate Emission Strength

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.


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][4][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][9][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.

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: where  is the mass pollutant concentration, mg/m 3 ; u 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: where  +1 and   are the column vectors of the contaminant concentration for each grid node at  + 1 and , 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: 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 where ) ,  = (  0 0)  . (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,  1 , to overcome the system uncertainty.Thus, the model turns into a new discretized stochastic model: where   represents the stochastic disturbance transition matrix;  1 is an uncorrelated Gaussian white sequence with a zero mean and [ 1 ⋅   1 ] = ;  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 .

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.

Preliminary Identification of Source Position and Emission
Strength.Assume that an unknown source releases CO 2 at a constant value in a short time, and we define a sensitivity coefficient as follows [16]: where 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 Since   is independent of the coordinate axes, we can rewrite (9) as Replacing /  by   , we get the direct well-posed equation as follows: 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 : 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 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  = ( 1 , . . .,   , . . .,   ) represent the predicted error at time , and then So source strength can be written as Δ  =   /   .Redefine the source strength as 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].

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  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.

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 where   represents the source strength  is the source position node.Assume that the source term is relatively unchanging in short time, and it satisfies where  is the uncertainty associated with the unknown source.Integrating and discretizing (18), we will obtain where   is an  ×  matrix.
The state equation including the unknown source is where is an × diagonal matrix and has the following structure: 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: where V +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  1 = [ 11 0], where  11 satisfies  11  = ; then The predicted estimation of  =  1  is given by The predicted error covariance matrix is given by The Kalman gain is given by has the structure  +1 = ( 1  ,  2  )  .
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 The actual estimation error covariance matrix  +1|+1 is then calculated from the predicted error covariance matrix by using From ( 27)-( 28), we can realize the identification of emission strength as well as the prediction of contaminant concentration distribution in the confined space.

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 2 concentration in the cabin is 0.5 g/m 3 .The CO 2 concentration in the inlet stream is also 0.5 g/m 3 .A sudden source appears at point A, and its emission strength versus time is Four sensors with 1% measurement noise are placed at B 1 , B 2 , B 3 , and B 4 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.
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.
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 3 .Therefore, the presented algorithm can predict the concentration in the enclosed space very well.

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.

4 Figure 2 :
Figure 2: Air flow, source, and sensors in the cabin.

Table 1 :
Emission strengths estimated by sensors.