1 Introduction↩︎

National Forest Inventories (NFIs) are designed to estimate and sometimes to monitor the change in forests characteristics at National scale with a given level of reliability. In order to downscale NFI estimates without increasing the cost of field surveys while preserving the reliability of the estimates, methods have been developed to combine NFI data with auxiliary data. Those so-called Multisource Forest Inventories (MSNFIs) usually rely on remote sensing data selected for their correlation with the forest attributes of interest, their low cost, and their availability at the national scale, together with their regular long-term renewal [1]. In terms of correlation, 3D data derived from airborne laser scanning, stereophotogrammetry, and, to a lesser extent, radar interferometry are the most suitable sources of auxiliary information for structural forest attributes such as volume and biomass [2]. With the development of national-scale LiDARdatabases, interest in their use within MSNFIs has been renewed. That said, data acquisition and processing are not standardized, making their integration into MSNFIs a challenge.

LiDARdata processing methods are established for point clouds acquired locally, and often with similar acquisition conditions. Nevertheless, a significant challenge arises with national acquisitions, which often rely on point clouds acquired under varying acquisition conditions [3]. Early studies have pointed out that LiDARacquisition parameters such as flying altitude, scan angle, pulse density, and sensor characteristics influence the estimation of forest structural attributes. [4] showed, from four Airborne Laser Scanning (ALS) acquisitions over 40 field plots in mature conifer forest, that instrument type, flying altitude, and pulse repetition frequency all affect ALS-derived height and density metrics, with systematic differences between acquisitions of up to 2.5% for mean tree height and 10.7% for timber volume, despite only minor differences in overall model precision. [5] further demonstrated, using linear regression on \(\sim\)​1,800 field plots across 76 heterogeneous ALS projects in the Southeastern US, that combining heterogeneous LiDARdatasets over large areas introduces additional variability, with factors such as point density, scan angle, and vegetation heterogeneity influencing model accuracy, though canopy height could still be predicted with a root mean squared error (RMSE) of 3.0\(\:\)m (R\(^2\) = 0.74), reduced to 2.4 m\(\:\)in structurally homogeneous stands. [6] specifically highlighted the scan angle as a critical parameter for foliage profile retrieval, identifying an optimal angle of 20\(^{\circ}\) regardless of flying altitude or pulse density, while showing that combining multiple scan angles can further improve estimation accuracy.

Differences between leaf-onand leaf-offacquisition conditions have been shown to substantially affect canopy structure representation and the estimation of forest attributes. In deciduous stands in particular, the absence of foliage leads to a greater penetration of the signal through the canopy, resulting in a higher number of single and last returns originating from the ground and a shift of the return distribution toward lower heights under leaf-offconditions [7]. These structural differences have direct consequences for model transferability: applying a k-nearest neighbors (kNN) algorithm, a non-parametric method that predicts values based on the most similar observations, to predict volume-related values by analyzing the volumes of the most similar data points, trained under leaf-onconditions to leaf-offdata (and vice versa) introduces substantial bias and should therefore be avoided [8]. This degradation in performance has been quantified in boreal forest contexts using Random Forest models across two inventory areas in Finland, where cross-using leaf-onand leaf-offmodels increased relative Root Mean Square Error (rRMSE) and led to both over and underestimation, particularly in deciduous-dominated plots. In these cases, error nearly doubled (from \(\sim\)​33% to \(\sim\)​56%) and bias reached up to 29%, although calibration via an empirical ratio estimator was shown to partially mitigate these errors [9]. Consistent with these findings, further evidence shows that cross-application of random forest models and acquisition conditions results in large increases in error for both forest types, with average absolute increases of rRMSE up to 16.8% for coniferous stands and 25.8% for deciduous stands, reinforcing that such mismatches should be avoided [10]. These results highlight the critical importance of accounting for phenological conditions when developing and applying LiDAR-based forest attribute models.

The reliability of estimates also vary with the forest type (even at a coarse level such as coniferous, deciduous, and mixed forests). For deciduous trees, leaf conditions influence the distribution of points, with skewed distributions toward the upper crown in leaf-onconditions with limited penetration to the ground, to shifted point distribution toward lower canopy levels with increased ground sampling in leaf-offconditions [11]. Pure coniferous stands tend to be minimally impacted by acquisition conditions [12], albeit the residual presence of broadleaf trees, lower vegetation and partial winter needle thinning [11] still introduces variability. One can conclude that most recent findings and models are heavily area- and acquisition-dependent [13].

Taken together, these acquisition-related and phenological sources of variability make it particularly challenging to develop forest variable estimation models that generalize beyond their local or regional scales, and the conditions under which they were trained. This motivates the development of more robust, large-scale approaches.
Estimating forest attributes at national scale with a single, generic model capable of handling a large diversity of forest types and LiDARacquisition settings remains, to the best of our knowledge, an open problem as to date, the authors are not aware of such investigations. In this paper, we address this challenge by leveraging the French national LiDARHD program. Such acquisition program provides a uniquely large and heterogeneous dataset of 3D point clouds, with minimum density of 10\(\:\)pts/m\(^2\) and acquisitions distributed across both leaf-offand leaf-onconditions. French forests are dominated at 68% by deciduous stands, a forest type whose structural signal is strongly modulated by phenology: such dataset captures a variability rarely encountered in previous studies. This heterogeneity is further compounded by the five-year acquisition window, which spans multiple campaigns conducted under different seasonal and instrumental conditions, making the dataset not only large but intrinsically variable.
Harmonization remains, to our knowledge, largely unexplored in this context, although several strategies have been proposed in the literature to address acquisition heterogeneity. One approach consists in developing separate estimation models for each acquisition campaign or season, as implemented by [14] within the Swedish NFI framework: block-specific (usually 25\(\:\)km by 50\(\:\)km in size) linear regression models were fitted using the nearest NFI plots to account for the variability introduced by different ALS sensors and acquisition settings across the country. Alternatively, linear mixed-effects models have been proposed to account for both forest type and acquisition variability, as demonstrated by [15] in the context of the Norwegian forest resource map, where data from 367 ALS projects acquired over ten years were combined with NFI plots through stratification by species and maturity class. Species-specific models yielded rRMSE of 35%, 34%, 31% and 12% for volume, aboveground biomass, basal area and Lorey’s height, respectively, outperforming general unstratified models by 2–7 percentage points.

However, explicitly stratifying models by acquisition condition and forest type remains challenging at national scale, where the compounding diversity of forest composition, climate, and topography, combined with the difficulty of aggregating small areas due to error propagation [16]–[18] limits the practical applicability of such approaches. The question of spatial transferability further illustrates these limitations: [19] proposed a Random Forest approach combining ALS and NFI data across Finland, yet reported a systematic increase in RMSE with geographic distance from the training set, highlighting the limited generalization capacity of such models under spatial heterogeneity. Moreover, phenology represents only one dimension of the problem: flight altitude, sensor type, and point density each influence point cloud properties in ways that interact with forest type and structure. This renders exhaustive stratification both impractical and potentially inadequate. A more robust solution requires a form of harmonization, whether applied at the input data level (making 3D point clouds similar in terms of geometry) or at the prediction level (estimating consistent forest variable whatever the above-mentioned fluctuations).

One straightforward approach would consist in homogenizing datasets by converting all 3D point clouds into a specific targeted condition (e.g., leaf-on). However, beyond phenology, many other acquisition characteristics like flight altitude, sensor type, density as well as their interaction with each forest type and structure have been shown to affect point cloud properties [20], [21]. Explicitly normalizing for every known effect is a complex reconstruction task, and likely impractical for large-scale operational use on heterogeneous datasets [22]. It is also based on the assumption that it is possible to retrieve the correct geometry of forests. However, the literature mainly focuses on point cloud completion and has only been addressed at the individual tree level or for opaque surfaces [23], [24].

