Jerome Villeneuve
email: jerome.villeneuve@imag.fr

Claude Cadoz
email: claude.cadoz@imag.fr

Understanding and Tuning Mass-Interaction Networks Through Their Modal Representation

Sound is all about vibration, and the GENESIS environment provides an efficient way for modeling and simulating complex vibrating structures, enabling to produce rich sounds.

In this paper, we propose an overview of tools recently developed and available within the GENESIS environment, allowing a better understanding on how massinteraction networks behave and introducing some enhanced tuning of their vibrating properties. All these tools try to address an inherent need of any creative process either in the physical world or in GENESIS, which is to create bidirectional connections between properties of a phenomenon, in our case, audible sounds, and properties of what produced it, here, mass-interaction networks.

For this purpose, we will introduce the topological and modal representations of such mass-interaction networks and appreciate how relevant it can be to switch between these different representations to really apprehend its inner properties and those of the sounds it produces.

INTRODUCTION

GENESIS [START_REF] Castagne | G3: GENESIS software environment update[END_REF] is a musician-oriented software environment for sound synthesis and musical composition. It implements the modular language CORDIS-ANIMA [START_REF] Cadoz | CORDIS-ANIMA : a modeling and simulation system for sound and image synthesis[END_REF] that provides a set of primitive elements, modules satisfying Newtonian physics.

The elementary bricks of CORDIS-ANIMA are distinguished in two categories: the <MAT> modules that designate "matter" being either moving or fixed and the <LIA> modules that designate interactions and allow interconnections between <MAT> modules. A GENESIS "topological" model is then constructed by connecting these modules to each other into a dedicated 2D space called "the bench", this by direct manipulation (Figure 1). Such a model will be defined by its "topology", meaning its network structuration, by the parameters of each one of the modules it is made of (M for inertia of MAT modules, K for stiffness, Z for damping for LIA modules), and finally by a set of initial conditions carried by MAT modules (position and/or velocity).

Once built and parameterized, the model will be set in motion through an outer excitation, and its movement will be « listened to » (Figure 1). Mass-interaction network of a "string" (indeed, its geometry is quite unusual, but it does not matter. Its topology does). Yellow dots are moving <MAT>, green ones stay fixed.

From there, one can build simple or very large massinteractions networks, composed with thousands of modules. Once simulated, these allow obtaining complex and rich sounds. But whatever the size of such models, when it comes to exploring the variety of sounds they are able to produce, or even to adjust their nature or structure with the aim of fine-tuning this sounds, a major problem arises. Its widest formulation expressing a need to establish two-ways relations between causes and effects. There are several complementary ways to go on with such concerns.

The first, empirical, is to build one's own knowledge of how objects, models, behave according to how we act on or modify them. One will have to learn how slight adjustment of model parameters will alter perceptive properties of the product of its simulation. Furthermore, in the case of the GENESIS environment, it is now possible to physically experiment virtual models by means of haptic interfaces [START_REF] Leonard | Interactive Physical Design and Haptic Playing of Virtual Musical Instruments[END_REF], which allow, by an additional gestural feedback, a very direct and instrumental way to investigate the properties and possibilities of models in terms of musical creation. In both cases, a closed loop is set that involves a model and its simulation, the perceptive outcomes of this simulation, and the adjustment of action or modification applied on this model.

A second one is to get to have a better understanding of the model and the sound it produces by using dedicated analytic technics and tools. These can be derived and adapted from the physical world to the virtual world proposed by GENESIS, or having a more mathematical essence. In our case, we will have to understand the matter at its very smallest scale, so the properties of large objects made of this matter can be easily deducted. For that we will introduce in the following a complementary representation to the topological mass-interaction networks. Based on both representations, and by switching between them, we will present some tools now implemented in the GENESIS environment and allowing to help users in their creative process.

MODEL ANALYSIS

Real objects, made out of real matter, and played as musical instruments are since quite recent times under the extensive survey of acoustic, fundamentally motivated in understanding how they produce sound [START_REF] Fletcher | The Physics of Musical Instruments[END_REF]. But even if the mechanisms of some musical instruments are still nowadays freshly discovered, instrument-makers haven't waited to develop tremendous know-how in their conception and neither have the musicians to learn how to play and to take the best out of them.

Sound synthesis by the mean of mass-interaction networks modeling and simulation is one of a few methods that, before producing sounds, require to handle matter and to build complex objects out of it. A virtual matter for virtual objects, of course, but still existing in itself, answering to its own "physics" and that can be studied as it is. And it is precisely by having this objects that one can learn a lot about the sound that it might ultimately produce. This, by developing and using a dedicated "acoustic" of mass-interaction networks

The concept is there to extract, directly by studying the properties of a model, its structure and parameters, without neither simulation nor signal analysis, relevant information regarding to the sound phenomena it will allow to produce. It is basically a transposition of what the Fourier analysis is for the sound signal to the model that enables its simulation.

The model modal representation

The model modal analysis relies on a mathematical description of the "mechanical network" that is intrinsically every topological CORDIS-ANIMA model. This mechanical network must be a linear vibrating structure moving with only one degree of freedom. Then, the analysis process first aims at switching in representation from a topological model to its equivalent modal model [START_REF] Djoharian | Generating Models for Modal Synthesis[END_REF]. The latter can be seen as a set of damped harmonic oscillators with unitary inertias (counting as much oscillators than the number of punctual moving <MAT> of the original topological models).

These oscillators represent the stationary vibration modes of the original topological structure. Ultimately, a proper study of the latters will allow a complete knowledge of the original model in terms of mechanical behaviour and acoustical properties. The following development applies to free vibrating structures without any outer perturbation.

From topological model to modal model

As a network, a linear topological CORDIS-ANIMA model can be mathematically described by the combination of 3 real symmetric matrices: M , K and Z , carrying respectively inertia, stiffness and damping parameters of its constitutive modules.

Obtaining the equivalent modal model of a topological one is done by diagonalizing two real symmetric matrices K ' et Z' in the same basis.

Z' = M -1/2 ZM 1/2 & K ' = M -1/2 KM 1/2 (1)
This is achievable only if the matrices K ' and Z' commute (ie : K'Z' = Z'K' ), which is always the case if they are linearly co-dependent [START_REF] Zhang | Matrix Theory: Basic Results and Techniques[END_REF], meaning that Z' is proportional to K ' :

Z' = a + bK ' or Z = aM + bK (2)
If it answers to relation (2), the viscosity matrix is then said "proportional". The diagonalization process will allow to obtain: two diagonal matrices K m and Z m , re- spectively containing n modal stiffness k i m and n modal damping z i m , and an additional Q m matrix such as :

Z m = Q m Z t Q m and K m = Q m K t Q m (3)
Mathematically, k i m and z i m are the respective eigenvalues of K ' and Z' , and the Q m matrix, representing the common basis of K m and Z m , contains their associ- ated eigenvectors (directly interpretable as modal shapes of the equivalent topological model).

The n harmonic oscillators of the modal model are then directly defined by unitary inertias, k i m and z i m as parameters.

Reading a model modal representation

Some of the mathematical tools that are used for describing and studying vibratory and acoustical properties of physical matter are adaptable to the virtual matter handled in GENESIS. Nevertheless, CORDIS-ANIMA relies on a discrete representation of time and matter, which implies well-quantified differences between continuous and discreet equations of movement.

Having the modal model, it is now possible to simply analyze each of its individual elementary harmonics oscillators [START_REF] Incerti | Time-Frequency Decomposition by Means of Physical Model[END_REF]. They are defined by 3 parameters, m (=1), k i m and z i m in direct relation with:

The vibration mode frequency:

f i = F e ⋅ arccos 2 -k i m + z i m ( ) 2 ⋅ 1-z i m # $ % % & ' ( ( (4)
Where F e is the sampling rate value of signals produced when simulating a GENESIS model (default value: 44100 Hz)

The vibration mode damping time:

1

τ i = - F e 2 ⋅ ln 1-z i m ( ) (5)
We can also learn about relative amplitudes of the different oscillators by a proper reading of the transfer matrix Q m . By choosing the topological structure moving <MAT> on which a listening module will be plugged (identified as !), and the topological model moving <MAT> that will be submitted to an impulsive excitation (identified as !). Then we can estimate the relative amplitude of each vibration mode impulse response according to:

A i = Q m (α ,i ) ⋅ Q m (β ,i ) (6)
The theory beyond complex excitations and sustained oscillations of mass-interaction networks is already greatly documented [START_REF] Poyer | Self-sustained vibrating structures physical modeling by means of massinteraction networks[END_REF]. From the modal representation perspective, calculation of each independent harmonic oscillator transfer function enables a first insight of vibrating model properties under such constraint. For instance, resonant frequencies can be determined and a whole frequency response of a topological model can be displayed.

Model analysis implementation

A specific tool for model modal analysis has been developed in GENESIS [START_REF] Castagne | GENESIS: A Friendly Musician-Oriented Environment for Mass-Interaction Physical Modelling[END_REF] and has now been updated. Follows a quick tour of its latest improvements and functionalities.

