Madura Pathirage

Danyang Tong

Flavien Thierry

Gianluca Cusatis

David Grégoire

Gilles Pijaudier-Cabot
email: gilles.pijaudier-cabot@univ-pau.fr

Discrete modeling of concrete failure and size-effect

Keywords: Concrete failure, Size-effect, Effect of geometry, Lattice Discrete Particle Model, Numerical modeling, Energy dissipation

  L'archive ouverte pluridisciplinaire

Introduction

Size-effect in quasi-brittle materials has been extensively investigated in the literature experimentally, theoretically, and numerically. In these materials, a reduction in strength is observed when the structural size increases in geometrically similar structures. This phenomenon has been confirmed for concrete (see for instance the work of [START_REF] Bažant | Determination of fracture energy from size effect and brittleness number[END_REF]; [START_REF] Hoover | Comprehensive concrete fracture tests: description and results[END_REF]; [START_REF] Grégoire | Failure and size effect for notched and unnotched concrete beams[END_REF]; C ¸aglar and S ¸ener ( 2016)) and its importance was acknowledged in the civil engineering community, as it was for the first time incorporated in the most recent ACI standard in 2019 (ACI-Committee, 2019). Two main reasons explain this size-effect [START_REF] Bažant | Fracture and Size Effect in Concrete and Other Quasibrittle Materials[END_REF][START_REF] Bažant | Scaling of Structural Strength[END_REF][START_REF] Bažant | Probabilistic Mechanics of Quasibrittle Structures: Strength, Lifetime, and Size Effect[END_REF]. The first one is the redistribution of stress due to stable crack propagation and release of stored energy into the fracture front. The second reason is the spatial randomness in material properties. In this study, only size-effect due to stress redistribution is studied, i.e. the energetic size-effect which is, in passing, purely deterministic. In concrete, the size of the fracture process zone (FPZ) is non-negligible as compared to the structural size and a complete non-linear theoretical fracture mechanics framework does not exist. Simplified analytical formulations such as Bažant's size-effect law [START_REF] Bažant | Scaling of Structural Strength[END_REF] only provide an approximate description of concrete scaling law. In order to account for the release of the stored energy in the FPZ, and for the development and propagation of the fracture front for different specimen sizes and shapes, one must carry out accurate numerical simulations.

Several types of models have been proposed over the years to describe concrete fracture and size-effect. One can mention for instance the cohesive model [START_REF] Elices | The cohesive zone model: advantages, limitations and challenges[END_REF], the crack-band model [START_REF] Bažant | Crack band theory for fracture of concrete[END_REF], non-local continuum damage models [START_REF] Bažant | Nonlocal integral formulations of plasticity and damage: survey of progress[END_REF], and discrete models [START_REF] Bolander | Discrete mechanical models of concrete fracture[END_REF]. In all the aforementioned formulations, two features are essential in capturing size-effect in strain softening materials such as concrete: (i) crack localization and (ii) existence of an internal characteristic length related to the size of the heterogeneity. In continuum models, these two ingredients are phenomenologically defined through constitutive laws. In this respect, random lattice or particle models such as the ones described by Cusatis et al. (2011a) or [START_REF] Eliáš | Stochastic discrete meso-scale simulations of concrete fracture: Comparison to experimental data[END_REF] are considered superior. The actual mix design and particle size distribution can be simulated to produce a realistic heterogeneous internal structure made of interacting aggregates.

In addition, the randomness in spatial distribution of particles reproduces the statistically isotropic nature of concrete and eliminates directional mesh bias during the fracturing process. Finally, a recent argument in favor of lattice particle models is their ability to capture the effect of stress parallel to cracks on the FPZ size and the induced size-effect, which was demonstrated to be significant in concrete (Nguyen et al., 2020a,b). In this study, one of such models, namely the Lattice Discrete Particle Model (LDPM) (Cusatis et al., 2011a,b) which simulates concrete at the coarse aggregate level is adopted.

In order to assess the capabilities of LDPM in simulating fracture and predicting size-effect, a comprehensive experimental data set on concrete fracturing [START_REF] Grégoire | Failure and size effect for notched and unnotched concrete beams[END_REF] was considered in this study. This data set is among one of the very few available in the literature that includes three-point bending tests on notched and unnotched beams and encompasses a large range of beam depths. It should be mentioned that this set of data was previously used in several studies involving other numerical models. More specifically, an integral-type non-local model was used but was found ineffective in capturing correct size and geometry effects [START_REF] Grégoire | Failure and size effect for notched and unnotched concrete beams[END_REF]. In another study, a lattice model [START_REF] Grassl | Meso-scale modelling of the size effect on the fracture process zone of concrete[END_REF] was able to simulate the data with a good accuracy. It must however be noted that the latter model falls into the miniscale category which implies computationally prohibitive simulations for more complex geometries that require an increase in structural size or account for the third dimension in the out-of-plane direction. It should also be mentioned that the model used in that study was bidimensional and that to the authors knowledge, there has been no attempt to extend the simulations to the 3D case. More recently, a study proposed to simulate concrete fracturing using a local isotropic damage constitutive model of the Rankine type through the crack band model [START_REF] Barbat | Structural size effect: Experimental, theoretical and accurate computational assessment[END_REF].

Although a good agreement with experimental data was found, this type of smeared crack model is limited in capturing complex cracking mechanisms and realistic crack tortuosity observed in concrete, which might have a non-negligible effect on energy dissipation during fracture. As a matter of fact, this effect is especially important in the case of unnotched beams where a large damage zone is generated before its collapse to a single propagating crack that is ultimately tortuous and not straight. Last but not least, a re-implementation of LDPM which includes stochasticity in material parameters was recently used to simulate the data used in this paper [START_REF] Eliáš | Fracture in random quasibrittle media: I. discrete mesoscale simulations of load capacity and fracture process zone[END_REF] with a very good accuracy. It should be however emphasized that the latter study did not assess the capability of the model to predict size-effect. Indeed, the study used all the load-crack mouth opening displacement (CMOD) curves for all the beam sizes and notch lengths in the calibration process, with the exception of the unnotched beams and the smallest sizes that were discarded.

