S Dupire

D Toe

J.-B Barré

F Bourrier

F Berger

Harmonized mapping of forests with a protection function against rockfalls over European Alpine countries

Keywords: Protection forests, Rockfall, Mapping, Ecosystem services, Alps, Natural hazard Highlights

Forest covers 40% of the European Alpine region and contributes to the protection of human beings and infrastructures against natural hazards such as rockfalls. However, despite the recognition of this ecosystem service, most mountain territories do not have a map of protection forests. When a map exists, it generally depends on data restricted to a limited extent, which prevents any replication or comparison on other areas.

The aim of this study is to develop a method using harmonized and open data to produce the first map of protection forests against rockfalls at the European Alpine region scale. Based on these data, we first identified potential rockfall release areas and calibrated the model according to 2812 real rockfall events located around the Alps. Second, 46.5 billion 3-D rockfall propagation simulations, taking into account topography, land use and human assets, were computed on the entire area. Protection forest is defined as being located on at least one rockfall trajectory that had impacted human assets.

Our results show that 14% of the forests have a potential protection function against rockfalls in the Alpine Space region. This proportion goes up to 21.5%, if we consider only the core of the Alpine area. 80% of the protection forest area contributes to mitigate rockfall hazard on road network, 55% on buildings and only 6% on railways. This work provides a robust, objective and reproducible method for locating protection forests on a large geographical scale. Such a map may serve as a basis for national and European risk management policies.

• We identified harmonized open data to model rockfall propagation over European Alps • 46.5 billions rockfall trajectories were computed to identify protection forests • 14% of protection forests in the Alpine region (21.5% in the core of the Alpine area) • 80% of protection forests area protect roads, 55% buildings and 6% railways • Reproducible method and result available for policy makers, foresters and risk managers

Introduction

Mountain forests provide many ecosystem services [START_REF] Briner | Trade-Offs between Ecosystem Services in a Mountain Region[END_REF], including the protection of human beings and infrastructures against natural hazards. In the Alpine region, the role of forests is often considered equivalent to civil engineering measures as it represents a cost-effective solution [START_REF] Notaro | The economic valuation of natural hazards in mountain forests: an approach based on the replacement cost method[END_REF] especially to mitigate avalanches, soil erosion or rockfalls [START_REF] Bebi | Assessing structures in mountain forests as a basis for investigating the forests' dynamics and protective function[END_REF][START_REF] Brang | Forests as Protection from Natural Hazards[END_REF].

The preservation and enhancement of the role of forests in mitigating natural hazards is essential in the strategies for the protection of inhabitants, users and economic activities of the mountain territories [START_REF] Accastello | A framework for the integration of nature-based solutions into environmental risk management strategies[END_REF]. The ongoing climate changes may promote disturbances such as fires [START_REF] Dupire | Vulnerability of forest ecosystems to fire in the french alps[END_REF][START_REF] Dupire | Spatio-temporal trends in fire weather in the french alps[END_REF][START_REF] Conedera | Characterizing alpine pyrogeography from fire statistics[END_REF], wind-storm [START_REF] Csilléry | Coupled effects of wind-storms and drought on tree mortality across 115 forest stands from the western alps and the jura mountains[END_REF] and drought [START_REF] Gobron | The state of vegetation in Europe following the 2003 drought[END_REF]. In order to avoid their catastrophic consequences related to the alteration, or even the disappearance, of the forest cover [START_REF] Gehring | Shallow landslide disposition in burnt European beech (Fagus sylvatica L.) forests[END_REF], it is first of all necessary to locate the forests that have a protective role, which requires a cartographic zoning of this ecosystem service [START_REF] Lepuschitz | Geographic information systems in mountain risk and disaster management[END_REF].

Different regional-scale maps of rockfall hazard have been proposed [START_REF] Toe | Regional Mapping of Forest with a Protection Function Against Rockfall[END_REF][START_REF] Acosta | Regional rockfall hazard assessment in the benasque valley (central pyrenees) using a 3d numerical approach[END_REF]), but they are not adapted to identify the forests with a potential effect on rockfall hazard mitigation. On the other hand, Switzerland [START_REF] Losey | Forêt protectrice en suisse. du projet SilvaProtect-CH à la forêt protectrice harmonisée[END_REF] and more locally certain administrative regions in the Alpine region [START_REF] Getzner | Gravitational natural hazards: Valuing the protective function of alpine forests[END_REF][START_REF] Meloni | Analisi della funzione protettiva delle foreste: l'esempio della "Carta delle foreste di protezione diretta della Valle d'Aosta[END_REF] have maps of protection forests but the natural hazard(s) against which forests protect are not always well informed. In addition, these maps are developed using different methods and input data from different sources making it nearly impossible to compare the results from one area to another. Thus, despite the recognition of the protection effect of forest against rockfall in both national and European forest policies, there is currently no harmonized method for carrying out such mapping.

In this context, our study shares the objective of the Alpine Space project "RockTheAlps"1 which aims to produce the first harmonized map of protection forests against rockfalls at the European Alpine Region scale. For this purpose, we first identified opensource and harmonized spatial data available for the entire study area. Second, 3-D rockfall simulations were conducted on 6 pilot areas in order to calibrate the rockfall model parameters. Third, rockfall simulations were computed on the whole study area and allowed to identify and locate the protection forests. Finally, we analysed the proportion of protection forests and human assets at risk at different scales.

Material and methods

Study area

This study focuses on the European Alps which geographical extent may be defined according to geographical, economical and cultural criteria (Fig. 1). Two areas of analysis were used in this work. The first one corresponds to the Alpine Region as defined by the European Alpine Space programme (www.alpine-space.eu). This area is based on administrative borders and is characterized by 188 Units for Territorial Statistics (referred as NUTS3) covering 390 227 km 2 over 8 countries (Austria, France, Germany, Italy, Liechtenstein, Monaco, Slovenia and Switzerland). NUTS3 refers to the hierarchical system used to divide up the economic territory of the European Union for the collection, development, and harmonization of socio-economic analyses of small regions (< 800 000 inhabitants).

The second area of analysis is defined by the Alpine Convention (www.alpconv.org) whose extent better fits to the orographic contours of the Alps. This area covers 190 717 km 2 and excludes the flattest regions present in the Alpine Space programme (e.g. Po Plain and Rhone Valleys).

General model description

This section describes the general workflow (Fig. 2) used to establish the harmonized mapping of protection forests against rockfalls and specifies the inputs and outputs of the model. Details on the process itself are given in section 2.3.

Input sources

Spatial input data were chosen to be harmonized over the entire study area and in open access. Therefore, they were taken from two main sources: the European Union's Earth Observation Programme Copernicus2 and Open-StreetMap3 .

Relief and topography were obtained from the EU-DEM4 v1.1 which is a digital terrain model (DTM) available in raster format for Europe at 25m spatial resolution with an overall vertical accuracy of 2.9 meters RMSE. EU-DEM was used as is to account for relief and topography in the modelling process. A slope raster was also produced from this source in order to identify the potential rock release areas (see section 2.3.1).

Land and forest cover were extracted from Corine Land Cover (CLC) 20185 [START_REF] Büttner | CORINE land cover 2018 technical guidelines[END_REF]. CLC is produced by visual interpretation of high resolution satellite imagery. It is available in shapefile (vector layer of polygons) format with a minimum mapping unit of 25 ha. The forest area input was derived from CLC by extracting entities with the CLC level 3 codes in [311,312,313,323,324]. CLC was also used to develop the classification of soil types (see section 2.3.3).

Water surfaces and rivers were obtained from © OpenStreetMap contributors in shapefile format and used also to develop the classification of soil types.

Human assets were taken from © OpenStreetMap contributors locating different topographical information over the world. The assets used were of three types (see Table 1) and were saved in four layers in shapefile format. No detailed information was available in OpenStreetMap in order to distinguish the different types of buildings (e.g. houses, schools, industry, etc.). Thus, we assumed the same vulnerability to rockfall for all buildings. Moreover, there was no information about linear assets that are already protected from rockfall with permanent civil engineering fences such as embankments. Therefore, linear assets were assumed to be vulnerable to rockfall on their entire length. Human asset shapefiles were rasterized at the EU-DEM resolution for the modelling process and were also used for the elaboration of soil classes. The rasterization process allowed to keep a trace of the type of asset present on each Highways + national or regional roads Line 100 raster cell. For example, if buildings and road are present on the same raster cell, it will have the code 110 (10 for Buildings + 100 for Roads).

Process

The propagation on the slope of 500 rocks with a random volume between [0.1,5] m 3 is modeled from each pixel identified as release area (see section 2.3.1). This range of volume corresponds to the optimal efficiency of forests against rockfalls [START_REF] Berger | Forest: a natural protection mean against rockfalls, but with which efficiency[END_REF]. The propagation of rocks down the slope is simulated on the rasterized digital terrain model (EU-DEM) by successive sequences of free flights through the air and rebounds on the slope surface (see section 2.3 for detailed description). For each rockfall simulation, the rock volume was randomly taken in the previous range and outputs relative to rockfall activity are systematically saved on the rock trajectory. If the rock trajectory intersects a human asset, all the pixels of forest located between the asset and the release area are mapped as protection forest.

Outputs

The simulation outputs are made up of two sets of results. The first set gathers Rockfall hazard outputs which consist in three 25 m resolution rasters. It is only used to calibrate (see section 2.4) and validate (see section 2.5) the modelling process. The first raster returns the number of passing rocks for each raster cell located on a simulated rock trajectory. The second refers to the raster cells where rocks ended their courses. Finally, the last raster gives the maximum energy line angle (see section 2.4 for definition) observed for each cell with a passing rock. It is important to underline that these data are adapted to identify forests with a protective function but not intended to be used for hazard zoning or risk assessment.

The second set concerns Protection forest outputs which are the final expected results of this study. They consists in three datasets. The first one is a shapefile locating the protection forests. The second is a raster at 25 m resolution returning the human assets type protected by the forest. Finally, a 25 m resolution raster inform on the minimum length of forest present between the release areas and the human assets for each pixel of protection forest.

Detailed model description 2.3.1. Rockfall source area identification

