Marie Clement-Rastello

A study on the size of snow particles in powder snow avalanches

In this work, we study the size of the particles involved in a powder snow avalanche phenomenon. To determine these sizes, we study all the phenomena the particles have to face before arriving in the "body" of the avalanche. We study the boundary layer which is at the bottom of the avalanche. We determine, with the help of experimental data, the range of size of the particles that can be entrained by the avalanche. We then examine the possibility for these particles to reach the top of the boundary layer, and so to take part in the avalanche. Our final result is that the more frequent particles suspended in a powder snow avalanche have a size lower than 200µm.

Introduction

Powder snow avalanches are very destructive and very ill-known phenomena. In this article we try to study them with a "microscopic" scale point of view. This kind of approach has been recently used on experimental sites where people try to know not only the macroscopic but also the inner properties such as pressure, velocity, etc. of avalanches (see Dent and others (1998), Nishimura and others (1997), Qiu and others (1997)). The purpose of this study is to know a bit better the size of the snow particles that are encountered in a powder snow avalanche. This has a real importance for someone who wants to study the interaction between particles and the turbulence of the flow. Indeed it is important to know the size of the particles that are involved in the avalanche so as to compare it with the characteristic scales of the turbulent flow, to determine the difference of velocity between air and particles. The study presented here is a theoretical study in which we try to find out which sizes of particles are en-1 countered more frequently in a powder snow avalanche. The way we proceed is the following: we study separately each of the mechanisms that are encountered by a snow particle before reaching the body of the avalanche. For each of them, we look at the range of sizes of the concerned particles. At the end we compare the different ranges found in each case so as to find the range that we must consider for a powder snow avalanche. The study is twodimensional and deals neither with the front part of the avalanche nor with its tail but with its body. A pure powder snow avalanche is studied going its way on the snow cover. The different phenomena that are considered are: the size of the particles in the snow cover, the picking up of the particles from the snow cover by air friction and by collisions, and finally the proportion of each type of particles that manage to reach the top of the boundary layer, which is at the bottom of the avalanche, and so that take part to the body of the avalanche. For this study, we need experimental data. We chose to use the results obtained by Nishimura and others (1997) while making measurements in an effective avalanche.

The particles of the snow cover

A first range for the size of the particles is given by the particles which are present in the snow cover. This range of size is not a fixed parameter, it depends on many param-eters such as temperature, wind velocity, etc. during and since the snow fall. [START_REF] Mellor | Properties of Snow[END_REF] studied the size of the particles that can be found in snow layers. We use his results concerning recent snow covers. From these we conclude that the size of the more frequent particles is smaller than 1 -2 mm.

3 The pick up of the particles from the ground

Mechanisms

A particle on the ground can be set into motion either due to air friction or to the impact of another particle which is projected into the ground.

Here, the aim is to obtain the size of the potentially mobilized snow grains. Let us first look at the collisions. What we show here is that the collisions between particles and the snow cover can increase the range of size with small radius but not really with large radius. Indeed let us first suppose the "best case" of energy transmission. One particle of radius a c arrives with a velocity u c . If it collides elastically with a particle of the snow cover of radius a l it gives all its energy to this particle. In this case the leaving particle will have a velocity

u l = ( a c a l ) 3 2 u c . (1)
One can see that even in this ideal case to have a significant velocity the leaving particle cannot have a size much larger than the incoming one. In reality it seems clear that the transmission of energy is not so perfect and that we do not have all the energy transferred from one particle to another. So we have

u l ≤ ( a c a l ) 3 2 u c (2)
and our previous conclusion remains valid.

On the other hand the collision phenomenon does not give any restriction on the lower limit of the sizes of the particles. So when we will compare the results of the picking-up by collisions and of the picking-up by air friction, because of the collisions we will lose the information on the lower limit. For the upper limit the collisions will have no influence. Let us now look at the range we obtain from pick up by air friction.

Threshold velocity for particle pick-up

In order to determine the radius of the particles which are entrained from the snow cover by the friction of the avalanche, we use results that have been established for sediment transport by the wind and that are also used in snow transport by the wind. Here we use the results given by Iversen and White (1982) (see Fig. 1). To collect their data they placed samples of material on a wind tunnel floor and determined optically the threshold friction velocity (u * t ) for each of them. The particles they used were of different kind and had different shapes. From their experimental re-

A =              0.129 1+ 6 * 10 -7 ρ p gd 2.5
1.928B 0.092 -1 , 0.03 ≤ B ≤ 10 0.120 1 + 6 * 10 -7 ρ p gd 2.5

(1 -0.0858e (-0.0617(B-10)) ), 10 ≤ B

(3)

with A = u * t ρ ρ p gd and B = u * t d
ν , where d is the diameter of a particle, A is the ratio of the threshold friction velocity to the shallow water velocity of a density current and B is a Reynolds number. Iversen and White give a lower limit (B = 0.03) but no upper limit for their law. The lower limit is outside the range of our study so does not limit us. Concerning a possible upper limit we plot on the figure the curve which separates the region in which Iversen and White made measurements and that where they did not make any. They say that for particles larger than 40 µm the above empirical law is generally within 5% of the experimental values. For smaller particles it is more difficult to perform the experiments and so the precision is lower. This is of no importance here because due to the collisions, information on the lower limit is lost. As usually done in snow transport studies we use this law here. One restriction would be that in their experiments Iversen and White had no adhesion between the particles, as we can have with snow. We suppose that for powder snow such bounds are not predominant. Another restriction is that the experiments have been performed each time with one given size of particles, so that the different phenomena that can occur when there are both small and large particles are not taken into account. A further major restriction is that Iversen and White's results were established in clean air on a horizontal floor, whereas in our case the air is loaded with particles and the ground can be inclined. From a recent study by [START_REF] Bintanja | The Interaction Between Drifting Snow and Atmospheric Turbulence[END_REF] it is inferred that the presence of particles modifies some characteristics of the boundary layer. However there are insufficient results that could be taken into account in a modified boundary layer model. So we use a clear air velocity profile essentially for lack of better knowledge.

Determination of the friction velocity in a powder snow avalanche

Now that a relation between the threshold friction velocity and the diameter of the potentially mobilized particles is available, the friction velocity in a powder snow avalanche must be determined. The avalanche studied by Nishimura and others was a "mixed avalanche": there was a dense part and above it a powder snow part. We suppose that the internal structure of the powder snow avalanche is not modified by the presence of the dense avalanche and use the results of the measurements for the pure powder snow. The velocity of the front part of the avalanche was from 40 m.s -1 to 60 m.s -1 . At a given point, the static pressure was measured and compared to the static pressure measured somewhere outside the influence of the flow. From this the velocity at the measuring point was measured and for the part of the avalanche behind the front a value of 10 m.s -1 was found. The point where the measurements were taken was located less than seventy centimeters above the bottom of the powder snow part; this is in the boundary layer of the powder snow avalanche. Let us suppose that the boundary layer is of the order of one meter thick. Later on we will see that the exact value of the size is of no importance. Thus, if z ≤ 25 cm (z ≤ 0.25δ with δ = 1 m the thickness of the boundary layer), the value of the friction velocity can be directly estimated from the logarithmic law which governs the velocity in this region :

