GEE-Based Spatial-Temporal Dynamics in a Ramsar Wetland, Honghe National Nature Reserve, Northeast China from 1985 to 2021

: Wetlands are vital to the human living environment, and with the degradation of wetland ecosystems, it is crucial to protect and restore them. Therefore, based on the long time-series Landsat images provided by Google Earth Engine (GEE), this study obtained the landscape spatial distribution maps of the Honghe National Nature Reserve (HNNR) from 1985 to 2021, analysed the spatial and temporal dynamics of the landscape patterns of the HNNR in the past 40 years, and explored the driving factors of the evolution of the HNNR wetland. The results show that from 1985 to 2008, the HNNR wetlands continued to degrade. The area of the wetland landscape declines and converts mainly into the meadow landscape, and the meadow landscape trends upwards and then downwards and converts mainly into woodland and arable land, with increased fragmentation of wetland and meadow. From 2008 to 2021, with the recovery in hydrological conditions, the area of the wetland landscape increased and fragmentation decreased. However, the meadow landscape continued to decline and fragmentation increased, with meadow converting mainly into wetland; changes in hydrological conditions were the main drivers of the evolution of the HNNR wetlands. The results of this study enable us to better understand the dynamics of the HNNR wetland landscape over the last 40 years and provide assistance for the management of the HNNR wetland ecosystems and the ecological restoration of degraded wetlands.


Introduction
Wetland ecosystems are the most biodiverse ecosystems with unique ecological structures and functions in nature [1], and as a transition zone between land and water [2], they are one of the most important environments for human survival [3]. Since the 20th century, in the context of climate change and rapid population growth, irrational economic development and human activities have had a non-negligible negative impact on the structure and function of wetland ecosystems [4], with declining wetland quantity and quality, degraded ecological functions, declining biodiversity, landscape fragmentation and a growing trend towards desiccation [5]. In ecology, ecological environment, landscape function and landscape structure are closely related, landscape pattern is the embodiment of landscape structure and refers to the arrangement of landscape elements of different sizes, shapes and properties in space [6]. The change of landscape pattern is the result of the interaction between natural and human factors [7]. With the in-depth study of wetland ecology, the study of dynamic changes in landscape patterns and the driving factors has become one of the hot spots [8,9]. Wetland landscape patterns depend on the distribution Land 2022, 11, 2137 2 of 28 and composition of wetland resources, and their changes will obviously affect the ecological service functions of wetland ecosystems such as disturbance resistance, resilience, productivity and biodiversity [10]. Therefore, the study of dynamic changes in wetland landscape patterns is important for revealing the intrinsic factors of wetland landscape changes, predicting future development directions and formulating scientific management strategies [11].
Compared with the traditional field survey methods used in wetland research, RS and GIS technology can quickly obtain real-time large-scale spatial data and can easily extract classification information from remote sensing images, which is an effective means to monitor and study wetlands [12,13], expanding the study of wetland landscapes at both temporal and spatial scales [14,15]. In recent years, Google Earth Engine (GEE) has emerged as a cloud-based computing platform dedicated to processing remote sensing images and other Earth observation data. It has distinctive advantages, especially in long time series and large scale remote sensing detection studies [16], and GEE provides free access to global time series satellite images and other ancillary data (e.g., vector data, digital elevation models, and meteorological data) as well as algorithms for processing large amounts of data with relative ease [17], greatly improving the efficiency of remote sensing scholars. Of these, based on GEE, Random Forest Algorithm is effective for accurate monitoring and mapping of wetlands [18,19].
Sanjiang Plain is home to the largest area of marsh wetlands in China, with marshes accounting for about 35% of the country's wetland area [20][21][22], and are important wetland types that play a fundamental role in purifying water, nourishing water, maintaining biodiversity, regulating regional microclimate, and serving as biological habitats [23][24][25], and hold great ecological value. The Honghe National Nature Reserve (HNNR) is an alluvial sedimentary plain of Sanjiang Plain, which contains almost all the biological species of Sanjiang Plain and plays an important role in protecting rare and endangered plant and animal resources and the wetland ecosystems on which they depend, as well as maintaining biodiversity. The HNNR reflects the full picture of the original wetlands of Sanjiang Plain and is regarded as a microcosm of the original wetlands of Sanjiang Plain [26]-combining typicality, rarity and diversity [27]-and is therefore of great research and conservation value. In fact, large-scale agricultural reclamation and cultivation was performed on Sanjiang Plain in the last century, and the excessive exploitation led to ecological problems such as wetland degradation, degradation of natural vegetation, expansion of soil erosion and sanding, reduction of soil fertility [28], and a sharp decline in the quantity and quality of wetland resources. Under the influence of human activities in the surrounding area, the HNNR wetland ecosystem has been severely degraded, posing a great threat to wetland organisms. Research on the dynamic changes in landscape patterns and the driving factors in the HNNR are therefore crucial to the restoration and conservation of the wetland ecosystems in the area. Fu et al. [29] used multi-frequency, fully polarized SAR image data to identify and classify vegetation in the HNNR using the Random Forest method, and the overall accuracy reached 86.77%. Xie et al. [30] selected multi-temporal Sentinel-1B and Sentinel-2A images of the HNNR and used the Random Forest method to classify the marsh wetlands by remote sensing, and the classification results achieved a maximum accuracy of 93.06% with a Kappa coefficient of 0.916 by parameter tuning. Wan et al. [31] selected remote sensing images from 1989 to 2010 and discussed the effects of land use/land cover change on the wetland landscape of the HNNR, and the results showed that the landscape diversity of the HNNR wetland had decreased. Yan et al. [32] found that the upstream cut-off of the Nongjiang River in 1989 led to a decrease in water level in the Nongjiang River section of the HNNR, and that the decrease in water level was most likely to affect the landscape pattern of the HNNR wetland. Wu et al. [33] suggested that water quantity changes were the main driver of wetland evolution in the HNNR. Yan et al. [34] modelled the plant distribution through GIS, and the results showed that water level fluctuations affected the function and structure of the wetland plant communities in the HNNR. Xue et al. [35] found that water level gradient was the primary environmental factor for Land 2022, 11, 2137 3 of 28 the type and distribution of wetland vegetation communities in the HNNR. Fu et al. [36] monitored the spatial and temporal dynamics of wetland vegetation and hydrological disturbance and restoration in the HNNR from 1985 to 2019 based on Landsat images provided by GEE. The results indicated that hydrological changes were the driving factor of wetland vegetation loss and restoration. Ning et al. [37] selected six Landsat images between 1987 and 2017 and predicted using the Markov model modified by linear regression methods, and the results indicated that the area of wetland, water and woodland in the reserve had recovered in 2023.
These research results have laid a good foundation for the HNNR regional research. However, by collating the research results on the HNNR landscape pattern changes, we found that they are mostly based on traditional RS and GIS technologies, which are limited by computer storage, computing power, and data quality. There are problems such as short research time series, early research data and lack of quantitative analysis of wetland landscape and its driving factors, which cannot effectively and comprehensively reflect the current status and change of the HNNR landscape pattern. Therefore, this study is based on Google Earth Engine to obtain Landsat images of the HNNR from 1985 to 2021, extract wetland information and classify the landscapes, analyse the dynamic change of wetland landscape pattern in the HNNR in the past 40 years, and further explore the driving factors such as the hydrological situation and meteorology that affect the evolution of wetlands. This study is of great theoretical value in enriching and improving the research on the dynamic changes of wetland ecological landscape and its driving factors. At the same time, it can provide a decision basis for the management of the HNNR wetland ecosystem and the ecological restoration of degraded wetlands, which is of great practical significance.

