Spatio-Temporal Land-Use/Cover Change Dynamics Using Spatiotemporal Data Fusion Model and Google Earth Engine in Jilin Province, China

Abstract

Jilin Province is located in the northeast of China, and has fragile ecosystems, and a vulnerable environment. Large-scale, long time series, high-precision land-use/cover change (LU/CC) data are important for spatial planning and environmental protection in areas with high surface heterogeneity. In this paper, based on the high temporal and spatial fusion data of Landsat and MODIS and the Google Earth Engine (GEE), long time series LU/CC mapping and spatio-temporal analysis for the period 2000–2023 were realized using the random forest remote sensing image classification method, which integrates remote sensing indices. The prediction results using the OL-STARFM method were very close to the real images and better contained the spatial image information, allowing its application to the subsequent classification. The average overall accuracy and kappa coefficient of the random forest classification products obtained using the fused remote sensing index were 95.11% and 0.9394, respectively. During the study period, the area of cultivated land and unused land decreased as a whole. The area of grassland, forest, and water fluctuated, while building land increased to 13,442.27 km² in 2023. In terms of land transfer, cultivated land was the most important source of transfers, and the total area share decreased from 42.98% to 38.39%. Cultivated land was mainly transferred to grassland, forest land, and building land, with transfer areas of 7682.48 km², 8374.11 km², and 7244.52 km², respectively. Grassland was the largest source of land transfer into cultivated land, and the land transfer among other feature types was relatively small, at less than 3300 km². This study provides data support for the scientific management of land resources in Jilin Province, and the resulting LU/CC dataset is of great significance for regional sustainable development.

1. Introduction

Land-use/cover change involves a process whereby humans analyze the natural characteristics of land, utilizing biological, technological, and legal policies, and subsequently modify socio-economic activities [1,2]. Land-use/cover change can directly affect the diversity of biological species, contribute to climate change and ecosystem transformation [3], and increase human–environment interactions [4]. With the development of science and technology and the rapid growth of the global population, human activities have caused distinct changes in global land cover [5,6]. Conducting land-use analysis is vital for understanding the intricate relationship between humans and nature, offering insights into current land-use patterns, and guiding the implementation of appropriate interventions. The study of human land-use and the resulting land-cover change allows for the assessment of current human land-use policies and the prediction of their consequences, in order to achieve sustainable land use [7,8].

LUCC is a critical subject in the realms of global climate change and sustainable development, emerging as a key factor in human–land relationship studies. In recent years, numerous researchers have applied remote sensing technology to investigate land-use changes. Recently, the development of domestic and foreign earth observation technology has entered a new stage, and a wide range of multi-temporal, multi-dimensional, and multi-scale remote sensing data have become available [9]. The traditional techniques for extracting information on land-use/cover changes are based on remote sensing data provided by a single sensor, such as MODIS, Landsat, etc. [10]. Due to the limitations of the physical properties of the remote sensing satellites themselves, the spatial and temporal resolutions of single remote sensing data cannot both be prioritized. Although high-spatial-resolution remote sensing images have obvious advantages in terms of the richness of the information expression of surface details, their spatio-temporal data are missing to a large extent due to the long revisit cycle of the data itself, the lower temporal resolution, and the greater influence of clouds, rain, and snow, etc. [11]. This makes it difficult to achieve continuous monitoring of the same region at the temporal or spatial scales, and, thus, leads to the inability to meet the requirements of continuous dynamic tracking and monitoring [12]. For high-temporal-resolution remote sensing data, due to its low spatial resolution in areas with fragmented or heterogeneous surface landscapes, one image element range may contain many different land-cover types, thus seriously affecting its remote sensing classification accuracy [13]. Therefore, the challenge of combining the advantages of different remote sensing satellite sensors to obtain remote sensing data with a high temporal and spatial resolution and to ensure the spatial and temporal continuity of 30 m remote sensing images is a significant problem.

Spatio-temporal fusion technology refers to the generation of high-spatio-temporal-resolution data through fusion algorithms based on high temporal–low spatial- and high spatial–low temporal-resolution remote sensing data [14,15]. Without changing the existing observation conditions, the method of fusion of remote sensing data with different spatial and temporal resolutions can improve the spatial and temporal resolution of remote sensing images [16,17,18]. In recent years, research on high-temporal- and high-spatial-resolution data fusion methods has mainly focused on constructing fusion models. The initial fusion models primarily utilized a linear hybrid algorithm, but the heterogeneity and complexity of the actual ground cover resulted in the algorithm no longer being applicable [19,20]. On this basis, ZhuKov et al. [21] and Lorenzo et al. [22] proposed an improved model by quantitatively analyzing the multiple differences that exist between the target and neighboring image elements. Gao et al. [16] proposed the STARFM fusion algorithm, which combines the advantages of Landsat’s high spatial resolution and MODIS’s high temporal resolution to generate high-spatial- and high-temporal-resolution synthetic data, and then verified the correlation between the fused data and Landsat data. The validation results show that the STARFM fusion algorithm can estimate the reflectance data with a high spatial and temporal resolution in the region. Since then, many researchers have improved the STARFM fusion algorithm. Hilker et al. [23] proposed the STAARCH algorithm based on the STARFM algorithm to complete the mapping of forest disturbances. Zhu et al. [18] developed the STARFM model and proposed the ESTARFM model, which is more suitable for monitoring vegetation dynamics in complex mountainous areas. Zhu et al. [24] also proposed the FSDAF model, which is a flexible spatio-temporal fusion model. It can detect rapid surface changes, i.e., it is able to capture reflectance changes caused by the conversion of land-cover types, and, thus, can effectively predict changes in land-cover types. Guo et al. [25] proposed a flexible object-level (OL) spatio-temporal adaptive reflectance fusion model (STARFM) processing strategy, which can better preserve the structural information of the image, and is 102.89 times faster than the original algorithm with high computational efficiency, and analyzed and verified the effectiveness of the OL processing strategy. Relatively few studies have been conducted on the application of multi-source spatial and temporal fusion techniques to long time series land-cover and its related products. For example, Hao et al. [26] used the ESTARFM fusion method to extract the lake water body area information of Siling Co Lake in Tibet during the period 2000–2014, in order to monitor the spatial and temporal changes of the lake with spatio-temporal fusion data. Liu Han and Gong Peng [27] proposed spatial and temporal fusion and reconstruction techniques based on virtual constellations to build a seamless 30 m day-by-day data cube SDC land-cover and land-use dynamic mapping of China in the 21st century.

