Spatiotemporal Estimation of Bamboo Forest Aboveground Carbon Storage Based on Landsat Data in Zhejiang, China

Abstract

China is one of the countries with the most abundant bamboo forest resources in the world, and Zhejiang province is among the top-3 Chinese provinces with richest bamboo forests. For rational bamboo forests management, it is of great significance to study the spatiotemporal dynamic changes of Aboveground Carbon (AGC) stocks of bamboo forest in Zhejiang. In this study, remote sensing variables, such as spectral, vegetation indices and texture features of bamboo forest in Zhejiang, were extracted from 32 Landsat TM and OLI images got from four different years (2000, 2004, 2008 and 2014). These variables were subsequently selected with stepwise regression method to build an estimation model of AGC of the bamboo forests. The results showed that (1) the accuracy of bamboo forest remote sensing information extracted from the four different years was high with a classification accuracy of >76.26% and an accuracy of users of >91.62%. The classification area of bamboo forest was highly consistent with the area from forest resource inventory, and the area accuracy was over 96.50%; (2) the estimation model performed well in predicting the AGC in Zhejiang for different years. The correlation coefficient for estimated and measured AGC was between 63% and 72% with low root mean square error; (3) the derived AGC of the bamboo forests in Zhejiang province increased gradually from 2000 to 2014, with the AGC density of 6.75 Mg·ha⁻¹, 10.95 Mg·ha⁻¹, 15.25 Mg·ha⁻¹ and 19.07 Mg·ha⁻¹ respectively, and the average annual growth of 0.88 Mg·ha⁻¹. The spatiotemporal evolution of bamboo forest AGC in Zhejiang province had a close relationship with the gradual expansion of bamboo forest in the province and the differentiation of management levels in different regions.

1. Introduction

Bamboos are naturally distributed in the tropical, subtropical and temperate regions of all the continents except Europe and western Asia, from lowland up to 4000 m in altitude [1], and mainly in Asia [2,3,4,5,6,7]. China locates in the center of world’s bamboo distribution area, and is the most important bamboo industry country in the world. Bamboo area in China reached 6.01 million hectares, accounting for approximately 20% of the total area of bamboo forests in the world [8]. Zhejiang, locating in the middle-east of China, has more than 0.9 million hectares of bamboo forest, being one of the top-3 Chinese provinces with richest bamboo forest resource. Such abundant bamboo forest resources in Zhejiang is mainly attributed to high level of management and developed bamboo industry. Therefore, it is said that “World bamboo forests focus on China, and Chinese bamboo forests focus on Zhejiang”.

Bamboo forest is an important part of subtropical ecosystem. Previous studies reported that bamboo forest, e.g., Moso Bamboo, has higher carbon dioxide sequestration than other subtropical forests [7,9]. With much effort focused on estimating carbon stocks of bamboo forest [10,11,12,13,14,15,16,17,18,19], it has showed that bamboo forest played an important role in coping with climate change [20,21,22,23,24,25,26].

Remote sensing data has been widely used for estimates of AGC in forests [27,28,29,30,31], and in recent years, high-resolution satellite imageries were used due to the development of technologies [32,33,34,35,36,37,38,39], e.g., the multi-resolution remote sensing imagery satellites [40,41,42,43,44,45]. Although the application of high spatial and spectral resolution sensors succeed in AGC estimation, there are still limitations such as high acquisition costs, small area coverage, multicollinearity, limited availability and a narrow bandwidth [46]. The above limitations lead to more workload and lower efficiency in estimating AGC of large-scale forest. Researches of the spatiotemporal evolution and the carbon storage estimation of bamboo forest by using medium-resolution remote sensing data has achieved plenteous results [47,48,49,50,51,52]. Meanwhile, algorithms for estimating carbon stocks/biomass based on remote sensing data have also been developed. In early studies, carbon stocks/biomass estimation was mostly conducted with linear or nonlinear regression models. For example, Anaya et al. [53] constructed a linear model using enhanced vegetation index (EVI) and related characteristics to estimate AGC in different vegetation types. Du et al. [49] used Partial Least Squares (PLS) method to estimate the carbon storage of bamboo forest in Anji County. Xu et al. [54] estimated the carbon storage of Phyllostachys praecox using the PLS regression. In recent years, machine learning algorithms have been widely used not only for inversion of land use changes but also for carbon storage estimation. Zhou et al. [26] used the K-Nearest Neighbor (KNN) method to estimate carbon storage in bamboo forests. Vafaei et al. [41] used Random Forests (RF), Support Vector Regression (SVR), Multi-Layer Perceptron Neural Networks (MPL Neural Net) and Gaussian Processes (GP) methods to estimate the forest AGC, and proved high accuracy of the methods. Gao et al. [55] combined Artificial Neural Network (ANN), SVR, RF, KNN, and Linear Regression (LR) methods to estimate forest AGC.

