domain: info math

text:
Matthieu Dolbeault

Olga Mula

Agustín Somacal

State estimation of urban air pollution with statistical, physical, and super-learning graph models

We consider the problem of real-time reconstruction of urban air pollution maps. The task is challenging due to the heterogeneous sources of available data, the scarcity of direct measurements, the presence of noise, and the large surfaces that need to be considered. In this work, we introduce different reconstruction methods based on posing the problem on city graphs. Our strategies can be classified as fully data-driven, physics-driven, or hybrid, and we combine them with super-learning models. The performance of the methods is tested in the case of the inner city of Paris, France.

Introduction

Background and motivation

Data-driven estimations are becoming increasingly relevant and widespread as the volume and heterogeneity of available data increases. A fundamental challenge is to build numerical methods for which one can estimate how optimally they exploit the given information. The present paper addresses some essential computational aspects connected to this question. More specifically, our goal is to reconstruct a state u of a physical process, for which we have at hand very heterogeneous sources of data coming from direct partial observations of u, from quantities related to u, and from the knowledge that the physics can be modelled by a Partial Differential Equations (PDE).

Assume that u belongs to some Banach space U of potentially infinite dimension, with associated norm • U , and that all the available information is given by an element x u from some abstract metric space X . Our goal is thus to build a mapping A : X → U such that A(x u ) approximates u at best, in the sense that the approximation error e(A, u) = u -A(x u ) U (1.1) is as small as possible, for any configuration (u, x u ) of the system. In practice, finding the optimal map A is not feasible, and various suboptimal reconstruction techniques have been proposed, each of them having its own virtues and drawbacks: statistical approaches such as BLUE [START_REF] Thomas | The design and analysis of computer experiments[END_REF] and kriging [START_REF] Cressie | Statistics for spatial data[END_REF], model order reduction of parametric PDEs [START_REF] Hesthaven | Certified reduced basis methods for parametrized partial differential equations[END_REF][START_REF] Mula | Inverse problems: A deterministic approach using physics-based reduced models[END_REF], or more recently approximations by neural networks and machine learning strategies [START_REF] Molinaro | Neural inverse operators for solving PDE inverse problems[END_REF][START_REF] Dahmen | Nonlinear reduced DNN models for state estimation[END_REF]. Since all of these strategies are sub-optimal, and each one is based on different a priori assumptions, one should not make the methods compete against each other, but rather collaborate with each other to enhance their respective strengths. This leads naturally to explore approaches based on ensemble super-learning [START_REF] Breiman | Stacked regressions[END_REF][START_REF] Van Der Laan | Super learner. Statistical applications in genetics and molecular biology[END_REF][START_REF] Polley | Super learner in prediction[END_REF] as we consider in the present work.

Urban air pollution modeling

There are numerous applications in which one is confronted with the above state estimation problem. As a guiding example, we consider in this paper the real-time reconstruction of urban pollution fields. Beyond the relevance of such a task to limit environmental and health risks in the city, pollution state estimation is an excellent example where collaborative, super-learning methods are required. This is because the problem accumulates several difficulties that make the reconstruction challenging for most common reconstruction methods. Among the issues, we may mention the following:

• Scarcity of pollution measurements: The amount of reliable sensor devices measuring pollutant concentrations is often limited, and the measurements are usually taken at fixed locations. As a result, reconstruction methods based solely on these measurements lack spatial resolution, and exhibit huge uncertainties in regions without sensors.

• Heterogeneous data: In addition to the pollutant measurements, other sources of relevant information are available such as traffic estimations in each street, wind speed, topography, temperature, etc. However, it is not obvious how to meaningfully combine this data to enhance the estimation. Some attempts have been tried in [START_REF] Criado | Data fusion uncertainty-enabled methods to map street-scale hourly NO 2 in Barcelona: a case study with CALIOPE-Urban v1.0[END_REF] through the use of Gradient Boosting Machines and Universal Krigging, with positive results for the estimation of PM10 particles in the city of Barcelona.

• Lack of training data: Even when incorporating other sources of information, the available data may be insufficient, noisy or hardly correlated to the pollution levels we wish to caracterize. Purely data-driven models greatly suffer from these impediments in their training phase.

• Complexity of the physical problem: The equations governing the dispersion of pollutants in the atmosphere are nonlinear, with turbulent effects at the street scale, thus imposing a fine spatial resolution, at least near the sensor stations [START_REF] Criado | Data fusion uncertainty-enabled methods to map street-scale hourly NO 2 in Barcelona: a case study with CALIOPE-Urban v1.0[END_REF]. On the other hand, the computational domain is of the size of a city, making it prohibitively expensive to solve a full model like 3D Navier-Stokes equations.

• Parameter calibration: Reduced models use effective parameters, which account for largescale averages of local effects, in order to alleviate the requirements on the resolution. However these parameters must be calibrated based on the available data or preliminary simulations, which is a hard task given the above issues.

The above obstructions advocate for collaborative strategies combining physics-driven and datadriven approaches such as the one that we develop in this paper. A similar idea has been explored in [START_REF] Niu | A novel hybrid decompositionand-ensemble model based on CEEMD and GWO for short-term PM2.5 concentration forecasting[END_REF], but for forecasting temporal series of pollutant, instead of performing state estimation on a large spatial domain.

It should be noted that the limited number of reliable measurements will still pose problems for validating and assessing the quality of each model, which is a crucial part in collaborative strategies. We will mitigate this defect by operating multiple leave-one out cross validations, which preserve as much data as possible for the training part of each model, while testing them on many instances.

Contributions and layout of the paper

Our main contributions are:

1. the construction of numerous physics-based and data-driven models for state estimation; 2. the construction of a very general ensemble super-learning method combining the above models;

3. its application to the task of recovering urban pollution maps at a city scale, together with a comparison of its constitutive submodels;

4. and the development of a routine extracting car emissions in each street from traffic maps. Moreover, we have created a dataset comprised of processed traffic data from Google Maps screenshots, which can be used for future research.

In our numerical experiments, we work with the inner city of Paris, which covers a surface of about 140 km2 . The pollutant we consider is NO 2 , which is monitored for its respiratory effects, while being mainly produced by vehicle emissions. We use concentration measurements from Airparif sensors1 , and real-time traffic data from Google Maps 2 . Compared to previous contributions and other existing reconstruction methods (see, e.g., [START_REF] Mallet | Ozone ensemble forecast with machine learning algorithms[END_REF][START_REF] Tilloy | BLUE-based NO 2 data assimilation at urban scale[END_REF]), the use of such online traffic data is rather novel. It gives a rough estimation of the spatial density of street traffic, benefits from a very fine spatial resolution, and can be freely updated as frequently as desired, in contrast to many existing approaches which only use time averages of traffic data.

Another distinctive aspect of our approach is the representation of the city by a graph, where nodes and edges correspond to crossroads and street segments, instead of considering an open subset of R 2 or R 3 as the spatial domain. This description immediately includes geometric specificities of the agglomeration under study, such as the orientation of each street or the configuration of each neighborhood. It is a natural framework for taking into account pollutant emissions caused by traffic, which are located on the graph. Moreover, it is in adequacy with our goal of estimating local variations in the concentration of pollutants close to the ground, since the streets are isolated from each other at this height.

Physical models can be solved on such domains thanks to the theory of quantum graphs, that is, metric graphs endowed with a differential operator acting on functions defined on the graph (see [START_REF] Berkolaiko | Introduction to quantum graphs[END_REF] for details and references). The metric graph structure leads to the definition of suitable and natural function spaces to pose the problem. Of course, several physical models of different complexity could be considered. In this paper, we work with simple elliptic operators, but the model could be refined by considering, for instance, advection-diffusion operators.

The rest of the paper is organized as follows. In Section 2 we present our guiding numerical example of the Parisian area, and the available data. Section 3 explains the different reconstruction methods we have used for our numerical experiments, including ensemble super-learning methods combining the previous ones. By construction, the super-learner has higher approximation power than each individual model. Section 4 discusses how to theoretically quantify performance and optimality of the numerical algorithms, and why leave-one-out is a good way to estimate this performance in practice. We summarize our numerical experiments and provide some illustrations in Section 5. Finally, Appendix A details the mathematical setting for the problem of pollution state estimation on graphs. 2 Available data and pre-processing

Pollution sensors

The main information we use consists of direct measurements of the NO 2 concentration field u at Airparif sensor stations. There are m = 13 such stations, which are placed at fixed locations

r obs := {r obs 1 , . . . , r obs m } ∈ (R 2 ) m ,
see Figure 1. Each of the stations provides hourly averages of the concentration of nitrogen dioxide, in µg/m3 , of the form

z i = u(r obs i ) + η i , i = 1, . . . , m,
where η i is some noise with nominal relative error |η i |/u(r obs i ) ≤ 15%. Note that, in principle, the equations governing the dispersion of pollutants are time dependent. However, as we only have measurements every hour, we opt for a static model, where the state at a given time is computed based on the data available at this time only. This essentially amounts to assuming that the emissions vary slowly over time, and that the system reaches an equilibrium state in less than one hour.

Meteorological conditions