With the development of geographic information, big data analytics, and other technological tools, machine learning methods are expanding in the field of land-cover change research, showing their ability to integrate and process training data [28]. Algorithms such as support vector machine (SVM), classification and regression tree (CART), and random forest (RF) have been widely used for land-cover classification [29,30,31]. However, there is still room for optimization of machine learning algorithms in terms of their performance in processing different satellite datasets and the improvement in the algorithmic classification accuracy [32,33]. Some commonly used remote sensing indices have good spatial resolution for vegetation, water, buildings, etc. Using these indices can address the limitation of relying solely on original image data, and can improve the classification accuracy of the remotely sensed images and the accuracy of information extraction [34,35].

Ensuring the continuous production of LU/CC products with high spatiotemporal resolution requires continuous monitoring on a large scale and with high accuracy, and traditional geographic information and remote sensing methods are challenged by their high computational complexity and processing costs [36]. The emergence and expansion of the GEE satellite remote sensing cloud storage and cloud computing platform integrates global remote sensing images from multiple sources and scales, and realizes functions such as full-waveband, high-intensity parallel processing, time series analyses, and mapping [37,38], enabling efficient and continuous remote sensing monitoring. Land-use dynamics changes have been extensively studied on the GEE platform, which has a very high computational efficiency [39].

In this study, in the area of Jilin Province, based on Landsat and MODIS fused long time series high-temporal- and high-spatial-resolution satellite images and the GEE platform, the random forest classification algorithm was used to fuse remote sensing indices to map the LU/CC dataset from 2000 to 2023, and to clarify the temporal trend and spatial characteristics of each land-use type in Jilin Province. This study provides technical support for the accurate acquisition of LU/CC information, and the resulting LU/CC land-use change data in Jilin Province are of great significance for regional sustainable development.

2. Study Area and Data Source

2.1. Study Area

Jilin Province is located in northeastern China, between longitude 121°38′–131°19′ E and latitude 40°50′–46°19′ N, in the geometric center of northeast Asia, which is composed of Japan, Russia, South Korea, North Korea, Mongolia, and northeastern China. It is adjacent to North Korea and Russia in the east, Heilongjiang Province in the northeast, and Inner Mongolia Autonomous Region and Liaoning Province in the southwest, with a total border length of 1438.7 km. Jilin Province has significant differences in geomorphology, with the terrain tilting from southeast to northwest, bounded by the Big Black Mountain in the center, mountains and low hills in the east, and plains in the northwest. This gives it the characteristics of the three major plates, and it is rich in arable land resources, with a vast area of forests. It is a key forestry province of the country, and the Changbai Mountainous Area in the northeast is an important ecological barrier in China. Since the 21st century, intense human activities have led to drastic changes in Jilin Province, which have had a large impact on resources and the environment. Jilin Province has eight prefectural-level cities, namely, Changchun, Jilin, Siping, Liaoyuan, Tonghua, Baicheng, Baishan, and Songyuan, and the Yanbian Korean Autonomous Prefecture, with Changchun as the provincial capital, covering a total area of 187.4 thousand square kilometers, as shown in Figure 1.

2.2. Data Sources

The GEE platform is a global-scale geospatial analysis tool that not only integrates massive geographic and remote sensing data resources, but also possesses powerful cloud computing capabilities, which provides an effective way to study geoscience-related issues [40,41].

2.2.1. Landsat Data

The satellite remote sensing data used in this paper were the Landsat Surface Reflectance SR dataset for the years 2000–2023. The 2005 and 2010 images are from Landsat 5 Thematic Mapper TM Imagery, the image from 2000 is from Landsat 7 Enhanced Thematic Mapper Plus (ETM+), and imagery from 2015, 2020, and 2023 is from Landsat 8 Operational Land Imager (OLI), as shown in

Table 1. Landsat satellite remote sensing data used in the study.

Data Set	Satellite	Sensor	Resolution	Year	Selected Bands
LANDSAT/LT05/C01/T1_SR	Landsat 5	TM	30 m	2005, 2010	Blue, green, red, NIR, SWIR1, SWIR2
LANDSAT/LE07/C01/T1_SR	Landsat 7	ETM+	30 m	2000	Blue, green, red, NIR, SWIR1, SWIR2
LANDSAT/LC08/C01/T1_SR	Landsat 8	OLI	30 m	2015, 2020, 2023	Blue, green, red, NIR, SWIR1, SWIR2

. The data products from the TM and ETM+ sensors were atmospherically corrected using the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS). The OLI sensor data products were atmospherically corrected using the Landsat 8 Surface Reflectance Code (LaSRC), which integrates USGS internal algorithms [42,43].

2.2.2. MODIS Data