In this work, the relevant LDPM parameters were calibrated on a single beam configuration, leaving the remaining fracture test results for blind predictions, including splitting tests. In addition and in contrast with most available studies, a rigorous quantification of the prediction accuracy was performed. It is the wish of the authors that this work would, to some extent, pave the way for more detailed and quantitative analyses of goodness of fits and prediction quality in future numerical studies within the concrete community, using the same set of experimental data but possibly other data. This would ultimately allow one to compare different model capabilities, develop the most effective models and abandon the less accurate ones. The proposed approach and preliminary results were presented in a recent conference (Pathirage et al., 2022b,a;[START_REF] Pijaudier-Cabot | Surface and size effects on elasticity and fracture[END_REF]. They are here largely extended and commented.

The Lattice Discrete Particle Model

In order to simulate the FPZ formation and propagation in concrete necessary to capture size-effect, the Lattice Discrete Particle Model is here adopted. This model originally proposed by Cusatis and coworkers (Cusatis et al., 2011a,b) simulates the mechanical interactions among major material heterogeneities, i.e. coarse aggregates in concrete. Over the years, this model has been used to simulate other granular quasi-brittle materials such as mortar (Pathirage et al., 2019b;[START_REF] Han | Lattice discrete particle modeling of size effect in slab scratch tests[END_REF], fiber reinforced concrete and engineered cementitious composites [START_REF] Schauffert | Lattice discrete particle model for fiber-reinforced concrete. i: Theory[END_REF][START_REF] Rezakhani | Influence of steel fiber size, shape, and strength on the quasi-static properties of ultra-high performance concrete: Experimental investigation and numerical modeling[END_REF][START_REF] Feng | Engineered cementitious composites using chinese local ingredients: Material preparation and numerical investigation[END_REF], unreinforced and reinforced stone masonry [START_REF] Mercuri | Computational modeling of the out-of-plane behavior of unreinforced irregular masonry[END_REF][START_REF] Angiolilli | Fiber reinforced cementitious matrix (frcm) for strengthening historical stone masonry structures: Experiments and computations[END_REF][START_REF] Angiolilli | Lattice discrete particle model for the simulation of irregular stone masonry[END_REF][START_REF] Mercuri | Masonry vaulted structures under spreading supports: Analyses of fracturing behavior and size effect[END_REF][START_REF] Mercuri | On the collapse of the masonry medici tower: An integrated discrete-analytical approach[END_REF], shale [START_REF] Li | A multiscale framework for the simulation of the anisotropic mechanical behavior of shale[END_REF], or cycling in concrete [START_REF] Zhu | Lattice discrete particle modeling of concrete under cyclic tension-compression with multi-axial confinement[END_REF]. LDPM was also coupled to multi-physics models describing cement hydration from microscale simulations, heat transfer and mositure diffusion, alkali-silica reaction, creep, aging [START_REF] Alnaggar | Lattice discrete particle modeling (ldpm) of alkali silica reaction (asr) deterioration of concrete structures[END_REF]Pathirage et al., 2019a;[START_REF] Yang | Computational modeling of temperature and relative humidity effects on concrete expansion due to alkali-silica reaction[END_REF][START_REF] Yang | Computational modeling of expansion and deterioration due to alkali-silica reaction: Effects of size range, size distribution, and content of reactive aggregate[END_REF], or more recently self-healing in concrete [START_REF] Cibelli | A discrete numerical model for the effects of crack healing on the behaviour of ordinary plain concrete: Implementation, calibration, and validation[END_REF].

Internal geometry

In order to generate the LDPM skeleton, spherical particles are placed in a volume of material from the largest to smallest size. This placement follows a prescribed particle size distribution that is based on the actual concrete mix design with the maximum and minimum in a prismatic sample. In order to simulate the interaction between particles, a lattice system 96 is generated by means of a Delaunay tetrahedralization with the centers of particles. A dual 97 tessallation is then performed which finally produces a system of polyhedral cells enclosing 98 the spherical particles. Figure 1

(a) (b) (c) (d) (e) (f) (g) (h) (i) (a) (b) (c) (d) (e) (f) (g) (h) (i) (a) (b) (c) (d) (e) (f) (g) (h) (i)

Elastic, tension, and tension-shear constitutive behaviors

If x i and x j denote the positions of nodes i and j, adjacent to the facet k, the facet strains are defined as:

e k = [e N k e M k e L k ] t = n t k [[u k ]] l k m t k [[u k ]] l k l t k [[u k ]] l k t (1)
where e N k is the normal strain component, and e M k and e L k are the tangential strain com-

ponents, [[u k ]] = u j -u i is the displacement jump corresponding to facet k, l k = ∥x j -x i ∥ is
the distance between the two nodes, n k = (x j -x i )/l k and m k and l k are two unit vectors mutually orthogonal in the facet plane projected orthogonally to the line connecting the adjacent nodes. The traction vector is defined as For the sake of readability, the subscript k that designates the facet is dropped. The elastic behavior is formulated through linear relations between the normal and shear stresses, and the corresponding strains as follows:

t k = [t N k t M k t L k ] t ,
t N = E N e N , t M = E T e M , t L = E T e L (2)
where E N = E 0 and E T = α 0 E 0 . E 0 ≈ E/(1 -2ν) and α 0 ≈ (1 -4ν)/(1 + ν) are the effective normal modulus and the shear-normal coupling parameter, respectively, and E is the macroscopic Young's modulus and ν is the macroscopic Poisson's ratio.

Because of the mesoscale nature of the model, concrete fracturing in mode I opening is always accompanied by shear at facets. This is a realistic feature since it is experimentally observed that most fracture paths are located at the interface between aggregates and cement paste. Therefore, the cohesive fracture behaviors in tension but also in tension-shear are important. This cohesive fracture occurs for e N > 0. One can define an effective strain as e = (e 2 N + α 0 (e 2 M + e 2 L )) where e N,max (τ ) = max

τ ′ <τ [e N (τ ′ )
] and e T,max (τ ) = max

τ ′ <τ
[e T (τ ′ )]. The strength limit of the effective stress that defines the transition between pure tension and pure shear is written as

σ 0 (ω sn ) = σ t -sin(ω sn ) + ( sin 2 (ω sn ) + 4α 0 cos 2 (ω sn )/r 2 st ) 1 2 2α 0 cos 2 (ω sn )/r 2 st ( 4
)
where r st = σ s /σ t is the ratio of the shear strength to the tensile strength, σ s is the shear strength and σ t is the tensile strength. The post-peak softening modulus is controlled by the effective softening modulus

H 0 (ω sn ) = H s /α 0 + (H t -H s /α 0 ) (2ω sn /π) nt , in which H t = 2E 0 /(l t /l -1), H s = r s E 0 and n t is the softening exponent.
Typically, the values of n t = 0.2 and r s = 0 are assumed and are fixed. l t is the tensile characteristic length defined as l t = 2E 0 G t /σ 2 t and G t is the mesoscale fracture energy.

Static equilibrium equations and numerical implementation

Finally, one can write the static linear and angular momentum equilibrium equations of each LDPM cell as follows:

k∈F I A p k t k = 0 , k∈F I A p k c k × t k = 0 (5)
where F I is the set containing all the facets of a generic polyhedral cell I, A p k = A k n t n k is the area of the projected facet k, n is the orientation of the tetrahedron edge associated to facet k and n k is the unit vector orthogonal to facet k of area A k (Cusatis et al., 2011a). c k is the vector that represents the distance between the center of facet k and the center of the cell. The model was implemented within a dynamic explicit scheme, with a central difference algorithm for time integration. Although the actual equations that are solved numerically are dynamic (see Cusatis et al. (2011a) for more details), the inertia terms are absent in the two expressions in Equation 5 because all the simulations presented next were performed under quasi-static conditions. In other words, loading rates were small enough to ensure the kinetic energy in the system would not exceed 5% of the internal energy throughout the analysis.

Regression statistics

Before describing the modeling and calibration process, it is important to introduce the regression statistics used in this work to perform the model parameter identification and list the different statistical indicators used to assess the quality of fit and the quality of predictions. In this study, a frequentist approach is proposed.

Experimental data and model response

Let y j k be the experimental results for replicate k for a given test j. For instance, y j 1 , y j 2 , and y j 3 can be the measured forces for three replicates of a three-point bending test denoted j. In order to compare the experimental data, one needs at some point to compute the arithmetic difference between experimental data and model response Y j k as

ε j k = y j k -Y j k .
This expression implies that a model response, which depends on the given spatial distribution of particles, corresponds to a specific replicate of the experimental test. This is not entirely correct since the internal structure of the experimentally tested material is only statistically replicated by the model.

To overcome this problem, a solution consists in: (i) averaging the experimental results of all replicates to obtain the mean data y j , (ii) running multiple simulations with different spatial distributions of particles and averaging the model responses to get Y j , and (iii) comparing the experimental and numerical results in terms of mean responses. In this process, there is a loss of information since the variance of the data is not known. Such a task would require the use a true stochastic model which appears unnecessary with respect to the very good results presented in Section 5. In addition, this loss remains minimal as the interest in this work is on central range statistics, viz. mean values for normally distributed regression errors, and is reduced as the number of replicates increases.

Parameter estimation

Let x j be the predictor variable, for instance the displacements corresponding to the forces y j or Y j . One can first discretize in equally spaced intervals the predictor variable over its range of values and define x j i , y j i , and Y j i at index i = 1, ..., n j where n j is the number of discretized points. Next, if one denotes f j the model response corresponding to the simulation of test j and θ = [θ 1 ... θ p ] t the parameter vector containing the p unknown model parameters to be estimated, then, Y j i = f j (x j i , θ). Next, one can formulate the general minimization problem as follows:

Find the least square estimate θ which minimizes

S j = n j i=1 y j i -f j (x j i , θ) 2 (6)
where S j is the residual sum of squares for curve j. Using the least-square method is here well justified as the conditional variance of the data used later on can be considered uniform and because the residuals of the regression ε j i = y j i -f j (x j i , θ) can be generally assumed independent or non-correlated.

In the case of a single test, i.e. j = 1, the formulas stated earlier can be directly applied and S 1 can be minimized to obtain θ. When multiple tests are considered simultaneously in the parameter identification process, one needs to normalize the residual sum of squares for each test and generate a global objective function to be minimized. This can be done by defining a new residual ϕ j for each test j as follows:

ϕ j = S j n j i=1 y j i 2 (7)
It is interesting to note that the normalization in Equation 7becomes the Mean Absolute Percentage Error (MAPE) in the case of a single point estimate (i = n j = 1).

The global residual is defined as a linear combination of ϕ j for different tests j, which is consistent with the additive property of the square of the coefficient of variation of errors defined in Section 3.3. The general minimization problem is written as follows:

Find the least square estimate θ which minimizes ψ = N j=1

w j ϕ j (8)
where w j can be seen as the weights for each test j and take the values between 0 and 1 such that N j=1 w j = 1. N is the total number of tests considered. In absence of specific information about the tests, one can assume equal weights, i.e. ∀j, w j = 1/N .

Fit and prediction quality

Once the estimate θ is obtained, one can appreciate the fit and prediction quality for each individual test. For this purpose, two statistical indicators are introduced and will be used throughout the paper.

The first indicator is the coefficient of variation of the regression errors ω j which characterizes the ratio of the scatter band width to the data mean

ω j = s j ȳj where s j = S j n j -p (9)
s j is the standard error of the regression, n jp is called degrees of freedom, and ȳj = n j i=1 y j i /n j is the data mean.

The second one is the adjusted (r 2 ) j , also called the coefficient of determination which characterizes the ratio of the scatter band width to the overall spread of data and indicates what percentage of data variation is accounted for by the model response. By defining the corrected total sum of squares S j c = n j i=1 y j i -ȳj 2 and the standard deviation of all data sj = S j c /(n j -1), this indicator is written as:

r 2 j = 1 - (s j ) 2 (s j ) 2 (10)
In the case of a single point estimate, i.e. i = n j = 1, the aforementioned indicators are not valid. One can instead use the mean absolute percentage error which is written as MAPE j = 100|ε j /y j | where ε j = y jf j (x j , θ), y j = y j 1 and x j = x j 1 .

In order to evaluate and quantify an overall quality of fit and predictions, one simply needs to regroup all data points for each test with their corresponding model responses into one global set and compute the overall coefficient of variation of errors, ω, and the overall coefficient of determination, r 2 .

Let us underline that the coefficient of variation and the coefficient of determination are mathematically derived for linear model regression only and are, in theory, not valid when one performs regression or assesses quality of fit for non-linear (parameter-wise) models.

Nevertheless, they are often used in practice because of the absence of clear indicator for quality of fit, aside from the standard error of regression which however is unit-dependent, is not expressed in percentage, and is not bounded. In addition, one might argue that when the model is already fitted and the problem is well-conditioned, the residual sum of squares which is originally non-linear can be linearized around the least square estimate. Therefore, close enough to the estimate and when the original and linearized residual sums of squares are close enough, the use of these statistical indicators can be justified.

LDPM internal geometry

The parameters required to construct the LDPM geometry were first identified based on the actual mix design used in the experiments. The particle size distribution was numerically

(c) (d) (e) (f) (g) (h) (i) (a) (b) (c) (d) (e) (f) (g) (h) (i) (a) (b) (c) (d) (e) (f) (g) (h) (i)
(c) (d) (e) (f) (g) (h) (i)

Modeling and calibration process

The identification of the parameters in the constitutive laws describing elastic, tension, and tension-shear behaviors followed a two-step procedure. First, the normal modulus E 0 and α 0 related to the elastic behavior were calculated using the approximated formulas listed in Section 2.2, based on the mean values of the macroscopic elastic modulus and Poisson's ratio reported by [START_REF] Grégoire | Failure and size effect for notched and unnotched concrete beams[END_REF]. The values of E 0 = 57180 MPa and α 0 = 0.25 were obtained. strength σ t , the mesoscale fracture energy G t , and the shear-to-tensile strength ratio r st were (e), (g) and (h). One can observe that the numerical results fit well the experimental data.

(c) (d) (e) (f) (g) (h) (i) (a) (b) (c) (d) (e) (f) (g) (h) (i) D =
In terms of quality of fit, a MAPE of 3.4% was reached for the compression test, whereas the three-point bending test fit was characterized by a coefficient of variation ω = 6.8% and a coefficient of determination r 2 = 0.995. It is interesting to note that the ratio between the macroscopic splitting tensile strength f st reported in Section 5.1 and the compressive strength f ′ c shown in Figure 3(a) is approximately 9%, which is consistent with the range of values reported in the literature. Note also that the mesoscale tensile strength σ t and the compressive yielding strength σ c0 are not macroscopic properties but only model parameters.

As explained in Section 2, concrete failure is characterized by multiple mechanisms that are different and LDPM is able to simulate all these mechanisms. For each mechanism, there is a set of relevant model parameters (less than 4) which makes the total number of parameters to be 16. The elastic behavior modeled by two parameters (E 0 , α 0 ) and the mesoscale mixed mode fracture governed by three parameters (σ t , G t , r st ) have been identified. The remaining parameters were assumed based on the actual mix design and Section 5.3 by Cusatis et al. (2011b), namely the softening exponent n t governing the interaction between shear and tensile behavior during softening, (σ c0 , H c0 , κ c0 , E d ) defining the behavior of the facet normal component under compression and affecting the macroscopic behavior in compression, (µ 0 , µ ∞ , σ N 0 ) contributing to the LDPM response in compression, mainly the triaxial compressive behavior at high-confinement, (κ c1 , κ c2 ) governing the nonlinear evolution of the normal facet stress in compression, and finally β controlling the coupling between the mesoscale compressive behavior and the macroscopic triaxial compressive behavior. All parameters are listed in Table 1 for the sake of clarity.

Prediction results and discussion

The capability of the model to simulate fracture and predict size-effect was assessed by carrying out blind simulations, i.e. without adjusting model parameters, on splitting and size-effect tests for different beam sizes and notch lengths. [START_REF] Rocco | Influencia del tamaño y mecanismos de rotura del ensayo de compresión diametral[END_REF][START_REF] Bažant | Fracture and Size Effect in Concrete and Other Quasibrittle Materials[END_REF].

D s = 113 mm L s = 215 mm (a) (b) (c) (d) (e) (f) (g) (h) (i) (a) (b) (c) (d) (e) (f) (g) (h) (i)

Size-effect tests

Load-CMOD curves

The model was next used to predict the load-CMOD curves for all sizes and notch lengths. The simulations were performed under displacement control and quasi-staticity was ensured by applying loads at a constant loading rate of 1 mm s -1 . Figures 5(a)-(c) show the predicted curves for the notch-to-depth ratios of 0.5, 0.2, and 0, respectively. The solid line is the mean response of the three individual simulations in dashed lines. The experimental scatter is represented with a gray area. One can observe that the numerical simulations predict well the mechanical behavior in the elastic, near-peak and post-peak regimes for the different geometries and sizes. The prediction in the post-peak of the smallest size beam with a notch-to-depth ratio of 0.5 deviates from the experimental results. The reason might be that the number of aggregates in the ligament is not enough and makes the model too coarse. In addition, boundary effects play a significant role for such small specimens. In the case of the two unnotched beams with largest sizes, the response stops at the peak due to snapback, similarly to what was observed in the experiments. For the notched specimens, the FPZ is localized and emanates at the crack tip. It develops for the increasing load and finally reaches an ultimate size at the peak load (Figures 5(d) and (e)). The FPZ of the constant size then propagates through the ligament (Figures 5(g) and (h)), which explains the strain softening behavior observed at the macroscale. (ii) For the unnotched specimens, the FPZ initiates at the bottom surface of the sample, where the stresses approach the material tensile strength, and is diffused on a zone much larger in size as compared to the notched-beam case. As the load increases up to the peak, the damaged zone Mesoscale crack opening (µm) at peak becomes larger (Figure 5(f)). At the peak load, the FPZ eventually collapses to a single crack that propagates, whereas the surrounding strained material unloads. One can also observe that the final crack does not necessarily originate at mid-span. This phenomena shows the direct effect of material heterogeneity, realistically captured by the mesoscale model.

Fracture process zone and dissipated energy

(a) (b) (c) (d) (e) (f) (g) (h) (i) (a) (b) (c) (d) (e) (f) (g) (h) (i) (a) (b) (c) (d) (e) (f) (g) (h) (i) D =
α = 0.5 α = 0.0 (a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l) (a) (b) (c) (d) (e) (f) (g) (h) (i) (a) (b) (c) (d) (e) (f) (g) (h) (i) D =
= 0.2 α = 0 (a) (b) (c) (d) (e) (f) (g) (h) (i) (a) (b) (c) (d) (e) (f) (g) (h) (i) (a) (b) (c) (d) (e) (f) (g) (h) (i) Mesoscale crack opening (µm) at displacement 0.2 mm (j) (k) (l) (j) (k) (l)
Although not quantified here, the evolution of the FPZ is consistent with the one described in [START_REF] Lefort | Analysis by ripley's function of the correlations involved during failure in quasi-brittle materials: Experimental and numerical investigations at the mesoscale[END_REF] based on Ripley's function analysis on the same experimental data.

As a matter of fact, the very nature of the model allows one to assess shear and tensile behaviors during fracture propagation in contrast with continuum models where smear cracks over a region and are unable to capture local oriented events. In general, the LDPM facets are subject to both tensile and shear strains even though the test configuration is designed for mode I opening. One can look at the evolution of the dissipated energy computed from the increment of the dissipated energy density ẇd = 3(t with the dissipated energies. As expected, the total energy dissipated at the end of the test/simulation is smaller for the larger notch length, i.e. for a smaller ligament length. This is consistent with the acoustic energy obtained from acoustic emission on the same type of experiments reported in the work of [START_REF] Grégoire | Mesoscale analysis of failure in quasi-brittle materials: comparison between lattice model and acoustic emission data[END_REF].

N ėine N + t M ėine M + t L ėine L ) for e N >
Moreover, one can observe a sharper increase in the dissipated energy as the notch length decreases, which is consistent with the increase in brittleness as α tends to zero. Up to the peak load, the energy dissipated in shear is negligible as compared to the one in tension, for all cases including the unnotched beam. However in the post-peak regime, the energy dissipated in shear becomes more than half of the energy dissipated in tension, which proves that the post-peak behavior involves both shear and tensile forces in the meso-structure. The quality of predictions was also assessed. Figure 6 (a) shows the scatter between experimental and numerical loads for all geometries and sizes, on the entire load-displacement curves. Figure 6(b) shows the same scatter but only for the peak values. In both figures, the cloud of points are close to the 1:1 lines. The set of load points that deviate the most to the 1:1 lines corresponds to the largest unnotched beam (Figure 6(a)). The prediction overestimates the peak load by about 15% (Figure 6(b)). One might be tempted to invoke statistical size-effect related to the randomness in material properties, which is not captured by the deterministic version of LDPM used in this study. Such an effect can be effectively simulated by introducing random fields on the mesoscale model parameters as performed in the work of [START_REF] Eliáš | Fracture in random quasibrittle media: I. discrete mesoscale simulations of load capacity and fracture process zone[END_REF]. Nevertheless, the deterministic predictions are inaccurate only for the unnotched beam of size D = 400 mm. Further testing on larger unnotched beam sizes would help validate or invalidate a possible statistical effect for beam with the size D = 400 mm. In addition, the use of random fields would add a new level of complexity by making the inverse parameter identification problem ill-conditioned with respect to the limited reported experimental data.

(c) (d) (e) (f) (g) (h) (i) (a) (b) (c) (d) (e) (f) (g) (h) (i)
Overall, a coefficient of determination of r 2 = 0.94 and a coefficient of variation of ω = 36.8% were found when all load points are considered. For the peak loads only, the values of r 2 = 0.95 and ω = 23.5% were obtained. The values of the coefficients of variation show that the numerical predictions are very reasonable and are within the typical scatter observed in concrete fracture testing. These results demonstrate quantitatively the capability of the model to predict the effects of size and geometry during concrete failure.

Universal size-effect law and fracture parameters

Macroscopic properties of concrete such as compressive or tensile strength are usually determined in laboratory using standardized sample dimensions. On the other hand, sizeeffect tests applied in a certain range of sizes were shown to provide an accurate estimation of fracture properties [START_REF] Rilem | Size-effect method for determining fracture energy and process zone size of concrete[END_REF]. This method is preferred to the work of fracture method that provides apparent fracture properties which are geometry and size dependent. The sizeeffect method only requires the knowledge of peak loads and sample geometry [START_REF] Bažant | Fracture and Size Effect in Concrete and Other Quasibrittle Materials[END_REF]). It appears therefore interesting to compare the fracture parameters obtained using experimental results as reported by [START_REF] Grégoire | Failure and size effect for notched and unnotched concrete beams[END_REF] and parameters one could identify using the numerical predictions, keeping in mind that only one load-CMOD curve on a single size notched beam was sufficient for model calibration.

For this purpose, many approximated formulae exist and can be used [START_REF] Morel | Size effect in quasibrittle fracture: derivation of the energetic size effect law from equivalent lefm and asymptotic analysis[END_REF][START_REF] Cusatis | Cohesive crack analysis of size effect[END_REF][START_REF] Bažant | Universal size effect law and effect of crack depth on quasi-brittle structure strength[END_REF][START_REF] Di Luzio | Cohesive crack analysis of size effect for samples with blunt notches and generalized size effect curve for quasi-brittle materials[END_REF][START_REF] Chen | On interchangeability and selection of size effect and boundary effect experiments for characterization and prediction of quasi-brittle fracture of concrete[END_REF]. Typically, two parameters are enough to capture size-effect in notched threepoint bending tests [START_REF] Rilem | Size-effect method for determining fracture energy and process zone size of concrete[END_REF][START_REF] Planas | Generalized size effect equation for quasibrittle materials[END_REF][START_REF] Cusatis | Cohesive crack analysis of size effect[END_REF][START_REF] Chen | On interchangeability and selection of size effect and boundary effect experiments for characterization and prediction of quasi-brittle fracture of concrete[END_REF]. Unnotched specimens can also be simulated but the size cannot be too large as statistical size-effect becomes non-negligible [START_REF] Eliáš | Stochastic discrete meso-scale simulations of concrete fracture: Comparison to experimental data[END_REF]. In this study, the so-called Universal Size Effect Law (USEL) in its deterministic version [START_REF] Bažant | Universal size effect law and effect of crack depth on quasi-brittle structure strength[END_REF] was considered. This fitting formula bridges Type I size-effect which occurs in structures that fail at crack initiation from a smooth surface and Type II size-effect occurring in notched structures. It also covers the two distinct asymptotic behaviors at large size in the typical double-logarithm nominal strength versus structural size representation: (i) 1/2 slope corresponding to linear elastic fracture mechanics for Type II size-effect and (ii) a straight horizontal line corresponding to the elastic limit for Type I size-effect. The formula is written as:

σ N u = EG f g ′ 0 c f + g 0 D 1 - rc 2 f g ′′ 0 e -kα 2 4(l p + D)(g 0 D + g ′ 0 c f ) 1/r (11)
where σ N u = (3P u S)/(2W D 2 ) is the nominal strength corresponding to the peak load P u , G f is the fracture energy , and c f is the effective length of the FPZ. g 0 = g(α 0 ) is the dimensionless energy release rate, and g ′ 0 = g ′ (α 0 ) and g ′′ 0 = g ′′ (α 0 ) are its first and second derivatives, respectively, evaluated at the initial notch-to-depth ratios α 0 = 0.5, 0.2, 0.

Finally, r, k, and l p are empirical constants. Based on the expression of g(α) reported for example by [START_REF] Bažant | Fracture and Size Effect in Concrete and Other Quasibrittle Materials[END_REF] or [START_REF] Grégoire | Failure and size effect for notched and unnotched concrete beams[END_REF], the following values were computed for the geometry studied in this paper: g 0 = 2.96, g ′ 0 = 18.95, and g ′ 0 = 153.88 for α = 0.5, g 0 = 0.57, g ′ 0 = 3.17, and g ′′ 0 = 10.77 for α = 0.2, and g 0 = 0, g ′ 0 = 3.41, and These predicted fracture parameters are to be compared with the ones identified using the experimental results, i.e. G f = 42.6 N m -1 and c f = 25.7 mm [START_REF] Grégoire | Failure and size effect for notched and unnotched concrete beams[END_REF]: both fracture energy and effective length can be qualitatively considered close enough with respect to the scatter of typical experiments on concrete. This result suggests that one could perform only one fracture test on a specific geometry in the laboratory, calibrate and use LDPM to simulate different sizes (and possibly geometries) to finally identify the fracture energy and effective length of the FPZ through a size-effect law.

g ′′ 0 =-
(c) (d) (e) (f) (g) (h) (i) 10 2 10 3 10 0 10 1 Experimental Numerical (a) (b) (c) (d) (e) (f) (g) (h) (i) Experimental Numerical

Conclusions

In this study, a large set of experimental results on fracture and size-effect including the effect of geometry was simulated using the Lattice Discrete Particle Model (LDPM).

The load-CMOD curve of a single size notched beam under three-point bending, and the compressive strength were used for model calibration. The remaining experimental results, namely one splitting test and eleven three-point bending tests of different beam sizes and notch lengths, were used for model validation without parameter adjustment. The fracturing process was discussed along with a rigorous quantification of quality of fit and quality of predictions. Fracture parameters using the numerical results were identified through the Universal Size Effect Law (USEL). Based on the obtained results the following conclusions can be drawn:

• The compressive strength and the load-CMOD curve of a single notched beam size are sufficient for a complete model calibration in elastic, tensile and tensile-shear behaviors.

• The predictions on splitting tests and on the different beam geometries and sizes are overall in excellent agreement with the experimental data.

• The fracturing process is well captured by LDPM for both notched and unnotched beams.

• The dissipated energy in shear constitutes a large part of the total dissipated energy in the post-peak.

• The fracture energy can be estimated through a size-effect law by using the simulation results.

• The use of a stochastic model does not seem to be justified with respect to the limited range of beam sizes.

Declarations of interest

The Authors declare that there is no conflict of interest.

Acknowledgment

Partial financial support from the investissement d'avenir French programme (ANR-16-IDEX-0002) under the framework of the E2S UPPA hub Newpores and from the Ezponda project under the FEDER program is gratefully acknowledged. This research was supported in part through the computational resources provided for the Quest high performance computing facility at Northwestern University which is jointly supported by the Office of the Provost, the Office for Research, and Northwestern University Information Technology. 9. Appendix 9.1. LDPM constitutive equations for pore collapse, compaction, and frictional behavior

The pore collapse and material compaction occur for e N < 0. Strain hardening plasticity behavior due to high compressive hydrostatic deformation is computed through a strain-dependent boundary σ bc (e D , e V ) at each facet which limits the normal compressive stress component via the inequality -σ bc (e D , e V ) ⩽ t N ⩽ 0, in which the volumetric strain is defined as e V = ∆V /(3V 0 ), computed as the change between the current and the initial volume of each LDPM tetrahedron. While e V is the same for all the facets of a given tetrahedron, the deviatoric strain defined as e D = e N -e V changes. Finally, the strain-dependent boundary is calculated for three different cases σ bc (e D , e V ) = σ c0 for -e V ≤ 0, σ bc (e D , e V ) = σ c0 + ⟨-e V -ε c0 ⟩H c (r DV ) for 0 ≤ -e V ≤ e c1 and σ bc (e D , e V ) = σ c1 (r DV ) exp [(-e V -e c1 )H c (r DV )/σ c1 (r DV )] otherwise, where r DV = |e D |/e V for e V > 0 and r DV = -|e D |/(e V -e V 0 ) for e V ≤ 0 in which e V 0 = κ c3 e c0 . e c0 = σ c0 /E 0 defines σ c0 the mesoscale yielding compressive stress, e c1 = κ c0 e c0 is the strain at which the rehardening starts, κ c0 and κ c3 are material constants and σ c1 (r DV ) = σ c0 + (e c1 -e c0 )H c (r DV ). The function H c (r DV ) reads as H c (r DV ) = H c1 + (H c0 -H c1 )/(1 + κ c2 ⟨r DV -k c1 ⟩) where H c0 , H c1 , κ c1 and κ c2 are material constants.

The frictional behavior due to compression-shear occurs also for e N < 0. In the presence of compressive stresses, the shear strength increases due to frictional effects. The frictional behavior is simulated by means of a nonlinear Mohr-Coulomb model in which the internal friction coefficient varies from an initial value µ 0 to zero. The formulation can be written as σ bs = σ s + µ 0 σ N 0 -µ 0 σ N 0 exp(σ N /σ N 0 ) where σ s is the cohesion stress and σ N 0 is the transitional stress.

  (b) shows an example of two adjacent polyhedral cells 99 enclosing the spherical particles. One can also refer to Figures 22(b)-(d) where cells at the 100 surface of prismatic samples are depicted. The surface of each polyhedral cell is composed of 101 triangular facets where failure can potentially occur. On each facet, stresses and strains are 102 formulated in a vectorial form through constitutive equations. Figure 1(c) shows the three 103 unit vectors defined at a generic facet colored in red, in the normal direction and in the two 104 tangential directions. LDPM incorporates specific constitutive equations to describe tensile 105 fracturing with strain softening, cohesive and frictional shearing, and compressive response 106 with strain-hardening. Since this study focuses on concrete fracturing, the corresponding 107 constitutive laws are detailed in the following. The reader is referred to the Appendix for 108 the complete set of constitutive equations.

  109

Figure 1 :

 1 Figure 1: LDPM internal geometry: (a) particle placement in a prismatic volume, (b) two adjacent LDPM cells, and (c) triangular facet and vector orientations.

  4. Modeling and calibration process4.1. Summary of experimental campaignIn the experimental work by[START_REF] Grégoire | Failure and size effect for notched and unnotched concrete beams[END_REF], four sizes of geometrically similar prismatic specimens with four depths D = 50 mm, 100 mm, 200 mm, and 400 mm, the span-to-depth ratio S/D of 2.5, and the out-of-plane thickness of 50 mm were tested in three-point bending. Unnotched and notched samples with a notch length a and the notchto-depth ratios α = a/D = 0.5, 0.2, 0 were tested under CMOD control to obtain a stable post-peak response. In the case of unnotched beams, the legs of the extensometer were attached to the bottom surface of the beams at a distance from mid-span of half the beam depth to ensure crack initiation between the legs. Splitting tests on cylinders were also conducted to estimate tensile strength. In complement to fracture tests, cylindrical samples were tested under unconfined compression and measurements of elastic parameters and compressive strength were obtained. More details on the experimental program are given by[START_REF] Grégoire | Failure and size effect for notched and unnotched concrete beams[END_REF].

  reproduced following the procedure described in Yang et al. (2022) with a cut-off size d 0 = 4 mm and a maximum size d a = 10 mm. Figure 2(a) shows the experimental and numerical (a) (b)

Figure 2 :

 2 Figure 2: (a) Simulated particle size distribution, LDPM cells, geometries, and dimensions of the simulated beams with the notch-to-depth ratios of (b) 0.5, (c) 0.2, and (d) 0 viz. unotched beam, for D = 50 mm, 100 mm, 200 mm, and 400 mm.

269

  identified simultaneously based on: (i) the compressive strength obtained from cylinders 270 with the diameter D c = 74 mm and the height H c = 142 mm, and (ii) the entire load-271 CMOD curve corresponding to the medium size beam with the depth D = 200 mm and 272 α = 0.2. For the compression test, rigid plates were used on the top and bottom of the 273 specimens. Friction between the plates and the sample was simulated through a simple 274 Coulomb friction law with a friction coefficient µ = 0.13. Concerning the bending test, the 275 loads were applied directly on the surface nodes. Both compression and fracture tests were 276 simulated under displacement control with a constant velocity of 1 mm s -1 to ensure quasi-277 static conditions. For each test, three simulations were performed with different spatial 278 distributions of particles. The least square estimate of the parameter vector θ = [σ t G t r st ] t 279 was obtained by minimizing the overall residual ψ = ϕ C + ϕ 3PBT where ϕ C and ϕ 3PBT are the 280 residuals computed through Equation 7 for the compression and three-point bending tests, 281 respectively. The values of σ t = 2.9 MPa, G t = 45.5 N m -1 , and r st = 3.276 were obtained.

Figure 3 : 1 2

 31 Figure 3: Calibration results: (a) stress-strain curve of the unconfined compression test on the cylinder with the diameter D c = 74 mm and height H c = 142 mm; the empty circle designates the mean peak value, (b) failure mode at the peak load (c) load-CMOD curve of the three-point bending test on the notched beam with D = 200 mm and α = 0.2.

Figure 4

 4 Figure 4(a) shows the mean predicted nominal stress versus displacement curve repre-

  crack opening (µm) at peak Mesoscale crack opening (µm) at peak

Figure 4 :

 4 Figure 4: Prediction results: (a) nominal stress-displacement curve of the splitting test; the empty circle designates the mean peak value, and (b) failure mode at the peak load.

  328using the formula f st = (2P st,u )/(πD s L s ) where P st,u is the splitting peak force. One can 329 observe that the numerical mean splitting tensile strength is within the scatter of the exper-330 iments. A MAPE of 7.5% was found when comparing the experimental and simulation peak 331 values, which indicates a good accuracy in the prediction with respect to the scatter observed 332 in experiments. Figure 4(b) shows the mode of failure at the peak load. As expected and 333 in accordance with experimental observation, fracture initiates at the center where tensile stresses are the highest. The main crack is tortuous and has a slight eccentricity with respect 335 to the vertical line passing through the center. This is due to the inherent ability of the 336 mesoscale model to reproduce heterogeneity in the material, and is by the way routinely observed in experiments

Figures 5

 5 Figures 5(d)-(i) show the typical failure modes for the beams with the size D = 100 mm for two different displacement values. Two types of failure can be distinguished. (i)

Figure 5 :

 5 Figure 5: Prediction results: load-CMOD curves for (a) α = 0.5, (b) α = 0.2, (c) α = 0, the empty circles designate mean peak values; failure modes at the peak load for the samples with the depth D = 100 mm for (d) α = 0.5, (e) α = 0.2, (f) α = 0; failure modes at a displacement of 0.2 mm for the samples with the depth D = 100 mm for (g) α = 0.5, (h) α = 0.2, (i) α = 0; dissipated energies for (j) α = 0.5, (k) α = 0.2, (l) α = 0.

  0 and the volume of the cell I containing the facet V I = ( k A p k l k )/3. Summing over the entire volume of the sample, one can obtain the evolution in time (or displacement here) of the total dissipated energy W d and dissociate the individual contributions of the normal component W d,N and the tangential components W d,T . In the present simulations, W d,N almost coincides with the energy dissipated in tension. Figures 5(j)-(l) show the load-displacement curves for the beams with the size D = 100 mm for the three different notch lengths, together

Figure 6 :

 6 Figure 6: Statistics of the model predictions: (a) scatter of the measured versus predicted values of loads, and (b) scatter of the measured versus predicted values of the peak loads only. The dashed lines correspond to the 1:1 lines.

  16.75 for α = 0. By keeping the empirical constants r = 0.11, k = 113, and l p = 12.9 identical to the ones reported in the work of[START_REF] Grégoire | Failure and size effect for notched and unnotched concrete beams[END_REF], the model was fitted using the simulation data. The values of the fracture energy and effective length, G f = 39 N m -1 and c f = 16.8 mm were obtained with a coefficient of determination r 2 = 0.97 and a coefficient of variation of ω = 1.82%. Figures7(a)-(b) show the fitted model plotted together with the experimental and numerical data.

Figure 7 :

 7 Figure 7: Fitting with the Universal Size Effect Law: nominal strength versus size for (a) experimental data and (b) simulation results.

  where t N k is the normal component, and t M k and t L k are the shear components.

Table 1 :

 1 Values of the material model parameters used in the numerical simulations

	Designation	Symbol	Units	Value	Source
	Density	ρ	kg m -3	2338 Grégoire et al. (2013)
	Water-to-cement ratio	w/c	-	0.626 Grégoire et al. (2013)
	Maximum aggregate size	d a	mm	10	Grégoire et al. (2013)
	Minimum aggregate size	d 0	mm	4	Fixed
	Effective normal modulus	E 0	MPa	57,180	Identified
	Shear-normal coupling parameter	α 0	-	0.25	Identified
	Tensile strength	σ t	MPa	2.9	Identified
	Fracture energy	G t	N mm -1	45.5	Identified
	Shear strength ratio	r st	-	3.276	Identified
	Softening exponent	n t	-	0.2	Cusatis et al. (2011b)
	Compressive yielding strength	σ c0	MPa	120	Cusatis et al. (2011b)
	Initial hardening modulus ratio	H c0 /E 0	-	0.4	Cusatis et al. (2011b)
	Transitional strain ratio	κ c0	-	2	Cusatis et al. (2011b)
	Deviatoric strain threshold ratio	κ c1	-	1	Cusatis et al. (2011b)
	Deviatoric damage parameter	κ c2	-	5	Cusatis et al. (2011b)
	Initial friction	µ 0	-	0.2	Cusatis et al. (2011b)
	Asymtotic friction	µ ∞	-	0.0	Cusatis et al. (2011b)
	Transitional stress	σ N 0	MPa	600	Cusatis et al. (2011b)
	Densification ratio	E d /E 0	-	1.0	Cusatis et al. (2011b)
	Volumetric deviatoric coupling	β	-	0.0	Cusatis et al. (2011b)

LDPM was first validated on splitting tests for which peak loads for nine replicates are 320 reported by [START_REF] Grégoire | Failure and size effect for notched and unnotched concrete beams[END_REF]. For this purpose, three cylinders with the height L s = 215 321 mm and the diameter D s = 113 mm with different spatial distributions of particles were 322 simulated. The load was applied directly on particles at the surface of the cylinder. In 323 addition, the simulations were performed under displacement control with a constant loading 324 rate of 1 mm s -1 to ensure quasi-static conditions.


Anthony Ung

Thierry Ranchin

Lucien Wald

Christiane Weber

Jacky Hirsch

Gilles Perron

Joseph Kleinpeter

Cartographie de la pollution de l'air : une nouvelle approche basée sur la télédétection et les bases de données géographiques. Application à la ville de Strasbourg

come

Introduction

L'•tude de la qualit• de l'air s'est largement diversifi•e pour une meilleure connaissance du ph•nom…ne de pollution et de lutte contre la pollution atmosph•rique (voir par exemple le programme PRIMEQUAL-PREDIT). La qualit• de l'air est aujourd'hui un probl…me multidisciplinaire qui mobilise autant les sp•cialistes •pid•miologiques, les sp•cialistes en mod•lisation des transports, en •mission et transformation des polluants, les sp•cialistes en syst…mes g•ographiques et en t•l•d•tection que les autorit•s locales et les industriels. L'application de la loi sur l'air et l'utilisation rationnelle de l'•nergie (LAURE, d•cembre 1996) a conduit ' un d•veloppement de la surveillance de la qualit• de l'air [START_REF] Mate | Minist…re de l'Am•nagement du Territoire et de l'Environnement[END_REF]. La pollution atmosph•rique est aujourd'hui sensiblement r•duite par rapport aux situations qui pr•valaient il y a quelques d•cennies (HCSP r•sum•, 2000). Cependant, du fait de la concentration des populations et des activit•s en ville, et ' cause du flux croissant des •changes et de l'accroissement constant du parc et de la circulation automobile, elle constitue toujours un enjeu et une pr•occupation particuli…rement sensible (IPTS, 2000). En l'•tat actuel des connaissances, la question de sant• public reste ouverte et des recherches sont n•cessaires afin de renforcer les connaissances en mati…re d'exposition r•elle de la population ' la pollution atmosph•rique. L'•valuation de cette exposition consiste ' croiser les donn•es de r•partition g•ographique de la population et les donn•es de pollution (AIRPARIF). Il est ainsi n•cessaire d'am•liorer nos connaissances th•oriques des comportements spatialis•s de d•placements mais aussi d'acqu•rir une estimation pr•cise de la pollution atmosph•rique dans l'espace et dans le temps. Plusieurs outils existent aujourd'hui qui permettent de d•river des informations sur la qualit• de l'air autant dans l'espace que dans le temps. Ils rencontrent aujourd'hui des probl…mes de validation et de limitation. La principale difficult• est l'insuffisance des stations de mesures. Les villes ne peuvent se permettre d'•quiper tous les quartiers, les rues et les places d'un tel dispositif de surveillance de la qualit• de l'air. Leur co"t •lev• limite le nombre de stations ' une dizaine par ville. Or la variabilit• de la distribution spatio-temporelle des concentrations des polluants est •lev•e et un suivi optimal requerrait des observations spatialement et temporellement denses.

Pour 1 Les m€thodes actuelles de cartographie de la pollution de l'air en milieu urbain L'objectif de la cartographie de la qualit• de l'air est d'acc•der ' une estimation des concentrations des polluants atmosph•riques dans l'espace et dans le temps. Divers outils sont couramment utilis•s par les agences de surveillance de la qualit• de l'air [START_REF] Maneux | Quelles sont les techniques pour surveiller la qualit• de l'air ? Document de[END_REF]). Il est par exemple possible de construire une carte de distribution spatiale et temporelle des polluants par mod•lisation de leur transport, de leur d•p‰t et/ou encore de leur transformation chimique dans l'atmosph…re [START_REF] Giambanis | Numerical simulation of pollutant dispersion and photochemical kinetics over complex terrain[END_REF][START_REF] Leuzzi | Particle trajectory simulation of dispersion around a building[END_REF][START_REF] Sall…s | Mobile source emission inventory model: application to Paris area[END_REF]. Ces mod…les se distinguent suivant les •chelles spatiales et temporelles prises en compte mais aussi suivant les •quations de transport utilis•es. Une description plus d•taill•e de tous ces mod…les peut ƒtre trouv•e dans le rapport du European Topic Center on Air Quality, 'Ambient Air Quality' [START_REF] Moussiopoulos | Ambient air quality, pollutant dispersion and transport models[END_REF] [START_REF] Lohmeyer | A comparison of street canyon concentration predictions by different modellers: final results now available from the Podbi-exercise[END_REF]. Les autres m•thodes pour estimer la distribution spatiale des polluants font appel ' l'interpolation et/ou extrapolation des donn•es in situ de pollution. Ici aussi, les mod…les et m•thodes d'interpolation et d'extrapolation sont vari•s. Certains travaux pr•conisent le krigeage [START_REF] Frangi | Suivi spatio-temporel de la concentration en NO x en ‹le de France[END_REF][START_REF] Carletti | Kriging and bilinear methods for estimating spatial pattern of atmospheric pollutants[END_REF], d'autres recommandent la m•thode des plaques minces [START_REF] Ionescu | M•thodes math•matiques pour estimer le champ de concentration d'une polluant gazeux ' partir des valeurs mesur•es aux points dispers•s[END_REF] (Basly 2000[START_REF] Brivio | Atmospheric pollution and satellite remotely sensed surface temperature in metropolitan areas[END_REF][START_REF] Finzi | Landsat images of urban air pollution in stable meteorological conditions[END_REF][START_REF] Mclellan | Remote sensing of atmospheric turbidity variations by satellite[END_REF][START_REF] Poli | Study of the heat island over the city of Rome from Landsat-TM satellite in relation with urban air pollution[END_REF][START_REF] Sifakis | Potentialit•s de l'imagerie ' haute r•solution spatiale pour le suivi de la repartition de pollutions atmosph•riques dans la basse troposph…re. Etude du cas d'Ath…nes[END_REF], Sifakis et al. 1992[START_REF] Sifakis | La t•l•d•tection des voiles de pollution atmosph•rique et de la d•gradation de l'environnement dans la r•gion d'Ath…nes[END_REF][START_REF] Sifakis | Quantitative mapping of air pollution density using Earth Observations : a New Processing Method and Application to an Urban Area[END_REF][START_REF] Wald | Observing air quality over the city of Nantes by means of Landsat thermal infrared data[END_REF], sugg…rent l'existence d'une relation qui pourrait lier les mesures de qualit• de l'air fournies par les stations de mesure et les observations satellitaires. Dans le domaine de l'infrarouge thermique, des corr•lations fortes ont •t• trouv•es entre le canal thermique TM6 du satellite Landsat et les mesures de concentration des particules noires en suspension. Dans le domaine du spectre visible et proche infrarouge, [START_REF] Carnahan | Evaluation of atmospheric particulate concentrations derived from analysis of ratio thematic mapper data[END_REF] ont recherch• des rapports de canaux Landsat permettant une discrimination grossi…re des diff•rentes concentrations de particules en suspension. [START_REF] Richards | Remote sensing digital image analysis: an introduction[END_REF]

Les observations de la qualit€ de l'air sur la ville de Strasbourg

L'association agr••e de surveillance de la qualit• de l'air ' Strasbourg est l'ASPA (Association pour la Surveillance et l'•tude de la Pollution atmosph•rique en Alsace). Elle est en charge de trente-deux stations de mesures r•parties sur toute l'Alsace dont quatorze sur la zone d'•tude de la CUS. Une pr•sentation compl…te de l'association est pr•sente sur le site : http://www.atmo-alsace.net. Ces stations de mesures sont situ•es soit dans des cabines am•nag•es par l'ASPA (2 m*3 m*4 m), soit dans des locaux oŽ une surveillance du mat•riel est possible (•cole, h‰pital). Les polluants effectivement mesur•s sont le dioxyde de soufre (SO 2 ), le monoxyde d'azote (NO), le dioxyde d'azote (SO 2 ), le monoxyde de carbone (CO), le dioxyde de carbone (CO 2 ), l'ozone (O 3 ), et les particules en suspension (PM). Les proc•dures d'•valuation de la qualit• de l'air •voluent au cours du temps et font aujourd'hui l'objet d'harmonisation (ADEME, 2000). Les recommandations actuelles de mesure de polluants sugg…rent de mesurer des particules plus fines. Ainsi les stations de mesures de particules, qui mesuraient les concentrations en PM 13 (c'est ' dire la concentration en particules de taille sup•rieure ' 13 m), mesurent d•sormais les PM 10 , voire mƒme les PM 2.5 . Les donn•es disponibles sont des concentrations de ces diff•rents polluants par quart horaire. A partir de ces mesures, il est possible de suivre en continu l'•volution des concentrations de polluants dans les lieux oŽ sont localis•es les stations de mesures.

Les observations satellitaires

•om•trique ont •t• calcul•s ; les images n'ont pas •t• re-projet•es dans ce r•f•rentiel SPOT.

La base de donn€es g€ographiques BD TOPO de l'IGN

Une caract•ristique •vidente de la ville r•side dans l'architecture mƒme des bOEtiments : les bOEtiments en centre ville poss…dent une toiture, et les maisons pr•sentes sont essentiellement des maisons ' colombages et aux toits pentus. La ville de Strasbourg pr•sente une configuration de bOEti plus resserr•e dans le centre ville (la petite France). Les •tudes de simulations num•riques montrent que l'architecture des bOEtiments et leurs configurations influencent l'•coulement des vents et a fortiori la pollution locale. A titre d'exemple, la place de la Cath•drale est c•l…bre pour les courants d'air tourbillonnaires. Notre •tude int…gre une base de donn•es g•ographiques BD TOPO de l'IGN [START_REF] Michel | Le territoire vu du ciel[END_REF][START_REF] Veillet | Evaluer la qualit• de la BD Topo : l'approche de l'IGN[END_REF]. Cette base, g•or•f•renc•e en Lambert I, contient une description de l'occupation du sol et de sa morphologie (r•seaux routiers et ferr•s, bOEti, v•g•tation, hydrographie, relief, limites administratives, etc.). De plus amples informations peuvent ƒtre trouv•es sur le serveur •ducatif d•di• ' l'information g•ographique : http://pse.ensg.ign.fr/

Une nouvelle approche basée sur la télédétection et les bases de données géographiques

Une fois •mis par une source, l'•volution d'un polluant dans l'atmosph…re d•pend du vent, de la stabilit• de l'atmosph…re, du rayonnement solaire, et de la topographie et morphologie du milieu influen•ant l'•coulement de l'air. La nature et le niveau de pollution ne sont, par exemple, pas les mƒmes dans un parc ou en ville. Il est n•cessaire de prendre non seulement en consid•ration les mesures de polluants mais aussi toutes les informations renseignant sur la surface terrestre et l'atmosph…re. La surveillance de la qualit• de l'air doit ƒtre coupl•e avec une description du milieu urbain et la perception de la morphologie urbaine [START_REF] Basly | ERS SAR imagery for urban climate studies[END_REF](Basly et al. , 1999(Basly et al. , 2000)). Une telle approche est d•j' amorc•e dans les travaux pr•c•dents couplant observations satellitaires et observations par les stations de mesures [START_REF] Wald | Observing air quality over the city of Nantes by means of Landsat thermal infrared data[END_REF], Basly et Wald 2000, Basly 2000[START_REF] Ung | Air pollution mapping over a city -virtual stations and morphological indicators[END_REF]. Elle vise ' augmenter et spatialiser le nombre de donn•es utilis•es pour la cartographie des polluants. Nous savons que la pr•sence des polluants modifie la composition de l'atmosph…re et par la mƒme ses propri•t•s optiques [START_REF] Mclellan | Remote sensing of atmospheric turbidity variations by satellite[END_REF][START_REF] Wald | Observing air quality over the city of Nantes by means of Landsat thermal infrared data[END_REF]. Basly (2000), en admettant que la pr•sence de particules noires en suspension soit la principale cause de modification de la valeur num•rique des pixels dans la bande thermique

Etape 3: construction des stations virtuelles ' l'aide d'images satellitaires

Nous appelons ‡ stations virtuelles ˆ le sous-ensemble des pseudostations sur lesquelles une estimation de la pollution est possible ' l'aide des images satellitaires et des observations de concentrations de polluants. L'estimation de la pollution ' partir des images satellitaires n•cessite l'•tablissement d'une loi reliant les mesures satellitaires ' la transmittance de l'atmosph…re puis aux concentrations des polluants. En premi…re approximation, la r•flectance mesur•e au niveau du satellite s'•crit : ρ sat = ρ atm + ρ surf * τ atm avec ρatm la r•flectance de l'atmosph…re ρsuf la r•flectance de la surface terrestre. τatm la transmittance de l'atmosph…re En principe, il est n•cessaire de s•parer la contribution de la surface du sol de la contribution de l'atmosph…re et d'isoler la transmittance de l'atmosph…re. En r•alit•, cette s•paration est difficile ' r•aliser. Dans certains cas, les propri•t•s optiques de l'atmosph…re peuvent ƒtre •valu•es et leurs effets mod•lis•s [START_REF] Tanr | Description of a computer code to simulate the satellite signal in the solar spectrum[END_REF][START_REF] Roujean | A bi-directional reflectance model of the Earth's surface for the correction of remote sensing data[END_REF]). On peut ainsi estimer puis •liminer les influences de l'atmosph…re sur la r

Application à la ville de Strasbourg

Etape 2 : localisation des pseudostations ' l'aide des cartes d'identit€

La localisation des pixels ayant une mƒme signature spectrale est effectu•e station de mesures par station de mesures sur l'ensemble des images. Une analyse en composante principale est effectu•e sur toutes les bandes spectrales pour l'ensemble des images. Les donn•es sont projet•es sur les trois axes principaux, dont la somme des variances repr•sente 99 % de la variance totale. La classification des pixels semblables aux pixels p…res prend en compte la dynamique du signal projet•. Un pixel est class• comme •tant semblable ' un pixel p…re si ses trois composantes sont similaires ' celles du pixel-p…re, prises une ' une. Une erreur de 20 % de la dynamique de chacune des trois composantes est permise. Ainsi, pour une date de prise de vue donn•e et ' partir de onze stations de mesures (onze pixels p…res), nous obtenons onze classes de pixels semblables dont la population varie suivant la classe. La position de ces pixels est ensuite projet•e dans la g•om•trie de la base topographique. La figure 3 est l'union de tous les pixels localis•s pour les huit dates.

Parmi ces pixels de signature spectrale identique ' la station STG Centre 2, certains sont des pseudostations, d'autres sont des art•facts dus ' l'utilisation d'une carte d'identit• non compl…te. En effet, la signature spectrale enregistr•e au niveau du capteur est un m•lange de toutes les signatures spectrales des objets pr•sents (v•hicule, arbre, …) dans la sc…ne. Ces objets ont une signature spectrale variable dans le temps et ind•pendante les unes des autres. Afin d'•liminer les art•facts et de d•terminer les pixels qui sont r•ellement des pseudostations, une intersection multidate des pixels est effectu•e en tenant compte de l'erreur due aux recalages des images. Les r•sultats de l'intersection sont des pixels pr•sentant une signature spectrale stable dans le temps et identique aux pixels p…res.

L'exemple pour la station STG Centre 2 est illustr• dans les figures suivantes. La figure 4a) est une superposition sur la BD TOPO de tous les pixels de signature spectrale identique ' la station STG Centre 2. La station STG Centre 2 est situ•e au centre de l'image (point rouge) et les diff•rences de couleurs uniformes correspondent aux huit diff•rentes dates. La figure 4b) est aussi une superposition des pixels mais ' un autre lieu de la ville. Les intersections des pixels sont repr•sent•es dans les figures 4c) et 4d) et correspondent respectivement ' une station de mesures et une pseudostation localis•e et stable dans le temps.

