Résumé-With the development of photogrammetry and computer vision, Unmanned Aerial Vehicle (UAV) embedded with digital compact camera plays an important role on gathering geographic information [START_REF] Sebastian Siebert | Mobile 3d mapping for surveying earthwork projects using an unmanned aerial vehicle (uav) system[END_REF]. In this context, a novel method is proposed in this article for monitoring displacement of natural environment by airborne multi-view photogrammetry. The method is implemented for calculating 3D displacement of a glacier area in one month and a half. We illustrate the application of Micmac, a powerful and free photogrammetry software, for point cloud reconstruction. With the GPS embedded on the UAV whose precision is not quiet high, we succeed in georeferencing point clouds. As we do not have any ground control point in land, two different ways of registration have been introduced to align point clouds. Afterwards, a set of tools is presented for calculating 2D displacement maps from point clouds. We introduce at last an inversion processing to calculate a 3D displacement map. The final result shows the 3D displacement of each pixel in the entire scene. Since final result is influenced by the precision of point cloud registration, this paper also highlights useful improvements to do in the future to obtain the 3D displacement with higher precision.

I. INTRODUCTION

Nowadays, the temporal evolution of glaciers is often monitored by image acquisition systems. Different observation techniques exist such as aerial photogrammetry using an airborne camera [START_REF] Bemis | Ground-based and uav-based photogrammetry : A multi-scale, high-resolution mapping tool for structural geology and paleoseismology[END_REF] and terrestrial photogrammetry using a digital camera [START_REF] Travelletti | Correlation of multi-temporal groundbased optical images for landslide monitoring : Application, potential and limitations[END_REF]. These two techniques provide 2D displacement measurement in diverse directions according to the geometry of each acquisition, while it is more useful to retrieve 3D displacement in order to interpret deformation fields.

With the aide of existing georeferenced 3D model and GPS, 3D displacement can be estimated by ground-based multi-view photogrammetry [START_REF] Travelletti | Correlation of multi-temporal groundbased optical images for landslide monitoring : Application, potential and limitations[END_REF]. However, to retrieve the 3D displacement, nearly all photogrammetric methods require Ground Control Points (GCPs) to georeference a point cloud and calculate 3D displacement. Because of the complexity of relief, it is not always easy to obtain GCPs. As a result, the retrieval of 3D displacement by photogrammetry still remains challenging so far.

This paper presents a low-cost methodology to monitor 3D displacement of glacial areas from airborne images. The method is able to detect deformation fields by using two series of images taken from different times. It is evaluated then used to calculate 2D displacement maps. The reason why we develop a virtual camera is that we can control camera acquisition positions and camera parameters including intrinsic and extrinsic parameters. As a result, ortho-images can be captured on two point clouds by using one camera at the same positions. Multiple 2D displacement maps are then calculated by multi-correlation. Finally, the 3D displacement is calculated by using inversion technique [START_REF] Yan | Fusion de mesures de déplacement issues d'imagerie sar : application aux modélisations séismo-volcaniques[END_REF] based on more than three different 2D displacement maps.

II. STUDY AREA

The Miage Glacier is one of the principal Italian glaciers of Mont Blanc massif. The lower glacier is a 8-km-long, gently sloping tongue, partly entrenched in a deep, straight valley [START_REF] Deline | Glacier fluctuations in the western alps during the neoglacial, as indicate by the miage morainic amphitheatre (mont blanc massif, italy)[END_REF]. Figure 1 shows the localization of study area and a point cloud we constructed. The point cloud of this glacier area can be manually divided into two parts : the northern part is the glacier area covered by supra-glacial debris and an ice-contact lake, the southern part is a fixed area covered by a great deal of vegetation. The whole scene is about 800 m from the east to the West and 900 m from the North to the South.

FIGURE 1: Localization of the study area and point cloud reconstruction zone

In this article, the fixed area is used to georeference and align point clouds. Several big rocks along the road were on two series of photos taken by a drone with integrated GPS. No GCP is available on the ground, only images and positions of the drone during the acquisition are available. The data are used for reconstruction of point clouds by using MicMac [START_REF] Deseilligny | MicMac, Apero and Other Beverages in a Nutshell[END_REF]. Based on two point clouds obtained from different image acquisitions, a virtual camera is created to take ortho-images from these point clouds, ortho-images are particularly selected as characteristic points. The moving area is deformation field where we calculate 3D displacement. Deformation types are not the same and depend on whether the field is glacial area or lake. Glaciers flows towards one direction, therefore its displacement is homogeneous. However, the type of change near the lake is not regular because the level of the lake brings down in summer, more edge areas show up. Since these areas appear only in the latter point cloud, the displacement is not able to calculate.

III. DATA ACQUISITION

Two series of images were taken by an UAV from the Miage glacier and lake. Between these two acquisitions, there is an interval of 7 weeks during which changes occurred. For example, the fall of ice pieces into the lake from the northern wall is a quite frequent event [START_REF] Tinti | The miage glacier in the valley of aosta (western alps, italy) and the extraordinary detachment which occurred on august 9[END_REF].

In this study the Swinglet CAM (SC) from SenseFly is used as an unmanned vehicle (UAV). The SC is a fixedwing UAV with 80 cm wingspan and a maximum take-off weight of 0,550 Kg. The nominal endurance is 30 minutes of flight time, the propeller is an electric brushless motor and the nominal cruise speed is 11 m/s (40 Km/h). The SC is mounted with a GPS receiver, wind meter and Canon IXUS 125HS 16 Mp digital compact camera, self-powered and electronically integrated and controlled. The camera is controlled by the autopilot system so as to have images at a correct position. The flights purpose was a piloting training and flight safety check up in operational conditions (altitude, wind, distance between take-off location and survey area). Therefore no GCP was collected.

The first acquisition was taken on 18th July 2014. 143 photos captured by the UAV at an altitude of 170 m above the ground were recorded in a linear way parallel to the land surface. All the images are in JPEG format with a size of 4608 * 3456 pixels. The second acquisition was carried out on 14th September 2014, we used the same UAV to take 250 pictures. Both flights have a side overlap of 80% and forward overlap of 85%. More UAV flight features is shown in table I The method of 3D displacement retrieval includes 3 parts : 3D dense point cloud reconstruction, 2D displacement detection and 3D displacement reconstruction. Details of each part can be find in figure 2.

For the surface reconstruction, Micmac is used to extract the tie-points (Tapioca), calibrate the camera (Tapas), georeference the point cloud (CenterBascule) and compute the dense matching (Malt). Micmac is a set of photogrammetric software that allows to compute a point cloud from a set of images [START_REF] Deseilligny | MicMac, Apero and Other Beverages in a Nutshell[END_REF]. By using Micmac, one point cloud was constructed for each series of images.

After the 3D reconstruction, we developed a set of tools for calculating 2D displacement maps. Firstly, georeferencing error of 3D reconstruction is corrected. Although the 3D reconstruction is georeferenced, an alignment step is required. This is mainly due to the imprecision of camera shooting position given by GPS embedded on UAV. Usually this step requires real GCPs or relative GCPs on the motion free area of one of the 3D reconstructions. Secondly, a virtual camera was created in order to take ortho-images from point clouds. As all intrinsic and extrinsic parameters can be manually set, the visual camera can be put in every location with any orientation. In this article, visual camera has been put in five different positions as shown in figure 3. Thirdly, at each position, the Once more than 3 maps of 2D displacement are calculated, a 3D displacement map can be extracted by using inversion technique [START_REF] Yan | Fusion de mesures de déplacement issues d'imagerie sar : application aux modélisations séismo-volcaniques[END_REF]. To clarify this point further, 2D displacement can be considered as projections of 3D displacement (with 3 components, East, North, Up) in different directions depending on the image acquisition geometry. Consequently, the 3D displacement can be retrieved from at least three 2D displacement maps in different directions by inversion in the least square sense. For this, projection vectors that link 2D displacement and 3D displacement are necessary. As we already have five 2D displacement maps, we are able to retrieve the 3D displacement field over Miage glacier.

V. RESULTS

A. 3D reconstruction and alignment

After the MicMac processing, we obtain two point clouds, one for each series of images. The first one is that of the Miage Glacier on 18th July 2014. For this cloud, we have about 1300 points/m 2 and we find gaps in areas covered by vegetation or water. The second one is Miage Glacier on 14th September 2014. Since we have more pictures for this series, the density is 1700 points/m 2 . As in the first cloud, there are also areas with gaps. This phenomenon is due to changes of vegetation and a lack of texture on the lake surface, both of these situations may cause bad images correlations. To georeference these two point clouds, the GPS embedded on the UAV is used. It records camera's position at each shooting time.

Then, we picked several rocks in the South part of the point cloud as reference points. These rocks are located in the fixed area and have not moved. These reference points are used to compare the difference of fixed area between both point clouds. After comparing 7 points in fixed area, differences are between 0 and 1 m in X axe, between 2 and 3 m in Y axe and between 6 and 10 m in Z axe. These results are not quite ideal because the fixed area has not moved so that differences should be close to zero. The reasons why we have a high difference between both point clouds are the following : firstly, as we have no GCPs on land surface, the georeferencing task is carried out by using embedded GPS of the UAV which has a low precision. Secondly, Micmac also produces errors in the point cloud construction procedure [START_REF] Strumpf | Ground-based multi-view photogrammetry for the monitoring of landslide deformation and erosion[END_REF].

Next, we tried two different ways to reduce differences in the fixed area between both point clouds. The first one is a point clouds alignment by using CloudCompare. To do so, we separated fixed and moving area for each point cloud and aligned only fixed areas. After comparison, differences in X, Y and Z axes are less than 0.5 m in fixed area. Then, we applied the transformation matrix onto the whole second point cloud and compared results. In the moving area, a few homologous points were picked to compare their coordinates. The result shows that the difference is less than 0.5 m in X and Y axes. However, it is about 5 m in Z axe at the border. The reason is that in moving area, no GCP is available so that the northern part of the point cloud has deviated when it was built. The second method consists in extracting GCP from fixed area of the first point cloud and to use these GCP to reconstruct the second point cloud. After reconstruction, we compared the new point cloud with the first one. The result obtained is very close to the aligned point cloud of the first method.

B. 2D displacement map

Once we have two well aligned point clouds, a virtual camera is placed into the 3D environment and takes pictures of both point clouds as shown in figure 3. Five ortho-images have been captured from each point cloud. For instance, figure 4 shows an ortho-image of the first point cloud taken in initial camera position and figure 5 shows an ortho-image of the second point cloud taken in the same camera position. Afterwards, the camera is moved to take pictures of the point cloud from four different points of view. These positions are all defined with 20 degrees of angle with the vertical axe Z (nadir to zenith) and 0, 90, 180 and 270 degrees around the Z axe. As a consequence, 5 pairs of orthoimages have been captured.