MODIS, on board the Terra and Aqua satellites, is the most widely used and important optical sensor in the U.S. Earth Observation Satellite System (EOS) series; these were launched by the National Aeronautics and Space Administration (NASA) in December 1999 and May 2000, respectively [44]. As a continuation of NOAA/AVHRR, MODIS data have been greatly improved in terms of performance and quality. Its scanning width is 2330 km, and the maximum spatial resolution can reach 250 m. It can transit scan up to four times a day, so it has a better time resolution. There are 36 spectral bands, realizing the full spectral coverage of visible to red light bands, which greatly improves the spectral resolution, and strengthens the monitoring and identifying of information regarding the change of land surface cover [45]. MODIS data are rapidly updated, easily accessible, and globally free of charge, which greatly improves the data quality and timeliness of remote sensing monitoring [46,47].

The data selected for this study are MOD09A1 products, which are the standard products of the terrestrial L2 level, and the acquired time range is from 2000 to 2023, with a spatial resolution of 500 m and a temporal resolution of 8 days. The data contain seven major bands, the latter six of which were mainly selected to participate in the subsequent spatial–temporal fusion calculation, and data were acquired against Landsat data at intervals of no more than 3 days. Since the image data were from different sensors, it was necessary to adjust the order of the image bands to match (

Table 2. The range of different sensor bands.

Landsat 5 (TM)		Landsat7 (ETM)		Landsat 8 (OLI)		MOD09A1
Band Name	Bandwidth/μm	Band Name	Bandwidth/μm	Band Name	Bandwidth/μm	Band Name	Bandwidth/μm
SR_B1 (blue)	0.45–0.52	SR_B1 (blue)	0.45–0.52	B2 (blue)	0.452–0.512	sur_refl_b03	0.459–0.479
SR_B2 (green)	0.52–0.60	SR_B2 (green)	0.52–0.60	B3 (green)	0.533–0.590	sur_refl_b04	0.545–0.565
SR_B3 (red)	0.63–0.69	SR_B3 (red)	0.63–0.69	B4 (red)	0.636–0.673	sur_refl_b01	0.620–0.670
SR_B4 (nir)	0.77–0.90	SR_B4 (nir)	0.77–0.90	B5 (nir)	0.851–0.879	sur_refl_b02	0.841–0.876
SR_B5 (SWIR1)	1.55–1.75	SR_B5 (SWIR1)	1.55–1.75	B6 (SWIR1)	1.566–1.651	sur_refl_b06	1.628–1.652
SR_B7 (SWIR2)	2.08–2.35	SR_B7 (SWIR2)	2.08–2.35	B7 (SWIR2)	2.107–2.294	sur_refl_b07	2.105–2.155

).

3. Methods

In this paper, oriented to regional characteristics, the OL-STARFM fusion model was selected to generate high temporal and spatial fusion data of Landsat and MODIS based on the platform of the Google Earth Engine (GEE). The long time series LU/CC mapping and spatial and temporal pattern analysis from 2000 to 2023 were realized by using the random forest remote sensing image classification method that fuses remote sensing indices.

3.1. Data Selection and Pre-Processing

After selecting the satellite remote sensing data of the corresponding time phase in the GEE platform, shadow removal and cloud mask operations were carried out using cloud mask [48]. After that, the images of the selected time were synthesized using median filtering, i.e., assigning a median value to each pixel of the entire image stack, so that all remote sensing images of the study year were synthesized into one, and remote sensing images that clearly and completely displayed the surface coverage information of the study area were obtained [49,50].

3.2. OL-STARFM High-Spatial-and-Temporal-Resolution Data Fusion Model

STARFM combines spectrally similar neighborhood pixel variation information for spatio-temporal fusion, which can be expressed as:

$${\widehat{\mathrm{F}}}_{\mathrm{p}}\left(x,y,B\right)={\sum }_{b=1}^M{\sum }_{i=1}^N{W}_b\left(x_i,y_i,B\right)\times \left[F_b\left(x_i,y_i,B\right)+C_P\left(x_i,y_i,B\right)-C_B\left(x_i,y_i,B\right)\right]$$

(1)

The coarse image was formed by resampling to a fine resolution by nearest neighbor interpolation. The principle of the OL processing version of STARFM (OL-STARFM) [25] after segmentation to process the auxiliary fine image can be expressed as:

$$F_p^{’}\left(s,B\right)=F_b\left(s,B\right)+M\left[C_P\left(s,B\right)-C_b\left(s,B\right)\right]$$

(2)

where M[·] refers to taking the median of all pixel values in the object. OL-STARFM uses the median value rather than the expected value to eliminate the impact of poor-quality pixels for the whole segmented object. Especially for fusing high-quality observations, a residual compensation (RC) step is added afterwards to avoid generating uncertainties and to enhance the initial fusion results:

$${\widehat{\mathrm{F}}}_{\mathrm{P}}={\widehat{\mathrm{F}}}_{\mathrm{P}}^{’}+\mathrm{R}$$

(3)

The residual R represents the differences between the temporal changes of coarse images and fine images. The bicubic interpolation was adopted to remove blocky artifacts in the residual. The pixel-level computation performed in (3) does not incorporate adjacent similar information and, therefore, takes less time.

3.3. LU/CC Information Extraction Based on Random Forests with Fused Remote Sensing Indices

According to the actual situation of land use in the study area, land-use types were divided into six categories: cultivated land, grassland, forest, water, building land, and unused land, as shown in

Table 3. Characteristics of land-use category and remote sensing image representation.

Category	Definition	Features of Remote Sensing Images (Standard False Color Composition)
Cultivated land	Farmland that can be cultivated normally in normal years
Grassland	Land dominated by herbaceous vegetation
Forest	Trees, shrubs, and other woodlands
Water	Natural land waters, land for water conservancy facilities and aquaculture ponds
Building	Urban and rural residential areas and industrial and mining transportation land
Unused land	The surface is soil or rock, the land is largely unvegetated

.

3.3.1. Integration of Remote Sensing Indices