A classical approach to determine the potential rockfall sources consists in identifying the areas with a slope angle greater than a given threshold value [START_REF] Frattini | Assessment of rockfall susceptibility by integrating statistical and physically-based approaches[END_REF][START_REF] Dorren | Comparison of three gis-based models for predicting rockfall runout zones at a regional scale[END_REF]. The value of this threshold angle highly depends on the geomorphology as well as the quality and resolution of the slope raster [START_REF] Michoud | Rockfall hazard and risk assessments along roads at a regional scale: example in swiss alps[END_REF]. For example, a vertical cliff of 10 m will be represented by slope values between 82°and 84°o n a 1 m grid cell (depending on the cliff aspect) and slopes values between 16°and 22°on a 25 m grid cell [START_REF] Loye | Identification of potential rockfall source areas at a regional scale using a dem-based geomorphometric analysis[END_REF]. Therefore, standard DEM (25 m resolution), such as the EU-DEM, considerably smooths the relief which makes difficult the identification of precise rockfall source areas. To overcome this issue, six calibrations areas spread over the Alpine Space (Table . 2) were used to identify the slope angle threshold to be used at the Alpine Space level. 1 m DEM derived from Lidar was available for each area and used to identify potential rockfall source pixels using a threshold slope angle of 50°. This value corresponds to the average obtained on the different geomorphological areas analysed at 1 m resolution by [START_REF] Loye | Identification of potential rockfall source areas at a regional scale using a dem-based geomorphometric analysis[END_REF]. These potential source areas were rasterized at the same resolution as slope raster derived from EU-DEM (25 m) in order to extract the slope angle distribution on their location. The distributions obtained considerably vary between the calibration zones (Table . 2) which conducted us to choose two thresholds to be applied at the Alpine Space scale that correspond to the 25 th and the 75 th percentiles of the distribution considering all the calibration zones. The first threshold of 42°is more likely to identify release areas in mountainous terrain and is in agreement with the values found in the literature for mountain areas which stand in the range [40-45°] at a 25 m resolution [START_REF] Dumperth | An easy approach to assess the susceptibility of a landslide by utilizing simple raster overlay methods: A case study on huangtupo landslide (p.r. china)[END_REF][START_REF] Loye | Identification of potential rockfall source areas at a regional scale using a dem-based geomorphometric analysis[END_REF][START_REF] Acosta | Regional rockfall hazard assessment in the benasque valley (central pyrenees) using a 3d numerical approach[END_REF][START_REF] Dorren | Comparison of three gis-based models for predicting rockfall runout zones at a regional scale[END_REF]. The second threshold of 28°is mostly interesting in mid-mountains and pre-Alps where the EU-DEM relief may be too smoothed to identify source areas as shown by the slope angle distribution of the FR-Bas-Rhin (Vosges) calibration area.

