Camille Couprie
email: c.couprie@esiee.fr

Leo Grady
email: leo.grady@siemens.com

Laurent Najman
email: l.najman@esiee.fr

Hugues Talbot
email: h.talbot@esiee.fr

Anisotropic diffusion using power watersheds

Keywords: Combinatorial optimization, image processing, denoising, mathematical morphology, watersheds

come

INTRODUCTION

The most common assumption about image noise is that the noise is high frequency. Therefore, a natural choice of image filter is a lowpass filter which dampens power in the high-frequency range. Unfortunately, a lowpass filter has the undesirable effect of blurring object boundaries, since these boundaries have a high-frequency. This dilemma was recognized early in the image processing literature and various approaches were proposed to solve it, with the goal of smoothing the image content internal to an object, but preserving image discontinuities across boundaries.

An approach to discontinuity-preserving filtering was proposed by Perona and Malik (PM) [START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF] who modeled image filtering as an anisotropic diffusion process that smoothed image intensities inside an object while preserving the intensity discontinuity between objects. Later, the anisotropic diffusion approach was shown by Black et al. [START_REF] Black | Robust anisotropic diffusion[END_REF] to be interpreted as a gradient descent method for optimizing a robust error func-tion model which had been proposed previously. Anisotropic diffusion has been a very successful filtering algorithm, in part because it is easy to implement. However, this algorithm does have the problem of setting two parameters, the robust estimator parameter and the diffusion time. Practical usage of anisotropic diffusion requires a choice between long computation times or blurry boundaries.

Recently, Couprie et al. introduced the power watershed (PW) method [START_REF] Couprie | Power watersheds: A new image segmentation framework extending graph cuts, random walker and optimal spanning forest[END_REF]. Although this technique was introduced in the context of image segmentation, the method can be used as an optimization method for some functionals. In this paper we show that PW is well-suited to address the robust estimator filtering model and therefore provide an alternative to anisotropic diffusion for the optimization of this model. There are several advantages of PW as compared to anisotropic diffusion in the context of optimizing the robust estimator filtering model. First, we remove the robust estimator parameter. Second, we preserve discontinuities without blurring while also allowing for a fast optimization based on PW. Third, no time step needs to be determined (i.e., no risk of a divergent solution). Therefore, our optimization via PW removes the false choice presented by anisotropic diffusion between fast optimization and sharp object boundaries.

Other models for discontinuity-preserving filtering are also well-known, including total variation [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], the Mumford-Shah functional [START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF] and the piecewise constant Mumford-Shah (Chan-Vese) model [START_REF] Vese | A multiphase level set framework for image segmentation using the Mumford and Shah model[END_REF]. The total variation model is convex, allowing for a relatively efficient solution [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF][START_REF] Darbon | Image restoration with discrete constrained total variation part I: Fast and exact optimization[END_REF] while the Mumford-Shah and Chan-Vese models are more difficult to optimize. However, recent work [START_REF] Bresson | Fast global minimization of the active contour/snake model[END_REF][START_REF] Grady | The piecewise smooth Mumford-Shah functional on an arbitrary graph[END_REF] has provided more efficient optimization methods. Despite this progress in the optimization of other discontinuity-preserving filtering methods, very little work has been done to improve the speed/blur tradeoff that is necessary in the anisotropic diffusion algorithm. Rather than promoting one discontinuity-preserving model over another, we simply show here how the PW may be used to very efficiently optimize the robust estimator model underlying anisotropic diffusion such that we no longer have to choose between fast optimization and blurred boundaries. Our optimization will be achieved at the cost of a few iterations of the efficient PW, which is simple to implement.

FORMULATION

Black et al. showed that anisotropic diffusion could be viewed as the minimization of a robust estimator. We employ this same robust estimator formulation of anisotropic diffusion, but show that the power watershed may be used to iteratively optimize the formulation.

Since the power watershed is defined on a graph, we begin by casting the anisotropic diffusion algorithm in discrete terms. A graph consists of a pair G = (V, E) with vertices v ∈ V and edges e ∈ E ⊆ V × V . An edge, e, spanning two vertices, v i and v j , is denoted by e ij . A weighted graph assigns a (typically non-negative and real) value to each edge called a weight. The weight of an edge e ij is denoted by w ij . We define the edge-node incidence matrix A of a graph G as

A eij v k =    +1 if the node i = k, -1 if the node j = k, 0 otherwise.
The incidence matrix A is used here as the combinatorial analogue of the continuous gradient operator. Given these definitions, the anisotropic diffusion algorithm is defined as

dx dt = A T g(Ax)Ax, (1)
where x is the image intensities of a filtered image, and x at time 0 equals the input, unfiltered image. The function g(x) is used to prevent blurring over edges and may be of any decreasing form. One possibility for g(x) suggested in [START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF] is

g(x) = exp -αx 2 , ( 2
)
where α is a free parameter.

The typical method for solving the anisotropic diffusion equation in ( 1) is via a forward Euler method in which the iteration

x k+1 = x k + dtA T g(Ax k )Ax k ,
is applied, using a time step dt which is designed to obey the CFL conditions and ensure stability of the solution. In practice, the number of iterations applied to produce the filtered image is a free parameter which is set to define the desired level of smoothing. Black et al. [START_REF] Black | Robust anisotropic diffusion[END_REF] showed that the anisotropic diffusion equation given by (1) could be viewed as the gradient of the energy

E(x) = σ(Ax), (3)
where σ(x) is a robust estimator or an M-estimator. A gradient of the Black et al. energy in (3) is given by

dE dx = A T σ ′ (Ax)Ax. (4)
Therefore, when σ ′ (Ax) = g(Ax), then the gradient of the Black et al. energy, (4) is the same as the anisotropic diffusion equation in [START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF]. In other words, anisotropic diffusion is gradient descent minimization of (3). Further, the σ(z) corresponding to the PM weighting function in ( 2) is given by the Welsch function σ(z) = 1 -exp -αz 2 . Black et al. explored the effect of different robust estimator forms of σ(z). However, a common characteristic of all robust estimators is that they exhibit linear or quadratic growth near z = 0, but provide a nearly constant output as z → ∞. The transition between quadratic growth and constant growth is controlled by the parameter α. If α → ∞, then the robust estimator ultimately achieves its constant growth phase for any z = 0. In this way, the σ(z) function may be viewed as a differentiable approximation to the || • || 0 norm in which ||0|| 0 = 0 and ||z|| 0 = 1 for any z = 0. Consequently, the anisotropic diffusion filtering algorithm may be viewed as gradient descent on an energy functional which is an approximation to the || • || 0 norm.

Previous filtering algorithms that explicitly formulated filtering from the standpoint of optimizing the || • || 0 norm took the form

E(x) = σ(Ax) + λh(x, f ), (5)
in which f represents the intensities of the input unfiltered image, h(x, f ) represents a loss function and λ is a free parameter. When h(x, f

) = ||x -f || 2 2
, the gradient of (5) be-

comes dE dx = A T σ ′ (Ax)Ax + λ (f -x)
, which achieves a stable point when

A T σ ′ (Ax)A + λI x = λf. (6)
Since ( 6) may be viewed as a backward Euler solution for the anisotropic diffusion equation ( 1) when λ = 1 dt , then the anisotropic diffusion algorithm may be seen as the optimization of a robust estimator of the image gradient balanced against a loss function of the form ||x-f || 2 2 where the tradeoff between gradient smoothness and data fidelity is governed by λ = 1 dt . Therefore, by setting a fixed time for the anisotropic diffusion, we may view the solution obtained for this time as a steady-state optimization of our second energy functional [START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF] with a corresponding λ parameter. Viewing the time parameter as a loss function for the fidelity of the filtered image with the original, we may freely alter the loss function.

Anisotropic diffusion and L 0 norm

As said in the previous section, it is possible to alter the loss function. For example, we could optimize the fully robust energy

E(x) = σ(Ax) + λσ (x -f ) ,
with gradient given by

dE dx = 2αA T σ ′ (Ax)Ax + 2αλσ ′ (x -f )x.
If, at any iteration, we fix the values of x inside the robust error function, then we have

dE dx k+1 = A T σ ′ (Ax k )Ax k+1 + λσ ′ (x k -f )x k+1 . (7)
This energy may be written as the steady-state optimization of the energy functional

E k+1 = x k+1 T A T σ ′ (Ax k )Ax k+1 + λx k+1 T σ ′ (x k -f )x k+1 ,
that may also possibly be written

E k+1 = eij σ ′ (Ax k ) x k+1 i -x k+1 j 2 +λ vi σ ′ (x k -f ) x k+1 -f 2 . ( 8
)
This expression for the energy to compute a minimum step is of a form that may be optimized by the power watershed (PW) [START_REF] Couprie | Power watersheds: A new image segmentation framework extending graph cuts, random walker and optimal spanning forest[END_REF] if the parameter α → ∞ for g(x).

The generalized PW energy is given by

min x eij ∈E w ij p |x i -x j | q + vi∈V w i p |x i -y i | q (9)
where y represents a measured configuration and x represents the target configuration. In this equation, the first term essentially forces x to vary smoothly within an object, and the second term enforces data fidelity. We can see that Eq. 8 is of the same form as Eq. 9, with q = 2, p = α. Therefore, we suggest the following filtering algorithm : Use PW with y = f to optimize (8) to obtain

