Long-Term Land Cover Dynamics (1986–2016) of Northeast China Derived from a Multi-Temporal Landsat Archive

Abstract

Northeast China is a major grain production area, an ecological important forest area, and the largest old industrial base which is now suffering from economic growth slowdown and brain drain. Accurate and long-term dynamic land cover maps are highly demanded for many regional applications. In this study, we developed a set of continuous annual land cover mapping product at 30 m resolution using multi-temporal Landsat images. The maps in year 2000 and 2015 were tested using another independent validation dataset and the overall accuracies were 80.69% and 88.38%, respectively. The accuracies of the maps were improved by the integration of multi-temporal Landsat images and post-classification strategies. We found a general trend that the total area of land that experienced a change in land cover each year increased over time. The area change of each land cover type is also detected. The area of forests was $3.92\times {10}^5$ km${}^2$ in 1986, fluctuated under fire disturbance, but declined in a quite high rate over the period of 1989 to 2006, and finally stayed relatively stable in area around $3.58\times {10}^5$ km${}^2$. The expansion of croplands was the leading land cover change from 1986 to 2000, and then the total area of croplands slightly declined under the Grain to Green Project of China, while shrublands, grasslands and wetlands began to increase. The area of impervious surfaces increased by more than 502% during the last three decades, and about 73% of the new built-up area was converted from croplands. We also demonstrated the our maps could capture the important land cover conversion processes, such as urbanization, forest logging activities, and agricultural expansion.

1. Introduction

Northeast China is a large geographical region, spanning a range of 17${}^{\circ }$ in latitude from 118${}^{\circ }$E to 135${}^{\circ }$E and 16${}^{\circ }$ in longitude from 38${}^{\circ }$N to 54${}^{\circ }$N. It consists of three provinces, including Liaoning, Jilin and Heilongjiang. Broadly speaking, Northeast China also encompasses the eastern part of Inner Mongolia, but, in this study, we used the former definition. It is a major grain producing area, including corns, soybeans and japonica rice, as it has fertile black soil and good condition for agricultural machinery, which thus plays a significant role in ensuring food security. Mapping the extent of croplands is important for spatially explicit crop modelling and yield estimation [1,2,3]. Occupying more than 30% of China’s total forest area [4], Northeast China is an ecologically important forest region. The natural vegetation of this area can be divided into four types: (1) cool temperate deciduous coniferous forest, (2) temperate mixed evergreen coniferous-deciduous broad-leaved forest, (3) warm temperate deciduous broad-leaved forest, and (4) temperate steppe [5]. According to a previous estimation, forests in Northeast China store 1.0–1.5 Pg C and contribute to 24% to 31% of forest C storage of the whole country [6]. Another study found that forests in this region tended to be a carbon source after a long time of overharvesting and degradation [7]. Forests in Northeast China are also significant habitats for many endemic wildlife, such as Siberian Tigers. The expansion of human settlements and croplands cuts off some of the wildlife corridors and breaks their habitats into fragments [8,9]. Thus, dynamic monitoring of forests provides essential input for carbon cycle modelling and background information for biodiversity conservation. Northeast China used to be the greatest industrial base before 1978, but it became a relatively backward region with a slowdown in economy growth [10]. The urban growth rate of this region is the lowest compared to other regions reflected by nighttime light data [11]. Long-term dynamics of urban development are also important for policy makers.

With better availability of higher spatial resolution satellite images, several global and regional land cover products are available at 30 m spatial resolution (e.g., [12,13,14,15]), but none of them are updated frequently to monitor land cover changes over the long-term. Some researchers produced map series to monitor the forest change at 30 m spatial resolution. Global forest extent change since 2000 was captured using Landsat TM/ETM+ images [16,17]. Similarly, Landsat TM/ETM+ images were also used to quantify the annual forest cover dynamics of Eastern Europe from 1985 to 2012 [18]. However, continuous monitoring of multiple land cover classes is uncommon. Most of the work produced land cover maps in sparser time series but not at annual frequency. For example, Ruelland et al. [19] used 20 Landsat images covering three periods of time around 1975, 1985 and 2000 to monitor the land cover change in the upper Niger watershed. There are a few studies produced annual land cover map series at local scale usually with mono-temporal image per year. Lyons et al. [20] mapped land cover of a coastal environment of South East Queensland for seagrass monitoring purpose using 27 winter or dry season images (1 per year) selected from Landsat archive from 1972 to 2010.

