Improvement of Spatio-Temporal Inconsistency of Time Series Land Cover Products

Abstract

In recent years, time series land cover products have been developed rapidly. However, the traditional classification strategy rarely considers time continuity and spatial consistency, which leads to the existence of unreasonable changes among the multi-period products. In order to solve the existing problems, this paper proposes a matrix decomposition model and an optimized hidden Markov model (HMM) to improve the consistency of the time series land cover maps. It also compares the results with the spatio-temporal window filtering model. The spatial weight information is introduced into the singular value decomposition (SVD) model, and the regression model is constructed by combining the eigenvalues and eigenvectors of the image to predict the unreasonable variable pixels and complete the construction of the matrix decomposition model. To solve the two problems of reliance on expert experience and lack of spatial relationships, this paper optimizes the model and proposes the HMM Land Cover Transition (HMM_LCT) model. The overall accuracy of the matrix decomposition model and the HMM_LCT model is 90.74% and 89.87%, respectively. It is found that the matrix decomposition model has a better effect on consistency adjustment than the HMM_LCT model. The matrix decomposition model can also adjust the land cover trajectory to better express the changing trend of surface objects. After consistent adjustment by the matrix decomposition model, the cumulative proportion of the first 15 types of land cover trajectories reached 99.47%, of which 83.01% were stable land classes that had not changed for three years.

1. Introduction

Land cover is a concept that appears with the development of remote sensing technology, which refers to the complex landform composed of natural elements and artificial objects. It contains different temporal and spatial features and therefore presents different states on various spatial and temporal scales [1,2,3]. Since the 1990s, several global land cover products based on remote sensing images acquired by the same sensors have introduced multi-temporal products with 30 m to 1 km resolution. These products provide more comprehensive surface information and are also suitable for long-term surface monitoring and change analysis [4,5,6,7]. One of the main tasks of time series land cover product production research is to extract land cover change data in the same region in different periods, understand the distribution and change trend of land cover classes, and effectively guide human adjustment activities through these changes. The more rational application of resource flow and social development also plays an important role in the sustainable development of natural resources [8,9,10,11,12].

The spatio-temporal inconsistencies of time series land cover products usually refer to the inconsistencies between land cover maps and the actual situation of the land surface, or the unreasonable changes between land cover maps in different periods [13]. Accurate and consistent land cover products are the key to continuous research related to land cover information [14,15,16]. However, due to the influence of many factors such as remote sensing image data quality, training sample selection, and phenological phenomena, the pixel radiation value of the same object in the remote sensing image at different times is often different [17,18]. In addition, there will be alignment errors between different images or between images and the ground truth, which causes the pixel class in the image to be offset from the surface objects. Therefore, these factors will lead to the wrong classification of the land cover map, resulting in the phenomenon of inconsistent spatio-temporal land cover goals [19,20,21].

In the production process of temporal land cover products, the method of single-phase classification is usually adopted, and only the accuracy and spatial correlation are considered [22,23]. However, there are relatively few studies on the continuity on the time scale, which leads to errors in the trend of surface feature changes and a large number of illogical changes when extracting the surface change in a period of time from land cover products [24,25]. In other words, the class changes in surface objects occurring in the current time period are often inconsistent with the actual situation [26]. Therefore, how to solve these problems and improve the accuracy and application value of temporal land cover products is one of the important directions of current research.

In the field of remote sensing, there are some studies on land cover classification that consider time consistency [27,28]. At present, time series land cover products are beginning to pay attention to time consistency, among which the algorithm update of the Moderate-resolution Imaging Spectroradiometer (MODIS) Collection series is the most typical representative [29,30]. Abercrombie and Friedl proposed a new framework for MODIS Collection 5 products. The HMM is used to help distinguish between real and false land cover changes in a classification time series. Taking the HMM as the post-processing step of supervised classification and using the existing HMM algorithm to solve the optimal label sequence, the stability presented is higher than that of MODIS Collection 5. The number of pixels with annual change decreases, the phenomenon of “lag” in the change in land cover labels is changed, and the number of labels monotonically decreases. Most of the changed pixels will only change once or twice [31]. Climate Change Initiative—Land Cover (CCI-LC) v2.0 also avoids the spurious change detection caused by the annual inconsistency of classification in the post-classification change detection method. Each change pixel must have a change class in the classification time series for two consecutive years or more before it is considered a change pixel [32].

In addition to the production of global temporal land cover products, there are many experts and scholars working on solving the time consistency problem of time series land classification in order to improve its accuracy and reliability. In order to more accurately reflect the dynamic change in land cover, some researchers try to incorporate the transition rule of surface objects over time into the process of time series classification or change detection [33,34]. However, the change law adopted does not objectively reflect the change in surface objects. For example, Cai et al. divided the class transformation of surface objects into two classes: reasonable transformation and unreasonable transformation, but the illogical transition matrix was completely based on expert experience [33]. Abercrombie and Friedl applied the same 0.1 transfer probability to all features [31]. In fact, the change law of ground features should be closely related to the ecological environment and land cover class of the study area [35,36,37]. Therefore, when studying the law of ground feature transfer, we should adopt more objective methods combined with the study area, such as the analysis of measured data and the modeling of environmental variables to determine the transfer probability matrix between different surface objects. In this way, the inconsistent and unreasonable changes in land cover products at different time points can be better solved, and the accuracy and availability of products can be improved.

