Rapid land cover classification using a 36-year time series of multi-source remote sensing data

: Long-time series land cover classification information is the basis for scientific research on


Introduction
Land cover classification is important to enable detailed studies of temporal and spatial environmental change, land resource management and sustainable development [1,2].Changes in land cover can affect carbon (C) balance; for example, a study in Shandong Province, China showed that between 2010 to 2020 land cover change resulted in the loss of 106 × 10 4 t C stored in vegetation [3].
Classification of land cover is usually based on natural geographic features such as vegetation type, climatic conditions and topographic features that enable the construc-tion of different types of thematic classification e.g., urban land [4], biogeoclimatic ecosystems [5], and forest types [6].Land classification methods traditionally rely on historical data of land classification and field observations than can require a large amount of time and resources to process as image-based land classification was mainly achieved through visual interpretation of photogrammetry.Subsequently, the availability of remotely sensed data enabled land classification based on statistically analysis of spectral features extracted from image pixels [7].As the availability and diversity of multi-source remote sensing data has increased there have been opportunities to greatly improve the accuracy of land classification.
. Remote sensing has an important role in determining land cover types as multi-sensor-derived waveband information can be used to classify land use cover quickly and reproducibly at different temporal and spatial scales [8].For example, Tadese et al. [9] used remote sensing data as a basis for analyzing and understanding the long-term dynamics of land use and land cover change in the Awash River Basin.Remote sensing imagery can also be used to generate macro-time-series land cover datasets for a region, country or even globally.An example is the Global Land Cover 30 series (GlobeLand30) dataset, which consists of ten primary land cover classes i.e., water bodies, wetland, artificial surfaces, cultivated land, permanent snow/ice, forest, shrubland, grassland, bareland and tundra [10].The release of GlobeLand30 provides a database for large-scale land cover change studies and has been used for large regional-scale studies [11].Whilst the above studies demonstrate the value of land classification at the spatial scale, the datasets are only available for specific years and are not regularly updated as the spectral characteristics of land cover or landscape features can vary interannually.As a result, sample points selected for analysis in one year are not optimal for other years, which can create issues related to training datasets and model migratability [12].To resolve this limitation, a sample point migration approach was developed that enables migration of classification thresholds for a feature from a single chronology to a long-time series dataset [13].
Google Earth Engine (GEE) has been recognized as a powerful tool for processing large-scale Earth observation data, with the ability to access and process large amounts of multi-source, multi-scale and time series remote sensing data via a cloud platform [14].GEE provides access to a variety of datasets in an integrated system, including various satellite image sources, geophysical data, climate data, and demographic data that facilitates the use of time series and multi-source datasets for land cover mapping [15,16].For example, Sidhu et al. [17] accessed the GEE platform's utility in processing raster and vector image manipulations for spatio-temporal analysis of urban and wetland land cover types in two subregions of Singapore, affirming the spatio-temporal analysis capabilities of GEE.However, most existing studies focus on one land cover type or generate land cover maps for certain areas at specific times of image collection.As a result, these studies often find it difficult to incorporate long-time series datasets.The utility of GEE for land cover detection using annual Landsat derived normalized difference vegetation index time-series data was demonstrated by Huang et al. [18] to create a dynamic map of land cover change in Beijing over a 30-year period with an overall accuracy of 86.61% The mulit-petabyte curated catalogue of geospatial datasets available in GEE permits and improves classification results by reducing the likelihood of dataset gaps and uncertainty through the provision of multiple sources of data [19].Multi-source remote sensing data is particularly effective at improving the efficiency of land cover classification as the data fusion and integration of spectral, spatio-temporal, and thermal information from multiple sensors can improves the accuracy of classification [20].For example, Li et al. [21] generated a land cover map of the entire African continent at a resolution of 10 m using a combination of Sentinel-2, Landsat-8, Nighttime Light and MODIS data.
Machine learning algorithms such as maximum likelihood [22], support vector machines [23], random forest (RF) [24] are recognized as accurate and effective methods of analyzing large-dimensional and complex spatio-temporal data when compared to traditional parametric algorithms [25].Selection of a good classification method is a key factor in the classification process that is dependent on the analysis objectives; for example, RF is one of the most frequently used supervised machine learning methods due to its high efficiency and accuracy in identifying single-class elements such as urban number space [26] in remotely sensed imagery as well as its ability to distinguish between multiple land types [27,28], time series data [29] and complex farming areas [30].The improvement of machine learning methods to achieve efficient, fast, and accurate land classification for long time series remains the focus of research.
In this study, we implemented the RF classifier in GEE to perform time series land classification at different spatial scales with Landsat-8 and Sentinel-2 datasets for the vegetation growing season in 2022.Our overarching aim was to use different land classification models constructed using multi-source remote sensing variables to establish an efficient, accurate, and general land classification model for time series datasets, and to identify land classification sample points and migration thresholds based on the differences in sample point image values without land classification changes.Our objectives were to: (1) determine the threshold value of sample point migration based on no change in land class; (2) analyse the accuracy of the land classification model produced using a 36-year time series of Landsat remote sensing imagery and high-precision Sentinel imagery based on threshold values; (3) determine the optimal RF land classification model based on different combinations of multi-source remote sensing variables and compare the impact of image resolution on the classification accuracy.