Several prior studies have produced land cover maps of different time periods for Northeast China. For example, two land cover maps of 1990 and 2000 were produced and compared to detect land cover change of 10 years to analyze the relationship between land cover changes and climate warming [21]. Gao et al. [22] developed three land cover maps using Landsat images acquired in the same season of 1990, 1995 and 2000 to detect land cover changes related to some land use policies. Chen et al. [23] produced four urban expansion maps by classifying 422 Landsat images of 1990, 2000, 2010 and 2015. However, none of these studies are continuous annual land cover mapping over a long time for this region.

The most challenging part of long-term land cover monitoring is the consistency among maps to avoid unreliable conclusions in land cover change caused by map inconsistency. Building a single model to be applied consistently over time was a proven approach in monitoring land cover change to improve the comparability of the map series generated using Landsat 5 images [24]. Since researchers have found the inconsistent spectral reflectance between sensors in Landsat series [25], the feasibility of land cover classification across space, time and sensors using a common model was not tested.

The objectives of this study are: (1) to develop a set of continuous annual land cover map at 30 m resolution for Northeast China; (2) to prove the accuracy improvements by the integration of multi-temporal Landsat images and post-classification strategies; and (3) to illustrate the potential of the map series in capturing the land cover conversion processes.

2. Data

We used Landsat surface reflectance images from 1986 to 2016 as our primary data source, including Landsat 5 Thematic Mapper (TM), Landsat 7 Enhanced Thematic Mapper Plus (ETM+) with Scan Line Corrector available (SLC-on) and Landsat 8 Operational Land Imager (OLI). In order to produce the land cover map series with multi-seasonal imagery, we defined four calendar-based seasons for each year, i.e., spring (March to May), summer (June to August), autumn (September to November), and winter (December to February). Preference for scene selection was based on (1) season of acquisition, (2) year of acquisition, and (3) cloud coverage. If there was no scene found for a specific season in a target year, we used the image with lowest cloud coverage acquired in the same season of recent year instead. The tolerance of time difference in acquisition year was denoted as $\Delta$Y. For example, for season i of year j, we searched for all scenes acquired in season i between year j −$\Delta$Y and year j + $\Delta$Y, sorted them by cloud coverage, and selected the first one in the list with the best image quality. We tested different $\Delta$Y ranging from 0 to 5, and decided to use 2 in order to have fewer substitutes and better image quality. In total, we collected 4834 Landsat scenes, of which 4832 scenes were ordered from United States Geological Survey (USGS), and 2 were collected from the Satellite Ground Station of China. The distribution of image acquisition time is shown in Figure 1. The images from USGS Landsat archive were all processed to Level 1T precision terrain correction, and scenes with cloud coverage higher than 90% were excluded. Since the two images from the Satellite Ground Station of China were without geometric correction, we did geometric correction for them using our own ground control points and the root mean square errors were under 60 m (2 pixels).

Data from 30 m Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM) were used as topographic data in this study. Slope and aspect were calculated from DEM.

Since the brightest areas at night are the most urbanized areas, the distribution and size of cities could be captured by the night light imagery. The Defense Meteorological Satellite Program (DMSP) Operational Linescan System (OLS) Nighttime Lights time series data [26] was used to reduce the errors in identifying impervious surfaces. To avoid the decline of spatial resolution in the mapping result caused by 1 km nighttime light data compared to 30 m Landsat data, we only use the nighttime light data as a mask to identify the potential areas of impervious surface.

3. Methods

3.1. Training Sample Collection

