Combination of UAV Photogrammetry and Field Inventories Enables Description of Height–Diameter Relationship within Semi-Arid Silvopastoral Systems

Abstract

Pollarding oak trees is a traditional silvopastoral technique practiced across wide areas of the northern Zagros mountains, a unique and vast semi-arid forest area with a strong cultural and ecological significance. So far, the effects of pollarding on tree structure in terms of DBH (diameter at breast height)~H (height) relationships within the typical pollarding cycle, which often lasts 4 years, has not been scientifically described. Here, we combine field inventories of DBH with H obtained from photogrammetric UAV flights for the first time to assess DBH~H relationships within this system. We conducted the research at six pollarded forest sites throughout the Northern Zagros. The sampling encompassed all three main species of coppice oak trees. In the case of multi-stem trees, we used the maximum DBH of each tree that formed a unique crown. A linear relationship between UAV and extracted H and the maximum DBH of pollarded trees explained a notable part of the variation in maximum DBH (R² = 0.56), and more complex and well-known nonlinear allometries were also evaluated, for which the accuracies were in the same range as the linear model. This relationship proved to be stable across oak species, and the pollarding stage had a notable effect on the DBH~H relationship. This finding is relevant for future attempts to inventory biomass using remote sensing approaches across larger areas in northern Zagros, as well as for general DBH estimations within stands dominated by pollarded, multi-stem coppice structures.

1. Introduction

Forests cover around 30% of the Earth’s land area [1], including semi-arid areas. Iran is climatically diverse and encompasses a considerable portion of the world’s climate zones [2]. However, the country is known to be dominated by semi-arid, arid, and fragile vegetation ecosystems [3]. The Zagros mountains are home to a significant portion of Iran’s semi-arid forests and woodlands, which encompass > 40% of Iran’s forests, accounting for ca. six million ha [4,5]. These forests are particularly important from the socioeconomic, soil conservation, and water-quality perspectives, as well as for providing considerable non-market ecosystem services [6]. However, their quantity and quality have been almost constantly deteriorating over the last few decades due to a variety of natural and human-induced factors [7,8,9,10,11], including, in particular, the high level of dependency of local livelihood on these forests [12,13].

This close connection between people and nature in northern Zagros (located in the West Azerbaijan and Kurdistan provinces, Iran) [14] has led to significant alterations in the vegetative composition [15] and the primary structure of both trees and forest stands [16]. Historically, the high dependence of local forest dwellers on natural resources has made inventive utilization of the forest to supply both human and livestock inevitable [17]. Most forest dwellers practice some sort of silvopastoralism and highly depend on oak forests to produce forage for their livestock throughout the year [18]. In particular, traditional animal husbandry in northern Zagros faces a lack of pastures and fodder in the cold seasons. To adapt their livelihoods to such limitations, the locals have developed a traditional silvopastoral method called pollarding (Kurdish: Galazani) [6,19]. Pollarding is a generic name for pruning tree branches to provide both fodder (for domestic animals) and fuelwood [20]. Pollarding has been reported to be practiced in parts of Europe, Asia, Africa, and America, especially for the benefit of using leaves and branches as animal fodder. Pollarding is hence an important factor in livestock management and relates to this important source of income in the farming economy [21,22,23,24,25]. Oak pollards are a part of traditional forest management in many ecosystems and have a very long history of about 4000 years [26]. In Iran, and as a part of the traditional management system in the northern Zagros, each rural family traditionally owns, i.e., tenures, a part of the forest (Kurdish: Galajar) [27]. Each Galajar is divided into three or four parcels with equal fodder productivity, termed Shanegala (in Kurdish) [28]. This is due to the fact that a given Shanegala is used to perform pollarding for winter fodder each year. Each Shanegala is pollarded in a three or four year pollarding cycle [4,28], which annually occurs from mid-September to early October (pre-fall). However, pollarding is not carried out in this way everywhere and is carried out irregularly [29]. During pollarding, forest dwellers cut leafy branches (Kurdish: Bakhe) of oak trees, mostly Quercus libani Olivier and Q. infectoria Olivier., and store them on large trees (Kurdish: Dargala), on the ground, or on rock fragments in a cone-shaped formation known as a Gala or a Loya Gala (in Kurdish). In winter, dried leafy Loya Gala branches are used to feed livestock (especially Markhoz goats and sheep) [16]. The locals estimate pollarding ages to conduct the above-mentioned pollarding cycle. Whereas a freshly pollarded tree is called Korpa, the tree is termed Kor, Khert, and 2 Khert following one, two, and three years from the time of its pollarding, respectively [30] (Figure 1).

One of the initial consequences of pollarding is reducing the tree canopy, which causes stress on the tree and changes the typical pattern of tree growth [31,32]. However, some studies [6,16,32,33,34] have shown that pollarding may result in a decrease in height, trunk, crown area, canopy vigor, health, seed and coppice regeneration, and density.

However, scientific knowledge on this traditional silvicultural approach is very sparse, and the exact effects of pollarding on tree structure and vitality are not yet fully understood (an overview of earlier studies focusing on this topic is summarized in Table 1). The first step to thoroughly examining the pollarding approach and its temporal trajectories is to adequately capture the tree structure of pollarded trees [35,36]. Field inventories typically involve measuring the diameter at breast height (DBH), tree height, and the tree crown diameter [37,38] but without considering the intrinsic features of pollarded trees, like their distorted DBH~height and DBH~crown area relationships. With respect to potential monitoring approaches, timely field measurements are challenging, as pollarding typically takes place on mountain slopes [16,30,39], a problem that is exacerbated not only by the subjectivity and required logistics for conducting field measurements (particularly for attributes like tree height) but also by the fact that measurements might be prevented or opposed by the local dwellers that assert a customary right to forest tenure over Galajars.