u = u * K ln( z z 0 ) ( 4
)
where K = 0.41 is the von Kàrmàn constant, z 0 the roughness height,

z 0 = C 0 (u * ) 2 2g
, with C 0 = 0.021 (see [START_REF] Owen | Saltation of Uniform Grains in Air[END_REF] and [START_REF] Rasmussen | Wind Tunnel Observations of Aeolian Transport Rates[END_REF]). There are, theoreti- cally two possible values for u * (see Fig. 2).

In comparison with the values of u * found in snow transport by wind, with velocities comparable to those of avalanches, it appears that the lower value is the physically relevant value. We deduce from Fig. 2 that in this region max(u * ) ≈ 1 m.s -1 .

In the upper region (25 cm ≤ z ≤ 70 cm) the boundary layer profile relating velocity and height is unknown. If the measurement has been done in this region, then u(zǫ[0; 0.25]) ≤ 10 m.s -1 . So with the formula relating u and u * in the logarithmic region, one necessarily has u * ≤ 1 m.s -1 . Thus here again max(u * ) = 1 m.s -1 .

Let us now suppose that δ = 1 m. With the same reasoning as before one can show that if u has been measured in the upper region the deduced value of u * is lower than it would be if u had been measured in the logarithmic region. Because in the logarithmic region u * is decreasing with z (see Fig. 2), for any δ, for z ≤ 0.25δ the value of max(u * ) remains 1 m.s -1 . Thus the size of the boundary layer does not change anything in this estimate. Two things should be improved in the future : first a better determination of u * with more field data, and second to account for the presence of the particles in the boundary layer in the determination of the velocity profile u(z).

Size of the entrained particles

We are now able to determine the maximum range of sizes of the particles which can be picked up from the snow cover by air friction : it is [5 * 10 -6 m; 10 -2 m] (see Fig. 1).

To determine these values of the size (d) of the particles, formulas given by Iversen and White were used (see section 3.2). The upper value of d still lies in the domain of validity of their formulas but is outside the region in which they performed the experiments. Because they did not put any upper limitation to their formula and because our range is not far from the range they studied experimentally we think their formula is applicable for the upper value of d.

The upper limit of this interval is not changed by the collisions of particles with the snow layer, but the lower limit can be changed by the impacts (see section 3.1). Considering both air friction and collisions, the size of the mobilized particles is therefore d ≤ 10 -2 m.

Sedimentation in the boundary layer

We now study the sedimentation of the different particles in the boundary layer. The purpose is to find out which particles "fall" sufficiently slowly to be present at the top of the boundary layer, and so to take part in the avalanche.

Law relating volumetric concentration and height

To find the governing law for the volumetric concentration in the boundary layer, we use the equations of mass conservation for the mean flow of air and the equation of the mean conservation of particles

∂u j ∂x j = 0 (5) ∂φ ∂t + ∂(φu pj ) ∂x j = ∂(D t ∂φ ∂x j ) ∂x j (6)
where u j are the components of the mean velocity of air, u pj the components of the mean velocity of the particles, φ the mean volumetric concentration, and D t the turbulent diffusivity that comes from the turbulent closure model:

φ ′ u ′ ip = -D t ∂φ ∂x j . () ′ are turbulent quantities.
We suppose that the mean flow is steady and that the velocity and concentration profiles are independent of x. With the boundary condition that near the ground u z = 0 (u z is the velocity perpendicular to the slope), we then have u z = 0 in the whole boundary layer. Boundary layer theories imply that D t = Ku * z σs , where σ s is the Schmidt number. The Schmidt number is the ratio σ s = Dt νt between the diffusivity of snow particle and the eddy diffusivity of momentum. Some more advanced theories [START_REF] Bintanja | The Interaction Between Drifting Snow and Atmospheric Turbulence[END_REF]) take into account that the suspended particles modify the turbulence of the air, so that the diffusive coefficient D t is modified :

D t = Ku * z σ s (1 + AR) (7) with R = -g( ρp ρ -1) ∂Φ ∂z σ s ( ∂ux ∂z ) 2 and A ≈ 5 -7
where R is the flux "snowdrift" Richardson number, ρ p is the density of particles, ρ is the density of air and g is the gravitational acceleration. The order of magnitude of the different quantities are the following:

[Φ] ≈ 10 -2 , [z] ≈ δ ≈ 1 m, [u] ≈ 50 m.s -1
and [σ s ] ≈ 0.5-1 (see section 4.2). From this we can deduce that [1 + AR] ≈ 1. So within the precision of our study, this modification need not to be taken into account.

To solve the second differential equation, we suppose that the particle volumetric concentration and its gradient vanish infinitely far from the boundary layer lim z→+∞ φ = 0 and lim z→+∞ ∂φ ∂z = 0 (8)

The steadiness of the flow, the independence of the concentration profiles of x and the fact that u z = 0 imply, using equation ( 6)

