Laetitia Tufféry

Hendrik Davi
email: hendrik.davi@inrae.fr

Noelia López-García
email: noelia.lopez.garcia@slu.se

Eric Rigolot
email: eric.rigolot@inrae.fr

Florence Jean
email: florence.jean@inrae.fr

Anne Stenger
email: anne.stenger-letheux@inrae.fr

Francois Lefèvre

Davi H López-García

Laetitia Tuffery
email: laetitia.tuffery@inrae.fr

INRAE Cee-M

Davi Hendrik

López-García Noelia

François Lefèvre
email: françois.lefevre.2@inrae.fr

Adaptive measures for mountain Mediterranean forest ecosystem services under climate and land cover change in the Mont-Ventoux

Keywords: mountain Mediterranean forests, ecosystem services, climate change, adaptive scenarios, socio-ecological approach

published or not. The documents may come

L'archive ouverte

Introduction

Climate change (CC) and land use land cover change (LULCC) have become the two biggest threats to Mediterranean forests, profoundly altering the future of forest ecosystem services (ES) [START_REF] Lamarque | Plant trait-based models identify direct and indirect effects of climate change on bundles of grassland ecosystem services[END_REF]Choe and Thornes, 2017). Awareness of CC and uncertainty modifies forest management and political decisions from sustainable management to the adaptive strategy perspective [START_REF] Von Detten | Sustainability as a guideline for strategic planning? The problem of long-term forest management in the face of uncertainty[END_REF][START_REF] Hoogstra-Klein | Rational versus adaptive forest management planning: Exploratory research on the strategic planning practices of Dutch forest management organizations[END_REF][START_REF] Lidskog | Managing uncertainty: Forest professionals' claim and epistemic authority in the face of societal and climate change[END_REF][START_REF] Pukkala | Which type of forest management provides most ES[END_REF]. In the adaptive management perspective, three main issues must be considered: (1) the forest ecosystems will change as driven by climate and human decisions, (2) the demand on ES will change as driven by societal evolution and (3) the changes in forest ecosystems will have multiple impacts on multiple ES. The exploration of future delivery of ES in socioecological systems (SES) must integrate complex interrelations and feedback among economic activities (e.g., forestry), land use, local environmental conditions and ecological dynamics. Currently considered as an interesting approach by the scientific and decision maker communities, the future delivery of ES has been increasingly studied through exploratory scenarios [START_REF] Harmáčková | Future uncertainty in scenarios of ecosystem services provision: Linking differences among narratives and outcomes[END_REF][START_REF] Berg | Building scenarios for ecosystem services tools: Developing a methodology for efficient engagement with expert stakeholders[END_REF][START_REF] Van Vuuren | The Representative Concentration Pathways: An Overview[END_REF]. Scenario construction relies on qualitative narratives (for example, the evolution of local management policies or social preferences) and their transcription into quantitative changes (input variations such as temperature, quality of soil, demand for renewable energy, etc.) [START_REF] Martinez-Harms | Scenarios for land use and ecosystem services under global change[END_REF]. Mediterranean forests provide good case studies of multifunctional SES [START_REF] Croitoru | Valuing the non-timber forest products in the Mediterranean[END_REF][START_REF] Masiero | From failure to value: economic valuation for a selected set of products and services from Mediterranean forests[END_REF], as they are biologically and ecologically diverse, rapidly changing, and heterogeneous in ownership and management systems [START_REF] Fao | A new dynamic for Mediterranean forests[END_REF][START_REF] Martinez De Araño | A forest-based circular bioeconomy for southern Europe: visions, opportunities and challenges[END_REF]. In the case of Mediterranean forests, the expected impacts of global changes on future ecosystem functioning and sustainability have already been listed [START_REF] Fao | A new dynamic for Mediterranean forests[END_REF]. The literature discusses several factors of changes (CC, LULCC, and ecosystem management policies) using different methods (statistical/correlation analysis and process-based modeling) to evaluate their impacts on ES production and value. A recent literature review on models and scenarios used to assess the future of the Mediterranean forest ES [START_REF] Morán-Ordóñez | The use of scenarios and models to evaluate the future of nature values and ecosystem services in Mediterranean forests[END_REF] showed that only a few studies used processbased models while the majority of studies evaluated only one driver of change, i.e., CC. In particular, integrated prospective analyses of future ES delivery can be improved by a combined approach addressing three major challenges:

(i) Implementing biophysical modeling to better estimate the balanced effects of CC (e.g., longer growing season or increasing drought) and CO 2 fertilization on ecosystem functions; (ii) Combining biophysical models and managers' expert assessments to elaborate scenarios of LULCC in response to climate and socioeconomic changes; and (iii) Combining biophysical and economic approaches to evaluate the impacts of forest ecosystem changes on multiple-ES delivery at a territorial scale for different temporal horizons. Such integrated assessment of ES is still largely missing due to integration complexity, implying different ecosystems, ES, values, temporal and spatial scales, stakeholders and disciplines [START_REF] Tardieu | The need for integrated spatial assessments in ecosystem service mapping[END_REF]. The biophysical evaluation of ES is commonly based on empirical linear models [START_REF] Lamarque | Plant trait-based models identify direct and indirect effects of climate change on bundles of grassland ecosystem services[END_REF][START_REF] Martinez-Harms | Scenarios for land use and ecosystem services under global change[END_REF][START_REF] Schirpke | Future impacts of changing land-use and climate on ecosystem services of mountain grassland and their resilience[END_REF] that estimate the determinants of ES, including ecological processes and human practices. These relationships are first documented in different climate or land use contexts, accounting for current spatial and past temporal variations, then applied in future conditions. These empirical models assume a stable relationship between eco-physiological ecosystem characteristics and ES provision. As these effects are empirically estimated using past and present observations of the relationships between climate and ES, they do not account for complex and dynamic changes in forest functioning in the future climate context. Thus, they are considered unable to incorporate the effects of changing environmental conditions on tree and stand growth [START_REF] Fontes | Models for supporting forest management in a changing environment[END_REF]. By contrast, process-based models have thus been considered particularly convenient for the investigation of forest dynamics under new environmental conditions [START_REF] Fontes | Models for supporting forest management in a changing environment[END_REF][START_REF] Oddou-Muratorio | Integrating evolutionary, demographic and ecophysiological processes to predict the adaptive dynamics of forest tree populations under global change[END_REF]. Only few publications consider a combined approach for the study of global change impact on the future provision of ecosystem services at the landscape scale [START_REF] Kubiszewski | The future value of ecosystem services: Global scenarios and national implications[END_REF][START_REF] Small | The challenge of valuing ecosystem services that have no material benefits[END_REF][START_REF] Costanza | Valuing ecological systems and services. F1000[END_REF] or at a local scale [START_REF] Ding | Valuing climate change impact on European forest ecosystem[END_REF][START_REF] Bottalico | Modeling the influence of alternative forest management scenarios on wood production and carbon storage: A case study in the Mediterranean region[END_REF][START_REF] Masiero | From failure to value: economic valuation for a selected set of products and services from Mediterranean forests[END_REF]. However, integrated study of ES considering both supply side (ES provided by ecosystems according to ecological functions) and demand side (beneficiaries of ES) in spatially explicit terms is currently recognized as a powerful tool to help environmental and management planning decisions [START_REF] Tallis | Mapping and Valuing Ecosystem Services as an Approach for Conservation and Natural-Resource Management[END_REF][START_REF] Groot | Estimates of monetary values of ecosystem services[END_REF]Garcia-Nieto et al., 2013). The economic valuation of CC and LULCC impacts on the production of ES is an important step to solving potential conflicts between beneficiaries or establishing guidelines for sustainable forest management considering the potential trade-offs or synergies between ES [START_REF] Kubiszewski | The future value of ecosystem services: Global scenarios and national implications[END_REF]. In this work, we studied Mont-Ventoux Natural Regional Park. We developed a methodology to explore the effects of different scenarios of forest policies and CC on changes in carbon sequestration (regulation and maintenance ES) and timber production (provisioning ES) in spatially sensitive and dynamic ways, combining biophysical and economic values. We designed different forest management options under climate scenarios to estimate the resilience of Mediterranean forest ES according to their ecological and economic trajectories in two time horizons, 2050 and 2100. We first used a stakeholderbased analysis to estimate the evolution of ES production according to the LULCC. This first step was based on the InVEST model and its stakeholder matrices tool (i.e., InVEST models; [START_REF] Kareiva | Natural capital: theory and practice of mapping ecosystem services[END_REF]. Then, we used a biophysical model to assess the future state and functions of different forest types (i.e., climate impact on ecosystem functions) using simulations on nine forest types at five altitude levels for three soil water depths, two canopy openness values and two aspects (i.e., CASTANEA model; [START_REF] Dufrêne | Modelling carbon and water cycles in a beech forest Part I: Model description and uncertainty analysis on modelled NEE[END_REF][START_REF] Davi | Assessing drought-driven tree mortality with a process-based model[END_REF]. This combined model-based approach (i.e., LULC modeling and physiological assessment of forest functioning) allowed us to jointly explore the future effect of LULC and CC on forest ecosystems. Moreover, this approach combined a spatially explicit ES analysis and an economic valuation. We used an integrated approach to investigate how climate and LULCC factors impact the different ES and to analyze whether the changes in ES are due to climate scenarios or LULCC scenarios (including adaptation) and whether some LULCC scenarios can offset the effect of climate.

Methodology

Figure 1 presents the four steps of our methodology.

Land cover map of the studied area

The case study is on Mont-Ventoux Natural Regional Park in France, a Mediterranean forest socioecological system. Located in the southern part of France (the Provence region), the current perimeter of the Mont-Ventoux park covers a total area of 987 km². Mont Ventoux peaks at 1,912 m a.s.l. The territory is occupied by 57.7% forest and semi-natural environments and hosts a great diversity of natural habitats linked to a unique bioclimatic and geomorphological context (Barbero and Quezel, 1976). Because of its altitudinal gradient and its intermediate geographical location between the Alpine and Mediterranean regions, the Ventoux area presents a high diversity of environments. We characterized our study territory according to a typology compatible with the ecological and land cover models, i.e., a mono-specific and spatially representative typology. We used two main geographic databases: the French forest database by the National Geographic Institute in 2012 (BDForêt © ) and the Corine Land Cover database in 2012. Finally, we obtained 37 types of land covers comprising (see Fig. 2): nine locally predominant tree species (Quercus pubescens, Quercus ilex, Fagus sylvativa, Pinus nigra, Pinus sylvestris, Pinus halepensis, Pinus unicinata, Abies alba, Cedrus atlantica) and five bioclimatic stages (plain, Mediterranean, supra-Mediterranean, lower mountain, and higher mountain), five groups of companion species, moors, fire areas, and three aggregated LULCs (artificial areas, agricultural areas and various vegetation). This typology allowed us to work on mono-specific process-based and fine-scale (i.e., non-aggregated) LULC models. The spatial resolution was 25 m X 25 m grid (i.e., 1 ha = 16 pixels).

Exploratory scenario narratives

Our exploratory scenarios are based on two main drivers, CC and forest policy, defined over two time horizons: short term (2050) and long term (2100).

Establishment of a local stakeholder panel

The scenario construction started with the identification of the effects of the drivers of change on local forest ecosystems, which we achieved in close relationship with local stakeholders and managers [START_REF] Mckenzie | Developing Scenarios to Assess Ecosystem Service Tradeoffs: Guidance and Case Studies for InVEST Users[END_REF]. To do so, the first step was to identify a panel of stakeholders and managers to define and improve our scenarios and collect their knowledge on local forest resources and management. As recommended by the NatCap team in their InVEST guide and in the existing literature on this topic [START_REF] Berg | Building scenarios for ecosystem services tools: Developing a methodology for efficient engagement with expert stakeholders[END_REF][START_REF] Tallis | InVEST 2.6.0 User's Guide[END_REF], our stakeholder panel participated in four iterative meetings: scenario narrative construction (e.g., local forest policy orientation, biodiversity protection, and timber market evolution), data and inputs for the assessment (e.g., forestry planning by species and stage, harvest mass, cost of production and prices), future LULC and the transition matrix construction, and reviewing and improving the ES outputs. The stakeholder panel is presented in Appendix A.