Study area
Shanxi Province is located within the Loess Plateau and the Yellow River Basin (N34°34′-40°44′, E110°14′-114°33′) and occupies a total area of 156,700 km 2 .Mountains account for more than 80% of the total surface of the region with the topography highest in the northeast and lowest in the southwest, with an average altitude of 1,500 m.Shanxi Province is an important coal energy base in China, with the retained reserves of coal resources reaching 270.9 billion metric tons.At the same time, Shanxi Province contains seven national nature reserves and is an important ecological barrier between mining activities and the Yellow River Basin.Within the Jinzhong Coal Base of Shanxi Province, the Huodong National Coal Planning Area covers a total area of 4,110 km 2 .The region is widely forested and includes the Taiyue Mountain National Forest Park that is an intimate mix of mining and forestry operations.The study area and land classification sample sites are shown in Figure 1.

Data sources
The Landsat series of satellites collect data at a resolution of 30 m and have been providing fundamental data for long-time series scientific research on a global scale since their launch in 1972.In this study remotely sensed data from 1 st June 2022 to 31 st August 2022 was used to capture the spectral reflectance of vegetation and assist in identification and extraction of information on land cover types, such as forests and grasslands, while effectively distinguishing bare ground and other landscape features.
Sentinel-2 satellite data offers 13 spectral bands, which include four 10 m, six 20 m, and three 60 m spatial resolution bands.MultiSpectral Instrument (MSI), Level-1C data is the standard of the Sentinel-2 archive and represents the top of the atmosphere (TOA) reflectance.Sentinel-2 imagery is commonly used to monitor land use and land cover change on a global scale and is designed to provide high-resolution, multispectral remote sensing data for monitoring surface change and environmental conditions.
In addition to the above images, we used the NASA digital elevation model (DEM) and Sentinel-1 synthetic aperture radar (SAR) as multi-source remote sensing images for land classification.All the multi-source remote sensing images involved in land classification are shown in Table 1.

Image pre-processing
The pre-processing of optical remote sensing images included image stitching, declouding, mosaicking, and cropping.In particular, the image de-clouding methods all remove clouds and cloud shadow elements by calling the QA quality bands of Landsat and Sentinel-2 data and operating the mask bit by bit.The mosaic processes of the images were both fused using the median method, which in turn resulted in the Landsat series of images from 1986-2022 and Sentinel-2 remote sensing images of vegetation growing seasons from 2019-2022, respectively.
The Sentinel-1 polarized data GEE has officially undergone ground range detection (GRD) boundary noise removal, thermal noise removal, radiometric calibration, and radiometric correction processes.In this study, the VV and VH polarization bands in the interferometric wide swath (IW) mode, which is suitable for remote sensing studies of land surfaces, were selected.The DEM data were reprojected and resampled to extract variables such as elevation, slope, and aspect as topographic factors to participate in the construction of the land classification model.