However, the previous studies mainly focused on small-scale area e.g., protected areas and countries. Estimation of spatiotemporal variation of bamboo forest carbon stocks on a large scale is of great significance for understanding the function of bamboo forest on climate change. In this study, we took Zhejiang province as a case study to establish AGC models of bamboo forest using the remote sensing information from Landsat5 TM and Landsat8 OLI data. The spatiotemporal dynamics of bamboo AGC was estimated and analyzed based on the extraction of spatial and temporal distribution of bamboo forests in four years (2000, 2004, 2008 and 2014) in Zhejiang Province. The study aimed at providing an insight into the spatiotemporal dynamics of bamboo forest carbon stocks in a long-time series at a national or global scale.

2. Materials and Methods

2.1. Study Area

Zhejiang province (Figure 1) is located in south of Yangtze River Delta in southeast coast, China (118–123°E and 27–31°N). Under subtropical monsoon climate, it has clearly demarcated seasons, suitable temperature with abundant rainfall. The annual mean temperature for Zhejiang is between 15 °C to 18 °C, and the annual precipitation varies in a range of 980–2000 mm. Zhejiang province has a rich forest resource, covering approximately 6.06 million hectares of forest land, 0.9 million hectares of which are occupied by bamboo forest [56]. Moso bamboo forest accounted for 87.22% of the total bamboo forest area, with 0.79 million hectares [57].

2.2. Dataset and Landsat TM Image Preprocessing

In this study, we selected 8 scene remote sensing images per year (2000, 2004, 2008 and 2014), and each of these 32 images covered the whole province (

Table 1. Acquisition date of the image data sets.

	Time Series	2000	2004	2008	2014
Row-Column Number		Landsat5 TM	Landsat5 TM	Landsat5 TM	Landsat8 OLI
118,039		16/01/2001	10/12/2004	24/03/2008	13/06/2014
118,040		16/01/2001	10/2/2004	24/03/2008	13/06/2014
118,041		29/12/1999	10/2/2004	24/03/2008	13/06/2014
119,039		17/09/2000	14/10/2004	05/07/2008	22/07/2014
119,040		04/11/2000	14/10/2004	06/06/2009	22/07/2014
119,041		29/06/2000	31/12/2003	28/02/2008	22/07/2014
120,039		10/10/2000	21/10/2004	16/10/2008	11/06/2014
120,040		10/10/2000	10/12/2004	16/10/2008	11/06/2014

). In order to reduce the influence of acquisition time, vegetation spectral reflectance, cloud thickness and other factors on classification results, we selected the images with cloud cover below 10%, smallest observational zenith angle, and consistent time interval.

We used the FLAASH method to make atmospheric correction for each image, with consideration of eliminating two major kinds of influence factors: (1) those such as water vapor, aerosol, bidirectional reflection and data transmission, which may influence the trend analysis and information extraction in time series [58,59]; and (2) the radiation difference between multi-temporal remote sensing data. Then, the corrected images were furtherly geometrically corrected using ground control points (GCPs) to splice the remote sensing data of Zhejiang Province [60].

2.3. Mapping of Bamboo

In this study, land use was classified into six types: bamboo, broad-leaved forest, coniferous forest, farmland, barren land, and water bodies [61]. Based on the fifth to eighth forest resource inventories (1994–1998, 1999–2003, 2004–2008 and 2009–2013) in Zhejiang Province, the spectral reflectance characteristics of samples were used to select Regions of Interest (ROIs) as the training samples for maximum likelihood classification based on visual interpretation, and the sampling data that derived from the continuous forest resource inventory data of Zhejiang Province, were used to validate the classification results.

Table 2. Validation samples in different years.

	Types	Barren Land	Water Bodies	Farmland	Broad	Coniferous	Bamboo
Years
2000		128	146	139	159	165	232
2004		128	146	139	142	152	215
2008		128	142	139	151	130	341
2014		128	137	139	153	156	139

shows the numbers of validation samples for bamboo, broad-leaved forest, coniferous forest, farmland, barren land and water bodies indifferent years.

2.4. AGC Estimation

As bamboo forest was dominated by Moso bamboo in Zhejiang province, the Moso bamboo sample plots were used for spatiotemporal estimation of carbon stocks on bamboo forest from 2000 to 2014 in this study. The numbers of bamboo forest plots were 137, 189, 203 and 139, respectively (Figure 1). Above Ground Biomass (AGB) of individual Moso bamboo was calculated based on the Equation (1) [62]:

$$\mathrm{M}\left(\mathrm{D},\mathrm{A}\right)=747.787{\mathrm{D}}^{2.771}{\left(\frac{0.148\mathrm{A}}{0.028+\mathrm{A}}\right)}^{5.555}+3.772,$$

(1)

where M, D, and A denote AGB (dry weight in Kg), DBH (cm), and age (du), respectively. For each plot, the AGB is a sum of all individual Moso bamboo AGB within the plot, and the expansion factor for the conversion from biomass to carbon for Moso bamboo forestis 0.5042 [9].