Nous obtenons ainsi vingt-huit pseudostations, ce qui est peu par rapport aux nombres de pixels de signature spectrale identique aux pixels p…res, mais largement sup•rieur aux onze stations de mesures. Cette proposition de pseudostations sera test•e durant la campagne de juin 2002. Une visite de chacun de ces sites sera effectu•e. L'environnement g•ographique sera compar• ' celui de la station de mesures correspondante (pixel-p…re) afin de valider la similarit•. Un traitement sera ' nouveau effectu• en prenant en compte les crit…res morphologiques ainsi que le cadastre d'•mission et les statistiques de pollution. Les diff•rences entre les deux propositions de carte de pseudostations seront analys•es. Par ailleurs, des camions mobiles de mesures positionn•s ' ces sites mesureront de mani…re continue les concentrations des polluants en juin 2002. On pourra ainsi analyser de mani…re comparative les mesures des camions et des stations fixes et valider l'approche. Basly L., Wald L., 2000. Remote sensing and air quality in urban area. In Proceedings of TeleGeo'2000, the 2nd international workshop on telegeoprocessing, Sophia Antipolis, France, May 10-12, 2000, pp. 213-219. Basly L., Cauneau F., Ranchin T., Wald L., 2000. SAR imagery in urban area. In: Proceedings of the 19th Symposium of EARSeL : Remote Sensing in the 21st Century, Casanova J-L. ed., Valladolid, Spain, 31 May -2 June 1999, pp. 563-568. Basly L., Baleynaud J.M., Wald L.,1999. Satellite data for the air pollution mapping. In : Proceedings of the 18th EARSeL Symposium on Operational Remote Sensing for Sustainable Development, Enschede, Netherlands, 11-14 May 1998, pp. 133-139.