Instead, another approach could focus on predicting similar outputs. It assumes a statistical model can capture such heterogeneity and can automatically learn to adjust its prediction to the acquisition and seasonal context. This hypothesis is supported by [10]. They showed, using Random Forest models trained on 787 ground plots in a lodgepole pine-dominated forest in Alberta, Canada, that pooling leaf-onand leaf-offdata into a single model results in only marginal performance loss (<2% rRMSE for coniferous, near-zero for deciduous stands). This suggests that a single model can implicitly accommodate both acquisition conditions without explicit stratification. At larger scales, [15] further demonstrated that such an approach remains feasible even under strong acquisition heterogeneity, successfully combining data from 367 ALS projects acquired over a decade across 17\(\:\)Mha of Norwegian forest within a single modeling framework.
Traditional machine learning approaches rely on pairing carefully designed statistical features derived from point clouds with ground-truth measurements from forest inventory plots to perform predictions. Regression models for forest attributes (volume, biomass, basal area, etc.), typically based on Random forest and kNN algorithms, often struggle to generalize when applied over heterogeneous conditions. Even if methods like kNN have the advantage of providing confidence intervals using empirical quantiles of the neighbors or bootstrapping [25], they are still limited by their ability to learn robust transferable representation [8], [26]. Recent advances in deep learning (DL) offer new perspectives to overcome these limitations. [27] highlight that DL regression of forest attributes can produce more accurate results compared to state-of the art machine learning methods. By operating directly on raw point clouds, DL models preserve the full three-dimensional structure of the data and automatically learn relevant spatial features, bypassing the need for manual feature selection, and fluctuation in the local geometry of point clouds.

Several DL architectures have been proposed to process 3D point clouds, broadly falling into two categories: point-based and voxel-based approaches. Point-based approaches rely on sampling or grouping strategies to capture local context [28], [29], a strategy that is not well suited to the structural irregularity of forest canopies, where vegetation density varies continuously and no clear object boundaries exist. Voxel-based discretizations, in contrast, preserve structural information while maintaining computational efficiency, and are better adapted to the sparse and heterogeneous nature of forest point clouds.
[30] evaluated two DL architectures for forest attribute regression in a temperate mixed forest: the point-based Dynamic Graph Convolutional Neural Network (DGCNN) and the voxel-based Octree Convolutional Neural Network (OCNN). Their results show that while both DL methods outperformed Random Forests, the voxel-based OCNN achieved the lowest prediction errors and reduced mean absolute percentage error. [31], in their study on airborne LiDARbiomass estimation over Danish national forests, adapted and compared three distinct DL architectures for regression: the point-based PointNet, the point-based Kernel Point Convolution (KPConv), and the voxel-based Minkowski convolutional neural network. Among these, the voxel-based Minkowski CNN consistently outperformed both point-based alternatives across all forest attributes and evaluation metrics. Both studies highlight that (i) voxel-based architectures exhibit greater robustness to variations in forest type, (ii) reliable predictions are conceivable by training a DL model with multi-annual NFI and LiDARdata that may not be fully temporally aligned.
In this paper, we propose a generic deep-learning solution, dubbed FLORA(Forest LiDAROctree Regression with Auxiliary data), for estimating forest attributes at a national-level, that can handle a large diversity of both forest and LiDARdata conditions usually encountered at national scales. FLORAis a voxel-based approach designed to estimate forest structural attributes from heterogeneous LiDARdata, validated over France where large-scale NFI and LiDARdata are jointly available. The proposed data does not have counterpart of this size so far [32]. Our model predicts six attributes commonly derived from LiDARdata [33]. We focused here on the following six : dominant height, total volume, deciduous stand volume, coniferous stand volume, stem density, and basal area. Model performance was evaluated using the rRMSE (%) and the coefficient of determination (\(R^{2}\)), which measure prediction error and explained variance, respectively. This work addresses three main objectives. We first investigate whether a single DL model trained on data acquired under contrasting phenological conditions (leaf-onand leaf-off) can generalize across seasons as compared to models trained with season-specific data , and whether this robustness holds across forest types. We then examine the specific contribution of overlapping plots, surveyed under both leaf-onand leaf-offconditions, as a training strategy that anchors the model to paired observations of the same forest under varying acquisition contexts. Finally, we assess the added value of auxiliary contextual variables in improving prediction accuracy across heterogeneous acquisition conditions.

2 Materials and methods↩︎

2.1 Data↩︎

Our approach perform at plot-level and is framed to handle a large diversity of context, has found in the French mainland territory. The advantages are two-fold: this ensures (i) a large amount of data and (ii) significant diversity, therefore capturing very heterogeneous behaviors necessary for our generalization task.

2.1.1 Study area↩︎

French forests cover approximately 17.5 million hectares, representing about 32% of the land area of mainland France [34]. They are predominantly composed of broadleaved species, with around 68% of the forest areas consisting of deciduous stands. Forest ownership is largely private, as nearly three quarters of forest land is privately owned. In terms of stand structure, more than half of French forests consist of mixed stands, combining either different broadleaved species or mixtures of broadleaved trees and conifers, while pure coniferous stands are less frequent. In metropolitan France, forests contain nearly 190 different tree species, distributed among seven main species: oak, beech, chestnut, maritime pine, Scots pine, spruce, and fir [35].

This structural and compositional diversity is closely linked to the high climatic and topographic heterogeneity of France, which includes oceanic, continental, Mediterranean and mountainous climates. As a result, French forests display a wide range of forest types and structures [36], [37].

2.1.2 LiDARdata↩︎

We benefit from the large scale 3D airborne LiDARdata acquired in the framework of the French National LiDARHD program, started in 2021 for a 5\(\:\)years duration. The objective of the program is to acquire and release in open access LiDARdata over the metropolitan area and ultramarine territories, with the exception of French Guyana, with a point density of at least 10\(\:\)pulses/m\(^2\). Such a pulse density was selected to fulfill spatial sampling requirements for both Digital Terrain Model (DTM) generation, land-cover mapping and statistical tasks [38], reaching sometimes 40\(\:\)pts/m\(^2\). To achieve such a large multi-year campaign, the territory was divided into 50\(\:\)km\(\times\)​50\(\:\)km patches, which were each attributed to one of the six selected private companies or to IGN. Acquisitions were conducted throughout the year, with survey parameters adapted to the relief and the LiDARdevice. Several LiDARsystems were employed, with varying emitted energy, pulse rate and operated under multiple flight heights, plane speed and scan angle ranges (Fig. 1). In addition, without inter- and intra-annual planning in the survey periods (e.g., based on forest management or landscape consistency), adjacent patches may have been acquired by different companies, under different seasons, and distinct years. Edge effects between patches are reduced by integrating an overlap region of 2-5\(\:\)km in which at least two surveys were performed.

Fig. 1 shows our dataset exhibits a significant heterogeneity in terms of season and LiDARsurveys, leading in very distinct samplings of the tree canopy and the vegetation underneath.

Raw data processing was mainly performed using automatic routines involving TerraScan software, and classified with Myria3D DL tool [39]. We used the classified point cloud as the input data for our approach, as it allowed us to directly retain only the ground and vegetation classes and discard the other points, thereby simplifying the preprocessing steps.

2.1.3 French national forest inventory data↩︎

The French NFI is based on a systematic sampling of semi-permanent plots distributed across the territory on a 1\(\:\)km grid. Each year, 1/10 of the grid is explored and around 14,000\(\:\)plots are surveyed on the field, half of them been visited for the first time, and the other half have been revisited 5\(\:\)years apart. Quantitative and qualitative data on standing forests and vegetation are collected in the field using concentric circular subplots of 6, 9, 15, and 25\(\:\)m radius. Stand description is done on the largest plots. Dendrometric data are collected on the three other radius, according to tree diameter: trees with a Diameter at Breast Height (DBH) in the range [7.5\(\:\)cm – 22.5\(\:\)cm] are measured on the 6\(\:\)m radius plot, those in the range [22.5\(\:\)cm \(-\) 47.5\(\:\)cm] are measured on the 9\(\:\)m radius plot, and larger trees on the 15\(\:\)m radius plot. Plot\(-\)level estimates of forest attributes are derived on a per \(\:\)ha basis using the probability of inclusion of trees, field measurements as well as allometric models for some variables like volume. In the framework of this research, we relied only on first visit plots and considered the stem density (stems ha\(^{-1}\)), the basal area (m\(^2\) ha\(^{-1}\)), the stem volume (m\(^3\) ha\(^{-1}\)) and its distribution between hardwoods and conifer, as well as dominant height (m).