Then, image correlation is applied for each pair and we can get 5 maps of 2D displacement. Figure 6 shows one of the five maps of displacement. For this map, the camera is in the initial position (0 degree with Z and 0 degree around Z). The 2D displacement computed in pixel unit by the correlation algorithm is converted in meters using the georeferenced point cloud. This result is filtered to eliminate most of the inconsistent errors. In this figure, arrows have been added in order to show the direction and distance of the displacement of each pixel. The direction of arrows is the direction of displacement and the length of arrows stands for the distance of displacement. We can see in figure 6 that in fixed area, the length of red arrows is very short -most of them are less than 0.5 m -which means that few movement occured. These red arrows illustrate the residual error after alignment.

In the moving area, directions of blue arrows are homogeneous and meet the direction of movement of the lower part of the Miage glacier (from West to East). For the entire area, we can find that red arrows are more than 0.5 m and inconsistent blue arrows correspond to correlation errors due to the vegetation or light effects.

C. 3D displacement map

In order to compute 3D displacement for each pixel, the least squares inversion algorithm is used to fuse 2D displacement maps. The forward equation is given by : R = PU Where P is projection matrices, U is the vector of 3D displacement and R is the 2D displacement depending on projection matrices. In our case, P and R are known. P is computed according to the positions of the virtual camera and R is the result of correlation algorithm. The inverse problem -to compute U depending on P and R -can be solved by the least squares inversion algorithm. In this case, U is given by : U = (P T P ) -1 P T R, where P T is the transposed matrix of P . As a result, a 3D displacement map is calculated as shown in figure 7.

FIGURE 7: 3D displacement amplitude mapped onto the ortho-image

From this figure, we can see the area where the 3D displacement was calculated. Displacement on the fixed area is close to 0. This result highlights a correct alignment of point clouds. On the glacier, displacement is between 3 and 5 m which is over estimated, the real displacement is close to 3 m. This error is due to an error on the Z axe displacement that corresponds to a deviation of alignment on the moving area where there is no GCP.

VI. CONCLUSION

The lack of GCPs led to a over estimated displacement of two 3D point clouds. In order to obtain a higher precision, higher accuracy point clouds are needed. In general terms, the accuracy of a photo-based point cloud depends on the scale and resolution of the input images, the distribution and accuracy of the control data (ground control points, scale measurements or camera positions), the precision and distribution of the matching image points, and the network geometry, which includes the number of photos, how much they overlap and how convergent the views are [START_REF] Bemis | Ground-based and uav-based photogrammetry : A multi-scale, high resolution mapping tool for structural geology and paleoseismology[END_REF]. If the surveyed area is easy to access, high precision GCP is recommended. Otherwise, it is better to raise the precision of GPS integrated on the UAV. In case neither of these two ways are possible, then at the moment of image acquisition, we can extract control points from an existing digital elevation model.

In conclusion, this paper offers the latest information regarding the 3D displacement detection by airborne multi-view photogrammetry. The method proposed in this paper is able to identify the 2D and 3D displacement of each pixel. Compared to others, this method requires only a GPS integrated drone for acquisition and offers a high precision 3D displacement map. During the procedures of data treatment, all used software and tools are open-source or free of charge, which strongly reduces the cost of the monitoring of glacial areas or landslides deformation. What's more, the method presented in this paper can not only be applied on airborne images but also be applied on optical satellite or terrestrial images with GPS.

  c

FIGURE 2 :

 2 FIGURE 2: Method diagram to retrieve the 3D displacement measurement

FIGURE 3 :

 3 FIGURE 3: Visual camera in five shooting positions

FIGURE 4 :

 4 FIGURE 4: Virtual ortho-image of the first point cloud.

FIGURE 5 :

 5 FIGURE 5: Virtual ortho-image of the second point cloud.

FIGURE 6 :

 6 FIGURE 6: Projection of the 2D displacement map onto the ortho-image

  .

		First Flight	Second Flight
	Start time	10 h 47	8 h 22
	Duration	18 min	24 min
	number of photo Area covered	143 0,36 Km 2	250 0,36 km 2
	Distance	14,0 Km	14,1 Km
	Altitude above the ground	170 m	170 m
	Resolution	5,0 cm/px	5,0 cm/px

TABLE I :

 I

UAV Acquisition features

IV. 3D DISPLACEMENT RETRIEVAL METHOD


Pyrène Larrey-Lassalle

Antoine Esnouf

Philippe Roux

Miguel Lopez-Ferber $

Ralph K Rosenbaum

Eléonore Loiseau

A methodology to assess habitat fragmentation effects through regional indexes: illustration with forest biodiversity hotspots

Keywords: Landscape spatial configuration, Forest loss, Metapopulation capacity, Large-scale indexes, Biological conservation, Global maps

A methodology to derive worldwide regionalised fragmentation indexes is developed.

It was applied to all forest ecoregions included in the biodiversity hotspots.

The results highlighted significant intra-and inter-ecoregions differences.

These indexes can be used for land use planning or macro-scale conservation planning.

Introduction

In December 2016, the Conference of Parties for the Convention on Biological Diversity (COP13) gathered together 196 countries in Cancún with the objective of stopping the dramatic decline of terrestrial and marine biodiversity. WWF recently estimated that by 2020, at the current biodiversity loss rate, the world will have witnessed a two-thirds decline in global wildlife populations in only half a century (WWF, 2016a). And yet, of the twenty Aichi biodiversity targets (Convention on Biological Diversity, 2016), only four are partly reached or in the process of being reached, while for all others the situation still stagnates or even worse, for five of them, deteriorates (e.g. damages to coral reefs, pollution from excess nutrients, or most endangered species protection) (Secretariat of the Convention on Biological Diversity, 2014).

Increasing international trade in our globalised economy accounts for a significant share of biodiversity threats [START_REF] Chaudhary | Land use biodiversity impacts embodied in international food trade[END_REF][START_REF] Lenzen | International trade drives biodiversity threats in developing nations[END_REF][START_REF] Moran | Identifying species threat hotspots from global supply chains[END_REF]. In particular, consumers in developed countries cause threats to species through their

Author-produced version of the article published in Ecological Indicators, 2018, N°89, p.543-551. The original publication is available at https://www. sciencedirect.com Doi: 10.1016sciencedirect.com Doi: 10. /j.ecolind.2018.01.068 .01.068 imported supplies from developing countries [START_REF] Moran | Identifying species threat hotspots from global supply chains[END_REF]. Similarly, numerous Life Cycle Assessment (LCA) studies showed that complex international value chains have environmental impacts, among others on biodiversity, all over the planet [START_REF] De Baan | Highresolution assessment of land use impacts on biodiversity in life cycle assessment using species habitat suitability models[END_REF][START_REF] Hellweg | Emerging approaches, challenges and opportunities in life cycle assessment[END_REF]. Land conversion from primary ecosystems to urban, industrial or agricultural land is a major biodiversity threat that has been studied for some years. The importance of subsequent environmental fragmentation, for which the situation is getting worse with respect to the corresponding Aichi biodiversity target (Convention on Biological Diversity, 2016), has been recognised more recently. However, to date most fragmentation studies have focused on the local scale whereas there are few analyses at regional scales or they are limited to the temperate zone or the tropics [START_REF] Bregman | Global patterns and predictors of bird species responses to forest fragmentation : Implications for ecosystem function and conservation[END_REF]. Except for a pioneering work that identified global fragmentation hotspots for mammalian carnivores [START_REF] Crooks | Global patterns of fragmentation and connectivity of mammalian carnivore habitat[END_REF], global patterns of habitat fragmentation and connectivity have not yet been examined, even if long-term experiments indicated that fragmentation effects are considerable and clearly consistent across a diverse range of terrestrial systems on five continents [START_REF] Haddad | Habitat fragmentation and its lasting impact on Earth's ecosystems[END_REF]. A study conducted at the biogeographical scale also showed that incorporating fragmentation metrics into largescale models may contribute for a better understanding of species distributions (Reino et al., 2013). With this in mind, global maps characterising the impacts of habitat fragmentation on species with relevant ecological indicators would be very useful for taking the appropriate decisions and actions for biodiversity conservation [START_REF] Bon | Global Biodiversity Change Indicators -Version 1[END_REF].

Fragmentation is commonly defined as a landscape-level process in which 'a large expanse of habitat is transformed into a number of smaller patches of smaller total area, isolated from each other by a matrix of habitats unlike the original' [START_REF] Wilcove | Habitat fragmentation in the temperate zone[END_REF]. Landscape structuration analysis is commonly conducted in parcels ranging in size from around 3 km 2 to 300 km 2 Author-produced version of the article published in Ecological Indicators, 2018, N°89, p.543-551. The original publication is available at https://www.sciencedirect.com Doi: 10.1016/j.ecolind.2018.01.068 [START_REF] Fischer | Landscape modification and habitat fragmentation : a synthesis[END_REF], and often set to 100 km 2 for different species in various studies [START_REF] Atauri | The role of landscape structure in species richness distribution of birds, amphibians, reptiles and lepidopterans in Mediterranean landscapes[END_REF][START_REF] Radford | Landscape-level thresholds of habitat cover for woodland-dependent birds[END_REF]. Fragmentation typically leads to the reduction of landscape connectivity, which is defined as the degree to which the landscape facilitates or hinders organisms movements among habitat patches [START_REF] Taylor | Connectivity is a vital element of landscape structure[END_REF].

