Quantifying Live Aboveground Biomass and Forest Disturbance of Mountainous Natural and Plantation Forests in Northern Guangdong , China , Based on Multi-Temporal Landsat , PALSAR and Field Plot Data

Spatially explicit knowledge of aboveground biomass (AGB) in large areas is important for accurate carbon accounting and quantifying the effect of forest disturbance on the terrestrial carbon cycle. We estimated AGB from 1990 to 2011 in northern Guangdong, China, based on a spatially explicit dataset derived from six years of national forest inventory (NFI) plots, Landsat time series imagery (1986–2011) and Advanced Land Observing Satellite (ALOS) Phased Array L-band Synthetic Aperture Radars (PALSAR) 25 m mosaic data (2007–2010). Four types of variables were derived for modeling and assessment. The random forest approach was used to seek the optimal variables for mapping and validation. The root mean square error (RMSE) of plot-level validation was between 6.44 and 39.49 (t/ha), the normalized root-mean-square error (NRMSE) was between 7.49% and 19.01% and mean absolute error (MAE) was between 5.06 and 23.84 t/ha. The highest coefficient of determination R2 of 0.8 and the lowest NRMSE of 7.49% were reported in 2006. A clear increasing trend of mean AGB from the lowest value of 13.58 t/ha to the highest value of 66.25 t/ha was witnessed between 1988 and 2000, while after 2000 there was a fluctuating ascending change, with a peak mean AGB of 67.13 t/ha in 2004. By integrating AGB change with forest disturbance, the trend in disturbance area closely corresponded with the trend in AGB decrease. To determine the driving forces of these changes, the correlation analysis was adopted and exploratory factor analysis (EFA) method was used to find a factor rotation that maximizes this variance and represents the dominant factors of nine climate elements and nine human activities elements affecting the AGB dynamics. Overall, human activities contributed more to short-term AGB dynamics than climate data. Harvesting and human-induced fire in combination with rock desertification and global warming made a strong contribution to AGB changes. This study provides valuable information for the relationships between forest AGB and climate as well as forest disturbance in subtropical zones.