∂(-Φv r cos(α) + D t ∂Φ ∂z ∂z = 0 (9)
with α the angle of the slope and v r = | -→ u p --→ u |. From this and with the boundary conditions (see equation 8) we deduce

φ = φ 1 ( z z 1 ) -γ with γ = v r cos(α)σ s Ku * . ( 10
)
φ 1 is the volumetric concentration at z = z 1 . Let z 1 be the roughness height, so that φ 1 is the volumetric concentration of the particles near the ground, i.e., the volumetric concentration of the mobilized particles. v r is a function of the size of the particles, so the variation in the boundary layer of the volumetric concentration will not be the same for different types of particles.

Presence of the particles

To know which particles effectively reach the top of the boundary layer, we study the ratio between the number of particles at the top of the boundary layer and the number of them near the ground; to this end we look at the repartition with the radius a defined as

R(a) = φ φ 1 = ( z z 1 ) -γ . ( 11
)
It defines the concentration of particles by which the body of the avalanche is fed from the boundary layer. To find it the relative velocity v r must be determined as a function of a. We assimilate it with the sedimentation velocity. The sedimentation velocity is obtained from the Lagrangian equilibrium state where the gravity force is balanced by the drag force. From this equilibrium we obtain

C d Re 2 p = g(ρ p -ρ) 32a 3 ρ 3µ 2 , (12)
in which C d is the drag coefficient and Re p = 2aρvr µ is the particle Reynolds number. For spherical particles, as we suppose here the snow particles in the flow to be, the drag coefficient can be given as follows, see Niño and Garcia (1994)

C d = 24 Re p (1+0.15 Re p +0.017Re p )- 0.208 1 + 10 4 √ Rep .
(13) Another parameter is the Schmidt number (σ s ) which takes different values in different configurations. [START_REF] Householder | Turbulent Diffusion and Schmidt Number of Particles[END_REF] found that for particles which are denser than the ambient fluid, the Schmidt number is smaller than 1. [START_REF] Naaim | Experimental and Theorical Determination of Concentration Profiles and Influence of Particle Characteristics in Blowing Snow[END_REF] performed wind tunnel experiments with PVC particles and found a Schmidt number of 0.5 to 0.6. We performed calculations both for σ s = 1 and σ s = 0.5. The last parameter to fix is u * . As before, because we look for the larger range of particle sizes, we take u * = 1 m.s -1 . The results are reported in Fig 3 . To demonstrate that the dependences on σ s and on the thickness of the boundary layer are weak we plotted the curves for either σ s = 1 and σ s = 0.5, and either δ = 50 cm and δ = 5 m. In all cases the particles which, in an important number, are

Conclusion

The purpose of our study was to estimate the size of the particles which are present in a powder snow avalanche. They correspond to the particles that reach the top of the boundary layer, between the ground and the body of the avalanche. Indeed in this way they can take part in the processes in the body of the avalanche. We analyzed the different phenomena which occur in this boundary layer and their influence on the size of the particles. Thus, we studied consecutively, the sizes of the particles of the snow cover, the pick-up of particles by the avalanche and the selection which is done between the particles while going up in the boundary layer. With this analyze, we were able to determine the maximum size of the particles which take part in the dynamics of a powder snow avalanche. These particles are those which have a size smaller than 200µm. This knowledge will allow to use the results, that have been established in diphasic flow studies, about the interaction between particles and turbulence. It will also help to improve numerical models which use as a parameter the size of the particles of snow. Finally for laboratory experiments, this information is necessary for the choice of particles. Another step, now, is to capture snow particles during the flow of a powder snow avalanche, and measure their size. This will help to improve our understanding of the processes and our theoretical model.

Figure 1 :

 1 Figure 1: Threshold friction velocity versus diameter of the particles

Figure 2 :

 2 Figure 2: Determination of u *

Figure 3 :

 3 Figure 3: Repartition for σ s = 1, σ s = 0.5, δ = 50 cm and δ = 5 m

* This work was supported by the Pole Grenoblois sur les Risques Naturels


Sijie Wu

Jianxi Huang

Xingquan Liu

Guannan Ma

Forest over classification from multi-temporal MODIS images in Southeast Asia using decision tree

Keywords: Multi-temporal, forest cover, decision tree

Abstract:MODIS data is of significant for the classification of regional forest cover due to its high temporal resolution and high spectral resolution. Forest cover is an important parameter for forest ecosystem. The objective of this preliminary study is to mapping forest cover from mutli-temporal MODIS data with decision tree. The classification forest samples were selected from four global land cover datasets with specific rules. The selected samples were used to generate rules of the decision tree for the classification of forest cover. The study results show that mutli-temporal remote sensing data with decision tree method have great potential to improve the regional forest cover mapping.

Introduction

Identification of types of forest has significant for forestry source monitoring and management. Because of the capability of acquiring regional surface information, remote sensing has become a reliable tool for identifying types of forest in regional and global scales. Currently, most of the applications of remote sensing classification are the traditional statistical pattern recognition method, such as minimum distance, parallelepiped, maximum likelihood, and mixed-distance method, cyclic cluster method and other supervision or unsupervised classification method. Because of the existence of spatial resolution of remote sensing image itself and "same object with different spectrum‖, -different objects with same image" phenomenon, misclassification and leakage of points occur more frequently. These factors lead to the low classification accuracy. New methods of pattern classification are as follows: fuzzy classification, classification based on texture description of Markov random field model, classification of wavelet analysis, fractal texture method, neural network and expert system classification, etc [START_REF] Huang | Remotely Sensed Image Classification Technology Review[END_REF][START_REF] Li | Study of Remotely Sensed Image Classification Technology[END_REF][START_REF] Wu | Purification of Training Samples in Supervised Classification of Remote Sensing Data[END_REF]. Currently, remote sensing information composite technique is widely used [START_REF] Chen | Remote Sensing Analyst[END_REF][START_REF] Zhao | Principles and methods of analysis of remote sensing applications[END_REF]. In recent years, researchers primarily utilize satellite remote sensing with vegetation surface temperature, terrain elements and other non-remote sensing of forest vegetation on the ground to identify sub-categories [START_REF] Liu | Research Advances and Countermeasures of Remote Sensing Classification of Forest Vegetation[END_REF].

In the past decades, applied research results on the large area forest cover mainly using AVHRR data has achieved great success [START_REF] Hansen | Global land cover classification at 1km resolution using a decision tree classifier[END_REF][START_REF] Kressler | Detecting land cover changes from NOAA 2AVHRR data by using spectral mixture analysis[END_REF]. Due to limitations AVHRR data for land cover mapping applications, there are still many uncertainties [START_REF] Loveland | An analysis of the IGBP global land cover characterization process[END_REF]. Thus, with the launch of Terra satellite, the use of MODIS data in regional scale studies of forest cover has been developed. Using MODIS global supervised classification model,Muchoney [START_REF] Muchoney | Application of the MODIS global supervised classification model to vegetation and land cover mapping of Central America[END_REF] classified the vegetation and land cover in central United States. On this basis, with IGBP classification system and the STEP global plots database and MLCCA (MODIS land cover classification algorithm method. Friedl [START_REF] Friedl | Global land cover mapping from MODIS: algorithms and early results[END_REF] carried out a global land cover classification with a total of five months of MODIS data in 2000. Based on pattern decomposition method (PDM), Cen [START_REF] Yi | Land cover classification in Kii Peninsula of Japan using Terra/MODIS data[END_REF] conducted a study of land cover classification using MODIS data of the Kii Peninsula of Japan in 2001. Using MODIS 8 day composite reflectivity products, Carrã o [START_REF] Carrã O | Contribution of multispectral and multi-temporal information from MODIS images to land cover classification[END_REF] evaluated the efficiency of MODIS hyper spectral data and relative land cover classification for a long time. On the basis of studying and comparing the abroad classification algorithm that having good application effects, Wu [START_REF] Wu | MODIS Data Land Use Classification[END_REF] achieved quantitative judgments on continuous coverage of MODIS data and generated the land use status classification. Liu [START_REF] Yonghong | Regional Land Cover Image Classification and Accuracy Evaluation Using MODIS Data[END_REF] proposed a classification method using MODIS data to select and extract classification feature and do large area land use/cover classification combined with the multi-temporal characteristics. Classification test was conducted in Shandong Province in China.

As for classification data, MODIS data can provide more data products. Visible data and near infrared data of the MODIS can response the growth characteristics of different vegetation types in different periods well, it is perfectly suited for forest classification. Taking into account of close contact of the distribution of vegetation with the climate and soil, coupled with influence of the climate by the altitude, slope and other terrain factors, it is necessary to add the soil and terrain data in the classification of forest types.

Study area and data

There are 11 countries in Southeast Asia: Vietnam, Laos, Cambodia, Thailand, Myanmar, Malaysia, Singapore, Indonesia, Brunei, and the Philippines, Timor-Leste. The special geographical location makes the Southeast Asia having hot and humid climate and lush tropical forests. Wet equatorial climate and a tropical monsoon climate are two types of the Southeast Asia. The main natural vegetation here is the tropical rain forest and tropical monsoon forest. Figure1 show the study area of Southeast Asia.

It can be divided into two sub-areas: In Southeast Asia, land surface is mainly composed by woodland, grassland, residential areas, etc., among which forest cover is the most typical one. In order to improve classification accuracy in the region, a lot of land surface sample data is needed. Because the large scope of the study area and too many type and the far distance, it is difficult to achieve the field sampling of vegetation types. Through data analysis on the existing global surface types, we found the land cover of the typical areas changes less, and many classification products can separate the type of surface area of typical areas very well. Therefore, in areas where vegetation integrated we reference map and interactively select the sample data, and using ROI to define the surface type. Decision tree is a tree structure similar with binary tree or multi-branches tree. The decision tree algorithm is one of the inductive learning algorithms. Induction is getting the regular conclusion by contrast, analysis, summary and summarization from characteristics and attributes of examples that known and seemingly with no specific order and rules. This conclusion was used to generate decision tree as the rule of interpretation, then classification was carried out on cases out of samples. Early famous decision tree algorithm was ID3 posed by Quinlan [START_REF] Quinlan | C4.5 Programs for Machine Learning[END_REF] in 1993. This is one of the most basic decision tree algorithms. C4.5 algorithm was posed by Quinlan [START_REF] Mehta | SLIQ: A fast sealable classifier for data mining[END_REF] in 1993 as Subsequence of ID3 and became basis of later decision tree algorithms. SLIQ (supervised Learning in Quest) algorithm is a fast scalable decision tree algorithm posed by IBM researchers that suitable for processing the data of large scale. CART decision tree, also called classification regression tree, was a data survey and prediction algorithm posed by Breiman,Friedman,Olshen and Stone in 1984.

QUEST is a new decision tree construction algorithm is formulated by Loh and Shih in 1997 The algorithm separates the variable selection and split point selection, not only suitable for continuous variables, but also for discrete variables, and it overcome the other decision tree algorithms that prefer to choose those predictor variables with more potential split point for classification. So basically it is unbiased on variable selection. And through hyper plane constituted by more than one predictor variable, this decision tree can distinguish between multiple members and non-class members of the class in the feature space.

Studies have shown that QUEST decision tree classification computation's speed and accuracy are better than other decision tree construction method.

Material and Method

Data

MODIS standard products which produced by this project were using to do forest classification. Taking into account the classification of forests and other land cover classification using vegetation status, MODIS 13Q products were downloaded, including visible light reflectivity. Considering the region's distribution of the forest influenced by the soil type and altitude, we collected global 10km soil data and STRM-DEM data and extracted the slope information. In order to maintain the consistency of the above data, the data project coordinate system and projection were normalized. Although the MODIS data was standardization, there are still clouds of pollution. Using the SG remove cloud algorithms, NDVI-SG Index data was generated. As the greenness of forests changed little during the growth in the 16 days, the index data was synthesized, that is, two 16-day month synthetic data was made into the index data. The 12th monthly indices data was instead by 23rd layer data. By integrating the above data, the input data for classification software was ready.

When carrying out the classification of the data by adding soil analysis, marked polygons were found which classification results were not the same size as same as grid soil data. After analysis, it was found that this is due to their resolution (the MODIS 1/10) was inconsistent with the image data. In the later experiments, the soil data were excluded. Taking into account the amount of data and storage limits, in order to support the Bit types of data files, when storing data, using altitude / 20 which altitude below 5,000 m in the region. Degree of slope data used (Degree) units. NDVI data was set to 0 when less than 0; Greater than 0, the draw ratio is 250 times.

Forest Classification Method

Based on four global land cover datasets, the samples with pure pixels were extracted. MERIS with 300 meters was used for validate the accuracy of percentage of each sample. And then the overall classification accuracy was discussed.

The choice of sample set of forest types

Five global Land Cover datasets, including UMD, IGBP-DIS Cover, MODIS, GLC2000 and MERIS 2005 Global Land Cover data were collected in our study. 1 km land cover map of the first four class of Southeast Asia was encrypted cut for classification of 500 meters. Combined with GLOBCOVER, a standard reference data for classification in this region was formed. Because these data in time is quite different and system of classification is different, classification system of forest mapping in Southeast Asia which was provided by Chinese Academy of Forestry was using for combining land cover types, encoding conversion and forming a comprehensive map of forest types. When producing the integrated forest type's graphs, the types with relatively large probability of appearance were on behalf of the forest types in the region, after the adoption of uniform codes. In order to analyze differences in forest types, using red color to represent the type inconsistent, the white for the same type. Types can be seen change little from on higher ground area; types change more in the flat terrain and more human activities places. Using ENVI-ROI function, the sample data were collected in each typical area. Based on the sample analysis and the decision-making algorithm, regression analysis of data for decision-making was conducted to form the classification decision in the region. Figure2, figure3, figure4 and figure5 and figure 6 represent IGBP-DIS, MODIS, GLC 2000, UMD and MERIS forest cover map respectively. And figure7 show the integrated forest cover results.

Decision Tree Classification Algorithm

The key of the decision tree algorithm is to establishment decision rules through the sample points. There are several decision-making algorithms to establish the rules, but algorithms for massive remote sensing data classification have deficiencies. Especially there are few algorithms suitable for application in forest cover classification. Through the relevant literature and data on foreign inquiries, Quest algorithm has been selected to generate the decision tree of classification. The basic principle is as follows:

1) First, predict variables are selected. Correlation between all predictor variables X and target variable Y is analyzed in turn. If X is discrete, the association strength between X and Y will be calculated using the chi-square test respectively, and find the P value; If X is ordered or continuous variables, P value will be calculated using analysis of variance.

2) All the P values are compared to the value preset boundary (the default is 0.05). If they are less than the critical value, the minimum value of the variable P is selected as a branch variable; If they are all larger than the cutoff value, when X is continuous or ordered variables, P value will be calculated using Levene homogeneity test of variances, and when P value is less than the cutoff value, the minimum value of the variable P will be chosen as a branch of variables. If P values of homogeneity of variance test are large than critical value, the minimum value of the variable P in the first step will be the direct selection as the branch variable.