Management to maintain or restore connectivity is crucial to ensure the survival of many species and preserve biological diversity. Therefore, many authors proposed metrics to measure habitat connectivity at the landscape scale [START_REF] Kindlmann | Connectivity measures : a review[END_REF][START_REF] Ortega | Indicator assessment for habitat fragmentation[END_REF][START_REF] Wang | Measuring habitat fragmentation: An evaluation of landscape pattern metrics[END_REF]. [START_REF] Calabrese | A comparison-shopper's guide to connectivity metrics[END_REF] proposed to distinguish between three main classes of connectivity metrics. First, the structural connectivity can be directly derived from landscape physical attributes (e.g. composition, shape or configuration metrics [START_REF] Rutledge | Landscape indices as measures of the effects of fragmentation: can pattern reflect process?[END_REF]). These types of indicators range from 'simple' landscape structure metrics (e.g. number of fragments) to more complex landscape indices (e.g. effective mesh size [START_REF] Jaeger | Landscape division, splitting index, and effective mesh size: new measures of landscape fragmentation[END_REF]), and dedicated software such as Fragstats [START_REF] Mcgarigal | FRAGSTATS: spatial pattern analysis program for quantifying landscape structure[END_REF] or, more recently, Conefor [START_REF] Saura | Environmental Modelling & Software Conefor Sensinode 2.2: A software package for quantifying the importance of habitat patches for landscape connectivity[END_REF] have been developed to easily compute them on different spatial levels. Data has also been collected to quantify structural fragmentation globally [START_REF] Haddad | Habitat fragmentation and its lasting impact on Earth's ecosystems[END_REF][START_REF] Riitters | Global-Scale Patterns of Forest Fragmentation[END_REF]. Second, the potential connectivity combines the landscape physical attributes, obtained from the structural connectivity analysis, with information on the focal species dispersal ability. Finally, the actual connectivity is based on observed or quantified (e.g. through species distribution models) movement pathways in the landscape. Structural connectivity is the easiest to collect, but it has a limited interest since it does not consider the species characteristics. The structural metrics ecological relevance (i.e. their relationship with the actual ecological processes taking place in the landscape) is often unproven and questionable [START_REF] Kupfer | Landscape ecology and biogeography : Rethinking landscape metrics in a post-FRAGSTATS landscape[END_REF]. Conversely, actual connectivity gives very detailed information on the landscape but requires a large amount of Author-produced version of the article published in Ecological Indicators, 2018, N°89, p.543-551. The original publication is available at https://www. sciencedirect.com Doi: 10.1016sciencedirect.com Doi: 10. /j.ecolind.2018.01.068 .01.068 data, and consequently hinders its implementation at larger scales. Thus, potential connectivity metrics seem to be a sound compromise for providing reasonable estimates of functional connectivity, i.e. species behavioural responses to landscape patterns [START_REF] Kindlmann | Connectivity measures : a review[END_REF], while being applicable at large-scale due to low data requirements.

One promising landscape potential connectivity metric at large-scale is the metapopulation capacity λ derived from metapopulation theory [START_REF] Hanski | The metapopulation capacity of a fragmented landscape[END_REF]. Hanski and Ovaskainen defined the metapopulation as a group of spatially separated populations of the same species interconnected by dispersal. In particular, the metapopulation capacity λ measures how a given spatial configuration of a set of fragments contributes to the long-term persistence of a particular species structured as metapopulation. This indicator has been widely used to estimate the fragmentation effects in various landscapes, e.g. agricultural landscapes [START_REF] Hietala-Koivu | Value of semi-natural areas as biodiversity indicators in agricultural landscapes[END_REF], forest landscapes [START_REF] Pardini | Beyond the Fragmentation Threshold Hypothesis : Regime Shifts in Biodiversity Across Fragmented Landscapes[END_REF]Schnell et al., 2013a), or fluvial landscapes [START_REF] Bertuzzo | Water Resources Research[END_REF], and it has also been tested and explored further on simulated landscapes [START_REF] Grilli | Metapopulation Persistence in Random Fragmented Landscapes[END_REF][START_REF] Rybicki | Species-area relationships and extinctions caused by habitat loss and fragmentation[END_REF]. Although the metapopulation capacity λ may be a promising indicator of landscape fragmentation, for now its use is limited to specific local landscapes. The aim of this paper is to develop an easily appropriable and parsimonious methodology to derive worldwide regionalised fragmentation indexes based on λ and to demonstrate the applicability and interest of this methodology with a dispersal distance of 1 km fitting a wide variety of species on all forest ecoregions in the biodiversity hotspots [START_REF] Myers | Biodiversity hotspots for conservation priorities[END_REF][START_REF] Mittermeier | Hotspots Revisited: Earth's Biologically Richest and Most Endangered Terrestrial Ecoregions[END_REF] as they are the richest and the most threatened plant and animal life reservoirs on Earth (Conservation International, 2016). This preliminary assessment will provide a global overview of the current abilities of birds or mammals with dispersal distances around 1 km to persist in the hotspots, which can be very useful in

Author-produced version of the article published in Ecological Indicators, 2018, N°89, p.543-551. The original publication is available at https://www. sciencedirect.com Doi: 10.1016sciencedirect.com Doi: 10. /j.ecolind.2018.01.068 .01.068 environmental assessment studies such as ones conducted in LCA or in footprint analyses. The innovative methodology itself also opens up vast prospects for biological conservation: with specific parameters adapted to a particular species in a specific region, its use could be extended to land use planning and areas prioritisation for both ecological protection and restoration.

Materials and methods

We used spatially explicit metapopulation models [START_REF] Hanski | Metapopulation dynamics[END_REF][START_REF] Hanski | The metapopulation capacity of a fragmented landscape[END_REF][START_REF] Ovaskainen | Spatially Structured Metapopulation Models : Global and Local Assessment of Metapopulation[END_REF] to describe habitat fragmentation at regional scales (e.g. ecoregion scale), since the metapopulation capacity λ can conveniently be used to rank different fragmented landscapes in terms of their capacity to support viable metapopulations. This section will present 1) the conventional metapopulation capacity calculation at the landscape scale and then 2) a methodology to extend it to large regions.

Metapopulation capacity of a species in a single landscape

Formally, in metapopulation theory, metapopulation capacity λ is given by: Leading eigenvalue λ of matrix M with elements , 0,

Where A i and A j are respectively the areas of fragments i and j, d ij is the Euclidian distance between the centroids of fragments i and j, and f(d ij ) is a dispersal function describing how the species arrival rate drops off with the distance d ij between fragments. The exponents x and y are scaling factors for the fragment area that depend on the extinction rate, the immigration rate and the emigration rate. The function f(d ij ) commonly presented is a negative exponential e -αd ij , in which α is the inverse of the species average dispersal distance. We adopted the following exponential dispersal kernel with a cut-off at 0.01 (equation 2) proposed by [START_REF] Hanski | Species-fragmented area relationship[END_REF]:

Author

Computing statistics on λ self calculated at the grid square scale for a large region

For the environmental assessment of human activities, baseline indicators on regional habitat fragmentation would be useful for decision-making along globalised supply chains (i.e. along product or service life cycles). However, the metapopulation capacity should be quantified at a scale that is relevant for the focal species. According to [START_REF] Olson | Terrestrial Ecoregions of the World : A New Map of Life on Earth[END_REF], the ecoregion average area is about 150,000 km 2 (median ≈ 60,000 km 2 ) while the biodiversity hotspots have a mean surface area around 800,000 km 2 (median ≈ 300,000 km 2 ). Given that the sizes of these biogeographic units are both much larger than the model species dispersal abilities, it would be ecologically meaningless to directly calculate the metapopulation capacity at such spatial scales. Thus, to compute λ self for larger areas, we divided the entire region of interest by a virtual grid and calculated λ self,GS for selected Grid Squares (GS) (see Figure A1 in Appendix A). The landscape structure was characterised by a nonhabitat matrix surrounding native habitat fragments. Using Q-GIS (QGIS Development Team, 2017), the habitat fragment centroids were identified, and the data necessary to calculate λ self,GS were extracted, i.e. each fragment size A i , and its centroid coordinates X i , Y i . The distances d ij between the fragment centroids within each grid square were calculated, then, by introducing the dispersal distances of the modelled species, λ self,GS was calculated according to equation 1 (see Appendix A, Section 1-E for programming details).

Deriving fragmentation indexes for large regions based on λ self,GS statistics

Statistical analysis of λ self,GS values calculated at the grid square scale can be used to obtain fragmentation indexes for a whole region.

Several studies support the hypothesis that fragmentation effects follow a threshold function related to the remaining available habitat. When habitat area is still large enough, fragmentation can have a positive impact, i.e. by increasing the local species richness, but which will be hidden by the superior, negative habitat loss effects. When the habitat area falls below a critical threshold, generally around 20-30% [START_REF] Andren | Effects of habitat fragmentation on birds and mammals of suitable habitat : a review landscapes with different proportions[END_REF][START_REF] Flather | Patchy reaction-diffusion and population abundance: the relative importance of habitat amount and arrangement[END_REF][START_REF] Hanski | Habitat fragmentation and species richness[END_REF][START_REF] Pardini | Beyond the Fragmentation Threshold Hypothesis : Regime Shifts in Biodiversity Across Fragmented Landscapes[END_REF][START_REF] Schneider | Habitat loss , fragmentation and predator impact : spatial implications for prey conservation[END_REF], the fragmentation effects become negative, affecting species survival, and should not be ignored. Consequently, λ self,GS were exclusively derived for grid squares with less than 30% of their original vegetation (defined as Conditions of λ self validity in Figure 1Figure 1). Original or native vegetation can refer to various ecosystems (among others, grasslands, savannas, forests, or deserts). In this study, only primary forest ecosystems were considered and λ self,GS was not calculated for grid squares with no primary forest at all. Based on λ self,GS values of all grid squares in which forest cover ranges from >0% to ≤30%thereafter referred to as 'highly converted forest'-, we can determine a median value λ self,j of the non-normally distributed λ self,GS (Figure A2 in Appendix A), that stands for the metapopulation capacity of the whole large region j most converted forest areas (grey box in Figure 1Figure 1). Once λ self,j has been computed for each large region j (e.g. ecoregions), the median values of all large regions are ranked from the 2.5 th percentile (P2.5) to the 97.5 th percentile (P97.5).

Author Subsequent winsorising reduces the effect of possible outliers by limiting extreme values in the data via setting outliers to a data specified percentile. Here, a 95% winsorising was applied, meaning that all data below the 2.5 th percentile are set equal to the 2.5 th percentile (=λ self,min ) and all data above the 97.5 th percentile are set equal to the 97.5 th percentile (=λ self,max ). A 95%

winsorising was found to deal adequately with extreme values without being too exclusive on the dataset. Normalisation from 0 (low fragmentation stress) to 1 (high fragmentation stress)

provides a commensurable and dimensionless index, i.e. the Forest Fragmentation Potential of highly converted forests FFP Forest≤30%,j of a large region j according to equation 3. In this way, the 5% extreme values of the global dataset, i.e. the least fragmented regions (λ self ≥ P97.5) and the most fragmented regions (λ self ≤ P2.5), have a Forest Fragmentation Potential set to 0 and 1, respectively. In between, the fragmented regions are classified according to their fragmentation degree.

FFP Forest≤30%,j indicates the fragmentation 'state' of the most critical forest areas within region j. The more region j is fragmented, the smaller is λ self,j and the higher is FFP Forest≤30%,j :

%, , , , , (3)
To characterise the Forest Fragmentation Potential of the whole region FFP Forest,j , we weight FFP Forest≤30%,,j by the proportion between the highly converted forest grid squares and the total number of forest grid squares in region j 30% FOREST , :

