Research Article  Open Access
Identifiability and Identification of Trace Continuous Pollutant Source
Abstract
Accidental pollution events often threaten people’s health and lives, and a pollutant source is very necessary so that prompt remedial actions can be taken. In this paper, a trace continuous pollutant source identification method is developed to identify a sudden continuous emission pollutant source in an enclosed space. The location probability model is set up firstly, and then the identification method is realized by searching a global optimal objective value of the location probability. In order to discuss the identifiability performance of the presented method, a conception of a synergy degree of velocity fields is presented in order to quantitatively analyze the impact of velocity field on the identification performance. Based on this conception, some simulation cases were conducted. The application conditions of this method are obtained according to the simulation studies. In order to verify the presented method, we designed an experiment and identified an unknown source appearing in the experimental space. The result showed that the method can identify a sudden trace continuous source when the studied situation satisfies the application conditions.
1. Introduction
A pollution incident has the potential to create significant harm on numbers of people, especially when it happens in an enclosed ventilation space such as a cabin in a manned spacecraft, submarine, or aircraft. It is urgent to develop related technology to efficiently locate a sudden pollutant source so that control actions can be taken rapidly.
Pollutant source identification is a process that searches a source term reversely by using limited information. It usually consists of detecting the initial emission time, estimating emission strength and finding the source location. However, inverse identification for a source is an extremely difficult challenge in the research field of inverse problems, namely, for the complexity in the coupling of source position, strength and emission time, and the incompleteness and uncertainty of observed data. With the improvement of computer processing capacity, the identification of an accidental pollutant source is becoming a more popular topic in some fields, such as research of pollution in atmospheric environment, water environment, enclosed spaces, and porous media. Many researchers attempt to use some potential method to detect contaminant source [1, 2].
Most of the methods used retrospective methods. These methods use the sensor observed data to estimate the unknown pollutant source. These methods can be further classified into analytical, optimization, probabilistic, and backward computational fluid dynamics (CFD) methods.
The analytical method can identify a pollutant source based on an analytical solution of a velocity field and a measured concentration distribution. It is often applied in an atmospheric environment or on surface water where a steady flow field is present [3–5].
The optimization method can search for the optimal value by comparing the measured pollutant concentration data with the corresponding calculated data. The linear or nonlinear optimization and the maximum likelihood method are widely adopted. This method requires a large amount of direct calculation and this may lead to a heavy computation load. This method is successfully used to identify the source in an atmospheric environment and in underground water [6–10]. Although this method is quick and efficient for simple situations, it is difficult to obtain the analytical solution for a more complicated flow field. In order to improve the robustness of this method, some researchers attempted to use sensitivity analysis in order to reduce the uncertainty [11, 12].
The probabilistic method uses probability theory to calculate the possibility that a pollution source appears at a certain position [13, 14]. By using this method, the illposed problem of source identification may be transformed into a wellposed problem in the extended statistical space [15, 16]. Many researchers have tried to apply this method to identify an indoor pollutant source and to optimize the placement of sensors [17–23]. Though there are some shortcomings of this method, such as the complexity of the solution of the adjoint state equations, and the fact that the uncertainty of the model cannot be overcome completely, it is indeed a promising method for identifying a source.
The backward CFD method can inversely evaluate the source information by using a negative time step and a reverse flow field. The inverse CFD model is illposed, so regularization technique and stabilization technique should be used to improve its solution stability [24]. Some researchers carried out several studies to improve this stability [25–28]. Currently, the combination of the probabilistic method and the backward CFD method is a significant advance in the identification of a pollutant source in an enclosed space [18–20]. They used backward PDF combined with the quasireversibility method and pseudoreversibility method, respectively, to identify an instantaneous source by using the accurate observed data of a single sensor. Although this research requires an accurate velocity field and exact sensor measured data, it has proved to be promising.
The above endeavors have been promoting the development of source identification. Actual pollutant source emission sometimes is a continuous and not instantaneous process, and the measured data is influenced by sensor noise. So, we develop a method to identify a sudden continuous emission source in a steady velocity field that remains unchanged or little changed during a source emission process by using single sensor information with noise. The application conditions for this method are presented and its assumptions are discussed in detail.
2. Model of Pollutant Transport Process
Assume that a source emits pollutant in a closed space with a constant temperature condition and a small viscous stress and then a gaseous pollutant transport process by air is governed by the ADE and momentum equations [29]: where is the pollutant mass concentration, mg/m^{3}; is the forward time, s; is the velocity vector, m/s; , , and are the velocity components in , , and directions of ; is the diffusion coefficient, m^{2}/s; is the pollutant source which is a function of (, , ); is the source position; is the emission strength, mg/(m^{3}·s); is the fluid density, kg/m^{3}; is the fluid pressure on the micro element; and , and are the body forces in , , and directions, N.
We can obtain an accurate solution of the concentration field by solving the ideal gas equation and (1)~(2). However, this simultaneous solution is very complex and timeconsuming. It will become even more difficult to identify a pollutant source reversely because it should solve and together by using the partially known .
In order to avoid this complex and timeconsuming calculation process, we assume the studied situation satisfies the following two conditions before and after a source appears: (1) is large enough; (2) is approximately uniform; and (3) they remain unchanged or little changed.
With these two assumptions, a pollutant source identification process can be transformed into a relatively simple inverse problem because (2) can be neglected and only (1) needs to be considered. We can estimate the three parameters (, , and ) of the source by solving (1) reversely when and a part of are known.
Equation (1) can be solved numerically. When the pollutant concentration distribution at time is known, the concentration distribution sequent at time can be obtained as follows: where and are the column vectors of the pollutant concentration for each grid node at time and , respectively; is the coefficient matrix of ; and is a diagonally dominant banded sparse matrix; hence, it must have the inverse matrix :
A can be uniquely determined when the computational zone, the mesh, and the discretization scheme are known. We chose the powerlaw scheme to solve the discretized equations of pollutant ADEs in this paper.
3. Source Identification Method
Based on the above assumptions about and , we will develop a method to identify a sudden trace pollutant source by using a single sensor with noise.
A sensor can detect a sudden pollution when its observed concentration rises upon a given threshold value, th. We define two time, and . is the initial emission time when a source begins to emit a pollutant. is the detecting time when a sensor takes to detect the existence of the pollutant. There is a delay time and .
An identification process can be conducted after and are known. We use the idea of multiple parameter multiple hypothesis optimization to identify the unknown three source parameters (, , and ) of source. The measured concentration sequence of sensor in the time interval is specially taken out from and defined by for the identification. There are four required steps to implement this process, as shown in Figure 1.
(A) Multiple Hypothesis of Source Position. A set of hypotheses of source position, () is established firstly. It includes all possible source positions. The second hypothesis assumes that a virtual pollutant source with a constant emission strength, , will appear on every possible hypothetical position , respectively. For a certain , a set of hypotheses, , will be established thirdly. If the backward searching time index is , then .
(B) Optimization Only for . The optimal three parameters of source only for can be searched with the objective function of minimum characteristic distance. It includes the following steps.(a)When and are known, the hypothetical measured concentration sequence, , at sensor position for can be calculated by forwardly solving (4) in with the assumption of ().(b)Assume that remains unchanged or little changed during the pollutant emission process and then the steadystate observed concentration is linear with respect to emission rate [23]; hence, the approximate relative emission strength can be obtained as follows: where is the constant emission strength of the virtual pollutant source, mg/(m^{3}·s); is the estimated emission strength, , mg/(m^{3}·s); is the estimated relative emission strength with the assumption of ; and is the initial concentration, mg/m^{3}.(c)After are calculated by using (5), the concentration at sensor position in the time interval can be calculated forwardly by using (4) and is denoted as the hypothetical measured concentration sequence of sensor : and can be obtained with (6).(d)We defined a characteristic distance to describe the similarity degree between and in : where is the characteristic distance between and .(e)After the set of has been calculated with (7), we define the minimum characteristic distance . will be an objective function to search the optimal three parameters of source only for the position . Let “” express the optimal one in the set of and then and can be optimized accordingly to . Now, the optimal result for the hypothetical position is .
(C) Minimum Characteristic Distance Set. Similarly, the set of minimum characteristic distance, , can be finally obtained for the other hypothetical positions in the set of .
(D) Source Identification Based on Location Probability. A location probability is established to clearly point out the most probable source positions after has been obtained. We use the location probability to show the probability of a source appearing at the th assumed position [14]: where is the measurement standard deviation of sensor in and is the number of the hypothetical source position.
A distribution map of location probability can be obtained after the set of location probabilities, , has been calculated by using (8). This map is directly able to show where the source appears with max location probability. We denoted as the global optimal objective function and let , where the superscript is the node corresponding to the .
After has been searched, the source can be located at ; that is, (8) is used for the source location. Hence, the global optimal source parameter is and this is the final identification result for all hypothetical positions, .
4. Identifiability
The presented identification method is relatively effective because only (1) is used for estimating a source reversely. But its performance of identification could be acceptable only when the assumptions about and are satisfied very well. If not, then an identification error will be generated. In this section, the impact of and on the performance of identification will be discussed in detail in order to obtain the application conditions.
4.1. Velocity Field
The velocity field, , has an important impact on the identification of pollutant source. In our analysis, is assumed to be large enough and remains unchanged or little changed in the whole process, so the source emission or other factors should not change velocity field greatly.
In order to test the change of velocity fields before and after a source happening, we define a synergy degree of velocity fields to quantitatively discuss how to satisfy the above velocity field assumption according to the field synergy principle [30].
For two scalar fields of the same variable in the same space, denoted by the field and , we define a nonnegative real number as the synergy degree of the field and the field for the same variable. is expressed as a percentage: where is the value of variable at position in the field ; is the value of at in the field ; and expresses the maximum absolute value of the scalar in the field .
If is a vector, take its components in every axis direction and obtain their corresponding synergy degrees. The largest component of synergy degrees will be the synergy degree of vector between the field and the field .
and are the velocity fields before and after a sudden source appears. We can obtain the synergy degrees of velocity fields with the definition of (9). If , then , which represents the velocity field, is not changed greatly; If , then , which represents the velocity field, is changed. The greater the value of is, the greater the velocity field will be changed.
One of the prerequisites to using our identification method is ; hence, should be required to be as small as possible. In order to discuss the impact of or on , we use a 2D simulation case, as shown in Figure 2, for the analysis.
It is a rectangular enclosed space, 560 mm long and 360 mm wide, with a 20 mm wide inlet and a 20 mm wide outlet. There is a pollutant source at , and it is a 10 mm × 10 mm square pollution source. There is a pollutant sampling point at . The direction of gravity is the direction of .
In our study, CO_{2} is the aim pollutant. Its molecular weight is larger than air. The source begins to release the CO_{2} pollutant with constant emission strength at 0 second after the velocity field has become steady. The source emits a trace amount of the pollutant continuously.
Four cases with different and were analyzed. According to the above definition, we can obtain the synergy degrees of velocity fields before () and after a sudden source appears in the gravity field space for the case in Figure 2 with different inlet velocity and source conditions. They are listed in Table 1.