The diagonalization process involved in switching from topological to modal representations is processed by using the Jacobi method [START_REF] Lascaux | Analyse numérique appliquée à l'art de l'ingénieur[END_REF][START_REF] Saas | Numerical Methods for Large Eigenvalue Problems[END_REF], which allows to obtain simultaneously and efficiently both eigenmodes and eigenvalues of a given matrix. This operation is quite costly but enable to compute the modal representation and acoustical properties of models counting more than 9000 moving <MAT> modules.

Regarding the proportionality constraint over the viscosity matrix (2), the modal analysis can be run even if the model is not proportional. However, we chose to implement a simple algorithm testing this proportionality prior to the diagonalization process, which alerts the user about the relative quality of the obtained result.

One can select a model or a portion of it directly from the bench and launch the Modal Analysis function. The displaying window of the result is organized in 3 distinct parts (Figure 3), which are:

• A general sheet section, gathering the whole information relative to each identified vibration mode of the structure, and by extension, to each partial that will compose a sound produced out of it. It also provides an explicit conversion from frequency to equivalent musical note, octave and cents (Figure 4). • A modal shapes visualization section, that allows to individually display an interactive 3D representation of each topological structure modal shapes. The latter can be animated and rendered with various static or dynamic color gradations. All this allowing for example to clearly identify nodes and antinodes of the vibrating structure, and ultimately determine where to act on a structure in accordance to modes/partials that should or should not be predominant (Figure 5). • A normalized impulse response graph section, providing a global picture of the whole frequency spectrum associated with relative amplitudes of model vibration modes. This supposes that both in and out points, respectively, excitation and listening points, are selected directly on the model displayed in the modal shape visualization section (Figure 6).

MODEL TUNING

Having a mass-interaction network and looking to have a better understanding of it is one thing. Trying to tune it, or in the most extreme case, having a "sound" in mind and trying to reproduce it by the mean of massinteraction networks suggest an other level of considerations. With the idea to provide tools to help the users, and to keep close to the illustrative concepts of the model modal analysis, we have tried to explore the way back from a modal model to its topological equivalent. First thing first, lets assume that one has built a topological model. How can they tune the latter by modifying properties of its modal equivalent?

From modal back to topological

In the very limited case of a topological model of which the modal equivalent carry a single harmonic oscillator, we can easily reverse the linear equations ( 3) and ( 4) and then tune this oscillator by recovering a proper set of stiffness and viscosity parameters with frequency and damping time as inputs.

Considering a modal model with more than one harmonic oscillator, it wouldn't be possible to tune each one of them independently because of K m and Z m that must be preserved for the relation (3) to stay true (necessary condition to get back from a modal representation to the original topological representation). Nevertheless it is possible to apply a single multiplicative coefficient to each of these matrices. This authorizes to set one vibration mode frequency and damping time. Of course all the other structure vibration modes will be proportionally tuned. This function exists in GENESIS and is built directly into the modal analysis window. Hence, from the general sheet section, one can simply edit the acoustical properties (frequency, or related musical note) of one vibration mode of the modal representation and tune the topological model.

Relative amplitudes of modes depend on the model properties and for the most part on how/where it is excited and listened. There is no trivial function enabling to tune a model in a way that the sound it produces has a precise amplitude distribution among its partials. For that, one will rely on the modal shapes visualization and on the impulsion response graph of the Modal Analysis window, which might help in finding the best set of input(s)/output(s).

It is not possible to have an extensive control over the physical parameters of a model, according to acoustical properties, without altering the equivalence/bridge standing between its modal and topological representations. But by breaking this relation that allows to keep causes and effects close together, advanced tuning possibilities can be foreseen.

From modal back to …?

Indeed, at a certain point, one could seek to obtain a GENESIS model given nothing else that the description of a sound they want it to produce.

A simple solution to that concern would be to recompose complex sounds by reproducing what's done in additive sound synthesis or modal synthesis [START_REF] Adrien | The Missing Link: Modal Synthesis[END_REF], which as already been specifically studied for CORDIS-ANIMA applications in [START_REF] Djoharian | Generating Models for Modal Synthesis[END_REF]: Knowing the modal characterization of a sound phenomenon, we could compose it by synthesizing as many sinusoidal signals as necessary, this by replicating a relevant amount of independents harmonic oscillators. The main issue of this approach is that one will have to artificially selectively distribute energy among the independent oscillators and selectively "listen" to each one of them. Thus, it prevents to obtain a topological representation of the proposed modal model and might considerably limit users in theirs exploration process.

In [START_REF] Kontogeorgakopoulos | Filtering within the Framework of Mass-Interaction Physical Modeling and of Haptic Gestural Interaction[END_REF], the authors already explore ways to tune a specific model under frequency constraints. Nevertheless a more general approach, aiming at obtaining a topological mass-interaction network out of multiple acoustical descriptors of sounds, is addressed in [START_REF] Villeneuve | Inverse Problem in Sound Synthesis using Mass-Interaction Networks[END_REF] and is referred as an "inverse problem". Follows the outcomes of preliminary investigations in simplified cases and their effective implementation in the GENESIS environment.

An inverse problem

A model of n moving <MAT> has an equivalent modal model of n independent oscillators and thus n eigenmodes. If we now consider the description of a sound composed of n partials, setting a modal model of n moving <MAT> as an input, then a generative topological model must have at least n moving <MAT>, and we can suppose it is enough to define it. Among the multiple solutions that this problem can have, it is necessary to arbitrary choose an alleged generative mass-interaction network. The simplest topological model of n moving <MAT> is the linear "chain", a "string" fixed at its extremities (such as in Figure 7), at one extremity or none. The reduced inverse problem is then to define, given 3 vectors of n frequencies, amplitudes and damping times, a generative mass-interaction string of n moving <MAT> connected by visco-elastic interactions. Different approaches have been followed to resolve it, all basically trying to recover the descripting matrices of the topological (M, K, Z) model given diagonal matrices of the tuned modal model ( K m and Z m ):

The first one, a numerical resolution, relies on an adapted optimization procedure. It partially addresses the reduced problem by tuning only the stiffness parameters of the generative model, allowing an individual control of each structure vibration mode frequencies (in this case, the moving <MAS> modules are set with unitary inertias).

The second is algebraic and allows to control frequencies as well as relative amplitudes of each vibration mode of the generative structure. Hence, the stiffness and inertia parameters of the latter are tuned. Unlike the numerical approach, this one is case specific and only allows to generate "string" models.

In this case the relative amplitudes are set by direct manipulation of the modal shapes matrix Q m . Further details about the previous can be found in [START_REF] Villeneuve | Inverse Problem in Sound Synthesis using Mass-Interaction Networks[END_REF] A third approach has been dedicated to enable a specific control over vibration modes damping time. We previously discussed the term of "proportional viscosity matrix" and expressed it with relation [START_REF] Cadoz | CORDIS-ANIMA : a modeling and simulation system for sound and image synthesis[END_REF]. Going from a modal representation to a topological one implies to stick close with this rule. Hence, having the inertia and stiffness parameters of a topological model already tuned regarding frequencies and relative amplitudes of its vibration modes, gives a pretty restricted control over the relation [START_REF] Cadoz | CORDIS-ANIMA : a modeling and simulation system for sound and image synthesis[END_REF]. Indeed, only two multiplicative values, a and b, can be edited, meaning that only two vibration modes damping time can be adjusted (the other ones will be proportionally balanced).

Actually, these two values are physically coherent and are a proper metaphor of external and internal viscosities of a mechanical object. a, applied to the inertias matrix, can be liken to an external viscosity and represent the dissipative effect of viscous environment in which a model is immersed. b, applied to the stiffness matrix, could represent an internal viscosity and stand for the dissipative effect inherent to the model internal physical matter. Thus, by renaming respectively a and b to z int and z ext in (2) we have a direct expression of the topological model viscosity matrix:

Z = z ext M + z int K (7)
Setting the damping-time of the i th and j th modes, with i < j, implies:

τ i ≥ τ j (8)
The internal and external viscosities constants are then defined by:

z int = exp - 2 τ i F e ! " # # $ % & & -exp - 2 τ j F e ! " # # $ % & & K ( j, j) m -K (i,i) m ( 9
)
z ext = K ( j, j) m 1-exp - 2 τ i F e ! " # # $ % & & ! " # # $ % & & -K (i,i) m 1-exp - 2 τ j F e ! " # # $ % & & ! " # # $ % & & K ( j, j) m -K (i,i) m ( 10
)
It is interesting to point out that the relation ( 8) guaranties positive internal viscosity values. On another hand, external viscosities might take negative values, which metaphorically would imply that the model environment provides energy instead of absorbing it. Also, viscosities do have an influence on the vibration mode frequencies, but, they are quite insignificant, and of course this influence has to be evaluated regarding perceptual matters.

This approach works on every linear mass-interaction network, whatever its topology and complexity.