x k+1 . k = k + 1; until ||x k+1 -x k || 2 < ǫ
This technique can be initialized by the image f , or, in the case of very noisy images by a smoother version of the image, for instance by the application of a Gaussian or a median filter.

In the initialization step of the algorithm, we need to build a graph in linking each pixel node with its neighbors, for example its 4 neighbors for 4-connectivity. The values of those nodes will be the x k . Noting N the number of pixels in the image, we add N nodes to this graph in linking each pixel node v i by an edge to the new node v i+N . Those additional nodes are set to constant values of f . A description of the graph on a example is shown in Fig. 1. The algorithm using power watershed for anisotropic diffusion is summarized in Alg. 1.

RESULTS

We now demonstrate the performance of the power watershed algorithm for anisotropic diffusion in presenting results on synthetic and real images. A first, the test in Fig. 2 is performed on a synthetic image corrupted with a Gaussian noise of standard deviation σ = 16. For each result, we may compute the peak signal-to-noise ratio (PSNR) relatively to the original image. Typical values for the PSNR in denoising lie between 20 and 40 dB where higher is better. In comparison with PM algorithm, the power watershed algorithm tends to produce piecewise constant results. The result obtained at Fig. 2(d) shows a good compromise between noise removing and edge preservation. However there is a choice to make for PM algorithm (Fig. 2(b) and Fig. 2(c)) between complete denoising and good contrast edge restoration.

Examples on real images in Fig. 3 and Fig. 4 show that the power watershed algorithm may be useful as a filtering step before segmentation. In terms of computation time, each iteration of the power watershed algorithm for anisotropic diffusion operates in quasi-linear time in practice. We noticed empirically that the convergence is very fast (fewer than 10 iterations). The result obtained at the first iteration is generally close to the final result as shown in Fig. 4. We give some timing with necessary number of iterations to reach convergence at Table 1.

CONCLUSION

The power watershed (PW) algorithm was originally presented as a technique applied to image segmentation. In this work, we showed that PW provides an optimization procedure that allows us to optimize robust error measures which are operated in a particular parameter range. Black et al. showed how the anisotropic diffusion method could be viewed as optimization of a robust error filtering model. We showed that PW could be used to optimize the same robust error filtering model, leading to an alternative optimization procedure.

An aspect of the PW optimization in this context is that it applies to the optimization of the robust error filtering model when operated in a particular parameter range which effectively models the image as being piecewise constant. For cases in which this model was reasonable (such as the synthetic images), the algorithm exhibited a strong denoising capability. When this denoising algorithm was applied to real images, the effect of imposing a piecewise constant model was to effectively quantize images into a small number of grayscale levels (Fig. 3 and4). This transformation of an image into a piecewise constant form may be helpful as a preprocessing step for segmentation or recognition.

We used the energy minimization aspect of PW to expand the traditional use of watersheds from a segmentation algorithm to an image filtering algorithm. Although our algorithm optimizes the filtering objective as a real-valued optimization problem, we showed that we still retain the fast speed of the watershed algorithm. However, the limit of PW as an energy optimization strategy remains unclear -What energy functions can be minimized via power watersheds? Future work will address this issue and continue to demonstrate applications which can benefit from the speed of power watershed.

Algorithm 1 :

 1 Anisotropic diffusion using PW Data: An image f , an initial solution x 0 , λ ∈ R * + Result: A filtered image x k Set k = 0. Build the graph described in Fig. 1. repeat Generate the pairwise weights exp -(Ax k ) 2 , and unary weights exp -(x kf ) 2 .

Fig. 1 .

 1 Fig. 1. Graph necessary to build for applying the power watershed diffusion algorithm to a 3 × 3 image. The unary weights in blue and pairwise weights in red are updated at each step.

  (a) Noisy image, PSNR = 24.24dB (b) PM, PSNR = 34.03dB (c) PM, PSNR = 30.46dB (d) PW, PSNR = 31.54dB

Fig. 2 .

 2 Fig. 2. Comparison of Perona-Malik(PM), and power watershed(PW) algorithms for denoising a synthetic image. (b) PM used with 80 iterations α = 0.0015, leading to a good PSNR but with remaining isolated noisy pixels. (c) PM, best compromise found for this image to remove the isolated pixels with 50 iterations and α = 0.0005. (d) A median filtered image as initialization and λ = 0.975 allows to obtain a better PSNR while removing isolated noisy pixels.

Fig. 3 .

 3 Fig.3. Filtering of a liver image by power watershed. Here noise and small vessels are both removed, leading to a result which may be used as a first step before segmentation.

Fig. 4 .

 4 Fig. 4. The segmentation of concrete images in three classes (bubbles in black, stones in light grey) is useful for the study of the material's mechanical properties. The segmentation is obtained by two thresholds of the filtered image by PW.

Table 1 .

 1 Number of iterations of Perona-Malik, and power watershed algorithms on the image of Fig.2, Fig.4, and Fig.3

	Perona-Malik	PW


Guyomarc'

Gilbert Guyomarc'
email: gilbert.guyomarch@meteo.fr

Measurement campaigns to investigate blowing snow and snow drift conditions at a high altitude site

Keywords:

In order to make progress in avalanche research and to improve avalanche risk forecasting, at least two elements must be taken into account: the investigation of the spatial variability of snow depth and the study of meteorological conditions during blowing snow events in high alpine areas. Numerical models (NEMO, Sytron 2 & 3) have been developed to simulate the redistribution of the snow pack during those events. More recently, we tested the ability of a meso-scale atmospheric model to explicitly simulate wind-induced snow transport in alpine terrain. The characteristics and mechanisms of blowing snow events must be verified and the results of the models can then be validated through the use of in-situ measurements. Over the last two years we have successfully performed intensive measurement campaigns at our experimental site (Col du Lac Blanc, French Alps, 2,700 m asl). We followed the evolution of snow depth over an area of 21 km² around our experimental site using data from a Terrestrial Laser Scanner. In addition, in-situ measurements were collected during blowing snow events:  Vertical profile of snow fluxes using mechanical traps and Snow Particles Counters (SPC),  Vertical profile of wind velocity, air temperature and relative humidity on a meteorological mast, 

Wind speed and air temperature at three automatic weather stations. These measurements represent a whole dataset suitable for the evaluation of models that simulate wind-induced snow transport in alpine terrain. We present in this paper the first results of this in-situ study.

INTRODUCTION *

While snow transport by the wind has a great influence on the snow depth distribution and on the evolution of the snow pack stability, we must take into consideration blowing snow events in order to improve the operational forecasting of avalanche risk [START_REF] Guyomarc'h | Integration of the snowdrift modeling into the French operational chain for avalanche hazard forecasting[END_REF]. The Snow Study Centre (specialized research centre of Météo-France) and IRSTEA (previously Cemagref) have been working for several years towards upgrading the knowledge on blowing snow and its numerical modeling in alpine areas. After numerous observations and measurements at an experimental site, some numerical models: NEMO [START_REF] Naaim | Numerical simulation of drifting snow: erosion and deposition models[END_REF], Sytron 2 & 3 [START_REF] Durand | Two-dimensional numerical modeling of surface wind velocity and associated snowdrift effects over complex mountainous topography[END_REF][START_REF] Durand | Improvement of a numerical snow drift model and field validation[END_REF] have been developed to simulate the occurrence of such phenomena and the redistribution of the snow pack during those events. More recently, we tested the ability of a meso-scale atmospheric model to explicitly simulate wind-induced snow transport in alpine terrain [START_REF] Vionnet | Simulations of snow-drift events in an alpine terrain[END_REF].

Other recent studies focus on fine scale process during blowing snow events [START_REF] Naaim-Bouvet | Back analysis of drifting snow measurements over an instrumented mountainous site[END_REF] and on the modeling of blowing snow occurrence using the Crocus snow pack model [START_REF] Vionnet | Occurrence of blowing snow events at an alpine site over a 10year period : observations and modeling[END_REF]. The study presented here takes place in this context and takes advantage of a high-altitude experimental site to investigate meteorological conditions during blowing snow events and their consequances on the snow pack.

EXPERIMENTAL SITE

In order to conduct observations on meteorological conditions of blowing snow events associated with snow features and wind velocity, we have used a high altitude site in the French Alps for around 20 years. This experimental site (Col du Lac Blanc) is situated at 2,700 m high close to the Alpe d'Huez ski resort. This pass could be considered as a "natural wind tunnel" due to its surrounding topography (figure 1). Since the beginning of the 1990's, blowing snow conditions were investigated by taking into account wind velocity thresholds according to snow grain characteristics at the snow surface. Three automatic weather stations (AWS) are located around the pass (figure 2) and record meteorological data every 15 minutes: wind velocity (mean, max and mini), wind direction, air temperature and snow depth by using an ultrasonic sensor. Two heating precipitation gauges provide information regarding the amount of precipitation (snow water equivalent). However, as usual at this altitude, they suffer of high uncertainty in the measurements when strong winds occur and that is the case at the most part of the time. The site was also fitted with an ultrasonic anemometer and six cup anemometers mounted on a 10-m vertical mast.

In parallel, blowing snow fluxes are also measured with Snow Particle Counters [START_REF] Sato | Field test of a new snow-particle counter (spc) system[END_REF]. For the winter season of 2011 three SPC were deployed. 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 could be raised manually when the height of snow increased. A snow depth sensor that measures the exact position of the SPC above the snow pack supplemented these devices.

DATA COLLECTED DURING MEASUREMENTS CAMPAIGNS.