The time window filtering model is widely used in identifying illogical changes. Sexton et al. developed a time filter that recalculates the probability of each class as the square root of the maximum probability of the current year and the joint probability of the previous year and the following year and assigns pixels to the class with the greatest probability [38]. Li et al. proposed an iterative time filtering method, defined a probability function, used a threshold to screen out unreasonable changes, and used the time consistency probability to eliminate the non-logical transformation from urban to non-urban in the land cover trajectory [39]. Liu et al. proposed a spatio-temporal land cover filter (STLCF), which “learns” the knowledge of non-logical land cover change events from the land cover map through the spatio-temporal transfer probability matrix. Using the maximum probability of land cover calculated by the naive Bayes equation to correct for illogical changes, the results show that the STLCF improves the average overall accuracy of annual change detection by about 6% [40].

The Markov random field model is widely used in pixel classification and post-classification processing in temporal land cover. According to the research of Tso and Mather, they developed a context classification model for remote sensing image classification based on the Maximum a Posteriori Markov Random Field (MAP-MRF) framework. The genetic algorithm is used to estimate the model parameters, and it has been verified that the simulated annealing algorithm has better performance than the conditional iterative model algorithm, and the accuracy of the simulated annealing algorithm is improved by about 2% [41]. Wang et al. applied the MAP-MRF, based on the 2010 Finer Resolution Observation and Monitoring-Global Land Cover (FROM-GLC) project, using a random forest classifier to obtain the original category label and improve the land cover classification for 3 consecutive years. The MRF outputs for 2001 and 2010 were then processed using a rule-based label adjustment method with MODIS/Terra Vegetation Continuous Fields (VCF), slope, and synthetic Enhanced Vegetation Index (EVI) series as auxiliary data. The labeling adjustment process re-labels overclassified forests, bodies of water, and bare land into alternative categories with the greatest probability [24].

The HMM is widely used to analyze time series data, where a series of unknown states is inferred from a series of observed states. The HMM is determined by three elements: initial probability distribution, state transition probability, and observation probability distribution. The HMM mainly solves three basic problems: (1) Probability calculation problem; (2) Learning problem, that is, estimating model parameters with known observation series; (3) Prediction problem, that is, to solve the optimal state sequence given the observation sequence [42]. Gong et al. proposed an algorithm to merge information from spatio-temporal neighboring observations into an HMM to improve the time series land cover map originally generated by a support vector machine (SVM). In order to study the effects of different initial distribution and transition probability matrices on HMM classification, experimental schemes with different input elements were designed, and the algorithm was verified by Landsat satellite images and China Environmental Disaster Reduction satellite (HJ-1) images. In order to make use of spatial information effectively, spatial weights are introduced into the HMM. Experimental results show that spatial weights can eliminate most unreasonable land cover changes that may occur in previous pixel classification [43]. Souverijns et al., based on the Google Earth Engine (GEE) platform, used Landsat data and random forest classification, combined with the HMM algorithm, to achieve time consistency and provided a 30-year land cover map with a resolution of up to 30 m for the Sudan-Sahel region [44].

In general, there are still some unsolved problems in the existing consistency adjustment methods. The time filter model can only adjust the classification results of the intermediate years and cannot adjust the consistency of the first and last years in the time series. Although the HMM has a good adjustment effect, it is more dependent on the previous time and cannot better capture the dynamic characteristics of the data, thus affecting its prediction effect [31,45]. The singular value decomposition (SVD) model is an optimal matrix decomposition in the sense of least squares, which can pack the most information into as few coefficients as possible. Huang et al., based on the SVD model, fused remote sensing images with different spatial and spectral characteristics through sparse matrix factorization. The experimental results demonstrate that the model outperforms other algorithms in generating fused imagery with both the well-preserved spectral properties of MODIS and the spatial properties of ETM+ [46]. Guo et al. proposed a low-rank approximate denoising algorithm based on the SVD model. For low-rank matrices, SVD can provide optimal energy compression in the least squares sense. The experimental results show that the method can effectively reduce noise and compete with the most advanced denoising algorithms in quantitative metrics and subjective visual quality at the time [47]. Wang et al. introduced the SVD model into the principal component analysis method, solved the problem of the large computational load of the principal component analysis, and verified the application of SVD in data analysis and latent semantic indexing [48]. Therefore, in the processing of remote sensing images, the image processing method based on SVD has been widely used and has achieved good results. With the advent of the era of big data, dimensionality reduction technology has become the primary task for high-dimensional large-scale data, and SVD can be applied in the field of data analysis and prediction [46,47,48,49]. Because of these characteristics and advantages of SVD, we proposed a matrix decomposition model based on SVD for consistency adjustment of the land cover classification map. This is a kind of classification post-processing method. For the pixels with unreasonable changes, their categories are predicted by the matrix decomposition method, and the improved land cover classification map is generated.

This paper uses the matrix decomposition model and HMM_LCT model for post-classification processing to adjust the consistency of the land cover maps and makes a comparative analysis of the results. The main contributions of this paper are as follows:

(1)

Two new methods of matrix decomposition model and HMM_LCT model are proposed to adjust the consistency of the time series land cover map.

(2)

The matrix decomposition model is applied to the consistency adjustment of the land classification map based on the transition probability matrix and SVD model.

(3)

To solve the problem of lack of spatial relationship in the HMM, the HMM_LCT model introduces a spatial factor.

(4)

The results of the two models are compared and analyzed. The experimental results show that the matrix decomposition model has a better effect on consistency adjustment than the HMM_LCT model.

2. Method