We sampled 2875 locations distributed across Northeast China for training. The land cover types of year 2000 and 2015 at each location were labelled separately by visually interpreting Landsat imagery (Landsat 5/7 images acquired in around year 2000, and Landsat 8 images acquired in around year 2015) and very high resolution images on Google Earth. Historical images on Google Earth were also used as reference. During image interpretation, we labelled not only the land cover type at each sample location, but also the number of pixels that had the same land cover type around the sample location (1 × 1, 3 × 3, 8 × 8 or 17 × 17 pixels). Then, we expanded each of the training sample from the centre to surrounding pixels to enrich the training sample dataset. All spatio-temporally distributed samples were used to train a common model to be applied consistently over time.

The classification scheme had to meet the basic need of mainly modelling communities and natural resource managers, so we included level 1 classes of FROM-GLC classification scheme (see

Table 1. The classification scheme used in this study.

Code	Land Cover Type
10	Croplands (CR)
20	Forests/Shrubs (FS)
30	Grasslands (GR)
50	Wetlands (WE)
60	Water bodies (WB)
80	Impervious surfaces (IS)
90	Barren lands (BL)
100	Snow/Ice (SI)
120	Clouds (CL)

). We merged forests and shrublands in the scheme because the land proportion of shrublands in Northeast China is small and they are found mainly in swamps and mountainous area. We categorized the shrubs in swamps as wetlands, while shrubs in mountainous area are always mixed with forests, so we categorized them as forests/shrublands.

3.2. Classification

Four seasonal land cover maps were produced by classification, and then integrated. For classification, we tested different feature sets and machine learning methods to get optimal performance. Random Forest was selected as the classifier from a trade-off between computation and accuracy. The random forest classifier used 86 features and 100 trees in this study. The feature set contains 3 parts: (1) Landsat spectral bands (band 1—blue, band 2—green, band 3—red, band 4—near infrared (NIR), band 5—short wave infrared (SWIR) 1 and band 7—SWIR 2 of Landsat 5/7; band 2—blue, band 3—green, band 4—red, band 5—near infrared, band 6—SWIR 1 and band 7—SWIR 2 of Landsat 8) of the target image to be classified and the corresponding indices, (2) Landsat spectral bands of the multi-temporal images in the time series (including spring, summer, autumn and winter image of the same year) at the same location and the corresponding indices, (3) topographic factors, including elevation, slope, and aspect. The indices used for classification include (1) Normalized Difference Vegetation Index (NDVI) [27,28]; (2) Enhanced Vegetation Index (EVI) [29]; (3) Atmospherically Resistant Vegetation Index (ARVI) [30]; (4) Soil Adjusted Vegetation Index (SAVI) [31]; (5) Modified Normalized Difference Water Index (MNDWI) [32]; (6) Normalized Difference Moisture Index (NDMI) [33]; and (7) Normalized Burn Ratio (NBR) [34,35]. Although clouds have been identified by CFMask, we still have training samples of clouds in the classification to reduce errors caused by thin cloud contamination. The two cloud layers (CFMask and clouds identified by our classification model) were combined and masked in post-processing, and the gap of missing pixels were filled according to the multi-temporal classification results.

$$\mathrm{NDVI}=\frac{R_{NIR}-R_{red}}{R_{NIR}+R_{red}},$$

(1)

$$\begin{array}{c}\hfill \mathrm{EVI}=G\times \frac{R_{NIR}-R_{red}}{R_{NIR}+C_1\times R_{red}-C_2\times R_{blue}+L};(G=2.5,C_1=6,C_2=7.5,L=1),\end{array}$$

(2)

$$\begin{array}{c}\hfill \mathrm{ARVI}=\frac{R_{NIR}-(R_{red}-\gamma (R_{red}-R_{blue}))}{R_{NIR}+(R_{red}-\gamma (R_{red}-R_{blue}))};(\gamma =1),\end{array}$$

(3)

$$\begin{array}{c}\hfill \mathrm{SAVI}=\frac{R_{NIR}-R_{red}}{R_{NIR}+R_{red}+L}\times (1+L);(L=0.5),\end{array}$$

(4)

$$\mathrm{MNDWI}=\frac{R_{green}-R_{SWIR1}}{R_{green}+R_{SWIR1}},$$

(5)

$$\mathrm{NDMI}=\frac{R_{NIR}-R_{SWIR1}}{R_{NIR}+R_{SWIR1}},$$