It doesn't seem possible to go further in controlling the damping-time of a modal model vibration modes without breaking the proportionality, altering the topology of the resulting model, or by staying physically coherent and thus being easy to handle by any user.

Model tuning implementation

PNSL Scripts

The different approaches previously proposed have been implemented thanks to an additional tool included within GENESIS: The Physics Networks Scripting Language (PNSL) [START_REF] Castagne | G3: GENESIS software environment update[END_REF]. This language has been developed to address mass-interaction networks modeling specific needs. For instance, by putting down some lines of script, it allows to create modules, to inter-connect them, to compute and set them with precise parameters, to automate repetitive tasks regarding topological or bench positions matters, to alter selected models, and so on. All this is fully integrated so one can edit and execute scripts directly into the GENESIS environment (Figure 8), and provides another modeling process, complementary to the bench direct manipulation of modules. One of the central considerations in building this tool and its syntax was to make it intelligible and easy to use, even for users non-acquainted with programming languages. Nevertheless, it is built over Tool Command Language (TCL) scripting language and therefore includes a lot of advanced libraries [START_REF] Ousterhout | Tcl and the Tk Toolkit[END_REF].

Inverse problems scripts

The numerical approach allows to efficiently obtain generative models defined by up to 20 frequencies. The time of resolution increases dramatically with the number of frequencies to deal with (in the case of 20 tuned vibration modes, the calculation might take about one hour). Its main advantage is that it doesn't depend on the kind of topology you may expect for the output model. It is also convenient to have all inertias set to 1 when bringing the generative model back to a common utilization in using GENESIS.

In its latest version, the algebraic approach is suitable for fixed, half-fixed or free linear chains. Any other topology for the output model would need to have its dedicated algorithm. However, with this approach we managed to obtain generative models defined by up to 50 frequencies and relative amplitudes. The global quality of this resolution, as well as the coherence of the resulting model, relies on the relative proximity of the wanted frequencies. For instance setting the modal model with very close frequencies (of about 1Hz) might introduce numerical imprecisions during the algorithm calculation and lead to partially tuned generative models.

Both of the previous methods can be executed in GENESIS and produce a topological network right on the bench. The damping-time tuning method needs the user to select a model directly on the bench. Running the script will then only tune its viscosity parameters.

With PNSL, a lot of new functions can be developed, easily used, and of course modified and enhanced, allowing easing prospective works in mass-interaction networks modeling. Inverse problem resolution scripts were the first ones to fully take advantage of this tool, and will be fully accessible as demonstrative scripts in the next update of GENESIS.

CONCLUSIONS

The mass-interaction networks modeling used in GENE-SIS rely on a very unique set of concepts. Its modularity makes it very generative and the exploration of what can be expected out of it seems far from its ending despite decades of developments, experimentations, and of course, of active utilization and artistic creations made by researchers and artists.

Luckily these concepts are close to the physics of real things, so on one hand, users can easily bring their physical instincts in the creative process involved in sound synthesis and musical composition and, on another hand, the tools developed with the aim of easing this creative process can themselves rely on physical metaphors and make direct sense to users.

Figure 1 .

 1 Figure1. Mass-interaction network of a "string" (indeed, its geometry is quite unusual, but it does not matter. Its topology does). Yellow dots are moving <MAT>, green ones stay fixed.

Figure 2 .

 2 Figure 2. Example of what can be an equivalent Modal Model of (Figure 1) Topological Model.

Figure 3 .

 3 Figure 3. Screenshot of the whole Modal Analysis Window.

Figure 4 .

 4 Figure 4. Modal Analysis Table, gathering modal model information. Purple columns collect the k i m and z i m modal parameters. The green columns display their associated acoustical properties.

Figure 5 .

 5 Figure 5. Collection of various topological models modal shapes. Left: we can identify the typical modal shapes of a circular membrane (here composed of 2148 moving elements). Middle: we represent some of the modal shapes of a "labyrinth" model (from which you can escape by looking closely at the first one). Right: various "punk smiley" model modal shapes.

Figure 6 .

 6 Figure 6. 3 impulse responses of a "string" model, composed of 40 moving <MAT> modules, excited and listened in different points. (a) Input: module 1, Output: module 40. (b) Input: module 4, Output : module 40. (c) Input: module 20, Output : module 20.

Figure 7 .

 7 Figure 7. CORDIS-ANIMA linear chain model.

Figure 8 .

 8 Figure 8. Screenshot of GENESIS PNSL Script Editor.

Proceedings ICMC|SMC|201414-20 September 2014, Athens, Greece


Jean-Matthieu Monnet
email: jean-matthieu.monnet@irstea.fr

Éric Mermin
email: eric.mermin@irstea.fr

Cross-Correlation of Diameter Measures for the Co-Registration of Forest Inventory Plots with Airborne Laser Scanning Data

Keywords: forest inventory, airborne laser scanning, co-registration, Global Navigation Satellite System

Continuous maps of forest parameters can be derived from airborne laser scanning (ALS) remote sensing data. A prediction model is calibrated between local point cloud statistics and forest parameters measured on field plots. Unfortunately, inaccurate positioning of field measures lead to a bad matching of forest measures with remote sensing data. The potential of using tree diameter and position measures in cross-correlation with ALS data to improve co-registration is evaluated. The influence of the correction on ALS models is assessed by comparing the accuracy of basal area prediction models calibrated or validated with or without the corrected positions. In a coniferous, uneven-aged forest with high density ALS data and low positioning precision, the algorithm co-registers 91% of plots within two meters from the operator location when at least the five largest trees are used in the analysis. The new coordinates slightly improve the prediction models and allow a better estimation of their accuracy. In a forest with various stand structures and species, lower ALS density and differential Global Navigation Satellite System measurements, position correction turns out to have only a limited impact on prediction models.

Introduction

Traditional forest inventory practices mainly rely on statistical descriptions obtained from field sample plots [START_REF] Duplat | L'inventaire par échantillonnage[END_REF]. This method provides reliable estimates at the forest level, while limiting fieldwork cost.

However, information at the management unit level (compartment) is often not usable, due to the low number of observations. In the past decade, airborne laser scanning (ALS) has demonstrated its ability to model the 3D geometrical structure of the canopy, even at the single tree level. The area-based approach aims at combining the high resolution of ALS data with the sample field plots in order to provide statistically calibrated, continuous maps of forest parameters [START_REF] Naesset | Predicting forest stand characteristics with airborne scanning laser using a practical two-stage procedure and field data[END_REF][START_REF] Hollaus | Operational wide-area stem volume estimation based on airborne laser scanning and national forest inventory data[END_REF]. Relationships between local descriptors of the ALS point cloud and forest parameters, such as basal area, stem volume or dominant height, are investigated based on the available field data. Once a prediction model is validated, it is applied to the whole ALS dataset in order to produce the map of the desired forest parameters.

Whereas the planimetric accuracy of an ALS point cloud is around 25 cm [START_REF] Baltsavias | Airborne laser scanning: Basic relations and formulas[END_REF], the position accuracy of field plots is lower. Indeed, plots are generally positioned relatively to elements that are visible on aerial pictures or with Global Navigation Satellite System (GNSS) receivers. Because of the canopy that reflects and attenuates the signal, the accuracy of measurements inside forest stands is lower than in open areas and of the order of a few meters [START_REF] Naesset | Assessing point accuracy of DGPS under forest canopy before data acquisition, in the field and after postprocessing[END_REF]. However, position accuracy is difficult to estimate, as it depends on acquisition parameters, satellite availability and canopy density. Moreover, precision values provided by post-processing software cannot be considered as reliable indicators [START_REF] Andersen | An accuracy assessment of positions obtained using survey-and recreational-grade global positioning system receivers across a range of forest conditions within the Tanana valley of interior Alaska[END_REF].

A proper co-registration is particularly important when field observations are linked to remote sensing data [START_REF] Korpela | Mapping forest plots: An efficient method combining photogrammetry and field triangulation[END_REF]. With ALS data, when the GNSS positions are shifted from the field positions, the predictive relationships are investigated between point clouds that do not exactly correspond to the inventoried sample plots. This error is likely to decrease the degree of fit of the final prediction models [START_REF] Gobakken | Assessing effects of positioning errors and sample plot size on biophysical stand properties derived from airborne laser scanner data[END_REF][START_REF] Frazer | Simulated impact of sample plot size and co-registration error on the accuracy and uncertainty of LiDAR-derived estimates of forest stand biomass[END_REF] and also result in a bad assessment of the prediction accuracy.

