Aboveground Biomass Estimation Using Structure from Motion Approach with Aerial Photographs in a Seasonal Tropical Forest

We investigated the capabilities of a canopy height model (CHM) derived from aerial photographs using the Structure from Motion (SfM) approach to estimate aboveground biomass (AGB) in a tropical forest. Aerial photographs and airborne Light Detection and Ranging (LiDAR) data were simultaneously acquired under leaf-on canopy conditions. A 3D point cloud was generated from aerial photographs using the SfM approach and converted to a digital surface model (DSMP). We also created a DSM from airborne LiDAR data (DSML). From each of DSMP and DSML, we constructed digital terrain models (DTM), which are DTMP and DTML, respectively. We created four CHMs, which were calculated from (1) DSMP and DTMP (CHMPP); (2) DSMP and DTML (CHMPL); (3) DSML and DTMP (CHMLP); and (4) DSML and DTML (CHMLL). Then, we estimated AGB using these CHMs. The model using CHMLL yielded the highest accuracy in four CHMs (R 2 = 0.94) and was comparable to the model using CHMPL (R 2 = 0.93). The model using CHMPP yielded the lowest accuracy (R 2 = 0.79). In conclusion, AGB can be estimated from CHM derived from aerial photographs using the SfM approach in the tropics. However, to accurately estimate AGB, we need a more accurate DTM than the DTM derived from aerial photographs using the SfM approach.


