Aboveground Tree Biomass Estimation of Sparse Subalpine Coniferous Forest with UAV Oblique Photography

: In tree Aboveground Biomass (AGB) estimation, the traditional harvest method is accurate but unsuitable for a large-scale forest. The airborne Light Detection And Ranging (LiDAR) is superior in obtaining the point cloud data of a dense forest and extracting tree heights for AGB estimation. However, the LiDAR has limitations such as high cost, low efﬁciency, and complicated operations. Alternatively, the overlapping oblique photographs taken by an Unmanned Aerial Vehicle (UAV)-loaded digital camera can also generate point cloud data using the Aerial Triangulation (AT) method. However, limited by the relatively poor penetrating capacity of natural light, the photographs captured by the digital camera on a UAV are more suitable for obtaining the point cloud data of a relatively sparse forest. In this paper, an electric ﬁxed-wing UAV loaded with a digital camera was employed to take oblique photographs of a sparse subalpine coniferous forest in the source region of the Minjiang River. Based on point cloud data obtained from the overlapping photographs, a Digital Terrain Model (DTM) was generated by ﬁltering non-ground points along with the acquisition of a Digital Surface Model (DSM) of Minjiang ﬁr trees by eliminating subalpine shrubs and meadows. Individual tree heights were extracted by overlaying individual tree outlines on Canopy Height Model (CHM) data computed by subtracting the Digital Elevation Model (DEM) from the rasterized DSM. The allometric equation with tree height (H) as the predictor variable was established by ﬁtting measured tree heights with tree AGBs, which were estimated using the allometric equation on H and Diameter at Breast Height (DBH) in sample tree plots. Finally, the AGBs of all of the trees in the test site were determined by inputting extracted individual tree heights into the established allometric equation. In accuracy assessment, the coefﬁcient of determination (R 2 ) and Root Mean Square Error (RMSE) of extracted individual tree heights were 0.92 and 1.77 m, and the R 2 and RMSE of the estimated AGBs of individual trees were 0.96 and 54.90 kg. The results demonstrated the feasibility and effectiveness of applying UAV-acquired oblique optical photographs to the tree AGB estimation of sparse subalpine coniferous forests.


Introduction
Biomass refers to the amount of material accumulated by plants in a unit area, and its essence is the organic matter and energy accumulated by the photosynthesis of green plants through their assimilation organs [1]. Forest biomass, especially tree biomass, is an important index to measure the carbon sequestration capacity of a forest and evaluate its carbon budget [2]. Accurate and quantitative estimations of forest biomass also form the data basis for study of the carbon cycle of a terrestrial ecosystem. As the underground portion of forest biomass is difficult to obtain, the Aboveground Biomass (AGB) of forest is usually used instead.
The tree AGB of a forest can be obtained through direct measurement or indirect estimation [3]. Direct measurement involves the widely used traditional harvest method. It can be divided into clear-felling, the average wood method, and the relative growth method [4]. Indirect estimations of tree AGB are mainly achieved by the biomass model, which can be classified into the single tree model and large-scale forest model [5]. However, direct measurement methods are time-consuming and prohibitively expensive, and have difficulties in the tree AGB estimation of a large-scale forest. Optical and thermal infrared remote sensing for tree AGB estimation has the advantage of large spatial coverage and multi-temporal observation. However, the accuracy of tree AGB estimation is relatively low, as there are problems of low image resolution, remote sensing signal response saturation, data conversion between different scales, and data heterogeneity [6]. Tree AGB estimation of a forest with Synthetic Aperture Radar (SAR) [7][8][9][10][11] and Light Detection And Ranging (LiDAR) [12][13][14][15] is of relative high accuracy, but the equipment involved is very expensive, and the data acquisition process is usually time-consuming and risky if the equipment is mounted on an airplane or a helicopter. Therefore, they are not suitable for wide deployments and routine operations within a large-scale forest.
In recent years, the Unmanned Aerial Vehicle (UAV) or drone has been gradually utilized in tree AGB estimations of forests, with the advantages of low cost, flexible take-off and landing, safety, under-cloud flying, and hyperspatial image resolution [16]. However, most of the current literature in this respect is focused on UAV-loaded LiDARs and the application of the acquired point cloud data to tree height extraction and AGB estimation [17][18][19][20][21][22]. Actually, the overlapping oblique photographs taken by a UAV-loaded optical digital camera can also generate point cloud data by the Aerial Triangulation (AT) method or Structure of Motion (SfM) algorithm [23][24][25][26][27]. In addition, the LiDAR must be mounted on a multi-rotor UAV of low flight speed, short duration, and small spatial coverage. However, the digital camera can instead be mounted on a fixed-wing UAV of high flight speed, long duration, and large spatial coverage. LiDAR has a good penetrating capability, which makes it superior in acquiring the point cloud data of a dense forest. The optical digital camera is more suitable for obtaining the point cloud data of a relatively sparse forest due to the relatively poor penetrating capacity of natural light [28,29].
In this study, an electric fixed-wing UAV loaded with a digital camera was employed to acquire oblique photographs of a sparse subalpine coniferous forest in the source region of the Minjiang River. Based on point cloud data produced from the overlapping photographs with the AT method, a Digital Surface Model (DSM) and Digital Terrain Model (DTM) of the test site with the low vegetation removed were obtained. Then, the AGBs of the Minjiang fir trees were estimated based on individual tree heights extracted from the computed Canopy Height Model (CHM) and the allometric equation using tree height as a predictor. The results indicated the feasibility and effectiveness of applying the fixed-wing UAV loaded with a digital camera to the tree AGB estimation of sparse subalpine coniferous forests.