2.1.4 Our dataset↩︎

The study focuses on the mainland territory including Corse. 80% of the territory was covered at the time of this research. We select NFI plots with less than a five-year difference between field measurements and the LiDARacquisition. This results in a total of 33,877 plots (see Table 1 and Fig. 2), which covers approximately 7.6 million m\(^2\), that sample most of the French sylvo-ecoregions (SER), French territory division based on forest ecosystems and management practices.The dataset includes 60 main species, with an average of 502.1 plots having a similar main species (\(\pm\) 901.8 SD). The three most represented species are pedunculate oak (4,166 plots), sessile oak (3,333 plots), and European beech (3,251 plots). We extract LiDARHD point clouds centered on NFI locations, and link them to the corresponding field measurements. LiDARpoint clouds were already classified and a 1\(\:\)m DTM made jointly available.

We used the LidR package [40] to generate Canopy Height Models (CHM) and to extract point clouds for each plot. To ensure data quality, we applied a three-step filtering procedure, retaining 32,052 plots (93.6% of the initial dataset).

• Plots with a CHM containing more than 10% of missing values were excluded: we discard plots only partially covered by the LiDARacquisition.

• Plots with large deviations occurring between the field campaign and the LiDARacquisition were removed. This is likely to reflect stand disturbances, e.g., harvesting or windthrow. We compared the canopy cover fraction derived from the CHM against the one recorded during field surveys as a qualitative estimate of the canopy cover by class of 10%. Plots exhibiting a discrepancy exceeding \(\pm 2\sigma\) of the observed distribution of differences were discarded. The canopy cover was defined in the LiDARdata as the proportion of pixels with height above 2\(\:\)m, a threshold widely used in the literature to separate low vegetation from actual canopy returns [41], [42], and then scaled to the definition of the NFI in 10% classes.

• Stand disturbances was also considered from another perspective. We compared the 95th percentile of CHM heights (\(H_{95}\)) to the dominant height measured in the field (\(H_d\)). The distribution of residuals \(H_{95} - H_d\) is positively skewed: we applied an asymmetric rejection criterion: plots were excluded if \(H_{95} - H_d < -2\sigma\) or \(H_{95} - H_d > +3\sigma\), where \(\sigma\) denotes the standard deviation of the residual distribution (\(\sigma = 3.89\:\)m). The stricter lower bound reflects the fact that a strong underestimation of \(H_{95}\) relative to \(H_d\) is most likely indicative of disturbances between the two acquisition dates, rendering the plot unsuitable for model training. Conversely, moderate overestimation of \(H_{95}\) relative to \(H_d\) can arise from legitimate physical sources: steep terrain slopes are known to introduce height errors in CHM normalization [43], and a sparse ground point density during DTM interpolation may lead to residual normalization biases, both resulting in inflated canopy height estimates. These artifacts do not necessarily reflect a mismatch between the field survey and the LiDARacquisition, and we wish to remain as faithful as possible to the field measurements rather than aggressively discarding potentially valid plots. Therefore, overestimation was tolerated up to a wider margin of \(+3\sigma\).

2.2 Methods↩︎

FLORAis designed to estimate various forest attributes at the plot level using a sufficient diversity of LiDARpoint clouds and that should be agnostic to the forest type. This is a regression task that is efficiently handled with machine learning solutions and that does not require any preliminary individual tree crown detection step [44].

2.2.1 Model architecture↩︎

2.2.1.1 Computational strategy

Our problem can be cast as an harmonization challenge, it should therefore exhibit a strong generalization ability over space, time, and plot composition. We selected an Octree-based Convolutional Neural Network combined with a High-Resolution Net (OCNN-HRNet) architecture developed by [45], and already proved effective for forest attribute regression by [30], while remaining simple. The proposed framework is designed to operate directly on raw point cloud data, as the conversion into octree structures is automatically performed as a preprocessing step prior to model training (Fig. 3). To improve training stability and accelerate convergence, each target variable is normalized to zero mean and unit variance prior to training using a z-score transformation [46]. This is motivated by the large differences in magnitude and variance across the six forest attributes (e.g., volumes are expressed in m\(^{3}\) ha\(^{-1}\) while basal area is in m\(^{2}\) ha\(^{-1}\)), which would otherwise cause the loss function to be dominated by the highest-variance attributes. Concretely, for each target \(k\), the normalized value is computed as \(\tilde{y}_k = (y_k - \mu_k) / \sigma_k\), where \(\mu_k\) and \(\sigma_k\) are the mean and standard deviation estimated over the training set. This ensures that all attributes contribute equally to the loss during training, regardless of their physical scale or absolute magnitude. Predictions are back-transformed to the original scale before evaluation.

2.2.1.2 Detailed design

The model architecture begins with feature extraction from octree-based inputs, transforming raw three-dimensional point cloud data into per-node feature representations (Fig. 3 part A). An initial stack of octree convolution and max-pooling layers reduces the spatial resolution while increasing feature dimensionality, producing the first set of feature maps at a coarser octree depth (Fig. 3 part B).

The core of the network follows a HRNet architecture adapted to octree-structured data (Fig. 3 part C). At each stage, parallel branches process feature maps at different octree depths, corresponding to different spatial resolutions, using residual blocks with bottleneck convolutions. The first stage operates on a single high-resolution branch (128 channels). In the second stage, a second branch is introduced at half the resolution (256 channels), and in the third stage a third branch is added at an even coarser resolution (512 channels). Between stages, a transition module simultaneously exchanges information across all branch pairs: each branch receives contributions from all other branches in parallel, computed from the same input state. These contributions are resampled to the target octree depth via successive max-pooling (for higher-resolution sources) or nearest-neighbour upsampling (for lower-resolution sources), with a \(1{\times}1\) convolution applied when channel dimensions differ, and summed element-wise to form the updated branch input.

After the HRNet stages, the multi-resolution feature maps are concatenated (Fig. 3 part D), resulting in a combined representation of 896 channels (128 + 256 + 512). A 1\(\times\)​1 convolution and global pooling layers then aggregate spatial information into a single 1024-dimensional feature vector, with dropout applied for regularization. Finally, in the original model from [30], a multilayer perceptron regressor composed of two linear layers with ReLU activations processes this vector and outputs the final regression estimates.
We propose several improvements with respect to existing works on deep-based regression of forest attributes : network optimization, increased versatility and context-aware estimation. First, we perform a systematic hyperparameter search to adapt the original architecture to our prediction task. Our six targets are continuous forest structural attributes of heterogeneous scale and ecological meaning: the optimal learning rate, weight decay and batch size may differ substantially from those reported in [30] that focuses on above-ground biomass estimation. Second, we design the model to rely exclusively on the 3D spatial coordinates (XYZ) of each point, deliberately excluding intensity and echo-related attributes such as return number, number of returns, and scan angle. Although these features were initially considered, hyperparameter optimization revealed no consistent performance gains from their inclusion, while intensity measurements were discarded due to the lack of radiometric sensor calibration. Restricting the input to XYZ coordinates also improves robustness across acquisition campaigns by reducing sensitivity to sensor-specific radiometric differences. Finally, we augment geometric information with contextual auxiliary data describing stand-level geographical, semantic, and temporal conditions. We hypothesize that this additional information enables the model to capture ecological and phenological significant variability that is not discernible from point cloud geometry alone. Specifically, we encode (see Fig. 3):