Introduction
Tropical forests have been recognized as an important type of ecosystem that can be used in mitigating climate change, because they sequester and store more carbon than any other vegetation types [1,2].However, tropical forests are being deforested and degraded dramatically through agricultural expansion, wood extraction, infrastructure development, and other natural and anthropogenic processes [3,4].As a result, carbon emissions from deforestation and forest degradation in tropical forests are a major issue in the global carbon budget [5][6][7][8].Reducing emissions from deforestation and forest degradation, and the role of conservation, sustainable management of forests and enhancement of forest carbon stocks in developing countries (REDD+) is one mitigation mechanism related to deforestation and forest degradation in tropical forests.For the implementation of REDD+, a scientifically robust method is required to quantify aboveground biomass (AGB) [9].
Remote sensing with ground-based inventories is expected to play an important role as the method that can be used to quantify AGB for the implementation of REDD+.
Airborne Light Detection and Ranging (LiDAR) is an active remote sensing system which directly estimates the vertical structure of objects by measuring the time of flight between the emitted laser pulses and their received reflectance [10].This system is well suited for measuring forest structural parameters [11][12][13].In the case of the tropics, a canopy height model (CHM) derived from airborne LiDAR can accurately estimate AGB [14,15] and aboveground carbon stocks [16][17][18].While airborne LiDAR is one of the most reliable tools used to estimate AGB in forests, acquisition costs of airborne LiDAR data are often prohibitive and hinder the monitoring of large forest areas [19].Thus, alternative approaches to estimate AGB are required.
Aerial photographs have been common tools used to retrieve forest information for many decades.They are well suited to measure forest cover [20], mean stand height [21], temporal dynamics of AGB [22] and individual tree parameters [23].It has long been demonstrated that tree height can be measured by manual methods of stereo photogrammetry using aerial photographs.We can also create digital surface models (DSMs), which are similar to those derived from airborne LiDAR, from a pair of stereo aerial photographs using stereo matching algorithms.The DSM derived from aerial photographs have significant potential to accurately estimate tree height, stem volume and basal area [24,25].Recent progress in computer science enables the production of DSMs using the Structure from Motion (SfM) approach with a much higher level of automation and much greater ease of use [26].SfM aims to generate 3D geometry from an unordered overlapping collection of photographs using standard automated techniques employing computer vision and photogrammetry [27].Applying the SfM approach enables us to produce a high spatial resolution 3D point cloud model similar to those derived from airborne LiDAR.Previous studies have demonstrated the usefulness of the SfM approach for topographic mapping and landslide monitoring [26,[28][29][30][31].A few studies have evaluated the quality of a canopy height model created by an SfM algorithm [32][33][34].Nevertheless, no known studies evaluate the utility of a CHM created by a SfM algorithm for the estimation of AGB in the tropics.
In general, CHMs are calculated as the difference between the height of the returned pulse or digital surface model (DSM) and ground elevation based on a digital terrain model (DTM).Because the pulses emitted from airborne LiDAR can penetrate below the canopy and reach the ground even in the forested areas, we can generate both a DTM and a DSM from airborne LiDAR alone.Compared with airborne LiDAR, obtaining under-canopy information using aerial photography is difficult.Thus, creating a DTM is the main limitation we face when calculating a CHM from aerial photography.The combined use of aerial photography derived DSM and airborne LiDAR derived DTM provides one solution.Previous studies have demonstrated that a CHM calculated from a DTM derived from airborne LiDAR and DSM derived from aerial photography have a significant potential for the successful monitoring of boreal forests [25,35,36].However, this type of research has not yet been conducted in unmanaged tropical forests.
In this study, we investigated the capabilities of a CHM derived from aerial photography using the SfM approach in estimating AGB in tropical forests.We created four types of CHMs, which were calculated from (1) aerial photography deriving a DSM and aerial photography deriving a DTM (CHMPP) (2) aerial photography deriving a DSM and airborne LiDAR deriving a DTM (CHMPL); (3) airborne LiDAR deriving a DSM and aerial photography deriving a DTM (CHMLP), and (4) airborne LiDAR deriving a DSM and airborne LiDAR deriving a DTM (CHMLL).The results of AGB estimation using an aerial photography-derived CHM are compared with the approach of using an airborne LiDAR-derived CHM for AGB estimation.Finally, we investigated the ability of a DSM derived from aerial photography to estimate AGB.

Study Area
The study area is located in Kampong Thom Province in central Cambodia (Figure 1).Covering a total land area of 1,244,764 ha, this province experiences a typical monsoon Asian climate with a distinct dry season from November to April and a rainy season with about 1700 mm annual precipitation [37].This lowland and nearly flat area has elevations from about 1 to 80 m above sea level.Four forest types occur in the study area: dense evergreen, deciduous, and degraded evergreen forests as well as an area of second-growth forest.While the dominant tree height is between 30m and 40 m in dense evergreen forest, the height is only 10-15 m in deciduous forest.An area of second-growth forest is defined as any area after clearcutting; these areas have shrubs with few residual trees left standing.Degraded evergreen forest is an evergreen forest with evidence of illegal logging based on field surveys.Degraded evergreen includes fewer large trees than evergreen forest, but the landscape still retains its forested nature.

Field Measurements
Field measurements in each forest type were conducted at previously established permanent plots that had been created as a part of a different series of ongoing field studies.The number of plots in dense evergreen, deciduous, degraded evergreen forest and regrowth were 10, 8, 4 and 8, respectively (Table 1).We used three sizes of rectangular permanent plots, which are 2500 m 2 (50 m × 50 m), 1200 m 2 (30 m × 40 m) and 900 m 2 (30 m × 30 m).Eight regrowth plots and four deciduous forest plots cover 2500 m 2 .Five evergreen and two deciduous forests plots cover 1200 m 2 .Five evergreen, four degraded forest and two deciduous forest plots cover 900 m 2 .Field measurements were collected under leaf-on canopy conditions between November 2011 and March 2012.Within each plot, the diameter at breast height (DBH) for all trees with DBH > 5 cm was measured.The coordinates of plot corners were collected using a Global Positioning System (GPS; GPSmap 62s, Garmin, Olathe, KS, USA).In addition, we measured the distance and the azimuth direction between each corner using a laser range finder (Trupulse 360, Laser Technology Inc., Centennial, CO, USA).Because the GPS instrument was not differentially corrected, the accuracy of GPS data was open to question.Thus, we selected the most reliable GPS corner coordinates for each plot through a verification process.Ota et al. (in press) will provide detailed descriptions of the verification process.From the selected corner, we determined the coordinate of other corners mathematically from the distance (i.e., 30 m, 40 m or 50 m) and azimuth directions between the corners.We calculated AGB for each measured tree using general allometric equations [38].AGB of each plot was then calculated by summing AGB of each tree and divided by the plot size.

Remote Sensing Data
Aerial photographs and airborne LiDAR data were simultaneously acquired under leaf-on canopy conditions from a helicopter with an Airborne GPS and inertial measuring unit on 18-21 January 2012.Table 2 shows details of aerial photograph data and airborne LiDAR data specifications.The average flight altitude was 500 m above ground level and the average flight speed was 25 m/s.In the case of aerial photograph data, camera focal length, image size, and pixel size inside the camera were, 51.2499 mm, 8984 × 6732 pixels and 6.0 μm, respectively.In the case of airborne LiDAR, pulse frequency was 100 kHz and average density of first returns-to-sensor was 26 point/m 2 .

Processing of Airborne LiDAR data
From the first and last returns of the airborne LiDAR data, we constructed a 1 m resolution grid of DSM (DSML) and DTM (DTML) (Figure 2).The last return data were used to create a DTM using ground echoes of reflected LiDAR pulses by filtering earlier returns.In the filtering process, the local minima with 10 × 10 m, assumed to represent the ground, were collected.Then, we produced a triangulated irregular network (TIN) that allowed for development of a DTM [39].Finally, the TIN was converted to a regular grid at 1 m resolution to create a DTM.The first return-to-sensor data were also used to create an airborne LiDAR derived DSM.The DSM was created from the highest first return value of pulses for each grid cell.

Processing of Aerial Photographs
We used Photoscan Professional (http://www.agisoft.com;Petersburg, Russia), which is a type of commercial computer vision software, to generate a 3D point cloud from the sets of aerial photographs.Dandois and Ellis [34] and Turner et al. [40] provide a detailed description of the software and workflow.Briefly, the software uses the SfM approach for 3D reconstruction from overlapping collection of photographs.The workflow starts with the "Align Photos" stage, which is the process used to find the camera position and orientation for each aerial photograph and build a sparse point cloud model [41].We selected "High accuracy" and "Ground control pre-selection" as settings.This stage was conducted in the real-world coordinate system, which was Universal Transverse Mercator projection (Zone 48N, WGS 84) based on the camera positions provided by Airborne GPS.We also manually identified ground control points (GCPs) within the aerial photographs to improve the accuracy of the align photos stage.The GCPs consisted of curbs, building corners and so on and were selected from a DSM created by airborne LiDAR data in 2012 (i.e., the DSM created in this study) and aerial orthophotography from 2014 (0.1 m pixel resolution, collected on 20 January 2014, developed by Asia Air Survey Co., Ltd.(Tokyo, Japan).We detected GCP positions in the aerial photographs by manual interpretation of 3D structures and RGB colors corresponding to GCP features identified in the DSM and the orthophotographs following Dandois and Ellis [33].Horizontal and vertical coordinates for GCPs positions were obtained from a DSM derived from airborne LiDAR.Extracting vertical coordinates from airborne LiDAR, we used "Extract Values to Points" function of ArcGIS desktop 10.2.The next workflow step is the "Build dense point cloud" stage which generates a 3D dense point cloud data based on the estimated camera position and orientation for each aerial photograph.We selected "medium quality" and "mild" as settings.Finally, a dense point cloud with 22 points/m 2 , in average, was generated.
From the generated 3D dense point cloud data, we constructed a 1 m resolution grid of DSM (DSMP) and DTM (DTMP) (Figure 2), similar to those derived from airborne LiDAR.We used the same approach when we generated DTML and DTMP.First, we selected points assumed to represent the ground using the local minima with 10 × 10 m.Then, the TIN was then generated from the selected points.The TIN was converted to a regular grid at 1 m resolution to create a DTM.In addition, the DSM was created from the highest value of the point cloud for each grid cell.

Calculation of CHM and CHM-Derived Variables
CHMs were calculated as the relative height between the DSM and DTM.In this study, we created four CHMs, which are the difference between 1) DSMP and DTMP (CHMPP), 2) DSMP and DTML (CHMPL), 3) DSML and DTMP (CHMLP) and 4) DSML and DTML (CHMLL).The pixels of each CHM were classified as either canopy or non-canopy pixels based on their relative height to avoid the potential inclusion of non-tree objects.The pixels of the CHM were classified as canopy if the normalized canopy height was ≥1 m.For each CHM, the variables derived from the CHM were then calculated within each permanent plot, including mean canopy height (the average value of the relative height Hmean), 50th percentile of height (H50), maximum height (H100), and canopy density (D).The canopy density was calculated as the ratio between the number of pixels representing the canopy and the total number of pixels in each plot.