Liste des illustrations

  pallier ' ce manque d'information, le pr•sent article pr•sente une m•thodologie pour la cartographie des concentrations des polluants bas•e sur une approche multi-sources. Apr…s une br…ve description des m•thodes actuelles de cartographie de la pollution, nous pr•sentons la zone d'•tude et les donn•es utilis•es. Dans le paragraphe 3, nous pr•sentons la d•marche g•n•rale qui a •t• mise en oeuvre. La m•thodologie est d•coup•e en trois •tapes et s'appuie sur la d•finition de ‡ carte d'identit• ˆ, de ‡ pseudostations ˆ et de ‡ stations virtuelles ˆ ' l'aide d'observation in situ, d'observations satellitaires et d'une base de donn•es g•ographiques. La notion de carte d'identit• d'un site est cr••e afin de faire une caract•risation aussi fine que possible de la ville. Les stations de mesures sont situ•es dans des lieux de la ville pr•sentant des caract•ristiques spatiales, structurales, morphologiques et des signatures spectrales qui leur sont propres. Des caract•ristiques semblables sont recherch•es dans la carte d'identit• des sites de la ville ' l'aide des donn•es disponibles et d•finissent les pseudostations. Une loi d'estimation est ensuite n•cessaire pour passer des pseudostations aux stations virtuelles avec une estimation de la concentration de polluants. Cette densification du nombre de points de mesures (r•els et virtuels) permet de calculer une nouvelle cartographie des concentrations des polluants. Enfin le paragraphe 4 est une application ' la ville de Strasbourg, illustrant partiellement chaque partie de la m•thodologie. Cette m•thodologie sera valid•e en juin 2002 par une campagne de mesures utilisant des moyens de mesures mobiles.

  . A l'•chelle locale, les simulations num•riques tiennent compte de la morphologie locale mais la zone d'•tude est limit•e ' celle d'une rue par les ressources de calcul n•cessaires. A l'•chelle urbaine, les mod…les se distinguent suivant la r•solution spatiale de travail, les param…tres d'entr•es utilis•s, l'initialisation du mod…le, les conditions aux fronti…res, les techniques num•riques utilis•es. La qualit• des r•sultats de chacun des mod…les d•pend de la qualit• des param…tres d'entr•es. Vingt et une institutions faisant appel ' vingt-quatre mod…les diff•rents ont r•cemment particip• au ‡ Podbi-Exercice ˆ, un exercice de mod•lisation de la pollution atmosph•rique. Les r•sultats montrent qu'il n'existe aujourd'hui aucune proc•dure standard pour ce type d'exercice et que les r•sultats diff…rent d'un facteur 4 suivant les mod…les et les param…tres d'entr•es utilis•s

  . Le principal probl…me dans ces m•thodes concerne le manque de donn•es de mesures fournies par les r•seaux d'observation de la qualit• de l'air. Les r•sultats sont d•pendants du nombre et de la r•partition spatiale des mesures, ainsi que de la m•thode d'interpolation. Quelques recherches ont •t• men•es r•cemment pour cartographier la pollution atmosph•rique par t•l•d•tection. Les donn•es satellitaires sont utilis•es afin d'extraire les param…tres n•cessaires aux diff•rents mod…les de transport des polluants : r•seau routier, occupation du sol, rugosit•, etc. Par ailleurs, quelques auteurs

  Les donn•es de t•l•d•tection utilis•es pour notre •tude sont principalement des images Thematic Mapper (TM) du satellite Landsat 5. Huit prises de vue ont •t• acquises les 31 mars 1998, 15 ao"t 1998, 10 septembre 1999, 8 juin 2000, 1 avril 2001, 22 juillet 2001 et 14 ao"t 2001. Chaque prise de vue est constitu•e de sept canaux correspondants ' sept bandes spectrales allant du domaine visible du spectre •lectromagn•tique ' celui de l'infrarouge thermique. La r•solution spatiale des canaux est de 30 m ' l'exception du canal thermique TM6 qui a une r•solution spatiale de 120 m. Nous disposons •galement d'une image SPOT panchromatique qui offre une r•solution spatiale de 10 m. Cette image a •t• soigneusement recal•e avec la base de donn•es g•ographiques pr•sent•e ci-dessous. Chaque image Landsat a ensuite •t• recal•e par rapport ' cette image. Toutefois, seuls les mod…les de d•formation g

  par diffusion et absorption, d•termine une relation lin•aire liant le contenu en particules ' la mesure satellitaire. Appliqu•e ' une image enti…re, cette relation a permis d'•tablir une premi…re cartographie de la distribution spatiale des polluants pour la date •tudi•e. Les structures apparaissant dans l'image finale sont identiques ' celles de l'image source, et ne refl…tent certainement pas la situation r•elle d'un •pisode de pollution. Le domaine de validit• de cette loi est donc restreint. La difficult• est de cerner ce domaine de validit• et d'isoler les sites pour lesquels les comptes num•riques sont r•ellement corr•l•s aux mesures de pollution.Nous proposons une approche multi-sources qui cherche ' identifier les sites pour lesquels les variations de luminance sont essentiellement dues ' la pr•sence de polluants, en dehors des variations optiques de l'atmosph…re ' m•so-•chelle (environ 100 km). La m•thodologie pr•sent•e ci-dessous s'articule autour des notions de ‡carte d'identit• ˆ, de ‡pseudostation ˆ et de ‡station virtuelle ˆ qui seront d•finies ult•rieurement. Elle est d•coup•e en trois •tapes : 1) caract•risation de la ville et construction des cartes d'identit•, 2) localisation des pseudostations ' l'aide des cartes d'identit•, et 3) construction des stations virtuelles ' l'aide des images satellitaires.Etape 1: caractérisation de la ville et construction des cartes d'identitéCette premi…re •tape a pour objectif de se doter de moyens pour la compr•hension et l'•valuation de l'environnement urbain en tant que facteurs agissant sur la pollution atmosph•rique. L'espace urbain est un milieu complexe constitu• ' proportions variables d'espaces construits, de voies de communication, de sols nus, d'espaces r•sidentiels, industriels, herbac•s et bois•s, etc. La pollution atmosph•rique n'est •videmment pas la mƒme dans chacun de ces espaces (Derbez 2001). Des •tudes ont montr• l'h•t•rog•n•it• de la pollution atmosph•rique. La position des bOEtiments et l'orientation des rues par rapport au vent influencent la dispersion des polluants[START_REF] Lamoureux | Etude de l'effet du vent sur la diffusion d'un polluant[END_REF]. Des concentrations en polluants tr…s diff•rentes peuvent ƒtre observ•es sur des rues voisines(Scaperdad et Colville 1999[START_REF] Croxford | Sitting consideration for urban pollution monitors[END_REF][START_REF] Croxford | Spatial distribution of urban pollution: civilizing urban traffic[END_REF]. C'est pourquoi il est n•cessaire de faire une caract•risation aussi fine que possible de la zone •tudi•e. Les informations n•cessaires peuvent ƒtre obtenues par des m•thodes d'analyse et traitement des donn•es, d'images et de la morphologie urbaine. Ces informations peuvent ƒtre de sources diverses, quantitatives ou qualitatives et d•finissent la carte d'identit• d'un •l•ment de la ville. Cette derni…re int…gre des variables en relation avec :  la position g•ographique du lieu consid•r•,  le type d'occupation du sol,  la proximit• par rapport aux sources d'•missions, notamment le trafic automobile,  la morphologie des bOEtiments entourant la zone,  les conditions climatiques et m•t•orologiques du lieu consid•r•. Chaque •l•ment de la ville, ainsi que les stations de mesures, est identifi• par une carte d'identit•. Il est alors possible d'•tudier les relations qui pourraient exister entre ces variables et la pollution aux niveaux des stations de mesures. Etape 2: localisation des pseudostations à l'aide des cartes d'identité L'implantation des stations de mesures est r•alis•e en fonction des objectifs de surveillance, de leur environnement imm•diat, des densit•s de population et des sources d'•missions de polluants. En raison du co"t d'une station de mesures, les associations en charge de la qualit• de l'air ne disposent pas plus de trois stations de mesures en ville et par classe de stations. Cependant l'urbanisation croissante des villes et le gigantisme de la surface des grandes villes laissent pr•sager l'existence de sites similaires aux stations de mesures. On se propose de densifier virtuellement le r•seau en d•tectant des lieux de la ville oŽ une station de mesures aurait pu ƒtre implant•e, station similaire ' l'une de celles d•j' implant•es. Il est fait l'hypoth…se qu'il sera possible de mod•liser la concentration des polluants en ces lieux ' partir d'images satellitaires et des relations •tablies entre ces images et les stations de mesures. La localisation de ces lieux doit prendre en compte les conditions d'implantations des stations, les objectifs de surveillance, l'environnement imm•diat, les densit•s de population et les sources d'•missions de polluants. Nous d•finissons les pseudostations comme •tant des lieux de la ville pr•sentant des attributs similaires aux stations de mesures, la notion de similitude •tant d•finie ' l'aide des cartes d'identit•. L'identification de ces pseudostations rel…ve alors d'une classification multi-vari•e des donn•es. A une station de mesures est alors associ•e une classe de pseudostations (fig. 1).

  •flectance r•fl•chie ' m•soou subm•so-•chelle (de 10 ' 100 km). Le principe retenu dans notre approche, est que la diff•rence de r•flectance d'une pseudostation par rapport ' la r•flectance de la station de la mƒme classe refl…te, d'une mani…re ou d'une autre, une variation d'•paisseur optique et donc de concentration de polluants par rapport ' cette station. Les difficult•s dans cette •tape consistent ' •liminer les effets purement atmosph•riques, g•n•ralement sans connaissance externe, puis ' estimer un mod…le convertissant des variations de r•flectance en concentration. Le travail dans cette •tape est peu avanc• dans son expression g•n•rale. Il est vraisemblable que le principe ne pourra s'appliquer r•ellement qu'' un sous-ensemble des pseudostations, d'oŽ cette notion de station virtuelle. Le faible nombre de mesures conduisant ' l'utilisation de plusieurs images acquises ' des instants diff•rents, voire par des capteurs diff•rents, en tant que mesure palliative, induit une difficult• suppl•mentaire. En effet, tous les termes de l'•quation pr•c•dente d•pendent du temps et de la longueur d'onde (et bien s"r, de l'espace).