2.5. Construction of Estimation Model for Carbon Storage of Bamboo Forest

2.5.1. Setting of Remote Sensing Variables

The variables of bamboo forest AGC model consisted of three types, i.e., original band combinations, vegetation indices and band texture (

Table 3. Information of remote sensing variables.

Type	Name	Calculate Model	Abbreviation	Remarks
Band combination	TM546	band5 * band4/band6	TM546	Suitable Landsat5 TM data (2000, 2004, 2008)
	TM543	band5 * band4/band3	TM543
	TM542	band5 * band4/band2	TM542
	TM432	band4 * band3/band2	TM432
	TM321	band3 * band2/band1	TM321
	TM754	Band7 * band5/band4	TM754	Suitable Landsat8 OLI data (2014)
	TM563	Band5 * band6/band3	TM563
	TM547	Band4 * band5/band7	TM547
	TM432	Band4 * band3/band2	TM432
	TM543	Band5 * band4/band3	TM543
Vegetation Index	Normalized Difference Vegetation Index	(NIR-R)/(NIR + R)	NDVI	NIR, R, and B represent Near-Infrared Reflectivity, Red reflectivity, Blue reflectivity, and L take value for 0.5
	Difference Vegetation Index	NIR-R	DVI
	Simple Ratio Index	NIR/R	SR
	Enhanced Vegetation Index	2.5(NIR-R)/(NIR + 6R − 7.5B + 1)	EVI
	Soil-Adjusted Vegetation Index	(NIR-R) * (1 + L)/(NIR + R + L)	SAVI
Texture	Mean	${\displaystyle \sum _{\mathrm{i}=0}^{\mathrm{N}-1}{\displaystyle \sum _{\mathrm{j}=0}^{\mathrm{N}-1}\mathrm{iP}(\mathrm{i},\mathrm{j})}}$	Mean	$\mathrm{P}(\mathrm{i},\mathrm{j})=\mathrm{V}(\mathrm{i},\mathrm{j})/{\displaystyle \sum _{\mathrm{i}=0}^{\mathrm{N}-1}{\displaystyle \sum _{\mathrm{j}=0}^{\mathrm{N}-1}\mathrm{V}(\mathrm{i},\mathrm{j})}}$ $\mathrm{V}(\mathrm{i},\mathrm{j})$ is the ith row of the jth column in the Nth moving window; ${\mathsf{\mu }}_{\mathrm{x}}={\displaystyle \sum _{\mathrm{j}=0}^{\mathrm{N}-1}\mathrm{j}{\displaystyle \sum _{\mathrm{i}=0}^{\mathrm{N}-1}\mathrm{P}(\mathrm{i},\mathrm{j})}}$ ${\mathsf{\mu }}_{\mathrm{y}}={\displaystyle \sum _{\mathrm{i}=0}^{\mathrm{N}-1}\mathrm{i}{\displaystyle \sum _{\mathrm{j}=0}^{\mathrm{N}-1}\mathrm{P}(\mathrm{i},\mathrm{j})}}$ ${\mathsf{\sigma }}_{\mathrm{x}}={\displaystyle \sum _{\mathrm{j}=0}^{\mathrm{N}-1}{{(\mathrm{j}-\mathsf{\mu }}_{\mathrm{i}})}^2{\displaystyle \sum _{\mathrm{i}=0}^{\mathrm{N}-1}\mathrm{P}(\mathrm{i},\mathrm{j})}}$ ${\mathsf{\sigma }}_{\mathrm{y}}={\displaystyle \sum _{\mathrm{i}=0}^{\mathrm{N}-1}{{(\mathrm{i}-\mathsf{\mu }}_{\mathrm{j}})}^2{\displaystyle \sum _{\mathrm{j}=0}^{\mathrm{N}-1}\mathrm{P}(\mathrm{i},\mathrm{j})}}$
	Variance	${\displaystyle \sum _{\mathrm{i}=0}^{\mathrm{N}-1}{\displaystyle \sum _{\mathrm{j}=0}^{\mathrm{N}-1}{(\mathrm{i}-\mathrm{mean})}^2\mathrm{P}(\mathrm{i},\mathrm{j})}}$	Var
	Homogeneity	${\displaystyle \sum _{\mathrm{i}=0}^{\mathrm{N}-1}{\displaystyle \sum _{\mathrm{j}=0}^{\mathrm{N}-1}\frac{\mathrm{P}(\mathrm{i},\mathrm{j})}{{1+(\mathrm{i}-\mathrm{j})}^2}}}$	Homo
	Contrast	${\displaystyle \sum _{\left|\mathrm{i}-\mathrm{j}\right|=0}^{\mathrm{N}-1}{\left|\mathrm{i}-\mathrm{j}\right|}^2\left\{{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{N}}{\displaystyle \sum _{\mathrm{j}=1}^{\mathrm{N}}\mathrm{P}(\mathrm{i},\mathrm{j})}}\right\}}$	Con
	Dissimilarity	${\displaystyle \sum _{\left|\mathrm{i}-\mathrm{j}\right|=0}^{\mathrm{N}-1}\left|\mathrm{i}-\mathrm{j}\right|\left\{{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{N}}{\displaystyle \sum _{\mathrm{j}=1}^{\mathrm{N}}\mathrm{P}(\mathrm{i},\mathrm{j})}}\right\}}$	Dissi
	Entropy	$-{\displaystyle \sum _{\mathrm{i}=0}^{\mathrm{N}-1}{\displaystyle \sum _{\mathrm{j}=0}^{\mathrm{N}-1}\mathrm{P}(\mathrm{i},\mathrm{j})}}\log (\mathrm{P}(\mathrm{i},\mathrm{j}))$	En
	Angular second moment	${\displaystyle \sum _{\mathrm{i}=0}^{\mathrm{N}-1}{\displaystyle \sum _{\mathrm{j}=0}^{\mathrm{N}-1}{\mathrm{P}(\mathrm{i},\mathrm{j})}^2}}$	Sec
	Correlation	$\frac{{\displaystyle \sum _{\mathrm{i}=0}^{\mathrm{N}-1}{\displaystyle \sum _{\mathrm{j}=0}^{\mathrm{N}-1}(\mathrm{i},\mathrm{j})}}{\mathrm{P}(\mathrm{i},\mathrm{j})-\mathsf{\mu }}_{\mathrm{x}}{\mathsf{\mu }}_{\mathrm{y}}}{{\mathsf{\sigma }}_{\mathrm{x}}{\mathsf{\sigma }}_{\mathrm{y}}}$	Corr

). Five different window sizes (3 × 3, 5 × 5, 7 × 7, 9 × 9, and 11 × 11) were set for calculation of texture variables. Due to different numbers between Landsat5 and Landsat8 bands, the number of variables was 250 in 2000, 2004 and 2008, while 290 in 2014. When the sample plot coordinates were matched with the pixel values of remote sensing variables, they might not fully matched with each other due to geometric correction and positioning error. In order to reduce the matching error, a window size of 3 × 3 pixels was used to extract the mean values of the selected remote sensing variables for each plot [63,64,65,66].