Current local climate and downscaled climate change scenarios

We used the Safran meteorological analysis system [START_REF] Durand | A computer-based system simulating snowpack structures as a tool for regional avalanche forecasting[END_REF], which provides past meteorological daily input data on a regular grid (8 × 8 km). Then, we used two points of the grid located in the Ventoux, one at the north and the other at the south, to calibrate the climatic clines used in the simulations. To account for altitudinal effects on climate, local measurements were taken at ten forest sites targeted in the study from local weather stations (2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015). The minimum, maximum, and average temperatures, rainfall and relative humidity were recorded with a Prosensor HOBO Pro (RH/Temp; Onset Computer Corporation, Bourne, MA 02532, USA). Regression coefficients were estimated between these local stations and the Safran climate data at the two grid points. These regression coefficients were then used to generate daily climate data from 1959 to 2015 based on the long-term SAFRAN outputs using the equations described by Oddou-Muratorio and [START_REF] Oddou-Muratorio | Simulating local adaptation to climate of forest trees with a Physio-Demo-Genetics model[END_REF]. For the future climate, we used the Hadley Centre Global Environmental Model version 2 (HadGEM2) and the Coupled CNRM -CM/AOGCMs -atmosphere-ocean general circulation models. From all existing climate models, CNRM -CM5 is the coldest model that can also predict an increase in precipitation, and HadGEM2 -ES is the warmest model that can predict a decrease in rainfall [START_REF] Mcsweeney | Selecting CMIP5 GCMs for downscaling over multiple regions[END_REF]. We considered two CC scenarios: RCP 4.5 (Representative Concentration Pathway 1 ) and RCP 8.5. RCP 4.5 is the more optimistic scenario, and RCP 8.5 is the more pessimistic scenario. RCP 8.5 is estimated using the Model for Energy Supply Strategy Alternatives and their General Environmental Impact [START_REF] Moss | The next generation of scenarios for climate change research and assessment[END_REF] predicting a rising up to 1,370 ppm CO 2 -eq in 2100. These climate predictions were downscaled at the two studied points of the Safran grid using the 1960-2015 period as a reference using the R-package 'meteoland' [START_REF] De Cáceres | Coupling a water balance model with forest inventory data to predict drought stress: the role of forest structural changes vs. climate changes[END_REF]. Both scenarios follow the same trend until 2050, but they strongly diverge between 2050 and 2100 (see Appendix B).

Global forest policy scenarios and local forest management options

The global forest policy scenarios are based on the plausible futures of the European forestry sector published by the UNECE FAO section (European Forest Sector Outlook Study II -EFSOS II; FAO, 2011). Among the four EFSOS scenarios, we selected the three most relevant forest policies in our study area: business as usual (BAU), promotion of wood energy (WE) and biomass and carbon sequestration maximization (BIO). The first corresponds to the evolution of land cover (i.e., LULCC) related to "as usual" forest management, without adaptation, whereas the latter two, BIO and WE, correspond to forest management policies for adaptation to global changes. We adapted the EFSOS scenarios to our case study by translating them into a short narrative description that was then validated by the stakeholders. A brief description of the scenarios used in this work is in Appendix C (Table C1). The three forest policy scenarios and their related forest management options affect forest management regimes (e.g., thinning, harvest period). The forestry planning is accurately described in Table 1. These scenarios resulted in 12 combinations (2 time horizons × 2 RCP × 3 forest policies), producing our 12 scenarios (see Appendix C, Fig. 2).

Models

As presented in Figure 1, our two global change drivers, climate change and forest policy, were then used as inputs for the LULC and process-based models. All data sources used in these methods are presented in Appendix D.

Land cover change model

The three forest policy scenarios and their related forest management options affect forest management regimes (e.g., thinning, harvest period). The forestry planning is accurately described in Table 1. 1 The modeled climate scenario data follow the Representative Concentration Pathways (RCP).

For the construction of different potential futures, we used InVEST, which allows incorporation of complex scenarios combining climate, policy factors and local stakeholder expertise [START_REF] Kareiva | Natural capital: theory and practice of mapping ecosystem services[END_REF][START_REF] Tallis | InVEST 2.6.0 User's Guide[END_REF]. Forest management options under different climate change contexts are translated first into forest planning options, then into transition matrices and finally into maps. Stakeholder and manager interviews provided local information on tree species choices according to the potential mortality or adaptation of each species in a global change context. The methodology to express the transition probabilities from the corpus of information provided by the panel of stakeholders follows the recommendations of the InVEST guide. The transition probabilities depend on the percentage of change from one land cover type to another and the two surface areas of the original and final land covers. Given our 37 land cover types and 12 scenarios, the required information represents a complex overview. Stakeholders provided the rate of area change in each land cover for each scenario and at the two time horizons, given in matrix 𝐴𝐴, as follows:

To 𝐹𝐹′ 𝑗𝑗 From 𝐹𝐹 𝑖𝑖 𝐹𝐹′ 1 … 𝐹𝐹′ 37
Current surface area

𝑆𝑆 𝑖𝑖 𝐹𝐹 1 𝑃𝑃 1,1 … 𝑃𝑃 1,37 𝑆𝑆 1 ⋮ ⋮ ⋱ ⋮ ⋮ 𝐹𝐹 37 𝑃𝑃 37,1 … 𝑃𝑃 37,37 𝑆𝑆 37
Future surface area

𝑆𝑆′ 𝑗𝑗 𝑆𝑆′ 1 … 𝑆𝑆′ 37 ∆ Surface area 𝑂𝑂 𝑗𝑗 𝑂𝑂 1 … 𝑂𝑂 37
with 𝐹𝐹 𝑖𝑖,𝑗𝑗 as the LULC categories, from type 𝑖𝑖 = 1, … ,37 to type 𝑗𝑗 = 1, … ,37.

𝑃𝑃 𝑖𝑖,𝑗𝑗 is the proportion, between 0 and 1, of the total surface of land cover type 𝑖𝑖 that becomes a land cover type 𝑗𝑗. The current surface related to each LULC type is 𝑆𝑆 𝑖𝑖 ; the future surface, 𝑆𝑆′ 𝑗𝑗 , based on probabilities in matrix A, is computed as presented by equation (1); and 𝑂𝑂 𝑗𝑗 is the difference between 𝑆𝑆 𝑖𝑖 and 𝑆𝑆′ 𝑗𝑗 :

𝑆𝑆′ 𝑗𝑗 = ∑ 𝑃𝑃 𝑖𝑖,𝑗𝑗 37 𝑗𝑗=1 𝑆𝑆 𝑖𝑖 (1)
We finally obtained 12 future land cover maps (one per scenario).

Process-based model

CASTANEA is an eco-physiological process-based model used to predict water and carbon fluxes in forest stands [START_REF] Dufrêne | Modelling carbon and water cycles in a beech forest Part I: Model description and uncertainty analysis on modelled NEE[END_REF]. The canopy is assumed to be horizontally homogeneous and is vertically subdivided into a variable number of layers, each with the same amount of leaf area. One averaged tree is assumed to be representative of the whole stand; therefore, each tree behaves as a dominant tree. In CASTANEA, the tree structure is assumed to be a combination of five different functional parts: stems, branches, leaves, coarse and fine roots [START_REF] Davi | Modelling carbon and water cycles in a beech forest: Part II.: Validation of the main processes from organ to stand scale[END_REF]. In addition, a carbohydrate storage section is included. The main simulated output variables are canopy photosynthesis, maintenance and growth respiration, tree carbon allocation, soil heterotrophic respiration, tree transpiration, and ecosystem evapotranspiration [START_REF] Davi | Modelling carbon and water cycles in a beech forest: Part II.: Validation of the main processes from organ to stand scale[END_REF].

Species-specific parameters were defined using a literature review, remote sensing and forest inventory data (see Appendix F). The model was first evaluated at four FLUXNET sites, where some of the studied species are present, using growth data on all of the studied species [START_REF] Davi | Modelling carbon and water cycles in a beech forest: Part II.: Validation of the main processes from organ to stand scale[END_REF][START_REF] Davi | Sensitivity of water and carbon fluxes to climate changes from 1960 to 2100 in European forest ecosystems[END_REF]Guillemot et al., 2017;[START_REF] Lopez-Garcia | Assessing ecosystem functioning on Mediterranean forests in the context of global change: the case study of Mont Ventoux[END_REF]. Then, the CASTANEA model was used to simulate the water and carbon fluxes across an elevation gradient for different types of stands (open vs. closed, three types of soil, and soil water content) under past and future climate conditions (two contrasted climate models, HADGEM2 and CNRM -CM/AOGCMs and two RCPs, 4.5 and 8.5) with or without silviculture (Appendix B).

Mortality is estimated in the model by accounting for both the mortality induced by hydraulic cavitation and the mortality related to carbon starvation [START_REF] Adams | A multispecies synthesis of physiological mechanisms in drought-induced tree mortality[END_REF]. Adult mortality is due to either carbon starvation or hydraulic failure [START_REF] Davi | Assessing drought-driven tree mortality with a process-based model[END_REF]. We assumed that conifers die when the percentage of conductance loss is greater than 50% [START_REF] Brodribb | Hydraulic Failure Defines the Recovery and Point of Death in Water-Stressed Conifers[END_REF] and that angiosperms die when the percentage of conductance loss is greater than 88% [START_REF] Urli | Xylem embolism threshold for catastrophic hydraulic failure in angiosperm trees[END_REF]. We then estimated the average concentration of non-structural carbohydrates ([NSC] threshold ) per species below which trees die from past mortality data recorded in the IGN data. For Holm oak, this method generated an excessively high value for [NSC] threshold ; we then used the maximum value obtained for the other species.

Ecosystem services assessment

We studied two categories of ES, provisioning for timber production and regulating for carbon sequestration, at every time horizon and for all scenarios. To respect the specificity of each ES and the forest multifunctionality outputs, we selected adapted methods for the economic valuation. Thus, we used the price in euros in the French market for the timber production service, and we used a socioeconomic value, i.e., the social cost of carbon (SCC), for the carbon sequestration. In these two cases, the demand is commonly equal to the flow of ES: the biophysical amount of ES is directly converted into monetary units through market prices and costs. These methods are part of the market valuation methods (versus nonmarket valuation, for example, revealed and stated preference techniques), allowing for a better comparison and joint analysis of our pair of ES.

Carbon sequestration

At every time horizon, for every scenario and for each forest land cover type (𝑖𝑖 = 1, … , 35), carbon sequestration, 𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆 𝑖𝑖 , was assessed yearly using the CASTANEA model (t C .ha -1 .year -1 ) using the following equation:

𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆 𝑖𝑖 = 𝐺𝐺𝑃𝑃𝑃𝑃 𝑖𝑖 -𝑅𝑅𝑅𝑅 𝑖𝑖 -𝑅𝑅ℎ 𝑖𝑖 (2)
where 𝐺𝐺𝑃𝑃𝑃𝑃 𝑖𝑖 is the gross primary production; 𝑅𝑅𝑅𝑅 𝑖𝑖 the sum of all autotrophic respiration of coarse and fine roots, branches, stem and leaves; and 𝑅𝑅ℎ 𝑖𝑖 is the heterotrophic respiration of soil.

The economic valuation is based on the social cost of carbon (SCC). It corresponds to the reduction of future damage, i.e., avoided damage, caused directly or indirectly by the emission of an additional ton of CO 2 , by reducing theses emissions. We computed the economic value of carbon sequestration as follows:

𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆_𝑁𝑁𝑃𝑃𝑃𝑃 𝑖𝑖 = 𝑆𝑆𝑆𝑆𝑆𝑆× 𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆 𝑖𝑖 (1+𝑟𝑟) 𝑡𝑡
(3)

𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆_𝑇𝑇𝑁𝑁𝑃𝑃𝑃𝑃 𝑖𝑖 = 𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆_𝑁𝑁𝑃𝑃𝑃𝑃 𝑖𝑖 × 𝑆𝑆 𝑖𝑖 (4)
where 𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆_𝑁𝑁𝑃𝑃𝑃𝑃 𝑖𝑖 and 𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆_𝑇𝑇𝑁𝑁𝑃𝑃𝑃𝑃 𝑖𝑖 are the net present value (NPV) and the total net present value (TPNV) of carbon sequestration in euros, respectively, SCC is the current and fixed social cost of carbon, 𝐶𝐶𝐶𝐶 𝑖𝑖 is the carbon sequestration (equation ( 2)), t is the total number of years during the studied period (to obtain the mean annual sequestration), and 𝑟𝑟 is the discount rate. The surface of each forest land cover type 𝑖𝑖 is given by 𝑆𝑆 𝑖𝑖 (or 𝑆𝑆′ 𝑖𝑖 for the future).

The literature has proposed a large range of SCC. Then, as recommended in the literature on this topic [START_REF] Bottalico | Modeling the influence of alternative forest management scenarios on wood production and carbon storage: A case study in the Mediterranean region[END_REF], we simplified the analysis and selected only two values (expressed in € 2017 ), 44€/tCO2 in 2017 and 57€ 2017 /tCO2 in 2050 [START_REF] Watkiss | The Social Cost of Carbon[END_REF], and we assumed that those values remain constant between 2050 and 2100.

Timber production

Timber production depends on two main pieces of information, the wood volume in m 3 .ha -1 and the forest planning (thinning, % of timber harvested, harvest frequency), which vary across time horizons, climate scenarios and policy scenarios.

In CASTANEA the wood volume per ha, 𝑃𝑃ℎ𝑅𝑅 𝑖𝑖 for forest land cover type 𝑖𝑖 (𝑖𝑖 = 1, … , 35), is estimated yearly using the following equation:

𝑃𝑃ℎ𝑅𝑅 𝑖𝑖 = 𝐵𝐵 𝑠𝑠𝑡𝑡𝑠𝑠𝑠𝑠 𝑖𝑖 𝜌𝜌 𝑤𝑤𝑤𝑤𝑤𝑤𝑤𝑤 𝑖𝑖 ×𝑆𝑆𝑆𝑆 𝑖𝑖 (5)
where 𝐵𝐵 𝑠𝑠𝑆𝑆𝑠𝑠𝑠𝑠 𝑖𝑖 is the stem biomass per ha, 𝜌𝜌 𝑤𝑤𝑆𝑆𝑆𝑆𝑤𝑤 𝑖𝑖 is the wood density, and 𝑆𝑆𝑆𝑆 is the carbon content of wood (the stem biomass in the model is in gC/m 2 of soil). We did not include new recruitment in our simulated stands. Then, from the forest planning, we obtained the mass of wood harvested (m 3 /ha), 𝑇𝑇ℎ𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑃𝑃ℎ𝑅𝑅 𝑖𝑖 , as follows:

𝑇𝑇ℎ𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑃𝑃ℎ𝑅𝑅 𝑖𝑖 = 𝑃𝑃ℎ𝑅𝑅 𝑖𝑖 × ℎ𝑅𝑅𝑟𝑟𝑎𝑎 𝑟𝑟𝑟𝑟𝑆𝑆𝑠𝑠 (6)
where ℎ𝑅𝑅𝑟𝑟𝑎𝑎 𝑟𝑟𝑟𝑟𝑆𝑆𝑠𝑠 is the percentage of wood harvested for each period studied (see Table 1). The model used for the timber production net present value, 𝑁𝑁𝑃𝑃𝑃𝑃 𝑖𝑖 , and the total net present value, 𝑇𝑇𝑁𝑁𝑃𝑃𝑃𝑃 𝑖𝑖 , in euros is the following:

𝑇𝑇ℎ𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑃𝑃ℎ𝑅𝑅_𝑁𝑁𝑃𝑃𝑃𝑃 𝑖𝑖 = 𝑇𝑇ℎ𝑖𝑖𝑖𝑖𝑖𝑖𝑠𝑠𝑤𝑤𝑖𝑖ℎ𝑟𝑟 𝑖𝑖 ×(𝑃𝑃𝑟𝑟𝑖𝑖𝑆𝑆𝑠𝑠 𝑖𝑖 -𝐻𝐻𝑟𝑟𝑟𝑟𝑣𝑣 𝑐𝑐𝑤𝑤𝑠𝑠𝑡𝑡 𝑖𝑖 )-𝑀𝑀𝑟𝑟𝑖𝑖𝑖𝑖𝑆𝑆 𝑐𝑐𝑤𝑤𝑠𝑠𝑡𝑡 𝑖𝑖 (1+𝑟𝑟) 𝑡𝑡 (7) 𝑇𝑇ℎ𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑃𝑃ℎ𝑅𝑅_𝑇𝑇𝑁𝑁𝑃𝑃𝑃𝑃 𝑖𝑖 = 𝑇𝑇ℎ𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑃𝑃ℎ𝑅𝑅_𝑁𝑁𝑃𝑃𝑃𝑃 𝑖𝑖 × 𝑆𝑆 𝑖𝑖 (8)
This model includes the mass of wood harvested (m 3 /ha), 𝑇𝑇ℎ𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑃𝑃ℎ𝑅𝑅 𝑖𝑖 ; the cost of wood harvesting (€/m 3 /ha), 𝐻𝐻𝑅𝑅𝑟𝑟𝑎𝑎 𝑆𝑆𝑆𝑆𝑠𝑠𝑆𝑆 𝑖𝑖 ; and the cost of maintenance of each plot (€/ha), 𝑀𝑀𝑅𝑅𝑖𝑖𝑖𝑖𝑆𝑆 𝑆𝑆𝑆𝑆𝑠𝑠𝑆𝑆 𝑖𝑖 . We then computed November 2020

