Can regional aerial images from orthophoto surveys produce high quality photogrammetric Canopy Height Model? A single tree approach in Western Europe

: 11 Forest monitoring tools are needed to promote effective and data driven forest management and 12 forest policies. Remote sensing techniques can increase the speed and the cost‐efficiency of the 13 forest monitoring as well as large scale mapping of forest attribute (wall‐to‐wall approach). Digital 14 Aerial Photogrammetry (DAP) is a common cost‐effective alternative to airborne laser scanning (ALS) 15 which can be based on aerial photos routinely acquired for general base maps. DAP based on such 16 pre‐existing dataset can be a cost effective source of large scale 3D data. In the context of forest 17 characterization, when a quality Digital Terrain Model (DTM) is available, DAP can produce 18 photogrammetric Canopy Height Model (pCHM) which describes the tree canopy height. While this 19 potential seems pretty obvious, few studies have investigated the quality of regional pCHM based on 20 aerial stereo images acquired by standard official aerial surveys. Our study proposes to evaluate the 21 quality of pCHM individual tree height estimates based on raw images acquired following such 22 protocol using a reference filed‐measured tree height database. To further ensure the replicability of 23 the approach, the pCHM tree height estimates benchmarking only relied on public forest inventory 24 (FI) information and the photogrammetric protocol was based on low‐cost and widely used 25 photogrammetric software. Moreover, our study investigates the relationship between the pCHM 26 tree height estimates based on the neighboring forest parameter provided by the FI program. 27 Our results highlight the good agreement of tree height estimates provided by pCHM using


Abstract: 11
Forest monitoring tools are needed to promote effective and data driven forest management and 12 forest policies. Remote sensing techniques can increase the speed and the cost-efficiency of the 13 forest monitoring as well as large scale mapping of forest attribute (wall-to-wall approach). Digital 14 Aerial Photogrammetry (DAP) is a common cost-effective alternative to airborne laser scanning (ALS) 15 which can be based on aerial photos routinely acquired for general base maps. DAP based on such 16 pre-existing dataset can be a cost effective source of large scale 3D data. In the context of forest 17 characterization, when a quality Digital Terrain Model (DTM) is available, DAP can produce 18 photogrammetric Canopy Height Model (pCHM) which describes the tree canopy height. While this 19 potential seems pretty obvious, few studies have investigated the quality of regional pCHM based on 20 aerial stereo images acquired by standard official aerial surveys. Our study proposes to evaluate the 21 quality of pCHM individual tree height estimates based on raw images acquired following such 22 protocol using a reference filed-measured tree height database. To further ensure the replicability of 23 the approach, the pCHM tree height estimates benchmarking only relied on public forest inventory 24 (FI) information and the photogrammetric protocol was based on low-cost and widely used 25 photogrammetric software. Moreover, our study investigates the relationship between the pCHM 26 tree height estimates based on the neighboring forest parameter provided by the FI program. 27 Our results highlight the good agreement of tree height estimates provided by pCHM using DAP with 28 both field measured and ALS tree height data. In terms of tree height modeling, our pCHM approach 29 reached similar results than the same modeling strategy applied to ALS tree height estimates. Our 30 study also identified some of the drivers of the pCHM tree height estimate error and found forest 31 parameters like tree size (diameter at breast height) and tree type (evergreenness/deciduousness) as 32 well as the terrain topography (slope) to be of higher importance than image survey parameters like 33 the variation of the overlap or the sunlight condition in our dataset. In combination with the pCHM 34 tree height estimate, the terrain slope, the DBH and the evergreenness factor were used to fit a 35 multivariate model predicting the field measured tree height. This model presented better 36 performance than the model linking the pCHM estimates to the field tree height estimates in terms 37 of r² (0.90 VS 0.87) and RMSE (1.78 VS 2.01 m). Such aspects are poorly addressed in literature and 38 further research should focus on how pCHM approaches could integrate them to improve forest 39 characterization using DAP and pCHM. Our promising results can be used to encourage the use of 40 1. Introduction 46 Forests cover almost a third of global land area (Keenan et al., 2015). They provide numerous 47 ecosystem services and are of major importance in public policies worldwide. DAP and pCHM can thus be used to timely update regional forest 3D structure as aerial images are 93 generally acquired on a regular basis by national or regional mapping agencies. This interest is 94 reinforced as countrywide ALS surveys providing accurate DTMs are occurring more and more 95 frequently throughout the world. Another major interest of using DAP approaches based on stereo-96 images acquired during regular national campaigns is to value potentially very dense time series 97 which can cover time period which may cover periods prior to the acquisition of data ALS. 98 While the potential of using pCHM build with aerial images regularly acquired by national or regional 99 mapping agencies seems pretty obvious, few studies have investigated the quality of such 100 regional/countrywide pCHM, especially on the single tree scale. Evaluating the 3D accuracy of DAP 101 derived from such regional/national aerial survey is nevertheless an essential topic as such surveys 102 are generally not designed to produced high density 3D point clouds but only orhtophotomosaics 103 which typically require less overlap. In this context, Ginzler  smaller reference trees set (respectively 51 and 356 trees) even if the lack of harmonized accuracy 112 metrics hampers real accuracy benchmarking. None of the pre-cited studies investigated the 113 relationship between pCHM error at single tree level and the forest characteristics around the 114 considered trees (e.g., stem density, volume, basal area, canopy roughness). 115 In this context, we propose to evaluate the accuracy of individual tree height estimates provided by 116 pCHM build using aerial images acquired in the specific context of countrywide orthophoto survey 117 protocols. To further ensure the replicability of the approach, the pCHM tree height estimates 118 benchmarking only relied on public FI information and the photogrammetric protocol was based on 119 low-cost and widely used photogrammetric software. Moreover, our study investigates the 120 relationship between the pCHM tree height estimates based on the neighboring forest parameter 121 provided by the FI program in a European temperate forest context. 122 are also regularly found. Needle-leaved forests are largely composed of spruce (Picea abies) and 131

