John Douglas

Thomas Ulrich

Caterina Negulescu

Risk-targeted seismic design maps for mainland France

Keywords: seismic hazard maps, earthquake risk, France, probabilistic seismic hazard assessment, fragility curves, risk-targeting, collapse rate, acceptable risk

In this article the recently proposed approach known as 'risk-targeting' for the development of national seismic design maps is investigated for mainland France. Risktargeting leads to ground-motion maps that, if used for design purposes, would lead to a uniform level of risk nationally. The Eurocode 8 design loads currently in force for France are used as the basis of this study. Because risk-targeting requires various choices on, for example, the level of acceptable risk to be made a priori and these choices are not solely engineering decisions but involve input from decision makers we undertake sensitivity tests to study their influence. It is found that, in contrast to applications of this methodology for US cities, risk-targeting does not lead to large modifications with respect to the national seismic hazard map nor to changes in the relative ranking of cities with respect to their design ground motions. This is because the hazard curves for French cities are almost parallel. In addition, we find that using a target annual collapse probability of about 10 -5 for seismically-designed buildings and a probability of collapse when subjected to the design PGA of 10 -5 leads to reasonable results. This is again in contrast to US studies that have adopted much higher values for both these probabilities.

Introduction

Like almost all current seismic design maps, the map underlying the French National Annex of the Eurocode 8 (EC8) design code is based on a 'constant hazard' assumption, i.e. the accelerations for design are those for a return period (or alternatively an annual exceedance probability) that is constant across the national territory. As with the majority of design maps [START_REF] Bommer | Adapting earthquake actions in Eurocode 8 for performancebased seismic design[END_REF] the return period for the French National Annex is thought to be 475 years (corresponding to a 10% probability of exceedance in 50 years), although the return period is not explicitly stated in the document. As discussed by [START_REF] Luco | Risk-targeted versus current seismic design maps for the conterminous United States[END_REF], however, it would be more consistent with the final use of seismic design maps to adopt a 'constant risk' assumption in which the design ground motions are defined to lead to a certain level of risk, e.g. annual probability of collapse. This approach was used by the United States Geological Survey (e.g. [START_REF] Luco | Preparation of new seismic design maps for building codes[END_REF] to develop the so-called Risk-Targeted Maximum Considered Earthquake (MCE R ) ground motions in the ASCE Standard 7-10 and the 2012 edition of the International Building Code (International Code Council, Inc., 2012). Elsewhere the 'risk-targeting' philosophy is called a 'performance-goal based' (United States Nuclear Regulatory Commission, 2007), 'risk-informed' [START_REF] Konno | A developing risk-informed design basis earthquake ground motion methodology[END_REF] or 'risk-consistent' (Hadjian, 2002) approach. To our knowledge, this idea has not been investigated for the development of national design maps outside the USA. Bommer et al (2005a) propose a more sophisticated approach in which the seismic design levels are fixed based on iterative loss estimation having defined the loss target, but they do not present any calculations using their method. Within the nuclear industry, analyses like that used to develop risk-targeted ground-motion maps for standard structures are becoming widely used, as shown by [START_REF] Kennedy | Performance-goal based (risk informed) approach for establishing the SSE site specific response spectrum for future nuclear power plants[END_REF] in relation to the ASCE Standard 43-05. The threshold risk level for nuclear installations, however, is generally lower than targeted for standard structures for the new generation of US seismic design maps.

The aim of this article is to present tests of the risk-targeting procedure of [START_REF] Luco | Risk-targeted versus current seismic design maps for the conterminous United States[END_REF] for mainland France (i.e. excluding French overseas departments and territories) using the hazard model developed by [START_REF] Martin | Revision du zonage sismique de la France: Etude probabiliste[END_REF] and the subsequent map that provides peak ground accelerations (PGA) to define the demand spectrum of the French seismic design code. We investigate whether the current PGA recommendations would require large changes if a risk-targeted approach was adopted. This study is concerned with the definition of earthquake loading for new structures rather than to define the loading for retrofitting or risk analysis of existing buildings. The following section briefly introduces: the French seismic design code, inputs to the subsequent analysis and the method used here to compute earthquake risk. The subsequent section presents the results of the sensitivity tests in which the free input parameters in the analysis are varied. Next, an example map using this approach is presented and compared to the current map. The article ends with a discussion of the results, some conclusions and suggestions for future research.

Computing risk-targeted ground motions

In 2002 GEO-TER (Géologie Tectonique Environnement et Risques) published a probabilistic seismic hazard assessment (PSHA) for France based on a classic Cornell-McGuire approach [START_REF] Cornell | Engineering seismic risk analysis[END_REF]McGuire, 1976) with logic trees (e.g. [START_REF] Kulkarni | Assessment of confidence intervals for results of seismic hazard analysis[END_REF]Bommer et al, 2005b) to account for epistemic uncertainty [START_REF] Martin | Revision du zonage sismique de la France: Etude probabiliste[END_REF]. They presented maps for return periods of 100, 475, 975 and 1975 years. The PGA map for a return period of 475 years was used as the basis for a new seismic zonation following the work of the Groupe d' Études et de Propositions pour la Prévention du risque sismique en France (GEPP) mainly during the period [START_REF] Grünthal | European Macroseismic Scale 1998, vol 19. Cahiers du Centre Européen de Géodynamique et de Séismologie, Luxembourg Hadjian AH[END_REF][START_REF] Duckett | Risk analysis and the acceptable probability of failure[END_REF](e.g. LeBrun et al, 2004)). This zonation has recently (1st May 2011) come into application as part of the French National Annex of EC8, the European seismic design code. The main difference between the PSHA maps of GEO-TER and the new French zonation is that the zonation classifies locations into five categories depending on the PGA for 475 years, rather than the continuous values of PGA provided by the original PSHA map. The five classes are: 'très faible' (very low) (PGA < 0.7 m/s 2 ), 'faible' (low) (0.7 ≤ PGA < 1.1 m/s 2 ), 'modéré' (moderate) (1.1 ≤ PGA < 1.6 m/s 2 ), 'moyen' (medium)