Sample point selection
The land classification of Shanxi Province was divided into six types: forest land, grassland, arable land, bare land, water bodies, and impervious surfaces.Additionally, a mining land type was added to the land classification system to account for the Huodong mining area and to assist in the differentiation of the mining and forest in the Taiyue Mountain National Forest Park complex area.
Fixed sample points for different land classifications were selected by importing the sample points into Google Earth to determine the accuracy of the sample points by comparing high-resolution remote sensing images.A total of 1507 sample points from Landsat imagery and 1235 sample points from Sentinel imagery were selected.70% of the sample points were used as training sample points and 30% as validation sample points in the classification process, the specific land classification sample points are shown in Table 2. Spectral features and indices are common methods used to analyze remotely sensed imagery.Spectral features are calculated from ratios or differences between reflectance or emissivity in different bands of the remotely sensed image.These features and indices can be used to extract feature information, monitor vegetation cover, and monitor water quality, among other things.In this study, the Normalized Difference Vegetation Index (NDVI), Normalized Difference Built-up Index (NDBI), Normalized Difference Water Index (NDWI), and Difference Vegetation Index (DVI) were used to calculate the difference in values between forest, grassland, and cropland, respectively, from year to year, and NDBI and DVI were used to calculate the difference between built-up (working) land and bare land.The spectral characteristics of unchanged land types are counted over a number of years so that a reasonable range of thresholds can be determined.In the GEE, the ee.spectralDistance() function was used for image difference statistics.The main purpose of this function is to compute the per-pixel spectral distance between two images.

Random Forest algorithm
Random Forest was used to train a decision tree with randomly selected samples and features from the data set, with the results of the decision trees are assessed to obtain a combined result.The advantage of using the RF algorithm is that it avoids the problem of overfitting and is reliable for handling data such as missing values and outliers.

Feature Model Construction
Comparison of single and multi-source remote sensing variables was conducted by combining different dimensions of remote sensing variables to investigate their influence on land classification results.Four remote sensing feature variables were selected: spectral band, spectral index, topographic features, and SAR data with the specific variable factors shown in Table 3.In the construction of the multi-source remote sensing variables, five combinations of spectral band, spectral band + spectral index, spectral band + SAR, spectral band + spectral index + SAR, and spectral band + spectral index + terrain features + SAR were used, respectively.Accuracy of the classification results was determined by calculating the Overall Accuracy (OA) as the ratio of correctly classified number of samples to the total number of samples, which is a common measure of classifier performance.The Kappa coefficient is a statistic used to measure the agreement between classifiers or evaluators.It can be used to assess the agreement between two evaluators on a classification task.Kappa coefficient values range from -1 to 1, with higher values indicating better agreement.

Determination of thresholds
A total of 180 sample points without land classification change were selected by comparing remote sensing images of the same period from 1986 to 2022 that used 30 sample points per land cover class and included data from each spectral band (Blue, Green, Red, Swir1, Swir2) and spectral indices (NDVI, NDBI, NDW) for each point year by year to obtain the maximum and minimum value range (Table 4).The results show that Landsat can vary somewhat in the image bands between bands and indices, but the fluctuation range is between 0.01 and 0.25.The variation between land classes indicated that water bodies are the most stable followed by grasslands; the bands associated with forests fluctuated more in the NDVI and NDWI indices.The final upper threshold value for land classification sample points was set at 0.25 for Landsat long-time series land classification.

