Use of vegetation change tracker, spatial analysis, and random forest regression to assess the evolution of plantation stand age in Southeast China

Abstract

Key message

By integrating vegetation change tracker (VCT), spatial analysis (SA), and random forest regression (RF), the spectral-temporal patterns of forest stand age were mapped for three typical plantations in Southern China. The spectral-temporal distribution of age structure indicated that the plantation stands in the study area were increasingly aging.

Context

Plantations play a major role in China for ecosystem restoration and carbon sequestration. Mapping plantation stand age distributions is essential for developing sustainable plantation forest management plans.

Aims

The purpose of this study was to propose two new remote sensing based models for mapping plantation ages and to test the model feasibility and accuracy in determining the spectral-temporal patterns of forest ages.

Methods

We first integrated vegetation change tracker (VCT) algorithm and spatial analysis (VCT-SA) for the pixels that were disturbed at least once from 1987 to 2017, and integrated VCT and random forest (VCT-RF) for the pixels were not disturbed during the study period. Then the forest age of these two parts were merged separately to generate annual forest age.

Results

The spectral-temporal (30-m resolution, from 1987 to 2017) of forest age for the three typical plantations in Lishui were generated. The results indicated that the plantation stands in our large study area were increasingly aging.

Conclusion

Our results revealed that it is reasonable to derive the distribution of plantation stand ages from combined remote sensing models. Besides, we confirmed that the stand ages of the plantations in our large study area of Lishui City are on the rise as the result of forest protection policies.

1 Introduction

Forests play vital roles in environmental improvement, landscape aesthetics, the mitigation of climate change, and determination of energy budgets (Goulden et al. 2011; Meijaard et al. 2013; Pan et al. 2011a; Poorter et al. 2015). As discussed by these authors, stable forests have the capacity to meet the social, economic, ecological, cultural, and spiritual needs of current and future generations in a sustainable manner. Currently, natural forests are decreasing, while plantation areas are increasing on a global scale. By the end of 2018, China had 80 million ha of plantations, which accounted for approximately 33% of the world’s total area of tree plantations and represented 36% of the forested area in China (SFA 2019). Plantation stands are anticipated to increasingly contribute to the world’s supply of wood, fibre, fuel, and non-timber forest products, as well as the provision of environmental and social services (Liu et al. 2014). Particularly, to meet the demands of timber production, short rotation and fast-growing plantations were principally created in Southern China due to its abundant precipitation, ample sunlight, and excellent soil fertility conditions (Zhou et al. 2011). Forest age, which reflects the past disturbance legacy, is an essential factor for quantifying harvesting potentials (Verkerk et al. 2011) and long-term forest carbon sequestration potentials (Poorter et al. 2016; Zhang et al. 2017). It can also be employed as an indicator for the assessment of biodiversity, recreation attractiveness, and disturbance risks (Payette and Frégeau 2019; Seidl et al. 2018). In Southern China, forest management activities have resulted in shifts of forest ages via frequent forest logging and post-harvesting recovery events. Therefore, reliable historical forest age structures are required to improve the modelling of forest carbon balances (Juutinen et al. 2018; Thom and Keeton 2019) and to quantify additional management measures.

Traditional methods for collecting forest parameters are based on a plot sampling design such as simple random or systematic sampling where the tree heights and diameters on the selected plots are measured using a variety of field instruments. Obviously, this manner of forest data collection is time-consuming, labour-intensive, and costly over large terrestrial regions (Almeida et al. 2019). In addition, this method is impossible or extremely costly to implement due to poor accessibility such as when measuring forests on very steep slopes in remote areas with poorly maintained or insufficient access roads (Wiggins et al. 2019). However, satellite remote sensing technique can overcome these limitations and provides a systematic means for the broad area spatially explicit estimation of biophysically essential variables. A considerable number of investigations have demonstrated that multispectral remote sensing imagery is a useful tool for discriminating the age of forests (Thom and Keeton 2019; Wunderle et al. 2007).

Conventional remote sensing based on forest age estimation approaches tend to examine spectral reflectance and texture analysis by employing a combination of shape, size, and spectral data to estimate forest as a class rather than as an exact forest age (Kayitakire et al. 2006; Nilson and Peterson 1994). The underlying physical mechanism for such approaches lies in the fact that different light reflectance and transmittance regimes of canopy structures and image texture vary from forest stands with different optical properties of species, tree branch, leaf, bark, trunk and cones (Zhou 2003). Canopy reflectance is a product of competing mechanisms of light absorption (Roberts et al. 1998). Chlorophyll concentrations increase with forest age which affect the multiple scattering of light in the canopy (Croft et al. 2014; Kokaly et al. 2003). Also, with increasing forest age, the image texture changes in forest canopies from a uniform style to a clumpy pattern. This enables the development of associations between the characteristics of canopies and age (Cohen et al. 1990; Cutler et al. 2012).