Statistical Analysis
The observed AGB based on field measurements was regressed against the variables derived from the CHM (Equation (1)): where B is the observed AGB (Mg/ha), β0 and β1 and β2 are the regression coefficients, h is the height variable (i.e., Hmean, H50, H100) and d is the canopy density.We used log transformation to simplify Equation (1) as a linear regression in Equation (2): In addition, the influence of forest type (Ftype) on the regression was assessed by expanding Equation ( 2) with dummy variables representing forest type (i.e., evergreen, degraded evergreen, deciduous, and regrowth) following Ota et al. [18].
where bi is the regression coefficient of ith class and zi is the dummy variable of ith class.The dummy variables' names z1, z2 and z3 take on values 1 for degraded evergreen, deciduous and regrowth, respectively and on value 0 for other types of forest in each respective group (z1, z2 and z3).
First, AGB was regressed against height variables and/or canopy density for each of CHMPP, CHMPL, CHMLP, and CHMLL.Additionally, the pertinence of considering forest type in the analysis was assessed.Then, the best model used to estimate AGB using each of four CHMs was selected.
Finally, we compared the accuracy of AGB estimation using the best model of four CHMs.The coefficient of determination (R 2 ), adjusted R 2 and the root mean square error (RMSE) expressed in MG/ha were calculated to express the accuracy of estimates of AGB.Because we used the log transformed equation, R 2 and adjusted R 2 were calculated using the observed and estimated values of AGB.Furthermore, R 2 , adjusted R 2 and RMSE were calculated using leave one-out cross-validation because of the limited number of field plots.