Land classification of Landsat imagery
Land cover classification using Landsat remote sensing images from 1986-2022 was conducted using the sample point migration threshold of 0.25 and the accuracy assessed using OA and kappa coefficient with the number of migrated sample points were counted (Figure 3).The results show that the classification accuracy of the images was highest in the years closer to the 2022 fixed land classification, while the difference between kappa coefficient and OA became larger as the number of years from the 2022 initial land classification sample points increased.However, the overall land classification accuracy remained high, with the lowest kappa coefficient being 0.60 and the lowest OA being 0.75 in 1999.The number of classification sample points decreases as the number of years from 2022 increases, with the migrated sample point data remaining stable at 900, which accounts for approximately 60% of the original number of sample points.It is noteworthy, that differences between Landsat TM/ETM and OIL sensor technology can explain the lower accuracy of results from the start of the study 1986 until 2012.

Land classification of sentinel-2 images
To verify the generality of this paper among different remote sensing images and its reproducibility under complex terrain conditions, we selected the Huodong national planning mining area in Shanxi Province with complex terrain conditions as the study area and added a mining class to the land classification system for Sentinel-2 high-resolution remote sensing images from 2019 to 2022.The land cover classification accuracies in different threshold ranges (0.1 -0.4) were assessed separately by counting each waveband for different years of the land class (Table 5).The results show that the land classification accuracy is higher when the threshold value of training sample point migration is set in the range of 0.20 -0.30 and the number of sample points for year-by-year land classification after threshold screening is maintained at about 70% of the original number, which can meet the number of sample points required for land classification to a greater extent.At the same time, the kappa coefficients between 2019 -2021 are stable around 0.90, while the OA is also all around 0.91.4), and the model accuracy is improved with an increase of different variables, especially the combination of spectral band + index + SAR + The combination model of feature variables and terrain has the best effect (Table 6 ; Figure 4).In 2019, for example, the kappa coefficient eventually increased from 0.863 for a single spectral band to 0.910 for Spectral band + index + Terrain + SAR, whilst the OA also increased from 0.888 to 0.927 for the sample variable combinations.In addition, compared to the 2022 participation in land classification accuracy, the sample points after threshold screening can be used to eliminate the misclassification of sample points in the selection process, so that the 2019-2021 land classification accuracy is better than the 2022 land classification accuracy.The land classification accuracy of Landsat-8 (Table 7; Figure 5) with various combinations of variables is lower than those of multi-source remote sensing land classification based on Sentinel-2 imagery.In 2022, for example, the highest land classification accuracy is achieved with the combination of spectral band + index + SAR, and the model combination of spectral band + SAR is better than spectral band + index.2019 and 2020 have the best accuracy for the full variable combination, while the best variable combination for 2021 and 2022 is spectral band + index + SAR.   7, in Huodong mining area forest and grassland area as a whole accounted for about 80% of the whole area, of which forest land accounted for about 40% of the whole study area, the coal mine area accounted for about 1.5% of the whole study area, and the water and bare ground area is only between 0.2%-0.3%.Three high-resolution land classification products from 2022 were obtained for the comparison of forested land classification results at resolutions ranging from 10 to 30 m (Table 8).However, notably, in the list of classification products shown in Table 8 the JAXA/ALOS/PALSAR/YEARLY/FNF4 products does contain both forest and non-forest land classes.Despite this subtle difference in classification procedure and imagery resolution, the area of forested land ranged from 1136.74 to 1418.27 km 2 with both the highest and lowest area estimates being produced at 10 m resolution.

