Extracting Frequent Sequential Patterns of Forest Landscape Dynamics in Fenhe River Basin, Northern China, from Landsat Time Series to Evaluate Landscape Stability

: The forest landscape pattern evolution can reveal the intensity and mode of action of human–land relationships at different times and in different spaces, providing scientiﬁc support for regional ecological security, human settlement health, and sustainable development. In this study, we proposed a novel method for analyzing the dynamics of landscape patterns. First, patch density (PD), largest patch index (LPI), landscape shape index (LSI), and contiguity index (CI) were used to identify the types of forest spatial patterns. The frequent sequential pattern mining method was used to detect the frequent subsequences from the time series of landscape pattern types from 1991 to 2020 and further evaluate the forest landscape stability of the Fenhe River Basin in China. The results show that different frequent sequence patterns have conspicuous spatial and temporal differences, which describe the evolution processes and stability changes during a certain period of forest evolution and play an important role in the analysis of forest dynamics. The proportion of the disturbed regions to the total forest area exhibited a downward trend. The long-term evolution pattern indicates that there are many evolution processes and trends in the forest at the same time, showing an aggregation distribution law. Compared with 2016, the forest landscape has become complete in 2020, and the overall stability of the Fenhe River Basin has improved. This study can provide scientiﬁc support to land managers and policy implementers and offer a new perspective for studying forest landscape pattern changes and evaluating landscape stability.


Introduction
The forest landscape is composed of forest ecosystems as the main body [1]. Due to the joint action of a variety of ecological processes (such as long-term mining activities, serious soil erosion, and a series of ecological restoration projects), which have a significant impact on landscape function, process, formation, and structure [2][3][4], the forest landscape of the Fenhe River Basin in China changes frequently and has a high degree of temporal and spatial heterogeneity [5][6][7]. The fragmentation, degradation, and loss of forest landscape caused by disturbance may threaten many ecosystems and processes [8] and even limit the diffusion of some species in the ecosystem [9]. The loss of biodiversity reduces the frequency and capacity of ecosystem services and hinders the cycle of the ecosystem [10]. Compared with connected forests, large-scale fragmented forests may be more vulnerable to stress and have lower resilience. These can destabilize the entire system, which can further weaken the resistance of the system to extreme environments and disaster events (floods and droughts) [11,12]. Specifically, the greater the degree of landscape fragmentation, the lower is its stability. Therefore, it is necessary to understand the landscape stability of the Fenhe River Basin by analyzing its dynamic characteristics. the first to identify the association between data record lists. Their method has been widely used in various fields [46][47][48], including remote sensing [49]. However, the potential of this method has not been used in the studies of forest landscape pattern dynamics. Frequent sequential pattern mining has certain advantages because (I) subsequences contain information from temporal and spatial dimensions, such as detailed state, spatiotemporal distribution, trend, and intensity. This information cannot be simultaneously described and quantified by other methods. (II) It can combine patterns and processes that are helpful in explaining and understanding the relationship between patterns and processes. (III) Noise can be removed to only focus on the main information. (IV) The most important evolutionary processes on different temporal scales can be mined without supervision.
This study proposes a new method to detect dynamic changes in landscape patterns by first identifying the types of landscape patterns. The frequent sequential pattern mining method was used to explore the evolution process that frequently occurs in the time series to find the information that users are interested in. Notably, the sequence pattern in the result is not given in advance but is determined by the method itself. This study has three objectives: (I) to identify the main types of forest landscape patterns based on multiple landscape indicators and analyze the spatial characteristics as well as inter-annual trends of landscape patterns; (II) to extract a typical sequence pattern that frequently occurs in the forest landscape patterns and analyze the disturbance, recovery characteristics and dynamic trends of the forest ecosystem in the study area from 1991 to 2020; and (III) to reveal the temporal and spatial distribution characteristics as well as the change characteristics of the stability of the forest landscape in the Fenhe River Basin based on the types of landscape patterns in the spatial scale and the frequent sequence patterns in the time scale.

Study Area
The Fenhe River Basin is located in northern China. Being the second-largest tributary of the Yellow River, the mainstream of the Fenhe River system is 713 km in length, with a basin area of 39,721 km 2 . It traverses across six cities: Xinzhou, Taiyuan, Jinzhong, Luliang, Linfen, and Yuncheng, and 45 counties (cities, districts). The Fenhe River Basin has a temperate continental climate and belongs to an inland semi-arid area. The annual rainfall is 360 to 700 mm, mainly concentrated between June, July, and August. The mean temperature is 4.2 to 14.2 • C [50]. The terrain of the entire Fenhe River Basin is high in the north and low in the south (Figure 1b). Additionally, the middle zone of the Fenhe River Basin is covered by loess soil (with uneven thickness, undulation, and crisscross ravines), which is a typical landform of the Loess Plateau, and has experienced serious soil and water loss. The Fenhe River Basin is one of the largest coal production zones of China [51]. Mining-related activities have damaged the ecological environment, such as causing forest degradation that is difficult to recovery [52]. Therefore, the government has adopted a series of ecological restoration projects to increase forest coverage. A large amount of agricultural land was converted into economic forest. In the past 30 years, the forest landscape pattern of the Fenhe River Basin has changed frequently with social and natural development [7].