2.5.2. Method of Model Construction

Stepwise regression screening variables method, one of the most widely used methods in regression models [67,68], was used to establish a remote sensing information model of AGC in bamboo forests. 70% of the sample plots were randomly selected for developing the model, and the others for evaluating the established model. AGC may have high or weak relationships with remote sensing variables. Because of the strong correlations among some explanatory variables, it was critical to eliminate the variables that have a high correlation between themselves and nonsignificant correlations between variables and AGC [55,69]. The advantage of stepwise regression is to determine the importance of explanatory variables and eliminate the influence of collinearity on accuracy of models. Correlation analysis can be used to examine the relationship between AGB and remote sensing variables. The basic idea of stepwise regression is to introduce the variables one by one into the model. After each of the explanatory variables was introduced, the F-test must be conducted and the explanatory variables that have been selected must be t-test one by one. When the originally introduced explanatory variable becomes less significant due to the introduction of later explanatory variables, it was deleted to ensure that the regression equation contained only significant variables before each new variable was introduced. The process did not stopped until there was no significant explanatory variable to choose the regression equation, and no significant explanatory variables were excluded from the regression equation. After applying stepwise regression, the obtained explanatory variables were optimal, and there was no serious collinearity among variables.

2.5.3. Model Evaluation

The model evaluation indexes mainly included Relative Error (RE), Mean Relative Error (MRE) and root mean square error (RMSE), as well as the analysis and evaluation of the extreme values of the model predictive value. Formulas were listed as follows, where i represents the ith sample:

$$\mathrm{R}{\mathrm{E}}_{\mathrm{i}}=\frac{\mathrm{Obs}\_\mathrm{AG}{\mathrm{C}}_{\mathrm{i}}-\mathrm{P}\mathrm{r}\mathrm{e}\_\mathrm{A}\mathrm{G}{\mathrm{C}}_{\mathrm{i}}}{\mathrm{O}\mathrm{b}\mathrm{s}\_\mathrm{A}\mathrm{G}{\mathrm{C}}_{\mathrm{i}}},$$

(2)

$$\mathrm{MRE}=\frac{1}{\mathrm{n}}{\displaystyle \sum \left|\mathrm{R}{\mathrm{E}}_{\mathrm{i}}\right|\times 100\%},$$

(3)

$$\mathrm{RMSE}=\sqrt{\frac{1}{\mathrm{n}}{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}{(\mathrm{O}\mathrm{b}\mathrm{s}\_\mathrm{A}\mathrm{G}{\mathrm{C}}_{\mathrm{i}}-\mathrm{P}\mathrm{r}\mathrm{e}\_\mathrm{A}\mathrm{G}{\mathrm{C}}_{\mathrm{i}})}^2}},$$

(4)