Several methods have been proposed to improve the co-registration of forest field plots with remote sensing data. Korpela et al. [START_REF] Korpela | Mapping forest plots: An efficient method combining photogrammetry and field triangulation[END_REF] used field triangulation and points derived from photogrammetry to improve the absolute mapping of large forest plots. Their method can be implemented by a single person with real-time estimation of positioning accuracy on the field. Hauglin et al. [START_REF] Hauglin | Geo-referencing forest field plots by co-registration of terrestrial and airborne laser scanning data[END_REF] co-register field plots by matching trees detected in ALS data with trees detected in terrestrial laser scanning data. When the relative tree position is recorded for the field plots and ALS with a sufficient pulse density is available, a visual comparison of the canopy height model (CHM) with the trees position and size is possible. However this task is time-consuming, in particular when the position error is large or when the forest canopy is closed and homogeneous. Some methods aim at automatically adjusting the tree position to ALS data. Olofsson et al. [START_REF] Olofsson | A method for linking field-surveyed and aerial-detected single trees using cross correlation of position images and the optimization of weighted tree list graphs[END_REF] first detected single trees in the CHM and then computed single trees position images with trees modeled as Gaussian surfaces. The image of the field trees is cross-correlated with the image of the detected trees in order to estimate the shift to be applied to the initial position. High co-registration accuracies are obtained, provided that the position error of single trees is small. This approach relies on single tree detection, which requires high density ALS data [START_REF] Hyyppä | A segmentation-based method to retrieve stem volume estimates from 3-D tree height models produced by laser scanners[END_REF]. Dorigo et al. [START_REF] Dorigo | An application-oriented automated approach for co-registration of forest inventory and airborne laser scanning data[END_REF] define a cost function that takes into account the difference in field and CHM heights, as well as the social status of a tree. Pascual et al. [START_REF] Pascual | Algorithm for improving the co-registration of LiDAR-derived digital canopy height models and field data[END_REF] also take into account the match between the ALS digital terrain model and topographic measurements. Some limitations are that these methods rely on single tree detection, which might be difficult in the case of low ALS density, or on height field measures, which are generally measured only for a sample of trees. An alternative would be to use the values of the diameter, which is the most common dendrometric measure and is highly correlated with tree height.

The main objective of this article is to investigate whether the diameter information suffices for co-registration and to evaluate the sensitivity to the number of georeferenced trees. Results obtained with diameter will also be compared to those obtained with height, either with cross-correlation or by adapting the cost function based on the difference between field and CHM heights [START_REF] Dorigo | An application-oriented automated approach for co-registration of forest inventory and airborne laser scanning data[END_REF]. The secondary objective is to evaluate to what extent the co-registration affects the assessment of the accuracy of ALS prediction models. The whole workflow is tested on two datasets with different forest structures and ALS pulse densities.

Material

Study Areas

The first study area is the public forest of Prénovel and Les Piards (Jura, France, Figure 1). The elevation range is narrow (900 to 1060 m), but the topography is rough. The forest is mainly constituted of uneven-aged, mixed stands of silver fir (Abies alba), Norway spruce (Picea abies) and beech (Fagus sylvatica). Permanent sample plots were installed in 2005, and the inventory was updated between 9 September and 6 October 2011. The sampling design is a grid pattern of a size of 200 m × 200 m, which results in 139 plots inside the forest (583 ha). Plot centers were positioned with a recreational-grade GNSS receiver. On each plot, all trees with a diameter at breast height (DBH) above 7.5 cm were inventoried up to a radius of 10 m and all trees with a DBH above 27.5 cm up to 17 m. For each tree, the azimuth and distance from the plot center, the DBH and species were recorded. Forest statistics are displayed in Table 1. The second study area is located in the department of Haut-Rhin (France), in the Vosges mountains (Figure 1). The elevation range is 200 to 1424 m. Ninety five plots were inventoried between 12 December 2011, and 19 April 2012. The sampling was designed to distribute plots into several environmental categories, including slope, site quality and forest stand type. The distribution of plots in stand development categories was: 14 young stands, 35 active growth stands, 29 mature stands and 17 uneven-aged stands. Plots are distributed equally into the poor and rich site quality categories and in the broadleaved, mixed and coniferous categories. Plot center positions were acquired with a survey-grade GNSS receiver during at least 15 min and post-processed for differential correction. All trees higher than 1.3 m and within 6 m from the plot center were inventoried. Trees with a DBH above 17.5 cm were inventoried up to a 15-m horizontal distance. For each tree, the DBH, azimuth and distance from the plot center were measured. The heights of the six trees with the largest diameters were also recorded. The main encountered species were silver fir (30% of basal area), beech (23%) and Norway spruce (18%).

Airborne Laser Scanning Data

In the Jura study area, the flight took place in leaf-on conditions, from 16-17 September 2012, with a full waveform RIEGL LMS-Q560 scanner on a fixed-wing aircraft. The flight speed was 180 km•h -1 at 500 m above the ground, with a strip overlap of 60%. The pulse frequency was 180 kHz with a scan angle of ±30 degrees. The obtained mean pulse density on the 9.2 km 2 was 9.3 m -2 . Pre-processing of the raw files was done by the contractor. Echoes were extracted and georeferenced with the RIEGL software suite. The resulting point cloud was classified into two classes (ground and vegetation) with TerraScan. The pulse density on the field plots ranges from 6.3 to 14.1 m -2 .

In the Vosges study area, airborne laser scanning data were acquired by the French Institute of Geographical and Forest Information during March and April 2011. Vegetation was mainly in leaf-off condition, and snow remained on the ground in some high areas. The scanner was an Optech ALTM3100 mounted on a fixed-wing aircraft flying at 290 km•h -1 at 1500 m above the ground. The pulse frequency was 71 kHz with a scan angle of ±16 degrees. The obtained mean pulse density on the 1362 km 2 was 2.6 m -2 . TerraScan was used for echo classification in ground and non-ground categories. The pulse density on the field plots ranges from 1.1 to 4.1 m -2 .

Methods

Co-Registration Algorithm

The workflow for the co-registration of a plot is presented in Figure 2. The required inputs are: trees position and value (diameter or height), plot center position and radius and the canopy height model (CHM). The CHM is the difference between the digital surface model, calculated by retaining the altitude of the highest ALS point in each pixel, and the digital terrain model, calculated by bilinear interpolation of ALS ground points at the center of pixels. The spatial resolution is one meter for all rasters. Even with this coarse resolution, the low density in the Vosges study area results in blank pixels with no ALS points. Besides, the CHM obtained in areas with leaf-off acquisition might display several low pixels, due to a deeper penetration of laser pulses inside the crown of broadleaved trees. A 3 × 3 median filter is applied to the CHM to remove this noise while retaining crown shapes [START_REF] Popescu | Measuring individual tree crown diameter with lidar and assessing its influence on estimating forest volume and biomass[END_REF]. The resulting filtered image is noted c. The trees map is rasterized into an image p by retaining for each pixel the largest value of the trees whose centers are contained in this pixel. Pixels without trees have their values set to zero. A plot mask m is computed by retaining only pixels whose centers are located inside the plot circle. A shift (dx, dy) is applied to the position of the mask and tree images. The correlation (cor) between the filtered canopy height model c and the shifted tree image is calculated over the mask according to Equation [START_REF] Duplat | L'inventaire par échantillonnage[END_REF]. When height measures are available, a weighted mean absolute error (wmae) between the tree image p and the filtered canopy height model c is computed according to Equation [START_REF] Naesset | Predicting forest stand characteristics with airborne scanning laser using a practical two-stage procedure and field data[END_REF]. 2 , where x is the mean of x in m

cor dx,dy = i∈m (p i -p)(c i -c) i∈m (p i -p) 2 i∈m (c i -c)
(1)

wmae dx,dy = i∈m |p i -c i | × p 2 i i∈m p 2 i ( 2
)
Figure 2. The workflow for the co-registration of a forest plot with the correlation (cor) method. For the weighted mean absolute error (wmae) method, the workflow is similar, except that the minimum value is used to select the correction offset. It is expected that the largest trees on a plot are likely to be the dominant ones, so that the height of their apices will be correctly represented in the ALS CHM. Smaller trees are likely to be suppressed, and the CHM values at their positions are likely to be higher than their heights due to the crown of neighboring dominant trees. The absolute error is thus weighted by the squared value of tree heights. This weight replaces in the equation proposed by Dorigo et al. [START_REF] Dorigo | An application-oriented automated approach for co-registration of forest inventory and airborne laser scanning data[END_REF] the tree class parameter, which is an "indicator of the visibility of a tree in the CHM", but not recorded in the Jura and Vosges datasets. The tree class parameter turned out to have little impact on their co-registration results, but the field data were based on angle-count plots where small trees account for a smaller proportion than on a fixed-radius plots. The correlation cor and weighted height error wmae are computed for dx and dy with one meter increments and within a search radius of 20 m from the initial position. Their values are then rasterized (Figure 3). The hypothesis is that when the co-registration is correct, the correlation between the trees image and the filtered canopy height model is maximal and that the weighted height error is minimal. Therefore, the coordinates of the maximum (resp. minimal) value in the cor (resp. wmae) image are retained as the shift to be applied to the trees' position. In order to estimate the reliability of the co-registration, the following values are recorded from the cor image: the position and value of the maximum (max1), the median value of a 3 × 3 window centered on the maximum (med1) and the value of the second local maximum (max2). In the case of a wmae image, the values correspond to the minimum instead of the maximum. The workflow is applied to all plots of both study areas, with the GNSS measurement as the initial plot center position. Corrected positions obtained after applying the shift are manually checked. Based on this correction and on his own appreciation, an operator proposes a shift for a better co-registration. In the case that the operator does not find any unambiguous correction, the plot is flagged as "doubtful".