Within such conventional estimation analyses, statistical models, such as multivariate (linear) regression (Chen et al. 2012a), artificial neural networks (Chen et al. 2012b), support vector machines (Kauffman and Prisley 2016), regression tree algorithms (Shataee et al. 2012), and recently used random forests (Healey et al. 2018), have been employed to establish relationships between the stand age and optical and/or microwave-based parameters (spectral reflectance and/or vegetation indices and backscatter coefficients) derived from remotely sensed images, including MODIS (Hovi et al. 2019; Kuusinen et al. 2014; Li and Fox 2012), Quick Bird (Dye et al. 2012), Landsat (Nelson et al. 2000), SPOT HRV (Wunderle et al. 2007), SAR (Cutler et al. 2012), LiDAR (Hovi et al. 2019; Wiggins et al. 2019), Radar (Imhoff 1995), etc. Although the methods listed above may be used to map the ages of plantation stands, there are several issues that should be reconsidered. Firstly, the spectral saturation of mature forests (Donaldson Hanna et al. 2019) and higher reflection of young forests, particularly in those locations that have undergone frequent land use changes (Chen et al. 2012a; Dong et al. 2013), might reduce the accuracy of models that estimate the ages of plantation stands. Secondly, forest reflections primarily emanate from particular tree species, canopy structures and background structures. Differences in the reflectance properties between tree species are attributable to the variable optical properties of foliage elements and canopy structures (Croft et al. 2013; Donaldson Hanna et al. 2019). Specific features should be selected, and dedicated models should be developed for different forest species. Thirdly, algorithms related to forest age invariably integrate an extensive feature set, which slows down the algorithms via the inclusion of excessive resources/data points. Many machine learning algorithms exhibit a decrease in accuracy when the number of variables is significantly higher than optimal (Nilsson et al. 2007). The selection of small (potentially minimal) feature sets that generate the best possible results is desirable for practical reasons. Lastly, existing forest age estimation products are mostly coarse, with a resolution of 250 m or 500 m (Hovi et al. 2019; Kuusinen et al. 2014; Li and Fox 2012), and they are created at limited time points, or even only one time point, which are unable to capture forest disturbances and recovery histories over adequate spatial or temporal domains.

A highly automated algorithm, referred to as vegetation change tracker (VCT), has been developed to reconstruct of recent disturbance histories using Landsat time series stacks (LTSS) (Huang et al. 2010). VCT is based on spectral-temporal characteristics of land cover and forest change processes and requires minimal, or no, fine tuning for most forests with closed, or nearly closed, canopy covers. In the VCT process, the forest is one of the darkest vegetated surfaces in the satellite images acquired during the leaf-on growing season via both visible and several shortwave infrared bands. During the growing season, forests have similar integrated forest z-score values, which are substantially more stable over time, and are mostly lower than those of non-forest land cover types (Huang et al. 2010). Previous studies have shown that VCT is effective for the mapping of forest disturbances and recovery histories, with overall accuracies ranging from 70% to 86%, which was achieved for disturbances that were mapped at the individual year level in the USA (Huang et al. 2010; Li et al. 2016; Shen et al. 2017) and Canada (Pickell et al. 2014; Stueve et al. 2011). Furthermore, the disturbance detection accuracy increases to over 92% when the disturbance data was aggregated into three classes: persisting non-forest, persistent forest and disturbed forest (Thomas et al. 2011).

By combining available historical inventory data with historical Landsat images, this study aimed to develop a reliable framework to track the trajectories of plantation stand age transitions in three types of plantation stands: pine (Pinus massoniana Lamb.), Chinese fir (Cunninghamia lanceolata (Lamb.) Hook.), and oak (Cyclobalanopsis glauca (Thunb.) Oerst.) from 1987 to 2017 in the Lishui City, Zhejiang Province, China. To obtain a robust framework, two methods are proposed to reconstruct the plantation stand age trajectories in Lishui City. The first is an integrated VCT and spatial analysis (VCT-SA) for the mapping of recovery forest ages, where pixels are disturbed at least once during the study period. The other was an integrated VCT and random forest (VCT-RF) model to calculate plantation stand ages for persisting forests in 1997 and then transfer to the remaining years. Specifically, this study addressed a critical question: can we demonstrate a robust approach for mapping age patterns and their changes for plantation forests in China?