where $\mathrm{O}\mathrm{b}\mathrm{s}\_\mathrm{A}\mathrm{G}{\mathrm{C}}_{\mathrm{i}}$ represent the observed AGC of the ith sample, $\mathrm{P}\mathrm{r}\mathrm{e}\_\mathrm{A}\mathrm{G}{\mathrm{C}}_{\mathrm{i}}$ represent the predicted AGC of the ith sample, n represent the number of sample, respectively.

3. Results

3.1. Spatiotemporal Distribution of Bamboo

The accuracies of bamboo forest remote sensing information extraction in 2000, 2004, 2008 and 2014 are shown in

Table 4. Accuracy of classification and bamboo forest extraction in Zhejiang Province.

Year	Overall Accuracy of Land Use Classification		Bamboo Forest Classification Accuracy		Bamboo Area Estimation Accuracy
	Accuracy (%)	Kappa Coefficient	Producer’s Accuracy (%)	User’s Accuracy (%)	(%)
2000	85.04	0.82	75.86	94.12	96.50%
2004	81.59	0.78	76.28	91.62	97.50%
2008	76.26	0.75	79.18	95.07	97.50%
2014	81.69	0.78	79.41	93.1	98.90%

. Table 4 shows that the overall accuracy of land use types at different times is above 76.26%, kappa coefficient is higher than 0.75. The overall classification accuracy is high. Producer’s accuracy of bamboo forest is above 75.86%, and the user’s accuracy is above 91.62%. In addition, the accuracy of area extraction is over 96.50%, which is a satisfying result according to the actual area of the forest management inventory. Figure 2 shows the spatial distribution of bamboo forests in Zhejiang Province. According to Figure 2, the area of bamboo forests in Zhejiang Province showed a gradually increasing trend in time and space from the year of 2000 to 2014.

Based on the time-series Landsat data, the distribution information of bamboo forests in Zhejiang Province was extracted, which exhibited high accuracies in terms of both classification results and area statistics. The spatiotemporal distribution characteristics of bamboo forests are consistent with the actual situation, which provide a more accurate data.

3.2. AGC Model of Bamboo

Estimation models of AGC from the year of 2000 to 2014 were obtained by the stepwise regression method (Formulas (5)–(8)):

AGC model of 2000:

$$\begin{array}{ll}\mathrm{A}\mathrm{G}{\mathrm{C}}_{2000}& =20.378-2.660\ast {\mathrm{W}}_9{\mathrm{b}}_1\mathrm{E}\mathrm{n}+1.495\ast {\mathrm{W}}_{11}{\mathrm{b}}_6\mathrm{M}\mathrm{e}\mathrm{a}\mathrm{n}-2.415\ast {\mathrm{W}}_5{\mathrm{b}}_2\mathrm{M}\mathrm{e}\mathrm{a}\mathrm{r}\\ & +2.153\ast {\mathrm{W}}_7{\mathrm{b}}_3\mathrm{C}\mathrm{o}\mathrm{n}-1.694\ast {\mathrm{W}}_7{\mathrm{b}}_4\mathrm{E}\mathrm{n}-1.402\ast {\mathrm{W}}_3{\mathrm{b}}_6\mathrm{C}\mathrm{o}\mathrm{r}\mathrm{r}\\ & +1.487\ast {\mathrm{W}}_7{\mathrm{b}}_2\mathrm{E}\mathrm{n}-1.107\ast {\mathrm{W}}_9{\mathrm{b}}_5\mathrm{E}\mathrm{n}\end{array},$$

(5)

AGC model of 2004:

$$\mathrm{A}\mathrm{G}{\mathrm{C}}_{2004}=-3.278+15.764\times \mathrm{S}\mathrm{A}\mathrm{V}\mathrm{I}-1.984\ast {\mathrm{W}}_{11}{\mathrm{b}}_3\mathrm{V}\mathrm{a}\mathrm{r},$$

(6)

AGC model of 2008:

$$\begin{array}{ll}\mathrm{A}\mathrm{G}{\mathrm{C}}_{2008}& =-8.112+14.312\ast \mathrm{S}\mathrm{A}\mathrm{V}\mathrm{I}+24.266\ast {\mathrm{W}}_{11}{\mathrm{b}}_5\mathrm{S}\mathrm{e}\mathrm{c}+6.455\ast {\mathrm{W}}_7{\mathrm{b}}_2\mathrm{V}\mathrm{a}\mathrm{r}\\ & +0.208\ast {\mathrm{W}}_3{\mathrm{b}}_5\mathrm{M}\mathrm{e}\mathrm{a}\mathrm{n}-4.788\ast {\mathrm{W}}_7{\mathrm{b}}_2\mathrm{C}\mathrm{o}\mathrm{n}\end{array},$$

(7)

AGC model of 2014:

$$\begin{array}{ll}\mathrm{A}\mathrm{G}{\mathrm{C}}_{2014}& =-4.436+10.597\ast \mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}+6.098\ast {\mathrm{W}}_3{\mathrm{b}}_7\mathrm{S}\mathrm{e}\mathrm{c}+1.99\ast {\mathrm{W}}_3{\mathrm{b}}_5\mathrm{E}\mathrm{n}+\\ & 0.016\ast \mathrm{T}{\mathrm{M}}_{547}+1.917\ast {\mathrm{W}}_{11}{\mathrm{b}}_2\mathrm{C}\mathrm{o}\mathrm{r}\mathrm{r}+3.373\ast {\mathrm{W}}_9{\mathrm{b}}_7\mathrm{C}\mathrm{o}\mathrm{r}\mathrm{r}\end{array}$$

(8)

Here, C is bamboo carbon storage; Wi for texture window, i = 3, 5, 7, 9, 11; b_i is i band, i = 1, 2, …, 7, and the band 7 was unique for 2014.

3.3. Accuracy Assessment of AGC Model

Figure 3 shows the correlation between the predicted AGC and the observed AGC for the model. All of the correlation coefficients R of training and testing data values range from 0.63 to 0.72 in different years (Figure 3; p < 0.01), and Mean Relative Error (MRE) is less than 0.377, and the highest accuracy of the AGB model was achieved in the year of 2014 (R = 0.72, RMSE = 2.9, MRE = 0.314). Both the model accuracy and the verifying accuracy pointed to good performance of the models.

Further analysis were statistics of model predicting ability (

Table 5. Statistical analysis of the model prediction and residual.

Year	Index	Minimum Value	Maximum Value	Average Value	STD
2000	Predicted value	3.611	12.950	6.933	2.253
	Residual	−3.395	3.600	0.369	1.888
2004	Predicted value	5.106	15.546	10.679	2.486
	Residual	−3.715	3.504	0.317	1.910
2008	Predicted value	5.493	18.035	11.493	2.813
	Residual	−4.570	5.009	0.059	2.421
2014	Predicted value	6.556	16.758	12.966	2.220
	Residual	−3.573	3.113	0.122	1.824

). According to the model predicted value, (1) the maximum, minimum and average value of model predicted value exhibited increasing trends over time from 2000 to 2014. As the economic benefit increasing of bamboo forest, and the management level of bamboo forests in Zhejiang Province has been improved year by year. Meanwhile, the biomass has also increased [70,71]. Therefore, the predicted value of AGC maximum, minimum and mean value are consistent with the actual situation; (2) Predicted STD and MRE value of bamboo forest AGC are relatively stable from 2.2 to 2.9, which indicates that the model is stable and has good performance in predicting the spatiotemporal distribution of AGC in bamboo forest.

3.4. AGC Spatiotemporal Evolution of Bamboo Forest

According to the AGC model of bamboo forest in Zhejiang Province, the spatial distribution of AGC in different periods from 2000 to 2014 is shown in Figure 4. As is known, the carbon density of bamboo Zhejiang province has been continuously increasing over time, especially in Hangzhou, Shaoxing, Quzhou, Ningbo, Lishui, Jinhua, and Wenzhou.

Table 6. The Spatial and Temporal Changes of Aboveground Carbon Storage in Bamboo Forest from 2000 to 2014 in Zhejiang Province.

City	2000			2004			2008			2014
	Bamboo Area (ha)	Carbon Density (Mg·ha−1)	Total Carbon Stock (Tg C)	Bamboo Area (ha)	Carbon Density (Mg·ha−1)	Total Carbon Stock (Tg C)	Bamboo Area (ha)	Carbon Density (Mg·ha−1)	Total Carbon Stock (Tg C)	Bamboo Area (ha)	Carbon Density (Mg·ha−1)	Total Carbon Stock (Tg C)
Hangzhou	106,772.76	7.62	0.81	136,085.36	11.62	1.58	175,642.38	17.94	3.15	157,954.05	21.98	3.47
Huzhou	88,271.10	7.99	0.70	97,067.25	10.92	1.06	104,735.63	19.10	2.00	104,331.98	22.44	2.34
Jiaxing	3.87	4.69	0.00	15.66	8.77	0.00	634.41	5.90	0.00	107.73	5.29	0.00
Taizhou	84,108.87	6.12	0.51	55,287.18	8.72	0.48	21,360.78	13.49	0.29	36,291.87	16.95	0.62
Shaoxing	30,607.29	6.91	0.21	21,675.96	9.98	0.22	41,277.42	15.33	0.63	75,569.90	17.14	1.30
Quzhou	72,458.82	7.19	0.52	50,822.19	13.09	0.67	74,815.83	15.70	1.17	72,991.44	21.37	1.56
Ningbo	91,198.08	5.70	0.52	58,383.90	10.15	0.59	62,439.17	13.89	0.87	88,918.11	17.29	1.54
Lishui	103,576.95	6.15	0.64	173,737.13	11.95	2.08	163,593.90	15.61	2.55	166,664.07	18.77	3.13
Jinhua	75,212.46	6.62	0.50	68,752.17	12.23	0.84	86,205.24	13.57	1.17	75,075.39	16.18	1.21
Zhoushan	12,782.16	4.37	0.06	3560.67	8.99	0.03	5068.85	11.89	0.06	5327.37	9.83	0.05
Wenzhou	104,261.31	6.42	0.67	114,684.30	7.59	0.87	76,519.65	10.64	0.81	97,825.57	17.60	1.72

