Robert E Alvarez
email: ralvarez@aprendtech.com

Conditions for the invertibility of dual energy data

Keywords: spectral x-ray, dual energy, energy selective, inverse

come

Introduction

In the Alvarez-Macovski [START_REF] Alvarez | Energy-selective reconstructions in X-ray computerized tomography[END_REF] method for extracting x-ray spectral data, the attenuation coefficient is expanded as a linear combination of functions of energy multiplied by basis set coefficients that depend on the composition of the material at points within the object. The line integrals of the basis set coefficients are then computed by inverting the transformation from the measurements with multiple spectra. Summarizing the measurements by the vector L and the line integrals of the basis set coefficients as the vector A, we need to show that the transformation L(A) is invertible. This question is examined in this paper. In particular, the transformation for a noise-free system using two measurement spectra with a two function basis set, i.e. a dual energy system, is studied.

Invertibility is of more than mathematical interest since spectra that are near non-invertible can lead to a large multiplication of noise in the output images. In the past, dual energy systems used measurements with two different x-ray tube voltages or a fixed tube voltage and filtering with two materials both with energy integrating detectors. Experience showed that these were invertible. However, the introduction of photon counting detectors into clinical systems [START_REF] Taguchi | Vision 20/20: Single photon counting x-ray detectors in medical imaging[END_REF] leads to potential invertibility problems due to the imperfections in these detectors [START_REF] Overdick | Towards direct conversion detectors for medical imaging with X-rays[END_REF] such as pileup, Kradiation escape, charge sharing, and others. My recent paper [START_REF] Alvarez | Non-invertibility of spectral x-ray photon counting data with pileup[END_REF] showed that a large pileup factor can lead to non-invertibility for a three dimension A-vector with three bin pulse height analysis and that this non-invertibility leads to large increases in noise.

I discussed invertibility in other papers. In early work [START_REF] Alvarez | Energy dependent information in x-ray imaging-part 1 the vector space description[END_REF][START_REF] Lehmann | Energy Selective Radiography: A Review[END_REF], I showed that a general, two dimensional inversion theorem [START_REF] Fulks | Advanced Calculus: An Introduction to Analysis[END_REF] could be applied to the dual energy transformation. The theorem requires that the transformation be invertible on a closed contour in the first quadrant of the A-vector plane and have a non-zero Jacobian determinant inside the contour. I showed that the transformation is invertible on the contour if the two spectra have different maximum energies. The proof was for ordinary energy integrating detectors as well as photon counting detectors without pileup. Later, I showed that the transformation is invertible on the contour for photon counting detectors with pileup [START_REF] Alvarez | Invertibility of spectral x-ray data with pileup-two dimension-two spectrum case[END_REF].

Recently, Levine [START_REF] Levine | Nonuniqueness in dual-energy CT[END_REF] described a different approach to determining invertibility. He studied a system with energy spectra with only three discrete energies, no scatter and a perfect detector. He found that, with appropriately chosen spectral values, the transformation is not one-to-one and is therefore not invertible. The spectra had the same maximum energy so my previous results were not applicable.

My papers did not address the proof that the Jacobian determinant is non-zero within the contour. In this paper, I derive analytical formulas are derived to test for the non-zero condition for a set of measurement spectra and detector. The formulas are applied to several systems including (1) the set proposed by Levine[9] that are not invertible as well as (2) similar spectra that are invertible and (3) invertible voltage switched x-ray tube spectra that have long been used in dual energy imaging. The results are verified by testing the errors with an iterative algorithm that inverts the transformation.

Methods

The transformation

This section reviews the Alvarez-Macovski method [START_REF] Alvarez | Energy-selective reconstructions in X-ray computerized tomography[END_REF], introduces notation, and specifies the transformation that is required to be inverted for the method.

For biological materials, we can approximate the x-ray attenuation coefficient µ(r, E) with a two function basis set [START_REF] Alvarez | Dimensionality and noise in energy selective x-ray imaging[END_REF]

µ(r, E) = a 1 (r)f 1 (E) + a 2 (r)f 2 (E). (1)
In this equation, the coefficients a i (r) are functions only of the material composition at a position r within the object and the functions f i (E) depend only on the x-ray energy E. If there is a high atomic number contrast agent, we can extend the basis set to higher dimensions. Neglecting scatter, the expected value of the number of transmitted photons λ k with an effective measurement spectrum S k (E) is

λ k = ˆSk (E)e -´L µ(r,E)dr dE (2)
where the line integral in the exponent is on a line L from the x-ray source to the detector. The effective measurement spectra include the effects of the source spectrum and the detector response. Using the attenuation coefficient decomposition, Eq. 1, the line integral is

