Impact of UAS Image Orientation on Accuracy of Forest Inventory Attributes

The quality and accuracy of Unmanned Aerial System (UAS) products greatly depend on the methods used to define image orientations before they are used to create 3D point clouds. While most studies were conducted in non- or partially-forested areas, a limited number of studies have evaluated the spatial accuracy of UAS products derived by using different image block orientation methods in forested areas. In this study, three image orientation methods were used and compared: (a) the Indirect Sensor Orientation (InSO) method with five irregularly distributed Ground Control Points (GCPs); (b) the Global Navigation Satellite System supported Sensor Orientation (GNSS-SO) method using non-Post-Processed Kinematic (PPK) single-frequency carrier-phase GNSS data (GNSS-SO1); and (c) using PPK dual-frequency carrier-phase GNSS data (GNSS-SO2). The effect of the three methods on the accuracy of plot-level estimates of Lorey’s mean height (HL) was tested over the mixed, even-aged pedunculate oak forests of Pokupsko basin located in Central Croatia, and validated using field validation across independent sample plots (HV), and leave-one-out cross-validation (LOOCV). The GNSS-SO2 method produced the HL estimates of the highest accuracy (RMSE%: HV = 5.18%, LOOCV = 4.06%), followed by the GNSS-SO1 method (RMSE%: HV = 5.34%, LOOCV = 4.37%), while the lowest accuracy was achieved by the InSO method (RMSE%: HV = 5.55%, LOOCV = 4.84%). The negligible differences in the performances of the regression models suggested that the selected image orientation methods had no considerable effect on the estimation of HL. The GCPs, as well as the high image overlaps, contributed considerably to the block stability and accuracy of image orientation in the InSO method. Additional slight improvements were achieved by replacing single-frequency GNSS measurements with dual-frequency GNSS measurements and by incorporating PPK into the GNSS-SO2 method.