Methodology
The framework of the research method used in this study is shown in Figure 2. First, landscape metrics were calculated according to the land cover map obtained by a random forest classifier on the Google Earth Engine (GEE) platform. Second, the four landscape indices were used to identify the three spatial heterogeneous patterns using the Self Organizing Maps (SOM)+K-means two-stage classification method, and the time series was constructed. Then, the operation rules of frequent sequence pattern mining were defined, frequent sequence patterns under different parameters were mined, and their spatiotemporal variation characteristics were analyzed. "Landscape stability" refers to the ability of a landscape to maintain its states and quickly recover after being disturbed by external forces [53]. By referring to a report by Mukherjee et al. [28], this study defines "landscape stability" from the spatial perspective as the weaker forest landscape stability with greater fragmentation as well as the stronger forest landscape stability with less fragmentation. From the time-scale perspective, the stability of forest landscapes with a high frequency of disturbance and recovery events is weaker, but the stability of those with a low disturbance and recovery frequency is stronger. Finally, based on the types of landscape patterns in the spatial scale and the frequent sequence patterns in the temporal scale, this study evaluates the stability of the forest landscape patterns in the Fenhe River Basin.

Methodology
The framework of the research method used in this study is shown in Figure 2. First, landscape metrics were calculated according to the land cover map obtained by a random forest classifier on the Google Earth Engine (GEE) platform. Second, the four landscape indices were used to identify the three spatial heterogeneous patterns using the Self Organizing Maps (SOM)+K-means two-stage classification method, and the time series was constructed. Then, the operation rules of frequent sequence pattern mining were defined, frequent sequence patterns under different parameters were mined, and their spatiotemporal variation characteristics were analyzed. "Landscape stability" refers to the ability of a landscape to maintain its states and quickly recover after being disturbed by external forces [53]. By referring to a report by Mukherjee et al. [28], this study defines "landscape stability" from the spatial perspective as the weaker forest landscape stability with greater fragmentation as well as the stronger forest landscape stability with less fragmentation. From the time-scale perspective, the stability of forest landscapes with a high frequency of disturbance and recovery events is weaker, but the stability of those with a low disturbance and recovery frequency is stronger. Finally, based on the types of landscape patterns in the spatial scale and the frequent sequence patterns in the temporal scale, this study evaluates the stability of the forest landscape patterns in the Fenhe River Basin.