(6)

$$\mathrm{NBR}=\frac{R_{NIR}-R_{SWIR2}}{R_{NIR}+R_{SWIR2}}.$$

(7)

To overcome the problem of coarser resolution of nighttime light than other images, we did not use nighttime light as a feature of the classifier. We used nighttime light to identify the potential area of urban using a digital number (DN) threshold of 7, and trained the classification model separately for potential urban area and non-urban area. For non-urban areas, we removed urban class from the classification scheme.

3.3. Post-Classification Processing

3.3.1. Seasonal Integration

For each path/row location, we produced four seasonal land cover maps using images acquired in different seasons. Normally, the seasonal land cover type is changing caused by the natural growing cycle, the cultivated production cycle, water level fluctuation or snowfall. For example, grasslands start to green up in spring and die back to ground in autumn, so the seasonal land cover type of grasslands could be grasslands, barren lands, or snow cover. However, atmospheric conditions, clouds contamination, and similar spectral characteristics of different land cover types may bring uncertainties to seasonal classification results.

We generated a yearly integrated land cover map by stacking four seasonal land cover maps for less errors and less cloud contaminations [15]. First, the land over type with the highest frequency for each pixel was captured as a base map. Then, the base map was updated pixel by pixel according to priority of different land cover types. The priority for different land cover types were categorized into six levels, and land cover type with lower priority was replaced by higher priority types in the classification result. The highest level of priority is vegetation, including croplands, forests, shrublands, grasslands and wetlands. The second level is impervious surfaces. The third level is water bodies. The fourth level is barren lands. The fifth level is snow and ice. The last level is clouds.

3.3.2. Temporal Consistency Correction

All land cover maps in the time series were produced independently, and the mapping error could be reduced by checking whether specific land cover change is reasonable under the temporal context [36,37]. Then, we developed a temporal consistency correction algorithm to distinguish the reasonable land cover change from the “pseudo-change” caused by mapping error. The algorithm was based on three basic assumptions: (1) Sudden land cover change in a continuous sequence is a very small probability event. For example, the land cover of a sample was continuously type A, but there was an outlier of type B in the time series, we consider it to be a “pseudo-change” caused by mapping error. (2) Conversion of cropland to other land cover types (except impervious surfaces) is a very small probability event before 2000. There are two reasons for this assumption. Firstly, the area of cropland was monotonically increasing before the year 2000 in China, according to the observation of different investigators [38,39,40]. After the launch of “Grain for Green” project in 2000 aiming at vegetation restoration, some cropland area started to be converted to forest or grassland. Secondly, Landsat imagery used in this study is more abundant after 2000, which should result in more reliable classification result, so we only do this correction for the results before 2000. (3) Conversion of the impervious surface to other land cover types is a very small probability event [41,42].

Therefore, there are three steps in the process of temporal consistency correction.

Step 1: For each spatial location, we calculated the relative frequency of each land cover type using a sliding window in multi-year time series of classification results. We tried window sizes between one year and seven years, and chose three years as window size to avoid over-smoothing. If the frequency of the most frequently occurring land cover type was greater than 60%, this pixel was considered “stable pixel”, and the most frequently occurring land cover type was assigned to the land cover type of the focal year (the center of the sliding window). If none of the relative frequency of the land cover types was over 60%, the land cover type of the focal year was kept unchanged.

Step 2: Since we had the assumption that cropland development was irreversible before year 2000, we used the cropland distribution of year 2000 as a benchmark, and inferred the presence of cropland of the previous year by the presence of cropland of the focal year.

Step 3: Similarly, according to the assumption that urban development was irreversible, we back-tracked the distribution of impervious surface from the year 2016.

3.4. Accuracy Assessment