, %, 30% FOREST , (4) With 30% FOREST , ∑ , % % ∑ , % %
Finally, to determine the whole region fragmentation state, we use the regional forest area proportion, defined as the ratio of forest grid squares and total number of grid squares (i.e.

With , ∑ GS,with 0% 100% ∑ GS

Implementing the methodology on forest ecoregions included in the biodiversity hotspots

For a global implementation, we first focus our assessment on the biodiversity hotspots. To qualify as a biodiversity hotspot, a region must meet two criteria, i.e. 1) holding high numbers of endemic vascular plants, and 2) having 30% or less of its original vegetation (Conservation International, 2016), thus corresponding to the 'fragmentation threshold' mentioned above, below which considering fragmentation is relevant. Furthermore, [START_REF] Sloan | Remaining natural vegetation in the global biodiversity hotspots[END_REF] recently provided updated estimates of natural, intact vegetation within these hotspots that show that their state could be even more critical than previously described. We considered all ecoregions included in the biodiversity hotspots with a focus on forest major habitat types (see Appendix A, Section 1-B for more details). Area loss in forests is rarely simply the removal of contiguous areas and these habitats are particularly experiencing fragmentation [START_REF] Whitmore | Tropical Deforestation and Species Extinction[END_REF], mainly because of industrial timber extraction, agricultural expansion, fire, and resource extraction [START_REF] Potapov | The last frontiers of wilderness : Tracking loss of intact forest landscapes from 2000 to 2013[END_REF]. Due to fragmentation, the number of intact forest landscapes (larger than 500 km 2 ) is decreasing drastically, around 7% in 13 years, particularly in tropical regions [START_REF] Potapov | The last frontiers of wilderness : Tracking loss of intact forest landscapes from 2000 to 2013[END_REF]. Thus, forests are under pressure like never before, while they host 90% of terrestrial species (WWF, 2016b).

Author-produced version of the article published in Ecological Indicators, 2018, N°89, p.543-551. The original publication is available at https://www.sciencedirect.com Doi: 10.1016/j.ecolind.2018.01.068

We applied the methodology developed (Figure 1Figure 1) to all biodiversity hotpots that contain forest ecoregions (i.e. 34 forest hotspots: six in Africa, fourteen in Asia-Pacific, four in Europe and Central Asia, five in North and Central America, and five in South America). 259 ecoregions met the conditions for λ self,j calculation (see Conditions of λ self validity in Figure 1Figure 1), and FFP j were provided for a total of 283 ecoregions, because fragmentation indexes of the 24 ecoregions for which all whole grid squares had a forest cover higher than 30% were set to 0 (see right side of Figure 1Figure 1). Details on spatial processing can be found in Appendix A, Section 1-D.

Input data for model parameters

Two kinds of data are required to implement the developed methodology, i.e. 1) species characteristics (dispersal ability and scaling factors), and 2) landscape inputs (original habitat choice, map resolution, and grid size for spatial analysis). For each of them, data constraints and availability are discussed to select the most appropriate data for the λ self global calculation (Table A1 in Appendix A).

Species data are essential for the metapopulation capacity model. Provided that 'forest' is its habitat, a species can be heavily or slightly affected depending on its ability to move (1/α, in connection with d ij ), and on its immigration and extinction abilities (x and y, related to A i ).

Various 'dispersal distances' are reported in literature, e.g. the median, mean and maximum dispersal distances, and the migratory distance. The differences between these variables are hard to capture because their calculation varies depending on the author. Furthermore, there is a large variability associated with species dispersal distances, depending on their size or migratory habits. Based on a literature review of available species dispersal distances [START_REF] Sutherland | Scaling of Natal Dispersal Distances in Terrestrial Birds and Mammals[END_REF], we arbitrarily considered a 1 km dispersal distance, which can apply to some birds and

Author-produced version of the article published in Ecological Indicators, 2018, N°89, p.543-551. The original publication is available at https://www.sciencedirect.com Doi: 10.1016/j.ecolind.2018.01.068 mammals, as a first example of the methodology application. Regarding scaling factors x and y, [START_REF] Hanski | Species-fragmented area relationship[END_REF] give a x realistic value for birds and mammals (1.5), while Schnell et al. (2013b) propose another value valid for many taxa (0.5). Both studies considered a y value of 1.

We used the values provided by Schnell et al. (2013b) because they offer a broader choice in terms of applications to different taxa.