For the purpose of this study, we organised intensive measurement campaigns over the last two winter seasons. Based on the forecasting of a blowing snow event (snow storm or blowing snow without concurrent precipitation) we defined a measurement protocol as following:

 before the event: we plan to perform a specific survey of snow height thanks to a terrestrial laser scanning (TLS) (figure 3) of the surrounding topography of the experimental site [START_REF] Prokop | Assessing the applicability of terrestrial laser scanning for spatial snow depth measurements[END_REF] and some measurements of snow characteristics at the snow pack surface.  During the event: measurements of the blowing snow flux are performed manually using snow nets along a vertical profile (figure 4) and at the same levels we record measurements at 1 min time step of a vertical profile of wind velocity and its direction, as well as automatic measurements of a three SPC profile (figure 5).  After the event, we perform again a terrestrial laser scanning of the surrounding topography and make measurements of snow features in accumulations.

Over the last two winter seasons, 6 blowing events have been investigated using the measurement protocol described above (figure 6).

RESULTS.

Measurements of snow fluxes

Measuring snow flux is always a key point when we have to validate numerical models with in-situ data. This explains why we focused our observation on these measurements. In winter season 2010-2011, we have compared profiles of blowing snow fluxes measured simultaneously on two poles separated from around 50 meters. Figure 9 shows the comparison for a blowing snow event during precipitations (16 profiles) and figure 10 without precipitation (4 profiles). The curve in red shows the mean values. The results show no bias between the two locations of measurements but a great variability when snow precipitations occur. We also noticed large differences near the snow surface when blowing snow event occurs without precipitation. The snow flux in saltation is greatly influenced by the micro-topography.

Comparison with SPC data

It appears interesting to compare now the data from the SPC and the measurements from "snow nets". The figures below show a comparison between the fluxes measured with the three SPC and the snow net fluxes for a given period (a) and a comparison with the mean values of snow fluxes that we have interpolated at the height of the intermediate SPC for three blowing snow events. We can observe that the SPC always measure larger values of flux that our manual measurements. The mean value of the flux collected with the nets is around 35% of the SPC flux values. This is probably due to the collection efficiency of the snow nets which measure the snow flux depending of the material porosity. The comparisons shown here are calculated by using 41 measurements and show a similar behavior for the three blowing snow event. Assuming an efficiency of 35% for the snow nets we can get the adjusted profile shown on figure 11(a).

Measurements of vertical profile of wind

For each event in winter season 2010-2011, the wind velocity and direction are recorded with a time step of 1 min along a vertical profile near the snow surface (figure 12). The friction velocity is calculated for each measurement period by using a logarithm fit:

Comparison wind velocity/snow mass

The figure 13 shows an example of the blowing snow density versus the friction velocity. The wind velocity is interpolated at the height of the corresponding net.

Measurements of snow distribution

The applicability of terrestrial laser scanning to measure the depth of the snow cover has been shown by [START_REF] Prokop | Assessing the applicability of terrestrial laser scanning for spatial snow depth measurements[END_REF]. According to meteorological and snow pack conditions, and to the distance to the target, the measurements may provide accurate results (<10 cm at a distance less than 500 m).

For this study, we have used the difference between two images calculated from the TLS which deliver a snow map around our experimental site (figure 14). These maps will be very useful to verify the results of blowing snow models.

CONCLUSION

Meteorological recordings and in-situ measurements were combined to describe as fine as possible the conditions of blowing snow events and their consequences on the snow pack distribution and stability. Six blowing snow events were investigated during these observation campaigns. include blowing snow fluxes from snow nets and SPC, near-surface meteorological conditions and maps of snow depth difference using TLS. A detailed analysis of the collected profiles will be conducted.

All of the data collected represent a complete set that we can use for the assesment of models of snow transport by the wind in Alpine terrain and notably for new works on Meso-NH/Crocus modeling [START_REF] Vionnet | Simulations of snow-drift events in an alpine terrain[END_REF].

Lots of works remain to be done to better understand the mechanisms involved in blowing snow events at a fine scale.

Figure 1 :

 1 Figure 1: View to South-East of the experimental site.

Figure 2 :

 2 Figure 2: View to North-East of the experimental site with the location of the 3 Automatic Weather Stations.

Figure 3 :

 3 Figure 3: The Terrestrial Laser Scan at the northerly point of measurement.

Figure 5 :

 5 Figure 5: Vertical profile of wind velocity (on the left) and vertical profile of SPC (on the right).

Figure 4 :

 4 Figure 4: Vertical profile of snow nets used for the measurement of snow concentration in the atmosphere ng a blowing snow event. duri

Figure 6 :

 6 Figure 6: These 2 graphs show the periods of measurement campaigns over 2010-2011 and 2011-2012 winter seasons.

Figures 7

 7 Figures 7 and 8 represent examples of vertical profiles of snow flux during blowing snow events respectively with and without concurrent snow falls. The three curves show different wind velocities during measurements and thus an increasing quantity of collected snow.

Figure 7 :

 7 Figure 7: three examples of snow flux profiles with precipitation.

Figure 8 :

 8 Figure 8: three examples of snow flux profiles without precipitation.

Figure 9 :

 9 Figure 9: Comparison of snow flux measured at two positions separated from 50 m (A & B) during a blowing snow event with precipitations.

Figure 10 :

 10 Figure 10: Same as figure 9 for a blowing snow event without precipitation

  Figure 11 : comparison between the fluxes measured with the SPC and the snow net fluxes for a given period and for 3 blowing snow periods.

  depending on the wind direction as shown by: Naaim-Bouvet et al, 2010.

Figure 12 :

 12 Figure 12: Average velocity of wind at 3.5 m versus the calculated friction velocity.

Figure 13 :

 13 Figure 13: measurements performed on 17 th and 18 th of march 2011. The lower net is 20 cm above the snow surface.

Figure 14 :

 14 Figure 14: example of difference map of snow depth between 12 to 23 February 2011.

Proceedings, 2012 International Snow Science Workshop, Anchorage, Alaska


C Mugnai

M Gerino

M Frignani

S Sauvage

L G Bellucci

Bioturbation experiments in the Venice Lagoon

Keywords: bioturbation, luminophores, sediment mixing, modeling, Venice Lagoon

Short experiments (14-21 days) were carried out during autumn 1998 and spring 1999 at one selected site of the Venice Lagoon to measure bioturbation activities and mixing rates, as well as to obtain quantitative information on benthos functionality. Fluorescent sediment particles (luminophores, 63-350 µm) were introduced as pulse inputs at the sediment surface. The concentration-depth profiles of the tracer were simulated with a new advectiondiffusion-non local model applied under non-steady state conditions. This allowed the quantification of the mixing parameters associated with different mechanisms: biodiffusion (D b ), bioadvection (W) and non-local mixing (Ke, z 1 , z 2 ). A parameter RS (removed sediment) was also calculated to account for the flux of sediment due to nonlocal transport. Results show that bioturbation was dominated by biodiffusion in autumn and by bioadvection in spring. Mean mixing parameters Db, W, and RS changed from 3.09 to 0.87 cm 2 y -1 , from 0.93 to 15.50 y -1 and from 5.85 to 7.79 g cm -2 y -1 , respectively.

Introduction

Bioturbation is defined as the process of sediment mixing that results from macrofauna burrowing, feeding, and reworking within the surficial sediment. This process has profound effects on both early diagenesis and preservation of sediment records. In fact, the rates of oxygen penetration and organic matter decomposition, the dissolution of biogenic components such as CaCO 3 and SiO 2 , and the pore water concentrations of nearly all dissolved species are affected by bioturbation as a function of its intensity [START_REF] Berner | Early Diagenesis. A Theoretical Approach[END_REF][START_REF] Aller | The effects of macrobenthos on chemical properties of marine sediment and overlying water[END_REF][START_REF] Wheatcroft | A mechanistic view of the particulate biodiffusion coefficient: step lenghts, rest period and transport directions[END_REF][START_REF] Wheatcroft | Experimental tests for particle sizedependent bioturbation in the deep ocean[END_REF]. Furthermore, the displacement of interstitial fluids, the increase of sediment surface in contact with bottom water, and irrigation enhance the exchange of dissolved species at the sediment-water interface. In addition, the displacement of particles and particle reactive species blurs or destroys the sediment record, i.e., the information that sediment can store in their strata. Therefore, a better understanding of bioturbation in its mechanisms, rates, and effects is important for a wide range of disciplines.