Study Area
Honghe National Nature Reserve (HNNR) (47 • 42 18" N~47 • 52 00" N, 133 • 34 38" Ẽ 133 • 46 29" E) is located in the hinterland of Sanjiang Plain in Heilongjiang Province (Figure 1), at the junction of Tongjiang City and Fuyuan City, bordering Qianfeng Farm in the east, Honghe Farm in the west, and Yalu River Farm in the north. In order to protect the wetlands from agricultural reclamation, the HNNR was established in 1984, with a total area of 251 km 2 . The region belongs to the temperate continental monsoon climate, with an annual mean temperature of 1.6 to 1.9 • C, annual precipitation of 400 to 900 mm, with more than 60% of the precipitation concentrated in July to September, average annual evaporation of 900 mm and an annual frost-free period of 114 to 150 d. The reserve has been divided into core, buffer and experimental zones. There are two swampy rivers in the reserve, one is the Nongjiang River, a primary tributary of the Heilongjiang River, with a total length of 116 km, of which 25.7 km in the middle and lower reaches flows through the buffer and experimental areas of the reserve, and the other is the Wolvlan River, the core water source of the reserve, with a total length of 7 km, mainly located in the core area of the reserve.

Remote Sensing Data Sources and Pre-Processing
The HNNR was established in 1984, but there were missing and irreparable images in 1984. Therefore, based on the Landsat series satellite data provided by Google Earth Engine (GEE), this study selected Landsat surface reflectance products during the vegetation growth period from 1985 to 2021, i.e., June to October, and with less than 20% cloud content. Finally, a total of 138 Landsat 5 TM (from 1985 to 1998, from 2003 to 2011), Landsat 7 ETM+ (from 1999 to 2002 and 2012) and Landsat 8 OLI (from 2013 to 2021) images with a spatial resolution of 30 m were obtained to construct the dataset, which has been subject to geometric correction and atmospheric correction. Landsat 5 satellite has stopped acquiring images since 2012 and Landsat 7 satellite lost data strips. We first used the updateMask function to remove the image edge defects and then used the fill function

Remote Sensing Data Sources and Pre-Processing
The HNNR was established in 1984, but there were missing and irreparable images in 1984. Therefore, based on the Landsat series satellite data provided by Google Earth Engine (GEE), this study selected Landsat surface reflectance products during the vegetation growth period from 1985 to 2021, i.e., June to October, and with less than 20% cloud content. Finally, a total of 138 Landsat 5 TM (from 1985 to 1998, from 2003 to 2011), Landsat 7 ETM+ (from 1999 to 2002 and 2012) and Landsat 8 OLI (from 2013 to 2021) images with a spatial resolution of 30 m were obtained to construct the dataset, which has been subject to geometric correction and atmospheric correction. Landsat 5 satellite has stopped acquiring images since 2012 and Landsat 7 satellite lost data strips. We first used the up-dateMask function to remove the image edge defects and then used the fill function provided by United States Geological Survey (USGS) to fill the images with the 2011 Landsat 5 TM image to obtain the filled Landsat 7 ETM+ images, so as to repair the 2012 Landsat 7 ETM+ images.
In order to obtain a high-quality dataset, a series of pre-processing was performed on the dataset using GEE. Firstly, the Clip and Mosaic functions were used to crop and mosaic the study area. Secondly, the pixel quality assessment bands (pixel_qa) generated by The C Function of Mask (CFMask) algorithm were used to identify clouds and cloud shadows for de-clouding. Then, the Normalized Difference Vegetation Index (NDVI), the Modified Normalized Water Index (MNDWI), and the Liquid Surface Water Index (LSWI) were calculated online using the normalized difference function and added them to the images in each collection, with NDVI used to help distinguish dense vegetation including woodland, MNDWI used to distinguish different water body types, and LSWI used to help distinguish out arable land [38]. The median function was applied to the datasets to achieve a median synthesis according to the actual conditions in the study area, and the synthesised images were used for subsequent supervised classification.