Wind, as well as stratification effects in the atmosphere due to variations in temperature, play a major role in the dispersion of pollutants [START_REF] Cai | The effects of thermal stratification on airborne transport within the urban roughness sublayer[END_REF]. Moreover, the chemical equilibrium between NO and NO 2 depends on the cloud cover [START_REF] Li | Quantifying the nitrogen isotope effects during photochemical equilibrium between NO and NO 2 : implications for δ 15 N in tropospheric reactive nitrogen[END_REF]. Therefore, we collect the temperature θ ∈ R and the wind speed w ∈ R 2 at every hour from a weather archive 3 . These two quantities are treated as global, that is, they are assumed to be constant over the spatial domain.

Traffic

Car traffic is responsible for more than half the emissions of NO 2 in urban environments [START_REF] Kurtenbach | Primary NO 2 emissions and their impact on air quality in traffic environments in germany[END_REF]. There is an increasing number of available sources that give access to traffic data. In our case, Figure 2: Raw data from Google Maps: the image contains the city with its main landmarks, and some streets are highlighted with one of the four colors corresponding to traffic.

we work with traffic information extracted from Google Maps. We have designed a script using the Python library Selenium to automatically take screenshots of Paris every 15 minutes over an area of 1253 × 1253 pixels with zoom level 13. An example of resulting raw image can be seen in Figure 2. Note that city landmarks could not be removed before taking the screenshot, nor even by substracting a background image, since each screenshot has slight color variations, rendering this approach impractical. Another issue is the absence of traffic data in the smallest streets of the city. In addition, linking it to the pollution field requires some calibration.

On the one hand, this kind of information is very rich because of its availability in real time, and its spatial coverage of the whole city at high resolution. On the other hand, it is very partial: it comes in the form of four colors, each representing a certain traffic intensity, which gives a qualitative estimate of the number of cars in each street: a street marked in red is, for instance, more congested than one marked in green. One main goal in our work is to examine the potential of incorporating such low-quality information for state estimation tasks.

Graph of Paris

In order to locate our sources and to express the spatial dependence of a state, we consider a graph domain. We use a metric graph G = (V, E) provided by Open Street Maps, together with the Python library osmnx. For the mathematical definition of metric graphs and their associated function spaces, we refer to Appendix A. The graph G covers the whole inner ring of the city, as shown in Figure 3. The full graph has |V| = 12963 vertices and |E| = 25476 edges, but we restrict it to the biggest connected component of the subgraph that remains after filtering out all the edges which have never been colored with traffic information. After this operation, our actual graph has |V| = 10116 vertices and |E| = 18713 edges. The street network is relatively dense, most nodes having 3 to 6 edges.

The vertices v ∈ V come with precise geographical coordinates. In the following, we assume that the graph is embedded in the two-dimensional plane, and do not take altitude into account. Each edge e ∈ E is a street or a portion of it, and we have access to its length e as well as its shape, number of lanes and speed limit. The information is so detailed that the streets are represented by paths that are not necessarily straight lines. However, in the following, we will work under this slightly simplified setting.

The location of the sensor stations does not exactly match with vertices of the graph. We therefore project their position r obs onto the nearest vertex, which yields observational nodes

v obs i := argmin v∈V |r obs i -v|, i = 1, . . . , m,
As Figure 1 illustrates, the projected locations are very close to the exact locations, with a maximal discrepancy of 165m, to be compared with the width of the domain, of about 12km. As a consequence, we will assume that the observations z i correspond to the values u(v obs i ), up to a slight increase in the measurement errors |η i |.

Pre-processing of traffic data

We also map the traffic information onto pollutant emissions on the graph edges, by implementing the following pipeline:

• Cropping: Starting from a raw image like Figure 2, we first crop it to the shape 800×1000, in order to eliminate toolbars and adapt it to the size of the graph. The background of Figure 1 is obtained by the same procedure.

• Traffic colors extraction: The colors associated with the four levels of traffic colors := {green, orange, red, dark-red} are easily identified 4 . They seem to be used exclusively for that purpose, hence it suffices to extract the pixels having one of these colors.

• Projection on graph edges: These pixels, once expressed in their geographical coordinates, almost perfectly overlap the metric graph from Open Street Maps. For each edge e ∈ E and each color c ∈ colors, we count the number p e c of pixels of color c that are closest to edge e. Note that the traffic color might change along an edge, in which case we give up on some local information by only considering the total traffic on the edge.

• Edge normalization: We then transform these pixel counts into proportions of traffic colors on each edge. As the edges may remain blank at times where there is no traffic, we take as a normalizing constant the maximal amount of pixels encountered over all times T for which we collect traffic data:

q e c (t) = p e c (t) max t ∈T c ∈colors p e c (t )
, t ∈ T.

In this way, q e c (t) ∈ [0, 1] indicates the proportion of edge e colored with c at time t, but remains null if no traffic is reported.

• Hourly averaging: As the pollution information is only available every hour, we take the average of the four values of q e c encountered every fifteen minutes, which we still denote q e c in the sequel.

• Projection on graph nodes: In our models, it is in fact simpler to localize emissions on the nodes of the graph. For this reason, we calculate the density of each traffic color c around a vertex v ∈ V as a weighted average on its neighboring edges E(v)

q v c (t) = e∈E(v) a e q e c (t) e∈E(v) a e ,
where a e stands for the area of the road associated to edge e, given by the product of its length e with the number of lanes.

Summary

While the history of sensor and weather data can be found on archives, our script for capturing traffic images only runs in real time, since Google Traffic only provides current information. We collected all types of data on an hourly basis for a period of time comprised between December 9, 2022 and March 19, 2023. After removing time stamps for which some data was missing, we end up with a set of acquisition times T , of cardinality |T | = 1712, which we divide into a set T train of 1338 training times, and a set T test of 374 testing times.

In the end, given the graph G = (V, E), the available information at any time t ∈ T is of the form

x = (v obs , z, θ, w, (q v c ) c,v ) ∈ X = V m × R m × R × R 2 × R 4|V| . (2.1)
In the next section, we present various models to estimate the pollution field from this data, using either statistical inference, linear mappings based on expert knowledge, or neural networks.

Reconstruction methods

We have implemented several methods of state estimation by leveraging the different information sources. Our methods give reconstructions on the metric graph G, that is, we consider mappings A : X → U where U is a space of functions defined on G. Typical examples are U = C(G), L 2 (G) or H 1 (G), as defined in Appendix A. As our main interest is in assessing the effect of incorporating indirect information like the real-time traffic data, we first consider models that take only a portion of x ∈ X as input.

Spatial average

If we give up on all the spatially-dependent data v obs and (q e c ) c,e , the reconstruction is necessarily constant over the whole domain G, which yields no better choice than the average of the observed concentration values

A avg (x)(r) = z = 1 m m i=1 z i , ∀r ∈ G.
This extremely simple reconstruction will serve as our baseline to compare more sophisticated reconstructions. In the sequel, we will add the spatially-dependent data and view the other models as corrections to the spatial average above. This will result in spatially unbiased estimators, provided that the locations of the sensors are representative of the whole pollution field. More precisely, assuming that the stations are randomly drawn according to the uniform probability distribution µ on G, the expectation over z 1 , . . . , z m of the spatially-averaged error is

E z G u(r) -A avg (x)(r) dµ(r) = G u dµ -E z 1 m m i=1 z i = 0,
and this remains true when adding to A avg a correction of vanishing expectation.

Best unbiased linear estimator

If we only want to estimate a missing measurement z i at a given station i ∈ {1, . . . , m}, we may also use statistical information stemming from the history (z t i ) of the station at previous times t, as well as the observations from other stations j = i in the present and the past, denoted respectively z j and (z t j ). For T train the set of training times, define the empirical average

z i := 1 |T train | t∈T train z t i and empirical covariance matrix K ∈ R m×m with entries K i,j := z i -z i z j -z j = z i z j -z i z j .
Any unbiased linear estimator zi of z i is of the form

zi = z i + j =i c j z j -z j ,
for some coefficients (c j ) j =i . Let c ∈ R m be the vector with coordinates c j for j = i and c i = -1.

Then the best linear unbiased estimator (BLUE) is obtained by optimizing the averaged squared error argmin

(c j ) j =i zi -z i 2 = argmin (c j ) j =i c Kc = (K j,k ) j,k =i -1 (K j,i ) j =i ,
where (K j,k ) j,k =i and (K j,i ) j =i are seen as a matrix in R (m-1)×(m-1) and a vector in R m-1 .

If the set of training times T train is large enough, we expect an ergodicity property of the form z i ≈ E(z i ) to hold. For this reason, BLUE should be a near minimizer of the expected squared error, given the available data. Therefore, in the numerical experiments, we will evaluate the different methods A by comparing A(x \ {z i })(z i ) and z i , and the error |z i -z i | 2 will act as an optimality benchmark.

It should be emphasized that BLUE itself is not a valid reconstruction method, since it requires statistical information which is accessible only at the locations of the stations v obs i . Hence this estimator cannot be computed at any point r ∈ G of the graph domain.

Kriging

In order to transform BLUE into a reconstruction method, one needs to propose a surrogate for the correlation between any two points in the graph. Moreover, as we don't know the average pollution at all points of the graph, we proceed without subtracting spatial averages z i in this subsection, in contrast to the previous one. Therefore, we consider the Gram matrix of normalized second-order moments

