Method for mapping the concentration of an analyte in an environment

Abstract

A method for estimating a mapping of the concentration of an analyte in an environment uses sensors distributed in the environment. Each sensor generates a measurement of the analyte concentration at various measurement instants, which measurements are carried out by each sensor at each measurement instant, forming an observation vector, each term of which corresponds to a measurement arising from a sensor. The environment is spatially meshed with a plurality of mesh cells. The analyte concentration at each mesh cell, at each measurement instant, forms a “state vector,” each term of which corresponds to an analyte concentration in a mesh cell. A “global bias” is determined and used to correct the state vector to obtain a “debiased state vector.” The state vector is also corrected by a local correction vector as a function of a correction vector.

Claims

1. A method for estimating a mapping of a concentration of an analyte in an environment, using sensors distributed in the environment: each sensor generating a measurement of an analyte concentration at various measurement instants, the measurements carried out by each sensor at each measurement instant forming an observation vector, each term of which corresponds to a measurement arising from a sensor of the sensors; the environment being spatially meshed with a plurality of mesh cells, the analyte concentration at each mesh cell, at each measurement instant, forming a state vector, each term of which corresponding to an analyte concentration in a mesh cell; the method comprising, using a processor: a) obtaining a measured observation vector at a measurement instant, using the measurements performed by each sensor; b) obtaining a state vector at the measurement instant and, using the state vector, estimating an observation vector at the measurement instant; c) comparing the estimation of the observation vector, obtained in b), with the measured observation vector resulting from a), and, on the basis of the comparison, determining a global bias at the measurement instant, the global bias being a scalar representative of a comparison between several terms of the estimated observation vector and of the measured observation vector respectively, d) correcting the state vector obtained in b) with the global bias obtained in step c), so as to obtain a debiased state vector at the measurement instant; e) on the basis of the debiased state vector obtained in d), obtaining a debiased estimation of observation vector at the measurement instant; f) comparing the debiased estimation of the observation vector resulting from e) with the measured observation vector resulting from a), and, on the basis of the comparison, determining of a local correction vector; and g) updating the state vector at the measurement instant, the latter being replaced with a sum of the debiased state vector resulting from d) and the local correction vector resulting from f), the updating of the state vector making it possible to estimate the mapping of the concentration of the analyte in various mesh cells.

2. The method according to claim 1, wherein b), the observation vector is estimated by applying an observation matrix to the state vector.

3. The method according to claim 1, wherein c) comprises: ci) establishing comparisons between various terms of the observation vector estimated in b) and of the observation vector measured in a); cii) calculating an average or median value of each comparison resulting from ci); and ciii) obtaining the global bias on the basis of an average or a median value resulting from cii).

4. The method according to claim 1, wherein d) comprises subtracting of the global bias from each term of the state vector.

5. The method according to claim 1, wherein in e), the debiased estimation of the observation vector is obtained by applying an observation matrix to the debiased state vector.

6. The method according to claim 1, wherein f) comprises: fi) establishing a comparison vector, resulting from a comparison, term by term, between the observation vector resulting from a) and the debiased estimation of the observation vector resulting from e); fii) taking into account of a gain matrix; and fiii) applying the gain matrix to the comparison vector so as to form the correction vector.

7. The method according to claim 2, wherein, following g), the method further comprises: h) iteratively updating the state vector, each iteration being associated to an iteration rank, the method further comprising: hi) taking into account a gain matrix corresponding to the iteration rank; hii) determining a comparison vector, associated with the iteration rank, by comparing the observation vector resulting from a) with a vector resulting from the application of the observation matrix to the state vector resulting from g), or to the state vector resulting from a previous iteration; hiii) applying the gain matrix of hi) to the comparison vector determined in hii), so as to obtain a local correction vector associated with the iteration rank; hiv) updating the state vector, the latter being replaced with a sum of the state vector resulting from g), or from a previous iteration, with the local correction vector resulting from hiii); and hv) repeating hi) to hiv) or stopping of the iteration.

8. The method according to claim 6, wherein each term of the gain matrix is associated with a mesh cell and with a sensor, a value of the term being all the higher the closer the mesh cell is to the sensor.

9. The method according to claim 7, wherein each term of the gain matrix is associated with a mesh cell and with a sensor, a value of the term being all the higher the closer the mesh cell is to the sensor.

10. The method according to claim 1, wherein b) comprises forming the state vector by using a model established on the basis of data relating to the road traffic in the environment, of the topography of the environment as well as of meteorological data relating to the environment, the model resulting in a concentration of the analyte at the level within each mesh cell.