3) If the elected Branch variables are of discrete categorical variables, through the transformation, the goal variable Y will be maximized when X values are different. And largest discriminate coordinate will be calculated.

4) If Y is more classified categories, the mean number of X will be calculated for each type of Y value. These categories will eventually merge into two categories using the cluster analysis algorithm. This will simplify the multiclass discrimination problem to two types of discrimination.

5) Split point position is finally decided using quadratic discriminate analysis. And the selected original value of predictor variable X is got for building classification rules.

Based on the QUEST algorithm, the system combines the selection of the sample, decision analysis and classification. Firstly, using the sample data, samples were selected based on ROI selection tools in ENVI software. Then sampled were changed into data formats supported by QUEST, and QUEST algorithm interface was called by automated programs developed by IDL for establishing the decision tree automatically (suitable for established rules for ENVI decision data). Finally, classification was carried out using ENVI decision tree analysis tools to achieve judgments of the types of classification of forests in Southeast Asia. Classification accuracy and speed are greatly improved. Using QUEST algorithm, regression accuracy analysis of sample points established by decision tree is showed as follows. The accuracy is better than 75%, and the average accuracy reaches 83.02%. The above results only show that the regression accuracy of samples is high, but cannot represent the actual accuracy. The classification results are needed to be analyzed.