G i,j = z i z j z 2 i z 2 j .
Taking the positions v obs of the stations into account, we observe that each entry G i,j partly depends on the distance |v obs i -v obs j | between the stations, see Figure 4. A typical choice of approximant is the Gaussian kernel Remark 3.1. The fact that C < 1 can be interpreted as the presence of random noise η i on the measurements z i = u(v obs i ) + η i . As a safety check, one may notice that the average relative error

G i,j ≈ Ĝi,j := C exp - |v obs i -v obs j | 2 2σ 2 + (1 -C)δ i,
η 2 i u(v obs i ) 2 = z 2 i u(v obs i ) 2 -1 ≈ 1 C -1 = 3.31%
is effectively much smaller than the uniform error guarantee η i /|u(r obs i ) L ∞ ≤ 15% that we discussed in Section 2.1. Adding the matrix (1 -C)I ensures that Ĝ has ones on its diagonal, as expected of a correlation matrix, and regularizes the system, by making the inversion of Ĝ stable. More practically, the reconstructed value A(x)(v obs i ) will not be exactly z i , but rather an average of the measurements close to v obs i .

Let r ∈ G, we again examine a linear model

A krig (x)(r) = m i=1 c r i z i ,
with coefficients c r ∈ R m to be determined. This estimator is unbiased if and only if m i=1 c i = 1, and in that case we can write it as a correction to the temporal or spatial average

A krig (x)(r) = A krig (x)(r) + m i=1 c r i (z i -z i ) = z + m i=1 c r i - 1 m (z i -z).
By analogy with BLUE, we thus define the weights as the renormalized solution of a system of correlation equation

c r = ĉr m i=1 ĉr i , ĉr = Ĝ-1 g r , g r i = C exp - |v obs i -r| 2 2σ 2 .
Although there are no optimality guarantees, we expect kriging to have intermediate performance when compared to the spatial average baseline, and to the ideal BLUE reconstruction. However, due to the important spacing between peripheral stations, we only observe an improvement in the central part of Paris. The insufficient density of pollution measurements calls for models involving other sources of information, such as traffic data. This is the objective of the next two subsections.

Source model

The simplest way to incorporate traffic data consists in using only local values q e c for estimating the pollution on an edge e ∈ E, or q v c for a node v ∈ V. As we projected the station locations on V, we focus on the latter case here. We call such a method a source model, since it directly maps the sources of emission to pollution values.

We opt for a linear model acting as a correction on the spatial average baseline:

A src (x)(v) = z + c∈colors α c (q v c -q), (3.1)
where we substracted the spatial average of traffic q for unbiasedness. The vector of coefficients α ∈ R 4 is found by solving a LASSO problem

min α∈R 4 t∈T train m i=1 |z t i -A src (x(t))(v obs i )| 2 + λ α 1 ,
and we perform a cross-validation to estimate the optimal parameter λ, in order to prevent overfitting.

Alternatively, we can also write nonlinear variants, of the form

A src (x)(v) = z + T α ((q v c ), θ, w) , (3.2)
which may take into account other sources of information like temperature θ and wind w. Here, T : R #α × R 4 × R × R 2 → R can be a polynomial combination of the inputs ((q v c ), θ, w), a neural network, or any other nonlinear mapping. The set of parameters α is no longer associated to the four traffic colors, but still needs to be learned via a LASSO regression.

All such models rely on the assumption that pollution depends on its sources in a very localized manner. However, the traffic charts (q v c ) v∈V exhibit sharp variations from one node to its neighbors, which incites to smooth the emissions before applying the above methods. This is the purpose of the following section, which attempts to model such a diffusive behavior.

Physical modeling

An important inspiration for reconstruction methods comes from the physical modeling of pollution dispersion. In our setting, we resort to building a quantum graph, that is, we endow our metric graph G with a differential operator acting on functions from functional spaces such as L 2 (G) or H 1 (G), as defined in Appendix A (we also refer to [START_REF] Berkolaiko | Introduction to quantum graphs[END_REF] for more details and references). Our approach can be summarized as follows:

Elliptic equation: One can first model pollution with a time-independent elliptic equation, by assuming that all time-dependent parameters have sufficiently slow variations, here over the course of an hour, for the pollution field to reach a steady state. For any point r ∈ G, the pollutant concentration is modelled by a function u : G → R solution to the diffusion-reaction equation

P(u) := - d dr • a(r) d dr u(t, r) + h(r)u(r) = q(r), r ∈ G (3.3)
which we choose to complement with "Newmann-Kirchoff" conditions on the vertices, that is,

e∈E(v) du dr e (v) = 0, v ∈ V, (3.4)
expressing the conservation of the quantity of pollutant at every crossroad v ∈ V. Here, E(v) denotes the edges having v as an endpoint, and the derivatives are assumed to be taken in the directions away from the vertex.

In equation (3.3), the function a ∈ L ∞ (G) is an effective diffusion coefficient, which takes into account turbulent dissipative effects. The absorption coefficient h ∈ L ∞ (G) models the leakage of pollutants from the streets to the higher atmosphere, as well as chemical reaction, in particular between NO 2 and other nitrogen oxides, which are not measured by the sensors. Lastly, the source term q ∈ L 2 (G) models all possible emissions of pollutant, which in the case of Paris essentially come from traffic, local heating, and industrial and urban activity outside the city. As we only have access to local traffic data, we assume that the other sources are spatially constant, and average them out by solving P(u -ū) = q -q, where ū is the spatial average of the pollutant concentration, estimated by A avg (x) = z, and q -q corresponds to the variations of traffic around its spatial average, computed through the procedure from Sections 2.5 and 3.4. Remark 3.2. Equation (3.4) has a similar effect as Newmann conditions at the borders of the spatial region under consideration, here the rectangle contour of Figure 1. Therefore the only exchanges with the exterior of this region are contained in the source term q. As this no-flux condition only give a very rough approximation of the solution close to the border, in the numerical experiments of Section 5, we will concentrate on the accurate prediction of the pollution in the central part of the city.

Variational formulation: The operator P(u) in (3.3) is defined for functions u ∈ H 2 (G), but the equation can be stated in a weak form, which only requires that u ∈ H 1 (G). Multiplying (3.3) by a sufficiently smooth test function v ∈ H 1 (G), and using the Kirchoff-Neumann boundary conditions, it follows that the corresponding weak formulation of the problem is to find u

∈ H 1 (G) such that b(u, v) = f(v), v ∈ H 1 (G) (3.5)
where b is the symmetric bilinear form defined as

H 1 (G) 2 → R b : (u, v) → e∈E e
a(r) du dr (r) dv dr (r)dr + e h(r)u(r)v(r)dr and f : v ∈ H 1 (G) → e∈E e q(r)v(r)dr is a continuous linear form.

Assuming that a(r) ≥ a 0 > 0 and h(r) ≥ h 0 > 0 for r ∈ G a.e., we see that b is continuous and coercive in H 1 (G) with coercivity constant min(a 0 , h 0 ), and continuity constant max( a L ∞ (G) , h L ∞ (G) ). By the Lax-Milgram theorem, problem (3.5) admits a unique solution u ∈ H 1 (G). Discretization: In our numerical tests, we discretize the equation with P 1 finite elements, that is, continuous functions whose restriction to any edge is affine. We describe below the main guidelines, and refer to [START_REF] Arioli | A finite element method for quantum graphs[END_REF] for further details and a complete analysis.

We define the set of hat functions {ϕ v } v∈V by ϕ v (v ) = δ v,v for any vertices v, v ∈ V, and

∀x e ∈ [0, e ], ϕ v (x e ) = 1 -xe e , if e ∈ E(v), 0, if e ∈ E(v),
for any edge e ∈ E. Fixing our finite element space

P 1 = span{ϕ v } v∈V ⊂ H 1 (G), we search for the Galerkin solution û = v∈V c v ϕ v ∈ P 1 such that b(û, v) = f(v), v ∈ P 1 .
Gathering the expansion coefficients of the solution in the vector c = {c v } v∈V , we obtain the linear system of equations

B c = f (3.6) with B = (b(ϕ v , ϕ v )) v,v ∈V and f = (f(ϕ v )) v∈V .
Again by Lax-Milgram theory, this system is invertible, which allows to compute the solution û.

Reduced models: Unfortunately, solving equation (3.6) is expensive, given the size |V| ≈ 10 4 of the graph, so we cannot afford to find û at each time step. In order to mitigate the computational cost, we rely on model order reduction techniques, which have received much attention in the context of parametrized elliptic PDEs [START_REF] Cohen | Approximation of high-dimensional parametric PDEs[END_REF][START_REF] Cohen | Analytic regularity and polynomial approximation of parametric and stochastic elliptic PDE's[END_REF][START_REF] Eftang | An "hp" certified reduced basis method for parametrized elliptic partial differential equations[END_REF][START_REF] Hesthaven | Certified reduced basis methods for parametrized partial differential equations[END_REF][START_REF] Rozza | Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations[END_REF][START_REF] Zou | An adaptive local reduced basis method for solving PDEs with uncertain inputs and evaluating risk[END_REF]. Here, the parameters would be the diffusion a, the reaction h, and the right-hand side q. We consider three reconstructions methods.

1. Eigenstates of the graph Laplacian: One option consists in taking as a reduced model the subspace V n ⊂ H 1 (G) spanned by the n first eigenfunctions of the Laplacian in P 1 . As this operator is self-adjoint and coercive, it admits a spectral decomposition with positive eigenvalues, and the coefficients of the eigenstates in the basis {ϕ v } v∈V are the eigenvectors of B. Assuming that the diffusion and reaction coefficients a and h are constants calibrated in a pre-processing phase, we define the reconstruction mapping A :