• Geographical-semantic variables (categorical, 20): canopy cover class, composition type, and sylvo-ecoregion using learnable embedding layers.

• Geographical knowledge (discrete, 4): we further include spatial context through a low-dimensional encoding of the plot center location.

• Temporal knowledge (discrete, 3): LiDARacquisition date is encoded with a cyclic sinusoidal representation (sine and cosine of the day of year), while the temporal gap between the LiDARacquisition and the field inventory date as a normalized scalar.

Land cover and type of forest formation are obtained from the BD Forêt®  v2 database. BD Forêt®  v2 is a vector database describing forest composition over France, and agregated into 32 main forest formations of pure and mixed compositions. Unlike field-based measurements (e.g., NFI plots), which are highly accurate but spatially sparse, BD Forêt®  offers wall-to-wall coverage suitable for large-scale modeling. The database includes a hierarchical classification that distinguishes closed and open forests according to canopy cover, as well as poplar plantations, heathlands, and herbaceous formations. Stand composition classes further differentiate pure broadleaf, pure conifer, and mixed forests.

All auxiliary data features are concatenated into a joint vector and passed through a dedicated multi-layer perceptron (MLP) to produce a learned auxiliary data representation in a shared latent space. Rather than modulating geometric features through feature-wise affine transformations like the FiLM method [47] (unsuccessful in our case), we adopt a late fusion strategy: the auxiliary data representation is concatenated with the features extracted by the HRNet backbone and jointly processed through subsequent layers (Fig. 3 part E). To allow the model to adaptively control the contribution of contextual information, we introduce a gating mechanism that learns, for each prediction target, how to interpolate between geometry-driven features and fused representations (Fig. 3 part F).

When the full set of auxiliary data is used, the resulting 27-dimensional representation is concatenated into a single feature vector. We adopt a late-fusion strategy, in which this auxiliary representation is combined with the geometric LiDARfeatures extracted by the HRNet backbone. To allow the model to adaptively modulate the contribution of contextual information depending on the target variable, we introduce a target-specific gating mechanism: \[\mathbf{g}_k = (1 - g_k) \odot \mathbf{L} + g_k \odot \mathbf{A}, \label{eq:gating}\tag{1}\] where \(\mathbf{L} \in \mathbb{R}^{128}\) denotes the latent representation derived solely from geometric LiDARfeatures, \(\mathbf{A} \in \mathbb{R}^{128}\) corresponds to the fused representation combining LiDARand auxiliary data. \(g_k = \sigma(\mathrm{MLP}([\mathbf{L}, \mathbf{A}])) \in [0,1]\) is the learned gating coefficient associated with target variable \(k\). When \(g_k \to 0\), the model relies exclusively on LiDAR-derived geometric information, whereas when \(g_k \to 1\), it fully incorporates auxiliary contextual information. This mechanism provides the model with partial interpretability, as the distribution of gate values on the test set enables quantification, for each forest attribute, of the relative contribution of contextual information to the final prediction. We hypothesize this design enables the network to selectively leverage auxiliary data when informative, while preserving robustness when such information noisy or less relevant.

2.2.1.3 Hyperparameter optimization

DL model performance is highly sensitive to the choice of hyperparameters, such as the learning rate, batch size, and regularization coefficients. Hyperparameter tuning is therefore formulated as an optimization problem in which hyperparameters are selected to minimize a validation loss. Common approaches include grid search, random search, and Bayesian optimization [48]. While grid search evaluates predefined combinations of hyperparameters, random search samples from predefined distributions and is generally more efficient in high-dimensional spaces [49]. Bayesian optimization further improves efficiency by modeling the relationship between hyperparameters and model performance using surrogate probabilistic models. Hyperparameter optimization was performed using Optuna [50], an open-source framework that supports dynamic construction of the search space during execution. Optuna relies on Bayesian optimization based on a Tree-structured Parzen Estimator, which adaptively samples promising regions of the hyperparameter space. In addition, Optuna integrates pruning mechanisms that enable early stopping of unpromising trials, making it particularly suitable for computationally expensive DL models [48], [51]. The optimized hyperparameters included the learning rate, batch size, dropout probability and weight decay coefficient. These hyperparameters were selected as they have the greatest influence on optimization dynamics (learning rate, batch size) and generalization (dropout, weight decay). The search space was defined using log-uniform sampling for the learning rate and categorical sampling for the remaining parameters.

The optimization procedure was implemented in two sequential phases. First, a fixed number of trials were executed. Due to the dataset size, 50 trials required 100 hours of computation. Although a more extensive exploration of the hyperparameter space would have been desirable, the high computational cost of this procedure led us to limit the search to 50 trials. Each trial consisted of training the model for a fixed number of epochs of 10 (preliminary results showed that models converge approximately at epoch 7). Candidate configurations were ranked according to their mean validation loss. The top three configurations were then selected for a second evaluation phase. In this phase, each configuration was retrained and validated using a five-fold cross-validation scheme and a larger number of training epochs (20). The final set of hyperparameters was chosen based on the lowest mean validation loss observed during this second phase. This two-stage strategy allows efficient exploration of the hyperparameter space while providing a robust and unbiased estimate of model generalization performance.

2.2.2 Study Design↩︎

The research questions outlined in Section 1 are addressed through four complementary experiments, each isolating a specific factor: the effect of phenological conditions on model generalization, the potential confound introduced by class imbalance, the role of overlapping multi-season acquisitions, and the contribution of auxiliary contextual variables together with the interpretability of their integration.

All models introduced below were built upon our FLORAarchitecture (see Section 2.2.1). They were trained in a supervised framework using paired LiDARpoint cloud samples and NFI ground truth measurement, for 20 epochs with optimized hyperparameters. Forest type was defined at the plot level from NFI records, distinguishing deciduous, coniferous, and mixed stands. Phenological classes were defined using fixed temporal boundaries applied to acquisition dates (leaf-start: Mar. 15 ; leaf-on: May 15 ; leaf-shut: Sep. 15 ; leaf-off: Nov. 15): leaf-oncorresponds to full canopy closure, and leaf-offto the leafless winter period. Transitional stages were excluded to ensure unambiguous class membership. leaf-onand leaf-offacquisitions represent 54% and 22% of the total dataset, respectively.

2.2.2.1 Experiment 1: All-seasons versus season-specific prediction

We assess whether a single model trained on both phenological seasons generalizes as well as season-specific models. Three independent models were trained: a leaf-onmodel (\(n = 14{,}142\)), a leaf-offmodel (\(n = 5{,}675\)), and a full model combining both seasons (\(n = 19{,}817\)). Each season-specific model was evaluated under two conditions: same-season testing (the test set matched the training season), and cross-season testing (the model was applied to acquisitions from the opposite season). The full model was evaluated on both seasons. To further disentangle the effect of class distribution from that of seasonal coverage, the full model was additionally trained on a balanced dataset enforcing equal leaf-onand leaf-offrepresentation through undersampling of the majority class; the balanced model and unbalanced original full model are compared to isolate the contribution of class imbalance (54% leaf-onvs.% leaf-off) from that of phenological diversity. For all configurations, samples were split into 70% training, 10% validation, and 20% test, stratified by phenological class and forest type.

2.2.2.2 Experiment 2: Benefiting from overlapping acquisitions

We quantify the impact of plots surveyed under both leaf-onand leaf-offconditions as a training strategy. Two datasets of equal size were compared: one composed exclusively of overlapping plot pairs (i.e., plots with acquisitions in both seasons), and one restricted to single-acquisition plots. Overlapping plots arise from two sources: spatial overlap between adjacent acquisition blocks, and flight-line sidelap within a single block (when point clouds contained data from multiple flights, points were separated using GPS timestamps). Among the \(1{,}239\) NFI plots acquired more than once across the full dataset, only \(38\) were observed in both phenological seasons, yielding \(76\) paired acquisitions (\(38\) leaf-onand \(38\) leaf-off). The single-acquisition dataset was constructed as \(K = 3\) disjoint subsamples of equal size, drawn without replacement from the pool of non-overlapping plots, ensuring no plot appeared in more than one subsample. Both conditions were evaluated on a fixed holdout test set composed exclusively of single-acquisition plots. Performance is reported as mean \(\pm\) standard deviation of rRMSE% across the \(K = 3\) runs.