In order to improve the accuracy of classification, in addition to the surface reflectance data, this study used various types of remote sensing index information. In this study, in order to better learn the characteristics of each category, the normalized difference vegetation index (NDVI), the normalized difference building index (NDBI), and the modified normalized difference water index (MNDWI) were imported into the fused raw images to enhance the distinguishability of each feature category and improve the classification accuracy. The formula for calculating each index is as follows:

$$\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}=\frac{NIR-RED}{NIR+RED}$$

(4)

$$\mathrm{N}\mathrm{D}\mathrm{B}\mathrm{I}=\frac{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}-\mathrm{N}\mathrm{I}\mathrm{R}}{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}+\mathrm{N}\mathrm{I}\mathrm{R}}$$

(5)

$$\mathrm{M}\mathrm{N}\mathrm{D}\mathrm{W}\mathrm{I}=\frac{\mathrm{G}\mathrm{R}\mathrm{E}\mathrm{E}\mathrm{N}-\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}}{\mathrm{G}\mathrm{R}\mathrm{E}\mathrm{E}\mathrm{N}+\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}}$$

(6)

where NIR, RED, GREEN, and SWIR denote the surface reflectance in the near-infrared, red, green, and short-wave infrared bands of Landsat images, respectively.

3.3.2. Classifier Construction Based on the Random Forest Algorithm

The random forest (RF) algorithm is a nonparametric regression method [51]. In contrast to other regression methods, the RF algorithm does not overfit, does not require variable selection, and includes built-in cross-validation methods that do not require a separate test dataset to evaluate performance [51]. In addition, the RF algorithm makes no distributional assumptions about predictor or response variables and can handle situations where the number of predictor variables greatly exceeds the number of observations. In this study, the number of decision trees was set to 100, the highest accuracy decision tree is selected for classification, and default values were selected for the remaining parameters such as the number of split variables, maximum leaf nodes, and randomized seeds. For sample points, 80% was used for classifier training and 20% for accuracy verification.

3.3.3. Evaluation of the Accuracy of Classification Results

Accuracy evaluation is an essential step in remote sensing information extraction and target recognition, which not only evaluates the accuracy of the results, but also serves as a reference for evaluating the performance of the method and optimizing the parameters [52,53]. Overall accuracy (OA) refers to the distribution of the number of correctly classified image elements along the diagonal of the confusion matrix, which directly reflects the proportion of correctly classified elements, and the kappa coefficient is used to test whether the model prediction results and the actual classification results are consistent. The formulae for the OA and the kappa coefficient are shown below:

$$\mathrm{O}\mathrm{A}=\left(\frac{{\mathrm{P}}_{\mathrm{C}}}{{\mathrm{P}}_{\mathrm{N}}}\right)·100$$

(7)

$$\mathrm{K}\mathrm{a}\mathrm{p}\mathrm{p}\mathrm{a}=\frac{\mathrm{N}{\sum }_{\mathrm{i}=1}^{\mathrm{r}}{x}_{ii}-{\sum }_{i=1}^{\mathrm{r}}\left(x_{i+}\times x_{+i}\right)}{{\mathrm{N}}^2-{\sum }_{i=1}^{\mathrm{r}}\left(x_{i+}\times x_{+i}\right)}$$

(8)

where P_c is the number of correctly categorized pixels, P_N is the total number of pixels. r is the number of rows and columns in the error matrix, x_ii is the number of observations in row i and column i, x_i+ is the marginal total number of rows i, x_+i is the marginal total number of columns i, and N is the total number of observations.

4. Results

4.1. Accuracy Test of the OL-STARFM Spatio-Temporal Fusion Data

In order to know whether the fusion-generated high-temporal-and-spatial-resolution surface reflectance data can accurately and reliably reflect the real surface cover changes and participate in the subsequent land-cover classification, we need to verify the accuracy of the prediction results. In this study, the Landsat_SR data image taken on 7 October 2020 and the MODIS_SR surface reflectance data images on taken 7 October 2020 and 23 October 2020 were selected as input for into the OL-STARFM model, and the predicted results were verified using the Landsat_SR data image taken on 23 October 2020.

In order to carry out the visual inspection accurately, we selected areas of the same size at typical land types such as cultivated land, forest land, water area, and urban area for detailed observation (

Table 4. Comparison of different ground conditions between real and fusion images.

Category	Real Image (True Color Composition)	Fusion Image (True Color Composition)
Cultivated land
Grassland
Forest
Water
Building
Unused land

). The visual inspection showed that the prediction results are closer to the real image, and also highlighted the spatial detail distribution characterized by the 30 m spatial resolution of the Landsat data.

The results from visual observation alone did not accurately represent the precision of the fusion results, so we needed to quantitatively analyze the predicted results. The correlation coefficients between the predicted data bands and Landsat_SR data bands were calculated by randomly selecting 1000 points in the study area with different corresponding bands, and the correlation coefficients between the different bands were positive, with Pearson correlation coefficients of more than 0.85, representing a high degree of correlation. After that, the Pearson correlation test was also performed on the calculated and fused remote sensing indices, and all three indices showed a positive correlation. The Pearson correlation coefficient of the NDVI was 0.8959, the Pearson correlation coefficient of the MNDWI was 0.8918, and the Pearson correlation coefficient of the NDBI was 0.6857, which were all relatively good correlations. The above results show that the fusion results better contain the high temporal–low spatial- and high spatial–low temporal-resolution image information, and could be applied to the subsequent land-cover classification study.

Figure 2 and Figure 3 show the fitting results between different bands and remote sensing indices in the real and predicted images.

4.2. Results of Remote Sensing Extraction of Land-Use/Cover Change Information in Jilin Province