X → H 1 (G) by taking û = A(q) the Galerkin projection of u onto V n , that is, by searching û ∈ V n solution to b(û, v) = f(v), v ∈ V n ,
which is simply a diagonal system in the eigenstate basis. We then plug û instead of q in equation (3.1) or (3.2), and learn the coefficients associated to each color, or the more general parameters α.

Principal components of traffic data:

Starting with the whole history of traffic data (q(t)) t∈T ∈ R |T |×|V| , we can also perform a singular value decomposition to find the n first modes q 1 , . . . , q n , compute the solutions to P(u k ) = q k , and assemble them in a reduced space V n = span{u 1 , . . . , u n }. In this way, we expect V n to better capture physical properties of the pollution field, such as strong correlations along a large avenue.

As full-order solves remain costly, we resort to a convolution with a gaussian kernel:

(u k ) v c = v ∈V e -d 2 vv /2δ 2 (q k ) v c φ v v ∈V e -d 2 vv /2δ 2
where d vv = |v -v | is the distance in R 2 (which is equivalent, up to constants, to the distance on the graph). We set δ = 400 m, after observing that pollution data is optimally correlated to regularized traffic information for δ close to this value.

In our experiments, we perform the smoothed projections q → û = n k=1 q, q k u k into a different reduced space for each of the four traffic colors. After this operation, we can apply any of the strategies described in section 3.4 to û instead of q.

These two methods regularize the traffic data, but they do not exploit the information from the pollution sensors, apart from the average value z. In order to assimilate data of both types, it is possible to use a combined least-squares fit of the form

A(x) = argmin v∈Vn z -v(v obs ) 2 2 + λ q -q 2 2 (V) ,
where λ > 0 balances the contributions of z and q, and q = P(v) = n k=1 ĉk q k for the coefficients ĉ ∈ R n such that v = n k=1 ĉk u k . However this approach did not perform well in practice, so we did not include it in the numerical experiments.

In the last method, if m ≥ n and λ tends to 0, the prediction û does a least squares fit of u at the available measurement points v obs i . In general, it is possible to enforce û(v obs ) = u(v obs ) by applying a correction to the prediction. This post-process, called Parameterized Background Data-Weak method, was originally introduced in [START_REF] Maday | A parameterized-background dataweak approach to variational data assimilation: formulation, analysis, and application to acoustics[END_REF] and has been analyzed and extended in a series of papers such as [START_REF] Binev | Data assimilation in reduced modeling[END_REF][START_REF] Cohen | Optimal reduced model algorithms for data-based state estimation[END_REF][START_REF] Cohen | Nonlinear reduced models for state and parameter estimation[END_REF][START_REF] Cohen | Nonlinear approximation spaces for inverse problems[END_REF]. The whole approach has found numerous applications, including pollution dispersion [START_REF] Hammond | PBDW: A non-intrusive reduced basis data assimilation method and its application to an urban dispersion modeling framework[END_REF].

It would of course be possible to gain in accuracy, by considering more refined equations for pollution dispersion, which capture additional physical properties, and thus by encoding these properties into the reduced space V n . One could for instance think of advection by wind, vertical fluxes or stratification of the atmosphere depending on the temperatures, changes in the chemical equilibrium between NO and NO 2 caused by cloud coverage and precipitations [START_REF] Li | Quantifying the nitrogen isotope effects during photochemical equilibrium between NO and NO 2 : implications for δ 15 N in tropospheric reactive nitrogen[END_REF], as well as local turbulent effects near the sensor stations. Note that, if nonlinear equations are involved, V n can be a nonlinear approximation space defined through a chart of n parameters, and approximation guarantees are more difficult to obtain [START_REF] Cohen | Nonlinear approximation spaces for inverse problems[END_REF].

Super-Learning as a collaborative approach

To gain in accuracy over each individual model, one can combine a set of p available mappings A 1 , . . . , A p coming from the previous methods, and build a super-learner

F(X , U ) p -→ F(X , U ) S : (A 1 , . . . , A p ) -→ S(A 1 , . . . , A p ),
where F(X , U ) denotes the set of functions from X to U . The most simple merger, usually called aggregator in statistics, amounts to taking a linear combination

S ω (A 1 , . . . , A p ) = p i=1 ω i A i ,
for some weights ω = (ω 1 , . . . , ω p ) expressing the confidence in each individual model. More sophisticated strategies involve nonlinear combinations and compositional structure. One could think of using a first model to obtain a rough estimation, and compose it with a second model performing refinements based on its output. This is already an underlying idea in our constructions, where we start with the spatial average, and add spatially-dependent corrections. The physical models involve one more compositional step, since they are of the form A src • û(q).

Neural networks constitute another prominent example of nonlinear super-learners: one could treat A 1 (x), . . . , A p (x) as inputs, and train the parameters ω by minimizing an empirical loss. We would like to emphasize here that properly training the super-learner requires to implement a nested leave-one-out strategy: one should first train each parametrized submodel by leave-oneout, and then optimize the neural network with another leave-one-out step, in order to avoid overfitting. As a consequence, at least two observation points are removed from the training set of the submodels, which may cause a loss of accuracy, especially when the number m of observations is small.

In our application, the neural network super-learner performed slightly worse than its linear counterpart, which should come as no surprise in view of the above observation.

Reconstruction benchmarks and Leave-One-Out

There are several ways to quantify the quality of a reconstruction map A : X → U . Ideally, given a state u ∈ U and the associated observations x ∈ X , one would like to find A such that the error u -A(x) U is as small as possible. Assuming that (u, x) is a random variable with distribution π ∈ Prob(U × X ), we define the performance of A as the L 2 norm of the error

e(A) 2 := U ×X u -A(x) 2 U dπ(u, x),
which acts as a good compromise between the worst-case error and the average error. In addition, although the state u has in principle H 2 regularity, we asses the spatial error also in L 2 , that is, we take U = L 2 (G). Unfortunately, finding A minimizing e(A) is out of reach for several reasons. First, we don't know the distribution π, nor even its support, which is the set of all possible states and observations. Second, given u ∈ U , we cannot evaluate u(r) at any point r ∈ G, making the computation of u -A(x u ) U intractable.

Concerning the first issue, as we have access to hourly data on a large period of time, we can replace the integral over π by an empirical average

e(A) 2 ≈ 1 T test t∈Ttest u(t) -A(x)(t) 2 U
over the set T test of 374 test times. Assuming that these states are independent, this approximation induces an error of order

E( u 2 U ) 1/2 |T test | ≈ 41 19.3 ≈ 2.1 µg/m 3 ,
which is totally acceptable in view of the noise level on the measurements. The second obstacle is more tricky, because we only know u at a very limited number m of fixed positions, and because these observations are also needed for constructing A. Ignoring the last issue leads to an systematic underestimation of e(A), as we detail below.

Assume that the observation points v obs for some linear space V n ⊂ U of dimension n. This setting is valid for most of our methods, with V n the set of constant functions in the case of A avg (of dimension n = 1), but also V n = span{r → c r i } 1≤i≤m in the case of A krig (of dimension n = m), and V n the reduced basis in the methods based on physical modeling.

Then, for A * the optimal map taking values in

V n A * = argmin A :X →Vn linear E( u -A (x) 2 U ) = E π (u|x),
we obtain, by applying Pythagoras theorem both for • U and • m ,

e(A) 2 = E( u -A(x) 2 U ) ≥ E( u -A * (x) 2 U ) = E( u -A * (x) 2 m ) ≥ E( u -A(x) 2 m ),
where the central equality comes from the assumption that the v obs i are random. This proves that u -A(x) 2 m is a biased estimator for e(A) 2 as soon as one of the inequalities is strict, that is, as soon as A(x) = A * (x). Note that separating the training and test data by splitting the set of time indices T is not sufficient, since the algorithm will then fit its prediction to the station locations, without generalization guarantees to the rest of the domain G.

As a consequence, we must separate the stations into a training set and test points. In order to compute an unbiased estimator of e(A) 2 with minimal variance, while keeping the maximal number of stations in the training set, we proceed to leave-one-out cross-validation. This procedure is very standard and has been used in other works on pollution reconstruction, such as [START_REF] Criado | Data fusion uncertainty-enabled methods to map street-scale hourly NO 2 in Barcelona: a case study with CALIOPE-Urban v1.0[END_REF]. For 1 ≤ i ≤ m, denote

u 2 m\i = 1 m -1 j =i |u(v obs j )| 2 and A i = argmin A:X →Vn linear 1 T train t∈T train u -A(x) 2 m\i .
Assuming that the station locations are independent random variables, the cross-validation estimator of the error

e CV (A) 2 := 1 T test t∈Ttest 1 m m i=1 |u(v obs i ) -A i (v obs i )| 2 is unbiased, since E(e CV (A) 2 ) = E(|u(v obs i ) -A i (v obs i )| 2 ) = E( u -A i (x) 2 U ) = e(A),
with the difference that x contains only m -1 direct evaluations of u this time.

Numerical results