Introduction
Forest biomass is a key biophysical parameter used for evaluating and modeling terrestrial carbon stocks and dynamics and in supporting forest disturbance responses to climate change modeling [1] under rising global temperatures [2].Traditional and frequently used approaches for estimating the spatial distribution of aboveground biomass (AGB) are through field plot measurements [3] or calculation using allometric regression equations or biomass expansion factors [4].Plot measurements that have inadequate spatial coverage and lengthy measurement intervals can limit the effectiveness of field plots in quantifying AGB.Allometric equations use forest structures like diameter at breast height (DBH) to obtain accurate AGB, but such details are only available in rich forest inventory data [5].Thus, both of these techniques are not well-suited for large area AGB spatial distribution measurements when used individually.
The lack of sufficient and high-quality field plots has been identified as a major barrier to develop robust AGB estimates and validation [6,7].The National Forest Inventory (NFI) provides highly detailed information about forest vegetation composition and structure, from which plot-based estimates of forest conditions can be calculated [8].China's NFI is constructed based on 5-year inventory periods including forest type, area, volume, growth, cutting and changes [9,10].
In recent years, remotely sensed data has provided spatially complete prediction information regarding forest cover and change across large areas [11][12][13][14][15]. Various estimation approaches could be used to derive AGB based on field plots' observations, including empirical models ranging from simple linear regression to machine-based modeling [5,14,16].Several studies have shown that the random forest (RF) approach provides the best accuracy among empirical models [17,18].Other physical models, need much more prior knowledge and high complexity and thus are rarely used [5].
Generally, lidar data and high-resolution images can obtain more accurate AGB prediction results [16,[19][20][21][22][23][24]; however, they are limited in spatial and temporal coverage.To date, Landsat images have been the most commonly used in AGB studies as multi-temporal medium-resolution data (30 m), leveraging the complete spatial coverage of imagery [14,17,25], while coarse resolution MODIS, AVHRR, SPOT-VGT are difficult to use because of the mismatch between field measurements and satellite observations [26,27].However, previous studies have commonly used a single image per location [28][29][30] or multiple images but with one in each year in peak growing season to estimate AGB based on empirical models [14,17,31], which would lead to a saturation issue involving underestimation of high AGB values in complex and mature forests.Use of a single year of multi-season NDVI results in more accurate and lower saturation than using a single NDVI [5].Synthetic Aperture Radars (SAR) data, especially ALOS PALSAR, has succeeded in retrieving AGB [32][33][34][35][36][37] (although the saturation effects are still present) [38].The 25 m resolution multi-seasonal mosaic image pixel of this instrument is similar to Landsat pixels and plot size, and the addition of L-band backscatter makes SAR data useful for accurate estimates in relatively homogeneous, young or sparse forests [39].The plantations in northern Guangdong are mostly young forest for timber resources.Therefore, further investigation is required of the combination between multi-date Landsat time-series image compositing [40,41] and PALSAR mosaic data for AGB estimation.Some vegetation indices (such as NDVI and EVI) coupled with field measurements have achieved moderate success for AGB estimation based on satellite imagery and promising results have recently been achieved using texture measurements (such as mean and homogeneity of texture) [17,19].
At present, there is a lack of information regarding AGB under different amounts of disturbance and varied climatic conditions.It is not feasible to periodically detect forest disturbances over large areas through field investigations.Therefore, most studies have used remote sensing data as the main data source for monitoring AGB.Such data includes multisource inventory data, lidar data, land cover classifications, climate and environmental variables are an efficient way to estimate the carbon stocks in forest AGB [6,18,42].Wall-to-wall lidar coverage for large areas is still costly and logistically prohibitive for many regional and NFI programs; however, 500 m resolution MODIS-based biomass linkage to climate and forest disturbance has demonstrated that forest development and climate controls were important contributors to the yearly AGB increase from 2001 to 2010 [43].The Landsat time series (LTS) change detection methods provide pixel-level characterization of forest disturbance and recovery [44][45][46][47] and thus LTS change metrics can improve and overcome some limitations in forest structure [38] and facilitate predictions of forest AGB dynamics over time [14,48,49].Forest disturbance data based on Vegetation Change Tracker (VCT) data can be used to integrate with AGB.Similarly, multiple factors such as social-economic data, rock desertification area and climate data can be used to build a relationship with AGB dynamics.Few studies have explored AGB distribution over large areas forest considering different forest disturbance levels and climate data based on image prediction and NFI data.The main obstacles to such an attempt include the limited annual forest disturbance data [18] and accurate field measurements [6].
In this study, we developed a forest AGB estimation method using a combination of NFI data, time-series passive and active remote sensing techniques , forest disturbance, climate data, topographic data and socio-economic data in Northern Guangdong, China.The methodology comprised several main steps: (a) calculating and examining AGB from national field inventory in 1988, 1992, 1997, 2002, 2007, 2012; (b) four types of explanatory variables (spectral indexes; texture measures and topography data, PALSAR components) combined with Landsat time series data  and PALSAR (2007PALSAR ( -2010) ) were assessed for their effectiveness to capture AGB characteristics; (c) identification of important predictor variables for RF modeling based on multi-temporal imagery, validation through an independent test set and identification of RF mapping quality; (d) integration of AGB changes with forest disturbance maps and (e) multiple exploratory factor analysis to identify potential factors of forest AGB dynamics.

Study Area
The study area (black region covering 24,275.5 km 2 ) is located in northern Guangdong Province, with a geographical coverage extending from 113.10 ˝E, 23.64 ˝N to 114.75 ˝E, 25.44 ˝N, including the cities of Shaoguan, Qingyuan and Heyuan (Figure 1).The local topography is undulating and its elevation is between 22 and 1353 m above sea level.The rocky desertification process is very typical in northern Guangdong, representing one of the serious rocky desertification areas in China.The climate is a mid-subtropical monsoon climate, with 1300 to 2400 mm of mean annual precipitation and of 18 to 21 ˝C mean annual temperature.The rainy season takes place from March to August, with approximately 53% of the annual rain falling between April and June.The vegetation includes natural forests and the large scale plantations.The tree cover is dominated by Pinus massoniana, Cunninghamia lanceolata, Pinuselliottii engelm, Eucalyptus, Pinus kwangtungensis, Castanopsis fissa, Acacia mangium and Phyllostachys edulis.Most of these species are considered evergreen and fast-growing.Other minor vegetation consists of deciduous trees and shrubs.The most common meteorological disasters in the region are chilling injury of plants, storms and flooding, and drought.

Field Data
The NFI in China is the first level of China's three-tiered inventory system, which, built upon permanent sample plots (PSPs) created from systematic sampling and statistical induction processes, is administered by the State Forestry Administration [10].The NFI has been carried out eight times from the 1970s to 2012 in China.The Guangdong Provincial Center for Forest Resources Monitoring has provided eight NFI datasets collected between 1979 and 2012.In this study, we used six years (1988,1992,1997,2002,2007,2012) of data from inventory plots located within the study area (Table 1).A total of 1355 plots with a size of 25.82 m ˆ25.82 m were distributed and established within the study area.Plot level AGB was derived from dominant trees biomass using tree species-specific allometric equations developed by the Guangdong Provincial Center for Forest Resources Monitoring.In cases where specific tree equations were not available, the allometric equations for a similar tree species were used (Appendix Table A1).2).The study area of Northern Guangdong is covered by the Landsat world reference system-2 path/row 122/043 (Figure 1).The best quality, cloud-free images as close to peak growing season were acquired to minimize the effects of vegetation phenology, sun angle differences and other exogenous factors.However, cloud contamination images (≤50% cloud cover) remained.Multiple images were composited to remove cloud and phenology effects.The image data were radiometrically and atmospherically calibrated to surface reflectance based on the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) algorithm [50].C correction in combination with DEM (elevation), solar elevation, solar azimuth and topographic kernel (based on the image resolution) was conducted in ENVI classic (topo correction express tool) to correct terrain effects.Iteratively reweighted multivariate alteration detection (IR-MAD) [51] was used to radiometrically normalize each image based on the high-quality reference image in 2004.2).The study area of Northern Guangdong is covered by the Landsat world reference system-2 path/row 122/043 (Figure 1).The best quality, cloud-free images as close to peak growing season were acquired to minimize the effects of vegetation phenology, sun angle differences and other exogenous factors.However, cloud contamination images (ď50% cloud cover) remained.Multiple images were composited to remove cloud and phenology effects.The image data were radiometrically and atmospherically calibrated to surface reflectance based on the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) algorithm [50].C correction in combination with DEM (elevation), solar elevation, solar azimuth and topographic kernel (based on the image resolution) was conducted in ENVI classic (topo correction express tool) to correct terrain effects.Iteratively reweighted multivariate alteration detection (IR-MAD) [51] was used to radiometrically normalize each image based on the high-quality reference image in 2004.1990206, 1993278, 1995284, 1996159, 1996191, 2001252, 2003290, 2004277, 2005199, 2006266, 2007205, 2008208, 2009290, 2010213, 2011232 L7 ETM+ 1999255, 1999287, 2000258, 2001260, 2002311 BJGS-China L5 TM 1986307, 1988313, 1992212, 1994313, 1997305 2).A tile product consists of two bands in HV and HH polarizations at 25 m spatial spacing, geometrically and radiometrically corrected and normalized for topography.The PALSAR data were processed to convert digital number (DN) to backscatter coefficient (σ ˝) in decibel (dB) for the HH and HV polarization components based on the following equation: where, CF is the calibration factor and its value for FBD 343 HH is ´83.2 and HV is ´80.2 before 9 January 2009, and ´83.0 is valid for data processed by JAXA after 9 January 2009 [53].The PALSAR data were co-registered with Landsat data.To reduce speckle effects, we later applied the enhanced Lee filter using a window size of 5 ˆ5 pixels [54].The variables derived incorporated both HH and HV polarizations and HH/HV (Figure 1).

Plot-Level Explanatory Variables
A large suite of available predictor variables (spectral, topography, texture, PALSAR backscatter coefficient) were used to predict AGB (Table 3).The spectral variables (spectral indices, tasseled cap transformations) were derived from six years of Landsat surface reflectance imagery coinciding with NFI acquisition years.The NFI plots in 2007 and 2012 were used for Landsat imagery from 2006 to 2011, respectively, as there was mass cloud contamination in 2007 and no available data for 2012.Eight texture variables (Table 3), chosen to represent spatial features in the imagery, were derived from the six Landsat bands.Texture variables were calculated with three window sizes of 3 ˆ3, 5 ˆ5, and 7 ˆ7 pixels at an offset ( [1,1]), and a 64 gray level quantization.Topography variables from the NASA Shuttle Radar Topographic Mission (STRM) were re-projected from 90 m resolution to 30 m spatial resolution, including elevation, slope and topographic solar radiation index.To examine the role of polarization in estimating AGB, PALSAR backscatter coefficient for HH, HV and HH/HV were derived.mean, variance, homogeneity, contrast, dissimilarity, entropy, second moment, correlation [68] GLCM texture measures PALSAR HH, HV, HH/HV PALSAR components

RF Modeling and Implementation
The RF model is a non-parametric algorithm with the following advantages: capacity to determine variable importance, robust to data reduction, not over fit, producing unbiased accuracy estimate and higher accuracy than decision trees and low sensitivity to tuning of the parameters.The weaknesses of the RF method are that the decision rules are unknown (black box), it is computationally intense and it requires input parameters [69,70].We used RF regression tree implemented through the R package ModelMap [71] by calling the R package random forest [72] to identify important predictor variables, to model the relationship between the predicted variables and AGB, and to apply the model over the study area for mapping AGB based on single or composited Landsat surface reflectance images and mosaic PALSAR data.There were 500 trees (ntree) used for modeling, and each tree was built from a four-fifths random sampling of the training data.For the parameter mtry, we used the default value of the square root of the total number of predictor variables.The parameter nodesize was set to the default value of 1.The result is an ensemble of low bias and high variance regression trees, where the final predictions are derived by averaging the predictions of the individual trees [69].We preformed RF modeling using predictor variables derived from multi-temporal imagery.The AGB time series trends were derived by extrapolating to the scene-level.

Variable Selection
All training data were used for internal model performance evaluation and for deriving two different variable importance measures (VIM) based on the response type.Importance type 1 was calculated by randomly permuting each predictor variable and computing the associated reduction in predictive performance using the 'out of bag' (OOB) error for RF models [69] and importance type 2 was calculated using the decrease in node impurities attributable to each predictor variable.The resulting VIM provides means to assess the contribution of each predictor variable to the modeling performance.The higher the percent increase in mean square error (MSE) (PercentIncMSE) and increase in NodePurity (IncNodePurity) indicated a stronger importance of these predictor variables.Mean squared prediction error (MSE) VIM was used to rank all predictor variables [17].The initial variable ranking (all predictor variables) is used throughout all the iterations.The less important predictor variables were then removed from the model and the small subset of predictor variables was selected to construct the final model and map AGB based on one RF model containing all plot measurements from six years.

Accuracy Assessment and Validation
The RF modeling was conducted using the reduced predictor variable datasets for multi-temporal Landsat imagery and PALSAR.The predictive ability of the model was assessed using an independent dataset (1/5), with OOB predictions on the training data.The independent validation datasets approach is based on reference data in Table 4. Four measures of model performance were calculated from the independent test set, including the coefficient of determination (R 2 ), the root mean square error (RMSE), the normalized root-mean-square error (NRMSE) and mean absolute error (MAE).In this study, the RMSE and MAE were both used to evaluate models by summarizing the differences between the observed and predicted values.The MAE gives equal weight to all errors, while RMSE gives extra weight to large errors.The NRMSE is often expressed as a percentage, where lower values indicate less residual variance.In this study, we validated the linear fit relationships between the predicted data based on RF AGB mapping results in 1988,1992,1997,2002,2006   Maps of the stochastic uncertainty remaining in the final RF model were created to determine the mapping quality.As the mean AGB increases, so does the standard deviation [71].High values of this standard deviation indicate a lack of agreement between the modeled trees.By calculating the standard deviation and the coefficient of variation (CV) for each pixel from the predictions of each of the independent randomly generated trees that compose the RF model, we estimated the RF mapping quality.The uncertainty is considered highest when the modeled trees were not in agreement, with some trees predicting low AGB and others predicting high AGB.

Integration of AGB Changes with Forest Disturbance Maps
Landsat time series stack (LTSS) ranging from 1986 to 2011 was used to input to the VCT algorithm to produce the disturbance components (e.g., disturbance year, annual disturbance maps, disturbance magnitude).Further information about the algorithm is available in Huang et al. [45].Thomas et al. 2011 indicated the accuracy of over 92% for change versus no-change independent based on validation of VCT maps across several of the sample scenes [73].We estimated the time since forest disturbance and quantified AGB decrease depending on the VCT disturbance pixels.The AGB time series stack was intersected with the annual map of forest disturbance derived from the VCT algorithm applied to the LTSS.The non-forest part of the predicted AGB map was masked based on the VCT annual non-forest type distribution map.The final annual AGB decrease map connected with forest disturbance map was clipped by the study area boundary.

Multiple Exploratory Factor Analysis about Potential Factors of Forest AGB Dynamics
The AGB decrease was strongly connected to the change in forest disturbance in northern Guangdong as a result of various factors including human activities, natural forest disturbance and climate change.Figure 2 records the time and cause of disturbances occurring in northern Guangdong.To determine the missing detailed statistical data of different forest disturbance types, VCT-based disturbance (which proven to be sensitive to harvests and fire [74] could be used as a proxy for these data) in northern Guangdong, because harvesting and slight human-induced fire were mainly disturbance types based on the national forestry yearbook, specifically clear cutting and prescribed fire of intentional clear cut (Figure 2b).The carbonatite rock topography was dominated by rocky desertification, northern Guangdong has the most area of the whole rocky desertification in Guangdong.The area of rocky desertification from 1974 to 2004 (based on Li et al. [75]) represented a clear linear decreasing trend.We calculated the approximate rock desertification area from 2005 to 2011 based on this linear relationship, constructing a dataset of yearly continuous rock desertification area change (Figure 2a).Furthermore, other human activity factors including population growth, industrial production, agricultural production, mining production, forestry production, per capita GDP and highway length were collected from the Statistics Bureau of Guangdong Province.

Multiple Exploratory Factor Analysis about Potential Factors of Forest AGB Dynamics
The AGB decrease was strongly connected to the change in forest disturbance in northern Guangdong as a result of various factors including human activities, natural forest disturbance and climate change.Figure 2 records the time and cause of disturbances occurring in northern Guangdong.To determine the missing detailed statistical data of different forest disturbance types, VCT-based disturbance (which proven to be sensitive to harvests and fire [74] could be used as a proxy for these data) in northern Guangdong, because harvesting and slight human-induced fire were mainly disturbance types based on the national forestry yearbook, specifically clear cutting and prescribed fire of intentional clear cut (Figure 2b).The carbonatite rock topography was dominated by rocky desertification, northern Guangdong has the most area of the whole rocky desertification in Guangdong.The area of rocky desertification from 1974 to 2004 (based on Li et al. [75]) represented a clear linear decreasing trend.We calculated the approximate rock desertification area from 2005 to 2011 based on this linear relationship, constructing a dataset of yearly continuous rock desertification area change (Figure 2a).Furthermore, other human activity factors including population growth, industrial production, agricultural production, mining production, forestry production, per capita GDP and highway length were collected from the Statistics Bureau of Guangdong Province.
Climate data were obtained from the China Meteorological Administration International Exchange Station Ground dataset.We used the annual mean temperature, annual mean minimum temperature, annual mean maximum temperature, extreme lowest temperature, extreme highest temperature, precipitation, daily maximum precipitation, mean relative humidity and minimum relative humidity for Shaoguan and Heyuan station in northern Guangdong Province from 1980 to 2013.In order to evaluate the impacts of natural and human disturbance on AGB, we first applied the z-score normalization method (still including all information of the original) based on original data's mean and standard deviation to eliminate the effect of these factors dimensions and data quantity.Next, the correlation among these factors (dimensionless, based on Spearman correlation) was Climate data were obtained from the China Meteorological Administration International Exchange Station Ground dataset.We used the annual mean temperature, annual mean minimum temperature, annual mean maximum temperature, extreme lowest temperature, extreme highest temperature, precipitation, daily maximum precipitation, mean relative humidity and minimum relative humidity for Shaoguan and Heyuan station in northern Guangdong Province from 1980 to 2013.
In order to evaluate the impacts of natural and human disturbance on AGB, we first applied the z-score normalization method (still including all information of the original) based on original data's mean and standard deviation to eliminate the effect of these factors dimensions and data quantity.Next, the correlation among these factors (dimensionless, based on Spearman correlation) was analyzed.We then adopted exploratory factor analysis based on the varimax orthogonal rotation method to calculate the principal component factor loading matrix of nine climate and nine human activities elements changing values.Varimax rotation is the most commonly used of the rotations that are available.This first involves scaling the loadings and maximizing the factor loadings, assuming no correlations between components [76].We scaled the loadings by dividing them by the corresponding communality as shown below: where, the loading of the ith variable on the jth factor after rotation, where ĥi is the communality for variable.What we want to do is to find the rotation which maximizes this quantity.The Varimax procedure, as defined below, selects the rotation to find this maximum quantity: which represents the sample variances of the standardized loadings for each factor, summed over the m factors.The final objective is then to find a factor rotation that maximizes this variance and represents the dominant factors affecting the AGB dynamics.

Variable Importance and Selection
The PercentIncMSE and IncNodePurity estimated from RF OOB data were used to rank all of the predictor variables by their capacity to predict AGB. Figure 3 shows the entire variable importance ranking for multi-temporal imagery.The progressive removal of the least important predictor variables generally resulted in reduced RMSE for the OOB data and the model with the lowest RMSE (2 t/pixel) and 13 predictor variables were selected for mapping AGB.The reduced variables(Figure 4) that are important for the prediction of AGB according to the ranking include Landsat band 1 (B1), EVI, TCA, DEM, Landsat band 2 mean texture feature (mean23: window size 3 ˆ3 mean, mean25: 5 ˆ5 mean, mean27: 7 ˆ7 mean), band 5 mean texture feature (mean57: 7 ˆ7 mean), band 6 mean texture feature (mean65: 5 ˆ5 mean, mean 67: 7 ˆ7 mean) (Table 3), and PALSAR HH, HV, HH/HV.The most important predictor was mean 67, this variable was highly correlated with several of the other variables (between-texture mean correlation up to 0.8) and moderately correlated with TCA of 0.55.Ultimately, ten variables were used for mapping AGB before 2006 because they did not include HH, HV and HH/HV.

Predictive Performance of the RF Regression Models
The results from independent test set with OOB predictions showed that the use of thirteen predictor variables and multi-temporal imagery based on 500 trees produces an accurate prediction of AGB.One pixel per plot for different years was included in the validation dataset, resulting in 100-200 reference pixels for AGB (Table 4).The range and detailed statistics of the observed and predicted AGB are shown in Table 4.The R 2 values for the linear fit of predicted variables and observed variables for 1988,1992,1997,2002,2006 and 2011 were all higher than 0.50.Plot-level RMSE values were between 6.44 and 39.49 t/ha, the NRMSE was between 7.49% and 19.01%and the MAE values for modeling AGB were between 5.06 and 23.84 t/ha.The highest R 2 of 0.8 and the lowest NRMSE of 7.49% were reported in 2006.The lowest RMSE of 6.44 t/ha, MAE of 5.06 t/ha and the second highest R 2 value of 0.71 occurred in 1988.Overall, although there were some prediction errors, the most accurate model for predicting and mapping AGB may be determined based on these reduction variables (Figure 5).
the predictor variables by their capacity to predict AGB. Figure 3 shows the entire variable importance ranking for multi-temporal imagery.The progressive removal of the least important predictor variables generally resulted in reduced RMSE for the OOB data and the model with the lowest RMSE (2 t/pixel) and 13 predictor variables were selected for mapping AGB.The reduced variables(Figure 4) that are important for the prediction of AGB according to the ranking include Landsat band 1 (B1), EVI, TCA, DEM, Landsat band 2 mean texture feature (mean23: window size 3 × 3 mean, mean25: 5 × 5 mean, mean27: 7 × 7 mean), band 5 mean texture feature (mean57: 7 × 7 mean), band 6 mean texture feature (mean65: 5 × 5 mean, mean 67: 7 × 7 mean) (Table 3), and PALSAR HH, HV, HH/HV.The most important predictor was mean 67, this variable was highly correlated with several of the other variables (between-texture mean correlation up to 0.8) and moderately correlated with TCA of 0.55.Ultimately, ten variables were used for mapping AGB before 2006 because they did not include HH, HV and HH/HV.

Predictive Performance of the RF Regression Models
The results from independent test set with OOB predictions showed that the use of thirteen predictor variables and multi-temporal imagery based on 500 trees produces an accurate prediction of AGB.One pixel per plot for different years was included in the validation dataset, resulting in 100-200 reference pixels for AGB (Table 4).The range and detailed statistics of the observed and predicted AGB are shown in Table 4.The R 2 values for the linear fit of predicted variables and observed variables for 1988,1992,1997,2002,2006 and 2011 were all higher than 0.50.Plot-level RMSE values were between 6.44 and 39.49 t/ha, the NRMSE was between 7.49% and 19.01%and the MAE values for modeling AGB were between 5.06 and 23.84 t/ha.The highest R 2 of 0.8 and the lowest NRMSE of 7.49% were reported in 2006.The lowest RMSE of 6.44 t/ha, MAE of 5.06 t/ha and the second highest R 2 value of 0.71 occurred in 1988.Overall, although there were some prediction errors, the most accurate model for predicting and mapping AGB may be determined based on these reduction  variables (Figure 5).
Figure 6 illustrates the Google Earth aerial photo with RF uncertainty regarding the mean, standard deviation and CV in 2011 for regions of northern Guangdong.When the CV value was high, the mean of AGB was low.The map of RF uncertainty highlighted small anomalous areas of high uncertainty.These are localized areas where the predictions from individual trees varied widely, with some of the trees in the RF predicting low AGB and other trees predicting high AGB.

Forest AGB Dynamics across Northern Guangdong
Based on the RF regression modeling results, the scene-level mean AGB time series change trends were derived (Figure 7).The modeled trajectories supported the ascending trend found by the NFI data between 1988 and 2011.In particular, the modeled mean AGB in NFI available years closely approximated the NFI mean AGB data.A strong rising trend of mean AGB from the lowest value of 13.58 t/ha to the highest value of 66.25 t/ha occurred between 1988 and 2000, while after 2000 there was a fluctuating ascending change.The peak value of mean AGB was 67.13 t/ha in 2004.

Forest AGB Dynamics across Northern Guangdong
Based on the RF regression modeling results, the scene-level mean AGB time series change trends were derived (Figure 7).The modeled trajectories supported the ascending trend found by the NFI data between 1988 and 2011.In particular, the modeled mean AGB in NFI available years closely approximated the NFI mean AGB data.A strong rising trend of mean AGB from the lowest value of 13.58 t/ha to the highest value of 66.25 t/ha occurred between 1988 and 2000, while after 2000 there was a fluctuating ascending change.The peak value of mean AGB was 67.13 t/ha in 2004.
By integrating the AGB change with forest disturbance, the trend in AGB decrease associated with forest disturbance is shown in Figure 8.The trend in disturbance area closely corresponds with the trend in AGB decrease, except for a couple of notable anomalies.Both the disturbance area and AGB decrease spiked in 1994, and a significant spike in AGB decrease occurred in 2008 with a much smaller associated spike in the disturbance area.This divergence presented a clear disagreement with the entire change trend.Although overall ascending trend of mean biomass was reported, forest disturbance deduced AGB decrease annually.

Regional Climate Change
To investigate climate change over northern Guangdong Province, we analyzed trends in the average values of temperature and precipitation of two stations from 1980 to 2013 in this region.During this period, an increasing linear trend in mean annual temperature was found, while a decreasing linear trend in annual precipitation was observed (Figure 9).The change trend predicted the increasing temperature and decreasing precipitation that occurred in this region.During this period, an increasing linear trend in mean annual temperature was found, while a decreasing linear trend in annual precipitation was observed (Figure 9).The change trend predicted the increasing temperature and decreasing precipitation that occurred in this region.
Correlation coefficient analysis of the following nine climate datasets from 1990 to 2011 was conducted in relation to the AGB decrease (Y) in the same trend: mean temperature (X1), mean maximum temperature (X2), mean minimum temperature (X3), annual precipitation (X4), extreme low temperature (X5), extreme high temperature (X6), mean humidity (X7), min humidity (X8) and maximum daily precipitation (X9).The correlation coefficients were showed in Table 5.Furthermore, the correlation coefficients between forest disturbance and mean temperature, mean minimum temperature, extreme low temperature, minimum humidity in relation to forest disturbance were −0.3, −0.33, −0.311 and −0.271, respectively.This showed that the above factors were slightly associated with AGB decrease, and minimum humidity and low temperature had a more notable influence on the AGB dynamics than other factors.Correlation coefficient analysis of the following nine climate datasets from 1990 to 2011 was conducted in relation to the AGB decrease (Y) in the same trend: mean temperature (X1), mean maximum temperature (X2), mean minimum temperature (X3), annual precipitation (X4), extreme low temperature (X5), extreme high temperature (X6), mean humidity (X7), min humidity (X8) and maximum daily precipitation (X9).The correlation coefficients were showed in Table 5.Furthermore, the correlation coefficients between forest disturbance and mean temperature, mean minimum temperature, extreme low temperature, minimum humidity in relation to forest disturbance were ´0.3, ´0.33, ´0.311 and ´0.271, respectively.This showed that the above factors were slightly associated with AGB decrease, and minimum humidity and low temperature had a more notable influence on the AGB dynamics than other factors.

Human Activities
We conducted correlation coefficient analysis of the following nine human activities from 1990 to 2011 in relation to the biomass decrease (Y) in the same period: forest disturbance (logging and fire) (X10), population (X11), industrial production (X12), agricultural production (X13), mining (X14), forestry production (X15), per capita GDP (X16), highway mileage (X17) and rock desertification area (X18) (Figure 10).The correlation coefficients were found in Table 6 (significant correlation at the 0.01 level).This showed that the above factors were significantly directly and indirectly associated with AGB decrease, and the annual forest disturbance and mining changes had a great influence on the AGB dynamics.The rotated component matrix (Table 7) showed the factor loadings for each variable.In Factor 1 (F1), the human activities variables X18, X17, X11, X13, X16, X12, X14 and X15 and climate variables X5 and X8 were loaded strongly.The rock desertification area (X18) was the strongest variable.Most strong loading variables were associated with human activities and the human-induced activities dominated Factor 1.In Factor 2 (F2), climate variables X1, X3, X2 and human activities variables X10 and X14 were strongly loaded.The climate variables performed better than human activities in Factor 2, especially temperature.In Factor 3 (F3), climate variable X4 was strongly loaded and variables X9,  The rotated component matrix (Table 7) showed the factor loadings for each variable.In Factor 1 (F1), the human activities variables X18, X17, X11, X13, X16, X12, X14 and X15 and climate variables X5 and X8 were loaded strongly.The rock desertification area (X18) was the strongest variable.Most strong loading variables were associated with human activities and the human-induced activities dominated Factor 1.In Factor 2 (F2), climate variables X1, X3, X2 and human activities variables X10 and X14 were strongly loaded.The climate variables performed better than human activities in Factor 2, especially temperature.In Factor 3 (F3), climate variable X4 was strongly loaded and variables X9, X7 and X6 were fairly strongly loaded, especially precipitation and humidity.The variance explanation of rotation sums of squared loadings showed Factor 1, 2 and, 3 account for 45.246%, 23.010% and 12.488% of the variability across all 18 variables, respectively.This indicates that Factor 1 explained much more variability in the AGB data than other factors, suggesting that human activities are more associated with AGB decrease than climate data.Six sets of NFI field data (1988 to 2012) covering northern Guangdong were derived to predict AGB based on remote sensing.We selected NFI data and filtered by removing non-forest plots (i.e., water and urban plots) and economic species (fruit trees, tea tree, mulberry), cloud contamination plots.We retained most plantations tree species and natural tree species such as masson pine, eucalyptus, Chinese fir and Pinus elliottii to construct a 'homogeneous' region.Six years of continuous recording field data, encompassing the sufficient and high quality plot-level AGB information, and the permanent location of field plots ensure an accurate geo-location, which was extremely important to link with spectral pixels.Su et al. collected over 8000 ground inventory records from published literature and stated that the huge plot location uncertainty can result in significantly different correspondence between field-measured AGB and predictors [6].To develop AGB, equations of different dominant trees per plot were used, which were developed by the local forestry agency.We need to determine their accuracy and use them to derive AGB; however, a bias originated from tree species' allometric equations, which were replaced by the similar tree species.These equations would also provide a reference for the same tree species to predict AGB in sub-tropical regions.

Image-Based Predicted Variables and Other Ancillary Data
Four types of variables were prepared as input parameters: spectral indices, topography variables, texture measures and PALSAR HH, HV, HH/HV.The best predictors were chosen based on the variable importance estimation.The relatively strong relationship between image texture, in particular the gray level co-occurrence matrix mean, and AGB agreed with previous research suggesting image texture variables have a stronger correlation with observed AGB than models using only physical and spectral variables [77].Although the single structure and regular distribution plantations are popular in northern Guangdong, the complex forests are co-existence among natural forests and the plantations.The relationship between image texture metrics and measurements of forest attributes can be used to help characterize complex forests, and enhance fine vegetation biophysical properties, a difficult challenge when using spectral vegetation indices especially in closed canopies [78].It is clear that image textural measures have the potential to provide an attractive opportunity for monitoring AGB [29].Texture metrics obtained from single Landsat 8 has also proved the plausible performance and strength of texture metrics in improving AGB estimates [29].Kelsey and Neff [77] stated that texture analysis was efficient in addressing saturation problems associated with vegetation indices when mapping biomass especially in dense canopies as it correlated more with AGB than spectral parameters.The L-band backscatter was strongly correlated with forest AGB apart from dense forests (up to 300-400 t/ha) AGB saturation [79], which can improve the predicted capability (Figure 5e,f).Some research have proved that the backscatter on HV polarisation gave better field-based AGB estimation than HH backscatter [80] and their derivates [52].Günlü et al. 2014 [81] reported that using Enhanced Vegetation Index (EVI)to estimate AGB presented better estimation than individual band reflectance values, and less susceptible to the saturation problem in temperate and tropical forests [82].Topographic information (DEM) was also included to analyze its influence on the biomass estimation in mountainous region.The STRM DEM data with 90 m spatial resolution was used, consequently, this issue requests future studies using a high resolution DEM (ASTER GDEM, ALOS-PRISM).

RF Regression and Validation
The RF regression approach has several advantages for modeling remote sensing data [69], but also is associated with some limitations.The R 2 between multi-temporal Landsat, PALSAR and the reference data reached 0.80 (RMSE = 24.02t/ha, NRMSE = 7.49%) in 2006 and at least R 2 of 0.5 in other years for AGB.Karlson et al. [17] reported that R 2 based on seven Landsat 8 images from 2013 to 2014 reached 0.57 (RMSE = 17.6 t/ha, NRMSE = 66%) for AGB.Powell et al. used RF based on multi-temporal Landsat images and reported an RMSE of 39.23Mg/ha of observed and predicted AGB in Minnesota and 32.19 Mg/ha in Arizona [14].The validation results of our study are better than other reported studies of modeling AGB using optical images with a combination of radar image.Topography correction and LEDPAS atmospheric correction mitigated the impacts of cloud shadow, mountain shadow and slight haze.C correction has been reported to be capable of reducing the topographic effects caused by the topographic relief [83].We found that RF modeling overestimated the low values and underestimated the high values (Table 4, Figure 5), which partly explained the absence of bias in the AGB prediction, indicating a saturation problem.This effect was most pronounced for AGB predictions and was due to both properties of the algorithm and characteristics of the reference data.In particular, the reference data need to cover the full range and represent the variability of the variable of interest in the specific study area.In general, the results of the RF method indicate that time-series derived from multi-temporal Landsat images and PALSAR data can improve the accuracy of AGB estimation and reduce the saturation problem.
In the RF uncertainty AGB maps (Figure 6), it is interesting to look at the predictions and uncertainty.Either particular predictors or particular training locations contributed to the high variation in the estimates of the response variable.The RF standard deviation map showed that although the mean prediction from the 500 trees in the model was for moderate AGB, there was a very high level of uncertainty for upper left (Figure 6) water areas.Because of the structure of the RF model, some of the trees would be constructed from subsets of the predictors containing neither NDWI nor MNDWI, leading to higher levels of between-tree RF uncertainty, while the right areas (Figure 6) were correctly identified as low AGB with low values of RF uncertainty.

Forest AGB and Forest Disturbance Dynamics Change
Extrapolating to the scene-level, we examined time-series AGB trends.The 'Green Guangdong' afforestation plan was established in the 1980s and resulted in mountain forests regrowth in the early 1990s.In 1994, mean AGB experienced a low value as a result of the implementation of the eucalyptus construction plan: a large number of forest was clear cut for planting fast-growing eucalyptus.After 2000, the permission for a low slope development plan and mountainous urbanization increase could explain the decreasing AGB trend.To protect water source forests and conservation forests, the eucalyptus plan was banned by local government, but some plantations remain.Furthermore, afforestation on the barren mountains, cutting-blank, fire-slash, inefficient forest lands was conducted, which together led to a slight AGB increase.The Chinese government carried out the first phase of a rocky desertification control project from 2006, which was beneficial for reducing cutting and protecting vegetation, water and soil [84].The overall rocky desertification area in Southwest China (containing northern Guangdong) has been reduced 7.4% from 2005 to 2011 [85].Figure 8 depicts a significant increase of AGB in 2008, which was likely to be the result of a VCT error in that time period associated with clouds providing a false-positive change detection [73].Stone (2008) estimated that freezing rain and snow disasters in 2008 caused the loss of one-tenth of China's forests and plantations, which was roughly equivalent to the number of hectares that were reforested between 2003 and 2006, according to China's State Forestry Administration [86].This highlights the importance of deriving spatially explicit maps of AGB dynamics within disturbed regions, because of the regional conditions of forest and disturbance.

Correlation Analysis
Although we found increasing temperature and decreasing precipitation, the correlation between climate and AGB decrease was not high compared with human-induced activities.Actually, climate effects were not a temporary and separate process, as many factors co-exist with human activities.Zhang et al. (2014) found that climate explained relatively little of the observed, stand-level variation in Alberta forest AGB [18], which is consistent with the finding of Stegen et al. [87] of AGB-climate relationships in temperate and tropical forests.
Harvesting (Figure 2b) and human-induced fire were significantly correlated with AGB, as the fast-growing and high-yield plantations have short rotation cycle and need to be cut for a new regeneration.The photos shown in Figure 2b indicate the harvest change from 2009 to 2014 and showed the fast cutting and regrowth during only 2 to 4 years.This was because the large-scale afforestation and reforestation efforts made eucalyptus plantations expand rapidly [88]; however, these afforested areas are primarily managed as successive short-rotation (around 6-8 years) plantations [89].This supports the considerable human disturbance leading to the AGB change.Asia Pulp and Paper (APP) planted a total plantation resource of 330,000 ha in Qingyuan and 400,000 ha in Shaoguan during 1995 to 1998, while prescribed burning is widely used for seeding or planting after harvesting, which is a common practice for forest plantations in southern China, especially for Chinese fir plantations [90].

Quantification and Qualitative Analysis
Rocky desertification in our study (Figure 2) occurs with karst topography as a result of human activities (illogical or excessive land use) and climatic changes, leading to vegetation degradation, soil erosion, surface water loss, bed-rock exposure.Increasing temperature and decreasing precipitation (Figure 9) has been witnessed, which has been proved to inevitably have profound effect on the reversal of rocky desertification [75].Local changes in temperature and precipitation can influence disturbance [91].The increasing temperature enhanced forest growth and the decreasing precipitation reduced the water erosion of rock desertification, slowing down the desertification speed [75], leading to an annual AGB increase.Also, extreme low temperature especially in 2008 affected vegetation growth [86].In the past 30 years, the area of rocky desertified land has begun to decrease (Figure 10), over-cultivation, over-cutting and over-grazing have gradually decreased or were banned completely under the guidance of national environmental and ecological policies.
Rock desertification was loaded strongly (0.989) among human activities factors.This has proven to be a complex process that has caused forest degradation and loss, soil erosion and water loss [75].In reality, it was mostly triggered by intensive human activities [84], such as slash-and-burn cultivation and destroying forest for land reclamation, which can be treated as a severe disturbance of the mountain region (Figure 2a).Thus, human activities played a major role while climate change played an indirect role in short term AGB dynamics, which is in agreement with the results presented in Table 7.

Uncertainty in Detection of AGB and Forest Disturbance
Overall, the findings of this study provide the necessary insight and motivation to policy makers and remote sensing community in forest protection, management and carbon storage, particularly for plantations regions.However, field-based AGB estimates themselves have many sources of sampling and measurement error, as well as model misspecifications [8].Future studies should pay more attention to more seasonal time-series data [5] and the integration of high resolution imagery or lidar data [92] to improve in accuracy of biomass estimation.Techniques need to be improved to include forest disturbance information in modeling to explain much more about AGB dynamics.Despite of the natural forests here, it is completely interesting to monitor the deforestation and afforestation activities of specific plantations species, which is provided with the characterization of fast-growing and high yield.

Conclusions
In this study, we developed a method to estimate the forest AGB distribution in northern Guangdong through the combination of multi-source datasets.Six years of NFI plot data in combination with spectral variables, topography data, texture measures and PALSAR HH, HV, HH/HV were collected for forest AGB mapping.The AGB map in 2006 shows highest R 2 of 0.8 and the lowest NRMSE of 7.49% with an independent test set.A strong increasing trend of mean AGB from the lowest value of 13.58 t/ha to the highest value of 66.25 t/ha occurred between 1988 and 2000, while after 2000 there was a fluctuating ascending change.The trend in disturbance area closely corresponds with the trend in AGB decrease.Human activities contributed more to short-term AGB dynamics than climate factors.Findings from this study will form the foundation for exploring the performance of forest AGB under human and climatic disturbance in subtropical regions, which can help policy makers and remote sensing community understand forest protection, management and carbon storage.

Figure 1 .
Figure 1.Landsat p122r043 footprint intersects with red polygon and color composite map of PALSAR images at the spatial resolution of 25 m for 2010 in a false-color combination of Red (HH), Green (HV), and Blue (HH/HV) in Northern Guangdong, showing the exact study site(black region).

Figure 1 .
Figure 1.Landsat p122r043 footprint intersects with red polygon and color composite map of PALSAR images at the spatial resolution of 25 m for 2010 in a false-color combination of Red (HH), Green (HV), and Blue (HH/HV) in Northern Guangdong, showing the exact study site(black region).

Figure 2 .
Figure 2. Example of Landsat time series (LTS) imagery (displayed in 5,4,3 false color composite; RGB uses band 5 images from 1984, 2008, and 2011) used to show rock desertification (a); harvesting (b) and spectral trajectories ((c), showing patterns of disturbance with arrows indicating year of onset) used to detect disturbance on random selected plots.A series of high resolution Google Earth photos (white color numbers) were also used to aid the interpretation of the Landsat data.

Figure 2 .
Figure 2. Example of Landsat time series (LTS) imagery (displayed in 5,4,3 false color composite; RGB uses band 5 images from 1984, 2008, and 2011) used to show rock desertification (a); harvesting (b) and spectral trajectories ((c), showing patterns of disturbance with arrows indicating year of onset) used to detect disturbance on random selected plots.A series of high resolution Google Earth photos (white color numbers) were also used to aid the interpretation of the Landsat data.

Figure 3 .
Figure 3. Variable importance for the whole predictor variables.Figure 3. Variable importance for the whole predictor variables.

Figure 4 .
Figure 4. Variable importance for stronger importance predicted variables.

Figure 4 .
Figure 4. Variable importance for stronger importance predicted variables.

Figure 6
Figure 6 illustrates the Google Earth aerial photo with RF uncertainty regarding the mean, standard deviation and CV in 2011 for regions of northern Guangdong.When the CV value was high, the mean of AGB was low.The map of RF uncertainty highlighted small anomalous areas of high uncertainty.These are localized areas where the predictions from individual trees varied widely, with some of the trees in the RF predicting low AGB and other trees predicting high AGB.

Figure 7 .
Figure 7. Scene-levelbiomass trajectories with NFI plot-based estimates shown for available years.Figure 7. Scene-levelbiomass trajectories with NFI plot-based estimates shown for available years.

Figure 7 .
Figure 7. Scene-levelbiomass trajectories with NFI plot-based estimates shown for available years.Figure 7. Scene-levelbiomass trajectories with NFI plot-based estimates shown for available years.

Figure 7 .
Figure 7. Scene-levelbiomass trajectories with NFI plot-based estimates shown for available years.

Figure 8 .
Figure 8. Annual aboveground biomass (AGB) decrease and VCT-based forest disturbance area.Although overall ascending trend of mean biomass was reported, forest disturbance deduced AGB decrease annually.

Figure 8 .
Figure 8. Annual aboveground biomass (AGB) decrease and VCT-based forest disturbance area.Although overall ascending trend of mean biomass was reported, forest disturbance deduced AGB decrease annually.

3. 4 .
Spatio-Temporal Multi-Scale Driving Factors of Forest AGB Dynamics 3.4.1.Regional Climate Change To investigate climate change over northern Guangdong Province, we analyzed trends in the average values of temperature and precipitation of two stations from 1980 to 2013 in this region.

Figure 9 .
Figure 9. Inter-annual variations in average values of mean temperature ( ˝C), precipitation (mm) based on Shaoguan and Heyuan meteorological stations.
Remote Sens. 2016, 8, 595 14 of 23 indirectly associated with AGB decrease, and the annual forest disturbance and mining changes had a great influence on the AGB dynamics.

Figure 10 .
Figure 10.Variations in biomass decrease and human activities.

Figure 10 .
Figure 10.Variations in biomass decrease and human activities.

Table 1 .Detail of
the field data showing plot level information for aboveground biomass (AGB).

Table 3 .
Summary of the four types of predictor variables of spectral, topography, texture and PALSAR components.
and 2011, and the observed data fromNFI in 1988NFI in  , 1992NFI in  , 1997NFI in  , 2002NFI in  , 2007NFI in   and 2012.   .

Table 4 .
Descriptive statistics of the validation dataset.

Table 7 .
The rotated component matrix of climate force, human activities and biomass decrease in northern Guangdong.