In this study, we developed a random forest remote sensing image classification method based on the fusion of Landsat and Modis data with long time series high-spatial- and high-temporal-resolution and fusion of remote sensing indices, and extracted LU/CC information with high spatial and temporal resolution of multiple time-phases in the study area, with the support of the GEE platform, and finally synthesized the annual average data, as shown in Figure 4.

The images, after fusing the remote sensing indices, show enhanced distinguishability of each feature class and improved classification accuracy.

Table 5. Comparison of accuracy before and after fusion of the remote sensing indices; 1 denotes cultivated land, 2 denotes grassland., 3 denotes forest, 4 denotes water, 5 denotes building land, and 6 denotes unused land, using the 23 October 2020 fusion data as an example.

Before Fusion Remote Sensing Index								After Fusion Remote Sensing Index
	1	2	3	4	5	6			1	2	3	4	5	6
1	17	0	0	0	1	0		1	20	2	0	0	0	0
2	4	10	0	0	0	0		2	2	13	0	0	0	0
3	0	1	22	0	1	0		3	0	0	19	0	0	0
4	0	0	0	14	0	0		4	0	0	0	20	0	0
5	0	0	2	0	17	2		5	0	1	1	0	17	1
6	0	0	0	0	0	17		6	0	0	0	0	0	20
Overall Accuracy = 89.81%								Overall Accuracy = 93.97%
Kappa = 0.8767								Kappa = 0.9274

demonstrates the significant improvement in the accuracy of the classification results of the OA and kappa coefficient evaluation of the confusion matrix before and after fusion of the remote sensing indices, when the training samples were consistent with all other parameters.

The OA and kappa coefficients were calculated based on the confusion matrix to evaluate the accuracy of the classification results, and the confusion matrix was normalized. The results are shown in

Table 6. Classification accuracy.

Year	2000	2005	2010	2015	2020	2023	Mean
Overall accuracy/%	93.09	94.55	92.96	96.48	97.24	96.31	95.11
Kappa	0.9148	0.9323	0.9134	0.9563	0.9658	0.9537	0.9394

and Figure 5. As can be seen from the table, the OA and kappa coefficients from 2000 to 2023 are greater than 92.96% and 0.9134, respectively, and the average values are 95.11% and 0.9394, respectively, indicating that the classification results are more accurate and could provide more accurate data support for the subsequent analysis of the evolution of spatio-temporal patterns.

4.3. Evolution of the Spatial and Temporal Patterns of Land-Use/Cover Change in Jilin Province

4.3.1. Changes in Time Trend

In this study, we counted indicators such as the area, the amount of area change, and the rate of area change of cultivated land, grassland, forest, water, building land, and unused land, and the results are shown in

Table 7. Land-use type area, area change, and rate of area change by period.

Year	Index	Cultivated Land	Grassland	Forest	Water	Building	Unused Land
2000	Area (km2)	80,538.54	17,443.88	76,389.85	2718.03	4718	5591.7
2005		91,559.24	4176.38	77,037.34	3009.07	6947.44	4670.53
2010		78,762.75	16,415.19	78,576.78	2526.83	8811.64	2306.81
2015		72,362.71	22,538.86	76,247.99	3346.58	9253.92	3649.94
2020		67,083.93	22,832.24	80,925.95	3410.7	10,841.82	2305.34
2023		71,951.81	11,773.05	82,513.3	4158.24	13,442.27	3561.3
2000–2005	Area change /km2	11,020.7	−13,267.5	647.49	291.04	2229.45	−921.17
2005–2010		−12,796.5	12,238.8	1539.45	−482.23	1864.19	−2363.73
2010–2015		−6400.04	6123.67	−2328.79	819.74	442.29	1343.13
2015–2020		−5278.78	293.38	4677.95	64.13	1587.9	−1344.6
2020–2023		4867.89	−11,059.2	1587.36	747.54	2600.45	1255.97
2000–2005	Rate of area change/%	13.68	−76.06	0.85	10.71	47.25	−16.47
2005–2010		−13.98	293.05	2	−16.03	26.83	−50.61
2010–2015		−8.13	37.30	−2.96	32.44	5.02	58.22
2015–2020		−7.29	1.30	6.14	1.92	17.16	−36.84
2020–2023		7.26	−48.44	1.96	21.92	23.99	54.48

.

(1)

Cultivated land. The area and the amount of change in the area of cultivated land in different study years were counted, as shown in Figure 6. Through Figure 6 and Table 7, it can be seen that the area of cultivated land increased from 80,538.54 km² in 2000 to 91,559.24 km² in 2005, and then decreased continuously to 67,083.93 km² in 2020. It then increased slightly to 71,951.81 km² in 2023. Overall, cultivated land decreased by an average of 357.78 km² per year, with rates of change of 13.68%, −13.98%, −8.13%, −7.29%, and 7.26% in each phase.

(2)

Grassland. The area and the amount of change in the area of grassland in different study years were counted, as shown in Figure 7. As can be seen from Figure 7 and Table 7, during the study period, contrary to cultivated land, the area of grassland as a whole showed a trend of decreasing, then increasing, and then decreasing, decreasing from 17,443.88 km² in 2000 to 4176.38 km² in 2005, then increasing to 22,832.24 km² in 2020, and decreasing to 11,773.05 km² in 2023. Overall, grassland decreased by an average of 236.28 km² per year, with rates of change of −76.06%, 293.05%, 37.30%, 1.30%, and −48.44% in each phase.

(3)

Forest. The area and the amount of change in the area of forest land in different study years were counted, as shown in Figure 8. As can be seen from Figure 8 and Table 7, the area of forest increased steadily, with an average annual growth of 255.14 km², and rates of change of 0.85%, 2%, −2.96%, 6.14%, and 1.96% in each phase.

(4)