11. The method according to claim 1, further comprising: taking into account of a later state vector with respect to a later instant following the measurement instant; and correcting the later state vector as a function of the state vector updated, at the measurement instant, during g).

12. The method according to claim 7, comprising the following steps: taking into account of a later state vector with respect to a later instant following the measurement instant; and correcting the later state vector as a function of the state vector updated, at the measurement instant, during step h).

13. The method according to claim 2, wherein each term of the observation matrix is associated with a mesh cell and with a sensor, a value of the term being all the higher the closer the mesh cell is to the sensor.

14. The method according to claim 5, wherein each term of the observation matrix is associated with a mesh cell and with a sensor, a value of the term being all the higher the closer the mesh cell is to the sensor.

15. The method according to claim 7, wherein each term of an observation matrix is associated with a mesh cell and with a sensor, a value of the term being all the higher the closer the mesh cell is to the sensor.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) FIG. 1A is the plan of an urban area under study, within which measurement sensors are distributed.

(2) FIG. 1B is a mapping of an analyte, in this instance nitrogen dioxide, this mapping being obtained by application of a model of OSPM type.

(3) FIG. 2 shows the main steps of a method according to the disclosure.

(4) FIG. 3A represents the mapping of FIG. 1B after taking a global bias into account.

(5) FIG. 3B shows a two-dimensional representation of the positive terms of a local correction vector.

(6) FIG. 3C shows a two-dimensional representation of the negative terms of a local correction vector.

DETAILED DESCRIPTION

(7) FIG. 1A represents a plan of an urban area in which a two-dimensional mapping of the concentration of an analyte is modelled according to a model known from the prior art, for example, the OSPM model mentioned previously. In this example, the analyte is a nitrogen dioxide molecule. In general, the analyte is a chemical molecule or a particle whose dispersion in the environment one wishes to know, that is to say a spatial distribution of its concentration or of its quantity. It may, in particular, be an analyte arising from road traffic. FIG. 1B represents a mapping of nitrogen dioxide in the streets of the urban area represented in FIG. 1A. In FIG. 1B, the grey levels correspond to a nitrogen dioxide concentration expressed in ppb.

(8) On the basis of the mapping modelled in FIG. 1B, it is possible to define a geographical meshing of the urban area, and to form, on the basis of the model, a vector, the so-called state vector M(t), each term M.sub.m(t) of which corresponds to a nitrogen dioxide concentration modelled at the level of a mesh cell 20.sub.m at an instant t, for example, at the level of each mesh cell centre. The index m is a strictly positive integer designating a mesh cell. The dimension of the state vector M(t) is (N.sub.m, 1), where N.sub.m represents the number of mesh cells considered. Each term of the state vector is obtained by the application of a predictive model, such as the OSPM model described with regard to the prior art, by taking into account data linked with the urban traffic, meteorological parameters, such as the temperature and/or the speed of the winds, as well as the three-dimensional topography of the environment, for example, the geometry of the streets as well as the height of the buildings between each street.

(9) Represented in FIG. 1A, in the form of dots, are simulated locations of sensors (10.sub.p), the index p being a strictly positive integer designating a sensor. In the example represented, each sensor is a nitrogen dioxide sensor, measuring a concentration c.sub.p(t) of this analyte at each measurement instant t. The measured concentrations c.sub.p(t) form a vector C(t), the so-called observation vector, at the measurement instant, each term of which is a concentration measured by a sensor at the measurement instant. The dimension of the vector C(t) is (N.sub.p, 1), where N.sub.p represents the number of sensors considered. The sensors are connected to a processor, for example, a microprocessor, the latter being programmed to execute instructions to implement the method described in this disclosure.

(10) On the basis of the state vector M(t) and of the observation vector C(t), the method described hereinbelow is aimed at updating the state vector, in such a way as to increase the precision of the mapping of the urban area considered, by taking into account the measurements performed by each sensor. Indeed, certain local features, not taken into account by the model, may have a local influence on the distribution of the analyte. This may, in particular, entail a bottleneck. This disclosure makes it possible to take them into account. The main steps of the method are represented in FIG. 2.

(11) Preferably, the distance between two adjacent sensors is less than 500 m, or indeed than 200 m. Indeed, the method described hereinbelow is all the more effective the higher the number of sensors. In terms of number of sensors per unit area, preferably, the number of sensors is greater than 2 or indeed 3 per km.sup.2. On the scale of a city, recourse to about ten or about twenty sensors is not sufficient to perform a sufficiently effective updating of the mapping.