The fragmentation assessment strongly depends on the land cover of interest definition [START_REF] Riitters | Global-Scale Patterns of Forest Fragmentation[END_REF]. In this work, we used Globcover 2009 maps [START_REF] Bontemps | GLOBCOVER 2009: Products Description and Validation Report[END_REF], produced by the European Space Agency, which count 22 land cover classes. Among them, several classes can be included in the forest definition, and considering only class 50 Closed (>40%) broadleaved deciduous forest (>5m), both classes 50 and 40 Closed to open (>15%) broadleaved evergreen or semi-deciduous forest ( >5m), or all eight potential forest classes combined, can have a large influence on the results, underestimating or overestimating the real forest cover. Furthermore, reports usually describe the remaining original habitat area in percentage 'of its original extent' (CEPF, 2016;[START_REF] Di Bitetti | A biodiversity vision for the Upper Parana Atlantic Forest ecoregion: designing a biodiversity conservation landscape and setting priorities for conservation action[END_REF], but all the ecoregion was not necessarily covered with this original habitat (other land covers, water bodies, mountains, etc.). To simplify our methodology, we considered all Globcover 2009 forest classes.

Map resolution is also a crucial parameter, since the notion of 'fragment' strongly depends on the scale considered. Many metrics indicated lower fragmentation at coarser spatial resolutions [START_REF] Saura | Effects of remote sensor spatial resolution and data aggregation on selected fragmentation indices[END_REF], and a recent study showed that metapopulation capacity is sensitive to the spatial scale at which the habitat units were defined (Blazquez-Cabrera, 2014). Maps with smaller resolutions identify more non-forest area where forest cover is dominant but not exclusive. [START_REF] Hanski | Species-fragmented area relationship[END_REF] used a map resolution of 30 m, but they used a bird species dispersal distance much smaller (300 m) than the median value found in literature, i.e. 4 km for birds

Author-produced version of the article published in Ecological Indicators, 2018, N°89, p.543-551. The original publication is available at https://www.sciencedirect.com Doi: 10.1016/j.ecolind.2018.01.068 [START_REF] Sutherland | Scaling of Natal Dispersal Distances in Terrestrial Birds and Mammals[END_REF]. A transposition based on ecological parameters would lead to a scaleup of a factor ten for map resolution (400 m). Thus, we took Globcover 2009 (300 m resolution) for land cover data.

Likewise, fragmentation can turn out greater when larger grid sizes are used [START_REF] Riitters | Global-Scale Patterns of Forest Fragmentation[END_REF]. Consequently, to perform a meaningful metapopulation capacity calculation, it is essential to define spatial scale(s) at which metapopulation processes predominantly take place. [START_REF] Olson | Terrestrial Ecoregions of the World : A New Map of Life on Earth[END_REF] examined a moth metapopulation with a measured dispersal distance of 100 m at four levels of spatial analysis and concluded that metapopulation processes were important at the regional and sub-regional scales (i.e. ≈35 km 2 and ≈10 km 2 , with 100 m and 50 m resolutions respectively). Consequently, for dispersal distances ranging from 500 m to 5 km, spatial scales between 50 and 2000 km 2 would be consistent, with associated resolutions ranging from 250 m to 5 km. Following the same logic as for map resolution, we could also have scaled-up the grid size for λ self calculations based on [START_REF] Hanski | Species-fragmented area relationship[END_REF], but we found 100 km 2 more adapted to the range of ecoregion sizes. 1000 km 2 grid squares were too large to adequately reflect some small ecoregions, e.g. for the Montane Fynbos And Renosterveld ecoregion (AT1203), only three grid squares were larger than 950 km 2 , representing only 6% of the whole ecoregion area. Thus, a 100 km 2 -grid was applied to each hotspot, and grid squares larger than 95 km 2 (5% margin) were selected for each ecoregion. Consequently, all ecoregions smaller than the grid square size (100 km 2 ) were excluded from the outset, representing less than 0. 2% in surface of all forest ecoregions belonging to the biodiversity hotspots.

Results

Spatial representation of the three FFPs

Author-produced version of the article published in Ecological Indicators, 2018, N°89, p.543-551. The original publication is available at https://www.sciencedirect.com Doi: 10.1016/j.ecolind.2018.01.068

The FFPs of all 30 ecoregions included in the Indo-Burma hotspot are spatially represented for each 100-km 2 grid square included in the hotspot (Figure 2Figure 4). FFP Forest≤30%,j is a median value of the most critical forest areas fragmentation 'state' within the ecoregion j, thus only grid squares in which forest cover is smaller than 30% are concerned (upper-map in Figure 2Figure 4). FFP Forest≤30% can be used to assess the potential impact on biodiversity of a current or future activity using land in a highly converted forest area. As FFP Forest weights the previous index by the total forest cover in the ecoregion, the index value is affected to all forest grid squares (middle-map in Figure 2Figure 4). FFP Forest can be used to compare activities for which land use surely occurs in a forest area, with no additional information on the forest cover. Lastly, FFP All applies to all grid squares in the ecoregion (lower-map in Figure 2Figure 4) and can be used to rank land use impacts in different ecoregions when the exact locations of activities in the ecoregions are unknown. The more spatially aggregated the fragmentation indexes (from FFP Forest≤30% , then FFP Forest , to finally FFP All ), the more uncertain they are.

More explanations on the different indexes are given with the example of the Indo-Burman ecoregions A, B and C (Figure 2Figure 4). All three ecoregions have a high FFP Forest≤30% (respectively 0.92, 0.82 and 0.99), which means that their most converted forest zones are highly fragmented and have a strong negative impact on species (compared to other ecoregions worldwide

Global maps of FFPs

The fragmentation indexes were calculated for each forest ecoregion included in the 34 biodiversity hotspots for an illustrative species dispersal distance of 1 km. Global map of FFP All is presented in Figure 3Figure The most converted forests within the ecoregions (i.e. 0% < forest cover ≤ 30%) globally show a high fragmentation potential (Figure A3 in Appendix A). Only 42 out of the 259 ecoregions for which λ self was calculated have a fragmentation potential FFP Forest≤30% smaller than 0.6, and 57 ecoregions (mainly in the Palearctic and Indomalayan realms) have it higher than 0.95. However, FFP Forest≤30% reflects the fragmentation potential of already altered forest areas, and does not reflect the forest 'state' of the whole ecoregion.

FFP Forest , i.e. FFP Forest≤30% weighted by the proportion of highly converted forest zones compared to all forest zones, provides more distributed results globally (Figure A4 in Appendix A). Only 22 ecoregions, mainly in Indo-Burma, Mediterranean Basin and Irano-Anatolian hotspots, present FFP Forest higher than 0.90. In addition, there are strong differences between FFP Forest≤30% and FFP Forest results (respective median values of 0.86 and 0.27), which indicate that fragmentation is most frequently a process confined to relatively specific geographical areas.

For instance, some ecoregions have highly fragmented forests (high FFP Forest≤30% ), but only over

Author-produced version of the article published in Ecological Indicators, 2018, N°89, p.543-551. The original publication is available at https://www.sciencedirect.com Doi: 10.1016/j.ecolind.2018.01.068 a small part of their total forest area (low FFP Forest ). The spatial aggregation at the ecoregion scale tends to smooth the results.

FFP All , weighted in proportion to ecoregion forest areas, has values very similar to FFP Forest (global median of 0.23). This is because there are few landscapes (100 km 2 grid squares) with no forest at all (i.e. the forest share FOREST TOTAL often approaches 1), which is expected for ecoregions with forest as major habitat types. Nevertheless, for some ecoregions, e.g. for the four North African ecoregions included in the Mediterranean Basin hotspot, FFP Forest can be high and FFP All much lower, due to their small forest cover. Conversely, e.g. in Indonesia or South America, both FFP Forest and FFP All can be high because of some ecoregions large forest share.

Only three ecoregions have FFP All higher than 0.9, and they all are in the Neotropical realm (ecoregions NT0233 Veracruz Dry Forests, NT0102 Atlantic Coast Restingas and NT0151

Pernambuco Coastal Forests).

Author-produced version of the article published in Ecological Indicators, 2018, N°89, p.543-551. The original publication is available at https://www.sciencedirect.com Doi: 10.1016/j.ecolind.2018.01.068

Discussion

An easily appropriable and performing methodology to derive worldwide regionalised λ-based fragmentation indexes at meaningful scales for fragmentation mechanisms was built, and applied to an illustrative, arbitrary example (by setting the species dispersal distance to 1 km, and using scaling factors adapted to birds or mammals). The objective was twofold: (1) to test the feasibility of the method at a worldwide scale, and ( 2) to evaluate the resulting indexes' ability to provide valuable information on fragmentation stresses. The methodology was successfully applied in about 300 ecoregions all over the world, and the resulting indexes clearly confirmed significant intra-and inter-ecoregions differences regarding landscapes' fragmentation stresses.

These results open up prospects for promising applications of the methodology to real ecological studies. However, the proposed methodology still faces some methodological and practical limitations and still has room for improvement.

Metapopulation model limitations

The interest of metapopulation theory to explain variations in species occurrence patterns in highly fragmented landscapes has been recognised. However, metapopulation theory implies that the populations exist in patchy distributions. When there is no evidence for metapopulation dynamics in a given species or region, this approach might not be valid. As a great range of dispersal models are available in the literature, the more adapted and realistic dispersal functions f(d ij ) should be incorporated in spatial models when possible [START_REF] Travis | Dispersal functions and spatial models: expanding our dispersal toolbox[END_REF]. In the proposed methodology, due to the metapopulation model flexibility, other more complex dispersal functions could easily replace the proposed one if needed. Furthermore, because

Author-produced version of the article published in Ecological Indicators, 2018, N°89, p.543-551. The original publication is available at https://www.sciencedirect.com Doi: 10.1016/j.ecolind.2018.01.068 species dispersal still needs to be better understood and assessed, new ecological developments and data could also be easily incorporated in the model as soon as they become available.

Modelling metapopulation dynamics is a very useful tool for conservation ecology. However, metapopulation models have historically treated a landscape as a collection of habitat patchesin our case, forest-separated by a homogenous, permeable but unsuitable habitat usually called matrix-in our case, non-forest-. In these models, movement between patches depends only on the distance between patches and the inherent species dispersal ability, and the matrix structure and quality effects on movement through the landscape are neglected. Yet, even if theoretical and empirical evidence show that matrix quality can be extremely important in determining metapopulation dynamics [START_REF] Vandermeer | Metapopulation Dynamics and the Quality of the Matrix[END_REF], these factors are rarely incorporated into metapopulation models. In particular, linear infrastructure such as road networks can considerably hinder species movements (often referred to as the 'barrier effect' in ecology [START_REF] Forman | Roads and Their Major Ecological Effects[END_REF]) and have a strong negative impact on species [START_REF] Loro | Ecological connectivity analysis to reduce the barrier effect of roads. An innovative graph-theory approach to define wildlife corridors with multiple paths and without bottlenecks[END_REF].

Moreover, this barrier effect varies among species, which are not similarly affected by a given barrier. [START_REF] Gebauer | Matrix Matters : Differences of Grand Skink Metapopulation Parameters in Native Tussock Grasslands and Exotic Pasture Grasslands[END_REF] investigated the importance of matrix type for metapopulation modelling and found that it had as much or sometimes more influence than patch sizes or distances between patches on metapopulation parameters. The use of mathematical tools enabled by the matrix-modelling framework, such as sensitivity and elasticity analyses [START_REF] Shima | Life history and matrix heterogeneity interact to shape metapopulation connectivity in spatially structured environments[END_REF], or the incorporation of a matrix resistance parameter, which is expected to vary among species [START_REF] Ricketts | The Matrix Matters: Effective Isolation in Fragmented Landscapes[END_REF], could definitively improve metapopulation modelling. If matrix quality was to be included in the metapopulation model, complementary road maps could be used in addition to Globcover 2009 maps to capture the barrier effects induced by linear infrastructure.

Author-produced version of the article published in Ecological Indicators, 2018, N°89, p.543-551. The original publication is available at https://www. sciencedirect.com Doi: 10.1016sciencedirect.com Doi: 10. /j.ecolind.2018.01.068 .01.068 Finally, the metapopulation capacity assesses species persistence in fragmented landscapes, which is a limited aspect of all possible ecological attributes and levels of biodiversity organisation, which encompasses not only species composition but also ecosystem structure and functions.

Data availability and uncertainty

There is inherent uncertainty in data (e.g. land cover maps, species dispersal distances), on choices associated with data treatment, and on the metapopulation model itself, i.e. to what extent does it reflect real species behaviour. The results should hence be interpreted in light of this information [START_REF] Langford | Map Misclassification Can Cause Large Errors in Landscape Pattern Indices : Examples from Habitat Fragmentation[END_REF]. Globcover 2009 maps are among the most detailed, reliable and up-to-date global land cover maps. Their overall accuracy weighted by the class area reaches 67.5%. The Globcover 2009 resolution (300 m) is consistent with the spatial scale considered (i.e. global statistics on 100 km 2 -landscapes), and execution times for the associated spatial processing are acceptable. Regarding species data, the selected dispersal distances were chosen to be of the same order of magnitude than estimates derived from empirical data from a global literature review [START_REF] Sutherland | Scaling of Natal Dispersal Distances in Terrestrial Birds and Mammals[END_REF]. The uncertainty associated with these estimates is very closely linked to the natural data variability. A way to overcome this shortcoming would be to generate more results to cover a more complete range of species dispersal distances. As more biodiversity monitoring data will become available, these dispersal distances estimates should be updated and results accuracy will be improved. Moreover, the indexes could be combined with data on specific species effectively inhabiting each ecoregion to enhance the results relevance. Another interesting development would be to weight the indexes by the threatened status and/or the level of species endemism, e.g. percent endemism by taxonomic group in each hotspot (CEPF, 2016).

Author-produced version of the article published in Ecological Indicators, 2018, N°89, p.543-551. The original publication is available at https://www. sciencedirect.com Doi: 10.1016sciencedirect.com Doi: 10. /j.ecolind.2018.01.068 .01.068 The λ self sensitivity to several crucial parameters related to 1) land cover data (original habitat and maps resolution), 2) grid size, and 3) species data (dispersal distance and scaling factor x)

was assessed and the results can be found in Appendix A, Section 3-A. Other data or methodological choices likely to affect the results were also discussed in the same section.

However, even if some parameters may have a great influence on λ self , no change in these will affect their ranking (i.e. for ecoregions ranking). The λ self final distribution and subsequent FFP indexes derived from λ self , might be affected though, given that the λ self dependency on the input parameters and choices is not necessarily linear.

Regarding the fragmentation index calculation based on λ self statistics, the λ self,GS values aggregation for large regions (e.g. ecoregions), by means of the median, is another important source of uncertainty, due to the intra-ecoregion spatial variability. As an example, the λ self,GS spatial variability of the larger forest ecoregion included in the biodiversity hotspots, i.e. Alto Paraná Atlantic Forest (ecoregion NT0150), is illustrated in Figure A5 in Appendix A. The results show that, rather than being distributed evenly throughout the ecoregion, fragmentation can be very critical in some specific places of the ecoregion. Similarly to other stress indexes (e.g. water stress), a global or regional fragmentation stress can strongly differ from a local stress, and too aggregated stress values will be unrepresentative of the local reality.

Likewise, the definition and spatial resolution of Forest Shares FS used for FFP calculation may influence the results and is discussed in Appendix A, Section 3-C.

Applicability

The fragmentation indexes can be used directly to rank all highly fragmented forest ecoregions included in the biodiversity hotspots. As they are based on metapopulation capacity, which is already a relative measure, and given their strong sensitivity to the input parameters, the Author-produced version of the article published in Ecological Indicators, 2018, N°89, p.543-551. The original publication is available at https://www.sciencedirect.com Doi: 10.1016/j.ecolind.2018.01.068 fragmentation indexes should not be taken as absolute values and should be used for comparative purposes only. Assuming that the species dispersal ability is known and similar in magnitude to the considered dispersal distance, conclusions can be drawn for this particular species. However, for now, they are limited to a specific range of species dispersal distances, based on animal dispersal abilities statistics. Plants should also be considered, especially as metapopulation concept has been largely applied to plants through seed dispersal [START_REF] Verheyen | Metapopulation Dynamics in Changing Landscapes : A New Spatially Realistic Model for Forest Plants[END_REF].

To overcome this problem, the proposed grid procedure and methodology for deriving metapopulation capacity statistics at regional scales was designed to be easily used and applied to other large datasets. Our analysis was restricted to biodiversity hotspots as defined by [START_REF] Myers | Biodiversity hotspots for conservation priorities[END_REF], but other prioritisation approaches exist for global biodiversity conservation, and other templates may be considered [START_REF] Brooks | Global Biodiversity Conservation Priorities[END_REF][START_REF] Moran | Identifying species threat hotspots from global supply chains[END_REF]. We also focused the scope of application on the ecoregions whose major habitat type was forest (i.e. approximately 70% in surface of all ecoregions belonging to the biodiversity hotspots), but other major habitat types, e.g. grasslands, could easily be incorporated. Conversely, the proposed methodology could also be applied to specific ecoregions or sub-ecoregions, with a grid size and a map resolution adapted to a particular species. To widen the range of model applications while ensuring reproducibility, and to ease new fragmentation indexes calculations, the matrix calculation program is available upon request.

Conclusions and perspectives

The global forest fragmentation is a critical issue facing major challenges for the biodiversity conservation and ecosystems management, particularly in tropical environments [START_REF] Bregman | Global patterns and predictors of bird species responses to forest fragmentation : Implications for ecosystem function and conservation[END_REF], and especially since these effects may be worsened by other global changes caused by human activity, e.g. climate change [START_REF] Haddad | Habitat fragmentation and its lasting impact on Earth's ecosystems[END_REF]. Macro-scale conservation planning

Author-produced version of the article published in Ecological Indicators, 2018, N°89, p.543-551. The original publication is available at https://www. sciencedirect.com Doi: 10.1016sciencedirect.com Doi: 10. /j.ecolind.2018.01.068 .01.068 gives a means of coping with the challenge arising from the environmental impacts of extensive land use changes for urban or agricultural activities [START_REF] Olsoy | 068 Quantifying the effects of deforestation and fragmentation on a range-wide conservation plan for jaguars[END_REF]. To this end, global information on areas sensitive to fragmentation is needed to improve the effectiveness of planning efforts.

Thus, an innovative methodology for modelling the fragmentation potential of large regions was proposed. The landscape spatial configuration is combined with a metapopulation model to obtain an ecologically relevant landscape fragmentation metric to compare any fragmented landscapes, producing an efficient tool for land management. The modelling approach applied in this work was developed to fit most species, provided that some adaptations are made to deal with species specificities. Indeed, the flexibility of the methodology allows the spatial parameters, e.g. the grid size and the map resolution, and the species parameters, e.g. the dispersal distance and the dispersal function of the metapopulation model, to be easily adapted to a specific species. Secondly, as an example, we provide a first set of fragmentation indexes for ranking all forest ecoregions (> 100 km 2 ) belonging to the biodiversity hotspots. These global maps are relevant to the subset of species corresponding to the types of habitat and to the dispersal capacities considered, and provide first-hand information on the ability of these species, i.e. birds or mammals with a dispersal distance of 1 km, to persist in highly fragmented habitats. These indexes can help decision makers, from producers, scientists, conservationists, to governments, to better target their actions. They can be used as a tool for land planners or industrial companies, for comparing different options regarding their activities locations at the ecoregion scale.

Conversely, they enable focusing on areas where offset measures would have the most impact [START_REF] Dalang | Trading connectivity improvement for area loss in patchbased biodiversity reserve networks[END_REF]. Finally, fragmentation metrics at large-scale can be used to

Author-produced version of the article published in Ecological Indicators, 2018, N°89, p.543-551. The original publication is available at https://www. sciencedirect.com Doi: 10.1016sciencedirect.com Doi: 10. /j.ecolind.2018.01.068 .01.068 improve predictive modelling of range shifts associated with land use (Reino et al., 2013). More particularly, they could contribute to better assess the environmental impacts of a product using land at different places around the globe along its life cycle. To date, several indicators are available in LCA to quantify a certain number of land use impacts (e.g. the impact of decreasing their habitat quantity on species), but, until now, habitat fragmentation impacts are poorly considered. In this perspective, the fragmentation indicators should be calculated for each taxon usually considered for land use impacts in LCA, i.e. birds, mammals, reptiles, amphibians and vascular plants, using median values of species or seed dispersal distances and scaling factors for which estimates are available in the literature (see for example [START_REF] Kharouba | Do ecological differences between taxonomic groups influence the relationship between species' distributions and climate? A global meta-analysis using species distribution models[END_REF]; [START_REF] Sutherland | Scaling of Natal Dispersal Distances in Terrestrial Birds and Mammals[END_REF] and [START_REF] Verheyen | Metapopulation Dynamics in Changing Landscapes : A New Spatially Realistic Model for Forest Plants[END_REF]). They should also be extended to other habitat types (e.g. grasslands, savannas, shrublands, tundra, mangroves and deserts). In addition, as they stand, the fragmentation indexes are stand-alone indicators which are not directly linked to a land use intervention. Incorporating them into analytical environmental assessment methods such as LCA will need further development.

Figure 1 :

 1 Figure 1: Methodology for calculating the fragmentation indexes (Forest Fragmentation Potential FFP) of a large region at three levels of spatial aggregation (highly converted forest, forest and the entire region). P: Percentile. NC: Not calculated.

Figure 2 :

 2 Figure 2: Fragmentation indexes for forest ecoregions included in the Indo-Burma biodiversity hotspot calculated for a species dispersal distance of 1 km; arrow A points to Meghalaya Subtropical Forests (ecoregion IM0126), arrow B points to Chin Hills-Arakan Yoma Montane Forests (ecoregion IM0109), and arrow C points to Irrawaddy Dry Forests (ecoregion IM0205)

  Figure 3: Forest Fragmentation Potential FFP All for all forest ecoregions included in the biodiversity hotspots (coloured areas) calculated for a species dispersal distance of 1 km

  Author-produced version of the article published in Ecological Indicators, 2018, N°89, p.543-551. The original publication is available at https://www.sciencedirect.com Doi: 10.1016sciencedirect.com Doi: 10. /j.ecolind.2018.01.068 .01.068

  Author-produced version of the article published in Ecological Indicators, 2018, N°89, p.543-551. The original publication is available at https://www.sciencedirect.com Doi: 10.1016sciencedirect.com Doi: 10. /j.ecolind.2018.01.068 .01.068

  -produced version of the article published in Ecological Indicators, 2018, N°89, p.543-551. The original publication is available at https://www.sciencedirect.com Doi: 10.1016sciencedirect.com Doi: 10. /j.ecolind.2018.01.068 .01.068

  Author-produced version of the article published in Ecological Indicators, 2018, N°89, p.543-551. The original publication is available at https://www.sciencedirect.com Doi: 10.1016sciencedirect.com Doi: 10. /j.ecolind.2018.01.068 .01.068

	including forest-free grid squares) in the region ( FOREST TOTAL	, in equation 5). The Forest
	Fragmentation Potential FFP All, j expresses the forest fragmentation potential for the whole large
	region j (e.g. an ecoregion):
	,	,	FOREST TOTAL	,

  ). Overall forest, however, is less damaged in ecoregion B than in ecoregions A or C, i.e. grid squares with less than 30% of forest are scarcer in ecoregion B, so FFP Forest,B (0.10) is much lower than FFP Forest,A (0.56) or FFP Forest,C (0.86). However, almost all 100 km 2 grid squares contain forest in both ecoregions A and B whereas ecoregion C contains many forestfree grid squares, so FFP All,A (0.54) and FFP All,B (0.10) are close to FFP Forest,A and FFP Forest,B whereas FFP All,C (0.18) decreases considerably compared to FFP Forest,C . -produced version of the article published in Ecological Indicators, 2018, N°89, p.543-551. The original publication is available at https://www.sciencedirect.com Doi: 10.1016/j.ecolind.2018.01.068

	Author