Water. The area and the amount of change in the area of the water body in different study years were counted, as shown in Figure 9. From Figure 9 and Table 7, it can be seen that the overall trend of water area was increasing, with an average annual increase of 60.01 km², but there was a fluctuating trend in the area of the water in individual study years, with a rate of change of 10.71%, −16.03%, 32.44%, 1.92%, and 21.92% for the five phases.

(5)

Building land. The area of construction land and the amount of area change in different study years were counted, as shown in Figure 10. It can be seen from Figure 10 and Table 7 that during the study period, the area of building land continuously increased from 4718 km² in 2000 to 13,442.27 km² in 2023, with an average annual increase of 363.51 km². The area of building land in the five change phases increased by 2229.45 km², 1864.19 km², 442.29 km², 1587.9 km², and 2600.45 km², corresponding to a rate of change of 47.25%, 26.83%, 5.02%, 17.16%, and 23.99%, respectively.

(6)

Unused land. The area of bare land and the amount of change in the area in different study years were counted, as shown in Figure 11. As can be seen from Figure 11 and Table 7, the unused land area showed an overall decreasing trend during the study period, with an average annual decrease of 84.6 km², but there were repeated increases and decreases in 2010–2023, with rates of change of −16.47%, −50.61%, 58.22%, −36.84%, and 54.48% in each phase.

4.3.2. Changes in Spatial Characteristics

The land-use transfer matrix can visualize the type and amount of land use that has changed, which has significance in clarifying the characteristics and laws of transformation between land-use types in each period. This study extracted the land-use change map and land-use transfer matrix for the first and last study years and the interval study years during the period 2000–2023, as shown in Figure 12 and

Table 8. Land-use transfer matrix for different phases.

Type of transfer-out	Type of transfer-in
		Land transfer matrix (2000–2005)							Land transfer matrix (2005–2010)
		Cultivated	Grassland	Forest	Water	Building	Unused		Cultivated	Grassland	Forest	Water	Building	Unused
	Cultivated	68,926.42	3260.86	4182.32	310.6	2624.38	395.86		68,244.91	13,182.89	5392.66	88.72	3725.98	92.32
	Grassland	13,368.64	654.9	2553.62	128.25	665.89	33.66		3445.67	468.72	117.58	5.34	64.67	1.33
	Forest	5582.88	113.83	71,152.05	72.69	523.05	69.79		3026.1	2555.99	72,257.95	100.93	201.58	3.78
	Water	236.08	9.21	179.73	2043.33	143.44	67.28		259.6	12.48	97.12	2207.91	291.96	95.57
	Building	1512.71	39.03	76.44	310.81	2193.76	513.48		2721.17	177.56	583.96	102.91	3152.45	155.21
	Unused	1102.4	25.22	1.99	97.82	742.42	3485.16		1061.3	18.12	130.23	22.76	1373.27	1959.3
		Land transfer matrix (2010–2015)							Land transfer matrix (2015–2020)
	Cultivated	57,420.24	10,969.68	4057.02	535.46	4573.56	1202.78		53,214.63	12,413.53	1877.7	147.49	3715.31	236.69
	Grassland	7826.63	5076.89	3241.43	43.51	211.54	15.75		6828.34	8478.7	6022.84	67.24	1110.28	6.81
	Forest	2753.04	5665.19	69,850.89	190.83	104.77	14.79		1359.45	1596.81	74,129.29	19.28	205.76	0.92
	Water	81.22	53.95	40.83	2152.89	112.55	87.12		284.7	40.19	47.46	2792.33	125.8	12.18
	Building	3365.56	741.96	120.92	317.86	3438.63	824.99		3439.04	268.88	19.94	213.38	4589.41	569.17
	Unused	158.65	6.54	0.43	62.1	658.77	1421.02		1041.09	14.14	2.69	120.96	961.57	1426
		Land transfer matrix (2020–2023)							Land transfer matrix (2000–2023)
	Cultivated	48,756.27	5294.46	4465.19	608.86	6015.31	1027.17		54,730.78	7682.48	8374.11	857.29	7244.52	811.25
	Grassland	14,038	4180.71	2902.48	71.8	1597.71	21.55		9942.33	1979.99	3261.63	359.51	1788.16	73.33
	Forest	5050.48	1161.81	75,706.55	18.62	158.86	3.6		4265.32	1190.8	71,663.49	76.78	303.06	14.86
	Water	22.76	75.56	78.87	2892.59	114.3	176.59		143.65	104.74	95.3	2021.27	141.08	173.01
	Building	3160.69	953.43	435	415.28	4730.36	1013.37		1004.21	447.4	173.7	495.34	2021.19	504.4
	Unused	250.91	22.95	4	80.31	657.77	1235.84		1194.46	283.16	24.72	275.88	1775.43	1901.37

.

As can be seen from Figure 12, the land-use pattern in Jilin Province underwent significant changes between 2000 and 2023. Cultivated land was the most important source of land transfer-out, decreasing from 80,538.54 km² in 2000 to 71,951.81 km² in 2023, with the proportion of the total area decreasing from 42.98% to 38.39%. During the 24-year period, cultivated land was mainly transferred to grassland, forest, and building land, with areas of transferring out of 7682.48 km², 8374.11 km², and 7244.52 km², respectively. Cultivated land was transferred out to other feature types and also received the transfer-in of other features, among which grassland and forest were the most significant, with an area of 9942.33 km² and 4265.32 km². The land transfer among the other feature types was relatively small, at less than 3300 km².

In order to more clearly describe the flow and diversity of LU/CC changes, the land-use transfer was categorized into three types, namely, transfer-in greater than transfer-out (grasslands, water, and building land), transfer-out is greater than transfer-in (cultivated land, and unused land), and transfer-in transfer-out change class (forest). The conversion relationship between the LU/CC types in different periods was quantitatively expressed using chordal diagrams, as shown in Figure 13.