Results
Table 3 shows the results of the regression model fit depicting the relationship between AGB and CHM derived variables when we use CHMPL and CHMPP.When we use CHMPP, the best single-variable model was Hmean (R 2 = 0.31, adjusted R 2 = 0.28, RMSE = 98.82Mg/ha), followed by H100 (R 2 = 0.24, adjusted R 2 = 0.22, RMSE = 103.06Mg/ha).The R 2 , the adjusted R 2 and the RMSE of the combined model of the height variable and canopy density were higher than the corresponding single height model.The best combined model was used Hmean + D (R 2 = 0.41, adjusted R 2 = 0.36, RMSE = 89.54Mg/ha), followed by H100 + D (R 2 = 0.36, adjusted R 2 = 0.31, RMSE = 93.68Mg/ha).Similarly, adding the forest type information improved the accurate estimation of AGB in terms of R 2 , adjusted R 2 and RMSE.The model using D + Ftype yielded the highest adjusted R 2 and the lowest RMSE in all models (R 2 = 0.79, adjusted R 2 = 0.76, RMSE = 51.79Mg/ha).Thus, we conclude that D + Ftype is the best model that can be used to explain AGB independent from forest type when we use CHMPP.
When we used CHMPL, the best single-variable model was Hmean (R 2 = 0.93, adjusted R 2 = 0.93, RMSE = 31.30Mg/ha), followed by H50 (R 2 = 0.92, adjusted R 2 = 0.92, RMSE = 32.55Mg/ha).The RMSE of the combined model of the height variable and canopy density was close to the corresponding single height model.In terms of the adjusted R 2 , canopy density was not informative because the adjusted R 2 of the combined model of the height variable and canopy density was equal to or smaller than the corresponding single height model.Similarly, while adding the forest type information improved the accurate estimation of AGB in terms of RMSE, forest type information is not informative in terms of the adjusted R 2 .The adjusted R 2 of single Hmean was the highest after considering forest type.Thus, we conclude that Hmean could explain AGB independent from forest type when we use CHMPL.Table 4 shows the results of the regression model fit depicting the relationship between AGB and CHM derived variables when we use CHMLP and CHMLL.When we use CHMLP, the best single-variable model was Hmean (R 2 = 0.25, adjusted R 2 = 0.22, RMSE = 104.53Mg/ha).As is the case with AGB estimation using CHMPP, adding the forest type information improved the accurate estimation of AGB in terms of R 2 , adjusted R 2 and RMSE.The model using D + Ftype yielded the highest adjusted R 2 and the lowest RMSE in all models (R 2 = 0.78, adjusted R 2 = 0.75, RMSE = 53.08Mg/ha).Thus, we conclude that D + Ftype is the best model that can be used to explain AGB independent from forest type when we use CHMLP.When we use CHMLL the best single-variable model was Hmean (R 2 = 0.94, adjusted R 2 = 0.94, RMSE = 30.73Mg/ha).As is the case with AGB estimation using CHMPL, canopy density and the forest type information did not improve the accuracy of AGB estimation.Thus, we conclude that Hmean is the best model that can be used to explain AGB independent from forest type when we use CHMLL.