As an alternative or complement, remote sensing- and photogrammetry-based approaches may partly replace or amend field measurements [40]. Remote sensing has been historically used in a variety of fields related to vegetation monitoring [41,42]. Well-known pattern recognition algorithms in combination with space-borne active [43,44,45] and passive [46,47,48] data have been largely employed yet face difficulties in small-scale applications. Estimating the structural variables of single trees can be accomplished using very-high-resolution (VHR) space-borne data [49], which are, however, often prohibitive due to their high price and the lack of accessibility for politically sensitive areas like Iran. As an alternative, three-dimensional ground sensors (e.g., terrestrial laser scanners) [50] or aerial sensors (e.g., aerial laser scanners) [51] can be used. However, these tools are often either unavailable, costly, or practically infeasible for large-area assessments or subsequent spatial upscaling to larger domains [52]. Consumer-grade unmanned aerial vehicles (UAVs) may offer an intermediate solution [53] at an affordable cost [54]. Compared with very-high-resolution spaceborne data, their application is associated with fewer limitations regarding atmospheric correction and cloud cover because of their relatively low flight altitudes [55]. Additionally, 3D point clouds can be created from raw UAV data by using stereo image-matching techniques [56]. This enhances the level of structural analysis from a planar space to a vertical space [57] and allows tree heights to be obtain with reasonable accuracy, particularly in open forest environments. UAV photogrammetry is not limited to forests when it comes to tree attribute estimation. Several studies have successfully employed UAV photogrammetry to estimate the height and diameter of orchard trees as well [58,59,60].

The overall workflow of UAV photogrammetry in structural tree attribute estimation involves mission planning, acquisition of aerial images using a UAV, preprocessing of the images to correct distortions, image matching to determine position and orientation, generation of a point cloud and digital surface model, subtraction of the digital surface model to create a canopy height model, and estimation of tree heights using the canopy height model [53,61]. The accuracy of the estimated tree heights depends on the quality of the images, the image-matching accuracy, and the effectiveness of segmentation and measurement algorithms.

Despite all of these benefits, it is important to also mention the limitations of UAV data. The very small coverage of most UAVs is the primary obstacle, which is a consequence of several factors, including (1) limited battery power/flight time, (2) limited flight area that can be used for stereo analysis, and (3) the significant data volumes needed for creating digital orthomosaics and digital terrain/surface models [62,63]. Despite these limitations, a few earlier studies applied UAVs in the Zagros region to assess various forest attributes (Table 1).

Table 1. The results of a systematic search on Web of Science^® on relevant UAV-based studies within Zagros. AB in the “Search Query” column stands for abstract.

Search Query	Result
	Reference	Year	Description	Region
AB = (UAV OR UAS) AND AB = (Zagros)	[61]	2022	Edge detection-based method for delineating overlapped crowns of coppice trees using UAV photogrammetry to monitor tree decline	Middle and South Zagros
	[53]	2021	Structural variable extraction of Zagros single trees	Middle and South Zagros
	[64]	2021	Intra- and interspecific interaction investigation and their contribution to growth inhibition in the patches of Pistacia trees and Amygdalus shrubs	South Zagros
	[65]	2017	Mapping part of a wild pistachio nature reserve in Zagros open woodlands	South Zagros
	[66]	2017	Challenges and quality of landslide assessment based on remotely sensed data	Zagros Mountains
AB = (pollarded) AND AB = (Zagros)	[25]	2018	Quantitative description of the pollarding process in the northern Zagros, with a focus on foliage stacks resulting from pollarding	North Zagros
	[67]	2018	Structural analysis of pollarded trees to evaluate the effect of aspect on biometric indices of oak trees	North Zagros
	[16]	2017	Investigation of the effect of pollarding on the increment in diameter of Lebanon oak trees by comparing them to less disturbed stands	North Zagros
	[6]	2014	Investigating the tree-pollarding process in Galajars	North Zagros
	[18]	2012	Allometric relationship introduction for estimating the crown and leaf biomass of Q. libani using the DBH, tree height, crown length, and crown width	North Zagros
	[27]	2010	Evaluation of the capability of IRS-P6 data to separate the pollarding areas	North Zagros

Table 1. The results of a systematic search on Web of Science^® on relevant UAV-based studies within Zagros. AB in the “Search Query” column stands for abstract.

Search Query	Result
	Reference	Year	Description	Region
AB = (UAV OR UAS) AND AB = (Zagros)	[61]	2022	Edge detection-based method for delineating overlapped crowns of coppice trees using UAV photogrammetry to monitor tree decline	Middle and South Zagros
	[53]	2021	Structural variable extraction of Zagros single trees	Middle and South Zagros
	[64]	2021	Intra- and interspecific interaction investigation and their contribution to growth inhibition in the patches of Pistacia trees and Amygdalus shrubs	South Zagros
	[65]	2017	Mapping part of a wild pistachio nature reserve in Zagros open woodlands	South Zagros
	[66]	2017	Challenges and quality of landslide assessment based on remotely sensed data	Zagros Mountains