2.2.2.3 Experiment 3: Assessing the contribution of auxiliary contextual information

To assess the contribution of auxiliary data, two models were compared under otherwise identical conditions: our full model incorporating ecological and spatiotemporal data, fused via the late-fusion gating mechanism (Aux), and a base model consisting of FLORAmodified to rely solely on LiDARpoint cloud features (No Aux). The gating coefficient \(g_k \in [0,1]\) learned for each target variable \(k\) quantifies the relative contribution of contextual information to the final prediction: values close to 0 indicate the model relies primarily on LiDAR-derived geometric features, while values close to 1 indicate a dominant contribution from auxiliary data. Since the gates are initialized at \(0.5\), values remaining close to this midpoint suggest that the auxiliary information neither contributed positively nor negatively in a substantial way. Analyzing the distribution of gate values on the test set thus provides a form of partial interpretability, enabling us to characterize, per forest attribute, how much the model draws on contextual information. This is complemented by an ablation protocol in which individual auxiliary variables are successively selected to quantify their marginal contribution. Both model variants were optimized independently using Optuna (50 trials), with the best hyperparameter configuration selected by five-fold cross-validation based on the lowest mean validation loss.

3 Results↩︎

3.1 Hyperparameter tuning↩︎

Hyperparameter optimization revealed that integrating auxiliary data fundamentally reshapes the optimization landscape, requiring a learning rate 3.5 times smaller, half the batch size, and double the dropout rate compared to the base model, while simultaneously demonstrating greater representational richness through resilience to aggressive regularization.

This part shows the results of the hyperparameter optimization using the aforementioned setup, conducted independently for the base model and the full model. Notably, the two variants converged to markedly different configurations, which reflects the impact of auxiliary data on the model’s learning dynamics. The learning rate converged to a value approximately 3.5 times lower in the auxiliary model (\(4.457 \times 10^{-4}\) vs.\(1.570 \times 10^{-3}\)), reflecting the more complex loss landscape induced by the fusion of heterogeneous information sources: the gating mechanism must learn to balance two complementary streams and is consequently more sensitive to large gradient steps. The optimal batch size likewise halved (16 vs.): smaller batches produce noisier gradient estimates that tend to converge towards flatter, more generalizable minima [52]. This is particularly beneficial here as the auxiliary metadata introduces a combinatorial diversity of input configurations (e.g., rare sylvo-ecoregion and acquisition period combinations), encouraging the model to find solutions that generalize across this expanded input space rather than overfitting to the most frequent covariate patterns. The most striking difference concerns the dropout rate, which doubled when auxiliary data was introduced (0.6 vs.): rather than viewing this solely as a defensive regularization measure, this result can be read as evidence of the representational richness brought by auxiliary data: the model remains performative despite a substantially higher fraction of masked activations, suggesting that the auxiliary features provide sufficiently redundant and complementary signals for the network to tolerate aggressive information dropout, while simultaneously preventing the model from short-circuiting point cloud processing by treating high-level contextual features, such as forest type, geographic region, or acquisition date, as direct proxies for the target variable. Conversely, weight decay was reduced by a factor of approximately three (\(3.6 \times 10^{-5}\) vs.\(1 \times 10^{-4}\)): since the stronger dropout already provides substantial regularization, a lower \(\ell_2\) penalty avoids over-constraining the learnable embeddings associated with categorical metadata, which require sufficient parametric freedom to encode meaningful distinctions between categories. Taken together, these results paint a coherent picture of how auxiliary data reshapes the optimization problem: it substantially enriches the representational space, as evidenced by the model’s resilience to aggressive dropout, while the added complexity of multimodal fusion calls for more conservative gradient updates, and the combination of smaller batches and reduced weight decay ensures that this richer representation generalizes effectively across the full diversity of input configurations.

3.2 All-seasons versus season-specific prediction↩︎

3.2.1 Results for full model↩︎

The full model achieves strong performance on height and total volume, while stem density remains the most challenging attribute (R\(^2\) < 0.47), with species-specific volume predictions exhibiting high variability across forest types.

Table 2 presents the global performance of the full model across all forest plots, evaluated on both leaf-onand leaf-offtest sets. Details for forest types are provided in 8.

Height is the best-predicted attribute by a large margin, with error below 13% and R\(^2\) \(\,>\,\)​0.89 on leaf-on, and remains stable across forest types. Total volume also yields strong results (error \(\approx\) 40%, R\(^2\) = 0.76 on leaf-on), with consistent performance across forest types. Basal area achieves moderate accuracy (error \(\approx\) 36%, R\(^2\) = 0.56 on leaf-on), while stem density is the most poorly predicted attribute, with R\(^2\) below 0.47 on both test sets.

Predictions for the proportion of deciduous and conifers exhibit the highest errors and the widest per-forest-type variability. Coniferous and deciduous volumes reach errors of 53.95% and 77.23% on leaf-onrespectively, with extreme rRMSE values on ecologically mismatched forest types (e.g., deciduous volume on coniferous plots), where target values are near zero. R\(^2\) values nonetheless remain relatively high for deciduous volume (0.82–0.83), reflecting good rank-order predictions despite high absolute errors.

With the exception of density, and volume of deciduous, other attributes leaf-offerrors are generally lower than or comparable to leaf-on, while R\(^2\) values tend to be higher on leaf-on. This trend indicates that leaf-oncapture better their overall variance, but is impacted by larger prediction errors.

When evaluated on the combined All test dataset, performance remains highly consistent with the separate leaf-onand leaf-off evaluations, with only marginal changes in both rRMSE and \(R^2\) across all targets. Overall trends observed in the separate test sets are preserved, with height and total volume remaining the most accurately predicted, while stem density and species-specific volume fractions continue to be the most challenging targets.

3.2.2 Balanced↩︎

Balancing the training set leads to moderate performance degradation on leaf-onwhile largely preserving leaf-offperformance, suggesting that the original training imbalance does not critically bias the full model.

Figure 4 compares the full and balanced models across all target attributes. Balancing the training set leads to a moderate degradation of leaf-onperformance across most attributes, while leaf-offperformance is largely preserved or marginally improved. The effect is most pronounced for coniferous volume and stem density, where leaf-onerror increases by 6 and 3 points respectively, and for height, where leaf-onR\(^2\) drops from 0.890 to 0.880. In contrast, height and basal area show the smallest sensitivity to the balancing strategy, with differences below 1 point of error on both test sets. Overall, the performance gap between balanced and unbalanced models remains limited, suggesting that the initial imbalance does not critically bias the full model. Detailed per-forest-type results are reported in 9.

3.2.3 Season-Specific Models↩︎

Season-specific models consistently underperform the full model even within their calibration domain, with cross-season generalization exhibiting strong asymmetry: the leaf-onmodel transfers reasonably well to leaf-off data, while the leaf-offmodel degrades severely when applied to leaf-onacquisitions.

Table 3 compares season-specific and full model performance under same-season and cross-season evaluation. The full model has better performance than the seasonal one when \(\Delta\)rRMSE is positive and \(\Delta\)R\(^2\) negative.

3.2.3.1 Same-season performance

Models trained on either leaf-onor leaf-offdata show slightly lower performances than the full one for both rRMSE and R\(^2\) (Tab. 3). This trend applies to both leaf conditions, indicating that the full model performs better than the specific ones, even in their calibration domain. For leaf-on, differences remain modest: volume error increases by less than 1 point (41.80% vs. 40.83%), and stem density and basal area stay within 1% of the full model. The degradation is more visible for coniferous volume (57.49% vs. 53.95%) and deciduous volume (84.10% vs. 77.23%).