(12) Step 100: acquisition of the data.

(13) This entails obtaining the state vector M(t) on the basis of the modelled mapping and of the observation vector C(t) on the basis of the sensors. FIG. 1B corresponds to a two-dimensional representation of the state vector M(t), which representation is obtained by establishing a correspondence between each term of this vector and a two-dimensional spatial coordinate M.sub.m(t) corresponding to a mesh cell 20.sub.m. In this example, the observation vector C(t) is obtained by simulation on the basis of concentrations established, at the level of each sensor 10.sub.p, on the basis of the model. A bias is added, as is an error term, the latter according to a Gaussian law.

(14) Step 110: estimation of the observation vector on the basis of the state vector.

(15) This entails estimating an observation vector, denoted Ĉ(t), on the basis of the state vector M(t). The estimation of the observation vector can be obtained by applying a matrix H, the so-called observation operator, to the state vector M(t), in the form of a matrix product. The matrix H makes it possible to spatially interpolate the measured data forming the observation vector C(t) so as to obtain, on the basis of the state vector, estimations (t) of the nitrogen dioxide concentration at the level of each sensor 10.sub.p. The matrix H is of dimension (N.sub.p, N.sub.m). With each row and with each column of the matrix H are respectively associated a sensor 10.sub.p and a mesh cell 20.sub.m. The estimation of the observation vector Ĉ(t) is obtained according to the expression: Ĉ(t)=H×M(t), (1) where x designates the matrix product.

(16) The terms of the matrix H(p, m) depend on the relative position of a sensor 10.sub.p with respect to the various mesh cells 20.sub.m. When the sensor 10.sub.p coincides with the centre of a mesh cell 20.sub.m, the row H(p, .) of the matrix H corresponding to the sensor 10.sub.p comprises only 0's, except at the level of the column corresponding to the mesh cell. In a general manner, the matrix is such that on a row H(p, .) corresponding to a sensor, the term of each column is all the higher when the column is associated with a mesh cell situated in proximity to the sensor. Preferably, the terms of the matrix H lie between 0 and 1.

(17) Step 120: comparison between the observation vector C(t) and its estimation Ĉ(t) and calculation of a global bias.

(18) In the course of this step, each term of the observation vector C(t) is compared with the term, corresponding to the same sensor 10.sub.p, of the estimation of the observation vector Ĉ(t) resulting from step 110. The comparison can take the form of a term-by-term subtraction or ratio.

(19) On the basis of the comparison, a global bias ε(t) is calculated, the bias representing, at the measurement instant t, a global comparison between the observation vector C(t) and its estimation Ĉ(t). The global bias is a scalar quantity. It may, in particular, be determined on the basis of an average or of a median of a comparison, term by term, of the vectors C(t) and Ĉ(t). For example,

(20) $\begin{array}{cc}\mathrm{.Math.}\left(t\right)=\frac{1}{N_p}\underset{p=1}{\overset{N_p}{.Math.}}\left(C_p\left(t\right)-{\hat{C}}_p\left(t\right)\right),& \left(2\right)\end{array}$
where C.sub.p(t) and Ĉ.sub.p(t) are respectively a term of rank p of the vectors C(t) and Ĉ(t), corresponding to one and the same sensor 10.sub.p.

(21) Step 130: debiasing of the state vector.

(22) In the course of this step, the state vector M(t) is corrected of the global bias ε(t), by subtracting the global bias from each term of the state vector. A so-called debiased state vector denoted M′(t) is then obtained.

(23) Thus, M′(t)=M(t)−E(t) (3) where E(t) is a bias vector, of dimension (N.sub.m, 1), each term of which is equal to the global bias ε(t). The term debiased signifies unbiased. The term debiasing signifies removal of the bias.

(24) This step forms a first correction of the state vector, on the basis of a global bias calculated on the basis of the observations obtained by the sensors 10.sub.p. Such a bias may be due to emissions affecting the whole urban area under study, and originating, for example, from urban heating, or a diffuse pollution. The inventors have observed that taking a global bias such as this into account allowed an appreciable improvement in the precision of the state vector M(t).

(25) FIG. 3A shows a two-dimensional representation of the debiased state vector M′(t). In this example, the value of the bias rises to 9.2 μg/m.sup.3. To each point of this figure there corresponds a mesh cell 20.sub.m and a term of the M′.sub.m(t) debiased state vector.