AB = (pollarded) AND AB = (Zagros)	[25]	2018	Quantitative description of the pollarding process in the northern Zagros, with a focus on foliage stacks resulting from pollarding	North Zagros
	[67]	2018	Structural analysis of pollarded trees to evaluate the effect of aspect on biometric indices of oak trees	North Zagros
	[16]	2017	Investigation of the effect of pollarding on the increment in diameter of Lebanon oak trees by comparing them to less disturbed stands	North Zagros
	[6]	2014	Investigating the tree-pollarding process in Galajars	North Zagros
	[18]	2012	Allometric relationship introduction for estimating the crown and leaf biomass of Q. libani using the DBH, tree height, crown length, and crown width	North Zagros
	[27]	2010	Evaluation of the capability of IRS-P6 data to separate the pollarding areas	North Zagros

However, the utility of UAVs as practical and affordable tools for structural assessments of pollarded trees has not yet been examined and was set as the main objective of this research. We were particularly interested in exploring the relationship between UAV-derived tree heights and DBHs of the trees. This is particularly interesting for potential future use cases to conduct simplified forest inventories exclusively based on UAV flights to estimate wood volume and biomass.

Therefore, the specific objectives of this research are:

• To test the feasibility of consumer-grade UAV data to estimate the height of pollarded trees and to use these heights to estimate the DBH of the trees;
• To understand how the essential DBH~H relationship is affected by the oak species and the stage of pollarding and whether there is a difference between multi- and mono-stemmed trees.

2. Materials and Methods

2.1. Subsection

The Zagros forests comprise 20% of Iran’s total vegetation cover and are a habitat to multiple and endemic woody species, in particular, three oak species, including Q. brantii Lindl., Q. infectoria Olivier, and Q. libani Olivier in both coppice and standard forms. The forests play a critical role in the water supply and the economic and social welfare of the local population [4]. The Zagros forests are biogeographically divided into northern, central, and southern zones [68], with all three oak species jointly occurring only in the northern part [4,68].

We analyzed the structural attributes of pollarded stands across 6 sites near Baneh county in Kurdistan province (Figure 2), one of the richest and most structurally diverse ecosystems in Zagros [18]. We considered the presence of Galajars of various pollarding ages when selecting the sites.

2.2. Materials

UAV photogrammetric products were applied to derive tree heights [37,38]. The derived heights were then used to estimate the field-measured DBH of the trees. Furthermore, the link between height and DBH was examined with respect to tree age, location, and species groups. Tree crown area was originally considered an additional variable [69], but we considered it to be less relevant since it is subject to continuous change during the pollarding process [16,70].

2.2.1. Field Measurements

The field-collected data included a general plot description (general morphology, aspect, presence of regeneration, and average slope), tree species, DBH (full callipering), and pollarding age (

Table 2. Summary of field-collected information.

Site	Area (m2)	Age	Average Slope (%)	Frequency of Tree Species
				Q. brantii	Q. libani	Q. infectoria
Site 1	6360	Korpa	37	14	52	165
Site 2	2850	Kor	30	0	122	81
Site 3	3050	Korpa	10	2	40	17
Site 4	5360	2 Khert	10	5	130	25
Site 5	4920	2 Khert	12	31	36	29
Site 6	1400	Korpa	25	19	22	6

). In addition, we selected 4 reference trees in each site with uniform distribution throughout the site and measured their height using a clinometer. Finally, we checked whether the trees were growing from seeds or from coppice. The precise UTM coordinates of each tree were recorded using a Raymand^® iRo-3 multi-frequency differential GNSS (Makanpardaz Raymand Inc., Tehran, Iran). The inventories were completed at the 6 study sites during July–September 2021, the closest possible date to the UAV flights (June 2021) and well before the start of the annual pollarding season (Section 2.2.2).

The fieldwork revealed that none of the trees in the 6 sites were of seed origin and were therefore all considered coppice trees. Many trees had multiple stems, for which we measured multiple DBH values. Figure 3 represents the structural variations to be expected in our site based on the field observations and according to an experienced local ecologist. We did not examine whether the root networks of individual trees were connected in the field.

The maximum DBH value within an individual coppice tree was considered the DBH to be estimated for each tree [53]. A visual representation of the DBH distribution at each study site can be seen in Figure S1, while

Table 3. Quantitative description of maximum DBHs measured at each site. N: number of trees, SD: standard deviation, Q1/Q3: 1st and 3rd quartiles.

Variable	Site	N	Mean	SD	Minimum	Q1	Median	Q3	Maximum
DBH	Site 1	231	20.0	4.5	11	17	20	23	33
	Site 2	203	17.2	3.5	10	15	17	19	28
	Site 3	59	18.9	4.6	10	16	18	21	34
	Site 4	160	18.7	4.4	11	15	18	21	34
	Site 5	96	13.6	5.1	5	9	13	17	27
	Site 6	47	15.8	5.1	5	13	15	19	29

shows the quantitative summary of DBH measurements.

2.2.2. UAV Imaging