ˆL µ (r, E) dr = A 1 f 1 (E) + A 2 f 2 (E). (3)
where A i = ´L a i (r) dr, i = 1 . . . 2 are the line integrals of the basis set coefficients. If the A i are summarized as the components of the A-vector, A, and the basis functions at energy E by the vector f (E), the expected values in Eq. 2 are

λ k (A) = ˆSk (E)e -A•f (E) dE. (4)
This paper assumes noise free data so the measurements are the expected values. These measurements are summarized by a vector, N, whose components are the expected value of the photon counts with each effective spectrum. Since the body transmission is exponential in A, we can approximately linearize the measurements by taking logarithms. The results is the log measurement vector

L(A) = -log( N(A) N 0 ) (5)
where N 0 is the expected value of the measurements with no object in the beam and the division means that corresponding members of the vectors are divided. Eq. 5 defines a transformation from the log measurements L(A) to the A-vector, A.

The Invertibility Theorem

The following theorem (Fulks [START_REF] Fulks | Advanced Calculus: An Introduction to Analysis[END_REF] page 284) gives conditions for any two dimension transformation to be invertible:

Let F be a continuously differentiable mapping defined on an open region D in E2, with range R in E2, and let its Jacobian be never zero in D. Suppose further that C is a simple closed curve that, together with its interior (recall the Jordan curve theorem), lies in D, and that F is one-to-one on C. Then the image Γ of C is a simple closed curve that, together with its interior, lies in R. Furthermore, F is one-to-one on the closed region consisting of C and its interior, so that the inverse transformation can be defined on the closed region consisting of Γ and its interior. Fig. 1: Closed contour used in proof of invertibility. The contour consists of the two axes, drawn slightly offset here to make them visible, and a quadrant of a circle centered at the origin joining the endpoints of the line segments along the axes.

In Fulks' notation, E2 is a two dimensional Euclidean space, so this theorem is applicable to the transformation L(A) for the two-measurement, two-basis function case. I paraphrase the theorem as

If the Jacobian of a continuously differentiable two dimensional mapping is nonzero throughout an open region D and if the mapping is one to one on a simple closed curve C which lies in D, then the mapping is one to one and invertible on C and its interior.

Fulks uses the term "Jacobian" to refer to the determinant of the Jacobian matrix

J = det ∂L i ∂A j .
I will assume that the converse of the theorem is also true. That is, if the Jacobian determinant is zero in a region, then the transformation is not invertible in that region. This is reasonable given the one-dimension case where a zero derivative implies the transformation is not one-to-one in a region around the zero value and is therefore not invertible.

I discussed the application of this theorem to L(A) in previous papers [START_REF] Alvarez | Energy dependent information in x-ray imaging-part 1 the vector space description[END_REF][START_REF] Lehmann | Energy Selective Radiography: A Review[END_REF][START_REF] Alvarez | Invertibility of spectral x-ray data with pileup-two dimension-two spectrum case[END_REF]. The proofs in these papers used the contour C, shown in Fig. 1, which consists of the A 1 and A 2 axes and a section of a circle centered on the origin and joining the end points of the segments on the axes. In the previous papers, I proved that L(A) is invertible on the A 1 and A 2 axes under very general conditions since objects whose A-vectors are on these axes are composed of only one material.

The previous invertibility proofs on the circular portion of C are not applicable to some of the spectra in this paper. The proof assumes that the maximum energies in the spectra are different but the spectra discussed below have the same three energies but with different weights. The previous argument was that, due to beam hardening, as the radius of the contour becomes larger the transmitted spectrum approaches a peak at the maximum energy. If the maximum energies of the measurement spectra are different, then the transformation will be invertible on the circle. However, if the maximum energies are the same, the transformation becomes non-invertible for large object thicknesses by the No-Linearize theorem (see Section 6.4 of my dissertation [START_REF]Extraction of Energy Dependent Information in Radiography[END_REF] and Chapter 44 of my book [START_REF]Topics in Energy-selective X-ray Imaging[END_REF]). However, the transformation may be invertible for moderate thicknesses if we can guarantee that the circular portion of the contour with a smaller radius is in a region where the Jacobian is non-zero.

Interpretation of the M Matrix and its derivative

For our case, the Jacobian matrix is the matrix M with coefficients

M ij = ∂L i ∂A j . (6)
Since L i = -log λi λi,0 , we can rewrite the elements of M as

M ij = - 1 λ i ∂λ i ∂A j . (7)
Substituting the expected value of measurements Eq. 4 in this equation

M ij = ´fj (E)S i (E)e -A•f (E) dE ´Si (E)e -A•f (E) dE . (8)
Defining the normalized spectrum