On leaf-off, the gap widens. While basal area and stem density remain comparable (35.79% vs. 34.03% and 58.67% vs. 57.44%), deciduous volume shows a pronounced degradation (118.77% vs. 96.25%) and height drops substantially in both error and explained variance (16.08%, R\(^2\) = 0.626 vs. 11.44%, R\(^2\) = 0.811), suggesting that training on a single season is particularly limiting for structurally sensitive attributes.

3.2.3.2 Cross-season performance

Cross-season generalization is strongly asymmetric between the two seasonal models. On one hand, the leaf-onmodel applied to leaf-offacquisitions maintains acceptable performance for most attributes: The rRMSE of volume increases by 3.22% and R\(^2\) decreases by 0.07, staying close to same-season performance, and height remains robust with less than 1% error difference. The main degradation in error is observed for deciduous volume (+33.15%) and coniferous volume (+7.24%). On the other hand, the leaf-offmodel applied to leaf-onacquisitions degrades severely across all attributes. Volume rRMSE increases by 24.32%, coniferous volume by 45.13%, and deciduous volume collapses by 119.06%. Basal area R\(^2\) drops by 0.4 and stem density R\(^2\) by 0.32, while height, though degraded in error (+5.79%), partially retains its R\(^2\) (-0.03). This asymmetry suggests that leaf-offacquisitions carry less transferable structural information than leaf-ondata when applied outside their training distribution.

3.2.3.3 Comparison with the full model

The full model consistently outperforms season-specific models in both settings. In same-season testing, its advantage is clearest on leaf-off, where it reduces volume error by more than 3 points and nearly halves the height error gap (11.44% vs. 16.08%, R\(^2\) = 0.811 vs. 0.626). In cross-season testing, the full model avoids the severe degradation of the leaf-off-specific model: it maintains volume error at 40.83% on leaf-onwhere the leaf-off-specific model reaches 65.15%, and preserves R\(^2\) above 0.75 for volume and height across both test sets. per-forest-type results are reported in 10.

3.3 Can we benefit from overlapping acquisitions ?↩︎

Spatial overlap between LiDARacquisitions produces contrasting effects depending on the target attribute: volumetric estimation improves substantially for total and deciduous volume, while coniferous volume, basal area, and density degrade, and height remains largely unaffected.

Overlapping acquisitions yield contrasting effects depending on the target attribute (Figure 5 and 11). Performance improves most clearly for total and deciduous volume, with rRMSE reductions of roughly 5 and 12 points respectively, suggesting that the complementary seasonal information captured by paired acquisitions benefits volumetric estimation of broadleaf stands. Coniferous volume shows the opposite trend, with rRMSE increasing from \(73.15\%\) to \(80.87\%\) and R² dropping from \(0.441\) to \(0.324\). Basal area and density also degrade under the overlap condition, while height remains largely unaffected in either direction. Compared to the full model, results in this experiment are strongly degraded for all target variables (e.g., \(\approx\)+10 points in volume error) probably due to the reduced train dataset size.

3.4 Do auxiliary contextual information bring additional knowledge ?↩︎

Auxiliary data yield limited improvements across most attributes, with species-specific volumes showing the highest gate activations and largest rRMSE reductions, while height demonstrates no improvement linked to auxiliary data and density systematically degrades.

3.4.0.1 Impact on predictive accuracy (Figure 6)

Figure 6 compares global rRMSE and R\(^2\) obtained with LiDAR-only inputs (grey bars) to the FLORAarchitecture including auxiliary data (black bars), both trained on the full training set. Detailed numerical results are provided in 8. Across attributes, changes remain moderate in absolute magnitude, but clear patterns emerge. For total volume, rRMSE decreases by approximately 0.7%, while R\(^2\) remains stable. Coniferous volume shows slightly larger effects, with rRMSE reductions of 2% and R\(^2\) gains reaching +0.02 while performance on deciduous volume are stable. Basal area (rRMSE -1.21, R\(^2\) +0.03) and Density (rRMSE -1.14, R\(^2\) +0.02) exhibit small improvements. In contrast, height show minimal sensitivity to auxiliary data: rRMSE variations remain under 0.2%, with R\(^2\) changes below 0.01. Overall, auxiliary data lead to limited improvements, but more noticeable for species-specific volume components.

Figure 7 reports the median gate values and associated performance differences for each auxiliary variable and target attribute. Full results are reported in 12.

Species-specific volumes benefit most from auxiliary information. Coniferous volume consistently exhibits the highest gate activations across all configurations (\(0.66\)–\(0.89\)), regardless of which variable is provided, with rRMSE reductions reaching \(-5.73\) for the time delta between LiDARand field acquisition alone and \(-6.20\) for the full auxiliary set. Deciduous volume follows a similar pattern, with gate values between \(0.63\) and \(0.78\) and consistent improvements up to \(-1.99\) points for center. Total volume and basal area show intermediate activations (\(\approx\) \(0.47\)–\(0.66\)) and more modest but generally positive contributions. Height displays the most heterogeneous gating behavior, ranging from 0.12 (date) to 0.79 (cover), yet without consistent performance benefit, high gate values do not systematically translate into lower error, suggesting the model gates toward auxiliary data without effectively exploiting it for this attribute. Density is the only attribute for which auxiliary information systematically degrades performance across configurations, with rRMSE increases up to \(+1.35\) for center, despite gate values remaining near 0.5.

Among individual variables, center and delta yield the strongest marginal contributions for volume-related attributes, while date and cover have negligible or inconsistent impact across most targets.

4 Discussion↩︎

4.1 Full model performance↩︎

The full model achieves competitive results for height and total volume, which are the best-predicted attributes across all configurations (rRMSE \(<\) 13% and R\(^2\) \(>\) 0.88 for height; rRMSE \(\approx\) 39% and R\(^2\) = 0.74 for total volume on the combined test set). For height, these results are consistent with recent studies using ALS data at national scale in Sweden: [14] report stand-level rRMSE of 9.8–11.2% using linear regression on 11,500 NFI plots, and [53] report approximately 10–12% using XGBoost and CNN on NFI plots. The differences on total volume is more pronounced: at the local scale, [54] report rRMSE of 24.16–31.69% over 360 stands within a single Polish forest district using ALS combined with Landsat composites and various machine learning methods, while at the national scale, [14] report rRMSE of 17.2–22.0% and [53] approximately 30%. Both papers were conducted on Swedish boreal forests where Norway spruce and Scots pine together account for approximately 80% of standing volume. Conversely, our dataset comprises 33,877 NFI plots spanning 60 main species with \(\approx\) 72% of plots dominated by deciduous stands, making national-scale prediction substantially more challenging.

This spatial contrast is further supported by [55], who show strong regional variability in Sentinel-2 kNN-based forest estimates over 151,000 NFI plots across 14 European countries divided in multiple ecological zones. In the area including Swedish boreal reference region, model performance is consistently higher than in most other areas, with rRMSE of 58.2% for volume, and moderate predictive skill (R\(^2\) = 0.54). In contrast, areas covering the French territory exhibit systematically higher errors, with volume rRMSE ranging from 61.7% to 74.6% confirming the impact of forest structure and composition on performance.

For basal area, our global results (rRMSE\(\approx\) 36.0%, R\(^2\) = 0.54 on the combined test set) align reasonably with recent literature. [56] in Romania report, at 10\(\:\)m resolution, an R\(^2\) of 0.28–0.35 and rRMSE of 21.7–32.3%. Over Sweden, [14] achieve rRMSE ranging between 20.4–26.7%. At an even larger scale, [57] mapped structural attributes across the  552 million hectares of Canadian boreal forest using 30\(\:\)m annual Landsat composites (1984–2012 time series) combined with LiDARplots from a national transect survey and a kNN-Random Forest imputation approach, reporting R\(^2\) = 0.544 and rRMSE of 50.7% for basal area. Our results, obtained with national scale high-density LiDAR, thus compare favourably with large-area mapping studies relying on coarser remotely sensed data, while remaining above the accuracy levels reported for local or compositionally simpler settings.