Influence of Co-Registration on the Accuracy of Prediction Models

In order to evaluate the effect of co-registration on the accuracy of ALS prediction models, the area-based method is applied to calibrate a prediction model for basal area. A Box-Cox transformation is first applied to the dependent variable (basal area) in order to normalize its distribution. Independent variables are selected among a list of ALS metrics by retaining the group that yields the highest adjusted-R 2 in ordinary least squares regression when testing all possible combinations with at most three variables. Candidate variables include 16 height metrics: the 20th, 40th, 60th, 80th and 99th height percentiles of points above two-meter height for the "single", "first of many" and "last" point groups; and the mean height of all points above two meters h mean . Percentiles are abbreviated as h g,x with g ∈ {s, f, l} the point group and x the percentile value. Three density metrics are considered: the proportion of points above two meters that are located below h 0.99 2 , 2h 0.99 3 and 5h 0.99 6 with h 0.99 the 99 th percentile of heights above two meters in each plot. They are abbreviated d x with x the ratio applied to h 0.99 . Two CHM metrics are added: the mean value of the CHM on the plot CHM mean and the percentage of CHM values above two meters on the plot CHM >2 .

Validation is based on the repetition of 1000 ten-fold cross-validations. For each repetition, ten groups of equal size are randomly constituted among the dataset. One group is discarded at a time when computing the coefficients of the regression for the previously selected dependent variables, and then, the prediction error is calculated for each plot of the discarded group. The global root mean square error (RMSE) for the repetition is recorded. This procedure is applied to the nine scenarios obtained by combining the three possibilities for choosing the positions of training data and validation data (GNSS, the result of cor method with all diameters and the operator), for both study areas. For "doubtful" plots where no operator position could be determined, the GNSS position is used.

Influence of Forest, Topography and ALS Data Parameters on Co-Registration

For each plot, the following forest parameters are computed: basal area, mean diameter, stem density, proportion of coniferous trees weighted by basal area and dominant height (mean height of the six largest trees, Vosges study area only). Two topographic parameters are recorded: slope, which is calculated by fitting a plane to the ALS ground points located inside the plot extent, and the altitude of the plot center. The following ALS parameters are computed: the pulse density, ground points density, vegetation points density and density of points below a 0.5-m height. Two indicators of the quality of the canopy height model are computed. Q b is the proportion of pixels inside the plot limits with no ALS points. It reflects the amount of smoothing required by the low pulse density. Q f is computed as the the mean absolute difference between the CHM and the canopy height model filtered by the median 3 × 3 window c, for the pixels located inside the plot limits, see Equation [START_REF] Hollaus | Operational wide-area stem volume estimation based on airborne laser scanning and national forest inventory data[END_REF]. It reflects the amount of pit-filling operated by the median filter, due to both blank pixels and low pixels inside the tree crowns.

Q f = 1 card(P ) i∈P |CHM i -c i | ,
with card(P ) the number of pixels inside the plot P (3)

Influence of the Number of Georeferenced Trees on Co-Registration

To assess the influence of the number of georeferenced trees on the possibility of automated co-registration, the workflow is applied to the plots with the GNSS positions and only the n t largest trees on the plot.

The tested numbers of reference trees are n t ∈ {2, [START_REF] Hollaus | Operational wide-area stem volume estimation based on airborne laser scanning and national forest inventory data[END_REF][START_REF] Baltsavias | Airborne laser scanning: Basic relations and formulas[END_REF][START_REF] Naesset | Assessing point accuracy of DGPS under forest canopy before data acquisition, in the field and after postprocessing[END_REF][START_REF] Andersen | An accuracy assessment of positions obtained using survey-and recreational-grade global positioning system receivers across a range of forest conditions within the Tanana valley of interior Alaska[END_REF][START_REF] Korpela | Mapping forest plots: An efficient method combining photogrammetry and field triangulation[END_REF][START_REF] Gobakken | Assessing effects of positioning errors and sample plot size on biophysical stand properties derived from airborne laser scanner data[END_REF][START_REF] Frazer | Simulated impact of sample plot size and co-registration error on the accuracy and uncertainty of LiDAR-derived estimates of forest stand biomass[END_REF][START_REF] Hauglin | Geo-referencing forest field plots by co-registration of terrestrial and airborne laser scanning data[END_REF][START_REF] Olofsson | A method for linking field-surveyed and aerial-detected single trees using cross correlation of position images and the optimization of weighted tree list graphs[END_REF][START_REF] Hyyppä | A segmentation-based method to retrieve stem volume estimates from 3-D tree height models produced by laser scanners[END_REF][START_REF] Dorigo | An application-oriented automated approach for co-registration of forest inventory and airborne laser scanning data[END_REF][START_REF] Pascual | Algorithm for improving the co-registration of LiDAR-derived digital canopy height models and field data[END_REF][START_REF] Popescu | Measuring individual tree crown diameter with lidar and assessing its influence on estimating forest volume and biomass[END_REF][START_REF] Hollaus | Accuracy of large-scale canopy heights derived from LiDAR data under operational constraints in a complex alpine environment[END_REF][START_REF] Vauhkonen | Comparative testing of single-tree detection algorithms under different types of forest[END_REF][START_REF] Monnet | Evaluation of a semi-automated approach for the co-registration of forest inventory plots and airborne laser scanning data[END_REF]19,20,22,24,26,28,30, all} for the cor method applied to diameters and n t ∈ {3, 4, 5, 6} for the cor and wmae methods applied to the height data of the Vosges study area. In the case that the number of trees in the plot is smaller than n t , all trees are used to compute the tree image.

Results

Co-Registration

Operator Validation

In the Jura study area, only one "doubtful" plot could not be unambiguously located by the operator (Plot 88). Seven plots were located outside of the 20-m buffer area investigated by the algorithm. In the Vosges area, ten plots were flagged as "doubtful". Two-tailed Wilcoxon rank-sum tests show that the proportion of coniferous trees weighted by basal area (p = 0.012), the altitude (p = 0.003) and, to a lesser extent, the amount of pit-filling of the CHM Q f (p = 0.078) have significantly different distributions for the "doubtful" and "non-doubtful" plots in the Vosges area. "Doubtful" plots have a lower coniferous proportion, a lower altitude and a noisier canopy height model.