Figure 1 is the overall flow chart of this study. The shared principle of the matrix decomposition model and HMM_LCT model adopted in this paper is introduced, and at the same time, the spatio-temporal window filtering method proposed by predecessors is used for comparative experiments. The general process of this study is as follows: Firstly, we use the initial classification results to calculate the land cover transfer probability matrix as the parameters of the consistency adjustment model. Secondly, based on the SVD method and spatial weights, a matrix decomposition model is constructed to adjust the image consistency after classification. The model mainly consists of using the land cover transfer probability matrix to determine the threshold value and screen out unreasonable variable pixels. Based on the SVD algorithm, the matrix decomposition model is constructed by combining spatially weighted interpolation weights and regression models to predict the unreasonable variable pixels. Then, based on the HMM previously studied, the HMM_LCT model is proposed. The model contains three elements: initial probability A, observation probability B, and state transition probability matrix $\mathsf{\pi }$. In this paper, the product of observation probability and spatial factor I is regarded as the observation probability of the model, and the land cover transfer probability matrix is calculated as the state transition probability matrix through the initial land cover classification maps, and the HMM_LCT model ($\lambda =(A,B\times I,\pi )$) is formed. The forward–backward algorithm is used to calculate the possibility of pixels belonging to various categories, and the initial classification results are adjusted. At last, the two models are compared and analyzed experimentally. The comparison of experimental results mainly relies on the confusion matrix and the trajectory of land category change, and the comparative analysis of local details is carried out. In addition, robustness experiments are conducted using data from different regions with different resolutions.

2.1. Study Area and Data

This paper selected several data sources and study areas: Huangpi District of Wuhan City as the study area of the HMM_LCT model, Shunyi District of Beijing City as the study area of the matrix decomposition model, and Changping District of Beijing City and Yaohai District of Hefei City, respectively, for robustness testing.

Huangpi District is located in the central part of Hubei Province of China; the total area is 2767 square kilometers, and its geographical coordinates are roughly 114°43′~115°31′ east longitude and 30°39′~31°23′ north latitude. The terrain is a typical hilly landform of middle and low mountains. Huangpi District belongs to the subtropical monsoon climate zone with four distinct seasons. Shunyi District is located in the northeast suburb of Beijing; the total area is 1021 square kilometers, and the geographical coordinates are roughly 116°28′~116°58′ east longitude and 40°00′~40°18′ north latitude. Shunyi District belongs to the temperate continental subhumid monsoon climate zone, with four distinct seasons and abundant rainfall. The land resources in this area are diverse and rich, including cultivated land, forest, grassland, towns and villages, industrial and mining land, and river and water conservancy facilities. These two areas have abundant surface elements, obvious seasonal alternations, and moderate area sizes, so we chose Huangpi District and Shunyi District as the research areas. The specific geographical location maps of the study areas are shown in Figure 2 and Figure 3.

In order to obtain a series of comparable data, it is usually necessary to select the same time point for remote sensing data acquisition and processing when making time series land cover products. In order to obtain temporal land cover products with high spatial resolution, four 10 m bands were used for this study. The original images of Huangpi District are shown in Figure 4. GF-1 is a remote sensing satellite successfully launched by the China National Satellite Application Center. The GF-1 satellite series consists of four satellites; the first satellite of the series was launched in 2013, followed by the launch of the same GF-1 02, 03, and 04 satellites in 2018. The satellites are equipped with two 2 m resolution panchromatic and 8 m resolution multispectral cameras and four 16 m resolution multispectral cameras. In this study, three phases of a multispectral 16 m resolution GF-1 image were selected. The original images are shown in Figure 5.

2.2. Land Cover Transfer Matrix

The probability matrix of land transfer refers to the probability matrix that a certain piece of land changes from one class to another class at a certain time. This matrix can reflect the trend and probability of land cover change, which is of great significance for land cover and planning management.

The land cover transfer matrix is used to predict the trend of land cover change by calculating the transfer probability between different classes of land. The purpose of this study is to adjust the land cover map coherently, and it is necessary to screen out some pixels with unreasonable changes through the transfer probability matrix, so it is crucial to calculate the transfer probability matrix and grasp the changing trend of land cover [50]. In this paper, the initial classification data set is selected to calculate the land cover transfer probability matrix. The transfer probability matrix of the initial land cover classification map of the two study areas is shown in

Table 1. Land cover transfer probability matrix of Huangpi District (%).

Land Cover Category	Impervious Surface	Water	Forest and Grass	Plowland
Impervious Surface	88.83	0.31	0. 80	9.90
Water	4.05	81.40	1.24	13.08
Forest and Grass	4.37	6.19	82.71	6.60
Plowland	8.73	2.24	2.00	87.00

and

Table 2. Land cover transfer probability matrix of Shunyi District (%).

Land Cover Category	Impervious Surface	Water	Forest and Grass	Plowland
Impervious Surface	66.66	3.39	2.72	27.23
Water	7.33	82.77	5.47	4.43
Forest and Grass	3.01	0.61	65.73	30.65
Plowland	16.66	2.87	12.47	68.00

, respectively.

2.3. Matrix Decomposition Model

2.3.1. Spatial Neighborhood Information

The land cover type usually has a certain continuity and correlation in space, and this correlation can be better obtained by considering the neighborhood information of pixels so as to improve the rationality of ground object classification. In the experiment, the spatial neighborhood information is introduced for consistency adjustment, and the constructed spatial weight matrix is used as an input in the matrix decomposition model to represent the spatial relationship between each pixel and its neighborhood pixel. Using the weight information in the spatial weight matrix, the eigenvalues and eigenvectors after matrix decomposition can be modified. By introducing the spatial weight matrix, the matrix decomposition model can better consider the spatial distribution relationship of ground objects and improve the expression ability of the model. One of the common methods to introduce spatial neighborhood information is distance-based weights [51]. We need to determine the spatial neighborhood of each pixel by defining a fixed-size window or a fixed range centered on the pixel. Let W be the spatial weight, d_ij represents the spatial distance between pixel i and pixel j, and w_ij represents the corresponding weight. Then, the element w_ij of the spatial weight matrix can be expressed as Formula (1):

$${\mathrm{W}}_{\mathrm{i}\mathrm{j}}={\displaystyle \frac{1}{{\mathrm{d}}_{\mathrm{i}\mathrm{j}}}}$$