shows the statistical results of bamboo forest AGC at different periods in 11 cities of Zhejiang. According to the statistics, the carbon density in Zhejiang Province increased from 6.75 Mg·ha⁻¹ in 2000 to 19.07 Mg·ha⁻¹ in 2014, with a growth rate of 182.52%. Carbon storage increased from 5.14 Tg in 2000 to 16.94 Tg in 2014.

From the spatiotemporal distribution of bamboo forest carbon stock in 11 cities, the cities with the highest and lowest increase in total carbon storage in four periods are Hangzhou and Jiaxing. The carbon storage in Hangzhou City increased from 0.81 Tg C in 2000 to 3.47 Tg C in 2014, more than quadrupled with an average annual increase of 0.19 Tg C; Jiaxing City has low carbon storage because of the small bamboo area; otherwise, Zhoushan City has maintained a relatively stable total carbon stock. At the same time, the carbon density of bamboo forests in all counties and cities maintained a continuous growth pattern at different periods (Table 6), but the difference in growth ranged greatly. Among them, the growth of carbon stock density in Huzhou and Hangzhou was the most obvious, from 7.99 And 7.62 Mg·ha⁻¹ in 2000 increased to 22.44 and 21.98 Mg·ha⁻¹ in 2014, a threefold increase. Jiaxing City, on the other hand, showed the smallest increase from 0.6 to 2.4 Mg·ha⁻¹ from 2000 to 2014.

The relationships between the Observed AGC and several factors established the models in different years are shown in Figure 5, Figure 6, Figure 7 and Figure 8. As is seen, W₉B₁En had the highest coefficient (R = 0.3608) with the Observed AGC in 2000, followed by W₉B₅En (R = 0.2209), W₅B₂Mean (R = 0.1794), W₇B₄En (R = 0.1676), W₇B₃Con (R = 0.1658), W₁₁B₆Mean (R = 0.1360), W₃B₆Corr (R = 0.1058) and W₇B₂En (R = 0.0824). Meanwhile, SAVI (R = 0.6368, R = 0.5790), NDVI (R = 0.4007) had the highest coefficient with the Observed AGC in 2004, 2008 and 2014, respectively.

The role of the predictive variable for the AGC estimation in this study of different years was assessed using factor analysis method in the SPSS.20 software (Statistical Product and Service Solutions, SPSS; International Business Machines Corporation, IBM; Chicago, America). Accordingly, the correlation coefficient was used to calculate the merit of the variables. The result is listed in

Table 7. The importance of the variables for the AGC estimation in different years.

Year	Merit Value	Variable	Ranking
2000	0.825	W7B3Con	1
	0.813	W9B1En	2
	0.81	W7B2En	3
	0.651	W9B5En	4
	0.468	W5B2Mean	5
	0.261	W11B6Mean	6
	0.532	W7B4En	7
	0.044	W3B6Corr	8
2004	0.869	SAVI	1
	0.766	W11b3Var	2
2008	0.823	W7B2Con	1
	0.764	W7B2Var	2
	0.643	SAVI	3
	0.611	W11B5Sec	4
	0.096	W3B5Mean	5
2014	0.663	NDVI	1
	0.656	W3b7Sec	2
	0.231	W11b2Corr	3
	0.23	TM547	4
	0.11	W3b5En	5
	0.016	W9b7Corr	6

. It can be seen that W₇B₃Con, SAVI, W₇B₂Con, and NDVI were the most important variables for predicting the bamboo forest AGC in different years. Meanwhile, we could realize that the most important variables in different years were derived from vegetable indexes and texture, and they could improve the accuracy. The results are inconsistent with the recent study reported by Vafaei et al. [41] and Eckert [37].

4. Discussion

The study shows that the stepwise regression method of Zhejiang province AGC spatiotemporal estimation of remote sensing information model has a good performance. The predicting ability of the model is strong. Errors such as RMSE and predicted error are small. Figure 5 shows the analysis of the residual distribution of the model prediction residual error when the standardization residual threshold of the test sample was 2. As is shown in Figure 9, the STD of all test samples were within the threshold range, which further illustrated that the model had good stability and reliability in predicting AGC of bamboo forests at the provincial scale.