We have implemented and tested numerous variants and combinations of the models from section 3. This was done thanks to a Python code we have developed, which can be found at https://github.com/agussomacal/CityPollutionModeling. The interested user could add its own models for further testing. In this section, we summarize the most important results that emerge from our tests. We report on the performance of the following reconstruction strategies:

• Spatial average: We take a simple spatial average, as in section 3.1. The resulting error serves as a baseline, which we expect to beat with the other more sophisticated models.

• BLUE: As explained in section 3.2, BLUE can be seen as an estimate of the optimal linear reconstruction method. It can be used as a benchmark of the best performance that we can expect of linear methods. Note that, in principle, nonlinear strategies could be more accurate than BLUE. However, we will see in our experiments that none of our methods achieves such accuracy.

• Kriging: We apply the kriging method depicted in section 3.3 with an exponential kernel.

The parameters σ and C are obtained by fitting an exponential to the correlation between training stations (that is, we do not into account correlations with the station that is set aside for testing) as a function of their distance (see fig. 4).

• Source: We apply a linear source model, as described in section 3.4, with temperature θ and wind w as extra regressor variables.

• Physical-PCA: We apply the second physical model from section 3.5, using a gaussian kernel to smooth the node traffic data. The four reduced spaces V n , associated to the four traffic colors, each consist of the first 10 principal components of the corresponding traffic data, as observed in the training set. After the smoothing and projection operations, we assemble the variables θ, w and the q v c on each node v into a vector s = (q v green , q v orange , q v red , q v dark-red , θ, w) ∈ R 6 and build a polynomial model of degree 2:

T α (s) := 6 j=1 α j s j + 6 j,k=1 α jk s j s k ,
where α ∈ R 42 is computed following the lines of section 3.4.

• Physical-Laplacian: We apply the first physical model from section 3.5, with the reduced space consisting of the first 5 eigenvectors of the graph laplacian. After projection into the subspace, we take the 4-colour traffic values on each node (q n c ) •,i and obtain their degree-3 polynomial combinations. Finally we apply a neural network consisting of two hidden layers of 20 neurons each and a ReLU activation function. The neural network is trained with ADAM optimizer with early stopping to prevent overfitting.

• Ensemble: We apply an ensemble model, as described in section 3.6, combining the Kriging method A krig , the Source model A src and the Physical-Laplacian model A lapl . We train each of them separately and compute the following linear combination:

A ens (x)(r) = ω(r)A krig (x)(r) + 1 -ω(r) 2 A src (x)(r) + 1 -ω(r) 2 A lapl (x)(r),
with a weight function ω(r) = exp(min 1≤i≤m |r -v obs i |/δ), where δ = 800 m. Essentially, we favour Kriging when r is close to one of the sensor stations, and average the predictions of models using local or global traffic information otherwise.

In fig. 5, we show the root mean square error (in µg/m 3 ) for each model's predictions on the test times T test and tested stations i:

e RMSE (A, i) := 1 T test t∈Ttest |u(v obs i ) -A i (v obs i )| 2 1/2 .
Note that the cross-validation error e CV (A) from section 4 is just the 2 -average of these errors over all stations. For the tests, we only keep the 10 stations located in the interior of Paris. The remaining 3 are set aside because they have a significant proportion of missing values (in average 10% of the data is lacking in each of these stations), and because they lay close to the border of the image, making the surrounding traffic information incomplete.

The shaded blue area is the region corresponding to errors smaller than the reference BLUE model. On the opposite side, the shaded red area marks situations in which the prediction is worse that the spatial average baseline. The white margin in between indicates the region where we expect feasible improvements.

We first notice that using a linear source model already yields reliable improvements with respect to the spatial average baseline. It fails, however, in HAUS and OPERA stations due to the absence of traffic information around the former, as the Google Maps symbol for the Paris Opera is drawn over the location of the sensor. This affects the predictions on both stations but most prominently on HAUS. This problem can be alleviated if we average traffic information over a bigger region, as done in Pysical-PCA thanks to the Gaussian smoothing, at the expense of losing precision on other stations like PA18.

The Kriging model manages to produce enhanced predictions in both HAUS and OPERA stations because of their proximity and correspondingly high correlation in pollution values. However, for other stations, especially those further from the center, the performance highly deteriorates.

With the Physical-Laplacian model, we get further improvements in 6 stations compared to the linear source model, while losing some advantage in the remaining 4. Finally, the ensemble method, by combining two traffic models and the Kriging method, manages to exploit the advantages of each in a pretty decent way. It yields predictions that beat or equal the spatial average baseline on all stations, and that reach the best average error among all our tested methods. In fig. 6, we show an example of pollution maps generated with this last model.

Conclusion and future works

In this work, we have showed that it is possible to leverage pollution sensor data, meteorological information and Google Traffic images to create pollution maps in real time. In particular, we explained how to build statistical, physics-based and ensemble reconstruction strategies by posing the problem of pollution state estimation on metric and quantum graphs. Furthermore, the right combination of these techniques produced systematic improvements over the proposed baselines, namely the Spatial average and Kriging.

Neither our linear reconstruction strategies nor our nonlinear ones could beat the BLUE benchmark that indicates the accuracy of the best linear estimator. We conjecture that this is due to the limited amount of stations giving us spatial information on the pollution field, and to the indirect nature of traffic data. Regarding the last point, even though the volume of traffic information is large, it still remains of reduced utility. This is due to the fact that we only measure the fluidity of traffic, instead of the actual amount of passing vehicles, which is the relevant variable directly impacting emissions. One can then hope to obtain further improvements by following a similar approach with better suited data.

Another limitation is the unavailability of local pollution averages for the city of Paris. Having access to this kind of data through intensive measurement campaigns lasting a few weeks, but employing hundreds to thousands of sensors, as done in [START_REF] Criado | Data fusion uncertainty-enabled methods to map street-scale hourly NO 2 in Barcelona: a case study with CALIOPE-Urban v1.0[END_REF] for the city of Barcelona, would give a much more precise baseline, and allow to learn corrections to the local average instead of the global one.

Finally, in setting the problem on the graph, we did not take into account the vicinity of open spaces like parks or rivers, nor the topology and variations in altitude. This is particularly visible in fig. 6, where the parks of Boulogne and Vincennes are colored in red because of the surrounding highways, and the absence of small internal streets. We leave the inclusion of such relevant features to future studies. continuous on [0, e ] for any edge e ∈ E, which implies in particular the continuity of u along any path in G. The space of square-integrable functions Finally, the Sobolev space

H 1 (G) = C(G) ∩ e∈E H 1 (r e ([0, e ]))
is also a Hilbert space, for the norm The restriction to C(G) in the definition of H 1 (G) stems from the fact that functions in H 1 (r e ([0, e ])) are continuous (because their domain is one dimensional), which automatically implies that functions in H 1 (G) must be continuous also at the vertices. In the same vein, one has to impose restrictions on the derivatives of u, such as Newmann-Kirchoff boundary conditions, for functions u ∈ H 2 (G).

Figure 1 :

 1 Figure 1: Cropped Google Map screenshot of Paris and the m = 13 available stations in the study: red dots represent the projection of the station locations to the nearest vertex in the graph of streets, while the blue crosses correspond to the exact position of the station.

Figure 3 :

 3 Figure 3: The metric graph downloaded from Open Street Maps, with the edges that never had Google Traffic activation in red, and the edges remaining after filtration in yellow

  j , with parameter values C = 0.968 and σ = 33.4km obtained by fitting the station data in our case.

Figure 4 :

 4 Figure 4: Correlation between stations as a function of the distance. The vertical slashed red line marks the maximal separation between vertex and station (165m) which still lays in the zone of high correlation.

i|u(v obs i )| 2 corresponding to an empirical version of u 2 UA:X →Vn linear 1 T

 221 are distributed uniformly at random on G, define the discrete semi-norm , and consider map A solution tomin train t∈T train u -A(x) 2 m

Figure 5 :

 5 Figure 5: Root mean square error on tested stations for the different proposed methods

Figure 6 :

 6 Figure 6: Pollution map for the ensemble model at 8am on March 1st, 2023. Based on the predictions on the node, one can linearly extrapolate pollution values even outside of the graph edges. Note that the fine variations in pollutant concentration (between 45 and 55 µg/m 3 ) seem to trace the main circulation axes.

L 2 (

 2 G) = e∈E L 2 (r e ([0, e ]))is a Hilbert space when endowed with the inner productu, v L 2 (G) := G u(r)v(r)dr =e∈E e 0 u(r e (x))v(r e (x))dx.

  e (x)) 2 dx + G d(u • r e ) dx 2 dx.

We extracted data from the Airparif database, which can be found at https://data-airparif-asso. opendata.arcgis.com

The permission to use Google Maps data for non-profit research is stated here: https://about.google/ brand-resource-center/products-and-services/geo-guidelines/#google-maps

See https://www.windguru.cz. For the wind, we combined the absolute wind speed with the wind direction to obtain a vector in R 2 .

The RGB value of each color is given by: green = (99, 214, 104), orange = (255, 151, 77), red = (242, 60,

50), dark-red = (129, 31, 31)

Acknowledgments and disclosure of funding: The authors would like to thank Prof. Albert Cohen, Prof. Joubine Aghili, Dr. Rachida Chakir, Dr. Vivien Mallet, and Dr. Fabien Brocheton for fruitful discussions, and for preliminary work on the extraction of applicable data.