Sample precision analysis

Table 2 accuracy analysis of selected sample points

Results

In order to fully use the data features in the data processing, four schemes have been taken into account: scheme 1 implement the SG cloud-remove processing, and monthsynthesis used maximum NDVI values (SGMaxLANDTYPE); scheme 2 carry out the SG cloud-remove processing, and month-synthesis used average NDVI values (MeanLANDTYPE); scheme 3 without implementing SG cloud-remove processing, and month-synthesis used maximum NDVI synthesis (SGMaxLANDTYPE); scheme 4 without implementing SG cloud-remove processing, and month-synthesis used average NDVI values (MeanLANDTYPE). Data samples were collected in the same position of sampling pointon the above to construct different trees, and carry out the forest cover classification with quest decision tree algorithm. Classification results are shown in figure 8 Since 2005 global land classification data (GLOBCover2005) is the update data and its spatial resolution is better than 500m, it can be used as reference data for the accuracy verification of classification results of this project. In order to evaluate the classification accuracy, the results of 4 categories of data were compared with the GLOB Cover 2005. The accuracy of MaxLANDTYPE, MeanLANDTYP, SGMaxLANDTYPE ,SGMeanLANDTYPE are 65.05%,66.19%,76.66%, 78.14% respectively. It can be seen that high precision can be reached using SG cloud-remove processing and month-synthesis using average NDVI values (more than 78%). Meanwhile, the proportion of each category and the proportion of surface vegetation in each country were statistically. Table 3 indicates the results. With analysis on the proportion of different vegetation in different countries, the results showed SGMean LANDTYPE classification has achieved the highest accuracy.

Discussions

MODIS provide the multi-temporal forest cover information. Decision tree has advantage of understandable in structure, interpreted by its rules, calculating and precise in getting result, which lead to its high development in the field of forest cover classification. Multi-temporal MODIS-NDVI data with SG filtering have been used to the forest cover classification. We compare forest cover classification accuracy with the four schemes for time series MODIS-NDVI processing using SG cloud remove algorithm using Quest decision tree method. The results show that SG Mean LANDTYPE classification has achieved the highest accuracy. It can be seen that the decision tree algorithm with multi-temporal MODIS has great potential in regional forest cover mapping.

Fig. 1

 1 Fig. 1 Study area (Southeast Asia) In the establishment of classification algorithms, including supervised and unsupervised classification are in common use. In order to improve the accuracy of classification, decision-making regression tree algorithm of supervised classification was used in this research. It is characterized by the better use of sample data and extracting decision to create classification algorithm according to the sample data. The user can understand the decision-making rules well.Decision tree is a tree structure similar with binary tree or multi-branches tree. The decision tree algorithm is one of the inductive learning algorithms. Induction is getting the regular conclusion by contrast, analysis, summary and summarization from characteristics and attributes of examples that known and

/Fig. 2

 2 Fig.2 IGBP-DIS Cover Global land coverdata

Fig. 8 Fig. 9

 89 Fig.8 Forest classification result with SG cloud-remove processing and months synthesis using maximum NDVI values

  1. Indochina area: The climate here is tropical continental monsoon climate. The climate of Malay Peninsula is wet equatorial climate. Annual rainfall of the Malay Peninsula and the rainy coast of Indochina are tropical rain forest landscape. Indochina with dry and wet season is Tropical monsoon forest landscape. Less rainfall Interior plains and valley are savannah landscape. Indochina base is mountain mixed forest. Coast of North Bay and the Gulf of Siam is filled with mangroves.

	rain forest (after Brazil). Papua new several second is in the western Pacific after
	Indonesia's second-largest forest resources garden.
	2. Southeast Asian island district is also called the Malay Archipelago area. It
	belongs to maritime equatorial rainy climate. The Philippine islands belong to
	maritime tropical monsoon climate, mainly for tropical rain forest landscape.
	Southeast Asia with the Indonesian forest area for first is the world's second-largest

  to figure11:

	Shrub lands	0	0	0	6	3	0	1	80	0	1
	Grasslands	7	1	0	0	1	6	0	2	11 9	4
	Cropland	6	3	0	0	6	1	0	1	6	1 60
	Classifica-tion accuracy	76.62%	91.46 %	93.33 %	88.89%	92.19%	88.57 %	97.73 %	88.89 %	88.15 %	91.4 3 %
	Overall
	Classifica-					83.02 %
	tion accuracy

  Table 3 the distribution table of the proportion of different vegetation in Southeast

				Asia
	No .	Types	SGMaxLA NDTYPE	SGMeanLA NDTYPE	MaxLAN DTYPE	MeanLA NDTYPE	GLOBCo ver2005
	1	Evergreen needle leaf forest	4.10	7.44	8.89	5.37	0.97
	2	Evergreen broadleaf forest	40.63	38.61	35.73	37.29	35.00
	3	Deciduous broadleaf forest	4.95	3.13	2.75	4.87	1.80
	4	Mixed forest	1.91	1.53	1.96	1.79	6.41
	5	Shrub lands	12.40	14.00	12.07	11.28	10.69
	6	Grasslands	14.42	14.43	30.34	24.76	0.33
	7	Cropland	16.22	15.14	0.00	0.00	39.56
	8	Urban and built-up	4.46	0.00	5.11	11.31	0.15
	9	Water body	0.31	0.31	0.20	0.21	2.06
	10	Snow and ice + Bare land	0.58	5.39	2.96	3.12	3.02

Acknowledgements

This work is supported by the National Science Foundation of China (NSFC) project (NO.40901161), and Chinese Universities Scientific Fund (Project No. 2011JS142).


Farida Manseur

Nadir Farhi

Haj-Salem Habib

Jean-Patrick Lebacque

Robust adaptive strategies for the guidance of users in road networks

Keywords: optimal guidance, travel time reliability, robustness, traffic modeling, traffic control

We present an algorithm for optimal guidance of users in road networks. It is a "stochastic-on-time-arrival (SOTA)"-like algorithm which calculates optimal guidance strategies with reliable paths, for road network origin-destination pairs. Our contribution consists here in extending an existing SOTA algorithm, in order to include robustness of the guidance strategy, towards path failures.