Fieldwork Data
The verification of landscape classification results is an essential part of classification. In order to obtain a high-quality dataset, a series of pre-processing was performed on the dataset using GEE. Firstly, the Clip and Mosaic functions were used to crop and mosaic the study area. Secondly, the pixel quality assessment bands (pixel_qa) generated by The C Function of Mask (CFMask) algorithm were used to identify clouds and cloud shadows for de-clouding. Then, the Normalized Difference Vegetation Index (NDVI), the Modified Normalized Water Index (MNDWI), and the Liquid Surface Water Index (LSWI) were calculated online using the normalized difference function and added them to the images in each collection, with NDVI used to help distinguish dense vegetation including woodland, MNDWI used to distinguish different water body types, and LSWI used to help distinguish out arable land [38]. The median function was applied to the datasets to achieve a median synthesis according to the actual conditions in the study area, and the synthesised images were used for subsequent supervised classification.

Fieldwork Data
The verification of landscape classification results is an essential part of classification. The condition of the sampling sites directly influences the judgement of the classification results and is a guarantee of an objective evaluation of accuracy. We have conducted ground surveys in August 2004 and September 2015, collecting 330 samples each year. More specifically, 119, 87, 72, 30 and 22 sample points were collected for wetland, meadow, woodland, water and arable land in 2004, and 93, 103, 90, 29 and 15 points for each category in 2015. In addition, 2/3 sample points were randomly selected as training samples, 1/3 of the sample points were used to verify the accuracy of the classification results. The published literature and Google Earth high-resolution images were also referred to, and the years 1988, 1996, 2002, 2014 and 2019 were additionally selected for accuracy validation to improve the accuracy of the classification results. The spatial distribution of the sample points for the selected years is shown in Appendix A.

Historical Hydrological Data
The Odyssey water level recording system was used to automatically monitor the hydrological changes in 2012, and the average water level in 2012 was calculated based on

Climatic Data
ERA5 global atmospheric reanalysis data from the European Centre for Medium-Range Weather Forecasts (ECMWF) is widely used for meteorological monitoring and the dataset provides the best estimate of the global atmosphere-land-ocean system and benefits from the absence of spatial and temporal gaps in the record. It has been shown [39] that the ERA5 dataset is reliable for meteorological monitoring in temperate regions. In this study, as ERA5 lacks data for 2020 and 2021, we selected the annual mean temperature and annual precipitation at 2 m height for the HNNR from 1985 to 2019, a dataset that is directly callable on GEE.

Random Forest Algorithm
Random Forest is a combination of a set of randomly trained decision trees. The key to its model is that the trees that make up the forest are random to each other. This allows the predictions between the decision trees to be de-correlated, which in turn improves generalisation and robustness [40]. Currently, the Random Forest algorithm is much less computationally intensive than other neural network computational algorithms and the accuracy of the results is guaranteed, so it is commonly used in large area remote sensing image classification, forest fire early warning, forest resource survey, etc. It has a more important role in geographic country monitoring and land classification coverage. In our study, image classification was conducted on GEE using a combination of the Random Forest algorithm and visual interpretation, and we selected years of fieldwork data, Google Earth high-resolution images and published literature for use in aiding visual interpretation. Five landscape types were identified: wetland landscape represented by Carex lasiocarpa, Carex pseudocuraica, meadow landscape represented by dwarf shrubs such as Deyeuxia angustifolia, Salix rosmarinifolia and Spiraea salicifolia, as well as shrubs and grasses, woodland consisting of Betula platyphylla and Quercus mongolica Fisch, etc., arable land containing paddy fields and dry fields, and open water including rivers and small lake bubbles. After several corrections, the final spatial distribution map of the HNNR landscape patterns from 1985 to 2021 was obtained. Based on the existing limited ground survey, published literature and high-resolution images, we assessed the accuracy of the images for 1988,1996,2002,2004,2014,2015 and 2019 only, and concluded that the average overall accuracy was greater than 90% and the Kappa coefficients were all greater than 0.83, indicating the usability of the corrected results.

Dynamic Change Transfer Matrix
Dynamic change transfer matrix is a common method used in land use/land cover and landscape patterns change studies. Its main function is to quantify the increase and decrease in change between landscape types at different times, as well as the source and direction of change and transfer. The overlay analysis was conducted using ArcGIS to produce landscape type shift maps for different landscape types in the HNNR during the study period. The formula of the transfer matrix is as follows: where A represents the area of each landscape type, i and j represent the landscape types at the beginning and end of the study period, respectively, and n refers to the total number of landscape types.

Landscape Indices Selection and Calculation
The analysis of wetland landscape patterns often requires the use of various quantitative indicators for the study and evaluation of landscape structure. A landscape index is a quantitative indicator that reflects its structural composition and spatial configuration. It better reflects the direction, rate and spatial variation of wetland landscape evolution and makes it easy to disentangle the driving factors that cause this evolution. The indices used for landscape patterns analysis are divided into three main types: patch level (P), class level (C), and landscape level (L), with the number of indices reaching more than 50. Currently, Fragstats 4.2 allows the calculation of these indicators. However, some indices are so similar to each other that it is not appropriate to use too many indices, especially those with a high degree of repetition. Based on various analyses and optimisations of landscape indicators selected by previous authors and taking into account the field conditions in the study area, four indices, Area_MN, COHESION, LPI and PLAND, were selected at class level, and another four indices, AI, IJI, SHAPE_MN and SHDI, were selected at landscape level for the analysis of landscape changes in the study area in each year. The landscape indices used for this study are shown briefly in Table 1 and the specific meanings are given in Appendix B. Table 1. Summary information about the landscape indices used in this study.

Name
Abbreviation Level Computing Equation Interspersion and juxtaposition index IJI

Correlation Analysis
The wetland landscape area in this study was analysed as a non-normally distributed continuous variable, so we used Spearman's correlation analysis in SPSS 26.0 to explore the correlation between the HNNR wetland landscape area and meteorological data.
The flow chart of the method used in this study is shown in Figure 2.

Correlation Analysis
The wetland landscape area in this study was analysed as a non-normally distributed continuous variable, so we used Spearman's correlation analysis in SPSS 26.0 to explore the correlation between the HNNR wetland landscape area and meteorological data.
The flow chart of the method used in this study is shown in Figure 2.