This work was done in the framework of the research project "Models and Measures" funded by the Parisian City Council (Emergences grant program). In addition, Matthieu Dolbeault acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) -Project number 442047500 through the Collaborative Research Center "Sparsity and Singular Structures" (SFB 1481).

A Metric graphs

Here we recall several notions about graphs that are necessary in our developments. The presentation is based on the book [START_REF] Berkolaiko | Introduction to quantum graphs[END_REF], which provides a comprehensive introduction to quantum graphs, and on the paper [START_REF] Arioli | A finite element method for quantum graphs[END_REF], which develops finite element discretizations of elliptic operators in quantum graphs. We sometimes narrow down the generality of certain notions for the purposes of the present paper.

A combinatorial graph G = (V, E) is a collection of a finite number of vertices V and of edges E ⊂ V × V connecting pairs of vertices. We restrict our attention to undirected graphs where no orientation is assigned to the edges, and denote |V| and |E| the number of vertices and edges, respectively.

We will work with connected graphs, where any two vertices v, w ∈ V are connected by at least one path (v, v 1 ), (v 1 , v 2 ), . . . , (v k , w) made by consecutive adjacent edges in E. A connected graph becomes a metric graph if we assign a length e > 0 and a local coordinate r e (x), for x ∈ [0, e ], to each edge e = (v, w) ∈ E, in such a way that r e (0) = v and r e ( e ) = w.

In our case, the crossroads V are embedded in R 2 through their geographical coordinates, and the streets r e ([0, e ]) ⊂ R 2 are differentiable curves with no loops. However, as done very often, we redefine them as simple straight lines joining the two vertices. Regardless of the choice of the edge curves, the points r in a metric graph G are thus not only its vertices but also all intermediate points on the edges as well, parametrized by the local coordinates r e :

As the name suggests, any metric graph can be endowed with a natural metric as follows. The distance between two vertices v, w ∈ V is as usual defined as the length of the shortest path connecting them. This notion of distance between vertices is then extended in a natural way to any two points possibly lying on different edges, by further adding the local coordinates along these edges.

We may now introduce function spaces and linear differential operators on a metric graph G. The space of continuous functions C(G) contains the functions u : G → R such that u • r e is


Raphaël Achddou

Yann Gousseau

Saïd Ladjal

Hybrid training of denoising networks to improve the texture acutance of digital cameras

Keywords: Image denoising, Deep learning, Image quality assessment

In order to evaluate the capacity of a camera to render textures properly, the standard practice, used by classical scoring protocols, is to compute the frequential response to a dead leaves image target, from which is built a texture acutance metric. In this work, we propose a mixed training procedure for image restoration neural networks, relying on both natural and synthetic images, that yields a strong improvement of this acutance metric without impairing fidelity terms. The feasibility of the approach is demonstrated both on the denoising of RGB images and the full development of RAW images, opening the path to a systematic improvement of the texture acutance of real imaging devices.

Introduction

In order to correctly visualize a photograph, its corresponding RAW image undergoes a complex sequence of development operations including white balancing, demosaicking, tone mapping, and image restoration operations such as deblurring and denoising. Camera manufacturers implement proprietary algorithms fine-tuned for each setting of each camera. As a result, the overall image quality is a combination of hardware characteristics (quality of the lens, size of the sensor) and software performances. In order to fairly assess the quality of an imaging device, independent agencies have defined standard tests and ISO protocols. Each of these tests focus on a specific characteristic such as chromatic aberrations, noise reduction, or texture rendering.

Recently, with the increase in computational power and the advent of deep learning for image processing, more and more digital image processing stages can be replaced by learned neural networks [START_REF] Gharbi | Deep joint demosaicking and denoising[END_REF]. Recent works already aim at completely replacing the full image development pipeline with a single neural network, producing impressive results in standard conditions [START_REF] Ignatov | Replacing mobile camera isp with a single deep learning model[END_REF] or extremely low-light conditions [START_REF] Chen | Learning to see in the dark[END_REF]. Moreover, light neural network architectures can now be integrated in embedded systems, e.g. on smartphone devices. Neural methods present another key advantage: one can easily optimize their response to specific test images, by including them in training databases.

For the specific task of texture rendering evaluation, Cao et al. [START_REF] Cao | Measuring texture sharpness of a digital camera[END_REF] first presented a protocol quantifying the ability of an imaging pipeline to preserve texture information. This is obtained through the frequential response of the system to dead leaves images with a specific perceptual metric called texture acutance. These images are known for their invariance properties, as well as statistical properties making them close to natural images (non Gaussianity, scaling property, distribution of the spectrum and gradient), as studied in [START_REF] Alvarez | The size of objects in natural and artificial images[END_REF][START_REF] Gousseau | Modeling occlusion and scaling in natural images[END_REF][START_REF] Lee | Occlusion models for natural images: A statistical study of a scale-invariant dead leaves model[END_REF]. This quality evaluation protocol later became an ISO standard to measure the preservation of textures [START_REF]Photography -Digital cameras -Part 2: Texture analysis using stochastic pattern[END_REF] and is now used by classical camera scoring protocols. In a different direction, Achddou et al. [START_REF] Achddou | Synthetic images as a regularity prior for image restoration neural networks[END_REF] showed that image restoration networks could be trained from synthetic images only, using databases of dead leaves images.

Inspired by these results, we propose, in this paper, to train a denoising neural network on natural and dead leaves images, to jointly optimize a new metric derived from the texture acutance and the classic data fidelity metrics on natural images. After presenting some related works on image restoration in Section 2, we first introduce in Section 3 the texture acutance metric and the corresponding perceptual loss for image restoration networks. We then show in Section 4 that we can strongly improve the texture acutance metric without impairing performances on natural images, first for the task of Additive White Gaussian Noise removal (AWGN) and then for the development of RAW images. These results open the path to an automatic improvement of standard quality evaluation tests.

Related Works

The goal of image restoration is to retrieve a clean image from distorted observations. In many cases, the distortion process can be modeled as follows: y = Ax + n, where x is the theoretically perfect image, y the distorted observation, A is a linear operator and n is some noise.

In order to solve this problem, a first class of methods are based on prior hypotheses on the distribution of natural images. These methods try to impose regularity properties on the restored solutions. For instance, wavelet shrinkage methods [START_REF] Donoho | Ideal Spatial Adaptation by Wavelet Shrinkage[END_REF][START_REF] Donoho | Minimax estimation via wavelet shrinkage[END_REF] or DCT-filtering methods [START_REF] Yu | DCT Image Denoising: a Simple and Effective Image Denoising Algorithm[END_REF] reconstruct an image assuming that the targeted images can be well approximated by a sparse decomposition. In turn, variational methods based on the total variation [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF][START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF] assume that the image gradient follows a Laplacian distribution. Based on the assumption of self-similarity, non-local methods leverage the redundancy in the image content. This is either done by weighted averaging (Non Local Means [START_REF] Buades | A non-local algorithm for image denoising[END_REF] Non Local Bayes [START_REF] Lebrun | A nonlocal bayesian image denoising algorithm[END_REF]) or by collaborative filtering (BM3D [START_REF] Dabov | Video denoising by sparse 3D transformdomain collaborative filtering[END_REF]).

Over the past decade, learning-based approaches for image restoration have developed drastically. After the success of neural networks for high-level computer vision tasks [START_REF] He | Deep residual learning for image recognition[END_REF], these methods have been adapted to image restoration through the use of generative models [START_REF] Zhang | Beyond a gaussian denoiser: Residual learning of deep cnn for image denoising[END_REF][START_REF] Zhang | Ffdnet: Toward a fast and flexible solution for cnnbased image denoising[END_REF]. Rather than using prior hypotheses, the parameters of the neural networks are tuned in a long optimization process to directly minimize the reconstruction error in a black box manner. For the training to succeed, these methods require large databases of pairs of distorted and clean images. Even though they are hard to interpret, they surpassed priorbased methods on most image restoration benchmarks by a large margin for a wide variety of tasks such as image denoising [START_REF] Zhang | Ffdnet: Toward a fast and flexible solution for cnnbased image denoising[END_REF], demosaicking [START_REF] Gharbi | Deep joint demosaicking and denoising[END_REF] etc.

Following these initial works, recent papers extended the use of deep learning methods to real-world problems of image restoration such as RAW image denoising [START_REF] Liu | Invertible denoising network: A light solution for real noise removal[END_REF][START_REF] Anwar | Real image denoising with feature attention[END_REF]. Ignatov et al. [START_REF] Ignatov | Replacing mobile camera isp with a single deep learning model[END_REF] and Chen et al. [START_REF] Chen | Learning to see in the dark[END_REF] also propose to fully replace the image development pipeline by a learned neural network, producing surprisingly good results. However, acquiring datasets of real-world pairs of distorted and clean RAW images is a cumbersome task [START_REF] Abdelhamed | A high-quality denoising dataset for smartphone cameras[END_REF][START_REF] Chen | Learning to see in the dark[END_REF], which often requires complex post-processing algorithms. In order to ease the training process, a more restrained approach consists in modeling the distortion process accurately, and to synthesize them accordingly [START_REF] Wei | Physics-based noise modeling for extreme low-light photography[END_REF].