(1.6 ≤ PGA < 3.0 m/s 2 ) and 'fort' (high) (PGA ≥ 3.0 m/s 2 )1 (Ministère de l' Écologie, de l' Énergie, du Développement durable et de l'Aménagement du territoire, 2005). The study of [START_REF] Martin | Revision du zonage sismique de la France: Etude probabiliste[END_REF] and the current seismic zonation are the hazard inputs to this article.

As the maps provided by [START_REF] Martin | Revision du zonage sismique de la France: Etude probabiliste[END_REF] are not sufficient for our purposes, we had to recalculate the hazard using their input parameters. The results of our recalculation are close to those of [START_REF] Martin | Revision du zonage sismique de la France: Etude probabiliste[END_REF], but not exactly the same, because we did not consider the smoothing branch of the logic tree (weight 20%) and because we used the 2003 version of the CRISIS PSHA software (rather than the version from 1999 used by GEO-TER).

Two performance levels are considered in EC8: 'no-collapse requirement', which should be met for a reference seismic action with a 10% probability of exceedance in 50 years (475year return period), and 'damage limitation requirement', which should be met for a seismic action with 50% probability of exceedance in 50 years (95-year return period). 'No-collapse' means that the structure should: withstand the seismic design action without collapse, retain its integrity and have residual load-bearing capacity after the earthquake. 'Damage limitation' means that the structure should withstand a more frequent seismic action without damage or limitations of use. The probability of collapse associated with designing for the 'no-collapse requirement' is not explicitly given in EC8 [based on the analysis of [START_REF] Luco | Risk-targeted versus current seismic design maps for the conterminous United States[END_REF] the new US code ASCE Standard 7-10 does assign a probability to this requirement].

In addition, the seismic actions of EC8 are expressed in terms of the PGA for a return period of 475 years (a gR ), which is used to scale the normalized spectral shape, rather than pseudo-spectral accelerations for 2475-year return period used by current American codes.

The PGA for the 'damage limitation requirement' is given as 0.4a gR . For structures of high importance (e.g. hospitals) an importance factor is used to increase the design PGA, which is equivalent to lengthening the return period of the design motions.

As shown by [START_REF] Kennedy | Performance-goal based (risk informed) approach for establishing the SSE site specific response spectrum for future nuclear power plants[END_REF], for example, the seismic risk, y(a 0 ), can be obtained by convolving the seismic hazard curve and the fragility curve using either of these equations:

y(a 0 ) = - ∞ 0 P a 0 (a)
dH(a) da da (1)

y(a 0 ) = ∞ 0 H(a) dP a 0 (a) da da ( 2
)
where P a 0 (a) is the fragility curve, i.e. the conditional probability of failure given the ground motion level a, for a design level of a 0 and H(a) is the seismic hazard curve, i.e. the exceedance frequency corresponding to ground motion level a. In the following, we consider the seismic risk y(a 0 ) to correspond to the annual probability of collapse of a given structure, as also considered by [START_REF] Luco | Risk-targeted versus current seismic design maps for the conterminous United States[END_REF]. As part of this study, we have computed y(a 0 ) using both equations and, as expected, they give the same result. We prefer Equation 2, as did [START_REF] Luco | Risk-targeted versus current seismic design maps for the conterminous United States[END_REF], which uses the derivative of the fragility curve, rather than Equation 1, using the derivative of the hazard curve. This is because the hazard is generally only computed in PSHA for a limited number of ground-motion (e.g. PGA) levels and, therefore, its derivative is less accurate than the derivative of a fragility curve, which is typically defined by an analytical equation. It is common to use the lognormal distribution, which is completely defined by its mean µ and standard deviation β , for fragility curves, i.e. P(a) = Φ{[ln(a)ln(µ)]/β }. [START_REF] Luco | Risk-targeted versus current seismic design maps for the conterminous United States[END_REF] use this equation for their analysis.

Using Equations 1 or 2, the hazard curve from PSHA and estimates of µ and β , the annual probability of collapse can be readily computed. The evaluation of this integral is the basis of risk targeting.

To derive a risk-targeted design map, firstly the hazard curves for each site are derived by PSHA in the standard way. Next the desired annual probability of collapse (Y ), which defines the acceptable risk threshold, is chosen. The design acceleration is then determined iteratively. At step i, the fragility curve for typical buildings is defined by the trial design PGA (a i ), conditional probability of collapse (X) for a i and the standard deviation β (a lognormal fragility curve is fully defined by the relation P a 0 (a 0 ) = X if the standard deviation is fixed). For this value of a i the risk integral is computed for the considered site to find y(a i ). Then, examining the value of Yy(a i ), a i+1 is generated for the next iteration. The procedure is then repeated until (a n ) converges. This iterative method is used for every location covered by the map. The risk integral takes into account the whole hazard curve rather than simply basing the design ground motions on the PGA for a given return period.

Consequently, the slopes of the hazard curves for each site have a large impact on the results: the relative ranking in design PGA for sites with parallel hazard curves would not change and only when hazard curves have different slopes could their ranking change, by adopting a risk-targeted approach.

The elegance of using knowledge of the level of risk to define the design ground motions comes at the cost of introducing various input parameters in place of a return period for a 'constant hazard' map. The evaluation of this integral requires that the acceptable risk to the population Y be quantified, which is obviously not a solely scientific question as compared to quantifying a constant hazard return period. Both are partly social questions, which might require input from politicians and sociologists, for example. The integral also requires the fragility curve for a building designed to the design code with, if assuming a lognormal distribution, its parameters µ (or alternatively X) and β , which are civil engineering questions.

In the documentation accompanying the French seismic regulations there is almost no guidance on what value to choose Y , or on the intrinsic values of X and β . The fragility curve (or curves) should, ideally, be based on analyses of typical structures designed using the EC8 regulations. Because the seismic design maps apply to all types of structures covered by the design code, an average curve should be used that takes account of the variability in the vulnerability of structures of different materials (e.g. reinforced concrete and steel) and geometries (e.g. different numbers of storeys and bays).

To develop risk-targeted maps for the USA, [START_REF] Luco | Risk-targeted versus current seismic design maps for the conterminous United States[END_REF] used fragility curves based on analyses of typical structures designed using the ASCE Standard 7-05 (Federal Emergency Management Agency, 2009). These analyses indicated that there is roughly a 10% chance that any structure will experience partial or total collapse as a result of its mapped design ground motion [START_REF] Luco | Risk-targeted versus current seismic design maps for the conterminous United States[END_REF]. Consequently [START_REF] Luco | Risk-targeted versus current seismic design maps for the conterminous United States[END_REF] used the equivalent of P(PGA) = 0.1. [START_REF] Luco | Risk-targeted versus current seismic design maps for the conterminous United States[END_REF] adopted a standard deviation β for their fragility curves of 0.8 and Chapter 21 of the ASCE Standard 7-10 specifies 0.6 for its sitespecific ground-motion procedures for seismic design.

EC8 does not provide nor discuss fragility curves for structures designed using this building code. The development of such curves for the purposes of developing risk-targeted ground motions requires the design of many structures for different levels of ground motion as indicated by the National Annex and different geometries and materials. Subsequently the fragility of these structures could be evaluated using finite-element codes and various levels of shaking (e.g. Seyedi et al, 2010). To our knowledge such calculations have not been made and fragility curves for EC8-designed structures have not been published. Therefore, we base our fragility curves on those used by [START_REF] Luco | Risk-targeted versus current seismic design maps for the conterminous United States[END_REF] but we acknowledge the uncertainty in this choice. Because of the uncertainty in all input variables, in the following section we undertake a number of sensitivity analyses.

Sensitivity tests

A series of sensitivity tests are conducted for these French cities with increasing level of seismic hazard: Paris (lowest zone), Lyon, Nantes, Nice, and Grenoble and Lourdes (highest zone in mainland France) (Figure 1), to assess the annual probabilities of collapse, y, resulting from the hazard map currently defined in the French National Annex of Eurocode 8. The seismic hazard curves for the six selected cities using the hazard model of [START_REF] Martin | Revision du zonage sismique de la France: Etude probabiliste[END_REF] are displayed in Figure 2. Due to the lack of fragility curves specifically for structures designed using Eurocode 8, various values of X and standard deviation of curves following the lognormal distribution were used. The results of these analyses are summarized in Table 1. From this table the sensitivity of y to the parameters X and β , for the 475-year design PGA, can be seen. This sensitivity is more clearly seen in Figure 3 where y is plotted against X for various values of β and two example French cities. In the following sections we try to choose appropriate values of the desired annual probability of collapse (Y ) based on our calculations, the available literature and engineering arguments.

Fragility curves

The value of β has a large effect on the fragility curves, since we are constraining the probability of collapse X for a design PGA with a 475-year return period, which should be a small value. Table 1 shows that when β is high [e.g. 0.8, the value adopted by [START_REF] Luco | Risk-targeted versus current seismic design maps for the conterminous United States[END_REF] when constraining for a design ground motion with a return period of roughly 2475 years] and X is low the corresponding fragility curve is unrealistic since it predicts low collapse probabilities even for very high shaking levels (high values of µ). Based on this table and the values of β reported by [START_REF] Crowley | D3.1 -Fragility functions for common RC building types in Europe[END_REF] for reinforced-concrete buildings classified into generic classes (e.g. 'mid-rise (4-to-7 storeys) seismically-designed'), which are around 0.5, we prefer to use 0.5 for this parameter, in contrast to 0.8 used by [START_REF] Luco | Risk-targeted versus current seismic design maps for the conterminous United States[END_REF] but relatively consistent with the 0.6 specified in the ASCE Standard 7-10. Even using β = 0.5 the values of µ for the 'medium' zone are still possibly unrealistically high when X is small.

This suggests that a high value of X (e.g. 10 -4 ) should be used. Figure 3 shows that the annual probability of collapse y and X are linked (if only small changes in the design ground motions based on a 475-year return period are desired) and, therefore, if X was large y would also need to increase to a unrealistically high value (see next section).

Target annual probability of collapse, Y

The choice of a target annual probability of collapse, Y , is discussed by [START_REF] Luco | Risk-targeted versus current seismic design maps for the conterminous United States[END_REF] for the USA. Analogous to what we have done, they computed the probabilities of collapse resulting from the design maps in ASCE Standard 7-05 using the methodology presented above, P(PGA) = 0.1, β = 0.8 and the current seismic hazard curves. They obtained an average value for western USA of 1% in 50 years (Y = 2 × 10 -4 ), which they then adopted as the risk target nationwide. [START_REF] Fajfar | A practice-oriented estimation of the failure probability of building structures[END_REF] and [START_REF] Goulet | Evaluation of the seismic performance of a code-conforming reinforced-concrete frame building -from seismic hazard to collapse safety and economic losses[END_REF] evaluate the failure probabilities of code-designed structures using two different approaches and find similar probabilities to that used by [START_REF] Luco | Risk-targeted versus current seismic design maps for the conterminous United States[END_REF].

Labbé (2010) tries to estimate the annual probability of collapse in mainland France with inputs from the French historical earthquake catalogue (BRGM-IRSN-EDF, 2011). Using a simple approach, he computes the mean areas affected annually by intensity higher than V on the European Macroseismic Scale (EMS98) [START_REF] Grünthal | European Macroseismic Scale 1998, vol 19. Cahiers du Centre Européen de Géodynamique et de Séismologie, Luxembourg Hadjian AH[END_REF]. From this he deduces the average probability of damage of grades 2 and 3 (moderate to significant damage) using the EMS98, which are about 10 -4 for grade 2 and about 10 -5 for grade 3. Partial collapse, which is of interest here, corresponds to grade 4, and can be estimated at about 10 -6 but caution should be used when interpreting this value because of the simplified methodology of [START_REF] Labbé | PSHA outputs versus historical seismicity: Example of France[END_REF]. [START_REF] Kennedy | Performance-goal based (risk informed) approach for establishing the SSE site specific response spectrum for future nuclear power plants[END_REF] computes the minimum seismic core damage frequency (SCDF) achieved at 28 nuclear power plants in the central and eastern USA when the site specific design response spectra for safe shutdown earthquakes are defined by state-of-the-art approaches.

The minimum SCDF lies in the range of less than 6 × 10 -6 to 0.6 × 10 -6 per annum with a median of about 3 × 10 -6 per annum. Based partly on the analysis of [START_REF] Kennedy | Performance-goal based (risk informed) approach for establishing the SSE site specific response spectrum for future nuclear power plants[END_REF], the target SCDF in United States Nuclear Regulatory Commission ( 2007) is 10 -5 . These values are of a similar order to those suggested for normal buildings by Labbé (2010) (and adopted here, see below), even though it could be argued that the consequences of an accident at a nuclear reactor are much higher than that associated with the collapse of a regular structure and, therefore, a lower risk target should be sought for nuclear facilities.

Based on a review of the literature and design codes, [START_REF] Duckett | Risk analysis and the acceptable probability of failure[END_REF] recommends that for accidental loads the engineer should aim for a risk of failure of roughly 10 -6 per annum but following a cost-benefit analysis and after consultation with the relevant authorities this could be increased up to 10 -4 per annum. He notes, however, that the public may accept higher risk levels for accidental loading, which cannot by definition be forecast, than for loading, such as earthquakes, that could be expected.

Based on simple reasoning, we believe the value used for risk-targeting in the USA is much too high. For example, assuming a 1% chance in 50 years and 1 million buildings (a rough number for areas within the 'moderate' and 'medium' hazard zones of France) we would have on average 200 collapsed buildings due to earthquakes per year, which seems much too high. As a compromise between roughly 10 -4 adopted by [START_REF] Luco | Risk-targeted versus current seismic design maps for the conterminous United States[END_REF] and computed by [START_REF] Fajfar | A practice-oriented estimation of the failure probability of building structures[END_REF], and Goulet et al ( 2007) and 10 -6 suggested by Labbé (2010), [START_REF] Duckett | Risk analysis and the acceptable probability of failure[END_REF] and [START_REF] Kennedy | Performance-goal based (risk informed) approach for establishing the SSE site specific response spectrum for future nuclear power plants[END_REF], we adopt Y = 10 -5 for subsequent analysis.

This value also means that the fragility curves for the selected value of β = 0.5 do not imply unrealistically strong buildings (i.e. the values of µ in Table 1 are quite small).

Probability of collapse at design PGA, X

The design PGAs used in Table 1 to define the fragility curves have been taken from the French seismic zoning, and are consequently based on ground motions computed for a 475year return period. Hence, we should not use the same value for X as the one used by [START_REF] Luco | Risk-targeted versus current seismic design maps for the conterminous United States[END_REF], which is based on accelerations at around a 2475-year return period, but instead we should use a lower value. After assuming β = 0.5 and Y = 10 -5 , Figure 3 and Table 1 show that X must be defined as around 10 -5 if large changes in the design PGAs (relative to 475-year accelerations) are not desired. Therefore, X = 10 -5 has been used here for subsequent analyses. We assume X = 10 -5 for a return period of 475 years because this leads to only minor modifications in the design PGAs for most locations in mainland France. This assumption means that application of the design levels defined by Eurocode 8 for France would lead to buildings with very small chances of collapse when subjected to their design loads.

Effect of using acceleration bands

Using ranges of acceleration rather than the exact PGA values for a 475-year return period computed by [START_REF] Martin | Revision du zonage sismique de la France: Etude probabiliste[END_REF] has a significant impact on the results. Table 3 illustrates that using the upper and lower limits of the EC8 acceleration levels can change the level of risk by one order of magnitude. Therefore, for the ground-motion maps associated with future design codes, banding perhaps should not be used but the original PGAs presented. This is especially true if a risk-targeting approach is used because the use of bands would reduce some of the benefits of the method.

4 Risk-targeted seismic design maps for mainland France

The risk-targeting approach can produce risk-constant acceleration maps, by back-calculating iteratively the accelerations constraining the fragility functions, after having fixed X, Y and β . Hence, if the annual probability of collapse computed with the trial design acceleration is higher than the targeted value, we have to increase the design acceleration (thereby, enforcing stronger structures) and vice versa. Table 4 summarizes the results obtained using this iterative procedure for the six cities.

It can be noticed that the accelerations computed with this approach are close to those on the current hazard map. This is what we expected because the slopes of the hazard curves (see values of k in Table 1) are quite similar over the national territory and the value of X was assumed such that it does not modify greatly the current design PGAs when Y = 10 -5

and β = 0.5. The PGAs from the risk-targeting approach are generally slightly lower than those given on the current hazard map since they generally correspond to a return period shorter than 475 years. As an example, the seismic design map of mainland France derived by targeting a risk of Y = 10 -5 with X = 10 -5 and β = 0.5 is shown in Figure 4.

[Fig. 4 about here.]

Comparing Figures 1 (the current seismic zonation) and 4 shows that the effect of risktargeting on the design PGAs is limited and the relative ranking of French locations would not greatly alter by adopting a risk-targeted approach. The ratio between these two maps is shown in Figure 5 -ASCE Standard 7-10 calls this ratio the risk coefficient. Locations of the highest hazard (e.g. Lourdes) would see their design PGAs reduced and the PGAs in areas with the lowest hazard (e.g. Paris) would increase. If a similar procedure to that followed in producing the current French National Annex, i.e. classifying locations into a handful of broad acceleration categories, was adopted for risk-targeted maps then most of the apparent differences between the maps shown on Figures 1 and4 would disappear. This is encouraging because it is important that the hazard levels do not change greatly from one generation of building code to the next so that buildings do not need significant retrofitting and that the loading conditions used by engineers are stable in time.

[Fig. 5 about here.]

k is roughly uniform across mainland France (about 4). This compares to 3-4 for western North America and to 1-2 for eastern North America (e.g. [START_REF] Ellingwood | Quantifying and communicating uncertainty in seismic risk assessment[END_REF] and a value of 3 recommended in the French National Annex of Eurocode 8 to modify the PGA for a 475-year return period to a different return period. This uniformity in k between locations in mainland France means that the effect of risk targeting on the relative ranking of French cities is limited. This contrasts with the findings for the western and eastern US, where the risk-targeted approach has been adopted to avoid surprising results when a uniform (and long) return period is used [START_REF] Luco | Risk-targeted versus current seismic design maps for the conterminous United States[END_REF].

Discussion

Rather than using the risk of collapse to define the acceptable risk level it may be more appropriate for countries of moderate seismicity, such as mainland France, if the risk target is defined in terms of the probability of slight damage, for example. Because of the nature of fragility curves and the observed distribution of earthquake damage, targeting the risk of a certain level of damage (e.g. collapse) roughly implies a risk target in terms of other damage levels too (e.g. if one building collapses, one hundred others may suffer slight damage).

The cost of each damage level could be quantified and then the annual cost of earthquake damage computed to facilitate the choice of what risk is acceptable. This is similar to the idea of Bommer et al (2005a) but less sophisticated.

As shown above, the acceptable risk target for collapse can be bracketed by arguments based on the observed rate of collapse in earthquakes and by comparisons to acceptable risks for other industries, for example. However, the other required value for risk targeting, X, the probability of collapse when a building designed to the seismic design code is subjected to the mapped ground motion, is a more difficult parameter to constrain since it is associated with the conservatism of the design code, which is almost never quantified. As also shown above, though, if large changes in the design ground motions are not desired, then once Y is chosen (or X is quantified) the necessary assumption for X (or choice for Y ) follows almost automatically.

For this article we have focused on providing a single ground-motion map for all types of structures covered by the EC8 design code (e.g. reinforced concrete and steel structures of various geometries), which is the current approach for expressing seismic loading in the majority of design codes. This means that for risk-targeting the fragility curves used must be sufficiently generic to cover the vulnerability of a large range of different buildings and β must capture this spread of vulnerabilities. This is in contrast to the fragility curves used for risk analysis of nuclear reactors, for example, where β is quite small since it does not need to account for the variability in response due to different building types or geometries. If future design codes adopted an approach where multiple ground-motion maps for different types of structures were acceptable, or their design procedures resulted in similar fragility curves for different types of structures, then more specific fragility curves could be used.

It also should be noted that this article is concerned with ground-motion maps to be used to define the earthquake loading for new structures rather than to define the loading for retrofitting of existing buildings or for risk assessments of such constructions. Like many European countries a large number of buildings date from before previous French seismic design codes and, therefore, they were not explicitly designed against earthquake loading, even if they will likely have some inherent resistance to minor shaking.

Conclusions

In this article we have applied the risk-targeting approach to the development of national seismic design maps for mainland France. We find that, by adopting reasonable values for the free input parameters, the maps produced are roughly in line with those currently in force. Before the risk-targeting approach can be routinely used, however, it is necessary to better constrain the free input parameters in this method. In particular, there is considerable uncertainty in appropriate fragility curves for a suite of buildings designed according to EC8, including the conditional probability of collapse when a building is subjected to its design load. Furthermore, what values should be used for the key parameter of desired annual probability of collapse is not simply a civil engineering issue. It should ideally involve input from the wider community (e.g. decision makers) as to what level of earthquake risk is acceptable.

Table 1 Annual probability of collapse of structures, y, based on the hazard map defined in the French National Annex of EC8. β is the standard deviation of the lognormal fragility curves used and X is the probability of collapse for the design PGA a gR , given as the range of PGA indicated on the map. PGA (GEO-TER 2002) is the PGA extracted from [START_REF] Martin | Revision du zonage sismique de la France: Etude probabiliste[END_REF], whereas PGA (recalculated) is the value recalculated for this study. k and k 0 are the shape parameter (slope) and constant, respectively, of the hazard curve when plotted on a log-log graph (e.g. [START_REF] Ellingwood | Quantifying and communicating uncertainty in seismic risk assessment[END_REF]. µ is the PGA corresponding to a probability of 0.5 of collapse for the fragility curve P a gR for the considered zone. The values of y[PGA(recalculated)] and X in the column 'Average for France' are the mean values for all locations in a grid covering mainland France. µ µ 10 -4 0.4 4.0 × 10 -5 5.7 × 10 -5 1.1 × 10 -4 3.4 × 10 -5 7.2 × 10 -5 7.6 × 10 -5 1.4 × 10 -4 1.9 7.8 " 0.5 3.0 × 10 -4 5.3 × 10 -5 1.1 × 10 -4 3.6 × 10 -5 7.4 × 10 -5 7.8 × 10 -5 1.4 × 10 -4 2.2 8.9 " 0.6 9.0 × 10 -4 4.9 × 10 -5 1.2 × 10 -4 3.8 × 10 -5 7.7 × 10 -5 8.0 × 10 -5 1.5 × 10 -4 2.6 10.4 " 0.8 2.5 × 10 -3 4.5 × 10 -5 1.2 × 10 -4 4.4 × 10 -5 8.6 × 10 -5 8.6 × 10 -5 1.7 × 10 -4 3.8 15.1 " 1 3.4 × 10 -3 4.1 × 10 -5 1.2 × 10 -4 4.8 × 10 -5 9.1 × 10 -5 8.7 × 10 -5 1.8 × 10 -4 6.0 24.0 10 -5 0.3 5.0 × 10 -13 7.9 × 10 -6 9.1 × 10 -6 1.9 × 10 -6 4.9 × 10 -6 3.1 × 10 -6 4.8 × 10 -6 3.3 13.2 " 0.4 8.0 × 10 -9 7.5 × 10 -6 9.5 × 10 -6 2.1 × 10 -6 5.3 × 10 -6 4.1 × 10 -6 6.6 × 10 -6 3.8 15.3 " 0.5 7.0 × 10 -7 7.0 × 10 -6 9.9 × 10 -6 2.4 × 10 -6 5.7 × 10 -6 5.0 × 10 -6 8.4 × 10 -6 4.5 17.9 " 0.6 8.0 × 10 -6 6.6 × 10 -6 1.0 × 10 -5 2.8 × 10 -6 6.2 × 10 -6 5.9 × 10 -6 1.0 × 10 2.0 × 10 -4 0.8 0.018 4.0 × 10 -5 3.9 × 10 -5 3.9 × 10 -5 4.0 × 10 -5 3.7 × 10 -5 4.2 × 10 -5 10 -5 0.4 5.0 × 10 -7 2.6 × 10 -6 5.2 × 10 -7 4.7 × 10 -7 4.2 × 10 -7 1.4 × 10 -7 5.3 × 10 -8 " 0.5 1.3 × 10 -5 2.6 × 10 -6 6.7 × 10 -7 6.2 × 10 -7 6.0 × 10 -7 3.1 × 10 -7 1.7 × 10 -7 " 0.6 7.0 × 10 -5 2.5 × 10 -6 8.2 × 10 -7 7.7 × 10 -7 7.6 × 10 -7 5.1 × 10 -7 3.4 × 10 -7 " 0.8 3.5 × 10 -4 2.4 × 10 -6 1.2 × 10 -6 1.2 × 10 -6 1.2 × 10 -6 9.8 × 10 -7 8.6 × 10 -7

Table 3 Influence of the use of ranges of acceleration instead of exact value for 475-year return period. The first part of this table shows the three acceleration levels used for constraining the fragility curve in m/s 2 : lower limit of the interval, exact value and upper limit. The second part shows y computed from these three acceleration levels. X = 10 -5 and β = 0.5 have been used. 2.4 × 10 -6 5.7 × 10 -6 5.0 × 10 -6 8.4 × 10 -6 1.0 × 10 -5

Paris

Exact 1.1 × 10 -5 2.1 × 10 -6 1.8 × 10 -6 1.7 × 10 -6 1.1 × 10 -6 4.9 × 10 -7 6.9 × 10 -5 a gR (upper limit) 1.0 × 10 -6 1.6 × 10 -6 4.4 × 10 -7 1.1 × 10 -6 1.8 × 10 -7 2.8 × 10 -7 1.5 × 10 -6

Table 4 Accelerations adjusted iteratively to obtain Y = 10 -5 with X = 10 -5 and β = 0.5. The accelerations from the French zoning map and the exact accelerations for a 475-year return period hazard are also given for comparison. All accelerations are expressed in m/s 2 . Also given are the return periods corresponding to these design PGAs.

Paris

[

  Fig. 1 about here.] [Fig. 2 about here.] [Fig. 3 about here.]

  It is interesting to compare the values of X and y computed for mainland France and for United States.[START_REF] Luco | Risk-targeted versus current seismic design maps for the conterminous United States[END_REF] based their fragility curves on pseudo-spectral accelerations (PSA) at 0.2 s computed for a 2475-year return period and P PSA(0.2s) [PSA(0.2 s)] = 0.1. Table 2 shows X(2475y) computed for various cities in mainland France and the corresponding y, for various values of β . Using the same β and targeting the same Y as[START_REF] Luco | Risk-targeted versus current seismic design maps for the conterminous United States[END_REF] (Y = 2 × 10 -4 and β = 0.8) leads to a lower value of X(2475y) (0.018) compared to Luco et al (2007) (0.1). The value that is assumed for X(2475y) is actually arbitrary because it simply defines which percentile of the fragility curve is fixed by the seismic design code.

Fig. 1 Fig. 2

 12 Fig. 1 French National Seismic Hazard Map (Ministère de l' Écologie, de l' Énergie, du Développement durable et de l'Aménagement du territoire, 2005) and the locations of the six cities chosen for this study.

Fig. 4

 4 Fig. 4 Seismic design map for mainland France derived by targeting a risk of collapse of Y = 10 -5 with X = 10 -5 and β = 0.5.

Fig. 5

 5 Fig. 5 Map of risk coefficients (as named by the ASCE Standard 7-10) defined by ratio of map shown in Figures 4 to current seismic zonation shown in Figure 1.

Table 2

 2 Annual probability of collapse of structures, y, based on PGA computed for a 2475-year return period using hazard curves computed by[START_REF] Martin | Revision du zonage sismique de la France: Etude probabiliste[END_REF].

	Average for France			Paris	Lyon	Nantes	Nice	Grenoble	Lourdes
	PGA(2475y)	( m/s 2 )		0.7	1.8	2.0	2.5	3.7	4.4
	y[PGA(2475y)]	β	X(2475y)			y[PGA(2475y)]
	3.5 × 10 -4	0.8	0.100	2.0 × 10 -4	2.6 × 10 -4	2.5 × 10 -4	2.7 × 10 -4	2.5 × 10 -4	3.2 × 10 -4

Only islands in the French Antilles (e.g. Guadeloupe and Martinique) are classified in the highest zone.

Acknowledgements

We thank Evelyne Foerster and Kushan Wijesundara for useful discussions on various aspects of this study and Nicolas Luco and Peter Fajfar for answering our questions on the targeted annual collapse rate. This article was supported by an internal BRGM research programme (VULNERISK/MULTIRISK 2011). We thank the anonymous reviewers for their comments on an earlier version of this article.

, β=1.0 Average, β=0.5 '' , β=0.6 '' , β=1.0 Fig. 3 y against X computed for two French cities and the average for France and various values of β (see Table 1).


Florence Naaim-Bouvet

Hervé Bellot

Kouichi Nishimura

Christophe Genthon

Cyril Palerme

Gilbert Guyomarc'

HOW TO DETECT SNOW FALL OCCURRENCE DURING BLOWING SNOW EVENT?

Keywords: Blowing snow, drifting snow, snow fall, precipitation, sensor, Alps, Antarctic

. In such way, it could be possible to evaluate precipitation. But depending on the height, the size of snow fall particle and the wind speed it is not so easy to distinguish between both types of events. In such case, the analysis of snow flux and mean diameter according to wind speed allows to separate blowing snow event with and without precipitation. It gives better results in predicting the timing of precipitation than done by precipitation gauges. A simpler photoelectric counter such as designed by Wenglor could also give some interesting results. These conclusions are supported by field campaigns conducted in French Alps and in Antarctica (coastal Adélie Land) during last winters.

INTRODUCTION

In mountainous areas, drifting snow influences the spatial distribution of the snow cover and consequently snowpack stability and avalanche danger. When comparing models with in situmeasurements, it is first necessary to identify snow-drift events with and without concurrent falling snow. In Antarctica, the mass balance of the Antarctic ice sheet is a key parameter of sea level rise. However, in the coastal areas where katabatic winds are strong and frequent, it is necessary to identify and separate blowing snow and precipitation. ______________________ * Corresponding author address: Florence Naaim-Bouvet, Snow Avalanches Engineering and Torrent Control research unit, IRSTEA, 2 rue de la papeterie, BP 76,email: florence.naaim@irstea.fr In mountainous areas but also in Arctic or Antarctic continent falling snow is rapidly redistributed by wind. In such conditions not only accurate measurement of precipitation but also detection of precipitation is a challenge. The blowing snow billow can be very high up to several hundred meters (Scarchilli et al., 2010). So the measured snow water equivalent suffers of high uncertainty whatever is the position of precipitation gauge. Present weather sensors such as VPF730 [START_REF] Bellot | Present weather sensor tests for measuring drifting snow[END_REF] are also unable, at the present time and without any additional treatment, to distinguish between falling and blowing snow. Empirical methods based on sensors (wind, precipitation and snow depth sensors, Webcam) were developed too but are not enough accurate (Vionnet et al., 2011). [START_REF] Nishimura | Blowing snow at Mizuho station, Antarctica[END_REF] showed that a Snow particle counter (SPC-S7) set at a specific height, preferably higher than 10 m, can be used to evaluate precipitation because in this case a bi-modal distribution could appear. But depending on the height, the size of snow fall particles and the wind speed it is not so easy to identify this bimodal distribution. That is why we turn to the concomitant analysis of snow flux, mean diameter and wind speed to detect blowing snow event with and without precipitation. First data stemming from SPC-S7 are analyzed. Then a simpler photoelectric counter designed by Wenglor Sensoric gmbh is also tested for this purpose. The results are supported by field campaigns conducted in French Alps and in Antarctica (coastal Adélie Land) during the last winters.

BLOWING SNOW EVENTS WITH AND WITHOUT SNOW FALL DETECTED BY A SNOW PARTICLES COUNTER (SPC-S7)

Snow particles Counter SPC-S7

The snow particle counter (SPC-S7, Niigata electric) (figure 1) is an optical device [START_REF] Nishimura | Blowing snow at Mizuho station, Antarctica[END_REF]. The diameter and the number of blowing snow particles are detected by their shadows on photosensitive semiconductors. Electric pulse signals of snow particles passing through a sampling area are sent to an analyzing logger. In this way the snow particle counter detects particles between 20 and 250 µm in mean radius size particle, divides them into 32 classes, and records the particle number every 1 s. Assuming spherical snow particles, the horizontal snow mass flux q is calculated as follows:

(1) Where q d is the horizontal snow mass flux for the diameter d, n d is the number of drifting snow particles, S the sample area, t the sample period, and ρ p the density of the drifting snow particles (917 kg m -3

).

Alps field experiment configuration

Lac Blanc Pass (French Alps) is located at the Alpe d'Huez ski resort near Grenoble, France. The large north-south-oriented pass, 2800 m high, consists of relatively flat terrain over a length of about 300 m where blowing snow has been studied for twenty years by IRSTEA (previously Cemagref) and Meteo France. This experimental site is presented in details in this proceeding [START_REF] Guyomarc | Measurement campaigns to investigate blowing snow and snow drift conditions at a high altitude site[END_REF]. The measurements presented here are done during winter 2010-2011 and 2011-2012. Three Snow Particles Counters (figure 3) were installed on mast. One SPC was set up at a fixed position (4.4 m above the ground). Two others were set up near the snow pack surface. A fixed distance of one meter separates them. These two SPC and two cup anometers set up at the same height could be raised manually when the snow depth increases and buries the sensors. A snow depth sensor (figure 2) that measures the exact position of the SPC above the snow pack and an ultrasonic anemometer supplements these devices.

Detection of snow fall during blowing snow event

During blowing snow event without precipitation it is known that the size distribution of snow particles at a specific height can be accurately approximated by the means of the two parameters gamma probability density function (in red in figures 4 and 5). As previously stated, during snow fall at a specific height, preferably higher than 10 m, we can observe a bi-modal distribution in Antarctica. But at Lac Blanc pass, it is not often the case even if blowing snow occurs with concurrent precipitation at least 37% of time [START_REF] Vionnet | Occurrence of blowing snow events at an alpine site over a 10year period: observations and modeling[END_REF]. It must be noticed that the highest SPC is set up at 4.4 m above the ground. Depending on the height, the wind speed and the snow particles size, the snowfall contributions can be hidden behind the large contribution of blowing snow particles (figure 5) and it is not always so easy to detect the bimodal distribution. That's why we choose to develop a method based on combined analysis of wind and drift sensors. distinguishing blowing snow events with and without snowfall. It is easier to detect in the highest layers (figures 6, 7, 8). Nevertheless the different types of events can be distinguished from 1, 2 m high. Exact transitions between events remain difficult to detect. A more detailed analysis of physical processes involved during blowing snow event without concurrent falling snow can be conducted and could improve the prediction.

-The snow particles diameter declines with height above the snow (figure 9) -Larger particles are being carried aloft at higher wind speed (figure 10) Previous example can be considered as a "textbook case" and sometimes the events are not easy to characterize (figure 11). In these case, the more detailed analysis based on snow particles diameter allow to conclude that snow falls occurred during this event. But precipitation gauges set up at Lac Blanc Pass did not detect any snow fall. Nevertheless other precipitation gauges recorded snow fall in a nearby area protected from wind at 1800 m a.s.l. and experimenters present at Lac Blanc Pass at this period confirmed the presence of snow fall. As conclusion this method gives better results in predicting the timing of precipitation than those obtained with precipitation gauges.

BLOWING SNOW EVENTS WITH AND WITHOUT SNOW FALL DETECTED BY A COMMERCIAL PARTICLE-COUNTING DEVICE (WENGLOR®)

3.1 The Wenglor® sensors (YH03PCT8 and YH08PCT8)

The Wenglor® sensors (YH03PCT8 and YH08PCT8), names Wenglor thereafter, consist of a laser (655 nm), photo sensor, and switching circuit (figure 12). When a particle passes through the 0.6 mm diameter, 30 mm (80 mm depending on the model) long laser beam, the sensor outputs a digital signal which is recorded on a data logger. This sensor is not capable of size discrimination but its main advantages are low cost and low power requirement which is essential in Antarctica.

As far as we know, Leonard (Leonard andCullather, 2008, Leonard et al., 2011) was one of the first researchers who proposes and tests the use of the Wenglor commercial particle-counters to quantify drifting snow.

Figure 12: Wenglor sensor deployed at Lac Blanc Pass

Alps field experiment configuration and results

The Wenglor and the SPC have been compared in a side-by-side deployment (figure 3) at Lac Blanc Pass in the winter of 2011-2012. SPC has a smaller particle size detection threshold than Wenglor. That's why it records far more particles. But despite this fact and for the studied blowing snow event (figure 13), there is a linear relationship between Wenglor particle count per 10 minutes and the average value of snow mass flux over 10 minutes calculated by the SPC-S7. In such way, the shapes of the curves linking wind speed to particle count and wind speed to snow mass flux are similar (figures 14 and 15). Similar analysis than those presented in section 2.3 can be done with Wenglor sensors. But we have to keep in mind that additional information concerning the snow particles diameter is no more available. Consequently it will be less accurate.

CONCLUSIONS AND PERSPECTIVES

The analysis of snow flux and mean diameter according to wind speed allows to separate blowing snow events with and without precipitation. It gives better results in predicting the timing of precipitation than those obtained with precipitation gauges. Even if it is less accurate a simpler photoelectric counter such as designed by Wenglor could also give some interesting results. Nevertheless this sensor still needs some tests and perhaps some improvements to be used in Antarctica under harsh conditions during several months. The next step will aim at evaluating precipitation [START_REF] Leonard | Snowfall measurements in the Amundsen and Bellingshausen Seas[END_REF].

Figure

  Figure 1: Snow Particles Counter

Figure 4 :

 4 Figure 4: Particle diameter distributions measured with the highest SPC (blowing snow event with snow fall -2011-02-23 4H00)

Figure 6 :Figure 7 :Figure 8 :

 678 Figure 6: Horizontal snow mass flux at 0.2 m high plotted against wind velocity from 23 to 26 February 2011

Figure 9 :

 9 Figure 9: Snow particles diameter plotted against height above surface during blowing snow events without falling snow from 23 to 26 February 2011

Figure 11 :

 11 Figure 11: Horizontal snow mass flux at 3.4 m high plotted against wind velocity from 20 to 21 February 2011

Figure 13 :

 13 Figure 13: Particle count measured by the Wenglor sensor at 3.4 m high plotted against snow mass flux from 10 to 11 April 2012

Figure 14 :

 14 Figure 14: Horizontal snow mass flux at approximately 1.2 m high plotted against wind

Proceedings, 2012 International Snow Science Workshop, Anchorage, Alaska


Lionel Bombrun
email: lionel.bombrun@ims-bordeaux.fr

Yannick Berthoumieu
email: yannick.berthoumieu@ims-bordeaux.fr

MULTIVARIATE TEXTURE RETRIEVAL USING THE KULLBACK-LEIBLER DIVERGENCE BETWEEN BIVARIATE GENERALIZED GAMMA TIMES AN UNIFORM DISTRIBUTION

Keywords: Texture, image retrieval, Kullback-Leibler divergence, Multivariate elliptical distribution

This paper presents a new multivariate elliptical distribution, namely the multivariate generalized Gamma times an Uniform (MGΓU) distribution. Because it generalizes the multivariate generalized Gaussian distribution (MGGD), the MGΓU distribution is able to fit a wider range of signals. For the bivariate case, we provide a closed-form of the Kullback-Leibler divergence (KLD). We propose the MGΓU distribution for modeling chrominance wavelet coefficients and exercise it in a texture retrieval experiment.

A comparative study between some multivariate models on the VisTex and Outex image database is conducted and reveals that the use of the MGΓU distribution of chromiance wavelet coefficient allows an indexing gain compared to other classical approaches such as MGGD and Copula based model).

INTRODUCTION

In the context of texture image recognition, the wavelet representation has been widely considered to characterize the texture. The multiscale recognition algorithm operates in two main steps. One is the feature extraction which consists in modeling each wavelet subband by a given probability density function (pdf). Its estimated parameters compose the signature of the texture. Then, a similarity measure based on a probabilistic metric is computed between the query image and all the other images of the database.

Many univariate models such as the generalized Gaussian [START_REF] Do | Wavelet-Based Texture Retrieval Using Generalized Gaussian Density and Kullback Leibler Distance[END_REF] and the Weibull [START_REF] Kwitt | Lightweight Probabilistic Texture Retrieval[END_REF] distributions have been proposed for modeling the magnitudes of wavelet detail coefficients. Those works have further been extended by the use of the generalized Gamma (GΓ) distribution which generalizes these two models [START_REF] Choy | Statistical Wavelet Subband Characterization Based on Generalized Gamma Density and its Application in Texture Retrieval[END_REF]. However, these approaches do not fully exploit the texture information in the image. Multivariate distributions have hence recently been introduced to model the spatial and/or color dependencies of the wavelet coefficients such as the Gaussian [START_REF] Stitou | Copulas based Multivariate Gamma Modeling for Texture Classification[END_REF] and Student-t [START_REF] Kwitt | A Joint Model of Complex Wavelet Coefficients for Texture Retrieval[END_REF] copulas. Other approaches based on the multivariate elliptical distributions have been proposed such as the multivariate generalized Gaussian distribution (MGGD) [START_REF] Verdoolaege | Multiscale Colour Texture Retrieval using the Geodesic Distance between Multivariate Generalized Gaussian Models[END_REF]. For this later one, the multivariate Laplace is a special case, and the multivariate G 0 distribution. For both models, a closed form of the geodesic distance (GD) [START_REF] Berkane | Geodesic Estimation in Elliptical Distributions[END_REF] has been established for fixed shape parameters and an approximation of the GD has been given for the general case by assuming the geodesic coordinate functions as straight lines [START_REF] Verdoolaege | Geodesics on the Manifold of Multivariate Generalized Gaussian Distributions with an Application to Multicomponent Texture Discrimination[END_REF] [START_REF] Bombrun | Mutlivariate Texture Retrieval Using the Geodesic Distance Between Elliptically Distributed Random Variables[END_REF]. Unfortunately, most of the time, only approximations of the GD can be derived for non trivial multivariate distributions (MGGD, G 0 , . . . ).

The main paper contribution is to propose a new multivariate model which generalizes the MGGD, namely the multivariate generalized Gamma times an Uniform (MGΓU) distribution. For the special case of a bivariate generalized Gamma times an Uniform (BGΓU) distribution, an analytical expression of the Kullback-Leibler divergence (KLD) [START_REF] Kullback | Information Theory and Statistics[END_REF] is derived. The paper is structured as follows. In section 2, we introduce the MGΓU distribution and its parameters estimation. Section 3 provides the derivation of the analytical expression of the KLD for the BGΓU distribution. The BGΓU distribution is then proposed for modeling the chrominance wavelet coefficients. Some indexing results on the VisTex and Outex databases are also presented. Conclusions are finally reported in Section 4.

THE MGΓU DISTRIBUTION

Definition

In this paper, we introduce a novel distribution, namely the MGΓU distribution, which belongs to the family of elliptically contoured distributions. This model admits the following stochastic representation [START_REF] Cambanis | On the Theory of Elliptically Contoured Distributions[END_REF]:

k = τ Σ 1 2 x, (1)
where x is an uniformly distributed random vector on the unit sphere R p , p being the dimension of vector k. Σ is the normalized scatter matrix (i.e. tr(Σ) = p). τ is a positive scalar random variable distributed according to a Generalized Gamma distribution, whose pdf p τ (τ ) is given by [START_REF] Stacy | A generalization of the gamma distribution[END_REF]:

p τ (τ |α, β, γ) = GΓ(τ |α, β, γ) = ατ αγ-1 e -( τ β ) α β αγ Γ(γ) , (2)
where α > 0, β > 0, γ > 0 and τ ∈ 0, +∞ . α, β and γ are respectively the power, scale an shape parameters of the distribution.

To derive the expression of the pdf of k = τ Σ 1 2 x, let first consider the scalar random variable u = k T Σ -1 k. As x is uniformly distributed on the unit sphere R p , it can be easily proved that u follows the same distribution as τ 2 . After the variable change u = τ 2 in (2), we obtain:

p u (u|α, β, γ) = GΓ u α 2 , β 2 , γ = αu αγ 2 -1 e -u α/2 β α 2β αγ Γ (γ) . (3)
By definition, k follows a multivariate elliptical distribution, then:

p k (k|α, β, γ, Σ) = 1 |Σ| 1 2 h p (k T Σ -1 k), (4)
where h p (•) is the density generator function. The change variable u = k T Σ -1 k in (4) leads to:

p u (u|α, β, γ) = π p 2 Γ(p/2) u p 2 -1 h p (u). (5)
By identification between (3) and ( 5), we obtain the analytical expression of h p (•) for the MGΓU distribution, it yields :

h p (u) = Γ(p/2) π p 2 α 2β αγ Γ (γ) u αγ-p 2 e -u α/2 β α . (6)
Finally, the expression of the pdf for an MGΓU distribution is obtained:

p k (k|α, β, γ, Σ) = 1 π p 2 |Σ| 1 2 αΓ(p/2) 2β αγ Γ(γ) k T Σ -1 k αγ-p 2 × e -( k T Σ -1 k ) α 2 β α . (7)
A particular case of the MGΓU distribution is the MGGD which is obtained when α = 2b, β = 2 1 2b and γ = p 2b [START_REF] Gómez | A Multivariate Generalization of the Power Exponential Family of Distributions[END_REF]. Note that the multivariate Gaussian distribution is obtained when α = 2, β = √ 2 and γ = p 2 , and the multivariate Laplace distribution is obtained when α = 1, β = 2 and γ = p.

Parameter estimation

Let (k 1 , . . . , k N ) be N vectors distributed according to an independent and identically distributed random MGΓU distribution. The expression of the maximum likelihood estimators (MLE) for the parameters of an MGΓU distribution can be easily derived by differentiating the joint distribution of (k 1 , . . . , k N ) with respect to α, β, γ and Σ. An iterative algorithm can be easily derived to solve the ML equations as in [START_REF] Bombrun | Performances of the Maximum Likelihood Estimators for the Parameters of Multivariate Generalized Gaussian Distributions[END_REF]. Unfortunately, this algorithm is very complex since the MLEs of the MGΓU distribution parameters are all dependent of each other. To circumvent this difficulty, we propose a simpler estimation algorithm based on a mixed moment and ML estimation procedure. Given the N vectors (k 1 , . . . , k N ), the normalized scatter matrix is estimated by using a momentbased approach, it yields:

Σ = 1 N N i=1 k i k T i with tr( Σ) = p. (8)
Then, the N scalar variables

u i = k T i Σ-1 k i are computed.
According to (3), the scalar random variable u follows a GΓ distribution of parameters α/2, β 2 and γ. Those parameters can hence be estimated by using the ML method.

TEXTURE INDEXING

In the following, only the bivariate case will be considered (i.e. p = 2). In this context, we can prove that the Kullback-Leibler divergence (KLD) can be analytically derived for a BGΓU distribution. We then propose this distribution as a suitable model for the chrominance wavelet coefficients in a texture retrieval experiment.

Kullback-Leibler divergence

After some cumbersome computations, for the bivariate case, one can derive the KLD between two populations distributed according to a BGΓU distributions of parameters

(α 1 , β 1 , γ 1 , Σ 1 ) and (α 2 , β 2 , γ 2 , Σ 2 ), it yields: KLD α1, β1, γ1, Σ1||α2, β2, γ2, Σ2 = ln α1 α2 + α2γ2 ln β2 -α1γ1 ln β1 + ln γ2 γ1 + 1 2 ln |Σ2| |Σ1| -γ1 + α1γ1 -α2γ2 2 2 ln β1 + 2 α1 Ψ(γ1) - α2γ2 2 -1 1 2 ln(g1g2) + 1 3 (g2 -g1) g2 2F1 1, 3 2 ; 5 2
; g2 -g1 g2

+ 1 3 (g1 -g2) g1 2F1 1, 3 2 ; 5 2
; g1 -g2 g1

+ β1 β2 α 2 Γ γ1 + α 2 α 1 2Γ(γ1) g α 2 /2 2 2F1 - α2 2 , 1 2 ; 3 2 ; g2 -g1 g2 + g α 2 /2 1 2F1 - α2 2 , 1 2 ; 3 2 ; g1 -g2 g1 , (9)
where 2 F 1 is the Gauss hypergeometric function and Ψ the digamma function. g i = 1 λi where λ i are the eigenvalues of Σ -1 1 Σ 2 , and λ 1 ≥ λ 2 . This expression of the KLD for the BGΓU distribution generalizes the KLD found for the BGGD [START_REF] Verdoolaege | Multiscale Colour Texture Retrieval using the Geodesic Distance between Multivariate Generalized Gaussian Models[END_REF], bivariate Laplace and bivariate Gaussian distributions when we replace the value of α, β and γ by their corresponding expressions given in Section 2.

Indexing results on the VisTex and Outex databases

As an analytical expression of the KLD exists for the BGΓU distribution, this new model is of potential interest for modeling the chrominance wavelet coefficients. An indexing experiment is hence conducted on the VisTex and Outex databases. For the VisTex database, 40 texture classes from the MIT Vision Texture database have been selected. From each of these texture images of size 512 × 512 pixels, 16 subimages (128 × 128) are created. Hence, the VisTex dataset contains 640 texture images [START_REF] Do | Wavelet-Based Texture Retrieval Using Generalized Gaussian Density and Kullback Leibler Distance[END_REF]. For the Outex database, the 319 color images in BMP format with 600dpi under inca lightning conditions from the Outex website have been used [START_REF]Outex Texture Database[END_REF]. Those images were then cropped to 128 × 128 pixels starting from the top-left hand corner, hence creating a dataset of 5109 images. In this experiment, the orthogonal wavelet decomposition (with 1 scale) with Daubechies' filter db4 is used for the wavelet decomposition.

To evaluate the potential of the BGΓU distribution for modeling the chrominance dependency of the wavelet coefficients, the empirical histogram of the scalar random variable u = k T Σ-1 k has been computed from the first scale and orientation of the grayscale Bark6 image of the VisTex database. This histogram is then modeled by the corresponding distributions of u (see ( 3)) for the BGΓU, the BGGD, the Laplace and the Gaussian distributions as shown in Fig. 1. As observed in Fig. 1, the BGΓU distribution allows a better fit of the scalar random variable u. As the BGΓU contains one more parameter compared to the BGGD, this new model is more flexible and seems promising for the modeling of the chrominance wavelet coefficients.

A texture retrieval experiment is now conducted on the VisTex database. The symmetric Kullback-Leibler divergence (SKLD) has been considered as similarity measure, i.e. half of the double-sided KLD. Moreover, as the subbands are assumed to be independent, the SKLD between two images I 1 and I 2 is obtained by computing the sum of the SKLDs evaluated for each subbands. Table 1 shows the average retrieval rate on the VisTex (top) and Outex (bottom) databases. The rows L and Ch means that an univariate generalized Gaussian (GG) distribution and a bivariate model have respectively been considered for the luminance and chrominance wavelets coefficients. The five proposed bivariate models are the Gaussian, Laplace, BGGD, BGΓU distributions and the Gaussian Copula with Weibull distributed margins. The row Ch + L shows the indexing results when combining the chrominance and luminance models. As the chrominance and luminance images can be considered as independent, the SKLD computed between two images I 1 and I 2 can be decomposed as follows: [START_REF] Kullback | Information Theory and Statistics[END_REF] where Ch i and L i are the chrominance and luminance images of I i (i = 1, 2). The line RGB corresponds to indexing results with a trivariate model for modeling the RGB wavelet coefficients dependencies. Note that for the trivariate case, only the Gaussian distribution and the Gaussian copula with Weibull distributed margins have been considered since their SKLD admits an analytical expression [START_REF] Lasmar | Gaussian Copula Multivariate Modeling for Texture Image Retrieval Using Wavelet Transforms[END_REF].

SKLD(I 1 ||I 2 ) = SKLD(Ch 1 ||Ch 2 ) + SKLD(L 1 ||L 2 ),
Note that here, indexing results on the Outex database seems to be very low, but this dataset is composed of very close texture classes (e.g. 11 classes of barleyrice). As shown in Table 1, the proposed BGΓU distribution seems promising for modeling the chrominance wavelet coefficients. A gain of 1% and 3% are respectively observed on the VisTex and Outex databases compared to the BGGD. Moreover, for both databases, combining the chrominance and luminance models (Ch+L) yields to a raise of the retrieval rate of about 10% compared to a separate indexing (Ch or L). For the Vis-Tex database, the combination of the BGΓU model for the chrominance and the univariate generalized Gaussian model for the luminance channel yields to the best indexing results. A gain of 2% is observed compared to a trivariate model applied directly on the RGB channels. Concerning the Outex, this gain is confirmed with respect to the trivariate Gaussian. If compared with the trivariate Gaussian copula, the indexing performance are approximately equivalent (32.39 vs. 33.15). Note also that the combination of chrominance and luminance models (Ch+L) allows parallelization of the code since the chrominance and luminance channels are treated independently.

Fig. 2 draws the recall/precision computed on the Vis-Tex database for the four best models in each category (L, Ch, Ch+L and RGB). In a retrieval context, the best model is the one which has the highest recall and precision values. This graph confirms that the best indexing results are obtained with a combination of the proposed BGΓU distribution (solid green line) for modeling chrominance wavelet coefficients and the univariate generalized Gamma distribution for luminance wavelet coefficients.

CONCLUSION

In this paper, a new stochastic model has been introduced i.e. the multivariate generalized Gamma times an Uniform (MGΓU) distribution. The multivariate Gaussian, Laplace and generalized Gaussian distribution (MGGD) can be considered as a special case. For bivariate data, a closed-form of the Kullback-Leibler divergence (KLD) has been established. The MGΓU distribution has hence be proposed for modeling chrominance wavelet coefficients, and implemented in a texture retrieval context. Some experiments on the VisTex and Outex databases have been conducted and revealed that the proposed MGΓU model allows a better retrieval rate compared to other classical approaches (MGGD, Copula based models).

Fig. 1 .

 1 Fig. 1. Modeling of the scalar random variable u for the first scale and orientation of the grayscale Bark6 image.

Fig. 2 .

 2 Fig. 2. Average retrieval rate as a function of the number of retrieved images for the VisTex database.

Table 1 .

 1 Average retrieval rate (in %) on the VisTex (top) and Outex (bottom) databases.

		Gaussian	Laplace	BGGD	BGΓU	Gaussian copula Weibull	GG
	L	-	-	-	-	-	69.54
	Ch	63.86	63.96	69.35	70.77	66.28	-
	Ch + L	77.73	79.46	80.48	81.46	78.48	-
	RGB	79.49	-	-	-	79.88	-
	L	-	-	-	-	-	21.79
	Ch	22.06	21.82	22.19	25.21	23.39	-
	Ch + L	30.70	30.76	29.80	32.39	31.33	-
	RGB	29.84	-	-	-	33.15	-