In the past, biogeochemists have approached bioturbation through simplified models based on the analogy with diffusion (e.g., [START_REF] Wheatcroft | A mechanistic view of the particulate biodiffusion coefficient: step lenghts, rest period and transport directions[END_REF]. The community-wide mixing rates were obtained by measuring the vertical concentration profiles of various tracers and fitting theoretical profiles to the observed gradients. The calculated biodiffusivity (D b ), which integrates all bioturbation activities over a sufficiently long time, cannot explain short time-scale effects. Among these, a tracer distribution with a maximum at the surface, followed by a sharp decrease with depth, is typically generated by biodiffusive mixing [START_REF] Guinasso | Quantitative estimates of biological mixing rates in abyssal sediments[END_REF][START_REF] Cochran | Particle mixing rates in sediments of the eastern equatorial Pacific: evidence from 210 Pb, 239,240 Pu and 137 Cs distributions at MANOP sites[END_REF][START_REF] Wheatcroft | A mechanistic view of the particulate biodiffusion coefficient: step lenghts, rest period and transport directions[END_REF]. This mixing is produced by organisms that move sediment particles in a random manner over short distances. On the contrary, subsurficial peaks are generated by vertical transfer of the tracers to depth caused by two independent processes: bioadvection [START_REF] Fisher | Vertical mixing of lake sediments by tubificid oligochaetes[END_REF] and non local transport [START_REF] Boudreau | Mathematics of tracer mixing in sediments: II Non local mixing and biological conveyor-belt phenomena[END_REF]. Bioadvective transport is generated by head-down oriented organisms called 'conveyor-belt species' that cause an active transport of sediment trough their gut from the ingestion zone located at depth in the sediment to the sediment surface [START_REF] Fisher | Vertical mixing of lake sediments by tubificid oligochaetes[END_REF][START_REF] Rice | Early diagenesis in bioadvective sediments: relationships between the diagenesis of 7 beryllium, sediment reworking rates, and the abundance of conveyor-belt depositfeeders[END_REF][START_REF] Gerino | Quantitative estimation of biodiffusive and bioadvective sediment mixing: in situ experimental approach[END_REF]. Their effect on a tracer profile is similar to that due to a high accumulation rate. In turn, non-local mixing is caused by two mechanisms: one is due to the egestion of fecal pellet at depth by surface deposit feeders denominated 'inverse conveyors', whereas the other, called 'regeneration', supplies fresh material that mechanically drops into burrows from the interface to depth (Smith et al, 1986/87;[START_REF] Gardner | A regeneration model for the effect of bioturbation by fiddler crabs on 210 Pb profiles in salt marsh sediments[END_REF][START_REF] François | A new approach for the modelling of sediment reworking induced by a macrobenthic community[END_REF].

Several models have been proposed to simulate bioturbation mechanisms different from diffusion [START_REF] Fisher | Vertical mixing of lake sediments by tubificid oligochaetes[END_REF]Smith et al., 1986/87;[START_REF] Gardner | A regeneration model for the effect of bioturbation by fiddler crabs on 210 Pb profiles in salt marsh sediments[END_REF][START_REF] Boudreau | Diagenetic models and their implementation[END_REF][START_REF] Gerino | Quantitative estimation of biodiffusive and bioadvective sediment mixing: in situ experimental approach[END_REF][START_REF] Gerino | Comparison of different tracers and methods used to quantify bioturbation during a bloom: 234 Thorium, luminophores and chlorophyll a[END_REF]. In particular, [START_REF] Gerino | Comparison of different tracers and methods used to quantify bioturbation during a bloom: 234 Thorium, luminophores and chlorophyll a[END_REF] used a particle transport model containing both diffusive and advective terms to simulate the depth distribution of the pulse input of a conservative tracer. An evolution of this model, which is able to account for non local mixing, is used in this paper. In the present study, an experiment was carried out to obtain information on bioturbation mechanisms and rates in a shallow environment of the Venice Lagoon. This site was selected because it is not disturbed by marine traffic and by the illegal clam fishing that affects a great part of the lagoon.

Materials and methods

Figure 1 shows the study area and the experimental site. G is located at 0.5 m depth in the central part of the Palude di Cona. It corresponds to a typical shallow water environment with estuarine characteristics that receives the discharge of the Dese river, one of the main freshwater suppliers to the lagoon [START_REF] Zonta | Heavy metal and grainsize distribution in estuarine shallow water sediments of the Cona Marsh (Venice Lagoon, Italy)[END_REF].

The experimental procedure, based on the use of dyed sediment particles called luminophores, was described in detail by [START_REF] Gerino | The effect of bioturbation on particle redistribution in Mediterranean coastal sediment. Preliminary results[END_REF] and [START_REF] Gerino | Quantitative estimation of biodiffusive and bioadvective sediment mixing: in situ experimental approach[END_REF]. Experiments were carried out in autumn 1998 and in spring 1999. At the experimental site, three open PVC tubes (25 cm long, 12.2 cm i.d.) were inserted into the sediment and a thin frozen mud cake (0.5 cm thick), containing 2 g of both pink and green tracer (63-125 µm), was set at the sediment surface thus simulating a pulse input. One control core with defaunated sediment was set in similar conditions. The tubes were carefully retrieved after 15 days in autumn and 21 days in spring, quickly transferred to the laboratory and subsampled at 0.5-5 cm intervals. Luminophores were counted with a microscope under UV light, and counts were converted into weight units using the number of particles per mass of lumino-phores in a given size fraction divided by the weight of sediment. Furthermore, the tracer profile in each core was normalized against its inventory.

Results and discussion

The tracer vertical distributions (Fig. 2) show evident penetration of luminophores in the experimental cores. On the contrary, the tracer remained confined at the surface in the control cores, even if in spring there is an expansion of the surficial peak down to a depth of 2 cm, thus excluding physical disturbance. However, the three depth-distributions obtained from each experiment are never identical, in that they exhibit multiple tracer peaks at different depth in the subsurface zone. In particular, the autumn profiles G(a) show an average penetration of 8 ± 3 cm (median ± range of variation, n = 3), with the highest concentration close to the surface. G(a) 1 shows three well-defined subsurface peaks (Fig. 2), whereas tracer accumulations are less marked and closer to the surface in G(a) 2 and G(a) 3 . In spring cores, especially G(s) 1 and G(s) 3 , the maximum concentrations are located immediately below the sediment surface. These distributions account for the presence of advective transport. The heterogeneity of biological transport processes in replicate cores is probably an effect of the non-homogeneous horizontal macrofauna distribution at the decimetric scale of the cores.

To estimate bioturbation parameters, the profiles were simulated using the bioadvection-biodiffusionnon local model in non-steady-state conditions. This model was developed by adding a non-local component to the regular diffusion advection model [START_REF] Officer | Interpretation procedures for the determination of sediment parameters from time-dependent flux inputs[END_REF][START_REF] Gerino | Quantitative estimation of biodiffusive and bioadvective sediment mixing: in situ experimental approach[END_REF][START_REF] Gerino | Comparison of different tracers and methods used to quantify bioturbation during a bloom: 234 Thorium, luminophores and chlorophyll a[END_REF]. The nonlocal mixing is modeled such as a removal function that determines the fraction of tracer removed from the surface within the particle collection zone, and an injection function that simulates the displacement of this material into the deposition zone. Since the experiments were conducted over a very short time period, and advection was not recorded in the control cores, the sedimentation rate is not taken into account in the model. The basic equation is:

∂C(z, t) ∂t = D b ∂ 2 C(z, t) ∂z 2 -W ∂C(z, t) ∂z (1) + K(z, t) -R(z, t),
where C is the normalized tracer concentration, t is the time (years), z is the depth (cm) positive down- ward, D b is the biodiffusive mixing rate (cm 2 y -1 ), W is the bioadvective transport rate (cm y -1 ), and R is the removal function that determines the fraction of tracer (g cm -3 y -1 ) removed from the surface. This quantity of tracer is collected at the sediment surface by inverse conveyor or lost from the sediment surface by passively falling into the burrow. K is the injection function of the non-local transport that simulates tracer inputs (g cm -3 y -1 ) into the injection zone of the sediment column and Ke is a constant parameter (y -1 ) estimated from the model. The depths z 1 and z 2 represent the upper and lower limits of the injection zone respectively. The non-local transport is thus quantified by a flux of sediment 'removed from the surface', called Removed Sediment (RS, g cm -2 y -1 ). Knowing that sediment was removed from a marked mud cake with a thickness of 0.5 cm, surface S (cm 2 ), density d (g cm -3 ) and volume V (cm 3 ), RS can be expressed as a function of Ke (y -1 ) that is RS = Ke (V/S) d. A parameter τ =0.5/W max is introduced to account for the time during which the tracer is still at the top 0.5 cm of surface sediment and limits the period available for non-local transport. In this case: and boundary conditions

R(z, t) = 0 for z > 0.5 or t > τ k(z, t) z 2 -z 1 0.5 for z < 0.5 and t < τ K(z, t) = Ke for z ∈ [z 1 , z 2 ] 0 for z / ∈ [z 1 , z 2 ] with model initial conditions C(z, t = 0) = 1 for z ∈ [0, 0.5] C(z, t = 0) = 0 for z > 0.5
C(z → +∞, t) = 0 W C (z=0,t ) -D b ∂C ∂z ] z=0,t = 0
an analytical solution is provided by Officier & Lynch (1982):

C(z, t) = 1 √ D b t exp - (z -W t) 2 4D b t - W 2D b exp( W z D b erf c z + W t √ 4D b t -Ret + Ket with Re = Ke Z 2 -Z 1 0.5 (2)
Since experimental time remains relatively short compared to advection velocities, and because the highest tracer concentrations remain very close to the surface, it is assumed that the initial tracer impulse has never migrated down to the limit depth where the organisms ingest food and sediment (Fischer et al., 1980) via bioadvective processes. The model allows the calculation of the theoretical tracer concentration given suitable values of the parameters D b , W , z 1 , z 2 , Ke, τ . These parameters are selected from profiles that produce the best fit with the experimental data using the least square method. When two or three subsurficial peaks were present, the upper tracer profile was assigned to bioadvective and/or biodiffusive processes, and other deeper tracers accumulations were taken into account as due to non local processes. The total theoretical profile is the sum of the theoretical sub-profile concentrations. Figure 3 shows an example through the bestfitting of the tracer depth-distribution in core G(a) 1 .

Table 1 shows the values of D b , W , Ke, z 1 , z 2 , and RS estimated from each experimental profile. Coefficients are comparable with those reported in the literature [START_REF] Aller | 234 Th/ 238 U Disequilibrium in near-shore sediment: particle reworking and diagenetic time scales[END_REF][START_REF] Benninger | Effects of biological sediment mixing on the 210 Pb chronology and trace metal distribution in a Long Island Sound sediment core[END_REF][START_REF] Clifton | The use of radionuclides in describing the mixing characteristics of estuarine sediments[END_REF][START_REF] Gerino | Quantitative estimation of biodiffusive and bioadvective sediment mixing: in situ experimental approach[END_REF][START_REF] Gerino | Comparison of different tracers and methods used to quantify bioturbation during a bloom: 234 Thorium, luminophores and chlorophyll a[END_REF] and in most cases obtained with the same method in environments with a morphology similar to that of the Venice Lagoon. D. Prevedelli (pers. com.) studied the macrofaunal density and composition at the experimental site in the two seasons. Results suggest that oligochaeta and the polychaeta Hediste diversicolor and Streblospio shrubsolii are the main responsible of non-local mixing. On the other hand, the seasonal differences of W could be attributed to several factors. In particular, the increased temperature could have enhanced biological activities or induced changes in the functional composition of the benthic community.

Conclusions

The following conclusions can be drawn: (1) the luminophores penetrated into the sediments only due to the activity of benthic organisms in that physical processes were not effective, as demonstrated by control cores; (2) the short time scale of the experiments and the use of the model allowed a clear discrimination and quantification of the diffusive, advective, and nonlocal processes under realistic in situ conditions. In spite of that, 2-3 weeks were not sufficient for the tracer to reach the depth of maximum ingestion, since the advective subsurficial peaks remain well defined;

(3) the polychaeta H. diversicolor and S. shrubsolii are probably the major responsible of the non local mixing associated with deep tracer peaks; and (4) repeated experiments evidence a greater importance of bioadvective processes during spring than during autumn. This change might be attributed to several factors like temperature, macroinvertebrate community composition, or organic matter content. Identification of environmental factors controlling bioadvection should be further considered in situ or in laboratory experiments. In a future prospective, the inclusion of sedimentation in the bioturbation model will let us to simulate the formation of the sediment record. This will enable to make provision for the effects of changes of the contaminant inputs, very important in deciding any monitoring program.

Table 1. Bioturbation coefficients and RS values from model simulations. D b ; biodiffusive mixing rate; W; bioadvective mixing rate; Ke; downward transport of tracer displaced by 'non-local mixing'; z 1 , z 2 ; upper and lower limits of the 'non-local mixing' zone; RS; removed sediment from the surface layer by non local transport, expressed per unit area of sediment-water interface. G(a) and G(s) represent average coefficients (mean ± confidence limit, n = 3, p < 0.05) for each season (a = autumn and s = spring). G(a) c and G(s) c are control core coefficients. These last values are given for information and do not enter in the calculation of average coefficients.

Site

Figure 1 .

 1 Figure 1. Study area and experimental site.

Figure 2 .

 2 Figure 2. Concentration-depth profiles of luminophores in sediments. G(a) 1 , G(a) 2 , G(a) 3 and G(s) 1 , G(s) 2 , G(s) 3 are the experimental replicates obtained in autumn and in spring, respectively; G(a) c and G(s) c are the defaunated control cores. Values are normalized against total inventories.

Figure 3 .

 3 Figure 3. Best-fitting of the experimental profile G(a) 1 (core with fauna at site G during autumn 1998) with a theoretical profile. The best-fitting of the corresponding control (core without fauna at site G during autumn 1998) is also shown.

Acknowledgements

Funds for this research were provided by the 'Project 2023' founded by the Ministry of Public Works through the Magistrato alle Acque of Venice and the Consorzio Venezia Nuova. This is contribution No. 1296 from the Istituto di Geologia Marina, Bologna (Italy).


Olivier Crespo
email: olivier.crespo@csag.uct.ac.za

J.-E É Bergez
email: jbergez@toulouse.inra.fr

Frederick Garcia
email: fgarcia@toulouse.inra.fr

P2q hierarchical decomposition algorithm for quantile optimization : application to irrigation strategies

Keywords: Simulation-based optimization, quantile optimization, irrigation management 1 Introduction Quantiles

risk management in finance, insurance, or banking [START_REF] Jorion | Value at risk[END_REF]. They are also largely involved in quantile regression [START_REF] Koenker | Quantile regression[END_REF] for optimal allocation design. [START_REF] Rostek | Quantile Maximization in Decision Theory[END_REF] lists advantages of using the quantile instead of the mean or standard deviation when faced with a decision problem under uncertainty. The main ones are the improved robustness to fat tails and the predictions not driven by outliers. The improved robustness is due to the uselessness of responses quantification. The user does not need indeed to prefer one solution 10 times more than another, he just needs to rank them. This implies that no a priori preferences are needed, which allows to resolve the problem without user dependence. This conclusion is also true for maxmax and maxmin optimization. The former consists in maximizing the maximum payoff of a random variable, since the latter consists in maximizing the minimum one. The first kind of optimization is said optimistic, and the decision maker target is to increase the best expected outcome without regard to the worst cases. The second kind of optimization is said pessimistic, and the decision maker try to limit the worst achievable outcome by maximizing it. These criteria have been applied to game theory, robust control, bargaining and others areas. However these criteria have been criticized for basing the choice on what may be extreme and unlikely outcomes. Quantile optimization is then intermediate between the extremes maximum and minimum maximization. Robustness is still improved compared to mean and standard deviation analysis, and predictions are not driven by unlikely outliers. Even more, no user preferences are needed but a risk acceptance, called anticipation level in [START_REF] Rostek | Reasoning in Strategic and Non-Strategic Interactions[END_REF].

The aim of this study is to design a simulation optimization algorithm adapted to quantiles optimization, and to analyze on some agricultural applications, the worth of optimizing quantiles instead of expected values. Due to the complex nature of our application, we assume that computing probability distributions or utility function is too expensive. We thus focus on optimization methods that do not need prior knowledge of the uncontrollable parameters. Simulation optimization concerns systems without an explicit analytical expression, but where instead only responses of sampled stochastic simulation runs L are available. The formulation of that kind of optimization involves input variables θ, constraints on these decision variables that should take value into a feasible region Θ and an objective function defined on these variables, J : θ → R. The common formulation of the optimization problem consists in defining J(θ) as the expected value E[L(θ, ω)] with ω the uncontrollable input variables vector of the stochastic system [START_REF] Andradóttir | A review of simulation optimization techniques[END_REF], [START_REF] Ólafsson | Simulation optimization[END_REF]Kim (2002) Fu et al (2005)). Quantile optimization consists then in replacing the expected value of L by the quantile of order α:

max θ∈Θ J(θ) = qα[L(θ, ω)]. (1)
We present in this paper P 2q, a new and original algorithm for solving this difficult simulation optimization problem. We applied the P 2q algorithm to the problem of optimizing irrigation strategies for corn, based on the crop simulation model MOD-ERATO [START_REF] Bergez | Moderato: an object-oriented decision tool for designing maize irrigation schedules[END_REF] in order to get the greatest simulated yield. This case application is a good illustration of the kind of simulation optimization problems that exists in crop management.

The plan of the paper is the following: first, we define quantiles and give limitation to the quantile optimization problem. Next, we describe the principle of the P 2q algorithm we propose and the risk interpretation of the quantile of order α. In section 4, we give some information on the biodecisional simulation model used to simulate the irrigation strategies and we set up the simulation experiments. We then present the

q α 1-α J ] E[ α J P(J)
Fig. 1 Quantile of order α separating the probability distribution of a random variable.

results of the quantile optimization for different quantile orders, and we discuss the pros and cons of this quantile optimization approach.

Quantile optimization

As extremes pessimistic and optimistic indices can be optimized with maxmin and maxmax approaches, a non extreme risk acceptance index can be improved using quantile optimization. Furthermore in most of the cases, quantile estimation will rely on more possible outcomes than the unique extreme taken into account for either maxmax or maxmin optimizations. In opposition to the usual mean expected value provided to decision makers, the quantile estimation can provide relative pessimistic or optimistic options in accordance to decision risk acceptance. This new optimization point of view is a tool that can provide relevant advices with both farmers that want to increase the X% worst possible outcomes or those who focus on the X% best.

Let α be a fixed probability in ]0, 1[. The quantile of order α of the random variable J is the threshold qα such that the probability of observing J lower than qα is equal to α: P (J < qα) = α. For instance, if we are interested in the 10% highest values of J, we will compute the quantile of order α = 0.9. If we take α = 0.5 then we get the median. Figure 1 shows the probability distribution of the random variable J. The quantile qα of order α is dividing the surface into two parts. The left hand side represents the probability of reaching a value lower than qα: it is equal to α. The right hand side represents the probability of reaching a value greater than qα: it is equal to 1 -α.

The optimization problem we consider consists in optimizing the α-quantile of the random variable L(θ, ω). The decision maker has thus to characterize his risk acceptance through the order α. This scalar α ∈ [0, 1] is a subjective belief up to the decision maker, such that α = 0 is equivalent to a maxmin (pessimistic) approach and α = 1 is a maxmax (optimistic) approach.

The common objective function is transformed from the expected value of simulation outcomes into the quantile of order α, which leads usually to different decisions. Thus the use of quantile in simulation-based optimization is not vain and we propose hereafter to evaluate its interest for irrigation strategies design.

3 The P 2q algorithm

A hierarchical algorithm

The P 2q algorithm is an extension of the P 2 and P 2p algorithms we introduced in [START_REF] Crespo | P2 hierarchical decomposition procedure: application to irrigation strategies design[END_REF] and [START_REF] Bergez | A hierarchical partitioning method for optimizing irrigation strategies[END_REF]. It is based on a hierarchical decomposition of the decision space into a binary tree and belongs to the family of stochastic branching methods, like stochastic branch-and-bound or nested partitions methods. It is dedicated to large continuous input variables problems. It encloses a single objective function, which can be the combination of weighted multiple objectives.

The decision space Θ is a hyper-rectangle or region. The P 2q optimization aims at finding small hyper-rectangles that we call R which include the decision vectors that maximize the quantile of order α of performance measures L(θ, ω). We assume that a hyper-rectangle is small enough, or unbreakable, when any decision vector included in this hyper-rectangle is indiscernible from the others. This is defined by the user for every dimension d of the D-dimensional problem as the width step of the dimension d ∈ D.

Let us call pending regions the regions that are still breakable (in opposition to the unbreakable ones). Initialization allocates the initial decision space as the single pending region of the list pendingRL of all the pending regions. The first step consists in selecting out from the pending regions list pendingRL, the region which is potentially optimal: we call it the promising region. This promising region is either selected in accordance to one or two indices evaluated from the decision vectors included in the pending regions, or it is randomly chosen following a low probability set for the whole process. The second step divides this promising region into two parts either following a trivial equal partitioning process or techniques looking at producing the most different sub regions. These two parts are offspring regions and are collected in the of f springRL list. During the third step, each of these of f springRL regions is sampled uniformly, simulated and quantiles are evaluated. This evaluation is the main concern of this study. Eventually, the pending regions list pendingRL is updated, and the three previous steps are repeated until stopping criteria are completed (usually time and/or simulation runs limits) or the pending regions list is empty.

The three main steps of selection, division and evaluation are discussed in [START_REF] Crespo | P2 hierarchical decomposition procedure: application to irrigation strategies design[END_REF]. The combination included in the P 2q algorithm involves the three following techniques. Selection considers only the quantile of order α of the performance measures simulated for each pending region R ∈ pendingRL. The division technique involves the determination of one included vector which translates as best as possible an ideal (defined by the user) separator vector called pivot. The promising region is then cut apart along the previously defined pivot vector in order to split the longest dimension d ∈ D into two parts. The quantile of performance measures is computed during the evaluation step as defined in the following section 3.2.

Quantile evaluation

Let us call R ⊂ Θ a decision subspace that includes decision vectors θ. The evaluation of R relies on N ×M performance measures L(θ i , ω j ) simulated

(i = [1..N ], j = [1..M ]).
The number N of θ simulated is proportionally related to the widest dimension of the considered subspace. These N decision vectors θ are uniformly sampled within the R decision subspace. Each of these decision is simulated coupled with M uncontrollable parameter ω randomly selected within a collection of either historical observed or computationally generated ω. The amount M of uncontrollable ω is set for the whole process.

Let α ∈ [0, 1] be a scalar such that the quantile q 0 [L(θ i , ω j )] is the smallest value of the performance measures and q 1 [L(θ i , ω j )] the largest. We use a weighted average method to estimate the quantile of order α. The first step consists in ordering the performance measures within a region such that L 0 is the smallest value and L (N ×M)-1 is the largest. Let p be the integer part of ((N × M ) -1) × α, and let g be its fractional part. The weighted average estimation of the α order quantile of performance measures is done as follow.

qα[L(θ i , ω j )] = Lp + g × (L p+1 -Lp) (2)
According to the previous quantile estimation process, the evaluation criterion J(R) of the decision subspace R is computed as follow.

J(R) θ∈R = qα[L(θ, ω)].
(3)

Application to irrigation management

We applied this new approach to design of irrigation strategies. The aim of this section is to present the problem we used to test the P 2q algorithm. [START_REF] Wallach | Parameter estimation for crop models a new approach and application to a corn model[END_REF]. In short, the model is based on the general captation / conversion of the solar radiation allowing for leaf area and biomass development.

Depending on phenological stages, efficiency of radiation conversion differs. Yield is calculated thank to the use of a dynamic harvest index. A water budget integrating transpiration and evaporation is computed on a 4-layer soil. Water stress modifies plant growth and harvest index.

The decision model consists of a set of decision rules for different management decisions. A decision rule can be seen as a function linking some indicators of system status, in this case the cropping system, to the action to be performed. It is written as a Boolean condition: "IF <indicator> <operator> <threshold> THEN <action1> ELSE <action2>". A crop management strategy is composed by three simple rules (Sowing, Fertilizing and Harvest) and by the Irrigation Strategy, which, in its turn, is defined by five elementary rules: i) a rule to determine if irrigation is to be used to facilitate plant emergence (Irrigation at Sowing Rule); ii) a rule to decide when to start the main irrigation period (Starting Rule); iii) a rule to determine when to start a new irrigation round (Returning Rule); iv) a rule to delay irrigation due to weather conditions (Delaying Rule); v) a rule to decide when to stop irrigation (Ending Rule). The crop model updates the state variables each day and passes their values to the decision model together with the explanatory variables of that day.