Discussion
When CHMPL and CHMLL were used, the single-variable models using Hmean yielded a higher R 2 and lower RMSE compared with other single-variable models (Tables 3 and 4).In the case of tropical forests, several studies using airborne LiDAR have shown how mean height can be used to estimate AGB or aboveground carbon density in the tropics [15,17,42,43].Our study has confirmed that mean height is particularly a well suited index to estimate AGB, even if we used a DSM derived from aerial photography using the SfM approach.
The R 2 , RMSE and the predicted values of the best model from CHMPL were close to the best model from CHMLL (Tables 3 and 4, Figure 3).These results imply that AGB estimation using metrics from a combination of aerial photography with the SfM approach and airborne LiDAR was comparable to that using metrics from only airborne LiDAR.Previous studies showed that height metrics from aerial photography allow the accurate estimation canopy height in temperate forests [32][33][34], and height, diameter, basal area, stem volume and AGB in boreal forests [36].Our study demonstrated the ability to estimate AGB in tropical forests using the combination of aerial photography with the SfM approach and airborne LiDAR.Acquiring under-canopy information from aerial photography is the difficult task.Thus, difficulty in creating a DTM is the main limitation when we calculate a CHM using only aerial photography.As a result, the estimation accuracy using CHMPP and CHMLP were lower than the estimation accuracy using CHMPL and CHMLL (Tables 3 and 4).These results imply that the accuracy of a DTM derived from aerial photography is not enough to estimate AGB in the tropics and that, for the accurate estimation of AGB, at least a single acquisition of airborne LiDAR is required to create a DTM.However, once a DTM was created by airborne LiDAR, aerial photographs can be used to accurately estimate AGB.This represents a strong advantage in the case of REDD+ implementation because AGB should be estimated repeatedly to monitor deforestation and forest degradation over extended time frames for REDD+ implementation.While this study used a helicopter to acquire the aerial photography, use of a low-cost Unmanned Aerial Vehicle (UAV) is becoming a popular method to collect aerial photographs [34,44].A UAV may contribute to a dramatic cost reduction of repeated monitoring of tropical forests.Thus, we need further research to evaluate the applicability of aerial photographs captured from a UAV, instead of manned aerial vehicle, with SfM approaches in tropical forests.CHMPL accurately estimated AGB without the consideration of forest types (Table 3).Similar results were obtained in previous studies that estimated AGB using airborne LiDAR in boreal coniferous forests [45], temperate coniferous forests [46] and tropical forests [18].Our results confirmed that the combination of aerial photography with the SfM approach and airborne LiDAR can also estimate AGB independently of forest type in tropical forests.These results imply that the calibration process of the equation between AGB and CHM for each individual forest type is not necessary.Only the calibration process for data sets gathering forest types all together is adequate.This may provide a strong advantage for the use of aerial photography and airborne LiDAR in areas where several forest types are intricately distributed and mixed, because it is possible to simplify the calibration process, including laborious and time-consuming forest inventory for each forest type [18].
Although several studies have demonstrated the power of using airborne LiDAR to estimate AGB, the acquisition of airborne LiDAR is costly.Thus, alternative approaches that can be used to estimate AGB are required.In this study, we investigated the capabilities of a DSM derived from aerial photograph using the SfM approach to estimate AGB in tropical forests and compared the results with the more promising airborne LiDAR approach.In conclusion, AGB can be estimated from metrics derived from aerial photographs using the SfM approach in the tropics.However, for the accurate estimation of AGB, we need a more accurate DTM than the DTM derived from aerial photography using SfM approach.Thus, the acquisition of airborne LiDAR that can be used to create a DTM is required.When we calibrate the parameters of the equation between AGB and related metrics, the calibration process for each individual forest type is not necessary.In the future, the applicability of UAV with SfM approaches in tropical forests should be evaluated.on tropical forests and related forest practices.Takio Sano contributed to acquiring remote sensing data in Cambodia.Heng Sokh and Vuthy Ma liaised to obtain permission to conduct research in the Cambodian forest and provided significant knowledge about forests in Cambodia.Eriko Ito, Jumpei Toriyama, Yukako Monda, Hideki Saito, Yoshiyuki Kiyono, Sophal Chann and Nang Ket contributed to conducting field measurements.

Figure 1 .
Figure 1.Study area in Kampong Thom Province, Cambodia with an inset map showing the study site location within extreme Southeast Asia.

Figure 3 Figure 3 .
Figure3shows the observed versus predicted values for AGB estimation using the best model of CHMPP, CHMPL, CHMLP and CHMLL.The predicted values using CHMPL model were close to the predicted values using the CHMLL model, while the predicted values using the CHMPP model were far from the predicted values using the CHM PL and CHMLL models.Both regression lines using the CHMPL and CHMLL model were almost in line with the 1:1 line.

Figure 3 .
Figure 3. Observed aboveground biomass (AGB) versus predicted AGB from the best model of each of the canopy height models (CHMs).(a) CHMPP; (b) CHMPL; (c) CHMLL.A diagonal dotted line and a cross solid line indicate the 1:1 line and regression line between predicted AGB and observed AGB, respectively.

Table 1 .
Summary of field measurements of aboveground biomass.

Table 2 .
Details of aerial photograph data and airborne light detection and ranging (LiDAR) data specifications.

Table 3 .
Results of estimation models each of canopy height model using canopy height model (CHM)PP and CHMPL.

Table 4 .
Results of estimation models each of canopy height model using CHMLP and CHMLL.