AUTOMATED CLASSIFICATION OF NATURAL FORESTS WITH LANDSAT TIME SERIES USING SIMPLIFIED SPECTRAL PATTERNS

Natural forests are a basic component of the earth ecology. It is essential for biodiversity, hydrological cycle regulation and environmental protection. Globally, natural forests are gradually degraded and reduced due to timber logging, conversion to cropland, production forest, commodity trees, and infrastructure development. Decreasing of natural forests results in loss of valuable habitats, land degradation, soil erosion and imbalance of water cycle in regional scale. Thus operational monitoring natural forest cover change, therefore, has been in interest of scientists for long time. Forest cover mapping methods are divided to two groups: field-based survey and remotely sensed image data based techniques. The field-based methods are conventional and they have been used widely in forestry management practice. Satellite-image-based methods were developed since beginning of earth observation. These methods, except visual image interpretation, can be grouped to supervised and unsupervised classification that rely on various algorithm as statistical, clustering or artificial intelligence. However, there is little report about method, which can extract natural forests from generic forest cover. Over the last couple of decades, natural forests have been over-exploited by various reasons. This practice led to urgent need of development of fast, reliable and automated method for mapping natural forests. In this study, a new method for mapping of natural forest by Landsat time series is presented. The new method is fully automated. It uses spectral patterns as principal classifier to recognize land cover classes. The proposed method was applied in study area consisted of Ratanakiri of Cambodia, Attapeu of Laos and Kon Tum of Vietnam. About 2000 Landsat images were used to generate land cover maps of the study area across years from 1989 to 2018.


INTRODUCTION
Natural forest is important for biodiversity and environmental protection as well as for socio-economy and the living conditions of forest-dependent populations (Stibig et al., 2014). Globally, natural forests are gradually degraded and reduced due to timber logging, conversion to cropland, production forest, commodity trees, and infrastructure development. Decreasing of natural forests results in loss of valuable habitats, land degradation, soil erosion and imbalance of water cycle in regional scale (Dennis et al., 2008). Thus operational monitoring natural forest cover change, therefore, has been in interest of scientists for long time. Methods for monitoring forest cover using Landsat time series (LTS) have been extensively studied during the last two decades. In general, application of LTS for forest mapping can be categorized into two types: image classification and trajectorybased analysis (Banskota et al., 2014). In image classificationbased analysis, individual images are separately classified, thereby minimizing the problem of radiometric calibration between dates (Coppin and Bauer, 1996). Trajectory based approaches utilize the temporal patterns of spectral variables to detect disturbance types and magnitude. Unlike classification, it involves analysis of a time series of single spectral variables and can be further divided into 4 categories: (i) threshold-based change detection, (ii) single curve fitting, (iii) hypothesized curve fitting, and (iv) trajectory segmentation (Banskota et al., 2014). Analysis of LTS requires careful image selection and preprocessing. Images of LTS are selected according to image quality and acquisition dates. The rules for selecting the best available pixel observation can be related to nearness to a target day of year and avoidance of atmospheric interference (be it cloud, haze or resultant shadows), or prioritization of a particular sensor . Radiometric correction and atmospheric correction also play important roles in LTS preparation. Radiometric correction includes calibration to radiance or at-satellite reflectance, while atmospheric correction is performed to achieve minimal radiometric differences between images observed in different time.
In this study, a new method for mapping of natural forest by LTS is presented. Landsat image data covering provinces Ratanakiri of Cambodia, Attapeu of Laos and Kon Tum of Vietnam since 1989 to 2018 were used to map natural forest cover. Natural forest is assumed as forest cover, which is not disturbed during the course of 30 years. More than 2000 Landsat scenes of sensors TM, ETM+ and OLI were used for testing the proposed algorithm. Automated classification was conducted with base land cover categories: Cloud, cloud shadow, water, wetland, closed forest, open forest, and area of intense human activities. The study showed that the proposed method for mapping natural forest is stable, independent and uniform for natural forest inventory in large areas. The advantage of the newly proposed method is that it provides not only statistic data, but also spatiotemporal patterns on how the natural forest was reduced across time.