The results shows that the bamboo forest AGC and carbon density both had an increasing trend from 2000 to 2014 in Zhejiang province. The study has great consistency with the previous researches by plot sampling [9] or ecological process model simulations [72] (

Table 8. The results were compared on the AGC of Bamboo Forest in Zhejiang Province.

Study Area	Year	Carbon Density (Mg·ha−1)		Carbon Stock (Tg C)
		The Study	References	The Study	References
Zhejiang province	2004	10.95	12.17; [72]	8.42	9.45; [72]
			13; [74]
	2008	15.72	13.85; [72]	12.72	11.33; [72]
	2014	19.07	18.15; [72]	16.84	16.41; [72]

, Except for the lack of data on related studies in 2000). Certainly, there are some discrepancies mainly coursed by different methods or models of estimation. The remote sensing information model is a kind of spatial information model [73]. There are some differences on spectrum, texture, and vegetation index in a same background object that may influence the estimation results; however, the remote sensing has obvious advantages in large-scale dynamic monitoring. The results of this study and previous studies provide a guarantee for using the remote sensing information to accurately monitor the spatiotemporal dynamics of bamboo carbon storage in a wide range.

The previous researches showed that there was obvious relationship between the trend of gradual increase of bamboo forest carbon stocks and carbon density, bamboo forest area gradual expansion and improving management level of Zhejiang province [51,71]. Based on Table 6, Figure 10 showed the relationship between spatiotemporal changes in carbon stocks and carbon densities and bamboo forest area in 2000–2004, 2004–2008, 2014–2008 and 2000–2014 in different counties and cities in Zhejiang Province. As we can see from Figure 10, there is a significant linear correlation between the change of bamboo forest carbon storage and bamboo forest area. The main reason why there was no obvious linear relationship between the change of carbon density and the area of bamboo forest was that average carbon storage per unit area was affected by management level. As is shown in Table 6, the general change of bamboo forest area in Huzhou was small, especially the year of 2008–2014, but the carbon density of bamboo forest remained relatively high, and the highest of 2014 was 22.44 Mg C/ha, which was related to the long-term intensive management in the region and the high biomass of bamboo forests [75,76,77,78]. Therefore, the management level of bamboo forests in different regions will also affect the spatial and temporal evolution of carbon storage to a huge extent.

In this paper, the estimation accuracy R of carbon storage models for different years range from 0.63 to 0.69 compared with previous studies. For instance, Shang et al. [22] used MODIS images to estimate bamboo forest AGC and reported an R of 0.68. Shang et al. [47] estimated the AGC of Moso bamboo forests in combination with Landsat and MODIS data, and the estimation accuracy R was 0.70. Zhu et al. [79] estimated forest biomass using time series Landsat data and the model accuracy was 0.69. Sandra Eckert [37] used WorldView-2 data estimated forest AGB and model accuracy R is up to 0.93. Vafaei et al. [41] used ALOS-2 PALSAR-2 and Sentinel-2A data to estimate above-ground biomass in forests, with a highest estimated accuracy R = 0.85. The estimation accuracy in this study is similar to those studies based on medium resolution images, and it was less accurate than models based on high-resolution images. In this research, the satellite remote sensing image used to estimate the carbon stocks were Landsat time series of medium resolution data with a spatial resolution of 30 × 30 m, and there were many disadvantages compared to the high-resolution satellite data used to estimate forest AGC [80]. Previous studies showed that the combination with multi-source remote sensing data could effectively improve the estimation accuracy of AGC [81]. Remote sensing data were affected by their own spectral resolution, resulting in differences in the extraction accuracy of the band spectrum, vegetation index, texture information compared with high-resolution images [37,38,42,81]. In addition, due to the large span of date and time acquisition of the eight scenes in the same period, although atmospheric correction was performed, radiation differences could not be completely eliminated, resulting in errors in bamboo forest information extraction [82]. This is what we might improve in the future study.

5. Conclusions

Based on the time series imageries of Landsat 5 TM and Landsat 8 OLI in 2000, 2004, 2008 and 2014, we took Zhejiang province as study area, as a precondition to extract bamboo forest distribution information in different periods of Zhejiang province. Model of remote sensing variables were constructed to estimate the spatiotemporal evolution of bamboo forest AGC in Zhejiang province. The results shows that:

(1)

The spatiotemporal distribution of bamboo forests in Zhejiang Province at different periods had a higher accuracy of information extraction, of which the classification accuracy reached above 76.36%, the user’s accuracy was above 91.62% and the area accuracy was above 96.50%.

(2)

Bamboo forest AGC spatiotemporal estimation model built by the stepwise regression method in Zhejiang Province has good performance and robustness. RMSE and prediction error are small. The estimated carbon storage results have a good consistency with the previous research.

(3)

Bamboo forest AGC storage shows gradually increased tend in Zhejiang province from 2000 to 2014, and the average carbon stock density at different years was 6.75 Mg·ha⁻¹, 10.95 Mg·ha⁻¹, 15.25 Mg·ha⁻¹ and 19.07 Mg·ha⁻¹, and an average annual growth was 0.88 Mg·ha⁻¹. Spatiotemporal evolution of bamboo forest carbon stocks has close relationships with the expansion of bamboo forest area and the differences in management level in various regions of Zhejiang province.