Spatial-Temporal Changes in the HNNR Landscape
The spatial distribution of the landscape for each five-year period from 1985 to 2021 ( Figure 3) is selected in this paper as a display, and the year-by-year spatial distribution of the landscape for different periods is shown in Appendix C, and the year-by-year area change is shown in Figure 4. More specifically, 2008 was the year when a dam was built on the Nongjiang River downstream in the reserve, together with the annual increase in precipitation from 2008 to 2021, resulting in more surface water storage, higher water level and a certain degree of hydrological recovery in the HNNR. The year 2008 is therefore taken as a turning point in the evolution of the HNNR landscape, and this paper analyses the evolution of the HNNR landscape in two phases: from 1985 to 2008 and from 2008 to 2021.
From 1985 to 2008, the wetland area showed a decreasing trend, dropping to a minimum of 114.25 km² in 2006, a decrease of 19.25 km² compared to 1985. From 2008 to 2021, the wetland showed an increasing trend, rising from 118.28 km² in 2008 to a maximum of 140.7 km² in 2012. After 2012, the area shares of the wetland decreased and remained stable at 52% to 54%. From the spatial distribution map, it is clear that the wetland was mainly located on both sides of the river floodplains formed by the Nongjiang River and the Wolvlan River, as well as in the southern part of the reserve. From 1985 to 2008, the wetland landscape was gradually replaced by other landscapes due to the degradation of the hydrological situation caused by the cut-off of the upper reaches of the Nongjiang River and the reduction of annual precipitation. After 2008, as mentioned above, the recovery of the hydrological situation led to the gradual restoration of the wetland landscape.
From 1985 to 2008, the meadow showed an upward and then downward trend, reaching a maximum of 113.97 km² in 1996, accounting for 45.29% of the HNNR. In 2008, the meadow area decreased by 1.32 km² compared to 1985; and from 2008 to 2021, the meadow showed a more obvious downward trend, with a minimum meadow area of

Spatial-Temporal Changes in the HNNR Landscape
The spatial distribution of the landscape for each five-year period from 1985 to 2021 ( Figure 3) is selected in this paper as a display, and the year-by-year spatial distribution of the landscape for different periods is shown in Appendix C, and the year-by-year area change is shown in Figure 4. More specifically, 2008 was the year when a dam was built on the Nongjiang River downstream in the reserve, together with the annual increase in precipitation from 2008 to 2021, resulting in more surface water storage, higher water level and a certain degree of hydrological recovery in the HNNR. The year 2008 is therefore taken as a turning point in the evolution of the HNNR landscape, and this paper analyses the evolution of the HNNR landscape in two phases: from 1985 to 2008 and from 2008 to 2021.
From 1985 to 2008, the wetland area showed a decreasing trend, dropping to a minimum of 114.25 km 2 in 2006, a decrease of 19.25 km 2 compared to 1985. From 2008 to 2021, the wetland showed an increasing trend, rising from 118.28 km 2 in 2008 to a maximum of 140.7 km 2 in 2012. After 2012, the area shares of the wetland decreased and remained stable at 52% to 54%. From the spatial distribution map, it is clear that the wetland was mainly located on both sides of the river floodplains formed by the Nongjiang River and the Wolvlan River, as well as in the southern part of the reserve. From 1985 to 2008, the wetland landscape was gradually replaced by other landscapes due to the degradation of the hydrological situation caused by the cut-off of the upper reaches of the Nongjiang River and the reduction of annual precipitation. After 2008, as mentioned above, the recovery of the hydrological situation led to the gradual restoration of the wetland landscape.
From 1985 to 2008, the meadow showed an upward and then downward trend, reaching a maximum of 113.97 km 2 in 1996, accounting for 45.29% of the HNNR. In 2008, the meadow area decreased by 1.32 km 2 compared to 1985; and from 2008 to 2021, the meadow showed a more obvious downward trend, with a minimum meadow area of 71.33 km 2 in 2021, accounting for only 28.4% of the HNNR-a decrease of 30.76 km 2 compared to 2008. According to the spatial distribution map, the meadow on both left and right sides of the Wolvlan River tended to expand towards the centre in the first period and then declined due to the expansion of other landscapes. The evolution of the meadow landscape is largely dependent on the decline or expansion of other landscapes. 71.33 km² in 2021, accounting for only 28.4% of the HNNR-a decrease of 30.76 km² compared to 2008. According to the spatial distribution map, the meadow on both left and right sides of the Wolvlan River tended to expand towards the centre in the first period and then declined due to the expansion of other landscapes. The evolution of the meadow landscape is largely dependent on the decline or expansion of other landscapes.   From 1985 to 2008, the woodland was basically surrounded by the meadow and showed a trend of first decreasing and then increasing. From 2008 to 2021, the woodland area continued to increase and reached a maximum of 20.37% in 2018. Over the years, the woodland area accounted for 3.51% to 10.6% of the total area of the HNNR. From the spatial distribution map, it can be seen that the woodland was distributed on the east side of the Nongjiang River and the northeast corner of the HNNR. The decline in the area of   From 1985 to 2008, the woodland was basically surrounded by the meadow and showed a trend of first decreasing and then increasing. From 2008 to 2021, the woodland area continued to increase and reached a maximum of 20.37% in 2018. Over the years, the woodland area accounted for 3.51% to 10.6% of the total area of the HNNR. From the spatial distribution map, it can be seen that the woodland was distributed on the east side of the Nongjiang River and the northeast corner of the HNNR. The decline in the area of side of the Nongjiang River and the northeast corner of the HNNR. The decline in the area of woodland was mainly caused by agricultural development activities in Sanjiang Plain. Then, a policy of returning farmland to forests was introduced in the 1990s for the protection of woodland, and as a result the area of woodland has since recovered. Arable land accounted for 0.23% to 5.21% of the total area of the HNNR. There was a significant upward trend in both phases, with the main distribution in the southwest corner of the HNNR in the first phase and then gradually spreading to the north. The water was minimally distributed before 2010, while after the completion of the dam construction in 2008, it was concentrated near the dam in the north, accounting for up to 2.36% of the total area of the HNNR during the study period. Water was negligible before 2008, while after 2008 it rose to a maximum of 5.94 km 2 in 2014, then declined and stabilised at around 4.3 to 4.6 km 2 after 2017. The creation of dams has a non-negligible positive impact on the wetland and water.