Etape 1 :

 1 caract€risation de la ville et construction des cartes d'identit€La caract•risation de la ville et la construction des cartes d'identit•s sont effectu•es ' l'aide de trois sources de donn•es : les stations de mesures, la base de donn•es g•ographiques et les images satellitaires. L'exploitation combin•e de ces donn•es n•cessite un r•f•rentiel g•ographique commun ' toutes les donn•es. Ce r•f•rentiel peut ƒtre soit un r•f•rentiel cartographique, soit directement une image de r•f•rence. Nous avons choisi pour notre •tude la BD TOPO de l'IGN. Un positionnement des stations sur la BD TOPO est effectu• en relatif sur le terrain. Selon l'appr•ciation de l'op•rateur, l'erreur effectu•e est approximativement de 5 m…tres. Ce positionnement permet de croiser les donn•es de pollution avec les donn•es de la BD TOPO et permet de d•river un ensemble d'informations renseignant sur la morphologique du site[START_REF] Chaignaud | Qualit• de l'air en milieu urbain: extraction d'indicateurs morphologiques[END_REF][START_REF] Durrenberger | Indicateurs morphologiques et pollution de l'air dans la CUS[END_REF][START_REF] Schnell | Mise en place et evaluation d'une base de donn•es sur l'emprise de la Communaut• Urbaine de Strasbourg, rapport de stage[END_REF] :-pourcentage de bOEti dans un secteur de rayon d•fini (fig.2a) -distance par rapport au premier bOEtiment dans une direction fix•e (fig.2b) -secteur de visibilit• dans une direction fix•e (fig. 2c) -surface ou volume de visibilit• en 2 ou 3 dimensions avec un rayon d•fini (fig. 2d) -indice de circularit• de la surface : indice de Miller (fig. 2e) -indice d'•talement de la surface : indice de Morton (fig. 2f) -indice de compacit• de la surface : indice de Grav•lius (fig. 2g) L'ensemble de ces informations est int•gr• dans la carte d'identit• et constitue une premi…re caract•risation des stations de mesures et de la ville. Quelques r•sultats sont pr•sent•s (tab. 1) pour les stations de mesures de l'ASPA et permettent par exemple de tirer les conclusions suivantes sur les stations de mesures. Les stations STG Amp…re, Rhin, et Reichstett •tant situ•es sur le toit d'un bOEtiment r•pertori• dans la base g•ographique, pr•sentent des indices morphologiques en 2D nuls. La station STG Hoerdt est localis•e ' l'ext•rieur de la zone couverte par la base g•ographique, le calcul des indices n'est donc pas possible. Pour les sept autres stations, les aires de visibilit• sont comprises entre 7000 et 22400 m 2 environ. Le maximum est atteint pour la station STG Ouest qui pr•sente effectivement une ouverture plus grande ' l'•coulement d'air, en accord avec la r•alit• du terrain. Quant au minimum, il est atteint pour la station STG Nord. La station est effectivement situ•e dans une configuration rectangulaire de bOEti identique ' STG Centre et Est, mais la proximit• de l'ensemble des bOEtiments ' la station STG Nord explique ce minimum. La station STG Cl•menceau pr•sente la deuxi…me plus grande aire de visibilit•. Sa localisation sur la place Cl•menceau, donc ' proximit• des grands axes routiers, offre une configuration de bOEti qui n'est plus rectangulaire mais ouverte vers les grands axes routiers. Tab. 1: les résultats sur le positionnement et les indices morphologiques Ajout•es ' ces indices morphologiques, des donn•es satellitaires viennent compl•ter les param…tres de la carte d'identit•. Un pr•traitement des images Landsat est effectu• pour les rendre superposable avec la base de donn•es g•ographiques et les donn•es des stations de mesures. Un recalage g•om•trique de chacune des images Landsat est effectu• avec le canal panchromatique (r•solution 10 m) d'une image SPOT ' l'aide d'un mod…le polynomial d'ordre 1. Cette image SPOT est elle-mƒme recal•e avec la BD TOPO'. Les r•sultats concernant les erreurs montrent que le recalage pr•sente une erreur quadratique de 0.6 pixels en colonne, de 0.6 pixels en ligne et de 1 pixel en circulaire. L'image de travail ayant une r•solution spatiale de 30 m, l'erreur est alors de 30 m et sera prise en compte dans la suite de notre •tude. Nous disposons d'une loi permettant d'associer ' un pixel de l'image SPOT, tous les pixels des images Landsat au pixel pr…s. Ce recalage des images Landsat avec la base de donn•es g•ographiques permet de localiser pr•cis•ment les pixels qui sont associ•s aux stations de mesures : les pixels ‡ p…res ˆ. La carte d'identit• d'une station de mesures contient alors : -la position de la station dans la BD TOPO ; -les donn•es relatives aux polluants mesur•s (concentration de polluants min, max, pic, …) ; -le pourcentage de bOEti dans un secteur de rayon d•fini ; -la distance par rapport au premier bOEtiment dans une direction fix•e ; -le secteur de visibilit• dans une direction fix•e ; -la surface ou le volume de visibilit• en 2 ou 3 dimensions avec un rayon d•fini ; -l'indice de circularit• de la surface de visibilit• : indice de Miller ; -l'indice d'•talement de la surface : indice de Morton ; -l'indice de compacit• de la surface : indice de Grav•lius ; -les contenus num•riques dans les sept bandes spectrales pour toutes les dates de prise de vue Landsat, -les indices d•riv•s des donn•es satellitaires. Notons que le calcul des indices morphologiques s'effectue ' l'aide du logiciel ESRI / ArcView et que le temps d'ex•cution n•cessaire ne nous a pas permis de disposer d'une carte compl…te de ces indices sur la CUS lors de cette application. De mƒme, l'exploitation du cadastre d'•mission et des statistiques d'observations de la pollution n'a pas pu ƒtre termin•e ' temps. Il en r•sulte que, pour cette application, nous nous limitons uniquement aux donn•es satellitaires provenant des huit prises de vue Landsat. Les donn•es relatives ' la morphologie de la ville et ' la pollution n'•tant pas prises en compte, la carte d'identit• utilis•e ici n'est pas compl…te et contient uniquement la signature spectrale des stations de mesures. Par manque d'exp•rience, nous avons d•cid• d'exploiter toutes les images disponibles simultan•ment. D'autres essais ont •t• effectu•s, par exemple, seulement avec les dates d'une mƒme saison, ou avec les images des diff•rentes saisons pour une mƒme ann•e. Ces essais n'ont pas abouti ' des r•sultats plus prometteurs et ne sont pas plus justifi•s scientifiquement. Il demeure que l'influence du choix des dates et du nombre d'images restera ' •valuer.