Stem density is the most poorly predicted attribute across all configurations (RMSE \(\approx\) 55%, R\(^2\) \(<\) 0.47), a result consistent with the broader literature. For example, in a subtropical planted forest in southern China, [58] found that stem density was among the most challenging forest structural parameters to estimate from airborne LiDARdata, despite using a hybrid deep-learning model (Deep-RBN) that outperformed conventional approaches. Stem density predictions exhibited substantially lower accuracy than estimates of tree height and volume, with an rRMSE =20.34%–31.03% and R\(^2\) = 0.52–0.72 under the optimal model configuration. In the two studies focusing on high-density ALS in New Zealand plantation forests of radiata pine, stem density consistently yields the lowest performance among all forest attributes considered. Using standard ALS metric-based regression models, [59] report a R\(^2\) of 0.48 and rRMSE of 34% for stand density, compared to 0.88 for mean top height, noting explicitly that ALS models produced poor estimates of stand density. Similarly, using Random Forest regression, [60], in the southeast of New Zealand, achieve an R\(^2\) of 0.47 and rRMSE of 32.3% for stem density, versus 0.86 for height. Unlike height or volume, stem density is strongly influenced by stand history, thinning regime, and understorey structure-information that is poorly encoded in plot-level point cloud geometry and that auxiliary data cannot fully compensate for, as further evidenced by [61], who attribute the comparatively poor density predictions of their k-NN model to the structural complexity of broadleaved-dominated forests in central France.

To our knowledge, no existing study has directly predicted absolute coniferous and deciduous volumes as separate target variables. The closest comparable work is that of [55], who produced Pan-European maps of deciduous-coniferous proportion (DCP). While their DCP metric provides the proportion of conifers (allowing broadleaf proportion to be derived as 100–DCP), the absolute volumes of each species group remain unquantified. In the French territory covered by their study, the rRMSE for conifer proportion ranged from 61.8% to 71.6% on the processing area corresponding to France. This finding underscores the inherent difficulty of species-specific attribution. Our results for coniferous volume (RMSE \(\approx\) 50%, R\(^2\) =  0.74) and deciduous volume (RMSE \(\approx\) 82%, R\(^2\) =  0.82) showcase the challenging task of predicting absolute species-specific volumes. Despite the high error values, the R\(^2\) metrics remain reasonably strong, suggesting that the models capture spatial variability and stand composition patterns well despite substantial prediction errors at the plot level.

4.2 Effect of balanced training↩︎

Balancing the training set to equalize the representation of leaf-onand leaf-offacquisitions leads to a moderate and selective degradation of performance. On leaf-ontest data, most attributes show slightly higher errors under the balanced configuration, most notably for coniferous volume (+6 points rRMSE) and stem density (+3 points), while leaf-offperformance is largely preserved or marginally improved (Table in  6). Height and basal area show the smallest sensitivity to balancing, with differences below 1 point on both test sets (Figure 4). This pattern is consistent with the composition of the original dataset: leaf-onacquisitions represent the majority of plots, so balancing effectively down-samples the dominant season and reduces the model’s exposure to the most frequently occurring acquisition context. The limited overall impact suggests that the original imbalance does not critically bias the full model, and that the unbalanced configuration is preferable for operational use given its better exploitation of available training data. This is consistent with findings in class-imbalanced learning more broadly, where down-sampling the majority class tends to improve minority-class performance but at the cost of overall accuracy [62]

4.3 Phenological robustness and season-specific generalization↩︎

The results presented in Tables 3 and 10 address two complementary questions: whether season-specific models outperform the full model on their own acquisition season, and whether they generalize adequately to the opposite season.

4.3.1 Same-season performance and the advantage of leaf-off↩︎

In same-season evaluation, the leaf-offmodel achieves lower rRMSE than the leaf-onmodel for most volumetric attributes (e.g., total volume: 39.67% vs. 41.80%), consistent with the literature showing that leaf-offacquisitions generally yield more accurate volume estimates due to greater canopy penetration and improved ground sampling [8], [10]. This advantage is however attribute-dependent: height estimation degrades substantially under leaf-offsame-season conditions (16.08% vs. 13.98% for the leaf-on model), likely because the absence of foliage reduces the contrast between canopy strata and lowers the upper return density that drives height predictions. This is consistent with findings by [11], who show that leaf-offconditions shift the return distribution toward lower canopy levels, particularly in deciduous stands.

The leaf-offmodel also shows marked degradation for deciduous volume in same-season evaluation, with rRMSE dropping from 118.77% to values exceeding 150% without auxiliary data. This reflects a well-known difficulty: under leaf-offconditions, mixed stands containing both deciduous and evergreen species present ambiguous return distributions, as retained conifer foliage introduces structural noise that complicates the separation of species-specific contributions [8]. [10] report similar patterns, noting that leaf-offmodels show larger errors for deciduous volume attributes, particularly for merchantable volume, where rRMSE can differ by up to 7% between leaf conditions.

4.3.2 Cross-season generalization and asymmetric transferability↩︎

Cross-season generalization is strongly asymmetric. The leaf-onmodel applied to leaf-offacquisitions maintains acceptable performance for most attributes: total volume reaches 39.43% rRMSE (R\(^2\) = 0.60), close to its same-season performance, and height remains robust at 12.22% (R\(^2\) = 0.78). By contrast, the leaf-offmodel applied to leaf-onacquisitions degrades severely: total volume increases to 65.15% (R\(^2\) = 0.38), coniferous volume to 85.96%, and deciduous volume collapses to 196.29% (R\(^2\) = \(-\)​0.19). This asymmetry is consistent with [10], who report that applying leaf-offmodels to leaf-ondata on eight forest attributes (top height, mean height, Lorey’s mean height, basal area, quadratic mean diameter, merchantable volume, total volume, and total aboveground biomass) results in average absolute rRMSE increases of 16.8% for coniferous stands and 25.8% for deciduous stands, and with [8], who observe that a leaf-offmodel applied to leaf-ondeciduous plots yields rRMSE exceeding 70% and bias approaching 50%. The physical explanation is straightforward: a model trained on leaf-offpoint clouds, where returns are concentrated in the lower canopy and on the ground, lacks the representational capacity to interpret the denser and higher return distributions, representative of leaf-onacquisitions [7].

4.3.3 Practical implications for acquisition planning↩︎

These results raise a practical question for operational forest inventory: when is it preferable to acquire leaf-onversus leaf-offdata? The answer depends on the intended application. Leaf-offacquisitions are generally preferred for terrain modeling and archaeological prospection, where reduced foliage improves ground return detection and the visibility of micro-topographic features beneath forest canopies [63], [64]. Conversely, leaf-onacquisitions provide a more realistic characterization of canopy fuels and vegetation continuity, making them more suitable for wildfire-related mapping applications [65]. For multi-attribute forest inventory, our results suggest that leaf-onacquisitions are preferable when cross-season generalization is required, given their superior transferability to leaf-offconditions. However, acquisition season is not always a free parameter: mountainous areas are often only accessible under leaf-onconditions due to snow cover precluding leaf-offsurveys. Regulatory constraints can further restrict the acquisition window in some regions. These operational realities motivate the development of models capable of handling heterogeneous acquisition conditions, precisely the challenge that the full FLORAmodel addresses.

4.3.4 The full model as a robust operational solution↩︎