In addition, 422 locations were sampled using the equal-area stratified sampling approach, similar to the sampling approach used in Finer Resolution Observation and Monitoring of Global Land Cover (FROM-GLC) project [14,43]. The land cover types of each sample location in year 2000 and 2015 were interpreted separately based on the very high resolution imagery on Google Earth. Unfortunately, we cannot get access to high resolution imagery in the past back to 2000 for all samples, so we abandoned some validation samples interpreted with high uncertainty due to the absence of very high resolution imagery or there is no strong evidence for us to infer the land cover type in the past. If the location sampled is a mixed pixel of several land cover types, the land cover type covering the largest area in the pixel was recorded. In total, we end up with 297 and 391 validation samples for year 2000 and 2015, respectively. Samples of year 2000 were used to validate the land cover map of 2000, and map of 2015 was validated using samples of 2015. Accuracy assessment was based on confusion matrices. Overall accuracy (OA), producer’s accuracy (PA) and user’s accuracy (UA) of both years were calculated.

Two controlled experiments were also conducted to compare (1) the OA of the maps with and without multi-temporal spectral bands in the feature set; and (2) the OA of the maps with and without post-classification processing.

4. Results

4.1. Accuracy of the Maps

The overall accuracy for the integrated land cover map is 80.69% for the year 2000, and 88.38% for the year 2015. The map accuracy of 2015 is higher than 2000, possibly due to the abundance of Landsat images acquired in 2015 and the sensor improvements of the Landsat 8.

From the confusion matrices presented in

Table 2. Confusion matrix for the integrated land cover map of 2000.

		Prediction
		Croplands	Forests/Shrubs	Grasslands	Wetlands	Water Bodies	Impervious	Bare Lands	Producer’s Accuracy (%)
Reference	Croplands	70	3	0	0	0	0	0	95.89
	Forests/Shrubs	7	87	0	0	0	0	0	92.55
	Grasslands	0	0	4	0	0	0	0	100.00
	Wetlands	10	3	1	1	1	0	0	6.25
	Water bodies	3	2	2	0	25	1	0	75.76
	Impervious	12	0	1	0	1	17	2	51.52
	Bare lands	1	0	0	0	0	0	5	83.33
	User’s accuracy (%)	67.96	91.58	50.00	100.00	92.59	94.44	71.43	80.69

and

Table 3. Confusion matrix for the integrated land cover map of 2015.

		Prediction
		Croplands	Forests/Shrubs	Grasslands	Wetlands	Water Bodies	Impervious	Bare Lands	Producer’s Accuracy (%)
Reference	Croplands	100	2	2	0	0	0	0	96.15
	Forests/Shrubs	9	106	0	0	0	0	0	92.17
	Grasslands	0	1	8	0	0	0	0	88.89
	Wetlands	5	0	0	13	4	0	0	59.09
	Water bodies	0	1	0	2	34	1	0	89.47
	Impervious	7	1	0	0	1	24	2	68.57
	Bare lands	0	0	0	0	0	0	4	100.00
	User’s accuracy (%)	82.64	95.50	80.00	86.67	87.18	96.00	66.67	88.38

, croplands, forests/shrubs, water bodies and barren lands were relatively well-classified in both years. The accuracy of wetlands was improved when the availability of multi-temporal imagery increased because wetlands were with high temporal dynamics. Wetlands are easily misclassified to croplands because the spectral characteristics of wetlands are similar to rice fields. The variation of mapping accuracy was also examined spatially, and we found that the misclassified samples were mainly distributed in the temperate forest steppe zone, the transition ecoregion between the temperate grassland and temperate broadleaf and mixed forest.

From the controlled experiments, we found that, without multi-temporal spectral bands in the feature set of the classifier, the average accuracy of the maps of 2000 and 2015 was 74.2%, while the average accuracy was 78.6% with multi-temporal spectral bands in the feature set but without post-classification process. The multi-temporal images and the post-classification processing were proved useful in this study.

4.2. Land Cover Changes of Northeast China

The land cover map series of the study area during the last three decades are given in Figure 2.

The area of different land cover classes for the last three decades (1986 to 2016) derived from our mapping results are presented in Figure 3. The first majority land cover of the study area in 1986 was forests, which were $3.92\times {10}^5$ km${}^2$, 49.7% of the total area. The area of forests was falling since 1986, with a high rate of decline over the period of 1989 to 2006, and slightly increased after 2006. The forest area in this region in 2016 was $3.58\times {10}^5$ km${}^2$, accounting for 45.3% of the total area. There was a sudden decline of forest cover in 1987 caused by the catastrophic wildfire in the northern part of Greater Xing’an Mountains, and approximately $1.1\times {10}^4$ km${}^2$ of forests were burnt. Our mapping results captured this wildfire disturbance to the landscape and the land cover change pattern on post-fire succession.