Transfer Direction of Landscape Types in the HNNR
As described above, this paper divides the change of the transfer direction of landscape types in the HNNR during the study period into two phases, from 1985 to 2008 and from 2008 to 2021.
The relationships between different landscape types in the HNNR from 1985 to 2008 and from 2008 to 2021 are shown in Tables 2 and 3, respectively. From 1985 to 2008, wetland was mainly converted into meadow, followed by arable land and woodland; meadow was converted into wetland, woodland and arable land. Woodland was mainly transformed into meadow. Arable land and water remain largely unchanged. From 2008 to 2021, wetland was mainly transformed into meadow and water; meadow was mainly transformed into wetland, followed by woodland, arable land and water; woodland was mainly transformed into wetland and meadow. Arable land and water remain largely unchanged. It can be seen that the two phases from 1985 to 2008 and from 2008 to 2021 are dominated by the conversion of wetland to meadow and meadow to wetland, respectively, with woodland in both phases being converted mainly from meadow, arable land likewise being converted mainly from meadow followed by wetland, and water being converted mainly from wetland followed by meadow.

Quantitative Analysis of Landscape Patterns of the HNNR
Analysed at the class level, as shown in Figure 5, the overall landscape type and dominance of the HNNR remained constant from 1985 to 2021, with the PLAND value of the wetland decreasing from 45.5% to 56%, showing a decreasing and then increasing trend, but always being the most dominant landscape type; the meadow decreasing from 45% to about 28%, proving that the dominance of meadow was decreasing; while the PLAND values of the woodland, arable land and water showed an increasing trend, indicating that the dominances of these three landscape types were expanding, mainly by reducing the dominance of meadow.

Quantitative Analysis of Landscape Patterns of the HNNR
Analysed at the class level, as shown in Figure 5, the overall landscape type and dominance of the HNNR remained constant from 1985 to 2021, with the PLAND value of the wetland decreasing from 45.5% to 56%, showing a decreasing and then increasing trend, but always being the most dominant landscape type; the meadow decreasing from 45% to about 28%, proving that the dominance of meadow was decreasing; while the PLAND values of the woodland, arable land and water showed an increasing trend, indicating that the dominances of these three landscape types were expanding, mainly by reducing the dominance of meadow. As shown in Figure 6, the LPI value of the wetland consistently far exceeded other landscapes, followed by the meadow, which also had higher LPI value with smaller percentage values for other landscape types. As with the PLAND value, it also shows that the wetland dominated the HNNR.
In addition, the LPI value of the wetland showed a decreasing trend from 1985 to 2008. The Area_MN and COHESION values (Figure 7) slightly decreased, indicating that the wetland became increasingly fragmented, more dispersed and less aggregated during this time period. The LPI, Area_MN and COHESION values showed a clear increasing trend from 2008 to 2021, implying that the degree of wetland fragmentation had improved with less dispersion and higher connectivity.
The LPI, Area_MN and COHESION values of the meadow consistently showed a decreasing trend, indicating continued fragmentation of the meadow and greater landscape dispersion. In contrast, the LPI, Area_MN and COHESION values of woodland, arable land and water were all on an upward trend, implying that fragmentation was not as severe in these three types of landscape and connectivity was becoming higher.
Analysed at the landscape level, as shown in Figure 8, the SHAPE_MN value showed an upward trend indicating that the landscape shape in the HNNR was becoming more complex, which was conducive to wetland stability. The IJI value showed a rising trend, indicating that the overall patch-to-patch connectivity was increasing due to the increased area of smaller occupied landscape types, and that the high degree of mixing and high connectivity of landscape types were conducive to the HNNR wetland stability. The AI As shown in Figure 6, the LPI value of the wetland consistently far exceeded other landscapes, followed by the meadow, which also had higher LPI value with smaller percentage values for other landscape types. As with the PLAND value, it also shows that the wetland dominated the HNNR. value showed an upward trend indicating a greater and then less fragmentation across the landscape in the HNNR. The overall trend was upward, ultimately leading to an increase in spatial homogeneity. The SHDI value showed a rising trend, indicating that the patches of each landscape type were distributed in an increasingly balanced trend in the HNNR. The heterogeneity of the landscape increased and the proportional structure tended to smooth out.     The LPI, Area_MN and COHESION values of the meadow consistently showed a decreasing trend, indicating continued fragmentation of the meadow and greater landscape dispersion. In contrast, the LPI, Area_MN and COHESION values of woodland, arable land and water were all on an upward trend, implying that fragmentation was not as severe in these three types of landscape and connectivity was becoming higher.
Analysed at the landscape level, as shown in Figure 8, the SHAPE_MN value showed an upward trend indicating that the landscape shape in the HNNR was becoming more complex, which was conducive to wetland stability. The IJI value showed a rising trend, indicating that the overall patch-to-patch connectivity was increasing due to the increased area of smaller occupied landscape types, and that the high degree of mixing and high connectivity of landscape types were conducive to the HNNR wetland stability. The AI value showed an upward trend indicating a greater and then less fragmentation across the landscape in the HNNR. The overall trend was upward, ultimately leading to an increase in spatial homogeneity. The SHDI value showed a rising trend, indicating that the patches of each landscape type were distributed in an increasingly balanced trend in the HNNR. The heterogeneity of the landscape increased and the proportional structure tended to smooth out. Land 2022, 11, 2137 13 of 28