(1)

2.3.2. Matrix Decomposition Model Prediction

The principle of the SVD model method is based on matrix eigenvalue decomposition. By calculating the eigenvalues and eigenvectors between matrix A and its transpose $A^T$, the three matrices U, ∑, and V required by the SVD model can be obtained. The eigenvalue represents the scaling factor of the transformation, and the eigenvector represents the transformation direction. In this experiment, the SVD method is used to decompose the land cover map C′(t) after removing unreasonable variable pixels, and its eigenvalues and eigenvectors are obtained. The specific formula is (2):

$${\mathrm{C}}^{\prime }\left(\mathrm{t}\right)=\mathrm{U}\sum {\mathrm{V}}^{\mathrm{T}}$$

(2)

where U is the left singular value vector matrix; ∑ represents a singular value matrix, which is a diagonal matrix, and the elements on the diagonal are the eigenvalues of the matrix; and V is the right singular value vector matrix.

The eigenvalues of the initial land cover classification results of the year to be predicted were obtained, and the eigenvalues and eigenvectors of other time series reference images were obtained using the SVD model. The relationship model between eigenvalues and time was established using the known eigenvalue data and the linear regression model. By selecting a series of time points and calculating the eigenvalues corresponding to each time point, a linear regression of Equation (3) is constructed:

$$\mathrm{X}\left({\mathrm{t}}_{\mathrm{i}}\right)={\mathsf{\beta }}_0+{\mathsf{\beta }}_1{\mathrm{X}}_1+{\mathsf{\beta }}_2{\mathrm{X}}_2+\cdots +{\mathsf{\beta }}_{\mathrm{n}}{\mathrm{X}}_{\mathrm{n}}+{\mathsf{\epsilon }}_{\mathrm{i}}$$

(3)

where $X\left(t_i\right)$ is the eigenvalue in year t, ${\beta }_n$ is the linear regression coefficient, and ${\epsilon }_i$ is the error. The regression model is used to forecast other years to be adjusted. The initial eigenvalue of the year to be predicted is obtained after completion. At the same time, the spatial weight is introduced, and the pixel $q_1,q_2,\cdots ,q_k$ in the spatial neighborhood of each pixel p with unreasonable change is found. The weighted average of the eigenvalues of these neighborhood pixels is calculated to obtain the eigenvalues closest to reality. The spatial weight matrix is shown in Formula (4):

$$\widehat{\sum }\mathrm{p}={\displaystyle \frac{{\sum }_{\mathrm{i}=1}^{\mathrm{k}}{\mathrm{W}}_{\mathrm{i}\mathrm{p}}{\sum }_{\mathrm{q}\mathrm{i}}}{{\sum }_{\mathrm{i}=1}^{\mathrm{k}}{\mathrm{W}}_{\mathrm{i}\mathrm{p}}}}$$

(4)

where $\widehat{\sum }\mathrm{p}$ is the eigenvalue of the pixel p after completion, and ${\sum }_{qi}$ is the eigenvalue of the neighborhood pixel q. The complete eigenvalue $\widehat{\sum }$ obtained after prediction is combined with the eigenvectors U and $V^T$ obtained by SVD, and the predicted complete matrix is obtained using the inverse matrix formula, Formula (5):

$${\mathrm{X}}_{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}\mathrm{l}\mathrm{e}\mathrm{t}\mathrm{e}}=\mathrm{U}·\widehat{\sum }·{\mathrm{V}}^{\mathrm{T}}$$

(5)

In the whole process, the time trend and spatial relationship are combined, and the pixels with unreasonable changes are predicted and completed by a linear regression model and calculating spatial weights [46,47,48,49,51,52,53,54]. After completion, the matrix is reconstructed to obtain the land classification result after consistency adjustment.

2.4. HMM_LCT Model

2.4.1. Hidden Markov Model (HMM)

The HMM_LCT model is proposed based on the HMM principle proposed by predecessors. The HMM is a probability-based statistical model that describes the process of randomly generating a sequence of observations from a hidden Markov chain. The model can be used to model time series data and infer unknown states from a series of observations. The HMM contains three elements: initial probability A, observation probability B, and state transition probability matrix $\pi$ [42]. The product of observation probability and spatial factor is regarded as the observation probability B of the model. In the process of initialization, recursion, and inverse algorithm of the HMM algorithm, the observation probability is necessary, that is, the spectral observation value generated by a certain class. The land cover transfer probability matrix is calculated as the state transition probability matrix $\pi$, and the HMM_LCT model is formed. The forward–backward algorithm is used to calculate the possibility of pixels belonging to various categories, and the initial classification results are adjusted.

2.4.2. Local Binary Pattern (LBP)

The LBP is a description method for local texture features that compares the neighborhood of each pixel, converts the comparison results into binary numbers, and then connects these binary numbers together to form LBP images. The main steps of the LBP algorithm include: (1) For each pixel in the image, the 8 pixels around it are compared with the center pixel, and the comparison result is converted to a binary number; (2) The obtained 8 binary numbers are sequentially connected together to form an 8-bit binary number, which is the LBP value of the pixel; (3) The above process is performed for each pixel in the image, and finally, an LBP image is obtained, and each pixel of the image is the LBP value corresponding to the pixel in the original image [55,56].

2.4.3. Spatial Factor

When dealing with image classification, local texture feature extraction is an important task. In this respect, the LBP operator is a classical method [56]. This paper aims to improve the LBP operator by introducing spatial confidence factor I so as to strengthen the spatial continuity of classification results. Specifically, we assume that A represents the category label of a pixel P in a certain period, and A1, A2, A3, A4, A5, A6, A7, and A8 are the category labels of 8 neighborhood pixels of P, respectively. We divide the pattern into “1” and “0” based on the class label difference between the neighborhood pixel and the center pixel, and then calculate n based on these binary codes. Where n is the number of pixels in the neighborhood that have the same label as P. Next, we calculate the ratio of the number of occurrences of a certain category in a 3 × 3 window to the size of the window as a spatial factor.