The largest area change was found in croplands. The expansion of croplands (+$6.30\times {10}^4$ km${}^2$) was the leading land cover change from 1986 to 2000, in a considerably high speed. Most of the new croplands were converted from forests and other natural land cover types like grasslands, water bodies, and wetlands. The total area of croplands declined under the Grain to Green Project of China, while shrublands, grasslands and wetlands began to increase.

The area of other vegetation types, i.e., grasslands and wetlands, has long been declining, especially from 1996 to 2006. Since 2006, the area for these land cover types started to increase, which is probably the result of conservation. The total area of water bodies fluctuates during these years, which was found obviously lower from 2001 to 2009. The area of water bodies and wetlands has similar variation trend generally—for example, they both reach a local peak in 2013. With the process of urbanization, impervious surfaces kept on expanding since 1986. The expansion of impervious surfaces started to accelerate since 1990 and the growing speed is relatively consistent, except for small leaps in 1998, 2001 and 2012.

We did further analysis to identify the main land cover transformation processes throughout the years. Based on our maps, we found a general trend that the total area of land that experienced a change in land cover each year increased over time, indicating an increasing spatial coverage of disturbance. The major land cover change processes (processes with the top 10 largest area of change) captured in this study are listed in Figure 4. The area of croplands converting to natural land cover classes is roughly increasing since 2000, especially from croplands to forests and wetlands. However, there is also a considerable amount of land converted from natural land cover classes to croplands at the same time. A large area of land converted between forests and barren lands only occurred in the period of 1986 to 1988 caused by the wildfire and forests recovery.

4.3. Application Highlights

The map series depict the continuous land cover change between 1986 and 2016, which can support wide application in regional planning, resource management, and academic studies. The potential of application will be presented through monitoring the change of two land cover classes: impervious surfaces and forests.

4.3.1. Impervious Surface Cover Changes

The urban growth of Northeast China in the last three decades is displayed in Figure 5. The total area of impervious surfaces in the study area expanded 7091 km${}^2$, increased by more than 502%. The map series indicated that about 73% of the new built-up area was contributed by cropland loss around the old built-up area from 1986 to 2016, while new cropland swallowed up other natural land cover classes (forests, grasslands, and wetlands). Besides croplands, bare lands, grasslands, water bodies, and forests are also sources to new built-up area, accounting for 12%, 6%, 4% and 4% of the observed new built-up area. A previous study produced the maps of urban built-up areas of China for the years 1990, 2000 and 2010 mainly by human interpretation [44]. The area of impervious surfaces of these three years in our maps is very close to the urban built-up area published in Wang’s paper, especially for Changchun and Shenyang. Since they defined built-up areas within five pixels as part of the urban areas in urban-rural transitional zone, except those that were connected by roads in Wang’s paper, the new district of Harbin, located to the north of Songhua River, including part of Songbei District and Hulan District, was not counted in the built-up area of Harbin. This is the reason why we have more impervious surfaces in Harbin than the built-up areas in the previous study. The urban dynamic maps produced in this study could provide an accurate base map for regional planning. It is reported that large labour outflow from Northeast China and low fertility rate has led to difficulties to recover from economic slowdown of this region. The three provinces in Northeast China are among the slowest GDP growth in the whole country. Better regional planning will promote economy growth and sustainable development by more efficient placement of infrastructure and zoning for reasonable human activities. Our maps could also help the users review the urban growth history in the last three decades and train the model for future prediction.

4.3.2. Forest Cover Changes

