Mapping Annual Forest Change Due to Afforestation in Guangdong Province of China Using Active and Passive Remote Sensing Data

Accurate acquisition of spatial distribution of afforestation in a large area is of great significance to contributing to the sustainable utilization of forest resources and the evaluation of the carbon accounting. Annual forest maps (1986–2016) of Guangdong, China were generated using time series Landsat images and PALSAR data. Initially, four PALSAR-based classifiers were used to classify land cover types. Then, the optimal mapping algorithm was determined. Next, an accurate identification of forest and non-forest was carried out by combining Landsat-based phenological variables and PALSAR-based land cover classifications. Finally, the spatio-temporal distribution of forest cover change due to afforestation was created and its forest biomass dynamics changes were detected. The results indicated that the overall accuracy of forest classification of the improved model based on the PALSAR-based stochastic gradient boosting (SGB) classification and the maximum value of normalized difference vegetation index (NDVI; SGB-NDVI) were approximately 75–85% in 2005, 2010, and 2016. Compared with the Japan Aerospace Exploration Agency (JAXA) PALSAR-forest/non-forest, the SGB-NDVI-based forest product showed great improvement, while the SGB-NDVI product was the same or slightly inferior to the Global Land Cover (GLC) and vegetation tracker change (VCT)-based land cover types, respectively. Although this combination of multiple sources contained some errors, the SGB-NDVI model effectively identified the distribution of forest cover changes by afforestation events. By integrating aboveground biomass dynamics (AGB) change with forest cover, the trend in afforestation area closely corresponded with the trend in forest AGB. This technique can provide an essential data baseline for carbon assessment in the planted forests of southern China.