As can be seen from Figure 12 and Figure 13, during the study period, construction land in Jilin Province increased rapidly with the expansion of urban agglomeration and constantly encroached on the surrounding cultivated land, and the land transfer between cultivated land and grassland, cultivated land and forest land, and cultivated land and building land was the most obvious.

(1)

Cultivated land. The total area of cultivated land decreased, mainly in the western part of Jilin Province, where agriculture predominates, and in the areas of Changchun City, Jilin City, and Yanbian Prefecture, which corresponded with the expansion of construction land and tourism development. During the study period, a total of 24,969.65 km² was transferred, of which 33.54% was converted to forest, 30.77% was converted to grassland, and 29.01% was converted to building land. In terms of transfer-in, water, building land, and unused land have always been the stable transfer-in types of cultivated land, and the largest areas transformed into cultivated land were 284.7 km² (2015–2020), 3439.04 km² (2015–2020), and 1102.4 km² (2000–2005). The area of grassland and forest transferred to cultivated land was unstable and fluctuated.

(2)

Grassland. Grassland makes up a relatively small proportion of Jilin Province, and is mainly distributed at in the eastern forest edge area of Jilin City and a small part of Changchun City. With the implementation of the policy of returning cultivated land to forest and grassland and the reconstruction of the ecological environment, the area of grassland increased between 2005 and 2020, but declined in 2020–2023. In the whole study period, under the influence of the policy of returning farmland to forest and grassland, the main transfer of grassland area occurred between cultivated land and grassland. In 24a, both the maximum area transferred in and the maximum area transferred out of grassland were cultivated land, at 7682.48 km² and 9942.33 km². In the interval study, the area of cultivated land transferred to grassland was the largest, with a maximum interval change of 13,182.89 km² (2005–2010), which was strongly influenced by institutional factors. The largest type of grassland transferred was also cultivated land, with an area of 13,368.64 km² in 2000–2005.

(3)

Forest. The forest land is mainly concentrated in the eastern part of the study area and in the green spaces of a few cities in the central and western regions, such as Changchun, Jilin, and Liaoyuan. The spatial pattern is relatively stable, and it is a land type with a low frequency of change. The main types of transfer-in and transfer-out are cultivated land. In terms of transfer-in, the main transfer-in source of forest is cultivated land, and the largest transfer-in area was 5392.66 km² (2005–2010). The largest transfer source was grassland from 2015 to 2020, at 6172.14 km². In terms of transfer-out, forest was mainly diverted to grassland and cultivated land, with maximum diversion areas of 5665.19 km² (2010–2015) and 5582.88 km² (2000–2005), respectively.

(4)

Water. The changes in the water body area in the study area are relatively stable, and mainly distributed in lakes in the western platform plain area of Jilin Province and the five river systems flowing through the province: Songhua River, Liao River, Yalu River, Tumen River, and Suifen River. The main reason for the changes in water was the mutual transformation between water and other land features such as cultivated land. In terms of transfer-in, the main sources were cultivated land and building land, and the largest transfer-in areas were 608.86 km² (2020–2023) and 415.28 km² (2020–2023). In terms of transfer-out, building land and cultivated land were the main transfer-out areas, and the largest transfer-out areas were 291.96 km² (2005–2010) and 284.7 km² (2015–2020), respectively.

(5)

Building land. Building land showed a trend of continuous expansion, and the development tended to start from a city cluster and spread around. The expansion area was mainly concentrated in central and western Jilin Province, with the urban expansion of Changchun City and the western platform plain cities being the most obvious. The area of building land increased year by year. In terms of transfer-in, the main sources of building land transfer-in were cultivated land, grassland, and unused land, and the largest transfer areas were 6015.315 km² (2020–2023), 1597.71 km² (2020–2023), and 1373.27 km² (2005–2010). In terms of transfer-out, the transfer-out area of building land was mainly cultivated land and unused land, and the maximum transfer-out area was 3439.04 km² (2015–2020) and 1013.37 km² (2020–2023).

(6)

Unused land. The unused land area decreased, occasionally fluctuated, and the spatial distribution was relatively small and mainly distributed in the west. Saline–alkali land was more likely to change. In terms of transfer-out, unused land was mainly transferred to building land and cultivated land, and the largest transfer areas were 1373.27 km² (2005–2010) and 1102.4 km² (2000–2005). In terms of transfer-in, the main sources of transfer were also cultivated land and building land, with the largest transfer areas of 1202.78 km² (2010–2015) and 1013.37 km² (2020–2023).

5. Discussion

5.1. Attitudes towards Land-Use Change Dynamics

Table 9. Single-movement attitudes to land-use change at different phases.

	Cultivated	Grassland	Forest	Water	Building	Unused
2000–2005	2.74%	−15.21%	0.17%	2.14%	9.45%	−3.29%
2005–2010	−2.80%	58.61%	0.40%	−3.21%	5.37%	−10.12%
2010–2015	−1.63%	7.46%	−0.59%	6.49%	1.00%	11.64%
2015–2020	−1.46%	0.26%	1.23%	0.38%	3.43%	−7.37%
2020–2023	2.42%	−16.15%	0.65%	7.31%	8.00%	18.16%
2000–2023	−0.46%	−1.41%	0.35%	2.30%	8.04%	−1.58%

shows the one-way dynamic attitudes of different land-use types in Jilin Province in each time period. Grassland had higher values of unidirectional motivation in 2000–2010, especially in 2005–2010, when the unidirectional motivation was as high as 58.61%; the next largest change was in unused land, which began to produce large recurring changes after 2005. The highest single-movement attitudes for building land were 9.45% in the period 2000–2005, followed by 8% in the period 2020–2023. Changes in cultivated land and forests were relatively small in all periods; the single dynamic attitude of water underwent repeated changes. The land-use type with the largest single-movement attitude over the 24a was building land at 8.04%, and the smallest was forest land at just 0.35%.