Discussion
In this study the utility of the GEE cloud computing platform for building land cover classification models using multiple sources of Landsat and Sentinel remote sensing imagery at different spatial resolution over a 36-year time series was assessed.High accuracy spatiotemporal land cover classification maps can help to reveal the impact of human activities such as coal mining and urban expansion on land use over time that could enhance our understanding of the impact of population growth, changes in demography and provide an evidence-base to facilitate future government policy decisions; for example accurate assessments to the spatiotemporal changes in forest C stocks in the context of C accounting and net zero targets [31].
Cloud computing platforms such as GEE, PIE-Engine and AI Earth have improved access to high-performance computing necessary to process large and complex datasets and facilitated an increase in both the speed and accuracy of land cover classification.The approach used in this study was to use the GEE platform to conduct classification based on sample point migration and determine the sample point threshold value required to detect land cover classification change.The selection of sample points migration method has the advantage of not needing to choose new sample points for each time period image and thereby the efficiency of the classification process is improved [32].
Fusion of multi-source remote sensing data into composite data products has been shown to improved accuracy of land cover classification [33].In this study, when assessing the classification of both Landsat and Sentinel multispectral images differences in the classification for crops and grassland were apparent because the imagery obtained for the vegetative growing season did not have substantial differences in the image spectra between grassland and crops.This finding supports the requirement for multi-sources of remotely sensed images obtained with different sensors (e.g., SAR and multispectral data available in Landsat and Sentinel series of images) to accurately classify land cover.
In our comparison of publicly available land classification products (Table 8), the land classification results for forested land ranged between 1136.74 and 1418.27km 2 for the Google and ESA products respectively, which is broadly consistent with our own classification of forest land area.The higher estimation of forested land area by the ESA product is likely due to the inclusion of sparse forest land in the classification of forested land, whereas the variance between the other four products is only 50.88 km 2 despite differences in image resolution.
Land cover classification based on the non-parametric RF algorithm is able to handle multi-dimensional and non-linear data sources whilst also removing the requirement for a balanced number of individual sample points [34] unlike the non-parametric minimum distance, maximum likelihood and Bayesian classification methods.The combination of multi-source remote sensing data and the RF method has been shown to perform land cover classification effectively and accurately.Random forest methods of land cover classification have generated higher accuracy outputs compared to other non-parametric machine learning methods such as support vector machine and artificial neural network [35,36].
Generation of accurate land cover maps over a 36-year period has several challenges relating to the detection of land cover change and technological advances.Sensor technology is continually evolving which has improved the diversity, quality, and quantity of remote sensing datasets available for analysis.The difference in satellite sensors between Landsat-8 Operational Land Imager (OLI) and Sentinel-2 MultiSpectral Instrument (MSI) did not have a large impact on the land cover classification results of the same area despite the higher resolution of the Sentinel-2 acquired datasets, which, in theory, should facilitate more accurate land cover classification and reduce the misclassification of features and reduce the necessity to filter imagery [37,38].However, the fusion of multi-source remote sensing datasets that incorporate textural features [39] has resulted in greater improvements in classifications than relying on the increased image resolution.For example, fusion of datasets from different sensors has been shown to improve accuracy of land classification [40], forest biomass estimation [41] and natural disaster monitoring [42].The complex topography and forest species composition and density in the typical mountainous mining area used in the study demonstrated that effective integration of topographic features such as elevation and slope can be more conducive to distinguishing forests from buildings and crops.

Conclusions
The GEE remote sensing cloud platform was used for rapid land cover classification using Landsat 5, 7, 8, and Sentinel-2 remotely sensed images with a time series spanning 36 years.Single sample point migration was used to produce a time series land cover classification map at both the provincial-regional scale and the scale of mining operations.The final sample point migration threshold value that corresponded to no change in classification was 0.25.The optimal combination of multi-source remote sensing variables used to parameterize the RF machine learning algorithm was the spectral band + index + terrain + SAR for both Landsat 8 and Sentinel-2 generated data.The RF model produced a classification map with highest accuracy for the year 2022 using the Landsat 8 data with an OA of 0.90 and Kappa coefficient of 0.919.Our analysis suggests that a higher accuracy can be achieved when imagery with higher spatial and temporal resolution is used.Further work, assessing the collation of low-resolution remotely sensed imagery and machine learning techniques will enable the assessment of a global-scale land cover classification map over a long time series.As sensor technology develops, we expect the accuracy of land cover classification will continue to improve enabling the future identification of land cover classes not yet considered.
To aid visualisation and interpretation, a GEE-based Land Classification based on Spectral Differences (1984 -present) application was developed and is available at the following URL: https://bqt2000204051.users.earthengine.app/view/land-classification-oflandsat-imagery.The main purpose of this land classification program is to allow users to input a predetermined set of land classification points for a specific year, choose a designated threshold, and utilize the RF algorithm to classify land images from 1984 to the current year.