Etape 3 :

 3 construction des stations virtuelles ' l'aide d'images satellitaires Cette partie est une illustration de construction des stations virtuelles ' l'aide d'images satellitaires appliqu•e ' la ville de Strasbourg et pour une date donn•e. Il ne s'agit pas d'une d•monstration, qui ne pourra ƒtre effectu•e qu'une fois l'approche des pseudostations valid•e. La figure 5 montre une carte de concentration en PM 10 pour la date du 31 mars 1998. Elle est obtenue par interpolation ' l'aide de la m•thode des ‡ plaques minces ˆ (Duchon 1976, Le M•haut• 1986) des trois donn•es de concentration en PM 10 fournies par le r•seau local de mesures (points noirs). La zone couverte par la carte est de 25,8*33,6 km 2 avec une r•solution spatiale de 30*30 m 2 . Pour une meilleure visualisation des sites pollu•s, le canal TM4 Landsat est incrust• en fond d'image, r•v•lant la structure des rues, les autoroutes et le Rhin. A cause du faible nombre de points de mesures disponibles, la carte de pollution obtenue est uniforme et ne repr•sente certainement pas la r•alit• de la situation de pollution pour la date •tudi•e. La figure 6 est une carte de pollution pour la mƒme date obtenue avec la notion de pseudostations. Trois cent une pseudostation sont localis•es ici en utilisant uniquement les signatures spectrales des pixels pour la date •tudi•e. En supposant que la loi qui relie les comptes num•riques TM6 aux concentrations de PM 10 est lin•aire, les stations virtuelles sont ensuite construites ' partir des pseudostations en estimant les niveaux de concentrations en PM 10 . Elles sont distribu•es dans tout l'espace. En comparaison avec la figure 5, la carte obtenue par interpolation des donn•es estim•es offre une repr•sentation plus r•aliste de la distribution des polluants en ville. Cependant, ce r•alisme n'est pas exactitude. Notre regard rep…re des structures li•es ' celles de la ville, mais il est pr•matur• de pr•tendre que la carte est exacte. En effectuant une op•ration similaire sur la ville de Nantes, Basly (2000) a trouv• une erreur relative d'estimation de la concentration de l'ordre de 100 %. Cette erreur repr•sente une am•lioration tr…s importante par rapport aux techniques actuellement employ•es mais reste extrƒmement •lev•e. L'analyse des images satellitaires acquises en juin 2002 permettra de confronter les comptes num•riques observ•s sur les pseudostations aux mesures des stations et des camions. On pourra tout d'abord d•montrer l'existence d'une relation puis tenter de mod•liser cette relation de mani…re suffisamment simple pour envisager une exploitation routini…re de la m•thodologie. 5 Conclusions et perspectives La strat•gie mise en place pour une cartographie des champs de polluants est complexe et d•coule des travaux publi•s ant•rieurement. Elle se compose de trois •tapes dont les principes ont •t• exprim•s mais dont les •tats d'avancement et de r•alisation diff…rent. Pour le moment, aucune d'entre elles n'a •t• valid•e. La premi…re •tape, intitul•e 'caract•risation de la ville et construction des cartes d'identit•' a permis effectivement d'int•grer les param…tres pertinents (indice de visibilit•, hauteur du bOEti …) pour la cartographie des polluants et ceci ' travers la d•finition de carte d'identit• d'un site de la ville. Cette carte d'identit• est construite ' l'aide de traitements multivari•s de bases de donn•es. Les comptes num•riques associ•s aux pixels de l'image contenant une station de mesures sont extraits. Des indicateurs morphologiques sont construits pour compl•ter la carte d'identit•. La carte d'identit• contient ainsi des param…tres qui sont en relation avec la pollution, les propri•t•s spectrales ou la nature du sol et les propri•t•s morphologiques du site. D'autres bases de donn•es, tel que le cadastre des •missions, seront int•gr•es afin d'adapter la m•thodologie pour une meilleure caract•risation de la ville. Nous avons d•fini les pseudostations comme •tant des sites pr•sentant des attributs semblables aux sites contenant les stations de mesures. La deuxi…me •tape est une •tape de classification des cartes d'identit• en vue de localiser ces pseudostations en int•grant des contraintes de stabilit• dans le temps. Une premi…re classification de ces cartes d'identit• est effectu•e pour l'application sur la ville de Strasbourg ' l'aide uniquement de donn•es satellitaires. Elle a fourni vingt-huit pseudostations localis•es sur la CUS. Des camions laboratoires seront envoy•s sur place lors de la campagne de juin 2002 dans l'optique de mesurer les concentrations en polluants afin de mieux appr•hender la r•alit• de ces stations. D'autres m•thodes de classification moins contraignantes feront l'objet de nos prochaines •tudes : elles pourraient fournir une plus grande densit• de pseudostations. Les avantages et les inconv•nients de diff•rentes m•thodes de classification seront discut•s.La troisi…me •tape concerne les relations existants entre les donn•es satellitaires utilis•es et les concentrations en polluants. L'identification de ces relations permet de construire des stations virtuelles ' partir des pseudostations. Une premi…re estimation de la concentration en particules est pr•sent•e ici en appliquant la mƒme m•thode queBasly (2000). Cependant ' travers divers calculs de corr•lation et de mod•lisation des voiles de polluants dans l'atmosph…re, de plus amples •tudes sont n•cessaires pour mettre en •vidence et mod•liser les relations entre les taux de polluants mesur•s au sol et les mesures de luminances effectu•es par les satellites.En d•finitive, l'approche d•velopp•e est un couplage synergique de toutes les sources d'information aboutissant ' la cr•ation des stations virtuelles. L'int•gration de ces points suppl•mentaires dans les m•thodes existantes de surveillance de la pollution permettra une cartographie plus fine des champs de concentration des polluants.BibliographieADEME et EDM, 2000. Classification et crit…res d'implantation des stations de surveillance de la qualit• de l'air, recommandation du groupe de travail ‡ caract•risation des sites ˆ, 70 pages. Basly L., 2000. T•l•d•tection pour la qualit• de l'air en milieu urbain. Th…se de doctorat, sciences de technologies de l'information et de la communication, Universit• de Nice Sophia Antipolis, France. 182 p.

