Spatiotemporal Pattern, Evolutionary Trend, and Driving Forces Analysis of Ecological Quality in the Irtysh River Basin (2000–2020)

Abstract

Considering climate change and increasing human impact, ecological quality and its assessment have also received increasing attention. Taking the Irtysh River Basin as an example, we utilize multi-period MODIS composite imagery to obtain five factors (greenness, humidity, heat, dryness, and salinity) to construct the model for the amended RSEI (ARSEI) based on the Google Earth Engine platform. We used the Otsu algorithm to generate dynamic thresholds to improve the accuracy of ARSEI results, performed spatiotemporal pattern and evolutionary trend analysis on the results, and explored the influencing factors of ecological quality. Results indicate that: (1) The ARSEI demonstrates a correlation exceeding 0.88 with each indicator, offering an efficient approach to characterizing ecological quality. The ecological quality of the Irtysh River Basin exhibits significant spatial heterogeneity, demonstrating a gradual enhancement from south to north. (2) To evaluate the ecological quality of the Irtysh River Basin, the ARSEI was utilized, exposing a stable condition with slight fluctuations. In the current research context, the ecological quality of the Irtysh River Basin watershed area is projected to continuously enhance in the future. This is due to the constant ecological protection and management initiatives carried out by countries within the basin. (3) Precipitation, soil pH, elevation, and human population are the main factors influencing ecological quality. Due to the spatial heterogeneity, the driving factors for different ecological quality classes vary. Overall, the ARSEI is an effective method for ecological quality assessment, and the research findings can provide references for watershed ecological environment protection, management, and sustainable development.

1. Introduction

The impact of human society’s rapid development on the Earth’s ecological environment is profound. For example, the global surface temperature has already risen by 1 °C, and further acceleration of human activities is expected to continue this trend, potentially exceeding 1.5 °C within the next two decades [1]. The rise in temperature has accelerated glacial melting, leading to an increase in sea levels. Additionally, this temperature rise has caused extreme weather events such as droughts, floods, heavy rainfall, and tropical cyclones, causing significant damage to biodiversity and vulnerable ecosystems [2]. Watershed ecology is a rich field for studying the spatiotemporal patterns, mechanisms, and regulations of the earth’s surface elements [3]. However, it is also a fragile ecosystem. The rise in global temperature has had a significant impact on watershed ecosystems, leading to the depletion of lakes, disrupted river flows, reduced biological diversity, diminished ecosystem function, and the disappearance of wetlands [4,5,6]. Therefore, it is important to investigate and assess the ecological quality (EQ) of watersheds, especially those located in harsh environments such as arid and semi-arid regions or cross-border areas, which are typically the focus of attention for countries related to the basin. Taking the Irtysh River as an example, its watershed ecology consists of a variety of ecosystems, such as arid, mountainous, piedmont plains, valley plains, and artificial oases, which is comprehensive and typical [7]. Therefore, more attention should be paid to the status and trends of the ecological environment of the Irtysh River Basin (IRB) under climate change. Existing research on the IRB has primarily focused on water resources [8,9,10], hydrology [11,12], climate change [13,14], and ecological species [15,16,17,18]. Nevertheless, there are relatively few works on evaluating EQ at the watershed scale [19,20]. It is important to note that a part of the watershed may not represent the ecological characteristics of the entire watershed. Additionally, characterizing ecological changes across the watershed over long periods of time can be challenging. Therefore, to fully assess the EQ and multi-year dynamic change characteristics of the watershed, it is necessary to consider the integrity of the geographical region on a large scale and over long time series.

The complexity and comprehensiveness of the environmental domain require the use of a combination of indicators for characterization. For example, the US Environmental Protection Agency used the Environmental Quality Index (EQI) to describe the overall environmental health status of all counties in the United States, and its domain-specific EQI loadings indicate which of the environmental domains account for most of the variability in the EQI environment [21]. The Forest Affinity Index, a newly proposed ecological index that is sensitive to qualitative changes in the specific composition of assemblages, was found to be the only parameter potentially useful for assessing the ecological complexity of poplar stands [22]. The current Ecological Index of the Chinese ecological environment standards can reflect the comprehensiveness, integrity, and hierarchy of ecological environment assessment, the component indices of which are difficult to obtain, and the visualization of this index is difficult [23].

With the rapid development of Earth observation remote sensing satellite missions, remote sensing (RS) has become increasingly important in the study of dynamic changes in the environment, mainly due to its advantages of surface information acquisition, large coverage, and high temporal continuity. In the past decades, various RS image band operation-based indices, such as the Normalized Difference Vegetation Index (NDVI), the Normalized Difference Water Index, and others, have been widely used in environmental research due to their direct, simple, and general characteristics [24,25]. It is difficult to comprehensively assess ecological conditions with a single indicator, so the Remote Sensing Ecological Index (RSEI) was proposed in 2013, which uses principal component analysis (PCA) to combine greenness, heat, dryness, and moisture to characterize EQ [26]. Benefiting from its obtainability without the need for artificial weighting or thresholding and its computability, the RSEI has been widely used as a comprehensive and detailed method to measure EQ [27,28,29].

In addition, to mitigate the limitations caused by the heterogeneity of different ecosystems, various scholars have modified the RSEI according to the characteristics of the study area. For example, aerosol optical depth as an indicator reflecting urban air quality was added to the RSEI, which can comprehensively assess urban ecological quality, and this improved method is mainly applicable to urban EQ assessment [30]. The modified RSEI combines the first to third principal components to express EQ, which ignores the characteristic representativeness of each principal component in detecting temporal changes, e.g., PC1 represents common information and PC3 represents changes in information [31]. The Comprehensive Salinity Index (CSI) and Water Network Density are included in the RSEI to assess the impact of soil salinity in arid regions and the role of rivers and lakes on the ecology, but the overall correlation of the ecological index may be reduced because the Water Network Density is a local characteristic of a watershed [32]. The moving window-based RSEI has been proposed to take into account the spatial distance between the object and the open pit mine, which provides more accurate information on spatial changes in the ecological environment caused by open pit mines than the RSEI [33]. The RSEI needs to mitigate the impact of significant water bodies on the moisture fraction representing soil and vegetation moisture to authentically reflect the ecological conditions of the study area [34]. Selecting an appropriate method to remove water is therefore crucial. The water index method, known for its high accuracy and low computational cost, is widely used for this purpose. The Modified Normalized Difference Water Index (MNDWI), widely employed due to its speed, simplicity, and precision in extracting water bodies, is extensively utilized [35]. When calculating the RSEI at different time intervals, the use of dynamic thresholds is essential to avoid overestimation or underestimation of water bodies. The Otsu thresholding algorithm, based on maximizing the inter-class variance and minimizing the intra-class variance in the image histogram, has been shown to be the most efficient method for obtaining thresholds in surface water extraction [36,37,38].

Google Earth Engine (GEE) provides PB-level RS data, massive cloud computing resources, and extensive processing and analysis algorithms to process large, long-term RS images [39]. It replaces the traditional desktop processing platforms that require significant time and hardware resources for pre-processing such as acquisition, stitching, atmospheric correction, and cloud and shadow removal of massive remote sensing data [40,41]. Researchers have combined GEE with the RSEI for EQ assessment and monitoring in regions such as the Jianghan Plain (Yi et al., 2023), the Loess Plateau (Gong et al., 2023), northern Anhui, China (Wang et al., 2022), the Erhai Lake Basin in Yunnan Province, China (Xiong et al., 2021), and the Yellow River Basin (Yang et al., 2022) [42,43,44,45,46]. This also demonstrates the convenience and efficiency of using the GEE platform for the RSEI assessment on a large scale and over extended time series of imagery. In previous research, the application of the RSEI has typically involved performing analyses based on the synthesis of images from two or more discrete time periods, but it cannot appropriately reflect the continuous long-term EQ changes [26,43,44,47]. Detailed changes in EQ over time scales can be better captured by analyzing time series trajectories [48]. The GEE platform provides appropriate time series data analysis methods, such as change vector analysis [49], continuous change detection and classification [50], and the LandTrendr algorithm for monitoring forest disturbances [51]. LandTrendr is a spectral–temporal segmentation algorithm suitable for detecting annual changes in the time series of medium-resolution satellite imagery. It can identify abrupt, short-term changes in target objects within images and distinguish long-term ecological recovery [51]. However, there has been limited research using LandTrendr to analyze EQ changes in continuous time series. As the RSEI is derived from synthetic imagery calculated for the vegetation growing season within a year, LandTrendr is particularly suited to analyzing continuous temporal changes in EQ based on RSEI assessments [52].

EQ is constrained by a combination of various natural and socio-ecological factors. Therefore, analyzing the impact of different factors on its changes is important for the evaluation, protection, and restoration of the ecological environment. The regression tree model [53], the regression PCA model [54], correlation analysis [55], and the Geographic Detector [56] are commonly used to study environmental drivers. However, the regression tree model, regression PCA model, and correlation analysis have varying degrees of uncertainty [57]. The Geographic Detector cannot analyze continuous data as it requires discrete data inputs [58]. The Patch-generating Land Use Simulation (PLUS) model, based on cellular automata, incorporates a new multi-type seed growth mechanism that can better simulate patch-level changes in multi-type land use. It has a strong advantage in studying the space–time dynamic evolution of simulated spatial complex systems. The land expansion analysis strategy (LEAS) module of the PLUS uses the random forest algorithm to mine the expansion and driving factors of different land use types one by one to obtain the evolution of the probability of different land uses and the contribution of driving factors to the expansion of different land uses in each period [59]. This method combines the advantages of the transformation analysis strategy (TAS) and the pattern analysis strategy (PAS). It avoids the analysis of transformation type, which increases exponentially with the number of categories, and retains the ability of the model to analyze the mechanism of land use change in each period, which has better interpretability. The LEAS supports the input of driving factor data with different resolutions and can realize automatic coordinate alignment, which greatly simplifies the operation and saves time in data processing.

Based on the above analysis, this paper aims to use Moderate Resolution Imaging Spectroradiometer (MODIS) synthetic images on the GEE platform to construct an amended RSEI (ARSEI). (1) Under the premise of considering long-term changes in the time series, dynamic thresholds are obtained by the Otsu algorithm to improve the accuracy of the ARSEI. The addition of the salinity index improves the accuracy of the ARSEI in describing EQ. (2) Taking the IRB as an example, the ARSEI is used to quantify its ecological environment over the past 20 years. (3) The paper reveals the inter-annual dynamics of regional EQ, emphasizing the value of regional EQ at a fine temporal scale. (4) The relative contributions of climate change and human activities to EQ are determined.

2. Materials and Methods

2.1. Study Area

The Irtysh River is an international river that originates in the Altay region of the Xinjiang Uygur Autonomous Region in China and flows through Kazakhstan and Russia (Figure 1). It spans 15° latitude, with a total length of 4248 kilometers. This basin boasts a vast drainage area and displays significant ecological diversity [60]. The IRB’s upper region exhibits a temperate semi-desert climate, notable for its relatively abundant rainfall, lower evaporation rates, and humid environment [61]. In contrast, the basin’s middle and lower sections, situated within the West Siberian Plain, are subject to a distinct continental climate, featuring extreme temperatures in winter and summer [62]. The prevalent westerly and southerly winds across the plain contribute to the dryness of the basin’s midstream areas, especially between Pavlodar and the Shagan River, where the annual rainfall is approximately 250 millimeters [63]. In the IRB, a systematic south-to-north reduction in soil pH is observed, marking a transition from weakly alkaline to weakly acidic conditions. The upper basin, particularly adjacent to the western or north-western slopes of the Altay Mountains in Kazakhstan, as well as the lower basin, are extensively forested with a mix of coniferous and deciduous trees. Predominant among these species are Pinus sibirica, Larix sibirica Ledeb., Picea obovata, and Abies sibirica [64,65,66]. Contrasting this, the middle reaches are characterized by diverse land use types, with farmland, grassland, and bare land being most prevalent. The IRB can connect the Arctic, Central Asia, and South Asia through roads, railways, and even pipelines to achieve vertical penetration of the Eurasian continent. It has obvious geographical advantages and provides a possible path for the integration of Asia and Europe. At the same time, as an important area of the Belt and Road, it has diverse ecological types and rich animal and plant resources. Monitoring and evaluating the EQ of this watershed can help cooperation in the field of ecological civilization under the framework of the Belt and Road, to the advantage of the citizens of all countries involved in the joint building of the Belt and Road [67].

2.2. Datasets

The data and its information used in this study can be found in

Table 1. Summary of the datasets.

Dataset	Time	Data Type	Spatial Resolution	Resample Method	Data Source
SPWI/CSI/NDBSI	2000–2020	Raster	500 m	-	MOD09A1 Dataset (https://lpdaac.usgs.gov/products/mod09a1v061/, accessed on 1 June 2023.)
LST	2000–2020	Raster	1000 m	bicubic	MOD11A2 Dataset (https://lpdaac.usgs.gov/products/mod11a2v061/, accessed on 1 June 2023.)
EVI	2000–2020	Raster	250 m	bicubic	MOD13Q1 Dataset (https://lpdaac.usgs.gov/products/mod13q1v061, accessed on 1 June 2023.)
DEM	-	Raster	90 m	bilinear	MERIT DEM (http://hydro.iis.u-tokyo.ac.jp/, accessed on 1 June 2023.)
Soil data	2020	Raster	250 m	bicubic	SoilGrids250m 2.0 data (https://data.isric.org, accessed on 1 June 2023.)
Total precipitation	2000–2020	Raster	0.1°	bicubic	Total precipitation (https://cds.climate.copernicus.eu, accessed on 1 June 2023.)
Population density	2000–2020	Raster	1000 m	Nearest neighbor	GPW v4 (https://sedac.ciesin.columbia.edu/, accessed on 1 June 2023.)
ARSEI	2000–2020	Raster	500 m	-	Data products for this article

Note: All data have been harmonized to 500 m resolution using the resampling methods listed in the table.

.

Remote sensing data: The MODIS images come from the GEE platform, including MOD09A1, MOD11A2, and MOD13Q1.

Natural factor data: The DEM data is a Multi-Error-Removed Improved-Terrain DEM, which is used by ArcGIS to generate slope and aspect [68]. The total precipitation per year data comes from the ERA5 MONTHLY data of the European Centre for Medium-Range Weather Forecasts [69]. The soil data comes from the soil information provided by the International Soil Reference and Information Center (ISRIC), including bulk density (bdod), cation exchange capacity at pH 7 (cec), coarse fragments (cfvo), clay (clay), total nitrogen (nitrogen), pH in H2O (phh2o), sand (sand), and silt (silt) [70].

Human factor data: The population density data comes from The Gridded Population of the World, Version 4 (GPWv4) of the United Nations World Population Prospects data [71]. The value in each cell represents the population density of that cell.

2.3. Methodology

This study proposes a research framework to explore the interaction between the natural environment and human activities on EQ (Figure 2). Human activities accelerate the degradation of EQ. At the same time, the natural environment, including climate, soil, and topography, is an important component of EQ. Using the IRB as an example, this study uses the ARSEI to assess and analyze EQ at a large spatial scale. The introduction of the LandTrendr and PLUS models allows for the analysis of EQ changes in both space and time at the pixel level, investigating the driving factors and their role in influencing EQ variations.

2.3.1. Construction of the ARSEI

In the watershed ecosystem, there is a close relationship between the different ecological factors, which promote and restrict each other and act together on the ecology of the watershed. How to determine the correlation between each ecological factor and the weight of its effect on the ecosystem becomes the key point when assessing EQ with multiple ecological factors. PCA can automatically and objectively determine its weight according to the contribution of each indicator to each principal component, thereby effectively reducing the range of data while retaining the maximum amount of information. To comprehensively monitor and evaluate the EQ of the watershed, this paper also selects PCA as the coupling method for five ecological factors, such as greenness, humidity, heat, dryness, and salinity, to construct the ARSEI (

Table 2. Summary of calculation methods for ecological factors.

Ecological Factor	Ecological Indicator	Calculation Method
Greenness	EVI	MODIS Product (MOD13Q1)
Humidity	SPWI	$\mathrm{S}\mathrm{P}\mathrm{W}\mathrm{I}=\frac{{\mathsf{\rho }}_{\mathrm{N}\mathrm{I}\mathrm{R}}-{\mathsf{\rho }}_{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}2}+{\mathsf{\rho }}_{\mathrm{B}}}{{\mathsf{\rho }}_{\mathrm{N}\mathrm{I}\mathrm{R}}+{\mathsf{\rho }}_{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}2}+{\mathsf{\rho }}_{\mathrm{B}}}$
Heat	LST	MODIS Product (MOD11A2)
Dryness	NDBSI	$\mathrm{N}\mathrm{D}\mathrm{B}\mathrm{S}\mathrm{I}=\left(\mathrm{S}\mathrm{I}+\mathrm{I}\mathrm{B}\mathrm{I}\right)÷2$ $\mathrm{S}\mathrm{I}=\frac{\left({\mathsf{\rho }}_{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1}+{\mathsf{\rho }}_{\mathrm{R}}\right)-\left({\mathsf{\rho }}_{\mathrm{N}\mathrm{I}\mathrm{R}}+{\mathsf{\rho }}_{\mathrm{B}}\right)}{\left({\mathsf{\rho }}_{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1}+{\mathsf{\rho }}_{\mathrm{R}}\right)+\left({\mathsf{\rho }}_{\mathrm{N}\mathrm{I}\mathrm{R}}+{\mathsf{\rho }}_{\mathrm{B}}\right)}$ $\mathrm{I}\mathrm{B}\mathrm{I}=\frac{2{\rho }_{SWIR1}/({\rho }_{SWIR1}+{\rho }_{NIR})-\left[\right({\rho }_{NIR}/({\rho }_{NIR}+{\rho }_R)+{\rho }_G/({\rho }_G+{\rho }_{SWIR1})]}{2{\rho }_{SWIR1}/({\rho }_{SWIR1}+{\rho }_{NIR})+\left[\right({\rho }_{NIR}/({\rho }_{NIR}+{\rho }_R)+{\rho }_G/({\rho }_G+{\rho }_{SWIR1})]}$
Salinity	CSI	$\mathrm{C}\mathrm{S}\mathrm{I}=\left(SI\_T+NDSI+SI3\right)÷3$ $\mathrm{S}\mathrm{I}\_\mathrm{T}=\left({\mathsf{\rho }}_{\mathrm{R}}/{\mathsf{\rho }}_{\mathrm{N}\mathrm{I}\mathrm{R}}\right)\times 100$ $NDSI=\left({\mathsf{\rho }}_R-{\mathsf{\rho }}_{NIR}\right)/\left({\mathsf{\rho }}_R+{\mathsf{\rho }}_{NIR}\right)$ $\mathrm{S}\mathrm{I}3=\sqrt{{\mathsf{\rho }}_G^2+{\mathsf{\rho }}_R^2}$

Note: ρ represents the reflectance of the remote sensing band; B is the blue band; R is the red band; G is the green band; NIR is the near infrared band; SWIR1 is the short-wave infrared 1 band; and SWIR2 is the short-wave infrared 2 band.

).

(1) Greenness: The spatial distribution pattern of vegetation in the study area is a high vegetation cover area with dense forest in the north and a low vegetation cover area with grassland and cultivated land in the south. Thus, EVI is chosen as the greenness factor to avoid the weaknesses of NDVI, which is quickly saturated in heavily vegetated areas and easily influenced by the soil background in sparsely vegetated areas [72].

(2) Humidity: This paper uses the Surface Potential Water Abundance Index (SPWI) as the humidity factor. It uses the normalized index method to calculate the NIR and SWIR-2 bands while adding the blue band to adjust the index value, which can correctly reflect the abundance of surface water resources compared to wetness (WET) [73].

(3) Heat: Land surface temperature (LST) is not only a measure of climate change and land desertification in ecological terms but also one of the parameters characterizing surface aridity.

(4) Dryness: The Normalized Differential Build-up and Bare Soil Index (NDBSI) expresses the dryness factor. It is obtained by averaging the Soil Index (SI) and the Index-based Build-up Index (IBI) [26].

(5) Salinity: Instead of a single salinity index, the Comprehensive Salinity Index (CSI) is used in this paper. It more accurately reflects the ecological impact of salinity on a large scale (soil heterogeneity and land use type differences are obvious). The CSI is based on the idea of integrated learning and integrates three universal salinity indices: the Salinity Index (SI-T), the Normalized Differential Salinity Index (NDSI), and the Salinity Index 3 (SI3) to improve the stability and reliability of the detection results [32].

This study implemented the dynamic thresholding MNDWI method for water extraction by introducing the Otsu algorithm and generating a water mask. The use of precise water masks for each period improved the accuracy of the ARSEI results [34,35,36]. To avoid weight imbalances due to dimensional inconsistencies, it is imperative that each ecological factor be traversed and normalized to the range [0, 1] prior to PCA. Then, according to the “$3\sigma$” principle of a normal distribution, the interval $\left(\mu -3\sigma ,\mu +3\sigma \right)$ can effectively eliminate outliers among pixels. Then, the processed ecological factors are used as bands to synthesize multi-band images for PCA. The formula for PCA is as follows:

$$\mathrm{P}{\mathrm{C}}_{\mathrm{i}}=\mathrm{f}\left(\mathrm{E}\mathrm{V}\mathrm{I},\mathrm{S}\mathrm{P}\mathrm{W}\mathrm{I},\mathrm{L}\mathrm{S}\mathrm{T},\mathrm{N}\mathrm{D}\mathrm{B}\mathrm{S}\mathrm{I},\mathrm{C}\mathrm{S}\mathrm{I}\right),$$

(1)

$$\mathrm{A}\mathrm{R}\mathrm{S}\mathrm{E}\mathrm{I}=1-\mathrm{P}{\mathrm{C}}_1,$$

(2)

where $\mathrm{f}\left(\cdot \right)$ represents the PCA operation; $\mathrm{P}{\mathrm{C}}_{\mathrm{i}}$ is the obtained principal component; $\mathrm{P}{\mathrm{C}}_1$ is the first principal component; and the ARSEI is the ultimate result of PC1, with its values representing the assessment outcomes of EQ. High values indicate better EQ, while low values suggest poorer EQ.

To support the measurement and comparison of ecological quality, the min-max normalization method adjusts the ARSEI to [0,1]. When the ARSEI is close to 1, it reflects a higher level of ecological quality, and when it is further away, the quality is lower. The formula is as follows:

$${\mathrm{X}}^{\mathrm{’}}=\left(\mathrm{X}-{\mathrm{X}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\right)/\left({\mathrm{X}}_{\mathrm{m}\mathrm{a}\mathrm{x}}-{\mathrm{X}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\right),$$

(3)

where ${\mathrm{X}}^{\prime }$ indicates the normalized ARSEI, ${\mathrm{X}}_{\mathrm{m}\mathrm{a}\mathrm{x}}$ is the maximum value of the ARSEI, and ${\mathrm{X}}_{\mathrm{m}\mathrm{i}\mathrm{n}}$ is the minimum value of the ARSEI.

2.3.2. Time Series Stability Analysis

The coefficient of variation (CV) is the ratio of the standard deviation to the mean, which is a measure of the variability relative to the mean. It is free of units of measurement and is often used to measure the dispersion of data and to compare them with each other [74,75,76]. In this paper, it is used to reflect the degree of dispersion of the ARSEI and to assess the inter-annual time series stability of the ARSEI. The larger the CV value, the more discrete the distribution of the ARSEI and the greater the inter-annual variation. Otherwise, the concentrated distribution of the ARSEI and the smaller inter-annual fluctuation are the results. The calculation is shown below:

$${\mathrm{C}\mathrm{V}}_{\mathrm{A}\mathrm{R}\mathrm{S}\mathrm{E}\mathrm{I}}={\mathsf{\sigma }}_{\mathrm{A}\mathrm{R}\mathrm{S}\mathrm{E}\mathrm{I}}/\overline{\mathrm{A}\mathrm{R}\mathrm{S}\mathrm{E}\mathrm{I}},$$

(4)

where ${\mathsf{\sigma }}_{\mathrm{A}\mathrm{R}\mathrm{S}\mathrm{E}\mathrm{I}}$ is the standard deviation of the ARSEI and $\overline{\mathrm{A}\mathrm{R}\mathrm{S}\mathrm{E}\mathrm{I}}$ is the mean of the ARSEI.

2.3.3. Spatiotemporal Change Detection Algorithm of the ARSEI

LandTrendr, utilizing an annual time series decomposition algorithm, identifies vegetation trends through spectral pixel trajectories (Figure 3). In this study, LandTrendr was used to investigate the degradation (loss) and improvement (gain) of EQ on an inter-annual scale. LandTrendr is now operational in Google Earth Engine (LT-GEE), significantly improving processing efficiency [51]. The parameters of LT-GEE can be found in the supplementary material (Table S1).

2.3.4. Future Trend Analysis

This paper uses Theil-Sen median trend analysis coupled with a Mann-Kendall test to analyze future ARSEI trends in the IRB. It has been successfully applied in various studies related to hydrological and meteorological trend changes [77,78]. Theil-Sen median has the advantages of not being affected by missing data, being insensitive to measurement errors and outlier data, and having high computational efficiency [79,80]. It was calculated as follows:

$${\mathrm{S}}_{\mathrm{A}\mathrm{R}\mathrm{S}\mathrm{E}\mathrm{I}}=\mathrm{M}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{n}\left(\left({\mathrm{X}}_{\mathrm{j}}-{\mathrm{X}}_{\mathrm{i}}\right)/\left(\mathrm{j}-\mathrm{i}\right)\right),\forall \mathrm{j}>\mathrm{i},$$

(5)

where ${\mathrm{S}}_{\mathrm{A}\mathrm{R}\mathrm{S}\mathrm{E}\mathrm{I}}$ is the Theil-Sen median slope; Median is the median function; the ARSEI inclines when ${\mathrm{S}}_{\mathrm{A}\mathrm{R}\mathrm{S}\mathrm{E}\mathrm{I}}>0$, and the ARSEI declines when ${\mathrm{S}}_{\mathrm{A}\mathrm{R}\mathrm{S}\mathrm{E}\mathrm{I}}<0$.

Mann-Kendall is a non-parametric test. This parametric test method is more appropriate for ordinal variables as it does not require samples to be drawn from a specified distribution and is less susceptible to outliers. The statistical significance of the Theil-Sen trend analysis is defined by the Z-value of MK. If the magnitude of Z is greater than 1.65 and 1.96, it means that the trend has passed the significance test with 90% and 95% confidence, respectively (

Table 3. Classification table of the Sen-MK trend and Hurst index changes.

SARSEI	Z	Trend	H	Types of EQ Changes
SARSEI > 0	1.96 < Z	Significant improvement	H > 0.5	Persistent and significant improvement
	1.65 < Z ≤ 1.96	Slight improvement		Persistent and slight improvement
SARSEI < 0	1.96 < Z	Significant deterioration	H < 0.5	Persistent and significant degradation → Persistent and significant improvement
	1.65 < Z ≤ 1.96	Slight deterioration		Persistent and slight degradation → Persistent and slight improvement
SARSEI = 0	0 < Z ≤ 1.65	Stable without change	-	Persistent and stable
SARSEI > 0	1.65 < Z ≤ 1.96	Slight improvement	H < 0.5	Persistent and slight improvement → Persistent and slight degradation
	1.96 < Z	Significant improvement		Persistent and significant improvement → Persistent and significant degradation
SARSEI < 0	1.65 < Z ≤ 1.96	Slight deterioration	H > 0.5	Persistent and slight degradation
	1.96 < Z	Significant deterioration		Persistent and significant degradation

). It is calculated as follows:

$$\mathrm{s}\mathrm{g}\mathrm{n}\left({\mathrm{x}}_{\mathrm{j}}-{\mathrm{x}}_{\mathrm{i}}\right)=\left\{\begin{array}{cc}1& \left({\mathrm{x}}_{\mathrm{j}}-{\mathrm{x}}_{\mathrm{i}}>0\right)\\ 0& \left({\mathrm{x}}_{\mathrm{j}}-{\mathrm{x}}_{\mathrm{i}}=0\right)\\ -1& \left({\mathrm{x}}_{\mathrm{j}}-{\mathrm{x}}_{\mathrm{i}}<0\right)\end{array}\right.,$$

(6)

$$\mathrm{S}={\sum }_{\mathrm{i}=1}^{\mathrm{n}-1}{\sum }_{\mathrm{j}=\mathrm{i}+1}^{\mathrm{n}}\mathrm{s}\mathrm{g}\mathrm{n}\left({\mathrm{x}}_{\mathrm{j}}-{\mathrm{x}}_{\mathrm{i}}\right),$$

(7)

$$\mathrm{V}\mathrm{A}\mathrm{R}\left(\mathrm{S}\right)=\mathrm{n}\left(\mathrm{n}-1\right)\left(2\mathrm{n}+5\right)/18,$$

(8)

$$\mathrm{Z}=\left\{\begin{array}{cc}\left(\mathrm{S}-1\right)/\sqrt{\mathrm{V}\mathrm{A}\mathrm{R}\left(\mathrm{S}\right)}& \left(\mathrm{S}>0\right)\\ 0& \left(\mathrm{S}=0\right)\\ \left(\mathrm{S}+1\right)/\sqrt{\mathrm{V}\mathrm{A}\mathrm{R}\left(\mathrm{S}\right)}& \left(\mathrm{S}<0\right)\end{array}\right.,$$

(9)

In the formula, $\mathrm{s}\mathrm{g}\mathrm{n}$ is a sign function, n is the amount of data in the sequence, and S is the test statistic when $\left|\mathrm{z}\right|<{\mathsf{\mu }}_{1-\mathsf{\alpha }/2}$ is the significant change in the ARSEI studied at a given significant level α.

The Hurst index is a numerical measure of the persistence of time series data used in hydrology, economics, and climatology [81,82]. In this paper, the Hurst index was calculated using the rescaled polar difference R/S analysis to analyze the persistence of future change trends of the ARSEI in the IRB, which can better express the inter-annual variability characteristics. For example, 0 < H < 0.5 represents the inverse persistence of the time series data, i.e., the future change trend is opposite to the past trend; H = 0.5 represents the randomness of the time series data, i.e., the future change trend is not correlated with the past; and 0.5 < H < 1 shows the time series data has persistence, i.e., the future trend is consistent with the past trend. The Theil-Sen median trend analysis and the Mann-Kendall test (Sen-MK) and Hurst indices were combined to analyze the persistence of the ARSEI change trends (Table 3).

2.3.5. Analysis of Driving Factors

The analysis of the drivers of EQ change is based on the PLUS model [59]. It integrates a rule mining framework based on the LEAS and a CA model based on multi-type random seeds (CARS). By analyzing different types of land use expansion and driving factors, the evolution of the probability of different types of land use and the contribution of driving factors to different types of land use expansion during this period can be obtained, and land use patches and land use scenario simulations can be predicted. In this paper, the LEAS explores the contribution of different driving factors to the hierarchical expansion of EQ and evaluates the role of each driving factor in the change in EQ in the IRB. The formula is as follows:

$${\mathrm{P}}_{\mathrm{i},\mathrm{k}\left(\mathrm{x}\right)}^{\mathrm{d}}={\sum }_{\mathrm{n}=1}^{\mathrm{M}}\mathrm{I}\left({\mathrm{h}}_{\mathrm{n}}\left(\mathrm{x}\right)=\mathrm{d}\right)/\mathrm{M},$$

(10)

where $P_{i,k\left(x\right)}^d$ is the development probability of land type k in unit i; x is a vector composed of multiple driving factors; I is the indicator function of the decision tree; $h_n\left(x\right)$ is the prediction type of the nth decision tree of the vector x; and M is the total number of decision trees. When d is 1, it means that it converted other land use types to type k, and when d is 0, it means other transitions.

3. Results

3.1. Performance Evaluation of the ARSEI in the IRB

3.1.1. Advantages and Applicability of the ARSEI

To visually demonstrate the suitability of the proposed ARSEI compared to the conventional RSEI, we first select several typical areas in the IRB for the experimental analysis in Figure 4. From the results, it can be easily seen that: (1) compared to the RSEI, the ARSEI can more accurately reflect the real surface conditions near the water; and (2) compared to the RSEI, the ARSEI can better reflect the real EQ.

The greenness index in the RSEI is represented using the NDVI, which tends to saturate in areas with high vegetation cover and is susceptible to soil influence in regions with low vegetation cover. The ARSEI addresses these issues by employing the EVI and introducing the CSI as a salinity index to better reflect conditions in the arid regions of the IRB in Central Asia. By comparing the ARSEI and the RSEI using Landsat 8 RS imagery and the ASTER Global Water Database as benchmarks, the following conclusions are deduced [83,84]. A comparison of Figure 4(7,10) reveals that the ARSEI, while accurately representing vegetation greenness, also more precisely reflects the ecological quality of agricultural lands under the influence of soil salinity. Further comparison using Figure 4(8,11) indicates that the RSEI tends to overestimate the ecological degradation around salt lakes and does not adequately capture the impact of soil salinity on EQ. In contrast, Figure 4(8) illustrates how the ARSEI, through the CSI indicator, effectively represents the effect of soil salinity on EQ, showcasing the superiority of the CSI index within the ARSEI. Additionally, in high-value vegetation areas, as depicted in (b) and (c), the ARSEI finely details the ecological textures near river channels, identifying low-quality ecological points within the channels. Comparisons of Figure 4(9,12) demonstrate the ARSEI’s capability to comprehensively reflect the overall state of forests in the upper right section of the (c) test area, contrary to the RSEI, which displays fragmented patterns. A comparison with Landsat 8 imagery confirms that the ARSEI’s outcomes align more closely with the actual ground conditions. This analysis substantiates the accuracy and superiority of the ARSEI in assessing EQ.

3.1.2. Correlation Evaluation of Ecological Factors

The PCA results and the water threshold for each year from 2000 to 2020 can be found in the supplementary material (Table S2). The PCA results of five selected periods showed that the contribution rate of PC1 from 2000 to 2020 exceeded 88%, integrating most of the information on each ecological factor (

Table 4. The PCA results and the water threshold of five indicators.

Year	Parameters	EVI	LST	SPWI	NDBSI	CSI
2000	Water threshold	−0.0898
	Eigenvalues	0.1227	0.0073	0.0047	0.0031	0.0004
	Percentage variance	88.79%	5.31%	3.37%	2.27%	0.26%
	PC1	0.4346	−0.4346	0.4610	−0.4622	−0.4428
	PC2	−0.6974	−0.5770	0.2405	−0.1680	0.3077
	PC3	−0.1105	0.6405	0.4869	−0.5034	0.2952
	PC4	0.5590	−0.2563	0.0241	0.0362	0.7874
	PC5	−0.0110	−0.0477	−0.7014	−0.7096	0.0464
2005	Water threshold	−0.0585
	Eigenvalues	0.1233	0.0081	0.0038	0.0028	0.0003
	Percentage variance	89.18%	5.85%	2.74%	1.99%	0.24%
	PC1	0.4322	−0.4326	0.4618	−0.4628	−0.4456
	PC2	−0.6611	−0.5772	0.2473	−0.1844	0.3669
	PC3	−0.3892	0.6761	0.4086	−0.4525	−0.1404
	PC4	0.4739	0.1436	0.2443	−0.2218	0.8037
	PC5	0.0101	0.0442	0.7064	0.7056	−0.0338
2010	Water threshold	−0.0743
	Eigenvalues	0.1266	0.0062	0.0030	0.0026	0.0003
	Percentage variance	91.3%	4.49%	2.13%	1.84%	0.24%
	PC1	0.4399	−0.4363	0.4585	−0.4583	−0.4426
	PC2	−0.5726	−0.6220	0.2332	−0.1613	0.4527
	PC3	−0.1337	0.6360	0.4637	−0.5356	0.2751
	PC4	0.6788	−0.1176	−0.0190	0.0460	0.7232
	PC5	−0.0093	−0.0665	−0.7211	−0.6892	0.0228
2015	Water threshold	−0.0820
	Eigenvalues	0.1236	0.0081	0.0033	0.0030	0.0004
	Percentage variance	89.26%	5.87%	2.41%	2.18%	0.28%
	PC1	0.4300	−0.4372	0.4623	−0.4620	−0.4436
	PC2	−0.6871	−0.5502	0.2397	−0.2142	0.3491
	PC3	0.4232	−0.6712	−0.3542	0.4135	0.2718
	PC4	0.4049	0.2329	0.3007	−0.2907	0.7790
	PC5	0.0083	−0.0365	−0.7162	−0.6965	0.0231
2020	Water threshold	−0.0509
	Eigenvalues	0.1248	0.0074	0.0035	0.0023	0.0004
	Percentage variance	90.21%	5.31%	2.54%	1.65%	0.28%
	PC1	0.4400	−0.4245	0.4622	−0.4607	−0.4476
	PC2	−0.5250	−0.7469	0.1286	−0.0564	0.3830
	PC3	0.5418	−0.4915	−0.4397	0.5210	0.0084
	PC4	0.4870	0.1249	0.2080	−0.2261	0.8080
	PC5	−0.0015	−0.0687	−0.7302	−0.6797	0.0093

). The loadings of each ecological factor in PC1 are stable and regular across years, with positive values for the EVI and the SPWI indicating that greenness and soil moisture had a positive effect on the ecology of the IRB and negative values for the LST, NDBSI, and CSI indicating that temperature, drought, and soil salinity had a negative effect on the ecology of the IRB. In contrast, the loadings of the other principal components were unstable and could not be used to characterize ecological importance.

To comprehensively assess the influence of the selected ecological factors on this catchment and validate their correlation with the PCA results, we performed a Pearson correlation analysis between the five ecological factors and the ARSEI (Figure 5). The correlation analysis results are illustrated in Figure 4. Notably, the EVI, SPWI, and ARSEI exhibited a positive correlation, while the LST, NDBSI, CSI, and ARSEI demonstrated a negative correlation, aligning with the PCA findings. The average correlation across the five periods for each ecological factor exceeded 0.69, and the correlation between each indicator and the ARSEI surpassed 0.79. These results affirm that the ARSEI effectively encapsulates the comprehensive information of each ecological factor, providing a holistic reflection of the overall EQ.

3.2. Ecological Environment Measurement Based on the ARSEI

The statistical results show that the ARSEI values in 2000–2020 are distributed in the range [0, 1], with an overall mean of around 0.54, a median in the range 0.48–0.58, and a standard deviation of 0.23, which is subject to some fluctuation (Figure 6). In terms of the density frequency distribution, the high density of the ARSEI values is concentrated in the range of 0.11–0.90, showing a bimodal distribution. This feature statistically highlights an imbalance in the EQ of the IRB.

To quantitatively depict the spatial EQ of the study area, the ARSEI results were categorized into five levels (I-V) for precise quantification and visualization: poor (0.0–0.2), fair (0.2–0.4), average (0.4–0.6), good (0.6–0.8), and excellent (0.8–1.0) (Figure 7).

In general, there is evident spatial heterogeneity in the ecological environment of the IRB. Using approximately 55°N as a boundary, EQ is generally better in the north, while EQ is generally worse in the south. The distribution pattern gradually improves from south-east to north-west. The areas with an overall better EQ are Russia and the northern slopes of the Altay Mountains in Xinjiang, China, where forests, herbaceous wetlands, and shrubs predominate. Areas of medium EQ are mainly located on the border between Russia and Kazakhstan, where grassland and cultivated land predominate. The EQ is poor in Kazakhstan and on the southern slopes of the Altay Mountains in the IRB. Combined with the land use classification, forests and shrubs are distributed in level V; shrubs, grasslands, and herbaceous wetlands in level Ⅳ; mainly grassland and cultivated land in levels Ⅲ and Ⅱ; and grassland, bare land, and alpine ice and snow in level Ⅰ (Figure 7f).

The study used a gravity transfer model to calculate the coordinates of the center of gravity for levels I–V between 2000 and 2020 and to plot the transfer path (Figure 8). The results show that the movement trends for levels I and II are from south-east to north-west. Conversely, the movement trends for levels III and IV are from south-west to north-east, while the transfer path of level V involves an initial movement to the south-east, followed by a shift to the north-west, and finally a downward movement from 2015 to 2020. The distances traveled for levels I to V are 1798 km, 199 km, 268 km, 165 km, and 98 km, respectively.

The transfer matrix further reflects the difference in the spatiotemporal distribution of the EQ pattern (Figure 9). Over the last 20 years, EQ in the IRB has shown an overall trend of deterioration followed by improvement. From 2000 to 2005, there were few transfers between different classes, but after 2005, transfers between classes became more frequent. The level IV and V areas showed relatively stable changes; the level III area showed a deteriorating trend; the level II area showed an increasing trend; and the level I area initially increased but then decreased.

3.3. Time Series Analysis of Ecological Quality in the Irtysh River Basin

3.3.1. Time Series Stability Analysis

To study the spatial variation characteristics of EQ from 2000 to 2020 in more detail, a pixel-scale spatial measurement and a temporal evolution stability simulation were carried out to obtain the ARSEI coefficient of variation (CV) (Figure 10). Its values are mainly concentrated at 0~0.2, accounting for 87.5%, indicating that the overall fluctuation of the ARSEI is small and in a stable state.

The area distribution across different variability classes is as follows, in descending order of proportion: low variability at 28.43%, medium variability at 20.16%, lower variability at 19.78%, higher variability at 19.13%, and high variability at 12.50%. The range of the coefficient of variation is 0.0104~1.0294, the mean is 0.1115, and the standard deviation is 0.0734. The spatial distribution of each fluctuation level shows that there is significant spatial heterogeneity in the EQ of the catchment. The low-variability areas are mainly distributed in the West-Siberian plains of Russia in the northern part of the catchment, and the high-variability areas are mainly distributed in the northern Kazakh hills and the Turgai depression in Kazakhstan. The forest ecosystem has a low degree of variability, while the grassland and cropland ecosystems have a high degree of variability. However, even forest areas with low EQ variability experience local fluctuations. For example, the area change in LUCC shows that the extent of evergreen needle-leaf forests decreases and the area of shrublands and savannas increases from 2001 to 2020, causing major changes in EQ in the area in the northern part of A.

3.3.2. Pixel-Based Analysis of Spatiotemporal Changes in Ecological Quality

This study uses LT-GEE to analyze the greatest loss and gain of spatiotemporal change in EQ on a pixel basis (Figure 11). The interannual loss analysis of EQ shows that the distribution areas of levels II (fair) and III (average) are the main areas where loss and gain occur. The EQ inter-annual loss analysis reveals that the south-western part of the basin, particularly the northern region of Kazakhstan, experiences more prolonged and substantial losses. In the time series stability analysis, region A in the northern part of the basin is identified as an area with severe and sustained losses over approximately 20 years (Figure 11c). The loss time plot shows that the decline in EQ is more pronounced in 2000–2005, and the end of the loss is more pronounced in 2005–2012 (Figure 11a–c). The analysis of the inter-annual growth of EQ shows that the area of gain growth coincides very closely with the area of loss. In the Kazakhstan-Russia border region, the period of gain is around 2000–2002, and within Kazakhstan, the period of gain is around 2012–2020, with most of the growth lasting only one year (Figure 11d–f).

3.3.3. Future Trend and Persistence Analysis

The Sen-MK and the Hurst index were used to quantify the spatial distribution characteristics of trend changes and the persistence of EQ in the study area (Figure 12). The results of future trend analysis, at a confidence level of α = 0.05, reveal distinct patterns: areas exhibiting an advancing EQ trend (composed of significant and slight improvements) constitute 14.71%, those with a consistent EQ trend (remaining stable without change) make up 76.15%, and regions witnessing a declining EQ trend (marked by significant and slight deteriorations) account for 9.14%. The Hurst index ranges from 0.053 to 1, with an average of 0.433. The portion with values below 0.5 represents 76.58% of the total area, while values above 0.5 represent only 23.42%. This suggests a robust negative persistence in the trend of EQ change within the IRB, implying that the trajectory of EQ change in most catchments will be opposite to that of the past. The regions with a Hurst index of less than 0.5 are mainly located in areas of levels Ⅰ-Ⅲ, indicating that the EQ of the catchment is gradually improving. It is worth noting that in the southern part of the Pavlodar region, the EQ classification typically falls within levels I or II, with a Hurst value less than 0.5. This suggests that the region is likely to encounter challenges in transitioning toward an improving trend.

The results of the Sen-MK analysis were combined with the Hurst index to analyze the persistence of future changes in the trend in EQ in the study area (Table 3). The analysis revealed the following distribution of results: 76.15% for persistent and stable conditions, 7.47% for persistent and significant improvement, and 4.36% for persistent and significant degradation. This was followed by 3.14% showing a trend from persistent and significant improvement to persistent and significant degradation, 3.08% for persistent and slight improvement, 2.55% for persistent and slight degradation, 1.53% for a shift from persistent and significant degradation to persistent and significant improvement, 1.02% for transitioning from persistent and slight improvement to persistent and slight degradation, and finally, 0.71% for persistent slight degradation evolving into persistent slight improvement (Figure 13). This indicates that the future trend of the ARSEI in the basin is dominated by persistent stability in the average, with an overall improvement in EQ. The spatial distribution of sustainable changes in future trends in EQ shows that the West Siberian Plain and the northern slopes of the Altay Mountains and East Kazakhstan Oblast continue to show significant improvements in EQ, while the northern part of Pavlodar Oblast and Kostanay Oblast continue to show significant degradation.

3.4. Analysis of Driving Factors for Ecological Quality in the Irtysh River Basin

This paper uses the PLUS model to analyze the driving forces of EQ changes at levels I–V). The PLUS model uses uniform sampling with 50 regression trees in Random Forest Regression (RFR) and a sampling rate of 0.1. The root mean square error (RMSE) values for levels I–V are 0.057, 0.091, 0.091, 0.097, and 0.079, respectively. In addition, out-of-bag (OOB) RMSE values of 0.165, 0.263, 0.264, 0.280, and 0.227, respectively, are also reported for the same categories. These data clearly demonstrate the effectiveness of the PLUS model, with minimal error in assessing the factors contributing to different types of environmental quality.

The poor area expansion is mainly driven by precipitation, elevation, and population density (Figure 14). Expansion occurs mainly in the southernmost part of the West-Siberian Plain, situated south of the Turgai Depression (Figure 15a). The main drivers of the increase in the fair are precipitation, population density, and soil pH (Figure 14). It increases mainly in the northern part of the territory of Kazakhstan, where the average annual precipitation from June to September is low (Figure 15b). The main drivers of average area expansion are population density, precipitation, and elevation (Figure 14). It increases mainly in the less densely populated areas, with an overall trend from northeast to southwest (Figure 15c). The main drivers of good area expansion are soil pH, precipitation, and elevation (Figure 14). The expansion area is in the north-central part of the basin and is characterized by slightly acidic soils with a pH range of 5.7–6.4. It shows a spatial pattern of increasing and then decreasing from south to north (Figure 15d). The predominant factors contributing to its excellent distribution are soil pH, precipitation, and elevation (Figure 14). It occurs mainly in the northern part of the watershed (Figure 15e).

4. Discussion

4.1. Rationality and Superiority of the ARSEI

In this paper, the GEE platform is used to calculate the ARSEI, which avoids the influence of pre-processing such as stitching and declouding on the results in software such as ENVI. The RSEI model has strong robustness and has been widely used in areas with different geographical conditions, with more than 1200 citations of related papers [85]. The ARSEI model improves the RSEI in three ways: (1) using the EVI to replace the NDVI; (2) adding a salinity factor that characterizes the soil environment for vegetation growth; and (3) using the Otsu algorithm to generate dynamic thresholds, which improve the accuracy of an EQ assessment using the IRB as an example. The ARSEI provides a more objective description of environmental realism than the RSEI results. The rationality of selected parameters has been clearly established with the Pearson correlation coefficient method. The analysis of the ARSEI results shows that it is not only consistent with the laws of ecological development (greenness, humidity, and productivity play a positive role in EQ; heat, dryness, and salinity play a negative role) but also with the results of various remote sensing ecological indices.

4.2. Spatiotemporal Pattern and Evolutionary Trend of Ecological Quality

The analysis of the spatial pattern and future trends of EQ in the IRB shows that the distribution of EQ exhibits pronounced spatial heterogeneity, with an overall spatial increase from south to north. The time series stability analysis shows that the EQ of the catchment is changing little and remains stable. The pixel-based analysis of spatiotemporal changes also indicates that within regions experiencing loss, there is a continuous process of gain. Simultaneously, the persistence of future trends suggests that EQ in the basin will show a slight trend towards improvement from a stable base.

The distribution areas of levels Ⅰ and Ⅱ are mainly located in arid areas of Central Asia, mainly grassland and cultivated land. Of these, EQ is very poor in the western part of Akmola Oblast, the north-western and western parts of Karaganda Oblast, and the areas around Zayan Lake in East Kazakhstan Oblast. The North Kazakhstan Oblast is the most stable and best-quality land use type in Kazakhstan, with an overall medium and above EQ. The East Kazakhstan Oblast and Pavlodar Oblast are typical industrial centers with numerous coal, non-ferrous metallurgy, oil, and chemical industry plants. These industries have a serious impact on soil, air, and rivers. The rapid development of these industries has led to serious environmental pollution [86]. A large amount of industrial wastewater is feared to be discharged into the Irtysh River every year, seriously affecting the EQ of the basin [87]. The Kostanay Oblast is the most typical area in the process of reclamation, abandonment, and then reclamation. According to the statistics, in the period 1988–2013, it had abandoned farmland, accounting for about 40% of the total area of its oblast [88]. Since 2000, exclusively land in the Kostanay Oblast with superior soil quality has been converted, resulting in further widespread land abandonment and degradation [88,89]. As human society encroaches on and alters the natural environment of this region, further ecological degradation is inevitable. The distribution area of level Ⅲ is mainly located on the border between Russia and Kazakhstan, where the main land use types are natural grassland and artificial cropland. The largest land type in North Kazakhstan Oblast is arable land, accounting for more than half of the total area. Grassland is another major land use type, accounting for about 30% of the total area. The conversion between these two types accelerates soil salinization, leading to ecological degradation. Abandoned arable land can accelerate the process of soil salinization through the growth of weeds and shrubs, leading to ecosystem degradation. In addition, the steps of southern West Siberia face a serious threat of water scarcity [90]. The distribution areas of levels Ⅳ and Ⅴ are mainly located in northern West Siberia. There are rich forests with high species richness that ensure the EQ of the area.

4.3. Analysis of Driving Factors for Ecological Quality

The driver factor analysis reveals that the ARSEI trend is consistent with the distribution of basin precipitation and soil pH and is negatively correlated with population distribution, while there is a clear gradient in the impact of altitude on EQ. Precipitation in the southern region indicates a decrease in EQ, while increased precipitation and slightly acidic soils favoring alpine forest vegetation in the northern region indicate an improvement in EQ. Human activities contribute to ecological degradation by converting land for agriculture and construction.

The reason for the decrease in EQ in the poor region is that the temperature increases and precipitation decreases, thereby exacerbating drought [91]. The region encompasses a mild gradient, large areas of cultivable land, and grassy terrain, rendering it a fitting ground for human settlement. A change in agricultural policy has markedly altered the landscape, transitioning cultivation grounds to grazing areas—a progression from cultivable to fallow, cut, pasture, and finally grazing lands, specifically within Akmola Oblast and Karaganda Oblast. Overutilization of grassland is a lasting challenge, as it is less productive than arable land. This results in alterations to the grassland ecosystem, causing a reduction in vegetation coverage, which subsequently becomes arid and saline. This leads to a significant decline in EQ [88,92].

The fair region has a typical continental climate, characterized not only by high radiation and large intra-annual temperature variations but also by considerable rainfall variability. It has many heavy rainfall events in the summer, sometimes reaching 24–65 mm or more in 24 hours. The heavy rainfall increases erosion and soil loss, while the high summer temperatures in the region lead to strong wind erosion due to the high evaporation of water after heavy rainfall [93]. According to Pavlodar Oblast statistics, 74% of arable land is affected by wind erosion. In Altayskiy Kray, 95% of the barren arable land is in the dry steppe zone bordering Kazakhstan [94]. This is because higher soil pH leads to lower species richness in the arid grasslands of northern Kazakhstan and south-western Siberia [95]. At the same time, the needs of people in Kazakhstan for economic and social development have led to overgrazing, land expansion, and mineral exploitation, resulting in increased soil erosion and pollution. Both adverse climatic and human influences are not conducive to improving EQ.

The impact of population density on EQ in average regions results from the development and use of urban and agricultural land. Kazakhstan’s urbanization efforts in the 21st century have exerted significant pressure on the ecological environment [96]. After independence, Kazakhstan began to restore natural pastures and rangeland that had been haphazardly cultivated. In the period 2000–2018, abandoned land and state-owned pastures were brought back to use. In 2006, Kazakhstan introduced the Sustainable Development Action Plan (2007–2024) and has continued to strengthen environmental protection efforts and improve relevant laws and regulations. It has achieved some useful results in ecological management and restoration, such as improving saline soils, seasonal grazing, and establishing nature reserves. As a result, the overall use of natural resources is better than in the past, and ecological conditions have improved.

The expansion of good regions is located in the north-central part of the basin, where increased precipitation has reduced soil pH to a slightly acidic range of 5.7–6.4, creating an environment for increased species richness [97]. The rate of forest restoration in Western Siberia is increasing because the President of the Russian Federation (No. 474 of 21 July 2020) established twelve national projects in 2018 [98]. In the Omsk region alone, about 8900 ha of forest will be restored and planted in 2020, increasing the forest cover and effectively restoring areas where forests have been degraded to shrubs or grassland because of overharvesting.

The excellent region is mainly distributed in the northern part of the IRB, which is the West Siberian Plain, with a large amount of coniferous and deciduous forests. The coniferous forests include spruce, pine, fir, and larch, while the deciduous forests include birch, poplar, willow, and kaiju. Research has shown that spruce and birch are suitable for soils with a pH of 5 to 7 and that this region provides a good ecological basis for forest regeneration [99]. The average precipitation from June to September every year is relatively high, and the abundant precipitation resources will further promote the improvement in EQ. It will further regulate deforestation in the region under Russia’s ecological protection policy. It will also find the best balance between EQ and economic benefits in the research area.

4.4. Limitations of the Study and Future Work

To complete an EQ assessment of the IRB, MODIS products were selected as the data source to reduce the computational requirements. However, it is important to recognize that the resolution of MODIS can affect the accuracy of EQ assessments. The selection of drivers is dominated by natural factors such as soil properties, rainfall, and topography. These indicators are key elements that influence EQ. The ARSEI shows that greenness and soil moisture are the most important contributors to EQ. There is a strong relationship between vegetation and soil. Soil and topography serve as the basis for vegetation growth, while precipitation affects its development. Therefore, we analyzed the driving mechanism of the above factors on the change in EQ. At the same time, these soil characteristics are easily obtained with GEE. Population density was chosen as a driving factor because the impact of human activities on the environment is becoming increasingly apparent. However, human activities are complex, and population density alone may not fully reflect the impact of human activities on the study area. Future work will focus on using higher-resolution datasets to improve the accuracy of the assessment and exploring the coupling relationship between human activities and the EQ of the IRB through additional human factors.

5. Conclusions

The ARSEI comprehensively, objectively, and quantitatively reflects the EQ of the IRB from 2000 to 2020 using remote sensing technology through the framework. Analysis of the temporal and spatial variations in its results leads to the following conclusions:

(1) The correlation between the ARSEI and each indicator reaches more than 0.88, which can effectively characterize EQ. Among them, EVI can express the vegetation cover of the catchment more accurately than NDVI. The EQ of the IRB displays significant spatial variation, with more impressive EQ in the northern areas and relatively inferior quality in the southern regions. The gradual improvement from south to north is exhibited by this spatial pattern.

(2) The ARSEI was used to assess the EQ of the IRB watershed, and the overall situation showed a stable state with small fluctuations. In the current research context, the continuity analysis of future trend changes in the EQ of the IRB watershed shows that it will continue to improve on a stable basis. Among these, arid grasslands are an area of significant and continuous improvement.

(3) Precipitation, soil pH, elevation, and human population are the main factors influencing the EQ of the IRB. Due to spatial heterogeneity, the driving factors for different EQ classes are different. Among them, precipitation has a more significant impact on EQ.

This study presents the ARSEI-based evaluation outcomes of EQ and its prospective alteration patterns, as well as investigating the drivers of diverse EQ categories. The findings of this investigation provide a theoretical framework that may be utilized for the formulation of ecological preservation, management, governance, and sustainable development strategies and practices.