The spatial factor is calculated as shown in Formula (6):

$$\mathrm{I}=\mathrm{n}/9$$

(6)

The values of n are 0,1,2,3,4,5,6,7,8,9, which are divided into 10 levels. But there are special cases, when n = 0, I = 0.5/9, because in the calculation of the HMM joint probability, this factor will be used as a multiplier; if I = 0, then the joint probability directly takes the lowest value, which will affect the final output result [57,58]. By introducing a spatial confidence factor I, we can more accurately consider the spatial relationship between different pixels in an image. This helps us improve classification accuracy, especially when dealing with images with spatial continuity.

2.4.4. HMM_LCT

Since the HMM in previous studies mostly uses expert knowledge as the state transition probability matrix, the parameters are mostly fixed values, and the subjectivity is strong. This paper introduces LCT and applies it to the state transition probability matrix in the HMM to build the HMM_LCT model. In the HMM, the initial probability is the probability of each state (class) occurring at the first time (year). In the method described in this paper, the initial probability is determined using the distribution of classification results for the initial year. Specifically, assuming that there are k categories, we use SVM classifiers to classify the study area and calculate the frequency $C\left(x_i\right)$ of each category at the first moment. Then, the initial probability of the category at the first moment can be expressed as $C\left(x_i\right)/C$, where C is the total number of pixels in the study area [59]. The initial probabilities of each category in this paper are shown in

Table 3. Initial probability distribution.

Initial Distribution	Impervious Surface	Water	Forest and Grass	Plowland
First Year Pixel Ratio	0.23	0.09	0.19	0.49

.

Observation probability is a necessary input in the prediction problem of the HMM, and it is an important part of calculating forward probability and backward probability. Observation probability refers to the probability that a specific symbol is observed in a given state, that is, the conditional probability $P\left(x|{\omega }_t\right)$ of spectral observations generated by a certain category in a given state, where ${\omega }_t$ represents the possible state of t at a given time and x represents the corresponding spectral vector. Since the initial time series classification diagram is generated by the spectral observation sequence and classifier, the observation probability can be obtained indirectly by using the classification process and output and applied to the calculation process in the forward and backward algorithms. Like the observed probability, the spatial factor I is obtained separately at each time, so the product of the observed probability and the spatial factor is regarded as the observed probability of the model. The three elements of the model are improved to $\mathsf{\lambda }=(\mathrm{A},\mathrm{B}\times \mathrm{I},\mathsf{\pi })$, and the HMM_LCT model is constructed [60,61].

$$P\left(x|{\omega }_t\right)={\displaystyle \frac{P\left({\omega }_t|x\right)P\left(x\right)}{P\left({\omega }_t\right)}}$$

(7)

where $P\left({\omega }_t|x\right)$ is the (equivalent) posterior probability, $P\left({\omega }_t\right)$ is the prior probability, and $P\left(x\right)$ is a fixed value for the pixel.

The transfer probability matrix in this model is calculated by the initial land cover classification map (Table 1). Finally, the forward–backward algorithm is used to predict the temporal land cover classes of all pixels in the study area, which is the reasoning process of the model.

2.5. Spatio-Temporal Window Filtering Model

Spatio-temporal window filtering is a method to improve land cover type by combining temporal and spatial information. The basic idea is to use the land cover transfer probability matrix to further refine and modify the preliminary classification results to improve the accuracy and reliability of land cover classification. The transfer probability matrix is obtained by the land cover map and 3 × 3 window and the threshold value is determined so as to screen out the unreasonable change in pixels. The method then needs to iteratively scan the land cover type for each pixel and modify it for unrealistic changes [40].

Specifically, the probability of land cover transfer from the beginning of year t to year t + 1 is calculated, as shown in Formula (8), and then compared to a predefined threshold, as shown in Formula (9). If the probability of land cover transfer from year t to year t + 1 is less than the predefined threshold, the land cover class for that year is modified to the one with the greatest joint transfer probability, as shown in Formula (10).

$$\mathrm{P}(\mathrm{C}(\mathrm{t}+1)={\mathsf{\alpha }}_{\mathrm{j}}|\mathrm{C}(\mathrm{t})={\mathsf{\alpha }}_{\mathrm{i}})={\displaystyle \frac{\sum (\mathrm{C}(\mathrm{t})={\mathsf{\alpha }}_{\mathrm{i}}\mathrm{A}\mathrm{N}\mathrm{D}\mathrm{C}(\mathrm{t}+1)={\mathsf{\alpha }}_{\mathrm{j}})}{\sum \mathrm{C}(\mathrm{t})={\mathsf{\alpha }}_{\mathrm{i}}}}$$

(8)

$$\mathrm{P}\left(\mathrm{C}\left(\mathrm{T}+1\right)={\mathsf{\alpha }}_{\mathrm{j}}|\mathrm{C}\left(\mathrm{t}\right)={\mathsf{\alpha }}_{\mathrm{i}}\right)<\mathrm{\varepsilon }$$

(9)