Study Area
The study area is composed of three neighbouring provinces: Ratanakiri of Cambodia, Attapeu of Laos and Kon Tum of Vietnam. These provinces situate in western site of the Truong Son range so they are of very similar climate characteristics ( Fig.  1). In the study area, there are two distinct seasons: dry and rainy. The rainy season starts in May and ends in October. The dry season lasts from November to April annually. Due to influences of monsoon climate there is low opportunity to obtain cloud free image, especially during rainy season. Prevalent vegetation cover of the study area is evergreen. However, there are also deciduous broad-leaved forests that lose their leaves during dry season.  The is not true TOA Reflectance, because it does not contain a correction for the solar elevation angle. This correction factor is left out of the Level 1 scaling at the users' request; some users are content with the scene-center solar elevation angle in the metadata, while others prefer to calculate their own per-pixel solar elevation angle across the entire scene. Once a solar elevation angle is chosen, the conversion to true TOA reflectance is as follows: (3)

Where
Planetary TOA reflectance [unitless] Local sun elevation angle; the scene center sun elevation angle in degrees is provided in the metadata. Local solar zenith angle; = 90° -

Cloud Masking
Cloud masking is done by using Landsat Collection 1 Level-1 quality assessment band (QA). The Landsat Collection 1 Level- The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B3-2020, 2020 XXIV ISPRS Congress (2020 edition) 1 Quality Assessment (QA) 16-bit band allows users to apply per pixel filters to all Landsat Collection 1 Level-1 data products. Each pixel in the QA band contains unsigned integers that represent bit-packed combinations of surface, atmospheric, and sensor conditions that can affect the overall usefulness of a given pixel. Landsat Collection 1 Level-1 QA bands (.TIF) are included in the Landsat Level-1 GeoTIFF Data Product downloaded from EarthExplorer (USGS, 2017). The QA band provides useful information on cloud, cloud confidence and cloud shadow confidence. The confidence levels of cloud and cloud shadow are divided to three levels: low, medium and high equivalent to 0%-33%, 34%-46% and 67%-100% respectively.

Methods
The goal of our proposed procedure for natural forest mapping is to detect contiguous forest, which is undisturbed during period from 1989 to 2018. We consider long-time undisturbed forest as natural one. This concept requires classification of study area to at least two classes: natural forest and non-natural forest areas. The procedure, therefore, is broken down to two major steps: classification of land cover and accumulation of bare land areas across time. We assume that natural forest changes are caused by human activities in form of timber logging or conversion to agricultural or developed land. Selective timber logging converts closed to open forest while clear cutting or conversion of forest land to other land use forms always experiencing a period when forest land changed to bare land. Hence, detection of bare land is one of steps for classification of non-natural forest areas.