The constraints to irrigation and available in MODERATO may be divided into three categories:

Equipment constraints: maximum flow available due to pumping capacity per hectare; maximum and minimum amount of irrigation allowed. Resource constraints: a limited total amount of water for irrigation and a decreased flow rate during the irrigation period. Regulatory or human constraints: irrigation bans on some days of the week, farmers reluctance to work for all or part of certain days, etc.

Constraints may vary during the irrigation period. For example river flow may decrease in summer, reducing the pumping flow. Regulations may be issued to reduce water consumption because of water shortage, possibly by prohibiting irrigation on certain days. These modifications of the general context over the irrigation period can be taken into account in the irrigation scheduling. However, as MODERATO is used at a strategic level, those modifications have to be anticipated to be tested.

An irrigation strategy is then defined as the set of rules (indicators -thresholds) activated for a given irrigation context.

MODERATO can run over a climatic series but every year is independent.

Case study

Tests have been run on an eight-parameter strategy as follows.

The main irrigation period starts from T1 as soon as the soil water deficit reaches D1. An amount I1 is applied. Once an irrigation cycle ends, a new cycle starts when the soil water deficit reaches D2. An amount I2 is applied.

For the irrigation cycle following T3, if the soil water deficit is greater than D3 before this irrigation cycle starts, a last irrigation cycle is performed; otherwise the irrigation campaign ends. An amount I3 is applied.