$$\mathrm{P}\{\mathrm{C}(\mathrm{t})=\mathsf{\alpha }\mathrm{k}|(\mathrm{C}(\mathrm{t}-1),\mathrm{C}(\mathrm{t}+1))\}={\displaystyle \frac{\mathrm{P}\{\mathrm{C}(\mathrm{t})=\mathsf{\alpha }\mathrm{k}|\mathrm{C}(\mathrm{t}-1)\left\}\mathrm{P}\right\{\mathrm{C}(\mathrm{t}+1)|\mathrm{C}(\mathrm{t})=\mathsf{\alpha }\mathrm{k}\}}{{\sum }_{\mathrm{k}=1}^{\mathrm{k}=4}\mathrm{P}\{\mathrm{C}(\mathrm{t})=\mathsf{\alpha }\mathrm{k}|\mathrm{C}(\mathrm{t}-1)\left\}\mathrm{P}\right\{\mathrm{C}(\mathrm{t}+1)|\mathrm{C}(\mathrm{t})=\mathsf{\alpha }\mathrm{k}\}}}$$

(10)

where C(t − 1), C(t), and C(t + 1), respectively, represent the land cover class of year t − 1, t, and t + 1. And ${\alpha }_k(k=\mathrm{1,2},\mathrm{3,4})$ represents the land cover class of the pixel. $\sum \left(C\right(t)={\alpha }_{i}ANDC(t+1)={\alpha }_j)$ is land cover class ${\alpha }_i(i=\mathrm{1,2},\mathrm{3,4})$ in year t, which is converted to the total number of pixels of ${\alpha }_j(i=\mathrm{1,2},\mathrm{3,4})$ in year t + 1; $\sum C\left(t\right)={\alpha }_i$ is the total number of pixels classified as land cover class ${\alpha }_j$ in year t [62].

3. Results Comparison and Analysis

3.1. Initial Classification Result

Since this study is a method of post-classification processing, a relatively simple and mature method was chosen for initial classification. The SVM supervised classification method was used to conduct the initial classification of the study area. In order to improve the classification accuracy, the normalized vegetation index (NDVI) band is superimposed to better assist the classification [63,64]. The initial classification results of Huangpi District and Yaohai District are shown in Figure 6 and Figure 7, respectively.

It can be seen from the classification results that the initial classification results are not ideal. Due to the impact of image quality differences, classification performance in different periods, and terrain fragmentation in non-urban areas of the study area, there is a lot of confusion between cultivated land and impervious surface, resulting in a lot of unreasonable changes between cultivated land and artificial cover. Therefore, there may be misclassification caused by factors such as classifiers and related year images, and irrational classification results generally exist in space.

3.2. Consistency Model Adjusting Results

The results of the HMM_LCT model are shown in Figure 6. Compared with the original classification results, the forest and grass area in 2019 decreased, while the cultivated land area increased. Compared with the original classification results, the HMM_LCT model can better fit the 2019 original image. In the original image, the cultivated land area in 2019 was significantly more than that in 2020 and 2021, so the adjusted results of the HMM_LCT model were more in line with the actual situation. Figure 7 shows the land cover classification map after consistency adjustment using the matrix decomposition model. It can be clearly seen that there is a large amount of plowland information in the initial classification, such as in 2020, but by observing the original image and the existing high-resolution image, it can be seen that there is less bare land area in 2020, and after consistency adjustment, it can be seen that the plowland area is significantly reduced, and the result after consistency adjustment is closer to the original image. The result of the spatio-temporal window filtering model for 2020 is shown in Figure 8. Since the category correction of spatio-temporal window filtering requires the classification results of the previous year and the following year, only the classification results of 2020 in the study area can be adjusted. By comparison, the results of the two models are not very different and are slightly higher than those of the spatio-temporal window filtering model. The comparison of the accuracy of the two models after adjustment is shown in

Table 4. Comparison of accuracy after consistency adjustment improvement.

		2019	2020	2021
SVM (Huangpi District)	Overall Accuracy	0.88	0.85	0.87
	Kappa Coefficient	0.82	0.78	0.81
HMM_LCT	Overall Accuracy	0.89	0.88	0.87
	Kappa Coefficient	0.83	0.82	0.81
SVM (Shunyi District)	Overall Accuracy	0.90	0.86	0.85
	Kappa Coefficient	0.87	0.81	0.80
Matrix decomposition model	Overall Accuracy	0.91	0.91	0.89
	Kappa Coefficient	0.86	0.87	0.85

.

3.3. Accuracy Evaluation and Comparison of Model Results

The confusion matrix of the matrix decomposition model, the HMM_LCT model, and the spatio-temporal window filtering model are shown in

Table 5. Matrix decomposition model result confusion matrix.

Overall Accuracy = 90.74%
Kappa Coefficient = 0.87
	Impervious Surface	Water	Forest and Grass	Plowland	Total	Producer Accuracy	User Accuracy
Impervious Surface	2257	13	20	104	2394	91.04%	94.28%
Water	0	1027	53	20	1100	90.17%	93.36%
Forest and Grass	6	69	2412	93	2580	95.68%	93.49%
Plowland	216	30	36	775	1057	78.13%	73.32%
Total	2479	1139	2521	992	7131

,

Table 6. HMM_LCT adjusts the result confusion matrix.

Overall Accuracy = 89.87%
Kappa Coefficient = 0.86
	Impervious Surface	Water	Forest and Grass	Plowland	Total	Producer Accuracy	User Accuracy
Impervious Surface	1195	4	0	486	1685	87.88%	76.95%
Water	8	1895	35	0	1943	99.53%	97.84%
Forest and Grass	13	10	1325	16	1340	95.35%	99.17%
Plowland	319	3	17	1216	1555	76.55%	86.49%
Total	1535	1912	1377	1723	6547

, and

Table 7. Spatio-temporal window filtering result confusion matrix.

Overall Accuracy = 86.56%
Kappa Coefficient = 0.82
	Impervious Surface	Water	Forest and Grass	Plowland	Total	Producer Accuracy	User Accuracy