The idea of SOTA algorithms is to calculate the maximum probability of reaching a destination node, starting from any node of a road network, and given a time budget. This calculus gives the optimal path for every origindestination pair of nodes in the network, with an associated optimal adaptive guidance strategy, and with respect to the considered time budget. The approach models very well the importance of the travel time variability in the route choice process. Indeed, by adopting this kind of guidance strategies, users may accept short deteriorations on the average travel time, if in return, they have a guarantee on the reliability of the travel time (in term of the maximum probability of reaching their destinations in given time budgets).

We propose here an extension of this approach in order to take into account the existence and the performance of alternative detours of the selected paths, in the calculus of the guidance strategy. We take into account the fact that one or many links of the selected optimal path may fail during the travel. We then consider that users may be sensitive to path changing. That is to say that they may prefer paths with efficient alternative detours, with respect to paths without, or with less efficient detours, even with a loss in the average travel time, and/or in its reliability. In order to take into account such behaviors, we propose a model that includes the existence as well as the performance of detours for selected paths, in the calculus of the travel time reliability (i.e. the maximum probability of reaching a destination node). This new way of calculating travel time reliability guarantees a kind of robustness of the guidance strategies. That is to say that, the travel time reliability associated to the obtained optimal guidance strategy is not likely to change, however associated adaptive paths change during the travel. The variation of the travel time reliability, with respect to a network structure changing, is thus improved.

We show the effectiveness of our algorithm on a grid network, where the vector of the travel times through the links of the network is assumed to have a multivariate gamma probability distribution, taking into account the correlations of the travel times. We present preliminary results of the application of the proposed guidance strategy, with comparisons to the existing algorithm. We show, in particular, the variance of the calculated optimal paths and strategies, with respect to the users' sensitivity towards the existence and the performance of alternative detours of the selected paths. The model is then able to calculate robust guidance strategies towards path failures.

Introduction

The optimal guidance problem of users in transportation networks knows a renewed interest in recent years. This is due mainly to the fast development of new technologies of information and communications. This development has made possible and necessary the development of real-time mobile and/or GPS computing applications, where adaptive algorithms give optimal solutions, updated according to the traffic conditions, and with better guarantees in terms of reliability, robustness and risk, compared to prior solutions. Optimal routing in transportation networks with highly varying traffic conditions is a challenging problem due to the stochastic nature of travel-time on links of the network. Depending on the users' preferences, an optimal route might need to consider notions of both the expected value the travel-time and its reliability (variance). This is a multi criterion optimization problem that is in general hard to solve. [START_REF] Loui | optimal paths in graphs with stochastic or multidimensional weights[END_REF] defined a formulation of routing optimality as a route with the Least Expected Time (LET).The LET problem has been well investigated and numerous efficient algorithms have been proposed to solve the different variants of the problem; see for example [START_REF] Fu | Expected shortest path in dynamic and stochastic traffic network[END_REF][START_REF] Miller-Hooks | Adaptive least expected time paths in stochastic, time varying transportation and data networks[END_REF][START_REF] Waller | On the online shortest path problem with limited arc cost dependencies[END_REF]. However, there are several frameworks in which the LET solution is not adequate, because it does not take into consideration the variance of travel-time distributions and does not give any guarantee on its reliability.

A very natural definition of a reliable optimal path is given in the Stochastic On Time Arrival (SOTA) approach, presented by Frank (1969). This formulation permits yet to obtain satisfactory adaptive solution in terms of travel time reliability. A formulation of the SOTA problem using stochastic optimal control is presented by Bertsekas (2005). [START_REF] Nie | Arriving on time problem[END_REF] formulated it as a stochastic dynamic programming problem, and solved it using a standard successive approximation (SA). However, in presence of cycles on the network, as is the case with all road networks, there is no finite bound on the maximum number of iterations required for the convergence of the algorithm. As an alternative, Nie et al. (2006) proposed a discrete approximation algorithm for the SOTA problem which converges in a finite number of steps and runs in pseudo-polynomial time. Samaranayake et al. (2012) presented a number of optimization techniques in order to speed up the computation time of the algorithm. The approach includes a label-setting algorithm based on the existence of a uniform strictly positive minimum link travel-time, advanced convolution methods centered on the Fast Fourier Transform and the approch of zero-delay convolution, and localization techniques for determining an optimal order of policy computation. However, the performance of all of these algorithms is limited by the large search space of the problem. Sabran et al. (2014) have shown how preprocessing methods can be used to further reduce the computation time of the SOTA problem. Unfortunately, the structure of the SOTA problem formulation limits the types of preprocessing methods that can be used for this problem, and prevents massive running time reductions in the deterministic case. Recently, [START_REF] Kobitzsch | Pruning techniques for the stochastic on-time arrival problem, An experimental study[END_REF] presented a novel approach to reduce the immense computational effort of stochastic routing based on existing techniques for alternative routes. In an extensive experimental study, they showed that the process of stochastic route planning can be speed-up immensely without sacrificing much in terms of accuracy.

In this paper, we focus on the reliable path problem that aims to find an adaptive optimal routing strategy, which takes into account the existence and the performance of alternative detours of the selected paths in the road network. That is to say, we assume that one or many links of the selected optimal path may fail during the travel, and that users may be sensitive to path changing. Some users may then prefer paths with efficient alternative detours, with respect to paths without or with less efficient detours.

Our approach is based on the idea of Frank (1969) considering that a reliable path from a given origin to a given destination maximizes the probability of realizing a travel time less than a given time budget. Thus, we propose an adaptation of this approach introducing robustness in the selection of the optimal path. We base here on the routing model of Samaranayak (2011). From the probability distributions of travel times through the links of the network, users evaluate their maximum probability to reach their destination in given time budgets, and through different possible routes. The model takes into account correlations between the travel times through the links of the network.

The remainder of this paper is organized as follows. In section 2, the SOTA problem and the SA method are summarized. The formulation of the problem and the description of our proposed algorithm are provided in section 3. Computational tests are conducted in section 4. In section 5 conclusions are drawn.

Summary of the SOTA problem formulation and SA algorithm

In this section, we briefly summarize the SOTA problem for the convenience of readers and the continuity of our discussion in later sections. The details of the approach summarized here are available for example in Samaranayak et al. (2011).

We consider the graph ( , ), where is the set of nodes, with | | = , and is the set of arcs, with | | = . The set of successor and predecessor nodes are denoted by Γ ( ) = { |( , ) ∈ } and Γ ( ) = { |( , ) ∈ } respectively. The SOTA problem is to find the best routing strategy from any starting node ( = 1, 2, … … , ) maximizing the probability of arriving at the destination node within a desired time . The routing procedure from a given origin to a given destination is considered as a multistage decision process, every decision stage of the routing process includes two states; the physical state which represents the location of the current decision point, and the information state which represents the time budget for completing the remaining process. The optimal decision for the remaining process is made on the basis of the two states of the current stage. This represents a typical feedback control process in which dynamic programming plays a primary role. Due to the randomness involved in the transition between steps, and contrary to deterministic processes, optimal strategy can not be identified a priori.