Other cultural operations are modeled with decision rules, such as the sowing or the harvest. Details of these rules are provided in the table 1. Soil properties and climatic series used are describing south west of France. The variation in rainfall during the two summer months as it ranges from 30 to 240 mm, underlining the unpredictable nature of rainfall in the area.

The irrigation equipment used for the study allows a 3.5 mm/day flow rate. A 180 mm limitation of available water is applied. No flow rate restrictions during summer (except those due to the equipment) are imposed.

The objective function to be maximized is the α-quantile of the direct margin (i.e. the gross margin minus specific costs for a given activity, here irrigation). The direct margin regarding irrigation can be written as a weighted sum of multiple criteria: where L(θ i , ω j ) is the direct margin for climate ω j and the strategy θ i , a(θ i , ω j ) is the grain yield obtained under climate ω j and using the strategy θ i , B is the selling price for corn, C is the operational costs for corn production, d(θ i , ω j ) is the amount of water used under climate ω j and using the strategy θ i , E is the cost of irrigation water, f (θ i , ω j ) is the number of irrigation cycles performed under climate ω j and using the strategy θ i and G is the cost of carrying out a new irrigation cycle. The average selling price for maize (grain) is assumed to be 106.71 ¤/Mg in the Toulouse area. Operational costs (seed, weeding, fertilizer, insurance) are assumed to be 327.77 ¤/ha. The cost of irrigation water is assumed to be 0.76 ¤/mm and the setting up of a new irrigation cycle is assumed to be 7.62 ¤/ha.

L(θ i , ω j ) = a(θ i , ω j ).B -[C + d(θ i , ω j ).E + f (θ i , ω j ).G] (4)

Experiments

We considered 10 replications for each alternative. The initial feasible region is defined in table 2 as the ranges of the different parameters from the strategy described in section 4.2. Some procedure parameters were set for all experiments. The maximum number of simulation was set to 2 million. The probability of randomly choosing the promising region is set to 20% and the sampling follows a uniform distribution. The density of samples which are use to evaluate a decision subspace is set to 1. According to the P 2 procedure, it means that the region R evaluation is based on at least as many θ as step number along the widest dimension. Finally the number M of uncontrollable ω is set to 49 and they are randomly chosen in weather data observed in the area of concern.

Each time the direct margin quantile increased, we stored it as well as the number of simulation runs used to achieve it. These stored values presented a strictly increasing Table 2 The eight parameters of the irrigation strategy to be optimized. Min and max show the range of each parameter within which the optimum is sought. A step is the minimum feasible range of the parameter.

Names