GEE-Based Spatial-Temporal Dynamics of the HNNR Landscape Patterns
By analysing studies published before 2008, we found that previous authors believed that the area of wetland declined with the rise of arable land, and that the wetland landscape gradually became more fragmented [31]. Based on the spatial distribution map of the HNNR landscape from 1985 to 2021, we obtained that the wetland and meadow landscapes were severely degraded before 2008, which is consistent with the results of previous studies before 2008 [41]. However, the study concluded that the wetland area decreased by 10% from 1998 to 2005. At the rate of degradation at that time, the wetland would be completely lost in the next 30 years [41], contrary to our result that the wetland area showed an increasing trend after 2008, mainly due to the establishment of the dam on the Nongjiang River in the reserve in 2008, which restored the hydrological situation and the wetland landscape followed by some recovery. The results of our study are consistent with the prediction of an increase in wetland area by 2023 using Markov models [37]. A study by Fu et al. [36] on the evolution of the spatial and temporal dynamics between wetland vegetation and hydrology in the HNNR from 1985 to 2019 concluded that there was a general trend of recovery of wetland vegetation driven by hydrological factors over the last 35 years, which can also provide some reference value for our study and is consistent with our belief that the HNNR wetland are recovering.

Drivers of the Evolution of the HNNR Wetland
Changes in hydrological conditions play a direct and important role in the formation, development, evolution and disappearance of wetlands [42], and are the main abiotic drivers of wetland evolution [43]. Climate change can also lead to wetland evolution, with rising temperatures increasing water evaporation from wetlands and decreasing precipitation leading to shrinkage of water within wetlands, which can lead to wetland degradation [44]. Therefore, climate change and changes in hydrological conditions are analysed in this paper as driving factors of changes in the HNNR landscape dynamics.
Yan et al. [45] showed that the annual mean temperature of Sanjiang Plain has been increasing and annual precipitation has generally decreased. Whereas changes in temperature affect evaporation, a sustained increase in temperature leads to an increase in

GEE-Based Spatial-Temporal Dynamics of the HNNR Landscape Patterns
By analysing studies published before 2008, we found that previous authors believed that the area of wetland declined with the rise of arable land, and that the wetland landscape gradually became more fragmented [31]. Based on the spatial distribution map of the HNNR landscape from 1985 to 2021, we obtained that the wetland and meadow landscapes were severely degraded before 2008, which is consistent with the results of previous studies before 2008 [41]. However, the study concluded that the wetland area decreased by 10% from 1998 to 2005. At the rate of degradation at that time, the wetland would be completely lost in the next 30 years [41], contrary to our result that the wetland area showed an increasing trend after 2008, mainly due to the establishment of the dam on the Nongjiang River in the reserve in 2008, which restored the hydrological situation and the wetland landscape followed by some recovery. The results of our study are consistent with the prediction of an increase in wetland area by 2023 using Markov models [37]. A study by Fu et al. [36] on the evolution of the spatial and temporal dynamics between wetland vegetation and hydrology in the HNNR from 1985 to 2019 concluded that there was a general trend of recovery of wetland vegetation driven by hydrological factors over the last 35 years, which can also provide some reference value for our study and is consistent with our belief that the HNNR wetland are recovering.