Going further, Achddou et al. [START_REF] Achddou | Synthetic images as a regularity prior for image restoration neural networks[END_REF][START_REF] Achddou | Fully synthetic training for image restoration tasks[END_REF] proposed to train image restoration neural networks on generated dead leaves images in order to completely circumvent the data acquisition process, reaching performances close to the networks trained on real images, for various image restoration tasks. These images indeed exhibit statistical properties close to those of natural images [START_REF] Alvarez | The size of objects in natural and artificial images[END_REF][START_REF] Gousseau | Modeling occlusion and scaling in natural images[END_REF][START_REF] Lee | Occlusion models for natural images: A statistical study of a scale-invariant dead leaves model[END_REF] even though they depend from few parameters. Following [START_REF] Achddou | Synthetic images as a regularity prior for image restoration neural networks[END_REF], similar synthetic databases were also used to pre-train image classification networks [START_REF] Baradad | Learning to see by looking at noise[END_REF] and disparity map estimators [START_REF] Madhusudana | Revisiting dead leaves model: Training with synthetic data[END_REF]. Prior to these works, dead leaves images were used to assess the capacity of cameras to render textures properly. This idea was first presented in 2009 by Cao et al. [START_REF] Cao | Measuring texture sharpness of a digital camera[END_REF], which was later improved in the following references [START_REF] Cao | Dead leaves model for measuring texture quality on a digital camera[END_REF][START_REF] Artmann | Image quality assessment using the dead leaves target: experience with the latest approach and further investigations[END_REF]. We will present in detail these works in the following section.

3 Texture acutance : a frequential loss assessing texture preservation Dead leaves images were first introduced by Matheron in 1975 [START_REF] Matheron | Random sets and integral geometry[END_REF], with the aim of modeling porous media. It was later shown that if object sizes fulfill some scaling property, this model accounts for many statistics of natural images [START_REF] Alvarez | The size of objects in natural and artificial images[END_REF][START_REF] Lee | Occlusion models for natural images: A statistical study of a scale-invariant dead leaves model[END_REF]. To generate such images, shapes of random size, color and position are superimposed on top of each other until the whole image plan is covered. In the simplest set-up, these shapes are disks of random radiuses. An example of a dead leaves image is given in Fig. 1, along with its spectrum. A precise mathematical formulation of dead leaves images is given in [START_REF] Bordenave | The dead leaves model: a general tessellation modeling occlusion[END_REF].

Dead leaves images were first used for camera evaluation in 2009 by Cao et al [START_REF] Cao | Measuring texture sharpness of a digital camera[END_REF]. The proposed idea is to measure the response of a camera to a specific image target. Because of their invariances and statistical properties, the dead leaves model was chosen by the authors as the generation algorithm for the target. Among the desired properties, scale invariance is achieved when the disks radii follows a power law with α = 3. The dead leaves target is therefore generated with this parameter. Note that to ensure the convergence of the algorithm, bounding parameters r min , r max are required [START_REF] Gousseau | Modeling occlusion and scaling in natural images[END_REF]. In [START_REF] Cao | Measuring texture sharpness of a digital camera[END_REF], the authors evaluate the response of a camera to the dead leaves target by computing the ratio of the power spectra, resulting in a Modulation Transfer Function (MTF). At each position (m, n) for an (N, N ) image:

Texture acutance

MTF 2D (m, n) = | Ŷ (m, n)| | X(m, n)| ,
where Ŷ is the spectrum of the obtained image and X is the ground truth spectrum. In all that follows, we compute the image spectra on a greyscale version of the color image, obtained by the standard linear combination Grey = 0.2126R + 0.7152G + 0.0722B. The classical idea behind the MTF is that the ratio of the power spectra corresponds to the Fourier coefficients of the blur kernel induced by the camera, excluding non linear transforms often involved in the image development, as well as the impact of noise.

In order to account for the impact of noise in the estimation of the MTF, Artmann first proposed a corrected version of the MTF, by subtracting an estimate of the noise spectrum. However, the latter was computed by taking a photograph of a uniform grey surface, assuming an additive and signal independent noise model, which is far from reality. Moreover, some image development pipelines include a nonlinear noise reduction operation, which affect the estimation of the real MTF.

In the same paper [START_REF] Artmann | Image quality assessment using the dead leaves target: experience with the latest approach and further investigations[END_REF], Artmann proposes a new computation trying to correct these issues. Here, we consider the complex spectrum of a reference digital dead leaves target X, rather than the estimate of the power spectrum | X| in the spatial domain. In the previous version, the phase information was lost. This is necessary in the context of camera calibration since phase information is reliable only if a registration algorithm is applied. In the context of training a denoising algorithm registration between noisy and restored image is supposed perfect and dealt with by the MSE-loss. Relying only on the amplitude of the spectrum meant that we could not differentiate frequencies which were already in the target and information that was added by the imaging device. Therefore noise and non linear functions had an impact in previous computations.

The proposed method, which we call MTF cross uses the cross power density between the target and the obtained image ϕ XY (m, n), and the auto power density ϕ XX (m, n). More precisely,

ϕ XY (m, n) = Ŷ (m, n) X * (m, n) and ϕ XX (m, n) = X(m, n) X * (m, n).
Given these quantities, the MTF becomes :

MTF cross (m, n) = ϕ XY (m, n) ϕ XX (m, n) . (1)
Since the dead leaves target is rotationally invariant, so is its spectrum. We therefore express the MTF as a 1D function by averaging it on concentric rings of width 1. The MTF becomes :

MTF 1D (k) = 1 C k #C k × MTF cross , where C k = (i, j) ∈ [-N/2, N/2] 2 |(k -1) 2 ≤ |i 2 + j 2 | < k 2 corresponds to a ring of radius k and #C k is its cardinal.
Though the full MTF 1D is a good indicator of the camera's capacity to render textures, it is more helpful to compute a single score. To that end, the texture acutance [START_REF] Cao | Measuring texture sharpness of a digital camera[END_REF] is defined as a weighted sum of the MTF 1D , with weights defined by a contrast sensitivity function (CSF), inspired by the slanted edge Spatial Frequency Response (SFR), used to evaluate the sharpness of a camera. Our visual system is indeed more sensitive to some frequencies than others. In that regard, the CSF models the sensitivity of the visual system to spatial frequencies expressed in cycle/degree.

Based on the physiological analysis of the contrast sensitivity of infants and monkeys led by Movshon and Kiorpes [START_REF] Movshon | Analysis of the development of spatial contrast sensitivity in monkey and human infants[END_REF], the chosen formula to model the CSF is :

CSF(ν) = a.ν c .e -bν , where ν is a spatial frequency expressed in cylces/degree, parameters are fixed as b = 0.2, c = 0.8, and a is a normalizing parameter so that Nyquist 0 CSF(ν)dν = 1, where Nyquist corresponds to the theoretical maximum spatial frequency of the device. Given this formula, the texture acutance score can be written as :

A = Nyquist 0 CSF(ν).MTF 1D (ν)dν.
Note that we need to convert spatial frequencies in cycles/degree to a digital frequency in cycles/pixel for homogeneity. To do so we use the following formula:

f spatial = 1 α f digital ,
where α is the viewing angle. The latter depends on viewing conditions with the equality α = 180 π arctan( P D ), where P is the pixel size and D is the viewing distance, assumed to be equal to 0.2 mm and 1m respectively. This corresponds to a maximal spatial frequency of 40 cycles/degree which is approximately the limit of the human visual system.

The perfect MTF corresponds to a constant function equal to 1, meaning that the frequential content has been perfectly restored by the camera for every frequency. This leads to an acutance A = 1. An acutance greater than 1 indicates that some frequential content was added to the image, probably because of noise or sharpening. An acutance lower than 1 indicates that some frequencies have been lost.

Acutance loss for image restoration CNNs

In [START_REF] Achddou | Synthetic images as a regularity prior for image restoration neural networks[END_REF], the authors showed that models trained on mixed databases (natural and synthetic images) perform on par with models trained on natural images only, while improving results on dead leaves image targets. We believe we can improve the frequential response of models trained on mixed sets, by using the acutance score in a loss function.

In the context of AWGN removal for color RGB images, the noisy image corresponds to Y = X + n where X is a ground truth dead leaves image of size (N, N, 3). The denoising network f θ produces an estimate of the clean image Z = f θ (Y ). For our restoration problem, we can consider that the denoising network is analogous to the camera which acquires the dead leaves target. We can compute MTF cross for the denoising network using Formula (1), based on the computation of the digital spectrum of both X and Z.

The obtained MTF cross is turned into a 1D signal as described above. For faster computation, concentric ring masks are stored in GPU so that the computation of MTF 1D can be accelerated with parallel computing. Since the best possible acutance is 1, we define the acutance loss function as :

L acutance (Y, X) = |1 -A(f θ (Y ), X)|,
which penalizes both adding or removing frequential information. In order to get a complete loss function, we add to it the L 2 loss (i.e. the Euclidian distance), the initial fidelity term of the network. Indeed, the acutance loss L acut is computed solely on an aggregation of the Fourier spectrum and is therefore blind to the spatial organisation of the image and can not replace an MSE-loss. When training on dead leaves images the loss is therefore

L = L 2 + λ.L acut ,
where λ is a weighting parameter.

Since we train the image denoiser on both natural images and dead leaves images, we compute the acutance loss only on the dead leaves images in a minibatch D of size K and the L 2 loss for all images. The formation of minibatches during training indeed randomly samples images from the mixed set. Thus, the loss in a batch becomes :