Author-produced version of the article published in Ecological Indicators, 2018, N°89, p.543-551. The original publication is available at https://www.sciencedirect.com Doi: 10.1016sciencedirect.com Doi: 10. /j.ecolind.2018.01.068 .01.068

Acknowledgments

The authors thank Jean-Louis Martin and Ana Rodrigues from CEFE, Montpellier for their inspiring and valuable insights on the fragmentation issue, Thomas Koellner and Asja Bernd from the University of Bayreuth, Germany for their useful inputs on the fragmentation assessment, Jean-Baptiste Mihoub from Pierre and Marie Curie University, Paris for fruitful discussions about global biodiversity assessment, and Samuel Alleaume from UMR TETIS, Irstea Montpellier for his good advice on spatial processing. The authors are members of the ELSA research group (Environmental Life Cycle and Sustainability Assessment, http://www.elsa-lca.org/) and thank all ELSA members for their advice. This work was supported by the French National Research Agency (ANR), the Occitanie Region, ONEMA, Ecole des mines d'Alès and the industrial partners (BRL, SCP, SUEZ, VINADEIS, and Compagnie Fruitière) of the Industrial Chair for Environmental and Social Sustainability Assessment 'ELSA-PACT' (ANR grant no. 13-CHIN-0005-01).

Appendices

Appendix A: Additional methods and results (PDF)

Appendix B: Metapopulation capacities λ self , forest shares FS, and FFP indexes for forest ecoregions included in the biodiversity hotspots for a 1 km species dispersal distance (XLSX)

Author-produced version of the article published in Ecological Indicators, 2018, N°89, p.543-551. The original publication is available at https://www.sciencedirect.com Doi: 10.1016Doi: 10. /j.ecolind.2018.01.068 .01.068


Nikolaus Hansen

Steffen Finck

Raymond Ros

Anne Auger

Real-Parameter Black-Box Optimization Benchmarking 2009: Noiseless Functions Definitions

Keywords: optimization, evolutionary algorithms, benchmarking, black-box

Quantifying and comparing performance of optimization algorithms is one important aspect of research in search and optimization. However, this task turns out to be tedious and difficult to realize even in the single-objective case -at least if one is willing to accomplish it in a scientifically decent and rigorous way. The BBOB 2009 workshop will furnish most of this tedious task for its participants: (1) choice and implementation of a well-motivated real-parameter benchmark function testbed, (2) design of an experimental set-up, (3) generation of data output for (4) post-processing and presentation of the results in graphs and tables. What remains to be done for the participants is to allocate CPU-time, run their favorite black-box real-parameter optimizer in a few dimensions a few hundreds of times and execute the provided post-processing script afterwards. In this report, the testbed of noise-free functions is defined and motivated.

Introduction

In the following, 24 noise-free real-parameter single-objective benchmark functions are defined (for a graphical presentation see [START_REF] Finck | Real-parameter black-box optimization benchmarking 2009: Presentation of the noiseless functions[END_REF]). 1 Our intention behind the selection of benchmark functions was to evaluate the performance of algorithms with regard to typical difficulties which we believe occur in continuous domain search. We hope that the function collection reflects, at least to a certain extend and with a few exceptions, a more difficult portion of the problem distribution that will be seen in practice (easy functions are evidently of lesser interest). We prefer benchmark functions that are comprehensible such that algorithm behaviours can be understood in the topological context. In this way, a desired search behaviour can be pictured and deficiencies of algorithms can be profoundly analysed. Last but not least, this can eventually lead to a systematic improvement of algorithms.

All benchmark functions are scalable with the dimension. Most functions have no specific value of their optimal solution (they are randomly shifted in x-space). All functions have an artificially chosen optimal function value (they are randomly shifted in f -space). Consequently, for each function different instances can be generated: for each instance the randomly chosen values are drawn anew2 . Apart from the first subgroup, the benchmarks are non-separable. Other specific properties are discussed in the appendix.

General Setup

Search Space All functions are defined and can be evaluated over R D , while the actual search domain is given as [-5, 5] D .

Location of the optimal x opt and of f opt = f (x opt ) All functions have their global optimum in [-5, 5] D . The majority of functions has the global optimum in [-4, 4] D and for many of them x opt is drawn uniformly from this compact. The value for f opt is drawn from a Cauchy distributed random variable, with zero median and with roughly 50% of the values between -100 and 100. The value is rounded after two decimal places and set to ±1000 if its absolute value exceeds 1000. In the function definitions a transformed variable vector z is often used instead of the argument x. The vector z has its optimum in z opt = 0, if not stated otherwise.

Boundary Handling On some functions a penalty boundary handling is applied as given with f pen (see next section).

Linear Transformations Linear transformations of the search space are applied to derive nonseparable functions from separable ones and to control the conditioning of the function.

Non-Linear Transformations and Symmetry Breaking

In order to make relatively simple, but well understood functions less regular, on some functions non-linear transformations are applied in x-or f -space. Both transformations T osz : R n → R n , n ∈ {1, D}, and T asy : R D → R D are defined coordinate-wise (see below). They are smooth and have, coordinate-wise, a strictly positive derivative. They are shown in Figure 1. T osz is oscillating about the identity, where the oscillation is scale invariant w.r.t. the origin. T asy is the identity for negative values. When T asy is applied, a portion of 1/2 D of the search space remains untransformed.

Symbols and Definitions

Used symbols and definitions of, e.g., auxiliary functions are given in the following. Vectors are typeset in bold and refer to column vectors.

⊗ indicates element-wise multiplication of two

D-dimensional vectors, ⊗ : R D ×R D → R D , (x, y) → diag(x) × y = (x i × y i ) i=1,...,D
. denotes the Euclidean norm, x 2 = i x 2 i .

[.] denotes the nearest integer value 0 = (0, . . . , 0) T all zero vector 1 = (1, . . . , 1) T all one vector Λ α is a diagonal matrix in D dimensions with the ith diagonal element as λ ii = α

1 2 i-1 D-1 , for i = 1, . . . , D. f pen : R D → R, x → D i=1 max(0, |x i | -5) 2
1 + -a D-dimensional vector with entries of -1 or 1 with equal probability independently drawn.

Q, R orthogonal (rotation) matrices. For one function in one dimension a different realization for respectively Q and R is used for each instantiation of the function. Orthogonal matrices are generated from standard normally distributed entries by Gram-Schmidt orthonormalization.

Columns and rows of an orthogonal matrix form an orthonormal basis.

R see Q T β asy : R D → R D , x i → x 1+β i-1 D-1 √ xi i if x i > 0 x i otherwise
, for i = 1, . . . , D. See Figure 1.

T osz : R n → R n , for any positive integer n (n = 1 and n = D are used in the following), maps element-wise x → sign(x) exp (x + 0.049 (sin(c

1 x) + sin(c 2 x))) with x = log(|x|) if x = 0 0 otherwise , sign(x) =      -1 if x < 0 0 if x = 0 1 otherwise , c 1 = 10 if x > 0 5.

otherwise and

c 2 = 7.9 if x > 0 3.1 otherwise . See Figure 1.

x opt optimal solution vector, such that f (x opt ) is minimal.

1 Separable functions 1.1 Sphere Function

f 1 (x) = z 2 + f opt (1)
• z = x -x opt
Properties Presumably the most easy continuous domain search problem, given the volume of the searched solution is small (i.e. where pure monte-carlo random search is too expensive).

• unimodal

• highly symmetric, in particular rotationally invariant, scale invariant Information gained from this function:

• What is the optimal convergence rate of an algorithm?

1.2 Ellipsoidal Function

f 2 (x) = D i=1 10 6 i-1 D-1 z 2 i + f opt (2)
• z = T osz (x -x opt )
Properties Globally quadratic and ill-conditioned function with smooth local irregularities.

• unimodal

• conditioning is about 10 6

Information gained from this function:

• In comparison to f10: Is separability exploited?

1.3 Rastrigin Function

f 3 (x) = 10 D - D i=1 cos(2πz i ) + z 2 + f opt (3)
• z = Λ 10 T 0.2 asy (T osz (x -x opt ))

Properties Highly multimodal function with a comparatively regular structure for the placement of the optima. The transformations T asy and T osz alleviate the symmetry and regularity of the original Rastrigin function

• roughly 10 D local optima

• conditioning is about 10

Information gained from this function:

• In comparison to f15: Is separability exploited?

1.4 Büche-Rastrigin Function

f 4 (x) = 10 D - D i=1 cos(2πz i ) + D i=1 z 2 i + 100 f pen (x) + f opt (4)
• z i = s i T osz (x i -x opt i ) for i = 1 . . . D • s i = 10 × 10 1 2 i-1 D-1
if z i > 0 and i = 1, 3, 5, . . .

10 1 2 i-1 D-1 otherwise for i = 1, . . . , D
Properties Highly multimodal function with a structured but highly asymmetric placement of the optima. Constructed as a deceptive function for symmetrically distributed search operators.

• roughly 10 D local optima, conditioning is about 10, skew factor is about 10 in x-space and 100 in f -space

Information gained from this function:

• In comparison to f3: What is the effect of asymmetry?

1.5 Linear Slope

f 5 (x) = D i=1 5 |s i | -s i z i + f opt (5)
• z i = x i if x opt i x i < 5 2 and z i = x opt i otherwise, for i = 1, . . . , D. That is, if x i exceeds x opt i it will mapped back into the domain and the function appears to be constant in this direction.

• s i = sign x opt i 10 i-1 D-1 for i = 1, . . . , D. • x opt = z opt = 5 × 1 + -
Properties Purely linear function testing whether the search can go outside the initial convex hull of solutions right into the domain boundary.

• x opt is on the domain boundary

Information gained from this function:

• Can the search go outside the initial convex hull of solutions into the domain boundary?

Can the step size be increased accordingly?

2 Functions with low or moderate conditioning 2.6 Attractive Sector Function

f 6 (x) = T osz D i=1 (s i z i ) 2 0.9 + f opt (6)
• z = QΛ 10 R(xx opt )

• s i = 10 2 if z i × x opt i > 0 1 otherwise
Properties Highly asymmetric function, where only one "hypercone" (with angular base area) with a volume of roughly 1/2 D yields low function values. The optimum is located at the tip of this cone.

• unimodal

Information gained from this function:

• In comparison to f1: What is the effect of a highly asymmetric landscape?

2.7

Step Ellipsoidal Function

f 7 (x) = 0.1 max |ẑ 1 |/10 4 , D i=1 10 2 i-1 D-1 z 2 i + f pen (x) + f opt (7)
• ẑ = Λ 10 R(xx opt )

• zi = 0.5 + ẑi if [ẑ i -→] |ẑ i | > 0.5 0.5 + 10 ẑi /10 otherwise for i = 1, . . . , D, denotes the rounding procedure in order to produce the plateaus.

• z = Qz Properties The function consists of many plateaus of different sizes. Apart from a small area close to the global optimum, the gradient is zero almost everywhere.

• unimodal, non-separable, conditioning is about 100

Information gained from this function:

• Does the search get stuck on plateaus?

2.8 Rosenbrock Function, original

f 8 (x) = D-1 i=1 100 z 2 i -z i+1 2 + (z i -1) 2 + f opt (8)
• z = max 1,

√ D 8 (x -x opt ) + 1 • z opt = 1
Properties So-called banana function due to its 2-D contour lines as a bent ridge (or valley) [START_REF] Rosenbrock | An automatic method for finding the greatest or least value of a function[END_REF].

In the beginning, the prominent first term of the function definition attracts to the point z = 0. Then, a long bending valley needs to be followed to reach the global optimum. The ridge changes its orientation D -1 times. Exceptionally, here x opt ∈ [-3, 3] D .

• tri-band dependency structure, in larger dimensions the function has a local optimum with an attraction volume of about 25%

Information gained from this function:

• Can the search follow a long path with D -1 changes in the direction?

2.9 Rosenbrock Function, rotated

f 9 (x) = D-1 i=1 100 z 2 i -z i+1 2 + (z i -1) 2 + f opt (9)
• z = max 1,

√ D 8
Rx + 1/2

• z opt = 1
Properties rotated version of the previously defined Rosenbrock function.

Information gained from this function:

• In comparison to f8: Can the search follow a long path with D -1 changes in the direction without exploiting partial separability?

3 Functions with high conditioning and unimodal 3.10 Ellipsoidal Function

f 10 (x) = D i=1 10 6 i-1 D-1 z 2 i + f opt (10) • z = T osz (R(x -x opt ))
Properties Globally quadratic ill-conditioned function with smooth local irregularities, nonseparable counterpart to f 2 .

• unimodal, conditioning is 10 6

Information gained from this function:

• In comparison to f2: What is the effect of rotation (non-separability)?

3.11 Discus Function

f 11 (x) = 10 6 z 2 1 + D i=2 z 2 i + f opt (11) • z = T osz (R(x -x opt ))
Properties Globally quadratic function with local irregularities. A single direction in search space is a thousand times more sensitive than all others.

• conditioning is about 10 6

Information gained from this function:

• In comparison to f1: What is the effect of constraint-like penalization?

3.12 Bent Cigar Function

f 12 (x) = z 2 1 + 10 6 D i=2 z 2 i + f opt (12)
• z = R T 0.5 asy (R(x -x opt ))

Properties A ridge defined as D i=2 z 2 i = 0 needs to be followed. The ridge is smooth but very narrow. Due to T 1/2 asy the overall shape deviates remarkably from being quadratic.

• conditioning is about 10 6 , rotated, unimodal Information gained from this function:

• Can the search continuously change its search direction?

Sharp Ridge Function

f 13 (x) = z 2 1 + 100 D i=2 z 2 i + f opt (13)
• z = QΛ 10 R(xx opt )

Properties As for the Bent Cigar function, a ridge defined as D i=2 z 2 i = 0 must be followed. The ridge is non-differentiable and the gradient is constant when the ridge is approached from any given point. Following the gradient becomes ineffective close to the ridge where the ridge needs to be followed in z 1 -direction to its optimum. The necessary change in "search behavior" close to the ridge is difficult to diagnose, because the gradient towards the ridge does not flatten out. Information gained from this function:

• In comparison to f12: What is the effect of non-smoothness, non-differentiabale ridge?

Different Powers Function

f 14 (x) = D i=1 |z i | 2+4 i-1 D-1 + f opt (14) • z = R(x -x opt )
Properties Due to the different exponents the sensitivies of the z i -variables become more and more different when approaching the optimum.

4 Multi-modal functions with adequate global structure 4.15 Rastrigin Function

f 15 (x) = 10 D - D i=1 cos(2πz i ) + z 2 + f opt (15)
• z = RΛ 10 Q T 0.2 asy (T osz (R(x -x opt )))
Properties Prototypical highly multimodal function which has originally a very regular and symmetric structure for the placement of the optima. The transformations T asy and T osz alleviate the symmetry and regularity of the original Rastrigin function.

• non-separable less regular counterpart of f 3

• roughly 10 D local optima

• conditioning is about 10

Information gained from this function:

• in comparison to f3: What is the effect of non-separability for a highly multimodal function?

4.16 Weierstrass Function

f 16 (x) = 10 1 D D i=1 11 k=0 1/2 k cos(2π3 k (z i + 1/2)) -f 0 3 + 10 D f pen (x) + f opt (16)
• z = RΛ 1/100 Q T osz (R(x -x opt ))

• f 0 = 11 k=0 1/2 k cos(2π3 k 1/2)
Properties Highly rugged and moderately repetitive landscape, where the global optimum is not unique.

• the term k 1/2 k cos(2π3 k . . . ) introduces the ruggedness, where lower frequencies have a larger weight 1/2 k .

• rotated, locally irregular, non-unique global optimum Information gained from this function:

• in comparison to f17: Does ruggedness or a repetitive landscape deter the search behavior?

4.17 Schaffers F7 Function

f 17 (x) = 1 D -1 D-1 i=1 √ s i + √ s i sin 2 50 s 1/5 i 2 + 10 f pen (x) + f opt (17)
• z = Λ 10 Q T 0.5 asy (R(x -x opt ))

• s i = z 2 i + z 2 i+1 for i = 1, . .

. , D

Properties A highly multimodal function where frequency and amplitude of the modulation vary.

• asymmetric, rotated

• conditioning is low 4.18 Schaffers F7 Function, moderately ill-conditioned

f 18 (x) = 1 D -1 D-1 i=1 √ s i + √ s i sin 2 50 s 1/5 i 2 + 10 f pen (x) + f opt (18)
• z = Λ 1000 Q T 0.5 asy (R(x -x opt ))

• s i = z 2 i + z 2 i+1 for i = 1, . . . , D
Properties Moderately ill-conditioned counterpart to f 17

• conditioning of about 1000

Information gained from this function:

• In comparison to f17: What is the effect of ill-conditioning?

4.19 Composite Griewank-Rosenbrock Function F8F2

f 19 (x) = 10 D -1 D-1 i=1 s i 4000 -cos(s i ) + 10 + f opt (19)
• z = max 1,

√ D 8
Rx + 0.5

• s i = 100 (z 2 i -z i+1 ) 2 + (z i -1) 2 for i = 1, . . . , D • z opt = 1
Properties Resembling the Rosenbrock function in a highly multimodal way.

5 Multi-modal functions with weak global structure 5.20 Schwefel Function

f 20 (x) = - 1 100D D i=1 z i sin( |z i |) + 4.189828872724339 + 100f pen (z/100) + f opt (20) • x = 2 × 1 + -⊗ x • ẑ1 = x1 , ẑi+1 = xi+1 + 0.25 xi -[x opt i -→] 2 x opt i for i = 1, . . . , D -1 • z = 100 (Λ 10 (ẑ -[x opt -→] 2 |x opt |) + [x opt -→] 2 |x opt |)
• x opt = 4.2096874633/2 1 + -, where 1 + -is the same realization as above Properties The most prominent 2 D minima are located comparatively close to the corners of the unpenalized search area, based on [START_REF] Schwefel | Numerical optimization of computer models[END_REF].

• the penalization is essential, as otherwise more and better minima occur further away from the search space origin 5.21 Gallagher's Gaussian 101-me Peaks Function

f 21 (x) = T osz 10 - 101 max i=1 w i exp - 1 2D (x -y i ) T R T C i R(x -y i ) 2 + f pen (x) + f opt (21) • w i =    1.1 + 8 × i -2 99 for i = 2, . . . , 101
10 for i = 1
, three optima have a value larger than 9

• C i = Λ αi /α 1/4 i
where Λ αi is defined as usual (see Section 0.2), but with randomly permuted diagonal elements. For i = 2, . . . , 101, α i is drawn uniformly randomly from the set 1000 2 j 99 | j = 0, . . . , 99 without replacement, and α i = 1000 for i = 1.

• the local optima y i are uniformly drawn from the domain [[-4.9, 4.9] D -→] [-5, 5] D for i = 2, . . . , 101 and y 1 ∈ [-4, 4] D . The global optimum is at x opt = y 1 .

Properties The function consists of 101 optima with position and height being unrelated and randomly chosen (different for each instantiation of the function), based on [START_REF] Gallagher | A general-purpose tunable landscape generator[END_REF].

• the conditioning around the global optimum is about 30

Information gained from this function:

• Is the search effective without any global structure? , two optima have a value larger than 9

• C i = Λ αi /α 1/4 i
where Λ αi is defined as usual (see Section 0.2), but with randomly permuted diagonal elements. For i = 2, . . . , 21, α i is drawn uniformly randomly from the set 1000 2 j 19 | j = 0, . . . , 19 without replacement, and α i = 1000 2 for i = 1.

• the local optima y i are uniformly drawn from the domain [-4.9, 4.9] D for i = 2, . . . Properties The function consists of 21 optima with position and height being unrelated and randomly chosen (different for each instantiation of the function), based on [START_REF] Gallagher | A general-purpose tunable landscape generator[END_REF].

• the conditioning around the global optimum is about 1000

Information gained from this function:

• In comparison to f21: What is the effect of higher condition?

• z = Q Λ 100 R(x -x opt )

Figure 1 :

 1 Figure 1: T osz (blue) and D-th coordinate of T asy for β = 0.1, 0.2, 0.5 (green)

5. 22 2 +

 222 Gallagher's Gaussian 21-hi Peaks Functionf 22 (x) = T osz 10 -y i ) T R T C i R(xy i ) f pen (x) + f opt

, 21 and y 1

 1 ∈[[-4, 4] D -→] [-3.92, 3.92] D ]. The global optimum is at x opt = y 1 .