The full model, trained on both leaf-onand leaf-offdata, consistently outperforms season-specific models regardless of the test set phenology. In same-season evaluation, its advantage is clearest on leaf-off, where it reduces total volume rRMSE from 39.67% to 36.21% and nearly halves the height error gap (11.44% vs. 16.08%, R\(^2\) = 0.81 vs. 0.63). In cross-season setting, it avoids the severe degradation of the leaf-off-specific model, maintaining volume rRMSE at 40.83% on leaf-on,where the leaf-offmodel reaches 65.15%. This is consistent with [10], who show that pooled models trained on both seasons achieve intermediate performance between leaf-onand leaf-offmodels in same-season evaluation, while substantially outperforming cross-season model-mixing. Our results extend this finding to a DL setting at the national scale, suggesting that training on heterogeneous phenological conditions is not only feasible but produces more robust predictions than any season-specific strategy, a practically important result given the uneven and fragmented nature of leaf-on/leaf-offcoverage across the French national LiDAR programme.

4.4 Contribution of overlapping acquisitions↩︎

Overlapping plots i.e.,  surveyed under both leaf-onand leaf-offconditions, provide a unique training signal by anchoring the model to paired observations of the same forest under contrasting contexts. The results in Table 9 and Figure 5 show that exploiting these overlaps yields contrasting effects, depending on the target attribute. Performance improves most clearly for total and deciduous volume, with rRMSE reductions of approximately 5 and 12 points respectively. This suggests that paired phenological observations help the model better disentangle broadleaf structural signals across seasons. Height also shows modest but consistent improvements under the overlap condition. These gains are coherent with the premise that co-registered multi-seasonal acquisitions provide complementary canopy penetration profiles that are informative for broadleaf volume estimation [7], [11].

Coniferous volume shows the opposite trend, with rRMSE increasing from 73.15% to 80.87% and R\(^2\) dropping from 0.441 to 0.324 under the overlap condition. Basal area and stem density also degrade slightly. This result is counterintuitive but may reflect the limited size and geographic scope of the overlap subset: overlapping areas represent a geographically constrained fraction of the territory, and the model trained on these plots may lack sufficient diversity in coniferous stand types to generalize. Furthermore, coniferous stands are structurally more stable across phenological conditions than deciduous stands [12], meaning that paired leaf-on/leaf-offobservations provide comparatively less complementary structural information for conifers. As a result, restricting training to overlap plots likely reduces effective sample diversity, without providing the compensating benefit of paired phenological contrast that drives improvements for broadleaf volume estimation. This suggests that exploiting dense optical time series (e.g., Landsat or Sentinel-2) can help mitigate these limitations by capturing phenological variability that single-date or acquisition-tied observations cannot resolve. Time series approaches improve the discrimination of forest types by separating phenology-driven spectral signals from structural variability, which is particularly relevant in mixed or conifer-dominated stands where seasonal contrast is weaker and less informative [66].

4.5 Added value of auxiliary variables↩︎

The relevance of contextual variables was evaluated both through direct Aux (FLORA)/No Aux comparisons (4 and 8) and through ablation gating analysis (Figure 7). Together, these results reveal that auxiliary metadata provides limited improvements in overall and same‑season conditions, but substantial gains in cross‑season generalization, with strong contrasts across forest attributes and forest‑type strata.

The improvements become much larger in cross‑season settings. Both season‑specific models exhibit strong error inflation when tested on the opposite season, and in these cases auxiliary metadata plays a major stabilizing role. For example, the leaf-offmodel tested on leaf-onexhibits rRMSE reductions from 97.7% to 65.1% for total volume, from 121.9% to 85.9% for coniferous volume, and from 367% to 196% for deciduous volume. Similar effects are observed for the leaf-onmodel cross‑tested on leaf-off, with deciduous volume rRMSE improving from 168.8% to 129.4%. These results confirm that temporal and ecological metadata compensate for seasonal structural shifts—especially in leaf-offacquisitions where canopy geometry is less informative.

Coniferous volume and deciduous volume benefit the most from auxiliary information, consistent with their low geometric discriminability in the point cloud. Ablation gates identify delta, center and SER as the most influential variables, and the Aux/No Aux comparisons confirm this through systematic rRMSE and R\(^2\) gains for pure coniferous and pure deciduous stands. This is coherent with the strong spatial segregation of species types (e.g., coniferous stands at higher elevation), which the contextual variables make explicit. In contrast, mixed stands often degrade when auxiliary data is added, suggesting that ecological metadata pushes the model towards over‑interpreting coarse species patterns that do not reflect the finer structural mosaics of mixed forests.

Total volume and basal area show moderate gate activations and correspondingly modest but stable improvements, particularly through the center and ser variables. These two variables capture regional trends in productivity and forest structure that LiDARalone cannot fully encode. Their effects remain positive across most models and test conditions.

Height is the only attribute essentially unaffected by auxiliary metadata. Gate activations are heterogeneous but do not translate into performance gains/ The Aux/No Aux comparisons across all tables show changes below 0.2% rRMSE. This confirms that height is almost entirely determined by LiDARcanopy geometry and largely independent of ecological descriptors.

Stem density exhibits a more complex pattern. In same‑season tests, auxiliary variables mildly degrade performance, consistent with shortcut learning dynamics where categorical metadata encodes coarse density priors that interfere with the structurally weak LiDARsignal. However, in cross‑season tests, the effect reverses: density rRMSE improves substantially (e.g., 82.1% → 67.8% for the leaf-offmodel tested on leaf-on). Thus, while metadata can be detrimental in homogeneous acquisition contexts, it provides valuable stabilization when acquisition conditions shift, mirroring the behaviour observed for species‑specific volumes.

Overall, center, ser, and delta stand out as the variables providing the most reliable improvements—particularly for volume‑related attributes and in cross‑season scenarios—while date and cover contribute little or inconsistently. The convergence of direct performance comparisons and gating analysis highlights that auxiliary variables are not universally beneficial but become crucial whenever LiDARgeometry alone is insufficient to capture ecological variability or to bridge strong shifts in acquisition conditions.

5 Conclusion↩︎

This study introduced FLORA, a voxel-based DL framework for predicting six forest structural attributes from heterogeneous airborne LiDARdata at national scale, trained and evaluated on 32,052 NFI plots across mainland France.

The main findings are as follows. A single model trained jointly on leaf-onand leaf-offacquisitions consistently outperforms season-specific models, including within their respective calibration domains, and avoids the strong cross-season degradation observed when leaf-offmodels are applied to leaf-ondata. The initial class imbalance between leaf-onand leaf-offsamples does not bias the full model, and using all available training data proves preferable to enforcing class balance. Plots acquired under both phenological conditions improve estimates of broadleaf and total volume, while their benefit for coniferous stands is limited. Auxiliary contextual variables provide modest but consistent gains, particularly for species-specific volumes and cross-season generalization, whereas dominant height is driven almost entirely by point cloud geometry.

Overall, these results show that a generic DL approach trained on sufficiently large and diverse national-scale data can handle acquisition heterogeneity such as that of the French LiDARHD without explicit stratification by season, sensor, or forest type. The model also exhibits strong generalization, an increasingly critical requirement for large-scale applications where LiDARacquisitions are spatially heterogeneous, calling for methods capable of robust performance over extensive areas. Species-specific volume prediction remains the main limitation, highlighting the need to integrate dense optical time series to better capture phenological variation and enable operational wall-to-wall mapping.

6 CRediT authorship contribution statement↩︎

Emilie Vautier: Writing – original draft, Writing – review and editing, Data Curation, Conceptualization, Investigation, Methodology, Formal Analysis, Software, Validation, Visualization
Clément Mallet: Writing – review and editing, Conceptualization, Methodology, Formal Analysis, Supervision
Cédric Vega: Writing – review and editing, Conceptualization, Methodology, Formal Analysis, Supervision

7 Declaration of Competing Interest↩︎

The authors declare that there are no conflicts of interest associated with this work.

8 Detailed results of FLORAtrained on the full training set↩︎

4pt

9 Detailed results of FLORAcomparing balanced and non-balanced full training sets↩︎

4pt

4pt

10 Results for season-specific model↩︎

4pt

4pt

11 Overlapping acquisitions↩︎

4pt

12 Auxiliary data↩︎

3pt

3pt

References↩︎