Northeast China is among the three most important forest distribution zones in China, which includes the largest zone of natural forests [45]. The forests in Northeast China are mainly distributed in three subregions, i.e., Greater Khingan Mountains, Lesser Khingan Mountains and Changbai Mountains. This forest zone experienced several periods of different forest management policies since the founding of the People’s Republic of China, for example, diameter-limit cutting, clearcutting and artificial regeneration [46], and unrestricted harvesting [47], under the part of ecological conservation [48]. The map series tracked the logging activities of Northeast China. As shown in Figure 6, clear-cut parcels are found increasing and decreasing in the forest farm of Dongfanghong in Jilin Province. Since in 1987 China introduced cutting quotas to limit annual timber production, the clear-cut areas do not appear to have grown much from 1987 to 1995, except for the new parcels in the pink circle. At the same time, many of the clear-cut areas have partly recovered. However, our maps also reveal that forests continued to be excessive harvested during this period, especially in Heilongjiang Province. As a result of more attention on ecological conservation and more strict control on forest harvesting since the end of 1990s, only a few more clear-cut parcels were found from 1995 to 2006, shown in the cyan circle. The logging activities were not as active as before in this area reflected in our maps from 2006 to 2016, and hopefully we will see better forest protection after the commercial logging of natural state-owned forests was banned in 2016 through continuous monitoring.

Another typical land cover convension type related to forest change is the expansion of croplands frontier into forests, which is the primary driver of deforestation. Figure 7 is an example of croplands expansion into forests in Liaoning province. The impervious surfaces (presented in pink in the figure) is Qingyuan County. There are a lot of hills and valleys in this area. We could see from our land cover maps that forests in the valley were replaced by croplands gradually.

5. Discussion and Conclusions

Accurate and long-term land cover dynamic maps are highly demanded for many regional applications of Northeast China. We developed a set of continuous annual land cover mapping product at 30 m resolution using a common classification model. The overall accuracies were 80.69% and 88.38% for year 2000 and 2015, respectively. Without multi-temporal spectral bands in the feature set and without post-classification processing, the average accuracy of the maps of 2000 and 2015 were 74.2% and 78.6%, respectively. Accuracy was improved by integrating multi-temporal Landsat images and adopting post-classification strategies. The total area of land that experienced land cover change in each year was found generally increasing over time. The resulting map series can capture the processes of major land cover conversion types, such as urbanization, forest logging activities, and agricultural expansion. The maps were consistent over time since they have the same classification scheme, image source, classification model, and post-classification strategies. Consistent and comparable land cover map series will support more robust land cover monitoring and analysis.

The validation of this study can be improved. Limited by the reference images, we only collected validation samples for 2000 and 2015. To evaluate the quality of the map series completely, we should validate more years by improving temporal sampling of the validation sample set. Some of the land cover classes were under-sampled compared to others, which resulted in higher uncertainty on the accuracies of the under-sampled classes. Land cover change areas are often biased because of the mapping errors, and there are approaches to calculate the confidence intervals of the area change [49]. The training sample set can also be enriched by signature extension methods, to cover the spectral variation across scenes, sensors, space and time [50].

Feature selection can help simplify the classification model and shorten training and classification time. We compared the map accuracies between the model with all features and the model with a subset of features, and found that even the least important features in the model still contributed to the overall accuracy, so we kept all features in this study. More experiments are needed to quantify the contribution of different features to the classification. In addition, many studies have proved that additional topographic metrics (such as topographic wetness index, upslope accumulated area) could contribute to the accuracy improvement [51,52], especially for identifying wetlands and water bodies. In the next stage, we will try to test and involve more predictor variables in the classification model, such as additional topographic metrics and phenological metrics.

The suitability of a classification-based (discrete) land cover product and index-based (continuous) product was compared in many applications (e.g., [53]). Both have their advantages: classification-based data can tell us clearly about which land cover type a pixel is for each year, while the index-based data can tell us quantitatively how much has changed for each year. Therefore, we plan to generate another layer of index-based continuous data as additional information to this land cover map series. After this improvement, our product will depict the process of forest loss and gain better. We will know how much a clear-cut parcel or a fire scar has recovered. Seasonal dynamics are also important for land cover monitoring [54]. The long-term map series in this study is annual, i.e., one integrated map per year, which can be improved by producing seasonal dynamic maps, i.e., multiple land cover maps of different seasons.