Data Availability Statement:
We encourage all authors of articles published in MDPI journals to share their research data.In this section, please provide details regarding where data supporting reported results can be found, including links to publicly archived datasets analyzed or generated during the study.Where no new data were created, or where data is unavailable due to privacy or ethical restrictions, a statement is still required.Suggested Data Availability Statements are available in section "MDPI Research Data Policies" at https://www.mdpi.com/ethics.

Figure 1 .
Figure 1.Overview of the study area.(a) Landsat 8 RGB image of Shanxi province in 2022, the red outline is the huodong mining area.(b) Sentinel-2 RGB image of Huodong mining area in 2022.

Figure 2 .
Figure 2. Flowchart of land classification method based on machine learning methods and multisource remote sensing variables.It consists of four parts: Data pre-processing, Sample migration, Land classification, and Accuracy assessment.

Figure 3 .
Figure 3. 1986-2022 Landsat land classification and sample sites.The y-axis on the left of the figure represents the accuracy of Kappa coefficient and Overall accuracy, and the y-axis on the right represents the number of sample points of land classification.

Figure 5 .
Figure 5. Landsat-8 Land Classification in 2019.(a) spectral band (b) Spectral Band + Index (3) Spectral band + SAR (d) Spectral band + Index + SAR (e) Spectral band + index + Terrain + SAR.3.4.3.Comparative analysis of the accuracy of land classification products Landsat 8 and Sentinel-2 remote sensing images of land classification results in 2020 are shown in Table7, in Huodong mining area forest and grassland area as a whole accounted for about 80% of the whole area, of which forest land accounted for about 40% of the whole study area, the coal mine area accounted for about 1.5% of the whole study area, and the water and bare ground area is only between 0.2%-0.3%.

Author
Contributions: Conceptualization, X.Y., J.L. and A.R.S.; methodology, X.Y.; software, X.Y., D.Y., and Y.T.; validation, D.Y. and A.R.S.; formal analysis, Y.S.; investigation, X.Y., Y.T. and Y.S.; data curation, X.Y.; writing-original draft preparation, X.Y.; writing-review and editing, X.Y., A.R.S.; visualization, X.Y.; supervision, A.R.S., J.L. and D.Y.; project administration, X.Y., J.L. and A.R.S.All authors have read and agreed to the published version of the manuscript.Funding: This research was support by the National Key Research and Development Program of China (Intergovernmental and international cooperation in science, technology and innovation) under Grant Number 2022YFE0127700; Royal Society International Exchanges 2022 Cost Share (NSFC) under Grant number IEC\NSFC\223567; University-Industry Collaborative Education Program under Grant Number 220702313180236; Alibaba Damo Academy 2023 AI Earth Joint Innovative Research.

Table 1 .
Multi-source remote sensing image data at two different resolutions used in this analysis. .

Table 2 .
Number of sample points for each land classification.

Table 3 .
Multi-source remote sensing variables feature variables.

Table 4 .
Threshold information of each band for sample points without land class change.

Table 5 .
Land classification accuracies for different thresholds in 2019-2021

Table 6 .
Land classification accuracy of sentinel-2 multi-source remote sensing variables in 2019-2022

Table 6 .
Landsat 8 multi-source remote sensing variable land classification accuracy for the years 2019-2022

Table 7 .
Land classification results of different remote sensing images in 2020

Table 8 .
Comparative analysis of forested land classification products in 2022