Data Preparation
The GEE is a free data cloud platform that can be used for planetary scale analysis (https://earthengine.google.org/ (accessed on 25 March, 2021)). Radiation correction with high accuracy was performed on the surface reflectance of the Landsat data archived on the platform. Atmospheric correction eliminated the atmospheric effects due to absorption and scattering. Terrain processing of level 1 accuracy was also performed. We downloaded the annual L1T Landsat remote sensing images of WRS-II numbered 124/34, 125/33-36, and 126/33-36 ( Figure 1c) from 1991 to 2020 on the GEE platform. Landsat-5 TM images from 1991-2000 and 2003-2011 were selected as supplements to reduce the impact of SLC-off gaps on Landsat-7 images. Landsat-7 ETM+ images from 2001 to 2002, 2012 and Landsat-8 OLI from 2013 to 2020 were also selected. To reduce the impact of phenological factors, images of the cutting period were selected from June to September for each year. The image quality was controlled by removing all pixels marked as clouds, shadows, or snow from the BQA (Quality Assessment Band) of the satellite images. Finally, the median value of the annual image collection was selected to obtain the annual image. The number of remote sensing image scenes used per year is shown in Table 1. The

Data Preparation
The GEE is a free data cloud platform that can be used for planetary scale analysis (https://earthengine.google.org/ (accessed on 25 March 2021)). Radiation correction with high accuracy was performed on the surface reflectance of the Landsat data archived on the platform. Atmospheric correction eliminated the atmospheric effects due to absorption and scattering. Terrain processing of level 1 accuracy was also performed. We downloaded the annual L1T Landsat remote sensing images of WRS-II numbered 124/34, 125/33-36, and 126/33-36 ( Figure 1c) from 1991 to 2020 on the GEE platform. Landsat-5 TM images from 1991-2000 and 2003-2011 were selected as supplements to reduce the impact of SLC-off gaps on Landsat-7 images. Landsat-7 ETM+ images from 2001 to 2002, 2012 and Landsat-8 OLI from 2013 to 2020 were also selected. To reduce the impact of phenological factors, images of the cutting period were selected from June to September for each year. The image quality was controlled by removing all pixels marked as clouds, shadows, or snow from the BQA (Quality Assessment Band) of the satellite images. Finally, the median value of the annual image collection was selected to obtain the annual image. The number of remote sensing image scenes used per year is shown in Table 1. The auxiliary data included the NASA Digital Elevation Model (NASADEM) and high-resolution historical images from Google. The download and processing of all Landsat data and digital elevation models (DEMs) were performed on the GEE platform. The forest distribution map is the basic data used to calculate the landscape index. The land cover types in the study area were divided into six categories: forest, water, cultivated land, construction land, grassland, and unused land [32,54]. This study defined the Fenhe River Basin forest as land area having a canopy density >30% [32]. The random forest classifier [55] encapsulated on the GEE platform was used to generate land cover maps from 1991 to 2020. All the bands of the Landsat images were input into a random forest classifier for prediction and training samples by visual interpretation based on Landsat images. Approximately 1000 points were marked each year, of which at least 500 were forest samples, and the remaining sample data were used to classify other land cover types. The accuracy of the classification was evaluated using the overall accuracy and kappa coefficient. The accurate criteria of the classification data were that the overall accuracy was >90%, and the Kappa coefficient was >85%. Finally, the forest mask corresponding to the land cover map from 1991 to 2020 was selected as the study area.

Calculating Forest Landscape Metrics
As the calculation of many landscape metrics is based on various statistical processing of patch circumference and area, redundancy is inevitable [56]. We selected four landscape metrics with less correlation: patch density (PD), largest patch index (LPI), landscape shape index (LSI), and contiguity index (CI). These four metrics can be referred to as the "FRAGSTATS Help" [57]. To select an appropriate spatial scale for the forest landscape pattern, we calculated the scale variance of patch density (from 11 × 11 pixels to 111 × 111 pixels (10 pixels between each scale)) and each pixel (with a grid cell size of 30 m × 30 m) [58][59][60]. The smallest scale variance was obtained from 21 × 21 pixels. Therefore, in this study, 21 × 21 pixels was the most appropriate grid cell size. FRAGSTATS 4.2 is a software that quantitatively analyzes landscape structure composition and spatial pattern and can calculate hundreds of landscape metrics [57]. We calculated four landscape metrics using the software; in total, 95712 effective grid cells were calculated.

Extracting Landscape Pattern Types
SOM can reveal the potential distribution patterns and relationships among multiple variables and transform the complex nonlinear statistical relationships between highdimensional data into simple geometric relationships and show them in a low-dimensional manner [61]. Therefore, we used SOM, which has been proven effective in extracting landscape pattern types in earlier studies [29,62]. Four landscape metrics (PD, LPI, LSI, and CI) were used as input data to construct a feature vector of x = [PD, LPI, LSI, CI]. Notably, the output layer of the SOM calculates the minimum distance of the feature vector, finds the rules, and classifies them. The SOM algorithm references the SOMPY library in Python (https://github.com/sevamoo/SOMPY (accessed on 25 March 2021)), and the SOM algorithm conducts preliminary clustering of massive sample data; feature vectors having similar features are classified into one category. Based on this, the output results of the SOM are taken as the input of the K-means algorithm for further clustering. In this study, the Davies-Bouldin index (DBI) was selected to evaluate the effect of clustering [62]. We selected the optimal cluster number by considering the minimum similarity and maximum heterogeneity within and between the regions, respectively [29,30]. In our study, different labels were assigned to landscape types according to the degree of fragmentation. Notably, we found that landscape fragmentation and the stability of the landscape are inversely proportional to each other [28].

Extracting Frequent Sequential Pattern
Frequent sequential pattern mining refers to the process of finding frequent subsequences (as patterns) from a time-series dataset composed of multiple time series [45]. Furthermore, a candidate set composed of different subsequences is generated, wherein each sequence is ordered by repeatable labels (the landscape pattern type in this study). When a threshold of the occurrence frequency of a subsequence is given, all subsequences that exceed this threshold are called frequent sequence patterns. In this study, we used frequent sequence patterns to identify the temporal relationship between landscape pattern types.
Before data mining, we built a dataset of n time series based on the landscape pattern types data. Frequent sequential pattern mining is divided into two processes: the first step is to generate a candidate database, and the second step is filtration ( Figure 3). First, we used the sliding window method to mine the candidate sequence patterns (subsequences). The sliding window is a time window model for common data pattern mining. In this method, the old data are eliminated as the new data enters the window; thus, the internal data of the window are constantly changing. Based on the landscape pattern data set of this study, a sliding window traversed each time series in the data matrix, by automatically adjusting the size of the sliding window. Therefore, each time the window slid, its window sequence formed a candidate pattern. By constantly changing the sliding window size from 1 to 30, all sequence patterns of different lengths were stored to generate a set of candidate sequences. Then, we counted the frequency of candidate sequences, calculated the relative support (Formula (1)) and confidence (Formula (2)) of each sequence pattern, and filtered the patterns that were greater than the minimum relative support (α rel ) and minimum confidence (β) thresholds. These sequential patterns were classified as frequent sequential patterns and are considered more valuable than other patterns. the rules, and classifies them. The SOM algorithm references the SOMPY library in Python (https://github.com/sevamoo/SOMPY (accessed on 25 March, 2021)), and the SOM algorithm conducts preliminary clustering of massive sample data; feature vectors having similar features are classified into one category. Based on this, the output results of the SOM are taken as the input of the K-means algorithm for further clustering. In this study, the Davies-Bouldin index (DBI) was selected to evaluate the effect of clustering [62]. We selected the optimal cluster number by considering the minimum similarity and maximum heterogeneity within and between the regions, respectively [29,30]. In our study, different labels were assigned to landscape types according to the degree of fragmentation. Notably, we found that landscape fragmentation and the stability of the landscape are inversely proportional to each other [28].

Extracting Frequent Sequential Pattern
Frequent sequential pattern mining refers to the process of finding frequent subsequences (as patterns) from a time-series dataset composed of multiple time series [45]. Furthermore, a candidate set composed of different subsequences is generated, wherein each sequence is ordered by repeatable labels (the landscape pattern type in this study). When a threshold of the occurrence frequency of a subsequence is given, all subsequences that exceed this threshold are called frequent sequence patterns. In this study, we used frequent sequence patterns to identify the temporal relationship between landscape pattern types.
Before data mining, we built a dataset of n time series based on the landscape pattern types data. Frequent sequential pattern mining is divided into two processes: the first step is to generate a candidate database, and the second step is filtration ( Figure 3). First, we used the sliding window method to mine the candidate sequence patterns (subsequences). The sliding window is a time window model for common data pattern mining. In this method, the old data are eliminated as the new data enters the window; thus, the internal data of the window are constantly changing. Based on the landscape pattern data set of this study, a sliding window traversed each time series in the data matrix, by automatically adjusting the size of the sliding window. Therefore, each time the window slid, its window sequence formed a candidate pattern. By constantly changing the sliding window size from 1 to 30, all sequence patterns of different lengths were stored to generate a set of candidate sequences. Then, we counted the frequency of candidate sequences, calculated the relative support (formula 1) and confidence (formula 2) of each sequence pattern, and filtered the patterns that were greater than the minimum relative support ( ) and minimum confidence (β) thresholds. These sequential patterns were classified as frequent sequential patterns and are considered more valuable than other patterns. . Schematic of frequent sequential pattern mining. ID is the location of the landscape cell; "T" is the year; "type" is the type of landscape pattern; and "Pattern" is sequence pattern classes. . Schematic of frequent sequential pattern mining. ID is the location of the landscape cell; "T" is the year; "type" is the type of landscape pattern; and "Pattern" is sequence pattern classes.
The support parameter represents the number of occurrences of sequential patterns in the candidate database. If the number of subsequences s 1 = <e 1 , e 2 , . . . . . . , e n > in the candidate database is m, and m is called support, which is recorded as sup(s 1 ), then, relative support is defined as the ratio of the number of s 1 to the total number of subsequences. We used |S| to denote the total number of subsequences in the candidate database. The relative Remote Sens. 2021, 13, 3963 8 of 21 support was recorded as Relsup(s 1 ). For sequential pattern s 1 , the degree of support can be expressed mathematically as follows: The confidence parameter is defined as the ratio of sup(s 1 ) to sup(s 2 ), which is recorded as con(s 1 ). Here, s 1 = <e 1 , e 2 , . . . . . . , e n > and s 2 = <e 1 , e 2 , . . . . . . , e n−1 >. The formula is as follows: The user specifies the minimum relative support and minimum confidence threshold. The principle of specifying a threshold is: the higher the threshold, the fewer patterns the algorithm generates, and lowering the threshold increases the number of patterns (Figure 4c,d). When |S| is certain for a dataset, the minimum confidence can be eliminated. To ensure the number of patterns generated and the effectiveness of this method, we used a strategy of repeated experiments to determine the threshold. Frequent sequential patterns (FSPs) of length "n" were recorded as "n-FSP". Assuming that A and B, respectively, represent the type of landscape pattern. Each frequent sequence pattern in the sequence set is expressed by a label, for example, "AAABB" represents the frequent sequence pattern "A→A→A→B→B". All the programs were implemented in Python. The support parameter represents the number of occurrences of sequential patterns in the candidate database. If the number of subsequences 1 = < 1 , 2 ,……, > in the candidate database is m, and m is called support, which is recorded as sup 1 , then, relative support is defined as the ratio of the number of 1 to the total number of subsequences. We used | | to denote the total number of subsequences in the candidate database. The relative support was recorded as Relsup 1 . For sequential pattern 1 , the degree of support can be expressed mathematically as follows: The confidence parameter is defined as the ratio of sup 1 to sup 2 , which is recorded as con 1 . Here, 1 = < 1 , 2 ,……, > and 2 = < 1 , 2 ,……, 1 >. The formula is as follows: The user specifies the minimum relative support and minimum confidence threshold. The principle of specifying a threshold is: the higher the threshold, the fewer patterns the algorithm generates, and lowering the threshold increases the number of patterns (Figure 4c,d). When | | is certain for a dataset, the minimum confidence can be eliminated. To ensure the number of patterns generated and the effectiveness of this method, we used a strategy of repeated experiments to determine the threshold. Frequent sequential patterns (FSPs) of length "n" were recorded as "n-FSP". Assuming that A and B, respectively, represent the type of landscape pattern. Each frequent sequence pattern in the sequence set is expressed by a label, for example, "AAABB" represents the frequent sequence pattern "A→A→A→B→B". All the programs were implemented in Python.  First, we mined frequent sequential patterns based on the time series of annual images. We tested the relationship between different parameters (the length of frequent sequential patterns, minimum relative support and minimum confidence threshold) and the number of frequent sequence patterns, as shown in Figure 4. Second, to analyze the spatiotemporal distribution of forest landscape pattern evolution, the following indicators were used to evaluate typical frequent sequence patterns: (1) frequency of pattern in the same area; (2) first occurrence time of the pattern; and (3) the proportion of the area that appears every year, that is, the ratio of the landscape cell area covered by the sequence pattern to the total forest landscape cell area.

Short Patterns about Fragmentation and De-Fragmentation
Each change in landscape pattern type is called an "event", and these changes were also observed in 2-FSP. These short patterns offer general information about fragmentation and de-fragmentation. Thus, after excluding stable 2-FSP, we categorized other 2-FSPs as fragmentation and de-fragmentation patterns. To reduce the missed detection rate, we set the relative support and confidence to 0. The fragmented and de-fragmented 2-FSPs were defined as disturbance and recovery events, respectively. The number and time of occurrence of the fragmentation and de-fragmentation events in the time series were recorded, and the area proportion of their occurrences each year was counted.

Long-Term Evolution Process and Trend of Forest Landscape Pattern
Due to the significant differences in the geographical location, development stage, and artificial or natural disturbance of the forests in the Fenhe River Basin, there may be potential differences in the forest evolution processes and distribution characteristics. Because there may be many potential sequential patterns in the 30-year-long time series that the coverage area of each pattern is very small, it is not conducive to interpretation and visualization. Thus, we used the method of selecting time intervals of four or five years, reconstructed the time series from the data of 1991,1996,2001,2006,2011,2016, and 2020, and extracted the frequent sequence patterns. Further, 7-FSP values larger than the minimum support and minimum confidence were visualized. Notably, the direction of landscape change summarized the historical trends of change and reflected the differences in different periods. According to the fragmentation trend, the long-term evolution process was defined as fragmentation, de-fragmentation, stability, and fluctuation.

Spatial Characteristics and Interannual Trends of Different Landscape Pattern Types
First, we used DBI to test the effect of k-means clustering, and the range of the number of clusters was 3-9. According to the minimum value of the DBI for different cluster numbers, we chose the optimal cluster number as 3 (Table 2). Second, landscape cells that are severely fragmented are labeled 1, partially fragmented are labeled 2, relatively complete are labeled 3, and non-forest are labeled 0. Figure 5 shows a comparison of forest landscape cells and high-resolution images from June to September 2020 with different labels. Third, we constructed a collection of time series, based on the labels of these landscape pattern types.   We observed the spatial characteristics and inter-annual trends of different landscape pattern types. As shown in Figure 6, according to each index of landscape in the three types of landscape patterns of the distribution, the severely fragmentized pattern had a larger PD and LSI and a smaller LPI and CI. This indicates that this type of forest landscape cell contains multiple scrawled forest patches, and patches of boundary shape are complicated, or patches with larger internal perforations. Generally, the forest area takes up less of the landscape cell. The average value of each landscape metric in the partially fragmented pattern was between a severely fragmented and a relatively complete pattern, and this category is primarily composed of large patches mixed with small patches. The relatively complete pattern had a larger LPI and CI and a smaller PD and LSI. The forest patches they contained were relatively complete, and almost all of them were covered by forests. As shown in Figure 7, the annual area of the landscape cell category (from 1991 to 2020) was determined. We observed the spatial characteristics and inter-annual trends of different landscape pattern types. As shown in Figure 6, according to each index of landscape in the three types of landscape patterns of the distribution, the severely fragmentized pattern had a larger PD and LSI and a smaller LPI and CI. This indicates that this type of forest landscape cell contains multiple scrawled forest patches, and patches of boundary shape are complicated, or patches with larger internal perforations. Generally, the forest area takes up less of the landscape cell. The average value of each landscape metric in the partially fragmented pattern was between a severely fragmented and a relatively complete pattern, and this category is primarily composed of large patches mixed with small patches. The relatively complete pattern had a larger LPI and CI and a smaller PD and LSI. The forest patches they contained were relatively complete, and almost all of them were covered by forests. As shown in Figure 7, the annual area of the landscape cell category (from 1991 to 2020) was determined.

Spatiotemporal Distribution of Typical Frequent Sequential Patterns
Under the constraint of minimum relative support and minimum confidence of 0.03% and 60%, respectively, we obtained 218 sequential patterns. We considered the sequential patterns "2233", "22111", and "110011" as examples because they correspond to the three most common and typical trends in the process of forest change. The relative support and confidence of pattern "2233" were 0.06% and 94.72%, respectively. The relative support and confidence of pattern "22111" were 0.06% and 82.78%, respectively. The relative support and confidence of pattern "110011" were 0.036% and 73.71%, respectively. Generally, the sequence patterns of the downward trend accounted for a small proportion. Notably, the fluctuant trend had the largest number of sequential patterns, and the relative support of each fluctuating sequential pattern was very small.
The pattern "2233" represents that the forest landscape pattern has changed from a partial fragmentation pattern to a relatively complete pattern in the third year, indicating that the local forest has become lush during this period, indicating forest restoration. This pattern has a large number of sequence patterns with time and a large coverage area (Figure 8a). The pattern "22111" changed from a partial to a serious fragmentation pattern in

Spatiotemporal Distribution of Typical Frequent Sequential Patterns
Under the constraint of minimum relative support and minimum confidence of 0.03% and 60%, respectively, we obtained 218 sequential patterns. We considered the sequential patterns "2233", "22111", and "110011" as examples because they correspond to the three most common and typical trends in the process of forest change. The relative support and confidence of pattern "2233" were 0.06% and 94.72%, respectively. The relative support and confidence of pattern "22111" were 0.06% and 82.78%, respectively. The relative support and confidence of pattern "110011" were 0.036% and 73.71%, respectively. Generally, the sequence patterns of the downward trend accounted for a small proportion. Notably, the fluctuant trend had the largest number of sequential patterns, and the relative support of each fluctuating sequential pattern was very small.
The pattern "2233" represents that the forest landscape pattern has changed from a partial fragmentation pattern to a relatively complete pattern in the third year, indicating that the local forest has become lush during this period, indicating forest restoration. This pattern has a large number of sequence patterns with time and a large coverage area (Figure 8a). The pattern "22111" changed from a partial to a serious fragmentation pattern in the third year, indicating forest degradation. The most frequent occurrence of this pattern was on the edge of the forest and near the road (Figure 8b). Pattern "110011" indicates that the level of the landscape pattern of forest decreased during a certain period, although the level returned to a seriously fragmented pattern in the fifth year. The pattern is primarily distributed in the middle of the Fenhe River Basin and the forest edge (Figure 8c).
The pattern "2233" showed a maximum of 2.52% in 2004, with an annual increase till 2011 (Figure 9a). The year with the highest frequency of the pattern "22111" was 2002, resulting from 3.49% (Figure 9b). The proportion of pattern "110011" was less in most years (Especially in 2010) and more in 1994, 1998, and 2004 (Figure 9c).
In 1991, pattern "2233" first appeared within the forest region and then, at the edge of the forest (Figure 10a). This indicates that the forest stability has improved during a certain period, and the process of improvement first occurred in the core area of the forest and then, at the edge of the forest. Pattern "22111" shows that forest stability declined during a certain period, most of which occurred before 2000 (Figure 10b). Pattern "110011" shows that the forest stability decreased first and then, increased in a period; for the first time, the pattern mostly occurred before 2000 (Figure 10c).
the third year, indicating forest degradation. The most frequent occurrence of this pattern was on the edge of the forest and near the road (Figure 8b). Pattern "110011" indicates that the level of the landscape pattern of forest decreased during a certain period, although the level returned to a seriously fragmented pattern in the fifth year. The pattern is primarily distributed in the middle of the Fenhe River Basin and the forest edge (Figure 8c).  In 1991, pattern "2233" first appeared within the forest region and then, at the edge of the forest (Figure 10a). This indicates that the forest stability has improved during a certain period, and the process of improvement first occurred in the core area of the forest and then, at the edge of the forest. Pattern "22111" shows that forest stability declined during a certain period, most of which occurred before 2000 (Figure 10b). Pattern "110011" shows that the forest stability decreased first and then, increased in a period; for the first time, the pattern mostly occurred before 2000 (Figure 10c). the level returned to a seriously fragmented pattern in the fifth year. The pattern is primarily distributed in the middle of the Fenhe River Basin and the forest edge (Figure 8c).  In 1991, pattern "2233" first appeared within the forest region and then, at the edge of the forest (Figure 10a). This indicates that the forest stability has improved during a certain period, and the process of improvement first occurred in the core area of the forest and then, at the edge of the forest. Pattern "22111" shows that forest stability declined during a certain period, most of which occurred before 2000 (Figure 10b). Pattern "110011" shows that the forest stability decreased first and then, increased in a period; for the first time, the pattern mostly occurred before 2000 (Figure 10c).

Short Patterns about Fragmentation and De-Fragmentation
As shown in Figure 11a, the places with more disturbance frequency are primarily distributed in the border area between the forest and non-forest areas and in the middle of the Fenhe River basin. With respect to time, as shown in Figure 11b, the largest area of disturbance occurred during 1991-2000 (0.11%) and the smallest was during 2016-2017 (0.005%). Overall, the proportion of disturbed forests showed a downward trend from 1991 to 2020. From a spatial perspective, as shown in Figure 11c, the ratio of 1-3 times was the highest (54.3%), and that of more than 6 times was the lowest (5.8%).
As shown in Figure 12, the places with higher frequency were primarily distributed in the area at the junction of the forest and non-forest areas and in the middle of the Fenhe River Basin. As shown in Figure 12b Figure 12c, the proportion of 1-3 recoveries were the highest (52.8%), and the lowest was more than 6 times (7.5%). Remote Sens. 2021, 13, x FOR PEER REVIEW 13 of 23 Figure 10. Gradient overlay map depicting spatial distribution of frequent sequential patterns that appeared for the first time: (a) pattern "2233"; (b) pattern "22111"; and (c) pattern "110011".

Short Patterns about Fragmentation and De-Fragmentation
As shown in Figure 11a, the places with more disturbance frequency are primarily distributed in the border area between the forest and non-forest areas and in the middle of the Fenhe River basin. With respect to time, as shown in Figure 11b, the largest area of disturbance occurred during 1991-2000 (0.11%) and the smallest was during 2016-2017 (0.005%). Overall, the proportion of disturbed forests showed a downward trend from 1991 to 2020. From a spatial perspective, as shown in Figure 11c, the ratio of 1-3 times was the highest (54.3%), and that of more than 6 times was the lowest (5.8%). As shown in Figure 12, the places with higher frequency were primarily distri in the area at the junction of the forest and non-forest areas and in the middle of the River Basin. As shown in Figure 12b  In the past 30 years, the total area of forest recovery events was 749 (10 ha), is larger than that of disturbance events of 621 (10 ha). The location of the recovery was similar to that of the disturbance events, and there was a negative correlation time, with a correlation coefficient of −0.4429.

Long-Term Evolution Process and Trend of Forest Landscape Pattern
For relative support of 0.10% and confidence of 50%, we obtained ten 7-FSP evo patterns ( Table 3). The distribution map of the pattern is visualized as shown in Figu where different colors are used to distinguish different sequential patterns, which sent the evolutionary trajectory of the landscape pattern. Notably, although all land types are forests, the evolutionary history within forests is significantly different tially, the same evolution pattern was gathered. The location was relatively close, may be due to the same artificial or natural interferences in the same period. Temp some patterns showed the transformation from disturbance to restoration. The chan the fragmentation degree of patterns 1 and 3 were relatively stable, and the elevat pattern 3 was higher and distributed centrally. Forests in areas covered by pattern 5, and 6 may disappear. Pattern 2 was primarily distributed in the central part of the River Basin, and it was observed where terraced fields were widely distributed. Pa 7 and 8 represented the areas having increased forest coverage. In the past 30 years, the total area of forest recovery events was 749 (10 4 ha), which is larger than that of disturbance events of 621 (10 4 ha). The location of the recovery events was similar to that of the disturbance events, and there was a negative correlation with time, with a correlation coefficient of −0.4429.

Long-Term Evolution Process and Trend of Forest Landscape Pattern
For relative support of 0.10% and confidence of 50%, we obtained ten 7-FSP evolution patterns ( Table 3). The distribution map of the pattern is visualized as shown in Figure 13, where different colors are used to distinguish different sequential patterns, which represent the evolutionary trajectory of the landscape pattern. Notably, although all land cover types are forests, the evolutionary history within forests is significantly different. Spatially, the same evolution pattern was gathered. The location was relatively close, which may be due to the same artificial or natural interferences in the same period. Temporally, some patterns showed the transformation from disturbance to restoration. The changes in the fragmentation degree of patterns 1 and 3 were relatively stable, and the elevation of pattern 3 was higher and distributed centrally. Forests in areas covered by patterns 2, 4, 5, and 6 may disappear. Pattern 2 was primarily distributed in the central part of the Fenhe River Basin, and it was observed where terraced fields were widely distributed. Patterns 7 and 8 represented the areas having increased forest coverage.  Figure 13. Spatial gradient map depicting the long-term forest evolution pattern in the Fenhe River Basin. As shown in Figure 14, each color represents an evolutionary trend. The area proportion of the fluctuant pattern was 12.9%, stable pattern was 12.1%, and de-fragmentation As shown in Figure 14, each color represents an evolutionary trend. The area proportion of the fluctuant pattern was 12.9%, stable pattern was 12.1%, and de-fragmentation pattern was 3.4%. The fluctuating pattern occupied a dominant position in the Fenhe River Basin and was primarily concentrated in the middle of the basin. The local people have developed many terraces in this area, which are also areas with serious soil and water loss. Notably, most of the stable ones are located in areas with higher elevations and less artificial disturbance. Contrastingly, the de-fragmentation regions are scattered in the transitional zone between the forest and non-forest areas.

Value and Applicability of Frequent Sequential Patterns
Landscape analysis aims to not describe the landscape but explain and understand the process [14]. The method proposed in this study can detect rising, falling, stable, and fluctuant frequent sequential patterns at different temporal scales. Sequential patterns can express information in both temporal and spatial dimensions simultaneously, which is difficult for other methods. The sequential pattern can provide scientific support for forest managers to formulate efficient forest management plans. Pattern "110011" may be the area of returning farmland to the forest because the forest may disappear, and it can be used as an area to prevent forest depletion. The pattern "2233" mostly appeared in the core areas within the forest, which can be the areas where the forest ecosystem continues to be optimized. Spatially, similar to the reports of Oikonomakis, Nikolaos, and Petros Ganatsas [63,64], pattern "22111" mostly occurred at the edge of the forest and near the road, and this process may be caused by artificial land-use change, leading to an increase in fragmentation. Temporally, this pattern occurred most frequently in 2002; the transformation of farmland to forests was implemented in this year [65]. The forest area increased rapidly, resulting in a rapid decline in fragmentation. However, many regions reemerged in 2004; this shows that the impact of returning farmland to forests is short-term, and the maintenance capacity of the forest ecosystem covered by the pattern "22111" was poor. Some research of Oikonomakis Nikolaos, and Petros Ganatsas [63,64] has considered the succession process between different tree species and the driving factors behind them. Elevation was regarded as a key driving factor. This study finds that high altitudes affect the evolution of forests because they hinder the fragmentation of forest landscapes. Forests in high-altitude areas maintain relatively complete and large patches with relatively strong landscape connectivity. Although with greater landscape connectivity and higher forest coverage, such areas may have a rich succession process of tree species within a

Value and Applicability of Frequent Sequential Patterns
Landscape analysis aims to not describe the landscape but explain and understand the process [14]. The method proposed in this study can detect rising, falling, stable, and fluctuant frequent sequential patterns at different temporal scales. Sequential patterns can express information in both temporal and spatial dimensions simultaneously, which is difficult for other methods. The sequential pattern can provide scientific support for forest managers to formulate efficient forest management plans. Pattern "110011" may be the area of returning farmland to the forest because the forest may disappear, and it can be used as an area to prevent forest depletion. The pattern "2233" mostly appeared in the core areas within the forest, which can be the areas where the forest ecosystem continues to be optimized. Spatially, similar to the reports of Oikonomakis, Nikolaos, and Petros Ganatsas [63,64], pattern "22111" mostly occurred at the edge of the forest and near the road, and this process may be caused by artificial land-use change, leading to an increase in fragmentation. Temporally, this pattern occurred most frequently in 2002; the transformation of farmland to forests was implemented in this year [65]. The forest area increased rapidly, resulting in a rapid decline in fragmentation. However, many regions reemerged in 2004; this shows that the impact of returning farmland to forests is short-term, and the maintenance capacity of the forest ecosystem covered by the pattern "22111" was poor. Some research of Oikonomakis Nikolaos, and Petros Ganatsas [63,64] has considered the succession process between different tree species and the driving factors behind them. Elevation was regarded as a key driving factor. This study finds that high altitudes affect the evolution of forests because they hinder the fragmentation of forest landscapes. Forests in high-altitude areas maintain relatively complete and large patches with relatively strong landscape connectivity. Although with greater landscape connectivity and higher forest coverage, such areas may have a rich succession process of tree species within a forest. This is worthy of continued research. Therefore, it is a very meaningful effort to detect the succession process of different forest stands from two dimensions of time and space simultaneously by using frequent sequence patterns in the future. This effort can help understand the environment-succession interactions. However, in the loess hilly area in the middle of the Fen River Basin, there are many small and often disappearing patches, because local people have developed the virgin forests into terraces, thus destroying the connectivity of the landscape. In addition, mining activities have also transformed large forest patches into small ones, thus weakening the landscape connectivity of forests. Although the forest coverage has been increased through ecological restoration projects and the forest landscape has become complete, the conversion of farmland into economic forests or artificially planted single stands may deliver certain negative impacts to the ecological environment [66]. All these driving factors affect the succession of forests as well as movements of species between patches. Jerzy Solon [67] and Lídia S. Bertolo [68] combined fragmentation and important driving factors to analyze landscape changes. Jerzy Solon pointed out that the distance from the traffic route had a significant impact on the landscape structure. They used a comparative analysis method, although they ignored the detailed process between different periods. Our method was able to mine this change process, describing the processes using two dimensions: time and space. In addition, we believe that different drive factors have different effects on the evolution process of the landscape.
Compared with the research of Hermosilla et al. [18]. When considering the overall dynamics of the forest by trend analysis, some instantaneous states are smoothed by "space filtration", making the dynamics of the entire forest look more stable [43,69]. However, in the history of the evolution of the forest, forest dynamics are not always stable. For example, as shown in Figure 10a, pattern "2233" indicates that the stability is enhanced in a certain period and may be affected by the policy. This process has appeared in the northern part of the Fenhe River Basin after 2000. These staged evolutionary processes show the direction of the forest changes. These typical sequence patterns represent landscape stability within a specific period, which helps to strengthen our understanding of the forest landscape.
Forest disturbance and restoration are usually carried out by detecting breakpoints and deviations in the time series composed of forest functions [70,71]. These detection methods pay attention to the time of disturbance occurrence and the trend of the subsequent recovery. Still, they cannot characterize the whole process of disturbance-recovery in both time and space dimensions simultaneously. Although the landscape pattern change has maintained a stable trend for several years, the forest function may still be disturbed before and after the pattern change. The changes in the forest landscape pattern may fall behind in the forest function. Landscape dynamics cannot completely characterize the evolution of forests. Therefore, there have been some studies on the interaction between forest landscape pattern change and function change recently [18,72,73].

Evaluation of Landscape Stability
Areas that are disturbed frequently or change rapidly may lead to the extinction of species and affect local biodiversity [72]. As shown in Figures 11a and 12a, the forests in most areas have been disturbed or restored 1-3 times, and the rest periods have kept a constant spatial pattern. This shows that once stable conditions are established, the forest will enter another state and last for several years. More than six times are distributed at the junction of forest and non-forest, which indicates that human activities have a greater driving force on landscape change than nature. We believe that the higher the interference frequency, the more unstable it is. The central area of the Fenhe River basin should also be the key area to prevent soil erosion. We used the proportion of the disturbed area to the total forest landscape area to express the overall stability of the forest in the Fenhe River basin [43,69]. As shown in Figure 11b, the area proportion of disturbance shows a downward trend, which indicates that the stability of the forest landscape in the basin has improved. Although we did not set a threshold to divide the landscape stability, the trend of the area proportion shows that the overall landscape stability of the basin is increasing.
The long-term evolution pattern expresses the difference in forest fragmentation in different periods, which can be used to analyze detailed process information and changes in stability. As shown in Figure 13, patterns 2, 4, 5, 7, and 8 have become landscapes that are severely fragmented forests from the non-forest area in 2001 or after 2001, which is consistent with the fact that China has fully implemented the conversion of farmland to forest to increase forest cover [65]. We take sequence pattern 2 as an example to explain the evolution process. The formation of pattern 2 may be due to excessive deforestation and land development (for construction) in the 1990s, coupled with frequent natural disasters (drought, pests, and diseases), which led to the deterioration of the ecological environment. The forest began to disappear in 1991. Since 2006, forests have reemerged in the area, indicating that local forest cover increased during 2001-2006 and has since been a stable sequence. With the implementation of environmental awareness and policies, forests have been effectively managed. As shown in Figure 14, similar to the report of Jerzy Solon, a specific feature of landscape pattern change is the simultaneous occurrence of numerous and frequently contradictory processes [67]. Although the fluctuation process is dominant, there are still many areas with de-fragmentation processes, which correspond to regions with increased stability. This asymmetrical development process also occurs within the forest, contributing to the comparison of forest landscapes in different regions.

Research Limitations and Future Work
(1) The frequent sequential pattern mining method is still feasible for detecting other non-forest landscapes and can further explore the pattern changes of different non-forest landscapes.
(2) Although we can adopt measures to reduce errors, such as the SLC-off gaps of Landsat7, there are still challenges in achieving accurate classification in cloudy areas (such as southern China).
(3) Indicators applied to frequent sequence pattern mining are not sufficiently rich. Although Julea et al. [49] detected the spatiotemporal changes in crop species, Normalized Difference Vegetation Index (NDVI), and subsidence/uplift movements related to water level fluctuations on the soil surface, other quantitative indicators or qualitative states are rarely combined with frequent sequential pattern mining methods. Therefore, further research can combine more quantitative indicators and states by exploring the evolution of forests.
(4) The evolution of landscape patterns is the result of the joint action of nature and society. It is necessary to determine the role of different driving factors in landscape change and establish a pattern-ecological process model. Therefore, future research can combine the driving factors (such as fire and logging) to construct time series to analyze the spatiotemporal association rules between driving forces and different indicators, for example, the spatiotemporal correlation between indicators and fire, precipitation, and logging. Mining these frequent event combinations <indicators, drivers> and their sequence patterns are particularly important for understanding the impact of ecological processes on the pattern or function.

Conclusions
In this study, we proposed a new method for detecting dynamic changes in landscape patterns. The temporal and spatial dimensions were considered to describe the evolution of forest landscape patterns simultaneously. Experiments were conducted in the Fenhe River Basin. Furthermore, we evaluated the landscape stability of forests in this area. The findings and conclusions drawn from this study can be summarized as follows.
(1) The frequent sequential pattern mining method is effective in mining the secondary evolution processes that frequently occur in the dynamics of landscape patterns. Shortterm sequential patterns can describe the process of disturbance, recovery, or interference recovery. Furthermore, long-term evolution patterns can describe the trend of forest evolution, state in the process, and landscape stability, which is difficult for other methods to mine simultaneously. The spatiotemporal distribution of sequential patterns is quite different, and its distribution characteristics are helpful for forest managers to design policies related to land degradation and potential environmental pressures.
(2) In the last 30 years, the total area of disturbed forest is less than the total area of forest recovery. There is a negative correlation between the extent of disturbance and recovery events every year, and the proportion of the disturbed regions show a downward trend overall.
(3) A variety of processes occur simultaneously in the same forest ecosystem. According to the evolution process of forest landscape patterns from 1991 to 2020, the same long-term evolution process is mostly aggregated, and the spatial position is relatively close.
(4) Generally, forest stability in the Fenhe River Basin has improved. Relatively complete forest landscape cells have better stability, although they account for a relatively small proportion. Since 2016, its proportion has increased, i.e., more forest landscapes have become complete and stable.