L batch = 1 K K i=0 ||x i -f θ (x i + n i )|| 2 2 + λ m T 1 K i=0 m i .L acut (f θ (x i + n i ), x i ), (2)
where m is a masking vector of size K such that m i = 1 if x i is a dead leaves image, or m i = 0 otherwise. In order to count the number of dead leaves images we sum this masking vector which is given by m T 1.

Image denoising results with FFDNet

We choose to train the FFDNet network [START_REF] Zhang | Ffdnet: Toward a fast and flexible solution for cnnbased image denoising[END_REF] to illustrate the impact of the perceptual loss we presented in the previous section. We adapt the training scheme of the network to the present problem as follows. First, we increase the size of the training patches from (50, 50, 3) to (100, 100, 3). The reason for this is that the estimation of the 1D-MTF on a small patch is not sufficiently accurate. Keeping the same rings' width would result in fewer estimates for the 1D-MTF.

On the other hand, decreasing the rings' width would lead to noisier estimates. Therefore, we perform the training with larger patches. Second, we reduce the batch size from 64 to 32 during training to decrease the memory footprint. We use 150000 samples, made of 100000 natural image patches and 50000 synthesized dead leaves patches. The other training hyper-parameters remain unchanged, such as the number of epochs or the learning rate decaying schedule.

Quantitative evaluation

In order to show that the proposed scheme indeed has the potential to improve the texture acutance without impairing the usual PSNR evaluation of the performances on natural images, we compute both these metrics for various values of λ, the weighting parameter in Equation [START_REF] Achddou | Fully synthetic training for image restoration tasks[END_REF]. We consider values of λ ∈ [0, [START_REF] Achddou | Fully synthetic training for image restoration tasks[END_REF][START_REF] Anwar | Real image denoising with feature attention[END_REF][START_REF] Cao | Measuring texture sharpness of a digital camera[END_REF][START_REF] Ignatov | Replacing mobile camera isp with a single deep learning model[END_REF]50,100,200,500]. Moreover, we also compute the classical SSIM metric and the perceptual metric PieAPP recently introduced in [START_REF] Prashnani | Pieapp: Perceptual image-error assessment through pairwise preference[END_REF].

The models are evaluated numerically on two datasets. First, we evaluate the data fidelity by computing the PSNR, SSIM and PieAPP metrics on the Kodak24 dataset, a benchmark test set of 24 natural images. Second, we evaluated the acutance metric on a test set of synthesized dead leaves images.

We report, in Table 1, the numerical evaluation of the trained models. We observe a similar behaviour for the tested noise levels σ = 25 and σ = 50. In both cases, we notice that the standard evaluation metrics, i.e., the PSNR and SSIM, are not affected by the increase of the weighting parameter λ until λ = 20. For values greater than λ = 100 these metrics decrease rapidly. On the other hand, the acutance metric keeps improving until λ = 100 and then reaches a plateau. This table shows that we can optimize the texture acutance without impairing classic denoising evaluation. The perceptual evaluation with the PieAPP metric suggests that, for high noise values, the perceptual image quality is slightly enhanced by the addition of the acutance loss. Some results can be visualized in Fig. 3 (please zoom in the electronic version of this document). The result with and without using the acutance loss appear quite close, despite the strong improvement of the texture acutance measurement. Nonetheless, we can notice some improvements in the preservation of low-contrast details in the first row. Moreover, the contrast is also rendered better when training with the acutance loss. Finally, on the third row, details on the dead leaves images are better preserved using the acutance loss. On the second and third row, we see that the network trained with natural images sometimes hallucinates details, which are removed when training with dead leaves images.

Spectral preservation

For mixed trainings of FFDNet, the texture acutance score is greatly improved when using the corresponding loss, which is expected. However, the acutance score only gives a partial information about the MTF of the trained network.

In order to further understand the impact of the acutance loss on the spectral preservation ability of the network, we compute its MTF as described next. We compute the 1D-MTF from the denoised image and the original image for each dead leaves image of the synthetic test set. Since the 1D-MTF depends on the image's content, which differs from image to image, we average the obtained MTF over the whole dataset. In Fig. 4, we report the MTF of FFD-Net trained with and without the acutance loss (with λ = 50) for a noise level σ = 25. Recall that a perfect MTF should be equal to one. We can observe that for low to medium frequency, the MTF of the model trained with the accutance loss is much closer to one. Actually, the values for low frequency exceed one which is one way the system can improve the acutance and which indeed is a limitation of the approach. For high frequency, the gap between the two MTF is smaller, probably as a result of the profile of the CSF function, which quickly decreases for high frequency, see Fig. 1. This behavior, as well as the addition of low frequency, could be modified by considering alternative CSF functions and can be easily integrated into our framework. This could also further improve the preservation of details on examples such as those of Fig. 3. In this paper, we have decided to keep the original definition of the acutance, since our main goal is to show that this standard measure of the ability to preserve texture can be greatly improved without impairing the other aspects of image quality.

RAW image denoising

As a proof of concept, we extended our experiments to real-world image denoising on the SIDD benchmark [START_REF] Abdelhamed | A high-quality denoising dataset for smartphone cameras[END_REF] for cameraphones denoising. To that end, we trained the same denoising network with a U-Net architecture to denoise real RAW images and synthesized RAW dead leaves images. This network produces a RAW denoised image from a RAW input noisy image. To simulate RAW noise for dead leaves images, we used a Poisson-Gaussian model with realistic noise parameters. Unlike Gaussian noise removal, the loss is here a combination of the L 1 loss and the acutance loss : L = L 1 +λL acut . For RAW images, the acutance computation differs slightly. In order to convert a RAW image to a grey-scale image, we first pack the (H, W ) image in a (H/2, W/2, 4) RGGB tensor, then we average them in a single (H/2, W/2) grey array, by weighting each channel with the white balance parameters. We ran the training for λ ∈ [0, 10, 100]. We report the numerical results obtained in Table 2. In comparison with λ = 0, the PSNR is still good for λ = 10, while the RAW acutance is largely improved. This improvement also translates in a better acutance in the RGB domain, which was not seen during training. This metric is computed on the denoised images developed with a standard ISP. This experiment shows that we can improve camera evaluation without impairing the image quality in the case of a full camera development pipeline. For λ = 100, the PSNR noticeably decreases while the RAW acutance reaches a plateau.

Conclusion

In this work, we have shown that a specific training of image restoration neural networks can greatly improve a standard evaluation metric quantifying the preservation of textures, without impairing classical performance evaluation criteria. As a proof of concept, we extended the use of the acutance loss for realworld image denoising networks, showing that the proposed framework can improve a complete RAW images development pipeline. Considering that the texture acutance metric is routinely used to evaluate digital camera, this founding has potential important practical applications.

Aknowledgements This work is supported by the project MISTIC (ANR-19-CE40-005).

3. 1

 1 Dead leaves images

Fig. 1 :

 1 Fig. 1: Grey level dead leaves image and its associated digital spectrum in logarthmic scales. The theoretical value is a straight line.

Fig. 2 :

 2 Fig. 2: Diagram explaining the computation of the acutance metric

Fig. 3 :

 3 Fig. 3: Comparison of FFDNet results on two natural images and on a dead leaves image. From left to right: original image, noisy image, image denoised with standard FFDNet, image denoised with FFDNet trained on a mixed database without the acutance loss, and finally with the acutance loss.

Fig. 4 :

 4 Fig. 4: Comparison of the MTF evaluated with FFDNet trained on a mixed database with or without the acutance loss, on the whole dead leaves image test set.

Table 1 :

 1 Denoising results of FFDNet trained with different weighting coefficients of the acutance loss for two noise levels. Each cell contains the PSNR, SSIM, PieAPP evaluated on Kodak 24 , and the Acutance metric evaluated on a test set of dead leaves images. Best results in blue, second results in red.

		λ	0	2	5	10	20	50	100	200	500
		PSNR ↑	31.88 31.87 31.88 31.87 31.88 31.85 31.77 31.65 31.56
	σ = 25	Acutance ↓ SSIM ↑	0.034 0.029 0.023 0.020 0.015 0.013 0.012 0.012 0.012 0.877 0.876 0.875 0.875 0.877 0.876 0.873 0.872 0.869
		PieAPP ↓	0.568 0.596 0.598 0.602 0.586 0.587 0.591	0.62	0.612
		PSNR ↑	28.81 28.79 28.80 28.80 28.80 28.76 28.66 28.58 28.42
	σ = 50	Acutance ↓ SSIM ↑	0.084 0.078 0.073 0.053 0.035 0.026 0.022 0.022 0.022 0.791 0.788 0.789 0.789 0.790 0.789 0.786 0.783 0.779
		PieAPP ↓	0.932 0.956 0.952 0.948 0.912 0.925 0.934 0.940 0.953

Table 2 :

 2 Denoising results obtained by training a denoising Unet for real RAW images evaluated on the SIDD test set of cameraphone images[START_REF] Abdelhamed | A high-quality denoising dataset for smartphone cameras[END_REF]. We report the PSNR, Acutance RAW and acutance RGB metrics. Best results in bold.

	λ	0	10	100
	PSNR RAW	51.31	51.24	50.61
	Acutance RAW	0.018	0.011	0.012
	Acutance RGB	0.086	0.061	0.049