Drivers of the Evolution of the HNNR Wetland
Changes in hydrological conditions play a direct and important role in the formation, development, evolution and disappearance of wetlands [42], and are the main abiotic drivers of wetland evolution [43]. Climate change can also lead to wetland evolution, with rising temperatures increasing water evaporation from wetlands and decreasing precipitation leading to shrinkage of water within wetlands, which can lead to wetland degradation [44]. Therefore, climate change and changes in hydrological conditions are analysed in this paper as driving factors of changes in the HNNR landscape dynamics.
Yan et al. [45] showed that the annual mean temperature of Sanjiang Plain has been increasing and annual precipitation has generally decreased. Whereas changes in temperature affect evaporation, a sustained increase in temperature leads to an increase in saturated vapour pressure at the water surface and a lack of saturation leads to increased evaporation, thus reducing its water capacity and degrading the wetlands [46]. The changes in temperature from 1985 to 2019 in the HNNR that we obtained are shown in Figure 9a, and a linear regression of the annual mean temperature indicates that the HNNR is continuing to warm at a rate of 0.32 • C per 5 years since 1985. The continued increase in temperature accelerates the evaporation process and reduces the water storage capacity of the HNNR wetland ecosystem, with the wetland area perhaps shrinking due to the reduction in water volume. The annual precipitation from 1985 to 2019 in the HNNR is shown in Figure 9b, where we can conclude that precipitation decreased at a rate of 2.23 mm per year from 1985 to 2008, reaching a minimum value of 432 mm in 2008, followed by an increasing trend after 2009. The low annual precipitation in the HNNR, coupled with the yearly decrease from 1985 to 2008, has reduced the water supply for the wetland vegetation, resulting in the degradation of the HNNR wetland habitat. Whereas annual precipitation has increased at a rate of 19.71 mm per year since 2009, the rise in annual precipitation will likewise go some way to replenishing the HNNR wetland vegetation and restoring the wetland habitat, which is consistent with the analysis of previous authors [33]. The HNNR wetland landscape area was analysed to be insignificantly correlated with temperature (r = −0.087, p > 0.05) and significantly correlated with precipitation (r = 0.368, p < 0.05).
Land 2022, 11, 2137 14 of 28 saturated vapour pressure at the water surface and a lack of saturation leads to increased evaporation, thus reducing its water capacity and degrading the wetlands [46]. The changes in temperature from 1985 to 2019 in the HNNR that we obtained are shown in Figure 9a, and a linear regression of the annual mean temperature indicates that the HNNR is continuing to warm at a rate of 0.32 °C per 5 years since 1985. The continued increase in temperature accelerates the evaporation process and reduces the water storage capacity of the HNNR wetland ecosystem, with the wetland area perhaps shrinking due to the reduction in water volume. The annual precipitation from 1985 to 2019 in the HNNR is shown in Figure 9b, where we can conclude that precipitation decreased at a rate of 2.23 mm per year from 1985 to 2008, reaching a minimum value of 432 mm in 2008, followed by an increasing trend after 2009. The low annual precipitation in the HNNR, coupled with the yearly decrease from 1985 to 2008, has reduced the water supply for the wetland vegetation, resulting in the degradation of the HNNR wetland habitat. Whereas annual precipitation has increased at a rate of 19.71 mm per year since 2009, the rise in annual precipitation will likewise go some way to replenishing the HNNR wetland vegetation and restoring the wetland habitat, which is consistent with the analysis of previous authors [33]. The HNNR wetland landscape area was analysed to be insignificantly correlated with temperature (r = −0.087, p > 0.05) and significantly correlated with precipitation (r = 0.368, p < 0.05). Sanjiang Plain has always been an important base for food production in China. Since the 1950s, it has undergone a number of land developments with a rapid rise in the area of arable land, all dominated by the reclamation of swampy wetlands and including deforestation and the construction of water conservancy projects. The HNNR is surrounded by Qianfeng Farm, Honghe Farm and Yalu River Farm, and the surrounding farms have established extensive irrigation systems. Since 1988, seven drainage trunk canals and a flood control dam have been constructed in the upper reaches of the Nongjiang River basin. This cuts off the Nongjiang River, resulting in the diversion of the Nongjiang River water that previously flowed through the HNNR to discharge directly into the Heilongjiang River [47], further cutting off the upper reaches of the Nongjiang River and depriving it of its fluvial regulating hydrological function, some studies have shown that the surface water resources of the HNNR from the Nongjiang River have been completely truncated [48]. Since then, the source of surface water within the HNNR has been limited to the annual rainfall and the flooding of the Heilongjiang and Ussuri rivers over the flood control works during the flood season, which directly affecting the surface hydrological Sanjiang Plain has always been an important base for food production in China. Since the 1950s, it has undergone a number of land developments with a rapid rise in the area of arable land, all dominated by the reclamation of swampy wetlands and including deforestation and the construction of water conservancy projects. The HNNR is surrounded by Qianfeng Farm, Honghe Farm and Yalu River Farm, and the surrounding farms have established extensive irrigation systems. Since 1988, seven drainage trunk canals and a flood control dam have been constructed in the upper reaches of the Nongjiang River basin. This cuts off the Nongjiang River, resulting in the diversion of the Nongjiang River water that previously flowed through the HNNR to discharge directly into the Heilongjiang River [47], further cutting off the upper reaches of the Nongjiang River and depriving it of its fluvial regulating hydrological function, some studies have shown that the surface water resources of the HNNR from the Nongjiang River have been completely truncated [48]. Since then, the source of surface water within the HNNR has been limited to the annual rainfall and the flooding of the Heilongjiang and Ussuri rivers over the flood control works during the flood season, which directly affecting the surface hydrological situation of the reserve. The construction of a dam on the lower reaches of the Nongjiang River within the HNNR in 2008 resulted in more surface water storage in the HNNR, higher water level and some recovery of the hydrological situation. According to the available information and our measured water level data, we know that the water level in the core area of the HNNR in 1993 was 52.0 m in the Wolvlan River, and decreased to 51.4 m in 2002, a decrease of 0.6 m, an average of 66 mm per year, while the water level in the HNNR rose to 51.7 m in 2012, an increase of 0.3 m. According to the spatial distribution and area of the HNNR landscape types we obtained in the previous section, it can be found that the wetland landscape area of the HNNR was 127.03 km 2 in 1993, decreasing to 125.02 km 2 in 2002, while rising to 140.7 km 2 in 2012, with the area of the water also decreasing before rising, thus showing that the change in water level is positively correlated with the change in wetland landscape.
In the mid-1990s, the structure of grain cultivation in Sanjiang Plain began to change significantly, from growing mainly soybeans and wheat to growing mainly soybeans, corn and rice. A rising temperature creates favourable conditions for paddy development, which is developing rapidly [49] with large amounts of groundwater being used to irrigate rice cultivation because of low annual precipitation and the low efficiency of crops in using natural precipitation. Although long-term hydrological data are lacking, available observations show that the groundwater level in the Wolvlan River was 1. In general, the evolution of the HNNR wetland landscape is significantly positively correlated with changes in hydrological conditions and precipitation. Changes in temperature also affect the evolution of the HNNR wetland to some extent, but the correlation is not significant. Precipitation changes ultimately lead to changes in the hydrological situation in the HNNR, which in turn leads to dynamic changes of the HNNR wetland landscape. Therefore, we believe that changes in the hydrological situation caused by natural and human factors are the main factors of landscape changes in the HNNR wetlands, which is consistent with the analysis of previous studies [34].

Deficiencies and Prospects
This paper invokes Landsat series data and uses the Random Forest algorithm for the HNNR wetland landscape classification based on GEE. Landsat series data has the longest time archived image data and high spatial resolution, which can obtain more accurate spatial and temporal changes in the study of long-term changes in wetland landscapes. GEE compensates for the disadvantages of difficult data acquisition, large data volume and inefficient interpretation in large-scale geological application studies, and is not limited by time and space to achieve rapid and batch processing of data [50][51][52], which provides us with great convenience in studying the dynamic changes of the HNNR landscapes from 1985 to 2021. In our study, better classification results were achieved according to the random forest algorithm. In addition, in our literature review, we did not find any analysis of long-term year-by-year landscape dynamics in the HNNR, so we are the only dataset to obtain annual dynamics of long-term HNNR landscapes, which could provide some reference value for subsequent studies.
However, while this study has advantages as a long time series study, it also has some problems. Due to the poor understanding of the conditions in the early study area, the insufficient resolution of remote sensing images and the inadequacy of field survey data, it is difficult to classify the landscape types. There is also a lack of accuracy verification. In addition, although we have generally concluded that hydrological conditions are the main drivers of the evolution of the HNNR wetland based on available data and literature, we cannot be completely sure using correlation analysis to fully determine due to the lack of hydrological data. In our subsequent study, we will try to simulate the hydrological data for the remaining years based on the available hydrological data, so as to analyse the driving factors of the evolution of the HNNR wetland more precisely and reasonably, and provide a more accurate basis for the ecological restoration of the HNNR wetland landscape.