Fig. 1 :

 1 Fig. 1 : Localisation des pseudostations par classification des cartes d'identit€.

Fig. 2 :Fig. 5 :

 25 Fig. 2 : Indicateurs et indices renseignant sur la morphologie du site

Fig. 2 :

 2 Fig. 2 : Indicateurs et indices renseignant sur la morphologie du site

gie de fusion de ces diff•rentes sources d'information pr•sent•e dans cet article utilise ces diff•rentes informations afin de mettre ' profit leur synergie et d'obtenir une meilleure cartographie de la pollution. 2 La ville d'€tude et les donn€es utilis€es

  propose une m•thode d'isolation des effets optiques de la pollution bas•e sur le principe de contamination des valeurs minimales radiom•triques trouv•es sur une image. Les r•sultats de Sifakis (2000) montrent que les trois premiers canaux Landsat (bleu, vert, rouge) dans le spectre du visible sont les plus sensibles ' des variations radiom•triques li•es ' la pr•sence de pollution. Ces travaux concourrent ' d•montrer les potentialit•s de l'imagerie satellitaire, notamment des images Landsat Thematic Mapper (TM), pour la cartographie de la pollution atmosph•rique. Cette discussion montre qu'il n'existe pas aujourd'hui de proc•dure standard permettant de conna‹tre la distribution des concentrations en polluants avec une r•solution spatiale et une couverture temporelle suffisante en milieu urbain. N•anmoins, les informations obtenues par chacune de ces m•thodes sont compl•mentaires et n•cessaires ' l'observation et la compr•hension du ph•nom…ne de pollution :

	 les mesures effectu•es par les stations de mesures permettent une surveillance en temps r•el de la
	pollution,
	 les mod…les num•riques prennent en compte les conditions m•t•orologiques et permettent une analyse et
	pr•diction de la pollution atmosph•rique ' grande •chelle spatiale,
	 ' petite •chelle, les simulations num•riques nous renseignent sur l'influence du bOEti sur l'•coulement de
	l'air et sur les ph•nom…nes de stockage des polluants.
	 enfin une cartographie de la pollution atmosph•rique ' l'aide d'images satellitaires permet de spatialiser les
	informations n•cessaires ' la compr•hension du ph•nom…ne.
	La strat•La zone d'•tude est la Communaut• Urbaine de Strasbourg (CUS). R•partie sur un territoire de 306 km 2
	(distance maximale du nord au sud : 28 km et d'est en ouest : 16 km) elle est s•par•e de l'Allemagne par le Rhin.
	Elle repr•sente 45 % de la population du Bas-Rhin alors qu'elle ne couvre que 6 % de sa superficie. La ville de
	Strasbourg est ' une altitude de 150 m…tres (au pied de la cath•drale) et se situe ' 48,33• latitude nord et 7,38•
	longitude ouest. Une description compl…te de la ville de Strasbourg (tourisme, •conomie, transport, patrimoine,
	tradition) peut ƒtre trouv•e sur le site Internet de la mairie de Strasbourg : http://www.mairie-strasbourg.fr . Les
	donn•es utilis•es pour la cartographie de la pollution dans la CUS sont des observations de la qualit• de l'air sur
	la ville de Strasbourg, des observations satellitaires et une base de donn•es g•ographiques BD TOPO de l'IGN.