Simplified Spectral Patterns
The simplified spectral pattern (SSP) concept has been used for automated classification of land cover and water body extraction with Landsat images (Duong, 2018(Duong, , 2016(Duong, , 2012Duong et al., 2017). An SSP is a transformation of the full spectral pattern into a simplified digital form, which allows direct incorporation of the spectral pattern into the classification. The SSP is constructed by non-repetitive pairwise comparison of reflectance values between two bands (Charalambides, 2002). Given a pixel value vector B6={b1,b2,b3,b4,b5,b6}, where b1, b2,...,b6 denote top of atmosphere (TOA) reflectance of bands 2, 3, 4, 5, 6, and 7 of OLI, and bands 1, 2, 3, 4, 5, and 7 of ETM+ sensor, respectively, the simplified spectral pattern is defined by a new 15 digit vector, shown in equation (1). (1) where mi,j is the result of comparison between the reflectance of bi and bj and has values of 0 (if bj<bi.), 1 (if bj=bi), or 2 (bj>bi). An example of SSP construction is presented in Fig. 2.   Table 1. Examples of SSPs for three land cover classes DGVdark green vegetation, BGV-bright green vegetation and DBL-dry bare land. In our method, every land cover class is recognized primarily by an SSP. However, to identify more specifically a particular class of land cover we need to use some supporting variables that are computed using TOA reflectance. There are seven variables, which help us to better describe a shape of a spectral pattern. List of those variables are given in Table 2 Figure 3. List of variables that support better description of a shape of a spectral pattern. B1 ---B3 denote TOA reflectance for bands blue, green, red, near infrared, short infrared 1and shortwave infrared 2 respectively.
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B3-2020, 2020 XXIV ISPRS Congress (2020 edition) Number of supporting variables needed for identification of a land cover class is not fixed for all classes. Some classes such as clear water does not need any supporting variables. Some other classes with high probability of mixture with other classes such as wetland or dry bare land require always supporting variables. The database of SSPs therefore contains both SSPs and supporting variables for every land cover classes. Example of using a software tool to identify SSP and supporting variables for a pixel in a image is given in figure 4. Figure 4. Example of using a software tool to identify SSP and supporting variables. The Pixel info windows shows SSP code, number of SSP codes available within the image, pixel vector in TOA reflectance and values of seven supporting variables.

Automated Classification of Land Cover
Algorithm for the automated classification of land cover is straightforward. Figure 4 shows major steps for the classification of land cover using the database of SSPs. The algorithm is a simple loop of two computation steps: conversion of pixel value vectors to TOA reflectance and comparison of SSPs and supporting variables of a pixel with SSPs and supporting variables. Computation of TOA reflectance is conducted differently for Landsat sensors using information provided in the metadata files. The complexity of spectral characteristics of surface materials results pixels with SSPs not yet defined in the databases. These pixels are classified to unknown class as usual in practice of image classification. If the number of unknown pixels is too large and those pixels threaten quality of classification result, we need to identify new SSP and to update the database of SSP to reduce unknown pixels in classification result. If there is still small number of unknown pixels after classification, spectral matching could be applied to fill those unknown pixels.

Stacks of Images within one year
Landsat images of tropical region are always disturbed by clouds. It is inevitable facts. Clouds and cloud shadows are required to be reliably identified if we need to achieve good result of classification of land cover with Landsat time series. To reach cloud free land cover map in certain time period (dry or wet season) we need to classify more Landsat images and stack them together. Information on clouds and cloud shadows can be extracted from band quality assessment (BQA) that associates with every Landsat image in Landsat collection 1 (USGS, 2019a(USGS, , 2017. By using the BQA we can identify cloud and cloud shadow pixels in one image and replace them by clear pixels from images of other date in Landsat time series. Figure 5 demonstrates stacking of classified images to create cloud free classified images.

Natural Forests Extraction
We assume that natural forest is a forest that is not substantially changed during its existence except consequence of hazard like landslides or fire. This assumption redirects detection of natural forest to detection of change of forest to bare land during its existence. In the study area, two prevalent forest covers evergreen and dry deciduous forest. Dry deciduous forest has leaf-off season from late December to April of the next year. During this period, dry deciduous forest looks like bare land on Landsat image. To avoid miss classification we do not collect Landsat image of area with distribution of dry deciduous forest. We accumulate classified images of one year to find changed and unchanged forest covers. When we stack multiple-year classified images, again, we can mark changed and unchanged forest covers. The long-time unchanged forest cover is then classified to natural forest cover. Figure 6. Displays flow chart of algorithm for automated classification of land cover and natural forest mapping.

RESULTS AND DISCUSSION
Algorithm described in the methods session was implement using C++ language for application in Windows 32-bit and 64-bit environment. The software running in command mode allows batch processing of large amount of Landsat images.  Figure 7 shows changes of natural forest of Ratanakiri from 1989 to 2018.   It is quite hard to validate reliability of classification for past years. We conducted field survey in December 2018 and March 2020 to collect ground data for classification accuracy analysis. For inaccessible area, high spatial resolution imagery of Google Earth was used for validation. Because spectral pattern was used as main classification measure in our proposed method so validation of classification results is equal to validation of spectral patterns and their simplified forms. Field survey in 2018 and the use of high spatial resolution of Google Earth showed very high accuracy of 2018 land cover maps. We assume that spectral patterns of ground materials observed by Landsat sensors are stable across time so we are of high confidence that classification results of Landsat images in the past years have high reliability.

CONCLUSION
In this paper, we introduced a new algorithm for automated classification of land cover with Landsat time series. Forest cover that remain unchanged during long period is considered as natural forest. The algorithm is composed of three processing steps. In the first step, Landsat image is automatically classified by using database of simplified spectral patterns developed for individual Landsat sensors. In the second step, single date classified land cover maps are stacked to create annual cloud free land cover maps. In the third step, by stacking all annual land cover maps we will extract natural forest and area of intense human activities. For each provinces in the study area, we generated series of land cover maps showing spatio-temporal changes of land cover across 30 years. Classification results The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B3-2020, 2020 XXIV ISPRS Congress (2020 edition)