2 Material and methods

2.1 Research site and sample data collection

Lishui is a mountainous city, situated in the Southwest of Zhejiang Province, with an administrative region that spans in latitude, from 27°25′ to 28°57′N, and longitude, from 118°41′ to 120°26′E. Lishui City encompasses 17,300 km², of which 80.79% is covered with forest (Fig. 1). The dominant tree species in this region are pine, Chinese fir, and oak, which constitutes 33.83%, 26.02%, and 16.99%, respectively, of its forested land, derived from the Forest Management Inventory (FMI) of Lishui City, conducted in summer 2013. Lishui City has undergone considerable land use changes in recent decades, which were largely due to deforestation and reforestation (Liu et al. 2017). Most deforestation events occurred prior to the 1970s, since at that time wood production was the sole aim of forest management.

Fig. 1

Maps showing location of study site and forest type distributions at Lishui City, Zhejiang Province, China. P119R040 indicates the path/row 119/040 of Landsat in WRS-2, and P119R041 indicates path/row 119/041 of Landsat in WRS-2

For this study, the Forest Resources Sampling Survey (FRSS) of Zhejiang Province was implemented during the summers of 1994, 2004, and 2014, and the survey data were used to construct and validate age estimation models. The FRSS is based on a systematic sampling design (the sample spacing is a 4 km × 6 km grid) and at each intersection of the grid a permanent field plot is established to perform a repetitive full calliper investigation to estimate the age of individual trees (Shen et al. 2018). The permanent sample plots were used to estimate and validate the forest age model in our study area (Table 1). Each plot covers 800 m² (28.28 m × 28.28 m), and it is permanently marked and periodically revisited. In the current work, forest tree species and their corresponding longevity data were compiled from the FMI.

Table 1 Number of training and validation plots

2.2 Datasets and preprocessing