Meaning unit min max step T1 Accumulated thermal unit to start the irrigation campaign curve in the simulation runs/direct margin criteria space. The general shape of this curve was made up of two distinct phases: one with a large improvement in direct margin quantile with a few simulation runs, and another with a tiny improvement with an infinite number of simulation runs. The large amount of simulation runs insured that every optimization processes reached the second phase of low improvement.

In order to highlight the advantages of quantile optimization, we first compare optimal regions reached for expected value and quantile optimizations. Then we compared quantile optimization of order 0.1 to order 0.9. The optimization process included in P 2q is a minimization. The translation of any maximization into a minimization allows any problem to be handled. However in this case, in opposition to what we said in the earlier quantile definition section, the quantile of order 0.1 divided the 10% best solutions from the 90% worst, and the quantile of order 0.9 divided the 90% best from the 10% worst.

Results

Quantile vs. expected value optimizations

Figure 2 shows 50 random strategies simulated in four regions. The + are strategies included in the initial research domain. The × are strategies included in the optimal region reached with expected value optimization. The * are strategies included in the optimal region reached with 0.1-quantile optimization. The are strategies included in the optimal region reached with 0.9-quantile optimization. Every strategy is plotted two times. The first is upside and is located on the graph by its average value and its 0.1-quantile. The second is downside and is located on the graph by its average value and its 0.9-quantile.

The optimal regions reached high α-quantile and expected value comparing to the strategies simulated within all the decision space. In one hand, the 50 strategies simulated in the 0.1-quantile optimal region are maximizing the threshold dividing the 10% best responses and the 90% lower ones. It reflects an optimistic decision indicator focusing on optimizing the best responses. In the other hand, the 50 strategies simulated in the 0.9-quantile optimal region are maximizing the threshold dividing the 90% best 100 200 300 400 500 600 700 250 300 350 400 450 500 550 direct margin 0.1 quantile (euros/ha) direct margin 0.9 quantile (euros/ha) direct margin average value (euros/ha) 50 strategies in expected value optimal region 50 strategies in initial domain 50 strategies in 0.1 quantile optimal region 50 strategies in 0.9 quantile optimal region direct margin (euros/ha) quantile of order 0.9 order 0.9 quantile optimization order 0.1 quantile optimization quantile of order 0.1 responses and the 10% lower ones. It reflects a pessimistic decision indicator focusing on optimizing the worst responses.

Different quantile orders optimizations

Figure 3 shows the cumulative distributions of the optimal regions respectively reached with 0.1-and 0.9-quantile optimization. These distributions rely on 2450 performance measures, so that the 0.1-quantile divided the 245 best measures from the 2205 others and so that the 0.9-quantile divided the 245 worst measures from the 2205 others. These thresholds are pictured as ⊓ upper bounded by 245 and 2205 measures from a distribution quantile to the other. So the width of these ⊓ translates the improvement of the α-quantile. For example, 0.9-quantile reached a value larger than 400 ¤/ha when the optimized quantile order is α = 0.9, and it reached a value close to 360 ¤/ha when the optimized quantile order is α = 0.1.

Latter figures characterized criteria responses of optimal regions. The next figures 4 and 5 are multiple axis charts representing the decision variable space. Couples of thin lines bound the parameters defining the optimal regions that were reached for the 10 replications of the P 2q algorithm. The couple of bold lines is an envelope of all of them. Each axis, leveled from 0 to 100, depicts each parameter from the minimum (0%) to the maximum (100%) values available and defined in table 2.

The largest 0.1-quantile direct margin reached over all replications is 691 ¤/ha and an average value of 537 ¤/ha (figures 3 and 4). The related decision vector depicted a first irrigation started after observing small accumulated temperature (below 725 o C.day) and water deficit (below 60 mm), and the application of a medium amount of water (between 18 and 32 mm 8 out of 10 times). The new irrigations are performed as soon as a small water deficit is observed (less than 55 mm) and a large amount of water is applied (more than 32 mm). irrigation campaign end are more variable. Yet we noticed that 6 out of 10 times the accumulated temperature T 3 leading to the irrigation campaign end is included in the 25% of low values (below 1440 o C.day). In the meanwhile 8 out of 10 times the water deficit D3 is included within the 25 and 50% of the proposed values (respectively 52 and 85 mm). The last irrigation amount does not show any of the previous consistencies.

The largest 0.9-quantile direct margin reached over all replications is 412 ¤/ha and an average value of 543 ¤/ha (figures 3 and 5). The related decision vector depicted a first irrigation started after observing a wide accumulated temperature (from 200 to 1000 o C.day) and small water deficit (below 40 mm), and the application of a variable amount of water. The new irrigations are performed as soon as a small water deficit is observed (less than 40 mm) and a large amount of water is applied (more than 32 mm). The accumulated temperature T 3, the water deficit D3 and the water amount I3 leading to the last irrigation are very variable and no consistency is ostensible.

Discussion

Quantile optimization efficiency

As we expected, P 2q as P 2, reached optimal regions with high objective function average value. The highest direct margin average value was achieved with the usual expected value optimization of the performance measures. However, α-quantile optimizations reached regions with greater α-quantiles (figure 2).

This means that the region distribution shapes were distinct. Figure 3 shows the cumulative distributions of the optimal regions reached with 0.1-and 0.9-quantile optimization. The deviation around the expected value which is about 540 ¤/ha in both cases, is very large. This highlights the significance of using quantile optimization. Indeed, optimization of the single expected value does not care about the shape of the region distribution. So, the decision maker can either reach optimal region with expected value probability of appearance close to one or zero. The perspective that the predicted decision will lead indifferently to the most likely or the most unlikely solution of optimal expected value, is not attractive. The optimization of quantile associated to the decision maker risk acceptance focus on more representative strategies and reduce the former misbehavior risk of appearance.

The results showed that the algorithm is as efficient to produce optimal region as the prior P 2. Furthermore it focuses on optimizing quantile thresholds which translate a new decision maker preference. Thus since the decision maker is able to characterize his risk acceptance, quantile optimization handled both with the efficiency of the predicted decision and the probability of its quantile of interest. That makes the quantile optimization more accurate to propose optimal decisions leading to optimal solutions according to the risk acceptance or aversion of the decision maker.

Optimal irrigation plans

Granting a large enough number of simulations, sensitive parameters are the same (figures 4 and 5). In one hand, the envelopes of optimal variables D1, D2 and I2 were thin, which underlines the large influence of these parameters. In the other hand, envelops of optimal variables D3 and I3 were wide, which underlines the small influence of these parameters. This observation was already true for the optimal decision space proposed with optimization of the expected value (figure 6). Main differences were related to T 1 and I1 parameters. The 0.1-quantile optimization led to relatively small intervals compared to the ones reached with 0.9-quantile optimization. We noticed two main limitations. The first one is that we simulated a water limited case, so that depending on the climate, the limit could have been reached or not and implied larger variability in the irrigation stopping rule decision parameters. The second one is that the decision parameter I2 was always high valued and may have been larger depending on the range we set to it. But I2 is an irrigation amount and its range variation should be discussed before any change.

Although the shapes of optimal regions were similar, the surface of the envelopes reached with quantile optimizations was larger than the one reached with expected value optimization. This observation corroborates that the quantile optimization is harder than expected value optimization.

Conclusion and perspectives

Granting that the decision maker is able to characterize his risk acceptance, results showed that optimizing quantiles leads to different optimal regions than optimizing expected value. As quantiles optimization allow the decision maker to take into account a new preference (the risk acceptance) which is not involved in expected value optimization, the new optimal regions reached by the P 2q algorithm address more accurately the decision maker preferences. Although we noticed that quantile optimization requires an extra sorting process and thus is more time demanding, it adds significant value to the design of irrigation strategies.

Though quantile optimization is commonly used in economic fields, it is unusual to use it to design agricultural strategies. This study emphasized the usefulness of quantile as an optimization criterion. This value should even be increased if combined with multiobjective optimization procedures. It seems indeed very attractive to be able to maximize simultaneously a small and a large quantile orders. It means that we could maximize the yield of poor crops and maximize the yield of efficient crops at the same time, such that the algorithm reaches decision subspace that avoid the small pay-off while maximizing the large ones.

4. 1

 1 MODERATO Simulator MODERATO (Bergez et al 2001) is a plot-based annual-based model aimed at evaluating current irrigation strategies for corn and at proposing improved strategies. It has been developed in close collaboration with irrigation advisers. It combines a dynamic and biophysical corn crop model with a dynamic decision model. The crop model is described in

Fig. 2

 2 Fig. 2 Two times 50 strategies of four different regions.

Fig. 3

 3 Fig. 3 Cumulative distributions based on 50 × 49 performance measures for each region.

Fig. 4

 4 Fig.4Input variables envelope of the 10 optimal regions reached for the 10 replications of the P 2q algorithm: 0.1-quantile optimization.

Fig. 5

 5 Fig.5Input variables envelope of the 10 optimal regions reached for the 10 replications of the P 2q algorithm: 0.9-quantile optimization.

Fig. 6

 6 Fig.6Input variables envelope of the 10 optimal regions reached for the 10 replications of the P 2 algorithm: expected value optimization.