2 j

 2 z i -[2 j z i ] 2 j   10/D 1.2 -10 D 2 + f pen (x)+f opt(23)

  General Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0.2 Symbols and Definitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Sphere Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Ellipsoidal Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Rastrigin Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Büche-Rastrigin Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Linear Slope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	Contents
	0 Introduction
	0.1 1 Separable functions
	1.1

2 Functions with low or moderate conditioning 2.6 Attractive Sector Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 Step Ellipsoidal Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.8 Rosenbrock Function, original . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.9 Rosenbrock Function, rotated . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Functions with high conditioning and unimodal 3.10 Ellipsoidal Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.11 Discus Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.12 Bent Cigar Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.13 Sharp Ridge Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.14 Different Powers Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Multi-modal functions with adequate global structure 4.15 Rastrigin Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.16 Weierstrass Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.17 Schaffers F7 Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.18 Schaffers F7 Function, moderately ill-conditioned . . . . . . . . . . . . . . . . . . 4.19 Composite Griewank-Rosenbrock Function F8F2 . . . . . . . . . . . . . . . . . . . 5 Multi-modal functions with weak global structure 5.20 Schwefel Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.21 Gallagher's Gaussian 101-me Peaks Function . . . . . . . . . . . . . . . . . . . . . 5.22 Gallagher's Gaussian 21-hi Peaks Function . . . . . . . . . . . . . . . . . . . . . . . 5.23 Katsuura Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.24 Lunacek bi-Rastrigin Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A Function Properties A.1 Deceptive Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.2 Ill-Conditioning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.3 Regularity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.4 Separability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.5 Symmetry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.6 Target function value to reach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Errata Errata (2019) are marked with colored text.