Introduction
Forest inventory provides critical information about the condition and status of forest resources, which is important for the implementation of sustainable forest management. Mapping forest resources is possible with the use of remote sensing technology through continuous data acquisition and advanced statistical analysis. Although widely used, remote sensing-based research remains challenging due to different properties of the sensors and the structural complexity of forests [1,2]. Over the last twenty number of studies have evaluated the spatial accuracy of UAS products derived by using different image block orientation methods in forested areas [38,39]. Gašparović et al. [38] evaluated the vertical accuracy of the Digital Surface Model (DSM) generated from UAS images collected with the low-cost UAS (DJI Phantom 4 Pro) over a dense forested area. When GNSS-SO approaches with no GCPs and with seven irregularly distributed GCPs were compared, a considerable improvement in the DSM vertical accuracy with GCPs was observed. Gašparović et al. [38] concluded that DSMs generated using GNSS-SO without GCPs can be used for visualization and monitoring purposes, whereas DSMs generated using GCPs have greater potential to be used in forest inventory. Another recent study that evaluated the spatial accuracy of UAS products using different image block orientation methods in a forested area was conducted by Tomaštik et al. [39]. Namely, GNSS-SO using the PPK dual-frequency carrier-phase GNSS data collected with eBee Plus and without GCPs was compared to two InSO approaches, one with four and another with nine GCPs. Contrary to the previously mentioned studies, Tomaštik et al. [39] obtained a significantly higher horizontal accuracy for the GNSS-SO (PPK) method compared to both InSO approaches, and a significantly higher vertical accuracy for the GNSS-SO (PPK) method compared to InSO with four GCPs. The vertical accuracy for GNSS-SO and for InSO with nine GCPs did not differ significantly.
Besides the studies of Gašparović et al. [38] and Tomaštik et al. [39], who focused on DSM and PCs, to the best of our knowledge no prior studies have examined the accuracy of forest inventory attributes estimated using UAS products generated by using different image block orientation methods. The main goal of this study is to investigate the impact of several image block orientation methods on the accuracy of estimated forest attributes, with a particular emphasis on the plot-level mean tree height (Lorey's mean height). Although the accuracy of the plot-level mean tree height estimates was evaluated in a number of recent UAS-based studies [14,[40][41][42][43][44][45], the influence of image block orientation methods on the reported results has not been examined nor discussed. In this study, three image block orientation methods are used and compared: (a) InSO method with five irregularly distributed GCPs; (b) GNSS-SO method using non-PPK single-frequency carrier-phase GNSS data (hereafter referred to as GNSS-SO 1 ); and (c) using PPK dual-frequency carrier-phase GNSS data (hereafter referred to as GNSS-SO 2 ). Since previous studies [32,35,36] have demonstrated an improvement in image orientation and vertical accuracy of generated products (PCs, DSMs) when GCPs were used, both GNSS-SO approaches are supported with five GCPs in the present study. In combination with field forest inventory data, plot-level metrics from UAS PCs obtained using the three different methods (InSO, GNSS-SO 1 , GNSS-SO 2 ) are used to generate the plot-level mean tree height models which are then compared and evaluated using field ground-truth data. Two validation approaches are used: (1) hold-out validation (HV), i.e., field validation across independent sample plots, and (2) the leave-one-out cross-validation (LOOCV) statistical technique. The latter is commonly used for smaller datasets [46,47] when there are not enough measurements to divide them into modeling and validation subsets.

Study Area
The study was conducted in part of the state-owned, actively managed forests of Pokupsko basin forest complex, located in Central Croatia (Figure 1). The main forest type of this area is even-aged pedunculate oak (Quercus robur L.) forests of different age classes ranging from 0 to 160 years. These oak stands are mainly mixed with other tree species (Carpinus betulus L., Alnus glutinosa (L.) Geartn., and Fraxinus angustifolia Vahl.), and with understory vegetation (Corylus avellana L. and Crataegus monogyna Jacq.). The study area is characterized by flat terrain with ground elevations ranging from 108 to 113 m a.s.l.

Field (Ground-Truth) Data
The field data from a total of 99 circular sample plots were collected in 2017 ( Figure 1). The plots were systemically distributed (100 m × 100 m, 100 m × 200 m, 200 m × 100 m, 200 m × 200 m) throughout 20 forest stands of different age classes and with a total area of 363.54 ha. Depending on the stand age and stand density, the sample plots had a radius of 8, 15, or 20 m. The coordinates of the plot centres were measured during leaf-off conditions using the GNSS RTK method. The receiver Stonex S9IIIN, connected to the Croatian network of GNSS reference stations called Croatian Positioning System (CROPOS) and Very Precise Positioning System (VPPS) service, was used. CROPOS VPPS service, based on the Virtual Reference Station (VRS) technique, provides real-time corrections with up to 2 cm horizontal and 4 cm vertical accuracies [48]. Despite the method and service used, only 40% of the plot centres had a fixed solution and 60% of the plot centres had a float solution. The average precisions of the positioning, i.e., standard deviations (SDs) provided by the receiver, were 0.038 m and 0.159 m for fixed and float solutions, respectively.
The plot data were collected by recording the tree species and measuring the diameter at breast height (dbh) for all trees with dbh ≥ 10 cm. Tree height was measured for 1908 trees or 73% of all sampled trees in the study area using a Vertex III hypsometer. In each plot, the trees sampled for height measurements included all available species and a wide dbh range. According to the extensive research of Stereńczak et al. [49], the relative errors for tree height measurements using a Vertex instrument are 2.52% and −1.25% for Q. robur and A. glutinosa, respectively. Other tree species (C. betulus, F. angustifolia) found in the study area of the present research were not in the scope of the research by Stereńczak et al. [49], but it can be assumed that expected errors do not significantly deviate from those reported for Q. robur and A. glutinosa. For trees without height measurements, tree

Field (Ground-Truth) Data
The field data from a total of 99 circular sample plots were collected in 2017 ( Figure 1). The plots were systemically distributed (100 m × 100 m, 100 m × 200 m, 200 m × 100 m, 200 m × 200 m) throughout 20 forest stands of different age classes and with a total area of 363.54 ha. Depending on the stand age and stand density, the sample plots had a radius of 8, 15, or 20 m. The coordinates of the plot centres were measured during leaf-off conditions using the GNSS RTK method. The receiver Stonex S9IIIN, connected to the Croatian network of GNSS reference stations called Croatian Positioning System (CROPOS) and Very Precise Positioning System (VPPS) service, was used. CROPOS VPPS service, based on the Virtual Reference Station (VRS) technique, provides real-time corrections with up to 2 cm horizontal and 4 cm vertical accuracies [48]. Despite the method and service used, only 40% of the plot centres had a fixed solution and 60% of the plot centres had a float solution. The average precisions of the positioning, i.e., standard deviations (SDs) provided by the receiver, were 0.038 m and 0.159 m for fixed and float solutions, respectively.
The plot data were collected by recording the tree species and measuring the diameter at breast height (dbh) for all trees with dbh ≥ 10 cm. Tree height was measured for 1908 trees or 73% of all sampled trees in the study area using a Vertex III hypsometer. In each plot, the trees sampled for height measurements included all available species and a wide dbh range. According to the extensive research of Stereńczak et al. [49], the relative errors for tree height measurements using a Vertex instrument are 2.52% and −1.25% for Q. robur and A. glutinosa, respectively. Other tree species (C. betulus, F. angustifolia) found in the study area of the present research were not in the scope of the research by Stereńczak et al. [49], but it can be assumed that expected errors do not significantly deviate from those reported for Q. robur and A. glutinosa. For trees without height measurements, tree heights were estimated using the developed local species-specific dbh-height models fitted with Michailloff's function [50]. Despite the fact that field tree height measurements and the use of dbh-height models produce certain errors, the field data were considered as ground-truth for the validation of UAS data, as the method was common in many other remote sensing studies, e.g., [14,15,21,[23][24][25][26][27][28]32]. Lorey's mean height (H L ) of each plot was calculated using Equation (1): where h i is the tree height of a single (i) tree at the plot (in m), g i is the basal area of an ith tree at the plot calculated from measured dbh (in m 2 ), and G plot the plot basal area (in m 2 ). A summary of the main forest attributes within the surveyed plots is presented in Table 1. Out of 99 sample plots, 44 were located in oak stands of age class 3 (41-60-year-old), 17 in stands of age class 4 (61-80-year-old), 8 in stands of age class 5 (81-100-year-old), and 30 plots in stands of age class 7 (>121-year-old). The irregular distribution of sample plots per age class is a consequence of the irregularity of age classes in the study area. Moreover, it can be noted that stands of age class 6 (101-120-year-old) are not present in the study area.

UAS Data
The UAS survey, which included GCPs measurement and image acquisition, was conducted during leaf-on conditions on 30 May 2017. Prior to image acquisition, 7 GCPs were signalized with 40 × 40 cm markers and surveyed across the study area. The GCPs were surveyed using a Trimble R7 GNSS receiver connected with the CROPOS VPPS service, which ensured horizontal and vertical measurement precision of SD ≤ 10 cm. The distribution of the GCPs was somewhat irregular (Figure 1) due to a very dense forest structure and lack of open areas. In the processing phase, gross errors were detected for 2 GCPs and they were excluded from further processing.
The UAS images were collected using a fixed-wing Trimble UX5 HP equipped with a Sony ILCE-7R camera ( Table 2). The main advantage of the UAS used is its ability to record dual-frequency (L1 and L2 frequency) GNSS carrier-phase data. Both single-frequency and dual-frequency GNSS data were recorded using the same on-board GNSS receiver. Single-frequency data were recorded during flight in metadata files of each image. Dual-frequency data were recorded in RINEX (Receiver Independent Exchange Format) files. In total, 458 and 475 images with high overlap (90% endlap, 80% sidelap) [51] were collected during the first and the second flights, respectively ( Figure 1). The average flying height was 600 m above ground level, resulting in a Ground Sampling Distance (GSD) of about 8 cm. The UAS survey was conducted by a private company and GSD was selected to fulfil the requirements of the project given the large study area and cost constraints.

Digital Terrain Model (DTM) Data
To normalize the PCs of the UAS, a combination of an airborne LiDAR Digital Terrain Model (DTM LiD ) and an improved version of the official Croatian DTM (DTM PHM ) with a spatial resolution of 0.5 m was used. Namely, the minor, northern part of the study area (31 plots) was not covered with the DTM LiD , and therefore, an improved version of DTM PHM was used to cover this part of the area.
DTM LiD was generated from airborne LiDAR data collected with an Optech ALTM Gemini 167 scanner during summer 2016. The resulting point density was 13.64 points·m −2 . The "ground" points were classified using TerraSolid (version 11) software [52]. From the classified ground points (0.91 points·m −2 ), a raster DTM LiD with a spatial resolution of 0.5 m was generated. It was provided by Hrvatske Vode Ltd. (Zagreb, Croatia), while both data acquisition and image processing were done by the Institute for Photogrammetry Inc. (Zagreb, Croatia) and Mensuras Ltd. (Maribor, Slovenia).
DTM PHM was generated from official national digital terrain data consisting of 3D line (breaklines, formlines) and point (spot heights, mass points) data. These data were primarily obtained from manual stereo photogrammetric methods using aerial images with a ground sampling distance of ≤ 30 cm, supported with the vectorization of existing topographic maps and field data collection. To eliminate eventual gross errors from the terrain data, a method proposed by Gašparović et al. [53] was used. Finally, DTM PHM in the raster format with a spatial resolution of 0.5 m was created from the 'error-free' digital terrain data using Global Mapper software (ver. 19, Blue Marble Geographics, Hallowell, Maine, USA). A detailed description regarding data characteristics and processing as well as accuracy analyses for both DTM LiD and DTM PHM can be found in the studies of Balenović et al. [54] and Gašparović et al. [53].

Methods
The simplified methodological workflow used in this research is presented in Figure 2.

UAS Image Orientation
Within this research, three different methods of UAS image block orientation were applied and evaluated ( Figure 2): 1. The InSO method based on image tie-points and 5 irregularly distributed GCPs; 2. The GNSS-SO1 method based on tie-points, 5 GCPs and non-PPK single-frequency carrier-phase GNSS data (absolute positioning); 3. The GNSS-SO2 method based on tie-points, 5 GCPs and using PPK dual-frequency carrier-phase GNSS data (relative positioning). The single-frequency GNSS data used for the GNSS-SO1 method were recorded in image metadata and did not require any additional pre-or post-processing. Dual-frequency GNSS data used for the GNSS-SO2 method were recorded in RINEX file, which further enabled the processing and application of relative GNSS positioning. Prior to image orientation, the PPK procedure was carried out in the Trimble Business Center v4.0 software. VRS, generated by the CROPOS Geodetic Precise Positioning Service (GPSS), was used for processing. All GNSS positions relating to the exposure moments (camera stations) of both flights were solved with the fixed solution. The differences between the single-and dual-frequency exposure position estimates for the first and second flights are presented in Figure 3, whereas Table 3 shows basic statistics for the above-mentioned differences. Both Figure 3 and Table 3 indicate that the offset of single-frequency measurements in respect to dualfrequency measurements is consistent for the N and h axes, while differences on the E axis show a significantly larger amount of both differences and dispersion. The main reasons for that could be the imaging method, navigation system synchronization, and flight direction. The direction of successive flight lines was from east to west and from west to east, and the ordinal numbers of the last and first images within each flight line correspond to the "jumps" on the Figure 3 E axis differences. Both datasets were collected with the same UAS; however, shutter synchronizations were carried out in a different way. Single-frequency measurements are primarily used for the navigation of the UAS, hence precise synchronization of the camera shutter is not mandatory. Therefore,

UAS Image Orientation
Within this research, three different methods of UAS image block orientation were applied and evaluated ( Figure 2):

1.
The InSO method based on image tie-points and 5 irregularly distributed GCPs; 2.
The GNSS-SO 2 method based on tie-points, 5 GCPs and using PPK dual-frequency carrier-phase GNSS data (relative positioning).
The single-frequency GNSS data used for the GNSS-SO 1 method were recorded in image metadata and did not require any additional pre-or post-processing. Dual-frequency GNSS data used for the GNSS-SO 2 method were recorded in RINEX file, which further enabled the processing and application of relative GNSS positioning. Prior to image orientation, the PPK procedure was carried out in the Trimble Business Center v4.0 software. VRS, generated by the CROPOS Geodetic Precise Positioning Service (GPSS), was used for processing. All GNSS positions relating to the exposure moments (camera stations) of both flights were solved with the fixed solution. The differences between the single-and dual-frequency exposure position estimates for the first and second flights are presented in Figure 3, whereas Table 3 shows basic statistics for the above-mentioned differences. Both Figure 3 and Table 3 indicate that the offset of single-frequency measurements in respect to dual-frequency measurements is consistent for the N and h axes, while differences on the E axis show a significantly larger amount of both differences and dispersion. The main reasons for that could be the imaging method, navigation system synchronization, and flight direction. The direction of successive flight lines was from east to west and from west to east, and the ordinal numbers of the last and first images within each flight line correspond to the "jumps" on the Figure 3 E axis differences. Both datasets were collected with the same UAS; however, shutter synchronizations were carried out in a different way. Single-frequency measurements are primarily used for the navigation of the UAS, hence precise synchronization of the camera shutter is not mandatory. Therefore, methods of synchronization associated with single-frequency measurements usually have a lower accuracy. One of the common methods of synchronization is logging the trigger signal message, but this method only logs the moment when the trigger command is sent to the camera. The exact moment of the sensor exposure is not logged, nor is the delay between the trigger signal and moment of exposure constant. In this specific case, considering that the cruise speed of the UAS used was 24 km·h −1 and the mean absolute error was ≈ 3 m on the E axis (Table 3), the average delay was ≈ 0.12 s. The method of synchronization used with the dual-frequency measurements is based on logging the exact moment of the sensor exposure via flash mount. Flash is very tightly synchronized with the shutter release and connecting it to the GNSS receiver allows us to directly log the sensor exposure moment in the GNSS observation RINEX file. This method of synchronization has a very small delay and it is stable in time [55]. Another possible cause or at least contributor to the detected error behavior is rolling shutter. The camera used employs a rolling shutter, which contributes to errors in the flying direction; however, rolling shutter can be, and was, fitted in the adjustment [56].
Remote Sens. 2020, 12,404 8 of 19 methods of synchronization associated with single-frequency measurements usually have a lower accuracy. One of the common methods of synchronization is logging the trigger signal message, but this method only logs the moment when the trigger command is sent to the camera. The exact moment of the sensor exposure is not logged, nor is the delay between the trigger signal and moment of exposure constant. In this specific case, considering that the cruise speed of the UAS used was 24 km•h −1 and the mean absolute error was ≈ 3 m on the E axis (Table 3), the average delay was ≈ 0.12 s. The method of synchronization used with the dual-frequency measurements is based on logging the exact moment of the sensor exposure via flash mount. Flash is very tightly synchronized with the shutter release and connecting it to the GNSS receiver allows us to directly log the sensor exposure moment in the GNSS observation RINEX file. This method of synchronization has a very small delay and it is stable in time [55]. Another possible cause or at least contributor to the detected error behavior is rolling shutter. The camera used employs a rolling shutter, which contributes to errors in the flying direction; however, rolling shutter can be, and was, fitted in the adjustment [56].   Table 3. Basic statistics of differences between the single-frequency and dual-frequency exposure (camera) position estimates for the first and second flights (E -easting, N -northing, h -altitude).  Photogrammetric processing, i.e., the image orientation for all three methods (InSO, GNSS-SO 1 , GNSS-SO 2 ), was conducted using Agisoft PhotoScan v1.4.3 software [57]. The a priori accuracies of the GCPs' were estimated by the GNSS receiver, and were set to 10 cm and 17 cm in the horizontal and vertical directions, respectively. The inaccuracy of the estimated positions, atypical for the RTK method, is mostly due to low-performance receiver hardware, accompanied with the unfavorable environment (leaf-on condition) for GNSS measurements. The accuracies of the exposure station positions obtained with single-frequency measurements were set to 3.8 m. The a priori accuracy was estimated in reference to the dual-frequency position measurements. Exposure station positions obtained with the dual-frequency measurements had a priori accuracy of 4 and 8 cm in the horizontal and vertical directions, respectively.

UAS Point Cloud Generation
After the image block orientation using three different methods (InSO, GNSS-SO 1 , GNSS-SO 2 ), the Agisoft PhotoScan (v1.4.3) DIM algorithm [57] was applied to generate three different dense point clouds (PCs). Preliminary tests [58] showed that for images with the technical characteristics used in this research, the best results in H L estimation were obtained with the level 1 image pyramid ("High" parameter selection). This means that images were firstly resampled to half of the original resolution on both axes, and thus GSD of 16 cm was used in the DIM procedure. Downsampling the images reduced the number of pixels by a factor of 4, which significantly increased the processing speed of the DIM algorithm. The resulting density of the generated PCs was 72.32 points·m −2 , on average, due to the flight plan, flying height, and image downsampling. This is in agreement with several studies [14,17,41,43,44], which showed that lower point density could be successfully used to estimate forest variables at plot levels.

Extraction and Calculation of Point Cloud (PC) Metrics
The PCs were normalized with the DTM combined from DTM LiD and improved DTM PHM . The normalized PCs were then clipped to the extent of the plots area, and various metrics were extracted for each plot using the FUSION LDV v3.80 open source software [59]. In accordance with other similar studies [3,5,7,43], the minimum height threshold of 2 m was applied to remove ground and understory vegetation (shrubs and trees with dbh≤10 cm). Furthermore, the height thresholds of 5, 10, 15, 20, and 25 m were applied to calculate the additional density (canopy) cover metrics. A total of 47 metrics were extracted and grouped into height metrics, height variability metrics, height percentiles, and canopy cover metrics (Table 4). For more details on extracted metrics, please refer to McGaughey [59]. Table 4. Metrics extracted from UAS point clouds (PCs), calculated for each plot and used as potential independent variables in statistical analysis (generation and validation of Lorey's mean height models).

Variable Group Variable (Abbreviation and Description)
Height metrics h min -minimum height; h max -maximum height; h mean -mean height; h mode -mode height Height variability metrics SD-standard deviation; VAR-variance; CV-coefficient of variation; IQ-interquartile distance; Skew-skewness; Kurt-kurtosis; AAD-average absolute deviation; MAD med -Median of the absolute deviations from the overall median; MAD mode -Median of the absolute deviations from the overall mode; CRR-Canopy relief ratio ((mean − min)/(max − min)); SQRT meanSQ -Generalized mean for the 2nd power (Elevation quadratic mean); CURT meanCUBE -Generalized mean for the 3nd power (Elevation cubic mean); L 1 , L 2 , L 3 , L 4 -L moments; L CV -moment coefficient of variation; L skew -moment skewness; L kurt -moment kurtosis

Generation and Validation of Lorey's Mean Height (H L ) Models
Out of a total of 99 circular sample plots, 66 plots were selected for the models' generation and 33 plots for the models' validation (HV validation approach). The plots were listed according to age, from youngest to oldest, and every third plot was selected in the validation dataset.
All extracted plot-level PC metrics were considered in statistical modelling (multivariate linear regression) as potential explanatory (independent) variables, while H L calculated from field ground-truth data (Equation (1)) was used as a response variable. Prior to the modelling, the large number of potential explanatory variables was reduced by applying a two-step pre-selection approach [43]. Firstly, the number of potential variables was reduced based on the Pearson correlation coefficient (r). All variables that had a low correlation (r < ±0.5) with the observed H L were excluded. The remaining variables were then included in the second step-the in-group collinearity analysis. Within each group separately, the correlation (r) between variables was calculated. Variables that showed high within-group collinearity (r ≥ ±0.7) and those that exhibited lower correlation with H L were excluded one by one, while the remaining variables (at least one per group) entered the multivariate linear regression. By applying a backward stepwise approach, the-best-fit H L model was selected for each PC, i.e., for PCs produced from UAS images, whose orientation was generated using the three proposed methods (InSO, GNSS-SO 1 , GNSS-SO 2 ).
The models were validated to evaluate the accuracy of each H L model by performing:

1.
The validation over the independent validation dataset of 33 plots (HV), which was not used to derive the models. The H L estimates from the developed models were compared with corresponding field data and evaluated by means of adjusted coefficients of determination (R 2 adj ), mean error (ME) (Equation (2)), relative mean error (ME%) (Equation (3)), root mean square error (RMSE) (Equation (4)), and relative root mean square error (RMSE%) (Equation (5)): Remote Sens. 2020, 12, 404 where H L i is the predicted (UAS estimated) Lorey's mean height of plot i, H L is the observed (from field data) Lorey's mean height of plot i, n is the number of plots, and H L is the mean of the observed values.

2.
The LOOCV statistical method [46,47], based on 66 sample plots, was used for the model's development. LOOCV is the iterative procedure of n iterations, where n is the number of all measurements (field plots). In n iterations (n = 66), one measurement was removed from the dataset and the selected model wasfitted using the remaining n-1 measurements. The model was then validated using the removed measurement. After the process of n-1 iterations was done, the model accuracy was estimated by averaging validation results (residuals) from all iterations.
A two-step pre-selection procedure, the backward stepwise regression, and the independent dataset validation were performed using STATISTICA 11 software [60], while MATLAB R2018a and freeware MATLAB RSR toolbox [61] were used to perform LOOCV.

Results
After generating the UAS images' orientation with three methods (InSO, GNSS-SO 1 , GNSS-SO 2 ), dense PCs were generated for each image orientation outcome. Vertical differences between InSo and GNSS-SO 2 PCs, as well as between GNSS-SO 1 and GNSS-SO 2 PCs for two selected (exemplary) plots, are presented in Figure 4. The GNSS-SO 2 PC was used as a reference since it was expected to have the strongest geometry of the image orientation in comparison to InSO and GNSS-SO 1 based PCs. Out of two selected exemplary plots, one plot (marginal plot) was located far from the GCPs (Figure 4a,b), while the second one (central plot) was located in the vicinity of the GCPs (Figure 4c,d). As expected, considerably larger vertical differences were observed for the marginal plot, due to the higher distance of the marginal plot from GCPs as well as to the irregular distribution of GCPs (Figure 1), which leads to vertical distortion after the image block orientation methods are applied, especially in the areas without GCPs (Figure 4a,b). Furthermore, for both marginal and central plots, larger vertical differences can be observed between InSO and GNSS-SO 2 PCs (Figure 4a,c) than between GNSS-SO 1 and GNSS-SO 2 PCs (Figure 4b,d). This confirms that despite their lower accuracy, the single-frequency GNSS measurements contribute to the block stability and accuracy of image orientation. Table 5. Plot-level H L models (equations with parameters and included predictors) developed using multivariate linear regression (backward stepwise) for each point cloud (orientation method).

Orientation Method Model
The plot-level metrics were extracted and calculated, as shown in Table 4, and the same statistical procedure was applied to generate H L estimation models (Table 5). For each PC, the best-fit plot-level H L linear models were developed based on the backward stepwise regression of pre-selected variables (plot-level PC metrics; independent variables) and field H L (dependent variable) in the selected 66 sample plots (Table 5). All models and their parameters show high statistical significance (p < 0.01). Both validation approaches (HV and LOOCV) show the consistent trend of models' performance (Table 6, Figure 5). The best performing is the GNSS-SO2-based model, followed by the GNSS-SO1and InSO-based models.
For HV (33 independent validation plots), the models exhibit just slightly lower R 2 adj values and slightly higher RMSE% and ME% values than for the modeling dataset (66 plots). The R 2 adj values between the modeling dataset and HV differ by 0.022, 0.028, and 0.021 for InSO-, GNSS-SO1-, and GNSS-SO2-based models, respectively. The differences between the RMSE% values for two datasets range from 1.063% to 1.295%, whereas differences between the ME% values range from 0.259% to 0.643%. A good agreement, i.e., a small difference between the validation results for the modeling and HV datasets, indicates that the derived models provide reliable HL estimates in the study area of lowland pedunculate oak forests.
As expected, LOOCV produces better validation results than HV, of which the results are very similar to the results of the modeling dataset. The R 2 adj values between the modeling dataset and The modeling results obtained from the three models proposed in Table 5 show small differences in their performance (Table 6, Figure 5). The differences between models are less than 0.02% and 0.6% for R 2 adj and RMSE%, respectively, whereas the ME% differences between all models are less than 0.001%. Furthermore, the high R 2 adj values (>0.943), low RMSE% (<4.488%), and ME% (<0.001%) values indicate a good model fit for all three models. However, the models differ by number and type of included independent variables (predictors) ( Table 5). The GNSS-SO 2 based model, which includes two predictors, exhibits the best fit to data when compared to InSO-and GNSS-SO 1 -based models, which include three and four predictors, respectively. Table 6. Results of multivariate linear regression (modeling dataset) and validation results (HV and LOOCV) for selected plot-level H L models described in Table 5. LOOCV differ by only 0.006, 0.004, and 0.002 for InSO-, GNSS-SO1, and GNSS-SO2-based models, respectively. The differences between the RMSE% values range from 0.169% to 0.356%, whereas the differences between the ME% values range from 0.004% to 0.056%. Table 6. Results of multivariate linear regression (modeling dataset) and validation results (HV and LOOCV) for selected plot-level HL models described in Table 5.  Both validation approaches (HV and LOOCV) show the consistent trend of models' performance (Table 6, Figure 5). The best performing is the GNSS-SO 2 -based model, followed by the GNSS-SO 1and InSO-based models.

Dataset
For HV (33 independent validation plots), the models exhibit just slightly lower R 2 adj values and slightly higher RMSE% and ME% values than for the modeling dataset (66 plots). The R 2 adj values between the modeling dataset and HV differ by 0.022, 0.028, and 0.021 for InSO-, GNSS-SO 1 -, and GNSS-SO 2 -based models, respectively. The differences between the RMSE% values for two datasets range from 1.063% to 1.295%, whereas differences between the ME% values range from 0.259% to 0.643%. A good agreement, i.e., a small difference between the validation results for the modeling and HV datasets, indicates that the derived models provide reliable H L estimates in the study area of lowland pedunculate oak forests.
As expected, LOOCV produces better validation results than HV, of which the results are very similar to the results of the modeling dataset. The R 2 adj values between the modeling dataset and LOOCV differ by only 0.006, 0.004, and 0.002 for InSO-, GNSS-SO 1 , and GNSS-SO 2 -based models, respectively. The differences between the RMSE% values range from 0.169% to 0.356%, whereas the differences between the ME% values range from 0.004% to 0.056%.

Discussion
To date, a number of recent studies have confirmed the great potential of UAS photogrammetry in forest inventory, at both individual tree- [16,18] and plot-levels [14,[40][41][42][43][44][45]. In general, the quality of the obtained forest inventory data depends on the quality and accuracy of the UAS products (e.g., PCs, DSMs) and DTM, while the quality of the UAS products depends on the image block orientation method applied to UAV data, among other factors such as flight characteristics and target properties [62].
The innovative component of the present study is that it focuses on the impact of widely-used image block orientation methods on not just UAS products but rather on the accuracy of forest attribute estimates. By adding this additional step to explore the performance of these methods on plot-level mean tree height (H L ), we explore their impact on the final remote sensing product directly used in forest inventory.
Overall, the best performing of the models is the GNSS-SO 2 -based model (RMSE%: HV = 5.18%, LOOCV = 4.06%), followed by the GNSS-SO 1 -based model (RMSE%: HV = 5.34%, LOOCV = 4.37%), while the worst performing is the InSO-based model (RMSE%: HV = 5.55%, LOOCV = 4.84%). In other words, the GNSS-SO 2 orientation method introduces the least errors in the estimation of the plot-level mean tree height (H L ). The GNSS-SO 1 method performs better than the InSO orientation method alone. The order of performance is as expected; however, the negligible differences in the performances of the regression models suggest that the image block orientation methods used in this study have no considerable effect on the estimation of forest attribute H L , which is somewhat surprising. This indicates that despite the low number of GCPs and their irregular distribution throughout the study area, the GCPs, as well as the high overlap (90%, 80%) of UAS images, contribute considerably to the block stability and accuracy of image orientation in the InSO method, and consequently have a positive impact on the geometric accuracy of the PCs. Despite their lower accuracy, the addition of single-frequency GNSS measurements to tie-points and GCPs slightly improves the accuracy of image orientation ( Figure 4) and H L estimates (Table 6, Figure 5). Additional slight improvements in image orientation and H L estimation are achieved by replacing single-frequency GNSS measurements with dual-frequency GNSS measurements and by incorporating PPK into the GNSS-SO 2 method (Table 6, Figure 5). It can be assumed that homogeneous forest structure and flat terrain in the study area also contribute to the small differences in regression models' performances. In areas with more heterogeneous forest structure and more complex terrain, greater differences among the results obtained with different image orientation methods could be expected. To achieve greater block stability and geometric accuracy of PCs, especially when applying the InSO method, additional cross-flights at different flying heights should be considered, particularly in more complex forest areas. Nevertheless, when using GNSS-SO (PPK) the observed EO positions contribute greatly to the block stability and parameter decorrelation, suggesting that the trade-off between the cost of additional cross-flights and accuracy improvements in forest variable estimation is justified.
In general, the results of this study are in agreement with the findings of recent studies [14,[40][41][42][43][44][45]] which reported that plot-level tree heights (e.g., dominant height; Lorey's height) can be estimated with high accuracy. Namely, the RMSE values for Lorey's mean height ranged from 4.24% to 17.87% [14,[40][41][42][43][44][45]. However, the estimates of Lorey's mean height in this study have attained higher accuracy than in most previous studies [14,[40][41][42][43][44]. Slightly higher accuracy (RMSE = 4.24%) was reported only by Shen et al. [45] compared to the InSO and GNSS-SO 1 results of this study, while the GNSS-SO 2 results of this study are of higher accuracy (RMSE = 4.06% for LOOCV) than reported in Shen et al. [45]. When applying the InSO orientation method, Shen et al. [45] used 30 GCPs for 45 plots. Within this study, only five GCPs were used to cover an area with 99 plots. Furthermore, the research of Shen et al. [45] was conducted in a ginkgo plantation, which is a considerably simpler forest type in terms of species composition and tree height variations than the mixed pedunculate oak forest. In addition to PC metrics derived from RGB images, Shen et al. [45] used spectral metrics derived from multispectral images to fit regression models and estimate H L . It should be emphasized that almost all previous studies [14,[40][41][42]44,45] applied the InSO approach. The only exception is the study of Balenović et al. [43], which preceded this research, where the GNSS-SO 2 method was applied on a similar dataset, but with a different research aim. Any direct comparison with other studies is difficult, as reported accuracies could depend on many other factors such as forest structure and site characteristics, UAS and camera specifications, flight conditions, photogrammetric software, type of used predictors, and modeling and validation strategies.
The present study considers the performance of two validation methods, HV and LOOCV. Ground sampling design, as part of every field-based validation process, plays an important role in remote sensing studies. It must include a careful selection, a sufficient number of points, and a proper sampling strategy [63]. The consistent results between the two validation approaches (HV and LOOCV) and modeling dataset suggest a minimal potential negative effect that may occur due to the sampling process, which is common for homogenous forests. By directly comparing the estimated H L information with field measurements over the independent dataset (HV), not used in the modeling process, the reliability of the models has been demonstrated. To avoid any potential uncertainties due to the number of sampling points, the LOCCV statistical approach additionally confirms the strong performance of the regression models, and ultimately reveals the performance of the image block orientation methods and their impact on the estimated H L information. The consistent results among the validation methods support the conclusion regarding the structural homogeneity of the study site, with no extreme values. Although the results are similar, it seems that the HV approach provides more realistic information than the LOOCV statistical approach. There are pros and cons for both methods, but the HV approach is preferable if there are enough measurements to be divided into modeling and validation subsets [64]. It is difficult to provide a general rule on the amount of data required for both modeling and validation, but it depends greatly on the complexity of the observed object [64,65] structure. However, as the collection of large amounts of field reference data for supporting remote sensing-based forest inventory is labor-and time-consuming, the LOOCV approach is a valid alternative. Thus, it is not a surprise that out of seven recently published studies, which explore the estimation of H L using UAS photogrammetry, six of them employ the LOOCV approach [14,[40][41][42][43][44], and one of them uses the 10-fold CV method [45]. No direct field-based validation on an independent dataset (HV approach) was used in the previous studies, which makes the present study more systematic.

Conclusions
In this study, the effect of different orientation methods (InSO, GNSS-SO 1 , GNSS-SO 2 ) applied to UAS images and PCs on the accuracy of plot-level H L estimates was tested on the example of mixed, even-aged pedunculate oak forests characterized by flat terrain. The results of this study confirmed the great potential of UAS photogrammetry in forest inventory, even when low-cost (consumer-grade) UASs were used. Such consumer-grade UASs are usually equipped with amateur cameras and low-quality GNSS sensors that record single-frequency GNSS data of lower accuracy. Thus, image orientation is usually performed with tie-points obtained from SfM and field-measured GCPs solely (InSO method) or with the addition of single-frequency GNSS measurements (GNSS-SO 1 method). Despite its lower accuracy, this study showed that the addition of single-frequency GNSS measurements to tie-points and GCPs was useful because it improved the accuracy of H L estimates (GNSS-SO 1 method). The highest accuracy of H L estimates was achieved by the GNSS-SO 2 method, i.e., by adding PPK dual-frequency GNSS measurements to tie-points and GCPs. Unlike single-frequency GNSS measurements obtainable by consumer-grade UASs, high-accuracy GNSS measurements are obtainable only with professional and more expensive UASs, which certainly limits their availability. Small differences among the results obtained by different image block orientation methods indicated that, despite the low number and non-optimal spatial distribution of GCPs, they contributed considerably to the block stability and accuracy of image orientation. This is especially important for the InSO method, which relies solely on GCPs for georeferencing, as well as for the GNSS-SO method when using low-accuracy GNSS data. It can be assumed that the more or less homogeneous forest structure and flat terrain in the study area also contributed to small differences among the results obtained by different orientation methods. Consequently, greater differences could be expected in more complex forest areas in terms of forest structure and terrain characteristics. The findings of this research should serve as a basis for further studies that should include different forest types and other important forest variables such as diameter at breast height, basal area, volume, and biomass.