Given a node ∈ and a time budget , ( , ) denotes the maximum probability of reaching destination node parting from origin node , within time , under the optimal policy, where is the predecessor node from which the user is arriving, and is the realized travel-time on the upstream link ( , ). At each node , the traveler should choose the link ( , ) that maximizes the probability of arriving on time at the destination. According to [START_REF] Samaranayake | A tractable class of algorithm for reliable routing in Stochastic networks[END_REF], the maximum probabilities ( , , ) satisfy:

( , ) = max ∈Γ = | = ( -, ) (1)
∀ ∈ \{ }, ∈ Γ ( ), ∈ Γ ( ), 0 ≤ ≤ , 0 ≤ ≤ - ( , ) = 1, ∀ ∈ Γ ( ), 0 ≤ ≤ , 0 ≤ ≤ -, (2)
where denotes the travel time on link ( , ), and ( | ) is the probability distribution function (pdf) of , conditioned on . The pdf (. ) is assumed to be known and can for example be obtained using historical data or real-time traffic information.

( -, ) is the maximum probability of arriving to destination node within time -, parting from node , conditioned that the travel time on link ( , ) is .

One approach to solve this problem is to use a successive approximation (SA) method (Algorithm 1) proposed by Fan and Nie (2006). Algorithm 1 solves the system equation (1) iteratively until convergence, and gives an optimal routing policy Given initial approximations, the iterative relationships for successive approximation are given by Algorithm1.

Algorithm 1. Successive approximation algorithm (Fan and Nie (2006)) Step 0. Initialization = 0 (iteration index)

( , ) = 0, ∀ ∈ \{ }, ∈ Γ ( ), 0 ≤ ≤ , 0 ≤ ≤ -, ( , ) = 1, ∀ ∈ Γ ( ), 0 ≤ ≤ , 0 ≤ ≤ -, Step 1. Update = + 1 ( , ) = 1, ∀ ∈ Γ ( ), 0 ≤ ≤ , 0 ≤ ≤ -, ( , ) = max ∈Γ = | = ( -, ) ∀ ∈ \{ }, ∈ Γ ( ), ∈ Γ ( ), 0 ≤ ≤ , 0 ≤ ≤ - Step 2. Convergence test If ∀( , ) ∈ × [0, ], max , ( , ) - ( ,
) = 0, then stop;

Otherwise go to step 1.

At each iteration , ( , ) gives the probability of reaching the destination node from origin node within a time budget ∈ [0, ], using a path that contains no more than links, and under the optimal policy. A convergence property is proposed by Fan and Nie (2006) using the bounded monotone convergence theorem. According to this theorem, the approximation error monotonically decreases with and the solution reaches an optimal value when is equal to the number of links in the optimal path. However, optimal paths in general stochastic networks may include loops. In this case, the number of iterations required to converge is undetermined.

Robust guidance

In this article we base on the routing model (equations ( 1) and ( 2)) presented by [START_REF] Samaranayake | A tractable class of algorithm for reliable routing in Stochastic networks[END_REF], where from the probability distributions of travel times through the links of the network, users evaluate their maximum probability to reach their destination in given time budget, and using different possible routes. The model takes into account correlations between the travel times through the links.

We propose here an extension of this approach in order to take into account the existence and the performance of alternative detours of the selected paths, in the calculus of the guidance strategy. We take into account the fact that one or many links of the selected optimal path may fail during the travel. We then consider that users may be sensitive to path changing. That is to say that they may prefer paths with efficient alternative detours, with respect to paths without, or with less efficient detours, even with a loss in the average travel time, and/or in its reliability. In order to take into account such behaviors, we propose a model that includes the existence as well as the performance of detours for selected paths, in the calculus of the travel time reliability (i.e. the maximum probability of reaching a destination node). This new way of calculating travel time reliability guarantees a kind of robustness of the guidance strategies. That is to say that, the travel time reliability associated to the obtained optimal guidance strategy is not likely to change, however associated adaptive paths change during the travel. The variation of the travel time reliability, with respect to a network structure changing, is thus improved. For that, we propose to calculate for each node the maximum probability ( , ) to reach the destination node, where we take into account the case where the selected path fails before the users who selected it reach the destination node; for witch case, alternative neighboring paths are used.

( , ) denotes, as above, the maximum probability to reach the destination node, departing from node , in a time budget , and knowing that the user comes from node upstream of , and that the realized travel time from to is . However, the mathematical definition of ( , ) is different from ( 1)-( 2), and is given below.

3.1. Calculation of the probabilities ( , )

*Corresponding author. E-mail address: farida.manseur@ifsttar.fr

We consider the following notations:

∀ ∈ \{ }, ∈ ( ), ∈ ( ), ∀ ≤ ≤ , ≤ ≤ -.

We then sort, for a given ( , ) pair, the quantities ( , ) on , in a decreasing order, and denote them ( , ) in that order, i.e. B ( , ) ≥ B ( , ) ≥ ⋯ ≥ B ( , ). We then have ( , ): = ( , ), ∀ ∈ Γ ( ) = { ( , ) ∈ ( , ), B ( , ) ≥ B ( , ) ≥ ⋯ ≥ B ( , )}

Then we write the maximum probability for a user to reach destination node from node in a time budget , knowing that the user comes from node and that the travel time from to is , as follows.

( , ) = Ψ B ( , ), ∀ ≠ , ∈ Γ ( ), 0 ≤ ≤ , 0 ≤ ≤ -

with Ψ are weighting coefficients satisfying.

Ψ = 1

Where, in case where > Γ ( ) , B ( , ) = 0 for > Γ ( ) . In this case, is the number of successors to be taken account in the network, independent of . Indeed, one can distinguishes the following two cases :  Case 1 : < Γ ( ) , in wich case, only some successors of are considered in the mean.  Case 2 : > Γ ( ) , in wich case, we have ∑

In the case 2 above, nodes that have a small number of successors are penalized; they get low values ( , ). Therefore, paths that pass through these nodes i.e. paths with small number of alternatives or detours shall have low probabilities to be selected as optimal paths. One way to choose would be to take the maximum over the cardinals of the sets Γ ( ) of successors of all the nodes of the network.

where | • | denotes the cardinal of a set. In order that formula (3) has meaning, Ψ have to be chosen such that Ψ ≥ Ψ ≥ ⋯ ≥ Ψ . That is to say that Ψ decrease as B ( , ) decrease on . This dependence of Ψ on A ( , ) makes the model nontrivial. Let us notice here that if = 1 or if > 1 and Ψ = 0, ∀ ≥ 2, then (3)-(4) coincides with (1)-(2).

Calculation of successor nodes

In the calculus of ( , ) with ( 3), instead of maximizing the quantities ( , ), we propose to take a wiegted mean of theses quantities, with weights , = 1,2, … . . The optimal guidance strategy is then determined by the sequence of successor nodes ( , ) as shown in the following formula.