From Table 1, we can see that(1)the synergy degree of velocity fields will increase with the decrease of . For example, when kg/(m^{2}·s), Φ is changed from 1.47% to 13.12% corresponding to changing from 1 m/s to 0.3 m/s.(2)If we assume that the application condition for the source identification is % (in fact, Φ should be as small one as possible). % in Case 2, in which m/s and kg/(m^{2}·s), and it satisfies the assumption condition of %. While for the same and a smaller inlet velocity ( m/s) in Case 3, its % and it does not satisfy %. But if becomes smaller, for example, kg/(m^{2}·s) in Case 4 with m/s, then its % and it satisfies % again. So, there is a close relationship between and . Therefore, a larger inlet velocity is more ideal to satisfy the assumption condition of application than a smaller one and will ensure that the source could be identified.
4.2. Momentum Equation
In this section, we will discuss in what situation the impact of the momentum equation can be ignored.
If only considering the axis direction (), the onedimensional flow velocity is taken as a fixed value and is very small for the air and then (2) can be written as follows:
According to the assumptions about and , we know that the absolute value of is smaller when is greater; thus, the influence of the second term on the right (gravity and air buoyancy) becomes weaker. If we assume that remains approximately uniform in the space, then the first term on the right of (10) can be ignored. Hence, the flow velocity and the amount of emitted pollutant are the key factors for effectively weakening the influence of the momentum equation.
We denote and as the true concentration and the calculated concentration at in Figure 2, respectively. can be obtained by solving the simultaneous equations of continuity, advectiondispersion, momentum, and energy equations. can be calculated by only solving (4). Because it is easier and more efficient to calculate than to calculate , we will compare with and determine in what situation is it more suitable to use instead of for the identification.
The above four cases in Table 1 are also used again for the comparison. is used to show the difference of and . Figures 3(a)~6(a) show the comparison results of for the four cases. We let stand for the virtual measured concentration at with an observed noise. We assume that the source is unknown and estimate it by using the presented method. Figures 3(b)~6(b) show the location results for the four cases. The deep red region indicates higher location probability and the blue region indicates that zero location probability.
(a) Concentration comparison
(b) Location result
We draw some conclusions from Figures 3~6.(1)The difference of and is very small in most time when the inlet velocity is m/s, as shown in Figures 3(a)~4(a). This implies that the influence of the momentum equation can be negligible. In this situation, we can use instead of for identification. The location results in Cases 1 and 2 can be acceptable because the deep red region is concentrated near the inlet and near the actual source position, as shown in Figures 3(b) and 4(b). But the difference of and are smaller in Case 1 than Case 2; then the location result in Case 1 is better than the one in Case 2.(2)When m/s, the difference between and becomes large as shown in Figure 5(a). This implies that the momentum equation has played a leading role in the CO_{2} transmission and can not be ignored. We can not obtain a reasonable location result, as shown in Figure 5(b). However, if the source strength becomes tracer, for example, kg/(m^{2}·s) in Case 4, then we can obtain a more reasonable location result, as shown in Figure 6(b). This is because in Case 4 is equal to 9.34% and it is less than the one in Case 3.
(a) Concentration comparison
(b) Location result
(a) Concentration comparison
(b) Location result
(a) Concentration comparison
(b) Location result
4.3. Application Conditions
According to the above analyses, we know the following.(1)The value of the inlet velocity or the velocity field is the key to satisfy the application assumptions of the presented method. An inlet velocity should be large enough to ensure that the source is effectively identified.(2)The inlet velocity and the source emission strength should be matching. The source emission process should not disturb the velocity field. So, a source with strong emission strength could be located in a large velocity field, but it cannot be located accurately in a small velocity field. For a small velocity field, a source needs to be one with small emission strength to be estimated accurately.
Therefore, the identification method is actually suitable for estimating a trace continuous pollutant source in a steady velocity field. In the actual application, we can refer to this idea to analyze a studied situation and determine its application conditions.
5. Experiment Study
5.1. Experimental System
An experimental device was set up to test the proposed method. CO_{2} is used as the aim pollutant. CO_{2} gas injected into the enclosed space from a point is supplied by a compressed gas cylinder through a pressure relief valve and a mass flow controller. A QIC20 mass spectrometer made by the UK company Hiden Analytical Ltd. was used to measure CO_{2} concentrations in real time at the sampling point.
The photo of the enclosed space is shown in Figure 7 and its size is given in Figure 2. It is a flat cavity topped and bottomed by two plexiglass flat panels and fourside sealing strips. Due to the thickness of the enclosed space being only 8 mm and relatively very thin, the flow velocity, which is formed by a small centrifugal fan, in the thickness direction can be considered approximately zero; thus, the flow field is close to a 2D case.
5.2. Source Identification
During the experiment, the inlet velocity was 1.68 m/s. CO_{2} concentration curve measured at is shown in Figure 8 when the CO_{2} was suddenly injected from . The source emission strength was constant and its actual emission strength was 16.37 mg/s.
For the experimental space in Figure 7, we know that(1) m/s and it is larger than 1.0 m/s of Case 1 in Table 1. At the same time, its is larger than 1.47% and is steady. Hence the assumption of is satisfied very well;(2) mg/s and it is less than 50 mg/s in Case 1. This means that the emission amount of source is small and the assumption of is satisfied.
Therefore, this experimental situation satisfies the application conditions in Section 4.3.
We apply the identification method to identify the source reversely by using the experimental data in Figure 8. The other parameters are s, s, mg/m^{3}, and s.
The estimated emission strength is 20.862 mg/s, and the final distribution map of location probability is shown in Figure 9. From Figure 9, we know that the position having the largest location probability is very near to the actual pollutant source position.
6. Conclusions
A method to identify a trace continuous source in a steady flow field was proposed in this paper. The identification method is realized by globally searching the optimal values with the objective function of maximum location probability.
In order to ensure the performance of identification, the assumptions about and are discussed in detail and the application conditions for the identification method are obtained. A synergy degree of velocity fields is introduced to analyze the application conditions of velocity field quantitatively. The value of should be as small as possible in order to satisfy this assumption of the velocity field. Based on this premise, we draw the following application conditions from the simulation studies:(1)the value of the velocity field (or the inlet velocity) is one key factor. A larger is more ideal in satisfying the application conditions;(2)The source emission strength is another key factor in addition to . The source emission process can not disturb greatly and and should be matching each other. Therefore, the identification method is suitable for estimating a trace continuous pollutant source in a steady velocity field.
In order to investigate the practical application effect of the above studies, the experimental system was set up and subsequently a sudden pollution experiment was conducted. The presented identification method was used. The studies show that the sudden trace continuous source can be successfully located when the experimental situation satisfies the application conditions.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
References
 H. Deng, Z. Dai, A. Wolfsberg, Z. Lu, M. Ye, and P. Reimus, “Upscaling of reactive mass transport in fractured rocks with multimodal reactive mineral facies,” Water Resources Research, vol. 46, no. 6, Article ID W06501, 2010. View at: Publisher Site  Google Scholar
 H. Deng, P. H. Stauffer, Z. Dai, Z. Jiao, and R. C. Surdam, “Simulation of industrialscale CO_{2} storage: multiscale heterogeneity and its impacts on storage capacity, injectivity and leakage,” International Journal of Greenhouse Gas Control, vol. 10, pp. 397–418, 2012. View at: Google Scholar
 L. C. Thomsona, B. Hirstb, G. Gibsona et al., “An improved algorithm for locating a gas source using inverse methods,” Atmospheric Environment, vol. 41, pp. 1128–1134, 2007. View at: Google Scholar
 T. K. Flesch, J. D. Wilson, L. A. Harper, and B. P. Crenna, “Estimating gas emissions from a farm with an inversedispersion technique,” Atmospheric Environment, vol. 39, no. 27, pp. 4863–4874, 2005. View at: Publisher Site  Google Scholar
 M. A. Islam, “Application of a Gaussian Plume model to determine the location of an unknown emission source,” Water, Air, and Soil Pollution, vol. 112, no. 34, pp. 241–245, 1999. View at: Publisher Site  Google Scholar
 A. Khemka, C. A. Bouman, and M. R. Bell, “Inverse problems in atmospheric dispersion with randomly scattered sensors,” Digital Signal Processing, vol. 16, no. 5, pp. 638–651, 2006. View at: Publisher Site  Google Scholar
 M. D. Sohn, P. Reynolds, N. Singh, and A. J. Gadgil, “Rapidly locating and characterizing pollutant releases in buildings,” Journal of the Air and Waste Management Association, vol. 52, no. 12, pp. 1422–1432, 2002. View at: Google Scholar
 A. C. Bagtzoglou and J. Atmadja, “Marchingjury backward beam equation and quasireversibility methods for hydrologic inversion: application to contaminant plume spatial distribution recovery,” Water Resources Research, vol. 39, no. 2, pp. 1038–1052, 2003. View at: Google Scholar
 S. M. Gorelick, B. Evans, and I. Remson, “Identifying sources of groundwater pollution: an optimization approach,” Water Resources Research, vol. 19, no. 3, pp. 779–790, 1983. View at: Google Scholar
 B. J. Wagner, “Simultaneous parameter estimation and contaminant source characterization for coupled groundwater flow and contaminant transport modelling,” Journal of Hydrology, vol. 135, no. 1–4, pp. 275–303, 1992. View at: Google Scholar
 A. C. Bagtzoglou and J. Atmadja, “Marchingjury backward beam equation and quasireversibility methods for hydrologic inversion: application to contaminant plume spatial distribution recovery,” Water Resources Research, vol. 39, no. 2, pp. 1001–1014, 2003. View at: Google Scholar
 M. Skliar and W. F. Ramirez, “Source identification of the distributed parameter transport processes,” in Proceedings of the American Control Conference, pp. 2228–2232, June 1997. View at: Google Scholar
 P. Kathirgamanathan, R. Mckibbin, and R. I. Mclachlan, “Source term estimation of pollution from an instantaneous point source,” Research Letters in Information Mathematic Science, vol. 3, pp. 59–67, 2002. View at: Google Scholar
 M. Lehning, D. R. Shonnard, D. P. Y. Chang, and R. L. Bell, “An inversion algorithm for determining areasource emissions from downwind concentration measurements,” Journal of the Air and Waste Management Association, vol. 44, no. 10, pp. 1204–1213, 1994. View at: Google Scholar
 L. P. Pang and H. Q. Qu, “Contaminant source location method based on multi hypothesis position,” in Proceedings of the International Conference on Indoor Air, p. pa260, 2011. View at: Google Scholar
 M. D. Sohn, P. Reynolds, N. Singh, and A. J. Gadgil, “Rapidly locating and characterizing pollutant releases in buildings,” Journal of the Air and Waste Management Association, vol. 52, no. 12, pp. 1422–1432, 2002. View at: Google Scholar
 A. Keats, E. Yee, and F.S. Lien, “Bayesian inference for source determination with applications to a complex urban environment,” Atmospheric Environment, vol. 41, no. 3, pp. 465–479, 2007. View at: Publisher Site  Google Scholar
 J. Arvelo, A. Brandt, R. P. Roger, and A. Saksena, “An enhancement multizone model and its application to optimum placement of CBW sensors,” ASHRAE Transactions, vol. 108, no. 2, pp. 818–825, 2002. View at: Google Scholar
 T. Zhang and Q. Chen, “Identification of contaminant sources in enclosed spaces by a single monitor,” in Proceedings of the IBPSA Building Simulation Conference, pp. 1477–1484, 2007. View at: Google Scholar
 T. F. Zhang and Q. Chen, “Identification of contaminant sources in enclosed environments by inverse CFD modeling,” Indoor Air, vol. 17, no. 3, pp. 167–177, 2007. View at: Publisher Site  Google Scholar
 T. Zhang, H. Li, and S. Wang, “Inversely tracking indoor airborne particles to locate their release sources,” Atmospheric Environment, vol. 55, pp. 328–338, 2012. View at: Publisher Site  Google Scholar
 M. Girault, D. Maillet, F. Bonthoux et al., “Estimation of timevarying pollutant emission rates in a ventilated enclosure: inversion of a reduced model obtained by experimental application of the modal identification method,” Inverse Problems, vol. 24, no. 1, Article ID 015021, 2008. View at: Publisher Site  Google Scholar
 X. Liu and Z. Zhai, “Probabilitybased inverse modeling algorithm for indoor pollutant source tracking,” in Proceedings of the IBPSA Building Simulation Conference, pp. 810–817, 2007. View at: Google Scholar
 M. F. Snodgrass and P. K. Kitanidis, “A geostatistical approach to contaminant source identification,” Water Resources Research, vol. 33, no. 4, pp. 537–546, 1997. View at: Google Scholar
 R. Lattes and J. L. Lions, The Method of QuasiReversibility: Applications to Partial Differential Equations, Elsevier, New York, NY, USA, 1969.
 A. Y. Sun, S. L. Painter, and G. W. Wittmeyer, “A robust approach for iterative contaminant source location and release history recovery,” Journal of Contaminant Hydrology, vol. 88, no. 34, pp. 181–196, 2006. View at: Publisher Site  Google Scholar
 T. H. Skaggs and Z. J. Kabala, “Recovering the history of a groundwater contaminant plume: method of quasireversibility,” Water Resources Research, vol. 31, no. 11, pp. 2669–2673, 1995. View at: Publisher Site  Google Scholar
 A. C. Bagtzoglou, D. E. Dougherty, and A. F. B. Tompson, “Application of particle methods to reliable identification of groundwater pollution sources,” Water Resources Management, vol. 6, no. 1, pp. 15–23, 1992. View at: Publisher Site  Google Scholar
 C.H. Huang, J.X. Li, and S. Kim, “An inverse problem in estimating the strength of contaminant source for groundwater systems,” Applied Mathematical Modelling, vol. 32, no. 4, pp. 417–431, 2008. View at: Publisher Site  Google Scholar
 D. John and J. Anderson, Computational Fluid Dynamics: The Basics with Applications, McGrawHill, New York, NY, USA, 1995.
Copyright
Copyright © 2014 Hongquan Qu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.