ŝi (E) = S i (E)e -A•f (E) ´Si (E)e -A•f (E) dE M ij = ˆfj (E)ŝ i (E)dE = f j (E) i . (9)
That is, M ij is the effective value of the f j (E) basis function in the transmitted measurement spectrum

S i (E)e -A•f (E) .
The derivative of M ij with respect to each of the A-vector components can be computed by differentiating Eq. 8 using the quotient rule

∂Mij ∂A k = -λi[ ´fj (E)f k (E)Si(E)e -A•f (E) dE]+[ ´fj (E)Si(E)e -A•f (E) dE][ ´fk (E)Si(E)e -A•f (E) dE] λ 2 i = -f j (E)f k (E) i -f j (E) i f k (E) i = M ij M ik -f j (E)f k (E) i . (10)
In this equation, λ i is given by Eq. 4 and i is defined in Eq. 9.

The effective values of M depend on the transmitted spectrum and are, for most spectra, a function of A. If the spectra are single delta functions at energy E 1 , E 2 such as from an isotope source or a synchrotron source, then

f j (E)f k (E) i,delta f unction = M ij M ik (11)
and ∂Mij /∂A k = 0. With single delta function spectra, there is no beam hardening and the M ij are constants for all values of A. Notice that Eq. 10 has the same form as a covariance.

Analytical formula for endpoints of zero determinant line

Since M is a 2 × 2 matrix, its determinant is

det (M) = M 11 M 22 -M 12 M 21 (12)
To apply the inversion theorem, we need to determine whether the determinant is zero in the first quadrant of the A plane. The simulations in Sec. 2.6 show that the product terms M ij M pq (A) are approximately planar and do not intersect if the transformation is invertible. If the transformation is not invertible, the M 11 M 22 (A) and M 12 M 21 (A) surfaces do intersect and the curve of intersection is approximately a straight line. One approach to develop analytical formulas for invertibility is to determine the line by points of zero determinant on the A 1 and A 2 axes. the intersections of the the product terms on the A 1 and A 2 axes. If the determinant is zero in the first quadrant so the system is invertible, both of the intersections will be not physically feasible, that is, the product terms at the intersection are less than or equal to zero. Using a linear approximation to the determinant along the

A k , k = 1, 2 axes det (M) (A k ) ≈ det (M) (0) + ∂ ∂A k [M 11 M 22 ] (0) -∂ ∂A k [M 12 M 21 ] (0) A k . ( 13
)
Setting the right hand side equal to zero and solving for

A k A k,0 = - det (M) (0) ∂ ∂A k [M 11 M 22 ] (0) -∂ ∂A k [M 12 M 21 ] (0) . ( 14
)
We can use the product rule to evaluate the derivatives in the denominator

∂ ∂A k [M ij M pq ] = M ij ∂Mpq ∂A k + M pq ∂Mij ∂A k (15)

Analytical condition for invertibility

Using a linear approximation, the value of the diagonal product term,

[M 11 M 22 ], at A k,0 is [M 11 M 22 ] lin (A k,0 ) ≈ M 11 M 22 (0) + ∂ ∂A 1 [M 11 M 22 ] (0)A k,0 . (16)
Note that the diagonal and off diagonal product terms are equal at the det(M)=0 points. In Sec. 2.6, I show numerically that the linear approximation to [M 11 M 22 ] lin (A k,0 ) is positive for the non-invertible spectrum and negative for at least one of the A k axes for the invertible spectra.

We can derive an analytical formula for [M 11 M 22 ] lin (A k,0 ) by substituting A k,0 from Eq. 14 into Eq. 16 and simplifying:

[M 11 M 22 ] lin (A k,0 ) = . . . M 11 M 22 M 12 M 21 f 2 f k 1 M 12 - f 1 f k 1 M 11 + f 1 f k 2 M 21 - f 2 f k 2 M 22 det(M)(0)(M 1k +M 2k )+M 21 f 2 f k 1 -M 22 f 1 f k 1 +M 12 f 1 f k 2 -M 11 f 2 f k 2 (17)
The simulations also show that the denominator of Eq. 17 is positive for the three spectra tested. Since we are only interested in whether the product value is negative or zero and the simulations show the denominator is always positive, we need consider the sign of the numerator of Eq. 17

N umerator = M 11 M 22 M 12 M 21   f 2 f k 1 M 12 - f 1 f k 1 M 11 + f 1 f k 2 M 21 - f 2 f k 2 M 22   . (18)
The elements of M are effective values of the basis functions and are always positive so the sign is determined by the expression in brackets of Eq. 18.

Notice that for delta function spectra the first and second and the third and fourth terms of the numerator cancel it is equal to zero. Similarly, the denominator is also zero for delta functions so the product term is indeterminate 0 0 and the system will be invertible if the Jacobian determinant at zero thickness is not equal to zero.