Conclusions
This study analysed remote sensing images from 1985 to 2021 to identify changes in the HNNR wetland landscape and also explored the drivers of change.
The main findings are as follows: (1) Based on Landsat series images provided by GEE, the HNNR landscape classification was conducted using the Random Forest algorithm. The verified average overall accuracies were all greater than 90%, and the Kappa coefficients were all greater than 0.83, achieving better classification results. (2) From 1985 to 2008, the area of wetland landscape in the HNNR showed a general downward trend, the area of meadow showed an upward and then downward trend, the area of woodland decreased and then increased, and the area of arable land continued to increase. The landscape change during this period was mainly from wetland to meadow, and meadow was transformed into woodland and arable land; both wetland and meadow landscapes were becoming increasingly fragmented. From 2008 to 2021, the area of wetland landscape in the HNNR gradually increased and reached a stable state, the area of meadow decreased significantly, the area of woodland and arable land both showed a steady rise, the area of water increased significantly and then stabilised, and the landscape change in this time period was mainly from meadow to wetland. Wetland landscape fragmentation has improved, but meadow landscape fragmentation continues to increase. At the overall level, the HNNR landscape is consistently dominated by wetland landscape, with woodland, arable land and water extending their dominance by encroaching on meadows; the HNNR landscape became complex in shape, favouring wetland stability and increasing landscape heterogeneity. (3) Changes in hydrological conditions due to natural and human factors are the main drivers of the evolution of the HNNR wetland. From 1985 to 2008, human activities such as the construction of water projects and irrigation of farmland, coupled with a declining trend in precipitation, led to a decline in the above-ground water and groundwater levels of the HNNR, resulting in the degradation of the HNNR wetland ecosystems. From 2008 to 2021, due to the dam constructed on the Nongjiang River downstream in the HNNR, the water level in the HNNR rose and both the area of the wetland and water increased, allowing the hydrological situation to recover and the wetland ecosystems to be restored.  = ( ) (100) / 0 < PLAND ≦100. PLAND characterises the total area of patches of a landscape type as a percentage of the total landscape area, quantifying the proportion of abundance of each patch type in the landscape and reflecting the degree of dominance of a landscape. A higher PLAND value indicates that one landscape type is dominant, while a lower PLAND value indicates that the various landscape types that make up the landscape are equally represented. 0 < COHESION <100. COHESION indicates the degree of aggregation of the different patches. Patch cohesion increases as the patch type becomes more clumped or aggregated in its distribution; hence, more physically connected.

1.
Largest Patch Index (LPI) LPI = max n j=1 a ij (100) /A 0 < LPI 100. LPI is the proportion of the entire landscape area occupied by the largest patches of a given patch type. It is an important indicator of the degree of fragmentation of an area at an overall scale. The change in value from largest to smallest is a reflection of the increasing fragmentation of the landscape.

2.
Mean patch size (Area MN ) a ij /N Area MN ≥ 0, without limit. Area MN represents the average condition of the patch area, with smaller values indicating greater fragmentation of the landscape type.

3.
Percentage of landscape (PLAND) PLAND = ( n ∑ j=1 a ij )(100) /A 0 < PLAND 100. PLAND characterises the total area of patches of a landscape type as a percentage of the total landscape area, quantifying the proportion of abundance of each patch type in the landscape and reflecting the degree of dominance of a landscape. A higher PLAND value indicates that one landscape type is dominant, while a lower PLAND value indicates that the various landscape types that make up the landscape are equally represented.

4.
Patch cohesion index (COHESION) Land 2022, 11, 2137 19 of 28 0 < COHESION < 100. COHESION indicates the degree of aggregation of the different patches. Patch COHESION increases as the patch type becomes more clumped or aggregated in its distribution; hence, more physically connected.

5.
Interspersion and juxtaposition index (IJI) 0 < IJI 100. IJI approaches 0 when the corresponding patch type is adjacent to only one other patch type and the number of patch types increases. IJI = 100 when the corresponding patch type is equally adjacent to all other patch types (i.e., maximally interspersed and juxtaposed to other patch types).

6.
Mean patch shape index (SHAPE_MN) 0.25p ij √ a ij /N SHAPE_MN ≥ 0, without limit. SHAPE_MN corrects for the size problem of the perimeter-area ratio index (see previous description) by adjusting for a square standard and, as a result, is the simplest and perhaps most straightforward measure of shape complexity.

7.
Shannon's Diversity Index (SHDI) The SHDI ≥ 0, without limit. The size of the SHDI reflects the number of landscape elements and the change in the proportion of each element The change in the size of the SHDI can characterise the richness and diversity of the various types of landscape structure in a region. The change in value from small to large reflects the increase in the richness of the landscape types and the tendency for the landscape to become more diverse. When a landscape is composed of a single element, the landscape is homogeneous and its diversity index is zero; a landscape composed of more than two elements has the highest diversity when the proportion of landscape types is equal; the diversity of the landscape decreases when the difference in the proportion of each landscape type increases.

8.
Aggregation index (AI) 0 AI 100. Given any P i , AI equals 0 when the patch types are maximally disaggregated (i.e., when there are no like adjacencies); AI increases as the landscape is increasingly aggregated and equals 100 when the landscape consists of a single patch.
The above equations all use the representation in the Fragstats 4.2 software, where the same symbols indicate the same meaning. Where n denotes the number of patches; a ij denotes the area of patch j of landscape category i in terms of number of cells; A denotes the total area of the total landscape; p ij denotes the perimeter of patch j of landscape category i in terms of number of cell surfaces; Z denotes the total number of cells in the landscape; e ik denotes the total length of the edge in landscape between patch types i and k; m denotes the number of patch types (classes) present in the landscape, including the