the TNPV, weighing the NPV by the total surface area (in ha). For the economic data, we used the marketplace value of the wood harvested from each plot, depending on the tree species (€/t) 𝑃𝑃𝑟𝑟𝑖𝑖𝑆𝑆𝑖𝑖 𝑖𝑖 , and the market discount rate 𝑟𝑟 (see Appendix E for all input details). Based on the literature (survey by [START_REF] Clark | The global wood market, Prices and plantation investment: an examination drawing on the Australian experience[END_REF][START_REF] Ding | Valuing climate change impact on European forest ecosystem[END_REF], we assumed that the timber selling price remains constant across the whole simulated period 2017-2100. We also justify this choice due to the very low impact of local shocks (supply or demand economic shocks resulting from the policy scenarios) on the timber marketplace.

Sensitivity analysis for the economic valuation

Economic forecasting analyses are characterized by a high level of uncertainty directly related to the selected level of discount rate applied to ES valuation. French legislation as well as many studies on land use planning projects recommend adopting a 4% discount rate (Lebègue Report, 2005;Quinet Report, 2013;Trivino et al.,2015). Other studies tested the effect of various discount rates: from 1% to 8% in [START_REF] Bottalico | Modeling the influence of alternative forest management scenarios on wood production and carbon storage: A case study in the Mediterranean region[END_REF] or from 1% to 5% in [START_REF] Pukkala | Which type of forest management provides most ES[END_REF]. However, to account for the scarcity of natural resources and the strong and positive correlation between ecological growth and economic growth [START_REF] Gosselin | Actualisation en forêt : pour quelles raisons et à quel taux ?[END_REF], the Stern Review and numerous studies recommend the use of a low discount rate of 1.5% [START_REF] Weitzman | Why the far-distant future should be discounted at its lowest possible rate[END_REF][START_REF] Stern | The Stern Review Report: The Economics of Climate Change[END_REF][START_REF] Weitzman | A review of the Stern Review on the economics of climate change[END_REF][START_REF] Gollier | Ecological discounting[END_REF][START_REF] Freeman | Biodiversity valuation and the discount rate problem[END_REF].

The selected discount rate can significantly modify the results. Therefore, we tested discount rates ranging from 1% to 4%. However, in this work we only compare the effect of the scenarios on the ES values based on the reference year, 2017. Thus, even if the results are modified in level, they are not modified in relative terms. We finally selected a discount rate of 1.5%, as recommended by the specific literature on economic valuation of biodiversity and ecosystems services [START_REF] Stern | The Stern Review Report: The Economics of Climate Change[END_REF][START_REF] Gollier | Ecological discounting[END_REF].

Results

To disentangle the effects of LULCC and forest adaptation policies on the socio-ecosystem evolution in various CC contexts, we present the results in three steps: (i) without LULCC and forest policy scenarios, i.e., the baseline; (ii) with LULCC but no adaptation policy, i.e., the BAU scenario; and (iii) with LULCC and adaptation policies, i.e., the BIO and WE scenarios (see Appendix C, Fig. C2). The analysis of the results was therefore conducted within time horizons, and between the different scenarios (CC, LULCC and adaptation) in 2050 and then in 2100. Carbon sequestration (t C .ha -1 ) and timber production (m 3 .ha - 1 ) in 2017 were computed from the past period, 1960-2017.

Baseline: climate change results

We defined the baseline changes as the CC effect alone under the climate scenarios used in the biophysical model. Figure 3 and Appendix I show the results in terms of biophysics and NPV for the climate change without LULCC and forest policy adaptation (BAU, BIO, and WE). Figure 3 presents the results of carbon sequestration and timber production for the two main species in terms of surface area in Mont Ventoux: Pinus halepensis and Quercus pubescens. We observed a weak effect of climate scenario on ES until 2050 but a slightly more noteworthy effect for the timber production starting in 2080 and the carbon sequestration of Quercus pubescens starting in 2050. Indeed, for this last species from 2050, the carbon sequestration related to RCP8.5 is approximately 5% higher than the carbon sequestration with RCP4.5. This is not the case for Pinus halepensis. Moreover, from 2080, the timber harvested is 5% to 10% higher with climate scenario 8.5. This is due to the fertilization effect of the CO 2 increase on photosynthesis and water use efficiency and to the higher increase in growing season under the 8.5 scenario. These effects counterbalance the higher mortality under the 8.5 scenario.

Assessing the mortality risk allows for more precisely including the effect of climate change on the future of species and therefore on the production of ES. Indeed, its variation is explained at 70-90% by climate scenarios and thus marginally explained by the LULCC and adaptation scenarios. More specifically, mortality on Mont Ventoux is on average 2% for RCP 4.5 and 4% for RCP 8.5 in 2050 and then 3% for RCP 4.5 and 5% for RCP 8.5 in 2100.

At the species level (Fig. 3 and Appendix I), the average difference between both RCPs is -/+ 15% for carbon sequestration and -/+ 5% for timber production. At the landscape level (Appendix J), the TNPV increases by 15% in 2050 from RCP4.5 to RCP8.5 and decreases by 15% in 2100. Finally, few effects of CC scenarios (RCP4.5 and RCP8.5) were found on ES production and their related TNPV both at the species and landscape levels.

Land cover changes

LULCC evolution differs according to CC and policy scenarios. Indeed, even if the signs of trends for all species remain stable across our scenarios, they strongly differ in terms of magnitude. For instance, the Pinus nigra surface increases from 22% to 134% according to the studied scenarios (Fig. 4).

Using the stakeholder transition matrices, we projected that the area covered by Quercus pubescens and the Pinus sylvestris will strongly decrease due to climate change. According to forest managers, these two species will not adapt to CC in their current location in the future; they cannot move to a higher altitude level, and the managers have already observed a strong drought-related mortality rate.

For Quercus ilex, Pinus halepensis, and Pinus nigra, forest managers anticipated different dynamics in the Ventoux. Due to the resilience of these species regarding CC, the area covered by these species is expected to increase because they will replace Quercus pubescens and Pinus sylvestris in different locations. Moreover, a policy scenarios effect occurs; Quercus ilex and Pinus nigra are the two species most favored by alternative forest policy scenarios: the BIO scenario supports Pinus nigra and Pinus halepensis development, and the WE scenario supports Quercus ilex and Pinus nigra.

Socio-ecosystem evolution with land cover changes and adaptation strategies

In this section, we present the ES variations, first including LULCC only (i.e., BAU scenario) and then including LULCC and adaptation strategies (i.e., BIO and WE scenarios). This allows us to observe the isolated and cumulative effects of the different drivers of change.

Figure 5 (and Appendix J) summarizes the changes in ES economic values according to these different factors. LULCC and adaptation positively affect the aggregated total value of ES in 2050 and 2100, increasing the TNPV by 30% to 120%, 50% on average, compared with the baseline. Moreover, in three cases, RCP4.5 and RCP8.5 in 2100 and RCP4.5 in 2050, the best scenario in terms of ecosystem services production aimed at promoting biomass; in the fourth case, RCP8.5 in 2050, the best scenario targeted the production of wood energy (see Fig. 5). The role of adaptation strategies on ES changes is highly variable according to the time horizon and the ES considered. The effects of LULCC are significant but opposite for the two ES (Appendix J). LULCC and the new species distribution at the landscape scale (see part 3.2) appear favorable for carbon sequestration but detrimental to timber production. The first effect is partly explained by the increase in Pinus halepensis and Pinus uncinata, and the second effect is partly explained by the decrease in Quercus pubescens. The effects of adaptation strategies, in addition to the LULCC, are significant and highly positive for the aggregated value of ES. However, each adaptation strategy affects ES differently. The biomass scenario is positive for carbon sequestration and negative for timber production, and the promotion of wood energy is positive for both ES but strongly favors timber production. As a result, the results vary according to the ES and the time horizons. The overall results are in line with the literature [START_REF] Bottalico | Modeling the influence of alternative forest management scenarios on wood production and carbon storage: A case study in the Mediterranean region[END_REF][START_REF] Chiabai | Economic Assessment of Forest Ecosystem Services Losses: Cost of Policy Inaction[END_REF][START_REF] Ding | Valuing climate change impact on European forest ecosystem[END_REF][START_REF] Ameztegui | Future trade-offs and synergies among ecosystem services in Mediterranean forests under global change scenarios[END_REF]. Finally, combined with LULCC and adaptation strategies, a climate scenario effect on the ES results exists. Figure 5 shows higher TNPV for RCP8.5 than for RCP4.5 for both horizons, 2050 and 2100. Indeed, comparing the TNPV in the case of LULCC and adaptation strategies with the baseline, we can observe that positive effects in RCP 8.5 are almost twice those in RCP 4.5.

All of these results assume a contrasting effect of the scenarios on the two ecosystem services, which is confirmed by Pearson's correlation analysis (see Appendix K). In the case of the spatial correlation analysis of the BAU, the correlation between the two ES is positive (around 0.7) in 2050 and then negative in 2100 (around -0.1). For adaptation scenarios (BIO and WE), the spatial correlations between ES vary more strongly for the WE in 2100, with correlations decreasing from 0.9 to 0.5, but the correlations remain stable for the BIO scenario.

Discussion

In this paper, we developed a socioeconomic and ecological methodology co-constructed with stakeholders to explore the effect of different scenarios of forest policies and CC on future ES provision (carbon sequestration and timber production) in a spatially sensitive and dynamic way, combining biophysical and economic values. This methodology, organized around the co-construction of change scenarios and transition matrices according to the stakeholders' local expertise, allows for the replicability of the study (as structured by the Wayfinder guideline2 ).

Our contribution is threefold. First, our predictions of the effect of CC scenarios on ES evolution are rather optimistic. Indeed, process-based predictions rely on physiological mechanisms to predict future ecosystem functioning, which corresponds to the fundamental niche of the species [START_REF] Keenan | Predicting the future of forests in the Mediterranean under climate change, with niche-and process-based models: CO2 matters![END_REF]. By contrast, predictions based on correlative species distribution models are usually more pessimistic because they derive their predictions from the realized bioclimatic niche of the species, which can be much more restricted than the fundamental niche for some species. Second, LULCC has the effect of increasing the changes and accentuating the difference between climate change scenarios. This result combines two effects: (1) a more pessimistic prediction by managers (probably the result of previous works based on correlative models and more pessimistic conclusions -particularly mapswhich have left their mark on foresters' management) and, as a consequence of this first point, (2) an anticipation action that tends to overanticipate and reinforce or even accelerate the expected changes. In this sense, the impact of adaptive actions, influenced by some models' predictions, can become greater than the impact of climate change itself. The action in adaptive strategy must be flexible and regularly reviewed according to the improvement of knowledge and scenario predictions. Third, even if the impact of economic choices is not major in our case study (we focus on a long-term management region rather than on short-term intensive management), our work shows that large-scale political orientations (i.e., adaptation strategies) impact ES. For the ES relationship, without adaptation scenarios, we can make wrong conclusions about synergy or conflict. Moreover, unlike the BIO scenario, which strengthens the synergies between ES, the WE scenario appears to weaken them considerably.

In line with several previous case studies, LULCC and adaptation strategies appear to be decisive drivers of the future ES evolution [START_REF] Ameztegui | Future trade-offs and synergies among ecosystem services in Mediterranean forests under global change scenarios[END_REF]Bréteau et al., 2019).

A few limitations and refinements of the study must be discussed to provide better modeling and policy orientations. First, rigorously comparing model and stakeholder results is difficult because simulation and forest management account for distinct elements. For example, managers not only anticipate the effect of CC on carbon sequestration but also the timber production and the future of the sectors. However, we can compare their global trends (upward/downward), and we can analyze what they are related to. Second, the NPV of carbon sequestration and the NPV of timber production (see Appendix J) differ by a hundred-fold. However, the "factor effect" is greatly attenuated in our analysis as we look at the results in terms of variations for each ES and then for the TNPV. Finally, regarding the use of process-based models, work still needs to be conducted in several directions. We may underestimate the tipping points in terms of hydraulic failure or the effects of multiple disturbances (drought, insects, and fire). In [START_REF] Davi | Sensitivity of water and carbon fluxes to climate changes from 1960 to 2100 in European forest ecosystems[END_REF], we showed that the main positive effect of climate change on the photosynthesis of evergreen species was due to the increase in atmospheric CO 2 , whereas the increase in vegetation duration mainly has an effect on deciduous species. We therefore looked at the relationship between GPP and CO 2 simulated for pines by CASTANEA in this study and compared this relationship with the measurements obtained in the FACE experiments on a Pinus ponderosa forest at a Duke site (Appendix G). CASTANEA reproduces the increasing water use efficiency (WUE) well for pine trees but underestimates the GPP fertilization (+23% instead of 30% at the Duke forest). However, we may still overestimate the fertilizing effect of CO 2 by the absence of acclimatization of photosynthesis to CO 2 and feedback of the nitrogen cycle. Moreover, huge uncertainties exist regarding soil organic carbon simulations (Smith et al., 2020), but CASTANEA reproduces the ranges of soil carbon stock and soil respiration, the relative variations between species and the effect of altitude well (Appendix H).

Conclusion

We conclude with directions that can be taken in future research. First, this work is valid for the public forest sector and management, but it does not include specific private forest silviculture. One obvious difference between private and public forests, which tends to disappear over time, is relying on the ES integration in forest management: if ES are normally integrated in public forest management, some incentives need to be implemented for private owners as these goods are public goods for most of them. Another possible extension addresses how some local impacts on demand or supply sides can diffuse at a more global scale.

Fig. 1 .Fig. 2 .Fig. 3 .

 123 Fig.1. Snapshot of the methodology in four steps.

Fig. 4 .

 4 Fig. 4. Land cover variations (in %) between 2017 and the two time horizons, 2050 (left side) and 2100 (right side) and scenarios (business as usual (BAU), biomass (BIO), wood energy (WE)), according to the expertise of the stakeholder panel for the five main tree species: Pinus halepensis, Pinus nigra, Pinus sylvestris, Quercus ilex, and Quercus pubescens.

Fig. 5 .

 5 Fig. 5. Evolutions (rate of change in %) of the total net present value (sum of the carbon and timber economic values in euros) between the scenarios (business as usual (BAU), biomass (BIO), wood-energy (WE)) and the baseline for two Representative Concentration Pathways, 4.5 and 8.5, and two time horizons, 2050 and 2100.

« Wayfinder represents a new generation of resilience practice that will guide development practitioners, policymakers and other change-makers navigating towards better futures. » https://graid.earth/projects/wayfinder/

Acknowledgements

This work was supported by a grant from the ERA-NET FORESTERRA project INFORMED. The authors wish to thank the following for their cooperation and constructive discussions regarding the issues raised in this paper: Olivier Delaprison, Jean Ladier, Anthony Roux, Philippe Dreyfus, the different members of the INFORMED workshop, and reviewers/participants at the Ecosystem Services Partnership (ESP) Regional Europe conference (San Sebastian, 2018) for their useful comments.


M-D S De Marco

T Hua

C Hamel

Samson

Position, Velocity, Attitude and Accelerometer-Bias Estimation from IMU and Bearing Measurements

This paper considers the problem of estimating the position, attitude and velocity of a rigid-body in a 3D space by fusing bearing measurements provided by a monocular camera with gyroscopic and accelerometer measurements provided by an Inertial Measurement Unit (IMU). The proposed deterministic observer is accompanied with an observability analysis, that points out the minimum number of image points (bearings) along with their configuration in the inertial frame, under which local exponential stability is guaranteed. The performance of the observer is demonstrated by performing experiments on a test-bed inertial-Visual sensor.

I. INTRODUCTION

The problem of full pose (position and attitude) estimation of a rigid body plays a major role for autonomous navigation of robotic vehicles. The orientation of a rigid body can be estimated from the output of low cost Inertial Measurement Unit (IMU) [START_REF] Mahony | Nonlinear complementary filters on the special orthogonal group[END_REF], whereas translational velocity and position can not be estimated from accelerometer measurement via sole integration due to unbounded growth in the estimation error caused by the sensor noise and bias [START_REF] Woodman | An introduction to inertial navigation[END_REF].

In order to circumvent this problem, IMUs are usually complemented by other sensors such as Global Positioning System (GPS) which allow to bound the error [START_REF] Farrell | Aided Navigation: GPS with High Rate Sensors[END_REF]. However, in many applications of interest such as indoors and urban canyon environments the GPS signal may be unreliable or heavily degraded. For this reason, the usage of other type of sensors such as laser scanners (LIDAR), Ultra-Wide Band sensors (UWB), monocular cameras and stereo cameras have been considered. Early approaches for state estimation were based on extended Kalman filters (EKF) [START_REF] Anderson | Optimal Filtering[END_REF], unscented Kalman filters (UKF) [START_REF] Julier | New extension of the Kalman filter to nonlinear systems[END_REF] and particle filters [START_REF] Samson | An Introduction to Sequential Monte Carlo Methods[END_REF]. The last fifteen years have witnessed an increased interest in deterministic non-linear observer design for mechanical systems with symmetry that can provide stronger stability results with large domain of attraction [START_REF] Aghannan | An intrinsic observer for a class of lagrangian systems[END_REF], [START_REF] Bonnabel | Non-linear symmetrypreserving observers on lie groups[END_REF], [START_REF] Mahony | Observers for kinematic systems with symmetry[END_REF]. Based on the type of sensors considered for the particular application at hand, several research directions have been investigated.

For instance, in the particular case of landmarks measurements (which can be retrieved by using depth cameras, range measurements and bearings, stereo cameras) an invariant extended Kalman filter (IEKF) on the group of double homogeneous matrices SE 2 (3) has been proposed in [START_REF] Barrau | The invariant extended kalman filter as a stable observer[END_REF]. In the same context, an observer on the Lie group of rigid body S. de Marco, M-D. Hua are with I3S, Université de Nice, Sophia-Antipolis, France sdemarco(hua)@i3s.unice.fr.

T. Hamel motions SE(3) based on the Lagrange-d'Alembert principle has been derived in [START_REF] Izadi | Rigid body pose estimation based on the lagrange-d'alembert principle[END_REF]. In case of velocity and landmark measurements a non-linear observer on SE(3) has been presented in [START_REF] Vasconcelos | A landmark based nonlinear observer for attitude and position estimation with bias compensation[END_REF]. Recently, a non-linear observer that estimates the position, velocity, attitude and gyro bias by fusing multiple measurements such as full position in the inertial frame (GPS), range (UWB), bearings in the inertial frame (MOCAP) and altimeter along with IMU measurements has been considered in [START_REF] Berkane | Position, velocity, attitude and gyro-bias estimation from imu and position information[END_REF].

In case of bearing measurements in the body frame, the estimation problem is referred to as Perspective-n-Point problem (PnP). Classical methods for the solution of the PnP problem are based on algebraic approaches [START_REF] Grunert | Das pothenotische problem in erweiterter gestalt nebst über seine anwendungen in der geodäsie[END_REF], [START_REF] Dhome | The inverse perspective problem from a single view for polyhedra location[END_REF] and iterative algorithms based on gradient search [START_REF] Haralick | Pose estimation from corresponding point data[END_REF], [START_REF] Lu | Fast and globally convergent pose estimation from video images[END_REF]. The PnP problem has been subsequently studied in the context of observer design. In particular, in [START_REF] Huster | Relative position estimation for auvs by fusing bearing and inertial rate sensor measurements[END_REF] an UKF like filter is derived considering only one source point and IMU measurements. The convergence of this observer depends on the motion of the system, since an implicit triangulation is required. In [START_REF] Aguiar | Minimum-energy state estimation for systems with perspective outputs and state constraints[END_REF], the authors present a minimum-energy state estimator with constraints considering bearing-only measurements. A non-linear observer evolving on SE(3) from bearing, angular and linear velocity measurements is derived in [START_REF] Baldwin | A nonlinear observer for 6 dof pose estimation from inertial and bearing measurements[END_REF]. In a prior work [START_REF] Hamel | Riccati observers for the nonstationary pnp problem[END_REF], a Riccati observer for the PnP problem based on angular and linear velocities measurements has been proposed.

In this paper, we propose a deterministic Riccati observer, based on the general framework presented in [START_REF] Hamel | Riccati observers for the nonstationary pnp problem[END_REF], providing pose, linear velocity and accelerometer-bias estimation by fusing IMU and bearings measurements, without relying on linear velocity measurements. The present work is motivated by a large number of applications in which linear velocity measurements are not available, and a reliable estimate of the linear velocity and of the pose is needed for control purposes. This is the case of low cost small-scale vehicles that operates in GPS denied environment for which the cost of body velocity sensors such as sensors based on Doppler effect or ad-hoc inertial velocity sensor systems are usually prohibitive.

II. PRELIMINARY MATERIAL A. Notation

We denote by R n the n-dimensional Euclidean space and by |x| the Euclidean norm of the vector x ∈ R n . Let B n r := {x ∈ R n : |x| ≤ r} denotes the closed ball in R n with radius r. The null matrix and the identity matrix of dimension n × n are denoted respectively by 0 n and I n .

The null matrix of dimension m × n is denoted by 0 m×n .

The Special Orthogonal group of order three is denoted by

SO(3) := {R ∈ R 3×3 : det(R) = 1, R R = RR = I 3 }.
For any x, y ∈ R 3 , we denote by x × the skew symmetric matrix associated to the cross product i.e. x × y = x × y, where × is the vector cross product in R 3 . Let S n := {x ∈ R n+1 : |x| = 1} denotes the n-dimensional sphere embedded in R n+1 with radius equal to one. For any p ∈ S 2 , we denote by Π p := (I -pp ) the projection onto the tangent space of the unit sphere S 2 at point p.

B. Riccati Observer for a Special class of Systems

The observer proposed in this paper is based on the Riccati Observer framework developed in [START_REF] Hamel | Riccati observers for the nonstationary pnp problem[END_REF]. In particular we consider the following class of nonlinear systems whose state

x := [x 1 , x 2 ], with x1 ∈ B n1
r and x 2 ∈ R n2 , and outputs y ∈ R m evolve according to

     ẋ = A(x 1 , t)x + u 1 u 2 + O(|x 1 | 2 ) + O(|x 1 ||u 1 |) y = C 1 (x 1 , x 2 , t)x 1 + C 2 (x 1 , x 2 , t)x 2 (1)
where A(x

A(x 1 , t) = A 11 (t) 0 n1×n2 A 21 (x 1 , t) A 22 (t) , (2)
and

C := [C 1 , C 2 ] ∈ R m×(n1+n2
) is a continuous matrixvalued function uniformly bounded w.r.t. t and uniformly continuous w.r.t. x.

In the sequel x1 will represent the attitude error between the camera orientation and the estimated one. We recall here that for a Riccati observer design one has to use a minimal parametrization of SO(3), however there are many ways to minimally parameterize SO(3) (e.g. Euler angles, vector part of Rodriques unit quaternion, Cardan angles). Consequently, system (1) can be seen as a handy way to take into account all the possible parameterizations, with the related nonlinear errors and first-order approximations, without choosing a priori one of them.

Let x2 be an estimate of x 2 and consider the following observer

     u 1 = -K 1 (y -C 2 x2 ) ẋ2 = A 22 x2 + u 2 + K 2 (y -C 2 x2 ) Ṗ = AP + P A -P C Q(t)CP + V (t) (3)
with P (0) a symmetric positive definite matrix, Q and V bounded continuous symmetric positive semidefinite matrixvalued functions, and

K := [K 1 , K 2 ] given by K = k(t)P C Q, with 0.5 ≤ k(t) ≤ k max < ∞.
Let us recall the definition of uniform observability of the following Linear Time Varying (LTV) system,

ẋ = Ā(t)x + B(t)u y = C(t)x (4)
which is instrumental for the stability of the observer in (3).

Given a continuous n × n-dimensional matrix-valued function Ā(t), a continuous m × n-dimensional matrixvalued function C(t), and a continuous m × m-dimensional matrix-valued function Q(t) > 0 (a positive definite matrix for all t ∈ R + ), the Riccati observability Gramian associated with the triplet ( Ā, C, Q) is the non-negative definite matrixvalued function defined by

W Ā, C Q (t, t + δ) := 1 δ t+δ t Φ (s, t) C (s)Q(s) C(s)Φ(s, t)ds (5)
where Φ(t, t 0 ) is the transition matrix associated with Ā(t), i.e. d dt Φ(t, t 0 ) = Ā(t)Φ(t, t 0 ), Φ(t 0 , t 0 ) = I. If Ā(t) and C(t) are bounded and if there exists δ > 0 and > 0 such that W Ā, C In (t, t + δ) > I n for all t ≥ 0, then we say that the pair ( Ā, C) is uniformly observable.

Define x2 := x 2 -x2 , then, from Corollary 3.2 in [START_REF] Hamel | Riccati observers for the nonstationary pnp problem[END_REF], (x 1 , x2 ) = (0, 0) is locally exponentially stable if Q(t) and V (t) are both larger than some positive matrix and the pair (A(0, t), C(0, x 2 , t)) is uniformly observable.

C. Measurements and System Equations

Let {I} denotes an inertial (fixed) frame of reference and {B} denotes a frame of reference attached to the camera. Let R ∈ SO(3) denote the orientation of the camera frame {B} w.r.t. the inertial frame {I}. Let ξ and v denote the position and the velocity of the camera frame {B} w.r.t. {I} expressed in {B}, respectively. Then, the camera second order kinematics equation are

     Ṙ = Rω × ξ = -ω × ξ + v v = -ω × v + a B + gR e g (6)
where a B denotes the specific acceleration expressed in the camera frame {B}, g denotes the gravitational acceleration and the unit vector e g denotes the gravitational acceleration direction expressed in the inertial frame.

The camera is equipped with a strap-down 6DOF IMU composed of a 3-axis gyroscope, which measures the angular velocity ω m ∈ R 3 , and a 3-axis accelerometer, which measures the specific acceleration a m B ∈ R 3 ; both measurements are expressed in the camera frame {B}. The rate gyros and accelerometer are usually modeled as

ω m = ω -b w -µ w a m B = a B -b a -µ a
where µ w , µ a denote stochastic additive noises and b w , b a denote constant or slowly time-varying biases.

In this paper we ignore the noise (µ w = µ a = 0) and for the sake of simplicity we assume that the gyro bias is zero (b w = 0) or already well calibrated by averaging the gyros measures when the IMU is at rest. The accelerometer calibration is more difficult because it cannot be performed off-line in a resting position and an incorrect calibration can yield drifts in the error (double integration). The following simplified sensor model is considered in this paper

ω m = ω a m B = a B -b a , ḃa = -ω × b a . ( 7
)
We assume to have a collection of n source points whose coordinates Pi w.r.t. the inertial frame {I} expressed in the basis of {I} are known, and whose calibrated projective coordinates p p i onto the camera plane expressed w.r.t. the camera frame in the basis of {B} are also known, see Fig. 1. Denoting by P i the coordinates of the ith source point w.r.t. the body-frame {B} expressed in {B}, and by z i its third component, one verifies p p i = P i /z i . Rather than using the perspective outputs typically used in computer vision algorithms, we use direction output representation of the measurement p i := p p i /|p p i | corresponding to projection onto a virtual spherical image plane and differing from perspective outputs only in the scaling. The unitary vector p i characterizes the bearing (direction) of P i in the basis of {B}. Let ξ represent the position of the camera frame with respect to the inertial frame expressed in {I}, then one verifies

p i := p p i |p p i | = -R ξ -Pi |ξ -Pi | .

III. RICCATI OBSERVER DESIGN

A. Problem Formulation and observer design

The problem at hand is to provide an estimate of ξ, v, b a and R given the measurement of the specific acceleration a m B , the gyro measurement ω m , and given a set of source points whose positions in the inertial frame Pi and bearings p i are known.

To this purpose, we proceed by showing that the system considered can be written in the form [START_REF] Samson | An Introduction to Sequential Monte Carlo Methods[END_REF]. In order to do so, as already mentioned in the previous section, we have to define the attitude error between the camera orientation and the estimated one, and then consider a minimal parametrization of the three-dimensional Lie Group SO(3).

Let R be an estimate of R and consider the following pre-observer attitude kinematics

Ṙ = Rσ × , (8)
where σ is an innovation term to be designed.

Denoting by R := R R the attitude error, one verifies that the error kinematics satisfy

Ṙ = R(ω -R σ) × ,
and by considering the gyros measurement model in [START_REF] Berkane | Position, velocity, attitude and gyro-bias estimation from imu and position information[END_REF] one has

Ṙ = R(ω m -R σ) × . (9)
The attitude parametrization here chosen is the vector part λ ∈ B 3 1 of the Rodriguez unit quaternion Λ = ( λ0 , λ), where λ0 ∈ B 1 1 denotes the scalar part. The rotation matrix error R is related to Λ by the Rodriguez formula

R = I 3 + 2 λ× ( λ0 I 3 + λ× ).
Then, in view of ( 9) and using the fact that λ0 = 1 -| λ| 2 , one has

R = I 3 + 2 λ× + O(| λ| 2 ) 2 λ = ω -ω m × 2 λ + O(| λ||ω|) + O(| λ| 2 ) ( 10
)
where ω := ω m -σ.

Concerning the dynamics of v one has

v = -ω × + gR R R e g + a B = -ω × + g R R e g + a B .
Substituting R with the first order approximation in ( 10)

yields v = -ω × v + g R e g + a B + g( R e g ) × 2 λ + O(| λ| 2 ),
and by taking into account the sensor model one obtains

v = -ω m × v + g R e g + a m B + b a + g( R e g ) × 2 λ + O(| λ| 2 ). (11
) Note that the term g( R e g ) × 2 λ that appears in the equation above represents a coupling term between the angular and the linear dynamics. This term will appear in the lower diagonal part of the A(x 1 , t) matrix in (1), namely A 21 (x 1 , t), making the problem of computing the transition matrix associated with A(x 1 , t) very difficult. In order to simplify the subsequent analysis we consider the following change of variable ba := b a + g( R e g ) × 2 λ.

(

) 12
Let ba denotes an estimate of ba , and denote by ba := baba the new bias estimation error. One verifies that in steady state ( λ, ba ) → (0, 0) implies ba = b a , and thus the new change of variable is a valid choice for the accelerometer bias estimation. We proceed, then, by computing the first order approximation of the time derivative of the new bias

ḃa = d dt b a + g(R e g ) × 2 λ + O(| λ| 2 ) = -ω × b a -ω × g( R e g ) × 2 λ + g( R e g )ω+ +O(| λ||ω|) + O(| λ| 2 ) + d dt O(| λ| 2 ), one verifies from (10) that d dt O(| λ| 2 ) = O(| λ||ω|)+O(| λ| 2 )
, bearing in mind the expression of ba in [START_REF] Hamel | Riccati observers for the nonstationary pnp problem[END_REF] and by taking into account the sensor model one has

ḃa = -ω m × ba + g( R e g ) × ω + O(| λ||ω|) + O(| λ| 2 ). ( 13
)
Concerning the output function for the ith source point one has

0 = -|ξ -Pi |Π pi p i = Π pi ξ -Π pi R Pi = Π pi ξ -Π pi R R Pi = Π pi ξ -Π pi (I + 2 λ× + O(| λ| 2 )) R Pi = Π pi ξ -Π pi R Pi -Π pi (R Pi ) × 2 λ + O(| λ| 2 )
and thus by setting y i := Π pi R Pi , one obtains 6), ( 10), [START_REF] Grunert | Das pothenotische problem in erweiterter gestalt nebst über seine anwendungen in der geodäsie[END_REF], ( 13) one obtains the system equation ( 1) with

y i = Π pi ξ -Π pi (R Pi ) × 2 λ + O(| λ| 2 ) ( 14
)
for i = 1, . . . , n. Finally, setting x1 = 2 λ, x 2 = [ ξ , v , b a ] , y = [y 1 , . . . , y n ] from (
                           A 11 = -ω m × , A 21 = 0, u 1 = ω A 22 =    -ω m × I 3 0 3 0 3 -ω m × I 3 0 3 0 3 -ω m ×    , u 2 =    0 3×1 a m B + g R e g g( R e g ) × w    C 1 =     -Π p1 ( R P1 ) × . . . -Π pn ( R Pn ) ×     , C 2 =     Π p1 0 3 0 3 . . . . . . . . . Π pn 0 3 0 3     (15)
Then, the observer is given by ( 3), ( 8) with σ = ω m -u 1 .

B. Observability Analysis

In this section we provide sufficient conditions that ensure good conditioning of the Continuous Riccati Equation (CRE) and the exponential stability of the origin of the observer error as well.

Proposition 1. Consider system (6) with measurement given by [START_REF] Berkane | Position, velocity, attitude and gyro-bias estimation from imu and position information[END_REF], and the observer given by (3) along with (8) and [START_REF] Izadi | Rigid body pose estimation based on the lagrange-d'alembert principle[END_REF]. Choose Q(t) and V (t) of (3) larger than some positive matrix. Assume that camera's position ξ, linear velocity v, angular velocity ω are bounded for all time. If the pair (A(0, t), C(0, x 2 , t)) is uniformly observable then, the origin of the observer error (x 1 , x2 ) = (0, 0) is locally exponentially stable.

Proof. The result is a direct application of Corollary 3.2 in [START_REF] Hamel | Riccati observers for the nonstationary pnp problem[END_REF].

As shown in [START_REF] Hamel | Riccati observers for the nonstationary pnp problem[END_REF], the uniform observability of pair (A(0, t), C(0, x 2 , t)) is required to ensure the local uniform exponential stability of the observer error. We show next the condition under which the uniform observability is violated.

Proposition 2. The pair (A(0, t), C(0, x 2 , t)) is not uniformly observable if the following conditions are fulfilled:

• The number of source points is less than or equal to two (n ≤ 2). • All source points are aligned independently of their number (n ≥ 3). • In case of three non-aligned source points (n = 3):

-Static Case: The camera lies on the circular cylinder generated by the circle passing through the three source points and whose axis is orthogonal to the plane containing the source points, known in literature as dangerous cylinder. -Moving Case: The camera moves along a straight line orthogonal to the plane containing the source points and passing through a source point. • In case of four and more non-aligned source points (n ≥ 4):

-Static Case: The source points are located on a horopter curve, which is uniquely determined by four points, and the camera frame {B} lies at the origin of this curve.

Proof. In view of ( 15) one verifies that the transition matrix associated with the state matrix A in ( 2) is of the form

Φ(s, t) =     R(s,t) 0 3 0 3 0 3 0 3 R(s,t) R(s,t) (s -t) R(s,t) (s-t) 2 2 0 3 0 3 R(s,t) R(s,t) (s -t) 0 3 0 3 0 3 R(s,t)     (16) where R(s,t) := R (s) R(t).
In view of (5) the pair (A, C) is not uniformly observable if there exists a non-zero constant vector ν ∈ R 12 such that

C(s)Φ(s, t)ν = 0, with s ∈ [t, t + δ] and t constant.
Substituting the expression of C(0, x 2 , s) in ( 15) and the one of Φ(s, t) in ( 16) on the equation above yields

Π p I i -Pi× I 3 I 3 (s -t) I 3 (s -t) 2 2 ν = 0 (17)
for ∀i ∈ {1, . . . , n} and s ∈ [t, t + δ]; where p I i := -(ξ(t) -P )/|ξ(t) -P | denotes the bearing of the ith source point expressed in the inertial frame {I}, and ν = R(t)ν.

Let ν := [ν 1 , ν 2 , ν 3 , ν 4 ] , with ν1 , ν2 , ν3 , ν4 ∈ R 3 . Equality (17) characterizes body position trajectories along which the body pose is not uniformly observable. For the case of only one bearing measurement, namely n = 1, it is trivial to verify that there exists an infinity of solutions for ν such that ( 17) is satisfied. For instance ν1 = P1 /| P1 |, ν2,3,4 = 0 is a non null solution to [START_REF] Lu | Fast and globally convergent pose estimation from video images[END_REF].

Consider now n ≥ 2. By choosing a particular value of s equal to t (i.e. s = t) one deduces from [START_REF] Lu | Fast and globally convergent pose estimation from video images[END_REF] that

Π p I i -Pi× I 3 ν1 ν2 = 0 (18)
for ∀i ∈ {1, . . . , n}, and subsequently (from ( 17) and ( 18))

(s -t)Π p I i ν3 + (s -t) 2 2 Π p I i ν4 = 0, ∀s ∈ [t, t + δ] This implies that Π p I i ν3 = 0 Π p I i ν4 = 0 (19)
∀i ∈ {1, . . . , n}, which in turn implies ν3 = ν4 = 0. Then from [START_REF] Mahony | Nonlinear complementary filters on the special orthogonal group[END_REF] the proof follows the same arguments as the proof in [START_REF] Hamel | Riccati observers for the nonstationary pnp problem[END_REF] for the case of mobile velocity measurements.

IV. EXPERIMENTAL RESULTS

In this section, we present experimental results to demonstrate the performance of the proposed observer. The experiments were performed on a custom platform composed of a low-cost IMU myAHRS+ and an oCam camera (see Figure 2). The camera mounts a standard M12 lens with 3.6mm focal length and 65 • field of view. Images are acquired with a resolution of 640 × 480 pixels at 30 frames per second. Inertial data from the IMU are acquired at rates of 100Hz. The point detection is performed with the ArUco library [START_REF] Francisco | Speeded up detection of squared fiducial markers[END_REF] which is included in OpenCV. Only one fiducial marker (size: 21 cm, dictionary: DICT 4X4 50) has been used during the experiments (Fig. 2). An Optitrack motion capture system (MOCAP) is used together with the markers mounted on the camera and on the fiducial marker in order to provide full pose ground truth measurements. Note that the imperfect positioning of the markers and the initialization process required in Optitrack introduce systematic errors in the ground truth pose.

A discrete version of observer (3) has been implemented on C++ and runs on an HP laptop with Intel R Core TM i7-6700HQ and Nvidia R Quadro M2000M graphic card.

The observer proposed in this paper is compared with the solvePnP() algorithm in OpenCV and with the ground truth data. The time evolutions of the roll, pitch and yaw Euler angles are depicted respectively in Figures 6-7-8 . The proposed observer shows better performances in the attitude estimate with respect to ArUco algebraic method, the difference between the two algorithm is evident when the camera performs complex dynamic motions (sec 145-185). Indeed, during secs 145-185 the camera is moving fast in the third component (Figure 5) and yawing (Figure 8) at the same time; the ArUco algorithm during this time window is switching between two possible solutions for the pitch (Figure 7) and roll (Figure 6) Euler angles.

V. CONCLUSIONS

In this paper, a novel deterministic observer based on the Riccati framework [START_REF] Hamel | Riccati observers for the nonstationary pnp problem[END_REF] for the pose and linear velocity estimation of a rigid-body along with accelerometer bias correction has been proposed. The conditions under which the state of the system is not uniformly observable, and thus local exponential stability not guaranteed, have been characterized in terms of the number and position of the source points. Experimental results have been provided as complement to the theoretical approach.

In this work, only one parametrization among the many possible local parametrizations of the pose error has been studied. The considered is for applications in which the measurements are given in the body-fixed frame.

Future work will focus on different types of measurements that could be available (e.g. barometer, altimeter) along with different pose error parametrizations, and their effects on the accuracy of the observer.

Fig. 1 .

 1 Fig. 1. Intuitive representation of inertial coordinates Pi , planar projective coordinates p p i and spherical projective coordinates p i of the ith source point.

Fig. 2 .

 2 Fig. 2. Experimental platform.

Fig. 3 .

 3 Fig. 3. First component of the camera position ξ from MOCAP ground truth, Riccati Observer and ArUco OpenCV algorithm.

Figures 3 -

 3  show the time behavior of the estimated position from the Riccati observer (red), ArUco OpenCV (yellow) and the ground truth data (blue). One notes that the position estimates obtained from the Riccati observer and

Fig. 4 .

 4 Fig. 4. Second component of the camera position ξ from MOCAP ground truth, Riccati Observer and ArUco OpenCV algorithm.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Third component of the camera position ξ from MOCAP ground truth, Riccati Observer and ArUco OpenCV algorithm.

Fig. 7 .

 7 Fig. 7. Pitch Euler angle from MOCAP ground truth, Riccati Observer and ArUco OpenCV algorithm.

Fig. 8 .

 8 Fig. 8. Yaw Euler angle from MOCAP ground truth, Riccati Observer and ArUco OpenCV algorithm.

ACKNOWLEDGMENT

This research was supported by the French DGA-Rapid project "Alcyon", the French EQUIPEX Robotex project and the French Astrid CONGRE project ANR-18-ASTR-0006.


Ninad Manerikar

Minh-Duc Hua

Simone De Marco

Tarek Hamel

Riccati observer design for homography decomposition

This paper presents a novel approach for tackling the classical problem of homography decomposition. The novelty of this paper lies in the design of a deterministic Riccati observer for addressing the homography decomposition problem instead of solving it on a frame-by-frame basis like traditional algebraic approaches. By exploiting the temporal correlation of image sequences and the natural low-pass response of the observer, we believe that the resulting estimates are robust and less prone to measurement noise. Moreover, an extensive observability analysis that points out a sufficient uniform observability condition under which local exponential stability is granted has been carried out. The large domain of convergence and good performance of the proposed observer have been demonstrated through simulation results.

I. INTRODUCTION

When a robotic vehicle equipped with a monocular camera navigates in an environment composed of (near) planar surfaces such as man-made structures or a ground terrain distantly viewed by an unmanned aerial vehicle, then the planar nature of the environment can be exploited in vision processing algorithms as well as vision-based control algorithms. Two different images of the same planar surface are typically related by the so-called homography mapping that essentially encodes the camera's pose (i.e. position and orientation), the distance between the camera and the scene and its normal vector in one single matrix [START_REF] Hartley | Multiple view geometry in computer vision[END_REF]. Homographies have been extensively used in robot control as a vision primitive, and one of the most successful visual servo control paradigms is the 2 1 2 D approach proposed by Malis et al. [START_REF] Malis | 2 1/2 D visual servoing[END_REF] that relies on the extraction of the camera displacement (i.e. orientation and translation up to a scale factor) and the scene's normal vector from homography (assuming that the camera's intrinsic parameters are known). This process of extraction is termed as Euclidean homography decomposition, also called Euclidean reconstruction from homography.

Classical approaches for homography decomposition such as Faugeras SVD-based [START_REF] Faugeras | and structure from motion in a piecewise planar environment[END_REF], Zhang SVD-based [START_REF] Zhang | 3D reconstruction based on homography mapping[END_REF] algorithms use the singular value decomposition to obtain numerical solutions. Malis and Vargas [START_REF] Malis | Deeper understanding of the homography decomposition for vision-based control[END_REF] lately solved the homography decomposition problem with an analytical approach, making it more suitable for real-time robot control applications. In that work, not only analytical expressions of the decomposition elements are provided, valuable insights regarding the robustness of vision-based control laws are also analyzed. In summary, to our knowledge all existing homography decomposition methods belong to algebraic category that only focuses on solving the homography decomposition N. Manerikar, M.-D. Hua, S. De Marco, T. Hamel are with I3S, Université Côte d'Azur, CNRS, Sophia Antipolis, France, manerika(hua, sdemarco, hamel)@i3s.unice.f r.

problem on a frame-by-frame basis, but not on filtering measurement noise. The precision of the decomposition elements (particularly of the estimated normal vector) is thus highly prone to noise, especially when the camera's translation is small [START_REF] Malis | Deeper understanding of the homography decomposition for vision-based control[END_REF]. The problem even becomes degenerate when the latter vanishes. Robotic vehicle applications, however, provide temporal sequences of images together with inertial measurements (from an embedded Inertial Measurement Unit). Therefore, it seems natural to exploit the temporal correlation rather than try to solve homography decomposition for each pair of image frames.

In this paper, a novel direction for solving the homography decomposition by exploiting the system dynamics is explored. The proposed solution is developed in the form of a nonlinear observer derived from the deterministic Riccati observer design framework recently proposed in [START_REF] Hamel | Riccati observers for the nonstationary PnP problem[END_REF]. We believe that the resulting estimated quantities would be less noisy since the noise can be filtered through a natural lowpass response of the observer. Moreover, the challenging theoretical issue related to the uniform observability, under which local exponential stability is granted, has been carefully addressed, and simulation results demonstrate a good performance and a large attraction domain of the proposed observer.

II. THEORETICAL BACKGROUND

A. Notation

• {e 1 , e 2 , e 3 } denotes the canonical basis of R 3 and the identity matrix of R n×n is denoted as I n .

• (•) × denotes the skew-symmetric matrix associated with the cross product, i.e.

u × v = u × v, ∀u, v ∈ R 3 . • π u I 3 -uu , ∀u ∈ S 2 ,
is the projection operator onto the plane orthogonal to u.

B. Brief overview of the Riccati observer design framework

The observer proposed in this paper is based on the Riccati design framework recently developed in [START_REF] Hamel | Riccati observers for the nonstationary PnP problem[END_REF]. The following nonlinear system (which is a particular case of the class of systems considered in [START_REF] Hamel | Riccati observers for the nonstationary PnP problem[END_REF]) is investigated:

Ẋ = A(t)X + U + O(|X| 2 ) + O(|X||U |) Y = C(X, t)X + O(|X| 2 ) (1) with state X = [X 1 , X 2 ] , X 1 ∈ B n1 r (the closed ball in R n1 of radius r), X 2 ∈ R n2 , output Y ∈ R m , A (
t) a continuous matrix-valued function uniformly bounded with respect to (w.r.t.) t in the form

A(t) = A 1,1 (t) 0 n1×n2 A 2,1 (t) A 2,2 (t) (2)
and C(X, t) ∈ R m×(n1+n2) a continuous matrix-valued function uniformly continuous w.r.t. X and uniformly bounded w.r.t. t. Then, apply the input

U = -P C DY (3)
with P ∈ R (n1+n2)×(n1+n2) a symmetric positive definite matrix solution to the following continuous Riccati equation (CRE):

Ṗ = AP + P A -P C D(t)CP + S(t) (4)
with P (0) ∈ R (n1+n2)×(n1+n2) a symmetric positive definite matrix, D(t) ∈ R m×m bounded continuous symmetric positive semi-definite, and S(t) ∈ R (n1+n2)×(n1+n2) bounded continuous symmetric positive definite. Then, from Theorem 3.1 and Corollary 3.2 in [START_REF] Hamel | Riccati observers for the nonstationary PnP problem[END_REF], X = 0 is locally exponentially stable (LES) when both matrices D(t) and S(t) are larger than some constant positive matrix and the pair (A (t), C (t)), with A (t) A(t), C (t) C(0, t), is uniformly observable.

C. Recalls of uniform observability

Definition 1 (uniform observability) The pair (A(t), C(t)) is called uniformly observable if there exist δ, µ > 0 such that ∀t ≥ 0

W (t, t + δ) 1 δ t+δ t Φ (t, τ )C (τ )C(τ )Φ(t, τ )dτ ≥ µI n
(5) with Φ(t, τ ) the transition matrix associated with A(t), i.e. such that d dt Φ(t, τ ) = A(t)Φ(t, τ ) with Φ(t, t) = I n . W (t, t + δ) is called the observability Gramian. The following lemma [START_REF] Scandaroli | Visuo-inertial data fusion for pose estimation and self-calibration[END_REF] establishes a sufficient condition of uniform observability.

Lemma 1 (see [START_REF] Scandaroli | Visuo-inertial data fusion for pose estimation and self-calibration[END_REF]) If there exists a matrix-valued function

M (•) of dimension (p × n) (p ≥ 1) composed of row vectors of N 0 = C, N k = N k-1 A + Ṅk-1 , k = 1, • • • such that for some positive numbers δ, μ and ∀t ≥ 0 1 δ t+ δ t det M (τ )M (τ ) dτ ≥ μ (6)
then W (t, t + δ) satisfies condition (5).

III. PROBLEM STATEMENT

A. Homography decomposition problem

Consider a moving camera observing a textured planar scene. Let {A} (resp. {B}) denote the reference (resp. current) right-hand frame attached to the camera at the reference (resp. current) view. Let ξ ∈ R 3 denote the position of frame {B} w.r.t. frame {A} expressed in frame {A}. The orientation of frame {B} w.r.t. frame {A} is represented by a rotation matrix R ∈ SO(3). The pose of the camera (R, ξ) determines a rigid body transformation from {B} to {A}.

The so-called Euclidean homography that maps Euclidean coordinates of the scene's points from {B} to {A} is given by (see Fig. 1) where d and η are, respectively, the normal distance from the origin of frame {B} to the observed scene and the normal vector pointing towards the scene expressed in {B}.

H = R + 1 d ξη (7)
Assume that the camera is well calibrated, then the Euclidean homography can be directly computed from the so-called projective homography or image homography estimated from image point correspondences (see [START_REF] Malis | Deeper understanding of the homography decomposition for vision-based control[END_REF] for details).

The well-known problem of Euclidean homography decomposition consists in decomposing the matrix H into the elements R, ξ d and η. Note that expression ( 7) of H can be rewritten as

H = R(I 3 + ξη ) (8)
with

ξ := R ξ d (9)
Since ξ d can be computed as R ξ, decomposing H into (R, ξ d , η) is equivalent to decomposing it into (R, ξ, η). In this paper, we propose to explore a novel direction of solving this homography decomposition problem dynamically by exploiting the differential equations guiding these variables (i.e. R, ξ, η). The underlying idea is to design an observer allowing for an exponential convergence rate of the state estimates ( R, ξ, η) to the real state variables (R, ξ, η).

B. System equations and measurements for observer design

The rigid body kinematics of (R, ξ) are given by

Ṙ = RΩ × ξ = RV (10)
with V ∈ R 3 and Ω ∈ R 3 denoting the linear and angular velocities expressed in {B}.

Since the scene is stationary the normal vector η ∈ S 2 expressed in the reference frame {A} is constant. Since η = R η one thus deduces

η = -Ω × η (11)
Using ( 9) and (10), one verifies that

ξ = -Ω × ξ + V d - ḋ d ξ = (-Ω × + φ ⊥ I 3 ) ξ + φ (12)
where φ := V d and φ ⊥ := -ḋ d = V η d are the so-called translational optical flow and optical flow divergence, respectively. Both the translational optical flow and the optical flow divergence are assumed to be measured (see our prior work on how these quantities can be estimated [START_REF] Manerikar | Homography observer design on special linear group SL(3) with application to optical flow estimation[END_REF]). Assume also that the angular velocity Ω is also measured using an embedded 3-axes gyrometers.

In this paper we focus on developing a nonlinear observer allowing for dynamically decomposing the homography matrix H into (R, ξ, η) by exploiting the system dynamics (10)-(12), i.e.

     η = -Ω × η Ṙ = RΩ × ξ = (-Ω × + φ ⊥ I 3 ) ξ + φ (13)
and using the homography measurement equation ( 8) together with the measurements of (Ω, φ, φ ⊥ ).

IV. OBSERVER DESIGN

A. Observer derivation

Inspired by [START_REF] Hua | Riccati observer design for pose, linear velocity and gravity direction estimation using landmark position and IMU measurements[END_REF], we avoid using minimal parametrization techniques such as spherical coordinate system to parametrize the normal unit vector η -an element of S 2 . Instead, an auxiliary rotation matrix

Q ∈ SO(3) is intro- duced such that η = Q e 3
The underlying idea is to over-parameterize an element of S 2 (dimension 2) by an element of SO(3) (dimension 3).

The advantage of such type of parameterization is that it reduces the complexities of the error system in first order approximations that arise due to minimal parameterization techniques for elements on S 2 (see [START_REF] Hua | Riccati observer design for pose, linear velocity and gravity direction estimation using landmark position and IMU measurements[END_REF] for more thorough discussions).

In view of the dynamics (11) of η, one deduces a possibility of the dynamics of Q as

Q = QΩ × (14)
In view of ( 13) and ( 14), the following general form of observer is proposed

         Q = QΩ × -σ Q× Q Ṙ = RΩ × -Rσ R× ξ = (-Ω × + φ ⊥ I 3 ) ξ + φ -σξ (15)
with initial conditions Q(0), R(0) ∈ SO(3), ξ(0) ∈ R 3 and with innovation terms σ Q , σ R , σξ ∈ R 3 to be designed thereafter. The estimated normal vector is then given by

η := Q e 3
The following error variables are defined:

Q := Q Q , R := R R, ξ := ξ -ξ (16)
Then the objective of observer design consists in stabilizing ( Qe 3 , R, ξ) about (e 3 , I 3 , 0).

From (10), ( 15) and ( 16), one verifies that the error system is given by

         Q = Qσ Q× Ṙ = RΩ × -Ω × R + σ R× R ξ = (-Ω × + φ ⊥ I 3 ) ξ + σξ (17)
For analysis purposes let us assume that ξ, Ω and φ ⊥ remain bounded for all time, which is a completely reasonable assumption.

The following step involves developing first order approximations of the error system (17) and of the measurement equation [START_REF] Scandaroli | Visuo-inertial data fusion for pose estimation and self-calibration[END_REF]. From the Rodrigues' formula, one deduces the following first order approximations of Q and R

Q = I + λ Q× + O(|λ Q| 2 ) R = I + λ R× + O(|λ R| 2 ) ( 18
)
with λ Q, λ R ∈ S 2 . One then deduces from the first two equations of ( 17) and ( 18) that in first order approximations

λ Q = σ Q + O(|λ Q||σ Q |) (19)
and

λ R× = λ R× Ω × -Ω × λ R× +σ R× +O(|λ R| 2 )+O(|λ R||σ R |) = (λ R × Ω) × + σ R× + O(|λ R| 2 ) + O(|λ R||σ R |) which yields λ R = -Ω × λ R + σ R + O(|λ R| 2 ) + O(|λ R||σ R |) (20)
As for the measurement equation ( 8), this homography expression can be developed in first order approximations as

R H -I 3 = ξη = ξ(Q e 3 ) ⇒ R R H -I 3 = ( ξ + ξ)( Q Q e 3 ) ⇒ (I 3 -λ R× ) R H -I 3 = ( ξ + ξ)( Q (I 3 -λ Q× )e 3 ) + O(|λ Q| 2 ) + O(|λ R| 2 ) ⇒( R H -I 3 )-λ R× ( R H) = ξ( Q e 3 ) + ξ( Q e 3 ) + ξ( Q e 3× λ Q) +O(|λ Q| 2 ) + O(|λ R| 2 ) ⇒ ( R H -I 3 ) -ξ( Q e 3 ) = λ R× ( R H)+ ξ( Q e 3 ) -λ Q,2 ξ( Q e 1 ) +λ Q,1 ξ( Q e 2 ) +O(|λ Q| 2 ) + O(|λ R| 2 )
Note that the last equality only involves the first two components of λ Q (i.e. λ Q,1 and λ Q,2 ) and can be equivalently written as

Y :=    ( R H -I 3 )( Q e 3 ) - ξ ( R H -I 3 )( Q e 2 ) ( R H -I 3 )( Q e 1 )    =     -( R H)( Q e 3 ) × λ R + ξ -( R H)( Q e 2 ) × λ R + λ Q,1 ξ -( R H)( Q e 1 ) × λ R -λ Q,2 ξ    +O(|λ Q| 2 )+O(|λ R| 2 ) (21)
From ( 19), (20), the third equation of (17), and (21), one obtains in first order approximations the system in compact form [START_REF] Faugeras | and structure from motion in a piecewise planar environment[END_REF] with output Y defined in (21) and

                                         X :=     λ Q,1 λ Q,2 λ R ξ     , U :=     σ Q,1 σ Q,2 σ R σξ     A(t) :=     0 0 0 1×3 0 1×3 0 0 0 1×3 0 1×3 0 3×1 0 3×1 -Ω × 0 3×3 0 3×1 0 3×1 0 3×3 -Ω × + φ ⊥ I 3     C(X, t) :=    0 3×1 0 3×1 -( R H)( Q e 3 ) × I 3 ξ 0 3×1 -( R H)( Q e 2 ) × 0 3×3 0 3×1 -ξ -( R H)( Q e 1 ) × 0 3×3    (22)
The fact that

X = [X 1 , X 2 ] with X 1 := [λ Q,1 , λ Q,2 , λ R ] ∈ B 5
2 (the closed ball in R 5 of radius equal to 2), X 2 := ξ ∈ R 3 , together with the particular form of the matrix A as (2), allows one to obtain the expression of the innovation terms from the input U calculated according to ( 3) and ( 4) where the matrices D and S (involved in ( 4)) are chosen larger than some constant positive matrix.

B. Observability and stability analysis

According to [START_REF] Hamel | Riccati observers for the nonstationary PnP problem[END_REF] the equilibrium X = 0 is locally exponentially stable, provided that the pair (A (t), C (t)) with A (t) = A(t) and C (t) := C(0, t) is uniformly observable. By setting X = 0 and looking at (22) one obtains

C =    0 3×1 0 3×1 -(I 3 + ξη )(Q e 3 ) × I 3 ξ 0 3×1 -(I 3 + ξη )(Q e 2 ) × 0 3×3 0 3×1 -ξ -(I 3 + ξη )(Q e 1 ) × 0 3×3    =    0 3×1 0 3×1 -η + ξ × I 3 ξ 0 3×1 -Q e 2 × 0 3×3 0 3×1 -ξ -Q e 1 × 0 3×3    with Q ∈ SO(3) satisfying Q e 3 = η and Q = Q Ω × .
For later use, let

q i ∈ S 2 , i = 1, 2, 3, denote the i-th row of Q (i.e. q i = Q e i ).
Theorem 1 Assume that there exists a positive number ν such that ∀t > 0 1

δ t+δ t | ξ(τ ) × η(τ )|dτ ≥ ν (23)
Assume also that Ω, φ, ξ remain uniformly bounded. Then, the pair (A * , C * ) is uniformly observable. Assume that the matrices D and S involved in (4) are chosen larger than some constant positive matrix. Then, the equilibrium ( Qe 3 , R, ξ) = (e 3 , I 3 , 0) of the error system is locally exponentially stable. Proof: According to Lemma 1 by choosing M = C , then the pair

(A * , C * ) is uniformly observable if ∃δ, µ > 0 such that 1 δ t+δ t det C (τ )C (τ ) dτ ≥ µ, ∀t > 0 (24)
We show thereafter that condition ( 23) is sufficient to guarantee (24). One verifies that

C C =     | ξ| 2 0 -ξ q 2× 0 1×3 0 | ξ| 2 ξ q 1× 0 1×3 q 2× ξ -q 1× ξ -η+ ξ 2 × -q 2 1× -q 2 2× η + ξ × 0 3×1 0 3×1 -η+ ξ × I 3    
which can be rewritten as

C C = | ξ| 2 I 2 1 | ξ| 2 B 1 | ξ| 2 B 1 | ξ| 2 G
where the expressions of B ∈ R 2×6 and G ∈ R 6×6 can be easily deduced. Thus, one deduces

det C C = | ξ| 4 det G - 1 | ξ| 2 B B = | ξ| 4 det -η+ ξ 2 × -q 1× π ξ | ξ| q 1× -q 2× π ξ | ξ| q 2× η + ξ × -η + ξ × I 3 = | ξ| 4 det -q 1× π ξ | ξ| q 1× -q 2× π ξ | ξ| q 2× = | ξ| 4 det Q -e 1× πQ ξ | ξ| e 1× -e 2× πQ ξ | ξ| e 2× Q = | ξ| 4 det -e 1× π η ξ e 1× -e 2× π η ξ e 2× = | ξ| 4 det e 3 e 3 + η ξ η ξ + (η ξ × e 3 )(η ξ × e 3 ) = | ξ| 4 |η ξ × e 3 | 4 = | ξ × η| 4 with η ξ := Q ξ | ξ| .
From the Cauchy-Schwarz integral inequality, one deduces (using (23))

t+δ t | ξ(τ ) × η(τ )| 4 dτ ≥ 1 δ t+δ t | ξ(τ ) × η(τ )| 2 dτ 2 ≥ 1 δ 3 t+δ t | ξ(τ ) × η(τ )|dτ 4 ≥ δν 4
One finally deduces (24) with µ = ν 4 . The remainder of the proof then directly follows by application of Theorem 3.1 and Corollary 3.2 in [START_REF] Hamel | Riccati observers for the nonstationary PnP problem[END_REF].

V. SIMULATION RESULTS

In this section the robustness and performance of the proposed observer is demonstrated through simulation results by considering two scenarios in particular. For the simulations, we consider that a monocular camera is attached to an aerial drone performing some specific trajectories and observing a planar target. In the first scenario, the drone realizes an aggressive periodical trajectory and noise is introduced to the measurements used in the proposed observer, namely H, φ , φ ⊥ and Ω. While in the second scenario, we consider the case where the camera trajectory passes through the reference position (i.e. the camera translation vanishes) again with the introduction of noise in the measurements.

With Matlab Simulink, we try to simulate the image noise by introducing white Gaussian noise of the level of about 10 percent of the real values on each individual component of the homography matrix H. We also introduce white The time evolutions of the estimated and real attitudes (represented by Euler angles), the attitude error estimate (represented by trace(I -R)) as well as the scaled position error estimate and the estimation error of the normal vector to the planar scene (represented by 1 -η η) are shown in Figs 234, respectively. From these figures, it can be clearly seen that all the estimated variables converge to the real ones after a short transition period of a few miliseconds despite the large initial estimation errors. The above results also show that the proposed observer is robust to noisy measurements as the latter marginally affects the overall performance of the proposed observer. Scenario 2: In this case the drone performs a rectilinear sinusoidal trajectory that passes through the reference position (i.e. the camera translation vanishes) with added noise on the measurements of H, φ, φ ⊥ , and Ω as in the previous simulation. The reference trajectory in the inertial frame is given by ξ = [5sin(πt/3), 0, 0] (m). It is worth noting that in such a degenerate case (i.e. when the camera translation vanishes) all traditional algebraic approaches fail to obtain a correct estimate of the normal vector (there exists in fact an infinity of solutions) whereas our approach still works since the uniform observability condition (23) is always satisfied. This is also because of the fact that the state estimates are propagated over time. From Figs. 5-7 it can be clearly seen that for this degenerate case and even with the addition of noise to the system measurements, the proposed observer is robust enough to provide a convincing performance in terms of convergence rate, smooth transient phase, and filtering of measurement noise.

VI. CONCLUSIONS

In this paper a novel approach for solving the classical problem of homography decomposition has been proposed by designing a nonlinear Riccati observer. To our knowledge this is the first study to address the homography decomposition problem by exploiting the system dynamics. By exploiting the temporal correlation of image sequences and the natural low-pass response of the observer, the resulting estimates are robust and less prone to measurement noise. Exhaustive observability and stability analysis has been carried out in order to support the proposed observer. At the end convincing simulation results have been presented to show the performance of the observer.

Fig. 1 .

 1 Fig. 1. Euclidean homography relating the camera's pose, the distance to the plane and the plane's normal vector by H = R + 1 d ξη .

Fig. 2 . 1 -Fig. 3 . 1 -Fig. 4 . 1 -

 213141 Fig. 2. Scenario 1 -Estimated and real attitudes represented by roll, pitch and yaw Euler angles (deg) versus time (s) (LEFT) and zooms (RIGHT)

Fig. 5 . 2 -

 52 Fig. 5. Scenario 2 -Estimated and real attitudes represented by roll, pitch and yaw Euler angles (deg) versus time (s) (LEFT) and zooms (RIGHT)

Fig. 6 .Fig. 7 . 2 -

 672 Fig. 6. Scenario 2 -Estimated and real scaled position versus time(s) (LEFT) and zooms (RIGHT)

Acknowledgment: This work was supported by the ANR Astrid CONGRE (ANR-18-ASTR-0006), FUI Green-Explorer and ANR EQUIPEX Robotex projects.