The calculated integrated land-use dynamics for the five time periods 2000–2005, 2005–2010, 2010–2015, 2015–2020, and 2020–2023 were 1.51%, 1.67%, 0.93%, 0.71%, and 1.97%, in that order (as shown in

Table 10. Integrated dynamic attitudes to land-use change at different phases.

	2000–2005	2005–2010	2010–2015	2015–2020	2020–2023	2000–2023
Integration of dynamic attitudes	1.51%	1.67%	0.93%	0.71%	1.97%	0.38%

). Although the individual land-use types showed different rates of change in different time periods, at the regional level, the combined dynamics of the five time periods showed an increasing, then decreasing, and then increasing trend. During the first 20 years, the integrated motivation was greatest during the period 2000–2005, indicating that land use in the region changed strongly and gradually intensified during the first period, while it gradually decreased during the latter period. Therefore, the rate of overall land-use change in the region has gradually slowed down since then, which suggests that the intensity of the combined effect of natural and anthropogenic factors on the land-use system in the region has been weakening.

5.2. Frequency Analysis of Land-Use Change

By calculating the number of changes in land-use types, land-use data for five time-phases were frequency-calculated and classified into four categories: no-change areas, low-change areas with 1–2 changes, middle-change areas with 3–4 changes, and high-change areas with 5 changes. Further frequency analysis of land-use changes was carried out.

Figure 14 shows that the frequency of land-use change in Jilin Province was dominated by no-change areas, accounting for 50.84% of the province’s land area, mainly distributed in the forests in the eastern part of Jilin Province, and Songyuan, which is dominated by cultivated land. Cultivated land and forests constitute two of the largest land-use types in Jilin Province, and their land-use types are stable. The low-change area account for 29.82% of the province’s area, dominating the total change area, and the main land-use changes were 1–2 type changes. These were mainly distributed in urban areas and the edge of cultivated land and forests, as well as part of the unused land in Baicheng City in western Jilin Province, which is mainly due to the expansion of agriculture and the process of urbanization leading to change to cultivated land and building land, and at the same time the improvement in unused land. The proportion of middle-change areas is relatively small, accounting for 17.54%, and the middle-change areas were mainly located in the northeast of Baicheng, the central part of Changchun, and the northwest of Yanbian. The areas of high change were the smallest, at just 1.8%. There were fewer areas where land-use changes occurred many times in 24a, and all land-use types were basically stable. Constant change leads to unstable land-use types, which may lead to damage to the land. The results show that frequent land type planning is being carried out in areas of high change, and that more attention should be paid to land in areas of high frequency change in order to carry out accurate planning assessments.

The statistic and comparison of different frequency areas across various cities (autonomous prefecture) in Jilin Province were illustrated in Figure 15. The no-change area is the largest in YB, covering 32,382.73 km². The low-change area is most widespread in BC and CC, while the middle-change area is largest in BC. The high-change area is most extensive in Jilin, indicating more frequent land-use planning activities. Additionally, it is observed that in three cities (autonomous prefecture) in the eastern part of Jilin Province (BS, TH, YB), the unchanged area is larger than the changed area. In contrast, the remaining cities have a larger changed area than the unchanged area.

By analyzing the frequency of land-use changes, we can better understand the causes and impacts of those changes and understand their patterns and trends, so as to provide a scientific basis for future land-use planning and policymaking, and predict future land use. The study of land-use change, especially its environmental and social impacts, needs to be strengthened in future research.

6. Conclusions

Based on long time series and high spatio-temporal resolution satellite images and the GEE platform fusion of Landsat and MODIS, this study developed an exponential random forest classification algorithm integrating remote sensing, carried out land-use classification research from 2000 to 2023, and analyzed the time trend and spatial characteristics of each land-use type. It provides basic information for the comprehensive management and sustainable development of Jilin Province. The main conclusions are as follows:

(1)

The prediction results using the OL-STARFM method were close to the real images, and also highlighted more spatially detailed features. The different bands obtained by spatio-temporal data fusion showed a positive correlation, and the Pearson correlation coefficient reached more than 0.85.

(2)

The land-use type maps with 30 m resolution for six time phases in Jilin Province from 2000 to 2023 were drawn by using the classification method of random forest remote sensing images based on the remote sensing index. The average overall accuracy and kappa coefficient of the products were 95.11% and 0.9394, respectively.

(3)

During the study period, the cultivated land area increased from 2000 to in 2005, then decreased continuously until 2020, and then increased slightly to 71,951.81 km² in 2023, with an average annual decrease of 357.78 km².The overall grassland area showed a trend of decreasing, then increasing, then decreasing, for an average annual decrease of 236.28 km². The area of forest increased steadily, with an average annual increase of 255.14 km². The overall water area showed an upward trend, with an average annual increase of 38.96 km², but it showed fluctuations in individual study years. The area of building land continuously increased, with an average annual growth of 363.51 km². The overall unused land area showed a downward trend during the study period, with an average annual decrease of 84.6 km², but there were fluctuating changes in 2010–2023.

(4)

During the study period, the pattern of land use in Jilin Province changed significantly. Cultivated land was the most important source of transfer-out, and the proportion of the total area decreased from 42.98% to 38.39%. Grassland, forest, and building land were mainly transferred out, and the transfer-out areas were 7682.48 km², 8374.11 km², and 7244.52 km², respectively. Grassland was the largest source of transfer-in, at 9942.33 km². The land transfer among other feature types was relatively small, at less than 3300 km² each.