Simulations of det(M(A))

In the following sections, I will study by computer simulation whether the Jacobian determinant is zero inside C for several spectra: (1) a set that are known to be non-invertible for two specific A-vectors, (2) a similar set that are better conditioned, and (3) x-ray tube spectra with different voltages.

3 energy non-invertible

Levine [START_REF] Levine | Nonuniqueness in dual-energy CT[END_REF] described a set of spectra that he found were not one-to-one. In this section, we analyze these spectra using the methods described in the previous sections. Fig. 2 (a Fig. 4 shows one dimensional plots of the diagonal and off-diagonal determinant terms on the A bone and A water axes. The curved lines are the actual surfaces with a circle at their intersection where det(M ) is zero. Also shown are straight line fits to the curved lines using the slope at the origin. The slope was computed with Eq. 15. The intersection of the fit lines is the diamond shape. A water axis Fig. 4: Plot of the diagonal and off-diagonal surfaces along the A 1 and A 2 axes. The curved lines are the actual surfaces with a circle at their intersection where det(M ) is zero. Also shown are straight line fits to the curved lines using the slope at the origin. The slope was computed with Eq. 15. The intersection of the fit lines is the diamond shape.

Once the zero determinant points are found on the A k axes, we can use them to determine the zero determinant region as the line joining the intersection points. This is the yellow line in part (b) of Fig. 3. The yellow line and the region of large inverse determinant coincide. A water axis Fig. 8: Plots of the diagonal and off-diagonal surfaces along the A 1 and A 2 axes for the voltage switched spectra. The actual product terms, shown by the curved lines, do not intersect. Also shown are the straight lines fit to the derivative at the origin. At their intersection, shown by the diamond icon, the product term values are negative, which is not physically possible.

3 energy invertible

L data for non-invertible spectra

Invert transform with iterative algorithm

The effect of the zero values of the Jacobian determinant on the invertibility of L(A) was tested with an iterative inverse transform algorithm. The algorithm iteratively found the value of A that minimized |L(A) -L input | 2 . It was implemented using the Matlab f minsearch function.

The norm of the errors of the iterative algorithm for A values on the green line in Fig. 10(a) were computed. They were plotted as a function of distance along the line in Fig. 10(b). The det(M ) was zero on the yellow line in Fig. 10(a

Iterative algorithm with Levine ambiguous points

The iterative algorithm was also used to test the ambiguous A vectors noted by Levine [START_REF] Levine | Nonuniqueness in dual-energy CT[END_REF], [2.5, 1] and [4.09, 0.13].

Using the L vector computed for the first A as input to the iterative algorithm, the solution depended on the initial point of the iteration. If the point was near one of the two values, the algorithm converged to that value. If the point was far from both, for example the origin, the algorithm converged to [2.5, 1].

Values of determinant terms with invertible and non-invertible spectra

As discussed in Sec. 2.5, the analytical formula for signs of values of the product terms at the intersections of the line of zero Jacobian values determines whether the spectra are invertible. The formula was tested with the simulation data described in Sec. 2.6 and the results for the product terms for the three spectra tested are shown in Table 1. The theoretical results indicated that the sign of the product terms is determined by the numerator of Eq. 17. Table 2 shows the numerator and denominator of the product terms equation for the tested spectra.

Discussion

The invertibility of spectral measurements is an important issue. Fig. 5 shows that the errors with the iterative algorithm show a sharp increase near the line of non-invertibility. My previous results with a three dimension system [START_REF] Alvarez | Non-invertibility of spectral x-ray photon counting data with pileup[END_REF] showed that the noise increased by a huge factor near the plane of non-invertibility. X-ray system designers need tools to determine whether their system may become non-invertible in some part of the A-vector operating range because the system may become unusable near the region of non-invertibility.

In the results shown here, the non-invertible region is nearly a straight line in the two dimension A-space. In my previous paper [START_REF] Alvarez | Non-invertibility of spectral x-ray photon counting data with pileup[END_REF] the non-invertible region for a three dimension basis set is a plane. The non-invertibility regions for both cases may be explained by the observation that the elements of the M matrix are nearly linear as a function of A. The results in Figures 2,2, and 5 show that the terms in the determinant of M are also close to linear so their intersection, which results in a zero determinant value, is close to a straight line in two dimensions or a plane in three dimensions.

Sections 2.6 and 2.5 derive analytical formulas for the line of non-invertibility in the two dimension A-space using a linear approximation to the terms in the determinant of M. The line is determined by its intersection with the A axes. For the non-invertible Levine spectrum, Fig. 2.6 shows that the actual and linear approximation intersection points are physically feasible. That is, at positive values and with positive values for the determinant terms at the intersection points. For the invertible three spike spectra, Fig. 2.6 shows that the actual determinant terms do not intersect. The linear approximation does intersect but at the intersection points the determinant product terms are negative which is not physically feasible. For the two voltage x-ray tube spectra, Fig. 2.6 shows that the actual determinant terms again do not intersect and, with the linear approximation, the determinant terms are again negative at the points of intersection.

Table 1 gives the values for the product terms at the intersections with the axes for the three spectra tested. As indicated by the theory, they are positive for the non-invertible spectra and at least one is negative for the three energy good condition and voltage switched spectra, which are both invertible. Table 2 shows that the denominator of Eq. 17 is always positive so the sign of the result for the product term is determined by the numerator and in particular by the expression in brackets in Eq. 18.

Therefore, for the cases tested, the linear approximation gives useful results and the formulas give insight into the conditions leading to non-invertibility. There may, of course, be cases where the linear approximation does not give accurate results. The approach taken was a trade-off between accuracy and simplicity of the formulas. Simple formulas allow us to gain physical insight into the invertibility conditions at the expense of lower accuracy.

Levine [START_REF] Levine | Nonuniqueness in dual-energy CT[END_REF] introduces an alternate approach to the dual energy inversion problem to that described here. The formulation focuses on the exponent of the transmission equation, Eq. 2. It puts the logarithm of the spectral density in the exponent and combines them with A • f (E). The combination is used to derive A-space surfaces specified by the measurement spectra. The A vectors that invert the equations are at the intersection of the surfaces. With a particular set of spectra, Levine found two A-vector points where the surfaces intersect. The transformation is therefore not one-to-one and is not invertible. Interestingly, the two ambiguous points are not on the line of zero determinant.

The results in Fig. 9 show that the L data on the line of non-invertibility are not constant. Instead, the L data vary but the two product terms of the determinant balance out resulting in a zero determinant value on this line.

Research is on-going into more accurate theoretical equations for invertibility and ways to apply the formulation described in this paper to clinical systems design.

Conclusion

An approach to studying the invertibility of spectral measurements is presented and applied to invertible and noninvertible spectra. Non-invertible dual energy spectra have zero Jacobian determinant in the first quadrant of the A vector plane. For the non-invertible spectra tested, the Jacobian determinant is equal to zero on a curve in the A vector plane that is close to a straight line. An analytical formula is derived for a linear approximation to the points where the line crosses the A-axes. If the spectra are invertible, the intersection points are not physically feasible while for the non-invertible spectra they are feasible indicating zero Jacobian determinant values in the first quadrant.

Fig. 2 :

 2 Fig. 2: Three dimension plot (a) of the products of the diagonal and off-diagonal elements of the determinant det(M (A)) = M 11 M 22 -M 12 M 21 as a function of A for the Levine non-invertible three energy spectrum. The determinant is zero if the surfaces intersect. Each product term is plotted as a separate color-coded surface as shown in the legend. Part (b) is a plot of the spectra as a function of x-ray energy.

Fig. 3 Fig. 3 :

 33 Fig. 3 shows views of the diagonal and off-diagonal M elements product data projected onto the A-plane. Part (a) is a view of the diagonal and off-diagonal determinant terms in Fig. 2 projected along the Z-axis. Part (b) is an image where the grayscale is proportional to 1 /|det(M(A))|. The bright line in the image shows the zero values. A straight line was fit to the line of maxima and is plotted on the image in yellow. The close fit of the straight line shows that the surfaces are close to planes so their intersection is a line. Fig. 3 also shows the ambiguous A-vectors found by Levine plotted as the green crosses. Notice that the points do not fall on the line of zero determinant values.

Fig. 5 Fig. 5 :

 55 Fig.5is a three dimension plot of the products of the diagonal and off-diagonal elements of det(M) as a function of A for well-conditioned three energy spectra. Part (b) is a plot of the spectra. Each product term is plotted as a separate color-coded surface as shown in the legend. Notice that the surfaces do not intersect so the determinant is not equal to zero in the first quadrant of the A-plane.

Fig. 6 Fig. 6 :Fig. 7 :

 667 Fig.6shows plots of the diagonal and off-diagonal surfaces along the A 1 and A 2 axes for invertible 3 energy spectra. Notice that the actual product terms, shown by the curved lines, do not intersect. Also shown are the straight lines fit to the derivative at the origin. Their intersection, shown by the diamond icon, occurs outside the region of interest and the product term values are equal to zero, which is mathematically not possible.

Fig. 8

 8 Fig.8shows plots of the diagonal and off-diagonal surfaces along the A 1 and A 2 axes for the voltage switched spectra. The actual product terms, shown by the curved lines, do not intersect. Also shown are the straight lines fit to the derivative at the origin. At their intersection, shown by the diamond icon, the product term values are negative, which is not physically possible.

Fig. 9 Fig. 9 :

 99 Fig. 9 shows the L data for non-invertible 3 energy spectra. Part (a) shows images where the gray scale is proportional to the L components. The images show the line where the determinant is zero and several lines parallel to it. Part (b) has plots of the L data on the lines. The values on the central point of the zero determinant line are subtracted for ease of display. Notice that the L data are not constant on the det(M(A)) = 0 line.

  ) and the large errors in Fig.10(b) are near the point of intersection of the green and yellow lines.

Fig. 10 :

 10 Fig. 10: Effect of det(M(A)) on errors with iterative inverse transform algorithm. The norm of the errors of the iterative algorithm for A values on the green line in panel (a) were computed. They were plotted as a function of distance along the line in panel (b). The det(M ) was zero on the yellow line in panel (a) and the large errors in (b) are near the point of intersection of the green and yellow lines.

Table 1 :

 1 Determinant product terms at intersections of zero Jacobian line with A 1 and A 2 axes.

	A 1	A 2

Table 2 :

 2 Numerator and denominator of Eq. 17 for the product terms of the tested spectra.

	Numerator M 11 M 22	A 1	A 2
	3 Energy Non-invertible	0.00035307 8.1111e-05
	3 Energy Good Condition 2.8091e-05 -4.3262e-06
	Tube Voltage Switched	-0.0055665 -0.0010313
	Denominator M 11 M 22	A 1	A 2
	3 Energy Non-invertible	0.0063432 0.0012686
	3 Energy Good Condition 0.020393 0.0036596
	Tube Voltage Switched	0.071413	0.012247


Léo Le Taro

Hervé Rivano

Simulating Collaborative Sensor Calibration: Convergence and Cost

Keywords: Wireless Sensor Networks, Air Pollution Monitoring, Distributed Algorithms, Sensor Calibration

Air pollution is an increasingly concerning issue in urban areas because of its impact on citizens' health. To tackle pollution effectively, accurate monitoring is a must. Precise stations managed by governmental or specialised institutions do exist, but they are both costly and bulky, which limits the potential to deploy them densely. However, recent progress in micro, connected sensors brings new alternative deployment schemes for dense monitoring by low-cost, imprecise sensors. For such a deployment to be relevant relative to urban air quality monitoring aspects, we are concerned with maintaining the system's properties over time. Indeed, one of the major drawbacks of cheap sensors is their drift: chemical properties degrade over time and alter the measurement accuracy. We challenge this issue by designing distributed, online recalibration procedures. We present a simulation framework modelling a mobile wireless sensor network (WSN) and we assess the system's measurement confidence using trust propagation paradigms. As WSN calibrations translate to information exchange between sensors, we also study means of limiting the number of such transmissions by skipping the calibrations deemed least profitable to the system.

I. INTRODUCTION

Urban pollution monitoring is traditionally carried out using governmental air quality stations, precise but costly and sparsely deployed. In large cities, this limits monitoring resolution to neighbourhood level, actual estimation being computed by numerical models. Recent advances in nanotechnologies is giving birth to small, affordable electrochemical sensors. One of the major research challenges is to achieve a higher spatial resolution with dense and mobile deployments of such lowcost sensors. To keep an exploitable measurement accuracy, sensor calibration must be considered [START_REF] Tolle | A macroscope in the redwoods[END_REF]. However, lowcost electrochemical sensors such as ones we use to measure NO 2 concentrations degrade over time [START_REF] Kumar | The rise of low-cost sensing for managing air pollution in cities[END_REF], hence recalibration is required if the system is to remain usable over a long deployment. Recalibrating nodes in-place is called non-blind calibration [START_REF] Ramanathan | Rapid deployment with confidence: Calibration and fault detection in environmental sensor networks[END_REF]. Ground truths, i.e., reference measurements, are needed to adjust the gain and offset of the low-cost nodes. However, such an approach does not scale because it is infeasible to move a high-quality reference to periodically visit hundreds of sensors.

Another paradigm, denoted blind calibration, assumes unknown ground truths and develops techniques to predict and compensate errors in measurements. Blind calibration methods rely on the underlying signal of interest being either bandlimited (i.e., varying smoothly in time and space) [START_REF] Bychkovskiy | A collaborative approach to in-place sensor calibration[END_REF][5] [START_REF] Lipor | Robust blind calibration via total least squares[END_REF] or sparse [START_REF] Ling | Self-calibration and biconvex compressive sensing[END_REF]. Yet NO 2 fields prove to be neither, exhibiting large spatiotemporal variations [START_REF] Moltchanov | On the feasibility of measuring urban air pollution by wireless distributed sensor networks[END_REF], thus negating the possibility to exploit such properties.

Mobile sensing is gaining more attention as recent studies found that a few mobile nodes on well-selected routes can reflect data as accurately as many static ones [START_REF] Calafate | On the use of mobile sensors for estimating city-wide pollution levels[END_REF]. On top of it, mobile sensor networks offer opportunities for multihop calibration, i.e., freshly calibrated sensors may in turn calibrate others. Work has been conducted to minimise error propagation in multihop calibration [START_REF] Saukh | Reducing multi-hop calibration errors in large-scale mobile sensor networks[END_REF].

II. PROBLEM STATEMENT

Given a heterogeneous set of sensing units, comprised of precise base stations and low quality sensors, assuming all low quality sensors are initially uncalibrated, we are wondering whether the system converges to an exploitable state, if so, we are concerned with the time required to reach a permanent regime.

Secondly, we are interested in limiting the number of energy-hungry transmissions between sensors, and wonder how saving data exchanges would affect the system's accuracy.

In Section III, we present the model upon which we based our simulation framework. In Section IV, we detail the algorithm that the framework implements. Section V presents mathematical analyses of binary calibration in our model. Simulation results are laid out in Section VI. We then conclude the paper in Section VII, suggesting future work ideas as well.

III. MODEL

We model our system as a discrete time and space process in which sensors follow a stochastic mobility pattern. There are S = {1..s} initially uncalibrated mobile sensors randomly moving within a space of P = {1..p} positions. For convenience, mobile sensors move following a uniform distribution: at any given time, a sensor moves to x ∈ P with probability 1/p. Several sensors may share the same position.

To simulate the presence of precise institutional or governmental air quality stations in a real-life urban scenario, let R ⊂ P be the subset of r positions featuring a static reference station, assumed perfectly calibrated and reliable.

Each mobile sensor a is assigned a trust value C a ∈ [0..1], 0 meaning completely inaccurate and 1 meaning completely accurate. Unless recalibrated, mobile sensors degrade and lower their trust following an exponential decay of rate γ ∈

[0..1]: C a (t) = C a (0) • e -γt .
Reference stations retain a constant trust of 1. Initially, all mobile sensors start with a trust of 0.

Mobile sensors are said to have a rendez-vous with a reference station or another mobile sensor when, at the same time instant, they stand at the same position as, respectively, a reference station or another mobile sensor.

Finally, we introduce the recalibration threshold q which conditions calibrations between mobile sensors. The rationale behind this parameter is detailed in Section IV.

IV. IMPLEMENTATION

Each sensor executes the following pseudo-code algorithm at every time step: We can observe that mobile sensors unconditionally recalibrate to reference stations, however calibration to a peer is conditioned by the variable q. In the real world, each calibration means wirelessly exchanging information. We wish to carry out only meaningful calibrations, with a high benefit. This threshold q allows skipping the calibrations deemed not worthwhile, as they would not lead to a noticeable increase of the trust.

Data: γ, q C ← C • e -γ ; foreach n in find_neighbours(p, t) do if n is a reference station then C ← 1; else if n is a mobile sensor and C + q < n.C then C ← n.C; end end
As a corollary, a q of 1 means mobile sensors do not cooperate between each other and only recalibrate to reference stations. Conversely, a q of 0 means that a sensor will always upgrade its trust to the maximum of its neighbors'.

V. BINARY CALIBRATION: ANALYTIC EVALUATION

Binary calibration is the concept of considering either fully uncalibrated or fully calibrated sensors, with trust values being respectively 0 or 1. In our framework, binary calibration is equivalent to letting γ be zero. Let the random variable T represent the date of a sensor's first calibration.

A. Sensors do not cooperate

Each sensor has probability r/p at each time step to encounter a reference station, independently of other sensors. The probability for such a sensor to remain uncalibrated at t follows a geometric progression of rate 1 -r/p: P (T > t) = (1 -r/p) t .

The system's global average trust can be derived:

C(t) = 1 -P (T > t) = 1 -(1 -r/p) t

B. Sensors cooperate

Being calibrated at t means:

• either having calibrated between 1 and t -1 (noted T ≤ t -1); • or not having calibrated between 1 and t -1 but having a rendez-vous at exactly t (noted T = t).

Both cases are mutually exclusive, let us sum their probabilities, keeping in mind sensors start the simulation uncalibrated (time step 0):

P (T = 0) = 0 P (T = t) = 1 -(1 - 1 p ) r+sP (T ≤t-1) P (T ≤ t) = P (T ≤ t -1) + (1 -P (T ≤ t -1))P (T = t)
The rest of this paper presents results of our simulation framework whose purposes are to validate our binary calibration analysis, as well as study cases of non-binary calibration, with γ > 0.

VI. NON-BINARY CALIBRATION: PERFORMANCE EVALUATION

In our system, individual sensors' trust occasionally jump when they recalibrate, then continuously variate to lower values because of decay. This behaviour is illustrated by the results of a simple simulation plotted in Figure 2 with a reduced number of sensors and grid positions to preserve readability. This figure hints at the possibility of a permanent regime, because after an initial period of 250 steps no sensor seems to fall below a trust of 0.8. How does the average trust of all sensors behave over a longer period of time and with different parameters?

Figure 3 shows that after a short hysteresis-shaped transient phase, the system's average trust converges to a permanent regime where the trust remains stable. Almost all sensors' trusts evolve between 0.5 and 0.95 despite a pessimistic decay rate γ = 10 -3 . Cooperative recalibration is therefore able to maintain the accuracy of the system.

The existence of a permanent regime raises the question of how much time is necessary to reach it. We define the time to converge t 0.9 as the time by which the average trust of the system C(t 0.9 ) exceeds 0.9. t 0.9 is plotted against r in Figure 4, which validates our framework against the analytic evaluations conducted in in Section V for binary calibration schemes (γ = 0). It shows that cooperation between sensors have a much stronger impact on the convergence speed than the number of reference stations.

We illustrate the impact of the decay γ on the system's average trust, as well as the maximal and minimal trusts achieved by non-reference sensors in Figure 5. We also consider the ratio of sensors above a trust of 0.5.

We can observe than for γ = 0.002, the average trust is close to 0.5, which is also the median trust with 50 percent of sensors lying above that value. Like all others, this simulation was run 100 times and results were aggregated to mitigate the stochastic nature of the results. Nevertheless, we notice that the strongest sensor's trust plot (in light blue) has a quite random aspect, explained by the fact that it depends a lot on its "luck" with rendez-vous rather than exclusively the values of simulation parameters. We deduce that while the average trust we can expect from our system is quite predictable, we cannot guarantee the reliability of a given sensor. Hence, in a real-world urban situation, areas where monitoring is critical (e.g., school or retirement home neighbourhoods) should be covered by many fixed sensors, and/or reference stations. Assuming sensors exchange information once for each recalibration, neglecting the cost of polling one's neighbour's trust, Figures 6 and7 depict the trade-off between the number of wireless transmissions and the accuracy of the system, i.e., the average trust. The recalibration threshold q is the parameter that restricts transmissions to the case of "useful" recalibrations. Numbers were gathered during the permanent regime. An increase of q degrades the system's trust as expected. However, using even a low q dramatically limits the number of wireless transmissions while sacrificing very little trust.

VII. CONCLUSION AND PERSPECTIVES

In this paper, we have presented a simulation framework capable of providing insight on what we can expect from autocalibrating low-cost electrochemical sensor networks. Cooperation between sensors make the system converge quickly and maintain a steady average accuracy despite the individual decay of the sensors, without requiring many reference stations. It is also shown that many wireless transmissions are spent on a very marginal improvement of the accuracy.

Adjusting certain parameters, we were able to make optimistic predictions of certain metrics like the measurement confidence of the system, its convergence time and the number of transmissions required to calibrate nodes.

Current and future work shall integrate more realistic mobility patterns into the model. One could expect that an urban mobility pattern induce correlation between subsets of sensors, hence a less stable "permanent" regime. The impact of the mobility model and the deployment of reference stations on the spatial mapping of trust values should then be investigated. Besides, realistic mobility models should take into account the distance between calibrating and calibrated nodes into the trust propagation function. Finally, calibration theory requires a non-correlated set of simultaneous calibrating and calibrated measurements. Such a set could be collected over a sequence of rendez-vous, at the cost of a more complex interpretation of the trust.

Figure 1 .

 1 Figure 1. Sensor trust updating process

Figure 2 .Figure 3 .

 23 Figure 2. Sensor trust overview (p = 400, s = 15, r = 1, γ = 10 -3 , q = 0)

Figure 4 .

 4 Figure 4. Time to converge (p = 10000, s = 150, γ = 0)

Figure 5 .

 5 Figure 5. System state versus Gamma (p = 10000, s = 150, r = 1, q = 0)

6 2.5x10 6 Figure 6 .Figure 7 .

 6667 Figure 6. Tradeoff between trust and transmission cost (p = 10000, s = 150, r = 1)

ACKNOWLEDGEMENT

This work has been supported by the "LABEX IMU" (ANR-10-LABX-0088) and the "Programme Avenir Lyon Saint-Etienne" of Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR).