Ecological Safety Assessment and Analysis of Regional Spatiotemporal Differences Based on Earth Observation Satellite Data in Support of SDGs: The Case of the Huaihe River Basin

: Terrestrial ecosystems provide a variety of beneﬁts for human life and production, and are a key link for achieving sustainable development goals (SDGs). The basin ecosystem is one type of terrestrial ecosystem. Ecological security (ES) assessments are an important component of the overall strategy to achieve regional sustainable development. The Huaihe River Basin (HRB) has the common characteristics of most basins, such as high population density, a rapidly developing economy, and many environmental problems. This study constructed an ES evaluation system by applying a pressure-state-response framework as an assessment method for the sustainable development of basins. Taking the HRB as an example, this study determined the ES status of the region from 2001 to 2019 and analyzed crucial factors for any variation observed by combining remote sensing and climate data, relevant policies, and spatial information technology. The results highlight the importance of reserves and the negative impact of urban expansion on ES. Additionally, the enactment of policies had a positive impact on ES, whereas precipitation had a negative effect on ES in most areas of the HRB. Based on these results, the government should strengthen the protection of forests, grasslands, and wetlands and improve water conservation facilities. This study provides guidance for the subsequent economic development, environmental protection, and the achievements of SDG 15 in the HRB.


Introduction
The 193 member states of the United Nations adopted the 2030 Agenda for Sustainable Development and its 17 Sustainable Development Goals (SDGs) in August 2015 [1]. The 17 SDGs aim to protect the Earth and ensuring peace and prosperity for humankind [2], and provide all countries with a common approach and an agenda to address the serious challenges facing the world today, such as poverty, climate change, and conflicts [3]. SDG 15 is focused on the terrestrial ecosystem, which is the basis of human development, and provides various benefits for human production and life. The aim of SDG 15 is to protect, restore, and promote sustainable use of land; sustainably manage forests; combat desertification; halt land degradation; and prevent biodiversity loss [4].
The basin ecosystem is considered an important component of the terrestrial ecosystem. It is a reasonably closed system with clear boundaries that are involved with material and energy exchange with the external environment [5]. The basin ecosystem provides multiple ecosystem goods and services to residents, including freshwater and oxygen, supplies, purified water, and carbon fixation [6]. With rapid urbanization, extensive economic development, and explosive population growth, basins currently face several environmental problems, such as 15, and combined with literature on the actual HRB situation, this study selected indicators based on the PSR framework. Then, the weights of the indicators were determined using the AHP and an ES indicator assessment system was established. By combining of satellite data and GIS technology, this study revealed the spatiotemporal variation characteristics of ES in the basin. Finally, this study analyzed the spatiotemporal changes in ES and explained the changes using climate and policy data. The results will be useful to determine the causes of eco-environment degradation and provide guidance for implementing SDG 15.

Study Area
The HRB is located in the eastern part of China and covers an area of approximately 2.7 × 10 5 km 2 crossing Jiangsu, Anhui, Shandong, and Henan Provinces [34]. The main stream of the Huaihe River originates from Tongbai Mountain in Henan and flows into the sea at Yangzhou and Yancheng in Jiangsu [35]. There are many tributaries in the middle and upper reaches of the Huaihe River. The Dabie Mountains and hills in the Jiang-Huaihe region are the sources of the southern tributaries. Tributaries are typically characterized by a rushing current and are close to the source [36]. Plains comprise the main landform type in the basin, and the elevation in the western region is higher than the elevation in the eastern region ( Figure 1a). The basin is situated in a climate transition zone, with a humid southern region and a semi-humid northern region [35]. The mean annual precipitation and temperature are approximately 883 mm and 11-16 • C, respectively [37]. The HRB is densely populated with 611 people/km 2 [38] and has experienced a high rate of population growth in recent decades. The utilization rate of land is very high; cropland accounts for the largest proportion of the total area (Figure 1b). The HRB is an important grain, cotton, and oilseed production region in China [39], resulting in high land utilization. The industries mainly consist of coal, power, and light textile industries,  The HRB is densely populated with 611 people/km 2 [38] and has experienced a high rate of population growth in recent decades. The utilization rate of land is very high; cropland accounts for the largest proportion of the total area ( Figure 1b). The HRB is an important grain, cotton, and oilseed production region in China [39], resulting in high land utilization. The industries mainly consist of coal, power, and light textile industries, which have developed rapidly in recent years. Hence, several industrial cities (e.g., Zhengzhou, Xuzhou, and Lianyungang) have emerged.

Workflow
The study selected indicators from ecological conditions and the literature, and determined the weights of the indicators using the AHP and expert judgment. After constructing the indicator evaluation system, the value of each indicator was calculated using satellite, statistical, and geospatial data. Then, the spatial distribution of ES was obtained in each year from 2001 to 2019 by normalizing the results. Finally, combining the precipitation, land use, and policy information, the spatiotemporal variations of ES were analyzed ( Figure 2).

Indicator Calculation
(1) Development intensity refers to the intensity of human construction activities. The overall economy of the HRB has been growing, but economic development does not always occur in a balanced manner [40]; excessive development and urbanization have harmed the environment [41].
Development intensity was calculated as follows: where is the development intensity in the grid, the size of the grids used in the study was 1 km × 1 km, is the built-up area in the grid, is the grid area, is the popu-  (1) Development intensity refers to the intensity of human construction activities. The overall economy of the HRB has been growing, but economic development does not always occur in a balanced manner [40]; excessive development and urbanization have harmed the environment [41].
Development intensity was calculated as follows: where DI is the development intensity in the grid, the size of the grids used in the study was 1 km × 1 km, B is the built-up area in the grid, S G is the grid area, PCC is the population concentration coefficient of the grid, L ave is the average light brightness in the grid, L min is the minimum light brightness in the grid, and L max is the maximum light brightness in the grid. The formula was improved based on the method proposed by Li et al. [42].
(2) The development speed is the area of building land expansion in the region each year. Rapid economic development in the HRB has resulted in damage to the ecoenvironment [41].
The development speed was calculated as follows: where DS is the development speed in the grid, B y is the developed area in the grid in year y, and B y−1 is the developed area in the grid in year y − 1.
(3) The pollution load is the total amount of point-source pollution from industries and urban residential areas and non-point-source pollution from pesticides and fertilizers into the water body in 1 year [43]; water pollution in the HRB is very serious [41].
The pollution load was calculated as follows: where PL is the pollutant source concentration in the grid, E i is the output coefficient of the pollutant for land use type i, C i is the area of land use type i, and PV is the pollutant point-source concentration. The validity of the calculation method and C i were confirmed by Chen et al. [44].
(4) The gradient difference is caused by the topography, as significant regional differences in elevation affects the precipitation, temperature, and other natural factors [45]; they have substantial impacts on human socioeconomic activities.
The gradient difference index was calculated as follows: where GDI is the gradient difference index for the grid, G var is the gradient variance in the grid, and G ave is the average gradient in the grid. The gradient was calculated using the DEM data. (5) The frequency of disasters refers to the number of disasters occurring in 1 year. Drought and flood disasters in the HRB strongly influence the local socioeconomic landscape and environment [46].
The frequency of disasters was calculated as follows: where FOD is the frequency of disasters occurring in the basin and T y is the number of disasters that occurred in year y. (6) Landscape fragmentation, in which habitats or areas of specific land use types are broken down into smaller plots, represents human disturbance of the landscape [47]. The extent of fragmentation can be preliminarily assessed by the number of patches in a particular landscape [47].
Landscape fragmentation was calculated as follows [42]: where LF is the extent of landscape fragmentation in the grid, N i is the number of land patches in the grid, and S G is the grid area.
Remote Sens. 2021, 13, 3942 6 of 21 (7) Ecological vitality refers to the integrity and health of the ecosystem supporting its ability to flourish and develop [48]. The quality of vegetation growth is a sign of ecosystem vitality [49].
Ecosystem vitality was calculated as follows: where EV is the ecosystem vitality of the grid, and N A y is the average NDVI in the grid. (8) Ecological elasticity refers to the ability of an ecosystem to maintain and adjust itself, while resisting various external pressures and disturbances [50].
Ecosystem elasticity was calculated as follows: where EEC is the ecosystem elasticity for the grid, f i is the ecosystem elasticity weight for land use i, C i is the area of land use type i, and S G is the grid area. The validity of the formula was confirmed by Li, Sun, Zhu, and Cao [42], and the ecosystem elasticity weight for each land use type was determined based on expert consultation as follows: forest, 1; grassland, 0.8; cropland, 0.6; building, 0.2; water and wetland, 1; and bare land, 0.4. (9) Reserves are established to conserve biodiversity and prevent habitat fragmentation. The area of reserves in the region indicates the importance of ecological protection by the local government [51].
The reserve area index was calculated as follows: where RI is the reserve area index of the grid, S R is the area of reserve in the grid, and S G is the grid area.

Determination of Weights
This study combined the AHP and the expert's judgment to determine the weights of the indicators ( Table 1). The AHP divides complex problems into a target layer, a criterion layer, and a parameter layer according to the relationship between internal factors, and forms a judgment matrix according to the importance of each factor in the same layer. The maximum eigenvector and random consistency ratio (CR) of the judgment matrix were calculated. The judgment matrix is effective when CR < 0.1. The weight is obtained by normalizing the maximum eigenvector [52].  The CR of the judgement matrix in the pressure layer was 0.0022, the CR of the judgement matrix in the state layer was 0.0462, and the CR of the judgement matrix in the target layer was 0.0079. All weights were obtained when CR values met the requirements.
The indicators were divided into two types: positive (+) and negative (−) depending on their effect on the environment. The ES value was calculated using the following formula: where V ES is the ES value, V P is the pressure value, V S is the state value, and V R is the response value. The range of possible ES values (0-1) can be divided into five equal intervals, such that each interval represents a different ES level ( Table 2).

Data Sources and Processing
The data used in this study were mainly satellite data. They were processed using GIS and RS-related software (e.g., ArcGIS and ENVI); data analysis was carried out using MATLAB software. All satellite data were processed by mosaic, clipping, projection, and format conversion, including nighttime light data, precipitation, surface temperature, digital elevation model (DEM), and MODIS product (i.e., MODIS/Terra Vegetation Indices 16-Day L3 Global 250 m SIN Grid (MOD13Q1 NDVI) and MODIS/Terra + Aqua Land Cover Type Yearly L3 Global 500 m SIN Grid (MCD12Q1 Land Cover Type)). The MOD13Q1 NDVI data was processed by S-G filtering and the maximum value composites (MVC) to eliminate the influence of noise and clouds. Nighttime light data were obtained from the Defense Meteorological Satellite Program Operational Linescan System and the National Polar Partnership's Visible Infrared Imaging Radiometer Suite (NPP/VIIRS). NPP/VIIRS data were derived from the Flint Earth Luminous product, which was provided by the Chinese Academy of Sciences and removed abnormal values. The nighttime light data also underwent supersaturation correction and continuity correction, in accordance with the methods used by Cao et al. [53] and Liang et al. [54]. Considering the difference in the spatial resolution of the data, the regional statistics tool in ArcGIS was used to calculate the mean value of each indicator data in 1 × 1 km grids and assigned the value to grids for subsequent normalization and the weighted calculation. The nighttime light data was used for calculating the development intensity. The MOD13Q1 NDVI data was used to calculate the ecosystem vitality. The DEM data was used for calculating the gradient difference index. The reserve area data was used to calculate the reserve area index. The statistical data was used for calculating the pollution load and frequency of disaster. The MCD12Q1 Land Cover Type was used to calculate the development intensity, development speed, pollution load, landscape fragmentation, and ecosystem elasticity.

Data Sources
This study combined satellite data from different sources and vector data to determine the ES status of the HRB, whereas climate data and policy information were used to analyze variations in the ES (Table 3).

Data Normalization
To eliminate differences (e.g., in terms of magnitude and units) among the indicators, the indicator values were transformed, such that all indicators were dimensionless with values that ranged between 0 and 1 [55]. The transformation formulas used were as follows: where x is the original value; max(x) and min(x) are the maximum and minimum values of indicator x, respectively; X s is the standard value of a positive indicator; and X s is the standard value of a negative indicator.

Land Change Analysis
The land-use transfer matrix is the main method used for investigating conversion among different land-use types and is written mathematically as follows [56]: where P represents the area of interest; n represents the number of land-use types; and i and j represent the land use types at the beginning and end of the study period, respectively.

Trend Analysis
Sen's slope estimator Q i can be used to assess trends in time series. This method was used to obtain the trends of changes in ES and annual average precipitation. Q i was calculated as follows [57]: where n is the number of datapoints in the time series; x j and x k are the data observed in year j and k (j > k), respectively; and Q med is the median of Q i .

Correlation Analysis
The Pearson's correlation coefficient is used to measure the strength of the correlation between two variables. Here, we calculated correlation coefficients using the following formula to assess the correlation between precipitation and the degree of ES [58]: where n is the number of time series datapoints; x i and y i are the precipitation amount and ES value, respectively; x and y are the average precipitation amount and ES value, respectively; and r x,y is the correlation coefficient.

Overall Spatiotemporal Variations in ES of the HRB
To more clearly observe the variations in the ES of the HRB, 1 km × 1 km grids were set up within the HRB area. Each grid contained the result of indicator normalization and weighting, which facilitated intuitive visualization of spatial variation across the region ( Figure 3). Then, proportions of areas with different levels of ES in the HRB were calculated ( Figure 4).      In summary, most areas of the HRB had good sustainable development ability. The better and poor sustainable development ability were associated with reserves and built-up areas, respectively. Shandong had a comparatively worse sustainable development ability, relative to the other provinces. Using 2001 as the starting point, the HRB exhibited its worst sustainable development ability in 2005 and its best sustainable development ability in 2015.

Variations of ES in the Four Provinces
The HRB area was divided into four province areas (Jiangsu, Anhui, Shandong, and Henan) based on provincial boundaries. The proportions of areas with different ES levels in each province were then calculated ( Figure 5). Overall, the sustainable development ability varied spatiotemporally across the four provinces. On a provincial level, Shandong had worse sustainable development ability, while Anhui had better sustainable development ability. The temporal variations in ES level according to the province were consistent with temporal variations for the whole basin.

Trends in ES value of the HRB
The ES ratings for all grids were calculated each year, and Sen's slope estimator was used to assess trends in the ES status across the basin from 2001 to 2019 ( Figure 6). Five trends were distinguished based on the natural points of discontinuity method using ArcGIS: significant decreasing trend, slight decreasing trend, inapparent change, slight increasing trend, and significant increasing trend (Table 4). Then, the proportion of area associated with each trend was calculated for the HRB and Jiangsu, Anhui, Shandong, and Henan provinces.
From 2001 to 2019, 16.33% of the HRB area exhibited an inapparent change trend. These areas were concentrated in the eastern part of the basin. Areas exhibiting significant and slight increasing trends comprised 5.66% and 74.86% of the basin area, respectively. Areas associated with a significant increasing trend were dispersed over the western, southwestern, and southeastern parts of the HRB, whereas areas associated with a slightly increasing trend were distributed across most of the basin. Areas exhibiting significant and slight decreasing trends comprised 0.59% and 2.56% of the total basin area, respectively. When the land use map was overlaid on the trend map, these regions were mainly located in and around built-up areas (Figure 1b).
In Jiangsu, 7.34% of the province's area was associated with a significant increasing trend, more than the proportions in the other provinces. In contrast, 57.23% of the province's area was associated with a slightly increasing trend, the lowest proportion among the four provinces. The proportions of areas associated with inapparent change and decreasing trends were 29.57% and 5.86%, respectively, both higher than the values for the other three provinces. In Anhui, 2.03% of the province's area was associated with decreasing trends, the lowest proportion among the four provinces. Shandong had the lowest proportion of area associated with a significant increasing trend (4.05%). In Henan, areas associated with increasing trends comprised 87.11% of the province's area, the largest proportion of the four provinces. These areas were concentrated in the eastern part of the basin. Areas exhibiting significant and slight increasing trends comprised 5.66% and 74.86% of the basin area, respectively. Areas associated with a significant increasing trend were dispersed over the western, southwestern, and southeastern parts of the HRB, whereas areas associated with a slightly increasing trend were distributed across most of the basin. Areas exhibiting significant and slight decreasing trends comprised 0.59% and 2.56% of the total basin area, respectively. When the land use map was overlaid on the trend map, these regions were mainly located in and around built-up areas (Figure 1b).
In Jiangsu, 7.34% of the province's area was associated with a significant increasing trend, more than the proportions in the other provinces. In contrast, 57.23% of the province's area was associated with a slightly increasing trend, the lowest proportion among the four provinces. The proportions of areas associated with inapparent change and decreasing trends were 29.57% and 5.86%, respectively, both higher than the values for the other three provinces. In Anhui, 2.03% of the province's area was associated with decreasing trends, the lowest proportion among the four provinces. Shandong had the lowest proportion of area associated with a significant increasing trend (4.05%). In Henan, areas associated with increasing trends comprised 87.11% of the province's area, the largest proportion of the four provinces.
In conclusion, the ES of the HRB exhibited an improving trend, namely the sustainable development ability of the HRB was getting better. Based on the proportion of area associated with increasing trends, the ES was mostly likely to improve in Henan. The ES variation trend in Jiangsu exhibited the largest regional variation.

Discussion
The achievement of SDG 15 is essential for human development. There is a lack of models for analyzing the sustainable development of basins. This study established an  In conclusion, the ES of the HRB exhibited an improving trend, namely the sustainable development ability of the HRB was getting better. Based on the proportion of area associated with increasing trends, the ES was mostly likely to improve in Henan. The ES variation trend in Jiangsu exhibited the largest regional variation.

Discussion
The achievement of SDG 15 is essential for human development. There is a lack of models for analyzing the sustainable development of basins. This study established an indicator evaluation system that makes ES analysis more convenient and offers a method to assess sustainable development. By analyzing the spatiotemporal variations in the ES of the HRB, this study contributes to a better understanding of the reasons for eco-environmental degeneration in the HRB, thus facilitating its ecological restoration and achievement of SDG 15.

The Indicator System to Evaluate Sustainable Development of Basins
The HRB has the commonalities of most basins. Industrial and agricultural development and accelerated urbanization have resulted in substantial pressure on the ecoenvironment in the HRB, and many ecological problems, such as water pollution, land degradation, and soil erosion, have been reported [59]. Thus, the HRB is well suited to be chosen as the case.
ES refers to the health and integrity of the ecosystem, and is the cornerstone of regional sustainable development [32]. ES assessment is usually carried out by selecting relevant indicators based on a specific framework. Indicators are selected from nature, society, and the economy (e.g., vegetation quality, land-use situation, and distribution of reserves) that inform ES consistent with the United Nations SDGs. Analysis of ES spatiotemporal variation helps to identify issues that need to be prioritized and the causes of ecological problems; thus, it promotes regional sustainable development. Therefore, this study selected ES as the basis for assessing the eco-environmental status of basins. The PSR model is more convenient and applicable than other models [60,61]; thus, this study chose evaluation indices and constructed an ES assessment system based on the PSR framework.
Due to the regional spatiotemporal heterogeneity of ES status, the importance of analyzing differences in spatial distribution is comparable with the importance of analyzing time-series changes. Combining temporal and spatial analyses allows clear and accurate assessment of the evolving regional ecological situation and the causes of such changes. With advancements in spatial information and RS technologies, the assessment of ES can currently include spatial characteristics, such as the combination of various spatial data to analyze changes in ES [26]. Therefore, based on the RS data and GIS techniques, this study determined the ecological situation and its spatiotemporal variations. The HRB covers a wide area, and considering the difficulty of data acquisition and the suitability of temporal and spatial resolution, the study mainly used MODIS data to calculate indicators. Additionally, MODIS products are widely used in ecological research [62] and they are reliable.
Therefore, a comprehensive ES evaluation system based on carefully selected indices from sufficient investigation of the HRB and SDG 15 clearly and precisely reflects the ecological situation and its spatiotemporal variations, hence facilitating follow-up analyses to elucidate the causes of the identified variations.

Analysis of the Factors Influencing Sustainable Development in the HRB
The results showed that the overall ES of the HRB exhibited an improving trend over time. In 2005, the ES level in most areas declined; most of these areas were concentrated in Shandong. However, the best overall ES was attained in 2015 in the HRB. There were obvious differences in the ES level and its change trend among the four provinces in the basin. Overall, the ES was better in Anhui and Henan.
In the 2001-2019 period, the average annual temperature did not extensively change, and ecological factors, such as streamflow and the frequency of natural disasters, did not appear to affect the temperature [63,64], while precipitation has more of an influence on ecological factors. Thus, Sen's slope estimator was used to calculate the change trends of average annual precipitation over 19 years (Figure 7b). The strength of the correlation between the ES level and the average annual precipitation in the HRB was assessed using Pearson correlation coefficients (Figure 7a). The results showed a negative correlation between the ES level and the average annual precipitation in most areas of the HRB, whereas a positive correlation was only observed in a few areas in the southern part of the basin. From 2001 to 2019, the average annual precipitation in most parts of the HRB exhibited a decreasing trend; an increasing trend was only observed in the eastern and southeastern parts of the basin. When these findings were compared with the change trends in the ES level over the same 19 years (Figure 6a), the change trends of ES levels were revealed to be roughly consistent with the change trends of precipitation. Precipitation in almost all areas in the northern, southwestern, and western parts of the HRB exhibited a decreasing trend. In those areas, the ES level was negatively correlated with precipitation; thus, the ES level exhibited an increasing trend. In eastern and southeastern areas, precipitation exhibited an increasing trend and was negatively correlated with the ES level. Therefore, the ES level mainly exhibited a decreasing trend in those areas. level over the same 19 years (Figure 6a), the change trends of ES levels were revealed to be roughly consistent with the change trends of precipitation. Precipitation in almost all areas in the northern, southwestern, and western parts of the HRB exhibited a decreasing trend. In those areas, the ES level was negatively correlated with precipitation; thus, the ES level exhibited an increasing trend. In eastern and southeastern areas, precipitation exhibited an increasing trend and was negatively correlated with the ES level. Therefore, the ES level mainly exhibited a decreasing trend in those areas.  Figure 8a,b). In the eastern HRB, some areas with a decreased ES rating were distributed in a band. Here, grassland was lost and cropland was gained between 2001 and 2005 (i.e., grassland was transformed into cropland; box 2 in Figure 8a,b). These changes in land use types likely caused the decrease in the ES rating for the HRB in 2005. Additionally, from 2001 to 2005, the precipitation level in the HRB was high and fluctuated greatly, which might have adversely impacted the environment. eastern HRB, some areas with a decreased ES rating were distributed in a band. Here, grassland was lost and cropland was gained between 2001 and 2005 (i.e., grassland was transformed into cropland; box 2 in Figure 8a,b). These changes in land use types likely caused the decrease in the ES rating for the HRB in 2005. Additionally, from 2001 to 2005, the precipitation level in the HRB was high and fluctuated greatly, which might have adversely impacted the environment.  (Figure 8c,d), the areas exhibiting a decreasing trend in ES roughly corresponded to areas that converted to built-up areas. Therefore, urban expansion is presumably a main cause of environmental degradation in the HRB. The largest increase in built-up area was observed between 2015 and 2019 (Figure 8f), which likely  (Figure 8c,d), the areas exhibiting a decreasing trend in ES roughly corresponded to areas that converted to built-up areas. Therefore, urban expansion is presumably a main cause of environmental degradation in the HRB. The largest increase in built-up area was observed between 2015 and 2019 (Figure 8f), which likely led to the increase in areas with a poor ES status. From 2005 to 2015, the annual precipitation did not extensively change (Figure 7c), the increase in built-up area was small, the area of grassland gained was enhanced, and the area of grassland lost was reduced (Figure 8e). During this period, the environmental status was steadily improving. In 2001-2019, areas exhibiting a significant increase in the ES were concentrated in the western and southeastern parts of the HRB. The land use transfer maps for the 2001-2019 period show that the cropland area in the southeast decreased, while grassland in the same area increased (box 4 in Figure 8c,d). The forest area in the west also increased (box 3 in Figure 8c,d) because of aerial seeding afforestation in Henan [65], thus improving the ES of this region. Furthermore, the waterbody and wetland area increased in 2019, compared with 2001, proving that the "Construction and management of key plains and depressions in the HRB" project was effectively implemented (Figure 9). led to the increase in areas with a poor ES status. From 2005 to 2015, the annual precipitation did not extensively change (Figure 7c), the increase in built-up area was small, the area of grassland gained was enhanced, and the area of grassland lost was reduced ( Figure  8e). During this period, the environmental status was steadily improving. In 2001-2019, areas exhibiting a significant increase in the ES were concentrated in the western and southeastern parts of the HRB. The land use transfer maps for the 2001-2019 period show that the cropland area in the southeast decreased, while grassland in the same area increased (box 4 in Figure 8c,d). The forest area in the west also increased (box 3 in Figure  8c,d) because of aerial seeding afforestation in Henan [65], thus improving the ES of this region. Furthermore, the waterbody and wetland area increased in 2019, compared with 2001, proving that the "Construction and management of key plains and depressions in the HRB" project was effectively implemented (Figure 9). At the national level, comprehensive management of the HRB has been underway for a few decades (Figure 9). In 1950, the Government of China promulgated its decision to harness the Huaihe River and clearly defined guidelines for regulating the use of the river; the goals were holding back flood waters and ensuring water conservation. By 2018, the Government of China had issued five major documents; the focus had shifted from flood storage to pollution control, ecological restoration, and sustainable development, consistent with changes in socioeconomic development. Each province issued a different number of documents emphasizing various aspects of Huaihe River management between 2001 and 2018 (Figure 9), which may explain the differences in ES status among the four provinces.

Moving toward Achieving SDGs of the HRB
The results of our analyses showed that overall, the ES status of the HRB is good and exhibits an improving trend. The ES of areas that underwent urban expansion declined; thus, although the overall environmental status of the HRB is good, socioeconomic development has exerted negative impacts in specific areas. Areas with improved ecological conditions in the basin experienced increases in forest, grassland, waterbody, and wetland areas, indicating that the measures taken by the government for ecological restoration had At the national level, comprehensive management of the HRB has been underway for a few decades (Figure 9). In 1950, the Government of China promulgated its decision to harness the Huaihe River and clearly defined guidelines for regulating the use of the river; the goals were holding back flood waters and ensuring water conservation. By 2018, the Government of China had issued five major documents; the focus had shifted from flood storage to pollution control, ecological restoration, and sustainable development, consistent with changes in socioeconomic development. Each province issued a different number of documents emphasizing various aspects of Huaihe River management between 2001 and 2018 (Figure 9), which may explain the differences in ES status among the four provinces.

Moving toward Achieving SDGs of the HRB
The results of our analyses showed that overall, the ES status of the HRB is good and exhibits an improving trend. The ES of areas that underwent urban expansion declined; thus, although the overall environmental status of the HRB is good, socioeconomic development has exerted negative impacts in specific areas. Areas with improved ecological conditions in the basin experienced increases in forest, grassland, waterbody, and wetland areas, indicating that the measures taken by the government for ecological restoration had positive effects. The afforestation activities in Henan Province have been of help in improving the eco-environment. Simultaneously, in most regions, there was a negative correlation between the change in annual precipitation and the ES value, indicating that the increase in precipitation may destroy the eco-environment.
The Government of China and the provinces within the basin have established governance policies focusing on pollution control and protection against disasters (Figure 9).
Although policies related to these areas are important and should be implemented, government departments should also carefully consider the rational utilization and saving of precipitation and river dredging. It is important to protect and restore land types with high ecosystem value (e.g., forests, grasslands, and wetlands) [66,67], thus improving overall ecosystem function in the HRB.

Conclusions
This study established a universal ES evaluation system based on the ecological characteristics of the HRB using the PSR model. The indicators are related to the SDG 15, to provide a method to evaluate the sustainable development of basins.
Combining GIS and RS technology, the ES of the HRB was calculated and factors affecting the sustainable development were determined. The results reveal that the spatiotemporal regional distribution of ES differed. Regions near the reserves and buildings are the areas with the best and worst ecological conditions, respectively. The annual precipitation had a negative impact on ES in most areas of the HRB. The addition of grasslands, forests, and wetlands would improve the eco-environment.
The results of these analyses can be used to guide future ecological restoration in the HRB, contributing to the realization of SDG 15. The government needs to focus on water conservancy and vegetation restoration. Some measures, such as aerial afforestation, should be implemented.