Mapping Deforestation in North Korea Using Phenology-Based Multi-Index and Random Forest

Phenology-based multi-index with the random forest (RF) algorithm can be used to overcome the shortcomings of traditional deforestation mapping that involves pixel-based classification, such as ISODATA or decision trees, and single images. The purpose of this study was to investigate methods to identify specific types of deforestation in North Korea, and to increase the accuracy of classification, using phenological characteristics extracted with multi-index and random forest algorithms. The mapping of deforestation area based on RF was carried out by merging phenology-based multi-indices (i.e., normalized difference vegetation index (NDVI), normalized difference water index (NDWI), and normalized difference soil index (NDSI)) derived from MODIS (Moderate Resolution Imaging Spectroradiometer) products and topographical variables. Our results showed overall classification accuracy of 89.38%, with corresponding kappa coefficients of 0.87. In particular, for forest and farm land categories with similar phenological characteristic (e.g., paddy, plateau vegetation, unstocked forest, hillside field), this approach improved the classification accuracy in comparison with pixel-based methods and other classes. The deforestation types were identified by incorporating point data from high-resolution imagery, outcomes of image classification, and slope data. Our study demonstrated that the proposed methodology could be used for deciding on the restoration priority and monitoring the expansion of deforestation areas.


Introduction
Forest ecosystems provide ecological benefits to both humans and wildlife [1].Growing global demand for food and fiber is accelerating the pressure on forest ecosystems around the world, from agriculture and logging [2].These pressures have strong impacts on forest ecosystems, which include direct damage to remaining trees, as well as changes in plant and animal species diversity through destruction of the soil cover vegetation [1,3].
North Korea (the Democratic People's Republic of Korea, DPRK) is known to have some of the most degraded forests in the world.Between 1990 and 2015, almost 40 percent of the forests in North Korea have either been converted to fields for food crop production, or cut down for fuel wood.This wide-spread deforestation has led to increased damage by natural disasters [4][5][6].
The characteristic forest landscape in North Korea is complex and heterogeneous; the major landscape types in the forest are hillside farm, unstocked forest, natural forest, and plateau vegetation.Deforestation of the forest landscape has been caused mostly by shifting cultivation practices and natural hazards in North Korea.This has affected nutrient loss from the soil in the forests and under other vegetation.To provide systematic, prioritized restoration efforts to address these problems, it is important to have a classification system for typical kinds of deforestation [7].
Remote sensing can be used for the mapping of forest degradation on a dynamic landscape, covering broad areas but retaining detail in the spatial distribution of target parameters.Most previous deforestation mapping studies were focused on the classification of broad land cover types, such as deforested areas or sand.In the classification algorithm, traditional classification methods, such as ISODATA [4,8] and maximum likelihood [9], were used for classification of degraded forests for mapping.However, these are not appropriate for the classification of combined multi-date images because of the heterogeneous spectral signature of land cover categories over large areas [10].Recently, machine learning algorithms, such as neural networks [11], support vector machines (SVM), [12] and random forest (RF), have been used to overcome these problems.Among these algorithms, RF classification accuracy has been shown to be better than with a support vector machine for soil mapping [13] and forest biomass mapping [14].Thus, RF could help overcome the shortcomings of using multi-parameters [12,[15][16][17][18]. RF can be used in mountain areas, where topographic variables tend to be highly collinear, thus RF could improve forest classification [19,20].Particularly, the RF combined with phenological information could improve the classification of farmland and semi-arid vegetation by capturing the specific seasonal patterns of each landscape type [18,21].Synthetically, RF and phenological information can be used to improve the accuracy of the mapping of specific deforestation types using various environmental variables such as the normalized difference vegetation index (NDVI), normalized difference water index (NDWI), normalized difference soil index (NDSI), and digital elevation models (DEM).
The objectives of this study were to determine a new approach to classify deforested areas in North Korea, and improve the accuracy of classification, using a phenology-based vegetation index derived from MODIS products.This index incorporates various environmental factors, such as vegetation, soil, and water, at a regional scale to improve its accuracy.Better classification of deforestation could provide essential information for decisions about forest management priorities and restoration of deforested areas.It would also improve the ability of the remote sensing science community to accurately monitor the condition of North Korean forests.

Study Area
The study site is North Korea as shown in Figure 1.North Korea encompasses approximately 123,138 km 2 of land located in far eastern Asia.The country is surrounded by rivers and ocean that make up its borders with China and Russia in the north, and the 38th parallel forms its southern border with South Korea.North Korea has four distinct seasons: spring, summer, autumn, and winter.
The majority of the land cover in North Korea is mountain forest.There are mountain ranges in the north and east regions of the country that have peaks that plateau above 2000 m.Examples of high plateau areas include the Gaema and Baeckmu Plateau, which are dominated by shrub and grass [22].The other stocked forests are temperate broadleaf, conifer and mixed forest.
In response to extreme pressure to provide adequate food and energy, the size of the forest area in North Korea has shown a clearly decreasing trend starting in the 1990s [23].This region is experiencing one of the highest deforestation rates worldwide.Moreover, there is high probability of continued forest degradation due to shifting, slash and burn fields for crop cultivation, and logging.Annually, about 127 thousand hectares are being deforested, and those in North Korea are already known to be among the most degraded forests in the world [24].The characteristic deforested land in North Korea includes areas in which crops are cultivated after logging off the trees, and to which cultivation is then shifted onto terrace fields with slopes >16 • to produce food and fuel [4][5][6].In the hillside fields, corn is the main crop.It is an annual plant with shallow roots, and is planted more densely on the hillsides than is usual in flat-land fields, which causes soil erosion and nutrient loss.
Unstocked forest (defined by the standard of South Korea, means areas where the crown cover has fallen to <20% because of slash and burn), is covered with grasses and shrubs.In the flatland the main land use types are: paddies for rice production, fields for producing other crops such as potatoes and other vegetables, and built-up areas.This means that each type of vegetation cover presents a different phenological characteristic.More specifically, the farmland on the hillsides and the unstocked forests have different phenology from the woody and plateau vegetation [25].
Remote Sens. 2016, 8, 997 3 of 14 Unstocked forest (defined by the standard of South Korea, means areas where the crown cover has fallen to <20% because of slash and burn), is covered with grasses and shrubs.In the flatland the main land use types are: paddies for rice production, fields for producing other crops such as potatoes and other vegetables, and built-up areas.This means that each type of vegetation cover presents a different phenological characteristic.More specifically, the farmland on the hillsides and the unstocked forests have different phenology from the woody and plateau vegetation [25].

MODIS Data Collection
The method we propose relies upon the use of multi-temporal and multi-spectral images for classifying deforested areas according to the vegetation, soil, and water phenological characteristics identified.For this, we used the MODIS MOD13Q1 product acquired through USGS Earth Explore [26].MOD13Q1 data were provided every 16 days at 250 m resolution, and included the normalized difference vegetation index (NDVI), enhanced vegetation index (EVI), blue, red, and near infrared (NIR), and short-wave infrared (SWIR) reflectance.NDVI, RED, NIR, and SWIR bands were used in this study to extract biomass and chlorophyll content related to biological characteristics [27] because most forest and crop plants have a relatively unique growth cycle [28].Therefore, we retrieved the data from the MODIS tiles H27V05, H27V04, and H28V05 (which cover North Korea), from March to October (eight months) considering the growth cycles of the vegetation found there.The elevation and slope were added to the classification scheme as predictor variables, which were derived from the SRTM (Shuttle Radar Topography Mission) 30 m digital elevation model (DEM).Both elevation and slope were presented using UTM projection with a bilinear interpolation method, and resampled to 231.7 m cell size to match the MODIS datasets.Thus, the input variables for classification are NDVI, NDSI, and NDWI from March-October, and the topographic variables, such as elevation and slope.

Processing of Time Series Data
The time-series raw data from satellite sensors had signal noise [15].We implemented smoothing of the time-series data and extraction of seasonality with TIMESAT v3.2 [29,30] to reduce the influence of signal noise in the raw NDVI, red, NIR, and SWIR data from March-October 2013.This required at least two steps to obtain the phenological dates.First, a spectral profile was built to fit the discrete data into a smoothing curve [31][32][33]; second, the phenological dates were extracted based on the curves, by applying certain thresholds or filtering models [34][35][36][37].The function-fitting

MODIS Data Collection
The method we propose relies upon the use of multi-temporal and multi-spectral images for classifying deforested areas according to the vegetation, soil, and water phenological characteristics identified.For this, we used the MODIS MOD13Q1 product acquired through USGS Earth Explore [26].MOD13Q1 data were provided every 16 days at 250 m resolution, and included the normalized difference vegetation index (NDVI), enhanced vegetation index (EVI), blue, red, and near infrared (NIR), and short-wave infrared (SWIR) reflectance.NDVI, RED, NIR, and SWIR bands were used in this study to extract biomass and chlorophyll content related to biological characteristics [27] because most forest and crop plants have a relatively unique growth cycle [28].Therefore, we retrieved the data from the MODIS tiles H27V05, H27V04, and H28V05 (which cover North Korea), from March to October (eight months) considering the growth cycles of the vegetation found there.The elevation and slope were added to the classification scheme as predictor variables, which were derived from the SRTM (Shuttle Radar Topography Mission) 30 m digital elevation model (DEM).Both elevation and slope were presented using UTM projection with a bilinear interpolation method, and resampled to 231.7 m cell size to match the MODIS datasets.Thus, the input variables for classification are NDVI, NDSI, and NDWI from March-October, and the topographic variables, such as elevation and slope.

Processing of Time Series Data
The time-series raw data from satellite sensors had signal noise [15].We implemented smoothing of the time-series data and extraction of seasonality with TIMESAT v3.2 [29,30] to reduce the influence of signal noise in the raw NDVI, red, NIR, and SWIR data from March-October 2013.This required at least two steps to obtain the phenological dates.First, a spectral profile was built to fit the discrete data into a smoothing curve [31][32][33]; second, the phenological dates were extracted based on the curves, by applying certain thresholds or filtering models [34][35][36][37].The function-fitting parameters used in TIMESAT were for the double logistic function, and the adaption strength was set to 2.0, season cutoff was set to 0, and the start and end of season threshold was set to 0.2.

Training Samples Collection
The representative samples for each class were collected by field survey, reference data from Google Earth and land cover map of North Korea which was provided by the Korean Forest Services in 2007.The direct collection of regions of interest (ROI) in the data-scarce, remote areas of North Korea is limited due to political conditions.In this study, we could only obtain field survey points from the China side of the North Korea-China border (Figure 3A).On the northern border along the Tumen River we were able to collect 112 GPS points and simultaneously observe the land cover types on the North Korean side of the river.We also interpreted the visual criteria of each land cover class (Table 1) in the satellite imagery of Quickbird and Google Earth [38].Since GE launched in 2005, many studies have been conducted using GE to explore the reference data [15,39].The hillside field and unstocked forest still have confusing visual characteristics in satellite images.Therefore, the land cover map of North Korea, classified using spot images, was referenced to distinguish deforestation types.

Training Samples Collection
The representative samples for each class were collected by field survey, reference data from Google Earth and land cover map of North Korea which was provided by the Korean Forest Services in 2007.The direct collection of regions of interest (ROI) in the data-scarce, remote areas of North Korea is limited due to political conditions.In this study, we could only obtain field survey points from the China side of the North Korea-China border (Figure 3A).On the northern border along the Tumen River we were able to collect 112 GPS points and simultaneously observe the land cover types on the North Korean side of the river.We also interpreted the visual criteria of each land cover class (Table 1) in the satellite imagery of Quickbird and Google Earth [38].Since GE launched in 2005, many studies have been conducted using GE to explore the reference data [15,39].The hillside field and unstocked forest still have confusing visual characteristics in satellite images.Therefore, the land cover map of North Korea, classified using spot images, was referenced to distinguish deforestation types.This sampling method was demonstrated to have high accuracy for use with RF [40].All training points were selected across North Korea with point centers at least 500 m apart [15].The land cover types of random points were identified by visual interpretation of each class (Table 1).The sample size per class was set to a minimum of 200 points for classification [41], and the total sample size was 1660 points containing the survey points (Figure 3).This sampling method was demonstrated to have high accuracy for use with RF [40].All training points were selected across North Korea with point centers at least 500 m apart [15].The land cover types of random points were identified by visual interpretation of each class (Table 1).The sample size per class was set to a minimum of 200 points for classification [41], and the total sample size was 1660 points containing the survey points (Figure 3).

Field
Crop field with annual crops on flatland.This can be distinguished from paddy or other landscapes by the characteristics of plow lines, rectilinear shapes.

Hillside field
Crop field with annual crops on the hillside.It can be detected by 3D view in Google Earth that offers elevation information.

Unstocked forest
This class is covered with shrubs or grasses, where crown cover of trees has fallen to <20% because of slash and burn.This class can be confused with hillside field, but it can be detected by the texture and whether there are plow lines or not.

Natural forest
Trees are the major components.

Plateau vegetation
This class is covered with shrubs or grasses.It can be detected by elevation information and the distribution information from advanced research.

Class Description
Built up Urban and buildings.

Waterbody
Lakes, reservoirs and rivers.

Paddy
Flooded field used to cultivate rice.

Field
Crop field with annual crops on flatland.This can be distinguished from paddy or other landscapes by the characteristics of plow lines, rectilinear shapes.

Hillside field
Crop field with annual crops on the hillside.It can be detected by 3D view in Google Earth that offers elevation information.

Unstocked forest
This class is covered with shrubs or grasses, where crown cover of trees has fallen to <20% because of slash and burn.This class can be confused with hillside field, but it can be detected by the texture and whether there are plow lines or not.

Natural forest
Trees are the major components.

Plateau vegetation
This class is covered with shrubs or grasses.It can be detected by elevation information and the distribution information from advanced research.

Test Samples Collection
Ground truth points used for assessing the accuracy of our results were developed using a random sampling method, which helps ensure that test samples are chosen in an unbiased manner and in proportion to the underlying, but as yet unmapped, land cover-land use categories [42].
The random points made in ArcGIS 10.1 were overlaid on the Google Earth images identifying land cover classes by the visual criteria described in Table 1.The random samples of natural forest were overabundant; in contrast, there were few samples for built-up, paddy, plateau vegetation, and water, as these areas covered only small fractions of the landscape.Thus, we implemented a stratified random sampling [15] for these classes.The polygons for built-up, paddy, and water categories were extracted from the land cover map provided by the Korean Forest Service as reference, and the random points for these landscape were made.Since there was no information about the plateau vegetation on the forest map, random points for plateau vegetation were made within the polygons with elevation >1800 m.Finally, 999 sample points, with a minimum of 50 points in each class were selected based on the recommendations of Congalton and Green [42] at intervals over 500 m.

Normalized Indices
Normalized indices show a much higher classification accuracy compared with original un-transformed spectral data.This is because normalization can be used to compensate for changing conditions of illumination, surface slope, and spectral variability due to reduction of the viewing angle [43].Such indices can effectively detect changes from early spring growth to late season maturity and senescence.In this study, we used multi-indices to enhance the classification accuracy with NDVI, NDSI, and NDWI, which indicated the existence of vegetation, soil, and water respectively [44].The equations of the indices are as follows: where red, NIR and SWIR are the surface reflectance for band 1, 2, and 7, respectively, of MODIS.
The NDVI is often used as an indicator to extract the seasonal data from the relationships between the spectral reflectance and vegetation canopy characteristics [45][46][47].In this case, R-band was used to represent spectral absorption by vegetation and the NIR-band to represent the reflectance value (Equation ( 1)).Therefore, NDVI can characterize the properties of vegetation well [48,49].NDVI is commonly used to make spatial and temporal comparisons of terrestrial photosynthetic activity; thus, it is one of the indicators of the state of land degradation and of the speed of increase or decrease of photosynthesis, which provides biophysical information [50-52].
The NDSI is more sensitive to canopy structure [53] because it identifies areas where soil is the dominant background or foreground material.NDSI uses the SWIR and NIR-bands (Equation ( 2)), where SWIR represents the difference in reflectance values between soil areas [54].Thus, NDSI is a good indicator for separating areas of vegetation from areas of soil.
The NDWI, using NIR and SWIR reflectance values [55] (Equation (3)), is proposed for use in remote sensing of the liquid water content of vegetation canopies.It reflects surface canopy moisture and forest humidity [55,56].The NDWI has been used for monitoring water stress [57], and for mapping burnt areas in boreal forest [58].Thus, the changing pattern of the NDWI with the seasons could distinguish areas of natural vegetation, farmland, and soil, in relation to moisture.

Classification Algorithm Using RF and Validation
Random forest was used for classifying the deforested areas in this study.RF is an ensemble of classification trees intended to improve classification accuracy [59].RF randomly selects samples of variables many times to produce a large number of classification trees.Each tree is individually trained on a bootstrapped sample of the training data and contributes a single vote for assignment of the most frequent class to the input data [59].RF can handle thousands of input features, and has been applied recently in remote sensing studies [60][61][62].
Multi-temporal satellite images are effective for reducing the uncertainty associated with land cover change [63].The traditional classification algorithms (such as maximum likelihood) may not be appropriate for the classification of combined multi-date images because of the heterogeneous spectral signature of land-cover categories over large areas.To overcome this problem, roles for the RF classifier have recently emerged in remote sensing, and RF has been applied successfully for classifying complex land-cover status and dynamics.
The RF was built using the software R 3.4.2[64] and the random forest package [65].The number of trees (ntree) was set to 500, which proved to be a sufficient number in previous experiments [62].At each node, a number of variables (mtry) were randomly sampled from a random subset of the features, and 100 runs were performed.To estimate the accuracy of trees, the out-of-bag prediction was used to estimate the accuracy of all trees, which represents an unbiased estimate of map accuracy, as long as the reference data were obtained via probability sampling [18].To verify the classification result, we used the confusion matrix to estimate overall, user's and producer's accuracies, and the kappa coefficient [66,67].

Temporal Indices of Land Cover Classes
Figure 4 shows that the phenological characteristics (via NDVI, NDSI, and NDWI) of the paddies, flatland fields, hillside field, unstocked forest, forest, and plateau vegetation fluctuated during the growing season.The NDVI consistently increased from March to June for forest, unstocked forest, and plateau vegetation during the growing season, and decreased during the period from August to October (end of the growing season).For paddy, field, and hillside field the growing season extended to July.The highest NDVI value of all land cover classes was observed in July and August.All landscapes showed slightly different trends of change in NDVI value.Natural forest had the highest NDVI value during all periods (mainly consisting of trees) while unstocked forest consisting of grasses and shrubs had slightly lower NDVI values, but the curves followed the same general form.In the hillside fields, the crops were planted more densely than in the flat fields [56] with resultant NDVI values higher overall in comparison to flatland fields.Natural forest showed the highest NDVI value, followed by unstocked forest, plateau vegetation, and hillside fields, followed by the flatland fields and paddies.
NDSI decreased from April to June or July, depending on the landscape type, as the NDVI increased.Hillside fields exhibited the highest NDSI value among the forest landscapes when farmers were preparing the fields early in the season creating conditions for maximum soil exposure.In March-May, at the start of the growing season, the area of soil exposure in the four types of forest landscapes varied because of the canopy.In this case, the hillside fields had the highest NDSI values, followed by unstocked forest and forest.The highest NDWI value was observed at the start of the growing season because of snow melt when the canopy water content was at maximum.When looking at the three individual graphs in Figure 4 one notices that the six landscape types followed similar trends by Index.However, the relative ranking of these landscape types varied according to the index being considered.The largest differences between forest landscape types occurred in different seasons based on the characteristics of each class.For example, hillside field and forest landscape types had similar trends of NDVI in the growing season and at the end of season.In these periods hillside fields had lower NDVI values of around 0.2, but in the middle of the season the difference narrowed to about 0.1 between hillside field and forest landscape types.It is important to note that the NDWI of plateau vegetation is higher than other classes because of the snow melt rate on plateaus was initially slower than for other forest landscape types, and this remained true until May.Additionally, the NDWI of paddies reached the summer maximum in June because of irrigation.

Temporal Indices of Land Cover Classes and Relevant Variables in the RF Classification
Using models with 500 trees, the importance of the contribution of each variable to the general classification model was calculated [59].The mean decrease in accuracy displayed in Figure 5 is the contribution of each variable to the classification model generated by considering the seasonal multiindices and topographic variables.Variables that have high mean decrease accuracy values are considered to be more important for overall or class-level classification [15].

Temporal Indices of Land Cover Classes and Relevant Variables in the RF Classification
Using models with 500 trees, the importance of the contribution of each variable to the general classification model was calculated [59].The mean decrease in accuracy displayed in Figure 5 is the contribution of each variable to the classification model generated by considering the seasonal multi-indices and topographic variables.Variables that have high mean decrease accuracy values are considered to be more important for overall or class-level classification [15].According to the mean decrease in accuracy, the relevant variables in the RF classification were elevation, NDSI during the growing season (March-May), NDVI in September, and NDWI in April.In general, elevation contributed significantly to increased classification accuracy for each category (Table 2).Since the plateau vegetation, hillside field, and (flatland) field categories had similar NDVI patterns, but different topography, the priority variable for classifying these vegetation types was elevation.NDSI had a high value for exposed soils, and showed a large difference between vegetation types in March.Since field preparation for crops began in March [68], the fields became bare soil and presented the highest NDSI values while the other land-cover types were covered with vegetation or snow at the start of the season.Therefore, in March, NDSI could effectively distinguish (flatland) fields from other vegetation types.At the end of the growing season (September), differences in NDVI values could be observed between farmland (field and paddy) and forest (forest and unstocked forest).The forests were covered with vegetation while farmland (field and paddy) were ready to harvest.The NDWI value in April contributed to a classification for plateau vegetation.In the plateau area, the rate of snow melt was slower than other landscape types and the NDWI value at the start of According to the mean decrease in accuracy, the relevant variables in the RF classification were elevation, NDSI during the growing season (March-May), NDVI in September, and NDWI in April.In general, elevation contributed significantly to increased classification accuracy for each category (Table 2).Since the plateau vegetation, hillside field, and (flatland) field categories had similar NDVI patterns, but different topography, the priority variable for classifying these vegetation types was elevation.NDSI had a high value for exposed soils, and showed a large difference between vegetation types in March.Since field preparation for crops began in March [68], the fields became bare soil and presented the highest NDSI values while the other land-cover types were covered with vegetation or snow at the start of the season.Therefore, in March, NDSI could effectively distinguish (flatland) fields from other vegetation types.At the end of the growing season (September), differences in NDVI values could be observed between farmland (field and paddy) and forest (forest and unstocked forest).The forests were covered with vegetation while farmland (field and paddy) were ready to harvest.The NDWI value in April contributed to a classification for plateau vegetation.In the plateau area, the rate of snow melt was slower than other landscape types and the NDWI value at the start of the growing season was higher than for the other types.For this reason, this high value clearly distinguished the plateau sites from the other vegetation types that occurred near them.
The other variables in each season also contributed to increased accuracy in the classification of each vegetation type.The contributions of these variables are specific for North Korea.For mapping regions where there is heterogeneous vegetation cover, it is necessary to find the variables most effective for classification.Each vegetation type in each season, was representative of different vegetation, soil, and water conditions [44], so the phenology-based indices were found effective for classifying these vegetation cover types [15,18,69,70].

Classification Result and Accuracy Assessment
Classification results in this study demonstrated that phenology-based indices using RF classification are effective in mapping deforested land and in understanding complex vegetation cover.High classification accuracies were achieved (Table 2), and it was possible to successfully distinguish areas covered by shrubs, crops, and plateau vegetation which, in turn, distinguished unstocked forest, hillside field, and plateaus from forest.The producer's accuracy in mapping the hillside field and unstocked forest was 80.3% and 87%, respectively, which indicates deforested areas.The corresponding user's accuracies were 94% and 91.2%.Consequently, the overall accuracy of the resultant map was 89.38%, with a kappa coefficient 0.87 (Table 2).Since the landscape types within forests have similar variation in spectral signatures [71], hillside field was misclassified as other classes, such as flatland field or unstocked forest.Some areas of water were misclassified as built-up for the reason that the sand on the riverside was recognized as bare soil.This is a limitation of classification using data with 250 m resolution.This is a common problem when classifying heterogeneous landscapes using per-pixel classification techniques [72][73][74].Figure 6A is the final classification map of the vegetation cover in North Korea, and Figure 5B shows the areas of deforestation in North Korea projected in this study.Approximately 4 million ha of forestland appeared deforested representing 34.2% of the total forestland in North Korea.In these areas forest was converted to hillside field (~2.7 million ha) and unstocked forest (~1.3 million ha) as shown in Figure 6B.Similar results were reported by Jeong [75] who also observed that about 30% of the forest was degraded in North Korea.
In North Korea, attempts to separate these cover types using traditional multispectral data or vegetation indices (Landsat, MODIS, or NDVI in one-time images) [9,75] have not been successful.The class accuracies indicated that phenology-based indices that also include vegetation, soil, and water information, along with the use of RF consistently increased the individual classification accuracies.Nevertheless, further research is needed to develop a technique to analyze detailed vegetation cover types at high spatial resolution.

Conclusions
The purpose of this study was to classify types of deforested areas caused by crop cultivation and logging and to increase the accuracy of classification over previous studies.We investigated a way to combine phenology-based multi-index distinctions derived from MODIS products and RF, as well as ways to distinguish complex, heterogeneous land cover in forests, such as hillside fields and unstocked forest, from plateau vegetation and natural forest.The outcomes of this study extend beyond those of most previous studies, which have usually been focused only on dryland forest.Previous work also involved the use of single image classifications based only on spectral data to distinguish the types of deforestation, and consequently had difficulty capturing the heterogeneous spectral signature of land cover categories over large areas.The outcomes of this study can be summarized as follows.First, the seasonal patterns indicated by three indices (NDVI, NDSI, and NDWI) showed differences typical of each type of vegetation cover in forest land, so it was possible to overcome the reflectance value confusion that occurs when using only one image, and to increase the classification accuracy.Second, to classify complex land cover and dynamics, RF proved itself a useful tool for classifying a variety of input features.Finally, the results highlighted the types of deforested land and their distribution in North Korea.That classification result showed overall accuracy of 89.38% when phenology-based multi-indices were combined with random forest.Our method greatly improved the classification accuracies for classification of heterogeneous vegetation cover and presented deforested areas more reliably.Therefore, it should be useful to continue monitoring variation of forest areas during forest restoration efforts in North Korea.The ecological impacts of forest degradation in the study area should also be urgently considered.

Conclusions
The purpose of this study was to classify types of deforested areas caused by crop cultivation and logging and to increase the accuracy of classification over previous studies.We investigated a way to combine phenology-based multi-index distinctions derived from MODIS products and RF, as well as ways to distinguish complex, heterogeneous land cover in forests, such as hillside fields and unstocked forest, from plateau vegetation and natural forest.The outcomes of this study extend beyond those of most previous studies, which have usually been focused only on dryland forest.Previous work also involved the use of single image classifications based only on spectral data to distinguish the types of deforestation, and consequently had difficulty capturing the heterogeneous spectral signature of land cover categories over large areas.The outcomes of this study can be summarized as follows.First, the seasonal patterns indicated by three indices (NDVI, NDSI, and NDWI) showed differences typical of each type of vegetation cover in forest land, so it was possible to overcome the reflectance value confusion that occurs when using only one image, and to increase the classification accuracy.Second, to classify complex land cover and dynamics, RF proved itself a useful tool for classifying a variety of input features.Finally, the results highlighted the types of deforested land and their distribution in North Korea.That classification result showed overall accuracy of 89.38% when phenology-based multi-indices were combined with random forest.Our method greatly improved the classification accuracies for classification of heterogeneous vegetation cover and presented deforested areas more reliably.Therefore, it should be useful to continue monitoring variation of forest areas during forest restoration efforts in North Korea.The ecological impacts of forest degradation in the study area should also be urgently considered.

Figure 2
shows an example fit to NDVI, red, NIR, and SWIR temporally smoothed from a forest sample.The twelve phenological variables extracted for each growing season included: length of the season, base level (average of the left and right minimum values), largest data value for the fitted function during the season, seasonal amplitude, rate of increase at the beginning of the season, rate of decrease at the end of the season, small seasonal integral, large seasonal integral, number of seasons in a calendar year, time of the start of the season, time of the end of the season, and time of the middle of the season.Use of these temporal factors allowed the program to have ample data to fit a full function to the growing seasons of North Korea.Remote Sens. 2016, 8, 997 4 of 14 parameters used in TIMESAT were for the double logistic function, and the adaption strength was set to 2.0, season cutoff was set to 0, and the start and end of season threshold was set to 0.2. Figure 2 shows an example fit to NDVI, red, NIR, and SWIR temporally smoothed from a forest sample.The twelve phenological variables extracted for each growing season included: length of the season, base level (average of the left and right minimum values), largest data value for the fitted function during the season, seasonal amplitude, rate of increase at the beginning of the season, rate of decrease at the end of the season, small seasonal integral, large seasonal integral, number of seasons in a calendar year, time of the start of the season, time of the end of the season, and time of the middle of the season.Use of these temporal factors allowed the program to have ample data to fit a full function to the growing seasons of North Korea.

Figure 2 .
Figure 2. Raw and fitted MODIS NDVI (A); red reflectance (B); NIR reflectance (C); and SWIR reflectance (D).TIMESAT seasonality variables derived from the functions are numbered on the figure for the interval from January 2012 to December 2013.

Figure 2 .
Figure 2. Raw and fitted MODIS NDVI (A); red reflectance (B); NIR reflectance (C); and SWIR reflectance (D).TIMESAT seasonality variables derived from the functions are numbered on the figure for the interval from January 2012 to December 2013.

Figure 3 .
Figure 3. Independent training samples (1660 points) (A) were collected regarding eight land cover classes in North Korea with a field survey; the 999 test points (B) used for the classification accuracy assessment were selected by random sampling.The displayed images are Landsat 8 OLE products.

Figure 3 .
Figure 3. Independent training samples (1660 points) (A) were collected regarding eight land cover classes in North Korea with a field survey; the 999 test points (B) used for the classification accuracy assessment were selected by random sampling.The displayed images are Landsat 8 OLE products.

Figure 4 .
Figure 4. Time series indices of MODIS NDVI, NDSI, and NDWI (derived from RED, NIR, and MIR reflectance) for paddy, field, hillside field, unstocked forest, forest, and plateau vegetation categories.Random samples of 100 points were individually extracted for each type.

Figure 4 .
Figure 4. Time series indices of MODIS NDVI, NDSI, and NDWI (derived from RED, NIR, and MIR reflectance) for paddy, field, hillside field, unstocked forest, forest, and plateau vegetation categories.Random samples of 100 points were individually extracted for each type.

Figure 5 .
Figure 5. Mean decrease accuracy values of overall and for each class from the RF classification using phenology-based multi-index and topographic variables.

Figure 5 .
Figure 5. Mean decrease accuracy values of overall and for each class from the RF classification using phenology-based multi-index and topographic variables.

Figure 6 .
Figure 6.Final forest classification map of North Korea: (A) land cover map in North Korea; and (B) distribution of deforested land in North Korea.

Figure 6 .
Figure 6.Final forest classification map of North Korea: (A) land cover map in North Korea; and (B) distribution of deforested land in North Korea.

Table 1 .
Classification types and visual interpretation used for identifying training samples in satellite images.

Table 1 .
Classification types and visual interpretation used for identifying training samples in satellite images.

Table 2 .
Confusion matrix for land cover classification of North Korea.