Impervious Surface	1231	0	2	505	1740	80.20%	70.75%
Water	0	1895	35	0	1930	99.11%	98.19%
Forest and Grass	3	10	1325	2	1340	96.22%	98.88%
Plowland	184	3	17	1216	1537	70.57%	79.12%
Total	1535	1912	1377	1723	6547

, respectively. First, for the confusion matrix results of the matrix decomposition model results, the overall accuracy is 90.74%, and the Kappa coefficient is 0.87, compared with the initial classification results, which increased by 6%. The model performs well in both producer accuracy and user accuracy in all categories. Especially in the category of impervious surfaces, both producer accuracy and user accuracy exceed 90%, which is a great improvement over the other two matrices. Second, for the confusion matrix results adjusted by the spatio-temporal filtering window model, the overall accuracy is 86.56%, and the Kappa coefficient is 0.82. Similar to the initial classification results, we find that the classifier has a low producer accuracy of 80.20% for impervious surfaces and 70.57% for cultivated land. The accuracy of the forest and grass category is higher, reaching 98.88%. Finally, for the confusion results adjusted by the HMM_LCT model, the overall accuracy is 89.87%, and the Kappa coefficient is 0.86. Especially in the category of plowland, both producer accuracy and user accuracy exceed 86%, which is a great improvement over the previous two matrices.

In general, from the perspective of overall accuracy and kappa coefficient, the adjustment effect of the matrix decomposition model and HMM_LCT model are higher than that of the spatio-temporal window filtering model, and the results of the matrix decomposition model are slightly better than the HMM_LCT model.

3.4. Local Comparison of Model Results

Figure 9 shows the local detail comparison between the initial classification map and the matrix decomposition model result. There is a clear road in the original image; the blurred boundary between roads and vegetation leads to the interruption of roads in the initial classification results. It can be seen from the adjustment results that the misdivided roads have been improved to a certain extent after the consistent adjustment by the matrix decomposition model. Figure 10 shows the local detail map of 2019. Compared with the original image, it can be seen that some pixels in the impervious surface are misclassified as cultivated land. After consistency adjustment by the matrix decomposition model, the misclassified impervious surface area is well improved, proving that the matrix decomposition model has a better effect on consistency adjustment of impervious surfaces.

In the inner part of the yellow box line in Figure 11, the similarity of spectral features caused the SVM classifier to misjudge the cultivated land as impervious water. The HMM_LCT model can make the classification result more accurate by adjusting the wrong cultivated land back to the correct category. In addition, from the results of cultivated land classification, the producer accuracy and user accuracy of the HMM_LCT model are higher than that of the spatio-temporal window filter model, which indicates that the HMM_LCT model has a better adjustment effect on the cultivated land category.

In summary, compared with the spatio-temporal window filter model, the HMM_LCT model performs better in the adjustment of cultivated land. From the local details, the matrix decomposition model can be judged to be relatively better in the consistency adjustment of vegetation, water, and impervious surfaces, and the adjusted results are more in line with the actual situation.

3.5. Land Cover Trajectory Analysis

In this paper, there are four classes of land cover in the annual land cover map, and a total of 64 land cover trajectories were generated [40,65]. Due to the limitation of space, this paper only shows the first 15 trajectories and arranges them in order of proportion. The change trajectory of the initial land cover shows that land cover changes frequently during the study period, especially between forest and cultivated land and cultivated land and impervious surface.

After adjustment by the HMM_LCT model, the unchanged change trajectories of surface object class accounted for the largest proportion in all the change trajectories, ranking the first four and accounting for more than 80% (including 40.62% of cultivated land, 22.76% of forest and grass, 12.40% of impervious surface, and 7.02% of water body). The total proportion of the first 15 change trajectories reached 96.79%, which was much higher than 90% of the initial classification result. For example, impervious surface–cultivated land–impervious surface disappeared in the first 15 change trajectories from 3.2% in the initial classification. This result shows that the HMM_LCT model has achieved a good effect in improving the time consistency of the initial time series land cover results and can improve the time consistency and accuracy of the land cover classification results, thus providing more reliable data support for land use planning and management. The land cover change trajectories adjusted by the HMM_LCT model in Huangpi District are shown in

Table 8. Top 15 land cover trajectories adjusted by the HMM_LCT model.

Number	Change Trajectory	Percentage (%)	Cumulative Percentage (%)
1	plowland–plowland–plowland	40.62	40.62
2	forest–forest–forest	22.77	63.39
3	impervious surface–impervious surface–impervious surface	12.40	75.79
4	water–water–water	7.02	82.81
5	plowland–forest–forest	4.09	86.90
6	impervious surface–plowland–plowland	2.48	89.37
7	plowland–plowland–forest	1.68	91.05
8	plowland–plowland–impervious surface	1.65	92.70
9	impervious surface–impervious surface–plowland	1.21	93.92
10	forest–plowland–forest	0.67	94.59
11	plowland–impervious surface–impervious surface	0.61	95.20
12	plowland–plowland–water	0.53	95.73
13	impervious surface–water–water	0.37	96.10
14	plowland–water–water	0.36	96.47
15	impervious surface–impervious surface–water	0.32	96.79

.

After consistent adjustment by the matrix decomposition model, the cumulative proportion of the first 15 types of land cover trajectories reached 99.47%, of which 83.01% were stable land classes that had not changed for three years, including 51.67% forest and grass, 24.91% impervious surface, 4.68% cultivated land, and 1.76% water body. After the consistent adjustment of the matrix decomposition model, the cumulative proportion of the first 15 types of land cover change trajectory increased and was the highest cumulative proportion. Meanwhile, some land cover classes with unreasonable changes were significantly improved after the consistent adjustment. If the impervious surface–forest–impervious surface trajectory, which accounted for 3% in the initial classification, disappears into the first 15 types of land cover change trajectories after consistent adjustment, it can be shown that the land cover classification change after consistent adjustment is more consistent with the geological logic and practical situation, making the classification results more accurate and reliable and thus improving the availability of the land cover map. The land cover change trajectories adjusted by the matrix decomposition model in Shunyi District are shown in