Study Area
The test site (Figure 1b (Figure 1a). This area is also the source region of the Minjiang River, which is one of the important tributaries of the upper Yangtse River. The Minjiang River has a length of about 735 km, with the majority of the river located in Sichuan Province. As shown in Figure 1c, the monoculture tree type of the dark subalpine coniferous forest is Minjiang fir (Abies faxoniana Rehder & E.H.Wilson) [30,31]. In addition, there are also subalpine shrubs and meadows. The dominant soil types are brown coniferous forest soil and subalpine meadow soil. This area belongs to the high mountain and canyon region between the Sichuan Basin and the Tibetan Plateau, with elevation variations from 3400 m to 3700 m, surface slopes up to 30 • , and annual accumulated temperatures of 3000 • C to 4500 • C [32,33].
Remote Sens. 2018, 10, x FOR PEER REVIEW 3 of 19 about 735 km, with the majority of the river located in Sichuan Province. As shown in Figure 1c, the monoculture tree type of the dark subalpine coniferous forest is Minjiang fir (Abies faxoniana Rehder & E.H.Wilson) [30,31]. In addition, there are also subalpine shrubs and meadows. The dominant soil types are brown coniferous forest soil and subalpine meadow soil. This area belongs to the high mountain and canyon region between the Sichuan Basin and the Tibetan Plateau, with elevation variations from 3400 m to 3700 m, surface slopes up to 30°, and annual accumulated temperatures of 3000 °C to 4500 °C [32,33].

UAV System
As shown in Figure 2, the employed UAV oblique photographing system consists of four major components: the electric fixed-wing UAV platform, the flight control system, the digital camera, and the ground control system. The UAV platform has a 1.2-m wingspan, 0.8-m fuselage length, and 4.2kg working weight. The model of both motors used in the UAV is SNNNYSKY V3508, with a KV value of 580. The flight control system is ANXIANG 2012, which mainly controls and stabilizes the takeoff/landing, flight altitude, and heading and attitude of the UAV platform. The POS data is recorded by an Inertial Measurement Unit (IMU) and Global Positioning System (GPS) embedded in the flight control system. The digital camera is a SONY ILCE-5100 with a resolution of 6000 × 4000 pixels, a size of 75 mm × 63 mm × 50 mm, and a weight of 192 g, including the lens. When the camera is mounted on the UAV, the angle between the lens axis and the line perpendicular to the ground is set to 20°. This aids in obtaining more side information about the trees, while also helping to improve the precision of three-dimensional (3D) point cloud generation and tree height extraction. The major functions of the ground control system include trajectory planning and uploading, flight attitude data downloading, and remote control commanding. The fixed-wing UAV has the advantages of a long duration, fast flight speed, strong anti-shear force, and stable and safe operation. Therefore, it is very suitable for aerial photographing in high altitudes and complex mountain environments.

UAV System
As shown in Figure 2, the employed UAV oblique photographing system consists of four major components: the electric fixed-wing UAV platform, the flight control system, the digital camera, and the ground control system. The UAV platform has a 1.2-m wingspan, 0.8-m fuselage length, and 4.2-kg working weight. The model of both motors used in the UAV is SNNNYSKY V3508, with a KV value of 580. The flight control system is ANXIANG 2012, which mainly controls and stabilizes the takeoff/landing, flight altitude, and heading and attitude of the UAV platform. The POS data is recorded by an Inertial Measurement Unit (IMU) and Global Positioning System (GPS) embedded in the flight control system. The digital camera is a SONY ILCE-5100 with a resolution of 6000 × 4000 pixels, a size of 75 mm × 63 mm × 50 mm, and a weight of 192 g, including the lens. When the camera is mounted on the UAV, the angle between the lens axis and the line perpendicular to the ground is set to 20 • . This aids in obtaining more side information about the trees, while also helping to improve the precision of three-dimensional (3D) point cloud generation and tree height extraction. The major functions of the ground control system include trajectory planning and uploading, flight attitude data downloading, and remote control commanding. The fixed-wing UAV has the advantages of a long

DOM of Test Site
As shown in Figure 2d, the overhead crossed trajectory for the UAV operation had an 80% longitudinal overlap and a 60% lateral overlap to meet the accuracy requirements of the later AT processing. The UAV flight was carried out at an average relative altitude of 400 m, with 476 oblique optical photographs and an average spatial resolution of 0.05 m being acquired. The AT processing of raw oblique photographs was completed using Bentley ContextCapture (CC) 4.3, with additional input data including 3D coordinates of Ground Control Points (GCPs), POS data, and a camera calibration model for the digital camera [34]. The raw oblique photographs were divided into rectangular blocks. For each block, the major steps included establishing GCPs; performing interior orientation; measuring and transferring all of the tie points, check points, and GCPs appearing on all of the photographs; and performing a least squares bundle adjustment. The AT processing ultimately produced exterior orientation parameters for the photographs and 3D coordinates for the measured object points, namely DSMs. Photographs of each block were then orthorectified and mosaicked to produce the resulting Digital Orthophoto Map (DOM) (Figure 3a).

DOM of Test Site
As shown in Figure 2d, the overhead crossed trajectory for the UAV operation had an 80% longitudinal overlap and a 60% lateral overlap to meet the accuracy requirements of the later AT processing. The UAV flight was carried out at an average relative altitude of 400 m, with 476 oblique optical photographs and an average spatial resolution of 0.05 m being acquired.
The AT processing of raw oblique photographs was completed using Bentley ContextCapture (CC) 4.3, with additional input data including 3D coordinates of Ground Control Points (GCPs), POS data, and a camera calibration model for the digital camera [34]. The raw oblique photographs were divided into rectangular blocks. For each block, the major steps included establishing GCPs; performing interior orientation; measuring and transferring all of the tie points, check points, and GCPs appearing on all of the photographs; and performing a least squares bundle adjustment. The AT processing ultimately produced exterior orientation parameters for the photographs and 3D coordinates for the measured object points, namely DSMs. Photographs of each block were then orthorectified and mosaicked to produce the resulting Digital Orthophoto Map (DOM) ( Figure 3a). As shown in Figure 3b, we used only the test site area of the DOM for subsequent analysis, as it was where sample tree plots were situated.

DOM of Test Site
As shown in Figure 2d, the overhead crossed trajectory for the UAV operation had an 80% longitudinal overlap and a 60% lateral overlap to meet the accuracy requirements of the later AT processing. The UAV flight was carried out at an average relative altitude of 400 m, with 476 oblique optical photographs and an average spatial resolution of 0.05 m being acquired. The AT processing of raw oblique photographs was completed using Bentley ContextCapture (CC) 4.3, with additional input data including 3D coordinates of Ground Control Points (GCPs), POS data, and a camera calibration model for the digital camera [34]. The raw oblique photographs were divided into rectangular blocks. For each block, the major steps included establishing GCPs; performing interior orientation; measuring and transferring all of the tie points, check points, and GCPs appearing on all of the photographs; and performing a least squares bundle adjustment. The AT processing ultimately produced exterior orientation parameters for the photographs and 3D

Field Data
The arboreal layer of the forest within the test site is characterized by Minjiang fir trees. There are also subalpine shrubs and meadows in the forest. Five 30 m × 30 m sample tree plots were selected in the test site, with plot 5 being used specially to assess the results. The heights and Diameters at Breast Height (DBH) of all of the Minjiang fir trees in all five plots were measured with a laser rangefinder and a DBH ruler, and are listed in Appendix A. The 3D coordinates of the 24 GCPs were measured with the Trimble R10 (8 mm H/15 mm V) for AT processing. The specific 3D coordinate data of the GCPs in Figure 3a are listed in Appendix B.

Methodology
In order to determine the tree AGBs of the sparse subalpine coniferous forest directly based on individual tree heights, a new allometric equation with tree height (H) as the predictor variable needed to be established. First, we obtained 3D point cloud data of the test site from UAV oblique photographs using the AT method, and calculated the CHM data by subtracting the DTM from the DSM without low vegetation. These were acquired by point cloud filtering and classification, respectively. Secondly, the DOM of the test site was segmented and classified to acquire the canopy outline of each individual tree, after which individual tree heights were extracted by overlaying tree canopy outlines on the CHM data. Thirdly, the AGB of each Minjiang fir tree in plots 1-4 was obtained by inputting the measured height and DBH into the allometric equation with variables of H and DBH; the coefficient values of the allometric equation with the unique variable of H were obtained through power function fitting. Finally, the AGBs of all of the trees in the test site were acquired by inputting the extracted individual tree heights into the newly established allometric equation on H. The flowchart is shown in Figure 4. As shown in Figure 3b, we used only the test site area of the DOM for subsequent analysis, as it was where sample tree plots were situated.

Field Data
The arboreal layer of the forest within the test site is characterized by Minjiang fir trees. There are also subalpine shrubs and meadows in the forest.   In order to determine the tree AGBs of the sparse subalpine coniferous forest directly based on individual tree heights, a new allometric equation with tree height (H) as the predictor variable needed to be established. First, we obtained 3D point cloud data of the test site from UAV oblique photographs using the AT method, and calculated the CHM data by subtracting the DTM from the DSM without low vegetation. These were acquired by point cloud filtering and classification, respectively. Secondly, the DOM of the test site was segmented and classified to acquire the canopy

Point Cloud Acquisition and Denoising
A point cloud is a set of points with both 3D coordinates and spatial sampling points of the object surface. The 3D point cloud data of the subalpine coniferous forest was obtained from the overlapping UAV oblique photographs. Based on the 3D model achieved by AT processing in Section 2.2.2, the point cloud data of the test site could be generated by undertaking the 3D reconstruction procedure in CC [34].
Point cloud noise can be classified as high or low noise according to the relative positions of noise points to normal points [35]. In the study, the point cloud obtained using the AT method had high noise. The high noise was possibly caused by false image matching, which was notably higher than the surrounding average elevation. The method for distinguishing this involved calculating the median and standard deviation of all of the point heights in the neighborhood determined by a certain searching radius for each point. If the point was higher than the specified height difference of the surrounding points, it was judged as a high noise point, and then removed [36].

CHM Computation
The point cloud data of a DSM usually includes ground points, as well as non-ground surface points of artificial buildings, vegetation, and vehicles [37]. Point cloud filtering serves to remove non-ground points and extract points that represent real ground. In the study, the gradual densification algorithm was implemented in Terrasolid, which is a professional point cloud processing software used for point cloud filtering. It first seeks the lowest point in each grid cell, takes them as surface points, and uses them as initial seed points to generate a sparse Triangulated Irregular Network (TIN) model. Then, the initial TIN model was raised to achieve the final Digital Terrain Model (DTM). The iterative angle (Figure 5a) is the maximum angle between the line that connects the target point and the nearest vertex of the TIN model, and the local triangular plane of the TIN model. The smaller the iterative angle, the smaller the surface fluctuation of the point cloud. The iterative angle usually has a small set value of about 4 • in flat areas, and a large value of about 10 • in mountainous areas. The iterative distance (Figure 5a) is the distance from the target point to the local triangular plane, which helps to exclude low height objects from the TIN model. Its typical value is 0.5-1.5 m, and it is also dependent on flat to mountainous terrain features. If the iterative angle and the iterative distance are less than the specified values, the target point will be added to the TIN model. Each added point brings the TIN model closer to the real ground, with the process repeating until no new ground points are found.
The vegetation of the test site consisted of Minjiang fir trees in the arboreal layer and low vegetation, including subalpine shrubs and meadows without any artifacts such as houses and wires ( Figure 5b). To obtain a DSM containing both the Minjiang fir trees and the ground, it was necessary to eliminate the disturbance of low vegetation. According to our survey of sample tree plots, low vegetation was usually not higher than three meters, while Minjiang fir trees were no less than three meters. Therefore, a height of three meters was determined as the demarcation line between Minjiang fir trees and low vegetation. The "classify by height from ground" algorithm was implemented in Terrasolid to classify and remove low vegetation according to the three-meter height between the target points and the DTM.
The CHM data could be obtained by subtracting the DTM from the corresponding DSM without low vegetation, which is illustrated in Figure 5c The CHM data could be obtained by subtracting the DTM from the corresponding DSM without low vegetation, which is illustrated in Figure 5c-e. Figure 5c was an example patch of the DSM as the result of eliminating low vegetation, and Figure 5d was the point cloud of the DTM. Figure 5e was the difference between the DSM and DTM, namely the CHM point cloud of the Minjiang fir trees in the patch.

Extraction of Individual Tree Canopy and Height
Individual tree heights needed to be extracted from CHM data before tree AGBs could be estimated in the test site. We first applied an object-based image analysis approach using the eCognition ® Developer 8.9 software package to extract individual trees canopies, and then overlaid the spatial outlines of the individual tree canopy on the CHM data to obtain individual tree heights.
Object-based image analysis includes two major steps: multiscale image segmentation and object classification. A segmentation scale is the key to image segmentation. The Estimation of Scale Parameter (ESP) tool [38] was used to segment an image with fixed scale parameter increments. The optimal segmentation scale was objectively selected according to the values of Local Variance (LV) and Rates of Change of LV (ROC) in the segmented results. ROC can assess the dynamics of LV from one object level to another. The scales at which ROCs reach local maxima will be regarded as candidate optimal segmentation scales. The candidate scale, which makes segmented objects the most separable, with the largest internal homogeneity and largest heterogeneity between objects, will be selected as the optimal segmentation scale. After the image was optimally segmented into objects, the nearest-neighbor classification method was utilized to classify the segmented objects into three classes: tree canopy, subalpine shrub, and subalpine meadow. The sample class images and their feature sets, which included features of spectrum, shape, and texture, were selected to participate in classification. Based on the sample images and the initial feature set of each class, the feature space was optimized to find the optimal feature set from more than 50 features, which made the most distinctive differences between the classes [39].
After segmentation and classification, we obtained the spatial outlines of individual Minjiang fir trees (Figure 6a). The spatial outlines were then overlaid on the rasterized CHM (Figure 6b), and the maximum heights in the outlines were extracted as individual tree heights. If there existed more than one pixel of maximum tree height, the pixel located in the center would be taken as the top of the tree

Extraction of Individual Tree Canopy and Height
Individual tree heights needed to be extracted from CHM data before tree AGBs could be estimated in the test site. We first applied an object-based image analysis approach using the eCognition ® Developer 8.9 software package to extract individual trees canopies, and then overlaid the spatial outlines of the individual tree canopy on the CHM data to obtain individual tree heights.
Object-based image analysis includes two major steps: multiscale image segmentation and object classification. A segmentation scale is the key to image segmentation. The Estimation of Scale Parameter (ESP) tool [38] was used to segment an image with fixed scale parameter increments. The optimal segmentation scale was objectively selected according to the values of Local Variance (LV) and Rates of Change of LV (ROC) in the segmented results. ROC can assess the dynamics of LV from one object level to another. The scales at which ROCs reach local maxima will be regarded as candidate optimal segmentation scales. The candidate scale, which makes segmented objects the most separable, with the largest internal homogeneity and largest heterogeneity between objects, will be selected as the optimal segmentation scale. After the image was optimally segmented into objects, the nearest-neighbor classification method was utilized to classify the segmented objects into three classes: tree canopy, subalpine shrub, and subalpine meadow. The sample class images and their feature sets, which included features of spectrum, shape, and texture, were selected to participate in classification. Based on the sample images and the initial feature set of each class, the feature space was optimized to find the optimal feature set from more than 50 features, which made the most distinctive differences between the classes [39].
After segmentation and classification, we obtained the spatial outlines of individual Minjiang fir trees (Figure 6a). The spatial outlines were then overlaid on the rasterized CHM (Figure 6b), and the maximum heights in the outlines were extracted as individual tree heights. If there existed more than one pixel of maximum tree height, the pixel located in the center would be taken as the top of the tree canopy. When the identified locations of individual tree tops were overlaid back, they matched the outlines of corresponding tree canopies (Figure 6c). In some cases, it was possible for tree shadows to be falsely classified as tree canopies; these could be conveniently removed, as the extracted heights of tree shadows were 0 m, and those of real trees were no less than 3 m. canopy. When the identified locations of individual tree tops were overlaid back, they matched the outlines of corresponding tree canopies (Figure 6c). In some cases, it was possible for tree shadows to be falsely classified as tree canopies; these could be conveniently removed, as the extracted heights of tree shadows were 0 m, and those of real trees were no less than 3 m.

Allometric Equation on Tree Height
According to the guidelines for measuring and monitoring the carbon sequestration of afforestation projects [40] , the allometric equation based on the H and DBH of Sichuan firs is as follows: where WT is tree AGB (kg), D represents DBH (cm), and H is tree height (m).
As Minjiang fir is one of the dominant subspecies of Sichuan fir, Formula 1 is applicable for the AGB estimation of Minjiang fir trees. The measured heights and DBHs of Minjiang fir trees in plots 1-4 were inputted into Formula 1 to obtain the AGBs of individual trees.
As tree heights were the only data that could be directly obtained by UAV oblique photography, the relationship between the tree height and its AGB needed to be established. The power function model (Formula (2)) was chosen as the allometric equation, which relied only on tree height for AGB estimation: = · (2) where W is the AGB of an individual Minjiang fir tree, H is the tree height, and α and β are coefficients of the allometric equation.
By fitting individual tree heights and the AGBs of individual Minjiang fir trees in plots 1-4 with Formula (2), the coefficient values of α and β were determined, and the allometric equation using tree height as a predictor for AGB estimation was established. Thus, the AGBs of all of the Minjiang fir trees in the test site could be determined by inputting their extracted tree heights.
The coefficient of determination (R 2 ) and the Root Mean Square Error (RMSE) were used to explore the correlation between extracted tree height and estimated tree AGB, and assess model fitness.

Computed CHM
The 3D point cloud was generated by sequentially applying AT processing and 3D reconstruction to UAV-acquired oblique photographs in CC. To remove high noise points, the point neighborhood was defined as a circle with a radius of 5 m, and the threshold of height difference was set to 0.5 m. A total of 25,391,900 points were ultimately obtained ( Figure 7a). As recommended by the user's guide of Terrascan [41], a module of Terrasolid, the iterative angle and iterative distance in mountainous area was set to 10° and 1.4 m respectively in non-ground point filtering. Then, the

Allometric Equation on Tree Height
According to the guidelines for measuring and monitoring the carbon sequestration of afforestation projects [40], the allometric equation based on the H and DBH of Sichuan firs is as follows: where W T is tree AGB (kg), D represents DBH (cm), and H is tree height (m).
As Minjiang fir is one of the dominant subspecies of Sichuan fir, Formula 1 is applicable for the AGB estimation of Minjiang fir trees. The measured heights and DBHs of Minjiang fir trees in plots 1-4 were inputted into Formula 1 to obtain the AGBs of individual trees.
As tree heights were the only data that could be directly obtained by UAV oblique photography, the relationship between the tree height and its AGB needed to be established. The power function model (Formula (2)) was chosen as the allometric equation, which relied only on tree height for AGB estimation: where W is the AGB of an individual Minjiang fir tree, H is the tree height, and α and β are coefficients of the allometric equation. By fitting individual tree heights and the AGBs of individual Minjiang fir trees in plots 1-4 with Formula (2), the coefficient values of α and β were determined, and the allometric equation using tree height as a predictor for AGB estimation was established. Thus, the AGBs of all of the Minjiang fir trees in the test site could be determined by inputting their extracted tree heights.
The coefficient of determination (R 2 ) and the Root Mean Square Error (RMSE) were used to explore the correlation between extracted tree height and estimated tree AGB, and assess model fitness.

Computed CHM
The 3D point cloud was generated by sequentially applying AT processing and 3D reconstruction to UAV-acquired oblique photographs in CC. To remove high noise points, the point neighborhood was defined as a circle with a radius of 5 m, and the threshold of height difference was set to 0.5 m. A total of 25,391,900 points were ultimately obtained ( Figure 7a). As recommended by the user's guide of Terrascan [41], a module of Terrasolid, the iterative angle and iterative distance in mountainous area was set to 10 • and 1.4 m respectively in non-ground point filtering. Then, the acquired DTM of the test site was rasterized, resulting in the Digital Elevation Model (DEM) with a spatial resolution of 0.05 m. Figure 7b, the elevation of the test site ranged from 3550 m to 3649 m, and the largest elevation difference was about 100 m. acquired DTM of the test site was rasterized, resulting in the Digital Elevation Model (DEM) with a spatial resolution of 0.05 m. As shown in Figure 7b, the elevation of the test site ranged from 3550 m to 3649 m, and the largest elevation difference was about 100 m. The DSM of the test site was rasterized (Figure 8a) from Figure 7a. Figure 8b was the rasterized result of filtering out the low vegetation with the height threshold of 3 m. As shown in Figure 8, two rectangular patches within the same location were clipped out of the two DSMs and enlarged, respectively. Patch I in Figure 8b only contained the surfaces of the ground and Minjiang fir trees, while patch I in Figure 8a still had some subalpine shrubs and meadows. According to the comparison of the patch IIs in Figure 8a,b, it could be clearly seen that the ground surface was covered with many subalpine shrubs and meadows before low vegetation removal. The disturbance generated by subalpine shrubs and meadows within the whole test site was almost completely eliminated after classifying and removing low vegetation.

As shown in
The CHM data of the test site (Figure 9a Figure 8b was the rasterized result of filtering out the low vegetation with the height threshold of 3 m. As shown in Figure 8, two rectangular patches within the same location were clipped out of the two DSMs and enlarged, respectively. Patch I in Figure 8b only contained the surfaces of the ground and Minjiang fir trees, while patch I in Figure 8a still had some subalpine shrubs and meadows. According to the comparison of the patch IIs in Figure 8a,b, it could be clearly seen that the ground surface was covered with many subalpine shrubs and meadows before low vegetation removal. The disturbance generated by subalpine shrubs and meadows within the whole test site was almost completely eliminated after classifying and removing low vegetation.

Segmented Tree Canopy Outlines
After ESP processing, the candidate optimal segmentation scales were 29, 45, 102, 112, 127, and 167. When the segmentation scales were set to 29 and 45, the resulting image was too fragmented (Figure 10a,b). When the segmentation scale was set to 102, it was slightly fragmented, but still did not achieve the ideal result (Figure 10c). When the segmentation scale was set to 112, the fragmented objects almost did not exist, and the tree canopy outlines were clearly identifiable ( Figure 10d). As the segmentation scale increased from 127 to 167, segmented objects were gradually merged into larger ones, which was also not an ideal result (Figure 10e,f). Therefore, the optimal segmentation scale for extracting the outlines of individual Minjiang fir tree canopies was 112.

Segmented Tree Canopy Outlines
After ESP processing, the candidate optimal segmentation scales were 29, 45, 102, 112, 127, and 167. When the segmentation scales were set to 29 and 45, the resulting image was too fragmented (Figure 10a,b). When the segmentation scale was set to 102, it was slightly fragmented, but still did not achieve the ideal result (Figure 10c). When the segmentation scale was set to 112, the fragmented objects almost did not exist, and the tree canopy outlines were clearly identifiable ( Figure 10d). As the segmentation scale increased from 127 to 167, segmented objects were gradually merged into larger ones, which was also not an ideal result (Figure 10e,f). Therefore, the optimal segmentation scale for extracting the outlines of individual Minjiang fir tree canopies was 112. The optimal feature set was composed of 15 features including the R-band ratio, the G-band ratio, the B-band ratio, the shape index, the aspect ratio, the R-band mean, the G-band mean, the Bband mean, the G-band mean, the G-band standard deviation, the gray level co-occurrence matrix entropy, the gray level co-occurrence matrix standard deviation, the gray level co-occurrence matrix homogeneity, the gray level co-occurrence matrix mean, the gray level co-occurrence matrix correlation, and the gray level co-occurrence matrix heterogeneity [39]. After object classification using the feature set, the spatial outlines of individual tree canopies were obtained (Figure 9b). The optimal feature set was composed of 15 features including the R-band ratio, the G-band ratio, the B-band ratio, the shape index, the aspect ratio, the R-band mean, the G-band mean, the B-band mean, the G-band mean, the G-band standard deviation, the gray level co-occurrence matrix entropy, the gray level co-occurrence matrix standard deviation, the gray level co-occurrence matrix homogeneity, the gray level co-occurrence matrix mean, the gray level co-occurrence matrix correlation, and the gray level co-occurrence matrix heterogeneity [39]. After object classification using the feature set, the spatial outlines of individual tree canopies were obtained (Figure 9b).

Extracted Individual Tree Heights
Individual tree heights were extracted by overlaying the spatial outlines of Minjiang fir trees in the test site (Figure 9b) with the raster CHM data (Figure 9a). There were about 4873 trees initially identified in the test site, of which 73 (about 1.5%) tree heights equaled 0 m. Removing those false trees, we finally obtained individual heights of Minjiang fir trees in the test site (Figure 11a). The tree height ranged from 3.0 m to 21.69 m. The optimal feature set was composed of 15 features including the R-band ratio, the G-band ratio, the B-band ratio, the shape index, the aspect ratio, the R-band mean, the G-band mean, the Bband mean, the G-band mean, the G-band standard deviation, the gray level co-occurrence matrix entropy, the gray level co-occurrence matrix standard deviation, the gray level co-occurrence matrix homogeneity, the gray level co-occurrence matrix mean, the gray level co-occurrence matrix correlation, and the gray level co-occurrence matrix heterogeneity [39]. After object classification using the feature set, the spatial outlines of individual tree canopies were obtained (Figure 9b).

Extracted Individual Tree Heights
Individual tree heights were extracted by overlaying the spatial outlines of Minjiang fir trees in the test site (Figure 9b) with the raster CHM data (Figure 9a). There were about 4873 trees initially identified in the test site, of which 73 (about 1.5%) tree heights equaled 0 m. Removing those false trees, we finally obtained individual heights of Minjiang fir trees in the test site (Figure 11a). The tree height ranged from 3.0 m to 21.69 m.

Estimated AGB
The AGBs of individual Minjiang fir trees in plots 1-4 (Table 1) were computed with Formula (1). Then, the measured tree heights in Appendix A and the corresponding AGBs in Table 1 were used to fit the power function model (Formula (2)), and the resulting coefficient values of α and β were 0.0391 and 3.4331, respectively ( Figure 12). The R 2 was 0.98, and the RMSE was 25.34 kg. As the R 2 was very close to one and the RMSE was smaller than the average tree AGB of 136.4634 kg, the measured tree heights and the estimated AGBs were well fitted.

Estimated AGB
The AGBs of individual Minjiang fir trees in plots 1-4 (Table 1) were computed with Formula (1). Then, the measured tree heights in Appendix A and the corresponding AGBs in Table 1 were used to fit the power function model (Formula (2)), and the resulting coefficient values of α and β were 0.0391 and 3.4331, respectively ( Figure 12). The R 2 was 0.98, and the RMSE was 25.34 kg. As the R 2 was very close to one and the RMSE was smaller than the average tree AGB of 136.4634 kg, the measured tree heights and the estimated AGBs were well fitted.    The allometric equation for the AGB estimation of individual Minjiang fir trees based only on tree height was obtained as Formula (3). Finally, the AGBs of individual Minjiang fir trees (Figure 11b) in the test site were determined by inputting their extracted tree heights (Figure 11a) into Formula (3). The total estimated AGB of Minjiang fir trees in the test site was about 328.5 t, and the average estimated AGB of individual Minjiang fir trees was about 68.44 kg. 3.4331 (3) where W is the AGB of an individual Minjiang fir tree, and H is the tree height.

Accuracy Assessment
We first assessed the accuracy of individual tree heights extracted from the CHM data. The extracted heights of individual Minjiang fir trees in plots 1-5 are listed in Table 2, and the corresponding measured heights are in Appendix A. As shown in Figure 13a, the R 2 was 0.92 and the RMSE was 1.77 m, which indicated that there was a high correlation between the extracted and measured tree heights in plots 1-5. In addition, the extracted tree heights were generally slightly shorter than the measured tree heights.

Discussion
In the tree AGB estimation of a sparse subalpine coniferous forest, many factors could have influenced the accuracy of the estimation. When acquiring optical photographs, the UAV attitudes and GPS accuracy would affect the georeferencing precision of the resulting data. While generating point cloud data, the image matching and AT processing might introduce some errors [42,43]. When producing DTM, the recommended iterative angle and iterative distance might cause some errors in filtering non-ground points. The three-meter height was chosen as the demarcation line of Minjiang  To assess the accuracy of AGB estimation with the allometric equation with respect to tree height only, the AGBs of individual Minjiang fir trees in plot 5 were calculated with Formulas (1) and (3), respectively ( Table 3).
As shown in Figure 13b, the R 2 was 0.96, and the RMSE was 54.90 kg, suggesting a fair correlation between the tree AGB estimated with the allometric equation on H and DBH, and that only on H. In most cases, the tree AGB estimated only on H was higher than that on H and DBH in plot 5.

Discussion
In the tree AGB estimation of a sparse subalpine coniferous forest, many factors could have influenced the accuracy of the estimation. When acquiring optical photographs, the UAV attitudes and GPS accuracy would affect the georeferencing precision of the resulting data. While generating point cloud data, the image matching and AT processing might introduce some errors [42,43]. When producing DTM, the recommended iterative angle and iterative distance might cause some errors in filtering non-ground points. The three-meter height was chosen as the demarcation line of Minjiang fir trees and low vegetation, but this could have excluded some Minjiang fir trees shorter than three meters. While segmenting individual tree canopies, it was possible to falsely merge multiple trees with overlapping branches as one single tree. The measured heights and DBHs of Minjiang fir trees in plots 1-5, which were regarded as true values, could introduce certain errors due to the nature of measuring instruments and their usage process. As Minjiang fir trees shrink sharply toward the top, it was difficult to ensure that the point cloud generated by the AT method would always exist at the real tops of tree canopies. Additionally, the AGBs of Minjiang fir trees in plots 1-5 might not be accurately estimated, as we directly adopted the allometric equation on H and DBH for Sichuan fir that was provided by the literature. Therefore, the accuracy of extracted individual tree heights and estimated individual tree AGBs depended on the cumulative result of all of the errors incurred throughout the whole process.
Nevertheless, some measures could be taken to further improve the tree AGB estimation accuracy for future applications. For example, the coefficient values of allometric equation on H and DBH could be determined through experiments instead of through the direct adoption of empirical values for Sichuan fir. As shown in Figure 13a, the extracted tree height was generally slightly shorter than the measured tree height, so we could take advantage of the linear relationship to inverse tree heights based on extracted tree heights, and use the inversed tree heights for subsequent AGB estimation. In addition, because the tree DBH cannot be directly obtained from UAV oblique photographs, the allometric equation with tree height as a predictor was used in the study, which may also have affected the AGB estimation accuracy. In future research, the tree AGB could be estimated by combining the tree height and some canopy structure parameters such as projected area and diameter, which could be acquired from UAV oblique photographs.
Accurate individual tree height extraction is the key to AGB estimation for individual trees. As the literature [19,44] shows, the R 2 of extracted tree heights from airborne LiDAR-acquired point cloud data is by far higher than 0.90. Even in some dense seasonal tropical forests, the R 2 is still as high as 0.94 [28]. On the other hand, as the studies [45,46] demonstrate, the R 2 of extracted tree heights from point cloud data generated with the AT method is usually higher than 0.80. In seasonal tropical forests, it may decrease to 0.79 [28]. In our study, due to the relatively low density of the subalpine coniferous forest, the R 2 of the extracted tree heights from overlapping oblique photographs is as high as 0.92, which is very close to the results from LiDAR technology. Compared with the LiDAR-based method for acquiring point cloud data and estimating tree AGB, the proposed AT method based on UAV-acquired oblique photographs is substantially cheaper, faster, and safer, and very suitable for tree AGB estimation of a relatively sparse large-scale forest. However, the 'sparse' or 'dense' forest is a qualitative phrase. The forest density can possibly be defined by the average tree gap in the forest as well as the average tree height and some structure characteristics of tree canopies. The point density [22] or photograph resolution [47] is another important factor in determining the accuracy of tree height extraction. The quantitative correlation between the forest density, the image resolution, and the accuracy of extracted individual tree heights can be further studied and established in future research. Based on the quantitative correlation, we could decide whether or not to adopt this cost-effective method for the tree AGB estimation for a large-scale forest, depending on specific application scenarios and accuracy requirements.

Conclusions
Accurate estimation of the tree AGB of a forest is fairly critical for measuring the carbon sequestration capacity of the forest and studying the carbon cycle of a terrestrial ecosystem. For a large-scale forest, the traditional harvest method is infeasible. The LiDAR is superior in obtaining the point cloud data of dense forest due to its sound penetrating capability, and tree heights can be extracted for AGB estimation. However, the LiDAR is limited by its expensive equipment and complex operations. The terrestrial LiDAR is also not suitable for the AGB estimation of a large-scale forest. Additionally, the airborne LiDAR mounted on a helicopter or a multi-rotor UAV is costly, slow, and risky for acquiring the point cloud data of a large-scale forest in mountainous areas. Instead, the overlapping oblique photographs acquired by a UAV-loaded optical digital camera can also generate point cloud data using the AT method. The digital camera can be mounted on a fixed-wing UAV of high flight speed, long duration, and large spatial coverage, so it can be used for the tree AGB estimation of a large-sale forest. In being restricted by its relatively poor penetrating capacity, the UAV-loaded digital camera is most suitable for obtaining the point cloud data of a relatively sparse forest. In this paper, a patch of sparse subalpine coniferous forest in the source region of the Minjiang River was chosen as the test site. An electric fixed-wing UAV loaded with a consumer-grade digital camera was employed to take oblique photographs of the test site. Based on point cloud data produced from the overlapping photographs, the DSM and DEM of the test site was generated. Then, the AGBs of Minjiang fir trees in the test site were estimated based on individual tree heights extracted from the CHM data and the allometric equation with respect to only tree height. In the accuracy assessment, the R 2 and RMSE of extracted individual tree heights from the CHM were 0.92 and 1.77 m, and the R 2 and RMSE of the estimated AGB of individual trees were 0.96 and 54.90 kg. According to the results of this accuracy assessment, the proposed method for the tree AGB estimation of a sparse subalpine coniferous forest with UAV oblique photography was feasible and effective. However, the quantitative correlation between the forest density, the photograph resolution, and the accuracies of extracted individual tree height and estimated tree AGBs need to be further studied in future research. Acknowledgments: The authors would like to thank Yangchun Wang and Xiaolin Du from Institute of Mountain Hazards and Environment, Chinese Academy of Sciences for help in acquiring oblique photographs with UAV. Furthermore, the authors would like to thank the reviewers for their helpful comments.

Conflicts of Interest:
The authors declare no conflict of interest.