For the "non-doubtful" plots, hereafter referred to as validated plots, the distance between the GNSS position and the operator position was 1.8 ± 2.1 m for the 85 remaining Vosges plots and 9.0 ± 8.7 m for the 138 Jura plots. In both study areas, this distance is not correlated with any of the forest or topographic variables (Spearman's ρ, p > 0.05 in two-tailed tests).

Algorithm Comparison

When the cor method is used with all diameters, the proportion of validated plots located no further than two meters away from the operator position is 91.3% for the Jura dataset and 65.8% for the Vosges dataset. These proportions increase respectively to 92.0% and 71.7% when the distance is below five meters. When the cor method is used with all available heights (Vosges study area), the proportions are smaller, respectively, 63.5% and 69.4% for two-and five-meter distances, but only six trees were sampled. With six diameters, the correlation results are lower, with, respectively, 61.2% and 68.2% for two-and five-meter distances. When the wmae method is used with six heights, the proportions decrease to 61.2% and 65.9%. Those results for the two-and five-meter distance thresholds are similar to the general tendency of the cumulative distribution (Figure 4). In the Vosges study area, the best co-registration is obtained with all diameters. With six values, using heights instead of diameters yields better results. The wmae method gives lower accuracy than the cor method.

Comparison of the GNSS, Algorithm and Operator Positions

Figure 4 also displays the cumulative proportion of distance between the GNSS and operator positions. Summary statistics are presented in Table 2. In the Jura area, the algorithm allows a better co-registration, as the GNSS curve is located below the curve of algorithm results. The mean distance decreases from 9.0 to 3.3 m. Meanwhile, the standard deviation slightly increases from 8.7 to 11.0 m. The improvement is also visible in the quartiles, as 75% of GNSS positions are more than 4.2 m away from the operator positions, whereas 75% of algorithm propositions lie within 1.1 m from the operator positions. In the Vosges area, the GNSS curve displays a plateau at 0.41 from zero to a 0.86-m distance. From a 1-to 2-m distance, the cumulative proportion is higher with the algorithm positions, but for larger distances, the GNSS proportion is higher. Indeed, the mean and standard deviation of errors are higher with the algorithm propositions than with the GNSS positions. The quartiles show that for more than half of the plots, the error is reduced, as the median error is 1.5 m for GNSS and 1.0 m for the cor method. However, for the other plots, the co-registration potentially leads to large errors. 3 displays Spearman's correlation coefficient ρ between the co-registration error, considered to be the distance between the algorithm and operator positions of the validated plots, and the forest, topographic and ALS variables. For the Jura study area, the only significant correlations are with the ratios of the maximum value in the correlation image to the second local maximum value max 1 /max 2 and to the median of its 3 × 3 neighborhood max 1 /med 1 . For the Vosges study area, the co-registration error is negatively correlated with basal area, mean diameter and coniferous proportion. When the wmae method is used, the error is positively correlated with stem density and negatively with the ALS vegetation point density. In all cases, the co-registration error is correlated with the ratio of the maximum value to the value of the second maximum in the cor image (in the wmae image, it corresponds to the first and second minima). Figure 5 displays the influence of the number of sample trees on the proportion of validated plots that are correctly co-registered. In the Jura study area, the cumulative proportion increases quickly from zero to a two-meter distance and then reaches a plateau. It is noteworthy that the proportion reached at the plateau is slightly lower when all trees are used than with a number between five and 26. The proportion of co-registration within two meters is between 91.3% and 93.5%, when more than five diameters are used. With three diameters, it is only 80.4%. The cumulative proportion of co-registered plots within rapidly as distance increases (not shown in the figure). Indeed, in the case of a failing algorithm, which randomly chooses a position inside the investigated area, it is expected that the probability to co-register a plot at a distance d is proportional to the prospected area at this distance, which is proportional to d 2 .

With height data in the Vosges area, trends for the cor method are similar to what is observed with diameter, but with slightly better results. The proportions of co-registered plots are 41.2% and 63.5% below two meters, 50.6% and 69.4% below four meters for three and six sampled trees, respectively. Trends for the wmae method are also similar, but with slightly poorer results. The proportions of co-registered plots are 37.6% and 61.2% below two meters, 48.2% and 65.9% below four meters for three and six sampled trees, respectively.

Prediction Models

With the Jura dataset, the selected independent variables are h s,0.2 , h f,0.6 and CHM >2 when the GNSS positions are used for the training data, h f,0.8 , h mean and CHM >2 with the algorithm positions and h f,0.4 , d 0.5 and CHM mean with the operator positions. With the Vosges dataset, the selected variables are h mean , d 0.83 and CHM >2 (GNSS), h l,0.8 , d 0.83 and CHM mean (algorithm) and h mean , d 0.83 and CHM >2 (operator). The mean values of the root mean square error (RMSE) of basal area prediction models obtained in 1000 ten-fold cross-validations are presented in Table 4. With the Jura dataset, the estimation of RMSE is mainly affected by the positions of the validation data. It is around 7.2 m 2 •ha -1 with the GNSS data, 6.4 with the algorithm data and 6.0 with the operator data. With the Vosges dataset, the value of the RMSE is lower when the GNSS or operator positions are used (around 9.45) than with the algorithm positions (around 9.59), but the tendency is reversed when the algorithm positions are used for the modeling. Table 4. Root mean square error of basal area prediction models (1000 repetitions of ten-fold cross-validations) depending on the plot positions used in the calibration and validation steps, for the two study areas. Units are in m

Discussion

Limitations of Reference Data

With a posteriori validation of the co-registration, a bias might exist as the operator would be tempted to consider a proposition as correct, even tough he would a priori have preferred a slightly different position. This explains the high proportion of operator positions located exactly at the GNSS positions and the low number of shifts in ]0, 1[, as the operator would rather keep the GNSS position rather than move it from a distance smaller than one pixel. The absolute precision of the manual co-registration is thus expected to be around one meter [START_REF] Dorigo | An application-oriented automated approach for co-registration of forest inventory and airborne laser scanning data[END_REF], while the relative position precision of a tree within a plot is around 0.5 m [START_REF] Hollaus | Accuracy of large-scale canopy heights derived from LiDAR data under operational constraints in a complex alpine environment[END_REF]. However, this value is based on the hypothesis that the positions of the maxima of the canopy height model are a correct indicator of the positions of stems. Broadleaved trees sometimes have multiple apices, and the highest one is not systematically located above the trunk. Coniferous trees, despite a stronger apical dominance, have not always a vertical stem. If those irregularities in tree growth result in random shifts of the apices relatively to the trunks for all trees in the plot, using a sufficient number of trees should result in an unbiased estimation of the plot center position. Meanwhile, in slope areas, normalizing the digital surface model with the digital terrain model leads to a deformation of the canopy shape and a possible virtual shift of the apices of round crowns (Figure 6), such as those of broadleaved trees. Moreover, in slope areas, trees tend to bend downslope or to develop a flag-shaped crown, so that the planimetric positions of apices are located downslope (Figure 6). In such cases, the co-registration is biased, because the crowns in the canopy height model are really or virtually shifted downslope compared to tree stems. The position error then depends on the slope, crown shape and tilt of trees. Methods that rely on the comparison of the positions of local maxima [START_REF] Olofsson | A method for linking field-surveyed and aerial-detected single trees using cross correlation of position images and the optimization of weighted tree list graphs[END_REF] avoid the deformation of tree crowns if local maxima are detected in the digital surface model, but fail to handle the case of tilted trees. Indeed, the co-registration of the CHM with the tree positions results in a co-registration of tree positions with the positions of their detected crowns. Anyway, if the final objective is to improve the analysis of remote sensing data with field plots, in our case, the modeling of stand parameters with point cloud data, then it is probably interesting to co-register the tree crowns rather than stems, as crowns are the most visible part of the trees in the remote sensing data.

GNSS Accuracy in Forests

Despite this potential error in the operator-validated ground truth, results show that the GNSS measurements made in the Jura area with a recreational-grade receiver are inaccurate. Seven plots are located further than 25 m, with a maximum of 73 m, so that it might be impossible to find those plots in the next field campaign. The mean error (9.0 m) is similar to the error of 98 manually checked plots of the Austrian National Forest Inventory (8.5 m, [START_REF] Dorigo | An application-oriented automated approach for co-registration of forest inventory and airborne laser scanning data[END_REF]). In contrast, the results obtained in the Vosges area with a survey-grade GNSS receiver and a long acquisition time are better. For an acquisition time of 20 min, Naesset and Jonmeister [START_REF] Naesset | Assessing point accuracy of DGPS under forest canopy before data acquisition, in the field and after postprocessing[END_REF] reported positioning errors of 0.49 ± 0.42, 0.75 ± 0.38 and 2.5 ± 2.31 m for plots whose basal areas are below 30, between 30 and 45 and above 45 m 2 •ha -1 , respectively. In the Vosges area, the values are 2.4 ± 2.5, 1.3 ± 2.0 and 2.2 ± 1.9. Surprisingly, higher errors are recorded for small basal areas, which might be due to inaccurate co-registration by the operator for stands with small trees or with a low number of trees. Errors for large basal area values are similar, but it is noteworthy that the "doubtful" plots, with potentially large errors, are not taken into account. The Scandinavian study also reported a correlation of the positioning error with basal area, whereas in both areas of the present study, no correlation was significant with any of the forest or topographic variables. This might be due to the fact that "doubtful" plots are representatives of the extreme tendencies, but missing from the analysis, and that for small distances, the error is not properly evaluated, due to approximations by the operator.

Possibility of Co-Registration

The co-registration by the operator is generally more difficult when the contrast in the CHM is low. This is particularly the case when the dominant trees are broadleaved. The crown of broadleaved trees is also less well described by ALS data in leaf-off acquisitions. With closed canopies, the absence of small gaps also makes the co-registration more difficult. Due to the spatial resolution of the CHM, the apices of small trees are also difficult to identify. In contrary, large trees, particularly broadleaved, may have a flat crown where the apex position is doubtful, especially when the CHM is distorted because of slope correction. Some plots with few trees also turned out to be difficult to co-register properly. Those observations confirm the results of previous studies where dense plots, and, more particularly, deciduous plots, had lower co-registration accuracies [START_REF] Olofsson | A method for linking field-surveyed and aerial-detected single trees using cross correlation of position images and the optimization of weighted tree list graphs[END_REF][START_REF] Dorigo | An application-oriented automated approach for co-registration of forest inventory and airborne laser scanning data[END_REF]. In the end, it is important to have enough trees for a better co-registration, and it is better if those trees are separated and of different sizes, which allows the operator to take advantage of the diameter information when comparing with the CHM heights. It is also crucial that the tree crowns are well represented in the CHM. If the ALS point density is too low, it is not possible to identify crowns, especially for small trees. Those empirical findings by the operator are partially confirmed by the differences between the "doubtful" and validated plots of the Vosges study area. "Doubtful" plots have significantly lower coniferous proportion and noisier CHM. They are also located lower in altitude. This link with a topographic variable might be explained by the fact that the coniferous proportion decreases with altitude (Spearman ρ = 0.39).

In the Jura study area, most trees are coniferous, and the forest stands are uneven-aged, so that on each plot, small and large trees are present, and the canopy is not closed. Moreover, the ALS point density is high, which makes it an easy case for co-registration. The main problems with this dataset are the presence of felled trees, as the ground data is one year older than the ALS data, and the very large positioning error, which requires exploring a large window.

The proposed algorithm does not rely on individual tree detection, which is implemented for pulse densities around 5-10 pulses•m -2 [START_REF] Vauhkonen | Comparative testing of single-tree detection algorithms under different types of forest[END_REF]. For such densities, results obtained on the Jura dataset show that co-registration performs properly. When the pulse density is lower (2-4 pulses•m -2 ), but the variations in the canopy height model are sufficient to produce a maximum in the correlation image, correct results can also be obtained, as shown with the Vosges dataset. With pulse densities below one pulse•m -2 , which can still be used for the implementation of the ALS area-based estimation method, one could expect that the modeling of the canopy is not sufficient to reflect the height variation and that the proposed approach will fail. In the absence of a field reference for the correct positions, the operator validation relied on the possibility to visually match the field trees with the CHM, so that, in the case of the Vosges dataset, some plots with a very low density had to be removed from the analysis.

The choice of a coarse resolution (one meter) for the rasters is driven by the objective to limit the amount of smoothing required to remove the noise of the CHM in the case of low pulse densities. Besides, the co-registration approach relies on correlation, which is mainly influenced by the largest variations linked to large trees and canopy gaps. Those features appear in the CHM, even at this resolution. In the case of stands with only small trees, co-registration might be improved when using finer spatial resolution if the ALS density allows it. In the Vosges area, the pulse density is a limiting factor. In the Jura area, there are no plots with only small trees. Preliminary results obtained with 0.5-m resolution yielded similar results [START_REF] Monnet | Evaluation of a semi-automated approach for the co-registration of forest inventory plots and airborne laser scanning data[END_REF]. The absence of improvement with a finer resolution might be due to the accuracy of the field tree positions.

Automated Co-Registration

Dorigo et al. [START_REF] Dorigo | An application-oriented automated approach for co-registration of forest inventory and airborne laser scanning data[END_REF] could automatically co-register 68% of 98 plots from the Austrian National Forest Inventory in Vorarlberg within five meters of the position previously determined by an operator. Those plots consisted in a combination of fixed-area (radius = 2.6 m) and angle count sampling. The proportion of correctly co-registered plots with the Jura data is higher (92.0%), which could be linked to several factors. First the plots are of a 17-m fixed-radius, so that more trees are inventoried inside smaller surfaces, leading to a better local fit between the stand described by the inventory and the canopy height model. As stated above, the forest stands are also particularly propitious to easy co-registration, despite a one-year gap between the ground and airborne acquisitions. In Vorarlberg, the time gap is up to 2-3 years. With the Vosges dataset, the obtained co-registration proportion is similar (72%), but the conditions are different: the time gap is smaller, plots are of a 15-m-fixed radius and the proportion of deciduous plots is higher.

With aerial images and azimuth-distances measures, Korpela et al. [START_REF] Korpela | Mapping forest plots: An efficient method combining photogrammetry and field triangulation[END_REF] were able to position reference points with decimeter accuracy in a stand dominated by Scots pine (Pinus sylvestris), Norway spruce and birch (Betula sp.). Their method can be applied to improve the co-registration of forest inventory plots when the azimuth and distance measures are available for a sample of trees. The obtained accuracy is better than with the methods based on ALS data, yet it might decrease in the case of stands dominated by large broadleaved trees, due to the uncertainty in the manual observations of treetops from multiple aerial images. Besides, results obtained with the fully-automated ALS-based methods can be further improved by a fast manual validation and correction of the algorithm propositions, especially if a procedure is designed to flag the uncertain cases [START_REF] Dorigo | An application-oriented automated approach for co-registration of forest inventory and airborne laser scanning data[END_REF]. In a boreal forest dominated by Scots pine and Norway spruce, Hauglin et al. [START_REF] Hauglin | Geo-referencing forest field plots by co-registration of terrestrial and airborne laser scanning data[END_REF] co-registered 82% (resp. 51%) of 71 plots within one (resp. 0.5) meter of the reference location, using terrestrial and low density airborne laser scanning data (0.7 pulses•m -2 ). This performance should yet be confirmed in the case of broadleaved forests. With the higher ALS density in the Jura area, 82% (resp. 59%) of plots are co-registered below one (resp. 0.5) meter with the correlation approach. In the case of existing field plots, the terrestrial laser scanning method requires additional fieldwork to perform the scans and also relies on tree detection in both the aerial and terrestrial laser data. A major advantage is that no information about tree positions are required.

The comparison of co-registration methods in the Vosges area shows that the cross-correlation yields slightly better results with heights than with diameters when the same number of tree measures are used. This could be expected, as there is some variability in the height-diameter relationship, which introduces some additional noise in the correlation between the field and ALS measures of trees. However, it is better to use all diameters than just a sample of heights. In the case that all tree positions are not available, results for the Jura dataset show that only with the five largest trees, the co-registration is almost as efficient as with all trees. For the Vosges dataset, ten trees are required to achieve similar results.

The correlation method searches for the best match between the relative variations of the CHM and those of the tree measures. The use of concentric fixed-radius plots with different DBH limits is sometimes problematic. Consider a plot covered with trees of 10 cm DBH, plus a few large trees on its western part. No trees are inventoried outside of the inner plot on the eastern part. If a canopy gap whose shape matches the area with no inventoried trees on the plot is present in the investigated area, the algorithm will erroneously match the plot to this location, because the relative variation between the eastern and western part of the plot would be apparently respected. The correlation algorithm is particularly sensitive to the presence of large gaps or stand borders, in relation to the proportion of trees inventoried in the different concentric plots with different DBH limits. This might explain why co-registration results are better with a lower number of trees. When small trees are added, they are sometimes all inventoried in the inner plot with a smaller DBH limit. Besides, small trees are likely to be suppressed trees, the values of which in the tree image will not correspond to variations in the CHM, where they are overtopped by the crown of neighboring, taller trees. To handle this issue, it is possible to model the tree crowns as Gaussian surfaces in the tree image [START_REF] Olofsson | A method for linking field-surveyed and aerial-detected single trees using cross correlation of position images and the optimization of weighted tree list graphs[END_REF], but this requires prior information about the relationship between tree size and crown shape. For a better performance, it seems important to find a trade-off between the number of considered trees and their repartition in the different concentric plots, but this optimal balance will probably depend on the forest structures (height distribution and spatial repartition of trees). In the Jura uneven-aged forests, large trees are present and spatially well distributed on all plots, so that good results are attained with only five trees and the best results with ten trees.

The weighted mean absolute error method is based on the absolute values of the CHM at the trees positions. It might also be affected by the different field sampling intensity in the different radii, even though it is not as flagrant as for the cor method. As it relies on the comparison of absolute values of height, it is probably affected by errors in the field and CHM heights, which are enhanced with broadleaved trees and in slope areas. Moreover, a shift in the positions of the tree apices relatively to the tree stems is also an issue. Unfortunately, due to the low number of sampled trees, it is not possible to investigate to what extent additional height measures could improve the co-registration.

Factors Affecting the Algorithm Co-Registration Error

In absence of high precision field data, co-registration error corresponds here to the distance to the operator proposition, which is based on visual analysis of the CHM and trees positions. In the Jura study area, the cumulative distribution of errors quickly rises and reaches a plateau at a two-meter distances, which shows that the distance is small for the correctly co-registered plots. No forest, topographic or ALS variable is correlated with the co-registration error. This could be linked with the low variability of forest stands and topography, and the high quality of the ALS data on the whole forest. The ratio of the maximum of the correlation image to the median of its neighborhood max 1 /med 1 indicates if the maximum is peak-shaped or round. Peak-shaped maxima correspond to situations where the shift associated with the global maximum is precise. In the Jura area, where the conditions are easy for co-registration, the question of error, hence, comes down to whether the proposed position is accurate or not.

In the Vosges area, the correlation method is affected by basal area and coniferous proportion. As mentioned before, broadleaved trees have less apical dominance and a more rounded crown, which makes the detection of the tree tops more difficult. The negative correlation with basal area is probably driven by the fact that for young stands or stands with a small number of trees, there is not enough information for a precise co-registration. The negative correlation with the mean diameter suggests that for big trees, there might be a large shift between the apices and trunks positions. In the Vosges area, both methods have a co-registration error linked with the ratio of the global maximum value to the second maximum value max 1 /max 2 . This ratio can be considered as a rough indicator of the signal-to-noise ratio. When its value is close to one, there are high chances that the global maxima is only an artifact caused by the irregular shapes of trees or to a low-quality CHM.

One would expect the ALS density to be a factor affecting co-registration; however, this is not significant in any dataset. In the Vosges area, the effect of ALS density might be hidden by the variability of forests characteristics, such as basal area, coniferous proportion and, possibly, the horizontal and vertical structures, which also influence the possibility and accuracy of co-registration. In the Jura area, the pulse density is high, so that it is probably not a limiting factor in this case.

Resulting Co-Registration and Prediction Models

In the Jura area, the co-registration results in better positions compared to the initial GNSS measurements, except for a few plots that are mainly located outside of the investigated radius. In the Vosges area, there is a small improvement for approximately half of the plots, but increased error for the remaining plots.

With the Jura data, the best model is obtained with the operator positions, with an RMSE of 5.88 m 2 •ha -1 . The models obtained with the GNSS or algorithm positions and validated with the operator positions have only a slightly lower accuracy (RMSE of respectively 6.09 and 6.03). However, their "apparent" error, which is obtained when calibration and validation steps are processed with the same positions, are respectively 7.13 and 6.44. Considering that the position error results in a random noise affecting the independent variables, the result of an ordinary least squares regression will be scarcely affected by the changes in ALS variables, provided that the number of plots is sufficient and that there are no outliers. Indeed, the noise will hardly modify the coefficients of the best fit, but will increase the This might explain why the calibration step is not affected by position error, whereas it is important to have precise positions for the evaluation of the prediction accuracy.

With the Vosges data, the unsupervised co-registration leads to large errors on some plots, so that the prediction models actually have a lower accuracy than with the GNSS positions. With the low level of positioning error (mean of 1.8 m), the correction of positions does not result in significantly better prediction models, which is consistent with the results obtained at the stand level [START_REF] Gobakken | Assessing effects of positioning errors and sample plot size on biophysical stand properties derived from airborne laser scanner data[END_REF][START_REF] Frazer | Simulated impact of sample plot size and co-registration error on the accuracy and uncertainty of LiDAR-derived estimates of forest stand biomass[END_REF].

Conclusions

This study demonstrated that the cross-correlation of diameter measures is a suitable method to co-register forest plots with ALS data. Moreover, co-registration results are almost as good as with height information and are obtained with only the largest trees. The cross-correlation method displayed some limitations when used with concentric plots with different diameter inventory limits. One solution might be to automatically adapt the number of trees considered in the analysis depending on the forest structure. Moreover, one could take advantage of the sampled heights in order to compare the results of the correlation and height error methods in order to provide the operator with concordant or complementary possibilities for a faster resolution of ambiguous situations.

A comparison of the results obtained on two different datasets showed that the accuracy of co-registration and its utility for the improvement of ALS prediction models highly depends on the GNSS positioning accuracy, the ALS data quality and the forest structure. With the dataset where the co-registration performed best and where the initial GNSS precision was low, using corrected positions produces slightly better models, but mainly allows a better estimation of their accuracy. With the dataset with already accurate GNSS positions, low-density ALS data and complex, mixed forests, the unsupervised algorithm is not able to improve the co-registration. Furthermore, when the positions are improved by an operator, there are no differences in the calibration and the validation of the prediction models.

Figure 1 .

 1 Figure 1. Location of the study areas and the positions of forest inventory plots.

Figure 3 .

 3 Figure 3. The workflow applied to Plot 8 of the Jura study area. (a) The plot position according to the GNSS positioning (• tree positions with a symbol size proportional to the tree diameter, × plot center); the background is the canopy height model; (b) The correlation image (× global maximum, its 3 × 3 neighborhood and + the second maximum); (c) The co-registered plot according to the (-10, -2) shift corresponding to the global maximum in the correlation image.

Figure 4 .

 4 Figure 4. Cumulative proportion of validated plots as a function of the distance in meters between the operator position and the Global Navigation Satellite System (GNSS) position or the algorithm position depending on the method used (cor correlation or wmae weighted mean absolute error), the input tree variable (diameter or height) and the number of considered trees (six or all). The curves are cropped to the distance interval 0-5 m.

Figure 6 .

 6 Figure 6. The difference between stem position and apex position in tilted trees (a); The shift of the apex position of a round crown when correcting a slope (b, c).

Table 1 .

 1 Forest statistics for the two study areas. The means ± standard deviations are on the first line, minimum-maximum values are on the second one.

	Study	Dominant Basal area Coniferous proportion	Stem	Mean
	area	height (m) (m

  2

  •ha -1 ) weighted by basal area density (ha -1 ) diameter (cm)

	Jura	31.8 ± 9.2	0 .87 ± 0.12	510 ± 150	25.8 ± 7.0
	n p = 139	13.3 -60.5	0 .49 -1.0	110 -1, 520	14.8 -54.5
	Vosges	23.9 ± 7.0 40.7 ± 9.2	0 .53 ± 0.38	1, 990 ± 2, 900	20.8 ± 11.5
	n p = 95 7.6 -38.4 10.8 -75.2	0-1	170 -17, 700	2.4 -51.5

Table 2 .

 2 Statistics on the distance in meters between the GNSS and operator positions and between the operator positions and various algorithm propositions for the validated plots. On the first line is the mean ± standard deviation; on the second line are the quartiles.

		GNSS	cor , Diameter, All cor, Diameter, 6 cor, Height, All wmae, Height, All
	Jura	9.0 ± 8.7 4.2, 6.6, 10.7	3 .3 ± 11.0 0.5, 0.5, 0.71	3 .1 ± 10.1 0.5, 0.5, 1.1
	Vosges	1.8 ± 2.1 0, 1.5, 2.8	5 .4 ± 7.5	5 .3 ± 7.0	5 .4 ± 7.3	6 .1 ± 7.5

0.71, 1.0, 11.3 0 .9, 1.1, 9.8 0.8, 1.1, 10.8 0 .9, 1.4, 12.8 4.1.4. Influence of Environmental Variables on Co-Registration Error Table

Table 3 .

 3 Spearman correlation coefficient ρ between forest, topographic, airborne laser scanning (ALS) and co-registration variables and co-registration error (the distance between the operator and algorithm positions for the validated plots) for different algorithm settings. Significance levels: * p < 0.05, ** p < 0.01, *** p < 0.001. CHM, canopy height model.

	Variable	Vosges cor, Diameter, All cor, Height, 6 wmae, Height, 6	Jura cor, Diameter, All
	Dominant height	-0.25 *	-0.11	-0.28 *
	Basal area	-0.43 * * *	-0.40 * * *	-0.29 * *	-0.05
	Stem density	0.2	0 .09	0.28 * *	0.03
	Mean diameter	-0.31 * *	-0.22 *	-0.31 * *	-0.04
	Coniferous proportion	-0.54 * * *	-0.53 * * *	-0.31 * *	-0.15
	Slope	-0.16	-0.08	-0.05	0.14
	Altitude	-0.12	-0.36 * * *	-0.37 * * *	0.04
	Pulse density	0.08	0.11	-0.06	-0.05
	Ground point density	0.18	0.17	-0.02	0.11
	Vegetation point density	-0.08	-0.21	-0.29 * *	-0.02
	Points below 0.5 m density	0.19	0.18	-0.02	0.05
	CHM pit-filling	-0.03	0.13	0.06	-0.11
	CHM empty pixels	0.06	0.05	0.08	-0.04
	max 1 /max 2	-0.41 * * *	-0.42 * * *	0.41 * * *	-0.20 *
	max 1 /med 1	0.04	0.02	-0.04	0.24 * *
	4.1.5. Influence of the Number of Sample Trees

  2

  •ha -1 .

	Calibration	Validation (Jura) GNSS cor, Diameter, All Operator	Validation (Vosges) GNSS cor, Diameter, All Operator
	GNSS	7.13	6.49	6.09	9.46	9.58	9.42
	cor, Diameter, All	7.17	6.44	6.03	9.72	9.31	9.62
	Operator	7.34	6.37	5.88	9.47	9.60	9.43

Acknowledgments

This work was funded by the European Commission (project Alpine Space 2-3-2-FR NEWFOR) and by the French National Research Agency (project ANR-2010-BIOE-008 FORESEE). We thank the French Forest Office (Office National des Forêts) and project BGF (Biodiversity and Forest Management) for the Prénovel-Les Piards field data.

one meter is also higher when only ten trees are used, compared to the proportion obtained with higher numbers of trees.

Figure 5. Cumulative proportion of validated plots as a function of the distance between the operator position and the algorithm position. Line colors refer to the sample size. Subfigure titles refer to the study area, the method used (cor correlation or wmae weighted mean absolute error) and the input tree variable (diameter or height). With the Vosges data, similar trends can be observed, except that the proportion increases more slowly and the plateau is reached around four meters. The proportion of co-registration within two meters is between 64.7% and 71.7%, when more than five diameters are used. With three (resp. five) diameters, it is only 40.0% (resp. 57.6). Whereas in the Jura area, the proportions reached at the plateau are similar when more than five trees are used, they are quite different in the Vosges area, with 64.7% for five trees, 71.7% for ten or all trees, 76.5% for 20 or 26 trees and 81.2% for 15 trees at a four-meter distance. With the Jura data, using ten trees yields a higher proportion of co-registration below one meter, but for larger distances, the proportions are similar whatever the sample size. With the Vosges data, using fifteen trees or more yields a higher proportion of co-registration, but curves then tend to diverge. In both study areas, the cumulative proportion of plots co-registered above five meters then increases slowly and more

Author Contributions

Jean-Matthieu Monnet: co-registration workflow, statistical analysis, writing. Éric Mermin: field data preparation, validation of co-registration results.

Conflicts of Interest

The authors declare no conflict of interest.