Ultimately, 9756 km 2 (15.6 million pixels) for release areas on slopes greater than or equal to 42°and 48388 km 2 (77.4 million pixels) for release areas on slope in the range [28-42°[ were detected on the entire Alpine Space region. As 500 rockfall simulations are launched by release area, it resulted in 46.5 billion rockfall trajectory simulations.

Model of rock propagation and rebounds on the slope surface

The block propagation is modeled by a succession of free flights and impacts on the slope surface. During the free flight phase, as the block is only subjected to gravity, its trajectory can be described as a parabola depending on the velocities at the beginning of the flight phase and on the acceleration due to gravity (g) [START_REF] Bourrier | Rockfall rebound: comparison of detailed field experiments and alternative modelling approaches[END_REF].

The computation of the rebound on the soil is made of two steps. The first step enables to calculate the block deviation after impact according to its incoming trajectory and the slope [START_REF] Dorren | Real-size experiments and 3-D simulation of rockfall on forested and non-forested slopes[END_REF]. The second step enables to calculate the translational and rotational velocities of the block after an impact on the soil depending on the incident block velocities and soil classes. It is based on the rebound model developed by [START_REF] Pfeiffer | Computer Simulation of Rockfalls[END_REF]. The normal and the tangential components of the velocity after impact are calculated as functions of the normal, tangential and rotational components of the velocity before impact, two restitution coefficients (R n and R t ), and the block characteristics (mass, diameter and inertia). R n values are directly associated with the 7 soil classes proposed. R t is calculated using the roughness of the soil (R g ) the block radius, and the penetration depth of the block in the soil [START_REF] Dorren | Real-size experiments and 3-D simulation of rockfall on forested and non-forested slopes[END_REF][START_REF] Dorren | Combining field and modelling techniques to assess rockfall dynamics on a protection forest hillslope in the european alps[END_REF]. Consequently, the block propagation is depending on two parameters related to soil properties: the normal restitution coefficient R n and the soil roughness R g .

Definition of soil classes

Soil parameter values used in rockfall propagation model are often derived from the type of land use [START_REF] Dorren | Rockyfor3D (v5.2) revealed -Transparent description of the complete 3D rockfall model[END_REF][START_REF] Volkwein | Rockfall characterisation and structural protection -a review[END_REF]. Corine Land Cover provides a homogeneous description of land use at the European scale and it is therefore the main input of the soil classes proposed in this study (Fig. 3). 8 soil classes (SC) were defined from CLC and precised according to elevation, slope or vector information on human assets, water surfaces and rivers. The SC (I) corresponds to agricultural area characterized by deep soils with a high capacity to reduce rockfall velocity. SC (II) are associated to forest area defined by a high diversity of soil structures in terms of depth, damping capacity and surface roughness. SC (III) represents mineral area, glacier and high elevation meadows. It is characterized by a wide range of soil parameters promoting rock propagation. SC (IV) gathers CLC codes linked to human activities and urban areas. SC (V) concerns permanent water area and rivers where rock propagation is stopped (R n =0). SC (VI) corresponds to temporary water area such as bog or mountain torrent with intermittent regime and are characterized by a low rock propagation. SC (VII) are associated to screes -i.e. CLC codes corresponding to mineral surface but located on slopes in the range [32,40[°(Francou and Manté, 1990) -and are defined by a high surface roughness. Finally SC (VIII) represents meadows, shrub-lands and grasslands from low to middle elevation (< 2000 m) characterized with soil depths and rock propagation in between those of forests and mineral areas.

Identification of potential protection forests

At the beginning of the simulation process, three rasters of the same size and resolution as the DEM are created. The first (PF) is initialized to zero and will be used to identify the protection forest pixels (which will be coded 1). The second raster (LMIN) will store the minimum length of forest on the slope between the source area and the human assets. The last one (PHA) is initialized to zero and will identify the types of human assets protected by the forest.

During the propagation simulation, the trajectory of each block (traj) is temporarily recorded as it progresses. Each time a raster cell is added to the trajectory, a test is performed to determine if human assets are located at that position. If yes, all forest pixels located on the path upstream of the human asset (trajFor) are identified as protective forest: PF[trajFor]=1. The length of forest between the starting area and the stake (L) is also calculated. The minimum values between L and those contained in LMIN for each forest cell on the trajectory are kept: LMIN[trajFor] = min(LMIN[trajFor], L). Finally, the code of the human asset type (see Table . 1) is added to each PHA forest raster cells on the trajectory: PHA[trajFor]=PHA[trajFor]+code. With this process, if a forest protects both a building (code 10) and a road (code 100), PHA will get the code 110.

Calibration of soil parameters for rockfall propagation

The calibration of rockfall propagation model consisted in defining the range of normal restitution coefficient R n and soil roughness R g for each soil class based on the distribution of the energy line angle. The energy line principle has been first described by [START_REF] Heim | Bergsturz und menschenleben[END_REF] before being used in several rockfall applications [START_REF] Toe | Regional Mapping of Forest with a Protection Function Against Rockfall[END_REF][START_REF] Berger | Principles of the tool Rockfor.net for quantifying the rockfall hazard below a protection forest[END_REF]. The energy line is characterized by the angle β between the horizontal plane and a fictive line drawn between the top of a rockfall source to the stopping point β is easy to measure on the field and can provide a rough statistical estimation of the rockfall run-out distance. In the framework of the Interreg -Alpine Space project RockTheAlps, a collection of rockfalls events across the Alps is carried out. This database (still in progress) gathers information on rockfalls events such as date of the event, rock volume, energy line angle β, location of source, location of stop and complete ground profile (not all the information is systematically known for each event). At the time of the study, β was available for 2812 rockfall events.

The model returns as output a raster of maximum energy line angle (β max ) which identifies for each cell on a block trajectory the maximum of all energy line angles β observed. Another output is a raster of stopped rocks (SR) returning the number of rocks that stopped their courses in each cell. Therefore, for each simulation, it is possible to obtain a distribution of the β angle by extracting the values of β max on the cells with SR >0.

The objective of the calibration step was to obtain a high level of similarity between observed distribution (β values observed on 2812 real rockfall events) and the modeled distribution calculated from β max and SR obtained with the simulations of rock propagation from release areas on slopes ≥ 42°(the most suitable to Alpine region, see section 2.3.1) on the 6 areas defined in Table .2. The similarity was assessed with two different criteria. First, we used a Wilcoxon rank sum test [START_REF] Fretz | Nonparametric statistical methods[END_REF][START_REF] Wilcoxon | Individual comparisons by ranking methods[END_REF] in order to center the modeled distribution on the observed values. Second, the Bhattacharyya coefficient (BC) [START_REF] Bhattacharyya | On a measure of divergence between two multinomial populations[END_REF] was applied to calculate the probability of distribution overlap. Its expression is given by:

BC(p, q) = p(x) • q(x) dx (1)
where p(x) and q(x) are the probability density functions to be compared. BC takes values between 0 (distributions do not overlap) and 1 (complete similarity).

In practice, different rockfall propagation simulations were carried out on the 6 calibration areas where soil parameters R n and R g were adjusted in order to have a conclusive Wilcoxon test (p-value ≫ 0.05) and a BC>0.95 for simulations from release areas located on slopes ≥ 42°.

Validation of soil parameters for rockfall propagation

The modeled distributions of the energy line angle β on the entire Alpine Space area were analysed and compared to the same observed distribution used for model calibration. Modeled distributions were split according to the release area source and their similarity with the observed values was assessed using Wilcoxon rank sum test and Bhattacharyya coefficient as described previously in section 2.4.

In addition to the previous point, we intersected the modeled propagation area (output raster of the number of passing rocks PR) with the locations of the deposits of real rockfalls events from the RockTheAlps database. 3685 events located in France, Italy and Slovenia were available at the time of this study. The distance between each real rockfall events to the nearest modeled propagation area was also computed in order to analyse the different plausible causes of non-intersection.

Results

Calibration and validation of the rockfall propagation model 3.1.1. Analyse of energy line angle distributions

Fig. 4 and Table . 3 return the results of the calibration and the validation steps regarding the distribution of the energy line angle β. .2 in dotted lines) and validation sets (entire Alpine Space region in full line). Cumulative probability curves (a.) and comparisons of density curves between observed and modeled β distributions for release areas in slopes ≥ 42°(b.) and slopes in [28-42°[ (c.) are displayed.

According to the statistical tests, observed and modeled distributions for release areas in slopes ≥ 42°have similar distributions. Non-significant difference are observed in median values (0.2°for calibration set -0.1°for validation set) as well as a high probability of distribution overlap (for both calibration and validation BC > 0.95). The standard deviation is however greater in the observed (8.15°) than in the modeled distributions (5.97°for calibration set -5.88°8 for validation set). These results suggests that the soil parameters used to model the propagation of the rocks are suitable for the Alpine region.

In contrast to the modeled distributions for release areas in slopes ≥ 42°, those for release areas in slopes in the range [28-42°[presented lower values of β angle than to the observed distribution (difference in median of -6.16°for calibration set and -7.21°for validation set). The modeled β angle values are also less scattered than observed values.

Final parameters of the soil classes

The ranges of R n and R g obtained after calibration process and used for the mapping of protection forests on the whole study area is shown in Table . 4. 3.1.3. Analysis of the run-out distance: comparison with real rockfall events 90% of the 3685 real rockfall events were located inside the modeled rockfall propagation area (Fig. 5). 65% of the real events intersected the propagation areas from release areas in slopes ≥ 42°and 25% in the propagation areas from release areas in slopes in [28-42°[. Only 10% of the real rockfall events did not intersect modeled rockfall propagation areas. Nonetheless, about half of them were located close to a propagation area, e.g. within a distance less or equal than the raster resolution of 25 m.

Results of the harmonized mapping of protection forest against rockfalls

All mapping results as well as thematic maps are freely available for downloading on the RockTheAlps project website.

When all release areas are considered, 92% of the protection forests of the Alpine Space region are located into the restricted perimeters of the Alpine Convention. This proportion is higher with 97% of the protection forests when only release areas on slopes greater than or equal to 42°are considered.

Results at the country level

Generic results at the country level are provided in Table . 5. The proportion of protection forest against rockfalls on the whole Alpine Space is in between 4.9% (while considering only potential release areas on slopes ≥ 42) and 14% (all potential release areas included). The proportion is greater on the Alpine Convention region with a proportion in between 7.9% and 21.5% (according to the potential release areas considered).

The most mountainous countries (Switzerland, Italy and Liechtenstein) present the highest rate of protection forest against rockfall (>25% considering all release areas) while the flattest (Germany) presents only 3.3% of protection forest in the Alpine Space region. The proportion of protection forests are higher for all countries while considering only the Alpine Convention region (except Liechtenstein which is unchanged because completely included in the Alpine Convention), however, the hierarchy is unchanged with the mountainous countries presenting the highest rates. Proportions are calculated at the NUTS3 levels for both regions.

Which human assets are most often protected by forests?

In the Alpine region, road network is the human asset the most often encountered underneath protection forests. Thus, about 80% of the total area of the protection forests are located above roads (Table. 6). Buildings are the second human asset type the most frequently protected by forests with about 55% of total area. Railways comes last with only 6% of the total protection forest area. If 63% of the protection forests concern a single human asset type, 34% of the them involve two different asset types and about 3% the three human asset types at the same time. The total proportion of human assets potentially exposed to rockfall hazard in the Alpine Space region is relatively low with only 3.4% of buildings and 6% of railways and roads. However, when considering only the Alpine Convention region, the proportion is much higher with about 10% of the buildings and 17% of roads and highways.

Discussion

The method developed in this study allowed to build up the first harmonized map of protection forests against rockfalls for the whole European Alpine region. Established from open source data, the process is easily reproducible and could be replicated for others mountainous region as long as harmonized data of relief, land use and human assets are available. Based on harmonized inputs, the results authorize for the first time the comparison of the importance of protection forests against rockfall between the different countries of the Alpine region or smaller administrative entities such as NUTS3.

The presented approach has certain advantages over previous regional mappings of protection forests against rockfalls. For example, Switzerland is the only country to have a nationwide map of protection forests [START_REF] Losey | Forêt protectrice en suisse. du projet SilvaProtect-CH à la forêt protectrice harmonisée[END_REF][START_REF] Zinggeler | Steinschlagsimulation in Gebirgswaldern. Modellierung der relevanten Teil-prozesse[END_REF] in the Alpine region. However, this map lies on a single rockfall simulation per source area. Our method uses a large number of rockfall simulations per source area (i.e. 500) in order to reflect the inherent randomness of any rockfall event [START_REF] Bourrier | Toward objective rockfall trajectory simulation using a stochastic impact model[END_REF] and to produce a more realistic covering of plausible rockfall trajectories. Another mapping method has been proposed by [START_REF] Toe | Regional Mapping of Forest with a Protection Function Against Rockfall[END_REF] and consists in extracting all straight terrain profiles from a digital terrain model from release area to the surrounding pixels that respect a threshold of energy line angle. Even if this allows fast calculation times, it does not account for topography and land use type in rock propagation and can lead to a wrong estimation of rockfall run out zones. In contrast, both topography and land use are essential to the process introduced in this paper. Land use participates to the calibration of soil parameters for the rockfall propagation, both topography and land use intervene in the validation of the simulations according to real rockfall events which ensures a robust modeling of potential rockfall trajectories.

© OSM Contributors

Real rockfall events

Outside propagation area

Propagation area

Release areas: slopes in [28-42°[ The accurate detection of rockfall release areas is one of the biggest challenges of a rockfall analysis on large areas [START_REF] Michoud | Rockfall hazard and risk assessments along roads at a regional scale: example in swiss alps[END_REF]. The approach presented here integrates several slope thresholds for the detection of release areas. Rockfalls from release areas on slopes greater than or equal to the first threshold of 42°are more susceptible to occur (Fig. 5: about 65% of the real rockfalls events) especially in the most mountainous regions. The second threshold of 28°is mostly interesting in order to detect the release areas corresponding to small rocky barriers on short slopes which are poorly represented with a raster resolution of 25 m. This type of configuration is mainly found in the pre-Alps or in the middle mountains. Despite this, some release areas might not be detected as illustrated in Fig. 8 which explains that less than 5% of the real rockfall events present in the database were at a distance of the modeled rockfall propagation areas greater than 25 m. At the opposite, a part of the release areas identified in this work might not be active due to stable bedrock. To overcome this issue, it is possible to use a self-produced spatial layer of release areas instead of the default layer produced during the process.

The overall calculation time required to produce this map at the Alpine Space scale is about 256 days without multi-threading (43 days when considering only release areas on slopes greater than or equal to 42°). This figure must be qualified by the fact that 46.5 billion single simulations of rockfall propagation were computed which finally returns only 475 µs per simulation. In order to overcome this issue, the computations were split on several threads on a server which allows us to get the results within two weeks of calculation. Moreover, it has to be said that the map does not need to be often updated, the same frequency as Corinne Land Cover updates should be sufficient (i.e every 6 years) in most of the Alpine region. However, in areas with very dynamic development, the appearance of new human assets might imply that some forests become protection forests.

This work emphasizes that a significant proportion of mountain forests protect against rockfalls (up to 21.5% in the Alpine Convention region). However, managing these forests represents a cost that forest owners should not bear alone in view of the service it represents for mountain societies [START_REF] Bianchi | The economic evaluation of forest protection service against rockfall: A review of experiences and approaches[END_REF]. By providing an objective and high resolution map of the forests having a potential to mitigate rockfall hazard at a large scale, this work could serve as a basis to implement a system of subsidies at the European scale to efficiently manage those forests and ensure that their protection capability will last or even improve in time. Such a system should rely on an objective evaluation of the protection capability of each forest identified here by using quantitative indicators (Dupire et al., 2016a,b). Remote sensing methods such as Lidar could also be of great interest in this topic [START_REF] Monnet | Suitability of airborne laser scanning for the assessment of forest protection effect against rockfall[END_REF] in order to map the protective effect of each forest identified and be able to ensure a sustainable service of protection in the most efficient forest ecosystems.

Finally, the results of this study could be used as a means of communication and mediation with foresters, risk managers, local stakeholders in order to raise awareness about this important ecosystem service in mountain areas and to find innovative and consensual ways to sustainably manage these forests in a changing world.

Conclusion

This paper proposes an original method based on opensource data to map forests with a potential protection function against rockfalls for the entire European Alpine region. The reproducible approach developed in this work lies on an objective and robust calibration and validation process as well as a high number of 3-D simulations of rockfall trajectories over the Alps. The results of the harmonized mapping allow a precise identification of the forests endorsing a protection function against rockfall and reflect the important role that mountain forests play in protecting people and human assets from natural hazards.

Figure 1 :

 1 Figure 1: Map of the Alpine Space area and Alpine convention perimeter (red line).

Figure 3 :

 3 Figure 3: Chart of assignment of soil classes.

Figure 4 :

 4 Figure4: Distributions of the energy line angle β of observed rockfall events and modeled rockfall propagations. Modeled distributions are split according to the slope of the release areas and between calibration (6 areas described in Table.2 in dotted lines) and validation sets (entire Alpine Space region in full line). Cumulative probability curves (a.) and comparisons of density curves between observed and modeled β distributions for release areas in slopes ≥ 42°(b.) and slopes in[28-42°[ (c.) are displayed.

Figure 5 :

 5 Figure5: Results of the spatial intersection between modeled propagation areas and location of real rockfall events (a.). Difference is made between propagation areas according to the slope threshold taken into account to identify the release zones. When outside the propagation area, difference is made between the events within a distance to the nearest modeled propagation area (D) less than or equal to the size of the working raster resolution (Res equal to 25 m) and those located beyond. An example of intersection in Isère (France) is also given (b.).

Figure 6 :

 6 Figure 6: Location of the protection forests against rockfalls in the European Alps. Difference is made between protection forests identified with rockfall simulations from release areas located on slopes ≥ 42°and those on slope in the range [28-42°[.

Fig. 7 .

 7 The highest rates are observed in NUTS3 located in the Central Alps (maximum for the Uri Canton in Switzerland with 61% of protection forest considering all release areas). The proportion decreases slowly from the Central Alps to the Pre-Alps. No protection forest was identified for some NUTS3 located on flat areas such as the Po Valley in Italy or the plateau of South Germany.

Figure 7 :

 7 Figure 7: Maps of the proportion of protection forests against rockfalls according to the region considered (Alpine Space: left column -Alpine Convention: right column) and to the slope of the release area (all slopes considered: top row -Only release areas on slope ≥ 42°: bottom row).Proportions are calculated at the NUTS3 levels for both regions.

Figure 8 :

 8 Figure 8: Examples of situation where real rockfall events are outside modeled propagation areas due to a non-detection of release areas in the South of France.

Table 1 :

 1 Description of the human assets used for mapping protection forests.

						OpenStreetMap
	INPUTS	EU-DEM Slope	Corine Land Cover	Seas, lakes & rivers	Human Assets
		Release area		Forest area	Soil classes
		Rockfall propagation simulation		Save rockfall outputs		Rock trajectory intersects human asset	yes
	PROCESS	next simulation 1	outputs Save protection forest	Identify forests on the trajectory area and human assets between release no
						1 n simulations for each release area
		Rockfall hazard outputs:	Protection forests outputs:
	OUTPUTS	Number of passing rocks (PR) Number of stopped rocks (SR) Maximum energy line angle (β max )	Protection Forest Area (PF) Protected human assets (PHA) Minimum length of forest (LMIN)
		Figure 2: General workflow used for the harmonized mapping of protection forests.
	Human asset type Description			Shape	Code
	Railways	National or regional railway lines	Line	1
	Buildings	Houses, socio-economic complexes, infrastructures Polygon 10
	Main roads

Table 2 :

 2 25 m DEM-based slope Angle obtained on potential rockfall sources identified from Lidar data on the different calibration areas. Locations of the different areas are displayed in Fig. 1.

	Zone	Area [km 2 ]	EU-DEM Slope angle (25m) [°] 25th percentile 75th percentile
	AT -Innsbruck	2096	30	43
	CH -Valais	6117	29	42
	IT -Valle d'Aosta	3262	29	43
	FR -Savoie (Bauges)	1614	27	40
	FR -Bas-Rhin (Vosges) 1364	19	31
	SI -Gorenjska	2168	23	39
	All zones	16621 28	42

Table 3 :

 3 Results of the statistical tests used to compare the energy line angle distributions. Observed distribution concerned the β angle of the 2812 rockfall events in the RockTheAlps database. The calibration set corresponds to the β angle distribution modeled on the 6 calibration areas (see Table.2). The validation set concerns the β angle distribution modeled on the entire Alpine Space area. Difference is made between rockfall simulations from release areas located on slopes ≥ 42°(suitable to the Alpine region) and those on slopes in the range[28-42°[.

		Release areas: slopes ≥ 42°Release areas: slopes in [28-42°[
		Calibration set Validation set	Calibration set Validation set
	Wilcoxon p-value*	0.48	0.99	<10 -16	<10 -16
	Bhattacharyya coefficient 0.98	0.96	0.85	0.81
	*Null hypothesis: True location shift is not equal to zero

Table 4 :

 4 Optimal values of R n and R g for each soil class obtained during the calibration process.

	Soil class Description	R n range[-] R g range [m]
	I	Agricultural area	0.23 -0.31 0.05 -0.18
	II	Forest area	0.30 -0.42 0.05 -0.23
	III	Mineral area	0.39 -0.58 0.05 -0.13
	IV	Urban area	0.32 -0.39 0.10 -0.58
	V	Permanent water area	0	0
	VI	Temporary water area	0.30 -0.42 0.40 -0.48
	VII	Screes	0.39 -0.58 0.15 -0.33
	VIII	Meadows and grassland 0.30 -0.42 0.05 -0.18

Table 5 :

 5 Results of the protection forest mapping at the country level for the entire Alpine Space region (a.) and the Alpine Convention region (b.). Difference is made between rockfall simulations from release areas located on slopes ≥ 42°(suitable to the Alpine region) and those on slope in the range[28-42°[.

	a. Alpine Space region
	Country	Total area [km²]	Release area density 28 -42°≥ 42°[%] Asset density	Forest area [km²]	Proportion of protection forest 28 -42°≥ 42°Total
	Austria	83 950	14.9%	2.5%	5.3%	37 280 (44.4%)	8.6%	3.7%	12.3%
	Switzerland	41 290	19.3%	5.2%	9.8%	12 650 (30.6%)	13.6%	12.5% 26.1%
	Germany	46 100	2.6%	0.4%	9.8%	16 410 (35.6%)	2.6%	0.7%	3.3%
	France	101 440	8.6%	1.5%	10.2%	42 765 (42.2%)	6.0%	2.6%	8.6%
	Italy	97 020	16.8%	3.6%	13.7%	35 165 (36.2%)	15.7%	9.0%	24.7%
	Liechtenstein 160	25.0%	6.3%	13.9%	75 (46.9%)	16.0%	11.7% 27.8%
	Slovenia	20 270	7.7%	1.2%	5.8%	11 810 (58.3%)	6.5%	1.9%	8.4%
	Total	390 230	12.4%	2.5%	9.7%	156 155 (40.0%)	9.1%	4.9%	14.0%
	b. Alpine Convention region
	Country	Total area [km²]	Release area density 28 -42°≥ 42°[%] Asset density	Forest area [km²]	Proportion of protection forest 28 -42°≥ 42°Total
	Austria	54 630	22.9%	3.9%	3.7%	28 830 (52.8%)	10.9%	4.7%	15.6%
	Switzerland	25 230	31.6%	8.5%	5.5%	7 610 (30.2%)	20.2%	20.4% 40.5%
	Germany	11 150	10.9%	1.7%	3.0%	4 655 (41.7%)	4.8%	2.2%	7.1%
	France	40 785	21.3%	3.7%	6.5%	19 365 (47.5%)	10.6%	5.4%	16.0%
	Italy	52 030	31.3%	6.7%	11.9%	28 460 (54.7%)	18.2%	11.0% 29.2%
	Liechtenstein 160	26.7%	6.5%	13.9%	75 (46.9%)	16.0%	11.7% 27.8%
	Slovenia	6 770	23.2%	3.7%	3.4%	4 630 (68.4%)	11.7%	4.2%	15.9%
	Total	190 755	25.3%	5.1%	6.8%	93 625 (49.1%)	13.5%	7.9%	21.5%
	3.2.2. Results at the Alpine Space NUTS3 level
	Detailed results at the NUTS3 level are provided in

Table 6 :

 6 Proportion of the human assets potentially affected by rockfall hazard for the entire Alpine Space region (a.) and the Alpine Convention region (b.). The proportion of protection forests that may mitigated rockfall hazard for each asset type is also returned. Difference is made between rockfall simulations from release areas located on slopes ≥ 42°(suitable to the Alpine region) and those on slope in the range[28-42°[.

	a. Alpine Space region
	Human asset type	Proportion of endangered asset Proportion of the protection forests 28 -42°42°Total
	Railways	4.2%	2.3% 6.4%	5.8%
	Buildings	2.4%	1.0% 3.4%	54.9%
	Roads	4.3%	1.9% 6.2%	78.7%
	b. Alpine Convention region
	Human asset type	Proportion of endangered asset Proportion of the protection forests 28 -42°42°Total
	Railways	10.6%	7.0% 17.6%	5.6%
	Buildings	7.2%	3.3% 10.6%	56.4%
	Roads	11.3%	5.8% 17.1%	78.2%

www.alpine-space.eu/projects/rockthealps/en/home

https://land.copernicus.eu/

https://www.openstreetmap.org/copyright/en

https://land.copernicus.eu/imagery-in-situ/eu-dem

https://land.copernicus.eu/pan-european/corine-land-cover

Funding

This work was supported by the project RockTheAlps [grant n°ASP462] from the European Union's Interreg Alpine Space Programme.


Pierre Auger

Olivier Pironneau

Multi-Site Fishing Models

Keywords: Fishery modelling, Marine Protected Area, Partial differential equations

come

Introduction

The logistic model introduced by Verhulst [START_REF] Verhulst | Notice sur la loi que la population poursuit dans son accroissement[END_REF] for population dynamics and controlling the fish biomass B(t) has been extended by Schaefer [START_REF] Schaefer | Some aspects of the dynamics of populations important to the management of commercial marine fisheries[END_REF] to fisheries with an effort function E(t) -conveniently representing the number of fishing boats at sea-and a catchability coefficient q. In [START_REF] Brochier | Can overexploited fisheries recover by self-organization? Reallocation of the fishing effort as an emergent form of governance[END_REF] an extension relating the fishing effort to the market price p of fish is analysed. It applies to a coastal area with a single fishing zone; the fish biomass is a function of the fish reproduction rate r (which includes natural death) and the fish mortality due to fishing; in absence of fishing, there is a natural limit to B(t ∞ ) which we denote by k (the carrying capacity of the site). Mortality due to fishing is proportional to fish biomass and fishing effort:

1 B dB dt = r(1 - B k ) -qE (1)
The fishing effort is a function of the price of fish at time t, p(t):

1 E dE dt = pqB -c. (2)
This says that the rate of the fishing effort is proportional to the difference between the revenue and the cost c of fishing per unit of fishing effort. When the market is liquid the price adjusts daily to balance supply qBE and demand D(p). The demand is a monotone decreasing function of price, for example D(p) = A 1+βp . Hence p = max{p min , (A -qBE)/(βqBE)}; it is necessary to introduce p min beyond which selling of fish cannot happen, but it may be convenient to take p min = 0. This model has 2 nonlinear ODE; the long limit solution can either be B(t ∞ ) = 0 or B ∞ > 0 depending on the parameters and the initial conditions B 0 , E 0 .

Extension of this model to multiple adjacent sites have been proposed in [START_REF] Brochier | Can overexploited fisheries recover by self-organization? Reallocation of the fishing effort as an emergent form of governance[END_REF] among others. It uses the following notations Notations.

• B i (t) the biomass at site i at time t.

• E i (t) the fishing effort at site i at time t.

• β i j is the rate of migration of boats from site i to site j. Assume β i j = β j i . • k i the carrying capacity of zone i.

• m i j is the rate of migration of fish from site i to site j. Assume : m i j = m j i . • q is the catchability.

• , φ time scales for population dynamics of fish and price adjustment respectively.

• r is the fish population growth rate.

• c is the cost per unit fishing effort.

• p market price of fish.

• A, β coefficients in the demand curve D(p).

• I i = 1 if i is a fishing site and 0 otherwise.

The coast is divided in adjacent sites where fishing is allowed or forbidden. I i is defined accordingly. The generalisation of the previous model, proposed in [START_REF] Brochier | Can overexploited fisheries recover by self-organization? Reallocation of the fishing effort as an emergent form of governance[END_REF], is

dB i dt -rB i (1 - B i k i -I i q r E i ) = m i i-1 B i-1 + m i i+1 B i+1 -(m i i-1 + m i i+1 )B i , dE i dt -(pqB i -c i )E i = β i i-1 E i-1 + β i i+1 E i+1 -(β i i-1 + β i i+1 )E i , dp dt = ( A 1 + βp - i q i I i B i E i )φ, i = 1.., I. (3)
In this article we look at these models from 4 angles:

1. What are the partial differential equations, limit of the model when the number of sites tends to infinity? 2. Can they assert the efficiency of marine protected zones? 3. Can it be generalised to two dimensions, i.e. not just the coast? 4. Conversely, starting from the PDE can we recover the above discrete model?

We will take these four questions one by one and simulate numerically all the derived models and make interesting conclusions at the end.

A one dimensional continuous model for coastal fishing

Given a small mesh size δx, for any function x → y(x) let y i denote the value at x i = iδx of y(x). Let us interpret m j i as m(jδx, iδx) and similarly for β. The system above can be rewritten as

dB i dt -rB i (1 - B i k i -I i q r E i ) = (m i B) i+1 + (m i B) i-1 -2(m i B) i -B i (m i i+1 + m i i-1 -2m i i ) = δ xx (m(x)B) -Bδ xx m(x) = δ x (mδ x B) + δ x Bδ x m (4) dE i dt -(pqB i -c i )E i = δ xx (β i E) -Eδ xx β i = δ x (βδ x E) + δ x Eδ x β (5) dp dt = ( A 1 + βp - i q i I i B i E i )φ (6)
where δ x and δ xx are finite difference operators:

(δ x y) i = y i+1 -y i-1 2δx , (δ xx y) i = y i+1 -2y i + y i-1 (δx) 2 .
Let δx → 0, then the system tends to

∂ t B -∂ x (m∂ x B) = ∂ x B∂ x m + rB 1 - B k -I(x) q r E ∂ t E -∂ xx (β(x)E) = ∂ x E∂ x β + (pqB -c)E, for all x such that I(x) = 1 ∂ t p = A 1 + βp - b a I(x)qBEdx φ
If the price adjusts instantaneously to the market then

p = 1 β A b a I(x)qnEdx -1
and, implicitly applying ∂ x (βE) = 0 at interfaces of fishing to non fishing zones. We summarize the results in the following proposition.

Proposition 1. When δx → 0 and fish price adjusts instantaneously, the fish biomass B and the fishing effort E are solution of the following integro-partial-differential system:

∂ t B -∂ x (m∂ x B) = ∂ x B∂ x m + B r 1 - B k -I q r E ∂ t E -∂ x (I∂ x (βE)) = I(x) ∂ x E∂ x β + E Bq β A b a IqBEdx -1 -c (7)
where I(x) is one if x is in a fishing zone and zero if x is in a marine protected area.

2.1. Stationary solution when m and β are small

B k + I q r E = 1, I Bq β A b a IBqEdx -1 = cI or equivalently BqI = cIβ λ A -λ , IE = r q (1 - B k ), λ = b a BqIE = b a r(B - B 2 k )
There is an autonomous equation for λ:

λ = b a I rcβ q λ A -λ (1 - cβ kq λ A -λ ) = b a I rcβ q λ A -λ (1 + cβ kq - cβA kq 1 A -λ ) ⇒ 1 = 1 A -λ b a I rcβ q (1 + cβ kq ) - 1 (A -λ) 2 b a I rA k ( cβ q ) 2 = βY + β 2 X A -λ - Aβ 2 X (A -λ) 2 with X = b a I rc 2 kq 2 and Y = b a I rc q . (8) So, equation a z 2 -b z + 1 = 0 gives z = 1/(A -λ) with b = βY + β 2 X and a = Aβ 2 X. If the discriminant b 2 -4a > 0, i.e. (Y + βX) 2 > 4AX, there is a positive root, z = 1 2a (b + √ b 2 -4a ), λ = A - 1 z , B ∞ = cβ q λ A -λ , E ∞ = r q (1 - B ∞ k ). If (7) is stable, its limit when t → ∞ is B ∞ , E ∞ . If B ∞ = 0, the fish disappear; it is possible only if A = cβ(b -a)E ∞ and E ∞ q = r, i.e. Aq = cβ(b -a)r.

Numerical Simulations

A semi-implicite finite difference scheme is used:

1 δt (B j+1 -B j ) -∂ x (m∂ x B j+1 ) + B j+1 r( B j k + I q r E j ) = ∂ x B j ∂ x m + B j r 1 δt (E j+1 -E j ) -∂ x I∂ x (βE j+1 ) + E j+1 I B j q β + c = I(x) ∂ x E j ∂ x β + E j AB j q β b a IqB j E j dx . ( 9
)
Then this time-discrete system is discretized in space by a finite element method of degree 1.

Numerical Example

When all parameters are equal to 1 and L = b a I > 1 the second degree equation analysed above has two roots:

X = Y = L, (Y + βX) 2 = 4L 2 > 4AX = 4L, b = L + L 2 , a = L, b + √ b 2 -4a = L + L 2 + (L + L 2 ) 2 -4L > 2a A -1 = 2L.
In the following simulation the coast has length L = 3; it is divided into a fishing zone

(0, L 2
) and a marine protected zone ( L 2 , L). The other parameters are m = 0.1 and β = 0.1. At time t = 0 we assume that B = 2B ∞ and E = E ∞ /2. At x = 0 we enforce n = 2B ∞ and at x = L, E = 0.01. The results are shown on Figures 1 and2, with a catchability is q = 1 and a fish and boat migration m = 0.1, β = 0.1. The biomass is decreasing in the middle of the fishing zone but reach an equilibrium in the no-fishing zone. To study the influence of the fish migration we ran another case with fish migration into the no-fishing zone m = 0.1 + 0.1x 2 ; Results are shown on Figure 3 and Figure 4. The migration stabilises the biomass because it diminishes the catch in the fishing zone; yet this diminution does not affect E very much.

To study the influence of boat migration we ran another case with β = 0.1+0.1x 2 ; Results are shown on Figure 5 and Figure 6.

To study the influence of the catchability we ran a case with q = 1 + x 2 . Results are shown on Figure 7 and Figure 8. These simulations show mild influence of migration coefficients and strong influence of the catchability. Futhermore while B is affected most, E doe not change much, indicating that the coupling between B and E is not strong enough.

Extension to two dimensional domains: fishing in the open sea

Consider a domain Ω with a subdomain O where fishing is forbidden. As before let I be the fishing zone indicator. By analogy with the one dimensional case, it is reasonable to consider , for all x ∈ Ω and all T > 0:

∂ t B -∇ • (m∇B) = ∇B • ∇m + B r 1 - B k -I q r E ∂ t E -∇ • (β∇E) = I(x) ∇E • ∇β + E Bq β A b a IqBEdx -1 -c (10)

A numerical test

The domain is a rectangle (-L, L)×(0, H) with a rectangular no fishing zone (-L 1 , L 1 )× (H 1 , H). Initially the number of fish is 2B ∞ on the left (i.e. x < 0) and B ∞ on the right. The number of boats is E ∞ on the left and 2E ∞ on the right. On the left vertical boundary we assume that n = 2B ∞ , i.e. many fish. On the right vertical boundary n = B ∞ . On the top boundary of the no fishing rectangle we impose E = 0.01 and for the second test m = 2B ∞ and nothing for the first test. All parameters are equal to 1, except β = 0.1 and m = 0.1 for the first test and m = 0.5 for the second. The system is simulated up to time t = 50 with 50 time steps. The mesh has 1483 vertices. and 2884 triangles. Results are shown on Figures 9 and10. It shows clearly how fishermen migrate from the right (where B is low initially) to the left and making the left zone harder to fish in. In the end there are no fish in the protected zone as well when the fish migration coefficient is small: m = 0.1; when m = 0.5 the left zone and the no fishing zone at the center are no longer empty at the end of the computation .

A continuous model for Plankton, Fish and Fishermen

The biggest difficulty of the previous model is the mobility function m which is probably very hard to guess in practice. In this section we will see that if we can give a more solid ground to the previous two dimensional model. To this end we need to introduce plankton as an incentive for the motion of the fish. The notations are as follow:

Notations Variables

• P (x, t) plankton density at point x and time t.

• B(x, t) fish population at x and time t.

• E i (t) fishing effort by fishermen i at time t

• X i (t) position of fisherman i at time t.

• π(x, t) probability density of fishing boats at x at time t.

• p(t) price of fish at time t.

Parameters

• Ω the domain occupied by the sea. (0, T ), the time interval for the study.

• time scale of the dynamics of plankton

• v(x, t) sea current at point x and time t.

• µ, ν, κ diffusion constant of plankton, fish position, boat position

• a reproduction rate of plankton

• b rate at which the plankton is eaten by fish

• α, β, γ, K = α T 2 β+γ behavioral parameters of a fisherman • r rate of reproduction of fish • d rate of death of fish
• q catchability of fish

• c cost of operation of a boat

• A, B parameters of the demande for fish for a given price.

Plankton

The plankton density P is regenerated at rate aP and eaten by the fish at rate b times B, the fish density. The logistic equation for P is:

∂ t P + v • ∇P -µ∆P = aP (1 -P ) -bBP, P (0) = P 0 in Ω × (0, T ) (11)
where P 0 (x) is the probability density function of plankton at initial time. The model assumes that in absence of fish the long time limit of P is one. Here v is the sea current velocity. The term -aP 2 models the natural death of the plankton. Having the rate of change of plankton density proportional to the concentration is convenient because of the following remark. More realistic models can be found in [START_REF] Peter | NPZ models of Planton dynamics[END_REF].

Remark 1. The right hand side is also aP (1 -b a B -P ); if B < a b , then P is positive and bounded by 1 -b a B if it is so at time zero. Otherwise P may become negative and the model is no longer meaningful. Remark 2. When c := 1 -b a B is constant and v = 0 and µ = 0, the solution of Ṗ = aP (c -P ) is P = c e ac t /(1 + e ac t ), tends to c when t → ∞. When µ > 0 and Ω is bounded there is no non-zero limit t → ∞ of ∂ t P -µ∆P = aP (c -P ), P (0) given in (0, c ) and P | ∂Ω = 0.

(12)

Modelling of a fishing boat

A single fisherman will want to maximize its catch while keeping his cost minimum. Its catch is proportional to the fish population B(x, t) and its cost is proportional to the square of the modulus of the velocity of the boat to operate the engine; he may also prefer not to go too faraway so X( T ) has a cost. So over a time interval (0, T ), max

U J(U ) := T 0 αB(X(t), t)dt - T 0 β 2 |U | 2 dt - γ 2 |X( T )| 2 : Ẋ = U (X(t)), X(0) = 0,
where α is a function of the price of the fish and β, γ include the price of fuel; they are constant over (0, T ) but vary from one day to the next. The solution is found by a calculus of variations. Let δU be a deviation from the optimal solution, then δ Ẋ = δU, δX(0) = 0,

δJ = T 0 [α∇BδX -βU • δU ] -γX( T ) • δX( T ) = T 0 α∇Bdτ -γX( T ) • δX( T ) - T 0 t 0 α∇Bdτ δ Ẋ + βU • δU (13)
We will have δJ = 0 if

t 0 α∇Bdτ + βU = T 0 α∇Bdτ -γX( T )
which implies

Ẍ(t) = - α 2β ∇B(X(t), t), X(0) = 0, β Ẋ( T ) + γX( T ) = 0. U = α β T t ∇Bdτ + αγ β + γ T 0 t 0 ∇B(X(τ ), τ )dτ -∇B(X(t), t) dt (14)
As a first approximation, we assume that ∇B is constant; it gives

U = α T 2 β + γ ∇B Interpretation.
In this formula T is the length of time over which the fisherman takes a decision; it is not the whole day because only a rough knowledge of the local value of ∇B near his boat is available to the fisherman. We retain from the above analysis that the velocity of the boat when it moves is proportional to ∇B. The coefficient of proportionality is a consequence of the psychology of the fisherman; it can only be measured.

Fish

Fish move so as to optimize their food intake for a given energy expenditure. The problem is similar to that of the fisherman who, in the previous paragraph optimized his catch for a given effort level. Hence by the same optimisation argument the optimal swimming velocity of the fish is proportional to the gradient of the plankton density.

Fish die naturally at rate d and reproduce at a rate r proportional, function of the food supply P . Fishing depletes the fish population; let N m (t) be the number of fish caught during the time interval (0, T ). Let n m (t) = Ṅm (t). Then the fish population B is driven by

∂ t B + ∇ • (B∇P ) -ν∆B = r(P )B -Bd - M 1 n m (t)δ Xm(t) B(•, 0) = B 0 (•) given, ( 15
)
where δ x denote the Dirac mass at x. In the sequel we assume that r(P ) = rP .

Modelling of the number of fish caught by a fisherman

In an ideal situation one would state that n m (t) = qB(X m (t), t), the coefficient of proportionality q being the catchability. However the fishermen have an operating cost c which must be kept less than the profit pqB where p is the market price of fish.

In [START_REF] Brochier | Can overexploited fisheries recover by self-organization? Reallocation of the fishing effort as an emergent form of governance[END_REF], n m (t) = B(X m (t), t)qE m (t), i.e. the coefficient of proportionality between n m and B(X m (t), t) has a new coefficient, the "fishing effort" E m of the fisherman m; it is given by

1 E dE dt = pqB -c (16)
The fishing effort is defined on the trajectory of the fisherman's boat, given by Ẋm = U = K∇B(X m (t), t), X(0) given.

The price of fish is given by balancing demand A 1 + Bp and supply qEB where EB =

∂ t B + ∇ • (B∇P ) -ν∆B = r(P )B -Bd - M 1 qE m (t)B(X m (t), t) B(•, 0) = B 0 (•) (18)
with (17) and

Ėm (t) = E m (t) max{0, A M 1 E m B(X m , t) -q B(X m (t), t) B -c }. ( 19
)
Proposition 2. A continuous model for P, B, {E m } M 1 , the plankton density, the fish biomass and the effort function of each fisherman, is given by (11),(18),(17), (19).

A large number of fishing boats

Let us assume that at t = 0 we know the distribution of fishing boats at sea. Let us add some randomness, with a 2D Brownian motion W , E[W 1 W 2 ] = 0, to the trajectories of each boat:

dX t = K∇B(X t , t)dt + κdW t . ( 20
)
The PDF of this process is given by the Kolmogorov equation

∂ t π + ∇ • (πK∇B) -κ∆π = 0, π(0) given in R 2 . ( 21
) So V := E[EB] = R 2
EBπdx with E given by ( 16) with p = 1 B ( A qV -1) and B given by

∂ t B + ∇ • (B∇P ) -ν∆B = r(P )B -Bd -qEB B(•, 0) = B 0 (•) (22)
Note that this is an approximation because it is not true that M 1 B(X m (t), t)qE m (t)δ Xm(t) converges to qEB when X m is random and M → ∞. Definition 1. When the number of fishing boats M is large, a continuous model for P, B, E, π, the plankton density, the fish biomass the global effort function and the PDF of the boat density, is given by ( 11),( 22),( 16), (21). This leads to an equation for B∞ ; with B = B/q and b = b/a it is

rb B2 ∞ -(r -d -B cb ) B∞ -(r -d)B c + A = 0. ( 23
)
With a marine protected zone. The interesting case is when parameters are such that the above equation has no solution (fish going extinct); in the no fishing zone the equations for P and B give an stationary average B∞ = (1 -d r ) a b , P∞ = d r .

4.8. Numerical Simulations 4.8.1. Plankton-Fish and no fishing boats

∂ t P -µ∆P = aP (1 -P ) -bBP, P (0) = P 0 (x) given, positive, bounded. ∂ t B + ∇ • (B∇P ) -ν∆B = r(P )B -dB B(•, 0) = B 0 (•) given. ( 24
)
The long time limit is B∞ = (1 18),( 20),( 16) with K = 80, q = 20, A = 10, B = 1, c = 1; P 0 and B 0 as above and 50 boats randomly distributed with a Gaussian distribution centered at x = -1 on the coast. Three Simulations are performed. One with fishing restriction in a rectangle in the middle of the domain, not touching its base, in which it is seen that the fish population does not go to zero near the fishing zone. In the second simulation fishing is allowed everywhere. Then it is seen that the fish population tends rapidly to zero everywhere. In the third simulation d = 2, forcing P∞ = 0. Even with a no fishing zone the fish population goes to zero everywhere. The fish do as above and a strip swims to the left; the fishermen hunt them and fish almost all the fish nearer to the coast. Finally figure 13 shows the fish price variation for the first 2 cases.

-d r ) a b , P∞ = d r . Hence it is important that d < r. The domain is a rectangle Ω = (-L, L) × (0, H) with L = 3, H = 2. Define II = Ω\(-L 1 , L 1 ) × (H 1 , H) with L 1 = 1, H 1 = 0.

Plankton -Fish -Fishing Effort -Fishermen density

We assume that there is a density of fishermen spread uniformly on the coast initially and spreading according to (21). As above when fishing is allowed everywhere, the fish density tends to zero everywhere. If fishing is not allowed in a rectangle then the fish density does not go to zero in the rectangle. The long time limit is B∞ = 0 because (23) is B2

∞ -79 80 B∞ + 9.925 = 0 which has no real root. Results are shown on figure 14. The conclusion of this section is that the two models are capable of asserting the benefit of a marine protected area. The model of Proposition 2 is more rigorous but harder to simulate. When the number of fishermen is large enough one may use the model defined in Definition 1, with some caution. Figure 14: Left Fish density at T=2, when fishing is allowed everywhere. Right, fishing is not allowed in the rectangle

Very coarse discretization

In this section we take the opposite path: going from the continuous system to a discrete system with a few fishing or no fishing zones.

Three zones

To connect this analysis with the discrete systems proposed earlier in [START_REF] Brochier | Can overexploited fisheries recover by self-organization? Reallocation of the fishing effort as an emergent form of governance[END_REF]; we consider 2 zones and an interface or, equivalently, 3 zones: fishing is allowed on the left, but may be forbidden on the right. The PDEs are discretized using finite differences on a very coarse mesh made of 5 × 3 points (see figure 15). On the boundary of Ω = (-L, L) × (0, H) normal derivatives of P and B are set to zero; referring to figure 15 this means that all quantities at vertices q n , n = 8, 7, 6 and q n , n = 11, 12, 13 are equal to the same quantities at q n , n = 1, 2, 3 respectively. Similarly quantities at q n , n = 0, 4 are equal to the same quantities at q n , n = 1, 3 respectively. Quantites at q n , n = 5, 9, 10, 14 do not appear in the discretisation. It implies that there are only 3 points of interest, of coordinates

q 1 = {- L 2 , H 2 }, q 2 = {0, H 2 }, q 3 = { L 2 , H 2 }.
q 0 q 1 q 2 q 3 q 4 q 9 q 8 q 7 q 6 q 5

q 10 q 11 q 12 q 13 q 14

Figure 15: Finite difference grid for the discretization of the domain Ω which can be viewed either as 2 zones and an interface in the center or 3 zones, left, right and the zone in red. On the right of the vertical arrowed line in the center we may impose no fishing.

The unknowns are {P n , B n , E n } n=1,2,3 . Let h = L/4 be the mesh size, then a finite element method of degree one gives

Ṗ1 - µ h 2 (P 2 -P 1 ) = aP 1 (1 -P 1 ) -bB 1 P 1 Ṗ2 - µ h 2 (P 3 + P 1 -2P 2 ) = aP 2 (1 -P 2 ) -bB 2 P 2 Ṗ3 - µ h 2 (P 2 -P 3 ) = aP 3 (1 -P 3 ) -bB 3 P 3 Ḃ1 - µ h 2 (B 2 -B 1 ) - K 2h [(B 2 + B 1 )(P 2 -P 1 )] = rB 1 P 1 -dB 1 -qB 1 E 1 Ḃ2 - µ h 2 (B 3 + B 1 -2B 2 ) + K 2h [(B 2 + B 1 )(P 2 -P 1 ) -(B 2 + B 3 )(P 3 -P 2 )] = rB 2 P 2 -dB 2 -qB 2 E 2 Ḃ3 - µ h 2 (B 2 -B 3 ) + K 2h [(B 2 + B 3 )(P 3 -P 2 )] = rB 3 P 3 -dB 3 -qB 3 E 3 Ėn = E n max 0, A E 1 B 1 + E 2 B 2 + E 3 B 3 -q B n B -c , n = 1, 2, 3. ( 25
)
Additionally P n (0), B n (0), E n (0), n = 1, 2, 3 must be given. Note that if E i (0) = 0 then E i (t) = 0 at all times t > 0. Hence a no-fishing zone is specified simply by putting E n (0) = 0 in the zone. w1,w2,w3 = 0.5,0.5,0.0 phi1,phi2,phi3 = 1,1,1 E1,E2,E3 = 0.01,0.,0.

Results: see Figure 16. The results are shown on figure 16 and show fish depletion on Figure 16: Evolution in time of the fish population at points q 1 , q 2 , q 3 , fishing is allowed everywhere the left where fishing occurs, and non zero limit of fish biomass in the center and on the right.

A Marine Protected Area surrounded by 2 fishing zones

The previous simulation has shown that a no-fishing zone surrounded by 2 fishing zones seems capable of keeping the fish biomass non-zero in the MPA. In this paragraph we further validate the concept with an asymptotic analysis. We consider a large fishing area in the middle of which a protected area where fishing is prohibited has been installed. We therefore consider a central zone 2 which is a Marine Protected Area (MPA) surrounded on both sides by fishing zones respectively 1 and 3.

The model takes into account the dynamics of the plankton population as well as that of the fish that feed on it and are caught. Let P 1 (t), P 2 (t) and P 3 (t) be the plankton biomass in zone 1, 2 and 3 at time t. Similarly, let B 1 (t), B 2 (t) and B 3 (t) be the fish biomass in zone 1, 2 and 3 at time t in Megatons. Let E(t) be the total fishing effort at time t. r is the growth rate of the plankton species and K its carrying capacity assumed to be the same in the three zones. Fish feed on plankton by predation. a is a predation rate. e is the conversion coefficient from plankton biomass to fish biomass. d is the fish mortality rate. Fishing boats catch fish. q is the catchability. The complete model reads:

dP 1 dτ = (kP 2 -kP 1 ) + rP 1 (1 - P 1 K ) -aP 1 B 1 (26) dP 2 dτ = (kP 1 + kP 3 -2kP 2 ) + rP 2 (1 - P 2 K ) -aP 2 B 2 (27)
dP 3 dτ = (kP 2 -kP 3 ) + rP 3 (1 - P 3 K ) -aP 3 B 3 (28) dB 1 dτ = (mB 2 -mB 1 ) + eaP 1 B 1 -dB 1 -qB 1 E (29) dB 2 dτ = (mB 1 + mB 3 -2mB 2 ) + eaP 2 B 2 -dB 2 (30)
dB 3 dτ = (mB 2 -mB 3 ) + eaP 3 B 3 -dB 3 -qB 3 E (31) dE dτ = A β -cE - q β B 1 E - q β B 3 E ( 32
)
Where τ is the fast time, t = ετ is the slow time and where ε << 1 is a small dimensionless parameter.

Remark 3. This model is identical to (25) except that notations have changed and the transfer rate K is zero.

Hence this complete model assumes that the displacement of the plankton and the fish between the zones is by diffusion. We assume that the three areas are large and that the migratory flows between the zones remain small. Thus, the diffusion of plankton and fish takes place on a slow time scale while population dynamics in each patch is fast in comparison. It is also assumed that the price varies according to the difference between the demand and the supply. The demand D(p) = A (1+βp) is assumed to be a monotone decreasing function of price. The supply is the instantaneous catch in the two fishing areas. We further assume that the price is fixed very rapidly.

We will now investigate whether it is possible to obtain an equilibrium corresponding to an overfishing situation leading to the extinction of the fish stock in the two fishing zones while the fish succeeds in surviving in the marine protected area (MPA ). Therefore, let us look for existence of a fast equilibrium with no fish in fishing zones and fish in MPA, (K, P * 2 , K, 0, B * 2 , 0, E * ). It is simple to check that such an equilibrium exists in the positive quadrant with:

P * 2 = d ea , B * 2 = r a (1 - d eaK ), E * = A cβ (33)
Furthermore it is locally asymptotically stable. Now let's look at the dynamics on the slow timescale. The reduced model is obtained by replacing the previous fast equilibrium in the complete model and using the slow time. The slow model is a diffusion model with plankton and fish exchanges between the zones: However, these flows between the different zones are small. Fish move from the central area to the two fishing zones and are harvested rapidly. As a consequence, there is no fish in the fishing areas and the plankton tends to its carrying capacity. This simple model shows that it is possible to obtain sites with fish and sites without fish in meta-population type models. For this, it is necessary to have a small flow of plankton and fish between the different sites. Cases of fast diffusion or going on at the same time scale as local population dynamics do not allow it.

dP 1 dt = (kP * 2 -kK),

Conclusion

Starting from the multi-site model of [START_REF] Brochier | Can overexploited fisheries recover by self-organization? Reallocation of the fishing effort as an emergent form of governance[END_REF] we have derived its continuous integro-differential limit. Numerical simulations have asserted the marine protected areas as an efficient tool for keeping a minimum fish biomass at sea, at least in simple cases. The continuous model can be extended in 2D by analogy, but without rigorous justification.

In an attempt to build similar 2D models with more mathematical justification we have turned to the optimisation of utility functions for the fish and the fishermen. Fish change place to gather more food and as the food density grows its feeding effort becomes less. It has been shown that the fish velocity is proportional to the gradient of its food density P . This requires the modelling of plankton (or smaller fish), which, in turn, moves with the sea currents. Similarly it is in the interest of the fishermen to move proportionally to the gradient of the fish biomass B. Such behavioural remarks lead to a model for P, B and the fishing effort E m of each fisherman m. The simulations have shown that it is again possible to design marine protected areas which will prevent the disappearance of fish. Unfortunately as the number of fishermen grows one would like to work with the probability density of the fishermen presence at point x at time t, but the consequent continuous model requires more mathematical studies to be asserted rigorously. Finally, in the last section, we have discretised the continuous model on a very coarse mesh in order to derive simpler multi-site 2D models consistent with the continuous models and for which all parameters have a clear meaning. Again they show that it is possible to design a marine protected zone in which the fish population has a non-zero asymptotic limit.

Let us recall that any model needs to be confronted to data for validity, a huge task for us for the future. Then all parameters must be guessed from the data, but we have shown in [START_REF] Auger | Parameter Identification by Statistical Learning of a Stochastic Dynamical System Modelling a Fishery with price variation[END_REF] that the techniques of AI are quite efficient for such a task.

Figure 1 :

 1 Figure1: The fish biomass B(x, t) computed by the continuous model[START_REF] Verhulst | Notice sur la loi que la population poursuit dans son accroissement[END_REF] versus costal position x ∈ (0, 3) at time (depth on the figure) t ∈ (0, 50); x > 1.5 is a no fishing zone. The catchability is q = 1. The fish and the boat migration are m = 0.1, β = 0.1.

Figure 2 :

 2 Figure 2: Fishing effort E near the coast x ∈ (0, 3) at time (depth on the figure) t ∈ (0, 50); x > 1.5 is a no fishing zone. The catchability is q = 1. The fish and the boat migration are m = 0.1, β = 0.1.

Figure 3 :

 3 Figure 3: Number of fish n near the coast x ∈ (0, 3) at time (depth on the figure) t ∈ (0, 50); x > 1.5 is a no fishing zone. The catchability is q = 1. The fish and the boat migration are m = 0.1 + 0.1x 2 , β = 0.1 respecctively.

Figure 4 :

 4 Figure4: Fishing effort E near the coast x ∈ (0, 3) at time (depth on the figure) t ∈ (0, 50); x > 1.5 is a no fishing zone. The catchability is q = 1. The fish and the boat migration are m = 0.1 + 0.1x 2 , β = 0.1 respecctively.

Figure 5 :

 5 Figure 5: Number of fish n near the coast x ∈ (0, 3) at time (depth on the figure) t ∈ (0, 50); x > 1.5 is a no fishing zone. The catchability is q = 1. The fish and the boat migration are m = 0.1, β = 0.1+0.1x 2 .

Figure 6 :

 6 Figure 6: Number of boats E near the coast x ∈ (0, 3) at time (depth on the figure) t ∈ (0, 50); x > 1.5 is a no fishing zone. The catchability is q = 1. The fish and the boat migration are m = 0.1, β = 0.1 + 0.1x 2 .

Figure 7 :Figure 8 :

 78 Figure 7: Number of fish n near the coast x ∈ (0, 3) at time (depth on the figure) t ∈ (0, 50); x > 1.5 is a no fishing zone.The catchability is q = 1 + x 2 . The fish and the boat migration are m = 0.1, β = 0.1.

Figure 9 :

 9 Figure 9: Case 1. Top 5 pictures: Evolution (one picture ever 10 time steps) of the fishing effort. The fish density is not prescribed at the top of the nofishing zone and m = 0.1. The last picture is the fish density at the final time.

Figure 10 :

 10 Figure 10: Case 2. Top 5 pictures: Evolution (one picture ever 10 time steps) of the fishing effort. The fish density is prescribed at the top of the no-fishing zone and m = 0.5. The last picture is the fish density at the final time.

M 1 E

 1 m (t)B(X m (t), t) is the sum of all fishing efforts: p = max{0,

4. 7 .

 7 Long time solution Fishing allowed everywhere. Let us assume that the time scale of the dynamics of P is small: << 1. Then P∞ := lim t→∞ w = 1 -b a B. Let B∞ = lim t→∞ Ω Bdx and similarly for Ē∞ . Integrating the equation of B over Ω and the equation of E gives B∞ (r P∞ -d -q Ē∞ ) = 0, pq B∞ -c = 0, Bp = A q Ē∞ B∞ -1.

  15. Starting from P 0 = 1 x<L 1 and B 0 = 0.01 + 0.51 -L 1 <x<L 1 , we obtained figure 11 at T = 2, when r = 1, a = 2, b = 1, d = 0.5. All diffusions, µ, κ, ν are set to 0.1. Results are shown on figure 11. The computation shows that most fish swim slightly to the left because that is where the plankton is. The long time limit is P∞ = 0.5, B∞ = 1.

Figure

  Figure 11: P

Figure 12 :Figure 13 :

 1213 Figure 12: Top 2 figures: Level lines of the fish population function and the boats positions at T = 2.Top: fishing is allowed everywhere. Below, fishing is not allowed in the rectangle. On the top image all fish have disappeared. On the second image the dynamics in the rectangle is sufficient to keep the fish population from extinction. The third image corresponds to d = 0; even with a no fishing zone fish disappear everywhere. Bottom: color map of the fishing effort without (left) and with (right) a non fishing zone

  = 10,1 # price constants c = 1 # cost for fisherman K = 3 # moving speed of fish to plankton Initial conditions.

pierre.auger@ird.fr, IRD UMI

209, UMMISCO, Sorbonne Université, Bondy, France 2 olivier.pironneau@sorbonne-universite.fr, LJLL, Sorbonne Université, Paris, France.

Acknowledgement

We thank C. Bardos and G. Pagès for their helpful comments; All 2D computations have been done with the public domain FreeFEM++ [4].


Milan Stepanov

Giuseppe Valenzise

Frédéric Dufaux

HYBRID LEARNING-BASED AND HEVC-BASED CODING OF LIGHT FIELDS

Keywords: light field, deep learning, HEVC, JPEG Pleno

Light fields have additional storage requirements compared to conventional image and video signals, and demand therefore an efficient representation. In order to improve coding efficiency, in this work we propose a hybrid coding scheme which combines a learning-based compression approach with a traditional video coding scheme. Their integration offers great gains at low/mid bitrates thanks to the efficient representation of the learning-based approach and is competitive at high bitrates compared to standard tools thanks to the encoding of the residual signal. The proposed approach achieves on average 38% and 31% BD rate saving compared to HEVC and JPEG Pleno transform-based codec, respectively.

INTRODUCTION

Light fields (LFs) capture light intensity of scene objects from various angles. The captured information offers novel applications such as refocusing and a perspective change at the cost of increased dimensionality and thus storage demand. The necessity of efficient compression methods was acknowledged by JPEG committee which started the JPEG Pleno initiative in order to provide a standard framework for the representation and coding of plenoptic data. In particular, some coding tools are proposed and optimized to encode 4D light fields [START_REF]JPEG PLENOlight field coding common test conditions[END_REF].

On the other side, some recent compression approaches based on learning focus on optimizing the whole coding pipeline in an end-to-end fashion [START_REF] Valenzise | Quality assessment of deeplearning-based image compression[END_REF]. Typically, these schemes employ deep auto-encoders to learn good representations of data, which are then quantized and entropy coded. The potential of these methods has been shown by their competitive performance in terms of objective and subjective metrics compared to traditional coding schemes [START_REF]Performance evaluation of learning based image coding solutions and quality metrics[END_REF]. However, it has been found that auto-encoder based approaches provide significant gains at low bitrates, while they generally fail to provide high quality and near-lossless reconstructions at high bitrates. This phenomenon is mainly due to the nature of autoencoders, which are intrinsically lossy. Conversely, traditional codecs are designed to span the full quality range, and in particular to provide near lossless performance at higher bitrates.

In order to incorporate the benefits of both approaches, we aim at overcoming the observed lack of scalability of deep learning approaches by adding an enhancement layer. We propose a hybrid scheme consisting of a base layer which provides high gains at low/mid bitrates and serves as an efficient predictor for high bitrates. This is complemented by an enhancement layer which allows the coding of the residual signal via a traditional coding scheme and provides improved performance at high bitrates. Furthermore, we explore various traditional coding schemes for the residual signal and show that even with a simple approach such as scalar quantization it is possible to achieve significant gains with respect to the base layer and to be competitive with state-of-the-art LF codecs.

RELATED WORK

We divide the related work into two categories: the LF compression using traditional approaches based on conventional prediction and transformation pipelines, and end-to-end learning-based image compression methods.

High Efficiency Video Coding (HEVC) and its extensions significantly influenced LF coding [START_REF] Vieira | Data formats for high efficiency coding of lytro-illum light fields[END_REF][5] [START_REF] Dricot | Integral images compression scheme based on view extraction[END_REF]. Typically, HEVC is used to encode a set of key views which is at the decoder used to reconstruct the complete LF. Zhao et al. [START_REF] Zhao | Light field image coding via linear approximation prior[END_REF] reconstruct the LF using a linear combination of decoded views. Jiang et al. [START_REF] Jiang | Light field compression using depth image based view synthesis[END_REF] recovers missing views via warping and inpainting of occluded regions. Viola et al. [START_REF] Viola | A graph learning approach for light field image compression[END_REF] use graph learning method to learn disparities between views and restore LF. These approaches are based on schemes including different functional blocks and therefore difficult to optimize together. Conversely, an end-to-end scheme learns a single function which jointly optimizes and integrates all needed operations. JPEG Pleno provided a coding tool for the light field modality which works in two modes: the prediction (WaSP) and the transformation (MuLE) [START_REF] Schelkens | JPEG pleno light field coding technologies[END_REF]. The prediction mode uses disparity-based warping and view merging to predict views, while the transformation mode exploits the redundancy in LF blocks using 4D Discrete Cosine Transform. Although, the solutions are tailored for LF compression, they are limited in the same way as the previous approaches.

End-to-end learning-based compression has been recently proposed in [START_REF] Ballé | Variational image compression with a scale hyperprior[END_REF] [START_REF] Toderici | Variable rate image compression with recurrent neural networks[END_REF] [13] [START_REF] Rippel | Real-time adaptive image compression[END_REF] and has gained huge popularity due to its ability to replace the whole traditional compression pipeline with a single function. Ballé et al. proposed an end-to-end compression approach which consists of analysis and synthesis functions corresponding to an encoder and a decoder in conventional pipelines, respectively, plus a uniform quantizer. In addition, they propose a differentiable quantization mechanism allowing to optimize the rate-distortion (RD) function directly [START_REF] Ballé | Variational image compression with a scale hyperprior[END_REF]. Theis et al. [START_REF] Theis | Lossy image compression with compressive autoencoders[END_REF] propose a similar approach but deals with quantization and bitrate estimation in different manner. Rippel et al. [START_REF] Rippel | Real-time adaptive image compression[END_REF] propose a real-time codec which applies a pyramidal analysis for the feature extraction and an adaptive coding module and regularization. Conversely to the previous approaches, Toderici et al. [START_REF] Toderici | Variable rate image compression with recurrent neural networks[END_REF] overcomes the necessity to train a separate model for each lambda value in the RD function by adopting encoding in a progressive manner.

Our approach extends the work of Ballé et al. [START_REF] Ballé | Variational image compression with a scale hyperprior[END_REF], and operates on data of a higher dimension, requiring a careful design of the network architecture to handle a particular filtering across different views. Furthermore, motivated by the limitations of auto-encoders in providing high-quality reconstructions, we also introduce an enhancement layer to encode the residual signal.

PROPOSED METHOD

Base Layer

As a base layer, we propose an end-to-end trained compression scheme based on the recent work of Ballé et al. [START_REF] Ballé | Variational image compression with a scale hyperprior[END_REF] with added modifications to adjust to the LF structure. More specifically, the scheme takes as input a LF image, reshapes it by extracting sub-aperture views and stacking the views along the third dimension following horizontal raster scan order, i.e.

row by row selection. The compression scheme is comprised of three functional blocks: an analysis function f a (θ) which creates a more compact representation of the input y = f a (x), a quantization block or an entropy bottleneck Q(η), which provides quantized version of y, ỹ, and a synthesis function f s (φ) whose goal is to reconstruct the input from the quantized compact representation x = f s (ỹ).

The analysis function comprises a set of sequential nonlinear, downsampling, and convolutional layers while the synthesis function is a symmetric counterpart of the encoding function with downsampling layers replaced by upsampling layers. More specifically, each layer utilizes 2D filters, which slide along spatial dimensions and accumulate the contribution of each feature map in the data signal in order to jointly learn structures in spatial and angular domains. The entropy bottleneck works in two modes depending on the inference phase, i.e. the training or testing phase. During the training phase the entropy bottleneck adds a uniform noise U(-0.5, 0.5) to the transformed representation of the input to approximate quantization in a differentiable manner. Furthermore, the entropy bottleneck learns the probability density function of each feature map of the compacted representation and utilizes them for entropy coding. A detailed description of the network parameters is summarized in Fig. 2. Weights of the analysis function, the entropy bottleneck and the synthesis function, θ, η and φ respectively, are learned by minimizing the RD function J(θ, η, φ; x) = R(ỹ) + λD(x, x) where the rate R(ỹ) is modeled with the entropy of the compressed bottleneck, the distortion D is the mean square error between the input x and the decoded LF x, and the parameter λ governs the trade-off between the rate and the distortion. We trained five models by selecting five different lambda values. For each model, the weights are learned using Adam optimizer with a learning rate of 1 × 10 -4 and 1 × 10 -3 , for θ and φ, and η respectively. The final bitstream is obtained by packing encoded coefficients, the LF size and the lambda parameter.

Enhancement layer

We observe that the auto-encoder based approach used in the base layer reaches a saturation in performance at higher bitrates. A possible solution to this problem is increasing the capacity of the network. However, this requires to use increasing model complexity when approaching higher quality points. Instead, we propose to introduce an enhancement layer to encode the residual signal between the original light field and the reconstruction from the base layer. The advantage of this hybrid approach is that the residual coding allows incorporating any traditional image coding scheme. We compute the residual signal by subtracting an input LF and its prediction obtained using the base scheme followed by cropping to a fixed range of values and translation to positive values.

For the purpose of an enhancement layer, we explore three coding schemes. We evaluate scalar quantization followed by sample entropy estimation due to its simplicity. Furthermore, we evaluate HEVC coder using Intra and Inter prediction to capture the remaining correlation in spatial and angular domains.

RESULTS

Testing conditions

In our test, we follow test conditions proposed by JPEG Pleno (CTC) [START_REF]JPEG PLENOlight field coding common test conditions[END_REF] with slight modifications due to the setup of our approach. Namely, as our approach operates on the luminance component only, we convert LF to YCbCr colour space and set chroma components to neutral colour. We use the EPFL LF image dataset [START_REF] Rerabek | New light field image dataset[END_REF] and select 4 test images as shown in Fig. 3 and divide the rest of the images into training and validation sets in the ratio 80/20. Each raw image is decoded using LFToolbox version 0.4 [16] [17], and a subset of 13 × 13 sub-aperture views, each having the spatial resolution of 625 × 434 pixels is selected to generate a LF image. The proposed approach is compared with HEVC (x265 implementation as defined in the JPEG Pleno CTC) and JPEG Pleno Verification Model 2.0 (VM 2.0) in terms of average PSNR and SSIM. Prior to x265 encoding, a LF is converted to YCbCr 4:4:4 10-bit colour space and arranged in a pseudovideo sequence following serpentine scanning order as defined in CTC. We evaluate VM 2.0 in the transformation mode (MuLE) as it is more efficient for the coding of lenselet data. Four bitrates are used to generate RD curve: 0.75 bpp (bits per pixel), 0.1 bpp, 0.02 bpp and 0.005 bpp.

Joint rate-distortion curve

The contribution of the enhancement layer is evaluated in terms of RD by generating a RD curve that corresponds to the joint coding of the base layer and the enhancement layer. The residual image is encoded using the proposed approaches for different values of the quantization parameter providing local RD curves around each RD point of the base layer. Fig. 4 (a) illustrates the RD curve of the base layer (red) and the points of local RD curves (asterisks). The final RD curve is obtained by computing the convex hull around all the points (cyan). Fig. 5. Performance with respect to PSNR (upper row) and SSIM (bottom row): proposed Hybrid approach using HEVC Inter for the enhancement layer, learning-based base layer, HEVC, MuLE.

The comparison of enhancement layers

Fig. 4 (b) compares the performances of the base layer against the proposed extensions for the LF Bikes. It can be noticed that at low bitrates enhancement layers do not improve performance suggesting the superiority of the base layer. On the contrary, at high bitrates we can notice the benefit of adding the enhancement layer and note that even a simple method such as scalar quantization provides ∼ 5 -15% savings comparing to the base layer. Further improvements are obtained with HEVC's Intra and Inter prediction modes which gain by exploiting the spatial and inter-view correlations in the residual signal. The quantification of the performance is presented in Table 1 for all test contents.

The comparison with respect to anchors

We compare the best performing proposed hybrid approach with the original base layer and two anchor methods proposed by JPEG Pleno, and provide RD comparison in terms of PSNR and SSIM in Fig. 5 and BD rate in terms of PSNR in Table 2. The hybrid approach gains at lower bitrates thanks to the learning-based base layer, while the introduction of the enhancement layer increases the performance and makes the approach competitive with the anchors at high bitrates. Nevertheless, overall performance suggests significant gains of the approach against the two anchors ranging from ∼ 16% to ∼ 50% saving with respect to HEVC and from ∼ 25% to ∼ 40% saving with respect to MuLE.

CONCLUSION

We propose a hybrid coding scheme for LF based on a learning-based approach and HEVC in order to overcome the saturation in performance at high bitrates of auto-encoders. Furthermore, we show that the proposed approach achieves better performances against state-of-the-art anchors. More precisely, the results show that it is possible to greatly improve the performance at high bitrates and moderately at mid bitrates compared to the learning-based base approach. At the same time, the latter provides highly competitive gains at low bitrates, demonstrating the superiority of learning-based representations also for the case of LF.

In future work we aim to optimize the base layer and the enhancement block jointly in an end-to-end fashion and to add the coding of chroma components.

ACKNOWLEDGEMENT

This project has received funding from the European Unions Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie Grant Agreement No. 765911 (Re-alVision).

Fig. 1 Fig. 1 .

 11 Fig.1shows our proposed scheme. It consists of a base layer illustrated with red blocks and an enhancement layer denoted by cyan blocks. The layers are presented in detail in the following subsections.

NFig. 2 .

 2 Fig.2. Neural network architecture. The parameters in each block denote the number of filters, spatial extent of the filter, stride, the usage of the bias and the activation function.

Fig. 3 .

 3 Fig. 3. Test images for the quality evaluation.

Fig. 4 .

 4 Fig. 4. RD curve of the base layer, local RD curves obtained for each lambda value and the RD curve obtained by computing the convex hull around all the point (a) and the comparison of the base layer and the improved versions (b) for the LF image Bikes.

Table 1 .

 1 BD rate savings of the proposed hybrid approach with three types of enhancement layer against the base layer.

			Base
	Im.	ScalQuant	HEVCIntra	HEVCInter
	I01	-13.2317%	-27.4263%	-31.9636%
	I02	-4.9793%	-13.8663%	-16.0852%
	I04	-8.3554%	-22.1310%	-22.6794%
	I09	-14.4142%	-22.2486%	-31.2664%

Table 2 .

 2 BD rate savings of the best proposed Hybrid approach using HEVC Inter for the enhancement layer, against three anchors: learning-based base layer, MuLE, HEVC.