Table 1

 1 General description of the strategy simulated.

	Operation	Rules
	Sowing	Sowing is between 20 April and 30 May as soon as the cumulative rainfall
		during the previous 3 days is less than 15 mm. Variety Cécilia is sown at 80
		000 plants/ha. Cécilia is a late growing variety requiring 1045 accumulated
		thermal units (ATU) from sowing to flowering and 1990 ATU from sowing
		to maturity (35% grain humidity).
	Fertilization	A single application of 200 kg/ha of nitrogen is made at sowing.
	Harvest The crop is harvested when grain moisture content reaches 20% or accu-
		mulated thermal units since sowing reach 2100 ATU and if the cumulative
		rainfall during the previous 3 days is less than 15 mm. In any case, the
		crop must be harvested before 15 October.
	Irrigation
		Sowing Irrigation to facilitate plant emergence (caused either by dryness
		or crust created by heavy rainfall on silty soil) is not taken into ac-
		count, nor irrigation to dissolve fertilizer.
		Starting irrigation Part of the optimization process.
		Next irrigation round Part of the optimization process.
		Delay irrigation Precipitation delays irrigation. When the cumulative
		rainfall over the 5 previous days is more than 10 mm, one day de-
		lay is applied for every 4 mm. The delay cannot exceed 7 consecutive
		days.
		Stopping irrigation Part of the optimization process.


Gabriel Sargent
email: gabriel.sargent@irisa.fr

Frédéric Bimbot
email: frederic.bimbot@irisa.fr

Emmanuel Vincent
email: emmanuel.vincent@inria.fr

A MUSIC STRUCTURE INFERENCE ALGORITHM BASED ON MORPHOLOGICAL ANALYSIS

Music structure refers to the description of the long term organization of a music piece through a sequence of structural segments. A structural segment can be defined by its structural borders (a start time, an end time) and a label reflecting the similarity of its music content compared to the other segments'. Its duration is typically around 16 s and more.

This document presents the music structure estimation system submitted to MIREX's structural segmentation task in 2012. It is composed of three steps : feature extraction, structural border estimation and segment labeling. First, the system produces a sequence of chroma vectors [6] expressed at the snap scale [1] (section 1). This sequence is used to calculate a segmentation criterion based on a morphological model of the structural segments [2] (section 2.1). The structural border estimation is performed by searching the segmentation with lowest cost, which combines this criterion and a regularity constraint (section 2.2). The segments are then labeled by clustering according to their similarity, through the minimization of an adaptive model selection criterion (section 3).

FEATURE EXTRACTION

The extraction of the sequence of chroma vectors of size 12 used to describe the music content of the piece is performed by means of the "Chroma Toolbox" by Muller and Ewert [START_REF] Müller | Chroma Toolbox: MATLAB implementations for extracting variants of chroma-based audio features[END_REF]. We use the CP features regularly and a hop of 0.1 s.

Then, they are expressed at the snap scale. The snap is here defined as the multiple of a beat whose period is closer to 1 s. The snap scale is synchronous to the downbeat scale, and they are often equal in practice. The beat and downbeat estimations are performed thanks to the MATLAB implementation by Davies et al. [START_REF] Davies | Towards Automatic Rhythmic Accompaniment[END_REF][START_REF] Stark | Real-Time Beat-Synchronous Analysis of Musical Audio[END_REF]. The downbeat estimator is tuned so as to consider 4 beats per bar.

We associate to each snap the mean of the CP features contained in the window centered on the snap that lasts the This document is licensed under the Creative Commons Attribution-Noncommercial-Share Alike 3.0 License. http://creativecommons.org/licenses/by-nc-sa/3.0/ c 2012 The Authors. duration of the snap period.

STRUCTURAL BORDER ESTIMATION

Morphological model

We assume that a structural segment can be characterized by its inner organization, according to its musical layers (timbre, harmony, melody ...). In this scope we consider the system and contrast model by Bimbot et al. [START_REF] Bimbot | Semiotic structure labeling of music pieces : concepts, methods and annotation conventions[END_REF]. It considers that each structural segment aimed is built from an group of typically four morphological elements of four snaps, we note {a 1 , a 2 , a 3 , a 4 }. The three first elements are related by simple transformations f and g so as a 2 = f (a 1 ) and a 3 = g(a 1 ). The fourth element can either follow the logic of the three elements and then form a system (a 4 = f (g(a 1 ))) or on the contrary contrasts with it (a 4 = δ(f (g(a 1 ))) where δ = id). Note that we assume that the relevant layers for structural analysis can vary from one structural segment to another.

However, in much cases, either f = id or g = id, or both. This leads to observe usual morphological motives like aaaa, abab, aabb in the case of systems with no contrast, or aaab, abac, aabc in the case of systems ending with a contrast. These motives can be extended to the case where the identity function id is replaced by "close to identity" functions id ′ (aaa ′ b, aba ′ c, aa ′ bc, ...). More information on this model can be found in [START_REF] Bimbot | COMPLEMENTARY REPORT TO THE ARTICLE "Semiotic structure labeling of music pieces : concepts, methods and annotation conventions[END_REF].

Segmentation criterion

The aim is to evaluate for each time unit considered the likelihood that it corresponds to the beginning of a system. We assume that at least one of the relations (f or/and g) between the elements of a system equals the identity function. For each snap t ∈ [1, T ] of a music piece, we consider the analysis window of size N = 16 snaps so as to consider three morphological elements starting from t (a 1 , a 2 , a 3 ), and one morphological element before this snap (a 0 ) as represented in figure 1. We consider that the size of each morphological element is N m = 4 snaps. The criterion Φ we consider in this work results from the linear combination of two quantities : with λ 1 and λ 2 ∈ R + learnt on a training database 1 .

Φ = λ 1 σ System + λ 2 σ Contrast (1)
σ System (t) quantifies the likelihood that t corresponds to the start time of a system through the analysis of the similarity of the morphological elements a 2 and a 3 with regard to a 1 . Be x = {X n } 1≤n≤N the sequence of features contained in the analysis window related to snap t. Let us define X = {X 0 , X 1 , X 2 , X 3 }, with vector X i = {x 1+iNm , ..., x Nm+iNm } for i = {0, 1, 2, 3}, and

Y = [Y 1 , Y 2 ] = [(X 2 -X 1 ) 2 , (X 3 -X 1 ) 2 ]
. X i contains the sequence of features of a i , and Y j is the squared distance between the features of a j+1 and a 1 , for each of its dimensions. We have :

σ System (t) = Nm j=1 min(Y 1 (j), Y 2 (j)) ||X 1 || 2 (2)
where ||X i || corresponds to the l 2 norm of vector X i .

σ System (t) allows to evaluate the contribution of X 1 to explain either X 2 or X 3 according to the various dimensions of the features. A high contribution implies that a system is likely to begin at snap t in the piece.

In the scope of the system and contrast model, a structural segment is likely to begin at t if a 2 and/or a 3 is similar to a 1 . This is considered through σ System (t). If the preceding structural segment is different from the current one, the third and the fourth morphological elements differ from a 1 . If the two segments are the same, then the third element may be similar to a 1 , but the fourth one generally differs from it. We therefore introduce σ Contrast (t) which evaluates the dissimilarity between the morphological element preceding snap t we note a 0 , and a 1 . We choose to formulate it as follows :

σ Contrast (t) = cotan(X 0 , X 1 ) (3)
where cotan(X i , X j ) corresponds to the cotangent of the angle between vectors X i and X j .

Regularity constraint

We consider a regularity constraint in the structural border estimation to favor segmentations with segments close to a typical segment size or structural pulse τ = 16 snaps . Let m be the size of a structural segment :

Ψ α (m) = | m τ -1| α (4)
with α ∈ R + a factor which controls the convexity of this function. The use of a non-convex function will favor segmentations with a majority of segments of size equal to τ and few segments whose size is far from this structural pulse. On the opposite, a convex function loosen this constraint.

Performing the structural border estimation

The segmentation criterion and the regularity constraint are combined through a linear combination to form a segmentation cost C :

C = (1 -λ 3 )Φ + λ 3 Ψ α (5) with λ 3 ∈ [0, 1].
We use the Viterbi algorithm described in [START_REF] Sargent | A regularityconstrained Viterbi algorithm and its application to the structural segmentation of songs[END_REF] to find the segmentation with the lowest segmentation cost.

SEGMENT LABELING

The labeling of the obtained segments is performed by the method described in [START_REF] Sargent | A music structure inference algorithm based on symbolic data analysis[END_REF]. Here, we transform the chroma sequence in a symbolic sequence by means of vector quantization, with the number of chroma clusters empirically fixed to 16. The Edit Distance on the symbolic features used to compare the content of the segments is replaced by a stripe distance [START_REF] Paulus | Music Structure Analysis Using a Probabilistic Fitness Measure and a Greedy Search Algorithm[END_REF] on the corresponding numeric chroma features. As we consider that the content of the fourth morphological element can be very variable as in section 2.1, we only consider the three fourth of the segment only when it lasts 16 snaps or more.

Figure 1 .

 1 Figure 1. Analysis window used for segmentation criterion calculation, containing the sequence of features x. It is composed of 4 small windows, each one related to the position of a morphological element a i of size 4 snaps, i = {0, 1, 2, 3}.

RWC Popular database with structural annotations from[START_REF] Bimbot | Decomposition into autonomous and comparable blocks : a structural description of music pieces[END_REF] were used for parameter tuning.

http://www.quaero.org/

ACKNOWLEDGEMENTS

The authors would like to thank Matthew Davies for sharing his beat and downbeat estimation scripts with us. This work was partly supported by the Quaero project 2 funded by OSEO.