We employed Landsat Surface Reflectance (SR) time series products due to their free accessibility, 30-m spatial resolution, and long consistent acquisition record (https://landsat.usgs.gov/landsat-surface-reflectance-data-products). In this study, we focused on the forested lands of Lishui, Zhejiang Province, which is jointly covered by the Landsat scene with World Reference System II (WRS2) path/rows (WPR) tiles 119/040 (hereafter P119R040) and 119/041 (hereafter P119R040). Those time series Landsat SR images (Band Blue, Green, Red, NIR, SWIR1, SWIR2, and Thermal) spanning from 1987 to 2017 (Fig. 2) were selected by two criteria image date and cloud cover percent. Images were selected to be in the growing season (from May to September for midlatitude regions) and with low cloud cover (defined as having < 10% cloud cover throughout this study). The first criterion was necessary to avoid the confusion of off-leaf deciduous trees with disturbance events outside this time window (Huang et al. 2010; Li et al. 2016). Also, images acquired in April or October served as substitutions of images acquired from May to September when they were heavily cloud contaminated. Due to heavy cloud contamination, the Landsat images for 1998 and 2012 were unavailable. The thermal band was resampled to 30-m resolution and used to detect the cloud and shadows due to their distinctive temperature difference compared to other land cover types (Huang et al. 2010). Then those pixels were replaced by clear pixels at the same location in other images of the same year. By carefully examining these two criteria but not very strictly following them due to the above-mentioned constraints, all of the selected surface reflectance images were combined to form an annual LTSS, defined as a temporal sequence of Landsat images, consisting of one image in each year. Basically, the stack consists of annual Landsat images or composites falling into the growing season during the entire research period. However, if these images of the growing season of a year are heavily contaminated by clouds, we will use a clear image acquired in May or October to replace these contaminated images to create this stack. Additionally, if these images of the growing season of a year are not heavily contaminated by clouds, and the clouds in these images are not spatially coincided, we will combine these images together to filter out all the clouds to form a clear composite image, feeding into the stack.

Fig. 2

Statistics of the number of Landsat images used by a P119R040 and b P119R041. Triangle dots represent images of Landsat TM. Square dots represent images of Landsat ETM+. Round dots represent images of Landsat OLI

2.3 VCT algorithm to detect disturbances

Based on the spectral-temporal characteristics of land cover and forest change, the VCT algorithm consists of two major steps: individual image analysis and time series analysis of the LTSS (Huang et al. 2010). Conceptually, VCT firstly performs cloud and shadow masking on the imagery and then detects changes in a synthesized forest index as a possible disturbance signal. For a given LTSS, this algorithm initially identifies forest samples in each image based on an integrated forest z-score (IFZ) that is calculated for each pixel. Pixels are persisting forest when their IFZ values are less than 3 (Huang et al. 2010). The algorithm then tracks the IFZ changes over time to detect where and when a forest disturbance(s) occurs and produces scene-level annual forest disturbance maps featuring a disturbance class (disturbed in this year) and five other classes (persisting forest, probable forest with recent disturbance, persisting non-forest, persisting water, and post-disturbance non-forest) (Huang et al. 2010). Forest disturbance is referred to as an event resulting in the sudden reduction or removal of forest canopy cover and woody biomass and is often manifested by abrupt spectral changes (Huang et al. 2010). The majority of forest disturbance is caused by frequent forest harvests in our study area, which could be captured by the VCT algorithm. After disturbance, the disturbed patch starts to recover naturally or artificially. A couple of years later, once this patch looks spectrally like a forest, the VCT algorithm will regard this patch as a successful forest regeneration event. At this point, the VCT declares that this patch spectrally resembles a forest. VCT comprises an automated algorithm that is designed for the evaluation of forest disturbances and post-disturbance recovery processes using spectral-temporal signals recorded during Landsat time series observations (Huang et al. 2010). We produced scene-level annual forest disturbance maps separately for the two Landsat scenes (P119R040 and P119R041) for our study forestlands. Next, the annual disturbance images of the two scenes were merged and then extracted by the administrative vector boundary of Lishui city to produce the annual forest disturbance map. Then we aggregated the “persisting forest”, “probable forest with recent disturbance”, as well as “disturbance in this year” into the “forest” class to derive estimates for the total annual forest cover. Subsequently, we extracted the disturbed pixels of each year from the disturbance map and calculated the forest disturbance rate by dividing the number of disturbed pixels by forest pixels.

2.4 Mapping plantation stand age

In accordance with the terminology of the VCT algorithm, forested pixels that were previously disturbed by forest management (primarily harvesting activities) represent the “VCT post-disturbance regrowing forest”. Forest pixels that did not experience recent disturbances during the time window of the Landsat LTSS are called “VCT undisturbed forest”. As shown in Fig. 3, we employed the VCT algorithm to initially detect these disturbances and then allocated the corresponding pixels into two classes to determine their ages separately: (1) integrating vegetation change tracker (VCT) algorithm and spatial analysis (VCT-SA) for the pixels that were disturbed at least once from 1987 to 2017 and (2) integrating VCT and random forest (VCT-RF) for the pixels were not disturbed during the study period. Lastly, the annual forest age of the two kinds of forest were merged as the annual forest age of Lishui city.

Fig. 3

Flow chart for the combined models

In the VCT-SA model, forest age was calculated by subtracting the disturbance year from its corresponding recovered year (when the disturbed pixels look spectrally like a forest on the image), which both could be identified from the LTSS. Estimated age is then modified by the addition of another 2 years (since reforestation generally occurs in the next year following disturbance by using 1-year-old nursery seedlings in Southern China). For example, when we map the age of post-disturbance regrowing forest for 2017 (Fig. 3), all the disturbed pixels from 1987 to 2016 were overlaid with the probable forest pixels in 2017. If a pixel was a disturbed pixel in 1987 and evolved to a probable forest pixel in 2017, then this pixel age should be 32 years (2017 minus 1987 and then plus 2). Subsequently, the forest pixels along with the estimated forest age were separately extracted from the pine, Chinese fir, and oak plantation samples and then validated for each type of plantation by forest age of FRSS sample plot.

The spectral characteristics such as canopy reflectance and texture change gradually with stand age at the early age (< 40 years) and become stable at a later ages (> 40 years) (Croft et al. 2014). Hence, the sample data with the forest age less than 40 years and the highest quality Landsat images preferably during the growing season that covered the study area were considered to build the VCT-RF model. Thus, we chose to map forest stand age of 1997 by the sample data and the Landsat images with the highest quality in the time series. Subsequently, corresponding indices of the Landsat image in 1997 (including vegetation index, textures, and tasselled cap transformations) were extracted to relate to the selected forest ages. These vegetation indices have been proven to be effective for mapping forest development stages and forest stand ages (Dye et al. 2012; Jensen et al. 1999; Johansen et al. 2007). Details of the spectral and textural measures used in the age model are summarized in Table 2. Based on these features, random forest-recursive feature elimination (RF-RFE) with tenfold cross-validation was applied for feature selection using training samples (Almeida et al. 2019). Following this, the best model was identified and applied to calculate the forest age for 1997 from which the remaining years in the LTSS were derived, by subtracting/adding the year(s) from/to this year.

Table 2 Spectral and textural features and their corresponding equations used in this study

2.5 Uncertainty of model

The Pearson correlation coefficient (COR) measures the strength of linear association between the predictions and observations (Sedgwick 2012). Thus, the COR was used to separately test the degree of uncertainty of VCT-SA and VCT-RF for the three plantations. Root mean square error (RMSE) was used to measure the differences between the predicted values of model and those observed. Modelling efficiency (MEF) reflects the ability of the model to project the data by measuring correlations between the predictions and observations, as well as their coincidence, such as the deviation from the 1:1 line (Bowman and Azzalini 1997). The bias (BIAS) of the model reveals whether the model has a systematic bias (Janssen and Heuberger 1995). It was applied to estimate the mean relative differences between the predictions and observations.

$$ \mathrm{COR}=\frac{\sum_{\mathrm{i}=1}^{\mathrm{N}}\left({\mathrm{P}}_{\mathrm{i}}-\overline{\mathrm{P}}\right)\left({\mathrm{O}}_{\mathrm{i}}-\overline{\mathrm{O}}\right)}{\sqrt{\sum_{\mathrm{i}=1}^{\mathrm{N}}{\left({\mathrm{P}}_{\mathrm{i}}-\overline{\mathrm{P}}\right)}^2{\sum}_{\mathrm{i}=1}^{\mathrm{N}}{\left({\mathrm{O}}_{\mathrm{i}}-\overline{\mathrm{O}}\right)}^2}} $$

(1)

$$ \mathrm{RMSE}=\sqrt{\sum_{\mathrm{i}=1}^{\mathrm{N}}\raisebox{1ex}{${\left({\mathrm{P}}_{\mathrm{i}}-\overline{\mathrm{O}}\right)}^2$}\!\left/ \!\raisebox{-1ex}{$\mathrm{N}$}\right.} $$

(2)

$$ \mathrm{MEF}=1-\frac{\sum_{\mathrm{i}=1}^{\mathrm{N}}{\left({\mathrm{P}}_{\mathrm{i}}-{\mathrm{O}}_{\mathrm{i}}\right)}^2}{\sum_{\mathrm{i}=1}^{\mathrm{N}}{\left({\mathrm{O}}_{\mathrm{i}}-\overline{\mathrm{O}}\right)}^2} $$

(3)

$$ \mathrm{BIAS}=\frac{1}{\mathrm{N}}{\sum}_{\mathrm{i}=1}^{\mathrm{N}}\raisebox{1ex}{$\left({\mathrm{P}}_{\mathrm{i}}-{\mathrm{O}}_{\mathrm{i}}\right)$}\!\left/ \!\raisebox{-1ex}{${\mathrm{O}}_{\mathrm{i}}$}\right. $$

(4)

where N is the number of paired values, P_I and O_i are the individual predictions and observations, whereas $\overline{\mathrm{P}}$ and $\overline{\mathrm{O}}$ represent the means of the predictions and observations, respectively.

3 Results

3.1 Annual VCT disturbance results

Figure 4 a illustrates an example of forest disturbance detected by the VCT algorithm. When calculating the percentage of two types of forest pixels, about 67.65% of the total forested area that was designated as VCT undisturbed forests, while the remaining 32.35% of the forested area was defined as VCT post-disturbance regrowing forests. Spatiotemporal patterns of forest disturbances detected by the VCT algorithm revealed where and how many forest areas were disturbed in earlier years in Lishui City (Fig. 4b). The annual forest disturbance rates were produced based on the annual forest disturbance map from 1987 to 2017. As shown in Fig. 5, the disturbance rate ranged from a low rate of 0.1% in 2002, to a high of rate 4.6% in 2008, except for the years 1998 and 2012. Most rates were ranging from 0.5% to 2.0%, which represented the normal forest disturbance levels for the study area.

Fig. 4

a Forest disturbance and restoration product developed by VCT algorithm in 2017. b Forest disturbance patterns mapped by the VCT algorithm during the 1987–2017 time period. For each year, the VCT algorithm was used to produce maps of disturbance. Dark green pixels representing “persisting forest” in the VCT algorithm were defined as “VCT undisturbed forest”; light green pixels representing “probable forest with recent disturbance” were defined as “VCT post-disturbance regrowing forest”

Fig. 5

Percentage of disturbance over time from the VCT algorithm

3.2 Results and accuracies of VCT-SA and VCT-RF

To examine the performance of these methods, we compared the VCT-SA predicted ages of the three plantations stand types against the actual ages derived from forest inventories, which were separately implemented in the summers of 1994, 2004, and 2014. For each forest type (Fig. 6), the top 21, 12, and 24 candidate variables were selected for pine, Chinese fir, and oak plantation, respectively (Table 3). The use of these variables produced the lowest RMSE compared to the use of all variables (n = 30) in the three type forests. Specifically, the decrease of RMSE was much higher in Chinese fir (0.33 years) than in the oak (0.05 years) and pine (0.07 years) when RF-RFE was applied in model building. Examining the variables selected for building the optimal model, we found that only image textures were selected for Chinese fir plantation modelling but spectral and textural variables were identified for pine and oak plantations.

Fig. 6

The RMSE obtained from the recursive variable selection (REF) method (X-axis represents number of candidates ranking by importance). The lowest RMSE (year) obtained from the lowest number of variables were 21, 12, and 24 for pine, Chinese fir, and oak forests, respectively

Table 3 Variables selected by random forest-recursive feature elimination (RFE) with tenfold cross-validation

The plots (Fig. 7) are linear regressions of predicted forest age values from VCT-SA and VCT-RF in three plantations. The observed plantation stand age revealed strong positive correlations with its predicted values. The correlation between predictions and observations of VCT-SA was higher than 0.80 (p < 0.001), while the correlation of VCT-RF (0.69~0.74) was lower. The RMSE of VCT-SA (5.95~6.84 years) was higher than that of VCT-RF (4.84~5.55 years). The MEF of VCT-SA was highest in Chinese fir plantation but lowest in pine plantation. The BIAS of all models was lower than 0.16 years, and the modelling efficiency (MEF) was higher than 0.4, except in the VCT-SA of pine plantation forest (Table 4).

Fig. 7

Relationships of observed versus predicted in VCT post-disturbance regrowing forest of pine plantation (a); VCT post-disturbance regrowing forest of Chinese fir plantation (b); VCT post-disturbance regrowing forest of oak plantation (c); VCT undisturbed forest of pine plantation (d); VCT undisturbed forest of Chinese fir (e); VCT undisturbed forest of oak (f). COR is Pearson’s correlation coefficient; RMSE is the root mean square error; MEF is modelling efficiency; BIAS is mean relative difference. Comparison of observed and predicted forest age. The diagonal line indicates the domain where predicted values equal field-measured values

Table 4 The performance validation of VCT-SA and VCT-RF in three plantations

3.3 Results of forest age mapping for Lishui City

Over the last 31 years, changes in the age structures of forests have revealed an increasing age structure (Fig. 8). During the period from 1987 to 2017, the dominant age class initially was 11~20, then 21~30, followed by 31~40 and most recently by 41~50 years in age. These four age classes accounted for about 85% of the trees in our study area in the year 2017. The oldest trees in age clasess 51~60 and 61~70 accounted for less than 15% of the trees in our study area in the year 2017. The spatial and temporal plantation stand ages from 1987 to 2017 are depicted in Fig. 9, where early in the 30-year window, the plantation stand age was dominated by trees under 20 years old. However, over time, the trees gradually grew, and the forest became dominated by an age bracket of 40~50 years. By analysing the age structure, we found that the plantation stands in our study area in 2017 were 26.55 years older than those in 1987 (the average ages of pine, Chinese fir, and oak were 15.65 years, 15.19 years, and 17.26 years in 1987 and then 42.78 years, 40.86 years, and 44.42 years in 2017, respectively), which might be attributed to the effects of forest protection policies.

Fig. 8

Annual distribution of plantation stand ages from 1987 to 2017 at Lishui City, Zhejiang Province, China. Different coloured bars represent different stand age classes (Y-axis) (refer to the legend across the top)

Fig. 9

Mapping of plantation stand age at Lishui City, Zhejiang Province, from 1987 to 2017. Seven classes of forest stand age (year) are shown in different colours (refer to the legend inset). From these smaller images, it was evident that, over time, the area covered by older forest stands increased

4 Discussion

With abundant rain and relatively fertile soils, the productivity of plantation stands in the Lishui City is among the highest in China. These plantation stands are expected to increasingly contribute to the supply of wood, fibre, fuel, and non-timber forest products, as well as the provision of environmental and social services (Liu et al. 2014). As environmental policies in China are increasingly focused on ecosystem-oriented management and protection (Zhou 2003), forest management in this area has significantly shifted from wood production only, to ecological protection. A snowstorm and freezing rain disaster struck the most populated and economically developed South Central region of China, from January 10 through February 8, 2008. Approximately 20 million ha of forest were damaged during this event through widespread crown and stem breakage, branch snapping and bending, and uprooting (Zhou et al. 2011). This can explain why there was a big jump of disturbance rate in 2008 of the LTSS (Fig. 5).

Based on the history of forest management and artificial and natural disturbances, we integrated remote sensing-based models (VCT, a spatial analysis model, random forest algorithm, as well as sample data) to map the resetting patterns of plantation stand ages resulting from forest disturbance and restoration events. In contrast to previous studies (Chen et al. 2012b; Nelson et al. 2000; Zhang et al. 2014; Zhang et al. 2017), the innovation of this study resided in the utilization of forest age products derived from the VCT algorithm. This study showed the spatial distribution of persisting forests and newly re-established forest stands (Fig. 4a). In this study, we applied VCT-SA in the post-disturbance regrowing plantation area to estimate forest ages, with the prerequisite of VCT-SA being that certain initial information on forest management should be available (Li et al. 2016). Thus, the age of seedlings transferred from nurseries and reforestation events that proceeded during the year following the disturbances should be known or reasonably presumed. Additionally, the number of observations could also affect the precision of the model. If we have observation with forest age less than 10 years in the model of VCT-SA for pine plantation, the model should perform consistent with the Oak and Chinese fir plantation types. Further, it should be noted that the initial information required may be site specific or latitude-dependent; thus, prior to transferring our VCT-SA to other regions, we needed to confirm this initial information with local authorities to ensure the accuracy of the age mapping.

Leaf chlorophyll values increase gradually with stand age and begin to plateau at a later age (40~50 years) which results in a change of canopy reflectance (Croft et al. 2014; Kokaly et al. 2003; Oester 1981). The optical properties of a tree not only depend on leaf optical properties but also on those of bark and cones especially for conifers (Guyot et al. 1989). The above studies provide the foundation for using three different VCT-RF models for specific forest stands by spectral characters for the earlier year (1997) when the majority of forest age was under 40 years in Lishui City. Multiple spectral data carries redundant information which might result in low accuracy and a high model calculation load if all relevant spectral variables were used in the modelling process. The RF-RFE is capable of selecting important spectral and texture features that are sensitive to forest age of specific forest species and improves the accuracy levels. Similarly, Pullanagari et al. (2018) and Granitto et al. (2006) used RF-RFE for analysing high-dimensional data and found it to be an efficient feature selection method, far better than traditional methods.

Li and Fox (2012) applied MODIS time series to extract young (< 4 years) and mature rubber trees (≥ 4 years) and gained an accuracy of over 97%. By integrating both Landsat time series and L-band PALSAR images, Kou et al. (2015) achieved accuracies between 80% and 90% for the three forest age groups mapping. The validation accuracies of our forest age were inferior to those results with the potential reasons as follows. Firstly, these studies only mapped age groups rather than quantitative annual ages. Conceptually, technical challenges of previous studies are not as high as our current work, i.e. mapping ages. Secondly, there was no L-band PALSAR information used in our current work, leading to the inadequate acquisition of forest structural information, although there is a strong connection between forest age and forest structure (Kou et al. 2015). Our subsequent research work may consider integrating available L-band PALSAR into plantation age modelling to improve the accuracy of mapping ages. Thirdly, there were not always sufficient available high-quality annual Landsat images acquired for Lishui City during the leaf-on season for this long time period (1987–2017). The biggest issue was the coexistence of abundant summer precipitation and cloud contamination in our study area. To include one imagery for each year, for second best, we had to use clear images acquired on April or October (Fig. 2), which is the leaf-off season heavily impacting the identification of deciduous trees (e.g. oak plantation in the current work), to replace the growing season images that had relatively high cloud contamination. Thus, the replacement may undermine the accuracy of mapping forest ages by using spectral and textural inputs from a different season (Fig. 7c and f, Table 4). Kuusinen et al. (2014) studied the relationships between forest age and albedo in the visible, near-infrared, and short-wave broadband based on MODIS and Landsat albedo retrievals for pine, spruce, and broadleaved forest and found that the stand level’s errors of the models applied to broadleaved forests were higher than those in coniferous forests. However, in our study, the RMSE of VCT-RF for oak plantation was lower than that of pine plantation. The major reason responsible for this difference may be the different instruments with different radiometric characteristics and different physical variables (Kuusinen used albedo considering semi-spherical reflections, whereas we used reflectance considering directional reflections). In Chinese fir plantation only texture variables were selected, which was different from the variables selected for the other two plantations types. The difference might be attributed to characteristics of Chinese fir trees. The first is that the forest canopy of the Chinese fir has a conical shape and often the crowns are not adjacent or overlapping (Guyot et al. 1989). Secondly, because Chinese fir is a fast-growing tree, the growth rate is higher than that of the other two tree species before 40 years of age, which results in rapid growth of DBH and tree height (Chen et al. 2015; Zhao et al. 2009). Chinese fir trees whose crowns typically do not overlap other trees and the different DBH and tree height for this species may result in the texture indices being more important than spectral and tasseled cap (TCT) indices in classifying plantation ages.

In contrast to our study, the age classes or groups maps with a spatial resolution of 250 m or 500 m derived from MODIS (Hovi et al. 2019; Kuusinen et al. 2014; Li and Fox 2012) were coarser in terms of age and spatial discrimination ability, which are less informative and undermine the practical usage of these maps in forest management practices. Besides, these studies focused on the mapping of forest ages for a specific year only, whereas our work mapped the ages of plantations and inferred the temporal transitions of ages in a time series manner. Thus, our mapping results showed good utility toward the development of reasonable forest management strategies. In addition, the premature tree growth stage (e.g. < 25 years) is typically challenging to estimate (Chen et al. 2012a; Dong et al. 2013); thus, a dedicated method should be applied for premature trees, which is distinct from mature trees. For the current work, we applied two different models to map the ages of VCT post-disturbance regrowing and VCT undisturbed forests, wherein this strategy fully considered various challenging issues when mapping the forest ages of premature forests (similar to post-disturbance regrowing forests here) and mature forests (similar to undisturbed forests).

Additionally, although Landsat data offers extensive spatial coverage, its average spatial resolution may result in errors related to the recognition of forest disturbances by VCT, in contrast to the SPOT HRV (high resolution visible) (Wunderle et al. 2007). Indeed, multiple mixed pixels in each scene may greatly complicate the identification of land features, while hindering the measurement of spectral indices. Atmospheric conditions were also of great concern as relates to object interpretation. Hence, some forests that were disturbed in the past years had incompletely recovered and do not readily present spectrally like a forest. In this case, the forest plantation age is underestimated by our model (Fig. 7). High-resolution remote sensing data with low cloud coverage is eagerly anticipated for applications in further research, such as mapping forest stand ages. Further, supplementary LiDAR derived data might be used in conjunction with Landsat data to gain further insights into forest stand age dynamics (Wiggins et al. 2019). Periodically updated information via satellite remote sensing technology could thus provide valuable data regarding changes in plantation stand age while assisting forest resource managers with the formulation of enhanced forest management plans.

For our study, we mapped plantation stand ages annually in the time series from 1987 to 2017 in Lishui City. These time-consistent forest age products by forest type can be employed by forest managers to inform forest harvesting quantities and locations in accordance with various management purposes for different tree species, as well as to evaluate the status of the sustainable utilization of forest resources by analysing age compositions against corresponding areas. Further, the stand ages of plantations coupled with allometric equations and biomass expansion factors can be utilized to accurately calculate the carbon sequestration capabilities and greenhouse gas fluxes of forests (Ball et al. 2010), which are critical for elucidating the carbon balances of forests across regional scales (Litvak et al. 2003). Additionally, the mosaic of forest stand ages has implications for the climate, water resources, habitat distributions, and ecosystem processes of landscapes (Helmer et al. 2008). Numerous studies (Fang et al. 2001; Pan et al. 2011b; Wang et al. 2007) have indicated that China’s forests have served as carbon sinks in recent decades. In this regard as well, it is clear that our developed high-quality forest stand age maps will be beneficial for improving carbon budget estimations and accurately projecting future trends for China’s forest plantation stands.

5 Conclusion

Based on an assembled LTSS consisting of annual Landsat time series observations and through the use of the VCT algorithm, spatial analysis model, and random forest regression algorithm, we mapped plantation stand ages for a VCT post-disturbance regrowing forest and VCT undisturbed forest in our study area. The results clearly revealed plantation stand age dynamics, both spatially and temporally, as well as the effects of forest management policy changes on forest age dynamics spanning 1987 to 2017. Since deforestation is contingent on the type, age, and management purposes of forest stands, our results demonstrated that mapping the distribution of plantation stand ages, in terms of forest type, can lay the foundations for improved and prudent forest management for large forest areas. Moreover, temporally consistent multi-temporal plantation stand age products may be considered as valuable data sources for multiple applications, including forest trajectory prediction and carbon sequestration dynamics. We believe that the combination of VCT, spatial analysis model, and random forest regression is robust for synergistically inferring historical plantation stand age distributions. This provides a critical new tool for assessing plantation stand management and carbon accounting under various climate change scenarios.