(26) Step 140: Debiased estimation of the observation vector.

(27) In the course of step 140, the debiased state vector M′(t), resulting from step 130, is fitted against the measurements arising from the sensors 10.sub.p. Accordingly, a so-called debiased estimation, denoted Ĉ′(t), of the observation vector is calculated by applying the observation operator H according to the expression Ĉ′(t)=H×M′(t) (4), in a manner analogous to step 110.

(28) Step 150: Fitting of the debiased state vector M′(t) to the measurements.

(29) In the course of step 150, a term-by-term comparison is performed between the debiased estimation of the observation vector Ĉ′(t), resulting from step 140, with the observation vector C(t) established during step 100. The comparison can take the form of a subtraction or of a ratio. From this comparison is formed a local comparison vector comp(t). The local comparison vector comp(t) is of dimension (N.sub.p, 1). In contradistinction to the debiasing step (steps 120 and 130), the comparison is a vector quantity. Thus, each term comp.sub.p(t) of the local comparison vector is such that comp.sub.p(t)=C′.sub.p(t)−Ĉ.sub.p(t) (5), the index p representing the rank of each term, p lying between 1 and N.sub.p. In contradistinction to the bias vector E(t), the terms of the local comparison vector comp(t) may differ from one another, and are mutually independent.

(30) Step 160: updating of the state vector.

(31) After having formed the subject of a debiasing, during step 140, the state vector forms the subject of a second, so-called local, correction based on the local comparison vector comp(t) formed during step 150. A matrix, the so-called gain matrix K, makes it possible to perform a weighting of the correction to be made as a function of the distance of a mesh cell 20.sub.m with respect to each sensor 10.sub.p. The gain matrix is of dimension (N.sub.m, N.sub.p). With each row and with each column of the gain matrix are respectively associated a mesh cell 20.sub.m and a sensor 10.sub.p. The terms of a row K(m, .), corresponding to a mesh cell 20.sub.m, are all the higher the closer a sensor 10.sub.p, corresponding to a column, is to the mesh cell. The terms K(m, p) of a gain matrix are preferably less than or equal to 1.

(32) The updating of the state vector is performed according to the following expression:
M*(t)=M′(t)+K×comp(t)=M′(t)+K×(C(t)−Ĉ′(t))  (6)
=M′(t)+K×(C(t)−H×M′(t))  (6′)
M*(t) corresponds to the updated state vector, making it possible to obtain a more realistic mapping of the pollutant.

(33) This operation is equivalent to applying a local correction vector corr(t) to the debiased state vector M′(t) so as to obtain a corrected (or updated) state vector M*(t). The local correction vector is of dimension (N.sub.p, 1) and corresponds to the application of the gain matrix to the local comparison vector comp(t), according to the expression:
corr(t)=K×(C(t)−Ĉ′(t))=K×comp(t)=K×(C(t)−H×M′(t))  (6″)

(34) In contradistinction to the debiasing step, the local correction vector is not uniform.

(35) In the course of this step, the correction of the state vector is not uniform, as during the debiasing, but differs from one term of the state vector to another. FIGS. 3B and 3C illustrate this aspect, and represent respectively the positive and negative terms of the correction vector corr(t) at the various mesh cells of the mapping. It is observed that the correction is local, the correction being more significant in certain parts than in others. It may be negative in certain parts, and positive in other parts.

(36) The method makes it possible, through a sufficiently high density of sensors, to obtain a mapping taking account of local features of the traffic, for example, the occurrence of a bottleneck. The combination between the taking into account of a global bias, followed by a local correction step, makes it possible to improve the spatial resolution of the mapping arising from the updated state vector. In particular, it makes it possible to take account of local evolutions, affecting only a few mesh cells 20.sub.m. The mapping obtained is thus more reactive in regard to the occurrence of local features.

(37) According to one embodiment, forming the subject of step 170 the updating of the state vector K is performed iteratively, by modifying the gain matrix at each iteration. Let n be the rank of each iteration and let K.sub.n be the gain matrix associated with each iteration, step 170 comprises an updating of the state vector resulting from step 160, or from a previous iteration n−1 in such a way that:
M.sub.n*(t)=M.sub.n−1(t)+K.sub.n×comp.sub.n(t)  (7)
With comp.sub.n(t)=C(t)−H×M.sub.n−1*(t)  (8), where: comp.sub.n(t) is a comparison vector associated with the iteration of rank n; M.sub.n*(t) is the state vector updated in the course of the iteration of rank n; C(t) is the previously defined measured observation vector; and H is the previously defined observation operator.