Material and Methods
Douglas fir (Pseudotsuga menziesii) and to a lesser extent, larches (Larix sp.) and pines (Pinus 132 sylvestris and P. nigra). The evergreen stands are mostly managed as even-aged stands and 133 harvested by the means of clear-cuts (Alderweireld et al., 2015). 134

Aerial surveys 135
The regional orthophoto surveys in the study area are achieved by private operators on a regular 136 basis, notably for the sake of controls related to European Union common agricultural policy. They 137 were initially acquired on a triennial basis and since 2015, on an annual basis. The timetable of the 138 regional surveys is driven by the objectives stated above and typically ranges from April to October 139 ( We also used a regional LiDAR survey performed from 12 December 2012 to 09 March 2014. This 147 LiDAR survey was used as a Digital Terrain Model (1m GSD) to compute the pCHMs (photo. DSMs -148 ALS DTM) as well an ALS CHM (1m GSD, ALS DSM -ALS DTM) to benchmark the 5 pCHMs in terms of 149 single tree height estimates. 150 We used Agisoft Metashape 1.5.4 in network mode for all the photogrammetric reconstruction steps 154 synthetized in Figure 1. Agisoft is one of the most used photogrammetric packages using a multiview 155 matching strategy (Smith et al., 2016). It also allows to handle large images (200 Mpx and more) 156 generated by large frame sensors like the UltraCam in an easy to set up network processing 157 interface. We choose this software for its relatively user friendly GUI as well as its rather low cost (ca. 158 549 $ for educational license and 3500 $ for commercial license) compared to other state-of-art 159 photogrammetric package like Trimble Inpho or Imagine Photogrammetry (LPS). These 160 characteristics will ease the reproducibility of the methodology. 161 We set up a processing network of 3 to 5 (depending on the resources available) computers 162 equipped with GPU processing (NVidia GTX) and 64 Go RAM. The very same photogrammetric 163 protocol was applied to process the raw images of the different regional orthophoto coverages. 164 Based on the GPS positions (metric coordinates system "Lambert 72", EPSG: 31370) and the camera 165 calibration information delivered by the service provider, we realized the tie point extraction in full 166 resolution using the "High" accuracy parameter, "Key points limit" and "Tie point limit" set to 40000 167 and 4000 respectively. To avoid overfitting, we followed the recommendations of James et al. (2017)  168 and used a rather conservative lens calibration strategy. The following lens calibration parameters 169 remained fixed: b1b2 (affinity and skew transformation coefficients), k4 and p3/p4 (additional 170 tangential and radial distortion coefficients). This lens calibration strategy was pursuit during the 171 entire photogrammetric processing. This first alignment process resulted in a first 3D reconstruction 172 (sparse point cloud) based on an initial bundle block adjustment (BA) as highlighted in Figure  The accuracy assessment process is based on a cross-validation using repeated k-fold technique (k=5, 186 repetition = 50) implemented in Agisoft Metashape through python scripting. For each of the 50 187 iterations, the GCP dataset is randomly divided into five folds of GCP and five BA processes are 188 successively run (Figure 2). Each BA process is run based on a set of GCP from one fold as checkpoint 189 (not used in the process, i.e. test set), the other GCP (associated to the 4 other folds) being used as 190 control points (to constrain the BA process, i.e. training set). The XYZ errors associated with the GCP 191 from the test fold are saved and another BA is run, using GCP's from one other fold as checkpoints. 192 Once all the GCP's from the k folds have been successively used as checkpoints (and thus as 193 independent test set), the entire process is repeated over 50 times to ensure the robustness of the 194 cross-validation. Indeed, all the GCPs of the network are used to assess the 3D accuracy of the 3D 195 reconstruction with 50 different neighboring conditions. As the accuracy of the 3D model is points. Once the accuracy assessment has been completed, all the GCP were used to run the final BA 203 process (Figure 1-2) with the GCP accuracy set to 0.01m in the Agisoft interface. 204 The dense matching process was run using aggregated images (aggregation factor of 2, half 205 resolution) as a compromise between spatial resolution and computing time. The "Aggressive" depth 206 filtering strategy was selected to limit the noise in the dense cloud. The final output is a DSM 207 ("Interpolation" enabled) which results in 5 regional DSM in raster format (Figure 1-4

pCHM Processing 228
The regional photogrammetric DSMs were combined with a regional ALS Digital Terrain Model 229 (photogrammetric DSM -ALS DTM) resampled according to the resolution of the DSM. As the ALS 230 survey occurred from 12 December 2012 to 09 March 2014 (see 2.2 section), we considered the 231 topography as constant during the entire study period. As our study is dedicated to trees located 232 inside forest landscapes, this assumption is reasonable considering the low erosion rate in forested 233 landscape as well as the infrequency in the study area of catastrophic events (e.g. landslide or 234 earthquakes). 235

Accuracy assessment of tree height estimates from pCHM 236
A selection of field reference plots from the regional FI program was performed using temporal and 237 spatial criteria. The FI plots collect various parameters such as tree height, diameter above breast 238 height (DBH), tree species, health condition, etc. The measurements of tree height were performed 239 with a vertex ultrasound instrument. Measured trees are spatially located using azimuth and distance 240 relative to the plot center which is georeferenced with off-the-shelf GNSS receivers. In terms of 241 absolute individual trees (Table 2) from 489 FI plots presenting an evenly spatial distribution in the study area. 264 To avoid 3D reconstruction issues linked to (partial) leaf-off conditions, FI plots were selected only 265 when they were in leaf-on condition during the associated aerial survey. In even-aged forests (mainly 266 spruce and Douglas fir), the field survey of the tree height is limited to few individuals in order to 267 compute the dominant height. As the evergreen forests are mainly managed as even-aged stands, 268 they occur in a rather low proportion in the reference tree database as compared to deciduous trees. 269 We performed the same matching between field measured tree heights from FI plots and tree height 270 estimates extracted from the ALS CHM. This database was composed of 1579 trees from 260 FI plots 271 (1144 deciduous, 435 evergreen). The Figure 3-III represents the FI plots network used in this study. 272 In order to test the reliability of tree height information provided by pCHM, we performed linear 273 modelling between field measured tree height (from the reference tree database) and tree height 274 estimates from local maxima of the 5 different pCHM's. We applied a similar method with ALS 275 reference dataset to benchmark tree estimates provided by the pCHM with a reference tree height 276 remote sensing data source. 277 278

Drivers of the pCHM tree height estimates error 280
To investigate the robustness of the pCHM tree height estimates, we compared the impact of various 281 parameters suspected to have an impact on the 3D reconstruction uncertainties or the tree height 282 estimate itself (see Table 3). Some parameters are assessed at the scale of the individual tree, others 283 at the forest inventory plot scale. The mean solar angle of the aerial images is a good proxy for the 284 light conditions during the surveys. Lower values can be associated to more important cast 285 shadowing and higher reconstruction uncertainties. The time difference between the aerial images 286 could also highlight artificial heterogeneity between the aerial images and hamper the 287 photogrammetric reconstruction. The number of overlapping images used for the photogrammetric 288 reconstruction is positively linked to the quality of the 3D reconstruction. As the 60% along-tracks 289 and 30% across-track overlap scheme induces varying image overlap condition, the number of 290 images used to reconstruct the forest canopy of the FI plot is an interesting parameter. The 291 uncertainties of the 3D photogrammetric reconstruction of tree canopy can also be linked to the 292 characteristic of the tree itself as well as its environment. The forest species and its evergreenness 293 were investigated as well as forest structure. Forest structure was here investigated in terms of stem 294 densities and tree size within the FI plot with the basal area. The relative DBH allows addressing the 295 size of the considered tree in relation with the size of the biggest trees in the FI plot (i.e. proxy of the 296 social status). Lastly, the terrain slope and altitude is also investigated as they can have a significant 297 impact on the 3D reconstruction but also on the accuracy of the ALS DTM itself. 298 To evaluate the significance of the linear relationship between the selected parameters and the tree 299 height absolute differences (abs( field tree height -pCHM tree height )), we used one-way analysis of 300 variance, (ANOVA) for qualitative factors, and linear regressions for quantitative factors. 301 We ran a best subset regression approach to build a multivariate linear model with the variables 302 previously highlighted using the regsubsets tools from the leaps package in R (Miller, 2017). We used 303 best subsets regression to fit all potential models of pCHM tree height estimate absolute error in 304 order to highlight the best combination of predictors (using bayesian information criterion, BIC) 305 Finally, we ran a second best subsets regression to test the potential of the highlighted parameters 306 from Table 3 to improve the tree height estimates with pCHM's data (using BIC). 307 308

325
The accuracy assessment of the tree height models highlights that the pCHM tree height estimates 326 agreed well with the field tree height estimates ( Figure 5) for the ALS CHM. The mean tree height estimate error is higher for the evergreen species as it ranges 339 from 2.5 to 3 m while it ranges from 0.92 to 2.02 m for deciduous species. The same trend can be 340 observed in the ALS models but again with lower mean error values (< 1 m). 341 342 Figure 5: biplots of tree height estimates from pCHM's and ALS CHM in comparison with FI tree height estimates.

344
Deciduous and evergreen trees are marked using " + " and " . " symbols respectively. The mean error (err. ) is computed 345 from the difference between the reference tree height (from FI) and the tree height estimates (provided by pCHM or ALS 346 CHM).

Drivers of the pCHM tree height estimates error 350
Among the 13 potential parameters listed in Table 3, our analysis highlighted 8 variables which  351 present a significant statistical link with the absolute error of the pCHM height estimates (Table 5). 352 Except the basal area and the number of trees in the FI plots, all the forest parameters were selected 353 as well as the topographic parameters (slope and altitude). In terms of correlation, the slope is the 354 only parameter which presents a negative correlation coefficient with the pCHM absolute error. 355 The reference year of the aerial survey is the only image parameter which was highlighted by our 356 analysis. As the number of levels in the factor is rather low for this variable, we ran a Tukey Honest 357 Significant

373
Within the set of variables which were proved to have a significant statistical relationship with the 374 pCHM tree height absolute error, we removed the reference year and the tree species factorial 375 variables before running the best subset regression process. This choice was made in order to ease 376 the interpretation (the * ( factor presents 29 levels in our dataset) and the replication of the 377 fitted pCHM error model. Indeed, the aerial survey reference gathers a bunch of environmental 378 parameters at the time of the flight survey (lights conditions, phenology …) with a rather low interest 379 for understanding the drivers of the pCHM tree height estimates error. 380 The best subsets regression model selected following the BIC a model with 2 variables: the DBH and 381 the evergreenness factor to predict pCHM tree height estimate error. This model presents a rather 382 low r² score (0.10) and a RMSE of 1.3 m. 383 Finally, we used the same initial set of variables used to fit the pCHM tree height estimate error 384 model to evaluate their potential income in terms of accuracy improvement of a global model linking 385 the tree height estimates provided with pCHM data (all pCHM used) and the FI tree height estimates. 386 The best subsets regression model selected following the BIC a model with 4 variables to predict the 387 field measured tree height: the pCHM tree height, the DBH, the evergreenness factor and the terrain 388 slope. This model is significantly different (ANOVA test, p-value < 0.001) of the linear model linking 389 pCHM tree height and field measured tree height. It presents a slightly higher r² score (0.90 VS 0.87) 390 and a smaller RMSE (1.78 VS 2.01 m). Therefore, the use of these variables in a multiple linear model