We used a consumer-grade DJI Phantom 4 pro multi-rotor UAV (DJI, 2016) for aerial imaging. The device comprised a three-axis stabilization gimbal, a 1” CMOS sensor camera, and an 8.8 mm/24 mm lens with a field of view of 84° (https://www.dji.com/phantom-4-pro/info/, accessed on 14 August 2023). The flights for different sites were designed according to site-specific topographic conditions and tree cover and were designed with an iOS version of Pix4DCapture installed on an iPad 2018 tablet (accessible at https://support.pix4d.com/hc/en-us/articles/204010419-iOS-Pix4Dcapture-Manual-and-Settings, accessed on 14 August 2023). The flights were conducted for all six sites between 14 and 16 June 2021 (

Table 4. UAV imaging technical specifications.

Group	Parameter	Value
Flight	Altitude (m)	80
	Overlap (%)	80
	Side overlap (%)	80
	Speed	Medium
	Aspect	Nadir, oblique (60°)
Georeferencing	3D control points	5 accurate 3D points were measured using PPK GNSS for each site
Images	Resolution (px)	5472 × 3648
Products	Models	DSM, DTM, CHM, orthomosaic
	Resolution (cm)	5

). Photogrammetric 3D modelling in forests is generally challenging [71]. To partially overcome the challenge, we additionally designed an oblique cross-flight in addition to the nadir cross-flights based on our previous experiences as well as suggestions by former studies [72,73]. The choice of a 60° oblique flight angle was based on preliminary tests and following the recommendations of earlier studies [56]. The nadir and oblique flights were conducted one minute apart.

Two approaches for processing the photogrammetric data were tested: (1) dense point cloud production for oblique and nadir flights as two separate projects, followed by merging the two models, and (2) merging the images of both flights and producing a 3D model using all of the images. Our visual and quantitative checks on the results of both scenarios showed that the first scenario added considerable noise during the merging of the two models despite being more efficient from a processing cost point of view. The second approach was more time-consuming due to the simultaneous processing of more images but resulted in denser and more accurate models. Therefore, we applied the second option for all sites.

2.3. Methodology

The applied methodology primarily involved extracting tree heights from the crown height model (CHM) [74], followed by examining DBH~H relationships using ordinary least squares (OLS) models as well as a set of common, more complex, nonlinear approaches tested and suggested by the literature (see Section 2.3.2).

2.3.1. Tree Height Estimation