Introduction
Historical land use, land cover changes, and forest management practices affect the exchange of carbon between sources and sinks in forests [1].Afforestation and reforestation programs in temperate regions have increased forest areas [2] and carbon accumulation since the 1970s, and this accumulation has significantly contributed to global terrestrial carbon sinks [3,4].China has the largest afforested area in the world (~62 million ha in 2008), and most of them are carbon sinks [5], while southern China accounts for 65% of the forest carbon sink in China, especially for the fast-growing tree species [6].However, historical time series of forest cover maps due to afforestation are still not available for generating the spatio-temporal dynamics of afforestation carbon storage or its biophysical mechanisms in response to climate change [7,8].Traditional methods that have been used to monitor forest change (e.g., afforestation) have relied on permanent sample plot (PSP) surveys at the provincial scale, and these have been used for the national forest resources inventory in China [9].However, the effectiveness of PSP for monitoring forest cover has been limited due to inadequate spatial coverage and a five-year survey rotation.
Time series remotely sensed data have been considered an effective spatial detection tool to monitor long-term forest cover changes at large scales [10][11][12][13].Past studies have identified forest cover primarily based on coarse resolution datasets (e.g., 300 m to 1 km) [14][15][16][17][18], thus leading to limited utility, especially at local scales [14].Consequently, 30 m resolution Landsat-like products (e.g., China's POK (pixel and object)-based 30 m GlobeLand30 (GLC30), Global Forest Change (GFC) data), or 25 m resolution PALSAR-based products (e.g., ALOS PALSAR-based global forest/non-forest mosaics) have been generated [11,13,[19][20][21], and a variety of time series forest disturbance detection techniques and products have been developed by using automated and semi-automated algorithms [22][23][24][25].However, bio-temporal or triple-temporal satellite images cannot capture the forest change spectrally [14,19,26].It is much more difficult and complex for time series algorithms to develop the above-mentioned medium to high resolution land cover products at multiple temporal and spatial scales due to many factors; for example, time series optical remote-sensing-based products need good quality observations without the limits of cloud or cloud shadow, and to save time and labor in big data processing, and to minimize the spectral confusion, etc. [13,23,24].Integration of optical sensors (e.g., Landsat and MODIS), radar sensor (Synthetic Aperture Radar (SAR), e.g., Phased Array L-band Synthetic Aperture Radar (PALSAR) on the Advanced Land Observing Satellite (ALOS)), and light detection and ranging (lidar) or high-resolution sensors can more accurately capture three-dimensional structures needed to delineate forest cover [27][28][29][30][31][32][33]; however, the latter two have limited spatial and temporal coverage due to the high costs [28,34].
Cloud-free L-band SAR has been shown to be advantageous for monitoring cloudy and rainy tropical or sub-tropical areas [35], and the optical remote sensing, e.g., Landsat, also has the potential to balance the deficiencies of radar data, e.g., PALSAR, in distinguishing between forest and other confusing types (e.g., rock, building, and urban) [36].Several successful published studies have investigated forest cover based on optical and radar data at different spatial scales [35,[37][38][39][40].However, previous studies have commonly been carried out in a single year [35,36], or multiple data comparisons (e.g., RapidEye, TM, PALSAR, Envisat ASAR)-based land cover (100 km by 100 km test site) monitoring by visual interpretation [41], or have taken multiple years to map forest cover  and to quantify forest encroachment into grasslands [42,43].So, the integration of multi-sensor and multi-temporal remote sensing systems, including Landsat-like optical sensors and SAR, shows great potential to develop dense time-series forest mapping projects and assist with dynamic monitoring endeavors.
Supervised classifications (e.g., support vector machine (SVM), boosting tree (stochastic gradient boosting (SGB)), decision tree, and random forest (RF)) [44][45][46] are more effective than unsupervised classification and object-oriented classification in terms of time series and large scales.The SVM classifier has been widely reported as an outstanding classifier in remote sensing [45].The RF classifier has been tested due to its reported performance in the machine learning community [44].The SGB algorithm usually outperformed traditional parameter or non-parameter methods (e.g., classification and regression, CART; RF) [47][48][49] in land use, land cover classification [48], and forest fuel type mapping [47].Furthermore, the SGB algorithm has been used to generate land cover types based on multispectral and hyperspectral images of individual years (e.g., IKONOS, Landsat ETM+, and Probe-1) [48,50].However, previous studies have proven that PALSAR-based machine learning algorithms in forest cover mapping had some commission or omission error when used alone [36,51,52].Temporally frequent Landsat data have long-term archives and free availability and similar image data [22], and the generated time-series spectral vegetation index can potentially provide vegetation phenology patterns, which are particularly useful in environments with limited accessibility and a lack of in situ measurements [53].The use of such an index can help us understand vegetation dynamics with regard to climate change impacts on vegetation identification, such as vegetation-greenness-related normalized difference vegetation index (NDVI) [54].Furthermore, in regard to the remote sensing change detection, forest cover produced by the plantation afforestation is defined as "forest stands that have been established artificially, either on land that has not supported forests in the last 50 years (i.e., afforestation), or on land that has supported forests in the past, but where the original vegetation has been replaced by forests (i.e., reforestation)."[55].Therefore, dense time series and high resolution free and open access data, e.g., Landsat or PALSAR, hold the ability to form dense time observations to generate the robust forest cover change due to afforestation.
The objective of this study is to extract annual forest change (1986-2016) due to afforestation in Guangdong, China.First, a novel procedure to identify and map annual forest cover caused by afforestation based on the integration of the PALSAR-based spectral and textural values and Landsat-based phenological variables is developed and tested.Then, the accuracy of the developed procedure is validated and compared using other forest/non-forest (FNF) products.Finally, forest aboveground biomass dynamics (AGB) under afforestation changes are investigated.

Study Area
The study area is the Guangdong Province (20 17.97 × 10 4 km 2 , Figure 1) in China.The local topography is undulating (elevation 22-1353 m above sea level).The climate varies from subtropical to tropical.The annual mean precipitation is 1300-2500 mm, and the average temperature ranges from 19 to 24 • C. The wet season occurs from April to September, and the dry season is from November to January (February, March, and October are transitional months).In March and April, the northern region is often wet, but the southern region is dry [56].In September, the pattern is reversed [56].Most of the forest species are considered evergreen and fast-growing [57].The most common extreme meteorological disaster includes chilling injury, storms and flooding, and drought [57].

Active-and Passive-Based Satellite Data
This study used radiometrically and geometrically corrected PALSAR mosaic data from Japan Aerospace Exploration Agency (JAXA) according to topography and atmospherically corrected Landsat data from USGS/EROS processes.There were six years of PALSAR mosaics used through raw-strips-based tiles from July to September (Table 1, https://www.eorc.jaxa.jp/ALOS/en/palsar_fnf/data/index.htm):tiles covering p120r043, p121043, p122r043, and p123r043 in September were dry season data, while the remaining paths/rows (p120r044, p121r044, p122r044, p122r045, p123r044, p123r045, p124r045, and p124r046) in September and all of tiles from July to August were wet season data.The proportion of Landsat images (e.g., by paths/rows (Figure S1a), months (Figure S1b), sensors (Figure S1c), and phenology (dry or wet season, Figure S1d)) was counted in Figure S1.
First, we converted the original PALSAR data to the backscatter coefficient in decibels, then implemented the enhanced Lee filter (window size, 5 × 5 pixels) to reduce speckles.Next, we produced some variables for land cover classification including HH, HV polarizations, HH/HV (ratio), HH-HV (difference), and HV texture measures (window size, 3 × 3 pixels, offset ([1,1]), and a 64 gray level quantization), which have been proven to distinguish well between forest and non-forest [21].Finally, the 25 m PALSAR mosaic data and their derivatives were re-projected using Landsat data to match the 30 m spatial resolution (Table 1).
The highest quality cloud-free images from the 12 Landsat path/row tiles (Figure 1, Figure S1) were used based on the Google Earth Engine cloud computing platform (https://earthengine.google.com).ETM+ data after the year 2003 were filled based on the USGS LS7 SLC-off gap-filling algorithm, which was recreated for the Google Earth Engine by Noel Gorelick (https://code.earthengine.google.com/20cba5268cbe117e2fc1c5fefc33f3) (Figure 2).

Extraction of PALSAR Backscatter Signatures for Land Cover Types
Ground truth samples in the regions of interest (ROIs) for forest (5841 polygons), cropland (5544 polygons), water (2267 polygons), urban (7036 polygons), and other types (short for others, 6474 polygons), were selected using 12 paths/rows of the Google Earth high resolution images (Figure 3) referring to National Forest Inventory (NFI) and sub-compartment data ("xiaoban" (XB) in the Forest Management Planning Inventory (FMPI).All of the ROIs were extracted in locations where only a single land cover type covered the area.They can be easily downloaded into different formats, such as Keyhole Markup Language files (.kml).Next, ArcGIS vector files (shapefile) were produced from ROIs in KML format.A series of land cover types were used (Figure 3), of which, a random 50% (13,581 polygons) of the total samples from six years was used as training data and overlaid on the PALSAR-based bands to classify the five land cover types by calculating the mean pixel value per polygon of the HH, HV, ratio, difference, and HV texture measures, and the chosen of the remaining samples were reserved as validation data for the classification accuracy assessment (Figure 2).Gaussian kernel density estimations [59] of the training ROIs (13,379) of land cover types (forest (2944 polygons, 21,800,050 pixels), other types (3153 polygons, 97,221 pixels, short for others), water (1114 polygons, 156,383 pixels), cropland (2739 polygons, 69,395 pixels), and urban areas (3429 polygons, 84,293 pixels)) in the bands of HH, HV, HV textural measures, Ratio, and Difference from all six years indicated the separability among these land cover types, especially water had lower HH and HV values, which can be easily identified (Figure 4).Urban has high HH values and lower ratio values, and can be identified based on these two indexes (Figure 4a,c); forest and urban have high HV values (Figure 4b), and lower values of forest can be seen in Difference (Figure 4d).Cropland can be identified in HV, while other types overlapped over the cropland (Figure 4b).Here, no obvious cost for HV-based texture measurements was found in distinguishing forest from the other non-forest types (Figure 4e-l).

Classification Algorithms
The support vector machine (SVM) classifier, RF classifier, stochastic gradient boosting (SGB) classifier, and C5.0 decision tree classifier were implemented to classify land cover types from the multi-temporal mosaic PALSAR and its derivatives (HH, HV, Ratio, and Difference, etc.) according to the above evaluation (Figure 2).The parameter settings used for each of the classifiers are listed in Table 2.
A parallel computing technique was performed for the SVM, RF, and C5.0 classifiers to improve the efficiency of large area image processing in R (R Development Core Team, 2008) [60].Gradient Boosting Machines (GBM) is an inherently sequential algorithm.The GBM package for SGB relies on a built-in parallel package [61].Each iteration depends on the results of the previous iteration.GBM creates an ensemble of decision trees that build on top of each other.Each tree predicts the error of the next tree.When combined this ensemble performs extremely well [61].The search for the best decision tree is done across the n.cores.GBM uses cross-validation to identify the best number of decision trees (either using the training or testing error).Contrasting this with random forest, where the algorithm is exceedingly parallel, every tree is independent of each other.This allows RF to be parallelized in ways GBM cannot.
"doParallel" and "foreach" [62] in R interface were used to perform collaborative parallel processing."foreach" allows for the creation of small trees, and they are then combined using the "combine" function.All of the images used for land cover classification were segmented based on the number of computer CPU cores along the latitude using "doParallel".Any of the available CPU cores were called to process the segmentation part using "foreach".Finally, all of the results from the segmented parts after processing were mosaicked.
A majority filter by calculating focal values for the neighborhood of the default moving window (3 × 3 pixels) based on "raster" package was applied to reduce the issue of "salt and pepper".Then a mathematical morphology opening operation (5 × 5 pixels) by "mmand" package was selected to eliminate the speckled and smooth boundaries, and to exclude the burrs and isolated pixels [63].  2 and 3).Also, confusion matrix plots that compared the actual and predicted classes for five items were produced.The tiles were colored according to the frequency of the intersection of the two classes, thus, the diagonal represented where the actual class was predicted correctly.The color represented the relative frequency of that observation in the data; given some classes occur more frequently, the values were normalized before plotting.Any row of tiles (save for diagonal) represented instances where items were falsely identified as belonging to the specified class.Finally, PALSAR-based land cover products (forest/non-forest, FNF) from the optimal classification algorithm were used for further study.

Further Forest Mapping Based on the Integration of PALSAR-Based FNF and Landsat Data
In this study, a similar pixel-based approach was used as that by Shen et al. (2018) to generate annual forest maps .The difference was that the newly produced PALSAR-based FNF data in Section 2.4 were used to integrate with Landsat-based phenological variables to map the forest, because PALSAR-based forest is often confused with other categories, for example, the commission error between forest and crops and grasses because of their different phenological patterns [58], or having some omission error associated with urban (buildings) and other features (Figure 4) in the performance of the similar PALSAR backscatter coefficients.The NDVI max Gaussian kernel densities for training ROIs of the forest, other types, water, cropland, and urban were plotted (Figure 5).
Figure 5 shows that forest has higher NDVI max values, follow by cropland, others, and urban, therefore, the highest Landsat-based NDVI vegetation greenness variables (NDVI max ) can be derived to differentiate similar high backscatter coefficients and different phenological patterns between forest and other non-forest (Figure 5) based on the previous studies [36,64].That is, statistical analysis based on the training ROIs among all of the types showed that the threshold value (greater than or equal to 0.72, a black dashed in Figure 5) of NDVI max was used to separate 80% of the forest pixels (>0.72) from 80% of the others' pixels (<0.72), 85% of the forest pixels (>0.72) from 85% of the cropland pixels (<0.72), 95% of the forest pixels (>0.72) from 99% of the water pixels (<0.72), and 95% of the forest pixels (>0.72) from 99% of the urban pixels (<0.72) (Figure 5).To reduce the error of the PALSAR-based forest (Figure 4) and further solve the mixed classified forest, a rule was built to eliminate those misclassified pixels in forest (commission error) and identify forest from other classified pixels (omission error) as follows: assume the PALSAR-based forest from Section 2.4 as 1, the PALSAR-based non-forest as 0; then, if PALSAR-based forest = 1 and ndvi max > 0.72, then a pixel is classified as forest to avoid the commission error; or if PALSAR-based non-forest = 0 and ndvi max > 0.72, then it is also classified as forest to avoid the omission error.Next, a median filter (window size 5 × 5) was adopted to solve speckle noise issues (e.g., salt-and-pepper noise) in the annual forest maps.Due to no long time-series PALSAR data, we used the PALSAR data in 2007 and 2015 to detect forest cover before 2007, and from 2011-2015, respectively (Figure 2).It was then deduced that the 2007 and 2015 PALSAR data produced the most accurate maximum forest area (8.27 × 10 6 ha, 9.06 × 10 6 ha, respectively) from the Chinese national forestry yearbook (2011)(2012)(2013)(2014)(2015), respectively) (Figure 2) [58].The land cover maps from 1986-2006 (2011-2015) were substituted by the PALSAR-based land cover map in 2007 (2015), then the above rule was also applicable to other years of forest/non-forest mapping (Figure 2).

Evaluation of PALSAR/Landsat-Based Forest Maps
A confusion matrix based on the validation plots was selected to assess the eventual forest maps.Half of the total plots for forests and non-forests were acquired based on data from Google Earth, NFI, and sub-compartment data from 2005 (520 polygons, 1641 polygons), 2010 (487 polygons, 1833 polygons), and 2016 (505 polygons, 1888 polygons) (Figure 2).

Evaluation of the PALSAR/Landsat-Based Forest Map with Mutlitple Forest Cover Products
Forest ROIs and non-forest ROIs in 2010 (4,871,833) were collected for validation of different forest cover products with PALSAR/Landsat-based forest map in 2010.We validated and compared the overall accuracy, kappa coefficient, user accuracy, and producer accuracy of forest classification and the total forest area among all of the forest cover products, including China's 30 m GlobeLand30 (GLC30) (Figure 1), JAXA PALSAR-FNF, vegetation change tracker (VCT) product (Table 3).PALSAR-FNF data was resampled from 25 m resolution to 30 m spatial resolution with nearest neighbor interpolation to make it consistent with other 30 m products.Annual forest maps derived from the PALSAR/Landsat-based FNF were used to provide forest change data to map the spatial pattern of afforestation distribution from 1986-2016 (Figure 2).Persisting forest (non-forest) indicated that the cover type of a pixel remained the forest (non-forest) during the entire observation period of the annual forest maps.Annual afforestation (e.g., 2016) was generated as the intersection between persisting non-forest from the year before the current year (e.g., 2015) to the starting year (e.g., 1986) of the entire annual forest maps and forest in the current year (e.g., 2016) (Figure 6).Per five (or four, or six) years afforestation (e.g., 2011-2016) was generated as the intersection between persisting non-forest from the start year (e.g., 2010) of the per time interval (e.g., 2010-2016) to the start year (e.g., 1986) of the entire annual forest maps and forest in the end year (e.g., 2016) of the per time interval (Figure 6).The forest AGB time-series stack was produced using ICESat/GLAS measurements, historical inventory data, and time-series optical and radar imagery.Further information about the algorithm is available in Shen et al. (2018).The combined remotely sensed algorithm for mapping AGB yielded a good accuracy (R 2 adj = 0.86, n = 558, p < 0.001, RMSE = 11.35t/ha).We estimated the forest cover area by afforestation and quantified AGB change depending on the "RF-based GLAS waveform-extrapolated footprint AGB model".The AGB time series stack was intersected with the above annual map of forest cover (Figure 6).The final annual AGB map with forest cover change map was clipped by the study area boundary.The trend of forest cover area changed due to afforestation associated with forest AGB ("afforestation-AGB") was counted based on individual year (e.g., 1990, 1995, 2000,..., 2010, 2016 ) or time intervals (e.g., 1987-1990, 1991-1995,..., 2006-2010, 2011-2016).

Analysis of Land Cover Types Classification from PALSAR
The parallel processing technique used on a regular single desktop computer increased the image computation efficiency per scene (referring to the Landsat footprint) approximately 9-10 times and had a shorter waiting time of approximately 20 minutes than the traditional per scene processing (more than 3 hours) (Code S1), and this was conducive to enhancing the efficiency of the classifiers for the PALSAR-based classification.
Table S1 shows the accuracy assessment of the PALSAR-based land cover classifications with ground-truth-based regions of interest (ROI) using four classifiers (SVM, RF, SGB, and C5.0).A total of 555 (546) ground truth forest polygon ROIs, 201 (186) water polygon ROIs, 467 (635) other polygon ROIs, 508 (588) cropland polygon ROIs, and 558 (631) urban polygon ROIs in 2007 (2016) were used for validation.The overall accuracy range of the four algorithms was 58.02-63.61%and 61.49-62.34% in 2007 and 2016, respectively.The Kappa coefficient range was 0.466-0.535and 0.502-0.513 in 2007 and 2016, respectively.In regard to the over accuracy, the SGB (RF) classifier was superior to the other classifiers, followed by the SVM (SGB) in 2007 (2016), respectively.Due to the consideration of user accuracy, the SGB classification results were the most effective among the four algorithms, especially for forest classification.However, in terms of producer accuracy, SVM classification results for forest classification were superior to other algorithms.Overall, any index may not separately determine which algorithm (SGB, SVM) was the best, while the SGB-based results showed a good balance of these indexes.
In the rendered plot (Figure 7), it can be observed that water and forest were identified as items belonging to all of the other classes in 2007 and 2016.The annual land cover maps with 30 m resolution from PALSAR were created using four classifiers.Figure 8 and Figure S2 show the land cover classification map in 2007 and 2016, respectively.The optimal SGB classification results were used to further distinguish forests and non-forests.

Assessment of PALSAR/Landsat-Based Forest/Non-Forest Mapping in Guangdong
The validation results of the PALSAR/Landsat-based forest maps with ROIs were demonstrated in Table 4 and Table S2.Ground truth forest polygon ROIs and non-forest ROIs were collected in 2005 (555, 1734), 2010 (518,1968), and 2016 (546, 2043) for validation.The overall accuracy was over 75% (95% CI: 75.11-78.6% in 2005), and up to approximately 85% (95% CI: 83.28-86.2% in 2010).The Kappa coefficient was over 0.45, and up to approximately 0.58.In regard to the producer accuracy, the PALSAR/Landsat-based forest mapping (85.5%) was superior to the PALSAR-based SGB forest mapping (66.48%,Table S1) in 2016.Apart from the user accuracy, the final results presented that the integration of the PALSAR-based SGB classification and the maximum value of NDVI ("SGB-NDVI")-based forest map had much better accuracy than that of the single PALSAR classification of FNF in Section 3.1.Originally, the overall accuracy of global land cover (GLC30) product for the year 2010 from Landsat TM/ETM+ and HJ-1 was 83.5% ± 0.18%, and the user accuracy of the forest classification was 89.00% [19].
Here, the differences between the overall accuracy of the forest classification from GLC30 (85.8%) and our results (SGB-NDVI-based FNF map) (84.8%) are less compared to the discrepancies between the JAXA PALSAR-FNF (80.7%) and our results for the entire Guangdong scale (Table 5 and Table S3).An assessment of forest and non-forest from the integrated forest z-score (IFZ)-based vegetation change tracker (VCT) product and our results in northern Guangdong (p122r043) was performed.This was used to show that the VCT-based forest product (90.3%)was superior to our results (86.1%) with regards to the overall accuracy and Kappa coefficient at a small scale (Table 5 and Table S3).The total forest area from the SGB-NDVI-based forest map in 2010 was calculated to be 8.53 × 10 6 ha in Guangdong, which was close to the results of the national forestry yearbook of China (8.74 × 10 6 ha), but lower than the calculation from the GLC30 map (9.59 × 10 6 ha) and higher than the calculation from the JAXA PALSAR-FNF map (7.83 × 10 6 ha).The areas of the PALSAR-based forest maps in Section 3.1 were about 8.33 × 10 6 ha, 8.13 × 10 6 ha, 8.22 × 10 6 ha, and 8.30 × 10 6 ha by SGB, SVM, RF, and C5.0 classifiers in 2010, respectively, which were lower than that of the SGB-NDVI-based forest map and the national forestry yearbook.

Relationships Between Forest Cover Change Dynamics (Afforestation) and Forest AGB
The annual forest cover change maps under afforestation (Figure 9) were created using the afforestation spatial pattern calculation based on the method in Figure 6.The trend of forest cover area changed due to afforestation associated with forest AGB was shown in Figure 10.Overall, the trend in afforestation area closely corresponded with the trend in forest AGB, except for a couple of notable anomalies.In Guangdong, the highest afforestation area and AGB value was observed during the period of 1991-1995, while the lowest was during the period from 2006-2010 (Figure 10a).In a single year, the afforestation area and forest AGB value tended to increase to the maximum, followed by a decrease up until 2010, then they continued to rise (Figure 10b).Northern Guangdong (p122r043/p121r043) (Figure 10c, d) was tested due to the highest forest AGB of the entirety of Guangdong Province [58].
The dramatic decline in the afforestation area was significant in 2000 (Figure 10d).Unlike Heyuan, Qingyuan City and Shaoguan City (p122r043) had no obvious fluctuation.

Choice of Mapping Algorithms
A variety of algorithms sensitive to land cover type classification were investigated.HH, HV, ratio, and difference contributed more in the separability evaluation of PALSAR-based classification, while a conclusion about the limited potential of the HV-based texture measures has been drawn (Figure 5), which was opposite to the well performance description in [21], because the training areas included pixels near from the edge, and texture measures are strongly influenced near edges due to the mixed pixels, especially in large analysis windows or multiscale analysis [65].
The overall accuracy (Kappa coefficients) of PALSAR-based classifications using SVM, RF, SGB, and C5.0 were not higher than 65% (0.54), showing the defects of the PALSAR-based land cover type mapping results directly because of the low accuracy, which can be explained as due to the PALSAR data lacking some regions because of the strong ionospheric distortion, especially near the image edge [66].Otherwise, there were color differences between two adjacent scenes in the mosaic images, which may have been caused by changes in the backscatter intensity induced by the freezing of trees in winter [67], which ultimately affected the PALSAR-based classification of forest and non-forest.However, the SGB classifier contributed more steadily, although there was a discrepancy from the previous study [68].Also, the area of the SGB-based forest classification in 2010 (8.33 × 10 6 ha) was closest to the true survey data (8.74 × 10 6 ha) from the national forestry yearbook of China among all of the classifiers.Usually, although specific parameter tuning of four classifiers need to be adjusted, the basic and default parameters used in classification can satisfy needs; for example, the SGB model is quite sophisticated since changing any setting can affect the optimal values of other settings [69].Moreover, high biomass crops were likely to be misclassified as forest, as the seamless PALSAR mosaic product was generated using the SAR image every summer from July to September, and the image data acquisition dates are equivalent to crop maturity dates with high biomass [70].Some of the raw strips comprising each tile were acquired during the wet season, the dielectric constant of moisture or water can affect radar backscatter [71] and may influence the results of the analysis.
After the integration of the Landsat-based NDVI max variable with PALSAR, the forest/non-forest classification led to accuracies (Kappa coefficients) ascension of up to 85% (0.6) in the current work.The area of the FNF in 2010 has gone up than that of the PALSAR-based FNF classification.Because the commission error and omission error between forest and non-forest has been improved, especially the urban (buildings) areas were included from the PALSAR/Landsat forest maps, while the croplands (or other types) were excluded led to the discrepancy of 0.2 × 10 6 ha between the true survey data and PALSAR/Landsat forest maps (SGB-NDVI-based forest map) (Figures 5 and 6).It was a converse result with Qing et al. 2016, where after including Landsat NDVI max , about 10% (~4000 km 2 ) of shrub, buildings, and rocky land were reduced in the area of the PALSAR/Landsat forest maps, and no obvious increase had been found.Moreover, Landsat observations during the vegetation growing season are limited.In the northern areas of Guangdong, the accumulated NDVI value during the dry season does not fully represent local phenological information because tree species are not completely evergreen, resulting in the underestimation of forest cover.So, when possible, wet season Landsat images were also included to produce the accumulated NDVI value for the full year.Generally, the number of dry season data that were used was more than the number of the wet season data.With sufficient quantity and superior quality without regard to the wet season, the maximum NDVI of dry season Landsat data can achieve good performance [58].The accuracy of forest cover maps during the years without PALSAR data also had acceptable accuracy (Table 4 and Table S2).Actually, many other time-series of vegetation indices (EVI, GNDVI, and NDWI) have been used to estimate vegetation phenology; however, most of them were derived from coarse resolution imagery from MODIS [53,72,73] or AVHRR observations [74].

Comparisons of Forest Cover Maps and the Existing Results
In view of the different definitions, data sources, and mapping methods that have been used to map forest cover [75], there are striking differences among land cover maps [36,76].The forest area from GLC30 was obviously greater in extent than that found in the national forestry yearbook of China and in our results, which may have resulted from the fact that the sparse woodland was also regarded as forest [19].The Landsat-based VCT algorithm is used to detect dense time-series forest changes and the VCT-based forest/non-forest product achieved outstanding performance; however, it must rely on images from the peak vegetation growing season [23], but these imageries cannot be guaranteed to be accurate in coastal or low-latitude regions.In a comparison of these results with forest maps created by an integration of the JAXA PALSAR global FNF map [21] and Landsat by Shen et al. (2018), we found the forest area of the latter was lower.The validation accuracy of the results generated using the "SGB-NDVI" algorithm was superior to that of the JAXA PALSAR FNF product [21].This is because the JAXA PALSAR-based forest from the FNF map is defined as areas with canopy cover of natural forests over 10%, and the area must be larger than 0.5 ha (http://www.eorc.jaxa.jp/ALOS/en/palsar_fnf/DatasetDescription_PALSAR2_Mosaic_FNF_revE.pdf).However, in addition to natural forests, there are a wide range of planted forests in southern China.The evaluation indexes (overall accuracy, Kappa coefficient, user accuracy, producer accuracy, and area) of the findings from this study were acceptable when time-series forest classification maps were produced that were based on an activeand passive-based improved algorithm.

Forest Cover Dynamics Change Due to Afforestation and Forest AGB
We examined the forest cover change dynamics.The forest area consisting of planted forests from 2006-2010 increased by 0.63 × 10 6 ha in Guangdong based on the national forestry year book of China, while that of afforestation was 0.59 × 10 6 ha.This difference is because the definition of planted forests contains new afforestation [55].The local government proposed a fast-growing eucalyptus plan in 1995 and the eucalyptus plan and slope improvement plans were discontinued in 2000.Furthermore, urban expansion possibly affected AGB changes under afforestation across Guangdong Province, which resulted in afforestation area combined with forest AGB increase and undulation, as Heyuan City in northern Guangdong (p121043) was the first to ban eucalyptus tree planting, and explains the sharp decline of afforestation area in 2000 (Figure 10d).Additionally, because the area covered by p121r043 is adjacent to the southwest area of Jiangxi Province, afforestation projects have increased in the past 20 years, so there is a large amount of afforestation area.

Uncertainties in the Detection of Forest Change Due to Afforestation
The uncertainties in the mapped historical forest distribution can be explained by poor data quality, inadequate data acquisition date, and errors in algorithm implementations [13,[77][78][79].To obtain greater spatial and temporal resolution observation capability, a combination of free and open access multi-source data (e.g., Landsat, PALSAR, and Sentinel)-based algorithms and high-performance computing systems for big data analysis [38,79], such as the NASA Earth Exchange (NEX) [80], are needed to provide better automatic extraction of seamless time-series forest change products.

Conclusions
In this study, a mapping method to detect changes in forest distribution under afforestation in Guangdong Province of China was developed using a combination of the PALSAR-based mosaic products and the dense time series Landsat-based phenology variable obtained from the Google Earth cloud platform.The final model was used to effectively construct an active-and passive-based forest cover detection framework.As validated by field measurements, the detection model generated reliable forest cover maps with some basic classification errors.By integrating spectral variables and phenology variables, the interannual and seasonal and spatio-temporal dynamics of changes in forest cover due to afforestation change were obtained.The combination of multiple sources and algorithms (advanced computing techniques, the optimal machine learning algorithms and remote sensing information) to develop models is a useful methodology, although inevitably there were some errors generated.This study shows that an integration of active and passive remote sensing data-based big data processing can fill in the lack of image data in low-latitude coastal areas and detect historical forest cover changes caused by afforestation.The findings from our study can improve the automatic identification of forest cover types.Future study could comprehensively incorporate multitemporal satellite observations and in situ measurements including lidar data, hyperspectral data, unmanned aerial vehicle (UAV), and forest structures (e.g., forest species, leaf area index, and forest age) to produce more accurate forest cover maps.This technique can provide a basis for understanding carbon dynamic related to forest biomass due to planted forests afforestation.

Figure 1 .
Figure 1.Twelve Landsat Paths/Rows covering the Guangdong Province of China, showing the exact study area (The background map is from China's 30 m GlobeLand30 (GLC30) data product in 2010 [19]).

Figure 2 .
Figure 2. The detailed flowchart for mapping annual PALSAR/Landsat-based forest/non-forest (1986-2016).First, support vector machine (SVM), random forest (RF), gradient boosting machines (GBM), and C5.0 based on PALSAR mosaic data was applied together with training and validation ROIs from Google Earth to generate five types of land cover maps.Second, the integration of PALSAR and Landsat-based maximum normalized difference vegetation index (NDVI) was used to generate PALSAR/Landsat-based forest/non-forest maps and improved the mapping accuracy.Finally, annual forest cover change due to afforestation was developed and to explore the relationship with forest aboveground biomass dynamics (AGB) distribution.

2. 4 .
Different Classification Algorithms for Mapping Forest and Non-Forest Based on Multi-Temporal PALSAR 2.4.1.Evaluation of the PALSAR Backscatter Signatures for Land Cover Types

Figure 5 .
Figure 5. Kernel density distribution plots of forest/non-forest (cropland, water, forest, and other types) from the dry and wet season maximum time series Landsat-based NDVI values over the corresponding six years, and the black dashed indicted the threshold value (0.72) of NDVI max to distinguish between forest and non-forest.

Figure 6 .
Figure 6.The method for the identification of annual or per five (or four, or six) years forest cover change due to afforestation.

Figure 10 .
Figure 10.Analysis of the relationship between afforestation and forest AGB change in Guangdong Province (a, b); including northern Guangdong: p122r043, p121r043 (c, d).

Supplementary Materials:
The following are available online at http://www.mdpi.com/2072-4292/11/5/490/s, Figure S1: Statistics of the number of Landsat images used by (a) 12paths/rows, (b) 12 months, (c) 3 sensors, and (d) dry season and wet season from Shen et al. (2018), Code S1: Taking a case of the PALSAR-based SVM land cove type classification by traditional classification and parallel processing classification techniques, Table S1: The accuracy assessment of land cover classification in 2007 and 2016, Figure S2: The 2007 land cover classifications based on SVM (a), RF (b), SGB (c), and C5.0(d) in the Guangdong province of China,

Funding:
This work was jointly supported by the National Natural Science Foundation of China [31670552], and the PAPD (Priority Academic Program Development) of Jiangsu provincial universities.Additionally, this work was performed while the corresponding author acted as an awardee of the 2017 Qinglan project sponsored by Jiangsu Province.

Table 1 .
Summary of satellite data used in this study.

Table 2 .
Classification algorithms and parameter settings.

Table 3 .
Summary of the forest cover products used.
max 2.7.Integration of Forest Aboveground Biomass Change with Annual Forest Cover Changes (Afforestation)

Table 4 .
The accuracy assessment of forest and non-forest in 2005, 2010, and 2016.

Table 5 .
Validation results of different forest and non-forest products in 2010.

Table S2 :
The accuracy assessment of forest and non-forest in 2005, 2010, and 2016, TableS3: Validation results of different forest and non-forest products in 2010.Author Contributions: W.S. designed the study, analyzed the data, and wrote the paper.M.L. and C.H. helped in project design, paper writing, and analysis.X.T. helped in paper review and editing.S.L. helped in the original data preparation.A.W. helped in field work and data analysis.