(38) Step 170 is repeated until an iteration criterion is attained. Such a criterion may be a predetermined number N.sub.n of iterations, or a sufficiently small disparity between two successive updates of the state vector M.sub.n*(t), M.sub.n+1*(t).

(39) Each gain matrix K.sub.n can be determined in the course of each iteration n, as a function of a weight w.sub.m,p.sup.n assigned to each iteration, the indices m and p representing respectively a row and a column of the gain matrix K.sub.n. The weight is defined according to the following expression:

(40) $\begin{array}{cc}{w}_{m,p}^n=\frac{R_{n,p}^2-r_{m,p}^2}{R_{n,p}^2+r_{m,p}^2},\mathrm{where}\text{:}& \left(9\right)\end{array}$ R.sub.n,p.sup.2 is a maximum radius of influence associated with each sensor 10.sub.p; for example, the radius of influence of a sensor disposed in the middle of a town square may be higher than the radius of influence of a sensor disposed in a narrow street; and r.sub.m,p is a distance between a sensor 10.sub.p and a mesh cell 20.sub.m.
The value of each term K.sub.n(m, p) is then such that:

(41) $\begin{array}{cc}{K}_n\left(m,p\right)=0\mathrm{if}{r}_{m,p}>R_{n,p}& \left(10\right)\\ K_n\left(m,p\right)=\frac{w_{m,p}^n}{\left({\alpha }^2-\underset{p}{.Math.}{w}_{m,p}^n\right)}\mathrm{if}{r}_{m,p}\le R_{n,p},& \left(11\right)\end{array}$
where:
α is a term representing an error between the observations and the model. α is a predefined scalar and may, for example, be equal to 0.1.

(42) Thus, with each sensor 10.sub.p is associated a neighbourhood V.sub.n,p, whose extent depends on the maximum radius of influence R.sub.n,p associated with the sensor 10.sub.p. It is considered that the concentrations of the analyte in the mesh cells 20.sub.m that are situated in this neighbourhood are impacted by the measurement arising from the sensor 10.sub.p. The more the iteration rank increases, the more the maximum radius of influence R.sub.n,p corresponding to one or more sensors 10.sub.p decreases. For example, it is assumed that the maximum radius of influence is identical at each sensor R.sub.n,p=R.sub.n. During the first iteration (n=1), R.sub.n=1 is fixed at 500 metres. During the second and third iterations R.sub.n=2 and R.sub.n=3 are fixed at 300 m and 100 m respectively.

(43) According to a variant, the neighbourhood V.sub.n,p associated with a sensor 10.sub.p, that is to say the mesh cells 20.sub.m at the level of which the concentration may be influenced by a measurement performed by the sensor, is not circular, but exhibits a predetermined shape, taking account of the topography, and, in particular, the presence of buildings around the sensor and/or of the dimensions of a street in which the sensor is placed. The neighbourhood of a sensor situated in a street may, for example, extend in a significant manner in a direction parallel to the axis of the street and in a lesser manner in a direction perpendicular to the axis of the street.

(44) According to one embodiment, on the basis of an updated state vector, whether it be the state vector M*(t) updated during step 160 or a state vector M.sub.n*(t) updated iteratively in the course of step 170, the method can comprise a step 200 of forecasting the state vector at a later instant t+dt following the measurement instant t. Accordingly, use is made of a state vector M(t+dt) provided by the model, in this instance the OSPM model. The time interval dt can be of the order of an hour. The state vector M(t+dt) can then be corrected using the state vector updated at the measurement instant t, according to the following expression:
M*(t+dt)=(M(t+dt)−M(t))+M.sub.n*(t)  (12), or
M*(t+dt)=(M(t+dt)−M(t))+M*(t)  (12′)

(45) Thus, the local correction performed on the model M(t+dt) depends on a variation between the state vectors at the respective measurement instants t and t+dt, and on the state vector updated at the measurement instant t, whether it be M.sub.n*(t) or M*(t). It is observed that the correction of the later state vector does not require any new measurements, and is performed with respect to the state vector updated at the measurement instant.

(46) Although described with regard to nitrogen dioxide, the present disclosure will be able to be implemented with other analytes, and, in particular, with polluting molecules or particles. Moreover, in the above example, the state vector established is of OSPM type, but other models known to the person skilled in the art may be applied to form the state vectors at each measurement instant.