Table 9. Top 15 land cover trajectories adjusted by the matrix decomposition model.

Number	Change Trajectory	Percentage (%)	Cumulative Percentage (%)
1	forest–forest–forest	51.67	51.67
2	impervious surface–impervious surface–impervious surface	24.91	76.58
3	forest–plowland–forest	5.07	81.65
4	plowland–plowland–plowland	4.68	86.32
5	impervious surface–plowland–impervious surface	4.47	90.79
6	plowland–forest–plowland	2.37	93.16
7	water–water–water	1.76	94.93
8	plowland–forest–forest	1.43	96.35
9	forest–forest–impervious surface	0.84	97.20
10	plowland–impervious surface–plowland	0.72	97.92
11	impervious surface–plowland–plowland	0.44	98.36
12	forest–forest–plowland	0.35	98.71
13	impervious surface–forest–forest	0.28	98.99
14	plowland–impervious surface–impervious surface	0.25	99.24
15	impervious surface–plowland–forest	0.23	99.47

.

3.6. Robustness Test of Two Models

In order to verify the validity of the two models, image sets from different data sources and regions were used in this paper to test the robustness of the model. The GF-1 image in Changping District of Beijing was selected as the experimental area of the HMM_LCT model. The Sentinel-2 image of Yaohai District of Hefei City was selected as the experimental area of the matrix decomposition model. Based on the classification results, the overall accuracy and Kappa coefficient after adjustment by the HMM_LCT model and matrix decomposition model are shown in

Table 10. Accuracy comparison before and after the HMM_LCT model improvement (Gf-1).

		2020	2021	2022
SVM	Overall Accuracy	0.83	0.85	0.84
	Kappa Coefficient	0.76	0.79	0.78
HMM_LCT	Overall Accuracy	0.85	0.87	0.86
	Kappa Coefficient	0.79	0.83	0.81

and

Table 11. Accuracy comparison before and after the matrix decomposition model improvement (Sentinel-2).

		2020	2021	2022
SVM	Overall Accuracy	0.89	0.91	0.90
	Kappa Coefficient	0.85	0.87	0.86
Matrix decomposition model	Overall Accuracy	0.93	0.93	0.91
	Kappa Coefficient	0.89	0.90	0.88

, respectively, and the results show that both overall accuracy and Kappa coefficient have been improved. It is proved that the two models are still applicable in different data sources and regions, and the universality and robustness of the method are proved. The robustness test results of the HMM_LCT model and matrix decomposition model are shown in Figure 12 and Figure 13, respectively.

In addition to the reason of mixed pixels, there are also image quality problems, which lead to differences in the location of surface objects in different years, and these differences are also reflected in the initial classification results. The improvement effect of the matrix decomposition model on spatial consistency is more obvious in the high-score data. As shown in Figure 14, a path was selected to demonstrate the improvement effect of the matrix decomposition model on surface object migration. It can be clearly seen from the perspective window of the original Sentinel-2 image that the path of two years has obvious spatial migration, as shown in Figure 14a. So the initial classification map also has corresponding migration, as shown in Figure 14b. After the adjustment of the matrix decomposition model, the corresponding migration of surface object classification results is significantly improved, as shown in Figure 14c. This result proves that the matrix decomposition model can still play a good role in improving the spatial consistency of the initial time series land cover results using different data sources.

In summary, both models have good adjustment effects, and the effect of the matrix decomposition model is slightly better than that of the HMM_LCT model, which provides certain reference values for the subsequent production and research of land cover products and can also provide accurate information for sustainable land use and resource management. This comprehensive analysis method contributes to a deeper understanding of the dynamic changes in the land system and provides key information for the sustainable development of human society.

4. Conclusions

At present, there are many large-scale land cover products serving users, but the availability of these products varies greatly in different regions. In the process of developing a time series land cover map, it is important to consider more accurate and specific spatio-temporal relationships to improve product quality. This requires the development of new algorithms to capture the spatio-temporal characteristics of land cover change, thereby improving the precision and accuracy of products.

Two consistency adjustment models are introduced and compared with those studied by predecessors. This paper draws the conclusion that the matrix decomposition model can improve the classification results better and demonstrates the effects of the models proposed in this paper from two aspects: land cover change trajectory and confusion matrix. The main novelty of this paper is as follows: In this study, the HMM_LCT model and the matrix decomposition model were used to improve the SVM classification results and compared them with the initial classification results and the spatio-time window filtering model. In this study, the matrix decomposition model is transferred to the consistency adjustment of the land cover classification map, and the spatial weight information is introduced on the basis of the matrix decomposition model so that the model can take into account the spatial distribution of land cover types and the correlation between adjacent pixels. In addition, the traditional HMM is improved, and the subjective influence is reduced by using the existing high-precision land cover products to calculate the transition probability matrix. Moreover, the spatial factor is added to the HMM_LCT model, which makes the model more suitable for consistency adjustment. By calculating the probability matrix of land cover type transfer and introducing spatial neighborhood information, the model can take into account the continuity of land cover type in time and the distribution law in space so that the classification results are more consistent with the law of actual land cover change. Both models improve the migration of surface objects in different years, making the position of surface objects in the temporal land cover classification map more consistent between different years and also making the change detection more accurate.

In general, the consistency adjustment effect of the two models is better than that of the spatio-temporal window filtering model, and the result of the matrix decomposition model is slightly better than that of the HMM_LCT model. At the same time, this study also provides a certain reference value for the research of land cover product production.