RR n°6829

For our experimental setup see[START_REF] Hansen | Real-parameter black-box optimization benchmarking[END_REF][START_REF] Hansen | COCO: the experimental procedure[END_REF] and for our performance assessment methodology see[START_REF] Hansen | Coco: performance assessment[END_REF].

The implementation provides an instance ID as input, such that a set of uniquely specified instances can be explicitly chosen.

Acknowledgments

The authors would like to thank Arnold Neumaier for his instructive comments, Verena Heidrich-Meisner for her contribution on information gain, and Dimo Brockhoff and Duc Manh Nguyen for relentlessly spotting mistakes. Steffen Finck was supported by the Austrian Science Fund (FWF) under grant P19069-N18.

Properties Highly rugged and highly repetitive function with more than 10 D global optima, based on the idea in [START_REF] Katsuura | Continuous nowhere-differentiable functions-an application of contraction mappings[END_REF]. [START_REF] Hansen | COCO: the experimental procedure[END_REF].24 Lunacek bi-Rastrigin Function

cos(2πz i ) +10 4 f pen (x)+f opt (24)

Properties Highly multimodal function based on [START_REF] Lunacek | The impact of global structure on search[END_REF] with two funnels around [µ

-being superimposed by the cosine. Presumably different approaches need to be used for "selecting the funnel" and for search the highly multimodal function "within" the funnel. The function was constructed to be deceptive for evolutionary algorithms with large population size.

• the funnel of the local optimum at [-µ All "deceptive" functions provide, beyond their deceptivity, a "structure" that can be exploited to solve them in a reasonable procedure.

A.2 Ill-Conditioning

Ill-conditioning is a typical challenge in real-parameter optimization and, besides multimodality, probably the most common one. Conditioning of a function can be rigorously formalized in the case of convex quadratic functions, f (x) = 1 2 x T Hx where H is a symmetric definite positive matrix, as the condition number of the Hessian matrix H. Since contour lines associated to a convex quadratic function are ellipsoids, the condition number corresponds to the square root of the ratio between the largest axis of the ellipsoid and the shortest axis. For more general functions, conditioning loosely refers to the square of the ratio between the largest direction and smallest of a contour line. The testbed contains ill-conditioned functions with a typical conditioning of 10 6 . We believe this a realistic requirement, while we have seen practical problems with conditioning as large as 10 10 .

A.3 Regularity

Functions from simple formulas are often highly regular. We have used a non-linear transformation, T osz , in order to introduce small, smooth but clearly visible irregularities. Furthermore, the testbed contains a few highly irregular functions.

A.4 Separability

In general, separable functions pose an essentially different search problem to solve, because the search process can be reduced to D one-dimensional search procedures. Consequently, nonseparable problems must be considered much more difficult and most benchmark functions are designed being non-separable. The typical well-established technique to generate non-separable functions from separable ones is the application of a rotation matrix R.

A.5 Symmetry

Stochastic search procedures often rely on Gaussian distributions to generate new solutions and it has been argued that symmetric benchmark functions could be in favor of these operators. To avoid a bias in favor of highly symmetric operators we have used a symmetry breaking transformation, T asy . We have also included some highly asymmetric functions.

A.6 Target function value to reach

The typical target function value for all functions is f opt + 10 -8 . On many functions a value of f opt + 1 is not very difficult to reach, but the difficulty versus function value is not uniform for all functions. These properties are not intrinsic, that is f opt + 10 -8 is not intrinsically "very good". The value mainly reflects a scalar multiplier in the function definition.