Production of 3D models from overlapping UAV imagery [75] is highly dependent on the surface color texture [63,71,76]. Due to the significant loss of leaved tree canopy within pollarded sites, trees exhibit a relatively poor texture, which results in technical obstacles during the image-matching process. Here, we produced optimal models by modifying the key parameters of 3D model production. We performed an exhaustive search for the optimal Agisoft Metashape 1.8 (https://www.agisoft.com, accessed on 14 August 2023) hyperparameters [77]. Then, we automatically filtered the key points according to a strategy that was based on an iterative process as follows: [78,79].

• Removing all points that are seen in three images or fewer;
• Removing the key points in such a way that the reprojection error is halved, followed by optimizing the camera parameters;
• Removing the points in such a way that the reconstruction uncertainty is halved, followed by optimizing the camera parameters;
• Removing points in such a way that the projection accuracy is halved, followed by optimizing the camera parameters;
• Returning to step 2 and repeating the steps until the stopping condition is reached.

It should be noted that our stop condition was reprojection error = 0.3 px, reconstruction uncertainty = 5 px, and projection accuracy = 3 px. Then, we created the digital surface model (DSM), filtered the trees, and interpolated the area of the tree locations to create the digital terrain model (DTM) [80] and created the CHM by subtracting the DTM from the DSM. The tree heights could then be extracted from the CHM [81]. To do so, the trees must first be segmented. Since the trees were generally sparse, we segmented the tree crowns using the straightforward and efficient marker-controlled watershed segmentation (MCWS) [82,83]. Using the maximum operator to extract the tree heights was infeasible due to the presence of noise in the UAV-based canopy height models [84]. To avoid a strong influence of corresponding outlier pixels, we assumed a normal distribution of the CHM values in each tree segment and identified the largest value as the tree height at a 95% confidence level. Thus, we considered heights that were >2SD (2 × standard deviation, the 95th percentile of all pixels in the tree crown segment) from the average. Figure S2 and

Table 5. Quantitative description of tree heights extracted from the CHM. H: height.

Variable	Site	N	Mean	StDev	Minimum	Q1	Median	Q3	Maximum
H	Site 1	231	6.7	1.2	3.4	5.7	6.5	7.6	9.6
	Site 2	203	4.8	0.9	2.4	4.1	4.9	5.3	7.5
	Site 3	59	3.7	0.5	2.9	3.4	3.5	3.7	5.7
	Site 4	160	5.5	1.2	2.2	4.5	5.6	6.3	9.7
	Site 5	96	4.5	1.2	1.7	3.8	4.0	5.5	6.9
	Site 6	47	3.5	0.7	2.1	3.1	3.4	3.6	5.7

show the visual and quantitative description of the tree height extraction results. To validate the UAV-derived tree heights, we compared them with the field-measured heights. The agreement with field-measured heights was reasonable, with an average difference of 0.3 m to 1.2 m (Figure S3).

2.3.2. DBH Estimation

Despite the inherent correlation of DBH and tree height [85,86,87], pollarded trees are not comparable to undisturbed (unpollarded) oak trees due to significant changes in the crown area during the pollarding cycle (Figure 1). Thus, the DBH~height relationship might also be influenced by other variables. We hypothesized that the link between tree height and DBH differs across (1) single- and multi-stemmed trees, (2) trees of different pollarding ages, (3) the number of trunks, (4) tree species, and (5) their interactions. Items 1, 3, 4, and 5 necessitate direct measurement of all individual trees and quantitative and qualitative measurements by an ecologist. Conversely, item 2 can be obtained exclusively through experience or by asking the Galajar’s local owner, thus negating the need for an examination of the individual trees. Thus, we examined all of these variables, along with the UAV-derived tree heights, to estimate the DBH to answer the question of whether field inventories for all single trees (items 1, 3, 4, and 5) are needed to estimate the DBH from the height obtained by the UAV. To identify the most important predictors, we applied a backward stepwise selection until the stopping condition was met. Then, the most effective independent variables were included in the final model [88]. Most ideally, we were interested in finding a linear relationship between H and DBH, which is the least dependent on additional field inventories and is easier to understand because of its linearity. However, we were also interested in examining the nonlinear relationship between these two variables. To consider possible nonlinear DBH~H relationships, we also used three well-known nonlinear relationships that were suggested by the relevant literature (

Table 6. List of well-known H~DBH nonlinear allometries used in this study.

Model	Equation	Source
1	$1.3+\mathrm{a}{\left(1-{\mathrm{e}}^{-\mathrm{b}\times \mathrm{DBH}}\right)}^{\mathrm{c}}$	[90]
2	$1.3+\mathrm{a}\left({ \mathrm{e}}^{-\mathrm{b}\times \exp \left(-\mathrm{c}\times \mathrm{DBH}\right)}\right)$	[91]
3	$\frac{\mathrm{a}}{1+\mathrm{b}\times { \mathrm{e}}^{-\mathrm{c}\times \mathrm{DBH}}}$	[92]

). We used the Levenberg–Marquardt algorithm [89] to solve these models, i.e., to find the optimal values of the model parameters. Also, since we required initial values to find the optimal parameters, 500 initial values were randomly considered for each model to ensure the global optimum was found. We chose the optimal value according to the evaluation criteria of root mean square error (RMSE) and R²_adj.

In OLS fits, one expects the model to improve as more variables are added but at the cost of overfitting. We used ANOVA (analysis of variance) to determine whether a relationship existed between two or more independent variables and a dependent variable. ANOVA can be used to determine whether the coefficients of the regression model are statistically significant and whether the differences among the independent variables are significant. It is also used to evaluate the overall effect of the independent variables on the dependent variable and to compare the relative importance of each of the independent variables on the dependent variable [93]. We used both the Akaike information criterion (AIC) and the Bayesian information criterion (BIC) to select variables [94]. Moreover, S, which is defined as $\sqrt{\frac{\mathrm{SSE}}{\mathrm{DF}}}$, denotes the variance of errors around the regression line, where SSE and DF are the sum of the squared residuals and the degree of freedom, respectively [95]. The R² from $1-\frac{\mathrm{SSE}}{\mathrm{SST}}$ demonstrates the regression model’s overall performance, where SST is the treatment sum of squares [96]. It is crucial to take a penalty for adding a variable to the model into account. Hence, we used the R²_adj and R²_pred coefficient and 10-fold R² and S to evaluate how well each step performed. The coefficient R²_adj is derived from 1-MSE/MST and R²_pred is derived from $1-\frac{\mathrm{PRESS}}{\mathrm{SST}},$ where PRESS is a measure of the difference between the predicted and observed values. To calculate PRESS, each observation is systematically removed from the data set, the regression equation is estimated, and the accuracy of the model in predicting the removed observation is determined. The 10-fold S measures the average difference between the actual values and the predicted values in the test data set. It provides an indication of how well the model fits the data. The 10-fold R² measures the proportion of variability in the response variable that is explained by the model. It indicates the strength of the relationship between the predictors and the response variable. The degree of freedom, which is a function of n (number of samples) and p (number of parameters), is included in the denominator of MSE and MST; thus, it is easy to see whether adding a variable to the model has a negative impact [95]. Therefore, it is necessary to use this coefficient in regression equations that have more than one variable [96]. Also, if there is a small difference between R²_adj and R², R² can be trusted, with a high probability that this difference is small in all stages.

2.3.3. Outlier Detection for OLS Models

We developed regressions with the independent variable H (height) for various Galajar ages. It was vital to identify and eliminate outliers prior to model evaluation, as noise is a constant companion of UAV-based photogrammetric products. The difference in fit (DFFITS) is an effective metric for detecting outlier data [94]; it is a mixture of the leverage and R-Student metrics and evaluates the effect of each observation on the fitted values of the regression model [97,98]. The ith data enters the model twice, once as test data and once as training data, followed by comparison of the predicted values for the two variants. Although this is somewhat time-consuming in practice, it was proven that DFFITS may be calculated using only the diagonal values of the hat matrix without re-modelling each data point [94,99]. Here, we identified outliers as data with $\mathrm{DFFITS}\ge 2\times \sqrt{\frac{\mathrm{p}}{\mathrm{n}}}$. The p and n represent the number of regression terms and samples, respectively, which were 6 and 796, respectively. Therefore, the data with $\mathrm{DFFITS}\ge 0.1736$ were discarded. Cook’s distance criterion [97,100] was the second strategy that was employed to filter outliers. The Cook’s distance criterion functions similarly to the DFITS criterion, with the exception that it evaluates the changes in the regression coefficients. Here, we calculated the aforementioned criterion for all of the data and identified outlier data as those with a $\mathrm{Cook}\prime \mathrm{s} \mathrm{distance}\ge 1$.

3. Results

3.1. DBH Estimation

3.1.1. Finding Effective Variables in DBH Estimation

We used the backward stepwise selection to eliminate redundant or ineffective variables, with the p-value used as the selection criterion, i.e., the variable with the greatest p was eliminated in each phase. The applied termination condition was alpha-to-remove value = 0.1. The results are summarized in

Table 7. The results of variable screening in the backward stepwise selection for inlier data. Height (H), number of trunks (nTr), trunk (Tr), species (Sp), and pollarding age (Age) were the candidate terms. nTr denotes the number of tree trunks, while the categorical variable Tr distinguishes between the single-stem and multi-stem trees.

Parameter	p-Value
	Step 1	Step 2	Step 3	Step 4	Step 5	Step 6	Step 7	Step 8	Step 9
Constant
H	0.000	0.000	0.000	0.000	0.000	0.000	0.000	0.000	0.000
nTr	0.305	0.396	0.166	0.266	0.204	----	----	----	----
Sp	0.028	0.013	0.000	0.000	0.000	0.000	0.000	0.000	0.000
Age	0.000	0.000	0.000	0.000	0.000	0.000	0.000	0.000	0.000
Tr	0.285	0.346	0.067	0.102	0.078	0.155	0.177	0.887	----
H × nTr	0.170	0.200	0.170	0.393	----	----	----	----	----
H × Sp	0.000	0.000	0.000	0.000	0.000	0.000	0.000	0.000	0.000
H × Age	0.000	0.000	0.000	0.000	0.000	0.000	0.000	0.000	0.000
H × Tr	0.049	0.056	0.039	0.079	0.045	0.040	0.174	----	----
nTr × Sp	0.348	0.408	0.309	----	----	----	----	----	----
nTr × Age	0.829	----	----	----	----	----	----	----	----
Sp × Tr	0.470	0.433	----	----	----	----	----	----	----
Age × Tr	0.769	0.012	0.011	0.023	0.028	0.041	----	----	----
S	3.1383	3.1350	3.1344	3.1351	3.1346	3.1358	3.1446	3.1463	3.1443
R2 (%)	58.9	58.9	58.8	58.6	58.6	58.5	58.2	58.1	58.1
R2adj (%)	57.8	57.8	57.9	57.8	57.9	57.8	57.6	57.5	57.6
AICc	4104	4100	4098	4096	4095	4094	4097	4097	4095
BIC	4211	4198	4186	4175	4169	4164	4157	4152	4146

.

Tree species and pollarding age were the most influential variables for OLS fitting, along with height (H) (Table 7). Furthermore, observing the effect of each of these two variables on the R² and R²_adj coefficients revealed that eliminating the species variable had a negligible impact on the model accuracy. Therefore, only the pollarding age variable was included in the final model. We also evaluated other nonlinear models (Table 6) to ensure that our linear model was the appropriate choice (

Table 8. Comparing the performance of nonlinear allometries with linear models for inlier data.

Model	Pollarding Stage	Estimated Parameter			Evaluation Criteria
		a	b	c	DF	SSE	RMSE	R2adj
$1$	Korpa	9.267	0.050	1.512	334	759.172	1.508	0.300
	Kor	6.234	0.099	2.701	198	39.224	0.445	0.729
	2 Khert	9.659	0.043	1.356	253	104.486	0.643	0.747
$2$	Korpa	7.972	2.872	0.084	334	759.004	1.507	0.300
	Kor	5.958	4.012	0.121	198	39.294	0.445	0.728
	2 Khert	7.724	2.842	0.084	253	101.395	0.633	0.754
$3$	Korpa	8.820	4.335	0.110	334	758.747	1.507	0.301
	Kor	7.082	5.493	0.146	198	39.356	0.446	0.728
	2 Khert	8.497	4.255	0.113	253	100.750	0.631	0.756
$\mathrm{a} +$ $\mathrm{b}\times \mathrm{DBH}$	Korpa	1.764	0.205		335	763.590	1.510	0.300
	Kor	1.328	0.203		199	41.646	0.457	0.715
	2 Khert	1.637	0.207		254	107.710	0.651	0.741

).

As shown in Table 8, the accuracies for different models were generally in the same range despite the model complexities applied. Therefore, we used the DBH~H linear model, which, in addition to its simplicity and comprehensibility, provides better interpretability than nonlinear models.

3.1.2. Outlier Detection/OLS Model

Figure 4 shows the outlier samples and residual errors identified by the criteria that were described in Section 2.3.3.

The final OLS was re-fitted after eliminating 35 trees (almost 4.4% of all trees) as outlier data. The comparison of residual errors to the fitted values is shown in Figure 5.

The estimated coefficients for this equation are listed in

Table 9. The results of the analysis of the coefficients of the linear model.

Term		Symbol	Coef	SE Coef.	95% CI	T-Value	p-Value	VIF
Constant		${\mathsf{\alpha }}_0$	10.852	0.581	(9.712, 11.992)	18.69	0.000
Age	2 Khert	${\mathsf{\alpha }}_1$	−12.36	1.01	(−14.35, −10.37)	−12.21	0.000	17.28
	Kor	${\mathsf{\alpha }}_2$	−9.83	1.40	(−12.57, −7.08)	−7.02	0.000	28.71
H		${\mathsf{\alpha }}_3$	1.4677	0.0972	(1.2769, 1.6585)	15.10	0.000	1.62
H × Age	2 Khert	${\mathsf{\alpha }}_4$	2.111	0.185	(1.748, 2.474)	11.41	0.000	16.46
	Kor	${\mathsf{\alpha }}_5$	1.896	0.278	(1.350, 2.442)	6.82	0.000	27.36

, while Figure 6 illustrates the relationship between the tree height and DBH by age. A 95% confidence interval for the coefficients was calculated to estimate whether the model might also be valid across other regions beyond the parametrization domain of the six test sites.

The general equation was:

$$\mathrm{DBHMax}={\mathsf{\alpha }}_0+{\mathsf{\alpha }}_1{\mathrm{x}}_1+{ \mathsf{\alpha }}_2{\mathrm{x}}_2+{\mathsf{\alpha }}_3\mathrm{H}+{\mathsf{\alpha }}_4{\mathrm{x}}_1\mathrm{H}+{\mathsf{\alpha }}_5{\mathrm{x}}_2\mathrm{H}.$$

(1)

The corresponding values for the variables of ${\mathrm{x}}_1$ and ${\mathrm{x}}_2$ are presented in

Table 10. The x₁ and x₂ coding and the final equation of linear regression for various pollarding ages.

Age	${\mathbf{x}}_1$	${\mathbf{x}}_2$	Equation (DBHMax=)
Korpa	0	0	−10.85 + 1.48 H
Kor	0	1	1.02 + 3.37 H
2 Khert	1	0	−1.51+ 3.58 H

.

Table 11. Summary of linear regression model evaluation.

S	R2	R2adj	PRESS	R2pred	AICc	BIC	10-Fold S	10-Fold R2
3.21095	56.06%	55.78%	8267.26	55.40%	4124.24	4156.86	3.22838	55.25%

provides a summary of the model evaluation based on multiple criteria, while

Table 12. Results of the ANOVA.

Source	DF	Seq SS	Contribution	Adj SS	Adj MS	F-Value	p-Value
Regression	5	10,393.00	56.06%	10,393.0	2078.60	201.61	0.000
H	1	8594.60	46.36%	2351.3	2351.28	228.05	0.000
Age	2	217.80	1.18%	1721.6	860.80	83.49	0.000
H × Age	2	1580.50	8.53%	1580.5	790.26	76.65	0.000
Error	790	8145.00	43.94%	8145.0	10.31
Lack of fit	785	8110.00	43.75%	8110.0	10.33	1.48	0.359
Pure error	5	35.00	0.19%	35.0	7.00
Total	795	18,538.00	100.00%

summarizes the analysis of variance (ANOVA).

Table 12 shows the analysis of variance of the model. The model included five parameters, in which high values of F and very low values of p emphasized their significance. The lack-of-fit error (indicated by the high value of p) increased due to the discrete nature of the DBH values (due to the cm measurement accuracy) and the continuous nature of the H values.

4. Discussion

The general objective of our study was to examine whether height derived from UAV photogrammetry can predict the DBH of pollarded oak trees, which is challenging due to the strongly altered canopy structures [27,67]. In the following, we first discuss the technical aspects of UAV-based tree height determination and then reflect on the identified DBH~H relationships.

Earlier studies discussed the challenges of producing realistic 3D models in forests using UAV photogrammetry [101,102]. One challenge may be that tree crowns are observed only with a limited number of pixels in UAV images. This complicates solving the collinearity equations of automatic image-matching techniques [102]. Increasing the number of sampled pixels has been suggested in such conditions [102]. Tu et al., (2021) showed that the combined use of oblique and nadir flights (as compared to a nadir flight alone) can improve 3D photogrammetry models by 35% [103]. However, this approach might also be associated with image-matching issues such as changes in wind speed/direction and the sun exposure that may occur due to the time lag between flights, hence complicating the image-matching procedure. Here, we coped with these two issues by taking images using continuous imaging under constant and calm wind conditions. One may also note that pollarded trees are less sensitive to wind because of the generally sparser canopy foliage compared with unpollarded trees. Additionally, we conducted the flights with the shortest possible time lag of approximately one minute between the two flights, due to which changes in sun exposure could be regarded as negligible. In order to overcome many common image-matching problems in forests, we processed images of two oblique (60°) and vertical flights simultaneously to increase the number of candidate pixels in the image-matching process [56].

All in all, our analysis returned accurate estimates of the tree height using UAV photogrammetric models, with an average difference of 0.3 m to 1.2 m compared to the values extracted from the CHM and the field measurement. This was in the range reported by other studies, which obtained differences of between 0.3 m to 2.9 m [104] and 0.1 m to 0.5 m [57] when estimating tree height from UAV data. The average heights of pollarded trees in our study were in the range of 3.5–6.7 m, which can be considered at the lower end of the height values reported in the literature for these silvopastoral systems [6,18,32]. The tree height of pollarded trees is strongly influenced by the local inhabitants who practice pollarding. The traditional owners of Galajars do not allow the crown of the tree to be out of their reach due to the increase in height [30]. Furthermore, another reason for the low height of these trees may be the low canopy cover of the stands. As a result, there is less competition between trees to receive light (a factor influencing height growth) in pollarded stands, and thus the trees tend to expand their crown horizontally rather than vertically (height growth) [105].

In terms of height–DBH relationships, previous studies on other forest types generally showed that DBH is correlated with height [106,107]. Ref. [29] showed that the DBHs and heights of pollarded trees of the species Quercus cerris L. were linearly correlated (correlation coefficient = 0.464). Furthermore, Niemczyk and Bruchwald (2017) presented relationships for estimating the DBH of Populus spp. trees and hybrids in coppice stands using height and the form factor [108]. Finally, Ref. [109] investigated and compared the relationship between the DBH and the height of Populus spp. in short rotations in northern Poland in two cycles of 5 and 6 years.

In field inventories, height measurement is generally more difficult than DBH measurement [110], often with lower accuracy of measured height than that of DBH [38,111]. Thus, allometries are often used to estimate height using the DBH in classical inventories of standard forests (see Table 6). We argue that the contrary is true in our case and other comparable cases in stands comprising both single- and multi-stem coppice structures. The use of photogrammetric UAV allows for the estimation of tree heights with reasonable accuracies, whereas DBH measurements may be more time-consuming or prohibited. Therefore, we estimated DBH using the UAV-derived tree heights for single- and multi-stem individuals in this study. In our special case, we used the maximum DBH of pollarded coppice groups (where trees are generally associated with more than one stem) as previously indicated in Section 2.2. We first checked whether distinguishing between single- and multi-stem trees (categorical variable of Tr in Table 7) improved the model quality. The high p-value for the Tr variable throughout the variable selection suggested that this categorization was not helping the model. As a result, we adopted the idea of using the maximum DBH.

The maximum DBH may also be suitable from an ecological perspective. In multi-stemmed coppice oak clumps, dominant stems benefit more from the well-developed parent-tree root system. Therefore, the dominant stem with maximum DBH typically grows faster both in diameter and height. Although multiple stems often persist within a clump for decades [112], as coppice stems grow larger, competition within the clump increases and reduces the diameter and height growth of stems. Finally, only one or two stems within each clump usually maintain their dominance in the crown position, and other stems are oppressed or excluded due to natural thinning or traditional silvicultural treatments [113]. Thus, we focused on establishing an allometry between the maximum DBH of the coppice oak trees and our UAV-derived H measurements.

In one of the few other studies focusing on DBH~H relationships in coppice stands, Niemczyk et al., (2016) used the allometric relationship suggested by [114] to estimate height from DBH. This relationship was nonlinear and included two fitted coefficients that were obtained by the samples measured in short rotations in northern Poland [109]. In northern Zagros, such allometric equations do not exist, especially not for different stages of pollarding. Therefore, we applied a linear model to explore the link between DBH and height (Figure 5), in which we used the SE coefficient, T-value, and VIF to evaluate the model coefficients (Table 9). The SE coefficient explains how using different sampling techniques might lead to changes in the estimated coefficients [95]. Here, we observed a higher SE coefficient for the pollarding age (Age) class 2 Khert, i.e., Galajar’s DBH at the 2 Khert age was estimated more accurately. These linear models are not ecologically realistic because the increase in tree height decreases with time, especially when the tree reaches the forest canopy [115]. But this is only true if the forests are minimally disturbed, whereas our studied forests are regularly changed by humans (see Figure 1). However, we additionally tested nonlinear allometries (Table 6). Although the findings demonstrated that nonlinear models exhibited slightly better accuracy than linear models, we chose to utilize the linear model. This decision was based on the simplicity, comprehensibility, and superior interpretability of linear models compared to nonlinear models.

In our study sites, pollarding results in the complete debranching of oak trees. Therefore, one year later, pollarding trees have almost no crown and their height decreases by 1 to 1.5 m. With time, the tree crown recovers by producing new twigs and branches both apically and laterally [33]. Just 4 years after pollarding (2 Khert class), the tree crown and height reach their typical shape and condition. Thus, the relationship between DBH and height is more reliable in such stands.

The linear model (Table 11) showed a reliable R² in the overall assessment due to the small difference between R² and R²_adj. The R²_pred is essentially a validation that denotes the degree of model reliability [116]. Furthermore, comparing k-fold S and k-fold R² with S and R², respectively, can warn of overfitting. We set k to 10, for which the small difference between the values rejected the hypothesis of model overfitting. In addition, the ANOVA (Table 12) showed relatively high SSError, i.e., changes in DBHMax that cannot be explained by the model. Since we chose DBHMax as the DBH value of the tree (except that the tree was multi-stem), it added an error to the model, which can be seen in SSError. We instead used DBHMax due to the fact that, in oak clumps, the dominant trees are significantly taller than the oppressed trees and stay dominant until all others are thinned or die.

The linear model suggested in this study allows for the prediction of DBH with reasonable accuracies in coppice stands at the age of Kor and 2 Khert, but the accuracy for stands in the Korpa age class is notably worse. Given that the models are also nearly unbiased, there seems to be a certain potential to inventory the DBH and hence also the biomass of oak coppice stands in Northern Zagros. Given the large areas covered by these coppice stands and given the frequently underestimated forest cover of semi-arid forests in global forest maps [117], this may be a relevant contribution to improve our understanding of global biomass distribution patterns and the corresponding carbon stocks. In future studies, it may be worthwhile to think about developing automated approaches to identify the pollarding age of coppice stands from remote sensing data. This may be possible with very-high-resolution UAV data but could also examined for other VHR data (if available).

5. Conclusions

In our research across silvopastoral systems (pollarded stands or Galajars) in Iran, we demonstrated that tree height can be effectively estimated at the scale of individual pollarded trees using UAV data. However, UAV photogrammetry is unable to directly derive DBH, which led us to demonstrate how DBH for Galajars of various ages may be estimated by tree height. We tested well-known nonlinear and linear allometries to model DBH~H and compared their performance, which suggested that DBH can be estimated with practically sufficient accuracy by solely using OLS with the independent variables of height and oak species across Galajar ages. We showed that removing species from the models was only associated with a marginal effect on the model accuracy; thus, we did not use species as a variable in the models. That is, we suggest eliminating species information when conducting pollarding-related field measurements, which would lead to significant savings in expert and logistics costs. We also suggest conducting a comparative analysis on pollarded and less-disturbed stands. Last but not least, one may note that we missed the Khert age group in our six inventoried sites, which calls for additional surveys on this age group to complement our results.