( , ) denotes here the optimal successor node of node for u user to reach the destination node, knowing that the user comes from node upstream of , and that the realized travel time on link ( , ) is . By taking a mean in (3) rather than the maximum (as in (1)), we do not only take into account the path maximizing the probabilities ( , ), but we also take into account the existence and the performance of alternative deviations at each node. The SA algorithm for the formulation (3)-( 5) is summarized as follow: Algorithm 2.

Step 0. Initialization = max ∈ Γ ( ) (maximum number of successor nodes to take into account) = 0

In the following, we study the efficiency of our algorithm (Algorithm 2) by applying it on the grid network of 25 nodes and 40 links, shown in Figure 1. We assume that link travel times on the network of Figure 1 are drawn from bi More precisely, we assume that the joint probability distribution of two successive links of the network is a bi variate Gamma. We base here on the bi joint PDF of two positively correlated random variables with shape and scale parameters follow.

where is the product-moment correlation coefficient of We assume that link travel times on the network of Figure 1 are drawn from bi-variate Gamma distribution. More precisely, we assume that the joint probability distribution of two successive links of the network is a bi variate Gamma. We base here on the bi-variate Gamma distribution given by Smith et al. (1982). Indeed, the joint PDF of two positively correlated random variables and , with Gamma marginal distributions with shape and scale parameters ( , ) and ( , ) respectively, was derived by Smith et al. (1982) as

if moment correlation coefficient of and , estimated from the sample data,

are the sample means and standard deviations of the variables and

time variate Gamma distribution. More precisely, we assume that the joint probability distribution of two successive links of the network is a bivariate Gamma distribution given by Smith et al. (1982). Indeed, the , with Gamma marginal distributions and ly, was derived by Smith et al. (1982) as

, estimated from the sample data, and respectively.

(7) (8) (9)

Numerical example

In the network of Figure 1, the maximum number of successors over all the nodes is equal to 2. Therefore, according to Section 3, we take here = 2. Then, we have two weighting coefficients Ψ and Ψ . We simply denote here Ψ = Ψ , and then Ψ = 1 -Ψ.

To reach destination node 25 departing from node 1, we have 154 elementary paths. We apply Algorithm 2 and derive the probabilities ( , ) for all origin nodes of the network. In order to illustrate our approach, let us consider the probability ( , ) of reaching the destination node 25 from node 3. We chose origin node 3 here because from node 3 we have two routing actions: go to successor node 4, which will give us only one routing action at the next step (go to successor node 5), or go to successor node 8, which will give us two routing actions at the next step (go to 7 or 13). Indeed, if we chose node 4 as successor of node 3, then at node 4, if link 6 fails, the whole routing strategy fails with it. However, if we chose node 8 as successor of node 3, then at node 8, if one of the links 9 or 19 fails, we have the possibility to change and take a detour using the other link; see Figure 1.

Then we consider the following scenario. Travel-times on two successor links of the network follow a bivariate gamma probability distribution, with given average travel time vector and variance-covariance matrix. We take the same average travel time (9 time units) for all the links, except link 4, for which we take a travel time of 5 time units. We take, for every pair of two successive links, the same covariance matrix (variance = 3, covariance = 1.5). We assume that we are at node 3 and that we have spent a time = 9 time units on link 3, upstream of node 3. The maximum time budget we consider here is 80 time units.

Figure 2 shows the maximum probabilities ( , 9) in function of time budget , and for different values of . Figure 3 gives the optimal successor nodes ( , 9) in function of time budget , and for different values of . The case = corresponds to the model ( 1)-( 2), while the cases < 1 correspond to the model ( 3)-( 4). We see from Figure 2 that the maximum probability ( , ) decreases as the values of decrease. That is to say that for lower values of , lower maximum probabilities ( , ) are obtained. This is because we replaced a maximum operator in (1)-( 2) by a mean value in (3)-( 4). Indeed, a user taking a lower value of , asks for more path-robustness or path-flexibility, and, in the counterpart, he loses in term of travel time reliability. The difference ( , , ) -( , , ) can then be interpreted as the price of path-robustness corresponding to a measure of it given by the difference -.

Figure 3 gives the optimal routing policy in term of the optimal successor node. We see from that figure, that for = 1, the optimal successor node is node 4, for all the time budgets up to 80 time units. Indeed, two successor nodes (4 and 8) are possible for node 3. However, since the average travel time on link 4 (from node 3 to node 4) is (5 time units) lower than the average travel time (9 time units) on link 11 (from node 3 to node 8), and since all the other parameters are symmetric in the network (average travel times, variances, etc.), and since we do not take into account path-robustness here ( = 1), it is trivial that paths passing from node 4 are better than those passing by node 8. With low values of = 0.9 (respectively 0.8, 0.7), where path-robustness is considered, we see that, with a time budget bigger than 57 (respectively 55, 53) time units, the optimal successor node of node 3 is node 8 rather than node 4, even though paths passing through node 4 have lower average travel time comparing to those passing through node 8. For = 0.9 for example, with a time budget of 57 time units, paths passing through node 8, with travel time reliability of (57,9) = 0.59 are preferable than those passing from node 4, with a travel time reliability of (57,9) = 0.97; see red and blue curves in Figure 2 and Figure 3. For = 0.9, in general, we see that with a time budget lower than 56 time units, the algorithm prefers paths passing through node 4 (see blue line in Figure 3). However, with a time budget higher than 56 time units, the optimal policy on node 3 changes, and node 8 becomes the optimal successor node (see blue line in Figure 3). That means that, node 4 which has only one successor is penalized, i.e. it gets low values ( , ). Therefore, paths that pass through that node i.e. paths with small number of alternatives or detours have low probability to be selected as optimal paths.

From Figure 2, if we want to reach the destination node within a time budget of 60 time units, then we get a travel time reliability of 0.99 with = 1, and a travel time reliability of O.63 with = 0.9. Then, the value 0.36 = 0.99 -0.63 can be interpreted as the price of robustness to pay in term of travel time reliability. On the other hand, if we want to reach the destination node with a travel time reliability of 0.60, we need a time budget of 49 time units with = 1 and a time budget of 58 time units with = 0.9. Then, the value 9 = 58 -49 can be interpreted as the price of robustness to pay in term of travel time budget. From these results we can conclude that:  If a user maximizes the travel time reliability of the paths, without taken into account their robustness, then he will choose the path that passes from node 4, because it is the one maximizing the probability of reaching the destination node 25 in the considered time budget.  If a user seeks a guarantee in terms of robustness even if he loses in terms of travel time budget, then he will choose the path that passes from node 8.

Conclusion

This article considers the routing problem of maximizing the probability of arriving on time with a new robust adaptive strategy for the guidance of users in road networks. An existing model is modified. In order to include the performance of alternative detours of the selected paths, we extended the concept of reliability by introducing a new reliability index. The modification that we made to the existing model allows the selection of an optimal path according to two criteria: the reliability of the path in term of travel time and the robustness of the path in term of flexibility (existence and performance of alternative detours).